in $T_{\rm K}$ and $\langle T_{\rm 21}\rangle$, around $f_{\rm*,III}\sim
0.02$.

We find that indeed direct injection of cosmic rays can heat up the IGM.
However, we do not find this particular model of cosmic ray heating to be as
efficient as in the original study by Sazonov & Sunyaev (2015). In part, this
can be attributed to differences in star formation prescription, the
relaxation of the global heating assumption that we introduced here, and the
truncation of cosmic ray heating by the LW feedback. Consequently, whereas
Sazonov & Sunyaev (2015) found cosmic rays to heat the IGM by $10$ to $100$ K
at $z=15$, in our simulations $10$ K of heating at $z=15$ is only achieved for
high values of star formation rates ($f_{\rm*,III}\geq 0.02$) or high escape
mass thresholds ($M_{\rm cr}^{\rm max}\geq 10^{7}$ M☉) with the largest $z=15$
temperature increase above the adiabatic cooling solution $\Delta T_{\rm
K}=20.3$ K achieved for $M_{\rm cr}^{\rm max}=10^{8}$ M☉. While the amount of
heat deposited into the IGM by cosmic rays is found to be smaller than in
Sazonov & Sunyaev (2015), the 21-cm signal is still sensitive to this
mechanism. Our results suggest that the 21-cm signal could provide a probe of
the nature of Pop III star supernovae through the small amount of cosmic ray
heating which they produce around $z=15-20$.

To recap, we find that the direct injection model can lead to a rich IGM
thermal history and unusual features in the 21-cm signal. Above we showed that
the shape of the global signal undergoes a qualitative transition as the
maximal halo mass changes from $M_{\rm cr}^{\rm max}=10^{6}$ to $10^{7}$ M☉. A
new inflexion point is developed as a manifestation of the complex heating
history owing to the contribution of cosmic rays. Now we consider this
transition in more detail by densely sampling the parameter $M_{\rm cr}^{\rm
max}$. We run the code 50 times with equal logarithmic spacing in $M_{\rm
cr}^{\rm max}$ in the range $10^{6}-10^{7}$. The results for $T_{\rm K}$ and
$\langle T_{\rm 21}\rangle$ are shown in Fig. 10.

Figure 10: Flattened global 21-cm signal minimum from the cosmic ray heating
cut off. Variation of the IGM kinetic temperature and the global 21-cm signal
is shown colour-coded with respect to the maximum halo mass that cosmic rays
can escape from. At the highest end of the range considered, $M_{\rm cr}^{\rm
max}=10^{7}$ M☉, $T_{\rm K}$ develops a local peak due to the contribution of
cosmic rays to heating. In turn, an inflexion point appears in the global
signal. As $M_{\rm cr}^{\rm max}$ is decreased the $T_{\rm K}$ maximum merges
with the first minimum, and the inflexion point in the global signal coincides
with the minimum of the absorption trough producing a flattened absorption
feature at $M_{\rm cr}^{\rm max}=4.57\times 10^{6}$ M☉(highlighted in black).
For lower values of the critical mass, the global signal minimum gets narrower
and deeper and is no longer flattened as cosmic ray heating becomes
insignificant.

Examining the resulting thermal histories we find a clear peak in temperature
for $M_{\rm cr}^{\rm max}=10^{7}$ M☉, the upper end of our considered range.
As $M_{\rm cr}^{\rm max}$ decreases, the peak shifts to higher redshifts and
lower temperature values. Eventually, the local maximum merges with the higher
redshift minimum transforming into two inflexion points that can be seen in
the highlighted signal (black line) corresponding to $M_{\rm cr}^{\rm
max}=4.57\times 10^{6}$ M☉. For lower values of the critical mass, thermal
evolution shows the classical behaviour driven by the adiabatic cooling and Ly
$\alpha$ heating, with no extra features and no clear sign of cosmic ray
heating.

Considering the effect on the global 21-cm signal, we see that as $M_{\rm
cr}^{\rm max}$ decreases the inflexion point shifts to higher redshifts.
Eventually, the inflexion point merges with the minimum of the absorption
trough producing a flattened profile (black line in the right panel of Fig.
10). As we reduce the value of $M_{\rm cr}^{\rm max}$ even further, the
absorption trough becomes deeper, narrower, and is no longer flattened.

Intriguingly, the flattened absorption feature that we find for $M_{\rm
cr}^{\rm max}=4.57\times 10^{6}$ M☉ has a similar shape to the disputed EDGES
Low Band detection of the global signal (Bowman et al., 2018) which has a
flattened minimum at a similar redshift. We note, however, that cosmic ray
heating cannot provide a self-consistent explanation of the EDGES best-fit
signal on its own as the resulting absorption trough does not have the
required depth or the steep sides observed by EDGES. However, it does
illustrate that such flattened global signals are not inherently pathological
and can be achieved in models with two separate heating mechanisms, one active
at higher redshifts and only mildly pre-heating the IGM and the other
representing a more sustained heating that becomes efficient at later times
(such as X-ray heating, Ly $\alpha$ or CMB heating). The high-redshift
preheating could be achieved also in other models such as e.g., fractional
annihilating dark matter with a cross-section that decreases with time (Liu &
Slatyer, 2018).

##  5 Discussion and Conclusions

In this paper, we have developed a model of cosmic ray heating for use in
semi-numerical 21-cm signal simulations. Unlike previous globally-averaged
models our approach includes realistic modelling of the spatial distribution
of cosmic ray heating. We bracket the range of possible cosmic ray heating
scenarios by considering free-streaming versus locally-confined models of
particle propagation, exploring the dependence of heating on the cosmic ray
energy spectrum, and considering different channels for escape from host
halos. Our method thus allows us to perform the first investigations of the
signatures of this heating mechanism in 21-cm tomographic maps and power
spectra, as well as the previously studied global 21-cm signal. It also
enables us to compare the contribution of cosmic rays to the thermal history
of the IGM to other heating mechanisms, such as X-ray heating, and thus to
propose novel methods to potentially distinguish an X-ray heated IGM from a
cosmic-ray heated IGM using upcoming observations.

We find the spatial range of cosmic ray heating to be much shorter than that
of X-ray heating, even when cosmic rays are assumed to travel along straight
lines. This short-range nature manifests in the IGM heating being clustered
around regions of efficient star formation and results in a sharp contrast in
the 21-cm signal between localized heated regions of emission and vast
unheated regions of absorption.

Comparing our inhomogeneous cosmic ray heating model to a uniform model
commonly used in the literature, we find that its localized nature results in
a slightly delayed evolution of the global 21-cm signal. This demonstrates the
existence of a small assumption-induced bias in the results of previous
globally-averaged studies. Much larger differences are seen in the power
spectra and tomographic maps, as with cosmic rays as a dominant heating source
the 21-cm signal traces matter fluctuations to smaller scales than in the case
of X-rays. Specifically, we find that the 21-cm power spectrum at higher $k$
is increasingly suppressed for the longer-ranged heating mechanisms allowing
us to constrain the propagation mechanism of cosmic ray particles and
distinguish cosmic ray from X-ray heating. Any clustering of cosmic ray
heating on scales smaller than our simulation resolution would enhance these
effects and thus strengthen our conclusions. The potential of the high-$k$ end
of the power spectrum to provide information about the spatial range of the
dominant heating mechanism is anticipated to apply more generally and be
robust to astrophysical uncertainties, making it a useful diagnostic tool. For
example, it has previously been suggested (Fialkov & Barkana, 2014; Fialkov et
al., 2014b) that the high-$k$ end of the power spectrum could be used to probe
the SED of early X-ray binaries.

Due to the uncertainties surrounding cosmic ray escape into the IGM, we
explored several plausible mechanisms. In models where cosmic rays can escape
from halos of any mass via diffusion or outflow advection, we found cosmic ray
heating of the IGM can be more efficient than X-ray heating. However,
similarly to Leite et al. (2017), we found that the efficiency of cosmic rays
as a heating source is strongly dependent on the spectrum of the particles
injected into the IGM. A change in the source kinetic energy spectral exponent
from $-2$ to $-1.7$ decreases the IGM temperature from $284$ to $19$ K at
$z=8$. This is because only the lower energy cosmic ray protons are able to
efficiently transfer their kinetic energy to the IGM as heat, as was
previously demonstrated in Sazonov & Sunyaev (2015).

For cosmic ray heating models similar to that of Sazonov & Sunyaev (2015),
wherein cosmic rays can only escape from lower mass Pop III star-forming
halos, we found cosmic ray heating to be relatively inefficient. Due to
Sazonov & Sunyaev (2015) not considering the LW feedback in their model, we
were only able to replicate the significant increase in temperature they
observed for unphysically weak LW feedback or if we raised the threshold halo
mass of cosmic ray sources. Instead at theoretically motivated parameter
values we find the interplay between this mass threshold and the LW feedback,
which raises $M_{\rm mol,crit}$, leads to unexpectedly complex thermal
histories. In such models, the IGM kinetic temperature was found to have a
local peak due to cosmic ray heating being cut off by the LW feedback leading
to a secondary episode of cosmic cooling before Ly $\alpha$ heating from Pop
II stars becomes efficient and re-heats the IGM. This unusual thermal
evolution results in an inflexion point in the heating arm of the global 21-cm
signal. By fine-tuning the mass threshold for cosmic ray escape, this
inflexion point can be merged with the minimum of the absorption trough
resulting in a flattened profile. Thus, we demonstrate that flat-based global
21-cm signals (like the one detected by EDGES) can in principle be achieved in
models where one heating mechanism becomes inefficient at a similar time to
the gas becoming fully coupled via the WF effect and another heating mechanism
turns on at a later time. We do not advocate for this model to be a self-
consistent explanation of the EDGES best-fit signal as it does not explain the
depth or the steep sides of the detected absorption trough (although a dual-
heating scenario could form a part of a more complex explanation which also
invokes a radio contribution to create the deep absorption).

There remains a major uncertainty in cosmic ray heating, namely the strength
of the IGM primordial magnetic field which scatters cosmic rays. This
uncertainty leads to a large uncertainty in the cosmic ray diffusivity within
the IGM. Furthermore, in the case of primordial magnetic fields at the higher
end of the experimentally allowed range ($>0.01$ nG comoving field strength),
the transfer of energy from cosmic rays to the IGM via Alfvén waves (Bera et
al., 2023) and synchrotron radio emission from cosmic ray electrons (Jana et
al., 2019) are anticipated to become efficient. Since these effects are only
efficient for the largest experimentally allowed primordial magnetic fields,
we do not include them in our modelling. However, should the primordial
magnetic field take such a high value, these effects would result in stronger
heating of the IGM from cosmic rays, clustered around star-forming halos, and
enhanced radio background, the former acting to diminish the 21-cm signal and
the latter enhancing it. This sensitivity of the 21-cm signal, in particular
the 21-cm power spectrum, to cosmic ray heating of the IGM provides an
indirect probe of the primordial IGM magnetic field.

Overall we found that the distinct nature of cosmic ray heating means it
leaves characteristic signatures in 21-cm tomographic maps, the high-$k$ end
of the power spectrum, and in some cases the 21-cm global signal. These
features are of sufficient magnitude to potentially be probed by the 21-cm
signal measurements from current and next-generation experiments. Thus, an
understanding of cosmic ray heating is necessary for the correct
interpretation of the 21-cm signal, and conversely, we may be able to
constrain early universe astrophysics that impacts cosmic rays (for example,
the primordial IGM magnetic field) through the 21-cm signal.

## Acknowledgements

The authors would like to thank Donghui Huang for helpful conversations
concerning the updating of the LW feedback prescription in our simulation
code, as well as to express our gratitude to John Cumner and Boyuan Liu for
their invaluable comments on earlier versions of this paper. We would
additionally like to thank the anonymous referee, for their careful reading of
this work and insightful comments, that have led to a much improved final
manuscript.

