# Seeding of the Self-Modulation in a Long Proton Bunch by Charge Cancellation
with a Short Electron Bunch

Mathias Hüther  Patric Muggli  Max Planck Institute for Physics, Munich,
Germany  Technical University of Munich, Garching b. München, Germany  CERN,
Geneva, Switzerland

###### Abstract

In plasma wakefield accelerators (e.g. AWAKE) the proton bunch self-modulation
is seeded by the ionization front of a high-power laser pulse ionizing a
vapour and by the resulting steep edge of the driving bunch profile inside the
created plasma.
In this paper, we present calculations in 2D linear theory for a concept of a
different self-modulation seeding mechanism based on electron injection. The
whole proton bunch propagates through a preformed plasma and the effective
beam current is modulated by the external injection of a short electron bunch
at the centre of the proton beam. The resulting sharp edge in the effective
beam current in the trailing part of the proton bunch is driving large
wakefields that can lead to a growth of the seeded self-modulation (SSM).
Furthermore, we discuss the feasibility for applications in AWAKE Run 2.

###### keywords:

AWAKE , Seeding of Proton Bunch Self-Modulation , Plasma Wakefield
Acceleration , Electron Injection

††journal: Nuclear Instruments and Methods in Physics Research A

##  1 Introduction

The proton bunch self-modulation in plasma wakefield accelerators (e.g. the
Advanced Wakefield Experiment (AWAKE) located at CERN [1, 2]) is usually
seeded by using a high-power laser pulse co-propagating in the centre of the
proton bunch (see Fig. 1(a)). The laser ionizes a vapour and creates a plasma,
resulting in a sharp relativistic ionization front separating vapour and
plasma. Due to this ionization front, the proton bunch self-modulation is
seeded and consequently growing within a distance of a few meters before
reaching its saturation level [3, 4]. The long proton bunch is split into a
train of micro-bunches with a period on the order of the plasma wavelength
[5]. This modulation is caused by the self-modulation instability (SMI), a
transverse beam-plasma instability. A witness electron bunch can thus
deterministically be injected into one of the buckets between the micro-
bunches and accelerated.
The seeding concept presented in this paper gives another approach to reach
this goal: A high-power laser pulse is ionizing the plasma some time ahead of
the proton bunch, so that the proton bunch is propagating through a preformed
plasma (see Fig. 1(b)). Hence, there is no seeding from the ionization
process. A short electron bunch is injected at the centre or somewhere along
the long proton bunch. The resulting local charge cancellation leads to a
distinct drop in the effective charge density distribution - and is therefore
driving wakefields in the rear part of the proton bunch. The front of the
bunch is still long with a smoothly increasing density and thus drives much
lower amplitude wakefields. Thanks to the seeding process [6], the wakefield
phase and amplitude is determined and weakly dependent on initial beam
parameters [7]. The witness electrons to be accelerated are again injected a
number of buckets behind the seeding, as in the usual injection scheme. An
overview of the different self-modulation seeding approaches using a short
electron bunch is given in Fig. 1. A variation of the approach described in
this paper, where a seed electron bunch is put ahead of the proton bunch, is
shown in Fig. 1(c) and will be explained later.
AWAKE is the world’s first proton-driven plasma wakefield accelerator aiming
for acceleration of externally injected electrons in gradients up to the GV/m
scale. It uses a 8-12 cm long proton bunch from CERN’s Super Proton Synchotron
(SPS) that propagates through a 10 m long Rubidium plasma (with a density of
$(1-10)\cdot 10^{14}\,$cm-3), induced by a high-power, short laser pulse
(duration $\tau=120\,$fs, energy $E=450\,$mJ).

(a)

(b)

(c)

Figure 1: Self-modulation seeding and electron injection schemes. (a) Current
seeding concept with wakefields seeded by a vapour-plasma transition front
created by a high-power laser pulse in the centre of the proton beam. (b)
Novel seeding scheme with an electron bunch injected in the centre of the long
proton bunch and resulting in a dip in the effective charge density. The laser
pulse is ionizing the vapour ahead of the proton bunch. (c) Laser pulse and
electron bunch ahead of the proton bunch. Seeding of SM process due to
wakefields driven by the electron bunch.

##  2 A Simplified Model in Linear Theory

For a more quantitative understanding of the approach towards seeding, we
introduce the following simplified model: The proton and electron bunch
longitudinal density distribution is described by a $cos$-profile, whereas
their radial profiles are Gaussian. We use 2D-linear plasma wakefield theory
[8] to describe the seeding effect by calculating the initial wakefields.

Figure 2: $\cos$-approximation of the longitudinal proton bunch density
$n_{b\parallel}$ as described by Eq. 1.

The total (proton and electron) longitudinal charge distribution
$n_{b\parallel}$ (see Fig. 2) is given by

$n_{b\parallel}(\xi^{\prime})=\begin{cases}A_{0}\cdot\sin\left(\frac{\pi\cdot(\xi^{\prime}+L)}{2(L-U)}\right)&\xi\in[-L,-U]\\\
\frac{A_{min}-A_{0}}{2}\cdot\cos\left(\frac{\pi\xi^{\prime}}{U}\right)+\frac{A_{min}+A_{0}}{2}&\xi\in[-U,U]\\\
A_{0}\cdot\cos\left(\frac{\pi\cdot(\xi^{\prime}-U)}{2(L-U)}\right)&\xi\in[U,L]\\\
0&~{}\text{elsewhere},\end{cases}$ (1)

where $\xi=ct-z$ is the spatial coordinate in the co-moving bunch frame, $2L$
the longitudinal size of the proton bunch and $2U$ the longitudinal size of
the electron bunch. Outside $[-L,L]$ the bunch charge density is assumed to be
zero. The amplitude of the proton charge density is $A_{0}$ and the amplitude
of the electron charge density is $A_{0}-A_{min}$ . The absolute number of
charges in one bunch for our model is given by
$N=8A_{0}Ln_{b}\sigma_{r,p}^{2}\cdot 0.683$. For comparison in case of a
Gaussian longitudinal bunch profile $N$ is given by
$N=(2\pi)^{\frac{3}{2}}n_{b}\sigma_{z,p}\sigma_{r,p}^{2}$. According to [8]
the parallel and perpendicular wakefields $W_{z}(\xi,r)$ and $W_{r}(\xi,r)$
are given by

$W_{z}(\xi,r)=\frac{e}{\epsilon_{0}}~{}~{}\int_{-\infty}^{\xi}n_{b\parallel}(\xi^{\prime})\cos\left(k_{pe}(\xi-\xi^{\prime})\right)d\xi^{\prime}\cdot~{}R(r)$
(2)

and

$W_{r}(\xi,r)=\frac{e}{k_{pe}\epsilon_{0}}\displaystyle\int_{-\infty}^{\xi}n_{b\parallel}(\xi^{\prime})\sin\left(k_{pe}(\xi-\xi^{\prime})\right)d\xi^{\prime}\cdot\frac{dR(r)}{dr},$
(3)

where $e$ is the elementary charge, $\epsilon_{0}$ the vacuum permittivity and
$k_{pe}=\omega_{pe}/v_{b}\approx\omega_{pe}/c$ the plasma wave number in case
of a relativistic particle bunch with bunch velocity $v_{b}\approx c$. The
plasma electron angular frequency $\omega_{pe}$ is

$\omega_{pe}=\sqrt{\frac{n_{pe}e^{2}}{\epsilon_{0}m_{e}}},$ (4)

where $m_{e}$ is the electron mass and $n_{pe}$ the plasma charge density. The
transverse component $R(r)$ in Eqn. 2 and 3 is given by

$\displaystyle R(r)=k_{pe}^{2}~{}K_{0}(k_{pe}r)$
$\displaystyle\int_{0}^{r}r^{\prime}n_{b\perp}(r^{\prime})~{}I_{0}(k_{pe}r^{\prime})~{}dr^{\prime}$
$\displaystyle+k_{pe}^{2}~{}I_{0}(k_{pe}r)$
$\displaystyle\int_{r}^{\infty}r^{\prime}n_{b\perp}(r^{\prime})~{}K_{0}(k_{pe}r^{\prime})~{}dr^{\prime},$
(5)

with $I_{0}$ and $K_{0}$ the modified Bessel functions of the first and second
kind, respectively. Assuming the radial profile and the radial size of the
electron and proton bunch are the same, $R(r)$ and $\frac{dR}{dr}$ are just
multiplying factors. Considering different radii is trivial.
Piecewise integration of the separate bunch parts of Eqs. 2 and 3 along $\xi$
(see Fig. 2) gives the contribution of each single part to the resulting
wakefields.
The wakefields driven by the leading part of the proton bunch $(-L<\xi<-U)$
are described by

$\footnotesize E_{1,\parallel~{}\,}(\xi)=\begin{cases}0\\\
\frac{A_{0}}{2k_{pe}}\left[-(C_{+}+C_{-})\left[\cos\left(\frac{\pi(\xi+L)}{2(L-U)}\right)+\cos(k_{pe}(\xi+L))\right]\right]\\\
\frac{A_{0}}{2k_{pe}}\left[-(C_{+}-C_{-})\sin(k_{pe}(\xi+U))+(C_{+}+C_{-})\cos(k_{pe}(\xi+L))\right]\\\
\end{cases}$ (6)

and

$\footnotesize E_{1,\perp~{}}(\xi)=\begin{cases}0\\\
\frac{A_{0}}{2}\left[(C_{+}-C_{-})\sin\left(\frac{\pi(\xi+L)}{2(L-U)}\right)+(C_{+}+C_{-})\sin\left(k_{pe}(\xi+L)\right)\right]\\\
\frac{A_{0}}{2}\left[(C_{+}-C_{-})\cos(k_{pe}(\xi+U))+(C_{+}+C_{-})\sin(k_{pe}(\xi+L))\right],\\\
\end{cases}$ (7)

where the first term is valid for $\xi<-L$, the second term for
-$L\leq\xi\leq-U$ and the third term for $\xi>-U$. The parallel and
perpendicular wakefields in the part where the electrons are injected and
therefore cause a drop in the effective charge density are given by

$\footnotesize E_{2,\parallel}(\xi)=\begin{cases}0\\\
\frac{A_{1}}{2k_{pe}}\left[(D_{+}+D_{-})\sin\left(\frac{\pi\xi}{U}\right)-(D_{+}-D_{-})\sin(k_{pe}(\xi+U))\right]\\\
~{}~{}~{}+\frac{A_{2}}{k_{pe}^{2}}\sin(k_{pe}(\xi+U))\\\
\frac{A_{1}}{2k_{pe}}\left[(D_{+}-D_{-})\sin(k_{pe}(\xi-U))-(D_{+}-D_{-})\sin(k_{pe}(\xi+U))\right]\\\
~{}~{}~{}-\frac{A_{2}}{k_{pe}^{2}}\left[\sin(k_{pe}(\xi-U))-\sin(k_{pe}(\xi+U))\right]\\\
\end{cases}$ (8)

