Large deviations in boundary-driven systems: Numerical evaluation and effective large-scale behavior

In many cases the typical size of fluctuations in a physical system with N degrees of freedom is of order $1/\sqrt {N}$. Larger fluctuations are rare, and their probability scales as $P\left [\rho \right ]\sim \exp \left (-N\phi \left [\rho \right ]\right)$, where ϕ is an intensive function of the state ρ. The function $\phi \left [\rho \right ]$ is known as the large-deviation function (LDF) and is of fundamental interest in statistical mechanics. In equilibrium systems, ϕ is equal to the free-energy density. Away from equilibrium, a simple expression for ϕ is in general not known, and it may be affected by details of the system's dynamics. Besides its fundamental interest for non-equilibrium physics, the function ϕ is important in various applications, e.g., for calculating escape rates from metastable states, with applications ranging from chemistry and population dynamics to cosmology [1–7].

In a non-equilibrium steady state, to compute the probability of a rare event, one must calculate the dynamics leading up to that event [8]. This is in general a difficult task, even more so for spatially extended systems, where only a handful of analytical solutions exist [9]. If a general understanding is to emerge, additional methods beyond exact solutions need to be considered. Indeed, recent years have seen a considerable effort to develop numerical techniques to calculate the LDF in a variety of systems [10–14].

In this letter, we study the LDF in bulk-conserving diffusive systems, which are driven out of equilibrium by the boundaries. These describe a broad range of transport phenomena, including electronic systems, ionic conductors, and heat conduction [15, 16]. We show that the LDF in such systems can be efficiently evaluated numerically for a general interacting system in one and two dimensions, giving us access to previously unavailable information. This is done by searching for the most probable history $\rho \left (x,t\right)$ of the conserved density function ρ leading to a rare state $\rho _{f}(x)$. Importantly, using the numerical technique we show that for many non-trivial cases, the LDF of a slowly varying configuration $\rho _{f}(x)$ can be calculated by considering only histories $\rho \left (x,t\right)$ which are slowly varying in space, i.e., which are given by the sum of only a few long-wavelength modes. This implies that the long-wavelength structure of the LDF can be understood using an effective finite-dimensional theory, instead of the full infinite-dimensional one. In addition, we find that a local-equilibrium assumption can capture much of the short-wavelength structure. This could suggest a simple framework to treat the LDF in these systems.

In bulk-conserving diffusive systems, the conserved density $\rho \left (\mathbf {x},t\right)$, representing, e.g., charge or energy density, is related to the current $\mathbf {J}\left (\mathbf {x},t\right)$ by

$$\begin{equation} \partial_{t}\rho+{\mathbf{\nabla}}\cdot{\mathbf{J}} = 0, \end{equation}\noindent \tag{ 1 }$$

where the current is given by

$$\begin{equation}{\mathbf{J=-}}D(\rho({\mathbf{x}},t)){\mathbf{\nabla}}\rho({\mathbf{x}},t)+\sqrt{\sigma(\rho({\mathbf{x}},t))}\boldsymbol{\eta}({\mathbf{x}},t). \end{equation}\noindent \tag{ 2 }$$

$D(\rho (x,t))$ is a density-dependent diffusivity function, while $\sigma (\rho ({\mathbf {x}},t))$ controls the size of the white noise $\boldsymbol {\eta }({\mathbf {x}},t)$, which in d dimensions satisfies $\langle \eta _{a}({\mathbf {x}},t)\rangle = 0$ and $\langle \eta _{a}({\mathbf {x}},t)\eta _{b}({\mathbf {x}}^{\prime},t^{\prime})\rangle = N^{-1}\delta _{ab}\delta ^{d}({\mathbf {x}}-{\mathbf {x}}^{\prime})\delta(t-t^{\prime})$. The prefactor N⁻¹ in the noise variance results from the fact that we have scaled distances by L = N^1/d the system size, and time by L². After this rescaling the noise is small as a consequence of the coarse graining. $D(\rho)$ and $\sigma (\rho)$ are related via a fluctuation-dissipation relation (Nyquist noise in electronic systems), which for particle systems reads $\sigma (\rho) =2k_{B}T\rho ^{2}\kappa (\rho)D(\rho)$, where $\kappa (\rho)$ is the compressibility [9]. Here we study a system on a domain A connected to reservoirs which fix the density at the boundary ∂A,

