# Coffee-stain growth dynamics on dry and wet surfaces

François Boulogne  Department of Mechanical and Aerospace Engineering,
Princeton University, Princeton, NJ 08544, USA  Laboratoire Matière et
Systèmes Complexes (MSC), UMR 7057 CNRS, Université Paris Diderot, Bâtiment
Condorcet, 10 rue Alice Domon et Leonie Duquet, Paris, France  Laboratoire de
Physique des Solides, CNRS, Univ. Paris-Sud, Université Paris-Saclay, Orsay
91405, France  François Ingremeau  Department of Mechanical and Aerospace
Engineering, Princeton University, Princeton, NJ 08544, USA  LIPhy, CNRS, and
Université Grenoble Alpes, 140 Rue de la Physique, 38402 Saint-Martin-d’Hères,
France  Howard A. Stone  Department of Mechanical and Aerospace Engineering,
Princeton University, Princeton, NJ 08544, USA

###### Abstract

The drying of a drop containing particles often results in the accumulation of
the particles at the contact line. In this work, we investigate the drying of
an aqueous colloidal drop surrounded by a hydrogel that is also evaporating.
We combine theoretical and experimental studies to understand how the
surrounding vapor concentration affects the particle deposit during the
constant radius evaporation mode. In addition to the common case of
evaporation on an otherwise dry surface, we show that in a configuration where
liquid is evaporating from a flat surface around the drop, the singularity of
the evaporative flux at the contact line is suppressed and the drop
evaporation is homogeneous. For both conditions, we derive the velocity field
and we establish the temporal evolution of the number of particles accumulated
at the contact line. We predict the growth dynamics of the stain and the
drying timescales. Thus, dry and wet conditions are compared with experimental
results and we highlight that only the dynamics is modified by the evaporation
conditions, not the final accumulation at the contact line.

##  1 Introduction

The investigation of the evaporation of a pure liquid drop started in the late
1800s with papers by Maxwell [1] and Stefan [2], and continued in 1918 with a
further contribution by Langmuir [3], for the specific case of a spherical
drop suspended in air. Also, the evaporation of disk-shaped liquid drops has
been motivated by the understanding of plant transpiration [4, 5, 6] and to
predict the lifetime of droplets [7, 8, 9]. The addition of small particles to
the liquid drop leaves, upon evaporation, a ring stain along the contact line,
as commonly observed for drops of tea or coffee. This so-called coffee-stain
effect was first investigated about 20 years ago by Deegan et al. [10, 11]. It
is now well-established that the accumulation of the particles at the contact
line is the result of an outward radial flow owing to evaporation [10, 12, 13,
14, 15]. This accumulation pins the contact line and the loss of the volatile
solvent decreases the drop height [16, 17]. Consequently, to replenish the
loss of solvent at the pinned contact line, the liquid flows radially. The
characteristics of this flow transporting the particles in the drop depends on
the flow field in the drop and on the details of the evaporation profile.

The laminar radial flow field commonly observed in simple volatile liquids can
be modified by changing the composition of the liquid. Thermal or solutal
Marangoni effects have been extensively studied in the literature, as they can
contribute to a chaotic flow in a binary solvent [18] or a reverse flow from
surfactant concentration gradients [18]. Also, instead of evaporating the
liquid from the liquid-vapor interface, the absorption by the substrate is
known to modify the particle deposition to achieve a spatially more uniform
coating [19].

For simple liquids characterized by a laminar radial flow, the transport of
particles suspended in a evaporating drop depends solely on the evaporative
flux. To describe the particle accumulation at the contact line, it is
necessary to establish first the evaporative flux, which is function of the
solvent diffusion in the gas phase, and the shape of the liquid-vapor
interface [10, 12, 20, 21, 22]. Mathematically, it is established that for a
sessile drop on a non-evaporating surface, the evaporative flux at the liquid-
vapor interface diverges at the contact line [10, 12, 20]. By using the
analogy between diffusive concentration fields and electrostatic potential
fields, the divergence of the evaporative flux can be understood as a tip-
shape effect [23, 11].

Therefore, the evaporation profile is well-established for a single drop on a
dry surface. However, in many common situations and industrial settings, drops
in the presence of other surrounding droplets, a liquid film or a porous
substrate imbibed with volatile liquid are also subject to evaporation. For
instance, the mutual influence of drops [24, 25, 26] can be crucial for
sprayed surface coatings where deposited droplets are closely spaced. As a
result, the vapor concentration field around a drop can vary substantially in
the presence of a drying environment. In particular, we can expect the
evaporative flux at the surface of such drops to be more homogeneous. Thus, as
the liquid flow is a function of the evaporation and as it transports
suspended particles [27, 28, 29, 30], the resulting number of particles
accumulated at the contact line may be expected to depend on the environmental
conditions. To the best of our knowledge, a detailed comparison between a
spatially homogeneous flux of solvent and a diverging flux has not been
reported to date.

In previous studies [19, 31], we focused our attention on the absorption of a
drop containing micrometer size particles by a swelling hydrogel. Instead of
absorption, our present study is motivated by the drying of a drop containing
microparticles on a wet hydrogel. In such a configuration, the drop and the
gel evaporate. Therefore, we expect that the evaporation of the drop is
modified by these new environmental conditions and that the transport of
colloidal particles is affected.

In this paper, we investigate theoretically the coffee stain effect of a drop
placed in the center of an evaporating surface that is much larger than the
drop diameter and we illustrate the concepts with experimental observations.
First, we argue theoretically that the divergence of the flux at the contact
line is strongly reduced in this geometry compared to the classic sessile drop
on a dry surface. Our study focuses on dilute suspensions, such that particles
are expected to act as tracers and do not disturb the velocity field. However,
their presence at the contact line helps to maintain a constant radius of the
drop and we ignore the second regime where the contact line eventually
recedes. The main result of this paper concerns the time evolution of the
particle density close to the contact line, which is different in the two
cases of a single drop on a non-evaporating surface and a drop on a drying
surface. This observation supports the suppression of the diverging flux at
the contact line. Nevertheless, we show that for both conditions, the
particles are transported toward the contact line such that the final patterns
are similar. Therefore, we conclude that the divergence of the evaporating
flux at the contact line for a drop on a dry surface is not crucial to obtain
coffee stain effect. Based on the model and the experimental results, the
coffee stain mechanism is discussed in section 3.2.

The outline of the paper is organized as follow. In Section 2, we first recall
the evaporation dynamics of a single drop on a non-evaporating surface. Then,
we derive the evaporation dynamics of a large circular hydrogel to show that a
small drop placed in the center of this gel evaporates with a nearly uniform
flux. Furthermore, we calculate the velocity fields in both cases as well as
the temporal evolution of the number of particles accumulated at the contact
line. In Section 3, we present our experimental setup for the preparation of
the dry and wet configurations. Then, we compare our experimental results to
the theoretical predictions regarding the drop lifetime and the time evolution
of the number of particles accumulated at the contact line.

##  2 Model

Figure 1:  Sketches of the (a) dry and (b) wet configurations that also
illustrate the materials used experimentally to achieve both cases. Each drop
is sitting on a thin glass disk of radius $R$ and the profile of the drop
interface is denoted $h(r,t)$. The thin glass disk prevents a solvent flux
between the drop and the substrate and ensures that the surface conditions are
the same in both cases for the particles in the drop. The vapor concentrations
at saturation and at infinity are denoted $c_{s}$ and $c_{\infty}$,
respectively. The thickness of arrows indicates qualitatively the strength of
the evaporative flux.

###  2.1 Drop shape

We consider two distinct situations depicted in Figure 1 where a liquid drop
evaporates on a substrate. We assume that the geometry is axisymmetric and we
use cylindrical $(r,z)$ coordinates. The capillary length is
$\ell_{c}=\sqrt{\gamma/(\rho g)}$, where $\gamma$ is the surface tension, $g$
is the gravitational acceleration and $\rho$ the liquid density. We assume
that the contact line at $r=R$ is pinned during the evaporation, which is
consistent with the experiments where we observe that the drop has a constant
radius. This defines the constant radius evaporation mode and the drop height
decreases in time. Nevertheless, the contact line eventually recedes at the
very end of the drying process, when the liquid film height is typically of
the order of a few times the particle size.

For a drop radius $R<\ell_{c}$, to a good approximation the drop shape is a
spherical cap. Thus, at each time, the drop profile $h(r,t)$ (Figure 1) is
described by a portion of a sphere,

$h(r,t)=\sqrt{\left(\frac{R^{2}+h_{0}(t)^{2}}{2h_{0}(t)}\right)^{2}-r^{2}}-\frac{(R^{2}-h_{0}(t)^{2})}{2h_{0}(t)},$
(1)

with $h_{0}(t)=h(0,t)$. This solution ensures that the curvature is only a
function of time, such that the capillary pressure is uniform in the drop and
capillary effects alone do not generate any flow.

We assume that the drop height is much smaller than the drop radius, $h_{0}\ll
R$, which corresponds to small contact angles. This assumption is necessary in
the later sections to apply the lubrication approximation to develop
analytical results. Therefore, the drop shape is approximated by

$h(r,t)\simeq h_{0}(t)\,\left(1-\frac{r^{2}}{R^{2}}\right).$ (2)

With the assumption of a pinned contact line ($\textrm{d}R/\textrm{d}t=0$),
the time derivative of the drop height $h(r,t)$ defined by equation (2) is
given by

$\frac{\partial h}{\partial t}=\left(1-\frac{r^{2}}{R^{2}}\right)\,\frac{{\rm
d}h_{0}}{{\rm d}t}.$ (3)

The drop has the volume of a spherical cap
$\Omega(t)=\frac{\pi}{2}h_{0}(t)\left(\frac{h_{0}(t)^{2}}{3}+R^{2}\right)\simeq\pi
R^{2}h_{0}(t)/2$. We define $Q_{e}(t)$ as the total evaporative flux given by
$Q_{e}(t)=\int_{S}v_{e}\mbox{d}S$, where $S$ designates the liquid-vapor
interface of the drop and $v_{e}$ denotes the local evaporation speed, i.e.
the volume of liquid that evaporates per unit area per unit time. The time
variation of this volume corresponds to the loss of liquid by evaporation,
i.e. ${\rm d}\Omega/{\rm d}t=-Q_{e}(t)$. Therefore, equation (3) simplifies to

$\frac{\partial h}{\partial
t}=-\left(1-\frac{r^{2}}{R^{2}}\right)\frac{2Q_{e}(t)}{\pi R^{2}}.$ (4)

In the next section, we present the evaporation dynamics in two different
environmental situations.

###  2.2 Diffusion-limited evaporation

We aim to analyze the difference between drops evaporating on a dry surface
and on a surface that presents an evaporative flux of the same liquid as the
drop. To do so, we need to calculate the evaporation speed in these distinct
cases. First, we present the well-known equations for diffusion-limited
evaporation. Then, we establish the vapor concentration surrounding a disk of
gel to define the conditions for which the evaporative flux can be considered
to a good approximation as homogeneous close to its center. Finally, we give
the expressions of the evaporative flux for a single drop on a dry surface and
for a drop surrounded by an evaporating surface.

####  2.2.1 Theoretical background

We denote by $c_{v}$ the vapor mass concentration in the gas surrounding the
evaporating liquid. For diffusion-limited evaporation [23, 32], the temporal
and spatial evolution of the vapor concentration follows

$\frac{\partial c_{v}}{\partial t}={\cal D}\nabla^{2}c_{v},$ (5)

where ${\cal D}$ is the diffusion coefficient of the vapor in the gas phase.

We consider the characteristic length scale ${\cal L}$ of the vapor
concentration gradient, which is related to the geometry of the system at long
timescales [33, 34]. We also denote ${\cal V}_{e}$ as the typical value of the
evaporation speed $v_{e}$. In particle-tracking experiments reported below,
the relative effect of the diffusive time of the vapor in the gas phase
$\tau_{D}={\cal L}^{2}/{\cal D}$ and the evaporation time $\tau_{e}={\cal
L}/{\cal V}_{e}$ defines the Péclet number:

$\textrm{Pe}=\frac{{\cal V}_{e}\,{\cal L}}{{\cal D}}.$ (6)

The appropriate choice of the length scale ${\cal L}$ is made in the next
sections for each geometry. The diffusion coefficient ${\cal D}$ for water
vapor in air at room temperature is ${\cal D}=2\times 10^{-5}$ m2/s. In our
experimental conditions, the typical evaporation speed for water is ${\cal
V}_{e}\approx 10^{-8}$ m/s (see Appendix). If the Péclet number is small, then
equation (5) can be simplified to the Laplace equation

$\nabla^{2}c_{v}=0.$ (7)

The evaporation velocity at the liquid-vapor interface is [35]

$v_{e}(r,t)=-\frac{{\cal D}}{\rho}{\bf
n}\cdot{\bm{\nabla}}c_{v}\qquad\mbox{at}\quad z=h(r,t)$ (8)