and

$\footnotesize E_{2,\perp}(\xi)=\begin{cases}0\\\
\frac{A_{1}}{2}\left[(D_{+}+D_{-})\cos\left(\frac{\pi\xi}{U}\right)-(D_{+}-D_{-})\cos(k_{pe}(\xi+U))\right]\\\
~{}~{}~{}+\frac{A_{2}}{k_{pe}}\left[1-\cos(k_{pe}(\xi+U))\right]\\\
\frac{A_{1}}{2}\left[(-D_{+}+D_{-})\cos(k_{pe}(\xi-U))+(D_{+}-D_{-})\cos(k_{pe}(\xi+U))\right]\\\
~{}~{}~{}-\frac{A_{2}}{k_{pe}}\left[\cos(k_{pe}(\xi+U))-\cos(k_{pe}(\xi-U))\right].\\\
\end{cases}$ (9)

Here, the first term is valid for $\xi<-U$, the second for $-U\leq\xi\leq U$
and the third for $\xi>U$.
The wakefields caused by bunch section between $U$ and $L$ are described by

$\footnotesize E_{3,\parallel~{}}(\xi)=\begin{cases}0\\\
\frac{A_{0}}{2k_{pe}}\left[(C_{+}+C_{-})\sin\left(\frac{\pi(\xi-U)}{2(L-U)}\right)+(C_{+}-C_{-})\sin(k_{pe}(\xi-U))\right]\\\
\frac{A_{0}}{2k_{pe}}\left[(C_{+}+C_{-})\cos(k_{pe}(\xi-L))+(C_{+}-C_{-})\sin(k_{pe}(\xi-U))\right]\\\
\end{cases}$ (10)

and

$\footnotesize E_{3,\perp~{}}(\xi)=\begin{cases}0\\\
\frac{A_{0}}{2}\left[(C_{+}-C_{-})\cos\left(\frac{\pi(\xi-U)}{2(L-U)}\right)+(-C_{+}+C_{-})\cos(k_{pe}(\xi-U))\right]\\\
\frac{A_{0}}{2}\left[(C_{+}+C_{-})\sin(k_{pe}(\xi-L))+(-C_{+}+C_{-})\cos(k_{pe}(\xi-U))\right].\\\
\end{cases}$ (11)

For this case, the first line is valid for $\xi<U$, the second for the
interval $U\leq\xi\leq L$ and the third line for $\xi>L$.
The constants in Eqn. 6 to 11 are given by $A_{1}=\frac{A_{min}-A_{0}}{2}$,
$A_{2}=\frac{A_{min}+A_{0}}{2}$, $C\pm=\frac{1}{\frac{\pi}{2(L-U)}\pm k_{pe}}$
and $D_{\pm}=\frac{1}{\frac{\pi}{U}\pm k_{pe}}$.
Making use of the superposition theorem of waves in linear theory, the final
wakefield distribution can be determined by adding the different contributions
of the bunch segments in every region. Hence, the resulting wakefields are
described by

$E(\xi)=\begin{cases}0&~{}\xi<-L\\\
E_{1}(\xi)\big{|}_{\xi\in\left[-L;-U\right]}&\xi\in\left[-L;-U\right]\\\
E_{1}(\xi)\big{|}_{\xi\in\left[-U;U\right]}+E_{2}(\xi)\big{|}_{\xi\in\left[-U;U\right]}&\xi\in\left[-U;U\right]\\\
E_{1}(\xi)\big{|}_{\xi\in\left[U;L\right]}+E_{2}(\xi)\big{|}_{\xi\in\left[U;L\right]}+E_{3}(\xi)\big{|}_{\xi\in\left[U;L\right]}&\xi\in\left[U;L\right]\\\
E_{1}(\xi)\big{|}_{\xi>L}+E_{2}(\xi)\big{|}_{\xi>L}+E_{3}(\xi)\big{|}_{\xi>L}&\xi>L.\\\
\end{cases}$ (12)

Fig. 3 shows the longitudinal and perpendicular wakefields (Eq. 12) for a
ratio $\frac{U}{L}=\frac{1}{30}$ and $\frac{U}{\lambda_{pe}}=1.6$. The front
of the proton bunch ($\xi<-U$) propagates in a preformed plasma where its
length is very long when compared to the plasma wavelength $\lambda_{pe}$ and
its profile is smooth. It is therefore not effective at driving wakefields,
but still drives a low amplitude fields that could lead to development of the
SMI process over long plasma distances. The adiabatic response of the plasma
generates a globally focussing force by charge neutralisation that, upon
propagation, can lead to an increase in bunch density and thus also of the
wakefield response. This effect could also lead to SMI growth (as opposed to
SSM) over long plasma distances. From the position of the electron bunch
($\xi>U$) and all along the second half of the proton bunch, much larger
wakefields are driven, which subsequently provide the strong seed for the SM
process.
For bunch parameters similar to AWAKE ($L=$ 6 cm, $U=$ 2 mm,
$\frac{U}{L}=\frac{1}{30}$, $\frac{U}{\lambda_{pe}}=1.6$, peak currents $\mid
Q_{p^{+},peak}\mid=\mid Q_{e^{-},peak}\mid$ (see Fig. 3) and
$k_{pe}=5000\,\text{m}^{-1}$, $n_{p}=7\cdot 10^{14}\,\text{cm}^{-3}$,
$\sigma_{r}=200\,\upmu\text{m}$) the wakefields in the trailing part of the
bunch could reach values up to the GV/m scale.

Figure 3: Parallel and perpendicular wakefields $E_{\parallel}(\xi)$ and
$E_{\perp}(\xi)$ for a bunch shape with $\sigma_{z,p}=$ 12 cm, $\sigma_{z,e}=$
4 mm, $\frac{U}{L}=\frac{1}{30}$, $\frac{U}{\lambda_{pe}}=1.6$, peak currents
$\mid Q_{p^{+},peak}\mid=\mid Q_{e^{-},peak}\mid$.

It can be easily seen from Eqn. 6 \- 11 that the amplitudes of the wakefields
in the trailing part of the bunch are scaling linearly with the electron bunch
charge $Q_{e}$, i.e. the depth of the gap in the effective charge density
$A_{0}-A_{min}$.
Furthermore, the wakefields are strongly depending on the width of the gap in
the effective charge density $2U$, corresponding to the longitudinal size of
the injected electron bunch $\sigma_{z,e}$. The amplitude of the wakefields
can be maximized for a given configuration of $A_{min}$ and $L$ and follow the
distribution shown in Fig. 4. For an electron bunch length $U$ on the order of
$\lambda_{pe}$, the peak wakefields are the highest (see Fig. 4).

Figure 4: Maximum wakefields for different ratios of $2U/\lambda_{pe}$.

##  3 Discussion of the Method

The seeding method described above could be an interesting approach for a
seeding concept for the planned extension of the AWAKE experiment after 2021,
called AWAKE Run 2 [9]. Current design studies for the plasma accelerator
consist of a split plasma, a short one ($\sim$ 4 m) for the seeding of the SSM
and a second, longer one ($\sim$ 10 m) driven by the modulated proton bunch
for the acceleration of the electrons (see Fig. 5). The injection point for
the witness electron bunch is foreseen to be between the two plasmas. The
short plasma source is similar to the one used in AWAKE Run 1, whereas the
long plasma cell could be either a Helicon or a discharge source [10].
By seeding the SMI with an electron bunch, as described here, the ionisation
front is not necessary for seeding and thus the first plasma could also be
preformed in one of these two types of sources. Hence, there would be no
longer need for a maintenance-intensive high-power laser system nor for the
complicated Rubidium handling and storage procedures required by a Rubidium
based vapour source.

(a)

(b)

Figure 5: Different schemes for a possible design of the split plasma cells in
AWAKE Run 2. (a) Design with laser seeding and (b) design with seeding by
electron injection.

The major uncertainty concerning the feasibility of this approach is the
effect of the first part of the proton bunch. Even though, according to the
simulations of Kumar et. al. [5] and according to the wakefields described in
Eq. 12, there is no strong seeding of the SMI expected by the leading part of
the bunch, the weak amplitude wakefields or those from noise may grow.
Subsequently, wakefields in the front of the bunch would interfere with the
seeded wakefields in the trailing part of the bunch. Hence, a phase-stable
seeding of the SSM would no longer be possible.
First preliminary experimental results from the AWAKE-experiment show some
evidence for a growth of wakefields along the bunch on a scale of a few meters
as well as micro-bunching even without seeding, i.e. propagation in a
preformed plasma.
In case the unseeded SM-growth in the leading part of the proton bunch is too
high for a seeding with an electron bunch in the centre of the beam, another
option would be to shift the electron bunch to an earlier position with
respect to the proton bunch or even ahead of it (see Fig. 1(c)). With this,
the full proton bunch charge could be used for driving wakefields, while
removing the risk of growth of SMI ahead of the seeding point.

##  4 Conclusions

We presented an alternative concept for the seeding of the SSM in plasma-
wakefield accelerators. Instead of seeding by an ionisation front created by a
high-power laser pulse, the seeding is archived by the injection of a short
electron bunch in the centre of a long proton bunch propagating through a
preformed plasma. Calculations in 2D linear theory show that the resulting
steep rising edge in the effective charge density of the proton bunch drives
large amplitude wakefields that seed the SSM. Although the approach is very
promising and could simplify the design for AWAKE Run 2, first experimental
results from the latest AWAKE Runs give evidence that wakefields do grow from
noise in the front of the proton bunch, when the bunch propagates in a
preformed plasma and no seeding is provided.

##  Appendix A References

## References

  * [1] P. Muggli, et al., Physics of the AWAKE project, in: IPAC 2013: Proceedings of the 4th International Particle Accelerator Conference, 2013, pp. 1179–1181.
  * [2] E. Gschwendtner, et al., AWAKE, The Advanced Proton Driven Plasma Wakefield Acceleration Experiment at CERN, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 829 (2016) 76–82.
  * [3] K. V. Lotov, A. P. Sosedkin, A. V. Petrenko, Long-term evolution of broken wakefields in finite-radius plasmas, Phys. Rev. Lett. 112 (2014) 194801.  doi:10.1103/PhysRevLett.112.194801.
  * [4] A. Caldwell, et al., Path to AWAKE: Evolution of the concept, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 829 (2016) 3–16.
  * [5] N. Kumar, A. Pukhov, K. Lotov, Self-modulation instability of a long proton bunch in plasmas, Physical Review Letters 104 (25) (2010) 255003.
  * [6] K. Lotov, et al., Natural noise and external wakefield seeding in a proton-driven plasma accelerator, Physical Review Special Topics-Accelerators and Beams 16 (4) (2013) 041301.
  * [7] N. Savard, P. Muggli, J. Vieira, Effect of proton bunch parameter variation on awake, in: North American Particle Accelerator Conf.(NAPAC’16), Chicago, IL, USA, October 9-14, 2016, JACOW, Geneva, Switzerland, 2017, pp. 684–686.
  * [8] P. Chen, et al., Plasma focusing for high-energy beams, IEEE transactions on plasma science 15 (2) (1987) 218–225.
  * [9] E. Adli, Towards AWAKE applications: Electron beam acceleration in a proton driven plasma wake, in: 7th Int. Particle Accelerator Conf.(IPAC’16), Busan, Korea, May 8-13, 2016, JACOW, Geneva, Switzerland, 2016, pp. 2557–2560.
  * [10] AWAKE Collaboration, AWAKE Status Report 2016, SPSC-SR-194.


