##  1 Methods

To fulfill the goals of this paper, we develop and make use of newly released
features within the Versatile Integrator for Chemical Evolution (;
Johnson2020), an open-source python package available for Unix system
architectures (for further details, see Appendix LABEL:sec:vice). These
features are designed to handle models such as these with a wide range of
flexibility. We reserve a description of ’s migration algorithm and our
simulation parameters for § 1.2, first describing our sample of star particles
from the hydrodynamic simulation.

Although we make use of a hydrodynamic simulation to drive stellar migration,
we do not make use of its SFH. While our study employs similar techniques as
Minchev2013, ours differs slightly in that they take a single SFH which is
similar, but not identical to that of their N-body+SPH simulation; in this
paper we present a handful of assumptions about the SFH, which do not
necessarily resemble that of . An alternative to modeling radial migration
based on a hydrodynamical simulation is to invoke dynamical arguments;
Schoenrich2009a and Sharma2020 take such an approach. This method, however,
introduces free parameters which then require fitting to data. An advantage of
our technique is that there are no free parameters describing radial migration
introduced to the model; it is unclear how much fitting radial migration
parameters could bias the model into agreement with parts of the data not
involved in the fitting process. Instead we adopt a physically motivated
migration model evolved from cosmological initial conditions. We rely on a
single realization of this numerical model and are thus subject to the
differences between the dynamical history of this simulated galaxy and that of
the Milky Way. However, one could compare the predictions made by our chemical
evolution models when applied to different hydrodynamical simulations, a
direction we plan to pursue in future work. We emphasize that the
hydrodynamical simulation only informs the mixing processes in our models, and
there is no N-body integration involved in our models aside from that which
was used to run the simulation in the first place.

###  1.1 The Hydrodynamical Simulation

In this paper we make use of star particles from the h277 simulation
(Christensen2012; Zolotov2012; Loebman2012; Loebman2014; Brooks2014). Recently
employed to study the stellar age-velocity relationship, a synopsis of its
detailed simulation parameters and cosmological model can be found in § 2 of
Bird2020. We do not repeat these details here, instead focusing on how we vet
the sample of star particles for use in our chemical evolution models.

The parameters of stars that we use in our analysis are the birth and final
radii and the final midplane distance. did not record the exact birth radius
of each star particle; however, each star particle does have an accurate age
at each snapshot. The orbital radii of stars that are sufficiently young in
their first snapshot should be good approximations of their birth radii. We
therefore restrict our sample to those star particles with an age at first
snapshot of $\leq$ 150 Myr, and we adopt their Galactocentric radius at that
time as their birth radius. We have repeated our analysis with a maximum age
at first snapshot of 50 Myr and found similar results, indicating that these
time intervals are short enough to not impact our conclusions. We adopt the
150 Myr interval because it provides a larger number of star particles to
sample from. Of the star particles that remain after imposing this cut, the
oldest has an age of 13.23 Gyr at the present day (i.e. at the simulation’s
final output). Our GCE model can only apply on timescales as long as or
shorter than the full range of ages of the sample of star particles; we
therefore subtract 0.5 Gyr from the formation times of all star particles,
allowing $T$ = 0 in our models to correspond to $T$ = 0.5 Gyr in . As a
consequence, our disc models trace the chemical evolution of the Galaxy out to
a lookback time of $\sim$13.2 Gyr, or a redshift of $z\approx$ 9, placing the
onset of star formation at that time.

We further restrict our sample of star particles to only those with both
formation and final radii of $\rgal\leq$ 20 kpc, and to have formed within
$\absz\leq$ 3 kpc of the disc midplane. These criteria are intended to
restrict our sample to star particles that formed within the disc and can
therefore be described by a disc GCE model. While it is possible that some
star particles formed in a dwarf galaxy happen to satisfy our geometrical cuts
at the star’s formation time, these particles are few in number, and are only
relevant at large $\rgal$ and high ages: a region of parameter space where few
stars form anyway.

Based on a kinematic decomposition performed on the present-day phase space
distribution of the h277 star particles conducted with the analysis.decomp
routine within the pynbody package (pynbody), we classify each star particle
as having thin disc, thick disc, bulge, pseudobulge, or halo-like kinematics.
Details on the decomposition process can be found in Brook2012 and Bird2013.
We include the entire sample in ’s public code base, but only make use of
those with disc-like kinematics in this paper. Our geometric selection yields
3,152,211 star particles in total, 1,751,765 of which have disc-like
kinematics and are included in our sample.

had a weak and transient bar during its evolution, but it does not have one at
$z$ = 0. This is a noteworthy difference between our model and that of
Minchev2013, and by extension the Minchev2014; Minchev2017 studies as well,
because they selected a hydrodynamic simulation of a galaxy specifically so
that it would have a strong bar at $z$ = 0. This could mean that the dynamical
history of our model Galaxy differs from that of Minchev2013 and perhaps the
Milky Way itself. However, the difference is likely within the uncertainties
of the current understanding of the Milky Way’s dynamical history. Although an
investigation of the impact of bar evolution on stellar migration and thus
chemical evolution is outside the scope of this paper, it is an interesting
question that can be probed by simply swapping the data within for another
simulation, then rerunning our numerical models and comparing the results.

Figure 1: Radial distributions of our sample of star particles from h277. In
the top row, we show distributions of final radius in bins of birth radius and
age, and in the bottom row, we show distributions of birth radius in bins of
final radius and age. Each bin in Galactocentric radius is shown in its own
panel, denoted in text at the top of each panel and by vertical black dashed
lines. We colour-code the distributions according to the age of the star
particles, denoted by the legend in the upper left panel. We smooth all
distributions with a box-car width of 0.5 kpc to improve visual clarity. We
omit the distributions for 8 - 10 Gyr old stars born in the 9 - 11 and 11 - 13
kpc bins due to an insufficient number of star particles with which to
calculate the distribution.

In the top row of Fig. 1, we plot the distributions of final radius in bins of
birth radius and age for our sample of star particles. Conversely, the bottom
row shows distributions of birth radii in bins of final radius and age.
Focusing on the top row of panels, we note that for star particles born at any
radius and time, the distribution of final radius is still peaked near the
birth radius, but the peak moves slightly inward with increasing age. The
tails of the distributions toward larger $\rgal$ are nearly age-independent,
while the tails toward smaller $\rgal$ are not. This suggests that radial
migration inward and outward occur on different timescales in h277,
specifically that inward migration is slower than outward migration. By
extension, this suggests that the two may be tied to different physical
processes. Alternatively, it may simply be that stars that migrate to the
outer Galaxy are no longer subject to strong dynamical perturbations while
stars that move inward can still experience strong orbital disturbances.

Focusing on the bottom row of panels in Fig. 1, we note that the modes of the
birth radius distributions show a much stronger dependence on age than the
modes of the final radius distributions. At any Galactocentric radius at the
present day, the youngest stars are overwhelmingly born at comparable radii,
while the oldest stars are overwhelmingly born at smaller radii. This trend is
most noticeable at large $\rgal$. The differences between the final radius and
birth radius distributions can be understood by considering the radial
gradient of stellar surface density: there are more stars at small radius to
move outward than vice versa, so roughly symmetric evolution of
$R_{\text}{final}$ produces strongly asymmetric evolution of
$R_{\text}{birth}$.

Taking $\left|\Delta\rgal\right|\geq$ 500 pc between birth and final radii as
the criterion for migration inward or outward, we find as global percentages
in our sample that 27% of star particles migrated inward, 29% migrated
outward, and the remaining 44% stayed near their birth radius. As one can see
from the top panels of Fig. 1, a large fraction of migration away from birth
radius has already occurred by the time stars are $\sim$2 Gyr old. If the SN
Ia delay-time distribution (DTD) is a $t^{-1.1}\approx t^{-1}$ power-law as
suggested by observational results (e.g. Maoz2012; Maoz2017), then we expect
similar numbers of SN Ia events to occur with delay times between 0.1 - 1 Gyr
and 1 - 10 Gyr. With an extended DTD and the timescales for migration implied
by Fig. 1, SN Ia progenitors can migrate significant distances before
exploding. This effect has largely been neglected by GCE studies to date on
the grounds that radial migration is a slow process, and thus the majority of
nucleosynthesis should occur near a star’s birth radius (e.g. as assumed in
Minchev2013, and the application of the Weinberg2017 analytic models in
Feuillet2018). However, we show below that radial migration within the
timescale of the SN Ia DTD can have an important impact on some aspects of
chemical evolution.

###  1.2 Radial Migration

As in previous studies (e.g. Matteucci1989; Schoenrich2009a; Minchev2013;
Sharma2020), in this paper we model the Milky Way as a series of concentric
rings111 For clarity, we use the term “ring” to refer to a computational zone
of our calculation (100 pc in radial range) and the term “annulus” to refer to
a larger radial range (typically 2 kpc).  with a uniform width $\Delta\rgal$.
To run numerical simulations of these models, we develop and make use of ’s
milkyway object, designed specifically for such an approach. The milkyway
object is a subclass of a more general object named multizone; at its core a
multizone object is an array of singlezone objects, which are designed to
handle one-zone models of GCE and were the focus of Johnson2020, ’s initial
release paper. The multizone object affords users full control over which zone
any individual stellar population is in at all timesteps following its
formation as well as the ability to move gas between any two zones with any
time dependence. In principle this should allow for arbitrarily complex zone
configurations and migration prescriptions. The milkyway object is a user-
friendly extension of the multizone base class, which enforces an annular zone
configuration as we take here. As defaults, it adopts the stellar migration
model detailed in this section, our star formation law discussed in § 1.6, and
the scaling of the outflow mass loading factor $\eta$ with radius
parameterized in § 1.4.

As in hydrodynamical simulations, star particles in are stand-ins for entire
stellar populations. They are said to be in a given zone if their radius is
between the inner and outer edges of the ring. At all times, their
nucleosynthetic products and returned envelopes are placed in the ISM of the
ring that they are in at that time. forms a fixed number of stellar
populations per zone per timestep, and it allows their masses to vary to
account for variations in the SFR. The total mass of stars formed in a given
zone and timestep is divided evenly among the corresponding stellar
populations, which can then experience different stellar migration histories.

The final radius of a stellar population is then determined based on the birth
and final radii of star particles in the hydrodynamical simulation. Describing
the Galaxy as a series of concentric rings, ’s milkyway object assumes stellar
populations are born at the centres of each ring. For a stellar population
born at a time $T$ and Galactocentric radius $\rgal$, it first searches for
star particles from h277 that formed at $T\pm$ 250 Myr and $\rgal\pm$ 250 pc.
It then randomly selects a star particle from this subsample to act as an
analogue. This stellar population then adopts the change in orbital radius
$\Delta\rgal$ of its analogue, and moves from its birth radius to the implied
final radius at $T$ = 13.2 Gyr with an assumed time dependence (see below). If
no candidate analogues are found, widens the search to $T\pm$ 500 Myr and
$\rgal\pm$ 500 pc. If still no analogue is found, then it finds the star
particle with the smallest difference in birth radius still within a birth
time of $T\pm$ 500 Myr, and assigns it as the analogue. While this
prescription allows stellar populations to be assigned analogues with
significantly different birth radii, this is only an issue for small $T$ and
large $\rgal$ where there are few star particles from , and where few stars
form in nature anyway due to the inside-out growth of galaxies (e.g.
Bird2013). Furthermore, due to the similarity of the histograms in the top row
of Fig. 1, we expect taking $\Delta\rgal$ from a star particle that formed at
a similar time but different birth radius in these instances to be accurate
enough for our purposes. When an star particle is assigned as an analogue, it
is not thrown out of the sample of candidate analogues, in theory allowing a
star particle to act as an analogue for multiple stellar populations. Because
these populations will have similar $R_{\text}{form}$ and $T_{\text}{form}$,
they will have similar but not identical abundances.

Figure 2: A diagram illustrating how Galactocentric radius changes with time
for two stellar populations under our migration schema: diffusion (crimson,
solid), linear (lime, dot-dashed), post-process (black, dotted), and sudden
(blue, dashed). Here the initial and final radii and birth times are chosen at
random for illustrative purposes. With the initial and final Galactocentric
radii of a stellar population, its birth time, and one of these assumptions
regarding the time-dependence of radial migration, the Galactocentric radius
at all times is known. Diffusion is our default model.

We neglect radial gas flows in the present paper (Bilitewski2012; Lacey1985),
instead focusing our comparisons on the time-dependence with which stars
migrate. We consider four models describing the evolution between a star
particle’s birth and final radius:

  * •

Post-Processing: Stars stay where they are born until the final timestep, at
which point they instantly migrate to their final radius. This retains the
assumption that stars do not contribute to nucleosynthesis beyond their birth
radius as employed in previous studies (e.g. Minchev2013). In this scenario,
the ISM of each ring is treated as a one-zone model independent of all other
zones. We illustrate this case with a black dotted line in Fig. 2.

  * •

Sudden: A random number is drawn from a uniform distribution between a stellar
population’s time of birth and the present day. That time is taken to be the
time of instantaneous migration to the present-day annulus. This emulates a
scenario in which a single dynamical interaction rapidly changes a star’s
orbital radius, and it can be thought of mathematically as a generalization of
the post-processing scenario. We illustrate this case with a blue dashed line
in Fig. 2.

  * •

Diffusion: Stars move to their final radii in a continuous, time-dependent
manner, with displacement $\propto\sqrt{\text{age}}$. This scenario
corresponds to diffusion of angular momentum by a random walk, similar to the
assumption used by Frankel2018; Frankel2020. We illustrate this case with a
red solid line in Fig. 2.

  * •