where $\rho$ is the liquid density, $h(r,t)$ is the vertical position of the
liquid-vapor interface and ${\bf n}$ the unit normal vector directed into the
gas phase.

####  2.2.2 Dry configuration

Under the assumption of diffusion-limited evaporation, we recall the
evaporation dynamics of a sessile drop of radius $R_{d}$ (Figure 1a). The
characteristic length scale ${\cal L}$ is the glass disk radius, thus
$\textrm{Pe}\ll 1$. For a contact angle $\theta_{c}$, the evaporation speed is
given by [10, 11]

$v_{e}^{d}(r)=\frac{2{\cal
D}(c_{s}-c_{\infty})}{\pi\rho}\left\\{\begin{array}[]{ll}\frac{f(\lambda)}{(R^{2}-r^{2})^{\lambda}}&\mbox{for
}r<R,\\\ 0&\mbox{for }r>R,\end{array}\right.$ (9)

where $\lambda=\frac{\pi-2\theta_{c}}{2\pi-2\theta_{c}}$ and $f(\lambda)$ is a
known function. For $\theta_{c}\ll 1$, we have $\lambda=1/2$ and
$f(\lambda)=1$ [11, 20], where we recover the particular solution of a flat
disk of radius $R$ given by equation (13). Consequently, the total evaporative
flux of the drop in this dry configuration is

$Q_{e}^{d}=\frac{4{\cal D}(c_{s}-c_{\infty})R}{\rho}.$ (10)

From equation (4) and $Q_{e}^{w}$ from equation (10), which is time
independent, the height at the center of the drop varies according to

$h_{0}(t)=h_{i}\left(1-\frac{t}{\tau_{e}^{d}}\right),$ (11)

where the initial drop height is $h_{i}=h(r=0,t=0)$ and
$\tau_{e}^{d}=\frac{\rho h_{i}R\pi}{8{\cal D}(c_{s}-c_{\infty})}$ is the
evaporation time.

####  2.2.3 Wet configuration

Next, the evaporation dynamics of the wet configuration represented in Figure
1(b) is investigated. We consider a circular slab of hydrogel of radius
$R_{g}$, without the glass disk and the drop. The typical length scale for the
vapor concentration gradient is the gel radius $R_{g}$ and, for example, for
${\cal L}=R_{g}=1.75$ cm, we have $\textrm{Pe}\ll 1$. This geometry
corresponds to the particular case $\theta_{c}=0$ in equation (9). The vapor
concentration field can be simply written as [36, 37]

${\cal
C}_{v}^{w}(\kappa,\sigma)=\frac{c_{v}^{w}(\kappa,\sigma)-c_{\infty}}{(c_{s}-c_{\infty})}=1-\frac{2}{\pi}\arctan(\sigma),$
(12)

where we use the oblate spheroidal coordinates ($\kappa,\sigma$) defined as
$r^{2}=R_{g}^{2}(1-\kappa^{2})(1+\sigma^{2})$ and $z=R_{g}\kappa\sigma$. The
evaporation velocity at $z=0$ is given by

$v_{e}^{w}(r)=\frac{2{\cal
D}(c_{s}-c_{\infty})}{\pi\rho}\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{R_{g}^{2}-r^{2}}}&\mbox{for
}r<R_{g},\\\ 0&\mbox{for }r>R_{g},\end{array}\right.$ (13)

which is equivalent to equation (9) with $R$ replaced by $R_{g}$ and
$\theta_{c}=0$.

The representation of the vapor concentration surrounding the hydrogel disk is
shown in Figure 2(a) and we focus our interest in the center of the circular
disk of radius $R_{g}$ where the drop will be placed, i.e. $r/R_{g}<0.1$. Near
the edge of the slab, we expect a divergence of the evaporative flux as
observed for a sessile drop. As is clear from Figures 2(b-c), the vapor
concentration and its vertical gradient above the center of a flat circular
slab only weakly depends on the radial and vertical coordinates. This makes
the hydrogel surface a good candidate to study the evaporation of a small
droplet exposed to a homogeneous evaporative flux near the center of the gel
disk.

Figure 2: Vapor concentration obtained from equation (12). (a) Vapor
concentration map surrounding a drying disk of radius $R_{g}$ obtained from
equation (12). (b) Radial vapor concentration profile at $z/R_{g}=0.1$. (c)
Vertical vapor concentration profiles at different radial positions
$\tilde{r}=r/R_{g}$.

As a consequence of the results shown in Figure 2, the drying velocity
$v_{e}^{w}$ has a weak variation for $r\in[0,R]$ if $R\ll R_{g}$. Therefore,
we will consider that in the configuration similar to the one represented in
Figure 1(b), but with $R<0.1R_{g}$, the drying velocity of the drop is uniform
and is equal to

$v_{e}^{w}=\frac{2{\cal D}(c_{s}-c_{\infty})}{\pi\rho R_{g}},$ (14)

where $\frac{{\cal D}(c_{s}-c_{\infty})}{\rho}=2\times 10^{-10}$ m2/s in our
experimental conditions. Thus, for an air-vapor interfacial area approximated
as $\pi R^{2}$, the evaporative flux from the drop is

$Q_{e}^{w}=v_{e}^{w}\,\pi R^{2}=\frac{2{\cal D}(c_{s}-c_{\infty})R^{2}}{\rho
R_{g}},$ (15)

considering that in the wet configuration, the local evaporative flux is
uniform. As expected on physical grounds, in the present experimental
conditions, the average evaporation velocity for the dry configuration
$\bar{v}_{e}^{d}=Q_{e}^{d}/(\pi R^{2})>v_{e}^{w}$ since $R<R_{g}$.

Similarly to equation (11), the height at the center of the drop follows

$h_{0}(t)=h_{i}\left(1-\frac{t}{\tau_{e}^{w}}\right),$ (16)

where $\tau_{e}^{w}=\frac{h_{i}}{2v_{e}^{w}}$ is the evaporation time.

###  2.3 Convective evaporation

Our previous derivations assume that the evaporation is limited by diffusion.
In this section, we briefly analyze the effect of convection on evaporation.
Recent studies have addressed the importance of convection on the evaporation
of drops [38, 39, 40, 41, 42, 43]. In the case of water, humid air is less
dense than dry air, which leads to natural convection in the surrounding air.

In order to close the problem, the density of air as to be expressed as a
function of the water concentration, which is related to the partial pressure
of vapor. We denote $p_{d}$ and $p_{v}$ the partial pressures of dry and humid
air, respectively. The density of air is then

$\rho_{m}=\frac{p_{d}M_{d}+p_{v}M_{v}}{{\cal R}T},$ (17)

where $M_{d}$ and $M_{v}$ are the molar masses of dry and wet air,
respectively, ${\cal R}$ is the ideal gas constant and $T$ is the absolute
temperature.

Therefore, we now consider the steady state convectively driven evaporation
above a circular disk of radius $R_{g}$ with a boundary-layer flow and use the
Boussinesq approximation (Figure 3) [44, 45]. The convective flow is driven by
the density difference between humid air located near the liquid-vapor
interface and the ambient air. We denotes $(u_{r},u_{z})$ as the velocity
field in the vapor phase in cylindrical coordinates. The continuity equation
is

$\frac{1}{r}\frac{\partial(ru_{r})}{\partial r}+\frac{\partial u_{z}}{\partial
z}=0$ (18)

and in the steady state, the radial momentum equation is

$u_{r}\frac{\partial u_{r}}{\partial r}+u_{z}\frac{\partial u_{r}}{\partial
z}=-\frac{1}{\rho_{m}}\frac{\partial p}{\partial r}++\nu\nabla^{2}u_{r},$ (19)

where $\nu$ is the kinematic viscosity of vapor and $\rho_{m}$, the vapor
phase density, considered as a constant in the framework of the Boussinesq
approximation. The vertical momentum equation is written,

$u_{r}\frac{\partial u_{z}}{\partial r}+u_{z}\frac{\partial u_{z}}{\partial
z}=-\frac{1}{\rho_{m}}\frac{\partial p}{\partial z}+\nu\nabla^{2}u_{z}.$ (20)

The vertical pressure gradient is given by

$\frac{\partial p}{\partial z}=(\rho_{m}-\rho_{\infty})g$ (21)

where $\rho_{\infty}$ is the air density far from the disk. This gradient is
the driving force for the convection, which induces a vertical motion
described by equation (20). The steady state vapor mass conservation satisfies

$u_{r}\frac{\partial c_{v}}{\partial r}+u_{z}\frac{\partial c_{v}}{\partial
z}={\cal D}\frac{\partial^{2}c_{v}}{\partial z^{2}}.$ (22)

On the surface of the disk, the boundary conditions are
$u_{r}(r,z=0)=u_{z}(r,z=0)=0$ and $c_{v}(r,z=0)=c_{s}$. Far from the disk,
$u_{r}(r,z\rightarrow\infty)=0$ and $c_{v}(r,z\rightarrow\infty)=c_{\infty}$.

The mass Grashof number compares the buoyancy forces to the viscous forces
represented, respectively, by the first and the second terms on the right-hand
side of both equations (19) and (20),

${\rm
Gc}=\left|\frac{\rho_{s}-\rho_{\infty}}{\rho_{\infty}}\right|\left(\frac{gR_{g}^{3}}{\nu^{2}}\right),$
(23)

where $\rho_{s}$ and $\rho_{\infty}$ are, respectively, the vapor density at
saturation and at infinity. For ${\rm Gc}\gg 1$, we expect that convection
takes place above the evaporating surface and that the diffusion-limited
solution given by equation (12) is not valid. As the mass Grashof number is
proportional to $R_{g}^{3}$, convection must be important for large
evaporating surfaces or very volatile liquids [46]. For instance, for
conditions typical of our experiments with a circular gel of radius
$R_{g}=1.75$ cm, setting the humidity of the glove box $c_{\infty}=0.5c_{s}$
leads to ${\rm Gc}=10^{3}$. The Schmidt number $Sc=\nu/{\cal D}$ describes the
relative thickness of the velocity and the vapor concentration boundary layers
and appears in the set of equations (18-22) to be solved. For air at room
temperature, the kinematic viscosity is $\nu=1.6\times 10^{-5}$ m2/s and ${\rm
Sc}=0.8$.

Figure 3:  Free-convection driven evaporative flux $\tilde{v}_{e}^{c}$ near
the center of the gel disk along the radius for ${\rm Sc}=0.8$ and ${\rm
Gc}=6\times 10^{4}$. The radial range corresponds to the droplet position
($R/R_{g}\approx 0.1$). The sketch introduces the notations with the
approximate flow field indicated.

Since we are interested in the evaporative flux of a drop placed at the center
of the gel, we focus on the evaporation rate of the gel at this location. To
solve equations (18-22), we follow the numerical scheme presented by [47] for
the free convection above the center of a heated horizontal circular disk. In
this numerical scheme, the absolute value of $\tilde{v}_{e}^{c}(\tilde{r})$
cannot be predicted by this numerical scheme. Nevertheless, we obtain the
typical variation of the evaporative flux for small radii where we place a
drop in our experiments. In Figure 3, we show that the evaporative flux
increases along the radius more significantly than for the diffusion-limited
situation.

From this analysis, we expect that the evaporative flux of a drop placed in
the center of the gel does not diverge at the contact line as a drop in a dry
configuration. In the following, as a first approximation, we will assume that
the evaporative flux is uniform in spite of the presence of convection above
the hydrogel disk. Thus, in this wet configuration, we will consider that
convection only increases the evaporation rate compared to the diffusion-
limited situation. So, the characteristics of the flow in the drop resulting
from the evaporation is independent of the Grashof number, which affects only
the evaporation rate and the flow speed.

###  2.4  Evaporation of dilute colloidal drops

In the previous section we characterized the velocity of the vapor at the
liquid-gas interface for both dry and wet configurations. Since we are
interested in the motion of the colloidal particles in the drop and more
specifically in the formation of particle deposits, we need to describe the
flow field in the drop during the constant radius evaporation mode.

In the following, we assume that the colloidal suspension is dilute so that
any accumulation does not affect the drop shape and the flow field; also,
hydrodynamic interactions are neglected and particles are assumed to follow
the streamlines.

####  2.4.1 Hydrodynamics

Let $(v_{r},v_{\theta},v_{z})$ denote the velocity field in the liquid in
cylindrical coordinates $(r,\theta,z)$. The problem is invariant by rotation
around the $z$ axis. Thus, $v_{\theta}=0$, $\partial/\partial\theta=0$, and we
only consider $v_{r}(r,z,t)$ and $v_{z}(r,z,t)$. Under this assumption, the
continuity equation is

$\frac{1}{r}\frac{\partial(rv_{r})}{\partial r}+\frac{\partial v_{z}}{\partial
z}=0.$ (24)

Also, the Navier-Stokes equation, analyzed within the lubrication
approximation, $h_{0}\ll R$ [48], provides a linear stress balance along the
$r$-direction

$\eta\frac{\partial^{2}v_{r}}{\partial z^{2}}=\frac{\partial p}{\partial r},$
(25)

