#  easyFermi: a graphical interface for performing Fermi-LAT data analyses

Raniere de Menezes1,2 1Universidade de São Paulo, Departamento de Astronomia,
Rua do Matão, 1226, São Paulo, SP 05508-090, Brazil
2Lehrstuhl für Astronomie, Universität Würzburg, Emil-Fischer-Strasse 31,
97074 Würzburg, Germany
Submitted to Astronomy and Computing on 07/March/2022, accepted on
22/June/2022

###### Abstract

Since its launch in 2008, the Fermi Large Area Telescope (LAT) allowed us to
peek into the extremely energetic side of the Universe with unprecedented
sensitivity and resolution. The tools available for analyzing Fermi-LAT data
are the Fermitools and Fermipy, both of which can be scripted in Python and
run via command lines in a terminal or in web-based interactive computing
platforms. In this work, we are providing the community with easyFermi, an
open-source user-friendly graphical interface for performing basic to
intermediate analyses of Fermi-LAT data in the framework of Fermipy. With
easyFermi, the user can quickly measure the $\gamma$-ray flux and photon
index, build spectral energy distributions, light curves, test statistic maps,
test for extended emission and even relocalize the coordinates of $\gamma$-ray
sources. The tutorials for easyFermi are available on YouTube and GitHub,
allowing the user to learn how to use Fermi-LAT data in about 10 min.

###### keywords:

methods: data analysis, methods: miscellaneous, gamma rays: observations

††journal: Astronomy and Computing

##  1 Introduction

The Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope
(Fermi) is a pair-conversion telescope with imaging and spectroscopic
capabilities, a large field-of-view of $\sim 2$ sr, and sensitive to
$\gamma$-rays in the energy range from $\sim 20$ MeV to $\sim 1$
TeV111https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm
[6]. The Fermi-LAT performs an all-sky survey every $\sim 3$ hours and, since
its launch in 2008, has collected more than 13 years of data.

Once the events detected by the Fermi-LAT are classified, the collected
$\gamma$-ray data are released to the public on the Fermi Science Support
Center (FSSC) data server222https://fermi.gsfc.nasa.gov/cgi-
bin/ssc/LAT/LATDataQuery.cgi. The users of these data can rely on
Fermitools333https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ and
Fermipy444https://fermipy.readthedocs.io/en/latest/ [11] to perform their
analyses, both tools relying on several online
tutorials555https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/,
https://fermipy.readthedocs.io/en/latest/quickstart.html. Fermipy, in
particular, is a high-level Python package that facilitates the analysis of
Fermi-LAT data in the framework of the Fermitools.

In this work we present easyFermi, an open-source graphical interface suited
to perform Fermi-LAT data analyses of point-like and extended $\gamma$-ray
sources. easyFermi facilitates the experience and can save a substantial
amount of time for those doing $\gamma$-ray astronomy, especially to
scientists just starting in the field of high-energy astrophysics. Here we
provide an introduction to what is possible to do with easyFermi and how to
use its graphical interface. Further tutorials can be found online on
GitHub666https://github.com/ranieremenezes/easyFermi and
YouTube777https://www.youtube.com/channel/UCeLCfEoWasUKky6CPNN_opQ. This paper
is organized as follows. We describe the installation and setup processes of
easyFermi in §2 and detail what is happening behind the scenes, especially the
dependency of easyFermi on Fermipy, in §3. In §4 we show some of the main data
products of easyFermi and, finally, we conclude in §6.

##  2 Installation and setup

The current (and upcoming) release of easyFermi (i.e. V1.0.7) is available in
the Python Package Index (PyPI) server888https://pypi.org/project/easyFermi/
and GitHub, together with detailed instructions for installation and usage.
The installation of easyFermi (V1.0.7) requires an existing installation of
the Fermitools V2.0.8 and Fermipy V1.0.1 (Python 3 version) and works on Linux
and Mac operational systems. The graphical interface will be maintained such
that its installation will be compatible with the upcoming releases of the
Fermitools and Fermipy.

In summary, once the Fermitools and Fermipy are installed (see the online
tutorial in the easyFermi GitHub webpage), one can use the terminal to open
the fermi environment with conda by typing:

$ conda activate fermi

And then simply type:

$ pip install easyFermi

To test if the installation is properly working, the user can type:

$ python

>> import easyFermi

Figure 1: The main window of easyFermi.

The graphical interface of easyFermi should appear at this point (see Figure
1) and the user can enter the desired configurations.

To exemplify the usage of the software, we will use the blazar PG 1553+113 as
a test target. The first step is to download the spacecraft and photon data
files from the FSSC server999https://fermi.gsfc.nasa.gov/cgi-
bin/ssc/LAT/LATDataQuery.cgi, as well as the Galactic (gll_iem_v07) and
isotropic (iso_P8R3_SOURCE_V3_v1) diffuse emission
models101010https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.
The parameters that we adopt to download the Fermi-LAT data for PG 1553+113
are:

  * 1.

Coordinates: 238.92935, 11.19010 (∘)

  * 2.

Search radius: 10 (∘)

  * 3.

Observation dates: 2008-08-04 15:43:36, 2009-08-04 00:00:00 (Gregorian)

  * 4.

Energy range: 1000, 500000 (MeV)

For a standard point-source analysis of PG 1553+113, the user only needs to
feed the empty boxes of easyFermi with the coordinates (i.e. 238.92935,
11.19010), desired energy range (1000, 500000), the directory where the
spacecraft, photon and background files were downloaded, and the adopted time
interval (from 2008-08-04 15:43:36 to 2009-08-04 00:00:00). For the users with
some experience on Fermi-LAT data analysis, there are a set of advanced
configurations that can also be controlled, but all of them are optional.
Since PG 1553+113 is listed in 4FGL, the configuration step is finished. For
targets not listed in 4FGL (or 3FGL), the user has to change the “Target
cataloged” entry to “No” and insert a nickname for the target in the box
“Target name/tag”.

The user can then simply choose what are the desired outputs by checking the
Light curve, Spectral Energy Distribution (SED), Extension, Re-localize and
test statistic (TS) map boxes. The results will then be saved as “.txt”,
“.npy”, “.fits”, and “.pdf” (or “.png”) files in the selected output
directory. If no output directory is chosen, easyFermi will save all files in
the current working directory. The analysis will start after pressing “Go!”,
however, the user can still make modifications on the fly, like activating or
deactivating one of the boxes in the Science panel. For more complex analysis,
we refer the user to the online tutorials.

##  3 Behind the scenes

The graphical interface and data analysis provided by easyFermi strongly
depend on the Python packages PyQt5 and Fermipy, although requiring minimum
maintenance. Once the user sets the desired configurations and start the
process, easyFermi automatically organizes all these information with PyQt5
and communicates them to Fermipy, thereafter starting a binned likelihood
analysis. This process generates several data subproducts, as a list of
selected events, a counts cube, an exposure map, a livetime cube, and a source
map111111More details on these data subproducts here:
https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html.

Before fitting the model to the data in the region of interest (ROI), we first
optimize the ROI model with the Fermipy function optimize() to ensure that all
parameters are close to their global likelihood maxima, and then look for
uncataloged sources with the function find_sources() (this one can be disabled
by the user). The model is then fitted to the data using the minimizer MINUIT
and the main results for the target are saved in the file
“Target_results.txt”. For the full set of results including all sources in the
ROI, the user can access the files “Results.fits” or “Results.npy”.

###  3.1 Current spectral models

In easyFermi we use the spectral models from the Fermi-LAT third and fourth
source catalogs [3FGL and 4FGL, respectively; 3, 2], namely i) the power-law
model, defined by

$\frac{dN}{dE}=N_{0}\left(\frac{E}{E_{0}}\right)^{\gamma},$

where $N_{0}$ is the normalization (in units of cm-2 s-1 MeV-1), $E$ is the
photon energy, $E_{0}$ is the pivot energy, and $\gamma$ is the photon index;
ii) the log-parabolic model [9] defined by

$\frac{dN}{dE}=N_{0}\left(\frac{E}{E_{0}}\right)^{-\alpha-\beta\log(E/E_{0})},$

where $\alpha$ and $\beta$ are indexes describing the hardness and curvature
of the spectrum; and iii) the power-law with a super exponential cutoff

$\frac{dN}{dE}=N_{0}\left(\frac{E}{E_{0}}\right)^{\gamma}\exp(-aE^{b}),$

where $a$ and $b$ are the exponential factor and index describing the shape of
the spectral cutoff. If the target is listed in one of the Fermi-LAT catalogs,
easyFermi automatically uses the cataloged spectral model in the analysis,
otherwise a new source with a power-law spectrum is added to the ROI model.
The user can always change the spectral model for the target under the
checkbox “Change model”. Depending on the feedback by the users of easyFermi,
we can add more spectral models in the next software releases.

###  3.2 Automatic configuration

For those using the Standard mode of easyFermi, the following set of
configurations apply.

  * 1.

The ROI is defined as an $L\times L$ square with size depending on the adopted
starting energy ($E_{min}$), i.e. $L=15^{\circ},12^{\circ}$ or $10^{\circ}$
for $E_{min}<500$ MeV, 500 MeV $\leq E_{min}<1000$ MeV and $E_{min}\geq 1000$
MeV, respectively.

  * 2.

For the same energy ranges defined in the previous item, we set the maximum
zenith angle to $z_{max}=90^{\circ},100^{\circ}$ or $105^{\circ}$.

  * 3.

The classification and point spread function type for each event are filtered
with evclass $=128$ and evtype $=3$, while the adopted instrument response
function is P8R3_SOURCE_V3 and the dataset is divided in 8 bins per energy
decade.

  * 4.

The radius from the ROI center in which the parameters of the sources are
allowed to vary during the fit is defined as $R_{free}=L/2$, where $L$ is the
size of the ROI.

The experienced user can change all of these configurations by passing a
customized configuration file to easyFermi under the selection of the “custom”
button and by modifying the entries in the “Advanced configurations” box (see
Figure 1).

###  3.3 Recovering the state of the analysis

A very useful characteristic of easyFermi is that it allows the user to
recover the latest state of the analysis, such that the user can quit
easyFermi and continue the exact same analysis some time later. Once the
analysis is done, the state of the graphical interface is automatically saved
in the file “GUI_status.npy” in the output directory, and can easily be
recovered by clicking on “Load GUI state” in the “Menu” button of the toolbar.

##  4 Main data products

The main results from easyFermi are the measurements of the flux and spectral
shape for all sources in the ROI, which are saved in the files
“Target_results.txt”, “Results.npy”, and “Results.fits”. Furthermore, the user
can relocalize the target (the ROI is then updated with the new location of
the target), compute a light curve, a $\gamma$-ray spectrum, look for extended
emission and compute a TS map.

In Figures 2, 3, and 4 we show the $\gamma$-ray spectrum and light curves for
PG 1553+113 built with the dataset described in §2. Upper-limits are displayed
whenever an energy or time bin has TS $\leq 9$. For the extension, since PG
1553+113 is a point-like source, the peak in the delta log-likelihood shown in
Figure 5 is centered at zero.

