Nonlinear optical frequency conversion was observed at the beginning of the laser era [1], and since then, it has been one of the most active areas of applied and fundamental research [2]. The development of high-power femtosecond and attosecond laser systems has opened numerous opportunities to use nonlinear processes for the generation of electromagnetic spectra between extreme UV (XUV) and THz frequency bands [3–5]. This has created a demand for an efficient theoretical and computational framework to model ultra-broad spectra, which could not be described by the mid-20th century theoretical tools of nonlinear optics. A variety of generalized envelope equations [6–10], and a class of carrier-resolved models, also known as unidirectional pulse-propagation equations (UPPE) [11–15], have met this demand.

Optical frequency combs and associated soliton pulses in microresonators have recently become a fast-growing area and a contender in the generation of temporally ultra-short and spectrally ultra-broadband waveforms [16,17]. The advantages of microresonators include a small footprint, small operational powers, and dispersion engineering used for spectral tuning of the generated combs. Microresonator comb spectra spanning across multiple octaves from mid-infrared to visible have been observed in [18–22]. Contrary to filamentation and supercontinuum experiments using high-energy ultra-short laser pulses, a typical microresonator experiment requires milliwatts of CW power, while short pulses are generated inside the resonator.

A numerical model well suited to describe comb generation in microresonators can be formulated using a coupled-mode theory [23–25]. This approach leads to a set of ordinary differential equations for amplitudes of the resonator modes evolving in time. Fourier transforming known coupled-mode equations to physical space, one finds the driven-damped nonlinear Schrödinger equation, which is known in the resonator context as the Lugiato–Lefever equation, or other envelope models [25–29]. Approximations underpinning the envelope models rely on the selection of a reference carrier frequency and wavenumber and their subtraction from the rest of the spectrum.

Our aim here is to bypass these approximations and present a carrier-resolved real-field theory for multi-octave frequency combs. We want this theory to capture desirable features of the carrier-resolved models known for propagation in bulk crystals, gases, and waveguides [12–15]. One such property is using the real-valued electric field, ${\cal E}$, oscillating with the required spectrum of electromagnetic frequencies. The real field can be conveniently used to express the second-, ${{\cal E}^2}$, and third-order, ${{\cal E}^3}$, effects without omitting any nonlinear terms, which is an unavoidable shortcoming of envelope based methods [28,29]. More generally, the carrier-resolved approach has not yet been formulated for a broad class of resonator systems where pump, gain, and dissipation are intrinsic for device operation [30–33], and where current trend towards the generation of increasingly extreme waveforms [3,4] requires development of novel theoretical methods beyond the envelope Haus-master and Ginzburg–Landau equations [30].

To develop the carrier-resolved theory of ring microresonators, we proceed with the initial steps of the formalism described in [29]. This formalism seeks the electric field components, ${{\cal E}_\alpha}$, obeying Maxwell equations as a superposition of the resonator modes ${F_{\alpha j}}(r,z){e^{ij\vartheta - i{\omega _j}t}}$, where $j$ is an azimuthal mode number, $\alpha$ can be $x$, $y$, or $z$, $\vartheta = [0,2\pi)$ is the polar angle, $r$ is the distance in the $(x,y)$ plane, $z$ axis is perpendicular to the resonator plane, and ${\omega _j}$ are the resonance frequencies. We account only per one mode for every $j$, and assume that this mode is dominated by one polarization direction, which allows us to omit the subscript $\alpha$,

In a resonator, every mode experiences damping with the rate ${\kappa _j}$, while the comb generation is achieved by pumping one of the modes, e.g., $j = M$, with laser light having frequency ${\omega _p}$ close to ${\omega _M}$. Substituting Eq. (1) into Maxwell equations, we find a set of equations for modal amplitudes:

${\cal N}$ is the nonlinear part of material response, which in our case includes second-, ${\tilde \chi ^{(2)}}$, and third-order, ${\chi ^{(3)}}$, nonlinear susceptibilities:

