# Fly-by galaxy encounters with multiple black holes produce star-forming
linear wakes

Nianyi Chen  McWilliams Center for Cosmology, Department of Physics, Carnegie
Mellon University, Pittsburgh, PA 15213  Patrick LaChance  McWilliams Center
for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh,
PA 15213  Yueying Ni  Center for Astrophysics $|$ Harvard & Smithsonian,
Cambridge, MA 02138, US  Tiziana Di Matteo  McWilliams Center for Cosmology,
Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213
Rupert Croft  McWilliams Center for Cosmology, Department of Physics, Carnegie
Mellon University, Pittsburgh, PA 15213  Priyamvada Natarajan  Black Hole
Initiative, Harvard University, Cambridge, MA 02138, USA  Department of
Astronomy, Yale University, New Haven, CT 06511, USA  Department of Physics,
Yale University, New Haven, CT 06520, USA  Simeon Bird  Department of Physics
& Astronomy, University of California, Riverside, 900 University Ave.,
Riverside, CA 92521, US Nianyi Chen<EMAIL_ADDRESS>
###### Abstract

We look for simulated star-forming linear wakes such as the one recently
discovered by van Dokkum et al. (2023) in the cosmological hydrodynamical
simulation ASTRID. Amongst the runaway black holes in ASTRID, none are able to
produce clear star-forming wakes. Meanwhile, fly-by encounters, typically
involving a compact galaxy (with a central black hole) and a star-forming
galaxy (with a duo of black holes) reproduce remarkably well many of the key
properties (its length and linearity; recent star formation, etc.) of the
observed star-forming linear feature. We predict the feature to persist for
approximately 100 Myr in such a system and hence constitute a rare event. The
feature contains a partly stripped galaxy (with $M_{\rm gal}=10^{9}\sim
10^{10}M_{\odot}$) and a dual BH system ($M_{\rm BH}=10^{5}\sim
10^{7}\,M_{\odot}$) in its brightest knot. X-ray emission from AGN in the knot
should be detectable in such systems. After $100\sim 200\,{\rm Myrs}$ from the
first fly-by, the galaxies merge leaving behind a triple black hole system in
an (still) actively star-forming early-type remnant of mass $\sim 5\times
10^{10}\,M_{\odot}$. Follow-up JWST observations may be key for revealing the
nature of these linear features by potentially detecting the older stellar
populations constituting the bright knot. Confirmation of such detections may
therefore help discriminate a fly-by encounter from a massive BH wake to
reveal the origin of such features.

Computational methods — Astronomy data analysis — Supermassive Black holes

##  1 Introduction

Given that nearly all galaxies host central black holes, it is likely that
during complex dynamical mergers involving multiple galaxies, one or more of
the BHs get ejected. Evidence from observations for such anticipated three-
body interactions has been sparse till recently. New observations have
reported serendipitous discoveries of extended linear star-forming features
with unclear origins that likely involve multiple BH interactions and
accompanying triggered star formation. van Dokkum et al. (2023) (VD23
hereafter) reports the discovery of a thin stellar streak extending $\sim
60\,{\rm kpc}$ from a nearby compact galaxy at $z\sim 1$. Based on the length
and linearity of the feature, as well as the emission lines associated with
dense star-forming gas, the authors interpret it as the stellar wake induced
by the passage of a massive black hole (MBH) kicked out from the compact
galaxy during a 3-body encounter of the MBHs. In another recent paper,
Zaritsky et al. (2023) reports finding an even longer, collimated galactic
tail $380\,{\rm kpc}$ in length, also associated with recent and ongoing star
formation. Having discussed the challenges of a more commonplace origin for
this feature such as ram pressure or tidal stripping, these authors also
propose that this feature is likely created by a three-body encounter between
three interacting galaxies. Other observed cases of potential runaway BHs
include the one detected in a peculiar source in the COSMOS survey at
$z=0.359$ (CXOCJ100043.1+020637 (CID-42)) that presents as two compact optical
sources embedded in the same galaxy with a large velocity offset (Civano et
al., 2010).

In this work, we focus on linear features with accompanying star formation
activity. If the observed linear features were indeed produced by a runaway BH
in a compact star cluster, or even a runaway galaxy itself, they open up a new
channel for finding ejected stars and MBHs undergoing complex dynamical
interaction. However, the runaway BH explanation itself also faces theoretical
challenges: it is not clear, for example, whether a BH with $M_{\rm BH}\sim
10^{7}\,M_{\odot}$ can trigger such a high level of star formation, and leave
a trail of stars weighing almost half the mass of the original galaxy
($>10^{9}\,M_{\odot}$ as per VD23). On the other hand, it is also
theoretically challenging to produce bright stellar wakes that are so thin and
straight, without having some ejected material associated with them. Further
clarification on the origin of these exotic systems demands work from the
theoretical side, and the first step would be to identify theoretical
counterparts of these observed systems in cosmological simulations, a path we
pursue here, so that we can better investigate their origins.

To this end, our work presents a first attempt to query and investigate the
origin of extended linear stellar wakes in cosmological simulations and their
possible association with runaway MBHs. We use the large volume ASTRID
simulation with a large galaxy population (Bird et al., 2022; Ni et al.,
2022), allowing us to find the very rare occurrences of extended star-forming
systems. Moreover, the updated BH dynamics model (e.g. Chen et al., 2022)
produces a vast number of wandering MBHs (Di Matteo et al., 2022) and even
potential runaway MBHs among MBH triplets (Hoffman et al., 2023), allowing us
to look for connections between stellar wakes and MBH tuples or runaway BHs.

Our letter is organized as follows: in Section 2 we first introduce the ASTRID
simulation, and the selection criteria we used to find runaway BHs in the
simulation as well as the theoretical counterparts of the currently observed
linear features. In Section 3, we present systems from the simulation that
morphologically resemble the observed linear star-forming features. We then
examine their origin, evolution, and detailed properties of the galaxies and
BHs associated with such systems.

##  2 Method

Figure 1: Target galaxies for the search of runaway BHs and linear stellar
features (orange points), plotted on top of the distribution of all galaxies
with $M_{\rm gal}>5\times 10^{8}\,M_{\odot}$ in ASTRID at $z=2$. The two
systems with the most prominent linear star-forming features are shown as
stars ($z\sim 2$ system in blue, $z\sim 1$ system in red). We focus our search
among the star-forming (SFR$>1\,M_{\odot}/{\rm yr}$), compact ($r_{\rm
half}<2{\rm kpc}$) galaxies with masses $5\times 10^{9}\,\,M_{\odot}<M_{\rm
gal}<5\times 10^{10}\,\,M_{\odot}$ that contains at least 3 BHs with a central
$M_{\rm BH}>10^{7}\,\,M_{\odot}$ BH, as indicated by the grey dashed lines.

###  2.1 Simulations