$$\begin{equation} \left.\rho({\mathbf{x}},t)\right.\vert_{\mathbf{x}}\in\partial A = \left.\rho_{B}({\mathbf{x}})\right.\vert_{{\mathbf{x}}\in\partial A}. \end{equation} \tag{ 3 }$$

If the boundary density is not constant, a current is induced through the system, driving it out of equilibrium.

The average, or most probable density profile for the system $\bar {\rho }$, is obtained by solving ${\mathbf {\nabla }}\cdot [D(\bar {\rho }) {\mathbf {\nabla }}\bar {\rho }] = 0$, with $\bar {\rho }\left ({\mathbf {x}} \right) = \rho _{B}({\mathbf {x}})$ at the boundaries. In equilibrium (i.e., when ρ_B is constant), the steady-state probability of any other density profile $\rho ({\mathbf {x}})$ is easy to obtain: the large-deviation functional $\phi \left [\rho \right ]$ is then given by the free energy, which is local in ρ. By contrast, the steady-state probability distribution away from equilibrium is notoriously hard to compute. In general, it is known that despite the local nature of the dynamics, $\phi \left [\rho \right ]$ is non-local in ρ, leading to generic long-range correlations [9]. Analytical results for $\phi \left [\rho \right ]$ are known only for a few models of interacting systems, corresponding to specific choices of $D(\rho)$ and $\sigma (\rho)$, and almost exclusively in one dimension (1d). The only known example in higher dimensions is the zero-range process, which admits a trivial product measure [17, 18]. Hence, despite the central role that $\phi \left [\rho \right ]$ plays in the understanding of non-equilibrium phenomena, little is known about its properties.

To compute the large deviation for the model described above, we first note that the probability of a noise realization $\boldsymbol {\eta }({\mathbf {x}},t)$ is Gaussian, $P\sim \exp \left (-N\frac {1}{2} \int \boldsymbol {\eta }^{2}{\mathrm { d}}{\mathbf {x}}{\rm d}t\right)$. Using this expression together with eq. (2), the probability of a history $\left \{ \rho ({\mathbf {x}},t),{\mathbf {J}}({\mathbf {x}},t) \right \}$ during time τ ⩽ t ⩽ 0 is $P\sim \exp (-NS)$, where the action S is given by

$$\begin{equation} S = \int_{\tau}^{0}{\rm d}t\int_{A}{\rm d}{\mathbf{x}}\frac{\left[{\mathbf{J}}({\mathbf{x}},t)+D(\rho({\mathbf{x}},t)) {\mathbf{\nabla}}\rho\left({\mathbf{x}},t\right)\right]^{2}}{2\sigma(\rho({\mathbf{x}},t))}. \end{equation} \tag{ 4 }$$

As the noise is small, the system spends most of the time close to $\bar {\rho }$, the unique fixed point of the diffusion equation at zero noise. To calculate the steady-state probability of a rare event, $\rho _{f}({\mathbf {x}})$, we consider trajectories $\rho \left ({\mathbf {x}} ,t\right)$ starting from $\bar {\rho }$ in the distant past $t\rightarrow -\infty$ and ending at $\rho _{f}({\mathbf {x}})$ at t = 0. For large N, its probability $P\sim \exp \left \{-N\phi \left [\rho _{f}\right ] \right \}$ is given by [8, 18–21]

$$\begin{equation} \phi\left[\rho_{f}\right] = \inf_{\rho,\mathbf{J}}S, \end{equation} \tag{ 5 }$$

where the infimum is over histories satisfying eq. (1), with initial and final conditions $\rho \left ({\mathbf {x}},t\rightarrow -\infty \right) = \bar {\rho }({\mathbf {x}})$, $\rho ({\mathbf {x}},t = 0) = \rho _{f}({\mathbf {x}})$, and the boundary conditions, eq. (3). We now describe a numerical method which utilizes this formulation to evaluate the large deviation $\phi \left [ \rho _{f}\right ]$.

