# Searching for universality of dineutron correlation at the surface of
Borromean nuclei

A. Corsi<EMAIL_ADDRESS>Y. Kubota  J. Casal  M. Gómez-Ramos  A. M. Moro  G.
Authelet  H. Baba  C. Caesar  D. Calvet  A. Delbart  M. Dozono  J. Feng  F.
Flavigny  J.-M. Gheller  J. Gibelin  A. Giganon  A. Gillibert  K. Hasegawa  T.
Isobe  Y. Kanaya  S. Kawakami  D. Kim  Y. Kiyokawa  M. Kobayashi  N. Kobayashi
T. Kobayashi  Y. Kondo  Z. Korkulu  S. Koyama  V. Lapoux  Y. Maeda  F. M.
Marqués  T. Motobayashi  T. Miyazaki  T. Nakamura  N. Nakatsuka  Y. Nishio  A.
Obertelli  A. Ohkura  N. A. Orr  S. Ota  H. Otsu  T. Ozaki  V. Panin  S.
Paschalis  E. C. Pollacco  S. Reichert  J.-Y. Rousse  A. T. Saito  S.
Sakaguchi  M. Sako  C. Santamaria  M. Sasano  H. Sato  M. Shikata  Y. Shimizu
Y. Shindo  L. Stuhl  T. Sumikama  Y.L. Sun  M. Tabata  Y. Togano  J. Tsubota
T. Uesaka  Z. H. Yang  J. Yasuda  K. Yoneda  J. Zenihiro  Département de
Physique Nucléaire, IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-
Yvette, France  RIKEN Nishina Center, Hirosawa 2-1, Wako, Saitama 351-0198,
Japan  Center for Nuclear Study, University of Tokyo, Hongo 7-3-1, Bunkyo,
Tokyo 113-0033, Japan  European Centre for Theoretical Studies in Nuclear
Physics and Related Areas (ECT∗), Villa Tambosi, Strada delle Tabarelle 286,
I-38123 Trento, Italy  Dipartimento di Fisica e Astronomia “G. Galilei” and
INFN - Sezione di Padova, Via Marzolo 8, I-35131 Padova, Italy  Departamento
de Física Atómica, Molecular y Nuclear, Facultad de Física, Universidad de
Sevilla, Apartado 1065, E-41080 Sevilla, Spain  Department of Physics,
Technische Universitat Darmstadt  Department of Physics, Peking University
Institut de Physique Nucleaire Orsay, IN2P3-CNRS, F-91406 Orsay Cedex, France
LPC Caen, ENSICAEN, Universite de Caen, CNRS/IN2P3, F-14050 Caen, France
Department of Physics, Tohoku University, Aramaki Aza-Aoba 6-3, Aoba, Sendai,
Miyagi 980-8578, Japan  Department of Applied Physics, University of Miyazaki,
Gakuen-Kibanadai-Nishi 1-1, Miyazaki 889-2192, Japan  Department of Physics,
Ehwa Womans University  Department of Physics, University of Tokyo, Hongo
7-3-1, Bunkyo, Tokyo 113-0033, Japan  Department of Physics, Tokyo Institute
of Technology, 2-12-1 O-Okayama, Meguro, Tokyo 152-8551, Japan  MTA Atomki,
P.O. Box 51, Debrecen H-4001, Hungary  Department of Physics, Kyoto
University, Kitashirakawa, Sakyo, Kyoto 606-8502, Japan  Department of
Physics, Kyushu University, Nishi, Fukuoka 819-0395, Japan  Research Center
for Nuclear Physics, Osaka University, 10-1 Mihogaoka, Ibaraki, Osaka
567-0047, Japan  Department of Physics, Technische Universitat Munchen  Center
for Exotic Nuclear Studies, Institute for Basic Science, Daejeon 34126,
Republic of Korea  School of Physics and State Key Laboratory of Nuclear
Physics and Technology, Peking University, Beijing 100871, China

###### Abstract

The dineutron correlation is systematically studied in three different
Borromean nuclei near the neutron dripline, 11Li, 14Be and 17B, via the
$(p,pn)$ knockout reaction measured at the RIBF facility in RIKEN. For the
three nuclei, the correlation angle between the valence neutrons is found to
be largest in the same range of intrinsic momenta, which can be associated to
the nuclear surface. This result reinforces the prediction that the formation
of the dineutron is universal in environments with low neutron density, such
as the surface of neutron-rich Borromean nuclei.

###### keywords:

quasi-free scattering , Borromean nuclei , three-body model , Jacobi
coordinates , dineutron

††journal: Physics Letters B

##  1 Introduction

Halo nuclei appear close to the neutron dripline and present a diffuse matter
distribution due to the reduced binding energy of the valence neutrons [1, 2].
Nuclei formed by a core and two loosely bound neutrons, such that the
subsystem formed by the core and the neutron is unbound, are called Borromean.
Most of them present nuclear halos near the neutron dripline. Some examples
are 6He, 11Li, 14Be and 17,19B. The correlation between the neutrons plays an
essential role to stabilize these nuclei and has been the subject of a number
of studies [3, 4, 5]. We discuss here a specific form of spatially localized
pairing correlation called dineutron [6]. The strength of the pairing
correlation evolves with density, going from the BCS
(Bardeen–Cooper–Schrieffer) regime of loosely spaced correlations to the BEC
(Bose-Einstein Condensate) regime of compact space correlation with decreasing
density. This regime is expected to appear at the surface of neutron-rich
nuclei, where neutron density is 10-4 to 0.5 of the saturation density [7, 8].
Its onset appears to be strongly linked to the admixture of different parities
in the wavefunction describing the valence neutrons [5]. Halo nuclei, with
their diffuse matter distribution, are an ideal probe to study this low-
density correlation. The dineutron was experimentally revealed in 6He [9, 10]
and 11Li [11, 12, 13, 14], however the experimental evidence is still scarce.
Typically, dineutron correlation is explored via the opening angle between the
two neutrons [11, 15]. An opening angle below 90∘ (90∘ corresponding to the
non-correlated case) in coordinate space or above 90∘ in momentum space [16,
17] points to a strong spatial correlation yielding a compact configuration.
Intuitively, the opening angle in coordinate space is related to the inter-
nucleon distance. References [11] measured the $E1$ strength after Coulomb
dissociation of 11Li, and extracted the opening angle in coordinate space
based on the cluster sum rule and assuming an inert core [18]. For 11Li the
average angle obtained was $\theta$ = 48${}^{+14}_{-18}$ ∘, corresponding to a
strong dineutron correlation. A more refined estimation can be obtained if the
average neutron-neutron distance is measured independently. In such a way the
authors of Ref. [19] deduced from the $B(E1)$ measurement of Ref. [11] a value
of the mean opening angle of $\sim$ 56.2∘. The average neutron-neutron
separation was estimated via measurements of the two-neutron correlation
function in dissociation reactions. This method has been applied to 6He, 11Li,
14Be [20, 21]. A combined analysis of the $B(E1)$ measurement [11] and the
correlation function [21] has given a somewhat larger value of the opening
angle in 11Li, $\sim 66^{\circ}$ [15], corresponding to a reduced dineutron
correlation. Two-neutron transfer reactions can also be used to study $nn$
correlations [4]. In Ref. [22], it was shown that the description of 11Li data
required a model with a large pairing correlation.

Nucleon removal reactions are another method to access the opening angle [23,
12, 21, 13]. The authors of Ref. [13] measured an opening angle in momentum
space of $103.4\pm 2.1^{\circ}$ for 11Li, and suggest that a dineutron
configuration exists also in 14Be, although less developed than in 11Li. In
both nuclei, the dineutron appears due to the mixing of different-parity
orbitals [5, 12]. In contrast, the structure of the halo in 17B is mainly of
$d$-wave character with small $s$-wave admixture [24], which should hinder the
development of the dineutron correlation as both orbitals have the same
parity. We note that the study of these correlations via breakup and knockout
reactions is complementary to that performed through 2$n$ decay [25, 26, 27],
which has been used to explore the properties of two-neutron unbound systems.

Recently, Ref. [14] introduced a new method based on quasi-free scattering
reactions to study dineutron correlation as a function of its peripherality,
i.e., distance from the baricenter of the system, and applied it to 11Li. The
observable related to the peripherality is the intrinsic momentum of the
removed nucleon in quasi-free scattering reactions. Being a fast removal
process, the impact of final-state interactions on the observable of interest
is assumed to be reduced, making its interpretation more straight-forward. In
this work we search for dineutron correlations applying this same method to
14Be and 17B, measured in the same experiment as Ref. [14]. Our goal is to
assess whether such correlation appears as a general feature at the surface of
neutron-rich Borromean nuclei. The data is compared to calculations using a
three-body model for the projectile and a quasi-free sudden model to describe
the knockout process [28].

##  2  Experimental results

###  2.1 Setup

The experiment was performed at the Radioactive Isotope Beam Factory operated
by the RIKEN Nishina Center and the Center for Nuclear Study (CNS) of the
University of Tokyo. Secondary beams were produced using projectile
fragmentation of a 48Ca primary beam at 345 MeV/nucleon with a typical
intensity of 400 particle nA on a Be target. Fragmentation products were
separated, detected and identified via the BigRIPS fragment separator [29].
The cocktail beam was composed by 11Li, 14Be, and 17B, with a percentage of
$\sim$ 80%, 12%, and 8%, respectively. It impinged on the secondary target
with an average energy of 246, 265 and 277 MeV/nucleon, respectively. The
secondary target was the 15-cm thick liquid hydrogen target from the MINOS
device [30], and was surrounded by a Time Projection Chamber acting as vertex
tracker together with the beam tracking MWDC detectors. The detection system
included the WINDS array of plastic scintillators for knockout neutron
detection, and a MWDC followed by an array of plastic scintillators for the
recoil proton detection. Those two detectors were key for the measurement of
the intrinsic momentum of the removed neutron and the opening angle in the
$(p,pn)$ reaction. The standard SAMURAI setup consisting of a set of drift
chambers, the SAMURAI dipole magnet and two hodoscope walls was used for
fragment analysis [31]. The neutrons emitted at forward angles were detected
by the NEBULA plastic scintillator array [32]. We evaluated the acceptance cut
they induce on the measurement of the intrinsic momentum and opening angle
distribution using a Geant4 simulation. No bias is introduced by the
experimental setup on the opening angle distribution, while the acceptance
decreases for increasing intrinsic momentum (leading to off-plane scattering),
as shown in the inset of Fig. 1. More details on the rest of the setup can be
found in Ref. [33, 14] and references therein.

###  2.2  Dineutron correlation

The measurement of the momenta of the outgoing proton and removed neutron
allows to reconstruct the intrinsic momentum of the neutron before removal
(within the quasi-free approximation):

$\vec{k_{y}}:=\vec{k_{n1}}=\vec{k_{n1}^{\prime}}+\vec{k_{p}^{\prime}}-\vec{k_{p}}\\\
$ (1)

where $k_{n1}$ ($k_{n1}^{{}^{\prime}}$) is the momentum of the neutron in the
initial (final) state and $k_{p}$ ($k_{p}^{{}^{\prime}}$) the one of the
target (recoil) proton. The correlation angle, or the opening angle $\theta$
in momentum space, is the angle between the Jacobi momenta $k_{x}$ and
$k_{y}$:

$cos(\theta)=\frac{\vec{k_{x}}\cdot\vec{k_{y}}}{|\vec{k_{x}}||\vec{k_{y}}|}$
(2)

with

$\vec{k_{x}}=\vec{k^{\prime}_{n2}}-\vec{k^{\prime}_{f}}$ (3)

where $k^{\prime}_{n2}$, $k^{\prime}_{f}$ are the momenta of the remaining
valence neutron and fragment in the final state. This representation of the
three-body system in terms of Y Jacobi coordinates is illustrated in the inset
of Fig. 2. In the following, we illustrate the intrinsic momentum and
correlation angle distribution for the case of 14Be. Both are compared with a
theoretical calculation performed within a quasi-free sudden model [28, 34],
using the three-body model for 14Be from [33].