Linear: A simple variation of the diffusion model in which the migration to
the final radius scales linearly with age rather than with
$\sqrt{\text{age}}$. We illustrate this case with a green dot-dashed line in
Fig. 2.

The diffusion model has the clearest physical motivation, and it is our
default assumption used in all cases unless otherwise noted. The other models
provide illustrative contrasts and allow us to investigate our model’s
sensitivity to the details of radial migration. We do not distinguish between
“blurring” and “churning” (Schoenrich2009a), terms frequenly used to refer to
a star’s epicyclic motions and changes in the guiding centre of its orbit,
respectively. Both effects, induced by a wide variety of underlying causes
such as molecular cloud scattering (Mihalas1981; Jenkins1990; Jenkins1992),
orbital resonances with spiral arms or bars (Sellwood2002; Minchev2011), and
satellite perturbations (Bird2012), are present in the simulation.

Our GCE model assumes that the star-forming ISM is vertically and azimuthally
mixed within each radial annulus. The abundances assigned to a stellar
population depend on its birth radius and time, but they do not depend on its
distance from the plane. Nonetheless, as shown below, the abundance patterns
of our simulation exhibit clear vertical gradients because older stellar
populations have thicker vertical distributions (Bird2020). Radial migration
is also coupled to vertical dynamics in complex ways (Solway2012;
Minchev2012b). We will see that these effects already suffice to explain many
of the observed vertical trends of Milky Way disc abundances.

###  1.3 Nucleosynthetic Yields

We focus our analysis on alpha and iron-peak elements, taking oxygen (O) and
iron (Fe) as the representative cases. The dominant enrichment channels of
interest in our models are thus CCSN and SN Ia (Johnson2019). We would expect
similar results for other alpha (e.g. Ne, Mg, Si) and iron-peak elements
(e.g., Cr, Ni), with quantitative differences reflective of their relative
yields. Odd-$Z$ iron-peak elements (e.g., V, Mn, Co) could behave somewhat
distinctly because of metallicity-dependent yields.

CCSN enrichment happens immediately following the formation of progenitor
stars in . This is an adequate approximation, because the lifetimes of massive
stars are short compared to the relevant timescales for galaxy evolution. For
the most massive stars, the lifetimes are comparable to the timestep size we
adopt in our numerical integrations. This assumption implies a linear
relationship between the CCSN enrichment rate and the SFR:

$\dot{M}_{x}^{\text}{CC}=y_{x}^{\text}{CC}\dot{M}_{\star}$ (1)

where $y_{x}^{\text}{CC}$ is the CCSN yield of some element $x$. Physically,
this quantity represents the mass of an element $x$ ejected to the ISM from
all CCSN events associated with a stellar population in units of the stellar
population’s initial mass. For example, if $y_{x}^{\text}{CC}=0.01$, a
hypothetical 100 stellar population would add 1 of $x$ to the ISM immediately.
In this paper, we adopt $y_{\text}{O}^{\text}{CC}$ = 0.015 and
$y_{\text}{Fe}^{\text}{CC}$ = 0.0012 from Johnson2020, who in turn adopt these
values from Weinberg2017.

SN Ia nucleosynthesis products are injected according to a $t^{-1.1}$ DTD with
a minimum delay time of $t_{\text}{D}$ = 150 Myr. This is the default DTD in ,
which was also adopted by Johnson2020, and is suggested by recent
observational results comparing the cosmic SN Ia rate to the cosmic SFH
(Maoz2012; Maoz2017). In a one-zone model at times $t>t_{\text}{D}$, the
enrichment rate of some element $x$ can be expressed as the product of some
yield $y_{x}^{\text}{Ia}$ and the SFH weighted by the DTD:

$\displaystyle\dot{M}_{x}^{\text{Ia}}$
$\displaystyle=y_{x}^{\text{Ia}}\langle\dot{M}_{\star}\rangle_{\text{Ia}}$
(2a)

where $R_{\text{Ia}}$ is the DTD itself, which has units of $M_{\odot}^{-1}$
Gyr-1. Like the CCSN yield, $y_{x}^{\text{Ia}}$ is the mass of some element
$x$ produced by SNe Ia over the time interval $t_{\text{D}}-t_{\text{max}}$,
in units of the stellar population’s initial mass. It can also be expressed as
an integral over the DTD:

$y_{x}^{\text{Ia}}=m_{x}^{\text{Ia}}\int_{t_{\text{D}}}^{t_{\text{max}}}R_{\text{Ia}}(t^{\prime})dt^{\prime}=m_{x}^{\text{Ia}}\frac{N_{\text{Ia}}}{M_{\star}}$
(2c)

where $m_{x}^{\text{Ia}}$ is the average mass of the element $x$ produced in a
single SN Ia event and the integral evaluates to the mean number of SN Ia
events $N_{\text{Ia}}$ per mass of stars formed $M_{\star}$. forces
$t_{\text{max}}$ = 15 Gyr always, though provided one is consistent with
equations eq:mdot_ia and eq:y_x_ia, the results are independent of
$t_{\text{max}}$ because the integrals cancel. Extending this formalism to
multi-zone models is simple; rather than an integral over the star formation
history of a given annulus, the rate becomes a summation over all stellar
populations that are in a given zone at some time:

$\dot{M}_{x}^{\text{Ia}}=y_{x}^{\text{Ia}}\ddfrac{\sum_{i}M_{i}R_{\text{Ia}}(\tau_{i})}{\int_{t_{\text{D}}}^{t_{\text{max}}}R_{\text{Ia}}(t^{\prime})dt^{\prime}}$
(2d)

where $M_{i}$ and $\tau_{i}$ are the mass and age of the $i$’th stellar
population, respectively.

Initially, we adopted $y_{\text{O}}^{\text{Ia}}$ = 0 and
$y_{\text{Fe}}^{\text{Ia}}$ = 0.0017 from Johnson2020, who in turn adopt these
values from Weinberg2017. However, we found that the e-folding timescales of
star formation in our models are sufficiently long (see discussion in § 1.5)
that our fiducial, inside-out SFH model predicted [O/Fe] $\approx$ +0.05 for
young stars. We therefore multiply $y_{\text{Fe}}^{\text{Ia}}$ by $10^{0.1}$,
adopting $y_{\text{Fe}}^{\text{Ia}}$ = 0.00214 so that our fiducial model
predicts a late-time [O/Fe] ratio in better agreement with observations.
Changes at this level are within the uncertainties of SN Ia rates and yields,
so we consider it reasonable to adjust the yields empirically to reproduce
observed abundances.

Our IMF-averaged O and Fe yields are based on a Kroupa2001 IMF combined with
supernova nucleosynthesis models in which most $M>8~{}\msun$ stars explode as
a CCSN (e.g. Chieffi2004; Chieffi2013). Recent studies have strongly suggested
that many high mass stars instead collapse directly to a black hole (see
theoretical discussion by, e.g., Pejcha2015; Sukhbold2016; Ertl2016, and
observational evidence from Gerke2015; Adams2017; Basinger2020). Our yields
would be lower in a scenario with extensive black hole formation and/or a
steeper high mass IMF (Griffith et al. 2021). These effects could also
introduce a metallicity dependence if the landscape of black hole formation
changes with metallicity. The strong increase of the specific SN Ia rate seen
at low galaxy masses (Brown2019) provides circumstantial evidence of a higher
$R_{\text{Ia}}$ normalization at low metallicity; the stellar close binary
fraction also depends on metallicity (Badenes2018; Moe2019). However, there is
presently no solid empirical basis for adopting metallicity dependent O and Fe
yields over the range relevant to this paper (roughly -0.8 $\leq$ $\leq$ 0.4),
and some empirical evidence that any metallicity trends in this range are weak
(Weinberg2019). has capabilities to compute IMF-averaged yields for a flexible
description of the massive star explodability landscape, as described by
Griffith et al. (2021) and the documentation, but we do not use this
methodology here. We expect that most of our results would be largely
unchanged if we were to lower all three yields ($y_{\text{O}}^{\text{CC}}$,
$y_{\text{Fe}}^{\text{CC}}$, $y_{\text{Fe}}^{\text{Ia}}$) by the same factor
and adjust our adopted outflow mass loading prescription to compensate (see
below).

Both AGB star enrichment and the recycling of previously produced metals in
this paper proceed as they did in Johnson2020, with the caveat that the mass
is added to the ring that a stellar population is in at a given time, which
may or may not be the ring it was born in. Recycling proceeds according to the
Kalirai2008 initial-remnant mass relation assuming a Kroupa2001 IMF and a
mass-lifetime relationship of $\tau=1.1\tau_{\odot}(M/M_{\odot})^{-3.5}$,
where $\tau_{\odot}$ = 10 Gyr is the main sequence lifetime of the sun and the
factor of 1.1 accounts for the post-main sequence lifetime. includes AGB
enrichment in all models; here we adopt the net yields sampled on a table of
stellar initial mass and metallicity from the FRUITY database (Cristallo2011).
However, the AGB star yields of O and Fe are tiny compared to their supernova
yields, so they have negligible impact on the results presented in this paper.

###  1.4 Outflows

Weinberg2017 demonstrate that, to first order, the nucleosynthetic yields of a
given element and the strength of outflowing winds determine the late-time
equilibrium abundance in the ISM, with a secondary dependence on the SFH. We
retain their characterization of outflows here, in terms of a mass-loading
factor $\eta$ describing the ratio of the mass outflow to the SFR:

$\eta\equiv\frac{\dot{M}_{\text{out}}}{\dot{M}_{\star}}.$ (2e)

We adopt a scaling of $\eta$ with $\rgal$ such that the late-time equilibrium
abundance as a function of radius describes a metallicity gradient in
agreement with observations. For a constant SFR, the equilibrium abundance of
an $\alpha$-element, produced by CCSNe with a metallicity-independent yield,
is given by

$Z_{\alpha,\text{eq}}=\frac{y_{\alpha}^{\text{CC}}}{1+\eta(\rgal)-r},$ (2f)

where $r$ is the recycling parameter ($\approx$0.4 for the sake of this
scaling with a Kroupa2001 IMF; see discussion in § 2.2 of Weinberg2017).
Solving for $\eta(\rgal)$ yields

$\eta(\rgal)=\frac{y_{\alpha}^{\text{CC}}}{Z_{\alpha,\text{eq}}}+r-1=\frac{y_{\alpha}^{\text{CC}}}{Z_{\alpha,\odot}}10^{-\text{mode([$\alpha$/H])}(\rgal)}+r-1,$
(2g)

where mode([$\alpha$/H])($\rgal$) denotes the mode of the stellar [$\alpha$/H]
distribution at a radius $\rgal$, which we assume to correspond to the
equilibrium abundance at that radius. Recent studies of the disc metallicity
gradient with APOGEE find values of -0.09 dex/kpc to -0.06 dex/kpc (e.g.
Frinchaboy2013; Hayden2014; Weinberg2019), consistent with earlier studies.
Here we adopt a slope of -0.08 dex/kpc and set mode([$\alpha$/H]) to be
$\sim$+0.3 at $\rgal$ = 4 kpc, producing mode([$\alpha$/H]) $\approx$ 0 at
$\rgal$ = 7 - 9 kpc. This results in the following form for $\eta$ as a
function of Galactocentric radius:

$\eta(\rgal)=\frac{y_{\alpha}^{\text{CC}}}{Z_{\alpha,\odot}}10^{(-0.08\text{
kpc}^{-1})(\rgal-4\text{ kpc})+0.3}+r-1,$ (2h)

where we adopt our CCSN yield of O for $y_{\alpha}^{\text{CC}}$ and the solar
photospheric abundance of O of $Z_{\text{O},\odot}$ = 0.00572 based on
Asplund2009. We plot this adopted scaling in the top panel of Fig. 3,
highlighting a value of $\sim$2.15 for the solar circle with a red dotted
line. In Fig. LABEL:fig:metallicity_gradient below, we show that our full
model including a time-dependent SFH and radial migration produces stellar and
gas phase gradients similar but not identical to those of equation
eq:eta_rgal. In Figs. LABEL:fig:mdf_3panel_fe and LABEL:fig:mdf_3panel_o
below, we show that the model achieves qualitative agreement with the
metallicity gradient observed by APOGEE.

Figure 3: Top: Our implemented scaling of the mass loading factor $\eta$ with
Galactocentric radius (black) as defined by equation eq:eta_rgal. Bottom: The
e-folding timescales of the star formation histories of our model galaxies
(black). These values come from a fit to the $\Sigma_{\star}$-age relation in
bins of $R/R_{\text{e}}$ for $10^{10.5}$ \- $10^{11}$ Sa/Sb Hubble type spiral
galaxies as reported by Sanchez2020. The horizontal and vertical red dashed
lines in both panels highlight a mass loading factor of $\eta\approx$ 2.15 and
a star formation timescale of $\tau_{\text{sfh}}\approx$ 15 Gyr at an assumed
radius of the sun of $R_{\odot}$ = 8 kpc.

###  1.5 Star Formation Histories

computes models in either “infall”, “star formation”, or “gas” mode, referring
to which component of the evolutionary history the user has specified. The
starburst models of Johnson2020 ran in infall mode, meaning that the gas
infall rate is specified and the SFR follows from the gas surface density and
adopted star formation law. Here we run in star formation mode so that we
achieve a specified form of the SFHs in our models. In Appendix
LABEL:sec:normalize_sfh, we explain how we normalize the parameters of our
SFHs to produce a realistic model Galaxy at the present day. In short, we take
a unitless description of the time-dependence of the SFH at a given
Galactocentric radius, denoted $f(t|\rgal)$, and a unitless description of the
present day stellar surface density gradient, denoted $g(\rgal)$. We integrate
$f(t|\rgal)$ with time for each annulus, assuming $\rgal$ to correspond to the
centre of the zone, and attach a prefactor to $f(t|\rgal)$ in each annulus
such that the desired gradient is achieved with a total stellar mass similar
to the Milky Way. This procedure neglects the impact of radial migration,
assuming that it does not significantly alter the form of $g(\rgal)$. We
demonstrate that this assumption holds in § 1.7, in which we also detail our
adopted form of $g(\rgal)$. The equation derived in Appendix
LABEL:sec:normalize_sfh can be used to calculate these prefactors for
alternative models of Milky Way-like galaxies.

