# Entanglement transition in rod packings

Yeonsu Jung1, Thomas Plumb-Reyes1, Hao-Yu Greg Lin2,
and L. Mahadevan1,3,4∗
1John A. Paulson School of Engineering and Applied Sciences,
2Center for Nanoscale Systems,
3Department of Physics,
4Department of Organismic and Evolutionary Biology, Harvard University,
Cambridge, MA 02138, USA
∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> Random packings of stiff rods are self-supporting mechanical structures
> stabilized by their geometrical and topological complexity. To understand
> why, we deploy X-ray computerized tomography to unveil the structure of the
> packing. This allows us to define and directly visualize the spatial
> variations in “entanglement,” a mesoscopic field that characterizes the
> local average crossing number, a measure of the topological complexity of
> the packing. We show that the entanglement field has information that is
> distinct from the density, orientational order, and contact distribution of
> the packing. We find that increasing the aspect ratio of the constituent
> rods in a packing leads to an abrupt change in the entanglement, correlated
> with a sharp transition in the mechanical response of the packing. This
> leads to an entanglement phase diagram for the mechanical response of dense
> rod packings that is likely relevant for a broad range of problems that goes
> beyond our specific study.
>
> One-sentence summary: Sufficiently long and thin rods entangle with each
> other to form mechanically stable structures.

##  1 Introduction

The structure and response of polydisperse granular systems consisting of
localized particles, such as spheres and ellipsoids, can be characterized via
their number density and positional order. At sufficiently large densities,
there is a transition to a disordered jammed packing with emergent mechanical
properties that are controlled by a topological order parameter, the average
coordination number (?, ?, ?, ?, ?). However, experimentally, spheroidal or
ellipsoidal packings need to be confined by boundaries to see these effects;
the packings are mechanically unstable in the absence of confining
constraints. This is in sharp contrast with the well known behavior of
disordered packings of long, thin rods seen in the familiar example of a bird
nest or a randomly packed felt that forms solid-like states by jamming and
densities that are much smaller than those associated with jammed spheres (?,
?, ?, ?).

The extended nature of individual filaments that can entangle naturally
implies the presence of topological constraints associated with contact in
addition to positional and orientational order/disorder. While these have been
well studied in thermalized systems exemplified by polymer melts (?, ?, ?),
their athermal macroscopic analogs seen in textiles (?), felts (?) etc. are
much less understood, although recent work has started to highlight the
importance of topological constraints (?, ?, ?) in rod packings of large
aspect ratio filaments (?, ?, ?) or curved filaments (?, ?) with applications
to soft robotics (?) and active systems (?, ?).

To understand the stability of dense rod packings, here we use X-ray
computerized tomography (X-ray CT) to examine the link between topology,
geometry and mechanics of strongly entangled structures.

##  2 Results

We consider a packing of steel rods with varying length-to-diameter ratios,
i.e. aspect ratios $\alpha\in[38,200]$, (see Table S1), imaged using X-ray CT
and processed to accurately locate the centerline of the rods and their
contacts (see Materials and Methods). This allows us to quantify the local
volume fraction, orientational order, contact density, and various measures of
entanglement (see Materials and Methods) and correlate them with the
mechanical stability of the packing.

###  2.1 Stability of topological cohesion

We prepared random packings of rods by pouring them gently into a cylindrical
container through a mesh. Figure 1A and Video S1 show that such a packing of
rods with $\alpha=200$ can resist its own weight while a packing of rods with
$\alpha=38$ leads to an untangled pile that cannot sustain its own weight. We
note that the presence of gravity breaks vertical symmetry, i.e. rods are
likely to be heterogeneously arranged in this direction, confirmed by our
later analysis of the packing geometry.

To understand the stability of these packings, we subject them to harmonic
excitation with varying frequencies and amplitudes. In Fig. 1B and C, we see
that a high aspect-ratio rod packing ($\alpha=200$) forms a dynamically stable
clump, while a low aspect-ratio rod packing ($\alpha=38$) collapses (data
shown for frequency, $f=5$ Hz, and amplitude, $A=25$ mm, so that the
dimensionless maximum acceleration $4\pi^{2}f^{2}A/g=2.5$ where $g=9.81$ m/s2,
although this behavior persists for a range of frequencies and amplitudes).

A method to quantify the stability of the entangled packings is afforded by
the height of the clump as a function of time. In Fig. 1D, we see that the
mechanical excitation can either destroy cohesion (for packings with
$\alpha\leq 100$), or enhance it (for packings with $\alpha\geq 100$). To
compare the dynamic stability of rod packings as a function of the aspect
ratio $\alpha$, we prepare samples that have a fixed average contact number
per rod ($\approx 6$ for randomly oriented rod packings) (?) with the volume
of minimum convex hull of the clump being the same. The response to harmonic
excitation of these packings allows us to define a scaled time
$\tau=t_{\mathrm{entangled}}/100~{}\mathrm{sec}$, the ratio of the duration
for which the entangled phase persists relative to the observation time. In
Fig. 1E, we show that $\tau\approx 0$ for $\alpha\ll 100$ but increases as
$\alpha\geq 100$, a dynamical signature of an entanglement transition.

###  2.2 Mesoscopic measures of volume fraction, orientational order, average
contact number and entanglement

X-ray CT and image analysis allow us to characterize the volume fraction
$\phi$, orientational order $S$, contact density $c$, and entanglement $e$
(see Fig. 2A and Materials and Methods) as a function of location and applied
loads. For high aspect packings which possess two well-separated length
scales, rod diameter $d$ and length $l$, a natural mesoscopic scale is the
geometric mean, $(dl)^{1/2}$ (see Materials and Methods for details). This
mesoscopic scale provides a natural scale over which to visualize the spatial
variations of $\phi(\mathbf{x})$, $S(\mathbf{x})$, $c(\mathbf{x})$, and
$e(\mathbf{x})$ defined in Fig. 2A (see Materials and Methods as well).

The spatially integrated probability densities $p(\phi)$, $p(S)$, $p(c)$, and
$p(e)$ for packings of rods with $\alpha=38$, $\alpha=100$ are shown in Fig.
2B-E and Fig. 2F-I, respectively, for different values of applied mechanical
strain $\epsilon\in[0,0.15]$, while the spatial distributions of these
mesoscopic fields is shown in Fig. 2K-N for the same values of $\alpha$ and
$\epsilon$.

The volume fraction of rods in a random packing are controlled by the aspect
ratio $\alpha$ and the average number of contacts per rod $Z$, via the
relation $\phi\alpha\sim Z$ (?). Therefore, $\langle\phi\rangle$ decreases as
$\alpha$ becomes larger, as shown in Fig. 2A,B where
$\langle\cdot\rangle=(1/V_{\mathrm{total}})\int(\cdot)\mathrm{d}^{3}\mathbf{x}$.

The distribution of the orientational order parameter $S$ is not completely
isotropic ($S=0$). Instead, they tend to be slightly aligned ($S>0$) with or
perpendicular to the local director ($-1/2<S<0$), although these variations
average out globally with $S_{\mathrm{total}}\approx 0$, as shown in Fig. 2C
and G. For long slender rod packings ($\alpha=100$), the positive side ($S>0$)
of $p(S)$ becomes larger than the negative side ($-1/2<S<0$), compared to
shorter rods ($\alpha=38$), which implies that the packing of rods with
$\alpha=100$ is more aligned than that with $\alpha=38$.

The contact density $p(c)$ and entanglement distribution $p(e)$ behave
similarly as a function of $\alpha$. For long slender rod packing
($\alpha=100$), $p(c)$ and $p(e)$ become fat-tailed distributions compared to
the distributions for shorter rods ($\alpha=38$), although $p(e)$ exhibits a
longer and fatter tail than $p(c)$ does. This implies that the entanglement
$e$ is more sensitive than contact $c$ to variations in the positional and
orientational order of the rods.

To reveal the spatial correlations between rod position, orientation, contact
density and entanglement, we visualize their 2D averaged displays in Fig. 2J-M
or 3D renderings (Video S4). We see that packings of rods with $\alpha=38$
show short-range fluctuations in the fields
$\phi(\mathbf{x}),S(\mathbf{x}),c(\mathbf{x}),~{}e(\mathbf{x})$. Additionally,
$c$ and $e$ show signs of heterogeneity in the bulk that increase with an
increase in the applied mechanical strain, but $S$ shows evidence of
heterogeneity primarily due to the ordering influence of boundaries.

In contrast, for rods with $\alpha=100$, we see the appearance of much more
strongly heterogeneous variations in the mesoscopic geometric fields. For
example, the orientational order $S$ shows the presence of a band of aligned
rods, $S>0.5$ near the lateral walls (see the arrows in Fig. 2K), a feature
not seen for rods with $\alpha=38$. Similarly, there are multiple hot spots of
$c$ and $e$ shown in Fig. 2L and M. These hot spots increase in number and
relative intensity in response to mechanical strain $\epsilon$, consistent
with an enhanced mechanical and dynamical stability described in Fig. 1 (see
Section S2.1, Fig. S5-S8, which shows the gradual change of the spatially
varying fields, $\phi(\mathbf{x})$, $S(\mathbf{x})$, $c(\mathbf{x})$, and
$e(\mathbf{x})$ for $\alpha=(38,66,76,100,200)$).