TGJ would like to thank the Science and Technology Facilities Council (UK) for
their continued support through grant number ST/V506606/1. AF is supported by
a Royal Society University Research Fellowship #180523. EdLA acknowledges the
support of the Science and Technology Facilities Council (UK) through a
Rutherford Fellowship. WJH thanks the Royal Society for their support through
a Royal Society University Research Fellowship. RB acknowledges the support of
the Israel Science Foundation (grant No. 2359/20), the Ambrose Monell
Foundation, the Institute for Advanced Study, the Vera Rubin Presidential
Chair in Astronomy, and the Packard Foundation.

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## Data Availability

The code to calculate the window functions used in this study to model cosmic
ray propagation is made available at
https://github.com/ThomasGesseyJones/CosmicRayHeatingFor21cm. All other data
used and generated in the writing of this article will be shared on reasonable
request to the corresponding author.

## References

  * Abdurashidova et al. (2022) Abdurashidova Z., et al., 2022, ApJ, 925, 221
  * Ackermann et al. (2013) Ackermann M., et al., 2013, Science, 339, 807
  * Alves Batista (2022) Alves Batista R., 2022, arXiv e-prints, p. arXiv:2210.12855
  * Barkana (2016) Barkana R., 2016, Phys. Rep., 645, 1
  * Barkana (2018) Barkana R., 2018, Nature, 555, 71
  * Barkana & Loeb (2004) Barkana R., Loeb A., 2004, ApJ, 609, 474
  * Barkana & Loeb (2005) Barkana R., Loeb A., 2005, ApJ, 626, 1
  * Barkana et al. (2023) Barkana R., Fialkov A., Liu H., Outmezguine N. J., 2023, Phys. Rev. D, 108, 063503
  * Bell (1978a) Bell A. R., 1978a, MNRAS, 182, 147
  * Bell (1978b) Bell A. R., 1978b, MNRAS, 182, 443
  * Bera et al. (2023) Bera A., Samui S., Datta K. K., 2023, MNRAS, 519, 4869
  * Berezinskii et al. (1990) Berezinskii V. S., Bulanov S. V., Dogiel V. A., Ptuskin V. S., 1990, Astrophysics of cosmic rays. Elsevier Science & Technology
  * Bowman et al. (2018) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67
  * Burns (2021) Burns J. O., 2021, Philosophical Transactions of the Royal Society of London Series A, 379, 20190564
  * Butler & Buckingham (1962) Butler S. T., Buckingham M. J., 1962, Physical Review, 126, 1
  * Caprioli & Spitkovsky (2014) Caprioli D., Spitkovsky A., 2014, ApJ, 783, 91
  * Caprioli et al. (2010) Caprioli D., Amato E., Blasi P., 2010, Astroparticle Physics, 33, 160
  * Chen et al. (2014) Chen K.-J., Heger A., Woosley S., Almgren A., Whalen D. J., 2014, ApJ, 792, 44
  * Chen et al. (2017) Chen K.-J., Heger A., Whalen D. J., Moriya T. J., Bromm V., Woosley S. E., 2017, MNRAS, 467, 4731
  * Choudhury (2022) Choudhury T. R., 2022, General Relativity and Gravitation, 54, 102
  * Cohen et al. (2016) Cohen A., Fialkov A., Barkana R., 2016, MNRAS, 459, L90
  * Drury et al. (1989) Drury L. O., Markiewicz W. J., Voelk H. J., 1989, A&A, 225, 179
  * Ellingson et al. (2009) Ellingson S. W., Clarke T. E., Cohen A., Craig J., Kassim N. E., Pihlstrom Y., Rickard L. J., Taylor G. B., 2009, IEEE Proceedings, 97, 1421
  * Ewall-Wice et al. (2016) Ewall-Wice A., et al., 2016, MNRAS, 460, 4320
  * Ewall-Wice et al. (2018) Ewall-Wice A., Chang T. C., Lazio J., Doré O., Seiffert M., Monsalve R. A., 2018, ApJ, 868, 63
  * Feng & Holder (2018) Feng C., Holder G., 2018, ApJ, 858, L17
  * Fialkov & Barkana (2014) Fialkov A., Barkana R., 2014, MNRAS, 445, 213
  * Fialkov & Barkana (2019) Fialkov A., Barkana R., 2019, MNRAS, 486, 1763
  * Fialkov et al. (2012) Fialkov A., Barkana R., Tseliakhovich D., Hirata C. M., 2012, MNRAS, 424, 1335
  * Fialkov et al. (2013) Fialkov A., Barkana R., Visbal E., Tseliakhovich D., Hirata C. M., 2013, MNRAS, 432, 2909
  * Fialkov et al. (2014a) Fialkov A., Barkana R., Pinhas A., Visbal E., 2014a, MNRAS, 437, L36
  * Fialkov et al. (2014b) Fialkov A., Barkana R., Visbal E., 2014b, Nature, 506, 197
  * Fialkov et al. (2018) Fialkov A., Barkana R., Cohen A., 2018, Phys. Rev. Lett., 121, 011101
  * Field (1958) Field G. B., 1958, Proceedings of the IRE, 46, 240
  * Fragos et al. (2013) Fragos T., et al., 2013, ApJ, 764, 41
  * Fraser et al. (2018) Fraser S., et al., 2018, Physics Letters B, 785, 159
  * Furlanetto & Oh (2005) Furlanetto S. R., Oh S. P., 2005, MNRAS, 363, 1031
  * Furlanetto et al. (2004) Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181
  * Gessey-Jones et al. (2022) Gessey-Jones T., et al., 2022, MNRAS, 516, 841
  * Ginzburg & Ozernoi (1966) Ginzburg V. L., Ozernoi L. M., 1966, Soviet Ast., 9, 726
  * Gould (1972) Gould R. J., 1972, Physica, 62, 555
  * Greif et al. (2007) Greif T. H., Johnson J. L., Bromm V., Klessen R. S., 2007, ApJ, 670, 1
  * Haiman et al. (2000) Haiman Z., Abel T., Rees M. J., 2000, ApJ, 534, 11
  * Hegde & Furlanetto (2023) Hegde S., Furlanetto S. R., 2023, MNRAS, 525, 428
  * Heger & Woosley (2002) Heger A., Woosley S. E., 2002, ApJ, 567, 532
  * Heger & Woosley (2010) Heger A., Woosley S. E., 2010, ApJ, 724, 341
  * Hibbard et al. (2022) Hibbard J. J., Mirocha J., Rapetti D., Bassett N., Burns J. O., Tauscher K., 2022, ApJ, 929, 151
  * Hills et al. (2018) Hills R., Kulkarni G., Meerburg P. D., Puchwein E., 2018, Nature, 564, E32
  * Hummel et al. (2016) Hummel J. A., Stacy A., Bromm V., 2016, MNRAS, 460, 2432
  * Ichiki et al. (2006) Ichiki K., Takahashi K., Ohno H., Hanayama H., Sugiyama N., 2006, Science, 311, 827
  * Jana & Nath (2018) Jana R., Nath B. B., 2018, MNRAS, 479, 153
  * Jana et al. (2019) Jana R., Nath B. B., Biermann P. L., 2019, MNRAS, 483, 5329
  * Jasche et al. (2007) Jasche J., Ciardi B., Enßlin T. A., 2007, MNRAS, 380, 417
  * Jones et al. (2021) Jones D., Palatnick S., Chen R., Beane A., Lidz A., 2021, ApJ, 913, 7
  * Kitayama & Yoshida (2005) Kitayama T., Yoshida N., 2005, ApJ, 630, 675
  * Koopmans et al. (2015) Koopmans L., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1 (arXiv:1505.07568), doi:10.22323/1.215.0001
  * Kulkarni et al. (2021) Kulkarni M., Visbal E., Bryan G. L., 2021, ApJ, 917, 40
  * Lacki (2015) Lacki B. C., 2015, MNRAS, 448, L20
  * Lagage & Cesarsky (1983) Lagage P. O., Cesarsky C. J., 1983, A&A, 125, 249
  * Leite et al. (2017) Leite N., Evoli C., D’Angelo M., Ciardi B., Sigl G., Ferrara A., 2017, MNRAS, 469, 416
  * Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
  * Lewis & Challinor (2011) Lewis A., Challinor A., 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library, record ascl:1102.026 (ascl:1102.026)
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Liu & Slatyer (2018) Liu H., Slatyer T. R., 2018, Phys. Rev. D, 98, 023501
  * Liu et al. (2019) Liu H., Outmezguine N. J., Redigolo D., Volansky T., 2019, Phys. Rev. D, 100, 123011
  * Longair (1994) Longair M. S., 1994, High energy astrophysics. Vol.2: Stars, the galaxy and the interstellar medium. Cambridge University Press
  * Ma et al. (2020) Ma X., Quataert E., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2020, MNRAS, 498, 2001
  * Madau et al. (1997) Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429
  * Magg et al. (2022a) Magg M., Hartwig T., Chen L.-H., Tarumi Y., 2022a, A-SLOTH: Semi-analytical model to connect first stars and galaxies to observables, Astrophysics Source Code Library, record ascl:2209.001 (ascl:2209.001)
  * Magg et al. (2022b) Magg M., et al., 2022b, MNRAS, 514, 4433
  * Mebane et al. (2018) Mebane R. H., Mirocha J., Furlanetto S. R., 2018, MNRAS, 479, 4544
  * Mebane et al. (2020) Mebane R. H., Mirocha J., Furlanetto S. R., 2020, MNRAS, 493, 1217
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Mertens et al. (2021) Mertens F. G., Semelin B., Koopmans L. V. E., 2021, in Siebert A., et al., eds, SF2A-2021: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 211–214 (arXiv:2109.10055)
  * Mesinger (2019) Mesinger A., 2019, The Cosmic 21-cm Revolution; Charting the first billion years of our universe. IOP Publishing, doi:10.1088/2514-3433/ab4a73
  * Mesinger et al. (2011) Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955
  * Mesinger et al. (2013) Mesinger A., Ferrara A., Spiegel D. S., 2013, MNRAS, 431, 621
  * Mirocha & Furlanetto (2019) Mirocha J., Furlanetto S. R., 2019, MNRAS, 483, 1980
  * Mirocha et al. (2018) Mirocha J., Mebane R. H., Furlanetto S. R., Singal K., Trinh D., 2018, MNRAS, 478, 5591
  * Muñoz et al. (2018) Muñoz J. B., Dvorkin C., Loeb A., 2018, Phys. Rev. Lett., 121, 121301
  * Muñoz et al. (2020) Muñoz J. B., Dvorkin C., Cyr-Racine F.-Y., 2020, Phys. Rev. D, 101, 063526
  * Muñoz et al. (2022) Muñoz J. B., Qin Y., Mesinger A., Murray S. G., Greig B., Mason C., 2022, MNRAS, 511, 3657
  * Nath & Biermann (1993) Nath B. B., Biermann P. L., 1993, MNRAS, 265, 241
  * Nebrin et al. (2023) Nebrin O., Giri S. K., Mellema G., 2023, MNRAS, 524, 2290
  * Ohira (2022) Ohira Y., 2022, ApJ, 929, 106
  * Ohira & Murase (2019) Ohira Y., Murase K., 2019, Phys. Rev. D, 100, 061301
  * Ohira et al. (2010) Ohira Y., Murase K., Yamazaki R., 2010, A&A, 513, A17
  * Pacucci et al. (2014) Pacucci F., Mesinger A., Mineo S., Ferrara A., 2014, MNRAS, 443, 678
  * Parsons et al. (2014) Parsons A. R., et al., 2014, ApJ, 788, 106
  * Philip et al. (2019) Philip L., et al., 2019, Journal of Astronomical Instrumentation, 8, 1950004
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A19
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Price et al. (2018) Price D. C., et al., 2018, MNRAS, 478, 4193
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Ptuskin (2012) Ptuskin V., 2012, Astroparticle Physics, 39, 44
  * Reis et al. (2020) Reis I., Fialkov A., Barkana R., 2020, MNRAS, 499, 5993
  * Reis et al. (2021) Reis I., Fialkov A., Barkana R., 2021, MNRAS, 506, 5479
  * Reis et al. (2022) Reis I., Barkana R., Fialkov A., 2022, MNRAS, 511, 5265
  * Samui et al. (2018) Samui S., Subramanian K., Srianand R., 2018, MNRAS, 476, 1680
  * Sartorio et al. (2023) Sartorio N. S., et al., 2023, MNRAS, 521, 4039
  * Sazonov & Sunyaev (2015) Sazonov S., Sunyaev R., 2015, MNRAS, 454, 3464
  * Schauer et al. (2019) Schauer A. T. P., Liu B., Bromm V., 2019, ApJ, 877, L5
  * Schauer et al. (2021) Schauer A. T. P., Glover S. C. O., Klessen R. S., Clark P., 2021, MNRAS, 507, 1775
  * Schlickeiser (2002) Schlickeiser R., 2002, Cosmic Ray Astrophysics. Springer Berlin
  * Schure et al. (2012) Schure K. M., Bell A. R., O’C Drury L., Bykov A. M., 2012, Space Sci. Rev., 173, 491
  * Scott & Rees (1990) Scott D., Rees M. J., 1990, MNRAS, 247, 510
  * Seager et al. (2011) Seager S., Sasselov D. D., Scott D., 2011, RECFAST: Calculate the Recombination History of the Universe, Astrophysics Source Code Library, record ascl:1106.026 (ascl:1106.026)
  * Semelin et al. (2023) Semelin B., et al., 2023, A&A, 672, A162
  * Sikder et al. (2023) Sikder S., Barkana R., Fialkov A., Reis I., 2023, arXiv e-prints, p. arXiv:2301.04585
  * Singh et al. (2022) Singh S., et al., 2022, Nature Astronomy, 6, 607
  * Sitwell et al. (2014) Sitwell M., Mesinger A., Ma Y.-Z., Sigurdson K., 2014, MNRAS, 438, 2664
  * Sivukhin (1965) Sivukhin D. V., 1965, Reviews of Plasma Physics, 1, 1
  * Sokolowski et al. (2015) Sokolowski M., et al., 2015, Publ. Astron. Soc. Australia, 32, e004
  * Spitzer & Scott (1969) Spitzer Lyman J., Scott E. H., 1969, ApJ, 158, 161
  * Stacy & Bromm (2007) Stacy A., Bromm V., 2007, MNRAS, 382, 229
  * Strong & Moskalenko (1998) Strong A. W., Moskalenko I. V., 1998, ApJ, 509, 212
  * Sur et al. (2010) Sur S., Schleicher D. R. G., Banerjee R., Federrath C., Klessen R. S., 2010, ApJ, 721, L134
  * Tanaka & Hasegawa (2021) Tanaka T., Hasegawa K., 2021, MNRAS, 502, 463
  * Tanaka et al. (2018) Tanaka T., Hasegawa K., Yajima H., Kobayashi M. I. N., Sugiyama N., 2018, MNRAS, 480, 1925
  * Tavecchio et al. (2011) Tavecchio F., Ghisellini G., Bonnoli G., Foschini L., 2011, MNRAS, 414, 3566
  * Tegmark et al. (1997) Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., Palla F., 1997, ApJ, 474, 1
  * Venumadhav et al. (2018) Venumadhav T., Dai L., Kaurov A., Zaldarriaga M., 2018, Phys. Rev. D, 98, 103513
  * Visbal et al. (2012) Visbal E., Barkana R., Fialkov A., Tseliakhovich D., Hirata C. M., 2012, Nature, 487, 70
  * Voytek et al. (2014) Voytek T. C., Natarajan A., Jáuregui García J. M., Peterson J. B., López-Cruz O., 2014, ApJ, 782, L9
  * Whalen et al. (2008) Whalen D., Van Veelen B., O’Shea B. W., Norman M. L., 2008, in Hunt L. K., Madden S. C., Schneider R., eds, Proceedings of the International Astronomical Union Vol. 255, Low-Metallicity Star Formation: From the First Stars to Dwarf Galaxies. pp 116–120 (arXiv:0808.0524), doi:10.1017/S1743921308024678
  * Whalen et al. (2013a) Whalen D. J., Joggerst C. C., Fryer C. L., Stiavelli M., Heger A., Holz D. E., 2013a, ApJ, 768, 95
  * Whalen et al. (2013b) Whalen D. J., et al., 2013b, ApJ, 777, 110
  * Wise et al. (2014) Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560
  * Woosley (2010) Woosley S. E., 2010, ApJ, 719, L204
  * Woosley & Heger (2007) Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269
  * Woosley & Weaver (1986) Woosley S. E., Weaver T. A., 1986, ARA&A, 24, 205
  * Wouthuysen (1952) Wouthuysen S. A., 1952, AJ, 57, 31
  * Wyithe & Loeb (2004) Wyithe J. S. B., Loeb A., 2004, Nature, 432, 194
  * Xia et al. (2022) Xia Z.-Q., Wang Y., Yuan Q., Fan Y.-Z., 2022, arXiv e-prints, p. arXiv:2210.13052
  * Yajima & Khochfar (2015) Yajima H., Khochfar S., 2015, MNRAS, 448, 654
  * Yajima et al. (2011) Yajima H., Choi J.-H., Nagamine K., 2011, MNRAS, 412, 411
  * Yajima et al. (2014) Yajima H., Li Y., Zhu Q., Abel T., Gronwall C., Ciardullo R., 2014, MNRAS, 440, 776
  * Yokoyama & Ohira (2022) Yokoyama S. L., Ohira Y., 2022, MNRAS, 515, 5467
  * Yokoyama & Ohira (2023) Yokoyama S. L., Ohira Y., 2023, MNRAS, 523, 3671
  * Zahn et al. (2011) Zahn O., Mesinger A., McQuinn M., Trac H., Cen R., Hernquist L. E., 2011, MNRAS, 414, 727
  * de Lera Acedo et al. (2022) de Lera Acedo E., et al., 2022, Nature Astronomy, 6, 984


