then we recover the original DMDc formulation [34].

These ideas are also valid for kernels that don’t take the form (96) but the
algebra is slightly more involved and we therefore do not report the results
here.

##  Appendix D Connection to exact dynamic mode decomposition

We now elucidate the connection between the present work and exact dynamic
mode decomposition of [4] which was introduced in section 2.1. In particular,
we demonstrate that exact DMD can be viewed as a special case of the present
work when there is no sparsity promotion and the kernel is linear. The linear
kernel is
$k(\boldsymbol{u},\boldsymbol{v})=\boldsymbol{u}^{\ast}\boldsymbol{v}$, so the
implicit feature space is simply
$\boldsymbol{\phi}(\boldsymbol{x})=\boldsymbol{x}$. As such, the full model is
the linear map $\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{L}\boldsymbol{x}$
where $\boldsymbol{L}$ is (60)

$\displaystyle\boldsymbol{L}=\tilde{\boldsymbol{W}}\tilde{\boldsymbol{X}}^{\ast}.$
(101)

In exact DMD there is no sparsity promotion so the dictionary used in our
algorithm is full: $\boldsymbol{\Pi}=\boldsymbol{I}$ and
$\tilde{\boldsymbol{X}}=\boldsymbol{X}$. Moreover, $\tilde{\boldsymbol{W}}$ is
given by (33) so

$\displaystyle\boldsymbol{L}=\boldsymbol{Y}\,\,k(\boldsymbol{X},\boldsymbol{X})^{\dagger}\,\,\boldsymbol{X}^{\ast}.$
(102)

Expanding the kernel yields

$\displaystyle\boldsymbol{L}=\boldsymbol{Y}\left(\boldsymbol{X}^{\ast}\boldsymbol{X}\right)^{\dagger}\boldsymbol{X}^{\ast}=\boldsymbol{Y}\boldsymbol{X}^{\dagger}(\boldsymbol{X}^{\ast})^{\dagger}\boldsymbol{X}^{\ast}=\boldsymbol{Y}\boldsymbol{X}^{\dagger}$
(103)

which is identical to the linear operator $\boldsymbol{A}$ from (13) defined
by exact DMD. In (103) we used the identities for the Moore–Penrose
pseudoinverse
$(\boldsymbol{M}^{\ast}\boldsymbol{M})^{\dagger}=\boldsymbol{M}^{\dagger}(\boldsymbol{M}^{\ast})^{\dagger}$
and
$\boldsymbol{M}^{\dagger}(\boldsymbol{M}^{\ast})^{\dagger}\boldsymbol{M}^{\ast}=\boldsymbol{M}^{\dagger}$
for any matrix $\boldsymbol{M}$.

Similarly, the eigenmodes computed by exact DMD are equivalent to those
defined in lemma 1 in the special case of a linear kernel without sparsity
promotion.

##  Appendix E Sensitivity to noise

All machine learning algorithms must be understood in the context of their
sensitivity to noise. To explore the effects of noise, we applied our learning
framework to noise-contaminated data generated by the Lorenz system from
section 5.1. The data are contaminated with Gaussian noise of magnitude $5\%$
of the variance of the original data; the noisy training data is visualised in
the left panel of figure 11. We use the same parameters as section 5.1, but
the system is now integrated to $t=50$. Unlike section 5.1, we assume that we
only have access to snapshot measurements of $\boldsymbol{x}$ and velocity
measurements $\dot{\boldsymbol{x}}$ are unavailable. Therefore, we approximate
the derivative $\dot{\boldsymbol{x}}$ from the noisy snapshot data with a
total-variation regularisation scheme [124]. Then, we use the algorithm of
section 3 to learn the Lorenz system with a quadratic kernel. The results of
the learned model are illustrated in the middle panel of figure 11, along with
the true local linear model, both evaluated at the equilibrium point
$\overline{\boldsymbol{x}}=\begin{bmatrix}-\sqrt{\beta(\rho-1)}&-\sqrt{\beta(\rho-1)}&\rho-1\end{bmatrix}^{T}$.
The reconstructed trajectory shows good qualitative agreement with the true
model, and the local linear model is a good approximation to the true
linearisation. The accuracy of these approximations usually improves as more
samples are added.

Figure 11:  Implicit learning of the Lorenz system in the presence of noise.
Initial measurements are corrupted by $5\%$ Gaussian noise. The velocity
vector $\dot{\boldsymbol{x}}$ is computed by differentiating the data
measurements using total variation regularised differentiation [124]. The
reconstruction captures the qualitative features of the original Lorenz
system, and also accurately reproduces the local linear model at
$\overline{\boldsymbol{x}}$. The trajectories are colored by the adaptive time
step, with red indicating a smaller time step.  Figure 12: Identifying the
spectrum of Burgers equation in the presence of noise over 20 trials. The
first three eigenvalues are plotted in black circles $\circ$ and their
approximations are plotted in blue crosses $\times$.

We also demonstrate the effect of noise on identifying the spectrum of the
viscous Burgers’ equation from section 5.3. Here, the kinematic viscosity is
$\nu=0.1$ and the equations are integrated to $t=4$. The data are snapshots of
the solution which are then corrupted by Gaussian noise of varying magnitude.
No velocity measurements $\dot{\boldsymbol{x}}$ are used and we do not de-
noise the data. Figure 12 plots the first three (repeated) eigenvalues learned
by the algorithm for 20 trials with a quadratic kernel. The first two
eigenvalues are recovered accurately for small values of noise but the
sensitivity to noise increases for eigenvalues of larger magnitude. Again,
these approximations can be improved by adding more samples, or by a suitable
de-noising.

Further experiments indicate that, in the absence of de-noising, the algorithm
is robust to noise in $\boldsymbol{Y}$ but relatively sensitive to noise in
$\boldsymbol{X}$. These observations can be explained through the linear
algebra and statistics underlying our kernel learning framework. The impact of
noise in $\boldsymbol{X}$ is felt in two forms. Firstly, noise may cause
elements to be included in the dictionary which may otherwise have been
excluded. The dictionary is independent of $\boldsymbol{Y}$ and is therefore
unaffected by noise in $\boldsymbol{Y}$. If the sparsification parameter is
not chosen carefully then noise in $\boldsymbol{X}$ may cause the dictionary
to become dense. Secondly, noise causes errors in the final regression,
whether performed online or in batch. As is typical of least-squares
regressions, our algorithm is an unbiased estimator when the noise is
restricted to $\boldsymbol{Y}$. This is because least-squares methods
implicitly assume that there are no “errors in variables” [77]. This
assumption becomes invalid when $\boldsymbol{X}$ is contaminated by noise. As
such, when noise is present, the naïve pseudoinverse (47) solution may prove
insufficient. For example, a similar issue arise in DMD when there is noise in
$\boldsymbol{X}$, and several approaches have been proposed to mitigate the
effects of noise [30, 31, 32, 33, 35, 36], for example solving the regression
problem with total least squares (TLS) [32]. Experiments with TLS in our
present setting were found to be unsuccessful, in part because the TLS problem
is unstable [125] but also because our problem is nonlinear. In particular,
TLS only guarantees the best solution to
$\operatorname*{arg\rm{}min}_{\boldsymbol{A}}\|\boldsymbol{A}\boldsymbol{X}-\boldsymbol{Y}\|_{F}$
only when the errors in $\boldsymbol{X}$ and $\boldsymbol{Y}$ are column-wise
independent and identically distributed with zero mean and covariance matrix
$\sigma^{2}\boldsymbol{I}$ [126, chapter 8]. Therefore, there are no
guarantees of the effectiveness of TLS for our problem, which is
$\operatorname*{arg\rm{}min}_{\tilde{\boldsymbol{W}}}\|\tilde{\boldsymbol{W}}\,k(\tilde{\boldsymbol{X}},\boldsymbol{X})-\boldsymbol{Y}\|_{F}$.
In the case of nonlinear dynamics, whether with SINDy or our kernel approach,
the noise in $\boldsymbol{X}$ is stretched and transformed through the
nonlinearity, adding nonlinear correlations.

In summary, we have demonstrated that our algorithm can remain effective in
the presence of noise. The Lorenz example in figure 11 indicates that the
algorithm is effective when applied to derivatives computed from noisy data by
total-variation regularisation. Additionally, experiments and theory suggest
that the present framework is insensitive to noise in $\boldsymbol{Y}$ but
moderately sensitive to noise in $\boldsymbol{X}$. In addition to total-
variation regularisation, there are several other methods that could be
deployed to limit the effects of noise. For example, one technique could
combine the KRLS algorithm with a Kalman filter, as explored in [127]. Another
filtering approach would be to use the optimal hard threshold criterion for
singular values of [128]. Fully addressing the challenge of noise is an area
of active ongoing research and will be the focus of future work.

