# Contact solitons as spin pistons in one-dimensional ferromagnetic channels

Medhanie Estiphanos  Ezio Iacocca  Center for Magnetism and Magnetic
Nanostructures, University of Colorado Colorado Springs, Colorado Springs, CO,
USA

###### Abstract

Ferromagnetic channels subject to spin injection have been theoretically shown
to sustain dissipative exchange flows (DEFs). In the strong injection regime,
a soliton is stabilized at the injection site, which has been termed a contact
soliton DEF or CS-DEF. Here, we investigate the modulation of CS-DEFs as a
mechanism to inject magnons into a DEF. By varying the injection current, the
parameters connecting the soliton with the DEF via a boundary layer are
varied. These lead to a modification in the DEF which is interpreted as a spin
piston. The injected magnons follow the expected dispersion relation of DEFs,
akin to the Bogoliubov - de Gennes dispersion relation. This work demonstrates
that changes in the injected current can pumping magnons along a ferromagnetic
channel via dissipative exchange flows.

##  I Introduction

One of the challenges of information storage, transport, and processing is
energy consumption. The field of magnonics [1] offers new possibilities for
remedying this problem by proposing low-dissipation information transport
mediated by the quanta of magnetic excitation, known as a magnon.
Fundamentally, each magnon has angular momentum $\hbar$ and thus can transport
spin with little Joule heating [2]. However, magnons decay exponentially with
the magnetic damping parameter $\alpha$ [3] so that any benefit from magnonic
devices is closely linked to the material used. Typically, permalloy (Py,
Ni80Fe20) can be patterned in films of $5$ nm to $10$ nm in thickness and has
a magnon decay length of about $1$ $\mu$m [4]; while Yttrium Iron garnet
(YIG), a magnetic insulator, has a much longer decay that can reach
centimeters for samples as thick as 200 $\mu$m [5]. Such long magnon decays
have made insulators an attractive material system for magnonics, including
amorphous Yttrium-Iron ferrites [6] and antiferromagnetic hematites [7].

To beat the magnon exponential decay, an alternative was theoretically
proposed as spin superfluids [8, 9, 10, 11]. These are large-amplitude
dynamical structures that appear as solutions to spin injection [12] in an in-
plane ferromagnet. While magnetic damping is still active, in contrast to
typical superfluids, spin superfluids exhibit an algebraic decay along the
full length of the ferromagnet. In other words, spin superfluids are solutions
to a boundary value problem (BVP) and their algebraic decay is enabled by the
spatial chirality of the structure.

The theoretical investigation of spin superfluids is conveniently achieved
through a fluid interpretation of magnetization dynamics [13, 10]. In fact,
this interpretation demonstrates the similarity between the equations of mass
and spin superfluidity. An exact dispersive spin hydrodynamics formulation was
derived later [14, 15, 16, 17, 18], which also allowed to generalize the BVP
solutions into dissipative exchange flows (DEFs). Several theoretical
investigations have further demonstrated that DEFs are general solutions
exhibiting nonlinear character [17, 19, 20, 21, 22] and can be stabilized in
the presence of anisotropy [16]. In addition, DEFs are theoretically stable
and can only be annihilated via vortex shedding [23, 15, 24]. Despite the
analytical and theoretical predictions, the experimental demonstration of DEFs
has been challenging [25, 26] and it can be argued that it remains
inconclusive. However, the investigation of hydrodynamic modes in magnetism
remains of theoretical interests and its use for interpretation of, e.g.,
ultrafast magnetic phenomena [27, 28] is still in its infancy.

In general, DEFs have been investigated in the low-injection regime due to
possible experimental constraints to inject spin via spin-Hall effect [29]. In
that regime, DEFs can be approximated precisely as spin superfluids [9, 10].
In the language of spin hydrodynamics, the injection conditions are subsonic
in this case, implying a laminar spin injection. In the high-injection regime,
the injection conditions become supersonic which would suggest the onset of
turbulence. However, the strong exchange in magnetic materials leads to the
formation of a stationary soliton at the injection site. This soliton prevents
turbulence and smoothly connects with a DEF in the remainder of the channel.
This solution was termed a contact-soliton DEF (CS-DEF) [17]. In contrast to
DEFs, the frequency exhibits a negative tunability to spin injection, which is
understood as the soliton becoming larger in amplitude. CS-DEFs and their
properties remain largely unexplored.

In this paper, we numerically investigate the time-dependent modulation of CS-
DEFs where the soliton acts as a spin piston. The modulation leads to
modifications of the parameters connecting the soliton and the DEF solution,
which are defined by a boundary layer approach. This implies that the soliton
per se does not “inject” magnons but its profile modification leads to a
varying spin density and fluid velocity of the DEF. Both continuum and pseudo-
spectral models are used to determine whether magnons are injected. We find
that the detailed soliton profile is different between the models due to their
different dispersion relation, but both models describe the injected magnons
well.

The remainder of the paper is organized as follows. In section II, we
summarize the analytical solution of a CS-DEF and explore its limit to
infinite injection. Modulation to achieve a spin piston is discussed in
section III. We consider the case where the injection is located in the middle
of the channel and with finite a width. The results are compared with
simulations where the magnon energy is accurately modeled by a pseudospectral
Landau-Lifshitz representation in section IV and demonstrate that the wave
injection is not dependent on the particular soliton profile. Concluding
remarks are provided in section V.

##  II Contact-soliton dissipative exchange flows

The equation of motion for a one-dimensional (1D) ferromagnetic chain is
described by the Landau-Lifshitz (LL) equation

$\partial_{t}\mathbf{m}=-[\mathbf{m}\times\mathbf{h}_{\mathrm{eff}}+\alpha\mathbf{m}\times(\mathbf{m}\times\mathbf{h}_{\mathrm{eff}})]$
(1)

expressed here in dimensionless form, such that the magnetization vector
$\mathbf{m}$ has unit magnitude and the effective field
$\mathbf{h}_{\mathrm{eff}}=\partial_{xx}\mathbf{m}-m_{z}$ contains,
respectively, exchange and uniaxial anisotropy (or local dipole) fields. To
obtain this representation [14], time is scaled by $\gamma\mu_{0}M_{s}$, where
$\gamma$ is the gyromagnetic ratio, $\mu_{0}$ is the vacuum permeability, and
$M_{s}$ is the saturation magnetization; and space is scaled by the exchange
length $\lambda_{\mathrm{ex}}$. Energy dissipation is quantified by the
Gilbert damping constant which can be used in Eq. (1) if $\alpha\ll 1$, which
will be assumed to hold for in-plane soft ferromagnets.

Equation (1) can be recast in dispersive hydrodynamic form [14, 18] by
introducing a spin density $n=m_{z}$ and a fluid velocity
$u=-\partial_{x}\Phi$, with $\Phi=\mathrm{atan}(m_{y}/m_{x})$ being the in-
plane magnetization’s phase. This results in the hydrodynamic equations

$\displaystyle\partial_{t}n$ $\displaystyle=$
$\displaystyle\partial_{x}[(1-n^{2})u]+\alpha(1-n^{2})\partial_{t}\Phi,$ (2a)
$\displaystyle\partial_{t}u$ $\displaystyle=$
$\displaystyle-\partial_{x}[(1+u^{2})n]-\left[\frac{\partial_{xx}n}{{1-n^{2}}}+\frac{(\partial_{x}n)^{2}}{(1-n^{2})^{2}}\right]$
(2b)
$\displaystyle+\alpha\partial_{x}\left[\frac{1}{1-n^{2}}\partial_{x}[(1-n^{2})u]\right].$

DEFs are solutions to Eqs. (2) in a channel of length $L$ when subject to
boundary conditions

$\displaystyle\partial_{t}n(0)=0,$ $\displaystyle\partial_{t}n(L)=0,$ (3a)
$\displaystyle u(0)=\bar{u},$ $\displaystyle u(L)=0,$ (3b)

so that $\bar{u}$ is a fluid injection strength.

The analytical solutions were presented in Ref. [17]. Here, we present the
contact soliton solution for completeness. This is a case where $\bar{u}$ is
sufficiently large and so the equations can be separated in two regions: an
inner supersonic injection region where exchange dominates and an outer
subsonic channel region where dissipation dominates. Using the boundary layer
method, a smooth solution is found by matching both regions, provided a
uniform precessional frequency $\Omega$ along the channel. The solution is

$\displaystyle n_{cs}$ $\displaystyle=$ $\displaystyle
n_{in}(x)+n_{out}(x/L)-n_{\infty},$ (4a) $\displaystyle u_{cs}$
$\displaystyle=$ $\displaystyle u_{in}(x)+u_{out}(x/L)-u_{\infty}.$ (4b)

In this solution, a soliton exist in the inner region, given by

$\displaystyle n_{in}$ $\displaystyle=$
$\displaystyle\frac{a\nu_{1}\tanh^{2}(\theta
x)+\nu_{2}(n_{\infty}-a)}{a\tanh^{2}(\theta x)+\nu_{2}},$ (5a) $\displaystyle
u_{in}$ $\displaystyle=$ $\displaystyle
u_{\infty}\frac{1-n_{\infty}^{2}}{1-n_{in}^{2}},$ (5b)

where ${\nu_{1}=a-n_{\infty}-2n_{\infty}u_{\infty}^{2}}$,
${\nu_{2}=a-2n_{\infty}-2n_{\infty}u_{\infty}^{2}}$,
${\theta=\sqrt{1-u_{\infty}^{2}-n_{\infty}^{2}(1+3u_{\infty}^{2})}}$, and
${a=n_{\infty}(1+u_{\infty}^{2})+\sqrt{(1-u_{\infty}^{2})(1-n_{\infty}^{2}u_{\infty}^{2})}}$.

The outer solution is, in general, a nonlinear DEF given by

$\displaystyle n_{out}$ $\displaystyle=$
$\displaystyle-\frac{\Omega}{1-u_{out}^{2}},$ (6a) $\displaystyle\alpha
L\Omega$ $\displaystyle=$ $\displaystyle u_{out}+4\tanh^{-1}(u_{out})$
$\displaystyle-2\left[\mathcal{N}^{-}(u_{out},\Omega)+\mathcal{N}^{+}(u_{out},\Omega)\right],$

where $u_{out}$ must be numerically solved from the transcendental Eq. (6)
given the functions

$\mathcal{N}^{\pm}(\kappa,\omega)=\sqrt{1\pm\omega}\tanh^{-1}\left(\frac{\kappa}{\sqrt{1\pm\omega}}\right).$
(7)

