Population boundary across an environmental gradient: Effects of quenched disorder

TL;DR

This study models how local heterogeneities affect ecological population boundaries across environmental gradients, revealing intermittent dynamics and universal scaling laws that could obscure climate change effects in observations.

Contribution

It introduces a disordered contact process model with strong-disorder renormalization to analyze population boundary dynamics under environmental gradients, highlighting the impact of quenched heterogeneities.

Findings

Population fronts advance intermittently with long quiescent periods.

Scaling laws relate the dynamics to the environmental gradient and correlation-length exponent.

Observations may underestimate long-term climate change effects due to intermittent boundary movements.

Abstract

Population boundary is a classic indicator of climatic response in ecology. In addition to known challenges, the spatial and dynamical characteristics of the boundary are not only affected by the spatial gradient in the environmental factors, but also by local heterogeneities in the regional characteristics. Here, we capture the effects of quenched heterogeneities on the ecological boundary with the disordered contact process in one and two dimensions with a linear spatial trend in the local control parameter. We apply the strong-disorder renormalization group method to calculate the sites occupied with an $O(1)$ probability in the stationary state, readily yielding the population front's position as the outermost site locally as well as globally for the entire boundary. We show that under a quasistatic change of the global environment, mimicking climate change, the front advances intermittently: long quiescent periods are interrupted by rare but long jumps. The characteristics of this intermittent dynamics are found to obey universal scaling laws in terms of the gradient, conjectured to be related to the correlation-length exponent of the model. Our results suggest that current observations might misleadingly show little to no climate response for an extended period of time, concealing the long-term effects of climate change.

Tables

Table 1. Table 1: Critical exponents ν ⟂ subscript 𝜈 perpendicular-to \nu_{\perp} and α 𝛼 \alpha for the clean (CP) and disordered contact process (DCP) in one and two dimensions.

	1d CP	1d DCP	2d CP	2d DCP
${\nu }_{⟂}$	$1.096854\left(4\right)$jensen	$2$ hiv	$0.733\left(4\right)$ gz ; voigt	$1.24\left(2\right)$ kovacs
$\alpha$	$0.523095\left(1\right)$	$2/3$	$0.423\left(1\right)$	$0.554\left(4\right)$

Equations

$$\mu (x)=p[1-g(x-x_c)],$$

$$\overline{\Delta }(x)=\overline{\ln [\mu (x)/\lambda ]}=\ln \mu (x)-\overline{\ln \lambda },$$

$$\overline{\Delta }(x)=\overline{\Delta }(x_c)-g(x-x_c)+O\{[g(x-x_c){]}^2\}.$$

$$\ln \tilde{\mu }=\ln {\mu }_i+\ln {\mu }_j-\ln {\lambda }_{ij}+\ln 2.$$

$$\ln {\tilde{\lambda }}_{jk}=\ln {\lambda }_{ij}+\ln {\lambda }_{ik}-\ln {\mu }_i.$$

$$\xi [\overline{\Delta }(x)]\sim |\overline{\Delta }(x)-{\Delta }_c{|}^{-{\nu }_{\perp }}\sim |g(x-x_c){|}^{-{\nu }_{\perp }},$$

$${\ell }_x\sim g^{-\alpha },$$

$$\alpha =\frac{{\nu }_{\perp }}{1+{\nu }_{\perp }}.$$

$$P(l)\sim e^{-\mathrm{const}\cdot g^{{\nu }_{\perp }}{l}^{1+{\nu }_{\perp }}}\sim e^{-\mathrm{const}\cdot (l/{\ell }_x{)}^{1+{\nu }_{\perp }}}.$$

$${\ell }_y\sim {\ell }_x\sim g^{-\alpha }.$$

$$Y_n=\sum _{m=n_0}^n{\beta }_m,$$

$$\overline{{\Delta }_x}=-gx$$

$$\overline{Y_{2x}}=-\frac{g}{2}{x}^2+O(|gx|)+\mathrm{const}.$$

$$P(x_f)\sim e^{-const\cdot g^2{x}_f^3}.$$

$${\ell }_x\sim g^{-2/3}.$$

$$\overline{{\Delta }_x}=-g(x-x_c),$$

$$\mathrm{Prob}(l>0)\sim g^{\kappa },$$

$$\mathrm{Prob}(l=0)=1-O(g^{\kappa }).$$

$$P(l,g)=g^{2\kappa }\tilde{P}(l{g}^{\kappa }),$$

$$\kappa =\alpha .$$

$$n(g)\sim \mathrm{Prob}(\mathrm{l}>0)/g\sim g^{\alpha -1}$$

$$P_{>}(l)=e^{-C{g}^{{\nu }_{\perp }}{l}^{1+{\nu }_{\perp }}},$$

$$\rho (z)\approx e^{-z-e^{-z}},$$

$$l_g\sim g^{-\alpha }(\ln n_y{)}^{1-\alpha }\sim g^{-\alpha }|\ln g{|}^{1-\alpha }.$$

$$p_{n_y}(l_g)=n_y(n_y-1){\int }_0^{\infty }[P_{<}(x){]}^{n_y-2}\rho (x+l_g)\rho (x)dx,$$

$$p_{n_y}(l_g)\sim \exp \{-C^{\prime }[g\ln (n_y/2){]}^{\alpha }{l}_g\},$$

$$l_g\sim (g|\ln g|{)}^{-\alpha }.$$

$$P(l)\sim e^{-l/{\xi }_i},$$

$$\ln P(l)\sim -{\int }^l{l}^{\prime }/\xi (l^{\prime })d{l}^{\prime },$$

$$P(x_f)=P(-x_f),$$

Full text

I Introduction

