\section{Introduction}
The Heisenberg model is arguably the most important construct in quantum magnetism. Experimentally, it is usually studied in transition metal oxides and metal-organic salts.
Unfortunately, symmetry-breaking magnetic anisotropy is unavoidable in such compounds. Even when weak compared to the dominant Heisenberg exchange, it can have a significant influence on the magnetic and dynamical properties.
In understanding these, existing microscopic theories are
usually disregarded in favor of minimalistic empirical models.
For example, magnetic exchange interactions in insulators are described within Anderson's theory of super exchange\cite{Anderson1959,Moriya1960PRL,Moriya1960}. In the presence of weak spin-orbit-coupling it predicts both symmetric (Ising) and antisymmetric (Dzyaloshinskii-Moriya, DM) exchange anisotropy always occurring in combination and at a fixed ratio\cite{Kaplan1983,SEWA1992,SEWA1993}.
Yet, experimental data is often successfully analyzed in comparison to a Heisenberg model with only an added DM term  \cite{Oshikawa1997,Affleck1999,Zvyagin2004,Coldea2002,Haelg2014} even though such a model lacks microscopic justification.
Even more often, signatures of weak magnetic anisotropy are detected in experiments, but their origin remains unknown altogether \cite{Kenzelmann2001,Jacobsen2018,Bettler2019}.
In short, anisotropy is ubiquitously present, yet pinpointing its microscopic origin is almost never possible.

Here we present a detailed case study: the effect of very weak magnetic anisotropy on spin excitations in
the almost ideal Heisenberg quantum spin ladder compound \BPCB{} (BPCB)\cite{Patyal1990,Savici2009,Ruegg2008BPCBThermodynamics,Klanjsec2008NMROrdering,ThielemannRuegg2009,CizmarBPCBAniso2010,Blosser2018}.
On this compound we have obtained exceptionally high resolution inelastic neutron scattering data showing a small but distinct anisotropy splitting of the triplet excitations as described in Sec.~\ref{sec:Experiment}.
At the same time, BPCB is an excellent realization of a strong rung spin ladder. It is described by only two exchange couplings: the rung exchange and an almost four times smaller leg coupling. For this model, the excitation spectrum in the presence of various possible exchange anisotropies can be precisely calculated as detailed in Sec.~\ref{sec:TripletDispersion}.
This combination of high resolution spectroscopic data and precise calculations of the excitation spectrum allows a detailed discussion of the anisotropy's microscopic origins in BPCB (Sec.~\ref{sec:Discussion}). These findings we compare to the microscopic theory of anisotropic superexchange and we discuss how the present case study may guide our understanding of magnetic anisotropy in related compounds.

\section{Experiment and results}\label{sec:Experiment}

\begin{figure}
\includegraphics{Fig_StructureAniso2.pdf}
\caption{\label{fig:structure}
Schematic of the crystal structure of BPCB. The spin ladders are formed by the magnetic Cu${^{2+}}$ cations and linking Br${^{-}}$ anions. Crystallographic centers of inversion symmetry are marked by green dots. Anisotropic Dzyaloshinskii-Moriya interactions, by symmetry are only allowed on the ladder legs and possible DM vectors $\mb{D}$ are sketched as blue arrows on these bonds.
Note that for this monoclinic structure, in both plots the third axis (projection axis) is not perpendicular to the plane of the figure.}
\end{figure}

\begin{figure*}[t]
\includegraphics{Fig_Spectra3.pdf}
\caption{\label{fig:spectra}False-color maps of the inelastic neutron scattering intensity measured using a) $E_\mr{i}=2.20$~meV and b) $E_\mr{i}=1.35$~meV incident energy neutrons at $T=0.35$~K. With the lower incident energy, a superior energy resolution is achieved but only the low energy parts of the dispersion can be probed.}
\end{figure*}

\subsection{The compound \BPCB: crystal structure and magnetic interaction pathways}
The compound \BPCB, bis-piperidinium copper bromide or BPCB for short crystallizes in a monoclinic crystal structure with space group P2$_1$/c and room temperature lattice parameters $a=8.49$, $b=17.22$, $c=12.38$~\AA{}, and $\beta=99.3^{\circ}$ \cite{Patyal1990}.
The spin ladders are formed by the magnetic Cu$^{2+}$ cations carrying a spin $S=1/2$ linked by super-exchange bridges via Br$^-$ anions \cite{Savici2009}.
The ladders run along the crystallographic $a$ axis and are well separated by non-magnetic organic piperidinium ions as depicted in Fig.~\ref{fig:structure}\footnote{The crystal structure plots shown in Fig.~\ref{fig:structure} are based on renderings generated by the VESTA software package\cite{VESTA2011}.}

There are two crystallographically equivalent spin ladders in BPCB, related by the glide plane symmetry of the P2$_1$/c space group. Furthermore, the center of each ladder rung corresponds to a crystallographic center of inversion symmetry. For this reason, on the ladder rungs, only symmetric exchange is possible\cite{DZYALOSHINSKY1958}, whilst for the ladder legs there are no symmetry restrictions and both symmetric and antisymmetric exchange are in principle allowed. On the ladder leg, the antisymmetric exchange contribution is parameterized by a Dzyaloshinskii-Moriya vector $\mb{D}$. These DM vectors are uniform within every ladder leg and anti-aligned between the two legs of each ladder. The DM vectors of the two crystallographically equivalent ladders are again related by the reflection of the glide plane symmetry as sketched in Fig.~\ref{fig:structure}. Besides these relations, the DM vectors may point in any direction and are not further constrained by symmetry.
Similarly, The $g$-tensors for the magnetic moments residing on the different Cu$^{2+}$ ions show different principal axes for the two types of ladders.
Thus, in an applied magnetic field, the two types of ladders become inequivalent, except for special orientations of the magnetic field \cite{ThielemannRuegg2009,CizmarBPCBAniso2010}. In the present study, however, we exclusively focus on zero-field properties where the two ladders remain equivalent.

\subsection{Inelastic neutron scattering}

Neutron scattering experiments were performed on four fully deuterated single crystals of total mass 2.07~g, co-aligned to better than 1$^\circ$ effective mosaic spread.
The measurements were carried out at the LET cold neutron time-of-flight spectrometer \cite{Bewley2011} at the ISIS facility, UK.
The sample was mounted on a $^3$He-$^4$He dilution refrigerator
with the crystallographic $b$ axis vertical.
Making use of repetition rate multiplication, data was collected simultaneously using neutrons of incident energies $E_\mr{i}=$ 1.35, 2.20, 4.20 and 11.0~meV. For these configurations we find an approximately Gaussian energy resolution of 19, 36, 97 and 410~$\mu$eV at FWHM, respectively, for elastic scattering ($\hbar\omega=0$) and improving towards higher energy transfer.

