# Pebble trapping in vortices: three-dimensional simulations

Natalie Raettig  Max-Planck-Institut für Astronomie, Königstuhl 17, 69117,
Heidelberg, Germany  Wladimir Lyra New Mexico State University, Department of
Astronomy, PO Box 30001 MSC 4500, Las Cruces, NM 88001, USA  Hubert Klahr Max-
Planck-Institut für Astronomie, Königstuhl 17, 69117, Heidelberg, Germany

(Received Dec 23, 2020; Accepted Apr 10, 2021)

###### Abstract

Disk vortices have been heralded as promising routes for planet formation due
to their ability to trap significant amounts of pebbles. While the gas motions
and trapping properties of two-dimensional vortices have been studied in
enough detail in the literature, pebble trapping in three dimensions has
received less attention, due to the higher computational demand. Here we use
the Pencil Code to study 3D vortices generated by convective overstability and
the trapping of solids within them. The gas is unstratified whereas the
pebbles settle to the midplane due to vertical gravity. We find that for
pebbles of normalized friction times of $\mathrm{St}=0.05$ and
$\mathrm{St}=1$, and dust-to-gas ratio $\varepsilon=0.01$, the vortex column
in the midplane is strongly perturbed. Yet, when the initial dust-to-gas ratio
is decreased the vortices remain stable and function as efficient pebble
traps. Streaming instability is triggered even for the lowest dust-to-gas
ratio ($\varepsilon_{0}=10^{-4}$) and smallest pebble sizes
($\mathrm{St}=0.05$) we assumed, showing a path for planetesimal formation in
vortex cores from even extremely subsolar metallicity. To estimate if the
reached overdensities can be held together solely by their own gravity we
estimate the Roche density at different radii. Depending on disk model and
radial location of the pebble clump we do reach concentrations higher than the
Roche density. We infer that if self-gravity was included for the pebbles then
gravitational collapse would likely occur.

††journal: ApJ

##  1 Introduction

Vortices, hydrodynamical flow in closed elliptic streamings, have long been
considered as a possible route for planet formation since Barge & Sommeria
(1995); Adams & Watkins (1995) and Tanga et al. (1996) independently suggested
that they would be exceptional sites for trapping solids aerodynamically. The
shear in a Keplerian disk means that only anticyclonic vortices survive; being
high pressure regions, these type of vortices will behave as dust traps
(Whipple, 1972; Bracco et al., 1999; Chavanis, 2000), concentrating the
marginally coupled pebbles that would otherwise drift toward the star
(Weidenschilling, 1980; Klahr & Bodenheimer, 2006).

Dust trapping in 2D was subsequently studied by Inaba & Barge (2006), finding
that solids rapidly sink toward the center of the vortices and increase the
solids-to-gas ratio by at least an order of magnitude. Lyra et al. (2008a,
2009a, 2009b) performed global models, also in 2D, and including selfgravity,
finding that the concentration of solids collapsed to objects of lunar and
Mars-mass, that subsequently undergo further collisions and pebble accretion,
building up Earth-mass planets.

Because these results were obtained in 2D, where the turbulent cascade is
inverse, the question remained whether a vortex could be formed and sustained
in 3D, where the vortices spontaneously develop internal turbulence (Lesur &
Papaloizou, 2009), leading to direct cascade and vortex destruction. Lyra
(2009) presented evidence that vortices would be sustained in 3D via the
Rossby wave instability (RWI, Lovelace et al., 1999) in special locations,
namely the boundary between viscosity transitions or the gaps opened by a
planet. Sustenance of Rossby vortices in 3D was subsequently demonstrated by
Meheut et al. (2010, 2012a, 2012b, 2012c); Lin (2012a, b); Lyra & Mac Low
(2012); Lin (2013) and Lyra et al. (2015). Beyond the RWI, Lesur & Papaloizou
(2010) and Lyra & Klahr (2011) also showed that in non-isothermal disks, a
mechanism for vortex formation exists in the bulk of the gas, a baroclinic
instability, first thought to be nonlinear, but then found to be linear, and
renamed Convective Overstability (COV, Klahr & Hubbard, 2014; Lyra, 2014).
Other non-magnetized hydrodynamical instabilities in the disk Ohmic zone, the
zombie vortex instability (Marcus et al., 2015, 2016; Barranco et al., 2018)
and the vertical shear instability (Nelson et al., 2013; Stoll & Kley, 2014;
Richard et al., 2016; Manger & Klahr, 2018; Flock et al., 2020; Manger et al.,
2020; Pfeil & Klahr, 2020) also show the formation of large scale vortices in
3D (see Lyra & Umurhan, 2019, for a review).

While these results from the last decade of research on gas dynamics in
protoplanetary disks decisively settled the question of vortex formation and
sustenance in 3D, two-dimensional simulations with dust cast doubt on the
survival of dusty vortices once the dust feedback becomes too strong and
disrupts the elliptic streamlines (Fu et al., 2014; Raettig et al., 2015;
Surville et al., 2016; Miranda et al., 2017). This is a phenomenon that had
already been noted by Inaba & Barge (2006), that the lifetime of vortices was
lowered as the dust accumulation increased. This has lead to the suggestion
that the asymmetries seen in ALMA dust continuum observations of transition
disks (van der Marel et al., 2013, 2020) are not vortices, despite the fact
that the observational results are generally consistent with the analytical
predictions for vortex-trapped dust in steady-state between drag and diffusion
(Klahr & Henning, 1997; Lyra & Lin, 2013; Sierra et al., 2017; Casassus et
al., 2019). The question therefore remains whether dusty vortices survive in
3D.

To answer this question, we conduct three-dimensional unstratified gas
simulation of the convective overstability, as in Lesur & Papaloizou (2010);
Lyra & Klahr (2011) and Lyra (2014), but with pebbles that are allowed to
sediment due to stellar gravity. As we will see, the dust vertical
stratification dramatically changes the result of two-dimensional models,
since pebbles sediment with a different scale height than the gas, undergo
vertical motion and are subject to vertical stirring. Another important factor
is how the pebbles affect the vortex structure, which we already saw in our 2D
simulations, most strongly for loosely coupled grains (Raettig et al., 2015).
In 3D with sedimentation, three different outcomes are possible:

  1. 1.

The vortex column is stable and pebbles accumulate inside the column.

  2. 2.

The pebbles disrupt the vortex column in the layer they reside around the
midplane, but do not affect the vortex column above and below the pebble
layer. The question is then how the remaining vortex column will affect the
pebbles concentration.

  3. 3.

The entire vortex column is destroyed by the pebbles. Pebbles settle to the
midplane with no significant difference with respect to non-turbulent runs.
The baroclinic vortex cannot be reestablished.

As shown in Lyra et al. (2018), the second option is the dominant outcome.
Here we detail and expand this result; the phenomenology depends strongly on
the pebble radius and the dust-to-gas ratio. First we will introduce in Sect.
2 the new physics that is needed for three-dimensional estimates. The
numerical setup is presented in Sect. 3. In Sect. 4 we will analyze the
evolution of pebbles in the simulations and, finally, discuss our results.

Figure 1: Vertical vorticity of unstratified gas in a three-dimensional box at
different times. The convective overstability amplifies seed noise into
growing vortices, that intensity and merge viscously, saturating at radial
length similar to the pressure scale height.

Figure 2: The vertical component of the vorticity ($\omega_{z}$) in the
midplane (left panels) and a vertical cut through the disk at the azimuthal
location of maximum negative vorticity (indicated by the dashed line). Also
shown is the vertical velocity component ($u_{z}$). The upper plots are shown
at 100 local orbits (top panels). At this time, there is no structure in the
vertical direction: the flow is essentially two-dimensional. After 275 local
orbits (bottom panels) there is now turbulent motion inside the vortex, which
is also visible in the vertical velocity. This is an indication of the growth
of elliptical instability in the vortex cores.

##  2 Model Equations

We model a local patch of the disk following the shearing box approximation.
The gas is rendered unstable to the COV by introducing baroclinicity in the
box, in the form of a linearized pressure gradient. The equations of motion
are

$\displaystyle\frac{{\mathcal{D}}\rho_{{g}}}{{\mathcal{D}}t}$ $\displaystyle=$
$\displaystyle-\rho_{{g}}{\bm{\nabla}}\cdot{\bm{u}}$ (1)
$\displaystyle\frac{{\mathcal{D}}{\bm{u}}}{{\mathcal{D}}t}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho_{{g}}}{\bm{\nabla}}{p}-2\varOmega\left(\hat{{\bm{z}}}\times{\bm{u}}\right)+\frac{3}{2}\varOmega
u_{x}\hat{{\bm{y}}}+\frac{\rho_{d}}{\rho_{{g}}}\frac{\left({\bm{v}}-{\bm{u}}\right)}{\tau}$
(2)
$\displaystyle+f(\beta)\left(\frac{1}{\rho_{{g}}}-\frac{1}{\rho_{{g}0}}\right)\hat{{\bm{x}}}$
$\displaystyle\frac{{\mathcal{D}}s}{{\mathcal{D}}t}$ $\displaystyle=$
$\displaystyle\frac{1}{\rho_{{g}}T}\left[{\bm{\nabla}}\cdot\left(K{\bm{\nabla}}{T}\right)-\rho_{{g}}c_{{}_{V}}\frac{\left(T-T_{0}\right)}{\tau_{c}}+f(\beta)\frac{u_{x}}{(\gamma-1)}\right]$

In the equations above, $\rho_{g}$ is the gas density ${\bm{u}}$ is the
velocity, ${\varOmega}$ is the Keplerian frequency, $\rho_{d}$ is the dust
density, ${\bm{v}}$ is the dust velocity, $\tau$ is the friction time between
gas and dust, $s$ is the entropy, $T$ is the temperature, $K$ is the radiative
conductivity, $c_{{}_{V}}$ is the heat capacity at constant volume,
$\gamma=c_{p}/c_{{}_{V}}$ is the adiabatic index, $c_{p}$ is the heat capacity
at constant pressure, $T_{0}$ is the gas equilibrium temperature, and
$\tau_{c}$ is the cooling time. The operator

$\frac{{\mathcal{D}}}{{\mathcal{D}}t}\equiv\frac{\partial{}}{\partial{t}}{+{\bm{u}}\cdot{\bm{\nabla}}}+u_{{k}}\frac{\partial{}}{\partial{y}}$
(4)

represents the Lagrangian derivative of a fluid parcel. It is the only place
where the (linearized) Keplerian flow
${\bm{u}}_{{k}}{=-\nicefrac{{3}}{{2}}\,\varOmega x\,\hat{y}}$ appears
explicitly. The function $f(\beta)$ is

$f(\beta)\equiv\frac{\beta p_{0}\varOmega}{c_{s0}}$ (5)

where the parameter $\beta\equiv h\xi$, (with $h=c_{s}/u_{k}$ the disk aspect
ratio, and $c_{s}$ the sound speed), is related to the power law of the
pressure gradient

$\bar{p}(r)\equiv p_{0}(r/r_{0})^{-\xi},$ (6)

where $r$ is the cylindrical radius and $r_{0}$ is a reference radius. The
overbar indicates that this quantity is time-independent. The total pressure
is $p_{\rm tot}$=$\bar{p}+p$, where $p$ is the local fluctuation. Linearizing
this pressure gradient contributes an extra term in the momentum equation and
in the energy equation, respectively. The derivation of this extra term and
the limitations of the approach are detailed in Lyra & Klahr (2011). All the
symbols used in this work are listed in

table:symbols.