Figure 1 shows the intrinsic momentum distribution of the removed nucleon for
two different relative-energy intervals in the 13Be system. The theoretical
distributions are already corrected for the experimental acceptance (see inset
of Fig. 1a) and convoluted with the experimental resolution of 0.17 fm-1. Each
relative-energy (${}^{12}\mathrm{Be}+n$) interval encompasses a peak in the
spectrum of 13Be [33]. The comparison to theoretical calculations show that
the 0-1.5 MeV interval is dominated by the $p_{1/2}$ component (72% of the
total in this energy range), while the 1.5-3 MeV interval is dominated by the
$d_{5/2}$ component (60% of the total). This is consistent with the
interpretation of the relative-energy spectrum of 13Be provided in Ref. [35,
33] as composed of a $p$-wave resonance centered at 0.5 MeV followed by a
broader $d$-wave resonance. The different lines in Fig. 1 are labeled as
$J^{\pi}[\ell_{j}\otimes S_{c}]$, where the single-particle angular momentum
$\ell_{j}$ couples with the spin of the core $S_{c}$ to give the total angular
momentum $J^{\pi}$ of the binary subsystem 13Be after knockout. Note that,
since the ground state of 14Be is a 0+ state, the angular momentum of the
knocked-out neutron has to match $J^{\pi}$, e.g., $5/2^{+}$ contributions
correspond to a $d$-wave knockout. It is worth noting that the calculations
presented in Fig. 1 are not a fit to the experimental data but the results of
the structure model (and corresponding partial-wave content) of Ref. [33],
therefore the agreement is not perfect. In particular the disagreement in the
peak in the lower energy range may suggest a larger $s$-wave component.
However, in [33], an increase in $s$-wave led to a worse description of the
low-energy distribution. Similarly, there is a slight disagreement for the
largest $k_{y}$ values that may be associated to missing components in the
wave function due to limitations of the model, as discussed in [33] for large
relative energies.

Figure 1:  Comparison between experimental and theoretical intrinsic momentum
distribution for the knocked-out neutron of 14Be for (a) $0<E_{\rm rel}<1.5$
MeV and (b) $1.5<E_{\rm rel}<3.0$ MeV in the 13Be system. The dashed lines
correspond to the contribution from different 13Be states (see the text for
details), and the solid line is the total. The inset shows the acceptance.

Fig. 2 shows the correlation angle distribution integrated over all intrinsic
momenta, and for intrinsic momenta between 0.2 fm-1 and 0.4 fm-1. One can see
that the inclusive distribution is rather symmetric, while an asymmetry
appears for some range of values of the intrinsic momentum. The calculations
are able to capture this behaviour. The range between 0.2 fm-1 and 0.4 fm-1 is
the one yielding the maximum asymmetry with an enhancement of values of the
correlation angle larger than 90∘. This points towards a geometrically compact
configuration of the two-neutron system (the dineutron) at low intrinsic
momenta, which can be associated to the nuclear surface.

Figure 2:  Correlation angle distribution for 14Be (integrated in
$E_{\mathrm{rel}}$), compared with our theoretical calculation. The lower
points and dashed line correspond to an intrinsic momentum window of
$0.2<k_{y}<0.4$ fm-1. A schematic representation of the angle in Jacobi
coordinates is embedded in the figure.

The average correlation angle $\theta$, obtained taking event by event the
arccos of the data plotted in Fig. 2, is plotted as a function of the
intrinsic momentum ranging from 0 to 1.8 fm-1 in Fig. 3. The data for 14Be are
compared to the ones for 11Li and 17B measured in the same experiment [14, 33,
24]. The nucleus of 11Li is considered as a reference case of well developed
dineutron correlation [11, 13, 14] so, as expected, it presents the largest
deviation from 90∘. It is however remarkable that for both 14Be and 17B the
data also show a significant deviation in the correlation angle distribution
in the same range of momenta. This deviation points to the appearance of a
dineutron correlation for intrinsic momenta smaller than 0.4 fm-1, which
corresponds to the nuclear surface [14], for all measured nuclei. We note that
the larger value of the correlation angle occurring around 0.2 fm-1 is clearly
above 90∘, even taking into account the errors. This constitutes the first
experimental evidence supporting universality of the dineutron correlation in
the low-density nuclear surface of Borromean nuclei, which had been previously
suggested [7, 8, 14]. It is worth noting that, depending on the probe, an
inclusive measurement of the correlation angle will be sensitive to a rather
large region of the nucleus (including the interior), and the dineutron
correlation signal may be damped, as shown in Fig. 2.

Figure 3:  Average correlation angle as a function the intrinsic momentum
$\vec{k}_{y}$ for $(p,pn)$ reaction on 11Li, 14Be and 17B. The points are
presented with their statistical error (bars) and systematic error (band).
Points for 14Be and 17B are slighty offsetted along the x-axis for clarity.

##  3 Theoretical analysis

To quantify and understand the mechanism behind the onset of the dineutron
correlation, we compare the experimental result of Fig. 3 to theoretical
predictions in Fig. 4. The theoretical description combines three-body models
within the hyperspherical framework to describe the structure of Borromean
nuclei [36, 37] and an eikonal sudden description of the $(p,pN)$ reaction
[34]. This description has been used to study 11Li$(p,pn)$ in [28]. We refer
the reader to [28] for a detailed overview of the model. The structure model
for 11Li was originally introduced in Ref. [38] to describe GSI $(p,pn)$ data
[39], and later revisited in [28] to incorporate $d$-wave contributions. This
model leads to a large admixture between $s$ and $p$ waves ($\sim$ 60% and
30%, respectively), and the computed core-$nn$ rms distance is 4.9 fm, which
compares well with the experimental value derived from Coulomb Dissociation
data [11]. For 14Be, we adopt the model in Ref. [33], which is dominated by a
low-lying $p$-wave resonance in 13Be ($\sim$ 60% of the wave function comes
from $p$ waves) and includes the effect of the first 2+ excited state of the
12Be core (which amounts to roughly 20% of the norm of the ground state of
14Be). For 17B, the three-body wave function was computed by fixing a simple
model neglecting the spin of the core, in the same spirit as the 19B
calculations in Ref. [40], with the low-lying $s$ and $d$ states adjusted to
reproduce the main features reported in the recent experimental work [24]. In
such a model, the wave function is mostly governed by the $d_{5/2}$ component
($\sim$ 80%), and the $p$-wave admixture is minimal ($\lesssim 2\%$) and comes
from the non-resonant continuum in 16B. The calculated matter radii for 14Be
and 17B are 3.0 and 2.8 fm, respectively, which compare well with the values
reported in Ref. [41] from interaction cross sections.

Using these structure inputs, the calculations capture the general trend of
the average correlation angle as a function of the intrinsic momentum, as
shown in Fig. 4. Figure 4a corresponds to 11Li and was already explored in
[28] with the same theoretical description. It should be remarked that for
missing momenta $k_{y}\gtrsim 0.5$ fm-1, the distribution is affected by the
core-proton interaction, so it is unreliable to extract nuclear structure
information from that region [28]. In the case of 14Be (Fig. 4b), the
calculated average correlation angle (blue solid line) follows the trend of
the experimental data but the results for intrinsic momenta smaller than 0.5
fm-1 are somewhat overestimated. Meanwhile, for 17B (Fig. 4c), the theoretical
model describes the maximum even with only a 2% $p$-wave admixture. This
remarkable sensitivity of the maximum to small opposite-parity components was
already noted in [28].

Only for 14Be there are significant differences between theoretical
calculation and experimental data. To understand these differences, we note
that in the analysis of the 13Be energy distribution in [33] the three-body
model used in this work was suggested to be missing some core-excited
components. Different components of the 12Be core can give opposing
contributions to the average correlation angle, as shown in Fig. 4b, where the
12Be$(0^{+}_{gs})$ component’s distribution (red dashed) goes over 90∘ at low
momenta, while the excited 12Be$(2^{+})$’s contribution (orange dashed) goes
under 90∘. Among the missing components in the used model, those where the
12Be core is in its first excited $0^{+}_{2}$ state are particularly
significant, since they are more likely to be populated, as its angular
momentum and parity are those of the 14Be ground state. To estimate their
effect, we note that the 12Be$(0^{+}_{2})$ state is usually described as an
orthogonal partner of the $0^{+}$ ground state [42, 43], with opposite
relative sign between its positive and negative-parity components when
compared to 12Be$(0^{+}_{gs})$. Therefore the components with
12Be($0^{+}_{2}$) should present a correlation angle smaller than 90∘
(opposite to 12Be$(0^{+}_{gs})$). Tentatively, for the correlation angle as a
function of missing momentum, we have assigned to the 12Be($0^{+}_{2}$)
components a distribution equal to that of 12Be$(0^{+}_{gs})$ but mirrored
around 90∘, and a weight of 16%, similar to the 20% obtained with the three-
body model for the similar-energy 12Be$(2^{+})$. This estimation produces the
magenta dot-dashed line, whose agreement with the data is much improved,
pointing to the excitation of the core having a significant effect in the
dineutron correlation, which was already indicated in [44]. Therefore, the
effect of the core may be responsible for 14Be and 17B showing a similar
correlation angle, despite their very distinct admixture of different-parity
components.

Figure 4:  Average correlation angle as a function of the intrinsic momentum
for (a) 11Li [28], (b) 14Be and (c) 17B, compared to the theoretical
calculations.

At this point, a natural question arises about how to compare the degree of
dineutron correlation among different nuclei. A possible criterium is based
solely on the experimental results, by comparing the maximum correlation
angle. The maximum correlation angle for 11Li, 14Be and 17B occurs at
$k_{y}=0.25$ fm-1 and corresponds to 100.0(2)${}^{+29}_{-29}$,
95.9(10)${}^{+29}_{-29}$ and 96.4(19)${}^{+29}_{-29}$ degrees, respectively.

A second criterion is to make use of the theoretical models employed. The
theoretical calculations in Fig. 4 give the average maximum values of 98.0
(11Li), 96.6 (14Be) and 95.4 (17B) degrees, which compare well to the
experimental results. It is worth noting that the corresponding three-body
models give rise to maximum of the two-neutrons wave function density around
the minimum of the average interneutron distance, as discussed in Ref. [8],
and this feature is directly linked to the present observations in momentum
space. The three-body model allows also to draw the ground-state probability
density as a function of the Jacobi-$T$ coordinates rnn and rc-nn, i.e. the
distance among the two neutrons and the two neutrons baricenter with respect
to the core. This is shown in Fig. 5 for the three cases considered and allows
to gain more insight on the configuration of the neutrons. In a purely non-
correlated scenario, the distributions would present equal weights at both
sides of the orange lines in the figure, which delimit two distinct regions
within the hyperspherical description of three-body nuclei [45]. Local maxima
above this line, i.e., for small $r_{nn}$, are usually associated to the
dineutron configuration, whereas the peaks below it correspond to the so-
called “cigar”-like structure. The dominance of one of these structures is
associated to correlations. We can see that a clear dineutron peak is obtained
for the three nuclei. For 11Li (Fig. 5a) the dineutron peak is clearly
dominant, with only a relatively small fraction of the probability exploring
larger neutron-neutron distances. For 14Be (Fig. 5b), the two configurations
are clearly separated, with the dineutron still being more pronounced. In the
case of 17B (Fig. 5c), three maxima appear (this is a consequence of the
dominant $d$-wave content of the ground state).

To quantify the degree of dineutron development for each nucleus, we may
define the quantity

$\chi=\frac{P_{d}-P_{c}}{P_{d}+P_{c}},$ (4)

where $P_{d}$ and $P_{c}$ are the integrated probabilities above and below the
symmetry lines in Fig. 5, i.e., $P_{d}$ is somehow a measure of the dineutron
component, while $P_{c}$ is related to the cigar component. Indeed, with this
definition $\chi=1$ ($-1$) would correspond to a “pure” dineutron (cigar). The
integration for 11Li, 14Be and 17B within the present calculations yields
$\chi=0.43$, $0.32$ and $0.19$, respectively. In this case, both criteria
agree and support the fact that the dineutron correlation is stronger for
11Li. The theoretical model also permits the extraction of the average opening
angle in configuration space, obtaining
$\left\langle\theta_{r}\right\rangle=66.9^{\circ}(^{11}\mathrm{Li})$,
$67.1^{\circ}(^{14}\mathrm{Be})$ and $77.4^{\circ}(^{17}\mathrm{B})$. The
results for ${}^{11}\mathrm{Li}$ and ${}^{14}\mathrm{Be}$ are consistent to
those presented in [15], while the angle for ${}^{17}\mathrm{B}$ is similar to
that presented for ${}^{6}\mathrm{He}$. Since both nuclei have very little
admixture of different-parity components, their opening angles should be
comparable.

Figure 5: Ground-state probability density, as a function of the distance
between the halo neutrons ($r_{nn}$) and that between the center of mass of
the $nn$ pair and the core ($r_{c\text{-}nn}$) for a) 11Li, b) 14Be and c)
17B. In a purely non-correlated case, the probability would be symmetric
around the orange lines.