To calculate numerically the large deviation, we directly minimize the action, eq. (4). We present a simple algorithm, which efficiently finds minima of the action for such systems. The algorithm is based on starting with a problem where the solution is known and gradually modifying it while maintaining a minimizing solution. Specifically, we start with a problem where the initial and final states are identical, $\rho \left ({\mathbf {x}},t\rightarrow \tau \right) =\rho ({\mathbf {x}},t=0) =\bar {\rho }({\mathbf {x}})$. In the exact solution τ = −∞; here we take a sufficiently early time τ, before any significant evolution has begun, and check for convergence. Equation (4) has a unique minimum to this problem, $\rho ({\mathbf {x}},t) = \bar {\rho }({\mathbf {x}})$ and ${\mathbf {J=-}}D(\bar {\rho })\mathbf {\nabla }\bar {\rho }$, for which the action vanishes, hence $\phi \left [\bar {\rho }\right ] = 0$. We now gradually change the final condition: define a series of gradually changing profiles $\rho^{(m)}({\mathbf {x}})$, m = 1,...,n, with $\rho^{(1)}={\bar{\rho}}$, and $\rho^{(n)}=\rho_{f}$. We call the series $\{\rho ^{(m)}\}$ the final-state trajectory. The solution for the problem with $\rho ^{(1)}({\mathbf {x}})$ is known. In the next iteration we solve the same problem as above, only with final conditions $\rho ({\mathbf {x}},t=0) = \rho ^{(2)}({\mathbf {x}})$, using as an initial guess $\rho \left ({\mathbf {x}} ,t\right)$ as obtained from the previous iteration. This procedure is iterated until we reach the final condition $\rho ^{(n)}({\mathbf {x}}) = \rho _{f}({\mathbf {x}})$. Standard algorithms can be used for the minimization at each step; we have experimented with different algorithms, and the final results did not depend on this choice, but computational efficiency did. Since S is a sum of squares, non-linear least-squares methods are applicable and were found to be efficient.

Within each iteration, we minimize the action $\inf _{\rho ,{\mathbf {J}}}S\left [ \rho ,{\mathbf {J}}\right ] = \inf _{\rho }\tilde {S}\left [\rho \right ]$, where $\tilde {S}\left [\rho \right ]\equiv \inf _{\mathbf {J}}S\left [\rho ,{\mathbf {J}}\right ]$. To evaluate $\tilde {S}\left [\rho \right ]$ we take into account the constraint given by eq. (1), by introducing a Lagrange multiplier ${\hat {\rho }}({\mathbf {x}},t)$, and optimizing $G=S+\int {\hat {\rho }}({\mathbf {x}},t)(\partial _{t} \rho +{\mathbf {\nabla }}\cdot {\mathbf {J}}){\mathrm { d}}{\mathbf {x}}\,{\rm d}t$ with respect to J and ${\hat {\rho }}$. δG/δJ = 0 gives ${\mathbf {J}} =-D(\rho){\mathbf {\nabla }}\rho +\sigma (\rho){\mathbf {\nabla }}{\hat {\rho }}$, which together with eq. (1) reads

$$\begin{equation} \partial_{t}\rho={\mathbf{\nabla}}\cdot(D(\rho){\mathbf{\nabla}}\rho) - {\mathbf{\nabla}}\cdot\left(\sigma(\rho){\mathbf{\nabla}}{\hat{\rho}}\right), \end{equation} \tag{ 6 }$$

with boundary conditions $\left .{\hat {\rho }}\right .\vert _{{\mathbf {x}} \in \partial A} = 0$ [17]. This is a linear equation for ${\hat {\rho }}({\mathbf {x}},t)$ in terms of $\rho ({\mathbf {x}},t)$, which can easily be solved numerically. Substituting the expression for J into eq. (4), we find that ${\tilde {S}}\left [\rho \right ] = \frac {1}{2}\int {\mathrm { d}}t\,{\rm d}{\mathbf {x}}\sigma (\rho)\left ({\mathbf {\nabla }}{\hat {\rho }}\right)^{2}$. In practice this is carried out numerically by discretizing the equations, as is described in detail in the appendix. For reference below, we note that on the minimal path $\delta G\left (\rho ,J\right)/\delta \rho = 0$ also holds, which yields an equation of motion for ${\hat {\rho }}$:

$$\begin{equation} \partial_{t}\hat{\rho} = -D(\rho)\nabla^{2}{\hat{\rho}}-\frac{1}{2}\sigma^{\prime}(\rho)\left({\mathbf{\nabla}}{\hat{\rho}}\right)^{2}.\end{equation} \tag{ 7 }$$

Note that the ${\hat {\rho }}$ field is the momentum conjugate to ρ in a Hamiltonian formulation of the problem [19].

In the context of boundary-driven diffusive systems numerical techniques were used to study current large deviations (and similar quantities) [10–13, 22, 23]. Algorithms have also been devised to calculate generating functions in related systems [24]. However, both of these quantities do not yield direct information on the probability density at a specific state. Our algorithm directly minimizes the action, as the works in ref. [14] do. A key feature of our algorithm is the gradual change of the final state, which makes it both stable and insensitive to the choice of optimization algorithm. The algorithm is easy to implement. As a further advantage, the algorithm is easily modified to use only a small number of modes, as explored below.

The numerics were tested against the known 1d and two-dimensional (2d) models for which analytical expressions for large deviation and the trajectory minimizing the action exist. As a first demonstration, we consider the simple symmetric exclusion process (SSEP) in 1d [19, 25, 26]. The model describes a lattice gas with hard-core exclusion, and in the continuum limit leads to D = 1, $\sigma (\rho) = 2\rho \left (1-\rho \right)$ with 0 ⩽ ρ ⩽ 1. We take the domain 0 < x < 1, and the boundary conditions are $\rho _{B}(0),\rho _{B}(1)$. When $\rho _{B}(0)\neq \rho _{B}(1)$, the system is driven out of equilibrium. The most probable state is ${\bar {\rho }}(x) = \left (1-x\right)\rho _{B}(0)+x\rho _{B}(1)$.

Figure 1 shows an example of a path $\rho \left (x,t\right)$ minimizing the action for the SSEP, with a given final state ρ_f. In this case it is known [26] that there is a unique local minimizer for the action, and indeed we find a single solution, independent of the different final-state trajectories tested. Shown is the trajectory from the initial state ${\bar {\rho }}$ to the final state ρ_f at different times, compared with the numerical solution, showing close agreement. The inset shows the contribution to the large deviation (the action integral, eq. (4)), integrated up to time t. The numerics where carried out with N_x = 48 space divisions and N_t = 30 time divisions. In order to capture the time evolution more exactly, the size of the time intervals Δt_i = t_i+1 − t_i was taken to be a geometric series, where the last division is 40 times smaller than the first. The initial time was τ = −0.5. As is clear from the inset, contributions from earlier time are negligible. The relative error in the large deviation is 4·10⁻³. By taking N_x = 150,N_t = 70 the error is reduced to 10⁻³.

Zoom In Zoom Out Reset image size

We now show that the dynamics of the large deviations of slowly varying configurations can be captured by a small number of variables, which describe the long-wavelength behavior. This provides an effective large-scale description of ϕ in terms of a small number of degrees of freedom. To this end we define a family of functions $\left \{\rho _{i}(\mathbf{x})\right\}$ which span the function space, ordered such that $\rho _{1}({\mathbf {x}})$ is the slowest varying in space, followed by $\rho _{2}({\mathbf {x}})$, etc. We then consider an approximation in which the configurations leading to ρ_f are restricted to be linear combinations of a finite number of the slowest-varying ρ_i,

$$\begin{equation} \rho_{M}(x,t) = {\bar{\rho}}(x)+\sum_{i=1}^{M}a_{i}(t)\rho_{i}(x).\end{equation}\noindent \tag{ 8 }$$

For $M\rightarrow \infty$ this recovers the exact extremal solutions. For finite M, when ρ_f is itself of the form of eq. (8) the solutions give upper bounds for the exact value of $\phi \left [\rho _{f}\right ]$, since the minimization is only on a subset all histories $\rho \left (x,t\right)$. Below we are interested in how well they approximate the exact solutions.