# Emergence of nonlinear dynamics from spatial structure
in tropical forest-grassland landscapes

Bert Wuyts  College of Engineering, Mathematics and Physical Sciences,
University of Exeter, EX4 4QF, UK  Jan Sieber  College of Engineering,
Mathematics and Physical Sciences, University of Exeter, EX4 4QF, UK
<EMAIL_ADDRESS>
###### Abstract

It is thought that tropical forests can exist as an alternative stable state
to savanna Staver  _et al._ (2011); Hirota  _et al._ (2011). Therefore, the
cumulative effects of perturbation by climate change or human impact may lead
to crossing of a tipping point beyond which there is rapid large-scale forest
dieback that is not easily reversed Lenton  _et al._ (2008); Scheffer  _et
al._ (2015). Empirical evidence for bistability due to fire-vegetation
feedbacks relies on tree cover bimodality in satellite-observed data Staver
_et al._ (2011); Hirota  _et al._ (2011), but this may also be explained by
spatial heterogeneity Wuyts  _et al._ (2017), or by biases in the data Hanan
_et al._ (2014); Gerard  _et al._ (2017). Most modelling studies of
alternative stable tree cover states have so far either relied on mean-field
assumptions Staver and Levin (2012); Schertzer  _et al._ (2015); Touboul  _et
al._ (2018); Goel  _et al._ (2020); Wuyts  _et al._ (2017, 2019) or not
included the spatiotemporal dynamics of fire Patterson  _et al._ (2021),
making it hard to compare model results to spatial data. In this work, we
analyse a microscopic model of tropical forest dynamics and fire spread,
proposed by Hébert-Dufresne  _et al._ (2018), to show how dynamics of forest
area emerge from spatial structure. We find that the relation between forest
perimeter and area determines the nonlinearity in forest growth while forest
perimeter weighted by adjacent grassland area determines the nonlinearity in
forest loss. Together with the linear changes, which are independent of
spatial structure, these two effects lead to an emergent relation between
forest area change and forest area, defining a single-variable ordinary
differential equation. Such a relation between pattern and dynamics offers
falsifiable theoretical predictions of the bistability hypothesis that are
more closely linked to the underlying mechanism than bimodality and provides a
criterion for forest vulnerability.

tropical forest $|$ bistability $|$ cellular automata $|$ percolation $|$
nonlinear dynamics $|$ mean field $|$ feedback control

##  I Introduction

Satellite observations of tropical tree cover show that in a wide range of
mean annual rainfall Staver  _et al._ (2011); Hirota  _et al._ (2011), both
savanna (low tree cover) and forest (high tree cover) are possible, making the
distribution of tree cover _bimodal_. On the one hand, fire exclusion
experiments show that fire can maintain low tree cover Bond (2008). On the
other hand, fire occurs almost exclusively below a tree cover threshold of
about 40% Staver  _et al._ (2011); Hoffmann  _et al._ (2012); Archibald  _et
al._ (2009); Wuyts  _et al._ (2017); Van Nes  _et al._ (2018), which is
consistent with fire being a _percolation_ process Stauffer and Aharony
(1994); Christensen and Moloney (2005) on grass patches Schertzer  _et al._
(2015); Cardoso  _et al._ (2022). Such a highly nonlinear response of fire to
grass together with an empirically consistent response of vegetation to fire
was shown to be sufficient for inducing _bistability_ in simple models Staver
and Levin (2012). Taken together, the bimodality, the two-way interaction
between fire and vegetation, and the availability of a plausible underlying
mechanism provide compelling evidence that forest and savanna are _alternative
stable states_ , maintained by a feedback between vegetation and fire Staver
_et al._ (2011), and between which transitions would neither be gradual nor
easily reversed Lenton  _et al._ (2008); Scheffer  _et al._ (2015). This
local-scale fire-vegetation feedback may exacerbate the effects of the
regional-scale feedback between vegetation and rainfall Spracklen  _et al._
(2018).

Nevertheless, subsequent analysis of the same data showed that at least part
of the bimodality can be explained by independent spatial heterogeneity of
parameters other than mean annual rainfall (Wuyts  _et al._ , 2017). Yet
another possibility is that the bimodality results from artificial
discontinuities in the data, introduced by its generative machine-learning
algorithm Hanan  _et al._ (2014); Gerard  _et al._ (2017). Therefore, the
fire-mediated alternative stable state hypothesis is still contested, with
some left unconvinced and in favour of other causal mechanisms that
distinguish between savanna and forest Veenendaal  _et al._ (2018). Taking
into consideration the possible alternative explanations for bimodality,
further tests of bistablity using remote-sensed data would have to rely on
more specific indicators than bimodality. For their specificity, they should
ideally be closely tied to the assumed underlying mechanism, while for
comparability with spatial data, they should relate to space.

Indicators are provided by theoretical models. Despite the fact that tropical
vegetation dynamics is a high-dimensional spatial stochastic process, its
conceptual analysis has mostly been concerned with models consisting of a
small number of ordinary differential equations (ODE) representing macroscale
dynamics (e.g. Van Langevelde  _et al._ , 2003; Accatino  _et al._ , 2010;
Staver and Levin, 2012), possibly including a noise D’Odorico  _et al._
(2006); Touboul  _et al._ (2018) or impulsive (e.g. Baudena  _et al._ , 2010;
Beckage  _et al._ , 2011; Hoyer-Leitzel and Iams, 2021) term. Such an approach
is justified under validity of _mean-field assumptions_ : (_i_) well
mixedness, allowing the probability of interaction between two species to be
written as the product of their overall occurrence probability, (_ii_) large
system size and validity of the law of large numbers in the deterministic case
or an appropriate limit theorem for the noise term in the stochastic case,
(_iii_) spatial homogeneity. E.g., the Staver-Levin model of tropical tree
cover bistability Staver and Levin (2012) is an ODE in which the variables
represent the system-scale fractions of area covered by grass, forest trees,
savanna trees or savanna saplings. Species interact through competition for
space or through the mediating effect of grass cover on fire occurrence. By
the well-mixedness assumption, competition terms involve products of these
macroscopic cover fractions and fire effects on tree cover are functions of
macroscale grass cover multiplied by macroscale tree cover. Taking account of
the threshold dependence of fire on grass cover, with fire affecting forest
trees, but not the fire-adapted savanna trees, except for their recruitment
from saplings, the Staver-Levin model exhibits tree cover bistability but also
oscillatory behaviour Staver and Levin (2012); Touboul  _et al._ (2018).
Assuming well-mixedness allows a greatly simplified description, but it comes
at the expense of disregard for spatial structure. Therefore, mean-field
models cannot provide predictions or indicators that relate to spatial
structure. Spatial mean-field models of tree cover bistability, such as
_reaction-diffusion_ Wuyts  _et al._ (2017, 2019); Goel  _et al._ (2020) and
_integro-partial-differential_ Patterson  _et al._ (2021) equations, partially
alleviate this problem by allowing for spatial variation on the mesoscale.
They show that, due to stationary fronts in heterogeneous environments,
bistability does not necessarily lead to bimodality, unless interfacial
effects are large Goel  _et al._ (2020) or the forest and savanna components
are sufficiently decoupled Wuyts  _et al._ (2019). However, they still rely on
mean-field assumptions between the meso- and microscale, thereby still
neglecting the effects of microscopic spatial structure. As trees or forest
patches are typically highly clustered in space, this neglect of structure by
the well-mixedness assumption is usually unjustified.

Contrary to mean-field models, _microscopic spatial simulation models_ , such
as interacting particle systems or cellular automata, capture detailed spatial
structure, thereby also allowing them to provide indicators that are testable
by spatial data. The _interacting particle system_ of Patterson et al.
Patterson  _et al._ (2021) reproduces many of the behaviours of the Staver-
Levin model, to which it was proven to converge under mean-field assumptions,
but it also shows front solutions, as in the spatial mean-field models, and
finite-size effects, such as extinction after a finite simulation time.
However, it did not include fire as a percolation process, even though
percolation was assumed to play the key role in inducing nonlinearity and
bistability Staver and Levin (2012); Schertzer  _et al._ (2015); Archibald
_et al._ (2012); Cardoso  _et al._ (2022). The cellular automaton of Hébert-
Dufresne et al. Hébert-Dufresne  _et al._ (2018), which can be considered as a
simplified version of Favier  _et al._ (2004), captures the joint dynamics of
forest cover and fire percolation. In their model, fires spread rapidly and
nearly perfectly on grass but can only burn forests at their edge. Besides
bistability of forest cover, they also found that isolated forest patches
obtain a steady state shape characterised by a perimeter-area power law
relation consistent with fits based on remote-sensed observations in 86
protected areas of the Brazilian Cerrado. Stochastic spatial simulation models
are, however, far harder to analyse than mean-field models. In particular,
simulations only give access to the stable parts of phase space, unless one is
prepared to wait prohibitively long.

In this work, we analyse the nonlinear dynamics of the cellular automaton
proposed by Hébert-Dufresne  _et al._ (2018) without relying on mean-field
assumptions. Instead, we start from the observation that fire spreads on grass
much faster than forest and assume that it spreads perfectly, reaching the
entire perimeter of the ignited grass patch, and damaging the forest at this
perimeter. This assumption makes two structural quantities relevant for the
long-term evolution of the forest-grass landscape: the _forest perimeter_ for
gain from forest spread, and, the _grassland-weighted forest perimeter_ ,
which counts each part of the forest perimeter weighted by the size of its
adjacent grass patch, for loss induced by fire. Via simulations, we test how
these perimeter quantities affect the dynamics. We are able to explore the
unstable parts of phase space by first obtaining unstable solutions from
feedback-controlled simulations. This will show that both perimeter quantities
are approximate functions of forest area throughout the dynamics, allowing the
evolution of forest area to be written as a single-variable ODE. In this
equation, the nonlinear terms emerge from competition of forest and grass
patches at their perimeter, while the linear terms arise from spatially
independent spontaneous conversions.

##  II Results

###  II.1 The FGBA probabilistic cellular automaton