The advection is made Galilean-invariant by means of the SAFI algorithm
(Johansen et al., 2009). The simulations were done with the Pencil Code 111See
http://pencil-code.nordita.org (Brandenburg & Dobler, 2002, 2010; Brandenburg
et al., 2020), a collocated finite-difference code which integrates the
partial differential equations with sixth-order spatial derivatives, and a
third-order Runge–Kutta time integrator. We add explicit sixth-order
hyperdiffusion, hyperviscosity, and hyper heat conductivity to the mass,
momentum, and entropy equations, respectively, explained in Lyra et al.
(2008b, 2009a, 2017). They are needed because the high-order scheme of the
Pencil Code has little overall numerical dissipation (McNally et al., 2012).
We use dimensionless code units; with all simulations using
$c_{p}=\varOmega=\rho_{{g}0}=1$, $\gamma=1.4$, $c_{s0}=h=0.1$, and $\xi=2$.
The models are scale-free and, thus, the translation to physical units depends
on specific parameter choices, as well as on the distance from the box to the
star. For densities of the Minimum-Mass Solar Nebula, (MMSN, Weidenschilling,
1977; Hayashi, 1981), $\Sigma_{g}=1700(r/{\rm AU})^{-1.5}{\rm g\,cm}^{-2}$ and
temperature of 100 K at 5AU and 30 K at 50 AU, the disk is stable against
gravitational instabillity (Safronov, 1960; Toomre, 1964), having a Toomre $Q$
value of 50 at 5 AU and 20 at 50AU.

###  2.1 Pebble evolution

The pebbles are represented by numerical particles that evolve Lagrangianly
according to

$\displaystyle\frac{d{\bm{x}}}{dt}$ $\displaystyle=$ $\displaystyle{\bm{v}}$
(7) $\displaystyle\frac{d{\bm{v}}}{dt}$ $\displaystyle=$
$\displaystyle-\frac{\left({\bm{v}}-{\bm{u}}\right)}{\tau}-2\varOmega\left(\hat{{\bm{z}}}\times{\bm{v}}\right)+\frac{3}{2}\varOmega
v_{x}\hat{{\bm{y}}}-\varOmega^{2}z\hat{{\bm{z}}}$

where ${\bm{x}}$ is the pebble position and ${\bm{v}}$ is the pebble velocity
corrected by the sub-Keplerian parameter

$\eta\equiv\frac{1}{2}h^{2}\xi{.}$ (9)

The extra velocity $\Delta v=\eta u_{k}$ is added to the pebbles, as explained
in Raettig et al. (2015). This term is needed because the frame of reference
is sub-Keplerian, as ensured by the extra $1/\rho_{{g}0}$ term in Eq. (2). In
the reference frame of the box, the sub-Keplerian gas is stationary, and the
pebbles flow with $\eta u_{k}$ azimuthal velocity. This is constrated with
simulations by Johansen et al. (2006) where the reference frame is Keplerian,
i.e. the pebbles are stationary and the gas flows with $-\eta u_{k}$ azimuthal
velocity. Here we use $\eta=0.01$, which corresponds to a pressure gradient
$\xi=2$ and disk aspect ratio $h=0.1$. The terminal velocities of the pebbles
for a given Stokes number relative to the Keplerian flow are

$\displaystyle v_{x}$ $\displaystyle=$ $\displaystyle-\frac{2\eta
u_{{k}}}{\mathrm{St}+\mathrm{St}^{-1}}$ (10) $\displaystyle v_{y}$
$\displaystyle=$ $\displaystyle-\frac{\eta u_{{k}}}{{1+}\mathrm{St}^{2}}$ (11)

In the vertical direction, the pebbles settle with a constant velocity so that
$\partial_{t}v_{z}\approx 0$. This leads to

$\frac{\partial{z}}{\partial{t}}=-\tau\varOmega^{2}z$ (12)

and

$z{(t)}=z_{0}e^{-t\tau\varOmega^{2}}.$ (13)

Here pebbles will approach the midplane exponentially but never reach it. For
pebbles with $\mathrm{St}\gg 1$ the term $\tau^{-1}v_{z}$ is negligible and
the equation of motion reduces to

$\frac{\partial{v_{z}}}{\partial{t}}=-\varOmega^{2}z$ (14)

where

$z{(t)}=\cos(\varOmega t)$ (15)

is a solution. Heavy solids will oscillate around the midplane, never settling
permanently. In a real protoplanetary disk there are additional effects acting
on the solids that will change this idealized picture. The sedimentation time
scale $t_{s}$ for small grains can be deduced from Eq. (13). After a time
$t_{s}=(\varOmega\,\mathrm{St})^{-1}$ a pebble will have sedimented by a
factor $e$ of its original height. The pebble stirring works as an effective
diffusion, so that the evolution of the pebble density can be described by

$\frac{\partial{\rho_{d}}}{\partial{t}}=D\frac{\partial{}}{\partial{z}}\left[\rho_{g}\frac{\partial{}}{\partial{z}}\left(\frac{\rho_{d}}{\rho_{g}}\right)\right]+\frac{\partial{}}{\partial{z}}\left(\varOmega^{2}\tau_{s}\rho_{d}z\right),$
(16)

where $D$ is the diffusion constant. For steady-state (Dubrulle et al., 1995;
Klahr & Henning, 1997; Lyra & Lin, 2013), the pebble distribution settles as a
Gaussian around the midplane, with pebble scale height

$H_{d}=H\sqrt{\frac{\delta}{\mathrm{St}+\delta}}$ (17)

where $H\equiv hr_{0}\equiv c_{s}/\varOmega$ is the gas scale height and
$\delta\equiv D/(c_{s}H)$ is a dimensionless diffusion coefficient similar to
the Shakura-Sunyaev $\alpha$-viscosity coefficient (Shakura & Sunyaev, 1973).
The ratio of $D$ and $\nu$ defines the dimensionless Schmidt number

$\mathrm{Sc}\equiv\frac{\nu}{D}=\frac{\alpha}{\delta}$ (18)

The pebble scale height then is

$H_{d}=H\sqrt{\frac{\alpha}{\mathrm{St}\,\mathrm{Sc}+\alpha}}$ (19)

Table 1: Symbols used in this work. Symbol | Definition | Description |  | Symbol | Definition | Description
---|---|---|---|---|---|---
$\rho_{g}$ |  | gas density |  | $\xi$ | $d\ln p/d\ln r$ | pressure gradient
${\bm{u}}$ |  | gas velocity |  | $H$ | $c_{s}/\varOmega$ | gas scale height
$p$ | $\rho_{g}c_{s}^{2}/\gamma$ | gas pressure |  | $h$ | $H/r_{0}$ | gas aspect ratio
$s$ | $c_{{}_{V}}[\ln(\nicefrac{{T}}{{T_{0}}})-(\gamma-1)\ln(\nicefrac{{\rho}}{{\rho_{0}}})]$ | gas entropy |  | $\beta$ | $h\xi$ | scaled pressure gradient
$T$ |  | gas temperature |  | ${\bm{\omega}}$ | ${\bm{\nabla}}\times{{\bm{u}}}$ | gas vorticity
$K$ |  | radiative conductivity |  | $\eta$ | $h^{2}\xi/2$ | sub-Keplerian parameter
$c_{s}$ | $[Tc_{p}(\gamma-1)]^{1/2}$ | sound speed |  | $\mathrm{St}$ | $\varOmega\tau$ | Stokes number
$c_{{}_{V}}$ |  | Specific heat at constant volume |  | $t_{s}$ |  | sedimentation time
$c_{p}$ |  | Specific heat at constant pressure |  | $\delta$ |  | dimensionless diffusion coefficient
$\gamma$ | $c_{p}/c_{{}_{V}}$ | Adiabatic index |  | $D$ | $\delta c_{s}H$ | diffusion coefficient
$r_{0}$ |  | Reference radius |  | $H_{d}$ | $H\sqrt{\delta/(\delta+\mathrm{St})}$ | dust scale height
$\varOmega$ |  | Keplerian frequency at $r_{0}$ |  | $\alpha$ |  | dimensionless viscosity coefficient
${\bm{u}}_{k}$ | $-\nicefrac{{3}}{{2}}\,\varOmega x\,\hat{y}$ | Keplerian velocity at $r_{0}$ |  | $\nu$ | $\alpha c_{s}H$ | viscosity coefficient
$r$ |  | cylindrical radius |  | $\mathrm{Sc}$ | $\nu/D$ | Schmidt number
$x$ | $r-r_{0}$ | Cartesian radial coordinate |  | $\varepsilon$ | $\rho_{d}/\rho_{g}$ | dust-to-gas ratio
$y$ |  | Cartesian azimuthal coordinate |  | $k$ |  | wavenumber
$z$ |  | vertical coordinate |  | $L$ |  | box length
$t$ |  | time |  | $\Sigma_{g}$ | $\int\rho_{g}dz$ | gas column density
$\rho_{d}$ |  | dust density |  | $\Sigma_{d}$ | $\int\rho_{d}dz$ | dust column density
${\bm{v}}$ |  | dust velocity |  | $\varphi$ | $[0,1]$ | Random phase
$\tau$ |  | drag time |  | $a_{\bullet}$ |  | pebble radius
$\tau_{c}$ |  | cooling time |  | $\rho_{\bullet}$ |  | pebble internal density
$\rho_{R}$ | $3\pi M_{\odot}/(2r^{3})$ | Roche density |  |  |  |
Table 2: Simulation setup and results Run | $\mathrm{St}$ | $\varepsilon_{0}$ | feedback | ${\xi}$ | $H_{d}$ | $D$ | $\delta$ | $\alpha$ | $\mathrm{Sc}_{z}$ | $t_{\rm end}$ | $\varepsilon_{\rm max}$ | ${\zeta}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  | $(\times 10^{3})$ | $(\times 10^{7})$ | $(\times 10^{5})$ | $(\times 10^{3})$ |  |  |  |
3DG | - | - | - | 2 | - | - | - | - | - | - | - | -
3DF05nt | 0.05 | $10^{-2}$ | yes | 0 | 1.20 | 0.72 | 0.72 | 0.1 | 13.87 | 43.4 | 41.6 | -
3DNF05 | 0.05 | $10^{-2}$ | no | 2 | 2.57 | 3.31 | 3.31 | 1.12 | 33.83 | 40.7 | 3540.0 | -
3DF05 | 0.05 | $10^{-2}$ | yes | 2 | 2.23 | 2.49 | 2.49 | 0.66 | 26.7 | 38.2 | 111.2 | 98.44
3DF05E-3 | 0.05 | $10^{-3}$ | yes | 2 | - | - | - | - | - | 25.3 | 8.7 | 25.67
3DF05E-4 | 0.05 | $10^{-4}$ | yes | 2 | - | - | - | - | - | 20.0 | 2.1 | 1.12
3DNF1 | 1.0 | $10^{-2}$ | no | 2 | 0.99 | 9.84 | 9.84 | 0.94 | 9.51 | 22.4 | 15630.0 | -
3DF1 | 1.0 | $10^{-2}$ | yes | 2 | 0.85 | 7.20 | 7.20 | 0.26 | 3.72 | 22.8 | 945.4 | 98.73
3DF1E-3 | 1.0 | $10^{-3}$ | yes | 2 | - | - | - | - | - | 7.7 | 227.7 | 26.41
3DF1E-4 | 1.0 | $10^{-4}$ | yes | 2 | - | - | - | - | - | 9.7 | 34.2 | 14.46

##  3 Numerical Setup

###  3.1 Gas Setup

Since we do not include vertical gravity for the gas simulations, the gas
setup is the same as in Raettig et al. (2015), with initial noise in the
density according to

$\rho_{{g}}=C\rho_{{g}0}e^{-(x/2\sigma)^{2}}\sum_{i=-k_{x}}^{k_{x}}\sum_{j=0}^{k_{y}}\sin\left[2\pi\left(i\frac{x}{L_{x}}+j\frac{y}{L_{y}}+\varphi_{ij}\right)\right]$
(20)

