# Simulating Chiral Magnetic and Separation Effects with Spin-Orbit Coupled
Atomic Gases

Xu-Guang Huang  Physics Department and Center for Particle Physics and Field
Theory, Fudan University, Shanghai 200433, China.

###### Abstract

The chiral magnetic and chiral separation effects—quantum-anomaly-induced
electric current and chiral current along an external magnetic field in
parity-odd quark-gluon plasma—have received intense studies in the community
of heavy-ion collision physics. We show that analogous effects occur in
rotating trapped Fermi gases with Weyl-Zeeman spin-orbit coupling where the
rotation plays the role of an external magnetic field. These effects can
induce a mass quadrupole in the atomic cloud along the rotation axis which may
be tested in future experiments. Our results suggest that the spin-orbit
coupled atomic gases are potential simulators of the chiral magnetic and
separation effects.

The recent experimental breakthroughs in generating synthetic spin-orbit
coupling (SOC) in both bosonic Lin:2011  and fermionic gases Wang:2012 ;
Cheuk:2012  have opened a new era for cold atomic physics. In these
experiments, a pair of Raman lasers induced an equal mixture of Rashba and
Dresselhaus SOCs between the (psuedo)spin-$1/2$ internal degree of freedom and
the orbital motion of the atoms. Very promisingly, other types of SOC, e.g.,
the Weyl SOC Anderson:2012 ; Anderson:2013 ; Li:2012 , could also be realized.
The presence of the SOC modifies the dynamics on both single-particle and
many-body levels and a variety of novel properties have been explored.
Furthermore, it is very suggestive that the spin-orbit coupled atomic gases
may provide ideal platforms to simulate intriguing phenomena that have
topological origins, e.g., the topological insulators or superfluid and
Majorana fermions Jiang:2011 ; Liu:2012 ; Zhang:2013 ; Ruhman:2015 , the spin
Hall effect Beeler:2013 ; Kennedy:2013 , and the Berezinskii-Kosterlitz-
Thouless transition He:2011 ; Xu:2015 . See Refs. review:2013 ; review:2014 ;
review:2015  for reviews.

In this article, we demonstrate that yet another topological phenomenon, the
quantum anomaly, can also be realized in a special setup for the spin-orbit
coupled atomic gases, namely the trapped rotating atomic gases with Weyl-
Zeeman SOC. As a consequence of this quantum anomaly, the currents of opposite
chiralities (see below for definition) are generated in parallel or anti-
parallel to the rotation axis. These currents mimic the chiral magnetic effect
(CME) Kharzeev:2008 ; Fukushima:2008  and chiral separation effect (CSE)
Son:2004 ; Met:2005  that are intensively studied in the context of quark-
gluon plasma (QGP) produced in heavy-ion collisions, with now the rotation
playing the role of an external magnetic field.

The QGP version of the CME and CSE is expressed as

$\displaystyle{\bf j}_{V}=\frac{e_{f}N_{c}\mu_{A}}{2\pi^{2}\hbar^{2}}{\bf
B},\;\;\;\;{\bf j}_{A}=\frac{e_{f}N_{c}\mu}{2\pi^{2}\hbar^{2}}{\bf B},$ (1)

for each flavor of light quarks, where $e_{f}$ is the electric charge of quark
with flavor $f$, ${\bf j}_{V}=\langle{\bar{\psi}}{\bm{\gamma}}\psi\rangle$ and
${\bf j}_{A}=\langle{\bar{\psi}}{\bm{\gamma}}\gamma_{5}\psi\rangle$ are
electric and chiral currents, $N_{c}=3$ is the color degeneracy, and $\mu$ and
$\mu_{A}$ are vector and chiral chemical potentials. Experimentally, signals
consistent with CME and CSE have been observed in heavy-ion collisions at
Relativistic Heavy Ion Collider (RHIC) STAR  and Large Hadron Collider (LHC)
ALICE . In these collisions, extremely strong magnetic fields arise due to the
fast motion of the ions Skokov ; Huang ; Huang2015 , and these magnetic fields
induce charge separation and chirality separation in the QGP via CME and CSE
which in turn lead to special azimuthal distributions of the charged hadrons
that are finally measured by the detectors.

It is worth noting that the CME was also discussed in astrophysical context
Vilenkin:1980  and more recently in Weyl and Dirac semimetals Burkov:2012 ;
Vaz:2013 ; Basar:2014 ; Land:2014 ; weyl:2014 . In Weyl and Dirac semimetals,
the low energy excitations are Weyl and Dirac fermions. These materials can
exhibit finite $\mu_{5}$ due to their special band structure or by applying
parallel electric and magnetic fields, and thus open the possibility of
realizing the CME. Comparing to these previously explored systems, the atomic
gases with Weyl-Zeeman SOC enable greater flexibility in controlling the
parameters and thus provide not only simulators but also a unique mean to
exploit new features (e.g., those raised by the presence of the Zeeman
splitting field) of the CME and CSE.

Semiclassical approach.— We begin by considering the following single-particle
Hamiltonian for spin-$1/2$ atoms (either bosons or fermions) in three
dimensions (3D),

$\displaystyle{\cal H}=\frac{{\mathbf{p}}^{2}}{2m}-\mu-{\bf
W}({\mathbf{p}})\cdot{\bm{\sigma}},$ (2)

where $m$ is the mass, $\mu$ is the chemical potential, and ${\mathbf{p}}$ is
the canonical momentum. The third term expresses a generic SOC where
${\bm{\sigma}}$ is the Pauli matrix and ${\bf W}$ represents a momentum-
dependent magnetic field. Its form will be given when necessary. The
Hamiltonian (2) represents a two-band system with band dispersions
$\varepsilon_{c}({\mathbf{p}})=p^{2}/(2m)-\mu-cW({\mathbf{p}}),c=\pm$.
Correspondingly, we define the chirality of band $c$ to be right-handed (left-
handed) if $c=+$ ($c=-$). We note that the chirality we defined here is
commonly called helicity in literature. We choose the term “chirality” to keep
the consistence with the terminologies “chiral anomaly”, “chiral magentic
effect”, etc.

Now consider that the atoms are trapped by an external potential
$V({\mathbf{x}})$ and at the meantime are subject to rotation with angular
velocity ${\bm{\omega}}$ which we assume to be a constant. The effect of the
trapping and rotation can be described by a gauge potential
$A^{\mu}=(A_{0},{\bf A})$ with
$A_{0}({\mathbf{x}})=V({\mathbf{x}})-(m/2)({\bm{\omega}}\times{\mathbf{x}})^{2}-\mu$
and ${\bf A}({\mathbf{x}})=m{\bm{\omega}}\times{\mathbf{x}}$, and then the
Hamiltonian becomes

$\displaystyle{\cal H}=\frac{[{\mathbf{p}}-{\bf
A}({\mathbf{x}})]^{2}}{2m}-{\bf W}[{\mathbf{p}}-{\bf
A}({\mathbf{x}})]\cdot{\bm{\sigma}}+A_{0}({\mathbf{x}}).$ (3)

Note that the Hamiltonian of atoms with SOC in the rotating frame is generally
time dependent and the minimal substitution adoptted in going from Eq. (2) to
Eq. (3) may not be applicable. However, for situations that, e.g., both the
trap and the lasers inducing the SOC are rotating and the SOC is linear in
momentum, one can find a time-dependent unitary transformation to eliminate
the time dependence from the Hamiltonian; see Ref. timedep  and the references
therein. Nevertheless, the present study is restricted to such a situation and
the use of the time-independent Halmitonian (3) is justified.

We now derive a set of semiclassical equations of motion (EOM) for the orbital
variables ${\mathbf{x}}$ and ${\mathbf{p}}$. (For non-rotating spin-orbit
coupled atomic gases, a similar derivation is given in Ref. bijl .) At
semiclassical level, we treat ${\mathbf{x}},{\mathbf{p}}$, and ${\bm{\sigma}}$
as classical variables, and their EOM are easily derived from Hamiltonian (3):

$\displaystyle\dot{{\mathbf{x}}}$ $\displaystyle=$
$\displaystyle\nabla_{\mathbf{p}}\xi-\nabla_{\mathbf{p}}{\bf
W}\cdot{\bm{\sigma}},$ (4) $\displaystyle\dot{{\mathbf{p}}}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbf{x}}\xi+\nabla_{\mathbf{x}}{\bf
W}\cdot{\bm{\sigma}},$ (5) $\displaystyle\dot{{\bm{\sigma}}}$ $\displaystyle=$
$\displaystyle\frac{2}{\hbar}{\bm{\sigma}}\times{\bf W},$ (6)

where $\xi({\mathbf{x}},{\mathbf{p}})=[{\mathbf{p}}-{\bf
A}({\mathbf{x}})]^{2}/(2m)+A_{0}({\mathbf{x}})$ and ${\bf
W}({\mathbf{x}},{\mathbf{p}})={\bf W}[{\mathbf{p}}-{\bf A}({\mathbf{x}})]$. To
proceed, we make an adiabatic approximation to the spin dynamics, that is, we
treat the orbital degrees of freedom ${\mathbf{x}},{\mathbf{p}}$ as slow
variables while the spin ${\bm{\sigma}}$ as fast variable and solve Eq. (6) up
to first order in time derivatives of the orbital variables. This gives

$\displaystyle{\bm{\sigma}}\approx c{\bf w}+c\frac{\hbar}{2W}{\bf
w}\times\dot{{\bf w}},$ (7)

where ${\bf w}={\bf W}/W$. This procedure is essentially equivalent to solving
${\bm{\sigma}}$ up to first order in $\hbar$. Inserting Eq. (7) to Eqs.
(4)-(5) we obtain

$\displaystyle\dot{{\mathbf{x}}}$ $\displaystyle=$
$\displaystyle\nabla_{\mathbf{p}}(\xi-
cW)+c\hbar{\bm{\Omega}}_{px}\cdot\dot{{\mathbf{x}}}+c\hbar{\bm{\Omega}}_{pp}\cdot\dot{{\mathbf{p}}},$
(8) $\displaystyle\dot{{\mathbf{p}}}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbf{x}}(\xi-
cW)-c\hbar{\bm{\Omega}}_{xx}\cdot\dot{{\mathbf{x}}}-c\hbar{\bm{\Omega}}_{xp}\cdot\dot{{\mathbf{p}}},$
(9)