## References

  * [1] S. L. Brunton and J. N. Kutz, Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press, 2019.
  * [2] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” J. Fluid Mech., vol. 656, pp. 5–28, 2010.
  * [3] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson, “Spectral analysis of nonlinear flows,” J. Fluid Mech., vol. 641, pp. 115–127, 2009.
  * [4] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” J. Comput. Dyn., vol. 1, no. 2, pp. 391–421, 2014.
  * [5] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems.  SIAM, 2016.
  * [6] M. O. Williams, I. G. Kevrekidis, C. W. Rowley, O. J. B. Matthew, and O. Williams, “A data-driven approximation of the Koopman operator: extending dynamic mode decomposition,” J. Nonlinear Sci., vol. 25, no. 6, pp. 1307–1346, 2015.
  * [7] M. O. Williams, C. W. Rowley, and I. G. Kevrekidis, “A kernel-based method for data-driven Koopman spectral analysis,” J. Comput. Dyn., vol. 2, no. 2, pp. 247–265, 2015.
  * [8] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proc. Natl. Acad. Sci., vol. 113, no. 15, pp. 3932–3937, 2016\.
  * [9] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Data-driven discovery of partial differential equations,” Sci. Adv., vol. 3, no. 4, p. e1602614, 2017.
  * [10] J. Bongard and H. Lipson, “Automated reverse engineering of nonlinear dynamical systems,” Proceedings of the National Academy of Sciences, vol. 104, no. 24, pp. 9943–9948, 2007.
  * [11] M. Schmidt and H. Lipson, “Distilling free-form natural laws from experimental data,” Science, vol. 324, no. 5923, pp. 81–85, 2009.
  * [12] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” J. Comput. Phys., vol. 378, pp. 686–707, 2019.
  * [13] M. Raissi, A. Yazdani, and G. E. Karniadakis, “Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations,” Science, vol. 367, no. 6481, pp. 1026–1030, 2020.
  * [14] M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho, “Lagrangian neural networks,” arXiv preprint arXiv:2003.04630, 2020\.
  * [15] C. Wehmeyer and F. Noé, “Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics,” J. Chem. Phys., vol. 148, no. 24, p. 241703, 2018.
  * [16] I. Mezić, “Analysis of fluid flows via spectral properties of the Koopman operator,” Ann. Rev. Fluid Mech., vol. 45, no. 1, pp. 357–378, 2013.
  * [17] E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G. Daviller, J. Osth, S. Krajnovic, and R. K. Niven, “Cluster-based reduced-order modelling of a mixing layer,” J. Fluid Mech., vol. 754, pp. 365–414, 2014\.
  * [18] S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C. Schütte, and F. Noé, “Data-Driven Model Reduction and Transfer Operator Approximation,” J. Nonlinear Sci., vol. 28, no. 3, pp. 985–1010, 2018\.
  * [19] S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, “Modern Koopman Theory for Dynamical Systems,” arXiv preprint arXiv:2102.12086, 2021.
  * [20] B. Peherstorfer and K. Willcox, “Data-driven operator inference for nonintrusive projection-based model reduction,” Comput. Methods Appl. Mech. Eng., vol. 306, pp. 196–215, 2016.
  * [21] E. Qian, B. Kramer, B. Peherstorfer, and K. Willcox, “Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems,” Phys. D, vol. 406, p. 132401, 2020.
  * [22] J. L. Lumley, Stochastic Tools in Turbulence.  Academic Press, 1970.
  * [23] K. Taira, S. L. Brunton, S. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: An overview,” AIAA Journal, vol. 55, no. 12, pp. 4013–4041, 2017.
  * [24] A. Towne, O. T. Schmidt, and T. Colonius, “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis,” J. Fluid Mech., vol. 847, pp. 821–867, 2018.
  * [25] K. Taira, M. S. Hemati, S. L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S. T. M. Dawson, and C.-A. Yeh, “Modal Analysis of Fluid Flows: Applications and Outlook,” AIAA Journal, vol. 0, pp. 1–36, 2019.
  * [26] N. Fabbiane, O. Semeraro, S. Bagheri, and D. S. Henningson, “Adaptive and model-based control theory applied to convectively unstable flows,” Applied Mechanics Reviews, vol. 66, no. 6, 2014.
  * [27] S. L. Brunton and B. R. Noack, “Closed-loop turbulence control: Progress and challenges,” Applied Mechanics Reviews, vol. 67, no. 5, 2015.
  * [28] D. Sipp and P. J. Schmid, “Linear closed-loop control of fluid instabilities and noise-induced perturbations: A review of approaches and tools,” Applied Mechanics Reviews, vol. 68, no. 2, 2016.
  * [29] C. W. Rowley and S. T. Dawson, “Model reduction for flow analysis and control,” Ann. Rev. Fluid Mech., vol. 49, pp. 387–417, 2017.
  * [30] S. Bagheri, “Effects of weak noise on oscillating flows: Linking quality factor, Floquet modes, and Koopman spectrum,” Physics of Fluids, vol. 26, no. 9, p. 094104, 2014.
  * [31] S. T. Dawson, M. S. Hemati, M. O. Williams, and C. W. Rowley, “Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition,” Exp. Fluids, vol. 57, no. 3, p. 42, 2016.
  * [32] M. S. Hemati, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, “De-biasing the dynamic mode decomposition for applied Koopman spectral analysis,” TCFD, vol. 31, no. 4, pp. 349–368, 2017.
  * [33] T. Askham and J. N. Kutz, “Variable projection methods for an optimized dynamic mode decomposition,” SIAM J. Appl. Dyn. Syst., vol. 17, no. 1, pp. 380–416, 2018.
  * [34] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” SIAM J. Appl. Dyn. Syst., vol. 15, no. 1, pp. 142–161, 2016.
  * [35] O. Azencot, W. Yin, and A. Bertozzi, “Consistent dynamic mode decomposition,” SIAM J. Appl. Dyn. Syst., vol. 18, no. 3, pp. 1565–1585, 2019.
  * [36] I. Scherl, B. Strom, J. K. Shang, O. Williams, B. L. Polagye, and S. L. Brunton, “Robust principal component analysis for particle image velocimetry,” Physical Review Fluids, vol. 5, no. 054401, 2020.
  * [37] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz, “Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control,” PLoS ONE, vol. 11, no. 2, p. e0150171, 2016.
  * [38] B. Schölkopf and A. J. Smola, Learning with kernels: Support vector machines, regularization, optimization, and beyond.  MIT Press, 2002.
  * [39] J. Shawe-Taylor and N. Cristianini, Kernel Methods for Pattern Analysis.  Cambridge University Press, 2004.
  * [40] J. Bouvrie and B. Hamzi, “Kernel methods for the approximation of nonlinear systems,” SIAM J. Control Optim., vol. 55, no. 4, pp. 2460–2492, 2017\.
  * [41] Y. Kawahara, “Dynamic mode decomposition with reproducing kernels for Koopman spectral analysis,” in Adv. Neural Inf. Process. Syst., pp. 919–927, 2016\.
  * [42] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis, “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 27, no. 10, p. 103111, 2017.
  * [43] H. Owhadi and G. R. Yoo, “Kernel Flows: From learning kernels from data into the abyss,” J. Comput. Phys., vol. 389, pp. 22–47, 2019.
  * [44] D. Giannakis, “Data-driven spectral decomposition and forecasting of ergodic dynamical systems,” Applied and Computational Harmonic Analysis, vol. 47, no. 2, pp. 338–396, 2019.
  * [45] S. Das and D. Giannakis, “Koopman spectra in reproducing kernel Hilbert spaces,” Applied and Computational Harmonic Analysis, vol. 49, no. 2, pp. 573–607, 2020.
  * [46] P. Gelß, S. Klus, I. Schuster, and C. Schütte, “Feature space approximation for kernel-based supervised learning,” tech. rep., 2020.
  * [47] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Neural operator: Graph kernel network for partial differential equations,” arXiv preprint arXiv:2003.03485, 2020.
  * [48] D. Burov, D. Giannakis, K. Manohar, and A. Stuart, “Kernel analog forecasting: Multiscale test problems,” arXiv:2005.06623, 2020.
  * [49] Y. Engel, S. Mannor, and R. Meir, “The kernel recursive least-squares algorithm,” IEEE Trans. Signal Process., vol. 52, no. 8, pp. 2275–2285, 2004.
  * [50] J. C. Loiseau and S. L. Brunton, “Constrained sparse Galerkin regression,” J. Fluid Mech., vol. 838, pp. 42–67, 2018.
  * [51] M. R. Jovanović and B. Bamieh, “Componentwise energy amplification in channel flows,” J. Fluid Mech., vol. 534, pp. 145–183, 2005.
  * [52] B. J. McKeon and A. S. Sharma, “A critical-layer framework for turbulent pipe flow,” J. Fluid Mech., vol. 658, pp. 336–382, 2010.
  * [53] B. J. McKeon, “The engine behind (wall) turbulence: Perspectives on scale interactions,” J. Fluid Mech., vol. 817, p. 1, 2017.
  * [54] M. R. Jovanović, “From bypass transition to flow control and data-driven turbulence modeling: An input–output viewpoint,” Ann. Rev. Fluid Mech., vol. 53, pp. 311–345, 2021.
  * [55] B. Herrmann, P. J. Baddoo, R. Semaan, S. L. Brunton, and B. J. McKeon, “Data-driven resolvent analysis,” J. Fluid Mech., vol. 918, p. A10, 2021\.
  * [56] K. K. Chen, J. H. Tu, and C. W. Rowley, “Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses,” J. Nonlinear Sci., vol. 22, pp. 887–915, dec 2012.
  * [57] J. L. Proctor and P. A. Eckhoff, “Discovering dynamic patterns from infectious disease data using dynamic mode decomposition,” International health, vol. 7, no. 2, pp. 139–145, 2015.
  * [58] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, “Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition,” Journal of Neuroscience Methods, vol. 258, pp. 1–15, 2016.
  * [59] J. Grosek and J. N. Kutz, “Dynamic mode decomposition for real-time background/foreground separation in video,” arXiv preprint arXiv:1404.7592, 2014.
  * [60] N. B. Erichson, S. L. Brunton, and J. N. Kutz, “Compressed dynamic mode decomposition for real-time object detection,” Journal of Real-Time Image Processing, 2016.
  * [61] E. Berger, M. Sastuba, D. Vogt, B. Jung, and H. B. Amor, “Estimation of perturbations in robotic behavior using dynamic mode decomposition,” Journal of Advanced Robotics, vol. 29, no. 5, pp. 331–343, 2015.
  * [62] I. Abraham and T. D. Murphey, “Active learning of dynamics for data-driven control using Koopman operators,” IEEE Transactions on Robotics, vol. 35, no. 5, pp. 1071–1083, 2019.
  * [63] D. Bruder, B. Gillespie, C. David Remy, and R. Vasudevan, “Modeling and control of soft robots using the Koopman operator and model predictive control,” in Robotics: Science and Systems, (Freiburg im Breisgau, June22-26,2019), 2019.
  * [64] X. T. G Mamakoukas, M Castano and T. Murphey, “Local Koopman operators for data-driven control of robotic systems,” in Proceedings of “Robotics: Science and Systems 2019”, Freiburg im Breisgau, June 22-26, 2019, IEEE, 2019\.
  * [65] J. Mann and J. N. Kutz, “Dynamic mode decomposition for financial trading strategies,” Quantitative Finance, pp. 1–13, 2016.
  * [66] Y. Susuki and I. Mezic, “Nonlinear Koopman Modes and Coherency Identification of Coupled Swing Dynamics,” IEEE Transactions on Power Systems, vol. 26, pp. 1894–1904, Nov. 2011.
  * [67] Y. Susuki, I. Mezić, and T. Hikihara, “Coherent Swing Instability of Power Grids,” Journal of Nonlinear Science, vol. 21, pp. 403–439, Feb. 2011.
  * [68] R. Taylor, J. N. Kutz, K. Morgan, and B. Nelson, “Dynamic mode decomposition for plasma diagnostics and validation,” arXiv preprint arXiv:1702.06871, 2017.
  * [69] A. A. Kaptanoglu, K. D. Morgan, C. J. Hansen, and S. L. Brunton, “Characterizing magnetized plasmas with dynamic mode decomposition,” Physics of Plasmas, vol. 27, p. 032108, 2020.
  * [70] M. R. Jovanović, P. J. Schmid, and J. W. Nichols, “Sparsity-promoting dynamic mode decomposition,” Phys. Fluids, vol. 26, no. 2, p. 024103, 2014\.
  * [71] B. R. Noack, W. Stankiewicz, M. Morzynski, and P. J. Schmid, “Recursive dynamic mode decomposition of a transient cylinder wake,” J. Fluid Mech., vol. 809, pp. 843–872, 2016.
  * [72] J. N. Kutz, X. Fu, and S. L. Brunton, “Multi-resolution dynamic mode decomposition,” SIAM J. Appl. Dyn. Syst., vol. 15, no. 2, pp. 713–735, 2016.
  * [73] M. S. Hemati, M. O. Williams, and C. W. Rowley, “Dynamic mode decomposition for large and streaming datasets,” Physics of Fluids (1994-present), vol. 26, no. 11, p. 111701, 2014.
  * [74] S. L. Brunton, J. L. Proctor, J. H. Tu, and J. N. Kutz, “Compressed sensing and dynamic mode decomposition,” Journal of Computational Dynamics, vol. 2, no. 2, pp. 165–191, 2015.
  * [75] F. Gueniat, L. Mathelin, and L. Pastur, “A dynamic mode decomposition approach for large and arbitrarily sampled systems,” Physics of Fluids, vol. 27, no. 2, p. 025113, 2015.
  * [76] J. H. Tu, C. W. Rowley, J. N. Kutz, and J. K. Shang, “Spectral analysis of fluid flows using sub-Nyquist-rate PIV data,” Experiments in Fluids, vol. 55, no. 9, pp. 1–13, 2014.
  * [77] G. H. Golub and C. F. Van Loan, Matrix Computations, vol. 3.  JHU Press, 2013.
  * [78] C. Eckart and G. Young, “The approximation of one matrix by another of lower rank,” Psychometrika, vol. 1, no. 3, pp. 211–218, 1936.
  * [79] W. Liu, P. P. Pokharel, and J. C. Principe, “The kernel least-mean-square algorithm,” IEEE Trans. Signal Process., vol. 56, no. 2, pp. 543–554, 2008.
  * [80] W. Liu, J. C. Príncipe, and S. Haykin, Kernel Adaptive Filtering: A Comprehensive Introduction.  Wiley, 2010.
  * [81] B. Schölkopf, A. Smola, and K. R. Müller, “Nonlinear component analysis as a kernel eigenvalue problem,” Neural Comput., vol. 10, no. 5, pp. 1299–1319, 1998.
  * [82] S. Mika, B. Schölkopf, A. Smola, K. R. Müller, M. Scholz, and G. Rätsch, “Kernel PCA and de-noising in feature spaces,” in Adv. Neural Inf. Process. Syst., pp. 536–542, 1999.
  * [83] R. Herbrich, Learning Kernel Classifiers: Theory and Algorithms.  MIT Press, 2002.
  * [84] G. Camps-Valls and L. Bruzzone, “Kernel-based methods for hyperspectral image classification,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 6, pp. 1351–1362, 2005.
  * [85] N. Cristianini and J. Shawe-Taylor, An introduction to support vector machines and other kernel-based learning methods.  Cambridge University Press, 2000.
  * [86] J. Mercer, “XVI. Functions of positive and negative type, and their connection the theory of integral equations,” Philos. Trans. R. Soc. London, vol. 209, no. 441-458, pp. 415–446, 1909.
  * [87] G. Kimeldorf and G. Wahba, “Some results on Tchebycheffian spline functions,” J. Math. Anal. Appl., vol. 33, no. 1, pp. 82–95, 1971.
  * [88] B. Schölkopf, R. Herbrich, and A. J. Smola, “A generalized representer theorem,” in Lect. Notes Comput. Sci., vol. 2111, pp. 416–426, 2001.
  * [89] J. T. Y. Kwok and I. W. H. Tsang, “The pre-image problem in kernel methods,” IEEE Trans. Neural Networks, vol. 15, no. 6, pp. 1517–1525, 2004.
  * [90] C. K. Williams and M. Seeger, “Using the nyström method to speed up kernel machines,” in Adv. Neural Inf. Process. Syst., 2001.
  * [91] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” NIPS, vol. 3, no. 4, p. 5, 2007.
  * [92] J. R. Bunch and C. P. Nielsen, “Updating the singular value decomposition,” Numer. Math., vol. 31, no. 2, pp. 111–129, 1978.
  * [93] M. Brand, “Fast low-rank modifications of the thin singular value decomposition,” Linear Algebra Appl., vol. 415, no. 1, pp. 20–30, 2006\.
  * [94] P. Benner, P. Goyal, B. Kramer, B. Peherstorfer, and K. Willcox, “Operator inference for non-intrusive model reduction of systems with non-polynomial nonlinear terms,” Comput. Methods Appl. Mech. Eng., vol. 372, 2020.
  * [95] I. Steinwart, D. Hush, and C. Scovel, “An explicit description of the reproducing Kernel Hilbert spaces of Gaussian RBF kernels,” IEEE Trans. Inf. Theory, vol. 52, no. 10, pp. 4635–4643, 2006.
  * [96] B. Haasdonk and H. Burkhardt, “Invariant kernel functions for pattern analysis and machine learning,” Mach. Learn., vol. 68, no. 1, pp. 35–61, 2007.
  * [97] S. Klus, P. Gelß, F. Nüske, and F. Noé, “Symmetric and antisymmetric kernels for machine learning problems in quantum physics and chemistry,” 2021.
  * [98] D. Decoste and B. Schölkopf, “Training invariant support vector machines,” Mach. Learn., vol. 46, no. 1-3, pp. 161–190, 2002.
  * [99] P. Reiterer, C. Lainscsek, F. Schürrer, C. Letellier, and J. Maquet, “A nine-dimensional Lorenz system to study high-dimensional chaos,” J. Phys. A. Math. Gen., vol. 31, no. 34, pp. 7121–7139, 1998.
  * [100] E. N. Lorenz, “Deterministic nonperiodic flow,” J. Atmos. Sci., vol. 20, no. 2, pp. 130–141, 1963.
  * [101] J. A. Acebrón, L. L. Bonilla, C. J. Vicente, F. Ritort, and R. Spigler, “The Kuramoto model: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys., vol. 77, no. 1, pp. 137–185, 2005.
  * [102] P. Gelß, S. Klus, J. Eisert, and C. Schütte, “Multidimensional Approximation of Nonlinear Dynamical Systems,” J. Comput. Nonlinear Dyn., vol. 14, no. 6, 2019.
  * [103] J. Snyder, J. L. Callaham, S. L. Brunton, and J. N. Kutz, “Data-driven stochastic modeling of coarse-grained dynamics with finite-size effects using Langevin regression,” arXiv preprint arXiv:2103.16791, 2021.
  * [104] T. A. Driscoll, N. Hale, and L. N. Trefethen, Chebfun Guide.  Oxford: Pafnuty Publications, 2014.
  * [105] S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems,” J. Comput. Phys., vol. 176, no. 2, pp. 430–455, 2002.
  * [106] H. Montanelli and N. Bootland, “Solving periodic semilinear stiff PDEs in 1D, 2D and 3D with exponential integrators,” Math. Comput. Simul., vol. 178, pp. 307–327, 2020.
  * [107] C. D. Brummitt and J. C. Sprott, “A search for the simplest chaotic partial differential equation,” Phys. Lett. A, vol. 373, no. 31, pp. 2717–2721, 2009.
  * [108] S. Filip, A. Javeed, and L. N. Trefethen, “Smooth random functions, random ODEs, and Gaussian processes,” SIAM Review, vol. 61, no. 1, pp. 185–205, 2019.
  * [109] A. Schwaighofer, V. Tresp, and K. Yu, “Learning Gaussian process kernels via hierarchical Bayes,” in Adv. Neural Inf. Process. Syst., 2005.
  * [110] P. C. Hansen, “The truncated SVD as a method for regularization,” BIT, vol. 27, no. 4, pp. 534–553, 1987.
  * [111] V. Roth, “The generalized LASSO,” IEEE Trans. Neural Networks, vol. 15, no. 1, pp. 16–28, 2004.
  * [112] G. Wang, D. Y. Yeung, and F. H. Lochovsky, “The kernel path in kernelized LASSO,” in J. Mach. Learn. Res., vol. 2, pp. 580–587, 2007.
  * [113] E. Åkervik, L. Brandt, D. S. Henningson, J. Hœpffner, O. Marxen, and P. Schlatter, “Steady solutions of the Navier-Stokes equations by selective frequency damping,” Physics of Fluids, vol. 18, no. 6, p. 068102, 2006.
  * [114] S. Bagheri, E. Åkervik, L. Brandt, and D. S. Henningson, “Matrix-free methods for the stability and control of boundary layers,” AIAA journal, vol. 47, no. 5, pp. 1057–1068, 2009.
  * [115] A. Monokrousos, E. Åkervik, L. Brandt, and D. S. Henningson, “Global three-dimensional optimal disturbances in the Blasius boundary-layer flow using time-steppers,” J. Fluid Mech., vol. 650, pp. 181–214, 2010.
  * [116] J.-C. Loiseau, M. A. Bucci, S. Cherubini, and J.-C. Robinet, “Time-stepping and Krylov methods for large-scale instability problems,” in Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics, pp. 33–73, Springer, 2019.
  * [117] J. H. M. Ribeiro, C.-A. Yeh, and K. Taira, “Randomized resolvent analysis,” Physical Review Fluids, vol. 5, no. 3, p. 033902, 2020.
  * [118] S. Bagheri, P. Schlatter, P. J. Schmid, and D. S. Henningson, “Global stability of a jet in crossflow,” J. Fluid Mech., vol. 624, no. 9, pp. 33–44, 2009.
  * [119] V. Theofilis, “Global linear instability,” Ann. Rev. Fluid Mech., vol. 43, pp. 319–352, 2011.
  * [120] P. Luchini and A. Bottaro, “Adjoint equations in stability analysis,” Ann. Rev. Fluid Mech., vol. 46, 2014.
  * [121] H. Zhang, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, “Online dynamic mode decomposition for time-varying systems,” SIAM J. Appl. Dyn. Syst., vol. 18, no. 3, pp. 1586–1609, 2019.
  * [122] C.-A. Yeh and K. Taira, “Resolvent-analysis-based design of airfoil separation control,” J. Fluid Mech., vol. 867, pp. 572–610, 2019.
  * [123] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit,” Proceedings of the Royal Society of London A, vol. 474, no. 2219, 2018.
  * [124] R. Chartrand, “Numerical Differentiation of Noisy, Nonsmooth Data,” ISRN Appl. Math., vol. 2011, pp. 1–11, 2011.
  * [125] G. H. Golub and C. F. van Loan, “An analysis of the total least squares problem,” SIAM J. Numer. Anal., vol. 17, no. 6, pp. 883–893, 1980.
  * [126] S. V. Huffel and J. Vandewalle, The total least squares problem: computational aspects and analysis.  SIAM, 1991.
  * [127] W. Liu, I. Park, Y. Wang, and J. C. Principe, “Extended kernel recursive least squares algorithm,” IEEE Trans. Signal Process., vol. 57, no. 10, pp. 3801–3814, 2009.
  * [128] M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is $4/\sqrt{3}$,” IEEE Trans. Inf. Theory, vol. 60, no. 8, pp. 5040–5053, 2014.


