# Comment on paper: Position: Rethinking Post-Hoc Search-Based Neural
Approaches for Solving Large-Scale Traveling Salesman Problems

Yimeng Min  Department of Computer Science
Cornell University
Ithaca, NY 14850
Email<EMAIL_ADDRESS>
###### Abstract

We identify two major issues in the SoftDist paper: (1) the failure to run all
steps of different baselines on the same hardware environment, and (2) the use
of inconsistent time measurements when comparing to other baselines. These
issues lead to flawed conclusions. When all steps are executed in the same
hardware environment, the primary claim made in SoftDist is no longer
supported.

##  1 Introduction

In recent years, machine learning (ML) has emerged as a promising avenue for
addressing optimization problems like the Travelling Salesman Problem (TSP).
ML techniques, particularly those involving neural networks and reinforcement
learning, have shown potential in learning heuristics and patterns that can
guide the search for optimal routes more efficiently. By leveraging data, ML
models can improve the quality of solutions and reduce computation time. One
notable approach is the use of heat map-based search, where ML models generate
heat maps that highlight promising regions of the solution space. These heat
maps are then used to focus the search process, potentially enhancing the
efficiency and effectiveness of finding optimal or near-optimal solutions [1].

Recently, the authors of paper [2] (referred to as SoftDist) discussed the
neural approach and claimed:

> Our theoretical and experimental analysis raises doubts about the
> effectiveness of ML-based heat map generation. In support of this, we
> demonstrate that a simple baseline method can outperform complex ML
> approaches in heat map generation.

Here, however, we show that the authors in SoftDist misconducted the
experiments, leading to an unfair comparison and a flawed conclusion. Using
the correct comparison, the experimental results on large-scale instances like
TSP-1000 contradict their claims. This seriously undermines the main assertion
in their abstract.

In fact, their experiments demonstrate the power of unsupervised learning
neural approaches on large instances like TSP-1000. Even after running a grid
search on their proposed SoftDist model, their performance still falls short
compared to unsupervised learning neural approaches.

Figure 1: The overall scheme of heat map-based search. In the SoftDist paper,
when comparing the running time of different models, the authors calculate
parts (a) and (c) for SoftDist, while the UTSP running time consists of parts
(a), (b), and (c), leading to an unfair comparison.

##  2 Background of Heat map-Based Search on TSP

Heat map-based search is a novel approach applied to the TSP, where heat maps
are utilized to guide the search process. In this method, heat maps are
generated to represent the solution space, highlighting areas of interest.
These heat maps help guide the search, improving the efficiency and
effectiveness of finding optimal or near-optimal solutions. This technique
contrasts with traditional TSP methods by incorporating data-driven methods to
direct the search, offering a unique way to tackle this well-known
combinatorial optimization problem.

There are three steps in heat map-based search:

  * •

Generating the heat maps.

  * •

Performing I/O operations to save or load the generated heat maps.

  * •

Running the search.

The overall scheme of heat map-based search is shown in Figure 1. When
comparing different methods, two factors are usually considered: the time
required to build the heat maps and perform the search, and the quality of the
TSP solutions. A method that provides better quality solutions within the same
time cost is considered superior.

In the SoftDist paper, the authors compare the SoftDist method with other
approaches using a search method proposed in UTSP [1], as referenced in Table
3 of the SoftDist paper. Here, we show that there are two major problems in
the SoftDist paper when comparing their proposed SoftDist method to the
performance of UTSP [1], a data-driven neural approach that leverages the
power of unsupervised learning to generate heat maps for TSP without requiring
labels.

###  2.1 Problem 1

In the SoftDist paper, the author states:

> Our SoftDist heat map generation is performed on an NVIDIA A100 GPU. Due to
> its simplicity, the inference time is negligible (e.g., $<$ 0.1 seconds),
> and the GPU type has minimal impact. For fairness in comparison, all MCTS
> computations are conducted on an AMD EPYC 7V13 64-Core CPU @ 2.45GHz. We use
> 64 threads for TSP-500 and TSP-1000.

However, as shown in Figure 1, there are three steps in the heat map-based
search. The authors only highlight running the search on the same AMD EPYC
CPU, ignoring the other steps. Specifically, they rerun the search using their
environment but do not account for the heat map generation part. In other
words, the comparison is incomplete and biased, as the full computational
process for heat map generation on the NVIDIA A100 GPU is not considered. This
omission leads to misleading conclusions about the effectiveness and
efficiency of the proposed method compared to other approaches.

|  UTSP (the measure in UTSP paper) |  SoftDist (the measure in SoftDist) |  UTSP (the measure in SoftDist)
---|---|---|---
Generating heat maps | Included | Included | Included
I/O operation | Included | Not Included | Not Included
Search | Included | Included | Included
Total time | 1s + 1.45m + I/O | 0s + 1.44m (w/o IO) | 1s + 1.45m (w/o IO)
Performance | 1.65% | 1.73% | 1.65%
Table 1: Results on TSP-1000. The I/O operation always takes the same amount
of time because we are saving or loading heat maps of the same size.
Performance is measured by the gap with respect to the optimal solution; a
lower performance value indicates a better solution. The performance values
are directly taken from Table 3 in the SoftDist paper [2]. SoftDist
underperforms compared to UTSP within nearly the same time budget.

###  2.2 Problem 2

The time measurements used in SoftDist and UTSP are not the same, leading to
an unfair comparison. In SoftDist, the total time considered by the authors
includes only two parts: generating the heat map and running the search.
However, when reporting UTSP’s total running time, the time includes three
parts: generating the heat maps, I/O, and running the search. This discrepancy
results in an unfair comparison between the methods.

To ensure a fair comparison, we recalculated the time cost using the SoftDist
measurement. The results are shown in Table 1. When using the SoftDist
measure, which only considers generating the heat maps and running the search,
the UTSP method outperforms SoftDist with almost the same time budget. This
suggests that the heat map generated using UTSP can better guide the search.
Furthermore, since the generated heat maps should always have the same size,
the I/O time for different heat map-based methods, when running in the same
environment, should be the same.

Consequently, the main claim of the SoftDist paper,

> we demonstrate that a simple baseline method can outperform complex ML
> approaches in heat map generation,

no longer holds.

It should also be highlighted that the authors of SoftDist conducted a grid
search on their method while not performing any hyperparameter search on other
methods. For example, they used a hyperparameter of 0.0066 for TSP-500 and
0.0051 for TSP-1000. This lack of consistent hyperparameter tuning across all
methods further contributes to the unfair comparison and may skew the results
in favor of SoftDist.

Overall, the performance on TSP-1000, using the same time measurement, shows
that SoftDist [2], even with hyperparameter grid search, still underperforms
compared to data-driven neural approaches within almost the same time budget.
Consequently, the claim in the SoftDist paper that their method can outperform
ML approaches in heat map generation does not hold.

##  3 Conclusion

We identify two major issues in the SoftDist paper: (1) not running all steps
of different baselines in the same hardware environment, and (2) using
different time measurements when comparing to other baselines. These issues
lead to flawed conclusions. When using the same hardware environment, the main
claims of SoftDist are no longer supported.

These findings emphasize the importance of using comprehensive and fair
evaluation metrics when comparing different methodologies. It also highlights
the potential of unsupervised learning neural approaches in solving large-
scale instances of TSP, demonstrating their superiority in both efficiency and
effectiveness. Future research should focus on further refining these data-
driven methods and exploring their applications in other complex combinatorial
optimization problems.

## References

  * [1] Min, Yimeng, Yiwei Bai, and Carla P. Gomes, Unsupervised learning for solving the travelling salesman problem. Advances in Neural Information Processing Systems 36 (2024).
  * [2] Xia, Yifan, Xianliang Yang, Zichuan Liu, Zhihao Liu, Lei Song and Jiang Bian, Position: Rethinking Post-Hoc Search-Based Neural Approaches for Solving Large-Scale Traveling Salesman Problems. arXiv:2406.03503 Accepted by International Conference on Machine Learning (2024)


# Evidence for a Dense, Inhomogeneous Circumstellar Medium in the Type Ia SNR
0519-69.0

Brian J. Williams NASA GSFC, Greenbelt, MD 20771  Parviz Ghavamian Department
of Physics, Astronomy and Geosciences, Towson, MD, 21252;
<EMAIL_ADDRESS>Ivo R. Seitenzahl School of Science, University of New
South Wales, Australian Defence Force Academy, Canberra, ACT 2600, Australia
Stephen J. Reynolds  Kazimierz J. Borkowski  Department of Physics, North
Carolina State University, Raleigh, NC, 27695-0001  Robert Petre  NASA GSFC,
Greenbelt, MD 20771

###### Abstract

We perform an expansion study of the Balmer dominated outer shock of the SNR
0519$-$69.0 in the LMC by using a combination of new HST WFC3 imagery obtained
in 2020 and archival ACS images from 2010 and 2011. Thanks to the very long
time baseline, our proper motion measurements are of unprecedented accuracy.
We find a wide range of shock velocities, with the fastest shocks averaging
5280 km s-1 and the slowest grouping of shocks averaging just 1670 km s-1. We
compare the H$\alpha$ images from HST with X-ray images from Chandra and mid-
IR images from Spitzer, finding a clear anti-correlation between the
brightness of the remnant in a particular location and the velocity of the
blast wave at that location, supporting the idea that the bright knots of
X-ray and IR emission result from an interaction with a dense inhomogeneous
circumstellar medium. We find no evidence for X-ray emission, thermal or
nonthermal, associated with the fastest shocks, as expected if the fastest
velocities are the result of the blast wave encountering the lower density
ambient medium of the LMC. We derive an age of the remnant of $\leq 670\pm 70$
yr, consistent with results derived from previous investigations.

##  1 Introduction

Supernova remnants (SNRs) in nearby galaxies with well-known distances, such
as the Large Magellanic Cloud (LMC), offer an excellent opportunity for
studying the energetics of supernova explosions with great accuracy. Of the
young remnants of Type Ia supernovae (SNe Ia), the four Balmer-dominated
remnants DEM L71, SNR 0509-67.5, SNR 0519-69.0 and SNR 0548-70.4, have offered
some of the best testing grounds for studying such aspects as Type Ia
progenitor models (Seitenzahl et al., 2019), collisionless shock physics
(Smith et al., 1991; Ghavamian et al., 2007; van Adelsberg et al., 2008; Heng,
2010) and the physics of cosmic ray acceleration (e.g. Morlino et al., 2013).
All four Balmer-dominated SNRs were discovered as X-ray sources by the
Einstein Observatory (Long et al., 1981) and confirmed as SNRs by Tuohy et al.
(1982), where their Balmer-dominated nature was also established. The optical
emission from Balmer-dominated shocks is produced by collisional excitation of
neutral hydrogen atoms as they cross the shock. There are two populations of
neutral atoms, a cold population excited by collisions with post-shock
electrons and protons, and a hot population generated by charge exchange
between the cold hydrogen atoms and hot protons. The two components produce
broad and narrow H$\alpha$ line emission which can be used to diagnose the
shock speed and degree of electron-ion equilibration (e.g. Chevalier et al.,
1980; Smith et al., 1991; Ghavamian et al., 2001, 2007; van Adelsberg et al.,
2008; Morlino et al., 2012).

In recent years, a small subclass of SNe Ia has been identified as interacting
with a circumstellar medium (CSM) in the months to years post-explosion
(Silverman et al., 2013; Graham et al., 2019). While the number of members in
this class is growing, it is still quite small, and systematic searches in the
radio (Chomiuk et al., 2016; Cendes et al., 2020) and the UV (Dubay et al.,
2022) have turned up only a few examples. In the SNR field, only Kepler’s SNR
(Reynolds et al., 2007) in the Galaxy and N103B in the LMC (Williams et al.,
2014) have shown evidence of substantial circumstellar interaction (at radii
of a few pc and ages of several centuries).

Of all the Balmer-dominated SNRs in the LMC, SNR 0519-69.0 (hereafter, 0519)
is one of the less well studied. The light-echo of the supernova that gave
rise to 0519 presents as that of a spectroscopically normal SN Ia and the
derived explosion age is $(600\pm 200)\,\mathrm{years}$ (Rest et al., 2005).
At least several solar masses of material have been swept-up by the blast wave
so far (Borkowski et al., 2006). Models taking the combined constraints of the
age and the forward and reverse shocks (as traced by the broad [Fe xiv]
emission) into account are consistent with a Chandrasekhar-mass progenitor
(Seitenzahl et al., 2019). Searches for surviving companions however have not
been able to identify a donor/companion star (Edwards et al., 2012), but see
Li et al. (2019) for a proposed candidate. A super-soft X-ray source has been
ruled out a progenitor for 0519 by Kuuttila et al. (2019), based on the
absence of a relic ionization nebular surrounding 0519. Li et al. (2021)
recently reported on the discovery of some higher density knots exhibiting low
surface-brightness forbidden line emission in 0519.