From Fig. 6, one can extract the correlation between the average rnn and rc-
nn. The minimum of rnn corresponds to a dineutron configuration, and its
position signals the region of the nucleus where the calculation predicts the
dineutron correlation to be stronger. We can notice that this occurs for
r${}_{c-nn}=3\text{-}4$ fm, corresponding to the nuclear periphery, again
supporting the results in [14] and generalizing them to 14Be and 17B. As
discussed in [8], this behaviour can be interpreted as a transition from BCS-
like correlations in the interior to a BEC-like one, the spatially compact
dineutron, around the surface. The probability density for the two valence
nucleons of 11Li is also displayed with a grey area. One can notice that the
probability maximum and the inter-nucleon distance minimum are attained around
the same value of r${}_{c-nn}\sim 3$ fm, which makes the dineutron
configuration dominant for the two valence nucleons. Within the adopted
theoretical framework, the maximum in rc-nn corresponding to the nuclear
surface can be associated to the maximum at low intrinsic momentum ky,
validating the use of ky as a proxy for peripherality.

Figure 6:  Root mean square inter-neutron distance rnn as a function of the
distance between the core and the center of mass of the two neutrons rc-nn for
11Li, 14Be and 17B. The shaded area represents the probability distribution of
the two neutrons in 11Li, obtained by direct integration in Fig. 5.

##  4  Conclusions

We have presented a comparative study of dineutron correlation in three
Borromean systems, 11Li, 14Be and 17B, based on the average correlation angle
as a function of the intrinsic momentum of the removed neutron. This work
follows the seminal work of Kubota _et al._ [14] who first proposed to use
this observable to probe the location of dineutron correlation inside the
nucleus, and extends the study to 14Be and 17B. A dineutron correlation
appears in the periphery of 14Be and 17B as well, but is damped compared to
11Li.
This study provides the first experimental hint of the universality of
dineutron correlation in the low-density surface of Borromean nuclei. Even
while fast nucleon removal induced by high-energy quasi-free scattering is the
tool of choice to reduce the effect of final-state interactions, consistent
measurements using different probes may help to confirm the universal
character of our observation. The damping of dineutron correlation in 14Be is
interpreted as due to the presence of configurations with an excited core,
that can be predicted within the three-body model. Higher statistics data
incorporating gamma-ray coincidences, which enable core excitations to be
probed, could be used to investigate this explanation.

## Acknowledgements

This work has been supported by the European Research Council through the ERC
Starting Grant No. MINOS-258567. J.C., M.G.R. and A.M.M. acknowledge financial
support by MCIN/AEI/10.13039/501100011033 under I+D+i project No.
PID2020-114687GB-I00 and under grant IJC2020-043878-I (also funded by
“European Union NextGenerationEU/PRTR”), by the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 101023609, by the Consejería de Economía, Conocimiento, Empresas
y Universidad, Junta de Andalucía (Spain) and “ERDF-A Way of Making Europe”
under PAIDI 2020 project No. P20_01247, and by the European Social Fund and
Junta de Andalucía (PAIDI 2020) under grant number DOC-01006. J.G., F.M.M. and
N.A.O. acknowledge partial support from the Franco-Japanese LIA-International
Associated Laboratory for Nuclear Structure Problems as well as the French
ANR14-CE33-0022-02 EXPAND. Z.K. and L.S. acknowledge partial support by the
Institute for Basic Science (IBS-R031-D1). S.P. acknowledges the support of
the UK STFC under contract numbers ST/L005727/1 and ST/P003885/1 and the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID
279384907 - SFB 1245.

## References

  * [1] P. Hansen, B. Jonson, Europhys. News 4 (1987) 409.
  * [2] I. Tanihata, H. Savajols, R. Kanungo, Prog. Part. Nucl. Phys. 68 (2013) 215–313.  doi:10.1016/j.ppnp.2012.07.001.
  * [3] K. Hagino, H. Sagawa, Phys. Rev. C 72 (2005) 044321.  doi:10.1103/PhysRevC.72.044321.
  * [4] K. Hagino, H. Sagawa, Few-Body Systems 57 (2016) 185.  doi:10.1007/s00601-015-1027-3.
  * [5] F. Catara, A. Insolia, E. Maglione, A. Vitturi, Phys. Rev. C 29 (3) (1984) 1091–1094.  doi:10.1103/PhysRevC.29.1091.
  * [6] A. B. Migdal, Yadern. Fiz. 16 (2) 427, English translation Sov. J. Nucl. Phys., 16 (1973) 238.
  * [7] M. Matsuo, Phys. Rev. C 73 (2006) 044309.  doi:10.1103/PhysRevC.73.044309.
  * [8] K. Hagino, H. Sagawa, J. Carbonell, P. Schuck, Phys. Rev. Lett. 99 (2007) 022506\.  doi:10.1103/PhysRevLett.99.022506.
  * [9] Y. Oganessian, V. Zagrebaev, J. Vaagen, Phys. Rev. Lett. 82 (1999) 4996.  doi:10.1103/PhysRevLett.82.4996.
  * [10] Y. Sun, et al., Phys. Lett. B 814 (2021) 136072.  doi:https://doi.org/10.1016/j.physletb.2021.136072.
  * [11] T. Nakamura, et al., Phys. Rev. Lett. 96 (2006) 252502.  doi:10.1103/PhysRevLett.96.252502.
  * [12] H. Simon, et al., Phys. Rev. Lett. 83 (1999) 496.  doi:10.1103/PhysRevLett.83.496.
  * [13] H. Simon, et al., Nucl. Phys. A 791 (2007) 267.  doi:10.1016/j.nuclphysa.2007.04.021.
  * [14] Y. Kubota, et al., Phys. Rev. Lett. 125 (2020) 252501.  doi:10.1103/PhysRevLett.125.252501.
  * [15] C. A. Bertulani, M. S. Hussein, Phys. Rev. C 76 (2007) 051602.  doi:10.1103/PhysRevC.76.051602.
  * [16] K. Hagino, H. Sagawa, Phys. Rev. C 89 (2014) 014331.  doi:10.1103/PhysRevC.89.014331.
  * [17] K. Hagino, H. Sagawa, Phys. Rev. C 93 (2016) 034330.  doi:10.1103/PhysRevC.93.034330.
  * [18] H. Esbensen, G. Bertsch, Nucl. Phys. A 542 (2) (1992) 310–340.  doi:10.1016/0375-9474(92)90219-A.
  * [19] K. Hagino, H. Sagawa, Phys. Rev. C 76 (2007) 047302.  doi:10.1103/PhysRevC.76.047302.
  * [20] F. M. Marqués, et al., Phys. Lett. B 476 (3) (2000) 219–225.  doi:https://doi.org/10.1016/S0370-2693(00)00141-6.
  * [21] F. M. Marqués, et al., Phys. Rev. C 64 (2001) 061301.  doi:10.1103/PhysRevC.64.061301.
  * [22] I. Tanihata, et al., Phys. Rev. Lett. 100 (2008) 192502.  doi:10.1103/PhysRevLett.100.192502.
  * [23] L. V. Chulkov, et al., Phys. Rev. Lett. 79 (1997) 201–204.  doi:10.1103/PhysRevLett.79.201.
  * [24] Z. H. Yang, et al., Phys. Rev. Lett. 126 (2021) 082501.  doi:10.1103/PhysRevLett.126.082501.
  * [25] K. Hagino, H. Sagawa, Phys. Rev. C 93 (2016) 034330.  doi:10.1103/PhysRevC.93.034330.
  * [26] S. M. Wang, W. Nazarewicz, Phys. Rev. Lett. 126 (2021) 142501.  doi:10.1103/PhysRevLett.126.142501.
  * [27] A. Spyrou, et al., Phys. Rev. Lett. 108 (2012) 102501.  doi:10.1103/PhysRevLett.108.102501.
  * [28] J. Casal, M. Gómez-Ramos, Phys. Rev. C 104 (2021) 024618.  doi:10.1103/PhysRevC.104.024618.
  * [29] N. Fukuda, et al. 317 (2013) 323–332.  doi:https://doi.org/10.1016/j.nimb.2013.08.048.
  * [30] A. Obertelli, et al., Eur. Phys. J. A 50 (2014) 8.  doi:10.1140/epja/i2014-14008-y.
  * [31] T. Kobayashi, et al., Nucl. Instrum. Methods Phys. Res. B 317 (2013) 294.  doi:10.1016/j.nimb.2013.05.089.
  * [32] T. Nakamura, Y. Kondo, Nucl. Instrum. Methods Phys. Res. B 376 (2016) 156.  doi:10.1016/j.nimb.2016.01.003.
  * [33] A. Corsi, et al., Phys. Lett. B 797 (2019) 134843.  doi:10.1016/j.physletb.2019.134843.
  * [34] Y. Kikuchi, K. Ogata, Y. Kubota, M. Sasano, T. Uesaka, Progr. Theor. Exp. Phys. 2016 (2016) 103D03.  doi:10.1093/ptep/ptw148.
  * [35] Y. Kondo, et al., Phys. Lett. B 690 (2010) 245.  doi:10.1016/j.physletb.2010.05.031.
  * [36] J. Casal, M. Rodríguez-Gallardo, J. M. Arias, Phys. Rev. C 88 (2013) 014327\.  doi:10.1103/PhysRevC.88.014327.
  * [37] J. Singh, J. Casal, W. Horiuchi, L. Fortunato, A. Vitturi, Phys. Rev. C 101 (2020) 024310.  doi:10.1103/PhysRevC.101.024310.
  * [38] M. Gómez-Ramos, J. Casal, A. Moro, Phys. Lett. B 772 (2017) 115.  doi:10.1016/j.physletb.2017.06.023.
  * [39] Y. Aksyutina, et al., Phys. Lett. B 718 (2013) 1309.  doi:10.1016/j.physletb.2012.12.028.
  * [40] J. Casal, E. Garrido, Phys. Rev. C 102 (2020) 051304.  doi:10.1103/PhysRevC.102.051304.
  * [41] T. Suzuki, et al., Nuclear radii of 17,19b and 14be, Nuclear Physics A 658 (4) (1999) 313.  doi:https://doi.org/10.1016/S0375-9474(99)00376-0.
  * [42] A. O. Macchiavelli, et al., Phys. Rev. C 97 (2018) 011302.  doi:10.1103/PhysRevC.97.011302.
  * [43] H. T. Fortune, R. Sherr, Phys. Rev. C 85 (2012) 051303.  doi:10.1103/PhysRevC.85.051303.
  * [44] F. Kobayashi, Y. Kanada-En’yo, Phys. Rev. C 93 (2016) 024310.  doi:10.1103/PhysRevC.93.024310.
  * [45] M. Zhukov, B. Danilin, D. Fedorov, J. Bang, I. Thompson, J. Vaagen, Phys. Rep. 231 (1993) 151.  doi:10.1016/0370-1573(93)90141-Y.


††thanks: Present address: Department of Physics, Yale University, New Haven,
Connecticut 06520, USA

# Twisted superfluid and supersolid phases of triplons in bilayer honeycomb
magnets

Dhiman Bhowmick1 Abhinava Chatterjee1,2,∗ Prasanta K. Panigrahi2 Pinaki
Sengupta1 1School of Physical and Mathematical Sciences, Nanyang Technological
University, 21 Nanyang Link, Singapore 637371, Singapore
2Department of Physical Sciences, Indian Institute of Science Education and
Research Kolkata, Mohanpur - 741246, India

###### Abstract

We demonstrate that low-lying triplon excitations in a bilayer Heisenberg
antiferromagnet provide a promising avenue to realise magnetic analogs of
twisted superfluid and supersolid phases that were recently reported for two-
component ultracold atomic condensate in an optical lattice. Using a cluster
Gutzwiller mean field theory, we establish that Dzyaloshinskii-Moriya
interactions (DMI), that are common in many quantum magnets, stabilize these
phases in magnetic system, in contrast to pair hopping process that is
necessary for ultracold atoms. The critical value of DMI for transition to the
twisted superfluid and twisted supersolid phases depends on the strength of
the (frustrated) interlayer interactions that can be tuned by applying
external pressure on and / or shearing force between the layers. Furthermore,
we show that the strength of DMI can be controllably varied by coupling to
tailored circularly polarized light. Our results provide crucial guidance for
the experimental search of twisted superfluid and supersolid phases of
triplons in real quantum magnets.

###### pacs:

Valid PACS appear here

##  I Introduction

The observation of twisted, multi-orbital superfluid in binary mixtures of
ultracold 87Rb atoms in two different hyperfine states on a honeycomb optical
lattice has attracted heightened interest in this novel quantum state of
matter Soltan-Panahi  _et al._ (2011). The twisted superfluid (or twisted
supersolid) state is characterised by a complex order parameter – the phase of
the local superfluid order parameter at each site changes continuously forming
a“twisting pattern”, thus breaking time reversal symmetry spontaneously.
Interestingly, complex order parameters have experimentally been shown to be
associated with other novel strongly correlated phases such as the
superconducting states of Sr_2RuO_4 Luke  _et al._ (1998); Wysokinski (2019)
and UPt_3 Luke  _et al._ (1993) and the pseudo-gap state in the cuprate high-
Tc superconductor, B-2212 Kaur and Agterberg (2003); Kaminski  _et al._
(2002). A detailed understanding of the twisted superfluid state can help gain
insight into these states as well. Subsequent theoretical studies have shown
that the extended Bose Hubbard model with an additional pair hopping term can
stabilize a twisted superfluid (TSF) ground state over a finite range of
parameters Jürgensen  _et al._ (2015).