In ecology much effort has recently been devoted to the study of population dynamics in the presence of an environmental gradient, which means that the environmental conditions that affect the reproduction rate or lifetime vary smoothly in space, along a given direction lennon ; holt ; hklmt ; oikos . This situation is typically realized for a plant species or vegetation living on a hillside, for which the environmental conditions become more and more favorable or unfavorable with increasing altitude. As a consequence, the density decreases in the unfavorable direction, and the connected patch of colonized area becomes, at sufficiently low densities, fragmented, i.e. it is typically composed of distinct islands. Moving further in the unfavorable direction the local density vanishes. The boundary (also known as the range margin) separating the connected and fragmented zones is the subject of intensive investigations oborny ; gastner , which is also motivated by the phenomenon of global warming. A global change of environmental conditions results in the a shift of the range margin, and this may help to predict the response of a species to climate change or, conversely, it could be used to monitor climate change, see Refs. mustin ; oborny and references therein.

Theoretical modeling is based on spatially explicit metapopulation models, the simplest case of which is the contact process cp . Here, habitat patches are represented by a two-dimensional grid, the sites of which can be in two states: either colonized or empty. Colonized sites stochastically colonize their neighboring sites, or go extinct with certain rates. The spatially homogeneous variant of the model exhibits a continuous phase transition separating a fluctuating active phase from the absorbing one as the control parameter, chosen as the relative rate of colonization vs. extinction, is varied liggett ; md ; odor ; hhl . In the presence of a weak gradient, i.e. a spatially slowly varying local control parameter, the bulk phase transition of the homogeneous model is transformed into a spatially explicit transition, with a vanishing local density below a critical coordinate (altitude, see Fig. 1). Due to the close relation with the critical behavior of the homogeneous model, the properties of the gradient model near the critical coordinate obey a scaling theory, which is characterized by critical exponents that are combinations of standard critical exponents of the homogeneous model sapoval ; oborny .

In addition to a gradient, in real systems there are also fine-scale heterogeneities in the environmental conditions, stemming from the depth of soil or the local surface relief, etc. Assuming these factors vary slowly in time compared to the dynamics of the model, these can be incorporated as quenched random reproduction or extinction rates. For the standard, gradient-free contact process, such a quenched disorder is known to alter the critical behavior, giving rise to singularities also in an extended region around the critical point noest . According to a real-space renormalization method, also known as strong-disorder renormalization group (SDRG) im , the dynamical scaling of the model in one hiv and two dimensions kovacs is ultra-slow moreira and static exponents differ from those of the clean system, related to the infinite-disorder fixed point of the transformation. Whether this type of behavior is valid for any weak randomness requires further investigations. So far, all critical exponents are found to be universal in the sense that they are independent of the specific distribution of random parameters, confirmed also by large scale simulations, even for relatively weak randomness vd ; vfm ; vojta_rev . In addition to the altered dynamical scaling and critical exponents, the quenched model has a number of qualitatively different hallmarks from the clean model. A key feature is that the disorder polarizes the density, i.e. typically anchors the activity to certain sites where the local density (i.e. occupation probability) is $O(1)$ while on the rest of sites it is negligibly small. In this work, we will consider the disordered contact process (DCP) with a gradient in the control parameter, and will use the aforementioned renormalization group method to identify the set of sites with $O(1)$ local density, called as colonized cluster (CC). We then study the width of the front of the CC in one and two spatial dimensions in the stationary state. In addition, mimicking climate change, we will study the evolution of the boundary as the global control parameter is slowly (quasistatically) tuned. According to the results, the boundary positions do not shift smoothly but exhibit significant jumps separated by extended quiescent periods. This phenomenon is reminiscent of punctuated equilibrium in evolution biology eldredge . Therefore, in monitoring climate change, no response could mean a quiescent period preceding a large jump. The distribution of jump lengths is found to have a scaling property in terms of the gradient, which is conjectured to involve the static correlation-length exponent of the model.

The rest of the paper is organized as follows. The model is introduced in section II, and its SDRG treatment is reviewed in section III. In section IV, general scaling considerations about the width of the front are presented. In section V, the distribution and the shift of the front position under a global change of the control parameter are studied in the one-dimensional model within the SDRG approach through a mapping to time-dependent random walks, while, in section VI, we address the same questions in the two-dimensional model by applying the SDRG method numerically. Some details of the calculations are presented in the Appendix. Finally, the results are summarized in section VII.

II The model

We consider a variant of the contact process cp ; liggett . The model is defined on lattice sites which can be in two states, either occupied (colonized) or empty (uncolonized). The dynamics is given by a continuous-time Markov process with two kinds of independent transitions. First, occupied sites can colonize the neighboring (empty) sites, second, occupied sites go spontaneously extinct. We study the simultaneous presence of quenched disorder in the transition rates and a constant gradient in one direction, leading to an average local control parameter that varies linearly in space. Due to the universality of the model around criticality, only the local ratio of the transition rates matters, irrespectively from the specific details of how the disorder and gradient are implemented. Therefore, as a convenient choice, we implement the gradient in the extinction rates and assign the quenched disorder to the colonization rates. The extinction rates vary linearly with the $x$ coordinate as

[TABLE]

where $g$ is the gradient strength, $p$ is a global control parameter, and $x_{c}$ denotes the critical coordinate. The colonization rate between sites $i$ and $j$, $\lambda_{ij}\in(0,1)$, is an independent, identically distributed random variable. Due to the gradient, the average local control parameter, which can be defined as usual for the disordered contact process in terms of the average logarithmic parameters hiv as

[TABLE]

varies with the $x$ coordinate. Here, and in the followings, the overbar denotes an average over quenched disorder. As a result, the mean local density $\overline{\rho}(x)$ also varies with $x$, and, above a critical coordinate $x_{c}$ where $\overline{\Delta}(x)$ hits the critical control parameter of the gradient-free model ($g=0$), $\overline{\Delta}(x_{c})=\ln p-\overline{\ln\lambda}\equiv\Delta_{c}$, the mean density is non-zero, whereas well below $x_{c}$ it tends rapidly to zero. Altogether, in the vicinity of the critical coordinate, $g|x-x_{c}|\ll 1$, the average control parameter varies linearly in leading order:

[TABLE]

