#  Viscous damping in weltering motion of trapped hydrodynamic dipolar Fermi
gases

Reuben R. W. Wang  John L. Bohn  JILA, NIST, and Department of Physics,
University of Colorado, Boulder, Colorado 80309, USA

###### Abstract

We consider collective motion and damping of dipolar Fermi gases in the
hydrodynamic regime. We investigate the trajectories of collective
oscillations – here dubbed “weltering” motions – in cross-dimensional
rethermalization experiments via Monte Carlo simulations, where we find stark
differences from the dilute regime. These observations are interpreted within
a semi-empirical theory of viscous hydrodynamics for gases confined to
anisotropic harmonic potentials. The derived equations of motion provide a
simple effective theory that show favorable agreement with full numerical
solutions. To do so, the theory must carefully account for the size and shape
of the effective volume within which the gas’ behavior is hydrodynamic.
Although formulated for dipolar molecules, our theoretical framework retains a
flexibility to accommodate arbitrary elastic cross sections.

††preprint: APS/123-QED

##  I  Introduction

Suppression of two-body collisional losses has been crucial for achieving
stable samples of molecular quantum gases. Within the last decade, theoretical
and experimental advances have brought to fruition the electric field
shielding of polar molecules against chemical reaction and complex formation
Karman and Hutson (2018); Karman (2020); Anderegg  _et al._ (2021); Bigagli
_et al._ (2023); Quéméner and Bohn (2016); González-Martínez  _et al._ (2017);
Matsuda  _et al._ (2020); Li  _et al._ (2021); Lassablière and Quéméner
(2022); Bigagli  _et al._ (2023); Lin  _et al._ (2023), permitting the
production of degenerate bulk molecular samples Valtolina  _et al._ (2020);
Schindewolf  _et al._ (2022). But even before the onset of quantum degeneracy,
these shielded molecules present a long-lived versatile platform for exploring
dipolar physics Baranov (2008); Lahaye  _et al._ (2009); Chomaz  _et al._
(2022). For instance, dipole-dipole interactions lead to highly anisotropic
two-body collision cross sections Bohn and Jin (2014) and observable
anisotropy in the collective dynamics of thermal gases Aikawa  _et al._
(2014); Tang  _et al._ (2016); Sykes and Bohn (2015); Wang  _et al._ (2020);
Wang and Bohn (2021). For these nondegenerate bulk gases, thermalization is an
essential mechanism with great utility in applications such as evaporative
cooling Anderson  _et al._ (1995); Davis  _et al._ (1995); DeMarco and Jin
(1999); Thomas  _et al._ (2003); Griesmaier  _et al._ (2005); Lu  _et al._
(2011); Aikawa  _et al._ (2012); Marco  _et al._ (2019) and scattering length
measurements Monroe  _et al._ (1993); Newbury  _et al._ (1995); Schmidt  _et
al._ (2003); Tang  _et al._ (2015); Patscheider  _et al._ (2022). The accuracy
and efficacy of both these applications, in turn, rely on a deep understanding
of thermalization in such systems.

The difference between dilute and hydrodynamic limits is revealed clearly in a
gas’ response to perturbation. In particular, in a cross-dimensional
rethermalization experiment, an initially equilibrated gas is preferentially
heated along a particular axis, then allowed to rethermalize back to
equilibrium Monroe  _et al._ (1993). Thermalization in the dilute regime is
closely related to the collision rate Monroe  _et al._ (1993); Guéry-Odelin
_et al._ (1999); Wang and Bohn (2021), while the hydrodynamic regime sees
similarly extracted relaxation rates close to the trapping frequency instead
Ma  _et al._ (2003); Thomas  _et al._ (2003); Schindewolf  _et al._ (2022).
The difference between the two regimes is illustrated in Fig. 1. In both
panels, a collection of 23Na40K molecules is subjected to the same harmonic
trapping potential

$\displaystyle V(\bm{r})$
$\displaystyle=\frac{1}{2}m\sum_{i}\omega_{i}^{2}r_{i}^{2}.$ (1)

and subsequently excited along the $z$ axis. The only difference is the
molecule number: for fewer molecules in the upper panel (a), the dynamics is
dilute, while for a greater number of molecules in the lower panel (b), it is
hydrodynamic.

Figure 1:  Pseudotemperatures (2) obtained from Monte Carlo simulations in the
dilute (upper panel, a) and hydrodynamic (lower panel, b) regimes. The gas
consists of microwave shielded 23Na40K molecules with dipole moment $d=0.75$
D, oriented along $\hat{\bm{x}}$, at temperature $T=700$ nK. The gas is
initially excited along $z$ by an instantaneous trap frequency ramp to
$\omega_{z}=2\pi\times 147$ Hz, while $\omega_{x}=\omega_{y}=2\pi\times 82.5$
Hz remain constant. The regimes are differentiated by the number of molecules
$N$, which are $N=10^{4}$ in panel (a), and $N=2\times 10^{5}$ in panel (b).

In both cases, the behavior is tracked using time trace plots of the
pseudotemperatures ${\cal T}_{i}(t)$, shown in Fig. 1. A pseudotemperature is
defined along axis $i$ as Sykes and Bohn (2015)

$\displaystyle k_{B}{\cal
T}_{i}(t)\equiv\frac{1}{2}m\omega_{i}^{2}\\{r_{i}^{2}\\}(t)+\frac{1}{2}m\\{v_{i}^{2}\\}(t).$
(2)

where $\\{\ldots\\}(t)$ denotes the time varying ensemble average over
molecular positions $\bm{r}$ and velocities $\bm{v}$, $m$ is the molecular
mass, $k_{B}$ is Boltzmann’s constant. Details of the calculation that
produced this figure are provided below.

The dilute regime is characterized by collision rates small compared to the
trap frequencies. Hence in this case, pseudotemperature in the warm, $z$
direction gradually diminishes, while that in the other, cooler directions
gradually increases, until the gas equilibrates on the time scale shown. The
hydrodynamic gas, by contrast, behaves like a somewhat compressible fluid;
excitation initially in the $z$ direction is distributed almost immediately
into the other directions, and the resulting dynamics is more like the
irregular flow to and fro of this liquid about its stationary center of mass.
The fluid expands sometimes in the radial direction, sometimes in the axial
direction, with irregularly varying amplitudes, reminiscent of waves on an
unquiet ocean. We therefore refer to this form of collective fluid excitation
as weltering. 111 Thus in a tense moment for Ahab aboard a whale boat,
Melville relates that the “… White Whale darted through the weltering ocean.”
(Moby Dick, Chap. 135.)

In the dilute gas case, the primary response of the gas is to come to thermal
equilibrium, whereby its dynamics is largely summarized in a single, density-
normalized equilibration rate, whose inverse defines the “number of collisions
per rethermalization” Monroe  _et al._ (1993). For dipolar gases, this
quantity can depend on the orientation of the dipoles relative to the
excitation axis Bohn and Jin (2014); Wang and Bohn (2021). Vice versa, the
complex dynamics of the hydrodynamic fluid requires a more complete
theoretical description.

The purpose of this paper is to provide such a description. We will base full
dynamics on a Monte Carlo simulation, to further elaborate the difference
between dilute and hydrodynamic regimes. Further, we will develop a simplified
formulation based on a Gaussian ansatz for the width of a gas, which semi-
empirically reproduces the numerics. Key to this model is the realization that
the periphery of a harmonically trapped gas is always dilute Schäfer (2014);
Bluhm and Schäfer (2015), which necessitates defining an effective volume
inside which hydrodynamics is a good idea. We identify the dependence of this
volume on the anisotropy of the trap and of the collision cross section among
polarized dipoles. Our theory is also presented in a manner that accommodates
arbitrary elastic cross sections, opening its applicability to a broader
variety of ultracold molecular gas experiments with far from threshold
collisions Sadeghpour  _et al._ (2000).

The remainder of this paper is organized as follows: In Sec. II, we describe
the numerical tools adopted to study trapped hydrodynamic gases, and present
notable differences from the dilute limit. We then introduce the equations of
motion employed to model a nondegenerate hydrodynamic dipolar gas in Sec. IV,
with the assumption of threshold scattering. A variational ansatz is employed
in Sec. IV.1, to derive effective dynamical equations governing weltering
oscillations in a harmonic trap. A comparison of our theory to full numerical
solutions is presented in Sec. IV.3, from which we purport several
considerations about the hydrodynamic extent of gases in traps. Finally,
conclusion are drawn in Sec. V, along with possible extensions of this current
work.

##  II  Numerical method

A gas is said to be hydrodynamic when the molecular mean-free path is much
smaller than the characteristic length over which fluid flow occurs Huang
(1963). The ratio of these scales is given by the Knudsen number ${\rm Kn}$.
For a harmonically trapped gas with mean density $\langle
n\rangle=\frac{1}{N}\int n^{2}(\bm{r})d^{3}r$ and molecules with total cross
section $\sigma_{\rm coll}$, the mean-free path is given by $L=(\langle
n\rangle\sigma_{\rm coll})^{-1}$. With a given geometric mean frequency
$\overline{\omega}$ and temperature $T$, the thermal width of the gas is
$R_{\rm th}=\sqrt{k_{B}T/m\overline{\omega}^{2}}$.

