Ab-Initio Calculations of Nonlinear Susceptibility and Multi-Phonon Mixing Processes in a 2DEG-Piezoelectric Heterostructure

Solid-state elastic-wave phonons are a promising platform for a wide range of quantum information applications. An outstanding challenge and enabling capability in harnessing phonons for quantum information processing is achieving strong nonlinear interactions between them. To this end, we propose a general architecture using piezoelectric-semiconductor heterostructures consisting of a piezoelectric acoustic material hosting phonon modes in direct proximity to a two-dimensional electron gas (2DEG). Each phonon in the piezoelectric material carries an electric field, which extends into the 2DEG. The fields induce polarization of 2DEG electrons, which in turn interact with other piezoelectric phononic electric fields. The net result is coupling between the various phonon modes. We derive, from first principles, the nonlinear phononic susceptibility of the system. We show that many nonlinear processes are strongly favored at high electron mobility, motivating the use of the 2DEG to mediate the nonlinearities. We derive in detail the first, second, and third-order susceptibilities and calculate them for the case of a lithium niobate surface acoustic wave interacting with a GaAs-AlGaAs heterostructure 2DEG. We show that, for this system, the strong third-order nonlinearity could enable single-phonon Kerr shift in an acoustic cavity that exceeds realistic cavity linewidths, potentially leading to a new class of acoustic qubit. We further show that the strong second-order nonlinearity could be used to produce a high-gain, traveling-wave parametric amplifier to amplify--and ultimately detect--the outputs of the acoustic cavity qubits. Assuming favorable losses in such a system, these capabilities, combined with the ability to efficiently transduce phonons from microwave electromagnetic fields in transmission lines, thus hold promise for creating all-acoustic quantum information processors.