All of these results are saved in files named “SOURCENAME_task.fits” and
“SOURCENAME_task.npy”, where “task” here can be “loc”, “lightcurve”, “SED”, or
“extension”. These files can easily be accessed with numpy [8], astropy [4, 5]
and TopCat [10]. Furthermore, easyFermi automatically plots the results as
.pdf or .png figures labeled as, e.g. “Quickplot_task.pdf”.

Figure 2: The $\gamma$-ray spectrum for PG 1553+113 fitted with a log-
parabolic model. Figure 3: Photon flux light curve for PG 1553+113. Figure 4:
Energy flux light curve for PG 1553+113. Figure 5: The extension of PG
1553+113 is compatible with zero, meaning that this source is point-like.

###  4.1 Goodness of fit

The detection, flux determination and spectral modeling of Fermi-LAT sources
with easyFermi is accomplished by means of a binned likelihood optimization
technique [1], using MINUIT as the minimizer. The multi-dimensional
minimization of parameters, however, is not an easy task and is susceptible to
fail. To be sure that the fit properly converged, the user can follow the
information in the Log box in the graphical interface (see Figure 1), where
the following information is displayed depending on the success of the fit:

  * 1.

Fit quality: 3. Excellent fit. Full accurate covariance matrix.

  * 2.

Fit quality: 2. Reasonable fit. Full matrix, but forced positive-definite
(i.e. not accurate).

  * 3.

Fit quality: 1. Poor fit. Diagonal approximation only, not accurate.

  * 4.

Fit quality: 0. Bad fit. Error matrix not calculated.

For any fit quality value other than 3, we recommend the user to rerun the
analysis changing the configurations in the panel Free source radius, in the
“Advanced configurations” box. In this panel, the user can select the radius
(from the ROI center) defining the circular region in which the parameters of
the $\gamma$-ray sources are left free to vary, can fix the spectral shape of
the target or of all sources in the ROI, and can freeze the Galactic and
Isotropic diffuse emission models. For a step-by-step guide on the
quality/goodness of fit, the user is invited to follow the tutorial on
YouTube121212https://www.youtube.com/channel/UCeLCfEoWasUKky6CPNN_opQ.

##  5 Validation, performance, and maintenance

The results obtained with easyFermi are exactly the same as those obtained
with Fermipy, since easyFermi is actually running Fermipy in the background.
We tested this by performing several analyses (e.g., the light curves, SEDs
and TS maps for ROIs centered on PG1551+113, 3C 273, 3C279, Omega Centauri,
and the extended lobes of Centaurus A) with easyFermi and comparing the
results with pure Fermipy analyses. In terms of performance, an analysis done
with easyFermi takes basically as much CPU time and RAM as an analysis done
directly with Fermipy on the terminal, however, easyFermi is slightly faster
and uses less RAM than an analysis performed with Fermipy in web-based
interactive computing platforms. To test if easyFermi is properly working, we
recommend the user to follow the tutorial on
GitHub/YouTube131313https://github.com/ranieremenezes/easyFermi. Furthermore,
the user can check if both Fermipy and easyFermi are working smoothly by
simply typing

$ python

>> from fermipy.gtanalysis import GTAnalysis

>> import easyFermi

If both modules are successfully imported, the installation is fine.

The maintenance of easyFermi will be done such that the users can ask for the
repair of possible issues via GitHub and in a way that we guarantee its
compatibility with the upcoming releases of the Fermitools and Fermipy.

##  6 Conclusions

The easyFermi graphical interface is a user-friendly tool that allows
astronomers from all niches to use Fermi-LAT $\gamma$-ray data. This tool is
especially indicated for those scientists just starting in the field of high-
energy astrophysics, and can be used for several goals like building a light
curve, computing a $\gamma$-ray spectrum, or looking for extended $\gamma$-ray
emission. Furthermore, easyFermi is meant to be simple. For more complex types
of analyses, we refer the user to Fermipy and Fermitools.

The tutorials and source code for easyFermi can be found
online141414https://github.com/ranieremenezes/easyFermi and allow the user to
learn how to use Fermi-LAT data in just a few minutes. Bug reports and
proposals for new functionality should be made through the GitHub issue
tracker. Although minimum maintenance is required, we aim to keep easyFermi
updated and in synergy with Fermipy, hence new versions of the tool can be
released in the future. We also plan to make a similar graphical interface for
the Cherenkov Telescope Array [7] and possibly for other Imaging Atmospheric
Cherenkov Telescopes in the near future.

## Acknowledgements

I would like to thank the anonymous referees for their suggestions and
comments, as well as Alessandra Azzollini, Clodomir Vianna, Douglas Carlos,
Fabio Cafardo, Kaori Nakashima, Lucas Costa Campos, Lucas Siconato, Raí
Menezes, Rodrigo Lang, and Romana Grossova for installing, testing, and
helping me with the development of easyFermi. Part of this project was also
supported by the European Research Council for the ERC Starting grant
MessMapp, under contract no. 949555.

## Data Availability

All data used to exemplify the usage of easyFermi in this work can be found
online in the Fermi-LAT data server151515https://fermi.gsfc.nasa.gov/cgi-
bin/ssc/LAT/LATDataQuery.cgi.

## References

  * [1] Abdo, Aous A., Ackermann, M., Ajello, M. et al. 2009, ApJS, 183, 46.
  * [2] Abdollahi, S., Acero, F., Ackermann, M. et al. 2020, Fermi large area telescope fourth source catalog, ApJS, 247, 33.
  * [3] Acero, F., Ackermann, M., Ajello, M. et al. 2015, Fermi large area telescope third source catalog, ApJS, 218 (2), 23.
  * [4] Astropy Collaboration, Robitaille, Thomas P., Tollerud, Erik J., Greenfield, P. et al. 2013, Astropy: A community Python package for astronomy, Astronomy & Astrophysics, 558, A33.
  * [5] Astropy Collaboration, Price-Whelan, A.M., Sipőcz, B.M., Günther, H. M. et al. 2018, The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package, The Astronomical Journal, 156 (3), 123.
  * [6] Atwood, W.B., Abdo, Aous A., Ackermann, M. et al. 2009, The large area telescope on the Fermi gamma-ray space telescope mission, ApJ, 697 (2), 1071.
  * [7] Bernlöhr, K., Barnacka, A., Becherini, Y. et al. 2013, Monte Carlo design studies for the Cherenkov telescope array, Astroparticle Physics, 43, 171.
  * [8] Harris, Charles R., Millman, K. Jarrod, van der Walt, Stéfan J. et al. 2020, Array programming with NumPy, Nature, 585, 7825.
  * [9] Massaro, E., Perri, M., Giommi, P., Nesci, R. 2004, Log-parabolic spectra and particle acceleration in the BL Lac object Mkn 421: Spectral analysis of the complete BeppoSAX wide band X-ray data set, Astronomy & Astrophysics, 413 (2), 489.
  * [10] Taylor, M.B. 2005, TOPCAT & STIL: Starlink Table/VOTable Processing Software, Astronomical Data Analysis Software and Systems XIV, Astronomical Society of the Pacific Conference Series, 347, 29.
  * [11] Wood, M., Caputo, R., Charles, E. et al. 2017, Fermipy: An open-source Python package for analysis of Fermi-LAT Data, 35th International Cosmic Ray Conference (ICRC2017), 301, 824.


\mathcal{V}\end{bmatrix}+\operatorname{\varepsilon}^{\frac{4}{3}}d\begin{bmatrix}{\mathcal{U}}_{xx}\\\
{\mathcal{V}}_{xx}\end{bmatrix}+\begin{bmatrix}\mu_{R}u_{\rm
Q}-(\mu_{I}+\omega_{0})v_{\rm Q}\\\ \mu_{R}v_{\rm
Q}+(\mu_{I}+\omega_{0})u_{\rm
Q}\end{bmatrix}+\operatorname{\varepsilon}\begin{bmatrix}{\rm Re}(d(A_{\rm
QSS})_{xx})\\\ {\rm Im}(d(A_{\rm QSS})_{xx})\end{bmatrix},$

where $\mu=\mu_{R}+i\mu_{I}$ and

$M=\begin{bmatrix}\operatorname{\varepsilon}^{\frac{1}{3}}\mu_{R}-3u_{\rm
Q}^{2}-v_{\rm Q}^{2}+2\alpha u_{\rm Q}v_{\rm
Q}&-\operatorname{\varepsilon}^{\frac{1}{3}}(\mu_{I}+\omega_{0})+\alpha u_{\rm
Q}^{2}+3\alpha v_{\rm Q}^{2}-2u_{\rm Q}v_{\rm Q}\\\
\operatorname{\varepsilon}^{\frac{1}{3}}(\mu_{I}+\omega_{0})-3\alpha u_{\rm
Q}^{2}-\alpha v_{\rm Q}^{2}-2u_{\rm Q}v_{\rm
Q}&\operatorname{\varepsilon}^{\frac{1}{3}}\mu_{R}-u_{\rm Q}^{2}-3v_{\rm
Q}^{2}-2\alpha u_{\rm Q}v_{\rm Q}\end{bmatrix}.$

The trace of $M$ is

${\rm tr}(M)=2\operatorname{\varepsilon}^{\frac{1}{3}}\mu_{R}-4(u_{\rm
Q}^{2}+v_{\rm Q}^{2}),$ (9.4)

and

${\rm
det}(M)=\operatorname{\varepsilon}^{\frac{2}{3}}(\mu_{R}^{2}+(\mu_{I}+\omega_{0})^{2})-4\mu_{R}\operatorname{\varepsilon}^{\frac{1}{3}}(u_{\rm
Q}^{2}+v_{\rm Q}^{2})(1+\mu_{I}+\omega_{0})+3(1+\alpha^{2})(u_{\rm
Q}^{2}+v_{\rm Q}^{2})^{2},$

so that ${\rm det}(M)>0$ for all $\mu_{R}<0$, as well as for $\mu_{R}>0$, at
least until $\mathcal{O}(\operatorname{\varepsilon}^{-\frac{1}{3}})$.

