# Effects of different discretisations of the Laplacian upon stochastic
simulations of reaction-diffusion systems
on both static and growing domains

Bartosz J. Bartmanski† and Ruth E. Baker† †Mathematical Institute, University
of Oxford,
Woodstock Road, Oxford, OX2 6GG, UK

###### Abstract

By discretising space into compartments and letting system dynamics be
governed by the reaction-diffusion master equation, it is possible to derive
and simulate a stochastic model of reaction and diffusion on an arbitrary
domain. However, there are many implementation choices involved in this
process, such as the choice of discretisation and method of derivation of the
diffusive jump rates, and it is not clear a priori how these affect model
predictions. To shed light on this issue, in this work we explore how a
variety of discretisations and method for derivation of the diffusive jump
rates affect the outputs of stochastic simulations of reaction-diffusion
models, in particular using Turing’s model of pattern formation as a key
example. We consider both static and uniformly growing domains and demonstrate
that, while only minor differences are observed for simple reaction-diffusion
systems, there can be vast differences in model predictions for systems that
include complicated reaction kinetics, such as Turing’s model of pattern
formation. Our work highlights that care must be taken in using the reaction-
diffusion master equation to make predictions as to the dynamics of stochastic
reaction-diffusion systems.

###### keywords:

stochastic simulation , diffusion , reaction , pattern formation.

###### MSC:

[2010] 35K05 , 68U20 , 65C35 92-08.

##  1 Introduction

Many biological processes hinge on the reaction and diffusion of constituent
components and most are, in addition, subject to the inherent stochasticity of
the natural world. Due to the nonlinearity of interacting components,
mathematical models provide a crucial means to interrogate how the combination
of reaction and diffusion drive a given process [4, 11, 14, 15, 19, 31].

Models of reaction-diffusion systems can be roughly divided into macroscopic,
mesoscopic and microscopic [10]. On the macroscale, partial differential
equations (PDEs) are often used to model reaction-diffusion systems in terms
of the evolution of species concentration in time and space. Such models can
be analysed using a range of existing mathematical tools, however they do not
take into account stochastic effects. At the other end of the spectrum we have
microscopic descriptions such as molecular dynamics models, where each
individual molecule is described by its position and velocity, and molecules
react when within a certain radius of each other [1]. Although they provide a
very detailed account of the reaction-diffusion process under consideration,
they incur extremely high computational cost and are, as such, unfeasible for
many real-world models [21, 29]. A mesoscopic description aims to strike a
balance between the two, including descriptions of individual molecule
dynamics without incurring prohibitive computational overheads. An example
mesoscopic description is that of the reaction-diffusion master equation
framework (RDME) [12, 13, 28]; in this model, the domain is discretised into
compartments, molecules can react with other molecules in the same
compartment, and diffusion is represented by jumps between neighbouring
compartments. The RDME framework has been used to explore a range of
biological phenomena, and it is this description of reaction and diffusion on
which we focus in this work.

To model diffusion, the RDME framework requires computation of the rate per
unit time of a molecule to jump from one compartment to a neighbouring one.
Previous research in this direction includes the work by Engblom et al. [5],
where the finite element method is used to calculate the diffusive jump rates,
and that of Hellander and Petzold [10] who developed a description of
diffusion for a more general compartment shapes. In addition, Meinecke and
Lötstedt explored the impact of different methods of diffusive jump rate
derivation upon simulation of diffusive systems [16, 18, 17]. To mitigate
issues with convergence of the RDME in the context of higher order reactions,
Isaacson developed a convergent RDME that allows molecules to react even when
they are not in the same compartment [13]. However, this approach it does not
consider how differences in the derivation of the diffusive jump rates affect
model predictions.

###  1.1 Aims and outline

Our aim in this work is to understand how the different methods of derivation
of diffusive jump rates, and different domain discretisations, affect model
predictions for reaction-diffusion systems, including those with bimolecular
reactions and reaction kinetics capable of pattern formation. We will use
Meinecke and Lötstedt’s [18] derivations of diffusive jump rates as a basis
for our work. These derivations are based on three different numerical
schemes: the finite volume method (FVM); the finite element method (FEM); and
the finite difference method (FDM). In addition we include the fourth approach
taken by Meinecke and Lötstedt that is based on first exit times (FETs) [18].
In Section 2 we outline the results of Meinecke and Lötstedt for a static
domain [18], extend them to a uniformly growing domain, and present the
stochastic simulation algorithm that we will use throughout this work. Our
results are outlined in Section 3, including the effects of different
diffusive jump rates on pattern formation in Section 3.1.4. A brief discussion
of the results outlined in this paper can be found Section 4.

##  2 The reaction-diffusion master equation framework

The chemical master equation (CME) framework describes the reaction of
chemical species under the assumption that the molecules are well-mixed [24].
In general, the CME is intractable and, as a result, Gillespie [8] developed a
stochastic simulation algorithm (SSA) to generate sample paths of the system
under consideration. The RDME framework is an extension of the CME where the
spatial domain is discretised into compartments, each of which is assumed
well-mixed, and diffusion is modelled using a set of first-order reactions
that describe the jumps of molecules between neighbouring compartments [24].
Again, the RDME is rarely tractable, and so stochastic simulation algorithms
have been developed to simulate reaction-diffusion systems modelled using the
RDME framework [2, 9].

In this section, we outline the RDME framework, and the stochastic simulation
algorithm that we will use throughout this work. For simplicity, we consider
the evolution of a reaction-diffusion system in the square domain
$\Omega=[0,L]\times[0,L]$, and partition $\Omega$ using a Cartesian mesh into
non-overlapping rectangular sub-domains $C_{i}$, called compartments, such
that $\Omega=\bigcup_{i=1}^{I}C_{i}$ (Figure 1) [18]. Each compartment is of
size $\kappa h\times h$, where $\kappa$ is the aspect ratio of a compartment,
and the total number of compartments, $I$, can be expressed as $I=n_{x}\times
n_{y}$, where $n_{x}$ and $n_{y}$ are the numbers of compartments in the
$x$-direction and $y$-direction, respectively. We number the compartments
starting at the bottom-left corner of the domain and then from left to right,
so that, for site $i$, the adjacent compartments have indices $i\pm 1$ in the
$x$-direction and indices $i\pm n_{x}$ in the $y$-direction.

We consider a reaction-diffusion system consisting of $L$ species, where
$U_{i}^{\ell}(t)$ is the number of molecules of species $\ell$ in compartment
$i$ at time $t$, for $\ell=1,\ldots,L$ and $i=1,\ldots,I$. We define
$\bm{U}^{\ell}(t)=(U_{1}^{\ell}(t),\ldots,U^{\ell}_{I}(t))$ to be the spatial
vector of species $\ell$ at time $t$,
$\bm{U}(t)=(\bm{U}^{1}(t),\ldots,\bm{U}^{L}(t))$ to be the state matrix of the
system and
$\mathbb{P}(\bm{u},t)=\mathbb{P}(\bm{U}(t)=\bm{u}|\bm{U}(0)=\bm{u}_{0})$.
Changes to the state matrix occur through either reactions or diffusive jumps.
We consider a system of $r=1,\ldots,R$ chemical reactions such that reaction
$r$ has propensity function $a_{i}^{r}$ in compartment $i$ and is defined such
that

$\displaystyle a_{i}^{r}(\bm{u}(t))\text{d}{t}$ $\displaystyle\coloneqq$
$\displaystyle\mathbb{P}\left(\text{reaction $r$ fires in compartment $i$
during $[t,t+\text{d}{t})$}\right.$ (1)
$\displaystyle\left.\qquad\quad\text{given the system is in state $\bm{u}$ at
time $t$}\right),$

where $\text{d}{t}$ is an infinitesimal time interval. For more details of the
specific forms of the propensity functions see [24]. Similarly, diffusive
jumps have propensity function $d_{i,j}^{\ell}$ where

$d_{i,j}^{\,\ell}(\bm{u}(t))\text{d}{t}=\lambda^{\ell}_{i,j}u_{i}^{\ell}(t)\text{d}{t},$
(2)

is the probability that a molecule of species $\ell$ jumps from compartment
$i$ to compartment $j$ in the infinitesimal time interval $[t,t+\text{d}{t})$
given the system is in state $\bm{u}$ at time $t$.

As such, the RDME can be written as

$\displaystyle\frac{\partial}{\partial{t}}\mathbb{P}\left(\bm{u},t\right)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{I}\sum_{r=1}^{R}\left[a_{i}^{r}\left(\textbf{u}-\bm{\nu}_{i,r}\right)\mathbb{P}\left(\bm{u}-\bm{\nu}_{i,r},t\right)-a_{i}^{r}\left(\textbf{u}\right)\mathbb{P}\left(\bm{u},t\right)\right]$
(3) $\displaystyle\quad+\sum_{i=1}^{I}\sum_{\begin{subarray}{c}j=1\\\
j\neq{i}\end{subarray}}^{I}\sum_{\ell=1}^{L}\left[\lambda^{\ell}_{i,j}\left(u_{i}^{\ell}(t)+1\right)\mathbb{P}\left(\bm{u}-\bm{e}_{i,j}^{\ell},t\right)-\lambda^{\ell}_{i,j}u_{i}^{\ell}(t)\mathbb{P}\left(\bm{u},t\right)\right],$

where $\bm{\nu}_{i,r}$ is the stochiometric matrix that describes the change
in state upon firing of reaction $r$ in compartment $i$, and
$\bm{e}_{i,j}^{\ell}$ is the stochiometric matrix that describes the change in
state upon a diffusive jump of a molecule of species $\ell$ from compartment
$i$ to compartment $j$ [12].

To generate sample paths for a reaction-diffusion system within the RDME
framework described in Equation (3), we use an SSA called the next sub-volume
method [2, 9]. Let $t_{i}$ be the time until the next reaction within
compartment $i$ for $i=1,\ldots,I$, and let the total propensity for
compartment $i$ be

$a_{i}^{0}(\bm{u}(t))=\sum_{k=1}^{K}a_{i}^{k}(\bm{u}(t))+\sum_{\begin{subarray}{c}j=1\\\
j\neq{i}\end{subarray}}^{I}\sum_{\ell=1}^{L}{d_{i,j}^{\,\ell}(\bm{u}(t))}.$
(4)

Then, given that $T$ is the end time of a simulation, we can simulate sample
paths consistent with the RDME in Equation (3) using the SSA given in
Algorithm 1.

Algorithm 1 Next sub-volume method SSA

Set $t=0$.

Initialise $\bm{u}$.

for $i$ in $\left\\{1,2,\ldots,I\right\\}$ do

Initialise the time until the next reaction for compartment $i$:
$t_{i}\sim\exp(a_{i}^{0}(\bm{u}(t)))$.

end for

while $t<T$ do

Find argmin of the set $\left\\{t_{1},t_{2},\ldots,t_{I}\right\\}$ and denote
it $m$.

Set $t=t_{m}$.

Generate a random number $r\sim\textit{U}(0,1)$.

Choose reaction / diffusion $q$ to fire. Note that reaction $k$ fires with
probability $a_{i}^{k}(\bm{u}(t))/a_{i}^{0}(\bm{u}(t))$ and diffusion of a
molecule of species $\ell$ from box $i$ to box $j$ occurs with probability
$d_{i,j}^{\,\ell}(\bm{u}(t))/a_{i}^{0}(\bm{u}(t))$.

Update molecule numbers: $\bm{u}:=\bm{u}+\bm{\nu}_{q}$, where $\bm{\nu}_{q}$
is the stoichiometric matrix of reaction / diffusion $q$.

Let $\Gamma$ be the set of indices of non-zero elements of $\bm{\nu}_{q}$.

for $\gamma$ in $\Gamma$ do

Update $a_{i}^{k}(\bm{u}(t))$, $d_{i,j}^{\,\ell}(\bm{u}(t))$ and
$a_{i}^{0}(\bm{u})$ accordingly.

Generate a random number $\Delta t\sim\exp(a_{\gamma}^{0}(\bm{u}))$.

Update the time until the next reaction: $t_{\gamma}=t+\Delta t$.

end for

end while

Although the RDME as outlined in Equation (3) and the accompanying text,
together with the SSA of Algorithm 1, holds in full generality, for the rest
of this work, we will assume for simplicity that all diffusion coefficients
and reaction rates are spatially homogeneous.

###  2.1 Derivations of the diffusive jump rates on a static domain

We first consider a diffusion-only system with a single molecular species, U.
As already outlined in Section 2, diffusion is modelled as a series of jumps
between neighbouring compartments where the propensity function for a jump
from compartment $i$ to compartment $j$ is
$d_{i,j}(\bm{u}(t))=\lambda_{i,j}u_{i}(t)$. Note that we have dropped the
superscript $\ell$ for clarity of notation since we are currently considering
only a single species. In this work we will assume that $\lambda_{i,j}=0$,
$i=1,\ldots,I$, unless $j$ is one of the eight nearest neighbour compartments
(indices $i\pm 1$, $i\pm{n_{x}}$ and $i\pm(n_{x}\pm 1)$), as depicted in
Figure 1. In this section we outline four methods for derivation of the
$\lambda_{i,j}$.

####  2.1.1 Diffusion of a single molecule

Consider a single molecule diffusing within the domain
$\Omega=[0,L]\times[0,L]$ with diffusion coefficient $D$. Then evolution of
the position, $\textbf{x}(t)=(x(t),y(t))$, of this molecule can be described
using the following stochastic differential equations (SDEs) [6]:

$\displaystyle x(t+\text{d}{t})$ $\displaystyle=$ $\displaystyle
x(t)+\sqrt{2D}\text{d}{W_{x}},$ (5) $\displaystyle y(t+\text{d}{t})$
$\displaystyle=$ $\displaystyle y(t)+\sqrt{2D}\text{d}{W}_{y},$ (6)

where $W_{x}$ and $W_{y}$ are the usual Wiener processes. The associated
Fokker-Planck equation (FPE) is

$\frac{\partial p(\textbf{x},t)}{\partial t}=D\Delta
p(\textbf{x},t)\quad\text{for}\quad\textbf{x}\in\Omega,$ (7)