In the present paper, we consider four forms of the SFH, which we dub
“constant”, “inside-out”, “late-burst”, and “outer-burst”. They are defined as
follows:

  * •

Constant: The SFH at a given radius is time-independent,

$f_{\text{C}}(t|\rgal)=1.$ (2i)

This case is of theoretical interest because it quantifies the effect of
stellar migration while removing the impact of a time-dependent SFH.

  * •

Inside-Out: This is our fiducial SFH,

$f_{\text{IO}}(t|\rgal)=(1-e^{-t/\tau_{\text{rise}}})e^{-t/\tau_{\text{sfh}}}.$
(2j)

We adopt this mathematical form over the somewhat more common
$te^{-t/\tau_{\text{sfh}}}$ because it allows separate control over the rising
and falling phase of the SFH. Equation eq:insideout_sfh has a maximum near
$\tau_{\text{rise}}$, which in this paper we set to 2 Gyr at all radii. This
form produces a peak in star formation at lookback times of $\sim$11 Gyr,
roughly corresponding to a redsfhit of $z\approx$ 2.5. We discuss the choice
of $\tau_{\text{sfh}}$ below.

  * •

Late-Burst: In this model, the inside-out SFH is modified to exhibit a recent,
slow “burst” in star formation described by a Gaussian in time:

$f_{\text{LB}}(t|\rgal)=f_{\text{IO}}(t|\rgal)\left(1+A_{b}e^{-(t-t_{b})^{2}/2\sigma_{b}^{2}}\right),$
(2k)

where $A_{b}$ is a dimensionless parameter describing the strength of the
starburst, $t_{b}$ is the time of the local maximum in the SFH during the
burst, and $\sigma_{b}$ is the width of the Gaussian describing it. To
approximate the findings of Mor2019 and Isern2019, we adopt $A_{b}$ = 1.5,
$t_{b}$ = 11.2 Gyr, and $\sigma_{b}$ = 1 Gyr, finding that $A_{b}$ = 1.5 with
a declining $f_{\text{IO}}(t|\rgal)$ produces a local maximum SFR that is a
factor of $\sim$2 larger than the preceding local minimum.

  * •

Outer-Burst: A variation of the late-burst model in which only zones at
$\rgal\geq$ 6 kpc experience the starburst, with inner regions following the
inside-out SFH. Because the empirical evidence for elevated recent star
formation comes from local observations, it is useful to investigate the case
where it is confined to the outer Galaxy. In their hydrodynamical simulation
of a Milky Way-like galaxy, Vincenzo2020 find that satellite perturbations
enhance gas accretion preferentially in the outer regions.

The inside-out model is our fiducial choice, and it is the model shown in
figures unless otherwise specified. The constant SFR model allows us to
investigate the impact of migration when the time-dependence of the SFH is
removed, and the two burst models allow us to explore the consequences of
elevated recent star formation supported by some recent data. More complex
scenarios with multiple bursts induced by repeated satellite passages (e.g.
Lian2020a; Lian2020b; RuizLara2020; Sysoliatina2021) can be modeled easily in
, but we do not explore them here.

We derive the $\tau_{\text{sfh}}-\rgal$ relation from the data of Sanchez2020,
who presents the stellar surface density $\Sigma_{\star}$ as a function of age
in bins of $R/R_{\text{e}}$ for MaNGA galaxies (Bundy2015), where
$R_{\text{e}}$ is the half-light radius. Here we take the
$M_{\star}=10^{10.5}-10^{11}~{}\msun$ bin for Sa/Sb spirals and simultaneously
fit the normalization and e-folding timescale $\tau_{\text{sfh}}$ of our
$f_{\text{IO}}(t|\rgal)$ form to the data. Although the normalization is
irrelevant to our models and determined via the method outlined in Appendix
LABEL:sec:normalize_sfh, we adopt the resulting $\tau_{\text{sfh}}-\rgal$
relation in our models. Our adopted stellar surface density gradient (see §
1.7) implies a present-day half-mass radius near 4 kpc. The findings of
Garcia-Benito2017 and GonzalezDelgado2014 suggest that half-light radii are
marginally larger than half-mass radii. Based on equation (4) of
GonzalezDelgado2014 relating the two for circular apertures, we expect our
model Galaxy to have a half-light radius near 5 kpc. We therefore adopt
$R_{\text{e}}$ = 5 kpc to convert the $\tau_{\text{sfh}}-\rgal/R_{\text{e}}$
relation resulting from our fit to the Sanchez2020 data into a
$\tau_{\text{sfh}}-\rgal$ relation.

We illustrate this relationship in the bottom panel of Fig. 3. The resulting
timescales are long, particularly for the outer Galaxy. With a red dotted
line, we highlight a value of $\tau_{\text{sfh}}\approx$ 15 Gyr at an assumed
orbital radius of the sun of $R_{\odot}$ = 8 kpc. The long timescales reflect
the fact that the $\Sigma_{\star}(\tau_{\text{lookback}})$ profiles in Fig. 11
of Sanchez2020 are fairly flat. For comparison, we have also considered the
assumption that the Galactic SFH may have resembled the cosmic SFH by running
models with e-folding timescales of a $\sim$few Gyr (e.g. Madau2014). We find
similar results in these cases, suggesting that our qualitative conclusions
are not sensitive to the exact values of $\tau_{\text{sfh}}$.

Figure 4: The surface densities of star formation $\dot{\Sigma}_{\star}$ (top
row), infall $\dot{\Sigma}_{\text{in}}$ (middle row), and gas
$\Sigma_{\text{g}}$ (bottom row) as functions of simulation time for our four
fiducial SFHs: constant (far left), inside-out (left middle), late-burst
(right middle), and outer-burst (far right). We plot curves for the rings
whose inner radii are 3 kpc (grey), 5 kpc (black), 7 kpc (red), 9 kpc
(yellow), 11 kpc (green), 13 kpc (blue), and 15 kpc (purple) (see equations
eq:constant_sfh, eq:insideout_sfh, and eq:lateburst_sfh for the mathematical
definition of each SFH).

We plot the resulting SFHs of our models in the top row of Fig. 4 for a
handful of radii. Because of the long $\tau_{\text{sfh}}$, the SFR at most
radii has fallen only modestly from its 2 Gyr peak in the inside-out model. At
all timesteps, the gas supply is known via the star formation efficiency
timescale $\tau_{\star}$ (see § 1.6), which is illustrated in the bottom row
of Fig. 4. automatically calculates the implied infall rate by comparing the
amount of gas lost to outflows and star formation in a given timestep to that
which is required to sustain the user-specified level of star formation at the
next timestep; we assume the infalling gas to be of zero metallicity at all
times. This quantity is shown in the middle row of Fig. 4.

###  1.6 Star Formation Efficiency

The term “star formation efficiency” (SFE) is somewhat overloaded in the
literature. In the star formation and feedback community, it usually refers to
the fraction of a molecular cloud’s mass that will eventually be converted
into stars. In the chemical evolution literature, however, it typically refers
to the inverse timescale relating the SFR within some star forming reservoir
to the mass of gas in that region:
$\tau_{\star}\equiv\Sigma_{\text{g}}/\dot{\Sigma}_{\star}$. High (Low) values
of $\tau_{\star}$ indicate slow (fast) conversion of gas and thus low (high)
SFE; when we refer to SFE here, we mean the definition based on this
timescale. In the star formation and feedback literature, $\tau_{\star}$ is
sometimes referred to as the “depletion time”, though here we follow the
terminology of Weinberg2017 who call it the “star formation efficiency
timescale.”

Based on the findings of Kennicutt1998, it is common practice in the chemical
evolution literature to adopt a single power-law describing the relationship
between the surface densities of gas and star formation $\Sigma_{\text{g}}$
and $\dot{\Sigma}_{\star}$, often referred to as the star formation law or the
Kennicutt-Schmidt relation:

$\dot{\Sigma}_{\star}\propto\Sigma_{\text{g}}^{N}.$ (2l)

Kennicutt1998 finds $N=1.4\pm 0.15$ relating the total $\dot{\Sigma}_{\star}$
and $\Sigma_{\text{g}}$ within the disc across a sample of quiescent spiral
galaxies and infrared and circumnuclear starbursts. However, recent studies
have found evidence that much of the observed scatter in this relation is
physical in origin (delosReyes2019) and that there are significant breaks in
both the power-law index and zero-points (Kennicutt2021). Some of the
uncertainty surrounding the details of the star formation law is a consequence
of the ongoing debate about the CO-to-H2 conversion factor (Kennicutt2012;
Liu2015). Although Ellison2021 demonstrate that there are significant galaxy-
to-galaxy variations in the star formation law, delosReyes2019 argue that the
mean trend is still a reasonable recipe for Galaxy evolution models. However,
for our purposes we also need the dependence of $\dot{\Sigma}_{\star}$ on
$\Sigma_{\text{g}}$ within a galaxy.

Krumholz2018a compare theoretically motivated star formation laws to the
observations of Bigiel2010 and Leroy2013 (see their Fig. 2). We find that the
following by-eye fit to the power-law index $N$ is a reasonable description of
the aggregate data:

$N=\begin{cases}1.0&(\Sigma_{\text{g}}\geq\Sigma_{\text{g},2})~{},\\\
3.6&(\Sigma_{\text{g},1}\leq\Sigma_{\text{g}}\leq\Sigma_{\text{g},2})~{},\\\
1.7&(\Sigma_{\text{g}}\leq\Sigma_{\text{g},1})~{},\end{cases}$ (2m)

where $\Sigma_{\text{g},1}=5\times 10^{6}$ and $\Sigma_{\text{g},2}=2\times
10^{7}$ . The apparent linearity of the relationship above $\sim 2\times
10^{7}$ suggests that in this regime, star formation proceeds at the highest
efficiency, and that
$\tau_{\star}\equiv\Sigma_{\text{g}}/\dot{\Sigma}_{\star}$ = constant. The
observational results of Leroy2013 and Kennicutt2021 would suggest that these
are the surface densities at which the molecular fraction
$f_{\text{mol}}=M_{\text{H}_{2}}/(M_{\text{H}_{2}}+M_{\text{HI}})\approx$ 1\.
We therefore adopt the assumption that above $\Sigma_{\text{g}}=2\times
10^{7}$ , $\tau_{\star}$ reaches its minimum value, and increases with
decreasing $f_{\text{mol}}$. We denote this value as $\tau_{\text{mol}}$, the
value of $\tau_{\star}$ for a gas reservoir with $f_{\text{mol}}$ = 1. This
identification, combined with our three-component power-law index $N$ results
in the following final form for our adopted star formation law:

$\dot{\Sigma}_{\star}=\begin{cases}\Sigma_{\text{g}}\tau_{\text{mol}}^{-1}&(\Sigma_{\text{g}}\geq\Sigma_{\text{g},2})~{},\\\
\Sigma_{\text{g}}\tau_{\text{mol}}^{-1}\left(\frac{\Sigma_{\text{g}}}{\Sigma_{\text{g},2}}\right)^{2.6}&(\Sigma_{\text{g},1}\leq\Sigma_{\text{g}}\leq\Sigma_{\text{g},2})~{},\\\
\Sigma_{\text{g}}\tau_{\text{mol}}^{-1}\left(\frac{\Sigma_{\text{g},1}}{\Sigma_{\text{g},2}}\right)^{2.6}\left(\frac{\Sigma_{\text{g}}}{\Sigma_{\text{g},1}}\right)^{0.7}&(\Sigma_{\text{g}}\leq\Sigma_{\text{g},1})~{}.\end{cases}$
(2n)

We choose the power-law indices such that this formalism is consistent with
equation eq:sf_law_indeces, and prefactors are added to ensure piece-wise
continuity. In implementation, requires the
$\tau_{\star}-\dot{\Sigma}_{\star}$ relation when running in star formation
mode and the $\tau_{\star}-\Sigma_{\text{g}}$ relation when running in infall
and gas modes. Both follow algebraically from this relationship given the
substitution $\tau_{\star}\equiv\Sigma_{\text{g}}/\dot{\Sigma}_{\star}$.

Based on the observed Kennicutt-Schmidt relation at different redshifts,
Tacconi2018 suggest that $\tau_{\text{mol}}$ should scale with redshift $z$
and with the deviation from the star forming main sequence $\delta$MS via
$\tau_{\text{mol}}\propto(1+z)^{-0.6}\delta\text{MS}^{-0.44}$. We do not
account for the effect of $\delta$MS in our models, but we do incorporate the
redshift dependence. For redshifts of $z\lesssim$ 3, encompassing most of the
timesteps in our models, a reasonable approximation to the $t-z$ relation
assuming typical $\Lambda$CDM cosmology is given by:

$\frac{t}{t_{0}}\approx(1+z)^{-5/4},$ (2o)

where $t_{0}$ is the present-day age of the universe, and $t$ is not
simulation time but the age of the universe. Plugging this relation into the
Tacconi2018 scaling yields

$\tau_{\text{mol}}=\tau_{\text{mol},0}(t/t_{0})^{12/25}\approx\tau_{\text{mol},0}(t/t_{0})^{1/2},$
(2p)

