# Fleeting Small-scale Surface Magnetic Fields Build the Quiet-Sun Corona

L. P. Chitta Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-
Weg 3, 37077 Göttingen, Germany  S. K. Solanki Max-Planck-Institut für
Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  J.
C. del Toro Iniesta Instituto de Astrofísica de Andalucía (IAA-CSIC), Apartado
de Correos 3004, E-18080 Granada, Spain  Spanish Space Solar Physics
Consortium (S3PC), Spain  J. Woch  Max-Planck-Institut für
Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  D.
Calchetti Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg
3, 37077 Göttingen, Germany  A. Gandorfer Max-Planck-Institut für
Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  J.
Hirzberger  Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-
Weg 3, 37077 Göttingen, Germany  F. Kahil Max-Planck-Institut für
Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  G.
Valori Max-Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3,
37077 Göttingen, Germany  D. Orozco Suárez Instituto de Astrofísica de
Andalucía (IAA-CSIC), Apartado de Correos 3004, E-18080 Granada, Spain
Spanish Space Solar Physics Consortium (S3PC), Spain  H. Strecker Instituto de
Astrofísica de Andalucía (IAA-CSIC), Apartado de Correos 3004, E-18080
Granada, Spain  Spanish Space Solar Physics Consortium (S3PC), Spain  T.
Appourchaux Univ. Paris-Sud, Institut d’Astrophysique Spatiale, UMR 8617,
CNRS, Bâtiment 121, 91405 Orsay Cedex, France  R. Volkmer  Leibniz-Institut
für Sonnenphysik, Schöneckstr. 6, D-79104 Freiburg, Germany  H. Peter Max-
Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077
Göttingen, Germany  S. Mandal Max-Planck-Institut für Sonnensystemforschung,
Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  R. Aznar Cuadrado Max-
Planck-Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077
Göttingen, Germany  L. Teriaca Max-Planck-Institut für Sonnensystemforschung,
Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  U. Schühle Max-Planck-
Institut für Sonnensystemforschung, Justus-von-Liebig-Weg 3, 37077 Göttingen,
Germany  D. Berghmans Solar-Terrestrial Centre of Excellence – SIDC, Royal
Observatory of Belgium, Ringlaan -3- Av. Circulaire, 1180 Brussels, Belgium
C. Verbeeck Solar-Terrestrial Centre of Excellence – SIDC, Royal Observatory
of Belgium, Ringlaan -3- Av. Circulaire, 1180 Brussels, Belgium  A. N. Zhukov
Solar-Terrestrial Centre of Excellence – SIDC, Royal Observatory of Belgium,
Ringlaan -3- Av. Circulaire, 1180 Brussels, Belgium  Skobeltsyn Institute of
Nuclear Physics, Moscow State University, 119991 Moscow, Russia  E. R. Priest
School of Mathematics and Statistics, University of St Andrews, St Andrews,
KY16 9SS, UK

(Received 2023 July 19; Accepted 2023 August 17)

###### Abstract

Arch-like loop structures filled with million Kelvin hot plasma form the
building blocks of the quiet-Sun corona. Both high-resolution observations and
magnetoconvection simulations show the ubiquitous presence of magnetic fields
on the solar surface on small spatial scales of $\sim$100 km. However, the
question of how exactly these quiet-Sun coronal loops originate from the
photosphere and how the magnetic energy from the surface is channeled to heat
the overlying atmosphere is a long-standing puzzle. Here we report high-
resolution photospheric magnetic field and coronal data acquired during the
second science perihelion of Solar Orbiter that reveal a highly dynamic
magnetic landscape underlying the observed quiet-Sun corona. We found that
coronal loops often connect to surface regions that harbor fleeting weaker,
mixed-polarity magnetic field patches structured on small spatial scales, and
that coronal disturbances could emerge from these areas. We suggest that
weaker magnetic fields with fluxes as low as $10^{15}$ Mx and/or those that
evolve on timescales less than 5 minutes, are crucial to understand the
coronal structuring and dynamics.

Solar extreme ultraviolet emission (1493), Solar photosphere (1518), Solar
coronal heating (1989), Solar magnetic fields (1503), Solar magnetic
reconnection (1504), Magnetohydrodynamics (1964), Solar coronal loops (1485)

††journal: ApJL

##  1 Introduction

The outer atmosphere of the Sun, the corona, is best observed from space in
the extreme ultraviolet (EUV) and X-rays. The corona is characterized by the
underlying magnetic activity (Golub et al., 1974; Vaiana & Rosner, 1978). The
so-called active Sun consists of large-scale transient phenomena such as
sunspots, active regions, prominences, flux ropes, flares and coronal mass
ejections (Priest, 2014). By comparison, outside the active Sun we have the
quiet Sun (QS), which varies generally on smaller scales, and whose corona
includes arch-like plasma loops that trace the closed magnetic field (Reale,
2014) and coronal holes, which appear darker in the EUV and X-rays and are
governed by open magnetic fields of the Sun that extend into interplanetary
space (Cranmer, 2009).

This magnetically-closed quiet-Sun corona (simply referred to as the QS corona
in the rest of the paper) is also the most persistent outer atmospheric
feature and is present throughout the 11-year solar activity cycle, with its
global mean temperatures varying from 1.4 MK to 1.8 MK with increasing
magnetic activity (Morgan & Taroyan, 2017). Plasma loops in the QS corona have
a wide range of lengths from the order of 10 Mm to more than 100 Mm (e.g.,
Wiegelmann & Solanki, 2004; Milanović et al., 2023). Some QS coronal loops
might lack a clear spatial structuring as they appear rather diffuse in the
EUV emission (Gorman et al., 2023). Although it is accepted that magnetic
fields govern coronal dynamics in general, the mechanism of coronal heating to
over a million Kelvin remains widely debated (Cranmer & Winebarger, 2019).

Surface magnetic fields, underlying the QS corona, generally reside on spatial
scales of less than 100 km, in the intergranular lanes of turbulent
magnetoconvection (Solanki, 1993). These QS magnetic fields are actually the
most dominant form of the solar surface magnetism, with three orders of
magnitude more magnetic flux emerging in the quiet regions compared to that in
active regions (Stein, 2012). Photospheric spectropolarimetric diagnostics
based on the Hanle effect indicate the existence of ubiquitous tangled
magnetic fields with an average strength over 100 G in the QS (Trujillo Bueno
et al., 2004; Shchukina & Trujillo Bueno, 2011; Bellot Rubio & Orozco Suárez,
2019). A local dynamo process could produce these QS magnetic fields (Vögler &
Schüssler, 2007; Danilovic et al., 2010; Lites, 2011; Rempel, 2014).

Three-dimensional magnetohydrodynamic models of the solar atmosphere that
self-consistently produce QS magnetic fields with average strengths over 100
G, do lead to an $\sim$1 MK hot corona (Amari et al., 2015; Rempel, 2017). In
these models coronal heating is facilitated via the dissipation of Alfvén
waves with lifetimes of 30–50 minutes, larger compared to the typical granular
lifetimes of 5–8 minutes (Amari et al., 2015) -or- through the stressing of
magnetic fields with motions on timescales of 20–50 minutes, fueling
intermittent energy release in the corona (Rempel, 2017). Such models also
reproduce features of some observed transient QS coronal brightenings (Chen et
al., 2021).

However, the correspondence between these simulations and observations in
terms of the spatial structuring and temporal evolution of the QS corona is
not fully established. Moreover, observational questions on how these small-
scale magnetic fields do couple to the solar corona, along with the
understanding of spatial, temporal scales and flux content of magnetic fields
relevant for QS coronal structuring and heating remain open. It is often
difficult and non-trivial to associate QS atmospheric loops to their surface
magnetic footpoints (e.g., Sánchez Almeida et al., 2007). This is mainly
because, until recently there is a lack of joint photospheric and coronal
observations at high-spatial resolutions of better than 500 km (i.e., sub-
granular scales where most of the magnetic fields reside in the photosphere).
Here we use unprecedented joint photospheric and coronal observations from
Solar Orbiter (Müller et al., 2020), during its second science perihelion
passage, to study the magnetic landscape of the QS corona.

Figure 1: Quiet-Sun photosphere and corona. Panel (a): Photospheric magnetic
flux density, close to the disk-center, obtained by SO/PHI on 2022 October 12.
The lighter and darker shaded regions represent positive and negative magnetic
polarity regions, respectively. The map is saturated at flux densities of
$\pm$50 Mx cm-2. The red (positive) and blue (negative) contours outline
magnetic concentrations with flux densities above 30 Mx cm-2 (i.e., 3 times
the noise level of $\sim$10 Mx cm-2). The yellow (positive) and cyan
(negative) contours mark regions with flux densities between 20–30 Mx cm-2
(i.e., 2–3 times the noise level). The green ellipse outlines a prominent
elongated magnetic network feature composed of predominantly negative magnetic
polarity concentrations. Boxes L1–L8 (side length of 10.55 Mm) identify
magnetic field patches underlying the loop features that we analyzed in this
study. Panel (b): QS coronal features captured by the
HRI${}_{\textrm{\scriptsize{EUV}}}$. The displayed snapshot is from the 30 s
effective cadence time series. The green ellipse, contours, and boxes all have
the same meaning as in panel (a). A 5 Mm scale is overlaid for reference. An
animation of panel (b), without the box and ellipse annotations (covering time
period UT 04:40–06:55), is available online. The real-time duration of the
animation is 9 s. See Sect. 2 and Appendix A for details. Figure 2: Small-
scale and weaker magnetic footpoints of QS coronal loops.
HRI${}_{\textrm{\scriptsize{EUV}}}$ sequences of sections of loops L1–L4,
marked in Fig. 1, and contours of the LOS component of the corresponding
magnetic field concentrations at their footpoints, as seen with SO/PHI are
displayed. The magnetic contours have the same meaning as in Fig. 1. Time runs
from left to right. See Sect. 2 and Appendix A for details.

##  2 Observations

We used photospheric line-of-sight magnetic field or flux density observations
and coronal EUV images covering a QS region close to the disk-center. The data
were recorded on 2022 October 12, at which time, Solar Orbiter was at a
distance of 0.293 astronomical units from the Sun. These observations are a
part of the Solar Orbiter Observing Plan, named R_SMALL_HRES_HCAD_RS-burst
(Zouganelis et al., 2020), and are available through the Solar Orbiter
Archive.111https://soar.esac.esa.int/soar

The magnetic field maps were obtained by the High Resolution Telescope (HRT)
of the Polarimetric and Helioseismic Imager on board the Solar Orbiter
spacecraft (SO/PHI; Solanki et al., 2020), which samples the four Stokes
parameters at five wavelength positions within the Fe i 617.3 nm line plus a
sixth one of its nearby continuum (reduction and calibration of SO/PHI data
are detailed in Kahil et al. 2022a, Sinjan et al. 2022, and Kahil et al.
2023).222https://www.mps.mpg.de/solar-physics/solar-orbiter-phi/data-releases
The data were recorded from UT 04:17 to UT 07:12, with a cadence of 300 s, an
image scale of 0.5″ pixel-1 (spatial resolution corresponding to
$2\,\textrm{pixels}\approx 212\,\textrm{km}$ on the Sun). The HRT field of
view (FOV) of 109 Mm $\times$ 109 Mm is sampled by 1024 $\times$ 1024 pixels.

The coronal EUV images were obtained using the High Resolution Imager
(HRI${}_{\textrm{\scriptsize{EUV}}}$) of the Extreme Ultraviolet Imager (EUI;
Rochus et al., 2020) onboard Solar Orbiter. The filter of
HRI${}_{\textrm{\scriptsize{EUV}}}$ is centered at 17.4 nm, and its thermal
response function has a peak at around 1 MK (e.g., Chen et al., 2021). The
passband has contributions from Fe ix (at 17.11 nm) and Fe x (at 17.45 nm and
17.72 nm). In this study we used level-2 data made available as a part of the
EUI Data Release 6 (Kraaikamp et al., 2023). The whole of the
HRI${}_{\textrm{\scriptsize{EUV}}}$ data sequence may be divided into three
sub-sequences. From UT 04:40 to UT 05:25, the images were recorded at a
cadence of 10 s with an exposure time of 5.6 s, then until UT 06:10 at a
cadence of 3 s with an exposure time of 1.65 s, and finally again at 10 s
cadence with exposure time of 5.6 s until about UT 06:55. The EUIPREP pipeline
that produces the level-2 images normalizes the data by the exposure time.
Using these three HRI${}_{\textrm{\scriptsize{EUV}}}$ sub-sequences, we
derived a new series with a uniform effective cadence of 30 s for display
purposes (see Appendix A for details).