Figure 1: The FGBA stochastic cellular automaton: (a) state transition diagram
(coloured rates: spread to neighbouring cell, black rates: spontaneous
conversion within cell), (b) example time series of a simulation starting at
zero tree cover, (c,d) $10$x and $10^{7}$x zoom of (b), (e-h) snapshots of a
simulation at indicated times for low fire ignition rate ($\phi N=0.075$). The
fire in (f) spreads throughout grassland in the whole domain whereas that in
(g) went extinct locally because forest splits grassland in clusters (notice
the area of ash near the top). Remaining parameters are shown in Table M1.
Domain size: 200x200 cells. Figure 2: Steady states and bistability of forest
area in the FGBA cellular automaton. (a) Bifurcation diagram of forest area
fraction $[F]$ versus fire ignition rate $\phi$ (shade: two-standard deviation
confidence interval of the mean). (b) Simulations initiated at two different
points on the saddle ($\phi N{=}0.257$ and $\phi N{=}1.32$). Remaining
parameters are shown in Table M1. Domain size: 100x100 cells.

The FGBA probabilistic cellular automaton (adapted from (Hébert-Dufresne  _et
al._ , 2018) – see Methods) models the stochastic dynamics of tropical
vegetation and fire on a square lattice and in continuous time. The key
empirical facts of tropical forest and fire dynamics captured by the FGBA
automaton are: (_i_) fires only naturally ignite in grasslands but they can
spread into forest, (_ii_) fires spread faster in grassland than in forest,
such that forests suppress fires, albeit imperfectly, (_iii_) forest dynamics
occurs on a _strongly separated time scale_ from fire spread and grass
regrowth.

This results in the following reaction rules in the cellular automaton. At any
time, each lattice cell can be in one of four states: $F$ – forest, $G$ –
grass, $B$ – burning, $A$ – ash. Conversions between these states can occur
spontaneously or due to spread to neighbouring cells (Fig. 1a and Table M1).
The spontaneous conversions are: forest recruitment on grass cells due to
long-distance seed dispersal ($G\to F$ at rate $\beta$), forest mortality
($F\to G$ at rate $\gamma$), fire ignition on grass cells ($G\to B$ at rate
$\phi$), and grass regrowth on ash cells ($A\to G$ at rate $\lambda$). The
conversions due to spread to neighbours are: forest recruitment due to short-
distance seed dispersal on grass or ash cells ($GF\to FF,\;AF\to FF$ at rate
$\alpha$), fire spread on grass ($GB\to BB$ at rate $\rho_{g}$) or on tree
cells ($FB\to BB$ at rate $\rho_{f}$). Chosen parameters are in the ranges
empirically justified by Hébert-Dufresne  _et al._ (2018) for a square domain
of $100$x$100$ cells, with cell size $\Delta x=\Delta y=30$m. The time scale
separation between fire and forest dynamics implies that the rates satisfy
$\rho_{g},\rho_{f},\mu,\lambda\gg\alpha,\beta,\gamma$. In particular, we
choose

$\displaystyle\rho_{g},\mu\sim 10^{6}>\rho_{f}\sim 10^{5}$ $\displaystyle\gg
1\text{y}^{-1},$ (1) $\displaystyle\alpha,\beta,\gamma\sim 10^{[-4,-2]}$
$\displaystyle\ll\lambda\sim 1\text{y}^{-1}.$ (2)

So, fire spreading and extinction $\rho_{g}$, $\rho_{f}$, $\mu$ occurs on the
scale of hours, while grass regrows on ash over months ($\lambda$) and forest
spread, growth and mortality $\alpha$, $\beta$, $\gamma$ occurs over decades.
We take fire ignition rate $\phi\sim 1/N$ such that fires spontaneously occur
about once per year in the modelled area. Figure 1b–d shows a typical time
profile for fractions of cells in each state. The initial state is all grass
and the ignition rate $\phi$ is low, such that nearly full forest cover is the
only steady state (reached after $300$ years, Fig. 1h). Brief events of
rapidly spreading fire counteract a gradual spread of forest. Time scale
separation of forest dynamics (Fig. 1b), grass regrowth (Fig. 1c), and fire
spread (Fig. 1d) shows clearly.

Figure 2a shows a bifurcation diagram of steady state forest cover in the FGBA
cellular automaton, denoted by $[F]$ (see Eq. M1), versus fire ignition rate
$\phi$. Unstable steady states (saddles) were obtained by applying feedback
control to the simulations (see Methods). _Bistability_ occurs above a
critical ignition rate $\phi$. Simulations initiated at the saddle will tip up
or down depending on initial stochasticity (Fig. 2b). Near the lower end of
the bistability range, the saddle solution is fairly homogeneous, but for
higher $\phi$ values, a single hole of grass in forest arises (insets in Fig.
2a).

###  II.2 Fast and slow subprocesses

The time scale separation (Eqs. 1 and 2) permits treatment of the joint
vegetation and fire dynamics as a fast-slow system. Fire spread occurs on the
fast time scale, where the vegetation landscape is treated as constant. Forest
dynamics occurs on the slow time scale, where the effects of fire are a steady
state function of the vegetation landscape.

#### Fast process: fires spreading in a given landscape

On the time scale of a single fire event, forest dynamics are negligible
($\alpha,\beta,\gamma\ll 1/$day) such that we can consider the total landscape
of forest patches as fixed. For each ignition event, this results in the
following dynamics. A fire ignites on a grass cell, then spreads across its
grassland cluster at a rate $\rho_{g}$ per $BG$ link, after which it reaches
the interface with adjacent forest, where it starts intruding the forest at a
rate $\rho_{f}$ per $BT$ link. At any time, burning cells can spontaneously
extinguish, converting to ash at a rate $\mu$. The probabilities of fire
spreading into a neighbouring grass or forest cell before the originating cell
extinguishes are given by

$p_{g}:=\frac{\rho_{g}}{\rho_{g}+\mu}=0.9,\qquad
p_{f}:=\frac{\rho_{f}}{\rho_{f}+\mu}=0.1,$ (3)

where we have shown the chosen values in our simulations (adopted from Hébert-
Dufresne  _et al._ (2018)). Since new fires arise and grass regrows only with
much lower rates than the fire spreading rate ($\phi
N\lesssim\lambda\ll\rho_{f},\mu,\rho_{g}$), we observe repeated fire waves
well separated in time (Fig. 1b).

When a fire in grassland cluster with index $j$ reaches its interface with
adjacent forest, the resulting forest loss due to this single fire event can
be approximated as a fraction $p_{f}$ of the forest cells at this interface.
This leads to the approximate loss of forest per fire event

$\Delta_{F,j}^{\mathrm{loss}}=p_{f}[FG]_{j},$ (4)

where $[FG]_{j}$ counts the number of forest cells adjacent to grassland
cluster $j$ (with both sides of the equation optionally normalised by $N$).
This approximation relies on the assumptions that the fire reaches the whole
interface with forest (i.e. $p_{g}\to 1$) and only once per fire (i.e.
$\rho_{g}\gg\lambda\gg\phi N$), and that burning forest does not pass on fires
further to adjacent forest (i.e. $p_{f}$ small).

#### Slow processes: forest demography and fire damage

Forest demography and loss due to repeated fires occur on the slow time scale.
Writing the number of neighbour pairs with $FG$ as $[FG]$ (divided by $N$,
equivalently the domain perimeter of forest or grass patches, see Eq. M1), the
dynamics for tree recruitment and mortality result in an expected rate of
change for $[F]$

$\Delta_{F}^{\mathrm{gain}}:=\beta[G]-\gamma[F]+\alpha[FG].$ (5)

In Eq. 5, the rates of change are $\beta[G]$ for spontaneous forest growth on
grass, $\gamma[F]$ for forest mortality, and $\alpha[FG]$ for spread of forest
into grass at its perimeter.

The rate of forest erosion at its perimeter due to fire damage over many fire
events is the weighted sum over all grass clusters $j=1,...,n_{c}$, i.e.

$\displaystyle\Delta_{F}^{\mathrm{loss}}$
$\displaystyle:=\sum_{j=1}^{n_{c}}\phi N[G]_{j}\Delta_{F,j}^{\mathrm{loss}}$
$\displaystyle=\phi Np_{f}\sum_{j=1}^{n_{c}}[G]_{j}[FG]_{j},$ (6)

where $[G]_{j}$ is the fraction of $G$ cells in grass cluster $j$ (so,
$[G]=\sum_{j=1}^{n_{c}}[G]_{j}$), $\phi N[G]_{j}$ is the rate at which fires
spontaneously ignite in grass cluster $j$ ($\phi$ is the rate per cell and
$N[G]_{j}$ is the area of the cluster), and
$\Delta_{F,j}^{\mathrm{loss}}=p_{f}[FG]_{j}$ is the conversion of forest to
ash caused by each fire event (see Eq. 4) (note also that
$[FG]=\sum_{j=1}^{n_{c}}[FG]_{j}$). By defining the _grassland-weighted forest
perimeter_ as

$\langle[FG]\rangle_{\mathrm{cg}}:=\sum_{j}^{n_{c}}\frac{[G]_{j}}{[G]}[FG]_{j},$
(7)

Eq. 6 becomes

$\Delta_{F}^{\mathrm{loss}}:=\phi p_{f}N[G]\langle[FG]\rangle_{\mathrm{cg}}.$
(8)

The grassland-weighted forest perimeter $\langle[FG]\rangle_{\mathrm{cg}}$ is
the average perimeter of forest clusters weighted by the relative size of
their adjacent grass cluster.

###  II.3 Emergent slow dynamics

Figure 3: Rate of change according to Eq. 9. (a,b) Spatial snapshots at
indicated times, (c,d) time series of $[F]$, (e,f) time series of
$\Delta_{F}^{\mathrm{gain}}$, $\Delta_{F}^{\mathrm{loss}}$, (g,h) time series
of the right-hand side of Eq. 9 (gain minus loss). At $t{=}0$, the simulation
is started on the saddle (on the left, $[F](0)\approx 0.2$ and on the right,
$[F](0)\approx 0.7$, see blue and red circles in Fig. 2). Towards the left
(along $t^{+}$), a simulation that tips up and towards the right (along
$t^{-}$) a simulation that tips down is shown. Parameters are shown in Table
M1. Columns correspond to leftmost and rightmost vertical dashed lines in Fig.
2 ($\phi N=0.257$ and $\phi N=1.32$). Domain size: 100x100 cells.

We now form the balance between the slow processes discussed above, assuming
fire converts trees immediately to grass (i.e. $\lambda\gg\phi N$). The
resulting expected rate of forest cover change during a short time interval is

$\displaystyle\langle\frac{\mathrm{d}[F]}{\mathrm{d}t}\rangle=$
$\displaystyle\Delta_{F}^{\mathrm{gain}}-\Delta_{F}^{\mathrm{loss}},$
$\displaystyle\frac{\mathrm{d}[F]}{\mathrm{d}t}=$
$\displaystyle\beta[G]-\gamma[F]+\alpha[FG]-\phi
p_{f}N[G]\langle[FG]\rangle_{\mathrm{cg}},$ (9)