where $\eta$ is the viscosity and $p$ is the pressure. A first integration of
this equation with respect to the variable $z$ gives

$\frac{\partial v_{r}}{\partial z}=\frac{1}{\eta}\frac{\partial p}{\partial
r}\left(z-h(r,t)\right)$ (26)

where we used the boundary condition of a zero stress at the liquid-vapor
interface $\frac{\partial v_{r}}{\partial z}(z=h,t)=0$. With the second
boundary condition of a zero tangential velocity on the substrate
$v_{r}(z=0,t)=0$, a second integration leads to a Poiseuille flow

$v_{r}(r,z,t)=\frac{1}{\eta}\frac{\partial p}{\partial
r}\left(\frac{z^{2}}{2}-h(r,t)z\right).$ (27)

We define the average velocity $\bar{v}_{r}(r,t)$ over the drop height
$h(r,t)$

$\displaystyle\bar{v}_{r}(r,t)=\frac{1}{h(r,t)}\int_{0}^{h(r,t)}v_{r}(r,z,t)\,{\rm
d}z=-\frac{h(r,t)^{2}}{3\eta}\frac{\partial p}{\partial r}.$ (28)

Substituting equation (28) in (27), we have

$v_{r}(r,z,t)=\frac{3\bar{v}_{r}(r,t)}{2h(r,t)^{2}}(2h(r,t)z-z^{2}).$ (29)

To obtain the average velocity $\bar{v}_{r}(r,t)$ as a function of the
evaporative flux of solvent, we use the local continuity equation obtained by
a mass conservation in a slice between $r$ and $r+{\rm d}r$, which yields

$\frac{\partial h}{\partial
t}+\frac{1}{r}\frac{\partial(hr\bar{v}_{r})}{\partial r}+v_{e}(r)=0,$ (30)

where the time derivative of the liquid height is given by equation (4). In
the next sections, we solve this equation for two distinct evaporating
conditions.

####  2.4.2 Uniform evaporation for a wet surface

Let us consider first the wet configuration. We assume that the drop dries
with a velocity $v_{e}^{w}$ independent of the radial coordinate $r$ and
constant in time. This assumption is supported by the analysis made in section
2.2.3 and 2.3 for the diffusive and the convective regimes respectively. By
integrating equation (30) with respect to $r$, we have

$\bar{v}_{r}^{w}(r,t)=\frac{r\,v_{e}^{w}}{2h(r,t)}\left(1-\frac{r^{2}}{R^{2}}\right)=\frac{rv_{e}^{w}}{2h_{0}(t)},$
(31)

where we used the condition $\bar{v}_{r}^{w}(r=0,t)=0$. Hence,
$\bar{v}_{r}^{w}$ is a linear function of $r$. Consequently, the radial
velocity $v_{r}^{w}(r,z,t)$ given by equation (29) is

$v_{r}^{w}(r,z,t)=\frac{3rv_{e}^{w}}{4h_{0}(t)h(r,t)^{2}}\left(2h(r,t)z-z^{2}\right).$
(32)

Now, we calculate the vertical component of the velocity field
$v_{z}^{w}(r,z,t)$. From equations (24) and (32), we have

$v_{z}^{w}(r,z,t)=\frac{v_{e}^{w}}{2}\left(\frac{z^{3}}{h(r,t)^{3}}\left(1+\frac{r^{2}}{R^{2}}\right)-\frac{3z^{2}}{h(r,t)^{2}}\right)$
(33)

where we used the boundary condition $v_{z}^{w}(r,z=0,t)=0$.

####  2.4.3 Evaporation with flux singularity at the contact line for a dry
surface

To compare this result with the configuration of a drop on a dry surface, we
now consider the flux singularity at the drop contact line. From equations (4)
and (10), we have

$\frac{\partial h}{\partial t}=-\frac{8{\cal D}(c_{s}-c_{\infty})}{\pi\rho
R}\left(1-\frac{r^{2}}{R^{2}}\right).$ (34)

and the mean radial velocity can be derived from (30),

$\bar{v}_{r}^{d}(r,t)=\frac{2{\cal
D}(c_{s}-c_{\infty})}{\pi\rho}\frac{R}{rh(r,t)}\\\
\left(\left(1-\frac{r^{2}}{R^{2}}\right)^{1/2}-\left(1-\frac{r^{2}}{R^{2}}\right)^{2}\right)$
(35)

with the boundary condition $\bar{v}_{r}^{d}(r,z=0,t)=0$. Therefore, we obtain
for the horizontal velocity the equation (36) below and, from the continuity
equation (24), the vertical velocity given by equation (37):

$v_{r}^{d}(r,z,t)=-\frac{3{\cal D}(c_{s}-c_{\infty})}{\pi\rho
h(r,t)^{3}}\frac{R}{r}\left(\left(1-\frac{r^{2}}{R^{2}}\right)^{1/2}-\left(1-\frac{r^{2}}{R^{2}}\right)^{2}\right)(z^{2}-2h(r,t)z).$
(36)

$v_{z}^{d}(r,z,t)=\frac{3{\cal D}(c_{s}-c_{\infty})}{\pi\rho
Rh(r,t)^{3}}\left[\frac{2h_{0}}{h}(z^{3}-2hz^{2})\left(\left(1-\frac{r^{2}}{R^{2}}\right)^{1/2}-\left(1-\frac{r^{2}}{R^{2}}\right)^{2}\right)\right.\\\
+\left.\left(\frac{z^{3}-3hz^{2}}{3}\right)\left(4\left(1-\frac{r^{2}}{R^{2}}\right)-\left(\frac{1}{1-\frac{r^{2}}{R^{2}}}\right)^{1/2}\right)\right].$
(37)

####  2.4.4 Typical flow fields

Figure 4:  Streamlines corresponding to the velocity fields (a-b) for the wet
configuration, given by equations (32) and (33), and (c-d) for the dry
configuration, given by equations (36) and (37). The drop aspect ratio is
$h_{0}(t)/R=0.1$. The black solid line represents the liquid-vapor interface.
(e) Mean radial dimensionless velocity $\bar{v}_{r}^{\star}$ defined as
$\bar{v}_{r}^{w}/(v_{e}^{w}R/(2h_{0}(t)))$ for the wet configuration and
$\bar{v}_{r}^{d}/(2{\cal D}(c_{s}-c_{\infty})/(\pi\rho h_{0}(t)))$ for the dry
configuration with $h_{0}(t)/R=0.1$. (f) Schematic that illustrates the
principle of the calculation for the number of particles accumulated at the
contact line. The particles contained in the dark blue area at $t=0$ are
transported to the contact line at time $t$.

To illustrate the results from the previous two sections, we present in Figure
4(a-d) the streamlines for both configurations. In Figure 4(e), we plot the
average radial velocity as a function of the radial position. The average
radial velocity diverges near the contact line in the dry configuration
because of the divergence of the evaporative flux. In the wet configuration,
this divergence does not exist and the radial velocity increases linearly
along the radial position.

###  2.5 Particle transport