ASTRID is a cosmological hydrodynamical simulation performed using a new
version of the Smoothed Particle Hydrodynamics code MP-Gadget. The simulation
evolves a cube of $250h^{-1}{\rm Mpc}$ per side with $2\times 5500^{3}$
initial tracer particles comprising dark matter and baryons, and has currently
reached $z=1.3$. ASTRID has a dark matter particle mass resolution of $M_{\rm
DM}=6.7\times 10^{6}{h^{-1}M_{\odot}}$ and $M_{\rm gas}=1.3\times
10^{6}{h^{-1}M_{\odot}}$ in the initial conditions. The gravitational
softening length is $\epsilon_{\rm g}=1.5h^{-1}\,{\rm kpc}$ for both DM and
gas particles. The simulation includes a full-physics sub-grid treatment for
modeling galaxy formation, SMBHs and their associated supernova and AGN
(Active Galactic Nuclei) feedback, as well as inhomogeneous hydrogen and
helium reionization. We refer the readers to the introductory papers Bird et
al. (2022); Ni et al. (2022) for detailed descriptions of physical models
deployed in the simulation.

Here we briefly summarize BH modeling, which is most relevant to our current
investigation of runaway MBHs. BHs are seeded in haloes with $M_{\rm
halo,FOF}>5\times 10^{9}{h^{-1}M_{\odot}}$ and $M_{\rm*,FOF}>2\times
10^{6}{h^{-1}M_{\odot}}$, with seed masses stochastically drawn between
$3\times 10^{4}{h^{-1}M_{\odot}}$ and $3\times 10^{5}{h^{-1}M_{\odot}}$,
motivated by the direct collapse scenario proposed in Lodato & Natarajan
(2007). The gas accretion rate onto the black hole is estimated via a Bondi-
Hoyle-Lyttleton-like prescription (Di Matteo et al., 2005). The black hole
radiates with a bolometric luminosity $L_{\rm bol}$ proportional to the
accretion rate $\dot{M}_{\bullet}$, with a mass-to-energy conversion
efficiency $\eta=0.1$ in an accretion disk according to Shakura & Sunyaev
(1973). 5% of the radiated energy is coupled to the surrounding gas as the AGN
feedback. Dynamics of the black holes are modeled with a sub-grid dynamical
friction model (Tremmel et al., 2015; Chen et al., 2022) to replace the
original implementation that directly repositioned BHs to the local minimum of
the potential. This gives well-defined black hole trajectories and velocities.
Per this implementation, two black holes can merge if their separation is
within two times the gravitational softening length $2\epsilon_{g}$, once
their kinetic energy is dissipated by dynamical friction and they are
gravitationally bound.

###  2.2 Galaxy Population and Candidate Selection

Galaxies in the simulation are identified with the SUBFIND algorithm (Springel
et al., 2001). The large volume of ASTRID provides a rich galaxy population
matching the host galaxy of the observed systems, from which we can conduct a
comprehensive search for runaway BHs and star-forming wakes. We perform the
search at two redshifts, $z=2$ and $z=1.3$ (the current lowest redshift of the
simulation), with the latter being closer in time to the observed streak in
VD23 at $z=0.964$.

Target galaxies: In Figure 1, we show the properties of $z=2$ galaxies in
ASTRID and our selection criteria. To match with the properties of the main
galaxy in van Dokkum et al. (2023), we focus on the star-forming
(SFR$>1\,M_{\odot}/{\rm yr}$), compact ($r_{\rm half}<2{\rm kpc}$), and non-
satellite galaxies with masses $5\times 10^{9}\,\,M_{\odot}<M_{\rm
gal}<5\times 10^{10}\,\,M_{\odot}$ that contain at least 3 BHs with a total BH
mass above $10^{7}\,\,M_{\odot}$. We apply the 3-BH criterion because it is
necessary to produce a runaway BH. After applying these cuts, we narrow down
our search to $29888$ (11399) galaxies at $z=1.3$ ($z=2$), and hereafter we
refer to these selected galaxies as our “target galaxies”.

Figure 2: An example runaway-BH involved in an MBH triplet (marked by red,
green, and blue crosses) identified in ASTRID at $z\sim 1.7$. The red cross is
the BH that got ejected from the host galaxy with a velocity of $\sim
700\,{\rm km/s}$, due to the complex dynamics within the halo. The host galaxy
mass is $M_{\rm gal}=1.3\times 10^{11}\,M_{\odot}$, and the runaway BH has a
mass of $M_{\rm BH}=1.3\times 10^{7}\,M_{\odot}$, and the central BH (blue
cross) has $M_{\rm BH}=6.6\times 10^{8}\,M_{\odot}$.

Runaway BH selection: We identify runaway BHs within the target galaxies above
(i.e. the orange population in Figure 1). To quantitatively define a “runaway”
MBH in these galaxies, we compute the ratio between consecutive apocentric
radii between each satellite BH in the galaxy and the central MBH (similar to
Hoffman et al., 2023). If we see a consistent and sudden increase in orbital
size by a factor of more than three in at least three full orbits, then we
categorize the MBH as a “runaway” MBH. Among the $\sim 11399$ target galaxies
at $z=2$, we found $\sim 200$ potential runaway BHs that satisfy these
criteria.

Linear star-forming wake selection: We look for radial linear star-forming
features with the following method: we divide all the young stars (stellar age
$<100\,{\rm Myrs}$) in each halo into 45 radial bins [10, 100] kpc from the
central galaxy center. To match a $\sim 50\,{\rm kpc}$ observable linear
feature, we select systems with at least $2\times 10^{6}\,M_{\odot}$ stars in
each bin for 25 continuous radial bins. This gives us 514 (435) candidates
with potentially observable stellar streaks at $z=1.3$ ($z=2$) respectively,
which makes up $2-4\%$ of the target galaxy population (we will refer to these
as “linear feature candidates”). We then visually inspect these potential
linear features to select the ones with large signal-to-noise, resulting in
$\sim 30$ visually-confirmed linear features at each redshift (note that this
number of viable candidates is likely a lower limit as we only look for
linearity along the $x$-$y$ and $x$-$z$ planes).

##  3 Results

###  3.1 Do runaway BHs produce linear stellar wakes?

van Dokkum et al. (2023) proposed that one of the most likely mechanisms for
producing the observed star-forming wake is a runaway BH induced by a three-
body encounter. Following this proposal and the works by Saslaw & De Young
(1972) and de la Fuente Marcos & de la Fuente Marcos (2008), we look for
possible associations between linear star-forming features and runaway BHs.
However, upon an initial search through the target galaxies (orange points in
Figure 1), we do not find any visible stellar features along the wakes of the
runaway BHs.

As the BH-induced star formation may be more prominent within a denser circum-
galactic medium (CGM) than that in our target galaxies, we relax the upper
limit on the galaxy mass and did another search among MBH triplets in larger
halos (the sample of MBH triplets from Hoffman et al. (2023)). In these halos,
the runaway BHs also tend to be more massive, and the formation of a stellar
wake is less affected by the resolution limit. Figure 2 shows an example
runaway BH candidate within a MBH triplet, found in a galaxy with mass
$1.3\times 10^{11}\,M_{\odot}$. In this system, a BH with mass $1.3\times
10^{7}\,M_{\odot}$ is ejected from the host galaxy $\sim 100\,{\rm Myrs}$
after the initial formation of the MBH triplet at a speed of $\sim 700\,{\rm
km/s}$, and it traverses a distance of $\sim 45\,{\rm kpc}$ from the galaxy
center in another $\sim 100\,{\rm Myrs}$. Even in this case, we still do not
see a star-forming signature associated with the wake of the runaway BH: we
found at most $<10^{6}M_{\odot}$ total newly formed stars within $10\,{\rm
kpc}$ of the ejected BH.