# Complex networks of interacting stochastic tipping elements:
cooperativity of phase separation in the large-system limit

Jan Kohler  Institute for Theoretical Physics, University of Leipzig, 04103
Leipzig, Germany, EU Earth System Analysis, Potsdam-Institute for Climate
Impact Research, Member of the Leibniz Association,14473 Potsdam, Germany, EU
Nico Wunderling<EMAIL_ADDRESS><EMAIL_ADDRESS>Earth System Analysis, Potsdam-Institute for Climate Impact Research, Member
of the Leibniz Association,14473 Potsdam, Germany, EU Institute of Physics and
Astronomy, University of Potsdam, 14476 Potsdam, Germany, EU Department of
Physics, Humboldt University of Berlin, 12489 Berlin, Germany, EU Jonathan F.
Donges  Earth System Analysis, Potsdam-Institute for Climate Impact Research,
Member of the Leibniz Association,14473 Potsdam, Germany, EU Stockholm
Resilience Centre, Stockholm University, 10691 Stockholm, Sweden, EU Jürgen
Vollmer<EMAIL_ADDRESS><EMAIL_ADDRESS>Institute
for Theoretical Physics, University of Leipzig, 04103 Leipzig, Germany, EU

###### Abstract

Tipping elements in the Earth System receive increased scientific attention
over the recent years due to their nonlinear behavior and the risks of abrupt
state changes. While being stable over a large range of parameters, a tipping
element undergoes a drastic shift in its state upon an additional small
parameter change when close to its tipping point. Recently, the focus of
research broadened towards emergent behavior in networks of tipping elements,
like global tipping cascades triggered by local perturbations. Here, we
analyze the response to the perturbation of a single node in a system that
initially resides in an unstable equilibrium. The evolution is described in
terms of coupled nonlinear equations for the cumulants of the distribution of
the elements. We show that drift terms acting on individual elements and
offsets in the coupling strength are sub-dominant in the limit of large
networks, and we derive an analytical prediction for the evolution of the
expectation (i.e., the first cumulant). It behaves like a single aggregated
tipping element characterized by a dimensionless parameter that accounts for
the network size, its overall connectivity, and the average coupling strength.
The resulting predictions are in excellent agreement with numerical data for
Erdös-Rényi, Barabási-Albert and Watts-Strogatz networks of different size and
with different coupling parameters.

tipping elements, nonlinear dynamics, complex networks, noise induced tipping

##  I Introduction

Hysteresis is a hallmark of first-order phase transitions. For thermodynamic
systems it leads to supercooling and superheating with subsequent explosive
phase changes [1, 2]. These rapid changes are well-understood for common
thermodynamic phase transitions [e.g. 3, 4, 5]. However, they still pose
challenges for systems with non-standard interaction rules, like the
Achlioptas Process [6, 7, 8, 9, 10], or non-standard interaction topologies
like processes on networks, where hysteresis can lead to cascading failure
[11, 12, 13, 14]. Hysteresis and explosive transitions between (meta-)stable
states are also commonly observed in other systems. In ecological and climate
systems [15] and in finance, economics and politics [16] they are commonly
denoted as tipping processes. Further applications are discussed in the recent
reviews [17, 12].

The present study is inspired by current models of climate change that are
formulated in terms of networks of interacting tipping elements [18, 19, 20,
21, 22, 23, 24, 25]. The interactions provide long-lived metastable states and
rapid cooperative transitions between the states. Often such a stabilization
and cooperation in a tipping event is caused by positive feedback effects,
i.e., coupling or interactions that tend to align individual elements. For
ferromagnetic systems the interaction tends to align individual spins with the
molecular field. An example in the Earth System is the surface albedo of sea
ice [26, 27]. A decrease in the ice-covered surface area due to an increase in
global mean temperature decreases the surface albedo. This, in turn, increases
the temperature and causes higher rates of melting [28, 29]. The interaction
between the elements leads to cascading behavior when the abrupt state shift
of an element causes the tipping of another [30]. The resulting rapid
nonlinear changes of the climate have been predicted almost 40 years ago [31].
Early on they were investigated regarding a Snowball-Earth/warm-Earth-state
transition [32, 33, 34, 35], and more recently towards a potential hothouse
state [23, 15].