Quantum magnets have long served as a versatile platform for realizing
magnonic analogs of complex bosonic phases, often under less extreme
conditions. For example, temperature needed for Bose Einstein condensation
(BEC) of magnons varies from a few Kelvins in many quantum magnets Tanaka  _et
al._ (2001); Coldea  _et al._ (2003); Singh and Johnston (2007) to room
temperature in Yttrium Iron Garnet (YIG) thin films Demokritov  _et al._
(2006); Dzyapko  _et al._ (2007), in contrast to nano-Kelvin temperature
scales required for BEC in ultracold atoms Bradley  _et al._ (1995); Davis
_et al._ (1995). In this work, we show that twisted superfluid and twisted
supersolid phases of magnons are realized in a bilayer honeycomb Heisenberg
model. Interestingly, pair hopping process of magnons is not essential for
stabilizing these phases Jürgensen  _et al._ (2015), in contrast to ultracold
atomic systems. Instead, a next nearest neighbor (NNN) Dzyaloshinskii-Moriya
interaction (DMI) – which is present in many quantum magnets – is sufficient
to yield field induced twisted superfluid (TSF) and twisted supersolid (TSS)
phases over wide ranges of parameters.

Our paper is structured as follows. In Sec. II, we introduce the microscopic
spin Hamiltonian and describe the physics in the non-interacting limit by
deriving the tight-binding triplon Hamiltonian and triplon band structure. The
cluster Gutzwiller mean field theory (CGMFT) is introduced in Sec. III which
is used to calculate order-parameters in interacting limit. In Sec. IV we
present the results of our study in the form of the order parameters and phase
diagram as a function of magnetic field, DMI and Heisenberg interactions. In
Sec.V, we propose the possible materials in search of TSS and TSF phases.
Finally, in Sec. V, we show that circularly polarized light can be used to
induce DMI greater than the cutoff DMI required to realize TSS (TSF), whereas
presence of frustration among interlayer and intralayer interactions is
helpful to lower the cutoff value of DMI (see also Appendix. C). The principal
findings are summarized in section Sec. VI.

Figure 1:  (Color online) (a) The spins on each lattice site interact via a
strong interlayer anti-ferromagnetic coupling, resulting in a ground state of
a honeycomb lattice of dimers on each interlayer nearest-neighbor bonds (shown
by dashed-black line). (b) The ferromagnetic honeycomb lattice. (c) Triplon
band structure at $D=0.1J$, (d) Triplon band structure at $D=0.8J$. The other
parameters for the band structure are $J_{\perp}=10J$, $B_{z}=0.0$,
$J_{z}=0.0$.

##  II The bilayer honeycomb magnet and effective triplon model

We start with a $S=1/2$ Heisenberg antiferromagnet on a bilayer honeycomb
lattice with out of plane exchange anisotropy and Dzyaloshinkii-Moriya
interaction (DMI), schematically shown in Fig. 1(a) and described by the
Hamiltonian,

$\displaystyle\pazocal{H}=$ $\displaystyle
J_{\perp}\sum_{\begin{subarray}{c}i,\,m\in A\\\ n\in
B\end{subarray}}\mathbb{S}_{i,m}\cdot\mathbb{S}_{i,n}-B_{z}\sum_{i,m}S^{z}_{i,m}$
$\displaystyle+\sum_{\left\langle
i,j\right\rangle,m}\left[J(S^{x}_{i,m}S^{x}_{j,m}+S^{y}_{i,m}S^{y}_{j,m})+J_{z}S^{z}_{i,m}S^{z}_{j,m}\right]$
$\displaystyle+D\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle,m}\nu_{ij}\hat{z}\cdot(\mathbb{S}_{i,m}\times\mathbb{S}_{j,m}).$
(1)

$\mathbb{S}_{i,m}$ denotes the spin operator at $i$-th interlayer bond at
layer $m$ ($m\in\\{A,B\\}$). $\langle\ldots\rangle$ and
$\langle\langle\ldots\rangle\rangle$ denote the nearest-neighbor (NN) and
next-nearest-neighbor (NNN) inter-dimer bonds respectively. $J_{\perp}(>0)$ is
the strength of (isotropic) interlayer Heisenberg interaction, while $J_{z}$
($>0$) and $J$ ($>0$) denote the Ising and XX type intralayer NN Heisenberg
exchange interactions respectively. $D$ is the Dzyaloshinskii–Moriya
interaction (DMI) which is constrained by the symmetry of the lattice to
intralayer NNN bonds; $\nu_{ij}=+1$, if $i\to j$ forms part of a
counterclockwise closed loop connecting the NNN sites in a hexagonal plaquette
in each layer (see counterclockwise circular arrows in Fig. 1(a)-(b)) and
$\nu_{ij}=-1$ otherwise. Finally $B_{z}S_{i,m}^{z}$ describes a Zeeman
coupling between a local spin-moment and an external longitudinal magnetic
field. It is noticeable that the change in sign of $J$ and $D$ do not alter
the magnetic ground state as well as excitations above the ground state,
whereas $J_{z}$ and $J_{\perp}$ are strictly positive in this study.

For $J_{\perp}\gg|J|$, the ground state of the system is a product of the
singlet dimers on each interlayer NN bond. In this limit, the lowest
excitations of the system are triplons, which are localized $S=1$
quasipartices. An out-of-plane magnetic field lowers the energy of the
$S_{z}=+1$ triplons and at a critical magnetic field, it crosses the energy of
the singlet state, populating the ground state with a finite density of
triplons. The other triplon branches ($S^{z}=0$ and $S^{z}=-1$) are separated
by a large energy gap. At low temperatures, one can restrict the local Hilbert
space of the dimers to the singlet and $S^{z}=+1$ triplon. By treating the
triplons as bosonic quasiparticles, one can formulate a description of the low
energy physics of the system in terms of hard core bosons which is known as
bond-operator formalism (see also Appendix.A). In this formalism, the zero
field ground state made of singlets on each interlayer bonds is considered an
empty lattice with number of triplons $n_{i}=0\;\forall i$. At the critical
field, a finite density of triplons is generated which increases with
increasing field. The inter-dimer exchange interactions induce an effective
hopping of the triplons. This delocalization induces a BEC of $S_{z}=+1$
triplons in the ground state. In the spin language, this corresponds to an
canted antiferromagnetic order with a spontaneously broken U(1) symmetry.
Considering singlets as a vacuum state in the system and triplon ($S_{z}=+1$)
as a hard-core bosonic quasi-particle excitations in vacuum of singlets, we
can use the bond operator formalism to express the bi-layer spin Hamiltonian
Eq. 1 as an effective triplon Hamiltonian on a single-layer honeycomb lattice
Zapf  _et al._ (2014); Sachdev and Bhatt (1990); Romhányi  _et al._ (2011)
(see Appendix. A),

$\displaystyle\pazocal{H}=$ $\displaystyle\frac{J}{2}\sum_{\left\langle
ij\right\rangle}\left[\hat{t}^{\dagger}_{i}\hat{t}_{j}+\text{H.c.}\right]+\frac{\mathfrak{i}D}{2}\sum_{\left\langle\left\langle
ij\right\rangle\right\rangle}\nu_{ij}\left[\hat{t}^{\dagger}_{i}\hat{t}_{j}-\text{H.c.}\right]$
$\displaystyle+\left(\frac{J_{\perp}}{4}-B_{z}\right)\sum_{i}\hat{t}^{\dagger}_{i}\hat{t}_{i}+\frac{J_{z}}{2}\sum_{\left\langle
ij\right\rangle}\hat{n}_{i}\hat{n}_{j}$ (2)

where, $\hat{t}_{i}$ is the triplon annihilation operator, $\hat{n}_{i}$ is
the triplon number operator and index-$i$ represents site index of effective
single-layer honeycomb lattice (which is equivalent to interlayer bond-index
of bi-layer honeycomb lattice). The first two terms represent hopping of
triplons between NN and NNN neighbor dimers respectively, the third term is an
on-site potential (effectively a chemical potential) and the last term
describes the effects of NN interaction between triplons. It is noticeable
that the NNN hopping has a complex weight which renders the Hamiltonian
unsuitable for quantum Monte Carlo simulations. When $J_{z}\approx 0$, Eq.2
reduces to a tight binding model of non-interacting triplons. The Bloch
Hamiltonian in the momentum basis, in terms of the momentum space triplon
operators, is determined via Fourier transformation as,

$\pazocal{H}=\sum_{\mathbb{k}}\Psi_{\mathbb{k}}^{\dagger}\left[g(\mathbb{k})\sigma_{0}+\mathbb{h}\cdot\bm{\sigma}\right]\Psi_{\mathbb{k}},$
(3)

where, $\Psi_{\mathbb{k}}=(\hat{a}_{\mathbb{k}},\hat{b}_{\mathbb{k}})^{T}$.
$\hat{a}_{\mathbb{k}}$ ($\hat{b}_{\mathbb{k}}$) denotes the ${\bf k}$-space
triplon annihilation operator on sublattice-a (b) as shown in Fig. 1(a)-(b)
and $\bm{\sigma}$ is the pseudo-vector of Pauli matrices and $\sigma_{0}$ is
the two-dimensional identity matrix. The coefficients of the $\sigma$-matrices
in the Bloch Hamiltonian are $g(\mathbb{k})=(J_{\perp}/4)-B_{z}$,
$h_{x}(\mathbb{k})=(J/2)\sum_{i}\cos(\mathbb{k}\cdot\bm{\alpha}_{i})$,
$h_{y}(\mathbb{k})=(J/2)\sum_{i}\sin(\mathbb{k}\cdot\bm{\alpha}_{i})$,
$h_{z}(\mathbb{k})=D\sum_{i}\sin(\mathbb{k}\cdot\bm{\beta}_{i})$, where
$\bm{\alpha}_{i}$ and $\bm{\beta}_{i}$ are the NN and NNN vectors for each
layer respectively. $\bm{\beta}_{i}$’s are chosen such that they form a
counter-clockwise triangular loop for sites in sublattice-a in a hexagonal
plaquette and clockwise triangular loop for sites in sublattice-b. The energy
eigenvalues are given by,

$E^{\pm}(\mathbb{k})=g(\mathbb{k})\pm|h(\mathbb{k})|.$ (4)

The band dispersion is shown in Fig. 1(c) and Fig. 1(d) for two values of DMI.
For ${\bf D}=0$, the energy spectrum is identical to that of graphene, with a
linear band crossing of the upper and lower bands at the Dirac points $K$ and
$K^{\prime}$. A finite DMI breaks time reversal symmetry and opens a band-gap
$6\sqrt{3}D$ at these points. The energy of the lower band at $\Gamma$ and $K$
($K^{\prime}$)-points are respectively given by (at $B_{z}=0$),

$E_{\Gamma}=\frac{J_{\perp}}{4}-\frac{3|J|}{2},\,\,E_{K}=\frac{J_{\perp}}{4}-\frac{3\sqrt{3}}{2}|D|.$
(5)

In the absence of DMI, the energy minimum is located at the center of the
Brillouin zone, the $\Gamma$ point. For a finite but small DMI, the band
minimum remains at $\Gamma$-point (Fig. 1(c)). Increasing DMI to
$|D|>|J|/\sqrt{3}$ shifts the band minimum from $\Gamma$-point to two
degenerate minima at the $K$ and $K^{\prime}$ (Fig. 1(d)). Thus with changing
DMI, the ground state changes from a one component BEC (condensation momentum
${\bf k}=0$) to a two component BEC (condensation momenta at k = $K$ and
$K^{\prime}$). The transition happens at $|D|=|J|/\sqrt{3}$ independent of
$J_{\perp}$ and $B_{z}$ for small $B_{z}$.

In presence of repulsive interaction between triplons ( the last term in Eq. 2
) co-existence of triplons at $K$ and $K^{\prime}$ points costs no additional
energy classically ( see Appendix. B ) , but the quantum-fluctuations around
$\mathbb{K}$ and $\mathbb{K}^{\prime}$ points introduces an energy cost Li
_et al._ (2014). Thus spontaneous breaking of valley-symmetry is energetically
favoured and the quantum fluctuation will lead to a ground state with only
single valley condensation which is known as “quantum order by disorder”
effect Li  _et al._ (2014); Adhikari  _et al._ (2010); Henley (1989). Thus, in
the limit of weak interaction, the super-fluid order parameter at lattice site
$\mathbb{r}$ is either $|b|e^{i\mathbb{K}\cdot\mathbb{r}+\phi}$ or
$|b|e^{i\mathbb{K}^{\prime}\cdot\mathbb{r}+\phi}$ depending on the valley for
the Bose-Einstein condensate, where $\phi$ is a global phase independent of
position of lattice site. This spontaneous breaking of valley symmetry
transforms the superfluid order-parameter from a real to a complex value and
the resultant BEC is known as twisted superfluid Jürgensen  _et al._ (2015).
Based on this, the ground state on the honeycomb lattice is no longer a
bipartite lattice in twisted superfluid phase, but become a lattice with six-
sublattices as shown in the figure Fig. 2. Results based on CGMFT also
supports this scenario as shown in section Sec. III.