The HRI${}_{\textrm{\scriptsize{EUV}}}$ observations have an image scale of
0.492″ pixel-1 (spatial resolution of $\sim$ 209 km), and FOV of $\sim$ 214 Mm
$\times$ 214 Mm, sampled by 2048 $\times$ 2048 pixels (Fig. 5), fully covering
the FOV of SO/PHI HRT. We aligned the HRI${}_{\textrm{\scriptsize{EUV}}}$ and
SO/PHI data, and selected an area of $\sim$ 84 Mm $\times$ 104 Mm for further
analysis. More details on the SO/PHI HRT and
HRI${}_{\textrm{\scriptsize{EUV}}}$ data processing and alignment are
described in the Appendix A.

The longitude of Solar Orbiter in the Heliocentric Earth Ecliptic coordinate
system at the time of observations was $\sim-118$°. Therefore, coordinated
observations from neither Earth-orbiting nor ground-based telescopes are
available.

The photospheric magnetic flux density and the overlying coronal 17.4 nm
emission maps, in the selected FOV are shown in Fig. 1. As expected, the
SO/PHI FOV is dominated by small-scale magnetic fields typical of the QS. The
standard deviation ($\sigma$) of the root-mean-square fluctuations of the
magnetic flux density is estimated to be about 10 Mx cm-2. We consider this to
be the 1-$\sigma$ noise level in the inferred magnetic flux densities. Based
on this, we classify all those pixels that possess magnetic flux densities
above 3-$\sigma$ level, including isolated ones, as regions where the fields
are well-detected. Clear magnetic network patches (known to host kilo Gauss
magnetic flux elements; Stenflo, 1973), distributed throughout the FOV are
distinguishable (blue and red colored contoured regions in Fig. 1). One such
elongated magnetic network feature composed of predominantly negative polarity
magnetic elements is highlighted by an ellipse in Fig. 1a.

Figure 3: Small-scale and weaker magnetic footpoints of QS coronal loops. Same
as Fig. 2, but plotted for sections of loops L5–L8. Ellipses in panels (f) and
(g) outline regions of footpoints of coronal loop bundle L6 without clear
underlying magnetic field structures.

We define weaker field patches as pixels with flux densities in the range of
2-$\sigma$ to 3-$\sigma$ level. In this case, we avoided isolated pixels and
considered only those patches which contain more than 4 contiguous pixels,
each possessing flux densities in the range of 2-$\sigma$ to 3-$\sigma$ level.
We found that these weaker field patches generally surround the network
patches (cyan and yellow colored contoured regions in Fig. 1), nonetheless,
there are also clear cases of isolated weaker field patches. In our
observations, the well-detected patches represent, 2.9% of the pixels in the
considered FOV on average, while the weaker field patches per our criterion,
represent about 1.9% of the pixels.

The overlying coronal emission is generally structured into arch-like loop
features. Many such loops are apparently rooted in the underlying network
patches (see loops outlined by the green ellipse in Fig. 1b). Together, these
SO/PHI and HRI${}_{\textrm{\scriptsize{EUV}}}$ observations reveal the QS
magnetic structures at almost exactly the same high spatial resolution of
$\sim$ 210 km. Using these unprecedented coordinated observations from Solar
Orbiter, we now examine the surface magnetic landscape of the QS corona.

Figure 4: Coronal dynamics emerging from weaker-field regions. Panels (a),
(c), (e), and (g) display examples of four jet-like events from the weak-field
sections of loops, L1, L4, L5, and L6. The green and purple colored boxes are
positioned along the parts of these jet-like activity. The magnetic flux
density contours have the same meaning as in Fig. 1. The green and purple
light curves in panels (b), (d), (f), and (h) show the average
HRI${}_{\textrm{\scriptsize{EUV}}}$ intensities from the corresponding boxes
overlaid on the left-side panels, as functions of time, normalized to their
respective minima and maxima within the plotted time period. These light
curves are derived from the images at their native cadence. Sect. 2 and
Appendix A for details.

##  3 Surface magnetic landscape of the quiet-Sun corona

Many coronal loops or bundles of loops in the considered FOV appear to have
one of their footpoints rooted in network magnetic patches. As such, the
central part of the HRI${}_{\textrm{\scriptsize{EUV}}}$ FOV displayed is
dominated by loops connected to the elongated network patch (see also the
time-averaged maps in Fig. 6). By following the animation associated with Fig.
1, we visually identified some coronal loop systems, labeled L1–L8. For
example, in the displayed instance in Fig. 1, loop features ending in boxes
L3, L5, and L7 are evident.

We present a closer look at one of the footpoint regions or segments of these
loop bundles and their evolution in Figs. 2 and 3. Alongside, we also plot the
contours of underlying magnetic field patches (as defined in Sect. 2) near the
considered footpoints of these loop segments. For instance, in case of loop
L1, only a part of the loop segment is initially seen to be connected to a
well-detected positive polarity magnetic patch of size less than 1 Mm (red
contoured region, over 3-$\sigma$ level, closer to the center of Fig. 2a).
With time, this magnetic patch is no longer distinguishable, while the coronal
loop continues to persist (Fig. 2b–d). This footpoint region is also bustling
with short-living weaker magnetic patches of both positive and negative
polarities.

Similar magnetic field and coronal behavior is observed in loops L2–L4. In
case of loop L2, the larger positive polarity magnetic field patches show some
persistence. But, there is also an increased presence of both well-detected
and weaker negative-polarity magnetic field patches, on small spatial scales,
at that location (Fig. 2e–h). Small-scale magnetic field patches at the bases
of loops L3 and L4 are rather scattered and show little persistence over time
(Fig. 2i–p).

Compared to all the other examples, the loop bundle L5 has one of its
footpoints rooted in magnetic field patches over 3-$\sigma$ level, that are
persistent. Still, these magnetic patches are flanked by fleeting small-scale
mixed-polarity magnetic fields (Fig. 3a–d). The segment of loop L6 is again
rooted in mixed-polarity footpoint region. It is striking to note that a
brighter portion of the footpoint region is devoid of clear underlying
magnetic structures (regions outlined by ellipses in Fig. 3f–g).

Loop features L7 and L8 are more dynamic. In case of loop bundle L7, it is
evident from Fig. 1 that its conjugate footpoint is rooted in the
predominantly negative polarity magnetic network patch. Therefore, it is
natural to expect that the region bounded by box L7 is composed of positive
polarity magnetic concentrations. However, a zoom into this region reveals the
presence of well-detected patches of negative polarity magnetic fields
directly at the footpoint regions along with well-detected positive polarity
fields, in addition to weaker mixed-polarity fields (Fig. 3i–l). Time-averaged
magnetic field map also reveals the persistence of negative polarity magnetic
elements in box L7 over the course of 2.25 hours of observations, along with
positive polarity magnetic elements surrounding them (see box L7 in Fig. 6).
Loop L8 is fully contained within box L8 and is related to an eruption or
larger-scale reconfiguration of magnetic fields in that location as inferred
from the rapidly varying coronal emission (see Fig. 1 and the associated
online animation). This loop bundle shows two clear footpoints on either side,
but the underlying magnetic field distribution is quite random without a clear
bipolar structure (Fig. 3m–p).

We find that while one footpoint of a QS coronal loop bundle might end in
network magnetic field patches, the opposite end or the conjugate footpoint
regions are complex, and are riddled with weaker and small-scale mixed
polarity magnetic fields that evolve on timescales of 5 minutes or less. But
the question is whether these small-scale, often weaker magnetic features play
any active role in the dynamics of the overlying corona.

##  4 Coronal dynamics emerging from weaker field regions

To demonstrate the active role of these weaker small-scale fields at QS
coronal loop footpoints, we present examples of jet-like eruptive activity
emerging from these regions in four of the studied loop bundles (Fig. 4). Near
the footpoints of loop L1, we observed a slender propagating coronal
disturbance that traversed a distance of at least 5 Mm from its source region
(green to purple box locations in Fig. 4a). The propagating nature of the
disturbance is evident from the localized intensity fluctuations: the source
region (green light curve) exhibits intensity enhancements before the
disturbance showing signatures farther away (purple light curve). The overall
disturbance, as inferred from the increase and decrease of localized intensity
fluctuations lasts for about 50–100 s. At this footpoint, we found
impersistent positive polarity magnetic elements. There are no well-detected
negative polarity magnetic patches, per our definition in Sect. 2, within at
least 1 Mm separation from the jet footpoint, at the time of the activity.

Events closer to the footpoint regions of loops L4 and L5 emerge over
locations with seemingly no well-detected underlying magnetic features (Fig.
4). The jet event from L6 region even exhibits a 1–2 Mm wide dome like
structure with an elongated spine feature and appears to emerge from weaker
magnetic field regions, if any.

We detect a continually changing magnetic landscape at the feet of coronal
loops. The data also reveal that the small-scale (weaker) magnetic fields act
not simply as passive photospheric endpoints of the corona, but might play a
role in generating disturbances that propagate into the higher layers.

##  5 Discussion

Our SO/PHI data show that almost all the weaker field patches including the
more isolated well-detected field patches are captured only once in any given
single frame of the 300 s cadence series. Obviously, lifetimes of these
patches must then be below this interval. Moreover, in our case, SO/PHI
scanned the Fe i 617.3 nm line in about 63 s, i.e., each magnetogram is a
temporal average over this period. Therefore, the shortest-lived structures
will appear in only parts of the raw data and they would not show up in the
magnetograms. Indeed, Sunrise observations with 33 s cadence and with nearly
two times better spatial resolution than SO/PHI (when close to the
perihelion), suggest that the median lifetime of weaker internetwork magnetic
features is only 66 s (Anusha et al., 2017). Given that this measurement in
itself might have been limited by the cadence of Sunrise observations, there
could be weaker surface magnetic features on the Sun with still shorter
lifetimes.

Interestingly, our HRI${}_{\textrm{\scriptsize{EUV}}}$ observations do reveal
clear spatial variations in the EUV emission patterns and apparent motions on
timescales as short as 30 s (see movie associated with Fig. 1). Except for the
case L8 that formed during an eruption, each of the other loop bundles as a
whole persist for longer timescales, ranging from a few 10 minutes to the
entire duration of our observations. But the individual strands within these
bundles exhibit intensity variations and evolve on shorter timescales of less
than 5 minutes as displayed in Figs. 2 and 3. While some of these variations
could be related to the general thermal and density changes in the EUV
emitting plasma, some could be directly in response to the rapidly changing
surface magnetic landscape.

Latter of these EUV Variations, that are potentially linked to the rapidly
changing magnetic landscape, can be seen as dynamic events333The high spatio-
temporal resolution offered by HRI${}_{\textrm{\scriptsize{EUV}}}$led to the
studies of a variety of small-scale dynamic and transient coronal events
(e.g., Berghmans et al., 2021; Chitta et al., 2021, 2023; Hou et al., 2021;
Kahil et al., 2022b; Mandal et al., 2021, 2022; Panesar et al., 2021, 2023;
Tiwari et al., 2022). emerging from the weaker field regions. The current
understanding is that the type of jet activity discussed in Fig. 4 are driven
by interaction and reconnection of a parasitic polarity magnetic element with
a dominant, opposite-polarity magnetic concentration, during surface flux
emergence or cancellation events (Shen, 2021). Here we cannot supply evidence
for such reconnection events, though. But the lack of obvious (well-detected)
small-scale parasitic polarity fields at the location of these jets suggests
that they are not being resolved by our magnetic field observations or that
they evolve on timescales shorter than 5 minutes. Given the median lifetime of
1.1 minutes of the internetwork magnetic features, the latter is rather
likely.

The higher-resolution Sunrise observations recorded smallest magnetic fluxes
of $9\times 10^{14}$ Mx (Anusha et al., 2017; Smitha et al., 2017). In our
data for comparison, the smallest value of the magnetic flux for well-detected
fields from a single pixel, that SO/PHI can measure is about $3.5\times
10^{15}$ Mx. Still weaker magnetic features will obviously remain undetected
in our data. Surface magnetic fields on the Sun follow a power-law
distribution over seven orders of magnitude in flux (1016–1023 Mx), with an
index of $-2.69$ (Thornton & Parnell, 2011, see also Anusha et al. 2017 and
Smitha et al. 2017). Assuming that the same power-law index continues to lower
fluxes of $10^{15}$ Mx, magnetic elements with flux contents in the range of
1015–1016 Mx will emerge at 4 times higher rate compared to those in the range
of 1016–1020 Mx (c.f. Eq. 3 in Thornton & Parnell, 2011). This implies that
magnetic elements with flux content below our detection limit (i.e.,
$<3.5\times 10^{15}$ Mx) still could be crucial for coronal loop structuring
and atmospheric dynamics.