where $L$ is the box size, and the phase $0<\varphi<1$ determines the
randomness. The subscripts underscore that the phase is the same for all grid
points, only changing with wavenumber. The constant $C$ sets the strength of
the perturbation. The width $\sigma$ is set to $L_{x}/5$. In this setup no
vertical modes are excited, but radial and azimuthal modes are excited in all
vertical layers of the box. Additionally, we add Gaussian noise to the
velocity. The physical size of our box is the same in radial ($\pm 2H$) and
azimuthal (${\pm 8}H$) direction as in the 2D runs in Raettig et al. (2013,
2015). We model $\pm 1H$ in the vertical direction. This is a compromise
between an adequately large physical domain to capture the elliptical
instability inside the vortex cores, but not high enough resolution to resolve
all modes that contribute to the streaming instability. The vertical
resolution is the same as in the radial direction. We restrict ourselves to an
entropy gradient of $\xi=2$, since we established in Raettig et al. (2013)
that the main difference between $\xi=2$ and lower entropy gradients was the
time to reach saturation. The size of the vortices and strength of relevant
parameters such as kinetic stresses depend only weakly on $\xi$. For $\xi=2$ a
resolution of $256\times 256\times 128$ is appropriate to resolve the COV to
convergence (Lyra & Klahr, 2011; Raettig et al., 2013; Lyra, 2014).

###  3.2 Pebble Setup

For the pebbles, we add the vertical equation of motion, also including linear
vertical gravity $g_{z}=-\varOmega^{2}z$. As in the 2D simulations we first
evolve the gas alone for 200 local orbits before adding the pebbles. In these
simulations, the pebbles are randomly distributed in the $x-y$ plane as
before, but now they follow a Gaussian profile in the vertical direction

$\rho_{d}(z)=C_{d}\ e^{-\frac{z^{2}}{2{H_{d0}}^{2}}}$ (21)

where $\rho_{d}$ is the dust density, $C_{d}$ a normalization constant and
${H_{d0}}$ gives the scale height of the initial distribution. We set
${H_{d0}}=0.01$, which gives a well settled profile, and speeds up the initial
relaxation. Like in the 2D simulations (Raettig et al., 2015) we want to have
about $4-5$ super-particles per grid point initially, to mitigate numerical
effects due to Poisson noise (Johansen et al., 2007). Because most pebbles
will settle into the midplane, we only need to take into account $\pm$10 grid
points around the midplane in this estimate, leading to a particle count of
$7\times 10^{6}$.

Pebble feedback onto the gas is included for all simulations except for the
control runs. We restrict ourselves to two pebble sizes, $\mathrm{St}=1$ for
loosely coupled pebbles and $\mathrm{St}=0.05$ for more strongly coupled
pebbles. The Stokes number scales with grain size as

$\mathrm{St}=\frac{a_{\bullet}\rho_{\bullet}}{H\rho_{g}}$ (22)

where $a_{\bullet}$ is the grain radius and $\rho_{\bullet}$ its internal
density. In the simulations we fix $\mathrm{St}$, not the grain radius, so
effectively $a_{\bullet}$ varies dynamically with $\rho_{g}$ and $H$ due to
compressibility. However, in practice, the vortices are not too compressible,
with density variations of 50% at most. For the MMSN and assuming silicate
density ($\rho_{\bullet}=3\,{\rm g\,cm}^{-3}$), the corresponding grain sizes
for $\mathrm{St}=0.05$ at 1, 5, and 50 AU are 10cm, 1cm, and 0.3 mm. For ice
internal density ($\rho_{\bullet}=1\,{\rm g\,cm}^{-3}$), the sizes would be
scaled by about 3, so about 30cm, 3cm, and 1 mm. For $\mathrm{St}=1.0$, the
grain radii would be 2m, 20cm, and 0.6cm at 1, 5, and 50 AU for silicate
density, and 6m, 60cm, and 2cm for ices at the same distances. For a nebula 5
times more massive than MMSN (Toomre $Q\sim 10$ at 5 AU and 5 at 50 AU,
marginally stable at these distances) the sizes are scaled up by 5. In short,
the $\mathrm{St}=0.05$ particles represent cm-sized pebbles in the inner disk
and mm-sized in the outer disk, whereas $\mathrm{St}=1.0$ particles represent
m-sized boulders in the inner disk and cm-sized pebbles in the outer disk.
While coagulation/fragmentation equilibrium does not predict the existence of
m-sized boulders (Güttler et al., 2010; Zsom et al., 2010), the
$\mathrm{St}=1$ particles representing cm-sized range in the outer disk is a
pebble size important for planet formation and captured by VLA observations
(Casassus et al., 2019).

We expect the setup to develop $\alpha\approx 10^{-3}$ (Lyra, 2014). If we
assume for now that the Schmidt number for our system is 1, then the pebble
scale height for $\mathrm{St}=0.05$ is $H_{d}=0.01$. Therefore, a pebble layer
of $2H_{d}$ around the midplane should be resolved with 12 grid points. The
layer of $\mathrm{St}=1$ will not be resolved. We set the initial dust-to-gas
ratio to $\varepsilon_{0}=0.01$. Additionally, we perform runs with
$\varepsilon_{0}=10^{-3}$ and $\varepsilon_{0}=10^{-4}$. We also perform one
simulation with $\xi=0$ for $\mathrm{St}=0.05$ for a control run without
turbulence. All setups are listed in Table 2. Particle block domain
decomposition is used for load balance (Johansen et al., 2011).

Figure 3: Gas velocity components averaged over the entire box. As the
elliptical instability sets in, the vertical velocity component (green line)
increases steeply, while staying well below the azimuthal and radial
components.

##  4 Results

Figure 4: Maximum dust-to-gas ratio for simulations with $\mathrm{St}=1$
(black line) and $\mathrm{St}=0.05$ (green and red line) pebbles. $\xi=0$
indicates that no baroclinic feedback was included in the simulation. Although
the high pebble concentration in the simulations with baroclinic effects
disrupts the vortex in the midplane, where pebbles are located, the maximum
dust-to-gas ratio is higher than without baroclinic effects.

Figure 5: Vertical vorticity component $\omega_{z}$ at the top of the box
(left) and in the midplane (middle) and dust-to-gas ratio $\varepsilon$
(right) in the midplane for $\mathrm{St}=0.05$ pebbles and
$\varepsilon_{0}=10^{-2}$. The vortex in the midplane is not apparent, and the
pebbles are spread out. Since there was an initial pebble accumulation inside
the vortex before it was disrupted, there is a residual accumulation at the
original vortex position. The maximum dust-to-gas-ratio in this snapshot is
$\varepsilon=31.7$.

Figure 6: Vertical vorticity component $\omega_{z}$ and dust-to-gas ratio
$\varepsilon$ in the midplane for $\mathrm{St}=0.05$ pebbles and
$\varepsilon_{0}=10^{-4}$. In contrast to simulations with
$\varepsilon_{0}=10^{-2}$ the vortices in the midplane are not disrupted and
the pebbles accumulate inside of the two vortices. The maximum dust-to-gas-
ratio reached in this snapshot is $\varepsilon=2.2$ which is one order of
magnitude lower than for $\varepsilon_{0}=10^{-2}$.

The evolution of the vortex in a 3D unstratified box is very similar to the
two-dimensional case. First, a number of small vortices emerge. These vortices
are then amplified by the convective overstability mechanism (Lesur &
Papaloizou, 2010; Lyra & Klahr, 2011; Klahr & Hubbard, 2014; Lyra, 2014) and
also merge. Fig. 1 shows snapshots of the vorticity in the fiducial model. At
some point the vertical gas motion sets in and generates turbulent features in
the vorticity structure (compare top and bottom panel of Fig. 2). After about
100 local orbits there is a steep increase in vertical gas velocity, due to
elliptical instability (Kerswell, 2002; Lesur & Papaloizou, 2009, 2010; Lyra &
Klahr, 2011; Lyra, 2013). Yet it remains $2-3$ orders of magnitude lower than
the radial and azimuthal velocity components (Fig. 3).

Figure 7: Vertical cut through the vortex after 15 local orbits for
$\mathrm{St}=0.05$ pebbles and $\varepsilon_{0}=10^{-2}$ (top) and
$\varepsilon_{0}=10^{-4}$ (bottom). Shown are the vertical component of
vorticity $\omega_{z}$ (left) and the dust-to-gas ratio $\rho_{d}/\rho_{g}$
(right). For $\varepsilon_{0}=10^{-2}$ the vortex located around
$x\approx-0.1H$ is strongly perturbed around the midplane. For
$\varepsilon_{0}=10^{-4}$ the vortex column remains roughly undisturbed.

Figure 8:  Maximum dust-to-gas ratio for simulations with $\mathrm{St}=0.05$
(to panel) and $\mathrm{St}=1$ (bottom panel) pebbles, and different initial
dust-to-gas ratios $\varepsilon_{0}$. The concentration $\rho_{d}/\rho_{d0}$
seem to be only weakly dependent of $\varepsilon_{0}$.

Figure 9:  The top row shows the dust-to-gas ratio at a meridional cut through
the box at the location of the highest dust-density. The bottom row shows the
time development of the dust surface density $\Sigma_{d}$ relative to its
initial value $\varSigma_{d0}$ and averaged over the entire azimuthal domain.
On the left side we show the simulation of $\mathrm{St}=1$ pebbles with
$\varepsilon_{0}=10^{-2}$ and on the right side $\mathrm{St}=1$ pebbles with
$\varepsilon_{0}=10^{-3}$. The location of the pebble concentration agrees
with the location of the vortices. The pebble concentration migrates inwards
with the vortex (even in the case where the vortex in the midplane is
disrupted). The oscillations in the pebble concentration correspond to
epicyclic motion of the pebbles within the vortex.

Figure 10: Same as Fig. 9, but for $\mathrm{St}=0.05$ pebbles. The vertical
column of pebbles is resolved. Pebbles migrate with the vortex column, showing
trapping is efficient even if the vortex column is disrupted around the
midplane.

Figure 11: Fraction of total dust mass with dust-to-gas ratio higher than a
certain value for $\mathrm{St}=0.05$ and $\mathrm{St}=1.0$ pebbles and
different initial dust-to-gas ratios. For all setups the threshold for
streaming instability ($\rho_{d}/\rho_{g}\geq 1$) is passed.

###  4.1 Pebble evolution

We first discuss the run without baroclinic driving. As expected, the pebbles
settle to the midplane, soon triggering the streaming instability. Although we
assumed an overall initial dust-to-gas ratio of $\varepsilon_{0}=0.01$, the
initial dust-to-gas ratio in the midplane is already $\varepsilon=1$ and
increases even further as the sedimentation progresses. The maximum dust-to-
gas ratio of the simulation with $\mathrm{St}=0.05$ pebbles and without
baroclinic driving can be seen in Fig. 4 (red line). The other two lines in
the plot show the maximum dust-to-gas ratios for the simulation with $\xi=2$,
$\mathrm{St}=0.05$ (green line) and $\mathrm{St}=1$ pebbles (black line).

The maximum concentration for $\mathrm{St}=0.05$ with baroclinic driving is
only a factor of two higher than in the case without it. Inspecting the
evolution of the dust density in the baroclinic simulation, we see that, as
expected, the pebbles settle to the midplane and get trapped in the vortices.
Yet, after only 3 orbits, the pebbles start to perturb the vortices in the
midplane via the streaming instability. Like in the 2D simulations of Raettig
et al. (2015), high pebble densities generate local perturbations of the
velocity field visible as steep vorticity gradients. Fig. 5 shows the
vorticity at the top of the box (left panel), in the midplane (middle panel)
and the dust-to-gas ratio in the midplane for $\mathrm{St}=0.05$. We see that
although the vortex in the midplane is heavily disrupted, the vortex column
above and below the midplane is still present. Since there are no pebbles in
the upper and lower areas of the box, there is no disruption, and we expect
these vortex columns to stay stable over long times.