where $\tau_{\text{mol},0}$ is simply $\tau_{\text{mol}}$ at the present day.
We generalize this formula to the following form:

$\tau_{\text{mol}}=\tau_{\text{mol},0}(t/t_{0})^{\gamma}$ (2q)

In this paper we present models which adopt $\tau_{\text{mol},0}$ = 2 Gyr
(Leroy2008; Leroy2013; Tacconi2018) and $\gamma$ = 1/2 based on this argument.
We have also run simulations witht $\tau_{\text{mol},0}$ = 1 Gyr and with
$\gamma$ = 0 (a time-independent $\tau_{\text{mol}}$), as well as combinations
of the two, and found similar results. In all timesteps and annnuli, infers
$\Sigma_{\text{g}}$ from $\dot{\Sigma}_{\star}$ given our equations eq:sf_law
and eq:tau_mol.

Figure 5: The star formation efficiency timescale $\tau_{\star}$ as a function
of Galactocentric radius at simulation times of 2 Gyr (red), 4 Gyr (yellow), 6
Gyr (green), 8 Gyr (blue), 10 Gyr (purple), and 13.2 Gyr (the present day,
black) predicted by our fiducial model.

In Fig. 5, we plot $\tau_{\star}$ as a function of $\rgal$ at six different
time stamps predicted by our fiducial, inside-out SFH model. At
$\rgal\lesssim$ 6 kpc, $\tau_{\star}$ is near $\tau_{\text{mol}}$ at all
times, implying a molecular fraction of unity at these radii. Although this
prediction is likely unrealistic because 21-cm line observations suggest the
presence of neutral hydrogen as close to the Galactic centre as several
hundred pc (Kalberla2009), we find in practice that changing our assumptions
about the star formation law does not impact our conclusions. In exploratory
work for this paper, we investigated purely linear, purely power-law, and
broken power-law characterizations, finding similar predictions in all cases.
In general we find that the detailed form of the SFH, and to some extent the
time-dependence of radial migration, exert much greater power in establishing
the model predictions than does the star formation law. Nonetheless it is an
interesting puzzle that a $\dot{\Sigma}_{\star}-\Sigma_{\text{g}}$ relation
informed by the observed population-averaged trends and the normalization of
SFHs implied by the stellar mass of the Milky Way predicts results in tension
with the observed HI distribution. Because the star formation law has minimal
impact on our results, we do not pursue this question further here.

###  1.7 Surface Density Gradient

As discussed in § 1.5, Appendix LABEL:sec:normalize_sfh presents the recipe by
which we select a unitless function describing the stellar surface density
gradient $g(\rgal)$. In setting the normalization, our model ensures that the
integral of $g(\rgal)$ over the surface area of the disc predicts a total
stellar mass in agreement with that observed for the Milky Way. For this
value, we adopt $M_{\star}^{\text{MW}}=5.17~{}\times~{}10^{10}~{}\msun$ from
Licquia2015. This is the total disc mass only; when the bulge is included, the
total becomes $(6.08~{}\pm~{}1.17)~{}\times~{}10^{10}~{}\msun$. Since we are
modeling only the disc populations here, we omit the contribution from the
bulge to the total mass budget.

We adopt a double exponential form for $g(\rgal)$, describing the thin and
thick disc components of the Galaxy:

$g(\rgal)=e^{-\rgal/R_{t}}+\frac{\Sigma_{T}}{\Sigma_{t}}e^{-\rgal/R_{T}}~{},$
(2r)

where $R_{t}$ and $R_{T}$ are the scale radii of the thin and thick discs,
respectively, and $\Sigma_{T}/\Sigma_{t}$ is the ratio of their surface
densities at $\rgal$ = 0. We adopt $R_{t}$ = 2.5 kpc, $R_{T}$ = 2.0 kpc, and
$\Sigma_{T}/\Sigma_{t}$ = 0.27 based on the findings of Bland-Hawthorn2016. We
plot the single exponential forms of each disc component as dotted black lines
in Fig. 6, with the solid black line denoting the sum of the two.

Figure 6: The surface density of gas (blue) and stars (red) as predicted by
our fiducial model. The dotted black lines denote thin and thick disc profiles
with scale lengths of $R_{t}$ = 2.5 kpc and $R_{T}$ = 2.0 kpc, respectively,
with a ratio of $\Sigma_{T}/\Sigma_{t}$ = 0.27 at = 0 (i.e. the thin disc
profile has the higher normalization). The solid black line denotes the sum of
the two; this is the stellar surface density gradient of the Milky Way as
reported by Bland-Hawthorn2016, renormalized according to our adopted total
stellar mass of $(5.17\pm 1.11)\times 10^{10}$ (Licquia2015).

We plot the resultant surface density gradients from our fiducial, inside-out
SFH model in Fig. 6 as well, with red denoting the stellar gradient and the
gas in blue. The stellar gradient very nearly follows our target gradient
(equation 2r) denoted by the solid black line. Stellar migration has not
altered the overall form of the gradient, simply introducing scatter around
the adopted trend. Although there are slight enhancements at small $\rgal$ and
deficits at large $\rgal$, the agreement is excellent in the regions of the
Galaxy where the gradient is best constrained observationally. The $\rgal$ ¿
15.5 kpc populations are composed entirely of stars that migrated there, since
that is the radius at which we shut off star formation. The gas gradient shows
a break in the scale radius near $\rgal\approx$ 6 kpc. This is a consequence
of our adopted star formation law; at $\Sigma_{\text{g}}=2\times 10^{7}$ the
relation changes from linear at higher densities to
$\dot{\Sigma}_{\star}\sim\Sigma_{\text{g}}^{3.6}$ at lower densities (see
discussion in § 1.6).

Although we adopt the Licquia2015 disc stellar mass of the Milky Way here, we
find similar results when taking a value which differs even by an order of
magnitude. If we used a linear star formation law, then our chemical evolution
would be independent of the mass normalization as it is in corresponding one-
zone models (Spitoni2017; Weinberg2017; Belfiore2019). The mass normalization
enters our calculation because it affects the transition between the linear
and non-linear regimes of our star formation law (equation 2n), but the impact
on abundance evolution is minimal.

###  1.8 Summary

In summary, our fiducial model has an inside-out SFH with e-folding timescales
derived from the observations of Sanchez2020. Radial migration proceeds in a
manner in which our model stellar populations have a change in radius
$\Delta\rgal$ informed from the h277 hydrodynamical simulation
(Christensen2012; Zolotov2012; Loebman2012; Loebman2014; Brooks2014, see
discussion in § 1.1). In the baseline model, stars move to their final radius
with a $\sqrt{\text{age}}$ dependence (Frankel2018; Frankel2020, see
discussion in 1.2). Using to calculate abundances for O and Fe in this paper,
our supernova yields are adopted from Johnson2020, who in turn take these
values from Weinberg2017 (see § 1.3). Outflows are characterized such that the
equilibrium abundance of oxygen under a constant SFH follows an abundance
gradient in agreement with observational results in the Milky Way (see § 1.4).
Our star formation law is based on the Bigiel2010 and Leroy2013 data presented
in comparison with theoretical models in Krumholz2018a. To describe the
stellar surface density gradient, we adopt the two-exponential form describing
the thin and thick discs from Bland-Hawthorn2016. We adopt the Kroupa2001 IMF
throughout this paper.

Our selection of star particles from yields a sample of 1,751,765 candidate
analogues with disc-like kinematics at the present day (see discussion in §
1.1). We take $\Delta\rgal$ = 100 pc as the width of each annulus from = 0 to
20 kpc and a timestep size of $\Delta T$ = 10 Myr from $T$ = 0 to 13.2 Gyr.
With the resulting 200 zones and 1,321 timesteps (one extra so that age = 0
stars are included), we let form $n$ = 8 stellar populations per zone per
timestep, resulting in 2,113,600 total stellar populations with predicted
masses and abundances. We set the SFR to zero beyond = 15.5 kpc; stellar
populations do form beyond this radius and are part of the computational
overhead, but they have zero mass and thus do not contribute to the chemical
evolution in our models. This results in 1,627,472 stellar populations with
non-zero masses and abundances, comparable to the total number of disc
particles in our sample from . These simulations run in $\sim$5 hours on a
single core with a 3 GHz processor and take up $\sim$235 MB of disc space per
output, including the extra data that we record for each stellar population’s
analogue star particle. We have also run variations with $n$ = 2 stellar
populations per zone per timestep, finding similar results in all cases.

We have run for all four of our SFHs, all four migration models, and all four
variations in $\tau_{\text{mol}}$ noted in § 1.6 — a total of 64 simulations,
as well as a variety of other test cases. For many of our results, only the
SFH variations have a substantial impact. We discuss the impact of migration
or $\tau_{\text{mol}}$ variations where they are relevant.


# Master Stability Functions of Networks of Izhikevich Neurons

Raul P. Aristides<EMAIL_ADDRESS>São Paulo State University (UNESP),
Instituto de Física Teórica, Rua Dr. Bento Teobaldo Ferraz 271, Bloco II,
Barra Funda, 01140-070 São Paulo, Brazil.  Hilda A. Cerdeira  São Paulo State
University (UNESP), Instituto de Física Teórica, Rua Dr. Bento Teobaldo Ferraz
271, Bloco II, Barra Funda, 01140-070 São Paulo, Brazil.  Epistemic, Gómez
$\&$ Gómez Ltda. ME, Rua Paulo Franco 520, Vila Leopoldina, 05305-031 São
Paulo, Brazil

###### Abstract

Synchronization has attracted the interest of many areas where the systems
under study can be described by complex networks. Among such areas is
neuroscience, where is hypothesized that synchronization plays a role in many
functions and dysfunctions of the brain. We study the linear stability of
synchronized states in networks of Izhikevich neurons using Master Stability
Functions, and to accomplish that, we exploit the formalism of saltation
matrices. Such a tool allows us to calculate the Lyapunov exponents of the
Master Stability Function (MSF) properly since the Izhikevich model displays a
discontinuity within its spikes. We consider both electrical and chemical
couplings, as well as total and partially synchronized states. The MSFs
calculations are compared with a measure of the synchronization error for
simulated networks. We give special attention to the case of electric and
chemical coupling, where a riddled basin of attraction makes the synchronized
solution more sensitive to perturbations.

##  I Introduction

Neurons are the building blocks of the brain. These structures can display a
rich variety of dynamics depending on their type and the inputs that it
receives from other neurons [1]. Although many features of neurons remain
unclear, models are capable of reproducing their observed dynamics. The
seminal work of Hodgkin-Huxley (HH) introduced a neuron model based on the
giant axon of the squid [2]. Being biophysically accurate, this model become a
standard for ODE-based neuron models. Naturally, the collective dynamics of
neurons are also of interest, and unfortunately, networks of complex models
like the HH can be difficult to deal with both analytically and
computationally. With this in mind, E. M. Izhikevich proposed a neuron model
that combines biological features of complex models, like the HH, with
computational efficiency [3].

Due to these characteristics, the Izhikevich model has been implemented in the
study of networks of spiking neurons, including large-scale models [4, 5].
Such models can be used to deepen our understanding of how the interplay
between synaptic and neuronal processes produces collective behaviors. Of
great interest is the emergence of rhythms and synchronization of neural
activity, as it suggests that the high-dimensional dynamics of neuronal
networks can collapse into low dimensional oscillatory modes [6]. Although
rhythmic activity and synchronization are often associated with brain
pathologies like Parkinson’s disease and epilepsy [7, 8], evidence suggests
that these two mechanics are involved in information processing, memory
formation, and cognition [9, 10, 11].

Since the mathematical description of synchronization is well established in
nonlinear dynamics, it can be applied to many natural phenomena, including
neural activity. A central question regarding synchronization theory is
whether a synchronized state is stable or not. Among the many contributions to
this topic over the last decades, here we focus on the Master Stability
Function (MSF) formalism [12]. The main result of such formalism is a
diagonalized variational equation, which allows us to calculate the Lyapunov
exponents associated with perturbations transverse to the synchronized
manifold. In the original work, it was made possible by assuming that the
system can be locally linearized around the synchronized solution, as well as
a diagonalizable coupling matrix. Fortunately, this formalism has been
extended to a variety of cases including near-identical systems, non-
diagonalizable, and multilayer networks [13, 14, 15]. Moreover, the stability
of partial synchronization, which includes cluster synchronized and chimera
states have also been studied in the context of master stability functions
[16, 17, 18, 19].

In this paper, we apply the MSF formalism to study the linear stability of
synchronized states in networks of Izhikevich neurons. Both global and partial
synchronization patterns are considered. To calculate the associated Lyapunov
exponents, we use the Saltation Matrices formalism, which takes into account
the discontinuities of the Izhikevich model. First, we consider the case of
global synchronization in networks with electrical coupling, followed by the
case with chemical coupling and then with both coupling schemes. For each
case, we study how the stability of the synchronized state depends on the
coupling strength. We verify that the MSF formalism is effective for probing
the threshold for synchronization, except for the last case, where the outcome
of the system depends on its initial conditions. This can be the case when the
system under study has a riddled basin of attraction [20]. As a result, the
system evolves towards a non synchronized solution even with initial
conditions in the neighborhood of the synchronization state [21].

##  II MSF applied to Izhikevich model

The Izhikevich model is a two-dimensional neuronal model, celebrated for its
biophysical accuracy, and its capacity of displaying many spiking patterns
without losing computational efficiency. Letting $\mathbf{x}=[x,y]^{T}$, the
model consists of a system of differential equations [3]

$\mathbf{\dot{x}}=F(\mathbf{x})=\begin{bmatrix}0.04x^{2}+5x+140-y+I\\\
a(bx-y)\end{bmatrix}$ (1)