where we used Eqs. 5 and 8, and assumed on the left-hand side that $N$ is
sufficiently large, such that, via the law of large numbers,
$\langle\mathrm{d}[F]/\mathrm{d}t\rangle\approx\mathrm{d}[F]/\mathrm{d}t$.
Equation 9 can be understood intuitively as forest and grass competing for
space within patches (spontaneous terms) and at their interface (interaction
terms). A larger interface $[FG]$ leads simultaneously to faster forest spread
(proportional to its perimeter $[FG]$), and to increased exposure to fires
(proportional to its grassland-weighted perimeter
$\langle[FG]\rangle_{\mathrm{cg}}$). Fires are most damaging to forest when
$[G]$ forms a single cluster, i.e. $\langle[FG]\rangle_{\mathrm{cg}}=[FG]$,
such that each fire reaches the whole interface. Conversely, when forest
patches break $[G]$ into several clusters $\langle[FG]\rangle_{\mathrm{cg}}$
is smaller than $[FG]$, such that several ignitions are required to have the
same effect, slowing forest erosion down. Additionally, the total amount of
grass $N[G]$ determines the number of ignitions and hence the rate at which
grass spreads into forest. The parameters determine the relative weight of
each of the discussed effects.

Figure 3 shows example simulations along trajectories starting from the saddle
equilibria of Fig. 2, showing forest area $[F]$ in space and time (a–d), the
gain/loss terms $\Delta_{F}^{\mathrm{gain}}$ and $\Delta_{F}^{\mathrm{loss}}$
defined in Eqs. 5 and 8 (e–f), and the right-hand side of Eq. 9 (gain minus
loss, g–h). The left column of Fig. 3 shows simulations for low fire ignition
rate $\phi$ and low $[F](0)$, and the right column for high $\phi$ and high
$[F](0)$. Each column shows two realisations, both starting from the same
saddle steady state. One realisation evolves toward high forest cover, shown
on axis $t^{+}$ (increasing to the left from $t=0$), the other realisation
evolves toward low forest cover, shown on axis $t^{-}$ (increasing to the
right from $t=0$). In the stable steady states, gain (green) and loss (red)
terms vary around the same mean. On the saddle (at time $0$) gain and loss
functions cross. This shows that the steady states and changes are accurately
captured by Eq. 9. The largest changes in forest cover $[F]$ occur when there
are large changes in forest loss due to fire. The snapshots in Fig. 3a,b show
that the high-cover state changes as an expanding/contracting hole in the
forest, whereas the low-cover state is more homogeneous.

###  II.4 Emergent nonlinear relations

Figure 4: Emergent relations between key quantities and forest area $[F]$:
(a–c) forest perimeter $[FG]$ and grassland-weighted forest perimeter
$\langle[FG]\rangle_{\mathrm{cg}}$, (d–f) forest gain terms and loss terms in
Eqs. 5 and 6, (g–i) forest area rate of change $(\mathrm{d}/\mathrm{d}t)[F]$
from Eq. 9. Columns correspond to vertical dashed lines in Fig. 2 ($\phi
N{=}0.257,\phi N{=}0.38,\phi N{=}1.32$). Domain size: 100x100 cells.

Equation 9 explains how the rate of change of $[F]$ depends on the perimeter
quantities $[FG]$ and $\langle[FG]\rangle_{\mathrm{cg}}$. Figure 4a–c shows a
scatterplot of $[FG](t)$ and $\langle[FG]\rangle_{\mathrm{cg}}(t)$ versus
$[F](t)$ for three different values of $\phi$, and for an ensemble of
simulations starting from the saddle in Fig. 2a, with each point being a value
observed at a discrete observation time. Remarkably, we observe that $[FG]$
and $\langle[FG]\rangle_{\mathrm{cg}}$ lie on a narrow band around some steady
state functions $[FG]^{*}$ and $\langle[FG]\rangle_{\mathrm{cg}}^{*}$ of $[F]$
(and $\phi$), which implies that $[FG],\langle[FG]\rangle_{c_{g}}$ are
changing on a much faster time scale, making them slaved to $[F]$. Figure 4d–f
shows the terms on the right-hand side of Eq. 9 depending on $[F]$, splitting
between gain and loss terms
$\Delta_{F}^{\mathrm{gain}},\Delta_{F}^{\mathrm{loss}}$, as defined in Eqs. 5
and 8. Steady states occur when gain equals loss
($\Delta_{F}^{\mathrm{gain}}{=}\Delta_{F}^{\mathrm{loss}}$). The resulting
graph of $\mathrm{d}[F]/\mathrm{d}t$ versus $[F]$ in Fig. 4g–i clearly shows
the bistability of $[F]$.

Replacing the quantities $[FG]$ and $\langle[FG]\rangle_{\mathrm{cg}}$ by
their steady-state functions $[FG]^{*}$ and
$\langle[FG]\rangle_{\mathrm{cg}}^{*}$ results in the single-variable ODE for
$[F]$,

$\frac{\mathrm{d}[F]}{\mathrm{d}t}=\beta[G]-\gamma[F]+\alpha[FG]^{*}-\phi
p_{f}N[G]\langle[FG]\rangle_{\mathrm{cg}}^{*},$ (10)

where $[FG]^{*},\langle[FG]\rangle_{\mathrm{cg}}^{*}$ are functions of $[F]$
and $\phi$ (as shown in Fig. 4a-c), and $[G]=1-[F]$. With these functions
$[FG]^{*}$ and $\langle[FG]\rangle_{\mathrm{cg}}^{*}$, the observed
bistability is caused by a classic scalar double-well potential of the
gradient system Eq. 10. In this ODE, nonlinearities emerge due to the
equilibrium dependence of the interface on forest area (affecting $[FG]^{*}$
and $\langle[FG]\rangle_{\mathrm{cg}}^{*}$), due to the segmentation of grass
cells near and below the percolation threshold (affecting
$\langle[FG]\rangle_{\mathrm{cg}}^{*}$) and due to dependence of the ignition
rate on grass patch size (multiplying $\langle[FG]\rangle_{\mathrm{cg}}^{*}$
with $[G]$). In Fig. S3, we show the roots of Eq. 10 using a nonparameteric
fit of $[FG]^{*}([F];\phi)$ and $\langle[FG]\rangle_{c_{g}}^{*}([F];\phi)$.
These match well with the steady states obtained via control (dot-dashed red).

If there is only one connected component of grass cells, we have
$\langle[FG]\rangle_{\mathrm{cg}}=[FG]$, such that Eq. 10 simplifies to

$\frac{\mathrm{d}[F]}{\mathrm{d}t}=\beta[G]-\gamma[F]+(\alpha-\phi
p_{f}N[G])[FG]^{*}.$ (11)

For homogeneous initial conditions (i.e. $[F]$ is about the same in different
large subsections of the domain), this approximation is expected to be valid
for small $[F]$, where most grass cells belong to the giant component. Figure
S3 shows the resulting steady states of Eq. 11 as a function of $\phi$ and
$[F]$ when only using the fit of $[FG]^{*}([F];\phi)$ (dashed blue). The
approximation is good for landscapes with low forest cover ($[F]\lesssim
0.2$). Above $[F]\approx 0.2$, it fails because the grassland breaks up into
multiple clusters and fires are smaller than in case of a single cluster.
Figure 4a–c already indicated that the single-cluster approximation is
accurate for low forest cover since
$\langle[FG]\rangle_{\mathrm{cg}}\approx[FG]$ for low $[F]$ in the
scatterplots.

##  III Discussion

Using as high-dimensional reference model a continuous-time stochastic
cellular automaton, adapted from ref. Hébert-Dufresne  _et al._ (2018), that
includes reaction rules for both forest dynamics and fire percolation, we have
shown how low-dimensional nonlinear dynamics of tropical forest emerge from
the spatial configuration of forest patches. Taking the assumption on the fast
time scale that fire spreads perfectly on grass, reaching the whole perimeter
of adjacent forest, makes two macroscopic structural quantities relevant for
the evolution of forest area on the slow time scale (Eq. 9): the forest
perimeter $[FG]$ for gain from forest spread, and, the grassland-weighted
forest perimeter $\langle[FG]\rangle_{\mathrm{cg}}$ for loss induced by fire.
Simulations initiated at the unstable steady states allowed us to explore the
whole phase space and revealed that these quantities are evolving on a much
faster time scale than forest area, such that on the slow time scale and for
given parameters, they are equilibrium functions of forest area. Hence,
replacing the perimeter quantities with their steady state functions of forest
area reduces the evolution of forest area on the slow time scale to a one-
variable ODE, Eq. 10. In this emergent relation between forest change and
forest area, nonlinear terms arise due to spatially structured interactions at
the forest-grass interface whereas linear terms arise from spontaneous
transitions, which are independent of spatial structure.

In the rest of our discussion, we will compare our result to what is
obtainable by mean-field models, illustrate the role of finite-size effects in
bistability, and consider the implications for empirical tests.

#### Inadequacy of mean-field models

Mean-field models fail to capture the role played by forest perimeter and
grassland clusters because they neglect spatial structure (by the well-
mixedness assumption, point (_i_) in the introduction). The mean-field
approximation for the FGBA process is given in Section S1 (Supporting
Information). Mean-field model Eq. S3 assumes absence of pairwise
correlations, which permits substitutions such as $[FG]=4[F][G]$ in Section S1
and closure of the system of equations. Such a substitution of counts of
larger units by an algebraic expression in terms of counts of smaller units in
order to close a system of equations is known as a _closure formula_ Wuyts and
Sieber (2021). Neglect of structure leads to a well-known bias in the location
of steady states compared to the simulations (e.g. Keeling, 1999). This may be
mitigated by mean-field models of higher order, which consider units of larger
size than single cells, but for dynamics on the square lattice in regimes with
long-range dependence (such as induced by fire or forest spreading) this would
require an impractically large number of equations Wuyts and Sieber (2021),
defeating its purpose of low-dimensional approximation.

By exploiting time scale separation of fire and forest dynamics, we arrived at
Eq. 10, which respects the spatially structured nature of the dynamics, in
contrast to mean-field approximations. In our description, the steady state
functions $[FG]^{*}$ and $\langle[FG]\rangle_{\mathrm{cg}}^{*}$ of $[F]$ are
closure formulas to obtain Eq. 10 from Eq. 9. Here the closure is observed
emerging in simulations instead of derived from an a priori imposed assumption
of statistical independence (see also Gross and Kevrekidis (2008)). Closure at
low dimension is feasible due to time-scale separation between forest
perimeter quantities and forest area.

#### Finite sizes and bistability