In the selected loop bundles, we have mainly focused on the footpoints rooted
in the generally weaker field regions. For instance, the conjugate footpoints,
i.e., the opposite ends of loop bundles L1–L7 (Fig. 1) end in network magnetic
field patches. There, too, the presence of dynamic, weaker mixed polarity
fields is clear. These features might be related to the internetwork elements
that emerge closer to or are swept into the network magnetic concentrations
(Gošić et al., 2014). Their interaction with the network magnetic features
might give rise to enhanced atmospheric activity (Samanta et al., 2019; Wang
et al., 2022; Bose et al., 2023).

Previously, the existence of coronal magnetic connections to photospheric weak
field regions could only be inferred and argued through magnetic field
extrapolations of synthetic magnetograms, representative of the QS (Schrijver
& Title, 2003; Jendersie & Peter, 2006) or observed magnetic fields
(Wiegelmann et al., 2010, 2013). Here we observationally demonstrated that QS
coronal loops do connect to surface regions with rapidly varying, weaker
magnetic fields. Moreover, we also have shown evidence for coronal dynamics
emerging from the underlying weaker field regions.

Here we implicitly assume that the EUV segments of a loop bundle connect
nearly vertically to the underlying magnetic fields. In addition, the EUV
segment of the loop bundle need not necessarily trace the full length/arc of
the magnetic connectivity. While the coronal segment of the loop around 1 MK
emitting in the EUV is detected by HRI${}_{\textrm{\scriptsize{EUV}}}$, the
instrument will not detect/capture the lower parts of the loop filled with
denser and cooler chromospheric material. In such a case, it could very well
be that the apparent EUV endpoints of a coronal loop do not actually point to
the surface magnetic footpoints. Then the magnetic fields threading a loop
could bifurcate closer to the surface and connect to other distinct flux
patches (Priest et al., 2002) or extend further along and connect to surface
regions that are laterally separated from the assumed footpoints. Even should
this be the case, our key interpretation of the presented observations still
stands, namely that the underlying magnetic patches in the vicinity of the
assumed footpoints of the loop bundles are still small-scale dynamic features
possessing generally weaker fields.

On a data-processing level, there could be residual offsets between the
coronal images and photospheric magnetic field maps. We assess that the
alignment accuracy is better than 1–2 Mm, but cannot be worse than 5 Mm in any
given point of the FoV (i.e. the scale shown in Fig. 1). Thus there could be
some offsets unaccounted for between the coronal ends of the loops and the
surface magnetic footpoints. Again, as discussed above, in their immediate
vicinity at the surface, even up to $\pm$5 Mm separation, the magnetic field
remains qualitatively similar. Therefore any minor offsets in the alignment
will not affect our inferences.

The morphology of the HRI${}_{\textrm{\scriptsize{EUV}}}$ loop bundles
generally resembles that of the QS coronal bright points around 1–2 MK
(Madjarska, 2019). Although the thermal response of
HRI${}_{\textrm{\scriptsize{EUV}}}$ has a peak around 1 MK, it also has a
secondary contribution at lower temperatures of $\sim$0.3 MK from O vi (Chen
et al., 2021), sampling emission from the transition region. This means that
the studied loop bundles could actually be the transition region counterparts
of coronal bright points (Tian et al., 2008). Irrespective of their exact
thermal characteristics, the evident dynamic nature of the magnetic fields at
their footpoints would still impact the evolution of any cooler/hotter
atmospheric sections of these loops.

##  6 Implications for the quiet-Sun coronal heating

Our observational findings have important implications for coronal heating
models. For instance, wave-based heating models (Van Doorsselaere et al.,
2020) do not yet consider the Poynting flux and wave energy injection into the
upper atmosphere due to the emergence of small-scale magnetic features, their
rapid variation, interaction with other small-scale features and possible
recycling.

Another traditional model of coronal heating is based on the energy dissipated
in the process of untangling of magnetic fields braided by slow stressing of
footpoints (Parker, 1988; Klimchuk, 2006; Cirtain et al., 2013; Chitta et al.,
2022). In this scenario, a typical QS coronal bright point with loop length of
10 Mm might reach a steady state between the energy injected into the corona
through slow footpoint motions and energy losses through radiation and
conduction, on timescales of 1.4 hours (Parker, 1988). Indeed, this timescale
is close to the coronal magnetic field recycling timescale, as inferred from
the analysis based on low-resolution SOHO/MDI magnetograms (Close et al.,
2004). But the rapid spatio-temporal variations we observed mean that magnetic
fields underlying even the QS loops are even more complex than previously
thought, further implying that coronal fields might be subjected to rapid
tangling and untangling, on shorter timescales of only a few 100 s or even
less. In essence, the coronal segment of a loop might not retain its surface
magnetic identity, for a slow build up of energy in the corona over a 1.4 hour
timescale, to begin with. Then how the processes of rapid untangling and the
slow build up of energy compete in generating coronal heating and
corresponding observational signatures in the braiding scenario, needs to be
evaluated.

At the same time, complex magnetic fields might reconnect in the lower
atmosphere and release the energy already in the chromosphere and transition
region at the coronal base, leading to hotter plasma filling the coronal loops
(Aschwanden et al., 2007). If these variations are sustained by frequent flux
emergence and cancellation events in the photosphere (Smitha et al., 2017),
they could power the overlying coronal loops through magnetic reconnection
(Priest et al., 2018, which is indeed observed in some active region loops,
e.g., Tiwari et al. 2014, Chitta et al. 2017, 2020).

As one can see in Figure 1, the correspondence between QS coronal loop bundles
and photospheric magnetic field patches is generally clear on larger scales.
However, such a correspondence for each individual small-scale loop is far
from certain. Even in cases with an obvious overall correspondence, like
regions L3 and L5 (see Figures 2i–l and 3a–d), it is not evident which small-
scale loop connects to which photospheric magnetic field patch. On the one
hand, this may indicate that the heating that leads to the EUV emission is
located directly in the corona and not in the transition region or even lower.
On the other hand, this could also point to the role of unresolved mixed
polarity magnetic elements at the base of coronal loops.

Indeed, high-resolution simulations of the surface magnetoconvection show the
existence of closely spaced mixed polarity magnetic fields in intergranular
lanes (Vögler & Schüssler, 2007; Rempel, 2014). Depending on their relative
strengths and spatial separation, their observable effect, within the
resolution element, could be diminished in our SO/PHI observations. A direct
indication of such an effect could be seen in our
HRI${}_{\textrm{\scriptsize{EUV}}}$ data where the loops apparently root in
regions with little surface magnetic flux (see Fig. 3). Therefore, it is
imperative that the role of (weaker) small-scale surface magnetic fields along
with their recycling, emergence and cancellation need to be carefully
evaluated in the processes of mass and energy transfer into the corona (see
also Wang, 2022). It is also interesting to note that despite the dynamism of
surface magnetic fields, there is a clear persistence in the overlying coronal
loops in the course of about 2.25 hours (Appendix B; Fig. 6). This kind of
temporal persistence could be an additional constraint on different heating
models.

Using unprecedented photospheric and coronal coordinated observations with
SO/PHI and EUI instruments onboard Solar Orbiter, we have demonstrated the
direct and crucial role of small-scale magnetic fields in the structuring of
the QS solar corona. The rapidly varying yet dominant weaker magnetic flux
patches could be important for the overlying atmospheric dynamics and
eventually for the powering of the solar corona itself.

The authors thank the anonymous referee for constructive comments that helped
improve the manuscript. L.P.C. gratefully acknowledges funding by the European
Union (ERC, ORIGIN, 101039844). Views and opinions expressed are however those
of the author(s) only and do not necessarily reflect those of the European
Union or the European Research Council. Neither the European Union nor the
granting authority can be held responsible for them. Solar Orbiter is a space
mission of international collaboration between ESA and NASA, operated by ESA.
We are grateful to the ESA SOC and MOC teams for their support. The EUI
instrument was built by CSL, IAS, MPS, MSSL/UCL, PMOD/WRC, ROB, LCF/IO with
funding from the Belgian Federal Science Policy Office (BELSPO/PRODEX PEA
4000112292 and 4000134088); the Centre National d’Etudes Spatiales (CNES); the
UK Space Agency (UKSA); the Bundesministerium für Wirtschaft und Energie
(BMWi) through the Deutsches Zentrum für Luft- und Raumfahrt (DLR); and the
Swiss Space Office (SSO). The German contribution to SO/PHI is funded by the
BMWi through DLR and by MPG central funds. The Spanish contribution is funded
by AEI/MCIN/10.13039/501100011033/ and European Union “NextGenerationEU”/PRTR”
(RTI2018-096886-C5, PID2021-125325OB-C5, PCI2022-135009-2, PCI2022-135029-2)
and ERDF “A way of making Europe”; “Center of Excellence Severo Ochoa” awards
to IAA-CSIC (SEV-2017-0709, CEX2021-001131-S); and a Ramón y Cajal fellowship
awarded to DOS. The French contribution is funded by CNES. A.N.Z. thanks the
Belgian Federal Science Policy Office (BELSPO) for the provision of financial
support in the framework of the PRODEX Programme of the European Space Agency
(ESA) under contract number 4000136424.

##  Appendix A Data processing and alignment

Figure 5: Quiet-Sun corona and plasma loops. The snapshot shows the full FOV
of HRI${}_{\textrm{\scriptsize{EUV}}}$. The white box marks a partial FOV that
we presented in Fig. 1. See Sect. 2 and Appendix A for details.

The EUI and SO/PHI data were further processed as follows. First we co-aligned
all the EUI data to the first image of the time sequence using a cross-
correlation technique described in Chitta et al. (2022). We used these aligned
level-2 EUI data at their native cadence to compute the light curves plotted
in Fig. 4. Next, the aligned EUI data (composed of sub-sequences with 10 s, 3
s, and 10 s cadences) were further processed to a uniform effective cadence of
30 s, by temporally averaging an appropriate number of images.

We also created an array of time-stamps with 30 s interval to match the 30 s
effective cadence sequence. Here the start time of the new time-stamps is the
same as the time-stamp of the first EUI image in the original L2 sequence.

The 10 s cadence sub-sequences were recorded with an Aluminum filter placed in
front of the detector, while the 3 s cadence sub-sequence images were recorded
with no filter. Consequently, per exposure time, roughly 1/0.6 more photons
were being recorded by the detector with no filter compared to the Aluminum
filter case. This means that the exposure normalized intensities in the 3 s
cadence images are brighter by a factor 1/0.6 compared to the images from the
10 s cadence sub-sequences. For display purposes, we account for these
intensity differences by further normalizing the 30 s effective cadence data
with the corresponding mean intensity at each time step, computed over the
field of view displayed in Fig. 1. Moreover, we applied an unsharp mask filter
to sharpen the images displayed in Figs. 1, 2, 3, and 4. The full FOV of
HRI${}_{\textrm{\scriptsize{EUV}}}$ displayed in Fig. 5 is also from the 30 s
effective cadence data, but with no unsharp mask filter applied.

The SO/PHI magnetic field vector is derived using the Milne-Eddington
approximation under the assumption that the Fe i 617.3 nm spectral line is
formed in local thermodynamic equilibrium (LTE; Solanki et al., 2020).
However, this line has been shown to be affected by non-LTE conditions in the
solar atmosphere. Ignoring these effects might result in an underestimation of
magnetic fields (Smitha et al., 2023). Such non-LTE effects, however, are not
considered in this study.

