# Superfluid Fraction in an Interacting Spatially Modulated Bose-Einstein
Condensate

G. Chauveau, C. Maury, F. Rabec, C. Heintze, G. Brochier, S. Nascimbene, J.
Dalibard, J. Beugnon<EMAIL_ADDRESS>Laboratoire Kastler Brossel, Collège
de France, CNRS, ENS-PSL University, Sorbonne Université, 11 Place Marcelin
Berthelot, 75005 Paris, France  S.M. Roccuzzo and S. Stringari
<EMAIL_ADDRESS>Pitaevskii BEC Center, CNR-INO and Dipartimento di
Fisica, Università di Trento, I-38123 Trento, Italy  Trento Institute for
Fundamental Physics and Applications, INFN, Italy

(August 28, 2024)

###### Abstract

At zero temperature, a Galilean-invariant Bose fluid is expected to be fully
superfluid. Here we investigate theoretically and experimentally the quenching
of the superfluid density of a dilute Bose-Einstein condensate due to the
breaking of translational (and thus Galilean) invariance by an external 1D
periodic potential. Both Leggett’s bound fixed by the knowledge of the total
density and the anisotropy of the sound velocity provide a consistent
determination of the superfluid fraction. The use of a large-period lattice
emphasizes the important role of two-body interactions on superfluidity.

Superfluidity is a unique state of matter exhibited by quantum many-body
systems in special conditions of temperature and interactions. It is
characterized by the absence of viscosity and by many other peculiar
phenomena, like the occurrence of quantized vortices, the reduction of the
moment of inertia, the propagation of second sound at finite temperature, and
Josephson effects. Superfluidity was first discovered in liquid helium Kapitza
(1938); Allen and Misener (1938). More recently, an impressive amount of
scientific activity has concerned the superfluid behavior of ultracold atomic
gases (for a review see Refs. Pitaevskii and Stringari (2016); Pethick and
Smith (2008); Leggett (2006)).

A key quantity characterizing superfluidity is the fraction of the total
density, the so-called superfluid fraction, which determines superfluid
transport phenomena. According to Landau’s theory of superfluidity Landau
(1941), at nonzero temperature the superfluid density does not coincide with
the total density. The thermal occupation of elementary excitations provides
the normal (non-superfluid) component responsible, for example, for the non-
vanishing moment of inertia and the propagation of second sound Pitaevskii and
Stringari (2016). The measurement of second sound velocity provides unique
information on the temperature dependence of the superfluid density, in both
liquid helium Peshkov (1944) and quantum gases Sidorenkov  _et al._ (2013);
Christodoulou  _et al._ (2021).

At zero temperature, superfluid and total densities still do not always
coincide as illustrated by the celebrated superfluid to Mott-insulator
transition Fisher  _et al._ (1989); Jaksch  _et al._ (1998); Greiner  _et al._
(2002). Even in the mean-field regime relevant for the present work,
consistent with the applicability of Gross-Pitaevskii theory, quenching of the
superfluid density can occur when translation or Galilean invariances are
broken, resulting in important consequences on the excitation spectrum. Such
effects have been already pointed out theoretically in the presence of
disorder Giorgini  _et al._ (1994); Gaul  _et al._ (2009), external periodic
potentials Javanainen (1999); Krämer  _et al._ (2002); Smerzi  _et al._
(2002); Machholm  _et al._ (2003); Taylor and Zaremba (2003); Watanabe  _et
al._ (2008), supersolidity Saccani  _et al._ (2012); Macrì  _et al._ (2013)
and spin-orbit coupling Li  _et al._ (2013); Martone  _et al._ (2012); Zhang
_et al._ (2016). Experimentally, the effects of the quenched superfluidity on
the collective frequencies of a harmonically trapped gas were investigated in
Ref. Cataliotti  _et al._ (2001). A similar situation emerges in astrophysics
in the context of neutron stars where the periodic lattice of nuclei
influences the superfluid density in the inner crust Chamel (2012); Watanabe
and Pethick (2017).

Here we provide a combined theoretical and experimental investigation of the
reduction of the superfluid fraction caused by the presence of a periodic
potential in a weakly interacting Bose-Einstein condensate (BEC) confined in a
box. We determine the superfluid density employing Leggett’s result Leggett
(1970, 1998), which is based on the knowledge of the in situ modulated total
density profile $\rho(\boldsymbol{r})$, experimentally available thanks to the
use of a large-period lattice. We also present an independent measurement of
the superfluid fraction taht exploits the anisotropic character of the sound
velocity.

#### Superfluid fraction in a modulated potential.

We consider a two-dimensional (2D) weakly interacting BEC confined in a box of
size $L\times L$, in the presence of the one-dimensional (1D) spatially
periodic potential $V(x)=V_{0}\cos\left(qx\right)$. This potential brings two
energy scales to the problem, its amplitude $V_{0}$ and the “recoil energy”
$\epsilon_{q}=\hbar^{2}q^{2}/2m$, where $m$ is the mass of an atom. Atomic
interactions provide the third energy scale relevant for the problem. They are
conveniently characterized by the chemical potential $\mu_{0}=g\rho_{0}$ of a
uniform condensate with a density equal to the average value
$\rho_{0}=\langle\rho(x)\rangle$, where $\rho(x)$ is the density profile and
the average is calculated over one period of the potential $V(x)$. The
interaction coupling constant between atoms $g=4\pi\hbar^{2}a_{s}/m$ is fixed
by the $s$-wave scattering length $a_{s}$.

In such a configuration, the superfluid fraction is an anisotropic rank-2
tensor with eigenaxes $x,y$ and with diagonal elements denoted $f_{s,\alpha}$
in the following ($\alpha=x,y$). They can be calculated by applying the
perturbation $-v\hat{P}_{\alpha}$ to the system, where
$\hat{P}_{\alpha}=\sum_{j=1}^{N}\hat{p}_{j,\alpha}$ is the momentum operator
along the axis $\alpha$ and $N$ the number of particles. This corresponds to
working in the frame moving with velocity $v$ with respect to the laboratory
frame. Only the normal part reacts to the perturbation, so that, by
calculating the average momentum $\langle\hat{P}_{\alpha}\rangle$ and imposing
periodic boundary conditions, one accesses to the superfluid fraction along
the axis $\alpha$

$f_{s,\alpha}=1-\lim_{v\to 0}\frac{\langle\hat{P}_{\alpha}\rangle}{Nmv}\;.$
(1)

A similar procedure, applied to the case of a rotating configuration,
employing the angular momentum operator rather than the linear momentum
operator, gives access to the moment of inertia, whose deviation from the
rigid value provides direct evidence of superfluid effects Stringari (1996).

In the presence of the periodic potential $V(x)$, the motion of the fluid is
slowed down along the $x$ direction, reflecting the quenching of the
superfluid density along this direction: $\rho_{s,x}<\rho_{0}$. The superfluid
density evaluated along the transverse $y$ direction is instead not modified:
$\rho_{s,y}=\rho_{0}$. The Gross-Pitaevskii equation (GPE) describing the
weakly interacting BEC, solved in the frame moving with velocity $v$, i.e.
subject to the constraint $-vP_{x}$, yields, according to the definition (1),
the result REF

$f_{s,x}=\frac{\rho_{s,x}}{\rho_{0}}=\frac{1}{\langle\rho(x)\rangle\
\langle\frac{1}{\rho(x)}\rangle}.$ (2)

According to the seminal work by Leggett Leggett (1970, 1998), the right-hand
side of (2) provides generally an upper bound to the superfluid density.
Remarkably, the bound reduces to an identity in the case of a weakly
interacting BEC subject to a 1D periodic potential 111Result (2) for the
superfluid density does not hold only in the presence of periodic potentials
of the form $V_{0}\cos\left(qx\right)$. For example, it provides the Josephson
energy in the presence of a junction separating two BECs Zapata  _et al._
(1998).

#### Effective mass and sound propagation.

Result (2) may be surprising because it relates a transport property (the
superfluid density) to a static quantity (the equilibrium density profile).
The concept of effective mass, commonly used in the context of interacting
Bose Krämer  _et al._ (2002); Zwerger (2003); Taylor and Zaremba (2003) and
Fermi gases Watanabe  _et al._ (2008) placed in a periodic potential,
elucidates this relation. In the present case, the superfluid fraction of the
BEC, defined according to (2), exactly coincides with the ratio

$\frac{m}{m_{x}^{*}}=f_{s,x}\;,$ (3)

where the effective mass $m_{x}^{*}$ fixes the curvature of the energy band
along the $x$ direction for small values of the quasi-momentum REF .

The relation (3) between the effective mass and $f_{s,\alpha}$ illustrates the
crucial role of the superfluid density in the propagation of sound. The
hydrodynamic formalism of superfluids indeed provides the following expression
for the velocity of a sound wave propagating along the $x$ direction in the
presence of $V(x)$ Krämer  _et al._ (2002); Pitaevskii and Stringari (2016);
Bensoussan  _et al._ (2011); Cataliotti  _et al._ (2001)

$c^{2}_{x}=\frac{1}{m_{x}^{*}\kappa}={f_{s,x}}\;\frac{1}{m\kappa}\;,$ (4)

where $\kappa=\left[\rho_{0}\,\partial_{\rho_{0}}\mu(\rho_{0})\right]^{-1}$ is
the compressibility of the gas.

Figure 1: Solid line: Superfluid fraction calculated for a 87Rb condensate by
injecting the numerical solution of the GPE inside Leggett’s formula (2) for
$\rho_{0}=60\,\rm$\mathrm{\SIUnitSymbolMicro m}$^{-2}$, $mg/\hbar^{2}=0.15$
and $2\pi/q=3.93\,\rm$\mathrm{\SIUnitSymbolMicro m}$$, corresponding to
$\mu_{0}/\epsilon_{q}=7.0$. Dotted line: LDA result. Red circles: Experimental
results obtained using Leggett’s formula. Violet squares: Experimental results
obtained from speeds of sound. The dashed line shows the prediction in the
opposite regime, $\mu_{0}\ll\epsilon_{q}$, of a very weakly-interacting
system. In all figures, the error bar represents the statistical uncertainties
of the measurements.