As we will discuss later, these vortex columns still have a minor influence on
the pebbles. The explanation for why we do not see a high increase of dust-to-
gas ratio for simulations with baroclinic vortices compared to simulations
without baroclinic driving lies in the disruption of vortices in the midplane.
To test whether there is a strong dust concentration when vortices are stable,
and also to account for the fact that not all pebbles will be of the same size
we perform simulations with lower initial dust-to-gas ratios,
$\varepsilon_{0}=10^{-3}$ and $\varepsilon_{0}=10^{-4}$.

Already for $\varepsilon_{0}=10^{-3}$ we see a different picture. The pebbles
settle to the midplane, but the vortices there are no longer disrupted: the
entire vortex column is a stable feature. Pebbles accumulate inside of the
vortices, as we saw in the 2D simulations, and then migrate slowly with them.
Fig. 6 shows the midplane vertical vorticity (left panel) and dust-to-gas
ratio (right panel) for $\mathrm{St}=0.05$ and $\varepsilon_{0}=10^{-4}$ after
15 orbits. Strong dust concentration is achieved in the vortex, yet no strong
vortex disruption in the midplane is observed. Fig. 7 shows a vertical cut of
the vertical vorticity (left panels) and dust-to-gas ratio (right panel),
again at 15 orbits, for the $\varepsilon_{0}=10^{-2}$ simulation (upper panel)
and the $\varepsilon_{0}=10^{-4}$ simulation (lower panel). For
$\varepsilon_{0}=10^{-2}$ the vortex located around $x\approx-0.1H$ is
strongly perturbed around the midplane. For $\varepsilon_{0}=10^{-4}$ the
vortex column remains roughly undisturbed.

Table 2 and Fig. 8 show the maximum values of $\varepsilon$. We see that in
all simulations $\varepsilon>1$ is reached and that streaming is triggered.
The baroclinic vortices prove to be an efficient pebble trap, even at extreme
subsolar metallicity. Even an initial dust-to-gas ratio of $10^{-4}$ is
sufficient to trigger the streaming instability inside a vortex.

Fig. 9 shows the azimuthally-averaged pebble column density as a function of
time (bottom panels), as well as the pebble density in the meridional plane of
the maximum pebble density (upper panels). The left hand side panels are for
$\varepsilon=10^{-2}$ and the right hand side panels for
$\varepsilon=10^{-3}$, both for $\mathrm{St}=1$ pebbles. The same graphs but
for $\mathrm{St}=0.05$ pebbles are shown in Fig. 10. The time evolution shows
bands of high pebble concentration that are slowly migrating radially inwards.
The migration speed is equal for both $\mathrm{St}=1$ and $\mathrm{St}=0.05$
pebbles, which evidences that this drift is not the pebble flux. Indeed, this
radial migration corresponds to the migration of the vortices (Paardekooper et
al., 2010). The oscillations in these bands are the epicyclic motion of the
trapped pebbles, induced by the internal vortex motion. We highlight that
these are simulations where the vortex motion was disrupted around the
midplane. The vortex column, above and below the disrupted midplane, keeps the
pebbles trapped.

We measured the intensity of pebble accumulation at different dust-to-gas
ratios. The result is shown in Fig. 11 for the two pebble sizes we consider,
and three different $\varepsilon_{0}$. This quantity, the mass fraction that
participates in the streaming, is also shown in the last column of Table 2 as
$\zeta$. As expected, simulations with higher $\varepsilon_{0}$ also have more
mass in large $\varepsilon$ bins. Promising for planet formation, even
simulations with $\varepsilon_{0}=10^{-4}$ pass the threshold of
$\varepsilon\geq 1$, where the streaming instability sets in. While only 1% of
$\mathrm{St}=0.05$ pebbles accumulate in regions with $\varepsilon\geq 1$ for
$\varepsilon_{0}=10^{-4}$, this is an important step towards higher
concentrations and may eventually lead to gravitational collapse.

###  4.2 Possible Gravitational Collapse

In 2D runs (Raettig et al., 2015) the pebble densities reached were lower than
the Roche density. But in 3D, sedimentation brings the density above this
threshold. We compare the achieved pebble densities with the Roche density to
assess if gravitational collapse would happen had we included self-gravity in
the simulation. We consider the same two disk models as in Raettig et al.
(2015), one with the MMSN with $\Sigma_{g}$(5.2 AU) = 150 g cm-2 (cases 1 and
2) and a more massive disk model with $\Sigma_{g}$(5.2 AU) = 800 g cm-2 (cases
3 and 4). The Roche density at the two radii that we consider is
$\rho_{R}(1AU)=2.83\times{10^{-7}}\ {\rm g\,cm}^{-3}$ (case 1 and 3) and
$\rho_{R}(5.2AU)=2.01\times{10^{-9}}\ {\rm g\,cm}^{-3}$ (case 2 and 4).

In Table 3 we show the maximum dust-to-gas ratios that were reached at the end
of the computational time. The last four columns show the ratio of $\rho_{R}$
and the gas density for the two disk profiles and the two radii. If
$\varepsilon_{\rm max}$ exceeds this ratio, it means that the local dust-to-
gas ratio exceeds $\rho_{R}$ at that radius, and the pebble accumulation would
collapse in the absence of diffusion (indicated by bold face in Table 3).

We see that at 1 AU the maximum dust densities are lower than the Roche
density for all runs except for $\mathrm{St}=1$ pebbles with
$\varepsilon_{0}=0.01$. Yet if we go to 5.2AU there are three more runs for
the more massive disk model that exceed the critical $\rho_{R}$. This means
that the clumps with these high overdensities should collapse and form bodies
that are held together by their own gravity. It is easier to form
planetesimals at larger radii, because $\rho_{R}$ falls off more steeply than
$\rho_{g}$. Therefore, the ratio $\rho_{R}/\rho_{g}$ decreases radially,
making it easier for dust densities to cross the critical density threshold.

Table 3: Dust density compared to Roche density. Run | $\varepsilon_{\rm max}$ | $\rho_{R}/\rho_{g}$ | $\rho_{R}/\rho_{g}$ | $\rho_{R}/\rho_{g}$ | $\rho_{R}/\rho_{g}$
---|---|---|---|---|---
|  | case 1 (5.2AU) | case 2 (1AU) | case 3 (5.2AU) | case 4 (1AU)
3DF05 | 111.2 | 262 | 1607 | 49 | 301
3DF05E-3 | 8.7 | 262 | 1607 | 49 | 301
3DF05E-4 | 2.1 | 262 | 1607 | 49 | 301
3DF1 | 945.4 | 262 | 1607 | 49 | 301
3DF1E-3 | 227.7 | 262 | 1607 | 49 | 301
3DF1E-4 | 34.2 | 262 | 1607 | 49 | 301

###  4.3 Limitations

One limitation of the work is that we include vertical gravity only for the
pebbles, treating the gas as unstratified. We can assess how important this
would be. Since the pebble scale height $H_{d}$ is much smaller than the gas
scale height $H$, the gas density is essentially constant over a distance
$H_{d}$. At one dust scale height, the gas density is 99.0% and 99.95% of the
original density for $\mathrm{St}=0.05$ and $\mathrm{St}=1$, respectively. We
do not expect this variation to be dynamically significant. An important
point, though, is the exclusion of vertical shear, given by (Lyra & Umurhan,
2019)

$\frac{d\varOmega}{dz}=\frac{1}{2}\varOmega h^{2}q_{T}\frac{z}{H^{2}}$ (23)

This vertical shear can in principle be of dynamical consequence (even though
the oscillations leading to vertical shear instability are stabilized by
buoyancy in our setup with finite cooling time). Yet, let us compare this
gravity-induced vertical shear to the shear the pebbles impart on the gas, due
to the backreaction of the drag force. In the absence of gas, the pebbles
would move at the Keplerian rate. Where the pebbles dominate dynamically (at
dust-to-gas ratios above unity), in the midplane, the drag force backreaction
accelerates the gas layer to move at the Keplerian rate. This establishes a
strong localized vertical shear in the gas that, in turn, breaks into Kelvin-
Helmholtz instability (Johansen et al., 2006) generating turbulence and
diffusion. This velocity difference is of magnitude $\eta u_{k}$, so the shear
over a distance $z$ is

$\frac{d\varOmega}{dz}=\frac{1}{2}\varOmega h^{2}\xi{z^{-1}}.$ (24)

So, for $\xi=q_{T}$ (ignoring the density gradient for the argument), the
ratio of the shear given by gravity and the one caused by the pebbles is

${f\equiv}\frac{\left(d\varOmega/dz\right)_{\rm
gravity}}{\left(d\varOmega/dz\right)_{\rm pebbles}}=\frac{z^{2}}{H^{2}}$ (25)

which is vanishingly small in the midplane. For $z=H_{d}$ it yields
$f=H_{d}^{2}/H^{2}=\delta/(\mathrm{St}+\delta)$. For $\delta\approx 10^{-3}$,
this translates into $f\approx 10^{-2}$ for $\mathrm{St}=0.05$ and $f\approx
10^{-3}$ for $\mathrm{St}=1$. That is, the pebbles dominate the vertical shear
as long as $\delta\ll\mathrm{St}$. Fundamentally, the approximation is
justified because the pebble layer is very thin in comparison to the gas
layer. Still, the effect of vertical gravity on the gas, even if small, should
be considered in future work.

##  5 Conclusion

In this paper we have performed three-dimensional simulations of self-
sustained baroclinic vortices, generated by the convective overstability,
along with Lagrangian particles to simulate $\mathrm{St}=0.05$ and
$\mathrm{St}=1$ pebbles in the disk, including the drag force and its
backreaction. These values of $\mathrm{St}$ correspond roughly to cm-sized
pebbles in the inner disk ($\sim$ 5 AU) and mm-sized pebbles in the outer disk
($\sim$ 50 AU) for $\mathrm{St}=0.05$, and m-sized boulders in the inner disk
and cm-sized pebbles in the outer disk for $\mathrm{St}=1.$ While the gas is
unstratified, we included vertical gravity for the pebbles, to allow for
sedimentation. Fu et al. (2014) report, via long evolution time 2D simulations
of Rossby vortices, a decreased vortex lifetime as the dust feedback
increases, either as a result of higher dust-to-gas ratio of larger pebble
sizes. In our 3D simulations, and consistent with these 2D results, for
$\varepsilon_{0}=10^{-3}$ or $\varepsilon_{0}=10^{-4}$ the vortex midplane is
not disrupted, remaining coherent over the simulation period. Also consistent
with the 2D results, for initial dust-to-gas ratio of
$\varepsilon_{0}=10^{-2}$ we see that the vortex column in the midplane is
disrupted after a very short time of pebble concentration. Raettig et al.
(2015) report a similar picture in the two-dimensional simulations for
$\mathrm{St}=1$ pebbles when the local dust-to-gas ratio exceeded a value of
$\approx$ 10\. Due to the sedimentation of dust, this dust-to-gas ratio is
reached after only a few orbital periods for both pebble sizes.

Yet, here we show with 3D models that because the pebble sedimentation occurs
only in a thin layer around the midplane, the vortex column remains coherent
above and below the midplane, and retains its ability to trap pebbles. The
pebbles disrupt the vortex flow within the thickness of the pebble layer, but
they do not destroy the vortex column. The initial concentration in the vortex
grows as the simulations progress. The pebble overdensity drifts not at the
pebble flux rate, but rather in the slower migration timescale of the vortex,
evidencing efficient trapping. Our results are important because based on the
previous 2D result suggesting eventual disruption, the vortex interpretation
of ALMA observations has been called into question. We show instead that the
vortex behaves like a Taylor column, and the pebbles as obstacles to the flow.