Finally, the matching conditions relate the fluid injection strength to the
solution

$\displaystyle n_{\infty}$ $\displaystyle=$
$\displaystyle-\frac{\Omega}{1-u_{\infty}^{2}},$ (8a) $\displaystyle\bar{u}$
$\displaystyle=$
$\displaystyle\frac{u_{\infty}(1-n_{\infty}^{2})}{1-(n_{\infty}-a)^{2}},$ (8b)

Details on the derivation of these equations are found in Ref. [17].

In the limit where $\bar{u}\rightarrow\infty$, the matching conditions
approach a linear DEF or spin superfluid solution so that
$n_{\infty}\rightarrow 0$. In turn, we can simplify some of the coefficients
in the soliton solution: $a\approx 1$, $\nu_{1}\approx 1$, $\nu_{2}\approx 1$,
and $\theta\approx 1$ as well as the precessional frequency $\Omega\approx
u_{\infty}/(\alpha L)=-n_{\infty}$. Notably, the soliton solution simplifies
to

$\lim_{\bar{u}\rightarrow\infty}n_{in}(x)=\frac{\tanh^{2}(x)-1}{\tanh^{2}(x)+1}.$
(9)

In other words, the soliton has a profile independent of $\bar{u}$ implying
that the injection is fully invested in the exchange energy. This can be also
appreciated from the boundary layer width $l$, Eq. (26) in Ref. [17], given by

$l=\frac{1}{\theta}\mathrm{atanh}{\left(\sqrt{\frac{\nu_{2}}{2(\nu_{1}-n_{\infty})+a}}\right)}\approx\mathrm{atanh}{\left(\sqrt{\frac{1}{3}}\right)}$
(10)

Figure 1:  Evolution of (a) the spin density at the injection site
$n_{in}(x=0)$, (b) $n_{\infty}$, (c), $u_{\infty}$, and (d) the boundary layer
width $l$ as a function of the spin injection strength $\bar{u}$. We observe
that the soliton solution is obtained when $n_{in}(x=0)=1$, roughly at
$\bar{u}=6$ in panel (a). At this point, the boundary layer width approach its
minimal value reported in Eq. (10) and shown by a black dashed line in (d).
While the soliton is apparently fully formed, there is still a DEF at
$\bar{u}\approx 10$. This DEF appears suppressed at about an order of
magnitude larger, $\bar{u}\approx 100$. The dotted gray vertical lines
indicate the equilibrium condition used for numerical calculations.

Several parameters of the CS-DEF solution are shown in Fig. 1 as a function of
$\bar{u}$. The spin density at the injection site, $n_{in}(x=0)$ is shown in
panel (a). In our geometry, injection leads to a negative tilt in the
magnetization so that $n_{in}(x=0)=-1$ in Eq. (9). This value is approximately
achieved at $\bar{u}\approx 10$. The matching conditions $n_{\infty}$ and
$u_{\infty}$ are shown in panels (b) and (c), respectively. Both approach zero
as $\bar{u}$ increases. Interestingly, the vanishing DEF is achieved at much
larger injection of $\bar{u}\approx 100$. This means that while spin density
at the injection site is in its vacuum state (out-of-plane magnetization), the
soliton has an amplitude smaller than one leading to a non-negligible
background. Finally, the boundary layer is shown in panel (d). Similar to (a),
the boundary layer reaches its smallest extent at $\bar{u}\approx 10$, shown
by a black dashed line calculated from Eq. (10). We note that
$\mathrm{atanh}\sqrt{1/3}\approx 0.66$, which is in units of exchange lengths.
For typical permalloy parameters, $\lambda_{\mathrm{ex}}=5$ nm, so that the
soliton is approximately 3.3 nm wide. Such a width is close to the region
where the micromagnetic approximation ceases to be accurate [30]. The
implications of such a soliton width is numerically investigated in section
IV.

Our goal is to modulate the spin injection to pump magnons into the channel
via the modification of the contact soliton. We set the injection to
$\bar{u}=1$ so that both $u_{\infty}$ and $n_{\infty}$ can be periodically
changed. The hypothesis is that such a change in DEF injection conditions can
give rise to a magnon propagating away from the injection region according to
the dispersion relation of magnons [14]

$\omega=2n_{\infty}u_{\infty}k\pm
k\sqrt{(1-n_{\infty}^{2})(1-u_{\infty}^{2})+k^{2}}$ (11)

where $k$ is the wavenumber and we assume that the DEF decays slowly and can
be approximated as a uniform solution, as done in Ref. [17].

##  III Spin piston

The numerical investigation of spin injection is performed by solving Eq. (1)
using a complex representation of the magnetization, $m_{z}$ and
$m_{p}=m_{x}-im_{y}$. This approach has been used in previous publications
[16, 15, 17]. We use a channel of length $L=2000$ discretized in cells of
$0.5$ and a material with $\alpha=0.01$. For permalloy parameters, $M_{s}=790$
kA/m and $\lambda_{\mathrm{ex}}=5$ nm, the channel would have a length of $10$
$\mu$m discretized at $2.5$ nm.

We consider injection in the middle of the channel, similar to Ref. [19] which
is a potentially more realistic geometry to investigate experimentally. The
spin-torque injection is implemented as [16]

$\tau=\bar{\mu}\mathbf{m}\times\left(\mathbf{m}\times\mathbf{p}\right),$ (12)

where $\bar{\mu}$ is a normalized injection strength and $\mathbf{p}=\hat{z}$
is the spin polarization direction. This spin polarization is possible, e.g.,
by a charge current passed through a perpendicular fixed layer. The injection
strengths $\bar{\mu}$ and $\bar{u}$ are related through geometric factors
[16]. This is because $\bar{\mu}$ is applied on a top surface while $\bar{u}$
considered lateral injection. We consider an injection of $\bar{\mu}=1$
applied in a region of size $11$. Again using permalloy parameters, this would
correspond to a lateral length of $27.5$ nm. The stabilized CS-DEF solution in
the vicinity of the injection region is shown in Fig. 2, while the complete
simulation domain is shown in the inset. The injection region is depicted by
the gray area. A single soliton is stabilized and the reminder of the channel
has effectively a linear DEF with $n_{\infty}\approx-0.02$ in the reminder of
the channel.

Figure 2:  CS-DEF solution stabilized with $\bar{\mu}=1$, shown by a solid
blue curve. The injection region is shown by a gray area. The analytical
solution of Eqs. (4) is shown by a dotted red curve. The full solution is
shown in the inset.

Comparison with the analytical solution Eqs. (4) is not directly possible
because those expressions are valid only for lateral injection. However, we
find reasonable agreement with the stabilized DEF when $\bar{u}\approx 4.4$,
shown in Fig. 2 by a dotted red curve. The soliton is clearly steeper for the
analytical solution because of the injection applied only at the boundary.
This condition is shown in Fig. 1 as vertical dotted gray lines. Because the
soliton is almost fully formed, this is a quasi linear regime for the
modulation of $n_{\infty}$ and $u_{\infty}$.

Modulation of the spin injection is implemented as a time-dependent torque

$\mu(t)=\bar{\mu}\left[1+A\sin{(2\pi ft)}\right].$ (13)

Both amplitude $A$ and frequency $f$ can be tuned. First, we vary the
amplitude to ensure that the resulting modulation is in a linear regime, i.e.
when the slope of $n_{\infty}$ with respect to $\bar{u}$ is approximately
linear in Fig. 1. For these simulations, we used a frequency of $f=0.015$
which corresponds to $420$ MHz. An example solution with $A=0.5$ is shown in
Fig. 3. The injection modulation indeed leads to the excitation of waves
propagating and decaying outwards from the soliton. The inset shows the
relationship between the amplitude and the maximum amplitude of $m_{z}$
relative to the background. This is a clearly linear relationship so that we
are operating in the linear regime.

Then, we set the amplitude $A=0.5$ and vary the modulation frequency from
0.00357 to 0.0357 (100 MHz to 1 GHz) in steps of 0.00178 ($50$ MHz). To
determine the frequency and wavenumber of the pumped magnon, we analyzed the
data in Fourier space. To estimate the frequency, we find the time trace for
the numerical cells in the regions 50 to 150 cells away from the injection
site, both to the right and left since the modulation must be symmetric. The
maximum peak of the Fourier transform for each cell was found and averaged to
estimate the resulting oscillation frequency. As expected, this frequency is
directly proportional to the modulation frequency, i.e., $\omega=2\pi f$,
estimated to an error of $1.5\times 10^{-4}$ (4.3 MHz for permalloy). The
wavenumber was estimated from a snapshot of the steady-state solution. We
consider a region 50 to 450 cells away from the injection point, both at its
right and left. First, a second-order polynomial fit is performed to subtract
the nonlinear DEF solution. The remaining signal is Fourier transformed and
its maximum peak corresponds to the wavenumber. Averaging the computed
wavenumber for both sides of the injection, we achieve an error of 0.013
($2.5~{}\mu$m-1 for permalloy).

Figure 3:  Modulated CS-DEF with $A=0.5$ and $f=0.015$. The modulation leads
to the excitation of propagating waves. The deviation of $m_{z}$ as a function
of the amplitude $A$ is shown in the inset. The linear behavior indicates that
we are operating in a linear regime.

The resulting dispersion relation for the excited magnons is shown in Fig. 4
by blue circles with errorbars. The analytical dispersion according Eq. (11)
is shown by a solid red curve. We find very good agreement between the
expected dispersion and the numerical simulations, within errorbars. Note that
the analytical solution assumes a uniform fluid mode rather than a DEF. In
this case, the agreement is due to the small $n_{\infty}$ so that the change
in fluid parameters along the DEF is small. This results demonstrates that we
are indeed pumping magnons into the DEF via the modulation of the contact
soliton, particularly of the matching conditions $n_{\infty}$ and
$u_{\infty}$.

Figure 4: Dispersion relation of the magnons pumped by the spin piston, shown
in blue circles with errorbars. The analytical dispersion estimated from a
uniform solution, Eq. (11) is shown by a solid red curve.

##  IV Comparison to pseudospectral LL model

The hydrodynamic formulation of magnetization dynamics is rooted in the
micromagnetic expansion of the Heisenberg exchange Hamiltonian. This means
that the exchange takes a Laplacian form which allows for the use of vector
identities to derive continuum solutions. For the average injection used,
above $\bar{\mu}=1$, the boundary layer width is approximately $1.5$ which
corresponds to a soliton with a width of $3$ exchange lengths. While it is
possible to increase the resolution numerically by using small cell-sizes, it
is important to note that the solution will obey the micromagnetic limit
regardless.