\begin{figure}
	\includegraphics{Fig_Cuts2.pdf}
	\caption{\label{fig:cuts}
		Cuts through the data shown in Fig.~\ref{fig:spectra} at the band extrema: a) Cut at the band minimum $q_\|/2\pi=-0.5$ through the $1.35$~meV data set. b,c) Cuts at the band extrema $q_\|/2\pi=-0.5$ or $q_\|/2\pi=0$, respectively, through the $2.20$~meV data set.
		Solid lines correspond to Gaussian fits. At the band minimum we find a slightly larger splitting of $\Delta_\mr{min}=50(1)$~$\mu$eV than at the band maximum where $\Delta_\mr{max}=40(2)$~$\mu$eV.			}
\end{figure}

\subsection{Experimental results}

An overview of the neutron scattering data\footnote{The raw data is available in the STFC Research Data Repository \cite{DataISIS_BPCB}.}
collected for $E_\mathrm{i}=2.2$ and 1.35 meV at a temperature of $0.35$~K is presented as false color intensity plots in Fig.~\ref{fig:spectra}.
For this quasi-one-dimensional spin ladder system, the measured scattering intensity is plotted versus reduced wave-vector transfer parallel to the ladder $q_\|=\mb{Q}\cdot\mb{a}$ and integrated fully along the non-dispersive $b^*$ and $c^*$ directions.

We observe a band of excitations with an excitation gap of $\Delta\sim 0.8$~meV and a bandwidth of $\Gamma\sim 0.6$~meV. The line-width of gapped excitations in 1D systems is known to show activation behavior\cite{Zheludev2008}. With $k_\mr{B}T\ll\Delta$, in the present experiment these excitations are long lived and their apparent width corresponds to the experimental resolution.
Cuts through both data sets at the band extrema are shown in Fig.~\ref{fig:cuts}. These correspond to the data of Fig.~\ref{fig:spectra} integrated in thin slices of $q_\|/2\pi=-0.5\pm 0.02$ and $q_\|/2\pi=0\pm 0.03$, respectively (the dispersion is less curved at the maximum allowing to integrate a wider slice for improved counting statistics).
To these cuts we fit two Gaussian peaks of equal width and a linear background shown as a solid lines in Fig.~\ref{fig:cuts}. We find a splitting of $\Delta_\mr{min}=50(1)$~$\mu$eV at the band minimum and $\Delta_\mr{max}=40(2)$~$\mu$eV at the band maximum.

\begin{widetext}
\section{Triplet dispersion for a spin ladder with exchange anisotropy}
\label{sec:TripletDispersion}

\subsection{Isotropic Heisenberg ladder}
\label{sec:LadderHamiltonian_Heis}
For an isotropic spin $S=1/2$ ladder we denote the rung and leg Heisenberg exchange constants as $J_\bot$ and $J_\|$, respectively. In the limit of $\lambda=J_\|/J_\bot\ll1$, the spin ladder can be understood as a 1D array of weakly interacting spin dimers. For $\lambda=0$ the groundstate corresponds to a product state of singlets and the low energy excitations are triplet excitations. In the presence of weak leg coupling, these excitations become mobile and for small $\lambda$ their dispersion can be calculated in a strong coupling expansion\cite{ReigrotzkiRice1994}.
For an isotropic ladder with Hamiltonian
\begin{equation}
\mc{H} = J_\bot \left( \mc{H}^\bot + \lambda \mc{H}^\| \right), \qquad
\mc{H}^\bot = \sum_{R=1}^{L}  \mb{S}_{R,1} \cdot \mb{S}_{R,2}, \qquad
\mc{H}^\| = \sum_{R=1}^{L} \sum_{\alpha=1}^{2}  \mb{S}_{R,\alpha} \cdot \mb{S}_{R+1,\alpha}, \qquad
\end{equation}
up to 3\textsuperscript{rd} order in $\lambda$, the following dispersion is obtained in Ref.~\onlinecite{ReigrotzkiRice1994},
degenerate for the three $\sigma=\{+,0,-\}$ triplet branches:
\begin{equation}
\epsilon_\sigma(k)/J_\bot =   1 + \gamma \cos(k)
+ \frac{\gamma^2}{4} \Big[3-\cos(2k)\Big]
+ \frac{\gamma^3}{8}\Big[3-2\cos(k)-2\cos(2k)+\cos(3k) \Big].
\label{eq:bpcb_TripletDispersion}
\end{equation}
{In the following we do the same calculation for a ladder with anisotropic rung or leg exchange, respectively. Some details of these computations are given in the appendix.}

