*6.1. Double-Frequency Tuning for Gain Factor Optimization*

Assume that a 250 THz (λ*f ree space* = 1.2 μm) infra-red stimulus wave *Est* and a high power pump wave *Ehp* that is composed of two high-intensity ultrashort pulses (frequencies are to be determined) are propagating inside a low-loss (high Q) cavity that has two reflecting walls (see Figure 5). The wall on the left side can be thought as an optical isolator, which fully transmits from its left side and almost fully reflects from its right side. The wall on the right side represents an optical band-pass filter with a frequency-dependent reflection coefficient Γ(*f*). Both waves are generated at x = 0 μm, at the time instant t = 0 s. The waves and the parameters of the gain medium are as given below:

$$E\_{hp}(\mathbf{x} = 0\mu m, t) = \sum\_{i=1}^{2} A\_i \cos(2\pi f\_i t + \psi\_i)(u(t) - u(t - \Delta T\_i)) \quad (u(t): \text{Unit step function})$$

$$\text{Where } A\_1 = 2 \times 10^8, \ A\_2 = 1.5 \times 10^8, \ \psi\_1 = 0, \ \psi\_2 = 0, \ \Delta T\_1 = 0.5\mu\text{s}, \ \Delta T\_2 = 1\text{ps}$$

$$E\_{xt}(\mathbf{x} = 0\mu m, t) = 1 \times \sin(2\pi (2.5 \times 10^{14})t) \text{ V/m}, \quad \text{for} \quad 0 \le t \le 10\text{ps}$$

$$\text{Dielectric constant of the gain medium} = \varepsilon\_{0\text{s}} = 1 + \chi = 12 \quad (\mu\_I = 1)$$

$$\text{Resonance frequency of the gain medium}: \quad f\_0 = 500T\hbar\text{Hz}$$

$$\text{Damping coefficient of the gain medium}: \quad \gamma = 1 \times 10^9 H\text{z}$$

$$\text{Time interval and duration of simulation}: \quad 0 \le t \le 10\mu\text{s}$$

$$\text{Spath range of the gain medium}: 0\mu\text{m} < \text{x} < 10\mu\text{m}$$

$$\text{Rigbt cavity wall location}: \quad \mathbf{x} = 10\mu\text{m}; \quad \text{Left cavity wall location}: \quad \mathbf{x} = 0\mu\text{m}$$

