High-Fidelity Harmonic Generation in Optical Micro-Resonators Using BFGS Algorithm

Abstract

Harmonic generation is an attractive research field that finds a variety of application areas. However, harmonic generation within a medium of micron-scale interaction length limits the magnitude of nonlinear coupling and leads to poor harmonic generation efficiency. In this study, we present a constrained non-linear programming approach based on the Quasi-Newton Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm to obtain high-fidelity harmonic generation in optical micro-resonators. Using this approach, one can achieve high-intensity harmonic generation in a simple Fabry–Perot type optical micro-resonator. The generation of super-intense harmonics at a typical ultraviolet (UV)-ablation frequency of 820 THz and at pure yellow-light (515 THz) is investigated in particular. Moreover, we achieved more than 98% accuracy compared to well-known theoretical results. Our approach enables the design of highly efficient microscale harmonic generators to be used in integrated photonic devices.

1. Introduction

Microscale harmonic generation is of high interest for many modern technologies, including optical antennas, lab-on-chip devices, and integrated optical devices. Currently, the conversion efficiency of microscale optical harmonic generators is very low due to the very small interaction length. Highly efficient targeted harmonic generation can be achieved with interaction mediums that are at least a few centimeters long. However, theoretically, the efficiency of generating a desired harmonic can be increased in the microscale by increasing the intensities of the waves that are intermixed, such an increase in intensity will definitely damage the micron-sized optical medium. Another solution may be to use an interaction material that has a superior nonlinearity, but materials with superior nonlinearity are usually not available. Even in the case of availability, super-nonlinear materials can be quite costly to employ in micron size as they are usually fabricated in thin film form and need to be arranged as an array of thin films to achieve the desired conversion efficiency. The discovery of novel super-nonlinear materials via experimentation is a hot-topic and some artificial materials do display an unusually high nonlinearity [1,2,3,4,5].

