Universality of noise-induced resilience restoration in spatially-extended ecological systems

Many systems may switch to an undesired state due to internal failures or external perturbations, of which critical transitions toward degraded ecosystem states are prominent examples. Resilience restoration focuses on the ability of spatially-extended systems and the required time to recover to their desired states under stochastic environmental conditions. The difficulty is rooted in the lack of mathematical tools to analyze systems with high dimensionality, nonlinearity, and stochastic effects. Here we show that nucleation theory can be employed to advance resilience restoration in spatially-embedded ecological systems. We find that systems may exhibit single-cluster or multi-cluster phases depending on their sizes and noise strengths. We also discover a scaling law governing the restoration time for arbitrary system sizes and noise strengths in two-dimensional systems. This approach is not limited to ecosystems and has applications in various dynamical systems, from biology to infrastructural systems. Studies on the role of noise in the recovery of a degraded ecosystem have so far been limited to single-variable systems. Here, the authors employ general concepts of nucleation theory in spatially-extended multi-variable systems and apply them to illustrative ecological models.

R esilience, a system's ability to retain its basic functionality when errors and failures occur, is a defining property of many complex systems [1][2][3] . Theoretical models of resilience loss and transitions between alternative states are used to study unanticipated and drastic changes. In real-world systems, many critical transitions are observed, including catastrophic shifts in ecological systems 4,5 , blackouts in power grids 6 , financial crises 7 , climate changes 8 , human depression 9 . These abrupt shifts may arise in the presence of alternative stable states. They can be modeled as critical transitions causing the resilience loss of the desired state, after which the systems switch from a functional to dysfunctional state 10-12 . Ecologists are particularly concerned about sudden critical switches between alternative stable states 4,5 .
As the system approaches the tipping point, its behavior becomes barely predictable when the gradual change of some environmental factors leads to the drastic change of the system state.
Once a system loses its resilience and becomes dysfunctional, restoring the environmental conditions only to those existing before the collapse is often insufficient. Instead, environmental conditions should be recovered to the critical point where the undesired state is destabilized and then resilience restoration would occur. Even though various restoration methods targeted at specific systems have been proposed, stochastic perturbations receive little attention for the recovery of multi-variable systems [13][14][15][16][17][18] . On the other hand, stability and resilience of stochastic dynamical systems containing only one variable has been extensively explored. In particular, it has been shown that noise can induce transitions between alternative stable states and the required time to transition has been established by computing the mean first passage time (MFPT) 2, [19][20][21] . Those studies indicate that single-variable systems can switch between alternative stable states in the presence of noise.
Despite advances in understanding the macroscopic characteristics of resilience restoration, previous research has mostly focused on single-variable or low-dimensional systems, which do not account for an exceptionally large number of variables that in reality are needed to control the state of a complex system. Indeed, many real-world systems consist of numerous components connected via a complex set of weighted, often directed, interactions 22 . The complicated interactions may lead to phenomena not arising in single-variable systems. For example, in the ecosystems, the recovery (or extinction) of a species in one location can impact the states of this species in the neighboring locations, leading to a recovery (or extinction) over the entire system 23 . Accordingly, a full understanding of the system evolution, stability, and resilience cannot be gained without considering interactions among a sufficiently large amount of components. However, hindered by the high-dimensionality of interaction topology and the nonlinear evolution dynamics, few analyses of critical transitions and resilience restoration had been done directly on high-dimensional systems consisting of a great number of participants until the effective reduction theory was recently developed by Gao et al. 3 . This mean-field theory can be used to effectively reduce a multi-dimensional complex system to a one-dimensional system by capturing the average activities of the original system. Furthermore, Liang et al. 24 designed a universal indicator for critical transitions in complex networks and concluded that noise compensates for the structural defects of complex networks, indicating that noise may alter the critical threshold. Jiang et al. 25 studied mutualistic networks through dimensional reduction and claimed that the tipping point can be predicted accurately even in the presence of noise. Tu et al. 26 developed an analytical framework for collapsing complex Ndimensional networked systems into an M + 1-dimensional manifold as a function of M effective control parameters with M ≪ N. While mean-field approaches may guide recovery strategies by indicating the conditions needed to destabilize undesired states, these approaches cannot accurately capture the transition pattern toward the desired state of spatially extended systems in stochastic environments. A recent study by Michaels et al. 27 combined nucleation theory with local-scale positive feedback to understand transitions and resilience in ecological systems. The classical approach to homogeneous nucleation theory was originally developed to describe phase transformation in materials [28][29][30][31][32] , in particular in ferromagnetic 33,34 and ferroelectric systems 35,36 . It was also applied to invasion phenomena in spatial ecological systems [37][38][39][40] . Korniss et al. 38 and O'Malley et al. 39 studied ecological invasion in spatially extended systems with competition. They discriminate between two fundamental modes of nucleating invasive clusters (single-cluster vs. multicluster) and their time-evolution and stochastic features.
In this work, we study the resilience restoration in spatially extended ecosystems, which focuses on a system's ability and required time to recover to its desired state from the undesired one. As noise can induce transitions in single-variable systems, and environmental stochasticity is an inherent property of realworld ecosystems, we investigate resilience restoration by introducing stochastic perturbations into multi-variable systems, providing a theoretical understanding of critical transitions in spatially extended systems subject to environmental stochasticity. Here, we show that nucleation theory builds an elegant bridge between the noise-induced transition and the spread of such transition over the entire system. Our study reveals the crucial effects of noise strength and system size on transition features. Generally, the stronger noise triggers the transition faster, and the larger system takes less time to recover than the smaller system under the same intensity of noise. For large systems or under intense noise, there are multiple clusters nucleated almost simultaneously in the beginning, and these clusters spread out to their neighbors and finally to the entire system. Only one cluster arises for small systems or under weak noise, which expands until the whole system finishes transitions. Another substantial quantity is the time required to recover a degraded system, and it is determined by many factors, including system sizes, noise strengths, and dynamical functions. Our numerical simulations reveal that the transition times are narrowly centered at an average value for the multi-cluster mode, signaling a deterministic feature. While for the single-cluster mode, they vary stochastically for different noise realizations and universally follow an exponential distribution. Also, our study discovers that noise eventually eliminates the deterministic critical threshold, and the recovery of the entire system from the dysfunctional state is possible in the presence of perturbations as long as noise can trigger the transition for just one component. This scenario is likely to occur when the system is close to the deterministic threshold where the undesired state loses stability. The farther away the system is from this point, the more difficult it is for noise to induce the transition.

