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 a prominent example. Resilience restoration focuses on the ability of spatially-extended systems and the required time to recover to their desired states under stochastic environmental conditions. While mean-field approaches may guide recovery strategies by indicating the conditions needed to destabilize undesired states, these approaches are not accurately capturing the transition process toward the desired state of spatially-extended systems in stochastic environments. The difficulty is rooted in the lack of mathematical tools to analyze systems with high dimensionality, nonlinearity, and stochastic effects. We bridge this gap by developing new mathematical tools that employ nucleation theory in spatially-embedded systems to advance resilience restoration. We examine our approach on systems following mutualistic dynamics and diffusion models, finding that systems may exhibit single-cluster or multi-cluster phases depending on their sizes and noise strengths, and also construct a new scaling law governing the restoration time for arbitrary system size and noise strength in two-dimensional systems. This approach is not limited to ecosystems and has applications in various dynamical systems, from biology to infrastructural systems.


I. INTRODUCTION
Resilience, 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 functional to a dysfunctional state [10][11][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. 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 after the transition to the undesired state. Even though various restoration methods targeted at specific systems have been proposed, stochastic perturbations receive little attention [13][14][15][16][17][18].
On the other hand, stability of stochastic dynamical systems containing only one variable has been extensively explored. Also, it has been known for a long time that noise can greatly affect the stability and resilience of the system with bistable states. 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 the single-variable systems can switch between alternative stable states in the presence of noise. One may expect that noise can also induce transitions in multi-variable systems. However, one should keep in mind that such noise-induced transition occurs only when the system is close to the bifurcation point where the basin of attraction of the current stable state vanishes. If the system is far from this bifurcation point, noise may drive the system back and forth between alternative stable states, or even not be able to trigger any transition. Since environmental stochasticity is an inherent property of real-world ecosystems, we investigate resilience restoration after introducing stochastic perturbations into multi-variable systems, providing a theoretical understanding of critical transitions in spatially-extended systems subject to environmental stochasticity.
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 the 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 number 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 theory enables us 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 N -dimensional networked systems into an M + 1-dimensional manifold as a function of M effective control parameters with M N .
Nonetheless, our study shows that noise 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 very 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.
The time required to recover a system is a quantity of great interest but determined by many factors, including system sizes, noise strengths, and dynamical functions. Nucleation theory provides an elegant bridge between the noise-induced transition and the spread of such transition over the entire system. The classical approach to homogeneous nucleation theory was originally developed to describe phase transformation in materials [27][28][29][30][31], in particular in ferromagnetic [32,33] and ferroelectric systems [34,35] [32].
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 compared with the smaller system under the same intensity of noise. For large systems or under strong noise, there are multiple clusters nucleated simultaneously in the beginning, and these clusters spread out to their neighbors and to the whole system.
Our numerical simulations reveal that the transition time is narrowly centered at an average value, signaling a deterministic feature. For small systems or under weak noise, one cluster is generated first, which expands until the entire system finishes transitions. The transition times vary stochastically for different realizations of noise, and they universally follow an exponential distribution.

II. MATHEMATICAL FRAMEWORK
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 singlevariable systems and the possibility of transitions between alternative stable states.

A. Single-variable system
For a single-dimensional system, the dynamics may follow the general form as In Eq. (1), 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 where σ is the standard deviation of the noise, which is referred to as noise strength. If the system has a fixed point, < 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 greater than β 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 less than β 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 manuscript should be applicable to the transition from ρ H to ρ L when ρ L has a larger basin of attraction.
(see Supplementary Information, Sec I for the detailed comparisons of one-way transitions and stochastic switching between stable states). For the case that we consider throughout the manuscript, the basin of attraction of the state ρ L is much smaller than that of ρ H .
In this case, 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 . Notably, 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 [37,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, 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 value, τ , presents interesting properties.
Following 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 calculated through the analysis of the mean first passage time (MFPT), which is given by Eq. (4).
From the quantitative predictions, 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 [42,48] for derivation).
With the knowledge of the average lifetime τ for single-variable 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.

B. Spatially-extended system
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 same interaction dynamics 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 η i (t)η j (t ) = σ 2 δ ij δ(t − t ) is introduced, which is the same as the noise applied to the single-variable 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 multivariable 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]. Hence, according to the analysis of single-variable systems, for the specific values of effective interaction strength (β ∈ (β 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 [32,37,38], and this theory can also predict the spatial-clustering 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 nucleated in the system. Homogeneous nucleation makes two assumptions [37]: 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 the 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, 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 weak-noise limit, the dominant term in the lifetime is the nucleation time, hence In contrast, for large systems, more than one independent cluster nucleates and expands separately, leading to the multi-cluster 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 [27-30, 34, 35], 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 τ = ( 3 ln 2 πv 2 I ) 1/3 . Then the evolution of ρ can be described by According to Avrami's Law and homogeneous nucleation in finite systems [32,37,38] , the average lifetime of two transition modes is summarized as R 0 (single-cluster mode) 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 [51] with the following asymptotic behavior, , 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 weaknoise case. In turn, the average lifetime τ for different system sizes and noise strength can be described by a universal scaling function: employing Eqs. (12) and (13)

III. RESILIENCE RESTORATION OF MUTUALISTIC SYSTEMS
To verify the predictions by nucleation theory of the noise-induced transition patterns in systems with alternative stable states, we first study resilience restoration of the mutualistic system.

A. Dynamical model
We use Eq. (14) 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 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 β c2 . 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 ).

B. Patterns of resilience restoration
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 success-fully show the snapshots of two cluster modes by numerical simulations, the single-cluster and multi-cluster mode. Notably, 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. It is expected since 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. Fig. 2b and 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 self-averaging. 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 multi-cluster mode (Fig. 2h) is not perfectly deterministic, but still much less random than the single-cluster mode (Fig. 2d).

C. The role of system size and noise strength
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 multi-cluster mode, while small systems under weak noise exhibit the singlecluster mode. However, low nucleation rates resulting from weak noise can induce the single-cluster mode even for a very large system (Fig. 3a, 3b). Also, the average nucleation time t n = (N I) −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. (10). For the given dynamics, τ entirely depends on the system size N and noise strength σ (Fig. 3c,   3d), and there is a decrease as N or σ increases. One can also notice that the slope of ln τ as a function of σ −2 for the single-cluster mode is larger than the slope for the multi-cluster mode.

D. Phase diagram and finite-size scaling
In accord with the scaling theory, Eqs. (10)-(12), 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 (redwhite-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

E. The universal scaling law
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. (11). 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 indicates that the initial state is not a 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 a new metastable state which is different from the 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, 5h), and the data from the two cluster modes follows the scaling function (Fig. 5e, 5f). It is expected that the average lifetime τ agrees better with Eq. (10) if the metastable state can be perfectly prepared.

IV. RESTORATION OF DIFFUSION DYNAMICS
To further validate the proposed theory, we adapt three well-studied ecological models exhibiting alternative stable states [54,55] with diffusive interaction and then apply the above theory to investigate the transition features.

A. Examples of some dynamics
The self-dynamics for three diffusion models are defined below.
The harvesting model in Eq. (15) describes the growing resource biomass with fixed graz-ing rate [10]. The first term on the right hand side of Eq. (15) 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 [56]. The system transitions from an underexploited state to overexploited state as the harvesting rate β exceeds a certain critical value.
The eutrophication model in Eq. (16)  in the location i of the lake. The first term a on the right hand side of Eq. (16) 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. (17) All three models can exhibit alternative stable states with the properly chosen bifurcation parameter β as illustrated in Fig. 7a, 7e, 7i. 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. (18), 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. (15) - (17) , and it is also influenced by dispersion to or immigration from neighbors and stochastic environmental fluctuations.
Once the system gets stuck in 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 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.

C. The phase diagram of Harvesting system
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 single-cluster mode. For the small system (N = 100 in the example), there is only one cluster formed during the transition (Fig. 6b), exhibiting the single-cluster pattern.
The individual lifetime τ in Fig. 6c 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. 6d, 6e. 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. 6e). For the large system (like N = 10000) with proper intensity of noise, more than one node in the separate location is recovered simultaneously (Fig. 6f), presenting a multi-cluster pattern. Spatial-averaging reduces randomness so that the individual lifetime τ is more centered about a certain value (Fig. 6g). The evolution is much more deterministic than the system in a single-cluster mode. The global state ρ evolves approximately as Eq. (9) predicts (Fig. 6h). If the applied noise is strong enough, the moderate-sized system can still exhibit the multi-cluster mode (Fig. 6i).