# Resilience in multiplex networks by addition of cross-repulsive links

Suman Saha, S. Saha is with National Brain Research Centre, NH-8, Manesar,
Gurugram-122052, India
E-mail<EMAIL_ADDRESS>
###### Abstract

A multiplex network of identical dynamical units becomes resilient against
parameter perturbation by adding selective linear diffusive cross-coupling
links. A parameter drift at any instant in one or multiple network nodes can
destroy synchrony, causing failure and even collapse in the network
performance. We introduced [Phys. Rev. E95, 062204(2017)] a recovery strategy
by selective addition of cross-coupling links to save synchrony in the network
from the edge of failure due to parameter mismatch (small or large) in any
nodes. This concept is extended to 2-layered multiplex networks when the
emergent synchrony becomes resilient against a small or large parameter
drifting. In addition, the stability of the synchronous state is enhanced from
local stability to global stability of synchrony. By the addition of cross-
coupling, the network revives complete synchrony in all the nodes except the
perturbed nodes, which emerges into a type of generalized synchrony with all
the unperturbed nodes. The generalized synchrony is manifested simply by a
linear amplitude response in the state variable(s) of the perturbed node(s) by
a scaling factor proportional to the mismatch. A set of systematic rules has
been derived from the linear flow matrix of the dynamical system representing
the nodes’ dynamics that helps find the connectivity matrix of the cross-
coupling links. Lyapunov function stability condition is used to determine the
cross-coupling link strength that, in turn, establishes global stability of
synchrony of the multiplex network. We verify the efficacy of our proposed
coupling scheme with analytical results and numerical simulations of two
examples of multiplex networks. In the first example, we use nonlocal
connectivity in each layer with nodal dynamics of the FitzHugh-Nagumo neuron
model. Furthermore, we use a second example of a two-layered multiplex network
of chaotic Rössler systems with random connectivity in both layers. We confirm
that our results are generic, independent of nodes’ dynamics and network
topology.

###### Index Terms:

Resilient network, synchrony, cross-coupling links, parameter perturbation,
globally stable synchrony, multiplex network .

##  1 Introduction

Multilayer networks [1, 2, 3, 4, 5] evolve in many real-world interactions
such as the transport systems consisting of many layers (road, air, and
railway) [6, 7], human disease network [8], ecological network [9], networks
of genes and tissues [10, 11], cortical and hippocampal neuronal network [12]
and human brain network [13, 14]. Multiplex network is a subclass of
multilayer networks, where one-to-one interactions between all the nodes in
different layers are maintained [15, 16]. In recent time, collective dynamics
in multi-layered networks has been extensively explored [15, 1, 17, 18, 19,
20, 21, 22] including investigations on chimera states [23, 24] and explosive
synchronization [25, 26]. In such multilayer dynamical networks [27, 28, 29],
synchrony is one of the most important desired states for all practical
purposes since it is essential for the efficient performance of the networks.
Stable synchrony in a network (single layer, multilayer) is always hindered by
internal or external perturbations when the network fails to maintain its
desired performance. The best example is the breakdown of a power grid (a
multilayer network) due to frequency drift in local nodes, generation, or
transmission centers. For efficient functioning, a spontaneous order should
sustain that demands resilience of the network synchrony against a disorder
(say, a parameter drift).

A network of identical dynamical units [30, 31] emerges into complete
synchrony (CS) for a coupling larger than a critical value [32, 33] when all
the units oscillate in a common rhythm with identical amplitude and phase for
all time. On the contrary, for heterogeneity in system parameters (induced or
drifting) of a network, CS is broken. At best, a phase synchrony (PS) [30, 34]
may be realized then. In the PS state, the oscillators are phase coherent with
almost no amplitude correlation when amplitude distortion of the individual
oscillators cannot be avoided. If the system parameters of the dynamical units
vary largely, hardly any synchrony can be realized in the network when the
individual units oscillate in complete incoherence; even a quenching of
oscillation [35, 36, 37, 38] can be evidenced, leading to an extreme situation
of a collapse of the network performance. The diversity of parameters as
heterogeneity in a system, in absence/presence of noise, may even produce a
resonance-like phenomenon with amplification of amplitude in an ensemble of
oscillators, bistable or excitable [39], and limit cycle or chaotic systems
[40, 41]. Similar situations of large parameter variation or drift may occur
in multiple nodes in a multilayer network when synchrony can never be
achieved, even for large coupling. Therefore, the stability of synchrony as a
resilience [42] against parameter drift in local nodes of a network is to be
established, i.e., the robustness of synchrony against parameter drift is the
most desirable. Recently, encouraging reports are coming that create optimism
about the constructive role of heterogeneity on synchrony in networks of non-
identical dynamical units [43, 44, 45, 46]. However, none of the reports
addressed how to restore synchrony against a parameter drift. We address the
question here to ensure the resilience of synchrony against parameter drifting
in vulnerable multiplex networks.

A common strategy of restoring synchrony in networks, in general, is to
introduce long-range interactions or links by randomizing the connectivity
between any pair of nodes or by rewiring of links [47]. A few good benefits of
rewiring have been reported in the literature, such as synchrony of dynamical
networks [47, 48, 49, 50], stopping infection progress, and enhancing recovery
in an epidemic network [51]. A recent study [52] has also proposed
optimization of synchrony in static multiplex networks by rewiring. The
resilience of synchrony was also realized against significant parameter
perturbation by varying the rewiring frequency [53] in time-varying multiplex
networks. However, the rewiring of links changes the topology of a network,
which is not always allowed or practically feasible in real-world networks of
fixed network connectivities. We search here for an alternative strategy to
realize the resilience of a multiplex network against significant parameter
perturbation in one/multiple nodes in any layer by keeping the original
network connectivity intact.

We address the problem in the following manner: Assume a two-layered multiplex
network of identical dynamical nodes performing in a stable CS state with its
stability conditions determined by the master stability function (MSF) [32],
which ensures local stability of synchrony. The critical question is how to
stop a breakdown of CS due to parameter drift in one or multiple nodes in any
layer? We propose a strategy to add directed cross-coupling links selectively
to the perturbed node(s) of one layer from the nearest/adjacent node of the
unperturbed layer and thereby try to prevent the breakdown of synchrony of the
multiplex network. The choice of the additional directed link(s) is not
arbitrary but made by following a set of systematic rules [54]. Then we define
a Lyapunov function of the error dynamics of the network and derive the
conditions for establishing the stability of synchrony. Lyapunov function-
based stability (LFS) conditions help derive the cross-coupling strength,
establishing globally stable network synchrony. The multiplex network then
becomes resilient against parameter perturbation (small or large) [54, 55].
The cross-coupling links from the unperturbed layer pass additional
information to the perturbed layer when the network performance starts
deteriorating due to parameter drifting, thereby preventing a network
breakdown.

