We introduce a finite-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 difficulty 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 find 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 algorithm should facilitate the design of nanophotonic on-chip devices for the harnessing of photon-phonon interactions.

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,17 and on-chip signal processing devices such as the microwave photonic filter18,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 105 optical wave cycles, and resolving an SBS resonance with a linewidth of around 1 MHz would require at least 109 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 equations26–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 as4,32

×μ01×𝐄(ω)ω2ε0εr𝐄(ω)ω2𝐏(ω)=iω𝐉(ω),
(1)

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

𝐔(t)=𝐔eiΩt+𝐔*eiΩt,
(2)

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

ρΩ2Ui+jklj(cijkl+iΩηijkl)kUl+Fi=0,
(3)

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 Ui and Fi 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 as30 

ρΩ2𝐔+(𝐜+iΩ𝜼¯):𝐔+𝐅=0,
(4)

where 𝐔 describes the tensor derivative of U as

𝐔=(xUxxUyxUzyUxyUyyUzzUxzUyzUz),
(5)

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

[𝐜:𝐔]ij=klcijklkUl,
(6)

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),

𝐏(t)=𝐏(b,PE)(t)+𝐏(s,MB)(t),
(7)

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,8 

χij(PE)(t)=εr2klpijklk(UleiΩt+Ul*eiΩt),
(8)

or written as a contracted tensor,

χ(PE)(t)=εr𝐩2:(𝐔eiΩt+𝐔*eiΩt),
(9)

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

𝐏(t)=ε0χ(PE)(t)𝐄(t),
(10)

which can be expressed in the frequency domain as

𝐏(b,PE)(ω)=ε0εr2[(𝐩:𝐔)𝐄(ωΩ)+(𝐩:𝐔*)𝐄(ω+Ω)].
(11)

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 form32 

𝐄(t)=m𝐄meiωmt+𝐜.𝐜.,
(12)

where Em 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

𝐏(b,PE)(ωm)=ε0εr2[(𝐩:𝐔)𝐄m1+(𝐩:𝐔*)𝐄m+1].
(13)

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,

Δ𝐄=(εb1εa1)ε01n^(n^𝐃),
(14a)
Δ𝐃=(εaεb)ε0(n^×n^×𝐄),
(14b)

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

𝐏(s,MB)(ωn)=ε0(εaεb)n^×[n^×𝐄m1(𝐔𝝈s)n^×𝐄m+1(𝐔*𝝈s)]+ε0(εb1εa1)n^[n^(εr𝐄m1)(𝐔𝝈s)+n^(εr𝐄m+1)(𝐔*𝝈s)],
(15)

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

×μ01×𝐄mωmε02εr𝐄mωmε02εr2[(𝐩:𝐔)𝐄m1+(𝐩:𝐔*)𝐄m+1]ωmε02(εaεb)n^×[n^×𝐄m1(𝐔𝝈s)n^×𝐄m+1(𝐔*𝝈s)]ωmε02(εb1εa1)n^[n^(εr𝐄m1)(𝐔𝝈s)+n^(εr𝐄m+1)(𝐔*𝝈s)]=iωm𝐉m(ωm).
(16)
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. 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.

Close modal

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,

𝐅=𝐅(b,ES)+𝐅(s,ES)+𝐅(s,MB),
(17)

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, as8 

Fl(b,ES)=ε0εr2ijkk(pijklm(𝐄m*)i(𝐄m+1)j),
(18a)
Fl(s,ES)=ε0εr2ijkσks(pijklm(𝐄m*)i(𝐄m+1)j).
(18b)

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

𝐅(b,ES)+𝐅(s,ES)=ε0εr2(𝝈s)𝐩:m𝐄m*𝐄m+1.
(19)

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

𝐅(s,MB)=𝝈s[ε0(εaεb)m(n^×𝐄m*)(n^×𝐄m+1)ε0(εb1εa1)m(n^εr𝐄m*)(n^εr𝐄m+1)].
(20)

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

ρΩ2𝐔+(𝐜+iΩ𝜼):𝐔+𝝈sεr2(σs)𝐩:m𝐄m*𝐄m+1+𝝈s[ε0(εaεb)m(n^×𝐄m*)(n^×𝐄m+1)ε0(εb1εa1)m(n^εr𝐄m*)(n^εr𝐄m+1)]=𝐅ext,
(21)

where we included an external force, Fext, 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 method34 to iteratively compute the self-consistent solution {𝐄m,𝐔}. The treatment below is similar in setup to the harmonic balance method for nonlinear circuit simulations35,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

