# Phase field models of active matter

Romain Mueller  Rudolf Peierls Centre for Theoretical Physics, University of
Oxford, Parks Road, OX1 3NP, Oxford, United Kingdom  Amin Doostmohammadi
<EMAIL_ADDRESS>The Niels Bohr Institute, University of Copenhagen,
Copenhagen, Denmark

###### Abstract

We present an overview of phase field modeling of active matter systems as a
tool for capturing various aspects of complex and active interfaces. We first
describe how interfaces between different phases are characterized in phase
field models and provide simple fundamental governing equations that describe
their evolution. For a simple model, we then show how physical properties of
the interface, such as surface tension and interface thickness, can be
recovered from these equations. We then explain how the phase field
formulation can be coupled to various active matter realizations and discuss
three particular examples of continuum biphasic active matter: active nematic-
isotropic interfaces, active matter in viscoelastic environments, and active
shells in fluid background. Finally, we describe how multiple phase fields can
be used to model active cellular monolayers and present a general framework
that can be applied to the study of tissue behaviour and collective
migration.111These are the lecture notes presented at the Initial Numerical
Training on Active Matter in a training school for the MSCA-ITN ActiveMatter,
(http://active-matter.eu)

##  I Introduction

Figure 1: Examples of natural and synthetic active interfaces. (a) Bacterial
competition during Pseudomonas ariginosa biofilm invasion. Fast and slow
moving bacterial strains are shown by yellow fluorescent protein (YFP) and
cyan fluorescent protein (CFP), respectively. Figure adapted from Meacock  _et
al._ (2020). (b) Epithelial cell progression during wound healing. Figure
adapted from Doostmohammadi  _et al._ (2015). (c) Synthetic active matter
composed of colloidal rollers forming fingering patterns at the progressive
front. Figure adapted from Driscoll  _et al._ (2017). (d) Self-deforming
active droplet composed of microtubule-motor protein mixtures in a fluid
background. Figure adapted from Keber  _et al._ (2014).

While understanding pattern formation within active matter has taken important
strides in the past decades Marchetti  _et al._ (2013a); Needleman and Dogic
(2017), the vast majority of active systems are characterized by dynamic
interfaces between an active phase and the surrounding medium, which can
itself be active or passive with significantly different properties (Fig. 1).
Whether it is a bacterial biofilm invading its surrounding environment
Burmølle  _et al._ (2006); Nadell  _et al._ (2015); Meacock  _et al._ (2020),
a cellular tissue closing a wound Brugués  _et al._ (2014), or a cluster of
human cancer cells migrating into extracellular matrices Friedl and Gilmour
(2009); Friedl  _et al._ (2012), the interaction between active matter and its
environment plays a crucial role in determining the behavior of prominent
physiological processes.

In addition, even within the active matter itself, material properties are not
always homogeneous and interfaces between active phases with different
physical and chemical properties can arise. Striking examples include
synthetic self-propelled particles that phase separate because of their
different self-propulsion Zhang  _et al._ (2020a), competition between
different phenotypes during bacterial collective invasion Meacock  _et al._
(2020), and cellular segregation in epithelial cells with distinct cell-cell
interactions Balasubramaniam  _et al._ (2020).

To correctly model such complex systems, it is then crucial to both capture
the dynamics of each phase independently but also to accurately track the
dynamic interface between them. This can represent, however, a major modelling
challenge. For example, there are many relevant biophysical scenarios where
active cellular materials interact with extracellular matrices with
viscoelastic properties Nadell  _et al._ (2015); Chaudhuri  _et al._ (2020).
Phase field methods provide a generic and versatile framework that is
adaptable to such complex situations and allow for the modelling of multiple
active or passive phases with their own complex dynamics and interactions.

##  II Phase field description of diffuse interfaces

The phase field method, originally introduced to model solidification
processes Hohenberg and Halperin (1977); Fix (1982); Langer , has since then
seen a wide range of applications, such as the description of fingering and
elastic surface instabilities, fluidization and crystallization in complex
media, and the modeling of soft vesicles, see Emmerich (2008) for a review.
The basic idea behind this method is to describe interfaces between two
distinct phases implicitly using an auxiliary field $\phi$ whose value
distinguishes the two phases. For instance, $\phi=0$ can denote the first
phase while $\phi=1$ the second, but this choice is arbitrary. The phase field
is then subject to a mixing free energy such that the balance between mixing
and demixing effects results in two distinct stable phases with a diffusive
interface of small, but finite, thickness. Under complete demixing, a sharp
interface separates these two phases, while the presence of a mixing free
energy results in diffusion of the order parameter in the interfacial region.
This method is extremely versatile and can be used to describe interfaces
between a variety of different phases, such as the solid and liquid phases of
a melting material, the nematic and homogeneous phases of a biphasic liquid
crystal, or simply the inside and outside of a vesicle or a cell. In what
follows, we present a simple phase field model with two phases and show how it
can be coupled to flow and stresses.

###  II.1 Mathematical basis

We present here the equations of a simple phase field model describing the
demixing of two distinct phases and show analytically how the model parameters
relate to the physical properties of the interface. For simplicity, we assume
that the phase field $\phi$ relaxes diffusively towards the minimum of the
free-energy $\mathcal{F}=\mathcal{F}[\phi,\nabla\phi]$ according to the
equation of motion

$\partial_{t}\phi=-\frac{\delta{\mathcal{F}}}{\delta\phi}=-\frac{\partial{\mathcal{F}}}{\partial\phi}+\nabla\cdot\frac{\partial{\mathcal{F}}}{\partial\nabla\phi}.$
(1)

This equation of motion, which was first introduced in Allen and Cahn (1972)
for the study of binary alloys, is arguably the simplest choice of evolution
for the phase field $\phi$ but others are possible and can be found throughout
the literature Emmerich (2008). A simple free-energy exhibiting two distinct
phases is the classical Cahn-Hiliard free energy first introduced in Cahn and
Hilliard (1958), which can be written as

$\mathcal{F}_{\text{CH}}=\int\mathrm{d}\vec{x}\,\left\\{\frac{A}{2}\>\phi^{2}\left(1-\phi\right)^{2}+\frac{\kappa}{2}(\nabla\phi)^{2}\right\\}.$
(2)

The first term is the demixing free energy and corresponds to a double well-
potential with minima at $\phi=0$ and $\phi=1$ corresponding to the first and
second phases, respectively, while the second term is the mixing part that
penalises gradients $\nabla\phi$ (see Fig. 2(a)). Note the phase field $\phi$
is not conserved in equation (1) but a conservation constraint can be easily
introduced using Lagrange multipliers, see for example Jacqmin (1999); Yue
_et al._ (2004).

Figure 2: Analysis of interfacial profile for $\phi$. (a) The double-well
potential corresponding to the demixing free energy in Eq. (2). The inset
shows the equilibrium phase field profile from Eq. (3) for different values of
the interface width $\lambda$. The interface becomes sharp as $\lambda$ goes
to $0$. (b) Simulation setup for measuring the interface between two phases
with $\phi=1$ (yellow) and $\phi=0$ (green) in a complex fluid with nematic
order parameter, see section III.1. The phase field profile is measured along
the black solid line and even in this it is very close to the theoretical
$\tanh$ form of $\phi_{\text{eq.}}$.

In order to get a better understanding of the parameters appearing in equation
(2), we search for steady state solutions of equation (1) which leads us to
solve

$\frac{\delta\mathcal{F}_{\text{CH}}}{\delta\phi}=A\phi(1-\phi)(1-2\phi)-\kappa\Delta\phi=0.$

It is easy to verify that

$\phi_{\text{eq.}}(x)=\frac{1}{2}\tanh\left(\frac{x_{0}-x}{\sqrt{4\kappa/A}}\right)+\frac{1}{2},$
(3)

is a solution of the above equation with boundary conditions
$\phi_{\text{eq.}}(-\infty)=1$ and $\phi_{\text{eq.}}(+\infty)=0$. This
corresponds to a one-dimensional interface along the $x$-axis located at
$x_{0}$ with interface width $\lambda=\sqrt{2\kappa/A}$, (see inset of Fig.
2(a)). Note that this solution is the absolute minimum of the free energy for
these boundary conditions, but that infinitely many unstable minima exist
where the phase field oscillates multiple times between $\phi=0$ and $\phi=1$
Mauri  _et al._ (1996).

Inserting $\phi_{\text{eq.}}$ into the free energy density, Eq. (2), and
integrating perpendicularly to the interface allows us to directly obtain the
surface tension between the two phases $\phi=1$ and $0$ as

$\displaystyle\gamma$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}x\,\left\\{\frac{A}{2}\>\phi_{\text{eq.}}^{2}\left(1-\phi_{\text{eq.}}\right)^{2}+\frac{\kappa}{2}(\nabla\phi_{\text{eq.}})^{2}\right\\}$
$\displaystyle=\sqrt{A\kappa}/6.$