The value of the sound velocity propagating along the transverse $y$ direction
is

$c^{2}_{y}=\frac{1}{m\kappa}\;,$ (5)

the effective mass $m_{y}^{*}$ being equal to the bare mass in this case. The
ratio between Eqs. 4 and 5 then provides

$f_{s,x}=\frac{c^{2}_{x}}{c^{2}_{y}}.$ (6)

The superfluid fraction can thus be determined either through the explicit
knowledge of the equilibrium density profile via (2) or through the
measurement of the ratio (6).

#### Limiting cases.

Results (2) and (6) hold for any values of the dimensionless parameters
$V_{0}/\mu_{0}$ and $\epsilon_{q}/\mu_{0}$, as long as the description of the
$T=0$ Bose gas by a macroscopic wave function is valid, i.e., as long as
quantum phase fluctuations between neighboring sites (a precursor of the
superfluid to Mott-insulator transition) can be ignored. We now examine some
limiting cases where $\rho(x)$ and $f_{s,x}$ take a simple expression.

We start with the very weakly-interacting regime where
$\mu_{0}\ll\epsilon_{q}$. In this case, the GPE approaches the Schrödinger
equation for a single particle subject to the periodic potential $V(x)$. In
this regime, the identity $m^{*}/m=\langle\rho\rangle\,\langle 1/\rho\rangle$
was already noticed in Ref. Leggett .

The opposite case $\epsilon_{q}\ll\mu_{0}$ is described by the Local Density
Approximation (LDA). The validity condition of the LDA is equivalent to
imposing that the period of the potential $2\pi/q$ be much larger than the
healing length $\hbar/\sqrt{2mg\rho_{0}}$. In the LDA, the equilibrium density
becomes $\rho^{\rm(LDA)}(x)=\rho_{0}-\rho_{1}\cos(qx)$ with
$\rho_{0}=\mu_{0}/g$ and $\rho_{1}=V_{0}/g$. When injected into (2), this
gives

$f^{\rm(LDA)}_{s,x}=\left(1-\frac{\rho_{1}^{2}}{\rho_{0}^{2}}\right)^{1/2}=\left(1-\frac{V_{0}^{2}}{\mu_{0}^{2}}\right)^{1/2}.$
(7)

For values $V_{0}>\mu_{0}$, the LDA density vanishes in a finite region within
each period of $V(x)$, with the consequent vanishing of the superfluid
fraction, (dotted line in Fig. 1).

We now consider the third limiting case of small $V_{0}$, which can be
addressed using the formalism of the static density response function. An
expansion of the solution of the GPE in powers of $V_{0}$ yields the amplitude
of the first Fourier component of $\rho(x)$

$\frac{\rho_{1}}{\rho_{0}}=\frac{2V_{0}}{2\mu_{0}+\epsilon_{q}}+{\cal
O}(V_{0}^{3}).$ (8)

Using (2), one finds the superfluid fraction Taylor and Zaremba (2003)

$f_{s,x}=1-\frac{2V^{2}_{0}}{(2\mu_{0}+\epsilon_{q})^{2}}+{\cal
O}(V_{0}^{4})\;,$ (9)

confirming the LDA result given above when we take the limit
$\epsilon_{q}/\mu_{0}\to 0$. Note that (8,9) also hold in the opposite limit
of large $\epsilon_{q}/\mu_{0}$, where the superfluid density no longer
depends on the interaction.

The expansion of the solution of the GPE in powers of $V_{0}$ also provides
the compressibility

$\kappa=\mu_{0}^{-1}\left[1-\frac{2V_{0}^{2}\epsilon_{q}}{(2\mu_{0}+\epsilon_{q})^{3}}\right]+{\cal
O}(V_{0}^{4})\;,$ (10)

showing that, different from the expansion (9) for the superfluid density, the
$V_{0}^{2}$ correction to the compressibility vanishes in the LDA limit
$\epsilon_{q}/\mu_{0}\to 0$. In this limit, we thus predict from 4 and 5 that,
at order 2 in $V_{0}$, the speed of sound $c_{y}$ in the direction
perpendicular to the lattice will not be affected by the presence of the
lattice, whereas $c_{x}$ will be reduced by an amount directly related to
$f_{s,x}$.

#### The ideal gas limit.

The addition of a lattice on a Bose gas sheds interesting light on the
controversial question of the possible superfluidity in the ideal case. The
fact that the Landau criterion is not satisfied points to a non-superfluid
character of the ideal gas, while the approach based on twisted boundary
conditions leads to $f_{s}=1$ for this system. To remove this ambiguity, we
take a gas with chemical potential $\mu_{0}$ placed in a lattice of large
spatial period ($\epsilon_{q}\ll V_{0}$) and consider the two limits (i)
$V_{0}\to 0$ and (ii) $\mu_{0}\to 0$. The order in which these limits are
taken is crucial. If we take limit (i) first (i.e. $\epsilon_{q}\ll
V_{0}\ll\mu_{0}$) and then limit (ii), we find $f_{s}=1$, see (9). Conversely,
taking first the limit (ii) (i.e. $\epsilon_{q},\mu_{0}\ll V_{0}$) leads to
$f_{s}\approx 0$ (see dashed line in Fig. 1). In our opinion, the latter
approach is more relevant as it implicitly takes into account the residual
(possibly disordered) modulated potentials acting on the gas.

#### Experimental setup.

We now describe the experimental determination of the superfluid fraction of a
planar BEC subjected to a sinusoidal potential along $x$. The setup has been
detailed in Refs. Ville  _et al._ (2017, 2018). We start from a single
quasi-2D Bose gas of 87Rb atoms confined in an optical dipole trap made of a
combination of repulsive laser beams at a wavelength $\lambda=532\,$nm. We
load all atoms around a single node of an optical lattice, which provides a
strong confinement along the vertical direction $z$. It leads to an
approximate harmonic confinement of frequency $\omega_{z}/2\pi\approx 3.7$
kHz. The associated characteristic length,
$\ell_{z}=\sqrt{\hbar/m\omega_{z}}\approx 180$ nm, is large compared to the
s-wave scattering length $a_{s}=5.3\,$nm, which leads to a quasi-2D regime
where collisions keep their three-dimensional (3D) character. The effective 2D
coupling constant describing the interactions in the cloud is
$mg/\hbar^{2}=\sqrt{8\pi}a_{s}/\ell_{z}\approx 0.15$. The 2D character of our
gas is not crucial for this experiment, which could also be performed in a 3D
box-like potential Gaunt  _et al._ (2013).

The in-plane confinement is created by spatially-shaped laser beams. A first
beam creates a square box potential of size
$L=40\,\rm$\mathrm{\SIUnitSymbolMicro m}$$. A second beam imposes the
sinusoidal potential modulated along the $x$ axis with a tunable amplitude
$V_{0}$ from 0 to 80 nK Zou  _et al._ (2021); REF . The lattice period
$d=3.93(4)\,\rm$\mathrm{\SIUnitSymbolMicro m}$$ and the average 2D density
$\rho_{0}=60(3)\,\rm$\mathrm{\SIUnitSymbolMicro m}$^{-2}$ are fixed. This
corresponds to $\mu_{0}/k_{B}\approx 50\,$nK and $\epsilon_{q}/k_{B}=7.1\,$nK.
The temperature of the gas is below the lowest measurable value in our setup,
i.e. $<20\,$nK.

Figure 2: (a) In situ absorption image of the 2D gas modulated with a lattice
along $x$ of period $3.93\,\rm$\mathrm{\SIUnitSymbolMicro m}$$ and amplitude
$V_{0}=54(5)$ nK. The length of the scale bar is 20
$\mathrm{\SIUnitSymbolMicro m}$. (b) Density profile integrated along $x$
(squares) and $y$ (circles) and the fit to a sinusoidal modulation (solid blue
line) for $\rho^{\rm(meas)}(x)$. The pixel size is 1.15
$\mathrm{\SIUnitSymbolMicro m}$. (c) Fourier components of the density
modulation $\rho_{n}/\rho_{0}$ versus the lattice depth $V_{0}$ for $n=1$
(circles) and $n=2$ (squares). Solid lines represent the corresponding
predictions from the GPE. The dotted line is the weak lattice limit of (8).
Figure 3: Speed of sound measurement. (a) Center of mass position of the cloud
after excitation as a function of time for an excitation along or
perpendicular to the lattice for $V_{0}/k_{B}=41(4)$ nK. The solid lines are
sinusoidal fits to the data giving $\nu_{x}=19(1)\,$Hz and
$\nu_{y}=27(1)\,$Hz. (b) Extracted speed of sound along $y$ (black squares)
and $x$ (blue circles) axes for different lattice depths. The solid lines are
the prediction from the GPE.

#### Superfluid fraction from Leggett’s formula.

To use Leggett’s result (2), we measure the in situ 2D density profile
$\rho^{\rm(meas)}(x,y)$ in the presence of the lattice using absorption
imaging, see Fig. 2(a). We integrate it along $y$ to obtain the 1D profile
$\rho^{\rm(meas)}(x)$ (Fig. 2(b)). For an ideal imaging system 222We use a
weakly-saturating imaging beam and partial transfer imaging to image low
density clouds and stay in the linear absorption regime.,
$\rho^{\rm(meas)}(x)=\rho(x)$ but finite optical resolution alters this
relation and has to be included in the analysis. The expected density
distribution can be expanded in Fourier series
$\rho_{0}-\sum_{n>0}\rho_{n}\cos(nqx)$, where the role of higher harmonics
becomes increasingly important for large $V_{0}$. We model our optical
resolution by multiplicative coefficients $\beta_{n}<1$:
$\rho^{\rm(meas)}(x)=\rho_{0}-\sum_{n>0}\beta_{n}\rho_{n}\cos(nqx)$.

We calibrate the first coefficients $\beta_{n}$ by studying the density
response to a lattice of wave number $q^{\prime}=nq$ for low lattice depths.
In this case, the density modulation $\rho(x)$ is dominated by its first
harmonic and we fit the measured profiles to a sinusoidal function whose
amplitude is adjusted to prediction (8). This adjustment provides
$\beta_{1}=0.73(2)$ and $\beta_{2}=0.27(6)$, while the values of the
coefficients $n\geq 3$ are below our experimental detectivity.

