# Leidenfrost Flows : instabilities and symmetry breakings

E. Yim1,†, A. Bouillant2,3,†, D. Quéré2,3, F. Gallaire1 1LFMI, École
Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland
2LadHyX, École polytechnique, 91128 Palaiseau, France
3PMMH, PSL-ESPCI, CNRS-UMR 7636, 75005 Paris, France
$\dagger$ Authors contributed equally to this work

(Received: XX 2021; Revised: XX XX 2021; Accepted: XX XX 2022)

###### Abstract

Leidenfrost drops were recently found to host strong dynamics. In the present
study, we investigate both experimentally and theoretically the flows
structures and stability inside a Leidenfrost water drop as it evaporates,
starting with a large puddle. As revealed by infrared mapping, the drop base
is warmer than its apex by typically 10∘C, which is likely to trigger bulk
thermobuoyant flows and Marangoni surface flows. Tracer particles unveil
complex and strong flows that undergo successive symmetry breakings as the
drop evaporates. We investigate the linear stability of the baseflows in a
non-deformable, quasi-static, levitating drop induced by thermobuoyancy and
effective thermocapillary surface stress, using only one adjustable parameter.
The stability analysis of nominally axisymmetric thermoconvective flows,
parametrized by the drop radius $R$, yields the most unstable, thus, dominant,
azimuthal modes (of wavenumber $m$). Our theory predicts well the radii $R$
for the mode transitions and cascade with decreasing wavenumbers from $m=3$,
$m=2$, down to $m=1$ (the eventual rolling mode that entails propulsion) as
the drop shrinks in size. The effect of the escaping vapor is not taken into
account here, which may further destabilize the inner flow and couple to the
liquid/vapor interface to give rise to motion (Bouillant et al., 2018b;
Brandão & Schnitzer, 2020).

Leidenfrost; Flow stability; Thermo-driven convection

##  I Introduction

A water drop can levitate above a hot surface provided the solid temperature
exceeds the boiling temperature of the liquid by typically 100∘C. This effect
was first reported in 1756 by J.G. Leidenfrost (Leidenfrost, 1756) and it has
been ever since a source of scientific curiosity. Evaporation produces a vapor
film underneath the drop with typical thickness 50 $\mu$m that thermally
insulates the liquid from its substrate and suppresses boiling. The drop can
thus survive a few minutes to plate temperatures as high as 300∘C. Levitation
also prevents the liquid from wetting the surface. As a consequence, a drop
adopts a quasi-spherical shape when its radius $R$ is smaller than the
capillary length $\ell_{c}=\sqrt{{\sigma_{0}}/{\rho_{0}g}}$ $-$ where
$\sigma_{0}$ is the surface tension, $\rho_{0}$ is the liquid density and $g$
is the acceleration of gravity, that is $\ell_{c}\approx 2.5$ mm for water at
100∘C $-$ while it gets flattened by gravity to a height 2$\ell_{c}$ when
$R>\ell_{c}$. The absence of contact also produces a (quasi-)frictionless
situation by which the liquid drop becomes highly mobile. Thermal and
mechanical insulation have contrasted consequences. On the one hand, above the
Leidenfrost temperature, vapor compromises the liquid cooling properties and
pilots the transition from nucleate to film boiling, unwelcome in nuclear
reactor or in metallurgy (loss of control in quenching processes). On the
other hand, it produces a purely non-wetting situation, which can serve as a
canonical model for superhydrophobicity. Recent studies have illustrated the
richness of spontaneous dynamics related to Leidenfrost drops (Quéré, 2013),
including oscillations (Garmett, 1878; Holter & Glasscock, 1952; Brunet &
Snoeijer, 2011; Bouillant et al., 2021a), bouncing (Celestini et al., 2012;
Waitukaitis et al., 2017) and directed propulsion on asymmetrically textured
(Linke et al., 2006; Cousins et al., 2012) or on non-uniformly heated surfaces
(Sobac et al., 2017; Bouillant et al., 2021b), which could be exploited in the
design of micro-reactors (Raufaste et al., 2016), of heat pipes and
exchangers, as well as for lab-on-a-chip technologies.
Among those dynamics is the ability of Leidenfrost droplets to self-rotate and
self-propel, even in the absence of external fields (Bouillant et al., 2018b).
Despite the early discovery of the Leidenfrost effect, the existence of strong
inner flows has only been reported recently, with velocities as high as a few
cm/s and whose morphology switches from four counter-rotating swirling cells,
preserving the axial symmetry, to an unique asymmetric swirl as evaporation
proceeds (Bouillant et al., 2018a). The eventual solid-like asymmetric
rotation comes with a tilt at the droplet base by an angle $\alpha$ of
typically a few milliradians (Bouillant et al., 2018b), producing
accelerations that scale as $a\sim\alpha g$ (Dupeux et al., 2013) of a few
tens of mm/s2. The internal rolling actually couples with the vapor cushion in
a feedback loop to sustain the self-rolling and self-propelling motions
(Brandão & Schnitzer, 2020), which thereby explains the intrinsic mobility of
Leidenfrost droplets. The origin of inner flows as well as their structure
remain however unclear. They may arise from: i) a thermal scenario: the base
of the drop is close to the hot plate while its apex is exposed to cooler air.
Temperature differences in the liquid may give birth to convective flows,
either driven by thermocapillary effects, known as Marangoni effect
(variations of the surface tension along the temperature gradient as in
Scriven & Sternling (1960)) or driven in the liquid bulk by thermobuoyant
effects, known as Rayleigh-Bénard convection (resulting from the thermal
expansion of the liquid); ii) a hydrodynamic scenario: vapor produced at the
drop base escapes through the subjacent film exerting a viscous radial stress
that may draw liquid with it. Both scenarios a priori preserve axisymmetry,
which is contradicted by experimental observations. Indeed, Particle induced
velocimetry (PIV) has revealed that internal flow structures evolve in time
and therefore with the drop geometry (Bouillant et al., 2018b), suggesting a
symmetry selection mechanism as the drop shrinks in size. Similar selection
mechanism has been reported in a non-levitating drop deposited on a warm plate
(Tam et al., 2009; Dash et al., 2014). In this close configuration, the vapor
cushion is suppressed and the origin of the internal dynamics is purely
thermal. Temperature gradients being greater in Leidenfrost drop than for
sessile drop on warm plates, we anticipate enhanced thermoconvection,
prompting us to focus on the thermal scenario (i) rather than on the
hydrodynamics scenario relying on the escaping vapor (ii). We therefore
address the stability and the symmetry of thermo-induced flows inside a
Leidenfrost drop as it evaporates, neglecting the hydrodynamic effect from the
escaping vapor and without explicitly modelling the surrounding gas. We aim at
predicting the radii at which are expected the successive symmetry breakings
starting with $R<4\ell_{c}$ (to prevent chimney formation), that is in the
regime where drops are flatten by gravity ($R>\ell_{c}$), while approaching
the regime where drops get quasi-spherical ($R<\ell_{c}$). Note that we do not
specifically explore the transition to the rolling mode $m=1$ as in Bouillant
et al. (2018b); Brandão & Schnitzer (2020) but we explore higher modes $m>1$
observed in large puddles $R>\ell_{c}$. In this regime, vapors carve a blister
underneath the liquid, whose amplitude increases with $R$ and reaches the
entire height of the puddle, $2\ell_{c}$, when $R\lesssim 4\ell_{c}$. Viscous
entrainment being markedly weakened along this vapor pocket, vapors would only
draw liquid along the very narrow, peripheral neck (Pomeau et al., 2012;
Burton et al., 2012; Sobac et al., 2014). The peculiar geometry of the vapor
cushion underneath large puddles undermines the hydrodynamic scenario,
corroborating our assumptions. Hydrodynamic effects of the vapor cushion may
however shift our prediction for the bifurcation radii, especially when
$R\rightarrow\ell_{c}$. As shown by Brandão & Schnitzer (2020) in the case of
quasi-spherical drops, the internal rolling couples with the vapor cushion in
a feedback loop, which sustains the self-rolling and self-propelling motions.
Such a coupling could also reshape the vapor cushion when $R>\ell_{c}$, with
consequences on the drop mobility that are not captured by our model.
The selection mechanism for the flow symmetry is a common feature in both
buoyancy and Marangoni instabilities. These instabilities are known to be very
sensitive to the geometry and confinement; the preferred discrete azimuthal
mode wavenumber decreases with the liquid domain size. The stability analysis
on the Marangoni-Rayleigh-Bénard convection has been extensively studied for a
rectangular and cylindrical domain with various boundary conditions Pearson
(1958); Nield (1964); Vrentas et al. (1981); Rosenblat et al. (1982); Kuhlmann
& Rath (1993); Johnson & Narayanan (1999). Yet, the stability of a free liquid
drop subjected to a vertical temperature gradient has heretofore not received
the same attention. Therefore, Leidenfrost drops constitute a toy model where
the evaporation-driven confinement enables to quasi-steadily sweep the states
in a non-wetting drop heated from below. We restrict our parametric study in
$R$ to the limit of $R<5$ mm to prevent Leidenfrost chimneys and pulsating
stars to appear (Quéré, 2013). Leaning on the experimental observations, we
develop a numerical model, which implicitly decouples the evaporation
timescale from the meanflow evolution timescale. We thus look for the
stability of the nominally axisymmetric thermo-convective baseflow in
Leidenfrost drops in order to explain the symmetry breaking from 4 to 1
convective cells reported in Bouillant et al. (2018a), as well as the prior
transition from 6 to 4 cells.
We first characterize the successive symmetry breakings in the internal flows
(§II) and extract from experiments physical quantities relevant to the problem
such as the temperature difference at the drop interface. The governing
hydrodynamic and thermal equations, as well as the linear stability analysis
are then presented in §III and the results are shown in §IV. We obtain the
baseflow generated by the stratification within the liquid and search for the
effective Marangoni number, which best captures our experimental observations,
such as the surface temperature and the velocity field. Then, a stability
analysis is carried out for different drop radii $R$. We compare for a given
$R$ the stability properties of each symmetry breaking mode, the mode with the
highest growth rate – the most unstable mode – being expected to dominate the
flow structure. Our study predicts the successive inner flow symmetries as
drops shrink in size. It also provides the critical radii for the modes
transition in quantitative agreement with observations. We eventually discuss
and compare the numerical outcomes to experiments, and add a few concluding
remarks in §LABEL:sec:discuss.

##  II Experimental observations

###  II.1 Quasi-steady state assumption

Leidenfrost drops levitate above a thin layer of vapor, of good insulating
properties since $k_{a}<<k$, where $k_{a}$ and $k$ are the air and water
thermal conductivities, provided in Table 1 of the SI. Evaporation, which
mainly takes place at the drop base, is thus markedly reduced and we verify
here that the drop is at quasi-static equilibrium. A water drop with initial
radius $R_{0}$ initially close to $\approx 4$ mm is deposited on a plate
brought to 300∘C. We use a slightly curved substrate to immobilize the highly
mobile liquid. The drop is observed using a top-view high speed camera, from
which we extract $R$, the drop equatorial radius as evaporation proceeds.
Figure 1 shows that $R$ decreases linearly with time $t$ at a rate
$\mathrm{d}R/\mathrm{d}t\approx-22\rm\;\mu m/s$, so that the drop survives
about $\tau_{0}\approx 3$ minutes.

Figure 1:  (a) Radius $R$ of a Leidenfrost drop levitating on a plate heated
at 350∘C as a function of time $t$. $R$ decreases as
$R(t)=R_{0}(1-t/\tau_{0})$ (eq.1), plotted as dotted line, denoting
$R_{0}=3.7\pm 0.1$ mm as the initial radius and $\tau_{0}=176.5$ s as the
lifetime. (b) Drop shape for some radii $R\in[0.9;4.5]$ mm readable in (a)
(see colored dots). Simulated shapes (full lines), obtained by numerically
integrating (2), are compared to experimental ones (dotted lines), obtained
for side-viewed drops pinned with a needle. $\ell_{c}$ is the capillary
length. Surface flows, viewed from the top, successively self-organize into
(c) 6 counter-rotating cells (mode $m=3$), (d) 4 counter-rotating cells (mode
$m=2$), and eventually a unique rolling cell (mode $m=1$). Drop keeps on
rolling but eventually stop. The consecutive snapshots are extracted from
movie SM1. (f-h) The horizontal median cut views of the most unstable mode
from the numerical stability analysis for some selected radii. The radii for
the inner flow symmetry transitions from experiments and from stability
analyses are plotted is (a) as black and orange lines, respectively.

Temporal variations in $R$ are best fitted by a linear law, plotted as the
blue dotted line and with equation:

$R(t)=R_{0}(1-t/\tau_{0}),$ (1)

