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

Eric Chatterjee  Wyant College of Optical Sciences, The University of Arizona,
Tucson, Arizona 85721, USA  Alexander Wendt  Wyant College of Optical
Sciences, The University of Arizona, Tucson, Arizona 85721, USA  Daniel Soh
<EMAIL_ADDRESS>Wyant College of Optical Sciences, The University of
Arizona, Tucson, Arizona 85721, USA  Matt Eichenfield<EMAIL_ADDRESS>Wyant College of Optical Sciences, The University of Arizona, Tucson, Arizona
85721, USA  Sandia National Laboratories, Albuquerque, New Mexico 87123, USA

###### Abstract

Solid-state elastic-wave phonons are a promising platform for a wide range of
quantum information applications including facilitating quantum transduction
from microwave to optical electromagnetic fields, long-lived quantum memories,
in addition to potentially acting as qubits themselves. An outstanding
challenge and enabling capability in harnessing phonons for quantum
information processing is achieving sufficiently strong nonlinear interactions
between them. To this end, we propose a general architecture for generating
strong quantum phononic nonlinearity using piezoelectric-semiconductor
heterostructures consisting of a piezoelectric acoustic material hosting
phononic modes that is 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 electric fields associated with the various
phonon modes. We derive, from first principles, the nonlinear phononic
susceptibility of a piezo-2DEG system and provide a prescription to calculate
all orders of the susceptibility in a perturbative expansion. 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 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
phononic nonlinearities generated 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, the combination of 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.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED111Distribution Statement A: Approved for Public
Release, Distribution Unlimited

##  I Introduction

Quantized elastic waves in solids, or phonons, have become increasingly
ubiquitous in quantum information processing applications. For example, they
can act as a nearly “universal” quantum bus Raniwala  _et al._ (2023); Neuman
_et al._ (2021), allowing coupling to many useful solid state qubit
modalities. Phonons naturally couple to optomechanical cavities Eichenfield
_et al._ (2009); Aspelmeyer  _et al._ (2014); Kippenberg and Vahala (2007);
Lemonde  _et al._ (2013); Metcalfe (2014), which provide a means for optical
generation of nonclassical phonon states Riedinger  _et al._ (2016); Hong  _et
al._ (2017); Li  _et al._ (2018); Zivari  _et al._ (2022) and optical
teleportation of quantum information Fiaschi  _et al._ (2021); Soh  _et al._
(2021). Moreover, piezoelectric elastic wave phonons can be used to
efficiently couple to traveling wave Hackett  _et al._ (2023a); Safavi-Naeini
and Painter (2011); Massel  _et al._ (2011); Wu  _et al._ (2020); Arnold  _et
al._ (2020) and standing wave Teufel  _et al._ (2008); Safavi-Naeini  _et al._
(2019) microwave fields in important quantum phononic materials, couple to
superconducting circuits Chamberland  _et al._ (2022); MacCabe  _et al._
(2020); Arrangoiz-Arriola  _et al._ (2018); Barzanjeh  _et al._ (2022);
Viennot  _et al._ (2018); Satzinger  _et al._ (2018), or potentially act as
quantum memories to store quantum information from superconducting circuits
due to their long lifetime Taylor  _et al._ (2022); Kilina  _et al._ (2011);
Hann  _et al._ (2019); Wallucks  _et al._ (2020); Chamberland  _et al._
(2022); MacCabe  _et al._ (2020); Chan  _et al._ (2011). If those
piezoelectric phonons are stored within optomechanical cavities, then they
provide a means for optical addressing of quantum information in
superconducting qubits Bochmann  _et al._ (2013); Mirhosseini  _et al._
(2020); Forsch  _et al._ (2020); Jiang  _et al._ (2020); Balram and Srinivasan
(2022). Phonons also can couple to solid state color centers and their spin
degrees of freedom Kepesidis  _et al._ (2013); Albrecht  _et al._ (2013a, b);
Bennett  _et al._ (2013); Teissier  _et al._ (2014); Goldman  _et al._ (2015);
Meesala  _et al._ (2016); Li  _et al._ (2016); Oeckinghaus  _et al._ (2020);
Liu and Cundiff (2020); Neuman  _et al._ (2021).

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

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 single-phonon 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 Hackett  _et al._ (2023b). Classical theories of such
nonlinearities Conwell and Ganguly (1971); Wu and Spector (1972); Johri and
Spector (1975) 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
Pfeiffer  _et al._ (1989); Umansky  _et al._ (2009).

To this end, we study a heterostructure consisting of a 2DEG bonded to a
strongly piezoelectric elastic slab waveguide substrate (e.g. thin-film LiNbO3
on Si, sapphire, or SiC or thin-film Al1-xScxN 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 Hackett  _et
al._ (2021, 2023c, 2023a, 2023d), phononic switches Storey  _et al._ (2021),
phononic circulators Hackett  _et al._ (2021), phononic mixing Hackett  _et
al._ (2023b), and acoustoelectrically Brillouin optomechanical processes
Otterstrom  _et al._ (2023). 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 stemming
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 Shao
_et al._ (2019). 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) $\omega_{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\omega_{0}$ ($2k_{0})$ via acoustic evanescent coupling. The
amplified phononic field at $\omega_{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 Siddiqui  _et al._ (2018); Liu  _et al._ (2019);
Eichenfield and Olsson (2013), phononic negative refraction in phononic
crystals Zhang and Liu (2004); Page (2016), and phononic crystal gradient-
index (GRIN) lenses Wu  _et al._ (2011). 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 Schuetz  _et al._ (2015):

$\displaystyle T_{ij}$ $\displaystyle=\sum_{\begin{subarray}{c}k,l\\\ l\geq
k\end{subarray}}\Big{(}c_{ijkl}S_{kl}-e_{kij}E_{k}\Big{)},$ (1) $\displaystyle
D_{i}$ $\displaystyle=\sum_{\begin{subarray}{c}j,k\\\ k\geq
j\end{subarray}}\Big{(}\epsilon_{ij}E_{j}+e_{ijk}S_{jk}\Big{)},$ (2)

where $c$, $\epsilon$, 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
$\bm{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):

$E_{i}=-\sum_{\begin{subarray}{c}j,k,l\\\ l\geq
k\end{subarray}}(\epsilon^{-1})_{ij}e_{jkl}S_{kl}.$ (3)

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:

$\displaystyle\begin{split}\bm{u}(\bm{r})&=\bm{u}_{n,q}f_{n,q}(z)\Big{(}e^{i(qx-\omega
t)}+e^{-i(qx-\omega t)}\Big{)}\\\ &=2\bm{u}_{n,q}f_{n,q}(z)\cos{(qx-\omega
t)}.\end{split}$ (4)

where the functions $f_{n,q}$ are orthogonal for a given wavevector $q$:

$\int dzf_{n,q}(z)f_{m,q}(z)=L_{z}\delta_{m,n}.$ (5)

The strain components corresponding to the displacement $\bm{u}(\bm{r})$ are
calculated as follows:

$\displaystyle\begin{split}S_{ij}(\bm{r})&=\frac{\partial u_{i}}{\partial
r_{j}}(\bm{r})+(1-\delta_{i,j})\frac{\partial u_{j}}{\partial
r_{i}}(\bm{r})\\\ &=u_{n,q}\Big{(}g_{ij,n,q}(z)e^{i(qx-\omega
t)}+g_{ij,n,q}^{*}e^{-i(qx-\omega t)}\Big{)},\end{split}$ (6)

where $g_{ij,n,q}$ is defined as follows:

$\displaystyle\begin{split}&g_{ij,n,q}(z)=\\\ &\quad
iqf_{n,q}(z)\Big{(}\cos{\theta_{n,q,i}}\delta_{j,x}+\cos{\theta_{n,q,j}}(1-\delta_{i,j})\delta_{i,x}\Big{)}\\\
&\quad+f^{\prime}_{n,q}(z)\Big{(}\cos{\theta_{n,q,i}}\delta_{j,z}+\cos{\theta_{n,q,j}}(1-\delta_{i,j})\delta_{i,z}\Big{)},\end{split}$
(7)

where $\theta_{n,q,m}$ is defined as the angle between the displacement vector
$\bm{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 $\bm{\dot{u}}$ to the kinetic energy $T$ and use
the fact that the total mode energy $N\hbar\omega$ is twice the kinetic energy
for a simple harmonic oscillator:

$\displaystyle\begin{split}&\bigg{(}N+\frac{1}{2}\bigg{)}\hbar\omega\\\
&=2T\\\ &=2\int d^{3}r\frac{1}{2}\rho\sum_{i}\dot{u_{i}}^{2}(\bm{r})\\\
&=4\rho\omega^{2}u_{n,q}^{2}\int dy\int dx\sin^{2}{(qx-\omega t)}\int
dzf_{n,q}^{2}(z)\\\ &=4\rho\omega^{2}u_{n,q}^{2}L_{y}\frac{L_{x}}{2}L_{z},\\\
u_{n,q}&=\sqrt{\frac{(N+1/2)\hbar}{2\rho\omega V}},\end{split}$ (8)

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 $\omega$ and wavevector
$q$ propagating in the $\hat{m}$-direction can be expanded in a manner
analogous to the displacement and strain amplitudes:

$\bm{E}(\bm{r})=\bm{E}_{n,q}(z)e^{i(qx_{m}-\omega
t)}+\bm{E^{*}}_{n,q}(z)e^{-i(qx_{m}-\omega t)}.$ (9)

We solve for the $z$-dependent electric field amplitude $\bm{E}_{n,q}$ as a
function of the phonon numbers in the respective modes using Eq. (3) as
follows:

$\displaystyle\begin{split}&E_{n,q,i}(z)\\\
&=-\sum_{\begin{subarray}{c}j,k,l\\\ l\geq
k\end{subarray}}(\epsilon^{-1})_{ij}e_{jkl}u_{n,q}g_{kl,n,q}(z)\\\
&=-\sqrt{\frac{(N+1/2)\hbar}{2\rho\omega V}}\sum_{\begin{subarray}{c}j,k,l\\\
l\geq
k\end{subarray}}(\epsilon^{-1})_{ij}e_{jkl}\bigg{(}iqf_{n,q}(z)\Big{(}\cos{\theta_{n,q,k}}\delta_{k,x}+\cos{\theta_{n,q,l}}(1-\delta_{k,l})\delta_{k,x}\Big{)}\\\
&\quad\quad+f^{\prime}_{n,q}(z)\Big{(}\cos{\theta_{n,q,k}}\delta_{l,z}+\cos{\theta_{n,q,l}}(1-\delta_{k,l})\delta_{k,z}\Big{)}\bigg{)}.\end{split}$
(10)

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)\propto
e^{-qz}$. Then, given a diagonal dielectric tensor (i.e.,
$\epsilon_{ij}=\epsilon_{ii}\delta_{i,j}$), the electric field amplitude at
the surface (which we define as $z=0$) reduces to the following:

$E=-\frac{C}{\epsilon}\sqrt{\frac{\hbar\omega}{v_{s}}}A_{\mathrm{ph}},$ (11)

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/\sqrt{\kappa}$, where $e$ and $\kappa$ are the
piezoelectric coupling and elasticity coefficients, respectively), $\epsilon$
is the effective dielectric constant, and $A_{\mathrm{ph}}$ is a measurement
of the square-root of the propagating strain field power density (i.e., power
per unit area) in units of $\hbar\omega$. Each of these parameters is defined
in the following manner:

$\displaystyle C$
$\displaystyle=\frac{f_{n,q}(0)}{\sqrt{2v_{s}^{2}\rho}}\sum_{\begin{subarray}{c}k,l\\\
l\geq
k\end{subarray}}e_{ikl}\bigg{(}i\Big{(}\cos{\theta_{n,q,k}}\delta_{k,x}+\cos{\theta_{n,q,l}}(1-\delta_{k,l})\delta_{k,x}\Big{)}-\Big{(}\cos{\theta_{n,q,k}}\delta_{l,z}+\cos{\theta_{n,q,l}}(1-\delta_{k,l})\delta_{k,z}\Big{)}\bigg{)},$
(12) $\displaystyle\epsilon$ $\displaystyle=\epsilon_{xx},$ (13)
$\displaystyle A_{\mathrm{ph}}$ $\displaystyle=\sqrt{v_{s}\frac{(N+1/2)}{V}}.$
(14)

For the electron-phonon interaction, we use the dielectric tensor for the 2DEG
as $\epsilon$, instead of the dielectric tensor for the piezoelectric
material. In Sec. V, we will use $A_{\mathrm{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^{\mathrm{th}}$-order susceptibility of the 2DEG electrons as a
function of the input electric field amplitudes. By definition, the
$N^{\mathrm{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^{\mathrm{th}}$-order susceptibility.

###  III.1 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_{\mathrm{eff}}$ and the electric field
amplitude as $E$, we apply the standard dipole interaction Hamiltonian:

$H_{\mathrm{int}}=-d_{\mathrm{eff}}E.$ (15)

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:

$H_{\mathrm{int}}=-\int d^{3}rE_{x}P_{x},$ (16)

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^{({\dagger})}$:

$E_{x}=E_{\mathrm{zpf}}be^{iqx}+E^{*}_{\mathrm{zpf}}b^{{\dagger}}e^{-iqx},$
(17)

where $E_{\mathrm{zpf}}$ is the zero-point piezoelectric field amplitude. Note
that the spatial dependence of the field is entirely isolated to the phases
$e^{\pm 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_{y}$-direction, 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}+\Delta k,k_{y})$ Todorov (2015):

$\displaystyle\begin{split}P_{x}&=\frac{|\phi(z)|^{2}}{L_{y}}\sum_{k_{i},\Delta
k,k_{y}}\frac{q_{e}\hbar}{2m\omega_{k_{i}+\Delta k,k_{i}}}\\\
&\quad\times\Big{(}\xi_{k_{i}+\Delta k,k_{i}}(x)c_{k_{i}+\Delta
k,k_{y}}^{{\dagger}}c_{k_{i},k_{y}}\\\ &\quad\quad+\xi^{*}_{k_{i}+\Delta
k,k_{i}}(x)c_{k_{i},k_{y}}^{{\dagger}}c_{k_{i}+\Delta
k,k_{y}}\Big{)},\end{split}$ (18)

where $c^{({\dagger})}_{k_{x},k_{y}}$ is the fermionic annihilation (creation)
operator for the electronic state at $(k_{x},k_{y})$, $\omega_{k_{i}+\Delta
k,k_{i}}=\omega_{k_{i}+\Delta k}-\omega_{k_{i}}$ is the change in carrier
energy (this is independent of $k_{y}$), $\phi(z)$ is the wavefunction along
the $z$ axis, $L_{y}$ is the span of the 2DEG along the $y$-axis, and $\xi$ is
defined as follows:

$\displaystyle\begin{split}&\xi_{k_{i}+\Delta k,k_{i}}(x)\\\
&=\frac{1}{L}\Big{(}e^{-i(k_{i}+\Delta
k)x}\partial_{x}e^{ik_{i}x}-(\partial_{x}e^{-i(k_{i}+\Delta
k)x})e^{ik_{i}x}\Big{)}\\\ &=i\frac{2k_{i}+\Delta k}{L}e^{-i(\Delta
k)x},\end{split}$ (19)

where $L$ is the waveguide length. It is worth noting that
$\omega_{k_{i}+\Delta k,k_{i}}$ can be solved as follows:

$\displaystyle\begin{split}\omega_{k_{i}+\Delta k,k_{i}}&=\omega_{k_{i}+\Delta
k}-\omega_{k_{i}}\\\ &=\frac{\hbar}{2m}\Big{(}(k_{i}+\Delta
k)^{2}-k_{i}^{2}\Big{)}\\\ &=\frac{\hbar\Delta k(2k_{i}+\Delta
k)}{2m}.\end{split}$ (20)

Substituting the above expressions into the Hamiltonian in Eq. (16), we
express the Hamiltonian in electron-phonon interaction form as follows:

$\displaystyle\begin{split}&H_{\mathrm{int}}\\\ &=\int
dx\Big{(}E_{\mathrm{zpf}}be^{iqx}+E^{*}_{\mathrm{zpf}}b^{{\dagger}}e^{-iqx}\Big{)}\sum_{k_{i},\Delta
k,k_{y}}i\frac{q_{e}}{L\Delta k}\Big{(}e^{-i(\Delta k)x}c_{k_{i}+\Delta
k,k_{y}}^{{\dagger}}c_{k_{i},k_{y}}-e^{i(\Delta
k)x}c_{k_{i},k_{y}}^{{\dagger}}c_{k_{i}+\Delta k,k_{y}}\Big{)}\int
dz|\phi(z)|^{2}\int\frac{dy}{L_{y}}\\\ &=i\sum_{k_{i},\Delta
k,k_{y}}\frac{q_{e}}{L\Delta k}\Big{(}E_{\mathrm{zpf}}bc_{k_{i}+\Delta
k,k_{y}}^{{\dagger}}c_{k_{i},k_{y}}\int dxe^{i(q-\Delta
k)x}+E^{*}_{\mathrm{zpf}}b^{{\dagger}}c_{k_{i},k_{y}}^{{\dagger}}c_{k_{i}+\Delta
k,k_{y}}\int dxe^{-i(q+\Delta k)x}\Big{)}\\\
&=\frac{q_{e}E_{\mathrm{zpf}}}{q}(c_{k_{i}+q,k_{i}}^{{\dagger}}b\delta_{q,\Delta
k}+b^{{\dagger}}c_{k_{i}+q,k_{i}}\delta_{q,-\Delta k}).\end{split}$ (21)

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_{{\mathrm{eff}},k_{x}+\Delta k,k_{x}}(q)=-\frac{q_{e}}{q}\delta_{q,\Delta
k},$ (22)

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.

Figure 2: Conceptual diagram of how the alternating electric field $\bm{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.

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.

###  III.2 Deriving the higher-order interaction probabilities

We now seek to determine the higher-order electron-phonon interaction
probabilities, which yields the polarization field as a function of the
acoustic (and thus piezoelectric) field strength, by deriving the generic
$N^{\textrm{th}}$-order susceptibility in terms of the general electronic
decay rate $\gamma=q_{e}/(m\mu)$, and from this deriving the conversion
efficiency for specific nonlinear processes. As in the case of an atom
interacting with an optical field Boyd (2008), 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 $\gamma$ as an anti-Hermitian effective Hamiltonian
term Chatterjee  _et al._ (2022), yielding the following unperturbed
Hamiltonian for the 2DEG states:

$H_{0}=\sum_{l>g}\hbar\bigg{(}\omega_{lg}-i\frac{\gamma_{l}}{2}\bigg{)}\ket{l}\bra{l},$
(23)

where $\omega_{lg}=\omega_{l}-\omega_{g}$ represents the energy of the excited
state $\ket{l}$ (with the energy of the ground state $\ket{g}$ set as the
baseline), and $\gamma_{l}=\gamma$ for all excited states $\ket{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 $+\hat{x}$-direction. The fields interact with the
electron dipole moment $\bm{d}$ via the following time-dependent interaction
Hamiltonian:

$\displaystyle\begin{split}&H_{\mathrm{int}}(t)\\\
&=-d_{x}\sum_{p(\omega_{p}>0)}\Big{(}E(\omega_{p})e^{i(q_{p}x-\omega_{p}t)}+E^{*}(\omega_{p})e^{-i(q_{p}x-\omega_{p}t)}\Big{)}\\\
&=-\sum_{p}d_{\mathrm{eff}}(q_{p})E(\omega_{p})e^{-i\omega_{p}t},\end{split}$
(24)

where $q_{p}=\omega_{p}/v_{s}$, and the spatial variation of the field is
incorporated into the effective dipole moment $d_{\mathrm{eff}}$ (which is
defined as in Eq. (22)). Note that $E(-\omega_{p})=E^{*}(\omega_{p})$,
ensuring the realness of the composite field at all points. We can also
express a generic field of frequency $\omega_{b}>0$ in terms of the
corresponding ladder operators $a_{b}^{({\dagger})}$ as
$E(\omega_{b})=E_{\mathrm{zpf}}(\omega_{b})a_{b}$ and
$E(-\omega_{b})=E^{*}(\omega_{b})=E_{\mathrm{zpf}}^{*}(\omega_{b})a_{b}^{{\dagger}}$.
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, $\left\langle d_{\mathrm{eff}}\right\rangle$, 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 $\left\langle
H_{\mathrm{int}}\right\rangle$ terms), we express this expectation value as
follows:

$\displaystyle\begin{split}&\left\langle
d_{\mathrm{eff}}(q_{p})\right\rangle\\\
&=\epsilon_{0}V_{\mathrm{2DEG}}\sum_{n}\sum_{\begin{subarray}{c}p_{1},...,p_{N}\\\
\omega_{p_{1}}+...+\omega_{p_{N}}=-\omega_{p}\end{subarray}}\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{N}})\\\
&\quad\quad\times
E(\omega_{p_{1}})...E(\omega_{p_{N}})e^{i\omega_{p}t},\end{split}$ (25)

where the $N^{\textrm{th}}$-order susceptibility $\chi^{(n)}$ is defined as a
function proportional to the $N^{\textrm{th}}$-order expansion of
$\left\langle d_{\mathrm{eff}}\right\rangle$, in the case where a unique set
of frequencies $\omega_{p_{1}},...,\omega_{p_{n}}$ add up to $-\omega_{p}$:

$\displaystyle\begin{split}\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{n}})&=\frac{\int_{-\infty}^{\infty}\left\langle
d_{\mathrm{eff}}^{(N)}(-q_{p_{1}}-...-q_{p_{N}})\right\rangle_{t}e^{i(\omega_{p_{1}}+...+\omega_{p_{N}})t}}{\epsilon_{0}V_{\mathrm{2DEG}}E(\omega_{p_{1}})...E(\omega_{p_{N}})}\\\
&=\frac{\int_{-\infty}^{\infty}e^{i(\omega_{p_{1}}+...+\omega_{p_{N}})t}\sum_{j=0}^{N}\braket{\Psi^{(j)}(t)}{d_{\mathrm{eff}}(-q_{p_{1}}-...-q_{p_{N}})}{\Psi^{(N-j)}(t)}}{\epsilon_{0}V_{\mathrm{2DEG}}E(\omega_{p_{1}})...E(\omega_{p_{N}})},\end{split}$
(26)

where $\ket{\psi^{(n)}(t)}$ is the $n^{\textrm{th}}$-order expansion of the
composite 2DEG wavefunction $\ket{\Psi(t)}$. Note that the variation of
$\left\langle d_{\mathrm{eff}}\right\rangle$ with the field amplitudes in Eq.
(LABEL:eq:_dipole_moment_expectation_value) 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,
$\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{n}})$ corresponds to the probability
that an electron becomes polarized upon interacting with $N$ fields of
wavevectors $q_{p_{1}},...,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 $\chi^{(N)}$ is calculated for a set of $N$ propagating
acoustic waves with wavevectors $q_{p_{1}},...,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 $\gamma$, the energy conservation requirement
is somewhat relaxed. This dynamic is depicted in Fig. 3.