We show in Fig. 2(c) the values of $\rho_{n}=\rho_{n}^{(\rm meas)}/\beta_{n}$
for $n=1,2$. Both measurements are in good agreement with the predictions of
the GPE (solid lines) over all the explored range of values of $V_{0}$. From
this measurement and restricting $\rho(x)$ to its two first Fourier
components, we calculate Leggett’s formula (2) and we plot the result as
circles in Fig. 1. We discuss in Ref. REF  the effect of the truncation of the
Fourier series on the solution of the GPE and confirm that, in our case,
restricting to the first two harmonics already gives a good estimate of
$f_{s,x}$.

#### Superfluid fraction from speed of sound.

We determine the speeds of sound along $x$ and $y$ by studying the response of
the cloud to an external perturbation of its density. Here, the perturbation
consists adding, during the preparation of the cloud, a weak linear magnetic
potential along $x$ or $y$, of amplitude $\approx$ 0.1
nK/$\mathrm{\SIUnitSymbolMicro m}$. At time $t=0$, we abruptly switch off this
potential and measure the evolution of the center of mass of the cloud, see
Fig. 3(a). For a perturbation along $x$, we observe a smaller frequency than
the one obtained for the same excitation along the $y$ axis, a clear signature
of the modification of the superfluid transport properties due to the presence
of the lattice.

From the frequency $\nu_{x,y}$ of the fitted oscillations, we determine the
speed of sound through the relation $c_{x,y}=2L\nu_{x,y}$ valid when the
lattice period and the healing length are both much smaller than the phonon
wavelength, equal here to $2L$. We show our measurements in Fig. 3(b) as a
function of the lattice depth. For a perturbation along $y$, we observe a
small increase of $c_{y}$ with $V_{0}$, which can be attributed to the
modification of the compressibility of the modulated gas with respect to the
uniform case (see Eqs. 5 and 10). Along the axis of the lattice, we note a
strong decrease of the speed of sound with the amplitude of the modulating
potential, which we associate with the decrease of the superfluid fraction of
the cloud. In addition, we plot with solid lines the result of a simulation of
the experimental protocol with the GPE. We observe an excellent agreement for
both $c_{x}$ and $c_{y}$. The superfluid fraction $f_{s,x}$ obtained from the
ratio (6) is plotted in Fig. 1.

Discussion and conclusion. The determination of the superfluid fraction
$f_{s,x}$ based on sound propagation is in excellent agreement with the
prediction from the GPE. The determination of $f_{s,x}$ based on Leggett’s
formula, although limited by the finite resolution of our optical system, also
agrees well with the prediction.

More generally, one may favor one of the two methods depending on the system
under study. For example, in a spin-orbit coupled BEC that violates Galilean
invariance Zhang  _et al._ (2016), the sound velocity measurement will give
access to the superfluid fraction, while the density may remain uniform, in
which case Leggett’s bound is not relevant. Conversely in supersolid BECs
Léonard  _et al._ (2017); Li  _et al._ (2017); Chomaz  _et al._ (2019);
Böttcher  _et al._ (2019); Tanzi  _et al._ (2019); Norcia  _et al._ (2021);
Tanzi  _et al._ (2021), the excitation spectrum is, in general, more complex
Josserand  _et al._ (2007); Hofmann and Zwerger (2021) and the sound velocity
measurement is not directly applicable to extract $f_{s}$, which, at least in
one-dimensional-like configurations, may be instead calculated using Leggett’s
formula Roccuzzo (2021). A challenging question concerns the determination of
the superfluid density in higher dimensions if the total density profile is
not factorizable along the various directions Leggett (1998), as in the case
of 2D dipolar supersolids Norcia  _et al._ (2021) and of a vortex lattice
Sonin (1987). Our work also paves the way for the investigation of the
superfluid fraction in other density modulated quantum gases, like Fermi
superfluids and 1D and disordered systems. It could be extended to study the
links between the quenching of the superfluid fraction and the emergence of
number squeezing and phase fluctuations effects Orzel  _et al._ (2001); Smerzi
_et al._ (2002) for deep periodic potentials as well as to investigate the
consequence of finite temperature effects.

While we were completing this work, a manuscript by Tao et al. explored
anisotropic superfluidity in a periodically modulated trapped BEC gas Tao  _et
al._ (2023). We employ here lattices with much larger periods, allowing for an
explicit measurement of Leggett’s formula and enhancing the role of two-body
interactions.

###### Acknowledgements.

Acknowledgments. G.C., C.M. and F.R. contributed equally to this work. We
acknowledge the support by ERC (Grant Agreement No 863880) and ANR (Grant
Agreement ANR-18-CE30-0010). We thank Brice Bakkali-Hassani for his
participation at the early stage of the project. We also acknowledge Tony
Leggett, Anne-Laure Dalibard, Stefano Giorgini, Augusto Smerzi, Alessio
Recati, Christopher Pethick, and Ian Spielman for fruitful discussions.