where $R_{0}=3.7\pm 0.1$ mm and $\tau_{0}=176.5$ s. This time is much larger
than the characteristic time of the internal motion $R/V$, since tracers
inside the liquid and at the drop surface reveal flow velocities $V$ as high
as a few cm/s. The separation of time-scales $\tau_{0}\gg R/V\approx 0.02$ s
suggests that the evaporation-driven dynamics can be decoupled from the inner
dynamics. As a result, a Leidenfrost drop can be considered in quasi-static
equilibrium at any time, a key assumption to discuss the stability of
Leidenfrost inner flows. Moreover, as will be discussed in §IV.3.2, the
instability develops faster than $\tau_{0}$, which supports the quasi-static
stability analysis of the Leidenfrost drop.

###  II.2 Leidenfrost drop shapes

A consequence of the timescale separation is that a given volume of liquid
adopts the static shape of a non-wetting drop (Roman et al., 2001). If we
denote $C(z)$ as the local curvature at a given height $z$ and $C_{0}$ as the
curvature at the apex, the balance of hydrostatic and Laplace pressures can be
written $C(z)=C_{0}+z/\ell_{c}^{2}$. We introduce the curvilinear abscissa
$s$, the horizontal radius at a given height $r(z)$ and the angle $\beta$
tangent to the interface, the previous equation can be recast into :

$\frac{\sin\beta}{r}+\frac{\mathrm{d}\beta}{\mathrm{d}s}=C_{0}+\frac{z}{\ell_{c}^{2}}.$
(2)

A numerical integration of eq.(2) for $\beta$ ranging from 0 to $\pi$ and
$\mathcal{C}_{0}\ell_{c}$ ranging from 0.5 to 10, provides the drop shape for
radii ranging from $R=0.9$ mm to 4.3 mm, plotted as full lines in Figure 1(b).
These theoretical shapes are found to match the shapes obtained for water
drops kept in place by a needle (see dotted lines), except at the drop north
pole, where the needle locally forms a meniscus. Increasing the drop size
tends to saturate the puddle height $H$ at its maximum value $H_{\max}\approx
2\ell_{c}$, that is roughly 5 mm for water at 100${}^{\circ}\mathrm{C}$. Drops
smaller than the capillary length $\ell_{c}$ are quasi-spherical while drops
larger than $\ell_{c}$ get flattened owing to gravity. As a consequence,
evaporation induces geometric changes, particularly on the drop aspect ratio
$2R/H$. The presence of strong internal flows could in principle deform the
liquid interface. The Reynolds number associated to the inner flows with
typical velocity $V\sim$ cm/s and kinematic viscosity $\nu$ writes as
$Re=RV/\nu$. For a millimetric drop, we have $Re\approx 10^{3}$, so that
inertia overcomes viscosity. The Weber number $We$, which compares inertia to
the resisting capillarity is expressed as $We=\rho_{0}V^{2}R/\sigma_{0}\sim
10^{-2}$. Capillary therefore outbalances inertia, which justifies that the
drop shape does not deviate from the static ones.

###  II.3 Drop internal flow structure

The apparent quietness of Leidenfrost drops does not reflect what really
happens inside the liquid. Side-viewed PIV measurements performed in a median
plane of a water drop containing tracer particles have revealed strong
internal flows, with velocities of a few cm/s. As the drop shrinks owing to
evaporation, Leidenfrost flows undergo a series of successive symmetry
breakings. This is further evidenced by focusing on the drop top surface. A
Leidenfrost drop, kept on a concave substrate, is seeded with surface
particles that have a greater affinity for the air interface. These hollow
glass beads are i) pre-dispersed in water,ii) skimmed from the interface
(where particles accumulate owing to a wetting or shape defect), and iii)
introduced in a pre-dispensed Leidenfrost drop. Despite the apparent
axisymmetry of the experiment, interfacial flow structures emerge, as
illustrated by the top-views in Figures 1(c-e) and visible in the movie SM1).
When the drop has a radius $R>2.5$ mm, it hosts multiple vortices, which can
be described by the azimuthal wavenumber $m\geq 3$ in a cylindrical coordinate
system, by denoting ($r$, $z$, $\theta$) and using a periodic wavenumber
expansion in $\theta$ as $e^{im\theta}$, $m\in\mathbb{N}^{*}$. Within the
range $R=[1.8;2.5]$ mm, a mode $m=2$ clearly appears, while for smaller radius
$R<1.5$ mm, the droplet rolls in an asymmetric fashion, corresponding to a
mode $m=1$. Inner flows thus successively self-organize into 6 counter-
rotating cells ($m=3$, Figure 1(c)); 4 counter-rotating cells ($m=2$, Figure
1(d)), and eventually a unique rolling cell ($m=1$, Figure 1(e)). We can
notice that at some instance of movie SM1, the convective cells loose their
organization and coherence. The flow structures seems to be transiently
perturbed by the drop oscillation in the slightly curved well. However, the
dominant modes reappear within a second as a hint of the robustness of the
unstable modes. The transition from $m=2$ to $m=1$ can be also visualized by
side views, using PIV techniques, as in Bouillant et al. (2018b) (images
reproduced in Figure 6(a) and 5(a)) or even indirectly measured, as the drop
acceleration $a$ (extracted from top-views, for water drops initially at rest)
suddenly increases with the mode switching onto $m=1$. The experiment proposed
by (Bouillant et al., 2018b) is reproduced in the SI, for plate temperatures
ranging from $250$ to $450^{\circ}\mathrm{C}$. The accelerations $a$ of about
80 drops as a function of their radii $R$ exhibit similar jumps from $\sim$1
mm/s2 (in the regime where drops are flattened by gravity), up to 60 mm/s2 as
$R$ decreases to $\sim 1$ mm. This indeed corresponds to entering the self-
rotation and propelling regime (Bouillant et al., 2018a). The onset for the
transition from $m=2$ to $m=1$, referred as $R_{2\rightarrow 1}$, is found to
weakly depend on the plate temperature and to be within the interval $[1;1.5]$
mm (see §S6 of the SI).

###  II.4 Temperature gradient at the drop interface

To test the aforementioned thermally-based instability scenario, we need to
specify the thermal boundary conditions at the drop surface. The drop base is
expected to be maintained at roughly the water boiling point
($T_{b}=100^{\circ}\mathrm{C}$), while its apex is cooled down by the ambient
air.

Figure 2: (a) Infrared side views of an evaporating Leidenfrost water drop
deposited on a slightly curved surface of brass heated at 350∘C. Images, taken
with a thermal camera give access to the surface temperature (calibration
range from $-40^{\circ}$C to $+150^{\circ}$C to focus on water surface).
Images are extracted from movie SM2. (b) Surface temperature $T$ of a given
water drop along its central vertical axis $z$ (white dotted line), showing
the change of temperature profile as the drop radius $R$ decreases.

The temperature field in a Leidenfrost drop is measured using an infrared
camera (FLIR A600 series), calibrated on the temperature range
$[-40;+150]^{\circ}\mathrm{C}$, only suitable to see water, and not the brass
substrate. Water being opaque to infra-red wavelengths, this measurement
provides the "skin" temperature $T$. Figure 2 shows at $t=0$, that is when
$R=3.5$ mm, $T$ linearly decreases with height $z$ (in millimeters) as
$T=-3.75z+97.0^{\circ}\mathrm{C}$, reaching a maximum
$T_{\max}=97.0^{\circ}\mathrm{C}$ at the drop base ($z=0$). A temperature
difference $\Delta T\approx 25^{\circ}\mathrm{C}$ thus develops along the
drop. At $t=16$ s, the gradient is slightly smaller with
$T=-3.53z+97.0^{\circ}\mathrm{C}$, and it keeps decreasing as $R$ decreases.
These observations are confirmed by introducing a thermocouple inside the
liquid (see Figure S1 in the (SI)). For $t=74$ s, when $R=1.8$ mm, $T$
suddenly becomes homogeneous with $T\approx 88^{\circ}$C. As best visible in
the supplementary movie SM1, this coincides with the moment where the flow
switches to a symmetry $m=1$. At this transition, the drop starts to vibrate,
which enhances mixing. For $t>100$ s, $T$ thus becomes roughly homogeneous,
with a minimum at the drop center (along the rolling axis), and a maximum at
its periphery since fluid are periodically brought close to the hot plate. In
this rolling state, the liquid temperature is $T\approx 80^{\circ}$C, a value
smaller than the boiling point of water, as consequences of i) the
intensifying evaporation-driven cooling; ii) a reduction of the flattened area
at the drop base from which water is heated, which scales as $R^{4}/l_{c}^{2}$
(Mahadevan & Pomeau, 1999); iii) the temperature homogenization due to the
rolling-enhanced mixing. We now try to link the existence of such temperature
distributions to the origin and structure of the internal flows.

##  III Theory

###  III.1 Problem formulation

Based on the experimental observations, we develop a minimal model assuming
that i) the liquid adopts at any time the static shape of a non-wetting drop;
ii) the baseflow inside a drop is steady; iii) the temperature at the bottom
is fixed to the liquid boiling temperature; iv) the interaction between the
liquid and the surrounding gas is not solved completely but modelled by heat
transfer correlation laws applied on the side boundary. We sketch in Figure
3(a) the problem, where we represent from the side a static drop provided by
eq.(2). We denote by $\Omega$ the liquid domain, and by $\partial\Omega$ the
boundaries, which are decomposed into $\partial\Omega_{F}$, the upper free
interface and $\partial\Omega_{S}$, the bottom interface of the drop. We also
introduce $\partial\Omega_{A}$ the vertical centerline of the drop, the axis
of symmetry. Hence, the numerical computation is done only on the half domain
$r=[0;R]$ as illustrated in Figure 3(b). The bottom interface
$\partial\Omega_{S}$ is assumed to be isothermal at the temperature
$T_{s}=100^{\circ}\mathrm{C}$ (the phase change of a pure body occurs at a
given constant temperature), while the temperature on the side of the drop
$\partial\Omega_{F}$ needs to be evaluated using the heat transfer balance.

Figure 3:  Sketch of the problem. (a) The drops (domain $\Omega$) presents
boundaries $\partial\Omega$ including the upper free surface
$\partial\Omega_{F}$, the bottom interface $\partial\Omega_{S}$ and the axis
of symmetry $\partial\Omega_{A}$. (b) Thermal conditions.

###  III.2 Governing equations

Under the Boussinesq approximation, the governing Navier-Stokes equation for
the velocity fields $\mathbf{u}=[u_{r},u_{\theta},u_{z}]^{\mathrm{T}}$ and the
temperature $T$ in the cylindrical coordinate $(r,\theta,z)$ defined in the
domain $\Omega$ of boundaries $\partial\Omega$ reads,

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}$ $\displaystyle=-\frac{1}{\rho_{0}}\nabla
p-\frac{\rho(T)}{\rho_{0}}\mathbf{g}+\nu\nabla^{2}\mathbf{u}\quad$
$\displaystyle\mathrm{in}\ \Omega,$ (3a) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=0\quad$ $\displaystyle\mathrm{in}\ \Omega,$ (3b)
$\displaystyle\frac{\partial T}{\partial t}+\mathbf{u}\cdot\nabla T$
$\displaystyle=\kappa\nabla^{2}T\quad$ $\displaystyle\mathrm{in}\ \Omega,$
(3c)

where $\mathbf{g}$ is the gravitational acceleration in $z$ direction, $\nu$
the liquid kinematic viscosity and $\kappa$ its thermal diffusivity. The
boundary conditions for $\partial\Omega_{F}$, $\partial\Omega_{S}$ and
$\partial\Omega_{A}$ are respectively:

$\displaystyle\mathbf{S}\cdot\mathbf{n}$
$\displaystyle=\mathbf{S}_{\mathrm{gas}}\cdot\mathbf{n}-\sigma(\nabla\cdot\mathbf{n})\mathbf{n}+(\mathbf{I}-\mathbf{nn})\cdot\nabla\sigma\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{F},$ (4a)
$\displaystyle\mathbf{u\cdot n}$ $\displaystyle=0\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{F},$ (4b) $\displaystyle
k\mathbf{n}\cdot\nabla T$
$\displaystyle=-h_{c}(T-T_{a})-n^{\prime}\mathcal{L}\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{F},$ (4c) $\displaystyle T$
$\displaystyle=T_{0}\ $ $\displaystyle\mathrm{on}\ \partial\Omega_{S},$ (4d)

with suitable boundary conditions on the flow axis of symmetry,
$\partial\Omega_{A}$, detailed in §III.3 for linear perturbations and reading
${u}_{r}={u}_{\theta}=0$ for the base axisymmetric case. $\mathbf{S}$ and
$\mathbf{S}_{\mathrm{gas}}$ are the stress tensors in the liquid and the gas,
respectively, $\sigma$ the surface tension, $\mathbf{n}$ the normal vector,
$k$ the conductivity, $h_{c}$ the convective heat transfer coefficient in air,
$n^{\prime}$ the evaporation rate and $\mathcal{L}$ the latent heat. Equation
4c expresses the heat flux balance at the interface stemming from the liquid
($j_{k}$) transferred into air ($j_{c}$) and into evaporative cooling
($j_{e}$), as schematized in Figure 3(b).
Both the liquid density $\rho$ and the surface tension $\sigma$ are assumed to
vary linearly with the temperature: $\rho(T)=\rho_{0}(1-\beta_{w}(T-T_{0})),\
\sigma=\sigma_{0}-\sigma_{1}(T-T_{0}),$ where $\beta_{w}$ is the water thermal
expansion coefficient, $\beta_{w}={\rho_{0}}^{-1}{\partial\rho}/{\partial T}$
and $\sigma_{1}$ the surface tension variation with the temperature,
$\sigma_{1}=-{\partial\sigma}/{\partial T}$. We decompose now the temperature
$T$ as:

$T=T_{0}-\Theta(r,z),$ (5)

where $T_{0}$ is the water boiling temperature taken as reference
($T_{0}=100^{\circ}\mathrm{C}$) at $z=0$ and $\Theta$ the deviation from
$T_{0}$. The reference density $\rho_{0}$ and the surface tension $\sigma_{0}$
are defined at $T_{0}$. With these definitions, the density and the surface
tension write:

$\displaystyle\rho(T)=\rho_{0}\left(1-\beta_{w}\Theta\right),\qquad\sigma=\sigma_{0}-\sigma_{1}\Theta.$
(6)

Denoting $p_{0}$ the hydrostatic pressure, the pressure $p$ inside the liquid
decomposes as $p=p_{0}(z)+p_{1}(r,z)$. The time, velocity and length are
scaled with $H^{2}/\kappa,\kappa/H$ and $H$, respectively where $H$ is the
drop height. The temperature is scaled with the temperature difference $\Delta
T(=|\Theta_{z=0}-\Theta_{z=H}|)$ and the pressure is scaled with
$\rho_{0}\nu\kappa/H^{2}$, which naturally appears when plugging eq. 7 in eq.
3a, and comparing the viscous term to the pressure term,

$r=H\hat{r},\quad z=H\hat{z},\quad
t=\frac{H^{2}}{\kappa}\hat{t},\quad\mathbf{u}=\frac{\kappa}{H}\hat{\mathbf{u}},\quad{\Theta}=\Delta
T\hat{\Theta},\quad
p=\frac{\rho_{0}\nu\kappa}{H^{2}}\hat{p},\quad\nabla=\frac{1}{H}\hat{\nabla}.$
(7)

We provide in the SI a Table S1 with the physical parameters relative to
water, vapor, air and their interface relevant to describe our problem. We
introduce the dimensionless numbers, Prandtl, Rayleigh, Marangoni, Biot and
Sherwood numbers as:

${Pr}=\frac{\nu}{\kappa},\quad{Ra}=\frac{\beta_{w}g\Delta
TH^{3}}{\nu\kappa},\quad Ma=\frac{\sigma_{1}\Delta
TH}{\rho_{0}\nu\kappa},\quad Bi=\frac{h_{c}H}{k},\quad
Sh=\frac{h_{m}H}{D_{va}}.$ (8)

The governing equation (3) can be now recast into

$\displaystyle
Pr^{-1}\left(\frac{\partial\hat{\mathbf{u}}}{\partial\hat{t}}+\hat{\mathbf{u}}\cdot\hat{\nabla}\hat{\mathbf{u}}\right)$
$\displaystyle=-\hat{\nabla}\hat{p}_{1}+Ra\hat{\Theta}\mathbf{e}_{z}+\hat{\nabla}^{2}\hat{\mathbf{u}}\quad$
$\displaystyle\mathrm{in}\ \Omega,$ (9a)
$\displaystyle\hat{\nabla}\cdot\hat{\mathbf{u}}$ $\displaystyle=0\quad$
$\displaystyle\mathrm{in}\ \Omega,$ (9b)
$\displaystyle\frac{\partial\hat{\Theta}}{\partial
t}+\hat{\mathbf{u}}\cdot\hat{\nabla}\hat{\Theta}$
$\displaystyle=\hat{\nabla}^{2}\hat{\Theta}\quad$ $\displaystyle\mathrm{in}\
\Omega,$ (9c)

with the boundary conditions:

$\displaystyle-\hat{p}_{1}\mathbf{n}+(\hat{\nabla}\hat{\mathbf{u}}+(\hat{\nabla}\hat{\mathbf{u}})^{T})\cdot\mathbf{n}$
$\displaystyle=-Ma(\mathbf{I}-\mathbf{nn})\cdot\hat{\nabla}\hat{\Theta}\ $
$\displaystyle\mathrm{on}\ \partial\Omega_{F},$ (10a)
$\displaystyle\hat{\mathbf{u}}\cdot\mathbf{n}$ $\displaystyle=0\ $
$\displaystyle\mathrm{on}\ \partial\Omega_{F},\partial\Omega_{S},$ (10b)
$\displaystyle\mathbf{n}\cdot\hat{\nabla}\hat{\Theta}$
$\displaystyle=-Bi\left(\hat{\Theta}+\frac{T_{0}-{T}_{a}}{\Delta
T}\right)-\frac{n^{\prime}\mathcal{L}H}{k}\ $ $\displaystyle\mathrm{on}\
\partial\Omega_{F},$ (10c) $\displaystyle\hat{\Theta}$ $\displaystyle=0\ $
$\displaystyle\mathrm{on}\ \partial\Omega_{S}.$ (10d)

and suitable symmetry conditions on $\partial\Omega_{A}$ (see 3.3 for
details). The evaporation rate $n^{\prime}$ is linked to the $Sh$ number,
which represents the evaporative heat transfer coefficient. The heat exchanges
on $\partial\Omega_{F}$, which determine $Bi$ and $Sh$ are modelled using the
Ranz-Marshall correlation (Ranz & Marshall, 1952; Bergman et al., 2011) as
detailed in the SI (§S2), where the input ambient air temperature $T_{a}$ is
also measured and provided in SI (SI) (§S4).

###  III.3 Stability analysis

The steady toroidal baseflow ($\mathbf{u}_{b},p_{b},\Theta_{b}$) is obtained
by solving the nonlinear steady state solution of (3) satisfying the boundary
condition (4) using the Newton method. The baseflow is solved using the
dimensional equations since some dimensionless numbers depend on $\Delta T$,
which is also an unknown of the thermal problem.111Note that the non-
dimensional equation can be resolved by defining $Ra$ and $Ma$ with the known
parameters, i.e. $T_{s}$. However, we kept the classical definition of $Ra$
and $Ma$ which are with $\Delta T$. We first compute the baseflows, estimate
$\Delta T$ and deduce the dimensionless numbers. This approach differs from
thermoconvective studies on a flat plate or in a cylinder. In these
situations, there is no radial pressure gradient and the vertical gradient is
solely balanced with the density, without inducing any velocity field. The
temperature difference $\Delta T$ and thermal dimensionless parameters are
control parameters. In contrast, when the radial pressure gradient is nonzero
(as in the Leidenfrost configuration, owing to the boundary conditions), it
induces a steady flow with nonzero velocity. Both $\Delta T$ and the
dimensionless parameters are solutions of the problem. Assuming infinitesimal
perturbations on the baseflow with the complex frequency $\omega$ and the
azimuthal wavenumber $m$, the dimensionless flow field, pressure and
temperature are decomposed using the normal mode expansion:

$[\hat{u}_{r},\hat{u}_{\theta},\hat{u}_{z},\hat{p}_{1},\hat{\Theta}](r,\theta,z,t)=[{u}_{r},{u}_{\theta},{u}_{z},{p},{\Theta}](r,z)\exp(\mathrm{i}m\theta-\mathrm{i}\omega
t)+c.c.,$ (11)

where $c.c.$ indicates complex conjugate. The linearized equation (9) becomes

$\displaystyle
Pr^{-1}\left(\mathrm{i}\omega\mathbf{u}+\mathbf{u}_{b}\cdot\nabla_{m}+\mathbf{u}\cdot\nabla_{0}\mathbf{u}_{b}\right)$
$\displaystyle=-\nabla_{m}{p}+Ra\Theta\mathbf{e}_{z}+\nabla^{2}_{m}\mathbf{u}\quad$
$\displaystyle\mathrm{in}\ \Omega,$ (12a)
$\displaystyle{\nabla}_{m}\cdot\mathbf{u}$ $\displaystyle=0\quad$
$\displaystyle\mathrm{in}\ \Omega,$ (12b)
$\displaystyle\mathrm{i}\omega\Theta+\mathbf{u}_{b}\cdot{\nabla_{m}}\Theta+\mathbf{u}\cdot{\nabla_{0}}\Theta_{b}$
$\displaystyle={\nabla}_{m}^{2}\Theta\quad$ $\displaystyle\mathrm{in}\
\Omega,$ (12c)

where $\nabla_{m}$ represents the derivative in $\theta$ is replaced by
$\mathrm{i}m$. The boundary conditions on $\partial\Omega_{F}$ are

$\displaystyle-{p}\mathbf{n}+({\nabla}_{m}\mathbf{u}+({\nabla}_{m}\mathbf{u})^{T})\cdot\mathbf{n}$
$\displaystyle=-Ma(\mathbf{I}-\mathbf{nn})\cdot{\nabla}_{m}{\Theta}\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{F},$ (13a)
$\displaystyle\mathbf{u}\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{F},$ (13b)
$\displaystyle\mathbf{n}\cdot{\nabla}_{m}{\Theta}$
$\displaystyle=-Bi{\Theta}\quad$ $\displaystyle\mathrm{on}\
\partial\Omega_{F}.$ (13c)

Note that the perturbation temperature is now only affected by the $Bi$ as the
other heat transfer coefficients are constant and contribute only in the zero
order baseflow equation. One could also linearize the last term in (4c) and
include in the stability analysis, but the effect on the stability results is
negligible (Yim et al., 2021). The condition prescribed on the drop axis of
symmetry $\partial\Omega_{A}$ (illustrated in figure 3) depends on $m$, in
particular on the mode symmetry. The following conditions are thus used as in
Batchelor & Gill (1962):

$\displaystyle m=0,$ $\displaystyle\quad u_{r}=u_{\theta}=\frac{\partial
u_{z}}{\partial r}=\frac{\partial\Theta}{\partial r}=0\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{A},$ (14a) $\displaystyle m=1,$
$\displaystyle\quad u_{z}=\Theta=p=\frac{\partial u_{r}}{\partial
r}=0=\frac{\partial u_{\theta}}{\partial r}=0\quad$ $\displaystyle\mathrm{on}\
\partial\Omega_{A},$ (14b) $\displaystyle m\geq 2,$ $\displaystyle\quad
u_{r}=u_{\theta}=u_{z}=\Theta=p=0\quad$ $\displaystyle\mathrm{on}\
\partial\Omega_{A}.$ (14c)

Finally, we prescribe on $\partial\Omega_{S}$, a Dirichlet condition for the
temperature and a free slip condition for the velocity since the bottom
surface is not in contact with the plate:

$\displaystyle\mathbf{u}\cdot\mathbf{n}$ $\displaystyle=0\quad$
$\displaystyle\mathrm{on}\ \partial\Omega_{S},$ (15a) $\displaystyle\Theta$
$\displaystyle=0$ $\displaystyle\mathrm{on}\ \partial\Omega_{S}.$ (15b)

###  III.4 Numerical methods

All the numerical analyses are performed using FreeFEM++ software (Hecht,
2012) for axisymmetric cylindrical coordinates $(r,z)$. The velocity, pressure
and temperature are discretized with Taylor-Hood P2, P1 and P2 elements,
respectively. The typical number of triangles is $\sim 10^{4}$. The linear
equations and the eigenvalue problem are solved using UMFPACK library and
ARPACK shift-invert method, respectively. Starting with the initial radius
$R=3.5$ mm, the nonlinear solution of (3) is solved for the given water
properties. Once the solution for one radius is found, it is used as the
initial guess for the smaller radius. Within 5 iterations, the L2 norm
residual becomes smaller than $1\cdot 10^{-8}$. The zero normal velocity
condition on the free surface is applied using the Lagrange multiplier method
(Babuška, 1973; Yim et al., 2021).

##  IV Numerical results

###  IV.1 Pure buoyancy induced flow ($Ma=0$)

####  IV.1.1 Baseflow ($Ma=0$)

Figure 4:  Baseflows for $Ma=0$: (a) $R=2$ mm, (b) $R=1.3$ mm and (c) $R=0.8$
mm. The temperature and velocity fields are shown in color and with arrows.
(d) The spectrum of growth rate $\omega_{i}$ with $R$.

Let us first consider the case of pure buoyant flows, neglecting
thermocapillary effects. The baseflow is thus computed following (3) while
setting the superficial stress on $\partial\Omega_{F}$ to zero. We restrict
this parametric study in $R$ to the limit $R\lesssim 2$ mm, since according to
IR measurements (Figure 2), the drop surface temperature tends to homogenize,
suppressing Marangoni surfaces flows. Figure 4 shows the baseflow obtained in
this limit of $Ma=0$, for drop radii $R=2$ mm, $R=1.3$ mm and $R=0.8$ mm,
respectively. In the absence of surface tension gradient, the baseflow
exhibits pure thermal convection: the flow rises along the center axis, warmer
since it is insulated from the drop interface, and descends along the side of
the drop, where the drop is cooled. The inner velocities for purely buoyant
flows ($\sim 1$ cm/s) underestimate the experimental observation ($\sim 5$
cm/s) of the similar size of drop (see Figure 5a,c for the experimental
measurements).