Figure 1:  Dynamics of the nodes in a network of tipping elements. (a)
Bifurcation diagram of Eq. (1) where solid and dashed red lines show the
positions of the stable and unstable steady states, respectively. In the range
between $r_{\text{crit,\ low}}$ and $r_{\text{crit,\ high}}$ there are two
stable and one unstable state. At the borders of the interval a stable state
merges with the unstable state, and they disappear in a saddle-node
bifurcation, and systems at that position will decay to the remaining stable
state. The trembling blue line illustrates the stochastic fluctuations around
the stable state that we are using here, see Eq. (2a). (b) Nodes in an
exemplary network of interacting tipping elements (see Eq. (2a)), where each
node is of the form shown in panel (a)).

The present paper adds to the ongoing research on the dynamic behavior of
coupled tipping elements by identifying universal behavior in the tipping of
large-scale systems. Tipping in a network comprising a large number of tipping
elements will be characterized by the first two cumulants (expectation and
variance) of the distribution of tipping elements. In large networks the
expectation behaves like a single aggregated tipping element. Moreover, the
variance remains small, and for small noise it rapidly decays except for an
initial transient growth early on. Hence, we adopted the expectation as order
parameter of the transition, and we will show that the tipping dynamics
depends only on two dimensionless control parameters. These parameters take an
analogous rise to the temperature and the external field in the thermodynamic
treatment of magnetic phase transitions [36].

The paper is organized as follows. In Sect. II we introduce the system of
nonlinear differential equations that are analyzed in the present study, and
we describe how they are coupled in network structures. Moreover, we also
describe how the systems are simulated. In Sect. III we present our main
results. First, the cumulants of the distribution of states of a network of
unbiased tipping elements are defined, and their time evolution is derived
(Sect. III.1). Then, we introduce rescaled coordinates that allow us to
discuss universal aspects of the late-time dynamics (Sect. III.2). In Sect. IV
the dynamics of the expectation is worked out and compared to numerical data.
Sect. V addresses the impact of additional forcing terms (bias) and offsets in
the coupling coefficients. In Sect. VI, we discuss the main findings of our
work, its relations to other work, and we suggest follow-up work. Sect. VII
concludes the paper with an exposition of our most important findings.

##  II Theory and Methods

The present study addresses the collective response of $N$ tipping elements
that are coupled linearly through network links distributed according to
paradigmatic network types (Erdös-Rényi, Barabási-Albert, Watts-Strogatz). The
tipping elements are described by a system of differential equations adapted
from Refs. [25, 24] using the software package PyCascades [37]. A Gaussian
white noise is applied to the tipping elements to model stochastic
fluctuations.

###  II.1 Equations of motion

Many natural systems show tipping paired with hysteresis-like behavior [16].
The behavior of such a tipping element is commonly modeled by cubic
differential equations [24, 25, 38] with normal form [39]

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}t}=-x^{3}+x+r$ (1)

where $r$ is the bifurcation parameter. The system is in a bistable state when
its right-hand side has three roots, i.e., for
$r\in(r_{\text{crit,low}},r_{\text{crit,high}})=(-\sqrt{4/27},\sqrt{4/27})$.
At $r_{\text{crit,low}}$ and $r_{\text{crit,high}}$ two of the roots disappear
in a saddle-node bifurcation. The bifurcation diagram of a single node
including noise is shown in Fig. 1a.

Commonly, tipping elements are not isolated. For instance, Krönke  _et al._
[25] suggested to study the dynamics of $N$ coupled tipping elements
$k\in\\{1,\dots,N\\}$ that are linearly coupled to other nodes $l\neq k$ (see
Fig. 1b),

$\displaystyle\frac{\mathrm{d}x_{k}}{\mathrm{d}t}=-x_{k}^{3}+x_{k}+{\sum_{l}}^{\prime}d_{kl}x_{l}+\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}\,.$
(2a) The prime at the sum indicates that $l$ takes values in
$l\in\\{1,\dots,N\\}\backslash\\{k\\}$. On average a node $k$ is coupled to
$pN$ other elements. $W_{k}$ denotes a Wiener process and $\xi$ the strength
of the noise.

Equation (2a) describes homogeneous systems where all nodes follow the same
dynamics and have the same critical parameter values for tipping. Such systems
are discussed in Sect. III and Sect. IV. In Sect. V we expand that discussion
to treat dynamics with additional drift terms, $d_{k}$, and offset terms,
$r_{kl}$,

$\displaystyle\frac{\mathrm{d}x_{k}}{\mathrm{d}t}=-x_{k}^{3}+x_{k}+d_{k}+{\sum_{l}}^{\prime}d_{kl}\bigl{(}x_{l}-r_{kl}\bigr{)}+\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}\,.$
(2b)

###  II.2 Network topologies

There is ample room for different types of interaction networks when one only
specifies that a node is connected to $p\,(N-1)$ other nodes on average. We
will therefore consider three types of networks to explore the impact of the
structure of the interaction network on the tipping dynamics:

ER

The connections form an Erdös-Rényi network [40] when they are assigned
randomly with probability $p$. The probability distribution to be connected to
$m$ nodes, i.e., the degree distribution, amounts to a binomial distribution.

BA

Barabási-Albert networks [41] are built by sequentially adding connections
under the constraint of preferential attachment: elements with a higher number
of connections have a larger chance to receive more links. This leads to
scale-free networks where the degree distribution has a power-law tail that
decays like $m^{-3}$.

WS

Watts-Strogatz [42] networks are built by arranging the elements as a one-
dimensional ring, where each node is connected to its $pN$ nearest neighbors.
Subsequently, a fraction $\beta$ of the links is reconnected at random to a
new element. This generates networks with the small-world property.

###  II.3 Numerical Methods

Tipping behavior in ER, BA, and WS networks has been analyzed previously by
Krönke  _et al._ [25]. Wunderling  _et al._ [37] developed the Python package
PyCascades to create the directed network of coupled tipping elements with the
python networkx2.3 package [43] and analyze the dynamics of tipping cascades
based on the package sdeint [44] for the integration of stochastic
differential equations. These programs are used here to generate the networks
and integrate the stochastic differential equations (2a). If not stated
otherwise the system is integrated up to $t=100$, using a step size of $\Delta
t=0.01$.

The coupling strengths $d_{kl}$ will be drawn from a uniform distribution with
$0<d_{kl}<2d$.

In the simulation runs, the nodes are initially placed at $x=0$, i.e., the
unstable fixed point, and the cascade is started by setting $\nu$ nodes to
$x=1$. Thus, we focus on the transient behavior of tipping.

At times, the noise turns the point $x=0$ into a stable fixed point of the
dynamics [45]. Then, we study cases where the coupling of the system is strong
enough such that all nodes follow the initial displacement of the displaced
node.

##  III Mean-Field Evolution

We characterize the state of the network of tipping elements by the
expectation of the state $x_{k}$ of the elements,

$\displaystyle\overline{x}=\left\langle
x_{k}\right\rangle=\frac{1}{N}\;\sum_{k=1}^{N}x_{k}$ (3a) and by their
variance
$\displaystyle\overline{v}=\left\langle(x_{k}-\overline{x})^{2}\right\rangle=\left\langle
y_{k}^{2}\right\rangle$ (3b)

where $y_{k}=x_{k}-\overline{x}$ denotes the deviation of $x_{k}$ from the
expectation $\overline{x}$. In Sect. III.1 we derive the equations for the
time evolution of $x$ and $v$ that derive from Eq. (2a). To this end we adopt
a closure where we suppress higher-order moments and correlations. Then we
point out universal aspects of the time evolution towards the asymptotic
steady state (Sect. III.2).

Figure 2:  Asymptotic Values of the Evolution. (a) Phase-space flow of the
dynamical system, Eq. (9), for $\sqrt{\alpha}\simeq 1/3$. The analytical
solutions of the nullclines are marked by gray lines, and fixed points are
marked by black circles, respectively. The red line shows the numerical
solution for $N=80$ tipping elements that are connected by an ER network with
$d=0.1$ and $p=0.25$. The numerical solutions for the blue vectors have been
computed using the function matplotlib.pyplot.streamplot in python. (b) The
red and the green line show the time evolution of the expectation and the
variance of the numerical solution shown in panel (a). The thinner blue line
shows the prediction, Eq. (13), for $x(0)=1/N$ (for explanation see Sect.
IV.2). The dotted line shows the same function shifted by $\Delta\tau=1.24$.
The shifted curve provides an excellent description of the late-time evolution
of the expectation. (c) The data points show the asymptotic values $x^{\ast}$
of numerical simulations of Eq. (2a) with $p=0.25$ and couplings $d_{kl}$
drawn from a uniform distribution $d_{kl}\in(0,0.2)$. Error bars indicate the
standard deviation derived from $50$ realizations of the network configuration
and the noise. The connectivity $\lambda_{d}$ is varied by increasing $N$ from
$20$ to $400$ (cf. Eq. (7)). The numerical values agree very well with the
prediction Eq. (10) that is shown by a solid blue line. The parameter value of
the trajectory shown in panels (a) and (b) is indicated by a vertical red
line. For small networks the system occasionally tips into the direction
opposite to the perturbation. Hence, we plot the absolute value of $x^{\ast}$.

###  III.1 Moment expansion

In order to evaluate the time evolution of the cumulants of the $x_{k}$
distribution we note that taking averages $\left\langle\\_\right\rangle$ is a
linear operation such that it commutes with taking the time derivative,

$\displaystyle\dot{}\overline{x}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\left\langle
x_{k}\right\rangle=\left\langle\dot{x}_{k}\right\rangle$
$\displaystyle\dot{}\overline{v}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\left\langle
y_{k}^{2}\right\rangle=2\>\left\langle y_{k}\>\dot{x}_{k}\right\rangle.$

In order to evaluate the averages we express Eq. (2a) as

$\displaystyle\dot{x}_{k}$
$\displaystyle=-(\overline{x}+y_{k})^{3}+\overline{x}+y_{k}+{\sum_{l}}^{\prime}d_{kl}\bigl{(}\overline{x}+y_{l}\bigr{)}+\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}$
$\displaystyle=\overline{x}\,\left(1-\overline{x}^{2}+{\sum_{l}}^{\prime}d_{kl}\right)+y_{k}\>(1-3\,\overline{x}^{2})-3\,\overline{x}\,y_{k}^{2}$
$\displaystyle-
y_{k}^{3}+{\sum_{l}}^{\prime}d_{kl}\,y_{l}+\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}.$
(4)

Taking the average of Eq. (4) and observing that $\left\langle
y_{k}\right\rangle=0$ provides

$\displaystyle\dot{}\overline{x}$
$\displaystyle=\overline{x}\,\left(1-\overline{x}^{2}-3\,\overline{v}+\left\langle{\sum_{l}}^{\prime}d_{kl}\right\rangle\right)$
(5) $\displaystyle-\left\langle
y_{k}^{3}\right\rangle+\left\langle{\sum_{l}}^{\prime}d_{kl}\,y_{l}\right\rangle+\left\langle\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}\right\rangle.$

The expressions in the second line of this equation account for the third
cumulant of the distribution of the state variables $x_{k}$, a biased average
of the deviations $y_{k}$ from the average of the state variables, and the
impact of noise on the dynamics, respectively. In the present study we
consider a closure of the dynamics where the third cumulant and the biased
average will be suppressed. Moreover, we focus on the average behavior of an
ensemble of systems, where we perform two averages at the same time, (i) we
average over the coupling strength in the network realizations, (ii) and we
average over the noise $W_{k}$. Hence, the noise term vanishes, and
$\left\langle{\sum_{l}}^{\prime}d_{kl}\right\rangle$ takes the value
$p\,(N-1)\,d$, since each node is connected to $p\,(N-1)$ other nodes with an
average coupling strength $d$. Potential subtleties involved in taking the
different averages will be addressed in a forthcoming paper. Altogether we
thus find