Alternatively, the Knudsen number can also be written as the ratio of mean
trapping frequency over the collision rate $\gamma_{\rm coll}=\langle
n\rangle\sigma_{\rm coll}\langle v_{\rm coll}\rangle$, where $\langle v_{\rm
coll}\rangle=\sqrt{16k_{B}T/(\pi m)}$ is the mean collision velocity.
Explicitly, these relations are summarized as

$\displaystyle{\rm Kn}$ $\displaystyle=\frac{L}{R_{\rm
th}}=\frac{4\>\overline{\omega}}{\pi^{\nicefrac{{1}}{{2}}}\gamma_{\rm
coll}}=\frac{8\pi^{\nicefrac{{3}}{{2}}}k_{B}T}{Nm\overline{\omega}^{2}\sigma_{\rm
coll}}.$ (3)

A trapped gas is said to be hydrodynamic if ${\rm Kn}\ll 1$. The relations
above provide an approximate mean Knudsen number. In practice, the thermal
width can differ in directions with different trap frequencies, while the
cross section, for dipolar scattering, can depend on the direction of the
collisions axis. Thus the boundary between hydrodynamic and dilute flow can be
anisotropic, a topic to be dealt with below.

To compute dynamics in either regime, we utilize the direct simulation Monte
Carlo (DSMC) method Bird (1970) to obtain numerical solutions to the Boltzmann
equation. In doing so, these numerical simulations allow for explorations of
hydrodynamic phenomena, while later also serving as a benchmark for our semi-
empirical theory.

The DSMC implementation we adopt for this work follows very closely that
described in Refs. Sykes and Bohn (2015); Wang  _et al._ (2020), which study
similar systems but in the dilute regime. Described briefly, the Boltzmann
equation is solved by approximating the phase space distribution with a
discrete ensemble of $N$ molecules

$\displaystyle
f(\bm{r},\bm{v})\approx\sum_{k=1}^{N}\delta^{3}(\bm{r}-\bm{r}_{k})\delta^{3}(\bm{v}-\bm{v}_{k}).$
(4)

Most crucial to an accurate hydrodynamic simulation is that collisions are
handled adequately. The DSMC does so by constructing a discrete spatial grid
within the simulation volume, binning particles into each grid cell based on
their positions, then sampling their collisional interactions from a
probability distribution derived from the differential cross section Sykes and
Bohn (2015).

Choosing a uniform grid that is appropriate for maintaining accuracy and
computational efficiency becomes tricky at large collision rates, so we
utilize a locally adaptive discretization scheme instead. At every numerical
time step, the locally adaptive grid is built in two phases. Phase one
constructs a master grid, consisting of uniform volume cells that span the
simulation volume. The resolution of the grid is then refined in phase two,
with an octree algorithm Franklin and Akman (1985). The octree algorithm
further discretizes the simulation volume by recursively subdividing cells
into eight octants, terminating when each cell has at most
$N_{\text{cell}}^{\max}$ particles. The parameter $N_{\text{cell}}^{\max}$, is
initialized at the start of the simulation, which we optimize for stochastic
convergence.

##  III  Numerical Results

Table 1:  Table of parameter values utilized in the Monte Carlo simulation for fermionic 23Na40K molecules. Da $=1.661\times 10^{-27}$ kg stands for Dalton (atomic mass unit) and ${\rm D}=3.33564\times 10^{-30}$ C$\cdot$m is a Debye.  Parameter | Symbol | Value | Unit
---|---|---|---
Relative molecular mass | $M_{r}$ | 63 | Da
Electric dipole moment | $d$ | 0.75 | D
Initial gas temperature | ${\cal T}(0)$ | 700 | nK
Trap frequency geometric mean | $\overline{\omega}$ | $2\pi 100$ | Hz

For our numerical experiments, we envision an ultracold gas of microwave
shielded 23Na40K molecules with the parameters in Tab. 1. The initial
temperature is chosen such that the gas remains nondegenerate with $T>T_{F}$
Butts and Rokhsar (1997) for all values of ${\rm Kn}$ in consideration, and
the trap is assumed cylindrically symmetric with
$\omega_{x}=\omega_{y}\equiv\omega_{\perp}$ but
$\omega_{\perp}\neq\omega_{z}$. Key variables of interest to this study will
be: a) the number of molecules $N$, which affects ${\rm Kn}$ and therefore how
hydrodynamic the gas is; b) the trap anisotropy
$\lambda=(\omega_{z}/\omega_{\perp})^{2}$; c) and the dipole orientation
$\hat{\bm{{\cal E}}}$. For the sake of illustration, collision cross sections
are described by the analytical formulas for point dipoles given in Ref. Bohn
and Jin (2014), although at sufficient temperature, realistic cross sections
may differ from these. For convenience, we only allow $\hat{\bm{{\cal E}}}$ to
tilt within the $x,z$-plane, allowing us to define a dipole tilt angle
$\Theta=\cos^{-1}\hat{\bm{{\cal E}}}\cdot\hat{\bm{z}}$, that parameterizes the
collisional anisotropy.

The behavior of the fluid after excitation in the $z$ direction is shown in
Fig. 2. This is done in a prolate (cigar) trap with $\lambda=0.2$, containing
$N=5\times 10^{5}$ molecules, with Knudsen number ${\rm Kn}\approx 0.04$. This
figure plots the separated position and momentum space pseudotemperatures
${\cal T}_{r_{i}}(t)=m\omega_{i}^{2}\\{r_{i}^{2}\\}(t)/k_{B}$ and ${\cal
T}_{v_{i}}(t)=m\\{v_{i}^{2}\\}(t)/k_{B}$ respectively. The position space time
trace shows the clear out-of-phase oscillations between the widths in the
radial and axial directions, expected for a weltering fluid. The momentum
space time trace has oscillations of considerably smaller magnitude than
${\cal T}_{r_{i}}$, and also shows a phasing in oscillations amongst the
different ${\cal T}_{v_{i}}$ traces. These observations showcase how large
collision rates diminish the effect of out-of-equilibrium thermodynamics on
the hydrodynamic welter of the gas.

Figure 2:  Plots of the ${\cal T}_{r_{i}}$ (upper panel a) and ${\cal
T}_{v_{i}}$ (lower panel b) vs time from a cross-dimensional rethermalization
experiment, with excitation along $z$. The gas is hydrodynamic with $N=5\times
10^{5}$ (${\rm Kn}\approx 0.04$), $\lambda=0.2$ and the parameters in Tab. 1.

The difference between dilute and hydrodynamic regimes is sharpened by
comparing the dependence of dynamics on the tilt angle $\Theta$ of the
dipoles. To this end, Fig. 3 plots the three components of pseudotemperature
${\cal T}_{i}$ for the dilute (upper row) and hydrodynamic (lower row) gases,
at the 3 different dipole tilt angles
$\Theta=0^{\circ},45^{\circ},90^{\circ}$.

As anticipated in Fig. 1, the dilute gas responds to the excitation primarily
by melting back to thermal equilibrium while the hydrodynamic gas exhibits
radial weltering motion, resulting from oscillating fluid flow toward and away
from the trap center. In Fig. 3 a second dramatic difference appears. For the
dilute gas, with the dipoles tilted away from the axis of trap symmetry ($z$),
the rates of warming of the gas in the $x$ and $y$ directions differ, as a
consequence of the anisotropic scattering cross section Bohn and Jin (2014);
Sykes and Bohn (2015); Wang and Bohn (2021). By contrast, the excitations in
the $x$ and $y$ directions in the hydrodynamic regime are nearly equal. In the
hydrodynamic regime, relatively rapid collisions scramble memory of the dipole
orientation. Note that a slight difference in $x$ and $y$ motions occurs, due
to a residual anisotropy of the viscosity tensor, described in the next
section. Nevertheless, this anisotropy is not a main driving force in the
dynamics. It is true, however, that the overall damping rate of the weltering
excitations does depend on the dipole tilt angle, as will be elaborated upon
in Sec. IV.3 below.

Figure 3:  Pseudotemperature times traces ${\cal T}_{x}(t)$ (solid green
curves), ${\cal T}_{y}(t)$ (dashed blue curves) and ${\cal T}_{z}(t)$ (dotted
red curves) for 3 values of $\Theta=0^{\circ},45^{\circ},90^{\circ}$, in
subplots (a, d), (b, e) and (c, f) respectively. The 2 rows are differentiated
by the number of molecules, with the upper row (subplots a, b, c) having
$N=2\times 10^{3}$ (${\rm Kn}\approx 11.10$), while the lower row (subplots d,
e, f) has $N=3\times 10^{5}$ (${\rm Kn}\approx 0.07$). The experimental
parameters are those in Tab. 1 with $\lambda=0.2$. Note that the simulation
times are different between the upper ($t=0$ to $0.1$s) and lower ($t=0$ to
$0.04$s) rows.

##  IV  Hydrodynamic formulation

The Monte Carlo simulation, while accurate, is nevertheless somewhat
cumbersome for calculating the response of the gas. For this reason, in the
hydrodynamic regime, it is useful to formulate the fluid’s motion directly in
terms of hydrodynamics. When hydrodynamic, a nondegenerate gas behaves as a
thermoviscous fluid Griffin  _et al._ (1997); Nikuni and Griffin (1998);
Kavoulakis  _et al._ (1998) with thermal conductivity $\kappa_{ij}$, and
viscosity $\mu_{ijk\ell}$, which are, in general, coordinate dependent and
formulated as rank-2 and rank-4 tensors respectively De Groot and Mazur
(2013). The equations of motion of the fluid are Fetter and Walecka (2003):