We note a few limitations of the simulation potentially responsible for this
null result: the resolution of ASTRID cannot fully resolve the ejection due to
three-body encounters, and BHs that “run away” are mostly disrupted by a blob
of incoming stars/gas. This implies that the maximum speed reached by any of
our runaway BH candidates is $\sim 800\,{\rm km/s}$ and decreases to $\sim
100\,{\rm km/s}$ at the target separation of $\sim 50\,{\rm kpc}$, barely
reaching the speed to shock the surrounding gas (e.g. Dopita & Sutherland,
1996; Allen et al., 2008). Also, the simulation may lack the mass resolution
to resolve gas shocked by a BH: the target runaway BH ($M_{\rm BH}\sim
10^{7}M_{\odot}$) is a few times larger than the gas particle mass, and so the
gas over-density drawn by the BH may not be directly resolvable. A follow-up
study would need to deploy higher-resolution simulations, to detect BHs
ejected at a higher velocity in environments similar to our target galaxies to
evaluate whether a similar star-forming wake can be produced.

###  3.2 Stellar wakes associated with fly-by galaxy encounters

Figure 3: Row 1/2: The stars (first row, color-coded by age) and gas (second
row, color-coded by temperature) associated with the linear stellar streaks at
$z\sim 2$ (left) and $z\sim 1$ (right). The linear features extend to $\sim
50\,{\rm kpc}$ from the compact galaxy on the left of each image. Row 3: Mock
image of the systems in the HST F606W and F814W filters. Row 4: F606W - F814W
color of the stellar population, modeled with with the Python-FSPS stellar
population model (Conroy & Gunn, 2010), and the mean age of the $5\%$ youngest
stars. Row 5: Mock image of the systems in the JWST filters will reveal the
older stellar population along the wake.

Although we do not find visible star-forming wakes associated with the passage
of runaway BHs, we do find a few tens of star-forming wakes in the simulation
following the method detailed in Section 2. The natural question to ask is
then: what mechanisms, if not runaway BHs, produce these wakes? By inspecting
the simulation counterparts to the observed systems selected in Section 2.2,
we find that most star-forming streaks originate from the fly-by encounter of
a recently merged, dual-BH young galaxy with a more massive compact galaxy. We
present some representative cases of such systems and their evolution in this
section.

Figure 3 shows star and gas visualizations of two systems with linear star-
forming features, and the mock images of the systems seen through the
different filters of HST and JWST. In both cases, the young stars extend
linearly to more than 50 kpc away from the main, compact galaxy (Galaxy 1, on
the left end of each image). The young galaxy (Galaxy 2, on the right) is
embedded in a stream of cold, dense gas with $T\sim 10^{4}-10^{5}\,{\rm K}$
clearly distinguishable from the background hot gas in the surrounding medium
(second row), along which we see ongoing star formation and young stars formed
within a few tens Myrs (fourth row).

Figure 3 also shows the mock HST and JWST observations, created by assigning
spectral energy distributions (SEDs) to each star particle according to the
Binary Population and Spectral Population Synthesis model (BPASS v2.2.1;
Stanway & Eldridge, 2018). These SEDs are convolved with the filter
transmission functions associated with the chosen filters. We smooth the stars
with an SPH kernel and make 2D projections with pixel sizes matching the
sensor for HST and JWST (0.049” for HST ACS, and 0.031”/0.063” for JWST NIRCAM
short/long wavelength).

The stellar age along the wake follows a bi-modal distribution: half of the
stars are relatively old with ages $\sim 1\,{\rm Gyr}$, and half are formed
within $<100\,{\rm Myrs}$ during the emergence of the streak. The young stars
form the linear feature when captured by the HST F606W and F814W bands (third
row), while the older stars are better seen with JWST at longer wavelengths
(fifth row). The redder stars may be a distinctive signature of the stripped
galaxy scenario, which are unlikely to be found in the runaway-BH case.

###  3.3 Evolution of the linear feature

Figure 4:  Time evolution of the two galaxies involved in producing the linear
stellar feature (projected face-on), and the associated BHs (central BH in the
Galaxy 1 in blue, two BHs in the tidally disrupted galaxy in green and
magenta, and the another BH nearby in red). The green line traces the
trajectory of BH2 relative to BH1, and the red line shows the orbit of BH3
relative to BH2.

Figure 5: Properties of the three BHs along the linear stellar feature. The
top two panels show the X-ray luminosities and BH masses of the three BHs. The
bottom panel shows the distance between the three BHs.

To better understand the formation of the star-forming streak, we trace System
1 across a Gyr before and after the time when the linear feature is most
prominent. Figure 4 shows the time evolution of the two galaxies involved in
the linear stellar feature in a face-on view. The feature is a by-product of a
fly-by encounter between two galaxies at a speed of $580\,{\rm km/s}$. The
galaxies barely touch each other during their first passage (frame 1), after
which Galaxy 1 remains unperturbed while Galaxy 2 develops elongated tidal
arms on both sides. The feature formed around $150\,{\rm Myrs}$ after the
initial passage between the two galaxies, and lasts for $\sim 200\,{\rm Myrs}$
(frame 2-4), after which the stars along the streak become old (frame 5), and
the two galaxies experience a head-on collision (frame 6-7).

Our simulated linear features are extreme cases of the tidal arms/tails found
in previous numerical studies of galaxy mergers (e.g. Toomre & Toomre, 1972;
Barnes & Hernquist, 1992; Ploeckinger et al., 2018). Idealized simulations
have shown that within these tails accumulations of gas and stars can be found
(e.g. Barnes & Hernquist, 1992; Baumschlager et al., 2019; Shin et al., 2020).
These clumps harbor molecular gas which provides a reservoir for new stars
(e.g. Mirabel et al., 1991, 1992), and could be the origin of the blue knots
along the VD23 system. Upon the close inspection of the $\sim 50$ visually-
confirmed linear features, we find that the long tails are usually associated
with dual-BH galaxies. $70\%$ of these tails contain close BH pairs,
indicating a recent galaxy merger within the tail. Isolated galaxy simulations
may not produce features like System 1 and System 2 if it takes two
consecutive galaxy mergers to produce very long star-forming tails.

We also show in Figure 4 the trajectories of the BHs along the tail. BH3 has
come within Galaxy 2 in a previous galaxy merger a few hundred Myrs before the
wake formation. Then it keeps orbiting around BH2 during and after the
production of the stellar wake, until the final merger between Galaxy 1 and
Galaxy 2 (frame 7), after which both BH2 and BH3 ended up in wide orbits
around BH1. We will examine the detailed evolution of the triple-BH system in
the next section.

###  3.4 BHs along the linear feature

A common feature shared by our linear feature candidates is a dual-BH system
along the streak merging into a central BH. Such MBH triplets are indicators
of consecutive galaxy mergers producing the stellar streak, and if active,
they may be also responsible for triggering the OIII emission along the
streak. Here we examine the evolution of BHs along the linear features and
their AGN activity.