where the tensors $\Omega_{AB}^{ij}$ with $A,B=x,p$ are defined as

$\displaystyle\Omega_{AB}^{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{\partial{\bf w}}{\partial
A_{i}}\times\frac{\partial{\bf w}}{\partial B_{j}}\right)\cdot{\bf w}.$ (10)

In terms of the kinetic momentum, ${\mathbf{k}}\equiv{\mathbf{p}}-{\bf
A}({\mathbf{x}})$, Eqs. (8)-(9) can be recast to more compact and
transparently gauge invariant forms Sundaram ; Xiao:2010 ,

$\displaystyle\dot{{\mathbf{x}}}$ $\displaystyle=$
$\displaystyle\nabla_{\mathbf{k}}\varepsilon_{c}+c\hbar\dot{{\mathbf{k}}}\times{\bm{\Omega}},$
(11) $\displaystyle\dot{{\mathbf{k}}}$ $\displaystyle=$ $\displaystyle{\bf
E}+\dot{{\mathbf{x}}}\times{\bf B},$ (12)

where $\varepsilon_{c}({\mathbf{k}})=k^{2}/(2m)-\mu-cW({\mathbf{k}})$ and

$\displaystyle E_{i}$ $\displaystyle=$ $\displaystyle-\frac{\partial
A_{0}}{\partial x_{i}},$ (13) $\displaystyle B_{i}$ $\displaystyle=$
$\displaystyle\epsilon_{ijk}\frac{\partial A_{k}}{\partial
x_{j}}=2m\omega_{i},$ (14) $\displaystyle\Omega_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{ijk}\Omega_{pp}^{jk},$ (15)

are the effective electric field, magnetic field, and Berry curvature,
respectively. Equation (14) exhibits the equivalence between ${\bf B}$ and
${\bm{\omega}}$. Note that for relativistic massless particles under rotation,
the effective magnetic field would be momentum-dependent Stephanov , ${\bf
B}\sim 2|{\mathbf{k}}|{\bm{\omega}}$, and thus the rotation is no longer
equivalent to a magnetic field. In that case, the rotation can induce an
independent current other than the CME/CSE which is called chiral vortical
effect (CVE) chialve ; chialve2 ; chialve3 ; Son2 ; cvw . In the non-
relativistic case, the CVE is equivalent to CME/CSE.

From Eq. (11) and Eq. (12) we obtain

$\displaystyle\sqrt{G_{c}}\dot{{\mathbf{x}}}$ $\displaystyle=$
$\displaystyle\nabla_{\mathbf{k}}\varepsilon_{c}+c\hbar{\bf
E}\times{\bm{\Omega}}+c\hbar({\bm{\Omega}}\cdot\nabla_{\mathbf{k}}\varepsilon_{c}){\bf
B},$ (16) $\displaystyle\sqrt{G_{c}}\dot{{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\bf E}+\nabla_{\mathbf{k}}\varepsilon_{c}\times{\bf
B}+c\hbar({\bf E}\cdot{\bf B}){\bm{\Omega}},$ (17)

where $\sqrt{G_{c}}=1+c\hbar{\bf B}\cdot{\bm{\Omega}}$, its physical meaning
will be clear soon. These are the semiclassical equations for the orbital
motion of atoms of chirality $c$. In these equations, the quantum effects are
reflected in the Berry curvature terms.

We will hereafter focus on Fermi gases and we will use the natural units
$\hbar=k_{B}=1$.

In the presence of the Berry curvature, the invariant measure of the phase
space integration for atoms of chirality $c$ needs to be modified to
$\sqrt{G_{c}}d^{3}{\mathbf{x}}d^{3}{\mathbf{k}}/(2\pi)^{3}$ Xiao ; dual:2006 .
With this notification, one is able to write down a kinetic equation for the
distribution function $f_{c}$ of chirality $c$,

$\displaystyle\partial_{t}f_{c}+\dot{{\mathbf{x}}}\cdot\nabla_{\mathbf{x}}f_{c}+\dot{{\mathbf{k}}}\cdot\nabla_{\mathbf{k}}f_{c}=I[f_{c}],$
(18)

where $\dot{{\mathbf{x}}}$ and $\dot{{\mathbf{k}}}$ are given by Eqs.
(16)-(17). In the context of relativistic chiral medium like QGP, similar
kinetic equation has been derived recently Stephanov ; Yamamoto ; Wangqun ;
Wangqun2 ; Mannuel . The density and current of chirality $c$ are defined as

$\displaystyle n_{c}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\sqrt{G_{c}}f_{c},$ (19)
$\displaystyle{\bf j}_{c}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\sqrt{G_{c}}\dot{{\mathbf{x}}}_{c}f_{c},$
(20)

and the continuity equation, following Eq. (18), reads

$\displaystyle\partial_{t}n_{c}+\nabla_{\mathbf{x}}\cdot{\bf j}_{c}$
$\displaystyle=$ $\displaystyle c({\bf E}\cdot{\bf
B})\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f_{c}\nabla_{\mathbf{k}}\cdot{\bm{\Omega}}$
(21) $\displaystyle=$ $\displaystyle
cf_{c}({\mathbf{k}}_{0})\frac{F}{4\pi^{2}}{\bf E}\cdot{\bf B},$

where we suppose that the collision kernel conserves the particle number for
each chirality. If not, there would be a term
$\int{d^{3}{\mathbf{k}}}/{(2\pi)^{3}}\sqrt{G_{c}}I[f_{c}]$ on the right-hand
side. The ${\mathbf{k}}_{0}$ in the second line specifies the location of the
Berry monopole which coincides with the band-crossing point determined by
${\bf W}({\mathbf{k}}_{0})={\bf 0}$. The $F$ is the total Berry curvature flux
associated with the Berry monopole. Its explicit expression is

$\displaystyle F$ $\displaystyle=$
$\displaystyle\frac{\epsilon_{ijm}}{8\pi}\int_{S^{2}}d^{2}\Sigma_{i}{\bf
w}\cdot\left(\frac{\partial{\bf w}}{\partial k_{j}}\times\frac{\partial{\bf
w}}{\partial k_{m}}\right).$ (22)

This counts the winding number of the map ${\bf w}:S^{2}\rightarrow S^{2}$ and
thus $F\in\pi_{2}(S^{2})=\mathbb{Z}$.

If the right-hand side does not vanish, Eq. (21) represents a quantum anomaly
for the current of chirality $c$ in the form analogous to the chiral anomaly
in gauge field theory. This can be seen more clearly if we consider a Fermi
gas at zero temperature with pure Weyl SOC, ${\bf W}=\lambda{\mathbf{k}}$,
where $\lambda$ is the strength of the SOC. In this case,
${\bm{\Omega}}={\mathbf{k}}/(2k^{3})$,
$\nabla_{\mathbf{k}}\cdot{\bm{\Omega}}=2\pi\delta^{(3)}({\mathbf{k}})$, $F=1$,
${\mathbf{k}}_{0}={\bf 0}$, and $f_{c}({\mathbf{k}}_{0})=1$. Thus the right-
hand side of Eq. (21) reads $c({\bf E}\cdot{\bf B})/(4\pi^{2})$, which
coincides exactly with the $U(1)$ chiral anomaly. In this case, the
conservation of particle numbers of chirality $c$ which is proportional to the
volume of its corresponding Fermi sphere is violated by the flux of the Berry
curvature across the Fermi surface. For a Weyl-Zeeman SOC, ${\bf
W}({\mathbf{k}})=\lambda{\mathbf{k}}+h\hat{\bf z}$ with $h>0$ being a constant
Zeeman field, the Berry curvature is

$\displaystyle{\bm{\Omega}}$ $\displaystyle=$
$\displaystyle\lambda^{2}\frac{{\bf W}({\mathbf{k}})}{2W^{3}},$ (23)

the Berry monopole locates at ${\mathbf{k}}_{0}=-h/\lambda\hat{\bf z}$, and
the winding number $F=1$. At zero temperature, we have
$f_{c}({\mathbf{k}}_{0})=\theta[\mu-h^{2}/(2m\lambda^{2})]$ which vanishes if
the Zeeman field $h>\lambda\sqrt{2m\mu}$. Thus the system has quantum anomaly
only when $h<\lambda\sqrt{2m\mu}$. The absence of the quantum anomaly when
$h>\lambda\sqrt{2m\mu}$ reflects a change of the Fermi-surface topology as
shown in Fig. 1.

We note here that there is no quantum anomaly for Rashba-Dresselhaus SOC ${\bf
W}({\mathbf{k}})=(\lambda k_{y}-\eta k_{x},-\lambda k_{x}+\eta k_{y},h)$ with
$\lambda,\eta,h>0$ being constants, because its Berry curvature,

$\displaystyle{\bm{\Omega}}$ $\displaystyle=$
$\displaystyle\frac{(\lambda^{2}-\eta^{2})h}{2[h^{2}+(\lambda^{2}+\eta^{2})(k_{x}^{2}+k_{y}^{2})-2\lambda\eta
k_{x}k_{y}]^{3/2}}\hat{\bf z},$ (24)

leads to zero winding number.

Figure 1: (Color online) Upper panels: The dispersion relations
$\varepsilon_{c}({\mathbf{k}}_{\perp}=0,k_{z})$. Lower panels: The Fermi-
surface topologies in $(k_{z},k_{x})$ plane. The green dashed line represents
the chemical potential. Blue (red) lines are for right-(left-)handed fermions.
When $h<\lambda\sqrt{2m\mu}$, the Berry monopole (the black point) locates
inside both the Fermi surfaces; At $h=\lambda\sqrt{2m\mu}$ the two Fermi
surfaces touch; when $h>\lambda\sqrt{2m\mu}$, the Berry monopole moves out
both the Fermi surfaces.

Before we proceed, let us comment on the validity regime of the semiclassical
approach. To arrive at Eqs. (16)-(17) the inter-band transition has been
neglected which means that the force acting on the atoms cannot be strong. In
particular, this require that $\sqrt{|{\bf B}|}\ll W({\mathbf{p}})/\lambda$;
see Refs. Sundaram ; Xiao:2010 ; Stephanov  for more discussions. Obviously,
this condition is violated if ${\mathbf{k}}$ is close to the Berry monopole
where $W({\mathbf{p}})$ is small. Thus the phase space integral in Eq. (21)
should be understood to exclude the region
$|{\mathbf{k}}-{\mathbf{k}}_{0}|<\Delta$ around the Berry monopole with
$\Delta$ large enough so that we can apply the classical description to
particles outside of it. The value of $\Delta$ depends on ${\bf B}$ and ${\bf
E}$. For example, for pure Weyl SOC, $\Delta$ should be larger than
$\sqrt{|{\bf B}|}$; this actually constrains the magnitude of ${\bf B}$:
because $\Delta$ should not exceed $k_{F}$, the maximum $|{\bf B}|$ should not
exceed $k_{F}^{2}$. This implies that the rotation frequency should not exceed
$\varepsilon_{F}$ in order to guarantee the validity of the semiclassical
approach. Our numerical simulations will always be within the validity regime
of the semiclassical approach.

Chiral magnetic and chiral separation effects.— A consequence of this quantum
anomaly is the appearance of CME and CSE. To see this, we substitute Eq. (16)
into Eq. (20),

$\displaystyle{\bf j}_{c}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}f_{c}\nabla_{\mathbf{k}}\varepsilon_{c}+c{\bf
E}\times\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}{\bm{\Omega}}f_{c}$ (25)
$\displaystyle+c{\bf
B}\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}({\bm{\Omega}}\cdot\nabla_{\mathbf{k}}\varepsilon_{c})f_{c}.$