Figure 18:  With the large-amplitude source term
$I_{\tilde{a}}(x)=\tilde{a}e^{-\frac{x^{2}}{4\sigma}}$, $\tilde{a}=0.2$, and
$\sigma=\frac{1}{4}$, the instantaneous Hopf bifurcation curve is $\mu_{\rm
Hopf}(x)=(2.8655...)e^{-\frac{2x^{2}}{3}}$ (blue curve). For each $x$, the
solution stays near the repelling QSS for some time (with ${\rm Re}(A)$ being
orange-yellow for $|x|\lesssim 2.24$, and green for $|x|\gtrsim 2.24$). Here,
the black dots mark the points $x=\pm 2.24$, where $\mu_{\rm
Hopf}=\sqrt{\operatorname{\varepsilon}}$. For $|x|\lesssim 2.24$, the duration
of the delayed onset of the oscillations is approximately $\omega_{0}$ time
units beyond the spatially-dependent curve $\mu_{\rm Hopf}(x)$. Then, for
$|x|\gtrsim 2.24$, where $\mu_{\rm
Hopf}\leq\sqrt{\operatorname{\varepsilon}}$, the oscillations set in just
before the spatio-temporal buffer curve (black curve), as expected since the
QSS has small-amplitude here and the curve is again given to leading order by
(3.13) for these $x$. Simulation performed on $[-50,50]$. For $|x|>5$ (not
shown), the oscillations also commence just before the space-time buffer
curve. The parameters are
$\operatorname{\varepsilon}=0.01,\omega_{0}=\frac{1}{2},d_{R}=1,d_{I}=0.5,\alpha=0.6$.

Therefore, for the solutions of (9.1), the Hopf bifurcation occurs to leading
order at the $x$-dependent time given by

$\mu_{\rm Hopf}(x)=\frac{2}{\operatorname{\varepsilon}^{\frac{1}{3}}}(u_{\rm
Q}^{2}+v_{\rm
Q}^{2})=\operatorname{\varepsilon}^{-\frac{1}{3}}\frac{2(I_{\tilde{a}}(x))^{\frac{2}{3}}}{(1+\alpha^{2})^{\frac{1}{3}}}.$
(9.5)

It is illustrated in Figures 18, 19, and 20. Also, we have checked (at several
points $x$) that the numerically observed duration of DHB in the full
nonlinear PDE scales as $\operatorname{\varepsilon}^{-\frac{1}{3}}$ (data not
shown).

In Figure 18, we compare the results of numerical simulations of (9.1) with a
Gaussian source term, $I_{\tilde{a}}(x)=\tilde{a}e^{-\frac{x^{2}}{4\sigma}}$,
${\tilde{a}}=0.2$, and $\sigma=\frac{1}{4}$, to the analytical results. The
Hopf bifurcation curve is $\mu_{\rm
Hopf}(x)=(2.8655...)e^{-\frac{2x^{2}}{3}}$, following (9.5). The narrow peak
of the repelling QSS manifests as the orange and yellow band about $x=0$. At
all points, the full solution stays near the repelling QSS at least until
$\mu$ reaches $\omega_{0}$. After that, the exit time from the neighbourhood
of the repelling QSS is $x$-dependent. In particular, for $|x|\geq
x_{*}=2.24...$, $\mu_{\rm Hopf}(x)\leq 0.1=\sqrt{\operatorname{\varepsilon}}$.
Hence, the use of the space-time buffer curve obtained from the linearisation
about $A=0$ (for which $\mu_{\rm Hopf}=0$ to leading order) is consistent
here, and we see that, for $|x|\gtrsim 2.24$, the time of exit from the
neighbourhood of the repelling QSS (and of the onset of oscillations) is
delayed beyond the Hopf bifurcation time essentially by the amount determined
by the space-time buffer curve, given by (3.13).

In contrast, for $|x|\lesssim x_{*}$, the solution continues to stay near the
repelling QSS for much longer, and there is a delay (beyond $\mu_{\rm
Hopf}(x)$) of approximately $\omega_{0}$ in duration before the oscillations
(rapid blue to red transitions) commence. For example, at $x=0$, the amplitude
crosses zero (transition from yellow to green) near $\mu=2.8$, and the
oscillations set in near $\mu=3.3$.

Figure 19:  For this simulation, $\mu_{\rm
Hopf}(x)=(2.9240...)e^{-\frac{x^{2}}{24}}$ (blue curve). The parameters are
the same as in Figure 18, except $\tilde{a}=0.5$, $\sigma=4$,
$\alpha=\sqrt{7}$, and $d_{I}=0$. Compared to the simulation in Figure 18, the
magnitude of $\mu_{\rm Hopf}(x)$ is larger due to the larger source term
amplitude $\tilde{a}$, the half-width $\sigma$ has increased for a larger
$\sigma$, and the space-time dynamics of the post-onset oscillations have
changed for larger values of $\alpha$. The space-time buffer curve (3.8)
(black curve) is shown only for those $x$ for which $\mu_{\rm Hopf}(x)\ll 1$,
i.e., where the linearisation about $A=0$ is valid.  Figure 20:  ${\rm
Re}(A(x,\mu))$ obtained from (9.1) with Gaussian source term. The QSS has
large amplitude in an interval about $x=0$. The space-dependent Hopf curve is
given by $\mu_{\rm Hopf}(x)=(8.3788...)e^{-\frac{2x^{2}}{3}}$ (blue curve).
The points $x$ where $\mu_{\rm Hopf}(x)=\sqrt{\operatorname{\varepsilon}}$ are
marked by the dots in the upper frame. The curve $\mu_{\rm Hopf}(x)$ extends
to $\mu=8.3788...$, as is shown in the lower frame. By the time $\mu$ reaches
the tip of the Hopf curve, the amplitude of the QSS has decreased (from red to
light green and blue). Inside the central interval, the oscillations set in
after $\mu$ crosses the blue curve. In contrast, outside this interval, the
oscillations set in at $\mu_{\rm stbc}(x)$, to leading order, since $\mu_{\rm
Hopf}(x)$ is negligibly small there (upper frame, where the space-time buffer
curve (black curve) is also super-imposed). The parameters are
$\tilde{a}=1.0$, $\sigma=\frac{1}{4}$, $\operatorname{\varepsilon}=0.01$,
$\omega_{0}=0.5$, $\alpha=0.6$, $d_{R}=3$, and $d_{I}=1$. The initial data at
$\mu_{0}=-1.5$ is
$A_{0}(x)=-\sqrt{\operatorname{\varepsilon}}\frac{I_{G}(x)}{\mu_{0}+i\omega_{0}}$.

Simulations with other values of $\tilde{a}$ (up to and including
$\tilde{a}=2$), $\sigma$, $\alpha$, and $d$ (complex) show similar results
(over all values simulated) for the delay in the onset of oscillations. In the
central regions, the delayed onset occurs beyond the Hopf curve, $\mu_{\rm
Hopf}(x)$, by an amount approximately equal to $\omega_{0}$. Then, outside the
central region, the DHB is given to leading order by the space-time buffer
curve (obtained by linearising the CGL about $A=0$). See Figure 19, where
$\tilde{a}=0.5$, $\sigma=4$, $\alpha=\sqrt{7}$,
$\operatorname{\varepsilon}=0.01$, and the Hopf bifurcation curve is $\mu_{\rm
Hopf}(x)=(2.9240...)e^{-\frac{x^{2}}{24}}$, by (9.5). See also Figure 20,
where $\tilde{a}=1.0$, $\sigma=\frac{1}{4}$,
$\operatorname{\varepsilon}=0.01$, $\alpha=0.6$, $d_{R}=3$, $d_{I}=1$, and the
Hopf bifurcation curve is $\mu_{\rm
Hopf}(x)=(8.3788...)e^{-\frac{x^{2}}{24}}$, by (9.5).

Remark. The first correction term $\mathcal{A}$ to the QSS in (9.2) is

$\mathcal{A}(x,\mu,\operatorname{\varepsilon})=\frac{1}{3(1+\alpha^{2})^{\frac{4}{3}}(I_{a}(x))^{\frac{1}{3}}}[(\mu_{R}-3\alpha^{2}\mu_{R}+4\alpha(\omega_{0}+\mu_{I}))+i(3(\omega_{0}+\mu_{I})-\alpha^{2}(\omega_{0}+\mu_{I})-4\alpha\mu_{R})]+\mathcal{O}(\operatorname{\varepsilon}^{\frac{1}{3}}).$

Remark. The linearised equation (9.3) for $\mathcal{A}$ has $x$-dependent
saddle points (or nilpotent points) in the complex $\mu$-plane. These are
located where ${\rm tr}(M)=0$ and ${\rm det}(M)=0$. Specifically, for values
of $x$ at which $I_{\tilde{a}}(x)$ is strictly of $\mathcal{O}(1)$, these
saddles occur at $\mu_{R}=\mu_{\rm Hopf}(x)$ and $\mu_{I}=({\mu_{I}})_{\pm}$,
where
$({\mu_{I}})_{+}=-\omega_{0}+\frac{8I_{\tilde{a}}(x)^{\frac{4}{3}}}{\operatorname{\varepsilon}^{\frac{2}{3}}(1+\alpha^{2})^{\frac{2}{3}}}-\frac{3\alpha^{2}-1}{8}+\mathcal{O}(\operatorname{\varepsilon}^{\frac{2}{3}}),$
and
$({\mu_{I}})_{-}=-\omega_{0}+\frac{3\alpha^{2}-1}{8}+\mathcal{O}(\operatorname{\varepsilon}^{\frac{2}{3}})$.
Hence, compared to the case of small-amplitude source terms
($\beta=\frac{1}{2}$) for which there is one saddle (at $\mu=-i\omega_{0}$),
the large-amplitude source term creates a second saddle, and the maximum and
spatial dependence of $I_{\tilde{a}}(x)$ determine the saddle locations.
Analysis of the Stokes lines through these saddles, especially where they
intersect the $\mu_{R}$-axis, would enable one to further quantify the DHB in
this region, though the analysis for (9.3) is more complex than it is for
(2.1), where the corresponding matrix $M$ is simpler, with
$M_{11},M_{22}=\mu_{R}$, $M_{12}=-(\mu_{I}+\omega_{0})$, and
$M_{21}=(\mu_{I}+\omega_{0})$.

##  10 An example with $\mathcal{O}(1)$ diffusivity and $\mathcal{O}(1)$
source term

In this section, we extend some of the results of the base case of the PDE
(1.1) to an example of $\mathcal{O}(1)$ diffusivity and $\mathcal{O}(1)$
amplitude source term in (1.1) (i.e., $\gamma=0$ and $\beta=0$). The PDE is

$\operatorname{\varepsilon}A_{\mu}=(\mu+i\omega_{0})A-(1+i\alpha)|A|^{2}A+I_{a}(x)+\hat{d}A_{xx},$
(10.1)

where we now use $\hat{d}$ to denote the diffusivity. The data
$A(x,\mu_{0})=A_{0}(x)$ is bounded with sufficiently many continuous
derivatives, with $\mu_{0}<-\omega_{0}$ again of primary interest.

###  10.1 The $\mathcal{O}(1)$ QSS and $\mu_{\rm Hopf}(x)$ for (10.1)

With $\mathcal{O}(1)$ diffusivity $\hat{d}$ and source terms $I_{a}(x)$, the
attracting and repelling QSS are $\mathcal{O}(1)$, which contrasts with the
$\mathcal{O}(\sqrt{\operatorname{\varepsilon}})$ amplitude of the QSS in the
base case, recall (1.2). For general source terms $I_{a}(x)$, one may use
variation of constants to find the QSS of the linearised equation, followed by
an iterative procedure on the mild form of the PDE to generate the QSS of the
nonlinear PDE.