To circumvent this issue, a pseudospectral LL (PS-LL) model has been recently
introduced [30]. In this model, the exchange interaction is written as a
convolution kernel where energy and momentum are described solely by the
magnon dispersion relation. The PS-LL model was found to be in excellent
agreement with atomistic spin dynamic simulations, yet utilizing a continuum
numerical and analytical formulation. Here, we use the PS-LL model extended to
include spin-transfer torques to compare with the continuum results.

Because spin-transfer torque is a local term, its implementation in the PS-LL
model is completely analogous to the LL equation, as described elsewhere [16].
For the PS-LL model, we use permalloy parameters directly. This implies that
the torque must be scaled accordingly.

A Slonczewski-type torque [31] is given by

$\tau=\frac{\gamma\hbar P}{M_{s}|e|\delta}J\epsilon,$ (14)

where $\hbar$ is Planck’s constant, $P$ is the spin polarization, $e$ is the
electron’s charge, $J$ is the current density, $\epsilon$ is the spin
asymmetry, and $\delta$ is the layer’s thickness, assumed to be identical to
the cell size. The torque of Eq. (14) is in units of Hz and so it can be
rendered dimensionless by scaling to $\gamma\mu_{0}M_{s}$. Using the same
geometric factor, we adjust for the use of a cell of $2.5$ nm to directly
compare with micromagnetic simulations. This results in a forcing term
$\tau\approx 1.34$ THz so that $J$ is on the order of $10^{14}$ A/m2 when
considering a symmetric torque $\epsilon=1/2$ and spin polarization $P=0.65$.
This is clearly a large current density that is unlikely to be delivered
without material damage. However, other forms of spin injection can be more
efficient, e.g., spin-Hall effect [29], and lead to the use of physically
realizable charge currents in the future.

Figure 5:  CS-DEF solution stabilized with the LL scaled to exchange length
(solid blue curve) and PS-LL (dotted black curve) models. The injection region
is shown by a gray area. The full solution is shown in the inset.

The relaxed CS-DEF is shown in Fig. 5. A immediate difference compared to Fig.
2 is that the soliton is wider. This can be understood as the exchange energy
is smaller in the PS-LL model because the magnon dispersion relation is, in
fact, bounded to the first-Brillouin zone. The inset shows the full solution,
with the reminder of the channel in good agreement between both models. The
difference originates from the exchange energy required in each model. The
exchange energy in the micromagnetic model (solid blue curve) scales as
$k^{2}$ so that sharp transitions require a significant torque. In the PS-LL
model (dashed black curve), the exchange energy scales with $1-\cos{(ka)}$,
with $a$ the atomic lattice, so that the energy is bounded. Even though we are
using a cell-size of $2.5$ nm, this is sufficient to account for the fact that
the magnetization is almost uniformly reversed at the injection site and the
soliton is defined outside of this region. The same soliton is stabilized with
atomic resolution, as shown in Appendix A. This result already indicates that
it is possible to sustain a soliton at lower currents, but we keep the large
current density here to compare with the continuum results.

Modulation in the current injection is applied with $A=0.5$ and the frequency
is varied in the same range as presented before. However, we raise the
frequency here in steps of $20$ MHz. This is because the PS-LL simulations are
much faster so that the time and storage resources needed are manageable. The
results are shown in Fig. 6 and are in good agreement with those shown in Fig.
4. In particular, the modulation of current leads to spin injection due to the
modulation of the $n_{\infty}$ and $u_{\infty}$ parameters and the dispersion
is in excellent agreement with the dispersion relation, corroborating that
these are magnons on top of a DEF. Use of the PS-LL model also confirms that
the modulation of the soliton is independent of the particularities of its
profile.

Figure 6: Dispersion relation of the magnons pumped by the spin piston, shown
in blue circles with errorbars. The analytical dispersion estimated from a
uniform solution, Eq. (11) is shown by a solid red curve.

##  V Conclusions

We investigated the possibility to inject magnons into a DEF by means of
modulating a contact soliton, which we consider an analog to a spin piston.
The injection is achieved by modulating the boundary layer parameters, i.e.,
the conditions linking the contact soliton solution at the injection site and
the DEF defined throughout the channel. Both LL and PS-LL approaches were used
to corroborate that injection occurs by this process. In addition, the process
is shown to be independent of the particularities of the soliton profile,
which are different between the LL and PS-LL descriptions. Nevertheless, it is
important to notice that the soliton profile depends on the short-wavelength
energies of the ferromagnetic channel which are accurately modeled by the PS-
LL model [30].

Magnon injection is in good agreement with the predicted dispersion relation
from a uniform hydrodynamic flow, which is a approximately applicable insofar
as the DEF is a negligible change in its spin density and fluid velocity. This
serves as proof that the injected waves are magnons on top of a DEF, following
a linear dispersion rather than the typical $k^{2}$ dispersion of magnons in
uniformly magnetized ferromagnets.

The currents needed to excite a CS-DEF via spin transfer torque are not
realistic, alternative spin injection methods could stabilize DEFs. Because
the magnons excited on a DEF follow a Bogoliubov - de Gennes dispersion
relation, detection of such magnons would provide an additional method to
verify the existence of a DEF.

## Acknowledgments

This work was supported by the U.S. Department of Energy, Office of Basic
Energy Sciences under Award Number DE-SC0024339.

Figure 7: CS-DEF solution stabilized with the PS-LL model using a
discretization of 2.5 nm (dotted black curve) and 0.355 nm (solid red curve)

.

##  Appendix A PS-LL simulations with atomic resolution

