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

We introduce a ﬁnite-difference frequency-domain algorithm for coupled acousto-optic simulations. First-principles acousto-optic simulation in time domain has been challenging due to the fact that the acoustic and optical frequencies differ by many orders of magnitude. We bypass this difﬁculty by formulating the interactions between the optical and acoustic waves rigorously as a system of coupled nonlinear equations in frequency domain. This approach is particularly suited for on-chip devices that are based on a variety of acousto-optic interactions such as the stimulated Brillouin scattering. We validate our algorithm by simulating a stimulated Brillouin scattering process in a suspended waveguide structure and ﬁnd excellent agreement with coupled-mode theory. We further provide an example of a simulation for a compact on-chip resonator device that greatly enhances the effect of stimulated Brillouin scattering. Our algo-rithm should facilitate the design of nanophotonic on-chip devices for the harnessing of photon-phonon interactions. © 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4975002]


I. INTRODUCTION
][22][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 5 optical wave cycles, and resolving an SBS resonance with a linewidth of around 1 MHz would require at least 10 9 optical cycles.Thus, accurately treating photon-phonon interactions in time-domain simulations becomes prohibitively numerically expensive.][10][11][12][13]16,[18][19][20][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 acoustooptic 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 acoustooptic 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.

II. OPTICAL AND ACOUSTIC WAVE EQUATIONS
To start, we first briefly review the physics of Maxwell's equations [26][27][28][29] and the acoustic wave equation in the context of acousto-optic interactions. 8,25,30,31In 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 ω.J (ω) is the external current density, and P (ω) 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 Ũ (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, c − 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 30 where ∇ ⊗ U describes the tensor derivative of U as 020801-3 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 P (ω) 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,33In what follows, we show that one can explicitly treat the effects of P (ω) 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 P (t) has two contributions (we ignore the moving polarization effect because it is a much weaker effect as noted in Ref. 8), P(t) = P(b, PE) (t) + P(s,MB) (t), (7)   where P(b,PE) (t) is the bulk polarization density from photo-elasticity (PE) and P(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 χ (PE) ij (t) caused by the acoustic wave, 8 or written as a contracted tensor, where p − 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 32 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 D = ε 0 ε r E. By treating ∆E and ∆D 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 F (b,ES) is the bulk electrostrictive (ES) force, F (s,ES) is the surface electrostrictive force, and F (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 8 With the tensor contraction notation, we can express, in the vectorial form, the sum of the forces as 020801-5 Shi, Cerjan, and Fan APL Photonics 2, 020801 (2017) 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.

III. ACOUSTO-OPTIC FDFD FORMALISM
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 selfconsistent 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 E m ≡ E (ω 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 34 to iteratively compute the self-consistent solution {E m , U}.8][39][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 Ô ∈ C 2(M+1)×2(M+1) is a block-diagonal operator that acts linearly on the fields 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 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 Ô 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 ∂C/∂v, 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 34 to iteratively solve for the self-consistent solution for which g(v) = 0. Given the initial condition v 0 = 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: 34 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) 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 g (v) = Ô 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. 41When 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 v ∈ R 2K(M+1) contains all discretized unknowns in Eq. (31).In doing so, the Jacobian in Eq. ( 27) is a well-defined R 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 43 for smaller problems or with various iterative methods such as the bi-conjugate gradient method 44 and the quasi-minimal residual method 45 for larger systems.

IV. SIMULATION VERIFICATION AND EXAMPLE
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 ŷ, and the infinite out-of-the-plane direction to be ẑ.2][3][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 vector q, 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 U = û (y) e −iqx+iΩ B t , 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,11The optical and acoustic parameters of the waveguide core material are chosen to be those of chalcogenide glass, As 2 S 3 . 46Optically, 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 = 3200 kg/m 3 , and in the Voigt notation, its stiffness tensor is [c 11 , c 12 ] = [18.7,6.1] GPa, and the viscosity tensor is η 11 , η 12 = [1.8,1.45] mPa s. 46 Here, we assume that As 2 S 3 is an isotropic material, which implies that c 66 = (c 11 − c 12 )/2 and η 66 = (η 11 − η 12 )/2. 30This 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,8The 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 12 .In the As 2 S 3 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 12 = 0.24. 46The core is surrounded by vacuum, which APL Photonics 2, 020801 (2017) has a stiffness tensor of 0 and enforces the traction-free boundary conditions. 30To 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 ŷ.
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 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, 41 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 = 25 nm, 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,47A 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]. 48In Fig. 3(d), the powers of the Stokes field P 1 and acoustic field P a from both FDFD and CMT results are plotted along the waveguide direction for various pump powers P 2 .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 100 W/µm is likely due to the breakdown of the slowly-varying envelope-approximation in CMT. 8,48ext, 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 acoustooptic 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, 22 and it may also be used as a microwave photonic filter. 17,18The material of the waveguides is again chosen to be chalcogenide glass As 2 S 3 , 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,47n the absence of the acousto-optic interactions, the ring resonator is critically coupled to the external waveguide at λ res = 1558.29 nm (ω res = 2π × 192.4 THz) with a quality factor of Q = 6.25 × 10 3 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,14o demonstrate the SBS gain from this structure, we apply a photo-elastic coefficient of p 12 = 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 (P out /P in ) from this acoustooptic 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 2 S 3 . 46For 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. 34here 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. 49This 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.

