Invited Article: Acousto-optic finite-difference frequency-domain algorithm for first-principles simulations of on-chip acousto-optic devices

In recent years, there has been increasing interest in acousto-optic devices, with an emphasis on the design of on-chip structures to efficiently harness various photon and phonon interaction mechanisms such as the stimulated Brillouin scattering (SBS).^1–4 As SBS yields an extremely strong nonlinear interaction with narrow resonances, it has found applications in many important areas of optics and acoustics.^1–3 Traditionally, SBS has been studied extensively in fiber-optic devices in order to inhibit undesired nonlinear effects induced by SBS.^1–4 More recently, SBS has been tailored for micron-scale on-chip devices,^2,3,5–10 where it is considered to be an attractive candidate in the creation of lasers with ultra-narrow bandwidths,^11–14 gigahertz frequency combs,^15,16 slow light,¹⁷ and on-chip signal processing devices such as the microwave photonic filter^18,19 and optical isolators.^20–23 

Given the wide range of devices that can harness SBS, it is important to develop general numerical techniques that can facilitate the device design process. However, direct simulations of these devices face an intrinsic challenge that arises from the enormous time scale difference between optical and acoustic waves, effectively rendering tradition time-domain simulation methods intractable. For instance, a typical optical wave has a frequency of around 200 THz, whereas SBS acoustic waves usually have frequencies of around 5 to 10 GHz. Thus, even though in principle, one could simulate acousto-optic interactions with a standard first-principles time-domain simulation technique such as the finite-difference time-domain (FDTD) algorithm,^24,25 a single acoustic wave cycle corresponds to around 10⁵ optical wave cycles, and resolving an SBS resonance with a linewidth of around 1 MHz would require at least 10⁹ optical cycles. Thus, accurately treating photon-phonon interactions in time-domain simulations becomes prohibitively numerically expensive. As such, when designing fiber-optic and on-chip acousto-optic devices, researchers typically adopt a mode-expansion technique, which calculates the optical and acoustic modes independently and then treats the acousto-optic coupling perturbatively using the coupled mode theory.^{5–3,16,18–21} Unfortunately, the application of coupled mode theory is not exact and becomes difficult for complex geometries, which may support a large number of interacting modes. Therefore, in order to accurately and realistically perform first-principles simulations of a general class of acousto-optic devices, there is an urgent need to develop a computational algorithm to efficiently and exactly simulate the interactions between optical and acoustic waves.

In this paper, we introduce an acousto-optic finite-difference frequency-domain (FDFD) technique in order to perform first-principles calculations of the photon-phonon interactions in acousto-optic devices. In the frequency domain, the physics of the acousto-optic system can be rigorously formulated as a system of coupled nonlinear equations, whose solution provides the steady-state dynamics of the acousto-optic systems. With such a frequency-domain solver, we bypass the need to compute field values at every time step and can therefore directly simulate a general class of acousto-optic devices without the limitations in time-domain simulations as imposed by the vastly differing time scales between optical and acoustic waves.

The remainder of this manuscript is structured as follows. In Section II, we review the physics of optical and acoustic waves along with their interactions. In Section III, we provide a general formalism of the acousto-optic FDFD algorithm based on the wave equations in Section II. In Section IV, we demonstrate two numerical examples of the acousto-optic FDFD algorithm. The first example is a verification of this algorithm, where we observe excellent agreement between its solutions and those of the coupled mode theory. The second example is a simulation of a realistic on-chip SBS resonator, where we capture features that are prominent to the SBS process. In Section V, we provide a summary of our work as well as a general discussion regarding to the application of our algorithm.

To start, we first briefly review the physics of Maxwell’s equations^26–29 and the acoustic wave equation in the context of acousto-optic interactions.^8,25,30,31 In frequency domain, Maxwell’s equations for the electric field E$(ω)$ at frequency ω can be written as^4,32

in which the spatial dependence of the source, field, and material parameters are implicitly defined. $εr$ is the relative permittivity at frequency ω. $𝐉(ω)$ is the external current density, and $𝐏(ω)$ is the nonlinear polarization density component at frequency ω.

The frequency-domain solution of the acoustic wave equation can be formulated in a similar manner. An acoustic (mechanical) displacement field $U∼(t)$ at frequency Ω can be expressed as