This approach may have implications in the design of resilient multiplex
networks, in general, with its broad applicability in various physical,
biological, and engineering [4, 5, 6, 7, 8, 10, 9, 11, 12, 13, 14] systems.
Although in a different context, the constructive role of cross-coupling has
been demonstrated in organic chemical reactions to process novel compounds for
use [56], ecological species for stabilization [57], and in the origin of
chimera states in complex Ginzburg-Landau systems [58, 59], other model
systems [60]. Here we explore a general procedure of constructing a resilient
multiplex network with an analytical description of the choice of self- and
cross-coupling and their strength of coupling. The efficacy of our proposition
is elaborated with numerical examples of two different network topologies and
two model systems as representative dynamical units, one periodic and another
chaotic model. It implies that the proposition is independent of nodes’
dynamical behavior, to some extent.

Before presenting the technical and analytical details, we brief our main
results. We consider first a multiplex network with nodal dynamics represented
by the slow-fast FitzHugh-Nagumo (FHN) neuronal model [61]. The intra-layer
coupling is assumed as nonlocal (say, as an example) and unidirectional
chemical synaptic type [62]. The interlayer coupling is diffusive
bidirectional, i.e., both the nodes send/receive information to each other
simultaneously. The local stability of CS in both the layers is first
established using the conditions derived by the MSF [32] and considering all
the dynamical units of the network as identical. That is, all the nodes follow
the same dynamical behavior. Next, a parameter of a single node in one layer
is detuned (increased or decreased), which destabilizes CS and destroys the
common rhythm. Then a selective, directed diffusive cross-coupling link is
added to the perturbed node from the adjacent neighbor of the unperturbed
layer. The LFS conditions determine the strength of the cross-coupling link.
The addition of the cross-coupling link restores synchrony as emergent
generalized synchrony (GS) of the perturbed node with all other nodes while
all the other unperturbed nodes remain in the CS state. We clarify later while
elaborating the examples on what we mean by a cross-coupling link. The
emergent GS [63, 54, 41, 64] as manifested here in the perturbed node is
explained here. The amplitude of a state variable of the perturbed node is
scaled up/down (amplified/attenuated) by the ratio between detuned parameter
and original parameter of the identical nodes, at the same time, phases remain
entirely coherent among all nodes. The frequency of the perturbed node thus
remains unchanged, which is an essential condition for many real-world
networks such as the power grid that does not allow a significant frequency
drift. Results are successfully verified for perturbation in multiple nodes in
one layer.

We briefly discuss the design procedure of selecting the coupling links in the
Supplementary material (Sec. 1.1, see for some details) and our previous
report [54]. The coupling profile or the complete set of coupling functions
(self- and cross-coupling) has been identified from a knowledge of the linear
flow matrix (LFM) [54] of a nonlinear dynamical system representing the nodal
dynamics. The flow of a nonlinear dynamical system, in general, can be
separated into its linear and nonlinear components. We extract information
from the linear flow matrix of the isolated system that defines the dynamics
of a node, how to choose the appropriate self-coupling of the nodes to realize
a stable CS in a network of identical systems. The LFM also tells us the
connectivity of the additional cross-coupling links to establish global
stability of synchrony. By adopting this procedure, the desired synchrony of
the multiplex network can be maintained even when all the nodes of one layer
show a distributed heterogeneity. However, each perturbed node needs an
additional directed cross-coupling link from the other layer. The perturbed
nodes then emerge into a GS state with the unperturbed nodes, which remain in
a stable CS state. This phenomenon of emergent GS in the multiplex network as
a whole and its resilience against parameter perturbation has further been
exemplified in a random multiplex network of chaotic Rössler oscillators. Thus
the efficacy of our strategy is not restricted to non-locally coupled networks
or a particular dynamical system. The coupling profile only changes with the
dynamics of the nodes of a network. We emphasize that the multiplex network
maintains globally stable CS after adding the cross-coupling links in an
identical parametric condition. The network emerges into a globally stable GS
under perturbed conditions.

##  2 Multiplex network of FHN systems

Figure 1a shows a multiplex network of two layers where each layer is a non-
locally coupled ring of oscillators. One-to-one inter-layer bidirectional
interactions are established between the two layers. The solid (black) and
dashed (blue) lines indicate the intra- and inter-layer links, respectively.
Intra-layer couplings are unidirectional chemical synaptic type (solid black
lines) and used for nonlocal self-coupling [65, 66, 67] between each node. The
coupling radius is $P=5$ (total neighbors of a node are $2P$) in a layer of
$N=50$ nodes (total number of nodes in the network, $2N=100$). The coupling
radius $P$ and the network size $N$ may be increased without any change in the
results. For this network topology, the dynamics of each node is represented
by the slow-fast FHN [61] system. The coupled systems of $i^{th}$ node in
layer-$l$ is,

$\displaystyle\dot{x}_{i}^{l}=$ $\displaystyle
F_{x}(x^{l}_{i},y^{l}_{i})+\frac{\varepsilon_{1}(x^{l}_{i}-V_{s})}{2P}\sum\limits_{j=i-P}^{i+P}H(x^{l}_{j})$
$\displaystyle+\varepsilon_{2}Q(x_{i}^{(n,l)}-x_{i}^{(l,n)})+\kappa~{}C_{i,i}^{(n,l)}(y_{i}^{n}-y_{i}^{l}),$
$\displaystyle\dot{y}_{i}^{l}=$
$\displaystyle\;F_{y}(x^{l}_{i},y^{l}_{i},r^{l}_{i}),$ (1)

where $x_{i}$ and $y_{i}$ are the state variables of the $i^{th}$ node,
$i=1,2,...,N$; $N$ is the number of nodes in both the layers with indices
$n,l=1,2$ ($n\neq l$). Subscript indicates the number of nodes in a layer and
superscript denotes the number of layers. $V_{s}$ is a constant bias. The flow
of intrinsic dynamics of the $i^{th}$ node in $l^{th}$-layer is described by,
$F_{x}(x^{l}_{i},y^{l}_{i})=x_{i}^{l}-x^{3l}_{i}-y^{l}_{i}+I,$ and
$F_{y}(x^{l}_{i},y^{l}_{i},r^{l}_{i})=r^{l}_{i}x^{l}_{i}-by^{l}_{i}$. The
intra-layer is nonlocal and the coupling function is represented by
$H(x^{l}_{j})=\frac{1}{1+e^{-10(x^{l}_{j}+0.25)}}$ as typically used for
defining chemical synapses in neurons. One-to-one inter-layer interactions are
bidirectional diffusive type, $Q(x_{i}^{(n,l)}-x_{i}^{(l,n)})$. The synaptic
coupling $H(x^{l}_{j})$ and bidirectional $Q(x^{(n,l)}_{i})$ are called as
self-coupling functions since they involve the similar variables $x^{n,l}$ and
added to the evolution equation of the same variables $x^{l}_{i}$. The
constants $\varepsilon_{1}$ and $\varepsilon_{2}$ represent their respective
coupling strengths. The self-coupling establishes locally stable CS in the
multiplex network when $\varepsilon_{1}$ and $\varepsilon_{2}$ are larger than
their critical values as determined by the MSF, see Figure 1b. On the other
hand, the coupling function involving $y^{n,l}$ variables is added to the
temporal evolution of the other state variable $x^{l}_{i}$ with a coupling
strength $\kappa$ and henceforth called as cross-coupling links.
$C^{(n,l)}_{i,i}$ represents the elements of the connectivity matrix of the
directed cross-coupling links from one layer to another layer. The cross-
coupling establishes a globally stable synchrony provided the LFS conditions
are satisfied. The multiplex network thereby becomes resilient against
parameter drifting. The choice of $C^{(n,l)}_{ii}$ is system specific, which
follows a set of systematic rules [54] [see Supplementary materials for
further details]. The parameters of the FHN system for all the nodes are
chosen as $I=0.4$, $b=0.1$, $V_{s}=2$. For identical case $r^{(l)}_{i}=r=1$.
To introduce heterogeneity, the parameter $r^{(l)}_{i}$ is perturbed from
$r=1$ to any other positive value, $r^{(l)}_{i}>1$; the perturbation is
randomly introduced within a range $r^{(1)}_{i}\in[1.5,4]$ as elaborated with
an example in Sec. 2.2. Results are also true for $r^{(l)}_{i}<1$, by the same
logic, hence not elaborated here.

For the realization of CS in an identical multiplex network of the FHN
systems, we first derive the critical values of $\varepsilon_{1}$ and
$\varepsilon_{2}$ (when no cross-coupling is present, $\kappa=0$) using the
MSF (see Supplementary material, Sec. 2). Information on stable CS is captured
on a $\varepsilon_{1}$-$\varepsilon_{2}$ plane in Figure 1(b) where the onset
of CS is marked by a boundary line (white dashed line) when the largest
Lyapunov exponent ($LLE$) crosses this line to a negative value. This dashed
line defines the critical self-coupling strengths ($\varepsilon_{1c}$ and
$\varepsilon_{2c}$): Right-hand side of the white dashed line ($LLE=0$) gives
us the range of $\varepsilon_{1}$ and $\varepsilon_{2}$ for the CS regime
($LLE<0$). The color bar indicates the finite values of $LLE$. Note that this
CS state based on MSF is locally stable (valid for small perturbation) when
the state variables $x_{i}^{l}$ and $y_{i}^{l}$ in all the nodes in Layer 1
and Layer 2 oscillate with identical amplitude and phase.

(a) (b)

Figure 1: (a) Schematic diagram of a two-layered multiplex network of FHN
systems. N=50 for both layers (Layer 1 and Layer 2). Each node in both layers
is connected to five neighbors ($P=5$) on both sides (black lines) via
chemical synaptic self-coupling function. Layer 1 interacts by one-to-one
connection to the immediate neighbouring node (dashed blue line) in Layer 2
via diffusive type self-coupling. (b) Phase diagram of $LLE$ in a
$\varepsilon_{1}-\varepsilon_{2}$ plane for identical oscillators when cross-
coupling links are absent ($\kappa=0$). Color bar shows the finite values of
$LLE$, where a white dashed line marks the $LLE=0$ line. The region on the
right hand side of the $LLE=0$ line belongs to locally stable CS state where
$LLE<0$. The system parameters are selected as $I=0.4$, $b=0.1$, $V_{s}=2$.
For identical case $r^{(1,2)}_{i}=r=1$, $i=1,\cdots,N$.

For a preliminary confirmation of the synchronous state, in the temporal
dynamics, we use an error function $e(t)$, considering only $x_{i}$ variable
of each node,

$\displaystyle e(t)=$
$\displaystyle\;\frac{1}{2N}\sqrt{\sum\limits_{i=1}^{2N}\left[x_{i}(t)-\bar{x}(t)\right]^{2}},$
(2)

where, $\bar{x}=\frac{1}{2N}\sum_{i=1}^{2N}x_{i}$ where $i=1,2,3,...2N$
includes all the nodes in two layers. We demonstrate the status of synchrony
with time evolution of the state variables of all the nodes and plots of
synchronization submanifold in the next section.

###  2.1 Generalized synchrony: Single node parameter mismatch