Although broad and narrow H$\alpha$ were first detected by Smith et al.
(1991), only the brightest filaments (corresponding to the slower shocks) were
observed in that study. Ghavamian et al. (2007) detected broad Ly $\beta$
emission in the UV with the Far Ultraviolet Spectroscopic Explorer, from which
they estimated a shock velocity of 3600$-$7100 km/s. More recently, Hovey et
al. (2018) used a combination of ground-based optical spectra and space-based
HST proper motions of 0519 to estimate the forward shock speed and cosmic ray
acceleration efficiency of the SNR. Hovey et al. (2018) used HST observations
separated by 1 year (2010 and 2011) to make their proper motion measurements.
In this work, we extend the temporal baseline by an order of magnitude, using
new observations obtained in 2020 to substantially reduce the uncertainties on
the proper motions of the various filaments.

##  2 Hubble Space Telescope Observations

The data considered here are a combination of existing archival HST images of
0519 acquired in 2010 and 2011 (Program ID GO-12017) with the ACS/F658N
instrument/filter combination, along with a new epoch of data acquired in 2020
with WFC3/F657N as part of our approved Chandra program (Program ID GO-15989,
which was part of a joint Chandra-HST program; B. J. Williams, PI). The
WFC3/UVIS data were obtained on 2020 Jun 21 and Aug 10 with the F657N
(H$\alpha$) filter, as shown in Table 1. The observations obtained in each
exposure were dithered to permit removal of artifacts from the data and (for
UVIS) to cover the chip gap, and the appropriate FLASH parameter was set for
UVIS to reduce the effects of charge transfer inefficiency. The combined
exposure was 8946 s for the WFC3 data. The corresponding earlier epoch ACS
exposures were 4757 s each for both the 2010 and 2011 observations.

Table 1: HST/ACS and WFC3 Data Used in Expansion Measurements Instrument | Dataset | Filter | Exposure Time (s) | Observation Date
---|---|---|---|---
ACS/WFC2 | JBDQ01010 | F658N | 3660 | 2010-04-17
ACS/WFC2 | JBDQ01E5Q | F658N | 1097 | 2010-04-17
ACS/WFC2 | JBDQ02010 | F658N | 3660 | 2011-04-21
ACS/WFC2 | JBDQ02Y7Q | F658N | 1097 | 2011-04-21
WFC3/UVIS2 | IE6M02020 | F657N | 3042 | 2020-06-21
WFC3/UVIS2 | IE6M52010 | F657N | 5914 | 2020-08-10
Figure 1: Difference image between images of 0519 taken at two epochs:
2010+2011 (ACS F657N images) minus 2020 (WFC3 F658N, rescaled to match ACS).
The black and white striping along the edge of the blast wave shows sharp,
distinct regions of positive and negative pixels, reflecting the substantial
expansion movement during the 10-year baseline.

The HST imagery provides exactly a decade in baseline (2010 to 2020), which,
combined with the high shock velocities in 0519 ($\sim$3000 km/s; Hovey et al.
(2018)), should allow a precise, easily measured expansion measurement. The
expansion rate in arcseconds per year, $\omega$, is given by

$\omega\,=(2.1\times 10^{-4})\frac{v_{s}}{D}$ (1)

where $v_{s}$ is the shock speed (km/s) and $D$ is the SNR distance in kpc.
Using a distance of 50 kpc for the LMC (Clementini et al. (2003)), shock
speeds of $2000-3000\,\mathrm{km/s}$ would result in filament displacements of
around 0.08″- 0.1″over a 10 year baseline, corresponding to 4-5 ACS/WFC3
pixels. Owing to the small expected displacement ($<$ half a pixel) of the
Balmer filaments expected between 2010 and 2011 (as compared to the 2020 WFC3
data), we combined the 2010 and 2011 observations into a single high signal-
to-noise image. Before combining, we utilized the tweakreg task in Drizzlepac
to align the WCS solutions of all WFC3 images and ACS images to one another
separately. The images include numerous well exposed LMC stars, with some near
the full-well saturation limit, so the PEAKMAX value was set to 70,000
electrons during tweakreg runs. A SEARCHRAD parameter was set to 3.0″to allow
a generous offset for search radius. We then used Astrodrizzle 3.2.1 to
combine the 2010 and 2011 ACS images. As per our calculation above, the Balmer
filaments in 0519 will only have moved less than 1/5 ACS pixels between the
2010 and 2011 data, so combining the 2010 and 2011 data will result in far
smaller smearing of the filament emission than expected between the ACS data
acquisition and the WFC3 data acquisition. Finally, we used Astrodrizzle to
rotate and align the 2020 WFC3 image with the combined 2010+2011 ACS image.
This procedure also resampled the WFC3 data to the ACS pixel scale (0.05′′
pixel-1).

Examination of the final WFC3 and ACS images showed no residual asymmetries in
the stellar PSFs, verifying that the image alignment had been successful. We
checked the absolute astrometric accuracy of the final combined images by
measuring centroids for a handful of stars in each field and comparing their
coordinates to those of of those stars in the Gaia DR2 catalogue (Gaia
Collaboration et al. (2018)). These coordinates matched to within
0.01′′-0.02′′, or less than a quarter of an ACS pixel. In Figure 1, we show a
difference image between the two epochs, where virtually all filaments can be
easily seen to have moved outward.

##  3 Measurements & Discussion

Our procedure for measuring the proper motions of the filaments is identical
to that used in our previous proper motion studies of various remnants in the
optical, radio, and X-ray bands (Williams et al., 2016; Winkler et al., 2014;
Coffin et al., 2021). Briefly, we extract the 1D radial profiles from both
epochs (where “radial” is defined as locally perpendicular to the shock front,
which is linear over the small scales that we examine), with uncertainties on
each data point, then shift epoch 1 relative to epoch 2, minimizing the value
of $\chi^{2}$ when subtracting one profile from the other. Profiles are
extracted in pixel space (using the ACS pixel scale size of 0.05′′), with
shifts calculated on a grid of 10000 points with a step size of 0.002 pixels.
We smooth the profiles very slightly, using a 1-pixel Gaussian. This has
virtually no effect on the profile shapes, but significantly decreases the
pixel-to-pixel Poisson noise level. Uncertainties reported are the 90%
confidence limits from the $\chi^{2}$ fitting procedure.

Figure 2: The WFC3 F658N image of 0519. The locations of our expansion
measurements around the rim are marked by the numbers. Proper motion values
for the 21 regions are given in Table 2.

Owing to the high angular resolution of HST and the high signal-to-noise ratio
of the filaments, we were able to choose very narrow regions (10 pixels wide,
or 0.5′′) over which to extract our 1D profiles. Such narrow regions ensure
that filament curvature is negligible, while also allowing us to easily avoid
stars. In Figure 2, we show the 21 regions we selected for analysis. This
choice of regions is of course somewhat arbitrary, but we followed two general
principles: 1) we chose enough regions around the remnant to get a
representative sample size of forward shock velocities; and 2) in places like
regions 3$-$6 or 17 and 18, where several distinct filamentary structures are
present, we measure each filament’s position to search for differences in the
expansion velocity. Region 1, 2, 8, and 10 are all located on faint “blowout”
regions, visible on the image. Results from our measurements are shown in
Table 2 and Figure 3. All values are also converted to velocities, using a
distance of 50 kpc to the LMC.

Figure 3: Forward shock velocities for 0519 derived from the expansion
measurements, shown for the 21 regions around the rim. Error bars are marked
by the vertical lines. In general, the fastest shocks correspond to the
faintest filaments, with the slowest shocks correspond to the brightest
filaments, consistent with greater deceleration in denser material. See the
discussion in the text.

It is apparent from the measurements that there is a spread in velocities
present throughout the remnant. In general, the slowest shock speeds (regions
5, 6, and 14) correspond to the brightest filaments, while the faster shock
speeds (regions 1, 2, 8, 18, 19) arise from the fainter filaments. Taking the
average numbers for those groupings of filaments, we find an average shock
velocity for the slower filament group of 1670 km s-1, and an average for the
faster group of 5280 km s-1, corresponding to a factor of around 3 range in
shock speed. Under the assumption of pressure equilibrium, where the product
of density, $\rho$, and the square of the shock velocity, $v$ (i.e., $\rho
v^{2}$) is constant, a factor of 3.1 difference in the shock velocities
implies about a factor of 10 difference in the ambient density.

The flux ratio between the groupings of faint filaments and bright filaments
measured above, as measured from the HST image, is about a factor of 5.
Because the H$\alpha$ emission behind the shock scales as the flux of incoming
atoms $\rho v$, a velocity contrast of 3 and a flux contrast of 5 imply a
density contrast of 15. Given the simplicity of these estimates, the rough
agreement between the factor of 10 inferred from pressure equilibrium
arguments and the factor of 15 inferred from H$\alpha$ brightness is
encouraging, and suggests that the remnant is not expanding into uniform
density material.

It is important to note that due to projection effects, the proper motion can
only measure the velocity of a filament in the plane of the sky. Getting the
true velocity requires having the line-of-sight velocity via spectroscopic
measurements. Geometrically, this is likely to be more important for the more
internal filaments than those on the rim. In a future work, we will explore
spectroscopy of these filaments, which can reveal both the bulk velocity and
the widths of the broad spectral components.

Figure 4: Left: An overlay of the 0.4-0.7 keV X-ray emission from 0519-69.0 in
red, the 0.7-0.8 keV emission in blue, and the H$\alpha$ emission in green,
from the 2000 observation. The 0.4-0.7 keV energy band is most sensitive to O
emission (ambient medium), while the 0.7-0.8 keV band is most sensitive to Fe
L-shell emission (primarily ejecta). While both are present throughout, there
is a clear enhancement of the O emission in red in the regions corresponding
to the brightest, slowest, optical filaments. Middle: The X-ray emission in
hard X-rays from 3.0-6.2 keV is shown in purple, with the H$\alpha$ emission
in green. In all energy bands, there is little to no correlation of X-ray
emission with the fastest shocks, consistent with those coming from only the
lowest density areas. Right: The de-convolved 24 $\mu$m Spitzer image is shown
in purple, with the H$\alpha$ in green. See discussion in the text.

One possible explanation for the large density variations on small scales is
the presence of a dense shell of circumstellar material which has, in some
locations, been penetrated by either the blast wave or clumps of ejecta.
Regions 8 and 10 are examples. We have examined in detail the correlation
between soft X-ray emission (0.4-0.8 keV) and H$\alpha$ emission. Due to the
degradation of Chandra’s soft X-ray response over the years, we use the
original Chandra observation of this remnant from June 2000 (Obs. ID 118, PI
S. Holt). In Figure 4, we show images of the remnant in the X-ray and optical
bands plotted together. Sharp edges in the X-ray emission match slower shocks
with brighter H$\alpha$ emission quite well (e.g., regions 5, 6, 14, and 20).
In the left panel of Figure 4, we separate the soft X-ray emission into two
components: the 0.4-0.7 keV emission, most sensitive to lines from O, and the
0.7-0.8 keV emission, most sensitive to emission from L-shell emission of Ne-
like Fe ions.

The O emission is most prominent in those regions where the shock speeds are
lowest. As was shown in Borkowski et al. (2006), these regions also correspond
to the brightest knots of emission seen in Spitzer 24 $\mu$m images. The
authors there concluded that the density in those knots is at least several
times higher than in other parts of the remnant. In Figure 4, we show an
overlay of the Spitzer 24 $\mu$m emission with the HST H$\alpha$ emission.
This image has been de-convolved to a resolution of $\sim 2.5$” (the native
resolution of Spitzer’s 24 $\mu$m camera is $\sim 7$”) using the iCore
software package (Masci & Fowler, 2009). Similarly to the X-ray emission, the
dust emission in IR is brightest where the H$\alpha$ emission is bright and
the shocks are slowest. As has been shown in previous studies of this and
other remnants with Spitzer (Borkowski et al., 2006; Williams et al., 2011),
it is possible to measure the post-shock gas density using emission from warm
dust grains, as grain temperature (and thus spectral shape) and overall IR
luminosity both scale with density. However, even with de-convolution,
Spitzer’s spatial resolution makes this difficult to do for LMC remnants,
where the angular sizes of SNRs is small ($\leq$1′). Future high-resolution
observations with JWST will be critical for gaining a detailed understanding
of the density structures encountered by 0519 and other Type Ia SNRs in the
LMC.