Although highly nonlinear natural materials and super-nonlinear artificial materials may help to increase the harmonic generation efficiency in the microscale, these materials usually lead to high dielectric, conduction, and scattering losses, along with being too expensive to afford or to fabricate. For these reasons, high-efficiency harmonic generation remains as an issue in the microscale. In this study, we used a computational approach to tackle this issue. A computational study that is based on maximizing the harmonic generation efficiency via a nonlinear programming algorithm is lacking [6,7,8]. Experimental studies usually focus on finding new techniques to increase the second and the third harmonic generation efficiency and there are relatively few experimental studies that have demonstrated a minor increase in the harmonic generation or conversion efficiency via certain experimental configurations and setups [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. Some computational studies have focused on increasing the computation efficiency of harmonic generation problems rather than proving that the harmonic generation efficiency itself can actually be increased. These include the studies mentioned in [6,7,8], which have managed to increase the efficiency of the nonlinear Finite-Difference-Time-Domain (FDTD) method by decreasing the dispersion error. The current literature of harmonic generation via nonlinear wave mixing is dominated by increasing the efficiency of the second/third harmonic generation rather than increasing the efficiency of an arbitrary harmonic [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. There are currently no known techniques in the literature that reported a high-efficiency in the microscale for an arbitrary harmonic generation (not necessarily the second harmonic) under monochromatic optical excitation [9,10,11,12,13,14,15,30,31,32,33,34,35,36,37,38,39]. In this study, we will show that for an arbitrary excitation frequency, one can maximize the harmonic generation efficiency at an arbitrary frequency. For example, if a microcavity is excited with an infra-red pump wave, and we want a super-high harmonic generation efficiency in the ultraviolet spectral range, we will show that this is possible by using constrained non-linear programming. The generation of a desired harmonic with an ultra-high efficiency is important for many applications in optics and biomedical engineering, such as for integrated optical devices that are used in spectroscopy and for Lab-On-Chip devices that are used for medical diagnostics. Another research field that can benefit from highly efficient harmonic generation at a desired frequency is the field of optical antennas. Optical antennas that can generate light at an arbitrary frequency would become ultra-wideband optical antennas through the embedding of a microscale controller. The microscale controller would have to do what we aim to do in this study, the adjustment of the micro-resonator parameters for maximizing the harmonic generation efficiency. Most importantly, ultra-efficient microscale harmonic generators can pave the way for macroscale high-efficiency harmonic generators that can generate powerful THz rays, ultraviolet (UV) rays, X-rays, and even gamma rays with much higher efficiencies than are currently available. High power yellow-light and ultraviolet-light generation is of particular interest as the power output of many yellow-light and ultraviolet-light generators is quite low. Especially high intensity UV-harmonic generation is an important issue as the harmonic conversion efficiency at this level is extremely low [40,41,42,43,44,45]. To provide a mathematical and computational proof that high-efficiency harmonic generation is feasible in a micro-resonator, we will present a full-fledged dispersion analysis. We will first investigate the nonlinear wave propagation concept in an arbitrary multi-resonant medium that is placed as an interaction medium in a plain Fabry–Perot type micro-resonator. Then, we will define the wave equation and its associated polarization density equations for each resonance frequency. Finally, we will define the pump wave as a combination of M ultrashort pulses and we will tune the pulses of these high-intensity pulses to maximize the frequency conversion efficiency at a target frequency that is to be generated. We will achieve this through an efficient constrained nonlinear programming algorithm that has a relatively lower computation cost and a relatively faster convergence rate such as the Quasi-Newton type Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. We will start our analysis by presenting the partial differential equations involved in nonlinear wave propagation in a multi-resonant interaction medium placed in a micro-resonator, and then we will present the formulations for the BFGS algorithm for the frequency tuning of the intense pulses that form the source wave which energizes the micro-resonator.

2. Nonlinear Wave Interaction in Optical Micro-Resonators

The interaction medium of an optical micro-resonator may have a single dominant resonance (emission) frequency or it may have multiple resonance frequencies associated with different resonance probabilities (oscillator strengths) [40,41,42,43,44,45]. Some materials, such as excitonic materials, have a single resonance frequency, while most materials have more than one resonance frequency. In this study, we will consider an arbitrary medium which has N resonance frequencies associated with N corresponding polarization damping rates. Quantum mechanics dictates that each electron that oscillates at a certain resonance frequency is associated with a resonance probability or oscillator strength and the sum of all resonance probabilities is equal to one. Assuming we have M different waves present in an arbitrary medium, if at least one of the waves has a sufficiently high intensity, then there will be a nonlinear coupling between the M waves and energy transfer between the waves will occur. For the total high-intensity wave that is present in a medium with N resonances, the wave equation that represents the time variation of the electric field $E$ of the high intensity wave and the associated equations that represent the components of the polarization density induced by the high intensity wave are given as

$${\nabla }^2\left(E\right)-{\mu }_0{\epsilon }_{\infty }\frac{{\partial }^2\left(E\right)}{\partial t^2}={\mu }_0\sigma \frac{\partial \left(E\right)}{\partial t}+{\mu }_0\frac{d^{2}P}{\partial t^2}$$

(1)

$$\frac{d^2{P}_1}{d{t}^2}+{\gamma }_1\frac{d{P}_1}{dt}+{\omega }_1{}^2{P}_1-\frac{{\omega }_1{}^2{P}_1{}^2}{Q_{1}ed}+\frac{{\omega }_1{}^2{P}_1{}^3}{Q_1{}^2{e}^2{d}^2}=\frac{Q_1{e}^{2}E}{m}$$

(2)

$$\frac{d^2{P}_2}{d{t}^2}+{\gamma }_2\frac{d{P}_2}{dt}+{\omega }_2{}^2{P}_2-\frac{{\omega }_2{}^2{P}_2{}^2}{Q_{2}ed}+\frac{{\omega }_2{}^2{P}_2{}^3}{Q_2{}^2{e}^2{d}^2}=\frac{Q_2{e}^{2}E}{m}$$

(3)

$$\frac{d^2{P}_N}{d{t}^2}+{\gamma }_N\frac{d{P}_N}{dt}+{\omega }_N{}^2{P}_N-\frac{{\omega }_N{}^2{P}_N{}^2}{Q_{N}ed}+\frac{{\omega }_N{}^2{P}_N{}^3}{Q_N{}^2{e}^2{d}^2}=\frac{Q_N{e}^{2}E}{m}$$

(4)

where P is the total polarization density, P_i is the polarization density component at the ith resonance frequency, γ_i is the polarization damping rate associated with the ith resonance frequency, ω_i is the ith angular resonance frequency, e is the electron charge, d is the atomic diameter, σ is the medium conductivity, and ${\epsilon }_{\infty }$ is the background permittivity. The number of electrons oscillating at each resonance is indicated by $Q_i$ and they are related to the number of electrons per unit volume (electron density) $Q$ via the oscillation strength ${\xi }_i$ such that

$${\displaystyle \sum }_{i=1}^N{\xi }_i=1, Q_i={\xi }_{i}Q, P={\displaystyle \sum }_{i=1}^N{P}_i={\displaystyle \sum }_{i=1}^N{Q}_i{p}_i=Q{\displaystyle \sum }_{i=1}^N{\xi }_i{p}_i$$

(5)

where, p_i is the dipole moment at the ith resonance.

Therefore, in order to determine the time variation of the electric field, we had to solve these N+1 equations. However, since our goal was to maximize the harmonic generation efficiency, we dealt with the nonlinear programming of these equations to maximize the spectral intensity at the desired frequency.

Assuming that we excite the optical microcavity with M high-intensity ultra-short pulses. The total wave that represents the overall excitation is stated as

$$E\left(x=0 \mathsf{\mu }\mathrm{m},t\right)={{\displaystyle \sum }}_{i=1}^M{A}_i\cos \left(2\pi {\nu }_{i}t+{\psi }_i\right)\left(u\left(t\right)-u\left(t-\Delta T_i\right)\right).$$

(6)

We wanted to tune the frequencies of the excitation, since there were N+1 equations and M optimization parameters involved in this configuration, we needed a cost-efficient optimization algorithm. The Quasi-Newton type Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm enabled a relatively fast convergence at a lower computational cost as the Hessin matrix was recursively updated and the computation of the second derivatives were not needed. Hence, we employed the BFGS algorithm for the optimization part.

Given the expression for excitation in Equation (6), and Equations (1)–(5) that modeled its propagation in the microcavity, we could then define the optimization problem. Note that this was a constrained optimization problem as the source device that generates the excitation could generate the ultrashort pulses at a certain frequency range and at a certain maximum pulse energy. Assuming that the we were solving the problem for a given pulse energy that was below the maximum available pulse energy, we had the allowable frequency range as the constraint of the problem. Most importantly, since we wanted a high-intensity output at a target frequency (desired harmonic), we had to define the cost function carefully and accurately. The intended description for the cost function can be stated as

$$F\left(\nu \right)=\left|E\left(\nu ={\nu }_{target}\right)\right|=\left|{{\displaystyle \int }}_{{\nu }_{target}-∆\nu }^{{\nu }_{target}+∆\nu }\left\{{{\displaystyle \int }}_0^{∆T}\left\{E\left(x=x^{\prime },t\right)e^{-i\left(2\pi \Omega \right)t}\right\}dt\right\}{e}^{i\left(2\pi \Omega \right)t}d\Omega \right|$$

(7)

where $∆\nu$ is the bandwidth around the target frequency. A certain target harmonic can be programmed to be maximized because the process of nonlinear wave mixing generates many new harmonics as a result of supercontinuum generation [20,21,22,23,24,25,26], which arise from high stored energy. What we intended to do is to program the microcavity to concentrate the supercontinuum spectral density around the target frequency. Therefore, the summary of the problem can be stated as follows: optimization parameters: $\nu =\left[{\nu }_1, {\nu }_2, \dots ,{\nu }_M\right]$, cost function to be maximized: $F\left(\nu \right)=\left|E\left(\nu ={\nu }_{target}\right)\right|$, constraints: ${\nu }_{min}\le \nu \le {\nu }_{max}$, equations:

$${\nabla }^2\left(E\left(\nu \right)\right)-{\mu }_0{\epsilon }_{\infty }\frac{{\partial }^2\left(E\left(\nu \right)\right)}{\partial t^2}={\mu }_0\sigma \frac{\partial \left(E\left(\nu \right)\right)}{\partial t}+{\mu }_0\frac{d^{2}P}{\partial t^2}$$

(8)

$$\frac{d^2{P}_1}{d{t}^2}+{\gamma }_1\frac{d{P}_1}{dt}+{\omega }_1{}^2{P}_1-\frac{{\omega }_1{}^2{P}_1{}^2}{N_{1}ed}+\frac{{\omega }_1{}^2{P}_1{}^3}{N_1{}^2{e}^2{d}^2}=\frac{N_1{e}^{2}E\left(\nu \right)}{m}$$

(9)

$$\frac{d^2{P}_2}{d{t}^2}+{\gamma }_2\frac{d{P}_2}{dt}+{\omega }_2{}^2{P}_2-\frac{{\omega }_2{}^2{P}_2{}^2}{N_{2}ed}+\frac{{\omega }_2{}^2{P}_2{}^3}{N_2{}^2{e}^2{d}^2}=\frac{N_2{e}^{2}E\left(\nu \right)}{m}$$

(10)

$$\frac{d^2{P}_N}{d{t}^2}+{\gamma }_N\frac{d{P}_N}{dt}+{\omega }_N{}^2{P}_N-\frac{{\omega }_N{}^2{P}_N{}^2}{N_{N}ed}+\frac{{\omega }_N{}^2{P}_N{}^3}{N_N{}^2{e}^2{d}^2}=\frac{N_N{e}^{2}E\left(\nu \right)}{m}$$

(11)

Note that even though the total excitation wave initially had M frequency components as seen in Figure 1, after the desired interaction duration, the spectrum of the excitation changed due to the nonlinear interaction. What we were trying to achieve here is to maximize the spectral density of the output spectrum around the target frequency, so that the output intensity at the desired harmonic was maximized.

Equations (1)–(4) can be discretized using the FDTD method at each iteration (k) of the optimization as follows:

$$\begin{array}{l}\frac{E_k\left(i+1,j\right)-2E_k\left(i,j\right)+E_k\left(i-1,j\right)}{∆x^2}-{\mu }_0{\epsilon }_{\infty }\left(i,j\right)\frac{E_k\left(i,j+1\right)-2E_k\left(i,j\right)+E_k\left(i,j-1\right)}{∆t^2}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\mu }_0\sigma \left(i,j\right)\frac{E_k\left(i,j\right)-E_k\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mu }_0\frac{P_k\left(i,j+1\right)-2P_k\left(i,j\right)+P_k\left(i,j-1\right)}{\Delta t^2}\end{array}$$

(12)

$$\begin{array}{l} \frac{P_{1,k}\left(i,j+1\right)-2P_{1,k}\left(i,j\right)+P_{1,k}\left(i,j-1\right)}{∆t^2}+{\gamma }_1\frac{P_{1,k}\left(i,j\right)-P_{1,k}\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+4{\pi }^2{f}_1{}^2\left(P_{1,k}\left(i,j\right)\right)-\frac{4{\pi }^2{f}_1{}^2}{Q_{1}ed}{\left(P_{1,k}\left(i,j\right)\right)}^2+\frac{4{\pi }^2{f}_1{}^2}{Q_1{}^2{e}^2{d}^2}{\left(P_{1,k}\left(i,j\right)\right)}^3\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{Q_1{e}^2}{m}\left(E_k\left(i,j\right)\right) \end{array}$$

(13)

$$\begin{array}{l} \frac{P_{2,k}\left(i,j+1\right)-2P_{2,k}\left(i,j\right)+P_{2,k}\left(i,j-1\right)}{∆t^2}+{\gamma }_2\frac{P_{2,k}\left(i,j\right)-P_{2,k}\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+4{\pi }^2{f}_2{}^2\left(P_{2,k}\left(i,j\right)\right)-\frac{4{\pi }^2{f}_2{}^2}{Q_{2}ed}{\left(P_{2,k}\left(i,j\right)\right)}^2+\frac{4{\pi }^2{f}_2{}^2}{Q_2{}^2{e}^2{d}^2}{\left(P_{2,k}\left(i,j\right)\right)}^3\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{Q_2{e}^2}{m}\left(E_k\left(i,j\right)\right) \end{array}$$

(14)

$$\begin{array}{l} \frac{P_{N,k}\left(i,j+1\right)-2P_{N,k}\left(i,j\right)+P_{N,k}\left(i,j-1\right)}{∆t^2}+{\gamma }_N\frac{P_{N,k}\left(i,j\right)-P_{N,k}\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+4{\pi }^2{f}_N{}^2\left(P_{N,k}\left(i,j\right)\right)-\frac{4{\pi }^2{f}_N{}^2}{Q_{N}ed}{\left(P_{N,k}\left(i,j\right)\right)}^2+\frac{4{\pi }^2{f}_N{}^2}{Q_N{}^2{e}^2{d}^2}{\left(P_{N,k}\left(i,j\right)\right)}^3\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{Q_N{e}^2}{m}\left(E_k\left(i,j\right)\right)\end{array}$$

(15)

$$P_k=P_{1,k}+P_{2,k}+\cdots +P_{N,k}, Q=Q_1+Q_2+\cdots +Q_N$$

(16)

Equations (12)–(15) were solved for each new update. For high precision, we chose $∆t$ and $∆x$ to be small.

3. Non-Linear Programming for Efficient Harmonic Generation

As already mentioned, we used the BFGS algorithm for its relatively low computational cost and high convergence rate. The BFGS algorithm is a Quasi-Newton algorithm that recursively computes the Hessian matrix instead of computing the second derivative of the cost function at every iteration. The BFGS algorithm is useful when the Hessian matrix is not available or when it is too costly to compute.

3.1. BFGS Algorithm-Based Optimization

Start with an initial estimate of the Hessian matrix

$${\mathit{H}}_0=\mathit{I}\left(\mathit{I}:\mathrm{Identity} \mathrm{matrix}\right).$$

(17)

Identify the cost function that includes the penalty terms for constraint violations

$$F\left(\nu \right)=\left|E\left(\nu ={\nu }_{target}\right)\right|=$$

(18)

$$\left|{{\displaystyle \int }}_{{\nu }_{target}-∆\nu }^{{\nu }_{target}+∆\nu }\left\{{{\displaystyle \int }}_0^{∆T}\left\{E\left(x=x^{\prime },t\right)e^{-i\left(2\pi \Omega \right)t}\right\}dt\right\}{e}^{i\left(2\pi \Omega \right)t}d\Omega \right|-L_1{\left\{{\displaystyle \sum }_{i=1}^N{\delta }_i\left(v_i-v_{max}\right)\right\}}^q-L_2{\left\{{\displaystyle \sum }_{i=1}^N{\zeta }_i\left(v_{min}-v_i\right)\right\}}^q$$

$${\delta }_i=\left\{\begin{array}{c}0 if v_i\le v_{max}\\ >0 if v_i>v_{max}\end{array}\right\}, {\zeta }_i=\left\{\begin{array}{c}0 if v_i\ge v_{min}\\ >0 if v_i<v_{min}\end{array}\right\}$$

q: exponent of the penalty (positive valued), $\left\{L_1,L_2\right\}$: penalty constants (positive valued), $\left\{{\delta }_i,{\zeta }_i\right\}$: penalty weights

3.2. Determining the Penalty Weights

In order to decrease $\left|E\left(\nu ={\nu }_{target}\right)\right|$ by a factor of $\left(100-\frac{100}{\mathsf{\sigma }}\right)\%$ for a deviation of $∆\mathsf{\nu }$ from ${\mathsf{\nu }}_{\max }$ or ${\mathsf{\nu }}_{\min }$, the penalty coefficients ${\delta }_1$ and ${\delta }_2$ were chosen as (assuming single-pulse tuning)

$${\delta }_1=\left\{\begin{array}{ll}0& if {\nu }_p<{\nu }_{max}\\ \frac{\left(1-\frac{1}{\sigma }\right)}{{\left(∆\nu \right)}^2}\left|E\left({\nu }_p\right)\right|& if {\nu }_p>{\nu }_{max}\end{array}\right\}, {\delta }_2=\left\{\begin{array}{ll}0& if {\nu }_p>{\nu }_{min}\\ \frac{\left(1-\frac{1}{\sigma }\right)}{{\left(∆\nu \right)}^2}\left|E\left({\nu }_p\right)\right|& if {\nu }_p<{\nu }_{min}\end{array}\right\}$$

(19)

$∆\nu$: deviation from max/min allowable frequency, $\sigma$: reduction factor.

$$\mathrm{Such} \mathrm{that} F=\left|E\left(\nu ={\nu }_{target}\right)\right|-{\delta }_1{\left(\nu -{\nu }_{max}\right)}^2-{\delta }_2{\left({\nu }_{min}-\nu \right)}^2$$

(20)

Identify the search direction, ${\mathit{p}}_k=-$ ${\mathit{H}}_k$ $\mathit{\nabla }{F}_k$ (Search direction). Determine a suitable step size (Backtracking approach). Select $\alpha >0,\rho ϵ \left(0,1\right),c ϵ \left(0,1\right).$ while $F(x_k+$ $\alpha p_k)\le$ $F\left(x_k\right)+$ $c\alpha \nabla F_k{}^T{p}_k \alpha \leftarrow \rho \alpha$ end.

Update the values of the optimization parameters: ${\mathit{\nu }}_{\mathit{k}+1}={\mathit{\nu }}_{\mathit{k}}+{\mathit{\alpha }}_{\mathit{k}}$ ${\mathit{p}}_{\mathit{k}}$ (Update equation).

Compute the parameter difference vector and the gradient difference vector ${\mathit{s}}_{\mathit{k}}={\mathit{\nu }}_{\mathit{k}+1}-$ ${\mathit{\nu }}_{\mathit{k}}$, ${\mathit{y}}_{\mathit{k}}=$ $\mathit{\nabla }{\mathit{F}}_{\mathit{k}+1}-$ $\mathit{\nabla }{\mathit{F}}_{\mathit{k}}$. Update the Hessian matrix for the next iteration ${\mathit{H}}_{\mathit{k}+1}=\left(\mathit{I}-{\rho }_k{\mathit{s}}_{\mathit{k}}{\mathit{y}}_{\mathit{k}}{}^{\mathit{T}}\right){\mathit{H}}_{\mathit{k}}\left(\mathit{I}-{\rho }_k{\mathit{y}}_{\mathit{k}}{\mathit{s}}_{\mathit{k}}{}^{\mathit{T}}\right)+{\rho }_k{\mathit{s}}_{\mathit{k}}{\mathit{s}}_{\mathit{k}}{}^{\mathit{T}}$ (BFGS update)

$$\mathit{\nabla }F=\left[\begin{array}{c}\frac{F\left({\nu }_1+ϵ,{\nu }_2,\dots ,{\nu }_N\right)-F\left({\nu }_1,{\nu }_2,\dots ,{\nu }_N\right)}{ϵ}\\ \frac{F\left({\nu }_1,{\nu }_2+ϵ,\dots ,{\nu }_N\right)-F\left({\nu }_1,{\nu }_2,\dots ,{\nu }_N\right)}{ϵ}\\ .\\ .\\ .\\ \frac{F\left({\nu }_1,{\nu }_2,\dots ,{\nu }_N+ϵ\right)-F\left({\nu }_1,{\nu }_2,\dots ,{\nu }_N\right)}{ϵ}\end{array}\right] , {\rho }_k=\frac{1}{{\mathit{y}}_{\mathit{k}}{}^{\mathit{T}}{\mathit{s}}_{\mathit{k}}}$$

(21)

Note that we solved the FDTD equations at every iteration of the BFGS algorithm until the desired intensity at a target frequency was obtained. The summary of the whole process is summarized in the flowchart presented in Figure 2.

4. Numerical Simulations

4.1. Simulation1: Intense Harmonic Generation in the Ultraviolet (UV) Frequency Range

The total excitation (pump) wave $E$ that was composed of two high-intensity ultrashort pulses, energized a Fabry–Perot type optical micro-resonator that had an optical isolator at the input (left) port and a band-pass filter at the output (right) port. These two ultrashort pulses were both initiated at x = 2.5 μm at time t = 0 sec (Figure 3). The total excitation wave at the input port can be expressed as

$$\begin{array}{c}E\left(x=0 \mathsf{\mu }\mathrm{m},t\right)={\displaystyle {\displaystyle \sum }_{i=1}^2}{A}_i\cos \left(2\pi {\nu }_{i}t+{\psi }_i\right)(u\left(t\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-u\left(t-\Delta T_i\right))\frac{V}{m}, u\left(t\right):Unit step function\end{array}$$

(22)

$$A_1=8\times {10}^7, A_2=1\times {10}^8, \Delta T_1=1.5 \mathrm{ps}, \Delta T_2=1 \mathrm{ps}$$

Our goal was to generate a desired monochromatic wave (harmonic) at f = 820 THz, a typical UV-ablation frequency. To achieve this, the excitation frequencies of the ultrashort pulses were tuned.

• $\mathrm{Spatial} \mathrm{and} \mathrm{temporal} \mathrm{intervals} \mathrm{of} \mathrm{the} \mathrm{simulation}: 0\le x\le 10 \mathsf{\mu }\mathrm{m}, 0\le t\le 10 \mathrm{ps}$
• $\mathrm{Resonance} \mathrm{frequencies} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: f_r=\left\{4\times {10}^{14} \mathrm{Hz}, 6.3\times {10}^{14} \mathrm{Hz}, 8.8\times {10}^{14} \mathrm{Hz}\right\}$
• $\mathrm{Polarization} \mathrm{damping} \mathrm{rates} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: \gamma =\left\{2\times {10}^9 \mathrm{Hz}, 3\times {10}^{10} \mathrm{Hz}, 1\times {10}^{11} \mathrm{Hz}\right\}$
• $\mathrm{Relative} \mathrm{permittivity} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: \left({\epsilon }_r\right)=10 ({\mu }_r=1)$
• $\mathrm{Location} \mathrm{of} \mathrm{the} \mathrm{optical} \mathrm{isolator}: x=0 \mathsf{\mu }\mathrm{m}, \mathrm{Bandpass} \mathrm{filter} \mathrm{location}: x=10 \mathsf{\mu }\mathrm{m}$
• $\mathrm{Spatial} \mathrm{range} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{material} :0 \mathsf{\mu }\mathrm{m}<x<10 \mathsf{\mu }\mathrm{m}$$\mathrm{Density} \mathrm{of} \mathrm{electrons}: Q=3.5\times \frac{{10}^{28}}{m^3}, \mathrm{Atom} \mathrm{diameter}: d=0.3 \mathrm{nanometers}$
• Resonance probabilities: $\xi =\left\{0.3,0.4,0.3\right\}, \mathrm{Cos}\mathrm{t} \mathrm{function} \mathrm{to} \mathrm{be} \mathrm{maximized}: C$

Problem statement: Identifying the excitation frequencies of the high-intensity ultrashort pulses for maximizing the intensity spectral density around the target frequency ($\left|E\left(\nu ={\nu }_{target}=820\mathrm{THz}\right)\right|$) inside the micro-resonator, for 50 THz < $\left\{{\nu }_1,{\nu }_2\right\}<500 \mathrm{THz}$, and for $0 \mathsf{\mu }\mathrm{m}<x<10 \mathsf{\mu }\mathrm{m}, 0\le \mathrm{t}\le 10 \mathrm{ps}$.

$$\begin{array}{c}C=\left|E\left(\nu ={\nu }_{target}=820\mathrm{THz}\right)\right|\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\left|{{\displaystyle \int }}_{8.2\times {10}^{14}-∆\nu }^{8.2\times {10}^{14}+∆\nu }\left\{{{\displaystyle \int }}_0^{∆T}\left\{E_{in}\left(x=x^{\prime },t\right)e^{-i\left(2\pi \Omega \right)t}\right\}dt\right\}{e}^{i\left(2\pi \Omega \right)t}d\Omega \right|\end{array}$$

(23)

$$∆T=10 \mathrm{ps}, \left(8.2\times {10}^{14}-∆\nu \right) \mathrm{Hz}<\Omega <\left(8.2\times {10}^{14}+∆\nu \right) \mathrm{Hz}, ∆\nu =10 \mathrm{THz}$$

Initial conditions of the electric field and polarization density: (Prime sign refers to the time derivative)

$$\begin{array}{c}{P}_2\left(x,0\right)=P_2{}^{’}\left(x,0\right)=E_2\left(x,0\right)=E_2{}^{’}\left(x,0\right)=P_1\left(x,0\right)=P_1{}^{’}\left(x,0\right)=E_1\left(x,0\right)\\ =E_1{}^{’}\left(x,0\right)=0\end{array}$$

(24)

Bandpass filtering: Frequency selective cavity wall (right port) is fixed at $\mathrm{x}=10 \mathsf{\mu }\mathrm{m}$

$$\begin{array}{l}\left|\Gamma \left(\nu \right)\right|\\ =1-e^{-{(\frac{\left(\nu -820\mathrm{THz}\right)}{\sqrt{200\mathrm{THz}}})}^2} \left(\mathrm{Magnitude} \mathrm{response} \mathrm{of} \mathrm{the} \mathrm{harmonic} \mathrm{selective} \mathrm{output} \mathrm{wall}\right)\end{array}$$

(25)

Cost function: $C\left({\nu }_1,{\nu }_2\right)=\left|E\left(\nu =820 \mathrm{THz}\right)\right|-{{\displaystyle \sum }}_{i=1}^2\{{\delta }_{i,1}{\left({\nu }_i-{\nu }_{\max }\right)}^2+{\delta }_{i,2}{\left({\nu }_{\min }-{\nu }_i\right)}^2\}$

$$\begin{array}{c}{\delta }_{i,1}=\left\{\begin{array}{ll}0& if {\nu }_i\le {\nu }_{\max }\\ \frac{\left|E\left(\nu =820 \mathrm{THz}\right)\right|}{{10}^{27}}& if {\nu }_i>{\nu }_{\max }\end{array}\right\},\\ {\delta }_{i,2}=\left\{\begin{array}{ll}0& if {\nu }_i\ge {\nu }_{\min }\\ \frac{\left|E\left(\nu =820 \mathrm{THz}\right)\right|}{{10}^{27}}& if {\nu }_i<{\nu }_{\min }\end{array}\right\}\end{array}$$

(26)

FDTD Equations: Discretization of Equations (1)–(4) via finite difference time domain is as follows [27,28,29,30]

$$\begin{array}{l}\frac{E_k\left(i+1,j\right)-2E_k\left(i,j\right)+E_k\left(i-1,j\right)}{∆x^2}-{\mu }_0{\epsilon }_{\infty }\left(i,j\right)\frac{E_k\left(i,j+1\right)-2E_k\left(i,j\right)+E_k\left(i,j-1\right)}{∆t^2}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\mu }_0\sigma \left(i,j\right)\frac{E_k\left(i,j\right)-E_k\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mu }_0\frac{P_k\left(i,j+1\right)-2P_k\left(i,j\right)+P_k\left(i,j-1\right)}{\Delta t^2}\end{array}$$

(27)

$$\begin{array}{l}\frac{P_{1,k}\left(i,j+1\right)-2P_{1,k}\left(i,j\right)+P_{1,k}\left(i,j-1\right)}{∆t^2}+{\gamma }_1\frac{P_{1,k}\left(i,j\right)-P_{1,k}\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+4{\pi }^2{f}_1{}^2\left(P_{1,k}\left(i,j\right)\right)-\frac{4{\pi }^2{f}_1{}^2}{Q_{1}ed}{\left(P_{1,k}\left(i,j\right)\right)}^2+\frac{4{\pi }^2{f}_1{}^2}{Q_1{}^2{e}^2{d}^2}{\left(P_{1,k}\left(i,j\right)\right)}^3\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{Q_1{e}^2}{m}\left(E_k\left(i,j\right)\right).\end{array}$$

(28)

$$\begin{array}{l}\frac{P_{2,k}\left(i,j+1\right)-2P_{2,k}\left(i,j\right)+P_{2,k}\left(i,j-1\right)}{∆t^2}+{\gamma }_2\frac{P_{2,k}\left(i,j\right)-P_{2,k}\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+4{\pi }^2{f}_2{}^2\left(P_{2,k}\left(i,j\right)\right)-\frac{4{\pi }^2{f}_2{}^2}{Q_{2}ed}{\left(P_{2,k}\left(i,j\right)\right)}^2+\frac{4{\pi }^2{f}_2{}^2}{Q_2{}^2{e}^2{d}^2}{\left(P_{2,k}\left(i,j\right)\right)}^3\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{Q_2{e}^2}{m}\left(E_k\left(i,j\right)\right).\end{array}$$

(29)

$$\begin{array}{l}\frac{P_{3,k}\left(i,j+1\right)-2P_{3,k}\left(i,j\right)+P_{3,k}\left(i,j-1\right)}{∆t^2}+{\gamma }_3\frac{P_{3,k}\left(i,j\right)-P_{3,k}\left(i,j-1\right)}{∆t}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+4{\pi }^2{f}_3{}^2\left(P_{3,k}\left(i,j\right)\right)-\frac{4{\pi }^2{f}_3{}^2}{Q_{3}ed}{\left(P_{3,k}\left(i,j\right)\right)}^2+\frac{4{\pi }^2{f}_3{}^2}{Q_3{}^2{e}^2{d}^2}{\left(P_{3,k}\left(i,j\right)\right)}^3\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{Q_3{e}^2}{m}\left(E_k\left(i,j\right)\right).\end{array}$$

(30)

• $P=$$P_1+P_2+P_3$
• x: Space coordinate, t: Time, k: Iteration, $E_k\left(x,t\right)=E_k\left(i∆x,j∆t\right) \to E_k\left(i,j\right)$

$$E_k: \mathrm{Total} \mathrm{wave} \mathrm{electric} \mathrm{field} \mathrm{at} \mathrm{the} \mathrm{kth} \mathrm{update}$$

Optimization via BFGS algorithm: Choose the identity matrix as the initial Hessian matrix

$${\mathit{H}}_0=\mathit{I} \left(\mathit{I}:2\times 2 \mathrm{identity} \mathrm{matrix}\right)$$

$$v_{p1,0}=250 \mathrm{THz}, v_{p1,1}=245 \mathrm{THz}, v_{p2,0}=225 \mathrm{THz}, v_{p2,1}=220 \mathrm{THz}, {\alpha }_1=0.4$$

$$\mathit{\nabla }{\mathit{C}}_k=\left[\begin{array}{c}\frac{C\left(f_{p1,k}, f_{p2,k}\right)-C\left(f_{p1,k-1}, f_{p2,k}\right)}{f_{p1,k}-f_{p1,k-1}}\\ \frac{C\left(f_{p1,k}, f_{p2,k}\right)-C\left(f_{p1,k}, f_{p2,k-1}\right)}{f_{p2,k}-f_{p2,k-1}}\end{array}\right]$$

(31)

• ${\mathit{p}}_k=-$${\mathit{H}}_k$$\mathit{\nabla }{\mathit{C}}_k$, ${\mathit{\nu }}_{ \mathit{k}+1}={\mathit{\nu }}_{\mathit{k}}+{\mathit{\alpha }}_{\mathit{k}}$${\mathit{p}}_k$, ${\mathit{s}}_k={\mathit{v}}_{k+1}-$${\mathit{v}}_k,$${\mathit{\nu }}_{\mathit{k}}=\left[\begin{array}{c}{\nu }_{1,k}\\ {\nu }_{2,k}\end{array}\right]$
• ${\mathit{y}}_k=$$\mathit{\nabla }{\mathit{C}}_{k+1}-$$\mathit{\nabla }{\mathit{C}}_k$${\rho }_k=\frac{1}{y_k{}^T{s}_k}$

BFGS recursion: ${\mathit{H}}_{k+1}=\left(\mathit{I}-{\rho }_k{\mathit{s}}_k{\mathit{y}}_k{}^T\right){\mathit{H}}_k\left(\mathit{I}-{\rho }_k{\mathit{y}}_k{\mathit{s}}_k{}^T\right)+{\rho }_k{\mathit{s}}_k{\mathit{s}}_k{}^T$$I:Identity matrix$

$$\mathit{\nabla }{\mathit{C}}_{k+1}=\left[\begin{array}{c}\frac{C\left(f_{p1,k+1}, f_{p2,k}\right)-C\left(f_{p1,k}, f_{p2,k}\right)}{f_{p1,k+1}-f_{p1,k}}\\ \frac{C\left(f_{p1,k}, f_{p2,k+1}\right)-C\left(f_{p1,k}, f_{p2,k}\right)}{f_{p2,k+1}-f_{p2,k}}\end{array}\right]$$

(32)

A simple formula for the step size that complies with the backtracking approach is determined as

$${\alpha }_k=c^{( \log \left|\frac{C\left({\nu }_k\right)}{C\left({\nu }_k\right)-C\left({\nu }_{k-1}\right)}\right| ) / ( \left|\frac{C\left({\nu }_k\right)}{C\left({\nu }_k\right)-C\left({\nu }_{k-1}\right)}\right| )}$$

(33)

The base of the step size ($c)$ was a simple constant ($1<c<1.5$) and ${\alpha }_k$ was the step size at iteration k. In this example, c was chosen as c = 1.45. According to these presented formulations, the maximum amplitude of the desired harmonic that was attained in the micro-resonator (for 0 < t < 10 ps) was obtained as $\left|E\left(\nu ={\nu }_{target}=820 \mathrm{THz}\right)\right|=7.6\times {10}^7 \mathrm{V}/\mathrm{m}$, and the resulting optimal excitation frequencies were identified as ${\nu }_1=273.2 \mathrm{THz}, {\nu }_2=284.7 \mathrm{THz}$(see

Table 1. Broyden–Fletcher–Goldfarb–Shanno (BFGS) updating process.

$\left|{\mathit{E}}_{\mathit{\nu }=820 \mathbf{THz}}\right|$	${\mathit{\nu }}_1$	${\mathit{\nu }}_2$	${\mathit{W}}_{\mathit{e},\mathit{p}} \left(\frac{\mathit{J}}{{\mathit{m}}^3}\right)$	${\mathit{P}}_{\mathit{p}\mathit{u}\mathit{m}\mathit{p}}\left(\frac{\mathit{C}}{{\mathit{m}}^2}\right)$	k (Iteration #)
5.8 $\times {10}^3 \mathrm{V}/\mathrm{m}$	250 THz	225 THz	2.9 $\times {10}^7$	0.09	1
2.9 $\times {10}^3 \mathrm{V}/\mathrm{m}$	245 THz	220 THz	1 $\times {10}^7$	0.05	2
2.7 $\times {10}^5 \mathrm{V}/\mathrm{m}$	249.4 THz	226.1 THz	8.66 $\times {10}^7$	0.14	3
2.0 $\times {10}^4 \mathrm{V}/\mathrm{m}$	257.2 THz	231.3 THz	5.49 $\times {10}^7$	0.13	4
1.3 $\times {10}^3 \mathrm{V}/\mathrm{m}$	264.9 THz	226.6 THz	1.31 $\times {10}^7$	0.06	5
7.6 $\times {10}^3 \mathrm{V}/\mathrm{m}$	276.0 THz	222.2 THz	1.76 $\times {10}^7$	0.07	6
1.4 $\times {10}^4 \mathrm{V}/\mathrm{m}$	286.7 THz	231.8 THz	1.73 $\times {10}^7$	0.07	7
1.1 $\times {10}^6 \mathrm{V}/\mathrm{m}$	279.7 THz	218.0 THz	3.06 $\times {10}^7$	0.10	8
9.8 $\times {10}^5 \mathrm{V}/\mathrm{m}$	272.1 THz	236.9 THz	2.49 $\times {10}^7$	0.10	9
5.1 $\times {10}^6 \mathrm{V}/\mathrm{m}$	263.5 THz	285.1 THz	3.30 $\times {10}^7$	0.10	10
1.2 $\times {10}^7 \mathrm{V}/\mathrm{m}$	273.4 THz	310.7 THz	8.21 $\times {10}^7$	0.16	11
3.2 $\times {10}^6 \mathrm{V}/\mathrm{m}$	275.1 THz	288.8 THz	3.48 $\times {10}^7$	0.12	12
7.6 $\times {10}^7 \mathrm{V}/\mathrm{m}$	273.2 THz	284.7 THz	5.75 $\times {10}^7$	0.17	13

).

$$\begin{array}{l}{W}_{e,p}=\mathrm{Stored} \mathrm{electric} \mathrm{energy} \mathrm{density}=\frac{1}{2}{\epsilon }_{\infty }{E}^2+\frac{1}{2}EP \left(\frac{Joules}{m^3}\right),\\ E :\mathrm{Electric} \mathrm{field} \mathrm{intensity}\end{array}$$

(34)

$P: \mathrm{Polarization} \mathrm{density} \mathrm{created} \mathrm{by} \mathrm{the} \mathrm{pump} \mathrm{wave}$$\left(\frac{Coulomb}{m^2}\right), {\epsilon }_{\infty }:\mathrm{Background} \mathrm{permittivity}.$

As indicated in Table 1, the optimal frequencies of the ultrashort excitation pulses gave rise to an enormous amount of stored electric energy and a corresponding large amount of stored polarization density (non-linear coupling coefficient). When we investigated Table 1, we realized that the induced electric energy density and the induced polarization density should be concurrently high for a super-intense harmonic generation at a desired frequency. The amount of stored electric energy was crucial in generating the high-intensity target harmonic and the induced polarization density was important in concentrating the available spectral energy around the desired frequency and in transferring the energy from one harmonic to another. It is important to note that in a multi-resonant interaction medium with many resonances, if the desired frequency (target harmonic) is near one of the resonance frequencies of the medium, generating the target harmonic with a high intensity is harder as the dielectric absorption is stronger near the resonance frequencies of the medium. Therefore, if the desired frequency of the harmonic is near any resonance, one should expect a relatively lower-intensity target harmonic after the optimization process is completed. To provide an increase in the intensity of the desired harmonic, one may use more ultrashort excitation pulses. At the end of 13th iteration, electric field at the output of the bandpass filter for 10 ps simulation time is seen in Figure 4. As seen in this figure, electric field at 820 THz reaches in a scale of 10⁸ V/m.

In Figure 5 we present spectral magnitude through the iterations. As seen in Figure 5a there is no generated target frequency (820 THz) at 6th iteration. When the iterations proceed (Figure 5b–d), wave mixing between two source waves happens, and we obtained generated target frequency at 820 THz. At the end of 13th iteration, intensity spectral density at the band-pass filter output for 822 THz is 1023 (W/m²Hz) as seen in Figure 6.

4.2. Simulation 2: Intense Quasi-Monochromatic Yellow-Light Generation Around 515THz

The total excitation (pump) wave $E$ that was composed of two high-intensity ultrashort pulses, energized a Fabry–Perot type optical micro-resonator that had an optical isolator at the input (left) port and a band-pass filter at the output (right) port as seen in Figure 7. All ultrashort pulses were initiated at x = 2.5 μm at time t = 0 sec. The total excitation wave at the input port can be expressed as

$$E\left(x=0 \mathsf{\mu }\mathrm{m},t\right)={\displaystyle \sum }_{i=1}^2{A}_i\cos \left(2\pi {\nu }_{i}t+{\psi }_i\right)\left(u\left(t\right)-u\left(t-\Delta T_i\right)\right) \mathrm{V}/\mathrm{m}$$

(35)

$$A_1=1\times {10}^8, A_2=1.2\times {10}^8, \Delta T_1=1\mathrm{ps}, \Delta T_2=0.7\mathrm{ps}$$

Our goal was to generate a desired monochromatic wave (harmonic) at f = 515 THz. To achieve this, the excitation frequencies of the ultrashort pulses were tuned.

• $\mathrm{Spatial} \mathrm{and} \mathrm{temporal} \mathrm{intervals} \mathrm{of} \mathrm{the} \mathrm{simulation}: 0\le x\le 10 \mathsf{\mu }\mathrm{m}, 0\le t\le 10\mathrm{ps}$
• $\mathrm{Resonance} \mathrm{frequencies} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: f_r=\left\{3\times {10}^{14}\mathrm{Hz}, 4.4\times {10}^{14}\mathrm{Hz}, 6.3\times {10}^{14}\mathrm{Hz}\right\}$
• $\mathrm{Polarization} \mathrm{damping} \mathrm{rates} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: \gamma =\left\{1\times {10}^{10}\mathrm{Hz}, 2.5\times {10}^{10}\mathrm{Hz}, 1\times {10}^{11}\mathrm{Hz}\right\}$
• Resonance probabilities: $\xi =\left\{\frac{1}{3},\frac{1}{3},\frac{1}{3}\right\}$, $\mathrm{Permittivity} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: \left({\epsilon }_r\right)=12 ({\mu }_r=1)$
• $\mathrm{Location} \mathrm{of} \mathrm{the} \mathrm{optical} \mathrm{isolator} \left(\mathrm{input} \mathrm{port}\right): x=0 \mathsf{\mu }\mathrm{m}, \mathrm{Filter} \left(\mathrm{output} \mathrm{port}\right) \mathrm{location}: x=10 \mathsf{\mu }\mathrm{m}$
• $\mathrm{Spatial} \mathrm{range} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{material} :0 \mathsf{\mu }\mathrm{m}<x<10 \mathsf{\mu }\mathrm{m}$, $\mathrm{Density} \mathrm{of} \mathrm{electrons}: N=3.5\times \frac{{10}^{28}}{m^3}$
• $\mathrm{Atom} \mathrm{diameter}: d=0.3 \mathrm{nanometers}, \mathrm{Cos}\mathrm{t} \mathrm{function} \mathrm{to} \mathrm{be} \mathrm{maximized}: C$

Problem statement: Identifying the excitation frequencies of the high-intensity ultrashort pulses for maximizing the intensity spectral density around a target frequency ($\left|E\left(\nu ={\nu }_{target}=515 \mathrm{THz}\right)\right|$) inside the micro-resonator, for 50THz < $\left\{{\nu }_1,{\nu }_2\right\}<500 \mathrm{THz}$, and for $0 \mathsf{\mu }\mathrm{m}<x<10 \mathsf{\mu }\mathrm{m}, 0\le t\le 10 \mathrm{ps}$.

$$\begin{array}{c}C=\left|E\left(\nu ={\nu }_{target}=515THz\right)\right|\hfill \\ \hfill =\left|{{\displaystyle \int }}_{5.15\times {10}^{14}-∆\nu }^{5.15\times {10}^{14}+∆\nu }\left\{{{\displaystyle \int }}_0^{∆T}\left\{E_{\mathrm{in}}\left(x=x^{\prime },t\right)e^{-i\left(2\pi \Omega \right)t}\right\}dt\right\}{e}^{i\left(2\pi \Omega \right)t}d\Omega \right|\end{array}$$

(36)

$$∆T=10\mathrm{ps}, \left(5.15\times {10}^{14}-∆\nu \right)\mathrm{Hz}<\Omega <\left(5.15\times {10}^{14}+∆\nu \right)\mathrm{Hz}, ∆\nu =10\mathrm{THz}$$

Initial conditions of the electric field and polarization density: (Prime sign refers to the time derivative)

$$P_2\left(x,0\right)=P_2{}^{\prime }\left(x,0\right)=E_2\left(x,0\right)=E_2{}^{\prime }\left(x,0\right)=P_1\left(x,0\right)=P_1{}^{\prime }\left(x,0\right)=E_1\left(x,0\right)=E_1{}^{\prime }\left(x,0\right)=0$$

Band-pass filtering: Frequency selective cavity wall (right port) is fixed at $\mathrm{x}=10\mathsf{\mu }\mathrm{m}$

$$\left|\Gamma \left(\nu \right)\right|=1-e^{-{(\frac{\left(\nu -515\mathrm{THz}\right)}{\sqrt{200}\mathrm{THz}})}^2}$$

(37)

Cost function:$C\left({\nu }_1,{\nu }_2\right)=\left|E\left(\nu =515\mathrm{THz}\right)\right|-{{\displaystyle \sum }}_{i=1}^2\{{\delta }_{i,1}{\left({\nu }_i-{\nu }_{max}\right)}^2+{\delta }_{i,2}{\left({\nu }_{min}-{\nu }_i\right)}^2\}$

$$\begin{array}{c}{\delta }_{i,1}=\left\{\begin{array}{ll}0& if {\nu }_i\le {\nu }_{max}\\ \frac{\left|E\left(\nu =515\mathrm{THz}\right)\right|}{{10}^{27}}& if {\nu }_i>{\nu }_{max}\end{array}\right\},\\ {\delta }_{i,2}=\left\{\begin{array}{ll}0& if {\nu }_i\ge {\nu }_{min}\\ \frac{\left|E\left(\nu =515\mathrm{THz}\right)\right|}{{10}^{27}}& if {\nu }_i<{\nu }_{min}\end{array}\right\}\end{array}$$

(38)

Optimization via BFGS algorithm: Choose the identity matrix as the initial Hessian matrix

$${\mathit{H}}_0=\mathit{I} \left(\mathit{I}:2\times 2 \mathrm{identity} \mathrm{matrix}\right)$$

$$f_{p1,0}=250 \mathrm{THz}, f_{p1,1}=245 \mathrm{THz}, f_{p2,0}=225 \mathrm{THz}, f_{p2,1}=220 \mathrm{THz}, {\alpha }_1=0.5$$

$$\nabla C_k=\left[\begin{array}{c}\frac{C\left(f_{p1,k}, f_{p2,k}\right)-C\left(f_{p1,k-1}, f_{p2,k}\right)}{f_{p1,k}-f_{p1,k-1}}\\ \frac{C\left(f_{p1,k}, f_{p2,k}\right)-C\left(f_{p1,k}, f_{p2,k-1}\right)}{f_{p2,k}-f_{p2,k-1}}\end{array}\right], \nabla C_{k+1}=\left[\begin{array}{c}\frac{C\left(f_{p1,k+1}, f_{p2,k}\right)-C\left(f_{p1,k}, f_{p2,k}\right)}{f_{p1,k+1}-f_{p1,k}}\\ \frac{C\left(f_{p1,k}, f_{p2,k+1}\right)-C\left(f_{p1,k}, f_{p2,k}\right)}{f_{p2,k+1}-f_{p2,k}}\end{array}\right]$$

(39)

• $p_k=-$${\mathit{H}}_k$$\mathit{\nabla }{\mathit{C}}_k$, ${\nu }_{ k+1}={\nu }_k+{\alpha }_k$${\mathit{p}}_k$, ${\mathit{s}}_k={\mathit{f}}_{p,k+1}-$${\mathit{f}}_{p,k},$${\mathit{\nu }}_{\mathit{k}}=\left[\begin{array}{c}{\nu }_{1,k}\\ {\nu }_{2,k}\end{array}\right]$, ${\mathit{y}}_k=$$\mathit{\nabla }{\mathit{C}}_{k+1}-$$\mathit{\nabla }{\mathit{C}}_k$, ${\rho }_k=\frac{1}{y_k{}^T{s}_k}$
• BFGS recursion:${\mathit{H}}_{k+1}=\left(\mathit{I}-{\rho }_k{\mathit{s}}_k{\mathit{y}}_k{}^T\right){\mathit{H}}_k\left(\mathit{I}-{\rho }_k{\mathit{y}}_k{\mathit{s}}_k{}^T\right)+{\rho }_k{\mathit{s}}_k{\mathit{s}}_k{}^T$$\mathit{I}:\mathrm{Identity} \mathrm{matrix}$

$${\alpha }_k=c^{( \log \left|\frac{C\left({\nu }_k\right)}{C\left({\nu }_k\right)-C\left({\nu }_{k-1}\right)}\right| ) / ( \left|\frac{C\left({\nu }_k\right)}{C\left({\nu }_k\right)-C\left({\nu }_{k-1}\right)}\right| )}$$

(40)

The base of the step size ($c)$ was a simple constant ($1<c<1.5$) and ${\alpha }_k$ was the step size at iteration k. In this example, c was chosen as c = 1.40. According to these presented formulations, the maximum amplitude of the desired harmonic that was attained in the micro-resonator (for 0 < t < 10 ps) was obtained as $\left|E\left(\nu ={\nu }_{target}=515 \mathrm{THz}\right)\right|=9.1\times {10}^7 \mathrm{V}/\mathrm{m}$, and the resulting optimal excitation frequencies were determined as ${\nu }_1=266.2 \mathrm{THz}, {\nu }_2=476.5 \mathrm{THz}$(

Table 2. BFGS updating process.

$\left|{\mathit{E}}_{\mathit{\nu }=515\mathbf{THz}}\right|$	${\mathit{\nu }}_1$	${\mathit{\nu }}_2$	${\mathit{W}}_{\mathit{e},\mathit{p}} \left(\frac{\mathit{J}}{{\mathit{m}}^3}\right)$	${\mathit{P}}_{\mathit{p}\mathit{u}\mathit{m}\mathit{p}}\left(\frac{\mathit{C}}{{\mathit{m}}^2}\right)$	k (iteration #)
1.9 $\times {10}^3 \mathrm{V}/\mathrm{m}$	250 THz	225 THz	2.4 $\times {10}^7$	0.20	1
2.6 $\times {10}^3 \mathrm{V}/\mathrm{m}$	245 THz	220 THz	$2.5\times {10}^7$	0.18	2
$3\times {10}^3 \mathrm{V}/\mathrm{m}$	248.6 THz	226.1 THz	2.2 $\times {10}^7$	0.19	4
$5\times {10}^3 \mathrm{V}/\mathrm{m}$	243.8 THz	258.5 THz	2.7 $\times {10}^7$	0.18	6
$7\times {10}^3 \mathrm{V}/\mathrm{m}$	234.3 THz	288.7 THz	2.3 $\times {10}^7$	0.18	9
1.3 $\times {10}^4 \mathrm{V}/\mathrm{m}$	249.8 THz	341.4 THz	3.5 $\times {10}^7$	0.18	12
3.4 $\times {10}^4 \mathrm{V}/\mathrm{m}$	259.9 THz	331.9 THz	5.2 $\times {10}^7$	0.21	15
2.2 $\times {10}^5 \mathrm{V}/\mathrm{m}$	272.8 THz	363.6 THz	4.6 $\times {10}^7$	0.21	18
9.4 $\times {10}^5 \mathrm{V}/\mathrm{m}$	277.7 THz	411.1 THz	3.9 $\times {10}^7$	0.18	21
4.3 $\times {10}^6 \mathrm{V}/\mathrm{m}$	270.4 THz	448.2 THz	6.5 $\times {10}^7$	0.18	24
1.1 $\times {10}^7 \mathrm{V}/\mathrm{m}$	261.5 THz	475.3 THz	8.8 $\times {10}^7$	0.20	27
2.7 $\times {10}^7 \mathrm{V}/\mathrm{m}$	265.3 THz	475.8 THz	7.1 $\times {10}^7$	0.22	30
9.1 $\times {10}^7 \mathrm{V}/\mathrm{m}$	266.2 THz	476.5 THz	$7.7\times {10}^7$	0.22	33

).

$$\begin{array}{l}{W}_{e,p}=\mathrm{Stored} \mathrm{electric} \mathrm{energy} \mathrm{density}=\frac{1}{2}{\epsilon }_{\infty }{E}^2+\frac{1}{2}EP \left(\frac{Joules}{m^3}\right), \\ E :\mathrm{Electric} \mathrm{field} \mathrm{intensity}\end{array}$$

(41)

$P: \mathrm{Polarization} \mathrm{density} \mathrm{created} \mathrm{by} \mathrm{the} \mathrm{pump} \mathrm{wave}$$\left(\frac{Coulomb}{m^2}\right), {\epsilon }_{\infty }:\mathrm{Background} \mathrm{permittivity}.$

At the end of 33th iteration, electric field at the output of the bandpass filter for 10 ps simulation time is seen in Figure 8. As seen in this figure, electric field at 515 THz reaches in a scale of 10⁸ V/m.

In Figure 9 we present spectral magnitude through the iterations. As seen in Figure 9a there is no generated target frequency (515 THz) at 6th iteration. When the iterations proceed (Figure 9b–d), wave mixing between two source waves happens, and we obtained generated target frequency at 515 THz. Compared to generation of a wave at 820 THz we obtained much more sharp spectral distribution around the target frequency of 515 THz. At the end of 33th iteration, intensity spectral density at the band-pass filter output for 515 THz is 234.1 (W/m²Hz) as seen in Figure 10.

5. Testing the Model Accuracy via Comparison with Experimental Results

In order to validate the accuracy of our model, we compared our numerical results with the well-established experimental formula of second harmonic generation efficiency for many different excitation amplitudes. The following example illustrates this comparison by computing the error percentage for different excitation amplitudes. Once the error percentage was below one percent for every excitation amplitude, we assessed the computational results to be in good agreement with the experimental results.

Example 5.1: Second Harmonic Generation by Nonlinear Wave Mixing

The goal of this example was to compute the second harmonic generation efficiency of a high intensity input wave after it had propagated through an interaction medium for a certain amount of time. The input wave was initially monochromatic with an angular frequency of ${\omega }_1$. The generated second harmonic of the input wave had an angular frequency of ${\mathsf{\omega }}_2=2{\mathsf{\omega }}_1$. Experimental and computational formulations were compared.

The 100 THz (monochromatic) high-amplitude input wave ${\mathit{E}}_1$ is excited at x = 2.4 μm (see Figure 11). The excitation amplitude of the input wave is $A_1$ (V/m).

• ${\mathit{E}}_1\left(x=2.4 \mathsf{\mu }\mathrm{m},t\right)=A_1\times \sin \left(2\pi \left(1\times {10}^{14}\right)t+{\phi }_1\right) \mathrm{V}/\mathrm{m}$$\left( {\phi }_1=0\right)$
• Spatial range and duration of the computation: $0\le x\le 10 \mathsf{\mu }\mathrm{m}, 0\le t\le 30 \mathrm{ps}$
• $\mathrm{Resonance} \mathrm{frequencies} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: {\mathit{f}}_{\mathit{r}}=\left\{7.8\times {10}^{14} \mathrm{Hz}, 9.5\times {10}^{14} \mathrm{Hz}, 1.4\times {10}^{15} \mathrm{Hz}\right\}$
• $\mathrm{Damping} \mathrm{coefficients} \mathrm{of} \mathrm{the} \mathrm{interaction} \mathrm{medium}: \mathit{\gamma }=\left\{4\times {10}^{12} \mathrm{Hz}, 3\times {10}^{12} \mathrm{Hz}, 1\times {10}^{12} \mathrm{Hz}\right\}$
• $\mathrm{Interaction} \mathrm{medium} \mathrm{background} \mathrm{permitttivity} \left({\epsilon }_{\infty }\right)=1+\chi =12 ({\mu }_r=1)$
• $\mathrm{Left} \mathrm{absorption} \mathrm{domain} \mathrm{ranges} \mathrm{from} x=0 to x=2.35 \mathsf{\mu }\mathrm{m} \left(\mathrm{absorbing} \mathrm{boundary}\right)$
• $\mathrm{Right} \mathrm{absorption} \mathrm{domain} \mathrm{ranges} \mathrm{from} x=7.65 \mathsf{\mu }\mathrm{m} to x=10 \mathsf{\mu }\mathrm{m} \left(\mathrm{absorbing} \mathrm{boundary}\right)$

The experimental formula for the efficiency of second harmonic generation is stated as [40,41,42]

$${\eta }_{\mathrm{experimental}}={(\tanh \sqrt{d^2{\eta }^3{\omega }^{2}cn{\epsilon }_0{A}_1{}^2{L}^2} )}^2$$

• d = Nonlinearity coefficient, $\eta$ = Medium impedance, n = Index of refraction
• $A_1$ = Input wave amplitude, L = Length of the medium

$${\omega }_2=\mathrm{Angular} \mathrm{frequency} \mathrm{of} \mathrm{the} \mathrm{generated} \sec \mathrm{ond} \mathrm{harmonic}$$

The computational results were obtained using Equations (12)–(16).

For a computation range of $0\le t\le t_{\max }$, the expression for the computational second harmonic generation efficiency is stated as [40,41,42]

$${\mathsf{\eta }}_{\mathrm{computational}}=\frac{\mathrm{Intensity} \mathrm{of} \mathrm{the} \sec \mathrm{ond} \mathrm{harmonic} \mathrm{of} \mathrm{the} \mathrm{source} \mathrm{wave} \mathrm{at} \mathrm{t}={\mathrm{t}}_{\max }}{ \mathrm{Intensity} \mathrm{of} \mathrm{the} \mathrm{first} \mathrm{harmonic} \mathrm{of} \mathrm{the} \mathrm{source} \mathrm{wave} \mathrm{at} \mathrm{t}=0 }$$

(42)

The parameters of the computation are as follows:

$${\omega }_1=\mathrm{First} \mathrm{harmonic} \mathrm{angular} \mathrm{frequency}=\left(2\pi \times 100\right) \mathrm{THz}$$

• L = Length of the interaction medium = 3.33 $\mathsf{\mu }\mathrm{m}$ (ranging from x = 3.33 $\mathsf{\mu }\mathrm{m}$ to 6.66 $\mathsf{\mu }\mathrm{m}$)
• ${\omega }_2=\mathrm{Second} \mathrm{harmonic} \mathrm{angular} \mathrm{frequency}=2\pi \times 200 \mathrm{THz},$n$=\sqrt{12}$
• $A_1$ = Input wave amplitude (Varied from ${10}^8\frac{\mathrm{V}}{\mathrm{m}}to 2.5\times {10}^9 \mathrm{V}/\mathrm{m}$, in increment of ${10}^8 \mathrm{V}/\mathrm{m}$)

Oscillator weights (Resonance probabilities) = ξ = {$0.3, 0.4, 0.3$}

Based on Equation (42), the computational second harmonic efficiency $\left({\eta }_{\mathrm{computational}}\right)$ was found as ${\eta }_{\mathrm{computational}}=1.74\times {10}^{-4}$ for an excitation amplitude of $A_1=3\times {10}^{8 }\mathrm{V}/\mathrm{m}$. Using this result, we tried to estimate the nonlinear coefficient of the medium by solving the following transcendental equation [40,41,42]

$${(\tanh \sqrt{d_{est}{}^2{\eta }^3{\omega }^{2}cn{\epsilon }_0{A}_1{}^2{L}^2} )}^2=1.21\times \frac{{10}^{-21}C}{V}$$

(43)

From which we solved for the estimated nonlinear coefficient as $d_{est}=1.21\times {10}^{-21}C/V$. However, in order to verify this estimated nonlinearity coefficient, we had to check whether this coefficient matched the theoretical results with the computational results for all input wave amplitudes. This verification is shown below in

Table 3. Comparison of the numerical and experimental results for different excitation amplitudes with the estimated nonlinearity coefficient of $d_{est}=1.21\times {10}^{-21}$ C/V.

Excitation Wave Amplitude (V/m)	Theoretical Efficiency	Numerical Efficiency	Error Percentage
$1\times {10}^8$	$1.93\times {10}^{-5}$	$1.92\times {10}^{-5}$	0.5
$3\times {10}^8$	$1.74\times {10}^{-4}$	$1.74\times {10}^{-4}$	0.4
$5\times {10}^8$	$4.82\times {10}^{-4}$	$4.88\times {10}^{-4}$	1.2
$1\times {10}^9$	$1.93\times {10}^{-3}$	$1.96\times {10}^{-3}$	1.55
$1.5\times {10}^9$	$4.33\times {10}^{-3}$	$4.39\times {10}^{-3}$	1.38
$2\times {10}^9$	$7.68\times {10}^{-3}$	$7.76\times {10}^{-3}$	1.04
$2.5\times {10}^9$	$1.20\times {10}^{-2}$	$1.21\times {10}^{-2}$	0.83

for some sample excitation amplitudes and illustrated in Figure 12 for a broad range of excitation amplitudes. The theoretical and computational results seem to be in good agreement, there is more than 98% accuracy between the theoretical and the numerical results.

6. Conclusions

For an arbitrary micro-resonator interaction medium of multiple resonances, the frequencies of the ultrashort excitation pulses must be tuned accordingly, based on a rapidly converging and computationally efficient Quasi-Newton algorithm (such as BFGS) in order to generate a desired harmonic at an ultra-high efficiency. The non-linear programming process is basically needed to maximize the stored spectral energy density around the desired frequency so that the conversion efficiency is maximized. For a micro-resonator medium with many resonance frequencies, the computational cost increases as the number of differential equations increases due to the occurrence of many polarization density components. However, the cost can be reduced by decreasing the number of excitation pulses and alleviating the burden on non-linear programming. Even a single parameter optimization via a single excitation pulse may allow us to achieve the target harmonic generation efficiency. However, for devices that have very limited output powers, multiple excitation pulses must be adjusted in order to achieve the optimal frequency combination that allows the highest efficiency to be attained and to compensate for the output-power limitations of the source device.