$\displaystyle\dot{}\overline{x}$
$\displaystyle=\overline{x}\,\Bigl{(}1+p\,(N-1)\,d-\overline{x}^{2}-3\,\overline{v}\Bigr{)}.$
(6a)

The evolution of $v$ is obtained by multiplying Eq. (4) by $2y_{k}$ and taking
the average,

$\displaystyle\dot{}\overline{v}=2\left\langle
y_{k}\,\dot{x}_{k}\right\rangle=$ $\displaystyle
2\,\overline{v}\,(1-3\,\overline{x}^{2})+\left\langle
y_{k}\>\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}\right\rangle$
$\displaystyle-6\,\overline{x}\,\left\langle
y_{k}^{3}\right\rangle-2\,\left\langle
y_{k}^{4}\right\rangle+2\,\left\langle{\sum_{l}}^{\prime}d_{kl}\,y_{k}\,y_{l}\right\rangle.$

The term involving noise vanishes on account of considering a small-noise
limit of the dynamics. The terms in the second row of the equation are higher-
order moments and correlations that are dropped due to the assumptions of our
closure of the moment hierarchy. Hence, we obtain

$\displaystyle\dot{}\overline{v}=2\,\overline{v}\,(1-3\,\overline{x}^{2}).$
(6b)

###  III.2 Asymptotics of the order parameter

The dynamics, Eq. (6), involves a single dimensionless parameter,

$\displaystyle\lambda_{d}=1+p\,(N-1)\,d\,,$ (7)

that depends on the average number of connections, $p\,(N-1)$, of a tipping
element and on the average coupling strength, $d$. Hence, we denote it as
_connectivity._

Rescaled variables that are based on the connectivity will allow us to discuss
the asymptotics of the order parameter in a more transparent form

$\displaystyle x(\tau)$ $\displaystyle=\overline{x}/\sqrt{\lambda_{d}}$
$\displaystyle y(\tau)$ $\displaystyle=\overline{v}/\lambda_{d}$ (8)
$\displaystyle\tau$ $\displaystyle=\lambda_{d}\,t$

Denoting the derivatives of $x(\tau)$ and $y(\tau)$ with respect to $\tau$ as
$\dot{x}$ and $\dot{y}$ provides

$\displaystyle\dot{x}$
$\displaystyle=\frac{\mathrm{d}x}{\mathrm{d}\tau}=x\,\bigl{(}1-x^{2}-3\,v\bigr{)}$
(9a) $\displaystyle\dot{v}$
$\displaystyle=\frac{\mathrm{d}v}{\mathrm{d}\tau}=6\,v\,\bigl{(}\alpha-x^{2}\bigr{)}$
(9b)

where $\alpha=(3\lambda_{d})^{-1}$.

In Fig. 2(a) we show the phase-space flow of Eq. (9). The evolution has
nullclines with $\dot{x}=0$ at $x=0$ and $v=(1-x^{2})/2$ that are marked by
solid gray lines. The nullclines with $\dot{v}=0$ lie at $v=0$ and
$x=\sqrt{\alpha}$. They are marked by dashed lines. The three fixed points,
$p_{0}=(0,0)$, $\left(\sqrt{\alpha},(1-\alpha)/3\right)$, and $(1,0)$, are
located at the intersections of the $\dot{v}=0$ and $\dot{x}=0$ nullclines.

By inspection of the flow crossing the nullclines one readily verifies that
the fixed point $p_{i}$ has $i$ stable directions. For all $\alpha>0$ the
trajectories approach a state with $x=\pm 1$ and vanishing variance. A
representative sample trajectory is given by a solid red line in Fig. 2(a).
Panel (b) shows the time evolution of $x(\tau)$ (red line) and $v(\tau)$
(green line) of this trajectory.

The approach towards the stable fixed point at $x=1$ implies that

$\displaystyle
x^{\ast}=\lim_{\tau\to\infty}\overline{x}=\sqrt{}\lambda_{d}=\sqrt{1+d\,p\,(N-1)}$
(10)

Fig. 2(c) shows that this dependence is indeed observed by our numerical data.
Only for small connectivities, $\lambda_{d}$, the data tend to lie closely
below the prediction. We will come back to this point when we discuss Fig. 4.

##  IV Late-time evolution

The $\dot{x}=0$ nullcline always takes the same form (solid gray line in Fig.
2(a)), while the $\dot{v}=0$ nullcline is a vertical (dashed) line whose
position, $x_{c}=\sqrt{\alpha}=(3\,\lambda_{d})^{-1/2}$, depends on the system
parameters. For trajectories that proceed to the left of $x_{c}$ the order
parameter $x$ remains small while $v$ diverges. This behavior is unphysical
and it lies out of the scope of our model because our closure assumption can
only be expected to work for small $v$. However, for large $N$ the fixed point
$p_{1}$ lies very close to $(0,1/3)$, and we only consider initial conditions
where $0<v<x\ll 1$. For these initial conditions one encounters trajectories
as shown by a solid red line. Its $x$ coordinates grows monotonically from
zero towards one. Initially, the variance grows rapidly. Subsequently, (after
crossing the nullcline at $x=\sqrt{\alpha}$) it decays towards zero. For
larger networks ($N\gg 1$) and for more strongly coupled networks (increasing
$pd$) the crossover arises very close to zero
($\sqrt{\alpha}=[3\,(N-1)\,p\,d]^{-1/2}\ll 1$), and the maximum value of the
variance decreases.

###  IV.1 Approach to the stable fixed point

The approach towards $p_{2}$ is governed by the linearized Eq. (9), of the
deviations $\epsilon=x-1$ and $v$ from the fixed point. We take into account
that $\alpha\simeq 1/3$ for large connectivity, $\lambda_{d}$, and find

$\displaystyle\begin{pmatrix}\dot{\epsilon}\\\
\dot{v}\end{pmatrix}\simeq\begin{pmatrix}-2&-3\\\
0&-4\end{pmatrix}\begin{pmatrix}\epsilon\\\ v\end{pmatrix}\,.$ (11)

The eigenvalue for the decay in $v$ is therefore twice as large as the one in
$x$, such that $v\sim(1-x)^{2}$ close to the fixed point; in line with the
flow lines of the phase-space plot shown in Fig. 2(a).

Figure 3:  Collective Tipping Behavior of Different Network Types. The left
and right panels show the evolution of the order parameter and relative
deviation of numerical data from Eq. (2a) and the prediction from Eq. (13).
The rows show data for (a,b) Erdös-Rényi, (c,d) Barabási-Albert, and (e,f)
Watts-Strogatz networks with reconnection probability $\beta=0.8$. The other
parameters of the simulations are fixed to $N=100$, $p=0.1$, $d=1.5$, and
$\xi=0.1$, such that $\lambda_{d}\simeq 16$. The data points and the error
bars represent the expectation and the standard deviation over $50$
realizations of the network and the noise, respectively. The prediction, Eq.
(13), is indicated by a solid line in the left panels. It is shifted by
$\Delta\tau$ to match the data at $x=0.5$ (see main text for details).

In order to discuss the evolution of $x$ we introduce the variable $w=x^{2}$
and observe that

$\displaystyle\dot{w}$ $\displaystyle=2\,w\>\bigl{(}1-w-3\,v\bigr{)}$ (12a)
$\displaystyle\dot{v}$ $\displaystyle=2\,v\>\bigl{(}1-3\,w\bigr{)}$ (12b)

The evolution of $w(\tau)$ for identically vanishing variance, $v\equiv 0$, is
readily obtained by variable separation and partial fraction decomposition

$\displaystyle 2\,\tau$
$\displaystyle=\int_{w(0)}^{w(\tau)}\frac{\mathrm{d}w}{w\>(1-w)}=\int_{w(0)}^{w(\tau)}\mathrm{d}w\>\left[\frac{1}{w}+\frac{1}{1-w}\right]$
$\displaystyle=\ln\left[\frac{w(\tau)}{w(0)}\>\frac{1-w(0)}{1-w(\tau)}\right]$

such that

$\displaystyle
x(\tau)=\sqrt{w(\tau)}=\left[1+\frac{1-x^{2}(0)}{x^{2}(0)}\>\mathrm{e}^{-2\,\tau}\right]^{-1/2}\,.$
(13)

For late times, $\tau\gg 1$, we hence obtain the expected exponential decay

$\displaystyle
1-x(\tau)\simeq\frac{1-x^{2}(0)}{2\,x^{2}(0)}\;\mathrm{e}^{-2\,\tau}\qquad\text{
for \ }\tau\gg 1\,.$ (14)

In order to discuss the impact of finite values of $v$ we look for the
solutions of Eq. (12) for $w=1-\varepsilon$ where

$\displaystyle\dot{\varepsilon}$
$\displaystyle=-2\,\varepsilon\,\bigl{(}1-\varepsilon\bigr{)}-6\,v\,\bigl{(}1-\varepsilon\bigr{)}$
$\displaystyle\dot{v}$
$\displaystyle=-4\,v\>\left(1+\frac{3}{2}\,\varepsilon\right)$

We approximate the equation for $\dot{\varepsilon}$ by inserting
$v=c\,\varepsilon^{2}$ (cf. the discussion of the relation of $v$ and
$\epsilon$ below Eq. (11)), and consider terms only till order
$\varepsilon^{3}$ and $v\,\varepsilon$. As a consequence

$\displaystyle(1+3c)\,\dot{\varepsilon}$
$\displaystyle=-2\,(1+3c)\,\varepsilon\,\bigl{(}1-(1+3c)\,\varepsilon\bigr{)}$

such that $\varepsilon(\tau)=1-x(\tau)$ remains unchanged up to multiplication
by the factor $1+3c$. However, in view of Eq. (14) such a factor can be
accounted for by a shift, $\Delta\tau$, of time, yielding the following
prediction for the time evolution of the order parameter,

$\displaystyle
x_{\text{theory}}(\tau)\simeq\left[1+\frac{1-x^{2}(0)}{x^{2}(0)}\>\mathrm{e}^{-2\,(\tau-\Delta\tau)}\right]^{-1/2}\,.$
(15)

This expression has a single free parameter $\Delta\tau$ that amounts to a
time delay to reach the half time $\tau_{1/2}$, where
$x_{\text{theory}}(\tau_{1/2})=1/2$,

$\displaystyle\tau_{1/2}=\Delta\tau-\frac{1}{2}\>\ln\frac{3\,x^{2}(0)}{1-x^{2}(0)}\,.$
(16)

###  IV.2 Numerical test of Eq. (15)

The left columns of Fig. 3 show the evolution of the order parameter for ER,
BA, and WS networks. As initial condition we consider a system where one node,
let it be $\ell$, is displaced to a value $x_{\ell}=\sqrt{\lambda_{d}}$, and
all other nodes, $m\neq\ell$, start at $x_{m}=0$. This amounts to initial
conditions with

$\displaystyle x(0)$ $\displaystyle=\frac{\left\langle
x_{k}\right\rangle_{\tau=0}}{\sqrt{}\lambda_{d}}=\frac{1}{\sqrt{}\lambda_{d}}\>\frac{1}{N}\sum_{k=1}^{N}x_{k}(0)$
$\displaystyle=\frac{0\times(N-1)+1\times\sqrt{\lambda_{d}}}{N\,\sqrt{}\lambda_{d}}=\frac{1}{N}\,.$

For each type of network we perform $50$ runs with different realizations of
the networks and noise. The expectation of the order parameter,
$x_{\text{sim}}$, for the resulting evolution is given by a circle, and the
error bars mark the standard deviation over the simulation runs.

The theoretical prediction, $x_{\text{theory}}$, is provided by a solid line.
The time offsets are determined such that the numerical data and theoretical
prediction match for $x=1/2$. In a logarithmic plot of the relative deviation
$(x_{\text{sim}}/x_{\text{theory}})/x_{\text{theory}}$ (Fig. 3(b,d,f)) this
results in a singularity at time $\tau_{1/2}$ that is marked by vertical
dotted lines. The lines align with the singularity for

$\displaystyle\Delta\tau_{\text{ER}}$ $\displaystyle=0.577\,,$
$\displaystyle\Delta\tau_{\text{BA}}$ $\displaystyle=-0.96\,,$ (17)
$\displaystyle\Delta\tau_{\text{WS}}$ $\displaystyle=0.527\,.$