## References

  * Kapitza (1938) P. Kapitza, “Viscosity of liquid helium below the $\lambda$-point,” Nature 141, 74 (1938).
  * Allen and Misener (1938) J.F. Allen and A.D. Misener, “Flow of liquid helium ii,” Nature 141, 75 (1938).
  * Pitaevskii and Stringari (2016) L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation and superfluidity_ (Oxford University Press, 2016).
  * Pethick and Smith (2008) C.J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * Leggett (2006) A.J. Leggett,  _Quantum liquids: Bose condensation and Cooper pairing in condensed-matter systems_ (Oxford university press, 2006).
  * Landau (1941) L. Landau, “Theory of the superfluidity of helium II,” Phys. Rev. 60, 356–358 (1941).
  * Peshkov (1944) V. P. Peshkov, “Second sound in He II,” J. Phys. Moscow 8, 131 (1944).
  * Sidorenkov  _et al._ (2013) L.A. Sidorenkov, M.K. Tey, R. Grimm, Y.-H. Hou, L. Pitaevskii, and S. Stringari, “Second sound and the superfluid fraction in a Fermi gas with resonant interactions,” Nature 498, 78–81 (2013).
  * Christodoulou  _et al._ (2021) P. Christodoulou, M. Gałka, N. Dogra, R. Lopes, J. Schmitt, and Z. Hadzibabic, “Observation of first and second sound in a Berezinskii-Kosterlitz-Thouless superfluid,” Nature 594, 191 (2021).
  * Fisher  _et al._ (1989) M.P.A. Fisher, Peter B. Weichman, G. Grinstein, and Daniel S. Fisher, “Boson localization and the superfluid-insulator transition,” Phys. Rev. B 40, 546 (1989).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C.W. Gardiner, and P. Zoller, “Cold bosonic atoms in optical lattices,” Phys. Rev. Lett. 81, 3108 (1998).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T.W. Hänsch, and I. Bloch, “Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms,” Nature 415, 39 (2002).
  * Giorgini  _et al._ (1994) S. Giorgini, L. Pitaevskii, and S. Stringari, “Effects of disorder in a dilute Bose gas,” Phys. Rev. B 49, 12938 (1994).
  * Gaul  _et al._ (2009) C. Gaul, N. Renner, and C.A. Müller, “Speed of sound in disordered Bose-Einstein condensates,” Phys. Rev. A 80, 053620 (2009).
  * Javanainen (1999) J. Javanainen, “Phonon approach to an array of traps containing Bose-Einstein condensates,” Phys. Rev. A 60, 4902 (1999).
  * Krämer  _et al._ (2002) M. Krämer, L. Pitaevskii, and S. Stringari, “Macroscopic dynamics of a trapped Bose-Einstein condensate in the presence of 1D and 2D optical lattices,” Phys. Rev. Lett. 88, 180404 (2002).
  * Smerzi  _et al._ (2002) A. Smerzi, A. Trombettoni, P. G. Kevrekidis, and A. R. Bishop, “Dynamical superfluid-insulator transition in a chain of weakly coupled Bose-Einstein condensates,” Phys. Rev. Lett. 89, 170402 (2002).
  * Machholm  _et al._ (2003) M. Machholm, C.J. Pethick, and H. Smith, “Band structure, elementary excitations, and stability of a Bose-Einstein condensate in a periodic potential,” Phys. Rev. A 67, 053613 (2003).
  * Taylor and Zaremba (2003) E. Taylor and E. Zaremba, “Bogoliubov sound speed in periodically modulated Bose-Einstein condensates,” Phys. Rev. A 68, 053611 (2003).
  * Watanabe  _et al._ (2008) G. Watanabe, G. Orso, F. Dalfovo, L.P. Pitaevskii, and S. Stringari, “Equation of state and effective mass of the unitary Fermi gas in a one-dimensional periodic potential,” Phys. Rev. A 78, 063619 (2008).
  * Saccani  _et al._ (2012) S. Saccani, S. Moroni, and M. Boninsegni, “Excitation spectrum of a supersolid,” Phys. Rev. Lett. 108, 175301 (2012).
  * Macrì  _et al._ (2013) T. Macrì, F. Maucher, F. Cinti, and T. Pohl, “Elementary excitations of ultracold soft-core bosons across the superfluid-supersolid phase transition,” Phys. Rev. A 87, 061602 (2013).
  * Li  _et al._ (2013) Y. Li, G.I. Martone, L.P. Pitaevskii, and S. Stringari, “Superstripes and the excitation spectrum of a spin-orbit-coupled Bose-Einstein condensate,” Phys. Rev. Lett. 110, 235302 (2013).
  * Martone  _et al._ (2012) G.I. Martone, Y. Li, L.P. Pitaevskii, and S. Stringari, “Anisotropic dynamics of a spin-orbit-coupled Bose-Einstein condensate,” Phys. Rev. A 86, 063621 (2012).
  * Zhang  _et al._ (2016) Y.-C. Zhang, Z.-Q. Yu, T.K. Ng, S. Zhang, L. Pitaevskii, and S. Stringari, “Superfluid density of a spin-orbit-coupled Bose gas,” Phys. Rev. A 94, 033635 (2016).
  * Cataliotti  _et al._ (2001) F.S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio, “Josephson junction arrays with Bose-Einstein condensates,” Science 293, 843–846 (2001).
  * Chamel (2012) N. Chamel, “Neutron conduction in the inner crust of a neutron star in the framework of the band theory of solids,” Phys. Rev. C 85, 035801 (2012).
  * Watanabe and Pethick (2017) G. Watanabe and C.J. Pethick, “Superfluid density of neutrons in the inner crust of neutron stars: New life for pulsar glitch models,” Phys. Rev. Lett. 119, 062701 (2017).
  * Leggett (1970) A.J. Leggett, “Can a solid be "superfluid"?” Phys. Rev. Lett. 25, 1543–1546 (1970).
  * Leggett (1998) A.J. Leggett, “On the superfluid fraction of an arbitrary many-body system at $T=0$,” J. Stat. Phys. 93, 927–941 (1998).
  * Stringari (1996) S. Stringari, “Moment of inertia and superfluidity of a trapped Bose gas,” Phys. Rev. Lett. 76, 1405–1408 (1996).
  * (32) “For details, see Supplemental Materials, which includes Refs. Dorrer and Zuegel (2007); Corman  _et al._ (2017),” .
  * Note (1) Result (2) for the superfluid density does not hold only in the presence of periodic potentials of the form $V_{0}\mathop{cos}\nolimits\left(qx\right)$. For example, it provides the Josephson energy in the presence of a junction separating two BECs Zapata  _et al._ (1998).
  * Zwerger (2003) W. Zwerger, “Mott–Hubbard transition of cold atoms in optical lattices,” J. Opt. B: Quantum Semiclass. Opt. 5, S9 (2003).
  * Bensoussan  _et al._ (2011) A. Bensoussan, J.-L. Lions, and G. Papanicolaou,  _Asymptotic analysis for periodic structures_, Vol. 374 (American Mathematical Soc., 2011).
  * (36) A.J. Leggett, Private communication (2022).
  * Ville  _et al._ (2017) J.L. Ville, T. Bienaimé, R. Saint-Jalm, L. Corman, M. Aidelsburger, L. Chomaz, K. Kleinlein, D. Perconte, S. Nascimbène, J. Dalibard, and J. Beugnon, “Loading and compression of a single two-dimensional Bose gas in an optical accordion,” Phys. Rev. A 95, 013632 (2017).
  * Ville  _et al._ (2018) J. L. Ville, R. Saint-Jalm, É. Le Cerf, M. Aidelsburger, S. Nascimbène, J. Dalibard, and J. Beugnon, “Sound propagation in a uniform superfluid two-dimensional Bose gas,” Phys. Rev. Lett. 121, 145301 (2018).
  * Gaunt  _et al._ (2013) A.L. Gaunt, T.F. Schmidutz, I. Gotlibovych, R.P. Smith, and Z. Hadzibabic, “Bose-Einstein condensation of atoms in a uniform potential.” Phys. Rev. Lett. 110, 200406 (2013).
  * Zou  _et al._ (2021) Y.-Q. Zou, É. Le Cerf, B. Bakkali-Hassani, C. Maury, G. Chauveau, P.C.M. Castilho, R. Saint-Jalm, S. Nascimbene, J. Dalibard, and J. Beugnon, “Optical control of the density and spin spatial profiles of a planar Bose gas,” J.Phys. B : At. Mol. Opt. 54, 08LT01 (2021).
  * Note (2) We use a weakly-saturating imaging beam and partial transfer imaging to image low density clouds and stay in the linear absorption regime.
  * Léonard  _et al._ (2017) J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, “Supersolid formation in a quantum gas breaking a continuous translational symmetry,” Nature 543, 87 (2017).
  * Li  _et al._ (2017) J.-R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F.C. Top, A.O. Jamison, and W. Ketterle, “A stripe phase with supersolid properties in spin–orbit-coupled Bose–Einstein condensates,” Nature 543, 91 (2017).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, “Long-lived and transient supersolid behaviors in dipolar quantum gases,” Phys. Rev. X 9, 021012 (2019).
  * Böttcher  _et al._ (2019) F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, “Transient supersolid properties in an array of dipolar quantum droplets,” Phys. Rev. X 9, 011051 (2019).
  * Tanzi  _et al._ (2019) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, “Observation of a dipolar quantum gas with metastable supersolid properties,” Phys. Rev. Lett. 122, 130405 (2019).
  * Norcia  _et al._ (2021) M.A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M.J. Mark, R.N. Bisset, L. Santos, and F. Ferlaino, “Two-dimensional supersolidity in a dipolar quantum gas,” Nature 596, 357 (2021).
  * Tanzi  _et al._ (2021) L. Tanzi, J.G. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G Modugno, “Evidence of superfluidity in a dipolar supersolid from nonclassical rotational inertia,” Science 371, 1162 (2021).
  * Josserand  _et al._ (2007) C. Josserand, Y. Pomeau, and S. Rica, “Patterns and supersolids,” Eur. Phys. J. Spec. Top. 146, 47 (2007).
  * Hofmann and Zwerger (2021) J. Hofmann and W. Zwerger, “Hydrodynamics of a superfluid smectic,” J. Stat. Mech. , 033104 (2021).
  * Roccuzzo (2021) S. Roccuzzo,  _Supersolidity in a dipolar quantum gas_ , Ph.D. thesis, Università di Trento (2021).
  * Sonin (1987) E.B. Sonin, “Vortex oscillations and hydrodynamics of rotating superfluids,” Rev. Mod. Phys. 59, 87 (1987).
  * Orzel  _et al._ (2001) C. Orzel, A.K. Tuchman, M.L. Fenselau, M. Yasuda, and M.A. Kasevich, “Squeezed states in a Bose-Einstein condensate,” Science 291, 2386 (2001).
  * Tao  _et al._ (2023) J. Tao, M. Zhao, and I. Spielman, “Observation of anisotropic superfluid density in an artificial crystal,” arXiv:2301.01258 (2023).
  * Dorrer and Zuegel (2007) C. Dorrer and J.D. Zuegel, “Design and analysis of binary beam shapers using error diffusion,” J. Opt. Soc. Am. B 24, 1268 (2007).
  * Corman  _et al._ (2017) L. Corman, J. L. Ville, R. Saint-Jalm, M. Aidelsburger, T. Bienaimé, S. Nascimbène, J. Dalibard, and J. Beugnon, “Transmission of near-resonant light through a dense slab of cold atoms,” Phys. Rev. A 96, 053629 (2017).
  * Zapata  _et al._ (1998) I. Zapata, F. Sols, and A.J Leggett, “Josephson effect between trapped Bose-Einstein condensates,” Phys. Rev. A 57, R28 (1998).
  * Note (3) We bypass the feedback loop technique reported in Ref.+.1667emZou  _et al._ (2021), which relies on LDA.

##  I Supplemental Material

###  I.1 Leggett’s bound and Gross-Pitaevskii equation

We consider a weakly interacting Bose-Einstein condensate (BEC), described by
Gross-Pitaevskii theory with the order parameter $\psi(\boldsymbol{r},t)$ and
the density $\rho(\boldsymbol{r},t)=|\psi(\boldsymbol{r},t)|^{2}$. We address
the case investigated in the main text, where the BEC is confined in a box of
length $L$ in the presence of an external potential $V(x)$. The equilibrium
density $\rho_{\rm eq}(x)$ thus depends only on the variable $x$. We assume
that the potential $V(x)$ satisfies the periodicity condition $V(L)=V(0)$ at
the boundary of the box.

We provide below a direct proof that Leggett’s upper bound

$f_{s,x}\leq\left(\langle\rho_{\rm eq}\rangle\;\langle\frac{1}{\rho_{\rm
eq}}\rangle\right)^{-1}$ (11)

for the superfluid fraction, relative to a fluid moving along the
$x$-direction, reduces to an identity. Here the average value of a quantity
$f(x)$ is defined by

$\langle f\rangle=\frac{1}{L}\int_{0}^{L}f(x)\,dx.$ (12)

In the main text, we take $V(x)=V_{0}\cos(qx)$ which is spatially periodic
with period $\ell=2\pi/q$. In this case, the density $\rho_{\rm eq}(x)$ is
also periodic with period $\ell$. Then, provided there is an integer number of
periods $\ell$ in the box of length $L$, the definition of the average value
(12) for a periodic function $f(x)$ coincides with the definition of the main
text

$\langle f\rangle=\frac{1}{\ell}\int_{0}^{\ell}f(x)\,dx.$ (13)

In the frame moving with velocity $v_{0}$, _i.e.,_ in the presence of the
perturbation $-v_{0}P_{x}$ (see main text), the Gross-Pitaevskii equation
(GPE) gives rise to the following expression for the equation of continuity
(only the motion along the $x$-axis is considered here)

$\partial_{t}\rho(x,t)+\partial_{x}\left\\{\rho(x,t)\left[v(x,t)-v_{0}\right]\right\\}=0\;,$
(14)

where $v(x)=(\hbar/m)\partial_{x}\phi(x)$ is the superfluid velocity field,
fixed by the gradient of the phase of the order parameter. In the limit of
small $v_{0}$, one can replace the density $\rho(x,t)$ with the equilibrium
value $\rho_{\rm eq}(x)$.

We look for a stationary solution in the moving frame by imposing a stationary
flow $\rho_{\rm eq}(x)(v(x)-v_{0})=j$, corresponding to a position and time
independent current $j$. We then obtain the result

$\frac{\hbar}{m}\partial_{x}\phi(x)=v(x)=v_{0}+\frac{j}{\rho_{\rm eq}(x)}$
(15)

for the gradient of the phase of the order parameter. Imposing that the phase
satisfies the periodic boundary condition $\phi(L)=\phi(0)$, the integration
of Eq.(15) gives

$0=v_{0}L+j\int_{0}^{L}\frac{dx}{\rho_{\rm eq}(x)}$ (16)

giving the current

$j=-\frac{v_{0}}{\langle\frac{1}{\rho_{\rm eq}}\rangle}.$ (17)