Figure 5 shows the time evolution of the X-ray luminosity, BH mass, and the
orbits of the three BHs associated with the merging galaxies in System 1 and
System 2. In both cases, the BH in the compact galaxy is fairly luminous,
maintaining an X-ray luminosity of nearly $10^{43}\,{\rm erg/s}$ for more than
a Gyr. The second brightest BH (BH2) embedded along the linear feature also
has a relatively high X-ray emission between $10^{42}-10^{43}\,{\rm erg/s}$.
If the observed wake is produced by the scenario we show here, then BH2 is
responsible for triggering the OIII emission along the wake. As mentioned in
VD23, the amount of OIII emission seen would correspond to an AGN with
$L_{X}>10^{43}\,{\rm erg/s}$ according to the scaling relation in Ueda et al.
(2015). Our BH2s are slightly fainter for a $1.9\times 10^{41}\,{\rm erg/s}$
OIII luminosity, but still fall well within the Ueda et al. (2015) relation
with the scatter and are capable of triggering the observed OIII signature.

The detectability of the third BH may depend on its orbits within Galaxy 2,
but it is possible that BH3 was captured at the pericenter of the orbit and at
it is brightest when the linear feature is produced (e.g. the case in System
1). In this (optimistic) case, follow-up observations may simultaneously see
an MBH triplet associated with the stellar wake.

##  4 Conclusions and Discussion

In this work, we identify runaway BHs ejected at $\sim 700\,{\rm km/s}$ from
host galaxies, as well as linear star-forming features extending $>50\,{\rm
kpc}$ from compact galaxies in the ASTRID simulation at $z=1\sim 2$. We look
for the association between the two types of systems following the proposal in
VD23, and find that runaway BHs do not produce visible star-forming wakes in
the simulation.

However, we find that linear star-forming features in fly-by galaxy encounters
resemble the recent observation by van Dokkum et al. (2023), with very bright,
thin, and straight young stars when seen in the HST F606W/F814W filters. Such
linear features potentially exist among $\sim 0.1\%$ compact, star-forming
galaxies with $5\times 10^{9}\,\,M_{\odot}<M_{\rm gal}<5\times
10^{10}\,\,M_{\odot}$ and multiple BHs, and the most prominent linear features
can last for $100-200\,{\rm Myrs}$.

We examine two representative systems with strong linear features in detail,
both with stellar masses between $10^{10}-5\times 10^{10}\,M_{\odot}$ and SFRs
between $10-50M_{\odot}/{\rm yr}$. Among the BHs associated with the two
galaxies producing the linear feature, the central BH has a bright X-ray
emission with $10^{42}\,{\rm erg/s}<L_{X}<10^{43}\,{\rm erg/s}$, which could
potentially be seen by Chandra. The brightest BH along the feature emits at
$10^{41}\,{\rm erg/s}<L_{X}<10^{43}\,{\rm erg/s}$, also with a chance to be
detected in the X-rays.

The production of the linear feature involves consecutive encounters between
three galaxies. The remnant of the first galaxy merger ($\sim
10^{9}-10^{10}\,M_{\odot}$ and star-forming) goes through a fly-by encounter
with a more massive, compact galaxy, during which its surrounding gas gets
tidally stripped, leaving a trace of star-forming gas and young stars. $\sim
70\%$ of the visually-confirmed systems with linear features harbor a dual-BH
star-forming galaxy along the feature, supporting the three-galaxy encounter
scenario. Our findings indicate that these linear features may offer a robust
signature for finding BH binaries and multiple-BH systems, but further work is
required to better establish the connection between linear long tidal tails
and binary-hosting galaxies. There have been catalogs of long tidal tails in
the COSMOS field (e.g. Civano et al., 2010; Wen & Zheng, 2016; Ren et al.,
2020), out of which two galaxy nuclei can be clearly identified in some
extended linear systems.

Finally, although our results show that only tidal features during fly-by
galaxy mergers can produce these recently observed features, we acknowledge
that our resolution limit leaves room for the possibility of a runaway-BH-
induced star-forming wake. Due to the similarity with many features of the
observed signatures (star-forming gas, OIII emission, and a population of
young stars) shared by these two scenarios, follow-up studies are needed to
potentially distinguish between them.

Follow-up observations using longer-wavelength bands with JWST can tell
whether there exists an underlying old stellar population along the linear
feature. The lack of old stars will greatly support the runaway BH origin of
the feature. Meanwhile, higher-resolution simulations can study the triggering
of such dense star formation by runaway BHs, and the mass/velocity of the BH
required. In both formation channels of the linear feature, follow-up
observations may be able to reveal multiple MBH signatures along the feature,
and such systems are very likely the sites for complex three-body encounters
and mergers between MBHs (e.g. Hoffman & Loeb, 2007; Bonetti et al., 2016;
Mannerkoski et al., 2021).

## Acknowledgements

We thank Mohit Bhardwaj, Pieter van Dokkum, Charlie Conroy and Qian Yang for
their helpful discussions. ASTRID was run on the Frontera facility at the
Texas Advanced Computing Center. TDM and RACC acknowledge funding from the NSF
AI Institute: Physics of the Future, NSF PHY-2020295, NASA ATP NNX17AK56G, and
NASA ATP 80NSSC18K101. TDM acknowledges additional support from NSF
ACI-1614853, NSF AST-1616168, NASA ATP 19-ATP19-0084, and NASA ATP
80NSSC20K0519, and RACC from NSF AST-1909193. SB acknowledges funding
supported by NASA-80NSSC22K1897.

## Data Availability

The code to reproduce the simulation is available at https://github.com/MP-
Gadget/MP-Gadget, and continues to be developed. Part of the ASTRID snapshots
are available at https://astrid-portal.psc.edu/.