\subsection{Anisotropic rung interactions}
\label{sec:LadderHamiltonian_RungAniso}
In addition to the isotropic rung and leg couplings, here Ising-type anisotropy on the ladder rungs is considered. We choose the $z$ axis parallel the Ising axis and consider the Hamiltonian
\begin{equation}
\mc{H} = J_\bot \left( \mc{H}^\bot + \lambda \mc{H}' \right), \qquad
\mc{H}' = \mc{H}^\| + C\, \mc{H}^\mr{Rung,Ising}, \qquad
\mc{H}^\text{Rung,Ising} = \sum_{R=1}^{L}  {S}_{R,1}^z {S}_{R,2}^z,
\end{equation}
where the Ising anisotropy on the ladder rungs is parameterized by $C$. Using the same strong coupling expansion as applied to the isotropic case in Ref.~\onlinecite{ReigrotzkiRice1994}, up to 3\textsuperscript{rd} order in $\lambda$ we find
\begin{align}
	\epsilon_0(k)/J_\bot &=
	1+ \lambda \cos(k)
	+\frac{\lambda^2}{4} \Big[3-\cos(2k)\Big]
	+\frac{\lambda^3}{8} \Big[3-2C  -2\cos(k)  -2\cos(2k)+\cos(3k)\Big]\\
	%
	\epsilon_\pm(k)/J_\bot &=
	1 +\lambda \Big[\frac{C}{2}+\cos(k)\Big]
	+\frac{\lambda^2}{4} \Big[3-\cos(2k)\Big]
	+\frac{\lambda^3}{8} \Big[3 -2C -2\cos(k) +\big(C-2\big)\cos(2k)+\cos(3k)\Big].
\end{align}

\subsection{Anisotropic leg interactions}
\label{sec:LadderHamiltonian_LegAniso}
Finally, we consider anisotropic exchange on the ladder legs with both symmetric and anti-symmetric contributions. Assuming a center of inversion symmetry on the ladder rungs, we consider a uniform DM vector $+\mb{D}$ on one leg and $-\mb{D}$ on the other leg. Furthermore we include Ising anisotropy of independent magnitude but pointing in the same direction as $\mb{D}$. We choose the $z$ axis parallel to $\mb{D}$ and consider
\begin{align}\label{eq:LadderHamiltonian_LegAniso}
	\mc{H} = J_\bot \left( \mc{H}^\bot + \lambda \mc{H}' \right), \quad\quad
	\mc{H}' = \mc{H}^\| + &B\, \mc{H}^\mr{Leg,DM} + A\, \mc{H}^\mr{Leg,Ising}\\
	%
	\mc{H}^\text{Leg,DM} = \sum_{R=1}^{L} \sum_{\alpha=1}^{2}  (-1)^\alpha \mb{e}_z\cdot(\mb{S}_{R,\alpha} \times \mb{S}_{R+1,\alpha}), \quad&\quad
	%
	\mc{H}^\text{Leg,Ising} =\sum_{R=1}^{L} \sum_{\alpha=1}^{2}
	S_{R,\alpha}^z S_{R+1,\alpha}^z
\end{align}
Up to 3\textsuperscript{rd} order in $\lambda$ we find the following triplet dispersion:
\begin{align}
	\epsilon_0(k)/J_\bot =&
	1 + \lambda (1+A)\cos(k)
	+
	\frac{\lambda^2}{4} \Big[3 -4B^2 + 2A + A^2 + 4B^2 \cos(k) - \Big(1 + 2A + A^2\Big)\cos(2 k) \Big] \notag\\
	&+
	\frac{\lambda^3}{8} \Big[3 +3A - \Big(2 +4B^2 +4A +4AB^2 +3A^2 + A^3\Big)\cos(k)   \notag\\&
	  \qquad\quad -
	\Big(2 +2A -4B^2 -4AB^2 \Big)\cos(2k)
	+ \Big(1 + 3A + 3A^2 + A^3\Big) \cos(3k) \Big]   \\
	%
	%
	%
	\epsilon_\pm(k)/J_\bot =& 1 + \lambda \cos(k)
	+
	\frac{\lambda^2}{4} \Big[3 -2B^2 +2A + A^2 - \cos(2k) \Big] \notag\\
	&+
	\frac{\lambda^3}{8} \Big[3 +3A - \Big(2 +12B^2 +A +\frac{A^2}{2} \Big)\cos(k)
	-
	\Big(2 + 2B^2 +2A +2AB^2\Big)\cos(2k)
	+ \frac{1}{8}\cos(3k) \Big]
\end{align}

\end{widetext}

\section{Discussion}\label{sec:Discussion}

\begin{figure*}
	\includegraphics{DispAniso.pdf}
	\caption{\label{fig:DispAniso}
		Triplet dispersion for a spin ladder with exchange anisotropy: a) Ising anisotropy on the ladder rungs, b) Ising anisotropy on the legs, c) DM anisotropy on the legs, and d) both Ising and DM anisotropy on the ladder legs, in the ratio as predicted for a simple super-exchange mechanism (Eq.~\ref{eq:SingleBondExchangeAniso}). For these plots, isotropic $J_\bot/k_\mr{B}=12.7$~K and $J_\|/k_\mr{B}=3.54$~K were used, approximately describing the ladders in BPCB. Further, in all four cases the anisotropy parameter was chosen such that the triplet splitting at the band minimum amounts to 50~$\mu$eV. In all plots, the blue line shows the doubly degenerate $\epsilon_\pm(k)$ dispersion and the orange line shows the $\epsilon_0(k)$ triplet. }
\end{figure*}

\subsection{Discussion of experimental results}
The experimental data shown in Fig.~\ref{fig:spectra} closely resembles previously published results on BPCB with dispersive triplet excitations \cite{Savici2009,ThielemannRuegg2009}. In these experiments the resolution was not sufficient to detect the small band splitting and those data were very well described using a model of an isotropic spin ladder with $J_\bot/k_\mr{B}=12.7(1)$~K and $J_\|/k_\mr{B}=3.54(3)$~K\cite{Blosser2018}.
Only with the superior resolution of the present experiment, the small splitting of the triplet band is observed.

In our data, perpendicular to the one-dimensional ladders, no dispersion is observed within the experimental resolution. This validates BPCB as an exceptionally one-dimensional ladder compound and it corroborates, that the observed band splitting is inherent to the spin ladders and is not related to residual 3D interactions. In addition, inter-ladder interactions were previously estimated on the order of a few 10s of mK\cite{Klanjsec2008NMROrdering} (a few $\mu$eV) -- much smaller than the observed band splitting.

A previous ESR study on BPCB\cite{CizmarBPCBAniso2010}
found clear signatures of magnetic anisotropy, of exactly the same magnitude as the band splitting observed in the present experiments.
By a careful study of the angular dependence of the ESR signal, an anisotropy axis tilted approx.~50$^\circ$ from the $b$ axis in the $(b,c^\star)$ plane could be identified as a special direction. However, the origin of this anisotropy remained unclear.

One source of magnetic anisotropy are dipolar interactions present in all materials\cite{Fazekas1999}. For classical magnetic moments at the shortest Cu--Cu distance found in BPCB, dipolar interactions are estimated as $J_\mr{dipolar}/k_\mr{B}\approx 6$~mK ($< 1$~$\mu$eV), much smaller than the observed band splitting. Dipolar interactions clearly can be ruled out as the primary cause of magnetic anisotropy in BPCB. Since for the $S=1/2$ Cu$^{2+}$ ions there cannot be any single-ion anisotropy, \emph{exchange anisotropy} must be the dominant cause of the observed small splitting of the triplet band. Before further analyzing our data in this regard, we first discuss the nature of anisotropic super-exchange interactions in the following section.

\subsection{Anisotropic super-exchange}
Most generally, the effective exchange interaction between two spins takes the form $\mc{H}_{1,2}=\mb{S}_1 \Gamma \mb{S}_2$ where $\Gamma$ is a tensor. Decomposing $\Gamma$ into its symmetric and anti-symmetric parts and choosing a basis such that the symmetric part is diagonal, this can be written as
\begin{equation}
\mc{H}_{1,2} =
J \,\mb{S}_1\!\cdot\mb{S}_2
+ \mb{D}\cdot\left(\mb{S}_1\!\times\mb{S}_2\right)
+ \!\! \sum_{\alpha=x,y,z}\!\!\! G^\alpha\, {S}_1^\alpha {S}_2^\alpha.
\end{equation}
Here, $J$ is the isotropic (Heisenberg) exchange. The so-called Dzyaloshinskii-Moriya vector $\mb{D}$ quantifies the anti-symmetric exchange and $G^\alpha$ ($\alpha=x,y,z$) determine the symmetric (Ising) exchange anisotropy where $\sum_\alpha G^\alpha=0$.

Including spin-orbit coupling (SOC) into Anderson's theory of super-exchange \cite{Anderson1959} indeed all these contributions are obtained \cite{Moriya1960,Moriya1960PRL}. However, in this setting, symmetric and anti-symmetric interactions are not independent and for a single bond, to a good approximation the super-exchange interaction Hamiltonian reads \cite{Kaplan1983,SEWA1992,SEWA1993}
\begin{equation}  \label{eq:SingleBondExchangeAniso}
\mc{H}_{1,2} =
\left(\! J \!-\frac{|\mb{D}|^2}{4J} \right) \! \mb{S}_1\!\cdot \mb{S}_2
+ \mb{D}\cdot(\mb{S}_1\!\times \mb{S}_2) + \frac{1}{2J} (\mb{D}\cdot \mb{S}_1 )(\mb{D}\cdot \mb{S}_2 ).
\end{equation}
Here a \emph{single} vector $\mb{D}$ defines both the symmetric and antisymmetric exchange.
This expression, although not obvious, is fully invariant under spin rotations\cite{SEWA1992,SEWA1993}: While the term $\mb{D}\cdot(\mb{S}_1\times \mb{S}_2)$ acts to confine the spins to a plane perpendicular to the vector $\mb{D}$, the Ising contribution acts to align the spins with $\mb{D}$. Thus, for a single bond, anisotropic super-exchange will never single out a particular direction.
Nonetheless, for multiple connected bonds, {\it frustration} of the $\mb{D}$ vectors may still beak spin rotation symmetry and lead to the appearance of anisotropy\cite{SEWA1992,SEWA1993}.