The first term on the right-hand side is the normal number current, the second
term is the (intrinsic) anomalous Hall effect and the last term represents a
${\bf B}$-induced current which we denote by ${\bf j}_{c}^{\rm{\bf B}-ind}$:

$\displaystyle{\bf j}^{\rm{\bf B}-ind}_{c}$ $\displaystyle=$
$\displaystyle\chi_{c}{\bf B},$ (26) $\displaystyle\chi_{c}$ $\displaystyle=$
$\displaystyle
c\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}({\bm{\Omega}}\cdot\nabla_{\mathbf{k}}\varepsilon_{c})f_{c},$
(27)

where $\chi_{c}$ is the ${\bf B}$-induced conductivity (BIC) of chirality $c$.
Let us consider $f_{c}$ to be the Fermi-Dirac distribution,

$\displaystyle f_{c}=f_{c}^{0}=\frac{1}{e^{\beta\varepsilon_{c}}+1}\;\;\;{\rm
for\;fermions}.$ (28)

In this case, one finds that

$\displaystyle\chi_{R}=-\chi_{L}=\frac{T}{4\pi^{2}}\ln\left[1+\exp{(-\beta\varepsilon_{0})}\right]$
(29)

with $\varepsilon_{0}=\mu[h^{2}/(2m\mu\lambda^{2})-1]$. Note that the
$\lambda$ and $h$ dependence of $\chi_{c}$ is through their ratio $h/\lambda$,
and thus when $h=0$ the BIC is independent of $\lambda$ as long as it is
nonzero. We present the numerical results for $\chi_{R}$ in Fig. 2.

Figure 2: (Color online) The ${\bf B}$-induced conductivity, $\chi_{R}$, for
Fermi gases with Weyl-Zeeman SOC as function of the temperature (left panel)
and $h/(\lambda\sqrt{2m\mu})$ (right panel). Note that $\chi_{L}=-\chi_{R}$.

It is seen that the BIC is enhanced at finite temperature and suppressed by
large Zeeman field or small SOC strength. The latter effect is more
transparent at zero temperature at which $f_{c}^{0}=\theta(-\varepsilon_{c})$
and $\chi_{c}$ can be analytically obtained:

$\displaystyle\chi_{c}(T=0)$ $\displaystyle=$
$\displaystyle\frac{c}{4\pi^{2}}\left(\mu-\frac{h^{2}}{2m\lambda^{2}}\right)\theta\left(\mu-\frac{h^{2}}{2m\lambda^{2}}\right).$
(30)

Thus, the ${\bf B}$-induced currents disappear when $h>\lambda\sqrt{2m\mu}$
which again reflects a Fermi-surface topology transition (see Fig. 1), in
parallel to the absence of the quantum anomaly.

The above result $\chi_{R}=-\chi_{L}$ reflects the fact that the right-handed
and left-handed atoms have the same chemical potential. If the atomic cloud
contains domains (a possible realization will be presented in next section) in
which the right-handed and left-handed chemical potentials differ, say,

$\displaystyle\mu_{R}=\mu+\mu_{A},\;\;\mu_{L}=\mu-\mu_{A},$ (31)

the two BICs, $\chi_{R}$ and $\chi_{L}$, in these domains will also differ in
magnitude and Eq. (30) becomes

$\displaystyle\chi_{c}(T=0)$ $\displaystyle=$
$\displaystyle\frac{c}{4\pi^{2}}\left(\mu_{c}-\frac{h^{2}}{2m\lambda^{2}}\right)\theta\left(\mu_{c}-\frac{h^{2}}{2m\lambda^{2}}\right).$
(32)

At zero Zeeman field, by inserting $\chi_{c}(T=0)$ into Eq. (26) we obtain the
${\bf B}$-induced vector and chiral currents:

$\displaystyle{\bf j}_{V}^{\rm{\bf B}-ind}$ $\displaystyle\equiv$
$\displaystyle{\bf j}_{R}^{\rm{\bf B}-ind}+{\bf j}_{L}^{\rm{\bf
B}-ind}=\frac{\mu_{A}}{2\pi^{2}}{\bf B},$ (33) $\displaystyle{\bf
j}_{A}^{\rm{\bf B}-ind}$ $\displaystyle\equiv$ $\displaystyle{\bf
j}_{R}^{\rm{\bf B}-ind}-{\bf j}_{L}^{\rm{\bf B}-ind}=\frac{\mu}{2\pi^{2}}{\bf
B}.$ (34)

These equations express the CME and the CSE in forms consistent with Eq. (1).

Finally, we note that the BICs do not receive perturbative corrections from
scattering Son2 ; Hou ; Banerjee ; Jensen ; Satow  (See the Method.). This
originates from the fact that the chiral anomaly is free of renormalization
(the Adler-Bardeen theorem).

Figure 3: (Color online) A schematic illustration of CME and CSE induced
chiral dipole and mass quadrupole in Fermi gases with Weyl SOC. First, the CSE
drives a chirality separation along the rotation axis with $\mu_{A}>0$ in the
upper tip and $\mu_{A}<0$ in the lower tip (the left panel). Then the CME in
turn drives particle number or equivalently the mass to flow away from the
center and the atomic cloud acquires a mass quadrupole along the rotation axis
(the right panel).

Chiral dipole and mass quadrupole.— Now we turn to the phenomenology of the
CME and CSE. We consider a Weyl-Zeeman spin-orbit coupled Fermi gas in normal
phase with $\mu>0$. Once it is rotating, the CSE (34) will drive the chirality
to move along ${\bf B}$ and causes a macroscopic separation of the chirality,
see Fig. 3 (left). This chiral dipolar distribution naturally forms two
separated domains, one with $\mu_{A}>0$ in the upper space while another with
$\mu_{A}<0$ in the lower space. Then in these domains the CME (33) in turn
drives particle number or equivalently the mass to flow away from the center
and the atomic cloud acquires a mass quadrupole along ${\bf B}$, see Fig. 3
(right). In the context of heavy-ion collisions, similar mechanism was
proposed to generate a charge quadrupole in QGP which may induce a difference
between the elliptic flows of $\pi^{+}$ and $\pi^{-}$ shov ; cmvexp  that has
been detected at RHIC Wangcmv .

The above argument is true only at the qualitative level, to reveal the real
dynamical process, one needs to solve the coupled evolution problem of the
right- and left-handed currents. Let us consider a pure Weyl SOC and the $T=0$
case. The general forms of the right-handed and left-handed currents should
contain diffusion and normal conducting terms, so they read (the anomalous
Hall terms vanish for pure Weyl SOC because of the time-reversal symmetry)

$\displaystyle{\bf j}_{R}$ $\displaystyle=$
$\displaystyle\frac{\mu_{R}}{4\pi^{2}}{\bf B}-D_{R}\nabla n_{R}+\sigma{\bf
E},$ (35) $\displaystyle{\bf j}_{L}$ $\displaystyle=$
$\displaystyle-\frac{\mu_{L}}{4\pi^{2}}{\bf B}-D_{L}\nabla n_{L}+\sigma{\bf
E},$ (36)

where $D_{c}$ is the diffusion constant and $\sigma$ is the normal
conductivity, they are linked by the Einstein relations $\sigma=D_{R}(\partial
n_{R}/\partial\mu_{R})=D_{L}(\partial n_{L}/\partial\mu_{L})$. Consider a
small fluctuation $\delta n_{c}$ in density $n_{c}$. This will induce a small
departure of the chemical potential from the background value,
$\mu_{c}=\mu_{\rm bg}+\delta\mu_{c}$, where $\mu_{\rm bg}=-A_{0}$ is linked to
${\bf E}$ by ${\bf E}=\nabla\mu_{\rm bg}$. Substituting ${\bf j}_{c}$ to Eq.
(21) and keeping linear order terms in $\delta n_{c}$ we obtain

$\displaystyle\frac{\partial\delta n_{R}}{\partial
t}+\frac{1}{4\pi^{2}}\frac{\partial\mu_{R}}{\partial n_{R}}{\bf
B}\cdot\nabla\delta n_{R}-D_{R}\nabla^{2}\delta n_{R}=0,$ (37)
$\displaystyle\frac{\partial\delta n_{L}}{\partial
t}-\frac{1}{4\pi^{2}}\frac{\partial\mu_{L}}{\partial n_{L}}{\bf
B}\cdot\nabla\delta n_{L}-D_{L}\nabla^{2}\delta n_{L}=0.$ (38)