To provide a comparative baseline for the nature of rod packings in our
physical experiments, we generated numerical simulations of rod packings by
adding rods of finite diameter at random positions and with random orientation
into a packing (as a numerical array of rod coordinates), accepting it only if
it is non-intersecting either with every other rod or the domain boundary (see
Section S2.2). In Fig. S9, we see that our numerically generated rod packings
do not show the experimentally observed variability of the mesosopic geometric
fields as a function of the rod aspect ratio $\alpha$ (see Video S5). To
understand this, we note that in experiments, rods enter a packing through the
boundary of a packing primarily via random excitation, so that topological
constraints associated with an increase in $\alpha$ have a strong effect on
the formation of entanglement tails in real rod packings shown in Fig. S9B. In
contrast, in numerically generated rod packings, the chosen protocol of random
insertion of rods is assumed to be uniform within a packing, but subject to
steric constraints. This makes the protocol much less dependent on the rod
aspect ratio, as seen in Fig. S9C, explaining the difference between the
numerical and experimental rod packings; however there are likely other
protocols where this is not necessarily true.

###  2.3 Analysis of vertical heterogeneity in rod packings

Gravity affects our rod packings owing to the sieving protocol that we use,
and we expect this to be visible in the vertical heterogeneity in the rod
packings. Indeed, in Fig. 2K, we see that a band of local orientational
ordering is strongly correlated with entanglement seen in Fig. 2M. By
averaging out the radial and azimuthal dependencies of the mesoscopic fields
(which are relatively homogeneous in these directions away from the
boundaries), we can evaluate $\phi(z),S(z),c(z)$ and $e(z)$ and extract the
vertical variation in them as a function of the height $z$. We then visualize
each of these fields over a thickness $h$, corresponding to the slice between
$z-h/2$ and $z+h/2$, see Materials and Methods.

In Fig. 3A, for elongated rods $\alpha\geq 66$, we see the appearance of a
strong peak in the normalized volume fraction $\tilde{\phi}$ (see definitions
in Materials and Methods) in the bulk of the packing (see arrows). Consistent
with this, the rods in this crowded layer are also aligned, see Fig. 3B.
However, the contact density $c(z)$ in Fig. 3C remains relatively homogeneous
away from the boundary, i.e. the rods are not overly constrained in these
oriented, dense regions. Most interestingly, the hot spots of $e(\mathbf{x})$
the 2D averaged images shown in Fig. 2M are now shown to be correlated with
the strong vertical heterogeneity in $e(z)$ which is largest in the zone where
$S(z)$ and $\phi(z)$ are large (see Fig. 3A,B, and D). Our observations
clearly indicate how our packing protocol (where gravity and steric effects
associated with the insertion of rods into a packing are relevant) affects
vertical heterogeneities in entanglement for $\alpha\geq 100$, even though the
contact density remains relatively homogeneous.

###  2.4 Entanglement phase diagram

Both mechanical testing and image analysis indicate that rod aspect ratio
controls the change in the mechanical response of the packing via changes in
the geometry and topology of the network. To further probe these aspects in
light of our packing protocol that relies on gravity, we calculate the
variations in entanglement as a function of depth in the packing, focusing on
the correlations in $e$ as a function of the depths $z$ in the packing. Using
a normalized version of Eq. S6, with
$\tilde{e}(z,z^{\prime})=e(z,z^{\prime})/e_{\mathrm{total}}$, where
$e_{\mathrm{total}}$ is the total entanglement of a rod packing, in Fig. 3E-H,
we show that the entanglement between two distinct layers can be displayed
either using matrix plots (E and F) or in the 3D rendering of the packing
color-coded with $\tilde{e}(z)=\sum_{z^{\prime}}\tilde{e}(z,z^{\prime})$ (G
and H).

For a low aspect ratio packing ($\alpha=38$), non-zero diagonal terms indicate
the presence of short-ranged entangling interactions between layers. However,
when $\alpha>100$, a clear entanglement band appears around the aligned layer
which mediates long-ranged entangling interactions throughout the entire
domain. Comparison with the entanglement between layers in simulated packings
confirm that this entangled band arises due to topological constraints as
shown in Fig. S10; this band persists under static and dynamic loading and is
responsible for the mechanical stability of the packing.

We note that the average entanglement for rod packings according to our
protocol is insensitive to change in the aspect ratio as shown in Fig. S9D. On
the other hand, higher moments of $p(e)$ is strongly and positively correlated
with the increasing aspect ratio (see Fig. S9). To quantify the heterogeneity
in $p(e)$, we compute the inverse of the coefficient of variation, i.e. the
ratio of the mean to the standard deviation of the packing $\sigma/\mu$ (see
Fig. S9E for a comparison across different rod packings). As expected, the
coefficient of variation is monotonically increasing with $\alpha$ as shown in
Fig. S9F.

This suggests a phase diagram (Fig. 3I) for the entanglement of extended
particles (i.e. rods) in terms of three natural variables, the entanglement
parameter, $\mu/\sigma$, and scaled measures of static and dynamic loading,
$F/\rho_{s}gd^{2}l$ and $a/g$ respectively, as shown in Fig. 3I, allowing us
to see a direct analogy to the classical jamming phase diagram for localized
particles (?), where the reciprocal of volume fraction ($1/\phi$), temperature
($T$), and the applied stress ($\Sigma$) are similar to the entanglement, the
dynamic load, and the static load, respectively. The pink shaded region in
Fig. 3I corresponds to the case of rod packings with $\alpha<50$ where the
experimentally observed distribution of $p(e)$ approaches that of the ideal
random packing of rods (see Fig. S11, along with Section S2.3), where stable
entanglement is impossible. As the entanglement parameter $\mu/\sigma$ varies
(due to a change in the aspect ratio of the rods), we expect to see a
transition between a stable entangled phase (green shade in Fig. 3I) and and
an untangled (white region in Fig. 3I) phase, suggesting its use as a natural
order parameter.

##  3 Summary

Inspired by the robust mechanical stability of packings of long stiff rods, we
used X-ray CT imaging to investigate and quantify a number of different
mesoscopic geometric fields that characterize the geometry and topology of
these systems. Our results show the emergence of a stable packing when the rod
aspect ratio is higher than a threshold ($\approx 100$), linking mechanical
rigidity and a mesocopic entanglement field that characterizes the geometry
and topology of the packing. Our results show that the coefficient of
variation of the entanglement probability distribution $p(e)$ increases
monotonically with the aspect ratio in the physical packings, while it does
not in ideal packings (?). This suggests its use as an order parameter to
define an entanglement transition, which we confirm using our experimental
data: the entanglement order parameter, $\mu/\sigma\approx 1$, is of order
unity for a packing of rods with $\alpha=100$. Understanding how the critical
value of this order parameter depends on the stiffness of the rods and the
physics of their interaction through friction, lubrication etc. awaits further
study, which we hope will be inspired by the first steps outlined here.

## Acknowledgments

Funding. This research was supported by National Research Foundation of Korea
(Grant No. 2021R1A6A3A03039239), and the Harvard NSF MRSEC 20-11754, NSF DMREF
19-22321, the Simons Foundation and the Henri Seydoux Fund, and performed in
part at the Harvard University Center for Nanoscale Systems (CNS); a member of
the National Nanotechnology Coordinated Infrastructure Network (NNCI), which
is supported by the National Science Foundation under NSF award no.
ECCS-2025158.

Author contributions. L.M. conceived of the question and approaches. Y.J.,
T.P-R. and L.M. designed the study. Y.J. and T.P-R. carried out the
experiments with the help of G.L. Y.J. analyzed the data. Y.J. and L.M. wrote
the paper. L.M. supervised the study. The authors declare no conflict of
interest.

Competing interests. The authors declare no competing interests.

Data availability. All data are available in the main text or the
supplementary materials.

## Supplementary materials

Materials and Methods
Figs. S1 to S11
Table S1
Sections S1 to S3
Videos S1 to S6
References (26 - 35)