D. The universal scaling in diffusion dynamics
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, 7d, 7g). The slope of ln τ 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, 7e, 7h). If the system starts from the prepared state, resembling the metastable configuration, the scaling agrees well with the designed scaling function (Fig. 7b, 7e, 7h). Overall, the conclusions of the single-cluster and multi-cluster transition are validated in three diffusion models.

V. 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 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 finitesize 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, 7d, 7g, 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

DATA AVAILABILITY STATEMENT
The authors declare that data supporting the findings of this study are available within the paper.    Fig. 4, 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 dashed curve corresponds to the center of the crossover region (provided by the above formula), separating the two cluster-growth modes. For (b) -(d), the system size N = 100 (10 × 10 lattice). (b) One snapshot shows the evolution of each node ρ i in the presence of noise of the standard deviation σ = 0.045. The initial states for all of the nodes are x L . At some time later, the transition to x H occurs to one node, which is treated as a single cluster, and it spreads out to its neighbors. (c) The evolution of the global state ρ for 100 realizations, which corresponds to the single-cluster mode in (a). (d) The distribution of waiting time P not for the fixed system size and the varied noise strengths σ = 0.04, 0.043, 0.045. (e) P not for the fixed weak noise σ = 0.035 and different system sizes N = 100, 900, 2500, 10000. In both (d) and (e), the single dots are simulation data, and different types of lines are obtained by linear fit according to Eq. (8). For (f) -(h), the system size N = 10000 (100 × 100 lattice). (f) The snapshot illustrates the evolution of the system starting from x L , where σ is also 0.045. Different from (b), the transition to x H occurs at several separate nodes, and they expand independently, forming the multiple cluster. (g) 100 realizations of the global state ρ, which are more centered around a certain value instead of being random in (c). (h) The evolution of ρ averaged over 100 realizations for σ = 0.04, 0.043, 0.045. The single dots are simulation data, and different types of lines are obtained by linear fit according to Eq. (9). (i) The evolution of the global state ρ averaged over 100 realizations for the fixed strong noise σ = 0.045 and large system sizes N = 900, 2500, 10000. 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.