These two equations represent two wave modes with dispersions
$E_{R}({\mathbf{k}})=v_{R}k-iD_{R}k^{2}$,
$E_{L}({\mathbf{k}})=v_{L}k-iD_{L}k^{2}$, one propagating along ${\bf B}$ with
velocity $v_{R}=B({\partial\mu_{R}}/{\partial n_{R}})/(4\pi^{2})$ and another
opposite to ${\bf B}$ with velocity $v_{L}=B({\partial\mu_{L}}/{\partial
n_{L}})/(4\pi^{2})$. We call them chiral magnetic waves (CMWs) in accordance
with the same wave modes found in the context of QGP cmv . It is the CMWs that
develops the chiral dipole and the mass quadrupole. Unlike the situation in
QGP, in our setup the trapping potential will finally balance the driving
force due to the CMWs and establish a new mechanical equilibrium. The new
equilibrium will be characterized by the position-dependent chemical potential

$\displaystyle\mu_{c}({\mathbf{x}})$ $\displaystyle=$
$\displaystyle\mu-\tilde{V}({\mathbf{x}})+\Delta\mu_{c}({\mathbf{x}}),$ (39)

where
$\tilde{V}=(m/2)\left(\tilde{\omega}_{\perp}^{2}(x^{2}+y^{2})+\omega_{z}^{2}z^{2}\right)$
is the effective trapping potential and
$\tilde{\omega}_{\perp}^{2}=\omega_{\perp}^{2}-\omega^{2}$. The chemical
potential shift $\Delta\mu_{c}({\mathbf{x}})$ is determined by the mechanical
equilibration condition ${\bf j}_{R}={\bf j}_{L}=0$. The numerical results are
shown in Fig. 4 where we simulate $5\times 10^{5}$ fermions in the harmonic
trap. Length is in unit of $1/\sqrt{m\omega_{z}}$, and we assume the
transverse effective trapping frequency $\tilde{\omega}_{\perp}$ being kept to
be $\omega_{z}$ when the ${\bf B}$-field changes. Experimentally, the chiral
dipole may be hard to detect but the mass quadrupole profile can be easily
detected by, e.g., light absorption images.

Figure 4: (Color online) Mass quadrupole (upper panels) and chiral dipole
(lower panels) induced by CME and CSE at different rotating frequencies.
Length is in unit of $1/\sqrt{m\omega_{z}}$.

Conclusion.— In summary, we have demonstrated that if the atomic gases with
Weyl-Zeeman SOC is 1) trapped by an external potential and 2) under rotation,
there can appear a quantum anomaly in the chiral currents. A consequence of
this chiral anomaly is the chiral magnetic effect and chiral separation
effect. The CME and CSE cause macroscopic separation of chirality and a mass
quadrupole along the rotation axis in the fermionic atomic cloud which may
possibly be detected in cold atomic experiments.

In the context of QGP, there has been found other transport phenomena that
stem from the quantum anomaly, e.g., the chiral electric separation effect
Huang2 ; Jiang , which may also be realized in spin-orbit coupled atomic gases
in the similar setup as we discussed in this Report. In addition, it is clear
from the derivation that the chiral anomaly may exist also in Bose gases. How
the quantum-anomaly-induced transports in bosonic gases also deserve detailed
exploration in future works.

Method.— Now we show the robustness of ${\bf B}$-induced currents against
scattering. Let us consider the collision kernel $I[f_{c}]$ in the relaxation
time approximation,

$\displaystyle I[f_{c}]=-\frac{\delta f_{c}}{\tau},$ (40)

where $\delta f_{c}=f_{c}-f_{c}^{0}\propto{\bf B}$ (which is assumed to be
small in this calculation) with $f_{c}^{0}$ the equilibrium distribution. We
now show that $\chi_{c}$ is independent of $\tau$. The relaxation time $\tau$
is assumed to be independent of the chirality, and we use it to characterize
the interaction strength among atoms. Concretely, the relaxation time can be
expressed as $\tau=1/(n\bar{v}\sigma_{T})$ with $n$ the density, $\bar{v}$ the
average velocity, $\sigma_{T}=4\pi a^{2}/(1+\bar{k}^{2}a^{2})$ the total cross
section, and $a$ the scattering length. For fermions, at low temperature,
$\tau\sim 3\pi m(1+(k_{F}a)^{2})/(4k_{F}^{2}(k_{F}a)^{2})$. At high
temperature, $\tau\sim m^{3/2}\sqrt{T}/n$.

To proceed, let us turn the ${\bf E}$-field off and assume a steady state with
no spatial dependence of $f_{c}$ and ${\bf B}$. Substituting $f_{c}$ into the
kinetic equation (18) in the main text, at linear order in ${\bf B}$, we
obtain

$\displaystyle\delta
f_{c}\approx-\frac{\tau}{\sqrt{G_{c}}}(\nabla_{\mathbf{k}}\varepsilon_{c}\times{\bf
B})\cdot\nabla_{\mathbf{k}}f_{c}^{0}.$ (41)

The ${\bf B}$-field induced current of chirality $c$ is given by

$\displaystyle{\bf j}_{c}^{{\bf B}-{\rm
ind}}=\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\sqrt{G_{c}}\dot{{\mathbf{x}}}(f_{c}^{0}+\delta
f_{c}),$ (42)

where the first term gives the result (26) in the main text and the second
term

$\displaystyle\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\sqrt{G_{c}}\dot{{\mathbf{x}}}\delta
f_{c}$ (43) $\displaystyle=$
$\displaystyle-\tau\int\frac{d^{3}{\mathbf{k}}}{(2\pi)^{3}}\dot{{\mathbf{x}}}(\nabla_{\mathbf{k}}\varepsilon_{c}\times{\bf
B})\cdot\nabla_{\mathbf{k}}\varepsilon_{c}\frac{\partial
f_{c}^{0}}{\partial\varepsilon_{c}}=0.$

Thus the ${\bf B}$-induced currents are solely given by CME and CSE and are
free of perturbative corrections.