𝐯=[𝐯1𝐯2𝐯M𝐯M+1𝐯M+2𝐯2M𝐯2M+1𝐯2M+2]T[𝐄1𝐄2𝐄M𝐄1*𝐄2*𝐄M*𝐔𝐔*]T.
(22)

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),

𝐠(𝐯)=O^𝐯+𝐂(𝐯)𝐛=0,
(23)

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

O^=diag[A1A2AMA1*A2*AM*BB*],
(24a)
Am=×μ01×()ωmε02εr,
(24b)
B=ρΩ2+(𝐜+iΩ𝜼):(),
(24c)

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

𝐛=[iω1𝐉1iωM𝐉Miω1𝐉1*iωM𝐉M*𝐅ext𝐅ext*]T,
(25)

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

𝐂(𝐯)=(𝐊(𝐯)𝐊(𝐯)𝐋(𝐯)𝐋(𝐯)).
(26a)

In Eq. (26a), K(v) describes M set of equations that govern the optical components {E1E2EM}, 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

𝐊m(𝐯)=ωm2(𝐏(b,PE)(𝐯)+𝐏(s,MB)(𝐯)),
(26b)

and L(v) can be written as

𝐋(𝐯)=𝐅(b,ES)(𝐯)+𝐅(s,ES)(𝐯)+𝐅(s,MB)(𝐯).
(26c)

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

D𝐠(𝐯)=𝐠(𝐯)𝐯=O^+𝐂(𝐯)𝐯,
(27)

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 Dg(v), we can apply the Newton-Raphson algorithm34 to iteratively solve for the self-consistent solution for which g(v) = 0. Given the initial condition v0 = 0, subsequent updates at the (k + 1)th step for vk+1 can be obtained as

𝐯k+1=𝐯k𝐬k,
(28)

where sk defines the step of the Newton-Raphson algorithm, computed by solving the following linear equation:34 

D𝐠(𝐯k)𝐬k=𝐠(𝐯k).
(29)

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

δ(k)=||𝐯k+1𝐯k||
(30)

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

𝐯=[Re{𝐄1}Im{𝐄1}Re{𝐄M}Im{𝐄M}Re{𝐔}Im{𝐔}]T
(31)

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.41 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 Dg(vk) is sparse, and Eq. (29) can thereby be very efficiently computed at each step with methods such as matrix factorization43 for smaller problems or with various iterative methods such as the bi-conjugate gradient method44 and the quasi-minimal residual method45 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 Ez, Hx, and Hy). 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 E1 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)eiqx+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

ρΩBu^i2+jkl(yy^iqx^)jcijkl(yy^iqx^)ku^l=0.
(32)

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, As2S3.46 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.46 Here, we assume that As2S3 is an isotropic material, which implies that c66=(c11c12)/2 and η66=(η11η12)/2.30 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.

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

Close modal

In this structure, the acousto-optic interaction occurs when the Ez field couples with the Ux and Uy 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 Ez, Ux, and Uy can be described by a single photo-elasticity element p12. In the As2S3 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 p12 = 0.24.46 The core is surrounded by vacuum, which has a stiffness tensor of 0 and enforces the traction-free boundary conditions.30 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μm1, there exists a backward acoustic mode at ΩB=2π×5.79 GHz and propagation constant q=15.08μm1 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 Ez at the origin of each Yee cell,41 and the Ux and Uy 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).

FIG. 3.

(a) Normalized plot for the Ez field in the back-propagating pump wave. (b) Electric field pattern of the forward-propagating Stokes waves at a pump power of P2=100W/μ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. 3.

(a) Normalized plot for the Ez field in the back-propagating pump wave. (b) Electric field pattern of the forward-propagating Stokes waves at a pump power of P2=100W/μ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.

Close modal

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].48 In Fig. 3(d), the powers of the Stokes field P1 and acoustic field Pa from both FDFD and CMT results are plotted along the waveguide direction for various pump powers P2. 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 lasers11,13,14 or as a nonreciprocal light storage unit,22 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 As2S3, 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

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.29nm with a Q factor of 6.25×103. (b) Left: a pump wave with a guided power of P2=90mW/μ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.

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.29nm with a Q factor of 6.25×103. (b) Left: a pump wave with a guided power of P2=90mW/μ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.

Close modal

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 As2S3.46 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.34 

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.49 This intriguing detail is typically neglected in the descriptions of experimental observations11,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.50 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 gradient44 or quasi-minimal residual methods45 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.52 Furthermore, this algorithm is not designed to handle thermally generated phonons with a broadband of frequency components.4 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.54 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.