For example, with the Gaussian $I_{a}(x)=e^{\frac{-x^{2}}{4\sigma}}$, the
attracting QSS of the linearised version of (10.1) on $\mu<-\delta$ is found
(using variation of constants) to be given to leading order by

$\begin{split}A_{\rm
QSS}(x,\mu)&=c_{1}(x,\mu)e^{-\sqrt{\frac{-(\mu+i\omega_{0})}{\hat{d}}}x}+c_{2}(x,\mu)e^{\sqrt{\frac{-(\mu+i\omega_{0})}{\hat{d}}}x},\quad{\rm
with}\\\
c_{1}(x,\mu)&=\frac{1}{2}\sqrt{\frac{\pi\sigma}{-{\hat{d}}(\mu+i\omega_{0})}}e^{-\frac{\sigma(\mu+i\omega_{0})}{\hat{d}}}\left[1+{\rm
erf}\left(\frac{x}{2\sqrt{\sigma}}-\sqrt{\frac{-\sigma(\mu+i\omega_{0})}{\hat{d}}}\right)\right]\\\
c_{2}(x,\mu)&=\frac{1}{2}\sqrt{\frac{\pi\sigma}{-\hat{d}(\mu+i\omega_{0})}}e^{-\frac{\sigma(\mu+i\omega_{0})}{\hat{d}}}\left[1-{\rm
erf}\left(\frac{x}{2\sqrt{\sigma}}+\sqrt{\frac{-\sigma(\mu+i\omega_{0})}{\hat{d}}}\right)\right].\end{split}$
(10.2)

We observe $\lim_{x\to-\infty}c_{1}(x,\mu)=0$ and
$\lim_{x\to\infty}c_{2}(x,\mu)=0$. A similar formula holds for the leading
order repelling QSS on $\mu>\delta$. Moreover, one can show using steepest
descents on the integral in $A_{p}$ (just as for the base case in Section 2),
that general solutions with $\mu_{0}<-\omega_{0}$ stay near the attracting QSS
for $\mu<-\delta$ and then near the repelling QSS at least until $\omega_{0}$
at each point $x$. However, the analysis is more involved, since the QSSs of
the CGL have $\mathcal{O}(1)$ amplitude.

We now find a formula for the $x$-dependent Hopf bifurcation curve $\mu_{\rm
Hopf}(x)$. Let

$A(x,\mu)=A_{\rm
QSS}(x,\mu)+\operatorname{\varepsilon}^{\frac{1}{3}}\mathcal{A}(x,\mu,\operatorname{\varepsilon}),$
(10.3)

where $\mathcal{A}$ is at most $\mathcal{O}(1)$ uniformly in
$\operatorname{\varepsilon}$ for all $x$ and $\mu$. The linearised equation
for $\mathcal{A}$ is

$\operatorname{\varepsilon}\mathcal{A}_{\mu}=(\mu+i\omega_{0})\mathcal{A}-(1+i\alpha)\left(2|A_{\rm
QSS}|^{2}\mathcal{A}+A_{\rm
QSS}^{2}\bar{\mathcal{A}}\right)+{\hat{d}}\mathcal{A}_{xx}.$

In terms of the real and imaginary parts, $A_{\rm QSS}(x)=u_{\rm QSS}+iv_{\rm
QSS}$ and $\mathcal{A}=\mathcal{U}+i\mathcal{V}$, the linearised equation for
$\mathcal{A}$ may be re-expressed as

$\operatorname{\varepsilon}\begin{bmatrix}{\mathcal{U}}_{\mu}\\\
{\mathcal{V}}_{\mu}\end{bmatrix}=M\begin{bmatrix}\mathcal{U}\\\
\mathcal{V}\end{bmatrix}+{\hat{d}}\begin{bmatrix}{\mathcal{U}}_{xx}\\\
{\mathcal{V}}_{xx}\end{bmatrix},$ (10.4)

where $\mu=\mu_{R}+i\mu_{I}$ and

$M=\begin{bmatrix}\mu_{R}-3u_{\rm QSS}^{2}-v_{\rm QSS}^{2}+2\alpha u_{\rm
QSS}v_{\rm QSS}&-(\mu_{I}+\omega_{0})+\alpha u_{\rm QSS}^{2}+3\alpha v_{\rm
QSS}^{2}-2u_{\rm QSS}v_{\rm QSS}\\\ (\mu_{I}+\omega_{0})-3\alpha u_{\rm
QSS}^{2}-\alpha v_{\rm QSS}^{2}-2u_{\rm QSS}v_{\rm QSS}&\mu_{R}-u_{\rm
QSS}^{2}-3v_{\rm QSS}^{2}-2\alpha u_{\rm QSS}v_{\rm QSS}\end{bmatrix}.$

Now, the trace of $M$ is

${\rm tr}(M)=2\mu_{R}-4(u_{\rm QSS}^{2}+v_{\rm QSS}^{2}).$ (10.5)

Hence, for the solutions of (10.1), the Hopf bifurcation is given implicitly
by

$\mu_{\rm Hopf}(x)=2|A_{\rm QSS}|^{2},$ (10.6)

where $A_{\rm QSS}=u_{\rm QSS}+iv_{\rm QSS}$ is evaluated at $(x,\mu_{\rm
Hopf}(x))$. This result for $\mathcal{O}(1)$ diffusivity and amplitude source
term shows that in the regions where $\mu_{\rm Hopf}(x)$ is $\mathcal{O}(1)$
the QSS changes from being attracting to repelling at a value of $\mu$
substantially to the right of zero, and the DHB needs to be determined from
(10.4). In contrast, for those $|x|$ at which $\mu_{\rm Hopf}(x)\ll 1$, this
formula shows that the changeover in the stability type of the QSS happens
again at $\mu=0$ to leading order, and the linearisation about $A=0$ is again
valid.

###  10.2 The space-time buffer curve for the linearised version of (10.1)

The general solution of the PDE obtained by linearising (10.1) about $A=0$ is
decomposed into homogeneous and particular components,
$A(x,\mu)=A_{h}(x,\mu)+A_{p}(x,\mu).$ These two components are derived using
Duhamel’s Principle (just as in Section 2.1) to solve the equation for $B$.
Here,

$A_{h}(x,\mu)=\frac{\sqrt{\operatorname{\varepsilon}}\ \ {\rm
exp}[\frac{1}{2\operatorname{\varepsilon}}((\mu+i\omega_{0})^{2}-(\mu_{0}+i\omega_{0})^{2})]}{\sqrt{4\pi\hat{d}(\mu-\mu_{0})}}\int_{R}{\rm
exp}\left[\frac{-\operatorname{\varepsilon}(x-y)^{2}}{4\hat{d}(\mu-\mu_{0})}\right]A_{0}(y)dy.$
(10.7)

Note that with $\hat{d}=\operatorname{\varepsilon}d$ in (10.7), one naturally
recovers the formula for the homogeneous solution (4.1) in the base case.

Next, in the scaled variable $\hat{x}=\sqrt{\operatorname{\varepsilon}}x$, the
equation for $B_{p}$ is

$\operatorname{\varepsilon}(B_{p})_{\mu}=\operatorname{\varepsilon}\hat{d}(B_{p})_{\hat{x}\hat{x}}+I_{a}(\hat{x}/\sqrt{\operatorname{\varepsilon}})e^{-\frac{1}{2\operatorname{\varepsilon}}(\mu+i\omega_{0})^{2}},$
(10.8)

with $B_{p}(\hat{x}/\sqrt{\operatorname{\varepsilon}},\mu_{0})=0$. The
solution is

$\begin{split}B_{p}\left(\frac{\hat{x}}{\sqrt{\operatorname{\varepsilon}}},\mu\right)&=\
\frac{1}{\operatorname{\varepsilon}}\int_{\mu_{0}}^{\mu}g(\hat{x},\mu-{\tilde{\mu}})e^{-\frac{1}{2\operatorname{\varepsilon}}({\tilde{\mu}}+i\omega_{0})^{2}}d{\tilde{\mu}},\quad{\rm
with}\\\
g(\hat{x}=\sqrt{\operatorname{\varepsilon}}x,\mu-\tilde{\mu})&=\frac{1}{\sqrt{4\pi\hat{d}(\mu-\tilde{\mu})}}\int_{R}e^{\frac{-(\hat{x}-\hat{y})^{2}}{4{\hat{d}}(\mu-\tilde{\mu})}}I_{a}\left(\frac{\hat{y}}{\sqrt{\operatorname{\varepsilon}}}\right)d{\hat{y}}=\frac{\sqrt{\operatorname{\varepsilon}}}{\sqrt{4\pi\hat{d}(\mu-\tilde{\mu})}}\int_{R}e^{\frac{-\operatorname{\varepsilon}(x-y)^{2}}{4\hat{d}(\mu-\tilde{\mu})}}I_{a}(y)dy.\end{split}$
(10.9)

Compare to (2.6).

Now, using the method of stationary phase along the contour $C=C_{1}\bigcup
C_{2}\bigcup C_{3}\bigcup C_{4}$ (recall Section 3), one finds that
$B_{p}=\frac{1}{\sqrt{\operatorname{\varepsilon}}}\sqrt{2\pi}g(\sqrt{\operatorname{\varepsilon}}x,\mu+i\omega_{0})+\mathcal{O}(1)$,
with the dominant contributions again coming from the segments of $C_{2}$ and
$C_{3}$ near the saddle. Hence, translating $B_{p}(x,\mu)$ via (2.2) to obtain
$A_{p}(x,\mu)$, one arrives at the following implicit formula for the space-
time buffer curve:

$\left\\{(x,\mu_{\rm stbc}(x))|\quad{\rm
Re}\left(\ln\left(\sqrt{\frac{2\pi}{\operatorname{\varepsilon}}}g(\sqrt{\operatorname{\varepsilon}}x,\mu_{\rm
stbc}(x)+i\omega_{0})\right)+\frac{1}{2\operatorname{\varepsilon}}(\mu_{\rm
stbc}(x)+i\omega_{0})^{2}\right)=0\right\\}.$ (10.10)

This space-time buffer curve is spatially flatter than (3.8), through the
argument $\sqrt{\operatorname{\varepsilon}}x$ of $g$.

The above analysis of the space-time buffer curve for the linearised version
of (10.1) with $\mathcal{O}(1)$ amplitude source term and $\mathcal{O}(1)$
diffusivity may be illustrated using the Gaussian source term
$I_{G}(x)=e^{-\frac{x^{2}}{4\sigma}}$. For $\mu$ on $C_{4}$, one finds

$g(x,\mu-{\tilde{\mu}})=\sqrt{\frac{\operatorname{\varepsilon}\sigma}{\hat{d}(\mu+i\omega_{0})+\operatorname{\varepsilon}\sigma}}e^{\frac{-\operatorname{\varepsilon}x^{2}}{4({\hat{d}}(\mu+i\omega_{0})+\operatorname{\varepsilon}\sigma)}}+\mathcal{O}(\operatorname{\varepsilon}).$
(10.11)