####  IV.1.2 Stability analysis ($Ma=0$)

The stability of purely thermobuoyant flows is herein considered. Figure 4(d)
shows the growth rate $\omega_{i}$ (imaginary part of the complex frequency
$\omega$) as a function of $R$ for the modes $m=0,1,2,3$. A positive growth
rate indicates the grow of the perturbations leading to the instability. As
shown in Figure 4(d), the $m=1$ mode is only unstable mode for $R>0.6$ mm. The
frequency $\omega_{r}$ of this mode (not shown) is zero, corresponding to a
steady unstable mode. The corresponding Rayleigh number $Ra$ decreases from
$20000$ to 10 as $R$ decreases from 2 to 0.5 mm (see Fig. S8b in SI), reaching
the value $Ra\sim 630$ when $R\sim 0.6$ mm, for which flows get stable.
Interestingly, this limit also corresponds to the radius where the droplet
ability to self-rotate and propel disappears, as noticeable in the last stage
of movie SM1. Both the measured propelling acceleration and the droplet base
asymmetry vanish below $R\lesssim 0.6$ mm until it stops (the measurement in
Figure 1(a) then ceases). This suggests that thermobuoyant effects become
stable to non-axisymmetric disturbances and flows stabilize as the drop size
reduces below a critical value. Figure 5 compares the $m=1$ unstable mode to
the flow fields obtained in a Leidenfrost drop with $R\sim 0.9$ mm. The mode
$m=1$ with a structure describes well the solid-like rolling motion in the
experimental observation.

Figure 5:  (a,c) Velocity fields within a droplet with $R\sim 0.9$ mm deduced
from PIV measurements. (b, d) Corresponding flow fields deduced from the
numerical stability analysis in the absence of Marangoni effects ($m=1$,
$Ma=0$, $Ra=6.1\cdot 10^{3}$). Velocity arrows are plotted within the lateral
(a,b) and horizontal (c,d) planes. The red dashed line in (c) indicates the
area of the flattened base of the drop on which the bottom camera focuses.
Color in (b,d) indicates the perturbed temperature field (normalized with its
maximum value).

Although the $m=1$ unstable mode with $Ma=0$ represents well the rolling
motion of small drops observed in the experiments, it fails to describe the
presence of higher azimuthal modes for large drops. This prompts us to look at
the $Ma\neq 0$ case for larger radii.

###  IV.2 Reduced Marangoni approximation

A major challenge for computing the baseflow in a Leidenfrost drop is to be
able to predict the surface tension distribution. When we take the exact
surface tension temperature dependence $\sigma_{1}$ provided in the Table S1
(SI), the Marangoni flows are largely overestimated compared to experimental
observations (see Figure S3 of the SI). It is known however that Marangoni
effects are very often markedly reduced (Hu & Larson, 2005a, 2006;
Dhavaleswarapu et al., 2010) or even almost absent (Marin et al., 2011;
Gelderblom et al., 2012; Dash et al., 2014).Based on measured quantities, we
first try to correct the surface tension temperature dependence $\sigma_{1}$.
To that end, we evaluate the Rayleigh $Ra$ and Marangoni $Ma$ numbers, which
compare the buoyancy and surface tension stresses in comparison to inertia,
respectively, as defined in (8). Using the parameters documented in Table S1
of the SI, they culminate to $Ra\sim 1.8\cdot 10^{5}$ and $Ma\sim 3.6\cdot
10^{5}$ for $\Delta T=25^{\circ}\mathrm{C}$, $H=4$ mm, $R=3.5$ mm. Both values
greatly exceed expected critical values for the onset for Rayleigh-Bénard and
Marangoni instabilities (see Table S2 in the SI) – typically $Ra_{c}\sim
O(10^{3})$ and $Ma\sim O(10^{2})$, for similar geometries (Chandrasekhar,
1961). As detailed in §S7 of the SI, the effective surface tension variation
$\sigma_{1,\rm eff}$ (and thus $Ma_{\rm\rm eff}$) is determined using
numerical analysis. We select $\sigma_{1,\rm eff}$ that best describes the
temperature difference $\Delta T$ within the liquid (Fig 2), as well as the
flow velocities, typically 5 cm/s. The temperature profiles reported in Figure
2 (left panels, $R=3.5$ mm and $R=3.0$ mm) are best represented by the curve
with $\sigma_{1,\rm eff}=4\cdot 10^{-5}\sigma_{1}$ as shown in Figure S7 of
the SI for the surface temperature and Figure 6 for velocity field. The
surface stress seems to be reduced by a few order of magnitude. This has also
been reported in Savino et al. (2002); Hu & Larson (2002, 2005b, 2006), where
it is ascribed to surface contamination, or to the fact that at large
Marangoni surface stress, it gets moderated by dissipation or by transport-
limited properties of the fluid, reducing the achievable velocities (see Fig.6
of the SI). In the following stability analysis, we thus use this reduced
surface tension variation by adjusting $\sigma_{1}$ to $\sigma_{1,\rm
eff}=4\cdot 10^{-5}\sigma_{1}$, which is the only tunable parameter of our
study (all other physical parameters are kept exact values of given thermal
properties, provided in the Table S1).

###  IV.3 Thermocapillary flow

####  IV.3.1 Baseflow with effective Ma

The baseflow in Leidenfrost drops is now computed as in IV.1.1, adding to
thermobuoyant effects reduced thermocapillary effects ($\sigma_{1,\rm
eff}=4\cdot 10^{-5}\sigma_{1}$ and $Ma_{\rm eff}$) and shown in Figure 6
compared to the experimental measurement (Figure 6a). As the surface tension
gradient induces Marangoni flow from low surface tension to the higher one
along the surface boundary, the flow direction is opposed than purely buoyant
flow: it sinks on the center-line and rises along the surface. The typical
velocity magnitude is about $\sim 5$ m/s which is similar to the experimental
observation,which corroborates the choice of effective thermocapillary
gradient done by tuning the temperature distribution in the previous section
IV.2.

Figure 6:  (a) PIV measuments in a drop with $R=2.5$ mm (Bouillant et al.,
2018b). (b) Baseflow for $Ma_{\rm eff}=11.3$ ($\sigma_{1,\rm eff}=4\cdot
10^{-5}\sigma_{1}$). The colormap and arrows give the inner temperature and
velocity.

####  IV.3.2 Stability analysis with effective Ma

Figure 7:  Growth rates and the corresponding frequencies for (a,b) $m=1$,
(c,d) $m=2$ and (e,f) $m=3$ for the two least unstable modes as a function of
decreasing radius. The dominant modes for each $m$: (g) growth rate and (h)
corresponding frequency.

Figure 7 shows the dominant eigenvalues as a function of decreasing radius
(with $\sigma_{1,\rm eff}=4\cdot 10^{-5}\sigma_{1}$). Compared to the pure
buoyant flow, there exist several unstable azimuthal modes ranging from $m=1$
to 3. Figures 7(a,c,e) show the two least stable growth rates $\omega_{i}$ and
their corresponding frequencies $\omega_{r}$ are shown in Figures 7(b,d,f) for
azimuthal wavenumbers $m=1,2$ and 3. For all modes $m=1,2,3$, there exist two
branches of unstable modes: one with $\omega_{r}=0$ (steady) and other with
$\omega_{r}\neq 0$ (unsteady). The $m=1$ mode (Figure 7a,b) is unstable both
at large radius $R>2.8$ mm and small radius $R<1.3$ mm, but it shows a window
of stability for intermediate values of the radius. The prevailing unstable
mode for $R>2.8$ mm is unsteady, with frequency $\omega_{r}\sim 10\
\mathrm{s^{-1}}$, while for $R<1.3$ mm, it becomes steady.
Mode $m=2$ (Figure 7c,d) is only unstable in the intermediate radius range
$R=[2.3;1.3]$ mm. At radius $R\sim 2$ mm, the unsteady branch dominates, yet,
with a small frequency, but the steady branch takes over at smaller $R$. Mode
$m=3$ (Figure 7e,f) is unstable when $R<2.8$ mm. Its steady branch dominates
at large $R$ but it becomes unsteady as the radius gets closer to $R\approx
2.8$ mm. The unsteady branch of mode $m=3$ reaches a maximum growth rate
around $R=2.4$ mm.
Finally, Figures 7g,h collect and retain only the most unstable mode $m$ and
their corresponding frequencies. We note in Figure 7g that the axisymmetric
mode $m=0$ is always stable ($\omega_{i}<0$). The dependence of the dominant
azimuthal mode on radius thus becomes explicit: for $R>2.1$ mm, the mode $m=3$
is the most unstable (first steady and then oscillatory). As the radius
decreases further, the $m=3$ mode becomes stable around $R=2$ mm. At $R=2.2$
mm, the mode $m=2$ starts to grow and becomes the only unstable mode in the
range of radii $R=[2;1.3]$ mm, with the maximum growth rate at $R=1.6$ mm. For
$R<1.3$ mm, the steady mode $m=1$ takes over as the dominant unstable mode.
This trend is very similar to the apparent mode transition in experiments,
which is $m\geq 3$ for $R\gtrsim 2.8$ mm, $m=2$ for $R\approx[2.5;1.8]$ mm and
$m=1$ for $R\approx 1.5$ mm, as shown in Figure 1.

Figure 8:  Eigenmodes of the most unstable modes at radii (a,d) $R=3$ mm
($Ma_{\rm eff}=12,Ra=1.5\cdot 10^{5}$), (b,e) $R=1.5$ mm ($Ma_{\rm
eff}=9.6,Ra=5\cdot 10^{4}$) and (c,f) $R=1$ mm ($Ma_{\rm eff}=8,Ra=2.4\cdot
10^{4}$). Colors indicate temperature perturbations and the arrows show the
in-plane velocity perturbations normalized with their maximum real values. The
top view is a plane cut at the maximum radius $R=R_{\max}$.

A growth rate $\omega_{i}\sim 1~{}\mathrm{s^{-1}}$ implies that within 2.3 s,
the perturbation amplitude becomes 10 times larger than its initial value. For
all modes, we verify that the growth rate is much faster than the drop
evaporation rate $\omega_{i}\gg(1/R)\mathrm{d}R/\mathrm{d}t$, underlying the
quasi-steady assumption in §II.1.
Figure 8 illustrates the eigenvectors of the most unstable mode at radii
$R=3,\ 1.5$ and $1$ mm. The top panel displays side-views while the right one
shows top-views: $xy$ plane cut at the maximum radius $R=R_{\max}$. The colors
indicate temperature perturbations and the arrows are the in-plane velocity
perturbations. For $R=3$ mm (Figure 8ad), the mode $m=3$ dominates.
Temperature perturbations are localized along the interface and the central
axis of the drop while velocity perturbations are localized where the
temperature perturbations are the weakest. The top view shows three counter-
rotating vortex pairs, very similar to our observations (Figure 1c). For
$R=1.5$ mm (Figure 8be), the mode $m=2$ becomes the most unstable thus
dominant mode. Temperature perturbations are maximum near the central axis,
while velocity perturbations remain localized where temperature perturbations
are weak. The top view shows 4 vortex cells (or 2 vortex pairs), in close
agreement to the experiments (Figure 1d). For $R=1$ mm (Figure 8cf), the
displacement mode $m=1$ is the most unstable. Temperature perturbations are
null on the central axis, where velocity perturbations are maximum. However,
its structure does not match the solid-like rolling motion reported in
Bouillant et al. (2018b) and illustrated in the top of Figure 1e. We interpret
this at a consequence of the temperature homogenization evidenced at small $R$
in Figure 2. The temperature difference at the liquid surface indeed tends to
vanish as $R$ becomes smaller than 1.8 mm. We thus expect thermocapillary
effects to weaken and eventually vanish. We provide in Figure S8 (SI), the
(reduced) Marangoni and Rayleigh numbers with $R$, showing that both effects
weaken as the drop shrinks in size, with a more abrupt decrease of
thermocapillary effects than thermobuoyant effects, prompting us to use our
predictions for the case $Ma=0$ when $R$ becomes sufficiently small. The
observed rolling motion is thus better captured by the eigenmode for $Ma=0$ as
shown in Figure 5 than the one with $Ma\sim O(10)$ in Figure 8e. Moreover, the
drop shape is, as yet, steady, but the internal flows could be coupled to the
drop envelop deformation. This extension of our model would enable the
exploration of the limit of even larger drops ($R\gtrsim 5$ mm), where star-
pulsations appear. Finally, it would be interesting to see how our model
applies to the inverse-Leidenfrost situation (Gauthier et al., 2019), for
which a rolling mode $m=1$ seems dominant in the millimetric, quasi-spherical
droplet. Their configuration is however essentially different since
temperature gradients are reversed, the drop shape remains quasi-spherical on
the deformable bath, and, as soon as the drop freezes, both the thermo-buoyant
and thermo-capillary flows should extinguish.