## References

  * Allen et al. (2008) Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20, doi: 10.1086/589652
  * Barnes & Hernquist (1992) Barnes, J. E., & Hernquist, L. 1992, Nature, 360, 715, doi: 10.1038/360715a0
  * Baumschlager et al. (2019) Baumschlager, B., Hensler, G., Steyrleithner, P., & Recchi, S. 2019, MNRAS, 483, 5315, doi: 10.1093/mnras/sty3418
  * Bird et al. (2022) Bird, S., Ni, Y., Di Matteo, T., et al. 2022, MNRAS, 512, 3703, doi: 10.1093/mnras/stac648
  * Bonetti et al. (2016) Bonetti, M., Haardt, F., Sesana, A., & Barausse, E. 2016, MNRAS, 461, 4419, doi: 10.1093/mnras/stw1590
  * Chen et al. (2022) Chen, N., Ni, Y., Tremmel, M., et al. 2022, MNRAS, 510, 531, doi: 10.1093/mnras/stab3411
  * Civano et al. (2010) Civano, F., Elvis, M., Lanzuisi, G., et al. 2010, ApJ, 717, 209, doi: 10.1088/0004-637X/717/1/209
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, FSPS: Flexible Stellar Population Synthesis, Astrophysics Source Code Library, record ascl:1010.043.  http://ascl.net/1010.043
  * de la Fuente Marcos & de la Fuente Marcos (2008) de la Fuente Marcos, R., & de la Fuente Marcos, C. 2008, ApJ, 677, L47, doi: 10.1086/587962
  * Di Matteo et al. (2022) Di Matteo, T., Ni, Y., Chen, N., et al. 2022, arXiv e-prints, arXiv:2210.14960, doi: 10.48550/arXiv.2210.14960
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604, doi: 10.1038/nature03335
  * Dopita & Sutherland (1996) Dopita, M. A., & Sutherland, R. S. 1996, ApJS, 102, 161, doi: 10.1086/192255
  * Hoffman et al. (2023) Hoffman, C., Chen, N., Di Matteo, T., et al. 2023, arXiv e-prints, arXiv:2303.04825, doi: 10.48550/arXiv.2303.04825
  * Hoffman & Loeb (2007) Hoffman, L., & Loeb, A. 2007, MNRAS, 377, 957, doi: 10.1111/j.1365-2966.2007.11694.x
  * Lodato & Natarajan (2007) Lodato, G., & Natarajan, P. 2007, MNRAS, 377, L64, doi: 10.1111/j.1745-3933.2007.00304.x
  * Mannerkoski et al. (2021) Mannerkoski, M., Johansson, P. H., Rantala, A., Naab, T., & Liao, S. 2021, ApJ, 912, L20, doi: 10.3847/2041-8213/abf9a5
  * Mirabel et al. (1992) Mirabel, I. F., Dottori, H., & Lutz, D. 1992, A&A, 256, L19
  * Mirabel et al. (1991) Mirabel, I. F., Lutz, D., & Maza, J. 1991, A&A, 243, 367
  * Ni et al. (2022) Ni, Y., Di Matteo, T., Bird, S., et al. 2022, MNRAS, 513, 670, doi: 10.1093/mnras/stac351
  * Ploeckinger et al. (2018) Ploeckinger, S., Sharma, K., Schaye, J., et al. 2018, MNRAS, 474, 580, doi: 10.1093/mnras/stx2787
  * Ren et al. (2020) Ren, J., Zheng, X. Z., Valls-Gabaud, D., et al. 2020, MNRAS, 499, 3399, doi: 10.1093/mnras/staa2985
  * Saslaw & De Young (1972) Saslaw, W. C., & De Young, D. S. 1972, Astrophys. Lett., 11, 87
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shin et al. (2020) Shin, E.-j., Jung, M., Kwon, G., et al. 2020, ApJ, 899, 25, doi: 10.3847/1538-4357/aba434
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726, doi: 10.1046/j.1365-8711.2001.04912.x
  * Stanway & Eldridge (2018) Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75, doi: 10.1093/mnras/sty1353
  * Toomre & Toomre (1972) Toomre, A., & Toomre, J. 1972, ApJ, 178, 623, doi: 10.1086/151823
  * Tremmel et al. (2015) Tremmel, M., Governato, F., Volonteri, M., & Quinn, T. R. 2015, MNRAS, 451, 1868, doi: 10.1093/mnras/stv1060
  * Ueda et al. (2015) Ueda, Y., Hashimoto, Y., Ichikawa, K., et al. 2015, ApJ, 815, 1, doi: 10.1088/0004-637X/815/1/1
  * van Dokkum et al. (2023) van Dokkum, P., Pasha, I., Buzzo, M. L., et al. 2023, ApJ, 946, L50, doi: 10.3847/2041-8213/acba86
  * Wen & Zheng (2016) Wen, Z. Z., & Zheng, X. Z. 2016, ApJ, 832, 90, doi: 10.3847/0004-637X/832/1/90
  * Zaritsky et al. (2023) Zaritsky, D., Crossett, J. P., Jaffé, Y. L., et al. 2023, arXiv e-prints, arXiv:2305.01335, doi: 10.48550/arXiv.2305.01335


# Quantized Quasi-Two Dimensional Bose-Einstein Condensates with Spatially
Modulated Nonlinearity

Deng-Shan Wang1, Xing-Hua Hu1, Jiangping Hu2 and W. M. Liu1 1Beijing National
Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy
of Sciences, Beijing $100190$, P.R. China  2Department of Physics, Purdue
University, West Lafayette, Indiana $47907$, U.S.A.

###### Abstract

We investigate the localized nonlinear matter waves of the quasi-two
dimensional Bose-Einstein condensates with spatially modulated nonlinearity in
harmonic potential. It is shown that the whole Bose-Einstein condensates,
similar to the linear harmonic oscillator, can have an arbitrary number of
localized nonlinear matter waves with discrete energies, which are
mathematically exact orthogonal solutions of the Gross-Pitaevskii equation.
Their novel properties are determined by the principle quantum number $n$ and
secondary quantum number $l$: the parity of the matter wave functions and the
corresponding energy levels depend only on $n$, and the numbers of density
packets for each quantum state depend on both $n$ and $l$ which describe the
topological properties of the atom packets. We also give an experimental
protocol to observe these novel phenomena in future experiments.

###### pacs:

03.75.Hh, 05.45.Yv, 67.85.Bc

_Introduction_.—Since the remarkable experimental realization Anderson ;
Hulet1 ; Davis  of Bose-Einstein condensations (BEC), there has been an
explosion of the experimental and theoretical activity devoted to the physics
of dilute ultracold bosonic gases. It is known that the properties of BEC
including their shape, collective nonlinear excitations are determined by the
sign and magnitude of the $s$-wave scattering length. A prominent way to
adjust scattering length is to tune an external magnetic field in the vicinity
of a Feshbach resonance Inouye1 . Alternatively, one can use a Feshbach
resonance induced by optical or electric field Theis . Since all quantities of
interest in the BEC crucially depend on scattering length, a tunable
interaction suggests very interesting studies of the many-body behavior of
condensate systems.

In the past years, techniques for adjusting the scattering length globally
have been crucial to many experimental achievements Herbig ; Bartenstein .
More recently, condensates with a spatially modulated nonlinearity by
manipulating scattering length locally have been proposed Rodas-Verde ;
Sakaguchi ; Konotop ; Qian ; Belmonte-Beitia12 . This is experimentally
feasible due to the flexible and precise control of the scattering length with
tunable interactions. The spatial dependence of scattering length can be
implemented by a spatially inhomogeneous external magnetic field in the
vicinity of a Feshbach resonance Xiong .

However, so far, the studies of BEC with spatially modulated nonlinearity are
limited in the quasi-one dimensional cases Rodas-Verde ; Sakaguchi ; Konotop ;
Qian ; Belmonte-Beitia12 . Moreover, in the study of nonlinear problems no one
discusses their quantum properties which are common in linear systems such as
the linear harmonic oscillator. In this Letter, we extend the similarity
transformation Belmonte-Beitia12  to the quasi-two dimensional (quasi-2D) BEC
with spatially modulated nonlinearity in harmonic potential, and find a family
of stable localized nonlinear matter wave solutions. Similar to the linear
harmonic oscillator, we discover that the whole BEC can be quantized which is
unexpected before. Their quantum and topological properties can be simply
described by two quantum numbers. We also formulate an experimental procedure
for the realization of these novel phenomena in 7Li condensate Hulet1 ; Hulet2
. This opens the door to the investigation of new matter waves in the high
dimensional BEC with spatially modulated nonlinearities.

Model and exact localized solutions.—The system considered here is a BEC
confined in a harmonic trap
$V(\textbf{r})=m(\omega_{\perp}^{2}r^{2}+\omega_{z}^{2}z^{2})/2$, where $m$ is
atomic mass, $r^{2}=x^{2}+y^{2}$, and $\omega_{\perp},\omega_{z}$ are the
confinement frequencies in the radial and axial directions, respectively. In
the mean-field theory, the BEC system at low temperature is described by the
Gross-Pitaevskii (GP) equation in three dimensions. If the trap is pancake-
shaped, i.e. $\omega_{z}\gg\omega_{\perp},$ it is reasonable to reduce the GP
equation for the condensate wave function to a quasi-2D equation Kivshar ;
Ueda1 ; Garcia-Ripoll

