\section{Instructions for RS \LaTeX\ Template}

This latex class file is available for authors to prepare the manuscript for the Royal Society journals.
It is assumed that the authors are familiar with either plain \TeX, \LaTeX,\ \AmS-\TeX\ or a standard latex set-up,
hence only the essential points are described in this document. For more details please see the \textit{\LaTeX\ User's Guide}
or \textit{The not so short introduction to \LaTeXe}.

\section{Installation}
Within the supplied set of files, \texttt{rsproca\_new.cls} need to be copied into a directory
where tex looks for input files. The other files need to be kept as a reference while preparing your manuscript.
Please use pre-defined commands from \texttt{RSPA\_Author\_tex.tex} for title, authors, address, abstract, keywords, body etc.

\section{How to start using rsproca\_new.cls}
Before you type anything that actually appears in the paper you need to include a \verb+\documentclass{rsproca_new}+ command
at the very beginning and then, the two commands that have to be part of any latex document, \verb+\begin{document}+
at the start and the \verb+\end{document}+ at the end of your paper. The main structure of your document should be as follows:

\begin{verbatim}
\documentclass{rsproca_new}
\begin{document}
\title{...}
\author{....}
\address{...}
\subject{...}
\keywords{...}
\corres{...}

\begin{abstract}
..........
..........
\end{abstract}

\begin{fmtext}
............
\end{fmtext}

\maketitle
....
\section{....}
...
\subsection{....}
....
\end{document}
\end{verbatim}

Note: In order to insert few lines on first page use the command \verb+\begin{fmtext}....\end{fmtext}+ before \verb+\maketitle+ command.
Please look into sample.tex for more clarification.

\section{Preamble}

In the preamble portion i.e. before \verb+\begin{document}+ insert the respective journal name:
\verb+\jname{rspa}+ and \verb+\Journal{Proc R Soc A\ }+ to get journal name in running heads.

\section{Packages used by the class file}

There are some packages that are essential when using the class file:
\begin{verbatim}
amsmath amssymb amsfonts amsthm graphicx endfloat endnotes setspace
verbatim geometry times helvet courier mathtime bm url babel dcolumn
\end{verbatim}

\noindent Some commonly-used packages are already used by this class file:
\begin{verbatim}
xspace amscd rotating latexsym multicol array algorithm subfigure
\end{verbatim}

\noindent Note:

For ethics section use the command \verb+\ethics{}+.

For dataccess section use the command \verb+\dataccess{}+.

For authors contribute section use the command \verb+\aucontribute{}+.

For competing section use the command \verb+\competing{}+.

For funding section use the command \verb+\funding{}+.

For acknowledgement section use the command \verb+\ack{}+.

For disclaimer section use the command \verb+\disclaimer{}+.

\vskip2.5pc

\noindent The rsproca\_new file is intended as a guide only. Article lengths based
on this estimation may be subject to change when the article is prepared for publication.
Authors whose article is estimated close to the page charge limit should contact the office regarding page charges.

\end{document}

\section{Instructions for RS \LaTeX\ Template}

This latex class file is available for authors to prepare the manuscript for the Royal Society journals.
It is assumed that the authors are familiar with either plain \TeX, \LaTeX,\ \AmS-\TeX\ or a standard latex set-up,
hence only the essential points are described in this document. For more details please see the \textit{\LaTeX\ User's Guide}
or \textit{The not so short introduction to \LaTeXe}.

\section{Installation}
Within the supplied set of files, \texttt{rsproca\_new.cls} need to be copied into a directory
where tex looks for input files. The other files need to be kept as a reference while preparing your manuscript.
Please use pre-defined commands from \texttt{RSPA\_Author\_tex.tex} for title, authors, address, abstract, keywords, body etc.

\section{How to start using rsproca\_new.cls}
Before you type anything that actually appears in the paper you need to include a \verb+\documentclass{rsproca_new}+ command
at the very beginning and then, the two commands that have to be part of any latex document, \verb+\begin{document}+
at the start and the \verb+\end{document}+ at the end of your paper. The main structure of your document should be as follows:

\begin{verbatim}
\documentclass{rsproca_new}
\begin{document}
\title{...}
\author{....}
\address{...}
\subject{...}
\keywords{...}
\corres{...}

\begin{abstract}
..........
..........
\end{abstract}

\begin{fmtext}
............
\end{fmtext}

\maketitle
....
\section{....}
...
\subsection{....}
....
\end{document}
\end{verbatim}

Note: In order to insert few lines on first page use the command \verb+\begin{fmtext}....\end{fmtext}+ before \verb+\maketitle+ command.
Please look into sample.tex for more clarification.

\section{Preamble}

In the preamble portion i.e. before \verb+\begin{document}+ insert the respective journal name:
\verb+\jname{rspa}+ and \verb+\Journal{Proc R Soc A\ }+ to get journal name in running heads.

\section{Packages used by the class file}

There are some packages that are essential when using the class file:
\begin{verbatim}
amsmath amssymb amsfonts amsthm graphicx endfloat endnotes setspace
verbatim geometry times helvet courier mathtime bm url babel dcolumn
\end{verbatim}

\noindent Some commonly-used packages are already used by this class file:
\begin{verbatim}
xspace amscd rotating latexsym multicol array algorithm subfigure
\end{verbatim}

\noindent Note:

For ethics section use the command \verb+\ethics{}+.

For dataccess section use the command \verb+\dataccess{}+.

For authors contribute section use the command \verb+\aucontribute{}+.

For competing section use the command \verb+\competing{}+.

For funding section use the command \verb+\funding{}+.

For acknowledgement section use the command \verb+\ack{}+.

For disclaimer section use the command \verb+\disclaimer{}+.

\vskip2.5pc

\noindent The rsproca\_new file is intended as a guide only. Article lengths based
on this estimation may be subject to change when the article is prepared for publication.
Authors whose article is estimated close to the page charge limit should contact the office regarding page charges.

\end{document}

\section{Introduction}

Tipping points, or critical transitions, are instabilities known to occur in natural and human systems subjected to changing external conditions, or external inputs. They may be explained in layman's terms as a sudden and large change in the state of a system in response to a small or slow change in the external input. The change in the state of the system may be permanent or transient.

Climate change presents an immediate and urgent challenge to the global community, and is an important factor for tipping points in natural systems. An immense amount of research is being conducted to predict and prevent its worst effects.
Hand in hand with this research effort,
interest in tipping
has accelerated due to theorised~\cite{Lenton1786,ritchie_overshooting_2021} and observed~\cite{caesar_rahmstorf_robinson_feulner_saba_2018, boers_rypdal2021} tipping points in  the  earth system caused by climate change.
Many of the tipping elements  identified in~\cite{Lenton1786}, for example the loss of  Arctic sea ice  or the shutdown of the Atlantic Meridional Overturning Circulation, can be captured by elaborate, high-resolution mathematical models referred to as General Circulation Models (GCMs). Additionally, there are tipping elements in the earth system that are not captured by GCMs, such as Amazon rainforest dieback~\cite{cox_betts_collins_harris_huntingford_jones_2004,boulton_lenton_boers_2022},
that may have an equally strong global impact. In this paper, we use a
conceptual
soil-carbon model
with time-varying climate as an external input to describe
one such tipping element, namely a release of gigatonnes of carbon from temperature-sensitive peat soils into the atmosphere.

Owing to the explicit time dependence of the external input, the ensuing dynamical system is nonautonomous.
This means that analysis of tipping points requires, in general, techniques beyond classical autonomous stability theory~\cite{wieczorek_excitability_2011,ashwin2017,alkhayuon_2018,kiers_jones_2019,okeeffe_tipping_2020,longo_2021,kuehn_longo_2022,alkhayuon2021,Wieczorek2021R}.
Nonetheless,
it is useful to consider the corresponding {\em autonomous frozen system} with fixed-in-time inputs. In the frozen system we identify a desired stable state, and refer to this state as the {\em base state}.
When the external input changes over time,  the shape and position of the base state may change too, and the nonautonomous system will try
to {\em track} the moving base state. However, in some cases tracking is not possible and tipping occurs.
For example, the base state may lose stability or disappear in a classical bifurcation at some {\em critical level} of the input. If this bifurcation is {\em dangerous}~\cite{thompson_sieber_2011},
the system tips to a different stable state, referred to as an {\em alternative stable state}~\cite{scheffer_carpenter_foley_folke_walker_2001,SCHEFFER1993275}.
We then say the system undergoes {\em bifurcation-induced tipping}, or in short {\em B-tipping}~\cite{thompson_sieber_2011,kuehn2011,ashwin2012}.
Another, arguably more interesting example is when the external input changes faster than some {\em critical rate}, the nonautonomous system deviates too far from the changing base state, crosses some {\em threshold}~\cite{fitzhugh_1955,Wieczorek2021R}
or some {\em quasithreshold}~\cite{fitzhugh_1955,wieczorek_excitability_2011,Mitry2013,perryman_adapting_2014}
and tips to an {\em alternative state}. If the base state never loses stability or disappears, such tipping is caused entirely by the {\em rate of change} of the external input. Therefore,  we say the system undergoes {\em rate-induced tipping}, or in short {\em R-tipping}~\cite{scheffer_2008,luke_soil_2011,wieczorek_excitability_2011,ashwin2012,perryman_adapting_2014,siteur_2014,siteur_2016,ashwin2017,alkhayuon_2018,alkhayuon_2019,Vanselow_2019,kiers_jones_2019,okeeffe_tipping_2020,kiers_2020,arumugam_2020,suchithra_2020,neijnens_2021,longo_2021,vanselow_2021,lohmann_2021,lohmann_ditlevsen_2021,pierini_ghil_2021,Merker2021,kuehn_longo_2022,kaur_sharathi_dutta_2022,Wieczorek2021R}.
Crucially, unlike B-tipping, R-tipping can occur to an {\em  alternative transient state} that lasts for a finite time\footnote{This is in contrast to an {\em alternative stable state} that lasts forever.}, after which the system returns to the base state~\cite{wieczorek_excitability_2011,Mitry2013,Vanselow_2019}. Such R-tipping is referred to as ``reversible" in~\cite{Wieczorek2021R}. Systems that exhibit reversible R-tipping are also known as {\em excitable systems}~\cite{Rinzel89,Izhikevich2006,wieczorek_excitability_2011,Wechselberger2013}.
The peat soil instability studied in this paper is an example of
{\em reversible R-tipping} to an {\em alternative metastable state} that lasts for a long but finite time, that occurs due to crossing an elusive quasithreshold.

The main motivation for our study is a combination of two environmental  features of the Arctic. First is the organic carbon content.
Estimates for the organic carbon contained in Northern and Arctic permafrost peat soil alone range from approximately 500Gt~\cite{bg-9-4071-2012} to approximately 1700Gt~\cite{tarnocai_soil_2009,Rapid_expansion_of_northern_peatlands,Spatial_heterogeneity_soil_carbon}. For comparison, the atmospheric carbon pool is estimated to be approximately 850 Gt~\cite{ClimateChange_ThePhysicalScienceBasis}.
Second is the rate of atmospheric warming and the increasing trends in summer heatwaves.
Due to so called ``Arctic Amplification''~\cite{screen_simmonds_2010},  both the so-far observed and future predicted warming for Arctic regions are approximately double the global mean~\cite{overland_urgency_2019}.
This means that the Arctic is home to massive deposits of ancient peat carbon and is the fastest warming region on the planet.
To visualise this combination of environmental features we combine in figure~\ref{fig:RateMap} (colourscale) recent rates of global warming\footnote{Note that these recent short-term rates of global warming exceed the long-term rates in the CMIP5 outputs.} obtained from the Berkeley Earth~\cite{berkely_earth_data,Berkely_Earth_Temp_generation} global temperature dataset and (greyscale) the global distribution of peat soils obtained from the PEATMAP dataset~\cite{PEATMAP2018}.
In addition to the increase in the mean global temperature, there is an increase in the intensity, frequency and duration of summer heatwaves in the Northern Hemisphere~\cite{met_office_2022,perkins-kirkpatrick_lewis_2020}, with the Arctic temperature record a scorching $38^\circ$C in 2020~\cite{united_nations_2021}.
Since peat soils are bioactive and thus temperature sensitive~\cite{Makiranta2009}, such conditions
can lead to a tipping point. This is the reason why northern latitude peat-soils were identified as a potential tipping element in~\cite{Lenton1786}.
To the best of our knowledge, the first example of R-tipping in peat soils
triggered by
atmospheric warming was reported, but not emphasised, by Khvorostyanov et al.~\cite[Fig.4(a)]{khvorostyanov_vulnerability_2008-1}.
Later, Luke and Cox~\cite{luke_soil_2011} proposed a conceptual soil-carbon model that exhibits a short-lived explosive release of soil carbon into the atmosphere above some critical rate of global warming, which they termed the ``compost-bomb instability"; see also~\cite{Clarke2021}. The dynamical mechanism responsible for this instability was explained by Wieczorek et al.~\cite{wieczorek_excitability_2011}.
Additionally, it has been known that spontaneous combustion is the main cause of fires at composting facilities~\cite{Buggeln2002}, which can then spread, e.g. the recent Wennington fire in London~\cite{sawer_2022,warburton_thrower_2022}.

\begin{figure}[t]
    \centering
    \includegraphics[]{RateMapv14.pdf}
    \caption{
     (Colourscale)  Historical global rates of warming from the period 1951-1980 to the period 2017-2021  obtained using  observations compiled by Berkeley Earth~\cite{berkely_earth_data,Berkely_Earth_Temp_generation}, together with (black and gray) the global distribution of peat soils obtained from PEATMAP~\cite{PEATMAP2018}. Note the largest warming rates at higher latitudes, in areas with a significant concentration of peat soils, such as (blue arrow) Cherskii in Northern Siberia.
    }
    \label{fig:RateMap}
\end{figure}

This paper is organised in two parts.
In the first, mathematical modelling part of the paper in Sections~\ref{sec: rtscm} and~\ref{sec:RHMSZF}, we:
\begin{itemize}
    \item
    Modify the  conceptual soil-carbon model introduced by Luke and Cox~\cite{luke_soil_2011} with a more realistic microbial soil respiration function.
    \item
    Show that the modified  conceptual  model has a new {\em alternative hot metastable state} and reproduces the key features of the medium-complexity model from~\cite[Fig.4(a)]{khvorostyanov_vulnerability_2008-1}.
    \item
    Demonstrate R-tipping to the hot metastable state for realistic climate change scenarios including global warming and summer heatwaves.
    \item
    Based on the above, propose an explanation for ``Zombie fires" in peatlands~\cite{BBC2021,Scholten2021,mccarty_arctic_2020,vaughan_arctic_2020} that disappear from the surface, smoulder underground
     during the winter, and ``come back to life" in the spring.
\end{itemize}
In the second, mathematical analysis part of the paper in Sections~\ref{sec: mtc}--\ref{sec:HSA}, we identify non-obvious dynamical mechanisms that are responsible for the R-tipping instability to the hot metastable state. The main obstacles to
analysis of this instability are twofold: the soil-carbon model is a non-autonomous dynamical system, so it does not
have any equilibria or compact invariant sets, and the R-tipping is due to crossing an elusive quasithreshold in the phase space~\cite{Fitzhugh_1961,wieczorek_excitability_2011,Mitry2013,perryman_adapting_2014,Vanselow_2019,vanselow_2021,kaur_sharathi_dutta_2022,arnscheidt_2022}.
To overcome these obstacles we combine three different strategies.
Firstly, we consider external inputs that decay to a constant at infinity. Secondly, we compactify the problem to include the equilibrium base states
for the autonomous limit systems from infinity~\cite{Wieczorek_2021}.
These two strategies alone work for R-tipping due to crossing  regular thresholds since these thresholds are anchored at infinity by  unstable compact invariant sets called {\em regular R-tipping edge states}~\cite{Wieczorek2021R}. However,  quasithresholds do not contain such edge states~\cite[Sec.8(d)]{Wieczorek2021R}. Therefore, thirdly, we build on~\cite{wieczorek_excitability_2011,perryman_adapting_2014}, take advantage of large timescale separation in the compactified soil-carbon  model and apply concepts from geometric singular perturbation theory~\cite{fenichel_1979,ckrt_1995,Wechselberger2013,kuehn_2016}.
Specifically, we:
\begin{itemize}
    \item
    Define R-tipping due to crossing a quasithreshold in terms of {\em slow manifolds} and {\em canard trajectories}~\cite{Wechselberger2013}
    in the autonomous compactified system.
    \item
    Identify {\em singular R-tipping edge states}: special points called {\em folded singularities}~\cite{SZMOLYAN2001419} that arise in the slow (reduced) system, and {\em new saddle equilibria} that arise in the fast (layer) system.
    \item
    Reduce an R-tipping problem due to crossing a quasithreshold to a heteroclinic orbit connecting the base state for the past limit system to a singular R-tipping edge state. We use this reduction to identify four different cases of R-tipping
    as follows:
    \begin{itemize}
        \item
        For global warming, we identify three
        %
        (slow) cases of R-tipping
        via unfolding of a  codimension-two {\em  heteroclinic folded saddle-node type-I singularity} that arises in the reduced system.
        \item
        For a summer heatwave, we identify a fourth (fast) case of R-tipping
        %
        via analysis of a  codimension-one {\em heteroclinic orbit} connecting the base state for the past limit system to a new saddle  equilibrium that arises in the layer system.
        \end{itemize}
    \item
    Show that a quasithreshold gives rise to {\em critical ranges of rate of change} of the external input rather than isolated critical rates.
    %
    Furthermore, we reveal new types of quasithresholds  through analysis of canard trajectories associated with singular R-tipping edge states.
\end{itemize}

\section{The Nonautonomous Soil-Carbon Model}
\label{sec: rtscm}

The starting point of our analysis is a
discussion of the soil-carbon model introduced
by Luke and Cox~\cite{luke_soil_2011}.
This model describes the time evolution of the soil temperature $T$ and soil carbon concentration $C$ in peat soils
\begin{align}
\label{eq: Intro Model}
    \mu\frac{dT}{dt} &= -\lambda(T-T_{a}(rt)) + A\, C\,R_s(T),\\
\label{eq: Intro Model2}
    \frac{dC}{dt} &= \Pi-C\,R_s(T).
\end{align}
The model incorporates three soil processes and one time-varying external input.
The first process describes equilibration of the soil temperature $T$ with atmospheric temperature $T_a$ according to Newton's Law of Cooling~\cite{maruyama_newtons_2021}, at a rate that depends on the soil-to-atmosphere heat transfer coefficient $\lambda$ and the specific heat capacity of the soil $\mu$.
The second process describes a linear increase in the soil carbon concentration $C$ over time at a
rate $\Pi$, due to carbon generated from decaying plant litter and other processes referred to as  ``Gross Primary Production"~\cite{roy_terrestrial_2001}.\footnote{Note that $dC/dt$ is the ``Net Primary Production" in the model.}
The third and only nonlinear process in the model describes temperature-sensitive microbial activity in the soil in terms of the soil respiration function $R_s(T)$. This
process couples the dynamics of $T$ and $C$, and is discussed in detail in section~\ref{sec:msrhms}.
Atmospheric temperature $T_a(rt)$ is a time-varying external input, which represents weather anomalies or climate variation.
The {\em rate parameter} $r > 0$ quantifies the time scale of climatic variability and is the key {\em input parameter} in the model.
An important aspect of our study is that we use realistic values of the soil parameters based on~\cite{luke_soil_2011}, which are given in Table~\ref{tab: Parameter Table} in the Appendix,
a realistic soil respiration function $R_s(T)$ introduced in
section~\ref{sec:msrhms}, and realistic climate-change scenarios $T_a(rt)$ based on real weather data from Cherskii in Siberia~\cite{noauthor_chersky_nodate}; see the arrow in figure~\ref{fig:RateMap}.

To facilitate analysis, we introduce  a {\em small parameter} $\epsilon = \mu/A \ll 1$
and rewrite {the soil-carbon model~\eqref{eq: Intro Model}--\eqref{eq: Intro Model2} as a fast-slow nonautonomous dynamical system}~\cite{wieczorek_excitability_2011}:
\begin{align} \label{eq: SmallParam}
    \epsilon\frac{dT}{dt} & = f_1(T,C,T_a(rt)) := -\frac{\lambda}{A}(T-T_{a}(rt)) + C\,R_s(T),\\
    \label{eq: SmallParam2}
    \frac{dC}{dt} &=  f_2(T,C)\quad\quad\quad\, := \Pi-C\,R_s(T),
\end{align}
where we define $f_1$ and $f_2$ for convenience.
Note that system~\eqref{eq: Intro Model}--\eqref{eq: Intro Model2} or~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} does not have any equilibria (stationary solutions) owing to the time-varying external input $T_a(rt)$.