## References

  * SI (SI) SI See Supplemental Material at [URL will be inserted by publisher] for supplementary experiments and details on the numerical model inputs. .
  * Babuška (1973) Babuška, Ivo 1973 The finite element method with lagrangian multipliers. Numerische Mathematik 20 (3), 179–192.
  * Batchelor & Gill (1962) Batchelor, G. K. & Gill, A. E. 1962 Analysis of the stability of axisymmetric jets. J. Fluid Mech. 14 (4), 529–551.
  * Bergman et al. (2011) Bergman, T.L., Incropera, F.P., DeWitt, D.P. & Lavine, A.S. 2011 Fundamentals of Heat and Mass Transfer. Wiley.
  * Bouillant et al. (2021a) Bouillant, A., Cohen, C., Clanet, C. & Quéré, D. 2021a Self-excitation of leidenfrost drops and consequences on their stability. PNAS 118 (26).
  * Bouillant et al. (2021b) Bouillant, A., Lafoux, B., Clanet, C. & Quéré, D. 2021b Thermophobic leidenfrost. Soft Matter pp. –.
  * Bouillant et al. (2018a) Bouillant, A., Mouterde, T., Bourrianne, P., Clanet, C. & Quéré, D. 2018a Symmetry breaking in leidenfrost flows. Phys. Rev. Fluids 3 (10), 100502\.
  * Bouillant et al. (2018b) Bouillant, A., Mouterde, T., Bourrianne, P., Lagarde, A., Clanet, C. & Quéré, D. 2018b Leidenfrost wheels. Nat. Phys. 14 (12), 1188–1192.
  * Brandão & Schnitzer (2020) Brandão, Rodolfo & Schnitzer, Ory 2020 Spontaneous dynamics of two-dimensional leidenfrost wheels. Physical Review Fluids 5 (9), 091601.
  * Brunet & Snoeijer (2011) Brunet, P & Snoeijer, J H 2011 Star-drops formed by periodic excitation and on an air cushion - A short review. European Physical Journal: Special Topics 192 (1), 207–226.
  * Burton et al. (2012) Burton, J. C., Sharpe, A. L., van der Veen, R. C. A., Franco, A. & Nagel, S. R. 2012 Geometry of the Vapor Layer Under a Leidenfrost Drop. Physical Review Letters 109 (7), 074301.
  * Celestini et al. (2012) Celestini, Franck, Frisch, Thomas & Pomeau, Yves 2012 Take off of small leidenfrost droplets. Physical Review Letters 109 (3), 1–5, arXiv: 1206.5932.
  * Chandrasekhar (1961) Chandrasekhar, Subrahmanyan 1961 Hydrodynamic and hydromagnetic stability. Courier Corporation.
  * Cousins et al. (2012) Cousins, Thomas R, Goldstein, Raymond E, Jaworski, Justin W & Pesci, Adriana I 2012 A ratchet trap for Leidenfrost drops. Journal of Fluid Mechanics 696, 215–227.
  * Dash et al. (2014) Dash, S., Chandramohan, A. a, Weibel, J. A. & Garimella, S. V. 2014 Buoyancy-induced on-the-spot mixing in droplets evaporating on nonwetting surfaces. Phys. Rev. E 90 (6), 062407.
  * Dhavaleswarapu et al. (2010) Dhavaleswarapu, H. K., Migliaccio, C. P., Garimella, S. V. & Murthy, J. Y. 2010 Experimental investigation of evaporation from low-contact-angle sessile droplets. Langmuir 26 (2), 880–888.
  * Dupeux et al. (2013) Dupeux, Guillaume, Baier, Tobias, Bacot, Vincent, Hardt, Steffen, Clanet, Christophe & Quéré, David 2013 Self-propelling uneven Leidenfrost solids. Physics of Fluids 25 (5), 1–7.
  * Garmett (1878) Garmett, WM. 1878 Leidenfrost’s Phenomenon. Nature 17 (441), 466–466.
  * Gauthier et al. (2019) Gauthier, Ana\̈mathrm{i}s, Diddens, Christian, Proville, Rémi, Lohse, Detlef & van der Meer, Devaraj 2019 Self-propulsion of inverse leidenfrost drops on a cryogenic bath. Proceedings of the National Academy of Sciences 116 (4), 1174–1179.
  * Gelderblom et al. (2012) Gelderblom, H., Bloemen, O. & Snoeijer, J. H. 2012 Stokes flow near the contact line of an evaporating drop. J. Fluid Mech. 709, 69–84.
  * Hecht (2012) Hecht, F. 2012 New development in freefem++. J. Numer. Math. 20 (3-4), 251–265.
  * Holter & Glasscock (1952) Holter, Norman J. & Glasscock, Wilford R. 1952 Vibrations of Evaporating Liquid Drops. Journal of the Acoustical Society of America 24 (6), 682–686.
  * Hu & Larson (2002) Hu, H. & Larson, R. G. 2002 Evaporation of a sessile droplet on a substrate. J. Phys. Chem. B 106 (6), 1334–1344.
  * Hu & Larson (2005a) Hu, H. & Larson, R. G. 2005a Analysis of the effects of marangoni stresses on the microflow in an evaporating sessile droplet. Langmuir 21 (9), 3972–3980.
  * Hu & Larson (2005b) Hu, H. & Larson, R. G. 2005b Analysis of the microfluid flow in an evaporating sessile droplet. Langmuir 21 (9), 3963–3971.
  * Hu & Larson (2006) Hu, H. & Larson, R. G. 2006 Marangoni effect reverses coffee-ring depositions. J. Phys. Chem. B 110 (14), 7090–7094.
  * Johnson & Narayanan (1999) Johnson, D. & Narayanan, R. 1999 A tutorial on the rayleigh–marangoni–bénard problem with multiple layers and side wall effects. Chaos 9 (1), 124–140.
  * Kuhlmann & Rath (1993) Kuhlmann, H. C. & Rath, H. J. 1993 Hydrodynamic instabilities in cylindrical thermocapillary liquid bridges. J. Fluid Mech. 247, 247–274.
  * Leidenfrost (1756) Leidenfrost, Johann Gottlob 1756 De Aquae Communis Nonnullis Qualitatibus Tractatus. Ovenius, Duisburg (1756); transl. Wares, C. On the fixation of water in diverse fire. International Journal of Heat and Mass Transfer 9, 1153–1166 (1966).
  * Linke et al. (2006) Linke, H., Aleman, B. J., Melling, L. D., Taormina, M. J., Francis, M. J., Dow-Hygelund, C. C., Narayanan, V., Taylor, R. P. & Stout, A. 2006 Self-propelled leidenfrost droplets. Physical Review Letters 96 (15), 2–5.
  * Mahadevan & Pomeau (1999) Mahadevan, L. & Pomeau, Y. 1999 Rolling droplets. Phys. Fluids 11 (9), 2449–2453.
  * Marin et al. (2011) Marin, A. G., Gelderblom, H., Lohse, D. & Snoeijer, J. H. 2011 Order-to-disorder transition in ring-shaped colloidal stains. Phys. Rev. Lett. 107 (8), 085502.
  * Nield (1964) Nield, D. A. 1964 Surface tension and buoyancy effects in cellular convection. J. Fluid Mech. 19 (3), 341–352.
  * Pearson (1958) Pearson, J. R. A. 1958 On convection cells induced by surface tension. J. Fluid Mech. 4 (5), 489–500.
  * Pomeau et al. (2012) Pomeau, Yves, Le Berre, Martine, Celestini, Franck & Frisch, Thomas 2012 The Leidenfrost effect: From quasi-spherical droplets to puddles. Comptes Rendus Mécanique 340 (11-12), 867–881.
  * Quéré (2013) Quéré, David 2013 Leidenfrost Dynamics. Annual Review of Fluid Mechanics 45 (1), 197–215.
  * Ranz & Marshall (1952) Ranz, W. E. & Marshall, W. R. 1952 Evaporation from drops. Chem. Eng. Progr 48 (3), 141–146.
  * Raufaste et al. (2016) Raufaste, C., Bouret, Y. & Celestini, F. 2016 Reactive Leidenfrost droplets. Europhysics Letters 114 (4), 46005.
  * Roman et al. (2001) Roman, B., Gay, C. & Clanet, C. 2001 Pendulum, drops and rods: a physical analogy. J. Fluid Mech.(submitted) .
  * Rosenblat et al. (1982) Rosenblat, S., Davis, S. H. & Homsy, G. M. 1982 Nonlinear marangoni convection in bounded layers. part 1. circular cylindrical containers. J. Fluid Mech. 120, 91–122.
  * Savino et al. (2002) Savino, R., Paterna, D. & Favaloro, N. 2002 Buoyancy and marangoni effects in an evaporating drop. Journal of Thermophysics and Heat Transfer 16 (4), 562–574.
  * Scriven & Sternling (1960) Scriven, L. E. & Sternling, C. V. 1960 The marangoni effects. Nature 187 (4733), 186–188.
  * Sobac et al. (2014) Sobac, B., Rednikov, A., Dorbolo, S. & Colinet, P. 2014 Leidenfrost effect: Accurate drop shape modeling and refined scaling laws. Physical Review E 90 (5), 053011.
  * Sobac et al. (2017) Sobac, B, Rednikov, A, Dorbolo, S & Colinet, P 2017 Self-propelled Leidenfrost drops on a thermal gradient: A theoretical study. Physics of Fluids 29.
  * Tam et al. (2009) Tam, D., von Arnim, V., McKinley, G. H. & Hosoi, A. E. 2009 Marangoni convection in droplets on superhydrophobic surfaces. J. Fluid Mech. 624, 101–123.
  * Vrentas et al. (1981) Vrentas, J. S., Narayanan, R. & Agrawal, S. S. 1981 Free surface convection in a bounded cylindrical geometry. Int. J. Heat Mass Transf. 24 (9), 1513–1529.
  * Waitukaitis et al. (2017) Waitukaitis, Scott R., Zuiderwijk, Antal, Souslov, Anton, Coulais, Corentin & Van Hecke, Martin 2017 Coupling the Leidenfrost effect and elastic deformations to power sustained bouncing. Nature Physics 13 (11), 1095–1099.
  * Yim et al. (2021) Yim, E., Bouillant, A. & Gallaire, F. 2021 Buoyancy-driven convection of droplets on hot nonwetting surfaces. Phys. Rev. E 103, 053105.


# The Dynamics of Quantum Fluids

Henri Godfrin and Eckhard Krotscheck

###### Abstract

We review experimental and theoretical progress in the physical description of
the dynamics of the quantum fluids 4He and 3He. Historically, the elementary
excitations in these systems have been identified as phonons and rotons and,
in 3He, collective zero sound and spin-fluctuations. Both recent high-
precision measurements and theoretical methods have shown that the dynamics of
these systems is actually very rich as will be discussed in detail in the body
of this contribution.

Keywords— Quantum Fluids, Liquid Helium, Fermions, Bosons, Neutron Scattering,
X-Ray scattering, Elementary Excitations, Phonons, Rotons, Ripplons

##  1 The Helium Liquids

Helium (4He and the less common isotope 3He) and hydrogen are simple condensed
matter systems, amenable to fundamental quantum mechanical descriptions.
Helium remains liquid even at the absolute zero of temperature, an effect
explained by Heisenberg’s uncertainty principle, which has the consequence
that light atoms confined in a small volume have a large kinetic energy, as is
the case in a liquid or a solid, where each atom is confined by several
surrounding atoms. For this reason, in helium, kinetic energies are larger
than the weak interatomic potential energy. Even at the lowest temperatures,
substantial atomic motion is present, and the system remains liquid, unless
pressures on the order of 25 MPa are applied, and solidification is achieved.
Another quantum effect plays an important role in the behavior of liquid
helium: as the temperature is reduced below a few degrees Kelvin, the thermal
de Broglie wave-length of the atoms becomes larger than the interatomic
distance. The waves describing the helium atoms overlap significantly, and a
description in terms of individual wave-packets of distinguishable particles
becomes inappropriate, it must be replaced by a quantum many-body description
in terms of indistinguishable particles, governed by quantum statistics.

Quantum Mechanics and Statistical Physics [1, 2] explain well the properties
of Quantum Gases, i.e. ensembles of non-interacting particles, which
necessarily belong to one of the two possible classes: Bosons having an
integer spin quantum number s), or fermions with half-integer spin. 4He atoms
($s$=0) are bosons, and 3He atoms ($s$=$\frac{1}{2}$) are fermions.

The wave-functions of particles confined in a volume $V$ are quantized plane
waves characterized by their wave-vector ${\bf k}$. As shown by A. Einstein,
the ground state of a Bose gas is the BEC (Bose-Einstein condensate), where
all particles occupy the same microscopic state with ${\bf k}$=0. The ground
state of a Fermi gas is very different, because Pauli’s principle forbids
double occupancy for Fermions: the particles form a Fermi sphere in wave-
vector space, by completing successive energy levels of increasing wave-
vector, up to the Fermi wave-vector $k_{\rm F}$. The energy $E_{\rm F}$ of the
last one, called the Fermi level, can be considerable.