Results
Transition in single-variable systems and the average lifetime. The current analytical framework for noise-induced transitions is specifically targeted at the low-dimensional system consisting of only a few interacting components (or only one variable) 12,19,[41][42][43] . It would be highly challenging to directly analyze the system consisting of many interacting components (for example, the square lattice). Before the investigation of multi-variable systems, it is helpful to understand the role of noise in single-variable systems and the possibility of transitions between alternative stable states.
For a single-dimensional system, the dynamics may follow the general form as Eq. (1), where f(x, β) governs the deterministic dynamics as the tunable parameter β captures the changing conditions, and η(t) is delta-correlated noise with zero mean and variance hηðtÞηðt 0 Þi ¼ σ 2 δðt À t 0 Þ, where σ is the standard deviation of the noise, which is referred to as noise strength. If the system has a fixed point, x 0 , it satisfies f(x 0 , β) = 0. This fixed point is stable if ∂f ðx;βÞ ∂x x 0 < 0, defined by the linear stability condition. These two equations enable us to derive the resilience function x(β), which represents all the possible steady states of the system as a function of β, as shown in Fig. 1a. This resilience function demonstrates that when the control parameter, β, is >β c 2 , the system only has one stable fixed point, denoted as x H , indicating a resilient state. The system will always recover to the fixed point, x H , for any state perturbations. When the parameter is <β c 1 , the system also has a single fixed point, x L , indicating a collapse of the system. The system is not restorable, unless we change the parameter β. In this work, we particularly focus on the parameter space where β is between β c 1 and β c 2 and close to one of the bifurcation points, such that the system has one unstable fixed point (the dashed line) and two stable fixed points, each with very different sizes of basin of attraction. In the absence of noise, to shift from x L to x H , the environmental conditions would need to be restored all the way to β c 2 , to allow for a transition back to x H . However, noise of the appropriate intensity can trigger the transition without restoring the condition to β c 2 . For example, for a system with β 1 and the system is in the low state, this system can recover to the high state if the noise has a chance to push the system across the unstable fixed point. Furthermore, for two systems (the system 1 and the system 2) with parameters, β 1 and β 2 (β 2 > β 1 ), respectively, it is easier to restore the system 2, because the undesired state of the system 2 is closer to the unstable fixed point. Here, we introduce a normalized variable ρ(t) in Eq. (2) to describe the state behavior during the restoration process. Thus, the normalized stable fixed points are ρ L = 0, ρ H = 1 and the normalized unstable fixed point is ρ u (ρ u ∈ [0, 1]).
In the following analysis, we focus on the case when the desired state ρ H has a much larger basin of attraction than the undesired state ρ L , suitable to examine the "one-way recovery" from ρ L to ρ H . Because of symmetry, the conclusions drawn from this study should be applicable to the transition from ρ H to ρ L when ρ L has a larger basin of attraction. (see Supplementary Note 1 for the detailed comparisons of one-way transitions and stochastic switching between stable states). For the case that we consider throughout the manuscript, the value of the unstable state ρ u is very close to ρ L , increasing the chance for perturbations to push the system from the state ρ L over the unstable fixed point ρ u until the system gets attracted to the state ρ H . It is noteworthy that the backward transition is highly improbable by the same level of perturbations because of the much stronger attraction of the state ρ H . To quantify the time needed for the transition from ρ L to ρ H , the half lifetime, τ (mathematically, τ = {t|ρ(t) = 0.5}), can be safely used as an indicator of the degree of inertia associated with the transition from the undesired state to the desired one. One can choose any value as the cutoff as long as it is sufficiently larger than the unstable state ρ u 38,44 . The lifetime τ is determined by the noise strength σ and the relative stability of alternative stable states controlled by β. Intuitively, larger noise strength indicates stronger fluctuations, which increases the chances of transition, making τ smaller (Fig. 1c). On the other hand, as β increases, the basin of attraction of the stable state ρ H gets larger, Fig. 1 The transition in the single-variable system. a The resilience function of a general bistable system. For the bifurcation parameter β 2 ðβ c 1 ; β c 2 Þ, there are two stable states (x L , x H ) and one unstable state (x u ). The initial stable state (x L ) evolves to the unstable state (x u ) by the aid of noise and is then naturally attracted to the other stable state (x H ) by its deterministic dynamics. b and c display the evolution of the rescaled state ρ in the presence of noise. The time to transition (measured by the lifetime of the initial state), τ, changes with the system parameter β and noise strength σ (defined by the standard deviation of Gaussian noise). b As β 2 is closer to the critical value β c 2 than β 1 , the unstable state ρ u 2 is lower, making the barrier in the landscape easier to cross in the presence of the same strength of fluctuations, and the lifetime τ 2 is thus smaller than τ 1 . c The parameter β is the same, leading to the same landscape. In the presence of stronger noise σ 2 , it is easier to drive the system to get over the barrier, causing τ 2 to be smaller than τ 1 . d shows the simulation results of the average lifetime 〈τ〉 under different values of β and different noise strengths. The dashed line is fitted based on Eq. (4). making the transition from ρ L to ρ H easier so that τ typically decreases (Fig. 1b). To better illustrate the transition process, the underlying landscape picture 45 is introduced, and the effective potential energy is provided by Eq. (3) with zero potential energy at x = 0. In the landscape representation, stable states are traditionally treated as valleys, and unstable states are pictured as hills 4,46,47 . The transition between stable states can be viewed as the transition from one shallow valley to a deeper valley by crossing the barrier.
The transition time from one stable state to the other is a random variable because of the stochastic fluctuations, but the average quantity, 〈τ〉, follows a certain law. Based on the derivations of the Kramers formula for the escape rate over a potential barrier by particles of Brownian motion 48,49 , the average transition time 〈τ〉 can be obtained through the analysis of the MFPT, which is given by Eq. (4). From the quantitative calculation, one can notice that 〈τ〉 increases exponentially with the potential energy difference ΔV = V(x u )−V(x L ) between the low stable state and the unstable state (also interpreted as the barrier height) and decreases with noise intensity σ 2 , which is numerically verified in Fig. 1d. The analysis of single-variable systems provides theoretical support for our intuitive assumption that the low barrier height and strong noise facilitate transitions, leading to resilience restoration (see refs. 42,48 for derivation).
With the knowledge of the average lifetime 〈τ〉 for singlevariable systems available, we are ready to investigate transitions in multi-variable (spatially extended) systems. Nucleation theory is utilized to analyze the generation and spread of the transition in multi-variable systems under external fluctuations.
Spatially extended system and nucleation theory. Generally, the evolution of a system that consists of N coupled components under external perturbations is described by Eq. (5). This study focuses on the spatially extended ecosystem, where one type of species is considered, and its density varies in two-dimensional space and is discretized according to the square lattice topology with periodic boundaries. The deterministic dynamics of each node follows the same self-dynamics F(x i ) and the interaction G(x i , x j ), and parameters of these functions are also set uniformly for all components. The element of the adjacency matrix A is either 0 or a positive value R, which decides the coupling strength between interacting elements. Additionally, to model the external fluctuations acting on the node i, delta-correlated Gaussian noise η i (t) with zero mean and variance hη i ðtÞη j ðt 0 Þi ¼ σ 2 δ ij δðt À t 0 Þ is introduced, which is the same as the noise applied to the singlevariable systems. This general framework can be used to describe a wide range of coupling systems under external perturbations.
Utilizing the dimensional reduction theory 3 , the deterministic evolution of multi-variable systems can be reduced to Eq. (6), containing only one variable, x eff . In the lattice, the fixed states and their stability for each component are proved to be the same as in the reduced one-dimensional system 50,51 . Hence, according to the analysis of single-variable systems, for the specific values of effective interaction strength (β 2 ðβ c 1 ; β c 2 Þ), each component in this system has two stable states and one unstable state. In the presence of noise, the transition from the low state x L to the high state x H is possible for every component. Once the first transition occurs to the node i, its neighbors will also recover through the interaction with it. To show the overall evolution properties of the entire system, the global state ρ(t) is defined in Eq. (7) by taking the average of the individual state value The spread of such transition can be well described by the theory of homogeneous nucleation and growth in finite systems 33,38,39 , and this theory can also predict the spatialclustering pattern formed during the spreading process. The transition from ρ L to ρ H occurs to some node at first, nucleating a cluster, and this cluster continues to expand until it fills the entire space, or the cluster's edge reaches the edge(s) of other clusters that have already nucleated in the system. Homogeneous nucleation makes two assumptions 38 : nucleation occurs in a Poisson process with a constant rate I both temporally and spatially; once a cluster nucleates, it grows homogeneously with a constant radial velocity v. Since the interaction environments for all components are identical, and the perturbations they receive come from the same distribution, each node has the same chance to nucleate a cluster, which satisfies the assumptions of homogeneous nucleation.
As predicted by homogeneous nucleation theory, the restoration process exhibits different patterns for small systems and large systems when the nucleation rate I is fixed. Small systems exhibit the single-cluster pattern because the number of candidates is so limited that the first cluster nucleates and spreads out to the rest of the system before the second possible cluster emerges. Since the nucleation for a specific individual follows a Poisson process, the global state ρ is expected to evolve distinctly for different noise realizations. Also, the time to nucleate the first cluster, t n , and the lifetime, τ, are inherently random. We introduce the waiting time to quantify the time to the system recovery. Because the underlying process is modeled as a random Poisson process, the complementary cumulative probability distribution of waiting time P not is derived as an exponential function, which represents the probability that the global state ρ has not exceeded 0.5 by time t. (Note that our chosen conventional cut-off value 0.5 does not affect the findings.) The distribution of waiting time P not is expressed as where 〈t n 〉 is the average time elapsing until the first cluster nucleates (i.e., the first transition occurs); t g represents the time required for the global state ρ to reach 0.5 after the first cluster emerges, and this time depends on the linear size of the system, t g~N 1/2 /v, where v is the constant radial velocity. Also, 〈t n 〉~(IN) −1 , where I is the nucleation rate per unit area. The average transition time from the initial undesired state to the desired state or the average lifetime of the initial state is expressed as 〈τ〉 = 〈t n 〉 + t g . For small system sizes or in the weaknoise limit, the dominant term in the lifetime is the nucleation time, hence 〈τ〉~(IN) −1 (see Supplementary Note 3 for details). In contrast, for large systems, more than one independent cluster nucleates and expands separately, leading to the multicluster mode. Spatial self-averaging reduces randomness of the global state ρ, making each individual τ closer to the average value and pushing P not closer to a step function. In the large system-size limit, the lifetime distribution is narrowly centered about the average. Based on Avrami's Law, for sufficiently large systems [28][29][30][31]35,36 , the evolution of ρ can be expressed in a deterministic form as ρ ¼ 1 À e À πv 2 I 3 t 3 . By setting ρ = 0.5 (without loss of generality), the average lifetime for the multi-cluster mode can be obtained as hτi ¼ 3ln 2 πv 2 I À Á 1=3 . Then the evolution of ρ can be described by Resilience restoration of mutualistic systems. To verify the predictions by nucleation theory about noise-induced transition patterns in systems with alternative stable states, we first study resilience restoration of the mutualistic system. We use Eq. (10) as the deterministic dynamics to track the abundance of one species distributed among a square lattice in the mutualistic system 52,53 . The self-dynamics F(x i ) describes that the growth of the species in each location follows the logistic law with the Allee effect, and the dynamics G(x i , x j ) accounts for the mutualistic interaction between the species in two neighboring locations i and j through the interaction strength A ij defined in Eq. (5).
Note that we use the same parameters as Ref. 3 . The parameters are node-uniform and set as , and the interaction strength R = 1 (leading to the effective interaction strength β eff = 4). At some moment, if the species in all locations are trapped in the low stable state x L , they will stay at this state forever if there is no action or perturbation. It is, for sure, not desired from the ecological viewpoint. If one would like to keep the system always in the high stable state x H , a straightforward approach is to increase the interaction strength to ensure that β eff is larger than the critical bifurcation value β c 2 . In this case, the system will be attracted to the high state no matter where it starts. Alternatively, noise has shown the ability to induce the transition between two stable states from the study of one-variable systems. We are particularly interested in how noise assists the resilience restoration (i.e., transferring the system from the undesirable state x L to the desired state x H ).
Let us consider the case when all nodes in the system start from the low state x L , indicating a collapsed state. Observed from simulations, the proper noise can excite some node to x H , and such transition spreads out to its neighbors via interaction until the rest of the system completes transitions. According to homogeneous nucleation theory, for different system sizes and nucleation rates, two possible transition patterns are present. We successfully show the snapshots of two cluster modes by numerical simulations, the single-cluster and multi-cluster mode. One can notice that these two modes possess radically distinct properties. For the single-cluster mode (Fig. 2a), there is only one node that switches from x L to x H in the beginning, and the rest of nodes shift to x H through the interactions with the neighbors which have already transitioned; while for the multi-cluster mode (Fig. 2e), there is more than one node in the separate location that transfers to x H simultaneously. This is expected as for the large system, there are enough candidates to receive fluctuations, increasing the chance to induce independent transitions.
As predicted by nucleation theory, the evolution of ρ(t) for the single-cluster mode and multi-cluster mode differs a lot, and numerical results confirm the difference. Figure 2b, f display 100 realizations for system sizes N = 100 and N = 10,000 under the same intensity of noise. For the single-cluster mode (Fig. 2b), the evolution varies for individual realizations, so the transition times are different, indicating the uncertain evolution feature. In contrast, for the multi-cluster mode (Fig. 2f), the evolution of ρ(t) is similar for different realizations, implying that the evolution is deterministic in the infinite system-size limit. Following this, the waiting time distribution P not for two cluster modes is also verified (Fig. 2c). For the single-cluster mode, P not is initially constant and then decreases exponentially with respect to t verifying Eq. (8). The slope of the distribution gets more negative as noise becomes stronger suggesting a larger nucleation rate. For the multi-cluster mode, P not gets closer to a step function (Fig. 2g) as noise strength increases. This is because the larger nucleation rate induces more separate clusters for a given system size and thus leads to the more deterministic evolution by selfaveraging. From the theoretical analysis, the global state ρ for the multi-cluster mode evolves predictably according to Eq. (9). Whereas, for the finite-size system, the evolution of the multicluster mode (Fig. 2h) is not perfectly deterministic, but still much less random than the single-cluster mode (Fig. 2d).
The cluster mode not only depends on the system size but also relies on the nucleation rate, which is decided by noise strength. Typically, large systems under strong noise belong to the multicluster mode, while small systems under weak noise exhibit the single-cluster mode. However, low nucleation rates resulting from weak noise can induce the single-cluster mode even for a very large system (Fig. 3a, b). Also, the average nucleation time 〈t n 〉 = (NI) −1 for the single-cluster mode is validated. The average lifetime 〈τ〉 behaves differently for the two cluster modes and exhibits two regimes (Fig. 3c), and its value increases exponentially as σ −2 increases for both cluster modes as predicted by Eq. (15). For the given dynamics, 〈τ〉 entirely depends on the system size N and noise strength σ (Fig. 3c, d), and there is a decrease as N or σ increases. One can also notice that the slope of ln hτi as a function of σ −2 for the single-cluster mode is larger than the slope for the multi-cluster mode.
In accord with the scaling theory (see the "Methods" section, Eqs. (15)-(17)), the two distinct cluster-growth modes are separated by a crossover region centered around N 1/2~R 0 (Fig. 4a). One should note that this crossover (centered around the dashed curve in Fig. 4a) is not a sharp transition separating the two cluster-growth modes. The gradual change of the background color (red-white-blue) is to qualitatively illustrate the continuous nature of the crossover from the single-cluster mode to the multi-cluster mode. The small system (Fig. 4b) or weak noise induces (Fig. 4c) the single-cluster mode, while the large system or relatively strong noise (Fig. 4d, e) produces the multi-cluster mode. According to the proposed scaling function, employing and plotting properly scaled variables, hτie À c 3σ 2 vs. e c 3σ 2 =N 1 2 , we expect that all numerical data would collapse onto a single curve, capturing general nucleation behavior 54 .
Nevertheless, observed in Fig. 5a, for the multi-cluster mode, 〈τ〉 is not precisely proportional to I −1/3 . Therefore, the two transition modes cannot be scaled in a satisfactory fashion as the scaled data for the multi-cluster regime is not constant (Fig. 5b), which contradicts the assumption of Eq. (16) (see the "Methods" section). The deviation from Avrami's Law suggests that the assumption(s) of homogeneous nucleation might be violated. For a large system subjected to relatively strong noise, which guarantees the multi-cluster mode, the nucleation rate changes in the beginning (Fig. 5c), defying the assumption of constant nucleation rate. The initial rise of the spatially distributed low state for the nodes which have not been driven to the high state  2), and it is unitless) from the undesired state (ρ i = 0) to the desired state (ρ i = 1). Initially, all of the nodes are at the low state. At some time later, the transition to the high state occurs to one node, which is treated as a single cluster. Such transition then spreads out to its neighbors. b The evolution of the global state ρ for 100 realizations of the single-cluster mode. c The probability distribution of waiting time P not for the fixed system size and various noise strengths σ = 0.08, 0.09, 0.1. The single dots are simulation data, and different types of lines are obtained by linear fit according to Eq. (8). d The evolution of the average global state ρ using the same data as in c. e One snapshot shows the state ρ i of each node for σ = 0.1. Different from a, the transition from the initial low state to the high state occurs at several separate nodes, and they expand independently, forming the multiple cluster. f The evolution of the global state ρ for 100 realizations, which are more centered around a certain value instead of being random in b. g The distribution P not for σ = 0.08, 0.09, 0.1, which approaches the step function as σ increases. h The evolution of the global state ρ averaged over 100 realizations using the same data as in g. The single dots are simulation data, and different types of lines are obtained by linear fit according to Eq. (9). indicates that the initial state is not the precise metastable state. In the presence of noise, the metastable state is a state distribution where all nodes are around the low state ρ L , but some nodes are closer to the unstable state ρ u and some are even below the state ρ L . The system in the presence of noise needs time to reach the actual metastable state which is different from the mean-field low-stable state. To keep the nucleation rate constant, one should use the metastable configuration as the initial condition. However, as seen in Fig. 5c, some cluster have already emerged before the system reaches the metastable configuration indicated by the increase of the average low state and the nucleation rate, so that the preprocessing of the initial state is needed to forbid any transition before the metastable configuration is constructed. One can artificially prepare the system close to the metastable state by reverting the node back to ρ L if any other sites nucleate during preprocessing.
To gain further insight into the source of this discrepancy, we carried out simulations with the initial configurations being very close to the metastable state. In such a case, the nucleation rate I stabilizes much faster (Fig. 5g) compared with the system without preprocessing. Afterward, the nucleation rate is roughly constant, so homogeneous nucleation theory can be reliably applied. This leads to a better agreement between the simulation and the theory, as the evolution of ρ is closer to what Eq. (9) postulates (Fig. 5d, h), and the data from the two cluster modes follows the scaling function (Fig. 5e, f). It is expected that the average lifetime 〈τ〉 agrees better with Eq. (15) if the metastable state can be perfectly prepared.
Restoration of diffusion dynamics. To further validate the proposed theory, we adapt three well-studied ecological models exhibiting alternative stable states 55,56 with diffusive interaction and then apply the above theory to investigate the transition features. The self-dynamics for three diffusion models are defined below.
The harvesting model in Eq. (11) describes the growing resource biomass with fixed grazing rate 10 . The first term on the right-hand side of Eq. (11) describes logistic growth, where r is the maximum growth rate, and K is the carrying capacity. The second term is the "Holling's type III" consumption function 57 . The system transitions from an underexploited state to overexploited state as the harvesting rate β exceeds a certain critical value.
The eutrophication model in Eq. (12) describes the dynamics of nutrient concentration in the eutrophic lake 58 . The variable x i represents the density of phosphorus mass (nutrient) in the location i of the lake. The first term a on the right hand side of Eq. (12) is the nutrient loading rate, the second term describes nutrient loss processes with the rate r, and the last term accounts for recycling processes following a sigmoid function. As the maximum recycling rate β increases to the critical point, the lake transfers from oligotrophic to eutrophic.
The vegetation-turbidity model in Eq. (13) describes the vegetation dynamics considering turbidity 59 . The variable x i represents the density of aquatic vegetation in the location i of the lake. The first term on the right-hand side of Eq. (13) characterizes the growth of vegetation with the maximum growth rate r v . The function E is an inverse Monod function and used to describe the vegetation effect on turbidity. Accordingly, the second term is a Hill function describing the sigmoidal decline of vegetation with turbidity. As the water becomes turbid, indicated by the background turbidity β, macrophytes suddenly decrease.
All three models can exhibit alternative stable states with the properly chosen bifurcation parameter β as illustrated in Fig. 7a, e, i. The transition from the low stable state to the high stable state is possible in the presence of noise for such single-variable systems. Likewise, the underlying topology is also a square lattice with periodic boundaries. The interaction dynamics is defined in Eq. (14), which represents the diffusive process between adjacent neighbors, and the interaction strength R determines the uniform diffusion rate. The species density at each location varies according to the internal dynamics as defined in Eqs. (11)-(13), and it is also influenced by dispersion to or immigration from neighbors and stochastic environmental fluctuations.
Once the system degrades to a malfunctioning state, resilience restoration is required. In the presence of noise, each component is likely to be driven from the undesired state to the functional stable state, and then the entire system undergoes substantial changes due to the transition of one or a few nodes. Nucleation theory can be employed as well to study the overall transition features. The results we collected from three diffusion models are is to qualitatively illustrate the continuous nature of the crossover from the single-cluster mode to the multi-cluster mode. The dashed curve corresponds to the center of the crossover region (provided by the above formula), separating the two cluster-growth modes. b and c describe the single-cluster mode (in red), while d and e display the multi-cluster mode (in blue). b The distribution P not for the fixed noise strength σ = 0.1 and different system sizes N = 9, 16, 25, 36, 49, 64, 81, 100. c P not for the fixed weak noise σ = 0.06 and N = 100, 900, 2500, 10,000. In both b and c, the single dots are simulation data, and different types of lines are obtained by linear fit according to Eq. (8). d P not for the fixed noise strength σ = 0.1 and different system sizes N = 900, 2500, 10,000. e The evolution of the global state ρ using the same data as in d. similar to those obtained from the mutualistic system and thus verify the predictions by nucleation theory. We consider the case when the undesired state has a much smaller basin of attraction than the desired state so that noise can drive the system from the initial undesired state to the functional state, resulting in resilience restoration. Take the harvesting model as an example to illustrate the successful application of nucleation theory to resilience restoration. Two cluster-growth modes are observed and separated by a crossover around N 1/2~R 0 (Fig. 6a). Similar to the mutualistic system, the large system or strong noise produces multi-cluster mode; conversely, the small system or weak noise induces singlecluster mode. For the small system (N = 100 in the example), there is only one cluster formed during the transition, exhibiting the single-cluster pattern. The individual lifetime τ in Fig. 6b varies a lot for different realizations, indicating a stochastic feature. As derived above, the distribution of waiting times P not follows an exponential function after a certain period t g , observed in Fig. 6c,d. If the system is exposed to weak noise which means a low nucleation rate, it is very likely to enter the single-cluster regime even the size is sufficiently large (Fig. 6d). For the large system (like N = 10,000) with proper intensity of noise, more than one node in the separate location is recovered simultaneously, presenting a multi-cluster pattern. Spatial-averaging reduces randomness so that the individual lifetime τ is more centered about a certain value (Fig. 6e). The evolution is much more deterministic than the system in a single-cluster mode. The global state ρ evolves approximately as Eq. (9) predicts (Fig. 6f). If the applied noise is strong enough, the moderate-sized system can still exhibit the multi-cluster mode (Fig. 6g).
As seen in Fig. 7, the system size and noise strength decide the recovery time and the cluster mode. Similarly, the average lifetime 〈τ〉 of three diffusion models displays two distinct regimes. One is the single-cluster and the other is the multi-cluster regime (Fig. 7a, d, g). The slope of ln hτi versus σ −2 reveals which mode is active. If the system starts from the undesired stable state, the scaling between two cluster modes deviates a little from the theoretical prediction (Fig. 7b, e, h). If the system starts from the prepared state, resembling the metastable configuration, the scaling agrees well with the designed scaling function (Fig. 7b, e, h). Overall, the conclusions of the single-cluster and multi-cluster transition are validated in three diffusion models.