Figure 3: Depiction of the ladder of intermediate states for the $N$-step
transition described by $\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{N}})$, with
the horizontal (vertical) axis representing the wavevector $q$ (frequency
$\omega$). 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
$\omega_{r_{1}}+...+\omega_{r_{N}}=\omega_{p_{1}}+...+\omega_{p_{N}}$. Also
note that each intermediate state is broadened in energy but fixed in
momentum.

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
$\omega_{r_{1}}+...+\omega_{r_{N}}=\omega_{p_{1}}+...+\omega_{p_{N}}$, and
hence $q_{r_{1}}+...+q_{r_{N}}=q_{p_{1}}+...+q_{p_{N}}$. The transition to
(from) $\ket{l_{m}}$ is mediated by the mode $r_{m}$ for $m\leq N-j$
($m>N-j$). The transition from $\ket{l_{N-j}}$ to $\ket{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
$\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{N}})$.

The most critical takeaway from Fig. 3 is that while the intermediate 2DEG
states are broadened in frequency $\omega$ (corresponding to a wide range of
possible real states), each intermediate state is fixed in wavevector $q$,
since the effective dipole moment $d_{\mathrm{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. (123) thus reduces to a sum over initial
states (which we index with the 2D wavevector $(k_{x},k_{y})$):

$\displaystyle\begin{split}&\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{N}})\\\
&=\frac{1}{\hbar^{N}\epsilon_{0}V_{\mathrm{2DEG}}}\frac{q_{e}^{N+1}}{|q_{p_{1}}...q_{p_{n}}(q_{p_{1}}+...+q_{p_{n}})|}\sum_{j=0}^{N}\sum_{k_{x},k_{y}}\mathcal{P}\sum_{p_{1},...,p_{N}}\\\
&\quad\frac{1}{(\omega_{k_{x}-q_{p_{N}},k_{x}}+\omega_{p_{N}}+i\frac{\gamma}{2})...(\omega_{k_{x}-q_{p_{N}}-...-q_{p_{N-j+1}},k_{x}}+(\omega_{p_{N-j+1}}+...+\omega_{p_{N}})+i\frac{\gamma}{2})}\\\
&\quad\times\frac{1}{(\omega_{k_{x}+q_{p_{1}}+...+q_{p_{N-j}},k_{x}}-(\omega_{p_{1}}+...+\omega_{p_{N-n}})-i\frac{\gamma}{2})...(\omega_{k_{x}+q_{p_{1}},k_{x}}-\omega_{p_{1}}-i\frac{\gamma}{2})},\end{split}$
(27)

where we have used the fact that
$d_{\mathrm{eff},k_{f},k_{i}}(q)=(q_{e}/|q|)\delta_{q,k_{f}-k_{i}}$, and the
symbol $\mathcal{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!$ permutations 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
Shao  _et al._ (2019)), the spectral broadening is dominated by the phonon
absorption by the 2DEG electrons, for which we need to determine the imaginary
part of $\chi^{(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\ll k_{F}$ (where
$q$ and $k_{F}$ are the acoustic and Fermi wavevectors, respectively).

First, we lay out the generic procedure for reducing the
$N^{\textrm{th}}$-order susceptibility to a calculable form. Given a 2DEG area
$A_{\mathrm{2DEG}}$, the following replacement can be used for any summation
over generic initial wavevectors $\bm{k}$:

$\sum_{\bm{k}}=\frac{2A_{\mathrm{2DEG}}}{(2\pi)^{2}}\int dk_{x}\int dk_{y},$
(28)

where the factor of 2 in the numerator is inserted to account for 2 spin
states per spatial state. The main challenge 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. 4(a)
depicts the initial and final states in phase space.

(a)

(b)

Figure 4: (a) Phase-space diagram of valid states for absorption of a phonon
with wavevector $\bm{\hat{x}}q$ 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 $+\hat{k_{x}}$-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.

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. 4(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:

$2\bigg{(}\sqrt{k_{F}^{2}-k_{x}^{2}}-\sqrt{k_{F}^{2}-(q+k_{x})^{2}}\bigg{)}.$
(29)

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:

$2\sqrt{k_{F}^{2}-k_{x}^{2}}.$ (30)

Consequently, the conversion of the summation over $\bm{k}$ into integral
form, laid out generically in Eq. (28), takes the following concrete form:

$\displaystyle\begin{split}\sum_{\bm{k}}&=\frac{2A_{\mathrm{2DEG}}}{(2\pi)^{2}}\Bigg{(}\int_{-q/2}^{k_{F}-q}dk_{x}2\bigg{(}\sqrt{k_{F}^{2}-k_{x}^{2}}-\sqrt{k_{F}^{2}-(q+k_{x})^{2}}\bigg{)}+\int_{k_{F}-q}^{k_{F}}dk_{x}2\sqrt{k_{F}^{2}-k_{x}^{2}}\Bigg{)}\\\
&=\frac{A_{\mathrm{2DEG}}}{2\pi^{2}}\Bigg{(}\int_{-q/2}^{k_{F}}dk_{x}2\sqrt{k_{F}^{2}-k_{x}^{2}}-\int_{-q/2}^{k_{F}-q}2\sqrt{k_{F}^{2}-(q+k_{x})^{2}}\Bigg{)}.\end{split}$
(31)

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
$\chi^{(3)}_{\mathrm{Kerr}}$ by applying Eq. (LABEL:eq:_chi(N)_reduced) to the
specific case of degenerate four wave mixing. The $\chi^{(3)}$ susceptibility
for this process features the input frequency triplet
$(\omega,\omega,-\omega_{1})$ (where $\omega$ and $\omega_{1}$ are positive
frequencies), yielding the following expression:

$\displaystyle\begin{split}&\chi^{(3)}(\omega,\omega,-\omega_{1})\\\
&=2\frac{q_{e}^{4}}{\hbar^{3}\epsilon_{0}V_{\mathrm{2DEG}}}\frac{v_{s}^{4}}{\omega^{2}\omega_{1}(2\omega-\omega_{1})}\sum_{k_{x},k_{y}}\bigg{(}\frac{1}{(\omega_{k_{x}+2q-q_{1},k_{x}}-(2\omega-\omega_{1})-i\frac{\gamma}{2})(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q,k_{x}}-\omega-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}+q_{1},k_{x}}-\omega_{1}+i\frac{\gamma}{2})(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q,k_{x}}-\omega-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}-q,k_{x}}+\omega+i\frac{\gamma}{2})(\omega_{k_{x}-2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}-q_{1},k_{x}}+\omega_{1}-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}-q,k_{x}}+\omega+i\frac{\gamma}{2})(\omega_{k_{x}-2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}-2q+q_{1},k_{x}}+(2\omega-\omega_{1})+i\frac{\gamma}{2})}\bigg{)},\end{split}$
(32)