$\displaystyle\frac{\partial\rho}{\partial t}+\sum_{j}\partial_{j}\left({\rho
U_{j}}\right)$ $\displaystyle=0,$ (5a) $\displaystyle\frac{\partial}{\partial
t}\left(\rho U_{i}\right)+\sum_{j}\partial_{j}\left(\rho U_{j}U_{i}\right)$
$\displaystyle=-\partial_{i}\left(nk_{B}T\right)-n\partial_{i}V(\bm{r})$
$\displaystyle+\sum_{j,k,\ell}\partial_{j}\left(\mu_{ijk\ell}\partial_{\ell}U_{k}\right),$
(5b) $\displaystyle\frac{\partial}{\partial t}(\rho
T)+\sum_{j}\partial_{j}\left(\rho TU_{j}\right)$
$\displaystyle=-\frac{2}{3}\rho T\sum_{i}\partial_{i}U_{i}$
$\displaystyle+\frac{2m}{3k_{B}}\sum_{i,j,k,\ell}(\partial_{j}U_{i})\mu_{ijk\ell}(\partial_{\ell}U_{k})$
$\displaystyle+\frac{2m}{3k_{B}}\sum_{i,j}\partial_{i}\left(\kappa_{ij}\partial_{j}T\right).$
(5c)

These equations govern the dynamics of the velocity averaged field variables
of mass density, flow velocity and temperature:

$\displaystyle\rho(\bm{r},t)$ $\displaystyle=mn(\bm{r},t)=\int
d^{3}vf(\bm{r},\bm{v},t)m,$ (6a) $\displaystyle\bm{U}(\bm{r},t)$
$\displaystyle=\frac{1}{n(\bm{r},t)}\int d^{3}vf(\bm{r},\bm{v},t)\bm{v},$ (6b)
$\displaystyle T(\bm{r},t)$ $\displaystyle=\frac{2}{3n(\bm{r},t)k_{B}}\int
d^{3}vf(\bm{r},\bm{v},t)\frac{1}{2}m\bm{u}^{2},$ (6c)

where $f(\bm{r},\bm{v},t)$ denotes the phase space distribution of the
molecules and $\bm{u}(\bm{r})=\bm{v}-\bm{U}(\bm{r})$ is the comoving molecular
velocity, relative to the frame of fluid flow.

It is worth pointing out that the local fluid kinetic temperature is related
to the flow velocity via

$\displaystyle\frac{3}{2}n(\bm{r},t)k_{B}T(\bm{r},t)$ $\displaystyle=\int
d^{3}vf(\bm{r},\bm{v},t)\frac{1}{2}m\bm{v}^{2}$
$\displaystyle\quad-\frac{1}{2}\rho\bm{U}(\bm{r},t)^{2},$ (7)

where the integral term is the local kinetic energy density. This relation
emphasizes a central difference between dilute and hydrodynamic trapped gases:
temperature, in the sense of equilibrium thermodynamics, is well defined
throughout the entire dynamical evolution when hydrodynamic, but only upon
global equilibration when dilute. Such a distinction identifies time-of-flight
imaging, common to ultracold gas experiments, as an indirect form of
thermometry to hydrodynamic gases, that probes an ensemble averaged sum of
both the fluid local temperature and mechanical energy from flow.

In this work, we assume that the transport tensors arise from two-body
collisions with elastic differential cross section $d\sigma/d\Omega$, as
derived with the first-order Chapman-Enskog method Chapman and Cowling (1990);
Wang and Bohn (2022a, b). We shall later see that only viscosity is relevant
to this work, so we omit further details of the thermal conductivity. At this
level of approximation, the anisotropic viscosity tensor for arbitrary
$d\sigma/d\Omega$ works out to be density independent, and is given as Wang
and Bohn (2022b, 2023)

$\displaystyle\bm{\mu}$
$\displaystyle=-\frac{2}{\beta}\left(\frac{n}{m\beta}\right)^{2}\left(\int
d^{3}u\bm{W}(\bm{u})\otimes C[f_{0}\bm{W}]\right)^{-1},$ (8)

where $\beta=(k_{B}T)^{-1}$ is the usual inverse temperature,

$\displaystyle\bm{W}$
$\displaystyle=\bm{u}\bm{u}^{T}-\frac{1}{3}\bm{u}^{2}\bm{I},$ (9)

is a rank-2 comoving velocity tensor, and $\bm{I}$ is the identity matrix. The
collision integrals

$\displaystyle C[f_{0}\bm{W}]$ $\displaystyle=\int
d^{3}u_{1}\absolutevalue{\bm{u}-\bm{u}_{1}}f_{0}(\bm{u})f_{0}(\bm{u}_{1})\int
d\Omega^{\prime}\frac{d\sigma}{d\Omega^{\prime}}\Delta\bm{W},$ (10)

with $\Delta\bm{W}=\bm{W}^{\prime}+\bm{W}^{\prime}_{1}-\bm{W}-\bm{W}_{1}$ and
primes denoting post-collision quantities, are evaluated with the Maxwell-
Boltzmann equilibrium phase space distribution function $f_{0}(\bm{u})$ Reif
(2009). The symbol $\otimes$ denotes a dyadic product which takes two tensors
of rank $N_{1}$ and $N_{2}$, and forms a tensor of rank $N_{1}+N_{2}$ (e.g.
$A_{ij}\otimes B_{k\ell}=C_{ijk\ell}$). Of interest here is the anisotropic
cross section resultant from close-to-threshold scattering 222 With the
parameters in Tab. 1, the typical collision energies are in fact much larger
than the dipole energy, which may place the assumption of threshold scattering
on shaky grounds. Nevertheless, we proceed with the threshold approximation in
this work. See App. B for further discussions.  between ultracold fermionic
polar molecules or dipolar atoms Matsuda  _et al._ (2020); Schindewolf  _et
al._ (2022); Patscheider  _et al._ (2022); Tang  _et al._ (2015). At low
enough temperatures with electric fields that align the dipoles along
$\hat{\bm{{\cal E}}}$, dipolar scattering is energy independent and permits
the viscosity tensor to be computed analytically Wang and Bohn (2023). It is
this analytic viscosity tensor that we use below.

###  IV.1  Viscous damping of a trapped fluid

The fluid equations in (5) are highly nonlinear and, in general, require
numerical methods to obtain solutions. For our purposes, we instead adopt a
variational ansatz approach to solving these partial differential equations
333 Our method has similarities to the scaling ansatz employed in Ref. Kagan
_et al._ (1997); Pedri  _et al._ (2003); Wächtler  _et al._ (2017), but
formulated slightly differently and with the inclusion of transport tensors. .
External confinement from a harmonic potential results in the equilibrium
(denoted by subscript $0$) density distribution following

$\displaystyle\rho_{0}(\bm{r})$
$\displaystyle=\frac{mN}{Z}\exp(-\frac{V(\bm{r})}{k_{B}T_{0}}),$ (11)

where $Z=\int d^{3}r{\rm e}^{-\frac{V(\bm{r})}{k_{B}T_{0}}}$ gives the
appropriate normalization and $N$ is the number of molecules. If we were then
only to consider collective oscillations and damping from long wavelength
excitations that do not induce center-of-mass sloshing, Eq. (11), motivates a
Gaussian variational ansatz for the local density:

$\displaystyle\rho(\bm{r},t)=mN\prod_{i=1}^{3}\frac{1}{\sqrt{2\pi\sigma_{i}^{2}(t)}}\exp\left(-\frac{r_{i}^{2}}{2\sigma_{i}^{2}(t)}\right),$
(12)

where $\sigma_{i}(t)$ is the distribution widths along each axis $i$ that we
allow to vary in time (depicted in Fig. 4).

Figure 4:  Cartoon of a density slice along axis $r_{i}$, through the Gaussian
ansatz for $\rho(\bm{r},t)$ with time varying widths $\sigma_{i}(t)$.

Plugging the ansatz of Eq. (12) into the continuity equation (5a) gives

$\displaystyle\sum_{i=1}^{3}\bigg{[}\partial_{i}U_{i}(\bm{r})$ $\displaystyle-
U_{i}(\bm{r})\left(\frac{r_{i}}{\sigma_{i}^{2}(t)}\right)$
$\displaystyle\quad+\left(\frac{r_{i}^{2}}{\sigma_{i}^{2}(t)}-1\right)\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\bigg{]}=0,$
(13)

which admits the velocity field solution

$\displaystyle
U_{i}(\bm{r})=\left(\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\right)r_{i}.$
(14)

Thus, as expected, the fluid flow vanishes in the trap’s center for the
excitations we consider. These functional forms for $\rho$ and $\bm{U}$ then
render the Navier-Stokes equation (5), of the form

$\displaystyle\ddot{\sigma}_{i}(t)+\omega_{i}^{2}\sigma_{i}(t)$
$\displaystyle=\frac{k_{B}}{m}\left(\frac{1}{\sigma_{i}(t)}-\frac{\sigma_{i}(t)}{r_{i}}\partial_{i}\right)T(\bm{r},t)$
$\displaystyle\quad+\sigma_{i}\sum_{j,k,\ell}\frac{\partial_{j}\mu_{ijk\ell}(T)}{r_{i}\rho(\bm{r})}\delta_{k,\ell}\frac{\dot{\sigma}_{k}}{\sigma_{\ell}},$
(15)