The excitation spectrum $\epsilon(k)$ of these quantum gases is simply given
by the kinetic energies of the individual particles,
$\epsilon(k)=\hbar^{2}k^{2}/2m$, where $m$ is the particle mass and $k$ the
the wave number. For bosons, at low temperatures, particles are excited from
the common ground state. The energy spectrum is the usual parabolic spectrum
of free particles described above. This behavior, however, is profoundly
modified by interparticle interactions, i.e., as is the case in Quantum Fluids
[1, 2, 3, 4].

The Helium quantum fluids are unique in the sense that they are very dense and
very quantum, making them a very challenging problem for theory. The issue is
made quantitative by looking at the simple Lennard-Jones model of the helium
fluids. In this model, the helium atoms interact via the interaction

$V(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]\equiv\epsilon
v\left(\frac{r}{\sigma}\right)$

where $\epsilon=10.22\,$K is the well depth of the potential and
$\sigma=2.556\,$Å is the core size [5].

The interaction is not the very best (quantitative calculations are in fact
made using the Aziz potential [6]), but it reproduces the ground state
properties of the liquids within a few percent. Measuring all energies in
units of the well depth, and all lengths in units of the core size, $x\equiv
r/\sigma$, the Hamiltonian of the system is

$\displaystyle\frac{1}{\epsilon}H({\bf x}_{1}\ldots,{\bf x}_{N})=$ (1)
$\displaystyle-\frac{\Lambda^{2}}{8\pi^{2}}\sum_{i}\nabla_{{\bf
x}_{i}}^{2}+\sum_{i<j}v(|{\bf x}_{i}-{\bf x}_{j}|)$

where

$\Lambda=\left(\frac{h}{m\epsilon\sigma^{2}}\right)^{\frac{1}{2}}$

is de Boer’s “quantum parameter” [7] $\Lambda$, giving the ratio between the
de Broglie wavelength and a typical diameter of the molecule. $\Lambda\approx
2.67$ for 4He, $\Lambda\approx 3.09$ for 3He, but $\Lambda\approx 1.3-1.7$ for
H2, HD, D2, and $\Lambda<0.1$ for heavy rare gases.

Moreover, the density of 4He at zero pressure is
$\rho_{0}=0.02185\,$Å${}^{-3}=0.365/\sigma^{-3}$ and for 3He
$\rho_{0}=0.0163\,$Å${}^{-3}=0.273/\sigma^{-3}$, in other words the average
particle distance is of the order of the core size.

Figure 1: Normalized Lennard-Jones “6-12” interatomic potential, as a function
of the normalized interatomic distance. The normalization factors for
distances ($\sigma$) and potential energies ($\epsilon$) are given in the
text. We also show, for comparison, the Aziz potential [6] which is the
currently accepted most accurate static 2-body interaction.

##  2 4He in 3D

Liquid helium (4He) is the archetype of a Bose Liquid. Immediately after the
discovery of the superfluidity of helium [8, 9] by Kapitza (Moscow) and Allen
and Misener (Toronto), F. London [10] understood that the system was
undergoing Bose-Einstein condensation, and L. Tisza proposed a “two-fluids
model” featuring a coexistence of a normal fluid and a superfluid [4]. A Bose-
Einstein condensed gas is not superfluid, however, and therefore the
superfluidity observed in helium is an additional effect, associated to the
interactions. A major step in the understanding of the physics of helium was
made by L.D. Landau, who found a deeper interpretation of Tisza’s two-fluids
model. Landau postulated that the “normal fraction”, which carries the
entropy, is in fact the ensemble of thermal excitations of the superfluid
ground state. This interpretation removed the unphysical consideration of two
categories within a set of indistinguishable atoms. Landau considered
elementary excitations of the fluid, in the form of quantized density
fluctuations (“phonons”), similar to the Debye phonon modes of a crystalline
solid. In order to calculate the thermodynamic properties by a Debye-like
approach, where the energy is obtained by a sum over all modes, weighted by
the Bose factor, Landau postulated that the dispersion relation $\epsilon(k)$
had a linear part $\epsilon=\hbar ck$, where $c$ is the speed of sound, and
another part displaying a minimum, the “roton gap” of energy $\Delta$.
Landau’s intuition was guided by specific heat measurements, showing a cubic
temperature dependence at low temperatures (below 0.5 K), as expected from the
linear part of the dispersion relation, as well a an exponential contribution
at higher temperatures (on the order of 1K), characteristic of an energy gap.

In a first article [11], Landau proposed a spectrum with two branches, with a
parabolic branch centered at k=0 in addition to the linear branch. This
spectrum was not found satisfactory by Landau himself, the presence of two
branches being reminiscent of Tisza’s “two kinds of atoms”.

In a second publication [12], he presented his celebrated “Landau dispersion
relation”, where a single branch of excitations characterizes the full
spectrum (Fig. 2). With two parameters, the sound velocity $c$ and the roton
gap $\Delta$, both obtained from the temperature dependence of the specific
heat, Landau’s theory provides an excellent account of the experimental
properties of superfluid helium at low temperatures. The existence of a
“second sound” collective mode, where the superfluid and normal components
oscillate with opposite phases while keeping the total density constant, as
proposed by Tisza, corresponds in Landau’s model to propagating entropy (and
hence temperature) waves.

Figure 2: The dispersion relation proposed by L. D. Landau in 1947 to describe
the elementary excitations in superfluid 4He at zero pressure.

Another fundamental observation due to Landau, is that the particular shape of
the dispersion relation leads to the superfluid behavior. A simple argument of
energy and momentum conservation shows that the creation of excitations by the
displacement of an object immersed in the fluid is not possible, unless the
object’s velocity exceeds a critical value. The linearity of the dispersion
relation at low wave vectors is essential for this to happen: in a Bose gas,
where the dispersion relations is parabolic, the critical velocity is zero,
the system remains normal even at the absolute zero of temperature. The
evolution of the dispersion relation from the parabolic shape characteristic
of a gas, to the linear behavior observed as the interaction is switched on,
has been explained by Bogoliubov (1947) [13] in the weak interaction limit.
The dispersion relation of the elementary excitations plays therefore a major
role in the determination of the thermal properties, and also in the
superfluid nature of the liquid’s ground state.

Landau’s phenomenological theory introduced the fruitful concept of
“quasiparticles”, describing in a simple way the elementary excitations of
complex many-body systems. The concept was of course already present in the
Debye theory of solids, but the generalization to helium, a non-periodic
system of strongly interacting particles, opened a new field of physics.

There were however two serious interrogations after Landau’s theory was
released. The first was obviously related to the necessity of an experimental
observation of the proposed dispersion curve. And the second, to the need of a
microscopic theory providing some insight in the superfluid helium behavior,
and able to test the validity of Landau’s model.

Experimentally, the evidence for a Landau-like dispersion relation was
missing: only indirect evidence obtained form thermodynamic measurement
(specific heat, thermal conductivity, viscosity, fountain pressure, etc.) was
available. Feynman and Cohen suggested the use of inelastic neutron
scattering, a technique already used to measure phonon dispersion curves in
crystals, to try to observe the excitation dispersion curve of superfluid
helium. The first proof of the existence of a roton minimum was reported by H.
Palevsky, K. Otnes, K. E. Larsson, R. Pauli and R. Stedman (1957) [16, 17],
clearly showing that the elementary excitations in superfluid helium are
extremely sharp, compared to those of normal helium. These exciting results
motivated extensive investigations of the dynamics of helium, covering the low
energy region, i.e. the Landau-like dispersion relation, the multi-excitations
region at higher energies, as well as the very high energy region of the
spectrum, where BEC can be investigated [18].

Figure 3: The time-of-flight neutron scattering spectrometer IN5 at the
Institut Laue-Langevin.

The principle of the inelastic neutron scattering is simple, but measurements
are demanding, and high flux neutron sources (reactors, spallation sources)
are needed. The helium sample is placed in an incident beam of “monochromatic”
(of a given energy) neutrons, which are scattered by the helium sample, and
then detected at some distance from the sample, as a function of angle and
time of arrival, in the so-called time-of-flight technique. Other methods
(triple-axis, back-scattering, spin-echo, etc.) have also been used, offering
a different coverage of the momentum-energy plane, and also different
resolutions in energy and momentum.

Energy and momentum conservation allow the determination of the energy
transfer and the momentum transfer from the neutron to the helium sample,
which correspond to the energy $\epsilon=\hbar\,\omega$ and momentum
$\hbar\vec{k}=\hbar\vec{Q}$ of the density fluctuations created in the helium.
The incident neutron has an energy $E_{i}$ and a wave-vector $\vec{k_{i}}$,
and, after scattering, a final energy $E_{f}$ and a wave-vector $\vec{k_{f}}$;
the wave-vector transfer is $\vec{Q}=\vec{k_{i}}-\vec{k_{f}}$, and the energy
transfer $\hbar\,\omega=E_{i}-E_{f}$. The first microscopic theory of
superfluid helium, due to R. P. Feynman (1954) [14], is a wonderful example of
physical intuition. He found a general form for a variational wave-function
which successfully explained the general shape of the dispersion relation,
including the roton minimum. It was soon completed by Feynman and Cohen [15],
who introduced, improving the variational wave-functions, the concept of
“backflow”: the flow of helium atoms around a moving helium atom.

Dealing with many-body correlations, including dynamic effects, has been
central to the development of microscopic theories, which have been
considerably refined in the last decades, reaching now a quantitative power of
prediction, as shown below.

Figure 4: Neutron scattering: an incident neutron of wavevector $\vec{k_{i}}$
creates an excitation of momentum $\hbar\vec{k}=\hbar\vec{Q}$ in the
superfluid. Q is the wavevector transfer (see text)

The theory of neutron scattering [19, 20, 18] relates the double differential
scattering cross section per target atom, which is the parameter measured by a
scattering experiment, to the dynamic structure factor $S(Q,\omega)$ and the
dynamic susceptibility $\chi(Q,\omega)$, the latter being the functions
calculated by microscopic theories. They are given by the following
expression, where $b_{c}$ is the coherent scattering length of helium:

$\frac{\partial^{2}\sigma}{\partial\Omega~{}\partial{E_{f}}}=\frac{b^{2}_{c}}{\hbar}\frac{k_{f}}{k_{i}}S(Q,\omega)$
(2)

Results are generally presented in terms of $S(Q,\omega)$, where sharp peaks
of high intensity are observed, as expected from the creation of Landau
“single-excitations”, as well as a more complex, broad landscape of multi-
excitations. A vast literature on the subject exists, as described in the
review articles by Woods and Cowley, Stirling, Glyde (see [18]), the book by
Glyde [20], and recent works [21, 22], which benefited from the substantial
progress of neutron scattering facilities and instruments.

The dispersion relation $\epsilon(k)$ of the sharp single-excitations, as
measured in recent inelastic neutron scattering experiments, is shown in Fig.
5. It closely resembles the curve predicted by Landau, with a linear part at
low wave-vectors (“phonons”), and a deep gap (“roton” minimum) at a wave-
vector $k$$\sim$2Å-1. The name roton is historical, no special rotation takes
place; the minimum in energy reflects an incipient localization (short-range
order), which eventually will lead to a liquid-solid transition under a
pressure $P$$\sim$2.5 MPa. The dispersion curve becomes flat for $k$$\geq$2.8
Å-1 as the energy reaches twice the roton gap, forming “Pitaevskii’s plateau”.

---
---

Figure 5: The dispersion relation $\epsilon(k)$ of excitations in superfluid
4He at zero pressure. Top: Experimental data, displaying the phonon, roton,
maxon, and plateau regions of the spectrum. The red dots are recent high-
accuracy, highly pixelized data (seen at this scale as a continuous line, the
dispersion is only visible at the highest wave-vectors) from Ref. [22], the
blue dots with error bars are from the original work by Cowley and Woods [23]
and the green dots from the high-momentum measurements by Glyde et al. [24].
Bottom: theoretical results, from Feynman’s model [14] to modern calculations
such as Brillouin-Wigner perturbation theory with correlated wave function
(BW, Ref. [25, 26]), diffusion (DMC [27, 28]) and path-integral (PIMC, [29])
Monte Carlo calculations and dynamic many-body theory (DMBT) [30, 21],
compared to the experimental curve [22] .

The calculation of the low temperature ($T<1.3\,K$) thermodynamic properties
of helium at low temperatures only involves the statistical physics of phonons
and rotons of low energy. That is, the properties of the strongly interacting
bosonic system are simply given by those of a collection of non-interacting
bosonic “quasi-particles” described by the phonon-roton curve [22]. At high
temperatures, where the number of rotons becomes very large, roton-roton
interactions start playing a role. At the “lambda point” (T=2.17 K), 4He
becomes a normal fluid.

Theoretical studies of the dynamic structure function in 4He began with the
work of Feynman [14] and Feynman and Cohen [15], the Feynman theory of
elementary excitations was developed in a systematic Brillouin-Wigner
perturbation theory by Jackson and Feenberg [31, 25, 32, 33] in terms of a
basis of Feynman excitation states

$\left\\{\left|{\bf k}\right\rangle\right\\}=\left\\{S^{-1}(k)\hat{\rho}_{{\bf
k}}\left|\Psi_{0}\right\rangle\right\\}$ (3)

