which models the instrument output resulting from observing a galaxy spectrum
on a patch of simulated sky, including instrumental and astrophysical noise.
(Euclid Collaboration: Lusso et al. 2023) have performed a similar analysis
assessing NISP performance for active galaxies, using mock AGN spectra created
from a library of empirical templates. We also refer the reader to ongoing
work by Mancini et al. (in prep.), which will further assess the impact of
Euclid’s instrumental effects on emission-line forecasts, with an empirical
approach based on the Millennium Mambo catalogues, complementary to ours.

Our approach further assumes that Euclid can recover separable line fluxes for
the H$\alpha$ and [N ii]$\lambda 6584$ line. Given the resolution of the NISP
spectrometer, the $\rm H\alpha+[{N\textsc{ii}}]\lambda
6584+[{N\textsc{ii}}]\lambda 6548$ complex will be blended in Euclid spectra
(see also Euclid Collaboration: Lusso et al. 2023). However, in the OU-SPE
spectroscopic pipeline (Le Brun et al., in prep) two different approaches are
implemented to measure the line fluxes: a direct integration (DI) and a
Gaussian-fit (GF) method. In the GF method, the $\rm H\alpha+[{N\textsc{ii}}]$
is modeled as three Gaussians, where the free parameters are the amplitude of
H$\alpha$ and of [N ii]$\lambda 6584$, the width of the line, the position of
H$\alpha$, and the value of the continuum. The flux ratio of the [N ii] lines
is set to 1/3 and their positions are fixed to that of the H$\alpha$ line. As
discussed in Le Brun et al. (in prep), the GF method provides an estimate of
the deconvolved H$\alpha$ and [N ii] fluxes.

##  9 Conclusion

In this work, we presented optical emission-line predictions at intermediate
redshifts for the upcoming Euclid Wide and Deep Surveys, addressing the lack
of comprehensive theoretical guidance in this regime. We followed a
methodology adapted from Hirschmann et al. (2017, 2019, 2023a, 2023b) to
construct emission-line catalogues by coupling galaxies from a mock light cone
based on the Gaea semi-analytic model to state-of-the-art photoionisation
models. This enabled us to self-consistently compute the emission lines of
galaxies at different cosmic epochs originating from young star clusters
(Gutkin et al. 2016), AGN narrow-line regions (Feltre et al. 2016), and post-
asymptotic giant branch stellar populations (Hirschmann et al. 2017). As a
last step, we validated the resulting emission-line catalogues by comparing
its predictions to observational data and well-calibrated theoretical
predictions for low redshifts (Sect. 3). This framework represents the first
emission-line catalogue based on a semi-analytic model which contains self-
consistent modelling of the ionising contributions not only from young star
clusters, but also from other ionising sources, such as AGN, and post-AGB
stellar populations.

We then focused on spectroscopic diagnostics in the redshift range 0.4–2.5
based on the seven optical emission lines: H$\alpha$, H$\beta$, [Sii], [Nii],
[Oi], [Oiii], and [Oii]. In order to make targeted predictions for their
observability with Euclid, we further computed observer-like fluxes based on
the location and redshift of each galaxy in the light cone and modelled
attenuation due to interstellar dust using the Calzetti et al. (2000)
relation, with attenuation $A_{V}$ modelled according to the mass-dependent
Garn & Best (2010) relation.

Our main results from the analysis can be summarised as follows:

  1. 1.

We tested how emission lines trace scaling relations in two observing
scenarios: observing H$\alpha$ in the EWS and the BPT lines, namely H$\alpha$,
H$\beta$, [O iii], and [N ii], in the EDS (Fig. 8). In both cases, the
resulting observable populations bias standard scaling relation towards high
stellar and halo masses, high specific SFR, and high metallicities, meaning
Euclid will predominantly observe line-emitting galaxies that are massive
($M_{\star}\gtrsim 10^{9}\,$\mathrm{\text{${M}_{\odot}$}}$$), star-forming
($\mathrm{sSFR}>10^{-10}\,\mathrm{yr^{-1}}$), and metal-rich
($\mathrm{\operatorname{log_{10}}(O/H)+12}>8$). We predict that both survey
configurations will recover galaxies containing AGN with black hole masses
between $10^{6}$–$10^{9.5}$\mathrm{\text{${M}_{\odot}$}}$$ and bolometric
luminosities of $10^{39}$–$10^{46}\rm\,erg\,s^{-1}$ (Fig. 9). We estimate that
the influence of interstellar dust could reduce observable percentages by an
additional 20–30% with respect to the intrinsic population, which may pose
challenges in measuring fainter lines, especially toward higher redshifts
(Fig. 10). If accounting for dust attenuation, we anticipate that at redshift
less than 1, Euclid will successfully capture approximately 30–70% of both SF
and active galaxies emitting H$\alpha$, [N ii], [S ii], and [O iii], given our
Gaea-lc galaxy sample with a mass resolution limit of
$10^{9}\,$\mathrm{\text{${M}_{\odot}$}}$$ and $H$-band magnitude cut of 25. At
higher redshifts, these percentages decline to below 10%. H$\beta$, [O ii],
and [O i] exhibit particularly faint signatures and, consequently, we expect
observable percentages to be limited to below 5% for both SF and active
galaxies at low redshifts, which are reduced to below 1% with increasing
redshift.

  2. 2.

For EDS-observable galaxies, we expect [O iii]/H$\beta$ versus [N
ii]/H$\alpha$ BPT diagrams to continue to distinguish between SF-dominated,
composite, and AGN-dominated galaxies up to at least redshift 1.8 (Fig. 11).
This can be attributed to the bias toward metal-rich systems, introduced by
requiring the observability of all four emission lines. After including the
impact of dust attenuation, we expect this to be the case for 1.4% of Gaea-lc
galaxies, with an upper limit of 11.8% in the no dust scenario.

  3. 3.

We find that relationships between H$\alpha$ and $\rm[{O\textsc{iii}}]+H\beta$
luminosities and the SFR show only a negligible, if any, evolution with
increasing redshift, when compared to the local calibrations from Kennicutt &
Evans (2012) and Osterbrock & Ferland (2006, Fig. 12). This indicates, that
they could be applied in the analysis of future Euclid data. The
$L_{\mathrm{AGN}}$-$L_{[{O\textsc{iii}}]}$ relationship shows a redshift
evolution of up to +0.5 dex with respect to the local Lamastra et al. (2009)
calibration, but generally retains the same slope. As a result, we find that
Lamastra et al. (2009) still provides a reasonable estimate for the AGN
luminosity, especially at the lower end of Euclid’s target redshift range. We
further indicate up to which luminosity thresholds these tracers appear
largely unbiased for EDS-observable and dust-attenuated galaxy populations at
different redshifts. It is important to note that these calibrations are only
valid for the dominant ionisation mechanism in the respective galaxy types,
meaning SF-dominated galaxies for the SFR relations and AGN-dominated galaxies
for the $L_{\mathrm{AGN}}$ relation. As a result, a pre-sorting of galaxies
according to their type is necessary.

  4. 4.

We find that [N ii] emission strongly depends on the $\rm BHAR/SFR$ ratio
(Fig. 13). As a result, we explored various [N ii]-based emission-line tracers
for the $\rm BHAR/SFR$ ratio and discover strong positive correlations with
$\rm[{N\textsc{ii}}]/H\alpha$, $\rm[{N\textsc{ii}}]/[{O\textsc{ii}}]$, and
$\rm[{N\textsc{ii}}]/[{S\textsc{ii}}]$ (Fig. 14). These relationships appear
to be largely redshift-invariant in the 0–2.5 range, are valid across all
galaxy types and show only slight biasing when applying Euclid flux cuts. We
propose these as novel tracers and provide linear best fits between
$\mathrm{-4\leq\operatorname{log_{10}}(BHAR/SFR)\leq-1}$, allowing for the
derivation of the relative strength of star-forming and accretion processes
from combinations of only two emission lines.

  5. 5.

We examined the potential evolution between redshift 0 and 2.5 of the strong
emission-line ratios commonly used to estimate oxygen abundance from observed
spectra (Fig. 15): $\rm[{N\textsc{ii}}]/H\alpha$ (N2),
$\rm[{S\textsc{ii}}]/H\alpha$ (S2), $\rm[{N\textsc{ii}}]/[{S\textsc{ii}}]$
(N2S2), $\rm[{N\textsc{ii}}]/[{O\textsc{ii}}]$ (N2O2),
$\rm[{O\textsc{ii}}]/H\alpha$ (R2), $\rm[{O\textsc{iii}}]/H\alpha$ (R3),
$\rm([{O\textsc{iii}}]+[{O\textsc{ii}}])/H\alpha$ (R23),
$\rm([{O\textsc{iii}}]/H\alpha)+([{S\textsc{ii}}]/H\alpha)$ (RS23),
$\rm[{O\textsc{iii}}]/[{O\textsc{ii}}]$ (O32),
$\rm([{O\textsc{iii}}]/H\alpha)/([{N\textsc{ii}}]/H\alpha)$ (O3N2), and
$\rm([{O\textsc{iii}}]/H\alpha)/([{S\textsc{ii}}]/H\alpha)$ (O3S2). We found
that, in general, they evolve away from their locally established
calibrations. This evolution manifests as a shift in normalisation, with
metallicity estimators gradually moving towards either lower or higher values,
while the overall shape remains consistent with local calibrations. This is in
tentative agreement with current, although sparse, JWST data. Notably, the
N2O2 and R23 ratios exhibit the weakest evolutions, suggesting that they are
the most robust to extend to higher redshifts without adjustments. We assessed
the robustness of various metallicity estimators for EDS-like samples and
found that [O ii]- and [O iii]-based estimators are reliable within the
observable $\rm 12+\operatorname{log_{10}}(O/H)$ range of approximately 8–9.
However differential dust extinction should be considered when using widely
spaced wavelength line combinations, like N2O2, O3N2, O3S2, O32, R23, and RS23
estimators.

In summary, the comprehensive predictions presented in this paper offer
valuable insights into emission-line properties of galaxy populations at
intermediate redshifts, and their relationship to the ionising source
properties and local gas conditions. This represents theoretical guidance for
a redshift range that has seen limited spectroscopic coverage thus far, and we
expect it to serve as a reference for interpreting results from the upcoming
Euclid surveys, as well as other spectroscopic surveys with instruments like
DESI and VLT/MOONs.

###### Acknowledgements.