which bears no dependence on the thermal conductivity. Since $\sigma_{i}(t)$
does not depend on spatial coordinates, consistency requires that we take a
spatial average to suppress local fluctuations of the temperature field in Eq.
(IV.1). This average is taken by multiplying Eq. (IV.1) and the temperature
balance equation (5), by $n(\bm{r},t)$, then integrating over $d^{3}r$. App. A
gives further details of the spatial averaging procedure, which results in

$\displaystyle\ddot{\sigma}_{i}(t)$
$\displaystyle+\omega_{i}^{2}\sigma_{i}(t)+\frac{1}{3\sigma_{i}(t)}\sum_{j}\left[\omega_{j}^{2}\sigma_{j}^{2}(t)+\dot{\sigma}_{j}^{2}(t)\right]-\frac{2k_{B}T_{0}}{m\sigma_{i}(t)}\approx-\frac{2}{5}\frac{{\cal
V}_{\rm
hy}}{Nm}\sum_{j}\frac{\mu_{iijj}(T(t))}{\sigma_{i}(t)}\frac{\dot{\sigma}_{j}(t)}{\sigma_{j}(t)}.$
(16)

The relevant viscosity matrix elements can be recast in terms of a unit-free
matrix

$\displaystyle M_{ij}(\Theta)$
$\displaystyle\equiv\frac{\mu_{iijj}(T;\Theta)}{\mu_{0}(T)}$
$\displaystyle=\frac{1}{512}\begin{pmatrix}117\cos(4\Theta)+84\cos(2\Theta)+415&-28(3\cos(2\Theta)+11)&-(117\cos(4\Theta)+107)\\\\[2.5pt]
-28(3\cos(2\Theta)+11)&616&28(3\cos(2\Theta)-11)\\\\[2.5pt]
-(117\cos(4\Theta)+107)&28(3\cos(2\Theta)-11)&117\cos(4\Theta)-84\cos(2\Theta)+415\end{pmatrix},$
(17)

as is taken from Ref. Wang and Bohn (2023), where the isotropic unit-full
viscosity coefficient is given by Chapman and Cowling (1990)

$\displaystyle\mu_{0}(T)$
$\displaystyle=\frac{5}{16a_{d}^{2}}\sqrt{\frac{mk_{B}T}{\pi}}.$ (18)

With the parameters in Tab. 1, the isotropic viscosity has a value of
$\mu_{0}\approx 2.5\times 10^{-15}$ Pa$\cdot$s, which is around $10^{10}$
times less than air at room temperature and pressure Kadoya  _et al._ (1985).
The $M_{ij}(\Theta)$ matrix elements are plotted in Fig. 5, with components
coupled to the $x$ and $z$ axes showcasing a significant variation with
$\Theta$. We see in Fig. 5 that the magnitude of off-diagonal matrix elements
$M_{13}=M_{xz}$ and $M_{23}=M_{yz}$ become maximally separated around
$\Theta\approx 45^{\circ}$, explaining the slight separation of ${\cal
T}_{x}(t)$ and ${\cal T}_{y}(t)$ in Fig. 3, otherwise negligible when
$\Theta=0^{\circ},90^{\circ}$.

Figure 5:  $M_{ij}$ matrix elements as a function of $\Theta$. The diagonal
elements are plotted on the left in subplot (a), whereas the negated
(multiplied by a minus sign) off-diagonal elements plotted on the right in
subplot (b).

Eq. (16) above treats the temperature field appearing in $\mu_{ijk\ell}(T)$ to
be spatially uniform over the region where the gas is hydrodynamic. Such an
approximation follows from the form of collective oscillations implied by the
density (12) and flow velocity fields (14) in an initially isothermal gas,
disallowing a spatial temperature variation on the order of the gas spatial
widths Kavoulakis  _et al._ (1998); Schäfer (2014). Hence, temperature as
appears in the viscosity is simply treated as $T\approx T(t)$. In doing so, we
were required to define an effective hydrodynamic volume ${\cal V}_{\rm
hy}=\int d^{3}r$ 444Refs. Kavoulakis  _et al._ (1998); Schäfer (2007) provides
another estimate for the effective hydrodynamic volume which, however, we find
systematically underestimates the results from our numerical Monte Carlo
simulations. . Proper identification of this volume, including its dependence
on aspect ratio, density, and dipole tilt, is essential to the performance of
the model, and is our main undertaking here. We define this volume to be the
spheroidal volume bounded by the outer classical turning radius of the trap,
multiplied by an empirical factor $\eta$. The outer turning radius is obtained
by equating $E_{\rm total}=V(R_{\rm HD},\theta,\phi)$, to give (see App. A)

$\displaystyle R_{\rm HD}^{2}(\theta)$
$\displaystyle=\frac{6k_{B}T(t)}{m\omega_{\perp}^{2}}\left[\sin^{2}\theta+\lambda\cos^{2}\theta\right]^{-1},$
(19)

where $\lambda=(\omega_{z}/\omega_{\perp})^{2}$ quantifies the trapping
anisotropy. The effective hydrodynamic volume is then computed as

$\displaystyle{\cal V}_{\rm hy}(\lambda,{\rm Kn})$
$\displaystyle=\frac{\eta(\lambda,{\rm Kn})}{3}\int R_{\rm
HD}^{3}(\Omega)d\Omega$
$\displaystyle=\frac{4\pi}{3}\left(\frac{6k_{B}T(t)}{m\omega_{\perp}^{2}}\right)^{3/2}\frac{\eta(\lambda,{\rm
Kn})}{\sqrt{\lambda}}.$ (20)

As written, we have assumed that $\eta$ could depend on the trapping geometry
through $\lambda$ and on the Knudsen number, which in turn, also implicitly
depends on $N$ and the dipole angle $\Theta$. These dependencies are addressed
later in the paper. Such generality allows $\eta$ to act as a coarse-graining
parameter which accounts for all non-hydrodynamic effects excluded from our
current theoretical treatment. Additionally, Eq. (18) implies the temperature
dependence of viscosity goes as $\mu_{iijj}(T)\propto\sqrt{T}$, for which we
will simply approximate as $T\approx T_{0}$ for all times 555 Incorporating
time-dependence in the temperature $T(t)$, requires more sophisticated
treatments such as second-order hydrodynamics Lewis and Romatschke (2017). We
leave this task to a subsequent study. .

For the relevance of time-of-flight imaging, we point out that the momentum
space temperature, which differs from the local temperature of Eq. (6c), can
also be obtained from solutions to Eq. (16) via the relation

$\displaystyle k_{B}T_{p}(t)$ $\displaystyle=\frac{1}{3N}\int
d^{3}rd^{3}vf(\bm{r},\bm{v},t)m\bm{v}^{2}$
$\displaystyle=2k_{B}T_{0}-\frac{1}{3}\sum_{i}m\omega_{i}^{2}\sigma_{i}^{2}(t),$
(21)

as follows from Eqs. (IV), (14) and (40).

###  IV.2  Linear analysis

Some proceeding discussions on collective dynamics are made more accessible in
the language of normal modes, motivating a linear analysis of Eq. (16). If
only taken perturbatively out-of-equilibrium, we can consider small deviations
away from the equilibrium widths by writing
$\sigma_{i}(t)=\sigma_{0,i}+\delta\sigma_{i}(t)$. Then expanding to first-
order in $\delta\sigma_{i}(t)$, Eq. (16) becomes

$\displaystyle\ddot{\delta\sigma}_{i}(t)+2\sum_{j}\Gamma_{ij}\dot{\delta\sigma}_{j}(t)+\sum_{j}O_{ij}\delta\sigma_{j}(t)$
$\displaystyle\approx 0,$ (22)

with squared-frequency and damping matrices

$\displaystyle O_{ij}$
$\displaystyle=2\omega_{i}^{2}\delta_{i,j}+\frac{2}{3}\omega_{i}\omega_{j},$
(23a) $\displaystyle\Gamma_{ij}$ $\displaystyle=\frac{\mu_{0}{\cal V}_{\rm
hy}}{5Nk_{B}T_{0}}\omega_{i}M_{ij}(\Theta)\omega_{j}.$ (23b)

The matrices above encode the anisotropies from both the trap and anisotropic
collisions. A factor 2 multiplies $\bm{\Gamma}$ in Eq. (22) as is convention
in damped harmonic oscillators. With $\bm{\Gamma}$ multiplying the first-order
time derivative terms $\dot{\delta\sigma}_{i}$, it is made clear that damping
of weltering oscillations results from the trap frequency weighted viscosities
within the hydrodynamic volume.

Diagonalizing the squared-frequency matrix $\bm{O}$ gives the eigenvalues

$\displaystyle\omega_{0}^{2}$ $\displaystyle=2\omega_{\perp}^{2},$ (24a)
$\displaystyle\omega^{2}_{\pm}$
$\displaystyle=\frac{1}{3}\left(4\lambda+5\pm\sqrt{16\lambda^{2}-32\lambda+25}\right)\omega_{\perp}^{2},$
(24b)

which are exactly those obtained for inviscid Euler flow in Refs. Griffin  _et
al._ (1997); Kavoulakis  _et al._ (1998), and correspond to the respective
eigenmodes (up to arbitrary normalization)

$\displaystyle\bm{o}_{0}$ $\displaystyle=\begin{pmatrix}1\\\ -1\\\
0\end{pmatrix},$ (25a) $\displaystyle\bm{o}_{\pm}$
$\displaystyle=\begin{pmatrix}5-4\lambda\pm\sqrt{25+16\lambda(\lambda-2)}\\\
5-4\lambda\pm\sqrt{25+16\lambda(\lambda-2)}\\\ 4\sqrt{\lambda}\end{pmatrix}.$
(25b)