Besides this mechanism, considering multiple orbitals of the ligands transmitting the super-exchange interactions, in the presence of SOC may result in {\it additional} sources of exchange anisotropy, beyond Eq.~\ref{eq:SingleBondExchangeAniso}. However, these additional effects have been argued to be much smaller in magnitude\cite{KoshibaeCommentSEWA,ReplySEWA}.

This leaves us with the following situation: When considering super-exchange interactions, to first approximation symmetric and antisymmetric anisotropy contributions always come together pointing in the same direction and at a fixed ratio. In this case, if the center of a bond corresponds to a center of inversion symmetry, not only is antisymmetric exchange anisotropy forbidden by symmetry as is commonly known\cite{DZYALOSHINSKY1958}, but we also expect symmetric exchange anisotropy to vanish. Only additional mechanisms may then lead to Ising-type anisotropy which is in principal allowed on a bond with inversion symmetry.

\subsection{Exchange anisotropy in BPCB}

For the compound BPCB with centers of inversion symmetry on the ladder rungs, we consider the following cases compatible with the crystallographic symmetry:
\begin{itemize}\vspace{-0.3\topsep}
	\setlength{\parskip}{0.5\topsep} \setlength{\itemsep}{0pt plus 1pt}
	\item [a)] Ising anisotropy on the ladder rungs

	\item [b)] Ising anisotropy on the ladder legs \setlength{\parskip}{0pt}
	\item [c)] DM anisotropy on the ladder legs
	\item [d)] Both DM and Ising anisotropy on the ladder legs as predicted for simple super-exchange (Eq.~\ref{eq:SingleBondExchangeAniso}).
\end{itemize}\vspace{-0.3\topsep}
Using the results of Sec.~\ref{sec:TripletDispersion}, for these four cases, the triplet dispersion is plotted in Fig.~\ref{fig:DispAniso}. For these plots, we have used the Heisenberg exchange constants $J_\bot/k_\mr{B}=12.7$~K and $J_\|/k_\mr{B}=3.54$~K approximately describing BPCB\cite{Blosser2018}. In all cases the anisotropy parameter was chosen such that the triplet splitting at the band minimum amounts to $50$~$\mu$eV.

Comparing the dispersion calculated for Ising-type exchange anisotropy on the ladder rungs (Fig.~\ref{fig:DispAniso}a) to the data of Fig.~\ref{fig:spectra}, we can clearly rule out this case as the dominant source of anisotropy in BPCB. The data qualitatively disagree with the calculated dispersion.
For the case of anisotropic leg exchange, all three types of anisotropy considered give qualitatively the same dispersion (Fig.~\ref{fig:DispAniso}b,c,d). There is a doubly degenerate band with a smaller bandwidth and a non-degenerate band with larger bandwidth. In all cases the calculated anisotropy splitting is slightly larger at the band minimum than at the band maximum, just as in our data%
\footnote{In Fig.~\ref{fig:DispAniso} the anisotropy parameter was chosen so as to ensure $\Delta_\mr{min}=50$~$\mu$eV. For the different models of leg anisotropy (Ising only, DM only, both), this gives $\Delta_\mr{max}=37,32$~and~$35$~$\mu$eV. However, since from 2\textsuperscript{nd} to 3\textsuperscript{rd} order perturbation expansion the dispersion still slightly changed, these small differences are probably not reliable enough to use as a criterion to select one of these three cases as the most suitable model describing the compound BPCB.}.

In Fig.~\ref{fig:NeutronFit} we show the triplet dispersion extracted from our data. These points were obtained from Gaussian fits to constant-$q_\|$-cuts, similar to the ones shown in Fig.~\ref{fig:cuts}.
The vertical bars denote the width of the observed peaks. Fitting the calculated dispersions with anisotropic leg exchange (cases b,c,d) to this data, we find excellent agreement for all three cases.  Indeed, the three calculated dispersions are so similar, that from our data, it is impossible to determine which type of exchange anisotropy on the ladder legs actually causes the observed band splitting in BPCB.
As an example, in Fig.~\ref{fig:NeutronFit} the solid line shows the dispersion calculated for anisotropic exchange on the ladder legs given by Eq.~\ref{eq:SingleBondExchangeAniso} (case d)
with fitted parameters
\begin{align}\label{eq:KSEA_FitParameters}
J_\bot/k_B=&12.77(1)\,\mr{K}, \quad \\
J_\|/k_B=3.55(1)\,\mr{K}, &\quad
D_\|/k_B=1.44(2)\,\mr{K}.\notag
\end{align}
In the parametrization of the Hamiltonian of Sec.~\ref{sec:LadderHamiltonian_LegAniso}, the coefficients are
$\tilde{J}_\|=J_\| - {D_\|^2}/{(4 J_\|)}$,
$\lambda = \tilde{J}_\|/{J_\bot}$,
$A = {D_\|^2}/{  (2 J_\| \tilde{J}_\|) }$ and
$B = {D_\|}/{\tilde{J}_\|}$.
Here, we stress that the numerical value of $D_\|/k_B=1.44$~K cannot serve as an estimate for the overall magnitude of the magnetic anisotropy which should be  estimated from the band splitting of 50~$\mu$eV, i.e. $\approx0.6$~K in units of Kelvin.

{
This leaves the question, how one might experimentally determine the type of leg exchange anisotropy dominant in BPCB.
The most promising route to answering this question would be to measure the triplet dispersion as a function of orientation and strength of a small applied magnetic field. In combination with calculations also including an applied magnetic field, such data might indeed allow to determine the precise type of exchange anisotropy on the ladder legs in BPCB. However, this is clearly beyond scope of the present work.
}

\begin{figure}
	\includegraphics{ExtremaAndFitKSEA.pdf}
	\caption{\label{fig:NeutronFit}
		Position of the triplet bands extracted from the neutron scattering data of Fig.~\ref{fig:spectra}. The dispersion drawn as a solid line represents a model with anisotropic interactions on the ladder legs (case d) as described in the text.
	}
\end{figure}

\subsection{Final remarks}
In BPCB the interactions on the ladder rungs are the strongest ones by far ($J_\bot\approx 3.6 \,J_\|$). On these bonds, antisymmetric exchange anisotropy is prohibited by symmetry and only Ising-type exchange anisotropy is allowed. It may thus seem plausible to expect this to be the dominant source of anisotropy.
Yet, experimentally we find that the anisotropy is predominantly due to the much weaker leg interactions of the ladder. Indeed, this is fully consistent with the theory of super-exchange interactions, where symmetric and antisymmetric anisotropy contributions are always coupled. If either is prohibited by symmetry, also the other will vanish. The negligible Ising anisotropy on the ladder rungs in BPCB we thus see as supporting evidence, for the theory of anisotropic super-exchange\cite{SEWA1992,SEWA1993}.