\begin{figure}[t]
    \centering
    \includegraphics[]{rsOut.eps}
    \caption{(Black) The modified non-monotone  microbial soil respiration function $R_s^{*}(T)$
    in~\eqref{eq:new_r} with $c=10$ and the die-off temperature $T^{*}=50$ and $70$. (Grey)
    In the limit $T^{*}\to\infty$, we recover the unmodified monotone  microbial soil respiration function
    $R_s^{\dagger}(T)$ in~\eqref{eq:old_r}. See Table~\ref{tab: Parameter Table} for other parameter values.
    }
    \label{fig:Rs Comparison}
\end{figure}

\subsection{ Modified Microbial Respiration}
\label{sec:msrhms}

Our contribution to the model introduced
by Luke and Cox is a modification of the
microbial soil respiration function from~\cite{luke_soil_2011}.
At low to moderate soil temperatures, microbial soil respiration can be described by a
$Q_{10}$ exponential function~\cite{kirschbaum_temperature_1995}:
\begin{equation}
\label{eq:old_r}
R_s^{\dagger}(T)= R_{s}(0)\,e^{\alpha T} \quad \mbox{with} \quad  \alpha = \ln{(Q_{10})}/10,
\end{equation}
where the dimensionless, soil-specific parameter $Q_{10}$ may be estimated from experimental data~\cite{kirschbaum_temperature_1995}.
While the monotone $R_s^{\dagger}(T)$ in~\eqref{eq:old_r} captures
the  trigger mechanism for the R-tipping instability,
it becomes unrealistic at high $T$.
Specifically, $T$ quickly increases to unrealistically high levels due to microbial soil respiration alone~\cite{luke_soil_2011,wieczorek_excitability_2011}.
To address this issue,
we account for an important limitation, that is, soil microbes die above some {\em die-off temperature} $T=T^*$. Specifically, we construct a modified non-monotone microbial soil respiration function
\begin{equation}
\label{eq:new_r}
R^*_{s}(T) = {R}_{s}(0)\,\frac{e^{\alpha b} + e^{-c\alpha b}}{e^{-\alpha (T-b)} + e^{c\alpha (T-b)}}
\quad \mbox{with} \quad
b = T^{*}+\frac{\ln{c}}{\alpha+c\alpha},
\end{equation}
that agrees with the $Q_{10}$ exponential growth~\eqref{eq:old_r} for $T < T^*$,  has a maximum at $T = b \approx T^*$, and decays exponentially for $T > T^*$; see figure~\ref{fig:Rs Comparison}.
Such $R^*_{s}(T)$ prevents soil temperature $T$ from reaching unrealistically high levels. Note that this construction introduces
two additional parameters, namely the die-off temperature $T^*$, and the ratio $c$ of the exponential decay rate for $T>T^*$ and exponential growth rate $\alpha$ for $T<T^*$.
In the remainder of the paper, we use
$$
R_s(T) = R_s^*(T),
$$
for the microbial soil respiration function.

\subsection{The  Autonomous Frozen System and the Moving Equilibrium}
\label{sec:qse}

To gain insight into the tipping mechanisms in the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}, or equivalently system~\eqref{eq: Intro Model}--\eqref{eq: Intro Model2},
we set $r=0$ and consider the resulting autonomous {\em frozen system} for different
but fixed-in-time values of the atmospheric temperature $T_a$.
The frozen system has just one equilibrium
\begin{align}
\label{eq:qse}
e(T_a) &=  (T^{e}(T_a),C^{e}(T_a)) =
\left(T_a+\frac{A\Pi}{\lambda},\frac{\Pi}{R_s(T_a+\frac{A\Pi}{\lambda})}\right),
\end{align}
which is the ``base state" described in the introduction.
Although the position of $e(T_a)$  in the  phase plane $(T,C)$ changes with $T_a$, the equilibrium
remains  linearly stable
and globally attractive\footnote{Meaning that $e(T_a)$ attracts all initial conditions.} within
the realistic range of the soil parameters and $T_a$ used in our study. Thus, we can exclude the possibility of B-tipping from $e(T_a)$~\cite{ashwin2012}.

To discuss R-tipping, we consider the stable equilibrium of the frozen system parameterised by time $t$ for a given input $T_a=T_a(rt)$, which we
denote
\begin{equation}
    \label{eq:me}
e(T_a(rt)),
\end{equation}
and refer to as the
{\em moving  stable equilibrium}~\cite{okeeffe_tipping_2020,Wieczorek2021R}.\footnote{ Note that $e(T_a(rt))$ is also referred to as a ``quasistatic equilibrium" or an ``instantaneous equilibrium"~\cite{ashwin2017}.}.
The moving stable equilibrium is not a solution to the nonautonomous
system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}. However, it can serve as a useful reference.
We follow the approach of~\cite{ashwin2017,okeeffe_tipping_2020,Wieczorek2021R} and relate solutions of the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} to
$e(T_a(rt))$ for different rates $r>0$.
For small $r$, solutions to~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} started near $e(T_a(rt))$
are guaranteed to stay near or track $e(T_a(rt))$~\cite[Th.7.1]{Wieczorek2021R}. However, a nonautonomous R-tipping instability
in the form of a large departure from $e(T_a(rt))$ to an  alternative {\em transient state} for a finite
period of time may appear for larger $r$; see~\cite{luke_soil_2011,wieczorek_excitability_2011,Clarke2021}.
In the next section we give an intuitive explanation of this instability and  propose it as a main cause of
recently observed Zombie fires in peatlands~\cite{BBC2021,Scholten2021,mccarty_arctic_2020,vaughan_arctic_2020} and, possibly, compost fires~\cite{Buggeln2002,NELSON20074612}.

\section{ R-tipping to a Hot Metastable State and Zombie Fires}
\label{sec:RHMSZF}

The  nonlinear process of microbial soil respiration gives rise to a {\em vicious cycle} that
underpins
R-tipping and can be understood in terms of {\em thermal runaway}~\cite{tr_wikipedia}
as follows.
As soil microbes release carbon from the soil into the atmosphere via organic matter decomposition, they generate heat. Typically, the cooler atmosphere allows the soil to achieve a balance between heat and carbon, which is represented by the  moving stable equilibrium $e(T_a(rt))$.
However, this balance may be disturbed
if the atmosphere warms up too fast.
The faster the atmosphere warms up, the faster the soil warms up,
and the quicker the increase in microbial activity, allowing  soil temperature $T$ to reach a level where it can no longer be equilibrated by the cooler atmosphere. This further accelerates  microbial activity, which in turn speeds up the growth of $T$ and the release of soil carbon $C$ into the atmosphere.
Once the soil carbon burns out, the soil temperature drops suddenly and carbon starts accumulating again.
Hence, despite having a linearly stable and globally attractive equilibrium $e(T_a)$ for any fixed value of the atmospheric temperature $T_a$,
the system is very sensitive to how fast the atmospheric temperature $T_a$ increases over time.
In other words, the system may evolve away from the moving stable equilibrium $e(T_a(rt))$ and visit an  alternative transient state for a finite period of time if $T_a(rt)$ rises too fast.

We now demonstrate that the soil respiration function $R^*_s(T)$ introduced in section~\ref{sec:msrmhs}  gives rise to a new alternative transient state at high $T$ that
could explain smouldering ``Zombie fires"~\cite{BBC2021}.
This is evidenced by two remarkable results of the soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with realistic soil parameters, realistic microbial soil respiration function $R_s(T)=R_s^*(T)$, and different climate-change scenarios $T_a(rt)$ based on real weather data from Cherskii in Siberia~\cite{noauthor_chersky_nodate}.

\begin{figure}[t]
    \centering
    \includegraphics[]{CherskiHSA021v2Out.eps}
    \caption{
    (a)
      Realistic seasonal variations of atmospheric temperature $T_a(rt)$ with a summer heatwave
      in year 10, based on observations at Cherskii in Northern Siberia \cite{noauthor_chersky_nodate}; see~\eqref{eq: Figure 2 Ta} for details of $T_a(rt)$.
     (b) Time evolution of (black) moving equilibrium soil temperature $T^{e}(T_a(rt))$ from~\eqref{eq:qse} and (red) the actual soil temperature $T(t)$ obtained by solving~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} in response to the external input $T_a(rt)$ from panel (a).
     Initial soil temperature is $T(0)= -10^\circ$C.
     Initial soil carbon is $C(0) = 120$\,kg\,m$^{-2}$, which corresponds to $\approx3.64$\,m soil depth
    assuming a volumetric carbon density of $33$\,kg\,m$^{-3}$~\cite{khvorostyanov_vulnerability_2008-1, zimov_permafrost_2006}.
    $R_s(T)=R_s^*(T)$ given in~\eqref{eq:new_r} with $T^* = 70^\circ$C.
    Other parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig: Periodic With Hot Summer}
\end{figure}

\begin{itemize}
    \item[(a)]
{\em R-tipping to a hot metastable state.} The more realistic $R_s(T)=R_s^*(T)$ in~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} breaks the vicious cycle (or thermal runaway) giving rise to a {hot metastable state}
    %
    at high soil temperatures $T \approx T^*$. This state is reminiscent of
    ``Zombie fires" observed in tropical and arctic peatlands~\cite{BBC2021,Scholten2021,mccarty_arctic_2020,vaughan_arctic_2020}:
    such fires appear to be extinguished, but smoulder underground throughout the winter and re-emerge the following year~\cite{Scholten2021}.
    Zombie fires are generally believed to happen as a result of surface wildfires~\cite{BBC2021}. Here,
    we suggest that R-tipping to a subsurface hot metastable state
    due to
    atmospheric warming is a main cause of
    %
    Zombie fires.
    Furthermore, we show that such R-tipping can be triggered by realistic climate patterns ranging from summer heatwaves to global warming scenarios.

    Figure~\ref{fig: Periodic With Hot Summer}
    shows the soil temperature change in response to seasonal variations of the atmospheric temperature. Rather amazingly, {\em a summer heatwave} in year 10 breaks the seasonal response pattern and  triggers a sudden transition to a hot metastable state
    that lasts for over a decade\footnote{The duration of the metastable state depends on its temperature $T \approx T^*$.} and releases most of the carbon from the soil into the atmosphere. Following this, the system settles to a lower than initial soil temperature pattern for a refractory period of over a century, during which time both soil carbon and soil temperature slowly increase back to their initial seasonal patterns.

    Figure~\ref{fig: Climate Change} shows the mean soil temperature change in response to  a slow increase in the mean atmospheric temperature of $4^\circ$C over 200 years. We use $4^\circ$C in this example because the global climate change mitigation target is an increase in the global mean temperature of $2 ^\circ$C by 2100 compared to pre-industrial levels (1850–1900)~\cite{IPCCsr15, unfccc2015adoption}.
    %
     Note that the $2 ^\circ$C target is within the range of uncertainty of CMIP5 outputs under two greenhouse gas concentration pathways: the ``very stringent'' RCP2.6 that gives a $1.5^\circ$C increase, and the ``intermediate" RCP4.5 that gives a $2.4^\circ$C increase~\cite{overland_urgency_2019}.
    %
    Due to the well documented phenomenon of ``Arctic amplification''~\cite{overland_urgency_2019,screen_simmonds_2010}, this target corresponds to an increase of roughly $4^\circ$C in the Arctic mean temperature over the same period.
    %
    In this scenario, a sudden transition to a hot metastable state that lasts for over a decade is triggered, and astonishingly occurs after only a modest mean temperature increase of $\approx1.1^\circ$C.
    Following this, the system settles to a lower than initial soil temperature for a refractory period of over a century, during which time both soil carbon and soil temperature slowly increase to their moving equilibrium levels.
    %
   \item[(b)]
   {\em Qualitative and quantitative agreement with PDE models.}
    The conceptual Ordinary Differential Equation (ODE) model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with the more realistic $R_s(T)=R_s^*(T)$
    captures the key nonlinearities of the soil-carbon system. This is evidenced by
    the ODE model
    reproducing both qualitatively and quantitatively the peat soil
    instability predicted by a medium-complexity
    Partial Differential Equation (PDE) model of Siberian permafrost carbon dynamics under climate change~\cite{khvorostyanov_vulnerability_2008-1,khvorostyanov_vulnerability_2008}.
    Specifically, figure~\ref{fig: Khvorostyanov Comparison}
    shows that the ODE model manages to capture the four key features of the PDE dynamics: permafrost thawing, R-tipping to the hot metastable state that lasts for half a century, followed by a sudden soil cooling to slightly above the air temperature.

\end{itemize}

\begin{figure}[t]
    \centering
    \includegraphics[]{CherskiGW2021v2Out.eps}
    \caption{
    (a)
    A realisation
    within the range of uncertainty of
    the ``very stringent" low-emissions RCP2.6 global warming scenario
    for Arctic regions~\cite{pachauri_climate_2015, overland_urgency_2019},
    based on observations at Cherskii in Northern Siberia~\cite{noauthor_chersky_nodate}; see~\eqref{eq: Figure 3 Ta} for details of $T_a(rt)$.
    (b)
    Time evolution of (black) the moving equilibrium soil temperature $T^{e}(T_a(rt))$ from~\eqref{eq:qse} and (red) the actual soil temperature $T(t)$ obtained by solving~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} in response to the external input $T_a(rt)$ from panel (a).
    %
    Initial soil temperature is $T(0)= T_a^-+\frac{A\Pi}{\lambda}\approx-1.95^\circ$C.
    %
    Initial soil carbon $C(0)$ and $R_s(T)$ are the same as in figure~\ref{fig: Periodic With Hot Summer}.
    %
    Other parameter values are given in Table~\ref{tab: Parameter Table}.
    %
    }
    \label{fig: Climate Change}
\end{figure}

The R-tipping instabilities in figures~\ref{fig: Periodic With Hot Summer} and~\ref{fig: Climate Change} arise when time-variation of the atmospheric temperature interacts with different {inherent} timescales of the soil-carbon system.
Our next aim is to understand this interaction and uncover the underlying dynamical mechanisms. Thus, the remainder of the paper is devoted to mathematical analysis of these instabilities in the soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with $R_s(T)=R_s^*(T)$ given in equation~\eqref{eq:new_r} and canonical inputs $T_a(rt)$
{introduced} in sections~\ref{sec:GW} and~\ref{sec:HSA}.