Second-order susceptibility can vary in $\vartheta$ if periodic poling is applied for quasi-phase-matching, i.e., ${\tilde \chi ^{(2)}}(\vartheta) = {\tilde \chi ^{(2)}}(\vartheta + 2\pi /Q)$, where $Q$ is the number of periods along the ring [34],

Detailed derivation of Eq. (2) can be found in [29], Eq. (3.2) and Eq. (6.1). Equation (2) is also structurally the same as Eq. (24) in [24], where ${\cal N}$ was the Kerr part of the refractive index. A standard way to simplify and then solve Eq. (2) is to introduce several groups of modes centered around the pump and its harmonics and simplify nonlinear terms by limiting nonlinear wave-mixing processes to the ones that resonantly excite selected modes, which was also the approach used in [29]. However, this method is compromised when the individual modal groups approach the octave width and omitted non-resonant nonlinear terms start playing a role.

Equation (2) is free from the above limitation. It treats all the modes on equal footing and accounts for all nonlinear wave-mixing processes by dealing with ${{\cal E}^2}$ and ${{\cal E}^3}$ without simplifications. However, substituting Eq. (1) into the right-hand side of Eq. (2) yields tens of thousands of nonlinear coefficients varying with $j$ within around one order of magnitude, and the equations appear to be entirely not practical even before numerical integration in time is approached. To derive a model with the practical number of nonlinear coefficients and, at the same time, to retain all terms inside ${{\cal E}^2}$ and ${{\cal E}^3}$, we make a simplification that for every $j$ there is the same coefficient for all nonlinear mixing terms entering the respective equation, but the coefficients vary with $j$. This is achieved by assuming $\iiint {{\cal E}^3}{F_j}{e^{- ij\vartheta}}r{\rm{d}}r{\rm{d}}z{\rm{d}}\vartheta \approx\def\LDeqbreak{}b_j^3\iint\! F_j^4r{\rm{d}}r{\rm{d}}z\int_0^{2\pi} {\cal E}_0^3{e^{- ij\vartheta}}{\rm{d}}\vartheta$ and $\iiint\! g(\vartheta){{\cal E}^2}\!{F_{\!j}}{e^{- ij\vartheta}}r{\rm{d}}r{\rm{d}}z{\rm{d}}\vartheta \def\LDeqbreak{}\approx b_j^2\iint F_j^3r{\rm{d}}r{\rm{d}}z\int_0^{2\pi} g(\vartheta){\cal E}_0^2{e^{- ij\vartheta}}{\rm{d}}\vartheta$, where

Equation (7) is our main result. Structurally, it can be matched with the carrier-resolved pulse propagation models; see, e.g., Eq. (85) in [12] and references therein. Apart from the pump and loss terms, a difference to be noted is that Eq. (7) describes the evolution of spatial harmonics (resonator modes) in time. At the same time, the aforementioned carrier-resolved [12] and multi-octave envelope models [9,10,28] deal with the evolution of discretized frequency spectrum in the propagation coordinate. Hence, unlike the above references, to reconstruct the frequency spectrum, we need to store the data in time and then Fourier transform them to frequency, $t \to \omega$, as post-processing. An advantage of our method for describing resonator systems generating narrow spectral features is that the frequency grid is not preset, and frequencies are captured as generated.

As a practical example, we consider a thin-film lithium niobate (LN) microresonator with a geometry close to that implemented in [36]. The ridge width at the top surface is 1.45 µm, ridge height is 360 nm, wall angle is 60°, total LN thickness is 560 nm, and the internal radius is 50 µm, so the repetition rate at ${\omega _M}/2\pi \approx 193\;{\rm{THz}}$ is $394\;{\rm{GHz}}$. Experimental measurements demonstrated that such geometry without periodic poling and pumped with ${\sim}100\;{\rm{mW}}$ of power produces 100 THz wide Kerr frequency combs with two strong dispersive waves emitted close to 150 and 250 THz [36]. For the data sets with ${\omega _j}$, ${\gamma _{2,j}}$, ${\gamma _{3,j}}$, and $g(\vartheta)$, see [35].