In this section, we establish the time evolution of the number of particles
accumulated at the contact line [12, 27, 29, 15]. These particles are advected
by the radial flow, which has a velocity $v_{r}(r,z,t)$, equations (32 and
(36). In the following, we assume that the particle concentration is
independent of the $z$-coordinate and so we consider the average radial
velocity $\bar{v}_{r}(r,t)$. Thus, the radial position $r(t)$ of a narrow
slice of liquid is given by

$\bar{v}_{r}(r(t),t)=\frac{\mbox{d}r}{\mbox{d}t}.$ (38)

We denote $r_{0}(t)$ the radial position for which particles initially present
at this position at $t=0$ reach the contact line at time $t$. This position
can be calculated by integrating equation (38).

The number of particles $N_{cum}(t)$ accumulating at the contact line is the
sum of the particles contained in the volume between $r_{0}(t)$ and $R$ (Fig.
4(f)).

$N_{cum}(t)=2\pi\int_{r_{0}(t)}^{R}c_{0}h(r^{\prime},t=0)\,r^{\prime}\mbox{d}r^{\prime},$
(39)

where $c_{0}$ is the initial particle concentration. Therefore, equation (39)
with the drop shape $h(r,t)$ given by equation (2) becomes

$N_{cum}(t)=2\pi
c_{0}h_{i}R^{2}\left(\frac{1}{4}-\frac{1}{2}\frac{r_{0}(t)^{2}}{R^{2}}+\frac{1}{4}\frac{r_{0}(t)^{4}}{R^{4}}\right),$
(40)

which depends on the radial velocity through $r_{0}(t)$, which has to be
determined, and thus on the evaporative conditions. The initial drop height is
$h_{i}=h(r=0,t=0)$.

We consider first the wet configuration. Using equations (31) and (16) in
equation (38), we have

$\frac{r_{0}^{w}(t)}{R}=\left(1-\frac{t}{\tau_{e}^{w}}\right)^{1/4}.$ (41)

Thus, the number of particles accumulated at the contact line can be simply
calculated by using equations (40) and (41).

Now, we consider the dry configuration. The time evolution of the drop height
follows equation (16) with an evaporation time $\tau_{e}^{d}=\frac{\rho
h_{i}R\pi}{8{\cal D}(c_{s}-c_{\infty})}$. With equations (35) and (38), we
have

$\frac{X}{\sqrt{X}-X^{2}}\,{\rm
d}X=\frac{1}{2}\frac{1}{(1-t/\tau_{e}^{d})}\,{\rm d}(t/\tau_{e}^{d}),$ (42)

where $X=1-(r_{0}^{d}/R)^{2}$. Thus, after integration, we obtain

$\frac{r_{0}^{d}(t)}{R}=\sqrt{1-\left(1-\left(1-\frac{t}{\tau_{e}^{d}}\right)^{3/4}\right)^{2/3}}.$
(43)

A more general solution taking into account the width of the ring has been
found previously by Popov [27]. We recover this solution in the limit of a
thin ring observed for dilute suspensions. From equation (40) with equations
(41) or (43), we have a prediction for the time evolution of the number of
particles $N_{cum}$ accumulated at the contact line for both evaporating
conditions. In the next section, we compare these predictions to experimental
observations.

##  3 Experiments

###  3.1 Materials and methods

Our experiments of controlled drying of drops are performed in a home-made
glove box ($30\times 30\times 30$ cm3), which has humidity regulation based on
an Arduino and a humidity sensor placed far from the drying drop (Figure 5).
Dry air is produced by circulating ambiant air in a container filled with
dessicant made of anhydrous calcium sulfate (Drierite) and moist air is
obtained by bubbling air in water. The relative humidity is set to $50$% in
all of our experiments. The glove box is placed on the top of an inverted
fluorescence microscope (Leica DMI4000 B) equipped with a 4$\times$ objective
(Plan achromatic, Olympus) and a Hamamatsu camera (digital camera ORCA-
Flash4.0 C11440). The acquisition is automated with the software Micromanager
[49].

Colloidal suspensions consist of fluorescent carboxylate-modified particles of
$1$ $\mu$m diameter (Lifetechnologies), which are suspended in deionized water
at a concentration of 2000 particles/mm3, which corresponds to a volume
fraction of $\phi=8\times 10^{-6}$. The suspension is surfactant free and no
Marangoni effect is observed during the experiment. The surface tension
measured with a pendant drop technique is $\gamma=70\pm 1$ mN/m. For all drop
evaporation experiments, the substrate consists of a circular glass disk of
diameter $2R=3$ mm and of thickness 0.15 mm (Thomas Scientific). Prior to each
experiment, a new glass disk is placed in a plasma chamber for 5 min to
enhance its wettability. A drop of $1.5$ $\mu\ell$ is dispensed with a
micropipette (Research Plus, Eppendorf) on the glass disk. This volume is
chosen to fully cover the glass disk with the colloidal suspension and it
ensures that all the drops have the same initial size and shape.

Two configurations are investigated, dry and wet surfaces, for the drop
evaporation experiments; the configurations effectively change the local
humidity of the drop’s environment between the two extremes. For experiments
with dry surfaces, to conveniently deposit a drop on the glass disk, the disk
is placed on a 3 mm thick PDMS slab bonded to a glass slide, which holds the
disk in place by adhesion. In contrast, a humid surface is made of a
polydimethyl acrylamide (PDMA) hydrogel [50]. The composition of the gel is
set by $m_{mono}/(m_{mono}+m_{w})=0.1$ and $[\rm{MBA}]/[\rm{mono}]=2\times
10^{-2}$, where $m_{mono}$ and $m_{w}$ are, respectively, monomer
(N,N-dimethylacrylamide) and water masses and $[\rm{MBA}]$, $[\rm{mono}]$ are,
respectively, molar concentrations of the crosslinker (N,N’-methylene-bis-
acrylamide) and the monomer. The quantity of initiator (potassium persulfate
and N,N,N’,N’-tetramethylethylenediamine) is set to a molar ratio of 1% of the
monomer quantity. All chemicals are purchased from Sigma-Aldrich, USA. The
solution is poured in a mold made of two glass plates ($75\times 50$ cm2, Dow
Corning) separated with a rubber spacer of $3$ mm thickness (McMaster-Carr).
After three days, the gel is removed from the mold, swollen in deionized water
overnight and its surface is gently wiped to remove any excess of water. A
hydrogel disk of $3.5$ cm diameter is prepared and a glass disk is placed at
its center.

Figure 5:  Schematic of the experimental setup. A PID controller injects dry
or moist air into the chamber based on the relative humidity measured by the
sensor.

Fluorescence images are taken in the plane of the glass disk surface every
$1.5$ s for the dry surface and every $3$ s for the wet surface. Particle
trajectories detected on fluorescence images are calculated by using the
library Trackpy [51]. These trajectories are calculated on sequences of $6$
and $20$ images for dry and wet surfaces, respectively.

The data analysis is performed with the programming language Python, the
scientific stack numpy-scipy [52, 53, 54] and Pandas [55]. The detection of
the position of the glass disk on a brightfield image is performed with
Scikit-image by using a Canny filter and the circular Hough transform [56].

Experimentally, we observe that the particles have an outward radial motion.
Thus, to evaluate the time evolution of the number of particles accumulated at
the contact line, we use the particle trajectories. For each time sequence, we
consider the trajectories that cross a circle centered on the glass disk with
a radius $0.92\,R$. We count as positive the trajectories that have an outward
direction and negative the inward trajectories. In both experimental
conditions studied in this work, the flow is outward and inward trajectories
can be observed on short timescales because of Brownian motion.

###  3.2 Experimental observations and discussion

We designed experiments to track the motion of particles, and their
corresponding accumulation at the contact line as described in the previous
section. From the particle trajectories calculated by particle tracking, we
are able to count the particles that are approaching the contact line, and to
measure their velocities.

Figure 6:  (a) Time evolution of the number of particles accumulated at the
contact line for dry and wet configurations. The time is non-dimensionalized
by the time $\tau_{e}$ at which the contact line recedes. The solid lines
represent the equation (40) for $N_{cum}(t)$ for $r_{0}^{d}(t)$ and
$r_{0}^{w}(t)$, respectively, given by equations (41) and (43). (b) Final
deposit for dry (left) and wet (right) configurations. The scale bar
represents 500 $\mu$m.

In Figure 6(a), we show for both dry and wet configurations, the time
evolution of the particles accumulated at the contact line. These values are
non-dimensionalized by the total number of particles known from the particle
concentration and the experimental evaporation time is normalized by the time
at which the contact line starts to recede.

In the theoretical model, all the particles are assumed to reach the contact
line since there is a radial flow that replenishes the contact line during the
constant radius evaporation mode (Figure 4). Indeed, the vertical velocity
$v_{z}$ is always zero at $z=0$, which means that the particles can reach the
solid surface only through Brownian motion, except in the vicinity of the
contact line where the liquid height is comparable to the particle diameter.
Experimentally, the particles can be trapped before the end of the evaporation
when the liquid film has a thickness of the order of magnitude of the particle
size. In fact, some particles do not reach the contact line and are deposited
in the inner part of the wetted area as illustrated in Figure 6(b).
Nevertheless, we observe that changing the evaporation conditions does not
affect significantly the quantity of particles forming the ring.

The percentage of particles that do not reach the drop edge is about $12$% and
$14$% for the dry and wet configurations, respectively. These particles can
originate either because they did not reach the drop edge before the contact
line recedes or because of sedimentation. However, since the number of counted
particles is not significantly different for the two drying times
$\tau_{e}^{d}$ and $\tau_{e}^{w}$, sedimentation can be neglected.

As displayed in Figure 6(a), the theoretical predictions based on equation
(40) are in good agreement with our experimental observations. The agreement
for the dry case confirms previous results obtained by Deegan [12], Popov [27]
and later by Berteloot et al. [15], among others, whereas the agreement of
theory and measurements for the wet configuration validates our choice of a
homogeneous evaporative flux at the liquid-vapor interface. A more accurate
prediction would probably require a more detailed description of the
evaporative flux at the surface of the drop. We conclude that the dynamics of
the particle accumulation at the contact line strongly depends on the spatial
variation of the evaporative flux but does not affect significantly the final
ring pattern.

In the dry configuration, equation (35) is combined with equation (11) and we
use the definition of the drying time $\tau_{e}^{d}$ to obtain

$\bar{v}_{r}^{d}(r,t)=\frac{r}{4\tau_{e}^{d}}f\left(\frac{r^{2}}{R^{2}}\right)\left(\frac{1}{1-t/\tau_{e}^{d}}\right),$
(44)

where $f$ is a function of $r^{2}/R^{2}$. This equation shows that the average
radial velocity $\bar{v}_{r}^{d}$ diverges in time for $t=\tau_{e}^{d}$. This
feature has been indicated by Gelderblom et al. [22] for drops in a dry
configuration and we show that the same mechanism is obtained in the wet
configuration. The corresponding equation to equation (44) for the wet
configuration is

$\bar{v}_{r}^{w}(r,t)=\frac{r}{4\tau_{e}^{w}}\left(\frac{1}{1-r^{2}/R^{2}}\right)\left(\frac{1}{1-t/\tau_{e}^{w}}\right).$
(45)

In figure 7, the comparison of the mean radial velocity deduced from the
particle tracking with the theoretical prediction given by equation (44) shows
good agreement.

Figure 7:  Mean radial velocity $\bar{v}_{r}^{w}(r,t)$ as a function of the
dimensionless time $t/\tau_{e}^{w}$. The radial position corresponds to
$r\approx 0.9R$.

We must also compare the predictions for the drying times with the
experimental observations. These times are $\tau_{e}^{d}\approx 1000$ s and
$\tau_{e}^{w}\approx 8000$ s for the dry and wet configurations, respectively.
In the dry configuration, the mass Grashof number (equation (23)) for the drop
is ${\rm Gc}=0.7$. From the diffusion-limited model, the prediction
$\tau_{e}^{d}=\frac{\rho h_{i}R\pi}{8{\cal D}(c_{s}-c_{\infty})}\approx 1500$
s is close to the observed value. The small discrepancy can be mainly
explained by the assumption of a small contact angle $\theta_{c}$ for the
evaporative flux in the model, which overestimates the evaporative flux at the
beginning of the experiment. Stauber et al. [8] provide an expression of the
constant radius duration including the effect of the time evolution of the
contact angle. We can numerically estimate this time to be $\approx 1150$ s.
For the wet configuration, the Grashof number for the hydrogel disk is ${\rm
Gc}=10^{3}$, which highlights the importance of buoyant forces and convection
in the vapor phase. Thus, as expected, the prediction of the drying time from
diffusion-limited model $\tau_{e}^{w}=\frac{h_{i}}{2v_{e}^{w}}\approx 34000$
s, is much larger than the experimental value.

Recent studies [40, 42] investigated the effect of convection on the
evaporation of highly volatile droplets and they rationalized their findings
with the mass Grashof number. For mass Grashof numbers between $5$ to
$10^{5}$, they observed experimentally that the total evaporation velocity
${\cal V}^{w}$ can be written in the form

${\cal V}^{w}=v_{e}^{w}(1+v_{c}^{\star}({\rm Gc})),$ (46)

where $v_{e}^{w}$ is the diffusion-limited evaporation velocity and
$v_{c}^{\star}({\rm Gc})$ is a dimensionless velocity, which is a power law of
the Grashof number, $v_{c}^{\star}\propto{\rm Gc}^{\beta}$ with $\beta\approx
0.2$ [40, 42]. Therefore, from expression (46) and their phenomenological
parameters, we can estimate ${\cal V}^{w}\in[$3.5, 4.5$]v_{e}^{w}$. Thus, the
corrected theoretical prediction for the drop lifetime in wet conditions
$\tau_{e}^{w}$ is between $7500$ and $9000$ s, which is in reasonable
agreement with our observations. Consequently, the drying of a hydrogel disk
cannot be considered to be limited by the diffusion of vapor in air. The
convection accelerates the drying velocity of the drop.

Finally, we would like to stress that the diffusion-limited evaporative
conditions can be relevant for some applications. Indeed, different solvent
properties or smaller systems such as a surface covered with water droplets of
few tens of micrometer of diameter can evaporate with negligible convection
effects.

##  4 Conclusion

In this paper, we have studied the drying dynamics of drops containing
colloidal microspheres in a dilute regime in two drying configurations. The
first configuration corresponds to a sessile drop dispensed on a dry flat
surface and is well-studied in the literature. As the evaporation proceeds,
particles migrate toward the contact line and form a solid ring of packed
particles that pins the contact line. The radial flux is the consequence of
the contact line pinning to compensate the loss of solvent. However, the
details of the radial velocity depends on the evaporation profile at the drop
interface. For the configuration of a drop on a dry surface, the evaporative
flux diverges at the contact line.

We explored the “wet” configuration of a drop placed at the center of a drying
surface. Theoretically, depending on the mass Grashof number, the evaporation
can be limited by the diffusion of vapor in the air or dominated by
convection. For diffusion-limited evaporation, we established that the
evaporative flux at the surface of the drop is homogeneous. When convection is
dominant, a simple analysis of the air flow above the center of the disk
indicates that the evaporative flux has a variation along the radius, which is
much smaller than the divergence observed for a droplet on a dry surface. We
concluded that, as a first approximation, a homogeneous evaporative flux can
be also applied as a boundary condition for convection-dominated evaporation.
Thus, we derived the velocity profiles in the drop for a homogeneous
evaporation speed and we calculated the time dependence of the number of
particles accumulated at the contact line.

We compared experimentally the effect of the dry and wet conditions on the
dynamics of the particle accumulation under these two evaporating conditions.
Our observations are in good agreements with the two respective models, which
validates our approach regarding the homogeneous evaporative flux for the wet
configuration. Also, as predicted by our model, we did not observe a
significant variation of the final particle number density at the ring
depending on the evaporating conditions. Finally, the evaporation timescales
confirm that a millimeter size sessile droplet evaporates according to a
diffusion-limited model. However, for a droplet in the center of a centimeter
size evaporating gel, convection above the disk dominates as suggested by the
large mass Grashof number. This significantly reduces the drop lifetime
compared to the prediction of pure-diffusive evaporation. As we indicated that
convection can be important for the drop evaporation surrounded by an
evaporative surface, future works should treat more precisely the effect of
convection to describe more accurately the drying conditions.

## Acknowledgments

F.B. acknowledges that the research leading to these results received funding
from the People Programme (Marie Curie Actions) of the European Union’s
Seventh Framework Programme (FP7/2007-2013) under REA grant agreement 623541.
We thank J. Dervaux and L. Limat for helpful discussions.

##  5 Appendix

In the main text, we assumed that a hydrogel dries at the same speed as pure
water. In order to validate this assumption, we measured the drying velocity
of water and gel in the same environmental conditions. The drying kinetics of
a water reservoir and a gel slab of circular shape (radius $R_{g}=1.75$ cm)
are recorded by placing them on a scale in a glove box at $50$% relative
humidity. The evolution of the mass versus time is plotted in Figure 8. The
weights decrease linearly with time at the same rate for the water and the
gel. These measurements validate the assumption that the drop deposited in the
center of the piece of gel is drying with the velocity given by equation (14).

Figure 8:  Time evolution of the weight of hydrogel disk of radius
$R_{g}=1.75$ cm and water in a short cylindrical container of the same radius.
The initial weight is arbitrarily shifted to distinguish the data points.
Black solid lines are linear fits with a slope $-0.032\pm 0.001$ mg/s.
Measurements are done in a controlled atmosphere at a relative humidity
$R_{H}=50$ %.

The total evaporative flux of water from the gel $\mathcal{Q}_{e}$ is given by
the integration of equation (13) over the total surface of the gel,

$\mathcal{Q}_{e}=\int_{0}^{R_{g}}v_{e}^{w}\pi r{\rm d}r=\frac{4{\cal
D}(c_{s}-c_{\infty})R_{g}}{\rho}.$ (47)

From the linear fit of data presented in Figure 8, we can deduce the total
evaporative flux

${\cal Q}_{e}=\frac{1}{\rho}\frac{\mbox{d}m}{\mbox{d}t}=3.2\times
10^{-11}~{}\textrm{m}^{3}/\textrm{s}.$ (48)

The theoretical value defined by equation (47) gives ${\cal Q}_{e}=1.4\times
10^{-11}~{}\textrm{m}^{3}/\textrm{s}$. Thus, the experimental average flux is
about a factor of two larger than the theoretical value, which supports a
convection enhanced evaporation.

## References

  * [1] J. C. Maxwell.  Diffusion, in Collected Scientific Papers.  Encyclopedia Britannica, Cambridge, 1877.
  * [2] J. Stefan.  Ueber die verdampfung aus einem kreisförmig oder elliptisch begrenzten becken.  Annalen der Physik, 253(11):550–560, 1882.
  * [3] I. Langmuir.  The evaporation of small spheres.  Phys. Rev., 12:368–370, 1918.
  * [4] N. Thomas and A. Ferguson.  On evaporation from a circular water surface.  Philosophical Magazine Series 6, 34:308–321, 1917.
  * [5] N. Thomas and A. Ferguson.  On the reduction of transpiration observations.  Annals of Botany, 31:241–255, 1917.
  * [6] H. Jeffreys.  Some problems of evaporation.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 35:270–280, 1918.
  * [7] H.G. Houghton.  A study of the evaporation of small water drops.  Journal of Applied Physics, 4:419–424, 1933.
  * [8] J.M. Stauber, S.K. Wilson, B.R. Duffy, and K. Sefiane.  On the lifetimes of evaporating droplets.  Journal of Fluid Mechanics, 744, 2014.
  * [9] J. M. Stauber, S. K. Wilson, B. R. Duffy, and K. Sefiane.  On the lifetimes of evaporating droplets with related initial and receding contact angles.  Physics of Fluids, 27(12), 2015.
  * [10] R.D. Deegan, O. Bakajin, T.F. Dupont, G. Huber, S.R. Nagel, and T.A. Witten.  Capillary flow as the cause of ring stains from dried liquid drops.  Nature, 389:827–829, 1997.
  * [11] R. Deegan, O. Bakajin, T. Dupont, G. Huber, S. Nagel, and T. Witten.  Contact line deposits in an evaporating drop.  Phys. Rev. E, 62:756–765, 2000.
  * [12] R. Deegan.  Pattern formation in drying drops.  Phys. Rev. E, 61:475–485, 2000.
  * [13] G. Berteloot, C.-T. Pham, A. Daerr, F. Lequeux, and L. Limat.  Evaporation-induced flow near a contact line: Consequences on coating and contact angle.  EPL, 83:14003, 2008.
  * [14] C.-T. Pham, G. Berteloot, F. Lequeux, and L. Limat.  Dynamics of complete wetting liquid under evaporation.  EPL, 92:54005, 2010.
  * [15] G. Berteloot, A. Hoang, A. Daerr, H.P Kavehpour, F. Lequeux, and L. Limat.  Evaporation of a sessile droplet: Inside the coffee stain.  Journal of Colloid and Interface Science, 370:155–161, 2012.
  * [16] H. Hu and R. Larson.  Evaporation of a sessile droplet on a substrate.  The Journal of Physical Chemistry B, 106:1334–1344, 2002.
  * [17] H. Hu and R. Larson.  Analysis of the microfluid flow in an evaporating sessile droplet.  Langmuir, 21:3963–3971, 2005.
  * [18] H. Kim, F. Boulogne, E. Um, I. Jacobi, E. Button, and H. A. Stone.  Controlled uniform coating from the interplay of Marangoni flows and surface-adsorbed macromolecules.  Phys. Rev. Lett., 116:124501, 2016.
  * [19] F. Boulogne, F. Ingremeau, J. Dervaux, L. Limat, and H. A. Stone.  Homogeneous deposition of particles by absorption on hydrogels.  EPL, 112(4):48004, 2015.
  * [20] J. Eggers and L.M. Pismen.  Nonlocal description of evaporating drops.  Physics of Fluids, 22, 2010.
  * [21] A. G. Marin, H. Gelderblom, D. Lohse, and J. H. Snoeijer.  Order-to-disorder transition in ring-shaped colloidal stains.  Phys. Rev. Lett., 107:085502, Aug 2011.
  * [22] H. Gelderblom, O. Bloemen, and J. Snoeijer.  Stokes flow near the contact line of an evaporating drop.  Journal of Fluid Mechanics, 709:69–84, 2012.
  * [23] R.G. Picknett and R. Bexon.  The evaporation of sessile or pendant drops in still air.  Journal of Colloid and Interface Science, 61:336–350, 1977.
  * [24] M. Sokuler, G. K. Auernhammer, C. J. Liu, E. Bonaccurso, and H.-J. Butt.  Dynamics of condensation and evaporation: Effect of inter-drop spacing.  EPL, 89(3):36004, 2010.
  * [25] G. Laghezza, E. Dietrich, J. M. Yeomans, R. A. Ledesma-Aguilar, S. Kooij, H. Zandvliet, and D. Lohse.  Collective and convective effects compete in patterns of dissolving surface droplets.  Soft Matter, 12:5787–5796, 2016.
  * [26] O. Carrier, N. Shahidzadeh-Bonn, R. Zargar, M. Aytouna, M. Habibi, J. Eggers, and D. Bonn.  Evaporation of water: evaporation rate and collective effects.  Journal of Fluid Mechanics, 798:774–786, 7 2016.
  * [27] Y. Popov.  Evaporative deposition patterns: Spatial dimensions of the deposit.  Phys. Rev. E, 71:036313, 2005.
  * [28] R. Zheng.  A study of the evaporative deposition process: Pipes and truncated transport dynamics.  The European Physical Journal E, 29(2):205–218, 2009.
  * [29] C. Monteux and F. Lequeux.  Packing and sorting colloids at the contact line of a drying drop.  Langmuir, 27:2917–2922, 2011.
  * [30] C. N. Kaplan and L. Mahadevan.  Evaporation-driven ring and film deposition from colloidal droplets.  Journal of Fluid Mechanics, 781, 10 2015.
  * [31] F. Boulogne, F. Ingremeau, L. Limat, and H. A. Stone.  Tuning the receding contact angle on hydrogels by addition of particles.  Langmuir, 2016.
  * [32] C. Poulard, G. Guéna, and A. M. Cazabat.  Diffusion-driven evaporation of sessile drops.  Journal of Physics: Condensed Matter, 17:S4213, 2005.
  * [33] E. Sultan, A. Boudaoud, and M. Ben Amar.  Evaporation of a thin film: diffusion of the vapour and Marangoni instabilities.  Journal of Fluid Mechanics, 543:183–202, 2005.
  * [34] F. Boulogne, A. Sauret, B. Soh, E. Dressaire, and H.A. Stone.  Mechanical tuning of the evaporation rate of liquid on crossed fibers.  Langmuir, 31:3094–3100, 2015.
  * [35] A.-M. Cazabat and G. Guena.  Evaporation of macroscopic sessile droplets.  Soft Matter, 6:2591–2612, 2010.
  * [36] N. Lebedev.  Special Functions and their Applications.  Courier Corporation, 1965.
  * [37] B.M. Marder and N.R. Keltner.  Heat flow from a disk by separation of variables.  Numerical Heat Transfer, 4:485–497, 1981.
  * [38] N. Shadizadeh-Bonn, S. Rafaï, A. Azouni, and D. Bonn.  Evaporating droplets.  Journal of Fluid Mechanics, 549:307–313, 2006.
  * [39] B.M. Weon, J.H. Je, and C. Poulard.  Convection-enhanced water evaporation.  AIP Advances, 1:012102, 2011.
  * [40] P.L. Kelly-Zion, C.J. Pursell, S. Vaidya, and J. Batra.  Evaporation of sessile drops under combined diffusion and natural convection.  Colloids and Surfaces A:, 381:31–36, 2011.
  * [41] P.L. Kelly-Zion, C.J. Pursell, N. Hasbamrer, B. Cardozo, K. Gaughan, and K. Nickels.  Vapor distribution above an evaporating sessile drop.  International Journal of Heat and Mass Transfer, 65:165 – 172, 2013\.
  * [42] F. Carle, B. Sobac, and D. Brutin.  Experimental evidence of the atmospheric convective transport contribution to sessile droplet evaporation.  Applied Physics Letters, 102, 2013.
  * [43] S. Somasundaram, T. N. C. Anand, and S. Bakshi.  Evaporation-induced flow around a pendant droplet and its influence on evaporation.  Physics of Fluids, 27(11), 2015.
  * [44] Z. Rotem and L. Claassen.  Natural convection above unconfined horizontal surfaces.  Journal of Fluid Mechanics, 39:173–192, 10 1969.
  * [45] S. C. Traugott.  Buoyant convection from horizontal surfaces at constant pressure.  Journal of Fluid Mechanics, 68:609–624, 4 1975.
  * [46] S. Dehaeck, A. Rednikov, and P. Colinet.  Vapor-based interferometric measurement of local evaporation rate and interfacial temperature of evaporating droplets.  Langmuir, 30(8):2002–2008, 2014.
  * [47] J.H. Merkin.  Free convection above a heated horizontal circular disk.  Zeitschrift für angewandte Mathematik und Physik, 34(5):596–608, 1983.
  * [48] G.K. Batchelor.  An Introduction to Fluid Dynamics.  Cambridge University Press, 2000.
  * [49] A. Edelstein, N. Amodaj, K. Hoover, R. Vale, and N. Stuurman.  Computer Control of Microscopes Using Micromanager.  John Wiley & Sons, Inc., 2010.
  * [50] G. Sudre.  Tunable Adhesion of Hydrogels.  PhD thesis, Université Pierre et Marie Curie - Paris VI, 2011.
  * [51] D. Allan, T. Caswell, and N. Keim.  Trackpy v0.2.  2014\.  doi:10.5281/zenodo.9971.
  * [52] T. Oliphant.  Python for scientific computing.  Computing in Science & Engineering, 9:10–20, 2007.
  * [53] K.J Millman and M. Aivazis.  Python for scientists and engineers.  Computing in Science & Engineering, 13:9–12, 2011.
  * [54] S. van der Walt, S.C. Colbert, and G. Varoquaux.  The NumPy array: A structure for efficient numerical computation.  Computing in Science & Engineering, 13:22–30, 2011.
  * [55] W. McKinney.  Data structures for statistical computing in Python.  In Proceedings of the 9th Python in Science Conference, 2010.
  * [56] S. van der Walt, J.L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J.D. Warner, N. Yager, E. Gouillart, and T. Yu.  scikit-image: Image processing in Python.  PeerJ, 2:e453, 2014.


††thanks: Equal contribution††thanks: Equal contribution

# Learning Lattice Quantum Field Theories with Equivariant Continuous Flows

Mathis Gerdes  Institute of Physics, University of Amsterdam, the Netherlands
Pim de Haan  QUVA Lab, University of Amsterdam  Qualcomm AI Research, Qualcomm
Technologies Netherlands B.V.  Qualcomm AI Research is an initiative of
Qualcomm Technologies, Inc.  Corrado Rainone  Qualcomm AI Research, Qualcomm
Technologies Netherlands B.V.  Qualcomm AI Research is an initiative of
Qualcomm Technologies, Inc.  Roberto Bondesan  Qualcomm AI Research, Qualcomm
Technologies Netherlands B.V.  Qualcomm AI Research is an initiative of
Qualcomm Technologies, Inc.  Miranda C. N. Cheng  Institute of Physics,
University of Amsterdam, the Netherlands  Korteweg-de Vries Institute for
Mathematics, University of Amsterdam, the Netherlands  Institute for
Mathematics, Academia Sinica, Taipei, Taiwan

###### Abstract

We propose a novel machine learning method for sampling from the high-
dimensional probability distributions of Lattice Quantum Field Theories.
Instead of the deep architectures used so far for this task, our proposal is
based on a single neural ODE layer and incorporates the full symmetries of the
problem. We test our model on the $\phi^{4}$ theory, showing that it
systematically outperforms previously proposed flow-based methods in sampling
efficiency, and the improvement is especially pronounced for larger lattices.
Compared to the previous baseline model, we improve a key metric, the
effective sample size, from 1% to 91% on a lattice of size $32\times 32$. We
also demonstrate that our model can successfully learn a continuous family of
theories at once, and the results of learning can be transferred to larger
lattices. Such generalization capacities further accentuate the potential
advantages of machine learning methods compared to traditional MCMC-based
methods.

## Introduction.

Lattice field theory (LFT) remains the only robust and universal method to
extract physical observables from a non-perturbative quantum field theory
(QFT). As a result, the applications of such lattice computations are
omnipresent in physics. Calculations for LFTs have traditionally relied on
Markov Chain Monte Carlo (MCMC) methods. A well-known challenge these sampling
methods face is the phenomenon of critical slowing down [1]. When moving
towards the critical point or the continuum limit, consecutively generated
samples become increasingly correlated as measured by a rapidly increasing
autocorrelation time (time between statistically independent samples). In
practice, this computational inefficiency constitutes a bottleneck for
achieving high-accuracy results. While efficient algorithms exist for certain
statistical mechanical models — e.g. the Swendsen-Wang algorithm for the Ising
and Potts models [2] — a general approach is still lacking. An additional
challenge arises from the fact that the entire sampling procedure typically
needs to be initiated anew each time the theory is deformed, for instance by
changing the UV cutoff (lattice spacing) or the coupling constants.

Figure 1: The effective sample sizes for lattice sizes $L\times L$ of our
model, compared to the real NVP baseline. Dashed lines show the learning
trajectory of transfer learning from $L=20$ to $L=64$.

Rapid progress in machine learning provides new tools to tackle this complex
sampling problem, see [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20] for earlier related works. Flow-based models, in particular, offer
a promising new framework for sampling which can mitigate the problem of
critical slowing down and amortizes the costs of sampling by learning to
approximate the target distribution. Moreover, since they estimate the
likelihood of configurations they can be used to estimate some observables
which are challenging to estimate with standard MCMC methods [21, 22]. In the
flow framework, a parametrized function $f_{\theta}$ is learned, which maps
from an auxiliary “latent” space to the space of field configurations. It is
optimized such that an easy-to-sample distribution $\rho$ (e.g. an independent
Gaussian) gets transformed, or pushed forward, to approximate the complex
distribution given by the action of the theory [23]. The capability to
incorporate physical symmetries of either the local or global type into the
flow can moreover circumvent difficulties of sampling caused by the
degeneracies of physical configurations which are related by symmetry actions
[24].

The idea of applying a flow model to lattice field theory has been tested in
[14] using the so-called real NVP model, and proof-of-concept experiments have
been successfully carried out for small lattices. However, significant hurdles
remain before a flow-based model can be applied in a physically interesting
setup. For instance, it has been argued that these models require a training
time that rises exponentially with the number of lattice sites [25]. In
practice, it becomes infeasible to attain high acceptance rates on large
lattices, while these are precisely the physically interesting regime.

To go beyond proof-of-concept studies, we devise a novel flow model111An
implementation can be found at https://github.com/mathisgerdes/continuous-
flow-lft. which improves the performance in the following three aspects, when
compared to the state-of-the-art baseline model [25, 14, 18]: 1) scalability,
2) efficiency, and 3) symmetries. In more details, we contribute the
following. 1) We propose a continuous flow model, where the inverse
trivializing map $f_{\theta}$ is given as the solution to an ordinary
differential equation (ODE). When applied to the two-dimensional $\phi^{4}$
theory, our model dramatically improves a key metric, the effective sample
size (ESS), to $91\%$ at $L^{2}=32^{2}$ after about 10 hours of training as
shown in Figure 1, in contrast with to the $1.4\%$ achieved by the real NVP
model after approximately 2 weeks of training. Going further beyond what has
been reported so far [25], we also obtain good results on lattices with
$L^{2}=64^{2}$ vertices. Moreover, as opposed to the real NVP flow tested in
[25], Figure 3 demonstrates that our continuous flow model is sample efficient
in training and in particular does not display an exponentially growing
training cost when lattice size increases. 2) We demonstrate the ability of
our model to learn many theories at once in the following sense. First, we
show that training efficiency can be gained by upscaling the network learned
for lattice size $L$ to a larger lattice of size $L^{\prime}$. See Figure 1.
Second, we train a family of flows $f_{\theta}^{(\psi)}$ parametrized by
theory parameters $\psi$ within a continuous domain, as used for Figure 3.
Once the training is done, sampling is extremely cheap and no particular
slowing down is encountered near the critical point. 3) In contrast to the
real NVP architecture (cf. [26]), our model is equivariant to all symmetries
of the LFT, including the full geometric symmetries of the lattice and the
global symmetry of the scalar theory. If not built into the architecture,
these symmetries are only approximately learned by the model, as exemplified
in Figure 6 in the appendix.