Finally, we calculate the average value of the momentum operator using Eqs.
(15,17):

$\displaystyle\langle\hat{P}_{x}\rangle$ $\displaystyle=$ $\displaystyle
m\int\rho_{\rm eq}(x)\,v(x)\,dx=mL\left(v_{0}\langle\rho_{\rm
eq}\rangle+j\right)$ (18) $\displaystyle=$ $\displaystyle
Nmv_{0}\left(1-\frac{1}{\langle\rho_{\rm eq}\rangle\;\langle\frac{1}{\rho_{\rm
eq}}\rangle}\right)$

where we introduced the number of particle $N=L\langle\rho_{\rm eq}\rangle$.
Using the definition of the main text:

$f_{s,x}=1-\lim_{v_{0}\to 0}\frac{\langle\hat{P}_{x}\rangle}{Nmv_{0}},$ (19)

we straightforwardly obtain the announced result for the superfluid fraction:

$f_{s,x}=\left(\langle\rho_{\rm eq}\rangle\;\langle\frac{1}{\rho_{\rm
eq}}\rangle\right)^{-1}.$ (20)

###  I.2 Superfluid density and effective mass

In this section, we explain why the two seemingly unconnected quantities,
superfluid fraction and effective mass, are connected when one considers a
zero-temperature Bose gas placed in a periodic potential and described by the
GPE for the order parameter $\psi(x)$. We assume that the potentiel $V(x)$ is
spatially periodic with period $\ell$ and that there are an integer number of
periods $n_{p}=L/\ell$ in the box of length $L$.

#### Superfluid density.

As explained in the previous section, the superfluid density can be calculated
using linear response theory for the perturbation $-v_{0}\hat{P}_{x}$.
Equivalently, one can apply twisted boundary conditions along the $x$ axis,
$\psi(L)=e^{i\theta}\,\psi(0)$ (with $\theta$ arbitrary small), and look for
the increase of energy of the ground state $\psi(x)$ at second order in
$\theta$ Leggett (1970). This energy increase directly provides $f_{s,x}$:

$E(\theta)\approx E(0)+Nf_{s,x}\frac{\hbar^{2}\theta^{2}}{2mL^{2}}.$ (21)

Since the potential $V(x)$ is periodic, it is natural to assume that the phase
twist $\theta$ of the wave function $\psi(x)$ minimizing the energy is
uniformly distributed over all lattice periods, hence

$\psi\left(n\ell\right)=e^{i\theta/n_{p}}\,\psi\left((n-1)\ell\right)\qquad
n=1,\cdots,n_{p}.$ (22)

To determine the ground state $\psi(x)$, it is convenient to perform the gauge
transform $\phi(x)=e^{-i\theta x/L}\,\psi(x)$, so that we recover periodic
boundary conditions for $\phi(x)$ with an additional vector potential. More
precisely, $\phi(x)$ is the solution of the equation

$\frac{\left(\hat{P}_{x}+\hbar\theta/L\right)^{2}}{2m}\phi(x)+V(x)\phi(x)+g|\phi(x)|^{2}\phi(x)=\mu\phi(x)$
(23)

minimizing the GP energy functional with the boundary condition deduced from
Eq. (22):

$\phi\left(n\ell\right)=\phi\left((n-1)\ell\right),$ (24)

_i.e.,_ a periodic boundary condition over each lattice period. This procedure
is directly connected to the one of the previous section by setting
$v_{0}=\hbar\theta/mL$ and leads to the result (20).

#### Effective mass.

We now turn to the determination of the effective mass for a BEC placed in the
periodic potential $V(x)$. We look for solutions that can be written as Bloch
functions, $\psi_{k}(x)=e^{ikx}\,u_{k}(x)$, where $u_{k}(x)$ is periodic over
the lattice and satisfies

$u_{k}\left(n\ell\right)=u_{k}\left((n-1)\ell\right).$ (25)

The equation satisfied by $u_{k}(x)$ is

$\frac{(\hat{P}_{x}+\hbar
k)^{2}}{2m}u_{k}(x)+V(x)u_{k}(x)+g|u_{k}(x)|^{2}u_{k}(x)=\mu\,u_{k}(x)$ (26)

and the effective mass $m_{x}^{*}$ is defined by the curvature of the energy
around $k=0$:

$E(k)\approx E(0)+N\frac{\hbar^{2}k^{2}}{2m_{x}^{*}}.$ (27)

Note that we define here the effective mass using the solutions of the non-
linear GPE. One can show that the same value for $m^{*}_{x}$ arises if one
considers the Bogoliubov excitations of the gas in the periodic potential
$V(x)$, i.e. the linearized version of the problem (see Eq. (28) of Taylor and
Zaremba (2003)).

The mathematical similarity between the two problems (23-24) and (25-26) is
clear, with the replacement $\theta/L\leftrightarrow k$, and the comparison
between Eqs. (21) and (27) provides the relation

$m_{x}^{*}=m/f_{s,x}$ (28)

for the zero-temperature Bose gas described by the GPE.

(a) (b) (c) (d)

Figure 4: (a) Typical dithered pattern displayed on the DMD. Each pixel is
either white (1) or black (0). (b) Cut of the dithered pattern (circles) along
with a fit (solid line) of the data to the function
$x\mapsto\sqrt{a_{0}+a_{1}\cos\left(kx\right)}$. (c) Measured image of the
projected lattice on the auxiliary camera. The residual intensity modulation
along $y$ is an artifact of the optical measurement due to parasitic fringes.
The axes are rescaled to coincide with the effective position in the atoms
plane. (d) Cut of the auxiliary camera image along the $x-$axis along with a
fit to the function $x\mapsto b_{0}+b_{1}\cos\left(kx\right)$.

###  I.3 Preparation and excitation of the gas

Our experimental setup has been described in Ref. Ville  _et al._ (2017). We
load an optical box potential made of a combination of laser beams at a
wavelength of 532 nm. The confinement along the vertical direction $z$ is
created by an optical lattice and we load a 3D Bose-Einstein condensate in a
single node of this lattice. The resulting trapping potential is well
approximated by an harmonic potential of frequency $\omega_{z}\approx
2\pi\times 3.7\,$kHz. The in-plane confinement is a square box potential of
size $L=40\,\rm$\mathrm{\SIUnitSymbolMicro m}$$ obtained by imaging the chip
of a digital micro-mirror device (DMD) on the atomic plane with a high-
resolution microscope objective. Atoms are cooled down thanks to evaporative
cooling to temperatures $T\lesssim 20\,$nK. For all experiments reported in
this work, the average surface density is fixed to $\rho_{0}=60(3)$
atoms$/\rm$\mathrm{\SIUnitSymbolMicro m}$^{2}$. Atoms are polarized in the
electronic ground state $\left|F=1,m=-1\right\rangle$. The in-situ density
distribution is measured using partial transfer to the
$\left|F=2\right\rangle$ state, followed by absorption imaging onto a CCD
camera with a typical resolution of $\lesssim
1\,\rm$\mathrm{\SIUnitSymbolMicro m}$$. The method of partial transfer allows
us to control the amount of atoms sensitive to the imaging beam so as to
maintain the optical depth of the imaged cloud to a low value. In this regime,
we avoid any collective effects in the light diffusion by the atomic cloud
Corman  _et al._ (2017).

The excitation of sound waves is performed by applying a magnetic field
gradient over the size of cloud, which creates a uniform force along either
the $x$ or the $y$ direction. In this work, we use a gradient of
$b^{\prime}=3.4\,\mu$G/$\mathrm{\SIUnitSymbolMicro m}$. For a chemical
potential of $\mu_{0}/k_{B}\sim 50\,$nK of the cloud, this corresponds to a
peak-to-peak potential of 0.09 $\mu_{0}$. The speed of sound is measured by
suddenly removing this gradient at time $t=0$ and then monitoring the time
evolution of the center-of-mass distribution of the cloud. The typical
amplitude of oscillation of the center-of-mass along the axis of the
excitation is $\approx 1\,\rm$\mathrm{\SIUnitSymbolMicro m}$$ and no clear
oscillation signal is observed along the orthogonal direction.

###  I.4 Optical lattice potential

We describe in this paragraph the preparation and the characterization of the
sinusoidal optical potential applied to the atomic cloud. This potential is
created by imaging the chip of a DMD onto the atomic plane with an imaging
system of resolution $\lesssim 1\,\rm$\mathrm{\SIUnitSymbolMicro m}$$. The
pixel size of the DMD is 13.7 $\mathrm{\SIUnitSymbolMicro m}$ and the
magnification is 1/70 leading to an effective pixel size on the atomic cloud
of $\approx\,$0.2 $\mathrm{\SIUnitSymbolMicro m}$, much smaller than the
optical resolution of the imaging system. The laser wavelength is 532 nm and
the beam waist in the atomic plane is $\sim 85\,\rm$\mathrm{\SIUnitSymbolMicro
m}$$ so that the intensity profile over the atomic cloud size is almost
uniform.

As described in Ref. Zou  _et al._ (2021), we use a dithering algorithm to
create effective “grey levels" of light intensity on the atomic plane 333We
bypass the feedback loop technique reported in Ref. Zou  _et al._ (2021),
which relies on LDA.. The lattice optical potential experienced by the atoms
is measured with an auxiliary imaging system with a similar resolution
$\lesssim 1\,\rm$\mathrm{\SIUnitSymbolMicro m}$$ and a magnification $\sim
3.2$. As the DMD operates as an amplitude modulator Dorrer and Zuegel (2007),
a sinusoidal modulation of the intensity on the atom plane will be obtained
using a target profile for the dithering algorithm given by

$f(x)=\sqrt{\alpha+\beta\cos(kx)}.$ (29)

The coefficients $\alpha$ and $\beta$ are chosen so that $\forall x$,
$0<f(x)$. For the data reported in this work, obtained with a lattice of
period $d\approx 4\,\rm$\mathrm{\SIUnitSymbolMicro m}$$, we used $\alpha=0.5$
and $\beta=0.4$.