We re-scaled the SO/PHI maps to the slightly smaller pixel scale of
HRI${}_{\textrm{\scriptsize{EUV}}}$, that resulted in 1041$\times$1041 pixels
data. These re-scaled SO/PHI maps were all co-aligned to the map recorded
closest to the middle of the time sequence (around UT 05:47). Then the first
image from the 30 s effective cadence HRI${}_{\textrm{\scriptsize{EUV}}}$ time
sequence was visually aligned to the closest-in-time re-scaled SO/PHI map. Our
visual alignment was primarily guided by a prominent elongated magnetic
network feature consisting of negative polarity magnetic patches (enclosed by
a green ellipse in Fig. 1a). Conspicuous coronal structures associated with
this magnetic network are seen in the HRI${}_{\textrm{\scriptsize{EUV}}}$ data
(green ellipse Fig. 1b). Keeping re-scaled aligned SO/PHI data as reference,
this visual alignment was achieved by extracting the
HRI${}_{\textrm{\scriptsize{EUV}}}$ sub-FOV with $(x,y)=(340,140)$ as origin
(Fig. 5). This alignment also naturally led to the alignment of other clearer
loop systems to their underlying magnetic fields (e.g., Loop bundle associated
with L5 and a bipolar loop structure at the bottom of the frame, adjacent to
L8). Overall, we assess that the alignment accuracy is better than 1–2 Mm.
Figure 1 then represents a still smaller sub-FOV covering $x=$100–899 pixels
and $y=$17–1016 pixels in the re-scaled SO/PHI data and the corresponding
extracted HRI${}_{\textrm{\scriptsize{EUV}}}$ sub-FOV. Finally, we padded the
SO/PHI time sequence by populating the temporal dimension with the
corresponding maps that are nearest in time to the 30 s effective cadence EUI
sequence.

##  Appendix B Average Maps

Figure 6: Persistent surface magnetic and coronal loop features. Time-averaged
SO/PHI magnetic flux density (panel a) and the
HRI${}_{\textrm{\scriptsize{EUV}}}$ coronal (panel b) maps are displayed.
Green ellipse and boxes L1–L8 have the same meaning as in Fig. 1. The red
(positive) and blue (negative) contours outline magnetic concentrations with
flux densities above 6 Mx cm-2 (i.e., 3-$\sigma_{\textrm{avg}}$), while the
yellow (positive) and cyan (negative) contours mark regions with flux
densities between 2-$\sigma_{\textrm{avg}}$ and 3-$\sigma_{\textrm{avg}}$
level, based on the time-averaged SO/PHI magnetic flux density map. See
Appendices A and B for details.

To check the persistence of surface magnetic fields and the overlying QS
coronal loops, we created temporal averages of SO/PHI magnetic flux density
and HRI${}_{\textrm{\scriptsize{EUV}}}$ coronal maps from the respective time
sequences spanning the length of the EUI observations (i.e. about 135
minutes). There are in total 27 SO/PHI magnetic field maps covering this
period, which we temporally average. The 1-$\sigma$ noise level in the
magnetic flux density in a single map is about 10 Mx cm-2. Assuming Poisson
statistics, the 1-$\sigma$ noise level in the time-averaged magnetic flux
density will reduce by a factor of $\sqrt{27}$, to $\sim$1.9 Mx cm-2 (referred
to as $\sigma_{\textrm{avg}}$ for brevity in the following). The averaged
SO/PHI map is displayed in Figure 6a. The prominent elongated magnetic network
is evident in this map (encircled by an ellipse; see also Fig. 1).

There are 1428 HRI${}_{\textrm{\scriptsize{EUV}}}$ images in the aligned
level-2 EUI native-cadence sequence, which are temporally averaged. We
contrast-enhance the time-averaged coronal map by a multi-scale Gaussian
normalization technique (Morgan & Druckmüller, 2014) to better visualize the
features (Fig. 6b). Persistent QS coronal loops rooted in the elongated
magnetic network are evident in this time-averaged map.

## References

  * Amari et al. (2015) Amari, T., Luciani, J.-F., & Aly, J.-J. 2015, Nature, 522, 188, doi: 10.1038/nature14478
  * Anusha et al. (2017) Anusha, L. S., Solanki, S. K., Hirzberger, J., & Feller, A. 2017, A&A, 598, A47, doi: 10.1051/0004-6361/201527738
  * Aschwanden et al. (2007) Aschwanden, M. J., Winebarger, A., Tsiklauri, D., & Peter, H. 2007, ApJ, 659, 1673, doi: 10.1086/513070
  * Bellot Rubio & Orozco Suárez (2019) Bellot Rubio, L., & Orozco Suárez, D. 2019, Living Reviews in Solar Physics, 16, 1, doi: 10.1007/s41116-018-0017-1
  * Berghmans et al. (2021) Berghmans, D., Auchère, F., Long, D. M., et al. 2021, A&A, 656, L4, doi: 10.1051/0004-6361/202140380
  * Bose et al. (2023) Bose, S., Nóbrega-Siverio, D., De Pontieu, B., & Rouppe van der Voort, L. 2023, ApJ, 944, 171, doi: 10.3847/1538-4357/acb544
  * Chen et al. (2021) Chen, Y., Przybylski, D., Peter, H., et al. 2021, A&A, 656, L7, doi: 10.1051/0004-6361/202140638
  * Chitta et al. (2020) Chitta, L. P., Peter, H., Priest, E. R., & Solanki, S. K. 2020, A&A, 644, A130, doi: 10.1051/0004-6361/202039099
  * Chitta et al. (2017) Chitta, L. P., Peter, H., Solanki, S. K., et al. 2017, ApJS, 229, 4, doi: 10.3847/1538-4365/229/1/4
  * Chitta et al. (2021) Chitta, L. P., Solanki, S. K., Peter, H., et al. 2021, A&A, 656, L13, doi: 10.1051/0004-6361/202141683
  * Chitta et al. (2022) Chitta, L. P., Peter, H., Parenti, S., et al. 2022, A&A, 667, A166, doi: 10.1051/0004-6361/202244170
  * Chitta et al. (2023) Chitta, L. P., Zhukov, A. N., Berghmans, D., et al. 2023, Science, 381, 867, doi: 10.1126/science.ade5801
  * Cirtain et al. (2013) Cirtain, J. W., Golub, L., Winebarger, A. R., et al. 2013, Nature, 493, 501, doi: 10.1038/nature11772
  * Close et al. (2004) Close, R. M., Parnell, C. E., Longcope, D. W., & Priest, E. R. 2004, ApJ, 612, L81, doi: 10.1086/424659
  * Cranmer (2009) Cranmer, S. R. 2009, Living Reviews in Solar Physics, 6, 3, doi: 10.12942/lrsp-2009-3
  * Cranmer & Winebarger (2019) Cranmer, S. R., & Winebarger, A. R. 2019, ARA&A, 57, 157, doi: 10.1146/annurev-astro-091918-104416
  * Danilovic et al. (2010) Danilovic, S., Schüssler, M., & Solanki, S. K. 2010, A&A, 513, A1, doi: 10.1051/0004-6361/200913379
  * Golub et al. (1974) Golub, L., Krieger, A. S., Silk, J. K., Timothy, A. F., & Vaiana, G. S. 1974, ApJ, 189, L93, doi: 10.1086/181472
  * Gorman et al. (2023) Gorman, J., Chitta, L. P., Peter, H., et al. 2023, arXiv e-prints, arXiv:2308.01640.  https://arxiv.org/abs/2308.01640
  * Gošić et al. (2014) Gošić, M., Bellot Rubio, L. R., Orozco Suárez, D., Katsukawa, Y., & del Toro Iniesta, J. C. 2014, ApJ, 797, 49, doi: 10.1088/0004-637X/797/1/49
  * Hou et al. (2021) Hou, Z., Tian, H., Berghmans, D., et al. 2021, ApJ, 918, L20, doi: 10.3847/2041-8213/ac1f30
  * Jendersie & Peter (2006) Jendersie, S., & Peter, H. 2006, A&A, 460, 901, doi: 10.1051/0004-6361:20066166
  * Kahil et al. (2022a) Kahil, F., Gandorfer, A., Hirzberger, J., et al. 2022a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12180, Space Telescopes and Instrumentation 2022: Optical, Infrared, and Millimeter Wave, ed. L. E. Coyle, S. Matsuura, & M. D. Perrin, 121803F, doi: 10.1117/12.2628942
  * Kahil et al. (2022b) Kahil, F., Hirzberger, J., Solanki, S. K., et al. 2022b, A&A, 660, A143, doi: 10.1051/0004-6361/202142873
  * Kahil et al. (2023) Kahil, F., Gandorfer, A., Hirzberger, J., et al. 2023, A&A, 675, A61, doi: 10.1051/0004-6361/202346033
  * Klimchuk (2006) Klimchuk, J. A. 2006, Sol. Phys., 234, 41, doi: 10.1007/s11207-006-0055-z
  * Kraaikamp et al. (2023) Kraaikamp, E., Gissot, S., Stegen, K., et al. 2023, SolO/EUI Data Release 6.0 2023-01, https://doi.org/10.24414/z818-4163
  * Lites (2011) Lites, B. W. 2011, ApJ, 737, 52, doi: 10.1088/0004-637X/737/2/52
  * Madjarska (2019) Madjarska, M. S. 2019, Living Reviews in Solar Physics, 16, 2, doi: 10.1007/s41116-019-0018-8
  * Mandal et al. (2021) Mandal, S., Peter, H., Chitta, L. P., et al. 2021, A&A, 656, L16, doi: 10.1051/0004-6361/202142041
  * Mandal et al. (2022) Mandal, S., Chitta, L. P., Peter, H., et al. 2022, A&A, 664, A28, doi: 10.1051/0004-6361/202243765
  * Milanović et al. (2023) Milanović, N., Chitta, L. P., & Peter, H. 2023, A&A, 673, A81, doi: 10.1051/0004-6361/202245544
  * Morgan & Druckmüller (2014) Morgan, H., & Druckmüller, M. 2014, Sol. Phys., 289, 2945, doi: 10.1007/s11207-014-0523-9
  * Morgan & Taroyan (2017) Morgan, H., & Taroyan, Y. 2017, Science Advances, 3, e1602056, doi: 10.1126/sciadv.1602056
  * Müller et al. (2020) Müller, D., St. Cyr, O. C., Zouganelis, I., et al. 2020, A&A, 642, A1, doi: 10.1051/0004-6361/202038467
  * Panesar et al. (2023) Panesar, N. K., Hansteen, V. H., Tiwari, S. K., et al. 2023, ApJ, 943, 24, doi: 10.3847/1538-4357/aca1c1
  * Panesar et al. (2021) Panesar, N. K., Tiwari, S. K., Berghmans, D., et al. 2021, ApJ, 921, L20, doi: 10.3847/2041-8213/ac3007
  * Parker (1988) Parker, E. N. 1988, ApJ, 330, 474, doi: 10.1086/166485
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge University Press, Cambridge), doi: 10.1017/CBO9781139020732
  * Priest et al. (2018) Priest, E. R., Chitta, L. P., & Syntelis, P. 2018, ApJ, 862, L24, doi: 10.3847/2041-8213/aad4fc
  * Priest et al. (2002) Priest, E. R., Heyvaerts, J. F., & Title, A. M. 2002, ApJ, 576, 533, doi: 10.1086/341539
  * Reale (2014) Reale, F. 2014, Living Reviews in Solar Physics, 11, 4, doi: 10.12942/lrsp-2014-4
  * Rempel (2014) Rempel, M. 2014, ApJ, 789, 132, doi: 10.1088/0004-637X/789/2/132
  * Rempel (2017) —. 2017, ApJ, 834, 10, doi: 10.3847/1538-4357/834/1/10
  * Rochus et al. (2020) Rochus, P., Auchère, F., Berghmans, D., et al. 2020, A&A, 642, A8, doi: 10.1051/0004-6361/201936663
  * Samanta et al. (2019) Samanta, T., Tian, H., Yurchyshyn, V., et al. 2019, Science, 366, 890, doi: 10.1126/science.aaw2796
  * Sánchez Almeida et al. (2007) Sánchez Almeida, J., Teriaca, L., Sütterlin, P., et al. 2007, A&A, 475, 1101, doi: 10.1051/0004-6361:20078124
  * Schrijver & Title (2003) Schrijver, C. J., & Title, A. M. 2003, ApJ, 597, L165, doi: 10.1086/379870
  * Shchukina & Trujillo Bueno (2011) Shchukina, N., & Trujillo Bueno, J. 2011, ApJ, 731, L21, doi: 10.1088/2041-8205/731/1/L21
  * Shen (2021) Shen, Y. 2021, Proceedings of the Royal Society of London Series A, 477, 217, doi: 10.1098/rspa.2020.0217
  * Sinjan et al. (2022) Sinjan, J., Calchetti, D., Hirzberger, J., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12189, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 121891J, doi: 10.1117/12.2629323
  * Smitha et al. (2017) Smitha, H. N., Anusha, L. S., Solanki, S. K., & Riethmüller, T. L. 2017, ApJS, 229, 17, doi: 10.3847/1538-4365/229/1/17
  * Smitha et al. (2023) Smitha, H. N., van Noort, M., Solanki, S. K., & Castellanos Durán, J. S. 2023, A&A, 669, A144, doi: 10.1051/0004-6361/202245130
  * Solanki (1993) Solanki, S. K. 1993, Space Sci. Rev., 63, 1, doi: 10.1007/BF00749277
  * Solanki et al. (2020) Solanki, S. K., del Toro Iniesta, J. C., Woch, J., et al. 2020, A&A, 642, A11, doi: 10.1051/0004-6361/201935325
  * Stein (2012) Stein, R. F. 2012, Living Reviews in Solar Physics, 9, 4, doi: 10.12942/lrsp-2012-4
  * Stenflo (1973) Stenflo, J. O. 1973, Sol. Phys., 32, 41, doi: 10.1007/BF00152728
  * Thornton & Parnell (2011) Thornton, L. M., & Parnell, C. E. 2011, Sol. Phys., 269, 13, doi: 10.1007/s11207-010-9656-7
  * Tian et al. (2008) Tian, H., Curdt, W., Marsch, E., & He, J. 2008, ApJ, 681, L121, doi: 10.1086/590410
  * Tiwari et al. (2014) Tiwari, S. K., Alexander, C. E., Winebarger, A. R., & Moore, R. L. 2014, ApJ, 795, L24, doi: 10.1088/2041-8205/795/1/L24
  * Tiwari et al. (2022) Tiwari, S. K., Hansteen, V. H., De Pontieu, B., Panesar, N. K., & Berghmans, D. 2022, ApJ, 929, 103, doi: 10.3847/1538-4357/ac5d46
  * Trujillo Bueno et al. (2004) Trujillo Bueno, J., Shchukina, N., & Asensio Ramos, A. 2004, Nature, 430, 326, doi: 10.1038/nature02669
  * Vaiana & Rosner (1978) Vaiana, G. S., & Rosner, R. 1978, ARA&A, 16, 393, doi: 10.1146/annurev.aa.16.090178.002141
  * Van Doorsselaere et al. (2020) Van Doorsselaere, T., Srivastava, A. K., Antolin, P., et al. 2020, Space Sci. Rev., 216, 140, doi: 10.1007/s11214-020-00770-y
  * Vögler & Schüssler (2007) Vögler, A., & Schüssler, M. 2007, A&A, 465, L43, doi: 10.1051/0004-6361:20077253
  * Wang et al. (2022) Wang, J., Lee, J., Liu, C., Cao, W., & Wang, H. 2022, ApJ, 924, 137, doi: 10.3847/1538-4357/ac374e
  * Wang (2022) Wang, Y. M. 2022, Sol. Phys., 297, 129, doi: 10.1007/s11207-022-02060-y
  * Wiegelmann & Solanki (2004) Wiegelmann, T., & Solanki, S. K. 2004, Sol. Phys., 225, 227, doi: 10.1007/s11207-004-3747-2
  * Wiegelmann et al. (2010) Wiegelmann, T., Solanki, S. K., Borrero, J. M., et al. 2010, ApJ, 723, L185, doi: 10.1088/2041-8205/723/2/L185
  * Wiegelmann et al. (2013) —. 2013, Sol. Phys., 283, 253, doi: 10.1007/s11207-013-0249-0
  * Zouganelis et al. (2020) Zouganelis, I., De Groof, A., Walsh, A. P., et al. 2020, A&A, 642, A3, doi: 10.1051/0004-6361/202038445


