††thanks: These authors contributed equally.††thanks: These authors
contributed equally.††thanks: These authors contributed equally.

# Non-Hermitian dynamics and nonreciprocity of optically coupled nanoparticles

Manuel Reisenbauer  Vienna Center for Quantum Science and Technology (VCQ),
Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna,
Austria  Henning Rudolph  Faculty of Physics, University of Duisburg-Essen,
Lotharstraße 1, 47048 Duisburg, Germany  Livia Egyed  Vienna Center for
Quantum Science and Technology (VCQ), Faculty of Physics, University of
Vienna, Boltzmanngasse 5, A-1090 Vienna, Austria  Klaus Hornberger  Faculty of
Physics, University of Duisburg-Essen, Lotharstraße 1, 47048 Duisburg, Germany
Anton V. Zasedatelev  Vienna Center for Quantum Science and Technology (VCQ),
Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna,
Austria  Murad Abuzarli  Vienna Center for Quantum Science and Technology
(VCQ), Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090
Vienna, Austria  Benjamin A. Stickler  Institute for Complex Quantum Systems,
Ulm University, Albert-Einstein-Allee 11, D-89069 Ulm, Germany  Uroš Delić
<EMAIL_ADDRESS>Vienna Center for Quantum Science and Technology
(VCQ), Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090
Vienna, Austria

###### Abstract

Non-Hermitian dynamics, as observed in photonic, atomic, electrical, and
optomechanical platforms, holds great potential for sensing applications and
signal processing. Recently, fully tunable nonreciprocal optical interaction
has been demonstrated between levitated nanoparticles. Here, we use this
tunability to investigate the collective non-Hermitian dynamics of two
nonreciprocally and nonlinearly interacting nanoparticles. We observe parity-
time symmetry breaking and, for sufficiently strong coupling, a collective
mechanical lasing transition, where the particles move along stable limit
cycles. This work opens up a research avenue of nonequilibrium multi-particle
collective effects, tailored by the dynamic control of individual sites in a
tweezer array.

A plethora of physical phenomena are well described by Hermitian dynamics,
such as the dynamics of closed quantum systems, Landau-type phase transitions,
or transport along oscillator chains. However, the growth of complexity in
quantum many-body systems, the occurrence of chiral transport properties, and
the emergence of experiments that strongly couple to the environment require
models that go beyond Hermitian descriptions. A particularly interesting
example of this so-called non-Hermitian dynamics are nonreciprocal
interactions, which seemingly break Newton’s third law of action equals
reaction. Besides being often encountered in biological systems Fruchart  _et
al._ (2021), nonreciprocity – in its broadest sense – has found various
applications in optics and photonics Feng  _et al._ (2017); Xiao  _et al._
(2020); Weidemann  _et al._ (2020); Wang  _et al._ (2021); Guo  _et al._
(2009); Zhen  _et al._ (2015); Miri and Alù (2019), ultracold atoms Li  _et
al._ (2019); Gou  _et al._ (2020); Takasu  _et al._ (2020); Öztürk  _et al._
(2021); Ferri  _et al._ (2021); Liang  _et al._ (2022); Gaikwad  _et al._
(2023), electrical circuits Helbig  _et al._ (2020); Zou  _et al._ (2021), and
metamaterials Brandenbourger  _et al._ (2019); Ghatak  _et al._ (2020); Chen
_et al._ (2021); Wang  _et al._ (2022). The size and sensitivity to
environmental perturbations make nonreciprocally interacting arrays of
mechanical objects an ideal ground for realizing unidirectional or topological
transport Metelmann and Clerk (2015); Fang  _et al._ (2017); Sanavio  _et al._
(2020), enhanced sensing due to nonreciprocity Metelmann and Clerk (2014); Lau
and Clerk (2018); McDonald and Clerk (2020), and topological states Peano  _et
al._ (2015); Rosenthal  _et al._ (2018); Wanjura  _et al._ (2020). So far,
there have been several experimental demonstrations of nonreciprocal or non-
Hermitian dynamics in various platforms in the classical regime Mathew  _et
al._ (2020); Xu  _et al._ (2016); Patil  _et al._ (2022); Ren  _et al._
(2022); Zhang  _et al._ (2022); Youssefi  _et al._ (2022); Doster  _et al._
(2022); Liu  _et al._ (2023).

Optically levitated nanoparticles have become a well-established system for
quantum physics with translational and rotational degrees of freedom Gonzalez-
Ballestero  _et al._ (2023); Stickler  _et al._ (2021). Recently, there has
been a surge of experiments that extend trapping and control from single
particles to particle arrays in a variety of geometries Rieser  _et al._
(2022); Vijayan  _et al._ (2023); Yan  _et al._ (2023); Liška  _et al._
(2023); Vijayan  _et al._ (2023); Bykov  _et al._ (2023), thus demonstrating
that this platform is highly versatile and scalable. In one of those
experiments, we have demonstrated direct, nonreciprocal, and nonlinear light-
induced dipole-dipole interactions between particles in a tweezer array Rieser
_et al._ (2022). Such optically interacting particle arrays offer several
benefits for investigations of (quantum) non-Hermitian physics. For example,
single-site readout enables the full reconstruction of the collective degrees
of freedom Penny  _et al._ (2023); Liška  _et al._ (2023); Bykov  _et al._
(2023). At the same time, the optically induced forces allow for a wide tuning
range from reciprocal to unidirectional to anti-reciprocal interactions
Rudolph  _et al._ (2023). Altogether, this system enables studies in
previously unexplored interaction regimes with an unprecedented level of
control.

In this work, we report on the experimental investigation of non-Hermitian
dynamics stemming from the anti-reciprocal interaction between the motion of
two trapped silica nanoparticles. We observe two exceptional points (EP) that
define a region where the collective motion is in the parity and time-reversal
($\mathcal{PT}$) symmetry-broken phase. In this phase, the interaction leads
to correlated particle motion, which we confirm by measuring a constant phase
delay between the oscillators. The system further exhibits a Hopf bifurcation
into the mechanical lasing phase, where the interaction-induced amplification
dominates over the intrinsic damping such that the motion becomes nonlinear.

Figure 1: Experimental setup. a. Each of the two orthogonal acousto-optical
deflectors (AOD-x/y) is driven by the same two radio-frequency (RF) tones at
frequencies $\omega_{1}$ and $\omega_{1}+\Delta\omega$, which creates $2\times
2$ laser beams. We use a slit to select the two beams that have equal optical
frequencies, while the distance between the beams can be tuned by changing the
RF tone frequency difference $\Delta\omega$. The Dove prism rotates the
optical plane to place the beams into the plane of the optical table. The
beams are subsequently focused in the vacuum chamber to form two traps at a
distance $d_{0}$. The light back-scattered by the particles is reflected with
a Faraday rotator (FR) and a polarizing beamsplitter (PBS) and sent to two
independent heterodyne detectors that monitor the particle motion. Inset: Two
particles are trapped and interact anti-reciprocally with the coupling rate
$\pm g_{a}$ tuned by the polarization angle $\theta$, which we set with a
half-wave plate (HWP) in front of the vacuum chamber. b. For polarization
along the $x$-axis ($\theta\approx\pi/2$), the interaction between the
particles is weak such that the two modes cross (green and red points:
measured eigenfrequencies and damping rates, black lines: fits). c. For
$\theta\neq\pi/2$ we observe degenerate eigenfrequencies between the two
exceptional points (EP) (top) and nondegenerate damping rates that are split
by $4g_{a}$ (bottom). Black lines are theory functions based on the measured
coupling rate. d. At the maximum splitting of the damping rates we can
reconstruct the power spectral densities (PSD) of the eigenmodes of the
particles’ positions as $z_{1}\pm iz_{2}$ (bottom) from the detected positions
$z_{1,2}$ (top).

#### Experimental setup.

We use a pair of orthogonal acousto-optical deflectors (AODs), both driven by
two radio-frequency (RF) tones at frequencies $\omega_{1}$ and
$\omega_{1}+\Delta\omega$, to create $2\times 2$ laser beams from a laser
source at a wavelength of $\lambda=1064~{}\text{nm}$ (Fig. 1a). We use a
narrow slit (width: $2~{}\text{mm}$) to select the laser beams on the diagonal
as they have equal laser frequencies, which is required to control the optical
interaction. A Dove prism tilted at an angle of $22.5^{\circ}$ rotates the
optical plane by $45^{\circ}$, thus placing the laser beams in the plane of
the optical table. A high numerical aperture lens ($\text{NA}=0.77$) inside
the vacuum chamber focuses the laser beams to two foci at a relative distance
$d_{0}$, in which we trap two silica nanoparticles of approximately equal
sizes (nominal radius $r=(105\pm 2)~{}\text{nm}$ SI ). We tune the particle
distance with the frequency difference of the RF tones as
$d_{0}\propto\Delta\omega$ on both AODs, which maintains equal laser
frequencies of the two tweezers SI . The optical phases at the traps
$\phi_{1,2}$ are controlled by the phases of the RF tones. For optical powers
of $P_{1,2}\approx 0.3~{}\text{W}$ the resulting oscillation frequencies of
the center-of-mass (CoM) motion along the tweezer axis ($z$-axis) are
$\Omega_{1,2}\propto\sqrt{P_{1,2}}\approx 2\pi\times 27.5~{}\text{kHz}$. To
sweep the mechanical frequency detuning $\Delta\Omega=\Omega_{2}-\Omega_{1}$,
we scan powers $P_{1,2}$ symmetrically such that the total power is conserved.
We collect the light backscattered from the particles with separate fiber-
based confocal microscopes, which allows us to independently detect the
particles’ motion with balanced heterodyne detections SI . To suppress
electrostatic coupling we place the particles at a large distance of
$d_{0}\sim 18.4~{}\mu\text{m}$ and discharge them. The intrinsic damping rate,
given by the gas pressure, is kept constant at $\gamma/2\pi=(0.46\pm
0.02)~{}\text{kHz}$ throughout the measurement.

Figure 2: Collective statistics of uncoupled (top row), linear (middle) and
nonlinear (bottom) particles’ motion $\mathbf{z_{1}}$ (blue) and
$\mathbf{z_{2}}$ (orange). In the middle and bottom rows the system is in the
$\mathcal{PT}$ symmetry-broken phase. a. The histograms of the particles’
motion show an increasing variance (top to middle) and eventually the
transition from linear into nonlinear motion (bottom). b. Uncoupled particles
move independently, which is confirmed by the uniform distribution of the
phase delay $\bar{\psi}$ between the oscillators. On the other hand, in the
$\mathcal{PT}$ symmetry-broken phase the histograms show a preferred phase
$\bar{\psi}_{\text{max}}$ as $z_{1}$ and $z_{2}$ are strongly correlated. The
black lines mark the most probable phase delay $\bar{\psi}_{\text{max}}$. c.
The joint $z_{1}$-$z_{2}$ distributions show the transition from a thermal
motion (top) to a correlated motion (middle) to a limit cycle (bottom).

#### Non-Hermitian Hamiltonian.

We model the particles’ $z$-motion in the frame rotating with the mean
oscillation frequency $\Omega_{0}=(\Omega_{1}+\Omega_{2})/2$. The linearized
equations of motion, averaged over one oscillation period $2\pi/\Omega_{0}$,
are:

$\frac{d}{dt}\begin{pmatrix}a_{1}\\\ a_{2}\\\
\end{pmatrix}=-iH_{\text{NH}}\begin{pmatrix}a_{1}\\\ a_{2}\\\
\end{pmatrix}+\sqrt{\gamma n_{T}}\begin{pmatrix}\xi_{1}\\\ \xi_{2}\\\
\end{pmatrix},$ (1)

where
$a_{j}=\left(z_{j}+ip_{j}/m\Omega_{j}\right)\exp(i\Omega_{0}t)/2z_{\text{zpf},j}$
is the complex amplitude of particle $j$ with position $z_{j}$, momentum
$p_{j}$, zero point fluctuation $z_{\text{zpf},j}=\sqrt{\hbar/2m\Omega_{j}}$
and mass $m$. Gas damping and diffusion enters via the gas damping rate
$\gamma$, the thermal occupation $n_{T}$ in equilibrium with the environment
and the complex white noises $\xi_{1,2}$ with correlations
$\langle\xi^{*}_{j^{\prime}}(t^{\prime})\xi_{j}(t)\rangle=\delta_{jj^{\prime}}\delta(t-t^{\prime})$
and $\langle\xi_{j^{\prime}}(t^{\prime})\xi_{j}(t)\rangle=0$. The effective
non-Hermitian Hamiltonian realized by the optical interaction is Rudolph  _et
al._ (2023):

$H_{\text{NH}}=\begin{pmatrix}-\frac{\Delta\Omega}{2}+g_{r}+g_{a}-i\frac{\gamma}{2}&-(g_{r}+g_{a})\\\
-(g_{r}-g_{a})&\frac{\Delta\Omega}{2}+g_{r}-g_{a}-i\frac{\gamma}{2}\\\
\end{pmatrix},$ (2)