A similar analysis can be carried out in two dimension for the case of a
radially symmetrical droplet even though no analytical solution for the
profile exists in this case. Considering a circular droplet whose radius $R$
is much larger that the interface width $\lambda$, it is easy to see that the
Laplacian in radial coordinates can be approximated by
$\Delta=\partial_{r}^{2}+1/r\partial_{r}\approx\partial_{r}^{2}+\mathcal{O}(\lambda/R)$,
in which regime the hyperbolic tangent profile from above is recovered. When
the radius of the droplet is somewhat similar to the interface width, this
term can not be neglected and the interface profile will only be approximated
by $\phi_{\text{eq.}}$. Also note that in the non-stationary case
$\partial_{t}\phi\neq 0$ the interface will also deviate from its equilibrium
profile, but these deviations are found to be small in practice Yue  _et al._
(2004).

Using a free energy-based description makes phase field models of interfaces
easily adaptable for modeling different types of complex microstructured
fluids, such as liquidc crystals or viscoelastic fluids, where the mixing free
energy for the phase field is simply complemented by the free energy of the
microstructured fluid Yue  _et al._ (2004). Another advantage of the phase
field models compared to the other approaches for modeling fluid interfaces is
that it can be proven that the energy is always conserved even in the case of
describing complex fluids such as liquid crystals Lin and Liu (1995).

###  II.2 Coupling to hydrodynamics

One attractive feature of phase models is that their governing equation for
the binary order parameter $\phi$ can be easily coupled to hydrodynamic
equations of fluid flow and other order parameters such as orientation field,
polarity field, and even viscoelasticity. In this part we introduce a general
framework through which such coupling is achieved, starting with a generic
coupling to hydrodynamics and discussing specific examples from active
nematics (where coupling to orientation field is introduced) and active matter
in viscoelastic surrounding (where coupling to viscoelasticity is described).
We further discuss how the generic framework for coupling to other order
parameters can be adapted for more complex setups such as active shells within
fluid backgrounds.

In order to couple phase field model to fluid flow, the governing equation for
the binary order parameter $\phi$ is simply supplemented with an advective
flux of $\phi$, $\vec{j}=\vec{u}\phi$, where $\vec{u}$ is the fluid velocity.
This gives:

$\partial_{t}\phi+\vec{\nabla}\cdot\vec{j}=-\frac{\delta{\mathcal{F}}}{\delta\phi}.$
(4)

As such, coupling to flow is easily achieved without introducing any
additional parameters in the phase field equation. In its most general form,
the fluid flow $\vec{u}$ is in turn determined from the Navier-Stokes
equation:

$\rho\left(\partial_{t}\vec{u}+\vec{u}\cdot\vec{\nabla}\vec{u}\right)=\vec{\nabla}\cdot\bm{\Pi},$
(5)

where $\rho$ is the fluid density and $\bm{\Pi}$ is the stress tensor. The
stress tensor, in general, comprises fluid pressure $p$ and viscous stresses
$\eta\bm{E}$, with $\eta$ the dynamic viscosity and
$\bm{E}=(\vec{\nabla}u)^{\text{Symm}}$ the strain rate tensor that
characterises the symmetric part of the velocity gradient tensor.

The presence of an interface between two phases gives rise to additional
stresses in the flow field. To capture these additional stresses, the back-
coupling from the binary order parameter $\phi$ to the fluid flow equation is
introduced through capillary stresses Sulaiman  _et al._ (2006); Blow  _et
al._ (2014); Cates and Tjhung (2018); Kempf  _et al._ (2019):

$\bm{\Pi}^{\text{capillary}}=\left(\mathcal{F}-\mu\phi\right)\mathds{1}-\vec{\nabla}\phi\frac{\partial\mathcal{F}}{\partial\vec{\nabla}\phi},$
(6)

where $\mu=-\delta\mathcal{F}/\delta\phi$ is the chemical potential and
$\mathds{1}$ is the identity matrix. As such, the first term describes
contributions to the fluid pressure (isotropic contributions to stress tensor)
due to particle exchange between two phases through chemical potential. The
second term on the right-hand-side of the Eq. (6) characterizes anisotropic
contributions to the fluid stress due to the presence of surface tension
between the binary phases. Calculating the divergence of the capillary stress,
it can be shown that the corresponding force field simplifies to Sulaiman  _et
al._ (2006)

$\vec{F}^{\text{capillary}}=\vec{\nabla}\cdot\bm{\Pi}^{\text{capillary}}=\phi\vec{\nabla}\mu.$
(7)

Therefore, from Eqs. (1)-(7), the coupled set of equations for a phase field
model in the presence of hydrodynamic fluid flow can be written as:

$\displaystyle\partial_{t}\phi+\vec{\nabla}\cdot\left(\vec{u}\phi\right)$
$\displaystyle=$ $\displaystyle\mu,$
$\displaystyle\rho\left(\partial_{t}\vec{u}+\vec{u}\cdot\vec{\nabla}\vec{u}\right)$
$\displaystyle=$ $\displaystyle\vec{\nabla}\cdot\bm{\Pi}+\phi\vec{\nabla}\mu,$
(8)

where, as before, the stress tensor $\bm{\Pi}$ includes pressure and viscous
contributions.

Eq. (II.2) can be solved numerically using hybrid schemes combining a finite-
difference integration method LeVeque (2007) for the scalar binary order
parameter with the desired fluid solver, e.g. lattice Boltzmann method
Marenduzzo  _et al._ (2007); Desplat  _et al._ (2001) (see Carenza  _et al._
(2019b) for a recent review of lattice Boltzmann method for simulating active
fluids).

The framework discussed above represents a general formulation for the
coupling between the binary order parameter and hydrodynamic fluid flow and
can be easily adapted to problems where the fluid flow is generated by active
particles. For active fluids, additional active stresses can easily be added
to the stress tensor Cates and Tjhung (2018); Jülicher  _et al._ (2018);
Doostmohammadi  _et al._ (2018) and their exact form depends on the type of
the active system under consideration. The important point, however, is that
the effect on the evolution of the binary order parameter can be described by
only considering the active flux ($\vec{j}=\vec{u}\phi$) that is solely
determined from the resulting velocity field. In what follows, we discuss
specific examples from active nematic models to make this more clear.

##  III Continuum biphasic active systems

Figure 3: Modeling active biphasic systems composed of active nematics and
isotropic phases. (a-d) Temporal evolution of the instability of an active
nematic band (white region) within an isotropic fluid background (blue
region). The black solid lines represent the director field of the active
nematic. An initially flat interface (a) is unstable to bend deformation of
the active nematics (b), which grow in amplitude (c) and generate a highly
deformed interface (d). Figure adapted form Blow  _et al._ (2014). (e-f) A
drop of active nematic elongates (e), forms a narrow neck (f), and eventually
divides (g). The black solid lines represent the director field of the active
nematic and white vectors indicate the flow streamlines. The active nematic
region is marked as red, while the isotropic fluid region is shown in blue.
Figure adapted form Giomi and DeSimone (2014). (h) Growing active nematic
(yellow region) invading a surrounding isotropic fluid (green region). I-III
show different patterns of invasion for increasing activities. Figure adapted
form Kempf  _et al._ (2019). (i) Active nematic fluid (yellow region) in a
background passive liquid crystal develops instabilities that do not grow
beyond a certain length set by the properties of the background liquid
crystal. Red solid lines represent the director field of the background liquid
crystal that are initially patterned to form a splay region in the middle.
Figure adapted form Turiv  _et al._ (2020).

###  III.1 Active nematic interfaces

Active nematics represent a class of active materials that are characterized
by having a local orientational order parameter Cates and Tjhung (2018);
Jülicher  _et al._ (2018); Doostmohammadi  _et al._ (2018). This comprises
solutions of elongated particles such as bacteria Volfson  _et al._ (2008);
Meacock  _et al._ (2020), microtubule filaments put in motion by kinesin motor
proteins Sanchez  _et al._ (2012); Guillamat  _et al._ (2016) or actin
filaments and myosin mixtures Kumar  _et al._ (2018), but also spindle shaped
cells such as mouse embryonic fibroblasts Duclos  _et al._ (2018) and neural
progenitor cells Kawaguchi  _et al._ (2017). More surprisingly, deformable
cells such as Madine-Derby Canine Kidney cells (MDCK) Saw  _et al._ (2017),
human bronchial cells (HBC) Blanch-Mercader  _et al._ (2018), and human
fibrosarcoma Yashunsky  _et al._ (2020) can also be described by active
nematic theories, where the orientational order is defined as the direction of
cell elongation. In nematic systems, the orientational order is characterized
by a second-rank nematic tensor $\bm{Q}$. The field associated with the
orientational order of elongated particles must be represented by a tensor
since it must show head-tail symmetry, i.e. be apolar. For this reason, it can
not be defined using a directional unit vector $\hat{n}$ only and instead
higher moments must be used to define the nematic tensor
${\bm{Q}}=q\frac{d}{d-1}(\hat{n}^{\mathsf{T}}\hat{n}-\mathds{1}/d)$, where $q$
is the magnitude of the nematic order that corresponds to the largest
eigenvalue of $\bm{Q}$ (see Doostmohammadi  _et al._ (2018) for a recent
review of active nematics).