with the after-spike reset given by

$\text{if}\quad x>30\text{mV},\quad\text{then}\begin{cases}x\rightarrow c\\\
y\rightarrow y+d\end{cases}$ (2)

The variables $x$ and $y$ represent the membrane potential and a membrane
recovery variable respectively. The parameter $I$ takes into account synaptic
currents or injected dc-currents. The role of the remaining parameters $a$,
$b$, $c$ and $d$ are discussed in [3]. For all calculations we used $a=0.2$,
$b=2$, $c=-56$, $d=-16$ and $I=-99$.

The Izhikevich model can be seen as a hybrid model, a continuous-time
evolution process interfaced with a logical process, in this case, the
resetting function when $x$ crosses the threshold. These discontinuities
interfere with the usual straightforward methods [24, 25, 26] to compute the
Lyapunov exponents of such a system.

Following Bizarri et al. [27] we resort to saltation matrices, which allow us
to calculate the Lyapunov exponents associated with the master stability
functions of networks of Izhikevich neurons. In this framework, these matrices
are used as correction factors at the instants of the discontinuities due to
the resetting function, yielding the correct Lyapunov exponents.

Between spikes the dynamics of the Izhikevich is smooth, and for a small
perturbation $\delta\mathbf{x}$ we have

$\delta\dot{\mathbf{x}}={DF\big{|}_{\mathbf{x}}}\delta\mathbf{x}\ ,$ (3)

where $DF\big{|}_{\mathbf{x}}$ is the Jacobian matrix applied to $F$ evaluated
at $\mathbf{x}$. If a spike occurs, the after-spike reset (2) induces a
discontinuity in the flow such that $DF\big{|}_{\mathbf{x}^{-}}\neq
DF\big{|}_{\mathbf{x}^{+}}$, where $-$ and $+$ indicates the before and after
the reset. Then, the perturbation $\delta\mathbf{x}^{+}$ after the reset event
is given by

$\delta\mathbf{x}^{+}=S({DF\big{|}_{\mathbf{x^{-}}}}\delta\mathbf{x}^{-})\ ,$
(4)

which allows us to write

$\dot{\delta\mathbf{x}^{+}}={DF\big{|}_{\mathbf{x}^{+}}}\delta\mathbf{x}^{+}\
,$ (5)

where $S$ is the saltation matrix, which for Izhikevich neurons is given by

$S=\begin{bmatrix}\dfrac{\dot{x}^{+}}{\dot{x}^{-}}&0\\\
\dfrac{\dot{y}^{+}-\dot{y}^{-}}{\dot{x}^{-}}&1\end{bmatrix}\ .$ (6)

The complete derivation of the procedure can be found at [28, 27].

###  II.1 Electrical coupling

First, we study the effect of electrical coupling, which represents electrical
synapses between nodes [29]. Consider a network with $N$ nodes, where the
isolated dynamics of each node is given by Eq. (1). An adjacency matrix $A$
indicates whether node $i$ is connected diffusively to node $j$. In this case,
the dynamics of each node is modified by

$W_{i}=-g_{e}\sum_{j=1}^{N}A_{ij}(x^{j}-x^{i})$ (7)

where $g_{e}$ is the coupling conductance (strength) and $A_{ij}$ is an
element of the adjacency matrix. Note that in this case, the coupling can be
written as a function of the Laplacian matrix of the network, given by
$L={D}-A$ where $D$ is the degree matrix defined by $D_{ii}=\sum_{j}A_{ij}$,
then $W_{i}=g_{e}\sum_{j=1}^{N}L_{ij}x^{j}$. If we define
$\mathbf{x}=[\mathbf{x}^{1},...,\mathbf{x}^{N}]^{T}$, where
$\mathbf{x}^{i}=[x^{i},y^{i}]^{T}$, and
$\mathbf{F}(\mathbf{x})=[\mathbf{F}(\mathbf{x}^{1}),...,\mathbf{F}(\mathbf{x}^{N})]^{T}$,
we can write the dynamical system of the network as

$\mathbf{\dot{x}}=\mathbf{F}(\mathbf{x})-g_{e}(L\otimes\mathbf{G})\mathbf{x}$
(8)

where $\otimes$ is the Kronecker product and

$\mathbf{G}=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix}$ (9)

is the matrix encoding the coupling scheme. The synchronization manifold
$\mathbf{x}^{s}$ is defined by the $N-1$ constraints
$\mathbf{x}^{1}=\mathbf{x}^{2}=...=\mathbf{x}^{N}$, and since the Laplacian
matrix is a zero row-sum matrix, $(L\otimes\mathbf{G})\mathbf{x^{s}}$ = 0, we
have $\mathbf{\dot{x}^{s}}=\mathbf{F}(\mathbf{x}^{s})$. The linear stability
of the synchronized solution can be investigated with the MSF of Eq. (8)

$\dot{\eta}_{k}=[{D\mathbf{F}\big{|}_{\mathbf{x}^{s}}}-g_{e}\gamma_{k}{D\mathbf{G}\big{|}_{\mathbf{x}^{s}}}]\eta_{k}\
,$ (10)

where $\gamma_{k}$ is the $k$-th eigenvalue of $L$. Thus, the analysis is
decomposed in eigenmodes, one being parallel to the synchronization manifold
$\gamma_{0}$, the others transverse to it. We are interested in the Lyapunov
exponents of transverse modes, which tell us whether or not perturbations to
the synchronized solution will decay. Equation (10) was first derived by
Pecora and Carroll [12] for undirected networks ($A=A^{T}$) but since then it
has been extended to a broad range of networks and synchronization patterns
[13, 14, 15, 16, 17, 18, 19]. We highlight that at the event of a spike, the
synchronized solution $\mathbf{x}^{s}$ encounters a discontinuity, and
therefore at this point, we need to apply the saltation matrix Eq. (6).

Figure 1: Maximum lyapunov exponent for Izhikevich model with electrical
coupling.

The maximum Lyapunov exponent (MLE) of Eq. (10) is presented in Fig. 1. The
expanded inset shows that $\Lambda_{max}$ becomes negative for
$g_{e}\gamma>0.2670$. To check our results, we simulate a network of $N=2$
nodes, as depicted in Fig. 2. The Laplacian matrix in such case is

Figure 2: Network with $N=2$ used to verify the result of Fig. 1.

$L=\begin{bmatrix}1&-1\\\ -1&1\end{bmatrix}$ (11)

with eigenvalues $0$ and $2$, and $\gamma_{2}=2$ is associated with the
eigenspace transverse to the synchronized manifold. Putting together that the
MLE is negative for $g_{e}\gamma>0.2670$, and that $2$ is the eigenvalue
associated with the transverse mode, we conclude that for $g_{e}>0.133$ the
synchronized solution is stable.

Fig. 3 depicts the synchronization error, which is computed as
$\langle|\mathbf{x}^{1}-\mathbf{x}^{2}|\rangle$, in good agreement with the
result given by the MSF.

Figure 3: Synchronization error for N=2 with electrical coupling. It goes to
zero around $\sigma\approx 0.133$, in conformity with the MSF.

###  II.2 Chemical coupling

We now study the synchronization of an Izhikevich network with chemical
synaptic coupling between nodes. The chemical synapses that node $i$ receives
from every other node $j$ are modeled by the sigmodal function $B_{i}$ [30]

$B_{i}=g_{c}(x^{i}-v_{s})\sum_{j=1}^{N}A_{ij}\zeta(x^{j}),$ (12)

where $v_{s}$ is the reversal potential,
$\zeta(x)=[1+\exp(-\epsilon(x-\theta))]^{-1}$, $\epsilon$ defines the slope of
the sigmoidal function and $\theta$ is the synaptic firing threshold.
Througout this study we set $v_{s}=20$, $\epsilon=7$ and $\theta=0$. This set
of parameters allows the chemical synapses to be both excitatory and
inhibitory, since $v_{s}$ lies inside the range of oscillation of the action
potential $x_{i}$. This is the case for some classes of neurons, as discussed
in [31, 32].

Note that this coupling term cannot be written as a function of the Laplacian
matrix, instead we can write it in terms of the adjacency matrix. To obtain an
equation analogous to (8), we introduce a few objects, which are defined in
the Appendix A, together with the complete derivation of the equation for the
dynamics of the network with chemical coupling, which is given by:

$\mathbf{\dot{x}}=\mathbf{F}(\mathbf{x})-g_{c}\mathbf{T}(\mathbf{x})[(A\otimes\mathbf{C})(\mathbf{x})]\
.$ (13)

The synchronized solution, $\mathbf{x}^{1}=\mathbf{x}^{2}=...=\mathbf{x}^{N}$,
derived in the Appendix A from Eq. (53) to (57), takes the form

$\mathbf{\dot{x}}^{s}=\mathbf{F}(\mathbf{x}^{s})-g_{c}k_{n}\mathbf{T}(\mathbf{x}^{s})\mathbf{C}(\mathbf{x}^{s})\
,$ (14)

and it exists only if the number of links $k_{n}$ is the same for all nodes
[33]. Following the Pecora-Carroll analysis, we obtain the Master Stability
Function for this case, and after we diagonalize $A$ we have

$\displaystyle\dot{\eta}=\\{[\mathbf{I}_{N}\otimes({D\mathbf{F}\big{|}_{\mathbf{x}^{s}}}-g_{c}k_{n}{D\mathbf{H}\big{|}_{\mathbf{x}^{s}}}\mathbf{K}(\mathbf{x}^{s}))]-g_{c}\Gamma\otimes\mathbf{T}(\mathbf{x}^{s}){D\mathbf{C}\big{|}_{\mathbf{x}^{s}}}\\}\eta\
.$ (15)

The complete derivation of this equation, which is analogous to Eq. (10), is
given in Appendix A and it was first derived by Checco et al. [33]. As an
example, we consider the network of $N=4$ nodes, depicted in Fig. 4, with an
adjacency matrix given by

$A=\begin{bmatrix}0&1&0&1\\\ 1&0&1&0\\\ 0&1&0&1\\\ 1&0&1&0\end{bmatrix}$ (16)

Figure 4: This network, with $N=4$, does not synchronize when the nodes
present only chemical coupling, this situation changes when the electrical
coupling is added.

The adjacency matrix Eq. (16) has an eigenvalue $\gamma_{1}=2$ associated with
the synchronized solution, so we need to calculate the MLE for the remaining
eigenvalues $\\{0,0,-2\\}$. Figure 5 shows the MLE of Eq. (65), which is
positive for the range studied. Our calculation for the synchronization error
in Fig. (6) shows good agreement with the MLE. The sudden increase in the
magnitude of the synchronization error at $g_{c}\approx 0.48$ will be
discussed in the context of partial synchronization on Section III.

Figure 5: Maximum Lyapunov exponent for Izhikevich model with chemical
coupling. Figure 6: Synchronization error for Izhikevich model with chemical
coupling.

###  II.3 Chemical and electrical coupling

Finally, we study the effect of both electrical and chemical coupling on
synchronization. In this case, the evolution of the network is given by

$\mathbf{\dot{x}}=\mathbf{F}(\mathbf{x})-g_{e}L\otimes\mathbf{G}(\mathbf{x}))-g_{c}\mathbf{T}(\mathbf{x})(A\otimes\mathbf{C})\mathbf{x}\
.$ (17)

For the sake of simplicity, we consider that both interactions depend on the
same diagonalizable adjacency matrix, i.e., $L=D-A$, and that every node has
the same number of neighbors $k_{n}$. Combining the results of previous
sections, we write the variational equation for the perturbations to the
synchronized state as

$\displaystyle\delta\mathbf{\dot{x}}=\\{[\mathbf{I}_{N}\otimes({D\mathbf{F}\big{|}_{\mathbf{x}^{s}}}-g_{c}k_{n}{D\mathbf{H}\big{|}_{\mathbf{x}^{s}}}\mathbf{K}(\mathbf{x}^{s}))]-g_{e}L\otimes{D\mathbf{G}\big{|}_{\mathbf{x}^{s}}}-g_{c}A\otimes\mathbf{T}(\mathbf{x}^{s}){D\mathbf{C}\big{|}_{\mathbf{x}^{s}}}\\}\delta\mathbf{x}\
.$ (18)

Since $L$ and $A$ commute, we can diagonalize both simultaneously, yielding

$\displaystyle\dot{\eta}=\\{[\mathbf{I}_{N}\otimes({D\mathbf{F}\big{|}_{\mathbf{x}^{s}}}-g_{c}k_{n}{D\mathbf{H}\big{|}_{\mathbf{x}^{s}}}\mathbf{K}(\mathbf{x}^{s}))]-[g_{e}\Gamma^{L}\otimes{D\mathbf{G}\big{|}_{\mathbf{x}^{s}}}-g_{c}\Gamma^{A}\otimes\mathbf{T}(\mathbf{x}^{s}){D\mathbf{C}\big{|}_{\mathbf{x}^{s}}}]\\}\eta\
.$ (19)

where $\Gamma^{L}$ and $\Gamma^{A}$ are diagonal matrices with eigenvalues of
$L$ and $A$ respectively as entries.

Once again, we take a network of $N=4$ neurons, with electrical and chemical
coupling given by $g=g_{e}=g_{c}$ and adjacency matrix (16). Thus, the
eigenvalues associated with transverse modes are
$\\{\gamma^{A}\\}=\\{0,0,\\-2\\}$ from the adjacency matrix and
$\\{\gamma^{L}\\}=\\{2,2,4\\}$, from the Laplacian matrix. We calculate the
MLE of Eq. (19) as a function of $g$, and as shown in Fig. (7), the global
synchronization pattern becomes stable for $g\approx 0.13$. To confirm this
result we calculate the synchronization error over $200$ simulations, which is
shown in Fig. 8 and find that it does not go to zero when the MLE becomes
negative.