## References

  * (1) Lin, Y.-J., Jimenez-Garcia, K. & Spielman, I. B. Spin-orbit-coupled Bose-Einstein condensates. Nature 471, 83 (2011).
  * (2) Wang, P. _et al._ Spin-orbit coupled degenerate Fermi gases. Phys. Rev. Lett. 109, 095301 (2012).
  * (3) Cheuk, L. W. _et al._ Spin-injection spectroscopy of a spin-orbit coupled Fermi gas. Phys. Rev. Lett. 109, 095302 (2012).
  * (4) Anderson, B. M., Juzeliūnas, G., Galitski, V. M. & Spielman, I. B. Synthetic 3D spin-orbit coupling. Phys. Rev. Lett. 108, 235301 (2012).
  * (5) Anderson, B. M., Spielman, I. B. & Juzeliūnas, G. Magnetically generated spin-orbit coupling for ultracold atoms. Phys. Rev. Lett. 111, 125301 (2013).
  * (6) Li, Y., Zhou, X. & Wu, C. Two- and three-dimensional topological insulators with isotropic and parity-breaking Landau levels. Phys. Rev. B 85, 125122 (2012).
  * (7) Jiang, L. _et al._ Majorana fermions in equilibrium and in driven cold-atom quantum wires. Phys. Rev. Lett. 106, 220402 (2011).
  * (8) Liu, X.-J. & Hu, H. Topological superfluid in one-dimensional spin-orbit-coupled atomic Fermi gases. Phys. Rev. A 85, 033622 (2012).
  * (9) Zhang, W. & Yi, W. Topological Fulde-Ferrell-Larkin-Ovchinnikov states in spin-orbit-coupled Fermi gases. Nature Communications 4, 2711 (2013).
  * (10) Ruhman, J., Berg, E., & Altman, E. Topological states in a one-dimensional Fermi gas with attractive interaction. Phys. Rev. Lett. 114, 100401 (2015).
  * (11) Beeler, M. C. _et al._ The spin Hall effect in a quantum gas. Nature 498, 201 (2013).
  * (12) Kennedy, C. J., Siviloglou, G. A., Miyake, H., Burton, W. C. & Ketterle, W. Spin-orbit coupling and quantum spin Hall effect for neutral atoms without spin flips. Phys. Rev. Lett. 111, 225301 (2013).
  * (13) He, L. & Huang, X.-G. BCS-BEC crossover in 2D Fermi gases with Rashba spin-orbit coupling. Phys. Rev. Lett. 108, 145302 (2012).
  * (14) Xu, Y. & Zhang, C. Berezinskii-Kosterlitz-Thouless phase transition in 2D spin-orbit-coupled Fulde-Ferrell superfluids. Phys. Rev. Lett. 114, 110401 (2015).
  * (15) Galitski, V. & Spielman, I. B. Spin-orbit coupling in quantum gases. Nature 494, 49 (2013).
  * (16) Goldman, N., Juzeliūnas, G., Öhberg, P. & Spielman, I. B. Light-induced gauge fields for ultracold atoms. Rep. Prog. Phys. 77, 126401 (2014).
  * (17) Zhai, H. Degenerate quantum gases with spin-orbit coupling: a review. Rep. Prog. Phys. 78, 026001 (2015).
  * (18) Kharzeev, D. E., McLerran, L. D. & Warringa, H. J. The effects of topological charge change in heavy ion collisions: “Event by event P and CP violation”. Nucl. Phys. A 803, 227 (2008).
  * (19) Fukushima, K., Kharzeev, D. E. & Warringa, H. J. Chiral magnetic effect. Phys. Rev. D 78, 074033 (2008).
  * (20) Son, D. T. & Zhitnitsky, A. R. Quantum anomalies in dense matter. Phys. Rev. D 70, 074018 (2004).
  * (21) Metlitski, M. A. & Zhitnitsky, A. R. Anomalous axion interactions and topological currents in dense matter. Phys. Rev. D 72, 045011 (2005).
  * (22) Vilenkin, A. Equilibrium parity-violating current in a magnetic field. Phys. Rev. D 22, 3080 (1980).
  * (23) Zyuzin, A. A. & Burkov, A. A. Topological response in Weyl semimetals and the chiral anomaly. Phys. Rev. B 86, 115133 (2012).
  * (24) Vazifeh, M. M. & Franz, M. Electromagnetic response of Weyl semimetals. Phys. Rev. Lett. 111, 027201 (2013).
  * (25) Basar, G., Kharzeev, D. E. & Yee, H.-U. Triangle anomaly in Weyl semimetals. Phys. Rev. B 89, 035142 (2014).
  * (26) Landsteiner, K. Anomalous transport of Weyl fermions in Weyl semimetals. Phys. Rev. B 89, 075124 (2014).
  * (27) Li, Q. _et al._ Observation of the chiral magnetic effect in ZrTe5. arXiv: 1412.6543.
  * (28) Abelev, B. I. _et al._(STAR Collaboration) Azimuthal charged-particle correlations and possible local strong parity violation. Phys. Rev. Lett. 103, 251601 (2009).
  * (29) Abelev, B. _et al._(ALICE Collaboration) Charge separation relative to the reaction plane in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. Phys. Rev. Lett. 110, 012301 (2013).
  * (30) Skokov, V., Illarionov, A. Yu. & Toneev, V. Estimate of the magnetic field strength in heavy-ion collisions. Int. J. Mod. Phys. A 24, 5925 (2009).
  * (31) Deng, W.-T. & Huang, X.-G. Event-by-event generation of electromagnetic fields in heavy-ion collisions. Phys. Rev. C 85, 044907 (2012).
  * (32) Deng, W.-T. & Huang, X.-G. Electric fields and chiral magnetic effect in Cu + Au collisions. Phys. Lett. B 742, 296 (2015).
  * (33) Zhou, X., Li, Y., Cai, Z. & Wu, C. Unconventional states of bosons with the synthetic spin Corbit coupling. J. Phys. B: At. Mol. Opt. Phys. 46, 134001 (2013).
  * (34) van der Bijl, E. & Duine, R. A. Anomalous Hall conductivity from the dipole mode of spin-orbit-coupled cold-atom systems. Phys. Rev. Lett. 107, 195302 (2011).
  * (35) Sundaram, G. & Niu, Q. Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and Berry-phase effects. Phys. Rev. B 59, 14915 (1999).
  * (36) Xiao, D., Chang, M.-C. & Niu, Q. Berry phase effects on electronic properties. Rev. Mod. Phys. 82, 1959 (2010).
  * (37) Stephanov, M. A. & Yin, Y. Chiral kinetic theory. Phys. Rev. Lett. 109, 162001 (2012).
  * (38) Kharzeev, D. E. Topologically induced local PandCP violation in QCD$\times$QED. Ann. Phys. (N.Y.) 325, 205 (2010).
  * (39) Landsteiner, K., Megias, E. & Pena-Benitez, F. Gravitational anomaly and transport phenomena. Phys. Rev. Lett. 107, 021601 (2011).
  * (40) Vilenkin, A. Parity nonconservation and rotating black holes. Phys. Rev. Lett. 41, 1575 (1978).
  * (41) Son, D. T. & Surowka, P. Hydrodynamics with triangle anomalies. Phys. Rev. Lett. 103, 191601 (2009).
  * (42) Jiang, Y., Huang, X.-G. & Liao, J. Chiral vortical wave and induced flavor charge transport in a rotating quark-gluon plasma. Phys. Rev. D 92, 071501 (2015).
  * (43) Xiao, D., Shi, J. & Niu, Q. Berry phase correction to electron density of states in solids. Phys. Rev. Lett. 95, 137204 (2005).
  * (44) Duval, C., Horvath, Z., Horvathy, P., Martina, L. & Stichel, P. Berry phase correction to electron density in solids and “exotic” dynamics. Mod. Phys. Lett. B 20, 373 (2006).
  * (45) Son, D. T. & Yamamoto, N. Berry curvature, triangle anomalies, and the chiral magnetic effect in Fermi liquids. Phys. Rev. Lett. 109, 181602 (2012).
  * (46) Gao, J.-H., Liang, Z.-T., Pu, S., Wang, Q. & Wang, X.-N. Chiral anomaly and local polarization effect from the quantum kinetic approach. Phys. Rev. Lett. 109, 232301 (2012).
  * (47) Chen, J.-W., Pu, S., Wang, Q. & Wang, X.-N. Berry curvature and four-dimensional monopoles in the relativistic chiral kinetic equation. Phys. Rev. Lett. 110, 262301 (2013).
  * (48) Manuel, C. & Torres-Rincon, J. M. Kinetic theory of chiral relativistic plasmas and energy density of their gauge collective excitations. Phys. Rev. D 89, 096002 (2014).
  * (49) Hou, D., Liu, H. & Ren, H.-C. Some field theoretic issues regarding the chiral magnetic effect. JHEP 1105, 046 (2011).
  * (50) Banerjee, N. _et al._ Constraints on fluid dynamics from equilibrium partition functions. JHEP 1209, 046 (2012).
  * (51) Jensen,K. Triangle anomalies, thermodynamics, and hydrodynamics. Phys. Rev. D 85, 125017 (2012).
  * (52) Satow, D. & Yee, H. U. Chiral magnetic effect at weak coupling with relaxation dynamics. Phys. Rev. D 90, 014027 (2014).
  * (53) Gorbar, E. V., Miransky, V. A. & Shovkovy, I. A. Normal ground state of dense relativistic matter in a magnetic field. Phys. Rev. D 83, 085003 (2011).
  * (54) Burnier, Y., Kharzeev, D. E., Liao, J. & Yee, H.-U. Chiral magnetic wave at finite baryon density and the electric quadrupole moment of the quark-gluon plasma. Phys. Rev. Lett. 107, 052303 (2011).
  * (55) Wang, G. _et al._(STAR Collaboration) Search for Chiral magnetic effects in high-energy nuclear collisions. Nucl.Phys. A 904, 248c(2013).
  * (56) Kharzeev, D. E. & Yee, H.-U. Chiral magnetic wave. Phys. Rev. D 83, 085007 (2011).
  * (57) Huang, X.-G. & Liao, J. Axial current generation from electric field: chiral electric separation effect. Phys. Rev. Lett. 110, 232302 (2013).
  * (58) Jiang, Y., Huang, X.-G. & Liao, J. Chiral electric separation effect in the quark-gluon plasma. Phys. Rev. D 91, 045001 (2015).

Acknowledgments We acknowledge the support from Fudan University Grant No.
EZH1512519, Shanghai Natural Science Foundation No. 14ZR1403000, the Key
Laboratory of Quark and Lepton Physics (MOE) of CCNU (Grant No. QLPL20122),
the Young 1000 Talents Program of China, and Scientific Research Foundation of
State Education Ministry for Returned Scholars.

Author Contributions
X.G.H conceived and conducted the research and wrote the manuscript.
Correspondence and requests for materials should be addressed to X.G.H.
(huangxuguang@fudan.edu.cn).

Competing Interests
The author declares that he has no competing financial interests.


# Thermal Creep on Mars: Visualizing a Soil Layer under Tension

Tetyana Bila  University of Duisburg-Essen, Faculty of Physics, Lotharstr.
1-21, 47057 Duisburg, Germany  Jonathan Kollmer  University of Duisburg-Essen,
Faculty of Physics, Lotharstr. 1-21, 47057 Duisburg, Germany  Jens Teiser
University of Duisburg-Essen, Faculty of Physics, Lotharstr. 1-21, 47057
Duisburg, Germany  Gerhard Wurm  University of Duisburg-Essen, Faculty of
Physics, Lotharstr. 1-21, 47057 Duisburg, Germany

###### Abstract

At low ambient pressure, temperature gradients in porous soil lead to a gas
flow, called thermal creep. With this regard, Mars is a unique as the
conditions for thermal creep to occur in natural soil only exist on this
planet in the solar system. Known as Knudsen compressor, thermal creep induces
pressure variations. In the case of Mars, there might be a pressure maximum
below the very top dust particle layers of the soil, which would support
particle lift and might decrease threshold wind velocities necessary to
trigger saltation or reduce angles of repose on certain slopes. In laboratory
experiments, we applied diffusing wave spectroscopy (DWS) to trace minute
motions of grains on the nm-scale in an illuminated simulated soil. This way,
DWS visualizes pressure variations. We observe a minimum of motion which we
attribute to the pressure maximum $\sim 2$ mm below the surface. The motion
above but especially below that depth characteristically depends on the
ambient pressure with a peak at an ambient pressure of about 3 mbar for our
sample. This is consistent with earlier work on ejection of particle layers
and is in agreement to a thermal creep origin. It underlines the supporting
nature of thermal creep for particle lift which might be especially important
on Mars.

Mars, thermal creep, particle lift, diffusing wave spectroscopy

††journal: Planetary Science Journal

##  1 Introduction

Sand and dust are mobile on the surface of Mars, readily observable from small
scale dust devils over sand dune motion to global dust storms (Fenton, 2020;
Viúdez-Moreiras et al., 2020; Toigo et al., 2018; Lorenz et al., 2021; Heyer
et al., 2020; Chojnacki et al., 2019). This requires some lifting mechanism
and gas drag related to wind and eolian transport is undoubtedly a main driver
(Rasmussen et al., 2015). It has a long tradition though to argue that winds
on Mars are just regularly slightly too weak to account for particle lift
(Greeley et al., 1980; Kok et al., 2012).

There are ways where rather low shear stress might suffice. E.g., the
continuation threshold for saltation is way lower than the initial threshold
(Kok, 2010), aggregates are easier to be picked up (Merrison et al., 2007),
settling at low gravity might make the soil more prone to wind drag (Musiolik
et al., 2018), sporadic motion by gusts might also contribute (Swann et al.,
2020), and electrostatic forces might play a role as well (Esposito et al.,
2016; Kruss et al., 2021). So, it might be debatable if the wind on Mars has a
general problem to initiate particle motion but, in any case, if at its best
the wind speed is always only close to the entrainment threshold velocity,
then every other effect might play a very important role in providing
conditions in favor of particle mobility (Neakrase et al., 2016).

Among the supporting mechanisms for particle lift or downhill flow might also
be sub-soil processes. Resurfacing water in selected locations has been
discussed as a sediment transport mechanism (Raack et al., 2017), also
pressure excursions coming along with the vortex of a dust devil might provide
lift to surface material (Balme & Hagermann, 2006; Koester & Wurm, 2017; Bila
et al., 2020). Especially in this latter case, the pressure within the soil is
higher than above the soil.