It is well-established that for an active nematic system, the prominent
contribution to active stresses is proportional to the nematic tensor
${\bm{Q}}$ such that the active stress takes the form $\zeta{\bm{Q}}$
Voituriez  _et al._ (2005); Ramaswamy (2010); Marchetti  _et al._ (2013b);
Doostmohammadi  _et al._ (2018), where $\zeta$ determines the strength of
activity. Therefore, in order to account for the effect of activity, an active
stress proportional to the nematic tensor ${\bm{Q}}$ is added to the stress
tensor is Eq. (II.2): $\bm{\Pi}^{\text{active}}=\phi\zeta{\bm{Q}}$. Since the
divergence of the stress determines the force acting on the fluid, this means
that any gradients in the nematic tensor ${\bm{Q}}$ result in additional
forces and actively generate fluid flow Ramaswamy (2010). Moreover, the
proportionality of the stress to the binary order parameter $\phi$ ensures
that the active stress is only applied in the phase with $\phi=1$ and is zero
in the other phase where $\phi=0$. Such a form for the active stress can be
used to model systems with a binary mixture of active and passive nematics
(where both phases have orientational order, but only one of them is active),
such as mixtures of live and dead bacteria Patteson  _et al._ (2018) or a
mixture of motor proteins-filaments in which only a part of the mixture is
activated through light-activated motor proteins Zhang  _et al._ (2019).

In several experimental conditions, however, an active phase of particles with
orientational order, forms an interface with an otherwise isotropic fluid.
This can be easily accommodated within the framework described here by
ensuring that the passive phase $\phi=0$ has no orientational order, i.e.,
${\bm{Q}}\rightarrow{\bm{0}}$ when $\phi\rightarrow 0$. To this end, the
Landau-De Gennes free energy for the nematic
$\mathcal{F}_{\text{L-G}}=\frac{1}{2}{\bm{Q}}^{2}+\frac{1}{4}{\bm{Q}}^{4}$ De
Gennes and Prost (1993); Beris and Edwards (1994), is coupled to the binary
order parameter $\phi$ as follows:

$\mathcal{F}_{\text{L-G}}=\frac{1}{2}\phi{\bm{Q}}^{2}+\frac{1}{4}{\bm{Q}}^{4}.$
(9)

This free energy ensures that in the active nematic phase $\phi=1$, the
orientational order is retained, while in the passive phase $\phi=0$, the
magnitude of orientational order that minimizes the free energy — i.e.
${\bm{Q}}$ for which $\delta\mathcal{F}/\delta{\bm{Q}}=0$ — diminishes to
zero, resulting in an isotropic fluid Blow  _et al._ (2014); Kempf  _et al._
(2019). In Figure 2b, we compare analytical results for the one-dimensional
$\phi$ profile to the results from numerical integration of Eq. (1).

An interesting feature of active nematic interfaces compared to their passive
counterpart, is the emergence of the phenomenon of active anchoring Blow  _et
al._ (2014): in the absence of any anchoring energies, the presence of
activity alone results in a certain alignment of active particles at the
interface. This effect has been shown for dividing bacteria, where cell
division provides a source of active force generation, that align parallel to
the interface of their colony Doostmohammadi  _et al._ (2016) and can also be
seen in Fig. 1(a). Considering the active force

$\displaystyle F_{\alpha}^{\text{active}}=$
$\displaystyle-\zeta\partial_{\beta}\left(\phi Q_{\alpha\beta}\right)$
$\displaystyle=$ $\displaystyle-\zeta\left(\partial_{\beta}\phi
S\right)\left(2n_{\alpha}n_{\beta}-\delta_{\alpha\beta}\right)$
$\displaystyle-2\;\zeta\phi
S\left(n_{\alpha}\left(\partial_{\beta}n_{\beta}\right)+\left(\partial_{\beta}n_{\alpha}\right)n_{\beta}\right),$
(10)

the active anchoring can be understood by decomposing the active force
parallel and perpendicular to the interface between active and passive phases
Blow  _et al._ (2014):

$\displaystyle f_{\perp}$
$\displaystyle=m_{\alpha}f_{\alpha}=\zeta\left\lvert\nabla\left(\phi
S\right)\right\rvert\left(2\left(m_{\alpha}n_{\alpha}\right)^{2}-1\right),$
(11) $\displaystyle f_{\parallel}$
$\displaystyle=l_{\alpha}f_{\alpha}=2\;\zeta\left\lvert\nabla\left(\phi
S\right)\right\rvert\left(l_{\alpha}n_{\alpha}\right)\left(m_{\beta}n_{\beta}\right),$
(12)

where
$\vec{m}\equiv-\vec{\nabla}S/\left\lvert\vec{\nabla}S\right\rvert\equiv-\vec{\nabla}\left(\phi
S\right)/\left\lvert\vec{\nabla}\left(\phi S\right)\right\rvert$ and $\vec{l}$
define the outward surface normal and the surface tangent directions,
respectively.

Depending on the sign of the activity parameter $\zeta$, the normal force on
the interface $f_{\perp}$, points outwards ($\zeta<0$) or inwards ($\zeta>0$)
for $\vec{n}\perp\vec{m}$ and hence induces active forces pushing the
interface outwards or inwards, respectively. Therefore, for $\zeta>0$ activity
favors alignment of active particles parallel to the interface, while
$\zeta<0$ favors perpendicular alignment. The case of $\zeta>0$ describes
extensile active systems such as bacteria, where particles pull the
surrounding fluid from their side and push it along their tail and head. The
$\zeta<0$ case, on the other hand, corresponds to contractile active systems
such as actomyosin networks, where active constituents contract along their
direction of elongation.

The method presented here has been extensively used to model fundamental
aspects of active interfaces, including the study of the instabilities of
active nematic interfaces Blow  _et al._ (2014, 2017) (Fig.3(a-d)), describing
dynamics of elongation, motility, and division of active nematic droplets as a
model of eukaryotic cells Giomi and DeSimone (2014); Carenza  _et al._ (2019a)
(Fig.3(e-g)), and wetting dynamics of active nematics on solid surfaces Coelho
_et al._ (2020); Neta  _et al._ (2020). Furthermore, the biphasic modeling of
active nematic interaction with isotropic fluid has been used to describe
morphologies of growing bacterial colonies Doostmohammadi  _et al._ (2016);
Dell’Arciprete  _et al._ (2018), active matter invasion into new territories
Kempf  _et al._ (2019) (Fig.3(h)), and bacterial streams within pre-patterned
medium in living liquid crystals Turiv  _et al._ (2020) (Fig. 3(i)).

###  III.2 Active matter in a viscoelastic medium

Figure 4: Modeling active biphasic system composed of active nematics and
passive polymeric phases. (a) A schematic of various interpretations of the
two-phase active-viscoelastic model. (b) Typical flow field of a drop of
active nematics elongating during the initial steps of cell division. Velocity
vectors are shown as red arrows, and polymer deformations are shown as blue
line segments. Inset: Zoom of the lower left part of the cell, showing
polymers aligning with the stretching direction. Figure adapted from Emmanuel
_et al._ (2020).

The exact same procedure used in the previous section can be followed to
introduce couplings to other order parameters. One important example is active
matter in a viscoelastic environment which has applications to the modelling
of several biophysical systems of interest. For example, cells in
extracellular matrices are characterized by the cross talk between an active
phase and a medium which has viscoelastic properties Chaudhuri  _et al._
(2020) (Fig.4(a)). Following the framework introduced in the previous section,
we extend our formalism to the case of an active nematic in contact with a
fluid phase that, rather than being isotropic, has now viscoelastic
properties.

To characterize a viscoelastic fluid, it is common to introduce a conformation
tensor ${\bm{C}}$ as the corresponding order parameter, which describes the
orientation and elongation of polymers that are present in the fluid Hemingway
_et al._ (2015); Emmanuel  _et al._ (2020). The conformation tensor $\bm{C}$
is defined such that: (i) it satisfies ${\bm{C}}=\mathds{1}$ at equilibrium,
(ii) its trace $\mathrm{tr}[{\bm{C}}]$ characterizes the square of the polymer
elongation, and (iii) the eigenvector corresponding to its largest eigenvalue
represents the local polymer orientation. Therefore, as before, a coupling
between the binary order parameter and the conformation tensor can be
introduced such that within the active phase ($\phi=1$) there is an
orientational order characterized by the nematic tensor ${\bm{Q}}$ but no
viscoelasticity ${\bm{C}}={\bm{0}}$, while in the passive phase ($\phi=0$)
there are no active nematic particles ${\bm{Q}}={\bm{0}}$ but polymers are
present ${\bm{C}}\neq{\bm{0}}$. To ensure this, the free energy of polymers is
coupled to the binary order parameter $\phi$ as follows:

$\mathcal{F}_{\text{polymer}}=\frac{\nu}{\tau}\frac{1}{2}\left(1-\phi\right)\left(\mathrm{tr}[{\bm{C}}-\mathds{1}]-\log\left(\det[{\bm{C}}]\right)\right),$
(13)

where $\nu$ characterizes the polymer viscosity and $\tau$ is the polymer
relaxation time. Note that the factor $1-\phi$ in front of the polymer free
energy ensures that the order parameter is retained within the passive phase
$\phi=0$, while it diminishes otherwise when $\phi=1$. Finally, since polymers
are within the passive phase, no extra addition to active stresses is
required. However, passive stresses need to be accounted for in Eq. (II.2) as
follows:

${\bm{\Pi}}^{\text{polymer}}=\frac{\nu}{\tau}(1-\phi)({\bm{C}}-\mathds{1}).$
(14)

It is noteworthy that the formulation described above uses an Oldroyd-B model
of polymer relaxation since it is one the simplest and most widely-used
polymer models Emmanuel  _et al._ (2020). The framework can be trivially
extended to other constitutive models of viscoelastic fluids in two and three
dimensions.