*Electron density o f the gain medium* : *<sup>N</sup>* = 3.5 <sup>×</sup> 1028/*m*3; *Atomic diameter* : *<sup>d</sup>* = 0.3 *nm*

**Figure 5.** Configuration of the cavity and the parameters for Section 6.1.

**Our problem:** find the optimal pump wave pulse frequencies *fp*1, *fp*<sup>2</sup> that maximize the magnitude of the monochromatic stimulus wave in the cavity ( *Est*( *fst* <sup>=</sup> <sup>250</sup> *THz*) ), for 10 THz <sup>&</sup>lt;  *fp*1, *fp*<sup>2</sup> < 1000*THz* (THz to UV), and for 0μ*m* < *x* < 10μ*m*, 0 ≤ *t* ≤ 10*ps*, such that

∇<sup>2</sup>*Ehp fp*1, *fp*<sup>2</sup> − μ0ε<sup>∞</sup> ∂<sup>2</sup>*Ehp fp*1, *fp*<sup>2</sup> <sup>∂</sup>*t*<sup>2</sup> <sup>=</sup> <sup>μ</sup>0<sup>σ</sup> ∂*Ehp fp*1, *fp*<sup>2</sup> <sup>∂</sup>*<sup>t</sup>* <sup>+</sup> <sup>μ</sup><sup>0</sup> ∂<sup>2</sup>*Php* ∂*t*<sup>2</sup> ∂<sup>2</sup>*Php* <sup>∂</sup>*t*<sup>2</sup> <sup>+</sup> <sup>γ</sup> ∂*Php* <sup>∂</sup>*<sup>t</sup>* <sup>+</sup> <sup>ω</sup><sup>0</sup> 2 *Php* <sup>−</sup> <sup>ω</sup><sup>0</sup> 2 *Ned Php*<sup>2</sup> <sup>−</sup> <sup>ω</sup><sup>0</sup> 2 *N*2*e*2*d*<sup>2</sup> *Php*<sup>3</sup> <sup>=</sup> *Ne*<sup>2</sup> *<sup>m</sup> Ehp fp*1, *fp*<sup>2</sup> . ∇2*Est fp*1, *fp*<sup>2</sup> − μ0ε<sup>∞</sup> ∂2*Est fp*1, *fp*<sup>2</sup> <sup>∂</sup>*t*<sup>2</sup> <sup>=</sup> <sup>μ</sup>0<sup>σ</sup> ∂*Est fp*1, *fp*<sup>2</sup> <sup>∂</sup>*<sup>t</sup>* <sup>+</sup> <sup>μ</sup><sup>0</sup> ∂2*Pst* ∂*t*<sup>2</sup> ∂2(*Pst*) <sup>∂</sup>*t*<sup>2</sup> + <sup>γ</sup> ∂(*Pst*) <sup>∂</sup>*<sup>t</sup>* + ω<sup>0</sup> <sup>2</sup>(*Pst*) <sup>−</sup> <sup>ω</sup><sup>0</sup> 2 *Ned Pst*<sup>2</sup> + 2*PstPhp* <sup>−</sup> <sup>ω</sup><sup>0</sup> 2 *N*2*e*2*d*<sup>2</sup> *Pst*<sup>3</sup> + 3*Pst*<sup>2</sup>*Php* + 3*PstPhp*<sup>2</sup> = *Ne*<sup>2</sup> *<sup>m</sup> Est fp*1, *fp*<sup>2</sup> 

Our aim is to maximize the magnitude of the stimulus wave at its original frequency *fst* = 250*THz* (monochromatic form). This is a precaution against any degree of spectral broadening that the stimulus wave may go through while being amplified. Therefore, our cost function is chosen as (at any spatial point *x* = *x* inside the cavity):

$$\begin{array}{lcl} Q = & \left| E\_{\rm st} (f\_{\rm st} = 250THz) \right| \\ &= \left| \int\_{2.5 \times 10^{14} - \Delta f}^{2.5 \times 10^{14} + \Delta f} \left\{ \int\_{0}^{\Delta T} \langle E\_{\rm st} (\mathbf{x} = \mathbf{x}', t) e^{-i (2\pi \Omega)t} | dt \right\} e^{i (2\pi \Omega)t} d\Omega \right| \end{array}$$

 

*where* Δ*T* = 10*ps*, 0 ≤ *t* ≤ 10*ps*, 2.5 <sup>×</sup> <sup>10</sup><sup>14</sup> <sup>−</sup> <sup>Δ</sup>*<sup>f</sup>* < Ω < 2.5 <sup>×</sup> 1014 + <sup>Δ</sup>*<sup>f</sup>* , Δ*f* = 1*THz* **Initial conditions:**

$$P\_{\rm lp}(\mathbf{x},0) = P\_{\rm lp} \, ^\prime(\mathbf{x},0) = E\_{\rm lp}(\mathbf{x},0) = E\_{\rm lp} \, ^\prime(\mathbf{x},0) = P\_{\rm sl}(\mathbf{x},0) = P\_{\rm sl} \, ^\prime(\mathbf{x},0) = E\_{\rm sl} \, (\mathbf{x},0) = E\_{\rm sl} \, ^\prime(\mathbf{x},0) = 0$$

**Boundary and excitation conditions:**

$$E\_{\rm hp}(\mathbf{x} = 0 \,\mu m, t) = \sum\_{i=1}^{2} A\_i \cos(2\pi f\_i t + \psi\_i) (u(t) - u(t - \Delta T\_i))^2$$

where *<sup>A</sup>*<sup>1</sup> = <sup>2</sup> <sup>×</sup> <sup>10</sup>14, *<sup>A</sup>*<sup>2</sup> = 1.5 <sup>×</sup> <sup>10</sup>14, <sup>ψ</sup><sup>1</sup> = 0, <sup>ψ</sup><sup>2</sup> = 0, <sup>Δ</sup>*T*<sup>1</sup> = 0.5*ns*, <sup>Δ</sup>*T*<sup>2</sup> = <sup>1</sup>*ns*

$$E\_{st}(\mathbf{x} = 0 \mu m, t) = 1 \times \frac{\sin \left( 2\pi \left( 2.5 \times 10^{14} \right) t \right)}{m}, \quad for \quad 0 \le t \le 10 ps$$

$$E\_{lp}(\mathbf{x} = 15\mu m, t) = E\_{st}(\mathbf{x} = 15\mu m, t) = 0 \qquad for \; 0 < t < 10ps$$

**Absorbing boundary condition (perfectly matched layer):**

$$\sigma(\mathbf{x}) = \{ (\mathbf{x} - (L - \Delta)) \sigma\_0 \, , \quad (L - \Delta) \le x < L \quad \}, \text{ for } L = 15 \mu m, \quad \Delta = 2.5 \mu m, \; \sigma\_0 = 4.5 \times 10^8 \text{S/m}^2$$

**Optical isolator condition:** full reflection at x = 0μm

$$
\Gamma(\mathbf{x} = 0 \mu \mathbf{m}, t) = 1 \quad \text{(Reflection coefficient is equal to 1)}.
$$

**Optical bandpass filter condition:** frequency dependent reflection at x = 10μm

$$\left|\Gamma(f')\right| = 1 - e^{-\left(\frac{\left(f'-2SOTx\right)}{\sqrt{2}Thx}\right)^2}$$

**Cost function to be maximized:**

$$\mathcal{Q}\left(f\_{p1}, f\_{p2}\right) = \left| \mathbb{E}\_{\mathcal{S}\mathcal{I}}\left(f\_{\mathcal{S}\mathcal{I}} = 250THz\right) \right| - \delta\_1 \left(f\_{p1} - f\_{\text{max}}\right)^2 - \delta\_2 \left(f\_{\text{min}} - f\_{p1}\right)^2 - \delta\_3 \left(f\_{p2} - f\_{\text{max}}\right)^2 - \delta\_4 \left(f\_{\text{min}} - f\_{p2}\right)^2$$

where

$$\delta\_{1} = \left\{ \begin{array}{c} 0 & \text{if } \begin{array}{l} f \\ \frac{\left[E\_{\text{af}}\left(\int\_{f} - 250Hz\right)\right]}{10^{2}} \end{array} \right\} \delta\_{1} = \left\{ \begin{array}{l} 0 & \text{if } \begin{array}{l} f \\ \frac{\left[E\_{\text{af}}\left(\int\_{f} - 250Hz\right)\right]}{10^{2}} \end{array} \right\} \end{array} \right\}$$

$$\delta\_{3} = \left\{ \begin{array}{l} 0 & \text{if } \begin{array}{l} f \\ \frac{\left[E\_{\text{af}}\left(\int\_{f} - 250Hz\right)\right]}{10^{2}} \end{array} \end{array} \right\} \delta\_{4} = \left\{ \begin{array}{l} 0 & \text{if } \begin{array}{l} f & f\_{p} \geq f\_{\text{min}} \\ \frac{\left[E\_{\text{af}}\left(\int\_{f} - 250Hz\right)\right]}{10^{2}} \end{array} \end{array} \right\}$$

**Optimization algorithm (BFGS):**

$$\text{Set } H\_1 = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix} \cdot f\_{p1,0} = 100THz, \quad f\_{p1,1} = 102THz,$$

$$f\_{p2,0} = 100THz, \qquad f\_{p2,1} = 103THz, \quad \alpha\_1 = 0.5$$

$$\nabla \mathbf{Q}\_k = \begin{bmatrix} \frac{Q(f\_{p1,k} \cdot f\_{p2,k}) - Q(f\_{p1,k-1} \cdot f\_{p2,k})}{f\_{p1,k} - f\_{p2,k-1}}\\ \frac{Q(f\_{p1,k} \cdot f\_{p2,k}) - Q(f\_{p1,k} \cdot f\_{p2,k-1})}{f\_{p2,k} - f\_{p2,k-1}} \end{bmatrix}$$

$$\mathbf{p}\_k = -\mathbf{H}\_k \nabla \mathbf{Q}\_k$$

$$f\_{p,\ k+1} = f\_{p,k} + \alpha\_k \mathbf{p}\_k, f\_{p,k} = \begin{bmatrix} f\_{p1,k} \\ f\_{p2,k} \end{bmatrix}$$

$$\mathbf{s}\_k = f\_{p,k+1} - f\_{p,k}$$

$$\nabla \mathbf{Q}\_{k+1} = \begin{bmatrix} \frac{Q(f\_{p1,k+1}, f\_{p2,k}) - Q(f\_{p1,k}, f\_{p2,k})}{f\_{p1,k+1} - f\_{p1,k}} \\ \frac{Q(f\_{p1,k}, f\_{p2,k}) - Q(f\_{p1,k}, f\_{p2,k})}{f\_{p2,k+1} - f\_{p2,k}} \end{bmatrix}$$

$$\mathbf{y}\_k = \nabla \mathbf{Q}\_{k+1} - \nabla \mathbf{Q}\_k$$

$$\rho\_k = \frac{1}{y\_k^T \mathbf{s}\_k}$$

$$\mathbf{H}\_{k+1} = \left(\mathbf{I} - \rho\_k \mathbf{s}\_k \mathbf{y}\_k^T\right) \mathbf{H}\_k \left(\mathbf{I} - \rho\_k \mathbf{y}\_k \mathbf{s}\_k^T\right) + \rho\_k \mathbf{s}\_k \mathbf{s}\_k^T \text{ (BFGS)}$$

$$\mathbf{I}: \text{Identity matrix}$$

In order to satisfy the Wolfe conditions, the step size at each iteration is chosen as:

$$\alpha\_k = \mathbb{C}^{(\lfloor \log \rfloor \frac{\mathbb{Q}(f\_{p,k})}{\mathbb{Q}(f\_{p,k}) - \mathbb{Q}(f\_{p,k-1})} \mid) \; / \; (\mid \frac{\mathbb{Q}(f\_{p,k})}{\mathbb{Q}(f\_{p,k}) - \mathbb{Q}(f\_{p,k-1})} \mid)} \; \tag{19}$$

where *C* is just a constant (1 < *C* < 1.5) and α*<sup>k</sup>* is the step size at iteration k. This formula (Equation (19)) was determined by trial and error and was found to satisfy Wolfe's conditions automatically at each iteration. This saves us from the huge computational cost of running another iteration loop to determine the step size at each iteration of the optimization process. In this simulation C = 1.445. Based on the above formulations, the maximum stimulus wave amplitude that has been reached in the cavity (for 0 < t < 10 ps) is determined as *Gainmax* = *Est*( *fst* <sup>=</sup> <sup>250</sup>*THz*) *max* <sup>=</sup> 4.67 <sup>×</sup> <sup>10</sup>8*V*/*m*, which corresponds to the frequencies *fp*<sup>1</sup> = 387.2 *THz*, *fp*<sup>2</sup> = 319.4 *THz* (see Table 1), or to the free space wavelengths λ*p*<sup>1</sup> = 0.939 μm, λ*p*<sup>2</sup> = 0.775 μ*m* as the ultrashort pulses can be practically generated outside the resonator. In this simulation, the gain factor of the amplification is defined as:

$$\begin{array}{lcl} \text{Gain}\_{\text{max}} &= \left| E\_{\text{st}} (f\_{\text{st}} = 250Thz) \right|\_{\text{max}} \\ &= \frac{\text{Amplitude of } f \text{ the stimulus curve at } t = t\_{\text{max}} \text{ for } x = x'}{Amplitude \ of \text{ the stimulus curve at } t = 0 \text{ } for \ x = x'}. \end{array}$$

where *x* is chosen as an intra-cavity point (*x* = 5.73μ*m*) and *tmax* = 10 *picoseconds*.

$$\mathcal{W}\_{\varepsilon,p} = \text{Stored electric energy density via pump wave} = \frac{1}{2}\varepsilon\_{\text{os}}E\_{pump}{}^2 + \frac{1}{2}E\_{pump}P\_{pump}\left(\frac{\text{Joules}}{m^3}\right)$$

*Ppump* : *Polarization density created by the pump wave Coulomb m*<sup>2</sup> , *Epump* : *Pump wave electric field intensity*



**Table 1.** BFGS algorithm-based optimization process.

As we can see from Table 1, the optimal ultrashort pulse frequencies correspond to very high stored electric energy density and high polarization density. The stored electric energy density and the polarization density must be simultaneously high for a significant stimulus wave amplification. The stored electric energy density indicates the achievable order of stimulus wave amplification [14,15], and the polarization density acts as a coupling coefficient, which is a measure of how much stored electric energy can be coupled to the stimulus wave.

The time variation of the spectrally broadened (polychromatic) stimulus wave between t = 6.6 picoseconds and t = 10 picoseconds is shown in Figure 6. From the figure, we can see that the polychromatic stimulus wave reaches an amplitude of approximately 8 <sup>×</sup> 108 *<sup>V</sup>*/*m*.

**Figure 6.** Stimulus wave amplification (in polychromatic form) inside the cavity at x = 5.73 μm.