Here, we stabilize a soliton using PS-LL with atomic resolution, $a=0.355$ for
permalloy. Because this is a 1D simulation, we also need to adjust the torque
to take into account the new geometric ratio. The stabilized soliton is shown
in Fig. 7 by a solid red curve. The soliton stabilized with a cell size of
$2.5$ nm is shown by a dashed black curve (also shown in Fig. 5. There are
clear discrepancies in the stabilized soliton but in general both
discretization show that the injection region is primarily saturated and the
soliton forms away from the injection region.

## References

  * Chumak  _et al._ [2022] A. V. Chumak, P. Kabos, M. Wu, C. Abert, C. Adelmann, A. O. Adeyeye, J. Åkerman, F. G. Aliev, A. Anane, A. Awad, C. H. Back, A. Barman, G. E. W. Bauer, M. Becherer, E. N. Beginin, V. A. S. V. Bittencourt, Y. M. Blanter, P. Bortolotti, I. Boventer, D. A. Bozhko, S. A. Bunyaev, J. J. Carmiggelt, R. R. Cheenikundil, F. Ciubotaru, S. Cotofana, G. Csaba, O. V. Dobrovolskiy, C. Dubs, M. Elyasi, K. G. Fripp, H. Fulara, I. A. Golovchanskiy, C. Gonzalez-Ballestero, P. Graczyk, D. Grundler, P. Gruszecki, G. Gubbiotti, K. Guslienko, A. Haldar, S. Hamdioui, R. Hertel, B. Hillebrands, T. Hioki, A. Houshang, C.-M. Hu, H. Huebl, M. Huth, E. Iacocca, M. B. Jungfleisch, G. N. Kakazei, A. Khitun, R. Khymyn, T. Kikkawa, M. Kläui, O. Klein, J. W. Kłlos, S. Knauer, S. Koraltan, M. Kostylev, M. Krawczyk, I. N. Krivorotov, V. V. Kruglyak, D. Lachance-Quirion, S. Ladak, R. Lebrun, Y. Li, M. Lindner, R. Macêdo, S. Mayr, G. A. Melkov, S. Mieszczak, Y. Nakamura, H. T. Nembach, A. A. Nikitin, S. A. Nikitov, V. Novosad, J. A. Otálora, Y. Otani, A. Papp, B. Pigeau, P. Pirro, W. Porod, F. Porrati, H. Qin, B. Rana, T. Reimann, F. Riente, O. Romero-Isart, A. Ross, A. V. Sadovnikov, A. R. Safin, E. Saitoh, G. Schmidt, H. Schultheiss, K. Schultheiss, A. A. Serga, S. Sharma, J. M. Shaw, D. Suess, O. Surzhenko, K. Szulc, T. Taniguchi, M. Urbánek, K. Usami, A. B. Ustinov, T. van der Sar, S. van Dijken, V. I. Vasyuchka, R. Verba, S. V. Kusminskiy, Q. Wang, M. Weides, M. Weiler, S. Wintz, S. P. Wolski, and X. Zhang, Advances in magnetics roadmap on spin-wave computing, IEEE Transactions on Magnetics 58, 1 (2022).
  * Kruglyak  _et al._ [2010] V. V. Kruglyak, S. O. Demokritov, and D. Grundler, Magnonics, J. Phys. D Appl. Phys. 43, 264001 (2010).
  * Madami  _et al._ [2011a] M. Madami, S. Bonetti, G. Consolo, S. Tacchi, G. Carlotti, G. Gubbiotti, F. Mancoff, M. A. Yar, and J. Åkerman, Direct observation of a propagating spin wave induced by spin-transfer torque, Nature nanotechnology 6, 635 (2011a).
  * Madami  _et al._ [2011b] M. Madami, S. Bonetti, G. Consolo, S. Tacchi, G. Carlotti, G. Gubbiotti, F. B. Mancoff, M. A. Yar, and J. Åkerman, Direct observation of a propagating spin wave induced by spin-transfer torque, Nat. Nanotechnol. 6, 635 (2011b).
  * Serga  _et al._ [2010] A. A. Serga, A. V. Chumak, and B. Hillebrands, Yig magnonics, J. Phys. D Appl. Phys. 43, 264002 (2010).
  * Wesenberg  _et al._ [2017] D. Wesenberg, T. Liu, D. Balzar, M. Wu, and B. L. Zink, Long-distance spin transport in a disordered magnetic insulator, Nature Physics 13, 987 (2017).
  * Lebrun  _et al._ [2018] R. Lebrun, A. Ross, S. A. Bender, A. Qaiumzadeh, L. Baldrati, J. Cramer, A. Brataas, R. A. Duine, and M. Kläui, Tunable long-distance spin transport in crystalline antiferromagnetic iron oxide, Nature 561, 222 (2018).
  * König  _et al._ [2001] J. König, M. C. Bønsager, and A. H. MacDonald, Dissipationless spin transport in thin film ferromagnets, Phys. Rev. Lett. 87, 187202 (2001).
  * Sonin [2010] E. Sonin, Spin currents and spin superfluidity, Advances in Physics 59, 181 (2010).
  * Takei and Tserkovnyak [2014] S. Takei and Y. Tserkovnyak, Superfluid spin transport through easy-plane ferromagnetic insulators, Physical review letters 112, 227201 (2014).
  * Chen  _et al._ [2014] H. Chen, A. D. Kent, A. H. MacDonald, and I. Sodemann, Nonlocal transport mediated by spin supercurrents, Phys. Rev. B 90, 220401 (2014).
  * Tserkovnyak  _et al._ [2002] Y. Tserkovnyak, A. Brataas, and G. E. Bauer, Enhanced gilbert damping in thin ferromagnetic films, Physical review letters 88, 117601 (2002).
  * Halperin and Hohenberg [1969] B. Halperin and P. Hohenberg, Hydrodynamic theory of spin waves, Physical Review 188, 898 (1969).
  * Iacocca  _et al._ [2017a] E. Iacocca, T. J. Silva, and M. A. Hoefer, Breaking of galilean invariance in the hydrodynamic formulation of ferromagnetic thin films, Physical Review Letters 118, 017203 (2017a).
  * Iacocca and Hoefer [2017] E. Iacocca and M. A. Hoefer, Vortex-antivortex proliferation from an obstacle in thin film ferromagnets, Physical Review B 95, 134409 (2017).
  * Iacocca  _et al._ [2017b] E. Iacocca, T. J. Silva, and M. A. Hoefer, Symmetry-broken dissipative exchange flows in thin-film ferromagnets with in-plane anisotropy, Physical Review B 96, 134434 (2017b).
  * Iacocca and Hoefer [2019a] E. Iacocca and M. A. Hoefer, Hydrodynamic description of long-distance spin transport through noncollinear magnetization states: Role of dispersion, nonlinearity, and damping, Physical Review B 99, 184402 (2019a).
  * Iacocca and Hoefer [2019b] E. Iacocca and M. A. Hoefer, Perspectives on spin hydrodynamics in ferromagnetic materials, Physics Letters A 383, 125858 (2019b).
  * Schneider  _et al._ [2021] T. Schneider, D. Hill, A. Kákay, K. Lenz, J. Lindner, J. Fassbender, P. Upadhyaya, Y. Liu, K. Wang, Y. Tserkovnyak, I. N. Krivorotov, and I. Barsukov, Self-stabilizing exchange-mediated spin transport, Phys. Rev. B 103, 144412 (2021).
  * Smith  _et al._ [2021] D. A. Smith, S. Takei, B. Brann, L. Compton, F. Ramos-Diaz, M. J. Simmers, and S. Emori, Diffusive and fluidlike motion of homochiral domain walls in easy-plane magnetic strips, Phys. Rev. Appl. 16, 054002 (2021).
  * Hu  _et al._ [2021] M. Hu, E. Iacocca, and M. A. Hoefer, Spin-injection-generated shock waves and solitons in a ferromagnetic thin film, IEEE Transactions on Magnetics 58, 1300105 (2021).
  * Hu  _et al._ [2022] M. Hu, E. Iacocca, and M. A. Hoefer, Spin-piston problem for a ferromagnetic thin film: Shock waves and solitons, Phys. Rev. B 105, 104419 (2022).
  * Kim and Tserkovnyak [2016] S. K. Kim and Y. Tserkovnyak, Topological effects on quantum phase slips in superfluid spin transport, Phys. Rev. Lett. 116, 127201 (2016).
  * Iacocca [2020] E. Iacocca, Controllable vortex shedding from dissipative exchange flows in ferromagnetic channels, Phys. Rev. B 102, 224403 (2020).
  * Stepanov  _et al._ [2018] P. Stepanov, S. Che, D. Shcherbakov, J. Yang, R. Chen, K. Thilahar, G. Voigt, M. W. Bockrath, D. Smirnov, K. Watanabe, T. Taniguchi, R. K. Lake, Y. Barlas, A. H. MacDonald, and C. N. Lau, Long-distance spin transport through a graphene quantum hall antiferromagnet, Nature Physics 14, 907 (2018).
  * Yuan  _et al._ [2018] W. Yuan, Q. Zhu, T. Su, Y. Yao, W. Xing, Y. Chen, Y. Ma, X. Lin, J. Shi, R. Shindou, X. C. Xie, and W. Han, Experimental signatures of spin superfluid ground state in canted antiferromagnet cr2o3 via nonlocal spin transport, Science Advances 4, eaat1098 (2018).
  * Iacocca  _et al._ [2019] E. Iacocca, T.-M. Liu, A. H. Reid, Z. Fu, S. Ruta, P. W. Granitzka, E. Jal, S. Bonetti, A. X. Gray, C. E. Graves, R. Kukreja, Z. Chen, D. J. Higley, T. Chase, L. Le Guyader, K. Hirsch, H. Ohldag, W. F. Schlotter, G. L. Dakovski, G. Coslovich, M. C. Hoffmann, S. Carron, A. Tsukamoto, A. Kirilyuk, A. V. Kimel, T. Rasing, J. Stöhr, R. F. L. Evans, T. Ostler, R. W. Chantrell, M. A. Hoefer, T. J. Silva, , and H. A. Dürr, Spin-current-mediated rapid magnon localization and coalescence after ultrafast optical pumping of ferrimagnetic alloys, Nat. Commun. 10, 1756 (2019).
  * Jangid  _et al._ [2023] R. Jangid, N. Z. Hagström, M. Madhavi, K. Rockwell, J. M. Shaw, J. A. Brock, M. Pancaldi, D. De Angelis, F. Capotondi, E. Pedersoli, H. T. Nembach, M. W. Keller, S. Bonetti, E. E. Fullerton, E. Iacocca, R. Kukreja, and T. J. Silva, Extreme domain wall speeds under ultrafast optical excitation, Phys. Rev. Lett. 131, 256702 (2023).
  * Hoffmann [2013] A. Hoffmann, Spin hall effects in metals, IEEE transactions on magnetics 49, 5172 (2013).
  * Rockwell  _et al._ [2024] K. Rockwell, J. Hirst, T. A. Ostler, and E. Iacocca, Pseudospectral landau-lifshitz description of magnetization dynamics, Phys. Rev. B 109, L180404 (2024).
  * Slonczewski [1996] J. C. Slonczewski, Current-driven excitation of magnetic multilayers, Journal of Magnetism and Magnetic Materials 159, L1 (1996).


reprint, superscriptaddress, amsmath, amssymb, aps, prl, floatfix revtex4-2

# Emerging supersolidity from a polariton condensate in a photonic crystal
waveguide

Dimitrios Trypogeorgos<EMAIL_ADDRESS>CNR Nanotec,
Institute of Nanotechnology, via Monteroni, 73100, Lecce, Italy  Antonio
Gianfrate  CNR Nanotec, Institute of Nanotechnology, via Monteroni, 73100,
Lecce, Italy  Manuele Landini  Institut für Experimentalphysik und Zentrum für
Quantenphysik, Universität Innsbruck, 6020 Innsbruck, Austria  Davide Nigro
Dario Gerace  Dipartimento di Fisica, Unversità degli Studi di Pavia, via
Bassi 6, 27100 Pavia, Italy  Iacopo Carusotto  Pitaevskii BEC Center, CNR-INO
and Dipartimento di Fisica, Università di Trento, I-38123 Trento, Italy
Fabrizio Riminucci  Molecular Foundry, Lawrence Berkeley National Laboratory,
One Cyclotron Road, Berkeley, California, 94720, USA  Kirk W. Baldwin  Loren
N. Pfeiffer  PRISM, Princeton Institute for the Science and Technology of
Materials, Princeton University, Princeton, New Jersey 08540, USA  Giovanni I.
Martone  INFN, Sezione di Lecce, 73100 Lecce, Italy  CNR Nanotec, Institute of
Nanotechnology, via Monteroni, 73100, Lecce, Italy  Milena De Giorgi  Dario
Ballarini  Daniele Sanvitto  CNR Nanotec, Institute of Nanotechnology, via
Monteroni, 73100, Lecce, Italy

###### Abstract

A supersolid is a counter-intuitive phase of matter where its constituent
particles are arranged into a crystalline structure, yet they are free to flow
without friction. This requires the particles to share a global macroscopic
phase while being able to reduce their total energy by spontaneous, spatial
self-organisation. This exotic state of matter has been achieved in different
systems using Bose-Einstein condensates coupled to cavities, possessing spin-
orbit coupling, or dipolar interactions. Here we provide experimental evidence
of a new implementation of the supersolid phase in a novel non-equilibrium
context based on exciton-polaritons condensed in a topologically non-trivial,
bound-in-the-continuum state with exceptionally low losses. We measure the
density modulation of the polaritonic state indicating the breaking of
translational symmetry with a remarkable precision of a few parts in a
thousand. Direct access to the phase of the wavefunction allows us to
additionally measure the local coherence of the superfluid component. We
demonstrate the potential of our synthetic photonic material to host phonon
dynamics and a multimode excitation spectrum.

The existence of a supersolid phase of matter, combining a crystalline
structure with superfluid properties, was speculated more than 50 years ago
PhysRev.106.161, 1969JETP...29.1107A, PhysRevLett.25.1543,
1971JETP...32.1191K, THOULESS1969403, GROSS195857 however only recently there
have been convincing experimental evidence, mainly using ultracold atomic
Bose-Einstein condensates (BECs) coupled to electromagnetic fields. There,
various guises of the supersolid were created using atoms coupled to high-
finesse cavities leonard.nature.2017, gopalakrishnan.nphys.2009,
guo.nature.2021, with large magnetic dipole moments tanzi.nature.2019,
tanzi.science.2021, PhysRevLett.122.130405, PhysRevX.9.021012,
bottcher.prx.2019, klaus.nphys.2022, guo.nature.2019, bottcher.prx.2019, and
spin-orbit-coupled, two-component systems showing stripe phases
li.nature.2017, putra.PhysRevLett.124.053605.

At the mean field level, supersolidity can be interpreted as two-mode
condensation; after condensation in the first, a second mode becomes
energetically available by tuning of interactions or an external
electromagnetic field and can then be dynamically populated. In a single-mode
BEC, the U(1) symmetry related to its phase is broken at the condensation
phase transition, giving rise to superfluidity in the system. The mean field
wavefunction has $k=0$ momentum, giving a constant density in space.
Supersolidity arises, in most cases, as an additional, macroscopic occupation
of one (or more) finite-momentum modes of the system leading to breaking of
translational symmetry. The occupation is preceded by the presence of a roton-
like minimum in the excitation spectrum that can be dynamically populated when
its energy is tuned to be equal to that of the stationary ground state
chomaz.nphys.2018, PhysRevA.90.063624, PhysRevLett.114.105301. The new mean
field wavefunction acquires therefore density modulations with a given length
scale and a macroscopic phase relation, giving rise to an extra U(1) symmetry
breaking. A change in the relative phase of the two condensates corresponds to
a rigid translation of the density pattern associated to a Goldstone mode of
the system wouters.goldstone.PhysRevA.76.043807. A gapped Higgs mode is also
present, corresponding to a change in the amplitude of the density modulation
leonard.science.aan2608.

However, despite all the progress in understanding the physics of the
supersolid phase, most current realisations are limited to ultracold atomic
systems. Here we demonstrate a new photonic platform that can give rise to a
supersolid phase in a novel driven-dissipative context and is poised to open
new directions for the study of these intriguing phases of matter. We use a
BEC of polaritons, i.e., strongly coupled light-matter excitations, formed in
a sub-wavelength, patterned waveguide that gives rise to a condensate at a
saddle point of the dispersion ardizzonePolaritonBoseEinstein2022,
kasprzak_boseeinstein_2006, carusotto_quantum_2013, amo_superfluidity_2009,
sanvitto_all-optical_2011. Under some mode-folding conditions and exciton-
photon detuning, the photonic crystal waveguide shows two extra modes
symmetrically distributed in $k$ with respect to the BEC at $k=0$, with the
same polarisation. This allows for the formation of energy-degenerate, optical
parametric oscillation (OPO) between the condensate and the adjacent modes,
which can be described in terms of polariton-polariton interactions arising
from a $\chi^{(3)}$ nonlinearity. Parametric processes are commonplace in
nonlinear photonics and can exhibit pattern formation close to threshold in
multimode systems pattern.formation.PhysRevLett.83.5278,
ardizzone2013.scientific.reports. However, differently from our results, they
are usually induced resonantly which imposes the coherence of the pump laser
on the wavefunction whittaker_polariton_2017, romanelli_four_2007,
lewandowski_polarization_2016. Recently, the onset of a non-resonantly-driven,
energy-degenerate parametric scattering was demonstrated in different systems
wu_nonlinear_2021, sawicki_polariton_2021, PhysRevB.87.155302,
doi:10.1021/acs.nanolett.1c04800, PhysRevB.107.L060303, PhysRevLett.108.166401
showing how these scattering processes are ubiquitous; however, no proof of
coherence nor pattern formation was reported. Here, instead, we used a
$\Gamma$-point, symmetry protected, bound-in-the-continuum (BiC) state formed
at the lower branch of the anticrossing ardizzonePolaritonBoseEinstein2022,
gianfrateReconfigurableQuantumFluid2024, which led to a few crucial
advantages: (1) a vanishing linewidth due to the supression of losses on the
photonic component, which reduces both the condensation threshold and the
incoherently scattered background, and does not overwhelm the emission from
other modes in the system, (2) a negative effective mass, which makes the
$k=0$ component self-localised, and (3) absence of accessible lower-energy
extrema in the dispersion where particles can accumulate; as such, the OPO is
adequately described by a three-branch iso-energetic model. These striking
properties of the BiC render the small density modulation of the BEC clearly
visible.

Figure 1: System in reciprocal and coordinate space. a. The mechanism that
leads to the formation of the supersolid. Linear and nonlinear scattering
processes combine and give rise to a density modulation in coordinate space.
b. The dispersion comprises of four fundamental propagating modes. The
measured dispersion closely matches the numerical dispersion obtained from
diagonalising an effective $8\times 8$ Hamiltonian (see SI). c. The BiC state
$\ket{0}$ at $k=0$ parametrically scatters to the modes $\pm k_{r}$, denoted
as $\ket{\pm 1}$. At powers above threshold, the BEC needs to be masked in
order for the $\ket{\pm 1}$ population to be measured. The parametric
scattering generates photon pairs that populate the $\ket{\pm 1}$ states as
indicated by the arrows. d. The density of the supersolid wavefunction
$\psi(x,y)$ shows a characteristic fast modulation due to the spontaneous
breaking of translational symmetry. The amplitude of the modulation along the
$y=0$ line is $\sim$2.6%.

##  I Formation of the density modulation

We prepared a polariton BEC, using off-resonant, pulsed excitation, at the
lower branch of an anticrossing formed by folding the propagating modes of the
one-dimensional waveguide to $k=0$ using a sub-wavelength, etched grating. The
corresponding wavevector sets the size of the irreducible Brillouin zone,
which is much larger than all other characteristic length scales of the
system, and lies outside our numerical aperture (see Methods). Aside from the
BiC, which corresponds to the anti-crossing of the left- and right-propagating
TE0 modes, the waveguide supports also TE1 modes with the mode spacing
determined by the planar waveguide thickness and index contrast, which are
responsible for the formation of the modulation. The off-resonant pumping
predominantly populates the BiC due to its negative effective mass
riminucciPolaritonCondensationGapConfined2023, nigro.PhysRevB.108.085305. It
is important to note that the TE1 modes are populated by some incoherent
population relaxing from the exciton reservoir, however they can also become
macroscopically populated via linear and in particular nonlinear coupling with
the TE0 (Fig. 1a).

The mode-spacing sets the natural unit of length for our system to be
$1/k_{r}=\qty{2.3}{\micro m}$, so we can write the four basis states as
$\ket{\pm,mk_{r}}$ for $m=0,1$, where the sign labels the propagation
direction of the relevant photonic dispersion branch. The TE0 modes are
coupled with a Rabi frequency of 1.5meV, and the BiC corresponds to the lowest
energy adiabatic state $\ket{0}=(\ket{+,0}+\ket{-,0})/\sqrt{2}$ that has zero
group velocity, with an energy $E=\qty{1.526}{eV}$. We will from now on limit
the description to the three states $\ket{0},\ket{\pm 1}\equiv\ket{\pm,k_{r}}$
to simplify notation. The single-particle dispersion shown in Fig. 1b is
measured below threshold since the formation of the BEC dominates the
luminescence above threshold and renders the rest of the dispersion scarcely
visible. The BEC is visible in the dispersion at the $\Gamma$-point, slightly
blueshifted from the bare dispersion due to intrinsic polariton-polariton
interactions, Fig. 1c ardizzonePolaritonBoseEinstein2022. The excitonic
Hopfield coefficient at the BiC energy is estimated to be 0.54 for $\ket{0}$
and 0.4 for $\ket{\pm 1}$, due to their different mode structure that leads to
them coupling to the exciton with different Rabi frequencies.

We will assume, in what follows, no linear coupling between TE0 and TE1 modes.
i.e. the off-diagonal element in the hamiltonian $U_{01}\to 0$ (see SI). When
the population of $\ket{0}$ becomes large enough, parametric scattering driven
by a four-polariton interaction term of the form
$a_{+1}^{\dagger}a_{-1}^{\dagger}a_{0}a_{0}+h.c.$, where $a_{m}$ are bosonic
operators, results in emission of polariton pairs that populate the $\ket{\pm
1}$ states at finite momentum according to energy and momentum conservation
criteria. The $\ket{0}$ and $\ket{\pm 1}$ states are commonly referred to as
pump, signal, and idler modes in optical parametric processes. Since the
interaction term is invariant to the transformation $a_{m}\to
a_{m}e^{mi\phi}$, with $m=0,\,\pm 1$ wouters.PhysRevA.76.043807, it has U(1)
symmetry which is spontaneously broken at OPO threshold. At the same time the
OPO breaks translational symmetry, due to it being energy-degenerate, which
leads to the appearance of a density modulation on top of the condensate
carusotto.PhysRevB.72.125335.

Thanks to the small radiative coupling of the BiC with the far field, the
luminescence of the BEC is strongly suppressed at $k=0$, showing the
characteristic two-lobe structure with a $\pi$ phase jump between them (Fig.
1d). This is essentially due to the co- and counter-propagating TE0 modes
being out of phase at $k=0$ as they Bragg scatter from the grating; the
propagation into the far field leads to the characteristic nodal line at the
centre of the BEC. The density modulation, indicating the translational
symmetry breaking of the supersolid, is visible atop both lobes of the BEC,
which are extended in reciprocal space and outside the laser spot in
coordinate space. Its amplitude is shown in Fig. 1e to be $A\sim$2.6%. The
modulation pattern extends over $2k_{r}L\simeq 30$ lattices sites, where
$L=\qty{30}{\mu m}$ is the effective $1/e^{2}$ waist of the BEC. The presence
of the modulation on top of the BEC is a clear sign of crystallisation: the
observation of a periodic pattern, similar to the formation of crystals,
despite the system otherwise still being in the superfluid regime.

The density in Fig. 1d is the integrated emission over approximately
$8\text{\times}{10}^{6}$ pulses exciting the sample at a rate of 80MHz, much
slower than the lifetime of the BEC (of the order of a few hundred
picoseconds); as such, each measurement is an average of multiple independent
BEC implementations. The dynamical behaviour of the optical parametric process
near- and above-threshold can be simply described by a multibranch Hamiltonian
involving three modes PhysRevB.87.155302, langbein.PhysRevB.70.205301,
ciuti.PhysRevB.63.041303. The single-mode solution becomes unstable with
increasing power and a three-mode solution becomes the ground state of the
system dunnett.PhysRevB.98.165307, PhysRevB.87.155302.

Figure 2: Parametric scattering. a. Growth of the BEC at $k=0$. Above the
condensation threshold power 250W both the blueshift and linewidth start
increasing to about 1meV, typical of pulsed excitation. The first-excited mode
appears for powers above 800W at the energy of the lower band. b. Emission
from $k=0$ and $k=k_{r}$ that are proportional to the poalriton BEC density.
Coherent and incoherent scattering processes can be separated since they have
different growth behaviour with the BEC density. The $n_{\pm 1}$ density
experiences sublinear growth above threshold that then shows a second
threshold in efficiency at about $\times 1.72$ above threshold power. The
insets show how the distribution of the density on the TE1 modes changes from
uniform at low powers to more localised along $k_{x}$ above the second
threshold.

Figure 2a shows that condensation takes place at around 0.3mW of pump power.
The luminescence and linewidth of the BEC dramatically increase above
threshold due to the pulsed excitation. Since our system is trapped, it
supports multiple confined modes; the first excited mode enters in the gap
around 0.8mW, and can be clearly seen in Fig. 2a at lower energies than the
main, blueshifted BEC riminucciPolaritonCondensationGapConfined2023. Figure 2b
shows the growth of the densities $n_{\pm 1}$ versus $n_{0}$, showing
approximate linear growth below BEC threshold. This suggests that at low
excitation power, the emission is dominated by incoherent photoluminescence
coming from the exciton reservoir. The $\ket{0}$ mode is populated much more
efficiently and $n_{\pm 1}$ starts to experience sublinear growth above
threshold. At roughly 430μW, the efficiency of the scattering process
increases indicating that nonlinear scattering is taking place. In this region
the $\ket{\pm 1}$ states, that are uniformly populated in reciprocal space for
low pump powers, become structured showing conjugate scattered pairs (see
insets of Fig. 2b) and the density modulation becomes even more evident in
coordinate space.

Our ability to observe the density modulation comes from the small linear
coupling present in the system. Since the $\ket{\pm 1}$ states are not empty
when the OPO process establishes, but are populated by the exciton reservoir
or via Rayleigh scattering, even a small amount of linear coupling contributes
towards fixing the phase of the density modulation and as such renders it
visible in time-integrated emission. The linear coupling is equivalent to
including two-operator exchange terms of the form
$a_{+1}^{\dagger}a_{0}+a_{-1}^{\dagger}a_{0}+h.c.$ in the hamiltonian, which
make it not invariant under the phase transformation discussed above; formally
the OPO tends to break the U(1) symmetry in a non fully spontaneous manner
because of this. Well above the OPO threshold the nonlinear scattering
dominates and a degree of phase freedom is recovered. To show this we consider
the ratio $K=(n_{1}+n_{-1})/2n_{0}$, which is effectively the Bragg signal
corresponding to a given modulation amplitude, and scales as $A^{2}/4$, where
$1/4$ is the structure factor for a sinusoidal modulation. For our observed
density modulation, this gives an expected Bragg signal of 0.02%, more than an
order of magnitude smaller than the measured 0.5%. Moreover, we note that (1)
the phase of the modulation varies in the transverse direction by about a
$1.1\pi$ rotation, as seen in Fig. 1d, and (2) the average phase in
continuous-wave excitation fluctuates predominantly at the density maxima (see
discussion in SI). These observations indicate underlying phase fluctuations
related to the OPO that are averaged out during an integrated measurement and
effectively reduce the visibility of the density modulation. Analogous effects
exist in atomic supersolids, where the phase of the modulation is fixed by the
presence of an external trapping potential, yet small phase variations could
still be observed, e.g. in dipolar supersolids PhysRevLett.122.130405. We note
the linear coupling term alone does not fix the phase of the modulation as
discussed for spin-orbit-coupled systems li.yun.PhysRevLett.108.225301,
geier.PhysRevLett.130.156001. Such effects are inherent in finite-size
systems, however they have not prevented the observation of well-behaved
collective excitation spectra tanzi.nature.2019, tanzi.science.2021.
Importantly, our system can recover full phase freedom in a much simpler
manner, i.e. by engineering the dispersion so that $U_{01}\to 0$.

##  II Local and global coherence

Having shown the behaviour of the populations undergoing parametric
scattering, it is interesting to see how the diagonal and off-diagonal phase
coherence is established. We use interferometric measurements to fully capture
the microscopic, spatial coherence of the system, which serves as a direct
proof that the fragile local and long-range order of the solid and superfluid
parts is preserved. The strength of using an all-optical system in this regard
is in the ease of access to the first-order spatial correlator,
$g^{(1)}(x-x^{\prime})$, due to our ability to measure the full wavefunction
of the BEC $\psi(x,y)=\sqrt{n(x,y)}\exp(-i\phi(x,y))$ using off-axis
holography (see Methods). Direct, local measurements of the first-order
coherence show the wavefunction to be fully coherent above threshold with a
modulated coherence amplitude locally correlated with the emerging supersolid
crystalline structure. The superfluid character of our bulk BEC on the BiC
state has already been characterised through the observation of a linear
Bogoliubov excitation spectrum nina_collective_2023, and as such, here we take
the observed increase in the $g^{(1)}(x-x^{\prime})$ coherence length (Fig.
3a) to be indicative of an increased superfluid fraction.

Figure 3a shows the increase with power of $g^{(1)}(x-x^{\prime})$ normalised
to its maximum. Below threshold only the the zero-time-space auto-correlator
$g^{(1)}(0)$ has appreciable amplitude. At ∼250μW, the onset of polariton
condensation, the coherence length starts increasing at a rate of 133μm/mW.
The fast spatial oscillations that appear in $g^{(1)}(x-x^{\prime})$ have an
amplitude of $\sim$1.5% and are in phase with the density modulations, showing
how coherence locally drops in regions of lower density of the crystalline
solid PhysRevLett.25.1543 (see Fig. 3a).

Figure 3: Spatial coherence through threshold. a. The first-order correlator
$g^{(1)}(x-x^{\prime})$ is a direct measurement of coherence of the superfluid
part. Since it is normalised to unity at each power, the increase of the width
of the distribution is proportional to the increase of the coherence length
(see right-hand panels) marked by the white diagonal lines with a slope
133μm/mW. The modulations on the $g^{(1)}(x-x^{\prime})$ are correlated with
the density modulations showing a local drop of coherence at the individual
crystal lattice sites (lower right panel at 754μW). b. Growth of $C_{mod}$
versus $C_{coh}$ over threshold. The trend is similar to Fig. 2b showing the
existence of two thresholds.

The emergence of the global, diagonal and off-diagonal long-range order can be
quantified by measuring
$\mathcal{F}_{x}\\{|\psi(x,y)e^{-ik_{1}x}+\psi(-x,-y)e^{-ik_{2}x}|^{2}\\}$ at
$y=0$, where $\mathcal{F}_{x}\\{\cdot\\}$ is the spatial Fourier transform
along $x$ and the quantity in the brackets is the output of the interferometer
(see SI). The magnitude of the $k=k_{r}$ component is a measure of the
amplitude of the density modulation arising from the phase coherence between
the three BEC components. Similarly, the component corresponding to the
wavevector of the interferogram fringes, $k=k_{2}-k_{1}$, captures the phase
coherence of the $k=0$ component. We denote these quantities as $C_{mod}$ and
$C_{coh}$ respectively. The relative increase of $C_{mod}$ against $C_{coh}$,
above threshold pump power, (see Fig. 3b) is similar to the behaviour of
$n_{\pm 1}$ against $n_{0}$ shown in Fig. 2b. Below the condensation
threshold, $C_{mod}$ is largely dominated by noise since the crystalline
pattern has not formed yet, then it slowly rises up to the OPO threshold where
it begins to increase sharply. In absence of linear coupling this would
indicate a two-step breaking of the superfluid U(1) and translational
symmetries; in our case the density modulation increasing is a clear signature
that the OPO process is established.

##  III Non-rigidity of the modulation

An important aspect of the study of quantum many-body systems is understanding
their collective excitations; they dictate the dynamics of the system and its
response to external perturbations. In order to be able to observe and
manipulate the collective modes, the quantum system needs to be allowed to
spontaneously break the relevant symmetry during a phase transition, without
external or intrinsic factors that influence this choice. For example, in
single-mode, high-finesse cavity realisations of the supersolid phase the
wavelength of the modulation is fixed by the wavelength of the relevant cavity
mode, implying that the phononic branch corresponding to the Goldstone mode
cannot be observed leonard.science.aan2608. In such systems, excitations at
finite $k$ are suppressed due to the infinite range of the coupling. This
means that the global translation mode at $k=0$ exists, but the local
compression/dilation motion of the fringes occurring at finite $k$ is not
observable; such phonon dynamics can be captured with multimode cavities or in
stripe phases of spin-orbit-coupled systems, as well as dipolar systems due to
the translational invariance of the interactions in those systems
guo.nature.2021, li.nature.2017, geier.PhysRevLett.130.156001. In our case,
the symmetry breaking results from short range polariton-polariton
interactions, we therefore expect similar behaviour to the latter systems.

The parametric scattering process is constrained to occur on the available
$\ket{\pm 1}$ modes where the density of states is different from zero. This
fixes the wavenumber of $\ket{\pm 1}$ to $k_{r}$. However the value of $k_{r}$
depends on polariton-polariton interactions and their resulting blueshift. As
the interactions blueshift the energy of the BEC, and since OPO is a iso-
energetic process, the energy of $\ket{\pm 1}$ has to shift as well, and this
is done by shifting their wavevector inwards. This shift is nonlinear due to
the strong-coupling of the photonic modes to the exciton, which leads to an
increased curvature at energies closer to the exciton resonance. The energy of
$\ket{\pm 1}$, associated with the breaking of translational symmetry,
increases by 0.35 meV as the pump power increases (Fig. 4a); the linewidth of
the TE1 mode is about 500μeV. For the same range of pump powers $k_{r}$ shifts
by 0.04m^-1 (Fig. 4b, c). As such, the wavevector of the density modulation is
not fixed, but depends on both the single-particle and manybody parameters of
the system, reminiscent of a non-rigid, roton-like behaviour.

Figure 4: Non-rigid scattering. a. The energy and momentum of $\ket{-1}$
changes with excitation power. For high powers the first-excited mode is
visible (right panel). The position of $\ket{-1}$ in the waveguide dispersion
defines the interference part of the supersolid. Since the modes are curved
due to them being strongly coupled to the exciton, the position moves both in
energy b. (left) and in wavenumber (right) giving rise to a non-rigid lattice
structure. The shift of the central peaks follows the dispersion curve shown
in the top panels. The secondary peak appearing at lower energy is the first
excited state that enters at $k=0$ for pump powers $>\qty{0.8}{mW}$ (marked by
the black solid line) and scatters to the TE1 modes. c. Energy dispersion
averaged over all excitation powers shwoing the movement of $\ket{-1}$ in
reciprocal space.

##  IV Discussion and outlook

We demonstrated evidence of a supersolid state of matter in a polaritonic
system which is a novel and flexible platform for the investigation of the
physics of supersolid systems. To place our work more generally in the context
of supersolids, as observed mainly in ultracold atomic systems, we
experimentally addressed its characteristics in coordinate and reciprocal
space, its diagonal and off-diagonal long range order, and its interaction-
dependent behaviour. Our supersolid state is formed using an iso-energetic OPO
process to spontaneously break the translation symmetry of a non-resonant
polaritonic superfluid. We stress that this is a novel mechanismg for the
creation of a supersolid, particular of the driven-dissipatve context of non-
equilibrium polariton systems, and not simply a photonic analogue of
mechanisms demonstrated in atomic platforms. The analysis of the
parametrically scattered light indicates that the system condenses at a given
threshold and OPO activates at higher powers. Theoretically, we also expect
the existence of two dictinct transition points.

By exploiting the richness of engineering opportunities in these photonic
crystal systems it is possible to eliminate completely the linear coupling
between the modes, for example using TE- and TM-polarised modes. This will
enable us to observe full phase freedom of the density modulation either in
single-shot experiments or higher-order coherence functions
carusotto.PhysRevB.72.125335, folling.spatial.coherence.nature.2005. Near and
above threshold the statistics of the emitted light changes to bunched leading
to an increase in density-density correlations that are peaked at threshold,
similar to a resonantly pumped OPO ardizzone_bunching_2012. A spatially
dependent $g^{(2)}(x-x^{\prime})$ will show signatures of density modulation
since it is only sensitive to how polaritons at relative distances are
correlated.

Since the spacing of the crystalline modulation is not rigid but influenced by
short range polariton interactions, the low-lying excitation spectrum of the
supersolid should reveal the two phononic Goldstone branches resulting from
the symmetries involved. The observation of such spectrum is a natural next
step in the investigation of the supersolid properties. As a further
development of this research, increasing the dimensionality of the system, we
expect additional breaking of discrete rotational symmetries to occur, as well
as the proliferation and dynamics of stochastic vortices
pomeau.PhysRevLett.72.2426. Moving beyond parametric scattering, similar
effects have been suggested utilising coupling to indirect excitons to
engineer long-range, dipolar interactions that can lead to crystallisation
under compression boning.PhysRevB.84.075130, or using Bose-Fermi mixtures
where excitons are coupled to an electron gas PhysRevLett.108.060401.

##  V Methods

Sample. The sample used in these experiments is an AlGaAs waweguide grown with
molecular beam epitaxy. It has twelve 20 nm embedded GaAs quantum wells with a
heavy hole excitonic transition at 1531.1(0.1). The propagating modes are
folded within the lightcone by inscribing a $50\,\times\,400\,\unit{\mu m}$
unidirectional grating within the structure with a lattice parameter
$a=\qty{242}{nm}$ corresponding to a wavevector translation of
$\pi/a=\qty{13}{\mu m^{-1}}$. The energy momentum dispersion therefore
consists of two counterpropagating TE0 modes forming an anticrossing at
$k_{x}=0$, 2.4meV blue detuned with respect to the excitonic transition and a
couple of symmetric TE1 modes located at $\approx\qty{2.3}{\mu m^{-1}}$ from
the TE0 modes. The filling factor of such grating is 76% in order to spatially
overlap the two $\ket{0}$ and ensure the coupling described by the $U$ term in
the Supplemental information responsible of the opening of the BiC gap. Due to
the confinement, the strong coupling regime with a Rabi frequency of 5.6 and
4.2 meV is achieved between the excitonic transition and the propagating TE0
and TE1 modes respectively. The fabrication details of the grating etching can
be found in riminucci_nanostructured_2022.

Laser excitation scheme. The sample is held in an attoDRY closed loop helium
cryostat at a temperature ≈4 in order to avoid exciton ionisation. The
waveguide is non-resonantly excited in reflection configuration using an 80MHz
repetition rate, 100 tunable laser. The polarisation is chosen to be parallel
to the grating corrugation direction, and the energy is set at 1.61. The
excitation beam is focused into a 3.90(0.03) $1/e^{2}$ spot through a NA =
0.68 aspheric objective, at the centre of the etched structure.

Interferometric measurements. The photoluminescence collected by the objective
is sent into a Michelson interferometer through the detection line. In one of
the two arms the image is centro-symmetrically flipped along $x$ and $y$ using
a backreflector and the coordinate-space photoluminescence distribution re-
focused onto the spectrometer slits through the same optical path. The total
magnification of the line, sample-to-camera, is approximately 72.6. The two
arms are synchronised and balanced, optimising the visibility of the
interference pattern of the excitation laser reflected from the sample
surface.

Reciprocal space measurements For the reciprocal space images, a similar
detection line is used, but the interferometer is bypassed. A 4f-system
collects the reciprocal plane and focuses it onto the monochromator slits,
resulting in a total magnification of the objective back focal plane equal to
1.66.

Coherence calculation. The two-point correlator
$g^{(1)}(x-x^{\prime})=\braket{E^{*}_{1}(x)E_{2}(x)}/\braket{}{E_{1}(x)}{{}^{2}}\braket{}{E_{2}(x)}{{}^{2}}^{1/2}$,
where $E_{1,\,2}(x)$ are the fields from the two interferometer arms, is
proportional to the visibility $\mathcal{V}(x)=(I_{u}-I_{l})/(I_{u}+I_{l})$ of
the interferometric fringes, where $I_{u}$ ($I_{l}$) is the upper (lower)
envelope of the interferogram; we omit the $x$ dependence for brevity. We
calculate $\mathcal{V}$ by first extracting the $k=k_{r}$ and $k=0$ components
of the interferogram by appropriate filtering in the reciprocal space and then
using the Hilbert transform to calculate the envelope of the ac component; the
visibility is then simply the envelope over the $k=0$ component. For slight
imbalance or not full reflection symmetry between the two images,
$\mathcal{V}$ needs to be rescaled by $2\sqrt{I_{1}I_{2}}/(I_{1}+I_{2})$,
where $I_{1,2}$ are the intensities in each arm of the interferometer, to give
a correctly normalised $g^{(1)}(x-x^{\prime})$. We forgo this normalisation
since small misalignments of the order of $1/k_{r}$ in the overlap of the two
images distorts the density modulation. This does not significantly alter the
form of $g^{(1)}(x-x^{\prime})$ since in our case this prefactor is $\simeq
1$.

Acknowledgements. We acknowledge fruitful discussions with V. Ardizzone, M.
Pieczarka, and A. Recati. We are thankful to G. Lerario for his valuable
feedback on the manuscript. This project was funded by the Italian Ministry of
University (MUR) PRIN project “Interacting Photons in Polariton Circuits” –
INPhoPOL (grant 2017P9FJBS); the project “Hardware implementation of a
polariton neural network for neuromorphic computing” – Joint Bilateral
Agreement CNR-RFBR (Russian Foundation for Basic Research) – Triennal Program
2021–2023; the MAECI project “Novel photonic platform for neuromorphic
computing”, Joint Bilateral Project Italia - Polonia 2022-2023; the PNRR MUR
project: “National Quantum Science and Technology Institute” - NQSTI
(PE0000023) co-funded by the European Union - NextGeneration EU; the Apulia
Region, project “Progetto Tecnopolo per la Medicina di precisione”, Tecnomed 2
(grant number: Deliberazione della Giunta Regionale n. 2117 del 21/11/2018);
the PRIN project “QNoRM: A quantum neuromorphic recognition machine of quantum
states” - (grant 20229J8Z4P). This research has been cofunded by the European
Union - NextGeneration EU, “Integrated infrastructure initiative in Photonic
and Quantum Sciences” - I-PHOQS [IR0000016, ID D2B8D520, CUP B53C22001750006].
IC acknowledges funding from the Provincia Autonoma di Trento, partly via the
Q@TN initiative. This research is funded in part by the Gordon and Betty Moore
Foundation’s EPiQS Initiative, grant GBMF9615 to L.P., and by the National
Science Foundation MRSEC grant DMR 2011750 to Princeton University. Work at
the Molecular Foundry is supported by the Office of Science, Office of Basic
Energy Sciences, of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231. We thank Scott Dhuey for assistance with electron beam
lithography and Paolo Cazzato for the technical support.

Competing interests. The authors declare no competing interests.

Data Availability. The data of this study is available from the corresponding
author upon reasonable request.

Author contributions. DT and ML conceived the experiment and convinced DS to
go ahead with it. AG took the data and together with DT performed the
analysis. DT, AG and ML wrote the manuscript. DN, DG, and IC provided
theoretical support. FR processed the sample; growth was performed by KB and
LP. DS supervised the work, discussed the data, and tirelessly explained to DT
how polariton OPOs should behave. DT tried to use the best photonics language
he could muster but in the end had to use kets. All authors contributed to
discussions and editing of the manuscript.

## References

Supplementary Information

##  I Dispersion of the waveguide

The particle-number-conserving part of the system Hamiltonian can be written
as a $8\times 8$ matrix

$\definecolor{light-
gray}{gray}{0.8}\mathcal{H}(k)=\left(\begin{array}[]{cc|cc|cc|cc}\color[rgb]{0,0,1}\delta_{0}&\color[rgb]{0,0,1}U_{00}&V_{0}&0&0&0&0&U_{01}\\\
\color[rgb]{0,0,1}U_{00}&\color[rgb]{0,0,1}-\delta_{0}&0&V_{0}&0&0&U_{01}&0\\\
\hline\cr V_{0}&0&0&0&0&0&0&0\\\ 0&V_{0}&0&0&0&0&0&0\\\ \hline\cr
0&0&0&0&0&0&V_{1}&0\\\ 0&0&0&0&0&0&0&V_{1}\\\ \hline\cr
0&U_{01}&0&0&V_{1}&0&\color[rgb]{1,0,0}\delta_{1}&\color[rgb]{1,0,0}U_{11}\\\
U_{01}&0&0&0&0&V_{1}&\color[rgb]{1,0,0}U_{11}&\color[rgb]{1,0,0}\delta_{-1}\end{array}\right),$
(S1)

which is represented on the basis $\\{\ket{\pm 0},\ket{X}^{4},\ket{\pm 1}\\}$
comprising the four TE0 and TE1 co- and counter-propagating modes and the
exciton which couples independently to left- and right-propagating modes of
different momenta. Here we take the energy of the bare exciton resonance as a
reference, i.e., we set it to zero. The linear dispersions of the modes
$\delta_{j}\equiv\delta_{j}(k)=\delta\pm u(k-jk_{r})$, where $j=0,\pm 1$,
$\delta$, act as effective detunings that are identical up to a translation in
$k$; $\delta$ is the photonic detuning, $u$ the photonic group velocity, and
the (plus) minus sign is for (right) left-propagating modes. The photonic
modes are coupled with Rabi frequencies $U_{ij}$, that depend on the mode
overlap with the etched grating, and then couple to the degenerate exciton
mode with Rabi frequencies $V_{0,1}$, with $V_{1}<V_{0}$ due to the different
mode distributions of $\ket{0}$ and $\ket{\pm 1}$. The block structure of the
Hamiltonian makes it straightforward to identify the contributing parts. The
photonic part of the BiC can be described with a $2\times 2$ matrix and is
analogous to a spin-$1/2$ system (top left block). The counter and co-
propagating TE1 modes exist at the same energy as the TE0 modes and define the
period of the supersolid (bottom right block).

The BiC state can be simply described as a dissipative spin-$1/2$ system

$\mathcal{H}_{B}(k)=\begin{pmatrix}\delta_{0}&U_{00}\\\
U_{00}&-\delta_{0}\end{pmatrix}+i\gamma_{0}\begin{pmatrix}1&-1\\\
-1&1\end{pmatrix},$ (S2)

with eigenvalues
$\lambda_{1,2}=-i\gamma_{0}\pm\sqrt{\delta_{0}^{2}+(U_{00}-i\gamma_{0})^{2}}$.
The near-zero linewidth of the BiC stems from the fact that at
$\delta_{0}=ak=0$ one of the eigenvalues is purely real (and corresponds to
the BiC state) while the other one has twice the dissipation $2i\gamma$.

##  II Resonant excitation of propagating modes

Figure S1: Resonant excitation of propagating modes. a. Two resonant spots
injected at the energy of the BiC interfere and give rise to an optical
lattice structure. b. The spots enter resonantly at $\ket{\pm 1}$ and
propagate towards $k=0$. The interference pattern, shown in c. has twice the
spatial frequency of that of the supersolid (off-resonant) since $\ket{\pm 0}$
is absent. The supersolid peak is split due to the non-connected spatial
structure of the BiC.

By resonantly exciting the TE±1 modes at the energy of the BiC we effectively
isolate the “solid” part of our polaritonic supersolid and measure its
characteristics. The interference that arises from the diffracted peaks has
$k=2k_{r}$. Note that this in contrast with results discussed in the main
text. Indeed, since in the present case the total emitted field is the
superposition of two plane-wave contributions counter propagating, without any
signal coming from TE0, the field intensity is $\cos(k_{r}\cdot
x)^{2}\simeq\cos(2k_{r}\cdot x)$, resulting in a $k=2k_{r}$ peak in the
emitted spectrum. To avoid interference coming directly from the laser we use
two distinct excitation spots for $\pm k_{r}$ at a distance of approximately
10μm and measure the interference of the propagating modes. Figure S1a shows
the polariton density in coordinate space normalised to the density of the
individual laser beams; the interference in the region between the spots is
evident. Both lasers enter at $\ket{\pm 1}$ (see Fig. S1b) and due to the
negative group velocity propagate towards $k=0$. Notice that the fringes do
not have an apparent curvature which is characteristic of the phase
fluctuations due to the spontaneous symmetry breaking of the supersolid.
Moreover, the Fourier transform has a single well-defined peak since the
gaussian beams have a slowly-varying spatial structure. Interference between
the BiC and $\ket{\pm 1}$ on the other hand, shows multiple peaks around
$k_{r}$ corresponding to the two BiC lobes (see Fig. S1c).

##  III Linear coupling and finite size effects

Since we are imaging an ensemble average of many supersolid realisations it
might seem surprising that we observe the breaking of translational symmetry
in the first place. Should this phase transition be truly spontaneous, the
phase of the density modulation would average to a uniform distribution in
$[0,2\pi]$ and significantly reduce the visibility of the fringes. As
explained in the main text, the observed residual modulation is due to a
symmetry breaking term resulting from the linear coupling of the modes in our
system so we would not recover full phase freedom even at the thermodynamic
limit unless $U_{01}\to 0$. We note that in the thermodynamic limit of very
large pump power, the nonlinear coupling scales more favourably with $n_{0}$
than the linear coupling and our system converges to an ideal supersolid. For
vanishing linear coupling, the next-order, finite-size effects are the
trapping potential from the laser and the negative effective mass
riminucciPolaritonCondensationGapConfined2023, and the length of the
waveguide. The length of the waveguides in this work is 300μm which is
comparable to the propagation length of the BEC and as such it reaches the
edges. However, we did not observe any differences using 1mm-long waveguides
with tapered edges that should minimise reflections, which means that such
standing waves do not appear in our system.

###  III.1 Spatial phase shift

When scattering from a moving grating the scattered orders pick up a phase
factor $\propto e^{imk_{G}(x-x_{0})}$, where $m$ is the diffraction order,
$k_{G}$ the grating wavevector, and $x_{0}$ the position of the grating with
respect to an arbitrary origin wise.PhysRevLett.95.013901. In our case the
origin is defined by the relative position of the laser spot and the grating.
Figure S2 shows the density modulation changing with the centre position of
the excitation spot $x_{0}$ along the waveguide axis. The density modulation
of our BEC can be considered as two plane waves impinging onto the grating. If
there is relative translation of the grating or the laser spot the modulation
acquires a phase shift proportional to $k_{r}x_{0}$. In our case, these small
spatial fluctuations are negligible since for a rigid translation of
$x_{0}=k_{r}^{-1}=\qty{2.4}{\mu m}$, much smaller than the extent of the BEC,
the phase winds over 2$\pi$.

Figure S2: Phase rolling. The phase of the density modulation changes smoothly
when varying the position $x_{0}$ of waveguide along $x$ keeping the pump
laser fixed. The phase winds over $w\pi$ for a rigid translation of
$k_{r}^{-1}$ indicating the pointing stability of the system.

###  III.2 Transverse phase fluctuations

A degree of phase freedom can be inferred by the transverse phase variations
of the density modulation. The two-dimensional density distribution is
proportional to $1+|a(y)|^{2}\cos{\big{(}k_{r}x+\theta(y)\big{)}}$, where
$a(y)$ and $\theta(y)$ are the transverse amplitude and phase of the density
modulation. Figure S3a shows the reciprocal amplitude and phase of a typical
BEC realisation. When $|y|<\qty{2.5}{\mu m}$ the peak maxima that define the
periodicity of the crystal are at $k_{r}$ as expected. Fluctuations increase
for larger values of $y$ as the density drops. Compare this with Fig. S1 where
$\ket{\pm 1}$ are excited resonantly and the interference pattern does not
show any transverse variation.

Figure S3: Transverse phase variations. a. Amplitude and phase of the spatial
frequencies of the BEC density. The periodicity of the crystal remains largely
constant when the intensity of the BiC is large enough, however the maxima of
the dislocations fluctuate for $|y|>\qty{2.5}{\mu m}$. b. The phase
distribution, at a fixed periodicity $k_{r}$, narrows with increasing power
spanning the full unit circle below threshold to about 135 at 500μW.

The transverse phase of the crystal structure $\theta(y)$ measured at $k_{r}$
is uniformly distributed in Fig. S3b for low powers (left panel) when the
system is still largely incoherent. The distribution narrows with increasing
pump power till it is mostly distributed over the first quadrant at 500μW
(right panel). Critical fluctuations at lower densities, close to threshold,
likely scramble the U(1) phase which is then pinned at higher densities.

###  III.3 Average phase fluctuations

Figure S4: Average phase fluctuations. The amplitude modulation for 50
different continuous-wave realisations of the supersolid (gray). As the pump
power goes above threshold the fluctuations are mainly visible in the minima
and maxima of the density. The standard deviation of the phase (right panel),
measured only at the maxima, reduces above threshold and seems to increase at
a power of about $\times 1.5$ the threshold power.

Here, we switched to continuous pumping in order to avoid averaging from the
pulsed laser as described in the main text. Figure S4 shows 50 BEC
realisations for powers around the condensation threshold which in this case
is around 14mW. Below this threshold the phase of the modulation is completely
random since the system is still incoherent, and its standard deviation is
about 0.7. The standard deviation is measured only in the region of 5μm to
10μm to avoid artefacts of fluctuations coming from near-zero density. The
standard deviation reduces monotonically till about $\times 1.5$ the threshold
power and saturates or slightly decreases again, indicating that the phase is
yet to be fully fixed in the high density regions. This measurement is not
sensitive to pointing fluctuations of the laser beam since it shows local
phase fluctuations only.

##  IV Holographic imaging

Figure S5: Holographic imaging. a. Fourier transform of the interferogram in
c. showing the main peak around 7.5μm^-1 and sidebands at $\pm k_{r}$ from it.
b. In energy resolved measurements the BiC is dark in the middle.c. An
interferogram as given by the Michelson interferometer described in Methods.
The bright spot in the middle is coming directly from the exciton reservoir
and its size is comparable to the laser spot. d. Retaining only the middle
frequencies leads to a density without the middle spot (bottom panel) since it
is incoherent, as opposed to retaining the low frequencies (top panel). e. In
both cases, the amplitude of the modulation and the density in the lobes of
the BiC is exactly the same.

The small excitation spot creates a hot exciton reservoir underneath it that
also emits photolumiscence covering the BiC notch at $k=0$ (see Fig. S5c and
d). This feature does not exist in energy-resolved measurements, as shown in
Fig. S5a. By taking advantage of interferometric techniques we can remove it
also in energy-averaged measurements and recreate a realistic representation
of the density in coordinate space. The reciprocal-space representation of the
interferogram (Fig. S5a) shows, apart from the low-lying wavevectors, a strong
peak at $k_{2}-k_{1}=\qty{7.5}{\mu m^{-1}}$ which corresponds to the fringe
spacing from the interferometer. The visible sidebands on the fringe spacing
are at $\pm k_{r}$ and arise from the beating of the density modulation with
the interferometer fringes. The fringe spacing has been chosen to be $\approx
3k_{r}$ so that the lowest sideband at $k_{2}-k_{1}-k_{r}$ is above all the
low-lying wavevectors. Using a bandpass filter before transforming back to
coordinate space allows us to selectively eliminate the incoherent central
spot while otherwise leaving the density unaffected (see Fig. S5d and e).

## References