# Quantum soliton-trains of strongly correlated impurities in Bose–Einstein
condensates

Hoshu Hiyane<EMAIL_ADDRESS>Quantum Systems Unit, Okinawa Institute of
Science and Technology Graduate University, Onna, Okinawa 904-0495, Japan
Thomas Busch  Quantum Systems Unit, Okinawa Institute of Science and
Technology Graduate University, Onna, Okinawa 904-0495, Japan  Thomás Fogarty
Quantum Systems Unit, Okinawa Institute of Science and Technology Graduate
University, Onna, Okinawa 904-0495, Japan

###### Abstract

Strongly correlated impurities immersed in a Bose–Einstein condensate (BEC)
can form a periodic structure of tightly localized single atoms due to
competing inter- and intra-species interactions, leading to a self-organized
pinned state. In this work, we show numerically that the impurities in the
self-pinned state form a soliton-train, as a consequence of a BEC-mediated
attractive self-interaction and ordering due to the exclusion principle. The
dynamics of the impurities possess similar characteristics to bright matter-
wave solitons as they appear in attractive BECs, however in the few impurities
case, the detailed nature of collisions is determined by their quantum
statistics.

Solitons, and more generally nondispersive solitary waves, appear in many
physical systems. The historical and paradigmatic example are shallow water
waves, whose dynamics can be modelled by the Korteweg–De Vries equation, which
is well-known to possess soliton solutions Lamb (1980). Over the last couple
of decades dark and bright solitons have been extensively studied in nonlinear
optical systems Zakharov and Shabat (1972); Emplit  _et al._ (1987); Krökel
_et al._ (1988); Kivshar and Malomed (1989); Kivshar (1993); Malomed  _et al._
(2005) and in atomic Bose–Einstein condensates (BECs) Burger  _et al._ (1999);
Busch and Anglin (2000); Carr  _et al._ (2000a, b); Denschlag  _et al._
(2000); Khaykovich  _et al._ (2002); Strecker  _et al._ (2002); Marchant  _et
al._ (2013); Kartashov  _et al._ (2019), which are both systems that are well
described by classical fields whose dynamics is governed by the nonlinear
Schrödinger equation (NLSE). As atomic BECs are amenable to clean and highly
controllable experiments many different realisations of solitons have been
explored, ranging from scalar dark Burger  _et al._ (1999); Busch and Anglin
(2000); Denschlag  _et al._ (2000) and bright Carr  _et al._ (2000b); Strecker
_et al._ (2002) solitons to vector dark-dark Hoefer  _et al._ (2011); Yan  _et
al._ (2012) and dark-bright Busch and Anglin (2001); Nistazakis  _et al._
(2008); Becker  _et al._ (2008); Hamner  _et al._ (2011) solitons. More
recently, complex soliton structures in two-dimensional two-component BECs,
called Townes solitons, have been experimentally observed Bakkali-Hassani  _et
al._ (2021).

Beyond single solitons, highly excited soliton arrays, known as soliton-
trains, have also been theoretically and experimentally investigated. While
the exact solutions of the free-space NLSE contains such states for both
repulsive (dark) Carr  _et al._ (2000a) and attractive (bright) Carr  _et al._
(2000b) interaction, their experimental realization is not easy. While more
controlled ways to create bright solitons have been suggested Edmonds  _et
al._ (2018), bright solitons are usually created by first preparing a stable
condensate with repulsive interactions and then suddenly quenching the
scattering length from positive to negative. However, this excites a
modulational instability that results in the cloud breaking up into bright
soliton-trains Strecker  _et al._ (2002); Carr and Brand (2004). This process
is therefore inherently uncontrollable and results in a collection of bright
solitons of differing widths and particle numbers, making single solitons and
soliton-train states hard to study systematically. In nonlinear optics, on the
other hand, the controllable generation of temporal soliton-trains is possible
using ultrashort laser pulses with high repetition rates Mollenauer  _et al._
(1990); Mamyshev  _et al._ (1991); Haus and Wong (1996); Kivshar and Luther-
Davies (1998); Buryak  _et al._ (2002); Chen  _et al._ (2012). The ability to
deterministically create and control matter-wave soliton trains could have
important applications in quantum engineering, such as in high-precision
sensing with atom interferometry Ramanathan  _et al._ (2011); Martin and
Ruostekoski (2012); Polo and Ahufinger (2013); Wales  _et al._ (2020); Szigeti
_et al._ (2020).

Recently it has been suggested that strongly repulsive bosonic impurities
immersed in a BEC can localize and self-organize into a periodic atomic array
due to the competition between impurity-impurity and impurity-BEC interactions
Keller  _et al._ (2022). In this work, we show that such self-localized
impurities form soliton-train states in the regime where their velocity is
below the speed of sound of the BEC. Such bright soliton-trains can be
understood as the spatial matter-wave counterpart of temporal optical solitons
in nonlinear media, with the nonlinear coupling mediated by the BEC matter-
wave. However, contrary to the similarities with optics, here only few atoms
are sufficient for these nonlinearities to appear Bruderer  _et al._ (2008);
Keller  _et al._ (2022) and the solitonic impurities can be treated fully
quantum mechanically. The system therefore supports quantum soliton-train
states, and we show that quantum properties, such as the statistics of the
impurity atoms (bosonic or fermionic), have a strong influence on their
nonequillibrium properties.

In the following we first briefly review the system of two coupled one-
dimensional quantum gases at ultralow temperatures where the first component
is a weakly correlated atomic BEC described in the mean-field limit, while the
second one is a minority component of $N\lll N_{\text{BEC}}$ bosonic
impurities Keller  _et al._ (2022). The impurities are assumed to strongly
interact with each other so that their many-body dynamics can be described
using the Tonks–Girardeau (TG) gas model Girardeau (1960). The inter-component
interaction is described by a repulsive density-density coupling of strength
$\gamma$, which we scale relative to the strength of the mean-field
interaction of the condensate $g$. We assume equal masses $m$ for both
components with the dynamics of the coupled system described by

$\displaystyle i\dot{\Psi}=$
$\displaystyle\left[-\frac{1}{2}\partialderivative[2]{x}+V_{\rm
BEC}+|\Psi|^{2}+\gamma\rho\right]\Psi,$ (1) $\displaystyle i\dot{\phi}_{n}=$
$\displaystyle\bqty{-\frac{1}{2}\partialderivative[2]{x}+V_{\rm
TG}+\gamma|\Psi|^{2}}\phi_{n}.$ (2)

Here $\Psi(x,t)$ is the condensate order parameter obeying mean-field
Gross–Pitaevskii equation (GPE) Pitaevskii and Stringari (2016) and
$\phi_{n}(x,t)$ are the single particle states needed to describe the TG gas
using the Bose–Fermi mapping theorem Girardeau (1960). All quantities are
scaled with respect to the characteristic length scale $x_{0}=\hbar^{2}/(mg)$,
which is related to the condensate-healing length, and time scale
$t_{0}=mg^{2}/\hbar^{3}$, so that all wavefunctions are in units of
$\sqrt{x}_{0}$. Since the density of the TG gas, $\rho(x,t)$, is equivalent to
the density of a gas of spin-polarised fermions, at zero temperature it can
simply be written as $\rho(x,t)=\sum^{N}_{n=1}|\phi_{n}(x,t)|^{2}$. We assume
the BEC to be trapped in a ring potential of length $L_{\rm BEC}$ so that it
has a flat density profile, while the TG gas is confined in a smaller box trap
$V_{\rm TG}(x)$ of length $L<L_{\rm BEC}$ Keller  _et al._ (2022). This setup
allows us to explore the motion of the impurities without considering edge
effects from the condensate.

The ground state of the coupled system can be found by solving the coupled
equations in a self-consistent manner by means of imaginary time evolution for
Eq. (1) and exact diagonalization for Eq. (2). When the systems are decoupled
for $\gamma=0$, the TG gas is delocalized in the box potential and may be
considered as a quasi-superfluid state. For sufficiently strong coupling
strength $\gamma$ the impurity atoms localize individually in the mean-field
potential provided by the BEC. The many-body groundstate then becomes a
regularly spaced array of pinned impurities (see Fig. 1 (a)). This array of
impurities therefore represents a self-pinned insulator state within a matter-
wave lattice, and the distance between impurities is set by the Fermi momentum
$k_{F}=N\pi/L$ Miyakawa  _et al._ (2004); Keller  _et al._ (2022).

In this self-pinned state, an effective model for the impurity eigenstates and
their dynamics can be constructed by taking the Thomas–Fermi approximation for
the BEC component. The solution to Eq. (1) is then given by Stringari (1996)

$\displaystyle\Psi(x,t)\approx\sqrt{\mu_{\rm TF}-\gamma\rho(x,t)},$ (3)

where the equilibrium chemical potential is $\mu_{\rm
TF}=\mu_{0}\quantity(1+\gamma N/N_{\rm BEC}),$ with $\mu_{0}=N_{\rm
BEC}/L_{\rm BEC}$ being the chemical potential evaluated at $\gamma=0$ Keller
_et al._ (2022). Under this approximation the Schrödinger equation for the
impurities Eq. (2) can be rewritten as