The ability of phonons to act as a universal quantum bus or even as a qubit in universal linear acoustics are undoubtedly useful, but these abilities would be far more useful if one could perform quantum logic operations on the phonons, themselves.This, of course, requires singlephonon nonlinearities that are large compared to system losses.No material currently possesses such strong nonlinearities and low losses; thus, it is imperative to explore other means for generating these properties.Recently, we have demonstrated that phononic nonlinearities can be increased by orders of magnitude by mediating them with electronic nonlinearities in a piezoelectric-semiconductor heterostructure [47].Classical theories of such nonlinearities [48][49][50] show that the nonlinearity depends directly on the mobility.Thus, one must wonder what would be possible in materials with ultra-high mobilities, such as two-dimensional electron gasses (2DEGs), which can attain mobilities multiple orders of magnitude higher than bulk semiconductors [51,52].
To this end, we study a heterostructure consisting of a 2DEG bonded to a strongly piezoelectric elastic FIG.1: Diagram of a system consisting of an acoustic Kerr cavity (left) coupled to a degenerate traveling-wave parametric amplifier (right).The coupling in the acoustic cavity between the 2DEG (red) and the piezoelectric material (blue) gives rise to a Kerr nonlinearity and in turn a qubit.The qubit is then amplified by the 3-wave-mixing in the parametric amplifier due to the coupling between the 2DEG there (red) and the piezoelectric material (blue), before being transduced into a microwave field at the system output.Note that the thin waveguide carries the pump phonons.slab waveguide substrate (e.g.thin-film LiNbO 3 on Si, sapphire, or SiC or thin-film Al 1-x Sc x N on SiC).Only the critical layers of the 2DEG remain on the surface of the piezoelectric substrate, which can be achieved with modern heteroepitaxy and chemical etching (e.g. in a GaAs/AlGaAs heteroepitaxial stack).Similar piezoelectric-semiconductor phononic heterostructure systems have been realized in the classical regime utilizing bulk semiconductors to achieve phonon amplification [14,[53][54][55], phononic switches [56], phononic circulators [53], phononic mixing [47], and acoustoelectrically Brillouin optomechanical processes [57].A 2DEGpiezoelectric heterostructure has itself been experimentally realized in the context of measuring a large acoustoelectric effect [58], while more recently, a quantumdot-piezoelectric heterostructure has been realized for optomechanical applications [59].In all these cases, including the one at-hand, the acoustic field in the piezoelectric material takes the form of a slab waveguide mode or surface acoustic wave.As in the case of the classical systems, an evanescent electric field associated with the elastic wave phonons extends into the 2DEG, where it induces the polarization of 2DEG electrons.The polarized electrons, in turn, interact with other electric fields stem-ming from other elastic wave phonon excitations in the system, yielding mixing of multiple phonons mediated by the 2DEG electrons.
Here, we theoretically analyze the electron-phonon interactions in the above system and perform ab-initio quantum calculations of the resulting acoustic nonlinearities and multi-phonon mixing processes that result from them.This allows for the exact calculation of phononic mixing processes such as second harmonic generation, parametric amplification, Kerr nonlinearity, etc.We find, for example, that these kinds of systems could produce high-gain and low-loss quantum-limited degenerate parametric amplification via three-wave mixing with only nanowatts of acoustic pump power; we also find that the effective Kerr nonlinearity in the system could be sufficient to produce an effective two-level system in a wavelength-scale cavity with sufficiently low losses, given a high-Q acoustic cavity, which has recently been achieved [60].Thus, our results allow one to imagine a new class of quantum information processing tools for phonons based on electronically mediated nonlinearity in heterostructures and calculate the exact details of these processes.
Such a system is illustrated in Fig. 1.Here, a 2DEG in a nearly single-mode cavity-the qubit cavity-creates a qubit at frequency (wavevector) ω 0 (k 0 ) via four-wave phononic Kerr nonlinearity.The qubit cavity is coupled to a tunable readout cavity that allows adiabatic transfer of the qubit state into a purely phononic state of the readout cavity.The readout cavity is strongly coupled to an output phononic bus waveguide that leads to a degenerate, traveling wave parametric three-wave mixing phononic amplifier.A pair of pump waveguides inject and extract the pump field at frequency (wavevector) wω 0 (2k 0 ) via acoustic evanescent coupling.The amplified phononic field at ω 0 can then be transduced into a microwave field at the system output, to be measured by room temperature electronics after a series of progressively higher temperature amplifiers.Note that elastic waves with highly restricted transverse dimensions can be efficiently excited (low microwave insertion loss and low elastic wave loss) by focusing transducers [61][62][63], phononic negative refraction in phononic crystals [64,65], and phononic crystal gradient-index (GRIN) lenses [66].Essentially, with the inclusion of 2DEGs playing the role of spatially selective sources of single-phonon nonlinearities, one can imagine circuit QED for phonons that mimic the functionalities of those of superconducting circuit QED.The manuscript is organized as follows.In Sec.II, we derive the electric field per phonon experienced by the 2DEG electrons based on the piezoelectric material's properties, as well as the 2DEG electric field screening.In Sec.III, we derive the susceptibility of the 2DEG to the electric fields associated with the elastic wave phonons.In Sec.IV, we present numerical and analytical results for the 2DEG susceptibility in certain relevant limits.In Sec.V, we combine the susceptibility values from Sec. IV with the electric-field-per-phonon values from Sec. II to calculate the resulting phononic susceptibility of the system and multi-phonon mixing dynamics.As a check on validity, in Sec.VI we derive the classical susceptibilities up to second order based on our quantum susceptibility model and compare to existing classical results in the literature.

II. DERIVING THE ELECTRIC FIELD PER PHONON
We start with the critical task of deriving the electric field per phonon from the elasticity, dielectric coupling, and piezoelectric coupling tensors given a shear-wave surface acoustic field.In general, the following relationships must hold for any piezoelectric material [67]: where c, ϵ, and e are the elasticity, dielectric, and piezoelectric coupling tensors, respectively, while T , S, E, and D represent the stress, strain, electric field, and electric displacement field, respectively.Our first goal is to determine the electric field for a given strain.In the absence of an external electric field, the electric displacement field D goes to 0 regardless of how the coordinate system is defined, yielding the following relationship between the electric field and the strain field based on Eq. ( 2): The electric field per phonon can thus be calculated from the strain amplitude per phonon.
We thus turn to solving the strain amplitude for a given mode as a function of phonon number in that mode.In general, for a surface acoustic wave propagating along the plane of the surface, the displacement amplitude will vary along the axis perpendicular to the surface.Defining the x-axis as the direction of propagation and the z-axis as the axis perpendicular to the surface, we can express the displacement field as a superposition of modes as follows: where the functions f n,q are orthogonal for a given wavevector q: dzf n,q (z)f m,q (z) = L z δ m,n . ( The strain components corresponding to the displacement u(r) are calculated as follows: where g ij,n,q is defined as follows: where θ n,q,m is defined as the angle between the displacement vector u n,q and the m-axis.
The goal is to determine the global amplitude u n,q as a function of the number of phonons N in the mode (n, q).To that end, we relate the rate of change of the displacement u to the kinetic energy T and use the fact that the total mode energy N ℏω is twice the kinetic energy for a simple harmonic oscillator: where V = L x L y L z is the mode volume in the piezoelectric material.Note that if the displacement consists of a superposition of modes, then there will be no direct coupling between the modes, due to two factors.First, if the modes are of different wavelengths, then the integral of the product of the amplitudes over the x-axis (i.e., the propagation axis) will go to zero.Second, even for two modes of the same wavelength, the fact that the z-varying functions f are orthogonal to one another ensures that the integral of the product of the amplitudes over the z-axis will go to zero.Finally, having solved for the strain amplitude per phonon, we calculate the electric field per phonon at the surface.In general, the electric field amplitude corresponding to a plane wave of frequency ω and wavevector q propagating in the m-direction can be expanded in a manner analogous to the displacement and strain amplitudes: We solve for the z-dependent electric field amplitude E n,q as a function of the phonon numbers in the respective modes using Eq. ( 3) as follows: We can assume that a surface acoustic wave exponentially decays into the bulk with a decay length similar to the wavelength, yielding f n,q (z) ∝ e −qz .Then, given a diagonal dielectric tensor (i.e., ϵ ij = ϵ ii δ i,j ), the electric field amplitude at the surface (which we define as z = 0) reduces to the following: where C is the effective ratio between the piezoelectric coupling coefficient and the square-root of the elasticity (i.e., in a linear piezoelectric material, C = e/ √ κ, where e and κ are the piezoelectric coupling and elasticity coefficients, respectively), ϵ is the effective dielectric constant, and A ph is a measurement of the square-root of the propagating strain field power density (i.e., power per unit area) in units of ℏω.Each of these parameters is defined in the following manner: For the electron-phonon interaction, we use the dielectric tensor for the 2DEG as ϵ, instead of the dielectric tensor for the piezoelectric material.In Sec.V, we will use A ph as a measurement of the overall field amplitude for each mode.

III. GENERAL DERIVATION OF HIGHER-ORDER SUSCEPTIBILITIES
Having determined the electric field amplitude per phonon, we now derive the general N th -order susceptibility of the 2DEG electrons as a function of the input electric field amplitudes.By definition, the N th -order susceptibility corresponds to the average electron dipole moment (normalized to material volume and amplitude of each input field) induced by N field modes simultaneously interacting with the electrons.We specifically consider the interaction at a quantum level: i.e., how N phonons simultaneously polarize a single electron via the electric field amplitudes associated with the respective phonons.First, we calculate the amplitude for a single electron-phonon interaction, which is proportional to the strength of an electron's dipole moment if it successfully interacts with N phonons.Then, we derive the probability (as a function of the field amplitudes) that the electron successfully interacts with N phonons, which yields the N th -order susceptibility.

A. Deriving the electron-phonon interaction amplitude
Here, we derive the electron-phonon interaction Hamiltonian from first principles, ultimately yielding the effective electron dipole moment as a function of the acoustic field wavelength.As previously discussed, the phonon acts on the electron through the piezoelectric field associated with the phonon.Generically defining the effective electron dipole moment along the axis of field polarization as d eff and the electric field amplitude as E, we apply the standard dipole interaction Hamiltonian: To solve for the dipole matrix elements, we thus need to calculate the Hamiltonian matrix elements.The reason for using this indirect method is that given the longitudinal nature of the electric field associated with the acoustic field, the dipole moment cannot be separately promoted to operator form in the basis of 2DEG electronic states, since the electric field also spatially varies along the axis of propagation.We therefore start with the fundamental piezo-2DEG interaction Hamiltonian, given an electric field polarized in the x-direction: where P is the 2DEG polarization.For an acoustic field (and thus an associated piezoelectric field) propagating in the x-direction with a wavevector q, the electric field E x takes the following form in terms of the acoustic field latter operators b ( †) : where E zpf is the zero-point piezoelectric field amplitude.Note that the spatial dependence of the field is entirely isolated to the phases e ±iqx .Next, we express the polarization P x in the basis of 2DEG electronic states.Here, we assume that the spacing between consecutive impurities is significantly smaller than the y-direction span of the 2DEG, yet much larger than the z-direction span (which represents the axis perpendicular to the piezo-2DEG interface).Consequently, the Fermi level cuts through several subbands in the k ydirection, while cutting through only the lowest subband in the k z -direction (with the next k z -subband being located a vast distance from the lowest subband in phase space).We can therefore ignore k z and index all states as (k x , k y ).The polarization P x can thus be expanded in terms of the carrier transitions between any initially occupied ground state (k i , k y ) and an initially unoccupied excited state (k i + ∆k, k y ) [68]: where c kx,ky is the fermionic annihilation (creation) operator for the electronic state at (k x , k y ), ω ki+∆k,ki = ω ki+∆k − ω ki is the change in carrier energy (this is independent of k y ), ϕ(z) is the wavefunction along the z axis, L y is the span of the 2DEG along the y-axis, and ξ is defined as follows: where L is the waveguide length.It is worth noting that ω ki+∆k,ki can be solved as follows: Substituting the above expressions into the Hamiltonian in Eq. ( 16), we express the Hamiltonian in electronphonon interaction form as follows: As desired, the spatially-varying phases in the field and the polarization combine to enforce momentum conservation.The effective dipole matrix element thus becomes a function of the phonon wavevector amplitude q: d eff,kx+∆k,kx (q) = − q e q δ q,∆k , for any initial wavevector k x .This implies that the matrix element amplitude is inversely proportional to the phonon momentum, and hence directly proportional to the acoustic wavelength.It is worth analyzing the reason why the dipole moment is proportional to the acoustic wavelength but independent of the material length.Intuitively, this results from the alternating nature of the piezoelectric field (as shown in Fig. 2), which limits the separation of charges to the order of a half-wavelength, regardless of the material dimensions.As we increase the number of phonons (corresponding to an increase in field intensity), the dipole moment per polarized electron will stay constant, but the number of polarized electrons will increase, as will be shown in the next subsection.

B. Deriving the higher-order interaction probabilities
We now seek to determine the higher-order electronphonon interaction probabilities, which yields the polarization field as a function of the acoustic (and thus piezo-electric) field strength, by deriving the generic N th -order susceptibility in terms of the general electronic decay rate γ = q e /(mµ), and from this deriving the conversion efficiency for specific nonlinear processes.As in the case of an atom interacting with an optical field [69], the energy conservation requirement for the intermediate states is relaxed due to the momentary nature of these intermediate states (per the energy-time uncertainty principle).However, unlike the atom-optical interaction, the much slower propagation speed of the acoustic signal ensures that each phonon features a significant momentum relative to its energy.The need for momentum conservation thus collapses the range of possible intermediate states, as the results in this section will show.
We start by considering the self-Hamiltonian of each excited 2DEG state.We incorporate the decay rate γ as an anti-Hermitian effective Hamiltonian term [70], yielding the following unperturbed Hamiltonian for the 2DEG states: where ω lg = ω l − ω g represents the energy of the excited state |l⟩ (with the energy of the ground state |g⟩ set as the baseline), and γ l = γ for all excited states |l⟩ in our system (though we express the decay rate independently for each state for the sake of generality).Turning to the piezoelectric field, we specifically consider a set of waves of varying frequencies polarized along the x-axis and traveling in the +x-direction.The fields interact FIG.2: Conceptual diagram of how the alternating electric field E associated with the acoustic field induces dipole formation in the 2DEG.Note that the dipole size is limited to the order of a half-wavelength.
with the electron dipole moment d via the following timedependent interaction Hamiltonian: where q p = ω p /v s , and the spatial variation of the field is incorporated into the effective dipole moment d eff (which is defined as in Eq. ( 22)).Note that E(−ω p ) = E * (ω p ), ensuring the realness of the composite field at all points.We can also express a generic field of frequency ω b > 0 in terms of the corresponding ladder operators a This will play an important role later in this section when we convert the semiclassical fields to quantum form for the purpose of deriving the frequency conversion rate from the Heisenberg equation of motion.
The interaction Hamiltonian enables us to derive the expectation value of the effective dipole moment, ⟨d eff ⟩, as a series expansion in the field amplitudes E. Reducing the expectation value of the interaction Hamiltonian using the rotating-wave approximation (such that we only consider the time-independent ⟨H int ⟩ terms), we express this expectation value as follows: where the N th -order susceptibility χ (n) is defined as a function proportional to the N th -order expansion of ⟨d eff ⟩, in the case where a unique set of frequencies ω p1 , ..., ω pn add up to −ω p : FIG. 3: Depiction of the ladder of intermediate states for the N -step transition described by χ (N ) (ω p1 , ..., ω p N ), with the horizontal (vertical) axis representing the wavevector q (frequency ω).Note that the set of generic modes (r 1 , ..., r N ) can be composed of any permutation of the set of input modes (p 1 , ..., p N ), with the requirement that ω r1 + ... + ω r N = ω p1 + ... + ω p N .Also note that each intermediate state is broadened in energy but fixed in momentum.
where |ψ (n) (t)⟩ is the n th -order expansion of the composite 2DEG wavefunction |Ψ(t)⟩.Note that the variation of ⟨d eff ⟩ with the field amplitudes in Eq. ( 25) is entirely due to the variation of the density of polarized electrons rather than the dipole strength corresponding to a given polarized electron (which is fixed by the acoustic wavelength, as discussed in the previous subsection).Intuitively, χ (N ) (ω p1 , ..., ω pn ) corresponds to the probability that an electron becomes polarized upon interacting with N fields of wavevectors q p1 , ..., q p N respectively (normalized to the amplitude of each field and to the material volume), multiplied by the strength of the dipole moment itself.
Appendix A shows how χ (N ) is calculated for a set of N propagating acoustic waves with wavevectors q p1 , ..., q p N .As for the 2DEG states, the well-defined momentum for each state ensures strict momentum conservation for each transition.However, due to the spectral broadening encapsulated in the decay rate γ, the energy conservation requirement is somewhat relaxed.This dynamic is depicted in Fig. 3. Due to the availability of permutations within the set of modes (p 1 , ..., p N ), we have left the mode indices in the generic form (r 1 , ..., r N ), where ω r1 + ... + ω r N = ω p1 + ... + ω p N , and hence q r1 + ... + q r N = q p1 + ... + q p N .The transition to (from) |l m ⟩ is mediated by the mode r m for m ≤ N − j (m > N − j).The transition from |l N −j ⟩ to |l N −j+1 ⟩ is unique in that it is not mediated by any individual mode in the set (p 1 , ..., p N ), but rather by the composite mode −p 1 −...−p N .Quantitatively, this mode will become relevant when considering the interaction Hamiltonian term associated with χ (N ) (ω p1 , ..., ω p N ).
The most critical takeaway from Fig. 3 is that while the intermediate 2DEG states are broadened in frequency ω (corresponding to a wide range of possible real states), each intermediate state is fixed in wavevector q, since the effective dipole moment d eff conserves momentum.As a result, for a given initial state, each intermediate state reduces to a single real state (rather than a superposition of multiple states).The summation over intermediate states in Eq. (A8) thus reduces to a sum over initial states (which we index with the 2D wavevector (k x , k y )): where we have used the fact that d eff,k f ,ki (q) = (q e /|q|)δ q,k f −ki , and the symbol P preceding the summation over p 1 , ..., p N denotes a permutation of the mode indices.The expression thus includes N !(N +1) = (N +1)! terms, due to the N possible values of j and the N !per-mutations of the field mode indices.

IV. SUSCEPTIBILITY RESULTS
Here, we specifically consider the susceptibility results up to third order.In calculating the nonlinearity produced by the third-order Kerr susceptibility, it is also important to compare its value to the spectral broadening.For a high-Q acoustic cavity (which has recently been achieved [60]), the spectral broadening is dominated by the phonon absorption by the 2DEG electrons, for which we need to determine the imaginary part of χ (1) .In performing these calculations, we make the assumption that the average spacing between consecutive charge carriers is much smaller than the acoustic wavelength, corresponding to the limit q ≪ k F (where q and k F are the acoustic and Fermi wavevectors, respectively).
The general procedure for calculating the susceptibilities from phase-space integrals is derived in Appendix B. We start by considering the third-order Kerr susceptibility.To this end, the degenerate four-wave-mixing suscep-tibility χ (3) (ω, ω, −ω 1 ) (corresponding to the processes where 2 phonons of frequency ω each are absorbed, and 2 phonons of frequency ω 1 and ω 2 = 2ω −ω 1 are emitted) reduces to the following form: where f kx (ω, ω, −ω 1 ) is approximately the following: where we have defined q = ω/v s and q 1 = ω 1 /v s .The first two terms are resonant, corresponding to the absorption of 2 phonons of frequency ω each, followed by the emission of a phonon of frequency ω 1 and another of frequency ω 2 = 2ω − ω 1 .On the other hand, the latter two terms are counter-resonant, corresponding to the emission of ω 1 and ω 2 phonons followed by absorption of two ω phonons.The process described by each term is depicted in Fig. 4. It is worth noting that if we ignored the initial occupation of the states, then there would be 12 separate processes contributing to χ (3) (ω, ω, −ω 1 ).However, assuming that q << k F , and given the fact that all states within the Fermi circle |k| < k F are occupied, we can make the approximation that the only valid processes are either two absorptions followed by two emissions (resonant) or two emissions followed by two absorptions (counter-resonant).The third-order Kerr susceptibility equals the degenerate four-wave-mixing suscepti-bility in the limit ω 1 → ω.
Next, we examine the second-order susceptibility, focusing first on the case of sum-frequency generation encapsulated by χ (2) (ω 1 , ω 2 ), where two phonons of frequency ω 1 and ω 2 are absorbed and a phonon of frequency 2ω = ω 1 + ω 2 is emitted: where f (2) ) is approximately the following: FIG. 4: Phase-space depiction of processes contributing to χ (3) (ω, ω, −ω 1 ), with (a), (b), (c), and (d) corresponding respectively to the first, second, third, and fourth terms of Eq. ( 29).The shaded (unshaded) regions represent initially occupied (unoccupied) electronic states.Note that the crescent region represents the range of allowed initial states, while electrons in the rest of the Fermi circle are prohibited from interacting with the given phonons.
where q 1 = ω 1 /v s and q 2 = ω 2 /v s .The first two terms are resonant and the latter two terms are counterresonant, with the processes depicted in Fig. 5.As with the case of χ (3) , we apply the high-carrier-density limit q ≪ k F , causing the requirement for unoccupied intermediate states to approximately reduce the 6 processes that would feed into a generic second-order process to 4 processes.The susceptibility governing the second-harmonic generation process equals the above susceptibility for the general sum-frequency generation in the limit ω 1 = ω 2 .
In the context of second-order susceptibility, it is also worth examining the case of parametric amplification  31).The shading scheme is identical to that in Fig. 4.
encapsulatd by χ (2) (2ω, −ω 1 ), where a pump phonon of frequency 2ω is absorbed, causing one phonon each of the signal frequency ω 1 and the idler frequency ω 2 to be emitted.Since the frequencies involved are the same as the sum-frequency generation, the coefficients for χ (2) (2ω, −ω 1 ) are the same as those for χ (2) (ω 1 , ω 2 ) in Eq. ( 30).However, f kx (2ω, −ω 1 ), defined as: Note the sign flip (relative to Eq. ( 31)) in the second imaginary term of the first and fourth processes.This has important implications for the relationship between the parametric-amplification susceptibility and the mobility in the case of degenerate parametric amplification (ω 1 = ω 2 ), as we will show in Sec.IV D. Finally, we calculate the linear susceptibility χ (1) (ω) corresponding to the absorption and emission of a phonon of frequency ω: where f (1) kx (ω) is defined as follows: where q = ω/v s .The first term is resonant, while the second is counter-resonant, as shown in Fig. 6.As previously discussed, these susceptibility expressions govern the corresponding processes in the high-carrier-density limit.The susceptibility values can be numerically evaluated from these expressions.However, it is useful to derive the analytical forms for the purpose of optimization.In the coming sections, we will do so in two limits: low-mobility and high-mobility.