where we have defined $q=\omega/v_{s}$ and $q_{1}=\omega_{1}/v_{s}$, and the
factor of 2 has been added in front to account for the permutations of
$\omega$. The first two terms are resonant, corresponding to the absorption of
2 phonons of frequency $\omega$ each, followed by the emission of a phonon of
frequency $\omega_{1}$ and another of frequency
$\omega_{2}=2\omega-\omega_{1}$. On the other hand, the latter two terms are
counter-resonant, corresponding to the emission of $\omega_{1}$ and
$\omega_{2}$ phonons followed by absorption of two $\omega$ phonons. The
process described by each term is depicted in Fig. 6.

(a)

(b)

(a)

(b)

Figure 6: Phase-space depiction of processes contributing to
$\chi^{(3)}(\omega,\omega,-\omega_{1})$, with (a), (b), (c), and (d)
corresponding respectively to the first, second, third, and fourth terms of
Eq. (LABEL:eq:_chi(3)(omega,omega,-omega1)_summation). 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.

It is worth noting that if we ignored the initial occupation of the states,
then there would be 12 separate processes contributing to
$\chi^{(3)}(\omega,\omega,-\omega_{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).

Finally, we convert the summation over $k_{x}$ and $k_{y}$ into integrals
using the procedure from Eq. (31). 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 $\omega_{k_{f},k_{i}}$, which is valid since the energy
spectrum and occupation probabilities are symmetric about the $k_{x}$-axis.
Then, $\chi^{(3)}$ reduces to the following form:

$\displaystyle\begin{split}&\chi^{(3)}(\omega,\omega,-\omega_{1})=\\\
&\quad\frac{q_{e}^{4}}{\pi^{2}\hbar^{3}\epsilon_{0}t_{\mathrm{2DEG}}}\frac{v_{s}^{4}}{\omega^{2}\omega_{1}(2\omega-\omega_{1})}\\\
&\quad\times\Bigg{(}\int_{-q/2}^{k_{F}}dk_{x}2\sqrt{k_{F}^{2}-k_{x}^{2}}f_{k_{x}}^{(3)}(\omega,\omega,-\omega_{1})\\\
&\quad-\int_{-q/2}^{k_{F}-q}dk_{x}2\sqrt{k_{F}^{2}-(q+k_{x})^{2}}f_{k_{x}}^{(3)}(\omega,\omega,-\omega_{1})\Bigg{)},\end{split}$
(33)

where $f_{k_{x}}^{(3)}(\omega,\omega,-\omega_{1})$ is approximately the
following:

$\displaystyle\begin{split}f_{k_{x}}^{(3)}(\omega,\omega,-\omega_{1})&\approx\frac{1}{(\omega_{k_{x}+q_{2},k_{x}}-\omega_{2}-i\frac{\gamma}{2})(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q,k_{x}}-\omega-i\frac{\gamma}{2})}\\\
&\quad+\frac{1}{(\omega_{k_{x}+q_{1},k_{x}}-\omega_{1}+i\frac{\gamma}{2})(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q,k_{x}}-\omega-i\frac{\gamma}{2})}\\\
&\quad+\frac{1}{(\omega_{k_{x}+q,k_{x}}+\omega+i\frac{\gamma}{2})(\omega_{k_{x}+2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}+\omega_{1}-i\frac{\gamma}{2})}\\\
&\quad+\frac{1}{(\omega_{k_{x}+q,k_{x}}+\omega+i\frac{\gamma}{2})(\omega_{k_{x}+2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}+\omega_{2}+i\frac{\gamma}{2})}.\end{split}$
(34)

The third-order Kerr susceptibility equals the degenerate four-wave-mixing
susceptibility in the limit $\omega_{1}\rightarrow\omega$.

Next, we examine the second-order susceptibility, focusing first on the case
of sum-frequency generation encapsulated by
$\chi^{(2)}(\omega_{1},\omega_{2})$, where two phonons of frequency
$\omega_{1}$ and $\omega_{2}$ are absorbed and a phonon of frequency
$2\omega=\omega_{1}+\omega_{2}$ is emitted. Applying Eq.
(LABEL:eq:_chi(N)_reduced) again, we obtain the following expression for the
second-order susceptibility:

$\displaystyle\begin{split}&\chi^{(2)}(\omega_{1},\omega_{2})\\\
&=\frac{q_{e}^{3}}{\hbar^{2}\epsilon_{0}V_{\mathrm{2DEG}}}\frac{v_{s}^{3}}{\omega_{1}\omega_{2}(\omega_{1}+\omega_{2})}\sum_{k_{x},k_{y}}\bigg{(}\frac{1}{(\omega_{k_{x}+q_{1}+q_{2},k_{x}}-(\omega_{1}+\omega_{2})-i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}-\omega_{1}-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}+q_{1}+q_{2},k_{x}}-(\omega_{1}+\omega_{2})-i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}-\omega_{2}-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}-q_{1}-q_{2},k_{x}}+(\omega_{1}+\omega_{2})+i\frac{\gamma}{2})(\omega_{k_{x}-q_{2},k_{x}}+\omega_{2}+i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}-q_{1}-q_{2},k_{x}}+(\omega_{1}+\omega_{2})+i\frac{\gamma}{2})(\omega_{k_{x}-q_{1},k_{x}}+\omega_{1}+i\frac{\gamma}{2})}\bigg{)},\end{split}$
(35)

where $q_{1}=\omega_{1}/v_{s}$ and $q_{2}=\omega_{2}/v_{s}$. The first two
terms are resonant and the latter two terms are counter-resonant, with the
processes depicted in Fig. 8.

(a)

(b)

(a)

(b)

Figure 8: Phase-space depiction of processes contributing to
$\chi^{(2)}(\omega_{1},\omega_{2})$, with (a), (b), (c), and (d) corresponding
respectively to the first, second, third, and fourth terms of Eq.
(LABEL:eq:_chi(2)(omega_1,omega_2)_summation). The shading scheme is identical
to that in Fig. 6.

As with the case of $\chi^{(3)}$, we apply the high-carrier-density limit
$|q|\ll 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. Converting the summation over wavevectors to
integrals using the same procedure and symmetry arguments as in the
$\chi^{(3)}$ calculation, we find that $\chi^{(2)}$ reduces to the following
form:

$\displaystyle\begin{split}&\chi^{(2)}(\omega_{1},\omega_{2})=\\\
&\quad\frac{q_{e}^{3}}{2\pi^{2}\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}}\frac{v_{s}^{3}}{\omega_{1}\omega_{2}(\omega_{1}+\omega_{2})}\\\
&\quad\times\Bigg{(}\int_{-q/2}^{k_{F}}dk_{x}2\sqrt{k_{F}^{2}-k_{x}^{2}}f_{k_{x}}^{(2)}(\omega_{1},\omega_{2})\\\
&\quad-\int_{-q/2}^{k_{F}-q}dk_{x}2\sqrt{k_{F}^{2}-(q+k_{x})^{2}}f_{k_{x}}^{(2)}(\omega_{1},\omega_{2})\Bigg{)},\end{split}$
(36)

where $f_{k_{x}}^{(2)}(\omega_{1},\omega_{2})$ is approximately the following:

$\displaystyle\begin{split}f_{k_{x}}^{(2)}(\omega_{1},\omega_{2})&\approx\bigg{(}\frac{1}{(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}-\omega_{1}-i\frac{\gamma}{2})}+\frac{1}{(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}-\omega_{2}-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}+2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}+\omega_{2}+i\frac{\gamma}{2})}+\frac{1}{(\omega_{k_{x}+2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}+\omega_{1}+i\frac{\gamma}{2})}\bigg{)}.\end{split}$
(37)

The susceptibility governing the second-harmonic generation process equals the
above susceptibility for the general sum-frequency generation in the limit
$\omega_{1}=\omega_{2}$.

In the context of second-order susceptibility, it is also worth examining the
case of parametric amplification encapsulatd by
$\chi^{(2)}(2\omega,-\omega_{1})$, where a pump phonon of frequency $2\omega$
is absorbed, causing one phonon each of the signal frequency $\omega_{1}$ and
the idler frequency $\omega_{2}$ to be emitted. Since the frequencies involved
are the same as the sum-frequency generation, the coefficients for
$\chi^{(2)}(2\omega,-\omega_{1})$ are the same as those for
$\chi^{(2)}(\omega_{1},\omega_{2})$ in Eq. (LABEL:eq:_chi(2)_integral).
However, $f_{k_{x}}^{(2)}(\omega_{1},\omega_{2})$ in the integral is replaced
by $f_{k_{x}}^{(2)}(2\omega,-\omega_{1})$, defined as:

$\displaystyle\begin{split}f_{k_{x}}^{(2)}(2\omega,-\omega_{1})&\approx\bigg{(}\frac{1}{(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}-\omega_{1}+i\frac{\gamma}{2})}+\frac{1}{(\omega_{k_{x}+2q,k_{x}}-2\omega-i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}-\omega_{2}-i\frac{\gamma}{2})}\\\
&\quad\quad+\frac{1}{(\omega_{k_{x}+2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}+\omega_{2}+i\frac{\gamma}{2})}+\frac{1}{(\omega_{k_{x}+2q,k_{x}}+2\omega+i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}+\omega_{1}-i\frac{\gamma}{2})}\bigg{)}.\end{split}$
(38)

Note the sign flip (relative to Eq. (37)) 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
($\omega_{1}=\omega_{2}$), as we will show in Sec. IV.4.

Finally, we use Eq. (LABEL:eq:_chi(N)_reduced) once more to briefly derive the
linear susceptibility $\chi^{(1)}(\omega)$ corresponding to the absorption and
emission of a phonon of frequency $\omega$:

$\displaystyle\begin{split}\chi^{(1)}(\omega)&=\frac{q_{e}^{2}}{\hbar\epsilon_{0}V_{\mathrm{2DEG}}}\frac{v_{s}^{2}}{\omega^{2}}\sum_{k_{x},k_{y}}\bigg{(}\frac{1}{\omega_{k_{x}+q,k_{x}}-\omega-i\frac{\gamma}{2}}\\\
&\quad\quad+\frac{1}{\omega_{k_{x}-q,k_{x}}+\omega+i\frac{\gamma}{2}}\bigg{)},\end{split}$
(39)

where $q=\omega/v_{s}$. The first term is resonant, while the second is
counter-resonant, as shown in Fig. 9.

(a)

(b)

Figure 9: Phase-space depiction of processes contributing to
$\chi^{(1)}(\omega)$, with (a) and (b) corresponding respectively to the first
and second terms of Eq. (39). The shading scheme is identical to that in Fig.
6.

Converting to an integral form over phase space in a manner analogous to the
$\chi^{(3)}$ and $\chi^{(2)}$ calculations, we find the following
$\chi^{(1)}$:

$\displaystyle\begin{split}&\chi^{(1)}(\omega)=\\\
&\quad\frac{q_{e}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}}\frac{v_{s}^{2}}{\omega^{2}}\Bigg{(}\int_{-q/2}^{k_{F}}dk_{x}2\sqrt{k_{F}^{2}-k_{x}^{2}}f_{k_{x}}^{(1)}(\omega)\\\
&\quad-\int_{-q/2}^{k_{F}-q}dk_{x}2\sqrt{k_{F}^{2}-(q+k_{x})^{2}}f_{k_{x}}^{(1)}(\omega)\Bigg{)},\end{split}$
(40)

where $f_{k_{x}}^{(1)}(\omega)$ is defined as follows:

$f_{k_{x}}^{(1)}(\omega)=\frac{1}{\omega_{k_{x}+q,k_{x}}-\omega-i\frac{\gamma}{2}}+\frac{1}{\omega_{k_{x}+q,k_{x}}+\omega+i\frac{\gamma}{2}}.$
(41)

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.

###  IV.1 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=\omega/v_{s}$ propagating in the $\hat{x}$-direction is governed by the
detuning $\Delta\omega$ between the phonon’s frequency and the transition
frequency corresponding to the electron jumping from $k_{x}\rightarrow
k_{x}+q$. Consequently, $\Delta\omega$ is a function of the electron’s initial
$x$-direction wavevector $k_{x}$ but is invariant in the initial $y$-direction
wavevector $k_{y}$:

$\Delta\omega=\omega_{k_{x}+q,k_{x}}-\omega_{0}=\frac{\hbar
q(2k_{x}+q)}{2m}-v_{s}q=\frac{\hbar q}{m}(k_{x}-k_{x,0}),$ (42)

where $\omega_{0}=v_{s}q$ denotes the phonon’s angular frequency, and
$k_{x,0}=mv_{s}/\hbar-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
$\Delta\omega\approx-\gamma/2$ to $\Delta\omega\approx\gamma/2$. In the high-
carrier-density/high-mobility regime, the $k_{x}$-span of this band is much
less than the diameter $2k_{F}$ of the Fermi circle. We solve for the
$k_{x}$-span of the near-resonance band from Eq. (42) as follows:

$k_{x,\textrm{span}}=\frac{m}{\hbar
q}(\Delta\omega)_{\textrm{span}}\approx\frac{2m\gamma}{\hbar q}.$ (43)

Conveniently, the $k_{x}$-span of the near-resonance band is proportional to
the electronic decay rate $\gamma$ (and hence inversely proportional to the
mobility $\mu$). Intuitively, this is due to the fact that the electron-phonon
detuning varies linearly with $k_{x}$. Applying the high-mobility requirement
$k_{x,\textrm{span}}\ll 2k_{F}$, we derive the following condition for
$\gamma$:

$\gamma\ll\frac{\hbar qk_{F}}{m},$ (44)

which corresponds to the following condition for the mobility $\mu$, upon
applying the definition $\mu=q_{e}/(m\gamma)$:

$\mu\gg\frac{q_{e}}{\hbar qk_{F}}.$ (45)

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}\ll k_{F}$), the $k_{y}$-span given in Eq. (29)
for $k_{x}$ values near $k_{x,0}$ can be approximately linearized in $k_{x}$
as follows:

$\displaystyle\begin{split}k_{y,\textrm{span}}&=2\Big{(}\sqrt{k_{F}^{2}-k_{x}^{2}}-\sqrt{k_{F}^{2}-(k_{x}+q)^{2}}\Big{)}\\\
&\approx
2k_{F}\bigg{(}-\frac{k_{x}^{2}}{2k_{F}^{2}}+\frac{(k_{x}+q)^{2}}{2k_{F}^{2}}\bigg{)}\\\
&=\frac{q(2k_{x}+q)}{k_{F}}.\end{split}$ (46)