We show in Fig. 4 a typical example of a dithered image that we display on the
DMD. We also report the corresponding intensity distribution over the size of
the cloud measured with the auxiliary imaging system. The obtained lattices
are well sinusoidal with an almost constant amplitude over the size of the
cloud. We also checked that the maximal amplitude of the lattice is linear
with the laser light intensity.

We now discuss the calibration of the lattice depth $V_{0}$. We shone large-
period lattices ($d>8\,\rm$\mathrm{\SIUnitSymbolMicro m}$$) on the atomic
cloud, we determined the dominant term of the density modulation $\rho_{1}$
and we deduced $V_{0}$ using the LDA prediction $\rho_{1}=V_{0}/g$. To make
this approach relevant, we restrict ourselves to low lattice depths
$V_{0}<\mu_{0}$ and we checked that, in this range and for a fixed laser
intensity, the contrast of the density modulation
$\mathcal{C}=\rho_{1}/\rho_{0}$ is independent of the lattice period. For such
large periods, our finite imaging resolution does not significantly affect the
measured signal, contrary to the case of shorter periods. Its measured value,
$\mathcal{C}\approx 0.75$, is consistent with the chosen values of $\alpha$
and $\beta$, leading to $\mathcal{C}=0.8$. Note that for large-period lattices
($\epsilon_{q}\ll\mu_{0}$), the correction to the LDA discussed in the main
text (see Eq. 8) is also negligible.

###  I.5 Role of finite optical resolution on the computation of Leggett’s
formula

For strong enough lattice depths, $V_{0}\gtrsim\mu_{0}$, the density
modulation $\rho_{n}(x)$ differs from a sinusoidal profile. While the
determination of the superfluid fraction through Leggett’s formula is still
exact (in the regime where the GPE is valid), its experimental determination
becomes more challenging. Deviations to a sinusoidal profile are described by
the emergence of higher-order harmonics $\rho_{n}$ in the density modulation,
which we have defined through $\rho_{0}-\sum_{n>0}\rho_{n}\cos(nqx)$. These
harmonics are subjected to spatial filtering by our imaging system. This
filtering is modeled by the $\beta_{n}$ coefficients in the main text.

We show in Fig. 5(a) the value of Leggett’s formula computed from the ground
state of the GPE taking into account this spatial filtering. In addition to
the LDA (thin solid line) and the GPE predictions (thick solid line), we show
the results obtained from the GPE keeping only the first $n=1$ harmonic
(dotted line), the two first harmonics (dot-dashed line) and the three first
harmonics (dashed line). With only two harmonics of the density modulation,
the estimated superfluid fraction is very close to the exact prediction for
$f_{s,x}\gtrsim 0.3$ and slightly deviates for lower superfluid fractions.
Including the third harmonic leads to an excellent approximation down to
$f_{s,x}\sim 0.1$. In Fig. 5(b), we compare our measurement of the density
modulation including only the first harmonic together with the corresponding
prediction of the GPE with filtering. The agreement is excellent. We also note
that for $V_{0}/\epsilon_{q}\gtrsim 8$ the computed superfluid fraction is
null because $\rho_{1}>\rho_{0}$, which leads to a vanishing density
$\rho_{0}-\rho_{1}\cos(qx)$ around the lattice maxima. Finally, we show in
Fig. 5(c) the data reported in the main text and compare them with the
prediction of the GPE when keeping the two first harmonics ($\rho_{1}$ and
$\rho_{2}$). Here again the agreement is very satisfactory. Note that the
observation of the third harmonic is not accessible with our setup. It
corresponds to a density modulation with a period $\sim
1.3\,\rm$\mathrm{\SIUnitSymbolMicro m}$$. To be detected, it would require a
numerical aperture (NA) of the optical system around 0.6 and is thus fully
filtered for our NA of 0.45. These results show that our choice of parameters,
and especially the lattice period, is suitable to obtain, within our
experimental uncertainties, a robust estimate of the superfluid fraction with
Leggett’s formula in a broad range of values of the lattice depth.

(a) (b) (c)

Figure 5: Influence of the finite optical resolution on the determination of
the superfluid fraction using Leggett’s formula. For all plots, we show the
LDA (thin solid line) and the GPE (thick solid line) predictions. In (a), the
predictions of the GPE with a filtering of the lowest harmonics are shown when
keeping only the first harmonic (blue dotted line), the two first harmonics
(red dot-dashed line) and the three first harmonics (green dashed line). (b)
Experimental determination of the Leggett’s formula keeping only the first
harmonic (blue circles). (c) Same when keeping only the two first harmonics,
as shown in the main text (red circles).


# Ice Core Dating using Probabilistic Programming

Aditya Ravuri
University of Cambridge
<EMAIL_ADDRESS>&Tom R. Andersson
British Antarctic Survey
<EMAIL_ADDRESS>&Ieva Kazlauskaite
University of Cambridge
British Antarctic Survey
<EMAIL_ADDRESS>&Will Tebbutt
University of Cambridge
<EMAIL_ADDRESS>&Richard E. Turner
University of Cambridge
<EMAIL_ADDRESS>&J. Scott Hosking
British Antarctic Survey
<EMAIL_ADDRESS>&Neil D. Lawrence
University of Cambridge
<EMAIL_ADDRESS>&Markus Kaiser
University of Cambridge
British Antarctic Survey
<EMAIL_ADDRESS>
###### Abstract

Ice cores record crucial information about past climate. However, before ice
core data can have scientific value, the chronology must be inferred by
estimating the age as a function of depth. Under certain conditions, chemicals
locked in the ice display quasi-periodic cycles that delineate annual layers.
Manually counting these noisy seasonal patterns to infer the chronology can be
an imperfect and time-consuming process, and does not capture uncertainty in a
principled fashion. In addition, several ice cores may be collected from a
region, introducing an aspect of spatial correlation between them. We present
an exploration of the use of probabilistic models for automatic dating of ice
cores, using probabilistic programming to showcase its use for prototyping,
automatic inference and maintainability, and demonstrate common failure modes
of these tools.

##  1 Introduction

Chemicals in the atmosphere are deposited onto ice sheets through
precipitation, with further deposition burying and eventually compacting the
snow into solid ice, recording the chemical composition of the atmosphere.
These chemicals provide evidence for the climate of the past and are known as
_proxy variables_. Annual cycles can be present in the data if 1) the
abundance of that proxy varies seasonally, 2) the precipitation rate at the
ice core site is large enough, and 3) the depth is not so great that the
annual layers have been excessively compressed. Given these conditions, annual
layer thickness is dictated by the amount of annual precipitation (a random
component) and compression of the ice with increasing depth (a systematic
component). A section of the data that we used, the Jurassic ice core from
Emanuelsson et al. [2022], is shown in Figure 1.

Constructing an ice core’s timescale manually through layer counting can be an
arduous process, lasting from days to years of person-time (Winstrup 2016).
Manual counting also has the disadvantage of poorly-quantified uncertainty in
the depth to time mapping, based on heuristics from expert disagreement or
uncertain layers. Probabilistic inference leads to a more principled treatment
of such uncertainties. A probabilistic approach also enables prior knowledge
to be incorporated during ice core timescale inference. For example, some
observations for time at certain depths may be available due to known volcanic
events depositing a layer of ash in the ice (these depth-time observations are
known as tie-points).

Our contribution is to present a case study in developing models for the ice
core dating problem using a Probabilistic Programming Language (PPL), mainly
Stan [Stan Development Team, 2022] alongside torchsde [Li et al., 2020], and
TensorFlow Probability [Dillon et al., 2017]. Our work highlights such tools
to practitioners at the intersection of probabilistic machine learning and
climate science. We aim to show the promise of PPLs in how they can enable the
composition of assumptions, easy experimentation, model extension and
maintainability. A core promise of widely used PPLs is that they automate and
abstract away inference details, thus ensuring that models are written at the
right level of abstraction without the need to write and maintain complex
inference algorithms. We also show current limitations of these tools, such as
inference only being computationally feasible over a limited set of models
that may be suitable for a task, and inference methods being limited in
practice depending on the specific functional requirements set out for each
PPL.

Figure 1:  Section of the data, showing annual seasonality in MSA (a proxy
exhibiting a strong seasonal pattern) plotted against depth. Manually counted
annual layers are overlaid as shaded bars.

##  2 Problem Setup

The proxy variables are measured in ice cores along a depth dimension. We
represent the depth series as a set of known random variables
$\boldsymbol{\delta}\coloneqq\\{\delta_{i}\\}_{i=1}^{n},\delta_{i}\in\mathbb{R}^{+}$
where $n$ is the number of sampled depths at which proxy readings are
available. The depth series $\boldsymbol{\delta}$ is mapped to latent time
(age) values associated with each observation, represented by the stochastic
process $\mathbf{t}\coloneqq\\{t_{\delta_{i}}\\}_{i=1}^{n}$ indexed by the
depth series. Proxy measurements are denoted as
$\mathbf{s}\coloneqq\\{s_{\delta_{i}}\\}_{i=1}^{n}$, with
$s_{\delta_{i}}\in\mathbb{R}$. The _ice core dating_ process can be stated as
an inference problem for time conditioned on the proxy and depth data,
$\mathbf{t}|\mathbf{s},\boldsymbol{\delta}$. Within this work, we assume that
proxies depend only on the latent time process, the depth to time mapping is
monotonic, and that there is only one proxy available with a clear seasonal
signal. This gives rise to the class of models shown in Fig. 4 (a), in the
Appendix.

##  3 Methodology

The setup of the problem gives rise to a class of models that vary in their
assumptions and complexity from discrete-index HMMs to continuous-index HMMs
to SDE-based models. Some of these models, particularly HMMs, were previously
studied by Winstrup [2011, 2016] who explored the the effect of batching, the
usage of other observation models and extensions to hidden semi-Markov models
(for allowing for priors to be set over lengths of year boundaries) and
inference therein.

###  3.1 A Hidden Markov Model

#### Model

Under the assumption that the depth-sampling is uniform, and that the latent
time process has a discrete domain, we can model the latent time process as a
Markov chain. If in addition, the proxies are conditionally independent given
the time periods they correspond to, the framework in Figure 4 (a) reduces to
a Hidden Markov Model shown in Figure 4 (b). We assume that $\boldsymbol{t}$
can occupy states $\forall
i:t_{\delta_{i}}\in\\{\nicefrac{{k}}{{n_{s}}}\\}_{k=1}^{m\cdot n_{s}}$, where
$n_{s}$ is the number of states within each yearly cycle and $m$ is an
arbitrarily large number of years. The transition matrix