where $\hat{\rho}_{\bf k}$ is the density operator and
$\left|\Psi_{0}\right\rangle$ is the ground state. An important contribution
was the identification of classes of theories for the dynamic structure
function [34] that satisfy the $\omega^{0}$ and $\omega^{1}$ sum rules
exactly. The most complete evaluation of the phonon–roton dispersion relation
in terms of Brillouin–Wigner perturbation theory was carried out by Lee and
Lee [26] who obtained an impressive agreement with the experimental
phonon–roton spectrum up the wave number of $Q=2.5$Å-1. The major drawback
with these old calculations was that the required input – pair- and three-body
distribution functions – were poorly known. Manousakis and Pandharipande [35]
tried to generalize the input states of the Brillouin-Wigner perturbation
theory to include “backflow” correlations in the spirit of Feynman and Cohen
[15]. Through the gradient operator acting on the wave function, in principle
dynamic correlations are introduced to all orders. The “backflow–function” is,
however, chosen per physical intuition rather than by fundamental principles,
and the evaluation of the perturbative series becomes very complicated.
Topologically, diagrams similar to those of Lee and Lee [26] were included.
While the accuracy of the theoretical roton energy is comparable to that of
Lee and Lee, one can clearly see an inconsistency since the energy of the
Pitaevskii-plateau [36] lies below twice the energy of the roton gap. More
recent progress [37, 30] used a hybrid approach of Brillouin–Wigner
perturbation theory and equations of motion for time–dependent multiparticle
correlation functions to derive a self- consistent theory of the dynamic
density–density response of 4He.

In a nutshell, the Feynman theory of excitations is generalized by writing the
dynamic wave function as

$\displaystyle\Psi({\bf r}_{1},\ldots,{\bf r}_{N};t)$ $\displaystyle=$
$\displaystyle\exp(U(t))\Psi_{0}({\bf r}_{1},\ldots,{\bf r}_{N})$ (4a)
$\displaystyle U(t)$ $\displaystyle=$ $\displaystyle\sum_{i}\delta u_{1}({\bf
r}_{i};t)+\sum_{i<j}\delta u_{2}({\bf r}_{i},{\bf r}_{j};t)\ +\ldots\,.$ (4b)

The Feynman form of the wave function is obtained by omitting all $n$-body
fluctuations $\delta u_{n}({\bf r}_{1},\ldots{\bf r}_{n};t)$ for $n\geq 2$.
The different theoretical descrptions basically differ by the way the $n$-body
fluctuations are determined.

Recent novel numerical methods [28, 27, 38, 39, 40] give access to dynamic
properties of quantum fluids. These are important algorithmic developments
that may ultimately aid in the demanding elimination of background and
multiple-scattering events from the raw data. Of course, it is generally
agreed upon that the model of static pair potentials like the Aziz interaction
[6] describes the helium liquids accurately. Hence, given sufficiently
elaborate algorithms and sufficient computing power, such calculations must
reproduce the experimental data. The aim of the works cited above is
different: The identification of physical effects like phonon-phonon, phonon-
roton, roton-roton, maxon-roton …couplings that lead to observable features in
the dynamic structure function is, from simulation data, only possible
a-posteriori whereas the semi-analytic methods permit a direct identification
of these effects, their physical mechanisms, their relationship to the ground
state structure, and their consequences on the analytic properties of the
dynamic structure function, directly from the theory.

The dynamic structure of 4He is actually quite rich as seen in Fig. 6 [30,
21]: “multi-excitations” are observed at energies above the single-excitation
dispersion curve discussed above. They can be related to combinations of
single-excitations, some of them are indicated in Fig. 6 by colored ellipses:
in blue, an extension of the linear phonon region (“ghost phonon”) deep into
the continuum; in red, an extension of the Pitaevskii-Plateau to small wave-
vectors; in green, above the roton minimum, a broadening of the dispersion
relation due to a Cherenkov effect; a red dashed ellipse at high energies
indicates a region of maxon-roton couplings. Multi-excitations are well
reproduced by recent theories. The mechanisms behind these are mostly the
same: It is kinematically possible that a perturbation of energy/momentum
$(E,{\bf k})$ can decay into two sharp quasiparticle excitations. For example,
the “ghost phonon” is caused by the effect that a perturbation can decay into
two almost collinear phonons, hence the effect disappears at about twice the
energy and wave number at which the dispersion relation $\epsilon(k)$ ceases
to be linear. Similarly, the “Pitaevskii-Plateau” is caused by the effect that
a perturbation of energy $2\Delta$ can decay into two rotons of energy
$\Delta$ and wave number $k_{\Delta}$. Hence, energy and momentum conservation
dictate that the plateau ends at a wave number of $2k_{\Delta}$, in this case
the two roton wave vectors are parallel. Perturbations with the same energy
and smaller wave numbers can decay into two rotons that are not collinear.
That means that the plateau can in principle extend to zero wave number; in
that case the perturbation decays into two anti-parallel rotons.

The Cherenkov effect is particularly interesting; it is caused by the fact
that the group velocity of a quasiparticle excitation beyond the roton
($R_{+}$ roton) is larger than the sound velocity. The effect is observed only
at low pressures.

Figure 6: Experimental (left) and theoretical (right) determinations of the
dynamic structure factor of superfluid 4He at zero pressure (upper graphs) and
near-solidification pressure (lower graphs) [30, 21]. In addition to the the
intense single-excitations dispersion seen in Fig. 5 (red lower curves), much
lower intensity multi-excitations (see color scale) can be observed here.
Colored ellipses indicate excitation couplings described in the text.

##  3 3He in 3D

The investigation of the dynamics of the 3He Fermi Liquid has naturally
followed that of the 4He Bose fluid [41, 20]. From the theoretical point of
view, the problem is notoriously difficult, due to the antisymmetry
requirement for the wave-functions of many-body fermionic systems [1, 2, 3].
Obtaining accurate results for the ground state is a challenge, even using
very sophisticated variational wave functions.

The dynamics is treated analogously to the boson case (4b); the local
excitation functions $\delta u_{n}({\bf r}_{1},\ldots,{\bf r}_{n};t)$ are
replaced by particle-hole operators

$\displaystyle\sum_{\genfrac{}{}{0.0pt}{1}{p_{1},\ldots,p_{n}}{h_{1}\ldots
h_{n}}}u_{p_{1},\ldots,p_{n};h_{1}\ldots h_{n}}(t)\times$
$\displaystyle{\vspace{-1cm}\quad\times a^{\dagger}_{p_{1}}\dots
a^{\dagger}_{p_{n}}a_{h_{n}}\dots a^{\dagger}_{h_{1}}}$

where the $p_{i}$ are the quantum numbers of unoccupied (“particle”) states
and the $h_{i}$ those of occupied (“hole”) states. Restricting the excitation
amplitudes to one-particle-one-hole components leads to a correlated version
of the random phase approximation (RPA) [42], the quantitative understanding
of the experiments requires at least 2-particle-2-hole amplitudes, see Refs.
[43, 44].

From the experimental point of view, even though the same inelastic neutron
techniques successfully applied to superfluid 4He can be used to investigate
liquid 3He, in practice the experiments are extremely difficult due to the
huge neutron absorption cross-section of the 3He nucleus [20]. In typical 4He
experiments, sample thicknesses are larger than 1 cm, while absorption limits
the maximum thickness in a liquid 3He measurement, to about 0.1 mm. In
addition, due to the lower excitation energies, experiments in 3He must be
performed at much lower temperatures ($\approx 0.1$ K) than in 4He ($\approx
1$ K), to remain in the low temperature limit, where few excitations are
present. For these reasons, in spite of the considerable interest in strongly
interacting fermions, accurate theoretical and experimental results have been
available only recently [45].

In order to understand the dynamics of a Fermi liquid, we begin by considering
the excitations of a Fermi gas. In the latter, excitations are created by
removing a particle from an occupied state inside the Fermi sphere, and
placing it outside the sphere, in a higher energy free state. The resulting
states are confined to a region of the energy vs wave-vector space called the
“particle-hole band”. The interacting fermionic system, i.e. the Fermi Liquid,
behaves in a similar way. A particle-hole band can be observed in bulk 3He
(see Fig 7), but in addition, collective modes are present. A density mode,
named “zero-sound” or “collision-less sound”, is observed. Its physical
origin, as noted by D. Pines [46], resides in the strong interactions, and
hence this mode is analogous to the phonon-roton mode described above for
bosons; statistics play here a minor role. Its substantial energy width
results from the possibility to decay into particle-hole excitations as it
enters the particle-hole band (Landau-damping mechanism).

The thermodynamic properties of bulk liquid 3He are often described, at very
low temperatures, by the phenomenological Landau Fermi Liquid model, where the
properties of the interacting system are related to those of the Fermi gas by
the introduction of an interaction function (“quasiparticle interaction”),
leading to a renormalization of parameters (effective mass, susceptibility
enhancement, etc.). Landau theory predicts, among other physical effects, the
existence of collective modes, the zero-sound and paramagnon modes.

Landau’s Fermi Liquid Theory is only applicable at very low temperatures. In
liquid 3He, where the Fermi energy, expressed in temperature units, is of the
order of several kelvins, the theory is valid below 0.1 K.

Figure 7: (Color online) Dynamic structure factor of bulk liquid 3He. Measured
density excitations [47] (triangles, upper branch) and spin-density
excitations [45] (squares, lower branch) are clearly visible (bars on the data
points indicate the observed width of the branch). Color intensity map, from
light blue to dark red: DMBT microscopic theory including exchange [44]. The
solid lines indicate the boundaries of the particle-hole band.

The reason becomes clear if one observes the actual excitation spectrum shown
in Fig. 7. The low energy mode, observed below the zero-sound excitation, is a
spin-density mode, completely immersed in the particle-hole band, and hence
severely Landau-damped. The origin of this mode, analogous to a dampened spin-
wave and often called “paramagnon”, can be traced back directly to the
antisymmetry of the many-body wave function.

There are also some interesting X-fay studies of the dynamic structure
function at high momentum transfers [48] that have led to a discussion about
the location of the particle-hole band [49, 50]. The authors of Ref. [48]
state that The obtained results show no evidence of such a decay: the zero-
sound mode remains well defined in the whole explored wave number range.
According to that, the particle-hole continuum should be much lower than the
continuum of the non-interacting Fermi system. An analysis of both the data
and the theoretical predictions within DMBT [51] shows that there is, when the
theoretical data are convoluted with the experimental resolution, actually no
contradiction to the conventional interpretation of the scenario.

Figure 8: The figure shows the normalized dynamic structure function
$S(k,\omega)/S(k)$ of 3He at a density of $\rho=0.0166$Å-3, for a sequence of
wave numbers. The theoretical results obtained by DMBT [51] are compared with
the X-ray scattering data [48] (squares) at saturated vapor pressure. Also
shown are the theoretical results convoluted with the experimental resolution
(dashed lines). The gray-shaded area shows the particle-hole continuum of a
non-interacting Fermi fluid. $t_{\rm F}=\hbar^{2}k_{\rm F}^{2}/2m$ is the
Fermi energy of the non-interacting system.

Details are shown in Fig. 8. To facilitate the comparison with experiments,
the theoretical spectra were convoluted with the experimental resolution.
Also, the results were scaled by $1/S(k)$ such that the integrated strength is
1 for all momentum transfers. The same scaling was applied to the experimental
data.

##  4 4He in reduced dimensions

The dynamics of quantum fluids has also been investigated in reduced
dimensions. Two-dimensional (2D) systems are obtained experimentally by
adsorption of gases onto solid substrates, usually graphite. One-dimensional
(1D) systems can also be created, they are obtained by confining the gases
inside silica or graphite nano-tubes [52, 53, 54]. There is also theoretical
interest [55, 56, 57, 58, 59] in the properties of 4He in one-dimension.

Extensive theoretical studies have clarified the nature of the excitations in
such systems. Superfluid 4He films of atomic thicknesses, for instance,
display in-plane density modes, as well as capillary waves.

Let us first discuss the dynamics of rigorously 2D 4He. Basically we see the
same effects as in 3D: a linear “phonon” branch as well as, at high momentum
transfers, a “roton minimum”.

Figure 9: (Color online) The figure shows contour plots of the dynamic
structure function of two-dimensional liquid 4He for a sequence of areal
densities as shown in the legends. The colors have been chosen to highlight
the prominent features, darker colors correspond to higher values of
$S(k,\hbar\omega)$. We also show the Feynman spectrum, and the simulation data
of Ref. [60]. From Ref. [61].