Then, substitution of (10.11) into (10.10) shows that the leading order space-
time buffer curve is

$\begin{split}\mu^{2}=\omega_{0}^{2}&+\frac{\operatorname{\varepsilon}^{2}x^{2}({\hat{d}}_{R}\mu-{\hat{d}}_{I}\omega_{0}+\operatorname{\varepsilon}\sigma)}{2\left(({\hat{d}}_{R}\mu-{\hat{d}}_{I}\omega_{0}+\operatorname{\varepsilon}\sigma)^{2}+({\hat{d}}_{R}\omega_{0}+{\hat{d}}_{I}\mu)^{2}\right)}\\\
&-\operatorname{\varepsilon}\ln(2\pi\sigma)+\frac{\operatorname{\varepsilon}}{2}\ln\left(({\hat{d}}_{R}\mu-{\hat{d}}_{I}\omega_{0}+\operatorname{\varepsilon}\sigma)^{2}+({\hat{d}}_{R}\omega_{0}+{\hat{d}}_{I}\mu)^{2}\right),\quad\mu\geq\omega_{0}.\end{split}$
(10.12)

An example of the space-time buffer curve obtained by solving (10.10)
numerically with $g$ given by (10.11) and a Gaussian source term is plotted in
Figure 21. (Note that with $\hat{d}=\operatorname{\varepsilon}d$ it reduces to
(3.13) obtained in the base case with
$\mathcal{O}(\operatorname{\varepsilon})$ diffusivity.) Compared to (3.13),
this space-time buffer curve for $\mathcal{O}(1)$ diffusivity defines a
spatially flatter space-time buffer curve $\mu_{\rm stbc}(x)$, because the
diffusivity is one order of magnitude larger here. Note that the simulation
presented in Figure 21 is on the domain $[-200,200]$, which is significantly
larger than that in Figure 5. Hence, the larger the modulus of the
diffusivity, the smaller the magnitude of the spatial contribution, and the
more uniform the delay time becomes. Also, the half-width $\sigma$ of the
Gaussian has a weaker impact.

Figure 21:  ${\rm Re}(A(x,\mu))$ for the solution of (10.1) on the spatial
interval $[-100,100]$. The black space-time buffer curve is obtained by
solving (10.10) numerically with $g$ given by (10.11), and the blue Hopf
bifurcation curve is obtained from formula (10.6) and (10.2). The
$\mathcal{O}(1)$ QSS is green in most of the region to the left of the space-
time buffer curve, where it has small amplitude; and the QSS is red in the
center, where it has maximum amplitude. The parameters are
$\operatorname{\varepsilon}=0.01,\omega_{0}=\frac{1}{2},{\hat{d}}_{R}=1,{\hat{d}}_{I}=0,\alpha=0.6$,
and $\mu_{0}=-1$, with $I_{G}(x)=e^{-\frac{x^{2}}{4\sigma}}$, and
$\sigma=\frac{1}{4}$.

##  11 Conclusions and discussion

###  11.1 Conclusions

Considering the prototypical CGL PDE (1.1) as an equation in its own right,
this article has presented a study of the phenomenon of delayed Hopf
bifurcation (DHB) as the parameter $\mu$ increases slowly in time through an
instantaneous Hopf bifurcation at $\mu=0$. It has been shown that solutions
with initial data given at $\mu_{0}\leq-\omega_{0}$ are not only near the
attracting QSS while $\mu<0$, but they remain near the QSS as $\mu$ continues
to evolve slowly until well after it has become repelling and at least until
$\mu$ reaches $\omega_{0}$. This analysis of the delay of the Hopf bifurcation
(DHB) was performed by directly using the classical methods of stationary
phase and steepest descents on the linear PDE, based on the topography induced
by the saddle point at $\mu=-i\omega_{0}$, and followed by using an iterative
method for solutions of the nonlinear PDE. Specifically, the nonlinear
analysis is based on an iterative method for the difference between the
solution $A$ of the full cubic PDE and the particular solution $A_{p}$ of the
linear PDE.

Then, with these explicit results, it was shown that there is a competition at
the heart of DHB between two exponentially small terms, one each from the
particular solution $A_{p}(x,\mu)$ and the homogeneous solution
$A_{h}(x,\mu)$, to see which component first ceases to be exponentially small.
The former stops being exponentially small and attains magnitude one along the
space-time buffer curve, $\mu_{\rm stbc}(x)$, and the latter along the
homogeneous exit time curve, $\mu_{h}(x)$. Explicit asymptotic formulas were
derived, and their properties were illustrated with different types of source
terms and initial data, including uni-modal, smoothed step function, and
spatially periodic. Furthermore, in some of the examples, it is possible to
calculate the curves from closed form solutions.

Based on an analysis of different outcomes of the competition between
$\mu_{\rm stbc}(x)$ and $\mu_{h}(x)$, i.e., between which comes first, several
primary cases of DHB were introduced and analyzed. The first threes cases of
DHB are for solutions $A(x,\mu)$ of (1.1) with initial data given at
$\mu_{0}\leq-\omega_{0}$. Here, Case 1 of DHB arises when the duration of the
bifurcation delay is determined at all points by $\mu_{\rm stbc}(x)$, i.e.,
when $\mu_{\rm stbc}(x)<\mu_{h}(x)$ for all $x$. Case 2 of DHB occurs when the
bifurcation delay is determined at some points by $\mu_{\rm stbc}(x)$ and at
others by $\mu_{h}(x)$. Case 3 of DHB arises when the duration of the delay is
determined at all points by $\mu_{h}(x)$, i.e., when $\mu_{h}(x)<\mu_{\rm
stbc}(x)$ for all $x$. Finally, Case 4 of DHB was introduced for solutions of
(1.1) with initial data given at $-\omega_{0}<\mu_{0}<-\delta$, where
$\delta>0$ is small but $\mathcal{O}(1)$. It was shown that, also for these
solutions, the exit time from a neighbourhood of the repelling QSS can be
spatially-dependent, as well.

Examples were presented of the different cases of DHB, and it was shown how to
classify the DHB for general source terms $I_{a}(x)$ and various initial data.
The local maxima of the source term and the initial data mark the sites at
which the solution of the full cubic PDE (1.1) first diverges from the
repelling QSS, and where the large-amplitude, post-DHB oscillations first set
in. The spatial dependence of the DHB and onset of oscillations was shown to
be quadratic in the case of Gaussians (uni-modal functions), a smoothed step
function in the case of source terms given by an error function, and
spatially-periodic in the case of spatially-periodic functions.

Finally, extensions of the main results were presented. Going beyond the main
DHB results established for bounded and positive source terms, it was shown
that DHB also occurs in the base case of (1.1) with algebraically-growing and
sign-changing source terms. The formulas for the space-time buffer and
homogeneous exit time curves (calculated either asymptotically or exactly)
also accurately predict when the oscillations set in at each point $x$, even
though the source terms are not bounded or positive. Next, for the PDE with
asymptotically large source terms ($\beta=-\frac{1}{2}$), it was shown that
the instantaneous Hopf bifurcation curve $\mu_{\rm Hopf}(x)$ can become large,
even asymptotically large, so that the duration of the DHB can be
asymptotically long. Combined with the information derived above about how the
properties of the source terms determine the space-time buffer curve, this
provides a high level of control or ability to design the spatial dependence
of when the oscillations set in. Moreover, with large-amplitude source terms,
it was found that there is more than one saddle point in the complex $\mu$
plane, and hence the topography of the Stokes and anti-Stokes lines is richer.
A final extension concerns the case of $\mathcal{O}(1)$ diffusivity, for which
the space-time buffer curve is also derived and found to be spatially flatter.
That the method also extends to $\mathcal{O}(1)$ diffusivity enables
application to a broader range of problems, in which the diffusivities are not
necessarily small.

There are important considerations about the stability of the numerical
simulations. As we have demonstrated throughout the article, the solutions are
rapidly oscillating with frequency on the order of
$\frac{\mu\omega_{0}}{\operatorname{\varepsilon}}$. The numerical stiffness
induced by the combination of rapid oscillations and slow drift in $\mu$
places an upper bound on the values of $\omega_{0}$ that can be used to
reliably compute the solution near the repelling QSS. On the other hand, our
space-time buffer curve predictions require that the leading order estimate of
the spatially-dependent Hopf bifurcation, $\mu_{\rm
Hopf}(x)=\frac{2\operatorname{\varepsilon}I_{a}(x)}{\omega_{0}^{2}}$, is
small. As such, the numerical values of $\omega_{0}$ that can be used are
bounded from below, for any fixed value of $\operatorname{\varepsilon}$. Thus,
to address the stiffness and satisfy the smallness of $\mu_{\rm Hopf}(x)$, we
have chosen $\omega_{0}$ values in the range
$C\operatorname{\varepsilon}^{1/4}<\omega_{0}<1$, where $C$ is an
$\mathcal{O}(1)$ constant. Many of our reported simulations use
$0.5\leq\omega_{0}\leq 0.75$ for $\operatorname{\varepsilon}=0.01$.

###  11.2 Post-DHB spatio-temporal patterns

Post DHB, several types of spatio-temporal patterns are observed depending on
the type of inhomogeneity, the specific case of DHB, and the dispersion
parameters $\alpha$ and $d_{I}$. For instance, in Case 1 of DHB, Gaussian
source terms break the $x$-translation symmetry in the system, and we find
that it leads to the bifurcation of large-amplitude periodic wave-trains which
organize into stationary, or “pinned,” defects, see Figures 5 and 22.

In the nonlinear CGL (1.1), with the parameter $\mu>0$ held constant (i.e.,
$\mu_{t}=0$), periodic waves have the explicit form $re^{i(kx-\omega t)}$ with
amplitude and nonlinear dispersion relations

$r^{2}=\mu-\operatorname{\varepsilon}d_{R}k^{2},\qquad\qquad\omega(k)=-\omega_{0}+\alpha\mu-\operatorname{\varepsilon}(\alpha
d_{R}-d_{I})k^{2}.$ (11.1)

The amplitude relation shows that periodic patterns exist for wavenumbers
$|k|<\sqrt{\mu/\operatorname{\varepsilon}d_{R}}$. The dispersion relation
shows that the frequency changes sign at $\mu=\frac{\omega_{0}}{\alpha}$,
provided $\alpha\neq 0$, to leading order (with
$0<\operatorname{\varepsilon}\ll 1$ and $k=\mathcal{O}(1)$). From these
relations, the phase and group velocities for patterns with wavenumber $k$
take the form

$c_{p}:=\omega(k)/k,\qquad\qquad c_{g}:=\omega^{\prime}(k)
=2\operatorname{\varepsilon}(d_{I}-\alpha d_{R})k.$