The ambient pressure range of mbar allows one more lifting mechanism on Mars
which also leads to a subsurface overpressure, eventually. This is based on a
Knudsen compressor induced by thermal creep within the soil as the surface is
insolated. We explain this in detail in the following section. The effect can
be large but occurs on a small scale and therefore is intrinsically difficult
to be measured directly. Therefore, a visualization is the focus of our work
here.

##  2 Thermal creep induced sub soil overpressure

The underlying physical principle of thermal creep and Knudsen compressors are
not widely spread and to some degree counter intuitive, so we review the
basics here along the lines of fig. 1. Part ”a” shows the general gas flow
situation in a single open tube with a temperature gradient along its wall. In
a thin layer on the order of the mean free path of the molecules
$\lambda_{g}$, gas molecules then get a net component of motion in the
direction from cold to warm. This is owed to the non-thermal equilibrium
situation and is called thermal creep. It is just that at a given constant
pressure, on the cold side the rate of molecules traveling to the warm side is
higher than the other way round. As one usually recalls heating a gas
reservoir with expansion or pressure increase and therefore a motion from warm
to cold, this flow is not very intuitive. In large tubes, it also essentially
goes unnoticed without larger effects outside the tube as there is ample space
in the center of the tube to lead to a backflow without much pressure
difference needed.

This changes for small tubes with an overall diameter on the order of the mean
free path as pictured in fig. 1 b. Thermal creep now dominates over the total
cross section. There can certainly still be a backflow but now this would need
to be driven by a significant pressure gradient. Therefore, if, e.g., the end
of the tubes would be closed, this would lead to a pressure increase on the
warm side until thermal creep flow and pressure-driven backflow would balance.
This would then be called a Knudsen compressor (s. fig. 1 d) according to
Knudsen who studied the pressure conditions at the end of small tubes in the
early twentieth century (Knudsen, 1910). At perfect conditions, he could
realize situations where the ratio between pressures at the warm to the cold
side was on the order of factor 10.

Figure 1:  The basic principle of thermal creep and Knudsen compressor. a) Gas
flow in an open tube of the large cross section with a temperature gradient
between the tube ends. In a tube with a large cross section, there are two
opposing gas flows: a gas flow from the colder end to the warmer end in a thin
layer on the order of the mean free path of molecules $\lambda_{g}$, called
thermal creep, and a backflow from the warmer to the colder end. b) In a tube
with a cross section on the order of $\lambda_{g}$, the thermal creep gas flow
dominates. c-d) A tube with a cross section on the order of $\lambda_{g}$ with
a temperature gradient between its ends. On the warmer side the tube is
followed by a smaller tube of the same cross section having a constant
temperature. The thermal creep gas flow cannot proceed into the smaller tube
and the pressure increases at the junction. This leads to a pressure-driven
gas flow through the smaller tube. The pressure adjusts itself until the
incoming temperature gradient-driven gas flow equals the pressure-driven gas
flow. e) In-depth basic pressure distribution in a dust bed heated from above.
$y_{0}$ denotes the surface of a particle bed. In the near-surface layer from
$y_{0}$ to $y_{1}$ the temperature is constant. Below $y_{1}$ down to $y_{2}$
the temperature gradient is present.

So at its extremes, thermal creep in a tube can either lead to a gas flow of
maximum rate, if the ends of the tube are open, or thermal creep leads to a
maximum pressure increase, if the flow cannot proceed, i.e. if it is blocked
by walls or if there are adjacent volumes but if these are negligible in
sizes. The situation relevant here is somewhere in between and pictured in
fig. 1 c. If we just add a (small) tube on the warm side which is at a
constant temperature, then there is no thermal creep that can transport the
incoming flow further on, so the pressure increases and the flow speed
decreases until the pressure gradient along the small tube transports the flow
onward and the temperature-driven thermal creep flow equals the pressure-
driven gas flow.

As outlined and entering in detail, thermal creep effects depend on the ratio
$Kn$ between the mean free path of the gas molecules $\lambda_{g}$ and tube
diameter $d$. This is known as Knudsen number. If the Knudsen number is small,
tubes are relatively large and no pressure can build up. If the Knudsen number
is high, the tube is very small and the flow rate eventually gets too small to
sustain pressure differences. Therefore, thermal creep effects are most
pronounced under conditions where the Knudsen number is around $Kn\sim 1$.

To make the connection to Martian soil, the question is how the tube picture
can be applied to a porous granular bed. While the idea would be the same,
namely that thermal creep transports gas through the pores, this situation is
more complex. There are pores of varying size. So for one thing, the Knudsen
number is not well defined a priori. A simple analogy is shown in fig. 2 from
Steinpilz et al. (2017), where the granular bed is idealized by a regular
cubic structure. Koester et al. (2017) studied various real particle beds
experimentally and showed that indeed a granular medium can very well be
described by a set of tubes where the diameter of the tube essentially equals
the average particle (pore) diameter, which would be slightly larger than the
sketch suggests.

Finally, as with the two tubes in fig. 1 c, we can have different parts of a
particle bed or soil. Thermal creep can pump gas without the need for a
pressure difference in a region with a temperature gradient. If this connects
i.e. to a top layer at a constant temperature, pressure builds up until gas is
driven further in the usual way by a pressure difference which is known as
Darcy flow. More detailed equations for calculating thermal creep for granular
matter can be found in Koester et al. (2017). The pressure gradient-driven top
layer under tension has a certain thickness, which we visualize below as the
overpressure will always lead to at least minute particle motions upwards and
downwards and there will be a minimum of particle motion close to the pressure
maximum.

Figure 2: Model of a granular medium as a collection of pore-size tubes (from
Steinpilz et al. (2017)).

There have been observations of particle motions before but at such large
levels of incident light flux that the pressure gradients were large enough to
eject grains directly against gravity and cohesion. First experiments on
particle lifting at low pressure when illuminated were carried out by Wurm &
Krauss (2006) but among the possible drivers, all related to non-equilibrium
low pressure gas physics, thermal creep generated overpressure was only
recognized by de Beule et al. (2014), supported by a number of works in
between (Wurm et al., 2008; Kelling & Wurm, 2009; Kelling et al., 2011;
Kocifaj et al., 2010, 2011).

Kelling et al. (2011) find aggregates on the order of 100 $\rm\mu m$ being
ejected in experiments. Numerical simulations in Kocifaj et al. (2011) show
that the temperature changes from flat to strongly decreasing at a few mm
depth. The layer under tension was probed by laboratory measurements by de
Beule et al. (2015) based on the thickness of an ejected layer to be up to a
few 100 $\rm\mu m$.

This overpressure might not be strong enough on its own under Martian sunlight
but can reduce the threshold shear velocity for particle lift by about 20 %
(Küpper & Wurm, 2015). The effect can be amplified under certain conditions of
illumination, i.e. at a shadow boundary (Kuepper & Wurm, 2016). This might
help explaining Recurring Slope Lineae on Mars (Schmidt et al., 2017).

In any case, the top layer that is active is rather thin. While this is
consistent with temperature profiles modeled, this is only on the order of the
grain size or slightly larger. Therefore, the pressure profile cannot be
probed directly.

In the spirit of evaluating different methods for tracing the pressure or at
least the maximum pressure line, we employed a new method here, called
diffusing wave spectroscopy and we present first measurements which qualify
this technique as suitable visualization tool for thermal creep pressure
variations.

Figure 3:  Top: The experimental setup used. 1) Vacuum chamber, 2) polarizing
filter, 3) camera lens, 4) camera, 5) red filter, 6) lens with f=300 mm, 7)
lens with f=10.8 mm, 8) polarizing filter, 9) laser, 10) optical fiber, 11)
digital thermometer, 12) analog-to-digital converter, 13) pressure display
unit and 14) vacuum gauge. Middle and bottom figures show the sketch of the
setup in top view and side view respectively. A granular soil simulant sample
within a vessel is placed in a vacuum chamber. Heating from the top (IR
radiation) simulates insolation. The speckle pattern of a laser backscattered
from the sample from the side is imaged by a camera.  Figure 4:  Volume size
distribution of the basalt sample used. From this distribution the mean
particle diameter is calculated to $51\pm 28\,\mathrm{\upmu m}$.

##  3 Experiment

Figure 3 middle and bottom presents a sketch of the experimental setup. A
vessel filled with sample soil material is placed inside a vacuum chamber. The
volume of the vessel holding the sample is $(9.5\times 9\times
1)\,\mathrm{cm^{3}}$, where the sample takes a volume of $(3.8\times 9\times
1)\,\mathrm{cm^{3}}$. One side of the vessel is transparent with a 2 mm flat
glass plate.

Heating of the surface by infrared radiation is induced by a heated wire
placed about $3\,\mathrm{cm}$ above the sample surface. We did not quantify
the radiative flux here as our focus is on the DWS technique and what it might
provide as information on particle motion. As sample we use basalt with a
broad size distribution (Figure 4) and a mean diameter of $51\pm
28\,\mathrm{\upmu m}$. The bulk density of the sample is about
$1.7\,\mathrm{g/cm^{3}}$, which results in a porosity of about 0.44.

A red (633 nm), 2 mW beam shaped ($\sim$ 3 cm diameter) laser illuminates a
large part of the soil from the side, especially including the surface region.
The scattered light is recorded by a video camera. The optical pathways of
laser and camera hold crossed linear polarizers to exclude geometrically
reflected light but limit the detection to light scattered by the sample which
changes the polarization of the incoming radiation. The camera is focused on
the front glass plate plane and has a resolution of $10\,\mathrm{\upmu
m/pixel}$. Due to constructive and destructive interference it records a
speckle image that changes with slight movements of the sample particles.