Since the theory of super-exchange predicts DM anisotropy $B\propto |\mb{D}|$ and Ising anisotropy $A\propto |\mb{D}|^2$, initially the latter was considered subdominant\cite{Moriya1960PRL,Moriya1960}. However, both contributions are now understood to be equally important\cite{Kaplan1983,SEWA1992,SEWA1993}. Whilst not immediately intuitive, the triplet dispersions calculated in Sec.~\ref{sec:LadderHamiltonian_LegAniso} illustrate this effect: In the obtained expressions $A$ appears linearly, while $B$ only appears to second order as $B^2$.

Finally, we mention an apparent contradiction: Super-exchange can never give rise to DM anisotropy exclusively. Yet, oftentimes pure DM interactions are used very successfully to explain anisotropy effects observed experimentally.
Here, the present case study offers some illustration. For the spin ladder, different types of exchange anisotropy on the ladder legs lead to almost identical triplet dispersions that all explain our data equally well. We speculate that also in other systems one might encounter similar situations and it is for this reason that pure DM anisotropy is used so successfully to explain experimental findings, even though (within the theory of super-exchange) it lacks any microscopic justification.
We are only aware of one experimental study on the helimagnet Ba$_2$CuGe$_2$O$_7$\xspace{}, where a model employing the full Hamiltonian for anisotropic super-exchange interactions (Eq.~\ref{eq:SingleBondExchangeAniso}) was compared to experimental data and indeed better agreement was found than for a model employing DM anisotropy only\cite{ZheludevPRL1998}.

Further case studies of the nature of magnetic exchange anisotropy would certainly be interesting. We also suggest, that whenever a model with pure DM interactions is considered, it would be enlightening to also consider a model with both Ising and DM anisotropy, compatible with super-exchange.

\section{Conclusion}
The prototypical spin ladder compound \BPCB{} (BPCB) has been studied by means of high resolution inelastic neutron scattering. We find a small splitting of the triplet band of $50(1)$~$\mu$eV at the band minimum and $40(2)$~$\mu$eV at the band maximum.
Further, for a spin ladder with exchange anisotropy, the triplet dispersion is calculated in a strong coupling expansion.

The Ising-type anisotropy allowed by crystallographic symmetry on the ladder rungs (by far the strongest bonds) we find to be negligible, in line with the theory of anisotropic super-exchange.
Three models with exchange anisotropy on the ladder legs, all describe the data equally well. Whilst we cannot distinguish these by comparison to our data, we note that only one is compatible with super-exchange interactions.

{\it Acknowledgement.}
This work is partially supported by the Swiss National Science Foundation under Division II. Experiments at the ISIS Neutron and Muon Source were supported by a beamtime allocation RB1720009 from the Science and Technology Facilities Council, UK.


\section{Introduction}
Interactions are ubiquitous in nature and understanding their effects is a primary challenge in physics, especially in many-body quantum systems. Generally, mean-field (MF) theories well reproduce a plethora of phenomena related to interparticle interactions. However, this approach fails when quantum fluctuations are non-negligible. The first-order correction to the MF energy, the so called Lee-Huang-Yang (LHY) term, was first calculated in 1957 \cite{LeeHuangYang} to describe bosons with hard-core repulsion. In ultracold gases the LHY term is normally negligible and only recently it has been experimentally investigated \cite{Altmeyer2007,Shin2008,Papp2008,Navon2010,Navon2011}.
Notwithstanding, there are situations where the LHY and MF contributions can be of the same order. For instance, in attractive bosonic mixtures the repulsive LHY term may stabilize the system against collapse, leading to the formation of self-bound droplets. Although very dilute, these states have a liquid-like behaviour, characterized by a core with uniform density in the large atom number limit \cite{Petrov}.

While originally predicted for Bose-Bose mixtures, quantum droplets have been first realized in dipolar gases, by exploiting the competition between contact repulsion and anisotropic dipole-dipole attraction \cite{kadau2016,FerrierBarbut2016,Schmitt2016,Ferlaino2016,FerrierBarbutJPB2016,Wenzel2017}.
Recently, tuning contact interactions through a Feshbach resonance, quantum droplets have been also observed in a spin mixture of $^{39}$K, both in the presence of an external potential \cite{Cheiney2018,Cabrera2018} and in free space \cite{Semeghini2018}.
By studying collisions between droplets it has been possible to observe the crossover between compressible and incompressible regimes \cite{Ferioli2019}. Free space $^{39}$K droplets have a short lifetime (only few ms), limited by three body losses \cite{Semeghini2018,Ferioli2019}, motivating the quest for stable long-lived droplets in different atomic mixtures. Longer lifetimes could indeed allow the investigation of many peculiar features of these states, like the characterization of the incompressible regime and the observation of self-evaporation \cite{Petrov}.