The eigenmode $\bm{o}_{0}$ is a strictly radial quadrupole mode, while
$\bm{o}_{-}$ and $\bm{o}_{+}$ are 3-dimensional quadrupole and breathing modes
respectively.

Similarly, $\bm{\Gamma}$ results in two nontrivial eigenvalues $\gamma_{\pm}$,
that constitute the principle damping rates of the system. Although it is
tempting to assign one of these principle rates as the overall relaxation
rate, the eigenmodes associated to each $\gamma_{\pm}$, are in general, not
the eigenmodes of $\bm{O}$. Consequently, coupling between the eigenmodes of
$\bm{\Gamma}$ is inevitable during dynamical evolution, enforcing that
accurate relaxation trajectories are best obtained from full solutions to Eq.
(22).

###  IV.3  The hydrodynamic volume

Returning to the main argument, Eq. (16) is expected to be a reasonable
representation of dynamics, provided the shape of the gas remains nearly
Gaussian. To employ these equations, we must establish the value of the
effective hydrodynamic volume. A first guess at this volume is given in Eq.
(IV.1), which left available a free parameter $\eta$, that may depend on
$\lambda$ and ${\rm Kn}$. As noted in Sec. IV.1, ${\rm Kn}$ is implicitly
dependent on $N$ and $\Theta$, which are taken as the relevant independent
variables for this study.

To extract $\eta$, we perform multiple DSMC runs while varying $\lambda$, $N$
and $\Theta$, which provides us time traces of $T_{p}(t)$ (IV.1) for each
combination of parameter values. We then fit $T_{p}(t)$ as computed from our
theory (16) to those from the DSMC simulations while floating $\eta$, such
that it minimizes the relative root-mean-squared error

$\displaystyle\varepsilon(\eta)=\sqrt{\sum_{t}\left(\frac{T_{p}^{\rm
DMSC}(t)-T_{p}^{\rm theory}(t;\eta)}{T_{p}^{\rm DMSC}(t)}\right)^{2}}.$ (26)

In these numerical experiments, we tune the trap anisotropy in a manner that
does not the affect ${\rm Kn}$, by setting
$\omega_{\perp}=\overline{\omega}/\lambda^{1/6}$ and
$\omega_{z}=\overline{\omega}\lambda^{1/3}$. This construction ensures that
$\overline{\omega}$, and therefore ${\rm Kn}$, both remain independent of
$\lambda$. The dipoles are taken to point along $\hat{\bm{x}}$ for the data
shown. Dependence on dipole orientation will be included below.

Results of several such fits are shown in Fig. 6, which compares the $T_{p}$
time traces for a series of cross-dimensional rethermalization experiments
with $N=5\times 10^{5}$ (${\rm Kn}\approx 0.04$) over a range of
$\lambda=0.13$ to $8.0$, as obtained from DSMC simulations (solid black
curves) and our fitted theory (dashed red curves). Noticeably, there is a
clear beating of various modes with different frequencies which our theory is
able to describe, showing favorable agreement in both the amplitude and phase
of oscillations. A representative comparison plot of ${\cal T}_{r}(t)$ as
obtained from DSMC and Eq. (16) is also provided in Fig. 7, with $N=5\times
10^{5}$ (${\rm Kn}\approx 0.04$) and $\lambda=0.32$. Good agreement is seen in
all ${\cal T}_{r_{i}}(t)$ time traces as well. We note that temperature time
traces tend to show better agreement to the DSMC ones for excitation along the
long axis of a prolate trap, even for larger Knudsen numbers (${\rm Kn}\approx
0.1$). So, we stick to this excitation geometry for a more focused study.

For a given orientation of the dipoles, it may be expected that $\eta$ depends
on both the trap aspect ratio $\lambda$ and the number of molecules $N$.
Increasing $N$, ceteris paribus, evidently increases the density and hence
likely the hydrodynamic volume. As for aspect ratio, a tentative $\lambda$
dependence of ${\cal V}_{\rm hy}$ is already taken into account by (IV.1),
whereby the scaling parameter $\eta$ may depend only weakly on $\lambda$. This
hypothesis is supported by the numerics as shown in Fig. 8, where we find that
$\eta$ is linearly dependent on $N$, but largely independent of $\lambda$ for
the range of these parameters we explore.

Finally, for a given $\lambda$ and $N$, it remains to resolve the dependence
of $\eta$ on the dipole orientation $\hat{\bm{{\cal E}}}$. In this context,
recall that the dilute and hydrodynamic regimes are distinguished by the
Knudsen number, which is inversely proportional to the collision cross
section, Eq. (3). We saw in Sec. IV.1, that this cross section results in
anisotropic viscosities, that work to bring local thermodynamic fluctuations
back to equilibrium. Having accounted for this aspect of differential
scattering, we posit that $\eta$ should only depend on the post-collision
averaged cross section $\sigma_{\rm coll}=\int
d\Omega^{\prime}\frac{d\sigma}{d\Omega^{\prime}}$, which still preserves an
incoming-collision angle dependence Bohn and Jin (2014). As to how so, we
present the following argument. Prolate traps have a weak trapping axis $z$,
along which the gas has a larger thermal width. As a result, the mean-free
path along that axis is relatively smaller compared to the sample size, and
consequently more hydrodynamic. Collisions that occur with relative momentum
directed along the long axis, are then most able to keep molecules behaving
collectively as hydrodynamic. The bulk total cross section is, therefore, most
simply taken as

$\displaystyle\sigma_{\rm
coll}=a_{d}^{2}\frac{\pi}{3}\big{[}3+18\cos^{2}(\hat{\bm{{\cal
E}}}\cdot\hat{\bm{e}}_{\rm hy})-13\cos^{4}(\hat{\bm{{\cal
E}}}\cdot\hat{\bm{e}}_{\rm hy})\big{]},$ (27)

where $\hat{\bm{e}}_{\rm hy}=\hat{\bm{z}}$ denotes the most hydrodynamic axis
(weakest trap frequency), so that $\hat{\bm{{\cal E}}}\cdot\hat{\bm{e}}_{\rm
hy}=\Theta$.

Figure 6:  Comparison of the momentum space temperature $T_{p}$ (IV.1) vs time
$t$, obtained from DSMC simulations (black solid curves) and our theory (red
dashed curves) with $N=5\times 10^{5}$ (${\rm Kn}\approx 0.04$) and parameters
in Tab. 1. The subplots (a) to (h) correspond to various values of trapping
anisotropy with $\lambda=0.13$ to $8.0$ as labeled in the subplot headers. The
fitted values of $\eta$ are also provided in the subplot headers with their
fitting standard uncertainties.  Figure 7:  Comparison of the position space
pseudotemperatures ${\cal T}_{r}$ vs time $t$, obtained from DSMC simulations
(upper subplot a) and our theory (lower subplot b) with the parameters in Tab.
1, $N=5\times 10^{5}$ (${\rm Kn}\approx 0.04$) and $\lambda=0.32$.  Figure 8:
Plot of $\eta$ vs $N$ for various values of $\lambda=0.13,0.20,0.32,0.50$, all
of which are prolate (cigar) geometries. Also plotted is a linear function
ansatz in Eq. (28) (gray dashed line), for comparison with data from DSMC
simulations (blue data). Error bars on the DSMC data points denote standard
fit uncertainties.  Figure 9:  Plot of $\eta$ vs $\Theta$ from a cross-
dimensional rethermalization experiment. The data points (points with error
bars) are obtained from DSMC simulations, which is compared to the fitting
function (dashed curves) in Eq. (28). The data is obtained with the parameters
in Tab. 1 and $\lambda=0.2$, for 3 values of $N=4\times 10^{5}$ (black data,
${\rm Kn}\approx 0.06$), $N=3\times 10^{5}$ (gray data, ${\rm Kn}\approx
0.07$) and $N=2\times 10^{5}$ (light gray data, ${\rm Kn}\approx 0.11$)

We indeed find that $\eta$ follows a $\Theta$ dependence very similar to that
of Eq. (27), when comparing $\eta$ as obtained from DSMC experiments, to a
fitting function of the form $(\sigma_{\rm coll}/\overline{\sigma}_{\rm
coll})\alpha+\beta$ in Fig. 9, where $\overline{\sigma}_{\rm
coll}=\int\sigma_{\rm coll}(\hat{\bm{e}}_{\rm hy})d\hat{\bm{e}}_{\rm
hy}={32\pi a_{d}^{2}/15}$ is the angular averaged total cross section. The
observations above motivate the functional form

$\displaystyle\eta\approx
a+b\left(\frac{N}{10^{5}}\right)\left[1+c\left(\frac{\sigma_{\rm
coll}}{\overline{\sigma}_{\rm coll}}\right)\right],$ (28)

for some constants $a,b$ and $c$, which we determine from fits to be $a\approx
2.21\pm 0.017$, $b\approx 0.67\pm 0.020$ and $c\approx 0.26\pm 0.015$. See
App. C for further details. Our functional guess for the hydrodynamic volume
is therefore