The Euclid Consortium acknowledges the European Space Agency and a number of
agencies and institutes that have supported the development of Euclid, in
particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian
Science Policy, the Canadian Euclid Consortium, the French Centre National
d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish
Space Research Institute, the Fundação para a Ciência e a Tecnologia, the
Ministerio de Ciencia e Innovación, the National Aeronautics and Space
Administration, the National Astronomical Observatory of Japan, the
Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the
Romanian Space Agency, the State Secretariat for Education, Research and
Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom
Space Agency. A complete and detailed list is available on the Euclid web site
(http://www.euclid-ec.org). L.S., M.H., and A.P. acknowledge funding from the
Swiss National Science Foundation (SNF) via a PRIMA Grant PR00P2 193577 “From
cosmic dawn to high noon: the role of black holes for young galaxies”. We also
acknowledge the computing centre of INAF-OATs, under the coordination of the
CHIPP project (Taffoni et al. 2020), for the availability of computing
resources and support. C.S. acknowledges support by NASA ROSES grant
12-EUCLID12-0004. A.F. acknowledges the support from project ”VLT-MOONS” CRAM
1.05.03.07, INAF Large Grant 2022 ”The metal circle: a new sharp view of the
baryon cycle up to Cosmic Dawn with the latest generation IFU facilities” and
INAF Large Grant 2022 ”Dual and binary SMBH in the multi-messenger era”. V.A.
acknowledges support from INAF-PRIN 1.05.01.85.08

## Data availability

An introduction to Gaea, a list of our recent work, as well as data file
containing published model predictions, can be found at
https://sites.google.com/inaf.it/gaea/home.

## References

  * Alarie & Morisset (2019) Alarie, A. & Morisset, C. 2019, Rev. Mex. Astron. Astrofis., 55, 377
  * Bagley et al. (2020) Bagley, M. B., Scarlata, C., Mehta, V., et al. 2020, ApJ, 897, 98
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Baugh et al. (2022) Baugh, C. M., Lacey, C. G., Gonzalez-Perez, V., & Manzoni, G. 2022, MNRAS, 510, 1880
  * Belfiore et al. (2016) Belfiore, F., Maiolino, R., Maraston, C., et al. 2016, MNRAS, 461, 3111
  * Bian et al. (2010) Bian, F., Fan, X., Bechtold, J., et al. 2010, ApJ, 725, 1877
  * Bian et al. (2017) Bian, F., Fan, X., Mcgreer, I., Cai, Z., & Jiang, L. 2017, ApJL, 837, L12
  * Bian et al. (2018) Bian, F., Kewley, L. J., & Dopita, M. A. 2018, ApJ, 859, 175
  * Blaizot et al. (2005) Blaizot, J., Wadadekar, Y., Guiderdoni, B., et al. 2005, MNRAS, 360, 159
  * Bresolin (2006) Bresolin, F. 2006, ApJ, 656, 186
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Cameron et al. (2023) Cameron, A. J., Katz, H., Rey, M. P., et al. 2023, MNRAS, 522, L89
  * Carton et al. (2017) Carton, D., Carton, D., Brinchmann, J., et al. 2017, MNRAS, 468, 2140
  * Carvalho et al. (2020) Carvalho, S. P., Dors, O. L., Cardaci, M. V., et al. 2020, MNRAS, 492, 5675
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Charlot & Longhetti (2001) Charlot, S. & Longhetti, M. 2001, MNRAS, 323, 887
  * Chevallard & Charlot (2016) Chevallard, J. & Charlot, S. 2016, MNRAS, 462, 1415
  * Chomiuk & Povich (2011) Chomiuk, L. & Povich, M. S. 2011, AJ, 142
  * Coil et al. (2014) Coil, A. L., Aird, J., Reddy, N., et al. 2014, ApJ, 801, 35
  * Colbert et al. (2013) Colbert, J. W., Teplitz, H., Atek, H., et al. 2013, ApJ, 779, 34
  * Cullen et al. (2017) Cullen, F., Mclure, R. J., Khochfar, S., Dunlop, J. S., & Dalla Vecchia, C. 2017, MNRAS, 470, 3006
  * Curti et al. (2022) Curti, M., D’Eugenio, F., Carniani, S., et al. 2022, MNRAS, 518, 425
  * Curti et al. (2019) Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2019, MNRAS, 491, 944
  * De Barros et al. (2019) De Barros, S., Oesch, P. A., Labbé, I., et al. 2019, MNRAS, 489, 2355
  * De Lucia & Blaizot (2007) De Lucia, G. & Blaizot, J. 2007, MNRAS, 375, 2
  * De Lucia et al. (2018) De Lucia, G., Hirschmann, M., & Fontanot, F. 2018, MNRAS, 482, 5041
  * De Lucia et al. (2014) De Lucia, G., Tornatore, L., Frenk, C. S., et al. 2014, MNRAS, 445, 970
  * Domínguez et al. (2012) Domínguez, A., Siana, B., Henry, A. L., et al. 2012, ApJ, 763, 145
  * Dors et al. (2017) Dors, O. L., Arellano-Córdova, K. Z., Cardaci, M. V., & Hägele, G. F. 2017, MNRAS, 468, L113
  * Dors et al. (2021) Dors, O. L., Contini, M., Riffel, R. A., et al. 2021, MNRAS, 501, 1370
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
  * Euclid Collaboration: Gabarra et al. (2023) Euclid Collaboration: Gabarra, L., Mancini, C., Rodriguez Munoz, L., et al. 2023, A&A, 676, A34
  * Euclid Collaboration: Lusso et al. (2023) Euclid Collaboration: Lusso, E., Fotopoulou, S., Selwood, M., et al. 2023, arXiv:2311.12096
  * Euclid Collaboration: Scaramella et al. (2022) Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112
  * Euclid Collaboration: Schirmer et al. (2022) Euclid Collaboration: Schirmer, M., Jahnke, K., Seidel, G., et al. 2022, 662, A92
  * Feltre et al. (2016) Feltre, A., Charlot, S., & Gutkin, J. 2016, MNRAS, 456, 3354
  * Feltre et al. (2023) Feltre, A., Gruppioni, C., Marchetti, L., et al. 2023, A&A, 675, A47
  * Ferland & Netzer (1983) Ferland, G. J. & Netzer, H. 1983, ApJ, 264, 105
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., Van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Fontanot et al. (2020) Fontanot, F., de Lucia, G., Hirschmann, M., et al. 2020, MNRAS, 496, 3943
  * Fontanot et al. (2017) Fontanot, F., Hirschmann, M., & De Lucia, G. 2017, ApJL, 842, L14
  * Fujita et al. (2003) Fujita, S. S., Ajiki, M., Shioya, Y., et al. 2003, ApJ, 586, L115
  * Gallego et al. (1995) Gallego, J., Zamorano, J., Aragón-Salamanca, A., & Rego, M. 1995, ApJ, 455, L1
  * Garg et al. (2022) Garg, P., Narayanan, D., Byler, N., et al. 2022, ApJ, 926, 80
  * Garn & Best (2010) Garn, T. & Best, P. 2010, MNRAS, 409, 421
  * Geach et al. (2008) Geach, J. E., Smail, I., Best, P. N., et al. 2008, MNRAS, 388, 1473
  * Genzel et al. (2011) Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101
  * Georgantopoulos & Akylas (2010) Georgantopoulos, I. & Akylas, A. 2010, A&A, 509, A38
  * Groves et al. (2006) Groves, B., Heckman, T., & Kauffmann, G. 2006, MNRAS, 371, 1559
  * Gutkin et al. (2016) Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757
  * Hainline et al. (2009) Hainline, K. N., Shapley, A. E., Kornei, K. A., et al. 2009, ApJ, 701, 52
  * Harikane et al. (2023) Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 18
  * Hasinger (2008) Hasinger, G. 2008, A&A, 490, 905
  * Hayes et al. (2009) Hayes, M., Schaerer, D., & Göran Östlin. 2009, A&A
  * Hayes et al. (2010) Hayes, M., Schaerer, D., & Östlin, G. 2010, A&A, 509, L5
  * Heckman et al. (1980) Heckman, T. M., Heckman, & M., T. 1980, A&A, 87, 152
  * Hirschmann et al. (2023a) Hirschmann, M., Charlot, S., Feltre, A., et al. 2023a, MNRAS, 526, 3619
  * Hirschmann et al. (2017) Hirschmann, M., Charlot, S., Feltre, A., et al. 2017, MNRAS, 472, 2468
  * Hirschmann et al. (2019) Hirschmann, M., Charlot, S., Feltre, A., et al. 2019, MNRAS, 487, 333
  * Hirschmann et al. (2023b) Hirschmann, M., Charlot, S., & Somerville, R. S. 2023b, MNRAS, 526, 3504
  * Hirschmann et al. (2016) Hirschmann, M., De Lucia, G., & Fontanot, F. 2016, MNRAS, 461, 1760
  * Hopkins et al. (2001) Hopkins, A. M., Connolly, A. J., Haarsma, D. B., & Cram, L. E. 2001, AJ, 122, 288
  * Hopkins et al. (2000) Hopkins, A. M., Connolly, A. J., & Szalay, A. S. 2000, AJ, 120, 2843
  * Hopkins et al. (2003) Hopkins, A. M., Miller, C. J., Nichol, R. C., et al. 2003, ApJ, 599, 971
  * Jaskot & Oey (2013) Jaskot, A. E. & Oey, M. S. 2013, ApJ, 766, 23
  * Jensen et al. (1976) Jensen, E. B., Strom, K. M., & Strom, S. E. 1976, ApJ, 209, 748
  * Jones et al. (2015) Jones, T., Martin, C., & Cooper, M. C. 2015, ApJ, 813, 126
  * Juneau et al. (2014) Juneau, S., Bournaud, F., Charlot, S., et al. 2014, ApJ, 788, 88
  * Kashino et al. (2014) Kashino, D., Silverman, J. D., Rodighiero, G., et al. 2014, ApJL, 777, L8
  * Kashino et al. (2017) Kashino, D., Silverman, J. D., Sanders, D., et al. 2017, ApJ, 835, 88
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kennicutt (1992) Kennicutt, R. C. 1992, ApJ, 388, 310
  * Kennicutt (1998) Kennicutt, R. C. 1998, ARA&A, 36, 189
  * Kennicutt & Evans (2012) Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531
  * Kewley & Dopita (2002) Kewley, L. J. & Dopita, M. A. 2002, ApJS, 142, 35
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961
  * Kewley et al. (2001) Kewley, L. J., Heisler, C. A., Dopita, M. A., et al. 2001, ApJS, 132, 37
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511
  * Khostovan et al. (2015) Khostovan, A. A., Sobral, D., Mobasher, B., et al. 2015, MNRAS, 452, 3948
  * Kocevski et al. (2022) Kocevski, D. D., Barro, G., McGrath, E. J., et al. 2022, ApJL, 946, L14
  * Kocevski et al. (2023) Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, ApJL, 954, L4
  * Kornei et al. (2012) Kornei, K. A., Shapley, A. E., Martin, C. L., et al. 2012, ApJ, 758, 135
  * Kroupa et al. (2003) Kroupa, P., Weidner, C., Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076
  * Lamastra et al. (2009) Lamastra, A., Bianchi, S., Matt, G., et al. 2009, A&A, 504, 73
  * Laseter et al. (2023) Laseter, I. H., Maseda, M. V., Curti, M., et al. 2023, arXiv:2306.03120
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193
  * Lehnert et al. (2009) Lehnert, M. D., Nesvadba, N. P. H., Le Tiran, L., et al. 2009, ApJ, 699, 1660
  * Liu et al. (2008) Liu, X., Shapley, A., Coil, A., Brinchmann, J., & Ma, C. 2008, ApJ, 678, 758
  * Ly et al. (2007) Ly, C., Malkan, M. A., Kashikawa, N., et al. 2007, ApJ, 657, 738
  * Maciaszek et al. (2022) Maciaszek, T., Ealet, A., Gillard, W., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12180, Space Telescopes and Instrumentation 2022: Optical, Infrared, and Millimeter Wave, ed. L. E. Coyle, S. Matsuura, & M. D. Perrin, arXiv:2210.10112
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, Annu. Rev. Astron. Astrophys., 52, 415
  * Maheson et al. (2024) Maheson, G., Maiolino, R., Curti, M., et al. 2024, MNRAS, 527, 8213
  * Maiolino & Mannucci (2019) Maiolino, R. & Mannucci, F. 2019, ARA&A, 27, 3
  * Maiolino et al. (2008) Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463
  * Marino et al. (2013) Marino, R. A., Rosales-Ortega, F. F., Sánchez, S. F., et al. 2013, A&A, 559, A114
  * Masters et al. (2014) Masters, D., McCarthy, P., Siana, B., et al. 2014, ApJ, 785, 153
  * McDonald et al. (2021) McDonald, M., McNamara, B. R., Calzadilla, M. S., et al. 2021, ApJ, 908, 85
  * Mclure et al. (2017) Mclure, R. J., Dunlop, J. S., Cullen, F., et al. 2017, MNRAS, 476, 3991
  * Mehta et al. (2015) Mehta, V., Scarlata, C., Colbert, J. W., et al. 2015, ApJ, 811, 141
  * Moustakas et al. (2006) Moustakas, J., Kennicutt, R. C., & Tremonti, C. A. 2006, ApJ, 642, 775
  * Murphy et al. (2011) Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, AJ, 737, 67
  * Nagao et al. (2006) Nagao, T., Maiolino, R., & Marconi, A. 2006, A&A, 459, 85
  * Netzer et al. (2006) Netzer, H., Mainieri, V., Rosati, P., & Trakhtenbrot, B. 2006, A&A, 453, 525
  * Newman et al. (2012) Newman, S. F., Genzel, R., Förster-Schreiber, N. M., et al. 2012, ApJ, 761, 43
  * Orsi et al. (2014) Orsi, A., Padilla, N., Groves, B., et al. 2014, MNRAS, 443, 799
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei (University Science Books)
  * Pagel et al. (1979) Pagel, B. E., Edmunds, M. G., Blackwell, D. E., et al. 1979, MNRAS, 189, 95
  * Panessa et al. (2006) Panessa, F., Bassani, L., Cappi, M., et al. 2006, A&A, 455, 173
  * Pellegrini et al. (2020) Pellegrini, E. W., Rahner, D., Reissl, S., et al. 2020, MNRAS, 496, 339
  * Perez-Gonzalez et al. (2003) Perez-Gonzalez, P. G., Zamorano, J., Gallego, J., Aragon-Salamanca, A., & de Paz, A. G. 2003, ApJ, 591, 827
  * Pérez-Montero & Contini (2009) Pérez-Montero, E. & Contini, T. 2009, MNRAS, 398, 949
  * Peterson (2006) Peterson, B. M. 2006, LNP, 693, 77
  * Pettini & Pagel (2004) Pettini, M. & Pagel, B. E. 2004, MNRAS, 348, L59
  * Plat et al. (2019) Plat, A., Charlot, S., Bruzual, G., et al. 2019, MNRAS, 490, 978
  * Popping et al. (2017) Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152
  * Pozzetti et al. (2016) Pozzetti, L., Hirata, C. M., Geach, J. E., et al. 2016, A&A, 590, A3
  * Price et al. (2014) Price, S. H., Kriek, M., Brammer, G. B., et al. 2014, ApJ, 788, 13
  * Racca et al. (2016) Racca, G. D., Laureijs, R., Stagnaro, L., et al. 2016, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, 9904, 99040O
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, 22
  * Rich et al. (2010) Rich, J. A., Dopita, M. A., Kewley, L. J., & Rupke, D. S. N. 2010, ApJ, 721, 505
  * Rich et al. (2011) Rich, J. A., Kewley, L. J., & Dopita, M. A. 2011, ApJ, 734, 87
  * Richards et al. (2006) Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, AJ, 131, 2766
  * Sanders et al. (2023a) Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2023a, ApJ, 955, 54
  * Sanders et al. (2023b) Sanders, R. L., Shapley, A. E., Topping, M. W., et al. 2023b, arXiv:2303.08149
  * Shapley et al. (2004) Shapley, A. E., Erb, D. K., Pettini, M., Steidel, C. C., & Adelberger, K. L. 2004, ApJ, 612, 108
  * Shapley et al. (2015) Shapley, A. E., Reddy, N. A., Kriek, M., et al. 2015, ApJ, 801, 88
  * Shapley et al. (2022) Shapley, A. E., Sanders, R. L., Salim, S., et al. 2022, ApJ, 926, 8
  * Shapley et al. (2016) Shapley, A. E., Steidel, C. C., Strom, A. L., et al. 2016, ApJL, 826, L24
  * Sharp & Bland-Hawthorn (2010) Sharp, R. G. & Bland-Hawthorn, J. 2010, ApJ, 711, 818
  * Shim et al. (2009) Shim, H., Colbert, J., Teplitz, H., et al. 2009, ApJ, 696, 785
  * Shimizu et al. (2016) Shimizu, I., Inoue, A. K., Okamoto, T., & Yoshida, N. 2016, MNRAS, 461, 3563
  * Shioya et al. (2008) Shioya, Y., Taniguchi, Y., Sasaki, S. S., et al. 2008, APJS, 175, 128
  * Singh et al. (2013) Singh, R., Van De Ven, G., Jahnke, K., et al. 2013, A&A, 558, A43
  * Sobral et al. (2009) Sobral, D., Best, P. N., Geach, J. E., et al. 2009, MNRAS, 398, 75
  * Sobral et al. (2015) Sobral, D., Matthee, J., Best, P. N., et al. 2015, MNRAS, 451, 2303
  * Sobral et al. (2013) Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128
  * Soto et al. (2012) Soto, K. T., Martin, C. L., Prescott, M. K., & Armus, L. 2012, ApJ, 757, 86
  * Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629
  * Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289
  * Steidel et al. (2014) Steidel, C. C., Rudie, G. C., Strom, A. L., et al. 2014, ApJ, 795, 165
  * Strom et al. (2017) Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164
  * Suzuki et al. (2016) Suzuki, T. L., Kodama, T., Sobral, D., et al. 2016, MNRAS, 462, 181
  * Taffoni et al. (2020) Taffoni, G., Becciani, U., Garilli, B., et al. 2020, ASPC, 527, 307
  * Teplitz et al. (2000) Teplitz, H. I., Malkan, M. A., Steidel, C. C., et al. 2000, ApJ, 542, 18
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898
  * Tresse & Maddox (1998) Tresse, L. & Maddox, S. 1998, ApJ, 495, 691
  * Tresse et al. (2002) Tresse, L., Maddox, S. J., Fevre, O. L., & Cuby, J. G. 2002, MNRAS, 337, 369
  * Urry & Padovani (1995) Urry, C. M. & Padovani, P. 1995, PASP, 107, 803
  * Veilleux & Osterbrock (1987) Veilleux, S. & Osterbrock, D. E. 1987, ApJS, 63, 295
  * Vidal-García et al. (2024) Vidal-García, A., Plat, A., Curtis-Lake, E., et al. 2024, MNRAS, 527, 7217
  * Villar et al. (2007) Villar, V., Gallego, J., Pérez-González, P. G., et al. 2007, ApJ, 677, 169
  * Weistrop et al. (2012) Weistrop, D., Nelson, C. H., Angione, R., et al. 2012, ApJ, 143, 11
  * Wilkins et al. (2020) Wilkins, S. M., Lovell, C. C., Fairhurst, C., et al. 2020, MNRAS, 493, 6079
  * Yabe et al. (2012) Yabe, K., Ohta, K., Iwamuro, F., et al. 2012, PASJ, 64, 601
  * Yan et al. (1999) Yan, L., McCarthy, P. J., Freudling, W., et al. 1999, ApJL, 519, L47
  * Zahid et al. (2013) Zahid, H. J., Yates, R. M., Kewley, L. J., & Kudritzki, R. P. 2013, ApJ, 763, 10
  * Zoldan et al. (2016) Zoldan, A., De Lucia, G., Xie, L., Fontanot, F., & Hirschmann, M. 2016, MNRAS, 465, 2236

##  Appendix A Redshift evolution of the ionisation parameter

Throughout this study, we have explained the evolution of emission-line
properties and their relationships to galaxy properties, like the gas-phase
metallicity and the location of star-forming galaxies in the BPT diagrams,
using the redshift evolution of the ionisation parameter. We demonstrate this
evolution at fixed metallicity and fixed stellar mass.

In Fig. 16, we show in different redshift bins (using the colour coding as in
Fig. 12) the average ionisation parameter $U_{\star}$ for SF photoionisation
models against global O/H abundance. We compare this against the empirical
relation from Carton et al. (2017, dashed line), which has been derived from
local star-forming SDSS galaxies. Across all redshift bins, our prediction
from Gaea presents an anti-correlation between the ionisation parameter and
the metallicity, in agreement with Carton et al. (2017). At fixed redshift,
the ionisation parameter decreases by around 1.5 dex between
$\mathrm{12+\operatorname{log_{10}}(O/H)=7}$ and 9.5. We further note that the
agreement is best in the lower redshift bins, while higher redshifts exhibit a
systematically increased average ionisation parameter at fixed metallicity.
This is caused by increased star-formation rates and gas densities at
increased redshift.

The increase of the ionisation parameter with redshift persists when keeping
the stellar mass fixed. In Fig. 17, we present the redshift evolution of
average SF ionisation parameter in different stellar mass bins (see colour
code in legend). Due to the anti-correlation of metallicity and ionisation
parameter, and the correlation of mass and metallicity (see Fig. 8 and
Maiolino et al. 2008), the ionisation parameter is most elevated at the lowest
stellar masses.

Figure 16: Average ionisation parameter $U_{\star}$ for SF models against O/H
abundance in different redshift bins (same colour coding as in Fig. 12).
Alongside, we show the relation from Carton et al. (2017). Figure 17: Redshift
evolution of average ionisation parameter $U_{\star}$ for SF models in
different stellar mass bins (colour coding as indicated in legend).


# Diffusion and pattern formation in spatial games

Alexandre Champagne-Ruel  Département de physique, Université de Montréal,
Canada<EMAIL_ADDRESS>Sascha Zakaib-Bernier
Département de physique, Université de Montréal, Canada  Paul Charbonneau
Département de physique, Université de Montréal, Canada

(July 1, 2024)

###### Abstract

Diffusion plays an important role in a wide variety of phenomena, from
bacterial quorum sensing to the dynamics of traffic flow. While it generally
tends to level out gradients and inhomogeneities, diffusion has nonetheless
been shown to promote pattern formation in certain classes of systems.
Formation of stable structures often serves as a key factor in promoting the
emergence and persistence of cooperative behavior in otherwise competitive
environments, however an in-depth analysis on the impact of diffusion on such
systems is lacking. We therefore investigate the effects of diffusion on
cooperative behavior using a cellular automaton (CA) model of the noisy
spatial iterated prisoner’s dilemma (IPD), physical extension and
stochasticity being unavoidable characteristics of several natural phenomena.
We further derive a mean-field (MF) model that captures the 3-species
predation dynamics from the CA model and highlight how pattern formation
arises in this new model, then characterize how including diffusion by
interchange similarly enables the emergence of large scale structures in the
CA model as well. We investigate how these emerging patterns favors
cooperative behavior for parameter space regions where IPD error rates
classically forbid such dynamics. We thus demonstrate how the coupling of
diffusion with non-linear dynamics can, counter-intuitively, promote large
scale structure formation and in return establish new grounds for cooperation
to take hold in stochastic spatial systems.

##  I Introduction

Cooperative phenomena are ubiquitous in nature. While biological evolution
hinges on the principle of selection, which also implies competition,
instances of cooperation are widespread both at small Turner and Chao (1999);
Bohl  _et al._ (2014); Hong  _et al._ (1992); Sievers and von Kiedrowski
(1994); Lee  _et al._ (1997) and larger Dugatkin (1997, 1988); Heinsohn and
Packer (1995); Kay  _et al._ (2020); Lombardo (1985); Milinski (1987) scales.
The evolution of mitochondrial cells is often cited as one prime example of
mutually beneficial behavior Sagan (1967); Margulis (1981), and cooperation is
believed to have been instrumental both in the emergence of life on Earth
Kauffman (1986); Hordijk and Steel (2004); Hordijk (2019) and, more generally,
in major evolutionary transitions throughout history that would have
necessitated the establishment of new forms of cooperation Queller (1997);
Szathmáry (2015); Massey and Mishra (2018).

Game theory is the framework of choice for examining the emergence of
cooperation in nature, a problem that has motivated extensive analysis across
various fields Colman and Colman (2013); Myerson (2013); Watson (2016), from
ecology Sigmund (1995) to computational neuroscience Keinan  _et al._
(2006)—and more recently in the fields of physics and biochemistry Schuster
_et al._ (2008). One way to study cooperative phenomena is through the
formalism of the Prisoner’s Dilemma (PD) Fudenberg and Tirole (1991): in this
game, two players must decide whether to cooperate or defect, without prior
knowledge of their opponent’s choice. The PD’s formalism can be defined as a
score matrix that highlights the tension between rational, selfish behaviors
and mutual cooperation that benefits the group as a whole (Table 1). A game is
classified as PD when the Reward $R$ (both players cooperate), the Punishment
$P$ (they both defect), the Temptation $T$ (one exploits its opponent) and the
Sucker’s payoff $S$ (one is being exploited) satisfy the constraint $T>R>P>S$
Nowak (2006).

Numerous competitive scenarios align with the structure of the PD. This
includes interactions among viruses Turner and Chao (1999) or fish Dugatkin
(1997) to the dynamics of WWI trench warfare Trivers (1985). As such, the PD
stands as the model of choice in studying what factors drive, influence, and
favor the emergence of cooperative behavior.

| Cooperation | Defection
---|---|---
Cooperation | 3 ($R$) | 0 ($S$)
Defection | 5 ($T$) | 1 ($P$)
Table 1: Canonical score matrix defining the Prisoner’s Dilemma (PD). Two
players, without knowing their opponent’s choice in advance, can either choose
to cooperate or defect. While exploitation (i.e., defecting while the opponent
cooperates) leads to the highest score for one player, the average score
between both players is higher when they simultaneously cooperate. The PD thus
highlights the tension between rational, selfish behavior, and an alternative
that benefits the group as a whole.

In real world scenarios, players rarely encounter one another only once.
Accordingly, the PD can be extended to account for repeated interaction
between players. The Iterated Prisoner’s Dilemma (IPD) defines a game of
length $M$ wherein two players will repeatedly choose between cooperation and
defection, guided by a specific strategy. Strategies range from
straightforward ones such as ”Always Cooperate” (ALLC) and ”Always Defect”
(ALLD) to reactive ones such as ”Tit-for-Tat” (TFT), which consistently
mirrors the opponent’s last move. Other well-known strategies include ”Playing
at Random”, ”Win-Stay-Lose-Shift”, ”Generous Tit-for-Tat”, among others
Delahaye and Mathieu (1995); Nowak (2006); Jurišić  _et al._ (2012).

Despite the abundant literature aiming to identify the optimal strategy for
the IPD Haider (2005); Sandholm and Crites (1996), TFT consistently emerges as
a successful contender against a wide array of opponents. The simple act of
mirroring the opponent’s last move remains unparalleled in versatility, even
when strategy makers know in advance that they would be playing against it
Axelrod and Hamilton (1981); Axelrod (1984). Just as the PD is considered the
cornerstone for studying cooperative phenomena, TFT is regarded as its pivot
Nowak and Sigmund (1992). Therefore, we explore in what follows the dynamics
between the most basic IPD strategies—ALLC, ALLD and TFT—as they are more
likely to emerge spontaneously in such a competitive setting.

In the context of the IPD, one important limitation of TFT is however its
susceptibility to the presence of errors. Given that statistical physics
asserts all physical systems experience fluctuations at non-zero temperatures
Reif (2009), this parameter can become crucial across a wide variety of
systems. If one of two cooperating TFT players mistakenly defects, they will
both lock onto a sequence of mutual defection. As illustrated on Fig. 1, the
average score of TFT playing the IPD against itself will quickly plummet with
an increase in the error rate Nowak (2006). Additionally, the absence of an
error-correction mechanism in TFT makes it susceptible to being superseded by
unconditional cooperators such as ALLC. Unlike TFT, ALLC doesn’t engage in
similar mutual defection cycles at non-zero error rates and can dominate once
their predators are removed.

Figure 1: Mean score for two TFT players competing against each other for
$M=1000$ moves, averaged over a statistical sample of $1000$ IPD simulations.
(a) As the error rate increases from a small value ($\log_{10}p=-5$, purple
bars) towards its maximum ($\log_{10}p=-1$, red bars), the score distribution
shifts from a narrow distribution with mean $\mu\sim 3.0$ to a broader
distribution where $\mu\sim 2.25$, reflecting the average score across all
four move combinations. (b) Two-dimensional histogram showing the distribution
of error rates for the same TFT players, with the vertical axis representing
variations in the error rate. The mean score peaks when $p$ is minimal (lower
region) and then shifts towards 2.25 (upper region) as the error rate
approaches the critical value $p_{\rm c}\sim M^{-1}=10^{-3}$, at which point
players commit at least one error per game.

Extensive research has revealed how several mechanisms can influence the
emergence of cooperation in challenging environments, such as payoffs Szabó
and Tőke (1998); Lim  _et al._ (2002); Szabó  _et al._ (2005), noise Chen  _et
al._ (2008); Szolnoki  _et al._ (2009); Fant  _et al._ (2023), or topology
Szabó  _et al._ (2000); Vukov  _et al._ (2006); Amaral  _et al._ (2016).
Likewise, the role of mobility as a phenotypical characteristic (i.e., under
direct control of the agents or their strategy) has been subjected to
comprehensive analyses Dugatkin (1997); Enquist and Leimar (1993); Ferrière
and Michod (1995); Aktipis (2004); Hamilton and Taborsky (2005); Le Galliard
_et al._ (2005); Helbing and Yu (2009); Meloni  _et al._ (2009); Smaldino and
Schank (2012); Yang and Zhang (2023). However, one related but prevalent
feature spanning a wide variety of physical systems, molecular diffusion—i.e.,
the Brownian-like (random) movement that characterizes the particles of any
physical system at non-zero temperature—, has been comparatively less explored
in this context. Systems that both rely on cooperative behaviors during their
evolution and simultaneously experience the physical effects of diffusion are
common in nature and across various scientific domains. Significant examples
include many biological systems, such as that of bacterial colonies engaging
in quorum sensing Ben-Jacob  _et al._ (2000) and where diffusion influences
the dynamics of signaling Bryant and Machta (2023), of ant populations that
act in coordinated actions and also rely on the diffusion of pheromones in the
environment to do so Hölldobler and Wilson (1990), or of bird flocking during
the formation of large patterns Cattivelli and Sayed (2011). Many other
instances highlight the pivotal influence of diffusion, from cooperative
chemical species in catalytic reaction systems Turing (1952) to the
coordinated action of vehicles Nagel and Paczuski (1995) or pedestrians von
Schantz and Ehtamo (2015); Yanagisawa (2016) in traffic flow that manifest
diffusion-like dynamics. Previous analyses were performed that included
molecular diffusion in the context of simulations of the IPD Ferriere and
Michod (1996); Vainstein  _et al._ (2007); Sicardi  _et al._ (2009); Droz  _et
al._ (2009); Cheng  _et al._ (2010); Suzuki and Kimura (2011); Yang and Wang
(2011); Gelimson  _et al._ (2013); Vainstein  _et al._ (2014), however an
extensive investigation of how pattern formation characterizing diffusion-like
processes affect the dynamics of the spatial IPD is still lacking.

In what follows, we aim to address this question about the influence of
diffusion on cooperation by characterizing its effect on the interplay between
cooperators and defectors. Using both a deterministic cellular automaton (CA)
model (Section II) and a mean-field (MF) model (Section III) of the spatial
iterated prisoner’s dilemma, we investigate how pattern formation affects the
dominance of strategies within the IPD. While diffusion usually tends to level
out gradients and inhomogeneities, it has long been known that coupling it to
nonlinear dynamics can give rise to pattern formation Turing (1952): our
results indicate that this phenomena further promotes cooperation (Section
IV), which could thus explain how cooperative behavior arises in otherwise
adverse environments where diffusion takes place.

##  II Cellular automaton model for a 3-strategy spatial game

Spatial constraints like topology or dimensionality exert significant
influence on the dynamics of evolutionary games Fu  _et al._ (2010); Lindgren
and Nordahl (1994); Nowak and Sigmund (1999). Given that these are inherent in
various biological, ecological, and social phenomena, we investigate the
iterated prisoner’s dilemma using numerical simulations on a 2D lattice of
sufficient size ($128\times 128$) to reduce finite-size effects such as
stochastic extinction of strategies.

Figure 2: Final lattice state for cellular automaton simulations with varying
error rates $p$. See also Appendix A for animated versions of Panels (a) and
(c). Simulations were carried out on a Cartesian lattice of size $128\times
128$, with periodic boundaries and random initial conditions, over $T=500$
model iterations. Games of the iterated prisoner’s dilemma of $M=2000$ moves
were played. (a) Fractal patterns are formed, across which strategies
propagate through Fisher waves (Fisher, 1937). (b) Filament-like structures of
defectors (ALLD, red) surviving among a TFT-filled lattice (blue). (c)
Clusters of cooperative strategies (ALLC, green) surrounded by parasites
(ALLD, red), in a TFT-filled lattice (blue). (d) Mottled patterns interweaving
cooperators (TFT, blue) and defectors (ALLD, red). (e-f) Close-up view of
stable ALLC clusters surrounded by ALLD strategies. The central strategies
(ALLC) survive through mutual cooperation, yielding the highest scores in
their Moore neighborhood, while the parasites can persist by defecting against
ALLC. The scores shown on panel (f) are averages taken over a statistical
ensemble of $10^{4}$ games.

In the deterministic cellular automaton (CA) model, simulations proceed by
distributing equal populations of each strategy randomly across the lattice.
At each iteration $t_{i}$ in the model, every player on the grid engages in
$M=2000$ rounds of the prisoner’s dilemma with the eight players in its Moore
neighborhood. After these rounds, the players will adopt the highest-scoring
player’s strategy in this neighborhood if it outperforms theirs. These
replacements are deterministic and sequential, but they proceed in a random
order to eliminate any spatial bias. The simulation runs for $T=500$
iterations, ensuring a relaxation of the system’s dependence on the initial
conditions for every realization of the model.

The repeated interaction between players, whereby they engage $M$ times in the
PD, is a key property of the model. These repeated interactions are a
reflection of the fact that systems relying both on cooperative behaviour and
experiencing diffusion are not well-mixed. Hence, in these systems the
characteristic diffusion time-scale is typically longer than the time of
interaction between the individuals, which increases the likelihood of
repeated encounters. While parasitic strategies like ALLD will inevitably
dominate cooperators like ALLC in single interactions, repeated encounters
enable reactive strategies such as TFT to reciprocate any potential
cooperative behaviour. Furthermore, selecting a larger value for $M$ allows
for the exploration of a broader range of error rates while maintaining the
statistical significance of the players’ mistakes over the total number of
games.

Integrating game playing errors in the CA model, which reflects the inevitable
stochasticity of natural processes, further expands the variety of behaviors
the simulation displays. Specifically, in each of the $M$ rounds of the IPD,
players have a probability $p$, drawn from a uniform distribution, of taking
the opposite action (i.e., defecting when supposed to cooperate, or vice-
versa). The spatial IPD has indeed been demonstrated to produce diverse
spatial patterns the error rate is varied: even for low values of $p$
($-5\lesssim\log_{10}p\lesssim-3$), novel behavior emerges, as depicted on
Fig. 2a: ALLC gains the ability to evade ALLD by propagating into areas
dominated by TFT. This sets off a series of wave-fronts where ALLC further
replaces TFT, with a second trailing front across which ALLD replaces ALLC.
Regions dominated by TFT finally expand into ALLD from behind (see Fig. 2a and
animations in Appendix A). Subsequent panels on Fig. 2 display the resulting
lattices from simulations conducted with varying error rates. Formed
structures thereafter include filaments of ALLD strategies interspersed within
a TFT-dominated lattice, clusters of ALLC strategies encircled by ALLD or
mottled patterns of defectors. These spatial constraints can foster
cooperation in conditions where it would otherwise be absent Flores  _et al._
(2022); Nowak  _et al._ (1994); Wang  _et al._ (2012). For example,
unconditional cooperators (ALLC) tend to be quickly eliminated by parasites
(ALLD) unless they can form stable clusters ensuring their survival (Fig. 2c)
Nowak (2006).

Examination of many such simulations indeed reveals that, barring very low or
very high error rates, the spatiotemporal evolution of strategies can—to good
approximation—be reduced to three distinct “predatory” interactions: (1) ALLD
rapidly replacing ALLC; (2) cooperating groups of TFT gradually replacing
ALLD; (3) ALLC slowly replacing TFT because of the latter’s inability to
correct mistakes when playing against itself at non-zero error rates. This
allows the design of a mean-field model Szabó  _et al._ (2000) capturing this
three-species predator-prey dynamics.

##  III Mean-field model

###  III.1 A 3-species predator-prey model

With the three primary predation interactions identified, we can construct a
model of the spatial IPD in the mean-field (MF) limit. The time evolution for
the population fractions $C$, $T$, and $D$ for the three strategies ALLC, TFT
and ALLD take the form:

$\displaystyle\frac{{\rm d}C}{{\rm d}t}$ $\displaystyle=$ $\displaystyle\alpha
CT-\gamma CD\leavevmode\nobreak\ ,$ (1) $\displaystyle\frac{{\rm d}T}{{\rm
d}t}$ $\displaystyle=$ $\displaystyle-\alpha CT+\beta TD\leavevmode\nobreak\
,$ (2) $\displaystyle\frac{{\rm d}D}{{\rm d}t}$ $\displaystyle=$
$\displaystyle\gamma CD-\beta TD\leavevmode\nobreak\ .$ (3)

with $0\leq C,T,D\leq 1$. Here $\alpha$ is the predation rate of TFT by ALLC,
$\beta$ the predation rate of ALLD by TFT, and $\gamma$ that of ALLC by ALLD.
This triangular predatory dynamics is in fact akin to that of the Rock-Paper-
Scissor game Szabo and Fath (2007); Szolnoki  _et al._ (2014).

Given the constraint $C+T+D=1$ one of the above ODEs is redundant and the
above system can be reduced to:

$\displaystyle\frac{{\rm d}C}{{\rm d}t}$ $\displaystyle=$ $\displaystyle\alpha
CT-\gamma C(1-C-T)\leavevmode\nobreak\ ,$ (4) $\displaystyle\frac{{\rm
d}T}{{\rm d}t}$ $\displaystyle=$ $\displaystyle-\alpha CT+\beta
T(1-C-T)\leavevmode\nobreak\ .$ (5)

The equilibrium solution (${\rm d}/{\rm d}t\equiv 0$) is:

$C_{\rm eq}={\beta\over s}\leavevmode\nobreak\ ,\qquad T_{\rm eq}={\gamma\over
s}\leavevmode\nobreak\ ,\qquad D_{\rm eq}={\alpha\over s}\leavevmode\nobreak\
.$ (6)

with $s\equiv\alpha+\beta+\gamma$. Figure 3 shows three representative
solutions to Eqs. (4)–(5), plotted as trajectories in the $[C,T]$ planes (the
corresponding 3D trajectories are restricted to the wedge defined by the
$C+T+D=1$ and $C=0,T=0$, and $D=0$ planes).The three solutions plotted differ
only in their initial conditions, as color-coded.

Figure 3: Representative solutions of Eqs. (4)–(5) in the $[C,T]$ plane, for
$\alpha=0.1$, $\beta=0.5$, et $\gamma=1$ and three different initial
conditions, as color-coded. The equilibrium solution is indicated by the black
dot. The trajectories describe temporal oscillations of the population
fractions qualitatively resembling oscillations observed in classical
predator-prey dynamics or autocatalytic sets of chemical reactions, although
both dynamics remain distinct (see text).

These trajectories described temporal oscillations of the population
fractions, qualitatively resembling oscillations observed in classical
predator-prey (e.g., Lotka-Volterra) dynamics (Song  _et al._ , 2015) or in
autocatalytic (e.g., Belousov-Zhabotinsky) sets of chemical reactions
(Kondepudi and Prigogine, 2015).

However, the resemblance is only superficial. Here we have no attractors in
the form of fixed points or limit cycles. Instead, for a given set of reaction
rates, an infinity of stable orbits are possible, and entirely determined by
the initial conditions through the invariant $C^{\beta}T^{\gamma}D^{\alpha}=$
constant, which can be obtained by algebraic manipulations of Eqs. (1)–(3).
The distinction hinges on the fact that predator-prey and autocatalytic
systems are open systems, with external reservoirs acting as sources and their
magnitude determining the equilibrium solutions and character of limit cycles;
while our system is closed and strictly conservative in its population size,
and its equilibrium solutions are entirely determined by the predations rates.

Further insight into these differences can be obtained through the
linearization of Eqs. (4)–(5) about the equilibrium solutions (6). This
yields, to first order:

$\displaystyle\frac{{\rm d}C_{1}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle\left[{\beta\gamma\over
s}\right]C_{1}+\left[{(\alpha+\gamma)\beta\over
s}\right]T_{1}\leavevmode\nobreak\ ,$ (7) $\displaystyle\frac{{\rm
d}T_{1}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle-\left[{(\alpha+\beta)\gamma\over
s}\right]C_{1}-\left[{\beta\gamma\over s}\right]T_{1}\leavevmode\nobreak\ ,$
(8)

with $C_{1}(t)$ and $T_{1}(t)$ the first order perturbations. These two
coupled linear ODEs can be combined in the decoupled harmonic equations:

$\left(\frac{{\rm d}^{2}}{{\rm
d}t^{2}}+\frac{\alpha\beta\gamma}{s}\right)\begin{pmatrix}C_{1}\\\
T_{1}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix}.$ (9)

These describe harmonic oscillation about the equilibrium solution, with
angular frequency $\omega=(\alpha\beta\gamma/s)^{1/2}$, and constant (small)
amplitudes set by the initial condition. Our system is thus dynamically closer
to an undamped nonlinear pendulum than it is to classical predator-prey
systems, or to autocatalytic chemical reaction systems.

Two of the three trajectories plotted on Fig. 3 spend a significant fraction
of their time running closely along the $C+T=1$ line, implying that $D\to 0$
there. Although stochastic extinction cannot take place in this model, we do
observe it in the discrete spatial IPD (CA model) described in §II. If we set
$D=0$, the system (4)–(5) reduces to the single ODE:

${{\rm d}C\over{\rm d}t}=\alpha C(1-C)\leavevmode\nobreak\ ,$ (10)

with $T=1-C$. This describes a slow drift towards the new equilibrium solution
$[C_{\rm eq},T_{\rm eq}]=[1,0]$ on the timescale $\alpha^{-1}$, as observed in
the IPD cellular automaton simulations when errors in game-playing occur,
following stochastic extinction of ALLD. Likewise, stochastic extinction of
ALLC would lead to a drift towards $[D_{\rm eq},T_{\rm eq}]=[0,1]$ on the
timescale $\beta^{-1}$.

###  III.2 Estimating predation rates from CA simulation

Using the inverse predation rate $\gamma^{-1}$ as a time unit, the equilibrium
solutions (6) can be recast as:

$\alpha={1\over T_{\rm eq}}(1-C_{\rm eq}-T_{\rm eq})\leavevmode\nobreak\
,\qquad\beta={C_{\rm eq}\over T_{\rm eq}}\leavevmode\nobreak\ ,$ (11)

allowing to estimate the dependence of predation rates on game playing error
rates in the cellular automaton IPD simulations. At low error rates
($-5\leq\log_{10}p\leq-2$), this yields the hierarchy $\alpha,\beta\leq 1$.
This no longer holds for regimes where the system enters absorbing states
(i.e., one or more strategies goes extinct), for instance regions where error
rates approach zero ($\log_{10}p\leq-5$) or where players have a high
probability of making mistakes ($\log_{10}p\geq-1$). The equivalency between
the error rate $p$ for the CA model and parameters $\alpha,\beta$ in the MF
model is shown on Fig. 4, calculated for 100 sampled values of the error rate
$\log_{10}p\in[-6,-1]$ with statistical ensembles over 10 realizations of the
simulation which were carried out with random initial conditions for each
value of $p$.

Figure 4: (a) Final population fractions for strategies ALLC (green), ALLD
(red) and TFT (blue) playing spatial IPD games of varying error rates which
were simulated using the cellular automaton model. A set of 100 values of the
error rate were sampled in the interval $p\in[10^{-6},0.5]$. Error bars
represent standard deviations over each statistical ensemble of 10
realizations of the model. Background hues indicate the strategy that
dominates for a given error rate. (b) Predation rates used as parameters for
the mean-field model, derived for each error rate from population fractions
shown in panel (a) using Equation (11).

###  III.3 Spatial extension and diffusion

To investigate the effects of diffusion, we introduce two-dimensional
spatiality and isotropic linear diffusion in the MF model; Equations (4)–(5)
are replaced by:

$\displaystyle\frac{\partial C}{\partial t}$ $\displaystyle=$
$\displaystyle\alpha CT-\gamma C(1-C-T)+D_{C}\nabla^{2}C\leavevmode\nobreak\
,$ (12) $\displaystyle\frac{\partial T}{\partial t}$ $\displaystyle=$
$\displaystyle-\alpha CT+\beta T(1-C-T)+D_{T}\nabla^{2}T\leavevmode\nobreak\
.$ (13)

with $C(x,y,t)$ and $T(x,y,t)$, and the diffusion coefficients $D_{C},D_{T}$
need not be identical.

We solve Eqs. (12)–(13) on the unit square with periodic boundary conditions
in $x$ and $y$, and using the inverse rate $\gamma^{-1}$ as a time unit. We
use an operator splitting scheme, whereby the nonlinear terms are advanced via
fixed-step fourth order Runge-Kutta, and the diffusion terms by the so-called
FTCS scheme, namely a time-explicit Euler method used jointly with a centered
finite difference discretization of the Laplacian terms. Each node of the
spatial mesh is initialized with random values of the two population
fractions, subject to the constraint $C+T<1$.

Figure 5 shows a representative set of results, obtained for predation rates
$\alpha=0.2$, $\beta=0.5$, and $\gamma=1.0$. The mosaic displays solutions for
increasing $D_{T}$ (horizontal) and $D_{C}$ (vertical, and increasing
downwards). The thick white dashed line marks the diagonal $D_{T}=D_{C}$. All
snapshots are extracted after an elapsed time $t=280\,\gamma^{-1}$, and are
rendered as a RGB coding of the population fraction, with blue$\equiv T$
(strategy TFT], green$\equiv C$ (ALLC), and red$\equiv 1-C-T$ (ALLD).

Figure 5: Top panels. Spatial patterns formed in numerical solutions of Eqs.
(12)–(13) in the mean-field model with predation rates $\alpha=0.2$,
$\beta=0.5$, and $\gamma=1$, which are equivalent to an error rate of
$\log_{10}p\simeq-5$. Solutions are integrated to $t=280\,\gamma^{-1}$,
starting from a random initial condition. The diffusion coefficients increases
logarithmically moving rightward and downward. Above the upper white diagonal
($D_{T}\geq 4.5D_{C}$) spatial patterns remain static after nonlinear
saturation. When $D_{\rm T}<D_{\rm C}$ (below the main diagonal), transient
spatial patterns are progressively smoothed by diffusion. Damped oscillatory
patterns materialize in the narrow region between the two diagonals, and no
patterns of significant amplitude build up along the main diagonal $D_{\rm
T}=D_{\rm C}$. An example of RGB color-coding scheme is shown below the top
panels. See also Appendix A for animated versions of simulations on the main
diagonal, and in the steady and damped regimes.

The striking spatial patterns located above the white dotted diagonal are all
steady, and set in very early in the evolution; Further examination of the
temporal evolution of solutions on and below the $D_{T}=D_{C}$ diagonal also
reveals global, system-wide oscillations of the population frequencies, with
faint spatial patterning damping out on a timescale decreasing as diffusion
becomes more important (moving towards the bottom right in the mosaic). The
solutions on the first superdiagonal, indicated by the white dotted line,
oscillate in time but with a spatially quasi-stationary planform. In all cases
the length scale of the spatial patterns, whether steady or oscillating,
decreases as diffusion diminishes (moving towards the upper left in the
mosaic), as one would intuitively expect.

These pattern-forming characteristics, and their dependence on the diffusion
coefficients, are similar for other values of the predation rates $\alpha$ and
$\beta$, but with the populations fractions of the three strategies (setting
the “color” in our RGB coding scheme) varying significantly. Steady patterns
only disappear in the extreme regime $\alpha\ll\beta\ll\gamma$, which occurs
only at very high error rates in the IPD simulations ($\log_{10}p\geq-1$, viz.
Fig. 4). It is well known from studies of pattern formation in classical
nonlinear reaction-diﬀusion systems à la Turing that the inhibition variable
must diﬀuse faster than the excitation variable for steady patterns to form.
Our necessary condition $D_{T}>D_{C}$ for stationary patterns similarly
implies that the strategy TFT acts here as the inhibitor in our 3-species
predator-prey-like system.

###  III.4 Stability analysis

The striking patterns on Fig. 5 can be understood through a linear stability
analysis. Working in Cartesian coordinates in two spatial dimensions, the
first order linearized versions of Eqs. (12)–(13) about the spatially uniform
equilibrium solutions become:

$\displaystyle\frac{\partial C_{1}}{\partial t}$ $\displaystyle=$
$\displaystyle\left[{\beta\gamma\over
s}\right]C_{1}+\left[{(\alpha+\gamma)\beta\over
s}\right]T_{1}+D_{C}\nabla^{2}C_{1}\leavevmode\nobreak\ ,$ (14)
$\displaystyle\frac{\partial T_{1}}{\partial t}$ $\displaystyle=$
$\displaystyle-\left[{(\alpha+\beta)\gamma\over
s}\right]C_{1}-\left[{\beta\gamma\over
s}\right]T_{1}+D_{T}\nabla^{2}T_{1}\leavevmode\nobreak\ .$ (15)

We are seeking solutions describing exponentially growing but spatially steady
planforms:

$\displaystyle\begin{pmatrix}C_{1}({\bf x},t)\\\ T_{1}({\bf
x},t)\end{pmatrix}=\begin{pmatrix}C_{\ast}\\\ T_{\ast}\end{pmatrix}\exp(i{\bf
k\cdot x}+\lambda t)\leavevmode\nobreak\ ,$ (16)

with $\lambda$ real and positive. Substitution in Eqs. (14)–(15) yields two
coupled nonlinear algebraic equations, whose solution yields the dispersion
relation:

$\displaystyle\lambda^{2}$ $\displaystyle+$
$\displaystyle\left[k^{2}(D_{C}+D_{T})\right]\lambda$ $\displaystyle+$
$\displaystyle\left[D_{C}D_{T}k^{4}+{\beta\gamma\over
s}(D_{C}-D_{T})k^{2}+{\alpha\beta\gamma\over s}\right]=0,$

with $k^{2}\equiv k_{x}^{2}+k_{y}^{2}$. Recalling that
$\alpha,\beta,\gamma,D_{T}$ and $D_{C}$ are all positive quantities, a
positive $\lambda$ demands that the positive root of this quadratic polynomial
be retained, and for the later to be purely real further requires

$D_{C}D_{T}k^{4}+{\beta\gamma\over
s}(D_{C}-D_{T})k^{2}+{\alpha\beta\gamma\over s}<0\leavevmode\nobreak\ .$ (18)

This is a convex quadratic polynomial in $k^{2}$, with its extremum at

$k^{2}_{\rm max}={\beta\gamma(D_{T}-D_{C})\over
2sD_{C}D_{T}}\leavevmode\nobreak\ .$ (19)

This corresponds to the wavenumber of the most rapidly growing mode. Note that
the requirement that $k_{\rm max}$ be a real quantity demands $D_{T}>D_{C}$.
Substituting $k_{\rm max}$ back into Eq. (18) leads to the further constraint:

$d^{2}-2(1+A)d+1>0\leavevmode\nobreak\ ,$ (20)

where we have defined $d=D_{T}/D_{C}$ and $A=2\alpha s/\beta\gamma$. The
positive root of this quadratic polynomial yields the minimal diffusion
coefficient ratio $d_{m}$ above which the formation of steady spatial pattern
is possible; for the predation rates used to generate Fig. 5, $d_{m}=4.5$
(dotted diagonal).

The dispersion relation (III.4) also allows oscillations ($\lambda$ picking up
an imaginary part) provided its discriminant is negative:

$(D_{C}-D_{T})^{2}k^{4}-{4\beta\gamma\over
s}(D_{C}-D_{T})k^{2}-{4\alpha\beta\gamma\over s}<0\leavevmode\nobreak\ .$ (21)

This is a quadratic in $k^{2}$, with roots:

$k^{2}(\pm)={2\beta\gamma\over s(D_{C}-D_{T})}\left[1\pm\sqrt{1+{\alpha
s\over\beta\gamma}}\right].$ (22)

Either positive $k^{2}(+)$ or $k^{2}(-)$ is possible, according to the sign of
$D_{C}-D_{T}$. However, these oscillations are always damped, because the real
part of $\lambda$ is now

${\rm Re}(\lambda)=-\frac{1}{2}(D_{C}+D_{T})k^{2}\leavevmode\nobreak\ .$ (23)

Damping is slowest for the largest length scales (small $k$), and weak
diffusion (small $D_{C}+D_{T}$), as expected for a diffusively damped
oscillation.

Note that the $k\to 0$ limit brings us back to uniform oscillations of the
population fractions, as in §III.1. The oscillation frequency
$\omega\equiv{\rm Im}(\lambda)=\sqrt{\alpha\beta\gamma/s}$ is indeed recovered
by setting $k=0$ in the dispersion relation (III.4). Such global oscillations
are not expected to persist in the nonlinear regime when $D_{T}/D_{C}>d(+)$,
because they then become overwhelmed by the exponentially growing spatially
stationary planform.

Overall, the linear stability analysis thus recovers quite well the pattern-
forming characteristics obtained from numerical solutions of the nonlinear
system (12)–(13). We now turn to the pattern-forming behavior of the spatial
IPD model of §II, upon the introduction of diffusion.

##  IV Diffusion in the cellular automaton model

###  IV.1 Pattern formation

Our next step is to investigate whether the structures identified in the mean-
field model (c.f., Fig. 5) can also arise at the microscale. Exploring the
effect of diffusion in the CA model can be carried out using diffusion by
interchange: for each iteration $t_{i}$ of the model, strategies perform a
stochastic random walk on the lattice, exchanging place with another strategy
in their von Neumann neighborhood (i.e., the four players in each cardinal
direction) which is the lattice-based equivalent process of continuous
classical linear diffusion (see, e.g., (Vallis, 2017)). While the CA score-
based update rules remain purely deterministic, the model now incorporates
stochastic elements in both game playing errors and through interchange
diffusion.

We analyzed the parameter space of the CA model under variation of interchange
frequency specified for individual species: Figure 6 shows the final lattice
state for CA simulations where each player on the grid has an interchange
frequency $\nu$. The horizontal axis represents the interchange frequency
defined for TFT ($\nu_{T}$), while the vertical axis displays the frequency
for ALLC ($\nu_{C}$). No distinctive behaviors were observed in the CA model
for regions of the parameter space where the interchange frequency for ALLD
($\nu_{D}$) was independently varied.

Figure 6: Final lattice state for simulations of the spatial IPD using the
cellular automaton model, which includes diffusion and an error rate set at
$\log_{10}p=-2$. Simulations were carried out over $T=500$ iterations on a
lattice of size $L=128$, with games of $M=2000$ moves. The horizontal axis
shows logarithmic variations of TFT’s interchange coefficient while the
vertical axis shows that of ALLC players, both ranging from $\log_{10}\nu=-4$
(negligible diffusion) to $\log_{10}\nu=2$ (close to the well-mixed limit).
The yellow contour lines shows the degree of non-stationarity, where values
approaching zero indicate simulations that become stationary over time whereas
higher values those remaining non-stationary. In the upper-left region, ALLC
players form stable cooperative clusters surrounded by ALLD strategies (see,
e.g., Fig. 2c). See also Appendix A for animated versions of systems on the
main ($\nu_{T}=\nu_{C}$), upper ($\nu_{T}\gg\nu_{C}$) and lower
($\nu_{T}\ll\nu_{C}$) diagonals.

The interchange frequencies progress logarithmically in both cases where the
minimum corresponds to negligible diffusion ($\log_{10}\nu=-4$), where each
player will be switching place with one of its neighbor once every $10^{4}$
iterations. Put differently, in this regime where diffusion is minimal, one
player will change place every iteration on a $128\times 128$ lattice on
average. Maximal values shown for the interchange frequencies approach the
limit where the system is well-mixed ($\log_{10}\nu=2$), as each player will
be switching place $10^{2}$ times at each iteration.

As both interchange coefficients are increased, (i.e., by moving rightward or
downward) larger structures are revealed, whose characteristic length
approaches that of the whole lattice. As interchange frequencies increase
further, the system transitions to an absorbing state where TFT occupies the
whole lattice. Generally, increasing $\nu$ results in the formation of
patterns at larger spatial scales, however moving along the
$\nu_{T}\gg\nu_{C}$ and $\nu_{T}\ll\nu_{C}$ diagonals do not reveal
symmetrical structures. Moreover, final population fractions between these two
cases show significant differences (see discussion in §IV.4).

Thus defined, diffusion by interchange necessarily implies coupled diffusion
frequencies, as interchange occurs with frequency $\nu_{i}$ between a player
belonging to the species $i$ and another random player $j$ whose strategy can
be different. This coupling of the diffusion frequencies between species also
implies that the effective $\nu_{T}$ and $\nu_{C}$ depend on population
frequencies for all three species. Since population frequencies also vary in
time, interchange frequencies are also affected when the system enters
absorbing states (e.g., when a strategy undergoes stochastic extinction).

###  IV.2 Spatiotemporal variability

Over time, simulations in the CA model can either stabilize or continue to
display diverse degrees of spatiotemporal variability. Letting
$\mathbf{G}^{t}$ be the matrix representing strategies on the lattice at
iteration $t$, and $g_{ij}^{t}$ its values, then the degree of non-
stationarity $n$ at time $t$ of a CA simulation corresponds to the number of
sites on the lattice where strategy replacement has taken place since
iteration $t-1$,

$n^{t}=\frac{1}{L^{2}}\sum_{\mathbf{G}^{t}}\Big{[}1-\delta\left(g_{ij^{\prime}}^{t},g_{ij}^{t-1}\right)\Big{]},$
(24)

with $\delta(\alpha,\alpha^{\prime})$ being the Kronecker delta and $n$
normalized by the size of the lattice. Lower non-stationarity thus implies
simulations that become stationary over time while higher non-stationarity
indicates simulations that continue to evolve, with $0\leq n^{t}\leq 1$. The
contour lines on Fig. 6 show the degree of non-stationarity $n^{t}$, averaged
over statistical ensembles of 10 realizations of the simulation with random
initial conditions for each set of parameters.

Varying $\nu_{T}$ and $\nu_{C}$ defines several regimes across which system
behavior ranges from the formation of small-scale clusters in static
equilibrium (top-left region, Fig. 6) to structures whose characteristic
length scale are of the same order of magnitude as the whole lattice (regions
bordering blue panels) which are highly dynamic. Increasing further the
interchange frequencies ultimately gives way to a sudden transition towards
absorbing states where TFT fills the whole lattice (rightmost and bottom
panels) after which the simulation becomes trivially stationary.

###  IV.3 Equivalence of the models

We now examine how pattern formation correlates between the mean-field model
and the cellular automaton model incorporating diffusion. Inspection of Fig. 5
and Fig. 6 suggests that moving along both diagonal axes towards increased
diffusion results in spatial patterns with increasing characteristic scale.
Here, we employ 2D Fourier transforms to formalize this relationship: by
calculating the power spectra of the lattices resulting from the execution of
simulations with both models, we can quantitatively assess how this increase
in wavelength is related to the increase in diffusion.

Utilizing this methodology, we have analyzed the spectra for simulations with
parameters identical to the first upper-diagonal of the MF model (Fig. 5) and
the main diagonal of the CA model (Fig. 6). The 2D Fourier transform of these
spectra, averaged over statistical ensembles of 10 independent simulations of
the models, and collapsed on a 1-dimensional axis, are displayed in Fig. 7a,
with the horizontal axis indicating the wavenumber, excluding the simulations
where TFT invades the entire lattice in the CA model (since the Fourier
spectra would be zero). The spectra for both the MF model (displayed in red)
and the dynamic regime of the CA model where simulations remain non-stationary
(displayed in green) is well fit by a lognormal distribution, while spectra
for the static regime of the CA model are found to align with Gaussian
distributions (spectra displayed in blue). Parameters for the best-fit
distribution are represented in Fig. 7b: the mean $\mu$ of the distributions
averaged over the statistical ensemble are shown, with error bars indicating
the standard deviation over the ensemble.

In Fig. 7c, we present the non-stationarity, computed using Equation (24), at
the final iteration of the model. This is averaged over the statistical
ensembles of each parameter set, effectively distinguishing between the static
regime (where simulations eventually become stationary) and the dynamical
regime (where simulations remain non-stationary).

Figure 7: Fourier analysis was performed on the lattices for both the MF and
CA models, using statistical ensembles of ten simulations conducted with
identical parameters as those used in the simulations on the first super-
diagonal (MF model) and the main diagonal (CA model). (a) The Fourier spectra,
averaged over the statistical ensemble and best-fit distributions, are plotted
against radial wavenumber $k$ ($k^{2}=k_{x}^{2}+k_{y}^{2}$) increasing to the
right. The spectra for the MF model are depicted in red, while those for the
CA model are shown in blue (static regime) and green (dynamic regime). (b) The
mean value of the best-fit distributions, averaged over the statistical
ensemble, is displayed, with error bars representing the standard deviation.
The left side (blue curve) represents the static regime, while the right side
(red and green curves) represents the dynamic regime. All three curves
demonstrate that for both systems the wavenumber decreases in every regime as
diffusion increases. (c) The non-stationarity $n$ is displayed in panel (c).
As diffusion increases, the final latice states gradually transit from static
to non-stationary ($n$ increasing).

This analysis of the underlying trends in pattern formation across both models
demonstrates that the behavior of the MF model reflects that of both the
static and dynamic regimes of the CA model: an increase in diffusion results
in a decrease of the wavenumber of the patterns, which mirrors what could be
inferred by examining Fig. 5 and Fig. 6 as non-stationarity increases.

###  IV.4 Dominance of strategies

Increasing either diffusion frequency ($\nu_{T}$ or $\nu_{C}$), which
corresponds to moving rightward or downward in Fig. 6, triggers regime changes
in the CA model’s behavior. Moving along both axis induces the formation of
large scale structures before well-mixed regimes take over (rightmost and
bottom panels in Fig. 6). In these well-mixed regimes, the stochastic
extinction of one strategy can occur: unlike in the MF model, where
populations are continuous values that can come close to but never actually
reach zero, the populations in the CA model are discrete. Consequently, if a
strategy’s population reaches zero, that species becomes extinct for the
remainder of the simulation, and one of the two remaining strategies will be
further eliminated, as outlined in §III.1. Careful inspection moreover reveals
that upon moving rightward or downward on Fig. 6, the system undergoes
configuration changes that are not symmetrical about the main diagonal: in
these panels, an increase in $\nu_{T}$ results in a more rapid rise of
cooperator populations (ALLC and TFT), contrasting with increases in $\nu_{C}$
where populations of parasites (ALLD) otherwise tend to dominate.

Investigating further this dissimilarity across error rates, Fig. 8 displays
resulting population fractions represented as RGB color coding, each value
calculated from statistical ensembles of 10 realizations with random initial
conditions, with the vertical axis showing the error rate and the horizontal
axis representing the diffusion frequency. In panel (a), the horizontal axis
represents $\nu_{C}$, while $\nu_{T}$ is three orders of magnitude smaller
(i.e., $\nu_{T}=\nu_{C}\times 10^{-3}$). In contrast, panel (b) represents
$\nu_{T}$ on its horizontal axis, with $\nu_{C}$ three orders of magnitude
smaller. Slicing horizontally across panels (a) and (b) of Fig. 8 at
$\log_{10}p=-2$ corresponds to moving along the lower and upper dotted white
diagonals in Fig. 6, respectively. Contour lines in Fig. 8 again represent the
degree of non-stationarity, as calculated via Equation (24) and averaged over
each 10-member ensemble.

Figure 8 illustrates the contrasting effects that varying either interchange
frequency and large-scale structure formation have on resulting population
fractions. Panels (a), where $\nu_{T}\ll\nu_{C}$, and (b), where
$\nu_{C}\gg\nu_{T}$, show crucial differences: while at low interchange
frequencies $\nu_{C}$-dominated diffusion favors TFT more explicitly, moving
rightward on both panels reveals that cooperative strategies (ALLC and TFT)
takes hold for much lower values of the interchange frequency—one exception to
this being when the error rate is very high ($\log_{10}p\gtrsim-1.0$). Simply
put, the system reaches an absorbing state where TFT invades the lattice much
earlier when $\nu_{T}$ is increased, compared to the scenario where $\nu_{C}$
is increased. Large-scale pattern formation regions on Fig. 6 yield higher
population fractions for cooperating strategies (ALLC and TFT) in the
corresponding regions on panel (b) for $\nu_{T}$-led diffusion. The
implications of TFT-driven diffusion ($\nu_{T}\gg\nu_{C}$) and ALLC-driven
diffusion ($\nu_{T}\ll\nu_{C}$) thus result in very different simulation
outcomes. In the latter scenario, parasites gain an advantage, whereas in the
former scenario, cooperators are favored at considerably lower interchange
frequencies.

Differences between the two regimes shown in Fig. 8 are however muddled due to
the coupling of interchange frequencies discussed in §IV.1. System dynamics
also differs significantly in parameter regions where the error rate $p$
approaches its maximum, which are shown in the upper regions of both panels on
Fig. 8. Importantly, in the regime where $\nu_{T}\gg\nu_{C}$, defectors (ALLD)
thrive at high error rates coupled with high interchange frequencies—a
behavior which is not observed in the $\nu_{T}\ll\nu_{C}$ scenario.

Figure 8: Final population fractions, shown as RGB color coding, for cellular
automaton simulations of the spatial IPD which includes diffusion with both
varying interchange frequencies and error rates. Population fractions were
averaged over statistical ensembles of 10 realizations of the simulation with
random initial conditions for each set of parameters. Other parameters such as
lattice size, IPD game length and number of iterations are identical to those
of simulations shown previously on Fig. 6. The horizontal axis represents
interchange frequencies, where panel (a) shows the one associated with
strategy ALLC (TFT’s interchange frequency being three orders of magnitude
smaller) and vice-versa for panel (b). The vertical axis displays the error
rate, ranging from a few errors per game ($p=10^{-3}$) to its maximum where
game dynamics is random ($p=0.5$). Contour lines show the degree of non-
stationarity, where values approaching zero indicate systems that become
stationary over time, whereas those with higher values continue to evolve.
Moving along the white dotted lines on both panels (a) and (b) amounts to
moving across lower and upper diagonals on Fig. 6 respectively.

##  V Conclusion

The emergence, promotion and maintenance of cooperation is an essential
feature of various natural phenomena, whether at small or large scale. Our
work sheds light on how the presence of diffusive processes in physical
systems affects the interplay of strategies competing in the spatial iterated
prisoner’s dilemma. We based our analysis on a deterministic cellular
automaton model (Section II), from which we derived a mean-field model in
which we introduced linear diffusion (Section III). We further investigated
the observed pattern formation through linear stability analysis, then used
diffusion by interchange to induce similar dynamics in the deterministic
cellular automaton model to characterize its effect on the dominance of
strategies (Section IV). Our results indicate that the emergence of patterns
whose characteristic length compares with the system’s size—a behavior
observed in other classes of systems such as autocatalytic
models—significantly influences the populations of strategies playing the IPD
and, importantly, is found to favor cooperation.

Our investigations motivates further research in analyzing how diffusion can
enable the rise of cooperation in various natural systems. While certain
classes of systems such as autocatalytic reaction-diffusion systems are known
to display pattern formation, we have shown that similar behavior can emerge
in other diffusive systems, with these dynamics fostering cooperative
behavior. Further work is needed in order to analyze how the various diffusive
processes at work in the context of specific natural systems—such as bacterial
colonies or bird flocks—influence coordinated actions and enable the emergence
of cooperation between individuals.

###### Acknowledgements.

Conceptualization, A.C.-R. and P.C.; methodology, A.C.-R. and P.C.; software,
A.C.-R. and S.Z.-B.; writing—original draft preparation, A.C.-R., S.Z.-B and
P.C.; writing—review and editing, P.C.; funding acquisition, P.C. All authors
have read and agreed to the published version of the manuscript. This research
was funded by the Natural Sciences and Engineering Research Council of Canada
(grant number RGPIN/05278-2018), the Fonds de recherche Nature et Technologies
of Québec (grant number 314488), the Fondation J. Armand Bombardier Excellence
Scholarship and the Google Cloud Research Credits Program (grant number
215904839).

##  Appendix A Supplementary Information

Animations for simulations presented in this article are available at the
following address:

https://github.com/champagnealexandre/DiffusionPattern

## References

  * Turner and Chao (1999) P. E. Turner and L. Chao, Nature 398, 441 (1999).
  * Bohl  _et al._ (2014) K. Bohl, S. Hummert, S. Werner, D. Basanta, A. Deutsch, S. Schuster, G. Theißen, and A. Schroeter, Molecular BioSystems 10, 3066 (2014).
  * Hong  _et al._ (1992) J. I. Hong, Q. Feng, V. Rotello, and J. Rebek, Science (New York, N.Y.) 255, 848 (1992).
  * Sievers and von Kiedrowski (1994) D. Sievers and G. von Kiedrowski, Nature 369, 221 (1994).
  * Lee  _et al._ (1997) D. H. Lee, K. Severin, and M. R. Ghadiri, Current Opinion in Chemical Biology 1, 491 (1997).
  * Dugatkin (1997) L. A. Dugatkin,  _Cooperation among animals: an evolutionary perspective_ (Oxford University Press, New York, 1997).
  * Dugatkin (1988) L. A. Dugatkin, Behavioral Ecology and Sociobiology 23, 395 (1988).
  * Heinsohn and Packer (1995) R. Heinsohn and C. Packer, Science 269, 1260 (1995).
  * Kay  _et al._ (2020) T. Kay, L. Keller, and L. Lehmann, Proceedings of the National Academy of Sciences 117, 28894 (2020).
  * Lombardo (1985) M. P. Lombardo, Science 227, 1363 (1985).
  * Milinski (1987) M. Milinski, Nature 325, 433 (1987).
  * Sagan (1967) L. Sagan, Journal of theoretical biology 14, 225 (1967).
  * Margulis (1981) L. Margulis,  _Symbiosis in Cell Evolution: Life and Its Environment on the Early Earth_ (W.H. Freeman & Company, San Francisco, 1981).
  * Kauffman (1986) S. A. Kauffman, Journal of Theoretical Biology 119, 1 (1986).
  * Hordijk and Steel (2004) W. Hordijk and M. Steel, Journal of Theoretical Biology 227, 451 (2004).
  * Hordijk (2019) W. Hordijk, Biological Theory 14, 224 (2019).
  * Queller (1997) D. C. Queller, The Quarterly Review of Biology 72, 184 (1997).
  * Szathmáry (2015) E. Szathmáry, Proceedings of the National Academy of Sciences 112, 10104 (2015).
  * Massey and Mishra (2018) S. E. Massey and B. Mishra, Journal of The Royal Society Interface 15, 20180429 (2018).
  * Colman and Colman (2013) A. M. Colman and A. M. Colman,  _Game theory and its applications in the social and biological sciences_ (2013).
  * Myerson (2013) R. B. Myerson,  _Game Theory Analysis of Conflict_ (Harvard University Press, Cumberland, 2013).
  * Watson (2016) J. Watson,  _Strategy: an introduction to game theory_ (Norton & Company, New York; London [etc., 2016).
  * Sigmund (1995) K. Sigmund,  _Games of life: explorations in ecology, evolution, and behaviour_ (1995).
  * Keinan  _et al._ (2006) A. Keinan, B. Sandbank, C. C. Hilgetag, I. Meilijson, and E. Ruppin, Artificial Life 12, 333 (2006).
  * Schuster  _et al._ (2008) S. Schuster, J.-U. Kreft, A. Schroeter, and T. Pfeiffer, Journal of Biological Physics 34, 1 (2008).
  * Fudenberg and Tirole (1991) D. Fudenberg and J. Tirole,  _Game theory_ (MIT Press, Cambridge, Mass., 1991).
  * Nowak (2006) M. A. Nowak,  _Evolutionary dynamics: exploring the equations of life_ (Belknap Press of Harvard University Press, Cambridge, Mass, 2006).
  * Trivers (1985) R. Trivers,  _Social Evolution_ (Benjamin-Cummings Pub Co, Menlo Park, Calif, 1985).
  * Delahaye and Mathieu (1995) J.-P. Delahaye and P. Mathieu, “Complex Strategies in the Iterated Prisoner’s Dilemma,” (1995).
  * Jurišić  _et al._ (2012) M. Jurišić, D. Kermek, and M. Konecki, in  _2012 Proceedings of the 35th International Convention MIPRO_ (2012) pp. 1093–1097.
  * Haider (2005) A. Haider, “Using Genetic Algorithms to Develop Strategies for the Prisoners Dilemma,”  (2005).
  * Sandholm and Crites (1996) T. W. Sandholm and R. H. Crites, Biosystems 37, 147 (1996).
  * Axelrod and Hamilton (1981) R. Axelrod and W. D. Hamilton, Science , 11 (1981).
  * Axelrod (1984) R. M. Axelrod,  _The evolution of cooperation_ , rev. ed ed. (Basic Books, New York, 1984).
  * Nowak and Sigmund (1992) M. A. Nowak and K. Sigmund, Nature 355, 250 (1992).
  * Reif (2009) F. Reif,  _Fundamentals of statistical and thermal physics_ (Waveland Press, Long Grove, Illinois, 2009).
  * Szabó and Tőke (1998) G. Szabó and C. Tőke, Physical Review E 58, 69 (1998).
  * Lim  _et al._ (2002) Y. F. Lim, K. Chen, and C. Jayaprakash, Physical Review E 65 (2002), 10.1103/PhysRevE.65.026134.
  * Szabó  _et al._ (2005) G. Szabó, J. Vukov, and A. Szolnoki, Physical Review E 72, 047107 (2005).
  * Chen  _et al._ (2008) X. Chen, F. Fu, and L. Wang, Physical Review E 78 (2008), 10.1103/PhysRevE.78.051120.
  * Szolnoki  _et al._ (2009) A. Szolnoki, J. Vukov, and G. Szabó, Physical Review E 80 (2009), 10.1103/PhysRevE.80.056112.
  * Fant  _et al._ (2023) L. Fant, O. Mazzarisi, E. Panizon, and J. Grilli, Physical Review E 108, L012401 (2023).
  * Szabó  _et al._ (2000) G. Szabó, T. Antal, P. Szabó, and M. Droz, Physical Review E 62, 1095 (2000).
  * Vukov  _et al._ (2006) J. Vukov, G. Szabó, and A. Szolnoki, Physical Review E 73, 067103 (2006).
  * Amaral  _et al._ (2016) M. A. Amaral, L. Wardil, M. Perc, and J. K. L. da Silva, Physical Review E 94, 032317 (2016).
  * Enquist and Leimar (1993) M. Enquist and O. Leimar, Animal Behaviour 45, 747 (1993).
  * Ferrière and Michod (1995) R. Ferrière and R. E. Michod, Proceedings of the Royal Society of London. Series B: Biological Sciences 259, 77 (1995).
  * Aktipis (2004) C. A. Aktipis, Journal of Theoretical Biology 231, 249 (2004).
  * Hamilton and Taborsky (2005) I. M. Hamilton and M. Taborsky, Proceedings: Biological Sciences 272, 2259 (2005).
  * Le Galliard  _et al._ (2005) J. Le Galliard, R. Ferrière, and U. Dieckmann, The American Naturalist 165, 206 (2005).
  * Helbing and Yu (2009) D. Helbing and W. Yu, Proceedings of the National Academy of Sciences 106, 3680 (2009).
  * Meloni  _et al._ (2009) S. Meloni, A. Buscarino, L. Fortuna, M. Frasca, J. Gomez-Gardenes, V. Latora, and Y. Moreno, Physical Review E 79, 067101 (2009).
  * Smaldino and Schank (2012) P. E. Smaldino and J. C. Schank, Theoretical Population Biology 82, 48 (2012).
  * Yang and Zhang (2023) Z. Yang and L. Zhang, Chaos: An Interdisciplinary Journal of Nonlinear Science 33, 043126 (2023).
  * Ben-Jacob  _et al._ (2000) E. Ben-Jacob, I. Cohen, and H. Levine, Advances in Physics 49, 395 (2000).
  * Bryant and Machta (2023) S. J. Bryant and B. B. Machta, Physical Review Letters 131, 068401 (2023).
  * Hölldobler and Wilson (1990) B. Hölldobler and E. O. Wilson,  _The ants_ (Springer, Berlin, 1990).
  * Cattivelli and Sayed (2011) F. S. Cattivelli and A. H. Sayed, IEEE Transactions on Signal Processing 59, 2038 (2011).
  * Turing (1952) A. M. Turing, Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237, 37 (1952).
  * Nagel and Paczuski (1995) K. Nagel and M. Paczuski, Physical Review E 51, 2909 (1995).
  * von Schantz and Ehtamo (2015) A. von Schantz and H. Ehtamo, Physical Review E 92, 052805 (2015).
  * Yanagisawa (2016) D. Yanagisawa, Collective Dynamics 1, 1 (2016).
  * Ferriere and Michod (1996) R. Ferriere and R. E. Michod, The American Naturalist 147, 692 (1996).
  * Vainstein  _et al._ (2007) M. H. Vainstein, A. T. C. Silva, and J. J. Arenzon, Journal of Theoretical Biology 244, 722 (2007).
  * Sicardi  _et al._ (2009) E. A. Sicardi, H. Fort, M. H. Vainstein, and J. J. Arenzon, Journal of Theoretical Biology 256, 240 (2009).
  * Droz  _et al._ (2009) M. Droz, J. Szwabiński, and G. Szabó, The European Physical Journal B 71, 579 (2009).
  * Cheng  _et al._ (2010) H. Cheng, H. Li, Q. Dai, Y. Zhu, and J. Yang, New Journal of Physics 12, 123014 (2010).
  * Suzuki and Kimura (2011) S. Suzuki and H. Kimura, Journal of Theoretical Biology 287, 42 (2011).
  * Yang and Wang (2011) H. Yang and B. Wang, Chinese Science Bulletin 56, 3693 (2011).
  * Gelimson  _et al._ (2013) A. Gelimson, J. Cremer, and E. Frey, Physical Review E 87, 042711 (2013).
  * Vainstein  _et al._ (2014) M. H. Vainstein, C. Brito, and J. J. Arenzon, Physical Review E 90, 022132 (2014).
  * Fu  _et al._ (2010) F. Fu, M. A. Nowak, and C. Hauert, Journal of Theoretical Biology 266, 358 (2010).
  * Lindgren and Nordahl (1994) K. Lindgren and M. G. Nordahl, Physica D: Nonlinear Phenomena 75, 292 (1994).
  * Nowak and Sigmund (1999) M. A. Nowak and K. Sigmund, in  _The Geometry of Ecological Interactions: Simplifying Spatial Complexity_ (1999).
  * Fisher (1937) R. A. Fisher, Annals of Eugenics 7, 355 (1937).
  * Flores  _et al._ (2022) L. S. Flores, M. A. Amaral, M. H. Vainstein, and H. C. M. Fernandes, “Cooperation in regular lattices,”  (2022).
  * Nowak  _et al._ (1994) M. A. Nowak, S. Bonhoeffer, and R. M. May, Proceedings of the National Academy of Sciences 91, 4877 (1994).
  * Wang  _et al._ (2012) J. Wang, C. Xia, Y. Wang, S. Ding, and J. Sun, Chinese Science Bulletin 57, 724 (2012).
  * Szabo and Fath (2007) G. Szabo and G. Fath, Physics Reports 446, 97 (2007).
  * Szolnoki  _et al._ (2014) A. Szolnoki, M. Mobilia, L.-L. Jiang, B. Szczesny, A. M. Rucklidge, and M. Perc, Journal of The Royal Society Interface 11, 20140735 (2014).
  * Song  _et al._ (2015) Y. Song, C. S. Gokhale, A. Papkou, H. Schulenburg, and A. Traulsen, BMC Evolutionary Biology 15 (2015), 10.1186/s12862-015-0462-6.
  * Kondepudi and Prigogine (2015) D. Kondepudi and I. Prigogine,  _Modern thermodynamics: from heat engines to dissipative structures_ (2015).
  * Vallis (2017) G. K. Vallis,  _Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation_, 2nd ed. (Cambridge University Press, Cambridge, 2017).