ER and WS networks have very similar time offsets of about $\Delta\tau\simeq
0.5$. In contrast BA networks decay significantly faster. The negative value
of $\Delta\tau_{\text{BA}}$ implies that these networks decay even faster than
predicted by the $v=0$ result, Eq. (13). Follow-up work will have to clarify
if and how this intriguing finding is related to the existence of hubs in BA
networks, nodes with a degree greatly exceeding the average.

In addition to the position of $\tau_{1/2}$, Fig. 3(b,d,f) demonstrates that
the prediction, Eq. (15), provides an excellent description of the data for
times beyond $\tau_{1/2}$. The time offset $\Delta\tau$ is the only fitting
parameter in this description.

###  IV.3 Influence of Initial Displacement

Figure 4:  Trajectories for Different Initial Displacements $x(0)={\nu}/{N}$.
Calculations using Eq. (2a) have been conducted for ER networks with $d=0.1$,
$p=0.25$, $\xi=0.1$, $N\in\\{100,200,300,800\\}$, and $\nu\in\\{1,2,3\\}$,
with couplings $d_{kl}$ drawn from a uniform distribution $d_{kl}\in(0,0.2)$.
This amounts to connectivities of $\lambda_{d}=3.475$, $5.975$, $10.975$, and
$20.975$, respectively. The data and their error bars in the upper panels
represent the expectation and standard deviation of $100$ simulations runs
with independent realizations of the network and the noise. The corresponding
lower panels shows the relative deviation of the trajectory from the
prediction Eq. (15) with $x_{0}=\nu/N$ and a time shift $\Delta\tau$ provided
in the respective legends.

In the previous sections one node has been perturbed by a fixed value
initially. To discuss the effect of larger initial perturbations we explore
now systems where initially $\nu$ nodes are tipped to the asymptotic
equilibrium value $x^{\ast}=\sqrt{}\lambda_{d}$. Consequently,

$\displaystyle x(0)$ $\displaystyle=\frac{\left\langle
x_{k}\right\rangle_{t=0}}{\sqrt{}\lambda_{d}}=\frac{1}{\sqrt{}\lambda_{d}}\>\frac{1}{N}\sum_{k=1}^{N}x_{k}(0)$
$\displaystyle=\frac{0\times(N-\nu)+\nu\times\sqrt{\lambda_{d}}}{N\,\sqrt{}\lambda_{d}}=\frac{\nu}{N}\,.$
(18)

According to Eq. (15) this change in $x(0)$ will be the only change in
$x_{\text{theory}}$. This prediction is scrutinized in Figure 4 where we show
data for ER networks with $100$, $200$, $400$, and $800$ nodes where we
excited initially $\nu=1$, $2$ or $3$ nodes:

  1. 1.

For a given value of $N$ the time offsets for different $\nu$ are indeed
constant to within our numerical accuracy. Moreover, the positions of the
dotted lines clearly show the decrease of $\tau_{1/2}$ that is predicted by
Eq. (16).

  2. 2.

For increasing $N$ the time offset decreases. This is in line with the
expectation that the variance reaches smaller values for increasing
$\lambda_{d}$ such that the dynamics more closely follows the $v=0$ dynamics.

  3. 3.

The relative deviation does not decay to zero. Indeed on a second sight one
recognizes that the data in Fig. 3 also tend to saturate at a $10^{-5}$ level.
Hence, the data do not quite saturate at $\sqrt{}\lambda_{d}$ but at a
slightly smaller value, as also observed for the small connectivity data in
Fig. 2. We attribute this effect to the small noise that enforces $v$ to take
a small positive value. According to Eq. (9a) this reduces the asymptotic
value of $x$ by a factor of $\sqrt{1-3\,v^{2}}$. Consequently, the relative
error decreases like $1-\sqrt{1-3v}\simeq 1.5\overline{v}/\lambda_{d}$. The
trend is clearly visible in the data.

Altogether, we established an analytic description for the time evolution of
the order parameter of coupled networks subject to small noise levels. For
different paradigmatic network types, the collective dynamics of the network
behaves just as a single aggregated tipping element with a new time scale and
equilibrium values selected by the connectivity $\lambda_{d}$, as provided by
the scaling in Eq. (8).

##  V  Dynamics and stability beyond the normal form

Up to now we considered systems that evolve according to Eq. (2a). In this
section we explore how the dynamics changes when there are additional drift
and offset terms, and how rescaling and shifting the variables $x_{k}$ affects
the parameters of the resulting description.

###  V.1 Influence of drift and offset terms

To discuss the impact of drift and offset terms we write Eq. (2b) in the form
of Eq. (4) with two additional terms,

$\displaystyle\dot{x}_{k}=$
$\displaystyle\overline{x}\,\left(1-\overline{x}^{2}+{\sum_{l}}^{\prime}d_{kl}\right)+y_{k}\>(1-3\,\overline{x}^{2})-3\,\overline{x}\,y_{k}^{2}$
$\displaystyle-
y_{k}^{3}+{\sum_{l}}^{\prime}d_{kl}\,y_{l}+\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}$
$\displaystyle+d_{k}-{\sum_{l}}^{\prime}d_{kl}\,r_{kl}$

The terms in the first row contributed to the average of $\dot{x}_{k}$, as
evaluated in Eq. (5). The terms in the second row amount to higher-order
cumulants, correlations, and noise that are dropped for the closure and the
weak-noise limit adopted here. The terms in the third row account for the
drift and offset terms. Let the expectation of the drift terms $d_{k}$ be
$r_{D}$, and for the offsets $r_{kl}$ in the coupling terms we introduce

$\displaystyle\left\langle{\sum_{l}}^{\prime}d_{kl}\,r_{kl}\right\rangle=p\,(N-1)\,d\,r_{O}=(\lambda_{d}-1)\,r_{O}\,.$

(Again, there may be correction terms that are subdominant for the considered
closure.) With these notations the average of $\dot{x}_{k}$ takes the form

$\displaystyle\dot{}\overline{x}$
$\displaystyle=\overline{x}\,\Bigl{(}\lambda_{d}-\overline{x}^{2}-3\,\overline{v}\Bigr{)}+r_{D}-(\lambda_{d}-1)\,r_{O}$
(19)

and $\dot{}\overline{v}$ still takes the form of Eq. (6b). After all, the
additional terms are constant such that they average to zero when multiplied
by $y_{k}$. In terms of the rescaled variables, Eq. (8) we hence find

$\displaystyle\dot{x}$
$\displaystyle=x\,\bigl{(}1-x^{2}-3\,v\bigr{)}+\frac{r_{D}+r_{O}}{\lambda_{d}^{3/2}}-\frac{r_{O}}{\sqrt{}\lambda_{d}}\,,$
(20a) $\displaystyle\dot{v}$
$\displaystyle=6\,v\,\bigl{(}\alpha-x^{2}\bigr{)}\,,$ (20b)

We note that for large connectivity $\lambda_{d}=1+p\,d\,(N-1)$ the drift and
offset terms are sub-dominant in the order-parameter dynamics, and they
disappear in the large-network limit. Hence, the order parameter of strongly
coupled, large networks evolves according to Eq. (9), as discussed in Sect.
III.

For large connectivity, $\lambda_{d}\gg 1$, every system that crosses a
tipping point relaxes to its new equilibrium value according to these
equations. The intrinsic order-parameter dynamics is symmetric. In terms of
the effective coupling

$\displaystyle
r_{\text{eff}}=\frac{r_{D}-r_{O}}{\lambda_{d}^{3/2}}-\frac{r_{O}}{\sqrt{}\lambda_{d}}$
(21)

one encounters saddle-node bifurcations at $r_{\text{eff}}=\pm\sqrt{4/27}$
(cf. Fig. 1 and the discussion below Eq. (1)).

###  V.2  Linear coordinate transformations

Scale changes can be expressed as a change of variables where the state of the
nodes is characterized by new variables $z_{k}$ that are linear functions of
$x_{k}$,

$\displaystyle z_{k}=m\,x_{k}+n\qquad\text{for fixed \ }m,\>n\,.$ (22)

Consequently, $\overline{z}=\left\langle z_{k}\right\rangle=m\,\overline{x}+n$
and $\dot{z}_{k}=m\,\dot{x}_{k}$. Averaging the latter equation and using Eq.
(19) provides

$\displaystyle\frac{\dot{}\overline{z}}{m}=\dot{}\overline{x}=\overline{x}\>\left(\lambda_{d}-\overline{x}^{2}-3\,\overline{v}\right)+r_{D}-(\lambda_{d}-1)\>r_{O}$

For $z=\overline{z}/(m\sqrt{}\lambda_{d})$, $\tau=\lambda_{d}\,t$, and
$z_{0}=n/(m\sqrt{}\lambda_{d})$ we have $x=z-z_{0}$, and we find

$\displaystyle\dot{z}(\tau)$
$\displaystyle=\frac{\mathrm{d}z}{\mathrm{d}\tau}=(z-z_{0})\>\left[1-(z-z_{0})^{2}-3\,v\right]+r_{\text{eff}}$
(23a) $\displaystyle\dot{v}(\tau)$
$\displaystyle=6\,v\>\left(\alpha-(z-z_{0})^{2}\right)$ (23b)

with a macroscopic bifurcation parameter, $r_{\text{eff}}$, given by Eq. (21).

The shift $z_{0}$ of $z$ can be interpreted as a vertical displacement of the
s-shaped red curve in Fig. 1. It affects the positions of the steady states,
but not the bifurcation parameters. The intrinsic dynamics remains symmetric
in the bifurcation parameter $r_{\text{eff}}$ with saddle-node bifurcations at
$r_{\text{crit}}=\pm\sqrt{4/27}$.

The dependence of the bifurcation parameter, Eq. (21), on connectivity
$\lambda_{d}$ provides another important insight into the stability of the
networks: The drift dominates the control parameter as long as the offset
terms $r_{kl}$ are very small as compared to the drift,
$(\lambda_{d}-1)\>|r_{O}|\ll|r_{D}|$, and the global system is bistable for
$-\lambda_{d}^{3/2}\lesssim 2.5\,r_{D}\lesssim\lambda_{d}^{3/2}$ (after all
$\sqrt{27/4}\approx 5/2$). In view of Eq. (7) the connectivity $\lambda_{d}$
increases linearly with the network size, such that the state of the network
in this parameter range is very robust against perturbations. However, for
strong connectivity the offset terms will always dominate the system response
eventually,

$\displaystyle(\lambda_{d}-1)\>|r_{O}|\gg|r_{D}|\quad\text{for \
}\lambda_{d}\ggg 1\,.$

Eventually, the coupling terms govern the strength of the bifurcation
parameter, and the range of stability increases more slowly with the
connectivity

$\displaystyle-\sqrt{{\textstyle\frac{4}{27}}\>\lambda_{d}}\;\lesssim\>r_{O}\>\lesssim\;\sqrt{{\textstyle\frac{4}{27}}\>\lambda_{d}}\quad\text{for
\ }\lambda_{d}\ggg 1\,.$

Still, the range of stability increases with network size and for increasing
network size the state will become ever more robust against perturbations.

Figure 5:  Asymptotic fixed point $z^{\ast}$ (a) and the critical points (b)
for networks with an offset $n$ and $r_{kl}=0$ in Eq. (24). All calculations
are conducted for an ER network with $N=100$ nodes, couplings $d_{kl}$ drawn
from a uniform distribution $d_{kl}\in(0,1.2)$, $\lambda_{d}=16$, and noise
strength $\xi=0.1$. (a) For fixed values $n=m=1$ and different drift $r_{D}$
we follow the dynamics for initial conditions $x_{\text{init}}=-1$ and one
excited node at $z=1$ (orange) and $x_{\text{init}}=2$ and one node at $z=-1$
(blue), respectively. The dots and the error bars indicate the average and the
standard deviation of the late-time asymptotic positions of fifty simulations
with different realizations of the network and noise. They lie nicely on the
prediction for the positions of the fixed points, Eq. (25a), that is provided
by the solid red line. The positions of the critical points, $r_{\text{crit}}$
are identified as having the largest distance between subsequent data points.
(b) The $n$-dependence of $r_{\text{crit}}$ determined from data as shown in
panel (a), $m=1$, and different values of $n$. The red lines provide the
predicted positions, Eq. (25b), of the saddle-node bifurcations.