Figure 2:  Samples used during training until the values of ESS is attained,
for lattice sizes $L=6,12,20,32,64$.

Figure 3: Susceptibility estimated by MCMC chains of length $10^{6}$, rescaled
using the Ising model critical exponents $\nu=1$, $\gamma=7/4$ to account for
finite-size scaling. The used value $\lambda_{c}=4.25$ of the critical
coupling was chosen based on previous numerical results [27].

## Scalar Quantum Field Theory.

Our goal is to improve the MCMC sampling performance for quantum field
theories with scalar fields which possess non-trivial symmetry properties. We
will now consider a scalar field theory on a two-dimensional periodic square
lattice $V_{L}\cong({\mathbb{Z}}/L{\mathbb{Z}})^{2}$. In particular, the field
configuration is a real function $\phi:V_{L}\to{\mathbb{R}}$ on the vertex set
$V_{L}$, and we denote by $\phi_{x}$ the value at the vertex $x\in V_{L}$. The
theory is described by a probability density $p(\phi)=\exp(-S(\phi))/Z$ with
action

$S(\phi)=\sum_{x,y\in V_{L}}\phi_{x}\Delta_{xy}\phi_{y}+\sum_{x\in
V_{L}}m^{2}\phi_{x}^{2}+\sum_{x\in V_{L}}V(\phi_{x})\,,$ (1)