## References

  * 1. A. Baule, F. Morone, H. J. Herrmann, H. A. Makse, Rev. Mod. Phys. 90, 15006 (2018).
  * 2. A. J. Liu, S. R. Nagel, Annual Review of Condensed Matter Physics 1, 347 (2010).
  * 3. A. J. Liu, S. R. Nagel, eds., Jamming and Rheology: Constrained Dynamics on Microscopic and Macroscopic Scales (CRC Press, London, 2014).
  * 4. C. S. O’Hern, S. A. Langer, A. J. Liu, S. R. Nagel, Physical Review Letters 88, 075507 (2002).
  * 5. C. S. O’Hern, L. E. Silbert, A. J. Liu, S. R. Nagel, Physical Review E 68, 011306 (2003).
  * 6. A. Ekman, A. Miettinen, T. Tallinen, J. Timonen, Phys. Rev. Lett. 113, 1 (2014).
  * 7. A. P. Philipse, Langmuir 12, 5971 (1996).
  * 8. N. Weiner, Y. Bhosale, M. Gazzola, H. King, Journal of Applied Physics 127, 050902 (2020).
  * 9. Y. Bhosale, et al., Phys. Rev. Lett. 128, 198003 (2022).
  * 10. L. Onsager, Annals of the New York Academy of Sciences 51, 627 (1949).
  * 11. S. F. Edwards, Proceedings of the Physical Society 91, 513 (1967).
  * 12. M. Doi, S. F. Edwards, J. Chem. Soc. Faraday Trans. 2 Mol. Chem. Phys. 74, 560 (1978).
  * 13. J. W. S. Hearle, W. E. Morton, Physical Properties of Textile Fibres (Elsevier Science, 2008).
  * 14. A. Kabla, L. Mahadevan, Journal of The Royal Society Interface 4, 99 (2007).
  * 15. G. Buck, J. Simon, Proc. R. Soc. A Math. Phys. Eng. Sci. 468, 4024 (2012).
  * 16. E. Panagiotou, L. H. Kauffman, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476, 20200124 (2020).
  * 17. V. P. Patil, J. D. Sandt, M. Kolle, J. Dunkel, Science 367, 71 (2020).
  * 18. M. Trepanier, S. V. Franklin, Physical Review E 82, 011308 (2010).
  * 19. S. V. Franklin, Physics Today 65, 70 (2012).
  * 20. C. Heussinger, Collapse of columns of granular rods (2023). ArXiv:2301.05596 [cond-mat].
  * 21. N. Gravish, S. V. Franklin, D. L. Hu, D. I. Goldman, Phys. Rev. Lett. 108, 208001 (2012).
  * 22. Y. Zhao, et al., Granular Matter 18, 24 (2016).
  * 23. K. Becker, et al., Proceedings of the National Academy of Sciences 119, e2209819119 (2022).
  * 24. V. P. Patil, et al., Science 380, 392 (2023).
  * 25. G. O. Bozdag, et al., Nature 617, 747 (2023).
  * 26. L. R. Gómez, N. A. García, T. Pöschel, Proceedings of the National Academy of Sciences 117, 3382 (2020).
  * 27. A. Rigort, et al., Journal of Structural Biology 177, 135 (2012).
  * 28. J. C. Lagarias, J. A. Reeds, M. H. Wright, P. E. Wright, SIAM Journal on Optimization 9, 112 (1998).
  * 29. J. Blouwolff, S. Fraden, Europhys. Lett. 76, 1095 (2006).
  * 30. C. Heussinger, Phys. Rev. E 103, 052903 (2021).
  * 31. M. A. Berger, C. Prior, J. Phys. A. Math. Gen. 39, 8321 (2006).
  * 32. Y. Diao, A. Dobay, R. B. Kusner, K. Millett, A. Stasiak, J. Phys. A. Math. Gen. 36, 11561 (2003).
  * 33. D. Baum, et al., Integrative and Comparative Biology 59, 1700 (2019).
  * 34. K. Klenin, J. Langowski, Biopolymers 54, 307 (2000).
  * 35. T. Banchoff, Indiana University Mathematics Journal 25, 1171 (1976).

Figure 1: Formation and stability of entangled rods. (A) Entangled steel rods
can support their own weight (See Video S1 for details). (B,C) Vertical
excitation with a fixed frequency (5 Hz) and amplitude (25 mm) of rod packing
with aspect ratios, $\alpha=38$ (B) and $\alpha=333$ (C). Inset shows that the
large aspect ratio rod packing survives, but the small aspect ratio packing
breaks down. (D) Normalized centroid height of rods,
$y_{\mathrm{c}}(t)/y_{\mathrm{c}}$, as a function of number of excitation
cycles, $N=t\times f$ with time $t$ and excitation frequency $f$, of packings
with different aspect ratios ranging from 38 to 200. Formation (i) and
destruction (ii) of an entangled unit during the excitation are denoted near
the onset of them for $\alpha=125$. (E) The ratio of entangled period to the
observation time, $\tau=t_{\mathrm{entangled}}/T$, ($T=100$ sec). $\tau$
serves as a measure of the persistence of the entangled phase under
excitation.

Figure 2: Experimental characterization of geometry and topology of a dense
rod packing. (A) Reconstructed 3D rendering of X-ray CT scan of a local
domain. For each rod $i$, we identify the centerline $\mathbf{r}_{i}(s_{i})$
and unit tangent vector
$\mathbf{t}_{i}=\dot{\mathbf{r}}_{i}/|\dot{\mathbf{r}_{i}}|$. Volume fraction
($\phi$), orientational order ($S$), contact density ($c$), and entanglement
($e$) fields can be calculated using the Eqs. S1-S5 where A(i) - A(iv) show
corresponding graphical description of those fields. (B-I) Statistics of local
volume fraction $\phi$ (B and F), orientational order $S$ (C and G), contact
density, $c$ (D and H), and entanglement $e$ in Eqs. S1-S5 (E and I) fields of
the rod packing when $\epsilon=0$ and $\epsilon=0.15$. (J-M) Visualizations of
horizontal- and vertical-averaging of $\phi$ (J), $S$ (K), $c$ (L), and $e$
(M) fields. The first two rows correspond to horizontal averaging and the last
two rows correspond to vertical averaging. For both vertical- and horizontal-
averaged visualizations include results for two different strain values,
$\epsilon=0$ and $\epsilon=0.15$. The centers of red and yellow circles in P
indicate two distinct points which have a large difference in $e$ but with
similar $\phi$, $e(\mathbf{x}_{\mathrm{red}})=142$ and
$e(\mathbf{x}_{\mathrm{yellow}})=29$ while
$\phi(\mathbf{x}_{\mathrm{red}})\approx\phi(\mathbf{x}_{\mathrm{yellow}})\approx
0.08$, as shown in E and F. (N and O) Zoomed-in views of the extremes of $e$
values, $e=29$ (N) and $e=142$ (O), for rods of aspect ratio $\alpha=100$ with
similar volume fractions $\phi=0.08$, from locations denoted by yellow and red
circles in the second and fourth rows of M.

Figure 3: Vertical variations in mesoscopic fields. (A) Normalized volume
fraction in the cylindrical slab at height $z$, $\tilde{\phi}(z)=\phi(z)/\Phi$
(see Materials and Methods for details), where $\Phi$ is the global volume
fraction domain. Inset: Two planes at $z-h/2$ and $z+h/2$ in a rod packing.
(B) Orientational order of rods within a slice at location $z$. (C) Normalized
contact density $\tilde{c}(z)$ within a slice at location $z$. (D) Normalized
entanglement at $z$-layer,
$\tilde{e}(z)=\sum_{z^{\prime}}\tilde{e}(z,z^{\prime})$ where
$\tilde{e}(z,z^{\prime})$ is given by Eq. S5 (E and F) Distance matrix plots
of $\tilde{e}(z,z^{\prime})$ for $\alpha=38$ (E) and $\alpha=100$ (F). (G and
H) A color-coded visualization of two different rod packings with aspect
ratios $\alpha=38$ (G) and $\alpha=100$ (H). (I) A tentative entanglement
phase diagram with respect to the reciprocal of the coefficient of variation
($\mu/\sigma$), dynamic load ($a/g$), and static load
($F/\rho_{\mathrm{s}}gd^{2}l$). Here $\rho_{\mathrm{s}}$ is the solid density.
In our experimental data, the critical $\mu/\sigma$ value was around unity,
for a packing of $\alpha\approx 100$ and a dynamic load $a/g\approx 2.5$.


# Cavity mediated level attraction and repulsion between magnons

Jayakrishnan M. P. Nair<EMAIL_ADDRESS>Institute for Quantum
Science and Engineering, Texas A$\&$M University, College Station, TX 77843,
USA  Department of Physics and Astronomy, Texas A$\&$M University, College
Station, TX 77843, USA  Debsuvra Mukhopadhyay<EMAIL_ADDRESS>Institute
for Quantum Science and Engineering, Texas A$\&$M University, College Station,
TX 77843, USA  Department of Physics and Astronomy, Texas A$\&$M University,
College Station, TX 77843, USA  Girish S. Agarwal<EMAIL_ADDRESS>Institute for Quantum Science and Engineering, Texas A$\&$M University,
College Station, TX 77843, USA  Department of Physics and Astronomy, Texas
A$\&$M University, College Station, TX 77843, USA  Department of Biological
and Agricultural Engineering, Texas A$\&$M University, College Station, TX
77843, USA

###### Abstract

We characterize some of the distinctive hallmarks of magnon-magnon interaction
mediated by the intracavity field of a microwave cavity, along with their
testable ramifications. In general, we foreground two widely dissimilar
parameter domains that bring forth the contrasting possibilities of level
splitting and level crossing. The former is observed in the regime of strong
magnon-photon couplings, particularly when the three modes bear comparable
relaxation rates. This character is marked by the appearance of three
distinguishable and non-converging polariton branches in the spectral response
to a cavity drive. However, when the bare modes are resonant and the couplings
perfectly symmetrical, one of the spectral peaks gets wiped out. This
anomalous extinction of polaritonic response can be traced down to the
existence of a conspicuous dark mode alongside two frequency-shifted bright
modes. In an alternate parameter regime, where the magnon modes are weakly
coupled to the cavity, features of level attraction unfold, subject to a large
relaxation rate for the cavity mode. Concurrently, for antisymmetric detunings
to the magnon modes, a transmission window springs into existence, exhibiting
transparency in the limit of negligible dissipation from the magnons. The
emergence of level attraction can be reconciled with a theoretical model that
embodies the dynamics of the magnon-magnon subsystem when the cavity field
decays rapidly into its steady state. In this limit, we identify a purely
dissipative coupling between the magnon modes.

##  I Introduction