The role of the global control parameter $p$ is merely to tune the critical coordinate $x_{c}$ to a prescribed position (typically to the middle of the sample in numerical calculations). We are interested in the properties of the boundary region for small gradients (i.e. small slopes of the hillside), especially in the scaling behavior in the limit $g\to 0$.

In the subsequent analytical calculations, we make restrictions neither on the form of disorder distributions nor on the way how quenched disorder and gradient is particularly realized in the transition rates. The only requirement is that the local control parameter is an independent random variable and its average varies in leading order linearly near $x_{c}$. Note that the latter requirement is generally fulfilled if the average control parameter varies smoothly with the coordinate $x$.

In the numerical calculations, we will consider finite systems with the coordinate $x$ in the range $[1,L]$, and similarly for the coordinate $y$ in two dimensions. For the latter, we consider periodic boundary condition in the $y$ direction. We use the extinction rate profile given Eq. (1), and fix the gradient to be inversely proportional to the system size: $g=a/L$. Due to this, the gradient is controlled by the system size, and we can apply standard finite-size scaling techniques. With the choice $a=2$, the extinction rate at the edge $x=L$ of the lattice is zero, $\mu(L)=0$, which is also known as active wall boundary condition hhl . This choice ensures for the finite systems to have a true non-trivial steady state.

III Determining the colonized cluster by renormalization

A basic question in ecology is to precisely define the range margin of a species or vegetation type exposed to an environmental gradient. On the modeling side, one usually samples a random configuration in the stationary state, where some of the sites are occupied while the rest is not (Fig. 1). Then, treating the sample as a percolation problem, the hull of the percolating cluster can serve as the boundary between the two zones gastner . The boundary defined this way is a stochastic object which fluctuates in time along with the configuration of the system. It is important to note that the mean $x$ coordinate of the hull, $x_{h}$ is different from $x_{c}$, as it is determined by the percolation threshold of typical configurations, which is different from the critical point gastner .

In presence of quenched disorder, the local densities, $\rho_{i}$, are site dependent. This enables an alternative classification of the sites as those having a local density larger than an arbitrary threshold $\rho_{c}$ and those having $\rho_{i}<\rho_{c}$. The so defined colonized sites have two important differences from the connected zone delineated by the hull. First, these zones do not fluctuate stochastically and are determined by the underlying random environment (and $\rho_{c}$). Second, since the local densities do not vary monotonically with the coordinates, the active zone is no longer necessarily connected in the percolation sense.

There is a well-established coarse-graining procedure of the disordered contact process to find the sites which are colonized with a high probability, without requiring to set an arbitrary threshold $\rho_{c}$. The method is known as strong-disorder renormalization group in the physics literature im and was originally designed for studying low-energy properties of disordered quantum spin chains. Its applicability and power to describe the long-time behavior of the disordered contact process was recognized later hiv . According to this method, the critical behavior of the model in one and two dimensions, at least for not too weak disorder, is governed by a so called infinite-disorder fixed point of the renormalization transformation, at which the procedure becomes asymptotically exact im .

The renormalization procedure consists of two kinds of steps, which are applied sequentially, depending on whether the momentarily largest parameter of the model is a colonization or an extinction rate. When the maximal rate is a colonization rate ($\lambda_{ij}$), the two sites form a highly correlated cluster, with an effective extinction rate obtained by perturbation calculation, since all other rates are smaller:

[TABLE]

On the other side, when the maximal rate is an extinction rate $\mu_{i}$, the site is most of the time extinct, therefore we can remove it from the system, apart from weak induced interactions between its neighboring sites, obtained perturbatively as

[TABLE]

At the critical point, the logarithmic rates increase in magnitude without limits during the renormalization procedure, therefore the term $\ln 2$ in Eq. (4) will not influence the asymptotic properties and can be safely omitted sm . At criticality, the distribution of the logarithmic rates also gets broader, increasing the accuracy of the perturbative approach as the process is performed. Beyond the one dimensional geometry, the above steps can generate a weak colonization rate between existing sites that have been already connected by a colonization rate. To resolve this problem we follow the standard maximum rule, according to which the smaller rate is omitted. An advantage of this approximation is that it is self-consistent and enables a very efficient numerical implementation kovacs . Then, in both steps the generated effective rates are smaller than the eliminated ones.

Applying the SDRG to a finite system, some of the sites will be eliminated during the procedure (i.e. rendered inactive), and some survive as constituents of the last remaining cluster the procedure ends up with. The sites of this cluster are precisely those which are colonized in the stationary state with an $O(1)$ probability, so they constitute the colonized cluster, while the local density on other (eliminated) sites becomes negligibly small. The active wall boundary condition $\mu(L)=0$, obtained with $a=2$, means that the CC contains the sites at $x=L$, as they are never eliminated.

The validity of the SDRG method in describing the criticality of the disordered contact process has been confirmed by Monte Carlo simulations in several works, see e.g. vd ; vfm ; im . In addition, a direct comparison of the density profile in the stationary state shows a good overlap with the CC of the SDRG in individual samples kj .

IV Scaling of the width of the front

In particle systems subject to a gradient, such as diffusion sapoval or contact process gastner , the hull of the percolating cluster fluctuates in time, and the amplitude of deviations in the gradient direction increases as the gradient decreases. Concerning the boundary of the colonized cluster in the disordered contact process, it does not have temporal fluctuations. Yet, it still has sample-to sample fluctuations, and in two dimensions, even in a fixed sample, the $x$ coordinate of the boundary varies along the $y$ direction. Therefore, we can define the width of the boundary even in quenched disordered models.

The scaling of the width of the diffusion frontier with a gradient was derived in Ref. sapoval . The argument, which will be briefly recapitulated, is, however, quite general and applies also to the boundary of the CC in the disordered contact process as follows. Away from the critical coordinate, in the unfavorable zone, the activity is typically restricted to distinct “islands”, while in the favorable zone, empty patches (“lakes”) appear in the occupied background (Fig. 1). The characteristic linear size of islands and lakes is given by the local correlation length, $\xi[\overline{\Delta}(x)]$, which varies with $x$ through the local control parameter:

[TABLE]

where $\nu_{\perp}$ is the correlation-length exponent of the model. The ruggedness of the population “shoreline” can be viewed as a result of large islands in the unfavorable zone reaching and joining to the occupied background and, similarly, large lakes in the favorable zone getting connected with the “sea”. The farthest position from $x_{c}$ at which this can happen is given by the condition that the characteristic size of islands and lakes at that position is comparable with the distance $|x-x_{c}|$. This gives the characteristic width of the margin and also the maximal distance, $\ell_{x}$ as $\ell_{x}\sim(g\ell_{x})^{-\nu_{\perp}}$, yielding

[TABLE]

where the exponent $\alpha$ is given by

[TABLE]

The above reasoning holds also for the disordered contact process, for which the characteristic size of islands and lakes is determined by the fluctuations of disorder through the correlation-length exponent of the disordered model. The numerical values of the exponents $\nu_{\perp}$ and $\alpha$ for the clean and disordered CP in one and two dimensions are summarized in Table 1.

The characteristic length $\ell_{x}$ gives the width of the fluctuations of the front in the $x$ direction, scaling as in Eq. (7). Besides, it can also be interpreted as a correlation length characterizing the spatial correlations in the gradient direction between $x_{c}$ and another $x$, across a medium where the local control parameter is continuously changing in space. This correlation function or the distribution of the $x$ coordinate of the front, as it is derived in section A, has a compressed exponential tail:

[TABLE]

The gradient thus gives rise to a finite correlation length $\ell_{x}$ in the gradient direction. We mention that this will also induce a finite correlation length $\ell_{y}$ in the $y$ direction at $x=x_{c}$, which, due to the intrinsic isotropy of the (gradient-free) model, is expected to scale with $g$ in the same way as $\ell_{x}$:

[TABLE]

V Disordered contact process in one dimension

V.1 Width of the boundary

We will show that, for the DCP in one dimension, the general result in Eq. (7) obtained by scaling considerations can be confirmed by a direct calculation within the frames of the SDRG. In this case, we will consider the coordinate $x_{f}$ of the outermost site of the CC, and define the width of the front by the standard deviation of $x_{f}$, $\ell_{x}=\sqrt{\overline{x_{f}^{2}}-\overline{x_{f}}^{2}}$, over different samples of disorder realizations. Again, the front position $x_{f}$ in a given realization of the random environment is sharp per definitionem, and we consider here its variation with the underlying random environment.

We start with a zig-zag path mapping of the random environment, see e. g. Ref. jkri for a similar model. Introducing the logarithmic variables $\beta_{2x-1}=\ln\mu_{x}$ and $\beta_{2x}=-\ln\lambda_{x,x+1}$, we define a sequence as

[TABLE]

which is a sum of random terms with alternating signs (Fig. 2). Here, we implicitly assumed that all rates lie in the interval $(0,1)$. Restricting ourselves to even values of the indices, $n=2x$, the alternation is eliminated, and $Y_{2x}$ can be regarded as a special, one-dimensional random walk in discrete time $2x$. Note that the space variable $x$ of the original problem appears here as the time variable of the random walk. The action of the SDRG transformation is particularly simple in terms of the sequence $Y_{n}$. The shortest in magnitude step of the walk $|\beta_{m}|=\min\{|\beta_{x}|\}$ is picked and eliminated together with its two neighbors and replaced by a single step $\tilde{\beta}=\beta_{m-1}+\beta_{m}+\beta_{m+1}$ of longer time. This corresponds to the coarse-graining of the path $Y_{n}$ as illustrated in Fig. 2a.

The local order parameter is given by $\Delta_{x}=\ln(\mu_{x}/\lambda_{x,x+1})=\beta_{2x-1}+\beta_{2x}$. Within our SDRG approach of the problem, the critical point is given by $\overline{\Delta}=0$, so the critical coordinate (or time) is simply $x_{c}=0$. Thus, in line with Eq. (3) we have

[TABLE]

at least for $|gx|\ll 1$. We assume that the left boundary $x_{0}=(n_{0}+1)/2$ of the system, as well as the right one is far from the origin compared to the yet unspecified width $\ell_{x}$ of the critical zone, so that the system can be regarded as practically infinite. At even values of the indices, $n=2x$, we have $Y_{2x}=\sum_{m=n_{0}}^{2x}\beta_{m}=\sum_{m=x_{0}}^{x}\Delta_{m}$. This means that the random walk $Y_{2x}$ has a time-dependent bias given in Eq. (12), which drives the walker to the positive (negative) direction for $x<0$ ($x>0$) (Fig. 2). The average path as the function of the discrete time is thus parabolic close to the critical time $x=0$ (Fig. 14):

[TABLE]

Now we make the following statement. The time $2m$ at which the path reaches its maximum, $Y_{2m}=\max_{n}\{Y_{n}\}$, determines the position of the front as $x_{f}=m+1$. The proof of this statement can be found in B. According to the approximate calculations based on properties of random walks in B, $P(x_{f})$ scales as

[TABLE]

This means that the width $\ell_{x}$ of the distribution of $x_{f}$ scales with $g$ according to

[TABLE]

These results are consistent with the general scaling results in Eqs. (7) and (9), using the correlation-length exponent $\nu_{\perp}=2$ of the DCP. Note that the coordinate $x$ of the original problem corresponds to the time in the equivalent random walk picture, therefore $\nu_{\perp}$ is given by the temporal correlation-length exponent of random walk, see Eq. (35).

At the characteristic distance $\ell_{x}$ from the critical coordinate, the local control parameter is $|\overline{\Delta_{\ell_{x}}}|\sim g\ell_{x}\sim g^{1/3}$, which tends to zero in the limit $g\to 0$. This further justifies the approximation made in the SDRG rules of the contact process, i.e. the neglection of the term $\ln 2$ in Eq. (4).

