**3. Wave Propagation in Non-Linear Dispersive Media**

As previously stated, non-linear interaction takes place when a wave that propagates in a medium have a very high intensity. This degree of non-linearity inducing high intensities is usually possible with ultra-short duration pulses and can be practically generated from a mode locked laser or a Q-switched laser, which have durations on the scale of picoseconds or femtoseconds. Since such high-intensity pulses have ultra-short durations, we must perform a dispersion analysis.

Impulse response of the polarization density of many dielectric media last much longer in duration than the pulse durations of such high intensity pulses [8,9].

The inclusion of dispersion in the analysis means that we need to solve for the polarization density of the interaction medium separately. In fact, it is much more accurate to solve the wave equation in parallel with the non-linear electron cloud motion equation for a non-linear dispersive medium, since parameters such as the resonance frequency, damping rate, atom density, and atomic diameter have typical values for most solid media, which makes the simulation results more realistic. To determine the time variation of the electric field in a non-linear dispersive medium, we need to solve the following equations (see Figure 1) [1]:

$$
\nabla^2 E - \mu\_0 \varepsilon\_{\alpha \nu} \frac{\partial^2 E}{\partial t^2} = \mu\_0 \sigma \frac{\partial E}{\partial t} + \mu\_0 \frac{\partial^2 P}{\partial t^2}.\tag{3}
$$

$$\frac{\partial^2 P}{\partial t^2} + \gamma \frac{\partial P}{\partial t} + \omega\_0 \,^2 P - \frac{\omega \nu\_0^2}{\text{N} \text{e} d} P^2 - \frac{\omega \nu\_0^2}{N^2 e^2 d^2} P^3 = \frac{\text{N} e^2}{m} E. \tag{4}$$