The expedient control of coherent magnon-photon coupling in hybrid cavity-
spintronic systems has laid bare a plethora of new avenues for magnon-based
quantum information transfer Huebl  _et al._ (2013); Tabuchi  _et al._ (2014);
Zhang  _et al._ (2014); Bai  _et al._ (2015); Zhang  _et al._ (2015); Tabuchi
_et al._ (2015); Chumak  _et al._ (2015); Zhang  _et al._ (2016); Osada  _et
al._ (2016); Bai  _et al._ (2017); Zhang  _et al._ (2017); Li  _et al._ (2019,
2018); Zhang  _et al._ (2019); Nair and Agarwal (2020); Wang  _et al._ (2018);
Lachance-Quirion  _et al._ (2019); Nair  _et al._ (2020, 2021a). The
collective spin-wave excitations called magnons in magnetic systems can
efficiently interface with microwave photons, thereby consolidating the
strength of dispersive interaction to produce well-separated hybridized states
demonstrating level repulsion and Rabi oscillations Khitrova  _et al._ (2006);
Agarwal (1984). Complementary to the phenomenon of coherent magnon-photon
coupling in a cavity, there exists a paradigm of dissipative coupling wherein
the magnon mode interacts with the modes of a waveguide Harder  _et al._
(2018); Yao  _et al._ (2019). In contrast to the mode anticrossing in
dispersively coupled systems, dissipative coupling induces level attraction
characterized by the coalescence of the polariton modes at the exceptional
points (EPs) Bhoi  _et al._ (2019); Yang  _et al._ (2019); Rao  _et al._
(2019); Xu  _et al._ (2019); Wang  _et al._ (2019). Considering the ubiquity
of both coherent and dissipative couplings in cavity-magnonics as well as the
wide disparity in the resulting physics, the competing effects of these
couplings have evoked significant interest of late. Consequently, cavity-
magnonics has emerged as a premier platform to translate some of the foremost
ideas of non-Hermitian optics into scalable quantum technologies with novel,
multitasking capabilities. By exploiting strongly coupled magnon-photon
systems, a plenitude of exotic effects have been brought to the fore, some of
which include single-spin detection Lachance-Quirion  _et al._ (2017) and
manipulation of spin current Bai  _et al._ (2015); Maier-Flaig  _et al._
(2016), polaritonic dark modes Zhang  _et al._ (2015); Lambert  _et al._
(2016) and magnonic quintuplets Yao  _et al._ (2017), observation of coherent
perfect absorption and exceptional point Zhang  _et al._ (2017),
multistability Wang  _et al._ (2018); Nair  _et al._ (2020), squeezing Li  _et
al._ (2019) and entanglement Li  _et al._ (2018); Zhang  _et al._ (2019); Nair
and Agarwal (2020), sensing of anharmonicities Nair  _et al._ (2021b), and
bidirectional microwave-to-optical transduction Hisatomi  _et al._ (2016); Ihn
_et al._ (2020); Mukhopadhyay  _et al._ (2021).

Owing to its potential applications in topological energy transfer Dembowski
_et al._ (2001); Doppler  _et al._ (2016); Xu  _et al._ (2016) and quantum
sensing Wiersig (2014); Chen  _et al._ (2017); El-Ganainy  _et al._ (2018);
Özdemir  _et al._ (2019); Miri and Alu (2019); Zhang  _et al._ (2020), there
has been a proliferation of efforts within the optics and condensed matter
community to realize level attraction in hybrid non-Hermitian systems. So far,
this effect has been successfully showcased in cavity-optomechanical systems
Bernier  _et al._ (2018), and cavity-magnonic configurations integrated with
waveguides Yang  _et al._ (2020) and split-ring resonators Bhoi  _et al._
(2019); Yang  _et al._ (2019). It is well known that two coherently coupled
harmonic systems naturally exhibit level repulsion while dissipatively coupled
systems are characterized by level attraction. Therefore, a customary
technique to implement level crossing lies in mitigating the effect of
dispersive interaction relative to its dissipative counterpart. Both
dispersive and dissipative couplings naturally occur in qubits interacting
with electromagnetic vacuum or with the modes of a bad cavity Agarwal (1974).

In the context of magnonics, magnon-photon level repulsion and level
attraction have both been experimentally realized Harder  _et al._ (2018); Yao
_et al._ (2019); Bhoi  _et al._ (2019); Yang  _et al._ (2019); Rao  _et al._
(2019); Xu  _et al._ (2019); Wang  _et al._ (2019). In this work, we explore
the distinctive qualities of magnon-magnon hybridization in a system of two
non-interacting magnetic systems coherently coupled to a microwave cavity. We
show that the mediating effect of the intracavity field brings to bear the
dual possibilities of effecting both level anticrossing and level crossing in
two distinct parameter regimes. Much to our anticipation, strong magnon-photon
couplings result in well-separated polaritonic states, or normal modes,
characterizing level repulsion. Curiously, however, we observe that the
central peak is washed out when the frequencies of the bare modes are matched
to resonance and the magnons are symmetrically coupled to the cavity mode. As
the symmetry is broken and the magnons are antisymmetrically detuned relative
to the cavity mode, the dark polariton gets reexcited on account of energy
transfer from the bright modes. The contrary feature of level attraction
transpires when the cavity leakage rate strongly dominates the dynamics, i.e.,
the magnon-photon couplings as well as the magnonic damping rates are much
weaker than the cavity relaxation rate. In the same parameter regime, on
examining the spectral properties of the hybrid system, we unveil a narrow
transmission window, with the possibility of transparent behavior in the limit
of negligible magnon dissipation. Both these features, i.e., the inception of
level crossing and the appearance of transparency, are borne out by a purely
dissipative form of interaction that ensues between the magnons when the
cavity field is adiabatically eliminated. In fact, the adiabatic limit
reproduces the characteristic dynamics of two spatially separated systems
coupled through an interceding reservoir, thereby unveiling an equivalence
with waveguide QED. It should be noted that qubits get coherently coupled with
each other in a dispersive cavity Agarwal  _et al._ (1997). In a cavity-
magnonic setup, such a coupling has been experimentally realized in the
dispersive limit, demonstrating level repulsion between the magnons Lambert
_et al._ (2016). In contrast, the level repulsion discussed in our paper makes
no reference to the dispersive regime, and the magnons could well be on
resonance with cavity mode.

This paper is organized as follows. In section II, we spell out the
theoretical model describing a pair of magnetic systems coupled to a microwave
cavity, and highlight the two disparate regimes in which level splitting and
level crossing could be observed. Following this, we study the induced
magnetization in response to a cavity drive in section III, to reaffirm the
contrasting spectral signatures to both these behaviors. In section IV, we
address the emergence of a dark polaritonic state and furnish the physical
mechanism governing the excitation of this mode. Finally, in section V, we
illustrate the generation of a transmission window as a consequence of
dissipative coupling between the two magnons.

Figure 1: Schematic of two YIGs coherently coupled to a single-mode microwave
cavity. The static magnetic field exciting the Kittel mode in both the YIGs
are aligned along the $z$-axis. The intracavity field mode is chosen to be
polarized along the $y$-axis, with the corresponding magnetic field directed
along the $x$-axis.

##  II Theoretical Model

Our model comprises a hybrid cavity-magnonic system, in which two macroscopic
yttrium iron garnet (YIG) samples are placed inside a microwave cavity, as
illustrated in Fig. 1. The YIG is a high-spin-density magnetic material, and
the low lying excitations of the collective angular momentum of these spins in
such magnetically ordered systems gives rise to quasi-particles, namely
magnons. In our model, we consider a spatially uniform Kittel mode of either
YIG as being dispersively coupled to the intracavity photons of a neighboring
frequency. Therefore, in essence, the magnetization associated with the YIG
sphere can simply be written as the magnetization of a ferromagnet with a
single large spin S, i.e., $\textbf{M}=\gamma_{e}\textbf{S}/V$, where
$\gamma_{e}=e/{m_{e}c}$ is the gyromagnetic ratio for electron spin and S
denotes the collective spin operator, and $V$ the volume of the YIG sphere.
The Hamiltonian of the hybrid magnon-cavity system is then provided by

$\begin{split}H/\hbar=-\gamma_{e}\sum_{i=1}^{2}&B_{0}S_{i,z}+\omega_{a}a^{\dagger}a-\gamma_{e}\sum_{i=1}^{2}S_{i,x}B_{i,x},\end{split}$
(1)

where $\textbf{B}_{0}$ is the applied bias magnetic field along the $z$
direction, $S_{i,x}$ $(S_{i,z})$ denotes the $x$\- ($z$-) component of the
collective spin S of the $i^{th}$ YIG, $\omega_{a}$ is the cavity resonance
frequency, and $a$ $(a^{\dagger})$ represents the annihilation (creation)
operator of the cavity field. The magnetic field of the cavity field
($B_{n,x}$) is assumed to be along the $x$ direction. A typical YIG sample
having diameter $d=1$ mm and an approximate spin density $\rho=4.22\times
10^{27}$ m-3 leads to the total number of spins $N\approx 10^{18}$. In the
limit of a high spin density, we invoke the Holstein-Primakoff transformation
to recast the raising and lowering operators as
$S_{i}^{+}=\sqrt{2S_{i}}m_{i},\ S_{i}^{-}=\sqrt{2S_{i}}m_{i}^{\dagger}$, with
$i=1,2$ labeling the two YIGs) and $m_{i}$ ($m_{i}^{\dagger}$) operators
satisfying the Bosonic algebra. With such simplifications, we can now tailor
the Hamiltonian in Eq. (1) into its corresponding operator representation as