In simple bistable chemical systems, it is known that bistability converges to
an abrupt phase transition in the thermodynamic limit ($N{\to}\infty$) Ge and
Qian (2009); Thiele  _et al._ (2019), at a value known as the Maxwell point
(e.g. Goldenfeld, 2018), making the macroscopic state of the system
deterministically dependent on the parameters (e.g. pressure or temperature).
With only forest and grass or provided that savanna and forest components are
sufficiently decoupled, the same behaviour occurs in spatial mean-field models
of tropical forest, with a front between forest and nonforest that is
deterministically dependent on environmental drivers Wuyts  _et al._ (2019).
We do not expect such determinism as $N{\to}\infty$ to arise in the FGBA
cellular automaton. The infinite FGBA cellular automaton possesses grass
clusters of arbitrary size, such that, even when assuming that fire spreads
instantly on grass, there will always be some parts of the forest perimeter
shielded from intruding fires by adjacent ash cells. Were it not for this
shielding effect, then there would be a deterministic dependence of the
dynamics on fire ignition rate away from the absorbing states: $\phi{=}0$
would lead to forest spread while $\phi{>}0$ would lead to forest extinction
(see Section S2, for the case $\beta{=}\gamma{=}0$). In reality, finite fire
spreading rates and, in particular, the effects of heterogeneity in time (e.g.
the influence of weather on spreading rates) or space impose stronger limits
on the reach of fires. Hence, while, due to the large size of a typical system
of molecules, the thermodynamic limit is usually of practical relevance in
context of materials, in case of ecosystems, _finite sizes_ and their effects
are the norm. In finite domains, both the FGBA system and scalar reaction-
diffusion equations show bistability due to critical patch sizes and
dependence on interfacial characteristics (e.g. Liehr, 2013; Goel  _et al._ ,
2020), but this correspondence requires further scrutiny. In realistic
scenarios, the amount of bistability then depends (besides the parameters) on
a ‘characteristic system size’, which in turn depends on e.g. fire spreading
and/or plant dispersal ranges. How this is affected by environmental
heterogeneity is subject of future work.

#### Implications for empirical tests

By taking the simplest possible model that was able to include fire
percolation and forest spread, we were able to obtain structural indicators of
bistable forests prone to fire percolation effects. These indicators are
theoretical predictions additional to bimodality that can be tested against
spatial data. We emphasise here that the perimeter quantities $[FG]$ and
$\langle[FG]\rangle_{\mathrm{cg}}$ do not have to be in equilibrium for Eq. 9
to be valid, such that it can also be tested in places where human impact
keeps them away from their natural equilibrium. Furthermore, they are in
principle easy to determine from observations. In the original treatment of
the FGBA cellular automaton Hébert-Dufresne  _et al._ (2018), a climate-
dependent equilibrium perimeter-area relation for _isolated_ forest patches in
low-cover simulations was found to be consistent with data of protected areas
in the Brazilian Cerrado. The perimeter-area relation at the patch level for
given parameters determines the perimeter-area relation at domain level
$[FG]^{*}$ versus $[F]$. Therefore, we expect this consistency with data to
hold at the domain level, at least for low cover areas. Inclusion of
additional realism, such as environmental heterogeneity, larger
neighbourhoods, non-lattice geometry, presence and dynamics of other tree
types and life stages, will lead to modified perimeter quantities, but these
can be derived and tested in a similar manner. We hope that the additional
opportunities that our theoretical predictions offer for falsification lead to
a faster rate of progress in the field. If this progress leads to theory with
sufficient reliability and validity, it may eventually function as a tool for
conservation of existing or restoration of lost tropical forests.

##  IV Methods

#### Details of the FGBA probabilistic cellular automaton

The FGBA probabilistic cellular automaton is a minimal spatial stochastic
process that models the joint dynamics of tropical vegetation and fire. It is
an adapted version of the BGT(A) model of ref. (Hébert-Dufresne  _et al._ ,
2018). The modifications compared to ref. (Hébert-Dufresne  _et al._ , 2018)
are: (i) it runs in continuous time, (ii) it includes a spontaneous forest
mortality rate $\gamma$, (iii) it optionally includes rainfall feedback (for
Fig. S1) (iv) species $T$ is labelled as $F$, consistent with other models of
tropical vegetation dynamics (Staver and Levin, 2012; Wuyts  _et al._ , 2019;
Patterson  _et al._ , 2021). Note that according to some definitions,
probabilistic continuous-time cellular automata are considered as interacting
particle systems. In general, when studying the stochastic dynamics of a
number $n$ of interacting species on a square lattice with $N$ cells, the
state of the system can be represented as

$X:=(X_{1},X_{2},...,X_{N}),$

where $X_{i}$ is the label of the species that occupies cell $i$. Each cell is
occupied by exactly one of four possible species: grass, forest, burning and
ash, with labels $G$, $F$, $B$ and $A$. Transitions between states (species)
occur in continuous time, resulting in a continuous-time Markov chain with
$n^{N}$ possible states. The reaction rules for transitions between states are
shown in Table M1, where spontaneous conversions are shown on the left and
conversions due to nearest-neighbour interactions on the right (see also Fig.
1a for the state transition diagram). The latter type of interaction occurs
over each four nearest neighbour connections of the indicated type. E.g., fire
will spread into a given grass cell with a rate $\rho_{g}$ for each burning
neighbour. For realistic time scales, our parameters satisfy Eqs. 1 and 2,
which were empirically justified in Hébert-Dufresne  _et al._ (2018). This
results in separated time scales for forest spread (slow) and fire spread
(fast). We borrow our notation from the moment closure literature (e.g.
Keeling, 1999; Kiss  _et al._ , 2017; Wuyts and Sieber, 2021), writing the
global fraction of species with label $x$ and the interface between species
with label $x$ with label $y$ respectively as

$[x]:=\frac{1}{N}\sum_{i}^{N}\delta_{x}(X_{i}),\quad[xy]:=\frac{1}{N}\sum_{i,j}^{N}\mathsf{A}_{ij}\delta_{x}(X_{i})\delta_{y}(X_{j}),$
(M1)

where both are normalised by $N$, $\delta$ is the Kronecker delta function
($\delta_{x}(y){=}1$ if $x{=}y$ and $0$ otherwise) and
$\mathsf{A}\in\\{0,1\\}^{N\times N}$ the adjacency matrix. We simulated the
cellular automaton via a Gillespie algorithm (Gillespie, 2007) and used a
domain of $N{=}100{\times}100$ ($N{=}200{\times}200$ in Fig. 1) cells with
periodic boundary conditions.

Table M1: Reactions and rates (y-1) in the FGBA cellular automaton. Spontaneous | Spread
---|---
$G\overset{\beta}{\to}F,$ | $\beta=2\cdot 10^{-4}$ | $GF\overset{\alpha}{\to}FF,$ | $\alpha_{\;}=3\cdot 10^{-2}$
$F\overset{\gamma}{\to}G,$ | $\gamma=2\cdot 10^{-2}$ | $AF\overset{\alpha}{\rightarrow}FF,$ | $\alpha_{\;}=3\cdot 10^{-2}$
$G\overset{\phi}{\to}B,$ | $\phi=[0,2]\cdot 10^{-4}$ | $GB\overset{\rho_{g}}{\to}BB,$ | $\rho_{g}=9\cdot 10^{6}$
$B\overset{\mu}{\to}A,$ | $\mu=10^{6}$ | $FB\overset{\rho_{f}}{\rightarrow}BB,$ | $\rho_{f}=1.11\cdot 10^{5}$
$A\overset{\lambda}{\to}G,$ | $\lambda=5$ |  |

#### Noninvasive feedback control

To study steady states regardless of their stability in a simulation, we apply
noninvasive feedback control, similar to what was used in (Sieber  _et al._ ,
2008; Schilder  _et al._ , 2015; Barton, 2017; Neville  _et al._ , 2018) for
tracking unstable vibrations and buckling states in mechanical experiments. To
study the dependence of equilibria of $[x]$ on one of the conversion rates
$r$, we treat the simulation like an ideal physical experiment and introduce
an artificial stabilizing feedback loop of the form

$r(t)=r_{0}+g([x](t)-[x]_{\mathrm{ref}})\mbox{.}$ (M2)

In Eq. M2, $r$ is one of the conversion rates of the lattice model. The
artificial feedback makes this rate time dependent by coupling it to the
species occurrence fraction $[x](t)$ through the relation (M2). The factor $g$
is called the feedback control gain and is problem specific. The property of
noninvasiveness means that the controlled simulations have the same equilibria
as regular simulations (Barton and Sieber, 2013; Sieber  _et al._ , 2014;
Renson  _et al._ , 2017). This implies that if one extracts the equilibrium
values of the controlled simulation $(r^{*},[x]^{*})$, one can use them to
plot a 1-parameter bifurcation diagram of the stochastic lattice model.

Figure M1: Feedback control applied to the cellular automaton without
spontaneous conversions ($\beta{=}\gamma{=}0$): (a) the unstable steady state
of the bifurcation diagram (dashed) was derived via feedback control by
letting $\phi$ be a function of $[F]$ (blue line) such that it is stabilised,
then obtaining $(\phi,[F])$ by averaging and repeating this for many
$[F]_{\mathrm{ref}}$ (with appropriate $g$), (b) a regular simulation with the
same $\phi$ value (solid grey in (a)) and starting from the final state of the
controlled simulation tips up or down depending on initial perturbations, (c)
snapshot of the domain at the saddle for the control indicated in (a) (black:
forest, white: grass). Other parameter values are shown in Table M1.

Figure M1 shows an example for which the controlled rate is the spontaneous
conversion rate $\phi$ for reaction $G\to B$ and the output variable $[x]$ is
total forest fraction $[F]$. The feedback Eq. M2, indicated in blue in Fig.
M1a, stabilises a steady state that is unstable in a regular simulation. This
can be seen in Fig. M1b, where the unstable steady state is first stabilised
via control, after which the control is removed and a regular simulation is
started with the effective rate (and initial conditions) obtained from the
controlled simulation. Depending on initial perturbations, the regular
simulation gets either attracted to the 100% forest state or to the low tree
cover state. When the controlled simulation is in equilibrium (steady part of
the blue curve in Fig. M1b), the steady state values of $[x]$ for
$x\in\\{F,G,B,A\\}$ are obtained via taking the time average, i.e.

$\overline{[x]}=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}[x](t)dt,$ (M3)

where $t_{0}$ is the time after which the dynamics have settled to a steady
state and $T$ the averaging time. If $n_{G\to B}$ is the number of ignition
events between $t=t_{0}$ and $t=t_{0}+T$, the steady state of $r$ is obtained
by calculating the mean ignition rate as $n_{G\to B}/T$ and dividing this by
the mean number of grass cells, such that

$\bar{\phi}=\frac{n_{G\to B}/T}{\overline{[G]}},$ (M4)

where $\overline{[G]}$ is obtained as in Eq. M3. When repeating this exercise
for many $[F]_{\mathrm{ref}}$ values (with for each an appropriately chosen
gain $g$), one can get multiple points on the unstable branch. Points on the
stable branches can be obtained with regular simulations. On the final
selection of points, we applied Gaussian process regression to obtain smooth
curves and used moving block bootstrapping (Kreiss and Lahiri, 2012) to obtain
confidence intervals. One of the advantages of applying control is that one
can obtain states for which one would have to wait prohibitively long in a
regular simulation due to their instability. The inset in Fig. M1b shows for
instance a snapshot of the unstable steady state obtained with the control
indicated in Fig. M1a.