We can express any given $k_{x}$ near $k_{x,0}$ as $k_{x}=k_{x,0}+\Delta
k_{x}$, where $\Delta k_{x}$ is the deviation of the initial wavevector from
that corresponding to a resonant electron-phonon interaction. Then, per Eq.
(42), the detuning $\Delta\omega$ varies linearly with the wavevector
deviation $\Delta k_{x}$ as $\Delta\omega=\hbar q(\Delta 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 $\Delta\omega$:

$\displaystyle\begin{split}k_{y,\textrm{span}}&\approx\frac{q(2k_{x,0}+q)}{k_{F}}+\frac{2q\Delta
k_{x}}{k_{F}}\\\ &=\frac{2m\omega_{0}}{\hbar
k_{F}}+\frac{2m\Delta\omega}{\hbar k_{F}},\end{split}$ (47)

where we have also substituted the phonon frequency $\omega_{0}=v_{s}q$. Using
the fact that $dk_{x}=md(\Delta\omega)/(\hbar
q)=mv_{s}d(\Delta\omega)/(\hbar\omega_{0})$, a generic susceptibility term
$f(\Delta\omega)$ can be integrated over phase space in the following manner:

$\displaystyle\begin{split}&\int dk_{x}\int dk_{y}f(\Delta\omega)\\\
&\approx\frac{2m}{\hbar
k_{F}}\int_{-\infty}^{\infty}\frac{mv_{s}}{\hbar\omega_{0}}d(\Delta\omega)(\omega_{0}+\Delta\omega)f(\Delta\omega)\\\
&=\frac{2m^{2}v_{s}}{\hbar^{2}k_{F}\omega_{0}}\int_{-\infty}^{\infty}d(\Delta\omega)(\omega_{0}+\Delta\omega)f(\Delta\omega)\\\
&=\frac{2m^{2}v_{s}}{\hbar^{2}k_{F}\omega_{0}}\begin{cases}\omega_{0}\int_{-\infty}^{\infty}d(\Delta\omega)f(\Delta\omega),&f(\Delta\omega)\textrm{
even},\\\
\int_{-\infty}^{\infty}d(\Delta\omega)(\Delta\omega)f(\Delta\omega),&f(\Delta\omega)\textrm{
odd}\end{cases}.\end{split}$ (48)

Note that we only need to consider the resonant terms for $f(\Delta\omega)$,
since the counter-resonant terms are negligible relative to the resonant terms
in the high-mobility limit. Intuitively, these results are explained by Fig.
10.

(a)

(b)

Figure 10: Phase-space diagram of near-resonance band spanning the electron-
phonon detuning range $\Delta\omega\approx-\gamma/2$ to
$\Delta\omega\approx\gamma/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).

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. 10(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. 10(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 $\sqrt{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 $\chi^{(1)}$ and $\chi^{(2)}$, discussed later in this section).

The second observation regards the scaling of the near-resonance band area
with the decay rate $\gamma$. For convenience, we define a new decay parameter
$\gamma^{\prime}=\gamma/2$. For each electron in the near-resonance band, the
average $N^{\textrm{th}}$-order susceptibility is proportional to
$1/\gamma^{\prime N}$. If a single-electron susceptibility term is even in the
detuning $\Delta\omega$, then the total susceptibility corresponding to that
term is simply determined by multiplying the average single-electron
susceptibility by the area of the near-resonance band (which is proportional
to $\gamma^{\prime}$, as shown in Fig. 10(b)). Therefore, the overall
susceptibility for the process represented by an even term is proportional to
$1/\gamma^{\prime N-1}$. However, if a single-electron susceptibility term is
odd in $\Delta\omega$, 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 positive-detuning cases.
As Fig. 10(b) shows, this difference is proportional to $\gamma^{\prime 2}$
rather than $\gamma^{\prime}$ when integrated over the band. As a result, the
overall susceptibility for the process represented by an odd term is
proportional to $1/\gamma^{\prime N-2}$.

It is worth addressing the special cases of $\textrm{Re}[\chi^{(1)}]$ and
$\textrm{Re}[\chi^{(2)}]$. In the case of $\textrm{Re}[\chi^{(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
$\textrm{Re}[\chi^{(1)}]$ are always approximately equal even in the high-
mobility limit. Consequently, $\textrm{Re}[\chi^{(1)}]$ is invariant in the
mobility. It is also important to note that since the contributions to
$\textrm{Re}[\chi^{(1)}]$ are far off-resonance, the initial electronic states
are primarily concentrated around $k_{x}=k_{F}$. Since $\Delta\omega\propto
k_{x}$ for $k_{x}\gg q$, and $\textrm{Re}[\chi^{(1)}]\propto 1/(\Delta\omega)$
for far-off-resonance 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
$\textrm{Re}[\chi^{(1)}]$ is invariant in the Fermi wavevector and thus
independent of carrier density.

The case of $\textrm{Re}[\chi^{(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
$\textrm{Re}[\chi^{(1)}]$. Consequently, in the high-mobility limit, the
electronic contributions to $\textrm{Re}[\chi^{(2)}]$ are neither dominantly
near-resonance or far-off-resonance, but rather in between the two limits. We
will thus tackle the problem of $\textrm{Re}[\chi^{(2)}]$ in the high-mobility
limit solely through numerical means.

###  IV.2 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 $\gamma$.
Following a derivation analogous to the high-mobility regime, the following
conditions on $\gamma$ and the mobility $\mu$ hold for the low-mobility
regime:

$\displaystyle\gamma$ $\displaystyle\gg\frac{\hbar qk_{F}}{m},$ (49)
$\displaystyle\mu$ $\displaystyle\ll\frac{q_{e}}{\hbar qk_{F}}.$ (50)

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 $\Delta\omega$ as follows:

$\textrm{Re}[f^{(N)}(\Delta\omega)]=\sum_{m=0}^{\frac{1}{2}(N-\frac{1}{2}+(-1)^{N}\frac{1}{2})}\frac{c_{N,m}\gamma^{\prime
2m}(\Delta\omega)^{N-2m}}{\prod_{n=1}^{N}a_{n}^{2}(\Delta\omega)^{2}+\gamma^{\prime
2}},$ (51)

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 $|\Delta\omega|\ll\gamma^{\prime}$. As
such, we can reduce each summation to the leading term in $\gamma^{\prime}$.
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):

$\textrm{Re}[f^{(N)}(\Delta\omega)]\approx\begin{cases}\frac{c_{N,N/2}}{\gamma^{\prime
N}},&N\textrm{ even},\\\ \frac{c_{N,(N-1)/2}\Delta\omega}{\gamma^{\prime
N+1}},&N\textrm{ odd}\end{cases}.$ (52)

Consequently, in the low-mobility limit, the real part of $\chi^{(N)}$ varies
with the mobility $\mu$ as $\mu^{N}$ ($\mu^{N+1}$) for even (odd) $N$,
assuming that the leading term does not cancel out. To determine the overall
$\chi^{(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 $q\ll k_{F}$ but not $k_{x}\ll k_{F}$:

$\displaystyle\begin{split}k_{y,\textrm{span}}(k_{x})&=2\Big{(}\sqrt{k_{F}^{2}-k_{x}^{2}}-\sqrt{k_{F}^{2}-(k_{x}+q)^{2}}\Big{)}\\\
&\approx
2\sqrt{k_{F}^{2}-k_{x}^{2}}-2\sqrt{k_{F}^{2}-k_{x}^{2}}\bigg{(}1-\frac{qk_{x}}{k_{F}^{2}-k_{x}^{2}}\bigg{)}\\\
&=\frac{2qk_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}.\end{split}$ (53)

The available phase-space area, which we label $A_{\mathrm{phase}}$, is solved
by simply integrating the $k_{y}$-span over $k_{x}$ as follows:

$\displaystyle\begin{split}A_{\textrm{phase}}=\int_{0}^{k_{F}}dk_{x}k_{y,\textrm{span}}(k_{x})&\approx\int_{0}^{k_{F}}dk_{x}\frac{2qk_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\\\
&=\frac{2\omega_{0}}{v_{s}}\int_{0}^{k_{F}}dk_{x}\frac{k_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\\\
&=\frac{2\omega_{0}k_{F}}{v_{s}}\int_{0}^{\pi/2}d\theta\cos{\theta}\\\
&=\frac{2\omega_{0}k_{F}}{v_{s}}.\end{split}$ (54)

Since the single-electron processes feeding into $\textrm{Re}[\chi^{(N)}]$ for
even $N$ are invariant in the detuning, the overall susceptibility can simply
be calculated by multiplying $\textrm{Re}[f^{(N)}]$ by the phase-space area
$A_{\mathrm{phase}}$. To determine the overall $\textrm{Re}[\chi^{(N)}]$ for
odd values of $N$, however, we need to integrate the detuning $\Delta\omega$
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}\gg q$ and thus
$\Delta\omega\approx\omega_{k_{x}+q,k_{x}}\approx\hbar qk_{x}/m$. We thus
integrate $\Delta\omega$ over the valid initial phase-space as follows:

$\displaystyle\begin{split}\int_{-q/2}^{k_{F}}dk_{x}k_{y,\textrm{span}}(k_{x})\Delta\omega&\approx\int_{0}^{k_{F}}dk_{x}\frac{2qk_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\frac{\hbar
qk_{x}}{m}\\\
&=\frac{2\hbar\omega_{0}^{2}}{mv_{s}^{2}}\int_{0}^{k_{F}}dk_{x}\frac{k_{x}^{2}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\\\
&=\frac{2\hbar\omega_{0}^{2}k_{F}^{2}}{mv_{s}^{2}}\int_{0}^{\pi/2}d\theta\cos^{2}{\theta}\\\
&=\frac{\pi\hbar\omega_{0}^{2}k_{F}^{2}}{2mv_{s}^{2}}.\end{split}$ (55)

Note that the result as $k_{F}^{2}$. 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 $\left\langle\Delta\omega\right\rangle\propto k_{F}$. As such,
$\textrm{Re}[\chi^{(N)}]$ scales linearly with the carrier density $n$ for odd
$N$. On the other hand, for even $N$, $\textrm{Re}[\chi^{(N)}]$ scales as
$n^{1/2}$, since the electron-phonon interaction probability in this case is
invariant in the detuning $\Delta\omega$.

It is also worth examining the integral of $(\Delta\omega)^{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 $(\Delta\omega)^{n}$ where $n>1$:

$\displaystyle\begin{split}&\int_{-q/2}^{k_{F}}dk_{x}k_{y,\textrm{span}}(k_{x})(\Delta\omega)^{n}\\\
&\approx\int_{0}^{k_{F}}dk_{x}\frac{2qk_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\bigg{(}\frac{\hbar
qk_{x}}{m}\bigg{)}^{n}\\\
&=\frac{2\hbar^{n}\omega_{0}^{n+1}}{m^{n}v_{s}^{n+1}}\int_{0}^{k_{F}}dk_{x}\frac{k_{x}^{n+1}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\\\
&=\frac{2\hbar^{n}\omega_{0}^{n+1}k_{F}^{n+1}}{m^{n}v_{s}^{n+1}}\int_{0}^{\pi/2}d\theta\cos^{n+1}{\theta}\\\
&=\frac{\sqrt{\pi}\hbar^{n}\omega_{0}^{n+1}k_{F}^{n+1}}{m^{n}v_{s}^{n+1}}\frac{\Gamma\Big{(}\frac{n+2}{2}\Big{)}}{\Gamma\Big{(}\frac{n+3}{2}\Big{)}}.\end{split}$
(56)

This result varies as $k_{F}^{n+1}$, since the number of available electrons
scales linearly with $k_{F}$, while the interaction probability for each
electron scales as $k_{F}^{n}$.

Finally, we comment on the imaginary parts of $\chi^{(N)}$. Separating these
into individual terms corresponding to specific processes, the imaginary parts
of the individual susceptibility terms $f^{(N)}$ take the following form:

$\displaystyle\begin{split}&\textrm{Im}[f^{(N)}(\Delta\omega)]=\\\
&\quad\sum_{m=1}^{\frac{1}{2}(N+\frac{1}{2}-(-1)^{N}\frac{1}{2})}\frac{c_{N,m}\gamma^{\prime
2m-1}(\Delta\omega)^{N-2m+1}}{\prod_{n=1}^{N}a_{n}^{2}(\Delta\omega)^{2}+\gamma^{\prime
2}}.\end{split}$ (57)

Using the same argument as with the real parts, the leading term reduces to
the following in the low-mobility limit:

$\textrm{Im}[f^{(N)}(\Delta\omega)]\approx\begin{cases}\frac{c_{N,N/2}\Delta\omega}{\gamma^{\prime
N+1}},&N\textrm{ even},\\\ \frac{c_{N,(N+1)/2}}{\gamma^{\prime N}},&N\textrm{
odd}\end{cases}.$ (58)

The leading term for $\textrm{Im}[\chi^{(N)}]$ therefore scales with mobility
as $\mu^{N+1}$ ($\mu^{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 $\textrm{Im}[\chi^{(N)}]$
for a generic $N$ must be calculated by making a higher-order expansion of the
terms to account for the slight difference between the absorption and emission
probabilities.

###  IV.3 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 $\chi^{(1)}$, which governs the dielectric constant for
the 2DEG. As discussed in the previous sections, $\textrm{Re}[\chi^{(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. (41) into Eq.
(51) and simplifying per Eq. (52), we find that $\textrm{Re}[f_{k_{x}}^{(1)}]$
approximately reduces to the following:

$\textrm{Re}[f^{(1)}(\Delta\omega)]\approx\frac{2\Delta\omega}{\gamma^{\prime
2}}.$ (59)

Note that the real part of $\chi^{(1)}$ thus scales with mobility as $\mu^{2}$
in the low-mobility limit. Integrating $\Delta\omega$ over the phase-space
area of the initial electronic states as in Eq. (55) and multiplying by the
constants shown in Eq. (LABEL:eq:_chi(1)_integral), we find the following
analytical expression for the real part of $\chi^{(1)}$ in the low-mobility
limit:

$\displaystyle\begin{split}\textrm{Re}[\chi^{(1)}(\omega_{0})]&\approx\bigg{(}\frac{q_{e}^{2}v_{s}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}\bigg{)}\frac{2}{\gamma^{\prime
2}}\bigg{(}\frac{\pi\hbar\omega_{0}^{2}k_{F}^{2}}{2mv_{s}^{2}}\bigg{)}\\\
&=\frac{q_{e}^{2}k_{F}^{2}}{2\pi\epsilon_{0}t_{\mathrm{2DEG}}m\gamma^{\prime
2}}\\\
&=\frac{2mk_{F}^{2}\mu^{2}}{\pi\epsilon_{0}t_{\mathrm{2DEG}}},\end{split}$
(60)

where we substituted the relationship $\gamma^{\prime}=q_{e}/(2m\mu)$ in the
last line. Note that this susceptibility is independent of the phonon
frequency $\omega_{0}$ in the limit $q\ll 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 $\omega_{0}$,
thus cancelling out the inverse dependence of the dipole interaction
probability with $\omega_{0}^{2}$.

Next, we examine the real part of $\chi^{(1)}$ in the high-mobility limit. As
discussed in Sec. IV.1, $\textrm{Re}[\chi^{(1)}]$ is dominated by the far-off-
resonance interactions, unlike the general case for the high-mobility limit.
Using the approximation $|\Delta\omega|\gg\gamma^{\prime}$, we find that
$\textrm{Re}[f_{k_{x}}^{(1)}]$ approximately reduces to:

$\textrm{Re}[f^{(1)}(\Delta\omega)]\approx\frac{2}{\Delta\omega}.$ (61)

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}\gg q$ and thus $\Delta\omega\approx\omega_{k_{x}+q,k_{x}}\approx\hbar
qk_{x}/m$. Integrating $(\Delta\omega)^{-1}$ over the valid initial phase-
space area yields the following result:

$\displaystyle\begin{split}\int_{-q/2}^{k_{F}}dk_{x}k_{y,\textrm{span}}(k_{x})(\Delta\omega)^{-1}&\approx\int_{0}^{k_{F}}dk_{x}\frac{2qk_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\frac{m}{\hbar
qk_{x}}\\\
&=\frac{2m}{\hbar}\int_{0}^{k_{F}}\frac{dk_{x}}{\sqrt{k_{F}^{2}-k_{x}^{2}}}\\\
&=\frac{2m}{\hbar}\int_{0}^{\pi/2}d\theta\\\ &=\frac{\pi
m}{\hbar}.\end{split}$ (62)

Multiplying this by 2 (per Eq. (61)) and by constants shown in Eq.
(LABEL:eq:_chi(1)_integral), we find the following analytical expression for
the real part of $\chi^{(1)}$ in the high-mobility limit:

$\displaystyle\begin{split}\textrm{Re}[\chi^{(1)}(\omega_{0})]&\approx\bigg{(}\frac{q_{e}^{2}v_{s}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}\bigg{)}2\bigg{(}\frac{\pi
m}{\hbar}\bigg{)}\\\
&=\frac{q_{e}^{2}v_{s}^{2}m}{\pi\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}.\end{split}$
(63)

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 far-off-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 $\sqrt{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/\sqrt{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 $\sqrt{n}$).

We now turn to the imaginary part of $\chi^{(1)}$, which governs the net
phonon absorption rate and thus the spectral broadening for the phonon energy
levels. As discussed in Sec. IV.2, 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_{k_{x}}^{(1)}$
in Eq. (41), yielding:

$\displaystyle f^{(1)}(\Delta\omega)$
$\displaystyle\approx\frac{1}{\Delta\omega-i\gamma^{\prime}},$ (64)
$\displaystyle\textrm{Im}[f^{(1)}(\Delta\omega)]$
$\displaystyle\approx\frac{\gamma^{\prime}}{(\Delta\omega)^{2}+\gamma^{\prime
2}}.$ (65)

Using the near-resonance approximation and the corresponding procedure in Eq.
(LABEL:eq:_chi(Delta_omega)_phase_space_integral), we integrate this over the
near-resonance band (noting that the integrand is even in $\Delta\omega$) and
multiply by the coefficients in Eq. (LABEL:eq:_chi(1)_integral), leading to
the following result for the imaginary part of $\chi^{(1)}$ in the high-
mobility limit:

$\displaystyle\begin{split}&\textrm{Im}[\chi^{(1)}(\omega_{0})]\\\
&\approx\bigg{(}\frac{q_{e}^{2}v_{s}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}\bigg{)}\frac{2m^{2}v_{s}}{\hbar^{2}k_{F}}\int_{-\infty}^{\infty}d(\Delta\omega)\frac{\gamma^{\prime}}{(\Delta\omega)^{2}+\gamma^{\prime
2}}\\\
&=\bigg{(}\frac{q_{e}^{2}v_{s}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}\bigg{)}\frac{2m^{2}v_{s}}{\hbar^{2}k_{F}}\pi\\\
&=\frac{q_{e}^{2}v_{s}^{3}m^{2}}{\pi\hbar^{3}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}k_{F}}.\end{split}$
(66)

As such, the imaginary part of $\chi^{(1)}$ is invariant in $\gamma^{\prime}$
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
$\gamma^{\prime}$, while the interaction probability for each electron scales
inversely with $\gamma^{\prime}$, causing the overall interaction rate to be
constant in $\gamma^{\prime}$.

Finally, we solve for the imaginary part of $\chi^{(1)}$ in the low-mobility
limit. If we take the resonant and counter-resonant terms (corresponding to
phonon absorption and emission, respectively) in Eq. (41) individually, then
the magnitude of the imaginary part of each of those terms should
approximately reduce to $1/\gamma^{\prime}$, since
$\gamma^{\prime}\gg|\Delta\omega|\gg\omega_{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:

$\displaystyle\begin{split}&\textrm{Im}[f^{(1)}_{k_{x}}](\omega_{0})\\\
&=\frac{\gamma^{\prime}}{(\omega_{k_{x}+q,k_{x}}-\omega_{0})^{2}+\gamma^{\prime
2}}-\frac{\gamma^{\prime}}{(\omega_{k_{x}+q,k_{x}}+\omega_{0})^{2}+\gamma^{\prime
2}}\\\ &\approx\frac{4\omega_{0}\omega_{k_{x}+q,k_{x}}}{\gamma^{\prime 3}}\\\
&\approx\frac{4\omega_{0}\Delta\omega}{\gamma^{\prime 3}}.\end{split}$ (67)

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_{k_{x}}^{(1)}$, the imaginary part is approximately
linear in the electron-phonon detuning $\Delta\omega$. We therefore integrate
over the phase-space area of initial electronic states as in Eq. (55) and
mutiply by the coefficients in Eq. (LABEL:eq:_chi(1)_integral):

$\displaystyle\begin{split}\textrm{Im}[\chi^{(1)}(\omega_{0})]&\approx\bigg{(}\frac{q_{e}^{2}v_{s}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}\bigg{)}\frac{4\omega_{0}}{\gamma^{\prime
3}}\bigg{(}\frac{\pi\hbar\omega_{0}^{2}k_{F}^{2}}{2mv_{s}^{2}}\bigg{)}\\\
&=\frac{q_{e}^{2}\omega_{0}k_{F}^{2}}{\pi\epsilon_{0}mt_{\mathrm{2DEG}}\gamma^{\prime
3}}\\\
&=\frac{8m^{2}\omega_{0}k_{F}^{2}\mu^{3}}{\pi\epsilon_{0}q_{e}t_{\mathrm{2DEG}}}.\end{split}$
(68)

Although the overall susceptibilities corresponding to the absorptive and
emissive processes individually scale linearly with the mobility $\mu$, the
net process scales as $\mu^{3}$. As with the real part of $\chi^{(1)}$ in the
low-mobility regime, this result also varies linearly with carrier density
$n\propto k_{F}^{2}$, 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 11 depicts the numerically-calculated results for the real and
imaginary parts of $\chi^{(1)}$ (given a phonon angular frequency of
$\omega_{0}=2\pi\times 10^{9}\textrm{ s}^{-1}$, speed of sound $v_{s}=4\times
10^{3}$ m/s, a carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier
effective mass $m=0.067m_{0}$ relative to the bare electron mass $m_{0}$, and
a 2DEG thickness $t_{\mathrm{2DEG}}=2\times 10^{-8}\textrm{ m}$), along with
the analytical results in the high-mobility and low-mobility limits.

(a)

(b)

Figure 11: Results for the real (a) and imaginary (b) parts of
$\chi^{(1)}(\omega_{0})$ as functions of the carrier mobility $\mu$, including
numerical (solid, green), low-mobility analytical (dotted, red), and high-
mobility analytical (dashed, blue) calculations. We use a phonon angular
frequency of $\omega_{0}=2\pi\times 10^{9}\textrm{ s}^{-1}$, speed of sound
$v_{s}=4\times 10^{3}$ m/s, a carrier density $n=2\times 10^{15}\textrm{
m}^{-2}$, a carrier effective mass $m=0.067m_{0}$, and a 2DEG thickness
$t_{\mathrm{2DEG}}=2\times 10^{-8}\textrm{ m}$.

As desired, the numerical and analytical results match in the low-mobility and
high-mobility regimes, with the real (imaginary) part varying as $\mu^{2}$
($\mu^{3}$) in the low-mobility regime, and with both parts plateauing in the
high-mobility regime. Note that for all mobility values,
$\textrm{Re}[\chi^{(1)}]$ vastly exceeds $\textrm{Im}[\chi^{(1)}]$. The
dielectric constant can thus be approximated as fully real rather than
complex, which features important implications when calculating the wave-
mixing dynamics.

###  IV.4 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 $\chi^{(2)}(\omega_{1},\omega_{2})$. In the low-
mobility limit, we can reduce $\textrm{Re}[f_{k_{x}}^{(2)}]$ in Eq. (37) to
the following using the approximation
$|\Delta\omega|\gg\omega_{k_{x}+q_{1}+q_{2},k_{x}},\omega_{1},\omega_{2}$:

$\textrm{Re}[f^{(2)}_{\mathrm{SF}}]\approx-\frac{4}{\gamma^{\prime 2}}.$ (69)

The sum-frequency-generation susceptibility thus varies with mobility as
$\mu^{2}$ in this regime, following the general low-mobility rule laid out in
Eq. (52). We can analytically solve for the overall susceptibility in the low-
mobility regime by multiplying by the available phase-space area in Eq. (54)
and the coefficients in Eq. (LABEL:eq:_chi(2)_integral). It is important to
consider, however, that for two of the terms in Eq. (37), the width of the
phase-space area is $q_{1}=\omega_{1}/v_{s}$ (see Figs. 8(a) and (d)), while
for the other two terms, the width is $q_{2}=\omega_{2}/v_{s}$ (see Figs. 8(b)
and (c)). Consequently, for the available phase-space area expression in Eq.
(54), we apply the replacement
$\omega_{0}\rightarrow(\omega_{1}+\omega_{2})/2$, yielding the following
overall second-order susceptibility in the low-mobility limit:

$\displaystyle\begin{split}\textrm{Re}[\chi^{(2)}(\omega_{1},\omega_{2})]&\approx\bigg{(}\frac{q_{e}^{3}v_{s}^{3}}{2\pi^{2}\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{1}\omega_{2}(\omega_{1}+\omega_{2})}\bigg{)}\\\
&\quad\times\bigg{(}-\frac{4}{\gamma^{\prime
2}}\bigg{)}\bigg{(}\frac{2(\omega_{1}+\omega_{2})k_{F}}{2v_{s}}\bigg{)}\\\
&=-\frac{2q_{e}^{3}v_{s}^{2}k_{F}}{\pi^{2}\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{1}\omega_{2}\gamma^{\prime
2}}\\\
&=-\frac{8q_{e}v_{s}^{2}m^{2}k_{F}\mu^{2}}{\pi^{2}\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{1}\omega_{2}}.\end{split}$
(70)

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 12 depicts the numerically-calculated results for the amplitude of the
real part of $\chi^{(2)}(\omega_{1},\omega_{2})$ in the degenerate limit
$\omega_{1}=\omega_{0}=\omega_{2}$ representing the second-harmonic-generation
process (we label this quantity $\chi^{(2)}_{\mathrm{SHG}}(\omega_{0})$),
along with the analytical result in the low-mobility limit, given the same
parameters as in the $\chi^{(1)}$ calculation.

Figure 12: Results for the amplitude of the real part of
$\chi^{(2)}(\omega_{0},\omega_{0})$ as a function of the carrier mobility
$\mu$, including numerical (solid, green) and low-mobility analytical (dotted,
red) calculations. We use a phonon angular frequency of $\omega_{0}=2\pi\times
10^{9}\textrm{ s}^{-1}$, speed of sound $v_{s}=4\times 10^{3}$ m/s, a carrier
density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier effective mass
$m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

As desired, the numerical and analytical results match in the low-mobility
regime, with the susceptibility varying as $\mu^{2}$. It is interesting to
note that the susceptibility plateaus in the high-mobility regime. As
discussed in Sec. IV.1, the near-resonance electronic contributions to the
real part of $\chi^{(2)}$ cancel out in the high-mobility regime, leaving
behind electronic contributions that are somewhat off-resonance (but not far-
off-resonance, since $\textrm{Re}[\chi^{(2)}]$ converges to zero for large
detuning values). Along with the very low value of $\gamma$, this causes the
effect of the decay rate $\gamma$ on the electron interaction probabilities to
be negligible, leading $\textrm{Re}[\chi^{(2)}]$ to be constant in mobility in
the high-mobility limit.

Next, we turn to parametric amplification, focusing specifically on the
degenerate case ($\omega_{1}=\omega_{0}=\omega_{2}$). Examining
$\textrm{Re}[f_{k_{x}}^{(2)}]$ from Eq. (38) in the low-mobility limit, we
note that the amplitude of each of the terms can be approximated as
$1/\gamma^{\prime 2}$ to first order, as with second-harmonic 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(\Delta\omega)^{2}/\gamma^{\prime 4}$ (where
$\Delta\omega=\omega_{k_{x}+q_{0},k_{x}}-\omega_{0}$). As a result, the
parametric amplification susceptibility varies with mobility as $\mu^{4}$
instead of $\mu^{2}$:

$\textrm{Re}[f^{(2)}_{\mathrm{PA}}(\Delta\omega)]\approx\frac{8(\Delta\omega)^{2}}{\gamma^{\prime
4}}.$ (71)

To analytically solve for the low-mobility susceptibility, we need to
integrate $(\Delta\omega)^{2}$ over the available phase-space area. Since the
electron-phonon transitions in the low-mobility limit are concentrated far-
off-resonance, we approximate $\Delta\omega\approx\hbar qk_{x}/m$, enabling us
to apply the generic higher-order integral shown in Eq.
(LABEL:eq:_generic_higher-order_integral_low-mobility) as follows:

$\int_{-q/2}^{k_{F}}dk_{x}k_{y,\mathrm{span}}(k_{x})(\Delta\omega)^{2}\approx\frac{4\hbar^{2}\omega_{0}^{3}k_{F}^{3}}{3m^{2}v_{s}^{3}}.$
(72)

Multiplying this $f^{(2)}_{\mathrm{PA}}$ and the coefficients in Eq.
(LABEL:eq:_chi(2)_integral), we find the following degenerate parametric
amplification susceptibility in the low-mobility limit:

$\displaystyle\begin{split}\textrm{Re}[\chi^{(2)}(2\omega_{0},-\omega_{0})]&\approx\bigg{(}\frac{q_{e}^{3}v_{s}^{3}}{2\pi^{2}\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}2\omega_{0}^{3}}\bigg{)}\\\
&\quad\times\bigg{(}\frac{8}{\gamma^{\prime
4}}\bigg{)}\bigg{(}\frac{4\hbar^{2}\omega_{0}^{3}k_{F}^{3}}{3m^{2}v_{s}^{3}}\bigg{)}\\\
&=\frac{8q_{e}^{3}k_{F}^{3}}{3\pi^{2}\epsilon_{0}t_{\mathrm{2DEG}}m^{2}\gamma^{\prime
4}}\\\
&=\frac{128m^{2}k_{F}^{3}\mu^{4}}{3\pi^{2}\epsilon_{0}q_{e}t_{\mathrm{2DEG}}}.\end{split}$
(73)

Note that this is independent of the frequency $\omega_{0}$. This is because
the number of available electrons varies linearly with $\omega_{0}$, while the
interaction probability for each electron varies inversely with $\omega_{0}$,
and the two variations cancel out. Furthermore, since the interaction
probability per electron now varies as $k_{F}^{2}$, the susceptibility scales
as $k_{F}^{3}$, implying a variation with carrier density as $n^{3/2}$.

Figure 13 depicts the numerically-calculated results for the amplitude of the
real part of $\chi^{(2)}(2\omega_{0},-\omega_{0})$ (we label this quantity
$\chi^{(2)}_{\mathrm{PA}}(\omega_{0})$).

Figure 13: Results for the amplitude of the real part of
$\chi^{(2)}(2\omega_{0},-\omega_{0})$ as a function of the carrier mobility
$\mu$, including numerical (solid, green) and low-mobility analytical (dotted,
yellow) calculations. We use a phonon angular frequency of
$\omega_{0}=2\pi\times 10^{9}\textrm{ s}^{-1}$, speed of sound $v_{s}=4\times
10^{3}$ m/s, a carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier
effective mass $m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

As desired, the numerical and analytical results match in the low-mobility
regime, with the susceptibility varying as $\mu^{4}$. The susceptibility
levels out in the high-mobility regime for the same reason as in the second-
harmonic-generation case.

###  IV.5 Third-Order Kerr Susceptibility

We conclude our susceptibility calculations by deriving the third-order Kerr
susceptibility $\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})$ governing the
Kerr shift and the degenerate four-wave-mixing process. In the low-mobility
limit, we simplify $\textrm{Re}[f_{k_{x}}^{(3)}]$ in Eq. (34) to the
following, since $|\Delta\omega|\gg\omega_{k_{x}+2q,k_{x}},\omega$:

$\textrm{Re}[f^{(3)}(\Delta\omega)]\approx(-4+2+2-4)\frac{\gamma^{\prime
2}\Delta\omega}{\gamma^{\prime 6}}=-\frac{4\Delta\omega}{\gamma^{\prime 4}},$
(74)

where $\Delta\omega=\omega_{k_{x}+q,k_{x}}-\omega_{0}$. Following the general
rule for the low-mobility regime as outlined in Eq. (52), the Kerr
susceptibility varies with mobility as $\mu^{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. (55) and multiplying by the coefficients in Eq.
(LABEL:eq:_chi(3)_integral):

$\displaystyle\begin{split}&\textrm{Re}[\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})]\\\
&\approx\bigg{(}\frac{q_{e}^{4}v_{s}^{4}}{\pi^{2}\hbar^{3}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{4}}\bigg{)}\bigg{(}-\frac{4}{\gamma^{\prime
4}}\bigg{)}\bigg{(}\frac{\pi\hbar\omega_{0}^{2}k_{F}^{2}}{2mv_{s}^{2}}\bigg{)}\\\
&=-\frac{2q_{e}^{4}v_{s}^{2}k_{F}^{2}}{\pi\hbar^{2}\epsilon_{0}mt_{\mathrm{2DEG}}\omega_{0}^{2}\gamma^{\prime
4}}\\\
&=-\frac{32v_{s}^{2}m^{3}k_{F}^{2}\mu^{4}}{\pi\hbar^{2}\epsilon_{0}mt_{\mathrm{2DEG}}\omega_{0}^{2}}.\end{split}$
(75)

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 $\omega_{k_{x}+2q,k_{x}+q}$ for the electron
upon absorbing a second phonon of wavevector $q$ exceeds that upon absorbing
the first phonon $\omega_{k_{x}+q,k_{x}}$ by a constant ($k_{x}$-independent)
offset frequency, which we label $\omega^{\prime}$:

$\displaystyle\begin{split}\omega^{\prime}&=\omega_{k_{x}+2q,k_{x}+q}-\omega_{k_{x}+q,k_{x}}\\\
&=\frac{\hbar q(2(k_{x}+q)+q)}{2m}-\frac{\hbar q(2k_{x}+q)}{2m}\\\
&=\frac{\hbar q^{2}}{m}.\end{split}$ (76)

In terms of this offset, $\textrm{Re}[f_{k_{x}}^{(3)}]$ in Eq. (34) takes the
following form, considering only the resonant (first and second) terms:

$\displaystyle\begin{split}&\textrm{Re}[f^{(3)}(\Delta\omega)]\\\
&\approx\sum_{\pm}\frac{1}{(\Delta\omega\mp
i\gamma^{\prime})(2\Delta\omega+\omega^{\prime}-i\gamma^{\prime})(\Delta\omega-i\gamma^{\prime})}\\\
&=\frac{(2\Delta\omega+\omega^{\prime})(\Delta\omega)^{2}-\gamma^{\prime
2}(\Delta\omega\pm
3\Delta\omega\pm\omega^{\prime})}{((\Delta\omega)^{2}+\gamma^{\prime
2})^{2}((2\Delta\omega+\omega^{\prime})^{2}+\gamma^{\prime 2})},\end{split}$
(77)

where the top and bottom indices in the $\pm$ and $\mp$ operations represent
the first and second terms of Eq. (34), respectively. We integrate this over
the near-resonance band (using the procedure in Eq.
(LABEL:eq:_chi(Delta_omega)_phase_space_integral)) and multiply by the
coefficients in Eq. (LABEL:eq:_chi(3)_integral), leading to the following
result for the overall Kerr susceptibility in the high-mobility limit:

$\displaystyle\begin{split}\textrm{Re}[\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})]&\approx\bigg{(}\frac{q_{e}^{4}v_{s}^{4}}{\pi^{2}\hbar^{3}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{4}}\bigg{)}\frac{2m^{2}v_{s}}{\hbar^{2}k_{F}}\int_{-\infty}^{\infty}d(\Delta\omega)(\omega_{0}+\Delta\omega)\frac{(2\Delta\omega+\omega^{\prime})(\Delta\omega)^{2}-\gamma^{\prime
2}(\Delta\omega\pm
3\Delta\omega\pm\omega^{\prime})}{((\Delta\omega)^{2}+\gamma^{\prime
2})^{2}((2\Delta\omega+\omega^{\prime})^{2}+\gamma^{\prime 2})}\\\
&=\beta\frac{\pi(3\gamma^{\prime
2}+\omega^{\prime}\omega_{0})}{\gamma^{\prime}(9\gamma^{\prime
2}+\omega^{\prime 2})}\\\ &=\frac{2\beta
m\mu}{q_{e}}\bigg{(}\frac{3q_{e}^{2}+4\omega^{\prime}\omega_{0}m^{2}\mu^{2}}{9q_{e}^{2}+4\omega^{\prime
2}m^{2}\mu^{2}}\bigg{)},\end{split}$ (78)

where $\beta$ is defined as follows:

$\beta=\frac{2q_{e}^{4}v_{s}^{5}m^{2}}{\pi^{2}\hbar^{5}\epsilon_{0}\omega_{0}^{5}k_{F}t_{\mathrm{2DEG}}}.$
(79)

Note that the Kerr susceptibility varies linearly with $\mu$ 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 $\textrm{Re}[f^{(3)}(\Delta\omega)]$ odd in
$\Delta\omega$ and yielding
$\left\langle\textrm{Re}[f^{(3)}(\Delta\omega)\right\rangle\propto
1/\gamma^{\prime 3}$. As such, integrating the product of this function and
$\Delta\omega$ over the near-resonance band results in
$\textrm{Re}[\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})]\propto
1/\gamma^{\prime}\propto\mu$. 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 $\textrm{Re}[f^{(3)}(\Delta\omega)]$
even in $\Delta\omega$ and yielding
$\left\langle\textrm{Re}[f^{(3)}(\Delta\omega)\right\rangle\propto
1/\gamma^{\prime 2}$. Integrating this over the near-resonance band also
results in $\textrm{Re}[\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})]\propto
1/\gamma^{\prime}\propto\mu$, but with a different slope than the medium-high-
mobility case.

Figure 14 depicts the numerically-calculated results for the amplitude of the
real part of $\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})$, which we label
$\chi^{(3)}_{\mathrm{Kerr}}(\omega_{0})$.

Figure 14: Results for the amplitude of the real part of
$\chi^{(3)}(\omega_{0},\omega_{0},-\omega_{0})$ as a function of the carrier
mobility $\mu$, including numerical (solid, green), low-mobility analytical
(dotted, red), and high-mobility analytical (dashed, blue) calculations. We
use a phonon angular frequency of $\omega_{0}=2\pi\times 10^{9}\textrm{
s}^{-1}$, speed of sound $v_{s}=4\times 10^{3}$ m/s, a carrier density
$n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier effective mass $m=0.067m_{0}$,
and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times 10^{-8}\textrm{ m}$.

We use the same parameters as in the $\chi^{(1)}$ and $\chi^{(2)}$
calculations. As desired, the numerical and analytical results match in the
low-mobility and high-mobility regimes, varying as $\mu^{4}$ in the low-
mobility regime and as $\mu^{1}$ in the high-mobility regime. The
aforementioned change in the proportionality coefficient with $\mu$ between
the medium-high-mobility and ultra-high-mobility regimes is apparent in the
intercept shift in the log-log plot around $10^{3}\textrm{
m}^{2}/(\textrm{V}\cdot\textrm{s})$.

##  V Calculation of Mixing Dynamics

The $N^{\textrm{th}}$-order susceptibility manifests itself in the
$(N+1)$-wave-mixing process. For a quantitative understanding, we can
substitute the expression for $\left\langle
d_{\mathrm{eff}}(q_{p})\right\rangle$ in Eq.
(LABEL:eq:_dipole_moment_expectation_value) into the interaction Hamiltonian
in Eq. (LABEL:eq:_interaction_Hamiltonian_multiple_modes), yielding an
expansion in product of field amplitudes:

$\displaystyle\begin{split}H_{\mathrm{int}}&=-\epsilon_{0}V_{\mathrm{2DEG}}\sum_{n}\sum_{p}\sum_{\begin{subarray}{c}p_{1},...,p_{N}\\\
\omega_{p_{1}}+...+\omega_{p_{N}}=\omega\end{subarray}}\\\
&\quad\quad\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{n}})E^{*}(\omega)E(\omega_{p_{1}})...E(\omega_{p_{N}}).\end{split}$
(80)

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
$\omega=\omega_{p_{1}}+...+\omega_{p_{N}}$ can be expressed in operator form
as follows:

$H=\hbar\Big{(}ga^{{\dagger}}a_{1}...a_{N}+g^{*}a_{1}^{{\dagger}}...a_{N}^{{\dagger}}a\Big{)},$
(81)

where, as evidenced by Eq. (80), the $(N+1)$-phonon interaction rate $g$ is
given by the following:

$\displaystyle\begin{split}g&=-\frac{\epsilon_{0}V_{\mathrm{2DEG}}}{\hbar}\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{n}})\\\
&\quad\times
E_{\mathrm{zpf}}^{*}(\omega)E_{\mathrm{zpf}}(\omega_{p_{1}})...E_{\mathrm{zpf}}(\omega_{p_{N}}),\end{split}$
(82)

where $E_{\mathrm{zpf}}$ represents the zero-point electric field experienced
by the 2DEG electrons. Intuitively, since $\chi^{(N)}$ represents the total
induced electron polarization per unit amplitude of the $N$ polarizing fields
$p_{1},...,p_{N}$, multiplying it by
$V_{\mathrm{2DEG}}E_{\mathrm{zpf}}(\omega_{p_{1}})...E_{\mathrm{zpf}}(\omega_{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 $\omega$) 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. (82).

In general, we can determine the time-evolution of the output electric field
$\omega$ in terms of the amplitudes of the input electric fields using the
Heisenberg equation of motion as follows:

$\displaystyle\begin{split}\dot{E}(\omega)&=E_{\mathrm{zpf}}(\omega)\dot{a}\\\
&=-\frac{i}{\hbar}E_{\mathrm{zpf}}(\omega)[a,H]\\\
&=-iE_{\mathrm{zpf}}(\omega)ga_{1}...a_{N}\\\
&=i\frac{\epsilon_{0}V_{\mathrm{2DEG}}}{\hbar}\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{n}})\Big{|}E_{\mathrm{zpf}}(\omega)\Big{|}^{2}\\\
&\quad\times E(\omega_{1})...E(\omega_{N}).\end{split}$ (83)

The zero-point electric field $E_{\mathrm{zpf}}$ is determined by substituting
$A_{\mathrm{ph}}=\sqrt{v_{s}/(2V)}$ (corresponding to a vacuum fluctuation
spanning the mode volume) into Eq. (11), yielding:

$E_{\mathrm{zpf}}(\omega)=-\frac{C}{\epsilon(\omega)}\sqrt{\frac{\hbar\omega}{2V}},$
(84)

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_{\mathrm{ph}}(\omega)$ in terms of the input amplitudes
$A_{\mathrm{ph}}(\omega_{1}),...A_{\mathrm{ph}}(\omega_{N})$:

$\displaystyle\begin{split}&\partial_{x}A_{\mathrm{ph}}(\omega)=\\\ &\quad
i\frac{\epsilon_{0}}{\hbar}\frac{\chi^{(N)}(\omega_{p_{1}},...,\omega_{p_{n}})\sqrt{\omega\omega_{1}...\omega_{N}}}{\epsilon^{*}(\omega)\epsilon(\omega_{1})...\epsilon(\omega_{N})}\bigg{(}-C\sqrt{\frac{\hbar}{v_{s}}}\bigg{)}^{N+1}\\\
&\quad\quad\times\frac{V_{\mathrm{2DEG}}}{2V}A_{\mathrm{ph}}(\omega_{1})...A_{\mathrm{ph}}(\omega_{N}).\end{split}$
(85)

Note that the conversion efficiency is proportional to the ratio between the
2DEG volume $V_{\mathrm{2DEG}}$ and the overall mode volume $V$ (which simply
reduces to the ratio between the 2DEG thickness $t_{\mathrm{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
$\epsilon^{N+1}$, where $\epsilon$ is the effective 2DEG dielectric constant,
which is in turn proportional to $\textrm{Re}[\chi^{(1)}]$. Naively, we would
assume that $\epsilon=\epsilon_{0}\textrm{Re}[\chi^{(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 B, this attenuates the dielectric constant by a factor $a$, which is
a function of the ratio between the 2DEG thickness $t_{\mathrm{2DEG}}$ and the
field’s half-wavelength $l$:

$a=\frac{t_{\mathrm{2DEG}}}{2\pi
l}\bigg{(}2\log{\bigg{(}\frac{l}{t_{\mathrm{2DEG}}}\bigg{)}+3}\bigg{)}.$ (86)

The second caveat is that when $\textrm{Re}[\chi^{(1)}]$ is low enough that
the corresponding dielectric constant is below the dielectric constant of the
piezoelectric material $\epsilon_{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:

$\epsilon(\omega_{0})\approx\begin{cases}a\epsilon_{0}\textrm{Re}[\chi^{(1)}(\omega_{0})],&a\epsilon_{0}\textrm{Re}[\chi^{(1)}]>\epsilon_{p},\\\
\epsilon_{p},&a\epsilon_{0}\textrm{Re}[\chi^{(1)}(\omega_{0})]<\epsilon_{p}\end{cases}.$
(87)

The 2DEG dielectric constant is plotted in Fig. 15 given
$\epsilon_{p}=43.6\epsilon_{0}$, with all the other parameters the same as in
the susceptibility plots.

Figure 15: Dielectric constant in units of $\epsilon_{0}$ as a function of the
carrier mobility $\mu$, given a piezoelectric material dielectric constant
$\epsilon_{p}=43.6\epsilon_{0}$, a phonon angular frequency of
$\omega_{0}=2\pi\times 10^{9}\textrm{ s}^{-1}$, speed of sound $v_{s}=4\times
10^{3}$ m/s, a carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier
effective mass $m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

Given the dielectric constant and the susceptibilities, we solve for the rates
of 3 processes: linear absorption (corresponding to
$\textrm{Im}[\chi^{(1)}]$), as well as second-harmonic generation and
parametric amplification (corresponding to $\textrm{Re}[\chi^{(2)}]$). We also
derive the figure-of-merit for the Kerr nonlinearity by calculating the ratio
between the phonon Kerr shift (corresponding to $\textrm{Re}[\chi^{(3)}]$) and
the phonon spectral broadening (corresponding to $\textrm{Im}[\chi^{(1)}]$,
given a high-Q acoustic cavity, which has recently been achieved Shao  _et
al._ (2019)).

###  V.1 Linear Absorption

We start with the case of linear absorption. Based on Eq.
(LABEL:eq:_spatial_evolution_of_output_field), the spatial evolution of a
single field propagating through the heterostructure and experiencing linear
absorption takes the following form:

$\partial_{x}A_{\mathrm{ph}}(\omega_{0})=-\alpha A_{\mathrm{ph}}(\omega_{0}),$
(88)

where $r$ is the absorption rate, given as:

$\displaystyle\begin{split}\alpha&=-i\frac{\epsilon_{0}}{\hbar}\frac{(i\textrm{Im}[\chi^{(1)}(\omega)])\sqrt{\omega_{0}^{2}}}{\epsilon^{*}(\omega_{0})\epsilon(\omega_{0})}\bigg{(}-C\sqrt{\frac{\hbar}{v_{s}}}\bigg{)}^{2}\frac{V_{\mathrm{2DEG}}}{2V}\\\
&=\frac{\epsilon_{0}C^{2}V_{\mathrm{2DEG}}\omega_{0}}{2v_{s}V}\frac{\textrm{Im}[\chi^{(1)}(\omega_{0})]}{|\epsilon(\omega_{0})|^{2}}.\end{split}$
(89)

Figure 16 depicts the absorption rate given $C^{2}=6.1\epsilon_{0}$ and
effective mode depth $L_{z}=2\times 10^{-6}\textrm{ m}$ along with the same
parameters used in the susceptibility calculations, assuming that the 2DEG
covers the piezoelectric material such that
$V_{\mathrm{2DEG}}/V=t_{\mathrm{2DEG}}/L_{z}$.

Figure 16: Linear absorption per unit length $\alpha$ as a function of the
carrier mobility $\mu$, given a piezoelectric material dielectric constant
$\epsilon_{p}=43.6\epsilon_{0}$, $C^{2}=6.1\epsilon_{0}$, a phonon angular
frequency of $\omega_{0}=2\pi\times 10^{9}\textrm{ s}^{-1}$, mode depth
$L_{z}=2\times 10^{-6}$ m, speed of sound $v_{s}=4\times 10^{3}$ m/s, a
carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier effective mass
$m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

At very low mobility, the imaginary part of $\chi^{(1)}$ varies as $\mu^{3}$,
while the dielectric constant is flat in mobility, causing the absorption rate
to vary as $\mu^{3}$. For moderately low mobility, though, the dielectric-
constant-squared varies as $\mu^{4}$, causing the absorption rate to decrease
with mobility as $\mu^{-1}$. The linear absorption coefficient reaches a
constant value of about $25\textrm{ m}^{-1}$ at the high-mobility limit (since
both the dielectric constant and the imaginary part of $\chi^{(1)}$ are
constant at high mobility), corresponding to an effective phonon coherence
length of 40 mm.

###  V.2 Second-Harmonic Generation

Next, we solve for the second-harmonic-generation coupling strength. From Eq.
(LABEL:eq:_spatial_evolution_of_output_field), the output field spatially
evolves as follows:

$\partial_{x}A_{\mathrm{ph}}(2\omega_{0})=iGA_{\mathrm{ph}}(\omega_{1})A_{\mathrm{ph}}(\omega_{2}),$
(90)

where $2\omega_{0}=\omega_{1}+\omega_{2}$, and $G$ is the coupling strength,
given as:

$G=\frac{\epsilon_{0}}{\hbar}\frac{\textrm{Re}[\chi^{(2)}(\omega_{1},\omega_{2})]\sqrt{2\omega_{0}\omega_{1}\omega_{2}}}{\epsilon^{*}(2\omega_{0})\epsilon(\omega_{1})\epsilon(\omega_{2})}\bigg{(}-C\sqrt{\frac{\hbar}{v_{s}}}\bigg{)}^{3}\frac{V_{\mathrm{2DEG}}}{2V}.$
(91)

For degenerate second-harmonic generation (i.e,
$\omega_{1}=\omega_{0}=\omega_{2}$), this becomes the following:

$G=-\frac{\epsilon_{0}C^{3}V_{\mathrm{2DEG}}}{2V}\sqrt{\frac{2\hbar\omega_{0}^{3}}{v_{s}^{3}}}\frac{\textrm{Re}[\chi^{(2)}(\omega_{0},\omega_{0})]}{\epsilon^{*}(2\omega_{0})|\epsilon(\omega_{0})|^{2}}.$
(92)

Figure 17 depicts the second-harmonic coupling strength using the same
parameters as in the linear absorption case.

Figure 17: Amplitude of coupling strength $G$ for second-harmonic generation
as a function of the carrier mobility $\mu$, given a piezoelectric material
dielectric constant $\epsilon_{p}=43.6\epsilon_{0}$, $C^{2}=6.1\epsilon_{0}$,
a phonon angular frequency of $\omega_{0}=2\pi\times 10^{9}\textrm{ s}^{-1}$,
mode depth $L_{z}=2\times 10^{-6}$ m, speed of sound $v_{s}=4\times 10^{3}$
m/s, a carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier
effective mass $m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

At very low mobility, $\chi^{(2)}$ varies as $\mu^{2}$, while the dielectric
constant is flat in mobility, causing $|G|$ to vary as $\mu^{2}$. Once the
dielectric constant starts rising with mobility, however, the strong inverse
dependence of $|G|$ with the dielectric constant (specifically, the
$\epsilon^{-3}$ variation) causes $|G|$ to sharply drop with increasing
mobility despite the fact that $\chi^{(2)}$ increases with mobility in this
regime. However, at the high-mobility limit, both $\chi^{(2)}$ and the
dielectric constant reach constant values, causing $|G|$ to level out at
$7.6\times 10^{-11}\textrm{ s}^{1/2}$.

###  V.3 Parametric Amplification

Figure 18: 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).

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 Caves (1982); Yurke  _et al._ (1989);
Levandovsky  _et al._ (1999); Castellanos-Beltran  _et al._ (2008); Karlsson
(2016). 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 phase-preserving amplifier can never achieve the same
additive quantum noise as a phase-sensitive (traveling-wave) amplifier Caves
(1982); Clerk  _et al._ (2010).

A diagram of a traveling-wave parametric amplifier is shown in Fig. 18. Given
input fields with frequencies $\omega_{1}$ and
$2\omega_{0}=\omega_{1}+\omega_{2}$, the output field $\omega_{2}$ evolves in
a manner analogous to the output field in Sec. V.2:

$\partial_{x}A_{\mathrm{ph}}(\omega_{2})=iG_{p,1}A_{\mathrm{ph}}(2\omega_{0})A_{\mathrm{ph}}^{*}(\omega_{1}),$
(93)

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:

$\displaystyle\begin{split}&G_{p,1}=\\\
&\quad\frac{\epsilon_{0}}{\hbar}\frac{\textrm{Re}[\chi^{(2)}(2\omega_{0},-\omega_{1})]\sqrt{2\omega_{0}\omega_{1}\omega_{2}}}{\epsilon(2\omega_{0})\epsilon^{*}(\omega_{1})\epsilon^{*}(\omega_{2})}\bigg{(}-C\sqrt{\frac{\hbar}{v_{s}}}\bigg{)}^{3}\frac{V_{\mathrm{2DEG}}}{2V}.\end{split}$
(94)

We now specifically consider the case of a degenerate traveling-wave
parametric amplifier (i.e., $\omega_{1}=\omega_{0}=\omega_{2}$). If we set the
phase of the pump field $A_{\mathrm{ph}}(2\omega_{0})$ such that it is in
phase with $-iA_{\mathrm{ph}}(\omega_{0})/A_{\mathrm{ph}}^{*}(\omega_{0})$,
then the time-evolution of the signal/idler field amplitude takes the
following form:

$\partial_{x}|A_{\mathrm{ph}}(\omega_{0})|=G_{p,1}\Big{|}A_{\mathrm{ph}}(2\omega_{0})\Big{|}|A_{\mathrm{ph}}(\omega_{0})|.$
(95)

As such, the amplitude $|A_{\mathrm{ph}}(\omega_{0})|$ grows exponentially at
a rate $G_{p,1}\Big{|}A_{\mathrm{ph}}(2\omega_{0})\Big{|}$. The signal/idler
field power $P(\omega_{0})$ is proportional to the field intensity
$|A_{\mathrm{ph}}(\omega_{0})|^{2}$ and thus evolves at twice the rate of the
amplitude:

$P(\omega_{0},x)=P(\omega_{0},0)e^{2G_{p,1}|A_{\mathrm{ph}}(2\omega_{0})|x}.$
(96)

The gain per unit length per pump phonon current density amplitude is thus
expressed in $dB$ form as:

$G_{p,2}=\frac{2G_{p,1}}{\ln{(10)}}.$ (97)

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:

$\displaystyle\begin{split}G_{p}&=\frac{G_{p,2}}{\sqrt{\hbar
2\omega_{0}L_{y}L_{z}}}\\\
&=\frac{\epsilon_{0}C^{3}V_{\mathrm{2DEG}}}{\ln{(10)}V}\sqrt{\frac{\omega_{0}^{2}}{v_{s}^{3}L_{y}L_{z}}}\frac{\textrm{Re}[\chi^{(2)}(2\omega_{0},-\omega_{0})]}{\epsilon^{*}(2\omega_{0})|\epsilon(\omega_{0})|^{2}}.\end{split}$
(98)

Note that $G_{p}$ features an inverse dependence on the cross-sectional area
$L_{y}L_{z}$ in addition to the standard $V_{\mathrm{2DEG}}/V$ 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 19 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 $\omega_{0}=\pi\times 10^{10}\textrm{ s}^{-1}$), a mode
depth of a single wavelength (i.e., $L_{z}=800\textrm{ nm}$), and a mode width
of 15 wavelengths (i.e., $L_{y}=12\textrm{ }\mu\textrm{m}$).

Figure 19: Degenerate parametric amplification gain per unit length per unit
power amplitude $G_{p}$ as a function of the carrier mobility $\mu$, given a
piezoelectric material dielectric constant $\epsilon_{p}=43.6\epsilon_{0}$,
$C^{2}=6.1\epsilon_{0}$, a signal/idler angular frequency of
$\omega_{0}=\pi\times 10^{10}\textrm{ s}^{-1}$, mode depth $L_{z}=8\times
10^{-7}$ m, width $L_{y}=1.2\times 10^{-5}$ m, speed of sound $v_{s}=4\times
10^{3}$ m/s, a carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier
effective mass $m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

At very low mobility, $\chi^{(2)}(2\omega_{0},-\omega_{0})$ varies as
$\mu^{4}$, while the dielectric constant is invariant in mobility, causing
$G_{p}$ to vary as $\mu^{4}$, peaking at $0.58\textrm{
dB}/(\mu\textrm{m}\sqrt{\mu\textrm{W}})$, corresponding to a 20-dB gain in 35
$\mu$m (about 17.5% of the Rayleigh length) given a 1-$\mu$W 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 for the intermediate-mobility case, as with the
second-harmonic case. Eventually, at the high-mobility limit, both
$\chi^{(2)}(2\omega_{0},-\omega_{0})$ and the dielectric constant become
constant, causing $G_{p}$ to level out at $0.13\textrm{
dB}/(\mu\textrm{m}\sqrt{\mu\textrm{W}})$, corresponding to a 20-dB gain in 150
$\mu$m (about 3/4 of the Rayleigh length) given a 1-$\mu$W pump power input.

###  V.4 Kerr Nonlinearity

Figure 20: 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.

Finally, we solve for the Kerr nonlinearity figure-of-merit (F.O.M.), which we
define as the ratio between the Kerr shift 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. 20. From Eq. (82), the anharmonicity induced by
the Kerr shift in a phonon ladder of fundamental angular frequency
$\omega_{0}$ is calculated from the third-order Kerr susceptibility
$\chi^{(3)}_{\mathrm{Kerr}}$ as follows:

$\displaystyle\begin{split}|g_{\mathrm{Kerr}}|&=\frac{\epsilon_{0}V_{\mathrm{2DEG}}}{\hbar}\chi^{(3)}_{\mathrm{Kerr}}(\omega_{0})\Big{|}E_{\mathrm{zpf}}(\omega_{0})\Big{|}^{4}.\end{split}$
(99)

On the other hand, given a high-Q acoustic cavity (which has recently been
achieved Shao  _et al._ (2019)), the spectral broadening is given by the
linear absorption per unit time, which corresponds to
$\textrm{Im}[\chi^{(1)}]$. We convert the result for the linear absorption per
unit propagation distance $\alpha$ from Sec. V.1 into per-unit-time absorption
by multiplying by the speed of sound $v_{s}$, yielding:

$\gamma_{ph}=\frac{\epsilon_{0}V_{\mathrm{2DEG}}}{\hbar}\textrm{Im}[\chi^{(1)(\omega_{0})}]\Big{|}E_{\mathrm{zpf}}(\omega_{0})\Big{|}^{2},$
(100)

where $\gamma_{ph}$ is the phononic spectral broadening. The Kerr-nonlinearity
figure-of-merit (F.O.M.) is therefore solved as follows:

$\displaystyle\begin{split}\textrm{F.O.M.}&=\frac{|g_{\mathrm{Kerr}}|}{\gamma_{ph}}\\\
&=\frac{\chi^{(3)}_{\mathrm{Kerr}}(\omega_{0})|E_{\mathrm{zpf}}(\omega_{0})|^{2}}{\textrm{Im}[\chi^{(1)}(\omega_{0})]}\\\
&=\frac{C^{2}\hbar\omega_{0}\chi^{(3)}_{\mathrm{Kerr}}(\omega_{0})}{2V\epsilon^{2}(\omega)\textrm{Im}[\chi^{(1)}(\omega_{0})]}.\end{split}$
(101)

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 21 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 $\omega_{0}$ is changed to $\pi\times
10^{10}\textrm{ s}^{-1}$ to align with the 5-GHz transmon frequency.

Figure 21: Kerr nonlinearity figure-of-merit (F.O.M.) as a function of the
carrier mobility $\mu$, given a piezoelectric material dielectric constant
$\epsilon_{p}=43.6\epsilon_{0}$, $C^{2}=6.1\epsilon_{0}$, a phonon angular
frequency of $\omega_{0}=\pi\times 10^{10}\textrm{ s}^{-1}$, mode length of
$4\times 10^{-7}$ m along each dimension, speed of sound $v_{s}=4\times
10^{3}$ m/s, a carrier density $n=2\times 10^{15}\textrm{ m}^{-2}$, a carrier
effective mass $m=0.067m_{0}$, and a 2DEG thickness $t_{\mathrm{2DEG}}=2\times
10^{-8}\textrm{ m}$.

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 $\mu$. This is because
$\chi^{(3)}_{\mathrm{Kerr}}$ also scales linearly with $\mu$ in this regime,
while both the real and imaginary parts of $\chi^{(1)}$ are constant in
mobility. At very high mobilies (above $10^{4}\textrm{
m}^{2}/(\textrm{V}\cdot\textrm{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 $\textrm{Im}[\chi^{(1)}]$) and the sum-
frequency generation (corresponding to $\textrm{Re}[\chi^{(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
Wherrett (1983); Lehmberg  _et al._ (1976); Yajima and Souma (1978); Khoo and
Code (1980); Panock and Temkin (1977).

###  VI.1 Linear Absorption in Fluid Limit

We start by considering the linear absorption. In the low-mobility regime, the
resonant part of $\textrm{Im}[f_{k_{x}}^{(1)}]$ from Eq. (41) reduces to:

$\textrm{Im}[f^{(1)}_{\mathrm{res}}]=\frac{\gamma^{\prime}}{(\omega_{k_{x}+q,k_{x}}-\omega)^{2}+\gamma^{\prime
2}}\approx\frac{1}{\gamma^{\prime}},$ (102)

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 $\chi^{(1)}$ simply by multiplying
this expression by the Fermi circle area and the coefficients in Eq.
(LABEL:eq:_chi(1)_integral), yielding:

$\displaystyle\begin{split}\textrm{Im}[\chi^{(1)}(\omega_{0})]&\approx\bigg{(}\frac{q_{e}^{2}v_{s}^{2}}{2\pi^{2}\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}}\bigg{)}\frac{1}{\gamma^{\prime}}\Big{(}\pi
k_{F}^{2}\Big{)}\\\
&=\frac{q_{e}^{2}v_{s}^{2}n}{\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}\gamma^{\prime}}\\\
&=\frac{2q_{e}v_{s}^{2}mn\mu}{\hbar\epsilon_{0}t_{\mathrm{2DEG}}\omega_{0}^{2}},\end{split}$
(103)

where we substituted $k_{F}^{2}=2\pi n$ in the second line (where $n$ is the
area carrier density) and $\gamma^{\prime}=q_{e}/(2m\mu)$ in the third line.

Now, we apply another consequence of the fluid treatment of electrons: that
all ground-state electrons feature zero momentum ($\bm{k}=0$). Then, a
momentum kick from absorbing a phonon of wavevector $q$ would elevate a
ground-state electron’s energy to $\hbar q^{2}/(2m)$. If we assume that the
electron-phonon interaction is resonant, then $\hbar q^{2}/(2m)$ must equal
$v_{s}q$ (i.e., $\hbar\omega_{0}=2mv_{s}^{2}$). Substituting this into the
above expression for the imaginary part of $\chi^{(1)}$ yields the following
result:

$\textrm{Im}[\chi^{(1)}](\omega_{0})\approx\frac{q_{e}n_{\mathrm{bulk}}\mu}{\epsilon_{0}\omega_{0}}=\frac{\omega_{c}\epsilon}{\epsilon_{0}\omega_{0}},$
(104)

where we defined a bulk volumetric electron density
$n_{\mathrm{bulk}}=n/t_{\mathrm{2DEG}}$, along with the frequency parameter
$\omega_{c}=n_{\mathrm{bulk}}q_{e}\mu/\epsilon$. It is worth noting that the
imaginary part of $\chi^{(1)}$ varies linearly with both carrier density and
mobility via $\omega_{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 $\alpha$ from the
imaginary part of $\chi^{(1)}$ by multiplying by
$\epsilon_{0}V_{\mathrm{2DEG}}|E(\omega_{0})|^{2}/\hbar$ 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 $\epsilon_{\mathrm{eff}}$, the elasticity $\kappa$, and
the frequency $\omega_{0}$:

$|E_{\mathrm{zpf}}(\omega_{0})|^{2}=\frac{e^{2}}{|\epsilon_{\mathrm{eff}}|^{2}}|S_{\mathrm{zpf}}|^{2}=\frac{K^{2}\epsilon\hbar\omega_{0}}{2|\epsilon_{\mathrm{eff}}|^{2}V_{\mathrm{piezo}}},$
(105)

where $\epsilon$ is the dielectric screening (i.e., the real part of
$\chi^{(1)}$) and $K^{2}=e^{2}/(\epsilon\kappa)$ represents the square of the
electromechanical coupling. This accords with Eq. (84), except for the fact
that $\epsilon_{\mathrm{eff}}$ incorporates both the real and imaginary parts
of $\chi^{(1)}$ for the sake of generality:

$\epsilon_{\mathrm{eff}}=\epsilon_{0}\Big{(}\textrm{Re}[\chi^{(1)}(\omega_{0})]+i\textrm{Im}[\chi^{(1)}(\omega_{0})]\Big{)}=\epsilon\bigg{(}1+i\frac{\omega_{c}}{\omega_{0}}\bigg{)}.$
(106)

The absorption per unit length $\alpha$ thus becomes the following:

$\displaystyle\begin{split}\alpha&=\frac{\epsilon_{0}V_{\mathrm{2DEG}}}{\hbar
v_{s}}\textrm{Im}[\chi^{(1)}](\omega_{0})\frac{K^{2}\epsilon\hbar\omega_{0}}{2|\epsilon_{\mathrm{eff}}|^{2}V_{\mathrm{piezo}}}\\\
&=\frac{\epsilon_{0}V_{\mathrm{2DEG}}}{\hbar
v_{s}}\frac{\omega_{c}\epsilon}{\epsilon_{0}\omega_{0}}\frac{K^{2}\epsilon\hbar\omega_{0}}{2\epsilon^{2}(1+\omega_{c}^{2}/\omega_{0}^{2})V_{\mathrm{2DEG}}}\\\
&=\frac{K^{2}}{2}\frac{\omega_{c}}{v_{s}}\bigg{(}1+\frac{\omega_{c}^{2}}{\omega_{0}^{2}}\bigg{)}^{-1}.\end{split}$
(107)

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_{\mathrm{piezo}}=V_{\mathrm{2DEG}}$ for a bulk
piezoelectric semiconductor. Note that this result matches that derived
previously in the classical limit White (1962).

###  VI.2 Sum-Frequency Generation in Fluid Limit

We now turn to the question of resolving the real part of $\chi^{(2)}$ in for
the purpose of deriving the dynamics of sum-frequency generation. In the low-
mobility regime, the resonant part of $\textrm{Re}[f_{k_{x}}^{(2)}]$ from Eq.
(37) reduces to the following:

$\displaystyle\begin{split}&\textrm{Re}[f^{(2)}_{\mathrm{res}}]\\\
&=\bigg{(}\frac{1}{(\omega_{k_{x}+q_{1}+q_{2},k_{x}}-(\omega_{1}+\omega_{2})-i\frac{\gamma}{2})(\omega_{k_{x}+q_{1},k_{x}}-\omega_{1}-i\frac{\gamma}{2})}+\frac{1}{(\omega_{k_{x}+q_{1}+q_{2},k_{x}}-(\omega_{1}+\omega_{2})-i\frac{\gamma}{2})(\omega_{k_{x}+q_{2},k_{x}}-\omega_{2}-i\frac{\gamma}{2})}\bigg{)}\\\
&\approx-\frac{2}{\gamma^{\prime 2}},\end{split}$ (108)

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 $\chi^{(2)}$ by multiplying this expression by the Fermi circle
area and the coefficients in Eq. (LABEL:eq:_chi(2)_integral), yielding the
following result:

$\displaystyle\begin{split}\textrm{Re}[\chi^{(2)}(\omega_{1},\omega_{2})]&\approx\bigg{(}\frac{q_{e}^{3}v_{s}^{3}}{2\pi^{2}\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{1}\omega_{2}(\omega_{1}+\omega_{2})}\bigg{)}\\\
&\quad\times\bigg{(}-\frac{2}{\gamma^{\prime 2}}\bigg{)}\Big{(}\pi
k_{F}^{2}\Big{)}\\\
&=-\frac{2q_{e}^{3}v_{s}^{3}n}{\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{1}\omega_{2}(\omega_{1}+\omega_{2})\gamma^{\prime
2}}\\\
&=-\frac{8q_{e}v_{s}^{3}m^{2}n\mu^{2}}{\hbar^{2}\epsilon_{0}t_{\mathrm{2DEG}}\omega_{1}\omega_{2}\omega_{3}},\end{split}$
(109)

where $\omega_{3}=\omega_{1}+\omega_{2}$, and we substituted $k_{F}^{2}=2\pi
n$ and $\gamma^{\prime 2}=q_{e}^{2}/(4m^{2}\mu^{2})$ in the second and third
lines, respectively.

Applying the ground-state approximation described in the linear absorption
derivation, along with the single-phonon resonance condition that $\hbar
q^{2}/(2m)$ must equal $v_{s}q$ (i.e., $\hbar^{2}\omega^{2}=4m^{2}v_{s}^{4}$,
where $\omega_{1}\approx\omega\approx\omega_{2}$), we find that the real part
of $\chi^{(2)}$ reduces to the following:

$\textrm{Re}[\chi^{(2)}(\omega_{1},\omega_{2})]\approx-\frac{2q_{e}n_{\mathrm{bulk}}\mu^{2}}{\epsilon_{0}v_{s}\omega_{3}}=-\frac{2\epsilon\omega_{c}\mu}{\epsilon_{0}v_{s}\omega_{3}},$
(110)

where $n_{\mathrm{bulk}}$ and $\omega_{c}$ are defined as in the linear
absorption calculation, and $\omega_{3}=\omega_{1}+\omega_{2}$.

Next, we use the real part of $\chi^{(2)}$ to determine how the $\omega_{3}$
field evolves in time given input fields $\omega_{1}$ and $\omega_{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_{\mathrm{zpf},n}a_{n}$, where $a_{n}$ is the lowering
operator for the mode $n$ and $E_{\mathrm{zpf},n}$ is the zero-point electric
field amplitude. Recall that the electric field amplitude is directly
proportional to the strain field amplitude:

$E_{n}=-\frac{e}{\epsilon_{n}}S_{n}=-\frac{e}{\epsilon\Gamma_{n}}S_{n},$ (111)

where $\Gamma_{n}$ is defined as follows for the given mode $n$:

$\Gamma_{n}=1+i\frac{\omega_{c}}{\omega_{n}}.$ (112)

Also, recall that the zero-point electric field intensity varies linearly with
the phonon frequency and inversely with the mode volume as follows:

$|E_{\mathrm{zpf},n}|^{2}=\frac{e^{2}}{|\epsilon_{n}|^{2}}|S_{\mathrm{zpf},n}|^{2}=\frac{K^{2}\epsilon\hbar\omega_{n}}{2|\epsilon\Gamma_{n}|^{2}V}=\frac{K^{2}\hbar\omega_{n}}{2\epsilon|\Gamma_{n}|^{2}V}.$
(113)

The time-evolution of $E_{3}$ due to sum-frequency generation is thus solved
by applying the 3-wave-mixing Hamiltonian:

$\displaystyle\begin{split}\dot{E_{3}}^{(2)}&=-i\frac{E_{\mathrm{zpf},3}}{\hbar}[a_{3},H_{\mathrm{int}}^{(2)}]\\\
&=i\frac{\epsilon_{0}V}{\hbar}\textrm{Re}[\chi^{(2)}](\omega_{1},\omega_{2})\Big{|}E_{\mathrm{zpf},3}\Big{|}^{3}E_{2}E_{1}[a_{3},a_{3}^{{\dagger}}]\\\
&\approx
i\frac{\epsilon_{0}V}{\hbar}\bigg{(}-\frac{2\epsilon\omega_{c}\mu}{\epsilon_{0}v_{s}\omega_{3}}\bigg{)}\bigg{(}\frac{K^{2}\hbar\omega_{3}}{2\epsilon|\Gamma_{3}|^{2}V}\bigg{)}\bigg{(}-\frac{e}{\epsilon\Gamma_{2}}S_{2}\bigg{)}\bigg{(}-\frac{e}{\epsilon\Gamma_{1}}S_{1}\bigg{)}\\\
&=-iK^{2}\frac{e^{2}\mu\omega_{c}}{\epsilon^{2}v_{s}}\frac{1}{\Gamma_{1}\Gamma_{2}|\Gamma_{3}|^{2}}S_{2}S_{1}.\end{split}$
(114)

Finally, we solve for $E_{3}^{(2)}$ by dividing both sides by $-i\omega_{3}$
(since $\dot{E_{3}}^{(2)}=-i\omega_{3}E_{3}^{(2)}$):

$E_{3}^{(2)}=K^{2}\frac{e^{2}\mu}{\epsilon^{2}v_{s}}\frac{\omega_{c}}{\omega_{3}}\frac{1}{\Gamma_{1}\Gamma_{2}|\Gamma_{3}|^{2}}S_{2}S_{1}.$
(115)

Note that this result matches that calculated in the classical limit Conwell
and Ganguly (1971) if $K^{2}\rightarrow 1$ and $|\Gamma_{3}|^{2}\rightarrow
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