where U is the complex amplitude of the displacement field. For such a wave, the fundamental equation of motion for U at frequency Ω can be expressed in the component form,^8,30

where ρ is the material density, $𝐜−$ is the stiffness tensor, $𝜼¯$ is the viscosity tensor, and F is the force acting on the acoustic wave. The frequency dependence of U_i and F_i is implicit since we are solving for the steady-state response at a single acoustic frequency $Ω$⁠. By adopting standard tensor-vector contraction notations, Eq. (3) can be written more compactly as³⁰ 

where $∇⊗𝐔$ describes the tensor derivative of U as

the operator describes a rank four tensor acting on a rank two tensor. For instance,

and $∇⋅()$ describes the divergence operator acting on a rank two tensor.

In order to treat acousto-optic phenomena, we need to incorporate both the nonlinear polarization density $𝐏(ω)$ and the force density F into Eqs. (1) and (4). In particular, one needs to include the optical and acoustic wave coupling both in the bulk of a material and at its boundaries.^7,8,33 In what follows, we show that one can explicitly treat the effects of $𝐏(ω)$ and F by developing a system of nonlinear equations that can be solved self-consistently.

When an acoustic wave $𝐔∼(t)$ with the form of Eq. (2) exists in optical media, the nonlinear polarization density $𝐏∼(t)$ has two contributions (we ignore the moving polarization effect because it is a much weaker effect as noted in Ref. 8),

where $𝐏∼(b,PE)(t)$ is the bulk polarization density from photo-elasticity (PE) and $𝐏∼(s,MB)(t)$ is the surface polarization density caused by moving boundaries (MBs). Subsequently, a superscript of (b) denotes a term acting on the bulk, whereas a superscript of (s) denotes a term acting on a surface. The photo-elastic effect can be described by an electromagnetic susceptibility $χij(PE)(t)$ caused by the acoustic wave,⁸ 

or written as a contracted tensor,

where $𝐩−$ is the rank four photo-elastic tensor. The polarization density induced by $χ(PE)(t)$ can be described as

which can be expressed in the frequency domain as

By inserting Eq. (11) into Eq. (1), we notice that the acoustic field density U causes the interaction between an optical field at frequency ω with its neighboring sideband frequency components $ω±Ω$⁠, which is a general property of acousto-optic interactions. Therefore, in the presence of an acoustic wave, we need to consider a time-domain electric field $𝐄∼(t)$ of the general form³² 

where E_m is the complex field component at frequency $ωm$⁠, and the neighboring frequencies are separated by $Ω$⁠, i.e., $ωm+1−ωm=Ω$⁠. For an optical wave equation at frequency $ωm$⁠, the exact form of the polarization density due to photo-elasticity can be expressed as

Next, we will discuss the polarization density induced by the movement of the boundary due to the acoustic field. When the acoustic field component normal to the surface causes the deformation of a structure as illustrated in Fig. 1, the electric field perturbation at the surface can be described by Ref. 8,

where the E and D fields are related by the constitutive relationship $𝐃=ε0εr𝐄$⁠. By treating $Δ𝐄$ and $Δ𝐃$ as perturbations caused by the acoustic field component normal to the material boundary, the polarization density on the surface at frequency $ωm$ as induced by moving boundaries can be calculated as

where we have defined $𝝈s$ as a one-dimensional delta function that lies on the material surface multiplied by the surface normal unit vector $n^$⁠. By substituting Eqs. (12), (13), and (15) into Eq. (1), we find the general form of an optical wave equation at frequency $ωm$ as

Having provided the general treatment for the effects of the acoustic wave on optical waves, we now describe how optical waves can produce forces that excite acoustic waves. In an acousto-optic medium, the optical waves in Eq. (12) can provide three types of mechanical forces,

where $𝐅(b,ES)$ is the bulk electrostrictive (ES) force, $𝐅(s,ES)$ is the surface electrostrictive force, and $𝐅(s,MB)$ is the surface force caused by radiation pressure (here, the superscript “MB” stands for “moving boundary”). The bulk and surface electrostrictive forces can be described in the component form, respectively, as⁸ 

With the tensor contraction notation, we can express, in the vectorial form, the sum of the forces as

For the radiation pressure force at the boundary, we use the analysis and results as derived in Ref. 8 to find