The important feature of the present framework that couples viscoelasticity,
fluid flow, and activity through the phase field formulation is its
versatility, which means that it can be applied to various biophysical
examples. For instance, similar formulation has been used to study a model of
cell division and motility within polymer gels with variable viscoelastic
properties (Fig.4(b)), and to study the hampering effect of polymers on
generic instabilities of active matter in a viscoelastic environment Emmanuel
_et al._ (2020) (Fig. 4).

###  III.3 Self-deforming active surfaces in a fluid background

So far, we have described applications of the phase field method to binary
mixtures, where $\phi=0$ is used to characterize one phase and $\phi=1$ to
describe the other, with $\phi=0.5$ marking the interface between the two
phases. It is, however, possible to extend this framework to model scenarios
in which thin shells of one phase are embedded within another phase. An
important example is the actomyosin cortex of eukaryotic cells that plays a
vital role in shaping cell deformation and various cell functions Naganathan
_et al._ (2014). Another interesting example is the active deformable shell
formed by stabilizing microtubule-motor protein mixtures at the oil-water
interface of a droplet Keber  _et al._ (2014); Zhang  _et al._ (2016);
Guillamat  _et al._ (2018). The resulting active shell continuously exerts
active stresses on the surrounding medium, generating dynamic patterns of
motion and exotic morphologies that can be considered as examples of self-
deforming and self-shaping materials (see Fig. 5(a)).

The distinguishing feature of these systems is that the shell itself is active
and as such activity is generated at a thin interface between the inside and
the outside of the shell. The phase field model of active interfaces described
in this section can be easily adapted to model such as active shells by
describing the coupling between the binary order parameter $\phi$ and the
desired order parameter (in the case of active nematic this is the nematic
tensor ${\bm{Q}}$ which is the source of active stresses) such that the active
phase is only stabilized at the interface, i.e., where $\phi\approx 0.5$. For
instance, starting from the active nematic system described in Sec. III.1, it
is only required to modify the free energy as follows:

$\mathcal{F}_{\text{L-G}}=\frac{1}{2}\left(1-2\lvert\phi-\frac{1}{2}\rvert\right){\bm{Q}}^{2}+\frac{1}{4}{\bm{Q}}^{4},$
(15)

where the term $\left(1-2\lvert\phi-\frac{1}{2}\rvert\right)$ ensures that
${\bm{Q}}\to{\bm{0}}$ as both $\phi\to 0$ and $\phi\to 1$, while the
orientational order is retained at the interface ($\phi\approx 0.5$). Again,
this provides only a simple example of how such coupling can be achieved. The
framework can be trivially extended to consider further complexities of the
interface, for example, by including Helfrich type free energies that describe
bending rigidity of the shell Metselaar  _et al._ (2019).

This method has been recently applied to model self-deforming active nematic
shells in fluid backgrounds, reproducing topological defect dynamics on the
shell surface observed in the experiments on microtubule-kinesin motor
mixtures Keber  _et al._ (2014); Zhang  _et al._ (2016); Metselaar  _et al._
(2019), and predicting exotic morphologies that can be obtained by tuning the
activity of the shell. More importantly, such modeling provides a tool to
probe the nature of cross talk between activity and self-induced surface
curvature. For example, it is predicted that protrusions are initiated at the
locations of $+1/2$ defects in the orientation field, creating finger-like
structures with a $+1$ topological defect at their tip, leaving behind $-1/2$
topological defects in regions with negative surface curvature (Fig. 5(b)).
Similar dynamics and cross talk between activity and morphology has been
recently reported in the development of regenerating animal hydra where actin
filaments form an active shell that is capable of self-deformation Maroudas-
Sacks  _et al._ (2020) (Fig. 5(c)).

Figure 5: Modeling three-dimensional active self-deforming surfaces. (a-c)
Temporal evolution of an active shell made of active nematics in an isotropic
fluid background. Active shell develops protrusions at the position of $+1/2$
defects. The protrusions eventuality develop $+1$ defects at their tip leaving
behind $-1/2$ defects. The solid lines illustrate the nematic director and are
colored by the magnitude of the nematic order. The shell is colored by the
curvature of the surface. Figure adapted from Metselaar  _et al._ (2019).
(d-e) Developments of protrusions and topological defects during the
morphogenesis of regenerating animal hydra. Figure adapted from Maroudas-Sacks
_et al._ (2020). (f) Temporal evolution of protrusions in a shell of active
nematics build from microtubule-motor protein mixtures. Figure adapted from
Keber  _et al._ (2014).

The examples described above provide an introduction to the applications of
phase field model for modeling active systems. We focused particularly on
active nematic interfaces, where dynamics of phase field is coupled to a
tensor order parameter and which describes orientational order of active
particles. Importantly, phase field models have been applied to various
realizations of active systems including scalar phase field models describing
the phenomenon of motility-induced phase separation Wittkowski  _et al._
(2014) (see Cates and Tjhung (2018) for a recent review), as well as polar
active matter, where phase field dynamics is coupled to a vector order
parameter describing the polarity of active particles Negro  _et al._ (2019);
Carenza  _et al._ (2019b). The latter has been particularly successful in
modeling cellular motility, where polar active matter is used to describes the
dynamics of actomyosin networks inside the cell, governing cell deformation
and its morphology Ziebert  _et al._ (2011); Tjhung  _et al._ (2012); Marth
and Voigt (2016); Carenza  _et al._ (2019a)

##  IV Phase field modeling of cell monolayers

In the previous sections we discussed applications of phase field models to
continuum representations of active interfaces. Interestingly, the same
framework can be adapted to agent-based modeling of active systems where each
agent is represented by an individual phase. An important example is the phase
field modeling of cell monolayers which allows to model the deformations of
individual cells and describe their physical properties, as well as
introducing cell-cell interactions easily through the interfaces of the
individual phase fields. In this section, we discuss how such models are
constructed and applied to the understanding of the dynamics of cellular
layers.

###  IV.1 Cells as active deformable droplets

Figure 6: Modeling cell monolayer composed of multiple active phase fields
representing each cell. (a) Schematic representation of the phase fields of
two cells (green and blue) as well as their overlap (red) seen from the side.
(b) A phase field representation of a cell monolayer as a collection of active
deformable drops seen from the top. The square around the cell represent the
sub-domain on which the equations of motion are solved, see section IV.3. (c)
A confluent monolayer in confinement where the overlap between cells and their
centres of mass are shown instead of the cell contours.

In silico models of cell motility have an important role to play in
unravelling the interplay between single cell properties and their collective
dynamics Bi  _et al._ (2015); Alert and Trepat (2020). There exist a breadth
of numerical approaches to tackle this problem, such as cellular Potts models
Graner and Glazier (1992), vertex models Farhadifar  _et al._ (2007); Bi  _et
al._ (2016); Sussman  _et al._ (2018), continuum models Prost  _et al._
(2015); Jülicher  _et al._ (2018), and phase field models Löber  _et al._
(2015), see the extensive reviews Alert and Trepat (2020); Camley and Rappel
(2017); Shaebani  _et al._ (2020). The phase field approach is particularly
interesting as it allows to describe the collective behaviour of cells while
modelling their physical properties at the individual level as well as their
interactions explicitly. It has been widely applied to problems involving
single cells Löber  _et al._ (2014); Ziebert  _et al._ (2011); Ziebert and
Aranson (2013); Tjhung  _et al._ (2012), few migrating cells Najem and Grant
(2016), as well as systems of colliding binary cells Camley  _et al._ (2014);
Löber  _et al._ (2015). More recent works Marth and Voigt (2016); Palmieri
_et al._ (2015); Peyret  _et al._ (2019); Wenzel  _et al._ (2019); Mueller
_et al._ (2019); Loewe  _et al._ (2019) have concentrated on the collective
dynamics of a large number of cells and have uncovered interesting phenomena
such as oscillatory patterns of cells under confinement Peyret  _et al._
(2019), as well as solid-liquid and flocking transitions Giavazzi  _et al._
(2018). However, since this approach is computationally quite intensive, it is
necessary to reduce the description of cells to their bare minimum by
abstracting away their internal machinery while retaining their relevant
physical properties. Here we concentrate on a simple model that represents
cells on a substrate in two-dimensions and assumes simple physical properties
for the cells.

A fruitful approach first pioneered by Palmieri  _et al._ (2015), is to
separate the description of the physical properties of the cells from their
dynamics. In this model, each cell is described as a deformable active
particle whose shape is defined using an individual phase field (Fig. 6(a-b)),
while its velocity is described separately and is given by a force-balance
equation. A cellular monolayer consisting of $N$ cells can then be described
by the equations of motion

$\partial_{t}\phi_{i}+\vec{v}_{i}\cdot\vec{\nabla}\phi_{i}=-\frac{\delta{\mathcal{F}}}{\delta\phi_{i}},\qquad
i=1,\ldots,N,$ (16)