In absence of interaction $J_{z}$, Eq. 2 transforms into well-known bosonic
topological Haldane model. At finite temperatures due to non-zero magnetic
excitations, this model is known to exhibit a finite thermal Hall effect due
to presence of non-zero Berry-curvature of the bands. On this basis it is
expected that Bose-Einstein condensate at $\Gamma$-point or $K$
($K^{\prime}$)-point would provide a non-zero Hall conductance even at zero
temperature Li  _et al._ (2014). However, we find that the Berry-curvature at
the condensation momenta – $\Gamma$-point when $|D|<|J|/\sqrt{3}$ and $K$
($K^{\prime}$)-point when $|D|>|J|/\sqrt{3}$ – is zero. Since the density of
magnons is concentrated around the condensation momentum, and is vanishingly
small away from it, the thermal Hall response of the Bose-Einstein condensate
is vanishingly small in our model.

An analogous bosonic model can also be observed in terms of Matsubara-Matsuda
bosons in magnetically ordered honeycomb ferromagnets (see Fig. 1(b) ) or
antiferromagnets. Thus the non-trivial phases like twisted-superfluid or
twisted-supersolid are also expected to emerge in magnetically ordered
systems. However, in practice, the presence of anisotropies in magnetically
ordered systems break the $U(1)$-symmetries destroying conservation of numbers
of particles which in-turn preclude a long time superfluidity response in this
systems Zapf  _et al._ (2014); Giamarchi  _et al._ (2008). Moreover due to
presence of strong exchange interactions ($J,\,J_{z}$ in Eq. 2) in general,
the phase transitions are also difficult to study varying the external
parameters in these systems Giamarchi  _et al._ (2008). That is why the
dimerized paramagnets are well suited to study superfluid phases and
superfluid-Mott transitions. Hence we focus on bilayer dimerized honeycomb
paramagnetic system to study the twisted superfluid and twisted supersolid
phases.

##  III  Cluster-Gutzwiller Mean Field Theory (CGMFT) and Observables

Figure 2:  (Color online) The cluster construction for CGMFT. There are 18
sites in the cluster which is located within the dashed black box. A periodic
boundary condition is applied along the horizontal direction and the mean-
field boundary condition is applied along the vertical direction. The
background of the cluster sites is denoted by a pink shade and the background
of the mean-field sites are denoted by blue shade.

The cluster Gutzwiller mean field theory or CGMFT Lühmann (2013); Jürgensen
_et al._ (2015) – equivalently, cluster mean field theory Chen and Yang
(2017); Gelfand  _et al._ (1989); Yamamoto  _et al._ (2014); Moreno-Cardoner
_et al._ (2014); Yamamoto  _et al._ (2012a, b); Singh  _et al._ (2014), self-
consistent cluster mean field theory Hassan  _et al._ (2007), multi-site mean
field theory McIntosh  _et al._ (2012), hierarchical mean field approach Isaev
_et al._ (2009a, b), composite boson mean field theory Huerga  _et al._ (2013)
– is a powerful technique to study superfluid phases in bosonic many body
systems with complex hopping terms. CGMFT improves over the conventional
single-site mean field approach by taking into account the short range
correlations present within a small lattice-cluster using exact-
diagonalization. Furthermore it is an alternative numerical method to study
the quantum systems like we described in section Sec. II, where sign problems
arises in quantum Monte-Carlo methods due to complex hopping terms or
geometric frustration Yamamoto  _et al._ (2017); Huerga  _et al._ (2014).
Whereas conventional mean-field theories Choudhury and Mueller (2013) and
exact-diagonalization of small systems Cao  _et al._ (2015) fail to show
existence of TSF and TSS phases, Density Matrix Renormalization Group for one
dimensional systems Lühmann (2016) and CGMFT for higher dimensional systems
Jürgensen  _et al._ (2015) are better alternatives for search of these non-
trivial phases.

We explore the ground state phases of the effective triplon Hamiltonian with
CGMFT by decomposing the system into clusters (pink shaded region) and mean-
field region (blue shaded region) as shown in figure Fig. 2. The effective
mean-field Hamiltonian of the cluster is given as,

$\pazocal{H}_{C}^{\text{eff}}=\pazocal{H}_{C}+\pazocal{H}_{\delta C},$ (6)

where, $\pazocal{H}_{C}$ is the Hamiltonian as in equation Eq.( 2) within the
cluster and $\pazocal{H}_{\delta C}$ is the Hamiltonian which takes into
account the interactions among the boundary sites of the cluster and the mean-
field region. The form of the boundary Hamiltonian is given by,

$\displaystyle\pazocal{H}_{\delta C}=\frac{J}{2}$
$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{\left\langle
i,j\right\rangle}\left[\hat{t}_{i}^{\dagger}\left\langle\hat{t}_{j}\right\rangle+\text{H.c.}\right]+\frac{iD}{2}\sideset{}{{}^{\prime}}{\sum}_{\left\langle\left\langle
i,j\right\rangle\right\rangle}\left[\hat{t}_{i}^{\dagger}\left\langle\hat{t}_{j}\right\rangle-\text{H.c.}\right]$
$\displaystyle+\frac{J_{z}}{2}\sideset{}{{}^{\prime}}{\sum}_{\left\langle
i,j\right\rangle}\hat{n}_{i}\left\langle\hat{n}_{j}\right\rangle,$ (7)

where the primed summations are over the boundary site-$i$ connected to the
mean-field site-$j$. $\left\langle\hat{t}_{j}\right\rangle$ and
$\left\langle\hat{n}_{j}\right\rangle$ are two mean-field parameters denoting
the superfluid order parameter and occupation number of triplons at site-$j$
respectively. We choose six inequivalent sites in each cluster (denoted by
different patterns in figure Fig. 2) to give a total of 12 mean-field
parameters. The ground state in the different parameter regimes are obtained
by evaluating these mean field parameters self-consistently in the following
manner,

  * (i)

Choose an initial set of mean field parameters
$\left\\{\left\langle\hat{t}_{j}\right\rangle,\,\left\langle\hat{n}_{j}\right\rangle\right\\}$,
$j=1,\dots,6$ and then exactly diagonalize the effective Hamiltonian of the
cluster $\pazocal{H}^{\text{eff}}_{C}$.

  * (ii)

Calculate new mean field parameters
$\left\langle\hat{t}^{\prime}_{j}\right\rangle$ and
$\left\langle\hat{n}^{\prime}_{j}\right\rangle$ from the sites within blue-
dashed rectangle in the figure Fig. 2 which reside within the cluster.
Periodic boundary condition is chosen along horizontal direction to eliminate
any boundary effect on the sites within the blue-dashed rectangle, so that the
mean field parameters obtained from those sites are free from boundary
effects.

  * (iii)

The initial and final set of mean field parameters are compared using the
tolerance

$\epsilon=\sum_{j}\left|\left\langle\hat{t}^{\prime}_{j}\right\rangle-\left\langle\hat{t}_{j}\right\rangle\right|+\sum_{j}\left|\left\langle\hat{n}^{\prime}_{j}\right\rangle-\left\langle\hat{n}_{j}\right\rangle\right|.$
(8)

If the tolerance $\epsilon$ is less than a certain cutoff then the obtained
mean-field parameters correspond to the ground state of the system. Otherwise
the step-(i) is repeated with new values of mean-field parameters
$\left\langle\hat{t}_{j}\right\rangle=\left\langle\hat{t}^{\prime}_{j}\right\rangle$
and
$\left\langle\hat{n}_{j}\right\rangle=\left\langle\hat{n}^{\prime}_{j}\right\rangle$.

We set the cutoff as $10^{-10}$ and start the simulation with different
initial mean-field parameter sets for a fixed set of parameters $J$, $D$,
$J_{\perp}$, $B_{z}$, and $J_{z}$. In general, the simulations with different
initial mean field parameter sets give different ground-states at the boundary
of two phases and we selected the phase with minimum energy as the ground
state.

After obtaining the ground state via self-consistent determination of the mean
field parameters, four order-parameters are calculated to identify the nature
of the ground state phase of the system. The magnitude of superfluid order
parameter is given by,

$|b|=\text{max}\left[|\left\langle\hat{t}_{j}\right\rangle|\right],$ (9)

where max denotes maximum value of the parameter obtained out of six-sites
within the dashed blue-border in Fig. 2. Additioally, the average number of
particles per site $\left\langle n\right\rangle_{av}$ and difference in number
of particles between NN sites $\Delta n$ are also enumerated. The superfluid
order-parameter is a complex quantity and for twisted super-fluid phases in
our study the phase difference of superfluid order parameter $b$ among NNN
sites is obtained to be $\theta=120^{o}$ and otherwise $\theta=0^{o}$.

##  IV  Numerical results

Figure 3:  (Color online) The order parameters are plotted for parameter
values (a) $D=0.2J,\,J_{z}=J$, (b) $D=0.9J,\,J_{z}=J$, (c)
$D=0.2J,\,J_{z}=4J$, (d) $D=0.9J,\,J_{z}=4.5J$. $J_{\perp}$ is fixed at value
$10J$. Order-parameters $|b|$, $\left\langle n\right\rangle_{av}$ and $\Delta
n$ are plotted as function of magnetic field $B_{z}$ and denoted by red, black
and blue dotted lines respectively. The dots on the lines denote the points in
parameter space where the CGMFT is performed and the lines just connect the
points. Moreover the phase difference of super-fluid order parameter $\theta$
is shown in the right-side vertical-axis and denoted in green colour.
Different phases are indicated by different coloured shades. “E” and “O”
denote empty and fully-occupied phase of the system respectively. All other
phases are described in the main text.

Using CGMFT, we determine the order parameters $|b|$, $\left\langle
b\right\rangle_{av}$, $\Delta n$ and $\theta$ as a function of magnetic field
$B_{z}$ for different sets of the parameters $(D,J_{z})$. The evolution of the
order parameters and the resulting field driven phases are shown in Figure
fig. 3 for four illustrative points of the $(D,J_{z})$ parameter space. In
fig. 3, the DMI increases from the left-column of figures to the right-column
of figures, whereas the interaction $J_{z}$ increases from upper-row of
figures towards the lower-row of figures.

For weak DMI ($D=0.2J$), the field driven phase diagram resembles that of the
canonical extended Bose Hubbard model for hard core bosons Gan  _et al._
(2007) (see Fig. 3(a), (c)). The zero field ($B_{z}=0$) ground state
corresponds to a singlet phase, or equivalently an empty lattice in the
bosonic language. All the order parameters vanish in this limit. This remains
true at small values of the applied field reflecting a finite gap to lowest
excitations due to the singlet-triplet gap of the local dimers. When the
applied field exceeds a critical value, the gap is closed and the ground state
acquires a finite density of triplons. These field induced triplons form a
superfluid (SF) driven by the NN triplon hopping and is characterised by a
finite SF order parameter, $|b|$. The mismatch between the occupancy of the
two sublattices ($\Delta n$) remains zero, reflecting the uniform nature of
the SF phase. A vanishing twist angle ($\theta=0$) completes the
characterization of the phase as a normal superfluid. In the weak interaction
limit ($J_{z}<2J$) with increasing magnetic field, the density of triplons
increases monotonically till full saturation is reached at an upper critical
field when each dimer is occupied by a triplon. At saturation, all the order
parameters (except average density, $\langle n\rangle_{av}$, of triplons)
vanish denoting a lattice fully occupied by triplons. In this weak interaction
limit, the physics is similar to the non-interacting limit as described in
section II and so the qualitative feature can be well described using band
structure as in figure Fig. 1(c).

