In recent years, there has been very significant interest in designing active nanophotonic devices. For instance, there are extensive works with electro-optical modulators for on-chip optical interconnects [1–15]. Such modulators have also been used in the construction of optical isolators [16–19] and nonreciprocal metasurfaces [20–22], as well as in the realization of photonic topological effects [17,23–25]. In designing such active nanophotonic devices, it is crucial to be able to simulate their performance from the first principle of Maxwell’s equations. Nevertheless, simulations of these devices present an intrinsic challenge that arises from the inherent time scale differences in the different physical processes in such systems. For example, in electro-optical modulation devices, the modulation frequencies are typically on the order of 10 GHz [1–15], while the optical carrier waves have frequencies around 200 THz. Although, in principle, one could simulate such a modulator with a standard first-principles time-domain simulation technique such as the finite-difference time-domain algorithm [6,16,19,20,26], a single modulation period corresponds to around ${10}^5$ optical cycles; obtaining device performance directly from solving time-domain Maxwell’s equations is thus prohibitively expensive. Due to this constraint, many previous active device performance instances are were only simulated approximately, such as by using a combination of coupled mode theory with passive device simulations or semi-analytical equations [3,7,8,11–13], deriving approximate circuit models for waveguide elements [14,15,21], or adopting the slowly varying envelope approximation [24,25]. None of the previous methods is exact, and they make approximations that inevitably neglect portions of the underlying device physics. In order to accurately and realistically simulate these active devices from first principles, there is an urgent need to develop a computational algorithm to efficiently perform multiple time scale simulations.

In this Letter, we introduce a multi-frequency finite-difference frequency-domain (MF-FDFD) technique in order to perform first-principles calculations of active devices that possess multiple time scales. In the frequency domain, the physics of active devices can be rigorously formulated in terms of interactions between waves at different frequencies. Therefore, by setting up frequency-domain simulations simultaneously at different frequencies and allowing interactions between wave components at these frequencies, one can directly simulate a large class of active devices without the limitation in time-domain simulations, as imposed by the existence of vastly differing time scales inherent in these problems.

To start, we first briefly review the formalism of frequency-domain simulation techniques for Maxwell’s equations [27–30]. In frequency domain, Maxwell’s equations for the electric field $\mathbf{E}(\omega )$ at frequency $\omega$ can be written as

In order to treat active optical devices, we incorporate the polarization density $\mathbf{P}(\omega )$ that arises from permittivity modulation into Eq. (1). In what follows, we show that one can explicitly treat the effect of $\mathbf{P}(\omega )$ by developing a matrix equation similar to Eq. (2). This treatment is similar in setup to those described in the harmonic balance method used for nonlinear circuit simulations [32,33], as well as frequency-domain algorithms developed to solve for the steady-state solutions of lasers under nonlinear effects due to saturation gain [34]. Consider a device whose relative permittivity function is modulated in time with the form

By Fourier transforming Eq. (4) in time, we obtain

Consider the case when an active optical device is excited at a frequency ${\omega }_0$ with current density $\mathbf{J}({\omega }_0)$. Modulation at a frequency $\mathrm{\Omega }$ induces a number of sideband components at frequencies ${\omega }_n={\omega }_0+n\mathrm{\Omega }$, where $n$ is an integer. Therefore, under modulation, the time-domain electric field $\tilde{\mathbf{E}}(t)$ in the system in general is of the form

By discretizing Eq. (7) and keeping a total of 2N + 1 frequency components, we obtain a linear matrix equation of the form of Eq. (2), which is the equation of the MF-FDFD algorithm for active devices. The solution $\mathbf{e}$ can be obtained through conventional linear algebra techniques such as matrix factorization [35] or various iterative methods such as the biconjugate gradient method [36] and quasi-minimal residual method [37]. Once $\mathbf{E}({\omega }_n)$ for $-N\le n\le N$ are calculated from Eq. (7), the time-domain electric field $\tilde{\mathbf{E}}(t)$ can be recovered from Eq. (6). Compared to other algorithms used in nonlinear circuit and laser simulations [32–34], which all require the solution of a system of nonlinear equations, the MF-FDFD technique requires the solution of a system of linear equations and is thus less computationally complex.