where ${\mathcal{F}}$ is a free energy describing the physical properties of
the cells, and $\phi_{i}$ and $\vec{v}_{i}$ are the phase field and the
center-of-mass velocity of cell $i$, respectively. The phase fields are
defined such that $\phi=0$ denotes the exterior of the cell $i$ while $\phi=1$
denotes its interior. This mirrors Eq. (1) with the difference that there is a
different phase field for each cell and that each individual cell and is
driven by a simple advection term. The velocity $\vec{v}_{i}$ for each cell
can be obtained from a force balance equation: since Reynolds numbers are
typically of the order of $\sim 10^{-4}$ for cell monolayers Kim  _et al._
(2013), we can safely assume overdamped dynamics and write

$\xi\vec{v}_{i}=\vec{F}^{\text{tot.}}_{i}=\vec{F}^{\text{pol.}}_{i}+\vec{F}^{\text{inter.}}_{i},$
(17)

where $\xi$ is a friction coefficient and $\vec{F}^{\text{tot.}}_{i}$ is the
total force exerted on cell $i$, which will be defined below. Note that
$\vec{v}_{i}$ is a single vector for each cell and does not depend on
$\vec{x}$. The definition of the total force $\vec{F}^{\text{tot.}}_{i}$
depends on the model and can include contributions from different sources such
as passive interface forces, interactions with other cells or walls, or active
forces such as polar or nematic driving. Note that we have assumed here for
simplicity that the free energy is the same for all cells but this does not
need to be the case in general and it is straightforward to extend the above
equation to cases where different cells have different physical properties,
see for example Balasubramaniam  _et al._ (2020). Equations (16) and (17) are
our master equations and define a general framework where only the free-energy
$\mathcal{F}$ and the total force $\vec{F}^{\text{tot.}}_{i}$ on each cell
need to be specified.

###  IV.2 Minimal model of cell monolayers

Let us now define a minimal model that allows the successful description of
cell monolayers following Palmieri  _et al._ (2015); Mueller  _et al._ (2019)
and write the free energy as a sum of three distinct contributions as
$\mathcal{F}=\mathcal{F}_{\text{CH}}+\mathcal{F}_{\text{area}}+\mathcal{F}_{\text{rep.}}$.
The first term is the Cahn-Hilliard free energy (2), which we rewrite as

$\mathcal{F}_{\text{CH}}=\sum_{i}\frac{\gamma}{\lambda}\int\mathrm{d}\vec{x}\left\\{4\phi_{i}^{2}(1-\phi_{i})^{2}+\lambda^{2}(\nabla\phi_{i})^{2}\right\\},$
(18)

using the definitions of the interface width $\lambda=2\sqrt{\kappa/A}$ and
surface tension $\gamma=\sqrt{A\kappa}/6$ from sec. II (up to an overall
factor of $3$). As shown in sec. II, this term creates and stabilises an
interface between the interior ($\phi_{i}=1$) and the exterior ($\phi_{i}=0$)
of each cell and results in a equilibrium shape for the interface that follows
an approximate hyperbolic tangent profile of size $\lambda$. The second term
is a soft area constraint defined by

$\mathcal{F}_{\text{area}}=\sum_{i}\mu\Big{(}1-\frac{1}{\pi
R^{2}}\int\mathrm{d}\vec{x}\,\phi_{i}^{2}\Big{)}^{2},$ (19)

which is a square potential of strength $\mu$ ensuring the cells areas
$A_{i}=\int\mathrm{d}x\,\phi_{i}^{2}$ are close to $\pi R^{2}$ 222The area is
proportional to the _square_ in order to ensure that it is always positive
even when the phase field is slightly negative which can happen during
simulations.. Note that the phase field $\phi_{i}$ is not conserved in Eq.
(16) and that even though cells are mostly incompressible in three dimensions,
the area that individual cells occupy within a monolayer can dramatically
change as they are squeezed by their neighbours and expand in the direction
perpendicular to the substrate Deforet  _et al._ (2014). The final term
discourages overlap between cells and is simply given by

$\mathcal{F}_{\text{rep.}}=\sum_{i}\sum_{j\neq
i}\frac{\kappa}{\lambda}\int\mathrm{d}\vec{x}\,\phi_{i}^{2}\phi_{j}^{2}\,,$
(20)

where $\kappa$ is the strength of the repulsion. Note that this term
implements a coupling between the cells and is proportional to the cell
overlaps, see an illustration on Fig. 6(c). The physical properties of the
cells are parameterized by $\lambda$, $\gamma$, $\mu$, and $\kappa$ which set
the interface width, the surface tension, the strength of the elastic
restoring force for the area, and the repulsive force between cells,
respectively. Normalisation is chosen such that properties of the cells are
roughly preserved when the interface width $\lambda$ is rescaled, see Mueller
_et al._ (2019).

Cellular deformations such as changes of shape or area should lead to forces
at the cell boundaries and contribute to the overall force balance. Such
interface forces can be constructed in a thermodynamically consistent way
using the free energy as

$\vec{F}_{i}^{\text{inter.}}=\int\mathrm{d}\vec{x}\sum_{j}\left(\frac{\delta\mathcal{F}_{\text{rep.}}}{\delta\phi_{j}}-\frac{\delta\mathcal{F}_{\text{CH}}}{\delta\phi_{j}}-\frac{\delta\mathcal{F}_{\text{area}}}{\delta\phi_{j}}\right)\vec{\nabla}\phi_{i},$
(21)

see Palmieri  _et al._ (2015); Mueller  _et al._ (2019) for more details.
Since $\vec{\nabla}\phi_{i}$ is only non-zero at the cell boundary and
pointing towards the cells center, this expression can be interpreted as the
integral over the interface of the cell $i$ of the total force density
generated by changes of the free-energy. Contributions from interactions with
other cells come with a positive sign (repulsion) while self contributions
come with a negative sign (restoring forces). With this and in the absence of
any active contribution, a confluent cell monolayer will relax to a minimum of
the total free energy and favour hexagonal lattice arrangement (Fig. 6(c)).

The great strength of this formulation is its versatility: it can be easily
extended to more complex models that can include forces such as adhesion
between cells Zhang  _et al._ (2020b) or viscous friction Peyret  _et al._
(2019). Moreover, as shown in Peyret  _et al._ (2019), it is also easy to
model non-trivial boundary conditions (such as walls) by introducing a static
phase field $\phi_{\text{walls}}$ that interacts with the cells with a simple
repulsion term.

###  IV.3 Efficient simulation algorithm

The coupled system of differential equations defined by (16) and (17) can be
solved using standard finite difference schemes LeVeque (2007). The main
difficulty is the large number of phase fields $\phi_{i}$ that are required
when simulating a monolayer with many cells. In fact, it is easy to see that a
naive simulation will have at least quadratic computational and memory
requirements with respect to the number of cells due to each phase field
extending over the whole domain. This can be mitigated by realising that the
individual phase fields are non-zero only in a well-defined region around the
center of the corresponding cell, which allows to simulate each phase field
only in a sub-domain centered around its center-of-mass Nonomura (2012);
Mueller  _et al._ (2019), see Fig. 6(b).

A common choice is to select a square of lattice points that tracks the center
of mass of each cell and use Dirichlet boundary conditions at the boundary of
the sub-domain Loewe  _et al._ (2019). It is important to note, however, that
because the centre of mass of the cell can lie between lattice points, the
location of the sub-domain changes abruptly every time the cell center moves
by a lattice length along any axis. This can lead to numerical artifacts when
the size of sub-domain is not large enough compared to the radius of the cell.
A more refined approach is to use periodic boundary conditions instead Mueller
_et al._ (2019), which has the advantage of being completely insensitive to
the discrete hopping of the sub-domain and allows smaller sub-domain sizes
compared to Dirichlet boundary conditions.

With such techniques, both the computational and memory costs of the algorithm
scale linearly with respect to the number of cells when simulating a
multicellular monolayer. Note that even though cells are only simulated around
their center of mass, the computation of the total force acting on a cell
interface will still require the computation of some quantities over the whole
domain, such as the sum of all phase fields $\sum_{i}\phi_{i}$ and their
squares $\sum_{i}\phi_{i}^{2}$, but the number of such quantities is fixed and
does not scale with the number of cells. An implementation by the authors of
the model presented in section IV.2 using these techniques can be found online
at https://github.com/rhomu/celadro.

###  IV.4 Active polar driving

Figure 7: Active forces and instabilities in cell monolayer. (a) Schematic
representation of different sources of active force. Figure adapted from Zhang
_et al._ (2020b). (b) Coarse-grained nematic field of the cells showing the
emergence of multicellular scale bend instability for extensile active forces
and splay instability for contractile active forces acting at the individual
cell level. The cells are driven out of equilibrium by an internal nematic
degree of freedom.

The driving forces behind the movements of single cells on a flat substrate
are well understood. Directional actin filaments, which continuously
polymerise and depolymerize, allow the cell to create a pushing force against
the substrate via focal adhesion, which are mechanical links between internal
actin bundles and the external surface Sarangi  _et al._ (2017). In doing so,
the cell polarizes its actin filaments and tends to elongate in the direction
of motion, creating protrusions (lamellipodia) Mitchison and Cramer (1996);
Alberts  _et al._ (2015). This suggests a minimal physical picture of cellular
motility where each cell experiences a net active force

$\vec{F}_{i}^{\text{pol.}}=\alpha\vec{p}_{i}$