\begin{figure}[]
   \centering
    \includegraphics[]{KHVabOut.eps}
    \caption{(a) The output of a PDE model~\cite{khvorostyanov_vulnerability_2008} parameterised by data from near Cherskii in Siberia to
    simulate a soil column with a depth of 12m and a volumetric carbon density of $33$\,kg\,m$^{-3}$, reproduced from~\cite[Fig.4(a)]{khvorostyanov_vulnerability_2008-1}.
    Shown are time evolution of (black) air temperature at 2 metres above the ground, and (red) soil temperature at a depth of 5m.
    %
    (b)
    Time evolution of (red) the soil temperature {$T(t)$} obtained by solving~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} in response to an increase in (black) atmospheric temperature $T_a(rt)$ given in~\eqref{eq: KhvTa}.
    For comparison with (a),
    we use $\Pi=0.09$ kg m$^{-2}$ y$^{-1}$,  $R_s(T)=0$ for $T\leq0$ as in~\cite{khvorostyanov_vulnerability_2008-1}, and $R_s(T)=R_s^*(T)$  given in~\eqref{eq:new_r} with $T^* = 30^\circ$C for $T>0$.
    %
    Initial soil temperature is $T(0) = -6^\circ$C. Initial soil carbon is $C(0) = 396$\,kg\,m$^{-2}$,
    which corresponds to a soil column with a depth of $12$\,m and a volumetric carbon density of $33$\,kg\,m$^{-3}$ matching~\cite{khvorostyanov_vulnerability_2008-1}.
    Other parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig: Khvorostyanov Comparison}
\end{figure}

\section{The Autonomous Compactified Soil-Carbon Model}
\label{sec: mtc}

To facilitate mathematical analysis of non-obvious dynamical mechanisms that are responsible for the R-tipping instabilities in figures~\ref{fig: Periodic With Hot Summer} and~\ref{fig: Climate Change}, we:
\begin{itemize}
    \item[$\bullet$]
    Consider canonical external inputs $T_a(rt)$ that
    %
    tend exponentially to a constant at infinity and generalise the linear ramp input used in~\cite{wieczorek_excitability_2011}.
    \item[$\bullet$]
    Use the compactification technique of~\cite{Wieczorek_2021,Wieczorek2021R}
    to reformulate the ensuing two-dimensional nonautonomous  soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} as a three-dimensional autonomous compactified system that contains equilibria of the limit systems from infinity.

\end{itemize}

\subsection{Asymptotically Constant External Inputs}
\label{sec: ei}

We choose to work with  external inputs $T_a(rt)$
that  decay exponentially to a constant $T_a^\pm$ as time $t$  tends to $\pm\infty$. To be precise, we follow~\cite[Def.6.1]{Wieczorek2021R} and define
\begin{defn}
\label{defn:expbiconst}
We say $T_a(rt)$ is {\em exponentially bi-asymptotically constant} if,
for all $r > 0$,
\begin{align*}
    \lim_{t\to\pm\infty}T_a(rt)=T_a^\pm\in\mathbb{R}\;\;\;
    \mbox{and}\;\;
    \lim_{t\to\pm\infty}\frac{dT_a(rt)/dt}
{e^{\mp\rho rt}}\in\mathbb{R},
\end{align*}
 for some {\em decay coefficient} $\rho>0$.
\end{defn}
Thus, we can define the autonomous {\em past limit system} with $T_a(rt)=T_a^-$,
\begin{align}
\label{eq:ls-}
    \epsilon\, \frac{dT}{dt} &= f_1(T,C,T_a^-),\\
\label{eq:ls-2}
    \frac{dC}{dt} &= f_2(T,C),
\end{align}
and the autonomous {\em future limit system} with $T_a(rt)=T_a^+$,
\begin{align}
    \label{eq:ls+}
    \epsilon\, \frac{dT}{dt} &= f_1(T,C,T_a^+),\\
        \label{eq:ls+2}
    \frac{dC}{dt} &= f_2(T,C),
\end{align}
which are examples of the frozen system.
We note that,
unlike the nonautonomous soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}, the autonomous past~\eqref{eq:ls-}--\eqref{eq:ls-2} and future~\eqref{eq:ls+}--\eqref{eq:ls+2} limit systems contain the equilibria
\begin{equation}
\label{eq:limeq}
e^- := e(T_a^-)\quad\mbox{and}\quad e^+ :=e(T_a^+),
\end{equation}
respectively, and
$e(rt)\to e^\pm$ as $t\to\pm\infty$ for any $r>0$.
The next step is to include these limit systems and their equilibria in the model.

\subsection{
 Compactification
}
\label{sec:cap}

In this section, we reformulate the two-dimensional nonautonomous  soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with exponentially bi-asymptotically constant input $T_a(rt)$ as a three-dimensional autonomous compactified system that includes the autonomous dynamics and equilibria { $e^-$ and $e^+$} of the limit systems{~\eqref{eq:ls-}--\eqref{eq:ls-2} and~\eqref{eq:ls+}--\eqref{eq:ls+2}} from infinity.

In the standard approach~\cite{kloeden2011}, the two-dimensional nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}
 is augmented with an unbounded dependent variable $u=rt$, which gives an extended autonomous system on $\mathbb{R}^3$ that does not have any equilibria.
Here, we use a different approach.
We follow~\cite[Sec.4]{Wieczorek_2021} and~\cite[Sec.6]{Wieczorek2021R}
and augment the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with {a bounded dependent variable}
\begin{align}
\label{eq:comptrans}
    s = g_{\nu}(rt) = \tanh\left(\frac{\nu}{2}\,rt\right).
\end{align}
We note that $s\in(-1,1)$ for $t\in\mathbb{R}$ and all $r>0$, use the definition of tanh$^{-1}$ to obtain the inverse
\begin{align}
\label{eq:comptransinv}
    rt = g^{-1}_{\nu}(s) = \frac{1}{\nu}\,\ln\frac{1+s}{1-s},
\end{align}
and  differentiate~\eqref{eq:comptrans} with respect to $t$ to derive
the differential equation for $s$:
$$
\frac{ds}{dt} = \frac{\nu r}{2}\,(1-s^2).
$$
Next, we continuously extend the new dependent variable $s$ to include the limits from $t=\pm\infty$, which correspond to $s=\pm 1$.
This gives the {\em autonomous compactified system}
\begin{align}
\label{eq: CompactEqns}
        \epsilon\; \frac{dT}{dt} &= f_1(T,C,T_a^{\nu}(s)),\\
    \label{eq: CompactEqns2}
        \frac{dC}{dt} &= f_2(T,C),\\
     \label{eq: CompactEqns3}
        \frac{1}{r}\;  \frac{ds}{dt} &= \frac{\nu}{2}(1-s^2),
\end{align}
with the continuously extended external input
\begin{equation}
\label{eq: LambdaGeneral}
    T_a^{\nu}(s) =
    \left\{
    \begin{array}{rcl}
    T_a\left( g^{-1}_{\nu}(s) \right) & \text{for} & s\in(-1,1),\\
    T_a^+ & \text{for} & s=1,\\
    T_a^- & \text{for} & s=-1.
    \end{array}
  \right.
\end{equation}
For $T_a(rt)$ with a decay coefficient $\rho$,
we choose the {\em compactification parameter} $\nu\in(0,\rho]$ so the compactified system is continuously differentiable
on the extended phase space $\mathbb{R}^2\times[-1,1]$; see~\cite[Cor.4.1]{Wieczorek_2021} and~\cite[Prop.6.1]{Wieczorek2021R}.
We discuss different choices of $\nu$ for the canonical external inputs in Sections~\ref{sec:GW} and~\ref{sec:HSA}, respectively.

\subsection{

 Equilibria for the Compactified System
}
\label{sec:SECS}

A particular advantage of compactification is that the flow-invariant planes
\begin{align}
    \label{eq:inv_planes}
\mathbb{R}^2\times\{-1\}
\quad\mbox{and}\quad
\mathbb{R}^2\times\{1\},
\end{align}
of the extended phase space $\mathbb{R}^2\times[-1,1]$
contain equilibria $e^{-}$ and $e^{+}$
of the autonomous past{~\eqref{eq:ls-}--\eqref{eq:ls-2}} and future{~\eqref{eq:ls+}--\eqref{eq:ls+2}} limit systems, respectively.
When embedded in the extended phase space,
$e^-$ gains one unstable eigendirection with positive eigenvalue $\nu\, r > 0$ and becomes a hyperbolic saddle
\begin{align}
\label{eq:limeqcomp-}
&\tilde{e}^- = {\left(e^-,-1 \right)},
\end{align}
whereas $e^+$ gains one additional stable eigendirection with negative eigenvalue $-\nu\, r < 0$
and becomes a hyperbolic sink
\begin{align}
\label{eq:limeqcomp+}
& \tilde{e}^+= { \left(e^+,1 \right)},
\end{align}
which is the only attractor for the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}; see ~\cite[Sec.3]{Wieczorek_2021} and~\cite[Prop.6.3]{Wieczorek2021R}.
 Furthermore,  we note that the moving equilibrium $e(T_a(rt))$ with $t\in\mathbb{R}$ corresponds to
$$
\tilde{e}(s) := \left(
e(T_a\left(g^{-1}_{\nu}(s)\right),s
\right)~\mbox{with}~s\in(-1,1),
$$
in the extended phase space of the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}, and $\tilde{e}(s)\to\tilde{e}^\pm$ as $s\to\pm 1$.

\section{Identifying the Model Timescales}
\label{sec:MT}

A key feature of many natural systems is that they evolve on multiple inherent timescales, which can be quantified by
one or more small parameters. The  nonautonomous soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} is one such system, and the left-hand side of~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} shows that it
may evolve on up to three different timescales, depending on the rate parameter $r$.
The {\em slow time} $t$ is the timescale of the {\em slow variable} $C$. The {\em fast time} $\tau=t/\epsilon$
is the timescale of the {\em fast variable} $T$. The {\em third time} $u = rt$ is the timescale of the external input $T_a(rt)$, which is represented by the augmented variable $s$.
The magnitude of the third time $u$ relative to $t$ and $\tau$ depends
on the magnitude of the rate parameter $r$.
Given $0<\epsilon\ll 1$, we consider three cases:
\begin{itemize}
    \item[$\bullet$]
    The {autonomous} {\em frozen system} is represented by $r=0$. This means that the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} becomes a 1-fast 1-slow system with the fast variable $T$, the slow variable $C$, and a fixed-in-time input parameter $T_a$.
    \item[$\bullet$]
    {\em A global warming scenario} is represented by
    $r \lesssim 1$. This means that $u=rt\lesssim t$ is another {\em slow time}, $s$ becomes the second {\em slow variable}, and the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} becomes a {\em 1-fast 2-slow} system with the fast variable $T$ and slow variables $C$ and $s$.
    \item[$\bullet$]
    {\em A summer heatwave} is represented by
    $r \lesssim 1/\epsilon$. This means that $u= \epsilon r \tau \lesssim \tau$ is another {\em fast time}, $s$ becomes the second {\em fast variable}, and the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} becomes a {\em 2-fast 1-slow} system with two fast variables, $T$ and $s$, and one slow variable $C$.
\end{itemize}
In other words, we consider two-timescale limits of a three-timescale problem.
Systems with multiple time scales
can be viewed as singular perturbation problems. In the following sections, we apply concepts and techniques of geometric singular perturbation theory~\cite{fenichel_1979,ckrt_1995,Wechselberger2013,kuehn_2016} first to the  autonomous frozen system and then to the autonomous compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}, which allows us to define and analyse R-tipping in the soil-carbon model.

\section{Defining R-tipping due to Crossing a Quasithreshold}
\label{sec:frozen}

In this section, we combine the compactification from section~\ref{sec: mtc} with concepts and techniques of geometric singular perturbation theory~\cite{fenichel_1979,ckrt_1995,Wechselberger2013,kuehn_2016} to:
\begin{itemize}
    \item
    Define R-tipping due to crossing a quasithreshold in the soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}.
    \item
    Give intuition for when to expect such R-tipping.
\end{itemize}

It is convenient to start the discussion with the frozen system  obtained by setting $r=0$ in system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}, so that $T_a$ becomes a fixed-in-time input parameter. The frozen system is a 1-fast 1-slow singular perturbation problem with a small parameter $\epsilon$. Taking the limit $\epsilon\to0$ in the slow time $t$ gives the singular {\em reduced problem}
\begin{equation}
\begin{split}
\label{eq: frozenreduced}
    \frac{dC}{dt} &= f_2(T,C),
\end{split}
\end{equation}
whose dynamics are restricted to the  one-dimensional critical manifold
\begin{align*}
    S(T_a) = \left\{(T,C)~:~f_1(T,C,T_a) = 0 \right\}\subset\mathbb{R}^2.
\end{align*}
Taking the limit $\epsilon\to0$ in the fast time $\tau=t/\epsilon$ gives the {\em layer problem}
\begin{align}
\begin{split}
\label{eq: frozenLayer}
    \frac{dT}{d\tau} &= f_1(T,C,T_a),
\end{split}
\end{align}
where the slow variable $C$ becomes an additional fixed-in-time parameter.
Note that, for a given $T_a$, the critical manifold $S(T_a)$ consists of all branches of equilibria of the layer problem~\eqref{eq: frozenLayer} parameterised by $C$.
Hence, stability analysis of these equilibria  gives stability of different branches of the critical manifold.
Specifically, $S(T_a)$ has two normally hyperbolic attracting branches, denoted $S_1(T_a)$ and $S_3(T_a)$, which are separated from a normally hyperbolic repelling branch $S_2(T_a)$ by two non-hyperbolic fold points, denoted $F_1(T_a)$ and $F_2(T_a)$.
In other words, the critical manifold can be partitioned as follows
$$
S(T_a) =S_1(T_a) \cup F_1(T_a) \cup S_2(T_a) \cup F_2(T_a) \cup S_3(T_a).
$$
The attracting branch $S_1(T_a)$ contains the equilibrium $e(T_a)$ of the frozen system. The other attracting branch, $S_3(T_a)$, is the hot metastable state
that arises from the non-monotone microbial respiration function~\eqref{eq:new_r}.\footnote{Note that this stable branch does not exist for the monotone respiration function in~\cite{luke_soil_2011}.}
In figure~\ref{fig:frozensystem}(a), we show the dynamics of the reduced and layer problems, and introduce the (blue) {\em special candidate trajectory} $\theta$.