In this work, we study the formation of quantum droplets in a heteronuclear Bose-Bose mixture of $^{41}$K and $^{87}$Rb where the two species experience different trapping potentials. We exploit Feshbach resonances for tuning the interspecies interaction. In particular, we explore the transition from the weakly to the strongly attractive regime, beyond the MF threshold for collapse, where self-bound quantum droplets are predicted to exist.
We observe that in the absence of any external trapping potential, the strongly attractive mixture remains localized on a timescale of several tens of milliseconds, consistently with the formation of long-lived self-bound droplets. Furthermore, we study the size evolution and the center-of-mass dynamics of the mixture in a horizontal optical waveguide, in the presence of a species-dependent magnetic force. Also in this geometry, where both the LHY correction and the radial confinement provide a mechanism to prevent the collapse \cite{Shabat72}, we observe the formation of a localized bound state, whose center-of-mass follows a trajectory intermediate between those of two weakly attractive $^{41}$K and $^{87}$Rb clouds.
By studying the dynamics, we extract the number ratio of $^{41}$K and $^{87}$Rb atoms forming the bound state, and we find it consistent with the value predicted by the theory \cite{Petrov}. Our experimental findings are well reproduced by numerical simulations performed by using two coupled generalized Gross-Pitaevskii (GP) equations  \cite{Gross,Pitaevskii}, including the LHY correction for heteronuclear mixtures \cite{Petrov,Ancilotto2018}.
\section{Experiment}
We produce a binary condensate of $^{41}$K and $^{87}$Rb both prepared in the $\left|F=2,m_{F}=2\right\rangle$ state in a crossed optical dipole trap \cite{Burchianti2018}.
To tune the interspecies interaction, two microwave pulses transfer both species in the $\left|F=1,m_{F}=1\right\rangle$ state, and a vertical homogeneous magnetic field is increased to $\simeq 71$~G, lying between two Feshbach resonances \cite{Thalhammer2008}. Due to their different masses, $^{41}$K and $^{87}$Rb atoms have a differential vertical gravitational sag and their equilibrium positions are separated by $\sim 15~\mu$m in the optical trap.
We compensate for this offset by compressing the dipole trap and by adding a vertical magnetic field gradient $b_z=-16$~G/cm \cite{Hansen2013,Lous2017}, exploiting the different magnetic moments of the two species ($\mu_\mathrm{K} = -0.83 \;\mu_B$ and $\mu_\mathrm{Rb}= -0.52\;\mu_{B}$, at $70$~G), as shown in Appendix \ref{AppendixA}.
At this stage, we normally have a binary condensate with $N_\mathrm{K} \simeq 1 - 2 \times 10^4$ and $N_\mathrm{Rb} \simeq 4 - 6 \times 10^4$ atoms, confined in an approximately harmonic potential with different frequencies for the two atomic species: ($\nu_{x},\nu_y, \nu_z) \simeq (130, 90, 180)$~Hz for $^{41}$K and
$(100, 70, 130)$~Hz for $^{87}$Rb.
The homonuclear and heteronuclear scattering lengths are respectively $a_\mathrm{K}=65 a_0$ \cite{DErrico2007}, $a_\mathrm{Rb}=100.4 a_0$ \cite{Marte2002}, and $a_\mathrm{KRb} \simeq 0 a_0$ \cite{Simoni2008,Thalhammer2009}.  To enter the attractive regime, we then change the magnetic field to the desired value $B_F$ in $30$~ms, thus tuning $a_\mathrm{KRb}$, as detailed in Appendix \ref{AppendixA}.
The mixture is completely characterized in terms of the intraspecies $g_{\mathrm{K}}=4\pi \hbar^2 a_{\mathrm{K}}/m_{\mathrm{K}}$, $g_{\mathrm{Rb}}=4\pi \hbar^2 a_{\mathrm{Rb}}/m_{\mathrm{Rb}}$ and interspecies $g_{\mathrm{KRb}}=2\pi \hbar^2 a_{\mathrm{KRb}}/m_r$ coupling constants,  with $m_{\mathrm{K}}$ and $m_{\mathrm{Rb}}$ the atomic masses and $m_r$  the reduced mass. The onset of the MF collapse regime corresponds to $\delta g = g_\mathrm{KRb}+\sqrt{g_\mathrm{K} g_\mathrm{Rb}}=0$ ($a_\mathrm{KRb}=-75.4a_0)$ \cite{Riboli2002}. We explore the crossover from the weak ($a_\mathrm{KRb}<0$, $\delta g >0$) to the strong ($\delta g <0$) attractive regime.
For the measurements in free space, the trap beams and the magnetic field gradient $b_z$ are simultaneously switched off and the atoms are imaged in time-of-flight (TOF). For the expansion in the horizontal waveguide (along the $\hat{y}$ direction) one of the two beams of the dipole trap (labeled as crossed beam in Fig.~\ref{fig:TOF}a)
is slowly turned off by linearly reducing its intensity in $10$~ms, and the atoms are imaged {\it in-situ} after a given evolution time. In the waveguide the atoms experience a radial trapping potential with average frequencies $\nu_{r} \simeq 160$~Hz ($120$~Hz)  and an axial anti-trapping potential, produced by the magnetic field gradient, with $\nu_{y} \simeq -1.8$~Hz ($-0.6$~Hz) for $^{41}$K ($^{87}$Rb).
\begin{figure}[t]
\includegraphics[width = 0.46 \textwidth]{TOF.pdf}
\hspace{5mm}
\caption{\label{fig:TOF} (a) Schematics of the geometry of the crossed trap. (b) Absorption images in the $xz$ plane ($190\;\mu$m x $190\;\mu$m) after $\mathrm{TOF}= 24.5$~ms ($27$ ms) for $^{41}$K ($^{87}$Rb) condensate in the weakly-attractive (top row) and in the strongly-attractive (bottom row) regimes. $^{41}$K atoms are imaged at the target Feshbach field, while $^{87}$Rb atoms are imaged at $B=0$~G after an additional TOF of $2.5$~ms. (c) Average size $\sigma = \sqrt{\sigma_x \sigma_z}$ of $^{41}$K cloud as a function of the TOF in free space for two values of the interspecies scattering length: $a_\mathrm{KRb} = (-17.5\pm1.4)  a_0$ (open circles), and $a_\mathrm{KRb} = (-84.5\pm1.6) a_0$ (closed circles). The error bars correspond to the standard deviation of, typically, five independent measurements.}
\end{figure}
\section{Results}
\subsection{Dynamics in free space}
We first consider the dynamics in free space.
In Fig.~\ref{fig:TOF}b we show typical absorption images of the two condensates in the $xz$ plane and $\mathrm{TOF}= 24.5$~ms for $\delta g>0$ (top) and $\delta g<0$ (bottom).
In the former case both $^{41}$K and $^{87}$Rb samples freely expand, while in the latter case we observe that a fraction of both species bounds and forms a small and dense component. The different shape of $^{41}$K and $^{87}$Rb clouds for $\delta g<0$ is a consequence of our imaging procedure.
To detect both atomic species, in fact, we first take an absorption image of $^{41}$K sample at the magnetic field $B_F$, then we switch off $B_F$ and after further $2.5$~ms of TOF we take the image of $^{87}$Rb atoms.
The actual size of the bound state can be measured only in the first image, after which the state is dissociated and $^{87}$Rb cloud expands.
In Fig. \ref{fig:TOF}c we plot the average width (rms) $\sigma = \sqrt{\sigma_x \sigma_z}$ of $^{41}$K sample as a function of TOF for both interaction regimes. Whereas for $\delta g>0$ (open circles) $\sigma$ increases in time, for $\delta g<0$ (closed circles) it does not exceed our imaging resolution ($5\;\mu$m) up to $28$~ms of TOF, indicating the formation of a droplet state.
In our mixture the droplet density is expected to be one order of magnitude smaller than in $^{39}$K at the same $\delta g$, leading to a substantial reduction of the three-body losses. This is consistent with a two orders of magnitude longer lifetime. Such difference is mainly due to the larger scattering lengths of $^{41}$K and $^{87}$Rb, since the density scales as $n \sim (\delta g/g)^2 a^{-3}$ for approximately equal scattering lengths, as discussed in Appendix \ref{AppendixB}.
In the experiment, for TOF $>28$~ms, we observe that the atomic clouds start to expand. We attribute this to the variation of $a_{\mathrm{KRb}}$, since for larger TOF the atoms fall off the region where the Feshbach field is spatially homogeneous. Nevertheless, this time is longer than the maximum observation time of $^{39}$K droplets in free space \cite{Semeghini2018}.
In order to increase the observation time, one could exploit a magnetic gradient $b_z$ to levitate the droplet, with an average magnetic moment per particle:
\begin{equation}
\mu_\mathrm{KRb} = \frac{N^d_\mathrm{K} \mu_\mathrm{K}+N^d_\mathrm{Rb} \mu_\mathrm{Rb}}{N^d_\mathrm{K}+N^d_\mathrm{Rb}},
\end{equation}
where $N^d_\mathrm{K}$ and $N^d_\mathrm{Rb}$ are the number of $^{41}$K and $^{87}$Rb atoms forming the bound state.
However, we observed that, due to the species-dependent magnetic force, the droplet breaks up for $b_z\gtrsim 5$~G/cm.
Thus, we apply a lower magnetic field gradient $b_z=1.1$~G/cm, that, even if too small to compensate for gravity, reduces the variation of the Feshbach field during TOF \footnote{After $25$~ms of TOF we estimate a variation of the Feshbach field of $0.17$~G, resulting in a $a_{\mathrm{KRb}}$ variation of $2.9 a_0$.}. In addition, it introduces a species-dependent acceleration sufficient to spatially discriminate the bound and the unbound atoms.
In the images at $\delta g<0$ (Fig. \ref{fig:TOF}b bottom), we note indeed the presence of a halo surrounding the localized clouds of both $^{41}$K and $^{87}$Rb. Due to the finite temperature, a fraction of atoms remains unbound. Moreover, for an efficient sympathetic cooling of the mixture, the atom number of $^{87}$Rb is typically three or four times larger than that of $^{41}$K, leading to an excess of unbound $^{87}$Rb atoms. For $^{41}$K the unbound halo is upshifted  ($|\mu_\mathrm{K}|>|\mu_\mathrm{KRb}|$), whereas for $^{87}$Rb is downshifted ($|\mu_\mathrm{Rb}|<|\mu_\mathrm{KRb}|$).
By measuring the accelerations of the different components one could extract $\mu_{\mathrm{KRb}}$ and therefore the ratio of the atom numbers $N^d_\mathrm{K}/N^d_\mathrm{Rb}$ in the droplet.
However, in free space, the acceleration is dominated by gravity (for the droplet the ratio between magnetic and gravitational forces is only $0.035$).
\subsection{Waveguide dynamics}
To minimize the effect of gravity, we hold the atoms in a horizontal waveguide and let them evolve under the dominant effect of a magnetic field gradient along the $\hat{y}$ direction.
In Fig. \ref{fig:ImmaginiEspansioneGuida} we show {\it in-situ} images of the two condensates as a function of the evolution time $t$ in the waveguide. In this configuration we are able to detect the bound states for a longer time, since the center-of-mass displacement is less than $300\;\mu$m in $50$~ms.
For $\delta g>0$ (top), both species expand as $t$ increases, while for $\delta g<0$ (bottom) we observe the formation of self-bound states. Also here, as well as in free space, a fraction of the atoms remains unbound and follows the same dynamics of the weakly-attractive gas.
\begin{figure}[t]
\includegraphics[width = 0.46\textwidth]{EspansioneGuida.pdf}
\caption{\label{fig:ImmaginiEspansioneGuida} \textit{In-situ} absorption images in the $yz$ plane of the atoms expanding in a horizontal waveguide ($760\;\mu$m x $266\;\mu$m) at increasing (from top to bottom) evolution time $t$ in the weakly attractive regime $a_\mathrm{KRb} = (-11.4\pm1.4) a_0$ (top row) and in the strongly attractive one $a_\mathrm{KRb} = (-82.4\pm1.6) a_0$ (bottom row). Both $^{41}$K and $^{87}$Rb atoms are imaged as described in Fig. \ref{fig:TOF}.}
\end{figure}
In Fig. \ref{fig:GraficiEspansioneGuida} we plot the size $\sigma_y$ and the center-of-mass position $y_\mathrm{cm}$ along the longitudinal direction. For $\delta g<0$ we fit the atomic density distribution with two Gaussians to distinguish the bound and unbound components.
For clarity, $\sigma_y$ and $y_\mathrm{cm}$ are shown only for the bound components, since the unbound fraction behaviour is very similar to the one observed for $\delta g>0$ (see Appendix \ref{AppendixC}).
The measured size for $\delta g<0$ (closed symbols) is constant and is comparable with the imaging resolution in the waveguide ($10\;\mu$m) \footnote{We use two different imaging systems: for measurements in free space (in the waveguide) the imaging beam propagates along the waveguide (crossed beam) and the optical resolution is $5\;\mu$m ($10\;\mu$m).}, while for $\delta g >0$ both species expand (open symbols).
In the waveguide, the center-of-mass moves under the combined effect of a magnetic force and a residual component of the gravity, since the guide is tilted of $\theta\simeq 0.1-0.2^{\circ}$ with respect to the horizontal plane.
Assuming a constant acceleration (see Appendix \ref{AppendixC}), we extract the ratio of atom number of the two species in the bound state as
\begin{equation}
\label{eq:ratio}
\frac{N^d_\mathrm{K}}{N^d_\mathrm{Rb}}=\frac{m_\mathrm{Rb}(\alpha_d-\alpha_\mathrm{Rb})}{m_\mathrm{K}(\alpha_\mathrm{K}-\alpha_d)},
\end{equation}
where $\alpha_\mathrm{K}$ and $\alpha_\mathrm{Rb}$ are the accelerations of the unbound $^{41}$K and $^{87}$Rb clouds, and $\alpha_d$ is the acceleration of the bound state.
From $\alpha_{\mathrm{K}}$ and $\alpha_{\mathrm{Rb}}$ measured at $\delta g >0$ we find $N^d_\mathrm{K}/N^d_\mathrm{Rb} = 0.8(3)$, which is consistent with the predicted value of $\sqrt{g_\mathrm{Rb}/g_\mathrm{K}}=0.85$ \cite{Petrov}. We note that the center-of-mass motion allows a measurement of the atom number ratio which is not affected by the usual systematic uncertainty of $N$ measurements.
\begin{figure}[t]
\includegraphics[width = 0.42\textwidth]{GraficiEspansioneGuida.pdf}
\hspace{5mm}
\caption{\label{fig:GraficiEspansioneGuida} Evolution of the longitudinal size $\sigma_y$ (a) and the center-of-mass position $y_\mathrm{cm}$ along the waveguide (b) of $^{41}$K (blue open/closed circles) and $^{87}$Rb (red open/closed triangles) clouds, for $a_\mathrm{KRb} = (-11.4\pm1.4)/(-82.5\pm1.6)$~$a_0$.
Solid lines in (b) are fits with constant acceleration. The measured accelerations of the expanding clouds are $\alpha_{\mathrm{K}}=(0.276\pm0.025)$~m$/$s$^2$ and $\alpha_{\mathrm{Rb}}=(0.104\pm0.010)$~m$/$s$^2$ for $a_\mathrm{KRb} = (-11.4\pm1.4) a_0$. While for the bound components we find $\alpha_d=(0.152\pm0.017)$~m$/$s$^2$ ($\alpha_d=(0.151\pm0.022)$~m$/$s$^2$) from images of $^{41}$K ($^{87}$Rb).  The error bars are the standard deviation of, typically, five independent measurements.}
\end{figure}
\subsection{Comparison with theory}
\label{subsec:Comparison}
Both in free space and in the waveguide, we study the full range of attractive interactions.
In Fig.~\ref{fig:Sigma_vs_a} we plot the width $\sigma_x$ of $^{41}$K in free space after TOF $=26.5$~ms (Fig.~\ref{fig:Sigma_vs_a}a) and $\sigma_y$ after $t = 45$~ms of expansion in the waveguide (Fig.~\ref{fig:Sigma_vs_a}b). In both cases we observe a transition from an expanding gas for $a_\mathrm{KRb} \gtrsim - 82 a_0$ to a localized self-bound state for strong attractive interactions ($a_\mathrm{KRb} \lesssim - 82 a_0$), whose width is below the optical imaging resolution. In Fig.~\ref{fig:Sigma_vs_a}, we also show the MF threshold for collapse (dash-dot gray line) and the critical scattering length to enter the droplet regime with our experimental atom number (dashed red line) \cite{Petrov}. Both lines are calculated for a homogeneous system. In the presence of the radial confinement a stable bound state may exist also in-between the two vertical lines. This is due to the combined effect of the LHY term and the dispersion along the guide giving rise to the formation of ``solitonic'' solutions \cite{Shabat72,Cappellaro2018}.