A stable CS state is realized in the identical multiplex network for a choice
of $\varepsilon_{1}>\varepsilon_{1_{c}}$=0.011 and
$\varepsilon_{2}>\varepsilon_{2_{c}}$=0.01 above the critical value as
estimated by the MSF; the broad range of $\varepsilon_{1}$ and
$\varepsilon_{2}$ is presented in Figure 1b. A plot of the error function
$e(t)$ in Figure 2 captures the collective temporal behavior of the network
before any perturbation of the nodes, which decays to zero after a transient
period. At an instant $t=200$ (marked by a vertical solid line), the parameter
$r^{(1)}_{2}$ of node-2 in Layer-1 is detuned manually keeping rest of the
oscillators identical and $\kappa=0$, which destabilizes the synchronous state
when $e(t)\neq 0$ is seen oscillatory. At the time $t\geq 400$ (vertical red
dashed line), one selective, directed cross-coupling link is added from node-2
in Layer-2 to the perturbed node-2 in Layer-1 with a coupling strength
$\kappa=-1$, resulting in restoration of synchrony of all the nodes of the
multiplex network after a transient time when $e(t)=0$. The strength of the
additional cross-coupling link ($\kappa=-1$) is analytically obtained by the
LFS analysis of two coupled FHN systems (node-2 of Layer-1 and Layer-2), and
it is found valid for the multiplex network (see sec. 2.1.1 for analytical
calculations) and confirmed by numerical results. The procedure to select the
particular type of the cross-coupling link is presented in the Supplementary
material. Now a question arises: What type of synchrony emerges in the
multiplex network? The state of synchrony of the multiplex network is
scrutinized using plots of the state variables of the perturbed node against
the respective state variables of unperturbed nodes in Figure 3. The network
becomes resilient to parameter drifting in the sense that any change in
$r^{(1)}_{2}$ cannot break the synchrony of the multiplex network.

Figure 2: Time evolution of the error function. In identical case
($r^{(1,2)}_{i}=r=1$), the error function becomes $e(t)$=0 after a transient
time, prior to any mismatch ($r^{(1)}_{2}=2$) introduced in node-2 of layer-1
at $t=200$ (solid vertical line), resulting in a loss of synchrony as
indicated by oscillations in $e(t)$ in the time interval $t=200$ to $400$, in
absence of cross-coupling ($\kappa=0$). The synchrony is re-established by the
addition of an unidirectional cross-coupling link ($\kappa=-1$) from node-2 in
Layer-2 to node-2 of Layer-1 at $t\geq 400$ (vertical red dashed line). The
cross-coupling link restores synchrony in the perturbed network. Self-coupling
strengths are $\varepsilon_{1}=0.1$ and $\varepsilon_{2}=0.1$, other
parameters are $I=0.4$, $b=0.1$, $V_{s}=2$. Figure 3: Intra- and inter-layer
synchrony in the multiplex network. Plot of $y_{1}^{(1)}$vs. $y_{i}^{(1)}$ for
intra-layer synchrony in (a, c, e) and $y_{i}^{(1)}$vs. $y_{i}^{(2)}$ for
inter-layer synchrony in (b, d, f). CS state (a, b) is realized in two layers
for $\varepsilon_{1}$=$\varepsilon_{2}$=$0.1$, determined by the MSF, in
absence of cross-coupling $\kappa$=$0$. (c, d) Loss of synchrony is observed
due to heterogeneity in node-2 induced by detuning $r^{(1)}_{2}$=$2$, where
$r^{(1)}_{i}$=$r^{(2)}_{i}$=$r$=$1$ for all other nodes, when $\kappa$=$0$.
(e, f) Synchrony is restored by the addition of a cross-coupling link
($\kappa$=$-1$) from node-2 of layer-2 to the perturbed node-2 of layer-1.
Addition of cross-coupling link enhances MSF based local stability of
synchrony to globally stable synchrony as established by the LFS conditions as
detailed in the main text. Node-2 of layer-1 transits to a GS state as
indicated by its rotation (red line) from the CS state (black lines)
maintained by all other nodes. Amplitude scaling of $y_{1}^{(1)}$ is indicated
by the enlargement of the GS manifold (red line). (g ,h) Actual nature of
synchrony of the perturbed node against all other nodes are shown by plotting
$y_{2}^{(1)}$vs. $2y_{i}^{(1)}$ and $y_{2}^{(1)}$vs. $2y_{i}^{(2)}$, which
also confirms the amplification by a factor of parameter mismatch ratio
($a$=$\frac{r^{(1)}_{2}}{r}=2$).

The status of intra-layer synchrony is shown in Figures 3(a, c, e), where
$y_{1}^{(1)}$ of node-1 ($i=1$) of Layer-1 is plotted against $y_{i}^{(1)}$ of
all the nodes in the same layer. The inter-layer synchrony is shown in Figures
3(b, d, f), where one-to-one correlation is seen between $y_{i}^{(1)}$ of all
nodes in Layer-1 against $y_{i}^{(2)}$ of all nodes in Layer-2. For identical
oscillators, CS is thus realized among all the intra- and inter-layer nodes as
shown in Figures 3(a) and 3(b), respectively, for $\kappa$=0 and when
$\varepsilon_{1}=0.1$ and $\varepsilon_{2}=0.1$, which are chosen larger than
their critical values as defined by the MSF plot in Fig. 1(b). Both the intra-
and inter-layer CS are destroyed by induced heterogeneity in node-2 (by our
choice) of Layer-1 as clearly seen in Figures 3(c) and 3(d), respectively. The
synchrony is then restored by adding a selective cross-coupling link as shown
in Figures 3(e) and 3(f), where $y_{1}^{(1)}$ vs. $y_{i}^{(1)}$ and
$y_{i}^{(1)}$ vs. $y_{i}^{(2)}$ plots indicate a strong correlation, but
rotate (red lines) away from the submanifold of the CS state (black lines) of
the identical network. This indicates a new emergent synchronous state in the
multiplex network, when the detuned node variable $y_{1}^{(1)}$ is no more in
a CS state with $y_{i}^{(1,2)}$ of all the nodes in two layers. The exact
nature of synchrony of the perturbed node against all other nodes is revealed
by the $y_{2}^{(1)}$vs. $2y_{i}^{(1)}$ and $y_{2}^{(1)}$vs. $2y_{i}^{(2)}$
plots (red lines) in Figures 3(g) and 3(h), respectively, where the detuned
parameter $r_{2}^{(1)}$=$ar_{i}^{(1,2)}$=$2r$ ($a=2$, amplification factor) is
involved in the rotation of the synchronization submanifold. This is
manifested by a rotation of the CS manifold
$y^{(1,2)}_{1}=y^{(2)}_{2}=y^{(1,2)}_{3}=\cdots=y^{(1,2)}_{N}$ (black lines)
to an emergent synchronization submanifold,
$y^{(1)}_{2}=2y^{(2)}_{2}=2y^{(1,2)}_{1}=2y^{(1,2)}_{3}=\cdots=2y^{(1,2)}_{N}$
(red lines). While all the state variables in the network remain completely
phase coherent, the amplitude of $y_{2}^{(1)}$ connected to the perturbed
node-2 is amplified/attenuated by a scaling factor
$a=\frac{r_{2}^{(1)}}{r_{i}^{(1,2)}}=2$, which is the ratio of the detuned
parameter and the original parameter. We define this emergent state of
synchrony as a class of GS state of the whole multiplex network. A most
exciting and novel part is that the state variables of all the intra- or
inter-layer nodes maintain the CS state except the perturbed node. The
synchrony quality is enhanced by adding the particular, directed cross-
coupling link establishing resilience against a parameter perturbation. The
phenomenon remains independent of the size of the perturbation (small or
large). The amount of perturbation cannot destroy GS in the multiplex network
when the amplification/attenuation of the particular state variable of the
perturbed node only varies linearly following the amount of perturbation with
a constant of proportionality $a$. Interestingly, the cross-coupling link does
not affect the CS state or sleeps when all the parameters are identical,
making the cross-coupling function zero yet making CS globally stable. The
cross-coupling starts working or becomes awake when a parameter drift appears
in any node, leading to GS of the network, which is globally stable as
elaborated in the next section.

####  2.1.1 Analytical Results: Cross-coupling strength in a reduced system

For a given multiplex network under parameter mismatch in one node (simplest
case), we only consider the detuned node and its immediate adjacent node from
the other layer, which is in CS with all other nodes in both layers of the
network. The entire multiplex network is then reduced to two representative
nodes, one detuned node, and one unperturbed node, and now we derive the
strength ($\kappa$) of the cross-coupling link using the LFS conditions,
between the two nodes, node-2 of Layer-1 and node-2 of Layer-2. For ease of
calculations, we replace the notations of all the variables of the two nodes
by $x_{2}^{(1)}$=$x_{1}$ and $x_{2}^{(2)}$=$x_{2}$, and $y_{2}^{(1)}$=$y_{1}$
and $y_{2}^{(2)}$=$y_{2}$, and the parameters by $r_{2}^{(1)}$=$r_{1}$ and
$r_{2}^{(2)}$=$r_{2}$. Similarly, discarding the layer indices, we can write
the dynamical equations of the two nodes as,

$\displaystyle\dot{x}_{1}=\;x_{1}-x_{1}^{3}-{y}_{1}+I+\varepsilon(x_{2}-x_{1})+\kappa(y_{2}-y_{1})$
$\displaystyle\dot{y}_{1}=\;r_{1}x_{1}-by_{1}$
$\displaystyle\dot{x}_{2}=\;x_{2}-x_{2}^{3}-{y}_{2}+I+\varepsilon(x_{1}-x_{2})$
$\displaystyle\dot{y}_{2}=\;r_{2}x_{2}-by_{2}$ (3)

where $(x_{1,2},y_{1,2})$ are state variables, $\varepsilon$ is self-coupling
strength. The conventional diffusive coupling function $(x_{2,1}-x_{1,2})$ is
defined as self-coupling. A parameter mismatch is introduced by taking $r_{1}$
and $r_{2}$ parameters as different. A directed cross-coupling link
$(y_{1}-y_{2})$ (from one node to the other) is then added to the dynamics of
$x_{1}$ variable with coupling strength $\kappa$. The choice of both the self-
and cross-coupling are made following systematic generic rules and detailed in
the Supplementary material.

Case I: For identical systems:

The LFS is analyzed first to establish globally stable synchrony in two nodes
for the identical case, i.e. $r_{1}=r_{2}=r$. The error function $\bf
e$=$[e_{x},e_{y}]^{T}=[x_{1}-x_{2},y_{1}-y_{2}]^{T}$ of system Eq. (3) evolves
as,

$\displaystyle\dot{e}_{x}=$
$\displaystyle\;e_{x}-\frac{e_{x}^{3}}{4}-\frac{3}{4}e_{x}e_{p}^{2}-e_{y}-2\varepsilon_{2}e_{x}-\kappa
e_{y}$ $\displaystyle\dot{e}_{y}=$ $\displaystyle\;re_{x}-be_{y}$ (4)

where, $e_{p}=x_{1}+x_{2}$ so that
$x_{1}^{3}-x_{2}^{3}=\frac{e_{x}}{4}(e_{x}^{2}+3e_{p}^{2})$. For a global
stability of $(e_{x}=0,e_{y}=0)$, we define a Lyapunov function,
$V(e)=\frac{1}{2}e_{x}^{2}+\frac{1}{2}e_{y}^{2}$. We first check the stability
of $x_{1}=x_{2}$, separately, by defining a partial Lyapunov function,
$V^{\prime}(e_{x})=\frac{1}{2}e_{x}^{2}$ when its time derivative is

$\displaystyle\dot{V^{\prime}}(e_{x})=$ $\displaystyle-
e_{x}^{2}\left(\frac{3}{4}e_{p}^{2}-1+2\varepsilon_{2}\right)-\frac{e_{x}^{4}}{4}-(1+\kappa)e_{y}e_{x}.$
(5)

$\dot{V}^{\prime}$($e_{x})<0$ provided $\kappa=-1$ and
$\frac{3}{4}e_{p}^{2}-1+2\varepsilon_{2}\geq 0$. To satisfy the condition, we
derive the roots of the equation,

$\displaystyle\frac{3}{4}e_{p}^{2}-1+2\varepsilon_{2}=0$ (6) $\displaystyle
e_{p}^{2}=\frac{4}{3}(1-2\varepsilon_{2})$ (7)

the positive root of $e_{p}$ is given by

$\displaystyle{e_{p}}=\sqrt{\frac{4}{3}(1-2\varepsilon_{2})}.$ (8)