For $0<\epsilon\ll1$,
the normally hyperbolic parts of $S(T_a)$ are guaranteed
by the `first' Fenichel theorem
to perturb smoothly to nearby normally hyperbolic and {\em locally invariant  slow manifolds}~\cite{fenichel_1971,fenichel_1979,ckrt_1995,Wechselberger2013,kuehn_2016}. To be specific,  $S_1(T_a)$ and $S_3(T_a)$ perturb to nearby attracting slow manifolds $S_{1,\epsilon}(T_a)$ and $S_{3,\epsilon}(T_a)$. Similarly, $S_{2}(T_a)$ perturbs to a nearby repelling slow manifold $S_{2,\epsilon}(T_a)$.
Each slow manifold is usually non-unique in the sense that there is a family of such slow manifolds that lie exponentially close in $\epsilon$ to each other; see~\cite[Th.3.1]{Wechselberger2013} and~\cite[Th.3.1.4]{kuehn_2016}. The strategy is to fix a representative
for each of these manifolds and work with these representatives.
Near the non-hyperbolic fold points $F_1(T_a)$ and $F_2(T_a)$,
the hyperbolic branches $S_{1,\epsilon}(T_a)$, $S_{3,\epsilon}(T_a)$ and $S_{2,\epsilon}(T_a)$ typically  split, meaning that they typically  become slow manifolds with  either one or two boundaries.
In particular, $S_{1,\epsilon}(T_a)$ has one inflow boundary near $F_1(T_a)$, $S_{3,\epsilon}(T_a)$ has one outflow boundary near $F_2(T_a)$, while $S_{2,\epsilon}(T_a)$ has one inflow boundary near $F_1(T_a)$ and one outflow boundary near $F_2(T_a)$.
Local invariance means that trajectories can enter or leave a slow manifold only through its boundary.
Furthermore, the splitting of $S_{1,\epsilon}(T_a)$ and $S_{2,\epsilon}(T_a)$ near $F_1(T_a)$
gives rise to
a narrow continuum of special solutions, referred to as {\em canards},
that move along $S_{2,\epsilon}(T_a)$ for some time~\cite{Benoit1983,Dumortier1993,krupaszmolyan2001,SZMOLYAN2001419,Wechselberger2013}{; see figure~\ref{fig:frozensystem}(b).}
In the remainder of the paper
we
follow~\cite{scholarpedia_canard,Benot1981ChasseAC} and~\cite[Sec.3.2.4]{Wechselberger2013}
and
identify three different types of canards.
\begin{defn}
\label{defn:canard}
In the autonomous frozen system and in the autonomous compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}:
\begin{itemize}
    \item[(i)]
    {\em Canards ``without head" } are solutions, or the corresponding trajectories, that move slowly along $S_{2,\epsilon}$ for time $ t= {\cal{O}}(1)$ before moving quickly and directly to $S_{1,\epsilon}$.
    \item[(ii)]
    {\em Canards ``with head"} are solutions, or the corresponding trajectories, that move slowly along $S_{2,\epsilon}$ for time $ t= {\cal{O}}(1)$, then move quickly and directly
    to $S_{3,\epsilon}$, then move slowly along $S_{3,\epsilon}$, before converging to $S_{1,\epsilon}$.
    \item[(iii)]
    A {\em maximal canard}
    is a special solution, or the corresponding trajectory,  that enters $S_{2,\epsilon}$ through its inflow boundary
    %
    and exits $S_{2,\epsilon}$ through its outflow boundary.
\end{itemize}
\end{defn}
Examples of a maximal canard,  denoted $\theta_\epsilon$, a canard ``without head" and a canard ``with head" in the autonomous frozen system are shown in blue, cyan and magenta, respectively, in figure~\ref{fig:frozensystem}(b).
Note that the (blue) maximal canard is a special  example of a canard without head, and a perturbation of the special candidate trajectory $\theta$.
In practice, a maximal canard is computed by choosing a suitable arclength and finding the solution along $S_2(T_a)$ that takes the longest time to trace out this arclength.
The above discussion of canards highlights two main obstacles to defining R-tipping to the hot metastable state. First, the hot metastable state is a transient and thus quantitative phenomenon. In the long term, the
system converges to the same stable state $e^+$ for any rate $r>0$. Second, there is no unique threshold for R-tipping to the hot metastable state~\cite{Wieczorek2021R}. Rather, one speaks of a ``quasithreshold" comprising a family of  canards.
As a consequence, there is no isolated critical rate for such R-tipping.
The fast-slow analysis together with compactification allow us to overcome both of these obstacles.
Our strategy is to relate nonautonomous and compactified system dynamics, and use these relations to define R-tipping for the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} in terms of
slow manifolds and
canards in the autonomous compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}.

\begin{figure}[t]
    \centering
   \includegraphics[]{NewFrozenCombined3Out.eps}
   \caption{(a) The critical manifold $S(T_a)$ with (grey) two normally hyperbolic attracting branches, $S_1(T_a)$ and $S_3(T_a)$,  and (dashed blue) a normally hyperbolic repelling branch $S_2(T_a)$, for the frozen system with $T_a = 0$ and $\epsilon = 0$. Shown also are {(single arrows) the slow dynamics of the reduced problem~\eqref{eq: frozenreduced} along $S(T_a)$ in the slow time $t$, and  (double arrows) the fast dynamics of the layer problem~\eqref{eq: frozenLayer}
   in the fast time $\tau$ for different values of $C$.}
   (b) Five solutions to the frozen system with $T_a = 0$, {$\epsilon \approx 0.064$, the same initial soil temperature $T(0)=0$, and different initial carbon levels: (green) $C(0) = 53$, (cyan) $C(0) = 54.287996178057110079$, (blue) $C(0) \approx 54.28799617805711008$, (magenta) $C(0) = 54.287996178057110088$, and
   (red) $C(0) = 55.5$.
   Included for reference is (grey) $S(T_a)$.
   Other parameter values are given in Table~\ref{tab: Parameter Table}.}
   }
    \label{fig:frozensystem}
\end{figure}

To be precise, we focus on R-tipping from $e^-$.
Using the notation $x=(T,C)\in\mathbb{R}^2$ for the state variable of the soil-carbon system,
we write
\begin{equation}
\label{eq:pullback}
x^{[r]}(t,e^-)\in\mathbb{R}^2,
\end{equation}
to denote the unique
solution to the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} at time $t$, that limits to $e^-$ as $t\to -\infty$ for a given rate parameter $r$.\footnote{This solution can be understood as a {\em local pullback attractor} for~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}~[Th.2.2]\cite{ashwin2017}.}
We also write
$W^{u,[r]}(\tilde{e}^-)$
to denote
the unique rate-dependent one-dimensional unstable invariant manifold of the hyperbolic saddle $\tilde{e}^-$ in the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}.\footnote{For convenience, we leave out the dependence on
the compactification parameter
$\nu$ in the notation for the unstable manifold.}  It follows from~\cite[Prop.6.4(a)]{Wieczorek2021R}
that  $W^{u,[r]}(\tilde{e}^-)$ contains ${x}^{[r]}(t,{e}^-)$ in the sense that
\begin{equation}
\label{eq:Wu}
W^{u,[r]}(\tilde{e}^-)
\supset
\left\{(x,s)~:~ x= x^{[r]}(t,e^-),~s=g_{\nu}(t) \right\}_{t\in\mathbb{R}}.
\end{equation}
This relation allows us to define tracking and R-tipping as follows:
\begin{defn}
\label{defn:track}
Consider the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with exponentially bi-asymptotically constant $T_a(rt)$ and decay coefficient $\rho>0$, and the corresponding autonomous compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}  with the compactification parameter $\nu\in(0,\rho]$. For a fixed $r>0$:
\begin{itemize}
    \item[(i)]
    We say $x^{[r]}(t,e^-)$
    {\em tracks} the moving stable equilibrium $e(T_a(rt))$ in system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} if $W^{u,[r]}(\tilde{e}^-)$
    connects to $\tilde{e}^+$ directly, that is without visiting $S_{2,\epsilon}$ or $S_{3,\epsilon}$,
    in system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}.
    \item[(ii)]
    We say $x^{[r]}(t,e^-)$ {\em R-tips} in system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}, or say system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}
    {\em R-tips}
    from $e^-$, if $W^{u,[r]}(\tilde{e}^-)$ visits
    $S_{3,\epsilon}$
    before
     connecting to
    $\tilde{e}^+$ in system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}.
\end{itemize}
\end{defn}

Next, we recognise that the
nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with $0<\epsilon\ll 1$ does not have unique R-tipping thresholds separating trajectories that track the moving stable equilibrium from those that R-tip,
and does not have isolated critical rates $r$.
More precisely:
\begin{defn}
\label{defn:critrng}
Consider the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with exponentially bi-asymptotically constant $T_a(rt)$ and decay coefficient $\rho>0$, and the corresponding autonomous compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with the compactification parameter $\nu\in(0,\rho]$.
\begin{itemize}
    \item[(i)]
    We define a {\em critical range of $r$} as an interval of $r$ for which $x^{[r]}(t,e^-)$ neither tracks the moving stable equilibrium {$e(T_a(rt))$} nor R-tips.
    \item[(ii)]
    For a fixed $r>0$, we define an {\em R-tipping quasithreshold}
    in system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} as a family of canards ``without head'' including a maximal canard, and in system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} as a family of solutions $x^{[r]}(t)$ or trajectories corresponding to this family of canards.
     \end{itemize}
\end{defn}
\begin{rmk}
There are important differences between {\em R-tipping quasithresholds} defined above and  {\em regular R-tipping thresholds} introduced in~\cite{Wieczorek2021R}.
To understand these differences, we consider the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} and recall that, for typical values of $r$, $W^{u,[r]}(\tilde{e}^-)$ connects to an attractor for the future limit system~\cite[Rmk.5.1]{Wieczorek2021R}.

A {\em regular R-tipping threshold} is a  codimension-one forward-invariant repelling manifold with inflow boundary.\footnote{``Regular" means that it is embedded, orientable and normally hyperbolic.} Such a threshold is contained in the stable manifold of a {\em regular R-tipping edge state} that, in turn, is a compact invariant set with one unstable direction (not an attractor) for the future limit system; see~\cite[Def.5.1]{Wieczorek2021R},~\cite[Def.5.3]{Wieczorek2021R} and~\cite[Prop.6.4(b)]{Wieczorek2021R}.
Crossing a regular threshold typically occurs at special isolated {\em critical rates} $r$ at which $W^{u,[r]}(\tilde{e}^-)$ enters a regular threshold and thus
connects to a regular R-tipping edge state. Note that this behaviour  is unusual for $W^{u,[r]}(\tilde{e}^-)$, which typically connects to an attractor.

In contrast, an {\em R-tipping quasithreshold} in Definition~\ref{defn:canard} is a codimension-zero family of canards ``without head". Thus, crossing a quasithreshold
occurs for {\em critical ranges of} $r$.
Specifically, we expect the following behaviour across a transition from tracking to R-tipping  due to crossing a quasithreshold. As $r$ traverses the critical range, $W^{u,[r]}(\tilde{e}^-)$ contains a canard  ``without head''. For different values of $r$ in the critical range, $W^{u,[r]}(\tilde{e}^-)$ typically contains a different canard  ``without head''. On the boundary of the critical range, $W^{u,[r]}(\tilde{e}^-)$ contains a maximal canard. Just past the critical range, $W^{u,[r]}(\tilde{e}^-)$ contains a canard ``with head'', giving R-tipping as per Definition~\ref{defn:track}(ii).
Since a quasithreshold does not contain any regular edge states, $W^{u,[r]}(\tilde{e}^-)$ always connects to the same attractor $\tilde{e}^+$.

To relate quasithresholds to regular thresholds, we show in Sections~\ref{sec:GW} and~\ref{sec:HSA} that if there is R-tipping due to crossing a quasithreshold in the soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with $0<\epsilon\ll 1$, then  $W^{u,[r]}(\tilde{e}^-)$ connects to a {\em singular R-tipping edge state} in the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with $\epsilon=0$. However, there is no guarantee that a critical range of $r$ for  $0<\epsilon\ll 1$ becomes an isolated critical rate for $\epsilon=0$.
For example, we identify the ``complicated case" in section~\ref{sec:GW3cases}, where the connection to a singular edge state for $\epsilon=0$ persists on a critical range of $r$.
\end{rmk}

To give intuition for when to expect R-tipping due to crossing a quasithreshold in system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}, we extend the concept
of threshold instability, introduced in~\cite[Def.4.5]{Wieczorek2021R} for regular thresholds, to
quasithresholds. Consider phase portraits for the frozen system with $\epsilon=0$ and three different values of $T_a$ in figure~\ref{fig:thr_inst}.
Suppose that $T_a=0^\circ$C and the system is settled at the base state $e(0)$;
see Figure~\ref{fig:thr_inst}(a).
Then, there is a special value
\begin{align}
\label{eq:thr_inst}
T_a = T_a^{inst} \approx -\frac{1}{\alpha}\left(1 + \log{\left(\alpha \frac{A\Pi}{\lambda}\right)} - \alpha \frac{A\Pi}{\lambda}\right),
\end{align}
such that the base state for $T_a=0$ crosses the {\em special candidate trajectory}
for $T_a=T_a^{inst}$, that is $e(0)\in \theta(T_a^{inst})$;
see Figure~\ref{fig:thr_inst}(b) and Appendix~\ref{sec:thr_inst}. If $T_a$ switches discontinuously from $0$ to $T_a^{max}$, then $e(0)$ becomes the initial condition for the frozen system with $T_a=T_a^{max}$. Thus, if $T_a^{max}> T_a^{inst}$, $e(0)$ will find itself on the other side of the { \em special candidate trajectory} $\theta(T_a^{max})$ and the system will undergo R-tipping; see Figure~\ref{fig:thr_inst}(c).  In Sections~\ref{sec:GW} and~\ref{sec:HSA}, we relate condition~\eqref{eq:thr_inst} to R-tipping conditions for $0<\epsilon\ll 1$ and continuously varying $T_a(rt)$.

\begin{figure}[t]
    \centering
   \includegraphics[]{ThresholdCombinedv2Out.eps}
   \caption{
    Phase portraits for the frozen system with $\epsilon=0$ and three different values of $T_a$ show when to expect R-tipping in the nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with $0 < \epsilon\ll 1$ and a time-varying $T_a$. Note the change in the position of the (blue) {\em special candidate trajectory} $\theta(T_a)$ for different values of $T_a$ relative to the fixed initial base state $e(0)$ at $T_a=0$.
   }
    \label{fig:thr_inst}
\end{figure}

In the remainder of the paper, we study the interplay between
the unstable manifold $W^{u,[r]}(\tilde{e}^-)$ and canards in the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} to provide numerical verification and geometric insight into the  different R-tipping mechanisms in figures~\ref{fig: Periodic With Hot Summer} and~\ref{fig: Climate Change}.
Furthermore, we exploit the fact that, in the singular limit $\epsilon=0$, one can work with a reduced problem that approximates the slow dynamics, and with a layer problem that approximates the fast dynamics. These
limit problems
are of lower dimension than the
compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}, making R-tipping analysis more tractable.
We use this approach for two  canonical external inputs $T_a(rt)$ from figure~\ref{fig: ExtInputs}, which are exponentially bi-asymptotically constant and approximate the (a) global warming scenario from figure~\ref{fig: Climate Change} and (b) summer heatwave from figure~\ref{fig: Periodic With Hot Summer}.

\section{R-tipping Mechanisms for Global Warming}
\label{sec:GW}

The R-tipping instability in figure~\ref{fig: Climate Change} arises because time-variation of the  mean atmospheric temperature $T_a$ interacts with the slow timescale of soil carbon $C$. Thus, we work with  the reduced problem for the ensuing 1-fast 2-slow compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} to uncover the underlying dynamical mechanisms.

As a model of the global warming scenario, we consider a slow nonlinear shift from $T_a^-=0^\circ$C to a given $T_a^+>0$ with the rate parameter $r\lesssim1$:
\begin{equation}
\label{eq: Tanh Ta}
    T_a(rt) = \frac{T_a^+}{2}\,(\tanh{(rt)}+1),
\end{equation}
that decays exponentially with the decay coefficient $\rho=2$ as per definition~\ref{defn:expbiconst}; see figure~\ref{fig: ExtInputs}(a).
We fix the compactification parameter $\nu=1$ and consider
system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with
\begin{align}
\label{eq: GWCompact}
T_a^\nu(s) =
    T_a^{1}(s) =
    \frac{T_a^+}{2}\, \frac{(1 + s)^2}{1 + s^2},
\end{align}
which is obtained by applying the inverse compactification transformation~\eqref{eq:comptransinv} to~\eqref{eq: Tanh Ta}. Our
choice of $\nu=1$ ensures that  {the corresponding} system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} is at least $C^1$-smooth on the extended phase space $\mathbb{R}^2\times[-1,1]$
; see section~\ref{sec:cap} and the references therein.

\begin{figure}[t]
   \centering
    \includegraphics[]{ExtInCombinedOut.eps}
    \caption{
    (a) A  nonlinear shift from $T_a^-$ to $T_a^+$ described by Eq.~\eqref{eq: Tanh Ta} approximating the global warming scenario in figure~\ref{fig: Climate Change}. (b) An impulse from $T_a^-=T_a^{min}$ to $T_a^{max}$ and then back to $T_a^+=T_a^{min}$ described by Eq.~\eqref{eq: Sech Ta} approximating the summer heatwave
    in figure~\ref{fig: Periodic With Hot Summer}.
   }
    \label{fig: ExtInputs}
\end{figure}
\begin{figure}[t]
    \centering
    \includegraphics{GWEpsv7Out.eps}
    \caption{
    R-tipping diagram for nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with
    global warming scenario~\eqref{eq: Tanh Ta}, in the plane of the shift amplitude $T_a^+$ and the rate parameter $r$, for $\epsilon\approx0.064$.
    Shown are regions of (green) tracking, (cyan, blue) critical range, and (magenta, red) R-tipping from $\tilde{e}^{-}$. System parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig: GWEps}
\end{figure}

To give an overview of the dynamics near transitions from tracking to R-tipping, we plot an R-tipping diagram in the plane $(T_a^+,r)$
of the input parameters; see figure~\ref{fig: GWEps}. The diagram was obtained by computing $W^{u,[r]}(\tilde{e}^-)$ in system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} for different values of $T_a^+$ and $r$, and  using
Definitions~\ref{defn:canard}--\ref{defn:critrng} to identify different dynamical regions for system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}.
There are two R-tipping regions, and one tracking-tipping transition found for a large enough shift magnitude $T_a^+$. The smaller R-tipping region is a curious R-tipping tongue.
This tongue gives rise to two (cyan) critical ranges of $r$  for a fixed $T_a^+$ (one of them being very narrow),
and suggests a possibility of different
R-tipping cases
for lower and higher magnitudes of the shift.

\begin{figure}[t]
   \centering
    \includegraphics[]{GWLayer135Out.eps}
    \caption{
    (Colour) The unstable invariant manifold $W^{u,[r]}(\tilde{e}^-)$ of the saddle $\tilde{e}^-$  for the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with $\epsilon\approx0.064$, global warming scenario~\eqref{eq: GWCompact} with $T_a^+=5^\circ$C, and three different values of the rate parameter: (green) $r_1=0.02$ , (blue) $r_2\approx0.0454218$, and (red) $r_3=0.1$;
    see the black dots in figure~\ref{fig: GWEps}.
    Included for reference is (grey) the critical manifold $S$ defined in~\eqref{eq:CritManGW}.
    System parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig: GWCompact}
\end{figure}

To gain geometric insight into the R-tipping instability caused by
global warming, we plot in figure~\ref{fig: GWCompact} the unstable manifold $W^{u,[r]}(\tilde{e}^-)$ for a fixed $T_a^+=5^\circ$C and three
different rates $0<r_1 <r_2<r_3$ (see the black dots in figure~\ref{fig: GWEps}), together with the critical manifold
\begin{align}
\label{eq:CritManGW}
    S = \left\{(T,C,s)~:~f_1(T,C,T_a^{1}(s)) = 0\right\}\subset \mathbb{R}^2\times[-1,1],
\end{align}
for reference. Owing to the additional slow $s$-direction, the critical manifold $S$ is now two-dimensional. It consists of two normally hyperbolic attracting sheets, $S_1$ containing $\tilde{e}(s)$, and $S_3$ being the hot metastable state, which are separated from a normally hyperbolic repelling sheet $S_2$ by two non-hyperbolic fold curves $F_1$ and $F_2$.
Figure~\ref{fig: GWCompact} shows that, as $r$ is increased,
tracking of $\tilde{e}(s)$ by (green) $W^{u,[r_1]}(\tilde{e}^-)$ is lost via canard trajectories, including the maximal canard contained in
(blue) $W^{u,[r_2]}(\tilde{e}^-)$ that crosses $F_1$ and moves along $S_2$ for the longest time.
This is followed by R-tipping at higher $r$ as shown by (red) $W^{u,[r_3]}(\tilde{e}^-)$ that visits the hot metastable state $S_{3,\epsilon}$ before connecting to $\tilde{e}^+$.
Thus, it should be possible to explain  R-tipping instabilities caused by
global warming in figure~\ref{fig: GWEps}, including the curious R-tipping tongue, in terms of the slow dynamics of the reduced problem on $S$~\cite{wieczorek_excitability_2011,perryman_adapting_2014,Vanselow_2019}.

\subsection{
The Desingularised Reduced Problem}
\label{sec:GWtm}
The reduced problem for the 1-fast 2-slow system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with global warming~\eqref{eq: GWCompact} is obtained by setting $\epsilon=0$ in~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}, which gives
\begin{align}
\label{eq: GWrp}
        0 &= f_1(T,C,T_a^{1}(s)),\\
        \label{eq: GWrp2}
        \frac{dC}{dt} &= f_2(T,C),\\
        \label{eq: GWrp3}
        \frac{ds}{dt} &=
        \frac{r}{2}\,
        (1-s^2).
\end{align}
This singular system describes the evolution of the slow variables $C$ and $s$ in the slow time $t$ on the two-dimensional critical manifold $S$ defined in~\eqref{eq:CritManGW}. However, the onset of R-tipping involves a sudden jump in the fast variable $T$, without any noticeable change in $C$. Thus, it is convenient to consider the evolution of the fast variable $T$ in the slow time $t$ on $S$, which is obtained  by differentiating the
critical-manifold condition in~\eqref{eq:CritManGW} or~\eqref{eq: GWrp} with respect to $t$. Furthermore,  we use this
condition to project the slow flow within $S$ onto the plane $(T,s)$, i.e. eliminate the dependence on $C$, and reformulate the reduced problem as
\begin{align}
\label{eq: ProjRed}
\frac{d T}{dt} &=
\left.
-
\frac{\partial f_1 /\partial C\cdot f_2 +
\partial f_1/\partial T_a^1 \cdot dT_a^1 /ds \cdot ds/dt}{\partial f_1/\partial T}
\right|_{S},\\
\label{eq: ProjReds}
\frac{ds}{dt} &=
 \frac{r}{2}\,
(1-s^2),
\end{align}
where $|_{S}$ denotes restriction to $S$.

In the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds},
the question of loss of tracking boils down to
whether $W^{u,[r]}(\tilde{e}^-)$
leaves $S_1$ through $F_1$.
To address this question, we note that the denominator in~\eqref{eq: ProjRed} changes sign at a fold
\begin{align}
\label{eq:foldsign}
    \left.\frac{\partial f_1}{\partial T}\right|_{S}
    \left\{
    \begin{array}{ll}
    <0 &\mbox{for}\quad (T,C,s)\in S_1 \cup S_3,\\
    =0 &\mbox{for}\quad (T,C,s)\in F_1 \cup F_2,\\
    >0 &\mbox{for}\quad (T,C,s)\in S_2.
    \end{array}
    \right.
\end{align}
It then follows that, depending on the numerator in~\eqref{eq: ProjRed}, there are two types of fold points:
\begin{itemize}
\item
A {\em jump point}~\cite{SZMOLYAN2001419} is a  point $p$ on a fold such that
\begin{equation}
\label{eq:jp}
\left.\left(\frac{\partial f_1 }{\partial C}\cdot f_2 + \frac{\partial f_1}{\partial T_a^1} \cdot \frac{dT_a^1}{ds} \cdot \frac{ds}{dt}\right)\right|_{p\in F_1 \cup F_2} \neq 0.
\end{equation}
If a solution of~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}
approaches a jump point, the denominator in~\eqref{eq: ProjRed} approaches zero while the numerator remains
finite, so that
$dT/dt$ tends to infinity, and $T(t)$ blows up in $t$.
This means that the solution reaches a jump point in finite time and ceases to exist within $S$.\footnote{For $0<\epsilon\ll 1$, the corresponding solution leaves the  slow manifold and moves along the fast $T$-direction.} Jump points are typically found on open subsets of a fold~\cite{SZMOLYAN2001419}.
\item
A {\em folded singularity}~\cite{SZMOLYAN2001419} is a point $q$ on a fold such that
\begin{equation}
\label{eq:fs}
\left.\left(\frac{\partial f_1 }{\partial C}\cdot f_2 + \frac{\partial f_1}{\partial T_a^\nu} \cdot \frac{dT_a^\nu}{ds} \cdot \frac{ds}{dt}\right)\right|_{q\in F_1\cup F_2} = 0.
\end{equation}
If a solution of~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}
approaches a folded singularity in a special direction where
the numerator in~\eqref{eq: ProjRed} approaches zero
faster than or at the same rate as
the denominator in~\eqref{eq: ProjRed},
then $dT/dt$ remains finite, and this solution either grazes or crosses  the fold at $q$ with a finite speed.
If the crossing is from an attracting to an unstable part of $S$, the solution, or the corresponding trajectory,
is termed a {\em singular canard}~\cite{Benot1981ChasseAC,Benoit1983}.
Folded singularities are typically found as isolated fold points.
\end{itemize}

Thus, in the singular limit $\epsilon=0$ for the slow time $t$, transitions from tracking to R-tipping caused by
global warming can be understood in terms of singular canards through folded singularities~\cite{wieczorek_excitability_2011,Mitry2013,perryman_adapting_2014}.
These folded singularities are examples of the ``singular R-tipping edge states" described in the introduction.
The obstacle to analysis of folded singularities and their singular canards is that the right-hand side of~\eqref{eq: ProjRed} is undefined on $F_1$ and $F_2$.
This obstacle is overcome by a {\em desingularisation}~\cite{SZMOLYAN2001419,dumortier1996canard} in the form of a state-dependent time rescaling:
\begin{align}
\label{eq:tresc}
  dt = \left.-d\hat{t}\; \frac{\partial f_1}{\partial T}\right|_{S},
\end{align}
which gives the {\em desingularised system}
\begin{align}
\label{eq: ProjDesing}
\frac{d T}{d\hat{t}} &=
\left.\left(\frac{\partial f_1}{\partial C}\cdot f_2 +
\frac{\partial f_1}{\partial T_a^1} \cdot \frac{dT_a^1}{ds} \cdot \frac{ds}{dt}\right) \right|_{S}\\
\label{eq: ProjDesing2}
\frac{ds}{d\hat{t}} &= -\left.\left(\frac{\partial f_1}{\partial T} \cdot \frac{ds}{dt}\right)\right|_{S},
\end{align}
defined everywhere on $S$. In detail we obtain:
 \begin{align}
\label{eq:DesingDetail}
\frac{d T}{d\hat{t}} &= R_s^*(T)\left(\Pi-\frac{\lambda}{A}\left(T-T_a^{1}(s)\right)\right) + \frac{r\lambda T_a^+}{2A}\left(\frac{1-s^2}{1+s^2}\right)^2,\\
\label{eq:DesingDetail2}
\frac{ds}{d\hat{t}} &= \frac{r\lambda}{2A}(1-s^2)\left(1+ \alpha\left(T-T_a^{1}(s)\right)\frac{ce^{c\alpha(T-b)}-e^{-\alpha(T-b)}}{e^{c\alpha(T- b)}+e^{-\alpha(T- b)}}\right).
\end{align}
The main advantages of
desingularisation are:
\begin{itemize}
    \item[(i)]
    Regular equilibria for the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} remain regular equilibria
    for the desingularised system~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
    Folded singularities for~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} become (new) regular equilibria for~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}. Hence the classification of folded singularities  into ``folded nodes",
    ``folded foci",
    ``folded saddles",
    and ``folded saddle-nodes",
    based on their classification into different types of equillibria in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}~\cite{SZMOLYAN2001419}.
    \item[(ii)]
    According to~\eqref{eq:foldsign} and~\eqref{eq:tresc}, the new time $\hat{t}$ in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2} flows in the same direction as $t$ on $S_1$ and $S_3$, passes infinitely faster on $F_1$ and $F_2$, and reverses direction on $S_2$. Thus, a phase portrait for~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}
    can be obtained by producing the corresponding phase portrait for~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}, reversing the direction of time (the arrows on trajectories) on $S_{2}$, and relabelling the (new) equilibria on $F_1$ and $F_2$ as folded singularities.
    \item[(iii)]
    It follows from points (i) and (ii) above that a singular canard through a folded singularity in~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}
    can be obtained by smoothly concatenating two trajectories tangent to a stable eigendirection of the corresponding equilibrium
    in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
\end{itemize}

\subsection{Three R-tipping Cases for $\epsilon=0$.
}
\label{sec:GW3cases}
The above relations (i)--(iii) between the desingularised and reduced system dynamics allow us to identify
three different cases of R-tipping  in the  singular reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} through analysis of the regular desingularised system~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
Two of these cases are typical and closely related to the two cases of non-obvious tipping thresholds
described in separate systems in~\cite{perryman_adapting_2014},
and termed as the ``simple case" and the ``complicated case". The third case is special and corresponds to a new {\em degenerate case} that separates the ``simple case" and the ``complicated case". Here, we identify all three cases in an unfolding of a {\em codimension-two global folded saddle-node type-I singularity}, which explains the problem at hand, consolidates the results of~\cite{perryman_adapting_2014},
and is of interest in its own right.

As the rate parameter $r$ is increased in the desingularised system~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}, there is a saddle-node bifurcation of equilibria on $F_1$ at some $r=r_{SN}$, giving rise to a saddle and a stable node. This corresponds to the following local bifurcation scenario in the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}.
When $r=r_{SN}$, a folded saddle-node type-I singularity (FSN-I) appears~\cite{KRUPA20102841,VoWechselberger2015}. As $r$ is increased past $r_{SN}$, FSN-I bifurcates into a folded saddle (FS) and a folded node (FN). For even higher $r$, FN becomes a folded degenerate node (FDN) and subsequently turns into a folded focus (FF).
Furthermore, it turns out that the saddle-node bifurcation in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2} can be a {\em a global  bifurcation} in the sense that it can involve a heteroclinic orbit
from $\tilde{e}^-$.
Thus, we speak of a  {\em global FSN-I}
in the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}.

To gain insight into the  global dynamics in the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds},
we recall that
FSN-I has one {\em folded-saddle-node singular canard} $\tilde{\gamma}^{SN}$, which
crosses
from $S_1$ to $S_2$~\cite{SZMOLYAN2001419}. $\tilde{\gamma}^{SN}$ is obtained by computing the unique stable invariant manifold of the corresponding saddle-node equilibrium in the desingularised system~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
FS has one {\em folded-saddle singular canard} $\tilde{\gamma}^S$, which
crosses
from $S_1$ to $S_2$, and one {\em folded-saddle singular faux canard} $\tilde{\gamma}^S_f$, which
crosses from $S_2$ to $S_1$~\cite{Benoit1983,SZMOLYAN2001419}.
$\tilde{\gamma}^S$ is obtained by computing the unique stable invariant manifold, and $\tilde{\gamma}^S_f$  by computing the unique unstable invariant manifold, of the corresponding saddle in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
In contrast,
FN has a family of singular canards that
cross
from $S_1$ to $S_2$~\cite{SZMOLYAN2001419,Wechselberger2005}. A {\em folded-node strong  singular canard} $\tilde{\gamma}^N_s$ is
obtained by concatenating the two special solutions that approach
the corresponding stable node in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2} along its
 strong eigendirection. In addition,
there is a continuum of {\em  folded-node weak singular canards}, each of which approaches FN along the weak eigendirection of the corresponding stable node in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}. The region of $S_1$ between $\tilde{\gamma}^N_s$ and $F_1$ containing the continuum of weak folded-node singular canards is referred to as the {\em singular funnel}~\cite[Def.3.11]{Wechselberger2013}.\footnote{The singular funnel is shown in yellow in figures~\ref{fig:GWSC},~\ref{fig:GWCC} and~\ref{fig:GWDC}.}
At
FDN, $\tilde{\gamma}^N_s$ turns into $\tilde{\gamma}^{DN}$ and joins the continuum of singular canards that approach FDN along the one (degenerate) eigendirection of the corresponding stable degenerate node in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
FF has no singular canards since every trajectory approaching the corresponding stable focus in~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2} reaches a jump point of $F_1$.

To be precise, we define tracking and R-tipping for slow external inputs and $\epsilon=0$ as follows,
\begin{defn}
In the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}:
\begin{itemize}
    \item [(i)]
    Tracking occurs when $W^{u,[r]}(\tilde{e}^-)$ connects to $\tilde{e}^+$ directly, that is without visiting $F_1$.
    \item[(ii)]
    R-tipping occurs when $W^{u,[r]}(\tilde{e}^-)$ reaches a jump point on $F_1$.
\end{itemize}
\end{defn}
We then note that, when $r<r_{SN}$,
tracking occurs
because the whole of $F_1$ is repelling,
so  $W^{u,[r]}(\tilde{e}^-)$ must remain on $S_1$ and connect directly to $\tilde{e}^+$.
As the rate parameter $r$ is increased past $r_{SN}$,
$W^{u,[r]}(\tilde{e}^-)$ may interact with singular canards to cross $F_1$ and cause loss of tracking.
Analysis of heteroclinic bifurcations in the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2} reveals
three cases of $W^{u,[r]}(\tilde{e}^-)$ crossing $F_1$, or equivalently three cases of loss of tracking:
\begin{itemize}
    \item [(i)] {\em Simple case:}
     $W^{u,[r]}(\tilde{e}^-)$ coalesces with the folded-saddle singular canard $\tilde{\gamma}^S$, and thus crosses from $S_1$ to $S_2$ via FS, at some critical rate $r=r_c>r_{SN}$.
     %
    \item [(ii)] {\em Complicated case:}
     $W^{u,[r]}(\tilde{e}^-)$
     grazes $F_1$ at FSN-I when $r=r_{SN}$,
     coalesces with a  folded-node weak singular canard when $r\in(r_{SN},r_{Ns})$,
     coalesces with the folded-node strong singular canard $\tilde{\gamma}^N_s$ at some $r = r_{Ns}>r_{SN}$,
     and thus crosses from $S_1$ to $S_2$ via FN when $r\in(r_{SN},r_{Ns}]$.
    %
    \item [(iii)] {\em Degenerate case:}
    $W^{u,[r]}(\tilde{e}^-)$ coalesces with the folded-saddle-node singular canard $\tilde{\gamma}^{SN}$,
    and thus crosses from $S_1$ to $S_2$ via FSN-I, at a critical rate $r_c=r_{SN}$.
\end{itemize}
In section~\ref{sec:HSA}, we refer to the three corresponding cases of R-tipping as the {\em slow R-tipping cases}.

We now discuss these three cases in terms of phase portraits for the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} on $S$ by fixing three representative values of the shift magnitude $T_a^+$ and increasing the rate parameter $r$.
To highlight  interactions of $W^{u,[r]}(\tilde{e}^-)$ with different singular canards, we colour $S_1$
in the phase portraits as follows.
Solutions initialised in the {\em green} region remain in this region and converge directly to $\tilde{e}^+$.
Solutions initialised in the {\em dark red}
region reach a jump point of $F_1$ and cease to exist within $S$.
Solutions initialised in the {\em yellow} (singular funnel) region are weak folded-node singular canards.
The green and red regions are separated by the  folded-saddle singular canard $\tilde{\gamma}^{S}$, whereas the red and yellow regions are separated by the  folded-node strong singular canard $\tilde{\gamma}^{N}_s$.
Finally, in the figures, we show projections of $S$ onto the plane $(T,s)$.

\subsubsection{The simple case}

As $r$ is increased past $r_{SN}$ for $T_a^+=1.5$, the appearance of FS and FN via FSN-I gives rise to folded-saddle and folded-node singular canards, and new dynamics. Of particular interest is
an attracting interval of jump points on $F_1$ between FS and FN, and
the corresponding  (dark red) region of solutions
between $\tilde{\gamma}^S$ and $\tilde{\gamma}^N_s$
that reach one of these jump points from $S_1$ and cease to exist within $S$; see figure~\ref{fig:GWSC}~(a)-(b).
In this case, $W^{u,[r]}(\tilde{e}^-)$ does not interact with FSN-I, meaning that the FSN-I is local.
$W^{u,[r]}(\tilde{e}^-)$ does not interact with
the folded-node singular canards either because it is separated from them by the folded-saddle singular canard $\tilde{\gamma}^S$.
Rather, when $r=r_c$, $W^{u,[r]}(\tilde{e}^-)$ and $\tilde{\gamma}^S$ coalesce; see figure~\ref{fig:GWSC}~(c). At this point, tracking is lost since $W^{u,[r]}(\tilde{e}^-)$ crosses $F_1$
from $S_1$ to $S_2$ via FS.
For $r>r_c$, $W^{u,[r]}(\tilde{e}^-)$ reaches a jump point of $F_1$ and ceases to exist within $S$, or R-tips; see figure~\ref{fig:GWSC}~(d).

Thus, in the  simple case for $\epsilon=0$, there is an isolated {\em critical rate} $r=r_c$. This critical rate is the value of
$r$ that gives a {\it codimension-one hetroclinic orbit} connecting $\tilde{e}^-$ to the corresponding saddle on $F_1$ in the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.

\begin{figure}[]
    \centering
    \includegraphics[width=13.5cm]{SimpPortraitShortCombinedOut.eps}
    \caption{
    Phase portraits of the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} with nonlinear shift~\eqref{eq: GWCompact}, obtained using the desingularised system~\eqref{eq:DesingDetail}.
    Shown is the simple case with a codimension-one local FSN-I in (a), and a codimension-one heteroclinic orbit connecting $\tilde{e}^-$ to FS in (c).
    $T_a^+=1.5$ and the rate parameter takes values $(a)\ r=r_{SN}\approx0.107194,\ (b)\ r\approx0.108119,\ (c)\ r = r_c\approx0.111459,$ and  $(d)\ r=0.12$; see the black dots in the left inset of figure~\ref{fig: GWTD}. System parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig:GWSC}
\end{figure}

\begin{figure}[]
    \centering
    \includegraphics[width=13.5cm]{CompPortraitShortCombinedOut.eps}
    \caption{
    Phase portraits of the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}   with nonlinear shift~\eqref{eq: GWCompact}, obtained using the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.  Shown is the complicated case with a codimension-one central heteroclinic FSN-I in (a), and a codimension-one heteroclinic orbit connecting $\tilde{e}^-$ to FN in (c).
    $T_a^+ = 3.5$ and the rate parameter takes values $ (a)\ r = r_{SN} \approx0.0050086,\ (b)\ r\approx0.050669,\ (c)\ r = r_c\approx0.052266,$ and  $(d)\ r=0.06$; see the black dots in the right inset of figure~\ref{fig: GWTD}. System parameter values are given in Table~\ref{tab: Parameter Table}.
    }
\label{fig:GWCC}
\end{figure}

\subsubsection{The complicated case}

As $r$ is increased past $r_{SN}$ for $T_a^+=3.5$, the local bifurcation scenario is the same as in the simple case. However, the global dynamics are different in that $W^{u,[r]}(\tilde{e}^-)$ interacts with various folded-node singular canards.

Similarly to the simple case, appearance of FS and FN via FSN-I gives rise to folded-saddle and folded-node singular canards,
an attracting interval of jump points on $F_1$ between FS and FN, and
the corresponding (dark red) region of solutions between $\tilde{\gamma}^S$ and $\tilde{\gamma}^N_s$ that reach one of these jump points from $S_1$ and cease to exist within $S$; see figure~\ref{fig:GWCC}.
In contrast to the simple case, when $r=r_{SN}$, $W^{u,[r]}(\tilde{e}^-)$
grazes $F_1$ at FSN-I, giving rise to a {\em codimension-one central heteroclinic FSN-I}; see figure~\ref{fig:GWCC}~(a).
In other words, the FSN-I is global.
At this point tracking is lost, but R-tipping cannot occur yet.
For $r$ just above $r_{SN}$, $W^{u,[r]}(\tilde{e}^-)$
passes through the yellow singular funnel and crosses $F_1$ from $S_1$ to $S_2$ via FN; see figure~\ref{fig:GWCC}~(b).  Thus, $W^{u,[r]}(\tilde{e}^-)$ contains a  folded-node weak singular canard.
Furthermore, it attracts all other singular canards within the yellow funnel; see
the relation~\eqref{eq:Wu} and~\cite[Th.2.2]{ashwin2017}.
Then, $W^{u,[r]}(\tilde{e}^-)$ crosses $F_1$ back to $S_1$ via FS and connects to $\tilde{e}^+$; see figure~\ref{fig:GWCC}~(b). Thus, $W^{u,[r]}(\tilde{e}^-)$ also contains the faux folded-saddle singular canard $\tilde{\gamma}^S_f$.
For some
$r=r_{Ns} > r_{SN}$, $W^{u,[r]}(\tilde{e}^-)$ and $\tilde{\gamma}^N_s$ coalesce. At this point, $W^{u,[r]}(\tilde{e}^-)$ crosses $F_1$ from $S_1$ to $S_2$ via FN and never returns to $S_1$; see figure~\ref{fig:GWCC}~(c).
It is only for $r>r_{Ns}$ that $W^{u,[r]}(\tilde{e}^-)$ reaches a jump point of $F_1$ and ceases to exist within $S$, or R-tips; see figure~\ref{fig:GWSC}~(d).

Thus, in the  complicated case for $\epsilon=0$,  there is a {\em critical range of}
$r\in[r_{SN},r_{Ns}]$.
The minimum of the critical range is the value of $r$ that gives a
a {\it codimension-one heteroclinic orbit} connecting $\tilde{e}^-$ to the corresponding saddle-node on $F_1$ along the centre eigendirection in the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
The maximum of the critical range is the value of $r$ that gives a
{\it codimension-one hetroclinic orbit} connecting $\tilde{e}^-$ to the corresponding stable node on $F_1$ along the strong eigendirection in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
For $r$ in the interior of the critical range, there is a {\em codimension-zero heteroclinic orbit} connecting $\tilde{e}^-$ to the corresponding stable node on $F_1$ along the weak eigendirection in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.

\begin{figure}[t]
    \centering
    \includegraphics[width=13.5cm]{DegenShortCombinedOut.eps}
    \caption{
    Phase portraits of the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} with nonlinear shift~\eqref{eq: GWCompact}, obtained using the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}. Shown is the degenerate case with a codimension-two non-central heteroclinic FSN-I in (a).  $T_{a}^+=2.15938$ and the rate parameter takes values $(a)\ r=r_{SN}=r_c \approx 0.0766488,$ and $(b)\ r\approx0.0776488$; see the black dots in the middle inset of figure~\ref{fig: GWTD}. System parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig:GWDC}
\end{figure}

\subsubsection{The degenerate case}

A natural question to ask, is what is the transition between the simple and complicated cases? It turns out that there is a special value $T_a^+=T_{a,c}^+$ at which the transition occurs. As $r$ is increased past $r_{SN}$ for this special value of $T_a^+$, the local bifurcation scenario is the same as in the simple and complicated cases. However, the global dynamics are different in that $W^{u,[r]}(\tilde{e}^-)$ interacts with the folded-saddle-node singular canard.
In contrast to the simple and complicated cases, if $r=r_{SN}$ and $T_a^+=T_{a,c}^+$,  $W^{u,[r]}(\tilde{e}^-)$ and $\tilde{\gamma}^{SN}$ coalesce,
giving rise to a {\em codimension-two non-central heteroclinic FSN-I}; see figure~\ref{fig:GWDC}~(a).
In other words, the FSN-I is global, but different from the complicated case.
At this point, tracking is lost since $W^{u,[r_c]}(\tilde{e}^-)$ crosses $F_1$ from $S_1$ to $S_2$ via FSN-I; see figure~\ref{fig:GWDC}~(a).
For $r>r_{SN}$, $W^{u,[r]}(\tilde{e}^-)$ reaches a jump point of $F_1$ and ceases to exist within $S$, or R-tips; see figure~\ref{fig:GWDC}~(b).

Thus, in the  degenerate case for $\epsilon=0$,  there is an isolated {\em critical pair} $(T_{a,c}^+,r_c)$.
This critical pair is the combination of $r$ and $T_{a}^+$ that give a
{\em codimension-two heteroclinic orbit} connecting $\tilde{e}^-$ to the corresponding saddle-node on $F_1$ along the stable eigendirection in the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.\footnote{Note the difference from typical codimension-one heteroclinic orbits connecting $\tilde{e}^-$ to a saddle-node along the centre eigendirection as in the complicated case.}

\subsection{Bringing it Together: Unfolding of Global FSN-I}
\label{sec:GWbit}

The aim of this section is twofold.
First, we summarise our results from section~\ref{sec:GW3cases} in the singular R-tipping diagram for $\epsilon=0$ in figure~\ref{fig: GWTD}.
Second, we use the singular R-tipping diagram to explain the regular R-tipping diagram for  $\epsilon\approx0.064$ in figure~\ref{fig: GWEps}.

The singular R-tipping diagram is obtained by the unfolding of a {\em codimension-two non-central heteroclinic
FSN-I} in the plane $(T_a^+,r)$ of the input parameters in the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}. This unfolding, in turn,
is obtained by the unfolding of a {\em codimension-two non-central saddle-node heteroclinic bifurcation} in the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.\footnote{This unfolding is reminiscent of the unfolding of a codimension-two non-central saddle-node homoclinic bifurcation in~\cite{chow1990}.}
The  diagram is partitioned into three regions by two curves that are tangent at the special point of  codimension-two non-central hetroclinic FSN-I, denoted $\tilde{e}$-to-FSN-I$_s$ in the middle inset.
Each curve has two branches emanating from the special point.
Both branches of the black curve of FSN-I were obtained by computing the saddle-node bifurcation of equilibria on $F_1$ in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}; we will return to this curve later.
The left branch of the blue curve,  denoted $\tilde{e}$-to-FS, was obtained by computing a codimension-one heteroclinic orbit connecting $\tilde{e}^-$ to the saddle on $F_1$ in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
This connection gives the simple case for
the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}, where
$W^{u,[r]}(\tilde{e}^-) \supset \tilde{\gamma}^S$. Thus, in the simple case, loss of (green) tracking
is the onset of (red) R-tipping.
Note that $\tilde{e}$-to-FS has a vertical asymptote
$T_a^+ = T_a^{inst}\approx 0.423364$, which is given by condition~\eqref{eq:thr_inst}  illustrated in figure~\ref{fig:thr_inst}.
The right branch of the blue curve, denoted $\tilde{e}$-to-FN$_s$,
was obtained by computing a codimension-one heteroclinic orbit connecting $\tilde{e}^-$ to the stable node on $F_1$ along the strong eigendirection in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}. This connection gives the upper boundary of the critical range in the complicated case for~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}, where $W^{u,[r]}(\tilde{e}^-) \supset \tilde{\gamma}^N_s$ and the onset of (red) R-tipping
does not correspond to loss of tracking.
Rather, loss of tracking occurs along the right branch of the black curve, denoted $\tilde{e}$-to-FSN-I$_{c}$. This branch was obtained by computing
a codimension-one heteroclinic orbit connecting
 $\tilde{e}^-$ to the saddle-node on $F_1$ along the centre eigendirection  in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
This connection gives the lower boundary of the critical range in the complicated case for~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}, where tracking is lost because $W^{u,[r]}(\tilde{e}^-)$ grazes $F_1$ at FSN-I; see the red trajectory in figure~\ref{fig:GWCC}(a).\footnote{The left branch of the black curve, denoted FSN-I, does not involve any global bifurcations and does not contribute to tracking-tipping transitions.}
In the (yellow) region between
$\tilde{e}$-to-FSN-I$_{c}$ and $\tilde{e}$-to-FN$_s$,  which is the interior of the critical range,
there is a codimension-zero heteroclinic orbit connecting $\tilde{e}^-$ to the stable node on $F_1$ along the weak eigendirection in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
This connection gives neither tracking nor R-tipping for~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}. Here, $W^{u,[r]}(\tilde{e}^-)$ contains a  folded-node weak singular canard, meaning that it crosses from $S_1$ to $S_2$ via FN, and also contains
the faux folded-saddle singular canard, meaning that it crosses back to $S_1$ via FS; see the red trajectory in figure~\ref{fig:GWCC}(b).
The special point $\tilde{e}$-to-FSN-I$_s$,
where the black and blue curves touch, corresponds to a codimension-two  non-central saddle-node heteroclinic bifurcation that involves a  heteroclinic orbit connecting $\tilde{e}^-$ to the saddle-node on $F_1$ along the stable eigendirection in~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
This connection gives
the degenerate case for~\eqref{eq: ProjRed}--\eqref{eq: ProjReds}, where
$W^{u,[r]}(\tilde{e}^-) \supset \tilde{\gamma}^{SN}$.

\begin{figure}[t]
    \centering
\includegraphics[]{GWTDExtendedv2Out.eps}
    \caption{
     Singular R-tipping diagram for the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} with global warming scenario~\eqref{eq: GWCompact}, obtained using the desingularised system~\eqref{eq:DesingDetail}--\eqref{eq:DesingDetail2}.
     (Green) region of tracking, (red) region of R-tipping, and (yellow) critical range of neither tracking nor R-tipping
     are separated by: (blue) $\tilde{e}^-$-to-FS, $\tilde{e}^-$-to-FSN-I$_s$ and $\tilde{e}^-$-to-FN$_s$, along which $W^{[u],r}(\tilde{e}^-)$ contains the singular canards $\tilde{\gamma}^S$, $\tilde{\gamma}^{SN}$ and $\tilde{\gamma}^{N}_s$, respectively, and (black) $\tilde{e}^-$-to-FSN-I$_c$, along which $W^{[u],r}(\tilde{e}^-)$ contains the weak folded saddle-node singular canard. System parameter values are given in Table~\ref{tab: Parameter Table}.
      }
    \label{fig: GWTD}
\end{figure}

Next, we use the singular R-tipping diagram in figure~\ref{fig: GWTD}, together with rigorous results on
persistence of singular canards as maximal canards for $0 < \epsilon \ll 1$~\cite{SZMOLYAN2001419,KRUPA20102841,krupaszmolyan2001,Wechselberger2005,VoWechselberger2015}, to explain  the regular R-tipping diagram in figure~\ref{fig: GWEps}.

In the simple case, the folded-saddle singular canard $\tilde{\gamma}^S$ perturbs to a family of canards when $0 < \epsilon \ll 1$. The family contains
one maximal canard, namely
the {\em folded-saddle maximal canard} $\tilde{\gamma}^S_\epsilon$, and  associated canards ``without head" and ``with head"~\cite{SZMOLYAN2001419}.
Thus,
these canards explain the simple part of the regular R-tipping diagram at lower values of $T_a^+$, including the (cyan) critical range and the vertical asymptote of the blue curve of maximal canards.
In the context of excitable systems, we use Definition~\ref{defn:critrng}(ii) to relate a family of canards ``without head" associated with a {\em codimension-one normally hyperbolic repelling slow manifold near FS} to a new type of excitability quasithreshold; see also~\cite{wieczorek_excitability_2011,Mitry2013}.

In the complicated case, the perturbed dynamics for $0 < \epsilon \ll 1$ are far less straightforward
due to the presence of
many folded-node maximal canards in addition to the folded-saddle maximal canard.
The number and type of folded-node maximal canards depend on
the distance between FN and FS, and on
the ratio of the weak and strong eigenvalues of FN, denoted $\mu\in(0,1)$. The folded-node strong singular canard $\tilde{\gamma}^N_s$ perturbs to a family of  canards for all $\mu\in(0,1)$. This family contains the {\em folded-node strong maximal canard} $\tilde{\gamma}^N_{s,\epsilon}$ and associated canards ``without head" and ``with head".
Depending on $\mu$, the family of folded-node weak singular canards from the singular funnel may perturb to a single {\em  folded-node weak maximal canard} $\tilde{\gamma}^{N}_{w,\epsilon}$~\cite{Wechselberger2005,Wechselberger2013}.
Additionally, there can be a number of {\em secondary folded-node maximal canards} that lie between $\tilde{\gamma}^N_{s,\epsilon}$ and $\tilde{\gamma}^N_{w,\epsilon}$~\cite{SZMOLYAN2001419,Wechselberger2005,Desroches_2010}.
Furthermore, the interplay between the folded-node and the folded-saddle can give rise to {\em real-faux canards} described in~\cite{perryman_adapting_2014,VoWechselberger2015}, that are a perturbation of the
red trajectory from figure~\ref{fig:GWCC}(b), and to {\em composite canards} identified in~\cite{perryman_adapting_2014}, that follow the folded-saddle maximal canard and then `switch' to and follow one of the folded-node maximal canards.
For example,  we have checked that as $r$ is increased for a fixed $T_a^+ \gtrsim  9.3 ^\circ$C
in figure~\ref{fig: GWEps}, $W^{u,[r]}(\tilde{e}^-)$ coalesces with a secondary folded-node canard at the bottom boundary of the R-tipping tongue, then with a composite canard at the upper boundary of the tongue, and finally with the folded-node strong maximal canard at
the bottom boundary of the main R-tipping region.\footnote{For higher $T_a^{max}$, we found additional R-tipping tongues, not shown in the figures, whose boundaries involve different secondary and composite canards; see also~\cite[Fig.4]{perryman_adapting_2014}.}
Thus, it is folded-node maximal canards and composite canards, together with the associated canards ``without head" and ``with head", that give rise to the complicated tracking-tipping transitions including R-tipping tongue(s), at higher values of $T_a^+$.
In the context of excitable systems, we use Definition~\ref{defn:critrng}(ii) to relate a family of canards ``without head" associated with a {\em codimension-one normally hyperbolic repelling slow manifold near FSN-I} to a new type of excitability quasithreshold. This quasithreshold is expected to have a very complicated shape as indicated by the computations of slow manifolds  in~\cite{desroches2008, perryman2015}.