Having presented the theoretical formalism for the MF-FDFD method, we first verify the MF-FDFD algorithm by demonstrating modal conversion in a waveguide structure with harmonically modulated permittivity. We consider a slab waveguide geometry, as shown in Fig. 1(a). This geometry has been previously considered as the basis for a dynamic isolator [18]. The waveguide core has a relative permittivity of ${ϵ}_{\mathrm{WG}}=4$ that is surrounded by vacuum, and it is 10 μm long and 750 nm wide. In the vicinity of $\lambda =1.5\mathrm{\mu m}$, the waveguide supports two transverse electric (TE) modes (with nonzero $E_z$, $H_x$ and $H_y$ field components): the even ${\mathrm{TE}}_0$ mode and the odd ${\mathrm{TE}}_1$ mode. In Figs. 1(b) and 1(c), we show the profiles of a ${\mathrm{TE}}_0$ mode at ${\lambda }_0=1.55\mathrm{\mu m}$ (${\omega }_0=1.215\times {10}^{15}\mathrm{rad}/\mathrm{s}$) and ${\mathrm{TE}}_1$ mode at ${\lambda }_1=1.50\mathrm{\mu m}$ (${\omega }_1=1.256\times {10}^{15}\mathrm{rad}/\mathrm{s}$) with propagation constants ${\beta }_0=7.47{\mathrm{\mu m}}^{-1}$ and ${\beta }_1=5.73{\mathrm{\mu m}}^{-1}$ respectively. To achieve modal conversion between the modes above, we apply a modulation frequency of $\mathrm{\Omega }={\omega }_1-{\omega }_0=40.5\times {10}^{12}\mathrm{rad}/\mathrm{s}$ with a strength of $|\delta |=0.1$ in the 5 μm long modulation region shown in Fig. 1(a), where the top half and bottom half of the waveguide are modulated with a relative phase of $\pi$ to maximize the coupling between the ${\mathrm{TE}}_0$ and ${\mathrm{TE}}_1$ modes.

 

Fig. 1. (a) Schematic of the waveguide structure with the source location and modulation region indicated. In the modulation region, the top and bottom halves of the waveguide are modulated with a relative phase difference of $\pi$. (b), (c) Modal profiles of the ${\mathrm{TE}}_0$ and ${\mathrm{TE}}_1$ modes of the waveguide at ${\lambda }_0=1.55\mathrm{\mu m}$ and ${\lambda }_1=1.50\mathrm{\mu m}$, respectively. Both modes carry a guided power normalized to $1\mathrm{W}/\mathrm{\mu m}$.

Download Full Size | PDF

In the simulation, we set up the MF-FDFD algorithm to analyze a total of $2N+1=9$ frequency components. The spatial discretization of the simulation domain is $\mathrm{\Delta }x=\mathrm{\Delta }y=25\mathrm{nm}$, and the simulation space is surrounded by 20 layers of stretched-coordinate perfectly matched layers (SC-PMLs) to suppress spurious reflection at the boundaries [38]. At the $x=1\mathrm{\mu m}$ position, we excite a continuous wave at ${\omega }_0$ (${\lambda }_0=1.55\mathrm{\mu m}$), whose guided power is normalized to $1\mathrm{W}/\mathrm{\mu m}$. As the wave propagates inside the modulation region, part of its amplitude is converted to the ${\omega }_n={\omega }_0+n\mathrm{\Omega }$ components. First, in Fig. 2(a) we plot the maximum amplitude of the electric field pattern at each ${\omega }_n$ (with $n\in [-4,4]$), and note that the maximum electric field amplitude decreases exponentially as we increase $|n|$. This indicates that convergence of the solution can be reached by considering a relatively small number of frequency components.

 

Fig. 2. (a) Plot of maximum field amplitude that exists in each frequency sideband as indexed with $n$. The field amplitude decreases exponentially with $n$. (b) (left) Field profiles of frequency components with $n\in [-2,2]$ and (right) comparison of the MF-FDFD results with coupled mode theory (CMT).

Download Full Size | PDF

Next, we examine the details of the MF-FDFD solution. On the left side of Fig. 2(b), we provide the field profiles for the $n\in [-2,2]$ sidebands, which have frequencies from ${\omega }_0-2\mathrm{\Omega }$ to ${\omega }_0+2\mathrm{\Omega }$. On the right side of Fig. 2(b), we plot the power amplitude along the propagation direction for each of the sidebands and compare them with coupled mode theory (CMT) [39] [for the CMT treatment, see Supplement 1]. We observe excellent agreement between the solution from the MF-FDFD method and CMT, which provides validation for our algorithm.

It is interesting to note that the MF-FDFD solution reveals additional features that are neglected in standard CMT calculations. On the right side of Fig. 2(b), for the MF-FDFD solution, there exist small oscillations in the power amplitudes within the modulation region that are absent in CMT, and we attribute this effect to the generation of backward propagating modes. For instance, for the field component at the frequency ${\omega }_0+\mathrm{\Omega }$, the periodicity of such small oscillations is ${\mathrm{\Lambda }}_1=2\pi /({\beta }_0+{\beta }_1)=0.476\mathrm{\mu m}$, which corresponds to the highly phase-mismatched coupling to a backward mode from modulation. The generation of such a backward wave can be further observed visually from the field profiles in Fig. 2(b) to the left of the modulation region ($x\le 3\mathrm{\mu m}$). Due to the large phase mismatch, backward coupling is typically ignored in standard CMT treatments. On the other hand, the MF-FDFD solution, being a first-principles method, can capture all possible wave interactions in a system.