where $g_{r}=G\cos{(kd_{0})}\cos{(\Delta\phi)}/kd_{0}$ and
$g_{a}=G\sin{(kd_{0})}\sin{(\Delta\phi)}/kd_{0}$, with $k=2\pi/\lambda$,
describe the reciprocal (conservative) and the anti-reciprocal (non-
conservative) coupling rate, respectively Rieser  _et al._ (2022); Rudolph
_et al._ (2023). The nonreciprocity of the effective interaction arises from
the interference of the tweezer and the light scattered off the particles,
which effectively carries away momentum. The interference and thus the
coupling rates can be tuned by the optical phase difference
$\Delta\phi=\phi_{2}-\phi_{1}$ and the distance $d_{0}$ between the particles.
In our work, we set $d_{0}$ and $\Delta\phi$ such that $g_{r}$ is negligible,
while $g_{a}$ takes its maximum value for a given distance. Its magnitude is
then determined by the coupling constant
$G\propto\sqrt{P_{1}P_{2}}\cos^{2}(\theta)$ as a function of the trapping
powers $P_{1,2}$ and the laser polarization angle $\theta$ (Fig. 1b,c), which
we modify with a half-wave plate (HWP) in front of the vacuum chamber SI . The
eigenvalues of the effective Hamiltonian $H_{\text{NH}}$ are in general
complex and are given by
$\lambda_{\pm}=-i\gamma/2\pm\sqrt{\Delta\Omega^{2}-4g_{a}\Delta\Omega}/2$,
such that the frequencies and damping rates of the eigenmodes are given by the
real and imaginary parts
$\Omega_{\pm}=\Omega_{0}+\operatorname{Re}(\lambda_{\pm})$ and
$\gamma_{\pm}=-2\operatorname{Im}(\lambda_{\pm})$, respectively. The
eigenvectors coalesce at the detunings of
$\Delta\Omega_{\text{EP}1,2}=2g_{a}\mp 2g_{a}$, which define the exceptional
points (EP). For $\Delta\Omega$ between the EPs, the frequencies are
degenerate and the damping rates become non-degenerate, resulting in the so-
called ”normal mode attraction” (Fig. 1c). The maximum splitting of the
damping rates $\gamma_{\pm}$ is achieved for $\Delta\Omega=2g_{a}$ and is
equal to $4g_{a}$, where the complex eigenmodes of the system can be
reconstructed as $a_{\pm}=a_{1}\pm ia_{2}$. Note that the sign in the
corresponding eigenmodes of the motion $z_{\pm}=z_{1}\mp iz_{2}$ is flipped
due to the definition of $a_{1,2}$. Therefore, the different damping rates
result in the suppression of the eigenmode $a_{-}$ ($z_{-}$), while $a_{+}$
($z_{+}$) is amplified (Fig. 1d). Note that the Hamiltonian $H_{\text{NH}}$ is
in general $\mathcal{PT}$-symmetric; however, the $\mathcal{PT}$ symmetry is
broken in the region between the EPs as the eigenmodes have different damping
rates.

#### Nonlinear anti-reciprocal interactions.

Once the effective damping rate $\gamma_{-}$ becomes negative, the linear
theory breaks down. An analytical model that includes the full nonlinear
dynamics – but no thermal fluctuations – yields a system of differential
equations for the amplitude of the collective motional state
$A=2k\sqrt{\hbar/m\Omega_{0}}\sqrt{|a_{2}|^{2}+|a_{1}|^{2}}$ and the phase
delay between the oscillators $\psi=\arg[a_{2}^{*}a_{1}]$ SI :

$\displaystyle\begin{aligned}
\dot{A}&=&-\frac{\gamma}{2}A+g_{a}A\sin(\psi)f\left(A\sin\frac{\psi}{2}\right),\\\
\dot{\psi}&=&\Delta\Omega-4g_{a}\sin^{2}\left(\frac{\psi}{2}\right)f\left(A\sin\frac{\psi}{2}\right).\end{aligned}$
(3)

Here, $f(x)=2J_{1}(x)/x$ depends on the first-order Bessel function of the
first kind $J_{1}$. The equations of motion (3) are valid only in the
$\mathcal{PT}$ symmetry-broken phase and have two steady-state solutions: (i)
a collective state with vanishing oscillations ($A=0$) but a stable phase
delay $\psi=2\arcsin(\sqrt{\Delta\Omega/4g_{a}})$, and (ii) a coherently
oscillating state with $A$ satisfying
$f\left(A\Delta\Omega/\sqrt{\gamma^{2}+\Delta\Omega^{2}}\right)=(\gamma^{2}+\Delta\Omega^{2})/4g_{a}\Delta\Omega$,
and phase delay $\psi=2\arctan(\Delta\Omega/\gamma)$. As $f(x)\leq 1$, the two
solutions exist in regions separated by the threshold defined by
$g_{a}=\left(\gamma^{2}+\Delta\Omega^{2}\right)/4\Delta\Omega$. Above this
threshold the first solution becomes unstable and the second, truly nonlinear
solution – a stable limit cycle – emerges, thus revealing a Hopf bifurcation.

The dynamics of the collective motion under conditions defined by the two
solutions of Eq. (3) are shown in the middle and bottom row of Fig. 2,
respectively, while the top row features the standard behavior of uncoupled
thermal oscillators for comparison. In the linear regime, the individual
particles’ motion follows a Gaussian distribution, and we observe an increased
motional amplitude as we transition into the $\mathcal{PT}$ symmetry-broken
phase. For higher coupling rates, the particles’ motion become nonlinear,
which is reflected in a modified motional statistics to a displaced Gaussian
distribution (Fig. 2a). We compute the phases $\bar{\psi}_{1,2}$ of the
individual particles’ motion via the Hilbert transform and calculate the
instantaneous phase delay between the oscillators as
$\bar{\psi}=\bar{\psi}_{2}-\bar{\psi}_{1}$ SI . The stable phase delay
$\bar{\psi}$ is resolved in the zoomed-in time trace of the nonlinear motion
(Fig. 2b). The interaction generates a preferred phase delay
$\bar{\psi}_{\text{max}}$ as observed in the histograms of $\bar{\psi}$, while
it is uniformly distributed in the case of uncoupled oscillators. The
$z_{1}$-$z_{2}$ distribution shows the collective dynamics of the two
particles (Fig. 2c). In the case of uncoupled particles, the distribution is
well-described by an uncorrelated two-dimensional Gaussian distribution.
However, in the case of solution (i), $z_{1}$ and $z_{2}$ are strongly
correlated, and thus the joint distribution is squashed under an angle that
depends on $\bar{\psi}$. For nonlinear motion, the $z_{1}$-$z_{2}$
distribution exhibits a stable path – the limit cycle.

Figure 3: Amplitude and phase delay of nonlinear motion. a. The spectrogram of
the particle motion reveals the mode eigenfrequencies as a function of
$\Delta\Omega/\gamma$, where from the fit (black lines) we obtain a coupling
of $g_{a}/\gamma=0.58\pm 0.05$. The frequencies are degenerate for
$0<\Delta\Omega/\gamma\lesssim 2.3$. b. Measured displacement amplitude
$\bar{A}$ (points) and theoretical limit cycle amplitude $A$ (line) as a
function of $\Delta\Omega/\gamma$. The displacement amplitude becomes nonzero
within the $\mathcal{PT}$ symmetry-broken phase whenever the particles’ motion
are limit cycles. c. Probability density of the phase delay $\bar{\psi}$ is
shown in color, where black points show the obtained $\bar{\psi}_{\text{max}}$
at each detuning $\Delta\Omega/\gamma$. For large detunings
$|\Delta\Omega/\gamma|\gg 0$, histograms of $\bar{\psi}$ follow an
approximately uniform distribution. Within the $\mathcal{PT}$ symmetry-broken
phase, $\bar{\psi}$ acquires a strongly preferred value that depends on
$\Delta\Omega$. The solid black line is the theory plot that combines the
linear and nonlinear models for $\psi$.

We compare the theoretical steady-state solutions $A$ and $\psi$ to the
measured displacement amplitude $\bar{A}$ and phase delay $\bar{\psi}$ as a
function of the mechanical detuning $\Delta\Omega/\gamma$ (Fig. 3). We obtain
the coupling rate of $g_{a}/\gamma=0.58\pm 0.05$ from the fit of the
eigenfrequencies $\Omega_{\pm}$ to the spectrogram of particle motion (Fig.
3a). The region between the EPs, given by $0\leq\Delta\Omega/\gamma\lesssim
2.3$, defines the $\mathcal{PT}$ symmetry-broken phase. To obtain the
displacement amplitude $\bar{A}$ from the particle motion $z_{1,2}$, we first
reconstruct the complex amplitudes $a_{1,2}$ SI . We model the histograms of
$|a_{1,2}|$ with Rice distributions and fit the individual displacements
$\bar{A}_{1,2}$, which we use to calculate
$\bar{A}=2k\sqrt{\hbar/m\Omega_{0}}\sqrt{|\bar{A}_{2}|^{2}+|\bar{A}_{1}|^{2}}$.
Although there is a good quantitative agreement between the predicted $A$
(line) and reconstructed amplitudes $\bar{A}$ (points) (Fig. 3b), we attribute
the discrepancy between them to thermal fluctuations that are not included in
our theory model. As the phase is periodic with $2\pi$, we plot the histograms
of $\bar{\psi}\mod 2\pi$ (colored density plot), centered around the preferred
phase $\bar{\psi}_{\text{max}}$ (black points), as a function of
$\Delta\Omega/\gamma$ (Fig. 3c). Throughout the scan, the phase difference
undergoes a $\pi$ shift, as predicted by our model (black line) that combines
the nonlinear theory in the range where $A\neq 0$ and the full linear theory
with included thermal fluctuations elsewhere.

Figure 4: Non-Hermitian phases. a. The motion of particles 1 (blue) and 2
(orange) along the optical axes generates Stokes and anti-Stokes sidebands at
$\omega_{\text{L}}-\Omega_{0}$ and $\omega_{\text{L}}+\Omega_{0}$,
respectively, from the intrinsic laser frequency $\omega_{\text{L}}$. The
optical interaction leads to modified amplitudes of Stokes and anti-Stokes
sidebands of the eigenmodes $a_{1}-ia_{2}$ (red) and $a_{1}+ia_{2}$ (green),
amplifying and damping the modes with suppressed anti-Stokes (red) and Stokes
sideband (green), respectively. b. Expected amplitude $A$ of the limit cycle
as a function of the detuning $\Delta\Omega/\gamma$ and the anti-reciprocal
coupling $g_{a}/\gamma$. The oscillators are in the $\mathcal{PT}$ symmetry-
broken phase (hatched region) between the exceptional points (dashed lines).
Outside of this region, the oscillators are in the $\mathcal{PT}$-symmetric
phase. The black line marks the Hopf bifurcation of the nonlinear model. The
collective motion follows the limit cycles in regions with nonzero amplitude
$A$ (shades of blue). c. Measured displacement amplitude $\bar{A}$, where each
pixel represents a single measurement. d. (e.) The autocorrelation function
$g^{(1)}(\Delta t)$ shows oscillatory (bi-exponential) behavior in the
$\mathcal{PT}$-symmetric (symmetry-broken) phase. The shaded region in b marks
the $\mathcal{PT}$ symmetry-broken phase if the measured autocorrelation
function $g^{(1)}(\Delta t)$ shows a bi-exponential decay. f. The region of
the limit cycle solutions ($\bar{A}>0$) occurs at a threshold coupling rate of
$g_{a}/\gamma>0.5$ for $\Delta\Omega/\gamma=0.9$. The points are measured
$\bar{A}$ and the line represents a theory fit.

#### Mechanical lasing transition.

The Hopf bifurcation and the increase of the limit cycle amplitude $A$ may be
interpreted as a mechanical lasing transition. In analogy to a laser, here the
two-level laser medium is represented by the highly populated tweezer (upper
level) and the Stokes sideband given by the particle motion (lower level).
Scattering of the tweezer mode into the Stokes sideband creates a phonon in
the mechanical degrees of freedom of the particles. In past works, one or more
cavities were required to suppress the opposite transition into the anti-
Stokes sideband Zhang  _et al._ (2018). Instead of enhancing the Stokes
scattering via the Purcell effect, here an interference between the mechanical
sidebands created by the particles’ motion leads to a suppression of the
Stokes or anti-Stokes scattering, depending on which mechanical mode is
excited (Fig. 4a). This interference leads to an amplified (suppressed anti-
Stokes scattering) and a decaying mechanical mode (suppressed Stokes
scattering). The lasing threshold is reached when the amplification exceeds
the net loss, given by the intrinsic damping of the mechanical modes.

To fully characterize the mechanical lasing transition, we repeat the
measurements of $\bar{A}$ as a function of $\Delta\Omega/\gamma$ for varying
couplings $g_{a}/\gamma$. Comparison of $A$ from the nonlinear theory model in
Eq. (3) (Fig. 4b) to the extracted $\bar{A}$ (Fig. 4c) shows a good agreement,
which demonstrates that the observed nonlinear dynamics is well-described by
our theory model. We distinguish three regions in both graphs. For
$0<\Delta\Omega/\gamma<4g_{a}/\gamma$ (hatched area) the oscillators are in
the $\mathcal{PT}$ symmetry-broken phase, while they are in the
$\mathcal{PT}$-symmetric phase elsewhere. Each particle motion is a
combination of the two eigenmodes, which have non-degenerate frequencies
outside and non-degenerate damping rates inside the $\mathcal{PT}$ symmetry-
broken phase. Therefore, the crossing between these phases is characterized by
the change of the amplitude correlation functions $g_{jj}^{(1)}(\Delta
t)=\langle a_{j}^{*}(t)a_{j}(t+\Delta t)\rangle$ from an oscillatory behavior
outside (Fig. 4d), to a bi-exponential decay inside the symmetry-broken phase
(Fig. 4e). Furthermore, above the threshold marked by the solid black line,
the oscillators exhibit the mechanical lasing transition of Eq. (3). To
clearly show the transition, we plot $\bar{A}$ as a function of the coupling
at a detuning of $\Delta\Omega/\gamma\approx 0.9$ (Fig. 4f). A nonzero
$\bar{A}$ emerges at the threshold coupling of $g_{a}/\gamma\approx 0.5$ and
follows a square-root dependence on the coupling, as predicted by the theory.
The mechanical lasing transition could be understood as a second-order (non-
Hermitian) phase transition Fruchart  _et al._ (2021), bearing in mind that
our system does not exhibit an apparent thermodynamic limit.

#### Conclusion.

We investigated the collective linear and nonlinear dynamics of two anti-
reciprocally coupled optically levitated nanoparticles, which interact through
light-induced dipole-dipole forces. The particles’ motion was monitored
independently, which allowed us to reconstruct the eigenmodes of the system
and to observe signatures of the $\mathcal{PT}$ symmetry breaking. Within the
$\mathcal{PT}$ symmetry-broken phase, two eigenmodes with different damping
rates emerge and the particles’ motion become strongly correlated, which we
confirm by measuring a stable phase delay between the oscillators. For a
sufficiently high coupling rate, the system shows a mechanical lasing
transition where the joint phase space distribution exhibits a limit cycle.