\section{R-Tipping Mechanism for Summer Heatwaves}
\label{sec:HSA}

The R-tipping instability in figure~\ref{fig: Periodic With Hot Summer} arises because time-variation of the atmospheric temperature $T_a$ interacts with the fast timescale of soil temperature $T$. Thus, we work with the layer problem for the ensuing 2-fast 1-slow compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} to uncover the underlying dynamical mechanisms.

As a model of a summer heatwave,
we consider a fast impulse rising from $T_a^{min}=0^\circ$C to a given $T_a^{max}>0^\circ$C and then dropping back to $0^\circ$C with the rate parameter $r\lesssim1/\epsilon$:
\begin{equation}
\label{eq: Sech Ta}
    T_a(rt) = T_a^{max}\sech{(rt)},
\end{equation}
that decays exponentially with the decay coefficient $\rho=1$ as per definition~\ref{defn:expbiconst}; see figure~\ref{fig: ExtInputs}(b).
We fix the compactification parameter $\nu= 1/2$ and consider
system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with
\begin{align}
\label{eq:HSACompact}
   T_a^\nu (s) =T_a^{\frac{1}{2}}(s) = \frac{2\,T_a^{max}\,(1-s)^{2}}{(1+s)^{4}+(1-s)^{4}},
\end{align}
which is obtained by applying the inverse compactification transformation~\eqref{eq:comptransinv} to~\eqref{eq: Sech Ta}.  Our
choice of $\nu$ ensures that the corresponding system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} is at least $C^1$-smooth on the extended phase space $\mathbb{R}^2\times[-1,1]$; see section~\ref{sec:cap} and the references therein.

