**Boundary and excitation conditions:**

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

where *<sup>A</sup>*<sup>1</sup> = 1.5 <sup>×</sup> <sup>10</sup>8, *<sup>A</sup>*<sup>2</sup> = <sup>2</sup> <sup>×</sup> <sup>10</sup>8, *<sup>A</sup>*<sup>3</sup> = 2.5 <sup>×</sup> 108, <sup>Δ</sup>*T*<sup>1</sup> = <sup>4</sup>*ps*, <sup>Δ</sup>*T*<sup>2</sup> = <sup>2</sup>*ps*, <sup>Δ</sup>*T*<sup>2</sup> = <sup>1</sup>*ps*

$$\begin{array}{rcl} E\_{\rm st}(\mathbf{x} = 0 \mu m, t) & = 1 \times \frac{\sin \left( 2 \pi \left( 3 \times 10^{14} \right) t \right)}{\text{m}}, \text{ for } \; 0 \le t \\ & \le 40 \text{ps}, \; E\_{\rm lp}(\mathbf{x} = 15 \mu m, t) = E\_{\rm st}(\mathbf{x} = 15 \mu m, t) = 0 \qquad \text{for } 0 < t \\ & < 40 \text{ps} \end{array}$$

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

$$
\sigma(\mathbf{x}) = \left\{ (\mathbf{x} - (L - \Delta)) \sigma\_0 \, , \quad (L - \Delta) \le x < L \right\},
\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; Γ(*x* = 0μ*m*, *t*) = 1(*Re flection coe f ficient* = 1)

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

$$\left|\Gamma(f')\right| = 1 - e^{-\left(\frac{(f'-300THz)}{\sqrt{2}THz}\right)^2} \quad , \quad for \quad \ge = 10 \mu m, \quad 0 < t < 40ps$$

Augmented cost function to be maximized:

$$\begin{array}{c} \mathsf{Q}[f\_{p1}, f\_{p2}, f\_{p3}] &= \left| \mathbb{E}\_{\mathsf{sl}}(f\_{d} = 300THz) \right| \\ & \quad \quad \quad \quad \quad \quad \quad \begin{cases} \delta\_{\mathsf{sl}}(f\_{d} = 300THz) \\ \delta\_{\mathsf{l}}(f\_{\mathsf{p},i} - f\_{\mathsf{max}})^{2} + \delta\_{\mathsf{l}}2\left(f\_{\mathsf{min}} - f\_{\mathsf{p},i}\right)^{2} \end{cases} \\\\ \delta\_{\mathsf{i},1} = \left\{ \begin{array}{cc} 0 & \mathrm{if} & f\_{\mathsf{p},i} \leq f\_{\mathsf{max}} \\ \frac{\left| \mathbb{E}\_{\mathsf{sl}}(f\_{d} = 300THz) \right|}{10^{2}} & \mathrm{if} & f\_{\mathsf{p},i} > f\_{\mathsf{max}} \end{array} \right\}, \delta\_{\mathsf{i},2} = \left\{ \begin{array}{cc} 0 & \mathrm{if} & f\_{\mathsf{p},i} \geq f\_{\mathsf{min}} \\ \frac{\left| \mathbb{E}\_{\mathsf{sl}}(f\_{d} = 300THz) \right|}{10^{2}} & \mathrm{if} & f\_{\mathsf{p},i} < f\_{\mathsf{min}} \end{array} \right\}. \end{array}$$

**Optimization algorithm (BFGS):**

$$Set\ \ H\_1 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix},\ f\_{p1,0} = 200THz,\ f\_{p1,1} = 205THz$$

*fp*2,0 = 150*THz*, *fp*2,1 = 153*THz*, *fp*3,0 = 100*THz*, *fp*3,1 = 102*THz*, α<sup>1</sup> = 0.5

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