For strong interactions ($J_{z}>2J$), an intervening charge density wave (CDW)
phase, driven by the strong NN-interaction between triplons, appears in
addition to the phases discussed above (see Fig. 3(c)). With increasing
magnetic field, when the density of triplons reaches $\langle
n\rangle_{av}=1/2$, the triplons form a staggered CDW pattern where one of the
sublattices is fully occupied, while the other remains empty. The potential
energy-cost due to nearest neighbor interaction is minimized as there are no
nearest neighbor pairs. This is accompanied by a complete quenching of
superfluidity, since any hopping of triplons will necessarily involve
configurations with energetically costly multiple nearest neighbor pairs. The
CDW phase has a finite gap to the addition of any more triplons and the
density remains constant at $\langle n\rangle_{av}=1/2$ over a finite rage of
applied field. This phase is characterized by a vanishing superfluid order,
and a non-zero density mismatch between the two sublattices ($\Delta n$),
reflecting the staggered order. When the increasing field strength reaches a
critical value where the Zeeman energy gain due to increasing magnetization
(equivalently, adding more triplons) exceeds the potential energy cost of
nearest neighbor repulsion, the density of triplons starts to increase again,
resulting in another normal SF phase. Finally, as the field is increased above
a saturation value, $B_{sat}$, the ground state enters the fully polarized
phase.

The above argument for the appearance of interaction driven CDW phase at half-
filling does not apply for weak to moderate interaction strengths
($J_{z}<2J$), as the kinetic energy gain due to the delocalization of triplons
exceeds the potential energy cost of NN-interactions.

The sequence of field-driven phase changes markedly for strong DMI. As shown
earlier in Sec. II, in the non-interacting limit, the triplon band minimum
shifts from the $\Gamma$ point to the $K$ and $K^{\prime}$ (Fig. 1(d)) and the
BEC of triplons occur at finite momentum. A local minimum persists at the
center of the Brillouin zone, and the energy gap between the triplon-sector
and the singlet dimer sector, $E_{K}$ in Eq. 5 decreases with increasing DMI.
This behavior persists in the presence of weak to moderate interaction
($J_{z}<2J$) and is reflected in Fig. 3(b). For the present choice of
parameters, the energy of the lowest triplon excitation is vanishingly small.
The triplon density acquires a finite value for an infinitesimally small
$B_{z}$, and increases monotonically with the strength of the applied field.
In this regime, the triplons form a superfluid ($|b|>0$). More interestingly,
the complex NNN hopping process imparts a complex phase to the superfluid
order parameter, as seen by a finite expectation value of the twist angle
($\theta\neq 0$). In other words, the ground state in this parameter range is
a twisted superfluid (TSF). The triplon density increases monotonically, with
the ground state remaining a TSF, till the fully polarized phase is reached at
a saturation field, $B_{sat}$.

Finally in the strong DMI and strong interaction limit (see Fig. 3(d)),the
twisted superfluid is replaced by a twisted supersolid phase, in addition to
the appearance of an interaction-driven CDW phase at $\langle
n\rangle_{av}=1/2$ over a finite range of applied field. In the twisted
supersolid phase (TSS), the ground state is characterized by a finite $\Delta
n$ (density mismatch between the two sublattices), in addition to a complex
superfluid order parameter ($|b|\neq 0,\,\theta\neq 0,$). The finite density
difference between two sublattices provide the diagonal order concurrently
with the finite (twisted) superfluid ordering. It is surprising that the
ground state exhibits TSS order even at low triplon densities. This is
understood by recalling that the primary delocalization process in this
parameter regime involves the DMI-induced intra-sublattice complex next-
nearest-neighbor hopping. The strong NN-repulsion between the triplons further
suppresses inter-sublattice hopping processes, resulting in a preferential
occupation of one of the two sublattices at small densities.

Figure 4:  (Color online) Phase diagram at two different interaction values
(a) $J_{z}=J$ and (b) $J_{z}=4.5J$. Each dot or circle denotes the parameter
point where the CGMFT is performed. Different color denotes different phases
in parameter space. The empty and dotted black circles represent the empty and
fully occupied lattice phases, respectively.

The phase diagram in $B_{z}-D$ parameter space is shown in the figure Fig. 4
for two different values of interaction $J_{z}=J$ and $J_{z}=4.5J$. A
comparison of the phase diagrams at moderate ($J_{z}<2J$) and strong
($J_{z}>2J$) interactions reveal,

  * (i)

Appearance of CDW phase at half-filling in the strong interaction limit. The
CDW phase appears at half filling dividing the SF-region which appear at weak
to moderate $J_{z}$ into two SF-regions.

  * (ii)

In the strong interaction limit the TSF phase is replaced by a TSS phase.

##  V Material Realization

There are two main ingredients to realize TSF or TSS phases; firstly, a
bilayer honeycomb valence bond state is required; secondly, a DMI greater than
a critical value is required. Although, to the best of our knowledge, there
are no materials described in the literature which satisfies our model, a
rough idea is sketched here to obtain TSS or TSF phases based on the real
materials Bi3Mn4O12(NO3) and CrI3.

The material Bi3Mn4O12(NO3) consists of Mn4+ ions carrying spin
$S=\frac{3}{2}$ arranged in a bilayer honeycomb lattice with A-A stacking. No
magnetic ordering is observed down to the lowest temperatures. Theoretical
studies indicate that the disordered magnetic state may be an interlayer dimer
phase which is adiabatically connected to the direct products of singlets
Zhang  _et al._ (2016); Kandpal and van den Brink (2011); Zhang  _et al._
(2014); Oitmaa and Singh (2012); Ganesh  _et al._ (2011a), although a spin
liquid phase cannot be completely ruled out Mulder  _et al._ (2010); Wang
(2010); Okumura  _et al._ (2010); Clark  _et al._ (2011); Cabra  _et al._
(2011); Ganesh  _et al._ (2011b); Albuquerque  _et al._ (2011); CABRA  _et
al._ (2011); Mezzacapo and Boninsegni (2012); Li  _et al._ (2012); Ganesh  _et
al._ (2013); Gong  _et al._ (2013). If the ground state phase is a spin
liquid, an interlayer dimer phase can be induced by applying a pressure along
the perpendicular direction Ganesh  _et al._ (2011b); Oitmaa and Singh (2012).

On the other hand, experimentally well studied van der Waals material CrI3 is
a honeycomb ferromagnet where Cr2+-ions carry spin-$\frac{1}{2}$ momentum
forming a honeycomb lattice. This material has been well studied by tuning the
number of layers Huang  _et al._ (2017) as well as varying the pressure Li
_et al._ (2019). It is shown that the interlayer antiferromagnetic Heisenberg
exchange interaction changes linearly with the distance between two layers by
application of pressure on the material and can be achieved a interlayer
coupling twice as compared with the initial interlayer coupling Li  _et al._
(2019). Although a valence bond state is not detected in CrI3 (and there are
indications of structural phase transition under pressure) Li  _et al._
(2019), the pressure induced high interlayer coupling in CrI3 motivates for
searching of pressure induced valance bond states in a family of Van der Walls
honeycomb magnets CrBr3 Samuelsen  _et al._ (1971), CrGeTe3 Gong  _et al._
(2017), CrSiTe3 Williams  _et al._ (2015), FePS3 Lee  _et al._ (2016), NiPS3
Kuo  _et al._ (2016).

Another challenge for achieving the TSF and TSS phases is to obtain a high DMI
($|D|>J/\sqrt{3}$). The materials Bi3Mn4O12(NO3) and CrI3 are experimentally
predicted to possess DMI Okubo  _et al._ (2012); Chen  _et al._ (2018).
However the magnitude of DMI is small for realizations TSS and TSF states.
This problem can be overcome by application of circularly polarized light. A
circularly polarized light couples either to the charge Vinas Boström  _et
al._ (2020) or to the magnetic degrees of freedom Owerre (2017) and it has
been shown theoretically that both kinds of mechanisms give rise to scalar
spin chiral interaction in a honeycomb magnetic insulator. The synthetic
scalar spin chiral interaction as in reference Ref. Vinas Boström  _et al._
(2020) given by,

$\displaystyle\pazocal{H}_{\chi}$
$\displaystyle=\chi\sum_{\left\langle\left\langle
j,k\right\rangle\right\rangle,m}\nu_{jk}\hat{S}_{i,m}\cdot\left(\hat{S}_{j,m}\times\hat{S}_{k,m}\right)$
$\displaystyle\approx\frac{3\mathfrak{i}\chi}{4}\sum_{\left\langle\left\langle
j,k\right\rangle\right\rangle}\nu_{jk}\left\langle\hat{n}_{i}\right\rangle\left[\hat{t}_{j}^{\dagger}\hat{t}_{k}-\text{H.c.}\right],$
(10)

where the $i$-th site is the NNN neighbor of both $j$-th and $k$-th sites,
$\chi$ is synthetic scalar spin chirality. Comparing the equations Eq. 2 and
Eq. 10, the effective DMI for bond $\left\langle\left\langle
jk\right\rangle\right\rangle$ due to circularly polarized light is
$D=\frac{3\chi\left\langle n_{i}\right\rangle}{2}$. The scalar spin chirality
$\chi$ is shown to have a resonance for the frequency of light near
$\omega_{n}=\frac{1}{\hbar}\frac{U}{n}$, where $U$ is onsite interaction
representing electron-electron repulsion and $n$ is a positive integer. Thus
the effective DMI can be tuned as high as possible tuning the frequency of the
light nearby $\omega_{n}$. However the results in Ref. Vinas Boström  _et al._
(2020) is based on single-band extended Hubbard model; for a magnetic material
which cannot be described by a one band Hubbard model, one needs a more
careful theoretical treatment. On the other hand, application of light may not
be required for a real material if the cutoff of DMI to realize TSF or TSS
states is reduced due to presence of frustration among interlayer and
intralayer interactions (see Appendix. C).

##  VI Conclusion

To summarize, we have shown that a magnetic analog of the novel twisted
superfluid (TSF) state reported in recent experiments with ultracold atoms in
an optical lattice Soltan-Panahi  _et al._ (2011) can be realized in a bilayer
quantum antiferromagnet with realistic interactions. We show that TSF phase is
induced by DMI greater than a critical value. For Ising-like anisotropy of the
intra-plane Heisenberg interactions, the TSF phase is replaced by a twisted
supersolid (TSS) phase. While the strength of DMI required for the
stabilization of TSF and TSS phases ($D/J\gtrsim 0.5$) is not observed
natively in most quantum magnets, recent experiments have shown that a strong
DMI can be induced in thin films of insulating magnets, by forming
heterostructures with heavy metals (with strong spin-orbit coupling) Tacchi
_et al._ (2017); Heide  _et al._ (2008); Di  _et al._ (2015); Gastaldo  _et
al._ (2019). Our results show that circularly polarized light can also induce
large DMI by varying the frequency of light Vinas Boström  _et al._ (2020).
Finally, the presence of frustration among interlayer and intralayer
Heisenberg exchange interactions in a material can lower the value of the
critical DMI required to realize TSF and TSS phases (Appendix. C), thus
facilitating their experimental observation. We propose that the material
Bi3Mn4O12(NO3) Zhang  _et al._ (2016); Kandpal and van den Brink (2011); Zhang
_et al._ (2014); Oitmaa and Singh (2012); Ganesh  _et al._ (2011a) and family
of honeycomb magnets CrBr3 Samuelsen  _et al._ (1971), CrGeTe3 Gong  _et al._
(2017), CrSiTe3 Williams  _et al._ (2015), FePS3 Lee  _et al._ (2016), NiPS3
Kuo  _et al._ (2016) are promising candidate materials to realize TSF and TSS
phases.

## Acknowledgement

It is a pleasure to thank Oleg Sushkov for useful discussions. Financial
support from the Ministry of Education, Singapore, in the form of grant
MOE2018-T1-001-021 is gratefully acknowledged. A. Chatterjee acknowledges the
financial support received from INSPIRE, Department of Science and Technology,
Govt. of India and the NTU-India Connect Program. A. Chatterjee would like to
thank the School of Physical and Mathematical Sciences at NTU for their kind
hospitality.

##  Appendix A Bond-operator formalism

The spin-operators in terms of bond-operators are given by Romhányi  _et al._
(2011),

$\displaystyle\hat{S}_{j,l}^{+}=\frac{\hat{t}_{1,j}^{\dagger}\hat{t}_{0,j}+\hat{t}_{0,j}^{\dagger}\hat{t}_{\bar{1},j}}{\sqrt{2}}\pm\frac{\hat{s}_{j}^{\dagger}\hat{t}_{\bar{1}j}-\hat{t}_{1,j}^{\dagger}\hat{s}_{j}}{\sqrt{2}}$
$\displaystyle\hat{S}_{j,l}^{-}=\frac{\hat{t}_{\bar{1},j}^{\dagger}\hat{t}_{0,j}+\hat{t}_{0,j}^{\dagger}\hat{t}_{1,j}}{\sqrt{2}}\mp\frac{\hat{s}_{j}^{\dagger}\hat{t}_{1,j}-\hat{t}_{\bar{1},j}^{\dagger}\hat{s}_{j}}{\sqrt{2}}$
$\displaystyle\hat{S}_{j,l}^{z}=\frac{\hat{t}_{1,j}^{\dagger}\hat{t}_{1,j}-\hat{t}_{\bar{1},j}^{\dagger}\hat{t}_{\bar{1},j}}{2}\pm\frac{\hat{s}_{j}^{\dagger}\hat{t}_{0,j}+\hat{t}_{0,j}^{\dagger}\hat{s}_{j}}{2},$
(11)