$i\dot{\phi}_{n}(x,t)=\bqty{-\frac{1}{2}\partialderivative[2]{x}-\gamma^{2}\rho(x,t)}\phi_{n}(x,t),$
(4)

with the impurities coupled to one another through the effective attraction
term $-\gamma^{2}\rho(x,t)$. Equation (4) resides in the large group of
integrable Manakov equations Manakov (1973), which possess (bright) soliton-
trains as stationary solutions that are thoroughly investigated in the field
of nonlinear optics Gordon (1983); Radhakrishnan  _et al._ (1997); Kanna  _et
al._ (2006). Since the stationary self-pinned solutions of Eqs. (1) and (2)
show that no overlap between neighbouring impurities exists (see Fig. 1 (a)),
one can treat them separately and therefore simply solve an effective single
impurity NLSE with attractive self-interaction (taking $N=1$ in Eq. (4)), for
which the solution is the celebrated bright soliton:
$\tilde{\rho}(x)=(\gamma/2)^{2}\sech^{2}(\gamma^{2}x/2)$. The many-body
density can then be approximated by a simple sum of the single impurity
densities located at individual sites in an effective BEC matter-wave lattice
potential, $\rho(x)\sim\sum^{N}_{n=1}\tilde{\rho}(x-x_{n})$, where $x_{n}$ are
the individual center-of-mass positions of the localised impurities which are
separated by the average distance $|x_{n}-x_{n+1}|\sim L/N$ Keller  _et al._
(2022).

Figure 1:  (a) Density profiles and (b) single particle eigenstates of the
ground state of the impurities (red solid lines) pinned in a BEC with density
profile given by the blue dashed line at $t=0$. Evolution of the density of
the impurities, $\rho(x,t)$, after receiving a momentum kick of (c) $v=0.15$
and (f) $v=1$. The momentum kick is applied to the left and right localised
parts of the wavefunction
$\psi_{n}(x)\to\psi_{n}(x<0)e^{-ivx},\psi_{n}(x>0)e^{ivx}$ for $n=1,2$.
Evolution of the respective single particle states $|\phi_{1}(x,t)|^{2}$ and
$|\phi_{2}(x,t)|^{2}$ for (d-e) $v=0.15$ and (g-h) $v=1$. The parameters
chosen are $N_{\rm BEC}=5000$, $L_{\rm BEC}=25$, $L=15$, and $\gamma=2$.

To go beyond the stationary state and explore dynamics where individual
impurities overlap or scatter, it is necessary to solve the full Eqs. (1) and
(2) numerically. However, in the following we will show that in certain
regimes it is possible to describe the behavior observed by just considering
the soliton-train solutions that Eq. (4) allows for. For this we will confirm
that the wave profile is maintained during the free evolution, and that the
shape and phase difference are maintained after a collision Lamb (1980);
Martin  _et al._ (2007, 2008); Malomed (2014). Our first step is therefore to
solve Eqs. (1) and (2) for the minimal system of $N=2$ impurities and to study
the collisional properties by applying opposite momentum kicks $e^{\pm ivx}$
to the initially localised parts of the wavefunction in Fig. 1 (a). The
resulting density evolutions of the impurities for two different relative
initial velocities is shown in Figs. 1 (c) and (f). One can immediately see
that in both cases the impurities stay localized and maintain their individual
shape after the kick and also after the collision process. However, while for
the collision with large relative velocity, $v=1$, the two localised solitonic
impurities strongly overlap (and possibly cross), for the collision with lower
velocity, $v=0.15$, a clear repulsion between the impurities is visible. To
investigate the impurity scattering we calculate the relative distance

$d(t)=\int dx_{1}dx_{2}|x_{1}-x_{2}||\Phi(x_{1},x_{2},t)|^{2}\;,$ (5)

which is shown for different velocities in Fig. 2 (a). One can see that with
increasing relative velocities the impurities approach and overlap with each
other more and more, and their trajectories approach that of noninteracting
classical particles represented by the dotted line. We also show that the
results obtained from the effective description Eq. (4) (dashed lines) agree
quite well with the result from the full coupled system of Eq. (1) and Eq.
(2). This shows that we remain in the regime of validity of the Thomas–Fermi
approximation as the impurity velocity never exceeds the Landau critical
velocity, which can be estimated by the Bogoliubov speed of sound for a
homogeneous BEC to be $v_{\rm B}\approx\sqrt{N_{\rm BEC}/L_{\rm BEC}}\sim 14$.
We note, however, that immediately after the instantaneous momentum kick, a
weak burst of phononic excitations with velocity $\pm v_{\rm B}$ appears in
the BEC (see Fig. 2 (b)) Andrews  _et al._ (1997), nevertheless they do not
alter the dynamics of the impurities as their amplitude is small. The
effective nonlinear coupled equations in Eq. (4) therefore well describe the
system and its dynamics in the regime where phononic excitations in the BEC
can be neglected, which in turn means that the soliton-train solutions provide
an accurate description of the dynamics of the $N$ impurities 111It is
important to note that even though these soliton-trains are embedded into the
background of a BEC, they are inherently different to the well-studied dark-
bright soliton, e.g. Refs. Busch and Anglin (2001); Mistakidis  _et al._
(2019), since the density inhomogeneities in the BEC are small and do not
carry the phase defects that define dark solitons. They are also different to
the Bose–Fermi composite solitons investigated previously in Refs. Karpiuk
_et al._ (2004); Salerno (2005); Karpiuk  _et al._ (2006); Santhanam  _et al._
(2006); Adhikari and Salasnich (2007); Rakshit  _et al._ (2019a, b); DeSalvo
_et al._ (2019), where strong inter-species attraction surpasses the intra-
species repulsion among bosons resulting in the effective attraction in the
BEC leading to a formation of a bright soliton in the BEC supported by
fermions..

This agreement then allows us to understand the velocity dependence of the
scattering process from the fact that these solitonic impurities can be seen
as composite vector-solitons formed by the single particle states $\phi_{1}$
and $\phi_{2}$ of Eq. (4), which are coupled by their total density $\rho$.
While the scattering behaviour of scalar many-body Bose-condensed bright
solitons is purely determined by their relative phase Lamb (1980); Martin  _et
al._ (2007, 2008); Malomed (2014), for the soliton-trains formed by the single
particle states $\phi_{1}$ and $\phi_{2}$, one has to carefully examine the
eigenstates’ individual dynamics. For this one can see in Fig. 1 (b) that the
phase of state $\phi_{1}$ is flat, whereas the phase for state $\phi_{2}$
possess a $\pi$-phase jump at $x=0$. This jump leads to an effective repulsion
between the localized parts on the left and the right hand side of $\phi_{2}$,
which will dominate over the kinetic energy at small distances. For low
relative velocities the impurity density therefore never shows a crossing
Aitchison  _et al._ (1991); Malomed (1991); Stegeman and Segev (1999); Hegde
_et al._ (2022), as the coupling between $\phi_{1}$ and $\phi_{2}$ ensures
that the state $\phi_{1}$ follows the trajectory given by $\phi_{2}$ as can be
seen from Fig. 1 (d) and (e). This results in the repulsive collision dynamics
observed in Fig. 1 (c). On the other hand, for large relative velocities the
kinetic energy can be high enough to lead to a closer approach for $\phi_{2}$
and therefore a significant overlap in the density can become
indistinguishable from a crossing (see Figs. 1 (f-h)).

Figure 2:  (a) Average distance between impurities as given by Eq. (5), for
$v=0.2$ (blue), $0.6$ (red), $1$ (green), $3$ (yellow) calculated by solving
the full two-component system of Eqs. (1) and (2). The dashed lines are
calculated from the effective model Eq. (4). The dotted line is the
corresponding quantity for two noninteracting classical particles:
$d(t)=|d(0)-2vt|$. (b) Density evolution of the BEC after an impurity momentum
kick of $v=1$. The overlaid white dotted line corresponds to $x=\pm v_{\rm
B}t-d(0)/2$, where $v_{\rm B}$ is the Bogoliubov speed of sound for a
homogeneous BEC.  Figure 3:  (left column) Density evolution of four
impurities coupled to a BEC after a sudden quench from an external box
potential to a harmonic trap with trapping frequency $\omega=0.15$. (right
column) Time evolution of the largest occupation number, $\lambda_{0}(t)$, of
the RSPDM. The solid line is obtained by evolving the coupled equations (1)
and (2), while the dashed lines in panels (e) and (f) are obtained by evolving
Eq. (4). For reference, we show the coherence $\lambda_{0}$ of the ground
state of the TG gas in the harmonic trap with frequency $\omega$ as a dotted
line. The inter-species interaction strengths are $\gamma=0$, $\gamma=1$, and
$\gamma=1.8$ from top to bottom and the other parameter choices are $N_{\rm
BEC}=10^{4}$, $L_{\rm BEC}=50$, and $L=30$.

To explore the dynamical behaviour of larger trains of impurity solitons we
quench the system by instantaneously changing the external trap from a box to
a harmonic oscillator potential, $V_{\rm TG}(x)=\omega^{2}x^{2}/2$, with
$\omega$ in units of $1/t_{0}$. The resulting density dynamics for $N=4$
impurities is shown in Fig. 3 (a-c) for different coupling strengths: for
$\gamma=0$ the impurities are not coupled to the BEC and undergo free
evolution, whereas for $\gamma=1.8$ the system is deep in the pinned regime
where the impurities behave as soliton-trains. In between ($\gamma=1$) the
system is in an intermediate regime where the impurities are quasi-localized
but still have finite overlap with one another. In all cases, the time
evolution of the density profile shows the harmonic trap-induced breathing
mode with periodic collisions and revivals of the impurities every
$t\sim\pi/\omega$, which is a realization of the quantum Newton’s cradle
Kinoshita  _et al._ (2006); Caux  _et al._ (2019) in a coupled two component
system. Furthermore, one can note that in the case of $\gamma=1.8$, the
impurities are localized even after all four impurities have collided multiple
times providing further evidence of being the soliton-train.

Since the evolution Eqs. (1) and (2) can be used to simulate both the bosonic
TG impurities as well as free fermionic impurities, it is interesting to
explore how their dynamics are affected by differences in their quantum
statistics. However, since the Bose–Fermi mapping theorem ensures that the
density evolution is identical for both exchange symmetries, $\rho_{\rm
F}(x,t)=\rho_{\rm TG}(x,t)$, it can not be used to discern any differences
between them Girardeau (1960). On the other hand, the reduced single particle
density matrices (RSPDMs) and the (experimentally measurable) momentum
distributions can be different for bosonic or fermionic impurities. The RSPDM
of the impurities is given by
$\varrho(x,x^{\prime},t)=\int\Phi^{*}(x,x_{2},\dots,x_{N},t)\Phi(x^{\prime},x_{2},\dots,x_{N},t)dx_{2},\dots,dx_{N}$,
where $\Phi$ is the respective many-body wavefunction Girardeau  _et al._
(2001); Pezer and Buljan (2007). This matrix can be diagonalised to find the
so-called natural orbital basis,
$\varrho(x,x^{\prime},t)=\sum_{n}\lambda_{n}(t)\varphi^{*}_{n}(x,t)\varphi_{n}(x^{\prime},t)$,
where $\lambda_{n}(t)$ are the occupation probabilities of the orbitals
$\varphi_{n}(x,t)$. The largest occupation number $\lambda_{0}(t)$
characterizes the degree of coherence in the system and takes a maximal value
of $\lambda_{0}\approx N$ for a weakly interacting coherent BEC, and is
$\lambda_{0}=1$ for noninteracting fermions.