For lower dust-to-gas ratio, although the concentration remains below the high
concentrations seen for $\varepsilon_{0}=0.01$, they still reach local
enhancements of $\varepsilon>1$ and, consequently, trigger the streaming
instability. To determine whether the local pebble density enhancements can
cause gravitational collapse were we to turn on pebble self-gravity, we
compare their density with the Roche density. Depending on the disk model and
radial position in the disk, we see that $\mathrm{St}=0.05$ pebbles with a
high initial dust-to-gas ratio, or $\mathrm{St}=1$ pebbles with
$\varepsilon_{0}=10^{-{3}}$ can exceed the Roche density. Had we considered
selfgravity, they would have formed gravitationally bound objects.

We conclude that baroclinic vortices in 3D are an efficient mechanism for
trapping and concentrating pebbles, allowing streaming instability and
planetesimal formation even at subsolar metallicity. In a future work we will
include selfgravity in order to assess the mass distribution formed in these
baroclinic vortices.

We acknowledge conversations with Nienke van der Marel. W. L. acknowledges
support from NSF through grant AST-2007422, and from NASA through grant
20-TCAN20-0011. The simulations presented in this paper utilized the Stampede
cluster of the Texas Advanced Computing Center (TACC) at The University of
Texas at Austin, through XSEDE grant TG-AST140014. This collaboration was made
possible through the support of the Annette Kade Graduate Student Fellowship
Program at the American Museum of Natural History.

## References

  * Adams & Watkins (1995) Adams, F. C., & Watkins, R. 1995, ApJ, 451, 314, doi: 10.1086/176221
  * Barge & Sommeria (1995) Barge, P., & Sommeria, J. 1995, A&A, 295, L1
  * Barranco et al. (2018) Barranco, J. A., Pei, S., & Marcus, P. S. 2018, ApJ, 869, 127, doi: 10.3847/1538-4357/aaec80
  * Bracco et al. (1999) Bracco, A., Chavanis, P. H., Provenzale, A., & Spiegel, E. A. 1999, Physics of Fluids, 11, 2280, doi: 10.1063/1.870090
  * Brandenburg & Dobler (2002) Brandenburg, A., & Dobler, W. 2002, Computer Physics Communications, 147, 471, doi: 10.1016/S0010-4655(02)00334-X
  * Brandenburg & Dobler (2010) —. 2010, Pencil: Finite-difference Code for Compressible Hydrodynamic Flows, Astrophysics Source Code Library.  http://ascl.net/1010.060
  * Brandenburg et al. (2020) Brandenburg, A., Johansen, A., Bourdin, P. A., et al. 2020, arXiv e-prints, arXiv:2009.08231.  https://arxiv.org/abs/2009.08231
  * Casassus et al. (2019) Casassus, S., Marino, S., Lyra, W., et al. 2019, MNRAS, 483, 3278, doi: 10.1093/mnras/sty3269
  * Chavanis (2000) Chavanis, P. H. 2000, A&A, 356, 1089
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237, doi: 10.1006/icar.1995.1058
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, arXiv e-prints, arXiv:2005.11974.  https://arxiv.org/abs/2005.11974
  * Fu et al. (2014) Fu, W., Li, H., Lubow, S., Li, S., & Liang, E. 2014, ApJL, 795, L39, doi: 10.1088/2041-8205/795/2/L39
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56, doi: 10.1051/0004-6361/200912852
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Inaba & Barge (2006) Inaba, S., & Barge, P. 2006, ApJ, 649, 415, doi: 10.1086/506427
  * Johansen et al. (2006) Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219, doi: 10.1086/502968
  * Johansen et al. (2011) Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62, doi: 10.1051/0004-6361/201015979
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022, doi: 10.1038/nature06086
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269, doi: 10.1088/0004-637X/697/2/1269
  * Kerswell (2002) Kerswell, R. R. 2002, Annual Review of Fluid Mechanics, 34, 83, doi: 10.1146/annurev.fluid.34.081701.171829
  * Klahr & Bodenheimer (2006) Klahr, H., & Bodenheimer, P. 2006, ApJ, 639, 432, doi: 10.1086/498928
  * Klahr & Hubbard (2014) Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21, doi: 10.1088/0004-637X/788/1/21
  * Klahr & Henning (1997) Klahr, H. H., & Henning, T. 1997, Icarus, 128, 213, doi: 10.1006/icar.1997.5720
  * Lesur & Papaloizou (2009) Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1, doi: 10.1051/0004-6361/200811577
  * Lesur & Papaloizou (2010) —. 2010, A&A, 513, A60, doi: 10.1051/0004-6361/200913594
  * Lin (2012a) Lin, M.-K. 2012a, ApJ, 754, 21, doi: 10.1088/0004-637X/754/1/21
  * Lin (2012b) —. 2012b, MNRAS, 426, 3211, doi: 10.1111/j.1365-2966.2012.21955.x
  * Lin (2013) —. 2013, ApJ, 765, 84, doi: 10.1088/0004-637X/765/2/84
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805, doi: 10.1086/306900
  * Lyra (2009) Lyra, W. 2009, PhD thesis, Uppsala University
  * Lyra (2013) Lyra, W. 2013, in European Physical Journal Web of Conferences, Vol. 46, European Physical Journal Web of Conferences, 04003, doi: 10.1051/epjconf/20134604003
  * Lyra (2014) —. 2014, ApJ, 789, 77, doi: 10.1088/0004-637X/789/1/77
  * Lyra et al. (2008a) Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008a, A&A, 491, L41, doi: 10.1051/0004-6361:200810626
  * Lyra et al. (2008b) —. 2008b, A&A, 479, 883, doi: 10.1051/0004-6361:20077948
  * Lyra et al. (2009a) —. 2009a, A&A, 493, 1125, doi: 10.1051/0004-6361:200810797
  * Lyra et al. (2009b) Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009b, A&A, 497, 869, doi: 10.1051/0004-6361/200811265
  * Lyra & Klahr (2011) Lyra, W., & Klahr, H. 2011, A&A, 527, A138, doi: 10.1051/0004-6361/201015568
  * Lyra & Lin (2013) Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17, doi: 10.1088/0004-637X/775/1/17
  * Lyra & Mac Low (2012) Lyra, W., & Mac Low, M.-M. 2012, ApJ, 756, 62, doi: 10.1088/0004-637X/756/1/62
  * Lyra et al. (2017) Lyra, W., McNally, C. P., Heinemann, T., & Masset, F. 2017, AJ, 154, 146, doi: 10.3847/1538-3881/aa8811
  * Lyra et al. (2018) Lyra, W., Raettig, N., & Klahr, H. 2018, Research Notes of the American Astronomical Society, 2, 195, doi: 10.3847/2515-5172/aaeac9
  * Lyra et al. (2015) Lyra, W., Turner, N. J., & McNally, C. P. 2015, A&A, 574, A10, doi: 10.1051/0004-6361/201424919
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001, doi: 10.1088/1538-3873/aaf5ff
  * Manger & Klahr (2018) Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125, doi: 10.1093/mnras/sty1909
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841, doi: 10.1093/mnras/staa2943
  * Marcus et al. (2016) Marcus, P. S., Pei, S., Jiang, C.-H., & Barranco, J. A. 2016, ApJ, 833, 148, doi: 10.3847/1538-4357/833/2/148
  * Marcus et al. (2015) Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87, doi: 10.1088/0004-637X/808/1/87
  * McNally et al. (2012) McNally, C. P., Maron, J. L., & Mac Low, M.-M. 2012, ApJS, 200, 7, doi: 10.1088/0067-0049/200/1/7
  * Meheut et al. (2010) Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31, doi: 10.1051/0004-6361/201014000
  * Meheut et al. (2012a) Meheut, H., Keppens, R., Casse, F., & Benz, W. 2012a, A&A, 542, A9, doi: 10.1051/0004-6361/201118500
  * Meheut et al. (2012b) Meheut, H., Meliani, Z., Varniere, P., & Benz, W. 2012b, A&A, 545, A134, doi: 10.1051/0004-6361/201219794
  * Meheut et al. (2012c) Meheut, H., Yu, C., & Lai, D. 2012c, MNRAS, 422, 2399, doi: 10.1111/j.1365-2966.2012.20789.x
  * Miranda et al. (2017) Miranda, R., Li, H., Li, S., & Jin, S. 2017, ApJ, 835, 118, doi: 10.3847/1538-4357/835/2/118
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610, doi: 10.1093/mnras/stt1475
  * Paardekooper et al. (2010) Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. B. 2010, ApJ, 725, 146, doi: 10.1088/0004-637X/725/1/146
  * Pfeil & Klahr (2020) Pfeil, T., & Klahr, H. 2020, arXiv e-prints, arXiv:2008.11195.  https://arxiv.org/abs/2008.11195
  * Raettig et al. (2015) Raettig, N., Klahr, H., & Lyra, W. 2015, ApJ, 804, 35, doi: 10.1088/0004-637X/804/1/35
  * Raettig et al. (2013) Raettig, N., Lyra, W., & Klahr, H. 2013, ApJ, 765, 115, doi: 10.1088/0004-637X/765/2/115
  * Richard et al. (2016) Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571, doi: 10.1093/mnras/stv2898
  * Safronov (1960) Safronov, V. S. 1960, Annales d’Astrophysique, 23, 979
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Sierra et al. (2017) Sierra, A., Lizano, S., & Barge, P. 2017, ApJ, 850, 115, doi: 10.3847/1538-4357/aa94c1
  * Stoll & Kley (2014) Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77, doi: 10.1051/0004-6361/201424114
  * Surville et al. (2016) Surville, C., Mayer, L., & Lin, D. N. C. 2016, ApJ, 831, 82, doi: 10.3847/0004-637X/831/1/82
  * Tanga et al. (1996) Tanga, P., Babiano, A., Dubrulle, B., & Provenzale, A. 1996, Icarus, 121, 158, doi: 10.1006/icar.1996.0076
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217, doi: 10.1086/147861
  * van der Marel et al. (2013) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199, doi: 10.1126/science.1236770
  * van der Marel et al. (2020) van der Marel, N., Birnstiel, T., Garufi, A., et al. 2020, arXiv e-prints, arXiv:2010.10568.  https://arxiv.org/abs/2010.10568
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, APSS, 51, 153, doi: 10.1007/BF00642464
  * Weidenschilling (1980) —. 1980, Icarus, 44, 172, doi: 10.1016/0019-1035(80)90064-0
  * Whipple (1972) Whipple, F. L. 1972, in IAU Symposium, Vol. 45, The Motion, Evolution of Orbits, and Origin of Comets, ed. G. A. Chebotarev, E. I. Kazimirchak-Polonskaia, & B. G. Marsden, 401
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57, doi: 10.1051/0004-6361/200912976


# A Minkowski space embedding to understand Markov models dynamics

David Andrieux

###### Abstract

To describe the behavior of Markov models as parameters are varied, I show how
to embed the space of Markov models within a Minkowski space. This embedding
maintains the inherent distance between different instances of the model. The
coordinates of this embedding emerge from the symmetrized Kullback-Leibler
divergence and are expressed in terms of thermodynamic quantities, organizing
the Minkowski space into equilibrium and nonequilibrium components. With this
approach, we can visualize models’ behaviors and gain a thermodynamic
interpretation of information geometric concepts, even in far-from-equilibrium
scenarios. I illustrate this approach using an analytically solvable molecular
motor model.

## Context and objectives

To gain insight into a system’s dynamical and thermodynamical properties, it
is crucial to understand their variations across the parameter space. For
Markov models, however, the space of parameters has no obvious intrinsic
structure. For example, multiple combinations of parameters can lead to the
same affinities and currents, creating ambiguity and arbitrariness in their
nonequilibrium response [1, 2]. Analysis of parameter sensitivity also shows
that most models are ’sloppy’, in the sense that they have a hierarchy of
eigenvalues spanning multiple scales [3, 4].