$i\psi_{t}=-\frac{1}{2}(\psi_{xx}+\psi_{yy})+\frac{1}{2}\omega^{2}(x^{2}+y^{2})\psi+g(x,y)|\psi|^{2}\psi,$
(1)

where $\omega=\omega_{\perp}/\omega_{z}$, the length, time and wave function
$\psi$ are measured in units of
$a_{h}=\sqrt{\hbar/m\omega_{z}},\omega_{z}^{-1},a_{h}^{-1}$ and $g(x,y)=4\pi
a_{s}(x,y)$ represents the strength of interatomic interaction characterized
by the $s$-wave scattering length $a_{s}(x,y)$, which can be spatially
inhomogeneous by magnetically tuning the Feshbach resonances Inouye1 ; Rodas-
Verde ; Sakaguchi ; Konotop ; Qian ; Belmonte-Beitia12 ; Xiong .

Figure 1: (color online). The interaction parameter $g(x,y)$ for two secondary
quantum numbers: $(a)~{}l=0$ and $(b)~{}l=1$ with $\omega=0.02,\nu=0.1$. It is
seen that $g(x,y)$ is a smooth function when $l=0$ and develops singularity
when $l$ gets large.

Now we consider the spatially localized stationary solution
$\psi(x,y,t)=\phi(x,y)e^{-i\mu t}$ of Eq. (1) with $\phi(x,y)$ being a real
function for $\lim_{|x|,|y|\rightarrow\infty}\phi(x,y)=0.$ This maps Eq. (1)
onto a stationary nonlinear Schrödinger equation
$\frac{1}{2}\phi_{xx}+\frac{1}{2}\phi_{yy}-\frac{1}{2}\omega^{2}(x^{2}+y^{2})\phi-g(x,y)\phi^{3}+\mu\phi=0$
Pethick . Here $\mu$ is the real chemical potential. Solving this stationary
equation by similarity transformation Belmonte-Beitia12 , we obtain a families
of exact localized nonlinear wave solutions for Eq. (1) as

$\psi_{n}={\frac{(n+1)K(k)\eta}{\sqrt{\nu}}}\,{\rm cn}(\theta,k)e^{-i\mu
t},n=0,2,4,\cdots$ (2) $\psi_{n}={\frac{(n+1)K(k)\eta}{\sqrt{2\nu}}}\,{\rm
sd}(\theta,k)e^{-i\mu t},n=1,3,5,\cdots$ (3)

where $k=\sqrt{2}/{2}$ is the modulus of elliptic function, $\nu$ is a
positive real constant,
$K(k)=\int_{0}^{\frac{\pi}{2}}[1-k^{2}\sin^{2}\varsigma]d\varsigma$ is
elliptic integral of the first kind, ${\rm sd}={\rm sn}/{\rm dn}$ with ${\rm
sn},{\rm cn}$ and ${\rm dn}$ being Jacobi elliptic functions, $\theta,\eta$
and $g$ are determined by

$\theta=(n+1)K(k){\rm erf}[\sqrt{2\omega}\left(x+y\right)/2],$
$\eta={e^{\omega\,xy}}{\rm
KummerU}[-\mu/(2\omega),1/2,\omega\,\left(x-y\right)^{2}/2],$ (4)
$g(x,y)=-2\omega\,\nu/(\pi\eta^{2}){e^{-\omega\,\left(x+y\right)^{2}}},$

here ${\rm erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-\tau^{2}}d\tau$ is
error function, and ${\rm KummerU}(a,c,s)$ Abramowitz  is Kummer function of
the second kind which is a solution of ordinary differential equation
$s\Lambda^{{}^{\prime\prime}}(s)+(c-s)\Lambda^{{}^{\prime}}(s)-a\Lambda(s)=0.$
It is easy to see that when $|x|,|y|\rightarrow\infty$ we have
$\psi_{n}\rightarrow 0$ for solutions $\psi_{n}$ in Eqs. (2)-(3) with Eq. (4),
thus they are localized bound state solutions.

In the above construction, it is observed that the number of zero points of
function $\eta$ in Eq. (4) is equal to that of function ${\rm
KummerU}[-\mu/(2\omega),1/2,\omega\,\left(x-y\right)^{2}/2],$ which strongly
depends on $\omega$ and the ratio $\mu/\omega.$ We assume the number of zero
points in $\eta$ along line $y=-x$ is $l.$ In the following, we will see that
integer $n$ is associated with the energy levels of the atoms and integers
$n,l$ determine the topological properties of atom packets, so $n$ and $l$ are
named the principal quantum number and secondary quantum number in quantum
mechanics. In addition, the three free parameters $\omega$, $\mu$ and $\nu$
are positive, so the dimensionless interaction function $g(x,y)$ is negative,
which indicates an attractive interaction between atoms. There are known
atomic gases with attractive interactions realized by modulating magnetic
Inouye1  technique, for examples, the 85Rb Cornish  and 7Li atoms Hulet1 ;
Hulet2 .

Figure 2: (color online). The density distributions of the quasi-2D BEC with
spatially modulated nonlinearities in harmonic potential, for different
principle quantum numbers $n$ in Eqs. (2)-(3) with Eq. (4), where the
secondary quantum number $l=0,$ the parameters $\omega,\nu$ are $0.02$ and
$0.1,$ respectively. The unit of space length $x,y$ is $1.69~{}\mu m.$ Figs.
2(a)-2(c) show the density profiles of the even parity wave function (2) for
$n=0,2$ and $4,$ respectively. Figs. 2(d)-2(f) demonstrate the density
profiles of the odd parity wave function (3) for $n=1,3$ and $5,$
respectively.

To translate our results into units relevant to the experiments Hulet1 ;
Hulet2 , we take the 7Li condensate containing $10^{3}\sim 10^{5}$ atoms in a
pancake-shaped trap with radial frequency $\omega_{\perp}=2\pi\times 10$ Hz
and axial frequency $\omega_{z}=2\pi\times 500$ Hz Rychtarik . In this case,
the ratio of trap frequency $\omega$ in Eq. (1) is $0.02$ which is determined
by $\omega_{\perp}/\omega_{z}.$ The unit of length is $1.69~{}\mu m$, the unit
of time is $0.32~{}ms$ and the unit of chemical potential is $nK$. The
spatially inhomogeneous interaction parameter $g(x,y)$ is independent of
principal quantum number $n$ but is strongly related to the secondary quantum
number $l$. In the Fig. 1, we show that for $\omega=0.02,\nu=0.1,$ function
$g(x,y)$ is smooth in space when $l=0$ and develops singularity when the $l$
gets large.

Quantized quasi-$2$D BEC.—In order to investigate the quantum and topological
properties of the localized nonlinear matter waves in quasi-2D BEC described
by Eqs. (2)-(3) with Eq. (4), we plot their density distributions by
manipulating the principal quantum number $n$ or secondary quantum number $l$.