$H/\hbar=\omega_{a}a^{\dagger}a+\sum_{i=1}^{2}\Big{[}\omega_{i}m_{i}^{\dagger}m_{i}+g_{i}(m_{i}^{\dagger}a+m_{i}a^{\dagger})\Big{]}$
(2)

where the frequency of Kittel mode $\omega_{i}=\gamma_{e}B_{0}$ with
$\gamma_{e}/2\pi=28$ $\text{GHz}/\text{T}$. The parameter
$g_{i}=\frac{\sqrt{5}}{2}\gamma_{e}\sqrt{N}B_{\text{vac}}$ quantifies the
coherent magnon-cavity coupling, with
$B_{\text{vac}}=\sqrt{2\pi\hbar\omega_{c}/\text{V}}$ denoting the magnetic
field of vacuum.

Figure 2: Level anticrossing and level crossing in two different regimes of
cavity QED; (a) Eigenfrequencies of the cavity-magnon system for $g=2$ and
$\gamma_{i}=1$ exhibiting level splitting; (b) Eigenfrequencies and (c)
Linewidths in the case of $g=0.2$ and $\gamma_{i}=0.01$, with the features of
level crossing. All the parameters are chosen in units of $\kappa$.

In the frame rotating at frequency $\omega_{a}$, the effective non-Hermitian
Hamiltonian describing the time evolution of the column matrix $X$, where
$X^{T}=$ [$\expectationvalue{a}$ $\expectationvalue{m_{1}}$
$\expectationvalue{m_{2}}$] is provided by

$\displaystyle\mathcal{H}=\begin{pmatrix}-i\kappa&g_{1}&g_{2}\\\
g_{1}&s-i\gamma_{1}&0\\\ g_{2}&0&-s-i\gamma_{2}\end{pmatrix},$ (3)

where $s=(\omega_{m_{1}}-\omega_{m_{2}})/2$, and we have set
$\omega_{a}=(\omega_{m_{1}}+\omega_{m_{2}})/2$. The parameters $\kappa$ and
$\gamma_{i}$ denote the rates of dissipation appearing in the master equation.
In what follows, we suppress the notation $\expectationvalue{.}$ in the mode
amplitudes, as we are interested in a semiclassical description of the system.
The three complex eigenvalues of $\mathcal{H}$, resulting from the
characteristic polynomial equation

$\displaystyle(\lambda+i\kappa)(\lambda+s+i\gamma_{1})(\lambda-s+i\gamma_{2})-g_{1}^{2}(\lambda+s+i\gamma_{2})$
$\displaystyle-g_{2}^{2}(\lambda-s+i\gamma_{1})=0,$ (4)

identify the normal modes of the hybridized system, also referred to as the
polaritons. Here, we present two distinct parametric regimes: (i) where
$g\gg\kappa,\gamma_{i}$ with the features of level repulsion and (ii), the
weak coupling regime $g,\gamma_{i}\ll\kappa$ where the eigenvalues relevant to
the system demonstrates level attraction. Before delving into the two cases,
we recall that a coherently coupled system typically demonstrates Rabi
splitting in the strong-coupling regime. That is to say, coherent coupling
opens up a gap between the normal mode frequencies of the system, and the
eigenfrequencies bend away from each other. In contrast, dissipatively coupled
systems exhibit level-attraction, which is characterized by smoothly
converging eigenfrequencies leading to mode coalescence over a finite
parameter domain. The points of convergence, known as EPs, are responsible for
a variety of exotic phenomena owing to the square-root singularity of the
modes at these points. We demonstrate in Fig. 2, the solutions of the
characteristic equation (4) in the two above mentioned parameter regimes where
level repulsion and level attraction come to light. Specifically, in case
(ii), when the cavity is weakly coupled to the magnon modes and the photons
decay over a much shorter time scale than the magnons, we can observe
signatures of level crossing between the magnon-like polaritons. The
contrasting features in Fig. 2(a, b) can be analytically understood by
inspecting the normal modes in these special regimes, as spelled out below.

For $\kappa=\gamma_{1}=\gamma_{2}$ and $g_{1}=g_{2}=g$, the three eigenvalues
of $\mathcal{H}$ are given by $\lambda_{0}=-i\kappa$,
$\lambda_{\pm}=-i\kappa\pm\sqrt{s^{2}+2g^{2}}$. The imaginary parts of the
eigenvalues, which quantify the corresponding linewidths, are all identical
and insensitive to the variations in $s$. The real parts of these eigenvalues
are plotted against $s$ in Fig. 2(a), where the minimum frequency gap between
$\lambda_{+}$ and $\lambda_{-}$, equaling $2\sqrt{2}g$, is a manifestation of
level repulsion in the system. Note that the minimum separation must be larger
than the linewidth $\kappa$ to observe level repulsion. However, in the limit
$\kappa\gg\gamma_{i}$ and $\kappa\gg g_{i}$, level crossing is introduced over
a small band of frequency detunings $s$, as portrayed in Fig. 2(b, c). We have
assumed that $\gamma_{1}=\gamma_{2}=\gamma$ and $g_{1}=g_{2}=g$. It is,
therefore, imperative to design a cavity with low quality factor to observe
this effect. Neglecting the effect of $\gamma$, which is weak in reference to
$\kappa$, leads to the characteristic equation

$\displaystyle(\lambda+i\kappa)(\lambda^{2}-s^{2})-2g^{2}\lambda=0.$ (5)

Since the eigenvalues in the absence of $g$ equal $\lambda_{0}=-i\kappa$, and
$\lambda_{\pm}=\pm s$, it can be argued that the perturbative correction to
each of them for sufficiently small $g$ would go as $\mathcal{O}(g^{2})$. This
consideration pins down the approximate form of $\lambda_{0}$ to be
$-i\kappa[1-2g^{2}/(s^{2}+\kappa^{2})]$. Next, in view of the fact that both
$s$ and $g^{2}$ are small compared to $\kappa$, we approximate
$\lambda_{\pm}+i\kappa\approx i\kappa$, which permits the reduction of the
cubic equation into a quadratic one
$\lambda_{\pm}^{2}+2\Gamma\lambda_{\pm}i-s^{2}\approx 0$, where
$\Gamma=g^{2}/\kappa$. This yields the remaining eigenvalues

$\displaystyle\lambda_{\pm}=-i\Gamma\pm i\sqrt{\Gamma^{2}-s^{2}}.$ (6)

The forms of $\lambda_{\pm}$ make apparent the level crossing between them,
with EPs located at $s=\Gamma$. The plots in Fig. 2(b, c) are found to be in
remarkable agreement with the expressions obtained above.

In the following section, we extract spectroscopic information about the
system by applying an external probe, which renders an empirical tool to
observe level crossing and level anticrossing.

Figure 3: (Panel A; a-e) Sumtotal of the spincurrent response
$\absolutevalue{m_{1}}^{2}+\absolutevalue{m_{2}}^{2}$ plotted as functions of
$\Delta_{a}$ in the case of strong magnon-photon coupling $g=2$ and equal
relaxation rates $\gamma_{i}=1$ (similar to Fig. 2(a)). Note the conspicuous
extinction of the central resonance in (A; c). (Panel B; a-e) Profiles for the
case of $g=0.2$, $\gamma_{i}=0.01$ (similar to Fig. 2(b,c)). Only two peaks
appear in this scenario. Drive power $P=1$ mW.

##  III Level repulsion and attraction between magnons

In the presence of an external electromagnetic field at frequency $\omega_{d}$
driving the cavity, the Hamiltonian in Eq. (2) is superseded by

$\begin{split}H/\hbar=\omega_{a}a^{\dagger}a+\sum_{i=1}^{2}\Big{[}\omega_{i}m_{i}^{\dagger}m_{i}+g_{i}(m_{i}^{\dagger}a+m_{i}a^{\dagger})\Big{]}\\\
+i\sqrt{\kappa}\mathcal{E}[a^{\dagger}e^{-i\omega_{d}t}-ae^{i\omega_{d}t}],\end{split}$
(7)

with $\mathcal{E}=\sqrt{\frac{P}{\hbar\omega_{d}}}$ being the strength of
external drive and $P$ the driving power. By switching to the rotating frame
of the drive, we deduce the mean-field equations for the cavity-magnon
dynamics to be

$\displaystyle\dot{X}=-i(\mathcal{H}-\Delta)X+F,$ (8)

where we have separated out the cavity detuning
$\Delta=\omega_{d}-\omega_{a}$, and $F^{T}=$$\sqrt{\kappa}$[$\mathcal{E}\hskip
2.84526pt0\hskip 2.84526pt0$]. The solution to the steady state can,
therefore, be expressed as

$\displaystyle X(\Delta)=-i(\mathcal{H}-\Delta)^{-1}F,$ (9)

With the poles of this expression being accorded by the eigenvalues of
$\mathcal{H}$, key information regarding the polaritonic properties would be
encoded in the empirically observable response, such as the steady-state
spincurrents $\absolutevalue{m_{1}}^{2}$ and $\absolutevalue{m_{2}}^{2}$. In
this section, we provide a numerical solution to the spectral features by
scanning the probe frequency.