In order to compare the experimental results with theoretical predictions \cite{Petrov,Ancilotto2018}, we have performed numerical simulations by using two coupled generalized time-dependent GP equations including the LHY term, taking into account the actual experimental preparation (see Appendix \ref{AppendixB} for details). In Fig. \ref{fig:Sigma_vs_a} the colored areas correspond to the numerical predictions including the experimental uncertainty on the atom number. Considering the finite imaging resolution, we find that the theory well reproduces the observed behavior.
\begin{figure}[t]
\subfigure{\includegraphics[width = 0.42 \textwidth]{Sigma_vs_a.pdf}}
\hspace{5mm}
\caption{\label{fig:Sigma_vs_a} (a) Size $\sigma_x$ of the $^{41}$K cloud expanding in free space as a function of the interspecies scattering length $a_\mathrm{KRb}$ after TOF $=26.5$~ms from the crossed trap. (b) Axial size $\sigma_y$ of $^{41}$K cloud as a function of the interspecies scattering length $a_\mathrm{KRb}$ after expansion in the waveguide of $t=45$~ms. The dashed red vertical line corresponds to the critical $a_\mathrm{KRb}$ for the existence of a droplet in the homogeneous case with $N_\mathrm{K}+N_\mathrm{Rb}=33000$, while the dash-dot gray vertical line corresponds to $\delta g = 0$. The solid lines are the numerical predictions for $N_\mathrm{K}+N_\mathrm{Rb}=33000$, where we fix $N_\mathrm{Rb}=\sqrt{g_\mathrm{K}/g_\mathrm{Rb}}\times N_\mathrm{K}$. The colored areas indicate the systematic relative uncertainty on the experimental atom number of 30 \%. We do not appreciate a substantial change of the average total atom number in the explored range of $a_\mathrm{KRb}$. The horizontal dotted lines show the imaging resolution of $5\;\mu$m and $10\;\mu$m for the measurements in free space and in the waveguide, respectively. Vertical error bars: standard deviation of, typically, five independent measurements.}
\end{figure}
From the simulations in free space (Fig.~\ref{fig:Sigma_vs_a}a) we find that, in the strongly attractive regime, almost spherical droplets form \cite{Petrov}. As expected, the critical scattering length to enter the droplet regime depends on the atom number. We find that this threshold  actually occurs at larger values of $|a_\mathrm{KRb}|$ with respect to the homogeneous case (see dashed red line in Fig.~\ref{fig:Sigma_vs_a}). This is due to the additional kinetic energy cost associated to the droplet surface, whose effects are more important for small droplets, as those formed in the experiment.