Figure 3: (color online). The density distributions of the quasi-2D BEC in
harmonic potential for different secondary quantum number $l$. Figs. 3(a)-3(d)
show the density profiles of the even parity wave function (2) for principle
quantum number $n=0,$ and Figs. 3(e)-3(h) show the density profiles of the odd
parity wave function (3) for $n=1,$ corresponding to $l=0,1,2,3$. The other
parameters are the same as that of Fig. 2.

Firstly, when the secondary quantum number $l$ is fixed, we can modulate the
principal quantum number $n$ to analyze the novel matter waves in quasi-2D
BEC. Fig. 2 shows the density profiles in quasi-2D BEC with spatially
modulated nonlinearities in harmonic potential for $l=0$. It is easy to see
that the matter wave functions in Eq. (2) satisfy
$\psi_{n}(-x,-y)=\psi_{n}(x,y),$ so they are even parity and are invariant
under space inversion. Figs. 2(a)-2(c) demonstrate the density profiles of the
even parity wave functions (2) with Eq. (4) for $n=0,2,4,$ which correspond to
a low energy state and two highly excited states. The numbers of atoms
$N_{n}=\int\int dxdy|\psi_{n}(x,y,t)|^{2}$ for the three states are
$N_{0}=3.76\times 10^{3},N_{2}=6.84\times 10^{4},N_{4}=2.633\times 10^{5}$,
respectively. The matter wave functions in Eq. (3) satisfy
$\psi_{n}(-x,-y)=-\psi_{n}(x,y),$ which denotes that they are odd parity.
Figs. 2(d)-2(f) demonstrate the density profiles of the odd parity wave
functions (3) with Eq. (4) for $n=1,3,5,$ which correspond to three highly
excited states. The numbers of atoms for the three states are
$N_{1}=4.016\times 10^{4},N_{3}=2.493\times 10^{5},N_{5}=7.28\times 10^{5},$
respectively. It is observed that when the secondary quantum number $l=0$, the
number of nodes along line $y=x$ for each quantum state is equal to the
corresponding principal quantum number $n$, i.e. the $n$th level quantum state
has $n$ nodes along $y=x$. And the number of density packets increases one by
one along line $y=x$ when the $n$ increases. This is similar to the quantum
properties in the linear harmonic oscillator.

Secondly, when the principal quantum number $n$ is fixed, we can tune the
secondary quantum number $l$ to observe the novel quantum phenomenon in
quasi-2D BEC. In Fig. 3 we demonstrate the density distributions of quasi-2D
BEC in harmonic potential for different secondary quantum number. Figs.
3(a)-3(d) show the density profiles of the even parity wave function (2) with
Eq. (4) for $n=0,$ and $l=0,1,2$ and $3,$ respectively. It is seen that the
number of nodes for the density packets along line $y=-x$ is equal to the
corresponding secondary quantum number $l$ which describes the topological
patterns of the atom packets, and the number of density packets increases one
by one when $l$ increases. Figs. 3(e)-3(h) show the density profiles of the
odd parity wave function (3) with Eq. (4) for $n=1$ and $l=0,1,2,3.$ We see
that the number of density packets increases pair by pair when $l$ increases.
The number of density packets for each quantum state is equal to
$(n+1)\times(l+1),$ and all the density packets are symmetrical with respect
to lines $y=\pm x,$ as shown in Figs. 2-3.

Normalization energy vs chemical potential. Next we calculate the
normalization energy of each quantum states numerically. The total energy of
the quasi-2D BEC is $E(\psi)=\int\int
dxdy[|\nabla\psi|^{2}+\frac{1}{2}\omega^{2}(x^{2}+y^{2})|\psi|^{2}+\frac{1}{2}g(x,y)|\psi|^{4}]$.
So the normalized energy is given by $E(\psi)/N=\mu-\frac{1}{2N}\int\int
dxdyg(x,y)|\psi|^{4}$ with $N=\int\int dxdy|\psi|^{2}$. Fig. 4 shows the
relations of the normalization energy $E(\psi)/N$ with chemistry potential for
different principle quantum numbers $n$. It is observed that for the fixed
$n$, the normalization energy is approximatively linear increase with respect
to chemistry potential, i.e., $d(E(\psi)/N)/d\mu>0.$ Fig. 4(a) demonstrates
that the normalization energy for the even parity wave function (2) increases
when the principal quantum number $n$ increases. So does the odd parity wave
function (3), as shown in Fig. 4(b). It is shown that the energy levels of the
atoms are only associated with the principle quantum number $n$. These are
similar to energy level distribution of the energy eigenvalue problem for the
linear harmonic oscillator described by linear Schrödinger equation.

Figure 4: (color online). The normalization energy $E(\psi)/N$ vs chemical
potential $\mu$, with $N=\int\int dxdy|\psi|^{2}$. (a) even parity wave
function (2) with principal quantum numbers $n=0,2,4$ and (b) odd parity wave
function (3) with $n=1,3,5$. Here the parameters $\omega=0.02$ and $\nu=0.1.$

Stability analysis.—Stability of exact solutions with respect to perturbation
is very important, because only stable localized nonlinear matter waves are
promising for experimental observations and physical applications. To study
the stability of our exact solutions (2)-(3) with Eq. (4), we consider a
perturbed solution $\psi(x,y,t)=[\phi_{n}(x,y)+\Psi(x,y,t)]e^{-i\mu t}$ of Eq.
(1). Here $\phi_{n}(x,y)$ are the exact solutions of the stationary nonlinear
Schrödinger equation
$\frac{1}{2}\phi_{xx}+\frac{1}{2}\phi_{yy}-\frac{1}{2}\omega^{2}(x^{2}+y^{2})\phi-g(x,y)\phi^{3}+\mu\phi=0$.
$\Psi(x,y,t)\ll 1$ is a small perturbation to the exact solutions and
$\Psi(x,y,t)=[R(x,y)+I(x,y)]e^{i\lambda t}$ is decomposed into its real and
imaginary parts Bronski . Substituting this perturbed solution to the quasi-2D
GP equation (1) and neglecting the higher-order terms in $(R,I)$, we obtain a
standard eigenvalue problem $L_{+}R=\lambda I,~{}L_{-}I=\lambda R,$ where
$\lambda$ is eigenvalue, $R,I$ are eigenfunctions with
$L_{+}=-\frac{1}{2}(\partial_{x}^{2}+\partial_{y}^{2})+3g(x,y)\phi_{n}(x,y)^{2}+\frac{1}{2}\omega^{2}(x^{2}+y^{2})-\mu$
and
$L_{-}=-\frac{1}{2}(\partial_{x}^{2}+\partial_{y}^{2})+g(x,y)\phi_{n}(x,y)^{2}+\frac{1}{2}\omega^{2}(x^{2}+y^{2})-\mu.$
Numerical experiments show that when $\omega=0.02$ and $\mu,\nu$ are arbitrary
non-negative constants, only for principle quantum number $n=0,1,2,3,4,5$ are
the eigenvalues $\lambda$ of this eigenvalue problem real. This suggests that
for $\omega=0.02$ the exact localized nonlinear matter wave solution (2) is
linear stability only for $n=0,2,4$ and solution (3) is linear stability only
for $n=1,3,5,$ see Fig. 5. It is seen that when the frequencies of pancake-
shaped trap is fixed, the stability of the exact solutions (2)-(3) with Eq.
(4) rests only on the principle quantum number $n.$