Figure 5: Simulated RGB image of Tycho at 50 kpc, as described in text, and
based on a 309-second portion of a 2004 Chandra observation (ObsID 3837), not
exposure-corrected. The chip gaps from the ACIS-I instrument are evident. Red:
0.5 – 1.8 keV; green, 1.8 – 4.0 keV; blue, 4.0 – 6.2 keV. Note the bluish rim
evident around most of the periphery; in the blue energy range, free of strong
lines, the synchrotron continuum is expected to be most evident.

In contrast, there is little or no X-ray emission present from the faster,
fainter filaments (e.g., 1, 2, 8, 15, 18), even in the hard X-rays, as Figure
4 shows. If the density contrast really is a factor of 10-15, then the thermal
X-ray emission (X-ray emissivity scales as $\rho^{2}$) would be fainter by a
factor of 100-200. In fact, this is almost exactly what we measure: when
comparing the X-ray count rate in a small area behind region 18 with an
identical area behind region 20, we find the X-ray brightness to be higher by
a factor of 105 in the latter. It appears that the X-ray emission in 0519 is
dominated by interaction with the densest material. We conclude that 0519, in
at least some places, is interacting with a dense CSM. This conclusion is
supported by the recent work of Li et al. (2021), who report that the dense
knots seen in forbidden line observations of 0519 are consistent with an
origin in a dense CSM. The dense knots seen in that paper correspond to the
brightest knots in H$\alpha$ seen in our observations.

All the shock velocities we obtain are well in excess of the $\sim 1,000$ km
s-1 value often taken as the rough lower limit for shock speeds capable of
accelerating electrons to X-ray-synchrotron-emitting energies $\gtrapprox 1$
TeV. However, X-ray images (shown in Figure 4, which use the recent deep
Chandra observations of $\sim 400$ ks) show no evidence of non-thermal
filamentary counterparts to the H$\alpha$ filaments that circumscribe 0519.
This is perhaps not surprising; the young Galactic SN Ia remnants Kepler,
Tycho, and SN 1006 all show Balmer-dominated shocks in some areas, but
typically without obvious X-ray counterparts. In SN 1006, dominated by
nonthermal X-ray emission, faint H$\alpha$ emission can be followed almost all
around the remnant periphery, and in a few locations, narrow X-ray rims become
evident as H$\alpha$ rim emission weakens, overlapping over tens of degrees in
azimuth. However, in the NW, where H$\alpha$ emission is strongest, there is
no nonthermal X-ray emission (see images in Winkler et al. (2014)). The
absence of nonthermal X-ray emission at the location of the H$\alpha$
filaments is somewhat expected on theoretical grounds, in that the partially
neutral upstream medium required for the Balmer emission is likely to suppress
electron acceleration to TeV energies, due to ion-neutral wave damping (Draine
& McKee, 1993; O’C Drury et al., 1996). However, a cosmic-ray precursor can be
inferred to be present ahead of Balmer-dominated shocks through detailed
optical spectroscopy (this and related topics are reviewed in Heng (2010) and
Ghavamian et al. (2013); also see Knežević et al. (2017) for more recent
observational results). The substantial Balmer line emission from 0519,
surrounding the entire remnant, is stronger and more extensive than in
Galactic remnants, perhaps indicating a higher neutral fraction upstream and
accounting for the absence of prominent nonthermal X-ray rims.

One may reasonably ask whether detection of thin rims of nonthermal X-ray
emission, such as those seen in Tycho’s SNR, is possible at a distance of 50
kpc. To simulate this, we took a single observation of Tycho from 2004 (Obs ID
3837), and assumed its distance to be 2.3 kpc. Translating this to 50 kpc
leads to angular scales smaller by a factor of 21.7 and fluxes smaller by a
factor of 472. We chose a random piece of the Obs ID 3837 observation (total
length 146 ks) that was only 0.309 ks in duration (146/472). We then binned
this image by a factor of 22, and generated an RGB image with red (0.5 - 1.8
keV), green (1.8-4.0 keV), and blue (4.0-6.2 keV). We show this image in
Figure 5, where one can see the narrow rims of nonthermal emission in blue.
While this is only qualitative, it shows that it should at least be possible
to discern nonthermal rims in 0519, if they were present. In a future work, we
will examine the X-ray observations in more detail.

At a distance of 50 kpc, the radius of the remnant of $\sim 15^{\prime\prime}$
corresponds to a physical radius of $\sim 3.6$ pc. Using the average shock
velocities for the fastest filaments reported above (including their
uncertainties), we can derive an undecelerated age of $\sim 670\pm 70$ years.
This assumes that the explosion site is in the center of the remnant, an
assumption that may not be valid, even for “circular” SNRs, though the effect
of this is relatively small (Williams et al., 2013). Because some amount of
deceleration has almost certainly occurred, the real age is somewhat younger
than this, which is completely consistent with the age reported in Rest et al.
(2005) of $600\pm 200$ years. Continued follow-up observations with HST will
further refine estimates of the shock velocity.

##  4 Conclusions

We re-observed 0519-69.0, a young Type Ia SNR in the LMC, using HST in 2020, a
decade after it was first observed in 2010. With this long baseline, we have
measured the proper motions of $\sim 20$ filaments in the remnant that
constitute a representative sample of the motions present. With the known
distance to the LMC, we convert these proper motions to absolute velocities,
finding a large range of about a factor of three. Combining our inferred
velocities with a morphological study of the remnant at optical, IR, and X-ray
wavelengths, we find that the slowest shock velocities result from an
interaction of the blast wave with substantially denser material than the
average ISM densities in the LMC ($\sim 0.1$ cm-3). This dense material likely
results from a circumstellar medium. We derive an undecelerated age (an upper
limit) for the remnant of 670 $\pm 70$ yr, consistent with previous studies.
Young remnants like this should be continuously monitored.

B.J.W. and P.G. acknowledge funding support through HST grant HST-GO-15989.

Table 2: Expansion Measurements of 0519 (2010-2020) Region Number | Proper MotionaaAll proper motions are reported in milliarcseconds over a 10 year baseline. | PM Lower Limit | PM Upper Limt | Shock VelocitybbAll shock velocities reported in km/s and assume a distance of 50 kpc for LMC | SV Lower Limit | SV Upper Limit
---|---|---|---|---|---|---
1 | 196 | 156 | 234 | 4890 | 3890 | 5840
2 | 162 | 124 | 212 | 4040 | 3090 | 5290
3 | 107 | 88 | 126 | 2670 | 2200 | 3140
4 | 104 | 93 | 115 | 2600 | 2320 | 2870
5 | 59 | 55 | 63 | 1470 | 1370 | 1570
6 | 79 | 77 | 81 | 1970 | 1920 | 2020
7 | 145 | 140 | 154 | 3670 | 3490 | 3840
8 | 200 | 159 | 243 | 4990 | 3970 | 6060
9 | 139 | 130 | 148 | 3470 | 3240 | 3690
10 | 96 | 60 | 144 | 2400 | 1500 | 3590
11 | 148 | 139 | 156 | 3690 | 3470 | 3890
12 | 71 | 63 | 80 | 1770 | 1570 | 2000
13 | 106 | 98 | 114 | 2650 | 2450 | 2840
14 | 63 | 59 | 66 | 1570 | 1470 | 1650
15 | 160 | 108 | 229 | 3990 | 2700 | 5710
16 | 171 | 155 | 185 | 4270 | 3870 | 4620
17 | 110 | 97 | 123 | 2750 | 2420 | 3070
18 | 244 | 208 | 299 | 6090 | 5190 | 7460
19 | 206 | 183 | 232 | 5140 | 4570 | 5790
20 | 121 | 117 | 126 | 3020 | 2920 | 3140
21 | 170 | 148 | 193 | 4240 | 3690 | 4820

## References

  * Borkowski et al. (2006) Borkowski, K. J., Williams, B. J., Reynolds, S. P., et al. 2006, ApJ, 642, L141, doi: 10.1086/504472
  * Cendes et al. (2020) Cendes, Y., Drout, M. R., Chomiuk, L., & Sarbadhicary, S. K. 2020, ApJ, 894, 39, doi: 10.3847/1538-4357/ab6b2a
  * Chevalier et al. (1980) Chevalier, R. A., Kirshner, R. P., & Raymond, J. C. 1980, ApJ, 235, 186, doi: 10.1086/157623
  * Chomiuk et al. (2016) Chomiuk, L., Soderberg, A. M., Chevalier, R. A., et al. 2016, ApJ, 821, 119, doi: 10.3847/0004-637X/821/2/119
  * Clementini et al. (2003) Clementini, G., Gratton, R., Bragaglia, A., et al. 2003, AJ, 125, 1309, doi: 10.1086/367773
  * Coffin et al. (2021) Coffin, S. C., Williams, B. J., & Katsuda, S. 2021, arXiv e-prints, arXiv:2112.08448.  https://arxiv.org/abs/2112.08448
  * Draine & McKee (1993) Draine, B. T., & McKee, C. F. 1993, ARA&A, 31, 373, doi: 10.1146/annurev.aa.31.090193.002105
  * Dubay et al. (2022) Dubay, L. O., Tucker, M. A., Do, A., Shappee, B. J., & Anand, G. S. 2022, ApJ, 926, 98, doi: 10.3847/1538-4357/ac3bb4
  * Edwards et al. (2012) Edwards, Z. I., Pagnotta, A., & Schaefer, B. E. 2012, ApJ, 747, L19, doi: 10.1088/2041-8205/747/2/L19
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Ghavamian et al. (2007) Ghavamian, P., Blair, W. P., Sankrit, R., Raymond, J. C., & Hughes, J. P. 2007, ApJ, 664, 304, doi: 10.1086/518686
  * Ghavamian et al. (2001) Ghavamian, P., Raymond, J., Smith, R. C., & Hartigan, P. 2001, ApJ, 547, 995, doi: 10.1086/318408
  * Ghavamian et al. (2013) Ghavamian, P., Schwartz, S. J., Mitchell, J., Masters, A., & Laming, J. M. 2013, Space Sci. Rev., 178, 633, doi: 10.1007/s11214-013-9999-0
  * Graham et al. (2019) Graham, M. L., Harris, C. E., Nugent, P. E., et al. 2019, ApJ, 871, 62, doi: 10.3847/1538-4357/aaf41e
  * Heng (2010) Heng, K. 2010, PASA, 27, 23, doi: 10.1071/AS09057
  * Hovey et al. (2018) Hovey, L., Hughes, J. P., McCully, C., Pandya, V., & Eriksen, K. 2018, ApJ, 862, 148, doi: 10.3847/1538-4357/aac94b
  * Knežević et al. (2017) Knežević, S., Läsker, R., van de Ven, G., et al. 2017, ApJ, 846, 167, doi: 10.3847/1538-4357/aa8323
  * Kuuttila et al. (2019) Kuuttila, J., Gilfanov, M., Seitenzahl, I. R., Woods, T. E., & Vogt, F. P. A. 2019, MNRAS, 484, 1317, doi: 10.1093/mnras/stz065
  * Li et al. (2021) Li, C.-J., Chu, Y.-H., Raymond, J. C., et al. 2021, ApJ, 923, 141, doi: 10.3847/1538-4357/ac2c04
  * Li et al. (2019) Li, C.-J., Kerzendorf, W. E., Chu, Y.-H., et al. 2019, ApJ, 886, 99, doi: 10.3847/1538-4357/ab4a03
  * Long et al. (1981) Long, K. S., Helfand, D. J., & Grabelsky, D. A. 1981, ApJ, 248, 925, doi: 10.1086/159222
  * Masci & Fowler (2009) Masci, F. J., & Fowler, J. W. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler, 67.  https://arxiv.org/abs/0812.4310
  * Morlino et al. (2012) Morlino, G., Bandiera, R., Blasi, P., & Amato, E. 2012, ApJ, 760, 137, doi: 10.1088/0004-637X/760/2/137
  * Morlino et al. (2013) Morlino, G., Blasi, P., Bandiera, R., & Amato, E. 2013, A&A, 557, A142, doi: 10.1051/0004-6361/201322161
  * O’C Drury et al. (1996) O’C Drury, L., Duffy, P., & Kirk, J. G. 1996, A&A, 309, 1002.  https://arxiv.org/abs/astro-ph/9510066
  * Rest et al. (2005) Rest, A., Suntzeff, N. B., Olsen, K., et al. 2005, Nature, 438, 1132, doi: 10.1038/nature04365
  * Reynolds et al. (2007) Reynolds, S. P., Borkowski, K. J., Hwang, U., et al. 2007, ApJ, 668, L135, doi: 10.1086/522830
  * Seitenzahl et al. (2019) Seitenzahl, I. R., Ghavamian, P., Laming, J. M., & Vogt, F. P. A. 2019, Phys. Rev. Lett., 123, 041101, doi: 10.1103/PhysRevLett.123.041101
  * Silverman et al. (2013) Silverman, J. M., Nugent, P. E., Gal-Yam, A., et al. 2013, ApJS, 207, 3, doi: 10.1088/0067-0049/207/1/3
  * Smith et al. (1991) Smith, R. C., Kirshner, R. P., Blair, W. P., & Winkler, P. F. 1991, ApJ, 375, 652, doi: 10.1086/170228
  * Tuohy et al. (1982) Tuohy, I. R., Dopita, M. A., Mathewson, D. S., Long, K. S., & Helfand, D. J. 1982, ApJ, 261, 473, doi: 10.1086/160358
  * van Adelsberg et al. (2008) van Adelsberg, M., Heng, K., McCray, R., & Raymond, J. C. 2008, ApJ, 689, 1089, doi: 10.1086/592680
  * Williams et al. (2013) Williams, B. J., Borkowski, K. J., Ghavamian, P., et al. 2013, ApJ, 770, 129, doi: 10.1088/0004-637X/770/2/129
  * Williams et al. (2016) Williams, B. J., Chomiuk, L., Hewitt, J. W., et al. 2016, ApJ, 823, L32, doi: 10.3847/2041-8205/823/2/L32
  * Williams et al. (2011) Williams, B. J., Blair, W. P., Blondin, J. M., et al. 2011, ApJ, 741, 96, doi: 10.1088/0004-637X/741/2/96
  * Williams et al. (2014) Williams, B. J., Borkowski, K. J., Reynolds, S. P., et al. 2014, ApJ, 790, 139, doi: 10.1088/0004-637X/790/2/139
  * Winkler et al. (2014) Winkler, P. F., Williams, B. J., Reynolds, S. P., et al. 2014, ApJ, 781, 65, doi: 10.1088/0004-637X/781/2/65