1.
B. J.
Eggleton
,
C. G.
Poulton
, and
R.
Pant
, “
Inducing and harnessing stimulated Brillouin scattering in photonic integrated circuits
,”
Adv. Opt. Photonics
5
,
536
587
(
2013
).
2.
A.
Kobyakov
,
M.
Sauer
, and
D.
Chowdhury
, “
Stimulated Brillouin scattering in optical fibers
,”
Adv. Opt. Photonics
2
,
1
59
(
2010
).
3.
G.
Agrawal
,
Nonlinear Fiber Optics
, 5th ed. (
Elsevier
,
2013
).
4.
R.
Boyd
,
Nonlinear Optics
, 3rd ed. (
Academic
,
2008
).
5.
E.
Peral
and
A.
Yariv
, “
Degradation of modulation and noise characteristics of semiconductor lasers after propagation in optical fiber due to a phase Shift induced by stimulated Brillouin scattering
,”
IEEE J. Quantum Electron.
35
,
1185
1195
(
1999
).
6.
V.
Laude
and
J. C.
Beugnot
, “
Generation of phonons from elecrostriction in small-core optical waveguides
,”
AIP Adv.
3
,
042109
(
2013
).
7.
P. T.
Rakich
,
C.
Reinke
,
R.
Camacho
,
P.
Davis
, and
Z.
Wang
, “
Giant enhancement of stimulated Brillouin scattering in the subwavelength limit
,”
Phys. Rev. X
2
,
011008
(
2012
).
8.
C.
Wolff
,
M. J.
Steel
,
B. J.
Eggleton
, and
C. G.
Poulton
, “
Stimulated Brillouin scattering in integrated photonic waveguides: Forces, scattering mechanisms, and coupled-mode analysis
,”
Phys. Rev. A
92
,
013836
(
2015
).
9.
W.
Qiu
,
P. T.
Rakich
,
H.
Shin
,
H.
Dong
,
M.
Soljacic
, and
Z.
Wang
, “
Stimulated Brillouin scattering in nanoscale silicon step-index waveguides: A general framework of selection rules and calculating SBS gain
,”
Opt. Express
21
,
31402
31419
(
2013
).
10.
M.
Santagiustina
,
S.
Chin
,
N.
Primerov
,
L.
Ursini
, and
L.
Thevenaz
, “
All-optical signal processing using dynamic Brillouin gratings
,”
Sci. Rep.
3
,
1594
(
2013
).
11.
H.
Lee
,
T.
Chen
,
J.
Li
,
K. Y.
Yang
,
S.
Jeon
,
O.
Painter
, and
K.
Vahala
, “
Chemically etched ultrahigh-Q wedge-resonator on a silicon chip
,”
Nat. Photonics
6
,
369
373
(
2012
).
12.
E. A.
Kittlaus
,
H.
Shin
, and
P. T.
Rakich
, “
Large Brillouin amplification in silicon
,”
Nat. Photonics
10
,
463
468
(
2016
).
13.
W.
Loh
,
J.
Becker
,
D. C.
Cole
,
A.
Coillet
,
F. N.
Baynes
,
S. B.
Papp
, and
S. A.
Diddams
, “
A microrod-resonator Brillouin laser with 240 Hz absolute linewidth
,”
New J. Phys.
18
,
045001
(
2016
).
14.
W.
Loh
,
A. A. S.
Green
,
F. N.
Baynes
,
D. C.
Cole
,
F. J.
Quinlan
,
H.
Lee
,
K. J.
Vahala
,
S. B.
Papp
, and
S. A.
Diddams
, “
Dual-microcavity narrow-linewidth Brillouin laser
,”
Optica
3
,
225
231
(
2015
).
15.
D.
Braje
,
L.
Hollberg
, and
S.
Diddams
, “
Brillouin-enhanced hyperparametric generation of an optical frequency comb in a monolithic highly nonlinear fiber cavity pumped by a cw laser
,”
Phys. Rev. Lett.
102
,
193902
(
2009
).
16.
M. S.
Kang
,
A.
Nazarkin
,
A.
Brenn
,
P.
St
, and
J.
Russel
, “
Tightly trapped acoustic phonons in photonic crystal fibres as highly nonlinear artificial Raman oscillators
,”
Nat. Phys.
5
,
276
280
(
2009
).
17.
Y.
Okawachi
,
M. S.
Bigelow
,
J. E.
Sharping
,
Z.
Zhu
,
A.
Schweinsberg
,
D. J.
Gauthier
,
R. W.
Boyd
, and
A. L.
Gaeta
, “
Tunable all-optical delays via Brillouin slow light in an optical fiber
,”
Phys. Rev. Lett.
94
,
153902
(
2005
).
18.
A.
Byrnes
,
R.
Pant
,
E.
Li
,
D.-Y.
Choi
,
C. G.
Poulton
,
S.
Fan
,
S.
Madden
,
B.
Luther-Davies
, and
B. J.
Eggleton
, “
Photonic chip based tunable and reconfigurable narrowband microwave photonic filter using stimulated Brillouin scattering
,”
Opt. Express
20
,
18836
18845
(
2012
).
19.
H.
Jiang
,
D.
Marpaung
,
M.
Pagani
,
K.
Vu
,
D. Y.
Choi
,
S. J.
Madden
,
L.
Yan
, and
B.
Eggleton
, “
Wide-range, high-precision multiple microwave frequency measurement using a chip-based photonic Brillouin filter
,”
Optica
3
,
30
34
(
2016
).
20.
X.
Huang
and
S.
Fan
, “
Complete all-optical silica fiber isolator via stimulated Brillouin scattering
,”
J. Lightwave Technol.
29
,
2267
2275
(
2011
).
21.
C. G.
Poulton
,
R.
Pant
,
A.
Byrnes
,
S.
Fan
,
M. J.
Steel
, and
B. J.
Eggleton
, “
Design for broadband on-chip isolator using stimulated Brillouin scattering in dispersion-engineered chalcogenide waveguides
,”
Opt. Express
20
,
21235
21246
(
2012
).
22.
C. H.
Dong
,
Z.
Shen
,
C. L.
Zou
,
Y. L.
Zhang
,
W.
Fu
, and
G. C.
Guo
, “
Brillouin-scattering-induced transparency and non-reciprocal light storage
,”
Nat. Commun.
6
,
6193
(
2015
).
23.
J.
Kim
,
M. C.
Kuzyk
,
K.
Han
,
H.
Wang
, and
G.
Bahl
, “
Non-reciprocal Brillouin scattering induced transparency
,”
Nat. Phys.
11
,
275
280
(
2015
).
24.
A.
Taflove
and
S.
Hagness
,
Computational Electrodynamics: The Finite-Difference Time-Domain Method
, 3rd ed. (
Artech House
,
2005
).
25.
X.
Yuan
,
D.
Borup
,
J. W.
Wiskin
,
M.
Berggren
,
R.
Eidens
, and
S. A.
Johnson
, “
Formulation and validation of Berenger’s PML absorbing boundary for the FDTD simulation of acoustic scattering
,”
IEEE Trans. Ultrason. Eng.
44
,
816
822
(
1997
).
26.
N. J.
Champagne
 II