V. DISCUSSION AND SUMMARY
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 firstprinciples frequency-domain techniques such as the finite element method (FEM), where a different discretization scheme is used. 50In 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,51n 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 biconjugate gradient 44 or quasi-minimal residual methods 45 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. 52Furthermore, this algorithm is not designed to handle thermally generated phonons with a broadband of frequency components. 4To 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,53t 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. 54Both 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.

SUPPLEMENTARY MATERIAL
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.

FIG. 1 .
FIG. 1. Schematics of the deformed waveguide geometry.The core has a relative permittivity of ε a , and it is surrounded by the cladding with a relative permittivity of ε b .

FIG. 2 .
FIG. 2. (a) Schematic of the waveguide structure with the electrostrictive region indicated.The Stokes wave is sent from left to right, and the pump is sent from right to left.(b) Dispersion relation of the optical waveguide.The blue and red points indicate the pump and the Stokes waves, respectively.The purple arrow represents the acoustic mode that couples the pump and Stokes waves.The insets on the right hand side show the acoustic and optical mode field patterns supported by the waveguide.

FIG. 3 .
FIG. 3. (a) Normalized plot for the E z field in the back-propagating pump wave.(b) Electric field pattern of the forwardpropagating Stokes waves at a pump power of P 2 = 100 W/µm.The optical field is clearly amplified as it propagates in x.(c) Field patterns of the acoustic mode inside the waveguide.(d) Comparison between the acousto-optic FDFD method with the solution from coupled-mode theory (CMT) under various pump powers.(e) Plot of maximum field amplitude that exists in each frequency sideband.The field amplitudes decrease drastically as we deviate from the pump and Stokes frequencies.(f) Plot of the relative error with number of iterations.The algorithm converges in only four iterations.

FIG. 4 .
FIG. 4. (a) Left: schematics of the ring resonator geometry, as well as the input locations of the pump (blue) and Stokes (red) waves.Right: without acousto-optic interactions, the ring resonator is critically coupled with the external waveguide at λ res = 1558.29 nm with a Q factor of 6.25 × 10 3 .(b) Left: a pump wave with a guided power of P 2 = 90 mW/µm is sent into the ring and is amplified inside the ring.Right: as a Stokes wave counter-propagates against the pump, it is amplified at the output.(c) Field plots of the acoustic modes generated from the interaction of the pump and Stokes waves.The fields show that the waveguide is highly multi-moded at the SBS frequency Ω B .(d) Plot of the power amplification experienced by the Stokes field as its frequency deviates from Ω B .The plot shows an ultra-narrow linewidth in the gain spectrum with a Fano-like line shape.