\begin{figure}[t]
    \centering
    \includegraphics[]{HSAEpsNewv2OutOut.eps}
    \caption{
    R-tipping diagram for nonautonomous system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2} with summer heatwave~\eqref{eq: Sech Ta}, in the plane of the impulse amplitude $T_a^{max}$ and the rate parameter $r$, for $\epsilon\approx0.064$.
    Shown are regions of (green) tracking, (cyan, blue) critical range, and (magenta, red) R-tipping from $\tilde{e}^{-}$. System parameter values are given in Table~\ref{tab: Parameter Table}.
    \label{fig:HSAEps}
    }
\end{figure}

To give a full overview of
transitions from tracking to R-tipping for impulse inputs, and to make connections to section~\ref{sec:GW}, we plot an R-tipping diagram in the plane $(T_a^{max},r)$
of the input parameters in figure~\ref{fig:HSAEps} for a wide range of the rate parameter $r$. The diagram was obtained by computing $W^{u,[r]}(\tilde{e}^-)$ in system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} for different values of $T_a^{max}$ and $r$, and  using Definitions~\ref{defn:canard}--\ref{defn:critrng} to identify different dynamical regions for system~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}. We note that the shape of the R-tipping region in figure~\ref{fig:HSAEps} is rather different to that obtained in figure~\ref{fig: GWEps}.
There are two R-tipping tongues, one large akin to that in~\cite[Sec.6]{okeeffe_tipping_2020} and one small akin to that in figure~\ref{fig: GWEps}, each enclosing a separate (magenta, red) region of R-tipping.
These two tongues give rise to  up to three critical ranges of $r$  for  a fixed  $T_a^{max}$. Furthermore they are reminiscent of resonance tongues for time-periodic inputs~\cite{marchionne2018} in the sense of a strongly enhanced response to external inputs with
`optimal' timing.
Most importantly, there are multiple tracking-tipping transitions for a large enough shift magnitude $T_a^{max}$.
The lower part of the diagram corresponds to slow impulses that last for decades.
The tracking-tipping transitions found in this part occur during the rise of the impulse $T_a(rt)$, and turn out to  be of the same type as the tracking-tipping transitions described in section~\ref{sec:GW}.
The upper part of the diagram corresponds
to fast impulses and is the focus of this section. Specifically, the rate parameter $r$ in the range of 10 to 15 yr$^{-1}$ gives a pulse duration in the range of 3 to 2 months, in line with realistic summer heatwaves~\cite{perkins-kirkpatrick_lewis_2020}.\footnote{The pulse duration is obtained using the formula
$2\ln(2 + \sqrt{3})/r \approx 2.6/r$ for the full width at half maximum of the hyperbolic secant $\sech(rt)$.} Thus, the single tipping-tracking transition found in the upper part of the diagram quantifies the intensity and duration of  summer heatwaves that trigger R-tipping to the hot metastable state in the soil-carbon model~\eqref{eq: SmallParam}--\eqref{eq: SmallParam2}.

To gain geometric insight into the R-tipping instability caused by
a summer heatwave, we consider~\eqref{eq: Sech Ta} with $r\lesssim 1/\epsilon$, and plot in figure~\ref{fig: HSALayer} the unstable manifold $W^{u,[r]}(\tilde{e}^-)$ for a fixed $T_a^+=15^\circ$C and three
different rates $r_1 >r_2>r_3 >0$ (see the black dots
in figure~\ref{fig:HSAEps}).
Additionally, for reference, we include useful geometric objects. One is the (light gray)
plane of constant $C$ fixed at the equilibrium soil-carbon concentration $C^e(0)$ for
the past and future limit systems. The other is the (dark gray) critical manifold
for~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3},
\begin{align}
\label{eq:critmanHSA}
S = S^- \cup S^+,
\end{align}
which now consists of two disconnected one-dimensional components
\begin{align}
\label{eq:CritManHSA}
S^\mp &=\left\{
(T,C,s)~:~f_1(T,C,0)=0,~s=\mp 1
\right\}.
\end{align}
Specifically,
$
S^- = S^-_1 \cup F_1^- \cup S^-_2 \cup F_2^- \cup S^-_3
$
is the critical manifold of the past limit system~\eqref{eq:ls-}--\eqref{eq:ls-2} with $T_a^-=0^\circ$C, embedded in the compactified phase space where it gains one unstable direction. Similarly,
$
S^+ = S^+_1 \cup F_1^+ \cup S^+_2 \cup F_2^+ \cup S^+_3
$
is the critical manifold of the future limit system~\eqref{eq:ls+}--\eqref{eq:ls+2} with $T_a^+=0^\circ$C, embedded in the compactified phase space where it gains one stable direction. Thus, the attracting branch $S_3^+$ of $S^+$is the hot metastable state. For tracking-tipping transitions, we will be interested in the hyperbolic saddle branch
$S_2^+$ of $S^+$.

Figure~\ref{fig: HSALayer} shows that,
as $r$ is {\em decreased} (meaning the duration of the heatwave is {\em increased}),
tracking by (green) $W^{u,[r_1]}(\tilde{e}^-)$ is lost via canard trajectories, including the maximal canard contained in (blue) $W^{u,[r_2]}(\tilde{e}^-)$ that  approaches
$S_2^+$
and moves along $S_2^+$ for the longest time.
This is followed by R-tipping at {\em lower} $r$ as shown by (red) $W^{u,[r_3]}(\tilde{e}^-)$ that visits the hot metastable state $S^+_{3,\epsilon}$ before connecting to $\tilde{e}^+\in S^+_{1,\epsilon}$.
Since $W^{u,[r]}(\tilde{e}^-)$ remains near $C=C^e(0)$  until $s(t)\approx 1$, it should be possible to explain the R-tipping instability caused by
a summer heatwave
in terms of the fast dynamics of a suitably chosen layer problem for~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}.

\begin{figure}[]
    \centering
     %
     \includegraphics[]{LayerHSAv2Out.eps}
    \caption{
    (Colour) The unstable invariant manifold $W^{u,[r]}_\nu(\tilde{e}^-)$ of the saddle $\tilde{e}^-$ for the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with $\epsilon\approx0.064$, summer heatwave~\eqref{eq:HSACompact} with $T_a^{max}=15^\circ$C, and three different values of the rate parameter: (green) $r= 20$, (blue) $r\approx12.9123$, and (red) $r=10$; see the black dots in figure~\ref{fig:HSAEps}. Included for reference are (dark grey) the two disconnected components $S^-$ and $S^+$ of the critical manifold, (light grey) the layer $\tilde{L}$ defined by $C = C^e(0)$, and (light blue) four new equilibria $\tilde{m}_2^\pm$ and $\tilde{m}_3^\pm$ for the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}.
    System parameter values are given in Table~\ref{tab: Parameter Table}.
    }
    \label{fig: HSALayer}
\end{figure}

\subsection{The Layer Problem}
\label{sec:HSAtlp}

The suitable layer problem for the 2-fast 1-slow system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} with a
summer heatwave~\eqref{eq:HSACompact} is obtained
by  rewriting~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3} in terms of the fast time $\tau = t/\epsilon$, setting $\epsilon= 0$ while
keeping $\epsilon r$ non-zero, and setting $C=C^e(0)$, which gives:
\begin{align}
\label{eq: HSALayer}
        \frac{dT}{d\tau} &= f_1(T,C^e(0),T_a^{\frac{1}{2}}(s)),\\
\label{eq: HSALayer2}
        \frac{ds}{d\tau} &=  \frac{\epsilon r\,}{4}\,(1-s^2).
\end{align}
This system describes the evolution of the fast variables $T$ and $s$ in the fast time $\tau$ on the two-dimensional {\em layer},
$$
    \tilde{L}=\left\{(T,C, s)~:~C = C^e(0) \right\},
$$
shown in figure~\ref{fig: HSALayer}.
Crucially, the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} has six equilibria.
Two of these equilibria, namely the saddle $\tilde{e}^- \in S^-_1$ and the sink $\tilde{e}^+ \in S^+_1$, are the equilibria of the compactified system~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}.
The four new equilibria, that do not exist in~\eqref{eq: CompactEqns}--\eqref{eq: CompactEqns3}, are found at the intersections of $\tilde{L}$ with $S^\pm_{2,3}$.
These are:
the source $\tilde{m}^-_2\in S_2^-$,
the saddle $\tilde{m}^-_3\in S_3^-$,
the saddle $\tilde{m}^+_2\in S_2^+$, and
the sink $\tilde{m}^+_3\in S_3^+$; see the (blue) dots in Figure~\ref{fig: HSALayer}.

\begin{figure}[t]
    \centering
     \includegraphics[]{ThreeTimeScalesCombinedv3Out.eps}
    \caption{
    Phase portraits of the reduced problem for the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} with impulse~\eqref{eq:HSACompact2}. Shown are (gray) different branches of the critical manifold $\tilde{M}$ defined in~\eqref{eq:Mtilde} for (a) $T_a^{max} = 0.2$, (b) $T_a^{max} = T_a^{inst} \approx 0.423364$ and (c) $T_a^{max} = 5$.
    Single arrows indicate the fast time $\tau=t/\epsilon$, and double arrows indicate the intermediate time $u=rt$.
    %
    System parameter values are given in Table~\ref{tab: Parameter Table}.
    %
    We use~\eqref{eq:HSACompact2} instead of~\eqref{eq:HSACompact}  to more easily see the separation of the branches of $\tilde{M}$.
    %
    }
    \label{fig:ReducedLayerProb}
\end{figure}

In the the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}, the question of loss of tracking boils down to whether or not $W^{u,[r]}(\tilde{e}^-)$ connects to the saddle $\tilde{m}_2^+$.
Thus, the saddle $\tilde{m}_2^+$ is another example of the ``singular R-tipping edge state" described in the introduction.
Before we address this question directly in Figure~\ref{fig:HSAfc}, we provide more
insight into R-tipping caused by
a summer heatwave.
We recognise that the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} itself can evolve on two different time scales. Specifically, we rewrite~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} in terms of the intermediate time $u=rt$, take the limit $\epsilon r \to 0$, and consider the critical manifold
\begin{align}
\label{eq:Mtilde}
\tilde{M} =
\left\{(T,s)~:~
0 = f_1(T,C^e(0),T_a^{\frac{1}{2}}(s))
\right\} \subset \tilde{L},
\end{align}
\begin{figure}[]
\centering
    \includegraphics[width=13.5cm]{HSACombinedOut.eps}
\caption{
Phase portraits of the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}
with impulse~\eqref{eq:HSACompact2}.
Shown is the fast case with codimension-one heteroclinic orbit connecting $\tilde{e}^-$ to $\tilde{m}_2^+$ in (b). Included for reference are different branches of (grey) $\tilde{M}$.
$T_a^+=15$, and the rate parameter takes values $(a)\ r=20$, $(b)\ r=r_c\approx15.1674$, and $(c)\ r =10$; see the black dots in the inset of figure~\ref{fig:HSAtd}. System parameter values are given in Table~\ref{tab: Parameter Table}.
We use~\eqref{eq:HSACompact2} instead of~\eqref{eq:HSACompact}  to more easily see the separation of the branches of $\tilde{M}$.
}
\label{fig:HSAfc}
\end{figure}
of the  reduced problem (not shown) for the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}.
It turns out that the number of branches\footnote{In a certain sense, a branch of $\tilde{M}$ corresponds to a moving equilibrium for the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}.} of $\tilde{M}$ depends on $T_a^{max}$
and changes at $T_a^{max} = T_a^{inst}\approx 0.423364$,  which is given by condition~\eqref{eq:thr_inst}  illustrated in figure~\ref{fig:thr_inst}.
If $T_a^{max} < T_a^{inst}$, there are three hyperbolic branches of $\tilde{M}$, two attracting, $\tilde{M}_1$ and $\tilde{M}_3$, and one repelling,
$\tilde{M}_2$. In figure~\ref{fig:ReducedLayerProb}(a), we show only $\tilde{M}_1$ and $\tilde{M}_2$, where the interesting change occurs.\footnote{$\tilde{M}_3$ is shown in figure~\ref{fig:HSAfc}. }
If $T_a^{max} = T_a^{inst}$, $\tilde{M}_1$ and $\tilde{M}_2$ become degenerate in a non-hyperbolic (transcritical) point at $s=0$; see figure~\ref{fig:ReducedLayerProb}(b).
If $T_a^{max} > T_a^{inst}$, $\tilde{M}_1$ and $\tilde{M}_2$ split up into four different branches, two attracting, $\tilde{M}_1^-$ and $\tilde{M}_1^+$,  and two repelling, $\tilde{M}_2^-$ and $\tilde{M}_2^+$. Furthermore, $\tilde{M}_1^\pm$ and $\tilde{M}_2^\pm$ meet at non-hyperbolic fold points $F^\pm$, respectively; see figure~\ref{fig:ReducedLayerProb}(c).
It is clear from figure~\ref{fig:ReducedLayerProb}(a) that tracking in the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} is guaranteed if $T_a^{max} < T_a^{inst}$, owing to the presence of the uninterrupted attracting manifold $\tilde{M}_{1}$ connecting $\tilde{e}^-$ to $\tilde{e}^+$. It is also clear from figure~\ref{fig:ReducedLayerProb}(c) that $T_a^{max} > T_a^{inst}$ is necessary for R-tipping to occur, owing to the presence of a `gap' in $\tilde{M}_{1}$, which allows trajectories to jump from $\tilde{M}_{1}^-$ via $F^-$ to $\tilde{M}_{3}$. Furthermore, it is intuitively clear that R-tipping occurs if $\epsilon r$ is small enough so that the external input $T_a(u)$, or equivalently $s(u)$, changes {\em too slowly} for the system to jump the `gap' in $\tilde{M}_{1}$ between $F^-$ and $F^+$.