In Fig. 3, on either of the panels (A) and (B), we plot the sum of the
spincurrents from the individual magnon modes, i.e.,
$\absolutevalue{m_{1}}^{2}+\absolutevalue{m_{2}}^{2}$ as a function of
$\Delta$ for five different values around $s=0$. The left panel (A) refers to
the barely strong-coupling regime, i.e., $g_{i}>\kappa,\gamma_{i}$, and
$\kappa=\gamma_{i}$. The total spincurrent, for $s\neq 0$, exhibits three
distinct peaks corresponding to the resonant frequencies of the system. The
central peak at $\Delta=0$ is flanked on either side by two symmetrical
polaritonic peaks. The disappearance of the central peak at $s=0$ is symbolic
of a dark polariton Yuen-Zhou and Menon (2019). The underlying physics of this
dark mode would be discussed at length in the subsequent section. From the
figures (a) to (e) in the panel (A), it is evident that the two magnon-like
polaritons move closer as $s\rightarrow 0$, leaving a minimum frequency gap
between the two levels at $s=0$. Since the levels never intersect, it is
referred to as mode anticrossing. However, this picture incurs a dramatic
inversion as we transition into the domain of weaker couplings and higher
relaxation rates of the cavity photons. Figure 3, panel (B), displays the
sumtotal of the spincurrents as a function of $\Delta$ in the limit
$\kappa\gg\gamma_{i}$ and $\kappa\gg g_{i}$, for five different values of $s$.
In this scenario, around $s=0$, the magnon-like polaritons coalesce together
into a single peak. Note that the cavity-like central resonance remains
strongly quenched on account of its large linewidth. The onset of level
crossing points to a hidden dissipative interaction between the magnon modes.
In Appendix A, we establish parity with these numerical results by uncovering
this dissipative magnon-magnon coupling through the adiabatic elimination of
the cavity mode.

##  IV Excitation of dark polariton via symmetry-breaking

As evidenced by the numerical illustration in Fig. 3, panel (A), when the
magnon-photon coupling is adequately strong, the tripartite cavity-magnon
system is endowed with three distinct spectral peaks characterizing the
polaritonic modes of the hybrid composite. The two peripheral peaks are pushed
inward as $s$ becomes smaller, settling for a minimum separation at resonance.
This limiting behavior is marked simultaneously by the complete extinction of
the central peak, laying bare a characteristic dark mode when the conditions
conform to perfect symmetry, i.e., when all three modes have identical
detunings, and the magnons are symmetrically coupled to the intracavity field.
We can delineate a physical argument supporting this observation by simply
inspecting the normal modes of the undamped system. In the rotating frame of
the cavity field, the Hamiltonian is given by

$\displaystyle H=g[a^{\dagger}(m_{1}+m_{2})+h.c].$ (10)

For this Hermitian system, the polaritonic modes turn out to be
$A=\dfrac{m_{1}-m_{2}}{\sqrt{2}}$,
$B=\dfrac{1}{\sqrt{2}}\bigg{[}a+\dfrac{m_{1}+m_{2}}{\sqrt{2}}\bigg{]}$, and
$C=\dfrac{1}{\sqrt{2}}\bigg{[}a-\dfrac{m_{1}+m_{2}}{\sqrt{2}}\bigg{]}$,
corresponding respectively to the eigenfrequencies $\lambda_{0}=0$,
$\lambda_{+}=\sqrt{2}g$, and $\lambda_{-}=-\sqrt{2}g$, which recasts the
Hamiltonian into the form

$\displaystyle H=\sqrt{2}g[B^{\dagger}B-C^{\dagger}C].$ (11)

Thus, $B$ and $C$ appear as two counteroscillating normal modes of the
symmetric system, with the frequency-splitting determined by the coupling
strength $g$. Note, however, the conspicuous absence of the mode $A$ in Eq.
(11). This mode, therefore, represents a dark polariton, whereas $B$ and $C$
identify the bright polaritons with detectable spectral signatures. We note
that dark polaritons have been extensively discussed in the context of
molecules interacting with cavities Yuen-Zhou and Menon (2019). This is what
is manifested in Fig. 3 (A; c). The picture obviously changes if we begin to
steer away from the symmetry, which are exemplified by the other subfigures in
the same panel. For non-zero value of $s$, the Hamiltonian transforms into
$s(m_{1}^{\dagger}m_{1}-m_{2}^{\dagger}m_{2})+g[a^{\dagger}(m_{1}+m_{2})+h.c]$.
This can be expressed in the original polaritonic representation as

$\displaystyle
H=\sqrt{2}g[B^{\dagger}B-C^{\dagger}C]+\frac{s}{\sqrt{2}}(A^{\dagger}B+AB^{\dagger})-\frac{s}{\sqrt{2}}(A^{\dagger}C+AC^{\dagger}).$
(12)

The additional terms in the new Hamiltonian, originating from the frequency
discord between the magnons, are clearly indicative of energy transitions from
$B$ and $C$ to $A$ and vice versa. These new channels for energy exchange
foreshadow how a mismatch in the magnon detunings would be instrumental to the
illumination of the non-radiative polariton. The stimulating effect of
symmetry-breaking can be experimentally observed via the spectroscopic
properties of the system when it is irradiated by a probe field.

To analytically clarify the polaritonic suppression at resonance, let us
examine the spectral signature of a real dissipative system. To that end, we
evaluate the expressions for the steady-state response functions

$\displaystyle
m_{j}=-i\sqrt{\kappa}\mathcal{E}[(\mathcal{H}-\Delta)^{-1}]_{j1},$ (13)

where $\mathcal{H}$ is the effective Hamiltonian for the $3\crossproduct 3$
system defined in Eq. (3), and $j=1\text{ or }2$ labels the magnon.
Substituting the expression for $\mathcal{H}$ from Eq. (3), we obtain,

$\displaystyle
m_{1}=\frac{-i\sqrt{\kappa}\mathcal{E}g_{1}(s+\Delta+i\gamma_{1})}{D}$ (14)
$\displaystyle
m_{2}=\frac{i\sqrt{\kappa}\mathcal{E}g_{2}(s-\Delta-i\gamma_{1})}{D},$ (15)

where $D=\det(\mathcal{H}-\Delta)$ given by

$\displaystyle\begin{split}D=(\Delta+i\kappa)(s-\Delta-i\gamma_{1})(s+\Delta+i\gamma_{2})+g_{1}^{2}(s+\Delta+i\gamma_{2})\\\
-g_{2}^{2}(s-\Delta-i\gamma_{1}).\end{split}$

In the event of symmetrical detunings ($s=0$) and equal damping rates
($\kappa=\gamma_{1}=\gamma_{2}$), the eigenvalues go as
$\lambda_{0}=-i\kappa$, $\lambda_{\pm}=\pm\sqrt{2}g-i\kappa$. The common
expression for the spincurrents then reduces to

$\displaystyle
m_{j}=\frac{i\sqrt{\kappa}g\mathcal{E}}{(\Delta+\sqrt{2}g+i\kappa)(\Delta-\sqrt{2}g+i\kappa)},$
(16)

for $j=1,2$. Clearly, the contribution from the pole at $\lambda_{0}$ remains
suppressed when the modes are on resonance and the magnon-photon couplings are
perfectly symmetrical. Given the equality of the two spincurrents, it is no
coincidence that the normal mode $A=\dfrac{m_{1}-m_{2}}{\sqrt{2}}$ remains
unpopulated in the steady state. On the other hand, if a small asymmetry is
introduced in the frequencies by letting $s$ become non-zero, the eigenvalues
get modified into $\lambda_{0}=-i\kappa$ and
$\lambda_{\pm}=\pm\sqrt{s^{2}+2g^{2}}-i\kappa$. In this case, the pole at
$\lambda_{0}$ leaves a detectable signature as underscored by the
corresponding spincurrents

$\displaystyle
m_{1}=\frac{i\sqrt{\kappa}g\mathcal{E}(\Delta+s+i\kappa)}{(\Delta+i\kappa)(\Delta+\sqrt{s^{2}+2g^{2}}+i\kappa)(\Delta-\sqrt{s^{2}+2g^{2}}+i\kappa)},$
$\displaystyle
m_{2}=\frac{i\sqrt{\kappa}g\mathcal{E}(\Delta-s+i\kappa)}{(\Delta+i\kappa)(\Delta+\sqrt{s^{2}+2g^{2}}+i\kappa)(\Delta-\sqrt{s^{2}+2g^{2}}+i\kappa)},$
(17)

which also remain unequal insofar as $s\neq 0$. This explains the emergence of
the central peak in Fig. 3 as a signature of population transfer to the mode
$A$. The corresponding peak height is provided by the expression

$\displaystyle\absolutevalue{m_{1}}^{2}+\absolutevalue{m_{2}}^{2}=\frac{2g^{2}\mathcal{E}^{2}(s^{2}+\kappa^{2})}{\kappa(s^{2}+2g^{2}+\kappa^{2})}.$
(18)

Clearly, the height reduces to a minimum at $s=0$.

##  V Dissipative-Coupling-Induced Transparency

Having explored some of the distinctive regimes of cavity QED, which uncovered
the contrasting features of level splitting and level crossing, it is
worthwhile to study the transmission characteristics of the cavity-magnon
system. For this purpose, we restrict our analysis to the regime where
$\kappa\gg\gamma_{i}$ and $\kappa\gg g_{i}$, and shed light on the possibility
of effecting transparent behavior. Note that this is precisely the parameter
domain that elicits level attraction.