We can now obtain a general acoustic wave equation with the interaction with optical waves by substituting Eqs. (20) and (19) into Eq. (3) and get

where we included an external force, F_ext, to capture any other driving forces that are non-optical. Together, Eqs. (16) and (21) fully capture the physics behind acousto-optic interactions.

Having presented the acousto-optic equations, in this section, we introduce the finite-difference treatment of these equations in order to construct the acousto-optic FDFD algorithm and reach a self-consistent solution for the fields. The formalism in this section is completely general; for a concrete example of a two-dimensional formalism, please refer to Part 1 of the supplementary material.

The coupled nonlinear equations (16) and (21) are the basis for the acousto-optic FDFD algorithm and must be solved simultaneously for the electric fields $𝐄m≡𝐄(ωm)$ at all frequency components $ωm$ and the acoustic field U. In a simulation, we keep a total of M frequency components and ensure that the solution converges as we increase M. By doing so, we need to solve a total of M + 1 complex nonlinear system of equations. To efficiently solve such a system of nonlinear equations, we adopt the Newton-Raphson method³⁴ to iteratively compute the self-consistent solution ${𝐄m,𝐔}$⁠. The treatment below is similar in setup to the harmonic balance method for nonlinear circuit simulations^35,36 as well as other frequency domain algorithms developed to solve for the steady-state solutions of lasers while accounting for the nonlinear effects due to gain saturation.^37–40 To start, we first define a vector v that contains 2(M + 1) complex field elements

With this definition, one can rewrite Eqs. (16) and (21) along with their complex conjugate counterparts into a set of 2(M + 1) functionals g(v),

where $O^∈ℂ2(M+1)×2(M+1)$ is a block-diagonal operator that acts linearly on the fields

$𝐛∈ℂ2(M+1)×1$ is the current sources for each field

and C(v) is a set of 2(M + 1) nonlinear functionals that captures the nonlinear coupling amongst the elements of v,

In Eq. (26a), K(v) describes M set of equations that govern the optical components {E₁ E₂ E_M}, and L(v) is the set of equations that governs the acoustic wave U. From Eqs. (16) and (21), the mth set of equations in K(v) can be identified as

and L(v) can be written as

Having explicitly written out g(v) and its constituents, the 2(M + 1) × 2(M + 1) Jacobian operator D_g(v) can be computed as

where $O^$ is given in Eq. (24). To derive the second term of the Jacobian, we simply apply partial derivatives with respect to each of the constituents of v, which consists of each of the field components as shown in Eq. (22). Because of the large volume of equations involved in calculating $∂𝐂/∂𝐯$⁠, we supply the details of this computation in Part 2 of the supplementary material.

With the Jacobian D_g(v), we can apply the Newton-Raphson algorithm³⁴ to iteratively solve for the self-consistent solution for which g(v) = 0. Given the initial condition v₀ = 0, subsequent updates at the (k + 1)th step for v_k+1 can be obtained as

where s_k defines the step of the Newton-Raphson algorithm, computed by solving the following linear equation:³⁴ 

The iterative solver is terminated when convergence is reached, defined by, when the Newton step

is sufficiently small.

In the presentation above, for simplicity and clarity in the formalism, we describe the Newton-Raphson method in terms of taking derivatives with respect to both the field and its complex conjugate. In the actual numerical implementations below, we alternatively treat v as 2(M + 1) real unknowns

and formulate Eq. (23) in terms of 2(M + 1) real set of equations. Then we solve for the real and imaginary parts of the fields using the Newton-Raphson method. In the limit where the acousto-optic coupling goes to zero, the Jacobian reduces to the linear operator $D𝐠(𝐯)=O^$ independent of v, and Eq. (28) converges in one iteration to the solution that corresponds to the uncoupled linear solutions at independent acoustics and optical frequencies. Since on-chip acousto-optic coupling is relatively weak, the Newton-Raphson algorithm converges in a relatively small number of iterations.