The coherence within the TG system strongly depends on the coupling strength
$\gamma$ (see Figs. 3 (d-f)). For $\gamma=0$ at time $t=0$, the TG impurities
overlap and strongly interact with each other, which reduces the coherence
compared to weakly interacting bosons, but still allows for a large degree of
coherence on the order of $\lambda_{0}\sim\sqrt{N}$ Girardeau  _et al._
(2001). Such a system may still be considered as a quasi-superfluid. When the
particles approach each other around the collision time, $t=\pi/2\omega$, the
coherence increases due to quasi-condensation Rigol and Muramatsu (2004).
Since for $\gamma=0$ the impurity system is integrable, the collisions and the
increase in coherence happen periodically. For $\gamma=1.8$ the pinned TG
state is initially fully incoherent due to the isolation of particles from
each other. However, during collisions, the coherence temporarily increases
rapidly as the impurities overlap, reaching about the same value as in the
uncoupled situation. The quench dynamics in the solitonic regime therefore
allows to realize a dynamical transition between different phases, with the
periodic oscillations between the pinned insulator and the quasi-superfluid
state being supported by the coupling to the BEC. Finally, in the intermediate
regime of $\gamma=1$, the impurities overlap throughout the dynamics, and
their coherence can be seen to decay after repeated collisions. This can be
attributed to the fact that the BEC-mediated coupling among the TG atoms
breaks the integrability of the system. We note that in all cases the maximum
velocity of the impurities is around $v\sim L\omega/2=2.25$, which is much
smaller than the Landau critical velocity estimated above. Thus, although we
do not show it here, the amplitudes of phonons excited in the BEC are
negligible, and the Thomas–Fermi approximation also holds in this case. The
dynamics observed from the effective nonlinear equation (dashed lines in Figs.
3 (e-f)) are therefore indistinguishable from the dynamics described by the
full system.

Figure 4: (top row) Momentum distribution $n(k,t)$ following a trap quench in
the free regime $\gamma=0$ for (a) TG and (b) Fermi impurities. (c) Peak of
the momentum distribution at $k=0$ for TG (solid line) and Fermi (dashed line)
impurities. (bottom row) Equivalent dynamics in the self-pinning regime
($\gamma=1.8$). The parameter choice is the same as in Fig. 3.

For a fermionic system the coherence is unchanged throughout the dynamics,
$\lambda_{n}(t)=1$ for $0\leq n\leq N-1$, due to the separable mean-field
ansatz between the noninteracting fermions and the condensate. However, the
dynamics of the fermionic momentum distribution $n(k,t)=\int
e^{-ik(x-x^{\prime})}\varrho(x,x^{\prime},t)dxdx^{\prime}$ can be compared to
the bosonic TG one, and distinct characteristics due to different particle
symmetries can be seen. The upper row of Fig. 4 shows the momentum
distributions for free evolution, $\gamma=0$, and one can see notable
differences. In particular, the TG impurities show a significant peak at $k=0$
whenever the atoms collide (see panel (c)). This is due to the impenetrable
hardcore character of their scattering, which reverses the impurity momenta
after each collision Atas  _et al._ (2017a, b). For fermions, this peak is
absent as they simply pass through each other. The lower row in Fig. 4 shows
the same quantities, but now for the situation where the impurities strongly
interact with the BEC, $\gamma=1.8$. One can see that the momentum
distribution of the TG and fermionic impurities are equivalent at all times
between collisions, highlighting the incoherent single particle nature of both
systems when deep in the pinned regime where we have $\lambda_{n}(t)=1$.
However, during the collisions the different scattering mechanisms are
enhanced by the BEC-mediated effective attraction. The $k=0$ mode for the
fermions is drastically depleted as the attractive interactions accelerate the
particles through one another when their densities overlap. Conversely, the
$k=0$ mode for TG impurities is increased as the particles scatter with larger
kinetic energy. Solitonic impurities in the BEC therefore act in different
ways according to their quantum statistics whenever the particles collide and
overlap, while at all other times they have an incoherent single particle
nature.

In conclusion, we have studied the dynamics of strongly correlated impurities,
namely bosons in the TG regime or fermions, which are coupled to a BEC in the
mean-field limit. Investigating the collision dynamics of two impurities in
the pinned limit, we provide strong evidence that the impurities satisfy the
characteristic scattering properties of soliton trains when the velocity of
the impurities is small enough such that the amplitude of the generated
phonons in the BEC are negligible. The system can then be accurately described
by a set of $N$-coupled nonlinear Schrödinger equations, which allows one to
identify the origin of the formation of the impurity brght soliton-train as a
consequence of exclusion principles and BEC-mediated self-attraction. Unlike
the soliton states formed by a (fully coherent) classical field such as in
BECs or nonlinear optical systems, the few-body nature of the impurity
soliton-train allows to study correlations among impurities and the quantum
statistical difference between TG bosons and noninteracting fermions. Indeed,
the coupling to the BEC can stabilize the impurities against dispersion,
allowing the localized pinned state to survive even after repeated collisions
and the induced breathing mode exhibiting dynamical transitions between
insulating and quasi-superfluid states. Moreover, TG bosons can be
distinguished from spinless fermions in the momentum distribution whenever the
impurities overlap, showing distinct collisional effects enhanced by the
presence of the BEC. An important extension would be the inclusion of
correlations between the BEC and impurities, allowing to explore beyond mean-
field effects on both the impurity localization and their soliton dynamics.
This would require describing the many-body state of the full coupled system
using techniques such as multi-configurational time-dependent Hartree methods
Streltsov  _et al._ (2011); Cao  _et al._ (2017), density-matrix
renormalization group Schollwöck (2005); Barbiero  _et al._ (2014) or exact
diagonalization Kanamoto  _et al._ (2005); Anh-Tai  _et al._ (2023), which we
leave as a future work.

###### Acknowledgements.

This work was supported by the Okinawa Institute of Science and Technology
Graduate University. The authors are grateful for the the Scientific Computing
and Data Analysis (SCDA) section of the Research Support Division at OIST. The
authors thank Tetsuro Nikuni and Tim Keller for insightful discussions. T.F.
acknowledges support from JSPS KAKENHI Grant No. JP23K03290. T.F. and T.B. are
also supported by JST Grant No. JPMJPF2221 and JSPS Bilateral Program No.
JPJSBP120227414.