###  V.3  Coordinate transformation that do not involve the coupling

It is by no means self-evident that a change of coordinates of individual
tipping elements affects the coupling terms in the same way. In particular,
studies of climate networks [46] address the situation where the individual
tipping elements are shifted at a _fixed_ symmetric coupling. In other words
one adopts Eq. (22) but considers the dynamics

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\frac{z_{k}}{m}=$
$\displaystyle-\left(\frac{z_{k}-n}{m}\right)^{3}+\frac{z_{k}-n}{m}+d_{k}$
$\displaystyle+{\sum_{l}}^{\prime}d_{kl}\bigl{(}z_{l}-r_{kl}\bigr{)}+\xi\frac{\mathrm{d}W_{k}}{\mathrm{d}t}$
(24)

with $\left\langle r_{kl}\right\rangle=0$. To identify the effects of this
asymmetry in the parameter change we observe that

$\displaystyle{\sum_{l}}^{\prime}d_{kl}\bigl{(}z_{l}-r_{kl}\bigr{)}={\sum_{l}}^{\prime}(m\,d_{kl})\>\left(\frac{z_{l}-n}{m}-\frac{r_{kl}-n}{m}\right)$

such that the substitution $d_{kl}\to m\,d_{kl}$ and $r_{kl}\to(r_{kl}-n)/m$
provides the system treated in Sect. V.2. On the level of the coarse-grained
system this entails that for the present system $\lambda_{d}=1+(N-1)\,p\,md$
and
$r_{O}=\left\langle(r_{kl}-n)/m\right\rangle=-n/m=-\sqrt{\lambda_{d}}\,z_{0}$.
Substitution into Eq. (23a) provides

$\displaystyle\dot{z}$
$\displaystyle=(z-z_{0})\>\left[1-(z-z_{0})^{2}-3\,v\right]+\frac{r_{D}+r_{O}}{\lambda_{d}^{3/2}}-\frac{r_{O}}{\sqrt{}\lambda_{d}}$
$\displaystyle=(z-z_{0})\>\left[1-(z-z_{0})^{2}-3\,v\right]+\frac{\lambda_{d}-1}{\lambda_{d}}\>z_{0}+r_{c}$

with control parameter

$\displaystyle r_{c}={r_{D}}/{\lambda_{d}^{3/2}}\,.$

Now, the fixed points, $\dot{z}=0$, are located on the line

$\displaystyle
r_{c}=\frac{r_{D}}{\lambda_{d}^{3/2}}=-(z^{\ast}-z_{0})\>\left[1-(z^{\ast}-z_{0})^{2}\right]-\frac{\lambda_{d}-1}{\lambda_{d}}\>z_{0}\,,$
(25a) and the saddle-node bifurcations arise when $\displaystyle
r_{\text{crit}}=\pm\sqrt{\frac{4}{27}}-\frac{\lambda_{d}-1}{\lambda_{d}}\>z_{0}\,.$
(25b)

Figure 5 demonstrates that this is indeed observed in simulations of the
networks. The right-hand side of Eq. (25a) mounts to the red s-shaped curve in
Fig. 5(a). The data points mark the positions of the fixed points $z^{\ast}$
adopted by the system. For each value of $z_{0}$ we thus find a hysteresis
loop that provides the critical driving parameters, $r_{\text{crit}}$, where
the saddle-node bifurcations arise. Figure 5(b) displays the $n$-dependence of
the positions of the saddle-node bifurcations. The data nicely follow the
expected behavior, Eq. (25b), with $z_{0}=n/(m\sqrt{}\lambda_{d})$. We
conclude that incongruent response of the tipping elements and their coupling
to parameter changes, as expressed by Eq. (24) and (22), gives rise to
dynamics where the system is bistable in a range of control parameters $r_{c}$
that is no longer symmetric around zero.

##  VI Discussion and Outlook

In the main part of the present paper we set up an analytical description of
the collective tipping of a large number, $N$, of tipping elements that are
coupled in a random network. The present section discusses these findings in
the context of related literature, and it highlights prospective extensions.

###  VI.1 Collective behavior

We characterized each tipping element $k\in\\{1\cdots N\\}$ by a scalar state
variable $x_{k}\in\mathbb{R}$. It evolves in a double-well potential with a
saddle at $x=0$ that separates the domains of attraction of two coexisting
stable states at $x_{k}=\pm 1$. The elements are pairwise connected in a
random network with mean coupling strength $d$, and a probability $p$ that
there is a connection for a given pair of elements.

The overall state of the system is characterized by the expectation
$\overline{x}$ and the variance $\overline{v}$ of the distribution of $x_{k}$.
The evolution of these collective coordinates is derived based on a closure
where we drop higher-order cumulants and correlations. In the weak-noise limit
they evolve according to Eq. (6) that involves a single dimensionless
parameter, the connectivity $\lambda_{d}=1+p\,(N-1)\,d$. The coarse-grained
dynamics has global fixed points at
$(\overline{x},\overline{v})=(\pm\sqrt{}\lambda_{d},0)$, and tipping proceeds
in terms of a dimensionless time $\lambda_{d}\,t$ (Fig. 2).

The dependence on $\lambda_{d}\,t$ implies that tipping occurs more abruptly
for increasing network degree $p\,(N-1)$ and coupling strength $d$. This
provides a rational for the observation of Eom [47] who also observed these
trends in numerical data.

The parameter dependence of the time scale and the position of the fixed
points of the collective coordinates can be absorbed by rescaling, Eq. (8). In
line with the observations of Wunderling  _et al._ [24] this implies that
systems with high average degree $p\,(N-1)$ only require low coupling
strengths $d$ to initiate a tipping cascade.

We have analyzed network-specific time delays in Eq. (17) and found very
similar offsets for ER and WS networks of about $\Delta\tau\simeq 0.5$.
Interestingly, BA networks decay with a significantly smaller time offset. We
expect this to result from the existence of hubs in the BA network that could
force the tipping process of the network. Follow-up work will have to clarify
if and how this finding is related to the existence of hubs in the BA network.

###  VI.2 Shape of the distribution

In terms of the rescaled collective coordinates
$x=\overline{x}/\sqrt{}\lambda_{d}$ and $v=\overline{v}/\lambda_{d}$ the
dynamics has a saddle point at position
$\Bigl{(}\alpha,(1-\alpha^{2})/3\Bigr{)}$ with $\alpha=(3\lambda_{d})^{-1/2}$.
Trajectories that proceed to the right of the saddle show collective behavior.
Up to initial transients the variance of the distribution of trajectory states
is small and the systems evolves according to a strongly-coupled collective
dynamics. Due to the assumptions taken to close the equations for $x$ and $v$
the present theory does not apply for trajectories that proceed to the left of
the saddle. Follow-up work will address higher-order closures in order to deal
with the broad and, at times, bimodal distributions of tipping elements in
this polydisperse parameter regime.

At this point we observe that $\alpha$ approaches zero for increasing
connectivity. Hence, a rapidly increasing fraction of initial conditions
proceed to the right of the saddle. In line with the findings of Brummitt  _et
al._ [16] this explains that increasing coupling strength $d$ entails
synchronous behavior of tipping elements for a larger set of initial
conditions.

###  VI.3 Noise

Careful inspection of our data in Fig. 4 revealed that the expectation
$\overline{x}$ saturates at values slightly below $\sqrt{}\lambda_{d}$. We
argued that the variance $\overline{v}$ of the distribution remains finite for
a system subjected to small noise, and that this will decrease the asymptotic
value by a factor $1-3\,\overline{v}/\lambda_{d}$. Moreover, preliminary data
indicate that noise has a non-trivial impact on the dynamics in the
polydisperse parameter regime. An expansion of our model that addresses the
impact of noise is in preparation.

###  VI.4 Negative feedback and multiple scales

The derivation of the collective dynamics, Eqs. (6), does not involve any
assumptions on the distribution of the coupling strength. The present study
describes the time evolution of systems where the feedback is positive on
average, $d>0$. In this case $\lambda_{d}$ increases for larger and more
strongly coupled networks, and the dynamics does not much depend on the
structure and realization of the network (cf. Fig. 3).

For negative overall feedback, $d<0$, and larger networks our connectivity as
defined in Eq. (7) becomes negative and thus no longer meaningful. The
nondimensionalization introduced in Eq. (8) to dimensionless parameters only
works for $\lambda_{d}>0$. In the case of negative overall feedback, the flow
described by Eqs. (6) will always enter the realm of large $\overline{v}$ that
is out of the scope of the present assumptions for the closure of the
equations. An appropriate expansion of the model will allow us to address
phenomena such as Kadyrov-oscillations [48, 14, 49].

We also expect that systems with noticeable negative feedback and vast
heterogeneity in the time scales of the tipping elements and their tipping
thresholds are much more susceptible to correlations of the characteristics of
tipping elements and their network environment. Expansions that account for
these correlations will make it possible to investigate systems like the
large-scale, Amazon rainforest, where spatially distant patches are dependent
on each other via the atmospheric moisture recycling feedback [e.g. 24, 50,
22]. Instead of a micro-scale (local-scale tipping elements) and a macro-scale
(entire forest), one will then also define an intermediary meso-scale for
strongly connected subcomponents of the forest.

###  VI.5 Time dependence of parameters

The present analysis provides a comprehensive framework for the description of
phase separation of networks of tipping elements for fixed-in-time control
parameters. However, the Earth climate system is subjected to a sustained
change of its parameters due to the release of vast amounts of carbon dioxide
into the atmosphere, and on historical scales its energy input varies
periodically due to Milankovitch cycles [51, 52]. Both effects can have a
severe impact on the time evolution. Parameter oscillations in a noisy
bistable system induce stochastic resonance [53, 54]. Parameter drift induces
size focusing [55] as discussed recently for applications in chemistry [56]
and soft biological matter [57]. Non-trivial new behavior emerges in the
latter systems due to overall constraints on the dynamics due to mass
conservation, or global pressure and elastic fields. It will therefore be
highly relevant to explore how the dynamics is impacted by the coupling to a
global parameter, like temperature in the Daisyworld models [31].

##  VII Conclusion

In the present paper we established a description of the collective tipping
dynamics of an assembly of coupled tipping elements that are set up in a
slightly perturbed unstable steady state. We described the assembly in terms
of the expectation and the variance of variables that describe the state of
the individual tipping elements. For a vast range of initial conditions of the
network the variance remains small throughout the evolution. In those cases
the expectation follows universal dynamics, as shown in Figs. 3 and 4 for
networks with vastly different numbers of nodes, $N$, network degrees,
$p\,(N-1)$, network types (Erdös-Rényi [40], Barabási-Albert [41], and Watts-
Strogatz [42]), and average coupling, $d$, between the elements. The dynamics
is characterized by a single dimensionless parameter, $\lambda_{d}$, that is
defined in Eq. (7). Based on the nondimensionalization, Eq. (8), it provides
the time scale of the dynamics and the parameter values of the expectation and
the variance of the target state of the dynamics (Fig. 2). In Sect. V we
showed that the description is robust in the large network and strong coupling
limit (see Fig. 5), as long as the coupling provides a positive feedback on
average. Expected changes for negative feedback, $d<0$, and other extensions
of the model were briefly discussed in Sect. VI.

In conclusion, in the present work we established an effective analysis of
cascading tipping behavior in strongly coupled networks. It provides a
comprehensive analytical description that is in excellent agreement with
numerical simulations, and it calls for extensions to address the tipping
dynamics in other parameter regimes and in settings with a global feedback or
other constraints.

###### Acknowledgements.

N.W. acknowledges support from the the IRTG 1740/TRP 2015/50122-0 funded by
DFG and FAPESP. N.W. is grateful for a scholarship from the Studienstiftung
des deutschen Volkes. J.F.D. is thankful for support by the Leibniz
Association (project DominoES) and the European Research Council project Earth
Resilience in the Anthropocene (743080 ERA). The authors gratefully
acknowledge the European Regional Development Fund (ERDF), the German Federal
Ministry of Education and Research and the Land Brandenburg for supporting
this project by providing resources on the high performance computer system at
the Potsdam Institute for Climate Impact Research.