## Acknowledgements

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) grants EP/N023544/1 and EP/V04687X/1. This project has
received funding from the European Union’s Horizon 2020 research and
innovation program under grant agreement 820970.

## References

  * Staver  _et al._ (2011) A. C. Staver, S. Archibald, and S. A. Levin, Science 334, 230 (2011).
  * Hirota  _et al._ (2011) M. Hirota, M. Holmgren, E. H. Van Nes, and M. Scheffer, Science 334, 232 (2011).
  * Lenton  _et al._ (2008) T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, and H. J. Schellnhuber, Proceedings of the National Academy of Sciences 105, 1786 (2008).
  * Scheffer  _et al._ (2015) M. Scheffer, S. R. Carpenter, V. Dakos, and E. H. Van Nes, Annual Review of Ecology, Evolution, and Systematics 46, 145 (2015).
  * Wuyts  _et al._ (2017) B. Wuyts, A. R. Champneys, and J. I. House, Nature Communications 8, 15519 (2017).
  * Hanan  _et al._ (2014) N. P. Hanan, A. T. Tredennick, L. Prihodko, G. Bucini, and J. Dohn, Global Ecology and Biogeography 23, 259 (2014).
  * Gerard  _et al._ (2017) F. Gerard, D. Hooftman, F. van Langevelde, E. Veenendaal, S. M. White, and J. Lloyd, Global Ecology and Biogeography 26, 854 (2017).
  * Staver and Levin (2012) A. C. Staver and S. A. Levin, The American Naturalist 180, 211 (2012).
  * Schertzer  _et al._ (2015) E. Schertzer, A. C. Staver, and S. A. Levin, Journal of Mathematical Biology 70, 329 (2015).
  * Touboul  _et al._ (2018) J. D. Touboul, A. C. Staver, and S. A. Levin, Proceedings of the National Academy of Sciences , 201712356 (2018).
  * Goel  _et al._ (2020) N. Goel, V. Guttal, S. A. Levin, and A. C. Staver, American Naturalist 195, 833 (2020).
  * Wuyts  _et al._ (2019) B. Wuyts, A. R. Champneys, N. Verschueren, and J. I. House, PLoS ONE 14, e0218151 (2019).
  * Patterson  _et al._ (2021) D. D. Patterson, S. A. Levin, C. Staver, and J. D. Touboul, SIAM Journal on Applied Dynamical Systems 19, 2682 (2021).
  * Hébert-Dufresne  _et al._ (2018) L. Hébert-Dufresne, A. F. Pellegrini, U. Bhat, S. Redner, S. W. Pacala, and A. M. Berdahl, Ecology Letters 21, 794 (2018).
  * Bond (2008) W. J. Bond, Annual Review of Ecology, Evolution, and Systematics 39, 641 (2008).
  * Hoffmann  _et al._ (2012) W. A. Hoffmann, E. L. Geiger, S. G. Gotsch, D. R. Rossatto, L. C. Silva, O. L. Lau, M. Haridasan, and A. C. Franco, Ecology Letters 15, 759 (2012).
  * Archibald  _et al._ (2009) S. Archibald, D. P. Roy, B. W. van Wilgen, and R. J. Scholes, Global Change Biology 15, 613 (2009).
  * Van Nes  _et al._ (2018) E. H. Van Nes, A. Staal, S. Hantson, M. Holmgren, S. Pueyo, R. E. Bernardi, B. M. Flores, C. Xu, and M. Scheffer, PloS One 13, e0191027 (2018).
  * Stauffer and Aharony (1994) D. Stauffer and A. Aharony,  _Introduction to percolation theory_ (CRC press, 1994).
  * Christensen and Moloney (2005) K. Christensen and N. R. Moloney,  _Complexity and Criticality_, Imperial College Press Advanced Physics Texts, Vol. 1 (Imperial College Press and Distributed, 2005).
  * Cardoso  _et al._ (2022) A. W. Cardoso, S. Archibald, W. J. Bond, C. Coetsee, M. Forrest, N. Govender, D. Lehmann, L. Makaga, N. Mpanza, J. E. Ndong, A. F. Koumba Pambo, T. Strydom, D. Tilman, P. D. Wragg, and A. C. Staver, Proceedings of the National Academy of Sciences of the United States of America 119 (2022), 10.1073/pnas.2110364119.
  * Spracklen  _et al._ (2018) D. V. Spracklen, J. C. Baker, L. Garcia-Carreras, and J. H. Marsham, “The effects of tropical vegetation on rainfall,”  (2018).
  * Veenendaal  _et al._ (2018) E. M. Veenendaal, M. Torello-Raventos, H. S. Miranda, N. M. Sato, I. Oliveras, F. van Langevelde, G. P. Asner, and J. Lloyd, New Phytologist 218, 153 (2018).
  * Van Langevelde  _et al._ (2003) F. Van Langevelde, C. A. Van De Vijver, L. Kumar, J. Van De Koppel, N. De Ridder, J. Van Andel, A. K. Skidmore, J. W. Hearne, L. Stroosnijder, W. J. Bond, H. H. Prins, and M. Rietkerk, Ecology 84, 337 (2003).
  * Accatino  _et al._ (2010) F. Accatino, C. De Michele, R. Vezzoli, D. Donzelli, and R. J. Scholes, Journal of Theoretical Biology 267, 235 (2010).
  * D’Odorico  _et al._ (2006) P. D’Odorico, F. Laio, and L. Ridolfi, The American naturalist 167 (2006), 10.1086/500617.
  * Baudena  _et al._ (2010) M. Baudena, F. D’Andrea, and A. Provenzale, Journal of Ecology 98, 74 (2010).
  * Beckage  _et al._ (2011) B. Beckage, L. J. Gross, and W. J. Platt, Ecological Modelling 222, 2227 (2011).
  * Hoyer-Leitzel and Iams (2021) A. Hoyer-Leitzel and S. Iams, Bulletin of Mathematical Biology 83 (2021), 10.1007/s11538-021-00944-x.
  * Archibald  _et al._ (2012) S. Archibald, A. C. Staver, and S. A. Levin, Proceedings of the National Academy of Sciences 109, 847 (2012).
  * Favier  _et al._ (2004) C. Favier, J. Chave, A. Fabing, D. Schwartz, and M. A. Dubois, Ecological Modelling 171, 85 (2004).
  * Wuyts and Sieber (2021) B. Wuyts and J. Sieber, “Automated derivation and closure of moment equations for multistate dynamics on networks,”  (2021), arXiv:2111.07643 .
  * Keeling (1999) M. J. Keeling, Proceedings of the Royal Society B: Biological Sciences 266, 859 (1999).
  * Gross and Kevrekidis (2008) T. Gross and I. G. Kevrekidis, Epl 82, 38004 (2008).
  * Ge and Qian (2009) H. Ge and H. Qian, Physical Review Letters 103, 148103 (2009).
  * Thiele  _et al._ (2019) U. Thiele, T. Frohoff-Hülsmann, S. Engelnkemper, E. Knobloch, and A. J. Archer, New Journal of Physics 21, 123021 (2019), arXiv:1908.11304 .
  * Goldenfeld (2018) N. Goldenfeld,  _Lectures on Phase Transitions and the Renormalization Group_ (CRC Press, 2018) pp. 1–394.
  * Liehr (2013) A. Liehr,  _Dissipative solitons in reaction diffusion systems_ , Vol. 70 (Springer, 2013).
  * Kiss  _et al._ (2017) I. Z. Kiss, J. C. Miller, and P. L. Simon,  _Interdisciplinary Applied Mathematics_, Interdisciplinary Applied Mathematics, Vol. 46 (Springer International Publishing, Cham, 2017) pp. 1–413.
  * Gillespie (2007) D. T. Gillespie, Annual Review of Physical Chemistry 58, 35 (2007).
  * Sieber  _et al._ (2008) J. Sieber, A. Gonzalez-Buelga, S. A. Neild, D. J. Wagg, and B. Krauskopf, Phys. Rev. Lett. 100, 244101 (2008).
  * Schilder  _et al._ (2015) F. Schilder, E. Bureau, I. F. Santos, J. J. Thomsen, and J. Starke, Journal of Sound and Vibration 358, 251 (2015).
  * Barton (2017) D. Barton, Mechanical Systems and Signal Processing 84, 54 (2017).
  * Neville  _et al._ (2018) R. M. Neville, R. M. Groh, A. Pirrera, and M. Schenk, Physical review letters 120, 254101 (2018).
  * Barton and Sieber (2013) D. A. W. Barton and J. Sieber, Phys. Rev. E 87, 052916 (2013).
  * Sieber  _et al._ (2014) J. Sieber, O. E. Omel’chenko, and M. Wolfrum, Phys. Rev. Lett. 112, 054102 (2014).
  * Renson  _et al._ (2017) L. Renson, D. Barton, and S. Neild, International Journal of Bifurcation and Chaos 27, 1730002 (2017).
  * Kreiss and Lahiri (2012) J. P. Kreiss and S. N. Lahiri, in  _Handbook of Statistics_, Vol. 30 (2012) pp. 3–26.
  * Tomé and De Oliveira (2015) T. Tomé and M. J. De Oliveira,  _Stochastic dynamics and irreversibility_ (Springer, 2015) p. 394.

Supplement

##  S1 Mean-field equations of the FGBA process on the square lattice

Here we show how to obtain the mean-field equations. Note that ref. (Hébert-
Dufresne  _et al._ , 2018) already analysed the dynamics of the mean-field
model (Eq. S3) with $\gamma=k=0$. Writing $x_{i}$ as shorthand for
$\delta_{x}(X_{i})$ (equalling $1$ if $X_{i}=x$ and $0$ otherwise) and taking
expectations in each cell $i$, we obtain equations for the rate of change of
the expectation that cell $i$ is occupied by species $x\in\\{G,F,B,A\\}$,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}G_{i}\bigr{\rangle}$
$\displaystyle=\lambda\bigl{\langle}A_{i}\bigr{\rangle}-\bigl{\langle}G_{i}\Bigl{[}\beta+\phi+\sum_{j\in{\cal
N}(i)}(\rho_{g}B_{j}+\alpha
F_{j})\Bigr{]}\bigr{\rangle}+\gamma\bigl{\langle}F_{i}\bigr{\rangle},$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}F_{i}\bigr{\rangle}$
$\displaystyle=\beta\bigl{\langle}G_{i}\bigr{\rangle}+\alpha\bigl{\langle}(A_{i}+G_{i})\sum_{j\in{\cal
N}(i)}F_{j}\bigr{\rangle}-\rho_{f}\bigl{\langle}F_{i}\sum_{j\in{\cal
N}(i)}B_{j}\bigr{\rangle}-\gamma\bigl{\langle}F_{i}\bigr{\rangle},$ (S1)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}B_{i}\bigr{\rangle}$
$\displaystyle=\phi\bigl{\langle}G_{i}\bigr{\rangle}-\mu\bigl{\langle}B_{i}\bigr{\rangle}+\bigl{\langle}(\rho_{f}F_{i}+\rho_{g}G_{i})\sum_{j\in{\cal
N}(i)}B_{j}\bigr{\rangle},$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}A_{i}\bigr{\rangle}$
$\displaystyle=-\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}G_{i}\bigr{\rangle}-\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}F_{i}\bigr{\rangle}-\frac{\mathrm{d}}{\mathrm{d}t}\bigl{\langle}B_{i}\bigr{\rangle}=\mu\bigl{\langle}B_{i}\bigr{\rangle}-\lambda\bigl{\langle}A_{i}\bigr{\rangle}-\alpha\bigl{\langle}A_{i}\sum_{j\in{\cal
N}(i)}F_{j}\bigr{\rangle}.$