in the direction of its polarity $\vec{p}_{i}$, where $\alpha$ is a parameter
denoting the strength of the focal adhesion (Fig. 7(a)). This force
contributes to the force-balance equation (17) and propels the phase field
$\phi_{i}$ uniformly in the direction of $\vec{p}_{i}$. A simple non-trivial
implementation can be obtained by introducing for each cell a new internal
degree of freedom $\theta_{i}^{\text{pol.}}$ describing the direction of its
polarity
$\vec{p}_{i}=(\cos\theta_{i}^{\text{pol.}},\sin\theta_{i}^{\text{pol.}})$ and
assume that the polarization angle relaxes diffusively towards an arbitrary
vector $\vec{d}_{i}$ Peyret  _et al._ (2019). This leads to the following
equation of motion

$\partial_{t}\theta_{i}^{\text{pol.}}=-J^{\text{pol.}}\,\lvert\vec{d}_{i}\rvert\,\measuredangle(\vec{p}_{i},\vec{d}_{i})+(2D^{\text{pol.}})^{\frac{1}{2}}\,\eta_{i}\,,$
(22)

where $\measuredangle(\vec{p}_{i},\vec{d}_{i})$ is the angle between
$\vec{p}_{i}$ and $\vec{d}$, $\eta_{i}$ is Gaussian white noise, and the
positive constants $J^{\text{pol.}}$ and $D^{\text{pol.}}$ are the strength of
the alignment torque and of the rotational diffusivity, respectively. This
defines a dynamics similar to an Ornstein-Uhlenbeck process and has been well
studied in the context of active systems Szabó  _et al._ (2006); Giavazzi  _et
al._ (2018).

There are multiple interesting choices for the definition of the aligning
direction $\vec{d}_{i}$ and there is as of today no clear microscopic or
dynamical basis favouring one of them Alert and Trepat (2020). A simple choice
that is consistent with the mechanism of contact inhibition of locomotion
(CIL) Smeets  _et al._ (2016) is to assume that the polarity aligns to the
direction of the total force exerted on a cell’s interface, namely
$\vec{d}_{i}=\vec{F}^{\text{int.}}_{i}$. Such a coupling was shown to describe
accurately patterns of sustained oscillations observed experimentally in
systems of confined MDCK cells Peyret  _et al._ (2019). Other interesting
possibilities are to assume that the polarity aligns either to the total
velocity $\vec{v}_{i}$ of each cell or to its main axis of elongation, and an
extensive evaluation of these two choices is presented in Alert and Trepat
(2020); Zhang  _et al._ (2020b). At high enough alignment strengths, the
former case shows a Vicsek-type phase transition where all cells move in the
same direction, while the latter exhibits unjamming to a liquid-like state
Zhang  _et al._ (2020b). Finally, the purely diffusive dynamics of the system
($J^{\text{pol.}}=0$) is interesting in its own right and it has been shown
that densely packed cells can show “bursts” of motions as they quickly relax
from their deformed shape in this case Palmieri  _et al._ (2015), while
confluent monolayers exhibit a solid-liquid phase transition as
$D^{\text{pol.}}$ is increased Loewe  _et al._ (2019); Zhang  _et al._
(2020b). Further research should investigate potential microscopic bases for
such alignment dynamics of the polarization and put on a firmer footing.

###  IV.5 Active nematic driving

Figure 8: Emergent features in a phase field model of cell monolayer. (a)
Emergence of active turbulence. Shown here are the cells in confluence, their
associated coarse-grained nematic field, and the corresponding velocity field.
The nematic field describes the local orientation of the cells defined as
their elongation direction. The singularities in the orientation field are
marked by green circles for $+1/2$ and blue triangles for $-1/2$ topological
defects. The velocity field is colored by the normalized magnitude of velocity
ranging from $0$ (purple) to $1$ (yellow). (b, c) Average velocity field,
isotropic stress, and pressure around $-1/2$ (b) and $+1/2$ topological
defects (c). Figures adapted from Mueller  _et al._ (2019).

While polar driving has been extensively studied theoretically in the context
of single cells, there is strong evidence that descriptions based on active
nematic liquid crystals offer a compelling framework to understand the
collective behaviour of cell monolayers Saw  _et al._ (2017); Blanch-Mercader
_et al._ (2018); Mueller  _et al._ (2019). Nematic driving can be introduced
following Mueller  _et al._ (2019) by rewriting the interface force (21) as

$\vec{F}^{\text{inter.}}_{i}=\int\mathrm{d}\vec{x}\,\phi_{i}\vec{\nabla}\bm{\Pi}^{\text{tissue}}=-\int\mathrm{d}\vec{x}\,\bm{\Pi}^{\text{tissue}}\vec{\nabla}\phi_{i},$
(23)

where $\bm{\Pi}^{\text{tissue}}$ is a tissue stress tensor that can be written
in the usual fashion as

$\bm{\Pi}^{\text{tissue}}=-p\mathds{1}-\zeta\bm{Q}.$ (24)

The pressure
$p=\delta(\mathcal{F}_{\text{rep.}}-\mathcal{F}_{\text{CH}}-\mathcal{F}_{\text{area}})/\delta\phi_{i}$
can be read directly from equation (21), while
$\bm{Q}=\sum_{i}\phi_{i}\bm{Q}_{i}$ is a newly introduced tissue nematic
tensor written as a weighted sum over the contributions of the individual
cells. Such a definition has two main advantages: (i) it bridges the gap
between local (cell level) and global (tissue level) properties and (ii)
allows for naturally introducing the usual active term $-\zeta\bm{Q}$ found in
continuum theories of active liquid crystals Marchetti  _et al._ (2013b);
Doostmohammadi  _et al._ (2018). It can be interpreted as creating a dipolar
force density distributed along the cells interfaces such that each cell
pushes or pulls its neighbours depending on the direction of their contact
area with respect to the nematic tensor (Fig. 7(a)).

The simplest choice for $\bm{Q}_{i}$ is to introduce an internal nematic
degree of freedom for each cell and define its dynamics similarly to the polar
case. In the continuum theory, most of phenomenology of active liquid crystals
depends on a balance between the restoring elastic forces and flow alignment
of the nematic tensor. Writing
$\bm{Q}_{i}=2\left(\hat{n}_{i}^{\mathsf{T}}\hat{n}_{i}-\mathds{1}/2\right)$
with
$\hat{n}_{i}=(\cos\theta_{i}^{\text{nem.}},\sin\theta_{i}^{\text{nem.}})$, we
can mirror these two components in our model by defining the following
dynamics of the angle $\theta_{i}$:

$\partial_{t}\theta^{\text{nem.}}_{i}=K^{\text{nem.}}\tau_{i}+J^{\text{nem.}}\omega_{i},$
(25)

where the torques are given by

$\displaystyle\tau_{i}=\frac{1}{\lambda}\int\mathrm{d}\vec{x}\;\phi_{i}\,\bm{Q}\wedge\bm{Q}_{i},\quad\omega_{i}=\int\mathrm{d}\vec{x}\;\vec{v}\wedge\nabla\phi_{i}.$
(26)

We have defined the tissue velocity as $\vec{v}=\sum_{i}\phi_{i}\vec{v}_{i}$
and $\bm{A}\wedge\bm{B}=A_{xx}B_{xy}-A_{xy}B_{xx}$ for symmetric traceless
matrices $\bm{A}$ and $\bm{B}$. The first torque $\tau_{i}$ aligns
$\bm{Q}_{i}$ to the tissue nematic tensor $\bm{Q}$ and induces an elastic
restoring force favouring the homogeneous state. The second torque
$\omega_{i}$ rotates $\bm{Q}_{i}$ with the local vorticity computed as the
integral of the neighbouring cells velocity projected on the cell interface.
Toghether with the active term, these torques are able to drive the cells out
of equilibrium and reproduce the bend and splay instabilities observed in
continuum theories of active nematic liquid crystals (Fig. 7(b)). In a
confluent epithelium, one can also check that the total force is approximately
zero and the system does not develop any total velocity under periodic
boundary conditions. In particular this means the transition to collective
movement in is different from the Vicsek-type phase transitions observed with
polar driving Giavazzi  _et al._ (2018).

Another interesting possibility explored in Mueller  _et al._ (2019) is to
avoid introducing supplementary degrees of freedom altogether and define an
active coupling that is directly proportional to the shape deformations of the
individual cells. This is motivated by recent studies indicating that the
local deformation of cells provides a suitable nematic order parameter that
allows the description of the dynamics of epithelial cells using theories of
active liquid crystals Doostmohammadi  _et al._ (2015); Duclos  _et al._
(2017); Saw  _et al._ (2017); Duclos  _et al._ (2018); Blanch-Mercader  _et
al._ (2018). Such a connection is quite surprising because individual
epithelial cells on a substrate are typically not elongated and have a well-
defined direction of movement, suggesting polar rather than nematic driving.
In this setting, the tissue nematic tensor is defined as
$\bm{Q}=\sum_{i}\phi_{i}\bm{S}_{i}$, where $\bm{S}_{i}$ is the deformation
tensor of cell $i$ given by

$\bm{S}_{i}=-\int\mathrm{d}\vec{x}\left((\vec{\nabla}\phi_{i})^{\mathsf{T}}\vec{\nabla}\phi_{i}-\frac{1}{2}\mathrm{tr}\left[(\vec{\nabla}\phi_{i})^{\mathsf{T}}\vec{\nabla}\phi_{i}\right]\right).$
(27)

