Large deviations in models of growing clusters with symmetry-breaking transitions

TL;DR

This paper investigates large deviations of magnetisation in models of growing clusters with symmetry-breaking transitions, revealing unexpected power-law decay of certain fluctuations and their relation to dynamical phase coexistence.

Contribution

It demonstrates that the rate function is zero for intermediate magnetisation values, leading to power-law decay of fluctuations, a novel insight into large deviations in such models.

Findings

Power-law decay of magnetisation fluctuations

Rate function is zero for intermediate magnetisation

Connection to dynamical phase coexistence phenomena

Abstract

We analyse large deviations of the magnetisation in two models of growing clusters. The models have symmetry-breaking transitions, so the typical magnetisation of a growing cluster may be either positive or negative, with equal probability. For large clusters, the magnetisation obeys a large deviation principle. We show that the corresponding rate function is zero for values of the magnetisation that are intermediate between the two steady state values, which means that fluctuations with these values of the magnetisation are much less unlikely than previously thought. We show that their probabilities decay as power laws in the cluster size, instead of the exponential scaling that would be expected from the large deviation principle. We discuss how this observation is related to dynamical phase coexistence phenomena. We also comment on the typical size of magnetisation fluctuations.

Figures

A of Elmatad *et al.* ([2010](#bib.bib50)) and Fig. 4B of Hedges *et al.* ([2009](#bib.bib15)).

A of Elmatad *et al.* ([2010](#bib.bib50)) and Fig. 4B of Hedges *et al.* ([2009](#bib.bib15)).

Equations

$$\overline{F}(\tau )=\frac{1}{\tau }{\int }_0^{\tau }F(x(t))\mathrm{d}t.$$

$$p\left(\overline{F}|\tau \right)\simeq \exp \left[-\tau I\left(\overline{F}\right)\right].$$

$$m_k=\frac{1}{k}\sum _{i=1}^k{s}_i.$$

$$s_{k+1}=\{\begin{array}{ll}+1,& with\; prob(1+\exp (-2J{m}_k){)}^{-1},\\ -1,& with\; prob(1+\exp (2J{m}_k){)}^{-1}.\end{array}$$

$$⟨m_{k+1}{⟩}_{m_k}\simeq m_k+\frac{\tanh (J{m}_k)-m_k}{k}.$$

$$p(m_K)\simeq \exp \left[-KI(m_K)\right]$$

$${\mathcal{I}}_K(m,ϵ)=-\frac{1}{K}\log \mathrm{Prob}\left(|m_K-m|\le ϵ\right).$$

$$\underset{K\to \infty }{\lim }{\mathcal{I}}_K(m,ϵ)=\underset{x\in [m-ϵ,m+ϵ]}{\inf }I(x).$$

$$p(m_K)\simeq \exp [-KI(m_K)].$$

$$(\Delta N_k,\Delta M_k)=\{\begin{array}{ll}(+1,+1),& w/prob(c/z_k)\\ (+1,-1),& w/prob(c/z_k)\\ (-1,+1),& w/prob{\mathrm{e}}^{J{m}_k}(1-m_k)/z_k\\ (-1,-1),& w/prob{\mathrm{e}}^{-J{m}_k}(1+m_k)/z_k\end{array}$$

$$z_k=2(c+\cosh J{m}_k-m\sinh J{m}_k)$$

$$⟨F⟩=⟨F\exp [-\mathcal{A}(\mathbf{m})]{⟩}_{\mathrm{con}}$$

$$\mathcal{A}(\mathbf{m})=\sum _{k=1}^K\log \frac{p_k^{\mathrm{con}}(s_k|m_{k-1})}{p_k(s_k|m_{k-1})}.$$

$${\mathcal{I}}_K(m,ϵ)\le {\mathcal{H}}_K(m,ϵ)$$

$$\begin{gathered} {\mathcal{H}}_K(m,ϵ)=-\frac{1}{K}\log ⟨{\chi }_{ϵ}(m_K-m){⟩}_{\mathrm{con}} \\ +\frac{⟨{\chi }_{ϵ}(m_K-m)\cdot \mathcal{A}[\mathbf{m}]/K{⟩}_{\mathrm{con}}}{⟨{\chi }_{ϵ}(m_K-m){⟩}_{\mathrm{con}}}. \end{gathered}$$

$$s_{k+1}=\{\begin{array}{ll}+1,& with\; prob(1+b)/2\\ -1,& with\; prob(1-b)/2\end{array}$$

$$⟨m_k{⟩}_{\mathrm{con}}=b,⟨(\Delta m_k{)}^2{⟩}_{\mathrm{con}}=(1-b^2)/k.$$

$$\begin{gathered} ⟨\mathcal{A}(\mathbf{m})/K{⟩}_{\mathrm{con}}\simeq \frac{1+b}{2}\log \frac{(1+b)(1+{\mathrm{e}}^{-2Jb})}{2} \\ +\frac{1-b}{2}\log \frac{(1-b)(1+{\mathrm{e}}^{2Jb})}{2} \end{gathered}$$

$$\begin{gathered} {\mathcal{H}}_K(m,ϵ)\simeq \frac{1}{2}\log (1-m^2)+\log \cosh (Jm) \\ +\frac{m}{2}\log \frac{1+m}{1-m}-J{m}^2 \end{gathered}$$

$$⟨\mathcal{A}(\mathbf{m})/K{⟩}_{\mathrm{con}}\simeq \frac{b^2(J-1{)}^2}{2}+\mathcal{O}(b^4)$$

$$\frac{\mathrm{d}m}{\mathrm{d}u}=\tanh Jm-m.$$

$$b^2{k}^{\ast }=a\sqrt{K}$$

$$⟨\mathcal{A}{⟩}_{\mathrm{con}}\approx k^{\ast }\frac{b^2(J-1{)}^2}{2}.$$

$${\mathcal{H}}_K(m,ϵ)\simeq \frac{a(J-1{)}^2}{2\sqrt{K}}.$$

$$I(m)=0,|m|<m_{\mathrm{S}}.$$

$${\rho }_K(m,ϵ)=\frac{1}{2ϵ}\mathrm{Prob}\left(|m_K-m|\le ϵ\right)$$

$$\underset{K\to \infty }{\lim }{K}^{-\beta }\log \mathrm{Prob}\left(|m_K-m|<ϵ\right)=0,$$

$${\mathcal{J}}_K(m,ϵ)=-\frac{\log {\rho }_K(m,ϵ)}{\log K}$$

$$k^{\ast }=\gamma K^{\frac{J-1}{J-(1/2)}}$$

$${\tilde{\rho }}_K(m,ϵ)=\frac{1}{2ϵ}{\mathrm{Prob}}_{\mathrm{con}}\left(|m_K-m|\le ϵ\right)$$

Full text

I Introduction

Large deviation theory is the mathematical framework for analysis of rare events den Hollander (2000). In recent years, this theory has been applied to a wide range of physical systems, in order to study dynamical fluctuations Lebowitz and Spohn (1999); Bodineau and Derrida (2004); Bertini et al. (2015); Mehl et al. (2008); Derrida (2007); Lecomte et al. (2007); Tailleur and Kurchan (2007); Garrahan et al. (2007); Hurtado et al. (2014); Klymko et al. (2018). In non-equilibrium systems, large-deviation theory can be used to analyse fluctuations of currents and of the entropy production Lebowitz and Spohn (1999); Bodineau and Derrida (2004); Mehl et al. (2008); Bertini et al. (2015); Derrida (2007), with implications for fluctuation theorems Lebowitz and Spohn (1999); Harris and Schütz (2007) and thermodynamic uncertainty principles Gingrich et al. (2016). In systems that exhibit metastability, including glasses Garrahan et al. (2007, 2009); Hedges et al. (2009) and biomolecules Weber et al. (2014), large deviation theory can be used to analyse deviations from ergodic behaviour, generating new insights into long-lived metastable states.

To illustrate the simplest case, consider a system with (time-dependent) state $x(t)$, which follows some stochastic dynamics. Select an observable quantity $F(x)$, and construct its time average as

[TABLE]

This quantity is a random variable. A simple statement of ergodicity is that $\overline{F}$ should converge to a corresponding ensemble average, in the limit where $\tau\to\infty$. Large deviation theory provides a framework in which this convergence can be analysed Derrida (2007); Lecomte et al. (2007); Maes and Netocny (2008); Garrahan et al. (2009). Definitions will be given below, here we give an outline of the physical picture. In simple cases, the probability density for $\overline{F}$ scales for $\tau\to\infty$ as

[TABLE]

This type of scaling relationship is called a large deviation principle (LDP) and $I$ is called the rate function den Hollander (2000); Touchette (2009). For physical systems that are finite and ergodic, one expects that the rate function is (strictly) convex, with a single zero when $\overline{F}$ is equal to the ensemble average.

However, there is a wide range of physical systems where this simple picture breaks down. For example, in systems with long-lived metastable states, it is natural for the rate function to develop singularities Garrahan et al. (2007); Jack and Garrahan (2010), which appear as the system size tends to infinity. Systems with long-ranged memory may also behave non-ergodically Harris and Touchette (2009); Harris (2015a, b), in which case the rate function may not be convex

In this work, we consider systems where clusters of particles grow, as a function of time Klymko et al. (2017); Morris and Rogers (2014). The clusters contain two kinds of particles, which can be distinguished either by their spins (up or down), or their colors (red or blue). Using the language of spins, one defines the (extensive) magnetisation, which is the difference between the number of spin-up particles and the number of spin-down particles. We focus here on models that are symmetric under interchange of spin-up and spin-down. In this case, Klymko, Garrahan, and Whitelam Klymko et al. (2017) showed that these systems exhibit a range of interesting behaviour, including symmetry-breaking transitions. That is, depending on the model parameters, the (intensive) magnetisation can converge to zero at large times, or it can converge to a non-zero value Klymko et al. (2017). In the latter case, positive and negative values are equally likely, but typical trajectories of the model involve spontaneous breaking of the up-down symmetry. More generally, these models are also members of a general class of (Polya) urn problems Pemantle (2007); Hill et al. (1980) that have been studied in mathematics; they can also be reformulated as random walks with memory, similar to the elephant random walk Schütz and Trimper (2004) (see also Usatenko and Yampol’skii (2003); Hod and Keshet (2004); Harris (2015b); Baur and Bertoin (2016)).

Recently Klymko et al. (2018), Klymko, Geissler, Garrahan and Whitelam (KGGW) analysed LDPs similar to (2), in these growth models. Since the clusters are always growing with time, the models are not ergodic. This means that one cannot apply standard methods from Derrida (2007); Lecomte et al. (2007); Chétrite and Touchette (2015a) in order to establish LDPs similar to (2). Nevertheless, KGGW showed that an LDP holds at large time, where the observable analogous to $\overline{F}$ is the magnetisation.

For parameters where the symmetry is spontaneously broken, Ref. Klymko et al. (2018) also predicted two unusual properties of the rate function. First, there are some values of the observable for which the rate function is concave (the second derivative is negative). Second, in cases where the steady states spontaneously break the symmetry, they predicted that the rate function should have local minima at unstable fixed points of the dynamics. From (2), it follows that these unstable points should be associated with local maxima in the probability, although these maxima were not observed in Klymko et al. (2018). These predictions were based on a method that provides upper bounds on the rate function Klymko et al. (2018); Whitelam (2018a). This is achieved by modifying the dynamics of the system of interest, in order to make the rare events of interest more likely. A recent preprint Jacobson and Whitelam discusses this method in more detail, and proposes a method for calculation of the difference between the upper bounds and the true rate functions.

In this work, we examine the LDPs of KGGW in more detail. We find that while their upper bounds on the rate function are valid, these bounds are not quantitative estimates of the rate functions themselves. In particular, we show that the second derivative of the rate function may be zero, but it is never negative. (The rate function is not strictly convex, but neither is it concave.) This also explains why the probability distribution of the magnetisation does not have local maxima at unstable fixed points of the dynamics.

Our analysis also reveals additional interesting and unusual behaviour in these systems. In cases where the symmetry is spontaneously broken, we find that the probability distribution of the magnetisation behaves similarly to (2) for some values of the magnetisation, while for some other values it behaves as a power law $p\!\left(\overline{F}|\tau\right)\simeq\tau^{-\alpha}$. We discuss the physical interpretation of this unusual scaling – it is associated with the fact that very large fluctuations can be observed for trajectories that diverge slowly from the unstable fixed point of the dynamics (a similar phenomenon has been discussed before Franchini (2017) in Polya urn models). In the light of these results, we discuss what general insights are available, including the strengths and weaknesses of numerical strategies for analysis of LDPs. In particular, we emphasise that while upper bounds on rate functions are useful, establishing quantitatively accurate bounds is likely to require a detailed understanding of the system of interest Jack and Sollich (2014, 2015, 2010), or a very flexible variational ansatz Bañuls and Garrahan .

In this paper, Sec. II describes the models and the main methods of Klymko et al. (2017, 2018); this includes a model of irreversible cluster growth, and one where the cluster grows reversibly, so that the number of particles can both increase or decrease. Sec. III has results for the irreversible process, and Sec. IV has results for the reversible process. In Sec. V we discuss what general conclusions can be drawn from our results.

II Models and methods

This Section defines the models considered in this work, which were introduced by Klymko, Garrahan and Whitelam Klymko et al. (2017). We also describe the method used by KGGW to obtain bounds on the probabilities of rare events in these models Klymko et al. (2018). The definitions of the models and methods are equivalent to those of Klymko et al. (2018). We also discuss the relationships between those works and other models and methods from the literature.

II.1 Irreversible model of growth

In the irreversible model of growth defined in Klymko et al. (2017), clusters grow irreversibly (see also Morris and Rogers (2014)). KGGW consider trajectories of this model that have a fixed number of events Klymko et al. (2018), which amounts to considering a stochastic process in discrete time. Let $s_{k}=\pm 1$ be the spin of the particle that is added on step $k$. The (intensive) magnetisation just after this step is

[TABLE]

The initial condition is that $s_{1}=\pm 1$, each with probability $1/2$. The dynamical rule is that

[TABLE]

The parameter $J>0$ describes a ferromagnetic interaction among the spins. We consider trajectories of this discrete-time model with $K$ steps in total: this set of trajectories is identical to the constant-event-number ensemble of KGGW Klymko et al. (2018), where the number of events is $K$.

There are several possible ways to analyse this growth process. It can be interpreted as a Markov process for $m_{k}$, in which the transition probabilities depend explicitly on the step $k$. (The process is not stationary.) This model can be mapped onto a (generalised) Polya urn problem Pemantle (2007); Hill et al. (1980) by interpreting the number of up/down spins in the growing cluster as the number of red/black balls in a container (similar to Klymko et al. (2017)). General results for large deviations in urn models have been obtained by Franchini Franchini (2017): the results presented here are consistent with that work. Alternatively, the system can be viewed as a non-Markovian process for $s_{k}$, in which the transition probabilities depend on the history, via $m_{k}$ Whitelam (2018b, a). This formulation is related to the elephant random walk Schütz and Trimper (2004), where $M_{k}=km_{k}$ is viewed as the position of a particle, which obeys a dynamical rule similar to (4) except that the probabilities of $s_{k+1}=\pm 1$ are linear in $m_{k}$. Within this non-Markovian formulation, the irreversible growth model falls in the broad class considered by Harris and Touchette Harris and Touchette (2009); Harris (2015a). They derived several general results for large deviations in models within this class. The results presented here are consistent with their theory.

The large-$k$ behaviour of this model depends on the value of $J$. For large $k$, the magnetisation $m_{k}$ converges to a limiting value. It is convenient to refer to this as a steady state for the model, even though the cluster is constantly growing so the system is not stationary. Given some $m_{k}$, the conditional average of $s_{k+1}$ is $\langle s_{k+1}\rangle_{m_{k}}=\tanh Jm_{k}.$ For large $k$ one finds

[TABLE]

This means that steady state values of the magnetisation must solve $m=\tanh Jm$ Klymko et al. (2017). This equation is familiar from mean-field models of ferromagnets: there is only one solution if $J\leq 1$ but for $J>1$ there are three solutions. In this case, the steady states have $m=\pm m_{\rm S}$ where $m_{\rm S}$ is the spontaneous magnetisation. The point $m=0$ is an unstable fixed point of (5) and is not a steady state of the growth model. This symmetry-breaking transition is similar to transitions in the non-Markovian random walk models of Harris (2015b), see in particular their (so-called) artificial model.

We consider trajectories with a total of $K$ steps, and we discuss an LDP that applies for large $K$. The analogue of the time-average in (1) is $m_{K}$, as defined in (3). The LDP discussed by KGGW is similar to (2): as $K\to\infty$ one has

[TABLE]

where $I$ is the rate function. This is an LDP with speed $K$ den Hollander (2000); Touchette (2009). We assume that $I$ is continuous, which is consistent with our results and those of KGGW. Define

[TABLE]

where $\mathrm{Prob}(\ldots)$ on the right hand side denotes the probability that $m_{K}$ is within $\epsilon$ of some value $m$. From the theory of LDPs den Hollander (2000), one has

[TABLE]

Hence, analysis of ${\cal I}$ provides information about the rate function. In particular, the assumed continuity of $I$ means that $I(m)=\lim_{\epsilon\to 0}\lim_{K\to\infty}{\cal I}_{K}(m,\epsilon)$.

II.2 Reversible model of growth

The reversible model considered by KGGW was introduced in Klymko et al. (2017). It is similar to the irreversible one, except that particles may leave the cluster as well as being added to it. (Note however, this is not a reversible Markov chain in the mathematical sense.) The number of particles in the cluster after step $k$ is $N_{k}$ and its (extensive) magnetisation is $M_{k}$. Also $m_{k}=M_{k}/N_{k}$ is the intensive magnetisation. There are several LDPs that could be considered, including the joint probability distribution for $M_{K}/K$ and and $N_{K}/K$, which was analysed by KGGW Klymko et al. (2018). Here we consider the behaviour of $m_{k}$, for which we expect that

[TABLE]

This is analogous to (6), although the functional form of $I$ will be different.

The dynamical rule for the model may be formulated in terms of increments $\Delta M_{k}=M_{k+1}-M_{k}$ and $\Delta N_{k}=N_{k+1}-N_{k}$. The physical idea Klymko et al. (2017) is that particles are added with rates that do not depend on their spin, but their unbinding rate is suppressed if they are aligned with the magnetisation of the cluster. Specifically

[TABLE]

where

[TABLE]

is a normalisation constant. Similar to the irreversible case, the set of trajectories of this discrete-time model is equivalent to the constant-event-number ensemble of the reversible growth model of Klymko et al. (2018).

Like the irreversible case, this model can also be mapped to generalised Polya urn problem, where particles can be removed from the container, as well as added. The model can also be interpreted as non-Markovian random walk in two dimensions, where the position is $(M_{k},N_{k})$. This would be a two-dimensional generalisation of the elephant random walk Schütz and Trimper (2004).

The steady states of this model are described in Klymko et al. (2017). There are symmetry-breaking transitions, similar to the irreversible model, but with some additional complexity. We give a brief summary of the relevant behaviour in Section IV, which also includes results for this model.

II.3 Bounds on probabilities via optimal control theory

One method for analysing large deviations is to establish upper bounds on the rate function, by considering processes where the rare events of interest become typical. This idea is common in the mathematical theory of large deviations den Hollander (2000); Dupuis and Ellis (1997), its application in physics is reviewed in Chétrite and Touchette (2015b), which also discusses its connection with optimal control theory Fleming (1985). This connection provides us with a useful terminology: we consider control forces which are added to the system, in order to enhance the probability of rare events (for a physical example, see Nemoto et al. (2019)).

In mathematical studies, one usually aims to prove bounds on rate functions. In KGGW Klymko et al. (2018), bounds were evaluated numerically, by direct simulation of the growth model, see also Whitelam (2018a); Jacobson and Whitelam . (This approach may be contrasted with other numerical methods Nemoto and Sasa (2014); Nemoto et al. (2016, 2017); Ray et al. (2018); Ferré and Touchette (2018); Nemoto et al. (2019), which use control forces to aid the computation of rate functions or other large-deviation properties, instead of computing bounds.)

We outline the derivation of the relevant bound, following Klymko et al. (2018); Whitelam (2018a). The method is very general, we focus here on the irreversible growth model. Let $p_{k}(s|m)$ be the probability that the $(k+1)$th particle has spin $s$, given that the magnetisation just after step $k$ is $m$. The probability of a trajectory $\bm{m}=(m_{1},m_{2},\dots,m_{K})$ for the original process is $P(\bm{m})=\prod_{k=1}^{K}p_{k}(s_{k}|m_{k-1})$, with $s_{k}=m_{k+1}-m_{k}$. In the irreversible model, (4) shows that $p_{k}(s|m)=[1+\exp(-2Jms)]^{-1}$, for $k\geq 2$. This is independent of $k$ because the update rule does not depend on $k$, except through the value of $m_{k}$. The initial condition is specified by taking $p_{1}(s|m)=(1/2)$ (for $s=\pm 1$).

Now introduce a controlled system in which the corresponding transition probabilities are $p^{\rm con}_{k}(s_{k+1}|m_{k})$. In Klymko et al. (2018); Whitelam (2018a); Jacobson and Whitelam , this process is called the reference system. The average of some observable quantity $F$ in the controlled process is denoted by $\langle F\rangle_{\rm con}$, and the corresponding average in the original system is $\langle F\rangle$. By considering the probabilities of individual trajectories, one sees that

[TABLE]

where the action $\cal A$ is

[TABLE]

We define an indicator function $\chi_{\epsilon}(x)$ that is equal to unity if $|x|<\epsilon$, and zero otherwise. Then use (7) to write ${\cal I}_{K}(m,\epsilon)=-K^{-1}\log\langle\chi_{\epsilon}(m_{K}-m)\rangle$. Using (12) to express the right hand side in terms of averages with respect to the controlled process, and noting that ${\rm e}^{x}$ is convex, one obtains by Jensen’s inequality Klymko et al. (2018); Whitelam (2018a) that

[TABLE]

with

[TABLE]

The first term on the right hand side of (15) is the log probability of the event that $|m_{K}-m|\leq\epsilon$, in the controlled process. The second term is a conditional average of the action, which is obtained by averaging over the trajectories that realise this event. If this event of interest is typical under the controlled dynamics, it is simple to evaluate ${\cal H}$ and hence to obtain an upper bound on ${\cal I}$. Obtaining accurate bounds (with ${\cal H}\approx{\cal I}$) requires a good choice of the controlled process.

For ergodic Markov processes in the classes considered by Lecomte et al. (2007); Jack and Sollich (2010); Chétrite and Touchette (2015a) (which include irreducible finite-state Markov chains and a large class of stochastic differential equations), the search for suitable controlled processes is somewhat simplified. In these cases, it can be shown Maes and Netocny (2008); Jack and Sollich (2010); Chétrite and Touchette (2015a) that for $K\to\infty$ one may achieve equality in (14) by using a controlled process where the transition probabilities do not depend on the step $k$ (except possibly for $k=1$, which controls the initial condition). However, the situation for growth models is more complicated because achieving equality in (14) may require a controlled process whose transition rates depend explicitly on $k$. The new results that we obtain here are obtained by considering this type of controlled process.

III Results : irreversible model

We analyse the LDP (6) in the irreversible model. KGGW derived bounds on the rate function $I$ that appears in (6). They argued that their bounds are equal to the rate function itself, because the controlled processes (reference processes) that they consider mimic the rare-event behaviour of the model. We show here that their bounds are valid but they are not equal to the rate function. Section III.1 focusses on the case $J>1$ in which the model spontaneously breaks up-down symmetry. We show that the rate function in (6) is zero throughout the range $|m|<m_{\rm S}$ (see also Corollary 5 of Franchini (2017)). Fig. 1 summarises this result. After that, Sec. III.2 discusses the behaviour of $p(m_{K})$ for $m_{K}$ in this range: we show that this probability decays as a power law in $K$. Then Sec. III.3 gives a brief discussion of the behaviour for $|m_{K}|\gtrsim m_{\rm S}$, focussing on the behaviour close to the peaks of the probability distribution.

III.1 LDP with speed $K$

This Section establishes the new bound shown in Fig. 1. Before starting the calculation, we recall the physical interpretation of this rate function. In many LDPs, zeros of the rate function $I(m)$ correspond to typical values of the observable $m$ Touchette (2009). The situation shown in Fig. 1 is different, in that the typical values of $m$ are $\pm m_{\rm S}$ but the rate function is zero throughout the range $|m|\leq m_{\rm S}$. This behaviour occurs when the probability density for $m$ decays to zero less fast than an exponential [for example $p(m_{K})\sim K^{-\alpha}$] so that $p(m_{K})\to 0$ at large $K$, but one still has ${\cal I}_{K}(m,\epsilon)\to 0$, and hence $I(m)=0$ by (8). This behaviour may be somewhat unusual but it is fully consistent with the existence of a large deviation principle with speed $K$ den Hollander (2000). It is also similar to the behaviour of the free energy in thermodynamic systems close to phase coexistence, see Sec. V.

III.1.1 Control forces that are independent of $k$

We first derive the bounds of KGGW, following their method, which uses (14). We use a slightly different controlled system in which the extensive magnetisation $M_{k}=km_{k}$ follows a biased random walk (see also Whitelam (2018a, b)). The dynamical rule for this controlled system is

[TABLE]

where $b$ is a numerical parameter with $|b|\leq 1$. In this case, the mean and variance of $m_{k}$ under the controlled dynamics are

[TABLE]

For large $k$, the variance tends to zero and $m_{k}$ becomes sharply peaked at $b$. This allows calculation of ${\cal H}_{K}$ in (15), and hence bounds on ${\cal I}_{K}$. For large $K$ then $\langle\chi_{\epsilon}(m_{K}-b)\rangle_{\rm con}\to 1$. Also, from (16) one has $p^{\rm con}_{k}(s|m)=(1+sb)/2$, independent of $m$, so (15) yields

[TABLE]

We used that the fraction of steps with $s_{k}=\pm 1$ is $(1\pm b)/2$ and the contribution to the action for each such hop is $\log[(1\pm b)(1+{\rm e}^{\mp 2Jm})/2]$; also $m_{k}$ is sharply peaked at $b$, so the average action for such a hop can be estimated as $\log[(1\pm b)(1+{\rm e}^{\mp 2Jb})/2]$. Since $m_{k}$ is sharply peaked, we note that this result for the action is somewhat insensitive to details of the controlled dynamics. For example, if the rates in (16) depended on also $m_{k}$, the action would only be sensitive to the values of the rates at the mean value of $m_{k}$. This insight is related to the temporal additivity principle of Harris and Touchette Harris and Touchette (2009); Harris (2015a). In our context, it means that adding extra complexity to the controlled process (16) does not yield improved bounds on ${\cal I}$.

Using again that $m_{K}$ is sharply peaked, the conditional average of the action in (15) can be replaced by the simple average in (18), and one obtains (after simplifying the right hand side of (18) and setting $b=m$)

[TABLE]

as in Klymko et al. (2018). This result is valid for large $K$. It is easily checked that ${\cal H}_{K}(m,\epsilon)$ is non-negative for all $m$ (as it must be, since it is a bound on ${\cal I}$). Using $2\tanh^{-1}m=\log(1+m)/(1-m)$, one also sees that ${\cal H}_{K}(m,\epsilon)=0$ if $m=\tanh(Jm)$. That means that ${\cal I}_{K}(m,\epsilon)=0$ if $m$ is fixed point of (5). For $b\ll 1$ we also obtain

[TABLE]

which determines the action of trajectories with $m\approx 0$. [Recall, we are considering $J>1$ so $m=0$ is an unstable fixed point of (5).]

At this fixed point, one sees that ${\cal I}_{K}(m,\epsilon)\simeq 0$ (for large $K$). This means that trajectories with $m_{K}=0$ have probabilities that do not decay exponentially in $K$. The physical interpretation of this fact is that if the growing cluster contains a symmetric mixture of up and down spins, there is no force that acts to increase or decrease the magnetisation. On the other hand, if $m$ is intermediate between [math] and $m_{\rm S}$, there is a force that drives the system towards the stable fixed point at $m_{\rm S}$. This is the intuition behind the LDP (6): the probability to remain for a long time at a non-typical value of $m$ is suppressed exponentially in $K$, because of the forces in the model that tend to drive $m$ towards a typical value. At $m=0$, the force is zero, so the probability to remain near this value is suppressed less strongly.

So far, all results are fully consistent with KGGW Klymko et al. (2018). We now show that for $m<m_{\rm S}$, the true value of ${\cal I}_{K}$ is much less than ${\cal H}_{K}$ in (19). We will obtain improved bounds by taking a controlled process in (14) in which the transition rates depend explicitly on the step $k$.

III.1.2 Control forces that depend on $k$

We take a controlled process that is a mixture of (16) and (4), as follows. We choose two parameters, which are the bias $b$ in (16) and a step $k^{*}$ at which the controlled dynamics changes its character. For the early part of the trajectory, which is $k\leq k^{*}$, the controlled system is an asymmetric random walk as in (16); for the later part ($k>k^{*}$) we revert to the original dynamics (4). One sees that the action ${\cal A}$ in (13) has no contribution from the later part of the path. Since smaller values of ${\cal A}$ lead to more accurate bounds, this is a desirable feature. We restrict to $|b|<m_{\rm S}$ which is sufficient for our purposes.

Typical trajectories of this controlled dynamics are illustrated in Fig. 2(a). If $b^{2}k^{*}\gg 1$ then one sees from (17) that the distribution of $m_{k^{*}}$ is sharply peaked, in the sense that its mean $b$ is much larger than its standard deviation, which is of order $(k^{*})^{-1/2}$. In this case, the distribution of $m_{k}$ also remains sharply peaked for the later part of the trajectory. For $k>k^{*}$, the system follows the original dynamics and the mean of $m_{k}$ can be obtained from (5) by solving $\mathrm{d}m/\mathrm{d}k=(\tanh Jm-m)/k$, as in Klymko et al. (2018). It is natural to change variables to $u=\log k$ so that

[TABLE]

This means that for $k>k^{*}$, trajectories will flow away from the unstable fixed point at $m=0$ and towards the stable point at $m=m_{\rm S}$, as in Fig. 2(a). Moreover, this evolution is very slow: the natural time variable is not the number of steps $k$ but the rescaled “time” $u=\log k$. (Physically, the slow variation with $k$ occurs because there are already many particles in the cluster, so making a significant change in its magnetisation requires the addition of many particles.)

In order to establish a bound on ${\cal I}_{K}(m,\epsilon)$, we treat $m$ as a target value for $m_{K}$: suitable controlled processes should hit this target with high probability. We choose $b$ and $k^{*}$ to achieve this, as follows.

We solve (numerically) the differential equation (21), going backwards in time. This yields a path which ends at $m_{K}=m$ and can be propagated back to any earlier time $k$, for example by Euler’s method. The magnetisation on this path is $\tilde{b}(k)$ with $k\leq K$ and $\tilde{b}(K)=m$. Both $k$ and $\tilde{b}(k)$ decrease as we solve the equation backwards in time. As $k\to 0$ then $\tilde{b}\to 0$. (Note that $m=0$ is an unstable fixed point of the forward equation, which corresponds to a stable fixed point of the backward equation.) We stop the solution at the point $(k,\tilde{b})=(k^{*},b)$ where

[TABLE]

where $a$ is a numerical parameter, of order unity (we take $a=2$, results depend weakly on this choice).

These $b,k^{*}$ are the parameters that we use for the controlled dynamics. As long as $K$ is reasonably large, the algorithm gives $b\ll 1$ and $k^{*}\gg 1$. Then the action for this controlled process can be estimated from (20), with $K$ replaced by $k^{*}$. This yields

[TABLE]

Since $(b,k^{*})$ solve (22), this means that $\langle{\cal A}\rangle_{\rm con}$ is of order $\sqrt{K}$. Combining this result with (15) and assuming that the controlled system hits the target with probability 1, one obtains

[TABLE]

Hence, the bound ${\cal H}_{K}$ tends to zero as $K\to\infty$. This result applies for large $K$ and is independent of the target chosen for $m_{K}$ (always assuming that this target is between $\epsilon$ and $m_{\rm S}-\epsilon$). Note however, that the controlled process depends on $K$, in that the parameters $b,k^{*}$ are chosen separately for each value of $K$. The assumption that the controlled system hits the target with probability 1 is valid as long as the magnetisation distribution at $k^{*}$ is sharply-peaked in the sense that its mean is much larger than its standard deviation. This requires $b^{2}k^{*}\gg 1$ which is true by (22) as long as $K$ is large.

Numerical results based on this construction are shown in Fig. 2. In particular, Fig. 2(a) shows that the parameters $(b,k^{*})$ obtained by this method are such that the controlled system hits the target $m$ with high probability. Also, Fig. 2(b) confirms that on increasing $K$, the bound ${\cal H}_{K}$ is quantitatively consistent with (24). This establishes that $\lim_{K\to\infty}{\cal H}_{K}(m,\epsilon)=0$, and hence from (14) one has $\lim_{K\to\infty}{\cal I}_{K}(m,\epsilon)=0$, for this value of $m$. Fig. 2(c) shows that the same behaviour occurs for several values of $m$ with $|m|<m_{\rm S}-\epsilon$. This is expected since the theoretical argument above is independent of the target for $m_{K}$. Hence $\lim_{K\to\infty}{\cal I}_{K}(m,\epsilon)=0$ throughout this range, which establishes that the rate function (6) obeys

[TABLE]

This was the result anticipated in Fig. 1. A similar result has been proven by a rigorous analysis of a general class of Polya urns, see Corollary 5 of Franchini (2017). We emphasise here that while the numerical results in Fig. 2 are a useful confirmation of our theoretical calculations, the bound in (24) is an analytical result.

III.2 Scaling of the probability that $|m_{K}|<m_{\rm S}$

As $K\to\infty$, we have shown that ${\cal I}_{K}(m,\epsilon)\to 0$ throughout the regime $|m|\leq m_{\rm S}-\epsilon$. The shows that $\mathrm{Prob}\!\left(|m_{K}-m|<\epsilon\right)$ does not decay exponentially with $K$. Nevertheless, we expect that this probability should vanish as $K\to\infty$, so the natural question is: how small is it?

To address this question we define an estimate of the probability density for $m_{K}$ as

[TABLE]

(In this Section, we emphasise that all large-$K$ limits are to be taken at fixed $\epsilon>0$, note also that we are assuming $J>1$.)

Fig. 3(a) shows the distribution of $m_{K}$, obtained by direct sampling of trajectories of the system, for the representative parameter value $J=1.3$. Two features are notable. First, the probability that $m_{K}\approx 0$ is small and decreases with $K$, but the decay is much slower than exponential in $K$, consistent with the arguments of section III.1. Second, there is no evidence for a local maximum in the probability at the unstable fixed point $m_{K}=0$. We have verified that the behaviour for $|m|<m_{\rm S}$ is similar for larger $J$, so this is a representative state point. However, the behaviour close to the peaks of $\rho_{K}$ has a more complex dependence on $J$, see Sec. III.3 below.

To estimate $\rho_{K}$ for $|m|<m_{\rm S}$, it is possible to repeat the argument of Sec. III.1, replacing (22) with $b^{2}k^{*}=aK^{\alpha}$ for any $\alpha\in(0,1)$. This can be used to show that for any $\beta>0$ one has

[TABLE]

for $|m|\leq(m_{\rm S}-\epsilon)$ as usual. In other words, the probability of a non-typical value of $m_{K}$ decays slower than $\exp(-cK^{\beta})$, for any $c,\beta>0$. Based on this observation, we propose that the probability decays as a power law in $K$. In that case

[TABLE]

should have a positive (non-zero) limit as $K\to\infty$. Fig. 3(b) shows results that are consistent with (28). We now present theoretical arguments that further support this conjecture, including bounds based on (14).

III.2.1 Accurate bounds for large $K$

Consider the same controlled dynamics as in Sec. III.1, but with $b=0$. The remaining parameter is $k^{*}$: this means that the extensive magnetisation $M_{k}$ in the controlled process is a simple (unbiased) random walk for $k<k^{*}$. We use the original growth dynamics (4) for $k>k^{*}$. In this case the distribution of $m_{k^{*}}$ has mean zero and its standard deviation is $(k^{*})^{-1/2}$, by (17). We will take $k^{*}\gg 1$ so the distribution of $m_{k^{*}}$ is sharply-peaked in the sense that its variance is small compared to unity. However, in contrast to Sec. III.1, this distribution does not remain sharply-peaked under the time evolution (see below).

To fix a suitable value for $k^{*}$, it is useful to consider the deterministic evolution of the mean of $m_{k}$ for $k>k^{*}$, assuming that $m_{k^{*}}$ has a typical value of order $(k^{*})^{-1/2}$. Since $m$ is small, one may linearise (21), leading to $(\mathrm{d}m/\mathrm{d}u)=(J-1)m$ and hence $m_{k}=m_{k^{*}}(k/k^{*})^{J-1}$. We choose $k^{*}$ in such a way that this deterministic equation gives $m_{K}=O(1)$. This leads to

[TABLE]

with $\gamma$ a constant of order unity. We take $\gamma=1$.

By analogy with (26), we define an estimate of the probability density for $m_{K}$ under the controlled dynamics as

[TABLE]

where $\mathrm{Prob}_{\rm con}$ indicates a probability under the controlled dynamics. Numerical results for $\tilde{\rho}$ are shown in Fig. 4(a), for several values of $K$, always with $k^{*}$ chosen according to (29). One sees that the distributions of $m_{K}$ are not sharply peaked. Instead $\tilde{\rho}$ is of order unity everywhere between $\pm m_{\rm S}$. Moreover, this distribution depends very weakly on $K$ (which varies over three decades). This is due to the choice proposed in (29): the value of $\gamma$ is not important but the correct power-law exponent is essential. (Different values of $\gamma$ lead to different distributions, but they are all similarly broad.)

Repeating the argument of Sec. II.3, one obtains

[TABLE]

with

[TABLE]

This bound is shown in Fig. 4(b), for the representative case $K=10^{5}$. It shows almost quantitative agreement over a range of $m$ that includes $m=0$. However, the agreement breaks down as $m$ gets close to $m_{\rm S}$. (Better bounds for larger $m$ might be obtained by using smaller $c$ in (29), but we have not explored this in detail. See also Sec. III.3.)

III.2.2 Asymptotic behaviour as $(\log K)\to\infty$

The numerical bound in Fig. 4(b) is valid at $K=10^{5}$ and gives a reasonable estimate of the probability density, but the regime of primary interest is $K\to\infty$. Based on Fig. 4(a), we expect that $\log\tilde{\rho}_{K}(m,\epsilon)$ remains of order unity as $K\to\infty$. In this case, its contribution to $\tilde{\cal H}_{K}$ will decay proportional to $1/(\log K)$ and the term involving ${\cal A}$ will dominate (32) at large $K$. It is useful to estimate this term directly, in order to predict the behaviour of $\tilde{\cal H}_{K}$ as $K\to\infty$. That is, assuming that $\log\tilde{\rho}$ remains finite as $K\to\infty$, one has from (31,32) that

[TABLE]

with

[TABLE]

We refer to this as an asymptotic bound since it does not bound ${\cal J}_{K}$ at any finite $K$, but only as $K\to\infty$: see (33). For finite $K$, the difference between the exact bound and the asymptotic bound scales as $(1/\log K)$, which decays very slowly with $K$. This means that the asymptotic behaviour is not easily accessible in numerics. However, the asymptotic bound can be evaluated numerically for finite $K$ and provides an estimate of the right hand side of (33).

Results are shown in Fig. 4(c). The estimates for the asymptotic bound are compared with the actual distribution for $K=10^{5}$. One sees that the asymptotic bound in Fig. 4(c) depends weakly on $m$, over a fairly wide range that includes $m=0$. This bound is evaluated numerically as an average of ${\cal A}$, which is restricted to trajectories whose final magnetisation $m_{K}$ is within the relevant bin of a suitably constructed histogram. We do not have a theoretical estimate of this conditional average. However, the unconditioned average of the action can be obtained from (13) as

[TABLE]

For $1\ll k\leq k^{*}$ then $m_{k}$ is small under the controlled dynamics and its variance is $1/k$. Hence, we expand (35) to second order in $m_{k}$, which yields $\langle{\cal A}(\bm{m})\rangle_{\rm con}\approx\sum_{k}J^{2}/(2k)$. Approximating the sum by an integral leads to

[TABLE]

where $A_{0}$ is a constant of order unity which depends on the behaviour of trajectories when $k$ is of order unity. (The behaviour of the system for $k=O(1)$ is not captured by our various approximations so we are not able to estimate this constant analytically.)

If we now assume that the the conditional average of the action in (15) is is the same as the corresponding unconditioned average then (29,32,36) together yield

[TABLE]

which is independent of $m$. Using the unconditioned average as an estimate of the conditioned one is an uncontrolled approximation, but (37) is consistent with the weak dependence of $\tilde{H}^{\infty}_{K}$ on $m$ in Fig. 4(c). For that case, (37) evaluates to 0.32, consistent with the data for small and moderate $m$ (given that we are still far from the limit where $\log K$ is large).

Based on this analysis, we summarise our conclusions for the probability distribution of $m_{K}$, at large $K$. We have argued that for $|m|<(m_{\rm S}-\epsilon)$ then ${\cal J}_{K}(m,\epsilon)$ has a finite limit as $K\to\infty$, which means that $\rho_{K}(m,\epsilon)$ decays as a power law. Based on the asymptotic bound in Fig. 4c, we offer two possibilities for the detailed behaviour of ${\cal J}_{K}$. The first is that $\lim_{K\to\infty}{\cal J}_{K}(m,\epsilon)$ is independent of $m$ within some range including zero, and that it takes a value $\alpha$ within that range. From (37), we expect $\alpha\approx J^{2}(J-1)/(2J-1)$ in this case. The result would be that

[TABLE]

within the relevant range of $m$, with $f(m,\epsilon)$ of order unity. The extreme version of this scenario would be that the relevant range is $|m|<m_{\rm S}-\epsilon$: this is not apparent from the finite-$K$ data presented here but we are still far from the limit where $\log K\to\infty$. In any case this limit is not expected to commute with a limit where $m\to(m_{\rm S}-\epsilon)$ so one may expect significant deviations from the asymptotic (large-$K$) behaviour in data obtained at finite $K$. The second scenario is that $\lim_{K\to\infty}{\cal J}_{K}(m,\epsilon)=\alpha(m,\epsilon)$ where $\alpha$ is now a function which takes values of order unity. In this case

[TABLE]

This might be interpreted as a large deviation principle with speed $(\log K)$, because the distribution of $m_{K}$ scales as $p(m_{K})\sim{\rm e}^{-(\log K)\beta(m_{K})}$, where $\beta$ would be interpreted as a rate function [with $\beta(m)=\lim_{\epsilon\to 0}\alpha(m,\epsilon)$]. However, establishing such an LDP would require a detailed mathematical analysis that is beyond the scope of this work.

Based on the available numerical data, we are not able to settle which (if either) of (38,39) is valid, because all results are limited by the fact that the numerical parameter $\log K$ governs the convergence to the large-$K$ regime, and this number is never very large. This might be addressed by noting that for large $k$, the change in $m_{k}$ on any single step is very small, so one might promote both $k$ and $m_{k}$ to continuous variables, and describe the discrete-time stochastic process (4) by a stochastic differential equation. Such a construction would be similar to the temporal addivity principle of Harris and Touchette (2009), but requires a detailed analysis that is beyond the scope of this work.

III.3 Scaling of the probability distribution of $m_{K}$ for $|m_{K}|\gtrsim m_{\rm S}$

We have shown how control forces that depend on $k$ can be used to derive bounds on the rate function, and that these bounds differ qualitatively from the (loose) upper bounds that are obtained using control forces that are independent of $k$ Klymko et al. (2018). All results so far are relevant for the probability distribution of $m_{K}$ when $|m_{K}|<m_{\rm S}$, which is the main focus of this paper.

This Section discusses the behaviour for $|m_{K}|\gtrsim m_{\rm S}$, based on the theory of Harris Harris (2015a). We show how $k$-dependent control forces can yield improved bounds on probability distributions in that case too. In particular, we analyse the behaviour of probability distributions such as the one in Fig. 3, for values of $m$ close to the peaks.

Before presenting the calculation, we first summarise the behaviour that would usually be expected Touchette (2009) based on the LDP (6). One would expect that a central limit theorem (CLT) holds, so that the sharp peaks of $\rho_{K}(m,\epsilon)$ have width (standard deviation) $\sigma K^{-1/2}$ with $\sigma=I^{\prime\prime}(m_{\rm S})^{-1/2}$. However, for $J>1$, one sees from Fig. 1 that $I^{\prime\prime}(m)$ is discontinuous at $m_{\rm S}$. In this situation, one expects in general that $\sigma=I^{\prime\prime}(m_{\rm S}^{+})^{-1/2}$ where the notation indicates that one should evaluate the derivative as $m\to m_{\rm S}^{+}$. The following analysis shows how the behaviour of the growth model differs from this simple picture.

III.3.1 Control forces with arbitrary time-dependence

We apply the theory of Harris (2015a); Harris and Touchette (2009) by connecting it to the controlled process introduced in Equ. (16). This requires that we consider $k$-dependent control forces where the parameter $b$ has an arbitrary dependence on $k$, in contrast to the simple choices considered so far. This leads to $\langle M_{k}\rangle_{\rm con}=\sum_{x=1}^{k}b_{x}$ where the sum runs over steps. Promoting $k$ to a continuous variable we define

[TABLE]

and one sees that $m^{\prime}(k)=[b(k)-m(k)]/k$ where the prime denotes a derivative with respect to $k$. Hence, for any (smooth) path $m(k)$, we can define a a controlled process that follows this path by taking

[TABLE]

Following (18,19), the action for such a time-dependent path can be estimated as

[TABLE]

which assumes that the distribution of $m_{k}$ is sharply peaked and is valid up to an additive constant of order $\log K$. To obtain the best available bound ${\cal H}_{K}(m,\epsilon)$ via (15), this action should now be minimised over the path $m(k)$, subject to (41) and $m(K)=m$.

Here, we take a simpler route. Let $m^{*}$ be a stable fixed point of the deterministic dynamics, so $m^{*}=0$ for $J<1$ and $m^{*}=\pm m_{\rm S}$ for $J>1$. We restrict our analysis to paths for which $m(k)-m^{*}$ and $b(k)-m^{*}$ are both small compared to unity. This allows us to apply the results of Harris (2015a). Under these assumptions the action (42) may be expanded as

[TABLE]

where we used (41); and $\Lambda=1-J[1-(m^{*})^{2}]$ and $\chi=1-(m^{*})^{2}$ are numerical parameters that correspond to $(1-A^{*})$ and $D^{*}$ in Equ. (7) of Ref. Harris (2015a). We neglected contributions at ${\cal O}(m-m^{*},b-m^{*})^{3}$ in the integrand. For forces that are independent of $k$ the action is $\langle{\cal A}\rangle_{\rm con}\approx{\cal A}_{\rm ind}$ with

[TABLE]

Using this result with (15,14) yields a bound on ${\cal I}_{K}$ which was also derived in KGGW Klymko et al. (2018).

The next step is to minimise the action $\langle{\cal A}\rangle_{\rm con}$ over paths $m(k)$ that satisfy $m(K)=m$. The resulting action yields a bound on ${\cal I}_{K}(m,\epsilon)$. We refer to Ref. Harris (2015a) for all details of the minimisation procedure. The important physical conclusion is that the solutions do not have constant forces, so one finds in general that $\langle{\cal A}\rangle_{\rm con}<{\cal A}_{\rm ind}$.

The reduction in the action depends strongly on the parameter $\Lambda$. One always has $0<\Lambda<1$, and the results depend qualitatively on whether $\Lambda$ is bigger or smaller than $\frac{1}{2}$ Harris (2015a). Small values of $\Lambda$ correspond to systems with large fluctuations.

III.3.2 Results for the growth model

In the irreversible growth model, we find that $\Lambda>\frac{1}{2}$ if the parameter $J$ is far from its critical value $J_{c}=1$. Specifically, $\Lambda>\frac{1}{2}$ if either $J<0.5$ or $J>1.37$. In these cases, the theory Harris (2015a) gives a minimal action

[TABLE]

Comparing with (44), one sees that introducing $k$-dependent control forces reduces the action by a factor of $(2\Lambda-1)/\Lambda^{2}$.

However, it is important to consider the optimal paths that are predicted by this theory. For convenience we restrict to the case where $m,m^{*}>0$. The optimal paths take $m(k)=m^{*}$ for $k<k_{*}$ where $k_{*}\ll K$ is a parameter. For $k>k_{*}$ then $m(k)$ has an excursion away from $m^{*}$, it grows to a value of order $(K/k_{*})^{1-\Lambda}(m-m^{*})$ before converging towards the target as $m(k)-m^{*}\approx(K/k)^{1-\Lambda}(m-m^{*})$. Note however that we required $m(k)-m^{*}\ll 1$ when deriving (43): hence one must have $(m-m^{*})={\cal O}(k_{*}/K)^{1-\Lambda}$. Since the derivation of (45) also requires that $k_{*}\ll K$, one sees that $m-m^{*}$ must be small.

For this reason, we focus on values of $m$ that are very close to the peaks of $\rho_{K}(m,\epsilon)$. We take $m-m^{*}=vK^{-x}$ for parameters $v,x>0$, so the difference between $m$ and $m^{*}$ vanishes as $K\to\infty$. In addition we assume that $x\geq(1-\Lambda)$: this includes the possibility that $x=\frac{1}{2}$ because $\Lambda>\frac{1}{2}$. In this case our various assumptions are all satisfied and (45) yields

[TABLE]

independent of $k_{*}$. Since all the assumptions of Harris (2015a) are satisfied, we expect that this action is optimal in the sense that substituting in (15) should yield a bound ${\cal H}_{K}$ that achieves equality in (14), as $K\to\infty$. The case $x=\frac{1}{2}$ is instructive because the resulting action is of order unity as $K\to\infty$: this indicates that fluctuations on this scale are typical and that the peaks of $\rho_{K}(m,\epsilon)$ can be approximated by Gaussian probability densities. The variances of these Gaussian peaks are predicted to be

[TABLE]

(Note, in contrast to previous results that described large deviations, this result is relevant for typical fluctuations.)

As an illustrative example of the behaviour near the peak of $\rho_{K}$, Fig. 5 shows results for $J=1.6$, for which $\Lambda\approx 0.67$. The data obtained by direct sampling are well-described by a Gaussian distribution with variance as in (47), which differs significantly from the corresponding prediction that was obtained using control forces that are independent of $k$. (This is the prediction obtained by using (44) as a bound and assuming that this bound coincides with the rate function).

Note that the dependence of (47) on $K$ resembles a CLT. However, it describes the behaviour of a single peak, while the distribution $\rho_{K}$ has two peaks; the true variance of $m_{K}$ is close to $m_{\rm S}^{2}$ which is of order unity. It is also possible to consider the distribution of $|m_{K}|$, which has a unimodal probability density function whose peak is well-described by a Gaussian with variance (47), recall Fig. 5. In this case we still find numerically (for $J=1.6$) that the true variance of $|m_{K}|$ is much larger than (47). This effect can be explained using the results of Sec. III.2: the distribution of $|m_{K}|$ has a heavy tail that extends all the way to $|m_{K}|=0$, and this has a significant contribution to the variance (recall Fig. 3). We conclude that (47) predicts the width of the peak of $\rho_{K}(m,\epsilon)$ but it is not sufficient to establish a CLT, neither for $m_{K}$ nor for $|m_{K}|$. On the other hand, for $J<0.5$ there is no heavy tail and we do expect a CLT for $m_{K}$. The results of Sec. III.2 also indicate that a CLT for $|m_{K}|$ might be expected for larger values of $J$, but this is beyond the scope of this work.

The result (47) raises the question of what happens as $\Lambda\to\frac{1}{2}$, where this variance seems to diverge. To understand this, we consider the case $\Lambda<\frac{1}{2}$, for which the physics changes qualitatively. Assuming as before that $K\gg k_{*}$, the analogue of (45) is Harris (2015a)

[TABLE]

In the optimal path that gives this result Harris (2015a), the maximal value of $m(k)$ is of order $(K/k_{*})^{\Lambda}(m-m^{*})$. Following the same procedure as before, this suggests that we take $m-m^{*}=vK^{-x}$ with the restriction $x\geq\Lambda$ (which again includes $x=\frac{1}{2}$ because we are now considering $\Lambda<\frac{1}{2}$). The result is that

[TABLE]

Note that $k_{*}$ is still a free parameter which we assume to be of order unity. For the case $x=\frac{1}{2}$ one sees that $\langle{\cal A}\rangle_{\rm con}\to 0$ for small $K$, in contrast to $\langle{\cal A}\rangle_{\rm con}={\cal O}(1)$ in (46). In fact $\langle{\cal A}\rangle_{\rm con}\to 0$ for all $x<\Lambda$, independent of $k_{*}$. For $x=\Lambda$ we expect an action of order unity. Obtaining its numerical value would require optimisation of the parameter $k_{*}$, but this depends on contributions to the action from steps with $k=O(1)$, which are not captured by our various approximations. However, since the action is of order unity for fluctuations of this size, so we expect that $\rho_{K}(m,\epsilon)$ has peaks with width of order $K^{-\Lambda}$. That is,

[TABLE]

which corresponds to peaks that are much broader than (47). This is consistent with (47): the apparent divergence of that variance as $\Lambda\to\frac{1}{2}$ signals a change in scaling behaviour from $K^{-1}$ to $K^{-2\Lambda}$. This means that the variance always tends to zero on taking $K\to\infty$ at fixed $\Lambda$. Equ. (50) should apply (for example) to the variances of the individual peaks in Fig. 3, since $\Lambda\approx 0.43<\frac{1}{2}$ in that case. We also predict that it should hold in the one-phase regime, for $0.5<J<1$: this behaviour was demonstrated for the (so-called) artificial model of Harris (2015b), whose behaviour is closely related to the model considered here. A more detailed analysis of (50) for the growth model would include numerical tests and determination of the prefactor in the scaling law, but this is beyond the scope of this paper.

In addition to the results derived so far, the arguments of Harris (2015a) also suggest that the behaviour in (50) should be associated with a new LDP whose speed is less than $K$. The calculations presented here are not sufficient to demonstrate this, because our analysis of optimal paths is restricted to very small $m-m^{*}$. However, we summarise the physical picture that may be expected, if the arguments of Harris (2015a) can be extended to situations where $m-m^{*}={\cal O}(1)$, for this model. For $\Lambda>\frac{1}{2}$ one expects the LDP (6) and also (47) to hold. For $\Lambda<\frac{1}{2}$ then one expects from (48) that $I(m)=0$ for all $m$ and one should have a different LDP with speed $K^{2\Lambda}$:

[TABLE]

In this case, (50) arises naturally, and the numerical value of the variance is related to the second derivative of $I_{\Lambda}$. For $J>1$ and $|m_{K}|<m_{\rm S}$, we recall from Sec. III.2 that $p(m_{K})$ always decays as a power law in $K$, which means that one still has $I_{\Lambda}(m_{K})=0$ in this range. Hence the rate function $I_{\Lambda}$ would be qualitatively similar to Fig. 1.

This is an appealing physical picture, which is consistent with the general expectations of Harris (2015a). However, we emphasise again that it would require significant further mathematical analysis to confirm it. For example, our results here are also consistent with a scenario where (6) holds for all values of $\Lambda$, with $I(m)>0$ for all $|m|>m_{\rm S}$. In this case (50) would imply that $I^{\prime\prime}(0)=0$ for $0.5<J<1$ and $I^{\prime\prime}(m_{\rm S}^{+})=0$ for $1<J<1.37$. For other values of $J$ then these second derivatives would be presumably be non-zero.

IV Results : Reversible model

In this section we analyse large deviations in the reversible model of growth, which was defined in Sec. II.2. This demonstrates that the use of time-dependent control forces to obtain improved bounds can be generalised beyond the single system considered so far, leading to results similar to Fig. 1. We focus on a representative state point where the deterministic dynamics has unstable fixed points at non-zero magnetisation. This is the case where the model has three steady states. (For parameters where the model has one steady state or two steady states, the behaviour is very similar to that of the irreversible model, so we do not discuss it in detail.)

IV.1 Summary of behaviour and definition of controlled dynamics

We summarise the behaviour of the model, as derived in Klymko et al. (2017). We assume throughout that the parameter $c$ in (10) is large enough that the steady state of the system involves a cluster that is grows with a constant rate. (The alternative is that there is a steady state with a cluster of finite size, the behaviour is quite different in that case and we do not discuss it here.) For the reversible model, the analogue of (5) is

[TABLE]

where the average is conditioned on the values of $m_{k}$ and $N_{k}$. As noted in Klymko et al. (2017), this means that steady state values for $m$ solve

[TABLE]

For $J>c$ there are three solutions to this equation, and the behaviour is similar to the irreversible model for $J>1$. For $J<c$ there are two possibilities – either a single solution, or five solutions. The first case corresponds to a single steady state with $m=0$. The latter case means that there are three possible steady states, which have magnetisations $0,\pm m_{\rm S}$. In this case there are also two unstable fixed points of (52) at $\pm m_{\rm U}$ (these are not steady states of the stochastic model, similar to the irreversible case). The situation with three steady states is only possible for $\sqrt{6}<J<c$ (this condition is necessary but not sufficient).

As a suitable controlled dynamics, we take

[TABLE]

with two parameters $\lambda,c^{\prime}$, also $z_{k}^{\prime}=2(c^{\prime}+1-\lambda m_{k})$ for normalisation. This controlled process is similar to that considered in Klymko et al. (2018), but not identical. In principle one may take different rates for the two cases with $\Delta N_{k}=+1$, which will lead to improved bounds in some cases. However, (54) is sufficient for our purpose, which is to show that $I(m)=0$ for $|m|\leq m_{\rm S}$.

In the steady state of this model, the rate of cluster growth is denoted by $r$ and the magnetisation is denoted by $b$ (by analogy with (17)). Considering the mean increments for $N_{k}$ and $M_{k}$, the rate of cluster growth can be verified to be $r=(c^{\prime}-1+\lambda b)/(c^{\prime}+1-\lambda b)$, which is assumed to be positive, as noted above. Also, $b$ is a solution of $\lambda(b^{2}-1)+bc^{\prime}=0$, specifically

[TABLE]

If $\lambda=0$ (equal probabilities for unbinding of up and down spins) then $b=0$ (no magnetisation in steady state).

Using that the cluster size and the magnetisation are both sharply peaked in the steady state, the analogue of (18) may be expressed as

[TABLE]

with $z^{\prime}_{b}=2(c^{\prime}+1-\lambda b)$. We emphasise that this formula is valid only if $b,c^{\prime},\lambda$ are related by (55), so that $b$ is the magnetisation of the steady state of the controlled process.

IV.2 Bounds on ${\cal I}_{K}(m,\epsilon)$

The main result of this section is that ${\cal I}_{K}(m,\epsilon)\to 0$ for $|m|<m_{\rm S}$ and $K\to\infty$, and hence that the rate function $I(m)$ in (9) is zero thoughout this range. This is illustrated in Fig. 6(a), which is similar to Fig. 1. As in the irreversible case, we prove this by obtaining bounds on ${\cal I}_{K}$.

As in Sec. III.1.1, we first consider controlled processes with constant rates (independent of $k$). This will recover results similar to KGGW. We use (56) with (14) to derive bounds on ${\cal I}_{K}(m,\epsilon)$, for $m=b$. Using

[TABLE]

from above and also substituting for $z^{\prime}_{b}$, the formula (56) for the action may be expressed as a function of $(b,\lambda)$. Minimising this function (numerically) over $\lambda$, one obtains a bound on ${\cal I}_{K}$, which is plotted with squares in Fig 6(a). Note that if $m=b$ is a solution of (53) then one may take $\lambda=\tanh Jb$ and $c^{\prime}=c/(\cosh Jb)$, and the action (56) evaluates to zero. Hence, the resulting bound on ${\cal I}_{K}$ is zero for fixed points of the deterministic dynamics (including the unstable fixed points). It is positive for other values of $m$.

If the growth process has multiple stable states, this bound (obtained with control forces independent of time) does not provide an accurate estimate of ${\cal I}_{K}$. The reason is similar to the irreversible model. To derive improved bounds we use the fact that the average action is very small if $b$ is close to $m_{\rm U}$. We then follow the method that was illustrated in Fig. 2. We introduce the additional parameter $k^{*}$ and we construct trajectories that have magnetisation $b$ for $k<k^{*}$, but then follow the natural dynamics of the model for $k>k^{*}$. The parameters $b,k^{*}$ are chosen such that the final magnetisation $m_{K}$ is close to its target value $m$. The action can then be minimised by taking $b$ close to $m_{\rm U}$ and $k^{*}/K$ to be small.

Compared to the irreversible model, the method for choosing $b,k^{*}$ is different; one must also fix $\lambda$, which determines $c^{\prime}$ via (57). Define $n_{k}=N_{k}/k$. Following Klymko et al. (2017), we derive analogues of (21), which are differential equations for the mean of $n_{k}$ and $m_{k}$ as a function of $u=\log k$. These are

[TABLE]

and

[TABLE]

These equations will be solved forwards in time, in contrast to the irreversible model. (The reason is that $m$ prescribes a target for $m_{K}$ but there is no target for $n_{K}$, so there are not enough boundary conditions to solve backwards in time.)

Suppose that we are given some $K$ and some target $m$ for $m_{K}$. We propose an initial guess for $b$ that is between $m_{\rm U}$ and $m$. For this $b$, we minimise (56) over $\lambda$ as above, to obtain a controlled process with $b$ as the steady state magnetisation. This minimisation fixes the parameter $\lambda$ and hence also $c^{\prime}$. We then use the steady state of this controlled process as an initial condition and solve (58,59) forwards in $u$, starting from an (arbitrary) initial value $u_{0}$. We stop the solution when the magnetisation hits the target. This happens at some $u=u_{1}$ and we set $u_{1}=\log K$ so that the magnetisation will hit the target at the required time. This requires that we identify $k^{*}=\log u_{0}$ as the point when we remove the control forces and allow the system to start evolving according to (58,59). Hence $k^{*}=K{\rm e}^{u_{0}-u_{1}}$. Given the initial choice $b$, this yields a value for $k^{*}$ such that the system with time-dependent control forces will hit the target $m$. The average action for this process is given by the product of (56) and $k^{*}$, which is straightforward to evaluate.

It remains to optimise the choice of $b$. By analogy with (22,23), we choose this parameter such that the average action for the controlled process is

[TABLE]

where $a_{A}$ is a parameter of order unity (we take $a_{A}=0.1$). Given a target $m$ (with $|m|<m_{\rm S}$) and a (sufficiently-large) value of $K$, it is possible to choose an initial guess for $b$ such that (i) the left hand side of (60) is larger than the right hand side, and (ii) the left hand side is reduced by moving $b$ closer to $m_{\rm U}$. Then, one may move $b$ towards $m_{\rm U}$ in suitably-chosen steps until one finds parameters $(b,k^{*})$ that solve (60). The method for computation of $\langle{\cal A}\rangle_{\rm con}$ also fixes the values for $(\lambda,c^{\prime})$, as described above.

In this procedure, there is only one pitfall, which is similar to the irreversible case: Eqs. (58,59) are only applicable if the distribution of $m_{k^{*}}$ is sharply-peaked. (Specifically, its mean $b$ should differ from $m_{\rm U}$ by an amount that is much larger than its standard deviation, which is of order $1/\sqrt{k^{*}}$). This requirement is always satisfied if $K$ is large enough.

Combining the ingredients, it follows that for large-$K$, we have established a bound on ${\cal I}_{K}$ that scales as

[TABLE]

This bound tends to zero at large $K$. We emphasise that while this procedure for fixing $b,k^{*}$ is numerical, it does not involve any simulation of the growth process, only minimisation of (56) and numerical solution of (58,59).

For cases where the model has two steady states, this method operates in the same way as the irreversible model and yields the same results. (The unstable fixed point has $m_{\rm U}=0$ in this case.) By (61) one has ${\cal H}_{K}(m,\epsilon)\to 0$ as $K\to\infty$ so the rate function in (9) reduces to $I(m)=0$ for $|m|\leq m_{\rm S}$. All details of this computation are very similar to the irreversible model: for reasons of brevity we do not show numerical results in this case.

The more interesting situation occurs when the model has three steady states. In this case we have performed simulations of the controlled process, in order to obtain bounds on ${\cal I}_{K}$. Results are shown in Fig. 6(b,c,d), following the procedure given above for determination of $(b,k^{*},\lambda,c^{\prime})$. The trajectories of the controlled process remain close to the unstable fixed point for $k<k^{*}$, and diverge from it at later times. Depending on the value of $b$, they may be attracted towards the fixed point at [math], or the one at $m_{\rm S}$. In either case, the method yields results that are consistent with (61) and sufficient to establish that ${\cal I}_{K}(m,\epsilon)\to 0$ as $K\to\infty$. Hence $I(m)=0$ whenever $|m|\leq m_{\rm S}$, as in the irreversible case.

It would be interesting to investigate further the scaling of the probability density $\rho_{K}$ for these values of $m$, as in Sec. III.2. We anticipate similar results to that section, but a detailed analysis is beyond the scope of this work.

V Discussion

These models of cluster growth show rich and interesting behaviour, both for typical trajectories Klymko et al. (2017) and for large deviations Klymko et al. (2018). They describe well-mixed clusters, in the sense that growth rates depend on the mean magnetisation and not, for example, the magnetisation near the boundary of the cluster. This often results in a self-averaging property, so that fluctuations are small when clusters are large. However, in cases where the deterministic dynamics has multiple fixed points (including unstable ones), large fluctuations are still possible, because trajectories may remain close to the unstable fixed point for large times, before eventually leaving it and converging (slowly) to a stable steady state. See Fig. 2(a) and Fig. 6(b).

At the level of large deviations, these trajectories are associated with large fluctuations and manifest in a rate function $I(m)$ that is zero whenever $|m|<m_{\rm S}$, recall Fig. 1. The probability to find a magnetisation in this range is not suppressed exponentially in $K$. For the irreversible model, we have shown in Sec. III.2 that these probabilities decay as power laws in $K$ and we expect similar behaviour for the reversible model too.

We have also emphasised that the models are not ergodic and do not fit the classes considered by Lecomte et al. (2007); Chétrite and Touchette (2015a). They can be expressed as non-Markovian processes and analysed using methods from Harris and Touchette (2009); Harris (2015a). As noted in those works, this can lead to complex behaviour even for typical fluctuations, see Sec. III.3.

It is also useful to recall the analogy between large deviation theory and equilibrium thermodynamics Lecomte et al. (2007); Garrahan et al. (2009); Touchette (2009). Within this framework, probabilities of individual trajectories (in $d$ dimensions of space and 1 dimension of time) are analogous to configurations of equilibrium systems in spatial $d+1$ dimensions. The growth models have no spatial degrees of freedom so $d=0$; this means that trajectories of the growth model are analogous to configurations of a one-dimensional Ising model, where $s_{k}$ is interpreted as the $k$th spin in the chain. The energy of a configuration of this Ising model is

[TABLE]

with $\bm{s}=(s_{1},s_{2},\dots,s_{K})$ and $m_{k-1}$ given by (3); also $E_{0}$ is an arbitrary additive constant and the analogy requires that the temperature $T=1$. Equ. (62) corresponds to an Ising model with long-ranged interactions, while the standard Markovian class of models would have only nearest-neighbour interactions.

Within this analogy, the rate function in the LDP corresponds to a free energy in the equilibrium system. This means that the results of Fig. 1 and Fig. 6(a) somewhat resemble a double-tangent construction, which would usually be associated with phase coexistence (it is equivalent to the Maxwell construction, see Sec. 4.7 of Huang (2010)). The analogy between dynamical large deviations and phase coexistence is discussed, for example, in Hedges et al. (2009); Elmatad et al. (2010); Nemoto et al. (2017). The physical analogue of phase coexistence in dynamical trajectories may depend on system details, but one possibility is that trajectories that realise the rare event of interest show different behaviour in early-time and late-time regimes, as in Fig. 2A of Elmatad et al. (2010) and Fig. 4B of Hedges et al. (2009).

In the growth models considered here, the analogy with equilibrium phase coexistence is not complete because of the long-ranged interactions in the Ising energy (62). From Fig. 2(a), one sees that the trajectories that realise the relevant rare events have qualitatively different behaviour in the early-time regime ($k<k^{*}$) and the late-time regime ($k>k^{*}$), similar to the behaviour for Markovian models Hedges et al. (2009); Elmatad et al. (2010). However the behaviour in the late-time regime is not at all stationary, for example the typical magnetisation $m_{k}$ depends on $k$ throughout the range $k^{*}<k<K$. This is contrary to the behaviour in Markovian models Hedges et al. (2009); Elmatad et al. (2010); Nemoto et al. (2017) where averages depend weakly on time within the late-time and early-time regimes, even if they differ strongly between these two regimes. For this reason, we prefer not to use the terminology of phases and phase coexistence to describe the behaviour shown in Fig. 1. Nevertheless, the behaviour of the rate function is the same as one would obtain from a double-tangent construction.

We also emphasise that while the rate functions for these growth models are never concave, the double-tangent construction does not hold generally in non-Markovian systems Duffy and Sapozhnikov (2008) nor even in Markovian systems on non-compact state spaces Nickelsen and Touchette (2018) – in such cases, the applicability of the double-tangent construction has to be tested on a case-by-case basis. This is similar to analysis of thermodynamic phase coexistence in systems with long-ranged interactions, where the applicability of the Maxwell construction depends on the decay of the interaction potential Campa et al. (2009). For the growth models considered here, we have shown that the construction is applicable.

Finally, we comment on the usefulness of the bound (14) for numerical estimation of small probabilities, as in Klymko et al. (2018); Whitelam (2018a). Our results here confirm that suitable choices of the controlled dynamics can make this bound accurate (see for example Figs. 4b and 5). However, construction of the relevant controlled dynamics in those cases required detailed understanding of the dynamical behaviour of the model (including analytical estimates of the action). Our conclusion is that this method is only reliable if one already has a precise understanding of the mechanism by which the relevant rare events (large deviations) will occur. In this case, one may design a controlled process with this mechanism in mind. However, experience with a range of model systems (see for example Sec. 3.4 of Jack and Sollich (2015)) indicates that it is difficult to predict suitable controlled dynamics, without prior theoretical analysis. If one evaluates the bound (14) using a controlled process does not fully account for the mechanism of the rare event, one may expect to obtain bounds that are not accurate estimates of the probabilities of interest.

Acknowledgements.

I would like to thank Steve Whitelam for many interesting discussions about growth models and numerical sampling methods. I am also grateful to Rosemary Harris, Hugo Touchette, and Juan P. Garrahan for useful discussions, including those related to LDPs in non-Markovian processes, and the possibility of optimal control forces that depend explicitly on time.