To be precise, we define tracking and R-tipping for fast external inputs and $\epsilon=0$ as follows,
\begin{defn}
In the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} with two attractors $\tilde{e}^+$ and $\tilde{m}_3^+$,
\begin{itemize}
    \item [(i)]
    Tracking occurs when $W^{u,[r]}(\tilde{e}^-)$ connects to $\tilde{e}^+$.
    \item[(ii)]
    R-tipping occurs when $W^{u,[r]}(\tilde{e}^-)$ connects to $\tilde{m}_3^+$.
\end{itemize}
\end{defn}
The corresponding dynamics of the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} on $\tilde{L}$
are shown  in a series of phase portraits in Figure~\ref{fig:HSAfc}, where we fix $T_a^{max} = 15 > T_a^{inst}$ and decrease $r$ across the critical rate.
The basin of attraction of $\tilde{e}^+$ is plotted in {\em green}, and the basin of attraction of $\tilde{m}_3^+$ is plotted in {\em dark red}.
(Blue) $W^s(\tilde{m}_2^+)$ is contained in the basin boundary separating these two basins of attraction.
For $r$ sufficiently large, (red) $W^{u,[r]}(\tilde{e}^-)$ lies in the basin of attraction of $\tilde{e}^+$, initially follows $\tilde{M}_1^-$, then jumps the `gap' between $F^-$ and $F^+$, and connects to $\tilde{e}^+$; see figure~\ref{fig:HSAfc}(a).
When $r=r_c$, $W^{u,[r]}(\tilde{e}^-)$ coalesces with $W^s(\tilde{m}_2^+)$ along the basin boundary. At this point, tracking is lost since $W^{u,[r]}(\tilde{e}^-)$ connects to the saddle $\tilde{m}_2^+$; see figure~\ref{fig:HSAfc}(b).
For $r<r_c$, $W^{u,[r]}(\tilde{e}^-)$ lies in the basin of attraction of $\tilde{m}_3^+$, initially follows $\tilde{M}_1^-$, then fails to jump the `gap' between $F^-$ and $F^+$, moves towards $\tilde{M}_3$ and connects to $\tilde{m}_3^+$ along $\tilde{M}_3$. In other words, the system R-tips for $r<r_c$; see figure~\ref{fig:HSAfc}(c).