A natural choice of the functions $\left \{\rho _{i}(x) \right \}$ are the normal modes of the linearized Hamilton evolution, eqs. (6) and (7), linearized around $\rho = \bar {\rho }$ and ${\hat {\rho }} = 0$ (see footnote ¹). These are obtained by substituting $\rho(x,t) = \bar {\rho}(x) +\sum_{i}\rho _{i}(x)e^{\lambda _{i}t}$, ${\hat {\rho }}(x,t) = \sum _{i}{\hat {\rho}}_{i}(x)e^{\lambda _{i}t}$ and keeping only linear terms in $\rho _{i}(x)$ and $\hat {\rho } _{i}(x)$,

$$\begin{equation} \lambda_{i}\rho_{i} = {\mathbf{\nabla}}\cdot\left(D(\bar{\rho}){\mathbf{\nabla}}\rho_{i}\right) - {\mathbf{\nabla}}\cdot\left(\sigma(\bar{\rho}){\mathbf{\nabla}}{\hat{\rho}}_{i}\right), \end{equation} \tag{ 9a }$$

$$\begin{equation} \lambda_{i}\hat{\rho}_{i} = -D(\bar{\rho}) \nabla^{2}\hat{\rho}_{i}, \end{equation} \tag{ 9b }$$

$\rho _{i} = 0 = {\hat {\rho }}_{i}$

$\rho \left (x,t\right)$

$\bar {\rho }$

${\hat {\rho }}_{i}$

$t\rightarrow -\infty$

$\rho _{f} = {\bar {\rho }}$

${\hat {\rho }}_{i}$

${\hat {\rho }}_{i}$

$\rho _{i}(x)$

$\int \rho _{i}^{2}{\mathrm { d}}x = 1$

As a first example we return to the profile ρ_f given in fig. 1, which is of the form $\rho _{f}=\bar {\rho }+\alpha _{1}\rho _{1}+\alpha _{2}\rho _{2}$, where λ_1,2 are the two lowest eigenvalues, $\lambda _{i}=\left (\pi i\right)^{2}$, and α₁ = −0.256, α₂ = −0.214. The modes ρ_i, for i = 1,2,3 are shown in fig. 2(a). We now minimize the action with histories constrained to be of the form of eq. (8), with M = 2,3,... . Figure 2(b) shows the histories for M = 2 (dashed lines). Even for M = 2, the large deviation is obtained to within 2%, see fig. 2(c), suggesting that even at this level the system's behavior can be captured by an effective model with only two degrees of freedom. Such a small error is striking considering the highly non-linear nature of the problem (as ρ_f is far from ${\bar {\rho }}$), which is generically expected to mix higher modes in significant amounts. By comparison, a local-equilibrium approximation (where a space-dependent chemical potential is set to reproduce $\bar {\rho }$) gives a relative error of 16%, whereas extending the linearized dynamics to the full evolution [23] gives an error in the large deviation of 67% (mostly due to the long-range Gaussian corrections). This highlights both the importance of non-linearities in this problem, and the success of the truncated approximation.

Zoom In Zoom Out Reset image size

The mode approximation can also be used as a high precision numerical method. As shown in fig. 2(c), for M = 15, the relative error is reduced to ∼10⁻⁵, well below the error obtained by a straightforward discretization of space. Indeed similar approaches have been used as numerical tools to improve accuracy in [14].

The low-mode approximation is also useful in cases where an exact analytical solution does not exist. For example, our method allows us to study two-dimensional (2d) systems. In 2d analytic solutions are only known for a set of models which exhibit no long-range correlations. We used one such model, with D = 1 and $\sigma (\rho) = 2\rho$, corresponding to a model of non-interacting particles, as a benchmark for our method, and found agreement between the numerics and the exact solution (see the appendix). In what follows we present results for an interacting system, the SSEP in 2d, for which an exact expression for the large deviation is not known. This exhibits the real power of the numerics. To show the generality of the low-mode approximation, we take a somewhat arbitrary choice of boundary conditions, $\rho \left (x,y\right) = 0.25\,{\mathrm { sin}}(3\,{\rm atan}\,(y/x)) + 0.55 + 0.5y$, on the square domain $A = \left [-1/2,1/2\right ]^{2}$, see fig. 3. The most probable density profile ${\bar {\rho }}$ is shown in fig. 3(a), and we present results for the profile ρ_f shown in fig. 3(b), which, as in the 1d discussion, is of the form $\rho _{f} = \bar {\rho }+\alpha _{1}\rho _{1}+\alpha _{2}\rho _{2}$, where ρ₁,ρ₂ are the lowest modes in this system, and α₁ = 0.097,α₂ = 0.128. Similar results were obtained for other profiles. Figure 3(c) shows the growth of the modes for M = 2,...,10. The first two modes give the exact large deviation to within 2·10⁻³, as estimated by $\left (\phi _{2}\left [\rho _{f}\right ]-\phi _{14}\left [\rho _{f}\right ]\right)/\phi _{14}\left [ \rho _{f}\right ]$, see fig. 3(d). Once more, as in 1d, this means that the evolution is well described by a two-parameter space $(a_{1}(t) ,a_{2}(t))$, despite the non-linear nature of the problem. Interestingly, a local-equilibrium approximation gives a relatively low error of 1.2%.