$\displaystyle{\cal V}_{\rm
hy}(\lambda,N,\Theta)\approx\frac{4\pi}{3}\left(\frac{6k_{B}T_{0}}{m\omega_{\perp}^{2}}\right)^{3/2}$
(29)
$\displaystyle\quad\quad\quad\times\frac{1}{\sqrt{\lambda}}\left[2.21+0.67\left(1+0.26\frac{\sigma_{\rm
coll}(\Theta)}{\overline{\sigma}_{\rm coll}}\right)\frac{N}{10^{5}}\right].$

This quasi-empirical formula constitutes the main result of the present paper.
Using this parametrization, the equations of motion (16) can be used to
reliably determine the evolution of a hydrodynamic dipolar Fermi gas in a
prolate trap, subject to excitation along the long axis of the trap.

##  V  Discussions and conclusions

A trapped gas transitions to one that is hydrodynamic when the molecular mean-
free path is far exceeded by the extent of its thermal cloud. Collisional
thermalization is then a local and rapid process, for which collective
dynamics becomes likened to that of a fluid. In this work, we have studied the
damping and oscillations of hydrodynamic welter in harmonically confined
dipolar gases, with cross-dimensional rethermalization experiments.

Unlike its dilute counterpart, a hydrodynamic dipolar gas has its distribution
width (second moment) dynamics closely follow the symmetries imposed by the
confining potential. This adherence to the extrinsic trap symmetry arises from
a high frequency of collisions, suppressing the intrinsic dipolar properties
from manifesting on macroscopic scales. But since local thermal equilibration
is not truly instantaneous, dipolar collisions still result in anisotropic
viscous shearing between fluid layers, damping the macroscopic fluid welter.
We have constructed a model to describe such damped weltering dynamics,
presented in Eq. (16). Embedded in this model is a semi-empirical quantity
${\cal V}_{\rm hy}$, which quantifies the hydrodynamic extent of the trapped
gas and its consequence to damping. Through use of numerical experiments, we
obtain a functional form for ${\cal V}_{\rm hy}$ in Eq. (29), expected to work
in the range of $\lambda$, $N$ and $\Theta$ explored here.

Larger Knudsen numbers and trap anisotropies will increase the dilute
fraction, requiring more nuanced treatments of the non-hydrodynamic regions.
Moreover, the approximation made in Sec. IV of threshold dipolar scattering,
may not be adequate in hydrodynamic samples of polar molecular gases.
Threshold scattering requires that the collision energies relative to the
dipole energy are sufficiently low Bohn  _et al._ (2009), but there be high
enough collision rates to remain hydrodynamic, as is detailed in App. B. This
raises issues for Bose gases within the presented formalism, since lowering
the temperature to achieve threshold scattering would result in a significant
condensate fraction. On the other hand, Fermi gases below $T_{F}$ still have
collective excitations well described by classical kinetic theories, if Pauli
blocking effects are included Nikuni and Griffin (1998). Lastly, dipolar mean-
field effects have been ignored, thermal energies being much larger than the
average dipolar mean-field energy per particle Lahaye  _et al._ (2009). All
these considerations, albeit important to current molecular ultracold
experiments, are not within the current scope of this work and will be
considered in future investigations.

###### Acknowledgements.

The authors would like to thank X. Y. Luo, A. Schindewolf and X. Y. Chen for
insightful and motivating discussions on ultracold molecular trapped gases in
the hydrodynamic regime. This work is supported by the National Science
Foundation under Grant Number PHY2110327.

##  Appendix A  Averaging out spatial coordinates

To obtain the spatially averaged equations of motion in Sec. IV.1, we start by
defining a notation for spatially averaged quantities:

$\displaystyle\langle\ldots\rangle=\frac{1}{N}\int
n(\bm{r},t)\left(\ldots\right)d^{3}r.$ (30)

This renders the density averaged equation for $\sigma_{i}(t)$ as

$\displaystyle\frac{\langle r_{i}^{2}T\rangle}{\sigma_{i}^{2}(t)}-\langle
r_{i}\partial_{i}T\rangle$
$\displaystyle=\frac{m}{k_{B}}\left(\frac{\ddot{\sigma}_{i}(t)}{\sigma_{i}(t)}+\omega_{i}^{2}\right)\langle
r_{i}^{2}\rangle$
$\displaystyle\quad-\sum_{j,k,\ell}\frac{\dot{\sigma}_{k}}{\sigma_{\ell}}\delta_{k,\ell}\int\frac{d^{3}r}{Nk_{B}}r_{i}\partial_{j}\mu_{ijk\ell}(T)$
$\displaystyle=\frac{m}{k_{B}}\left(\frac{\ddot{\sigma}_{i}(t)}{\sigma_{i}(t)}+\omega_{i}^{2}\right)\sigma_{i}^{2}(t)$
(31)
$\displaystyle\quad-\sum_{j,k,\ell}\frac{\dot{\sigma}_{k}}{\sigma_{\ell}}\int\frac{d^{3}r}{Nk_{B}}r_{i}\partial_{j}\mu_{ijk\ell}(T)\delta_{k,\ell}.$

As for the temperature balance equation:

$\displaystyle\frac{\partial T(\bm{r},t)}{\partial
t}+\sum_{i}U_{i}\partial_{i}T(\bm{r},t)+\frac{2}{3}\sum_{i}\partial_{i}U_{i}T(\bm{r},t)$
$\displaystyle\quad\quad\quad=\frac{2}{3n(\bm{r},t)k_{B}}\sum_{i,j,k,\ell}(\partial_{j}U_{i})(\partial_{\ell}U_{k})\mu_{ijk\ell}(T)$
$\displaystyle\quad\quad\quad\quad+\frac{2}{3n(\bm{r},t)k_{B}}\sum_{i,j}\partial_{i}\left[\kappa_{ij}\partial_{j}T(\bm{r},t)\right],$
(32)

we first note the relation

$\displaystyle\frac{d\langle T\rangle}{dt}$
$\displaystyle=\int\frac{d^{3}r}{N}\left[n(\bm{r},t)\frac{\partial
T(\bm{r},t)}{\partial t}+T(\bm{r},t)\frac{\partial n(\bm{r},t)}{\partial
t}\right]$ $\displaystyle=\left\langle\frac{\partial T}{\partial
t}\right\rangle+\sum_{i}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\left(\frac{\langle
r_{i}^{2}T\rangle}{\sigma_{i}^{2}(t)}-\langle T\rangle\right),$ (33)

where we utilized the continuity equation. Then multiplying the temperature
balance equation by $n(\bm{r},t)/N$ and integrating over $d^{3}r$ gives

$\displaystyle\frac{d\langle T\rangle}{dt}$
$\displaystyle+\frac{5}{3}\sum_{i}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\langle
T\rangle-\sum_{i}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\left(\frac{\langle
r_{i}^{2}T\rangle}{\sigma_{i}^{2}(t)}-\langle
r_{i}\partial_{i}T\rangle\right)$
$\displaystyle\quad=\frac{2}{3Nk_{B}}\sum_{i,j,k,\ell}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\delta_{i,j}\left(\int
d^{3}r\mu_{ijk\ell}\right)\delta_{k,\ell}\frac{\dot{\sigma}_{\ell}(t)}{\sigma_{\ell}(t)}$
$\displaystyle\quad\quad\quad\quad+\frac{2}{3Nk_{B}}\sum_{i,j}\int
d^{3}r\left[\partial_{i}(\kappa_{ij}\partial_{j}T)\right].$ (34)

Combining equations (A) and (A), we get

$\displaystyle\frac{d\langle T\rangle}{dt}$
$\displaystyle+\frac{5}{3}\sum_{i}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\langle
T\rangle$
$\displaystyle\quad-\frac{m}{k_{B}}\sum_{i}\dot{\sigma}_{i}(t)\left[\ddot{\sigma}_{i}(t)+\omega_{i}^{2}\sigma_{i}(t)\right]$
$\displaystyle\quad\approx\frac{2}{3Nk_{B}}\sum_{i,j,k,\ell}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\delta_{i,j}\left(\int
d^{3}r\mu_{ijk\ell}\right)\delta_{k,\ell}\frac{\dot{\sigma}_{\ell}(t)}{\sigma_{\ell}(t)}$
$\displaystyle\quad\quad-\frac{1}{Nk_{B}}\sum_{i,j,k,\ell}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\left(\int
d^{3}rr_{i}\partial_{j}\mu_{ijk\ell}\right)\delta_{k,\ell}\frac{\dot{\sigma}_{k}}{\sigma_{\ell}}$
$\displaystyle\quad\quad+\frac{2}{3Nk_{B}}\sum_{i,j}\int
d^{3}r\left[\partial_{i}(\kappa_{ij}\partial_{j}T)\right].$ (35)

At this point, conservation of energy has that

$\displaystyle E_{\rm total}$
$\displaystyle=\frac{m}{2}\sum_{i}\left(\omega_{i}^{2}\langle
r_{i}^{2}\rangle+\int\frac{d^{3}rd^{3}v}{N}f(\bm{r},\bm{v},t)v_{i}^{2}\right)$
$\displaystyle=\frac{m}{2}\sum_{i}\left(\omega_{i}^{2}\sigma_{i}^{2}+\int\frac{d^{3}rd^{3}v}{N}f(\bm{r},\bm{v},t)v_{i}^{2}\right),$
(36)

where $E_{\rm total}$ is the total energy of the hydrodynamic system.
Therefore, the relation above along with Eqs. (IV) and (14) motivates the form
for $\langle T\rangle$ as