## References

  * Lamb (1980) G. L. Lamb,  _Elements of soliton theory_ (Wiley New York, 1980).
  * Zakharov and Shabat (1972) V. E. Zakharov and A. B. Shabat, “Exact Theory of Two-dimensional Self-focusing and One-dimensional Self-modulation of Waves in Nonlinear Media,” Soviet Journal of Experimental and Theoretical Physics 34, 62 (1972).
  * Emplit  _et al._ (1987) P. Emplit, J.P. Hamaide, F. Reynaud, C. Froehly, and A. Barthelemy, “Picosecond steps and dark pulses through nonlinear single mode fibers,” Optics Communications 62, 374–379 (1987).
  * Krökel  _et al._ (1988) D. Krökel, N. J. Halas, G. Giuliani, and D. Grischkowsky, “Dark-pulse propagation in optical fibers,” Phys. Rev. Lett. 60, 29–32 (1988).
  * Kivshar and Malomed (1989) Yuri S. Kivshar and Boris A. Malomed, “Dynamics of solitons in nearly integrable systems,” Rev. Mod. Phys. 61, 763–915 (1989).
  * Kivshar (1993) Y.S. Kivshar, “Dark solitons in nonlinear optics,” IEEE Journal of Quantum Electronics 29, 250–264 (1993).
  * Malomed  _et al._ (2005) Boris A Malomed, Dumitru Mihalache, Frank Wise, and Lluis Torner, “Spatiotemporal optical solitons,” Journal of Optics B: Quantum and Semiclassical Optics 7, R53 (2005).
  * Burger  _et al._ (1999) S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G. V. Shlyapnikov, and M. Lewenstein, “Dark solitons in bose-einstein condensates,” Phys. Rev. Lett. 83, 5198–5201 (1999).
  * Busch and Anglin (2000) Th. Busch and J. R. Anglin, “Motion of Dark Solitons in Trapped Bose-Einstein Condensates,” Phys. Rev. Lett. 84, 2298–2301 (2000).
  * Carr  _et al._ (2000a) L. D. Carr, Charles W. Clark, and W. P. Reinhardt, “Stationary solutions of the one-dimensional nonlinear schrödinger equation. i. case of repulsive nonlinearity,” Phys. Rev. A 62, 063610 (2000a).
  * Carr  _et al._ (2000b) L. D. Carr, Charles W. Clark, and W. P. Reinhardt, “Stationary solutions of the one-dimensional nonlinear schrödinger equation. ii. case of attractive nonlinearity,” Phys. Rev. A 62, 063611 (2000b).
  * Denschlag  _et al._ (2000) J. Denschlag, J. E. Simsarian, D. L. Feder, Charles W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, S. L. Rolston, B. I. Schneider, and W. D. Phillips, “Generating solitons by phase engineering of a bose-einstein condensate,” Science 287, 97–101 (2000).
  * Khaykovich  _et al._ (2002) L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, “Formation of a matter-wave bright soliton,” Science 296, 1290–1293 (2002).
  * Strecker  _et al._ (2002) Kevin E. Strecker, Guthrie B. Partridge, Andrew G. Truscott, and Randall G. Hulet, “Formation and propagation of matter-wave soliton trains,” Nature 417, 150–153 (2002).
  * Marchant  _et al._ (2013) A. L. Marchant, T. P. Billam, T. P. Wiles, M. M. H. Yu, S. A. Gardiner, and S. L. Cornish, “Controlled formation and reflection of a bright solitary matter-wave,” Nature Communications 4, 1865 (2013).
  * Kartashov  _et al._ (2019) Yaroslav V. Kartashov, Gregory E. Astrakharchik, Boris A. Malomed, and Lluis Torner, “Frontiers in multidimensional self-trapping of nonlinear fields and matter,” Nature Reviews Physics 1, 185–197 (2019).
  * Hoefer  _et al._ (2011) M. A. Hoefer, J. J. Chang, C. Hamner, and P. Engels, “Dark-dark solitons and modulational instability in miscible two-component bose-einstein condensates,” Phys. Rev. A 84, 041605 (2011).
  * Yan  _et al._ (2012) D Yan, J J Chang, C Hamner, M Hoefer, P G Kevrekidis, P Engels, V Achilleos, D J Frantzeskakis, and J Cuevas, “Beating dark–dark solitons in bose–einstein condensates,” Journal of Physics B: Atomic, Molecular and Optical Physics 45, 115301 (2012).
  * Busch and Anglin (2001) Th. Busch and J. R. Anglin, “Dark-bright solitons in inhomogeneous bose-einstein condensates,” Phys. Rev. Lett. 87, 010401 (2001).
  * Nistazakis  _et al._ (2008) H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and R. Carretero-González, “Bright-dark soliton complexes in spinor bose-einstein condensates,” Phys. Rev. A 77, 033612 (2008).
  * Becker  _et al._ (2008) Christoph Becker, Simon Stellmer, Parvis Soltan-Panahi, Sören Dörscher, Mathis Baumert, Eva-Maria Richter, Jochen Kronjäger, Kai Bongs, and Klaus Sengstock, “Oscillations and interactions of dark and dark–bright solitons in bose–einstein condensates,” Nature Physics 4, 496–501 (2008).
  * Hamner  _et al._ (2011) C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer, “Generation of dark-bright soliton trains in superfluid-superfluid counterflow,” Phys. Rev. Lett. 106, 065302 (2011).
  * Bakkali-Hassani  _et al._ (2021) B. Bakkali-Hassani, C. Maury, Y.-Q. Zou, É. Le Cerf, R. Saint-Jalm, P. C. M. Castilho, S. Nascimbene, J. Dalibard, and J. Beugnon, “Realization of a townes soliton in a two-component planar bose gas,” Phys. Rev. Lett. 127, 023603 (2021).
  * Edmonds  _et al._ (2018) M. J. Edmonds, T. P. Billam, S. A. Gardiner, and Th. Busch, “Noise-free generation of bright matter-wave solitons,” Phys. Rev. A 98, 063626 (2018).
  * Carr and Brand (2004) L. D. Carr and J. Brand, “Spontaneous soliton formation and modulational instability in bose-einstein condensates,” Phys. Rev. Lett. 92, 040401 (2004).
  * Mollenauer  _et al._ (1990) L. F. Mollenauer, M. J. Neubelt, S. G. Evangelides, J. P. Gordon, J. R. Simpson, and L. G. Cohen, “Experimental study of soliton transmission over more than 10,000 km in dispersion-shifted fiber,” Opt. Lett. 15, 1203–1205 (1990).
  * Mamyshev  _et al._ (1991) P.V. Mamyshev, S.V. Chernikov, and E.M. Dianov, “Generation of fundamental soliton trains for high-bit-rate optical fiber communication lines,” IEEE Journal of Quantum Electronics 27, 2347–2355 (1991).
  * Haus and Wong (1996) Hermann A. Haus and William S. Wong, “Solitons in optical communications,” Rev. Mod. Phys. 68, 423–444 (1996).
  * Kivshar and Luther-Davies (1998) Yuri S. Kivshar and Barry Luther-Davies, “Dark optical solitons: physics and applications,” Physics Reports 298, 81–197 (1998).
  * Buryak  _et al._ (2002) Alexander V. Buryak, Paolo Di Trapani, Dmitry V. Skryabin, and Stefano Trillo, “Optical solitons due to quadratic nonlinearities: from basic physics to futuristic applications,” Physics Reports 370, 63–235 (2002).
  * Chen  _et al._ (2012) Zhigang Chen, Mordechai Segev, and Demetrios N Christodoulides, “Optical spatial solitons: historical overview and recent advances,” Reports on Progress in Physics 75, 086401 (2012).
  * Ramanathan  _et al._ (2011) A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K. Campbell, “Superflow in a toroidal bose-einstein condensate: An atom circuit with a tunable weak link,” Phys. Rev. Lett. 106, 130401 (2011).
  * Martin and Ruostekoski (2012) A D Martin and J Ruostekoski, “Quantum dynamics of atomic bright solitons under splitting and recollision, and implications for interferometry,” New Journal of Physics 14, 043040 (2012).
  * Polo and Ahufinger (2013) J. Polo and V. Ahufinger, “Soliton-based matter-wave interferometer,” Phys. Rev. A 88, 053628 (2013).
  * Wales  _et al._ (2020) Oliver J. Wales, Ana Rakonjac, Thomas P. Billam, John L. Helm, Simon A. Gardiner, and Simon L. Cornish, “Splitting and recombination of bright-solitary-matter waves,” Communications Physics 3, 51 (2020).
  * Szigeti  _et al._ (2020) Stuart S. Szigeti, Samuel P. Nolan, John D. Close, and Simon A. Haine, “High-precision quantum-enhanced gravimetry with a bose-einstein condensate,” Phys. Rev. Lett. 125, 100402 (2020).
  * Keller  _et al._ (2022) Tim Keller, Thomás Fogarty, and Thomas Busch, “Self-pinning transition of a tonks-girardeau gas in a bose-einstein condensate,” Phys. Rev. Lett. 128, 053401 (2022).
  * Bruderer  _et al._ (2008) M. Bruderer, W. Bao, and D. Jaksch, “Self-trapping of impurities in bose-einstein condensates: Strong attractive and repulsive coupling,” Europhysics Letters 82, 30004 (2008).
  * Girardeau (1960) M. Girardeau, “Relationship between systems of impenetrable bosons and fermions in one dimension,” Journal of Mathematical Physics 1, 516–523 (1960).
  * Pitaevskii and Stringari (2016) Lev Pitaevskii and Sandro Stringari,  _Bose-Einstein condensation and superfluidity_ , Vol. 164 (Oxford University Press, 2016).
  * Miyakawa  _et al._ (2004) Takahiko Miyakawa, Hiroyuki Yabu, and Toru Suzuki, “Peierls instability, periodic bose-einstein condensates, and density waves in quasi-one-dimensional boson-fermion mixtures of atomic gases,” Phys. Rev. A 70, 013612 (2004).
  * Stringari (1996) S. Stringari, “Collective excitations of a trapped bose-condensed gas,” Phys. Rev. Lett. 77, 2360–2363 (1996).
  * Manakov (1973) S. V. Manakov, “On the theory of two-dimensional stationary self-focusing of electromagnetic waves,” Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki 65, 505–516 (1973).
  * Gordon (1983) J. P. Gordon, “Interaction forces among solitons in optical fibers,” Opt. Lett. 8, 596–598 (1983).
  * Radhakrishnan  _et al._ (1997) R. Radhakrishnan, M. Lakshmanan, and J. Hietarinta, “Inelastic collision and switching of coupled bright solitons in optical fibers,” Phys. Rev. E 56, 2213–2216 (1997).
  * Kanna  _et al._ (2006) T. Kanna, M. Lakshmanan, P. Tchofo Dinda, and Nail Akhmediev, “Soliton collisions with shape change by intensity redistribution in mixed coupled nonlinear schrödinger equations,” Phys. Rev. E 73, 026604 (2006).
  * Martin  _et al._ (2007) A. D. Martin, C. S. Adams, and S. A. Gardiner, “Bright matter-wave soliton collisions in a harmonic trap: Regular and chaotic dynamics,” Phys. Rev. Lett. 98, 020402 (2007).
  * Martin  _et al._ (2008) A. D. Martin, C. S. Adams, and S. A. Gardiner, “Bright solitary-matter-wave collisions in a harmonic trap: Regimes of solitonlike behavior,” Phys. Rev. A 77, 013620 (2008).
  * Malomed (2014) B.A. Malomed,  _Spontaneous Symmetry Breaking, Self-Trapping, and Josephson Oscillations_, Progress in Optical Science and Photonics (Springer Berlin Heidelberg, 2014).
  * Andrews  _et al._ (1997) M. R. Andrews, D. M. Kurn, H.-J. Miesner, D. S. Durfee, C. G. Townsend, S. Inouye, and W. Ketterle, “Propagation of sound in a bose-einstein condensate,” Phys. Rev. Lett. 79, 553–556 (1997).
  * Note (1) It is important to note that even though these soliton-trains are embedded into the background of a BEC, they are inherently different to the well-studied dark-bright soliton, e.g. Refs. Busch and Anglin (2001); Mistakidis  _et al._ (2019), since the density inhomogeneities in the BEC are small and do not carry the phase defects that define dark solitons. They are also different to the Bose–Fermi composite solitons investigated previously in Refs. Karpiuk  _et al._ (2004); Salerno (2005); Karpiuk  _et al._ (2006); Santhanam  _et al._ (2006); Adhikari and Salasnich (2007); Rakshit  _et al._ (2019a, b); DeSalvo  _et al._ (2019), where strong inter-species attraction surpasses the intra-species repulsion among bosons resulting in the effective attraction in the BEC leading to a formation of a bright soliton in the BEC supported by fermions.
  * Aitchison  _et al._ (1991) J. S. Aitchison, A. M. Weiner, Y. Silberberg, D. E. Leaird, M. K. Oliver, J. L. Jackel, and P. W. E. Smith, “Experimental observation of spatial soliton interactions,” Opt. Lett. 16, 15–17 (1991).
  * Malomed (1991) Boris A. Malomed, “Bound solitons in the nonlinear schrödinger–ginzburg-landau equation,” Phys. Rev. A 44, 6954–6957 (1991).
  * Stegeman and Segev (1999) George I. Stegeman and Mordechai Segev, “Optical spatial solitons and their interactions: Universality and diversity,” Science 286, 1518–1523 (1999).
  * Hegde  _et al._ (2022) Gautam Hegde, Sandra M. Jose, and Rejish Nath, “Dynamics of a pair of overlapping polar bright solitons in spin-1 bose-einstein condensates,” Phys. Rev. A 106, 043307 (2022).
  * Kinoshita  _et al._ (2006) Toshiya Kinoshita, Trevor Wenger, and David S. Weiss, “A quantum newton’s cradle,” Nature 440, 900–903 (2006).
  * Caux  _et al._ (2019) Jean-Sébastien Caux, Benjamin Doyon, Jérôme Dubail, Robert Konik, and Takato Yoshimura, “Hydrodynamics of the interacting Bose gas in the Quantum Newton Cradle setup,” SciPost Phys. 6, 070 (2019).
  * Girardeau  _et al._ (2001) M. D. Girardeau, E. M. Wright, and J. M. Triscari, “Ground-state properties of a one-dimensional system of hard-core bosons in a harmonic trap,” Phys. Rev. A 63, 033601 (2001).
  * Pezer and Buljan (2007) R. Pezer and H. Buljan, “Momentum distribution dynamics of a tonks-girardeau gas: Bragg reflections of a quantum many-body wave packet,” Phys. Rev. Lett. 98, 240403 (2007).
  * Rigol and Muramatsu (2004) Marcos Rigol and Alejandro Muramatsu, “Emergence of quasicondensates of hard-core bosons at finite momentum,” Phys. Rev. Lett. 93, 230404 (2004).
  * Atas  _et al._ (2017a) Y. Y. Atas, I. Bouchoule, D. M. Gangardt, and K. V. Kheruntsyan, “Collective many-body bounce in the breathing-mode oscillations of a tonks-girardeau gas,” Phys. Rev. A 96, 041605 (2017a).
  * Atas  _et al._ (2017b) Y. Y. Atas, D. M. Gangardt, I. Bouchoule, and K. V. Kheruntsyan, “Exact nonequilibrium dynamics of finite-temperature tonks-girardeau gases,” Phys. Rev. A 95, 043622 (2017b).
  * Streltsov  _et al._ (2011) Alexej I. Streltsov, Ofir E. Alon, and Lorenz S. Cederbaum, “Swift loss of coherence of soliton trains in attractive bose-einstein condensates,” Phys. Rev. Lett. 106, 240401 (2011).
  * Cao  _et al._ (2017) L. Cao, V. Bolsinger, S. I. Mistakidis, G. M. Koutentakis, S. Krönke, J. M. Schurer, and P. Schmelcher, “A unified ab initio approach to the correlated quantum dynamics of ultracold fermionic and bosonic mixtures,” The Journal of Chemical Physics 147, 044106 (2017), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4993512/15529588/044106_1_online.pdf .
  * Schollwöck (2005) U. Schollwöck, “The density-matrix renormalization group,” Rev. Mod. Phys. 77, 259–315 (2005).
  * Barbiero  _et al._ (2014) L. Barbiero, B. A. Malomed, and L. Salasnich, “Quantum bright solitons in the bose-hubbard model with site-dependent repulsive interactions,” Phys. Rev. A 90, 063611 (2014).
  * Kanamoto  _et al._ (2005) Rina Kanamoto, Hiroki Saito, and Masahito Ueda, “Symmetry breaking and enhanced condensate fraction in a matter-wave bright soliton,” Phys. Rev. Lett. 94, 090404 (2005).
  * Anh-Tai  _et al._ (2023) Tran Duong Anh-Tai, Mathias Mikkelsen, Thomas Busch, and Thomás Fogarty, “Quantum chaos in interacting Bose-Bose mixtures,” SciPost Phys. 15, 048 (2023).
  * Mistakidis  _et al._ (2019) S. I. Mistakidis, L. Hilbig, and P. Schmelcher, “Correlated quantum dynamics of two quenched fermionic impurities immersed in a bose-einstein condensate,” Phys. Rev. A 100, 023620 (2019).
  * Karpiuk  _et al._ (2004) T. Karpiuk, M. Brewczyk, S. Ospelkaus-Schwarzer, K. Bongs, M. Gajda, and K. Rza¸żewski, “Soliton trains in bose-fermi mixtures,” Phys. Rev. Lett. 93, 100401 (2004).
  * Salerno (2005) Mario Salerno, “Matter-wave quantum dots and antidots in ultracold atomic bose-fermi mixtures,” Phys. Rev. A 72, 063602 (2005).
  * Karpiuk  _et al._ (2006) Tomasz Karpiuk, Mirosław Brewczyk, and Kazimierz Rza¸żewski, “Bright solitons in bose-fermi mixtures,” Phys. Rev. A 73, 053602 (2006).
  * Santhanam  _et al._ (2006) J. Santhanam, V. M. Kenkre, and V. V. Konotop, “Solitons of bose-fermi mixtures in a strongly elongated trap,” Phys. Rev. A 73, 013612 (2006).
  * Adhikari and Salasnich (2007) Sadhan K. Adhikari and Luca Salasnich, “One-dimensional superfluid bose-fermi mixture: Mixing, demixing, and bright solitons,” Phys. Rev. A 76, 023612 (2007).
  * Rakshit  _et al._ (2019a) Debraj Rakshit, Tomasz Karpiuk, Paweł Zin, Mirosław Brewczyk, Maciej Lewenstein, and Mariusz Gajda, “Self-bound bose–fermi liquids in lower dimensions,” New Journal of Physics 21, 073027 (2019a).
  * Rakshit  _et al._ (2019b) Debraj Rakshit, Tomasz Karpiuk, Mirosław Brewczyk, and Mariusz Gajda, “Quantum Bose-Fermi droplets,” SciPost Phys. 6, 079 (2019b).
  * DeSalvo  _et al._ (2019) B. J. DeSalvo, Krutik Patel, Geyue Cai, and Cheng Chin, “Observation of fermion-mediated interactions between bosonic atoms,” Nature 568, 61–64 (2019).