We notice that the synchronization error only vanishes around $g\approx 0.2$,
and for values below that, the system can either synchronize or not. Moreover,
for $g\approx 0.23$, we see that some simulations yielded a synchronization
error greater than expected. To better understand these results, we calculate
the synchronization error for different initial conditions, for $g=0.18$ and
$g=0.23$. As we see in Fig. 9 (a), for $g=0.18$, depending on the initial
conditions, the system can end up in different solutions, a characteristic of
multistable systems [34]. Moreover, the different solutions are intertwined in
a complex way. This type of basin of attraction is often called a riddled
basin, and can be the cause for the divergence between the MLE and the
synchronization error results [21].

Figure 7: Maximum Lyapunov exponent for Izhikevich model with electrical and
chemical coupling. Figure 8: Synchronization error for N=4 with electrical and
chemical coupling. The grey lines represent maximum and minimum errors and the
circles the average over $200$ simulations.

For $g=0.23$, although the basin of attraction is dominated by synchronized
solutions (which represent $96.33\%$ the cases), it is studded with non-
synchronized ones, as shown in Fig. 9 (b). Among these cases, few ones yield a
synchronization error greater than $15$ ($0.095\%$ of the cases). We verified
that in such cases, the network switches from the synchronized solution to a
cluster synchronized state with two clusters $(1,3)$ and $(2,4)$ in such a way
that these clusters are synchronized but out of phase, resulting in a greater
synchronization error. It is worth mentioning that the dynamics of the nodes
also changed, as we will discuss later in the text. These results lead us to
conclude that the electrical coupling induces synchronization in our network.

Besides that, it is usually assumed that the presence of riddled basins of
attraction is related to the MLE approaching zero from below, as discussed in
[22, 23]. A similar result, where the system fails to exhibit synchronization
when the MLE approaches zeros from above was reported by [21] but the presence
of the riddled basin was not investigated.

Figure 9: Riddled basin of attraction in $(x_{1},x_{3})$ plane for different
values of $g$, calculated with $10.24\times 10^{4}$ initial conditions in
$x_{1},x_{3}\in[1,-1]$. (a) is calculated with $g=0.18$. (b) $g=0.23$

##  III Izhikevich model: Partial Synchronization

###  III.1 Electrical coupling

As mentioned in the first section, the electrical coupling can be written in
terms of the Laplacian matrix. In this case, the variational equation
concerning the linear stability around the cluster synchronized state (CSS) is
given by [16]

$\delta\dot{\mathbf{x}}=[\sum_{m=1}^{M}(E^{m}\otimes{D\mathbf{F}\big{|}_{\mathbf{s}^{m}}}-g_{e}(LE^{m}\otimes{D\mathbf{G}\big{|}_{\mathbf{s}^{m}}}))]\delta\mathbf{x}$
(20)

where $M$ is the number of clusters and $E^{m}$ is an $N\times N$ diagonal
indicator matrix for each cluster such that $E^{m}_{ii}=1$ if $i$ belongs to
cluster $m$, and it is equal to zero otherwise. If we consider a CSS with $M$
clusters, the dynamical evolution of such state is given by [17, 18]

$\displaystyle\dot{\mathbf{x}}$
$\displaystyle=(Z\otimes\mathbf{I}_{q})\dot{\mathbf{s}},\quad\text{where}$
(21) $\displaystyle\dot{\mathbf{s}}$
$\displaystyle=\mathbf{F}(\mathbf{s})-g_{e}(L^{\pi}\otimes\mathbf{I}_{q})\mathbf{G}(\mathbf{s})$
(22)

and $Z$ is the $N\times M$ indicator matrix which encodes the cluster state
configuration $Z:Z_{ij}=1$ if node $i$ is part of cluster $C_{j}$ and
$Z_{ij}=0$ otherwise. $\mathbf{I}_{q}$ is an $q\times q$ identity matrix,
where $q$ is the dimension of the isolated dynamical system, $q=2$ in our
case. The state of each cluster is encoded in the $M\times q$-dimensional
vector $\mathbf{s}$. If the quotient network Laplacian $L^{\pi}$ satisfies the
condition [35]

$LZ=ZL^{\pi}\ ,$ (23)

where $L^{\pi}=(Z^{T})^{-1}Z^{T}LZ=Z^{+}LZ$ and $Z^{+}$ is the left Moore-
Penrose of pseudoinverse $Z$ [36, 37], then the partition encoded by $Z$ is
said to be an Equitable External Partition (EEP) [35]. In this case, the
network is divided into $M$ clusters in such way that the number of
connections from nodes in a cluster $C_{i}$ to nodes in a cluster $C_{j}$
depends only on $i,j$ with $i\neq j$. The quotient network dynamics is a
coarse-grained version of the original network, where each cluster becomes a
node and the weights between these nodes are the out-degrees between clusters
in the original network [38, 18].

Again, we start our example with a network of Izhikevich neurons coupled
through electrical synapses. Consider a network with the following Laplacian
matrix

$L=\begin{bmatrix}2&\\-1&0&\\-1&0&0&0&0\\\ \\-1&2&0&\\-1&0&0&0&0\\\
0&0&1&\\-1&0&0&0&0\\\ \\-1&\\-1&\\-1&5&\\-1&\\-1&0&0\\\
0&0&0&\\-1&3&0&\\-1&\\-1\\\ 0&0&0&\\-1&0&3&\\-1&\\-1\\\
0&0&0&0&\\-1&\\-1&3&\\-1\\\ 0&0&0&0&\\-1&\\-1&\\-1&3\end{bmatrix}\ $ (24)

and a CSS with indicator matrix $Z$ given by

$Z=\begin{bmatrix}1&0&0&0&0\\\ 1&0&0&0&0\\\ 0&1&0&0&0\\\ 0&0&1&0&0\\\
0&0&0&1&0\\\ 0&0&0&1&0\\\ 0&0&0&0&1\\\ 0&0&0&0&1\\\ \end{bmatrix}\ .$ (25)

One can verify that this choice of $L$ and $Z$ satisfies (23) and therefore
this cluster configuration is an EEP. A representation of such network is
given in Fig. 10, where nodes of the same color represent a cluster. Putting
this together, Eq. (22) reads

$[\dot{\mathbf{x}}]^{T}=[\dot{\mathbf{s}}^{1},\dot{\mathbf{s}}^{1},\dot{\mathbf{s}}^{3},\dot{\mathbf{s}}^{4},\dot{\mathbf{s}}^{5},\dot{\mathbf{s}}^{5},\dot{\mathbf{s}}^{7},\dot{\mathbf{s}}^{7}]^{T}\
.$ (26)

Figure 10: Network studied in this section, where the clusters are indicated
with colors.

After block-diagonalizing the variational equation (20), we investigated the
linear stability of the CSS with the following equation

$\mathbf{\dot{\eta}}=\Big{[}\sum_{m}Q^{m}\otimes{D\mathbf{F}\big{|}_{\mathbf{s}^{m}}}-g_{e}\Gamma^{L}Q^{m}\otimes{D\mathbf{G}\big{|}_{\mathbf{s}^{m}}}\Big{]}\mathbf{\eta}\
.$ (27)

Where $Q^{m}=V^{-1}E^{m}V$ is a block diagonal matrix and
$\Gamma^{L}=V^{-1}LV$ is a diagonal matrix filled with eigenvalues of $L$. The
orthogonal matrix that block-diagonalizes both ${E^{m}}$ and $L$ can be found
using group theory [17], External Equitable Partitions (EEP’s) [18] and
Simultaneous Block Diagonalization (SBD) [19] approaches, resulting in the
decoupling between modes within the cluster synchronized manifold and
transverse to it.

To check if all transverse modes are damped, we calculated the maximum
Lyapunov exponent of Eq. (27), which is given in Fig. 11 as a function of
$g_{e}$.

Figure 11: Maximum Lyapunov exponent for Izhikevich model with electrical
coupling under partial synchronization.

Fig. 12 shows the synchronization error of the cluster synchronization state
as a function of $g_{e}$. The error is computed as
$\langle|x^{i}(t)-x_{m}^{i}(t)|\rangle$, where we take an average both in time
and over the neurons $i$ in a cluster $m$. $x_{m}^{i}$ is the average state
for the neurons in the cluster to which node $i$ belongs.

Fig. 11 and Fig. 12 show excellent agreement between the stability analysis
and the synchronization error: the MLE associated with the transverse modes
goes negative around $g_{e}\approx 0.10$, and the synchronization error
vanishes around the same value.

Figure 12: Synchronization error for Izhikevich model with electrical
coupling, partial synchronization.

###  III.2 Chemical coupling

We now extend the linear stability study of cluster synchronized states to
networks with chemical coupling, which is mediated by the adjacency matrix of
the network. Consider a network with $N$ nodes and adjacency matrix $A$, given
a CSS with $M$ clusters encoded by an indicator matrix $Z$, the coarse-grained
dynamics can be written as follows

$\displaystyle\dot{\mathbf{x}}$
$\displaystyle=(Z\otimes\mathbf{I}_{q})\dot{\mathbf{s}},\quad\text{where}$
(28) $\displaystyle\dot{\mathbf{s}}$
$\displaystyle=\mathbf{F}(\mathbf{s})+g_{e}\mathbf{T}(\mathbf{s})(A^{\pi}\otimes\mathbf{C})(\mathbf{s})\
.$ (29)

If the adjacency matrix of the quotient network is given by $A^{\pi}=Z^{+}AZ$,
its diagonal entries corresponds to self-loops in the quotient network
whenever we have connections between nodes in the same cluster, and the
partition enconded by $Z$ is an External Partition (EP). For EPs the number of
connections from nodes in a cluster $C_{i}$ to nodes in a cluster $C_{j}$
depends only in $i,j$ [18].

The variational equation of Eq.(28) is

${\delta\dot{\mathbf{x}}}=\big{\\{}\sum_{m=1}^{M}[E^{m}\otimes({D\mathbf{F}\big{|}_{\mathbf{s}^{m}}}-g_{c}PR^{m})]-[\sum_{m=1}^{M}(E^{m}A\otimes\mathbf{T}(\mathbf{s}^{m}))(\sum_{m=1}^{M}E^{m}\otimes{D\mathbf{K}\big{|}_{\mathbf{s}^{m}}})]\big{\\}}\delta\mathbf{x}\
.$ (30)

And after block-diagonalizing Eq. (30), we obtain