Equation (27) defines a $2\times 2$ matrix whose eigenvalues and eigenvectors
describe the strength and orientation of the main deformation axes of each
cell, see Mueller  _et al._ (2019); Asipauskas  _et al._ (2003) for details.
This corresponds to a nematic tensor with order parameter and director given
by the largest eigenvalue and its associated eigenvector. With this simple
definition, the system shows an activity-driven transition to non-zero nematic
order and flows for high enough activity strengths $\zeta$. The spontaneous
creation of defects in the coarse grained nematic field as well as observed
patterns of flows and mechanical stresses around topological defects are
accurately predicted, see Fig. 8. This points to a strong connection between
the shape deformation of cells and their active behaviour, irrespective of the
corresponding microscopic mechanism.

##  V Summary

We have presented an introduction to phase field modeling of active systems,
providing examples from both continuum biphasic active materials and agent-
based models of cellular monolayers. A general framework that allows the
introduction of various forms of coupling to the phase field was developed. In
particular, it was shown how hydrodynamic effects can be included in phase
field models and how tensor order parameters such as active particles
orientation field or polymer conformation (describing viscoelastic effects)
can be coupled to the equations of motions of the phase field in order to
construct complex models of active interfaces. Furthermore, we have described
how using multiple phase fields can be used to construct a versatile model of
cellular monolayers, where each cell is described as an active deformable
droplet using an individual phase field. Applications of this approach to
modeling collective cell motion and capturing cell-cell interaction forces
were described.

As we gain more knowledge of complex spatiotemporal features of active
systems, understanding the dynamics of physical interfaces between different
active materials as well as between active matter and its surrounding media
becomes more and more important. Such understanding will be of prime
importance to tackle biophysical problems such as bacterial biofilm
development, cellular invasion, and morphogenesis. Building versatile and
predictive computational models of active interfaces will be, together with
experimental advances, the determining factor in advancing our understanding
of such systems, both by helping to explain the mechanisms behind experimental
observations and, even more importantly, by providing predictive tools for
experimentally inaccessible scenarios.

##  VI Acknowledgement

The notes here are part of the “Initial Training on Numerical Methods of
Active Matter” organised by MSCA-ITN ActiveMatter (Grant: 812780) funded by
the EU H2020 programme. AD acknowledges support from the Novo Nordisk
Foundation (grant no. NNF18SA0035142), Villum Fonden (grant no. 29476), Danish
Council for Independent Research, Natural Sciences (DFF-117155-1001), and
funding from the European Union’s Horizon 2020 research and innovation program
under the Marie Sklodowska-Curie grant agreement no. 847523 (INTERACTIONS).