One approach to tackle these challenges is through the lens of stochastic
thermodynamics. In stochastic thermodynamics, the (thermo)dynamical properties
of the system are examined as parameters are varied [6, 5]. In the linear
regime, the parameter space has a natural structure based on the Fisher
information metric, which is linked to equilibrium fluctuations [7, 8].
Outside the linear regime, however, the structure of the parameter space and
its connections to thermodynamics has remained largely unexplored. One notable
exception are dynamical equivalence classes, which provide a ’hidden’
structure in the parameter space induced by their transport properties [9].
These equivalence classes uniquely characterize the nonquilibrium conditions
and parameters of the system, providing a well-defined and symmetric
nonequilibrium theory [10].

Another important angle to study the space of parameters is through the lens
of information geometry. In information geometry, the space of parameters
(e.g., for probability distributions or Markov dynamics) is analyzed by
introducing a divergence as a measure of the separation between models [12,
11]. In this context, Teoh et al. demonstrated that the symmetrized Kullback-
Leibler divergence plays a special role in vizualizing probabilistic models
[13, 14]. The symmetrized KL divergence induces a Minkowksi structure in the
space of probabilistic models which preserves the ’distance’ between models.

This paper brings these different concepts and results together. Following the
approach of Teoh et al. [14], I show that the space of Markov models possesses
a Minkoswki structure induced by the symmetrized Kullback-Leibler divergence.
In this Minkowski space, the embedding coordinates disantangle the system’s
equilibrium and nonequilibrium components. The former is expressed in terms of
kinetic activities and forces while the latter coincides with the dynamical
equivalence classes and are expressed in terms of thermodynamic currents and
affinities [6]. These findings provide a unified and systematic way to study
nonequilibrium Markov systems and their dynamical properties.

##  I Thermodynamic decomposition of the space of Markov chains

Let’s consider the space $\Lambda$ of primitive Markov chains on a finite
state space of size $N$. A Markov chain $P\in\Lambda$ admits a unique
stationary state $\pi$ and can be depicted as a graph with $N$ states, $S$
self-transitions, and $E$ edges, each edge representing the allowed transition
between states.

The space $\Lambda$ has $M=2E+S-N$ dimensions while the subspace $\Sigma$ of
equilibrium systems has dimension $E+S-1$. A chain $P$ then has $E-N+1$
independent currents $J_{\alpha}$ and affinities $A_{\alpha}$ [15, 1]. Each
independent current is associated with a fundamental cycle in the graph and a
specific edge $\alpha$ [15, 1].

Multiple representations of the spaces $\Lambda$ and $\Sigma$ exist, e.g.
using exponential families [12, 11, 16] or the rotational representation of
Markov chains [17, 18]. Here I will represent the space $\Lambda$ using a
thermodynamic decomposition [9]. In this decomposition, a chain $P$ is
determined by an equilibrium component $E$ and by it affinities $A_{\alpha}$
as

$\displaystyle P=s\left[E\circ Z(\boldsymbol{A})\right]\,.$ (1)

Here the mapping $s$ transforms a positive matrix into a stochastic one:

$\displaystyle s[K]=\frac{1}{\rho}{\rm diag}(\beta)^{-1}\,K\,{\rm
diag}(\beta)\,,$

where $\rho$ is the largest eigenvalue of $K$ and $\beta$ its largest right
eigenvector. The matrix

$\displaystyle Z_{ij}\equiv\begin{cases}\exp\left(\pm
A_{\alpha}/2\right)&\text{if $i\rightarrow j=\alpha$ in the $\pm$
direction,}\\\ 1&\text{otherwise}\end{cases}\ $

captures the nonequilibrium (asymmetric) components of $P$. The equilibrium
component $E$ of a given $P$ can be recovered as $E=s[(P\circ
P^{T})^{(1/2)}]$, where $P^{(1/2)}$ denotes the Hadamard or elementwise
exponentiation ($E$ is often denoted by $P^{(e)}$ in information geometry).
The full space of Markov models is then obtained by varying the elements
$E\in\Sigma$ and the affinities $A_{\alpha}$.

The dynamics $P$ sharing the same equilibrium component $E$ in the
decomposition (1) are expected to display similar properties [9]. These
dynamical equivalence classes indeed play a special role in nonequilibrium
transport [20]. For example, the equilibrium and nonequilibrium fluctuations
are identical within an equivalence class [19]. In addition, the currents
fluctuations can be obtained from the current response curve itself and all
response coefficients are fully symmetric and are expressed in terms of
equilibrium coefficients [10]. The equivalence classes also capture the
minimal distance, in the information-theoric sense, between models under fixed
nonequilibrium conditions [2].

The decomposition (1) offers a clear thermodynamic interpretation and provides
important features when studying transport properties [19, 10]. However, it
still lacks a natural geometric structure or notion of distance. In the next
section, I show how the space $\Lambda$ can be embedded in a Minkowski space
with a clear physical interpretation.

##  II Intensive embedding in Minkowski space

Teoh et al. showed that the symmetrized Kullback-Leibler divergence $D_{sKL}$
offers a natural distance in the space of exponential families of probability
distributions [14]. Using that distance [21], an exponential family of
probability distributions $f$ with $N$ parameters can be embedded
isometrically in a $N+N$-dimensional Minkowski space. That is, $D_{sKL}$ can
be decomposed into $N+N$ components $S^{\pm}_{i}(f)$ such that

$\displaystyle D_{sKL}(f,f^{\prime})=\sum_{i=1}^{N}\left[\left(\Delta
S^{+}_{i}\right)^{2}-\left(\Delta S^{-}_{i}\right)^{2}\right]\,,$

with $\Delta S_{i}^{\pm}=S^{\pm}_{i}(f)-S^{\pm}_{i}(f^{\prime})$. Using these
$N+N$ coordinates, the models can be visualized as control parameters are
varied, preserving the distance between probability distributions. All
pairwise distances between model instances are given by the symmetrized
Kullback-Leibler divergence [14].

In this section I adapt this approach to the context of Markov chains and show
that the resulting coordinates have a clear physical interpretation in terms
of kinetic and thermodynamical quantities.

The symmetrized Kullback-Leibler divergence between two Markov chains reads

$\displaystyle D_{sKL}(P||G)$ $\displaystyle=$ $\displaystyle
D_{KL}(P||G)+D_{KL}(G||P)$ (2) $\displaystyle=$
$\displaystyle\sum_{ij}\left[\pi_{i}P_{ij}-\nu_{i}G_{ij}\right]\ln\frac{P_{ij}}{G_{ij}}\geq
0\,.\quad$

Here $D_{KL}(P||G)=\sum_{ij}\pi_{i}P_{ij}\ln P_{ij}/G_{ij}$ is the Kullback-
Leibler divergence and $\pi$ and $\nu$ are the stationary distributions of $P$
and $G$, respectively.

To calculate and recast the divergence (2) in a more useful form, I now use
the thermodynamic decomposition $P=s[E\circ Z]$ and $G=s[E^{\prime}\circ
Z^{\prime}]$. As demonstated in appendix A, this leads to

$\displaystyle D_{sKL}(P||G)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\alpha}(A_{\alpha}-A_{\alpha}^{\prime})(J_{\alpha}-J_{\alpha}^{\prime})$
(3) $\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{e}(X_{e}-X_{e}^{\prime})(Y_{e}-Y_{e}^{\prime})\,.$

Here $J_{\alpha}$ is the average current with corresponding affinity
$A_{\alpha}$ while $Y_{e}=\pi_{i}P_{ij}+\pi_{j}P_{ji}$ is the activity along
the edge $e=i\longleftrightarrow j$ and $X_{e}=\ln P_{ij}P_{ji}$. By analogy
with the affinities, the quantities $X_{e}$ play the role of ’kinetic forces’
driving the activity $Y_{e}$ along the edge [22].

Not all variables $Y_{e}$, however, are independent since $\sum_{e}Y_{e}=1$.
One variable $r$ can thus be eliminated as $Y_{r}=1-\sum_{e\neq r}Y_{e}$. With
this substitution, Eq. (3) retains the same form but with effective forces

$\displaystyle\bar{X}_{e}=X_{e}-X_{r}\,.$

For notational simplicity, in the rest of the paper I assume that one of the
dependent edge $e$ has been eliminated and write $X_{e}=\bar{X}_{e}$.

The divergence (3) between two models $P$ and $G$ comprises two types of
contributions. The first contribution captures the nonequilibrium (or
antisymmetric) aspects and is expressed in terms of the thermodynamic current
and affinities, while the second contribution reflects the changes in the
equilibrium (or symmetric) part of the dynamics. As I will argue below, the
divergence (3) measures both the time and ’spatial’ asymmetry between two
models.

Let’s now express the divergence (3) in terms of space-time coordinates. As
shown in appendix A, the coordinates

$\displaystyle S_{\alpha}^{\pm}=\frac{1}{4}(A_{\alpha}\pm J_{\alpha})\quad{\rm
and}\quad S_{e}^{\pm}=\frac{1}{4}(X_{e}\pm Y_{e})$ (4)

provide an intensive embedding in a $M+M$ Minkowski space:

$\displaystyle D_{sKL}(P||G)=\sum_{l=1}^{M}\left[(\Delta
S^{+}_{l})^{2}-(\Delta S^{-}_{l})^{2}\right]\geq 0\,,$ (5)

where the sum runs over $l=(\alpha,e)$ with $e\neq r$.

The Minkowski space is $M+M$ dimensional, where $M=2E+S-N$ corresponds to the
number of independent parameters of the Markov models. The spacelike
coordinates $S^{+}$ are greater than the timelike coordinates $S^{-}$ for each
variable,
$[S^{+}_{l}(P)-S^{+}_{l}(G)]^{2}\geq[S^{-}_{l}(P)-S^{-}_{l}(G)]^{2}$. It is
also possible to center and rescale the coordinates by factors $\lambda_{l}$,
e.g. to minimize the sum of the square of the coordinates or adjust the units
[14].

Using coordinates (4), the set of Markov chains can thus be vizualized as a
manifold as control parameters are varied, preserving the symmetric divergence
between systems and with clear thermodynamic interpretation of the
coordinates. This constitues our main result.

##  III Kinetic and thermodynamics properties of the model manifold

I now explore the kinetic and thermodynamic properties of the embedding (5)
before analyzing an example of molecular motor in Section IV. Taken together,
the results below show that the symmetrized KL divergence measures both the
time asymmetry (through the currents and affinities) and the spatial asymmetry
(through the activities and kinetic forces) of the dynamics.

###  III.1 Equilibrium and iso-affinity manifolds

Let’s first consider the case where the affinities between models are
identical, $A_{\alpha}=A^{\prime}_{\alpha}$. In these sub-manifolds, the
separations $\Delta S_{\alpha}^{\pm}$ vanish and the divergence (3) simplifies
to