The experiment was carried out at 8 different pressures $1.9\pm 0.1$ mbar,
$2.5\pm 0.1$ mbar, $3.6\pm 0.1$ mbar, $5.2\pm 0.1$ mbar, $8.2\pm 0.1$ mbar,
$148.9\pm 0.4$ mbar, $293.9\pm 0.6$ mbar and $780\pm 2$ mbar. The initial
pump-down took tens of hours (a weekend) in order to assure that most of the
water within the pore space is evaporated and not influencing the
measurements. Final adjustments were made before each measurement using a
precision valve. The temperature before the heater was switched on was
$296.1\pm 0.5$ K for all measurements. The maximum temperature at the end of
the measurement depends on the pressure. Figure 9 (bottom) shows the pressure
dependence of the absolute temperature increase ($\Delta T$) of the top soil
layers.

Figure 5:  Example of a time evolution of the correlation for the measurement
at $5.2\,\mathrm{mbar}$; times after the heater is switched on are imprinted
in each image. The dashed green line indicates the position of the sample
surface. $y$ denotes downward direction and $x$ is the horizontal. The size of
one metapixel is 20 pixel.

##  4 Diffusing wave spectroscopy

A speckle pattern is the essential part of diffusing wave spectroscopy, which
is a technique which visualizes very small displacements in granular materials
(Crassous, 2007; Erpelding et al., 2010; Amon et al., 2017; Weitz & Pine,
1993). In the used backscattering geometry, the laser beam illuminates the
granular bed from the side. The light is then entering the granular medium, is
scattered several times on various grains and part of the radiation leaves the
sample again on the same side. As mentioned, due to the coherent nature of the
light source and interference of the various scattered beams, a focused image
results in a speckle image with bright spots for constructive interference and
dark ones if pathways result in destructive intereference.

The scattering medium, in our case the soil simulant, is characterized by its
transport mean free path $\lambda$. On their way through the sample, the
photons scan a volume on average in the order of $\lambda^{3}$ (Erpelding et
al., 2008). Hence, the scattering process limits the spatial resolution of the
method to the transport mean free path. Several pixels of an image are grouped
into a _metapixel_. In our case, the speckles have a size of about $2$ pixel
$(20\rm\upmu m)$ and we take 20 by 20 pixels as metapixel. $\langle I\rangle$
is the average of the brightness $I$ over all pixels inside a metapixel.

If grains shift on the size scale of the wavelength, the speckle pattern
changes. Therefore, areas where displacements take place can be visualised in
a correlation plot comparing corresponding metapixels from images taken at
different times $t$. The intensity correlation function $G_{l}$ is calculated
for the corresponding $l$-th metapixel between the initial frame at $t=0\,$s
and one of the following frames as

$G_{l}=\frac{\langle I_{\textit{0l}}\cdot I_{\textit{tl}}\rangle-\langle
I_{\textit{0l}}\rangle\langle I_{\textit{tl}}\rangle}{\sqrt{\langle
I_{\textit{0l}}^{2}\rangle-\langle
I_{\textit{0l}}\rangle^{2}}\cdot\sqrt{\langle
I_{\textit{tl}}^{2}\rangle-\langle I_{\textit{tl}}\rangle^{2}}}.$ (1)

In eq. 1, each metapixel $l$ of a frame is represented by its brightness
$I_{\textit{tl}}$ where t denotes the frame at time $t$. $\langle
I_{\textit{0l}}\cdot I_{\textit{tl}}\rangle$ denotes the average of the
element wise multiplication of respective pixels within the two metapixels.

Figure 6:  Same as fig. 5 but without heating at $\sim 1190$ mbar ambient
pressure.

At $t=0\,$s the sample is undeformed. The deformation develops with time and
therefore can be seen as a decrease of the correlation value. A correlation
value close to 1 indicates that no or very weak displacement of scattering
particles has taken place. The correlation values can be plotted spatially
resolved, building a correlation map where each correlation value is a measure
for the deformation occurring in a cell of $\lambda^{3}$-size.

Fig. 5 shows an example of the correlation plot for the heated sample in the
illuminated part of the vessel over time. Fig. 6 shows the same as fig. 5, but
for the sample at room temperature without heating.

##  5 Results and discussion

To characterize the depth dependence of the particle motion in more detail, we
average the data horizontally. Fig. 7 then shows the time evolution of the
depth profile in absolute spatial units now. We reversed the colors here for
visibility, i.e. bright colors refer to large motions. Especially deeper
layers down to 10 mm loose correlations over time, i.e. show large particle
motion on timescales of minutes. This is consistent with heat conduction
timescales. We did not follow the evolution longer than 400 s here. While heat
is transported further downwards, the correlations do no longer change in the
top layers, eventually.

To discuss the depth profiles we therefore then choose the profile of the
latest time at 395 s for each of the 8 pressures sampled. These are shown in
fig. 8.

To show the characteristics of the depth profiles we divide the data into
different pressure ranges. Fig. 8 top shows the depth profile for the three
highest pressures from about 149 mbar to 780 mbar. No thermal creep should be
present at this high pressures and, indeed, the curves are all essentially the
same. Here we focus on the depth down to $5\,$mm. These motion depth profiles
at higher pressure can clearly be split into two parts. Down to 2 mm, the
motion decreases roughly like a power law (linear in log-log plot). Below 2 mm
the motion profile is flat. As thermal creep is not acting, the motion close
to the surface has to have a different origin. We cannot pin down the reason
but in general, the topmost grains are bound the least, are not compressed by
the grains close by and therefore are the most mobile. So they are moved the
easiest by any disturbance. The further down, the more confined are the grains
until disturbances are too faint to result in any motion. Disturbances might
be simple gas flow or expansion as the sample heats up but, again, we do not
know the exact reason. In any case, these profiles provide a firm base for
comparison to the lower pressure range.

Figure 7:  Strength of the motion of grains with depth y over time at 5.2 mbar
ambient pressure. Dark color marks little motion (reversed color from
correlation plot fig. 5).

To make the difference clear, we now overplot one example of a depth profile
for the lower pressure range (5.2 mbar) in fig. 8 middle where we expect
thermal creep to act. For completeness, we plot all low pressure data in fig.
8 bottom. In comparison to the higher pressure range, the low pressure data
have a simple but very distinguished shape with 3 points to note.

  * •

The topmost decrease of motion is still a power law and essentially one with
the same power as for higher pressures but the absolute value is shifted
upwards. This implies more motion in the top layer. This is consistent with a
subsurface overpressure with a pressure gradient above that moves all grains
up to the very top somewhat more.

  * •

There is a clear minimum for all low pressure data at about 2 mm. The position
of this minimum does not change significantly.

  * •

Finally, the motion now increases again deeper within the particle bed. This
is also consistent with a subsurface overpressure that works in both
directions upwards but also downwards.

The motion of grains requires forces acting on them. These can be provided by
the sub-soil pressure gradients induced by thermal creep. As indicated in fig.
1 e, there is a pressure maximum close to the surface, as the upward thermal
creep gas flow can only be maintained by a pressure-driven gas flow in the
upper layers with little temperature gradient. Further down below the surface,
along the flat temperature distribution in the near-surface layer, there will
be another transition to a flat temperature distribution, which might result
in pressure variations, though this is not shown in fig. 1 d which was focused
on the top layer.

Figure 8:  Motion depth profiles at the latest time (395 s); top: high
pressure data with clear split in two parts, a downward decreasing and then
nearly constant part; middle: overplotted example of low pressure data with an
additional increase at larger depths. The orange dashed line is a fit with
$r_{high}=const.$ on the flat parts of all three of the higher pressure depth
profiles. $r_{low}$ is the value of the depth profile of lowest pressure at
$y=5\,$mm. $\Delta r$ is the difference between $r_{low}$ and $r_{high}$;
bottom: all low pressure data.

It seems reasonable to assume that these pressure variations drive the motion
of grains seen in fig. 7. It has to be noted though that the location of
particle motion and the location of gas pressure extrema might not necessarily
coincide 1 to 1 a priori. The gas pressure gradient implies certain local
forces on the grains but forces on grains will in addition be commuted to the
other grains by force chains. In other words, the minimum of particle motion
traces an equilibrium spot where upward and downward directed forces would
balance but even a bulk motion, i.e. large scale lift or compression would be
seen as motion in DWS as it always comes with small variations. Large scale
motions are not present in the data though. The sample further below the
surface only moves significantly once the temperature rises in that region as
the heat is conducted downwards. Any motion further down in the particle bed
occurring then is not significantly influencing the layers above, so pressure
gradients and particle motion can be considered to act only locally.

Figure 9:  Top: pressure dependence of $\Delta r$; Center: dependence of
$\Delta r$ on Knudsen number; Bottom: Temperature difference at the end of the
measurement over the ambient pressure at the beginning of the measurement.

The most pronounced differences with ambient pressure are the strengths of the
motion increase below the surface. We therefore studied the differences
$\Delta r=r_{low}-r_{high}$ over pressure and, indeed, we observe a systematic
variation as shown in fig. 9 (top). The effect is strongest at an ambient
pressure of about 3 mbar rapidly decreasing to the offset value as the
pressure is increased by a few mbar. In terms of the Knudsen number
$\rm{Kn}=\lambda_{g}/d$ with the mean free path of the gas molecules
$\lambda_{g}$ and an efficient pore diameter of $d=20\upmu\mathrm{m}$ the
maximum is located at about Kn = 1 (fig. 9 (center)).

This trend is consistent with earlier work by de Beule et al. (2015). They
find that the layer removed in their experiments at three pressures of 0.1, 1,
and 10 mbar is largest at 1 mbar. So we clearly trace the same effect. More
specific, we find that the maximum is situated at a Knudsen number of 1. As
thermal creep is most efficient at Kn=1 (Chambers, 2004) and does not work at
large ambient pressure and less at lower ambient pressure (Koester et al.,
2017), this strongly supports the idea that thermal creep is the underlying
mechanism for the particle motion.

##  6 Caveats

###  6.1 Boundary effects

The vessel holding the sample has a relatively small thickness. This is
intended as it e.g. allows a faster cooling after one measurement before the
next is started. However, it also means that the boundaries, i.e. the glass
windows on the side have a significant influence on the temperature
distribution. This cannot be avoided as we need a transparent wall that allows
DWS at all and the particle columns need a stabilizing wall. With an order of
magnitude of 1 W/(m K), the thermal conductivity of the walls is orders of
magnitudes larger than the thermal conductivity of the sample. Therefore, the
temperature of the window next to the relevant top few mm can be considered as
being constant. So this is a heat source for the sample but overlaying as
constant temperature it will not change the principle vertical stratification
of the temperature differences and therefore, the position of the pressure
maximum will only change slightly.