Discussion
We have utilized nucleation theory to analyze the noise-induced resilience restoration in ecosystems where the desired state has a much larger basin of attraction than the undesired state has. This is a general theory, and we successfully apply it to four ecological models, revealing the transition features. During the restoration process, homogeneous nucleation theory distinguishes two different cluster modes: the single-cluster and multi-cluster transition modes. We also derive the formulas for the recovery time under different conditions and propose a scaling function that collapses all the data onto one universal line.
The two cluster modes possess quite distinct features. The individual lifetime is random for the single-cluster phase, and the waiting time before the transition follows an exponential distribution. In contrast, for the multi-cluster mode, the lifetime is less random and centered about its average value so that the evolution of the global state ρ is more deterministic. Which cluster mode the system follows is decided by its size and noise strength, and the crossover region is theoretically derived and can separate two phases. Generally, the large system subjected to strong noise presents the multi-cluster mode, and the small Fig. 5 Scaling between two cluster modes in the mutualistic system. For a-c, the system starts from the low state x L , and the time to reach the high state is characterized by τ. a The relationship of 〈τ〉 and e À c σ 2 differs between two cluster modes, where c is a constant specific to the dynamics and σ is the standard deviation of noise. a, b, e, and f share the same legend. b The finite-size scaling is drawn by assuming the slope of multi-cluster mode in a is À 1 3 . For c and d, the system size is set as N = 10,000 and the noise strength is σ = 0.08. c The nucleation rate increases before the average low-state nodes stabilizes. d The system starts to evolve from x L and the prepared state, respectively. The single dots are simulation data, and different types of lines are obtained by linear fit according to Eq. (9). The related uncertainty can be calculated by the standard deviation std ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 nÀ2 ∑ n i¼1 ðlog ð1 À ρ fitted Þ À log ð1 À ρ simulation ÞÞ 2 q between the fitted state and the simulated state, where n is the number of data for the linear fit. The uncertainties are 0.42 and 0.09 for the initial state being x L and the prepared state, respectively. With the slope being fitted, the global state evolution for the latter case agrees better with Eq. (9) than the former case. For e-h, the system starts from the prepared metastable state, and the time to reach the state x H is recorded as τ. e The average lifetime 〈τ〉 for two cluster modes. f The finite-size scaling is consistent with the theoretical prediction in Eq. (17). g The nucleation rate takes less time to stabilize compared with c. h The evolution of the global state ρ for the multi-cluster mode when N = 10,000 and σ = 0.08, 0.085, 0.09, 0.095, 0.1. The uncertainties of linear fit are std = 0.09, 0.14, 0.13, 0.13, 0.15, respectively.
system with weak noise displays the single-cluster regime. One quantity of interest for resilience restoration is the recovery time, which also depends on the system size and the noise strength. The rise of noise strength increases the nucleation rate, diminishing 〈τ〉. There is no finite-size effect when noise is strong enough; that is, 〈τ〉 is the same for different system sizes. The decrease of noise strength reveals the size effect, where larger systems take less time to complete transitions than smaller systems on average. Due to the distinct evolution features, the relationship between 〈τ〉 and σ varies for two cluster modes. What can be clearly seen in Fig. 3 and Fig. 7a, d, g, is that 〈τ〉 exhibits two distinctive regimes, corresponding to two cluster modes. The scaling between two cluster modes is proposed according to Avrami's Law. The deviation observed in numerical simulation can be corrected by preparing the initial state close to metastable configurations to satisfy two homogeneous nucleation assumptions. Employing the nucleation theory, we successfully extend the noise-induced transition in single-variable systems to spatially extended multi-variable systems. Our framework is useful to predict critical transitions and guide the resilience restoration in the general dynamical systems presenting alternative stable states. One may wonder how the mean-field theory works in noisy environments as the noise-induced transition in single-variable systems has been well explored. The mean-field approach makes a satisfactory prediction of the system state when the system is close to the stable states. Unfortunately, it cannot capture the transition features observed in spatially-extended systems, including the average lifetime, the dependence on the system size, and the spatial-clustering patterns (see Supplementary Note 2 for the results of mean-field theory and the comparisons with the nucleation approach.). Therefore, it is still an open question whether the mean-field theory can be used to study transitions in noisy environments. One may need to develop other dimension reduction approaches for this problem. Also, there are further questions to be addressed. For example, we analyze the transition in the lattice model. In reality, the interaction relationship in ecosystems exhibits various structures. In addition, the Fig. 6 Single-cluster and multi-cluster modes in the harvesting system. The parameters are set as r = 1 (defined as the maximum growth rate), K = 10 (defined as the carrying capacity), the diffusion rate R = 0.02, and the bifurcation parameter β = 1.8. Initially, all of the nodes are at the low stable state ρ L , and they are driven to the high stable state ρ H in the presence of noise. a Two insets describe the snapshots of individual normalized state values ρ i for the single-cluster and multi-cluster mode, respectively. They are separated by a crossover region centered around the curve N 1=2 ¼ e c 3σ 2 , where c is a constant decided by the dynamics, and σ is the standard deviation of noise (noise strength). As in Fig. 4, the gradual change of the background color (red-whiteblue) is to qualitatively illustrate the continuous nature of the crossover from the single-cluster mode to the multi-cluster mode. b-d are about the singlecluster mode. b 100 realizations of the global state evolution (state ρ changes with time t) for the system size N = 100 and the noise strength σ = 0.045. c The distribution of waiting time P not for the fixed system size N = 100 and the varied noise strengths σ = 0.04, 0.043, 0.045. d P not for the fixed weak noise σ = 0.035 and different system sizes N = 100, 900, 2500, 10,000. e-g correspond to the multi-cluster mode. e 100 realizations of the global state evolution for the system size N = 10,000 and the noise strength σ = 0.045, which are more centered around a certain value instead of being random in b. f The evolution of ρ averaged over 100 realizations for different noise strengths σ = 0.04, 0.043, 0.045. g The evolution of the global state ρ averaged over 100 realizations for the fixed strong noise σ = 0.045 and large system sizes N = 900, 2500, 10,000.
interaction strength between components in the real-world complex systems varies, which may add to the difficulty of analysis. Furthermore, we focus on the one-way transition, which requires that the system be close to the theoretical bifurcation point where the basin of attraction of the undesired state vanishes. For the case when two alternative stable states have basins of attraction of similar sizes, stochastic switching (namely back-and-forth switching) between states may arise, which is beyond the scope of this study. In summary, network topology (in addition to spatial structure), coupling strength, and stochastic switching are of great interest to be investigated in future studies.