Then, for a dynamic Hopf parameter, we expect amplitudes to vary adiabatically
as $\mu$ is increased, unless a stability boundary is reached, after which we
expect a secondary dynamic bifurcation.

In Figure 5 (where $\alpha=0$, $d_{R}=3,$ and $d_{I}=1$), we observe that the
heterogeneity induces a symmetric defect which connects wavetrains of
wavenumber $k_{+}=k$ and $k_{-}=-k$ on $x>0$ and $x<0$, respectively, with
$k\approx 1.15$. Inward pointing phase curves indicate the arrangement
$c_{p,-}>0>c_{p,+}$ which, since $\omega(k)<0$ for all $\mu$ and
$\operatorname{\varepsilon}\ll 1$, indicates a wavenumber arrangement
$k_{-}<0<k_{+}$. This then implies that the group velocities take the form
$c_{g,-}<0<c_{g,+}$, so that the defect is a source.

To leading order, a fixed wavenumber is selected at the space-time buffer
curve which persists for increasing $\mu\in[0.5,2]$ with fixed oscillation
frequency $\omega(k)$ (note that, with $\alpha=0$, the frequency $\omega(k)$
has no direct dependence on $\mu$) and increasing amplitude $r(k)$ as $\mu$
increases. As the space time buffer curve expands in $x$, additional periods
of the wavetrains are added, leading to a persistent defect solution. We also
remark that in Case 1, where the DHB is governed by the growth of the
particular solution, the large amplitude pattern is independent of small white
noise perturbations of the initial data. This indicates that such dynamic Hopf
heterogeneities could be used in applications to _select_ organised patterns.

In Figure 22 (where $\alpha>0$, $d_{R}=1,d_{I}=0$), the defect dynamics are
somewhat different. Both for $\alpha=0.1$ and $\alpha=0.6$ (see the top
frames), the defect at $x=0$ is a symmetric sink. This may be seen as follows.
For $\alpha=0.1$, we find that $\omega(k)<0$ for
$\mu\in[\mu_{\mathrm{stbc}}(x),\frac{\omega_{0}}{\alpha}]$ (which lies beyond
the edge of the figure), and the phase velocities point inward. Hence, one
readily calculates that, in this case, $c_{g,-}>0>c_{g,+}$, so that the defect
is a symmetric sink. Then, for $\alpha=0.6$, a sink defect is also formed as
$\mu$ increases past $\mu_{\mathrm{stbc}}(x)$ for $x\sim 0$. As $\mu$
increases further, the frequency $\omega(k)$ changes sign at
$\mu\approx\frac{\omega_{0}}{\alpha}=\frac{5}{6}$ to leading order in
$\operatorname{\varepsilon}\ll 1$ so that phase velocities switch direction,
$c_{p,-}<0<c_{p,+}$. However, the defect is still a sink with inward pointing
group velocities $c_{g,-}>0>c_{g,+}$.

We remark that symmetric sink defects are in general transverse heteroclinic
orbits when viewed in the appropriate spatial dynamics formulation, and come
in one-parameter families parameterised by the wavenumber $k$ [57, 64]. Hence,
we expect them to persist under the slow drift of $\mu$. We also mention that
source defects, which select the asymptotic wavenumbers $(k_{-},k_{+})$, are
not transverse heteroclinics and only exist for isolated asymptotic wavenumber
pairs. These solutions are in general robust under the introduction of a small
localised heterogeneity, and it is believed that they become pinned at the
site of the heterogeneity, as we numerically observe here. It thus would be
interesting to more fully understand how source and sink defects are selected
in the presence of a spatio-temporal dynamic Hopf bifurcation.

For $\alpha$ yet larger, say $\alpha=2.5$ and $\alpha=5$, we find that core
and wave-train instabilities begin to arise. Since $d_{I}=0$, the Benjamin-
Feir condition $1+\alpha d_{I}>0$ holds and the band of stable wavenumbers is
given, for stationary $\mu$, by

$ |k|^{2}<\frac{\mu(1+\alpha d_{I})}{3+\alpha
d_{I}+2\alpha^{2}}=\frac{\mu}{3+2\alpha^{2}}.$

Hence, an increase in $\alpha$ brings a decrease in the range of stable
wavenumbers. See [1, 64] for more information on this topic. The effects of
adiabatic changes in $\mu$ on these large-amplitude spatio-temporal patterns
is a topic of future research.

Figure 22:  Brief illustration of how, beyond the space-time buffer curve,
$\alpha$ affects the direction of propagation of the maxima and minima of the
oscillations in the base case with a Gaussian source. (a) $\alpha=0.1$, (b)
$\alpha=0.6$, (c) $\alpha=2.5$, (d) $\alpha=5.0$. In all frames, the space-
time buffer curve (black curve) given by (3.13) is superimposed. The
parameters are
$\operatorname{\varepsilon}=0.01,\omega_{0}=\frac{1}{2},d_{R}=1,d_{I}=0$, and
$\sigma=\frac{1}{4}$. The initial data at $\mu_{0}=-1$ is
$A_{0}(x)=-\sqrt{\operatorname{\varepsilon}}\frac{I_{G}(x)}{\mu_{0}+i\omega_{0}}$.

###  11.3 Open questions

The analysis and numerical simulations presented in this article raise a
number of open questions and suggest avenues of additional research, beyond
that on the effect of adiabatic variation in $\mu$ on the post-DHB spatio-
temporal patterns, just discussed in Section 11.2.

The formal calculations in Sections 2 and 6 of the attracting and repelling
QSS of the linear and nonlinear systems for $\mu<0$ and $\mu>0$, respectively,
show that they are unique to all orders, but differ in the exponentially small
terms. This suggests studying the Gevrey regularity properties of the QSS.
Moreover, this question is motivated directly by the recent demonstration that
slow manifolds for analytic fast-slow systems of ODEs are Gevrey regular, in
the absence of singularities in the slow flow, see [17], which is a major
extension of the classical results in [24].

Next, we expect that the nonlinear analysis of Section 6, which is only formal
in this article, could be made rigorous and used to prove the convergence of
the sequence of functions $\\{b_{j}\\}$ in an appropriate function space. In
particular, one would try to show that

$e^{\frac{1}{2\operatorname{\varepsilon}}(\mu+i\omega)^{2}}(b_{j+1}-b_{j})\sim
\operatorname{\varepsilon}^{\frac{2j+3}{2}},$

uniformly in $x$ to conclude that $\\{b_{j}\\}_{j\in\mathbb{N}}$ is Cauchy and
hence convergent to a fixed point $b_{*}$ of the mapping $H$ as
$j\rightarrow+\infty$. This limiting fixed point would, by construction, yield
a solution of the mild formulation (6.5) with the desired dynamics for
$\mu\in[\mu_{0},\omega_{0}-\tilde{\delta})$. Standard bootstrapping and energy
estimate methods [42] could be used to lift this mild solution to a classical
solution.

The results in the base case of (1.1) with unbounded and sign-changing source
terms (Section 8) raise the question about how to generalize the analysis of
DHB beyond the assumptions made here. Good agreement is found between the
analytically calculated DHB and numerically calculated exit times from the
neighbourhood of the repelling QSS and the onset of the oscillations, even
though the general conditions in Section 2 are not satisfied with these source
terms.

In the study of DHB for solutions of (1.1), the analyses of $A_{h}$, the
homogeneous component, in Section 4 and of the solution $A$ of the full cubic
PDE (1.1) in Section 6 were done working with $\mu$ on the real line.
Analyticity was used in the study of the particular solution, $A_{p}$, of the
linear PDE, to deform the contour in the complex $\mu$ plane to derive the
space-time buffer curve. It is of interest to determine whether analyticity is
strictly necessary for DHB, or instead if some slightly weaker assumption
might suffice to observe DHB. We refer to sections 2 and 3 of [4], for
example, for aspects of the general theory of Gevrey regularity of dissipative
PDEs with analytic nonlinearities. Also, there is a well-developed theory of
Gevrey regularity for solutions of nonlinear R-D equations on bounded domains,
see for example [53], especially the results for the GL equation in section
2.2 there. This could be useful for developing a better understanding of DHB
in nonlinear R-D equations, the minimal hypotheses necessary for it, and the
Gevrey regularity of the solutions.

There are a number of possible avenues to explore to rigorously determine
whether or not slow invariant manifolds exist on $\mu\in[\mu,-\delta)$ and
$\mu\in(-\delta,-\mu_{0}]$, respectively. We are presently studying this
question for (1.1) using scaling variables on the real line and previous
results on slow manifolds; see for example [34] and references therein.

Finally, the type of competition analyzed herein between the homogeneous and
particular solutions of the linearised problem, in terms of which stops being
exponentially small first and transitions through modulus one to becoming
exponentially large, arises in other PDEs. We are currently exploring it in a
number of other reaction-diffusion systems that exhibit DHB, see [36].

## Acknowledgements

We thank Antonella Cucchetti, Arjen Doelman, Edgar Knobloch, and Gene Wayne,
for useful comments and questions. T.K. and T.V. thank the Lorentz Center
(Leiden University, NL) for their hospitality and stimulating research
environment during the 2017 Workshop on Singular Perturbation Theory. T.K.
thanks the organisers of the February 2019 conference on Advances in Pattern
Formation (Ben-Gurion University, Israel), where the initial results about DHB
in the CGL equation were presented. T.K. is grateful to Prof. G.B. Whitham,
from whom he learned the methods of stationary phase and steepest descents.
This work was supported by the US National Science Foundation [DMS-2006887 to
R.G., DMS-1616064 to T.J.K., and DMS-1853342 to T.V.].