${e_{p}}$ will now be positive real if $\frac{4}{3}(1-2\varepsilon_{2})\geq 0$
that implies ${\varepsilon_{2}}\leq 1/2$. Thus $\dot{V^{\prime}}(e_{x})<0$ is
satisfied when

$\displaystyle\kappa=-1,\quad\text{and}\quad\varepsilon_{2}\leq 1/2,$ (9)

this ensures partial stability of $x_{1}=x_{2}$. Substituting the conditions
in (9) into Eq. (4) and assuming identical systems ($r_{1}=r_{2}$), we obtain
$\dot{e}_{y}=-be_{y}$, when

$\dot{V}(e_{x},e_{y})=-\frac{e_{x}^{4}}{4}-be_{y}^{2}<0,\quad\varepsilon_{2}\leq\frac{1}{2}\quad\text{and}\quad\kappa=-1.$
(10)

Once this LFS condition (10) is satisfied, the CS state $x_{1}=x_{2}$ and
$y_{1}=y_{2}$ becomes globally stable in presence of the selective self- and
cross-coupling links.

Case II: For non-identical systems:

Now we study the effect of heterogeneity on stability of CS by detuning one of
the parameters when $r_{1}\neq r_{2}$. The induced heterogeneity parameters
$r_{1}$ and $r_{2}$ are not involved in Eq.(5), hence stability of the
$x_{1}=x_{2}$ is still preserved. After detuning $r_{1,2}$, the equation of
$\dot{e}_{y}$ from Eq. (4) becomes,

$\begin{array}[]{l}\dot{e}_{y}=r_{1}x_{1}-r_{2}x_{2}-be_{y}=(r_{1}-r_{2})x_{1}-be_{y}\\\
\quad\;=\frac{r_{1}-r_{2}}{r_{1}}(\dot{y}_{1}+by_{1})-be_{y},\end{array}$ (11)

and this leads to

$\displaystyle\dot{y}_{1}$
$\displaystyle(1-\frac{r_{1}-r_{2}}{r_{1}})-\dot{y}_{2}=-by_{1}(1-\frac{r_{1}-r_{2}}{r_{1}})+by_{2}$
$\displaystyle\dot{y}_{1}$
$\displaystyle\frac{r_{2}}{r_{1}}-\dot{y}_{2}=-b(y_{1}\frac{r_{2}}{r_{1}}-y_{2}).$
(12)

From Eq.(12), the revised error dynamics is $\dot{e}_{y}^{*}=-be_{y}^{*}$,
where the modified error function becomes

$\displaystyle e_{y}^{*}=y_{1}\frac{r_{2}}{r_{1}}-y_{2}.$ (13)

Rewriting Eqn. (13) in actual terms, we have

$\displaystyle e_{y}^{*}=y_{2}^{(1)}\frac{r}{r_{2}^{(1)}}-y_{2}^{(2)}.$ (14)

Accordingly, the Lyapunov function is redefined in terms of the modified error
functions whose time derivative is

$\dot{V}^{*}(e_{x},e_{y}^{*})=-\frac{e_{x}^{4}}{4}-be_{y}^{*2}<0$ (15)

which ensures a globally stable synchrony, which basically represents an
emergent globally stable GS state as cliamed in the previous section. The
manifestation of the GS state is found in the perturbed node,

$y_{2}^{(1)}\frac{r}{r_{2}^{(1)}}=y_{2}^{(2)},\;\text{or}\;y_{2}^{(1)}=\frac{r_{2}^{(1)}}{r}y_{2}^{(2)}$

while all other state variables of all the nodes in the multiplex network
maintain CS ($x_{1}=x_{2}$ and $y_{1}=y_{2}$).

Under this stability condition, all other identical nodes emerge into CS and
can be treated as if isolated since the self-coupling terms eventually vanish.
However, the cross-coupling link in the perturbed node remains active for any
instant of time and plays a constructive role for both identical and non-
identical cases. Global stability of CS is maintained, in an identical case,
beyond the MSF-based local stability, when the cross-coupling function becomes
effectively zero or sleeping since $y_{2}^{1}$=$y_{2}^{2}$. Under the
perturbed conditions, the cross-coupling function is awake, and the network
emerges into a globally stable GS state while MSF fails to ensure the
stability of synchrony of any form. We find that this particular strength of
$\kappa$ as analytically derived for two nodes here works perfectly well for
an extensive network and for multi-node perturbation to ensure a resilient
multiplex network of FHN systems with an emergent globally stable GS. We must
mention here that the amplification factor ($a=\frac{r_{2}}{r_{1}}$) is a very
general effect for any pair of nodes between the two layers for any topology
of the multiplex network. The question remains, how this particular choice of
cross-coupling link is to be made? The choice of the cross-coupling function
is systematically made by a semi-analytic approach as detailed in [54] and
elaborated briefly in the Supplementary material. For a particular choice of
node dynamics, an appropriate cross-coupling link can always be found to
realize globally stable GS in perturbed multiplex networks of various
topologies. Most importantly, the original topology of the multiplex network
remains unchanged. Only one sleeping directed cross-coupling link is added for
each perturbed node to the identical network that awakes only when a parameter
drift appears. The next section presents numerical results on perturbation in
multiple nodes in a layer.

###  2.2 State of coherence: Heterogeneity in a layer

We numerically demonstrate the efficacy of the cross-coupling strategy in the
case of multiple perturbed nodes. All the nodes of one layer (say, Layer-1)
are made heterogeneous by a normal distribution of $r^{(1)}_{i}\in(1.5,4]$
with a mean 3 and standard deviation 0.5 where $i=1,2,\cdots,50$. This range
of $r_{i}$ and its distribution are arbitrarily chosen.

Figure 4: (Color online) Synchrony in the multiplex network for distributed
heterogeneity in one layer. Parameters $r^{(1)}_{i}$ in layer-1 are detuned
randomly within a range of $r^{(1)}_{i}\in(1.5,4)$ with a normal distribution
of mean 3 and standard deviation 0.5; for unperturbed nodes ($i=51$ to $100$)
in Layer-2, $r_{i}^{(2)}$=$r$=1. (a) Snapshot of temporal and (b) spatio-
temporal dynamics of $x_{i}$ of all the nodes, show CS in amplitude and phase.
(c) Snapshot of temporal evolution and (d) spatio-temporal plots of
$y_{i}^{(1,2)}$ confirms amplified response all perturbed nodes in Layer-1
($i=1$ to $50$) while all the identical nodes ($i=51$ to $100$) in Layer-2
maintain CS. (e) Time evolution of $y_{i}^{(1)}$ corresponding to the
perturbed nodes in Layer-1 are shown in various colors, which shows phase
coherent evolution with varying amplitude. A plot of $y_{i}^{(2)}$ of all the
unperturbed identical nodes in Layer-2 is shown (red line) for comparison and
amplitude response of the perturbed nodes. The scaling up of amplitude in
$y_{i}^{(1)}$ is due to distributed heterogeneity in the parameter
$r^{(1)}_{i}$. The amplified response in the perturbed nodes is dictated by
the heterogeneity of the nodes (various colors) quantified by the ratio of the
parameter mismatch, $r^{(1)}_{i}/r$. (f) Mean phase velocity ($\Omega_{j}$)
confirms complete phase coherence in all the nodes in the network. The
frequency of all the nodes remains identical and is preserved against
parameter drifting.

The parameter values $r^{(2)}_{i}$ of all the nodes in Layer-2 are assumed
identical i.e., $r^{(2)}_{i}=r=1$, other system parameters are kept identical
too. Unidirectional cross-coupling links from Layer-2 to all adjacent
perturbed nodes in Layer-1 are added, keeping other coupling functions and the
original topology of the network unchanged. The selection process of the
coupling profile (self- and cross-coupling) remains the same as done for
single node perturbation. The self-coupling strengths are chosen once again as
$\varepsilon_{1}=0.1$ and $\varepsilon_{2}=0.1$ larger than the critical
values defined by the MSF. The strength of the cross-coupling links is chosen
$\kappa=-1$, as derived by the LFS condition for two coupled systems as shown
in sec. 2.1.1. We show that our strategy successfully works for perturbation
in multiple nodes in the network.

The emergent behavior of the network under perturbation is described in Figure
4. The CS state in $x_{i}$ variable of both the layers is not disturbed as
shown by a snapshot of $x_{i}^{(1,2)}$ and their spatio-temporal dynamics in
Figs. 4(a) and 4(b), respectively. The amplitude of the $y_{i}$-variables of
all the nodes in Layer-1 are amplified by the ratio,
$a_{i}=\left(\frac{r^{(1)}_{i}}{r}\right)$, where $a_{i}$ follows the
distribution of $r_{i}^{(1)}$ of the perturbed nodes and $r=r^{(2)}_{i}=1$,
$i=1,2,\cdots,50$. A snapshot of $y^{(1,2)}_{i}$-variables in Figure 4(c)
shows a pattern of distributed amplitude in the nodes of Layer-1 (1 to 50),
while the rest of the nodes (51 to 100) of Layer-2 have identical amplitude.
The distribution of amplitude in the first 50 nodes reflects the distribution
of $a_{i}$ due to distribution in $r_{i}^{(1)}$ in layer-1. The spatio-
temporal dynamics in Figure 4(d) corroborates the fact that the first 50 nodes
in Layer-1 are phase coherent but with a distribution of amplitude (color bar
indicates the distribution in amplitude) while the rest of the 50 nodes in
Layer-2 are both amplitude and phase coherent (coherent pattern in cyan
color). The amplitude response of the nodes in Layer-1 reveals a distortion-
free amplification in the heterogeneous nodes, as confirmed by the temporal
dynamics of $y_{i}^{(1,2)}$ (color lines) in Figure 4(e). To estimate the
local phase correlation of all the nodes in Layer-1 and layer-2, a long time
average of phase velocity [68] of $j^{th}$ oscillators is taken,

$\displaystyle\Omega_{j}=\frac{2\pi M_{j}}{\Delta
t},~{}~{}~{}~{}j=1,2,\cdots,2N,$ (16)

where $M_{j}$ is the average number of periods of the $j^{th}$ node for a long
time interval $\Delta t$. All the nodes in the two layers are seen completely
phase coherent as seen in Figure 4(f), where they all have identical mean
frequency $\Omega_{j}$. It is an essential requirement for many real-world
multiplex networks, where any frequency drift is prevented by these additional
cross-coupling links against parameter perturbation in any node. The design of
cross-coupling links follows a systematic semi-analytical technique as usual,
and hence it is quite a novel proposition for realizing a resilient multiplex
network.

##  3 Random multiplex network: Rössler system

Now we confirm that our coupling strategy is independent of both the
connectivity matrices of the layers of the multiplex network and their nodal
dynamics. A random multiplex network of two layers is considered with nodal
dynamics of a chaotic Rössler oscillator and a distributed parameter in one
layer. The governing dynamics of the $i^{th}$ node,

$\displaystyle\dot{x}_{i}^{(1)}=-y^{(1)}_{i}-z^{(1)}_{i}+\kappa(z^{(2)}_{i}-z^{(1)}_{i})$
$\displaystyle\dot{x}_{i}^{(2)}=-y^{(2)}_{i}-z^{(2)}_{i}$
$\displaystyle\dot{y}_{i}^{(1,2)}=\;x^{(1,2)}_{i}+ay^{(1,2)}_{i}+\varepsilon_{1}\sum\limits_{j=1}^{N}A^{(1,2)}_{ij}(y^{(1,2)}_{j}-y^{(1,2)}_{i})$
$\displaystyle\quad\quad\quad+\varepsilon_{2}(y^{(1,2)}_{i}-y^{(2,1)}_{i})$
$\displaystyle\dot{z}_{i}^{(1,2)}=\;b^{(1,2)}_{i}-cz^{(1,2)}_{i}+x^{(1,2)}_{i}z^{(1,2)}_{i},$
(17)