$\dot{\eta}=\big{\\{}\sum_{m=1}^{M}Q^{m}\otimes[{D\mathbf{F}\big{|}_{\mathbf{s}^{m}}}-g_{c}\mathbf{G}R^{m}]-(Q^{m}\Gamma^{A}\otimes\mathbf{T}(\mathbf{s}^{m})(\sum_{m=1}^{M}Q^{m}\otimes{D\mathbf{K}\big{|}_{\mathbf{s}^{m}}})\big{\\}}\eta\
.$ (31)

where we defined

$\mathbf{G}=\begin{bmatrix}1&0\\\ 0&0\end{bmatrix}\quad\text{and}\quad
R^{m}=\sum_{n=1}^{M}A^{\pi}_{mn}\zeta(x_{n})\ .$ (32)

Once again, we consider the network represented in Fig. 4, but now we look for
the stability of a CSS, which can explain the sudden increase in the
synchronization error around $g_{c}\approx 0.48$ in Fig. 6. The CSS is encoded
by the indicator and quotient adjacency matrices, which are given by

$Z=\begin{bmatrix}1&0\\\ 0&1\\\ 1&0\\\ 0&1\\\
\end{bmatrix}\quad\text{and}\quad A^{\pi}=\begin{bmatrix}0&2\\\ 2&0\\\
\end{bmatrix}\ .$ (33)

The linear stability analysis of such CSS is given in Fig. 13, where we notice
that the CSS becomes stable for $g_{c}>0.53$, which is around the value for
which the synchronization error suddenly increases in Fig. 6. This can be
explained because the network dynamics are similar to one composed of $N=2$
neurons out-of-phase. Figure 14 shows the synchronization error of the CSS in
question, which confirms our stability analysis.

Figure 13: Maximum Lyapunov exponent for Izhikevich model with chemical
coupling partial synchronization (N=4).  Figure 14: Synchronization error for
Izhikevich model with chemical coupling, partial synchronization (N=4).

We now consider the network of the previous subsection, depicted in Fig. 10,
and the same indicator matrix $Z$ (Eq. (25)), with chemical coupling between
neurons. In this case, the adjacency matrix is given by

$A=\begin{bmatrix}0&1&0&1&0&0&0&0\\\ 1&0&0&1&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
1&1&1&0&1&1&0&0\\\ 0&0&0&1&0&0&1&1\\\ 0&0&0&1&0&0&1&1\\\ 0&0&0&0&1&1&0&1\\\
0&0&0&0&1&1&1&0\\\ \end{bmatrix}\ .$ (34)

Unlike the case with electrical coupling, here the CSS is not stable, as seen
in Fig. (15), where the maximum Lyapunov exponent is positive in the range of
$g_{c}$ considered. This analysis is supported by the synchronization error of
this case (Fig. 16), which is greater than $0$ in the range of $g_{c}$
studied.

Figure 15: Maximum Lyapunov exponent for Izhikevich model with chemical
coupling partial synchronization (N=8).  Figure 16: Synchronization error for
Izhikevich model with chemical coupling, partial synchronization (N=8).

###  III.3 Electrical and chemical coupling

In the case where the network admits both eletrical and chemical coupling, the
dynamics of the CSS is governed by

$\displaystyle\dot{\mathbf{x}}=$
$\displaystyle(Z\otimes\mathbf{I}_{q})\dot{\mathbf{s}},\quad\text{where}$ (35)
$\displaystyle\dot{\mathbf{s}}=\mathbf{F}(\mathbf{s})-g_{c}\mathbf{T}(\mathbf{s})$
$\displaystyle(A^{\pi}\otimes\mathbf{C})(\mathbf{s})-g_{e}(L^{\pi}\otimes\mathbf{I}_{q})\mathbf{G}(\mathbf{s})\
.$ (36)

where $Z$ is the indicator matrix of the CSS state, $A^{\pi}$ and $L^{\pi}$
are the quotient adjacency and Laplacian matrices, and for the sake of
simplicity we take $g_{c}=g_{e}=g$. Putting together the results of the later
subsections, the master stability function of a given CSS state can be written
as

$\dot{\eta}=\big{\\{}\sum_{m=1}^{M}Q^{m}\otimes[{D\mathbf{F}\big{|}_{\mathbf{x}^{m}}}-g_{c}\mathbf{G}R^{m}]-g_{c}(\sum_{m=1}^{M}Q^{m}\Gamma^{A}\otimes\mathbf{T}(\mathbf{x}^{m})(\sum_{m=1}^{M}Q^{m}\otimes{D\mathbf{K}\big{|}_{\mathbf{x}^{m}}})-g_{e}\sum_{m=1}^{M}\Gamma^{L}Q^{m}\otimes{D\mathbf{H}\big{|}_{\mathbf{x}^{m}}}\big{\\}}\eta\
.$ (37)

We proceed with the analysis of the network in Fig. 10, with indicator matrix
(25). The MLE of Eq. (37) is shown in Fig. 17, where we notice that the
introduction of electrical coupling in the network seems to induce the CSS,
since the MLE becomes negative for $g>0.172$.

Figure 17: Maximum Lyapunov exponent for the Izhikevich model with chemical
and electric coupling, partial synchronization (N=8).

The initial conditions, used to make Fig. 18 (nearly identical but random),
were equivalent to adding perturbations to all three directions transverse to
the synchronized solution, plus the five directions parallel to the
synchronized solution. The dotted line represents the average over $100$
simulations and the shaded area the error bars of these simulations. Although
the lower bound of the synchronization error goes to zero around $g\approx
0.172$, we see that on average the system does not reach the synchronized
solution. This sensibility to perturbations is characteristic of multistable
systems [34], and as discussed in [21], a riddled basin of attraction can be
the cause for such behavior, which in turn may be a consequence of the
discontinuity of the differential equation. A thorough understanding of this
phenomenae needs an intense research which is outside the scope of this work.

Figure 18: Synchronization error as a function of $g$ for Izhikevich model
with electrical and chemical coupling. The CSS considered is the one Fig. 10.
The average over $150$ runs, with nearly identical initial conditions, is
depicted with its bounds.

To verify that, we consider adding a disturbance in only one direction
transverse to the synchronized solution, which means perturbing only one
cluster. For example, if we take the solution

$\mathbf{x}=\begin{bmatrix}\mathbf{s}_{1}\\\ \mathbf{s}_{1}\\\
\mathbf{s}_{3}\\\ \mathbf{s}_{4}\\\ \mathbf{s}_{5}\\\ \mathbf{s}_{5}\\\
\mathbf{s}_{7}\\\ \mathbf{s}_{7}\end{bmatrix}\ ,\quad\text{and
add}\quad\mathbf{p}_{\perp}=\alpha\begin{bmatrix}1\\\ \\-1\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\ 0\end{bmatrix}\ \otimes\begin{bmatrix}1\\\ 1\end{bmatrix},$ (38)

where $\alpha$ is a constant used to control the magnitude of the perturbation
$\mathbf{p}_{\perp}$, we are pulling the first cluster out of synchronicity.
Fig. 19 (a) shows how the synchronization error, represented by colors, is
sensible to initial conditions, for $g=0.18$. In this case, we consider
initial conditions slightly perturbed as in Eq. (38). The initial conditions
for $\mathbf{x}_{(1,2)}=[x^{(1,2)},y^{(1,2)}]^{T}$ vary from $[1,1]$ to
$[\\-1,\\-1]$, which is equivalent to $\alpha$ varying from $1$ to $-1$ in Eq.
(38). Along the diagonal $\mathbf{x}_{1}=\mathbf{x}_{2}$ and the system is
always synchronized. In the off-diagonal region, we have both synchronized and
non-synchronized solutions. The basin is riddled with these two kinds of
solutions, which receive different colors depending on the value of the
synchronization error. The same calculation is shown in Fig. 19 (b), which
shows a zoom on the center of Fig. 19 (a).

This broad error range is due to the different types of solutions that the
system allows. In Fig. 20 we show the time-series of the two first nodes in
red and blue lines. In Fig. 20(a), a synchronized solution, which yields Err =
$0$ (taking into account only $x_{1}$ and $x_{2}$), for (b) and (c) the
synchronization error is equal to $2.059$ and $7.964$. The larger
synchronization error comes from the last type of solution, depicted in (d),
where the two nodes are out of phase, resulting in an Err $=20.884$. In
addition to that, in the last case we notice a hard increase in the frequency
of spikes.

Fig. 21 depicts the attractors of five different solutions color-coded in the
same fashion of Fig. 19. The attractors corresponding to solutions with Err
$\leq 10$ are superposed in the same region, while the yellow and blue
attractor, which yields an Err between $10$ and $20$, are far apart and have
different shape.

Figure 19: Riddled basin of attraction in $(x_{1},x_{2})$ plane for $g=0.18$.
(a) is calculated with $10.24\times 10^{4}$ initial conditions in
$x_{1},y_{1}\in[1,-1]$. (b) Inset of (a), calculated with $2.56\times 10^{4}$
initial conditions in $x_{1},y_{1}\in[0.01,-0.01]$. Figure 20: Time series of
the $x$-variable from nodes $1$ (red) and $2$ (blue) for different initial
conditions. (a) Complete synchronization with Err $=0$, (b) and (c) burst
synchronization with Err $=2.05$ and Err $=7.96$, and (d) out of phase
synchronization with Err $=20.88$. Figure 21: Depending on the initial
conditions the system settles into different attractors. Each attractor is
color-coded according to the synchronization error in Fig. 19.

Figure 22: Riddled basin of attraction in $(x_{1},x_{2})$ plane for $g=0.26$.
(a) is calculated with $2.56\times 10^{4}$ initial conditions in
$x_{1},y_{1}\in[1,-1]$. (b) Inset of (a), calculated with $2.56\times 10^{4}$
initial conditions in $x_{1},y_{1}\in[0.01,-0.01]$.

It is worth mentioning that the transition to synchronization observed in Fig.
18 is due to the nearly identical initial conditions. Taking $g=0.26$, and
using initial conditions as in Fig. 22, yields a bistable system, as we see in
Fig. 22 (a) and (b). Instead of complete synchronization (Fig. 18) or multiple
solutions Fig. (19), in this case, the system is either in a completely
synchronized (Err $\approx 0$) or in an out of phase synchronized (Err
$\approx 20$) state. The inset in Fig. 22(b) shows that points from these two
attractors are interwoven in a riddled fashion.

##  IV Conclusions

We have calculated the MSF of networks of Izhikevich neurons, considering
global and cluster synchronized states with electrical and chemical coupling
schemes. To do so, we combined the MSF formalism with the use of saltation
matrices, which allow the calculation of Lyapunov exponents of systems with
discontinuities such as the Izhikevich model. For the networks studied, we
found that the addition of electrical coupling to a network with only chemical
coupling induces synchronization for both global and cluster synchronized
states. However, the presence of both coupling schemes yielded a riddled basin
of attraction. Due to this, even in the range where the MSF is negative and
with nearly identical initial conditions concerning the CSS, the system can
end up in a not synchronized state and with a different kind of dynamics.

It will be interesting, for future work, to study not only networks with more
elements but also how the stability of synchronized solutions if we consider
non-identical neurons and the presence of noise in the network. It would also
be appealing to investigate why the riddled basins displayed by these systems
appear as the MLE approaches zero from above, and not from below as it is
usually reported in the literature.

##  V Acknowledgements

We acknowledge useful discussions with Federico Bizzarri, Bruno Boaretto and
Eleonora Catsigeras. R.P.A. acknowledges financing from Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES), Finance Code
001. H.A.C. thanks ICTP-SAIFR and FAPESP Grant No. 2016/01343-7 for partial
support.

##  Appendix A Derivation of MSF for the chemical coupling

Starting from Eq. 12 we define the $2\times 2$ matrices

$\mathbf{H}=\begin{bmatrix}h&0\\\
0&0\end{bmatrix}\quad\text{and}\quad\mathbf{C}=\begin{bmatrix}\zeta&0\\\
0&0\end{bmatrix}\ ,$ (39)

where $h$ and $\zeta$ are nonlinear operators

$\displaystyle hx$ $\displaystyle=h(x)=x-v_{s}\ ,$ (40) $\displaystyle\zeta
x=\zeta(x)$ $\displaystyle=[1+\exp(-\epsilon(x-\theta))]^{-1}\ ,$ (41)

and an operator $\mathbf{Dg}$, that returns a diagonal matrix whose elements
are the components of the vector on which it acts. For example, for
$\mathbf{x}$:

$\mathbf{Dg}(\mathbf{x})=\mathbf{Dg}\begin{bmatrix}x^{1}\\\ y^{1}\\\ \vdots\\\
x^{N}\\\ y^{N}\end{bmatrix}=\begin{bmatrix}x^{1}&0&...&0&0\\\
0&y^{1}&...&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&...&x^{N}&0\\\
0&0&...&0&y^{N}\end{bmatrix}\ ,$ (42)

where in this notation $\mathbf{Dg}(\mathbf{x})$ reads $\mathbf{Dg}$ acts on
$\mathbf{x}$. We can write the vector encoding the synaptic coupling of a
network of $N$ nodes as

$\mathbf{Dg}[(\mathbf{I}_{N}\otimes\mathbf{H})(\mathbf{x})](A\otimes\mathbf{C})(\mathbf{x})\
.$ (43)

Moreover, it is convenient to introduce the matrices

$\mathbf{T}=\mathbf{Dg}(\mathbf{I}_{N}\otimes\mathbf{H})$ (44)
$\mathbf{K}=\mathbf{Dg}(\mathbf{I}_{N}\otimes\mathbf{C})\ .$ (45)

With this notation, the dynamics of a network coupled through chemical
synapses reads

$\mathbf{\dot{x}}=\mathbf{F}(\mathbf{x})-g_{c}\mathbf{T}\mathbf{x}[(A\otimes\mathbf{C})(\mathbf{x})]\
.$ (46)

For example, if $N=1$,

$\mathbf{T}\mathbf{x}=\mathbf{Dg}\Bigg{[}\Bigg{(}\mathbf{I}_{1}\otimes\begin{bmatrix}h&0\\\
0&0\end{bmatrix}\Bigg{)}\begin{bmatrix}x^{1}\\\
y^{1}\end{bmatrix}\Bigg{]}=\mathbf{Dg}\begin{bmatrix}h(x^{1})\\\
0\end{bmatrix}$ $=\begin{bmatrix}x^{1}-v_{s}&0\\\ 0&0\end{bmatrix}$ (47)

If we consider a network with $N=2$ neurons and adjacency matrix given by

$A=\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}\ ,$ (48)

the coupling term reads

$\displaystyle\mathbf{T}\mathbf{x}$
$\displaystyle=\mathbf{Dg}\Bigg{(}\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\otimes\begin{bmatrix}h&0\\\
0&0\end{bmatrix}\Bigg{)}\begin{bmatrix}x^{1}\\\ y^{1}\\\ x^{2}\\\ y^{2}\\\
\end{bmatrix}=\mathbf{Dg}\begin{bmatrix}h&0&0&0\\\ 0&0&0&0\\\ 0&0&h&0\\\
0&0&0&0\\\ \end{bmatrix}\begin{bmatrix}x^{1}\\\ y^{1}\\\ x^{2}\\\ y^{2}\\\
\end{bmatrix}$ $\displaystyle\mathbf{Tx}$
$\displaystyle=\mathbf{Dg}\begin{bmatrix}h(x^{1})\\\ 0\\\ h(x^{2})\\\ 0\\\
\end{bmatrix}=\begin{bmatrix}h(x^{1})&0&0&0\\\ 0&0&0&0\\\ 0&0&h(x^{2})&0\\\
0&0&0&0\\\ \end{bmatrix}$ (49)

times

$(A\otimes\mathbf{C})(\mathbf{x})=\begin{bmatrix}0&0&\zeta&0\\\ 0&0&0&0\\\
\zeta&0&0&0\\\ 0&0&0&0\\\ \end{bmatrix}\begin{bmatrix}x^{1}\\\ y^{1}\\\
x^{2}\\\ y^{2}\\\ \end{bmatrix}=\begin{bmatrix}\zeta(x^{2})\\\ 0\\\
\zeta(x^{1})\\\ 0\\\ \end{bmatrix}$ (50)

yielding

$\displaystyle\mathbf{T}(\mathbf{x})[(A\otimes\mathbf{C})(\mathbf{x})]$
$\displaystyle=\begin{bmatrix}h(x^{1})&0&0&0\\\ 0&0&0&0\\\ 0&0&h(x^{2})&0\\\
0&0&0&0\\\ \end{bmatrix}\begin{bmatrix}\zeta(x^{2})\\\ 0\\\ \zeta(x^{1})\\\
0\\\ \end{bmatrix}$ (51)
$\displaystyle\mathbf{T}(\mathbf{x})[(A\otimes\mathbf{C})(\mathbf{x})]$
$\displaystyle=\begin{bmatrix}h(x^{1})\zeta(x^{2})\\\ 0\\\
h(x^{2})\zeta(x^{1})\\\ 0\\\ \end{bmatrix}\ .$ (52)