In practice, to obtain the acousto-optic FDFD numerical solution of the system as described by Eq. (23), one can discretize the simulation domain on, for instance, the Yee lattice.⁴¹ When discretizing both the optical and acoustic parameters in the same cell, the optical fields are discretized in accordance with the standard Yee method, whereas the acoustic displacement fields are co-located with their electric field counterparts. The optical and acoustic material parameters are located at the center of each Yee cell. To construct the boundary operator $𝝈s$ as a one-dimensional delta function at the material boundary, we adopt the numerical treatment as done in Ref. 42, where we locate the boundary pixels within a material, derive its surface-normal unit vector, and assign a value of 1/$Δh$ to the $𝝈s$ term in Eqs. (16) and (21) associated with that pixel, where $Δh$ is the pixel size.

Regarding to the numerical properties of our algorithm, for concreteness, let us consider a simulation space that is discretized into K field components, and thus the discretized g(v) contains 2K(M + 1) real nonlinear equations, where $𝐯∈ℝ2K(M+1)$ contains all discretized unknowns in Eq. (31). In doing so, the Jacobian in Eq. (27) is a well-defined $ℝ2K(M+1)× 2K(M+1)$ matrix, and the Newton-Raphson update equation in Eqs. (28) and (29) can be computed at each iteration. In particular, the Jacobian D_g(v_k) is sparse, and Eq. (29) can thereby be very efficiently computed at each step with methods such as matrix factorization⁴³ for smaller problems or with various iterative methods such as the bi-conjugate gradient method⁴⁴ and the quasi-minimal residual method⁴⁵ for larger systems.

In this section, we will use the acousto-optic FDFD algorithm that we developed in Sec. III and dedicate the rest of this paper to provide validation for our algorithm as well as an example of its application to a realistic acousto-optic ring resonator. For simplicity, we restrict our analysis to two-dimensional transverse-electric (TE) optical fields (where the nonzero field components are E_z, H_x, and H_y). The details of the formalism for this two-dimensional algorithm are provided in Part 1 of the supplementary material.

In the first example, we verify the acousto-optic FDFD solution by applying it to a waveguide with SBS gain. To be concrete, we assume the propagation direction to be $x^$⁠, the transverse direction to be $y^$⁠, and the infinite out-of-the-plane direction to be $z^$⁠. In a backwards SBS gain process, one typically considers the interaction of three modes: a backward-propagating optical pump $E2$ at $ω2$ with propagation constant $β2$⁠, a forward optical Stokes wave E₁ at $ω1$ with propagation constant $β1$⁠, and an acoustic wave $𝐔=Uxx^+Uyy^$ at frequency $Ω=ω2−ω1$ and wave vector $q=β2−β1$⁠.^1–4 Upon the generation of and the mixing with the acoustic wave, the Stokes wave experiences exponential growth along its propagation direction. To maximize the gain of the Stokes field, given the generated acoustic wave vectorq, the frequency of a generated acoustic wave Ω must be approximately equal to the frequency of an acoustic guided mode of the waveguide at this wavevector q, as denoted by $ΩB$⁠. For an acoustic mode of the form $𝐔=𝐮^(y)e−iqx+iΩBt$⁠, where $𝐮^(y)$ describes the acoustic modal profile, its dispersion relation can be obtained from the solution of the following eigenvalue equations:^7,8,30

As a numerical demonstration, we consider a slab waveguide geometry that is shown in Fig. 2(a). Such a structure can be used to model the suspended waveguide geometry that has been widely used for achieving efficient on-chip SBS gain processes.^9,11 The optical and acoustic parameters of the waveguide core material are chosen to be those of chalcogenide glass, As₂S₃.⁴⁶ Optically, the core has a relative permittivity of $εWG=5.6169$⁠, and it is 20 $μm$ long and 275 nm wide. Acoustically, the core has a density of $ρ0=3200kg/m3$⁠, and in the Voigt notation, its stiffness tensor is $[c11,c12]=[18.7,6.1]GPa$⁠, and the viscosity tensor is$[η11,η12]=[1.8,1.45]mPa s$⁠.⁴⁶ Here, we assume that As₂S₃ is an isotropic material, which implies that $c66=(c11−c12)/2$ and $η66=(η11−η12)/2$⁠.³⁰ This choice is purely to make it possible to derive the analytical solution of the waveguide as to compare with the acousto-optic FDFD algorithm. However, the acousto-optic FDFD algorithm is completely general and can be applied to solids with any form of stiffness and viscosity tensors.