Figure 5: (Color online) Multiplex network of Rössler oscillators. Node
indices 1-50 and 51-100, respectively, belong to Layer-1 and Layer-2. Self-
coupling strengths $\varepsilon_{1}=1$ and $\varepsilon_{2}=1$, and the
strength of cross-coupling link is $\kappa=-1$ [54, 55]. For a choice of
parameters $a=0.2$, $c=4.8$ and $b^{(2)}_{i}=b=0.2$, nodes in isolation
exhibit chaotic dynamics. We perturb layer-1 by taking a distribution in
$b^{(1)}_{i}\in(0.4,0.8)$, $i=1,\cdots,50$. The entire network perfectly
synchronized even under the induced heterogeneity due to the presence of
cross-coupling links from Layer-2 to Layer-1, when the $x$-variable shows
complete coherence in amplitude and phase as shown in (a) and (b). The
amplitude of $z$-variables in layer-1 are linearly scaled up by the amount of
parameter mismatches, although the phase coherence is not disturbed as
confirmed by a snap shot of $Z_{i}^{1,2}$ of all the nodes in two layers (c)
and a spatio-temporal plot (d). (e) Time evolution of $Z_{i}^{1,2}$ shows
amplified replica (multiple colors) of all the perturbed nodes against the
identical nodes (red line), when only the amplitude is scaled up with complete
phase coherence. (f) $Z_{i}^{1}$ vs $Z_{i}^{2}$ plot. Synchronization
manifolds of the detuned nodes are rotated along the transverse direction
against the CS manifold (red line) of the identical nodes. The rotation of the
individual plot depends upon the ratio of mismatched parameters to the
original (unperturbed) parameter.

where $x^{(1,2)}_{i},y^{(1,2)}_{i}$ and $z^{(1,2)}_{i}$ are the state
variables; superscripts $1$ and $2$ denote Layer-1 and Layer-2, respectively.
$A^{(1,2)}_{ij}$ is the adjacency matrix of an individual layer. We use a
rewiring strategy developed by Watts and Strogatz [47] to construct a random
network for each layer independently, hence making two different network
structures for the two layers where the inter-layer connections are one-to-
one. The randomness of the two layers is defined by their degree distributions
($\in$ [1,8]). The average degrees of Layer-1 and layer-2 are 6 and 5,
respectively.

The intra-layer coupling for both the layers is defined by a diffusive self-
coupling function $(y^{(1,2)}_{j}-y^{(1,2)}_{i})$ with a coupling strength
$\varepsilon_{1}$. One-to-one diffusive self-coupling links
$(y^{(1,2)}_{i}-y^{(2,1)}_{i})$ establishes the inter-layer connectivity
between the adjacent nodes of the two layers with a strength
$\varepsilon_{2}$. An isolated node exhibits chaotic dynamics by the choice of
parameters as $a=0.2,c=4.8$ and $b^{(1,2)}_{i}=b=0.2$ when all the nodes are
identical. The identical multiplex network maintains CS with local stability
for coupling strength $\epsilon_{1}$ and $\epsilon_{2}$ decided by the MSF
condition. Next, we introduce heterogeneity in all the nodes of Layer-1 with
uniform distribution in $b^{(1)}_{i}$ in a range $b^{(1)}_{i}\in(0.4,0.8)$
while all the nodes in Layer-2 are identical $b^{(2)}_{i}=b=0.2$, thereby
destroying CS state in the multiplex network. Then, we apply our strategy of
selective addition of directed cross-coupling links to each perturbed node of
Layer-1 from the adjacent nodes of Layer-2. The directed cross-couplings are
linear and diffusive as defined by $(z^{(2)}_{i}-z^{(1)}_{i})$ and added to
the dynamics of the $x_{i}^{(1)}$-variable with a strength $\kappa$. The
choice of the coupling links is made by our proposed systematic rules [54],
and the strength of $\kappa$ is defined by the LFS condition [see
Supplementary material [supple] for details]. Globally stable synchrony of the
whole multiplex network is then reestablished. Figure 5 presents numerical
results.

Figures for CS in the identical multiplex network are not presented here. The
$x_{i}^{(1,2)}$-variable for the two layers in the perturbed network shows CS
in amplitude and phase in a snapshot in Figure 5(a) and a spatio-temporal plot
in Figure 5(b). It is also found true for $y_{i}^{(1,2)}$ dynamics in all the
nodes in two layers in the perturbed multiplex network (not shown here).
However, the amplitudes $z_{i}$ of 50 nodes in Layer-1 are seen distributed in
Figure 5(c) while the other 50 nodes of layer-2 remain in CS. The spatio-
temporal plot in Figure 5(d) shows phase coherence among all the nodes.
Looking at the temporal dynamics of the nodes in both the layers in Figure
5(e), it is clear that the amplitudes of the perturbed 50 nodes in Layer-1 are
linearly scaled up by $\frac{r^{(1)}_{i}}{r}$, where $r_{i}^{(1)}$ follows a
distribution, which is reflected in the amplitude distribution of the first
$50$ nodes of Layer-1 in Figure 5(c). Due to induced distributed
heterogeneity, the GS manifold of the detuned nodes, $z_{i}^{(2)}$ vs.
$z_{i}^{(2)}$, in Figure 5(f), shows a variation in rotation angle of
individual plots along the transverse direction of the CS manifold (red line).
However, the phase coherence is not disturbed as usual by the distribution of
parameter $b_{i}^{(1)}$ as revealed in Figures 5(e)-(f).

##  4 Conclusion

An essential question in multiplex networks is how to prevent loss of
synchrony against parameter perturbation (or drifting)? We addressed the
question with two examples of two-layered multiplex networks with a suggestion
of adding selective linear diffusive cross-coupling links to the perturbed
nodes in one layer from the adjacent nodes of the unperturbed layer without
deleting any link from the original network. The selection of the coupling
functions (self- and cross-coupling links) is most important for the multiplex
network that follows a systematic rule formulated from the LFM of the dynamics
of the nodes and an analytical procedure based on the LFS measure. This
selection procedure is primarily elaborated in our earlier work [54] and
briefly explained in the Supplementary material.

An identical multiplex network with self-coupling maintains locally stable CS
following the MSF conditions. Under parameter perturbation in any node of one
layer, CS is broken when selective cross-coupling links restore synchrony as a
globally stable GS state in the multilayer network. More categorically, all
the nodes attain a globally stable CS except the perturbed nodes that emerge
into a GS state. The amplitude of some of the state variables of the perturbed
nodes shows a linear amplitude response to the amount of parameter
perturbation in the GS state; the perturbed nodes remain phase coherent
preventing a frequency drift, and it is independent of the size of the
perturbation. This approach enables us to design resilient multiplex networks.

To demonstrate the efficacy of our proposal, we consider two exemplary
multiplex networks with different topologies and different dynamical systems.
As a first example, we use a multiplex network of two layers with the slow-
fast FHN system representing the nodal dynamics and consider unidirectional
chemical synapses that interact between the nodes of both layers, representing
the nonlinear coupling between the nodes. We use the nonlocal coupling
topology of the individual layers. CS is established with an assumption of
identical nodes and deriving the critical coupling strength using the MSF
formalism in the multiplex network. One of the parameters in one node in a
layer is then perturbed, and we define the appropriate cross-coupling function
as decided from the FHN dynamics and add to the perturbed node from the
unperturbed layer. A globally stable GS state is realized when the stability
conditions are derived using the LFS measure.

Most importantly, none of the original links of the network is deleted in
contrast to what usually happens during the rewiring of links as used by
others for enhancing synchrony in networks. In the unperturbed state, the
cross-coupling links sleep without affecting the original dynamics and the CS
state and awake or become active when a drift in parameters of the nodes
appears. We successfully tested the results as true for induced heterogeneity
in all the nodes in one layer when one cross-coupling link is to be added to
each node from the immediate adjacent node of the other layer. Such a design
of cross-coupling links can be adopted a priori for realizing a resilient
multiplex network. We confirm our results with the second example of a
multiplex network, where two layers of random networks with multiplexing are
used, and a chaotic Rössler oscillator represents the nodal dynamics. Our
proposed framework of additional cross-coupling links is thus independent of
the dynamics of the nodes and the topology of the layers in a multiplex
network.

## Additional information

Selection of coupling profile from the linear flow of FitzHugh Nagumo and
Rössler system are provided in the Supplementary material, Secs. 1.1 and 1.2,
respectively. The stability conditions of CS of the FHN system (1) in the
multiplex network are derived using the MSF in the supplementary material,
Sec. 2. For Matlab codes, please visit, https://github.com/ecesuman06/cross-
coupling-multiplex.

## Acknowledgments

My immense gratitude to Dr. Syamal K. Dana for fruitful discussions,
suggestions and help preparing the article. I also acknowledge Dr. Ritwika
Mondal for fruitful discussions and help in MSF derivation. I am very grateful
to the anonymous reviewers for their critical comments and suggestion that
helped me improve the quality of the work.