Following the Pecora-Carroll analysis, the master stability function can be
obtained via the variational equation for
$\delta\mathbf{x}^{i}=\mathbf{x}^{i}-\mathbf{x}^{s}$. So first, we derive the
synchronized solution $\mathbf{x}^{s}$, where $\mathbf{x}^{i}=\mathbf{x}^{j}$,
$\forall$ $i,j$. If we take an adjacency matrix equal to Eq. 16, we have

$\mathbf{\dot{x}}^{s}=\mathbf{F}(\mathbf{x}^{s})-g_{c}\mathbf{T}(\mathbf{x}^{s})[(A\otimes\mathbf{C})(\mathbf{x}^{s})]\
,$ (53)

the second term reads

$\mathbf{Dg}\Bigg{\\{}\Bigg{(}\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{bmatrix}\otimes\begin{bmatrix}h&0\\\
0&0\end{bmatrix}\Bigg{)}\begin{bmatrix}x^{s}\\\ y^{s}\\\ \vdots\\\ x^{s}\\\
y^{s}\\\ \end{bmatrix}\Bigg{\\}}\Bigg{\\{}\Bigg{(}\begin{bmatrix}0&1&0&1\\\
1&0&1&0\\\ 0&1&0&1\\\ 1&0&1&0\end{bmatrix}\otimes\begin{bmatrix}\zeta&0\\\
0&0\end{bmatrix}\Bigg{)}\begin{bmatrix}x^{s}\\\ y^{s}\\\ \vdots\\\ x^{s}\\\
y^{s}\end{bmatrix}\Bigg{\\}}$ (54)

which, after application of $\mathbf{Dg}$ becomes:

$g_{c}\begin{bmatrix}h(x^{s})&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&h(x^{s})&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&h(x^{s})&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&h(x^{s})&0\\\ 0&0&0&0&0&0&0&0\\\
\end{bmatrix}\Bigg{(}\begin{bmatrix}0&1&0&1\\\ 1&0&1&0\\\ 0&1&0&1\\\
1&0&1&0\end{bmatrix}\otimes\begin{bmatrix}\zeta&0\\\
0&0\end{bmatrix}\begin{bmatrix}x^{s}\\\ y^{s}\\\ \vdots\\\ x^{s}\\\
y^{s}\end{bmatrix}\Bigg{)}$ (55)

simplifying we obtain:

$g_{c}\begin{bmatrix}h(x^{s})&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&h(x^{s})&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&h(x^{s})&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&h(x^{s})&0\\\ 0&0&0&0&0&0&0&0\\\
\end{bmatrix}\begin{bmatrix}\zeta(x^{s})+\zeta(x^{s})\\\ 0\\\
\zeta(x^{s})+\zeta(x^{s})\\\ 0\\\ \zeta(x^{s})+\zeta(x^{s})\\\ 0\\\
\zeta(x^{s})+\zeta(x^{s})\\\ 0\\\ \end{bmatrix}$ (56)

which is symply

$g_{c}k_{n}\mathbf{T}(\mathbf{x}^{s})\mathbf{C}(\mathbf{x}^{s})$ (57)

where in this case, the number of links that each node has is $k_{n}=2$. Now,
we linearize $\mathbf{x}^{i}$ around the synchronized state $\mathbf{x}^{s}$

$\mathbf{x}^{i}\approx\mathbf{F}(\mathbf{x}^{s})+D\mathbf{F}\big{|}_{\mathbf{x}^{s}}\delta\mathbf{x}^{i}-g_{c}D\Big{[}\mathbf{T}(\mathbf{x}^{i})\sum_{j}a_{ij}\mathbf{C}(\mathbf{x}^{j})\Big{]}\delta\mathbf{x}^{k}\
,$ (58)

where $k=(i,j)$ depending on whether the Jacobian matrix $D$ acts on a
function of $\mathbf{x}^{i}$ or $\mathbf{x}^{j}$. Evaluating the last term
yields

$D[\mathbf{T}(\mathbf{x}^{i})\sum_{j}a_{ij}\mathbf{C}(\mathbf{x}^{j})]\delta\mathbf{x}^{k}=D\begin{bmatrix}h(x^{i})\sum_{j}a_{ij}\zeta(x^{j})\\\
0\end{bmatrix}\delta\mathbf{x}_{k}$ (59) $=\begin{bmatrix}\frac{\partial
h(x^{i})}{\partial x^{i}}\sum_{j}a_{ij}\zeta(x^{j})&0\\\
0&0\end{bmatrix}\delta\mathbf{x}^{i}+\begin{bmatrix}h(x^{i})\sum_{j}a_{ij}\frac{\partial\zeta(x^{j})}{\partial
x^{i}}&0\\\ 0&0\end{bmatrix}\delta\mathbf{x}^{j}$ (60)
$=\begin{bmatrix}\frac{\partial h(x^{i})}{\partial
x^{i}}k_{n}\zeta(x^{j})&0\\\
0&0\end{bmatrix}\Bigg{|}_{\mathbf{x}^{s}}\delta\mathbf{x}^{i}+\begin{bmatrix}h(x^{i})\sum_{j}a_{ij}\frac{\partial\zeta(x^{j})}{\partial
x^{i}}&0\\\ 0&0\end{bmatrix}\Bigg{|}_{\mathbf{x}^{s}}\delta\mathbf{x}^{j}$
(61)

which is

$g_{c}k_{n}D\mathbf{H}\big{|}_{\mathbf{x}^{s}}\mathbf{K}(\mathbf{x}^{s})\delta\mathbf{x}^{i}-g_{c}\mathbf{T}(\mathbf{x}^{s})\sum_{j}a_{ij}D\mathbf{C}\big{|}_{\mathbf{x}^{s}}\delta\mathbf{x}^{j}\
.$ (62)

Then, we can write $\delta\mathbf{x}^{i}$ as

$\displaystyle\delta\mathbf{x}^{i}=D\mathbf{F}\big{|}_{\mathbf{x}^{s}}-g_{c}k_{n}D\mathbf{H}\big{|}_{\mathbf{x}^{s}}\mathbf{K}(\mathbf{x}^{s})\delta\mathbf{x}^{i}-g_{c}\mathbf{T}(\mathbf{x}^{s})\sum_{j}a_{ij}D\mathbf{C}\big{|}_{\mathbf{x}^{s}}\delta\mathbf{x}^{j}\
.$ (63)

and finally, we write
$\delta\mathbf{x}=[\delta\mathbf{x}^{1},...,\delta\mathbf{x}^{N}]^{T}$ as

$\displaystyle\delta\dot{\mathbf{x}}=\\{[\mathbf{I}_{N}\otimes
D\mathbf{F}\big{|}_{\mathbf{x}^{s}}-g_{c}k_{n}D\mathbf{H}\big{|}_{\mathbf{x}^{s}}\mathbf{K}(\mathbf{x}^{s}))]-g_{c}A\otimes\mathbf{T}(\mathbf{x}^{s})D\mathbf{C}\big{|}_{\mathbf{x}^{s}}\\}\delta\dot{\mathbf{x}}\
.$ (64)

The diagonalization of $A$ leads to the desired block diagonalized variational
equation:

$\displaystyle\dot{\eta}=\\{[\mathbf{I}_{N}\otimes(D\mathbf{F}\big{|}_{\mathbf{x}^{s}}-g_{c}k_{n}D\mathbf{H}\big{|}_{\mathbf{x}^{s}}\mathbf{K}(\mathbf{x}^{s}))]-g_{c}\Gamma\otimes\mathbf{T}(\mathbf{x}^{s})D\mathbf{C}\big{|}_{\mathbf{x}^{s}}\\}\eta\
.$ (65)

## References

  * Purves  _et al._ [2009] D. Purves, G. J. Augustine, D. Fitzpatrick, W. C. Hall, A.-S. LaMantia, and R. D. Mooney,  _Neuroscience_ (Oxford, 2009).
  * Hodgkin and Huxley [1952] A. L. Hodgkin and A. F. Huxley, J. Physiol. 117, 500 (1952).
  * Izhikevich [2003] E. M. Izhikevich, IEEE Trans. Neural Net. 14, 1569 (2003).
  * Izhikevich and Edelman [2008] E. M. Izhikevich and G. M. Edelman, Proc. Natl. Acad. Sci. U.S.A. 105, 3593 (2008).
  * Modolo  _et al._ [2007] J. Modolo, E. Mosekilde, and A. Beuter, J. of Physiol. Paris 101, 56 (2007).
  * Buzsaki [2006] G. Buzsaki,  _Rhythms of the Brain_ (Oxford University Press, 2006).
  * Jiruska  _et al._ [2013] P. Jiruska, M. Curtis, J. G. R. Jefferys, C. A. Schevon, S. J. Schiff, and K. Schindler, J. Physiol. 594, 787 (2013).
  * Uhlhaas and Singer [2006] P. Uhlhaas and W. Singer, Neuron 52, 155 (2006).
  * Cohen and D’Esposito [2016] J. R. Cohen and M. D’Esposito, J. Neurosci. 36, 12083 (2016).
  * Axmacher [2006] N. Axmacher, and F. Mormann, and G. Fernández, and C. E. Elger, and J. Fell Brain. Res. Rev. 52, 170-182 (2006).
  * Fries [2009] P. Fries, Annu. Rev. of Neurosci. 32, 209 (2009).
  * Pecora and Carroll [1998] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
  * Nishikawa and Motter [2006] T. Nishikawa and A. E. Motter, Phys. Rev. E 73, 065106(R) (2006).
  * Sun  _et al._ [2009] J. Sun, E. M. Bollt, and T. Nishikawa, Europhys. Lett. 73, 60011 (2009).
  * Rossa  _et al._ [2020] F. D. Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, and F. Sorrentino, Nat. Commun. 11 (2020).
  * Pecora  _et al._ [2014] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, Nat. Commun. 5 (2014).
  * Pecora  _et al._ [2016] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, Sci. Adv. 2 (2016).
  * Schaub  _et al._ [2016] M. Schaub, N. O’Clery, Y. N. Billeh, J.-C. Delvenne, R. Lambiotte, and M. Barahona, Chaos 26, 094821 (2016).
  * Zhang  _et al._ [2021] Y. Zhang, V. Latora, and A. E. Motter, Commun. Phys. 4 (2021).
  * Alexander  _et al._ [1992] J. C. Alexander, J. A. Yorke, Y. Z., and I. Kan, Int. J. Bifurc. Chaos 2, 795 (1992).
  * Huang  _et al._ [2009] L. Huang, Q. Chen, Y.-C. Lai, and L. M. Pecora, Phys. Rev. E 80, 036204 (2009).
  * Heagy [1994] J. F. Heagy, T. L. Carroll, and L. M. Pecora, Phys. Rev. Lett. 73, (1994).
  * Ott [1994] E. Ott, J. C. Sommerer, Phys. Lett. A 188, 39-47 (1994).
  * Shimada and Nagashima [1979] I. Shimada and T. Nagashima, Prog. Theor. Phys. 61, 1605 (1979).
  * Benettin  _et al._ [1980] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 9 (1980).
  * Wolf  _et al._ [1985] A. Wolf, J. . B. Swift, H. L. Swinney, and J. A. Vastano, Phys. D. 16, 285 (1985).
  * Bizzari  _et al._ [2013] F. Bizzari, A. Brambilla, and G. S. Gajani, J. Comput. Neurosci. 35, 201 (2013).
  * di Bernardo  _et al._ [2008] M. di Bernardo, C. J. Budd, A. R. Champneys, and P. Kowalczyk,  _Piecewise-smooth dynamical systems: theory and applications_ (Springer, 2008) p. 106.
  * Perkel  _et al._ [1981] D. H. Perkel, b. Mulloney, and R. W. Budelli, Neuroscience 6, 823 (1981).
  * Somers and Kopell [1993] D. Somers and N. Kopell, Biol. Cybern. 68, 393 (1993).
  * Root  _et al._ [2014] D. H. Root, C. Mejias-Aponte, S. Zhang, H. Wang, A. F. Hoff, C. R. Lupica, and M. Morales, Nat. Neurosci. 17, 1543 (2014).
  * Pelkey  _et al._ [2020] K. A. Pelkey, D. Calvigioni, C. Fang, G. Vargish, T. Ekins, K. Auville, J. C. Wester, M. Lai, C. M.-C. Scott, X. Yuan, S. Hunt, D. Abebe, Q. Xu, J. Dimidschstein, G. Fishell, R. Chittajallu, and C. J. McBain, eLife 9, e51995 (2020).
  * Checco  _et al._ [2008] P. Checco, M. Righero, M. Biey, and L. Kocarev, IEEE Trans. Circuits Syst. II 55, 1274 (2008).
  * Feudel [2008] U. Feudel, Int. J. Bifurc. Chaos 18, 1607 (2008).
  * Cardoso  _et al._ [2007] D. M. Cardoso, C. Delorme, and P. Rama, Eur. J. Comb. 28, 665 (2007).
  * Moore [1920] E. H. Moore, Bull. Amer. Math. Soc 26, 394 (1920).
  * Penrose [1955] R. Penrose, Math. Proc. Camb. Philos. Soc. 51, 406 (1955).
  * O’Clery  _et al._ [2007] N. O’Clery, Y. Yuan, G.-B. Stan, and M. Barahona, Phys. Rev. E 88, 042805 (2007).