Figure 4: Refection spectra for (a) $s=0.01$ with $\gamma_{i}=0$, and (b)
$s=0.04$ with $\gamma_{i}=0.01$. Other parameters are $g=0.2$ and $\kappa=1$,
which yields $\Gamma=g^{2}/\kappa=0.04$.

For the cavity-driven system described by Eq. (7), we can obtain the the
reflected and transmitted output fields from the input-output formalism, with
the relevant input-output relations being

$\displaystyle b_{out}$ $\displaystyle=\sqrt{\kappa}a,$ $\displaystyle
a_{out}+\mathcal{E}$ $\displaystyle=\sqrt{\kappa}a,$ (19)

where $a=-i\sqrt{\kappa}[(\mathcal{H}-\Delta)^{-1}]_{11}\mathcal{E}$ is the
steady-state intracavity field obtained by setting $\dot{X}=0$ in Eq. (8).
Substituting this into the aforementioned relations, we extract the
transmission and reflection coefficents to be

$\displaystyle t$
$\displaystyle=\frac{i\kappa(\Delta-s+i\gamma)(\Delta+s+i\gamma)}{(\Delta+i\kappa)(\Delta-s+i\gamma)(\Delta+s+i\gamma)-2g^{2}(\Delta+i\gamma)},$
$\displaystyle r$
$\displaystyle=-1+\frac{i\kappa(\Delta-s+i\gamma)(\Delta+s+i\gamma)}{(\Delta+i\kappa)(\Delta-s+i\gamma)(\Delta+s+i\gamma)-2g^{2}(\Delta+i\gamma)}.$
(20)

In Figs. 4(a, b), we sketch the reflection profiles for a couple of values of
$s$, vindicating the existence of a transmission window. In the limit
$\gamma/\kappa\rightarrow 0$, perfect transparency is observed at $\Delta=0$,
as demonstrated in Fig. 4(a). Similarly, Fig. 4(b) depicts the reflection
spectra at the exceptional point $s=\Gamma=0.04$ for $\gamma_{i}=0.01$,
exhibiting a transmission window. The transparency owes its origin to a
dissipative form of magnon-magnon coupling surfacing in this limit of
negligible magnon dissipation. The details on this origin are elaborated in
Appendix A. At $\Delta=0$, the expressions in Eq. (20) reduce to

$\displaystyle t$
$\displaystyle=\frac{s^{2}+\gamma^{2}}{s^{2}+\gamma^{2}+2\gamma\Gamma},$
$\displaystyle r$
$\displaystyle=-\frac{2\gamma\Gamma}{s^{2}+\gamma^{2}+2\gamma\Gamma},$ (21)

where $\Gamma=g^{2}/\kappa$. As long as $\gamma\Gamma\ll s^{2}$, the
approximate behavior would be given by $t\approx
1-\mathcal{O}(2\gamma\Gamma/s^{2})$ and
$r\approx\mathcal{O}(2\gamma\Gamma/s^{2})$.

As $s$ becomes smaller, in the limit $\gamma/\kappa\approx 0$, the EIT-like
window gets progressively narrower leading to a sharpening of the central dip.
This can be understood from the behavior of the eigenvalues derived earlier in
Sec. II. Through pertinent approximations, we had obtained them to be
$\lambda_{0}=-i\kappa[1-2g^{2}/(s^{2}+\kappa^{2})]$ and
$\lambda_{\pm}=-i\Gamma\pm i\sqrt{\Gamma^{2}-s^{2}}$. While $\lambda_{0}$ had
rather pedestrian features, the other two normal modes
$\lambda_{\pm}\approx\lambda_{\pm}=-i\Gamma\pm i\sqrt{\Gamma^{2}-s^{2}}$ were
the source of level attraction. In the same vein, the unfolding of the sharp,
precipitous dip in the reflection lineshapes can be traced down to the
linewidths of the modes $\lambda_{\pm}$. Apropos of our parameter choices,
when $s/\Gamma$ is tiny compared to $1$, as is the case in Fig. 4(a), they can
be further approximated as

$\displaystyle\lambda_{+}$ $\displaystyle\approx-\frac{is^{2}}{2\Gamma},$
$\displaystyle\lambda_{-}$ $\displaystyle\approx-2i\Gamma.$ (22)

While both of them have zero real parts, $\imaginary(\lambda_{+})$ pales in
comparison to $\imaginary(\lambda_{-})$. In fact, the pole at $\lambda_{+}$
has a vanishingly small linewidth, which explains the sharp and narrow trough
in the reflection spectrum. The effect of $\lambda_{-}$, on the other hand, is
embodied in the Lorentzian envelope circumscribing both the peaks flanking the
central minimum.

##  VI Conclusions

In conclusion, we have discussed the feasibility of engineering both level
repulsion and level attraction between two magnons coupled to a single-mode
cavity, by manipulating the relaxation rates and the magnon-photon coupling
strengths. In the weak-coupling regime, when the cavity mode relaxes
considerably faster than the magnon modes, an effective dissipative
interaction is engendered between the two magnons, which is directly
responsible for the unfolding of level attraction as well as a narrow
transmission window. The level crossing is manifested as the coalescence of
polaritonic peaks in the spincurrent response of the hybridized system. On the
contrary, the strong-coupling regime, for comparable relaxation rates, results
in distinguishable polaritonic peaks with the earmarks of level repulsion.
Interestingly, when the magnons are tuned in resonance to the cavity together
with symmetrical couplings, the spectral response is characterized by the
appearance of a dark polaritonic state. The central dip designating the dark
polariton transitions into a peak as we deviate from resonance, indicating the
transfer of population into the dark mode. It could be noted that the bulk of
the precursive literature on level attraction in cavity-magnonics was grounded
in the study of magnon-photon coupling and the resulting polariton
characteristics. On the other hand, the level attraction highlighted in our
paper, is rooted in an effective magnon-magnon coupling infused via the
photonic interaction with the individual magnon modes.

##  VII Acknowledgements

The authors acknowledge the support of The Air Force Office of Scientific
Research [AFOSR award no FA9550-20-1-0366], The Robert A. Welch Foundation
[grant no A-1943] and the Herman F. Heep and Minnie Belle Heep Texas A&M
University endowed fund.

##  Appendix A Adiabatic Theory of Level Attraction

To explain the physics of our observation, in the limit $\kappa\gg\gamma_{i}$
and $\kappa\gg g_{i}$, we employ an adiabatic model to advance an effective
two-mode description for the magnon modes, since the cavity mode relaxes
rapidly into its steady state. This treatment reinforces forthwith the
emergence of a dissipative interaction between the magnons in the adiabatic
limit, which underpins the phenomenon of level attraction. In the rotating
frame of the cavity mode, the mean field equations of the hybrid cavity magnon
system are given by

$\displaystyle\dot{X}=-i\mathcal{H}X,$ (23)

In the adiabatic limit, we can eliminate the cavity mode by setting
$\dot{a}=0$, yielding

$\displaystyle a=\frac{-i(g_{1}m_{1}+g_{2}m_{2})}{\kappa}.$ (24)

Substituting this into Eq. (A1), we obtain a coupled system of differential
equations encompassing the dynamics of the magnon modes,

$\displaystyle\dot{Y}=-i\tilde{\mathcal{H}}Y,$ (25)

with ${Y}^{T}=$[$m_{1}$ $m_{2}$] and

$\displaystyle\tilde{\mathcal{H}}=\begin{pmatrix}{s}-i\tilde{\gamma_{1}}&\frac{-ig_{1}g_{2}}{\kappa}\\\
\frac{-ig_{1}g_{2}}{\kappa}&{-s}-i\tilde{\gamma_{2}}\end{pmatrix},$ (26)

where $\tilde{\gamma_{i}}=\gamma_{i}+\frac{g_{i}^{2}}{\kappa}$. For identical
couplings $g_{1}=g_{2}=g$ and magnon relaxation rates
$\gamma_{1}=\gamma_{2}=\gamma$, we can simplify $\tilde{\mathcal{H}}$ further
into

$\displaystyle\tilde{\mathcal{H}}=\begin{pmatrix}s-i(\gamma+\frac{g^{2}}{\kappa})&\frac{-ig^{2}}{\kappa}\\\
\frac{-ig^{2}}{\kappa}&-s-i(\gamma+\frac{g^{2}}{\kappa})\end{pmatrix}.$ (27)

The eigenvalues of the Hamiltonian are given by

$\displaystyle\lambda_{\pm}=-i\bigg{(}\gamma+\frac{g^{2}}{\kappa}\bigg{)}\pm\sqrt{s^{2}-(g^{2}/\kappa)^{2}},$
(28)

which bear the earmarks of level attraction, and also reproduce Eq. (6) in the
limit of vanishingly small $\gamma$. For,
$\absolutevalue{s}\leq\frac{g^{2}}{\kappa}$, the real parts remain identical,
yielding the phase transition points $s=\pm\frac{g^{2}}{\kappa}$. So it is no
coincidence that we had observed the onset of mode coalescence at this exact
same point in the preceding section. Thus, under the purview of the adiabatic
approximation, in which the short-lived nature of the intracavity field mode
is accounted for, the full tripartite system reduces to an effective two-mode
subsystem with conformable properties in the linear response. The above matrix
exactly replicates the dynamics of a two-mode system with density matrix
$\rho$, coupled through an interposing reservoir, and described by the master
equation Agarwal (2012)