where $V:\mathbb{R}\to\mathbb{R}$ is a potential function, $\Delta$ is the
Laplacian matrix of the periodic square lattice, and $m^{2}$ is a numerical
parameter. The partition function is defined as $Z=\int\prod_{x\in
V_{L}}d\phi_{x}e^{-S(\phi)}$. For polynomials $V(\phi)$ of degree greater than
two, the normalization factor $Z$ and the moments of $p(\phi)$ are not known
analytically, and are typically estimated numerically. The statistical
correlations between spatially-separated degrees of freedom of the theory are
measured by the _correlation length_ $\xi$. As a result, it also controls the
difficulty of directly sampling from $p(\phi)$.

We focus on the case where the potential is invariant to all spatial
symmetries of the lattice and has a $\mathbb{Z}_{2}$ symmetry
$V(\phi)=V(-\phi)$. This family of theories is also known as the Landau-
Ginzburg model, and possesses a rich phase diagram with multi-critical points
described by the unitary minimal model conformal field theories [28].

## Samplers Based on Normalizing Flows.

To correct for the approximate nature of the learned map, a Metropolis-
Hastings step is applied to reject or accept the pushed-forward samples. This
guarantees asymptotic exactness, while, assuming the acceptance rate is
sufficiently high, the computational cost is not significantly increased.
Recall that the Independent Metropolis-Hastings algorithm is a flexible method
to sample from a target density $p$ given a proposal distribution $q$ [29]. At
each time $i$, a proposed sample $\phi^{\prime}$ drawn from the proposed
distribution $q$ is accepted with probability
$\min\left(1,{q(\phi^{(i-1)})}{p(\phi^{\prime})}/{q(\phi^{\prime})}{p(\phi^{(i-1)})}\right)$,
which depends on the previous sample $\phi^{(i-1)}$. Accepting the proposal
means setting $\phi^{(i)}=\phi^{\prime}$, else we repeat the previous value
$\phi^{(i)}=\phi^{(i-1)}$. Traditionally, the distribution $q$ is
painstakingly handcrafted to maximize the acceptance rate. In contrast, the ML
approach [30] to this problem is to learn a proposal distribution $q(\phi)$,
for example by using a normalizing flow [14]. As mentioned before, in this
framework $q(\phi)$ is given by the push-forward under a learned map
$f_{\theta}$ of some simple distribution $\rho$ such as a normal distribution.
Since the proposals are generated independently, rejections in the ensuing
Metropolis-Hastings step are the only source of autocorrelation.

When samples from $p$ are available, $q$ can be matched to $p$ by optimizing
the parameters $\theta$ to maximize the log-likelihood of those samples:
$\mathbb{E}_{\phi\sim p}[\log q_{\theta}(\phi)]$. However, obtaining such
samples requires expensive MCMC sampling which is precisely what we would like
to avoid. Instead, we aim to minimize the reverse Kullback-Leibler (KL)
divergence $\text{KL}(q|p)$, also known as the reverse relative entropy, which
uses samples from the model distribution $q$:

$\begin{split}\text{KL}(q|p)&=\mathbb{E}_{\phi\sim
q}[\log\tfrac{q(\phi)}{p(\phi)}]\\\ &=\mathbb{E}_{z\sim\rho}[\log
q(f_{\theta}(z))+S(f_{\theta}(z))]+\log Z\,.\end{split}$ (2)

In the second equality, we have used the flow $f_{\theta}$ to replace samples
from $q$ by $f_{\theta}(z)$, where $z$ denote samples of $\rho$. Note that the
last term does not contribute to the gradient.

In our continuous normalizing flow model, we transform a sample $z\sim\rho$
with an invertible map
$f_{\theta}:\mathbb{R}^{L^{2}}\to\mathbb{R}^{L^{2}},z\mapsto\phi$, defined as
the solution to a neural ODE [31] for time $t\in[0,T]$:

$\frac{\mathrm{d}\phi(t)_{x}}{\mathrm{d}t}=g_{\theta}(\phi(t),t)_{x}~{}{\rm
with}~{}z\equiv\phi(0),~{}\phi\equiv\phi(T)\,.$ (3)

Here, the vector field $g_{\theta}(\phi(t),t)_{x}$ is a neural network with
weights $\theta$. The probability $q(\phi)\equiv p(\phi(T))$ of samples
transformed by the flow, necessary for computing the KL divergence, can be
obtained by solving another ODE [31] which contains the divergence of
$g_{\theta}$:

$\begin{split}&\frac{\mathrm{d}\log
p(\phi(t))}{\mathrm{d}t}=-(\nabla_{\phi}\cdot
g_{\theta})(\phi(t),t)~{}\end{split}$ (4)

with the boundary condition $p(\phi(0))=\rho(z)$. For the network
architectures of the vector field $g_{\theta}$ described below, the divergence
can be computed analytically. The gradients of the KL divergence are computed
with another ODE derived using the adjoint sensitivity method [31].

As mentioned in [32], if the vector field $g_{\theta}$ is equivariant to the
symmetries of the theory and the chosen latent space distribution is invariant
under the symmetry action, then the resulting distribution on $\phi$ is
automatically invariant. In other words, for a group element $h\in G$ acting
on the space of field configurations $\mathbb{R}^{L^{2}}$, given
$\rho(hz)=\rho(z)$ and $h(g_{\theta}(\phi(t),t))=g_{\theta}(h\phi(t),t)$ it
follows that $p(h\phi(t))=p(\phi(t))$. In the next section we show how to
construct a fully equivariant vector field $g_{\theta}$.

## Method.

Inspired by equivariant flows used for molecular modelling [32], we propose to
construct a time-dependent vector field for the neural ODE featuring pair-wise
interactions between the lattice sites and a tensor contraction with a chosen
set of basis functions:

$\displaystyle\frac{\mathrm{d}\phi(t)_{x}}{\mathrm{d}t}=\sum_{y,d,f}W_{xydf}K(t)_{d}H(\phi(t)_{y})_{f}\,.$
(5)

With chosen dimensions $D$ and $F$, the sums of $d$ and $f$ are taken over
$1,\dots,D$ and $1,\dots,F$, respectively. Here,
$H:\mathbb{R}\to\mathbb{R}^{F}$ is a basis expansion function for the local
field values $\phi_{y}$, $K:[0,T]\to\mathbb{R}^{D}$ is our chosen time kernel,
and $W\in\mathbb{R}^{L^{2}\times L^{2}\times D\times F}$ is a learnable weight
tensor.

Our choice for the basis function is $H(\phi)_{1}=\phi$ and
$H(\phi)_{f}=\sin(\omega_{f}\phi)$ for $f\in\\{2,...,F\\}$, namely a linear
term combined with a sine expansion where $\omega\in\mathbb{R}^{F-1}$ are
learnable frequencies, inspired by Fourier Features [33]. This choice of basis
functions is made in order to respect the global symmetry $\phi\mapsto-\phi$,
as is done in [21]. In particular, the cosine function is prohibited by the
symmetry. The functional dependence of (5) on $\phi$ is also chosen so that
the divergence of the vector field, necessary for computing the density of the
resulting distribution $q$, can be easily computed analytically.

The time basis functions $K:[0,T]\to\mathbb{R}^{D}$ are chosen to be the first
$D$ terms of a Fourier expansion on the interval $[0,T]$. This is chosen to
allow for simple time-dependent dynamics.

The learnable weight tensor $W$ is initialized to $0$, so that the initial
flow is the identity map. Linear constraints on the tensor imposed by
equivariance with respect to spatial symmetries reduces the number of
independent entries of the weight tensor. The spatial symmetry group of the
periodic lattice $(\mathbb{Z}/L\mathbb{Z})^{\times 2}$ is $G=C_{L}^{2}\rtimes
D_{4}$, the semi-direct product of two cyclic groups $C_{L}$ acting as
translations, and the Dihedral group $D_{4}$ generated by 90-degree planar
rotations and the mirror reflection. $G$-equivariance dictates that $W_{xydf}$
for a given $d$ and $f$ only depends on the orbit of the pair $(x,y)$.
Specifically, for any element $g$ of the symmetry group $G$ acting on a pair
of lattice points as $(x,y)\mapsto(gx,gy)$, we have $W_{gx,gy,df}=W_{xydf}$
for all values of $x$, $y$, $d$, and $f$. Explicitly, since we can use
$C_{L}^{2}$-equivalences to move $x$ to the (arbitrarily chosen) origin given
any pair $(x,y)$, and for most of these pairs there are eight $y^{\prime}$
satisfying $(x,y)\sim(x,y^{\prime})$ under $D_{4}$ equivalance, the number of
free parameters per $d$ and $f$ grows like $L^{2}/8+O(L)$ as opposed to
$L^{4}$ without symmetry constraints.

Finally, we observe that the following factorization of the matrix leads to
better training results