$\displaystyle D_{sKL}(P||G)$ $\displaystyle=$
$\displaystyle\sum_{e}\left[(\Delta S^{+}_{e})^{2}-(\Delta
S^{-}_{e})^{2}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{e}(X_{e}-X_{e}^{\prime})(Y_{e}-Y_{e}^{\prime})\,.$

The embedding coordinates and the separation between two models are directly
expressed in terms of the systems’ activities and kinetic forces $X_{e}$ and
$X^{\prime}_{e}$, which in turn measure the inhomogeneity (spatial asymmetry)
of the dynamics across the graph.

Note, however, that the coordinates $S^{\pm}_{\alpha}$ still vary along the
iso-affinity manifold even if $\Delta S^{\pm}_{\alpha}=0$. Only when
$A_{\alpha}=0$ are the coordinates $S^{\pm}_{\alpha}$ independent of
$S^{\pm}_{e}$ (more precisely, $S^{\pm}_{\alpha}=0$ in that case). The
equilibrium manifold is thus special in the sense that it is entirely
described by the coordinates $S^{\pm}_{e}$.

###  III.2 Dynamical equivalence classes

Dynamical equivalence classes play a special role in nonequilibrium response
theory [10]. Remarkably, they emerge naturally in this framework.

A dynamical equivalence class is obtained when all the model instances share
the same equilibrium component $E$, i.e. when $X_{e}=X_{e}^{\prime}$. In that
case the separations $\Delta S^{\pm}_{e}$ vanish and the symmetrized KL
divergence (3) simplifies to

$\displaystyle D_{sKL}(P||G)$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\left[(\Delta S^{+}_{\alpha})^{2}-(\Delta
S^{-}_{\alpha})^{2}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\alpha}(A_{\alpha}-A_{\alpha}^{\prime})(J_{\alpha}-J_{\alpha}^{\prime})\,.$

The embedding coordinates and the separation between two models in an
equivalence class are directly expressed in terms of the system’s thermodyamic
currents and affinities, which in turn measure the time asymmetry of the
dynamics.

As in the case of the iso-affinity manifolds, the coordinates $S^{\pm}_{e}$
vary along the dynamical equivalence class even if $\Delta S^{\pm}_{e}=0$. In
contrast to the equilibrium case however, even when $X_{e}=0$ the coordinates
$S^{\pm}_{e}=\pm Y_{e}$ still depend on $S^{\pm}_{\alpha}$. Nonetheless, the
class $X_{e}=0$ can act as a distinct reference dynamics: it is obtained from
$E=s[T]$, where $T$ is the adjacency matrix of $P$, and thus corresponds to
Markov models derived from a ’spatially homogeneous’ equilibrium component
(this will be illustrated in the molecular motor example in section IV).

###  III.3 Entropy production

The dissipation or entropy production
$\Delta_{i}S=\sum_{\alpha}J_{\alpha}A_{\alpha}\geq 0$ plays an important role
in thermodynamics. In the present framework, the entropy production emerges as
the symmetrized KL divergence with respect to the equilibrium component
$E=P^{(e)}$ [23]:

$\displaystyle\Delta_{i}S=\sum_{\alpha}J_{\alpha}\times
A_{\alpha}=2\,D_{sKL}\left(P||E\right)$

since $J_{\alpha}^{\prime}=A_{\alpha}^{\prime}=0$ and $X_{e}^{\prime}=X_{e}$
as $E$ belong to the equivalence class of $P$.

The entropy production can also be obtained as $D_{sKL}(P||P^{(m)})$ where
$P^{(m)}=(P+P^{*})/2$ with $P^{*}$ the time-reversed dynamics [23]. This
occurs because $J_{\alpha}^{\prime}=A_{\alpha}^{\prime}=0$ and
$Y_{e}^{\prime}=Y_{e}$ in that case.

Finally, the entropy production can also be expressed as the divergence
between a dynamics and its time-reserval [24, 23]:
$\Delta_{i}S=2D_{sKL}(P||P^{*})$, as seen from Eq. (3) and using that the
currents and affinities of $P^{*}$ take opposite values compared to $P$ while
having the same equilibrium components [25].

##  IV A molecular motor example

Consider a Markov chain representing a molecular motor with $2\ell$ states
corresponding to different conformations of the protein complex. These states
form a cycle of periodicity $2\ell$ corresponding to a revolution by 360° for
a rotary motor or a reinitialization step for a linear motor. The motor
alternates between two types of states according to the transition matrix [26]

$\displaystyle P=\begin{pmatrix}0&p_{1}&&&&1-p_{1}\\\ 1-p_{2}&&p_{2}&&&\\\
&1-p_{1}&&p_{1}&&\\\ &&\ddots&0&\ddots&\\\ &&&1-p_{1}&0&p_{1}\\\
p_{2}&&&&1-p_{2}&0\end{pmatrix}_{2\ell\times 2\ell}.$

The matrix $P$ is doubly stochastic, so that its stationary state is given by
the uniform distribution ${\bf\pi}=(1,1,\cdots,1)/2\ell$ for all parameters
$(p_{1},p_{2})$. The average current $J$ and affinity $A$ take the form

$\displaystyle J=\frac{1}{2\ell}\left(p_{1}+p_{2}-1\right)$ (6)

and

$\displaystyle A=\ell\,\ln\frac{p_{1}p_{2}}{(1-p_{1})(1-p_{2})}\,.$ (7)

The system is at equilibrium when $J=A=0$ (i.e., when $p_{2}=1-p_{1}$).

The activities $Y_{e}=J_{e}^{+}+J_{e}^{-}$ along the edges read
$Y_{1}=(p_{1}-p_{2}+1)/2\ell$ and $Y_{2}=(p_{2}-p_{1}+1)/2\ell$ depending if
the transition is odd or even, respectively. However, only one edge is
independent due to the symmetry of the model and the constraint
$Y_{1}+Y_{2}=1/\ell$. Taking

$\displaystyle Y=Y_{1}=\frac{1}{2\ell}(p_{1}-p_{2}+1)$ (8)

as the independent activity leads to an effective force

$\displaystyle
X=\ell\,(X_{1}-X_{2})=\ell\,\ln\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})}\,.$ (9)

Using expressions (7) and (9), the model parameters can be expressed as

$\displaystyle p_{1}=\frac{{\rm e}^{(A+X)/2\ell}}{1+{\rm
e}^{(A+X)/2\ell}}\,,\quad\quad p_{2}=\frac{{\rm e}^{(A-X)/2\ell}}{1+{\rm
e}^{(A-X)/2\ell}}\,.$ (10)

The iso-affinity classes are obtained by fixing $A$ and varying $X$ while the
dynamical equivalence classes are obtained by varying the affinity $A$ at
fixed value of $X$ [27]. The equilibrium manifold $A=0$ is given by the line
$p_{2}=1-p_{1}$ (i.e., reversible motors) and the dynamical class $X=0$ by the
line $p_{1}=p_{2}$ in parameter space (i.e., spatially homogeneous motors).

Figure 1:  Current and activity response curves for the molecular motor along
equivalence classes. (a) Current response curves for different dynamical
equivalence classes $X$. Each response curve is anti-symmetric and invariant
under the transformation $X\rightarrow-X$ (not shown). (b) Activity response
curves for different equivalence classes $A$. Each response curve is anti-
symmetric (after translation) and invariant under the transformation
$A\rightarrow-A$. In this system, the current and activity response curves
have the same functional form within their respective (dynamical and iso-
affinity) equivalence classes.

Within these manifolds or equivalence classes, the response curves are well-
defined (Fig. 1). In this system, the current and activity have the same
functional form within their respective equivalence classes. The current
response curves are anti-symmetric, satisfy a fluctuation theorem, and capture
fluctuation at all orders (i.e., deriving $J_{X}(A)$ with respect to $A$
generates the various cumulants for any value of the affinity) [27].

The activity response curves are antisymmetric (after translation to the
origin) in the iso-affinity classes and are invariant under the transformation
$A\rightarrow-A$ due to the time-reversal symmetry. As demonstrated in
Appendix B, the activity satisfies a fluctuation theorem and, similar to the
current, the response curves encode information about its fluctuations [10,
27]. The generalization of these observations to an arbitrary Markov system
will be explored elsewhere.

The symmatrized KL divergence between two models $P=(p_{1},p_{2})$ and
$G=(g_{1},g_{2})$ reads

$\displaystyle 2\,D_{sKL}$ $\displaystyle=$
$\displaystyle(p_{1}-g_{1})\ln\frac{p_{1}}{g_{1}}+(p_{2}-g_{2})\ln\frac{p_{2}}{g_{2}}$
$\displaystyle+$
$\displaystyle(g_{1}-p_{1})\ln\frac{1-p_{1}}{1-g_{1}}+(g_{2}-p_{2})\ln\frac{1-p_{2}}{1-g_{2}}\,.$

Inserting the parametrization (10) into this expression leads to

$\displaystyle D_{sKL}(P||G)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(A_{P}-A_{G})(J_{P}-J_{G})$ $\displaystyle+$
$\displaystyle\frac{1}{2}(X_{P}-X_{G})(Y_{P}-Y_{G})\,.\quad\quad$

As expected from our general result, the divergence is expressed in terms of
the system’s current and affinity on the one hand, and the activity and
kinetic force on the other hand.

As a result, the four embedding coordinates take the form

$\displaystyle S_{\alpha}^{\pm}=\frac{1}{4}(A\pm J)\,,\quad\quad S_{e}^{\pm}$
$\displaystyle=$ $\displaystyle\frac{1}{4}(X\pm Y)\,.$

In these coordinates, the molecular motor models form a two-dimensional
manifold (Fig. 2) satisfying the Minkowski structure (5).

The embedding coordinates have similar shapes due to the similarity of the
functional forms of $J_{X}(A)$ and $Y_{A}(X)$ (Fig. 1). Note that the number
of embedding coordinates is independent of the size $\ell$ of the motor.

Figure 2: Model manifold for the molecular motor. The molecular model manifold
is embedded into $2+2$ dimensions. (a) A three-dimensional projection of the
manifold, which is colored based on the current $J$. (b) The different two-
dimensional projections of the manifold. The spacelike coordinates are denoted
as $S^{+}$ and the timelike by $S^{-}$. For visual clarity, the embedding
coordinates are scaled as $S_{\alpha}^{\pm}=(A/\lambda\pm\lambda J)/4$ and
$S_{e}^{\pm}=(X/\lambda\pm\lambda Y)/4$ with $\lambda=L$. All 2-d projections
have the same scale.

The entropy production is obtained as the divergence from the equilibrium
component $E=P(X,0)=P^{(e)}$ [23]:

$\displaystyle\Delta_{i}S=J\times A=2\,D_{sKL}\left(P||E\right)\,.$

Similarly, the entropy production is also measured as the divergence with
respect to $P^{(m)}=(P+P^{*})/2=(1/2)P(p_{1}-p_{2}+1,p_{2}-p_{1}+1)$ [23]. In
that case $Y_{P}=Y_{P^{(m)}}$ even if $X_{P}\neq X_{P^{(m)}}$, so that
$\Delta_{i}S=2D_{sKL}(P||P^{(m)})$.

##  V Summary and outlook

This work brings together thermodynamic and information geometric concepts,
showing that thermodynamic quantities emerge naturally from the Minkoswki
structure induced by the symmetrized KL divergence. While other divergences
could be used, the symmetrized KL divergence bypasses the curse of
dimensionality and provides the lowest embedding dimension, at least for
exponential families of probability distributions [14, 28].

The embedding coordinates structure the space of Markov models in terms of
equilibrium and nonequilibrium (i.e. symmetric and antisymmetric) components.
The equilibrium coordinates are expressed in terms of the system’s activities
and kinetic forces while the nonequilibrium coordinates coincide with the
dynamical equivalence classes and are expressed in terms of the system’s
currents and affinities. The symmetrized Kullback-Leibler divergence thus
measures the difference in temporal and spatial structures between models.

This framework provides a systematic way to analyze the space of Markov models
and to gain insights into the behavior and design of mesocopic devices. For
example, the Minkowski embedding allows to explore a system’s response in the
directions ’orthogonal’ to the dynamical equivalence classes, where additional
symmetries are expected to appear. The embedding can be used to determine the
minimum parameters needed to describe model data or to identify the governing
parameters for model prediction and coarse graining [4, 28]. The space-time
structure could also help untangle important phenomena such as dynamical or
nonequilibrium phase transitions.

## Appendix A: Derivation of the embedding coordinates