Methods
Scaling law of single-cluster and multi-cluster modes. According to Avrami's Law and homogeneous nucleation in finite systems 33,38,39 , the average lifetime of two transition modes is summarized as hτi $ where R 0 $ ð v I Þ 1=3 $ I À1=3 is the typical distance between separate clusters (and N 1/ 2 is the linear size of the two-dimensional lattice). Transition patterns for different system sizes and nucleation rates are classified into two distinct cluster nucleation modes, separated by a crossover region centered around the curve, N 1/2~R 0 . Small systems or low nucleation rates induce the single-cluster mode, while large systems or high nucleation rates exhibit the multi-cluster mode.
By constructing a scaling function 54 with the following asymptotic behavior: SðuÞ $ u 2 ; u ) 1 ð single-cluster mode Þ const:; where u = R 0 /N 1/2~v1/3 /(I 1/3 N 1/2 )~I −1/3 N −1/2 , one can capture the average lifetime of any system size and nucleation rate (including the crossover between the single-cluster and multi-cluster regimes), Motivated by the study on the average transition time 〈τ〉 for single-variable systems [Eq. (4)], the relationship between local nucleation rate and noise strength in spatially extended systems is expected to scale as where c is a constant specific to the given dynamics and can be empirically fitted in the weak-noise case. In turn, the average lifetime 〈τ〉 for different system sizes and noise strength can be described by a universal scaling function: employing Eqs. (17) and (18), when plotting hτie À c 3σ 2 vs. e c 3σ 2 =N 1 2 , all curves should collapse on the same curve S(u).

Data availability
The authors declare that data supporting the findings of this study are available within the paper.
Code availability Fig. 7 Average lifetime 〈τ〉 and the scaling results for three diffusion models. The diffusion rate is set as R = 0.02. a-d harvesting model, r = 1 and K = 10; e-h eutrophication model, a = 0.5 and r = 1; i-l vegetation model, r = 1, r v = 0.5, h v = 0.2. a, e, and i are the resilience diagrams. These models exhibit alternative stale states for a β ∈ (1.79, 2.60), e β ∈ (0.86, 6.35), and i β ∈ (2.59, 3.64). The arrows represent transitions from the low stable state to the high stable state. In simulation, the bifurcation parameter β takes values of 1.80, 6.00, 2.60 for a, e, and i, respectively. For b, f, and j, the system starts from the undesired states x L , and the average lifetime 〈τ〉 changes with N and σ. For e, g, and k, the system starts from the undesired states x L , the lifetimes of two clusters and the crossover are shown in the scaling form. For d, h, and l, the system starts from the prepared metastable states, the lifetimes of two clusters and the crossover are shown in the scaling fashion.