*P* : *Polarization density* , *e* : *Electron charge*, *m* : *Electron mass*, *E* : *Electric field* (*V*/*m*)

**Figure 1.** An interaction medium placed inside a cavity.

In Equation (4) we have accounted up to the third order of the polarization density as higher order terms are negligibly small. For a dielectric medium, the electric conductivity is almost zero (<sup>σ</sup> <sup>≈</sup> <sup>0</sup>). Some typical values for solid dielectric media are [1]; *Resonance f requency* : *<sup>f</sup>*<sup>0</sup> = 1.1 <sup>×</sup> <sup>10</sup>15*Hz, Damping rate* : <sup>γ</sup> = 1 <sup>×</sup> 109 *Hz*, *Electron density* : *<sup>N</sup>* = 3.5 <sup>×</sup> 1028/*m*3*, Atomic diameter* : *<sup>d</sup>* = 0.3 *nanometer*

Electromagnetic wave amplification by non-linear wave mixing involves two waves, namely the low intensity stimulus (input) wave to be amplified, and the high-intensity pump wave. Let us first consider the high-intensity pump wave *E*<sup>2</sup> without the presence of the low intensity stimulus wave *E*1. In this case, the pair of equations that model the propagation of *E*<sup>2</sup> is given as:

$$
\nabla^2(E\_2) - \mu\_0 \varepsilon\_{\infty} \frac{\partial^2(E\_2)}{\partial t^2} = \mu\_0 \sigma \frac{\partial(E\_2)}{\partial t} + \mu\_0 \frac{\partial^2 P\_2}{\partial t^2}.\tag{5a}
$$

$$\frac{\partial^2 P\_2}{\partial t^2} + \gamma \frac{\partial P\_2}{\partial t} + a\nu\_0^2 (P\_2) - \frac{a\nu\_0^2}{N\epsilon d} (P\_2)^2 - \frac{a\nu\_0^2}{N^2 \epsilon^2 d^2} (P\_2)^3 = \frac{N\epsilon^2}{m} (E\_2) \tag{5b}$$

*P*<sup>2</sup> : *Polarization density due to the electric field E*2, *P*<sup>1</sup> : *Polarization density due to the electric field E*<sup>1</sup>

Now assume that *E*<sup>1</sup> and *E*<sup>2</sup> are simultaneously propagating in the same medium, in that case the pair of equations that describe the total electric field propagation is written as:

$$
\nabla^2 (E\_1 + E\_2) - \mu\_0 \varepsilon\_{\rm os} \frac{\partial^2 (E\_1 + E\_2)}{\partial t^2} = \mu\_0 \sigma \frac{\partial (E\_1 + E\_2)}{\partial t} + \mu\_0 \frac{\partial^2 (P\_1 + P\_2)}{\partial t^2}.\tag{6a}
$$

$$\begin{split} \frac{\partial^2 (P\_1 + P\_2)}{\partial t^2} + \gamma \frac{\partial (P\_1 + P\_2)}{\partial t} + a \rho\_0^2 (P\_1 + P\_2) - \frac{\omega \eta^2}{\text{NaCl}} (P\_1 + P\_2)^2 - \frac{\omega \eta^2}{N^2 c^2 d^2} (P\_1 + P\_2)^3 \\ = \frac{Ne^2}{m} (E\_1 + E\_2) .\end{split} \tag{6b}$$

Our aim is to solve for the low amplitude wave *E*<sup>1</sup> in the presence of the high-amplitude wave *E*2, i.e., we want to solve for the stimulus wave *E*<sup>1</sup> when there is an energy coupling from *E*2. To model this problem, we subtract Equation (5a,b) from Equation (6a,b) respectively [1], which yields:

$$
\nabla^2(E\_1) - \mu\_0 \varepsilon\_{\alpha \nu} \frac{\partial^2(E\_1)}{\partial t^2} = \mu\_0 \sigma \frac{\partial(E\_1)}{\partial t} + \mu\_0 \frac{\partial^2(P\_1)}{\partial t^2}.\tag{7a}
$$

$$\frac{\partial^2(P\_1)}{\partial t^2} + \gamma \frac{\partial(P\_1)}{\partial t} + a\nu\_0^2(P\_1) - \frac{a\nu\_0^2}{\text{N}cd} \left\{ P\_1^2 + 2P\_1P\_2 \right\} - \frac{a\nu\_0^2}{\text{N}^2c^2d^2} \left\{ P\_1^3 + 3P\_1^2P\_2 + 3P\_1P\_2^2 \right\} = \frac{\text{N}e^2}{m}(E\_1) \tag{7b}$$

Equations (5a,b) and (7a,b) show that *E*<sup>2</sup> and *E*<sup>1</sup> are coupled to each other. The pump wave *E*<sup>2</sup> acts as a source field distribution for electromagnetic wave propagation, with the stimulus wave *E*<sup>1</sup> acting as the carrier distribution in reference to the model presented in [10]. Based on Equation (7a,b), we will examine if it is feasible to amplify the low-intensity electric field *E*1, via energy coupling from the microresonator that is energized by the high-intensity electric field *E*2(see Figure 2).

**Figure 2.** Concurrent propagation of the stimulus and the pump waves in a micro-resonator.

From Equation (7a), we can see that the stimulus wave *E*<sup>1</sup> has a source term *P*1, which is the polarization density of the stimulus wave. The source term *P*<sup>1</sup> is also dependent on the source term (polarization density) of the high-intensity pump wave *P*2. This means that in order to solve for the stimulus wave *E*1, we need to solve for the pump wave *E*<sup>2</sup> as the equations for the stimulus wave and the pump wave are coupled to each other via the coupling term *P*2. Therefore, the stimulus wave can be solved by simultaneously solving all these four equations (Equations (5a,b) and (7a,b)).

Our goal is to maximize the stimulus wave magnitude at a given stimulus wave frequency; max( *E*1(*fst* <sup>=</sup> *<sup>f</sup>* ) ). In order to achieve this, we will make a dispersion analysis that is based on the high-intensity pump wave frequency. By sweeping the excitation frequency of the pump wave (*fp*) in a large spectral interval that extends from the far-infrared region to the near-ultraviolet region, we can investigate the pump wave frequency dependent variation of the maximum stimulus wave magnitude for *fst* = *f* , and we can select the optimal pump wave frequency value that maximizes the magnitude of the stimulus wave for *fst* = *f* . Mathematically our optimization problem can be stated as:

> *Given that fmin* < *fp* < *fmax* ; *max E*1(*fst* <sup>=</sup> *<sup>f</sup>* ) *based on the f ollowing equations*

$$
\nabla^2 \left( E\_2 \left( f\_p \right) \right) - \mu\_0 \varepsilon\_{\infty} \frac{\partial^2 \left( E\_2 \left( f\_p \right) \right)}{\partial t^2} = \mu\_0 \sigma \frac{\partial \left( E\_2 \left( f\_p \right) \right)}{\partial t} + \mu\_0 \frac{\partial^2 P\_2}{\partial t^2} \tag{8a}
$$

$$\frac{\partial^2 P\_2}{\partial t^2} + \gamma \frac{\partial P\_2}{\partial t} + \omega\_0 \,^2(P\_2) - \frac{\omega\_0}{N \epsilon d}^2 (P\_2)^2 - \frac{\omega\_0^2}{N^2 e^2 d^2} (P\_2)^3 = \frac{N \epsilon^2}{m} \Big( E\_2 \Big(f\_p\Big) \Big) \tag{8b}$$

$$
\nabla^2 \left( E\_1(f\_p) \right) - \mu\_0 \varepsilon\_{\rm so} \frac{\partial^2 \left( E\_1(f\_p) \right)}{\partial t^2} = \mu\_0 \sigma \frac{\partial \left( E\_1(f\_p) \right)}{\partial t} + \mu\_0 \frac{\partial^2 (P\_1)}{\partial t^2} \tag{8c}
$$

$$\begin{split} \frac{\partial^2 (P\_1)}{\partial t^2} + \gamma \frac{\partial (P\_1)}{\partial t} + a \nu\_0^2 (P\_1) - \frac{a \nu\_0^2}{N \text{cl}} \Big\{ P\_1 \,^2 + 2P\_1 P\_2 \Big\} - \frac{\omega\_0^2}{N^2 c^2 d^2} \Big\{ P\_1 \,^3 + 3P\_1 \,^2 P\_2 + 3P\_1 P\_2 \Big\} \\ = \frac{\hbar \omega^2}{m} \Big( E\_1 (f\_p) \Big) \end{split} \tag{8d}$$

where

$$E\_2\{\mathbf{x} = \mathbf{x}\_{\text{input}}, t\} = A\_p \cos\{2\pi f\_p t + \psi\_p\} (\mathbf{u}(t) - \mathbf{u}(t - \Delta T\_p)) \, V/m \quad (\mathbf{u}(t) : \text{Unit step function})$$

$$E\_1\{\mathbf{x} = \mathbf{x}\_{\text{input}}, t\} = A\_{st} \cos\left(2\pi f\_{st} t + \psi\_{st}\right) (\mathbf{u}(t) - \mathbf{u}(t - \Delta T\_{st})) \frac{V}{m} \, , \, A\_p \gg A\_{st} \quad \Delta T\_p \ll \Delta T\_{st}$$

If we are using *N* ultrashort high-intensity pulse excitations to amplify the low-intensity stimulus wave,

$$E\_2\{\mathbf{x} = \mathbf{x}\_{\text{input}}, t\} = \sum\_{i=1}^{N} A\_i \cos(2\pi f\_i t + \psi\_i) \left( u(t) - u(t - \Delta T\_i) \right) V/m \quad \text{( $u(t)$  : llinite step function)}$$

$$E\_1\{\mathbf{x} = \mathbf{x}\_{\text{input}}, t\} = A\_{st} \cos(2\pi (f\_{st})t + \psi\_{st}) \left( u(t) - u(t - \Delta T\_{st}) \right) V/m$$

$$\text{where } A\_i \gg A\_{st}, \ \Delta T\_i \ll \Delta T\_{st}, \ i = 1, 2, \dots, N$$

Then the dispersion analysis-based optimization problem is stated as follows:

$$f\_p = \{f\_1, f\_2, \dots, f\_N\}$$

$$f\_{\min} = \left\{f\_{\min, 1}, f\_{\min, 2}, \dots, f\_{\min, N}\right\}, \ f\_{\max} = \left\{f\_{\max, 1}, f\_{\max, 2}, \dots, f\_{\max, N}\right\}$$

$$\begin{array}{c|c} \text{Given that } f\_{\min} &< f\_p \\ &< f\_{\max}; \max \left| \mathcal{E}\_1(f\_{st}) \\ &= f' ) \right| \text{ based on the following equations} \end{array}$$

*Appl. Sci.* **2020**, *10*, 1770

$$\nabla^2 \left( \mathbf{E}\_2 \left( \mathbf{f}\_p \right) \right) - \mu\_0 \varepsilon\_{\infty} \frac{\partial^2 \left( \mathbf{E}\_2 \left( \mathbf{f}\_p \right) \right)}{\partial t^2} = \mu\_0 \sigma \frac{\partial \left( \mathbf{E}\_2 \left( \mathbf{f}\_p \right) \right)}{\partial t} + \mu\_0 \frac{\partial^2 P\_2}{\partial t^2} \tag{9a}$$

$$\frac{\partial^2 P\_2}{\partial t^2} + \gamma \frac{\partial P\_2}{\partial t} + \omega\_0^2 (P\_2) - \frac{\omega\_0^2}{N \epsilon d} (P\_2)^2 - \frac{\omega\_0^2}{N^2 \epsilon^2 d^2} (P\_2)^3 = \frac{N \epsilon^2}{m} \left( E\_2 \{ f\_p \} \right) \tag{9b}$$

$$\nabla^2 \left( E\_1(f\_p) \right) - \mu\_0 \varepsilon\_{\infty} \frac{\partial^2 \left( E\_1(f\_p) \right)}{\partial t^2} = \mu\_0 \sigma \frac{\partial \left( E\_1(f\_p) \right)}{\partial t} + \mu\_0 \frac{\partial^2 (P\_1)}{\partial t^2} \tag{9c}$$

$$\begin{split} \frac{\partial^2 \langle P\_1 \rangle}{\partial t^2} + \gamma^\flat \frac{\partial \langle P\_1 \rangle}{\partial t} + a \nu^2 \langle P\_1 \rangle - \frac{a \nu\_0^2}{N \text{ad}} \Big[ P\_1 \,^2 + 2P\_1 P\_2 \Big] - \frac{a \nu\_0^2}{N^2 c^2 d^2} \Big[ P\_1 \,^3 + 3P\_1 \,^2 P\_2 + 3P\_1 P\_2 \Big] \\ = \frac{N c^2}{m} \Big( E\_1 \{ f\_p \} \Big) \end{split} \tag{9d}$$

Note that in this case the pump wave comprises *N* ultrashort pulses instead of a single ultrashort pulse. Therefore, we have a multiparameter optimization problem.

The physics behind the efficient amplification of the stimulus wave involves the simultaneous maximization of the stored electric energy density and the polarization density originated by the pump wave in the resonator (see Figure 3). This can be explained in two steps, first we need to maximize the stored energy in the cavity in order to transfer a high amount of energy to the stimulus wave. Second, we need to maximize the polarization density of the pump wave, which acts as the energy coupling coefficient based on Equation (7b).

**Figure 3.** Configuration of the cavity.

Even if we maximize the stored electric energy density in a resonator, if the non-linear coupling coefficient *P*<sup>2</sup> is not high, then we cannot efficiently transfer the accumulated energy into the stimulus wave and high-gain amplification of the stimulus wave does not occur.

In order to perform the optimization of the stimulus wave magnitude at a given stimulus wave frequency *max E*1(*fst* <sup>=</sup> *<sup>f</sup>* ) , we need an efficient optimization algorithm. In the next section, we will use the computationally efficient quasi-Newton Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm to perform the maximization of the stimulus wave magnitude at a desired frequency.