In this structure, the acousto-optic interaction occurs when the E_z field couples with the U_x and U_y fields through both moving boundary effects and electrostriction/photo-elasticity.^7,8 The moving boundary phenomena can be incorporated through the interplay of radiation pressure and boundary perturbation, and the electrostriction/photo-elasticity processes are effected by the photo-elastic tensor as described in Sec. II. In this 2D example, the bulk electrostrictive/photo-elastic coupling between E_z, U_x, and U_y can be described by a single photo-elasticity element p₁₂. In the As₂S₃ waveguide in Fig. 2(a), the acousto-optic interaction happens between the x = 1 $μm$ and the x = 19 $μm$ region of the waveguide core, where we assign p₁₂ = 0.24.⁴⁶ The core is surrounded by vacuum, which has a stiffness tensor of 0 and enforces the traction-free boundary conditions.³⁰ To treat the surface coupling, we use the discretized version of Eqs. (16) and (21), where we set $𝝈s$ to be $1/Δy$ on the waveguide boundary, which is normal to $y^$⁠.

To determine the parameters for the efficient excitation of the SBS process, we first calculate the waveguide’s optical and acoustic dispersion relations. This waveguide supports only one optical mode, but acoustically it is multi-moded. The optical dispersion curve is shown in Fig. 2(b), and we find that for a Stokes mode with frequency $ω1=2π×193.4$ THz (corresponding to an optical wavelength of $1.55μm$⁠) and propagation constant $β1=7.543μm−1$⁠, there exists a backward acoustic mode at $ΩB=2π×5.79$ GHz and propagation constant $q=−15.08μm−1$ that is phase- and frequency-matched with a pump mode with $β2=q+β1$ and $ω2=ω1+ΩB$⁠.

In constructing the acousto-optic FDFD simulation, we keep a total of M = 6 frequency components equally spaced at the Brillouin frequency $ΩB$⁠. For this two-dimensional TE simulation, we place E_z at the origin of each Yee cell,⁴¹ and the U_x and U_y are located halfway along the x and y edges of the cell, respectively. The optical and acoustic material parameters are placed at the center of that cell, and the boundaries of the waveguide region are defined by the pixels inside the waveguide that are immediately adjacent to vacuum (see Part 1 of the supplementary material). The spatial discretization of the simulation domain is $Δx=Δy=25nm$⁠, and the simulation domain is surrounded by 15 layers of stretched-coordinate perfectly matched layers (SC-PMLs) on all four edges to suppress spurious reflections.^25,47 A forward Stokes wave with a guided power of $1μW/μm$ is excited from the $x=0.8μm$ position. At the $x=19.2μm$ position, we inject a backward continuous pump wave at $ω2$⁠, whose normalized field profile is shown in Fig. 3(a). Under this backward SBS configuration, the generated sideband frequency components alternate between propagating forward and backward. However, our simulations do not make a priori assumptions about the propagation direction of these sidebands. Instead, the directionality is inferred from analyzing the simulation results. In Fig. 3(b), we plot the field profile of the Stokes wave when the pump power is chosen to be 100 $W/μm$⁠; such a high pump power is used in order to observe an appreciable SBS process of a relatively short waveguide. Visually, we see that the Stokes field is amplified along the propagation direction. Furthermore, the field profiles of the generated acoustic field are shown in Fig. 3(c).

We now compare the results from the acousto-optic FDFD solutions with those from the coupled mode theory (CMT) [see Part 3 of the supplementary material].⁴⁸ In Fig. 3(d), the powers of the Stokes field P₁ and acoustic field P_a from both FDFD and CMT results are plotted along the waveguide direction for various pump powers P₂. For all the pump powers analyzed, we observe remarkable agreement between the FDFD and CMT solutions, which provide a validation that our algorithm can accurately predict the physics of acousto-optic interactions. The slight disagreement for the acoustic power at a pump power of $100W/μm$ is likely due to the breakdown of the slowly-varying envelope-approximation in CMT.^8,48

