A Hamiltonian treatment of stimulated Brillouin scattering in nanoscale integrated waveguides

We present a multimode Hamiltonian formulation for the problem of opto-acoustic interactions in optical waveguides. We establish a Hamiltonian representation of the acoustic field and then introduce a full system with a simple opto-acoustic coupling that includes both photoelastic/electrostrictive and radiation pressure/moving boundary effects. The Heisenberg equations of motion are used to obtain coupled mode equations for quantized envelope operators for the optical and acoustic fields. We show that the coupling coefficients obtained coincide with those established earlier, but our formalism provides a much simpler demonstration of the connection between radiation pressure and moving boundary effects than in previous work [C. Wolff et al, Physical Review A 92, 013836 (2015)].

We present a multimode Hamiltonian formulation for the problem of opto-acoustic interactions in optical waveguides. We establish a Hamiltonian representation of the acoustic field and then introduce a full system with a simple opto-acoustic coupling that includes both photoelastic/electrostrictive and radiation pressure/moving boundary effects. The Heisenberg equations of motion are used to obtain coupled mode equations for quantized envelope operators for the optical and acoustic fields. We show that the coupling coefficients obtained coincide with those established earlier, but our formalism provides a much simpler demonstration of the connection between radiation pressure and moving boundary effects than in previous work [C. Wolff et al., Physical Review A 92, 013836 (2015)].