## References

  * [1] I.S. Aranson and L. Kramer [2002], The world of the Complex Ginzburg-Landau equation, Rev. Mod. Phys., 74, 99-143
  * [2] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox [2012], Tipping points in open systems: Bifurcation, noise-induced and rate-dependent examples in the climate system, Phil. Trans. R. Soc. A, 370, 1166-–1184
  * [3] D. Avitabile, M. Desroches, R. Veltz, and M. Wechselberger [2020], Local theory for spatio-temporal canards and delayed bifurcations, SIAM J. Math. An., 52, 5703-5747
  * [4] H. Bae and A. Biswas [2015], Gevrey regularity for a class of dissipative equations with analytic nonlinearity, Meth. Applic. An., 22, 377-408
  * [5] S.M. Baer, T. Erneux, and J. Rinzel [1989], The slow passage through a Hopf bifurcation: delay, memory effects, and resonance, SIAM J. Appl. Math., 49, 55–71
  * [6] W. Balser [2005], Power series solutions of the inhomogeneous heat equation, Recent Trends in Microlocal Analysis, Kokyuroku RIMS, 1412, 151-159
  * [7] E. Barreto and J.R. Cressman [2011], Ion concentration dynamics as a mechanism for neuronal bursting, J. Biol. Phys. 37, 361-–373
  * [8] C.M. Bender and S.A. Orszag [1999], Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory, Springer-Verlag, New York
  * [9] F. Berthier, J.-P. Diard, and S. Nugues [1997], On the nature of the spontaneous oscillations observed for the Koper-Sluyters electrocatalytic reaction, J. Electroanalytical Chem., 436, 35–42
  * [10] R. Bertram, M. Butte, T. Kiemel, and A. Sherman [1995], Topological and phenomenological classification of bursting oscillations, Bull. Math. Biol., 57, 413–-439
  * [11] L.M. Bilinsky and S. M. Baer [2018], Slow passage through a Hopf bifurcation in excitable nerve cables: Spatial delays and spatial memory effects, Bull. Math. Bio., 80, 130-–150
  * [12] C. Borgers, M. Krupa, and S. Gielen [2010], The response of a classical Hodgkin-Huxley neuron to an inhibitory input pulse, J. Comp. Neurosci., 28, 509–526.
  * [13] B. Braaksma [1993], Critical phenomena in dynamical systems of van der Pol type, Ph.D. Thesis, Utrecht University, NL
  * [14] S.J. Chapman, J.R. King, and K.L. Adams [1998], Exponential asymptotics and Stokes lines in nonlinear ordinary differential equations, Proc. Roy. Soc. Lond. A, 454, 2733-2755
  * [15] S.J. Chapman and D.B. Mortimer [20015], Exponential asymptotics and Stokes lines in partial differential equations, Proc. Roy. Soc. Lond. A, 461, 2385-2421
  * [16] C.A. Del Negro, C.-F. Hsiao, S.H. Chandler, and A. Garfinkel [1998], Evidence for a novel bursting mechanism in rodent trigeminal neurons, Biophys. J., 75, 174-–182
  * [17] P. De Maesschalck and K. Kenens [2020], Gevrey properties of slow manifolds, Nonlinearity, 33, 341-387
  * [18] R.B. Dingle [1973], Asymptotic expansions: Their derivation and interpretation, Academic Press, London, UK
  * [19] M. Dolnik, I. Berenstein, A.M. Zhabotinsky, and I.R. Epstein [2001], Spatial periodic forcing of Turing structures, Phys. Rev. Lett., 87, 238301
  * [20] W. Eckhaus [1983], Relaxation oscillations including a standard chase on French ducks, in Asymptotic Analysis II, F. Verhulst (ed.) LNM, 985, Springer, Berlin, 449-494
  * [21] H. Engler, H.G. Kaper, T.J. Kaper, and T. Vo [2017], Dynamical systems analysis of the Maasch-Saltzman model for glacial cycles, Physica D, 359, 1-–20
  * [22] T. Erneux and E. Reiss [1990], Delaying the Transition to Hopf Bifurcation by Slowly Varying the Bifurcation Parameter,” in Spatial inhomogeneities and transient behaviour in chemical kinetics, P. Gray, G. Nicolis, F. Baras, P. Borckmans, and S. Scott, eds., Proc. Nonlin. Sci., Manchester University Press, Manchester, Chap. 18, 267–-278
  * [23] T. Erneux, E. Reiss, L. Holden, and M. Georgiou [1991], Slow passage through bifurcations and limit points: Asymptotic theory and applications, in Dynamic Bifurcations, Lecture Notes in Math., 1493, E. Benoit, ed., Springer, Berlin, 14–-28
  * [24] N. Fenichel [1979], Geometric singular perturbation theory for ordinary differential equations, J. Diff. Eq., 31, 53–98
  * [25] D. Golomb, K. Donner, L. Shachan, D. Shlosberg, Y. Amitai, and D. Hansel, Mechanisms of firing patterns in fast-spiking cortical interneurons, PLOS Comp. Bio., 3, e156
  * [26] M. Golubitsky, K. Josic, and T.J. Kaper [2001], An unfolding theory approach to bursting in fast-slow systems, in Global Analysis of Dynamical Systems, H.W. Broer, B. Krauskopf, and G. Vegter (eds.), Institute of Physics, Bristol, UK, 277–308
  * [27] J. Grasman and J.J. Wentzel [1994], Co-existence of a limit cycle and an equilibrium in Kaldor’s business cycle model and its consequences, J. Econ. Behav. Organization, 24, 369–377
  * [28] L. Haim, A. Hagbereg, and E. Meron [2015], Non-monotonic resonance in a spatially-forced Lengyel-Epstein model, Chaos, 25, 064307
  * [29] X. Han, F. Xia, P. Ji, Q. Bi, and J. Kurths [2016], Hopf-bifurcation-delay-induced bursting patterns in a modified circuit system, Comm. Nonlin. Sci. Num. Sim., 36, 517–527
  * [30] M. Haragus and G. Iooss [2010], Local Bifurcations, Center Manifolds, and Normal Forms in Infinite-Dimensional Dynamical Systems, Springer, London
  * [31] M.G. Hayes, T.J. Kaper, P. Szmolyan, and M. Wechselberger [2016], Geometric desingularization of degenerate singularities in the presence of fast rotation: A new proof of known results for slow passage through Hopf bifurcations, Indag. Math., $\bf 27$, 1184–1203
  * [32] M. Higuera, J. Porter, and E. Knobloch [2008], Faraday waves, streaming flow, and relaxation oscillations in nearly circular containers, Chaos, 18, 015104
  * [33] L. Holden and T. Erneux [1993], Slow passage through a Hopf bifurcation: From oscillatory to steady state solutions, SIAM J. Appl. Math., 53, 1045-–1058
  * [34] F. Hummel and C. Kuehn [2020], Slow manifolds for infinite-dimensional evolution equations, preprint, arxiv:2008.10700
  * [35] E. Izhikevich [2001], Synchronization of elliptic bursters, SIAM Rev., 43, 315-–344
  * [36] T.J. Kaper and T. Vo [2018], Delayed loss of stability due to the slow passage through Hopf bifurcations in reaction-diffusion equations, Chaos, $\bf 28$, 091103
  * [37] J. Kevorkian and J.D. Cole [1981], Perturbation Methods in Applied Mathematics, Applied Mathematical Sciences series, vol. 34, Springer, New York
  * [38] M.T.M. Koper and B.D. Aguda [1996], Experimental demonstration of delay and memory effects in the bifurcations of nickel electrodissolution, Phys. Rev. E, 54, 960–963
  * [39] M.T.M. Koper [1998], Non-linear phenomena in electrochemical systems, J. Chem. Soc. Faraday Trans., 94, 1369–-1378
  * [40] P. Kügler [2016], Early Afterdepolarizations with growing amplitudes via delayed sub-critical Hopf bifurcations and unstable manifolds of saddle-foci in cardiac action potential dynamics, PLOS One, 11, e0151178
  * [41] C. Kuehn [2015], Multiple Time Scale Dynamics, Applied Mathematical Sciences series, vol. 191, Springer, New York
  * [42] C.D. Levermore and M. Oliver [1996]. The complex Ginzburg-Landau equation as a model problem, Lectures in Applied Mathematics 3. 141-190, AMS, Providence, Rhode Island.
  * [43] D.A. Lutz, M. Miyake, and R. Schäfke [1999], On the Borel summability of divergent solutions of the heat equation, Nagoya Math. J., 154, 1-29
  * [44] K.A. Maasch and B. Saltzman [1990], A low-order dynamic model of global climate variability over the full Pleistocene, J. Geophys. Res., 95, D2:1955–1963
  * [45] P.D. Miller [2006], Applied Asymptotic Analysis, AMS Graduate Studies in Applied Mathematics, 75, Providence
  * [46] J.D. Murray [1984], Asymptotic Analysis, 2nd ed. Applied Mathematical Sciences series, vol. 48, Springer, New York
  * [47] A.I. Neishtadt [1987], Persistence of stability loss for dynamical bifurcations. I, Diff. Urav., 23, 2060–2067 (Russian). [English translation: Diff. Eq., 23 (1988) 1385–1391, Plenum Pub. Corp.]
  * [48] A.I. Neishtadt [1988], Persistence of stability loss for dynamical bifurcations. II, Diff. Urav., 24, 226–233 (Russian). [English translation: Diff. Eq., 24 (1988) 171–176, Plenum Pub. Corp.]
  * [49] E.W. Ng and M. Geller [1969], A table of integrals of the error functions, J. Res. Nat. Bur. Standards – B. Mathematical Sciences, 73B (1), 1–20
  * [50] F.W.J. Olver [1974], Asymptotics and special functions, AK Peters, Wellesley, MA.
  * [51] Y. Park, Y. Do, and J.M. Lopez [2011], Slow passage through resonance, Phys. Rev. E, 84, 056604
  * [52] D. Premraj, K. Suresh, T. Banerjee, and K. Thamilmaran [2016], An experimental study of slow passage through Hopf and pitchfork bifurcations in a parametrically-driven nonlinear oscillator, Comm. Nonlin. Sci. Num. Sim., 37, 212-221
  * [53] K. Promislow [1991], Time analyticity and Gevrey regularity for solutions of a class of dissipative partial differential equations, Nonlin. An., Theory, Meth., Applic., 16, 959-980
  * [54] P. Remy [2019], Gevrey index theorem for the inhomogeneous $n$-dimensional heat equation with a power-law nonlinearity and variable coefficients, HAL archives-ouvertes, hal-02117418, preprint.
  * [55] J. Rinzel and S. Baer [1988], Firing threshold of the Hodgkin-Huxley model for a slow current ramp: A memory effect and its dependence on fluctuations, Biophys. J., 54, 551–-555
  * [56] J. Rubin and D. Terman [2002], Geometric singular perturbation analysis of neuronal dynamics, in Handbook of Dynamical Systems, volume 2, B. Fiedler, ed., Elsevier, chapter 3, 93-–146
  * [57] B. Sandstede and A. Scheel [2004], Defects in oscillatory media: toward a classification, SIAM Journal on Applied Dynamical Systems 3, 1-68
  * [58] M.A. Shishkova [1973], A discussion of a certain system of differential equations with a small parameter multiplying the highest derivatives, Dokl. Akad. Nauk SSSR, 209, 576–579
  * [59] P. Shorten and D. Wall [2000], A Hodgkin-Huxley model exhibiting bursting oscillations, Bull. Math. Biol., 62, 695-–715
  * [60] R. Speth, W. Green, S. MacNamara, and G. Strang [2013], Balanced splitting and rebalanced splitting, SIAM J. Numer. Anal., 51, 3084-3105
  * [61] P. Strizhak and M. Menzinger [1996], Slow passage through a supercritical Hopf bifurcation: Time-delayed response in the Belousov-Zhabotinsky reaction in a batch reactor, J. Chem. Phys., 105, 10905
  * [62] J. Su [1993], Delayed oscillation phenomena in the FitzHugh Nagumo equation, J. Diff. Eq., 105, 180–215
  * [63] G. Ullah, Y. Wei, M.A. Dahlem, M. Wechselberger, and S.J. Schiff [2015], The role of cell volume in the dynamics of seizure, spreading depression, and anoxic depolarization, PLOS Comp. Bio., 11, e1004414
  * [64] W. van Saarloos [1992], Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations, Physica D, 56, 302-367
  * [65] T. Vo, R. Bertram, and T.J. Kaper [2020], Multi-mode attractors and spatio-temporal canards, Physica D, 411, 132544
  * [66] M. Wechselberger [2020], Geometric Singular Perturbation Theory Beyond the Standard Form, Frontiers in Applied Dynamical Systems: Reviews and Tutorials, 6, Springer, New York
  * [67] H. Wu, Y. Ye, M. Chen, A. Xu, and B. Bao [2019], Extremely slow passages in a low-pass filter-based memristive oscillator, Nonlinear Dynamics, 97, 2339–2353