Next, we analyze the convergence of the acousto-optic FDFD algorithm. In Fig. 3(e), we plot the maximum electric field amplitudes at each of the sideband frequencies, and we note that the field amplitudes decrease rapidly as the sideband frequencies deviate farther away from $ω1$ and $ω2$⁠. This justifies our choice of keeping only a relatively small number of frequency sidebands. Furthermore, in Fig. 3(f), we plot the error at each Newton step, defined as $δ(k)$ in Eq. (30), relative to the error at the first Newton step, $δrelative=δ/δ(0)$⁠. We see that the solution converged in just four iterations with the update equation in Eq. (28).

We now demonstrate the application of the acousto-optic FDFD algorithm to a realistic acousto-optic device. The structure we consider here is shown in Fig. 4(a), and it consists of an external waveguide coupled to a ring resonator. Such a resonator device has been previously demonstrated experimentally as low pump-threshold SBS lasers^11,13,14 or as a nonreciprocal light storage unit,²² and it may also be used as a microwave photonic filter.^17,18 The material of the waveguides is again chosen to be chalcogenide glass As₂S₃, whose optical and acoustic material parameters are provided before in the simulation of a straight waveguide. In the ring structure, both the external and the ring waveguides are surrounded by vacuum, and they both have widths of 275 nm—the same as the previous straight waveguide example. Because this waveguide supports numerous acoustic modes, it is difficult to have a comprehensive analytical description of the acoustic field patterns inside the ring waveguide. The ring has a center-to-center diameter of 4.506 $μm$⁠, and it is separated from the external waveguide by an edge-to-edge distance of 475 nm. In the simulation, the discretization of space is chosen as $Δx=Δy$ = 25 nm, and the simulation domain is surrounded by 15 layers of SC-PML on each boundary.^25,47

In the absence of the acousto-optic interactions, the ring resonator is critically coupled to the external waveguide at $λres=1558.29nm$ (⁠$ωres=2π×192.4THz$⁠) with a quality factor of $Q=6.25×103$ as shown in the right panel of Fig. 4(a). When operating on resonance, the optical power inside the ring waveguide is strongly enhanced over that of the external waveguide, which drastically increases the acousto-optic interaction and reduces the pump power required to observe the SBS gain.^2,3,11,13,14

To demonstrate the SBS gain from this structure, we apply a photo-elastic coefficient of $p12=0.24$ inside the ring and compute the surface term $𝝈s$ at the boundary pixels of the ring according to the surface-normal direction at each pixel. We then send in an acoustic pump wave from the right-hand side at $ω2=ωres$ with a guided power of 90 $mW/μm$ [Fig. 4(b), left]. Meanwhile, we inject a Stokes wave from the left-hand side at a wide range of the Stokes frequencies $ω1$ with a power of 1 $μW/μm$ [Fig. 4(b), right], which, together with the pump wave, generates a large number of acoustic modes inside the ring waveguide [Fig. 4(c)]. As we sweep the Stokes frequency $ω1$⁠, we find that $ΩB=ω2−ω1=2π×5.88$ GHz, the Stokes field becomes resonant with an acoustic mode, where the Stokes field is highly amplified inside the ring waveguide, resulting in an amplified transmission at the Stokes frequency $ω1$ [Fig. 4(b), right]. For the chosen pump power and ring geometry, the transmitted power at $ω1$ is 12 times stronger than the input Stokes power. In Fig. 4(d), we plot the power amplification (⁠$Pout/Pin$⁠) from this acousto-optic interaction as we vary the Stokes frequency by $ΔΩ$ around $ω2−ΩB$⁠. From this plot, we observe an SBS linewidth of approximately 13 MHz, which is congruent with the SBS linewidth of As₂S₃.⁴⁶ For the simulations above, we achieve convergence by keeping a total of $M=6$ frequency components, and the acousto-optic algorithm converges in four steps of the Newton-Raphson algorithm.³⁴ 

There are several interesting observations that we can make from the simulations above. First, we note that although we are using the same waveguide geometry, the Brillouin frequency for the ring waveguide differs from that of the straight waveguide by 90 MHz, which captures the change in the acoustic wave vector due to the bending of the ring structure. In addition, another interesting feature is that the gain spectrum exhibits an asymmetric Fano-resonance line shape as seen in Fig. 4(d). This arises from the interference between the sharp acousto-optic resonance mode and the much broader mode of the ring resonator.⁴⁹ This intriguing detail is typically neglected in the descriptions of experimental observations^11,13,22,23 but can be captured through first-principles calculations via the acousto-optic FDFD method.