To test the assumptions behind Eq. (14), we performed numerical simulations of random walks with a time-dependent bias. We have drawn $\Delta_{n}$ from uniform distributions of unit width with an offset of the support increasing linearly in time, so that the average control parameter changes in time as $\overline{\Delta_{x}}=gx$. Distributions of $x_{f}$ for different values of the gradient $g$ are shown in Fig. 3. The data show a good collapse in terms of the scaling variable $x_{f}g^{2/3}$. Furthermore, as it is demonstrated in Fig. 3b, the tail of the distribution is in accordance with the compressed exponential form in Eq. (14).

V.2 Evolution of the front under a global change of the environment

Motivated by the potential role of range margins in monitoring climate change mentioned in the Introduction, in the following we will examine the question how the front $x_{f}$ shifts in a fixed (large) random environment when the control parameter is globally changed. In particular, we examine the effects of quenched disorder in such a scenario. To start, we assume a constant average gradient, as before,

[TABLE]

but in a random environment with fixed colonization rates. Then, we tune the critical coordinate $x_{c}$ by shifting the extinction rate profile in Eq. (1). The front position $x_{f}$ is calculated for a range of $x_{c}$, in the corresponding stationary states for each value of $x_{c}$ by the SDRG method, and we are interested in the change of $x_{f}$ under a unit increase in $x_{c}$. This approach corresponds to the quasistatic change of the global environment, i.e. sufficiently slow that the system is able to relax to the instantaneous stationary state. The validity of this approach is quantitatively discussed in Appendix C.

First, to gain a general impression of the change of the density profile under shifting $x_{c}$, we performed numerical simulations of the DCP in samples of size $L=200$, when $x_{c}$ is swept from $x_{c}=51$ to $x_{c}=150$ in unit steps, leaving sufficient time for the system to relax and for a measurement of the density profile for each value of $x_{c}$. Density profiles obtained in three different samples are shown in Fig. 4.

In comparison, the SDRG approach classifies the lattice sites to be either practically empty $\rho_{i}\approx 0$ or fully occupied $\rho_{i}\approx 1$, with the binary assumption being valid for asymptotically large scales, while the local densities are blurred on a microscopic scale in the simulations. Yet, with increasing system size, the profiles, when viewed on a macroscopic scale, are expected to become more and more contrasted and to approach to the outcome of the SDRG method. In spite of these difficulties in small systems, we can clearly observe that for increasing $x_{c}$, there are idle periods in which the front of the high-density region practically does not shift (just a weak overall decrease of the density occurs), and then abruptly a cluster of sites fades away, resulting in a considerable advance of the front.

Let us continue now with the SDRG treatment of the model. We have seen in the previous section that, for a given value of $x_{c}$, $x_{f}$ is determined by the maximum point of the path $Y_{x}$, or, in other words, by the coordinate at which the path $Y_{x}$ is touched from above by a straight, horizontal line. The path has a parabolic trend in leading order around $x_{c}$, $\overline{Y_{x}}=-\frac{g}{2}(x-x_{c})^{2}$. To get a clearer picture of the problem, this trend can be transferred from the path to the touching line. Considering an unbiased random walk path, $\tilde{Y}_{x}\equiv Y_{x}-\overline{Y_{x}}$, for which $\overline{\tilde{Y}_{x}}=0$ for all $x$, $x_{f}$ is equivalently determined by the touching point of $\tilde{Y}_{x}$ with a parabola $Y(x)=\frac{g}{2}(x-x_{c})^{2}$. Going further, for small $g$, the parabola is not much different from a circle of diameter $d=1/g$ in the vicinity of its bottom point, so we may think of the following picture, illustrated in Fig. 5.

We have a wheel of diameter $d=1/g$, the horizontal coordinate of its center is $x_{c}$, and it is rolling on a rough (scale-free) ’surface’ formed by the unbiased random walk path $\tilde{Y}_{x}$. The horizontal coordinate of the contact point gives $x_{f}$. It is now intuitive that the contact point will change intermittently, i.e. for most steps of $x_{c}$ it remains pinned, while there are certain rare steps at which it suddenly makes a large jump. This is also supported by Fig. 6, in which $x_{f}$ is plotted as a function of $x_{c}$ for different values of the gradient.

To obtain quantitative results, $x_{f}$ was calculated numerically in long samples under unit changes of $x_{c}$ ($10^{7}$ steps for each value of $g$). We denote the change in $x_{f}$ by $l$. According to the results, the probability that a non-zero jump occurs scales as

[TABLE]

while the probability that the front remains pinned under a unit change of $x_{c}$ is

[TABLE]

This means that the shift of the front is a rare event for small gradients, but once it happens, the front makes a long jump of typical length $O(g^{-\kappa})$. Accordingly, the distribution of jump lengths has the scaling property

[TABLE]

see Fig. 7. The exponent $\kappa$ is found by the scaling collapse to be $\sim 2/3$.

Obviously, the typical non-zero jump lengths can not exceed the width of the critical zone, $\ell_{x}$, therefore $\kappa\leq 2/3$. Here we conjecture that the jump lengths are in the order of the correlation length $\ell_{x}$, and thus

[TABLE]

So far we considered an elementary, unit shift in $x_{c}$, corresponding to a change of $g$ in the order parameter $\overline{\Delta_{x}}$, see Eq. (3). If the order parameter $\overline{\Delta_{x}}$ is changed by $O(1)$, the front moves forward by $\Delta x_{f}\sim g^{-1}$ and, on average, this $O(g^{-1})$ front displacement is realized in distinct jumps,

[TABLE]

in number. This is in accordance with numerical results shown in Fig. 8, where the average number of non-zero jumps of the front, ${\rm Prob(l>0)}/g$, during a change $\Delta x_{c}=1/g$ is plotted against $g^{\alpha-1}$.

VI Disordered contact process in two dimensions

VI.1 Width of the front