###  6.2 Other motion

The maximum in motion at around 10 mm (s. fig. 7) is also interesting but is
not the focus of this work. We currently have no sophisticated model to
explain these data.

###  6.3 Temperature dependence on pressure

Fig. 9 (bottom) shows that the overall temperature increases at low pressure
but since it reaches a rather constant plateau, we do not consider this to
influence our thermal creep measurements.

##  7 Conclusion

Thermal creep gas flow can significantly support the lifting of grains on
Mars. However, the pressure variations are small and change on sub-mm spatial
scales. Implementing arrays of commercial pressure sensors on sub-mm scale
within a dust bed and resolutions on the Pa level without influencing the
temperature profile if the dust bed is illuminated is currently virtually
impossible. Therefore, the pressure gradient itself is hard to be measured
directly. With this lack of direct confirmation, even if the effect of thermal
creep on larger scales is well known, every indirect method can be an
important verification. After all, thermal creep is not something encountered
on Earth in natural settings and it does not come intuitively (we think).

Tracing the gas motion in and out of a dust bed was a milestone to attribute
particle lifting to thermal creep gas flow (de Beule et al., 2014). Observing
the direct lifting of whole layers of grains in a work by de Beule et al.
(2015) could pin down an active layer of $\sim 0.2$ mm within a dust bed,
where we refer to activity as being under tension.

Here, this pressure maximum is visualized by subtle particle motions made
visible by diffusing wave spectroscopy. The depth of about 2 mm is slightly
larger than these earlier results but noting the different methodology and
observing the same typical dependence of particle motion on ambient pressure,
these findings are in agreement. Our results therefore provide one more
confirmation of the effect of thermal creep under Martian conditions making
thermal creep an ever more likely effect to influence the motion of grains on
the Martian surface quite generally.

## acknowledgements

This project is supported by DLR Space Administration with funds provided by
the Federal Ministry for Economic Affairs and Climate Action (BMWK) under
grant numbers 50WM1943 and 50WM2049. This project also has received funding
from the European Union’s Horizon 2020 research and innovation program under
grant agreement No 101004052. The manuscript gained significantly by the
reviews of Frédéric Schmidt and an anonymous reviewer.

## References

  * Amon et al. (2017) Amon, A., Mikhailovskaya, A., & Crassous, J. 2017, Review of Scientific Instruments, 88, 051804, doi: 10.1063/1.4983048
  * Balme & Hagermann (2006) Balme, M., & Hagermann, A. 2006, Geophys. Res. Lett., 33, L19S01, doi: 10.1029/2006GL026819
  * Bila et al. (2020) Bila, T., Wurm, G., Onyeagusi, F. C., & Teiser, J. 2020, Icarus, 339, 113569, doi: 10.1016/j.icarus.2019.113569
  * Chambers (2004) Chambers, A. 2004, Modern Vacuum Physics, Masters Series in Physics and Astronomy (CRC Press), 25–48
  * Chojnacki et al. (2019) Chojnacki, M., Banks, M. E., Fenton, L. K., & Urso, A. C. 2019, Geology, 47, 427, doi: 10.1130/G45793.1
  * Crassous (2007) Crassous, J. 2007, The European Physical Journal E, 23, 145, doi: 10.1140/epje/i2006-10079-y
  * de Beule et al. (2015) de Beule, C., Wurm, G., Kelling, T., Koester, M., & Kocifaj, M. 2015, Icarus, 260, 23, doi: 10.1016/j.icarus.2015.06.002
  * de Beule et al. (2014) de Beule, C., Wurm, G., Kelling, T., et al. 2014, Nature Physics, 10, 17, doi: 10.1038/nphys2821
  * Erpelding et al. (2008) Erpelding, M., Amon, A., & Crassous, J. 2008, Physical Review E, 78, doi: 10.1103/physreve.78.046104
  * Erpelding et al. (2010) —. 2010, EPL (Europhysics Letters), 91, 18002, doi: 10.1209/0295-5075/91/18002
  * Esposito et al. (2016) Esposito, F., Molinaro, R., Popa, C. I., et al. 2016, Geophysical Research Letters, 43, 5501, doi: 10.1002/2016GL068463
  * Fenton (2020) Fenton, L. K. 2020, Icarus, 352, 114018, doi: 10.1016/j.icarus.2020.114018
  * Greeley et al. (1980) Greeley, R., Leach, R., White, B., Iversen, J., & Pollack, J. B. 1980, Geophys. Res. Lett., 7, 121, doi: 10.1029/GL007i002p00121
  * Heyer et al. (2020) Heyer, T., Raack, J., Hiesinger, H., & Jaumann, R. 2020, Icarus, 351, 113951, doi: 10.1016/j.icarus.2020.113951
  * Kelling & Wurm (2009) Kelling, T., & Wurm, G. 2009, Phys. Rev. Lett., 103, 215502, doi: 10.1103/PhysRevLett.103.215502
  * Kelling et al. (2011) Kelling, T., Wurm, G., Kocifaj, M., Klačka, J., & Reiss, D. 2011, Icarus, 212, 935, doi: 10.1016/j.icarus.2011.01.010
  * Knudsen (1910) Knudsen, M. 1910, Annalen der Physik, 336, 633, doi: 10.1002/andp.19103360310
  * Kocifaj et al. (2011) Kocifaj, M., Klačka, J., Kelling, T., & Wurm, G. 2011, Icarus, 211, 832, doi: 10.1016/j.icarus.2010.10.006
  * Kocifaj et al. (2010) Kocifaj, M., Klačka, J., Wurm, G., Kelling, T., & Kohút, I. 2010, MNRAS, 404, 1512, doi: 10.1111/j.1365-2966.2010.16370.x
  * Koester et al. (2017) Koester, M., Kelling, T., Teiser, J., & Wurm, G. 2017, Ap&SS, 362, 171, doi: 10.1007/s10509-017-3154-4
  * Koester & Wurm (2017) Koester, M., & Wurm, G. 2017, Planet. Space Sci., 145, 9, doi: 10.1016/j.pss.2017.07.005
  * Kok (2010) Kok, J. F. 2010, Phys. Rev. Lett., 104, 074502, doi: 10.1103/PhysRevLett.104.074502
  * Kok et al. (2012) Kok, J. F., Parteli, E. J. R., Michaels, T. I., & Karam, D. B. 2012, Reports on Progress in Physics, 75, 106901, doi: 10.1088/0034-4885/75/10/106901
  * Kruss et al. (2021) Kruss, M., Salzmann, T., Parteli, E., et al. 2021, Planetary Science Journal, 2, 238, doi: 10.3847/PSJ/ac38a4
  * Kuepper & Wurm (2016) Kuepper, M., & Wurm, G. 2016, Icarus, 274, 249, doi: 10.1016/j.icarus.2016.02.049
  * Küpper & Wurm (2015) Küpper, M., & Wurm, G. 2015, Journal of Geophysical Research (Planets), 120, 1346, doi: 10.1002/2015JE004848
  * Lorenz et al. (2021) Lorenz, R. D., Lemmon, M. T., & Maki, J. 2021, Icarus, 364, 114468, doi: 10.1016/j.icarus.2021.114468
  * Merrison et al. (2007) Merrison, J. P., Gunnlaugsson, H. P., Nørnberg, P., Jensen, A. E., & Rasmussen, K. R. 2007, Icarus, 191, 568, doi: 10.1016/j.icarus.2007.04.035
  * Musiolik et al. (2018) Musiolik, G., Kruss, M., Demirci, T., et al. 2018, Icarus, 306, 25, doi: 10.1016/j.icarus.2018.01.007
  * Neakrase et al. (2016) Neakrase, L. D. V., Balme, M. R., Esposito, F., et al. 2016, Space Sci. Rev., 203, 347, doi: 10.1007/s11214-016-0296-6
  * Raack et al. (2017) Raack, J., Conway, S. J., Herny, C., et al. 2017, Nature Communications, 8, 1151, doi: 10.1038/s41467-017-01213-z
  * Rasmussen et al. (2015) Rasmussen, K. R., Valance, A., & Merrison, J. 2015, Geomorphology, 244, 74, doi: 10.1016/j.geomorph.2015.03.041
  * Schmidt et al. (2017) Schmidt, F., Andrieu, F., Costard, F., Kocifaj, M., & Meresescu, A. G. 2017, Nature Geoscience, 10, 270, doi: 10.1038/ngeo2917
  * Steinpilz et al. (2017) Steinpilz, T., Teiser, J., Koester, M., Schywek, M., & Wurm, G. 2017, MicST, 29, 325, doi: 10.1007/s12217-017-9550-0
  * Swann et al. (2020) Swann, C., Sherman, D. J., & Ewing, R. C. 2020, Geophys. Res. Lett., 47, e84484, doi: 10.1029/2019GL084484
  * Toigo et al. (2018) Toigo, A. D., Richardson, M. I., Wang, H., Guzewich, S. D., & Newman, C. E. 2018, Icarus, 302, 514, doi: 10.1016/j.icarus.2017.11.032
  * Viúdez-Moreiras et al. (2020) Viúdez-Moreiras, D., Newman, C. E., Forget, F., et al. 2020, Journal of Geophysical Research (Planets), 125, e06493, doi: 10.1029/2020JE006493
  * Weitz & Pine (1993) Weitz, D. A., & Pine, D. J. 1993, in Monographs on the physics and chemistry of materials, Vol. 49, Dynamic Light Scattering: The Method and Some Applications, ed. W. Brown (Clarendon Press), 652 – 720
  * Wurm & Krauss (2006) Wurm, G., & Krauss, O. 2006, Phys. Rev. Lett., 96, 134301, doi: 10.1103/PhysRevLett.96.134301
  * Wurm et al. (2008) Wurm, G., Teiser, J., & Reiss, D. 2008, Geophys. Res. Lett., 35, L10201, doi: 10.1029/2008GL033799