A. High-Mobility Regime
Here, we seek to derive a general procedure for analytically determining the susceptibilities in the high-mobility (low-decay) regime.We start by considering an electron in its ground state inside the Fermi circle.The probability that this electron interacts with an phonon of wavevector q = ω/v s propagating in the x-direction is governed by the detuning ∆ω between the phonon's frequency and the transition frequency corresponding to the electron jumping from k x → k x + q.Consequently, ∆ω is a function of the electron's initial x-direction wavevector k x but is invariant in the initial y-direction wavevector k y : where ω 0 = v s q denotes the phonon's angular frequency, and k x,0 = mv s /ℏ − q/2 represents the initial wavevector corresponding to a fully resonant electron-phonon interaction.
We therefore focus on the near-resonance region, i.e., the band of electrons around the initial wavevector k x,0 ranging in detuning from ∆ω ≈ −γ/2 to ∆ω ≈ γ/2.In the high-carrier-density/high-mobility regime, the k xspan of this band is much less than the diameter 2k F of the Fermi circle.We solve for the k x -span of the nearresonance band from Eq. ( 35) as follows: Conveniently, the k x -span of the near-resonance band is proportional to the electronic decay rate γ (and hence inversely proportional to the mobility µ).Intuitively, this is due to the fact that the electron-phonon detuning varies linearly with k x .Applying the high-mobility requirement k x,span ≪ 2k F , we derive the following condition for γ: which corresponds to the following condition for the mobility µ, upon applying the definition µ = q e /(mγ): Since the susceptibility corresponding to each initial electron is invariant in k y , the total susceptibility corresponding to all electrons with a particular initial k x is proportional to the k y -span of valid initial electronic states with the given k x -value.In the limit of high carrier density (i.e., k x,0 ≪ k F ), the k y -span given in Eq. (B2) for k x values near k x,0 can be approximately linearized in k x as follows: We can express any given k x near k x,0 as k x = k x,0 +∆k x , where ∆k x is the deviation of the initial wavevector from that corresponding to a resonant electron-phonon interaction.Then, per Eq. ( 35), the detuning ∆ω varies linearly with the wavevector deviation ∆k x as ∆ω = ℏq(∆k x )/m.Substituting this relationship as well as the definition of k x,0 into the above expression, we find that to first-order, the k y -span varies linearly with the detuning ∆ω: where we have also substituted the phonon frequency ω 0 = v s q.Using the fact that dk x = md(∆ω)/(ℏq) = mv s d(∆ω)/(ℏω 0 ), a generic susceptibility term f (∆ω) can be integrated over phase space in the following man- FIG.6: Phase-space depiction of processes contributing to χ (1) (ω), with (a) and (b) corresponding respectively to the first and second terms of Eq. ( 34).The shading scheme is identical to that in Fig. 4. ner: Note that we only need to consider the resonant terms for f (∆ω), since the counter-resonant terms are negligible relative to the resonant terms in the high-mobility limit.
Intuitively, these results are explained by Fig. 7. Two observations are worth making here.First, we note that the k y -span of valid initial states in the near-resonance band varies inversely with the Fermi wavevector k F (as shown in Fig. 7(b)), ensuring that the overall susceptibility also varies inversely with k F as well.It might seem counterintuitive that a larger Fermi circle would yield a smaller span of valid initial states.This can be visually explained based on Fig. 7(a) by the fact that for the given value of k x , as the Fermi wavevector is increased, the total range of k y included inside the Fermi circle also increases, but the part of that range inside the forbidden (purple/circle-minus-crescent) zone grows even faster due to a reduced circle curvature.Therefore, the range of k y included inside the valid (green/crescent) region actually decreases as the Fermi wavevector increases.Since k F varies as √ n (where n is the carrier density), this implies that any susceptibility term varies as n −1/2 (with the exceptions of the real parts of χ (1) and χ (2) , discussed later in this section).
The second observation regards the scaling of the nearresonance band area with the decay rate γ.For convenience, we define a new decay parameter γ ′ = γ/2.For each electron in the near-resonance band, the average N th -order susceptibility is proportional to 1/γ ′N .If a single-electron susceptibility term is even in the detuning ∆ω, then the total susceptibility corresponding to that term is simply determined by multiplying the average single-electron susceptibility by the area of the nearresonance band (which is proportional to γ ′ , as shown in Fig. 7(b)).Therefore, the overall susceptibility for the process represented by an even term is proportional to 1/γ ′N −1 .However, if a single-electron susceptibility term is odd in ∆ω, then the contributions from the electrons featuring positive and negative electron-phonon detuning values largely cancel each other out.In this case, the net contribution comes from the difference in the number of electrons between the negative-detuning and positivedetuning cases.As Fig. 7(b) shows, this difference is proportional to γ ′2 rather than γ ′ when integrated over the band.As a result, the overall susceptibility for the process represented by an odd term is proportional to 1/γ ′N −2 .
It is worth addressing the special cases of Re[χ (1) ] and Re[χ (2) ].In the case of Re[χ (1) ], the integral does not converge.Therefore, the contributing electron-phonon interactions are concentrated far off-resonance rather than near-resonance.This is evidenced by the fact that the resonant and counter-resonant terms for Re[χ (1) ] are always approximately equal even in the high-mobility limit.Consequently, Re[χ (1) ] is invariant in the mobility.It is also important to note that since the contributions to Re[χ (1) ] are far off-resonance, the initial electronic states are primarily concentrated around k x = k F .Since ∆ω ∝ k x for k x ≫ q, and Re[χ (1) ] ∝ 1/(∆ω) for far-offresonance electrons, the susceptibility for each electron varies inversely with k F .On the other hand, the size of the region of valid initial states varies directly with k F (i.e., the radius of the Fermi circle).As a result, these two variations cancel out, and Re[χ (1) ] is invariant in the Fermi wavevector and thus independent of carrier density.
The case of Re[χ (2) ] is complicated by the fact that the near-resonance contributions cancel out when integrated.On the other hand, it also converges to zero for high detuning values, unlike Re[χ (1) ].Consequently, in the high-mobility limit, the electronic contributions to Re[χ (2) ] are neither dominantly near-resonance or far-offresonance, but rather in between the two limits.We will thus tackle the problem of Re[χ (2) ] in the high-mobility limit solely through numerical means.