I. INTRODUCTION
Almost a century after it was first proposed [1,2] and fifty years since the invention of the laser allowed its first observation [3], the phenomenon of stimulated Brillouin scattering (SBS) may only now be entering its golden age. At its simplest, SBS refers to the stimulated interaction between a pair of coherent optical waves and a resonant hypersonic acoustic wave. SBS has traditionally been encountered as the scattering of an optical pump beam into a backward traveling Stokes beam of slightly lower frequency by an acoustic wave oscillating at the optical beat frequency. The acoustic wave is generated by the process of electrostriction [4,5], and both the Stokes and acoustic wave grow by a process of positive feedback. In optical fiber, this process can be highly efficient and is often described as the "strongest" fiber nonlinearity. This can be problematic, as SBS prevents the propagation of high power narrow-bandwidth pumps. Nevertheless, SBS in fibers has long provided a mechanism for producing narrow linewidth lasers and amplifiers [6][7][8], filters and other spectral components for microwave photonics [9][10][11][12], as well as various sensors [13,14].
Like most optical nonlinearities, the emergence of sub-wavelength scale waveguides with strong confinement and tunable dispersion has greatly increased the efficiency, utility and reach of stimulated Brillouin processes. This includes the generation of frequency combs by cascaded Brillouin generation in microstructured small-core fibers in both backward [15] and forward [16] configurations, slow and fast light effects [17], and novel microstructured fiber lasers [18,19].
Motivated by such studies, the development of on-chip SBS in highly-nonlinear integrated waveguides has recently been pursued aggressively [20]. Attaining efficient on-chip SBS is complicated by the requirement of simultaneous confinement of both the optical and acoustic fields. This is non-trivial because optically dense materials suitable for optical waveguide cores are commonly mechanically stiff and therefore susceptible to leakage of the acoustic wave into softer substrates; the silicon on silica system is an important example. Consequently, on-chip SBS was first achieved [21] in rib waveguides made from nonlinear chalcogenide glasses, which combine high refractive index and nonlinearity with relative mechanical softness. Subsequently, SBS in silicon waveguides has been observed in Si/SiN membranes [22] and elevated rails [23,24], which both exploit physical isolation of the waveguide to minimize acoustic losses. Considerable development will be needed to reach designs suitable for mass-fabrication, but a practical platform for on-chip SBS would enable numerous applications [25] in microwave photonics [10,12,[26][27][28][29], sensing, isolators [30,31] and chip-based lasers [32,33].
A key driver for developing SBS in sub-micron waveguides was the realization by Rakich et al. [34,35] that at small scales, there are new contributions to SBS associated with radiation pressure of light on the waveguide boundaries, and the back-action of "moving boundaries" on the optical field [36]. Depending on the particular waveguide configuration and combination of optical modes, these contributions can either reinforce or counteract the more familiar bulk contributions from electrostriction and photoelasticity [34]. As well, waveguides in which the acoustic fields are strongly confined can enhance the scattering efficiency of near-stationary quasi-transverse acoustic waves, a requirement for efficient forward SBS where the pump and Stokes wave co-propagate [16,35].
Following this realization there was some variation in the literature as to how best to incorporate the new effects into a coupled mode theory self-consistently. In conventional SBS, electrostriction (the mechanical stress induced by the optical field) is accompanied by the complementary process of photoelasticity (the change in the dielectric response induced by the acoustic strain) [4]. The two processes are captured by identical coupling terms in the coupled mode theory, as required by the Manley-Rowe relations [23,37]. One should expect the same symmetry between the effect of radiation pressure on the waveguide boundaries driving the acoustic field, and the reverse effects of moving boundaries on the optical field. However, the initial formulations in terms of optical forces and the Maxwell stress tensor led to some confusion about whether certain additional coupling terms arise or whether they are essentially "double-counting". These include concepts of an electrostrictive "boundary pressure" [22] and bulk contributions to the radiation pressure term [34,35].
Recently, our group provided a new derivation of the coupled mode equations [37] that avoids the formalism of optical forces. Instead we used thermodynamic arguments to unambiguously identify the correct coupling term describing both radiation pressure and moving boundaries effects. We found that this term indeed matches expressions suggested in several papers [34,38] and the matter seems to be resolved. It turns out, for instance, that the appearance of an electrostrictive boundary pressure depends on whether electrostriction is viewed in terms of a stress or a force density [37]. Nevertheless, the argument establishing the correct form for the moving boundary coupling was quite involved, and a simpler, more direct derivation would be desirable.
In this work, we provide such a derivation. Rather than the standard approach of applying slowly-varying envelope approximations to the wave equation, we extend a quantized multimode Hamiltonian formalism of integrated optical waveguides [39] to include opto-acoustic interactions. The photoelastic and radiation pressure couplings are introduced through a single interaction energy term in the Hamiltonian, and fully quantum equations of motion are obtained from the Heisenberg equations. In the classical limit, the coupled mode equations of Wolff et al. [37] emerge naturally and simply, with no ambiguity about double-counting.
There are number of other advantages to our approach. The fundamental quantum process underlying SBS-the stimulated decay of a pump photon into a lower energy Stokes photon and an acoustic phonon-is manifestly visible in the interaction term. Further, obtaining quantum equations that respect the appropriate operator commutation relations provides an important starting point for the investigation of effects at the boundary of quantum and classical opto-acoustics. This is likely to become more important as the distinction between cavity optomechanics and guided wave opto-acoustics becomes increasingly blurred [40,41], and phonon confinement strategies are improved. The propagation of guided wave acoustic fields with strongly modified phonon density of states is likely not far off, which raises the prospect of Brillouin interactions in the quantum regime. Finally, the Hamiltonian formalism has proved very powerful for the description of other integrated quantum nonlinear processes such as spontaneous four wave mixing [42,43] and spontaneous parametric downconversion [44].
We should note that quantum or analytical dynamics approaches to guided wave opto-acoustics themselves have some pedigree. Hamiltonian approaches to SBS date back to the first rigorous treatment by Shen and Bloembergen [5] who gave an analysis for plane waves in the semi-classical limit. Drummond and Corney incorporate Raman gain into their quantized theory of nonlinear fiber propagation [45]. In that case, the Raman response, which depends on the detailed glass composition and network, is introduced through a phenomenological measured response function. In contrast, for the Brillouin couplings we consider here, the phonon response is entirely determined by the bulk elastic properties and waveguide geometry, and so can be calculated exactly using acoustic mode solvers. Finally, van Laer et al. [40] have recently discussed the connections between quantum optomechanics for single or few resonator systems and classical SBS [40]. They identify an elegant connection between the opto-acoustic coupling in waveguide SBS and the corresponding coupling in quantum optomechanical systems. The latter is treated with a single mode Hamiltonian approach which is appropriate for the optomechanics of a resonator consisting of a single cavity, but limits its application to longer structures with continuous phonon spectra. In contrast, for the waveguide problem that is our focus, a full multi-mode treatment is appropriate in order to treat arbitrary input optical fields.
The paper is structured as follows. In Section II we construct a Hamiltonian description of acoustics, including useful expressions for group velocity and power flow. In Section III we review some necessary results from guided wave electromagnetic quantization. In Section IV, which is the core of the paper, we provide the full optoacoustic Hamiltonian, find expressions for the coupling terms, and show directly how the symmetry between radiation pressure and moving boundary effects emerges. In Section V we derive quantum coupled mode equations for the system, make connections to prior expressions for the coupling strength, and recover the classical coupled mode equations for SBS. Finally, in Section VI we discuss directions for future examination including the important issue of phonon dissipation. A comprehensive supplementary materials document provides detailed derivations of many of the results.

II. HAMILTONIAN FORMULATION OF GUIDED WAVE ACOUSTICS
The classical theory of guided elastic waves is of course very mature and can be formulated in many guises. Auld [46] provides an excellent introduction for readers with an optics background. Since our goal is a Hamiltonian operator representing the complete opto-acoustic system we begin by re-framing guided acoustic wave propagation in a quantum Hamiltonian picture; we have not found such a formulation in the literature.
For purely classical applications, one might build a Hamiltonian from which the dynamics are determined by Hamilton's equations. For generality, we construct the theory in a quantum form, using canonical quantization with commutators that follow from the standard association with the Poisson brackets of the classical formulation:

A. Hamiltonian operator
To identify a classical theory suitable for canonical quantization we should begin with canonical variables, which will become the canonical operators in the quantum theory, and a classical Hamiltonian that both yields the standard equations of motion in the form of the elastic wave equation, and is numerically equal to the classical energy of the system. To that end we introduce vector field variables u(r) describing the displacement and π(r) as their conjugate momenta; these will become operators in the quantum theory, although we will not explicitly include "hats" in our notation.
Following standard quantum mechanics, we naturally choose the commutation relations where the integration is over all space. Here ρ(r) and c ijkl (r) are position-dependent c-number quantities representing the density and stiffness tensor respectively, while is the strain tensor operator [46], and again we emphasise that unless stated otherwise, u(r) and π(r) are to be read as operators. In (3) and throughout, repeated indices are to be summed over the Cartesian coordinates x, y, z. We will formally assume these quantities are continuous functions of position, although they may change drastically as one moves from solid to air, for example. The limit to these functions changing discontinuously can be taken at the end of the calculation when matrix elements and the like are evaluated, but in all dynamical equations and derivations ρ(r) and c ijkl (r) should be taken as continuous functions. Let us comment on the precise meaning of the operators u and π. In elastic theory, the displacement field applies to "volume elements" much larger than the atomic scale, but much smaller than the characteristic wavelength of any acoustic excitation. Therefore, the displacement and momentum operators do not correspond directly to any individual physical oscillator but describe collective excitations of the mesoscopic bulk medium. Nevertheless the low-energy phonon excitations that emerge from the theory are very real and their quantization is physical. For example, completeness relations are correct up to an appropriate wavevector cutoff, far above the typical wavevector of any Brillouin-induced excitation.
Classically these equations indicate, as expected, that the conjugate momentum is the product of the mass density and the velocity fieldu n , and that the momentum evolves according to Newton's laws. From (8), we recover the standard acoustic wave equation [46]: Thus we have the desired equations of motion in both the classical and quantum descriptions. Moreover, substituting the relation (8a) into (3) and taking the classical limit by dropping the operator nature of u and S ij , we recover the classical Hamiltonian This is clearly numerically equal to the sum of kinetic and potential energy [46], and so (3) meets all the criteria to be an appropriate Hamiltonian for quantization. Note that at this point, no acoustic dissipation (viscosity) has been included. Direct incorporation of losses is less straightforward in Hamiltonian approaches than in Lagrangian approaches, in which a Rayleigh dissipation function can be introduced. The phonon loss is important to the physics of SBS, but can be incorporated later perturbatively. We discuss this in Section VI.

B. Modes and new fields
We now seek to reduce the Hamiltonian to standard harmonic oscillator form; the low energy states will describe the quasi-particle or phonon excitations of the medium. To diagonalise the acoustic Hamiltonian we look for classical solutions of (9) of the form u(r, t) = U Λ (r)e −iΩΛt + c.c., (10) where Λ is a mode index. This implies where again we have used the fact that c njkl (r) is symmetric in its last two indices. It is convenient however, if the linear operator generating mode functions is Hermitian, which the form in Eq. (11) is not. To obtain a Hermitian system, we introduce new fieldsũ which clearly preserve the commutation relations. That is, We now look for mode solutions in terms of the new fields in the form π(r, t) =Π Λ (r)e −iΩΛt + c.c.
We introduce the operator M nk which acts on a general vector function C(r) as Section S.II of the supplementary material shows that M nk (r) is indeed Hermitian. Then, using integration by parts we can write the Hamiltonian (3) in terms of the new fields (12) in the simple form As the basis of our acoustic modal expansion, Eq. (15) plays an analogous role to that which the vector Helmholtz equation (commonly known as the "master equation" in the photonic crystal literature [48],) plays in the quantization of the electromagnetic field (see (43) below). Its Hermitian form is useful both for developing the formalism but also for formulating mode-solving algorithms based on energy functional minimization [48]. Now if we construct eigenfunctions F Λ (r) of M nk such that it follows from (11) and (12) that we can takeŨ where the second equation follows from substituting (14) in (8a). From the Hermiticity of M nk , the eigenfunctions F Λ (r) with different eigenvalues will automatically be orthogonal, and if there are eigenfunctions with the same eigenvalue we can construct them to be orthogonal. Thus we can write Here the integration is over a finite volume that, in the end, can be allowed to pass to infinity if we wish to generate a continuum of eigenfunctions. We take the set of eigenfunctions with positive eigenvalues Ω 2 Λ as complete, at least for the problems of interest. For each Ω 2 Λ we choose a positive Ω Λ and take it as the frequency of the eigenfunction. Some additional properties of the mode functions that are required for reduction of the acoustic Hamiltonian to canonical harmonic oscillator form are developed in section S.III of the supplementary material.
Following similar arguments to those used in quantization of the vacuum electromagnetic field [49] (or the electromagnetic field in nondispersive media [39]), one can show using (20) that introducing new operators b Λ and b † Λ with commutation relations we can preserve the commutation relations (13), by expanding the fields as where h.c. is the Hermitian conjugate. Then, through a series of manipulations using properties of the F Λ (r), the Hamiltonian can be reduced to the canonical harmonic oscillator form where we have dropped the zero point energy which has no dynamical effect. Derivations of Eqs. (22) and (23) are provided in section S.IV of the supplementary material. Finally, returning to the physical modes by introducing we can write the full displacement and momentum field operators as with prefactors in the expansions reminiscent of simple harmonic oscillator physics.

C. Waveguide acoustics
The results to this point apply to any acoustic structure. We now specialize to waveguides running along the z direction, and choose a box of length L in that direction, which includes all x and y. Focusing on acoustic modes confined to the waveguide, we label the modes by a wavenumber q = 2πn/L, for integer n, and a band α which identifies the transverse spatial mode structure in the xy plane. Then the eigenfunctions of the previous section can be written with the normalization (19) guaranteed by requiring where the integration is over the whole x-y plane. The derived mode amplitudes (24) are then of the form Note that the normalization condition (27) can be written, using (28a), as If we now let L → ∞, moving to a continuous distribution of modes, the commutation relations, Hamiltonian and field operator expansions respectively become Here the q integrals are over the range of q for which each mode α exists, taking account of any modal cutoffs. Normally we work in the Heisenberg regime so the b αq are time-dependent. Note our convention that while u and π are operators, the mode functions u αq and π αq , which carry modal index subscripts, are c-number quantities.

D. Envelope functions
To make the connection to the more familiar waveguide representation of slowly-varying envelopes, we now introduce envelope functions for each type of acoustic mode. We assume the excitation is centered at some wavenumber q o and factor out that dependence to produce a function varying slowly in space. However we retain the full time-dependence in the operators b αq (t) writing If the integral in (32) is taken to range over all q it is easily checked that we obtain the canonical equal-time continuous commutators In reality, the range of integration in (32) is restricted by modal cutoffs, which will temper the Dirac delta function in (33). However, assuming spectrally narrow envelopes far from any cut-offs, Eq. (33) should normally be an excellent approximation. From (31) we see that if Ω αq , u αq (x, y), and π αq (x, y) vary little over the range of significant q we can approximate the field expansions in terms of envelope functions centered at q = q o : Corrections can be included by expanding the prefactors Ω αq /2 u αq (x, y) and Ω αq /2 π αq (x, y) about q = q o , and combining the resulting powers of (q − q o ) with the exp(i(q − q o )z) in the integrals over q to yield expressions involving derivatives of φ α (z, t), but we neglect them here.
Returning to Eq. (32), we can derive approximate equations of motion for the φ α (z, t) in the Heisenberg picture by expanding the dispersion relation Ω αq of mode α about q o : where etc. An expression for the group velocity v αq in terms of the modal fields is worked out in section S.V of the supplementary material.

E. Acoustic powers
Finally, we establish some expressions concerning the power carried by the acoustic modes in the envelope representation. Classically, the acoustic power density at a point in the medium in a directionn is given by which has the natural interpretation of power being the dot product of an applied force and the velocity of the point of application [46]. Using (8), the power carried by a waveguide mode in theẑ direction is thus To construct a quantum operator corresponding to the power density we use the symmetrized form of the noncommuting operators π(r) and S kl (r), in the usual way, as shown in section S.VI of the supplementary material. This leads to the result that the operator for the power carried by the acoustic field is Here the contribution from the pair of modes α and α at wavenumbers q and q is given by where we have used (28) to introduce the modified mode functions f α,q .
For a phonon field involving only one mode α and assuming we can neglect the q dependence of p αα (q, q ) over the pulse spectrum we obtain the slowly-varying power operator as Finally, it may be shown (see section S.VI of the supplementary material), that p A αα (q o , q o ) = Ω αqo v αqo , so that the power carried by the acoustic envelope is and it follows that in this limit, φ † α (z)φ α (z) has the natural interpretation of a phonon number density operator.

III. QUANTIZATION OF THE ELECTROMAGNETIC FIELDS
To construct the full opto-acoustic Hamiltonian we will need similar results for the quantization of the electromagnetic field in integrated structures. The procedure is well known and we simply summarize some essential results [39].

A. Hamiltonian and modes
As the fundamental quantum fields we take the electric displacement D (r) and magnetic B (r) fields with commutation relations This choice, which dates back to Born and Infeld [50] has the advantage that the transversality of the two fields is easily imposed. The quantization procedure has been discussed at length elsewhere [39], and we simply quote the necessary. The parallels to the acoustic problem in the previous section are very apparent. The electromagnetic Hamiltonian operator is taken as where describes the "background" dielectric response of the waveguide structure in terms of the relative dielectric constant ref (r), without any acoustic effects. Note that it is straightforward to include a tensor response in ref (r) but to reduce cluttering the tensor notation, here we treat the material as optically isotropic. In principle, we could also extend the treatment to include dispersion of the dielectric [51] at the expense of considerably more complexity. For Brillouin processes, the linewidths of the interacting optical waves are usually narrow, and we can safely neglect dispersion within each optical field. As with the acoustic case, we are interested in optical waveguide modes with translational invariance along z. We find these modes {B Λ (r), ω Λ } by solution of the vector Helmholtz equation together with Ampere's law and then introduce waveguide mode functions d γk defined by where γ indexes the transverse spatial bands of the waveguide. We choose the normalization and expand the displacement field operator as The mode operators a γk satisfy the standard commutation relations and neglecting the vacuum energy, the electromagnetic Hamiltonian reduces to

B. Envelope operators
As with the acoustic modes we can introduce envelope function operators associated with a mode γ and a range of wavenumbers in the neighborhood of some k j : The integration is to be taken over the range of wavenumbers k that we wish to associate with the center value k j . This allows the introduction of distinct envelope operators for fields that occupy the same spatial mode but occupy distinct frequency ranges (such as a pump and Stokes wave in the same mode). If the integrals in (49) were to extend over all k we would obtain In principle, the existence of cutoffs and the possible partitioning of each channel γ into separate bands for pump and Stokes waves means that the k integrals have restricted range, but as with the acoustic fields, we assume the excitations are sufficiently narrow band and away from cutoff that the integrals leading to the Dirac delta function in (50) can be safely extended to infinity. Assuming in fact that only values of k close to k j are important for each mode γj, as we now label them, we can write where we have put ω j γ ≡ ω γkj and neglected the variation in ω j γ and the k dependence of d γk (x, y); as in our treatment of acoustic fields, corrections to these expressions can be easily identified.
Again in analogy with the treatment of the acoustic fields, an operator for the slowly-varying part of the power in the waveguide can be constructed from the Poynting vector (section S.VII of the supplementary material) which takes the form For ranges of k and k close enough to k j , and assuming that for different γ the corresponding k j ranges are distinct, we can write this as and since we find that where v j γ is the group velocity of electromagnetic mode type γ centered at k j , we have so that ψ † γj (z, t)ψ γj (z, t) behaves as a photon number density operator; note again the similarity with the treatment of the acoustic field, for which (40) is the corresponding result. The form in Eqs. (54) and (55) does not seem to have been presented earlier, and is derived in section S.VII of the supplementary material.
Finally, in similar fashion to Eqs. (34), we find that the envelope operator obeys the dynamical equation where are respectively the group velocity and group velocity dispersion of mode γj at its reference wavenumber k j . For the narrow bandwidths involved in SBS physics, higher dispersive terms are unlikely to be needed.

IV. THE COMPLETE OPTO-ACOUSTIC HAMILTONIAN
At last, we can now assemble the complete opto-acoustic Hamiltonian Being composed of different classes of oscillators, the electromagnetic and acoustic fields commute with each other. All quantities are taken as position dependent, varying continuously (if rapidly) across any material boundaries. Only at the end of various calculations we will allow them to acquire step-wise discontinuities. The opto-acoustic coupling is captured by the new quantity β ij (r). This is the total inverse (relative) dielectric tensor, which includes both the purely electromagnetic properties (the background waveguide structure) in β ref , and the photoelastic and radiation pressure couplings in the correctionβ ij . Naturally, we have β ij (r) jk (r) = δ ik , where jk (r) is the complete relative dielectric tensor. The coupling between sound and light enters because we assume that β ij (r) depends on the displacement field u (r), through both its dependence on the strain in the material and the motion of the interfaces. Thus we can write where the opto-acoustic coupling is In principle, an additional coupling arises from the dependence of the mechanical density and stiffness on the electromagnetic field variables. These effects lead to terms quadratic in u(r) or π(r) (corresponding to two-phonon-singlephoton interactions) which are of higher order than we consider here. They are also of much lower energy and the processes are very unlikely to be phase-matched, so they are safely neglected.
To proceed to a set of coupled mode equations, we seek to expand the interaction Hamiltonian (59) in terms of the mode operators constructed earlier. The physics of the opto-acoustic interaction is introduced by writing keeping only linear terms in the strain in keeping with our neglect of two-photon interactions. The photoelastic tensor p ijlm (x, y) accounts for the conventional electrostrictive/photoelastic contribution to SBS. The effect of moving boundaries and radiation pressure enters through the second term's dependence on the displacement u(r). Note that while we focus below on the effects associated with material discontinuities, this expression also accounts for a bulk contribution to the radiation pressure in graded index materials for which β ref (r) varies smoothly in space. Equation (60) is a key expression because the symmetric relationship of radiation pressure and moving boundary effects follow directly from its form. This identification is what allows us to avoid the rather lengthy thermodynamic arguments which were required in the previous rigorous derivation of the classical coupled mode equations [37]. Expanding β ref (r − u(r)) for small displacements we take and using (57) we have for the opto-acoustic correctioñ Then interaction (59) becomes Now since the phonon energy is much smaller than the photon energy the only significant terms in V will involve the creation and annihilation of photons. On substituting (47) for D, we normal order the photon mode operators that arise and neglect the resulting constant terms corresponding to vacuum fluctuation corrections to both the photoelastic tensor and the displacement-induced change in the dielectric properties. We thus find From (4) and the second of (31) we can write where S lm αq (r) is of the form After some manipulation (see section S.VIII of the supplementary material), the interaction can be reduced to the form where the coupling parameter is Note that if there is a slow variation of the nonlinear properties, due to longitudinal variation in, say, the composition or waveguide dimensions, then Γ(γk; γ k ; αq) will acquire this variation too, and the integration over z in (68) would capture this effect.
For an infinite homogeneous waveguide we can do the remaining integral over all z in (68) to obtain a delta function δ(k − k + q). Using this to eliminate the q integral in (68), the total Hamiltonian (58) becomes Observe that the final two terms explicitly display momentum conservation, and are clearly identified as describing anti-Stokes and Stokes processes respectively. This infinite structure form is a good starting point for investigating the enhancement and suppression of Brillouin scattering by adjusting the matrix elements or the density of states of optical or acoustic modes [52]; a simple Fermi's Golden Rule calculation reveals much of the underlying physics, as we will show in a subsequent contribution [53]. Technically, the finite phonon lifetime requires the δ(k − k + q) function to be broadened into a linewidth function. In practice the δ function is a reasonable approximation, with the impact of loss entering through the linear properties as discussed later.

V. QUANTUM COUPLED MODE EQUATIONS
To derive coupled mode equations for the envelope function operators, we assume that the process of interest destroys pump photons with transverse mode and center wavenumber (γ P , k P ), creating Stokes photons with (γ S , k S ), and phonons with mode-wavenumber pair (α, q) = (α o , q o ). Additional processes could be included to describe cascaded SBS phenomena [16,54]. We assume that the center wavenumbers and frequencies involved satisfy energy conservation and phase matching; that is, where we note the wavevectors can be positive or negative, and we have defined ω P ≡ ω PkP , ω S ≡ ω SkS , and Ω o = Ω αoqo . Figure 2 indicates the significance of these quantities for the processes of backward, forward intermodal and backward intermodal SBS. Then we can write (68) as V = dz dkdk dq (2π) 3/2 a † Pk e −i(k−kP)z a Sk e i(k −kS)z b αoq e i(q−qo)z Γ(γ P k; γ S k ; α o q) + dz dkdk dq (2π) 3/2 b † αoq e −i(q−qo)z a † Sk e −i(k −kS)z a Pk e i(k−kP)z Γ * (γ P k; γ S k ; α o q). Given the narrow bandwidths associated with SBS, we treat the coupling strengths Γ(γ P k; γ S k ; α o q) as constant over the range of wavenumbers integration: Pulling this term out, moving to the Heisenberg picture and using the definitions of the envelope functions in Eqs. (32) and (49), we reach where ψ S , ψ S , φ are the pump and Stokes photon fields and phonon fields respectively. Using (34) and (56), we find and similarly To remove the fast time-dependence, we define so that , and Φ † (z, t)Φ(z, t) can be identified as the power flowing in the waveguide in each mode at position z; this is a positive quantity with the direction along z given by the sign of the velocity.
Defining a reduced coupling constant and evaluating it at the center frequencies of our equations: we obtain the Heisenberg evolution equations for the field operators We have thus established a fully quantum form of the opto-acoustic coupled evolution equations and shown that the correspondence of coupling constants for photoelasticity plus moving boundaries and the reverse processes of electrostriction plus radiation pressure emerge naturally from the starting point of (60). We next evaluate the coupling constants (or "matrix elements") to connect them to familiar forms in the literature.

A. The coupling constants
We separate the coupling constantsΓ(γk; γ k ; αq) in Eq. (74) into the contributions Γ(γk; γ k ; αq) =Γ bulk (γk; γ k ; αq) + Γ surf (γk; γ k ; αq), The label "surface" for Eq. (79) is perhaps too restrictive since the derivative ∂β ref /∂r j is non-zero in graded index materials, however this contribution is typically weak, and in current experiments it is the surface contribution that is of most interest. We seek expressions in terms of the optical modes d γk and acoustic modes u αq that can be obtained from numerical solvers.

The bulk coupling constant
The bulk term in (78) describes photoelasticity (in (76a) and (76b)) and electrostriction (in (76c)). Its evaluation is straightforward. Even in the limit of a step discontinuity in material parameters across a waveguide interface, the d γk (x, y) will suffer at most a step discontinuity, as will s lm αq (x, y) and, by assumption, p ijlm (x, y). Hence (78) remains well-defined for step discontinuities.
Using the definition of s lm αq (x, y), one can show that for the standard interaction with modes (P, S, o), the bulk element reduces to the form which is consistent with the results in Wolff et al. [37] for the photoelastic coupling, there termed Γ ePE .

The surface matrix element
The surface term in (79) describes driving of the optical fields by moving boundaries (in (76a) and (76b)) and radiation pressure (in (76c)). Its evaluation is more subtle, since in the limit of a step discontinuity in material parameters β ref (x, y) will change in a step-like fashion, but its derivative can be Dirac delta-function-like. Consequently, a smoothing operation is required to make sense of this term. Following Johnson et al. [55] one can show (see section S.IX of the supplementary material) that for a boundary contour R c (s) parameterized by arc length s that separates two materials with dielectric constants + and − the surface contribution reduces tō Here the unit normaln(s) points in the direction from − to + . The full expression for Γ surf (γk; γ k ; αq) then involves a sum over all such curves separating distinct dielectrics. Note there is no ambiguity in evaluating these terms, since d ⊥ γk (r) is continuous across a step discontinuity in β ref (x, y), as is e γk (r). For the SBS combination of modes (P, S, o), the surface term can also be written as which coincides with the expression of Wolff et al. [37]. A normalized form of these expressions convenient for working with numerical mode solvers is provided in section S.X of the supplementary material.

B. Recovery of classical coupled mode equations
Finally, the standard classical coupled mode equations [37] can be recovered from (76) by dropping the dispersion terms and taking mean values for the operators: We have introduced the acoustic loss α by hand following the expression of Wolff et al. [37]: where η ijkl is the viscosity tensor.

VI. DISCUSSION
Equations (76) represent a full quantum description of guided-wave optoacoustic interactions. Although it requires some preliminary derivations to identify the effective fields that are involved, our approach provides a direct derivation of both the photoelastic/electrostrictive and moving-boundaries/radiation pressure components of the SBS interaction. The latter has clear contributions from both surface effects at material boundaries and bulk effects due to smooth variation in dielectric properties. By avoiding any discussion of forces and stress tensors, the ambiguities and challenges of prior treatments do not arise.
The equations of motion have a number of potential applications. In the classical limit they provide a rigorous confirmation of the earlier treatment [37]. In the quantum regime, we have a theory of opto-acoustic interactions that faithfully represents the photon and number statistics including any non-classical behavior. To date, quantum acoustic effects in guided wave systems have not been observed due to the overwhelming thermal contribution to the phonon field, though this may change in the near future. However, we can certainly envisage mixed systems in which a classical coherent state phonon field interacts with non-classical photon states in order to transfer quantum information between different channels, and our treatment is ideal for studies at this quantum-classical boundary.
An obvious and significant extension to our work would be a complete treatment of the acoustic dissipation, which we plan to present in the future. To incorporate dissipation into a Hamiltonian picture, a natural approach would be to introduce a thermal bath of phonon oscillators rather than just the Brillouin excited mode. These modes would couple with the coherent modes of interest and with each other through a three-phonon collision term associated with anharmonicity in the phonon Hamiltonian. Tracing over the additional oscillators would lead to an effective dissipation on the preferred phonon mode. According to need or preference, one could derive a dissipative master equation for the reduced phonon density operator, or perhaps more usefully, a set of Heisenberg equations with Langevin noise terms associated with the loss. This kind of treatment would be particularly important in understanding the impact of phonon loss on the photon quantum noise.
As mentioned earlier, another avenue is the description of enhancement and inhibition of SBS through density of states engineering [52]. In the spontaneous regime, this is well handled by a Fermi Golden Rule calculation of phonon generation rates [53], as we will demonstrate in a subsequent work.
and ∂π n (r, t) ∂t where in the second last line we have used (5).

S.II. OPERATOR M nk (r) IS HERMITIAN
Here we show that the operator M nk (r) of (15) is Hermitian. Consider an integral over an appropriate volume and assume that fields are either periodic over the volume or vanish at the surface of the volume. Then for vector functions C(r) and D(r), integrating by parts twice gives (D n (r)) * M nk (r)C k (r) dr = − (D n (r)) * Now using (5) we put c njkl (r) = c klnj (r) and switching the dummy indices j ↔ l, we can write this as and so the differential operator M nk (r) is Hermitian.

S.III. PROPERTIES OF MODE FUNCTIONS AND PARTNER FUNCTIONS
Here we establish some useful properties of the mode functions F Λ (r) introduced in (17) to (20), that are required to reduce the acoustic Hamiltonian to canonical harmonic oscillator form in S.IV.
Note that since M nk (r) is real, if F Λ (r) is an eigenfunction then F * Λ (r) is also an eigenfunction with the same eigenvalue ω Λ . This may happen simply because F Λ (r) is purely real. In fact, as we show in S.XI, it is always possible to choose the set of eigenfunctions {F Λ (r)} such that each of them is purely real. But it is often more convenient to work with complex eigenfunctions (traveling waves rather than standing waves, for example). Section S.XI establishes that if we include complex eigenfunctions in the set {F Λ (r)}, the set can be chosen so that each eigenfunction F Λ (r) is either real or, if not, there is another eigenfunction FΛ(r) in the set such that FΛ(r) = F * Λ (r). Typically the naturally chosen set of eigenfunctions will make this so; for example, if F Λ (r) is a traveling wave to the right, then FΛ(r) is a traveling wave to the left. We refer to FΛ(r) as the "partner" of F Λ (r). That is, each complex eigenfunction in the set has a partner that is also in the set. If there are purely real eigenfunctions in the set {F Λ (r)}, we take them to be their own partners. Then the set of eigenfunctions {F Λ (r)} is equivalent to the set of {FΛ(r)} of partner eigenfunctions, and FΛ(r) = F * Λ (r) for each Λ. Since the set of partners is equivalent to the original set, then from Eq. (20) of the main paper we can also write Similarly, from Eq. (18) of the main paper we see that we havẽ

S.IV. ACOUSTIC MODE EXPANSION
Here we show that through use of the partner functions introduced in S.III, the acoustic field operators and Hamiltonian can be expanded in terms of the mode functions as expressed in Eqs. (22) and (23). Recall that the {F Λ (r)} are the eigenfunctions of (15) and that the Ũ Λ (r) and Π Λ (r) are defined as in Eq. (18). Both sets of functions are proportional to the {F Λ (r)}, so we can take each to constitute a complete set of states. We can then expandũ Λ and C Λ are operators and the factors /2Ω Λ are added for later convenience.
In the Heisenberg picture, C Λ and C Λ are time-dependent, and therefore so areũ(r) andπ(r). However, sincẽ u(r) andπ(r) are Hermitian the C so that (S.4) may also be written as a sum over partner modes, where we have used ΩΛ = Ω Λ . Then from the Hermiticity of the canonical fields we see that we require Λ , which may be satisfied by setting For real eigenfunctions F Λ (r) we haveΛ = Λ and this just says that C Λ is (proportional to) a coordinate operator, and C Λ is (proportional to) a momentum operator (as we will see). For partners, b Λ and bΛ are independent operators (or independent amplitudes in the classical case). Using (S.5) in (S.4) we havẽ where in the third line we use the fact that the sum is over the whole set of partner functions. Similarlỹ in accordance with Eq. (22). Postulating the commutation relations where we have used (20) and (S.2). Since [ũ m (r),π m (r)] = [u m (r), π m (r)] we recover the starting commutation relations (2). It is possible but more complicated to show that demanding the result (S.10) one can find that the b Λ and b † Λ must satisfy (S.9). Now we look at the Hamiltonian in terms of the b Λ and b † Λ . From (16) we have From the above we haveũ In the last two terms, orthogonality gives Λ = Λ. In the first two, we replace the sum over Λ by a sum overΛ and use the fact that F Λ (r) = F * Λ (r); then orthogonality demands that Λ =Λ. Recalling that ΩΛ = Ω Λ we then have Then since we have Using the same strategy as above this gives Combining (S.11,S.11) we have

S.V. THE GROUP VELOCITY
Here we work out the group velocity of the acoustic modes in terms of the modal field providing an explicit expression for Eq. (35a), a result which is needed in the next section. We take the continuous limit of Eq. (26), writing It is helpful to re-express M nk in terms of an operator L nk q operating on the transverse spatial variables only.
Applying M nk to the mode (S. 11) gives where the last line defines the action of the operator L nk q on f αq (x, y). It follows from the Hermiticity of M nk that L nk q is Hermitian with respect to integration over the transverse plane. We can then write so that the f n αq are eigenfunctions of L nk q and may be taken as orthogonal. Taking the inner product with (f n αq ) * and using the orthogonality of the f αq , we have Differentiating with respect to q gives 2Ω αq dΩ αq dq = d dq dxdy (f n αq (x, y)) * L nk q f k αq (x, y). (S. 16) Since L nk q is Hermitian, we may invoke the Hellmann-Feynman theorem to simplify the right hand side: Then the group velocity Swapping the dummy indices n ↔ k in the second term and using (5) gives where the final line follows from (28).

S.VI. ACOUSTIC POWER FLOW
Even in the presence of coupling the displacement to the electromagnetic fields, or other forces, we expect the first of (8a) still to hold, Since in general the power density at a point in the medium in a directionn is classically given by the power in the waveguide in theẑ direction, integrated over the xy plane, is where the second line follows from the symmetry properties of the stiffness tensor. We form the operator corresponding to the classical P cl (z) by a usual procedure. Since P cl (z) involves the product of the classical fields π i (r) and ∂u l (r)/∂r m , we obtain the operator P (z) by using the symmetrized version of the operators corresponding to π i (r) and ∂u l (r)/∂r m : where we put Using (31), we see that K ilm (r) has the form dqdq Ω αq Ω α q κ ilm αα (q, q ), (S. 27) where κ ilm αα (q, q ) =κ ilm αα (q, q ) +κ ilm αα (q, q ). (S.28) The first term contains parts rapidly-varying in space and time: andκ ilm αα (q, q ) contains the slowly-varying terms, We write where P rv (z) contains the contributions fromκ ilm αα (q, q ) and P sv (z) those fromκ ilm αα (q, q ). Our interest is in the latter. Since the sums and integrals in (S.27) are over all α, α , q, and q , we may switch the dummy indices in the second term on the right-hand-side of (S.29): Moving to normal-ordering with where L ilm α α (q , q; x, y) = π i α q (x, y)e iq z * ∂ ∂r m u l αq (x, y)e iqz + π i αq (x, y)e iqz ∂ ∂r m u l α q (x, y)e iq z * , (S.32) The term involving T ilm α (q) in Eq. (S.31) represents vacuum zero-point contributions and should give no net contribution to P sv (z), which is a directed quantity. Indeed, using Eqs. (28) and the property f k αq (x, y) = (f k α(−q) (x, y)) * , which follows from the Hermiticity of M nk (see S.III), it can be shown that its contribution to Eq. (S.27) vanishes.
The remaining contribution to P sv can be written where the pairwise term p A α a (q , q) = − 2 Ω αq Ω α q dxdy c izkl ρ L ikl α α (q , q; x, y).
Evaluating the derivatives in Eq. (S.32), Using (28) we have Consequently, In the second term we may exchange i and k because the other elements of the stiffness tensor are both the same to obtain which by Eq. (S.19) is simply the group velocity of the acoustic mode. The desired result (40) then follows from (39).

S.VII. ELECTROMAGNETIC POWER FLOW
Here we justify the relations (53) to (55) in the main paper for the optical power transport in terms of the optical field envelope operators.
The operator for the power carried by the field is given by the Poynting vector which we write in the symmetrized form The temporally slowly-varying part of this expression is Since the k integral is over all wavenumbers including all partner modes, it is easy to show using (S.41) that the second term in this expression, associated with vacuum contributions, vanishes, as we would expect for a signed quantity. For the remaining non-vacuum contribution, since the sums and integrals are over all values we may swap the indices γ, γ and k, k in the second term in square brackets to give where we have introduced the quantity Finally, if the different modes γ have very different center wavenumbers k j , then only the γ = γ terms will contribute significantly to (S.45) and we may approximate with p EM γγ (k, k) the power carried by the normalized mode functions γ at center wavenumber k.

S.1. Interpretation as the photon number density operator
To convert the result in (S.47) to a simple expression involving the photon envelope operators we require the group velocity in terms of the fields.
Noting that the basis functions B γk (r) = b γk (x, y)e ikz are eigenmodes of the vector Helmholtz equation (43) where the k-dependent operator O k operates on a vector function f as and where ∇ t = [∂ x , ∂ y , 0]. It can be shown that O k is Hermitian such that From Ampere's law, we also have that We now take the inner product with b * γk in (S.48) and differentiate both sides with respect to k: where we used the normalization dxdy b * γk · b γk /µ 0 = 1 which follows from (46) and Maxwell's equations. By the Hermiticity of O k , we can invoke the Hellmann-Feynman theorem to write the left hand side as Since the inverse dielectric tensor is symmetric even under strain, we have p ijlm (x, y) = p jilm (x, y), and swapping the dummy indices k, k in the second term gives × dxdy d i γk (x, y) * d j γ k (x, y) p ijlm (x, y)s lm αq (x, y) − δ ij ∂β ref (x, y) ∂r l u l αq (x, y) .

S.IX. SMOOTHING THE SURFACE MATRIX ELEMENT
In the (x, y) plane we can in general identify a number of curves C that indicate where β ref (x, y) will change discontinuously from one value to another. These may be straight or curved lines. We only contemplate discontinuous changes in β ref (x, y), adding up the neighborhoods of all such curves identifies the regions where β ref (x, y) is assumed to vary in the xy plane. We write R = (x, y), and parameterize such a curve by R c (s) = (x c (s), y c (s)), and for a given curve let s range from 0 to 1. As s increases along the curve we have x = x c (s) + ζ (x ·n(s)) , y = y c (s) + ζ (ŷ ·n(s)) .
For fixed s, β ref changes as ζ passes from < 0 to > 0; that is, it is only a function of ζ. We assume now that the change in β ref (x, y) occurs only in such a small region about R = R c (s) (we will eventually take that change to be a Dirac delta function there) that the mapping from (ζ, s) to (x, y) is one-to-one. Then we can write where in the second line we understand x = x(s, ζ) and y = y(s, ζ); that is, we have switched integration variables from x and y to s and ζ. where d ⊥ γ k (x, y) =n(s) · d γ k (x, y), d γ k (x, y) =ẑẑ · d γ k (x, y) +û (s)û (s) · d γ k (x, y) = d γ k (x, y) −n (s)n (s) · d γ k (x, y),