##  Appendix A Second parametrisation for the integral $\mathcal{I}_{a2}$ in
Section 2.2

In this appendix, we briefly show that $\mathcal{I}_{a2}$ (2.13) may be
derived also using a second, implicit parametrisation (by $\sigma$) of the
contour $C_{a2}$,

$C_{a2}:\quad-\frac{1}{2}({\tilde{\mu}}+i\omega_{0})^{2}=-\frac{1}{2}(\mu+i\omega_{0})^{2}+\sigma,\quad\sigma\leq
0.$ (A.1)

Here, $\sigma$ increases monotonically along $C_{a2}$, starting from
$-\infty$, hitting $-\frac{1}{2}(\omega_{0}^{2}-\mu^{2})$ at $q_{a}$, and
continuing up to zero at the point $\mu$. Therefore, the integral for
$\mathcal{I}_{a2}$ may be written as

$\mathcal{I}_{a2}=\frac{e^{-\frac{1}{2\operatorname{\varepsilon}}(\mu+i\omega_{0})^{2}}}{\sqrt{\operatorname{\varepsilon}}}\int_{-\infty}^{0}g(x,\mu-{\tilde{\mu}}(\sigma))e^{\frac{\sigma}{\operatorname{\varepsilon}}}\frac{d{\tilde{\mu}}}{d\sigma}d\sigma.$
(A.2)

Now, we fix $\Delta>0$ sufficiently small, independent of
$\operatorname{\varepsilon}$. The dominant contribution to the integral comes
from $\sigma\in(-\Delta,0]$, since the integrand in $\mathcal{I}_{a2}$ is
exponentially small for all $\sigma<-\Delta$, due to the exponential term,
where we recall $\mu\in[-\omega_{0},-\delta]$ is fixed, $g$ is analytic and
hence bounded, and the differential is bounded along $C_{a2}$. Taylor
expanding $g$ and the differential about $\sigma=0$, and using
${\tilde{\mu}}(\sigma)=-i\omega_{0}+(\mu+i\omega_{0})\left[1-\frac{2\sigma}{(\mu+i\omega_{0})^{2}}\right]^{\frac{1}{2}},$
we find

$g(x,\mu-{\tilde{\mu}}(\sigma))=g\left(x,\mu+i\omega_{0}-(\mu+i\omega_{0})\left[1-\frac{2\sigma}{(\mu+i\omega_{0})^{2}}\right]^{\frac{1}{2}}\right)=g\left(x,\frac{\sigma}{\mu+i\omega_{0}}+\mathcal{O}\left(\frac{\sigma^{2}}{(\mu+i\omega_{0})^{3}}\right)\right),$
(A.3)
$\frac{d{\tilde{\mu}}}{d\sigma}=\frac{-1}{(\mu+i\omega_{0})}\left[1-\frac{2\sigma}{(\mu+i\omega_{0})^{2}}\right]^{-\frac{1}{2}}=\frac{-1}{(\mu+i\omega_{0})}\left[1+\frac{\sigma}{(\mu+i\omega_{0})^{2}}+\mathcal{O}\left(\frac{\sigma^{2}}{(\mu+i\omega_{0})^{4}}\right)\right].$
(A.4)

Hence,

$\mathcal{I}_{a2}=\frac{-e^{-\frac{1}{2\operatorname{\varepsilon}}(\mu+i\omega_{0})^{2}}}{\sqrt{\operatorname{\varepsilon}}(\mu+i\omega_{0})}\int_{-\Delta}^{0}g\left(x,\frac{\sigma}{\mu+i\omega_{0}}+\mathcal{O}\left(\frac{\sigma^{2}}{(\mu+i\omega_{0})^{3}}\right)\right)e^{\frac{\sigma}{\operatorname{\varepsilon}}}\left[1+\mathcal{O}\left(\frac{\sigma}{(\mu+i\omega_{0})^{2}}\right)\right]d\sigma.$
(A.5)

Then, by (2.6)(b) and standard properties of the 1-D heat kernel, we observe
that $\lim_{\chi\to 0^{+}}g(x,\chi)=I_{a}(x).$ Therefore, we directly obtain
(2.13), and hence the same formula for $A_{p}$.

Remark. The two parametrisations (2.10) and (A.1) are consistent. For example,
expressing (2.10) in terms of small $|\sigma|$, we find
${\tilde{\mu}}_{R}(\sigma)=\mu-\frac{\mu\sigma}{\mu^{2}+\omega_{0}^{2}}+\mathcal{O}(\sigma^{2})$,
and hence
${\tilde{\mu}}(\sigma)=\mu-\frac{\mu-i\omega_{0}}{\mu^{2}+\omega_{0}^{2}}\sigma+\mathcal{O}(\sigma^{2})$,
via (2.10). This is identical to the expansion of (A.1),
${\tilde{\mu}}(\sigma)=\mu-\frac{\sigma}{(\mu+i\omega_{0})}+\mathcal{O}\left(\frac{\sigma^{2}}{(\mu+i\omega_{0})^{3}}\right).$

##  Appendix B The solutions of (2.1) are also near the QSS for
$\mu\in(-\delta,\delta)$

In this appendix, we briefly show how the asymptotic analysis of the solutions
$A_{p}(x,\mu)$ of the linear CGL (2.1) in the base case, which was carried out
separately for $\mu\in[-\omega_{0},-\delta]$ in Section 2.2 and for
$\mu\in[\delta,\omega_{0}]$, in Section 2.3 extends to $\mu\in(0,\delta)$, and
to $\mu\in(-\delta,0]$, where we recall that $\delta>0$ is $\mathcal{O}(1)$
and small.

First, we consider the case of $\mu\in(0,\delta)$. The contour is again
$C_{r}=C_{r1}\bigcup C_{r2}\bigcup C_{r3}\bigcup C_{r4}$, as in the previous
subsection. Here, the point $q_{r}$ is close to the saddle at $-i\omega_{0}$
since $\mu\in(0,\delta)$, recall Figure 3. The contribution from $C_{r1}$ and
$C_{r2}$ is again

$A_{p}(x,\mu)=\sqrt{\frac{\pi}{2}}\left(g(x,\mu+i\omega_{0})+\mathcal{O}(\sqrt{\operatorname{\varepsilon}})\right)e^{\frac{1}{2\operatorname{\varepsilon}}(\mu+i\omega_{0})^{2}},$
(B.1)

where the dominant contribution comes from the final portion of the rise along
$C_{r2}$, the steepest ascents path up to the saddle. Then, the contributions
along $C_{r3}$ and $C_{r4}$ yield

$\begin{split}A_{p}(x,\mu)&=-\sqrt{\operatorname{\varepsilon}}\frac{I_{a}(x)}{\mu+i\omega_{0}}+\operatorname{\varepsilon}^{\frac{3}{2}}\left(\frac{I_{a}(x)+d(\mu+i\omega_{0}){I_{a}}^{\prime\prime}(x)}{(\mu+i\omega_{0})^{3}}\right)+\mathcal{O}\left(\frac{\operatorname{\varepsilon}^{\frac{5}{2}}}{(\mu+i\omega_{0})^{5}}\right)\\\
&+\left(\sqrt{\frac{\pi}{2}}+c(\mu)\right)\left(g(x,\mu+i\omega_{0})+\mathcal{O}(\sqrt{\operatorname{\varepsilon}})\right)e^{\frac{1}{2\operatorname{\varepsilon}}(\mu+i\omega_{0})^{2}}.\end{split}$
(B.2)

Here,
$c(\mu)=(1+i)\int_{0}^{\sqrt{\frac{\omega_{0}\mu}{\operatorname{\varepsilon}}}}e^{-i{\tilde{\sigma}}^{2}}d{\tilde{\sigma}}$,
which comes from the contribution along $C_{r3}$, the stationary phase path
from emerging from the saddle; recall (2.19). The function $c(\mu)$ increases
monotonically from zero in the limit $\mu\to 0^{+}$ to $\sqrt{\frac{\pi}{2}}$
for $\delta>0$ and for other ${\cal O}(1)$ values of $\mu>0$.

Next, we consider the case of $\mu\in(-\delta,0)$. We again use the contour
$C_{a}$ from Section 2.2, recall Figure 2. Then, $A_{p}(x,\mu)$ is given by
(2.15), since also for $\mu\in(-\delta,0)$ the dominant part comes from the
final segment up to $\mu$ on the real axis.

Finally, we consider the case of $\mu=0$. The contour of integration is the
union of the line segments

$\mathcal{L}_{1}=\\{{\tilde{\mu}}_{I}=-\omega_{0},{\tilde{\mu}}_{R}\leq
0\\}\quad{\rm
and}\quad\mathcal{L}_{2}=\\{{\tilde{\mu}}_{R}=0,-\omega_{0}\leq{\tilde{\mu}}_{I}\leq
0\\}.$ (B.3)

These line segments lie on the horizontal and vertical anti-Stokes lines
$\psi=0$, respectively, through the saddle. Application of the method of
steepest descents along $\mathcal{L}_{1}$ yields
$\left(\sqrt{\frac{\pi}{2}}g(x,i\omega_{0})+\mathcal{O}(\sqrt{\operatorname{\varepsilon}})\right)e^{-\frac{\omega_{0}^{2}}{2\operatorname{\varepsilon}}}$.

Then, the contribution from the steepest ascents path $\mathcal{L}_{2}$
consists of two pieces, one from the lower limit of integration and a second
piece from the upper limit of integration at ${\tilde{\mu}}=0$. The former
cancels out the contribution from $\mathcal{L}_{1}$. Hence, for $\mu=0$, we
find

$A_{p}(x,0)=\sqrt{\operatorname{\varepsilon}}\frac{iI_{a}(x)}{\omega_{0}}+\operatorname{\varepsilon}^{\frac{3}{2}}\left(\frac{iI_{a}(x)-d\omega_{0}I_{a}^{\prime\prime}(x)}{\omega_{0}^{3}}\right)+\mathcal{O}(\operatorname{\varepsilon}^{\frac{5}{2}}).$
(B.4)

Therefore, the solution is also the QSS to all orders at $\mu=0$, just as for
$\mu\in(-\delta,0)$ and for $\mu\in(0,\delta)$.