Thus, in the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}, there is an isolated {\em critical rate} $r=r_c$. This critical rate is the value of $r$ that gives a {\it codimension-one hetroclinic orbit} connecting $\tilde{e}^-$ to the new saddle $\tilde{m}_2^+$ for the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}. This gives just one case of R-tipping, which we refer to as the {\em fast R-tipping case} to distinguish it from the three {\em slow R-tipping cases} identified in section~\ref{sec:GW3cases}.

\subsection{Bringing it Together: Saddle-to-Saddle Heteroclinic Orbit}
\label{sec:HSAbiat}

The aim of this section is twofold. First, we combine our results from section~\ref{sec:HSAtlp} with the approach of section~\ref{sec:GW} to produce the singular R-tipping diagram for $\epsilon = 0$ in figure~\ref{fig:HSAtd}.
Second, we use the singular R-tipping diagram
to explain the regular R-tipping diagram for $\epsilon\approx0.064$ in figure~\ref{fig:HSAEps}.

In the singular R-tipping diagram, the lower-$r$ tracking-tipping transitions were obtained in the same manner as outlined in section~\ref{sec:GW}, that is by computing
different heteroclinic orbits in system~\eqref{eq: ProjDesing}--\eqref{eq: ProjDesing2}.
To avoid repetition,  we skip the derivations and refer to sections~\ref{sec:GW3cases} and~\ref{sec:GWbit} for more details.

The higher-$r$ tipping-tracking transition, which is the focus of this section, is given by the (blue) curve denoted $\tilde{e}^-$-to-$\tilde{m}_2^+$. This curve was obtained by computing a
codimension-one heteroclinic orbit connecting the saddle $\tilde{e}^-$ to the saddle $\tilde{m}_2^+$ in the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}; see figure~\ref{fig:HSAfc}(b).  This connection gives the only case for~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}, namely the {\em fast R-tipping case}, where $\overline{W^{u,[r]}(\tilde{e}^-)} = \overline{W^{s}(\tilde{m}_2^+)}$.\footnote{We write $\overline{A}$ to denote the closure of $A$, that is the smallest closed set containing $A$.}
Thus, in the fast R-tipping case, (red) R-tipping is lost via $W^{s}(\tilde{m}_2^+)$, and this corresponds to the onset of
(green) tracking.
Note that curves $\tilde{e}^-$-to-$\tilde{m}_2^+$ and $\tilde{e}^-$-to-FS have a common vertical asymptote $T_a^{max} = T_a^{inst} \approx 0.423364$, which is given by condition~\eqref{eq:thr_inst}  illustrated in figure~\ref{fig:thr_inst}.

\begin{figure}[t]
    \centering
    \includegraphics[]{HSATDv6Out.eps}
    \caption{
    Singular R-tipping diagram obtained using a combination of (the lower curves) the reduced problem~\eqref{eq: ProjRed}--\eqref{eq: ProjReds} with slow impulse~\eqref{eq:HSACompact} and (the upper curve $\tilde{e}^-$-to-$\tilde{m}^+_2$) the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} with summer heatwave~\eqref{eq:HSACompact}.
    At higher $r$, the (red) region of R-tipping is separated from the (green) region of tracking by (blue) $\tilde{e}^-$-to-$\tilde{m}^+_2$, along which
    $W^{[u],r}(\tilde{e}^-)$ intersects $W^{s}(\tilde{m}^+_2)$ in~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2}.
    %
    For lower $r$, see the caption of figure~\ref{fig: GWTD}.
    %
    System parameter values are given in Table~\ref{tab: Parameter Table}.
    %
    }
    \label{fig:HSAtd}
\end{figure}

Next, in addition to the `first' Fenichel theorem on persistence of normally hyperbolic critical manifolds as slow manifolds for $0 < \epsilon \ll 1$, we recall the `second' Fenichel theorem that guarantees persistence of stable and unstable invariant manifolds of normally hyperbolic saddle critical manifolds when $0 < \epsilon \ll 1$~\cite{fenichel_1971,fenichel_1979,ckrt_1995,Wechselberger2013,kuehn_2016}.
We then use the singular R-tipping diagram in figure~\ref{fig:HSAtd}, together with both Fenichel theorems and rigorous results on
persistence of singular canards as maximal canards for $0 < \epsilon \ll 1$~\cite{SZMOLYAN2001419,KRUPA20102841,krupaszmolyan2001,VoWechselberger2015,Wechselberger2005}, to explain  the regular R-tipping diagram in figure~\ref{fig:HSAEps}.
For the relevant results in the  three {\em slow R-tipping cases} we refer to the discussion in section~\ref{sec:GWbit}, which explains the presence of the small R-tipping tongue and the lower part of the large R-tipping tongue in figure~\ref{fig:HSAEps}.

In the {\em fast R-tipping case}, the heteroclinic $\tilde{e}^-$-to-$\tilde{m}_2^+$ connection for the layer problem~\eqref{eq: HSALayer}--\eqref{eq: HSALayer2} perturbs to an intersection of
$W^{u,[r]}(\tilde{e}^-)$ with the (fast) stable manifold of the normally hyperbolic {\em saddle slow manifold} $S_{2,\epsilon}^+$ in the compactified system~\eqref{eq: CompactEqns}, giving rise to the (blue) maximal canard shown in figure~\ref{fig: HSALayer}.  In other words,
$W^{u,[r]}(\tilde{e}^-)$ approaches the normally hyperbolic saddle slow manifold
$S^+_{2,\epsilon}$ along its stable manifold, thus moves along $S^+_{2,\epsilon}$ for as long $S^+_{2,\epsilon}$ exists, then approaches the stable slow manifold $S^+_{1,\epsilon}$
and eventually connects to $\tilde{e}^+$. In addition, for nearby values of $r$, we
expect associated canards ``without head'' and canards``with head''. These occur when
$W^{u,[r]}(\tilde{e}^-)$ near misses the stable manifold of  $S_{2,\epsilon}^+$ to the left or right, respectively, and diverges away from $S_{2,\epsilon}^+$ before $S_{2,\epsilon}^+$ ceases to exist.
Thus, it is the (blue) maximal canard shown in figure~\ref{fig: HSALayer}, together with the associated canards ``without head",  that give rise to the higher-$r$ tipping-tracking transition
in the regular R-tipping diagram in figure~\ref{fig:HSAEps}.
In the context of excitable systems, we use Definition~\ref{defn:critrng}(ii) to relate a family of canards ``without head" associated with a {\em codimension-one stable manifold of a normally hyperbolic saddle slow manifold} to yet another new type of excitability  quasithreshold.

\section{Conclusion and Outlook}
In this paper, we obtain two kinds of results.
First, we demonstrate that sufficiently fast atmospheric warming can cause  R-tipping to a subsurface hot metastable state in bioactive peat soils, using a conceptual process-based ODE model with realistic soil parameter values and contemporary climate patterns. This gives rise to the hypothesis that such R-tipping is a main cause of ``Zombie fires" observed in peatlands, that disappear from the surface, smoulder underground during the winter, and ``come back to life" in the spring.
Second, we recognise that such R-tipping is a nonautonomous instability that occurs due to crossing an elusive quasithreshold in the phase space of a multiple timescale dynamical system, and thus cannot be explained by traditional autonomous stability theory.  Therefore, to understand the underlying dynamical mechanisms, we provide a mathematical framework that is underpinned by a compactification technique for asymptotically autonomous dynamical systems and concepts from geometric singular perturbation theory, such as slow manifolds, canard trajectories and folded singularities, and use the conceptual soil-carbon model as an illustrative example. This
 explains R-tipping to the hot metastable state in the soil-carbon model and, more generally, identifies generic cases of  R-tipping due to crossing a quasithreshold. Furthermore, it shows that a quasithreshold  gives rise to critical ranges of the rate of change of the external input rather than isolated critical rates, and reveals new types of quasithresholds. These results pose two types of challenges for future research.

First, to strengthen our hypothesis, it would be interesting to build on the excellent agreement with the medium-complexity PDE model in figure~\ref{fig: Khvorostyanov Comparison} and extend the conceptual ODE model in different directions. For example,
    include terms describing peat soil ignition processes at higher temperatures~\cite{yuan_restuccia_rein_2021}. This would allow investigation of existence of additional (quasi)thresholds for the onset of flames.
    %
    Include other physical processes, in addition to soil temperature, that contribute to microbial soil respiration. For example, microbial soil respiration is also strongly moisture dependent~\cite{Sierra_2015,luke_soil_2011,khvorostyanov_vulnerability_2008}, and the inclusion of soil hydrology would extend the applicability of the model; see for example DigiBog\_Hydro~\cite{Digibog_hydro_2022} and~\cite{Waddington2014}.
    Include random fluctuations to account for the fact that climate and weather patterns are inherently noisy. This would give a more accurate description of the tipping being a combination of R-tipping and noise-induced tipping (N-tipping). Understanding how the two tipping mechanisms interact is non-trivial~\cite{Ritchie2017} and an area of ongoing research for  quasithresholds~\cite{Slyman2022}.
    %
    Account for spatio-temporal dynamics by extending the conceptual ODE model to a PDE model. Vertical diffusion has already been studied in the extended Luke and Cox model with interesting results confirming validity of the conceptual ODE model~\cite{Clarke2021}. Further to this, the inclusion of lateral diffusion on a heterogeneous spatial domain would allow investigation of lateral heat and fire patterns, and how they evolve over time, in the ensuing reaction-diffusion equation. This, in turn, would enable comparisons with satellite data of Zombie fires.
    %
    Finally, it might be interesting to couple the conceptual model to spatially extended land-surface models such as JULES~\cite{JULES_2022} or DigiBog~\cite{Digibog_2022}. Currently, such models often neglect heat produced by microbial decomposition which is a key factor for the R-tipping instability to the hot metastable state described in this paper.

Second, there are challenges related to extending the mathematical framework and obtaining rigorous results.
The definitions of R-tipping and a quasithreshold introduced here for the soil-carbon model are a good starting point towards a general theory of R-tipping due to crossing quasithresholds in multiple timescale systems.
Of particular interest to scientists and mathematical modelers are rigorous yet easily testable criteria for the occurrence of such R-tipping, akin to those derived in~\cite{Wieczorek2021R} for R-tipping due to crossing regular thresholds.

\section{Appendix}
\subsection{Details for Figure 1}
The Berkeley Earth dataset~\cite{berkely_earth_data} gives a field of the land temperature anomalies $\Phi(\mathbf{x},m)$ for a given month $m$ compared to the average temperature for that month for the period 1951-1980. We take the mean temperature anomaly field $\bar{\Phi}(\mathbf{x})$ for the 60 months from January 2017 to December 2021 inclusive. A field of the rates of warming, as shown in figure~\ref{fig:RateMap}, is then obtained as follows:
\begin{align}
    r(\mathbf{x}) &= \frac{\bar{\Phi}(\mathbf{x})}{\frac{1}{2}(2021+2017)-\frac{1}{2}(1980+1951)}.
\end{align}

\subsection{External Input and Parameter Details for Figure 2}

In figure~\ref{fig: Periodic With Hot Summer}, the external input is given by
\begin{equation}
\label{eq: Figure 2 Ta}
    T_a(rt) = \frac{T_a^{Jul}-T_a^{Jan}}{2}\left(\sin{(2\pi t)}+1\right)+T_a^{Jan}+T_a^+\sech{\left(r(t-10.25)\right)}.
\end{equation}
where $r=10$,  $T_a^+ = 20$, and $T_a^{Jan} = -33.1^\circ$C , $T_a^{July} = 12.9^\circ$C are the mean January, July temperatures recorded at Cherskii in Northern Siberia \cite{noauthor_chersky_nodate}.

\subsection{External Input and Parameter Details for Figure 3}

In figure~\ref{fig: Climate Change}, the external input is given by
\begin{equation}
\label{eq: Figure 3 Ta}
  T_a(rt) = \frac{T_a^+-T_a^-}{2}\left(\tanh{(r(t-140))}+1\right)+T_a^-,
\end{equation}
where $r=0.025$,  $T_a^+ = -6 $, and $T_a^- = -10.1$ is the mean value of seasonal component of~\eqref{eq: Figure 2 Ta}, which is based on temperatures recorded at Cherskii in Northern Siberia \cite{noauthor_chersky_nodate}.

\subsection{External input and Parameter Details for Figure 4}
In figure ~\ref{fig: Khvorostyanov Comparison}, atmospheric temperature external input is given by
\begin{equation}
\label{eq: KhvTa}
   T_a(rt) = \frac{T_a^+-T_a^-}{2}\left(\tanh{(r(t-2050))}+1\right)+T_a^-.
\end{equation} Parameter values used are $T_a^-=-5^\circ$C,\ $T_a^+=2^\circ$C,$\ \Pi = 0.09,\ r = 0.025$ and $T^{*}=35$, which correspond to values used by \cite{khvorostyanov_vulnerability_2008}.

\subsection{Threshold Instability Condition}
\label{sec:thr_inst}
For $T < T^*$, the modified non-monotone respiration function $R_s^*(T)$ in~\eqref{eq:new_r} is well approximated by the monotone $R_s^{\dagger}(T)$~\eqref{eq:old_r},
which may be used to a derive a closed-form approximation for $F_1$. We recall the fold condition for a fixed value of $T_a$ with $R_s(T)=R_s^{\dagger}(T)$:
\begin{align*}
   \left.\partial f_1/\partial T \right|_{S} & = -\frac{\lambda}{A} \left(1 + (T - T_a\frac{r^\prime_s(T)}{R_s^{\dagger}(T)}\right) = 0.
\end{align*}
This gives an approximation of the $T$-component of the fold $F_1$, which we denote $T^{F_1}$. We use $\approxeq$ to indicate that the expression is exact when $R_s(T)=R_s^{\dagger}(T)$, and near exact (has an exponentially small error) when $R_s(T) = R^*_s(T)$. Applying the critical manifold condition gives
\begin{align*}
   \left(T^{F_1}(T_a),C^{F_1}(T_a)\right) \approxeq \left(\frac{1}{\alpha} + T_a,\frac{\lambda}{\alpha A R_s^{\dagger}(T^{F_1}(T_a))}\right).
\end{align*}
We note from section~\ref{sec:frozen} that there is a special value $T_a^{inst}$ for which $e(0)\in\theta(T_a^{inst})$. This condition is equivalent to $C^{e}(0) = C^{F_1}(T_a^{inst})$, which gives
\begin{align*}
    \frac{\Pi}{R_s^{\dagger}(\frac{A\Pi}{\lambda})} = \frac{\lambda}{\alpha A R_s^{\dagger}(T^{F_1}(T_a^{inst}))}.
\end{align*}
Solving for $T_a^{inst}$ gives the following expression
\begin{align}
    T_a^{inst} \approxeq -\frac{1}{\alpha}\left(1 + \log{\left(\alpha \frac{A\Pi}{\lambda}\right)} - \alpha \frac{A\Pi}{\lambda}\right).
\end{align}

\subsection{Compactified Summer Heatwave}
In figure~\ref{fig:HSAfc},
to more easily see the separation of the branches of $\tilde{M}$
rather than $T_a^{0.5}(s)$~\eqref{eq:HSACompact} we use $T_a^{\nu}(s) = T_a^{0.1}(s)$ which is given by:
\begin{align}
\label{eq:HSACompact2}
   T_a^\nu (s) =T_a^{\frac{1}{10}}(s) = \frac{2\,T_a^{max}\,(1-s)^{10}}{(1+s)^{20}+(1-s)^{20}}.
\end{align}

\begin{table}[t]
    \centering
\begin{tabular}{ | m{1.5cm} || p{6.5cm} | m{4cm} | }
 \hline
 \multicolumn{3}{|c|}{\textbf{Parameters}} \\
 \hline
Symbol & Description & Value\\
 \hline
 $\mu$ & Areal Soil Heat Capacity & $2.5 \times 10^6$\,J\,m$^{-2}$\,$^\circ$C$^-1$\\
 $\lambda$ & Soil-to-Atmosphere Heat Transfer Coefficient &$5.049\times10^6$\,J\,y$^{-1}$\,m$^{-2}$\,$^\circ$C$^{-1}$\\
 $r$ & Rate Parameter & $r\in(0,\infty)$ \,y$^{-1}$\\
 $A$ & Specific Heat of Microbial Respiration & $3.9 \times 10^7$\,J\,kg$^{-1}$\\
 $\Pi$ & Gross Primary Production & $1.055$\,kg\,m$^{-2}$\,y$^{-1}$\\
 $\alpha$ & Respiratory Growth Rate &  $\ln{(2.5)}/10$\\
 $R_s(0)$ & Reference respiration at $0^\circ$C &  $0.01$\\
 $c$ & Decay : Growth Ratio &  $10$\\
 $T^*$ & Respiration Die-Off Temperature &  $30,50 \text{ or }70^\circ$C\\
 \hline
\end{tabular}
\caption{}
\label{tab: Parameter Table}
\end{table}
\newpage

\enlargethispage{20pt}

\bibliographystyle{rs}