$\mathbb{P}(t_{\delta_{i}}|t_{\delta_{i-1}})=\begin{bmatrix}p_{1/n_{s}}&1-p_{1/n_{s}}&0&\cdots&0\\\
0&p_{2/n_{s}}&1-p_{2/n_{s}}&\cdots&0\\\ 0&0&0&0&1\\\ \end{bmatrix}$ (1)

is bidiagonal, with $p_{i}$ denoting the probability that the chain stays in
state $i$ given that it was in the same state at the last depth measurement.
In other words, as depth increases, the time process either advances or stays
the same. This enforces monotonicity of the latent time sequence w.r.t. depth
and also allows one to track what year an observation corresponds to (given by
the floor of the state, $\lfloor i\rfloor$). The parameters for each state
within a year are repeated across years to constrain the model (i.e.
$p_{c+k/n_{s}}=p_{k/n_{s}}$ for $c\in\mathbb{N}$). The following observation
model is used as part of the HMM,

$\forall i:s_{\delta_{i}}|t_{\delta_{i}}\sim\mathcal{N}(a\cos(2\pi
t_{\delta_{i}})+b,\sigma^{2}).$

Note that, as the transition matrix parameters are not constant within each
annual layer, the model’s annual layer shapes can be a warping of the mean
cosine function in the observation model, enabling some flexibility in
modelling the real proxy cycle shapes.

#### Inference

Using the Stan language, we perform maximum likelihood inference for the
parameters given data with time marginalised
($a,b,\sigma,\\{p\\}_{j}|\mathbf{s}$), with the posterior over times
$\mathbf{t}|\mathbf{s},a,b,\sigma,\\{p\\}_{j}$ estimated using the forward-
backward algorithm in Stan. The probabilistic program is shown in Appendix B.
Two problems are encountered with standard tooling:

  * •

The runtimes increase quadratically w.r.t. the state space (due to the forward
algorithm, which computes the log likelihood with hidden states marginalised).
We remedy this by rewriting the forward algorithm exploiting the sparsity of
our transition matrix, providing a direct replacement of the native Stan
function (shown in Appendix B). In a PPL, users can typically implement their
own efficient functions for such likelihood calculations.

  * •

The model is misspecified as the data is expected to be non-stationary (due to
compression of the ice core and temporal variation in precipitation). This can
pragmatically be remedied by processing the data in batches, allowing the
parameters such as $a$ and $b$ to change between different sections of the ice
core.

Inference using this implementation takes 2.5 minutes using a single-thread
run, without GPU utilization, for the entire ice core (about 2.5k
observations), without requiring any special initialization strategies.

###  3.2 An extension allowing for tie-point specification

To allow for tie-points to be integrated, which constrain the depth-to-time
mapping, the observation model is extended to account for non-stationarity in
the signal, because accounting for tie-points would require processing data in
batches large enough to cover the tie-points. This is done by changing the
observation model to allow for parameters $a,b$ to change with each data point
along the ice core,

$\displaystyle\forall i:s_{\delta_{i}}|t_{\delta_{i}}$
$\displaystyle\sim\mathcal{N}(a_{i}\cos(2\pi
t_{\delta_{i}})+b_{i},\sigma^{2}),$

with a prior (as part of a hierarchical model) placed over $a_{i},b_{i}$
(thus, the change is “slow”/constrained). A similar hierarchical treatment is
given to parameters of the transition matrix, using the prior

$\forall j\in 1/n_{s},...,m:p_{j}\sim\text{Beta}(\alpha_{j*n_{s}\text{ mod
}n_{s}},\beta_{j*n_{s}\text{ mod }n_{s}}).$

This allows the transition probabilities to change over years (and hence
depths). The hierarchical distribution of states at the same point in any
yearly cycle however share the same parameters, providing some constraint.

Having extended the model to different observation models per data-point, we
specify volcanic tie-points using an alternative observation model using the
state directly rather than proxy information. We use a categorical
distribution ensuring that the time states must reach the volcanic tie-points,

$p(s^{\prime}_{\delta_{\text{tie}}}|t_{\delta_{\text{tie}}})=\begin{cases}1/n_{s}&\text{if
}\lfloor t_{\delta_{\text{tie}}}\rfloor=t_{\text{tie}}\\\
0&\text{otherwise}\end{cases},$

to enforce that the volcanic ash observed at depth $\delta_{\text{tie}}$
corresponds to a known year of the eruption $t_{\text{tie}}$.

Maximum likelihood inference for parameters such as $a_{i},b_{i}$ in such
models produces subpar results, due to posterior modes of such unidentifiable
models lying outside their typical sets, raising a need to integrate out the
parameters. As MCMC is computationally infeasible due to the number of
parameters and data, we use mean-field variational inference (a seamless
change in the Stan API). This mode of inference greatly increases runtimes
w.r.t. the model in Section 3.1 however, to about a few hours due to the
increased number of parameters. Results of inference in the two models
presented thus far is shown in Figure 3.

Figure 2: A comparison between the inference results, i.e. sample paths of
$\mathbf{t}|\mathbf{s}$, obtained using Stan. These correspond to models
presented in sections 3.1 and 3.2, with expert labels overlaid.

Figure 3: On top, a short data series with a missing section. At the bottom,
posterior paths of time inferred using the data above - the missing section
induces multimodality in the depth to time posterior.

###  3.3 An extension to continuous-index Hidden Markov Models

Depth measurements may not be regularly spaced, due to missing sections of an
ice core as a result of the extraction process, or due to the sampling
frequency changing over its depth. In this case, modelling the latent time
process as a continuous-index Markov chain leads to the times observed at
irregularly spaced points being described by an index-(depth) inhomogeneous
Markov chain. Here, the transition matrix associated with a step from
$\delta_{i-1}$ and $\delta_{i}$ can be computed via the transition rate matrix
$\mathbf{Q}$ with the same sparsity structure as Equation 1,

$\mathbb{P}(t_{\delta_{i}}|t_{\delta_{i-1}})=\exp_{\text{matrix}}((\delta_{i}-\delta_{i-1})\mathbf{Q}).$
(2)

Such models are termed _continuous-time HMMs_ [Liu et al., 2015] (although we
use the terminology continuous-index to avoid confusion as time is not the
index in our application). They allow a better representation of posterior
uncertainty arising from missing observations, as seen in Fig. 3. Inference
was performed on small datasets by using tensorflow-probability for computing
log likelihoods of HMMs in a differentiable manner. High-level pseudocode for
MLE/VI when using such tools is shown in Appendix E. Inference involving
larger datasets would involve bespoke inference code due to a lack of
functionality around efficient computation of the forward algorithm involving
matrix exponentials within the PPLs considered.

###  3.4 An exploration of Stochastic Differential Equations for ice core
dating

Assuming instead that we are working with a continuous index and a continuous
state space for $\mathbf{t}$, the class of models under consideration
naturally extends to state space/ SDE models. Given an index $\delta$, a prior
over a stochastic process $t_{\delta}$ can be formulated as,

$\begin{bmatrix}d\mathbf{z}_{\delta}\\\
dt_{\delta}\end{bmatrix}=\begin{bmatrix}\mu(\mathbf{z}_{\delta},\delta)\\\
-\exp(\mu^{\prime}(\mathbf{z}_{\delta},t_{\delta},\delta))\end{bmatrix}d\delta+\begin{bmatrix}\Sigma\\\
\epsilon\end{bmatrix}d\mathbf{W}_{\delta},$

with $\epsilon\rightarrow 0$ enforcing monotonicity of sample paths and where
$\mathbf{z}_{\delta}$ is a latent process that, for example, can have a GP
prior represented as an SDE [Särkkä and Solin, 2019]. Such a prior is similar
to the one used in Ustyuzhaninov et al. [2020]. Inference in this class of
models can be performed by specifying a variational SDE using the same
diffusion but with a different drift [Li et al., 2020]. As in the case above,
due to customised variational inference not being supported in Stan, we use
functionality from torchsde (differentiable solvers) that enable computation
of the objective.

Inference in this class of models was particularly difficult, mirroring the
findings of GPCore [Andersson, 2019], where a GP prior is placed on
$\mathbf{t}|\boldsymbol{\delta}$, and where maximum likelihood inference is
performed for $\mathbf{t}|\mathbf{s},\mathbf{d}$, in a constrained manner to
ensure monotonicity in the depth-time sample paths. The inference is very
dependent on good initialisation; using existing estimates of the chronology
in the case of GPCore, and using inverse Lomb-Scargle spectrograms for
$\mathbf{z}_{\delta}$ in our SDE models, following ideas from Gay et al.
[2014]. However, we had greater difficulties due to local minima than in
GPCore. More discussion on this topic can be found in Appendix F.

##  4 Conclusion

In this paper we presented a taxonomy of models for ice core dating, showing
how simplifying assumptions lead naturally to HMMs and how these models, and
corresponding probabilistic programs can be extended. We exemplified the
shortcomings of PPLs as the models are extended towards SDEs and discussed
difficulties with inference in such models. Future work would involve
extension of the models presented to utilize multiple proxies, account for
spatial correlation between ice cores, and an exploration to determine how to
aid inference in latent SDE models for such modelling tasks.

### Data and Code

Code to reproduce key results in this paper can be found at
https://github.com/infprobscix/icecores.

### Acknowledgements

AR is supported by a studentship from the Accelerate Programme for Scientific
Discovery. IK was funded by a Biometrika Fellowship awarded by the Biometrika
Trust. The authors would also like to thank Liz Thomas and Dieter Tetzner from
the British Antarctic Survey for providing the datasets we used, for their
insight and helpful discussions on ice core dating task.