Having observed excellent agreement between the MF-FDFD algorithm and CMT, we now demonstrate the application of the MF-FDFD algorithm to a realistic system with small modulation strength and slow modulation frequency compared with the optical frequency. The structure we consider consists of an external waveguide coupled to a modulated ring resonator, as shown in Fig. 3(a). Such a device has previously been demonstrated as a miniaturized on-chip electro-optical modulator [1,2]. In our structure, both the external and ring waveguides have a width of 280 nm. Both waveguides are single mode. The external waveguide has a relative permittivity of ${ϵ}_{\mathrm{WG}}=4$. The ring has a diameter of 5.80 μm as measured at the middle of the ring, and it is separated from the external waveguide by an edge-to-edge distance of 480 nm. In the simulation, the discretization of space is chosen as $\mathrm{\Delta }x=\mathrm{\Delta }y=40\mathrm{nm}$, and the simulation domain is surrounded by 20 SC-PMLs on each boundary [38].

 

Fig. 3. (a) Schematic of the ring resonator geometry. (b) At $\lambda =1564.55\mathrm{nm}$, the ${\mathrm{TE}}_0$ mode from the waveguide is critically coupled to the ring resonator. (c), (d) Field profiles for waves at the $n=1,2$ sidebands, respectively. The power inside the ring at these sidebands is coupled to the waveguide and can be detected at the probe. (e) Power amplitude at the probe at each frequency. (f) Plot of power at the probe as a function of time.

Download Full Size | PDF

In the absence of modulation, the relative permittivity of the ring is set at ${ϵ}_{\text{ring}}=4-5\times {10}^{-4}i$ so that at ${\lambda }_0=1564.55\mathrm{nm}$, the ring resonator is critically coupled to the external waveguide—the transmission reaches zero due to a combination of material and radiation losses of the ring, as shown in Fig. 3(b).

To include modulation, we assume that the permittivity of the ring has the form described in Eq. (3), with $\mathrm{\Omega }=10\times {10}^9\mathrm{rad}/\mathrm{s}$ and $\delta /{ϵ}_{\mathrm{WG}}=5\times {10}^{-4}$, both of which are achievable in realistic electro-optical modulators [2,3,40]. For the input wave with the same wavelength of ${\lambda }_0=1564.55\mathrm{nm}$ as above, the generated sideband components can now couple from the ring to the waveguide, resulting in nonzero transmission. The field profiles of the first two sidebands are shown in Figs. 3(c) and 3(d), and the transmitted power amplitudes at all nine frequency components are shown in Fig. 3(e). The transmitted power decreases exponentially with $n$, and the dominant transmitted power occurs at the first sidebands where $\omega ={\omega }_0\pm \mathrm{\Omega }$. In Fig. 3(f), the transmitted power upon modulation is plotted in time, and the power follows the waveform of the sinusoidal modulation. The calculation here provides a first-principles analysis of a ring resonator electro-optic modulator with no further approximations.

In the examples presented above, as an illustration we performed two-dimensional MF-FDFD simulations under transverse electric polarization. Nevertheless, because the theory behind our algorithm is completely general, this algorithm can be applied to more complicated simulation setups, such as three-dimensional, nonlinear, or anisotropic simulations [41]. Furthermore, we are free to prescribe a spatially varying modulation phase profile in the simulation domain, which makes it possible to simulate a traveling wave index modulation system such as those induced by acousto-optics. Finally, our algorithm is not restricted by the choice of spatial discretization methods, and therefore it can be formulated for other first-principles frequency-domain simulation techniques such as finite element analysis [42].

In obtaining the two-dimensional MF-FDFD solutions above, since the discretized system equation in Eq. (7) is sparse and not symmetric, we use the UMFPACK package built into MATLAB to efficiently solve Eq. (7) [35,43]. For three-dimensional MF-FDFD simulations, one can use iterative techniques for solving large linear systems, such as the biconjugate gradient [36] or quasi-minimal residual [37] method, in which case one may need to precondition the system matrix as detailed in Ref. [28] to obtain solutions with accelerated convergence.

In summary, we have presented the formulation and applications of the MF-FDFD algorithm that can efficiently perform simulations where a large time scale difference exists. When incorporating the polarization density into Maxwell’s equations, we can construct system equations that capture the interactions between different frequency components in modulated optical systems. Using this general technique, we demonstrate that one can effectively simulate, from first principles, devices that have frequency components with arbitrarily large discrepancies in their time scales.