In this section, we consider the more realistic, two-dimensional, DCP in the presence of a gradient. As the colonized cluster is generally not connected, it is not possible to delineate it by a single connected boundary. Of course, one could still consider the percolation hull of the colonized sites, in the same way as for the clean model. However, by this definition, some “outposts” which are physically part of the CC but not connected with it in the percolation sense, would fall in the outer side of the hull. It is thus an improper choice for the characterization of the disordered problem. Here, we study the set of $x$ coordinates of outermost sites for each value of the coordinate $y$.

We performed numerical SDRG calculations in samples of size $L\times L$, with $g=1/L$, and collected the front coordinates $x_{f}$ at all transversal coordinate $y$. The distributions are shown in Fig. 9.

As can be seen, using the correlation-length exponent of the two-dimensional DCP quoted in Table 1, the distributions follow the general scaling law in Eq. (7). Furthermore, the tails of the distribution are compatible with the compressed exponential form in Eq. (9). In contrast to the one-dimensional model, the distributions are not symmetric around $x_{c}$ and the most likely position is in the active domain ($x>x_{c}$).

In a system of finite transversal size $L_{y}$, one can also consider the $x$ coordinate of the outermost site of the CC in the whole sample, $x_{g}=\min_{1\leq y\leq L_{y}}\{x_{f}(y)\}$, corresponding to the global boundary. Then, the approximate distribution of $x_{g}$ can be obtained by extreme-value analysis. For this, we have to take into account that the front coordinates $\{x_{f}(y)\}$ in a given sample are correlated on a scale $\ell_{y}$, see Eq. (10). The effectively independent number of the data is therefore $n_{y}\sim L_{y}/\ell_{y}$. Starting with a compressed exponential parent distribution of $l=|x_{f}-x_{c}|$,

[TABLE]

it is straightforward to show that the maximal, global deviation $l_{g}=\max_{1\leq n\leq n_{y}}\{l_{n}\}$ follows, for $Cg^{\nu_{\perp}}l_{g}^{1+\nu_{\perp}}\gg 1$, the Gumbel distribution

[TABLE]

given in terms of the reduced variable $z=Cg^{\nu_{\perp}}l_{g}^{1+\nu_{\perp}}-\ln n_{y}$. As $z$ is an $O(1)$ random variable, for $\ln n_{y}\gg 1$ we obtain that the typical value of the global deviation scales as

[TABLE]

Here, the second relation is valid for a sample of size $L\times L$, for which $n_{y}\sim g^{\alpha-1}$, see Eq. (10).

VI.2 Evolution of the front

We now turn to the question of how the front advances under a global change of the environment. We follow the same protocol as for the one-dimensional model, i.e. the critical coordinate $x_{c}$ is changed in unit steps, and the set of local front coordinates $\{x_{f}(y)\}$ is determined by the SDRG method for each $x_{c}$. Typical colonized clusters at subsequent steps of $x_{c}$ are shown in Fig. 10.

We monitored the local shift $l$ of $x_{f}(y)$ for fixed transversal coordinates $y$, as well as the global shift $l_{g}$ of the extremal coordinate $x_{g}$. Note that in the latter case the $y$ coordinate of the globally outermost site is not fixed. The local shift distributions calculated for each $y$ are shown in Fig. 11 in a few thousands of samples for each $g$.

Similarly to the one-dimensional model, the local shifts follow the intermittent dynamics described by Eqs. (18-19), and a good scaling collapse of the distributions is achieved by $\kappa=0.554$, in accordance with the generalization of our conjecture $\kappa=\alpha$ formulated in one dimension. We can also study the number of jumps to achieve a total displacement $\Delta x_{f}\sim g^{-1}$ of the local front position, corresponding to an $O(1)$ change in the control paremeter. Numerical results for the average number of non-zero jumps of the local front coordinate shown in Fig. 8 are in accordance with the expectation $O(g^{\alpha-1})$, see Eq. (21).

Concerning the global jumps $l_{g}$, we make the following considerations. When the outermost site disappears (together with a bunch of sites in the correlation volume $\ell_{x}\times\ell_{y}$) under the change of $x_{c}$, the next outermost site will be typically in a different correlation volume of the front region. The length $l_{g}$ of a global jump (once it occurs) is thus typically given by the difference of the smallest and the second smallest front coordinate $x_{f}(y)$ in the sample. Assuming again independent variables, $n_{y}$ in number, one obtains in straightforward calculations that the distribution of the “gap” $l_{g}$ between the first and second extrema has the distribution

[TABLE]

where $P_{<}(x)=1-P_{>}(x)$ is the parent distribution given in Eq. (22) and $\rho(x)$ is the corresponding probability density. Here, we assumed for the sake of simplicity that the parent distribution is continuous. Although the integral cannot be evaluated with the distribution in Eq. (22), the scaling of $l_{g}$ can be determined by the saddle-point approximation. It is straightforward to show that, under the condition $\ln n_{y}\gg 1$, the saddle point of the integrand is at $x^{*}\approx\left[\frac{g^{-\nu_{\perp}}}{C}\ln(n_{y}/2)\right]^{1-\alpha}$, and the tail of the distribution behaves as

[TABLE]

up to exponential precision. In square samples, for which $n_{y}\sim g^{\alpha-1}$ we thus expect the rare, non-zero global jumps to scale as

[TABLE]

Accordingly, the global shift $l_{g}$ obeys similar scaling laws as in Eqs. (18-19) with the difference that $g$ is replaced by $g|\ln g|$. This is supported by the numerical distributions shown in Fig. 12, with a good scaling collapse.

The average number of global jumps under an $O(1)$ change of the global control parameter is expected to scale as $n(g)\sim{\rm Prob(l>0)}/g\sim(g|\ln g|)^{\alpha}/g$, in agreement with numerical data presented in Fig. 13.

VII Conclusion

Motivated by the problem of range margin in an environmental gradient, studied intensively in ecology, we considered in this work the effects of quenched disorder by studying the disordered contact process with a linear spatial trend in the local control parameter. We applied the SDRG approach to the model, enabling us to construct the colonized cluster, i.e. the set of sites occupied with a high probability. The CC is a disconnected set, and we identified its front as the outermost constituent of the CC.