## References

  * [1] S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, “The structure and dynamics of multilayer networks,” Physics Reports, vol. 544, no. 1, pp. 1–122, 2014.
  * [2] D. Y. Kenett, M. Perc, and S. Boccaletti, “Networks of networks–an introduction,” Chaos, Solitons & Fractals, vol. 80, pp. 1–6, 2015.
  * [3] Y. Moreno and M. Perc, “Focus on multilayer networks,” New Journal of Physics, vol. 22, no. 1, p. 010201, 2019.
  * [4] R. L. Breiger and P. E. Pattison, “Cumulated social roles: The duality of persons and their algebras,” Social Networks, vol. 8, no. 3, pp. 215–256, 1986.
  * [5] M. Jalili, Y. Orouskhani, M. Asgari, N. Alipourfard, and M. Perc, “Link prediction in multiplex online social networks,” Royal Society Open Science, vol. 4, no. 2, p. 160863, 2017.
  * [6] A. Cardillo, M. Zanin, J. Gómez-Gardenes, M. Romance, A. J. G. del Amo, and S. Boccaletti, “Modeling the multi-layer nature of the european air transport network: Resilience and passengers re-scheduling under random failures,” The European Physical Journal Special Topics, vol. 215, no. 1, pp. 23–33, 2013.
  * [7] J. Wu, C. Pu, L. Li, and G. Cao, “Traffic dynamics on multilayer networks,” Digital Communications and Networks, vol. 6, no. 1, pp. 58–63, 2020.
  * [8] A. Halu, M. De Domenico, A. Arenas, and A. Sharma, “The multiplex network of human diseases,” NPJ Systems Biology and Applications, vol. 5, no. 1, pp. 1–12, 2019.
  * [9] S. Pilosof, M. A. Porter, M. Pascual, and S. Kéfi, “The multilayer nature of ecological networks,” Nature Ecology & Evolution, vol. 1, no. 4, pp. 1–9, 2017.
  * [10] L. Liang, V. Chen, K. Zhu, X. Fan, X. Lu, and S. Lu, “Integrating data and knowledge to identify functional modules of genes: a multilayer approach,” BMC Bioinformatics, vol. 20, no. 1, pp. 1–15, 2019.
  * [11] M. Zitnik and J. Leskovec, “Predicting multicellular function through multi-layer tissue networks,” Bioinformatics, vol. 33, no. 14, pp. i190–i198, 2017.
  * [12] N. Timme, S. Ito, M. Myroshnychenko, F.-C. Yeh, E. Hiolski, P. Hottowy, and J. M. Beggs, “Multiplex networks of cortical and hippocampal neurons revealed at different timescales,” PloS one, vol. 9, no. 12, p. e115764, 2014.
  * [13] D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson, and S. T. Grafton, “Dynamic reconfiguration of human brain networks during learning,” Proceedings of the National Academy of Sciences, vol. 108, no. 18, pp. 7641–7646, 2011.
  * [14] J. J. Crofts, M. Forrester, and R. D. O’Dea, “Structure-function clustering in multiplex brain networks,” EPL (Europhysics Letters), vol. 116, no. 1, p. 18003, 2016.
  * [15] P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, “Community structure in time-dependent, multiscale, and multiplex networks,” Science, vol. 328, no. 5980, pp. 876–878, 2010.
  * [16] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, “Multilayer networks,” Journal of Complex Networks, vol. 2, no. 3, pp. 203–271, 2014.
  * [17] F. Sorrentino, “Synchronization of hypernetworks of coupled dynamical systems,” New Journal of Physics, vol. 14, no. 3, p. 033035, 2012.
  * [18] R. Sevilla-Escoboza, I. Sendiña-Nadal, I. Leyva, R. Gutiérrez, J. Buldú, and S. Boccaletti, “Inter-layer synchronization in multiplex networks of identical layers,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 26, no. 6, p. 065304, 2016.
  * [19] V. Nicosia, P. S. Skardal, A. Arenas, and V. Latora, “Collective phenomena emerging from the interactions between dynamical processes in multiplex networks,” Physical Review Letters, vol. 118, no. 13, p. 138302, 2017.
  * [20] I. Leyva, I. Sendiña-Nadal, R. Sevilla-Escoboza, V. Vera-Avila, P. Chholak, and S. Boccaletti, “Relay synchronization in multiplex networks,” Scientific Reports, vol. 8, no. 1, p. 8629, 2018.
  * [21] A. Bergner, M. Frasca, G. Sciuto, A. Buscarino, E. J. Ngamga, L. Fortuna, and J. Kurths, “Remote synchronization in star networks,” Physical Review E, vol. 85, no. 2, p. 026208, 2012.
  * [22] S. Rakshit, B. K. Bera, J. Kurths, and D. Ghosh, “Enhancing synchrony in multiplex network due to rewiring frequency,” Proceedings of the Royal Society A, vol. 475, no. 2230, p. 20190460, 2019.
  * [23] V. A. Maksimenko, V. V. Makarov, B. K. Bera, D. Ghosh, S. K. Dana, M. V. Goremyko, N. S. Frolov, A. A. Koronovskii, and A. E. Hramov, “Excitation and suppression of chimera states by multiplexing,” Physical Review E, vol. 94, no. 5, p. 052205, 2016.
  * [24] S. Majhi, M. Perc, and D. Ghosh, “Chimera states in uncoupled neurons induced by a multilayer structure,” Scientific Reports, vol. 6, no. 1, pp. 1–11, 2016.
  * [25] P. Khanra and P. Pal, “Explosive synchronization in multilayer networks through partial adaptation,” Chaos, Solitons & Fractals, vol. 143, p. 110621, 2021.
  * [26] S. Jalan, A. Kumar, and I. Leyva, “Explosive synchronization in frequency displaced multiplex networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 29, no. 4, p. 041102, 2019.
  * [27] A. Babloyantz and A. Destexhe, “Low-dimensional chaos in an instance of epilepsy,” Proceedings of the National Academy of Sciences, vol. 83, no. 10, pp. 3513–3517, 1986.
  * [28] W. W. Lytton, “Computer modelling of epilepsy,” Nature Reviews Neuroscience, vol. 9, no. 8, p. 626, 2008.
  * [29] J. Machowski, J. W. Bialek, and J. Bumby, Power system dynamics: stability and control.  John Wiley & Sons, 2011.
  * [30] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscillators,” Physical Review Letters, vol. 76, no. 11, p. 1804, 1996.
  * [31] S. Rakshit, B. K. Bera, and D. Ghosh, “Synchronization in a temporal multiplex neuronal hypernetwork,” Physical Review E, vol. 98, no. 3, p. 032305, 2018\.
  * [32] L. M. Pecora and T. L. Carroll, “Master stability functions for synchronized coupled systems,” Physical Review Letters, vol. 80, no. 10, p. 2109, 1998\.
  * [33] L. M. Pecora, T. L. Carroll, G. A. Johnson, D. J. Mar, and J. F. Heagy, “Fundamentals of synchronization in chaotic systems, concepts, and applications,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 7, no. 4, pp. 520–543, 1997.
  * [34] P. K. Roy, S. Chakraborty, and S. K. Dana, “Experimental observation on the effect of coupling on different synchronization phenomena in coupled nonidentical chua’s oscillators,” Chaos: An interdisciplinary Journal of Nonlinear Science, vol. 13, no. 1, pp. 342–355, 2003.
  * [35] R. E. Mirollo and S. H. Strogatz, “Amplitude death in an array of limit-cycle oscillators,” Journal of Statistical Physics, vol. 60, no. 1-2, pp. 245–262, 1990.
  * [36] G. Saxena, A. Prasad, and R. Ramaswamy, “Amplitude death: The emergence of stationarity in coupled nonlinear systems,” Physics Reports, vol. 521, no. 5, pp. 205–228, 2012.
  * [37] C. Hens, P. Pal, S. K. Bhowmick, P. K. Roy, A. Sen, and S. K. Dana, “Diverse routes of transition from amplitude to oscillation death in coupled oscillators under additional repulsive links,” Physical Review E, vol. 89, no. 3, p. 032901, 2014.
  * [38] M. Nandan, C. Hens, P. Pal, and S. K. Dana, “Transition from amplitude to oscillation death in a network of oscillators,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 24, no. 4, p. 043103, 2014\.
  * [39] C. J. Tessone, C. R. Mirasso, R. Toral, and J. D. Gunton, “Diversity-induced resonance,” Physical Review Letters, vol. 97, no. 19, p. 194101, 2006.
  * [40] H. Daido, “Aging transition and disorder-induced coherence in locally coupled oscillators,” Europhysics Letters, vol. 84, no. 1, p. 10002, 2008.
  * [41] E. Padmanaban, S. Saha, M. Vigneshwaran, and S. K. Dana, “Amplified response in coupled chaotic oscillators by induced heterogeneity,” Physical Review E, vol. 92, no. 6, p. 062916, 2015.
  * [42] J. Gao, B. Barzel, and A.-L. Barabási, “Universal resilience patterns in complex networks,” Nature, vol. 530, no. 7590, pp. 307–312, 2016.
  * [43] Y. Zhang, T. Nishikawa, and A. E. Motter, “Asymmetry-induced synchronization in oscillator networks,” Physical Review E, vol. 95, no. 6, p. 062215, 2017\.
  * [44] Y. Sugitani, Y. Zhang, and A. E. Motter, “Synchronizing chaos with imperfections,” Physical Review Letters, vol. 126, no. 16, p. 164101, 2021\.
  * [45] E. S. Medeiros, U. Feudel, and A. Zakharova, “Asymmetry-induced order in multilayer networks,” Physical Review E, vol. 104, no. 1, p. 024302, 2021\.
  * [46] R. Banerjee, D. Ghosh, E. Padmanaban, R. Ramaswamy, L. Pecora, and S. K. Dana, “Enhancing synchrony in chaotic oscillators by dynamic relaying,” Physical Review E, vol. 85, no. 2, p. 027201, 2012.
  * [47] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’networks,” Nature, vol. 393, no. 6684, pp. 440–442, 1998\.
  * [48] I. Leyva, R. Sevilla-Escoboza, I. Sendiña-Nadal, R. Gutiérrez, J. Buldú, and S. Boccaletti, “Inter-layer synchronization in non-identical multi-layer networks,” Scientific Reports, vol. 7, no. 1, pp. 1–9, 2017.
  * [49] P. Schultz, T. Peron, D. Eroglu, T. Stemler, G. M. R. Ávila, F. A. Rodrigues, and J. Kurths, “Tweaking synchronization by connectivity modifications,” Physical Review E, vol. 93, no. 6, p. 062211, 2016.
  * [50] J. P. Pade and T. Pereira, “Improving network structure can lead to functional failures,” Scientific Reports, vol. 5, p. 9968, 2015.
  * [51] S. Saha, A. Mishra, S. K. Dana, C. Hens, and N. Bairagi, “Infection spreading and recovery in a square lattice,” Physical Review E, vol. 102, no. 5, p. 052307, 2020.
  * [52] S. K. Dwivedi, M. S. Baptista, and S. Jalan, “Optimization of synchronizability in multiplex networks by rewiring one layer,” Physical Review E, vol. 95, no. 4, p. 040301, 2017.
  * [53] S. Rakshit, S. Majhi, B. K. Bera, S. Sinha, and D. Ghosh, “Time-varying multiplex network: Intralayer and interlayer synchronization,” Physical Review E, vol. 96, no. 6, p. 062308, 2017.
  * [54] S. Saha, A. Mishra, E. Padmanaban, S. K. Bhowmick, P. K. Roy, B. Dam, and S. K. Dana, “Coupling conditions for globally stable and robust synchrony of chaotic systems,” Physical Review E, vol. 95, no. 6, p. 062204, 2017.
  * [55] S. Saha, A. Mishra, P. K. Roy, and S. K. Dana, “Cross-coupling plays constructive role on global stability of synchrony in neuronal networks,” Opera Medica et Physiologica, no. 3-4, pp. 93–98, 2017.
  * [56] A. Suzuki, “Recent advances in the cross-coupling reactions of organoboron derivatives with organic electrophiles, 1995–1998,” Journal of Organometallic Chemistry, vol. 576, no. 1-2, pp. 147–168, 1999.
  * [57] R. Karnatak, R. Ramaswamy, and U. Feudel, “Conjugate coupling in ecosystems: Cross-predation stabilizes food webs,” Chaos, Solitons & Fractals, vol. 68, pp. 48–57, 2014.
  * [58] G. C. Sethia and A. Sen, “Chimera states: the existence criteria revisited,” Physical Review Letters, vol. 112, no. 14, p. 144101, 2014.
  * [59] C. Hens, A. Mishra, P. Roy, A. Sen, and S. Dana, “Chimera states in a population of identical oscillators under planar cross-coupling,” Pramana, vol. 84, no. 2, pp. 229–235, 2015.
  * [60] A. Mishra, C. Hens, M. Bose, P. K. Roy, and S. K. Dana, “Chimeralike states in a network of oscillators under attractive and repulsive global coupling,” Physical Review E, vol. 92, no. 6, p. 062920, 2015.
  * [61] R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical Journal, vol. 1, no. 6, p. 445, 1961.
  * [62] M. Bazhenov, R. Huerta, M. Rabinovich, and T. Sejnowski, “Cooperative behavior of a chain of synaptically coupled chaotic neurons,” Physica D: Nonlinear Phenomena, vol. 116, no. 3-4, pp. 392–400, 1998.
  * [63] N. F. Rulkov, M. M. Sushchik, L. S. Tsimring, and H. D. Abarbanel, “Generalized synchronization of chaos in directionally coupled chaotic systems,” Physical Review E, vol. 51, no. 2, p. 980, 1995.
  * [64] S. Boccaletti, D. Valladares, J. Kurths, D. Maza, and H. Mancini, “Synchronization of chaotic structurally nonequivalent systems,” Physical Review E, vol. 61, no. 4, p. 3712, 2000.
  * [65] Y. Kuramoto and D. Battogtokh, “Coexistence of coherence and incoherence in nonlocally coupled phase oscillators,” Nonlinear Phenomena in Complex Systems, vol. 5, p. 380, 2002.
  * [66] I. Omelchenko, E. Omel’chenko, P. Hövel, and E. Schöll, “When nonlocal coupling between oscillators becomes stronger: patched synchrony or multichimera states,” Physical Review Letters, vol. 110, no. 22, p. 224101, 2013.
  * [67] H. Sakaguchi, “Instability of synchronized motion in nonlocally coupled neural oscillators,” Physical Review E, vol. 73, no. 3, p. 031907, 2006.
  * [68] S. Saha, N. Bairagi, and S. K. Dana, “Chimera states in ecological network under weighted mean-field dispersal of species,” Frontiers in Applied Mathematics and Statistics, vol. 5, p. 15, 2019.

|  Suman Saha received the B.TECH., and M.TECH. degrees in Electronics and
Communication Engineering both from West Bengal University of Technology,
Kolkata, India, in 2006 and 2010, respectively. He completed his Ph.D. at
Jadavpur University, Kolkata-700032, India, in 2020. His research interests
include nonlinear dynamics, complex network, synchronization in dynamical
systems and collective behaviors- cluster and chimera states, extreme events
in neuronal systems, computational biology, the role of network structure and
initial infectivity pattern in progression of infectious disease, predicting
complex signal using machine learning, and neuroscience: cognition and brain
dynamics.
---|---