# Nanomechanically-induced nonequilibrium quantum phase transition
to a self-organized density wave of a Bose-Einstein condensate

Milan Radonjić<EMAIL_ADDRESS>I. Institut für
Theoretische Physik, Universität Hamburg, Notkestr. 9, 22607 Hamburg, Germany
Institute of Physics Belgrade, University of Belgrade, Pregrevica 118, 11080
Belgrade, Serbia  Leon Mixa<EMAIL_ADDRESS>I. Institut für
Theoretische Physik, Universität Hamburg, Notkestr. 9, 22607 Hamburg, Germany
The Hamburg Center for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg,
Germany  Axel Pelster<EMAIL_ADDRESS>Physics Department and Research
Center OPTIMAS, University Kaiserslautern-Landau, Erwin-Schrödinger Str. 46,
67663 Kaiserslautern, Germany  Michael Thorwart<EMAIL_ADDRESS>hamburg.de I. Institut für Theoretische Physik, Universität Hamburg, Notkestr.
9, 22607 Hamburg, Germany  The Hamburg Center for Ultrafast Imaging, Luruper
Chaussee 149, 22761 Hamburg, Germany

###### Abstract

We report on a nonequilibrium quantum phase transition (NQPT) in a hybrid
quantum many-body system consisting of a vibrational mode of a damped
nanomembrane interacting optomechanically with a cavity, whose output light
couples to two internal states of an ultracold Bose gas held in an external
quasi-one-dimensional box potential. For small effective membrane-atom
couplings, the system is in a homogeneous Bose-Einstein condensate (BEC)
steady state, with no membrane displacement. Depending on the transition
frequency between the two internal atomic states, either one or both internal
states are occupied. By increasing the atom-membrane couplings, the system
transitions to a symmetry-broken self-organized BEC phase, which is
characterized by a considerably displaced membrane steady-state and density
wave-like BEC profiles. This NQPT can be both discontinuous and continuous for
a certain interval of transition frequencies, and is purely discontinuous
outside of it.

##  I Introduction

Many physical systems can be classified in terms of their emergent collective
behavior when they undergo phase transitions [1, 2, 3]. While equilibrium
phase transitions are fairly well understood, nonequilibrium quantum phase
transitions (NQPTs) are a relatively new and exciting field of research [4,
5]. Such phase transitions occur in a nonequilibrium steady state of driven
quantum systems in contact with an environment, so that their energy is not
conserved due to external driving and dissipation [6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16]. Among other cases, NQPTs have been studied in ultracold atoms in
an optical cavity [17, 18, 19, 20, 21, 22, 23, 24], in driven-dissipative
Bose-Einstein condensates [25, 26], in microcavity-polariton systems [27, 28,
29], and in photon Bose-Einstein condensates [30, 31, 32]. Additionally,
dynamical quantum phase transitions and universal scaling have been observed
both in the nonequilibrium dynamics of isolated quantum systems after a
quench, with time playing the role of the control parameter [33], and recently
in open quantum systems, whose dynamics is driven by the dissipative contact
to an environment [34].

Hybrid atom-optomechanical systems, which combine optomechanics with atom
optics, represent one recent promising platform for studying NQPTs [35, 36,
37, 38, 39, 40, 41, 42, 43]. In these systems, a nanomembrane in an optical
cavity is coupled to the motional [35] or to the internal [36] degrees of
freedom of a distant cloud of cold atoms that are trapped in the optical
lattice of the outcoupled light field. By employing different cooling
mechanisms, such as sympathetic cooling [35, 37, 38, 39] and optical feedback
cooling [39, 44, 45], the nanomembrane can be cooled down almost to its
quantum ground state. Moreover, hybrid atom-optomechanical systems are a
versatile playground for various quantum phenomena, such as indirect quantum
measurement, atom-membrane entanglement and coherent state transfer [46, 47,
48, 49, 50].

The internal state coupling scheme [36] is realized by a nanomembrane that is
coupled to transitions between internal states of the atoms via translating
the phase shift of the light, caused by the membrane displacement, into a
polarization rotation using a polarizing beam splitter. This scheme allows for
resonant coupling and mitigates the drawbacks of the motional coupling scheme
[35], like strong frequency mismatch between the nanooscillator and the atomic
motion in the optical trap. It has been utilized for membrane cooling [36,
51], displacement squeezing [52], the realization of a effective negative mass
oscillator [53], and quantum back-action evading measurements [54]. In
addition, a peculiar NQPT, whose order can be controlled by changing a
directly accessible experimental parameter, has recently been proposed in this
scheme [55]. It has been found that the system can undergo both a first- and a
second-order phase transition in the same physical setup, by tuning the atomic
transition frequency.

In this work, we examine NQPTs in a hybrid atom-optomechanical system in the
internal state coupling scheme. We consider an ultracold atomic gas trapped in
a large quasi-1D box, where the lattice potentials generated by the light
fields have been canceled. Since the membrane is naturally damped by its
suspension, the atoms relax into a steady state whose features critically
depend on the state of the membrane. For large enough trapping box and weak
effective coupling between the atoms and the nanomembrane, we find that the
atomic condensates are essentially uniform and the membrane is not displaced.
When the atom-membrane couplings become strong enough, the system transitions
into a symmetry-broken self-organized phase that exhibits density wave-like
condensate profiles and significant nanomembrane displacement. Depending on
the values of the readily controllable experimental parameters, this NQPT can
be of either first or second order. We map out the relevant parts of the
parameter space and study the hallmarks of the phase transition. Hence, the
rest of the paper is organized as follows. In Sec. II we outline the basic
mean-field theory. In the following Sec. III we examine the homogeneous steady
states. Next, in Sec. IV we focus on the inhomogeneous steady state and its
features. We also address and examine the occurring 1st and 2nd order phase
transitions. We discuss our results in Sec. V and provide an outlook.

##  II Model and equations of motion

The system under consideration [36] consists of $N$ identical ultracold
bosonic atoms held in an external quasi-1D trap (see Fig. 1). The relevant
atomic states $\\{|-\rangle,|+\rangle,|e\rangle\\}$ constitute a $\Lambda$
level scheme, where the two lower states are energetically separated by the
atomic transition frequency $\Omega_{\rm a}$. An applied $\sigma_{-}$
circularly polarized laser of frequency $\omega_{\rm L}$ drives the internal
transition $|+\rangle\leftrightarrow|e\rangle$ at a finite large detuning
$\Delta$. The passing beam is sent to a polarizing beam splitter (PBS), which
splits the circularly polarized light into linearly polarized $\pi_{x}$ and
$\pi_{y}$ light beams in two perpendicular directions. In the vertical arm, a
fixed mirror simply reflects light with conserved polarization $\pi_{x}$. The
horizontal beam enters a low-finesse cavity which hosts a semitransparent
nanomembrane with the mechanical resonance frequency $\Omega_{\rm m}$. The
lengths of the two arms are made equal when the membrane is undisplaced. In
that case, the cavity reflects $\pi_{y}$ light without any relative phase
shift. In a quasistatic picture, any finite displacement of the membrane
induces a relative phase shift between the $\pi_{x}$ and $\pi_{y}$ light beams
from the two arms. This relative phase is converted to a polarization rotation
after the two beams have passed backwards through the PBS. The emerging
$\sigma_{+}$ polarized light then drives the internal atomic transition
$|-\rangle\leftrightarrow|e\rangle$ and may induce a two-photon Raman
transition between the two lower states of the $\Lambda$ scheme via the
excited state $|e\rangle$. On the other hand, an internal transition between
the two lower atomic states leads to the emission of a $\sigma_{+}$ polarized
photon which alters the radiation pressure on the membrane after reaching it,
which leads to an effective atom-membrane coupling.

Figure 1: A semitransparent nanomembrane in an optical cavity is coupled to
the internal states of a distant atomic ensemble that is trapped in an
external quasi-1D trap. The internal states of the atoms constitute a
$\Lambda$ level scheme according to the inset. The figure is adapted from Ref.
[55].

An adiabatic elimination of the excited state and the light field results [55]
in an effective Hamiltonian involving the membrane and the low-energy atomic
degrees of freedom ($\hbar=c=1$)