$\displaystyle\langle T\rangle$ $\displaystyle=\frac{2E_{\rm
total}}{3k_{B}}-\frac{m}{3k_{B}}\sum_{i}\left[\omega_{i}^{2}\sigma_{i}^{2}(t)+\dot{\sigma}_{i}^{2}(t)\right],$
(37)

and its time-derivative

$\displaystyle\frac{d\langle T\rangle}{dt}$
$\displaystyle=-\frac{2m}{3k_{B}}\sum_{i}\left[\omega_{i}^{2}\dot{\sigma}_{i}(t)\sigma_{i}(t)+\ddot{\sigma}_{i}(t)\dot{\sigma}_{i}(t)\right].$
(38)

Plugging these relations into Eq. (A) and assuming each axis can be solved
independently, we obtain

$\displaystyle\dot{\sigma}_{i}(t)\left[\ddot{\sigma}_{i}(t)+\omega_{i}^{2}\sigma_{i}(t)\right]$
$\displaystyle+\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\left[\frac{1}{3}\sum_{j}\left(\omega_{j}^{2}\sigma_{j}^{2}(t)+\dot{\sigma}_{j}^{2}(t)\right)-\frac{2E_{\rm
total}}{3m}\right]$
$\displaystyle\approx\frac{3}{5Nm}\sum_{j,k,\ell}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\left(\int
d^{3}rr_{i}\partial_{j}\mu_{ijk\ell}\right)\delta_{k,\ell}\frac{\dot{\sigma}_{k}}{\sigma_{\ell}}$
$\displaystyle\quad-\frac{2}{5Nm}\sum_{j,k,\ell}\frac{\dot{\sigma}_{i}(t)}{\sigma_{i}(t)}\delta_{i,j}\left(\int
d^{3}r\mu_{ijk\ell}\right)\delta_{k,\ell}\frac{\dot{\sigma}_{\ell}(t)}{\sigma_{\ell}(t)}$
$\displaystyle\quad-\frac{2}{5Nm}\sum_{j}\int
d^{3}r\left[\partial_{i}(\kappa_{ij}\partial_{j}T)\right].$ (39)

Finally, the conserved total energy $E_{\rm total}$, is made up of the
potential energy and thermal equilibrium temperature $T_{0}$:

$\displaystyle E_{\rm total}$
$\displaystyle=\frac{3}{2}k_{B}T_{0}+\frac{m}{2}\sum_{i}\omega_{i}^{2}\sigma_{0,i}^{2}=3k_{B}T_{0},$
(40)

where we utilized that $\sigma_{0,i}=\sqrt{{k_{B}T_{0}/m\omega_{i}^{2}}}$.

##  Appendix B  Considerations for threshold scattering

The analytic results obtained for the viscosities in Sec. IV.1 are applicable
for close to threshold dipolar scattering, which is energy independent Bohn
and Jin (2014). However this assumption is only appropriate when the collision
energy is much smaller than the characteristic dipole energy $E_{\rm
dd}=16\pi^{2}\epsilon_{0}^{2}\hbar^{6}/m^{3}d^{4}$, where $d$ is the electric
dipole moment Bohn  _et al._ (2009). At the same time, the transport
coefficients are derived with classical kinetic theory that assumes a
nondegenerate sample. Implicit in this formulation is, therefore, that the gas
temperature remains well above the Fermi temperature
$T_{F}=\hbar\overline{\omega}(6N)^{1/3}/k_{B}$ Butts and Rokhsar (1997). The
applicability of our current theory requires that temperature lies in the
range $T_{F}<T\ll E_{\rm dd}/k_{B}$.

Furthermore, the derivation above relies on the gas being hydrodynamic, as is
characterized by the Knudsen number ${\rm Kn}$. The requirements to remain in
the regime of validity as formulated in Sec. IV.1 are summarized as

$\displaystyle\frac{\hbar^{2}}{4ma_{d}^{2}}$ $\displaystyle\gg
k_{B}T>\hbar\overline{\omega}(6N)^{1/3},$ (41a) $\displaystyle N$
$\displaystyle\gg\frac{15\sqrt{\pi}}{4}\frac{k_{B}T}{m\overline{\omega}^{2}a_{d}^{2}},$
(41b)

which is only ever possible if $a_{d}/a_{\rm HO}\ll 0.04$, where
$a_{d}=md^{2}/(8\pi\epsilon_{0}\hbar^{2})$ is the dipole length and $a_{\rm
HO}=\sqrt{\hbar/(m\overline{\omega})}$ is the harmonic oscillator length. In
heteronuclear alkali dimers, these microwave shielded molecules with $d\sim 1$
D and $m\sim 50$ amu have dipole lengths on the order of $a_{d}\sim 5000a_{0}$
to $10,000a_{0}$, in units of Bohr radius $a_{0}$. The necessary trap
frequencies to permit threshold scattering above $T_{F}$ would thus need to be
of order $\omega\ll 10$ Hz, which is very weak compared to typical ultracold
experiments.

For the parameters in Tab. 1, we find that $k_{B}T/E_{\rm dd}\approx 28$,
implying a more accurate cross section would be that obtained from the semi-
classical Eikonal approximation Glauber (1955, 1959); Bohn  _et al._ (2009).
We opt to proceed with the effective cross section obtained with threshold
energy scattering as it still serves to illustrates the effectiveness of our
theory, as formulated for arbitrary cross sections.

##  Appendix C  A simple functional form for the hydrodynamic volume

From Fig. 8, we saw that $\eta$ is mostly independent of $\lambda$, which
leaves us with $\eta=\eta(N,\Theta)$. Then assuming that $\eta$ is separable
in its 2 arguments, this allows us to write
$\eta(N,\Theta)=\eta_{N}(N)\eta_{\Theta}(\Theta)$. Within the range of $N$ we
explore, we could Taylor expand $\eta_{N}$ around a number of molecules that
is sure to be hydrodynamic $N_{0}$, so that

$\displaystyle\eta(N,\Theta)$
$\displaystyle\approx\left(\eta_{N}(N_{0})+(N-N_{0})\left.\frac{\partial\eta_{N}}{\partial
N}\right|_{N_{0}}\right)\eta_{\Theta}(\Theta).$ (42)

Then also assuming that the dependence of $\eta_{\Theta}$ on $\Theta$ purely
arises through $\sigma_{\rm coll}(\Theta)$ (i.e.
$\eta_{\Theta}=\eta_{\Theta}(\sigma_{\rm coll})$), we then treat
$\xi=\sigma_{\rm coll}/\overline{\sigma}_{\rm coll}$ as a small parameters and
Taylor expand $\eta_{\Theta}$ to give

$\displaystyle\eta(N,\Theta)$ $\displaystyle\approx
a+b\left(\frac{N}{10^{5}}\right)\left[1+c\left(\frac{\sigma_{\rm
coll}(\Theta)}{\overline{\sigma}_{\rm coll}}\right)\right],$ (43)

as in Eq. (28), where

$\displaystyle a$
$\displaystyle=\eta_{\Theta}(0)\left(\eta_{N}(N_{0})-N_{0}\left.\frac{\partial\eta_{N}}{\partial
N}\right|_{N_{0}}\right),$ (44a) $\displaystyle b$
$\displaystyle=10^{5}\times\eta_{\Theta}(0)\left.\frac{\partial\eta_{N}}{\partial
N}\right|_{N_{0}},$ (44b) $\displaystyle c$
$\displaystyle=\frac{1}{\eta_{\Theta}(0)}\left.\frac{\partial\eta_{\Theta}}{\partial\xi}\right|_{\xi=0},$
(44c)

having used the notation $\eta_{\Theta}(0)=\eta_{\Theta}(\xi=0)$.