Zoom In Zoom Out Reset image size

We have shown so far that the large-deviation function is well approximated using only a few modes, provided that ρ_f itself is a slowly varying function of space. Here we discuss how to extend these results to cases where ρ_f is not necessarily slowly varying, but has some high-mode content. Recalling that the bulk behavior is governed by equilibrium dynamics, one might expect that small, high-mode perturbations around the low-mode behavior would be captured by a local-equilibrium theory. In particular, define a local free-energy density $f_{eq,\bar {\rho }}(\rho)\equiv \int _{\bar {\rho }}^{\rho }{\mathrm { d}}\rho ^{\prime }\int _{\bar {\rho }}^{\rho ^{\prime }} \frac {2D\left (\rho ^{\prime \prime }\right)}{\sigma \left (\rho ^{\prime \prime }\right)}{\rm d}\rho ^{\prime \prime }$. In equilibrium, when all boundary densities are equal, this is precisely the free-energy density. We then expect

$$\begin{equation} \phi\left[\rho_{f}\right]-\phi\left[\rho_{f,M}\right]\simeq\int\left[ f_{eq,\bar{\rho}}\left(\rho_{f}\right)-f_{eq,\bar{\rho}}\left(\rho _{f,M}\right)\right]{\rm d}x, \end{equation} \tag{ 10 }$$

where ρ_f,M is ρ_f projected to the subspace spanned by $\{\bar {\rho },\rho _{1},\ldots ,\rho _{M}\}$. In other words, the error due to truncation of the high modes is approximately accounted for by a local-equilibrium theory. We now show that this is indeed the case in an example on the Kipnis-Marchioro-Presutti (KMP) model for heat transfer [16], whose continuum limit [27] gives D = 1 and $\sigma (\rho)=2\rho ^{2}$ (similar results are obtained for other models). In fig. 4, we take ρ_f of the form ${\bar {\rho }} +0.463\rho _{1}+0.507\rho _{2}-0.1\left (\rho _{5}+\rho _{6}+\rho _{7}\right)$. For this profile, the LHS of eq. (10) for M = 2 equals −0.03667 and the RHS equals −0.03623. Thus, eq. (10) is satisfied with relative accuracy of 1.2%. Hence, for profiles with high-mode content the effective low-mode description can be corrected for using eq. (10).

Zoom In Zoom Out Reset image size

In summary, our findings suggest that the evolution leading to a "smooth" rare event is smooth: the continuous diffusion $\partial _{t}\rho ={\mathbf {\nabla }}\cdot \left (D(\rho){\mathbf {\nabla }}\rho \right) + \left \langle {\mathrm{small \ noise}}\right \rangle$, makes an enduring high-frequency perturbation highly unlikely. Indeed we show that even when the evolution leading to the rare event is restricted to profiles involving just a few modes, good quantitative agreement with the exact large deviation may be found. This means that the large deviation is well described in a space involving just a few variables.

Finally, we note that in general the action may have more than a single local minimum, an effect which is well-known in finite-dimensional systems [4–6, 28, 29]. While this does not happen in the well-known models used above in our comparisons, it does in fact exist in other models, and can be studied using the numerical method described here by considering various different final-state trajectories. For example, when two solutions are present, different choices of final-state trajectories will lead to the different local minima². A detailed account of this issue and its physical implications is beyond the scope of the present work, and will be disscussed elsewhere [30].