As for the waveguide dynamics, we observe expanding clouds in-between the two vertical lines instead of the predicted ``solitonic'' bound states. We attribute this to a twofold effect: \textit{(i)} we adjust the initial size of the mixture to match the droplet size (smaller than the one of the solitonic ground state); \textit{(ii)} the dynamics is triggered by a non adiabatic process. However, droplets formed for very large attractive scattering lengths ($|a_\mathrm{KRb}| \gtrsim 95 a_0$), as already observed in previous experiments \cite{Cheiney2018}, eventually decay into ``solitonic'' solutions due to the decrease of the atom number induced by three body losses. In the waveguide, indeed, atomic losses do not affect the existence of self-bound states allowing their observation for a wider range of $a_\mathrm{KRb}$. In free space, instead, as confirmed by numerical simulations with three body losses (Appendix \ref{AppendixB}), we are limited to $|a_\mathrm{KRb}|\lesssim 95 a_0$. This corresponds to a broader range of $\delta g$ with respect to the case of $^{39}$K droplets, pushing the investigation far away from the MF collapse.
\section{Conclusions}
In conclusion, we have observed quantum droplets in a mixture of two atomic species experiencing different potentials, providing a clear picture of their formation and evolution both in free space and under radial confinement. This proves that the presence of identical traps for the two components is not essential to the droplets production. We have found  evidence that the lifetime of our droplets is substantially longer than those in $^{39}$K, as expected. The droplet density obtained by our numerical simulations results in a lifetime up to $ 400$~ms.
This opens new perspectives for future studies of the droplets peculiar properties \cite{Petrov,Ancilotto2018} and of their collective modes \cite{Cappellaro2018}. Indeed, the lower density of our sample allows for a larger ratio $\tau_\mathrm{life}/\tau \sim n^{-1}$, with $\tau_\mathrm{life}$ and $\tau$ being respectively the lifetime and the characteristic time scale of the droplet dynamics (Appendix \ref{AppendixB}). Other interesting directions could be the creation of droplets in reduced dimensionality \cite{Petrov2016,Santos2017,Astrakharchik2018,Zin2018,Parisi2019}, where quantum fluctuations become easily dominant, and the investigation of beyond LHY corrections \cite{Cikojevic2019} and other stabilization mechanisms \cite{Adhikari2017,Gautam_2019}. Finally, the nucleation of vortex-antivortex ring pairs \cite{Kartashov2018,Ancilotto_He,Escartin2019} or more exotic scenarios, like droplets clusters \cite{Kartashov2019}, could be studied with rotating droplets.

\begin{acknowledgments}
We acknowledge inspiring discussions with M. Fattori. We gratefully thank M. Inguscio for support and A. Simoni for collisional calculations.
This work was supported by the European Commission through FET Flagship on Quantum Technologies - Qombs Project (Grant No. 820419) and by Fondazione Cassa di Risparmio Firenze through project "SUPERACI-Superfluid Atomic Circuits".
M. M. acknowledges support by the Spanish Ministry of Science, Innovation and Universities and the European Regional Development Fund FEDER through Grant No. PGC2018-101355-B-I00 (MCIU/AEI/FEDER,UE), and the Basque Government through Grant No. IT986-16.
\end{acknowledgments}