The presented steady-state measurements are a first step toward probing
effects arising from the dynamical operation of nonreciprocally coupled
particle chains induced by, for example, encircling the exceptional point for
topological energy transfer Xu  _et al._ (2016); Patil  _et al._ (2022);
Doppler  _et al._ (2016). Larger arrays of tunable, nonreciprocally
interacting optically levitated particles have a great potential for studies
of nonreciprocal phase transitions Fruchart  _et al._ (2021), nonequilibrium
physics Loos and Klapp (2020); Xu  _et al._ (2023); Liu  _et al._ (2020), and
may have sensing applications McDonald and Clerk (2020); Malz  _et al._
(2019); Porras and Fernández-Lorenzo (2019); Sharma  _et al._ (2022); Wiersig
(2020); Li  _et al._ (2023). A combination of long-range nonreciprocal
interactions and the already established quantum control of particle motion
Delić  _et al._ (2020) holds promise to be a game-changer for studying
collective quantum phenomena, such as non-Hermitian quantum physics with long-
range interactions Lourenço  _et al._ (2022), entanglement in an optical
cavity, and observation of quantum optical binding Rudolph  _et al._ (2023).

###### Acknowledgements.

Acknowledgments. We thank Tobias Donner for insightful discussions, Markus
Aspelmeyer for his support, Helmut Ritsch for his idea of how to calibrate the
particle distance, and Baptiste Courme and Iurie Coroli for their help with
programming the acousto-optical deflectors. This research was funded in whole
or in part by the Austrian Science Fund (FWF, Project No. I 5111-N). A.V.Z.
acknowledges support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant LOREN (grant
agreement ID: 101030987). M.A. acknowledges support from the European Union’s
Horizon 2022 research and innovation programme under the Marie Sklodowska-
Curie grant NEOVITA (grant agreement ID: 101109773). H.R., K.H., B.A.S.
acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Grant No. 439339706. B.A.S. acknowledges support by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Grant No. 510794108. For the purpose of Open Access, the author has applied a
CC BY public copyright license to any Author Accepted Manuscript (AAM) version
arising from this submission.

Note. We are aware of related work by Liška et al Liška  _et al._ (2023).