*p<sup>k</sup>* = − *Hk*∇*Qk*, *fp*, *<sup>k</sup>*+<sup>1</sup> = *fp*,*<sup>k</sup>* + α*kpk*, *fp*,*<sup>k</sup>* = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎣ *fp*1,*<sup>k</sup> fp*2,*<sup>k</sup> fp*3,*<sup>k</sup>* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎦ , *s<sup>k</sup>* = *fp*,*k*+1−*fp*,*<sup>k</sup>* ∇*Qk*+<sup>1</sup> = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *<sup>Q</sup>*(*fp*1,*k*+1, *fp*2,*k*, *fp*3,*k*)−*Q*(*fp*1,*k*, *fp*2,*k*, *fp*3,*k*) *fp*1,*k*+1−*fp*1,*<sup>k</sup> <sup>Q</sup>*(*fp*1,*k*, *fp*2,*k*+1, *fp*3,*k*)−*Q*(*fp*1,*k*, *fp*2,*k*, *fp*3,*k*) *fp*2,*k*+1−*fp*2,*<sup>k</sup> <sup>Q</sup>*(*fp*1,*k*, *fp*2, *fp*3,*k*+1)−*Q*(*fp*1,*k*, *fp*2,*k*, *fp*3,*k*) *fp*3,*k*+1−*fp*3,*<sup>k</sup>* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

$$\mathbf{y}\_k = \nabla \mathbf{Q}\_{k+1} - \nabla \mathbf{Q}\_k \quad \rho\_k = \frac{1}{\mathbf{y}\_k \mathbf{^T s}\_k}$$

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

$$\mathbf{I}: \mathbf{Id} \\mathbf{x}$$

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

$$\alpha\_k = \mathcal{C}^{(\lceil \log \rceil \frac{\mathbb{Q}(f\_{p,k})}{\mathbb{Q}(f\_{p,k}) - \mathbb{Q}(f\_{p,k-1})} \big| \ ) \ / \ (\! \lfloor \frac{\mathbb{Q}(f\_{p,k})}{\mathbb{Q}(f\_{p,k}) - \mathbb{Q}(f\_{p,k-1})} \big| \ ))$$

where *C* is just a constant (1 < *C* < 1.5) and α*<sup>k</sup>* is the step size at iteration k. In this simulation C = 1.45. Based on the above formulations, the maximum stimulus wave amplitude that has been reached in the cavity (for 0 < t < 40 ps) is determined as *Gainmax* = *Est*(*fst* <sup>=</sup> <sup>300</sup>*THz*) *max* <sup>=</sup> 6.34 <sup>×</sup> <sup>10</sup>8*V*/*m*, which corresponds to the frequencies *fp*<sup>1</sup> = 290.8*THz*, *fp*<sup>2</sup> = 410.6*THz*, *fp*<sup>3</sup> = 209.7*THz* (see Table 2), or to the free-space wavelengths λ*p*<sup>1</sup> = 1.032μm, λ*p*<sup>2</sup> = 0.731μ*m*, λ*p*<sup>3</sup> = 1.431μ*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} \textbf{Gain}\_{\text{max}} &= \left| \mathrm{E}\_{\text{st}} \left( f\_{\text{st}} = \text{300}THz \right) \right|\_{\text{max}} \\ &= \frac{\text{Amplitude of } f \text{ the stimulus move at } t \to t\_{\text{max}} \text{ for } x \coloneqq \textbf{x'}}{\text{Amplitude of } f \text{ the stimulus move at } t \to \textbf{0} \text{ for } x \coloneqq \textbf{x'}} \end{array}$$

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

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

$$P\_{pump}: \text{ Polarization density created by the pump wave} \left(\frac{Coulomb}{m^2}\right)$$

*Epump* : *Pump wave electric field intensity*, ε<sup>∞</sup> : *Background permittivity*


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

As we can see from Table 2, the optimal ultrashort pulse frequencies correspond to a very high stored electric energy density and a high polarization density. If we have a look at Table 2, we can see that 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. The polarization density serves as an energy-coupling coefficient and is a measure of the stored electric energy that can be coupled to the stimulus wave.