## References

  * Safari  _et al._ [2017] A. Safari, R. Saidur, F. A. Sulaiman, Y. Xu, and J. Dong, A review on supercooling of phase change materials in thermal energy storage systems, Renewable and Sustainable Energy Reviews 70, 905 (2017).
  * Wunderlich [2007] B. Wunderlich, One hundred years research on supercooling and superheating, Thermochimica Acta 461, 4 (2007).
  * Binder [1987] K. Binder, Theory of first-order phase transitions, Reports on Progress in Physics 50, 783 (1987).
  * Kuwahara  _et al._ [1995] H. Kuwahara, Y. Tomioka, A. Asamitsu, Y. Moritomo, and Y. Tokura, A first-order phase transition induced by a magnetic field, Science 270, 961 (1995).
  * Onuki [2002] A. Onuki,  _Phase Transition Dynamics_ (Cambridge University Press, Cambridge, 2002).
  * da Costa  _et al._ [2010] R. A. da Costa, S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Explosive percolation transition is actually continuous, Physical Review Letters 105, 255701 (2010).
  * D’Souza and Mitzenmacher [2010] R. M. D’Souza and M. Mitzenmacher, Local cluster aggregation models of explosive percolation, Physical Review Letters 104, 195702 (2010).
  * D’Souza and Nagler [2015] R. M. D’Souza and J. Nagler, Anomalous critical and supercritical phenomena in explosive percolation, Nature Physics 11, 531 (2015).
  * Grassberger  _et al._ [2011] P. Grassberger, C. Christensen, G. Bizhani, S.-W. Son, and M. Paczuski, Explosive percolation is continuous, but with unusual finite size behavior, Physical Review Letters 106, 225701 (2011).
  * Riordan and Warnke [2011] O. Riordan and L. Warnke, Explosive percolation is continuous, Science 333, 322 (2011).
  * Buldyrev  _et al._ [2010] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin, Catastrophic cascade of failures in interdependent networks, Nature 464, 1025 (2010).
  * Motter and Timme [2018] A. E. Motter and M. Timme, Antagonistic phenomena in network dynamics, Annual Review of Condensed Matter Physics 9, 463 (2018).
  * Watts [2002] D. J. Watts, A simple model of global cascades on random networks, Proceedings of the National Academy of Sciences 99, 5766 (2002).
  * Abraham  _et al._ [1991] R. Abraham, A. Keith, M. Koebbe, and G. Mayer-Kress, Computational unfolding of double-cusp models of opinion formation, International Journal of Bifurcation and Chaos 1, 417 (1991).
  * Lenton  _et al._ [2019] T. M. Lenton, J. Rockström, O. Gaffney, S. Rahmstorf, K. Richardson, W. Steffen, and H. J. Schellnhuber, Climate tipping points—too risky to bet against, Nature 575, 592 (2019).
  * Brummitt  _et al._ [2015] C. D. Brummitt, G. Barnett, and R. M. D’Souza, Coupled catastrophes: sudden shifts cascade and hop among interdependent systems, Journal of the Royal Society Interface 12, 20150712 (2015).
  * Feudel  _et al._ [2018] U. Feudel, A. N. Pisarchik, and K. Showalter, Multistability and tipping: From mathematicsand physics to climate and brain—minireview and preface to the focus issue, Chaos 28, 033501 (2018).
  * Scheffer [1989] M. Scheffer, Alternative stable states in eutrophic, shallow freshwater systems: a minimal model, Hydrobiological Bulletin 23, 73 (1989).
  * Hughes  _et al._ [2013] T. P. Hughes, S. Carpenter, J. Rockström, M. Scheffer, and B. Walker, Multiscale regime shifts and planetary boundaries, Trends in Ecology & Evolution 28, 389 (2013).
  * Lenton and Williams [2013] T. M. Lenton and H. T. P. Williams, On the origin of planetary-scale tipping points, Trends in Ecology & Evolution 28, 380 (2013).
  * Rocha  _et al._ [2018] J. C. Rocha, G. Peterson, Ö. Bodin, and S. Levin, Cascading regime shifts within and across scales, Science 362, 1379 (2018).
  * Zemp  _et al._ [2017] D. C. Zemp, C.-F. Schleussner, H. M. J. Barbosa, M. Hirota, V. Montade, G. Sampaio, A. Staal, L. Wang-Erlandsson, and A. Rammig, Self-amplified Amazon forest loss due to vegetation-atmosphere feedbacks, Nature Communications 8, 1 (2017).
  * Steffen  _et al._ [2018] W. Steffen, J. Rockström, K. Richardson, T. M. Lenton, C. Folke, D. Liverman, C. P. Summerhayes, A. D. Barnosky, S. E. Cornell, M. Crucifix, _et al._ , Trajectories of the Earth system in the Anthropocene, Proceedings of the National Academy of Sciences 115, 8252 (2018).
  * Wunderling  _et al._ [2020a] N. Wunderling, B. Stumpf, J. Krönke, A. Staal, O. A. Tuinenburg, R. Winkelmann, and J. F. Donges, How motifs condition critical thresholds for tipping cascades in complex networks: Linking micro-to macro-scales, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 043129 (2020a).
  * Krönke  _et al._ [2020] J. Krönke, N. Wunderling, R. Winkelmann, A. Staal, B. Stumpf, O. A. Tuinenburg, and J. F. Donges, Dynamics of tipping cascades on complex networks, Physical Review E 101, 042311 (2020).
  * Garbe  _et al._ [2020] J. Garbe, T. Albrecht, A. Levermann, J. F. Donges, and R. Winkelmann, The hysteresis of the Antarctic ice sheet, Nature 585, 538 (2020).
  * Robinson  _et al._ [2012] A. Robinson, R. Calov, and A. Ganopolski, Multistability and critical thresholds of the greenland ice sheet, Nature Climate Change 2, 429 (2012).
  * Curry  _et al._ [1995] J. A. Curry, J. L. Schramm, and E. E. Ebert, Sea ice-albedo climate feedback mechanism, Journal of Climate 8, 240 (1995).
  * Müller-Stoffels and Wackerbauer [2011] M. Müller-Stoffels and R. Wackerbauer, Regular network model for the sea ice-albedo feedback in the arctic, Chaos 21, 013111 (2011).
  * Kriegler  _et al._ [2009] E. Kriegler, J. W. Hall, H. Held, R. Dawson, and H. J. Schellnhuber, Imprecise probability assessment of tipping points in the climate system, Proceedings of the National Academy of Sciences 106, 5041 (2009).
  * Watson and Lovelock [1983] A. J. Watson and J. E. Lovelock, Biological homeostasis of the global environment: the parable of Daisyworld, Tellus B: Chemical and Physical Meteorology 35, 284 (1983).
  * Budyko [1969] M. I. Budyko, The effect of solar radiation variations on the climate of the Earth, Tellus 21, 611 (1969).
  * Sellers [1969] W. D. Sellers, A global climatic model based on the energy balance of the Earth-atmosphere system, Journal of Applied Meteorology and Climatology 8, 392 (1969).
  * Hoffman  _et al._ [1998] P. F. Hoffman, A. J. Kaufman, G. P. Halverson, and D. P. Schrag, A Neoproterozoic Snowball Earth, Science 281, 1342 (1998).
  * Lucarini and Bódai [2019] V. Lucarini and T. Bódai, Transitions across melancholia states in a climate model: Reconciling the deterministic and stochastic points of view, Physical Review Letters 122, 158701 (2019).
  * Callen [1985] H. B. Callen,  _Thermodynamics and an Introduction to Thermostatistics_ (Wiley, New York, 1985).
  * Wunderling  _et al._ [2021a] N. Wunderling, J. Krönke, V. Wohlfarth, J. Kohler, J. Heitzig, A. Staal, S. Willner, R. Winkelmann, and J. F. Donges, Modelling nonlinear dynamics of interacting tipping elements on complex networks: the pycascades package, The European Physical Journal Special Topics , 1 (2021a).
  * Klose  _et al._ [2020] A. K. Klose, V. Karle, R. Winkelmann, and J. F. Donges, Emergence of cascading dynamics in interacting tipping elements of ecology and climate, Royal Society Open Science 7, 200599 (2020).
  * Murdock [2006] J. Murdock, Normal forms, Scholarpedia 1, 1902 (2006).
  * Erdös and Rényi [1959] P. Erdös and A. Rényi, On random graphs, Publicationes Mathematicae 6, 290 (1959).
  * Barabási and Albert [1999] A.-L. Barabási and R. Albert, Emergence of scaling in random networks, Science 286, 509 (1999).
  * Watts and Strogatz [1998] D. J. Watts and S. H. Strogatz, Collective dynamics of ‘small-world’ networks, Nature 393, 440 (1998).
  * Hagberg  _et al._ [2008] A. Hagberg, P. Swart, and D. S Chult,  _Exploring network structure, dynamics, and function using NetworkX_ , Tech. Rep. (Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008).
  * Aburn [2017] M. J. Aburn, sdeint: Numerical integration of stochastic differential equations (SDE), software available under GNU general public license (2017), accessed: 2020-07-18.
  * Arnold [1992] V. I. Arnold,  _Catastrophe theory_ , 3rd ed. (Springer, Berlin, 1992).
  * Wunderling  _et al._ [2021b] N. Wunderling, J. F. Donges, J. Kurths, and R. Winkelmann, Interacting tipping elements increase risk of climate domino effects under global warming, Earth System Dynamics 12, 601 (2021b).
  * Eom [2018] Y.-H. Eom, Resilience of networks to environmental stress: From regular to random networks, Physical Review E 97, 042313 (2018).
  * Kadyrov [1984] M. N. Kadyrov, A mathematical model of the relations between two states, in  _Global Development Processes: Modeling and Analysis_ , Collected Papers, Vol. 3 (Institute for System Studies, Moscow, 1984) pp. 87–99.
  * Wunderling  _et al._ [2020b] N. Wunderling, M. Gelbrecht, R. Winkelmann, J. Kurths, and J. F. Donges, Basin stability and limit cycles in a conceptual model for climate tipping cascades, New Journal of Physics 22, 123031 (2020b).
  * Staal  _et al._ [2018] A. Staal, O. A. Tuinenburg, J. H. C. Bosmans, M. Holmgren, E. H. van Nes, M. Scheffer, D. C. Zemp, and S. C. Dekker, Forest-rainfall cascades buffer against drought across the Amazon, Nature Climate Change 8, 539 (2018).
  * Hays  _et al._ [1976] J. D. Hays, J. Imbrie, and N. J. Shackleton, Variations in the earth’s orbit: Pacemaker of the ice ages, Science 194, 1121 (1976).
  * Sørensen  _et al._ [2020] A. L. Sørensen, A. T. Nielsen, N. Thibault, Z. Zhao, N. H. Schovsbo, and T. W.Dahl, Astronomically forced climate change in the late cambrian, Earth and Planetary Science Letters 548, 116475 (2020).
  * Benzi  _et al._ [1982] R. Benzi, G. Parisi, A. Sutera, and A. Vulpiani, Stochastic resonance in climatic change, Tellus 34, 10 (1982).
  * Gammaitoni  _et al._ [1998] L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Stochastic resonance, Reviews of Modern Physics 70, 223 (1998).
  * Vollmer  _et al._ [2014] J. Vollmer, A. Papke, and M. Rohloff, Ripening and focusing of aggregate size distributions with overall volume growth, Frontiers in Physics 2, 18 (2014).
  * Clark  _et al._ [2011] M. D. Clark, K. S. K, O. J. S, and E. M. Chan, Focusing nanocrystal size distributions via production control, Nano Letters 10.1021/nl200286j (2011).
  * Rosowski  _et al._ [2020] K. A. Rosowski, T. Sai, E. Vidal-Henriquez, D. Zwicker, R. W. Style, and E. R. Dufresne, Elastic ripening and inhibition of liquid–liquid phase separation, Nature Physics 16, 422 (2020).