B. Low-Mobility Regime
Here, we focus on deriving a procedure for analytically determining the susceptibilities in the low-mobility (high-decay) limit.Conceptually, this corresponds to the regime in which the k x -span of the near-resonance band is much greater than the diameter 2k F of the Fermi circle, such that all electrons are included in the band and the electron-phonon detuning for any electron is negligible compared to the electronic decay rate γ.Following a derivation analogous to the high-mobility regime, the following conditions on γ and the mobility µ hold for the low-mobility regime: In this regime, the far-off-resonance electrons will dominate the susceptibility terms.We start by considering how the real-valued susceptibility terms are reduced.In general, the real part of f (N ) can be expressed as a power series in the electron-phonon detuning ∆ω as follows: where each a n is a positive integer and each c N,m represents a coefficient that derives from expanding the numerator.In the low-mobility limit, we can make the approximation |∆ω| ≪ γ ′ .As such, we can reduce each summation to the leading term in γ ′ .This yields a real part of f (N ) that varies linearly with the detuning for odd N , while being invariant in the detuning for even N (assuming that this result does not cancel out): Consequently, in the low-mobility limit, the real part of χ (N ) varies with the mobility µ as µ N (µ N +1 ) for even (odd) N , assuming that the leading term does not cancel out.To determine the overall χ (N ) values, we devise a scheme to integrate the single-electron values over the entire phase space of valid initial states rather than just the near-resonance region.To this end, we approximate the k y -span by applying the approximation The available phase-space area, which we label A phase , is solved by simply integrating the k y -span over k x as follows: Since the single-electron processes feeding into Re[χ (N ) ] for even N are invariant in the detuning, the overall susceptibility can simply be calculated by multiplying Re[f (N ) ] by the phase-space area A phase .To determine the overall Re[χ (N ) ] for odd values of N , however, we need to integrate the detuning ∆ω over this phase-space area.Since the electrons in the overall range of valid initial states are concentrated near k x = k F , we can make the approximation k x ≫ q and thus ∆ω ≈ ω kx+q,kx ≈ ℏqk x /m.We thus integrate ∆ω over the valid initial phase-space as follows: Note that the result as k 2 F .This is due to two factors.The first is that the available phase space area (corresponding to the number of available electrons) scales as k F .The second is that the probability that an electron interacts with the phonon fields also scales as k F in the given regime, since ⟨∆ω⟩ ∝ k F .As such, Re[χ (N ) ] scales linearly with the carrier density n for odd N .On the other hand, for even N , Re[χ (N ) ] scales as n 1/2 , since the electron-phonon interaction probability in this case is invariant in the detuning ∆ω.
It is also worth examining the integral of (∆ω) n more generally over the valid initial phase space.This is especially useful for the cases where the leading term cancels out, leaving behind secondary terms that vary as (∆ω) n where n > 1: This result varies as k n+1 F , since the number of available electrons scales linearly with k F , while the interaction probability for each electron scales as k n F .Finally, we comment on the imaginary parts of χ (N ) .Separating these into individual terms corresponding to specific processes, the imaginary parts of the individual susceptibility terms f (N ) take the following form: Using the same argument as with the real parts, the leading term reduces to the following in the low-mobility limit: The leading term for Im[χ (N ) ] therefore scales with mobility as µ N +1 (µ N ) for even (odd) N .However, it is critical to note that in the low-mobility regime, the resonant and counter-resonant terms largely cancel out, corresponding to the fact that stimulated emission largely balances out absorption.As a result, the overall Im[χ (N ) ] for a generic N must be calculated by making a higherorder expansion of the terms to account for the slight difference between the absorption and emission probabilities.