One can also see the somewhat finer details discussed above: At low densities,
the “ghost-phonon”, and at high densities the plateau coming down and
extending to long wave lengths. A second striking feature is the appearance of
a sharp mode below the plateau. We stress the difference: normally, the
plateau is a threshold above which a wave of energy/momentum $(\hbar\omega,k)$
can decay into two rotons. This has the consequence that the imaginary part of
the self-energy $\Sigma(k,\hbar\omega)$ is a step function and the real part
has a logarithmic singularity [36]. A collective mode is, on the other hand,
characterized by a singularity of the $S(k,\hbar\omega)$. Figs. 9 show, for
the two highest densities, the appearance of a sharp discrete mode below the
plateau. At a wave number of $k\approx 2.6\,$Å-1, the collective mode is still
merged into the continuum. With increasing wave number, we see, however, a
clearly distinguishable but rather weak ($Z(k)\approx 0.02$) mode about 0.3 K
below the plateau. We have mentioned already above that the roton should be
seen as an emergence of short–ranged order, the “ghost of a Bragg Spot” [62,
63]. If this is the case, then one might see a second Bragg spot. Such a thing
is not seen, but location of the secondary minimum in 2D corresponds indeed to
the position of a second Bragg spot of a triangular lattice.

Density modes are analogous to the phonons and rotons already mentioned in the
case of bulk helium, but due to the adsorption potential of the substrate, the
fluid is inhomogeneous, and solidification of the first two or three atomic
layers is observed in thick films. Layered-excitations are observed in this
regime.

Capillary waves, also called “ripplons”, have a different nature: these are
surface waves, which do not (in first approximation) correspond to
compression/expansion of the fluid, but to a change of height in the external
potential (gravitational in the case of a bulk helium surface, or substrate
potential in adsorbed films).

A typical example is shown in Fig. 10. It corresponds to a multilayer film of
about 6 atomic layers, where both density excitations and ripplons are visible
with comparable intensity [64]. The experimental curve is broader than the
theoretical one, an effect largely due to the mosaic spread (angular
distribution) of the graphite substrate, a powder made of microscopic flat
crystals.

Figure 10: Experimental [64] (top) and theoretical (bottom) [65]
determinations of the dynamic structure factor of a multilayer film of 4He,
displaying density and ripplon excitations.

Finally, in the sub-monolayer coverage regime, one can observe gaseous, fluid
and solid phases, as well as “commensurate phases” where the atoms adopt a
periodic arrangement commensurate with the substrate atomic lattice.

##  5 3He in reduced dimensions

Figure 11: Experimental (top) and theoretical (left) determinations of the
dynamic structure factor of two-dimensional liquid 3He [66]. The theoretical
calculation shows the density and spin-density modes; the latter is not
visible experimentally, it is masked by a strong elastic signal due to the
substrate. A comparison with Fig. 7 shows that the roton-like excitation in
the 2D system is found outside the particle-hole band (the region limited by
the blue solid lines), thus avoiding Landau damping.

The excitations in 3He films are shown in Fig. 11. They are very similar to
those discussed above for bulk liquid, with the remarkable exception [66] that
the roton-like excitations are at the edge of the particle-hole band and hence
not strongly affected by Landau damping. These excitations can therefore
propagate, suggesting interesting physical effects specific to 2D Fermi
liquids.

##  6 Summary

The investigation of the dynamics of quantum fluids has been performed
experimentally at very low temperatures using various techniques, mainly
specific heat, compressibility, neutron and X-ray scattering and, for 3He,
nuclear magnetic susceptibility. Theoretical calculations are based on semi-
phenomenological approaches [67, 68, 69], dynamic many body theory (DMBT) [43,
30, 44] as well as different versions of Quantum Monte Carlo numerical
calculations (PIMC, DMC) [27, 29, 60, 57] to cite only a few. These have
brought a very detailed understanding of these systems. It is sometimes very
difficult to perform experimental measurements on some particular systems (for
instance in reduced dimensionality), or even impossible to create quite
arbitrary systems (”mathematical models” of interaction or external
potentials, etc.). In this case the combination of numerical and microscopic
techniques is essential.

The results find a natural application to many different physical systems,
like nuclear matter, neutron stars, particle physics and cosmology, where
strongly interacting particles are investigated. Particles themselves, in
fact, can be interpreted as quantized excitations of an underlying field. The
helium liquids have been, in some sense, the Ariadne thread in these
investigations: Their interactions are well established [70, 71] and, compared
to nucleon interactions, quite simple. The density of these systems is very
high and they are very “quantum”, hence they do not permit the simple
treatments that have become popular by the experimental success to generate
and investigate cold quantum gases. The experimental challenge is equally high
and only the very best of the experimental equipment and the experimental
techniques can stand up to the challenge. Nevertheless, these challenges have
been met on both the theoretical and experimental side.

## References

  * [1] Fetter, A. L. and Walecka, J. D.,  Quantum Theory of Many-Particle Systems,  McGraw-Hill, New York, 1971.
  * [2] Thouless, D. J.,  The quantum mechanics of many-body systems,  Academic Press, New York, 2 edition, 1972.
  * [3] Pines, D. and Nozieres, P.,  The Theory of Quantum Liquids, volume I,  Benjamin, New York, 1966.
  * [4] Wilks, J.,  The Properties of Liquid and Solid Helium,  Clarendon, New York, Oxford, 1967.
  * [5] de Boer, J. and Michels, A.,  Physica 6 (1938) 945.
  * [6] Aziz, R. A., Slaman, M. J., Koide, A., Allnatt, A. R., and Meath, W. J.,  Mol. Phys. 77 (1992) 321337.
  * [7] de Boer, J.,  Physica 14 (1948) 139.
  * [8] Kapitza, P.,  Nature 141 (1938) 74.
  * [9] Misener, J. F. A. . A. D.,  Nature 141 (1938) 75.
  * [10] London, F.,  Nature 141 (1938) 643.
  * [11] Landau, L.,  J. Phys. U.S.S.R. 5 (1941) 71.
  * [12] Landau, L.,  J. Phys. U.S.S.R. 11 (1947) 91.
  * [13] Bogoliubov, N. N.,  J. Phys. U.S.S.R 11 (1947) 23.
  * [14] Feynman, R. P.,  Phys. Rev. 94 (1954) 262.
  * [15] Feynman, R. P. and Cohen, M.,  Phys. Rev. 102 (1956) 1189.
  * [16] Palevsky, H., Otnes, K., Larsson, K. E., Pauli, R., and Stedman, R.,  Phys. Rev. 108 (1957) 1346.
  * [17] Palevsky, H., Otnes, K., and Larsson, K. E.,  Phys. Rev. 112 (1958) 11.
  * [18] Glyde, H. R.,  Rep. Prog. Phys. 81 (2018) 014501.
  * [19] Lovesey, S. W.,  Theory of Neutron Scattering from Condensed Matter,  Clarendon Press, Oxford, 1984.
  * [20] Glyde, H. R.,  Excitations in liquid and solid helium,  Oxford University Press, Oxford, 1994.
  * [21] Beauvois, K. et al.,  Phys. Rev. B 97 (2018) 184520.
  * [22] Godfrin, H. et al.,  Phys. Rev. B 103 (2021) 104516.
  * [23] Cowley, R. A. and Woods, A. D. B.,  Can. J. Phys. 49 (1971) 177.
  * [24] Glyde, H. R., Gibbs, M. R., Stirling, W. G., and Adams, M. A.,  Europhysics Letters 43 (1998) 422.
  * [25] Jackson, H. W. and Feenberg, E.,  Rev. Mod. Phys. 34 (1962) 686.
  * [26] Lee, D. K. and Lee, F. J.,  Phys. Rev. B 11 (1975) 4318.
  * [27] Boronat, J. and Casulleras, J.,  Europhys. Lett. 38 (1997) 291.
  * [28] Boronat, J. and Casulleras, J.,  J. Low Temp. Phys. 110 (1998) 443.
  * [29] Ferré, G. and Boronat, J.,  Phys. Rev. B 93 (2016) 104510.
  * [30] Campbell, C. E., Krotscheck, E., and Lichtenegger, T.,  Phys. Rev. B 91 (2015) 184510/1.
  * [31] Jackson, H. W. and Feenberg, E.,  Ann. Phys. (NY) 15 (1961) 266.
  * [32] Jackson, H. W.,  Phys. Rev. 185 (1969) 186.
  * [33] Jackson, H. W.,  Phys. Rev. A 8 (1973) 1529.
  * [34] Jackson, H. W.,  Phys. Rev. A 9 (1974) 964.
  * [35] Manousakis, E. and Pandharipande, V. R.,  Phys. Rev. B 33 (1986) 150.
  * [36] Pitaevskii, L. P.,  Zh. Eksp. Theor. Fiz. 36 (1959) 1168,  [Sov. Phys. JETP 9, 830 (1959)].
  * [37] Campbell, C. E. and Krotscheck, E.,  Phys. Rev. B 80 (2009) 174501.
  * [38] Vitali, E., Rossi, M., Reatto, L., and Galli, D. E.,  Phys. Rev. B 82 (2010) 174510/1.
  * [39] Roggero, A., Pederiva, F., and Orlandini, G.,  Phys. Rev. B 88 (2013) 094302.
  * [40] Nava, M., Galli, D. E., Moroni, S., and Vitali, E.,  Phys. Rev. B 87 (2013) 144506/1.
  * [41] Dobbs, E. R.,  Helium Three,  Oxford University Press, Oxford, 2001.
  * [42] Krotscheck, E.,  Phys. Rev. A 26 (1982) 3536.
  * [43] Böhm, H. M., Holler, R., Krotscheck, E., and Panholzer, M.,  Phys. Rev. B 82 (2010) 224505.
  * [44] Krotscheck, E. and Lichtenegger, T.,  Dynamic many body theory IV: (Spin-)density fluctuations and exchange in normal-liquid 3He, 2022,  in preparation.
  * [45] Glyde, H. R. et al.,  Phys. Rev. B 61 (2000) 1421.
  * [46] Pines, D.,  Physics Today 34 (Nov. 1981) 106.
  * [47] Fåk, B., Guckelsberger, K., Körfer, M., Scherm, R., and Dianoux, A. J.,  Phys. Rev. B 41 (1990) 8732.
  * [48] Albergamo, F., Verbeni, R., Huotari, S., Vankó, G., and Monaco, G.,  Phys. Rev. Lett. 99 (2007) 205301/1.
  * [49] Schmets, A. J. M. and Montfrooij, W.,  Phys. Rev. Lett. 100 (2008) 239601/1.
  * [50] Albergamo, F., Verbeni, R., Huotari, S., Vankó, G., and Monaco, G.,  Phys. Rev. Lett. 100 (2008) 239602/1.
  * [51] Krotscheck, E. and Panholzer, M.,  J. Low Temp. Phys. 163 (2011) 1.
  * [52] Yano, H., Yoshizaki, S., Inagaki, S., Fukushima, Y., and Wada, N.,  J. Low Temp. Phys. 110 (1998) 573–578.
  * [53] Kato, H., Ishioh, K., Wada, N., Ito, T., and Watanabe, T.,  J. Low Temp. Phys. 68 (1987) 321.
  * [54] Maestro, A. D., Nichols, N. S., Prisk, T. R., Warren, G., and Sokol, P. E.,  Experimental realization of one dimensional helium, 2022,  arXiv:2203.11984.
  * [55] Stan, G. and Cole, M. W.,  Surface Science 395 (1998) 280.
  * [56] Stan, G., Crespi, V. H., Cole, M. W., and Boninsegni, M.,  J. Low Temp. Phys. 113 (1998) 447–452.
  * [57] Bertaina, G., Motta, M., Rossi, M., Vitali, E., and Galli, D. E.,  Phys. Rev. Lett. 116 (2016) 135302.
  * [58] Krotscheck, E. and Miller, M. D.,  Phys. Rev. B 60 (1999) 13038.
  * [59] Gordillo, M. C., Boronat, J., and Casulleras, J.,  Phys. Rev. B 61 (2000) R878.
  * [60] Arrigoni, F., Vitali, E., Galli, D. E., and Reatto, L.,  Fizika Nizkikh Temperatur 39 (2013) 1021.
  * [61] Krotscheck, E. and Lichtenegger, T.,  J. Low Temp. Phys. 178 (2015) 61.
  * [62] Nozières, P.,  J. Low Temp. Phys. 137 (2004) 45.
  * [63] Nozières, P.,  J. Low Temp. Phys. 142 (2006) 83.
  * [64] Lauter, H. J., Godfrin, H., Frank, V. L. P., and Leiderer, P.,  Phys. Rev. Lett. 68 (1992) 2484.
  * [65] Clements, B. E., Krotscheck, E., and Tymczak, C. J.,  Phys. Rev. B 53 (1996) 12253.
  * [66] Godfrin, H. et al.,  Nature 483 (2012) 576.
  * [67] Aldrich, C. H. and Pines, D.,  J. Low Temp. Phys. 25 (1976) 677.
  * [68] Aldrich III, C. H. and Pines, D.,  J. Low Temp. Phys. 31 (1978) 689.
  * [69] Hsu, W. and Pines, D.,  J. Stat. Phys. 38 (1985) 273.
  * [70] Aziz, R. A., Nain, V. P. S., Carley, J. C., Taylor, W. J., and McConville, G. T.,  J. Chem. Phys. 70 (1979) 4330.
  * [71] Aziz, R. A., McCourt, F. R. W., and Wong, C. C. K.,  Molec. Phys. 61 (1987) 1487.