## References

  * Meacock  _et al._ (2020) O. J. Meacock, A. Doostmohammadi, K. R. Foster, J. M. Yeomans, and W. M. Durham, Nature Physics , 1 (2020).
  * Doostmohammadi  _et al._ (2015) A. Doostmohammadi, S. P. Thampi, T. B. Saw, C. T. Lim, B. Ladoux, and J. M. Yeomans, Soft Matter 11, 7328 (2015).
  * Driscoll  _et al._ (2017) M. Driscoll, B. Delmotte, M. Youssef, S. Sacanna, A. Donev, and P. Chaikin, Nature Physics 13, 375 (2017).
  * Keber  _et al._ (2014) F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, Science 345, 1135 (2014).
  * Marchetti  _et al._ (2013a) M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013a).
  * Needleman and Dogic (2017) D. Needleman and Z. Dogic, Nature Reviews Materials 2, 1 (2017).
  * Burmølle  _et al._ (2006) M. Burmølle, J. S. Webb, D. Rao, L. H. Hansen, S. J. Sørensen, and S. Kjelleberg, Applied and environmental microbiology 72, 3916 (2006).
  * Nadell  _et al._ (2015) C. D. Nadell, K. Drescher, N. S. Wingreen, and B. L. Bassler, The ISME journal 9, 1700 (2015).
  * Brugués  _et al._ (2014) A. Brugués, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Muñoz, G. W. Brodland, B. Ladoux, and X. Trepat, Nature Physics 10, 684 (2014).
  * Friedl and Gilmour (2009) P. Friedl and D. Gilmour, Nature Reviews Molecular Cell Biology 10, 445 (2009).
  * Friedl  _et al._ (2012) P. Friedl, J. Locker, E. Sahai, and J. E. Segall, Nature Cell Biology 14, 777 (2012).
  * Zhang  _et al._ (2020a) J. Zhang, R. Alert, J. Yan, N. S. Wingreen, and S. Granick, arXiv preprint arXiv:2011.03175 (2020a).
  * Balasubramaniam  _et al._ (2020) L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. Sankara, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, Y. Toyama, _et al._ , bioRxiv (2020).
  * Chaudhuri  _et al._ (2020) O. Chaudhuri, J. Cooper-White, P. A. Janmey, D. J. Mooney, and V. B. Shenoy, Nature 584, 535 (2020).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Fix (1982) G. Fix (1982).
  * (17) J. S. Langer, “Models of pattern formation in first-order phase transitions,” in  _Directions in Condensed Matter Physics_, pp. 165–186.
  * Emmerich (2008) E. Emmerich, Advances in Physics 57, 1 (2008).
  * Allen and Cahn (1972) S. Allen and J. Cahn, Acta Metallurgica 20, 423 (1972).
  * Cahn and Hilliard (1958) J. W. Cahn and J. E. Hilliard, The Journal of Chemical Physics 28, 258 (1958), https://doi.org/10.1063/1.1744102 .
  * Jacqmin (1999) D. Jacqmin, Journal of Computational Physics 155, 96 (1999).
  * Yue  _et al._ (2004) P. Yue, J. J. Feng, C. Liu, and J. Shen, Journal of Fluid Mechanics 515, 293 (2004).
  * Mauri  _et al._ (1996) R. Mauri, R. Shinnar, and G. Triantafyllou, Physical Review E 53, 2613 (1996).
  * Lin and Liu (1995) F.-H. Lin and C. Liu, Communications on Pure and Applied Mathematics 48, 501 (1995).
  * Sulaiman  _et al._ (2006) N. Sulaiman, D. Marenduzzo, and J. Yeomans, Physical Review E 74, 041708 (2006).
  * Blow  _et al._ (2014) M. L. Blow, S. P. Thampi, and J. M. Yeomans, Physical review letters 113, 248303 (2014).
  * Cates and Tjhung (2018) M. E. Cates and E. Tjhung, Journal of Fluid Mechanics 836, P1 (2018).
  * Kempf  _et al._ (2019) F. Kempf, R. Mueller, E. Frey, J. M. Yeomans, and A. Doostmohammadi, Soft matter 15, 7538 (2019).
  * LeVeque (2007) R. J. LeVeque, SIAM (2007).
  * Marenduzzo  _et al._ (2007) D. Marenduzzo, E. Orlandini, M. E. Cates, and J. M. Yeomans, Phys. Rev. E 76, 031921 (2007).
  * Desplat  _et al._ (2001) J.-C. Desplat, I. Pagonabarraga, and P. Bladon, Computer Physics Communications 134, 273 (2001).
  * Jülicher  _et al._ (2018) F. Jülicher, S. W. Grill, and G. Salbreux, Reports on Progress in Physics 81, 076601 (2018).
  * Doostmohammadi  _et al._ (2018) A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Nature Communications 9, 1 (2018).
  * Giomi and DeSimone (2014) L. Giomi and A. DeSimone, Physical review letters 112, 147802 (2014).
  * Turiv  _et al._ (2020) T. Turiv, R. Koizumi, K. Thijssen, M. M. Genkin, H. Yu, C. Peng, Q.-H. Wei, J. M. Yeomans, I. S. Aranson, A. Doostmohammadi, _et al._ , Nature Physics 16, 481 (2020).
  * Volfson  _et al._ (2008) D. Volfson, S. Cookson, J. Hasty, and L. S. Tsimring, Proceedings of the National Academy of Sciences 105, 15346 (2008).
  * Sanchez  _et al._ (2012) T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012).
  * Guillamat  _et al._ (2016) P. Guillamat, J. Ignés-Mullol, and F. Sagués, Proc. Natl. Acad. Sci. U.S.A 113, 5498 (2016).
  * Kumar  _et al._ (2018) N. Kumar, R. Zhang, J. J. de Pablo, and M. L. Gardel, Sci. Adv. 4, eaat7779 (2018).
  * Duclos  _et al._ (2018) G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost, and P. Silberzan, Nat. Phys. 14, 728 (2018).
  * Kawaguchi  _et al._ (2017) K. Kawaguchi, R. Kageyama, and M. Sano, Nature 545, 327 (2017).
  * Saw  _et al._ (2017) T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Nature 544, 212 (2017).
  * Blanch-Mercader  _et al._ (2018) C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi, and P. Silberzan, Physical Review Letters 120, 208101 (2018).
  * Yashunsky  _et al._ (2020) V. Yashunsky, D. Pearce, C. Blanch-Mercader, F. Ascione, L. Giomi, and P. Silberzan, arXiv preprint arXiv:2010.15555 (2020).
  * Voituriez  _et al._ (2005) R. Voituriez, J.-F. Joanny, and J. Prost, Europhy. Lett. (EPL) 70, 404 (2005).
  * Ramaswamy (2010) S. Ramaswamy, Annu. Rev. Condens. Matter Phys. 1, 323 (2010).
  * Marchetti  _et al._ (2013b) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013b).
  * Patteson  _et al._ (2018) A. E. Patteson, A. Gopinath, and P. E. Arratia, Nature communications 9, 1 (2018).
  * Zhang  _et al._ (2019) R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant, M. L. Gardel, _et al._ , arXiv preprint arXiv:1912.01630 (2019).
  * De Gennes and Prost (1993) P.-G. De Gennes and J. Prost,  _The physics of liquid crystals_ , Vol. 83 (Oxford university press, 1993).
  * Beris and Edwards (1994) A. N. Beris and B. J. Edwards,  _Thermodynamics of Flowing Systems_ (Oxford University Press, 1994).
  * Doostmohammadi  _et al._ (2016) A. Doostmohammadi, S. P. Thampi, and J. M. Yeomans, Phys. Rev. Lett. 117, 048102 (2016).
  * Blow  _et al._ (2017) M. L. Blow, M. Aqil, B. Liebchen, and D. Marenduzzo, Soft matter 13, 6137 (2017).
  * Carenza  _et al._ (2019a) L. N. Carenza, G. Gonnella, D. Marenduzzo, and G. Negro, Proceedings of the National Academy of Sciences 116, 22065 (2019a).
  * Coelho  _et al._ (2020) R. C. Coelho, N. A. Araújo, and M. M. T. da Gama, Soft Matter 16, 4256 (2020).
  * Neta  _et al._ (2020) P. Neta, M. Tasinkevych, M. T. da Gama, and C. Dias, arXiv preprint arXiv:2011.05963 (2020).
  * Dell’Arciprete  _et al._ (2018) Dell’Arciprete, M. Blow, A. Brown, F. Farrell, J. S. Lintuvuori, A. McVey, D. Marenduzzo, and W. C. Poon, Nat. Commun. 9, 4190 (2018).
  * Emmanuel  _et al._ (2020) L. Emmanuel, J. M. Yeomans, and A. Doostmohammadi, Physical Review Fluids 5, 023102 (2020).
  * Hemingway  _et al._ (2015) E. Hemingway, A. Maitra, S. Banerjee, M. C. Marchetti, S. Ramaswamy, S. M. Fielding, and M. E. Cates, Physical review letters 114, 098302 (2015).
  * Naganathan  _et al._ (2014) S. R. Naganathan, S. Fürthauer, M. Nishikawa, F. Jülicher, and S. W. Grill, Elife 3, e04165 (2014).
  * Zhang  _et al._ (2016) R. Zhang, Y. Zhou, M. Rahimi, and J. J. De Pablo, Nature communications 7, 1 (2016).
  * Guillamat  _et al._ (2018) P. Guillamat, Ž. Kos, J. Hardoüin, J. Ignés-Mullol, M. Ravnik, and F. Sagués, Science advances 4, eaao1470 (2018).
  * Metselaar  _et al._ (2019) L. Metselaar, J. M. Yeomans, and A. Doostmohammadi, Physical Review Letters 123, 208001 (2019).
  * Maroudas-Sacks  _et al._ (2020) Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, bioRxiv (2020).
  * Wittkowski  _et al._ (2014) R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo, and M. E. Cates, Nature communications 5, 1 (2014).
  * Negro  _et al._ (2019) G. Negro, L. N. Carenza, A. Lamura, A. Tiribocchi, and G. Gonnella, Soft matter 15, 8251 (2019).
  * Carenza  _et al._ (2019b) L. N. Carenza, G. Gonnella, A. Lamura, G. Negro, and A. Tiribocchi, The European Physical Journal E 42, 81 (2019b).
  * Ziebert  _et al._ (2011) F. Ziebert, S. Swaminathan, and I. S. Aranson, Journal of the Royal Society Interface 9, 1084 (2011).
  * Tjhung  _et al._ (2012) E. Tjhung, D. Marenduzzo, and M. E. Cates, Proceedings of the National Academy of Sciences 109, 12381 (2012).
  * Marth and Voigt (2016) W. Marth and A. Voigt, Interface Focus 6, 20160037 (2016).
  * Bi  _et al._ (2015) D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning, Nature Physics 11, 1074 (2015).
  * Alert and Trepat (2020) R. Alert and X. Trepat, Annual Review of Condensed Matter Physics 11, 77 (2020).
  * Graner and Glazier (1992) F. Graner and J. A. Glazier, Physical Review Letters 69, 2013 (1992).
  * Farhadifar  _et al._ (2007) R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, Current Biology 17, 2095 (2007).
  * Bi  _et al._ (2016) D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Physical Review X 6, 021011 (2016).
  * Sussman  _et al._ (2018) D. M. Sussman, M. Paoluzzi, M. C. Marchetti, and M. L. Manning, EPL (Europhysics Letters) 121, 36001 (2018).
  * Prost  _et al._ (2015) J. Prost, F. Jülicher, and J.-F. Joanny, Nature Physics 11, 111 (2015).
  * Löber  _et al._ (2015) J. Löber, F. Ziebert, and I. S. Aranson, Scientific Reports 5, 9172 (2015).
  * Camley and Rappel (2017) B. A. Camley and W.-J. Rappel, Journal of Physics D: Applied Physics 50, 113002 (2017).
  * Shaebani  _et al._ (2020) M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Nature Reviews Physics 2, 1 (2020).
  * Löber  _et al._ (2014) J. Löber, F. Ziebert, and I. S. Aranson, Soft Matter 10, 1365 (2014).
  * Ziebert and Aranson (2013) F. Ziebert and I. S. Aranson, PLoS one 8, e64511 (2013).
  * Najem and Grant (2016) S. Najem and M. Grant, Phys. Rev. E 93, 052405 (2016).
  * Camley  _et al._ (2014) B. A. Camley, Y. Zhang, Y. Zhao, B. Li, E. Ben-Jacob, H. Levine, and W.-J. Rappel, Proceedings of the National Academy of Sciences 111, 14770 (2014).
  * Palmieri  _et al._ (2015) B. Palmieri, Y. Bresler, D. Wirtz, and M. Grant, Scientific Reports 5, 11745 (2015).
  * Peyret  _et al._ (2019) G. Peyret, R. Mueller, J. d’Alessandro, S. Begnaud, P. Marcq, R.-M. Mège, J. M. Yeomans, A. Doostmohammadi, and B. Ladoux, Biophysical Journal 117, 464 (2019).
  * Wenzel  _et al._ (2019) D. Wenzel, S. Praetorius, and A. Voigt, The Journal of Chemical Physics 150, 164108 (2019).
  * Mueller  _et al._ (2019) R. Mueller, J. M. Yeomans, and A. Doostmohammadi, Physical Review Letters 122, 048004 (2019).
  * Loewe  _et al._ (2019) B. Loewe, M. Chiang, D. Marenduzzo, and M. C. Marchetti, arXiv preprint arXiv:1912.10549 (2019).
  * Giavazzi  _et al._ (2018) F. Giavazzi, M. Paoluzzi, M. Macchi, D. Bi, G. Scita, M. L. Manning, R. Cerbino, and M. C. Marchetti, Soft matter 14, 3471 (2018).
  * Kim  _et al._ (2013) J. H. Kim, X. Serra-Picamal, D. T. Tambe, E. H. Zhou, C. Y. Park, M. Sadati, J.-A. Park, R. Krishnan, B. Gweon, E. Millet, _et al._ , Nature Materials 12, 856 (2013).
  * Note (1) The area is proportional to the _square_ in order to ensure that it is always positive even when the phase field is slightly negative which can happen during simulations.
  * Deforet  _et al._ (2014) M. Deforet, V. Hakim, H. Yevick, G. Duclos, and P. Silberzan, Nature Communications 5, 3747 (2014).
  * Zhang  _et al._ (2020b) G. Zhang, R. Mueller, A. Doostmohammadi, and J. M. Yeomans, Journal of The Royal Society Interface 17, 20200312 (2020b).
  * Nonomura (2012) M. Nonomura, PLOS ONE 7, 1 (2012).
  * Sarangi  _et al._ (2017) B. R. Sarangi, M. Gupta, B. L. Doss, N. Tissot, F. Lam, R.-M. Mège, N. Borghi, and B. Ladoux, Nano Letters 17, 399 (2017).
  * Mitchison and Cramer (1996) T. Mitchison and L. Cramer, Cell 84, 371 (1996).
  * Alberts  _et al._ (2015) B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter, (2015).
  * Szabó  _et al._ (2006) B. Szabó, G. Szöllösi, B. Gönci, Z. Jurányi, D. Selmeczi, and T. Vicsek, Physical Review E 74, 061908 (2006).
  * Smeets  _et al._ (2016) B. Smeets, R. Alert, J. Pešek, I. Pagonabarraga, H. Ramon, and R. Vincent, Proceedings of the National Academy of Sciences 113, 14621 (2016).
  * Duclos  _et al._ (2017) G. Duclos, C. Erlenkämper, J.-F. Joanny, and P. Silberzan, Nat. Phys. 13, 58 (2017).
  * Asipauskas  _et al._ (2003) M. Asipauskas, M. Aubouy, J. A. Glazier, F. Graner, and Y. Jiang, Granular Matter 5, 71 (2003).