The nonlinear part of Eq. (7) is formulated in a pseudo-spectral form [12], i.e., it requires going back to physical space to compute ${{\cal E}_0}$. Fully spectral representation is also possible but computationally less efficient. In solving Eq. (7) numerically, the integrals and transformation from the real field, ${{\cal E}_0}$, to the mode amplitudes, ${B_j}$, and back, $\vartheta \leftrightarrow j$, were implemented using fast Fourier transform (FFT). The evolution in time was computed using the third-order Runge–Kutta method. The number of modes was set to embrace the third harmonic. In that case, the maximal frequency of oscillations produced by the ${\gamma _{3,j}}$ terms is well estimated by ${\varepsilon _3} = 3{\omega _M} - {\omega _{3M}}$. ${\varepsilon _3}$ is zero in the absence of dispersion (equidistant resonator modes), revealing that the maximal suitable time step is determined not by the optical frequency but by the resonator dispersion, which is enhanced in integrated photonics due to tight structural confinement. In our case, $|{\varepsilon _3}|/2\pi \approx 75\;{\rm{THz}}$. In practice, stability and convergence of numerical code were achieved for steps $\le \pi /4|{\varepsilon _3}| \approx 1.7\;{\rm{fs}}$. If the third-harmonic modes are set to zero, then the near-phase-matched resonators can be integrated using time steps ${\sim}100\;{\rm{fs}}$, determined by the minimal frequencies produced by the ${\gamma _{2,j}}$ terms, $|{\varepsilon _2}|/2\pi$; see Eq. (5). Setting second and third harmonics to zero increases permitted steps to above 1 ps.

In our resonator, modes $M = 404,405$ (${\omega _M}/2\pi \approx 193\;{\rm{THz}}$) are quasi-phase-matched with $2M + Q$ for $Q = 87,88,89$, so that $|{\varepsilon _2}|$ reaches its minimum, where it becomes an order of $100\;{\rm{GHz}}$. We start presenting our simulation results from the strongly mismatched case of $Q = 10$. Bistable response of the pumped resonance for power ${\cal W} = 200\;{\rm{mW}}$, which includes (black line) and excludes (red line) ${\chi ^{(2)}}$ nonlinearity, is shown in Fig. 1(a). One can see that this arrangement corresponds to the quasi-Kerr regime since second-order nonlinearity makes little impact. Therefore, we expect to find Kerr-soliton combs with two dispersive waves as were observed in [36]. The shaded area in Fig. 1(a) shows an interval of the pump detuning where we have found solitons and breathers by solving Eq. (7). A typical spectrum and carrier-resolved soliton pulse are shown in Figs. 2(a) and 2(b), respectively.

 

Fig. 1. Black (red) lines show power of the pump mode, $|{B_M}{|^2}$, versus detuning, $\delta$, when ${\tilde \chi ^{(2)}} + {\chi ^{(3)}}$ (only ${\chi ^{(3)}}$) nonlinearities are accounted for. $Q$ is the number of poling periods, and ${\cal W}$ is laser power. Shading shows intervals investigated for frequency comb generation using the carrier-resolved model, Eq. (7). Nonlinear parameters are ${\gamma _{2,2M + Q}}/2\pi = 8.39\;{\rm{GHz}}/\sqrt {\rm{W}}$ and ${\gamma _{3,M}}/2\pi = 3.79\;{\rm{MHz/W}}$, linewidth is $\kappa /2\pi = 100\;{\rm{MHz}}$, and coupling parameter and finesse are $\eta = 0.5$ and ${\cal F} = 3940$, respectively.

Download Full Size | PDF

 

