Optimal paths of non-equilibrium stochastic fields: the Kardar-Parisi-Zhang interface as a test case

Atypically large fluctuations in macroscopic non-equilibrium systems continue to attract interest. Their probability can often be determined by the optimal fluctuation method (OFM). The OFM brings about a conditional variational problem, the solution of which describes the"optimal path"of the system which dominates the contribution of different stochastic paths to the desired statistics. The OFM proved efficient in evaluating the probabilities of rare events in a host of systems. However, theoretically predicted optimal paths were observed in stochastic simulations only in diffusive lattice gases, where the predicted optimal density patterns are either stationary, or travel with constant speed. Here we focus on the one-point height distribution of the paradigmatic Kardar-Parisi-Zhang interface. Here the optimal paths, corresponding to the distribution tails at short times, are intrinsically non-stationary and can be predicted analytically. Using the mapping to the directed polymer in a random potential at high temperature, we obtain"snapshots"of the optimal paths in Monte-Carlo simulations which probe the tails with an importance sampling algorithm. For each tail we observe a very narrow"tube"of height profiles around a single optimal path which agrees with the analytical prediction. The agreement holds even at long times, supporting earlier assertions of the validity of the OFM in the tails well beyond the weak-noise limit.

Non-equilibrium behaviors of stochastic macroscopic systems are ubiquitous in nature. One interesting group of questions about such systems concern atypically large fluctuations, which manifest themselves as tails of probability distributions of different fluctuating quantities. Although a universal description of the tails is not likely to emerge, a great progress has been achieved in the last two decades in particular systems. Much of the progress has been due to a method which emerged in different areas of physics under different names: the optimal fluctuation method (OFM), the instanton method, the weak noise theory, the macroscopic fluctuation theory, etc. A key ingredient of this method is a saddle-point evaluation of the path integral of the stochastic process in question, conditioned on the specified large deviation. This evaluation, which employs a model-dependent small parameter (colloquially called "weak noise"), brings about a conditional variational problem. Its solution -a deterministic field, evolving in time -describes the "optimal path": the most probable history of the system which dominates the contribution to the desired statistics. The OFM has been proved efficient in a host of non-equilibrium problems: turbulence and turbulent transport [1][2][3][4][5][6], stochastic surface growth and related systems [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], diffusive transport in the absence  and in the presence [62][63][64][65][66][67][68][69][70][71] of particle reactions or energy dissipation, etc. However, we are aware only of one class of systems -diffusive lattice gases -where theoretically predicted optimal paths were actually observed in stochastic simulations [36,38,42,47,54,55]. Furthermore, in all these systems the optimal density patterns are either stationary, or travel with constant speed.
In the last decade the studies of the KPZ equation and related systems have seen spectacular developments [76][77][78][79]. They shifted from the interface roughness and correlation function to more detailed characteristics such as the full probability distribution of the interface height at a point, P (H, t). Remarkably, exact representations were derived for P (H, t) for the KPZ equation on the line |x| < ∞ for several basic initial conditions: "droplet" [80][81][82][83], flat [84], stationary [85,86], and some of their combinations. At long times and for typical height fluctuations, P (H, t) converges to the Tracy-Widom (TW) distribution for the Gaussian unitary ensemble (GUE) of random matrices [87] for the droplet, to the TW distribution for the Gaussian orthogonal ensemble (GOE) [88] for the flat interface, and to the Baik-Rains (BR) distribution [89] for the stationary interface. These theoretical arXiv:1907.05677v2 [cond-mat.stat-mech] 10 Oct 2019 predictions were confirmed in ingenious experiments with liquid-crystal turbulent fronts [90], see also Ref. [91]. Atypically large fluctuations -the tails of P (H, t) -behave differently and are unrelated to the TW and BR distributions. These tails were determined, for the three basic initial conditions and for several other initial conditions, by the OFM [9-11, 14-16, 18, 20]. The OFM predictions were verified in all cases where the corresponding tails were also extracted from exact representations.
The most thoroughly studied is the droplet case, for which the OFM predicted the following tails [15]: The higher tail (2) coincides with the corresponding tail of the GUE TW distribution, while the lower tail (3) is different from the GUE TW tail (the latter scales as |H| 3 /t). Using the exact representations, it was shown that the tails (2) and (3) are observed both at short times [92] and at long times [93][94][95][96][97][98][99]. Furthermore, these tails have been recently verified in rare-event stochastic simulations [100] which employed a standard mapping to (a discrete version of) the directed polymer in a random potential at high temperature [72,73,75,76,101,102]. The problem of the KPZ one-point height statistics in 1+1 dimension for the droplet initial condition can therefore serve as a good test case for the studies of optimal paths in spatially explicit and non-stationary stochastic systems conditioned on large deviations. Here we use the directed polymer mapping and the importance sampling algorithm [100] to obtain "snapshots" of the optimal paths. As in the previous work [100], our simulations probe the distribution tails with reaching the probability densities as small as 10 −1000 . For each of the two tails we observe a very narrow tube of height profiles around a single optimal path which agrees with an analytical prediction. The agreement holds even at relatively large times, supporting earlier assertions of the validity of the OFM in the tails well beyond the weak-noise limit. 1. Optimal path: theoretical predictions for h(x, t/2). Let us recall the OFM formulation of the problem of one-point KPZ height statistics [9,14,15]. For typical fluctuations, the OFM relies on the smallness of the dimensionless parameter = (t/t NL ) 1/2 1, where t NL = ν 5 /(D 2 λ 4 ) is the characteristic nonlinear time of the KPZ equation. A saddle-point evaluation of the path integral, corresponding to Eq. (1), leads to a variational problem for the action. The resulting Euler-Lagrange equations can be recast in the Hamiltonian form. In the rescaled variables τ /t → τ , x √ νt → x, and λh/ν → h, the Hamilton's equations are where ρ(x, τ ) is the "momentum" density field, canonically conjugate to h(x, τ ); it describes the optimal realization of the noise ξ(x, t). In Ref. [15] a parabolic initial condition h(x, 0) = −x 2 /L was considered. The droplet case follows from these results in the limit of ; the Lagrange multiplier Λ is ultimately expressed through H. Once the OFM problem is solved, one can calculate the (rescaled) action s = s(H) = (1/2) t 0 dτ dx ρ 2 (x, τ ) and obtain, in the original variables, As the noise magnitude D drops out from the "classical mechanics" formulation, the optimal path -the solution of the OFM problem for h(x, τ ) -is independent of D.
As of present, the exact optimal path is unknown, but asymptotic solutions for very large positive and negative H [which determine the tails (2) and (3)] are available [15]. Here are "snapshots" of these asymptotic solutions at see Fig. 1 a. In the leading order in λH/ν 1, h(x, τ ) does not depend on the diffusivity ν. Equation (7) describes a snapshot of two outgoing "ramps" of h(x, τ ), which correspond to a stationary "antishock" of the interface slope v(x, τ ) = ∂ x h(x, τ ) at x = 0, sustained by a very narrow stationary pulse of ρ [8,9,14,15]. The solution (8) describes the noiseless (ξ = 0) evolution of the KPZ equation (1) starting from the droplet initial condition. The solutions (7) and (8) are continuous with their first derivatives at the moving boundaries between them [15]. At τ = t/2 these boundaries are at |x| = λHt/2. If one accounts for the diffusion, the corner singularity of h(x, t/2) at x = 0 is smoothed [14]: For very large negative H the optimal path is quite different. Using the results of Ref. [15], we obtain where the function h(z) is defined parametrically, z(u) = 1 + u arctan u, h(u) = 1 + 2 π u + (u 2 − 1) arctan u , 0 ≤ u < ∞. Here too the optimal path is independent of ν in  the leading order. Noticeable is an extended flat region of h(x, t/2), see Fig. 1 b. As one can see, the optimal paths for large positive and negative H are strikingly different. 2. Directed polymer mapping. Now let us recall the discrete version of the mapping between the KPZ height h(x, t) and the free energy of a directed polymer in a two-dimensional random potential at high temperature T [81]. Consider all directed polymers which start at the point (0, 0) and end at the point (L, L) of a square lattice indexed by integers (i, j) which run from 0 to L, as shown in Fig. 2. The value of the potential V at each lattice point is normally distributed with zero mean and unit variance. Introduce the new variablesτ = i + j andx = (i − j)/2. The partition function Z(i, j) of a given random configuration of the potential, Z(i, j) =Ẑ [(i − j)/2, i + j], with the end point (i, j), obeys the exact recursive equation [81] Z (x,τ + 1) = Ẑ where V (x,τ ) V (x ,τ ) = δxx δττ , and δxx and δττ are Kronecker deltas. Let Z * (x,τ ) = 2 −τẐ (x,τ ). The mapping to the continuous Eq. (1) is achieved via a truncated Taylor expansion of all the discrete quantities and exp(−V /T ) in Eq. (12). For example, To justify these expansions, we use strong inequalities L 1 and T 1. This procedure approximates the discrete Eq. (12) by the stochastic heat equation where we dropped the hats ofx andτ . The (now continuous) noise term ξ is delta-correlated in x and τ . The initial condition is Z * (x, 0) = δ(x). The 3. Importance sampling algorithm. The rare-event simulation method used here is based on the idea to bias the underlying distribution towards the outcomes of interest. It was introduced in computer science under the names importance sampling or variance reduction [103] and has been frequently used in physics, in particular when studying explicitly time-dependent processes with transition-path sampling [104], cloning approaches [105] or density-matrix renormalization group [106]. Our approach is more straightforward and somewhat more general, as it applies to a very broad range of stochastic systems, both in and out of equilibrium [107].
To obtain realizations of the random field V which correspond to extreme values of H = H(V ), we do not sample the random field according to its natural Gaussian product weight G(V ), but according to the modified weight P (V ) ∼ G(V ) exp(−H/Θ). Θ is an auxiliary "temperature" parameter, which allows us to shift the peak of distributions to smaller than (Θ > 0 close to zero) or larger than (Θ < 0 close to zero) typical values. The random realizations cannot be generated directy according to the weight P (V ). Instead we use a standard Markov-chain approach with the Metropolis-Hastings algorithm [108], where the configurations of the Markov chain are the random realizations. By sampling for different values of Θ and combining and normalizing the results, H can be sampled over a larger range of the support [100], down to probabilities like 10 −1000 . For details of the approach, see Refs. [100,107,109]. We extend the approach by storing regularly during the simulation (after equilibration) configurations sampled at various values of Θ. We also extend the simulations to even smaller values of H. For configurations exhibiting corresponding values of H of interest we then can analyze h(x, t/2).
4. Simulation results and comparison with theory. In our first series of simulations we chose L = 2 7 and T 4 = 2 13 . Although = (t/t NL ) 1/2 = 2 is not small, the simulations of Ref. [100] confirmed the validity of the OFM for the whole distribution P(H, t). One can expect, therefore, that the OFM predictions of the optimal paths will be accurate. This is indeed what our simulations show. Figure 3 is a result of processing of 10 simulated configurations at τ = t/2, conditioned on reaching a large positive H (close to 21.9) at τ = t. They appear in the natural ensemble with a probability near 10 −200 [100]. Instead of showing the 10 actual profiles, we only presented the profile averaged over the realizations and the error bars. As one can see, the error bars are strikingly small, implying a very narrow tube of stochastic trajectories around the optimal path. Furthermore, this optimal path agrees very well with the leading-order (zero-diffusion) analytical prediction (7) and (8), and remarkably well with the finite-diffusion expression (9). Figure 4 a shows 10 sampled configurations at τ = t/2, conditioned on a large negative H (close to -31.5) at τ = t. They appear in the Gaussian random ensemble with a probability near 10 −1000 . Here too we only presented the average profile with error bars. Again, the error bars are very small, so that the tube of stochastic paths around the optimal path is very narrow, clearly indicating the presence of a well-defined optimal path. The profile exhibits an extended flat region, as predicted by Eq. (10). The quantitative agreement with the leadingorder zero-diffusion prediction (10) and (11) is fairly good [110]. It should further improve if one solves Eqs. (4) and (5) for this H numerically, rather than rely on the leading-order asymptotic at |H| 1. As we already mentioned, the OFM correctly predicts the leading-order tails (2) and (3) at arbitrary long times [93][94][95]98]. Do the optimal paths still dominate the contribution to the height probability at long times? To answer this question we simulated the directed polymer at much lower temperature, T = 2, when = (t/t NL ) 1/2 = 10 5.5 1. At this temperature the continuous stochastic heat equation (14) is not expected to accurately approximate the exact recursive equation (12). However, we can still expect an emergence of a well-defined optimal path if we condition the discrete polymer on a large deviation of its free energy. This is indeed what our simulations clearly show, see Fig. 4 b. These configurations appear in the Gaussian random ensemble with a probability near 10 −160 . Again, the error bars are very small, implying a well-defined optimal path. The quantitative agreement with Eqs. (10) and (11) is still fair, although h decreases with |x| at large |x| faster than predicted. This latter effect reflects a systematic difference in the behaviors of the deterministic solutions of the discrete equation (12) and of its continuous approximation (14).
5. In conclusion, the problem of one-point height statistics of the KPZ equation clearly demonstrates the versatility and robustness of the OFM. Importantly, the optimal path can persist, and the OFM can be applicable, in the distribution tails, well beyond the weaknoise limit. From a broader perspective, our findings extend previous observations of optimal paths in lowdimensional noisy dynamical system out of equilibrium [111][112][113][114][115] to the realm of non-stationary stochastic fields.
sions. The simulations were performed in Oldenburg on the HPC cluster CARL which is funded by the DFG through its Major Research Instrumentation Programme (INST 184/157-1 FUGG) and the Ministry of Science and Culture (MWK) of the Lower Saxony State. B.M. acknowledges financial support from the Israel Science Foundation (grant No. 807/16). P.S.'s work is supported by the project "High Field Initiative" (CZ.02.1.01/0.0/0.0/15 003/0000449) of the European Regional Development Fund.