where the upper-sign is for layer $l=$A and the lower sign is for the layer
$l=$B. $\hat{t}_{1,j}^{\dagger}$, $\hat{t}_{\bar{1},j}^{\dagger}$ and
$\hat{t}_{0,j}^{\dagger}$ create states $\ket{\uparrow\uparrow}$,
$\ket{\downarrow\downarrow}$ and
$(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})/\sqrt{2}$ on $j$-th bond
respectively. After the bond operator transformation, the bilayer honeycomb
spin system transforms into an effective single layer honeycomb lattice system
with triplon and singlet operators on each site. At low temperature, in the
limit $|J_{\perp}|\gg|J|$, it can be assumed that the ground state is product
of singlets on the interlayer NN-bonds. Assuming singlets form the empty or
vacuum state and taking advantage of the hardcore nature of the bosons, we can
simply transform the following quadratic operators into a single triplon
operators in equation Eq. 11,

$\hat{s}_{j}^{\dagger}\hat{t}_{\alpha j}\rightarrow\hat{t}_{\alpha
j},\,\,\hat{t}_{\alpha j}^{\dagger}\hat{s}_{j}\rightarrow\hat{t}_{\alpha
j}^{\dagger}.$

Moreover application of magnetic field in the $z$-direction lowers the energy
of triplon correspond to the operator $\hat{t}^{\dagger}_{1,j}$ and so
neglecting all other triplon operators in equation Eq. 11, we get the
following bond operator transformation Zapf  _et al._ (2014),

$\hat{S}_{j,l}^{+}=\mp\frac{1}{\sqrt{2}}\hat{t}_{j}^{\dagger},\,\hat{S}_{j,l}^{-}=\mp\frac{1}{\sqrt{2}}\hat{t}_{j},\,\hat{S}_{j,l}^{z}=\frac{1}{2}\hat{t}_{j}^{\dagger}\hat{t}_{j},$

where the subscript-$1$ is omitted from the operator
$\hat{t}_{1,j}^{\dagger}$.

##  Appendix B Interaction term in reciprocal space

The last term in equation Eq. 2 represents the NN repulsive interaction and in
reciprocal space it is given by,

$\pazocal{H}_{V}=\frac{J_{z}}{2}\sum_{\bm{\alpha}}\sum_{\mathbb{k},\mathbb{k}^{\prime},\mathbb{k}^{\prime\prime}}e^{i(\mathbb{k}-\mathbb{k}^{\prime})\cdot\bm{\alpha}}\,\,\hat{t}^{\dagger}_{\mathbb{k}}\hat{t}_{\mathbb{k}^{\prime}}\hat{t}^{\dagger}_{\mathbb{k}^{\prime\prime}}\hat{t}_{\mathbb{k}+\mathbb{k}^{\prime\prime}-\mathbb{k}^{\prime}}.$
(12)

Considering triplons are only present at $K$ and $K^{\prime}$ points as well
as neglecting quantum fluctuations around these points, the interaction
Hamiltonian can be explicitly written as,

$\displaystyle\pazocal{H}_{V}$
$\displaystyle=\frac{3J_{z}}{2}\left[\hat{t}^{\dagger}_{\mathbb{K}}\hat{t}_{\mathbb{K}}\hat{t}^{\dagger}_{\mathbb{K}}\hat{t}_{\mathbb{K}}+\hat{t}^{\dagger}_{\mathbb{K}^{\prime}}\hat{t}_{\mathbb{K}^{\prime}}\hat{t}^{\dagger}_{\mathbb{K}^{\prime}}\hat{t}_{\mathbb{K}^{\prime}}+\hat{t}^{\dagger}_{\mathbb{K}}\hat{t}_{\mathbb{K}}\hat{t}^{\dagger}_{\mathbb{K}^{\prime}}\hat{t}_{\mathbb{K}^{\prime}}+\hat{t}^{\dagger}_{\mathbb{K}^{\prime}}\hat{t}_{\mathbb{K}^{\prime}}\hat{t}^{\dagger}_{\mathbb{K}}\hat{t}_{\mathbb{K}}\right]$
$\displaystyle+\frac{J_{z}}{2}\sum_{\bm{\alpha}}\left[e^{i(\mathbb{K}-\mathbb{K}^{\prime})\cdot\bm{\alpha}}\hat{t}^{\dagger}_{\mathbb{K}}\hat{t}_{\mathbb{K}^{\prime}}\hat{t}^{\dagger}_{\mathbb{K}^{\prime}}\hat{t}_{\mathbb{K}}+e^{i(\mathbb{K}^{\prime}-\mathbb{K})\cdot\bm{\alpha}}\hat{t}^{\dagger}_{\mathbb{K}^{\prime}}\hat{t}_{\mathbb{K}}\hat{t}^{\dagger}_{\mathbb{K}}\hat{t}_{\mathbb{K}^{\prime}}\right]$
$\displaystyle=\frac{3J_{z}}{2}\left[\hat{n}_{\mathbb{K}}^{2}+\hat{n}_{\mathbb{K}^{\prime}}^{2}+2\hat{n}_{\mathbb{K}}\hat{n}_{\mathbb{K}^{\prime}}\right]$
$\displaystyle=\frac{3J_{z}}{2}\left(\hat{n}_{\mathbb{K}}+\hat{n}_{\mathbb{K}^{\prime}}\right)^{2},$
(13)

thus coexistence of particles at $\mathbb{K}$ and $\mathbb{K}^{\prime}$ points
does not seem to increase the energy of the system. However, quantum
fluctuations around the points $\mathbb{K}$ and $\mathbb{K}^{\prime}$
increases the energy due to presence of triplons at both $\mathbb{K}$ and
$\mathbb{K}^{\prime}$ points as shown in the reference Ref. Li  _et al._
(2014).

##  Appendix C Introducing more terms in the spin Hamiltonian

In this appendix, we have taken into account additional interlayer bonds (see
Fig. 5(a), (b)). We introduce intralayer NNN Heisenberg-exchange interaction
$J_{1}$ as well as interlayer NN Heisenberg exchange interaction $J_{2\perp}$
in the spin Hamiltonian in equation Eq. 1,

$\displaystyle\pazocal{H}=$ $\displaystyle
J_{\perp}\sum_{\begin{subarray}{c}i,\,m\in A\\\ n\in
B\end{subarray}}\mathbb{S}_{i,m}\cdot\mathbb{S}_{i,n}-B_{z}\sum_{i,m}S^{z}_{i,m}+\sum_{\left\langle
i,j\right\rangle,m}\left[J(\mathbb{S}^{x}_{i,m}\mathbb{S}^{x}_{j,m}+\mathbb{S}^{y}_{i,m}\mathbb{S}^{y}_{j,m})+J_{z}\mathbb{S}^{z}_{i,m}\mathbb{S}^{z}_{j,m}\right]$
$\displaystyle+D\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle,m}\nu_{ij}\hat{z}\cdot(\mathbb{S}_{i,m}\times\mathbb{S}_{j,m})+D_{\perp}\sum_{\begin{subarray}{c}\left\langle\left\langle
i,j\right\rangle\right\rangle\\\ m\in A,B,\,n\neq
m\end{subarray}}\nu_{ij}\hat{z}\cdot(\mathbb{S}_{i,m}\times\mathbb{S}_{j,n})$
$\displaystyle+J_{1}\sum_{\left\langle\left\langle
i,j\right\rangle\right\rangle,m}\mathbb{S}_{i,m}\cdot\mathbb{S}_{i,m}+J_{2,\perp}\sum_{\begin{subarray}{c}\left\langle
i,j\right\rangle\\\ m\in A,B,\,n\neq
m\end{subarray}}\mathbb{S}_{i,m}\cdot\mathbb{S}_{i,n}+J_{3,\perp}\sum_{\begin{subarray}{c}\left\langle\left\langle
i,j\right\rangle\right\rangle\\\ m\in A,B,\,n\neq
m\end{subarray}}\mathbb{S}_{i,m}\cdot\mathbb{S}_{i,n}.$ (14)

Moreover we have added symmetry allowed interlayer NNN DMI $D_{\perp}$ and
Heisenberg interaction $J_{3\perp}$. The DMI on dimer-bond and NN interlayer
bonds are zero due to presence of inversion center at the middle of the bonds.
The corresponding real space triplon Hamiltonian is given as,

$\displaystyle\pazocal{H}=$
$\displaystyle\frac{J-J_{2\perp}}{2}\sum_{\left\langle
ij\right\rangle}\left[\hat{t}^{\dagger}_{i}\hat{t}_{j}+\text{H.c.}\right]+\frac{J_{1}-J_{3\perp}}{4}\sum_{\left\langle\left\langle
ij\right\rangle\right\rangle}\left[\hat{t}^{\dagger}_{i}\hat{t}_{j}+\text{H.c.}\right]+\frac{\mathfrak{i}(D+D_{\perp})}{2}\sum_{\left\langle\left\langle
ij\right\rangle\right\rangle}\nu_{ij}\left[\hat{t}^{\dagger}_{i}\hat{t}_{j}-\text{H.c.}\right]$
$\displaystyle+\left(\frac{J_{\perp}}{4}-B_{z}\right)\sum_{i}\hat{t}^{\dagger}_{i}\hat{t}_{i}+\frac{J_{2\perp}+J_{z}}{2}\sum_{\left\langle
ij\right\rangle}\hat{n}_{i}\hat{n}_{j}+\frac{J_{1}+J_{3\perp}}{4}\sum_{\left\langle\left\langle
ij\right\rangle\right\rangle}\hat{n}_{i}\hat{n}_{j}$ (15)

It is noticeable that the Heisenberg interaction $J_{2\perp}$ renormalizes the
interactions $J$ and $J_{z}$, whereas $J_{1}$ and $J_{3\perp}$ add additional
NNN hopping term and NNN interaction term. Moreover $D_{\perp}$ renormalizes
$D$ by simply adding up and so $D_{\perp}$ is absorbed into $D$ in the rest of
this section. In this appendix, we investigate the triplon Hamiltonian,
neglecting the interaction terms. The non-interacting Hamiltonian in $k$-space
is similar to equation Eq. 3 and given by,

$\pazocal{H}_{0}=\sum_{\mathbb{k}}\Psi_{\mathbb{k}}^{\dagger}\left[g^{\prime}(\mathbb{k})\sigma_{0}+\mathbb{h^{\prime}}\cdot\bm{\sigma}\right]\Psi_{\mathbb{k}},$
(16)

where,

$\displaystyle g^{\prime}(\mathbb{k})$
$\displaystyle=(J_{\perp}/4)-B_{z}+((J_{1}-J_{3\perp})/2)\sum_{i}\cos({\mathbb{k}\cdot\bm{\beta}_{i}}),$
$\displaystyle h^{\prime}_{x}(\mathbb{k})$
$\displaystyle=((J-J_{2\perp})/2)\sum_{i}\cos(\mathbb{k}\cdot\bm{\alpha}_{i}),$
$\displaystyle h^{\prime}_{y}(\mathbb{k})$
$\displaystyle=((J-J_{2\perp})/2)\sum_{i}\sin(\mathbb{k}\cdot\bm{\alpha}_{i}),$
$\displaystyle h^{\prime}_{z}(\mathbb{k})$
$\displaystyle=D\sum_{i}\sin(\mathbb{k}\cdot\bm{\beta}_{i})$

, where $\bm{\alpha}_{i}$ and $\bm{\beta}_{i}$ are the NN and NNN vectors for
each layer respectively. The band structure for two different $D$-values are
plotted in the figure Fig. 5(c) and (d). The minima of the bands are at the
$\Gamma$-point and $K$ (or $K^{\prime}$)-point for $D=0.1J$ and $D=0.8J$
respectively. The results are same when $J_{2\perp}=0$ and $J_{1}=0$ as in the
figures Fig. 1(c) and (d). The condition for the band minima at $K$ or
$K^{\prime}$-point is given as,

$|D|>\frac{2|J-J_{2\perp}|-3(J_{1}-J_{3\perp})}{2\sqrt{3}}.$ (17)

We note that the critical DMI required to realise TSF or TSS phases depends on
various Heisenberg exchange interactions. The conditions $J\approx J_{2\perp}$
and $J_{1}\approx J_{3\perp}$ lowers the critical value of the DMI to realize
TSF or TSS phases. Thus we can conclude that frustration among interlayer and
intralayer Heisenberg interactions can lower the critical DMI required to
realize TSF and TSS phases.