In one dimension, we have shown that the problem of finding the front position is equivalent to an extremal problem of a random walk with a time-dependent bias. We confirmed by explicit calculations that the distribution of the front position for small gradients is in agreement with the general scaling considerations, which predict a compressed exponential tail of the distribution, and in which the model-dependence enters through the correlation-length exponent. Concerning the shift of the front under a unit change of the critical coordinate, we demonstrated that the front advances intermittently, meaning that it remains unchanged in almost all steps, while it makes long jumps in rare steps, the fraction of which vanishes with decreasing gradient. These characteristics of the intermittent motion of the front are found to follow universal scaling laws, conjectured to be related to the correlation-length exponent of the model.

In two dimensions, we studied both the coordinates of outermost sites of the CC for a fixed transversal coordinate, as well as its globally extremal value, by applying the SDRG method numerically. The distribution of local front coordinates was found to follow the general scaling law. Implementing the same protocol as in one dimension, the corresponding local and global shifts were measured. These are found to show the same intermittent motion as in the one-dimensional model, and to obey the same scaling laws (including a logarithmic correction for the global shift) with the corresponding correlation-length exponent.

It is worth comparing the motion of the boundary under the change of control parameter in two dimensions to the thoroughly studied problem of driven interfaces in random media fisher ; giamarchi . The basic ingredients of that model at zero temperature, written as an overdamped dynamical equation, are the elastic interactions between the constituents of the interface, which tend to keep the interface flat, and the static random pinning forces of the medium. As a response to an external force, provided it exceeds a critical depinning force, the interface performs a persistent but jerky motion, of which the intermittent motion of the boundary in our model is reminiscent. Beside the similar appearances, there are, however, several differences between the two models. A formal one is that the description of the dynamics of the boundary in our model by an autonomous dynamical equation in terms of the degrees of freedom of the boundary is unresolved. Furthermore, the boundary, as it is defined, is not necessarily a connected object. Finally, the mean velocity of the boundary (or the change rate of the control parameter) is an external parameter in our model, as opposed to driven interfaces where it is a response to the applied force.

It would be interesting to examine whether quenched disorder has observable effects in real ecological systems. As the difference between the characteristics of clean and quenched disordered systems is most prominent at criticality, an ideal testing ground for disorder effects would be provided by ecological systems subject to an environmental gradient, since here a critical region naturally appears without any fine-tuning of a control parameter. In particular, it would be interesting to validate our key result, the intermittent evolution of the front under global climate change by empirical data. Our results indicate that climate change effects might stay concealed for an extended period of time, leading to sudden large changes in the observations.

Acknowledgements.

We thank B. Oborny, M. Gastner, G. Ódor, and F. Iglói for useful discussions. This work was supported by the National Research, Development and Innovation Office NKFIH under grant No. K128989. I.A.K. was supported by the Domus Hungary Scholarship of the Hungarian Academy of Sciences. This publication was made possible through the support of a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation.

Appendix A Distribution of the front

Let us first assume that the local control parameter is discontinuous at $x_{c}$, while constant otherwise: $\overline{\Delta}(x)=\Delta_{1}<\Delta_{c}$ for $x<x_{c}$ and $\overline{\Delta}(x)=\Delta_{2}>\Delta_{c}$ for $x>x_{c}$. The front is then pinned on average to $x_{c}$ and its penetration into either phase in a distance $l=|x-x_{c}|$ is exponentially improbable:

[TABLE]

where $\xi_{i}$ is the correlation length $\xi_{i}\sim|\Delta_{i}-\Delta_{c}|^{-\nu_{\perp}}$, for $i=1,2$. In a medium with a continuously varying control parameter this can be generalized to

[TABLE]

which, substituting $\xi(l)\sim[\Delta(l)-\Delta_{c}]^{-\nu_{\perp}}\sim(gl)^{-\nu_{\perp}}$, results in a compressed exponential tail of the distribution given in Eq. (9). This result can also be confirmed by formulating lower and upper bounds on $P(l)$ as will be done for a similar function in section B.

Appendix B The front position in one dimension and time-dependent random walks

The proof of the statement formulated in section V.1 is simple. Clearly, the maximum point must be at an even index, since the terms $\beta_{n}$ are alternately positive and negative for even and odd $n$, respectively. Assume that the maximum point is at $2m$. It is obvious from Fig. 2a, that the SDRG procedure always just eliminates some $Y_{n}$, so, until $Y_{2m}$ is not removed it will remain the maximal value at any stage of the procedure. But its removal can never happen, which can be seen as follows. The extinction rate of site $m+1$ is represented by the descending segment of path between $2m$ and $2m+1$, see Fig. 2b. The removal of site $m+1$ could only happen if both of its neighboring ascending segments were higher in magnitude. This, however, cannot be fulfilled for its right neighbor since then $Y_{2m+2}$ would be greater than $Y_{2m}$, contradicting our assumption. Therefore site $m+1$ is never removed during the SDRG procedure. A fusion of site $m+1$ (or the cluster containing site $m+1$) with a cluster on its left would require that the height of the ascending segment on the left of the maximum be smaller in magnitude than those of its neighboring descending segments. However, this cannot be true for its left descending neighbor since then $Y_{2m-2}$ would be greater than $Y_{2m}$. This is again in contradiction with our assumption. Consequently no sites with $n<m+1$ are fused to site $m+1$. We thus conclude that site $m+1$ is the outermost site of the colonized cluster, i.e. $x_{f}=m+1$.

In the sequel, we will make two simplifications. First, as we are interested in the scaling of the front position for small $g$, where $x_{f}$ is typically large, we can write $x_{f}\approx m$. Second, we will consider the series $Y_{n}$ restricted to even indices, $n=2x$, and replace the indices $2x$ with $x$. The problem of determining $x_{f}$ can thus be summarized as follows. We have a random walk with a time-dependent bias given in Eq. (12) which changes sign from positive to negative at time $x=0$, and we are looking for the time $x_{f}$ at which the walker reaches the farthest position in the positive direction.