## References

  * Karman and Hutson (2018) T. Karman and J. M. Hutson, Phys. Rev. Lett. 121, 163401 (2018).
  * Karman (2020) T. Karman, Phys. Rev. A 101, 042702 (2020).
  * Anderegg  _et al._ (2021) L. Anderegg, S. Burchesky, Y. Bao, S. S. Yu, T. Karman, E. Chae, K.-K. Ni, W. Ketterle, and J. M. Doyle, Science 373, 779 (2021).
  * Bigagli  _et al._ (2023) N. Bigagli, C. Warner, W. Yuan, S. Zhang, I. Stevenson, T. Karman, and S. Will, Collisionally stable gas of bosonic dipolar ground state molecules (2023), arXiv:2303.16845 [cond-mat.quant-gas] .
  * Quéméner and Bohn (2016) G. Quéméner and J. L. Bohn, Phys. Rev. A 93, 012704 (2016).
  * González-Martínez  _et al._ (2017) M. L. González-Martínez, J. L. Bohn, and G. Quéméner, Phys. Rev. A 96, 032718 (2017).
  * Matsuda  _et al._ (2020) K. Matsuda, L. D. Marco, J.-R. Li, W. G. Tobias, G. Valtolina, G. Quéméner, and J. Ye, Science 370, 1324 (2020).
  * Li  _et al._ (2021) J.-R. Li, W. G. Tobias, K. Matsuda, C. Miller, G. Valtolina, L. De Marco, R. R. Wang, L. Lassablière, G. Quéméner, J. L. Bohn, _et al._ , Nature Physics 17, 1144 (2021).
  * Lassablière and Quéméner (2022) L. Lassablière and G. Quéméner, Phys. Rev. A 106, 033311 (2022).
  * Lin  _et al._ (2023) J. Lin, G. Chen, M. Jin, Z. Shi, F. Deng, W. Zhang, G. Quéméner, T. Shi, S. Yi, and D. Wang, Microwave shielding of bosonic narb molecules (2023), arXiv:2304.08312 [physics.atom-ph] .
  * Valtolina  _et al._ (2020) G. Valtolina, K. Matsuda, W. G. Tobias, J.-R. Li, L. De Marco, and J. Ye, Nature 588, 239 (2020).
  * Schindewolf  _et al._ (2022) A. Schindewolf, R. Bause, X.-Y. Chen, M. Duda, T. Karman, I. Bloch, and X.-Y. Luo, Nature 607, 677 (2022).
  * Baranov (2008) M. Baranov, Physics Reports 464, 71 (2008).
  * Lahaye  _et al._ (2009) T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Reports on Progress in Physics 72, 126401 (2009).
  * Chomaz  _et al._ (2022) L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Reports on Progress in Physics (2022).
  * Bohn and Jin (2014) J. L. Bohn and D. S. Jin, Phys. Rev. A 89, 022702 (2014).
  * Aikawa  _et al._ (2014) K. Aikawa, A. Frisch, M. Mark, S. Baier, R. Grimm, J. L. Bohn, D. S. Jin, G. M. Bruun, and F. Ferlaino, Phys. Rev. Lett. 113, 263201 (2014).
  * Tang  _et al._ (2016) Y. Tang, A. G. Sykes, N. Q. Burdick, J. M. DiSciacca, D. S. Petrov, and B. L. Lev, Phys. Rev. Lett. 117, 155301 (2016).
  * Sykes and Bohn (2015) A. G. Sykes and J. L. Bohn, Phys. Rev. A 91, 013625 (2015).
  * Wang  _et al._ (2020) R. R. W. Wang, A. G. Sykes, and J. L. Bohn, Phys. Rev. A 102, 033336 (2020).
  * Wang and Bohn (2021) R. R. W. Wang and J. L. Bohn, Phys. Rev. A 103, 063320 (2021).
  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * DeMarco and Jin (1999) B. DeMarco and D. S. Jin, Science 285, 1703 (1999).
  * Thomas  _et al._ (2003) A. M. Thomas, S. Hopkins, S. L. Cornish, and C. J. Foot, Journal of Optics B: Quantum and Semiclassical Optics 5, S107 (2003).
  * Griesmaier  _et al._ (2005) A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94, 160401 (2005).
  * Lu  _et al._ (2011) M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011).
  * Aikawa  _et al._ (2012) K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012).
  * Marco  _et al._ (2019) L. D. Marco, G. Valtolina, K. Matsuda, W. G. Tobias, J. P. Covey, and J. Ye, Science 363, 853 (2019).
  * Monroe  _et al._ (1993) C. R. Monroe, E. A. Cornell, C. A. Sackett, C. J. Myatt, and C. E. Wieman, Phys. Rev. Lett. 70, 414 (1993).
  * Newbury  _et al._ (1995) N. R. Newbury, C. J. Myatt, and C. E. Wieman, Phys. Rev. A 51, R2680 (1995).
  * Schmidt  _et al._ (2003) P. O. Schmidt, S. Hensler, J. Werner, A. Griesmaier, A. Görlitz, T. Pfau, and A. Simoni, Phys. Rev. Lett. 91, 193201 (2003).
  * Tang  _et al._ (2015) Y. Tang, A. Sykes, N. Q. Burdick, J. L. Bohn, and B. L. Lev, Phys. Rev. A 92, 022703 (2015).
  * Patscheider  _et al._ (2022) A. Patscheider, L. Chomaz, G. Natale, D. Petter, M. J. Mark, S. Baier, B. Yang, R. R. W. Wang, J. L. Bohn, and F. Ferlaino, Phys. Rev. A 105, 063307 (2022).
  * Guéry-Odelin  _et al._ (1999) D. Guéry-Odelin, F. Zambelli, J. Dalibard, and S. Stringari, Phys. Rev. A 60, 4851 (1999).
  * Ma  _et al._ (2003) Z.-Y. Ma, A. M. Thomas, C. J. Foot, and S. L. Cornish, Journal of Physics B: Atomic, Molecular and Optical Physics 36, 3533 (2003).
  * Note (1) Thus in a tense moment for Ahab aboard a whale boat, Melville relates that the “… White Whale darted through the weltering ocean.” (Moby Dick, Chap. 135.).
  * Schäfer (2014) T. Schäfer, Phys. Rev. A 90, 043633 (2014).
  * Bluhm and Schäfer (2015) M. Bluhm and T. Schäfer, Phys. Rev. A 92, 043602 (2015).
  * Sadeghpour  _et al._ (2000) H. R. Sadeghpour, J. L. Bohn, M. J. Cavagnero, B. D. Esry, I. I. Fabrikant, J. H. Macek, and A. R. P. Rau, Journal of Physics B: Atomic, Molecular and Optical Physics 33, R93 (2000).
  * Huang (1963) K. Huang,  _Statistical mechanics, john wily & sons_ (John Wiley & Sons, Inc New York, 1963) p. 10.
  * Bird (1970) G. A. Bird, The Physics of Fluids 13, 2676 (1970), https://aip.scitation.org/doi/pdf/10.1063/1.1692849 .
  * Franklin and Akman (1985) W. R. Franklin and V. Akman, in  _Computer-Generated Images_ , edited by N. Magnenat-Thalmann and D. Thalmann (Springer Japan, Tokyo, 1985) pp. 176–185.
  * Butts and Rokhsar (1997) D. A. Butts and D. S. Rokhsar, Phys. Rev. A 55, 4346 (1997).
  * Griffin  _et al._ (1997) A. Griffin, W.-C. Wu, and S. Stringari, Phys. Rev. Lett. 78, 1838 (1997).
  * Nikuni and Griffin (1998) T. Nikuni and A. Griffin, Journal of Low Temperature Physics 111, 793 (1998).
  * Kavoulakis  _et al._ (1998) G. M. Kavoulakis, C. J. Pethick, and H. Smith, Phys. Rev. A 57, 2938 (1998).
  * De Groot and Mazur (2013) S. De Groot and P. Mazur,  _Non-Equilibrium Thermodynamics_, Dover Books on Physics (Dover Publications, New York, 2013).
  * Fetter and Walecka (2003) A. Fetter and J. Walecka,  _Theoretical Mechanics of Particles and Continua_, Dover Books on Physics (Dover Publications, 2003).
  * Chapman and Cowling (1990) S. Chapman and T. G. Cowling,  _The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_ (Cambridge university press, New York, 1990).
  * Wang and Bohn (2022a) R. R. W. Wang and J. L. Bohn, Phys. Rev. A 106, 023319 (2022a).
  * Wang and Bohn (2022b) R. R. W. Wang and J. L. Bohn, Phys. Rev. A 106, 053307 (2022b).
  * Wang and Bohn (2023) R. R. W. Wang and J. L. Bohn, Phys. Rev. A 107, 033321 (2023).
  * Reif (2009) F. Reif,  _Fundamentals of statistical and thermal physics_ (Waveland Press, Long Grove, Illinois, 2009).
  * Note (2) With the parameters in Tab. 1, the typical collision energies are in fact much larger than the dipole energy, which may place the assumption of threshold scattering on shaky grounds. Nevertheless, we proceed with the threshold approximation in this work. See App. B for further discussions.
  * Note (3) Our method has similarities to the scaling ansatz employed in Ref. Kagan  _et al._ (1997); Pedri  _et al._ (2003); Wächtler  _et al._ (2017), but formulated slightly differently and with the inclusion of transport tensors.
  * Kadoya  _et al._ (1985) K. Kadoya, N. Matsunaga, and A. Nagashima, Journal of Physical and Chemical Reference Data 14, 947 (1985), https://pubs.aip.org/aip/jpr/article-pdf/14/4/947/11457831/947_1_online.pdf .
  * Note (4) Refs. Kavoulakis  _et al._ (1998); Schäfer (2007) provides another estimate for the effective hydrodynamic volume which, however, we find systematically underestimates the results from our numerical Monte Carlo simulations.
  * Note (5) Incorporating time-dependence in the temperature $T(t)$, requires more sophisticated treatments such as second-order hydrodynamics Lewis and Romatschke (2017). We leave this task to a subsequent study.
  * Bohn  _et al._ (2009) J. L. Bohn, M. Cavagnero, and C. Ticknor, New Journal of Physics 11, 055039 (2009).
  * Glauber (1955) R. J. Glauber, Phys. Rev. 100, 242 (1955).
  * Glauber (1959) R. Glauber, Interscience, New York 1, 315 (1959).
  * Schäfer (2007) T. Schäfer, Phys. Rev. A 76, 063618 (2007).
  * Bluhm and Schäfer (2016) M. Bluhm and T. Schäfer, Phys. Rev. Lett. 116, 115301 (2016).
  * Lewis and Romatschke (2017) W. E. Lewis and P. Romatschke, New Journal of Physics 19, 023042 (2017).
  * Kagan  _et al._ (1997) Y. Kagan, E. L. Surkov, and G. V. Shlyapnikov, Phys. Rev. A 55, R18 (1997).
  * Pedri  _et al._ (2003) P. Pedri, D. Guéry-Odelin, and S. Stringari, Phys. Rev. A 68, 043608 (2003).
  * Wächtler  _et al._ (2017) F. Wächtler, A. R. P. Lima, and A. Pelster, Phys. Rev. A 96, 043608 (2017).