$\displaystyle\hat{H}_{\rm ef{}f}$ $\displaystyle=\Omega_{\rm
m}\hat{a}^{\dagger}\hat{a}+\sum_{\nu=\pm}\int
dz\;\\!\hat{\psi}_{\nu}^{\dagger}(z)\bigg{[}\nu\frac{\Omega_{\rm
a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{\nu}(z)+\frac{1}{2}\sum_{\nu^{\prime}=\pm}g_{\nu\nu^{\prime}}\hat{\psi}_{\nu^{\prime}}^{\dagger}(z)\hat{\psi}_{\nu^{\prime}}(z)\bigg{]}\hat{\psi}_{\nu}(z)$
$\displaystyle-(\hat{a}+\hat{a}^{\dagger})\int
dz\>\\!\sin(2z)\bigg{[}\frac{\lambda}{2}\big{(}\hat{\psi}_{+}^{\dagger}(z)\hat{\psi}_{-}(z)+\hat{\psi}_{-}^{\dagger}(z)\hat{\psi}_{+}(z)\big{)}+\lambda_{\rm
ex}\hat{\psi}_{+}^{\dagger}(z)\hat{\psi}_{+}(z)\bigg{]},$ (1)

where the position coordinate has been re-scaled as $z\to z/\omega_{\rm L}$.
The bosonic operator $\hat{a}$ denotes the annihilation operator of the
excitations of one single mechanical mode of the nanomenbrane, while
$\hat{\psi}_{\nu}(z)$ denotes the field operator of the atoms in the state
$\nu$. The atomic recoil frequency reads $\omega_{\rm R}=\omega_{\rm
L}^{2}/(2m)$, the external trapping potentials are denoted by $V_{\pm}(z)$,
while $g_{\nu\nu^{\prime}}$ are atomic state-dependent contact interaction
strengths, and the effective membrane-atom coupling strengths are $\lambda$
and $\lambda_{\rm ex}$. Note that the effective cavity-mediated global atom-
atom interaction [56] has been neglected. Taking into account that the
membrane is damped at the rate $\Gamma_{\rm m}$, we derive the following
equations of motion for the membrane and the atomic operators

$\displaystyle i\partial_{t}\hat{a}$ $\displaystyle=\big{(}\Omega_{\rm
m}-i\Gamma_{\rm m}\big{)}\hat{a}-\int
dz\>\\!\sin(2z)\bigg{[}\frac{\lambda}{2}\big{(}\hat{\psi}_{+}^{\dagger}(z)\hat{\psi}_{-}(z)+\hat{\psi}_{-}^{\dagger}(z)\hat{\psi}_{+}(z)\big{)}+\lambda_{\rm
ex}\hat{\psi}_{+}^{\dagger}(z)\hat{\psi}_{+}(z)\bigg{]}-\hat{\xi}_{\rm m},$
(2a) $\displaystyle i\partial_{t}\hat{\psi}_{-}(z)$
$\displaystyle=\bigg{[}-\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{-}(z)+\sum_{\nu=\pm}g_{-\nu}\hat{\psi}_{\nu}^{\dagger}(z)\hat{\psi}_{\nu}(z)\bigg{]}\hat{\psi}_{-}(z)-\frac{\lambda}{2}\sin(2z)\big{(}\hat{a}+\hat{a}^{\dagger}\big{)}\hat{\psi}_{+}(z),$
(2b) $\displaystyle i\partial_{t}\hat{\psi}_{+}(z)$
$\displaystyle=\bigg{[}\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{+}(z)+\sum_{\nu=\pm}g_{+\nu}\hat{\psi}_{\nu}^{\dagger}(z)\hat{\psi}_{\nu}(z)\bigg{]}\hat{\psi}_{+}(z)-\frac{\lambda}{2}\sin(2z)\big{(}\hat{a}+\hat{a}^{\dagger}\big{)}\hat{\psi}_{-}(z)$
$\displaystyle-\lambda_{\rm
ex}\sin(2z)\big{(}\hat{a}+\hat{a}^{\dagger}\big{)}\hat{\psi}_{+}(z),$ (2c)

where we have omitted the time argument for brevity.
In accordance with the damping of the membrane mode, we have introduced the
corresponding bosonic noise operator $\hat{\xi}_{\rm m}$ that is characterized
by the zero mean fluctuations $\langle{\hat{\xi}_{\rm
m}}\rangle=\langle{\hat{\xi}_{\rm m}^{\dagger}}\rangle=0$, as well as by the
auto-correlations

$\displaystyle\langle{\hat{\xi}_{\rm m}(t)\hat{\xi}_{\rm
m}^{\dagger}(0)}\rangle$ $\displaystyle=2\Gamma_{\rm m}(n_{\rm
m}+1)\delta(t),$ (3a) $\displaystyle\langle{\hat{\xi}_{\rm
m}^{\dagger}(t)\hat{\xi}_{\rm m}(0)}\rangle$ $\displaystyle=2\Gamma_{\rm
m}n_{\rm m}\delta(t).$ (3b)

The environment occupation number $n_{\rm m}$ determines the steady-state
occupation of the membrane mode when it is affected only by its environment
[57].

We assume that the atoms are at ultra-low temperature and that the atom-
membrane coupling is such that a large fraction of atoms condenses. Thus, the
system dynamics can be described by the following coupled mean-field equations
of motion

$\displaystyle i\partial_{t}\alpha$ $\displaystyle=\big{(}\Omega_{\rm
m}-i\Gamma_{\rm m}\big{)}\alpha-\int
dz\>\\!\sin(2z)\bigg{[}\frac{\sqrt{N}\lambda}{2}\big{(}\psi_{+}^{*}(z)\psi_{-}(z)+\psi_{-}^{*}(z)\psi_{+}(z)\big{)}+\sqrt{N}\lambda_{\rm
ex}|\psi_{+}(z)|^{2}\bigg{]},$ (4a) $\displaystyle i\partial_{t}\psi_{-}(z)$
$\displaystyle=\bigg{[}-\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{-}(z)+N\sum_{\nu=\pm}g_{-\nu}|\psi_{\nu}(z)|^{2}\bigg{]}\psi_{-}(z)-\frac{\sqrt{N}\lambda}{2}\sin(2z)\big{(}\alpha+\alpha^{*}\big{)}\psi_{+}(z),$
(4b) $\displaystyle i\partial_{t}\psi_{+}(z)$
$\displaystyle=\bigg{[}\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{+}(z)+N\sum_{\nu=\pm}g_{+\nu}|\psi_{\nu}(z)|^{2}\bigg{]}\psi_{+}(z)-\frac{\sqrt{N}\lambda}{2}\sin(2z)\big{(}\alpha+\alpha^{*}\big{)}\psi_{-}(z)$
$\displaystyle-\sqrt{N}\lambda_{\rm
ex}\sin(2z)\big{(}\alpha+\alpha^{*}\big{)}\psi_{+}(z),$ (4c)

where we have neglected all mutual correlations and introduced the complex
membrane amplitude $\alpha=\langle{\hat{a}}\rangle/\sqrt{N}$, and the wave
functions of the two condensates
$\psi_{\pm}(z)=\langle{\hat{\psi}_{\pm}(z)}\rangle/\sqrt{N}$. We assume that
the membrane damping is so fast that the membrane almost instantaneously
relaxes to its steady state determined by the condensate wave functions

$\displaystyle\bar{\alpha}$ $\displaystyle=\frac{1}{\Omega_{\rm
m}-i\Gamma_{\rm m}}\int dz\>\\!\sin(2z)\bigg{[}\sqrt{N}\lambda_{\rm
ex}|\psi_{+}(z)|^{2}+\frac{\sqrt{N}\lambda}{2}\big{(}\psi_{+}^{*}(z)\psi_{-}(z)+\psi_{-}^{*}(z)\psi_{+}(z)\big{)}\bigg{]},$
(5)

so that we obtain the following coupled equations for the two condensates

$\displaystyle i\partial_{t}\psi_{-}(z)$
$\displaystyle=\bigg{[}-\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{-}(z)+N\sum_{\nu=\pm}g_{-\nu}|\psi_{\nu}(z)|^{2}\bigg{]}\psi_{-}(z)-\frac{\sqrt{N}\lambda}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\psi_{+}(z),$
(6a) $\displaystyle i\partial_{t}\psi_{+}(z)$
$\displaystyle=\bigg{[}\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{+}(z)+N\sum_{\nu=\pm}g_{+\nu}|\psi_{\nu}(z)|^{2}\bigg{]}\psi_{+}(z)-\frac{\sqrt{N}\lambda}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\psi_{-}(z)$
$\displaystyle-\sqrt{N}\lambda_{\rm
ex}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\psi_{+}(z).$ (6b)

Due to the damping of the membrane, the atoms also eventually reach the
minimal-energy steady state, i.e., the ground steady state (GSS), described by
the ansatz $\psi_{\pm}(z,t)=e^{-i\mu t}\psi_{\pm}(z)$. It involves the
chemical potential $\mu$ that is fixed by the normalization condition

$\displaystyle\int
dz\;\\!\big{(}|\psi_{-}(z)|^{2}+|\psi_{+}(z)|^{2}\big{)}=1.$ (7)

It follows that the GSS solutions of the equations of motion (6) can be
obtained either by solving the stationary equations

$\displaystyle 0$ $\displaystyle=\bigg{[}-\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}-\mu+V_{-}(z)+N\sum_{\nu=\pm}g_{-\nu}|\psi_{\nu}(z)|^{2}\bigg{]}\psi_{-}(z)-\frac{\sqrt{N}\lambda}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\psi_{+}(z),$
(8a) $\displaystyle 0$ $\displaystyle=\bigg{[}\frac{\Omega_{\rm
a}}{2}-\omega_{\rm
R}\partial_{z}^{2}-\mu+V_{+}(z)+N\sum_{\nu=\pm}g_{+\nu}|\psi_{\nu}(z)|^{2}\bigg{]}\psi_{+}(z)-\frac{\sqrt{N}\lambda}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\psi_{-}(z)$
$\displaystyle-\sqrt{N}\lambda_{\rm
ex}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\psi_{+}(z),$ (8b)

or, equivalently, by minimizing the energy functional

$\displaystyle\mathcal{E}_{N}[\psi_{\pm}]$ $\displaystyle=N\sum_{\nu=\pm}\int
dz\;\\!\psi_{\nu}^{*}(z)\bigg{[}\nu\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{\nu}(z)+\frac{N}{2}\sum_{\nu^{\prime}=\pm}g_{\nu\nu^{\prime}}|\psi_{\nu^{\prime}}(z)|^{2}\bigg{]}\psi_{\nu}(z)$
$\displaystyle-\frac{N^{2}\Omega_{\rm m}}{\Omega_{\rm m}^{2}+\Gamma_{\rm
m}^{2}}\bigg{\\{}\int
dz\>\\!\sin(2z)\bigg{[}\frac{\lambda}{2}\big{(}\psi_{+}^{*}(z)\psi_{-}(z)+\psi_{-}^{*}(z)\psi_{+}(z)\big{)}+\lambda_{\rm
ex}|\psi_{+}(z)|^{2}\bigg{]}\bigg{\\}}^{2},$ (9)

with respect to $\psi_{\pm}^{*}(z)$ and under the normalization constraint
(7). This minimization justifies our notion of the GSS. The chemical potential
stems from the relation $\mu=\partial\mathcal{E}_{N}[\psi_{\pm}]/\partial N$
and reads

$\displaystyle\mu$ $\displaystyle=\sum_{\nu=\pm}\int
dz\;\\!\psi_{\nu}^{*}(z)\bigg{[}\nu\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+V_{\nu}(z)+N\sum_{\nu^{\prime}=\pm}g_{\nu\nu^{\prime}}|\psi_{\nu^{\prime}}(z)|^{2}\bigg{]}\psi_{\nu}(z)$
$\displaystyle-\frac{2N\Omega_{\rm m}}{\Omega_{\rm m}^{2}+\Gamma_{\rm
m}^{2}}\bigg{\\{}\int
dz\>\\!\sin(2z)\bigg{[}\frac{\lambda}{2}\big{(}\psi_{+}^{*}(z)\psi_{-}(z)+\psi_{-}^{*}(z)\psi_{+}(z)\big{)}+\lambda_{\rm
ex}|\psi_{+}(z)|^{2}\bigg{]}\bigg{\\}}^{2}.$ (10)

Let us next examine the ground steady state of the system. The external
trapping potentials are designed to cancel the lattice potentials due to the
driving laser and result in the box-like traps $V_{\pm}(z)$ with the hard
walls at $z\in\\{z_{1},z_{2}\\}$, where $z_{1}=-\pi/4-\ell\pi$,
$z_{2}=-\pi/4+\ell\pi$ for some large integer $\ell\gtrsim 100$. That imposes
the Dirichlet boundary conditions $\psi_{\pm}(z_{1})=\psi_{\pm}(z_{2})=0$.
According to our assumption $\ell\gg 1$, we can safely neglect any edge
effects that significantly affect the condensates only within a few healing
lengths of the walls. Since the atom-membrane interaction is periodic with the
period of $\pi$, we can reduce our attention onto the single period
$[-\pi/4,3\pi/4]$. We impose the periodic boundary conditions
$\psi_{-}(-\pi/4)=\psi_{-}(3\pi/4)$, $\psi_{+}(-\pi/4)=\psi_{+}(3\pi/4)$ and
take $V_{\pm}(z)=0$ henceforward. Note that the spatial dependence of the
atom-membrane interaction, $\propto\sin(2z)$, is even with respect to the
center of the chosen period. This brings us to the energy per period

$\displaystyle\tilde{\mathcal{E}}_{\tilde{N}}[\tilde{\psi}_{\pm}]$
$\displaystyle\equiv\frac{\mathcal{E}_{N}[\psi_{\pm}]}{2\ell}=\tilde{N}\sum_{\nu=\pm}\int_{\pi}dz\;\\!\tilde{\psi}_{\nu}^{*}(z)\bigg{[}\nu\frac{\Omega_{\rm
a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+\frac{\tilde{N}}{2}\sum_{\nu^{\prime}=\pm}g_{\nu\nu^{\prime}}|\tilde{\psi}_{\nu^{\prime}}(z)|^{2}\bigg{]}\tilde{\psi}_{\nu}(z)$
$\displaystyle-\frac{\tilde{N}^{2}\Omega_{\rm m}}{\Omega_{\rm
m}^{2}+\Gamma_{\rm
m}^{2}}\bigg{\\{}\int_{\pi}dz\>\\!\sin(2z)\bigg{[}\frac{\tilde{\lambda}}{2}\big{(}\tilde{\psi}_{+}^{*}(z)\tilde{\psi}_{-}(z)+\tilde{\psi}_{-}^{*}(z)\tilde{\psi}_{+}(z)\big{)}+\tilde{\lambda}_{\rm
ex}|\tilde{\psi}_{+}(z)|^{2}\bigg{]}\bigg{\\}}^{2},$ (11)

where we introduced the rescaled quantities

$\displaystyle
N=2\ell\tilde{N},\quad\psi_{\pm}(z)=\tilde{\psi}_{\pm}(z)/\sqrt{2\ell},\quad\lambda=\tilde{\lambda}/\sqrt{2\ell},\quad\lambda_{\rm
ex}=\tilde{\lambda}_{\rm ex}/\sqrt{2\ell}.$ (12)

The integrals in Eq. (II) run over one single period of length $\pi$. In terms
of the above, the membrane amplitude becomes

$\displaystyle\bar{\alpha}$ $\displaystyle=\frac{1}{\Omega_{\rm
m}-i\Gamma_{\rm
m}}\int_{\pi}dz\>\\!\sin(2z)\bigg{[}\sqrt{\tilde{N}}\tilde{\lambda}_{\rm
ex}|\tilde{\psi}_{+}(z)|^{2}+\frac{\sqrt{\tilde{N}}\tilde{\lambda}}{2}\big{(}\tilde{\psi}_{+}^{*}(z)\tilde{\psi}_{-}(z)+\tilde{\psi}_{-}^{*}(z)\tilde{\psi}_{+}(z)\big{)}\bigg{]},$
(13)

while the rescaled wave functions evolve in accordance with

$\displaystyle i\partial_{t}\tilde{\psi}_{-}(z)$
$\displaystyle=\bigg{[}-\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+\tilde{N}\sum_{\nu=\pm}g_{-\nu}|\tilde{\psi}_{\nu}(z)|^{2}\bigg{]}\tilde{\psi}_{-}(z)-\frac{\sqrt{\tilde{N}}\tilde{\lambda}}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\tilde{\psi}_{+}(z),$
(14a) $\displaystyle i\partial_{t}\tilde{\psi}_{+}(z)$
$\displaystyle=\bigg{[}\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+\tilde{N}\sum_{\nu=\pm}g_{+\nu}|\tilde{\psi}_{\nu}(z)|^{2}\bigg{]}\tilde{\psi}_{+}(z)-\frac{\sqrt{\tilde{N}}\tilde{\lambda}}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\tilde{\psi}_{-}(z)$
$\displaystyle-\sqrt{\tilde{N}}\tilde{\lambda}_{\rm
ex}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\tilde{\psi}_{+}(z),$
(14b)

and are normalized according to

$\displaystyle\int_{\pi}dz\;\\!\big{(}|\tilde{\psi}_{-}(z)|^{2}+|\tilde{\psi}_{+}(z)|^{2}\big{)}=1.$
(15)

The stationary wave function equations, which determine the GSS, turn into

$\displaystyle 0$ $\displaystyle=\bigg{[}-\frac{\Omega_{\rm a}}{2}-\omega_{\rm
R}\partial_{z}^{2}-\mu+\tilde{N}\sum_{\nu=\pm}g_{-\nu}|\tilde{\psi}_{\nu}(z)|^{2}\bigg{]}\tilde{\psi}_{-}(z)-\frac{\sqrt{\tilde{N}}\tilde{\lambda}}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\tilde{\psi}_{+}(z),$
(16a) $\displaystyle 0$ $\displaystyle=\bigg{[}\frac{\Omega_{\rm
a}}{2}-\omega_{\rm
R}\partial_{z}^{2}-\mu+\tilde{N}\sum_{\nu=\pm}g_{+\nu}|\tilde{\psi}_{\nu}(z)|^{2}\bigg{]}\tilde{\psi}_{+}(z)-\frac{\sqrt{\tilde{N}}\tilde{\lambda}}{2}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\tilde{\psi}_{-}(z)$
$\displaystyle-\sqrt{\tilde{N}}\tilde{\lambda}_{\rm
ex}\sin(2z)\big{(}\bar{\alpha}+\bar{\alpha}^{*}\big{)}\tilde{\psi}_{+}(z),$
(16b)

where the chemical potential is given by

$\displaystyle\mu$
$\displaystyle=\sum_{\nu=\pm}\int_{\pi}dz\;\\!\tilde{\psi}_{\nu}^{*}(z)\bigg{[}\nu\frac{\Omega_{\rm
a}}{2}-\omega_{\rm
R}\partial_{z}^{2}+\tilde{N}\sum_{\nu^{\prime}=\pm}g_{\nu\nu^{\prime}}|\tilde{\psi}_{\nu^{\prime}}(z)|^{2}\bigg{]}\tilde{\psi}_{\nu}(z)$
$\displaystyle-\frac{2\tilde{N}\Omega_{\rm m}}{\Omega_{\rm m}^{2}+\Gamma_{\rm
m}^{2}}\bigg{\\{}\int_{\pi}dz\>\\!\sin(2z)\bigg{[}\frac{\tilde{\lambda}}{2}\big{(}\tilde{\psi}_{+}^{*}(z)\tilde{\psi}_{-}(z)+\tilde{\psi}_{-}^{*}(z)\tilde{\psi}_{+}(z)\big{)}+\tilde{\lambda}_{\rm
ex}|\tilde{\psi}_{+}(z)|^{2}\bigg{]}\bigg{\\}}^{2}.$ (17)

In the following, we will assume that the rescalings (12) have been performed,
and, for the sake of brevity, we will omit the tildes. All extensive
quantities addressed below are therefore per period.

##  III Homogeneous ground steady states

We first analyze the homogeneous solutions using the ansatz
$\psi_{\pm}(z)=\sqrt{n_{\pm}/\pi}$, where $n_{\pm}$ are the occupation
fractions of the two atomic states. In such a case, the normalization
condition simply reads $n_{-}+n_{+}=1$ and the membrane is not excited, i.e.,
$\bar{\alpha}=0$. The stationary equations (16) become

$\displaystyle 0$ $\displaystyle=\bigg{[}-\frac{\Omega_{\rm
a}}{2}-\mu+\frac{N}{\pi}\sum_{\nu=\pm}g_{-\nu}n_{\nu}\bigg{]}\sqrt{n_{-}^{\vphantom{*}}}\>\\!,$
(18a) $\displaystyle 0$ $\displaystyle=\bigg{[}\frac{\Omega_{\rm
a}}{2}-\mu+\frac{N}{\pi}\sum_{\nu=\pm}g_{+\nu}n_{\nu}\bigg{]}\sqrt{n_{+}^{\vphantom{*}}}\>\\!.$
(18b)

We will be considering the on average repulsive regime where
$g_{--},g_{++}>0>g_{+-}$ and $g_{--}g_{++}>g_{+-}^{2}$. Three solutions are
possible, which read

$\displaystyle n_{-}$ $\displaystyle=$ $\displaystyle 1,\quad
n_{+}=0,\quad\mu=-\frac{\Omega_{\rm a}}{2}+\frac{N}{\pi}g_{--},$ (19a)
$\displaystyle n_{-}$ $\displaystyle=$ $\displaystyle 0,\quad
n_{+}=1,\quad\mu=\frac{\Omega_{\rm a}}{2}+\frac{N}{\pi}g_{++},$ (19b) as well
as $\displaystyle n_{-}={}$ $\displaystyle\frac{g_{++}-g_{+-}+\pi\Omega_{\rm
a}/N}{g_{--}+g_{++}-2g_{+-}},$ $\displaystyle n_{+}={}$
$\displaystyle\frac{g_{--}-g_{+-}-\pi\Omega_{\rm
a}/N}{g_{--}+g_{++}-2g_{+-}},$ $\displaystyle\mu={}$
$\displaystyle\frac{\Omega_{\rm
a}}{2}\frac{g_{--}-g_{++}}{g_{--}+g_{++}-2g_{+-}}$
$\displaystyle+\frac{N}{\pi}\frac{g_{--}g_{++}-g_{+-}^{2}}{g_{--}+g_{++}-2g_{+-}}.$
(19c)

The first solution is the lowest in energy when $\pi\Omega_{\rm
a}/N>g_{--}-g_{+-}$ holds, while for $\pi\Omega_{\rm a}/N<g_{+-}-g_{++}$, the
second one takes over. The third solution exists only if
$g_{+-}-g_{++}<\pi\Omega_{\rm a}/N<g_{--}-g_{+-}$ is satisfied and has the
lowest energy of the three. For $\pi\Omega_{\rm a}/N<(g_{--}-g_{++})/2$ we
have $n_{-}<n_{+}$, while for $\pi\Omega_{\rm a}/N\geq(g_{--}-g_{++})/2$ we
get $n_{-}\geq n_{+}$. Thus, it is convenient to introduce the reference
transition frequency $\Omega_{\rm a,0}=N(g_{--}-g_{++})/(2\pi)$ and study the
behavior of the system as a function of the relative detuning

$\displaystyle\delta\Omega_{\rm a}\equiv\Omega_{\rm a}-\Omega_{\rm a,0},$ (20)

of the atomic transition frequency from it. For strong enough couplings
$\lambda$ and $\lambda_{\rm ex}$ the system transitions into a self-organized
GSS. The membrane becomes excited and the condensates become inhomogeneous.
The nature of this phase transition is what we are going to investigate next.

##  IV Phase transitions to and from inhomogeneous ground steady state

Let us now explore the self-organized phase. To this end, we evolve the
equations of motion (14) in imaginary time starting with a suitable initial
state until the ground steady state is reached. As can be seen from these
equations, the membrane couples to the condensates via the combination ${\rm
Re}[\bar{\alpha}]\sin(2z)$. The inhomogeneous GSS wave functions have to reach
their maximum together with $\sin(2z)$, i.e., in the middle of the chosen
period. Also, $\psi_{\pm}(z)$ have to be even functions with respect to the
same point, in the same manner as $\sin(2z)$. In this way, the atom-membrane
interaction term in Eq. (II) contributes most to the minimization of the
energy of the system. Thus, the combination ${\rm Re}[\bar{\alpha}]\sin(2z)$
captures the essence of the membrane’s influence on the atoms, both in terms
of its spatial behavior and its relative strength. Figure 2 shows two
exemplary cases of spatial profiles of ${\rm Re}[\bar{\alpha}]\sin(2z)$,
together with the two GSS wave functions $\psi_{\pm}(z)$ over the chosen
period. In the panel (a) we present the case of a negative relative detuning
$\delta\Omega_{\rm a}$, where the $|+\rangle$ state usually has higher
occupation. The panel (b) displays the opposite case with the $|-\rangle$
state being more populated. As can be seen, larger spatial variances of the
wave functions correspond to larger membrane amplitude.

Figure 2: Spatial profiles of ${\rm Re}[\bar{\alpha}]\sin(2z)$ (green) and the
ground steady state condensate wave functions $\psi_{-}(z)$ (red) and
$\psi_{+}(z)$ (blue) over one period for large negative (a) and large positive
(b) relative detuning. Parameters: $\Omega_{\rm m}=100\omega_{\rm R}$,
$\Gamma_{\rm m}=10\omega_{\rm R}$, $Ng_{--}=Ng_{++}=100\omega_{\rm R}$,
$Ng_{+-}=-90\omega_{\rm R}$, $\sqrt{N}\lambda_{\rm ex}=35\omega_{\rm R}$,
while in (a) $\delta\Omega_{\rm a}=-40\omega_{\rm R}$,
$\sqrt{N}\lambda=5\omega_{\rm R}$ and in (b) $\delta\Omega_{\rm
a}=70\omega_{\rm R}$, $\sqrt{N}\lambda=70\omega_{\rm R}$.

Figure 3: Dependence of the real part of the membrane amplitude ${\rm
Re}[\bar{\alpha}]$ (a) and the population imbalance
$d_{-+}\equiv(N_{-}-N_{+})/N$ (b) between the $|-\rangle$ and the $|+\rangle$
state on the relative detuning $\delta\Omega_{\rm a}$ for a fixed value of the
coupling $\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ and for three values of
the coupling $\sqrt{N}\lambda\in\\{20,70,120\\}\;\\!\omega_{\rm R}$ (green,
red, blue curves). Other parameters: $\Omega_{\rm m}=100\omega_{\rm R}$,
$\Gamma_{\rm m}=10\omega_{\rm R}$, $Ng_{--}=Ng_{++}=100\omega_{\rm R}$,
$Ng_{+-}=-90\omega_{\rm R}$.

For large negative relative detunings, the condensate is homogeneous and the
atoms occupy exclusively the $|+\rangle$ state. Likewise, for large positive
$\delta\Omega_{\rm a}$ the atoms homogeneously condense into the $|-\rangle$
state. In between, both atomic states exhibit a macroscopic population, while
the system may transition from the homogeneous to the self-organized phase.
Such transitions occur for sufficiently large couplings to the membrane
$\lambda$ and $\lambda_{\rm ex}$. Their values determine whether the
transitions are continuous or discontinuous.

In Fig. 3 we showcase the behavior of ${\rm Re}[\bar{\alpha}]$ and the
population imbalance $d_{-+}\equiv(N_{-}-N_{+})/N$ when the relative detuning
is scanned across zero, for a fixed value of the coupling
$\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ and for three values of the
coupling $\sqrt{N}\lambda\in\\{20,70,120\\}\;\\!\omega_{\rm R}$ (green, red,
blue curves). Here, we defined $N_{\pm}$ as the population per period of the
state $|\pm\rangle$. The green curves depict a case of a continuous second-
order phase transition on the negative side, followed by a continuous second-
order transition on the positive side as well. This feature is manifest in
both ${\rm Re}[\bar{\alpha}]$ and $d_{-+}$. For a larger value of $\lambda$
(red curves), the transition on the negative side becomes discontinuous and of
first order, while on the positive side it is still continuous. For even
larger value of $\lambda$ (blue curves), the latter transition also becomes a
discontinuous first-order one.

Figure 4: Spatial profiles of the condensate wave functions $\psi_{+}(z)$
(left panels) and $\psi_{-}(z)$ (right panels) depending on the relative
detuning $\delta\Omega_{\rm a}$ for a fixed value of the coupling
$\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ and for three values of the
coupling $\sqrt{N}\lambda\in\\{20,70,120\\}\;\\!\omega_{\rm R}$ (top, middle,
bottom row). The remaining parameters are $\Omega_{\rm m}=100\omega_{\rm R}$,
$\Gamma_{\rm m}=10\omega_{\rm R}$, $Ng_{--}=Ng_{++}=100\omega_{\rm R}$,
$Ng_{+-}=-90\omega_{\rm R}$.

In order to acquire additional insight into the properties of the system, in
Fig. 4 we show the behavior of the condensate spatial profiles $\psi_{\pm}(z)$
along the three curves shown in Fig. 3 as we sweep the detuning from negative
to positive values. One can visually distinguish the homogeneous and the self-
organized phases. The wave functions $\psi_{\pm}(z)$ are symmetric and maximal
in the center of the period. The sharpness of the color gradients around the
transition points reflects the (dis)continuity of the respective transition.
The width of $\psi_{+}(z)$ decreases until shortly before the second
transition, while the width of $\psi_{-}(z)$ increases monotonically. Stronger
couplings $\lambda$ and $\lambda_{\rm ex}$ typically lead to narrower wave
function profiles $\psi_{\pm}(z)$ due to the membrane-induced localization of
the atoms.

Figure 5: Dependence of the real part of the membrane amplitude ${\rm
Re}[\bar{\alpha}]$ (left panels) and the population imbalance
$d_{-+}=(N_{-}-N_{+})/N$ (right panels) on the relative detuning
$\delta\Omega_{\rm a}$ and the coupling strength $\sqrt{N}\lambda$, for
$\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ (top row) and
$\sqrt{N}\lambda_{\rm ex}=25\omega_{\rm R}$ (bottom row). The horizontal
dashed lines in the top left panel mark those cases shown in Figs. 3 and 4.
The black solid lines in the right panels correspond to $N_{-}=N_{+}$. Other
parameters are: $\Omega_{\rm m}=100\omega_{\rm R}$, $\Gamma_{\rm
m}=10\omega_{\rm R}$, $Ng_{--}=Ng_{++}=100\omega_{\rm R}$,
$Ng_{+-}=-90\omega_{\rm R}$.

We extend the study of the transitions discussed above to a range of values of
the coupling strength $\lambda$. Figure 5 displays the dependence of ${\rm
Re}[\bar{\alpha}]$ (left panels) and $d_{-+}$ (right panels) on
$\delta\Omega_{\rm a}$ and $\lambda$ for two values of the coupling
$\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ (top row) and
$\sqrt{N}\lambda_{\rm ex}=25\omega_{\rm R}$ (bottom row). The horizontal
dashed lines in the top left panel mark explicitly those cases shown in Figs.
3 and 4. The single-colored blue (red) regions in the right panels correspond
to the homogeneous phase where the atoms condense into the $|+\rangle$
($|-\rangle$) state. Careful inspection of the left panels reveals that there
are blue intermediate regions where the condensates are homogeneous and both
atomic states have significant population. These appear for
$\sqrt{N}\lambda\lesssim 20\omega_{\rm R}$ and only around $|\delta\Omega_{\rm
a}|\approx 50\omega_{\rm R}$ for the larger coupling $\lambda_{\rm ex}$, as
opposed to the entire range $|\delta\Omega_{\rm a}|\lesssim 50\omega_{\rm R}$
for the smaller coupling $\lambda_{\rm ex}$. Therefore, to discriminate
experimentally between the homogeneous and the inhomogeneous phase, one has
first to measure the membrane displacement. If it turns out to be essentially
zero, further characterization of the homogeneous phase can be done by
measuring the population imbalance $d_{-+}$.

The continuous nature of the transitions to the self-organized phase is
manifest in the gradual changes of the color-coding on both the negative and
the positive detuning side for $\sqrt{N}\lambda\lesssim 40\omega_{\rm R}$ in
the former case and for $\sqrt{N}\lambda\lesssim 110\omega_{\rm R}$ in the
latter case. In the top row, the transition on the negative side becomes sharp
and discontinuous within the strip $40\omega_{\rm
R}\lesssim\sqrt{N}\lambda\lesssim 110\omega_{\rm R}$, while the one on the
positive side is still continuous. Above that, both transitions become
discontinuous. For the smaller coupling $\lambda_{\rm ex}$ in the bottom row,
the corresponding coupling strengths are $\sqrt{N}\lambda\approx
110\omega_{\rm R}$ and $\sqrt{N}\lambda\approx 160\omega_{\rm R}$,
respectively.

Figure 6: Dependence of the energy per period, Eq. (II), on the relative
detuning $\delta\Omega_{\rm a}$ and the coupling strength $\sqrt{N}\lambda$,
for $\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ (top panel) and
$\sqrt{N}\lambda_{\rm ex}=25\omega_{\rm R}$ (bottom panel). The black curves
represent the contour lines. Used parameters: $\Omega_{\rm m}=100\omega_{\rm
R}$, $\Gamma_{\rm m}=10\omega_{\rm R}$, $Ng_{--}=Ng_{++}=100\omega_{\rm R}$,
$Ng_{+-}=-90\omega_{\rm R}$.

Figure 7: Dependence of the real part of the membrane amplitude ${\rm
Re}[\bar{\alpha}]$ on the coupling strengths $\sqrt{N}\lambda$ and
$\sqrt{N}\lambda_{\rm ex}$, for $\delta\Omega_{\rm a}=-120\omega_{\rm R}$ (top
panel) and $\delta\Omega_{\rm a}=120\omega_{\rm R}$ (bottom panel). Other
parameters are: $\Omega_{\rm m}=100\omega_{\rm R}$, $\Gamma_{\rm
m}=10\omega_{\rm R}$, $Ng_{--}=50\omega_{\rm R}$, $Ng_{++}=200\omega_{\rm R}$,
$Ng_{+-}=-90\omega_{\rm R}$.

Further understanding can be gained by investigating the energy behavior.
Figure 6 portrays the energy per period, given by Eq. (II), as a function of
the relative detuning $\delta\Omega_{\rm a}$ and the coupling strength
$\sqrt{N}\lambda$, for $\sqrt{N}\lambda_{\rm ex}=40\omega_{\rm R}$ (top panel)
and $\sqrt{N}\lambda_{\rm ex}=25\omega_{\rm R}$ (bottom panel). Noticeably,
there are no discontinuities within the energy landscape. However, the energy
features gradient jumps along the curves that correspond to discontinuous
phase transitions observed in Fig. 5. This is illustrated both by the shading
and the non-smoothness of the contour lines. The energy is maximal when both
the relative detuning and the coupling $\lambda$ are zero. As expected, the
larger either of the couplings $\lambda$ or $\lambda_{\rm ex}$ is, the lower
is the energy due to the stabilizing atom-membrane interaction. We also
observe that for the fixed couplings the energy of the homogeneous phases
decreases with increasing magnitude of $\delta\Omega_{\rm a}$, while the self-
ordered phase exhibits a non-monotonic dependence in this respect.

It is worth illustrating the asymmetric role played by the couplings $\lambda$
and $\lambda_{\rm ex}$, which is obvious from the equations of motion (14).
Figure 7 shows ${\rm Re}[\bar{\alpha}]$ as a function of the two coupling
strengths, for $\delta\Omega_{\rm a}=-120\omega_{\rm R}$ (top panel) and
$\delta\Omega_{\rm a}=120\omega_{\rm R}$ (bottom panel). For presentation
purposes, and different than in other plots, we have chosen in this context
the imbalanced contact interaction strengths $Ng_{--}=50\omega_{\rm R}$ and
$Ng_{++}=200\omega_{\rm R}$. The regions that correspond to the homogeneous
phase are not symmetric with respect to $\lambda$ and $\lambda_{\rm ex}$. For
large enough relative detuning, increasing $\lambda_{\rm ex}$ leads to a
first-order transition for small $\lambda$ and potentially to a second-order
transition for larger $\lambda$. Conversely, increasing $\lambda$ leads to a
second-order transition for small $\lambda_{\rm ex}$ and to a first-order
transition for larger $\lambda_{\rm ex}$.

##  V Discussion and outlook

The above examples demonstrate that the considered system is versatile and
offers multiple tuning knobs for changing the order of the phase transition
between the homogeneous and the self-organized density wave-like state. These
include, but are not limited to, the relative detuning $\delta\Omega_{\rm a}$
of the atomic transition frequency $\Omega_{\rm a}$ and the coupling strengths
$\lambda$ and $\lambda_{\rm ex}$. The two states can be experimentally
distinguished by measuring the membrane displacement. To further characterize
the homogeneous state, one can resort to the population imbalance $d_{-+}$. We
have not explored the corresponding role of the atomic contact interaction
strengths $g_{\nu\nu^{\prime}}$ in this work. They represent, of course,
further means of controlling the behavior of the system. The asymmetric nature
of the atom-membrane coupling can either be further enhanced or effectively
suppressed by an appropriate asymmetric choice of $g_{\nu\nu^{\prime}}$.
Depending on the path taken through the experimental parameter space, one can
encounter either first- or second-order transition hyper-surfaces. A possible
future work could involve dynamics beyond the studied nonequilibrium steady
state. In such a case, the system parameters could be dynamically controlled
along different open or closed paths. In certain scenarios, hysteresis should
appear and can be studied.

It would be interesting to investigate how the above mean-field analysis may
be amended once quantum fluctuations have been taken into account. The spatial
profiles of the two condensates critically determine which atomic normal modes
are coupled to the membrane. Those that are coupled should exhibit enhanced
occupation and may significantly alter primarily the homogeneous state.
Therefore, the order of the observed phase transitions should be reanalyzed.
Moreover, quantum fluctuations may lead to the emergence of new phases, akin
to the formation of self-bound quantum droplets in Bose-Bose mixtures and
dipolar Bose gases [58, 59]. In the considered system the membrane may have an
essential role to play. Beyond mean-field quantum effects may also reveal the
imprint of NQPTs and their order on the statistics of the outcoupled light.

We find it important to note the differences between the NQPT reported here
and the one studied in Ref. [55]. First, in the cited work it was essential to
have a lattice potential present for the atoms. This led to exclusively
density wave-like condensate steady states. Second, only the special case of
atomic contact interaction strengths $g_{\nu\nu^{\prime}}=g$ was considered.
Moreover, the phase transition studied therein would occur even in the absence
of atom-atom interactions within the BEC, i.e., for $g_{\nu\nu^{\prime}}=0$.
Obviously, our present scenario requires the overall repulsive interactions.

In conclusion, in this paper we have presented a flexible hybrid atom-
optomechanical experimental platform for the study of both continuous and
discontinuous NQPTs. Importantly, it is possible to choose between the two by
tuning readily accessible experimental parameters. The present work lays the
foundation for further exploration of possible quantum phases of ultracold
atomic gases along the lines of self-bound quantum droplets.

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) via the Collaborative Research Center SFB/TR185 (Project
No. 277625399) (A.P.) and via the Research Grant 274978739 (M.R. and M.T.). We
also acknowledge the support from the DFG Cluster of Excellence “CUI: Advanced
Imaging of Matter” $-$ EXC 2056 (Project ID 390715994).

## References

  * [1] S. Sachdev, Quantum Phase Transitions, 2nd ed. (Cambridge University Press, Cambridge, 2011).
  * [2] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena, 4th ed (Oxford University Press, Oxford, 2002).
  * [3] H. Kleinert and V. Schulte-Frohlinde, Critical Properties of $\phi^{4}$-Theories, (World Scientific, Singapore, 2001).
  * [4] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Rep. Prog. Phys. 79, 096001 (2016).
  * [5] L. M. Sieberer, M. Buchhold, J. Marino, and S. Diehl, Universality in driven open quantum matter, arXiv:2312.03073 (2023).
  * [6] A. Mitra, S. Takei, Y. B. Kim, and A. J. Millis, Nonequilibrium Quantum Criticality in Open Electronic Systems, Phys. Rev. Lett. 97, 236808 (2006).
  * [7] A. Mitra and A. J. Millis, Coulomb gas on the Keldysh contour: Anderson-Yuval-Hamann representation of the nonequilibrium two-level system, Phys. Rev. B 76, 085342 (2007).
  * [8] V. M. Bastidas, C. Emary, G. Schaller, and T. Brandes, Nonequilibrium quantum phase transitions in the Ising model, Phys. Rev. A 86, 063627 (2012).
  * [9] F. Brennecke, R. Mottl, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Real-time observation of fluctuations at the driven-dissipative Dicke phase transition, Proc. Natl. Acad. Sci. USA 110, 11763 (2013).
  * [10] J. Marino and S. Diehl, Driven Markovian Quantum Criticality, Phys. Rev. Lett. 116, 070407 (2016).
  * [11] M. F. Maghrebi and A. V. Gorshkov, Nonequilibrium many- body steady states via Keldysh formalism, Phys. Rev. B 93, 014307 (2016).
  * [12] J. Hannukainen and J. Larson, Dissipation-driven quantum phase transitions and symmetry breaking, Phys. Rev. A 98, 042113 (2018).
  * [13] R. Rota, F. Minganti, C. Ciuti, and V. Savona, Quantum Critical Regime in a Quadratically Driven Nonlinear Photonic Lattice, Phys. Rev. Lett. 122, 110405 (2019).
  * [14] T. O. Puel, S. Chesi, S. Kirchner, and P. Ribeiro, Mixed-Order Symmetry-Breaking Quantum Phase Transition Far from Equilibrium, Phys. Rev. Lett. 122, 235701 (2019).
  * [15] C. Klöckner, C. Karrasch, and D. M. Kennes, Nonequilibrium Properties of Berezinskii-Kosterlitz-Thouless Phase Transitions, Phys. Rev. Lett. 125, 147601 (2020).
  * [16] T. Behrle, T. L. Nguyen, F. Reiter, D. Baur, B. de Neeve, M. Stadler, M. Marinelli, F. Lancellotti, S. F. Yelin, and J. P. Home, Phonon Laser in the Quantum Regime, Phys. Rev. Lett. 131, 043605 (2023).
  * [17] D. Nagy, G. Szirmai and P. Domokos, Self-organization of a Bose-Einstein condensate in an optical cavity, Eur. Phys. J. D 48, 127 (2008).
  * [18] C. Maschler, I. B. Mekhov, and H. Ritsch, Ultracold atoms in optical lattices generated by quantized light fields, Eur. Phys. J. D 46, 545 (2008).
  * [19] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Dicke quantum phase transition with a superfluid gas in an optical cavity, Nature 464, 1301 (2010).
  * [20] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Cold atoms in cavity-generated dynamical optical potentials, Rev. Mod. Phys. 85, 553 (2013).
  * [21] J. Klinder, H. Keßler, M. Reza Bakhtiari, M. Thorwart, and A. Hemmerich, Observation of a Superradiant Mott Insulator in the Dicke-Hubbard Model, Phys. Rev. Lett. 115, 230403 (2015).
  * [22] M. Reza Bakhtiari, A. Hemmerich, H. Ritsch, and M. Thorwart, Nonequilibrium Phase Transition of Interacting Bosons in an Intra-Cavity Optical Lattice, Phys. Rev. Lett. 114, 123601 (2015).
  * [23] F. Mivehvar, F. Piazza, T. Donner, and H. Ritsch, Cavity QED with quantum gases: new paradigms in many-body physics, Adv. Phys. 70, 1 (2021).
  * [24] N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni, Long-range interacting quantum systems, Rev. Mod. Phys. 95, 035002 (2023).
  * [25] C. D. Mink, A. Pelster, J. Benary, H. Ott, and M. Fleischhauer, Variational truncated Wigner approximation for weakly interacting Bose fields: Dynamics of coupled condensates, SciPost Phys. 12, 051 (2022).
  * [26] J. Benary, C. Baals, E. Bernhart, J. Jiang, M. Röhrle, and H. Ott, Experimental observation of a dissipative phase transition in a multi-mode many-body quantum system, New J. Phys. 24, 103034 (2022).
  * [27] G. Dagvadorj, J. M. Fellows, S. Matyjaskiewicz, F. M. Marchetti, I. Carusotto, and M. H. Szymanska, Tuning across Universalities with a Driven Open Condensate, Phys. Rev. X 5, 041028 (2015).
  * [28] A. Zamora, L. M. Sieberer, K. Dunnett, S. Diehl, and M. H. Szymanska, Tuning across Universalities with a Driven Open Condensate, Phys. Rev. X 7, 041006 (2017).
  * [29] P. Comaron, G. Dagvadorj, A. Zamora, I. Carusotto, N. P. Proukakis, and M. H. Szymanska, Dynamical Critical Exponents in Driven-Dissipative Quantum Systems, Phys. Rev. Lett. 121, 095302 (2018).
  * [30] H. J. Hesten, R. A. Nyman, and F. Mintert, Decondensation in Nonequilibrium Photonic Condensates: When Less Is More, Phys. Rev. Lett. 120, 040601 (2018).
  * [31] J. D. Rodrigues, H. S. Dhar, B. T. Walker, J. M. Smith, R. F. Oulton, F. Mintert, and R. A. Nyman, Learning the Fuzzy Phases of Small Photonic Condensates, Phys. Rev. Lett. 126, 150602 (2021).
  * [32] F. Öztürk, T. Lappe, G. Hellmann, J. Schmitt, J. Klaers, F. Vewinger, J. Kroha, and M. Weitz, Observation of a non-Hermitian phase transition in an optical quantum gas, Science 372, 6537 (2021).
  * [33] M. Heyl, Dynamical quantum phase transitions, Rep. Prog. Phys. 81, 054001 (2018).
  * [34] L.-N. Wu, J. Nettersheim, J. Feß, A. Schnell, S. Burgardt, S. Hiebel, D. Adam, A. Eckardt, and A. Widera, Indication of critical scaling in time during the relaxation of an open quantum system, arXiv:2208.05164 (2023).
  * [35] B. Vogell, K. Stannigel, P. Zoller, K. Hammerer, M. T. Rakher, M. Korppi, A. Jöckel, and P. Treutlein, Cavity-enhanced long-distance coupling of an atomic ensemble to a micromechanical membrane, Phys. Rev. A 87, 023816 (2013).
  * [36] B. Vogell, T. Kampschulte, M. T. Rakher, A. Faber, P. Treutlein, K. Hammerer, and P. Zoller, Long distance coupling of a quantum mechanical oscillator to the internal states of an atomic ensemble, New J. Phys. 17, 043044 (2015).
  * [37] A. Jöckel, A. Faber, T. Kampschulte, M. Korppi, M. T. Rakher, and P. Treutlein, Sympathetic cooling of a membrane oscillator in a hybrid mechanical–atomic system, Nat. Nanotechnol. 10, 55 (2015).
  * [38] H. Zhong, G. Fläschner, A. Schwarz, R. Wiesendanger, P. Christoph, T. Wagner, A. Bick, C. Staarmann, B. Abeln, K. Sengstock, and C. Becker A millikelvin all-fiber cavity optomechanical apparatus for merging with ultra-cold atoms in a hybrid quantum system, Rev. Sci. Instrum. 88, 023115 (2017).
  * [39] P. Christoph, T. Wagner, H. Zhong, R. Wiesendanger, K. Sengstock, A. Schwarz, and C. Becker, Combined feedback and sympathetic cooling of a mechanical oscillator coupled to ultracold atoms, New J. Phys. 20, 093020 (2018).
  * [40] A. Vochezer, T. Kampschulte, K. Hammerer, and P. Treutlein, Light-Mediated Collective Atomic Motion in an Optical Lattice Coupled to a Membrane, Phys. Rev. Lett. 120, 073602 (2018).
  * [41] N. Mann, M. Reza Bakhtiari, A. Pelster, and M. Thorwart, Nonequilibrium Quantum Phase Transition in a Hybrid Atom-Optomechanical System, Phys. Rev. Lett. 120, 063605 (2018).
  * [42] C. Gao and Z. Liang, Steady-state phase diagram of quantum gases in a lattice coupled to a membrane, Phys. Rev. A 99, 013629 (2019).
  * [43] T. M. Karg, B. Gouraud, C. T. Ngai, G.-L. Schmid, K. Hammerer, and P. Treutlein, Light-mediated strong coupling between a mechanical oscillator and atomic spins 1 meter apart, Science 369, 174 (2020).
  * [44] G.-L. Schmid, C. T. Ngai, M. Ernzer, M. B. Aguilera, T. M. Karg, and P. Treutlein, Coherent Feedback Cooling of a Nanomechanical Membrane with Atomic Spins, Phys. Rev. X 12, 011020 (2022).
  * [45] M. Ernzer, M. B. Aguilera, M. Brunelli, G.-L. Schmid, T. M. Karg, C. Bruder, P. P. Potts, and P. Treutlein, Optical Coherent Feedback Control of a Mechanical Oscillator, Phys. Rev. X 13, 021023 (2023).
  * [46] K. Hammerer, M. Aspelmeyer, E. S. Polzik, and P. Zoller, Establishing Einstein-Poldosky-Rosen Channels between Nanomechanics and Atomic Ensembles, Phys. Rev. Lett. 102, 020501 (2009).
  * [47] K. Hammerer, M. Wallquist, C. Genes, M. Ludwig, F. Marquardt, P. Treutlein, P. Zoller, J. Ye, and H. J. Kimble, Strong Coupling of a Mechanical Oscillator and a Single Atom, Phys. Rev. Lett. 103, 063005 (2009).
  * [48] M. Wallquist, K. Hammerer, P. Rabl, M. Lukin, and P. Zoller, Hybrid quantum devices and quantum engineering, Phys. Scr. T137, 014001 (2009).
  * [49] M. Paternostro, G. De Chiara, and G. M. Palma, Cold-Atom-Induced Control of an Optomechanical Device, Phys. Rev. Lett. 104, 243602 (2010).
  * [50] C. Genes, H. Ritsch, M. Drewsen, and A. Dantan, Atom-membrane cooling and entanglement using cavity electromagnetically induced transparency, Phys. Rev. A 84, 051801(R) (2011).
  * [51] H. K. Lau, A. Eisfeld, and J. M. Rost, Cavity-free quantum optomechanical cooling by atom-modulated radiation, Phys. Rev. A 98, 043827 (2018).
  * [52] C. F. Ockeloen, R. Schmied, M. F. Riedel, and P. Treutlein, Quantum Metrology with a Scanning Probe Atom Interferometer, Phys. Rev. Let. 111, 143001 (2013).
  * [53] E. S. Polzik and K. Hammerer, Trajectories without quantum uncertainties, Ann. Phys. (Leipzig) 527, A15 (2015).
  * [54] C. B. Møller, R. A. Thomas, G. Vasilakis, E. Zeuthen, Y. Tsaturyan, M. Balabas, K. Jensen, A. Schliesser, K. Hammerer, and E. S. Polzik, Quantum back-action-evading measurement of motion in a negative mass reference frame, Nature 547, 191 (2017).
  * [55] N. Mann, A. Pelster, and M. Thorwart, Tuning the order of the nonequilibrium quantum phase transition in a hybrid atom-optomechanical system, New J. Phys. 21, 113037 (2019).
  * [56] J.-M. Cheng, Z.-W. Zhou, G.-C. Guo, H. Pu, and X.-F. Zhou, Bose-Einstein condensates in an atom-optomechanical system with effective global nonuniform interaction, Phys. Rev. A 103, 023328 (2021).
  * [57] C. W. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics, Springer Series in Synergetics (Springer, Berlin, 2004).
  * [58] Z.-H. Luo, W. Pang, B. Liu, Y.-Y. Li, B. A. Malomed, A new form of liquid matter: Quantum droplets, Front. Phys. 16, 32201 (2021).
  * [59] F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, New states of matter with fine-tuned interactions: quantum droplets and dipolar supersolids, Rep. Prog. Phys. 84, 012403 (2021).