$W_{xydf}=\sum_{d^{\prime}f^{\prime}}\tilde{W}_{xyd^{\prime}f^{\prime}}W^{K}_{d^{\prime}d}W^{H}_{f^{\prime}f}$
(6)

with matrices $W^{K}\in\mathbb{R}^{D^{\prime}\times
D},W^{H}\in\mathbb{R}^{F^{\prime}\times F}$. The bond dimensions $D^{\prime}$
and $F^{\prime}$ are hyperparameters of the model.

## A Theory-Conditional Model.

A small change in the parameters of the theory typically leads to a small
deformation of its associated distribution, which is reflected in our
observation that the network trained for a specific theory still results in a
relatively high value of ESS even when applied to a nearby theory with
slightly altered parameters. This motivates the extended network architecture
and the transfer learning between lattice sizes, as discussed in the current
and the following paragraphs respectively.

Instead of a single action $S$, here we consider a set of theories ${\cal M}$
on the lattice. Each theory $\psi\in{\cal M}$ is specified by the action
$S_{\psi}:\mathbb{R}^{L^{2}}\to\mathbb{R}$, and we can then learn a single
neural network to generate samples for all these theories. We do this via an
additional basis expansion function $J:{\cal M}\to\mathbb{R}^{A}$ and a new
index to each of the terms in the factorized $W$ tensor: instead of (6), we
let

$W_{xydf}(\psi)=\sum_{d^{\prime}f^{\prime}abc}\tilde{W}_{xyd^{\prime}f^{\prime}a}J(\psi)_{a}W^{K}_{d^{\prime}db}J(\psi)_{b}W^{H}_{f^{\prime}fc}J(\psi)_{c}\,.$
(7)

By choosing a distribution $r(\psi)$ over the theories, the parameters of the
model are optimized by maximising the reverse KL divergence averaged over the
theories: $\mathbb{E}_{\psi\sim r}\text{KL}(q_{\psi}|p_{\psi})$, where
$q_{\psi}$ is the push-forward distribution using to the map solving the ODE
(5) with the conditional weight matrix (7) and $p_{\psi}$ is the Boltzmann
distribution of action $S_{\psi}$.

Table 1: ESS, MCMC acceptance rate and observables $\xi$ and $\chi_{2}^{(1)}$ for the flow models trained with different lattice sizes $L$ and the corresponding $\lambda$ tuned to give $L/\xi=4$ [27].  $L$ | $\lambda$ | ESS | MCMC | $L/\xi$ | $\chi_{2}^{(1)}$
---|---|---|---|---|---
6 | $6.975$ | $99.759(7)$ | $97.503(4)$ | $3.968(5)$ | $1.064(2)$
12 | $5.276$ | $98.54(2)$ | $93.367(7)$ | $3.981(5)$ | $4.129(6)$
20 | $4.807$ | $96.63(4)$ | $89.77(1)$ | $4.039(5)$ | $10.59(2)$
32 | $4.572$ | $91.07(8)$ | $83.15(2)$ | $4.031(5)$ | $25.53(4)$
64 | $4.398$ | $66(5)$ | $64.96(2)$ | $4.012(4)$ | $91.6(2)$

## Transfer-Learning Between Lattice Sizes.

The parameters $W$ of a model that is trained on a square lattice with $L^{2}$
sites may be used as the initial parameters $W^{\prime}$ for training a model
on a larger lattice with length $L^{\prime}$. This is an example of transfer
learning or curriculum learning [34] and may lead to a speedup compared to
training directly on the lattice of length $L^{\prime}>L$. Roughly speaking,
to transfer we embed the kernel $W$ into the kernel $W^{\prime}$ for the
larger lattice. We refer the reader to the appendix for details. The
frequencies $\omega_{f}$ are kept fixed during the transfer. After
initializing the parameters this way, the model can be further trained on the
lattice of length $L^{\prime}$.

## Experiments.

For our experiments, we consider the $\phi^{4}$ theory, which has an action as
in equation (1) with potential $V(\phi)=\mathit{\lambda}\sum_{x}\phi_{x}^{4}$,
specified by the coupling constant $\lambda$. In the first experiment, we
evaluate the ability of our model as defined in equation (5) to learn to
generate samples from a single theory, i.e. for a fixed value of $\lambda$. We
test this on lattices with length $L$ varying between 6 and 64, where the
theory parameter $m^{2}$ is held fixed to be $m^{2}=-4$, and $\lambda$ is
tuned so that the correlation length $\xi$ approximately equals to $1/4$ of
the lattice side length. We choose $F=50$ for the number of field basis
functions, $D=21$ for the number of time kernels, and $F^{\prime}=20$,
$D^{\prime}=20$ for the bond dimensions of the respective factorization
matrices. The parameters are optimized by gradient descent, where the
expectation in the loss of equation (2) is approximated by averaging over
$256$ generated samples. More algorithmic details can be found in the
appendix.

To assess model quality, we use two different but related metrics, the ESS and
the MCMC acceptance rate. We remind the reader that the ESS is a quantity
closely related to the acceptance rate. It controls the variance of the mean
of $N$ correlated random variables, which scales as $(\text{ESS}\times
N)^{-1/2}$. In particular, an ESS of $100\%$ indicates perfectly independent
samples and is attained when the target and the proposal distribution
coincide. In Figure 1 we report the ESS values against the training time,
computed as a moving average over 100 training steps, with the real NVP
baseline for comparison, which is based on [18] and was trained for two weeks.
For each solid line, the corresponding coupling $\lambda$ is chosen such that
the correlation length is approximately $L/\xi=4$, as shown in Table 1. Our
model attains much larger values of the ESS and acceptance rates in much
shorter times, for all sizes $L$, despite the fact that each of its training
steps involves integration of an ODE. Table 1 shows the final ESS and MCMC
acceptance rates as well as $L/\xi$ and the two-point susceptibility
$\chi^{(1)}_{2}$ computed with MCMC using the trained flows. Note that the
correlation length computed with our flow-assisted MCMC is indeed comparable
with that computed using the traditional MCMC method. The description of how
$\xi$, ESS, and two-point susceptibility are estimated is recorded in the
appendix. In addition, we demonstrate that transfer learning as described
above can be used to reduce the training time for larger lattices. This is
shown by dashed lines in Figure 1. All three dashed lines represent a single
run of transfer learning between lattice sizes with the target $L=64$ and the
correspondingly fixed value $\lambda=4.398$. Initially, we train on a theory
with $L=20$ and the $\lambda$ from Table 1 listed for $L=64$. After the ESS
approximately saturates during training, the kernel is scaled to size
$32\times 32$ and training is continued with the same value of the coupling
constant $\lambda$ but at $L=32$. Finally, this is repeated to move to $L=64$.
We observe that, when using transfer learning, it takes a much shorter time to
attain a similar value of the ESS when compared to training directly with
$L=64$.

In the final experiment, we train the theory-conditional models as described
above and summarized in equation (7) over a range of coupling constants from
$\lambda=4$ to $\lambda=6$. For the model hyperparameters, we now choose
$F=300$, $F^{\prime}=20$ for the field basis functions, $D=21$,
$D^{\prime}=20$ for the time kernels, and $A=50$ for the number of $\lambda$
basis expansion functions. The latter are chosen to be of the form

$\displaystyle J(\lambda)_{i}=\frac{\exp(-w\cdot(\lambda-
c_{i})^{2})}{\sum_{j}\exp(-w\cdot(\lambda-c_{j})^{2})}\,,$ (8)

where $w$ is trainable and the centers $c_{i}$ are fixed and uniformly spaced
in the range of $\lambda$. In order to speed up training, we initially train
at $L=12$ and then transfer to $L=20$ and finally $L=32$ as described above.
For each size, training was continued until the ESS approximately stabilized.
Our model achieves high MCMC acceptance rates between $60\%$ and $90\%$ over
the whole range of considered theory space. The robustness of the model was
tested by performing the training multiple times with different random seeds,
which renders similar performance each time as shown in more detail in the
appendix.

To demonstrate the physical usefulness of this approach, Figure 3 shows the
susceptibility, computed using the estimator $\chi_{1}^{(2)}$ via MCMC with
the trained model as proposal distribution. With the finite-size scaling taken
into account, our numerical results display the expected collapse of the
curves for different lattice sizes. To obtain these results, a theory
conditional flow was trained three times for each lattice size with different
random seeds over the range of $\lambda$. Using each of the three trained
models as proposal, $\chi_{1}^{(2)}$ was estimated by MCMC. The values shown
are the mean while the corresponding standard deviations averaged over
$\lambda$ are $0.3\%$, $0.5\%$, $1.6\%$ for $L=12$, $20$, $32$, respectively,
and are thus too small to be visible in the figure.

## Conclusions.

In this work we propose a fully equivariant continuous flow model as an
effective tool for sampling in lattice field theory. Our experiments indicate
that such neurally-augmented sampling is a viable method even for large
lattices. Moreover, we have shown that the learned parameters can be
meaningfully transferred to nearby theories with different correlation
lengths, coupling constants, and lattice sizes, even across the critical
point. This makes clear that our approach can be employed for tasks where
computations need to be done for families of theories, such as phase
transition detection and parameter scans. Finally, we have demonstrated that
our flow model can easily incorporate the geometric and global symmetries of
the physical system, resulting in in a visible computational advantage which
increases with the lattice size. These developments make our flow-model
approach a practical tool ready to be employed in lattice field theory
computations.

We will end with a comment on a potential interpretation of the learned ODE
flow. Note that the flow of the density described by the ODE (4) can readily
be interpreted as a flow in the space ${\cal M}$ of coupling constants by
identifying the probability distribution $p(t)$ as the Boltzmann distribution
of a Hamiltonian $H_{t}$. As such, it is tantalizing to try to understand the
flow in terms of the physical renormalization group of the theory, which also
renders a flow in ${\cal M}$. We will report on this in the future.

## Acknowledgements.

We wish to thank Max Welling for inspiring discussions, Vassilis Anagiannis
for his contribution in the early stages of this work, and Jonas Köhler for
his helpful suggestions. We also thank Joe Marsh Rossney, Luigi Del Debbio and
Uroš Seljak for helpful correspondence.

Appendix

##  Appendix A Observables

In the MCMC setting, observables are estimated using a set of $N$ samples
$\quantity{\phi^{(i)}}_{i=1}^{N}$ from the target distribution. These are
generated here by applying the Metropolis-Hasting step on samples proposed by
the trained continuous normalizing flows. Here we consider scalar fields on a
two-dimensional periodic square lattice of length $L$,
$V_{L}\cong({\mathbb{Z}}/L{\mathbb{Z}})^{2}$.

The effective sample size (ESS) is computed as in [18]:

$\displaystyle\text{ESS}=\frac{\quantity[\frac{1}{N}\sum_{i=1}^{N}p\quantity(\phi^{(i)})/q\quantity(\phi^{(i)})]^{2}}{\frac{1}{N}\sum_{i=1}^{N}\quantity[p(\phi^{(i)})/q\quantity(\phi^{(i)})]^{2}}\,.$
(9)

To compute the correlation length $\xi$ we first introduce the estimator for
the connected two-point Greens function

$\displaystyle G(x)=\frac{1}{L^{2}}\sum_{y\in
V_{L}}\frac{1}{N}\sum_{i=1}^{N}\quantity[\phi_{y}^{(i)}\phi_{x+y}^{(i)}-\phi_{y}^{(i)}\frac{1}{N}\sum_{j=1}^{N}\phi_{x+y}^{(j)}]\,.$
(10)

Summing over one of the lattice dimensions, we have

$\displaystyle G_{s}(x_{2})=\sum_{x_{1}=1}^{L}G(x_{1},x_{2})\,.$ (11)

The correlation length can then be estimated via

$\displaystyle\frac{1}{\xi}=\frac{1}{L-1}\sum_{x_{2}=1}^{L-1}\mathrm{arcosh}\quantity(\frac{G_{s}(x_{2}+1)+G_{s}(x_{2}-1)}{2G_{s}(x_{2})})\,.$
(12)

The two-point susceptibility is estimated as

$\displaystyle\chi^{(1)}_{2}=\sum_{x\in
V_{L}}G(x)=L^{2}\quantity(\frac{1}{N}\sum_{i=1}^{N}\overline{\phi^{(i)}}^{2}-\quantity(\frac{1}{N}\sum_{i=1}^{N}\overline{\phi^{(i)}})^{2})\,,\quad\text{where}\quad\overline{\phi^{(i)}}=\frac{1}{L^{2}}\sum_{x\in
V_{L}}\phi^{(i)}_{x}\,.$ (13)

In the broken phase, the distribution of $\overline{\phi}$ becomes bimodal. To
address this, we introduce an alternative estimator for the two-point
susceptibility [35]

$\displaystyle\chi^{(2)}_{2}=L^{2}\quantity(\frac{1}{N}\sum_{i=1}^{N}\overline{\phi^{(i)}}^{2}-\quantity(\frac{1}{N}\sum_{i=1}^{N}\quantity|\overline{\phi^{(i)}}|)^{2})\,.$
(14)