,
J. G.
Berryman
, and
H. M.
Buettner
, “
FDFD: A 3D finite difference frequency domain code for electromagnetic induction tomography
,”
J. Comput. Phys.
170
,
830
848
(
2001
).
27.
W.
Shin
and
S.
Fan
, “
Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell’s equations solvers
,”
J. Comput. Phys.
231
,
3406
3431
(
2012
).
28.
W.
Shin
and
S.
Fan
, “
Accelerated solution of the frequency-domain Maxwell’s equations by engineering the eigenvalue distribution of the operator
,”
Opt. Express
21
,
22578
22595
(
2013
).
29.
.
G.
Veronis
,
R. W.
Dutton
, and
S.
Fan
, “
Method for sensitivity analysis of photonic crystal devices
,”
Opt. Lett.
29
,
2288
2290
(
2004
).
30.
B. A.
Auld
,
Acoustic Fields and Waves in Solids
, 2nd ed. (
Krieger
,
1990
).
31.
J.-C.
Beugnot
and
V.
Laude
, “
Electrostriction and guidance of acoustic phonons in optical fibers
,”
Phys. Rev. B
86
,
224304
(
2012
).
32.
Y.
Shi
,
W.
Shin
, and
S.
Fan
, “
A multi-frequency finite-difference frequency-domain algorithm for active nanophotonic device simulations
,”
Optica
3
,
1156
1159
(
2016
).
33.
J. E.
Sipe
and
M. J.
Steel
, “
A Hamiltonian treatment of stimulated Brillouin scattering in nanoscale integrated waveguides
,”
New J. Phys.
18
,
045004
(
2016
).
34.
W. H.
Press
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
Numerical Recipes in C: The Art of Scientific Computing
, 3rd ed. (
Cambridge University Press
,
Cambridge, England
,
2007
).
35.
K. S.
Kundert
and
A.
Sangiovanni-Vincentelli
, “
Simulation of nonlinear circuits in the frequency domain
,”
IEEE Trans. Comput. Aided Des. Integr. Circuits Syst.
5
,
521
535
(
1986
).
36.
B.
Troyanovsky
,
Z.
Yu
, and
R. W.
Dutton
, “
Physics-based simulation of nonlinear distortion in semiconductor devices using the harmonic balance method
,”
J. Comput. Methods Appl. Mech. Eng.
181
,
467
482
(
2000
).
37.
H. E.
Tureci
,
A. D.
Stone
, and
B.
Collier
, “
Self-consistent multimode lasing theory for complex or random lasing media
,”
Phys. Rev. A
74
,
043822
(
2006
).
38.
L.
Ge
,
Y. D.
Chong
, and
A. D.
Stone
, “
Steady-state ab initio laser theory: Generalizations and analytic results
,”
Phys. Rev. A
82
,
063824
(
2010
).
39.
S.
Esterhazy
,
D.
Liu
,
M.
Liertzer
,
A.
Cerjan
,
L.
Ge
,
K. G.
Makris
,
A. D.
Stone
,
J. M.
Melenk
,
S. G.
Johnson
, and
S.
Rotter
, “
Scalable numerical approach for the steady-state ab initio laser theory
,”
Phys. Rev. A
90
,
023816
(
2014
).
40.
A.
Cerjan
,
Y. D.
Chong
, and
A. D.
Stone
, “
Steady-state ab initio laser theory for complex gain media
,”
Opt. Express
23
,
6455
6477
(
2015
).
41.
K.
Yee
, “
Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media
,”
IEEE Trans. Antennas Propag.
14
,
302
307
(
1966
).
42.
A.
Cerjan
and
A. D.
Stone
, “
Why the laser linewidth is so narrow: A modern perspective
,”
Phys. Scr.
91
,
013003
(
2016
).
43.
J.
Bunch
and
J.
Hopcroft
, “
Triangular factorization and inversion by fast matrix multiplication
,”
Math. Comput.
28
,
231
236
(
1974
).
44.
R.
Fletcher
, “
Conjugate gradient methods for indefinite systems
,” in
Numerical Analysis
, Lecture Notes in Mathematics Vol. 506 (
Springer-Verlag, Berlin
,
1976
), pp.
73
89
.
45.
R.
Freund
and
N.
Nachtigal
, “
QMR: A quasi-minimal residual method for non-Hermitian linear systems
,”
Numerische Math.
60
,
315
339
(
1991
).
46.
M. J. A.
Smith
,
B. T.
Kuhlmey
,
C. M.
de Sterke
,
C.
Wolff
,
M.
Lapine
, and
C. G.
Poulton
, “
Metamaterial control of stimulated Brillouin scattering
,”
Opt. Lett.
41
,
2338
--
2341
(
2016
).
47.
J. P.
Berenger
, “
A perfectly matched layer for the absorption of electromagnetic waves
,”
J. Comput. Phys.
114
,
185
200
(
1994
).
48.
H.
Haus
,
Waves and Fields in Optoelectronics
(
Prentice-Hall
,
1984
).
49.
S.
Fan
,
W.
Suh
, and
J. D.
Joannopoulos
, “
Temporal coupled-mode theory for the Fano resonance in optical resonators
,”
J. Opt. Soc. Am. A
20
,
569
572
(
2003
).
50.
J. M.
Jin
,
The Finite Element Method in Electromagnetics
, 2nd ed. (
Wiley-IEEE Press
,
2002
).
51.
T. A.
Davis
, “
Algorithm 832
,”
ACM Trans. Math. Software
30
,
196
199
(
2004
).
52.
T. F.
Buttner
,
M.
Merklein
,
I. V.
Kabakova
,
D. D.
Hudson
,
D. Y.
Choi
,
B.
Luther-Davies
,
S. J.
Madden
, and
B. J.
Eggleton
, “
Phase-locked, chip-based, cascaded stimulated Brillouin scattering
,”
Optica
1
,
311
314
(
2014
).
53.
A.
Pick
,
A.
Cerjan
,
D.
Liu
,
A. W.
Rodriguez
,
A. D.
Stone
,
Y. D.
Chong
, and
S. G.
Johnson
, “
Ab-initio multimode linewidth theory for arbitrary inhomogeneous laser cavities
,”
Phys. Rev. A
91
,
063806
(
2015
).
54.
R.
Zecca
,
P. T.
Bowen
,
D. R.
Smith
, and
S.
Larouche
, “
Transformation optics simulation method for stimulated Brillouin scattering
,”
Phys. Rev. A
94
,
063818
(
2016
).

Supplementary Material