Appendix. Action evaluation

The action evaluation, including the Lagrange multipliers, is implemented directly in the discrete setting, which gives discrete variants of eq. (6), along with the boundary conditions. In the numerical implementation, time and space are discretized, and ρ is kept at points $\rho _{i,k} = \rho \left (x_{i},t_{k}\right)$ in 1d, and $\rho _{i,j,k}=\rho \left (x_{i},y_{j},t_{k}\right)$ in 2d. We start by describing the method in 2d, and then discuss the simplifications which occur in 1d.

The action, eq. (4), is discretized as $S = \sum _{k}\Delta S_{k}$, where ΔS_k is the value of S associated with the time interval $\left [t_{k},t_{k+1}\right ]$. This allows for the time resolution to vary. For each k separately, ΔS_k is evaluated as

$$\begin{equation} \Delta S_{k} = \Delta t\sum_{i,j}\left(S_{i+1/2,j}^{x}+S_{i,j+1/2}^{y}\right), \end{equation} \tag{ A.1 }$$

where Δt = t_k+1 − t_k. S^x_i+1/2,j corresponds to the bond connecting $\left (i,j\right)$ to $\left (i+1,j\right)$ (and similarly for other half-integer indices).

Let $\rho _{i,j}\equiv \frac {1}{2}\left (\rho _{i,j,k+1}+\rho _{i,j,k}\right)$. Then S^x_i+1/2,j is given by

$$\begin{equation} S_{i+1/2,j}^{x}=\frac{\left(J_{i+1/2,j}^{x}+D\left(\rho_{i+1/2,j}\right) \left(\partial_{x}\rho\right)_{i+1/2,j}\right)^{2}}{2\sigma\left( \rho_{i+1/2,j}\right)} \end{equation} \tag{ A.2 }$$

with $\rho _{i+1/2,j}\equiv \frac {1}{2}\left (\rho _{i,j}+\rho _{i+1,j}\right)$, $\left (\partial _{x}\rho \right)_{i+1/2,j}\equiv (\rho _{i+1,j} - \rho _{i,j})/\Delta x$, and a similar expression for S^y_i,j+1/2. The currents J^x,J^y are constrained to satisfy a discretized version of the continuity equation, eq. (1),

$$\begin{eqnarray} 0 & =&R_{ij}\equiv\left(\partial_{t}\rho\right) _{i,j}\nonumber\\ &&+\frac{J_{i+1/2,j}^{x}-J_{i-1/2,j}^{x}}{\Delta x}+\frac{J_{i,j+1/2}^{y}-J_{i,j-1/2}^{y}}{\Delta y},\qquad \end{eqnarray} \tag{ A.3 }$$

where $\left (\partial _{t}\rho \right)_{i,j}\equiv \left (\rho _{i,j,k+1}-\rho _{i,j,k}\right)/\Delta t$, and Δx,Δy are the (constant) spacings in the x- and y-directions. To minimize the currents subject to eq. (A.3), we define $\Delta G_{k} = \Delta S_{k} -\sum _{ij}\hat {\rho }_{ij}R_{ij}$. Differentiating ΔG_k with respect to the currents gives

$$\begin{eqnarray} J_{i+1/2,j}^{x}& =&-D\left(\rho_{i+1/2,j}\right)\left(\partial_{x}\rho\right)_{i+1/2,j}\nonumber\\ && +\sigma\left(\rho_{i+1/2,j}\right)\displaystyle\frac{\hat{\rho}_{i+1,j}-\hat{\rho}_{i,j}}{\Delta x} \end{eqnarray} \tag{ A.4 }$$

with a similar expression for J^y_i,j+1/2. This is a discrete variant of J = $-D(\rho){\mathbf {\nabla }}\rho +\sigma (\rho){\mathbf {\nabla }}{\hat {\rho }}$. Substituting eq. (A.4) into eq. (A.3), one obtains a linear set of equations for the ${\hat {\rho }}$-variables, which corresponds to eq. (6). These are solved to find the $\hat {\rho }$-variables. Note that on boundary sites eq. (A.3) involves only three currents (or two at the corners of the lattice), which is equivalent to setting ${\hat {\rho }}_{i,j} = 0$ for i,j outside the lattice. This corresponds to the boundary conditions $\left . {\hat {\rho }}\left (x,t\right)\right \vert _{x\in \partial A} = 0$ in the continuum.