In summary, we have presented a numerically efficient, first-principles method for simulating SBS in optical devices. Although both devices simulated here consisted of two-dimensional structures with a transverse electric polarization operating using the backwards SBS configuration, the theory underlying the acousto-optic FDFD algorithm is completely general to acousto-optic and optomechanic wave phenomena, and thus this algorithm can be extended to three dimensions, as well as to other forms of acousto-optic interactions, such as the forward SBS process or an on-chip acousto-optic modulator. Furthermore, the concept behind this algorithm is not restricted by the method with which we discretize the simulation domain, so it can also be formulated for other first-principles frequency-domain techniques such as the finite element method (FEM), where a different discretization scheme is used.⁵⁰ In fact, since FEM is formulated as the solution to boundary value problems and is superior to FDFD in modeling curved surfaces, we should expect the equivalent FEM formalism of the acousto-optic interaction to be more accurate in handling moving boundary effects.

When applying the acousto-optic FDFD algorithm, the most computationally expensive step lies in solving the linear equation in Eq. (29). In the two-dimensional acousto-optic simulations above, since the resulting Jacobian in Eq. (27) is sparse and not symmetric, we use the UMFPACK package built within MATLAB to efficiently factorize the Jacobian matrix and solve Eq. (29).^43,51 On a computer cluster, using 12 cores of CPU, where each core is an AMD Opteron 6386SE (2.8 GHz/16 MB/140 W) processor, each of the first simulation examples converged in 40 min, whereas each of the second examples converged in 13 min. For larger acousto-optic FDFD simulations such as those performed in three dimensions, one needs to resort to iterative techniques, such as bi-conjugate gradient⁴⁴ or quasi-minimal residual methods⁴⁵ for solving a larger system of equations. In using an iterative solver, we can expect a similar convergence property as described by Ref. 27. Depending on the choice of PML and conditioning of the Jacobian matrix in Eq. (19), one may need to precondition the Jacobian matrix as detailed in Refs. 27 and 28 to obtain solutions with accelerated convergence.

We should note that despite the versatility of the acousto-optic FDFD algorithm, there are some limitations. First, from a computation point of view, the algorithm may be incapable of simulating three-dimensional devices that are larger than several hundred microns in length. This algorithm is much better suited for compact micron-scale acousto-optic devices with a complex geometry, where many optical and acoustic modes would interact in nontrivial ways such that a modal description is difficult or intractable. Second, we made the underlying assumption that there exists only one acoustic frequency. While this is true for the vast majority of SBS devices, there are other structures that harness a cascaded SBS process that produce acoustic waves at various frequencies.⁵² Furthermore, this algorithm is not designed to handle thermally generated phonons with a broadband of frequency components.⁴ To include the generation of other acoustic frequencies, we may use the same concept developed in Sec. II to construct a larger system of nonlinear equations and capture the interactions amongst all of the frequency components. However, if the number of optical and acoustic frequency components becomes too large, the solution to the discretized system could become computationally infeasible. Lastly, in understanding practical SBS devices, it is important to treat the effect of pump fluctuation on the linewidth of the device. The formalism presented in the paper does not directly treat such an effect of pump fluctuation. Nevertheless, one can imagine a treatment where one calculates the response of a structure to a pump at a given frequency, and then determine the effect of pump frequency fluctuation by a perturbation approach.^14,40,42,53

At the final stage of the revision, it was brought to our attention that concurrent to our work, there is another proposal for performing first-principles simulations of the SBS process using a transformation optics approach.⁵⁴ Both our work and the work in Ref. 54 point to the emerging importance of performing first-principles simulations of photon-phonon interactions for the design and characterization of acousto-optic devices.

See supplementary material for an implementation of the acousto-optic FDFD algorithm in two dimensions, the derivation of the Jacobian operator, and the coupled mode theory calculation of SBS in a waveguide.

We gratefully acknowledge fruitful discussions with Martin M. Fejer, Michael J. Steel, and Christopher Sarabalis. This work is supported by United States Air Force Office of Scientific Research (USAFOSR) grants (Grant Nos. FA9550-12-1-0024, FA9550-15-1-0335, and FA9550-17-1-0002). Y. Shi in addition acknowledges the support of a Stanford Graduate Fellowship.