Figure 5:  (Color online) (a) Intradimer and interlayer NN bonds are indicated
by red solid and red dotted lines, respectively. (b) The intralayer and
interlayer NNN bonds are indicated by green solid and green dotted lines,
respectively. (c) Triplon band structure at $D=0.1J$, (b) Triplon band
structure at $D=0.8J$. The other parameters for the band structures are
$J_{\perp}=10J$, $J_{1}=0.8J$, $J_{2\perp}=3J$, $B_{z}=0.0$, $J_{z}=0.0$,
$J_{3\perp}=0$.

## References

  * Soltan-Panahi  _et al._ (2011) P. Soltan-Panahi, D.-S. Lühmann, J. Struck, P. Windpassinger, and K. Sengstock, Nature Physics 8, 71–75 (2011).
  * Luke  _et al._ (1998) G. M. Luke, Y. Fudamoto, K. M. Kojima, M. I. Larkin, J. Merrin, B. Nachumi, Y. J. Uemura, Y. Maeno, Z. Q. Mao, Y. Mori, and et al., Nature 394, 558–561 (1998).
  * Wysokinski (2019) K. Wysokinski, Condensed Matter 4, 47 (2019).
  * Luke  _et al._ (1993) G. M. Luke, A. Keren, L. P. Le, W. D. Wu, Y. J. Uemura, D. A. Bonn, L. Taillefer, and J. D. Garrett, Phys. Rev. Lett. 71, 1466 (1993).
  * Kaur and Agterberg (2003) R. P. Kaur and D. F. Agterberg, Phys. Rev. B 68, 100506(R) (2003).
  * Kaminski  _et al._ (2002) A. Kaminski, S. Rosenkranz, H. Fretwell, J. C. Campuzano, Z. Li, H. Raffy, W. Cullen, H. You, C. Olson, C. Varma, and H. Höchst, Nature 416, 610 (2002).
  * Jürgensen  _et al._ (2015) O. Jürgensen, K. Sengstock, and D.-S. Lühmann, Scientific Reports 5, 12912 (2015).
  * Tanaka  _et al._ (2001) H. Tanaka, A. Oosawa, T. Kato, H. Uekusa, Y. Ohashi, K. Kakurai, and A. Hoser, Journal of the Physical Society of Japan 70, 939 (2001).
  * Coldea  _et al._ (2003) R. Coldea, D. A. Tennant, and Z. Tylczynski, Phys. Rev. B 68, 134424 (2003).
  * Singh and Johnston (2007) Y. Singh and D. C. Johnston, Phys. Rev. B 76, 012407 (2007).
  * Demokritov  _et al._ (2006) S. Demokritov, V. Demidov, O. Dzyapko, G. Melkov, A. Serga, B. Hillebrands, and A. Slavin, Nature 443, 430 (2006).
  * Dzyapko  _et al._ (2007) O. Dzyapko, V. Demidov, S. Demokritov, G. Melkov, and A. Slavin, New Journal of Physics 9, 64 (2007).
  * Bradley  _et al._ (1995) C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * Zapf  _et al._ (2014) V. Zapf, M. Jaime, and C. D. Batista, Rev. Mod. Phys. 86, 563 (2014).
  * Sachdev and Bhatt (1990) S. Sachdev and R. N. Bhatt, Phys. Rev. B 41, 9323 (1990).
  * Romhányi  _et al._ (2011) J. Romhányi, K. Totsuka, and K. Penc, Phys. Rev. B 83, 024413 (2011).
  * Li  _et al._ (2014) X. Li, S. S. Natu, A. Paramekanti, and S. D. Sarma, Nature Communications 5, 5174 (2014).
  * Adhikari  _et al._ (2010) S. K. Adhikari, B. A. Malomed, L. Salasnich, and F. Toigo, Phys. Rev. A 81, 053630 (2010).
  * Henley (1989) C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * Giamarchi  _et al._ (2008) T. Giamarchi, C. Rüegg, and O. Tchernyshyov, Nature Physics 4, 198–204 (2008).
  * Lühmann (2013) D.-S. Lühmann, Phys. Rev. A 87, 043619 (2013).
  * Chen and Yang (2017) Y.-C. Chen and M.-F. Yang, Journal of Physics Communications 1, 035009 (2017).
  * Gelfand  _et al._ (1989) M. P. Gelfand, R. R. P. Singh, and D. A. Huse, Phys. Rev. B 40, 10801 (1989).
  * Yamamoto  _et al._ (2014) D. Yamamoto, G. Marmorini, and I. Danshita, Phys. Rev. Lett. 112, 127203 (2014).
  * Moreno-Cardoner  _et al._ (2014) M. Moreno-Cardoner, H. Perrin, S. Paganelli, G. De Chiara, and A. Sanpera, Phys. Rev. B 90, 144409 (2014).
  * Yamamoto  _et al._ (2012a) D. Yamamoto, I. Danshita, and C. A. R. Sá de Melo, Phys. Rev. A 85, 021601(R) (2012a).
  * Yamamoto  _et al._ (2012b) D. Yamamoto, A. Masaki, and I. Danshita, Phys. Rev. B 86, 054516 (2012b).
  * Singh  _et al._ (2014) M. Singh, T. Mishra, R. V. Pai, and B. P. Das, Phys. Rev. A 90, 013625 (2014).
  * Hassan  _et al._ (2007) S. R. Hassan, L. de Medici, and A. M. S. Tremblay, Phys. Rev. B 76, 144420 (2007).
  * McIntosh  _et al._ (2012) T. McIntosh, P. Pisarski, R. J. Gooding, and E. Zaremba, Phys. Rev. A 86, 013623 (2012).
  * Isaev  _et al._ (2009a) L. Isaev, G. Ortiz, and J. Dukelsky, Phys. Rev. B 79, 024409 (2009a).
  * Isaev  _et al._ (2009b) L. Isaev, G. Ortiz, and J. Dukelsky, Phys. Rev. Lett. 103, 177201 (2009b).
  * Huerga  _et al._ (2013) D. Huerga, J. Dukelsky, and G. E. Scuseria, Phys. Rev. Lett. 111, 045701 (2013).
  * Yamamoto  _et al._ (2017) D. Yamamoto, H. Ueda, I. Danshita, G. Marmorini, T. Momoi, and T. Shimokawa, Phys. Rev. B 96, 014431 (2017).
  * Huerga  _et al._ (2014) D. Huerga, J. Dukelsky, N. Laflorencie, and G. Ortiz, Phys. Rev. B 89, 094401 (2014).
  * Choudhury and Mueller (2013) S. Choudhury and E. J. Mueller, Phys. Rev. A 87, 033621 (2013).
  * Cao  _et al._ (2015) L. Cao, S. Krönke, J. Stockhofe, J. Simonet, K. Sengstock, D.-S. Lühmann, and P. Schmelcher, Phys. Rev. A 91, 043639 (2015).
  * Lühmann (2016) D.-S. Lühmann, Phys. Rev. A 94, 011603(R) (2016).
  * Gan  _et al._ (2007) J. Y. Gan, Y. C. Wen, J. Ye, T. Li, S.-J. Yang, and Y. Yu, Phys. Rev. B 75, 214509 (2007).
  * Zhang  _et al._ (2016) H. Zhang, C. A. Lamas, M. Arlego, and W. Brenig, Physical Review B 93, 235150 (2016).
  * Kandpal and van den Brink (2011) H. C. Kandpal and J. van den Brink, Phys. Rev. B 83, 140412(R) (2011).
  * Zhang  _et al._ (2014) H. Zhang, M. Arlego, and C. A. Lamas, Phys. Rev. B 89, 024403 (2014).
  * Oitmaa and Singh (2012) J. Oitmaa and R. R. P. Singh, Phys. Rev. B 85, 014428 (2012).
  * Ganesh  _et al._ (2011a) R. Ganesh, S. V. Isakov, and A. Paramekanti, Phys. Rev. B 84, 214412 (2011a).
  * Mulder  _et al._ (2010) A. Mulder, R. Ganesh, L. Capriotti, and A. Paramekanti, Phys. Rev. B 81, 214419 (2010).
  * Wang (2010) F. Wang, Phys. Rev. B 82, 024419 (2010).
  * Okumura  _et al._ (2010) S. Okumura, H. Kawamura, T. Okubo, and Y. Motome, Journal of the Physical Society of Japan 79, 114705 (2010).
  * Clark  _et al._ (2011) B. K. Clark, D. A. Abanin, and S. L. Sondhi, Phys. Rev. Lett. 107, 087204 (2011).
  * Cabra  _et al._ (2011) D. C. Cabra, C. A. Lamas, and H. D. Rosales, Phys. Rev. B 83, 094506 (2011).
  * Ganesh  _et al._ (2011b) R. Ganesh, D. N. Sheng, Y.-J. Kim, and A. Paramekanti, Phys. Rev. B 83, 144414 (2011b).
  * Albuquerque  _et al._ (2011) A. F. Albuquerque, D. Schwandt, B. Hetényi, S. Capponi, M. Mambrini, and A. M. Läuchli, Phys. Rev. B 84, 024406 (2011).
  * CABRA  _et al._ (2011) D. C. CABRA, C. A. LAMAS, and H. D. ROSALES, Modern Physics Letters B 25, 891 (2011).
  * Mezzacapo and Boninsegni (2012) F. Mezzacapo and M. Boninsegni, Phys. Rev. B 85, 060402(R) (2012).
  * Li  _et al._ (2012) P. H. Y. Li, R. F. Bishop, D. J. J. Farnell, and C. E. Campbell, Phys. Rev. B 86, 144404 (2012).
  * Ganesh  _et al._ (2013) R. Ganesh, J. van den Brink, and S. Nishimoto, Phys. Rev. Lett. 110, 127203 (2013).
  * Gong  _et al._ (2013) S.-S. Gong, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Phys. Rev. B 88, 165138 (2013).
  * Huang  _et al._ (2017) B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, and et al., Nature 546, 270–273 (2017).
  * Li  _et al._ (2019) T. Li, S. Jiang, N. Sivadas, Z. Wang, Y. Xu, D. Weber, J. E. Goldberger, K. Watanabe, T. Taniguchi, C. J. Fennie, and et al., Nature Materials 18, 1303–1308 (2019).
  * Samuelsen  _et al._ (1971) E. J. Samuelsen, R. Silberglitt, G. Shirane, and J. P. Remeika, Phys. Rev. B 3, 157 (1971).
  * Gong  _et al._ (2017) C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, and et al., Nature 546, 265–269 (2017).
  * Williams  _et al._ (2015) T. J. Williams, A. A. Aczel, M. D. Lumsden, S. E. Nagler, M. B. Stone, J.-Q. Yan, and D. Mandrus, Phys. Rev. B 92, 144404 (2015).
  * Lee  _et al._ (2016) J.-U. Lee, S. Lee, J. H. Ryoo, S. Kang, T. Y. Kim, P. Kim, C.-H. Park, J.-G. Park, and H. Cheong, Nano Letters 16, 7433 (2016).
  * Kuo  _et al._ (2016) C.-T. Kuo, M. Neumann, B. Karuppannan, H. Park, S. Kang, H. Shiu, J. Kang, B. Hong, M. Han, T. Noh, and J.-G. Park, Scientific Reports 6, 20904 (2016).
  * Okubo  _et al._ (2012) S. Okubo, T. Ueda, H. Ohta, W. Zhang, T. Sakurai, N. Onishi, M. Azuma, Y. Shimakawa, H. Nakano, and T. Sakai, Phys. Rev. B 86, 140401(R) (2012).
  * Chen  _et al._ (2018) L. Chen, J.-H. Chung, B. Gao, T. Chen, M. B. Stone, A. I. Kolesnikov, Q. Huang, and P. Dai, Phys. Rev. X 8, 041028 (2018).
  * Vinas Boström  _et al._ (2020) E. Vinas Boström, M. Claassen, J. McIver, G. Jotzu, A. Rubio, and M. Sentef, SciPost Physics 9, 061 (2020).
  * Owerre (2017) S. A. Owerre, Journal of Physics Communications 1, 021002 (2017).
  * Tacchi  _et al._ (2017) S. Tacchi, R. E. Troncoso, M. Ahlberg, G. Gubbiotti, M. Madami, J. Åkerman, and P. Landeros, Phys. Rev. Lett. 118, 147201 (2017).
  * Heide  _et al._ (2008) M. Heide, G. Bihlmayer, and S. Blügel, Phys. Rev. B 78, 140403(R) (2008).
  * Di  _et al._ (2015) K. Di, V. L. Zhang, H. S. Lim, S. C. Ng, M. H. Kuok, J. Yu, J. Yoon, X. Qiu, and H. Yang, Phys. Rev. Lett. 114, 047201 (2015).
  * Gastaldo  _et al._ (2019) D. Gastaldo, N. Strelkov, L. D. Buda-Prejbeanu, B. Dieny, O. Boulle, P. Allia, and P. Tiberto, Journal of Applied Physics 126, 103905 (2019).