C. First-Order Susceptibility
We now apply our findings for the low-mobility and high-mobility limits to derive the first-order susceptibility in these limits.In parallel, we numerically determine the susceptibility as a function of mobility.We start with the real part of χ (1) , which governs the dielectric constant for the 2DEG.As discussed in the previous sections, Re[χ (1) ] follows a special rule in the high-mobility limit, while following the generic behavior in the low-mobility limit.
We thus focus first on the low-mobility limit.Substituting Eq. ( 34) into Eq.( 44) and simplifying per Eq. ( 45), we find that Re[f (1) kx ] approximately reduces to the following: Note that the real part of χ (1) thus scales with mobility as µ 2 in the low-mobility limit.Integrating ∆ω over the phase-space area of the initial electronic states as in Eq. ( 48) and multiplying by the constants shown in Eq. ( 33), we find the following analytical expression for the real part of χ (1) in the low-mobility limit: Re[χ (1) where we substituted the relationship γ ′ = q e /(2mµ) in the last line.Note that this susceptibility is independent of the phonon frequency ω 0 in the limit q ≪ k F .This is because the total number of electrons (encapsulated by the phase space area) and the real part of the inverse detuning per electron each varies linearly with ω 0 , thus cancelling out the inverse dependence of the dipole interaction probability with ω 2 0 .Next, we examine the real part of χ (1) in the highmobility limit.As discussed in Sec.IV A, Re[χ (1) ] is dominated by the far-off-resonance interactions, unlike the general case for the high-mobility limit.Using the approximation |∆ω| ≫ γ ′ , we find that Re[f (1) kx ] approximately reduces to: As with the generic low-mobility case, we use the fact that the valid initial states are concentrated near k x = k F , leading to the approximation k x ≫ q and thus ∆ω ≈ ω kx+q,kx ≈ ℏqk x /m.Integrating (∆ω) −1 over the valid initial phase-space area yields the following result: Multiplying this by 2 (per Eq. ( 54)) and by constants shown in Eq. ( 33), we find the following analytical expression for the real part of χ (1) in the high-mobility limit: Re[χ (1) Note that this is independent of both mobility and carrier density.The invariance in mobility can be explained by the fact that the dielectric screening is dominated by faroff-resonance electron-phonon interactions, causing the detuning to dominate over the decay rate, as explained earlier.The invariance in carrier density n is explained by the fact that the number of electrons available to interact with phonons varies as √ n (since the corresponding phase-space area varies as k F ), while the probability that any given electron actually interacts with a phonon varies as 1/ √ n (since in the far-off-resonance limit, this probability scales inversely with the average detuning, which in turn varies as k F and thus as √ n).We now turn to the imaginary part of χ (1) , which governs the net phonon absorption rate and thus the spectral broadening for the phonon energy levels.As discussed in Sec.IV B, the imaginary parts follow a special rule in the low-mobility limit due to cancellation between absorptive and emissive processes.We thus start by analyzing the high-mobility limit.Here, we can approximately ignore the second (counter-resonant) term in f (1) kx in Eq. ( 34), yielding: Using the near-resonance approximation and the corresponding procedure in Eq. ( 41), we integrate this over the near-resonance band (noting that the integrand is even in ∆ω) and multiply by the coefficients in Eq. ( 33), leading to the following result for the imaginary part of χ (1) in the high-mobility limit: As such, the imaginary part of χ (1) is invariant in γ ′ and hence independent of mobility in the high-mobility limit.This is because the number of electrons in the near-resonance region scales linearly with γ ′ , while the interaction probability for each electron scales inversely with γ ′ , causing the overall interaction rate to be constant in γ ′ .Finally, we solve for the imaginary part of χ (1) in the low-mobility limit.If we take the resonant and counter-resonant terms (corresponding to phonon absorption and emission, respectively) in Eq. ( 34) individually, then the magnitude of the imaginary part of each of those terms should approximately reduce to 1/γ ′ , since γ ′ ≫ |∆ω| ≫ ω 0 in the low-mobility limit.However, this dynamic is fundamentally altered by the fact that the two terms cancel each other out to first-order.This requires us to make a higher-order approximation, altering the mobility-dependence of the net susceptibility: Intuitively, this result can be conceptualized as follows: the probability of the 2DEG electrons absorbing a phonon from the field (represented by the resonant term) is almost (though not fully) cancelled out by the probability of 2DEG electrons emitting a phonon into the field.The former (deamplifying) process narrowly edges out the latter (amplifying) process, yielding a positive imaginary part of the electron-phonon interaction probability.As with the real part of f kx , the imaginary part is approximately linear in the electron-phonon detuning ∆ω.We therefore integrate over the phase-space area of initial electronic states as in Eq. ( 48) and mutiply by the coefficients in Eq. ( 33): Although the overall susceptibilities corresponding to the absorptive and emissive processes individually scale linearly with the mobility µ, the net process scales as µ 3 .As with the real part of χ (1) in the low-mobility regime, this result also varies linearly with carrier density n ∝ k 2 F , due to a combination of the overall number of available electrons varying and the single-electron interaction probability each varying as n 1/2 .Figure 8 depicts the numerically-calculated results for the real and imaginary parts of χ (1) (given a phonon angular frequency of ω 0 = 2π×10 9 s −1 , speed of sound v s = 4 × 10 3 m/s, a carrier density n = 2 × 10 15 m −2 , a carrier effective mass m = 0.067m 0 relative to the bare electron mass m 0 , and a 2DEG thickness t 2DEG = 2 × 10 −8 m), along with the analytical results in the high-mobility and low-mobility limits.As desired, the numerical and analytical results match in the low-mobility and highmobility regimes, with the real (imaginary) part varying as µ 2 (µ 3 ) in the low-mobility regime, and with both parts plateauing in the high-mobility regime.Note that for all mobility values, Re[χ (1) ] vastly exceeds Im[χ (1) ].The dielectric constant can thus be approximated as fully real rather than complex, which features important implications when calculating the wave-mixing dynamics.

D. Second-Order Susceptibility
We now derive the second-order susceptibility, using numerical means to calculate the susceptibility for all mobility values and analytical means to specifically derive the low-mobility limit.We start with sum-frequency generation, governed by χ (2) (ω 1 , ω 2 ).In the low-mobility limit, we can reduce Re[f (2) kx ] in Eq. ( 31) to the following using the approximation |∆ω| ≫ ω kx+q1+q2,kx , ω 1 , ω 2 : The sum-frequency-generation susceptibility thus varies with mobility as µ 2 in this regime, following the general low-mobility rule laid out in Eq. (45).We can analytically solve for the overall susceptibility in the lowmobility regime by multiplying by the available phasespace area in Eq. ( 47) and the coefficients in Eq. (30).It is important to consider, however, that for two of the terms in Eq. ( 31), the width of the phase-space area is q 1 = ω 1 /v s (see Figs.  other two terms, the width is q 2 = ω 2 /v s (see Figs. 5(b) and (c)).Consequently, for the available phase-space area expression in Eq. ( 47), we apply the replacement ω 0 → (ω 1 + ω 2 )/2, yielding the following overall secondorder susceptibility in the low-mobility limit: Note the linear variation in k F , implying that the susceptibility scales with the carrier density as n 1/2 .This is because the number of available electrons (corresponding to the available phase-space area) varies linearly with k F , while the interaction probability per electron is independent of k F . Figure 9 depicts the numerically-calculated results for the amplitude of the real part of χ (2) (ω 1 , ω 2 ) in the degenerate limit ω 1 = ω 0 = ω 2 representing the second-harmonic-generation process (we label this quantity χ (2) SHG (ω 0 )), along with the analytical result in the low-mobility limit, given the same parameters as in the χ (1) calculation.As desired, the numerical and analytical results match in the low-mobility regime, with the susceptibility varying as µ 2 .It is interesting to note that the susceptibility plateaus in the high-mobility regime.As discussed in Sec.IV A, the near-resonance electronic contributions to the real part of χ (2) cancel out in the high-mobility regime, leaving behind electronic contributions that are somewhat off-resonance (but not far-offresonance, since Re[χ (2) ] converges to zero for large detuning values).Along with the very low value of γ, this FIG.9: Results for the amplitude of the real part of χ (2) (ω 0 , ω 0 ) as a function of the carrier mobility µ, including numerical (solid, green) and low-mobility analytical (dotted, red) calculations.We use a phonon angular frequency of ω 0 = 2π × 10 9 s −1 , speed of sound v s = 4 × 10 3 m/s, a carrier density n = 2 × 10 15 m −2 , a carrier effective mass m = 0.067m 0 , and a 2DEG thickness t 2DEG = 2 × 10 −8 m. causes the effect of the decay rate γ on the electron interaction probabilities to be negligible, leading Re[χ (2) ] to be constant in mobility in the high-mobility limit.
Next, we turn to parametric amplification, focusing specifically on the degenerate case (ω kx ] from Eq. ( 32) in the low-mobility limit, we note that the amplitude of each of the terms can be approximated as 1/γ ′2 to first order, as with secondharmonic generation.However, the key difference here is that the leading terms cancel out for both the resonant case and the counter-resonant case, due to destructive interference of probability waves in the degenerate traveling-wave parametric amplifier.Consequently, we shift to the secondary term, which equals 2(∆ω) 2 /γ ′4 (where ∆ω = ω kx+q0,kx − ω 0 ).As a result, the parametric amplification susceptibility varies with mobility as µ 4 instead of µ 2 : To analytically solve for the low-mobility susceptibility, we need to integrate (∆ω) 2 over the available phasespace area.Since the electron-phonon transitions in the low-mobility limit are concentrated far-off-resonance, we approximate ∆ω ≈ ℏqk x /m, enabling us to apply the generic higher-order integral shown in Eq. ( 49) as follows: Multiplying this f (2) PA and the coefficients in Eq. ( 30), we find the following degenerate parametric amplification susceptibility in the low-mobility limit: Note that this is independent of the frequency ω 0 .This is because the number of available electrons varies linearly with ω 0 , while the interaction probability for each electron varies inversely with ω 0 , and the two variations cancel out.Furthermore, since the interaction probability per electron now varies as k 2 F , the susceptibility scales as k 3  F , implying a variation with carrier density as n 3/2 .Figure 10 depicts the numerically-calculated results for the amplitude of the real part of χ (2) (2ω 0 , −ω 0 ) (we label this quantity χ (2) PA (ω 0 )).As desired, the numerical and analytical results match in the low-mobility regime, with the susceptibility varying as µ 4 .The susceptibility levels out in the high-mobility regime for the same reason as in the second-harmonic-generation case.

E. Third-Order Kerr Susceptibility
We conclude our susceptibility calculations by deriving the third-order Kerr susceptibility χ (3) (ω 0 , ω 0 , −ω 0 ) governing the Kerr shift and the degenerate four-wavemixing process.In the low-mobility limit, we simplify Re[f (3) kx ] in Eq. ( 29) to the following, since |∆ω| ≫ ω kx+2q,kx , ω: where ∆ω = ω kx+q,kx −ω 0 .Following the general rule for the low-mobility regime as outlined in Eq. ( 45), the Kerr susceptibility varies with mobility as µ 4 in this regime.We analytically solve for the overall susceptibility in the low-mobility limit by integrating over the phase-space area of initial electronic states using the procedure in Eq. ( 48) and multiplying by the coefficients in Eq. ( 28): The quadratic variation in k F implies a linear scaling with the carrier density n.This is because the number of available electrons and the interaction probability per electron each scales as n 1/2 .
We now turn to the high-mobility limit.Here, we have to account for the fact that the transition frequency ω kx+2q,kx+q for the electron upon absorbing a second phonon of wavevector q exceeds that upon absorbing the first phonon ω kx+q,kx by a constant (k x -independent) offset frequency, which we label ω ′ : In terms of this offset, Re[f kx ] in Eq. ( 29) takes the following form, considering only the resonant (first and sec-ond) terms: where the top and bottom indices in the ± and ∓ operations represent the first and second terms of Eq. ( 29), respectively.We integrate this over the near-resonance band (using the procedure in Eq. ( 41)) and multiply by the coefficients in Eq. ( 28), leading to the following result for the overall Kerr susceptibility in the high-mobility limit: Re[χ (3) where β is defined as follows: Note that the Kerr susceptibility varies linearly with µ in the high-mobility limit.However, the slope of this variation is 3 times greater in the ultra-high-mobility limit than in the medium-high-mobility limit.This is because the linear variation originates from different sources for the two cases.In the medium-high-mobility case, the electronic spectrum effectively becomes harmonic, thus making Re[f (3) (∆ω)] odd in ∆ω and yielding Re[f (3) (∆ω) ∝ 1/γ ′3 .As such, integrating the product of this function and ∆ω over the near-resonance band results in Re[χ (3) (ω 0 , ω 0 , −ω 0 )] ∝ 1/γ ′ ∝ µ.On the other hand, in the ultra-high-mobility limit, the size of the electronic anharmonicity in phase space exceeds the width of the near-resonance band, making Re[f (3) (∆ω)] even in ∆ω and yielding Re[f (3) (∆ω) ∝ 1/γ ′2 .Integrating this over the near-resonance band also results in Re[χ (3) (ω 0 , ω 0 , −ω 0 )] ∝ 1/γ ′ ∝ µ, but with a different slope than the medium-high-mobility case.
Figure 11 depicts the numerically-calculated results for the amplitude of the real part of χ (3) (ω 0 , ω 0 , −ω 0 ), which we label χ (3) Kerr (ω 0 ).We use the same parameters as in the χ (1) and χ (2) calculations.As desired, the numerical and analytical results match in the low-mobility and high-mobility regimes, varying as µ 4 in the lowmobility regime and as µ 1 in the high-mobility regime.The aforementioned change in the proportionality coefficient with µ between the medium-high-mobility and ultra-high-mobility regimes is apparent in the intercept shift in the log-log plot around 10 3 m 2 /(V • s).FIG.11: Results for the amplitude of the real part of χ (3) (ω 0 , ω 0 , −ω 0 ) as a function of the carrier mobility µ, including numerical (solid, green), low-mobility analytical (dotted, red), and high-mobility analytical (dashed, blue) calculations.We use a phonon angular frequency of ω 0 = 2π × 10 9 s −1 , speed of sound v s = 4 × 10 3 m/s, a carrier density n = 2 × 10 15 m −2 , a carrier effective mass m = 0.067m 0 , and a 2DEG thickness t 2DEG = 2 × 10 −8 m.

V. CALCULATION OF MIXING DYNAMICS
The N th -order susceptibility manifests itself in the (N + 1)-wave-mixing process.For a quantitative understanding, we can substitute the expression for ⟨d eff (q p )⟩ in Eq. ( 25) into the interaction Hamiltonian in Eq. ( 24), yielding an expansion in product of field amplitudes: Intuitively, this corresponds to the fact that the different field modes are coupled to one another via interactions with the 2DEG electrons.At the elemental level, the field modes p 1 through p N simultaneously polarize the electrons, which in turn interact with the field mode p through field-dipole coupling.It is apparent that the real part of the susceptibility corresponds to the conversion efficiency between modes (or to an energy shift in the acoustic states, for the case of mode self-coupling), while the imaginary part yields a broadening of the acoustic states caused by loss due to phonon absorption by the electrons.
The Hamiltonian corresponding to the interaction of a phonon from each of the modes p 1 , ..., p N and a phonon of frequency ω = ω p1 + ... + ω p N can be expressed in operator form as follows: where, as evidenced by Eq. ( 73), the (N + 1)-phonon interaction rate g is given by the following: where E zpf represents the zero-point electric field experienced by the 2DEG electrons.Intuitively, since χ (N ) represents the total induced electron polarization per unit amplitude of the N polarizing fields p 1 , ..., p N , multiplying it by V 2DEG E zpf (ω p1 )...E zpf (ω p N ) yields the product of the probability that an electron becomes polarized by absorbing N phonons and the resulting dipole moment of the polarized electron.This polarized electron then interacts with the (N + 1)-st field (represented here by the frequency ω) through the dipole interaction.The net process is electron-mediated coupling of N + 1 phonons from the respective fields, with the interaction rate given by Eq. ( 75).
In general, we can determine the time-evolution of the output electric field ω in terms of the amplitudes of the input electric fields using the Heisenberg equation of motion as follows: The zero-point electric field E zpf is determined by substituting A ph = v s /(2V ω ) (corresponding to a vacuum fluctuation spanning the mode volume) into Eq.( 11), yielding: where V ω is the overall mode volume for the acoustic wave.From this, we can determine the spatial evolution of the output field amplitude A ph (ω) in terms of the input amplitudes A ph (ω 1 ), ...A ph (ω N ): Note that the conversion efficiency is proportional to the ratio between the 2DEG volume V 2DEG and the overall mode volume V ω (which simply reduces to the ratio between the 2DEG thickness t 2DEG and the mode depth L z,ω if the 2DEG covers the piezoelectric material's surface).This is because the fraction of the total mode energy stored in the 2DEG field is proportional to the ratio of the two volumes.
It is also worth noting that the conversion efficiency varies inversely with ϵ N +1 , where ϵ is the effective 2DEG dielectric constant, which is in turn proportional to Re[χ (1) ].Naively, we would assume that ϵ = ϵ 0 Re[χ (1) ].However, there are 2 critical caveats.The first is that due to the thin-film nature of the 2DEG (i.e., the fact that the 2DEG thickness is much smaller than the field wavelength), the dielectric screening is strongly suppressed.As explained in Appendix C, this attenuates the dielectric constant by a factor a, which is a function of the ratio between the 2DEG thickness t 2DEG and the field's half-wavelength l: The second caveat is that when Re[χ (1) ] is low enough that the corresponding dielectric constant is below the dielectric constant of the piezoelectric material ϵ p , the latter provides the effective dielectric constant for the 2DEG, since the fact that the 2DEG is much thinner than the piezoelectric material ensures that the fringe induced field from the piezoelectric material fully penetrates the 2DEG.As such, the dielectric constant can be approximated as follows: The 2DEG dielectric constant is plotted in Fig. 12 given ϵ p = 43.6ϵ0 , with all the other parameters the same as in the susceptibility plots.Given the dielectric constant and the susceptibilities, we solve for the rates of 3 processes: linear absorption (corresponding to Im[χ (1) ]), as well as second-harmonic generation and parametric amplification (corresponding to Re[χ (2) ]).We also derive the figure-of-merit for the Kerr nonlinearity by calculating the ratio between the phonon Kerr shift (correspond- ing to Re[χ (3) ]) and the phonon spectral broadening (corresponding to Im[χ (1) ], given a high-Q acoustic cavity, which has recently been achieved [60]).

A. Linear Absorption
We start with the case of linear absorption.Based on Eq. ( 78), the spatial evolution of a single field propagating through the heterostructure and experiencing linear absorption takes the following form: where α is the amplitude attenuation rate, given as: Figure 13 depicts the absorption rate given C 2 = 6.1ϵ 0 and effective mode depth L z = 2 × 10 −6 m (i.e., the half-wavelength of the ω 0 field) along with the same parameters used in the susceptibility calculations, assuming that the 2DEG covers the piezoelectric material such that V 2DEG /V = t 2DEG /L z .At very low mobility, the imaginary part of χ (1) varies as µ 3 , while the dielectric constant is flat in mobility, causing the absorption rate to vary as µ 3 .For moderately low mobility, though, the dielectric-constant-squared varies as µ 4 , causing the absorption rate to decrease with mobility as µ −1 .The linear absorption coefficient reaches a constant value of about 25 m −1 at the high-mobility limit (since both the dielectric constant and the imaginary part of χ (1) are constant at high mobility), corresponding to an effective phonon coherence length of 40 mm.

B. Second-Harmonic Generation
Next, we solve for the second-harmonic-generation coupling strength.From Eq. ( 78), the output field spatially evolves as follows: where 2ω 0 = ω 1 + ω 2 , and G is the coupling strength, given as: where V 2ω0 is the mode volume corresponding to the 2ω 0 field.Here, the half-wavelength yields an effective mode depth of L z,2ω0 = 1 × 10 −6 m.For degenerate secondharmonic generation (i.e, ω 1 = ω 0 = ω 2 ), the coupling strength becomes the following: Figure 14 depicts the second-harmonic coupling strength using the same parameters as in the linear absorption case (except for the mode depth, as previously mentioned): At very low mobility, χ (2) varies as µ 2 , while the dielectric constant is flat in mobility, causing |G| to vary as µ 2 .Once the dielectric constant starts rising with mobility, however, the strong inverse dependence of |G| with the dielectric constant (specifically, the ϵ −3 variation) causes |G| to sharply drop with increasing mobility despite the fact that χ (2) increases with mobility in this regime.However, at the high-mobility limit, both χ (2)  and the dielectric constant reach constant values, causing |G| to level out at 7.6 × 10 −11 s 1/2 .

C. Parametric Amplification
Here, we solve for the parametric amplification gain per unit length.We specifically design a phase-sensitive, degenerate, traveling-wave parametric amplifier due to its ability to achieve the minimum possible additive quantum noise for an amplification process [71][72][73][74][75].It is worth noting that resonant, degenerate parametric amplifiers could alternatively be designed to minimize the pump power required for a given gain, though a phasepreserving amplifier can never achieve the same additive quantum noise as a phase-sensitive (traveling-wave) amplifier [71,76].
A diagram of a traveling-wave parametric amplifier is shown in Fig. 15.Given an input signal field with frequency ω 1 and a pump field with frequency 2ω 0 = ω 1 +ω 2 , the output field ω 2 evolves in a manner analogous to the output field in Sec.V B: where G p,1 is the parametric conversion efficiency per unit length per pump phonon current density amplitude (defined as the square-root of the number of pump phonons per unit time per unit cross-sectional area), defined as: Conversely, if we consider ω 2 as the input signal field, then the time-evolution of ω 1 is modeled by switching ω 1 and ω 2 in Eq. ( 86), yielding: We now specifically consider the case of a degenerate traveling-wave parametric amplifier (i.e., ω 1 = ω 0 = ω 2 ).
If we set the phase of the pump field A ph (2ω 0 ) such that it is in phase with −iA ph (ω 0 )/A * ph (ω 0 ), then the timeevolution of the signal/idler field amplitude is solved by summing Eqs. ( 86) and (88): As such, the amplitude |A ph (ω 0 )| grows exponentially at a rate 2|G p,1 | A ph (2ω 0 ) .The signal/idler field power and thus evolves at twice the rate of the amplitude: The gain per unit length per pump phonon current density amplitude is thus expressed in dB form as: An alternative metric for parametric amplification is the gain per unit length per pump power amplitude (defined as the square-root of the pump power).This is solved from G p,2 in the following manner: Note that G p features an inverse dependence on the pump mode's cross-sectional area L y L z,2ω0 in addition to the standard V 2DEG /V ω0 variation for G p,1 .This is because for a given pump power, the pump field amplitude increases with decreasing cross-sectional area.Since the gain is proportional to the pump field amplitude rather than the power, this causes the gain to increase for lower cross-sectional area values.Figure 16 depicts the parametric amplification gain per unit length per pump power amplitude, given the same parameters as in the second-harmonic-generation case, except that we use a signal/idler frequency of 5 GHz (corresponding to ω 0 = π × 10 10 s −1 ) and a mode width of 15 wavelengths (i.e., L y = 12 µm).The mode depths corresponding to the signal and pump fields, respectively, are the half-wavelengths L z,ω0 = 400 nm and L z,2ω0 = 200 nm.At very low mobility, χ (2) (2ω 0 , −ω 0 ) varies as µ 4 , while the dielectric constant is invariant in mobility, causing G p to vary as µ 4 , peaking at 43 dB/(µm √ µW), corresponding to a 20-dB gain in 4.7 µm (about 2.4% of the Rayleigh length) given a 10-nW pump power input.As the mobility is increased from this point, though, the strong inverse dependence of the gain with the dielectric constant causes G p to decline with increasing mobility FIG.15: Diagram of the degenerate traveling-wave parametric amplifier system.Note that the S-curve on the left carries the pump input, while the curve on the right drains the pump phonons remaining after the 3-wave-mixing process.The region in between represents the core amplifier region, with the 2DEG (red) on top of the piezoelectric material (blue).for the intermediate-mobility case, as with the secondharmonic case.Eventually, at the high-mobility limit, both χ (2) (2ω 0 , −ω 0 ) and the dielectric constant become constant, causing G p to level out at 11 dB/(µm √ µW), corresponding to a 20-dB gain in 19 µm (about 9.5% of the Rayleigh length) given a 10-nW pump power input.

D. Kerr Nonlinearity
Finally, we solve for the Kerr nonlinearity figure-ofmerit (F.O.M.), which we define as the ratio between the Kerr-induced anharmonicity and the spectral broadening.This will determine whether the electron-phonon interaction induces enough anharmonicity in the phonon spectrum such that the system can be used as an artificial atom.A diagram of the acoustic cavity giving rise to the Kerr nonlinearity is shown in Fig. 17.From Eq. ( 75), the anharmonicity induced by the Kerr shift in a phonon ladder of fundamental angular frequency ω 0 is calculated by doubling the Kerr coupling coefficient, which we calculate from the third-order Kerr susceptibility χ On the other hand, given a high-Q acoustic cavity (which has recently been achieved [60]), the spectral broadening is given by the linear absorption per unit time, which corresponds to Im[χ (1) ].We convert the result for the linear amplitude attenuation coefficient per unit propagation distance α from Sec. V A into per-unit-time absorption by doubling it (to convert from amplitude attenuation to energy absorption) and multiplying by the speed of sound v s , yielding: where γ ph is the phononic spectral broadening.The Kerrnonlinearity figure-of-merit (F.O.M.) is therefore solved as follows: As this expression shows, the Kerr nonlinearity is maximized by minimizing the acoustic field's mode volume.This is because the Kerr shift varies with the electric field intensity per phonon, which varies inversely with mode volume.
Figure 18 depicts the Kerr nonlinearity figure-of-merit using the same parameters as in the linear absorption and second-harmonic generation cases, except that the angular frequency ω 0 is changed to π × 10 10 s −1 to align with the 5-GHz transmon frequency.We minimize the mode volume by setting the mode length equal to the half-wavelength of 400 nm along each dimension.At high mobilities, the Kerr figure-of-merit scales linearly with µ.

This is because χ
(3) Kerr also scales linearly with µ in this regime, while both the real and imaginary parts of χ (1)  are constant in mobility.At very high mobilies (above 10 4 m 2 /(V•s)), the figure-of-merit exceeds unity, indicating that the phononic system can be used as an artificial atom.

VI. COMPARISON OF QUANTUM AND CLASSICAL RESULTS IN FLUID LIMIT
Although the previous sections treated the carrier states as coherent fermionic states, it is interesting to analyze the system in the incoherent (fluid) carrier limit for the purpose of comparing the quantum results for the linear absorption (corresponding to Im[χ (1) ]) and the sum-frequency generation (corresponding to Re[χ (2) ]) processes to classical results from previous studies.Here, we operate in the low-mobility limit and effectively treat the electrons as bosonic.Due to the latter, all electrons are capable of interacting with phonons, and due to the former, the probability of such an interaction is approximately uniform for all of the electrons.Furthermore, we only consider the resonant interaction processes, in line with previous theoretical studies of various types of interactions [77][78][79][80][81].

A. Linear Absorption in Fluid Limit
We start by considering the linear absorption.In the low-mobility regime, the resonant part of Im[f (1) kx ] from Eq. ( 34) reduces to: Im[f (1)  res ] = which is uniform for all electrons, as expected.Then, if we make the assumption that all electrons in the Fermi circle are capable of absorbing phonons, we solve for the imaginary part of χ (1) simply by multiplying this expression by the Fermi circle area and the coefficients in Eq. ( 33), yielding: where we substituted k 2 F = 2πn in the second line (where n is the area carrier density) and γ ′ = q e /(2mµ) in the third line.Now, we apply another consequence of the fluid treatment of electrons: that all ground-state electrons feature zero momentum (k = 0).Then, a momentum kick from absorbing a phonon of wavevector q would elevate a ground-state electron's energy to ℏq 2 /(2m).If we assume that the electron-phonon interaction is resonant, then ℏq 2 /(2m) must equal v s q (i.e., ℏω 0 = 2mv 2 s ).Substituting this into the above expression for the imaginary part of χ (1) yields the following result: Im[χ (1) ](ω 0 ) ≈ q e n bulk µ ϵ where we defined a bulk volumetric electron density n bulk = n/t 2DEG , along with the frequency parameter ω c = n bulk q e µ/ϵ.It is worth noting that the imaginary part of χ (1) varies linearly with both carrier density and mobility via ω c .The former is due to the fact that a higher carrier density leads to a larger number of electrons available for interacting with phonons, while the latter is because for a given electron, the electron-phonon interaction probability scales inversely with the spectral broadening in the limit in which the detuning is much smaller than the broadening.Next, we determine the phonon absorption per unit length α from the imaginary part of χ (1) by multiplying by ϵ 0 V 2DEG |E(ω 0 )| 2 /ℏ to determine the absorption per unit time and dividing by the propagation speed v s .Recall that the zero-point electric field intensity is a function of the effective dielectric constant ϵ eff , the elasticity κ, and the frequency ω 0 : where ϵ is the dielectric screening (i.e., the real part of χ (1) ) and K 2 = e 2 /(ϵκ) represents the square of the electromechanical coupling.This accords with Eq. ( 77), except for the fact that ϵ eff incorporates both the real and imaginary parts of χ (1) for the sake of generality: (100) The absorption per unit length α thus becomes the following: where in the first line, we inserted v s in the denominator to convert from per-unit-time to per-unit-length absorption, while in the second line, we used the fact that V piezo = V 2DEG for a bulk piezoelectric semiconductor.Note that this result matches that derived previously in the classical limit [82].

B. Sum-Frequency Generation in Fluid Limit
We now turn to the question of resolving the real part of χ (2) in for the purpose of deriving the dynamics of sum-frequency generation.In the low-mobility regime, the resonant part of Re[f (2) kx ] from Eq. ( 31) reduces to the following: Re[f (2)  res ] which is also uniform for all electrons.Then, assuming that all electrons are capable of undergoing second-order interaction with phonons, we solve for the real part of χ (2)  by multiplying this expression by the Fermi circle area and the coefficients in Eq. ( 30), yielding the following result: where ω 3 = ω 1 + ω 2 , and we substituted k 2 F = 2πn and γ ′2 = q 2 e /(4m 2 µ 2 ) in the second and third lines, respectively.
Applying the ground-state approximation described in the linear absorption derivation, along with the singlephonon resonance condition that ℏq 2 /(2m) must equal v s q (i.e., ℏ 2 ω 2 = 4m 2 v 4 s , where ω 1 ≈ ω ≈ ω 2 ), we find that the real part of χ (2) reduces to the following: where n bulk and ω c are defined as in the linear absorption calculation, and Next, we use the real part of χ (2) to determine how the ω 3 field evolves in time given input fields ω 1 and ω 2 .To that end, the rate of change of the field E 3 due to the sum-frequency-generation process can be calculated by promoting the fields to operators, such that E n = E zpf,n a n , where a n is the lowering operator for the mode n and E zpf,n is the zero-point electric field amplitude.Recall that the electric field amplitude is directly proportional to the strain field amplitude: where Γ n is defined as follows for the given mode n: Also, recall that the zero-point electric field intensity varies linearly with the phonon frequency and inversely with the mode volume as follows: (107) The time-evolution of E 3 due to sum-frequency generation is thus solved by applying the 3-wave-mixing Hamiltonian: Finally, we solve for E (2) = −iω 3 E (2) 3 ): Note that this result matches that calculated in the classical limit [48] if K 2 → 1 and |Γ 3 | 2 → 1.

VII. CONCLUSION
Using time-dependent perturbation theory, we have theoretically derived the nonlinear acoustic susceptibilities up to an arbitrary order for a technologically feasible heterostructure consisting of a 2DEG stacked on top of a piezoelectric material.We have also used this susceptibility to derive the equations of motion for various nonlinear quantum phononic processes.We presented the first, second, and third order susceptibilites in detail, both analytically and numerically, and their variations with important parameters such as carrier concentration and mobility.Our results demonstrate that both the second-order and third-order susceptibilities are maximized at high mobility, demonstrating the advantage that a 2DEG creates by virtue of its ultra-high mobility.We note that in the case of classical acoustoelectric mediation of phonon-mixing nonlinearities [47], the application of bias fields causes profound modification of these nonlinear processes.These modifications range from amplification or deamplification in space of the waves involved in the mixing processes, to relaxation of phase-matching conditions, to direct modification of the phononic susceptibility.Though the treatment of the interaction of these quantum phononic mixing processes with bias fields is beyond the scope of this paper, we expect that in this regime, the effect of bias fields will be similarly profound, and we plan to analyze it in future work.
It is particularly noteworthy that the Kerr nonlinearity continues to increase with mobility, even at high mobility values.If we are able to achieve a sufficiently highmobility in a 2DEG and fabricate a device such that a heterostructure island mostly fills the inside of an acoustic cavity, then the phenomena such as phononic cavity blockades (analogous to photon blockades [83,84]) and other types of single-phonon level quantum interactions or even logic gates are achievable.The relevant metric to enter this regime is the Kerr shift per phonon inside the cavity relative to the cavity linewidth, given a high-Q acoustic cavity (which has recently been achieved [60]).If we can produce such a system with a high enough mobility and small enough cavity linewidth such that only the transition from the ground state to a single phonon excitation is effectively resonant with the cavity, then this system could provide a new type of artificial atom, analogous to a superconducting circuit qubit and at the same microwave frequencies but with a length scale of microns instead of centimeters.Furthermore, the large second-order nonlinearity of our heterostructure system at high mobilities could be used to create a parametric phonon amplifier.This combination of artificial atom cavities and traveling-wave parametric amplifiers, together with bus waveguides, make a complete analog to the primary hardware components of superconducting circuit quantum computers.Together with a method for microwave-to-phonon transduction (wellestablished interdigital electrode transducers) and for tuning the coupling rate between cavities [27,85], one has the pieces necessary to perform state preparation, single and multi-qubit gates via tunable cavity-cavity coupling, and room temperature readout after parametric amplification.Thus, it could be possible to make microwave frequency phononic quantum processors based on these heterostructures that are highly analogous to superconducting circuit quantum computers and operating at the same frequency and temperature but with many orders of magnitude higher qubit density.We note these qubits would require no shielding from the environment, as there are no phononic modes in the vacuum to exclude and their size compared to a microwave wavelength gives them a vanishingly small electromagnetic dipole moment.Finally, their reduced size and wholly different mechanism for generating nonlinearity may lead to reduced decoherence from sources such as charge fluctuations and two-level systems.
where in the second line, we use the fact that ∞ −∞ dte iωt = δ ω,0 .It is also worth noting the symbol P in front of the summation over the modes p 1 , ..., p N , which denotes a permutation over the mode indices.This is due to the fact that there is no specific requirement that each generic mode frequency ω rm equal the corresponding input frequency ω pm , but rather that the sum of the generic mode frequencies (i.e., ω r1 + ... + ω r N ) equal the sum of the input frequencies (i.e., ω p1 + ... + ω pn ).As a result, the set of transitions encapsulated by the real part of χ (N ) conserves the energy of the acoustic field, which is required in order for an electron starting in the ground state |g⟩ to return to that state after transitioning through the intermediate states |l 1 ⟩ , ..., |l N ⟩.Since each mode features a well-defined momentum as well, the net momentum of the field is conserved too.

Appendix B: Calculating Susceptibility Through Phase-Space Integrals
Here, we convert the summation over states in the susceptibility expressions to phase-space integrals.Given a 2DEG area A 2DEG , the following replacement can be used for any summation over generic initial wavevectors k: where the factor of 2 in the numerator is inserted to account for 2 spin states per spatial state.The main chal-lenge is to determine the boundaries of integration for k x and k y .The key requirement is an occupied initial state paired with an unoccupied final state.To this end, Fig. 19(a) depicts the initial and final states in phase space.As explained in the caption, the possible final states for a single excitation from the ground state is represented by the area enclosed by the red (right) circle, minus the shaded (overlapping) area.To represent the zone of forbidden initial states, the shaded zone representing the forbidden final states can be shifted leftward by q to cover the corresponding part of the blue (left) circle.The resulting range of valid initial states is depicted in Fig. 19(b).The states in the purple region (Zone 1) are forbidden from serving as initial states, while the rest of the states within the circle (green region) are allowed to act as initial states.We separate the range of allowed states into Zone 2 (left of dotted line) and Zone 3 (right of dotted line) in order to streamline the process of establishing the integral boundaries.In Zone 2, k x ranges from −q/2 to k F − q.For each given k x in this range, the span of k y is determined as follows: In Zone 3, k x ranges from k F − q to k F .The corresponding span of k y for each value of k x is the following: Consequently, the conversion of the summation over k into integral form, laid out generically in Eq. (B1), takes the following concrete form: (B4) We will evaluate this integral both numerically and analytically (in the high-mobility and low-mobility limits) to determine the susceptibility results.We start by examining the third-order Kerr suscepbility χ Kerr by applying Eq. ( 27) to the specific case of degenerate four wave mixing.The χ (3) susceptibility for this process features the input frequency triplet (ω, ω, −ω 1 ) (where ω and ω 1 are positive frequencies), yielding the following expression: where we have defined q = ω/v s and q 1 = ω 1 /v s , and the factor of 2 has been added in front to account for the permutations of ω.We convert the summation over k x and k y into integrals using the procedure from Eq. (B4).Note that the integral boundaries for all terms can be aligned if we flip the signs of all wavevector arguments for the frequency shift terms ω k f ,ki , which is valid since the energy spectrum and occupation probabilities are symmetric about the k x -axis.Then, χ (3) reduces to the following form: kx (ω, ω, −ω 1 ) where f (3) kx (ω, ω, −ω 1 ) is approximately the following: kx (ω, ω, −ω 1 ) ≈ . (B7) The third-order Kerr susceptibility equals the degenerate four-wave-mixing susceptibility in the limit ω 1 → ω.

(B10)
The susceptibility governing the second-harmonic generation process equals the above susceptibility for the general sum-frequency generation in the limit ω 1 = ω 2 .Finally, we use Eq. ( 27) once more to briefly derive the linear susceptibility χ (1) (ω) corresponding to the absorption and emission of a phonon of frequency ω: χ (1) (ω) = q 2 e ℏϵ 0 V 2DEG where q = ω/v s .Converting to an integral form over phase space in a manner analogous to the χ (3) and χ (2)   calculations, we find the following χ (1) : where f  where we have used the approximation t 2DEG ≪ l in the fourth line.The attenuation factor k takes the following form in terms of the dipole length and 2DEG thickness: The overall field E inside the 2DEG can be solved by substituting the relationship P = ϵ 0 χ (1) E into Eq.(C3): ϵ 0 E ext = ϵ 0 E + aϵ 0 χ (1) E, E = E ext 1 + aχ (1) ≈ E ext aχ (1) , (C7) where in the last line, we have used the approximation aχ (1) ≫ 1.

FIG. 7 :
FIG. 7: Phase-space diagram of near-resonance band spanning the electron-phonon detuning range ∆ω ≈ −γ/2 to ∆ω ≈ γ/2, zoomed out to show the position of the orange/light-colored band relative to the overall Fermi circle (a), and zoomed in to show the dimensions of each half of the band (b).
5(a)  and (d)), while for the

FIG. 17 :
FIG. 17: Diagram of the acoustic cavity, with the 2DEG (red) stacked on top of the piezoelectric material (blue) in the central region.The coupling between the two gives rise to a Kerr nonlinearity through 4-wave-mixing.

FIG. 19 :
FIG. 19: (a) Phase-space diagram of valid states for absorption of a phonon with wavevector xq by an electron.The blue (left) circle is the Fermi circle, which in the ground state is completely filled.The red (right) circle is the Fermi circle shifted in the + kx -direction by q, with the enclosed area representing potential destinations for the excited electrons.However, the shaded (overlapping) area represents states that are banned from serving as final states, since they are already occupied.Consequently, the range of possible destination states are represented by the area within the red (right) circle, less the shaded (overlapping) area; (b) Phase-space diagram of the range of possible initial states.The purple region (Zone 1) represents the forbidden initial states, while the green region (Zones 2 and 3) represents the allowed initial states. k

( 1 ) 2 + 1 ω 2 .FIG. 20 :
FIG.20: Depiction of induced electric field in the 2DEG due to charge polarization effectively giving rise to 2 plates of charge density ±P separated along the x-axis by the dipole length in the material, which we label as l.