so that $p(\textbf{x},t)\text{d}{\textbf{x}}$ is the probability to find the
molecule in the region $[x,x+\text{d}{x}]\times[y,y+\text{d}{y}]$ at time $t$.

For simplicity, throughout this work we will use Neumann boundary conditions

$\left(\textbf{n}\cdot\nabla\right)p=0\quad\text{for}\quad\textbf{x}\in\partial\Omega,$
(8)

and we take the initial distribution for the molecule to be

$p(\textbf{x},0)=p_{0}(\textbf{x}).$ (9)

Using the Cartesian mesh shown in Figure 1, we discretise the Laplacian as
follows:

$D\Delta
p(\textbf{x}_{i},t)\approx\sum_{j=1}^{8}\lambda_{i,j}\,\tilde{p}_{j}(t)-\sum_{j=1}^{8}\lambda_{j,i}\,\tilde{p}_{i}(t),$
(10)

where $\tilde{p}_{i}(t)$ is the approximate solution to Equation (7) at point
$\textbf{x}_{i}$, i.e. at the centre of compartment $i$ at time $t$. The
compartments are all of the same size, $\kappa h\times h$, with area $A=\kappa
h^{2}$ (Figure 1). Hence, exploiting translational symmetry and noting that
the diffusion coefficient, $D$, is spatially homogeneous, we refer to the jump
rates, $\lambda_{i,j}$, from now on as $\lambda_{n}$ for
$n\in\left\\{1,2,\ldots,8\right\\}$ (Figure 1), with

$\lambda_{0}=\sum_{n=1}^{8}\lambda_{n}.$ (11)

Figure 1: Schematic diagram of the discretisation of the two-dimensional
domain $\Omega$, with the dimensions of the compartments and the indexing of
the jump rates, $\lambda_{n}$, for $n\in\left\\{1,\ldots,8\right\\}$ for a
given compartment $i$ shown.

We now make a few remarks. First, note that this description allows for jumps
between diagonally neighbouring compartments. Second, the indexing of the jump
rates starts at one for a jump to the right and continues in an anti-clockwise
fashion, as shown in Figure 1. Third, due to the reflective symmetries of a
Cartesian mesh and since $D$ is spatially homogeneous, we have the following
identities

$\lambda_{1}=\lambda_{5}\,,\quad\lambda_{3}=\lambda_{7}\,,\quad\lambda_{2}=\lambda_{4}=\lambda_{6}=\lambda_{8}.$
(12)

Finally, note that in what follows we present derivations of the jump rates
for a two-dimensional domain but that all derivations naturally extend to the
three-dimensional case.

####  2.1.2 The finite difference method

The first method we consider to derive the diffusive jump rates is the FDM.
Meinecke and Lötsedt [18] chose a 9-point stencil which allows one to include
a parameter $\alpha$ to study the effects of diagonal jumps, while maintaining
a discretisation of the Laplacian that is still second order accurate [26].
The standard 5-point stencil can be recovered by setting $\alpha=0$ and hence
neglecting diagonal jumps. We have

$\Delta
p\left(\textbf{x}_{i},t\right)\approx\delta^{2}_{x}\tilde{p}_{i}+\delta^{2}_{y}\tilde{p}_{i}+\frac{1}{2}\alpha\kappa
h^{2}\delta^{2}_{x}\delta^{2}_{y}\tilde{p}_{i},$ (13)

where $\alpha$ is a free parameter controlling the frequency of diagonal
jumps, and the difference operators are defined as follows:

$\displaystyle\delta^{2}_{x}\tilde{p}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{\kappa^{2}h^{2}}\left(\tilde{p}_{i+1}+\tilde{p}_{i-1}-2\tilde{p}_{i}\right);$
$\displaystyle\delta^{2}_{y}\tilde{p}_{i,j}$ $\displaystyle=$
$\displaystyle\frac{1}{h^{2}}\left(\tilde{p}_{i+n_{x}}+\tilde{p}_{i-n_{x}}-2\tilde{p}_{i}\right).$
(14)

As such, the right-hand side of Equation (7) can be approximated as

$\displaystyle D\Delta p$ $\displaystyle\approx$
$\displaystyle\frac{D-D\alpha\kappa}{\kappa^{2}h^{2}}\left(\tilde{p}_{i+1}+\tilde{p}_{i-1}\right)+\frac{D\kappa-D\alpha}{\kappa
h^{2}}\left(\tilde{p}_{i+n_{x}}+\tilde{p}_{i-n_{x}}\right)$ (15)
$\displaystyle\quad+\frac{D\alpha}{2\kappa
h^{2}}\left(\tilde{p}_{i+n_{x}+1}+\tilde{p}_{i-n_{x}+1}+\tilde{p}_{i+n_{x}-1}+\tilde{p}_{i-n_{x}-1}\right)$
$\displaystyle\quad-2D\frac{\kappa^{2}-\alpha\kappa+1}{\kappa^{2}h^{2}}\tilde{p}_{i}\,.$

Comparing the general expression for the discretised Laplacian in Equation
(10) with the discretised Laplacian in Equation (15) we have

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{D-D\alpha\kappa}{\kappa^{2}h^{2}},$ (16)
$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{D\alpha}{2\kappa h^{2}},$ (17) $\displaystyle\lambda_{3}$
$\displaystyle=$ $\displaystyle\frac{D\kappa-D\alpha}{\kappa h^{2}},$ (18)
$\displaystyle\lambda_{0}$ $\displaystyle=$
$\displaystyle\frac{2D(\kappa^{2}-\alpha\kappa+1)}{\kappa^{2}h^{2}}.$ (19)

####  2.1.3 The finite volume method

We now consider derivation of the jump rates using the FVM. We start by
averaging the Laplacian at $\textbf{x}_{i}$ by integrating over compartment
$i$

$D\Delta p(\textbf{x}_{i})\approx\frac{D}{A}\int_{C_{i}}\Delta
p\,\text{d}\Omega,$ (20)

where $A$ is the area of the compartment $i$. Using the divergence theorem, we
rewrite Equation (20) as

$D\Delta p(\textbf{x}_{i})\approx\frac{D}{A}\int_{\partial
C_{i}}\left(\textbf{n}\cdot\nabla\right)p\,\text{d}s,$ (21)

where $\text{d}s$ is a line element of the compartment boundary, $\partial
C_{i}$, and n is the normal vector to the compartment boundary at the point x.

Next, we split the integral in Equation (21) into integrals along segments of
compartment boundary $\partial C_{i}$, as shown in Figure 2:

$D\Delta p(\textbf{x}_{i})\approx\frac{D}{A}\sum_{k=1}^{4}\int_{\partial
C_{i}^{k}}\left(\textbf{n}\cdot\nabla\right)p\,\text{d}s,$ (22)

and we approximate the gradient of $p$ across the different parts of the
boundary of the compartment using a finite difference method:

$\displaystyle\textbf{n}\cdot\nabla p\left\rvert{}_{\partial
C_{i}^{1}}\right.$ $\displaystyle\approx$
$\displaystyle\frac{\tilde{p}_{i+1}-\tilde{p}_{i}}{\kappa h};$ (23)
$\displaystyle\textbf{n}\cdot\nabla p\left\rvert{}_{\partial
C_{i}^{2}}\right.$ $\displaystyle\approx$
$\displaystyle\frac{\tilde{p}_{i+n_{x}}-\tilde{p}_{i}}{h};$ (24)
$\displaystyle\textbf{n}\cdot\nabla p\left\rvert{}_{\partial
C_{i}^{3}}\right.$ $\displaystyle\approx$
$\displaystyle\frac{\tilde{p}_{i-1}-\tilde{p}_{i}}{\kappa h};$ (25)
$\displaystyle\textbf{n}\cdot\nabla p\left\rvert{}_{\partial
C_{i}^{4}}\right.$ $\displaystyle\approx$
$\displaystyle\frac{\tilde{p}_{i-n_{x}}-\tilde{p}_{i}}{h}.$ (26)

Substituting Equations (23)–(26) into Equation (22) gives

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{D}{\kappa^{2}h^{2}},$ (27) $\displaystyle\lambda_{2}$
$\displaystyle=$ $\displaystyle 0,$ (28) $\displaystyle\lambda_{3}$
$\displaystyle=$ $\displaystyle\frac{D}{h^{2}},$ (29)
$\displaystyle\lambda_{0}$ $\displaystyle=$
$\displaystyle\frac{2D(\kappa^{2}+1)}{\kappa^{2}h^{2}}.$ (30)

Figure 2: Schematic diagram of how compartment $i$ boundary, $\partial C_{i}$,
is split into four segments $\partial C_{i}^{k}$ for $k=1,2,3,4$.

####  2.1.4 The finite element method

The last numerical method we will use to derive the jump coefficients is the
FEM. We multiply Equation (7) by a test function $v(\textbf{x})$, which is a
function in Sobolev space [3] (i.e. is bounded and has bounded first partial
derivatives), and integrate over the whole domain $\Omega$:

$\displaystyle\int_{\Omega}\frac{\partial p}{\partial
t}v\,\text{d}x\,\text{d}y$ $\displaystyle=$ $\displaystyle
D\int_{\Omega}v\Delta p\,\text{d}x\,\text{d}y$ (31) $\displaystyle=$
$\displaystyle D\int_{\Omega}\left[\nabla\cdot\left(v\nabla p\right)-\nabla
p\cdot\nabla v\right]\,\text{d}x\,\text{d}y.$

Now we use the divergence theorem to give

$\int_{\Omega}\frac{\partial p}{\partial
t}v\,\text{d}x\,\text{d}y=D\int_{\partial\Omega}\textbf{n}\cdot\left[v\nabla
p\right]\text{d}s-D\int_{\Omega}\nabla p\cdot\nabla v\,\text{d}x\,\text{d}y,$
(32)

where $\text{d}s$ is the line element of the domain boundary,
$\partial\Omega$, and n is the normal vector to the domain boundary at the
point x. Due to the Neumann boundary conditions, the first term in Equation
(32) vanishes and we are left with

$\int_{\Omega}\frac{\partial p}{\partial
t}v\,\text{d}x\,\text{d}y=-D\int_{\Omega}\nabla p\cdot\nabla
v\,\text{d}x\,\text{d}y.$ (33)

We express the functions $p$, $\partial p/\partial t$ and $v$ in terms of the
basis functions $\varphi_{j}$,

$\displaystyle v(x,y)$ $\displaystyle=$
$\displaystyle\varphi_{j}(x,y),\vphantom{\sum_{i=1}^{I}}$ (34) $\displaystyle
p(x,y,t)$ $\displaystyle\approx$
$\displaystyle\sum_{i=1}^{I}\tilde{p}_{i}(t)\varphi_{i}(x,y),$ (35)
$\displaystyle\frac{\partial p}{\partial t}$ $\displaystyle\approx$
$\displaystyle\sum_{i=1}^{I}\frac{\text{d}\tilde{p}_{i}}{\text{d}t}(t)\varphi_{i}(x,y),$
(36)

where the $\varphi_{i}$ are defined as