Clearly, as the bias is symmetric in time, the distribution of $x_{f}$, $P(x_{f})$ is an even function, i.e.

[TABLE]

and the mean value $\overline{x_{f}}=0$. Furthermore, the distribution $P(x_{f})$ is expected to have a maximum at $x_{f}=0$, since, for $x_{f}\neq 0$, the walker moves against the bias during the time interval $|x_{f}|$. Due to the symmetry of $P(x_{f})$, it is enough to deal with the positive part $x_{f}>0$ in the following. Now observe that the probability that the maximum is reached at time $x_{f}$, is equivalent with the probability that the path $Y_{x}$ remains below $Y_{x_{f}}$ in both time directions starting from $x_{f}$, see Fig. 14.

That means, $Y_{x_{f}\pm x}<Y_{x_{f}}$ for all $x=1,2,\cdots$. Thus $P(x_{f})$ is a product of two survival probabilities:

[TABLE]

Here, $P_{\rightarrow}(x_{f})$ is the survival probability of a random walk in a time-dependent bias starting at time $x_{f}$, whereas $P_{\leftarrow}(x_{f})$ is the survival probability of a random walk starting at time $x_{f}$ and moving in reversed time, toward $x=0$ and beyond. The second factor, $P_{\rightarrow}(x_{f})$, describes a situation where the walker is biased away from the absorbing wall it has to survive. For the first factor, $P_{\leftarrow}(x_{f})$, the walker is initially attracted toward the wall from time $x_{f}$ to [math] (keep in mind that time is reversed), and subsequently, for times $x<0$ it is repelled away from the wall. As the more rapidly varying factor is the second one, we assume that the scaling properties of $P(x_{f})$ are determined by those of $P_{\leftarrow}(x_{f})$ up to an exponential precision. $P_{\leftarrow}(x_{f})$ can be decomposed as $P_{\leftarrow}(x_{f})=P_{+}(x_{f})P_{-}(x_{f})$, where $P_{+}(x_{f})$ is the probability that the walk starting at time $x_{f}$ survives up to time [math], whereas $P_{-}(x_{f})$ is the conditional probability that under the survival down to time [math], the walk will survive throughout in the regime $x<0$. Clearly, $P_{-}(x_{f})$ will also depend on $x_{f}$ through the position of the walker at time $x=0$ which depends on $x_{f}$. But, as the walker arrives at time $x=0$ through a time domain in which it is attracted toward the wall, the above dependence is expected to be weak, and remaining within the exponential precision, we may write

[TABLE]

The calculation of $P_{+}(x_{f})$, which is the survival probability of a time-dependent random walk is still a hard problem for which, up to our knowledge, no general results exist. Nevertheless, we can formulate lower and upper bounds of $P_{+}(x_{f})$ and will see that both have the same scaling property.

This can be done as follows. $P_{+}(x_{f})$ describes the survival probability in a time-dependent bias given by the series $gx_{f},g(x_{f}-1),\dots,0$. A lower bound is provided by the survival probability in $x_{f}$ steps in the presence of a constant bias $gx_{f}$, $P_{s}(\Delta=gx_{f},t=x_{f})$. An upper bound can be obtained by a survival probability in a shorter time interval $x_{f}/2$ in the presence of a smaller constant bias $gx_{f}/2$, $P_{s}(\Delta=gx_{f}/2,t=x_{f}/2)$. So we have

[TABLE]

Now we make use of the result that in the presence of a constant bias toward the absorbing wall, the survival probability decays exponentially in time

[TABLE]

where the correlation time diverges according to

[TABLE]

as $\Delta\to 0$. This result was derived in Ref. ri for random walks with a particular (bimodal) distribution of step lengths but, due to universality, it is expected to hold generally. We have thus for the lower bound $P_{s}(\Delta=gx_{f},t=x_{f})\sim e^{-Cg^{2}x_{f}^{3}}$, and the upper bound $P_{s}(\Delta=gx_{f}/2,t=x_{f}/2)\sim e^{-C^{\prime}g^{2}x_{f}^{3}}$, where $C$ and $C^{\prime}=C/8$ are constants. As a consequence, $P_{+}(x_{f})$, as well as $P(x_{f})$ must have the scaling form given in Eq. (14).

Appendix C Validity of the quasistatic approximation

In this section, we derive a condition for the change rate of the control parameter, $v_{\Delta}=\frac{d\overline{\Delta_{x}}(t)}{dt}$, valid for both one and two dimensions, under which the process can be regarded as quasistatic. We have seen that the jumps of the front are rare; if the critical coordinate advances by $\ell_{x}\sim g^{-\alpha}$, only one jump [of length $O(\ell_{x})$] occurs typically. This event is realized by the vanishing of the cluster extending from the old to the new position of front, having a linear size $O(\ell_{x})$. This cluster lies in the critical zone around $x_{c}$, therefore the time scale of this event is well approximated by the lifetime $\tau$ of the critical DCP in a finite sample of linear size $\ell_{x}$. This is known to scale as $\ln\tau\sim\ell_{x}^{\Psi}$, where the barrier exponent $\Psi$ is $1/2$ in one dimension hiv and close to $1/2$ in two dimensions kovacs . A quasistatic change in a weak sense is then achieved if the typical time between subsequent jumps of the front, $t\sim g\ell_{x}/|v_{\Delta}|$, is much longer than the relaxation time: $t\gg\tau$. This leads to the condition for the change rate of the control parameter

[TABLE]

If this condition is fulfilled, the system has sufficient time to relax before the next jump of the front occurs, although, there will be a time lag of $O(\tau)$ in the evolution of the front position, $x_{f}(t)$, compared to the instantaneous stationary value. A more stringent condition of a quasistatic change is obtained by the requirement that the time $t$ under which the critical coordinate makes a unit displacement, $\Delta x_{c}=1$, is much longer than the relaxation time. This imposes the following condition for the change rate of the control parameter:

[TABLE]