## References

  * Fruchart  _et al._ (2021) M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Nature 592, 363 (2021).
  * Feng  _et al._ (2017) L. Feng, R. El-Ganainy, and L. Ge, Nature Photonics 11, 752 (2017).
  * Xiao  _et al._ (2020) L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Nature Physics 16, 761 (2020).
  * Weidemann  _et al._ (2020) S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Science 368, 311 (2020).
  * Wang  _et al._ (2021) K. Wang, A. Dutt, K. Y. Yang, C. C. Wojcik, J. Vučković, and S. Fan, Science 371, 1240 (2021).
  * Guo  _et al._ (2009) A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Phys. Rev. Lett. 103, 093902 (2009).
  * Zhen  _et al._ (2015) B. Zhen, C. W. Hsu, Y. Igarashi, L. Lu, I. Kaminer, A. Pick, S.-L. Chua, J. D. Joannopoulos, and M. Soljačić, Nature 525, 354 (2015).
  * Miri and Alù (2019) M.-A. Miri and A. Alù, Science 363, eaar7709 (2019).
  * Li  _et al._ (2019) J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Nature Communications 10, 855 (2019).
  * Gou  _et al._ (2020) W. Gou, T. Chen, D. Xie, T. Xiao, T.-S. Deng, B. Gadway, W. Yi, and B. Yan, Phys. Rev. Lett. 124, 070402 (2020).
  * Takasu  _et al._ (2020) Y. Takasu, T. Yagami, Y. Ashida, R. Hamazaki, Y. Kuno, and Y. Takahashi, Prog Theor Exp Phys 2020, 12A110 (2020).
  * Öztürk  _et al._ (2021) F. E. Öztürk, T. Lappe, G. Hellmann, J. Schmitt, J. Klaers, F. Vewinger, J. Kroha, and M. Weitz, Science 372, 88 (2021).
  * Ferri  _et al._ (2021) F. Ferri, R. Rosa-Medina, F. Finger, N. Dogra, M. Soriente, O. Zilberberg, T. Donner, and T. Esslinger, Physical Review X 11, 041046 (2021).
  * Liang  _et al._ (2022) Q. Liang, D. Xie, Z. Dong, H. Li, H. Li, B. Gadway, W. Yi, and B. Yan, Phys. Rev. Lett. 129, 070401 (2022).
  * Gaikwad  _et al._ (2023) C. Gaikwad, D. Kowsari, W. Chen, and K. W. Murch, (2023), arXiv:2306.14980 [quant-ph] .
  * Helbig  _et al._ (2020) T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. W. Molenkamp, C. H. Lee, A. Szameit, M. Greiter, and R. Thomale, Nature Physics 16, 747 (2020).
  * Zou  _et al._ (2021) D. Zou, T. Chen, W. He, J. Bao, C. H. Lee, H. Sun, and X. Zhang, Nature Communications 12, 7201 (2021).
  * Brandenbourger  _et al._ (2019) M. Brandenbourger, X. Locsin, E. Lerner, and C. Coulais, Nature Communications 10, 4608 (2019).
  * Ghatak  _et al._ (2020) A. Ghatak, M. Brandenbourger, J. van Wezel, and C. Coulais, Proceedings of the National Academy of Sciences 117, 29561 (2020).
  * Chen  _et al._ (2021) Y. Chen, X. Li, C. Scheibner, V. Vitelli, and G. Huang, Nature Communications 12, 5935 (2021).
  * Wang  _et al._ (2022) W. Wang, X. Wang, and G. Ma, Nature 608, 50 (2022).
  * Metelmann and Clerk (2015) A. Metelmann and A. A. Clerk, Physical Review X 5, 021025 (2015).
  * Fang  _et al._ (2017) K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Nature Physics 13, 465 (2017).
  * Sanavio  _et al._ (2020) C. Sanavio, V. Peano, and A. Xuereb, Phys. Rev. B 101, 085108 (2020).
  * Metelmann and Clerk (2014) A. Metelmann and A. A. Clerk, Phys. Rev. Lett. 112, 133904 (2014).
  * Lau and Clerk (2018) H.-K. Lau and A. A. Clerk, Nature Communications 9, 4320 (2018).
  * McDonald and Clerk (2020) A. McDonald and A. A. Clerk, Nature Communications 11, 5382 (2020).
  * Peano  _et al._ (2015) V. Peano, C. Brendel, M. Schmidt, and F. Marquardt, Physical Review X 5, 031011 (2015).
  * Rosenthal  _et al._ (2018) E. I. Rosenthal, N. K. Ehrlich, M. S. Rudner, A. P. Higginbotham, and K. W. Lehnert, Phys. Rev. B 97, 220301 (2018).
  * Wanjura  _et al._ (2020) C. C. Wanjura, M. Brunelli, and A. Nunnenkamp, Nature Communications 11, 3149 (2020).
  * Mathew  _et al._ (2020) J. P. Mathew, J. d. Pino, and E. Verhagen, Nature Nanotechnology 15, 198 (2020).
  * Xu  _et al._ (2016) H. Xu, D. Mason, L. Jiang, and J. G. E. Harris, Nature (London) 537, 80 (2016).
  * Patil  _et al._ (2022) Y. S. S. Patil, J. Höller, P. A. Henry, C. Guria, Y. Zhang, L. Jiang, N. Kralj, N. Read, and J. G. E. Harris, Nature (London) 607, 271 (2022).
  * Ren  _et al._ (2022) H. Ren, T. Shah, H. Pfeifer, C. Brendel, V. Peano, F. Marquardt, and O. Painter, Nature Communications 13, 3476 (2022).
  * Zhang  _et al._ (2022) Q. Zhang, C. Yang, J. Sheng, and H. Wu, Proceedings of the National Academy of Sciences 119, e2207543119 (2022).
  * Youssefi  _et al._ (2022) A. Youssefi, S. Kono, A. Bancora, M. Chegnizadeh, J. Pan, T. Vovk, and T. J. Kippenberg, Nature (London) 612, 666 (2022).
  * Doster  _et al._ (2022) J. Doster, T. Shah, T. Fösel, P. Paulitschke, F. Marquardt, and E. M. Weig, Nature Communications 13, 2478 (2022).
  * Liu  _et al._ (2023) T. Liu, J.-Y. Ou, K. F. MacDonald, and N. I. Zheludev, Nature Physics 19, 986 (2023).
  * Gonzalez-Ballestero  _et al._ (2023) C. Gonzalez-Ballestero, M. Aspelmeyer, L. Novotny, R. Quidant, and O. Romero-Isart, Science 374, eabg3027 (2023).
  * Stickler  _et al._ (2021) B. A. Stickler, K. Hornberger, and M. S. Kim, Nature Reviews Physics 3, 589 (2021).
  * Rieser  _et al._ (2022) J. Rieser, M. A. Ciampini, H. Rudolph, N. Kiesel, K. Hornberger, B. A. Stickler, M. Aspelmeyer, and U. Delić, Science 377, 987 (2022).
  * Vijayan  _et al._ (2023) J. Vijayan, Z. Zhang, J. Piotrowski, D. Windey, F. van der Laan, M. Frimmer, and L. Novotny, Nature Nanotechnology 18, 49 (2023).
  * Yan  _et al._ (2023) J. Yan, X. Yu, Z. V. Han, T. Li, and J. Zhang, Photon. Res. 11, 600 (2023).
  * Liška  _et al._ (2023) V. Liška, T. Zemánková, V. Svak, P. Jákl, J. Ježek, M. Bránecký, S. H. Simpson, P. Zemánek, and O. Brzobohatý, Optica 10, 1203 (2023).
  * Vijayan  _et al._ (2023) J. Vijayan, J. Piotrowski, C. Gonzalez-Ballestero, K. Weber, O. Romero-Isart, and L. Novotny, (2023), arXiv:2308.14721 [quant-ph] .
  * Bykov  _et al._ (2023) D. S. Bykov, L. Dania, F. Goschin, and T. E. Northup, Optica 10, 438 (2023).
  * Penny  _et al._ (2023) T. W. Penny, A. Pontin, and P. F. Barker, Physical Review Research 5, 013070 (2023).
  * Rudolph  _et al._ (2023) H. Rudolph, U. Delić, K. Hornberger, and B. A. Stickler, (2023), arXiv:2306.11893 [quant-ph] .
  * (49) See supplementary materials for the detailed experimental setup, comparison of the particle sizes, linear and nonlinear theory models, particle discharging, calibration of particle distance, and the heterodyne detection.
  * Zhang  _et al._ (2018) J. Zhang, B. Peng, c. K. Özdemir, K. Pichler, D. O. Krimer, G. Zhao, F. Nori, Y.-x. Liu, S. Rotter, and L. Yang, Nature Photonics 12, 479 (2018).
  * Doppler  _et al._ (2016) J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Nature (London) 537, 76 (2016).
  * Loos and Klapp (2020) S. A. M. Loos and S. H. L. Klapp, New Journal of Physics 22, 123051 (2020).
  * Xu  _et al._ (2023) H. Xu, U. Delić, G. Wang, C. Li, P. Cappellaro, and J. Li, (2023), arXiv:2309.07731 [quant-ph] .
  * Liu  _et al._ (2020) S. Liu, Z.-q. Yin, and T. Li, Advanced Quantum Technologies 3, 1900099 (2020).
  * Malz  _et al._ (2019) D. Malz, J. Knolle, and A. Nunnenkamp, Nature Communications 10, 3937 (2019).
  * Porras and Fernández-Lorenzo (2019) D. Porras and S. Fernández-Lorenzo, Phys. Rev. Lett. 122, 143901 (2019).
  * Sharma  _et al._ (2022) S. Sharma, A. Kani, and M. Bhattacharya, Phys. Rev. A 105, 043505 (2022).
  * Wiersig (2020) J. Wiersig, Nature Communications 11, 2454 (2020).
  * Li  _et al._ (2023) A. Li, H. Wei, M. Cotrufo, W. Chen, S. Mann, X. Ni, B. Xu, J. Chen, J. Wang, S. Fan, C.-W. Qiu, A. Alù, and L. Chen, Nature Nanotechnology 18, 706 (2023).
  * Delić  _et al._ (2020) U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Science 367, 892 (2020).
  * Lourenço  _et al._ (2022) J. A. S. Lourenço, G. Higgins, C. Zhang, M. Hennrich, and T. Macrı, Phys. Rev. A 106, 023309 (2022).
  * Liška  _et al._ (2023) V. Liška, T. Zemánková, P. Jákl, M. Šiler, S. H. Simpson, P. Zemánek, and O. Brzobohatý, (2023), arXiv:2310.03701 [physics.optics] .
  * Dholakia and Zemánek (2010) K. Dholakia and P. Zemánek, Rev. Mod. Phys. 82, 1767 (2010).
  * Tebbenjohanns  _et al._ (2019) F. Tebbenjohanns, M. Frimmer, and L. Novotny, Phys. Rev. A 100, 043821 (2019).
  * Magrini  _et al._ (2021) L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Nature 595, 373 (2021).
  * Tebbenjohanns  _et al._ (2021) F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Nature 595, 378 (2021).
  * Gieseler  _et al._ (2013) J. Gieseler, L. Novotny, and R. Quidant, Nature Physics 9, 806 (2013).
  * Aspelmeyer  _et al._ (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * Kepesidis  _et al._ (2016) K. V. Kepesidis, T. J. Milburn, J. Huber, K. G. Makris, S. Rotter, and P. Rabl, New Journal of Physics 18, 095003 (2016).
  * Šantić  _et al._ (2018) N. Šantić, A. Fusaro, S. Salem, J. Garnier, A. Picozzi, and R. Kaiser, Phys. Rev. Lett. 120, 055301 (2018).

##  Appendix A Supplementary Materials

###  A.1 Theoretical model

In this section, we discuss the derivation of the linear and nonlinear models
(1) and (3) for the collective motion of the two trapped nanoparticles from
the theory of optical interactions Dholakia and Zemánek (2010); Rudolph  _et
al._ (2023). The starting point is the quantum Langevin equations from Ref.
Rudolph  _et al._ (2023) for light-induced dipole-dipole interactions between
arbitrary mechanical degrees of freedom. We consider two identical spherical
particles in their respective tweezer traps (see Eq. (69) Rudolph  _et al._
(2023)), ignore quantum noise, treat all operators as classical variables, and
add gas damping and diffusion with damping constant $\gamma$ and gas
temperature $T_{\rm g}$. The tweezers are assumed to have identical (linear)
polarizations and waists but may be driven with a different power. The
mechanical frequencies along transverse directions are significantly higher
than the frequency along the optical axis, therefore we restrict the dynamics
to the particle motion along the tweezer propagation direction.

To determine the main nonlinearity of the two-particle dynamics around the
tweezer foci, we expand the optical potential and the non-conservative
radiation pressure force to the first-order nonlinearity, while in the dipole-
dipole interaction, the tweezer fields are well-approximated by plane waves.
Furthermore, the equations of motion are expanded to leading order in the
inverse distance $1/d_{0}$ between the tweezer foci, since far-field coupling
dominates.

The resulting equations of motion for the particle coordinates $z_{1,2}$ and
momenta $p_{1,2}$ along the optical axis can be written in terms of the
complex amplitudes $a_{1,2}$, as defined below of Eq. (1). For the experiment
described in this article, the influences of gas collisions, nonlinearities in
the optical potential, the non-conservative radiation pressure forces, dipole-
dipole interactions as well as the mechanical frequency difference
$\Delta\Omega$ are weak during one cycle of the mean harmonic motion
$\Omega_{0}$. This allows us to average the equations of motion for $a_{1,2}$
over one mechanical period $2\pi/\Omega_{0}$, resulting in

$\dot{a}_{j}=\pm
i\frac{\Delta\Omega}{2}a_{j}-\frac{\gamma}{2}a_{j}+i\beta|a_{j}|^{2}a_{j}+\sqrt{\gamma
n_{T}}\xi_{j}(t)\pm i(g_{r}\pm g_{a})(a_{2}-a_{1})f(2kz_{\rm
zpf,0}|a_{2}-a_{1}|),$ (S1)

where the upper (lower) sign holds for $j=1$ $(j=2)$. Here, we define $z_{\rm
zpf,0}=\sqrt{\hbar/2m\Omega_{0}}$, the thermal occupation $n_{T}=k_{\rm
B}T_{\rm g}/\hbar\Omega_{0}$, the optical potential nonlinearity
$\beta=3\Omega_{0}z_{\rm zpf,0}^{2}/z_{\rm R}^{2}$, depending on the Rayleigh
range $z_{\rm R}$ of the tweezers, and the reciprocal and anti-reciprocal
coupling rates $g_{r}=G\cos{(kd_{0})}\cos{(\Delta\phi)}/kd_{0}$ and
$g_{a}=G\sin{(kd_{0})}\sin{(\Delta\phi)}/kd_{0}$ as in the main text. They
depend on the constant

$\displaystyle G=\frac{\varepsilon_{0}\chi^{2}V^{2}k^{5}E_{0}^{2}}{16\pi
m\Omega_{0}}\cos^{2}\theta,$ (S2)

with the particle volume $V$ and electric susceptibilty
$\chi=3(\varepsilon_{r}-1)/(\varepsilon_{r}+2)$ with relative electric
permeability $\varepsilon_{r}$, the local tweezer field strength $E_{0}$ and
the angle ($\pi/2-\theta$) between tweezer polarization and particle-
connecting axis. Importantly, the relevant nonlinearity in Eqs. (3), (S1)
causing the limit cycle is not the optical potential nonlinearity $\beta$, but
the nonlinearity in the dipole-dipole interaction described by the function
$f$. As a result of the average over one mechanical cycle, (i) the non-
conservative radiation pressure forces cancel, and (ii) the gas diffusion
noises $\xi_{j}(t)$ turn complex-valued.

First order correlation functions.– To obtain the linear model (1) for deeply
trapped particles, we harmonically expand Eq. (S1) around $a_{1,2}=0$, using
$f(0)=1$. Then, the autocorrelation functions $g^{(1)}_{jj}=\langle
a^{*}_{j}(t)a_{j}(t+\Delta t)\rangle$ for $j\in\\{1,2\\}$, as well as the
cross-correlation function $g^{(1)}_{12}=\langle a^{*}_{2}(t)a_{1}(t+\Delta
t)\rangle$ of deeply trapped particles follow from a straightforward
calculation, yielding

$\displaystyle g_{jj}^{(1)}(\Delta
t)=n_{T}\Bigg{[}\Bigg{(}1+\frac{4g_{a}(g_{a}\pm
g_{r})}{\gamma^{2}+\kappa^{2}}\Bigg{)}\cos\left(\frac{\kappa}{2}\Delta
t\right)+\frac{4\gamma g_{a}(g_{a}\pm
g_{r})}{\kappa(\gamma^{2}+\kappa^{2})}\sin\left(\frac{\kappa}{2}|\Delta
t|\right)\pm
i\frac{\Delta\Omega-2g_{a}}{\kappa}\sin\left(\frac{\kappa}{2}\Delta
t\right)\Bigg{]}e^{-ig_{r}\Delta t}e^{-\frac{\gamma}{2}|\Delta t|},$ (S3)
$\displaystyle g_{12}^{(1)}(\Delta
t)=n_{T}\Bigg{[}\frac{2g_{a}(i\gamma+2g_{a}-\Delta\Omega)}{\gamma^{2}+\kappa^{2}}\Bigg{\\{}\cos\left(\frac{\kappa}{2}\Delta
t\right)+\frac{\gamma}{\kappa}\sin\left(\frac{\kappa}{2}|\Delta
t|\right)\Bigg{\\}}+i\frac{2g_{r}}{\kappa}\sin\left(\frac{\kappa}{2}\Delta
t\right)\Bigg{]}e^{-ig_{r}\Delta t}e^{-\frac{\gamma}{2}|\Delta t|},$ (S4)

with $\kappa=\sqrt{\Delta\Omega^{2}-4g_{a}\Delta\Omega+g_{r}^{2}}$. In the
$\mathcal{PT}$ symmetry-broken regime, where $\kappa$ becomes imaginary, the
oscillatory behavior in Eq. (S3) changes to a bi-exponential decay, in
excellent agreement with the observed correlation functions (Fig. 4d-e).

Limit cycle model.– The limit cycle model (3) follows from Eq. (S1) by
ignoring gas diffusion, choosing the anti-reciprocal interaction as $g_{r}=0$
and writing the complex amplitudes in terms of occupations and phases as
$a_{1,2}=\sqrt{n_{1,2}}\exp{(-i\psi_{1,2})}$. It follows that the occupation
difference $n_{2}-n_{1}$ is always damped to zero on the timescale of the gas
damping constant $\gamma$, while both occupations $n_{1,2}$ and the mechanical
phase delay $\psi=\psi_{2}-\psi_{1}$ decouple from the mean mechanical phase
$(\psi_{1}+\psi_{2})/2$. Assuming that $n_{2}-n_{1}$ has already settled at
zero, and defining the dimensionless effective oscillation amplitude by
$A=2\sqrt{2}kz_{\rm zpf,0}\sqrt{n_{1}+n_{2}}$, we finally arrive at Eq. (3).
Note that the optical potential nonlinearity $\beta$ cancels from the
effective limit cycle model (3), as its influence on the phase delay $\psi$
vanishes if the particles oscillate with identical amplitude.

###  A.2 Detection of the particle motion

The particles scatter light and imprint the information about their position
onto the phase of the scattered light Tebbenjohanns  _et al._ (2019). In our
case, the optical phase of the backscattered light optimally maps the
particles’ motion along the optical axis (Fig. S1a). The optical phase
$\varphi_{j}$ can be approximated as the sum of contributions from the three-
dimensional motion:

$\varphi_{j}(r_{j}(t))=\varphi_{j}(z_{j})+\varphi_{j}(y_{j})+\varphi_{j}(x_{j}).$
(S5)

However, the $z$-motion is well-separated in frequency from the $x$\- and
$y$-motion due to the trap asymmetry
($\Omega_{z,j}\ll\Omega_{x,j},\Omega_{y,j}$). Therefore, we can isolate the
$z$-motion contribution by applying a bandpass filter with the bandwidth
$[5,50]$ kHz around the average $z$ frequency $\Omega_{0}/2\pi=27.5$ kHz.
Furthermore, the transverse motion influences the optical phase around $100$
times less than $z_{j}$ in the backplane detection.

To split the backscattered light from the trapping light, we use a polarizing
beamsplitter and a Faraday rotator (Fig. S1b). The image of the particle plane
is magnified $160$-fold onto a prism mirror that separates the individual
particle images into opposite directions. The split particle signals are
collected with two single-mode fibers such that the confocal microscope
magnification ($\times 7$) between the trapping plane and the fiber input
plane maximizes the collection efficiency Magrini  _et al._ (2021);
Tebbenjohanns  _et al._ (2021). The collected scattered light (optical power
of $5$-$10$ $\mu$W) is mixed on a beamsplitter with a reference beam (local
oscillator) (optical power of $1$-$2$ mW, frequency-shifted by
$\omega_{\text{LO}}/2\pi=1.1$ MHz) to implement a balanced heterodyne
detection.

Figure S1: Detection of the particles’ motion. a. A particle is trapped in a
tightly focused Gaussian beam at a position $z_{0,j}$ away from the focus. The
phase of the backscattered light changes with the particle motion $z_{j}$
around the trap position as $\varphi_{j}\approx 2kz_{j}$. b. Sketch of the
balanced heterodyne detection scheme. The backscattered light is separated
from the incoming beams by a Faraday rotator (FR) and a polarizing
beamsplitter (PBS) and subsequently reimaged onto a prism mirror to spatially
separate the light scattered by the two particles. Each particle’s image is
collected with a single-mode fiber in a confocal microscope configuration for
further spatial filtering from the background reflections. The particles’
signals are combined with a local oscillator at a frequency of
$\omega_{\text{LO}}$ with a 50:50 beamsplitter and are measured with balanced
photodetectors.

The heterodyne signal to the leading order in the particle motion $z_{j}$
yields:

$I_{-,j}(t)\propto\left|E_{LO}^{*}E_{sc}\right|\cos\left(\omega_{LO}t+\varphi_{j}[z_{j}(t)]\right),$
(S6)

where $E_{LO}$ and $E_{sc}$ are the strengths of the local oscillator and
scattered fields, respectively. We retrieve the optical phases $\varphi_{j}$
by calculating the Hilbert transform of $I_{-,j}$, which allows us to obtain
the unwrapped argument of $I_{-,j}$. Finally, we subtract the known
accumulated phase of the local oscillator to obtain:

$\varphi_{j}(z_{j}(t))=\textrm{unwr}\left\\{\textrm{arg}\left(\mathcal{H}\left[I_{-,j}(t)\right]\right)\right\\}-\omega_{LO}t.$
(S7)

The obtained phase can now be converted into the position $z_{j}$ as
Tebbenjohanns  _et al._ (2019):

$\varphi_{j}(z_{j})=2kz_{j}-\textrm{arctan}\left((z_{j}+z_{0,j})/z_{\rm
R}\right),$ (S8)

where $\textrm{arctan}\left((z_{j}+z_{0,j})/z_{\rm R}\right)$ is the
modification due to the Gouy phase with the constant offset from the focus
$z_{0,j}$, the Rayleigh length $z_{\rm R}=w_{x}w_{y}\pi/\lambda$ and waists
$w_{x}$ and $w_{y}$ along the $x$\- and $y$-axis, respectively. We numerically
calculate the tweezer waists $w_{x}=0.678\ \mu$m and $w_{y}=0.775\ \mu$m from
the measured mechanical frequencies and $z_{0,j}=0.8\ \mu$m from the particle
radius to estimate the correction due to the Gouy phase. The obtained Rayleigh
length of $z_{\rm R}=1.55\ \mu$m results in the modification of $\sim 3\%$ to
the particle motion detection. Therefore, we safely omit this correction and
use $\varphi_{j}(z_{j})=2kz_{j}$ instead in the main text.

###  A.3 Comparison of particle sizes

We assume identical spherical particles in size and mass in the theoretical
models. The particles used in the experiment have a nominal radius of
$r=(105\pm 2)$ nm (microParticles GmbH), therefore they could be of slightly
different sizes. In the regime dominated by gas damping, the damping rate of
the particle CoM motion depends on the particle radius as $\gamma\propto 1/r$
Gieseler  _et al._ (2013). Therefore, we measure and compare the damping rates
of the two particles as a function of gas pressure to determine the ratio of
the particle radii Yan  _et al._ (2023). At each pressure, we calculate the
PSDs of the CoM motion and extract the gas damping rates of particles 1 (blue)
and 2 (purple) from the fit of the driven damped harmonic oscillator (Fig.
S2a) Aspelmeyer  _et al._ (2014). We assume equal mass densities and, due to
the absence of rotational signatures in the PSDs, spherical particles. This
relates the ratio of the damping rates of the two particles simply to the
inverse ratio of the particle radii $\gamma_{1}/\gamma_{2}=r_{2}/r_{1}$. We
obtain the average ratio of the damping rates $\gamma_{1}/\gamma_{2}=0.99\pm
0.01$, and determine the size ratio to be $r_{1}/r_{2}=1.01\pm 0.01$ (Fig.
S2b). We conclude that the particles can be considered of approximately equal
sizes.

Figure S2: Ratios of particle damping rates at different pressures. a.
Calculated particle PSDs at three different pressures for particles 1 (blue)
and 2 (purple) illustrate how the damping rates $\gamma_{1,2}$ decrease as the
gas pressure decreases. The black lines correspond to the fitted lineshapes.
b. Ratios of the fitted damping rates at different pressures (blue points and
error bars) yields the weighted average (orange line) and its standard
deviation (orange shaded region) of $\gamma_{1}/\gamma_{2}=0.99\pm 0.01$.

###  A.4 Calibration of the particle distance

The mean particle distance $d_{0}$, given by the positions of the foci of
their respective tweezers, influences the reciprocal and anti-reciprocal
coupling rates. In our experiment, we tune the frequency difference
$\Delta\omega$ of the RF tones to modify the distance as
$d_{0}=C_{0}\Delta\omega$, where $C_{0}$ is the conversion factor that can be
calibrated by comparing the distance to a known etalon, e.g. the laser
wavelength. Therefore, we perform an interferometric measurement of the dipole
radiation from the two particles with a detector placed in the far field,
along the axis connecting the particles (Fig. S3). The tweezer polarizations
are set perpendicular to this axis, such that the dipole radiation in the
direction of the detector is maximal. As the images of the two particles are
indistinguishable at the detector, the overlap of the dipole radiations leads
to an interference that depends on $d_{0}$. In our measurement, we sweep
$\Delta\omega/2\pi$ in the range of $4-12$ MHz and record the intensity of the
interference pattern. We note that the magnitude of the interference pattern
changes as a function of $\Delta\omega$. As the measurement is fully
repeatable, we attribute this to the varying AOD diffraction efficiencies as a
function of the RF tone frequencies. We fit the data with a sine function and
obtain the conversion factor $C_{0}=(2.3845\pm 0.0005)\>\mu$m/MHz. All
measurements in the main text were performed at a frequency difference of
$\Delta\omega/2\pi=7.72$ MHz, which leads to an estimated distance of
$d_{0}=(18.408\pm 0.005)\>\mu$m.

Figure S3: Distance calibration. The scattered light from particles 1 and 2
interferes in the far field. We measure the intensity of the interference
pattern while scanning the particle distance (black line). The fringe spacing
is given by the wavelength of the trapping light. We fit the pattern (red
line) to determine the calibration constant from the frequency difference to
particle distance.

###  A.5 Measurement of the coupling rates

In this section, we detail the experimental methods to optimize the anti-
reciprocal coupling rate $g_{a}$ and to reduce the reciprocal coupling rate
$g_{r}$ and the coupling rate $g_{C}$ due to electrostatic interactions.

Optical coupling— The optical coupling rates $g_{r}$ and $g_{a}$ depend on the
particle distance $d_{0}$ and the optical phase difference $\Delta\phi$. To
maximize $g_{a}$ and minimize $g_{r}$, we scan both parameters while the trap
polarization is set to $\sim 0^{\circ}$ to maximize the optical interaction.
For a general combination of parameters where $g_{a}\neq 0$ and $g_{r}\neq 0$,
the total coupling rates $g_{r}\pm g_{a}$ are different in opposite
directions. In this case, the response to the motion of particle 2 in the PSD
of particle 1, with a peak amplitude $\propto(g_{r}+g_{a})^{2}$ at frequency
$\Omega_{2}$ in the PSD, is different from the peak amplitude
$\propto(g_{r}-g_{a})^{2}$ of the particle 1 at frequency $\Omega_{1}$ in the
PSD of particle 2. For any combination of $d_{0}$ and $\Delta\phi$, we monitor
the PSD of both particles and we choose the parameters such that the two
contributions have matching amplitudes, which is satisfied either for
$g_{r}=0$ or $g_{a}=0$. As the final step, we make sure that there is no
normal mode splitting but normal mode attraction to confirm that $g_{r}\equiv
0$. For such a configuration, we estimate that $\Delta\phi=(0.5\pm 0.08)\>\pi$
due to the step size of the scan. With the previously estimated distance of
$d_{0}=(18.408\pm 0.005)\>\mu$m, the reciprocal coupling
$g_{r}\propto\cos(kd_{0})\cos(\Delta\phi)=0.00\pm 0.08$ is effectively zero.

Figure S4: Anti-reciprocal coupling rate. a. We fit the eigenfrequency
difference $\Omega_{+}-\Omega_{-}$ as a function of the intrinsic frequency
detuning $\Delta\Omega$ (points) with linear theory (line) to extract the
coupling rate $g_{a}$. b. The anti-reciprocal coupling rate follows a
$\cos^{2}\theta$ dependence of the tweezer polarization angle $\theta$. The
red point represents the detuning scan in a.

To determine the magnitude of the anti-reciprocal coupling rate $g_{a}$, we
fit the difference of the eigenfrequencies
$\Omega_{+}-\Omega_{-}=\operatorname{Re}(\sqrt{\Delta\Omega^{2}-4g_{a}\Delta\Omega})$,
obtained from the PSDs, as a function of the intrinsic frequency detuning
$\Delta\Omega$ (Fig. S4a). We attribute the nonzero values in the
$\mathcal{PT}$ symmetry-broken phase to small differences in the particle
masses and the thermal fluctuations, which are not accounted for in the linear
model of the eigenfrequencies. We determine $g_{a}$ as a function of the
tweezer polarization $\theta$ (Fig. S4b), which follows a $\cos^{2}\theta$
function, as predicted from the light-induced dipole-dipole interactions.

Electrostatic coupling— To minimize the electrostatic interaction, we
discharge the particles by igniting a plasma at a high-voltage pin inside the
vacuum chamber (Fig. S5a). We attach wires to the holders of the trapping lens
and the lens that collects light after the tweezer foci, which forms a
capacitor for the charged particle with holders acting as the electrodes. We
continuously drive the particle motion with a sinusoidal signal applied to two
electrodes, which yields a sharp peak at the drive frequency in the PSD of the
particle motion. When plasma is ignited, the number of charges on each
particle fluctuates, which we monitor as changes in the peak amplitude in
real-time. We switch off the high voltage when both particles have just a few
charges, which is reflected in low amplitudes of the drive peak in the PSD.

Figure S5: Charge control and calibration. a. To change the amount of net
charges attached to the particles we ignite a plasma at a high-voltage pin
inside the vacuum chamber. We monitor the number of charges by driving the
particles via an AC field, applied to electrodes mounted to the trapping lens
(signal, red line) and a lens tube (ground, dark blue line), at a relative
distance of $D_{e}=(5\pm 1)$ mm. b. The electric drive can be seen as a sharp
peak (red) in the power spectral densities if the particles are charged. We
estimate the remaining charges to be $q_{1}=(2.5\pm 0.7)\>e$ and
$q_{2}=(1.5\pm 0.5)\>e$, where $e$ is the elementary charge.

To estimate the electrostatic interaction after the discharging procedure, we
follow the procedure in the supplementary information of Rieser  _et al._
(2022). To calibrate the number of charges on each particle, we drive the
motion of both particles by applying a sinusoidal voltage with an amplitude of
$V_{AC}=10$ V and at a frequency $\Omega_{d}/2\pi=31$ kHz to the electrodes
and record the particle motion. Assuming a model of a point charge in a plate
capacitor, we estimate the number of charges $N_{1,2}$ to be:

$N_{1,2}=\frac{\sqrt{2}mD_{e}}{eV_{AC}}\left|\Omega^{2}_{1,2}-\Omega^{2}_{d}\right|\sqrt{\left\langle
z^{2}_{d1,2}\right\rangle},$ (S9)

where $m$ is the particle mass, $D_{e}=(5\pm 1)$ mm is the distance between
the electrodes, $e$ is the elementary charge, $\Omega_{1,2}$ are the
oscillation frequencies of particles 1 and 2, and $\langle
z^{2}_{d1,2}\rangle$ is the time-averaged displacement of particle 1,2 along
the $z$-axis due to the electronic drive. The PSDs of the driven motion can be
seen in Fig. S5b. The peak due to the electric drive is marked in red. We
estimate the number of charges to be $N_{1}=2.5\pm 0.7$ and $N_{2}=1.5\pm
0.5$, which gives an electrostatic coupling rate of:

$g_{C}=-\frac{N_{1}N_{2}e^{2}}{8\pi\varepsilon_{0}m\Omega^{\prime}d^{3}_{0}}=(-0.077\pm
0.030)\>\text{Hz}.$ (S10)

Here, $\varepsilon_{0}$ is the vacuum permittivity and
$\Omega^{\prime}=\sqrt{\Omega^{2}-N_{1}N_{2}e^{2}/(4\pi\varepsilon_{0}md^{3}_{0})}$.
This coupling rate is more than two orders of magnitude smaller than the
intrinsic damping rate and the measured anti-reciprocal coupling rates for
polarization angles in the range $[0^{\circ},80^{\circ}]$. Therefore, the
electrostatic interaction does not affect the non-Hermitian dynamics in our
experiment.

###  A.6 Evaluation of the phase delay and the displacement amplitude

Phase delay.– To compute the oscillation phases $\bar{\psi}_{1,2}$, we first
apply the Hilbert transform to the individual particle motion $z_{1,2}$ to
obtain the complex amplitude
$\tilde{z}_{1,2}(t)=\mathcal{H}\left[z_{1,2}\right]$ such that
$z_{1,2}(t)=\operatorname{Re}(\tilde{z}_{1,2}(t))$. The oscillation phases are
recovered from the unwrapped phases of the complex amplitudes:

$\bar{\psi}_{j}(z_{j}(t))=\textrm{unwr}\left\\{\textrm{arg}\left(\mathcal{H}\left[z_{j}(t)\right]\right)\right\\}.$
(S11)

We define the instantaneous phase delay between the oscillators as
$\bar{\psi}=\bar{\psi}_{2}-\bar{\psi}_{1}$. The interaction generates a
preferred phase delay $\bar{\psi}_{max}$ as observed in the histograms of
$\bar{\psi}$, while it is uniformly distributed in the case of uncoupled
oscillators (Fig. 2b).

Displacement amplitude.– The Nonlinear Interaction Regime (NIR) arises for the
negative effective damping $\gamma_{-}<0$ as predicted by the linearized
theory. This results in the modified statistics of the particles’ motion and
the emergence of a limit cycle (Fig. 2a). To analyze this regime
quantitatively and to compare the results to the theoretical model (3), we
calculate the histograms of the oscillation envelopes $|\tilde{z}_{j}|$ for
each particle and fit them with a model accounting for the possible emergence
of a limit cycle with a displacement amplitude $\bar{A}_{j}$ (Fig. S6a)
Kepesidis  _et al._ (2016).

Outside the NIR, there is no limit cycle, therefore $\bar{A}_{j}=0$. There,
the positions (as well as velocities) assume Gaussian distributions centered
around zero with a width $\sigma$ that depends on the effective thermal
occupation. The histograms of $|\tilde{z}_{j}(t)|$ are described by the
Rayleigh distribution $P_{\rm
Ray}(x,\sigma)=(x/\sigma^{2})\exp(-x^{2}/2\sigma^{2})$ as a function of the
position $x$. As soon as the limit cycle trajectories emerge, the displaced
amplitude $\bar{A}_{j}$ becomes nonzero. There, the oscillation envelope
follows a Rice distribution that generalizes the Rayleigh distribution for the
case of non-centered Gaussian random variables. The Rice distribution has
previously been used for the statistical description of precondensed light
Šantić  _et al._ (2018). In our case, it is defined as:

$P_{\rm
Rice}(x,\sigma,\bar{A}_{j})=\frac{x}{\sigma^{2}}e^{-\frac{x^{2}+\bar{A}_{j}^{2}}{2\sigma^{2}}}\mathcal{I}_{0}\left(\frac{x\bar{A}_{j}}{\sigma^{2}}\right),$
(S12)

where $\mathcal{I}_{0}$ is the modified Bessel function of the first kind with
order zero. Note that the thermal (coherent) statistics are recovered from the
Rice distribution in the limit of $\bar{A}_{j}\rightarrow 0$
($\bar{A}_{j}\gg\sqrt{\sigma}$).

Figure S6: Evaluation of the displacement amplitude. a. Short time trace of
the particle motion (blue), velocity (orange), and the reconstructed
oscillation envelope (black). b. Phase space distributions of the motion of
particle 1 sampled over the whole time trace for four different
configurations, labeled as (i)-(iv). The limit cycle is observed in (ii) and
(iii), while no limit cycle is observed in (iv). c. The probability
distribution functions (pdf) of the oscillation envelopes for configurations
(i)-(iii). Solid black lines represent the fits of the Rice distribution,
where (ii) and (iii) clearly show the increase of the displaced amplitude
$\bar{A}_{1}$ from the case of $\bar{A}_{1}=0$ in (i). The presented
experimental data corresponds to the following detuning and coupling rate
$(\Delta\Omega/2\pi,g_{a}/2\pi)$ and extracted displacement amplitudes
$\bar{A}_{1}$: (i) $(4.68,0.01)$ kHz and $0$, (ii) $(0.45,0.29)$ kHz and
$1.7$, (iii) $(0.79,0.51)$ kHz and $2.7$. d. Histogram of
$\tilde{z}_{1}/\lambda$ in the vicinity of the NIR region for
$(\Delta\Omega/2\pi,g_{a}/2\pi)=(0.56,0.51)$ kHz, corresponding to the
distribution (iv), exhibits two peaks. The solid black line is the bi-modal
Rice distribution fit, while the dashed and dash-dotted lines show the
contributions of two Rice distributions with fitted $\bar{A}_{1}=0$ and
$\bar{A}_{1}=3.1$, respectively.

We fit the Rice distribution to the histograms of the envelope
$|\tilde{z}_{j}|$ with $\bar{A}_{j}$ and $\sigma$ as free parameters. The
fitted amplitudes are then used to calculate the collective displaced
amplitude
$\bar{A}=2k\sqrt{\hbar/(m\Omega_{0})}\sqrt{\bar{A}_{1}^{2}+\bar{A}_{2}^{2}}$
and compared to the theoretical amplitude $A$ from Eq. (3). Fig. S6c shows the
histogram of the oscillation envelope for particle 1 fitted with a Rice
distribution outside the NIR ((i)) as well as within the NIR for increasing
interaction strength ((ii) and (iii)). In Fig. S6b we plot the corresponding
position-velocity distributions of particle 1. In the vicinity of the NIR
region, we observe a bi-modal distribution (Fig. S6d), which can be
interpreted as a time average of the collective state occasionally jumping
between the linear and nonlinear states. We attribute that to the thermal
fluctuations-induced drifts of the mechanical frequencies around their average
values. As a result, if the average frequencies are at the boundary of the
limit cycle region, the particles spend some time in the limit cycle phase
($\bar{A}_{j}>0$) and the rest in the correlated thermal state
($\bar{A}_{j}=0$) if otherwise. In these cases, we fit the histograms of the
oscillation envelopes with a weighted sum of two different Rice distributions:

$P_{j}(x)=\alpha P_{\rm Rice}(x,\sigma_{j,1},\bar{A}_{j,1})+(1-\alpha)P_{\rm
Rice}(x,\sigma_{j,2},\bar{A}_{j,2}).$ (S13)

For example, at a detuning of $\Delta\Omega/2\pi=(0.56\pm 0.15)$ kHz and a
coupling rate of $g_{a}/2\pi=(0.51\pm 0.03)$ kHz we obtain the amplitudes
$\bar{A}_{1,1}=0$ with the weight $\alpha=0.56$ and $\bar{A}_{1,2}=3.1$ with
the weight $1-\alpha=0.44$, thus confirming that our model explains the
average collective dynamics well. In the case of a bimodal distribution for
the individual oscillation envelope, we use the fitted displaced amplitude
with the highest weight for the calculation of the collective displaced
amplitude, shown in Fig. 4 of the main text.

###  A.7 Experimental evaluation of the correlations

The time traces of the complex envelopes $a_{1/2}$ were used to calculate the
first order correlations $g^{(1)}_{jj}/n_{T}=\langle
a^{*}_{j}(t)a_{j}(t+\Delta t)\rangle/n_{T}$ for $j\in\\{1,2\\}$, and
$g^{(1)}_{12}/n_{T}=\langle a^{*}_{2}(t)a_{1}(t+\Delta t)\rangle/n_{T}$ by
using python’s ”scipy.signal.correlate” function. In Fig.S7a, we show the
measured autocorrelation functions of particles 1 and 2 (left and middle
column) compared to the theoretical model (S3) (right column) for increasing
coupling rates $g_{a}$ from top to bottom. The autocorrelation is highest at
zero time delay, where it corresponds to the motion variance. In the
oscillatory phase, it exhibits decaying oscillations as a function of the time
delay with a decay rate of half of the damping rate, in excellent agreement
with the theory. The oscillation frequency depends on the intrinsic detuning,
therefore for uncoupled particles (upper row), the autocorrelations show an
oscillating decay everywhere except at zero detuning. The variances increase
with the coupling rate in the $\mathcal{PT}$-symmetry-broken phase and
eventually lead to the limit cycle phase, where the linearized theory is not
valid anymore (gray region in the theory plot, bottom row). Note that the
autocorrelations look very similar for both particles, confirming the
negligible contribution of the reciprocal coupling component ($g_{r}\approx
0$), according to the Eq. (S3). The crosscorrelation color plots as a function
of the mechanical detuning and the time delay are compared to the linearized
model (S4) for increasing coupling rates from top to bottom in Fig.S7b. The
magnitude of the crosscorrelation increases with the coupling rate for both
theory and experiment and is maximal at $\Delta t=0$, where it corresponds to
the amplitude covariance. Within the $\mathcal{PT}$ symmetry-broken regime,
the sign of the correlations changes due to the change in the relative phase
delay between the oscillations of the two particles.

Figure S7: First-order correlations: comparison between the experiment and
theory. Correlation functions are plotted as a function of the mechanical
detuning and the time delay for increasing coupling rate $g_{a}$ from top to
bottom. a. The autocorrelations in the left and middle columns correspond to
the experimental data from particles 1 and 2, respectively, and the right
column corresponds to the theoretical model. b. The experimental and
theoretical crosscorrelations are shown in the left and right columns,
respectively.


# Acoustic Rashba–Edelstein effect

Takumi Funato  Center for Spintronics Research Network, Keio University,
Yokohama 223-8522, Japan  Kavli Institute for Theoretical Sciences, University
of Chinese Academy of Sciences, Beijing, 100190, China  Mamoru Matsuo  Kavli
Institute for Theoretical Sciences, University of Chinese Academy of Sciences,
Beijing, 100190, China  CAS Center for Excellence in Topological Quantum
Computation, University of Chinese Academy of Sciences, Beijing 100190, China
RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama 351-0198, Japan
Advanced Science Research Center, Japan Atomic Energy Agency, Tokai, 319-1195,
Japan

###### Abstract

We theoretically study the mechanical induction of the spin density via the
Rashba spin–orbit interaction (SOI). The spin density in the linear response
to lattice distortion dynamics is calculated based on the microscopic theory.
We reveal that there are two mechanisms of spin induction: one is the acoustic
Edelstein effect (AEE) from the acceleration of the lattice dynamics and the
other is caused by the Rashba spin–vorticity coupling (RSVC). We find that the
AEE induces a more efficient spin-to-charge conversion in comparison with the
conventional electric Edelstein effect. The induced spin density due to the
RSVC is expressed as a Berry curvature-like quantity; therefore, it can be
attributed to the spatial symmetry breaking due to the Rashba SOI. Our work
demonstrates high-efficiency spin generation in Rashba systems.

The spin-orbit interaction (SOI), which represents coupling of orbital motion
and spin angular momentum, plays a significant role in the interconversion
between spin and charge in the field of spintronics. In particular, the Rashba
SOI, which arises from spatial inversion symmetry breaking, is expected to
provide a highly efficient spin–charge conversion. In Rashba systems, the spin
degeneracy of electronic energy bands is split, and the electron momentum
direction and spin direction are fixed in each band, a phenomenon known as the
spin–momentum lockingbychkov1984 ; rashba1 ; rashba2 ; ast2007 ; landolt2012 ;
eremeev2012 ; review2004 ; manchon2015 . When an external electric field is
applied, the energy bands shift, and a non-equilibrium spin density is
induced, which is known asthe Edelstein effect (EE)lyanda-geller1989 ;
pikus1989 ; edelstein1990 ; sanchez2013 ; fan2014 ; nomura2015 ; nakayama2016
; isasa2016 ; karube2016 . The direct/inverse EE provides high efficient
interconversion between charge and spin and is expected to substitute the
direct/inverse spin Hall effectd'yakonov1971 ; hirsch1999 ; zhang2000 ;
rashba2003 ; dimitrova2005 ; shi2006 ; saitoh2006 ; sinova2015 . Investigation
of the Rashba SOI is essential for the realization of efficient and versatile
spintronics devices.

Exploring different methods to generate a spin current is also an important
topic in the field of spintronics. Recently, the generation of spin currents
through the conversion from mechanical angular momentum to spin angular
momentum has attracted considerable attention. As one of its underlying
mechanisms, the spin–vorticity coupling (SVC), the coupling between the
electron spin and the effective magnetic field associated with the rigid
rotation, has been proposedbarnett ; ed ; matsuo1 ; matsuo2 ; matsuo3 ;
matsuo4 . Spin-current generation via the SVC has been observed in a shear
flow of liquid metals (spin hydrodynamic generation)svc_expt1 ; svc_expt2 ;
svc_expt3 ; svc_expt4 ; svc_theo1  and surface acoustic waveskobayashi2017 ;
kurimune1 ; kurimune2 ; tateno2020 ; tateno2 . Alternatively, mechanical spin-
current generation via the SOI has been investigated both theoretically and
experimentallyfunato2018 ; helicity ; kawada2021 . It is suggested that the
SOI mechanism can be dominant in strong SOI systems, and it provides spin
currens with a characteristic symmetry. Spin currents are expected to be
generated efficiently via mechanical means in Rashba systems; however, the
conversion of mechanical rotation into spin via the Rashba SOI still remains
an open problem.

In this work, we theoretically studied the induction of spin density from
lattice distortion dynamics in two dimensional (2D) and three-dimensional (3D)
Rashba electron systems. To treat a couple of lattice distortion dynamics and
the electrons, we performed the local coordinate transformation and calculated
the spin density in linear response to the lattice displacement. The results
indicate that there are two mechanisms through which the spin density
induction occurs: one is the acoustic Edelstein effect (AEE) due to the
lattice acceleration and the other is the Rashba spin–vorticity coupling
(RSVC) due to the lattice vorticity. The AEE was found to provide a more
efficient charge-to-spin conversion than the conventional EE. The spin density
due to the RSVC includes terms that are independent of the magnitude of the
Rashba SOI. This spin density can be attributed to the spatial symmetry
breaking due to the Rashba SOI and is expressed in terms of the Berry
curvature-like quantity. It is suggested that the spin density induced by the
RSVC is comparable to or larger than that induced by the conventional SVC.
Therefore, the combination of the RSVC and the conventional SVC provides a
more efficient spin generation in Rashba systems.

We model 2D and 3D electrons scattered by impurities in the crystal according
to the following equation:

$\displaystyle\hat{H}=\int d^{3}r^{\prime}$
$\displaystyle\hat{\psi}^{\prime\dagger}(\bm{r}^{\prime},t)\biggl{[}\frac{p^{2}}{2m}+V(\bm{r}^{\prime})$
$\displaystyle+\lambda_{\text{so}}\bm{\sigma}\cdot[\nabla^{\prime}V_{\text{p}}(\bm{r}^{\prime})\times\bm{p}\biggr{]}\hat{\psi}^{\prime}(\bm{r}^{\prime},t),$
(1)

where $\hat{\psi}^{\prime\dagger}$ and $\hat{\psi}^{\prime}$ are the electron
field operators, $V=V_{\text{p}}+V_{\text{i}}$ is the total potential which
combines the lattice periodic potential $V_{\text{p}}$ and the impurity
potential $V_{\text{i}}$, $\lambda_{\text{so}}$ is the SOI strength and
$\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ are the Pauli matrices.
Assuming that an acoustic wave is applied, the Hamiltonian is then modulated
according to the following equation:

$\displaystyle\hat{H}_{\text{lab}}=$ $\displaystyle\int
d^{3}r^{\prime}\hat{\psi}^{\prime\dagger}(\bm{r}^{\prime},t)\biggl{[}\frac{p^{2}}{2m}+V(\bm{r}^{\prime}-\delta\bm{r})$
$\displaystyle+\lambda_{\text{so}}\bm{\sigma}\cdot[\nabla^{\prime}V_{\text{p}}(\bm{r}^{\prime}-\delta\bm{r})\times\bm{p}\biggr{]}\hat{\psi}^{\prime}(\bm{r}^{\prime},t),$
(2)

where $\delta\bm{r}(\bm{r}^{\prime},t)$ is the lattice displacement vector.
The effects of the lattice distortion are taken into account by performing a
local coordinate transformation from the laboratory frame to the “material”
frame, i.e. $\bm{r}=\bm{r}^{\prime}-\delta\bm{r}(\bm{r}^{\prime},t)$; the
“material frame” is locally fixed to the moving latticefunato2018 ; helicity ;
tsuneto . Specifically, we consider the connection
$\nabla_{i}^{\prime}=(\partial r_{i}/\partial r^{\prime}_{j})\nabla_{j}$ and
the renormalization of the wave functions, i.e.
$\hat{\psi}^{\prime}(\bm{r}^{\prime},t)=\left|\partial\bm{r}^{\prime}/\partial\bm{r}\right|\hat{\psi}(\bm{r},t)$,
where $\left|\partial\bm{r}^{\prime}/\partial\bm{r}\right|$ is the Jacobian
and $\hat{\psi}$ is the wave function of the material frame. The Lagrangian is
given by the following equation:

$\displaystyle\hat{L}$
$\displaystyle=\int\psi^{\dagger}(\bm{r},t)i\partial_{t}\psi(\bm{r},t)d^{3}r-\hat{H}_{\text{mat}},$
(3)

where $\hat{H}_{\text{mat}}$ is the Hamiltonian given in the material frame.

Assuming a free-electron system with spatial symmetry breaking along
$z$-direction, the second quantized total Hamiltonian up to the first order in
$\delta\bm{r}$ is given by $\hat{H}_{\text{mat}}=\hat{H}+\hat{H}^{\prime}(t)$,
where the first term is the time-independent Hamiltonian, and the second term
is the perturbation Hamiltonian due to the acoustic waves. The time-
independent term $\hat{H}=\hat{H}_{0}+\hat{H}_{\text{i}}$ consists of the
unperturbed term (kinetic energy and Rashba SOI),
$\hat{H}_{0}=\sum_{\bm{k}}c_{\bm{k}}^{\dagger}\left[k^{2}/2m+\alpha_{\text{R}}\hat{z}\cdot(\bm{k}\times\bm{\sigma})\right]c_{\bm{k}}$,
and the impurity potential term,
$\hat{H}_{\text{i}}=\sum_{\bm{k},\bm{k}^{\prime}}V_{\bm{k}^{\prime}-\bm{k}}c_{\bm{k}^{\prime}}^{\dagger}c_{\bm{k}}$,
where $c^{\dagger}_{\bm{k}}(c_{\bm{k}})$ is the electron creation
(annihilation) operator, $\alpha_{\text{R}}$ is the strength of the Rashba
SOI, and $V_{\bm{k}}$ is the Fourier component of the impurity potential
$V_{\text{i}}(\bm{r})$. The perturbation Hamiltonian consists of two parts,
namely, $H^{\prime}(t)=H^{\prime}_{\text{k}}(t)+H^{\prime}_{\text{r}}(t)$,
where $H^{\prime}_{\text{k}}$ is the modulation of the kinetic energy and the
time derivativetsuneto ; funato2018 , and $H^{\prime}_{\text{r}}$ is the
modulation of the Rashba SOIhelicity :

$\displaystyle\hat{H}^{\prime}_{\text{k}}(\omega)$
$\displaystyle=-\sum_{\bm{k},\bm{q}}c^{\dagger}_{\bm{k}_{+}}w^{\text{k}}_{j}c_{\bm{k}_{-}}u^{j}_{\bm{q},\omega},$
(4) $\displaystyle\hat{H}^{\prime}_{\text{r}}(\omega)$
$\displaystyle=-\sum_{\bm{k},\bm{q}}c^{\dagger}_{\bm{k}_{+}}w^{\text{r}}_{j}c_{\bm{k}_{-}}u^{j}_{\bm{q},\omega},$
(5)

where $\bm{k}_{\pm}=\bm{k}\pm\frac{\bm{q}}{2}$ and $\bm{u}_{\bm{q},\omega}$
are the Fourier components of the lattice velocity,
$\bm{u}(\bm{r},t)=\partial_{t}\delta\bm{r}(\bm{r},t)$. $w^{\text{k}}_{j}$ and
$w^{\text{r}}_{j}$ are the vertex parts of the perturbation Hamiltonian:

$\displaystyle
w^{\text{k}}_{j}=\left(1-\frac{\bm{v}\cdot\bm{q}}{\omega}\right)k_{j},$ (6)
$\displaystyle
w^{\text{r}}_{j}=\frac{\alpha_{\text{R}}}{\omega}[\bm{k}\cdot(\bm{\sigma}\times\bm{q})\delta_{jz}-\hat{z}\cdot(\bm{\sigma}\times\bm{q})k_{j}].$
(7)

Throughout this study, we assume that the wavenumber $\bm{q}$ and the
frequency $\omega$ are much smaller than the inverse of the mean-free path $l$
and the relaxation time of the electrons, i.e. $q\ll l^{-1}$ and
$\omega\ll\tau^{-1}$, respectively.

Assuming a uniform random distribution of the impurities and a delta-function
like an impurity potential,
$V_{\text{i}}(\bm{r})=u_{\text{i}}\sum_{j}\delta(\bm{r}-\bm{R}_{j})$, where
$\bm{R}_{j}$ is the position of the $j$-th impurity, the impurity averaged
equilibrium retarded/advanced Green functions with impurity scattering
included by the Born approximation are given by the following equation:

$\displaystyle
G^{\text{R/A}}_{\bm{k}}(\epsilon)=\frac{1}{2}\sum_{s=\pm}\frac{1+s\Gamma_{\bm{k}}}{\mu+\epsilon-\epsilon_{s}\pm
i\gamma},$ (8)

where $\Gamma_{\bm{k}}=\epsilon_{z\beta l}\sigma^{\beta}k_{l}/k_{\text{t}}$
with the Levi-Civita symbol $\epsilon_{z\beta l}$,
$\epsilon_{s}=\frac{k^{2}}{2m}+s\alpha_{\text{R}}k_{\text{t}}$ is the
eigenenergy with $k_{\text{t}}=(k_{x}^{2}+k_{y}^{2})^{1/2}$ , $\mu$ is the
chemical potential, $\gamma=\frac{\pi}{2}n_{\text{i}}u_{\text{i}}^{2}\nu$ is
the damping rate with the impurity concentration $n_{\text{i}}$ and the Fermi-
level density of state $\nu=2a_{1}\nu_{0}$. The parameter $a_{1}$ and the
density of states for free-electrons $\nu_{0}$ are defined below Eq. (15).
Here, we use the ensemble average for impurity positions $\langle
V_{\bm{k}}V_{\bm{k}^{\prime}}\rangle_{\text{av}}=n_{\text{i}}u_{\text{i}}^{2}\delta_{\bm{k},\bm{k}^{\prime}}$.

The spin density and spin-current density operators are given by the following
equation:

$\displaystyle\hat{j}^{\alpha}_{\text{s},0}(\bm{q})=\hat{\sigma}^{\alpha}(\bm{q})=\sum_{\bm{k}}c^{\dagger}_{\bm{k}_{-}}\sigma^{\alpha}c_{\bm{k}_{+}},$
(9)
$\displaystyle\hat{j}^{\alpha}_{\text{s},i}(\bm{q})=\sum_{\bm{k}}c^{\dagger}_{\bm{k}_{-}}j^{\alpha}_{\text{s},i}(\bm{k})c_{\bm{k}_{+}},$
(10)

where the Greek indices $\alpha(=x,y,z)$ specify the spin direction and the
Roman indices $i(=x,y,z)$ specify the flow direction.
$j^{\alpha}_{\text{s},i}=-e(\sigma^{\alpha}v_{i}+\alpha_{\text{R}}\epsilon_{z\alpha
i})$ are the vertex parts of the spin-current density operators with
$v_{i}=k_{i}/m$. The current density operators are given by the following
equation:

$\displaystyle\hat{j}_{\text{e},i}(\bm{q})=\sum_{\bm{k}}c^{\dagger}_{\bm{k}_{-}}j_{\text{e},i}(\bm{k})c_{\bm{k}_{+}},$
(11)

where $j_{\text{e},i}(\bm{k})=-e(v_{i}+\alpha_{\text{R}}\epsilon_{z\alpha
i}\sigma^{\alpha})$ are the vertex parts of the current density operators.
Note that the contribution of the anomalous velocity due to
$\bm{u}_{\bm{q},\omega}$ is negligible in this study.

The expectation value in non-equilibrium states is given by the following
equation:

$\displaystyle\langle\hat{j}_{\text{s},\mu}^{\alpha}(\bm{q},\omega)\rangle_{\text{ne}}=\int^{\infty}_{-\infty}\frac{d\epsilon}{2\pi
i}\sum_{\bm{k}}\text{tr}\left[j^{\alpha}_{\text{s},\mu}G^{<}_{\bm{k}_{+},\bm{k}_{-}}(\epsilon_{+},\epsilon_{-})\right],$
(12)

where $\mu=(0,x,y,z)$ and $\epsilon_{\pm}=\epsilon\pm\frac{\omega}{2}$. Here,
the trace is taken for the spin space.
$G^{<}_{\bm{k},\bm{k}^{\prime}}(\epsilon,\epsilon^{\prime})$ represents the
lesser component of the non-equilibrium path-ordered Green function,
$G_{\bm{k},\bm{k}^{\prime}}(t,t^{\prime})=-i\langle
T_{\text{C}}c_{\bm{k}}(t)c_{\bm{k}^{\prime}}^{\dagger}(t^{\prime})\rangle_{\text{mat}}$,
where $T_{\text{C}}$ is a path-ordering operator, and
$\langle\cdots\rangle_{\text{mat}}$ represents the expectation value in the
non-equilibrium state

. We expand the non-equilibrium Green function for each of the perturbation
Hamiltonian $H^{\prime}_{\text{k}}(t)$ and $H^{\prime}_{\text{r}}(t)$ and
calculate the linear response to the lattice velocity field
$\bm{u}_{\bm{q},\omega}$ up to the first order in $\bm{q}$ and $\omega$.

First, we consider the so-called Fermi surface term, which can be dominant in
clean metals. The non-equilibrium spin density and the spin-current density
are given by

$\displaystyle\langle\hat{j}^{\alpha}_{\text{s},\mu}\rangle^{\text{k/r}}_{\text{surf}}=\frac{i\omega}{2\pi}\sum_{\bm{k}}\text{tr}\left[j^{\alpha}_{\text{s},\mu}G^{\text{R}}_{\bm{k}}\tilde{w}^{\text{k/r}}_{j}G^{\text{A}}_{\bm{k}}\right]u^{j}_{\bm{q},\omega},$
(13)

where $\tilde{w}^{\text{k/r}}_{j}=w^{\text{k/r}}_{j}+\Lambda^{\text{k/r}}_{j}$
with the three-point vertices $\Lambda^{\text{k/r}}_{j}$ presented in Fig. 1.

Figure 1:  (a) Linear response of the spin-current density ($\mu=x,y,z$) and
spin density ($\mu=0$) to the lattice velocity field $u^{j}_{\bm{q},\omega}$.
(b) Three-point vertices with ladder vertex corrections. The crosses and
dashed lines indicate the impurities and impurity potential, respectively.

The spin density and spin-current density are given by the following equation:

$\displaystyle\langle\hat{\sigma}^{\alpha}\rangle^{\text{k}}_{\text{surf}}=i\omega\alpha_{\text{R}}m\nu_{0}\tilde{a}_{2}\tau[\hat{z}\times\bm{u}_{\bm{q},\omega}]^{\alpha},$
(14)
$\displaystyle\langle\hat{j}^{z}_{\text{s},i}\rangle^{\text{k}}_{\text{surf}}=i\omega\frac{e\nu_{0}}{8}\tilde{a}_{2}\epsilon_{zij}u^{j}_{\bm{q},\omega},$
(15)

where
$a_{n}=\nu_{0}^{-1}\sum_{s,\bm{k}}(s\tilde{k}_{\text{t}})^{n-1}\delta(\mu-\epsilon_{(s)})$
are dimensionless parameters with $\tilde{k}=k/k_{\text{F}}$, and
$\nu_{0}=\frac{m}{2\pi}$ (for 2D), $\nu_{0}=\frac{mk_{\text{F}}}{2\pi^{2}}$
(for 3D) are the density of states in free-electron system without the Rashba
SOI with $k_{\text{F}}=\sqrt{2m\mu}$ being the Fermi wavenumber. For the
calculations, we also defined
$\tilde{a}_{2n}=-\tilde{\alpha}_{\text{R}}a_{2n}$. Note that a spin current
proportional to the first order in $\bm{q}$ is also generated, which can be
found in Ref. helicity . Eqs. (14) and (15) show that the spin density and
spin currents are generated by the acceleration of the lattice and are
considered to be the AEE and the spin Hall effect, respectively. Eq. (14) can
be put in the form of
$\langle\hat{\bm{\sigma}}\rangle=\lambda_{\text{A}}[\hat{z}\times\langle\hat{\bm{j}}_{\text{e}}\rangle]$
with the current density $\langle\hat{j}_{\text{e},i}\rangle$. Here,
$\lambda_{\text{A}}$ is the current–spin conversion efficiency:

$\displaystyle\lambda_{\text{A}}=-\tilde{\alpha}_{\text{R}}\frac{k_{\text{F}}}{e\mu}\left[\frac{a_{3}}{\tilde{a}_{2}}-2\tilde{\alpha}_{\text{R}}^{2}\left(1-\frac{\tilde{a}_{2}}{2a_{1}}\right)\right],$
(16)

where $\tilde{\alpha}_{\text{R}}=m\alpha_{\text{R}}/k_{\text{F}}$ is the
dimensionless Rashba parameter. The dependence of the conversion efficiency
$\lambda_{\text{A}}$ on the Rashba parameter $\tilde{\alpha}_{\text{R}}$ is
presented in Fig. 2 together with the case of the (electric) conventional EE,
$\lambda_{\text{e}}=-\tilde{\alpha}_{\text{R}}\frac{k_{\text{F}}}{e\mu}(1-\frac{a_{1}}{\tilde{a}_{2}})[\frac{a_{3}}{\tilde{a}_{2}}-2\tilde{\alpha}_{\text{R}}(2-\frac{\tilde{a}_{2}}{a_{1}}-\frac{a_{1}}{\tilde{a}_{2}})]$.
If Eq. (14) originates from the EE, $\lambda_{\text{A}}$ should match
$\lambda_{\text{e}}$ because EE is due to the band shift of the systems with
spin-momentum locking. However, $\lambda_{\text{A}}$ and $\lambda_{\text{e}}$
show different $\alpha_{\text{R}}$-dependencies, which indicates that the AEE
is a different mechanism from the conventional EE due to the electric field.
The current density generated by the lattice distortion is given by the
following equation:

$\displaystyle\langle\hat{j}_{\text{e},i}\rangle=-i\omega
e\mu\nu_{0}\tau\bigg{\\{}2c_{1}\delta^{z}_{ij}$
$\displaystyle+\biggl{[}a_{3}-2\tilde{\alpha}_{\text{R}}^{2}\tilde{a}_{2}\left(1-\frac{\tilde{a}_{2}}{2a_{1}}\right)\biggr{]}\delta^{2\text{d}}_{ij}\biggr{\\}}u_{\bm{q},\omega}^{j},$
(17)

where $\delta^{z}_{ij}=\delta_{iz}\delta_{jz}$,
$\delta^{\text{2d}}_{ij}=\delta_{ix}\delta_{jx}+\delta_{iy}\delta_{jy}$, and
$c_{n}=\nu_{0}\sum_{s,\bm{k}}\tilde{k}_{z}^{2}(s\tilde{k}_{\text{t}})^{n-1}\delta(\mu-\epsilon_{(s)})$.
The term independent of $\alpha_{\text{R}}$ in the current density is
consistent with the results obtained for the system without the Rashba SOI:
$\langle\hat{\bm{j}}_{\text{e}}\rangle=-i\omega
en_{\text{e}0}\tau\bm{u}_{\bm{q},\omega}$ for the $\alpha_{\text{R}}=0$ limit
with $n_{\text{e}0}=\frac{2k_{\text{F}}^{2}}{3m}\nu_{0}$ being the electron
number density for systems without the Rashba SOI.

Figure 2:  $\tilde{\alpha}_{\text{R}}$-dependence of the dimensionless
current–spin conversion efficiency,
$\tilde{\lambda}_{\text{A/e}}\equiv-\frac{e\mu}{k_{\text{F}}}\lambda_{\text{A/e}}$.
The red lines indicate the case in which the conversion occurs via mechanical
means and the blue lines indicate the (electric) conventional EE. The rigid
lines represent 3D systems and the dashed lines represent 2D systems. The
results show that the conversion efficiency due to mechanical means is
different from that due to the EE.

It is well-known that the spin Hall effect due to the uniform and static
external field vanishesrashba2003 ; dimitrova2005 ; shi2006 . Conversely, the
present spin Hall currents remain as presented in Eq. (15), since it is caused
by the lattice acceleration associated with the surface acoustic waves.

Next, we consider the so-called Fermi sea term, which contributes only to the
spin density response.

$\displaystyle\langle\hat{\sigma}^{\alpha}\rangle^{\text{r}}_{\text{sea}}=\frac{i\omega}{4\pi}\sum_{\bm{k}}\int^{0}_{-\infty}d\epsilon\text{tr}\biggl{[}\sigma^{\alpha}\partial_{\epsilon}G^{\text{R}}_{\bm{k}}(\epsilon)w^{\text{r}}_{j}G^{\text{R}}_{\bm{k}}(\epsilon)$
$\displaystyle-\sigma^{\alpha}G^{\text{R}}_{\bm{k}}(\epsilon)w^{\text{r}}_{j}\partial_{\epsilon}G^{\text{R}}_{\bm{k}}(\epsilon)-\langle\text{R}\leftrightarrow\text{A}\rangle\biggr{]}u^{j}_{\bm{q},\omega}.$
(18)

This term is calculated using the following equation:

$\displaystyle\langle\hat{\sigma}^{\alpha}\rangle^{\text{r}}_{\text{sea}}=-\frac{\nu}{4}\biggl{[}\left(1-\frac{\tilde{a}_{2}}{a_{1}}\right)\omega_{\bm{q},\omega}^{z}\delta_{\alpha
z}$
$\displaystyle-\left(\frac{15}{4}-\frac{7\tilde{a}_{2}}{8a_{1}}-\frac{\tilde{\alpha}_{\text{R}}^{2}\mu}{4\nu}\frac{d}{d\mu}\tilde{a}_{2}\nu_{0}\right)\omega_{\bm{q},\omega}^{\alpha}\delta_{\alpha\neq
z}\biggr{]},$ (19)

where $\bm{\omega}_{\bm{q},\omega}=i\bm{q}\times\bm{u}_{\bm{q},\omega}$ is the
vorticity of the lattice distortion. Eq. (19) suggests that the spin density
is induced by the vorticity motion of the lattice via the Rashba SOI, which
can be referred to as the RSVC. This includes the terms that do not depend on
the Rashba SOI, i.e.
$\langle\hat{\sigma}^{\alpha}\rangle^{\text{r}}_{\text{sea}}=\nu_{0}(\frac{1}{2}\omega^{z}_{\bm{q},\omega}+\omega_{\bm{q},\omega}^{\alpha}\delta_{\alpha\neq
z})$ for the non-Rahsba limit $\alpha_{\text{R}}\rightarrow 0$, however, the
RSVC requires the Rashba SOI. The dependence of the spin density induced via
the RSVC on the Rashba parameter is presented in Fig. 3 for the 3D system.
Here,
$\langle\hat{\sigma}^{\alpha}\rangle^{\text{tot}}=\langle\hat{\sigma}^{\alpha}\rangle^{\text{r}}+\langle\hat{\sigma}^{\alpha}\rangle^{\text{sv}}$
is the total spin density due to the RSVC and the conventional SVC. As the
RSVC does not occur in the absence of the Rashba SOI, the RSVC can be
attributed to the spatial symmetry breaking due to the Rashba SOI. The spin
density in response to the SVC
$H_{\text{sv}}=-\frac{1}{4}\bm{\sigma}\cdot\omega_{\bm{q},\omega}$ is given by
$\langle\hat{\sigma}^{\alpha}\rangle^{\text{sv}}=\frac{1}{4}\chi^{\text{s}}_{\alpha\beta}\omega^{\beta}$,
where
$\chi^{\text{s}}_{\alpha\beta}=\frac{\nu_{0}}{2}[\delta^{\text{2d}}_{\alpha\beta}\tilde{a}_{2}+2\delta^{z}_{\alpha\beta}(\tilde{a}_{2}-a_{1})]$
is the spin-spin correlation function. Note that the spin-spin correlation
function approach to the Pauli paramagnetism is
$\chi^{\text{s}}_{\alpha\beta}=2\nu_{0}\delta_{\alpha\beta}$ for
$\alpha_{\text{R}}\rightarrow 0$. In the 2D system, the spin density is
constant and independent of the Rashba parameter. As presented in Eq. 19, the
spin density induced by the RSVC is comparable to or larger than that caused
by the conventional SVC. Therefore, it is expected that the combination of the
RSVC and the conventional SVC provides more efficient spin generation in the
Rashba systems.

Figure 3:  $\tilde{\alpha}_{\text{R}}$-dependence of the dimensionless
current-spin conversion efficiency,
$\tilde{\lambda}_{\text{A/e}}\equiv\frac{\mu}{k_{\text{F}}}\lambda_{\text{A/e}}$.
The red lines represent the case of mechanical means and the blue lines
represent (electric) conventional EE. The rigid lines represent 3D systems and
the dashed lines represent 2D systems. The results show that the conversion
efficiency due to mechanical means is different from that due to EE.

Finally, we show that Eq. (19) can be expressed in terms of the Berry
curvature-like quantity according to the following equation:

$\displaystyle\langle\sigma^{\alpha}\rangle_{\text{sea}}^{\text{r}}=\frac{i\omega}{4}\sum_{\bm{k},s}f(\epsilon_{s})\Omega^{\alpha}_{s,j}(\bm{k})u^{j}_{\bm{q},\omega},$
(20)

where $f(\epsilon_{s})$ is the Fermi distribution function.
$\Omega^{\alpha}_{s,j}(\bm{k})$ is similar to the Berry curvature:

$\displaystyle\Omega^{\alpha}_{s,j}(\bm{k})=\frac{2\text{Im}[\sigma^{\alpha}_{s\bar{s}}w^{\text{r}}_{j,\bar{s}s}]}{(\epsilon_{s}-\epsilon_{\bar{s}})^{2}},$
(21)

where $\bar{s}=-s$. Here, $\sigma^{\alpha}_{s\bar{s}}$ and
$w^{\text{r}}_{j,s\bar{s}}$ are the matrix elements of the spin density and
the vertex part given in Eq. (7) for the eigenstates of $\epsilon_{s}$ and
$\epsilon_{\bar{s}}$, respectively.

In summary, we theoretically studied the generation of spin density from the
lattice distortion dynamics in 2D and 3D Rashba systems. The results indicate
that there are two mechanisms for the spin density induction: one is the AEE
due to the lattice acceleration and the other is the RSVC due to the lattice
vorticity. It is suggested that the AEE allows for more efficient spin-charge
conversion than the conventional EE. The spin density due to the RSVC remained
nonzero even in the non-Rashba limit $\alpha_{\text{R}}\rightarrow 0$. This
induced spin dnsity could be attributed to the spatial symmetry breaking due
to the Rashba SOI and was expressed in terms of the Berry curvature-like
quantity. The spin density induced by the RSVC is suggested to be comparable
to or larger than that induced by the conventional SVC. Since the RSVC and the
conventional SVC are cooperative, a more efficient spin generation is
envisaged in the Rashba systems.

We would like to thank Y. Nozaki, K. Yamanoi, and T. Horaguchi for
enlightening discussions. We also thank H. Kohno, A. Yamakage, Y. Imai, J. J.
Nakane. Y. Yamazaki, Y. Ogawa, Y. Ozu, and Y. Hayakawa for daily discussions.
This work was partially supported by JST CREST Grant Number JPMJCR19J4, Japan.
TF is supported by Grant-in Aid for JSPS Fellows Grant Number 19J15369, and by
a Program for Leading Graduate Schools ”Integrative Graduate Education and
Research in Green Natural Sciences”. MM is supported by JSPS KAKENHI for
Grants (Nos. 20H01863 and 21H04565) and the Priority Program of the Chinese
Academy of Sciences, Grant No. XDB28000000.

## References

  * (1) E. I. Rashba and V. I. Sheka, Fiz. Tverd. Tela: Collected Paper 2, 162 (1959).
  * (2) E. I. Rashba, Sov. Phys. Solid State 2, 1109 (1960).
  * (3) A. G. Aronov and Y. B. Lyanda-Geller, J. Phys. C: Solid State Phys., 17, 6039 (1984).
  * (4) C. R. Ast, J. Henk, A. Ernst, L. Moreschini, M. C. Falub, D. Pacilé, P. Bruno, K. Kern, and M. Grioni, Phys. Rev. Lett. 98, 186807 (2007).
  * (5) G. Landolt, S. V. Eremeev, Y. M. Koroteev, B. Slomski, S. Muff, T. Neupert, M. Kobayashi, V. N. Strocov, T. Schmitt, Z. S. Aliev, Phys. Rev. Lett. 109, 116403 (2012).
  * (6) S. V. Eremeev, I. A. Nechaev, Yu. M. Koroteev, P. M. Echenique, and E. V. Chulkov, Phys. Rev. Lett. 108, 246802 (2012).
  * (7) I. Žutić, J. Fabian, and S. D. Sarma, Rev. Mod. Phys. 76, 323 (2004).
  * (8) A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A. Duine, Nature Materials 14, 871 (2015).
  * (9) A. G. Aronov and Y. B. Lyanda-Geller, Pis’ma Zh. Eksp. Teor. Fiz. 50, 398 (1989) [JETP Lett. 50, 431 (1989)].
  * (10) E. L. Ivchenko, Y. B. Lyanda-Geller, and G. E. Pikus, Pis’ma Zh. Eksp. Teor. Fiz. 50, 156 (1989) [JETP Lett. 50, 175 (1989)].
  * (11) V. M. Edelstein, Solid State Commun. 73, 233 (1990).
  * (12) J. C. R. Sánchez, L. Vila1, G. Desfonds, S. Gambarelli, J. P. Attane, J. M. D. Teresa, C. Magén, and A. Fert, Nat. Commun. 4, 2944 (2013).
  * (13) Y. Fan, P. Upadhyaya, X. Kou, M. Lang, S. Takei, Z. Wang, J. Tang, L. He, L.-T. Chang, M. Montazeri, G. Yu, W. Jiang, T. Nie, R. N. Schwartz, Y. Tserkovnyak. and K. L. Wang, Nat. Mat. 13, 699 (2014).
  * (14) A. Nomura, T. Tashiro, H. Nakayama, and K. Ando, Appl. Phys. Lett. 106, 212403 (2015).
  * (15) H. Nakayama, Y. Kanno, H. An, T. Tashiro, S. Haku, A. Nomura, and K. Ando, Phys. Rev. Lett. 117, 116602 (2016).
  * (16) M. Isasa, M. C. M.-Velarte, E. Villamor, C. Magén, L. Morellón, J. M. D. Teresa, M. R. Ibarra, G. Vignale, E. V. Chulkov, E. E. Krasovskii, L. E. Hueso, and F. Casanova, Phys. Rev. B 93, 014420 (2016).
  * (17) S. Karube, K. Kondou, and Y. Otani, Appl. Phys. Exp. 9, 033001 (2016).
  * (18) M. I. D’yakonov and V. I. Perel’, JETP Lett. 13, 467 (1971).
  * (19) J. E. Hirsch, Phys. Rev. Lett. 83, 1834 (1999).
  * (20) S. Zhang, Phys. Rev. Lett. 85, 393 (2000).
  * (21) E. I. Rashba, Phys. Rev. B 68, 241315(R) (2003).
  * (22) O. V. Dimitrova, Phys. Rev. B 71, 245327 (2005).
  * (23) J. Shi, P. Zhang, D. Xiao, and Q. Niu, Phys. Rev. Lett. 96, 076604 (2006).
  * (24) E. Saitoh, M. Ueda, H. Miyajima, and G. Tatara, Appl. Phys. Lett. 88, 182509 (2006).
  * (25) J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, Rev. Mod. Phys. 87, 1213 (2015).
  * (26) S. J. Barnett, Phys. Rev. 6, 239 (1915).
  * (27) A. Einstein and W. J. de Haas, Verh. Dtsch. Phys. Ges. 17, 152 (1915).
  * (28) M. Matsuo, J. Ieda, E. Saitoh, and S. Maekawa, Phys. Rev. Lett. 106, 076601 (2011).
  * (29) M. Matsuo, J. Ieda, E. Saitoh, and S. Maekawa, Phys. Rev. B 84, 104410 (2011).
  * (30) M. Matsuo, J. Ieda, and S. Maekawa, Phys. Rev. B 87, 115301 (2013).
  * (31) M. Matsuo, J. Ieda, K. Harii, E. Saitoh, and S. Maekawa, Phys. Rev. B 87, 180402(R) (2013).
  * (32) R. Takahashi, M. Matsuo, M. Ono, K. Harii, H. Chudo, S. Okayasu, J. Ieda, S. Takahashi, S. Maekawa, and E. Saitoh, Nat. Phys. 12, 52 (2016).
  * (33) R. Takahashi, H. Chudo, M. Matsuo, K. Harii, Y. Ohnuma, S. Maekawa, and E. Saitoh, Nat. Commun. 11, 3009 (2020).
  * (34) H. T. Kazerooni, A. Thieme, J. Schumacher, and C. Cierpka, Phys. Rev. Appl. 14, 014002 (2020).
  * (35) H. T. Kazerooni, G. Zinchenko, J. Schumacher, and C. Cierpka, Phys. Rev. Fluids 6, 043703 (2021).
  * (36) M. Matsuo, Y. Ohnuma, and S. Maekawa, Phys. Rev. B 96, 020401(R) (2017).
  * (37) D. Kobayashi, T. Yoshikawa, M. Matsuo, R. Iguchi, S. Maekawa, E. Saitoh, and Y. Nozaki, Phys. Rev. Lett. 119, 077202 (2017).
  * (38) Y. Kurimune, M. Matsuo, S. Maekawa, and Y. Nozaki, Phys. Rev. B 102, 174413 (2020).
  * (39) Y. Kurimune, M. Matsuo, and Y. Nozaki, Phys. Rev. Lett. 124, 217205 (2020).
  * (40) S. Tateno, G. Okano, M. Matsuo, and Y. Nozaki, Phys. Rev. B 102, 104406 (2020).
  * (41) S. Tateno, Y. Kurimune, M. Matsuo, K. Yamanoi, and Y. Nozaki, Phys. Rev. B, to appear.
  * (42) T. Funato and H. Kohno, J. Phys. Soc. Jpn. 87, 073706 (2018).
  * (43) T. Funato and M. Matsuo, arXiv:2105.04804 (2021).
  * (44) T. Kawada, T. Funato, H. Kohno, and M. Hayashi, Sci. Adv. 7, eabd9697 (2021).
  * (45) T. Tsuneto, Phys. Rev. 121, 402 (1961).