One may think of $\langle\cdot\rangle$ as ensemble averages here. This set of
equations can be derived rigorously from the master equation (e.g. Tomé and De
Oliveira, 2015). Summing over $i$, dividing by $N$ and using Eq. M1, we obtain
equations for total frequencies of all four species in the domain (relative to
$N$):

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle[G]\rangle$
$\displaystyle=\lambda\langle[A]\rangle-(\beta+\phi)\langle[G]\rangle-\rho_{g}\langle[GB]\rangle-\alpha\langle[FG]\rangle+\gamma\langle[F]\rangle,$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle[F]\rangle$
$\displaystyle=\beta\langle[G]\rangle+\alpha\langle[FA]\rangle+\alpha\langle[FG]\rangle-\rho_{f}\langle[FB]\rangle-\gamma\langle[F]\rangle,$
(S2) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle[B]\rangle$
$\displaystyle=\phi\langle[G]\rangle-\mu\langle[B]\rangle+\rho_{f}\langle[FB]\rangle+\rho_{g}\langle[GB]\rangle,$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle[A]\rangle$
$\displaystyle=\mu\langle[B]\rangle-\lambda\langle[A]\rangle-\alpha\langle[FA]\rangle,$

where $[xy]$ is the total number of neighbouring $xy$ pairs divided by $N$,
which will later be referred to as the $xy$ interface or $xy$ perimeter.
Section S1 is not a closed system. To close the system, we need to determine
all undetermined terms $[xy]$ on the right-hand side without creating new
unknowns. The simplest way to do this is to assume absence of pairwise
correlations, i.e. $\langle[xy]\rangle=4\langle[x]\rangle\langle[y]\rangle$.
We take the additional assumption of $N\to\infty$, such that the law of large
numbers applies and $[x]\to\langle[x]\rangle$. Taking these assumptions leads
to the _simple mean-field approximation_ of the cellular automaton model:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}[G]$
$\displaystyle=\lambda[A]-(\beta+\phi)[G]-4\rho_{g}[G][B]-4\alpha[F][G]+\gamma[F],$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}[F]$
$\displaystyle=\beta[G]+4\alpha[F][A]+4\alpha[F][G]-4\rho_{f}[F][B]-\gamma[F],$
(S3) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}[B]$
$\displaystyle=\phi[G]-\mu[B]+4\rho_{f}[F][B]+4\rho_{g}[G][B],$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}[A]$
$\displaystyle=\mu[B]-\lambda[A]-4\alpha[F][A],$

In the main text we focus exclusively on the fire-vegetation feedback due to
percolation effects but here we also test the effect of a rainfall feedback
that works globally, by taking a linear dependence of forest spreading rate on
total forest area, i.e.

$\alpha([F])=\alpha_{0}+k[F].$ (S4)

The mean field assumes spatial homogeneity, well-mixedness and large system
size. These assumptions are valid when the dynamics occur on a large complete
graph and interactions occur at contact rates of order $1/N$.

Figure S1: FGBA process (Table M1 for simulations and Eq. S3 for mean field):
(a) mean-field steady states of forest cover $[F]$ as a function of fire
ignition rate $\phi$ and rainfall feedback strength $k$ (blue: steady state
manifold, green: contours at fixed $k$ values, red: saddle-node bifurcations),
(b) steady states of $[F]$ in the simulations versus $\phi$ (see legend for
parameter values, shading: two-standard deviation confidence interval of the
mean). The curves in (b) were obtained via a Gaussian process regression on
the output of a finite number of simulations (regular for the stable and
controlled for the unstable branches).

The FGBA cellular automaton exhibits bistability of tree cover (first shown by
ref. (Hébert-Dufresne  _et al._ , 2018)) in ranges of all parameters that are
expected to show considerable spatial heterogeneity in a given ecosystem:
$\alpha,\beta,\gamma,\phi$. The bistability range is widened by the rainfall
feedback, as shown in Fig. S1b (red versus blue curve). Yet, for the mean-
field model at low rainfall feedback gain $k$, bistability can also be
decreased (Fig. S1a). Hence, the effect may not be monotonic. That bistability
of tree cover occurs even without rainfall feedback shows that both fire and
rainfall feedbacks can cause it independently. Note that the mean field shows
significant bias.

##  S2 Evolution of fronts — heterogeneous states

Figure S2: Perimeter quantities and change rates according to Eq. 9 when
initial conditions are heterogeneous and there are no spontaneous transitions
($\beta{=}\gamma{=}0$) for $\phi{=}6.98\cdot 10^{-5}$, displayed as in Figs. 3
and 4. (a) Snapshots of the domain at indicated times. (b,d,f) Cover,
interface and rate of change from Eq. 9 versus time. (c,e,g) Interface, grass-
cluster weighted average of the interface, and rate of change from Eq. 9
versus cover. Remaining parameters are
$\alpha{=}0.03,\rho_{g}{=}10^{10},\rho_{f}{=}1.11\cdot
10^{5},\mu{=}10^{6},\lambda{=}5$.

Here, we illustrate the case where grass and forest are initially separated
into two contiguous areas with their interface extending along a straight line
Because for this type of initial conditions the single-cluster approximation
(Eq. 11) is valid, we can focus on the evolution of the interface. As
spontaneous conversion between forest and grass (with rates $\beta$ and
$\gamma$) increases mixing and promotes homogeneity at large scale, we expect
the effects of heterogeneous initial conditions to be most persistent when the
spontaneous conversion rates $\beta$ and $\gamma$ are small. Therefore, we
will set $\beta=\gamma=0$, for which Eq. 9 becomes

$\frac{\mathrm{d}[F]}{\mathrm{d}t}=(\alpha-\phi p_{f}N[G])[FG].$ (S5)

Hence, the precise shape of the interface $[FG]$ does not affect the location
of the steady states, only the rate at which they are approached or receded
from. The trivial steady states of Eq. S5 (for $[FG]=0$) are $[F]=0$ and
$[F]=1$, which are stable, and between them, there is the (unstable) saddle

$[F]^{*}=1-\frac{\alpha}{\phi Np_{f}}.$ (S6)

As seen in Fig. S4, this analytical prediction (solid black) matches the
controlled simulations with $p_{g}=0.9999$ (shaded blue). For $p_{g}=0.9$
(shaded red), which we used before, there is a small bias. In the limit of
$N\to\infty$, Eq. S6 converges to $[F]^{*}_{\infty}=1$, implying that in an
infinite domain, any positive fire rate leads to extinction of forest below
$[F]^{*}=1$. When ignoring the effect of ash, this would also occur for
heterogeneous initial conditions. That is, considering an infinite domain with
many grass clusters of which the size is a random variable (with support
$[0,\infty)$), there will be initial grass clusters of arbitrarily large size,
which will expand and eventually drive forest extinct. However, such
determinism does not occur in the simulations because at high fire rates,
patches with ash start to block fires, and the rate of exposure of the forest
interface to fire becomes limited by the rate at which ash is converted back
into grass. As a first correction for this, one can multiply $p_{f}$ with the
average proportion of grass sites that are in the ash state after the expected
waiting time between fires $1-\exp(-\lambda/\phi N)$, such that
$[F]^{*}_{\infty}=1-\alpha/\lambda p_{f}$. Keeping in mind that we are
focusing on heterogeneous states, the analysis here implies that for
$\gamma=\beta=0$, there is a critical patch size above which the forest patch
expands and below which it contracts. The intuition is that above the critical
forest patch size, there is not enough grass area to reach the minimum number
of ignitions required to erode the forest.

In Fig. S2, we show how the dynamics and steady states arise from $[FG]$, as
we did in Fig. 3 and Fig. 4, but now starting with separated patches of grass
and forest that interface on a line on both sides (showing for $[F](0)=0.57$).
From $[FG](0)=2L/N=2N^{-1/2}=0.02$, the interface quickly gains roughness due
to the dynamics, until it reaches a steady state of $[FG]^{*}\approx
6L/N=6N^{-1/2}=0.06$ (see Fig. S2a–c), which appears independent of $\phi$
(not shown). Figure S2b confirms that
$\langle[FG]\rangle_{\mathrm{cg}}\approx[FG]$ (except when forest cover
approaches $[F]=0$ or $[F]=1$), confirming that the single-cluster
approximation is valid. Away from $[F]=0$ and $[F]=1$, $[FG]$ stays about
constant as the forest interface moves (Fig. S2b–c). Therefore the gain and
loss terms (as defined in Eqs. 5 and 8) are now, respectively, constant and
linearly decreasing with $[F]$ (Fig. S2e), such that
$\mathrm{d}[F]/\mathrm{d}t$ increases linearly with $[F]$ (Fig. S2g), except
near $[F]=0$ and $[F]=1$, where it connects to 0, as here $[FG]=0$.

##  S3 Supplementary figures

Figure S3: Steady states of Eq. 10 (multicluster – dot-dashed red) and of Eq.
11 (single cluster – dashed blue) compared to controlled simulations (solid
black with shading indicating two-standard deviation confidence interval of
the mean). Parameters in the cellular automaton:
$\gamma{=}0.02,\alpha{=}0.03,\beta{=}2\cdot 10^{-4},\rho_{g}{=}9\cdot
10^{6},\rho_{f}{=}1.11\cdot 10^{5},\mu{=}10^{6},\lambda{=}5$. Figure S4:
Saddle of Eq. S5 via Eq. S6 (single cluster – solid black) compared to
controlled simulations (shaded red: $p_{g}=\rho_{g}/(\rho_{g}+\mu)=0.9$,
shaded blue: $p_{g}=0.9999$) in case of heterogeneous initial conditions and
without spontaneous interactions ($\beta{=}\gamma{=}0$). Other parameters:
$\alpha{=}0.03,\rho_{f}{=}1.11\cdot 10^{5},\mu{=}10^{6},\lambda{=}5$.