$\varphi_{i}\left(x,y\right)=\left\\{\begin{array}[]{cc}\left(1-\cfrac{x-x_{i}}{\kappa
h}\right)\left(1-\cfrac{y-y_{i}}{h}\right)&\text{for}\quad\textbf{x}\in\left[x_{i},x_{i}+\kappa
h\right]\times\left[y_{i},y_{i}+h\right],\\\ \left(1+\cfrac{x-x_{i}}{\kappa
h}\right)\left(1-\cfrac{y-y_{i}}{h}\right)&\text{for}\quad\textbf{x}\in\left[x_{i}-\kappa
h,x_{i}\right]\times\left[y_{i},y_{i}+h\right],\\\
\left(1+\cfrac{x-x_{i}}{\kappa
h}\right)\left(1+\cfrac{y-y_{i}}{h}\right)&\text{for}\quad\textbf{x}\in\left[x_{i}-\kappa
h,x_{i}\right]\times\left[y_{i}-h,y_{i}\right],\\\
\left(1-\cfrac{x-x_{i}}{\kappa
h}\right)\left(1+\cfrac{y-y_{i}}{h}\right)&\text{for}\quad\textbf{x}\in\left[x_{i},x_{i}+\kappa
h\right]\times\left[y_{i}-h,y_{i}\right],\\\
0&\text{otherwise},\end{array}\right.$ (37)

for $i\in\left\\{1,\ldots,I\right\\}$, and $\varphi_{j}(\textbf{x})=0$ for
$\textbf{x}\notin\Omega$. We choose linear basis functions for simplicity,
however, any form will suffice, as long as the following condition is met:

$\varphi_{j}\left(x_{i},y_{i}\right)=\begin{cases}1\quad\text{if}\quad i=j;\\\
0\quad\text{if}\quad i\neq j.\end{cases}$ (38)

Substituting Equations (34), (35) and (36) into Equation (33) gives, for each
$j$,

$\sum_{i=1}^{I}\left(\int_{\Omega}\varphi_{i}\varphi_{j}\,\text{d}x\,\text{d}y\right)\frac{\text{d}\tilde{p}_{i}}{\text{d}t}=-\sum_{i=1}^{I}\left(\int_{\Omega}D\nabla\varphi_{i}\cdot\nabla\varphi_{j}\,\text{d}x\,\text{d}y\right)\tilde{p}_{i},$
(39)

which can be expressed as a matrix equation

$M\frac{\text{d}\tilde{\textbf{p}}}{\text{d}t}=-N\tilde{\textbf{p}}\,,$ (40)

where

$\displaystyle\tilde{\textbf{p}}$ $\displaystyle=$
$\displaystyle(\tilde{p}_{1},\tilde{p}_{2},\ldots,\tilde{p}_{I}),\vphantom{\int_{\Omega}}$
(41) $\displaystyle M_{ij}$ $\displaystyle=$
$\displaystyle\int_{\Omega}\varphi_{i}\,\varphi_{j}\,\text{d}x\,\text{d}y,$
(42) $\displaystyle N_{ij}$ $\displaystyle=$
$\displaystyle\int_{\Omega}D\nabla\varphi_{i}\cdot\nabla\varphi_{j}\,\text{d}x\,\text{d}y.$
(43)

Evaluating the elements of the matrices $M$ and $N$, and using the lumped mass
matrix approach (see Supplementary Information), we finally arrive at the
following expressions for the jump coefficients:

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{D(2-\kappa^{2})}{3\kappa^{2}h^{2}};$ (44)
$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{D(\kappa^{2}+1)}{6\kappa^{2}h^{2}};$ (45)
$\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{D(2\kappa^{2}-1)}{3\kappa^{2}h^{2}};$ (46)
$\displaystyle\lambda_{0}$ $\displaystyle=$
$\displaystyle\frac{4D(\kappa^{2}+1)}{3\kappa^{2}h^{2}}.$ (47)

(r,0cm,0.1cm)[0.5cm]—c—c—c—c——c—c—c—c—Jump rates & FDM FVM FEM
$\lambda_{1}$ $\displaystyle\frac{D-D\alpha\kappa}{\kappa^{2}h^{2}}$
$\displaystyle\frac{D}{\kappa^{2}h^{2}}$
$\displaystyle\frac{D(2-\kappa^{2})}{3\kappa^{2}h^{2}}$
$\lambda_{2}$ $\displaystyle\frac{D\alpha}{2\kappa h^{2}}$ $\displaystyle 0$
$\displaystyle\frac{D(\kappa^{2}+1)}{6\kappa^{2}h^{2}}$
$\lambda_{3}$ $\displaystyle\frac{D\kappa-D\alpha}{\kappa h^{2}}$
$\displaystyle\frac{D}{h^{2}}$
$\displaystyle\frac{D(2\kappa^{2}-1)}{3\kappa^{2}h^{2}}$

Table 1: Comparison of jump rates from different derivations.

####  2.1.5 The first exit time method

Meinecke and Löstedt proposed an additional method of calculating the jump
rates [16, 18]. This method involves considering the FET distribution from a
domain $\omega_{i}$ (centered on $C_{i}$) of a single molecule undergoing
Brownian motion that is placed initially at the centre of $C_{i}$. Note that
to model diffusion in a RDME framework using FET-derived jump rates we
approximate the probability distribution function (PDF) for the FET by an
exponential distribution with the same mean. We follow the choice of domain
$\omega_{i}$ from [18], where $\omega_{i}$ is defined by its boundary,
$\partial\omega_{i}$, and taken to be a rectangle that intersects the centres
of neighbouring compartments of compartment $i$, as shown in Figure 3, so that
$\omega_{i}$ completely contains the compartment $i$111The choice of
$\omega_{i}$ has been briefly discussed in [16] and so we do not explore it
further here..

Figure 3: Diagram of the domain $\omega_{i}$, which has been shaded in dark
blue. The boundary of the domain $\omega_{i}$, indicated by the dashed line,
is split into eight sections coloured as follows: $\partial\omega_{1}$ is the
rightmost light blue; $\partial\omega_{2}$ is the upper-right dark blue;
$\partial\omega_{3}$ is the uppermost green; $\partial\omega_{4}$ is the
upper-left dark blue; $\partial\omega_{5}$ is the leftmost light blue;
$\partial\omega_{6}$ is the lower-left dark blue; $\partial\omega_{7}$ is the
lower green; and $\partial\omega_{8}$ is the lower-right dark blue.

We follow Meinecke and Löstedt [18] in choosing jump rates of the form

$\lambda_{n}=\lambda_{0}\,\theta_{n},$ (48)

where $\theta_{n}$ is the probability that a molecule exits $\omega_{i}$ into
the $n$th neighbouring compartment (see Figure 1). We define a molecule to
exit $\omega_{i}$ into the $n$th neighbouring compartment if it exits through
the segment of the domain boundary denoted by $\partial\omega_{i}^{n}$ in
Figure 3. There is flexibility in the choice of discretisation of the boundary
$\partial\omega_{i}$ into segments $\partial\omega_{i}^{n}$, $n=1,\ldots,8$,
and we define the discretisation using parameter $\beta$ as follows:

$\begin{split}&\left|\partial\omega_{1}\right|=\left|\partial\omega_{5}\right|=2\beta
h;\\\
&\left|\partial\omega_{3}\right|=\left|\partial\omega_{1}\right|=2\beta\kappa
h;\\\
&\left|\partial\omega_{2}\right|=\left|\partial\omega_{4}\right|=\left|\partial\omega_{6}\right|=\left|\partial\omega_{8}\right|=(1-\beta)(1+\kappa)h.\end{split}$
(49)

We will explore the effect of changing the parameter $\beta$ on the jump rates
in Section 2.1.6.

We proceed by considering two random variables associated with a molecule
governed by Brownian motion exiting from the domain $\omega_{i}$: $\tau$ and
$\eta$. The random variable $\tau$ denotes the FET out of $\omega_{i}$,
whereas $\eta$ denotes the exit segment. We use $\tau$ to calculate
$\lambda_{0}$ by, as explained earlier, equating the mean of the PDF of $\tau$
with the rate parameter of the exponential distribution used in Algorithm 1.
The random variable $\eta$, on the other hand, is used to calculate
$\theta_{n}$.

We calculate the PDF of $\tau$ using the FPE (as outlined in Section 2.1) that
describes the probability for the molecule to be at position x at time $t$:

$\frac{\partial p}{\partial t}=D\Delta
p\quad\text{for}\quad\textbf{x}\in\omega_{i},$ (50)

where $D$ is the diffusion coefficient. To calculate the PDF of the FET, we
use Dirichlet (absorbing) boundary conditions

$p(\textbf{x},t)=0\quad\text{for}\quad\textbf{x}\in\partial\omega_{i},$ (51)

and assume that initially the molecule at the centre of $\omega_{i}$ so that

$p(\textbf{x},0)=\delta(\textbf{x}-\textbf{x}_{i}).$ (52)

Using $p(\textbf{x},t)$ we calculate the survival probability, $S(t)$, which
is the probability that a molecule remains in the domain $\omega_{i}$ at least
until time $t$, as [22]

$S(t)=\mathbb{P}(\tau\geqslant
t)=\int_{\omega_{i}}p(\textbf{x},t)\,\text{d}\omega,$ (53)

where $\text{d}\omega$ is a domain element of the domain $\omega_{i}$. Now,
given $S(t)$, we derive the probability that a molecule exits $\omega_{i}$ at
time $t$, i.e. the PDF of the random variable $\tau$, denoted by
$\mathbb{P}(\tau=t)$, as

$\mathbb{P}(\tau=t)=-\frac{\partial S}{\partial t}\,.$ (54)

Substituting Equation (53) into Equation (54), results in

$\displaystyle\mathbb{P}(\tau=t)$ $\displaystyle=$
$\displaystyle-\int_{\omega_{i}}\frac{\partial p}{\partial t}\,\text{d}\omega$
(55) $\displaystyle=$ $\displaystyle-D\int_{\omega_{i}}\Delta
p\,\text{d}\omega$ $\displaystyle=$
$\displaystyle-D\int_{\partial\omega_{i}}\textbf{n}\cdot\nabla p\,\text{d}s,$

where $\text{d}s$ is a line element of $\partial\omega_{i}$. Finally, we
calculate the mean of the PDF for $\tau$ as

$E_{\tau}=\mathbb{E}\left[\tau\right]=\int_{0}^{\infty}t\,\mathbb{P}(\tau=t)\,\text{d}t=\int_{0}^{\infty}S(t)\,\text{d}t,$
(56)

and set

$\lambda_{0}=\frac{1}{E_{\tau}}.$ (57)

To calculate the PDF for $\eta$, we use the law of total probability

$\mathbb{P}(\eta=n)=\int_{0}^{\infty}\mathbb{P}(\eta=n|\tau=t)\mathbb{P}(\tau=t)\,\text{d}t.$
(58)

Following [18], we have the probability of a molecule exiting through segment
$\partial\omega_{i}^{n}$ given that a molecule exits at time $t$, i.e. the
probability $\mathbb{P}(\eta=n|\tau=t)$, as follows

$\mathbb{P}(\eta=n|\tau=t)=\frac{\int_{\partial\omega_{i}^{n}}\textbf{n}\cdot\nabla
p\,\text{d}s}{\int_{\partial\omega_{i}}\textbf{n}\cdot\nabla p\,\text{d}s}.$
(59)

Therefore, substituting Equation (55) and Equation (59) into Equation (58), we
have

$\displaystyle\theta_{n}$ $\displaystyle=$
$\displaystyle\mathbb{P}(\eta=n)\vphantom{\int_{0}^{\infty}}$ (60)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{\int_{\partial\omega_{i}^{n}}\textbf{n}\cdot\nabla
p\,\text{d}s}{\int_{\partial\omega_{i}}\textbf{n}\cdot\nabla
p\text{d}s}\left(-D\int_{\partial\omega_{i}}\textbf{n}\cdot\nabla
p\,\text{d}s\right)\text{d}t$ $\displaystyle=$
$\displaystyle-D\int_{0}^{\infty}\int_{\partial\omega_{i}^{n}}\textbf{n}\cdot\nabla
p\,\text{d}s\,\text{d}t.$

Therefore, given solutions to Equation (48) and Equation (50), we can
calculate the jump rates $\lambda_{n}$ for $n=1,\ldots,8$. To see the full
expressions for the jump coefficients derived using the FET method see the
Supplementary Information.

####  2.1.6 Comparison of the diffusive jump rates

Figure 4: Plot of the jump rates as a function of the implementation
parameters, $\alpha$ and $\beta$, for $h=0.025\,\mu m$, $\kappa=1.0$ and
$D=1.0\,\mu m^{2}s^{-1}$ for the total jump rate, $\lambda_{0}$, as well as
$\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$. Since the derivations using
the FEM and the FVM do not have any implementation parameters, the jump rates
are constant for the given simulation parameters.

Figure 4 shows how the four methods of derivation compare to one another as
the parameters $\alpha$ and $\beta$, introduced during the derivations of the
FDM and FET jump rates, respectively, are varied. Both $\alpha$ and $\beta$
control the frequency of diagonal jumps. For the FDM-derived jump rates,
increasing the parameter $\alpha$ increases $\lambda_{2}$, the rate of
diagonal jumps. On the other hand, increasing $\beta$ has the opposite effect
on $\lambda_{2}$ in the FET case. As the FEM- and the FVM-derived jump rates
do not have any free parameters, they are constant in each of the plots.

###  2.2 Derivation of the diffusive jump rates on a growing domain

The SDE description of a molecule undergoing Brownian motion on a two-
dimensional growing domain $\Omega_{t}=[0,L(t)]\times[0,L(t)]$ we use in this
work is

$\displaystyle x(t+\text{d}t)$ $\displaystyle=$ $\displaystyle
x(t)+v_{x}(x,y,t)\text{d}t+\sqrt{2D}\text{d}{W}_{x},$ (61) $\displaystyle
y(t+\text{d}t)$ $\displaystyle=$ $\displaystyle
y(t)+v_{y}(x,y,t)\text{d}t+\sqrt{2D}\text{d}{W}_{y},$ (62)

where the drift terms, $v_{x}(x,y,t)$ and $v_{y}(x,y,t)$, describe the effects
of domain growth. The corresponding FPE is

$\frac{\partial p}{\partial
t}=D\nabla^{2}p-\nabla\cdot(\textbf{v}p)\quad\text{for}\quad\textbf{x}\in\Omega_{t},$
(63)

with Neumann boundary conditions

$\left.\nabla p\right|_{\partial\Omega_{t}}=0,$ (64)

and initial condition

$p(\textbf{x},0)=p_{0}(\textbf{x})\quad\text{for}\quad\textbf{x}\in\Omega_{0}.$
(65)

Assuming growth is isotropic, the drift terms,
$\textbf{v}(\textbf{x},t)=(v_{x}(x,y,t),v_{y}(x,y,t))$, are given by

$\textbf{v}(\textbf{x},t)=\left(\frac{x\dot{L}(t)}{L(t)},\frac{y\dot{L}(t)}{L(t)}\right).$
(66)

Care now must be taken when considering the equation that is to be used to
calculate the diffusive jump rates, since on the growing domain the
compartments are also changing in size, with $x$ and $y$ dimensions
$\kappa{h(t)}=L(t)/n_{x}$ and $h(t)=L(t)/n_{y}$, respectively. This means that
the probability to find the particle in compartment $C_{i}(t)$ is given by
$\int_{C_{i}(t)}p(\textbf{x},t)\text{d}\textbf{x}$ with

$\displaystyle\frac{\text{d}}{\text{d}{t}}\int_{C_{i}(t)}p(\textbf{x},t)\,\text{d}\textbf{x}$
$\displaystyle=$ $\displaystyle\int_{C_{i}(t)}\frac{\partial}{\partial
t}p(\textbf{x},t)\,\text{d}\textbf{x}+\int_{\partial
C_{i}(t)}(\textbf{v}\cdot\textbf{n})\,p(\textbf{x},t)\,\text{d}\textbf{x}.$
(67)

Using equation (63) together with the divergence theorem, we see that the
terms involving the drift, v, cancel so that we should consider the particle
only as diffusing on the growing domain. We employ a change of coordinate to
transform to a static domain,

$(x,y,t)\rightarrow(\xi,\eta,\tau)=\left(\frac{x}{L(t)},\frac{y}{L(t)},t\right).$
(68)

so that we can derive the diffusive jump rates using the diffusion equation
with time-varying diffusion coefficient:

$\frac{\partial
p}{\partial\tau}=\frac{D}{L(t)^{2}}\nabla_{\xi}\,p\quad\text{for}\quad\bm{\xi}\in\Omega_{0},$
(69)

where $\bm{\xi}=(\xi,\eta)$ and
$\nabla_{\xi}=(\partial/\partial\xi,\partial/\partial\eta)$.

The key advantage is that we can now simply use the diffusive jump rate
derivations for the static case for the FDM, FVM and FEM (see Table 2). Note
that we no longer consider the FET approach to deriving the diffusive jump
rates. Simpson et al. [25] have previously shown that on a growing domain the
survival probability does not approach zero in the limit $t\to\infty$. As a
result, the FET distribution has infinite mean so we cannot use it to
parameterise the exponential distribution as required in Equation (57).

(r,0cm,0.1cm)[0.5cm]—c—c—c—c——c—c—c—c—Jump rates & FDM FVM FEM
$\lambda_{1}(t)$ $\displaystyle\frac{D-D\alpha\kappa}{\kappa^{2}h^{2}(t)}$
$\displaystyle\frac{D}{\kappa^{2}h^{2}(t)}$
$\displaystyle\frac{D(2-\kappa^{2})}{3\kappa^{2}h^{2}(t)}$
$\lambda_{2}(t)$ $\displaystyle\frac{D\alpha}{2\kappa h^{2}(t)}$
$\displaystyle 0$ $\displaystyle\frac{D(\kappa^{2}+1)}{6\kappa^{2}h^{2}(t)}$
$\lambda_{3}(t)$ $\displaystyle\frac{D\kappa-D\alpha}{\kappa h^{2}(t)}$
$\displaystyle\frac{D}{h^{2}(t)}$
$\displaystyle\frac{D(2\kappa^{2}-1)}{3\kappa^{2}h^{2}(t)}$

Table 2: Comparison of jump rates from different derivations on a growing
domain.

##  3 Results

The SSA outlined in Algorithm 1 of Section 2 allows us to simulate a range of
reaction-diffusion systems and test the impact of the different diffusive jump
derivations on model predictions. We first explore a range of systems on a
static domain, before moving to a growing domain.

###  3.1 Static domain

For simple systems, e.g. those that contain only zeroth and first order
reactions, it is possible to compare results with analytic solutions of the
corresponding macroscopic partial differential equation (PDE) models (see
Supplementary Information for more details of the models and their analytic
solution). The PDE represents the evolution of molecule concentration in time,
so to compare with the results of stochastic simulation, we use the following
error measure:

$e(t)=\sqrt{\sum_{i=1}^{I}A(t)\left|\frac{U_{i}(t)}{A(t)}-u(\textbf{x}_{i},t)\right|^{2}},$
(70)

where $A(t)$ is the compartment area at time $t$, and $u(\textbf{x}_{i},t)$ is
the analytical solution of the corresponding macroscale PDE at
$\textbf{x}_{i}$ at time $t$.

Figure 5: Plot of the error, as defined in Equation (70), as a function of
$\kappa$, the compartment aspect ratio, for simulations with (a) diffusion
only, (b) diffusion and production, (c) diffusion and decay and (d) diffusion,
decay and production. The domain $\Omega=[0,20]\times[0,20]$ is discretised
into $21\times 21\kappa$ compartments and $5\times 10^{6}$ molecules are
placed initially in the central compartment. Parameters are: $D=1\,\mu
m^{2}s^{-1}$, $\alpha=0.7$, $\beta=0.5$, $T=5\,s$ and $k_{1}=100\,\mu
m^{-2}s^{-1}$, $k_{-1}=0.1\,s^{-1}$ for the systems involving production and
decay, respectively.

####  3.1.1 Diffusion

We first investigate the effects of varying the compartment aspect ratio,
$\kappa$. The macroscale PDE that describes the evolution of molecule
concentration is

$\frac{\partial u(\textbf{x},t)}{\partial t}=D\Delta
u(\textbf{x},t)\quad\text{for}\quad\textbf{x}\in\Omega,$ (71)

with Neumann boundary conditions

$\left(\textbf{n}\cdot\nabla\right)p=0\quad\text{for}\quad\textbf{x}\in\partial\Omega,$
(72)

where n is the normal vector to the domain boundary $\partial\Omega$. We
choose the initial condition to be

$u(\textbf{x},0)=\left\\{\begin{array}[]{cc}5\times 10^{6}/A&\text{for
}0<x<\kappa h,0<y<h,\\\ 0&\text{otherwise}.\end{array}\right.$ (73)

We use the analytic solution of this PDE (see Supplementary Information) to
calculate the error, $e$, using Equation (70). Figure 5(a) shows that the
error is relatively unaffected by variation in the compartment aspect ratio,
$\kappa$. For the values of $\alpha$ and $\beta$ used, the FVM-derived jump
rates give rise to the smallest error and all other methods (FDM, FEM and FET)
give rise to errors that are very similar to each other. Note that for
$\alpha=0$ the diffusive jump rates for the FDM and the FVM are identical,
indicating that increasing $\alpha$ increases the error.

####  3.1.2 Zeroth and first order reactions

We now consider reaction-diffusion systems involving the decay of molecules,
given by the chemical equation

$\text{U}\xrightarrow{k_{-1}}\emptyset,$ (74)

and the production of molecules, given by the chemical equation

$\emptyset\xrightarrow{k_{1}}\text{U},$ (75)

alongside diffusion. The corresponding propensity functions take the form

$\displaystyle a_{i}^{\text{decay}}(\textbf{u}(t))$ $\displaystyle=$
$\displaystyle k_{-1}u_{i}(t),$ (76) $\displaystyle
a_{i}^{\text{prod}}(\textbf{u}(t))$ $\displaystyle=$ $\displaystyle
k_{1}A(t),$ (77)

for the decay reaction and the production reaction in compartment $i$,
respectively. As before, $A(t)$ is the area of a compartment at time $t$, and
$u_{i}(t)$ is the number of molecules in compartment $i$ at time $t$. We
consider production-only, decay-only and production-and-decay systems, and the
macroscale PDEs that describe the evolution of molecule concentration are,
respectively,

$\displaystyle\frac{\partial u(\textbf{x},t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\Delta u(\textbf{x},t)+k_{1},$ (78)
$\displaystyle\frac{\partial u(\textbf{x},t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\Delta u(\textbf{x},t)-k_{-1}u,$ (79)
$\displaystyle\frac{\partial u(\textbf{x},t)}{\partial t}$ $\displaystyle=$
$\displaystyle D\Delta u(\textbf{x},t)+k_{1}-k_{-1}u.$ (80)

Again, we use Neumann boundary conditions

$\left(\textbf{n}\cdot\nabla\right)p=0\quad\text{for}\quad\textbf{x}\in\partial\Omega,$
(81)

where n is the normal vector to the domain boundary $\partial\Omega$, and
initially $5\times 10^{6}$ molecules are placed in the central compartment.

Figures 5(b)–(d) show how the error, $e$, varies with $\kappa$ for these three
systems (solutions to the PDEs in Equations (78)–(80) are given in the
Supplementary Information). Again, the FVM out-performs the other three
methods (FDM, FEM and FET) for the chosen values of $\alpha$ and $\beta$. In
addition, for systems involving decay (a first order reaction) the error
increases as the aspect ratio of the compartments, $\kappa$, increases.

####  3.1.3 Second-order reactions

To investigate the effects of different derivations of the jump rates in
systems involving second-order reactions, we now consider a system composed of
two species, U and V, which undergo diffusion (with rates $D_{u}$ and $D_{v}$,
respectively) and the following reactions

$\text{U}+\text{V}\xrightarrow{k_{1}}\text{V},\quad\emptyset\xrightarrow{k_{2}}\text{U},$
(82)

i.e. U is consumed upon contact with V at rate $k_{1}$, and U is produced at
rate $k_{2}$. The corresponding propensity functions take the form

$\displaystyle a_{i}^{1}(\textbf{u}(t))$ $\displaystyle=$
$\displaystyle\frac{k_{1}}{A}u_{i}(t)v_{i}(t),$ (83) $\displaystyle
a_{i}^{2}(\textbf{u}(t))$ $\displaystyle=$ $\displaystyle k_{2}A(t).$ (84)

We cannot write down equivalent closed-form macroscale PDEs for the evolution
of U and V concentration over time because the system contains second order
reactions. As such, we investigate how the stationary distribution, estimated
from stochastic simulations, is affected by the choice of compartment size,
for the four different methods of deriving the diffusive jump rates. In Figure
6(a) we see that as the number of compartments increases the stationary
distribution shifts to the right, towards higher molecule numbers. This is
because the rate of bimolecular reactions decreases as the number of
compartments increases [13]. Figure 6(b) shows how the stationary distribution
changes with the method of derivation of the diffusive jump rates. All methods
give slightly different results; in particular, using $\beta=0.1$ or
$\beta=0.9$ in the FEM results in stationary distributions that are shifted
noticeably to the left.

Figure 6: Stationary distributions as calculated from simulations run until
$T=1\times 10^{6}s$, on the domain $\Omega=[0,1]\times[0,1]$, with the
following reaction rates: $k_{1}=0.2\,\mu m^{2}s^{-1}$ and $k_{2}=1.0\,\mu
m^{-2}s^{-1}$. We initialised the simulations with five molecules of U and one
molecule of V distributed at random across the compartments. (a) Comparison of
the stationary distributions as the compartment number is varied, where we
used FVM-derived jump rates. (b) Comparison of the stationary distributions
for different methods of derivation of the diffusive jump rates. The domain
has been discretised into $1000\times 1000$ compartments.

####  3.1.4 Schnakenberg kinetics

Finally, we investigate the effects of different derivations of diffusive jump
rates on Turing’s reaction-diffusion model of pattern formation [27]. We focus
on Schnakenberg kinetics [7, 23] as they are simple to simulate. The system is
composed of two species, U and V, which undergo diffusion (with rates $D_{u}$
and $D_{v}$, respectively) and the following reactions

$\emptyset\xrightarrow{k_{1}}\text{U},\quad\text{U}\xrightarrow{k_{2}}\emptyset,\quad
2\text{U}+\text{V}\xrightarrow{k_{3}}3\text{U},\quad\emptyset\xrightarrow{k_{4}}\text{V}.$
(85)

The corresponding propensity functions take the form

$\displaystyle a_{i}^{1}(\textbf{u}(t),\textbf{v}(t))$ $\displaystyle=$
$\displaystyle k_{1}A(t),$ (86) $\displaystyle
a_{i}^{2}(\textbf{u}(t),\textbf{v}(t))$ $\displaystyle=$ $\displaystyle
k_{2}u_{i}(t),$ (87) $\displaystyle a_{i}^{3}(\textbf{u}(t),\textbf{v}(t))$
$\displaystyle=$
$\displaystyle\frac{k_{3}}{A(t)^{2}}u_{i}(t)(u_{i}(t)-1)v_{i}(t),$ (88)
$\displaystyle a_{i}^{4}(\textbf{u}(t),\textbf{v}(t))$ $\displaystyle=$
$\displaystyle k_{4}A(t).$ (89)

We quantify the patterns formed using the discrete Fourier cosine transform.
Let $f(x,y):\mathbb{R}^{2}\rightarrow\mathbb{R}$, then the discrete Fourier
cosine transform is defined as

$\hat{f}(k_{x},k_{y})=\Delta_{x}\Delta_{y}\sum_{i=1}^{n_{x}}\sum_{j=1}^{n_{y}}\cos(k_{x}\Delta_{x}(i-1))\cos(k_{y}\Delta_{y}(j-1))f(x,y),$
(90)

where $\Delta_{x}$ is the spacing between points $\textbf{x}_{i}$ and
$\textbf{x}_{i+1}$ and $\Delta_{y}$ is the spacing between points
$\textbf{x}_{i}$ and $\textbf{x}_{i+n_{x}}$ (here, $\Delta_{x}=\kappa h$ and
$\Delta_{y}=h$). The wavenumbers, $k_{x}$ and $k_{y}$, give the spatial
frequency of a pattern, with wavemodes $m_{x}$ and $m_{y}$ such that
$k_{x}=m_{x}\pi/L$ and $k_{y}=m_{y}\pi/L$, where $L$ is the side length of the
square domain $\Omega$ [33]. To showcase results, we plot the power spectrum,
defined as

$P_{s}=\left|\widehat{f}(k_{x},k_{y})\right|^{2},$ (91)

where $\widehat{f}(\cdot,\cdot)$ is the discrete Fourier cosine transform
defined in Equation (90).

We can predict which patterns are possible in equivalent PDE models of pattern
formation using linear stability analysis (LSA) [20] (see the Supplementary
Information for more details). Woolley et al. [32] developed a method using a
weak noise expansion (WNE) of the RDME [28] to determine the wavemodes that
can evolve in stochastic simulations (see the Supplementary Information for
more details). The advantage of the WNE approach over LSA in determining
whether a particular wavemode is possible for a given parameter set is that
the WNE takes into account the details of the diffusive jump rates. We compare
the simulation results to the predictions of the WNE and LSA to see whether
the patterns formed fall within the ranges of wavemodes predicted by those
methods.

All the simulations are run until $T=1800\,s$ on the domain
$\Omega=[0,1]\times[0,1]$, which is discretised into $40\times\lfloor
40\kappa\rfloor$ compartments. We use the following reaction rates:
$D_{u}=1\times 10^{-5}\,\mu m^{2}s^{-1}$, $D_{v}=1\times 10^{-3}\,\mu
m^{2}s^{-1}$, $k_{1}=1\,s^{-1}$, $k_{2}=0.02\,s^{-1}$, $k_{3}=1\times
10^{-6}\,s^{-1}$ and $k_{4}=3\,s^{-1}$. We initialise the system at the
spatially uniform steady state, placing 200 molecules of species U and 75
molecules of species V in each compartment.

##### The finite difference method

The left-most column of Figure 7 shows examples of patterns formed with FDM-
derived diffusive jump rates, while the second-from-left column of Figure 7
shows corresponding power spectra averaged over 100 simulations. There is good
agreement with the possible wavemodes predicted using the WNE (right-most
column of Figure 7) and LSA (second-from-right column of Figure 7). Our
results indicate that varying the parameter $\alpha$, which controls the rate
of diagonal jumping, has no discernible effect on pattern formation, except
when $\alpha=1.0$ (only diagonal jumps are possible) where checkerboard
patterns arise.

Figure 7: Results from simulations with FDM-derived diffusive jump rates with
$\kappa=1.0$. (a), (e), (i) An example pattern. (b), (f), (j) Averaged power
spectrum over 100 simulations. (c), (g), (k) Wavemodes predicted using LSA.
(d), (h), (l) Wavemodes predicted using the WNE.

##### The finite element method

The left-most plot in Figure 8 shows an example pattern formed with FEM-
derived diffusive jump rates, while the second-from-left plot shows the
corresponding power spectra averaged over 100 simulations. There is good
agreement with the possible wavemodes predicted using the WNE (right-most plot
of Figure 8) and LSA (second-from-right plot of Figure 8).

Figure 8: Results from simulations with FEM-derived diffusive jump rates with
$\kappa=1.0$. (a) An example pattern. (b) Averaged power spectrum over 100
simulations. (c) Wavemodes predicted using LSA. (d) Wavemodes predicted using
the WNE.

##### The finite volume method

For the case of FVM-derived diffusive jump rates, the simulations produced
very similar results to the FEM-derived diffusive jump rates. The averaged
power spectrum, as well as the predicted wavemodes using LSA and the WNE, are
very similar for the FVM and FEM cases, as can be seen in Figure 9.

Figure 9: Results from simulations with FVM-derived diffusive jump rates with
$\kappa=1.0$. (a) An example pattern. (b) Averaged power spectrum over 100
simulations. (c) Wavemodes predicted using LSA. (d) Wavemodes predicted using
the WNE.

##### The first exit time method

Lastly, we investigate the effects of using FET-derived diffusive jump rates
on pattern formation. Based on the averaged power spectrum shown in Figure
10(b), we see that changing the value of $\beta$ shifts the range of excited
wavemodes towards larger values. This observation is consistent with the
predictions made using the WNE (Figure 10(c)). In this respect, the
predictions of LSA fall down as they do not take into account details of the
diffusive jump rates (Figure 10(d)).

Figure 10: Results from simulations with FET-derived diffusive jump rates with
$\kappa=1.0$. (a), (e), (i) An example pattern. (b), (f), (j) Averaged power
spectrum over 100 simulations. (c), (g), (k) Wavemodes predicted using LSA.
(d), (h), (l) Wavemodes predicted using the WNE. The chequerboard pattern seen
in (a) is due to the jump rates being non-zero only in the diagonal directions
when $\beta=0$.

To explore whether changing the compartment aspect ratio has any bearing on
the patterns formed, we also simulated the system with $\kappa=1.4$. As in the
case of $\kappa=1.0$, larger wavemodes appear as the value of $\beta$
increases (Figure 11(b)). However, the range of wavemodes observed in the
$x$-direction is more significantly affected than the range in the
$y$-direction. This observation is, again, consistent with the predictions
made using the WNE, but cannot be predicted using LSA.

Figure 11: Results from simulations with FET-derived diffusive jump rates with
$\kappa=1.4$. (a), (e), (i) An example pattern. (b), (f), (j) Averaged power
spectrum over 100 simulations. (c), (g), (k) Wavemodes predicted using LSA.
(d), (h), (l) Wavemodes predicted using the WNE. Again, we see checkerboard
patterns for $\beta=0$ where there are no diagonal jumps.

In both cases, $\kappa=1.0$ and $\kappa=1.4$, changing the parameter $\beta$,
which controls the rate of molecules jumping diagonally in the FET derivation,
has a significant effect on the type of patterns formed, as shown in Figure 10
and Figure 11.

###  3.2 Growing domain

In this section, we explore the effects of using the different diffusive jump
rate derivations on a growing domain. As for the static domain case, we first
consider cases where we can compare with the analytic solution of the
corresponding macroscale PDE and we then consider models for pattern
formation. We assume domain growth to be isotropic, such that
$\Omega_{t}=[0,L(t)]^{2}$ and, for simplicity, we consider the domain to be
growing exponentially in time so that

$L(t)=L_{0}e^{rt},$ (92)

where $L_{0}$ is the initial domain size and $r$ the growth rate.

As we now have time-dependent propensity functions, we cannot use the
algorithm as stated in Section 2. We adapt Algorithm 1 using the Extrande
method [30]. As each of the propensity functions decreases with time, we use
the propensity function calculated at time $t$ as the upper bound in the time
interval, $[t,t+\delta t]$ in the Extrande method. Hence, the modified
algorithm for SSA is as stated in Algorithm 2.

Algorithm 2 Modified next sub-volume method SSA

Set $t=0$.

Initialise $\bm{u}$.

for $i$ in $\left\\{1,2,\ldots,I\right\\}$ do

Calculate the total propensity in compartment $i$, $a_{i}^{0}(\bm{u}(t),t)$.

Set the upper bound $\tilde{a}_{i}(\bm{u}(t),t)=a_{i}^{0}(\bm{u}(t),t)$

Initialise the time until the next reaction for compartment $i$:
$t_{i}\sim\exp(\tilde{a}_{i}(\bm{u}(t),t))$.

end for

while $t<T$ do

Find argmin of the set $\left\\{t_{1},t_{2},\ldots,t_{I}\right\\}$ and denote
it $m$.

Set $t=t_{m}$.

Generate a random number $r\sim\textit{U}(0,1)$.

Choose reaction / diffusion $q$ to fire. Note that reaction $k$ fires with
probability $a_{i}^{k}(\bm{u}(t),t)/\tilde{a}_{i}(\bm{u}(t),t)$, diffusion of
a molecule of species $\ell$ from box $i$ to box $j$ occurs with probability
$d_{i,j}^{\,\ell}(\bm{u}(t),t)/\tilde{a}_{i}(\bm{u}(t),t)$ and no reaction
happens with probability
$(\tilde{a}_{i}(\bm{u}(t),t)-a_{i}^{0}(\bm{u}(t),t))/\tilde{a}_{i}(\bm{u}(t),t)$.

Update molecule numbers: $\bm{u}:=\bm{u}+\bm{\nu}_{q}$, where $\bm{\nu}_{q}$
is the stoichiometric matrix of reaction / diffusion $q$.

Let $\Gamma$ be the set of indices of non-zero elements of $\bm{\nu}_{q}$.

for $\gamma$ in $\Gamma$ do

Update $a_{\gamma}^{k}(\bm{u}(t),t)$, $d_{\gamma,j}^{\,\ell}(\bm{u}(t),t)$,
$a_{\gamma}^{0}(\bm{u}(t),t)$ and $\tilde{a}_{i}(\bm{u}(t),t)$ accordingly.

Generate a random number $\Delta t\sim\exp(\tilde{a}_{\gamma}(\bm{u}(t),t))$.

Update the time until the next reaction: $t_{\gamma}=t+\Delta t$.

end for

end while

####  3.2.1 Diffusion

As in the static case, there is an analytic solution to the macroscale PDE
(see the Supplementary Information for details). Figure 12 shows how the
error, as defined in Section 3.1, evolves in time. We can see that in all the
cases shown in Figure 12, whether for a static or a growing case, the error
decreases with time and that the simulations with FVM-derived diffusive jump
rates have the lowest error.

Figure 12: Plots of the error, $e$ as a function of time for simulations with
diffusion coefficient $D=1\,\mu m^{2}s^{-1}$ and growth rate $r=0.1s^{-1}$. In
each plot we compare the errors from simulations for a static domain and a
growing domain where the jump rates were derived using: (a) the FDM
($\alpha=0.9$), (b) the FEM, and (c) the FVM. The domain, of initial size
$[0,5]\times[0,5]$, has been discretised into $21\times 21$ compartments and
all $5\times 10^{6}$ molecules are initially placed into the bottom-right
compartment.

Next, we explore the effects of changing the voxel aspect ratio, $\kappa$, on
the error, presenting results in Figure 13. The FVM-derived diffusive jump
rates provide the smallest error across the range of compartment aspect ratios
considered, increasing slightly as the aspect ratio increases. On the other
hand, for both the FEM- and FVM-derived diffusive jump rates, there is a
slight decrease in the error as the compartment aspect ratio is increased.

Figure 13: Plots of the error, $e$, as a function of voxel aspect ratio,
$\kappa$, for simulations with diffusion coefficient $D=1\,\mu m^{2}s^{-1}$
and growth rate $r=0.1\,s^{-1}$. We compare different derivations of the
diffusive jump rates as we vary $\kappa$ For the FDM, $\alpha=0.7$. The
domain, of initial size $[0,5]\times[0,5]$, is discretised into $21\times
21\kappa$ compartments and all $5\times 10^{6}$ molecules are initially placed
into the bottom-right compartment.

####  3.2.2 Schnakenberg kinetics

We now consider pattern formation on a growing domain to determine to what
extent the different derivations of diffusive jump rates affect patterns
formed. We again use the Schnakenberg system, as in Section 3.1.4. Figure 14
showcases the range of patterns that form as the domain growth rate is varied.

Figure 14: Comparison of different growth rates for simulations of the
Schnakenberg system with the following parameters: $D_{u}=1\times 10^{-5}\,\mu
m^{2}s^{-1}$, $D_{v}=0.001\,\mu m^{2}s^{-1}$, $k_{1}=1.0\,\mu m^{-2}s^{-1}$,
$k_{2}=0.02\,s^{-1}$, $k_{3}=1\times 10^{-6}\,\mu m^{4}s^{-1}$,
$k_{4}=3.0\,\mu m^{-2}s^{-1}$ run until $T=1800\,s$ and FVM-derived jump rates
used. The domain has been discretised into $40\times 40$ compartments.

We see in Figure 15 that the patterns formed are different compared to the
static case. However, using different derivations for the diffusive jump rates
does not seem to markedly affect the patterns produced. The reason we see
different patterns when the domain grows is because higher wavemodes become
available, and the inherent stochasticity in the system causes clusters of
molecules (peaks) to split and these higher wavemodes are then realised. This
change in the wavemode is highlighted in the centre and left-hand columns in
Figure 15.

Figure 15: Results from simulations of the Schnakenberg system on a growing
domain with the following parameters: $D_{u}=1.0\times 10^{-5}\,\mu
m^{2}s^{-1}$, $D_{v}=0.001\,\mu m^{2}s^{-1}$, $k_{1}=1.0\,\mu m^{-2}s^{-1}$,
$k_{2}=0.02\,s^{-1}$, $k_{3}=1\times 10^{-6}\,\mu m^{4}s^{-1}$,
$k_{4}=3.0\,\mu m^{-2}s^{-1}$ run until $T=1800\,s$. The domain has been
discretised into $40\times 40$ compartments. Left-hand column – example
patterns; centre column – averaged power spectrum over 100 simulations; right-
hand column – predictions from LSA.

##  4 Conclusions

Our aim in this work was to investigate the effects of different methods of
derivation of diffusive jump rates on stochastic reaction-diffusion systems
modelled using the RDME. To do so, we replicated, extended and applied the
results of Meinecke and Lötstedt [18] to a number of reaction-diffusion
systems, that included a range of reaction types, on both static and growing
domains.

For the production-decay examples, which consisted solely of zeroth and first
order reactions, the FVM generally gives the most reliable results in terms of
the error between stochastic simulations and solution of the corresponding
macroscale PDE (as does the FDM for $\alpha=0$ where there are no diagonal
jumps and we have the same diffusive jump rates as for the FVM) (Figure 5).
When bimolecular reactions are included, we showed that the spatial
discretisation affects results, for example shifting the stationary
distribution (Figure 6).

We then moved to study the effects of different methods of derivation of
diffusive jump rates on pattern formation, using the Turing reaction-diffusion
model with Schnakenberg kinetics. We compared the averaged power spectra from
simulations with the wavemodes predicted using LSA, which analyses the
corresponding PDE model, and the WNE, which takes into account the details of
the diffusive jump rates (Figures 7–10). Our results show that the WNE
approach provides good predictions of the range of possible wave modes, as
expected from the work of Woolley et al. [32]. In particular, we showed that
the aspect ratio of the compartments can significantly affect the patterns
formed (see Figure 11).

We also considered the effects of uniform domain growth on derivations of the
diffusive jump rates, and simulations of pattern formation. We considered the
FDM, FVM and FEM showing, in line with the results of Woolley et al. [32],
that the diffusive jump rates become time-dependent as the compartment
dimensions now change in time. It is not possible to use the FET approach in
the context of growing domains because the mean of the FET distribution
becomes infinite, and it is difficult to approximate the FET distribution with
an exponential distribution.

Taken together, the results presented in this work demonstrate that care
should be taken when choosing how to model diffusion in the context of RDME
models of stochastic reaction-diffusion systems. Future work will consider how
to extend these results to unstructured meshes and more complicated forms of
domain growth.

Acknowledgements

R.E.B. is a Royal Society Wolfson Research Merit Award holder, would like to
thank the Leverhulme Trust for a Research Fellowship. B.J.B. would like to
thank the EPSRC for supporting this research through grant EP/G03706X/1.

## References

  * [1] D. F. Anderson.  A modified next reaction method for simulating chemical systems with time dependent propensities and delays.  The Journal of Chemical Physics, 127(21):214107, 2007.
  * [2] J. Elf and M. Ehrenberg.  Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases.  Systems Biology, 1(2):230–236, 2004.
  * [3] H. C. Elman, D. J. Silvester, and A. J. Wathen.  Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics.  Oxford University Press, 2010.
  * [4] M. B. Elowitz, A. J. Levine, E. D. Siggia, et al.  Stochastic gene expression in a single cell.  Science, 297(5584):1183–1186, 2002.
  * [5] S. Engblom, L. Ferm, A. Hellander, and P. Lötstedt.  Simulation of stochastic reaction-diffusion processes on unstructured meshes.  SIAM Journal on Scientific Computing, 30(6):2709–2733, 2008.
  * [6] R. Erban, S. J. Chapman, and P. K. Maini.  A practical guide to stochastic simulations of reaction-diffusion processes.  arXiv preprint, 2007.
  * [7] A. Gierer and H. Meinhardt.  A theory of biological pattern formation.  Kybernetik, 12(1):30–39, 1972.
  * [8] D. T. Gillespie.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  Journal of Computational Physics, 22(4):403–434, 1976.
  * [9] J. Hattne, D. Fange, and J. Elf.  Stochastic reaction-diffusion simulation with MesoRD.  Bioinformatics, 21(12):2923–2924, 2005.
  * [10] S. Hellander and L. Petzold.  Reaction rates for reaction-diffusion kinetics on unstructured meshes.  Journal of Chemical Physics, 146(6), 2017.
  * [11] M. Howard and A. D. Rutenberg.  Pattern formation inside bacteria: fluctuations due to the low copy number of proteins.  Physical Review Letters, 90(12):128102, 2003.
  * [12] S. A. Isaacson.  The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target.  SIAM Journal on Applied Mathematics, 70(1):77–111, 01 2009.
  * [13] S. A. Isaacson.  A convergent reaction-diffusion master equation.  Journal of Chemical Physics, 139(5), 2013.
  * [14] S. A. Isaacson, D. M. McQueen, and C. S. Peskin.  The influence of volume exclusion by chromatin on the time required to find specific DNA binding sites by diffusion.  Proceedings of the National Academy of Sciences, 2011.
  * [15] D. Longo and J. Hasty.  Dynamics of single-cell gene expression.  Molecular Systems Biology, 2(1):64, 2006.
  * [16] P. Lötstedt and L. Meinecke.  Simulation of stochastic diffusion via first exit times.  Journal of Computational Physics, 300:862–886, 2015.
  * [17] L. Meinecke, S. Engblom, A. Hellander, and P. Lötstedt.  Analysis and design of jump coefficients in discrete stochastic diffusion models.  SIAM Journal on Scientific Computing, 38(1):A55–A83, 2016.
  * [18] L. Meinecke and P. Lötstedt.  Stochastic diffusion processes on Cartesian meshes.  Journal of Computational and Applied Mathematics, 294:1–11, 2016\.
  * [19] R. L. Mort, R. J. H. Ross, K. J. Hainey, O. J. Harrison, M. A. Keighren, G. Landini, R. E. Baker, K. J. Painter, I. J. Jackson, and C. A. Yates.  Reconciling diverse mammalian pigmentation patterns with a fundamental mathematical model.  Nature Communications, 7(1):10288, 2016.
  * [20] J. D. Murray.  Mathematical Biology II: Spatial Models and Biomedical Applications.  Springer-Verlag New York Incorporated, 2001.
  * [21] J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini, and D. J. Gavaghan.  Comparing individual-based approaches to modelling the self-organization of multicellular tissues.  PLoS Computational Biology, 13(2):1–34, 2017.
  * [22] S. Redner and J. R. Dorfman.  A guide to first-passage processes.  American Journal of Physics, 70(11):1166–1166, 2002.
  * [23] J. Schnakenberg.  Simple chemical reaction systems with limit cycle behaviour.  Journal of Theoretical Biology, 81(3):389–400, 1979.
  * [24] D. Schnoerr, G. Sanguinetti, and R. Grima.  Approximation and inference methods for stochastic biochemical kinetics—a tutorial review.  Journal of Physics A: Mathematical and Theoretical, 50(9):093001, 2017.
  * [25] M. J. Simpson and R. E. Baker.  Exact calculations of survival probability for diffusion on growing lines, disks, and spheres: The role of dimension.  Journal of Chemical Physics, 143(9):94109, 2015.
  * [26] J. C. Strikwerda.  Finite Difference Schemes and Partial Differential Equations, volume 88.  SIAM, 2004.
  * [27] A. M. Turing.  The chemical basis of morphogenesis.  Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 237(641):37–72, 1952.
  * [28] N. G. van Kampen.  Stochastic Processes in Physics and Chemistry.  North Holland Personal Library. Elsevier, Oxford, 3rd edition, 2007.
  * [29] P. Van Liedekerke, M. M. Palm, N. Jagiella, and D. Drasdo.  Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results.  Computational Particle Mechanics, 2(4):401–444, 2015.
  * [30] M. Voliotis, P. Thomas, R. Grima, and C. G. Bowsher.  Stochastic Simulation of Biomolecular Networks in Dynamic Environments.  PLoS Computational Biology, 2016.
  * [31] A. Volkening and B. Sandstede.  Modelling stripe formation in zebrafish: an agent-based approach.  Journal of The Royal Society Interface, 12(112):20150812, 2015.
  * [32] T. E. Woolley, R. E. Baker, E. A. Gaffney, and P. K. Maini.  Power spectra methods for a stochastic description of diffusion on deterministically growing domains.  Physical Review E, 84(2), 2011.
  * [33] T. E. Woolley, R. E. Baker, E. A. Gaffney, and P. K. Maini.  Stochastic reaction and diffusion on growing domains: understanding the breakdown of robust pattern formation.  Physical Review E, 84(4):46216, 2011.


# Spin triplet exciton condensations in LaCoO3 at ultrahigh magnetic fields up
to 600 T

Akihiko Ikeda<EMAIL_ADDRESS>Institute for Solid State Physics, University
of Tokyo, Kashiwa, Chiba 277-8581, Japan  Department of Engineering Science,
University of Electro-Communications, Chofu, Tokyo 182-8585, Japan  Yasuhiro
H. Matsuda  Institute for Solid State Physics, University of Tokyo, Kashiwa,
Chiba 277-8581, Japan  Keisuke Sato  National Institute of Technology, Ibaraki
College, Hitachinaka, Ibaraki 312-0011, Japan  Yuto Ishii  Hironobu Sawabe
Daisuke Nakamura  Shojiro Takeyama  Institute for Solid State Physics,
University of Tokyo, Kashiwa, Chiba 277-8581, Japan  Joji Nasu  Department of
Physics, Tohoku University, Sendai, Miyagi 980-8578, Japan  PRESTO, Japan
Science and Technology Agency, Honcho Kawaguchi, Saitama 332-0012, Japan

Figure 1: Schematic drawings of (a) crystal structure, (b) spin-state degrees
of freedom, (c) and corresponding excitonic degrees of freedom in LaCoO3. LS
and a pair of shaded electron and hole in (b) correspond to vacuum and an
exciton in (c), respectively.

Bose-Einstein condensation of electron hole pairs, exciton condensation, has
been effortfully investigated since predicted 60 years ago [1, 2]. Irrefutable
evidence has still been lucking due to experimental difficulties in verifying
the condensation of the charge neutral and non-magnetic spin-singlet excitons
[3, 4, 5, 6, 7]. Whilst, condensation of spin-triplet excitons is a promising
frontier because the spin supercurrent [8, 9], spin-Seebeck effects [10] will
be observable. A canonical cobaltite LaCoO3 under very high magnetic fields is
a propitious candidate, yet to be verified [11, 12, 13, 14, 15]. Here, we
unveil the exotic phase diagram of LaCoO3 up to 600 T generated using the
electromagnetic flux compression method [16] and the state-of-the-art
magnetostriction gauge [17]. We found successive emergence of two spin-triplet
exciton condensates identified by the gapless magnetostriction curves and
model calculations. The spin-triplet exciton condensation in a cobaltite, that
is three dimensional and thermally equilibrated, opens up a novel venue for
spintronics technologies with spin-supercurrent such as a spin Josephson
junction.

Figure 2: Schematics drawings of (a) the interactions and dynamics of the
excitons, and (b) the magnetic field effect on the exciton states in LaCoO3.
(c) Expected emergence of the magnetic field induced phases in LaCoO3, and the
resultant magnetostriction curves.

Novel phases in condensed matters often emerge at the boundaries of
incompatible phases aided by the competing interactions, fluctuations and
spontaneous symmetry breakings. Insulating state by a Bose-Einstein
condensation of electron-hole pairs, excitonic insulator or exciton
condensation, has been eagerly investigated for over 50 years since the
prediction [2] at the boundary region between semi-metals and semi-conductors
and their photo-excited states. However, compelling evidence has been lucking,
due to the difficulty in verifying the condensation of charge neutral and
spin-singlet excitons [3, 4, 5, 6, 7].

Figure 3:  Magnetostriction data of LaCoO3 (a) as a function of time at
various initial temperatures along with the magnetic field profile generated
using the electromagnetic flux compression method, and (b) those as a function
of magnetic field at various temperatures, along with the previous data (Green
curves) using STC in Ref. [18]. (c) Summary of the magnetostriction
measurement of LaCoO3 on a magnetic field-temperature plane up to 600 T along
with the previously obtained data using magnetization [19] and
magnetostriction measurements [18].

On the other hand, the condensation of the spin-triplet excitons is more
beneficial to investigate thanks to observable magnetic phenomena such as
spin-supercurrent [8, 9] and spin Seebeck effect [10], which are attractive
for potential applications in spintronics and quantum computing technologies.
Condensation of spin-triplet bosons is of a fundamental interest in analogy
with the spin-triplet superconductivity [20, 21] and superfluidity of 3He
[22]. The condensation of spin-triplet excitons has been rarely considered in
conventional studies of exciton condensations [3, 4, 5, 6, 7].

Recently, theories have predicted that spin-triplet exciton condensation is
stable in perovskite cobaltites with strongly correlated electrons [11, 12,
13]. This is based on the fact that the spin-state degrees of freedom are
inherent in the cobaltites [see Fig. 1(a) and 1(b)] due to the competing
Hund’s coupling and the crystal field splitting, which can also be viewed as
the novel degrees of freedom of atomic size excitons [see Fig. 1(c)]. The
spin-triplet exciton condensation is essentially the quantum hybridization of
the distinct exciton states. The hybridization is due to the delocalization,
interaction, and duality of excitons [see Fig. 2(a)] that are promoted when
the energies of the vacuum and exciton states are in proximity by such
external fields as magnetic fields. The LS and HS states correspond to the
Mott insulator and the band insulator phases, respectively [11]. Thereby, the
spin-triplet exciton condensation occurs at the boundary region between the
incompatible phases [See Fig. 2(b)]. Note that it is also a multipole order of
spin and orbital degrees of freedom [23].

A doped perovskite cobaltite Pr0.5Ca0.5CoO3 and its family compounds are the
first candidate for the spin-triplet excitonic condensation, exhibiting an
insulating and paramagnetic ground phase at low temperatures [24, 25]. The
origin of the insulator phase is, however, not completely uncovered because
the significant valence transition of Pr complicates the phase transition.
Next, the most well-known cobaltite LaCoO3 has been claimed to be a promising
candidate when placed under very high magnetic fields exceeding 100 T [14,
15]. LaCoO3 has a variety of spin-states, such as low-spin (LS: $S=0$,
$t_{2g}^{6}e_{g}^{0}$), intermediate spin (IS: $S=1$, $t_{2g}^{5}e_{g}^{1}$)
and high spin (IS: $S=2$, $t_{2g}^{4}e_{g}^{2}$) states [See Fig. 1(b)], which
are viewed as exciton states of vacuum, single exciton and bi-exciton,
respectively [See Fig. 1(c)]. The spin-triplet excitonic condensation will
emerge when the magnetic field eliminates the spin gap and changes the vacuum
state (LS) to a magnetic single exciton state (IS) or bi-exciton states (HS)
[See Fig. 2(b)].

Figure 4: The calculated results of the proposed model of the excitonic states
of $\rm{LaCoO_{3}}$ under ultrahigh magnetic field. The parameters are chosen
to be
$(E_{1},E_{2},J_{11},J_{02},t,V,g_{1},g_{2})=(1.0,0.7,0.2,0.05,0.07,0.17,1,2)$.
See main text and the Supplemental Material for more details.

However, the experimental investigation has been challenging due to the need
for very high magnetic fields beyond 100 T necessitated by the large spin-gap
of $\sim 100$ K, and also the need for the probe of the exciton states at such
high fields. Recently, we reported the generation of 1200 T with
electromagnetic flux compression (EMFC) method using a newly installed
discharge capacitor banks in ISSP [16]. Besides, we succeeded in implementing
a state-of-the-art high-speed 100 MHz strain gauge using fiber Bragg grating
(FBG) and optical filter method [17], which have enabled us to measure
magnetostriction in the $\mu$s-pulsed 1000 T environment. Magnetostriction is
a crucial probe of the exciton states because the density of excitons is
coupled to the lattice volume of LaCoO3. As indicated in Refs. [14, 15, 26],
the plateaux and slopes in the magnetostriction correspond to solidifications
and Bose-Einstein condensation of excitons in LaCoO3, respectively [See Fig.
2(c)]. Note that the correspondences are analogues to the magnetization
plateaux and slopes in magnon solids and superfluids, respectively [27, 28].

Previously in LaCoO3, the field induced phase transitions are uncovered below
30 K ($\alpha$ phase) [29, 30, 31] and above 30 K ($\beta$ phase), which
extend beyond 100 T, using magnetization and magnetostriction measurements up
to 200 T generated by single turn coil (STC) method [19, 18]. The two phases
are considered superlattice formation of the single excitons, bi-excitons, and
vacuum based on the plateaux of magnetostrictions. Considering that the spin-
state and magnetization below 200 T being far from the full polarization,
further orderings of excitons are expected to appear at above 200 T. Now, we
are fully equipped for the investigation of the evolution of the excitonic
states in LaCoO3 beyond 200 T.

Figs. 3(a)-3(c) and shows the results of longitudinal magnetostriction
measurements of polycrystalline LaCoO3 up to 600 T with the initial sample
temperatures $T_{\rm{ini}}=5$, 78, and 108 K. Fig. 3(a) and 3(b) shows the
magnetostriction data as functions of time and magnetic fields, respectively.
The magnetostriction data up to 150 T in Ref. [18] are shown for comparison
using green curves in Fig. 3(b). The transition points are indicated in 3(b)
and summarized using star shaped symbols in Fig. 3(c) along with the previous
results up to 150 T [19, 18] shown using filled and open triangle and squares.
The transitions and features below 200 T well reproduce the previous
observations up to 150 T [See Fig. 3(b)] where the transition fields are
denoted as $B_{\rm{C1}}$ and $B_{\rm{C2}}$. In Fig. 3(a), the propagation of a
shock wave inside the sample initiated by the strong transitions of
$B_{\rm{C1}}$ and $B_{\rm{C2}}$ are present [18, 32, 33]. The temperature
changes of the sample during the application of the $\mu$-second pulsed
magnetic field are considered significant at $T_{\rm{ini}}=5$ K, while it is
negligible at $T_{\rm{ini}}>30$ K based on the sufficiently large heat
capacity above 30 K [34] and a preliminary measurement of magnetocaloric
effects up to 65 T [35].

We found that the magnetostriction plateau of $\alpha$ phase at
$T_{\rm{ini}}=5$ K that appears at 70 T persists up to 600 T, revealing the
significant stability under magnetic fields. On the other hand, $\beta$ phase
has been found to show a negatve slope of magnetostriction [See slope 1 in
Fig. 3(b)] between 170 T and 380 T. Further, we found a new phase transition
at 380 T denoted as $B_{\rm{C3}}$ from $\beta$ phase to a novel state in the
data at $T_{\rm{ini}}=78$ K, which we term $\gamma$ phase, being characterized
with a positive slope beyond 380 T [See slope 2 in Fig. 3(b)]. Previously in
Ref. [18], $\beta$ phase was falsely considered a plateau due to the limited
measurement regions below 200 T. We note that the sharpness of the transition
at $B_{\rm{C3}}$ is lost with increasing temperature at $T_{\rm{ini}}=108$ K
data, implying that the upper boundary of $\gamma$ phase is located at $\sim
100$ K.

Here, we discuss the origins of $\alpha$, $\beta$, and $\gamma$ phases in
LaCoO3 induced at ultrahigh magnetic fields. We claim that $\alpha$ phase is a
superlattice of spin-triplet excitons and vacuum, based on the present
observation of the magnetostriction plateau at 5 K shown in Fig. 3(a) and that
the magnetization of $\alpha$ phase is only 1/8 or 1/4 of the saturation
excitons [29, 30]. In the exciton superlattice, the number of exciton is
discrete [13], leading to the magnetostriction plateau as is shown in the data
at 5 K in Fig. 3(a). This is in good agreement with the previous observations
and considerations on $\alpha$ phase [29, 30, 31, 18]. The remarkable
stability of $\alpha$ phase that sustains up to 600 T suggests the localized
nature of the spin-triplet excitions in $\alpha$ phase of LaCoO3.

On the other hand, $\beta$ and $\gamma$ phases show gapless magnetostriction
slopes denoted as slope 1 and slope 2 in Figs. 3(d), which indicate that
exciton density varies continuously as a function of magnetic fields. These
gapless behaviors evoke the Bose-Einstein condensation of spin-triplet
excitons. In a condensate of excitons, the single site wave function of Co is
described by $\ket{\phi_{\rm{EC}}}=a\ket{0}+b\ket{1}+c\ket{2}$, where
$\ket{0}$, $\ket{1}$, $\ket{2}$ represents the vacuum, a single-exciton and a
bi-exciton occupied states, respectively, with a global wave-function of
$\ket{\psi_{\rm{EC}}}=\prod_{i}\ket{\phi_{\rm{EC}}}_{i}$. The continuous
change of magnetostriction and magnetization is realized because the
coefficients $a$, $b$, $c$ change continuously in exciton condensates as a
function of magnetic field [15, 14]. A possible origin of the continuous
magetostriction is thermal activation, which is ruled out by the negative
magnetostriction slopes as shown in Fig. 3(d). Another possibility, orbital
order, is also ruled out as the main cause of the present result because
polycrystalline samples are used and LaCoO3 is a mostly isotoropic pseudo-
cubic system [36].

For considering the condensation of spin-triple exciton condensation in LaCoO3
at ultrahigh magnetic fields, we tentatively construct a phenomenological
model and analyze it using the mean-field theory with the $A$, $B$
sublattices. In the presence of strong magnetic field, the model Hamiltonian
is constructed as,

$\mathcal{H}=\mathcal{H}_{\rm{loc}}+\mathcal{H}_{\rm{int}}+\mathcal{H}_{\rm{trans}}+\mathcal{H}_{\rm{dual}}.$
(1)

Each terms are expressed as

$\displaystyle\mathcal{H}_{\rm{loc}}=$
$\displaystyle\left(E_{\rm{2}}-g_{2}h\right)\sum_{i}n_{i}^{\rm{(2)}}+\left(E_{\rm{1}}-g_{\rm{1}}h\right)\sum_{i}n_{i}^{\rm{(1)}},$
(2) $\displaystyle\mathcal{H}_{\rm{int}}=$ $\displaystyle-
J_{\rm{11}}\sum_{\braket{ij}}n_{i}^{\rm{(1)}}n_{j}^{\rm{(1)}}-J_{\rm{02}}\sum_{\braket{ij}}\big{(}n_{i}^{\rm{(2)}}n_{j}^{\rm{(0)}}+n_{i}^{\rm{(0)}}n_{j}^{\rm{(2)}}\big{)},$
(3) $\displaystyle\mathcal{H}_{\rm{trans}}=$
$\displaystyle-t\sum_{\braket{ij}}\left(\ket{0,1}_{ij}\bra{1,0}_{ij}+\ket{1,0}_{ij}\bra{0,1}_{ij}\right),$
(4) $\displaystyle\mathcal{H}_{\rm{dual}}=$
$\displaystyle-V\sum_{\braket{ij}}\bigg{(}\ket{1,1}_{ij}\left(\bra{2,0}_{ij}+\bra{0,2}_{ij}\right)$
(5)
$\displaystyle+\left(\ket{2,0}_{ij}+\ket{0,2}_{ij}\right)\bra{1,1}_{ij}\bigg{)},$
(6)

where, for simplicity, only three states {$\ket{0}$, $\ket{1}$, $\ket{2}$} are
considered as the local bases [37]. The local term $\mathcal{H}_{\rm{loc}}$
defines the energy of each exciton state at magnetic fields $h$. $E_{k}$,
$g_{k}$, and $n^{(k)}_{i}$ with $k=0,1,2$ are the energy, effective $g$-factor
and number operator of $\ket{k}_{i}$ excitonic states. The interaction term
$\mathcal{H}_{\rm{int}}$ describes the attractive interactions $J_{11}$,
$J_{02}$ between neighboring excitons as $\ket{1,1}_{ij}$ and
$\ket{0,2}_{ij}$, respectively, where $\braket{ij}$ represents the nearest
neighbor pair. The transfer term $\mathcal{H}_{\rm{trans}}$ describes the
delocalization of single exciton states as
$\ket{0,1}_{ij}\rightleftharpoons\ket{1,0}_{ij}$. The duality term
$\mathcal{H}_{\rm{dual}}$ hybridizes the neighboring pair $\ket{2,0}_{ij}$ and
the pair $\ket{1,1}_{ij}$ as $\ket{2,0}_{ij}\rightleftharpoons\ket{1,1}_{ij}$,
which describes the duality of a $\ket{2}$ state and two $\ket{1}$ states [38,
39]. The order parameters of the excitonic condensations are defined as
$\braket{\bm{\tau}}$ and $\braket{\bm{\rho}}$, where
$\tau_{i}^{x}=\ket{0}_{i}\bra{1}_{i}+{\rm H.c.}$,
$\tau_{i}^{y}=i\ket{0}_{i}\bra{1}_{i}+{\rm H.c.}$,
$\rho_{i}^{x}=\ket{1}_{i}\bra{2}_{i}+{\rm H.c.}$, and
$\rho_{i}^{y}=i\ket{1}_{i}\bra{2}_{i}+{\rm H.c.}$, respectively. We further
define
$\bar{n}_{k}^{\pm}=\left|\braket{n^{(k)}}_{A}\pm\braket{n^{(k)}}_{B}\right|$,
$\bar{\tau}^{\pm}=\sqrt{\left(\braket{\tau^{x}}_{A}\pm\braket{\tau^{x}}_{B}\right)^{2}+\left(\braket{\tau^{y}}_{A}\pm\braket{\tau^{y}}_{B}\right)^{2}}$,
and
$\bar{\rho}^{\pm}=\sqrt{\left(\braket{\rho^{x}}_{A}\pm\braket{\rho^{x}}_{B}\right)^{2}+\left(\braket{\rho^{y}}_{A}\pm\braket{\rho^{y}}_{B}\right)^{2}}$,
with $\langle\cdots\rangle_{A(B)}$ being the expectation value at the $A$
($B$) sublattice. We note that nonzero $\bar{\tau}^{\pm}$, $\bar{\rho}^{\pm}$
indicate the hybridizations of
$\ket{0,1}_{ij}\rightleftharpoons\ket{1,0}_{ij}$, and
$\ket{2,0}_{ij}\rightleftharpoons\ket{1,1}_{ij}$, respectively, that are
excitonic condensations. Note also that nonzero $\bar{n}^{-}_{k}$,
$\bar{\tau}^{-}$, $\bar{\rho}^{-}$ indicate a staggered order.

The calculated results are summarized in Fig. 4. Two kinds of exciton
condensates actually appear, indicated by the quantum hybridization of the
single-site exciton states. From $h=0$ to $h_{\rm{C1}}$, the spatially uniform
$\psi_{\rm{Vacuum}}=\prod_{i}\ket{0}_{i}$ state is the ground state. With a
first order transition, the first exciton condensate (EC1) emerges between
$h_{\rm{C1}}$ and $h_{\rm{C2}}$, where the magnetic excitonic condensation of
the type $\phi_{\rm{EC1}}=p_{0}\ket{0}+p_{1}\ket{1}+p_{2}\ket{2}$ appears as
indicated by $n_{k}^{+}>0$ for any $k$. EC1 is a spatially uniform
condensation of single exciton $\ket{1}$ and also of bi-exiciton $\ket{2}$,
which is indicated by $n_{k}^{-}=0$ for any $k$. The emergence of EC1 is
assisted by the terms $\mathcal{H}_{\rm{trans}}$ and $\mathcal{H}_{\rm{dual}}$
that induces quantum hybridizations of $\ket{0}$, $\ket{1}$, and $\ket{2}$
states. Next, the second exciton condensate (EC2) emerges between
$h_{\rm{C2}}$ and $h_{\rm{C3}}$ with a distinct character from EC1, where the
$A$ and $B$ sublattices hosts different excitonic states of the type
$\phi^{A}_{\rm{EC2}}\simeq q_{0}\ket{0}+q_{1}\ket{1}$ and
$\phi^{B}_{\rm{EC2}}\simeq q^{\prime}_{2}\ket{2}+q^{\prime}_{1}\ket{1}$,
respectively. As shown in Fig. 4, $n_{k}^{-}$ for $k=0$ and 2 take large
values, suggesting the spatially periodic arrangement of $\ket{0}$ and
$\ket{2}$. The superlattice formation coexists with the excitonic condensation
in EC2, which thus corresponds to a supersolidity of excitons. The quantum
terms $\mathcal{H}_{\rm{trans}}$ and $\mathcal{H}_{\rm{dual}}$ are
indispensable for the appearance of EC1 and EC2, respectively, which are
separated by a second-order phase transition. Lastly, above $h_{\rm{C3}}$, the
saturated phase of $\psi_{\rm{Bi-exciton}}=\prod_{i}\ket{2}_{i}$ appears due
to the largest magnetization and Zeeman energy.

With only the classical terms $\mathcal{H}_{\rm{loc}}$ and
$\mathcal{H}_{\rm{int}}$, the phase transitions occur as
$\psi_{\rm{Vacuum}}\leftrightarrow$ superlattice $\ket{0,2,0,2\cdots}$
$\leftrightarrow\psi_{\rm{Bi-exciton}}$ with increasing magnetic fields, as is
discussed in Ref. [30]. With the additional quantum effects
$\mathcal{H}_{\rm{trans}}$ and $\mathcal{H}_{\rm{dual}}$, the nontrivial
excitonic phases EC1 and EC2 emerge over the classical phases. In fact, EC2
has a root in the classical solid $\ket{0,2,0,2\cdots}$ and the quantum
hybridization $\mathcal{H}_{\rm{dual}}$. When $t\gg V$, only EC1 appears, and
when $t\ll V$, only EC2 appears, on the contrary. When $t$ and $V$ balance,
both EC1 and EC2 emerge as a result of the competing energy. We note that the
experimentally observed transitions from LS $\leftrightarrow$ $\beta$
$\leftrightarrow$ $\gamma$ may correspond to the successive phase transitions
in the calculated results, $\psi_{\rm{Vacuum}}\leftrightarrow$ EC1
$\leftrightarrow$ EC2, considering the features that two gapless phases appear
in high magnetic fields and that those phases are connected with a continuous
phase transition. The correspondence of the two kinds of gapless phases
$\beta$ and $\gamma$ to EC1 and EC2 infers that considering the duality of the
excitons $\mathcal{H}_{\rm{dual}}$ and the situation $t\sim V$ may actually be
important in LaCoO3. See the caption of Fig. 4 and the Supplemental Material
for the details of the parameters in the calculation [37]. The negative slope
in $\beta$ may result from the non-trivial spin-lattice coupling found in a
high-resolution x-ray diffraction study, where the lattice volume of $\ket{1}$
is claimed to be smaller than a pair of neighboring site $\ket{0}$ and
$\ket{2}$ [40]. Increasing the weight of $\ket{1}$ might contract the lattice
volume in EC1 with increasing magnetic fields. While, increasing the weight of
the alignment with $\ket{0}$ and $\ket{2}$ might expand the lattice volume in
EC2. No phase corresponding to $\alpha$ appear in the calculated result. We
speculate that $\alpha$ phase emerges with modified interaction parameters by
the strong lattice contraction as compared to $\beta$ phase where the
localization nature of excitons is enhanced resulting in exciton superlattice
formation in $\alpha$ phase. Such lattice contraction in the course of
$\beta$-$\alpha$ phase transition is actually observed in Ref. [18].

We comment on the spin-crossover physics in LaCoO3, which has been a long-
standing controversy over half a century regarding the temperature evolution
of the microscopic understanding of the electronic and magnetic state in terms
of the spin-states [41]. As depicted in Fig. 1, there is a correspondence of
the exciton occupation states {$\ket{0}$ , $\ket{1}$, $\ket{2}$} and spin-
states {LS, IS, HS}. HS-LS order have been anticipated by many studies in
LaCoO3 [42, 30]. EC2, that may correspond to $\beta$ phase, possesses the
component of such HS-LS order in the present experiment. Besides, $\alpha$
phase is claimed to originate in a IS-LS superlattice [18]. Recently, the HS-
IS duality has been claimed to be at the heart of the controversy over HS and
IS states in LaCoO3 [38, 39], in addition to the delocalization of IS. The
delocalization of IS to LS sites, and the HS-IS duality correspond to the
itineracy of single exciton and the exciton duality
$\ket{2,0}_{ij}\rightleftharpoons\ket{1,1}_{ij}$, respectively, that give rise
to the exciton condensation EC1 and EC2, respectively. Thus, we argue that the
gapless features of $\beta$ and $\gamma$ phases may originate in the coherent
hybridization of LS-IS-HS states due to the HS-IS duality. Incoherent spin-
state hybridizations and superlattices are also argued in Ref. [43]. With an
influence of the hybridization effects of LS, IS and HS states, the variety of
ordered states of LaCoO3 at high magnetic fields are realized as shown in Fig.
3(c) and Fig. 4. To fully understand and characterize the possible exciton
condensates in $\beta$, and $\gamma$ phases of LaCoO3, experimental
investigations such as spin transport, neutron scattering and theoretical
examinations of more realistic models are mandatory.

## References

  * Mott [1961] N. F. Mott, The transition to the metallic state, Philos. Mag. 6, 287 (1961).
  * Halperin and Rice [1968] B. I. Halperin and T. M. Rice, Possible anomalies at a semimetal-semiconductor transistion, Rev. Mod. Phys. 40, 755 (1968).
  * Kogar  _et al._ [2017] A. Kogar, M. S. Rak, S. Vig, A. A. Husain, F. Flicker, Y. I. Joe, L. Venema, G. J. MacDougall, T. C. Chiang, E. Fradkin, J. van Wezel, and P. Abbamonte, Signatures of exciton condensation in a transition metal dichalcogenide, Science 358, 1314 (2017).
  * Lu  _et al._ [2017] Y. F. Lu, H. Kono, T. I. Larkin, A. W. Rost, T. Takayama, A. V. Boris, B. Keimer, and H. Takagi, Zero-gap semiconductor to excitonic insulator transition in Ta2NiSe5, Nat. Commun. , 14408 (2017).
  * Lai  _et al._ [2007] C. W. Lai, N. Y. Kim, S. Utsunomiya, G. Roumpos, H. Deng, M. D. Fraser, T. Byrnes, P. Recher, N. Kumada, T. Fujisawa, and Y. Yamamoto, Coherent zero-state and $\pi$-state in an exciton-polariton condensate array, Nature 450, 529 (2007).
  * Jia  _et al._ [2021] Y. Jia, P. Wang, C.-L. Chiu, Z. Song, G. Yu, B. J’́ack, S. Lei, S. Klemenz, F. A. Cevallos, M. Onyszczak, N. Fishchenko, X. Liu, G. Farahi, F. Xie, Y. Xu, K. Watanabe, T. Taniguchi, B. A. Bernevig, R. J. Cava, L. M. Schoop, A. Yazdani, and S. Wu, Evidence for a monolayer excitonic insulator, Nat. Phys. 10.1038/s41567-021-01422-w (2021).
  * Bretscher  _et al._ [2021] H. M. Bretscher, P. Andrich, Y. Murakami, D. Golez, B. Remez, P. Telang, A. Singh, L. Harnagea, N. R. Cooper, A. J. Millis, P. Werner, A. K. Sood, and A. Rao, Imaging the coherent propagation of collective modes in the excitonic insulator Ta2NiSe5 at room temperature, Sci. Adv. 7, eabd6147 (2021).
  * Yuan  _et al._ [2018] W. Yuan, Q. Zhu, T. Su, Y. Yao, W. Xing, Y. Chen, Y. Ma, X. Lin, J. Shi, R. Shindou, X. C. Xie, and W. Han, Experimental signatures of spin superfluid ground state in canted antiferromagnet Cr2O3 via nonlocal spin transport, Sci. Adv. 4, eaat1098 (2018).
  * Jiang  _et al._ [2020] Z. Jiang, W. Lou, Y. Liu, Y. Li, H. Song, K. Chang, W. Duan, and S. Zhang, Spin-Triplet Excitonic Insulator: The Case of Semihydrogenated Graphene, Phys. Rev. Lett. 124, 166401 (2020).
  * Nasu and Naka [2021] J. Nasu and M. Naka, Spin seebeck effect in nonmagnetic excitonic insulators, Phys. Rev. B 103, L121104 (2021).
  * Kaneko  _et al._ [2015] T. Kaneko, B. Zenker, H. Fehske, and Y. Ohta, Competition between excitonic charge and spin density waves: Influence of electron-phonon and Hund’s rule couplings, Phys. Rev. B 92, 115106 (2015).
  * Kuneš [2015] J. Kuneš, Excitonic condensation in systems of strongly correlated electrons, J. Phys. Cond. Mat. 27, 333201 (2015).
  * Nasu  _et al._ [2016] J. Nasu, T. Watanabe, M. Naka, and S. Ishihara, Phase diagram and collective excitations in an excitonic insulator from an orbital physics viewpoint, Phys. Rev. B 93, 205136 (2016).
  * Sotnikov and Kuneš [2016] A. Sotnikov and J. Kuneš, Field-induced exciton condensation in LaCoO3, Sci. Rep. 6, 30510 (2016).
  * Tatsuno  _et al._ [2016] T. Tatsuno, E. Mizoguchi, J. Nasu, M. Naka, and S. Ishihara, Magnetic field effects in a correlated electron system with spin-state degree of freedom ― implications for an excitonic insulator ―, J. Phys. Soc. Jpn. 85, 083706 (2016).
  * Nakamura  _et al._ [2018] D. Nakamura, A. Ikeda, H. Sawabe, Y. H. Matsuda, and S. Takeyama, Record indoor magnetic field of 1200 T generated by electromagnetic flux-compression, Rev. Sci. Instrum. 89, 095106 (2018).
  * Ikeda  _et al._ [2017] A. Ikeda, T. Nomura, Y. H. Matsuda, S. Tani, Y. Kobayashi, H. Watanabe, and K. Sato, High-speed 100 MHz strain monitor using fiber Bragg grating and optical filter for magnetostriction measurements under ultrahigh magnetic fields, Rev. Sci. Instrum. 88, 083906 (2017).
  * Ikeda  _et al._ [2020] A. Ikeda, Y. H. Matsuda, and K. Sato, Two spin-state crystallizations in LaCoO3, Phys. Rev. Lett. 125, 177202 (2020).
  * Ikeda  _et al._ [2016] A. Ikeda, T. Nomura, Y. H. Matsuda, A. Matsuo, K. Kindo, and K. Sato, Spin state ordering of strongly correlating $\rm{LaCoO}_{3}$ induced at ultrahigh magnetic fields, Phys. Rev. B 93, 220401(R) (2016).
  * Manago  _et al._ [2019] M. Manago, S. Kitagawa, K. Ishida, K. Deguchi, N. K. Sato, and T. Yamamura, Spin-triplet superconductivity in the paramagnetic UCoGe under pressure studied by 59Co NMR, Phys. Rev. B 100, 035203 (2019).
  * Ran  _et al._ [2019] S. Ran, C. Eckberg, Q. P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I. L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Nearly ferromagnetic spin-triplet superconductivity, Science 365, 684 (2019).
  * Halperin [2019] W. P. Halperin, Superfluid 3He in Aerogel, Annu. Rev. Condens. Matter Phys. 10, 155 (2019).
  * Kaneko and Ohta [2016] T. Kaneko and Y. Ohta, Electric and magnetic multipoles and bond orders in excitonic insulators, Phys. Rev. B 94, 125127 (2016).
  * Hejtmánek  _et al._ [2013] J. Hejtmánek, Z. Jirák, O. Kaman, K. Knížek, E. Šantavá, K. Nitta, T. Naito, and H. Fujishiro, Phase transition in Pr0.5Ca0.5CoO3 and related cobaltites, Eur. Phys. J. B 86, 305 (2013).
  * Kuneš and Augustinský [2014] J. Kuneš and P. Augustinský, Excitonic instability at the spin-state transition in the two-band hubbard model, Phys. Rev. B 89, 115134 (2014).
  * Sotnikov and Kuneš [2017] A. Sotnikov and J. Kuneš, Competing phases in a model of Pr-based cobaltites, Phys. Rev. B 96, 245102 (2017).
  * Rice [2002] T. M. Rice, Quantum mechanics. To condense or not to condense, Science 298, 760 (2002).
  * Zapf  _et al._ [2014] V. Zapf, M. Jaime, and C. U. Batista, Bose-Einstein condensation in quantum magnets, Rev. Mod. Phys. 86, 563 (2014).
  * Sato  _et al._ [2009] K. Sato, A. Matsuo, K. Kindo, Y. Kobayashi, and K. Asai, Field induced spin-state transition in $\rm{LaCoO}_{3}$, J. Phys. Soc. Jpn. 78, 093702 (2009).
  * Altarawneh  _et al._ [2012] M. M. Altarawneh, G. W. Chern, N. Harrison, C. D. Batista, A. Uchida, M. Jaime, D. G. Rickel, S. A. Crooker, C. H. Mielke, J. B. Betts, J. F. Mitchell, and M. J. R. Hoch, Cascade of magnetic field induced spin transitions in $\rm{LaCoO}_{3}$, Phys. Rev. Lett. 109, 037201 (2012).
  * Rotter  _et al._ [2014] M. Rotter, Z. S. Wang, A. T. Boothroyd, D. Prabhakaran, A. Tanaka, and M. Doerr, Mechanism of spin crossover in LaCoO3 resolved by shape magnetostriction in pulsed magnetic fields, Sci. Rep. 4, 7003 (2014).
  * Schönemann  _et al._ [2021] R. Schönemann, G. Rodriguez, D. Rickel, F. Balakirev, R. D. McDonald, J. A. Evans, B. Maiorov, C. Paillard, L. Bellaiche, A. V. Stier, M. B. Salamon, K. Gofryk, and M. Jaime, Magnetoelastic standing waves induced in UO2 by microsecond magnetic field pulses, Proc. Natl. Acad. Sci. U. S. A. 118, e2110555118 (2021).
  * sup [a] (a), see section I in the Supplemental Material.
  * Kyômen  _et al._ [2003] T. Kyômen, Y. Asaka, and M. Itoh, Negative cooperative effect on the spin-state excitation in $\rm{LaCoO}_{3}$, Phys. Rev. B 67, 144424 (2003).
  * [35] T. Kihara, Y. Kohama, M. Tokunaga, and K. Sato, Unpublished .
  * Noguchi  _et al._ [2002] S. Noguchi, S. Kawamata, K. Okuda, H. Nojiri, and M. Motokawa, Evidence for the excited triplet of LaCoO3, Phys. Rev. B 66, 094404 (2002).
  * sup [b] (b), see section II in the Supplemental Material.
  * [38] K. Tomiyasu, T. Nomura, Y. Kobayashi, S. Ishihara, S. Ohira, and M. Kofu, Emergent units of itinerant spin-state excitations in LaCoO3, arXiv:1808.05888 .
  * Hariki  _et al._ [2020] A. Hariki, R.-P. Wang, A. Sotnikov, K. Tomiyasu, D. Betto, N. B. Brookes, Y. Uemura, M. Ghiasi, F. M. F. d. Groot, and J. Kuněs, Melting of excitonic dispersion in LaCoO3: theory and experiment, Phys. Rev. B 101, 245162 (2020).
  * Radaelli and Cheong [2002] P. G. Radaelli and S. W. Cheong, Structural phenomena associated with the spin-state transition in $\rm{LaCoO}_{3}$, Phys. Rev. B 66, 094408 (2002).
  * Okimoto  _et al._ [2021] Y. Okimoto, T. Saitoh, Y. Kobayashi, and S. Ishihara, eds.,  _Spin-Crossover Cobaltite (Review and Outlook)_, Springer Series in Materials Science (Springer Singapore, 2021).
  * Raccah and Goodenough [1967] P. M. Raccah and J. B. Goodenough, First-Order Localized-Electron $\rightleftharpoons$ Collective-Electron Transition in LaCoO3, Phys. Rev. 155, 932 (1967).
  * Knížek  _et al._ [2009] K. Knížek, Z. Jirák, J. Hejtmánek, P. Novák, and W. Ku, Gga+u calculations of correlated spin excitations in $\rm{LaCoO}_{3}$, Phys. Rev. B 79, 014430 (2009).

## Methods

The magnetic fields up to 600 T are generated using electromagnetic flux
compression method [16] in the Institute for Solid State Physics, the
University of Tokyo, Japan. The magnetostriction measurement is performed
using our state-of-the-art high-speed strain gauge utilizing fiber Bragg
grating and optical filter method [17]. Polycrystalline samples of LaCoO3 were
prepared by the solid-state-reaction of predried La2O3 and CoO at 1300 C∘ in
air. The sample is cut out in $1\times 3\times 0.5$ mm. FBG fiber is glued
tightly onto a specimen of the sample using Stycast 1266. See Supplemental
Material for details [33].

## Awknowledgements

The work is financially supported by JSPS KAKENHI Grant-in-Aid for Early-
Career Scientists Grant No. 18K13493 and Basic Science Program No. 18-001 of
Tokyo Electric Power Company (TEPCO) memorial foundation.

## Author contributions

A.I. designed the research. S.K prepared the sample. A.I., Y.H.M., H.S., D.N.,
S.T. performed the experiments. A.I. analysed the data. J.N. and A.I.
developed the model. J.N. performed the calculations. A.I Y.H.M., J.N.
discussed the results. A.I. wrote the manuscript with input from all co-
authors.

## Competing interests

The authors declare no competing interests.