## References

  * Andersson [2019] Tom Andersson.  GPCore: A Gaussian process approach for inferring ice core chronologies.  Master’s thesis, University of Cambridge, 2019.
  * Barber [2012] David Barber.  _Bayesian Reasoning and Machine Learning_.  Cambridge University Press, 2012.
  * Dillon et al. [2017] Joshua V. Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A. Saurous.  Tensorflow distributions, 2017.
  * Emanuelsson et al. [2022] B. Daniel Emanuelsson, Elizabeth R. Thomas, Dieter R. Tetzner, Jack D. Humby, and Diana O. Vladimirova.  Ice core chronologies from the antarctic peninsula: The palmer, jurassic, and rendezvous age-scales.  _Geosciences_ , 12(2), 2022.
  * Gay et al. [2014] M. Gay, M. De Angelis, and J.-L. Lacoume.  Dating a tropical ice core by time–frequency analysis of ion concentration depth profiles.  _Climate of the Past_ , 10(5), 2014.
  * Hensman et al. [2016] James Hensman, Nicolas Durrande, and Arno Solin.  Variational fourier features for gaussian processes, 2016.
  * Li et al. [2020] Xuechen Li, Ting-Kam Leonard Wong, Ricky T. Q. Chen, and David Duvenaud.  Scalable gradients for stochastic differential equations, 2020.
  * Liu et al. [2015] Yu-Ying Liu, Shuang Li, Fuxin Li, Le Song, and James M. Rehg.  Efficient learning of continuous-time hidden markov models for disease progression.  _Advances in Neural Information Processing Systems_ , 28:3599–3607, 2015.
  * Murphy [2012] Kevin P Murphy.  _Machine learning: a probabilistic perspective_.  2012\.
  * Stan Development Team [2022] Stan Development Team.  Stan modeling language users guide and reference manual 2.29, 2022.  URL https://mc-stan.org.
  * Särkkä and Solin [2019] Simo Särkkä and Arno Solin.  _Applied Stochastic Differential Equations_.  Institute of Mathematical Statistics Textbooks. 2019.
  * Ustyuzhaninov et al. [2020] Ivan Ustyuzhaninov, Ieva Kazlauskaite, Carl Henrik Ek, and Neill Campbell.  Monotonic Gaussian process flows.  In _International Conference on Artificial Intelligence and Statistics_ , pages 3057–3067. PMLR, 2020.
  * Winstrup [2011] Mai Winstrup.  _An Automated Method for Annual Layer Counting in Ice Cores_.  PhD thesis, University of Copenhagen, 2011.
  * Winstrup [2016] Mai Winstrup.  A hidden markov model approach to infer timescales for high-resolution climate archives.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence_ , AAAI’16, pages 4053–4060. AAAI Press, 2016.

##  Appendix A Appendix: Supplementary Figures

$\boldsymbol{\delta}$$\mathbf{t}$$\mathbf{s}$

$s_{\delta_{1}}$$s_{\delta_{2}}$$...$$t_{\delta_{1}}$$t_{\delta_{2}}$$...$

Figure 4:  (a.) On the left, a graph that summarises the class of models under
consideration. (b.) On the right, a graph representing a hidden Markov model
obtained when assumptions of the form in Section 3.1 are made.

##  Appendix B Appendix: Discrete-index HMMs in Stan

Probabilistic programs in Stan are composed as blocks corresponding to,

  * •

data and transformed data: where users can input / specify random variables
for which observations have been made and other fixed values,

  * •

parameters and transformed parameters: where users specify random variables on
which inference is performed,

  * •

model: where users specify distributional assumptions on and between the
random variables defined,

  * •

generated quantities: where downstream analysis can be computed as a function
of posterior draws of the parameters conditioned on data.

Extensive documentation on example models in the Stan language, on functions
and the language can be found here

We define the data blocks in Stan as follows, for all use cases.

A simple HMM in Stan, corresponding to the model defined in Section 3.1, would
consist of parameters,

on which inference will occur. Following the stan syntax on how HMMs are
defined, we define log_omega corresponding to our observation model, and other
variables as follows,

Note that we create these variables in the transformed parameters section
instead of the model section (which would be more Stan like, as it would be
easier to read the observation model in the model section) to be able to
access these variables in the generated quantities block without having to
redefine them. It’s an unusual case syntactical case as that the parameters
governing the likelihood of the HMM with hidden states marginalised is a
function of the observation model likelihood. Nevertheless, the observation
model is specified in the line,

The model then simply specifies the log posterior, which is a sum of priors
over our parameters (implicitly assumed to be uniform as they’re not
specified) and the likelihood of

$\mathbf{s}|\mu,\sigma,\text{scale and transition matrix diagonal}$

which is specified as,

where target corresponds to the log posterior.

As we’re interested in the posterior probabilities (for analysis and to
initialize the inital state probabilities in an iterative manner if this model
is used on sequential batches of data) and sample paths, we generate the
following quantities,

for post-hoc analysis.

As the forward algorithm here, i.e. the function hmm_marginal is inefficient
for our use case, we can define a more efficient function for bidiagonal
transition matrices as follows. A discussion on the computation of the forward
algorithm can be found in Barber [2012], Murphy [2012].

Then, the model changes to,

Other convenience functions we define are as follows,

##  Appendix C Appendix: discrete-index HMMs in Stan, hierarchical case

To change the model to allow for changing parameters as described in 3.2, we
simply change the following aspects of the Stan model.

The parameters are extended such that there is one for every data point, with
further parameters corresponding to hierarchical distribution parameters.

The observation model in the transformed parameters section changes to account
for the varying parameters,

The model block changes to account for the hierarchical priors,

We do not change the generated quantities block (and do not code up a more
efficient version of the forward backward algorithm as the generated
quantities block is only run once at the end of the inference process unlike
the forward algorithm, which would be run multiple times during the inference
process as the parameters are changed).

We change the inference method to variational inference instead of maximum
likelihood due to the need to integrate out the hierarchical parameters during
inference (as a maximum likelihood estimate of parameters such as
$a_{i},b_{i}$ in such an overparameterised model is perhaps unlikely to lie in
the posterior’s typical set). This is done by replacing model.optimize with
model.variational. We also increased (doubled with respect to Stan’s defaults)
the number of gradient samples for the ELBO calculation, which was needed to
obtain reasonable outputs.

##  Appendix D Appendix: cts-HMM forward algorithms in Stan

A computationally inefficient example allowing for time varying transition
matrices in the forward algorithm is shown below. In this example, the time
varying transition matrix is created using Equation 2.

##  Appendix E Appendix: other models

We implemented the continuous time Markov chain model in tensorflow-
probability due to out of the box support for time varying transition
matrices. We implemented our SDE models using torchsde due to support for
differentiable SDE solvers.

The basic algorithm followed in both cases (as VI and maximum likelihood
estimation are both optimization problems) is,

$\texttt{log\\_posterior\\_or\\_lower\\_bound(params)}\leftarrow\texttt{function(params,
data) ...}$

$\texttt{parameters}\leftarrow...$ $\triangleright$ variational params or
params on which inference is done

$\texttt{optimizer}\leftarrow\texttt{Optimizer(parameters, lr)}$

while loss not converged do

$\texttt{loss}\leftarrow-\texttt{log\\_posterior\\_or\\_lower\\_bound(parameters)}$

$\texttt{grad}\leftarrow\texttt{loss.grad()}$

$\texttt{params}\leftarrow\texttt{optimiser.step(params, grad)}$

end while

posthoc analysis

##  Appendix F Appendix: Results of SDE models

Our SDE models ran into severe difficulties with inference. As an example, we
assume the prior,

$\displaystyle\begin{bmatrix}dz^{a}_{\delta}\\\ dz^{b}_{\delta}\\\
dt_{\delta}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}z^{b}_{\delta}\\\
-\lambda^{2}z^{a}_{\delta}-2\lambda z^{b}_{\delta}\\\
\alpha*\sigma^{+}(z^{a}_{\delta})\end{bmatrix}d\delta+\begin{bmatrix}0\\\ 1\\\
10^{-2}\end{bmatrix}\odot d\mathbf{W}_{\delta},$

where $\alpha$ and $\lambda$ are set to be constant, and $\sigma^{+}$
corresponds to the softplus operation. Note that the prior over $\mathbf{z}$
is a Matérn-3/2 Gaussian process. We use the following observation model,

$s_{\delta_{i}}|t_{\delta_{i}}\sim\text{Laplace}(\sin(\pi
t_{\delta_{i}}),0.05).$

We tried to fit a variational posterior over $t_{\delta_{i}}$, using the
variational SDE,

$\displaystyle\begin{bmatrix}dz^{a}_{\delta}\\\ dz^{b}_{\delta}\\\
dt_{\delta}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}f^{a}_{nn}(\mathbf{z}_{\delta};t_{\delta})\\\
f^{b}_{nn}(\mathbf{z}_{\delta};t_{\delta})\\\
\alpha*\sigma^{+}(f^{c}_{nn}(\mathbf{z}_{\delta};t_{\delta}))\end{bmatrix}d\delta+\begin{bmatrix}0\\\
1\\\ 10^{-2}\end{bmatrix}\odot d\mathbf{W}_{\delta},$ (3)

where $\mathbf{f}_{nn}$ was parameterized using a neural network, to some data
simulated from the prior. The expected mean of the observation model for a few
samples from the variational posterior are shown below in Figure 5.

Figure 5:  Results of our latent SDE models on synthetic data, showing a poor
local maximum in the ELBO reached by VI.

We also tried to use sparse Gaussian processes for placing priors and
variational approximations on $\mathbf{z}_{\delta}$, utilizing random Fourier
features [Hensman et al., 2016] to sample functions from these GPs to work
seamlessly with torchsde, finding no improvement in results. Pre-fitting a
neural-network based SDE prior (using spectra derived from the data) however
results in a better fit to the ice core data (illustrated below in Figure 6),
however, results are still very (impractically) sensitive to initialisation.

Figure 6:  Results of our latent SDE model using spectra to initialize the
prior, showing a relatively good fit to the ice core data. The blue line shows
ground truth data, while the grey lines show the mean of the observation model
using different samples of the posterior over time.

Future work can also involve exploration of Kalman filtering algorithms for
this problem, as our SDE priors are valid models that may be assumed for state
estimation in extended Kalman filtering.