The first step to obtain Eq. (3) is to calculate the log ratios $\ln
P_{ij}/G_{ij}$ that appear in Eq. (2). Using the thermodynamic decomposition
(1) with $P=s[E\circ Z]$ and $G=s[E^{\prime}\circ Z^{\prime}]$, the log ratios
take the form

$\displaystyle\ln\frac{P_{ij}}{G_{ij}}=\ln\frac{\rho^{\prime}}{\rho}+\ln\left(\frac{\alpha_{j}^{\prime}}{\alpha_{i}^{\prime}}\frac{\alpha_{i}}{\alpha_{j}}\right)+\ln\frac{Z_{ij}}{Z^{\prime}_{ij}}+\ln\frac{E_{ij}}{E^{\prime}_{ij}}\,.$
(11)

The first term $\ln\rho^{\prime}/\rho$ vanishes when inserted into expression
(2) since $\sum\pi_{i}P_{ij}=\sum\nu_{i}G_{ij}=1$. The second term vanishes as
well when averaged over a stochastic dynamics (see for example Lemma 4.3 (iii)
in ref. [16]). The third term in Eq. (11) reads

$\displaystyle(1/2)\sum_{\alpha}(A_{\alpha}-A_{\alpha}^{\prime})(J_{\alpha}-J_{\alpha}^{\prime})$
(12)

when inserted into (2), where I used that $\ln Z_{ij}=\pm A_{\alpha}/2$ if the
transition $i\rightarrow j$ corresponds to $\alpha$ in the positive (negative)
direction, and $0$ otherwise.

To calculate the last term $\ln E_{ij}/E^{\prime}_{ij}$, remember that
$E_{ij}=(1/\lambda)(\gamma_{i}/\gamma_{j})\sqrt{P_{ij}P_{ji}}$ where $\lambda$
and $\gamma$ are the Perron eigenvalue and eigenvector of
$\sqrt{P_{ij}P_{ji}}$, and similarly for $E^{\prime}$. Here also the terms
$\ln\lambda^{\prime}/\lambda$ and
$\ln(\gamma^{\prime}_{j}\gamma_{i}/\gamma^{\prime}_{i}\gamma_{j})$ will vanish
when averaged over $\sum\pi_{i}P_{ij}-\nu_{i}G_{ij}$. The remaining factor
then becomes

$\displaystyle(1/2)\sum_{e}(X_{e}-X^{\prime}_{e})\left(Y_{e}-Y^{\prime}_{e}\right)\,,$
(13)

where the sum runs over the $E+S$ edges and self-transitions
$e=i\longleftrightarrow j$, $Y_{e}=J_{e}^{+}+J_{e}^{-}$ is the average
activity (sum of forward and backward currents) along an edge, and $X_{e}=\ln
P_{ij}P_{ji}$. Note that only $E+S-1$ factors $Y_{e}$ are independent since
$\sum_{e}Y_{e}=1$. Summing the different terms gives the decomposition (3).

Finally, to obtain the embedding coordinates (4), each term in expression (12)
can be written as [14]

$\displaystyle(1/2)(A_{\alpha}-A_{\alpha}^{\prime})(J_{\alpha}-J_{\alpha}^{\prime})=(\Delta
S^{+}_{\alpha})^{2}-(\Delta S^{-}_{\alpha})^{2}$

with the two coordinates

$\displaystyle S^{\pm}_{\alpha}=(1/4)(A_{\alpha}\pm J_{\alpha})\,.$

The same factorization applies to the terms
$(X_{e}-X^{\prime}_{e})(Y_{e}-Y^{\prime}_{e})$ in Eq. (13). $\square$

## Appendix B: Fluctuation theorem for the molecular motor activity

In this section I show that the activity (8) of the molecular motor model
satisfies a fluctuation theorem similar to the one for the currents (analytic
expressions for the current fluctuations of the molecular motor are given in
Ref. [27]).

In a given realization of the Markov chain, the activity will fluctuate as
$y(n)=1$ if the transition $i_{n-1}\rightarrow i_{n}$ matches one of the edges
$i\longleftrightarrow i+1$ in any direction and with $i$ odd (here denoted as
an ’odd edge’), and $0$ otherwise. The fluctuations of the activity are
captured by the generating function

$\displaystyle
q(\lambda)=\lim_{n\rightarrow\infty}-\frac{1}{n}\left\langle{\rm
e}^{-\lambda\sum_{n}y(n)}\right\rangle_{P}\,.$

The mean activity and its fluctuations are obtained by successive derivations
of $q$ with respect to $\lambda$. The generating function $q$ is in turn
obtained as

$\displaystyle q(\lambda)=-\ln\rho[P\circ G(\lambda)]\,,$

where $\rho$ denotes the largest eigenvalue of $P\circ G$ and

$\displaystyle
G_{ij}(\lambda)\equiv\begin{cases}\exp\left(-\lambda/\ell\right)&\text{if
$i\rightarrow j$ is an odd edge,}\\\ 1&\text{otherwise .}\end{cases}$

However, the activity $y$ and its generating function $q(\lambda)$ do not
display any symmetry. Rather, the centered activity

$\displaystyle\bar{y}=y-1/2\ell$

such that $\left\langle\bar{y}\right\rangle(X=0)=0$ satisfies a fluctuation
theorem.

The fluctuation theorem for the centered activity $\bar{y}$ can be
demonstrated by calculating its generating function explicitly:

$\displaystyle\bar{q}(\lambda)=-\frac{1}{2}\ln\Big{[}\left(p_{1}+(1-p_{1})\text{e}^{\lambda/\ell}\right)\left(p_{2}+(1-p_{2})\text{e}^{-\lambda/\ell}\right)\Big{]}$

which satisfies the fluctuation theorem

$\displaystyle\bar{q}(\lambda)=\bar{q}(X-\lambda)\,.$

The kinetic force $X$ plays the same role as the affinity $A$ in the current
fluctuation theorem. Note that generating function of $\bar{y}$ and $y$ are
related as $\bar{q}(\lambda)=q(\lambda)-\lambda/2\ell$.

Similar to the current fluctuations within dynamical equivalence classes, the
activity fluctuations within iso-affinity classes display additional
symmetries beyond the fluctuation theorem. Expressing $\bar{q}$ in terms of
the kinetic force and affinity of $P$,

$\displaystyle\bar{q}(\lambda)=-\frac{1}{2}\ln\frac{\left(\text{e}^{(A-X)/2\ell}+\text{e}^{-\lambda/\ell}\right)\left(\text{e}^{(A+X)/2\ell}+\text{e}^{\lambda/\ell}\right)}{\left(\text{e}^{(A-X)/2\ell}+1\right)\left(\text{e}^{(A+X)2\ell}+1\right)},$

shows that two dynamics $X$ and $X^{\prime}$ within an iso-affinity class
share the same generating functions modulo the shifts
$\lambda\rightarrow\lambda+(X^{\prime}-X)/2$ and
$q^{\prime}=q+(\bar{q}^{\prime}(0)-\bar{q}(0))$. This shared functional form
allows to obtain the activity fluctuations at all orders and for any value of
the affinity $A$ by deriving the activity response curve [10].

The activity fluctuation theorem and the role of iso-affinity classes for a
general system will be explored elsewhere.

## References

  * [1] T. L. Hill, Free Energy Transduction and Biochemical Cycle Kinetics (Dover, 2005).
  * [2] D. Andrieux, Revealing hidden structures and symmetries in nonequilibrium transport, arXiv:2401.14496 (2024).
  * [3] RN Gutenkun, JJ Waterfall, FP Casey, et al., Universally sloppy parameter sensitivities in systems biology models, PLoS Comput Biol 3: e189 (2007).
  * [4] B. B. Machta, R. Chachra, M. K. Transtrum, J. P. Sethna, Parameter Space Compression Underlies Emergent Theories and Predictive Models, Science 342, 604 (2013).
  * [5] G. Falasco and M. Esposito, Macroscopic Stochastic Thermodynamics, arXiv:2307.12406 (2023).
  * [6] T. Aslyamov and M. Esposito, Nonequilibrium Response for Markov Jump Processes: Exact Results and Tight Bounds, Phys. Rev. Lett. 132, 037101 (2024).
  * [7] D. Sivak and G. E. Crooks, Thermodynamic metrics and optimal paths, Phys. Rev. Lett. 108, 190602 (2012).
  * [8] P. R. Zulkowski, D. A. Sivak, G. E. Crooks, and M. R. DeWeese, The geometry of thermodynamic control, Phys. Rev. E. 86, 0141148 (2012).
  * [9] D. Andrieux, Equivalence classes for large deviations, arXiv:1208.5699 (2012).
  * [10] D. Andrieux, Fully symmetric nonequilibrium response of stochastic systems, arXiv:2205.10784 (2022).
  * [11] Shun-ichi Amari, Information Geometry and Its Applications (Springer, 2016).
  * [12] H. Nagaoka, The exponential family of Markov chains and its information geometry, The proceedings of the Symposium on Information Theory and Its Applications 28, 601 (2005). Also available on arXiv:1701.06119.
  * [13] K. Quinn et al., Visualizing probabilistic models and data with Intensive Principal Component Analysis, Proc. Nat. Acad. Sci. 116, 13762 (2019).
  * [14] H. K. Teoh et al., Visualizing probabilistic models in Minkowski space with intensive symmetrized Kullback-Leibler embedding, Phys. Rev. Res. 2, 033221 (2020).
  * [15] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Rev. Mod. Phys 48, 571 (1976).
  * [16] G. Wolfer and S. Watanabe, Information Geometry of Reversible Markov Chains, Information Geometry 4, 393 (2021).
  * [17] S. Alpern, Rotational representations of stochastic matrices, The Annals of Probability 11, 789 (1983).
  * [18] S. Kalpazidou, Cycle Representations of Markov Processes (Springer, 2006).
  * [19] D. Andrieux, Nonequilibrium large deviations are determined by equilibrium dynamics, arXiv:1212.1807 (2012).
  * [20] Two dynamics $P$ and $Q$ belong to the same class if there exists a scalar $\rho\geq 0$ such that: $(P\circ P^{T})^{(1/2)}=\rho\,\,(Q\circ Q^{T})^{(1/2)}$.
  * [21] While $D_{sKL}(P||G)$ is positive, symmetric and vanishes when $P=G$, in general it doesn’t define a distance since the triangle inequality is not always respected.
  * [22] The role of the symmetric component in nonequilibrium response was highlighted by C. Maes and colleagues, see e.g. C. Maes and M. H. van Wieren, Time-symmetric fluctuations in nonequilibrium systems, Phys. Rev. Lett. 96, 240601 (2006).
  * [23] D. Andrieux, Irreversibility as divergence from equilibrium, arxiv:2404.01978 (2024).
  * [24] P. Gaspard, Time-reversed dynamical entropy and irreversibility in Markovian random processe, J. Stat. Phys. 117, 599 (2004).
  * [25] In analogy with the entropy production, the quantity $D_{sKL}(P||R)=(1/2)\sum_{e}X_{e}(Y_{e}-Y^{\prime}_{e})\geq 0$, where $R=s[T\circ Z(\boldsymbol{A})]$ with $\boldsymbol{A}$ the affinities of $P$ and $T$ its adjacency matrix, can be interpreted as measuring the asymmetry or ’inhomogeneity’ in the graph’s dynamics. Its properties remain to be explored.
  * [26] D. Andrieux and P. Gaspard, Fluctuation theorems and the nonequilibrium thermodynamics of molecular motors, Phys. Rev. E 74, 011906 (2006).
  * [27] D. Andrieux, Making sense of nonequilibrium current fluctuations: A molecular motor example, arxiv:2306.01445 (2023).
  * [28] K. N. Quinn, M. C. Abbott, M. K. Transtrum, B. B. Machta and J. P. Sethna, Information geometry for multiparameter models: new perspectives on the origin of simplicity, Rep. Prog. Phys. 86, 035901 (2023).