Fig. 2. (a), (b) Spectrum [Fourier transform, $t \to \omega$, of ${{\cal E}_0}(t,\vartheta = 0)$] and real field profile, ${{\cal E}_0}(t,\vartheta = 0)$, of the Kerr soliton generated in strongly mismatched resonator: $Q = 10$, $M = 404$, ${\varepsilon _2}/2\pi = 30\;{\rm{THz}}$. (c), (d) Three-octave spectrum corresponding to the randomly fluctuating electric field generated in quasi-phase-matched resonator: $Q = 89$, $M = 404$, ${\varepsilon _2}/2\pi = - 445\;{\rm{GHz}}$, ${j_{{\min}}} = 147$, ${j_{{\max}}} = 1478$. Pump is ${\cal W} = 200\;{\rm{mW}}$ in both cases. Total integration time is 191 ns in (a) and 143 ns in (b). Spectrum was computed over last 40 ns in (a) and 95 ns in (b). Round trip time is 2.538 ps.

Download Full Size | PDF

Selecting $Q = 89$, we bring the resonator to quasi-phase-matching. Figures 1(b)–1(d) show the respective bistable curves for ${\cal W} = 200\;{\rm{mW}}$, $20\;{\rm{mW}}$, and $7\;{\rm{mW}}$ of pump power. For all three, tilts of the resonance happen dominantly through the ${\chi ^{(2)}}$ nonlinearity, and the soliton and comb generations in the shaded intervals of detunings were found outside the Kerr-only resonance (red line). Spectra computed in the $200\;{\rm{mW}}$ case continuously extend across three octaves; see Fig. 2(c). They correspond to the incoherent randomly varying fields, whose carrier-resolved form is shown in Fig. 2(d).

Reducing power to $20\;{\rm{mW}}$ and then to $7\;{\rm{mW}}$, we have observed formation of ${\chi ^{(2)}}$ solitons and breathers; see Fig. 3. For higher powers, these solitons are surrounded by strong dispersive waves emitted around the pump and second harmonic. Dispersive wave emission is greatly reduced for powers less than $10\;{\rm{mW}}$, where the soliton combs become spectrally narrow. Spectral resolution in our approach is set by inverse integration time, which is 88 ns in Fig. 3, providing 11 MHz resolution. A section of the recorded pulse train, zoomed-in spectra, and real-field profiles of individual solitons are also shown.

 

Fig. 3. ${\chi ^{(2)}}$ solitons generated in quasi-phase-matched resonator: $Q = 89$, $M = 405$, ${\varepsilon _2}/2\pi = - 445\;{\rm{GHz}}$, for relatively small powers. (a)–(c) Spectrum, section of the soliton train, and real-field profile of a single soliton for ${\cal W} = 20\;{\rm{mW}}$, respectively. (d)–(f) Spectrum, its zoomed interval showing the resolved comb lines, and real-field profile of a single soliton for ${\cal W} = 7\;{\rm{mW}}$, respectively. Total integration time is 190 ns in (a) and 143 ns in (b). Spectrum was computed over last 80 ns in (a) and 88 ns in (d). Round trip time is 2.538 ps.

Download Full Size | PDF

In summary, we developed the real-field carrier-resolved model to describe optical frequency comb generation in ring microresonators. Our model is capable of reproducing frequency combs with multi-octave spectral width, while simultaneously accounting for every mode physically present in the resonator and all nonlinear wave-mixing terms entering the nonlinear polarization specified via the real-valued electric field; see Eq. (3). For example, if we consider ${\chi ^{(2)}}$ nonlinearity and two envelope functions ${A_1}$ and ${A_2}$ centered around the pump and its second harmonic, the opening of the second-order nonlinear polarization ${\chi ^{(2)}}{({A_1}{e^{i{\phi _1}}} + {A_2}{e^{i{\phi _2}}} + \text{c.c.})^2}$ yields a combination of 10 terms. Here, ${\phi _{1,2}}$ are the phases oscillating with optical frequencies. Omission of any term and other common approximations fail to be justifiable if the spectra of ${A _{(1,2)}}$ approach an octave width, while our method is free from these constraints and solves the problem for nonlinearities of arbitrary orders. Overall, our results open an avenue for theoretical research into a growing family of resonator systems requiring carrier-resolved multi-octave modeling of ultra-short pulses [32,33,37].