Experimental protocol.—We now provide an experimental protocol for creating
the quasi-2D localized nonlinear matter waves. To do so, we take the
attractive 7Li condensate Hulet1 ; Hulet2 , containing about $10^{3}\sim
10^{5}$ atoms, confined in a pancake-shaped trap with radial frequency
$\omega_{\perp}=2\pi\times 10$ Hz and axial frequency $\omega_{z}=2\pi\times
500$ Hz Rychtarik . This trap can be determined by combination of
spectroscopic observations, direct magnetic field measurement, and the
observed spatial cylindrical symmetry of the trapped atom cloud Rychtarik .
The next step is to realize the spatial variation of the scattering length.
Near the Feshbach resonance Inouye1 ; Xiong ; Ueda1 ; Khaykovich , the
scattering length $a_{s}(B)$ varies dispersively as a function of magnetic
field $B,$ i.e. $a_{s}(B)=a[1+\Delta/(B_{0}-B)],$ with $a$ being the
asymptotic value of the scattering length far from the resonance, $B_{0}$
being the resonant value of the magnetic field, and $\Delta$ being the width
of the resonance. For the magnetic field in $z$ direction with gradient
$\alpha$ along $x$-$y$ direction, we have $\vec{B}=[B_{0}+\alpha
B_{1}(x,y)]\vec{e_{z}}$. In this case, the scattering length is dependent on
$x$ and $y$. In real experiments, the spatially dependent magnetic field may
be generated by a microfabricated ferromagnetic structure integrated on an
atom chip Vengalattore1 ; Vengalattore2 , such that interaction in Fig. 1 can
be realized. In order to observe the density distributions in Figs. 2-3
clearly in experiment, the 7Li atoms should be evaporatively cooled to low
temperatures, say in the range of 10 to 100 $nK$. After the interaction
parameter in Fig. 1(a) is realized by modulating magnetic field properly, the
density distributions in Fig. 2 can be observed for different numbers of atoms
by evaporative cooling, for example, the numbers of atoms in Fig. 2(a)-2(c)
are $3.76\times 10^{3},6.84\times 10^{4},2.633\times 10^{5}$, respectively.
The density distributions in Fig. 3 can also be observed by changing the
scattering lengths through magnetic field for various atom numbers.

Figure 5: (color online). Eigenvalue for different principal quantum number
$n$ with parameters $\omega=0.02,\mu=0.001$ and $\nu=0.1.$ It is shown that
only for $n=0,1,2,3,4,5$ are the localized nonlinear matter wave solutions
(2)-(3) with Eq. (4) linear stability.

Conclusion.—In summary, we have discovered a new family of stable exact
localized nonlinear matter wave solutions of the quasi-2D BEC with spatially
modulated nonlinearities in harmonic potential. Similar to the linear harmonic
oscillator, we introduce two classes of quantum numbers: the principle quantum
number $n$ and secondary quantum number $l$. The matter wave functions have
even parity for the even principle quantum number and odd parity for the odd
one, the energy levels of the atoms are only associated with the principle
quantum number, and the number of density packets for each quantum state is
equal to $(n+1)\times(l+1)$. We also provide an experimental scheme to observe
these novel phenomena in future experiments. Our results are of particular
significance to matter wave management in high dimensional BEC.

This work was supported by NSFC under Grants No. 10874235, No. 10934010, No.
60978019 and by NKBRSFC under Grants No. 2006CB921400, No. 2009CB930704 and
No. 2010CB922904.

## References

  * (1) M. H. Anderson et al., Science 269, 198 (1995).
  * (2) C. C. Bradley et al., Phys. Rev. Lett. 75, 1687 (1995).
  * (3) K. B. Davis et al., Phys. Rev. Lett. 75, 3969 (1995).
  * (4) S. Inouye et al., Nature 392, 151 (1998); J. L. Roberts et al., Phys. Rev. Lett. 81, 5109 (1998).
  * (5) M. Marinescu and L. You, Phys. Rev. Lett. 81, 4596 (1998); M. Theis et al., Phys. Rev. Lett. 93, 123001 (2004).
  * (6) L. Khaykovich et al., Science 296, 1290 (2002); J. Herbig et al., Science 301, 1510 (2003).
  * (7) M. Bartenstein et al., Phys. Rev. Lett. 92, 203201 (2004).
  * (8) M. I. Rodas-Verde et al., Phys. Rev. Lett. 95, 153903 (2005); A. V. Carpentier et al., Phys. Rev. A 74, 013619 (2006).
  * (9) H. Sakaguchi and B. A. Malomed, Phys. Rev. E 72, 046610 (2005); G. Theocharis et al., Phys. Rev. A 72, 033614 (2005); F. K. Abdullaev and J. Garnier, Phys. Rev. A 72, 061605(R) (2005); J. Garnier and F. K. Abdullaev, Phys. Rev. A 74, 013604 (2006); M. A. Porter et al., Physica D 229, 104 (2007); D. A. Zezyulin et al., Phys. Rev. A 76, 013621 (2007).
  * (10) M. Salerno, V. V. Konotop, Y. V. Bludov, Phys. Rev. Lett. 101, 030405 (2008); V. A. Brazhnyi et al., Phys. Rev. Lett. 102, 144101 (2009).
  * (11) G. Dong and B. Hu, Phys. Rev. A 75, 013625 (2007); L. C. Qian et al., Phys. Rev. A 77, 013611 (2008).
  * (12) J. Belmonte-Beitia et al., Phys. Rev. Lett. 98, 064102 (2007); J. Belmonte-Beitia et al., Phys. Rev. Lett. 100, 164102 (2008); A. T. Avelar, D. Bazeia, and W. B. Cardoso, Phys. Rev. E 79, 025602(R) (2009).
  * (13) H. Xiong et al., Phys. Rev. Lett. 95, 120401 (2005).
  * (14) J. M. Gerton et al., Nature London 408, 692 (2000); M. Junker et al., Phys. Rev. Lett. 101, 060406 (2008); S. E. Pollack et al., Phys. Rev. Lett. 102, 090402 (2009).
  * (15) G. Theocharis et al., Phys. Rev. Lett. 90, 120403 (2003).
  * (16) H. Saito, M. Ueda, Phys. Rev. Lett. 90, 040403 (2003).
  * (17) J. J. Garcia-Ripoll et al., Phys. Rev. Lett. 87, 140403 (2001).
  * (18) G. Baym, C. J. Pethick, Phys. Rev. Lett. 76, 6 (1996).
  * (19) M. Abramowitz, I. Stegun, Handbook of Mathematical Functions, Dover, New York (1965).
  * (20) S. L. Cornish et al., Phys. Rev. Lett. 85, 1795 (2000).
  * (21) D. Rychtarik et al., Phys. Rev. Lett. 92, 173003 (2004).
  * (22) P. A. Ruprecht et al., Phys. Rev. A 51, 4704 (1995).
  * (23) J. C. Bronski et al., Phys. Rev. Lett. 86, 1402 (2001).
  * (24) F. K. Abdullaev et al., Phys. Rev. Lett. 90, 230402 (2003).
  * (25) M. Vengalattore et al., J. Appl. Phys. 95, 4404 (2004).
  * (26) M. Vengalattore et al., Eur. Phys. J. D 35, 69 (2005).