Given the ${\hat {\rho }}$ values, the final expression for S is obtained by combining eqs. (A.1), (A.2) and (A.4), and reads

$$\begin{equation}\Delta S_{k}=\frac{1}{2}\Delta t\Delta x\Delta y\sum_{i,j}\left[ \begin{array} [c]{c} \sigma\left(\rho_{i+1/2,j}\right)\left(\frac{\hat{\rho}_{i+1,j} -\hat{\rho}_{i,j}}{\Delta x}\right)^{2}\\ +\sigma\left(\rho_{i,j+1/2}\right)\left(\frac{\hat{\rho}_{i,j+1} -\hat{\rho}_{i,j}}{\Delta y}\right)^{2} \end{array} \right], \end{equation} \tag{ A.5 }$$

which serves as the discrete analog of $\tilde {S}\left [\rho \right ] =\frac {1}{2}\int {\mathrm { d}}t{\rm d}{\mathbf {x}}\sigma (\rho)\left ({\mathbf {\nabla }}{\hat {\rho }}\right)^{2}$. This concludes the evaluation of the action S for a given ρ. This procedure is used as a building block in the optimization algorithm, where S is evaluated for different histories $\rho \left (x,t\right)$, see the main text.

In 1d the above scheme is somewhat simplified. Of course, only terms in the x-direction appear. The continuity equation (A.3) is now $J_{i+1/2}^{x}=J_{i-1/2}^{x}-\Delta x\left (\partial _{t}\rho \right)_{i}$, so $J_{i+1/2}^{x}=J_{c}-\Delta x\sum _{n=1}^{i}\,\left (\partial _{t} \rho \right)_{n}$, where J_c is independent of the position i (but may depend on time). Summing over eq. (A.4), and using $\hat {\rho }_{0} = 0$ we find

$$\begin{eqnarray} &&\hat{\rho}_{i}=\nonumber\\ &&\Delta x\sum_{m=0}^{i-1}\frac{D\left(\rho_{m+1/2}\right) \left(\partial_{x}\rho\right)_{m+1/2}+J_{c}-\Delta x\sum_{n=1} ^{m}\,\left(\partial_{t}\rho\right)_{n}}{\sigma\left(\rho_{m+1/2} \right)},\nonumber\\ \end{eqnarray} \tag{ A.6 }$$

where J_c is fixed by requiring that the boundary condition ${\hat {\rho }}_{N_{x}+1} = 0$ holds.

As an additional tool to improve accuracy, it is possible to interpolate $\rho \left (x,t\right)$ onto a finer grid in $\left (x,t\right)$ before evaluating the action. This simple step improves accuracy and stability at low resolutions. In the example presented below, we use this technique to double the time resolution.

Computations involving modes use exactly the same action evaluation scheme, and only differ in the profiles $\rho \left (x,t\right)$ allowed in the density optimization process.

Appendix. Benchmark for 2d numerical method

The algorithm was tested in 2d against the model σ = 2ρ and D = 1. This model is a particular case of the open-boundary zero-range process [17, 18], and its large deviation is given by

$$\begin{equation} \phi\left[\rho_{f}\right]=\int {\rm d}^{2}x\left(\rho_{f}\ln\frac{\rho_{f} }{\bar{\rho}}+\bar{\rho}-\rho_{f}\right). \end{equation} \tag{ A.7 }$$

$\phi \left [\rho _{f}\right ]$ was calculated for ρ_f in fig. 5(a). Figure 5(b) shows a comparison of the numerical method with the exact result at a relatively low resolution, with N_x = N_y = 8 divisions in each space dimension and N_t = 25 divisions in time, starting from τ = −0.3. The profiles were interpolated onto a grid with twice the time resolution before the action evaluation. The relative error in $\phi \left [\rho _{f}\right ]$ was 5·10⁻⁴.

Zoom In Zoom Out Reset image size