In the table of estimated observables at different lattice sizes, the ESS was
computed $100$ times using different random seeds, each time using $500$
samples. MCMC of length $10^{6}$ was run $10$ times with different random
seeds. In both cases, the average and its estimated error are shown. The
observables were computed with a single MCMC run of length $10^{6}$ and errors
estimated using bootstrap drawing $100$ times and with bin size 4. A burn-in
phase of length $10^{5}$ was discarded.

##  Appendix B Experimental details

## Network architecture.

The width factor in the chosen form of the $\lambda$ basis expansion functions
$J(\lambda)_{i}$ should not become negative during training. It is thus
parametrized by the trainable parameter $\tilde{w}$ as
$w=(A-1)\log(1+\exp(\tilde{w}))$, which is initialized to
$\tilde{w}=\log(e-1)$. The frequencies $\omega_{f}$ are initialised by random
uniform values between 0 and 5. Orthogonal initialization was used for the
matrices $W^{K}$ and $W^{H}$, while $\tilde{W}$ was initialized with zeros.

## Loss computation.

The KL-divergence between the network’s output distribution and the theory
distribution is computed by a Monte Carlo approximation. For the initial flow
model, $256$ samples are used in each training step to approximate the
integral. To train the theory-conditional flows, a sampling scheme for the
considered theory parameters must be chosen. Here, in each training step $8$
different values of $\lambda$ are sampled uniformly in the chosen range. For
each, the KL-divergence is computed by a Monte Carlo approximation to the
integral with $128$ samples. The loss is then the average of the $8$ different
KL-divergences. The training performance was not observed to strongly depend
on the chosen sample sizes. The above values are chosen as a compromise based
on available computational resources and are not the result of a comprehensive
hyperparameter search.

## Optimizer.

The flow networks are optimized by minimizing the KL-divergence using
stochastic gradient descent. For this, the Adam optimizer is used with a
learning rate of $0.005$, decaying exponentially by a factor of $0.01$ every
$8000$ steps. After transferring to a larger lattice, the initial learning
rate was reduced to $10^{-5}$. The decay parameters $\beta_{1}=.8$, and
$\beta_{2}=.9$ of the Adam optimizer were found to improve training speed. The
same values were used for all lattice sizes. Training efficiency can likely be
improved by tuning these hyperparameters.

## Training theory-conditional models.

Figure 4 shows the MCMC acceptance rate achieved by training a theory-
conditional model over a range of $\lambda$ couplings from $4$ to $6$. To
demonstrate the robustness of the model, it was independently trained three
times with different random seeds. In all runs a similar performance was
achieved after training.

Figure 4: Acceptance rates of MCMC chains of length $10^{6}$ for a theory-
conditional flow trained over a range of $\lambda$ for $L=32$.

## Symmetries.

Figure 6 shows the result of the ablation study we conducted by limiting the
symmetries manifestly preserved by the model. We observe a noticeable loss of
ESS during training. This becomes more significant as the lattice size is
increased. We also observed that models without the full symmetries built-in
are more prone to training instabilities when the lattice size increases.

To investigate the effect of having manifest equivariance compared to
approximate equivariance learned during training, we look at whether
symmetries are broken by the trained real NVP network. The result, shown in
Figure 6, shows the quality of the approximation when equivariance of the
network is not imposed.

Figure 5: Training performance of the continuous normalizing flow at $L=20$
with and without enforcing $D4$-symmetry. Shown are the mean and standard
deviation (shaded area) over $8$ training runs with different random seeds,
where the shown value of the ESS is the moving average over $100$ training
steps.

Figure 6: For 4 samples from an MCMC chain with $L=6$, we show 1) the true log
likelihood given by the action, 2) the model log likelihood of our model and
3) real NVP. Log likelihood of 4 samples of an MCMC chain with $L=6$. The
$y$-axis shows the likelihoods when the sample is transformed by all $8\times
6^{2}$ symmetries of the lattice.

## Transfer Learning.

To embed the kernel $W$ into the kernel $W^{\prime}$ for the larger lattice,
we consider a pair of points $(x,y)$ in $C_{L}^{2}$ and a pair of points
$(x^{\prime},y^{\prime})$in $C_{L^{\prime}}^{2}$ and identify $t=y-x\in
C_{L}^{2}$ with $t=[t_{1},t_{2}]$ with
$t_{i}\in\\{-\lfloor{\frac{L}{2}}\rfloor,\dots,\lfloor{\frac{L-1}{2}}\rfloor\\}$,
and similarly
$t^{\prime}=y^{\prime}-x^{\prime}=[t_{1}^{\prime},t_{2}^{\prime}]$ with
$t_{i}\in\\{-\lfloor{\frac{L^{\prime}}{2}}\rfloor,\dots,\lfloor{\frac{L^{\prime}-1}{2}}\rfloor\\}$.
To transfer, we set as initial values
$W^{\prime}_{x^{\prime}y^{\prime}af}=W_{xyaf}N_{L}(t)/N_{L^{\prime}}(t^{\prime})$
if $t=t^{\prime}$, and zero otherwise. The scaling factor
$N_{L}(t)=|\\{RtR^{-1}\in C_{L}^{2}\mid R\in D_{4}\\}|$, where $RtR^{-1}$
denotes the action of $R$ on $t$, is the number of sites reached by
rotating/mirroring $y$ around $x$ in a periodic lattice of length $L$. It
accounts for the fact that weights may be used for more pairs in the larger
lattice, as illustrated in Figure 7.

Figure 7: Illustration of transferring weights $W$ between lattices of size
$L=4$ and $L^{\prime}=6$.

## References

  * Binder and Heermann [2019] K. Binder and D. W. Heermann,  _Monte Carlo Simulation in Statistical Physics: An Introduction_ (Springer, 2019).
  * Sokal [1997] A. Sokal, Monte carlo methods in statistical mechanics: Foundations and new algorithms, in  _Functional Integration: Basics and Applications_, edited by C. DeWitt-Morette, P. Cartier, and A. Folacci (Springer US, Boston, MA, 1997) pp. 131–192.
  * Tomiya and Nagai [2021] A. Tomiya and Y. Nagai, Gauge covariant neural network for 4 dimensional non-abelian gauge theory (2021), arXiv:2103.11965 [hep-lat] .
  * Wetzel and Scherzer [2017] S. J. Wetzel and M. Scherzer, Phys. Rev. B 96, 184410 (2017).
  * Blücher  _et al._ [2020] S. Blücher, L. Kades, J. M. Pawlowski, N. Strodthoff, and J. M. Urban, Phys. Rev. D 101, 094507 (2020).
  * Boyda  _et al._ [2021a] D. L. Boyda, M. N. Chernodub, N. V. Gerasimeniuk, V. A. Goy, S. D. Liubimov, and A. V. Molochkov, Phys. Rev. D 103, 014509 (2021a).
  * Yoon  _et al._ [2019] B. Yoon, T. Bhattacharya, and R. Gupta, Phys. Rev. D 100, 014504 (2019).
  * Zhang  _et al._ [2020] R. Zhang, Z. Fan, R. Li, H.-W. Lin, and B. Yoon, Phys. Rev. D 101, 034516 (2020).
  * Matsumoto  _et al._ [2020] T. Matsumoto, M. Kitazawa, and Y. Kohno, Progress of Theoretical and Experimental Physics 2021, 10.1093/ptep/ptaa138 (2020), 023D01.
  * Tanaka and Tomiya [2017] A. Tanaka and A. Tomiya, Towards reduction of autocorrelation in hmc by machine learning (2017).
  * Zhou  _et al._ [2019] K. Zhou, G. Endrődi, L.-G. Pang, and H. Stöcker, Phys. Rev. D 100, 011501 (2019).
  * Pawlowski and Urban [2020] J. M. Pawlowski and J. M. Urban, Mach. Learn. Sci. Tech. 1, 045011 (2020), arXiv:1811.03533 [hep-lat] .
  * Nagai  _et al._ [2020] Y. Nagai, A. Tanaka, and A. Tomiya, Self-learning monte-carlo for non-abelian gauge theory with dynamical fermions (2020).
  * Albergo  _et al._ [2019] M. Albergo, G. Kanwar, and P. Shanahan, Physical Review D 100, 034515 (2019).
  * Rezende  _et al._ [2020] D. J. Rezende, G. Papamakarios, S. Racaniere, M. Albergo, G. Kanwar, P. Shanahan, and K. Cranmer, in  _Proceedings of the 37th International Conference on Machine Learning_ , Proceedings of Machine Learning Research, Vol. 119, edited by H. D. III and A. Singh (PMLR, 2020) pp. 8083–8092.
  * Kanwar  _et al._ [2020] G. Kanwar, M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, S. Racanière, D. J. Rezende, and P. E. Shanahan, Phys. Rev. Lett. 125, 121601 (2020).
  * Boyda  _et al._ [2021b] D. Boyda, G. Kanwar, S. Racanière, D. J. Rezende, M. S. Albergo, K. Cranmer, D. C. Hackett, and P. E. Shanahan, Phys. Rev. D 103, 074504 (2021b).
  * Albergo  _et al._ [2021] M. S. Albergo, D. Boyda, D. C. Hackett, G. Kanwar, K. Cranmer, S. Racanière, D. J. Rezende, and P. E. Shanahan, Introduction to normalizing flows for lattice field theory (2021), arXiv:2101.08176 [hep-lat] .
  * Vaitl  _et al._ [2022] L. Vaitl, K. A. Nicoli, S. Nakajima, and P. Kessel, in  _International Conference on Machine Learning_ (PMLR, 2022) pp. 21945–21959.
  * Abbott  _et al._ [2022] R. Abbott, M. S. Albergo, A. Botev, D. Boyda, K. Cranmer, D. C. Hackett, G. Kanwar, A. G. D. G. Matthews, S. Racanière, A. Razavi, D. J. Rezende, F. Romero-López, P. E. Shanahan, and J. M. Urban, Sampling qcd field configurations with gauge-equivariant flow models (2022).
  * Nicoli  _et al._ [2021a] K. A. Nicoli, C. J. Anders, L. Funcke, T. Hartung, K. Jansen, P. Kessel, S. Nakajima, and P. Stornati, Phys. Rev. Lett. 126, 032001 (2021a), arXiv:2007.07115 [hep-lat] .
  * Nicoli  _et al._ [2020] K. A. Nicoli, S. Nakajima, N. Strodthoff, W. Samek, K.-R. Müller, and P. Kessel, Phys. Rev. E 101, 023304 (2020), arXiv:1910.13496 [cond-mat.stat-mech] .
  * Dinh  _et al._ [2017] L. Dinh, J. Sohl-Dickstein, and S. Bengio, Density estimation using real nvp (2017), arXiv:1605.08803 [cs.LG] .
  * Nicoli  _et al._ [2021b] K. A. Nicoli, C. Anders, L. Funcke, T. Hartung, K. Jansen, P. Kessel, S. Nakajima, and P. Stornati, arXiv preprint arXiv:2111.11303 (2021b).
  * Del Debbio  _et al._ [2021] L. Del Debbio, J. M. Rossney, and M. Wilson, Phys. Rev. D 104, 094507 (2021), arXiv:2105.12481 [hep-lat] .
  * Hackett  _et al._ [2021] D. C. Hackett, C.-C. Hsieh, M. S. Albergo, D. Boyda, J.-W. Chen, K.-F. Chen, K. Cranmer, G. Kanwar, and P. E. Shanahan, Flow-based sampling for multimodal distributions in lattice field theory (2021), arXiv:2107.00734 [hep-lat] .
  * Vierhaus [2010] I. Vierhaus,  _Simulation of $\varphi^{4}$ theory in the strong coupling expansion beyond the Ising limit_, Ph.D. thesis, Humboldt University (2010).
  * Mussardo [2010] G. Mussardo,  _Statistical Field Theory: An Introduction to Exactly Solved Models in Statistical Physics_, Oxford Graduate Texts (OUP Oxford, 2010).
  * Montvay and Münster [1997] I. Montvay and G. Münster,  _Quantum fields on a lattice_ (Cambridge University Press, 1997).
  * Noé  _et al._ [2019] F. Noé, S. Olsson, J. Köhler, and H. Wu, Science 365 (2019).
  * Chen  _et al._ [2018] T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, in  _Advances in neural information processing systems_ (2018) pp. 6571–6583.
  * Köhler  _et al._ [2020] J. Köhler, L. Klein, and F. Noe, in  _Proceedings of the 37th International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 119, edited by H. D. III and A. Singh (PMLR, 2020) pp. 5361–5370.
  * Tancik  _et al._ [2020] M. Tancik, P. P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. T. Barron, and R. Ng, arXiv preprint arXiv:2006.10739 (2020).
  * Bengio  _et al._ [2009] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, in  _Proceedings of the 26th annual international conference on machine learning_ (2009) pp. 41–48.
  * Sandvik [2010] A. W. Sandvik, AIP Conference Proceedings 1297, 135 (2010).