$\frac{\differential\rho}{\differential
t}=-\frac{i}{\hbar}[{H},\rho]+\gamma_{1}\mathcal{L}(m_{1})\rho+\gamma_{2}\mathcal{L}(m_{2})\rho+2\Gamma\mathcal{L}(c)\rho,$
(29)

with $H=\Delta_{1}m_{1}^{\dagger}m_{1}+\Delta_{2}m_{2}^{\dagger}m_{2}$,
$\Gamma=g^{2}/\kappa$, $c=\frac{1}{\sqrt{2}}(m_{1}+m_{2})$, and the
Liouvillian $\mathcal{L}$ is defined by
$\mathcal{L}(\sigma)\rho=2\sigma\rho\sigma^{\dagger}-\sigma^{\dagger}\sigma\rho-\rho\sigma^{\dagger}\sigma$.
More precisely, the mean-field equations resulting from the above master
equation are identical to Eq. (A3) with an effective Hamiltonian given by
$\tilde{\mathcal{H}}$. To put it into perspective, the cavity, when weakly
coupled to the magnons, merely acts as a common reservoir for the magnon modes
in the adiabatic limit, thereby precipitating in a dissipative form of
interaction between them.

## References

  * Huebl  _et al._ (2013) H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, Phys. Rev. Lett. 111, 127003 (2013).
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 113, 083603 (2014).
  * Zhang  _et al._ (2014) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * Bai  _et al._ (2015) L. Bai, M. Harder, Y. P. Chen, X. Fan, J. Q. Xiao, and C.-M. Hu, Phys. Rev. Lett. 114, 227201 (2015).
  * Zhang  _et al._ (2015) X. Zhang, C.-L. Zou, N. Zhu, F. Marquardt, L. Jiang, and H. X. Tang, Nature communications 6, 1 (2015).
  * Tabuchi  _et al._ (2015) Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Science 349, 405 (2015).
  * Chumak  _et al._ (2015) A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Nature Physics 11, 453 (2015).
  * Zhang  _et al._ (2016) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Science advances 2, e1501286 (2016).
  * Osada  _et al._ (2016) A. Osada, R. Hisatomi, A. Noguchi, Y. Tabuchi, R. Yamazaki, K. Usami, M. Sadgrove, R. Yalla, M. Nomura, and Y. Nakamura, Phys. Rev. Lett. 116, 223601 (2016).
  * Bai  _et al._ (2017) L. Bai, M. Harder, P. Hyde, Z. Zhang, C.-M. Hu, Y. P. Chen, and J. Q. Xiao, Phys. Rev. Lett. 118, 217201 (2017).
  * Zhang  _et al._ (2017) D. Zhang, X.-Q. Luo, Y.-P. Wang, T.-F. Li, and J. You, Nature communications 8, 1 (2017).
  * Li  _et al._ (2019) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. A 99, 021801 (2019).
  * Li  _et al._ (2018) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. Lett. 121, 203601 (2018).
  * Zhang  _et al._ (2019) Z. Zhang, M. O. Scully, and G. S. Agarwal, Phys. Rev. Research 1, 023021 (2019).
  * Nair and Agarwal (2020) J. M. Nair and G. Agarwal, Applied Physics Letters 117, 084001 (2020).
  * Wang  _et al._ (2018) Y.-P. Wang, G.-Q. Zhang, D. Zhang, T.-F. Li, C.-M. Hu, and J. Q. You, Phys. Rev. Lett. 120, 057202 (2018).
  * Lachance-Quirion  _et al._ (2019) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Applied Physics Express 12, 070101 (2019).
  * Nair  _et al._ (2020) J. M. P. Nair, Z. Zhang, M. O. Scully, and G. S. Agarwal, Phys. Rev. B 102, 104415 (2020).
  * Nair  _et al._ (2021a) J. M. P. Nair, D. Mukhopadhyay, and G. S. Agarwal, Phys. Rev. B 103, 224401 (2021a).
  * Khitrova  _et al._ (2006) G. Khitrova, H. Gibbs, M. Kira, S. W. Koch, and A. Scherer, Nature physics 2, 81 (2006).
  * Agarwal (1984) G. S. Agarwal, Phys. Rev. Lett. 53, 1732 (1984).
  * Harder  _et al._ (2018) M. Harder, Y. Yang, B. M. Yao, C. H. Yu, J. W. Rao, Y. S. Gui, R. L. Stamps, and C.-M. Hu, Phys. Rev. Lett. 121, 137203 (2018).
  * Yao  _et al._ (2019) B. Yao, T. Yu, X. Zhang, W. Lu, Y. Gui, C.-M. Hu, and Y. M. Blanter, Phys. Rev. B 100, 214426 (2019).
  * Bhoi  _et al._ (2019) B. Bhoi, B. Kim, S.-H. Jang, J. Kim, J. Yang, Y.-J. Cho, and S.-K. Kim, Phys. Rev. B 99, 134426 (2019).
  * Yang  _et al._ (2019) Y. Yang, J. Rao, Y. Gui, B. Yao, W. Lu, and C.-M. Hu, Phys. Rev. Applied 11, 054023 (2019).
  * Rao  _et al._ (2019) J. Rao, C. Yu, Y. Zhao, Y. Gui, X. Fan, D. Xue, and C. Hu, New Journal of Physics 21, 065001 (2019).
  * Xu  _et al._ (2019) P.-C. Xu, J. W. Rao, Y. S. Gui, X. Jin, and C.-M. Hu, Phys. Rev. B 100, 094415 (2019).
  * Wang  _et al._ (2019) Y.-P. Wang, J. W. Rao, Y. Yang, P.-C. Xu, Y. S. Gui, B. M. Yao, J. Q. You, and C.-M. Hu, Phys. Rev. Lett. 123, 127202 (2019).
  * Lachance-Quirion  _et al._ (2017) D. Lachance-Quirion, Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, and Y. Nakamura, Science advances 3, e1603150 (2017).
  * Maier-Flaig  _et al._ (2016) H. Maier-Flaig, M. Harder, R. Gross, H. Huebl, and S. T. B. Goennenwein, Phys. Rev. B 94, 054433 (2016).
  * Lambert  _et al._ (2016) N. J. Lambert, J. A. Haigh, S. Langenfeld, A. C. Doherty, and A. J. Ferguson, Phys. Rev. A 93, 021803 (2016).
  * Yao  _et al._ (2017) B. Yao, Y. Gui, J. Rao, S. Kaur, X. Chen, W. Lu, Y. Xiao, H. Guo, K.-P. Marzlin, and C.-M. Hu, Nature communications 8, 1 (2017).
  * Nair  _et al._ (2021b) J. M. P. Nair, D. Mukhopadhyay, and G. S. Agarwal, Phys. Rev. Lett. 126, 180401 (2021b).
  * Hisatomi  _et al._ (2016) R. Hisatomi, A. Osada, Y. Tabuchi, T. Ishikawa, A. Noguchi, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. B 93, 174427 (2016).
  * Ihn  _et al._ (2020) Y. S. Ihn, S.-Y. Lee, D. Kim, S. H. Yim, and Z. Kim, Phys. Rev. B 102, 064418 (2020).
  * Mukhopadhyay  _et al._ (2021) D. Mukhopadhyay, J. M. Nair, and G. S. Agarwal, arXiv preprint arXiv:2111.01335 (2021).
  * Dembowski  _et al._ (2001) C. Dembowski, H.-D. Gräf, H. Harney, A. Heine, W. Heiss, H. Rehfeld, and A. Richter, Physical review letters 86, 787 (2001).
  * 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 537, 76 (2016).
  * Xu  _et al._ (2016) H. Xu, D. Mason, L. Jiang, and J. Harris, Nature 537, 80 (2016).
  * Wiersig (2014) J. Wiersig, Phys. Rev. Lett. 112, 203901 (2014).
  * Chen  _et al._ (2017) W. Chen, Ş. K. Özdemir, G. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Physics 14, 11 (2018).
  * Özdemir  _et al._ (2019) Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nature materials 18, 783 (2019).
  * Miri and Alu (2019) M.-A. Miri and A. Alu, Science 363 (2019).
  * Zhang  _et al._ (2020) H. Zhang, R. Huang, S.-D. Zhang, Y. Li, C.-W. Qiu, F. Nori, and H. Jing, Nano Letters 20, 7594 (2020).
  * Bernier  _et al._ (2018) N. R. Bernier, L. D. Tóth, A. K. Feofanov, and T. J. Kippenberg, Phys. Rev. A 98, 023841 (2018).
  * Yang  _et al._ (2020) Y. Yang, Y.-P. Wang, J. W. Rao, Y. S. Gui, B. M. Yao, W. Lu, and C.-M. Hu, Phys. Rev. Lett. 125, 147202 (2020).
  * Agarwal (1974) G. S. Agarwal, in  _Quantum Optics_ (Springer, 1974) pp. 1–128.
  * Agarwal  _et al._ (1997) G. S. Agarwal, R. R. Puri, and R. P. Singh, Phys. Rev. A 56, 2249 (1997).
  * Yuen-Zhou and Menon (2019) J. Yuen-Zhou and V. M. Menon, Proceedings of the National Academy of Sciences 116, 5214 (2019).
  * Agarwal (2012) G. S. Agarwal,  _Quantum optics_ (Cambridge University Press, 2012).