Tail of the two-time height distribution for KPZ growth in one dimension

Obtaining the exact multi-time correlations for one-dimensional growth models described by the Kardar-Parisi-Zhang (KPZ) universality class is presently an outstanding open problem. Here, we study the joint probability distribution function (JPDF) of the height of the KPZ equation with droplet initial conditions, at two different times $t_1<t_2$, in the limit where both times are large and their ratio $t_2/t_1$ is fixed. This maps to the JPDF of the free energies of two directed polymers with two different lengths and in the same random potential. Using the replica Bethe ansatz (RBA) method, we obtain the exact tail of the JPDF when one of its argument (the KPZ height at the earlier time $t_1$) is large and positive. Our formula interpolates between two limits where the JPDF decouples: (i) for $t_2/t_1 \to +\infty$ into a product of two GUE Tracy-Widom (TW) distributions, and (ii) for $t_2/t_1 \to 1^+$ into a product of a GUE-TW distribution and a Baik-Rains distribution (associated to stationary KPZ evolution). The lowest cumulants of the height at time $t_2$, conditioned on the one at time $t_1$, are expressed analytically as expansions around these limits, and computed numerically for arbitrary $t_2/t_1$. Moreover we compute the connected two-time correlation, conditioned to a large enough value at $t_1$, providing a quantitative prediction for the so-called persistence of correlations (or ergodicity breaking) in the time evolution from the droplet initial condition. Our RBA results are then compared with arguments based on Airy processes, with satisfactory agreement. These predictions are universal for all models in the KPZ class and should be testable in experiments and numerical simulations.


Introduction
Numerous exact results have been obtained in recent years for discrete growth models in the so-called Kardar-Parisi-Zhang (KPZ) universality class in one dimension. This includes discrete models for the stochastic growth of an interface, such as the PNG model [1,2], but also, more broadly, particle transport models, such as asymmetric exclusion processes (TASEP, ASEP and q-TASEP) [3][4][5][6], directed polymers (DP) in random potentials (last passage percolation [12]), combinatorial problems, and others, where the analogous of a height field evolving under a local stochastic rule can be defined. The KPZ behavior was also observed in non-linear fluctuating hydrodynamics [7,8] and equilibrium and non-equilibrium properties of Bose-Einstein condensates [9][10][11]. At the center of this universality class lies the continuum KPZ equation [13], for which exact results were obtained more recently, for the three main types of initial conditions (droplet, flat, stationary) [14][15][16][17][18][19][20][21][22][23][24][25]. The remarkable aspect of these results was to show how the probability distribution function (PDF) of the fluctuations of the growing height function at one space point, h(X, t), is described at large time by the Tracy-Widom distributions, describing also the statistical behavior of the largest eigenvalue of the Gaussian ensembles of random matrices [26]. In particular for droplet initial conditions the PDF of the height converges to the Tracy-Widom distribution associated to the Gaussian Unitary Ensemble (GUE), while for flat initial conditions it converges to the Tracy-Widom distribution associated to the Gaussian orthogonal ensemble (GOE) of random matrices. Remarkably, the PDF of the height with Brownian initial conditions (stationary class) converges to the so called Baik-Rains distribution [1], which has no equivalent in the random matrices theory. Furthermore it was found that the multi-space point joint PDF (JPDF) identifies at large time [27][28][29][30][31], with the one of a well characterized determinantal process, the Airy processes (Ai 2 /Ai 1 ) [2,32]. These predictions were confirmed in remarkable experiments on liquid crystals [33,34].
All these results and progresses however deal with one time quantities, i.e. observing the KPZ interface at some fixed time t after the initial condition at t = 0. However the two-time observables are also of great interest, for instance the JPDF of the couple of random variables (h(0, t 1 ), h(X, t 2 )) at two very large times t 1 , t 2 , the ratio t 2 /t 1 being fixed, a limit believed to be universal across the KPZ class. Remarkably little is known theoretically about such multi-time correlations, which encode the full time evolution of the KPZ equation. On the other hand the numerical results [35] and experiments are more advanced on the topic and there is a strong need for exact results to be compared with the recent remarkable experimental results [36][37][38].
In order to make substantial theoretical advances on the KPZ equation and, more broadly, on the whole KPZ class, the study of directed polymer models [39] (also called last passage percolation [12]) has often proved useful in the past. Indeed, in the continuum the two models are identical, i.e. the continuum KPZ equation is equivalent, under the Cole-Hopf mapping [40][41][42][43], to a problem of a continuum directed polymer (DP) in a random potential, with h(x, t) = ln Z(x, t) being minus the free-energy and the time t being the polymer length.
In the present context, Johansson recently obtained a rigorous formula for the twotime JPDF of the free energies of the zero-temperature Brownian semi-discrete DP model [44], known to be in the KPZ class. However this formula is quite involved and until now it has not been possible to analyze it analytically or even numerically, in order to produce predictions that can be compared with experiments. Using Airy processes, some results were also recently obtained by Ferrari and Spohn, with an explicit formula valid only for the two-time height correlation (i.e. second moment), and only in the stationary case [45]. Another notable exception is the pioneering work of Dotsenko [46][47][48], who proposed a formula for the two-time JPDF of the free energies (i.e. of the KPZ heights). However we believe that this formula is not correct and we briefly explain our reasons in the text.
For the continuum directed polymer/KPZ equation, the replica Bethe ansatz method allows to write aesthetically appealing, formal equations for the joint integer moments of the partition sum of the directed polymer Z(X, t) = e h(X,t) at fixed temperature and for two different and finite lengths. It involves summations over the known eigenstates of the attractive one-dimensional delta-Bose gas, the so-called Lieb-Liniger model [49]. Hence one could attempt a summation over these states, to obtain an explicit formula for these moments, and from them, extract the two-time JPDF. This program was indeed successfully carried through for the one-time problems [14,15,19].
Unfortunately, when more than one time is involved, the summand contains the socalled form factors of the bosonic creation and annihilation operators. Those form factors are notoriously difficult to calculate in full generality [50,51]. In Refs. [46][47][48] the above described summation was indeed attempted with a guess for the form factors which, in our opinion, is not correct, although it leads to a simplified summation. As will be confirmed here, many terms of the exact sum are thus missing, leading then to an incomplete result for the JPDF. Although one could hope that that these terms are negligible in the large time limit, we show here explicitly that this is not the case. Hence the problem is still mostly unsolved.
The aim of the present paper is to study the joint probability distribution function (JPDF) of the two heights (h(0, t 1 ), h(X, t 2 )) of the Kardar-Parisi-Zhang equation at two different times t 1 < t 2 and with droplet initial conditions. We do not obtain a full solution of the problem, namely the full JPDF, which still remains an open problem. However, we are able to obtain an exact expression for the JPDF in the limit where h(0, t 1 ) (i.e. its fluctuations) is positive and large, i.e. the exact tail of the JPDF when one of the two arguments is large. The dependence in the second argument, describing the fluctuations of h(X, t 2 ), remains non-trivial. Although we do also obtain an expression valid for finite time, we focus here on the limit where both times are large while their ratio t 2 /t 1 is fixed. As this ratio t 2 /t 1 is varied, our formula for the JPDF exhibits a crossover between two limits where it decouples into a product of two distributions, each of them describing the fluctuations at one time. In the limit where the two times become equal t 2 /t 1 → 1 + we recover the product of a GUE Tracy-Widom PDF and a Baik-Rains distribution. By independent arguments, based on Airy processes, we can show that this decoupling should hold exactly in that limit (beyond the tail), which provides an important check of our method. The opposite limit, where the ratio t 2 /t 1 is very large, leads to a factorized JPDF, into a product of two GUE Tracy-Widom distributions. The subleading corrections decay in time as (t 1 /t 2 ) 1/3 , which leads to a non-trivial two-time correlations also in the limit where t 2 t 1 , a phenomenon known as persistence of correlations [34,37,45,52] and that we show here show from our exact expression for the tail of the JPDF. We obtain exact expressions for the cumulants of h(X, t 2 ) conditioned to a fixed large value of h(0, t 1 ), and their expansion in the limit of large t 1 /t 2 . Finally, our predictions are universal for all models in the KPZ class and should be testable in experiments and numerical simulations.
As is detailed below, our method is based on a partial summation over the Bethe ansatz states: we include only a single string state (bound state of all the particles describing the different replicas) propagating in the first time interval [0, t 1 ], which can be written explicitly, and, as we argue, gives the exact tail of the JPDF when h(0, t 1 ) is large. The resulting non-trivial summation over the other many-string states related to the second variable can then be performed exactly at arbitrary finite times and with no approximation.
The outline of the paper is as follows. In Section 2 we review the known results for the one-time PDF of the height of the KPZ equation, introduce units and notations, and define the tail approximation. In Section 3 we define the JPDF that we aim to calculate, and we present the main results of the paper, some in explicit form, and for others we refer to formula later in the text. In Section 4 we present the method of the replica Bethe ansatz (RBA) for the two-time problem. We recall the quantum mechanical approach based on the attractive delta-Bose gas, i.e. the Lieb Liniger (LL) model. We use it to express the two-time joint moments of the DP partition sums as sums over two sets of eigenstates of the LL model, the so-called strings. We define and present a formula for the generating function. We define the tail approximation used in this paper, were the summation is restricted to 1-string states for the first set, and unrestricted for the second. In Section 5 we derive the expression of the 1-string form factor, namely the form factor of a generic power of the annihilation operator between a state containing only one single string (bound state) and a generic Bethe state. The explicit formula allows to perform the summation explicitly in Section 6. After a number of manipulations, we obtain our main result for the tail generating function, and the JPDF. We analyze this JPDF in three limits, and perform expansions around these limits (the t 2 /t 1 → 1 limit, the large t 2 /t 1 limit and the large h 1 limit). In Appendix 7 we analyze the problem using Airy processes, and compare with the results of the RBA method. Finally Section 8 contains the Conclusion.
Further technical details are given in the Appendices. In Appendix A we detail the comparison between our work and the papers [46][47][48]. In Appendix B we present the calculation of the form factor used in the text. In Appendix C we present an independent calculation of the "double tail" in both height variables. In Appendix D we show a more general result for the t 2 /t 1 → 1 + limit. In Appendix E we give more details about extended Baik-Rains distributions and derive some new expansion formula. In Appendix F we explain a contour integral method for constrained summations needed in the text. In Appendix G we give useful identities, such as expressions of derivatives of the GUE-TW distribution. Finally, in Appendix H we calculate higher orders in the large t 2 /t 1 expansion.

Units and the Cole-Hopf mapping
In this paper we study the one-dimensional Kardar-Parisi-Zhang (KPZ) equation [13]. It describes, in the continuum, the stochastic growth of an interface, of height h(x, t) at point x ∈ R, as a function of time t driven by a unit white noise ξ(x, t)ξ(x , t ) = δ(x − x )δ(t − t ). We will use the following scales as units i.e. we set x → x * , t → t * and h → h * so that from now on x, t, h are in dimensionless units and the KPZ equation becomes Eq. (1) with ν = 1, λ 0 = 2 and D = 2. The Cole-Hopf mapping solves the KPZ equation from an arbitrary initial condition h(x, 0) as follows: the solution at time t can be written as: Here Z η (x, t|y, 0) is the partition function of the continuum directed polymer in the random potential − √ 2 η(x, t) with fixed endpoints at (x, t) and (y, 0): which solves the (multiplicative) stochastic heat equation (SHE): with Ito convention and initial condition Z η (x = 0, t|y, 0) = δ(x − y). Equivalently, Z(x, t) is the solution of (5) with initial conditions Z(x, t = 0) = e h(x,t=0) . We will sometimes omit the "environment" index η. Here and below overbars denote averages over the white noise η.

Large time results
Let us briefly review the known results on the one-time fluctuations of the KPZ field. We focus here on the droplet initial condition. To define it for arbitrary time one usually considers the wedge initial condition defined as and the "hard wedge" limit w 0 → +∞, so that exp(h w 0 (x, 0)) → δ(x). In that limit h w 0 (x, t) → h(x, t) for any t > 0, and one has where Z η (x, t|0, 0) is the partition function of the DP defined in the previous Section ‡. Hence the droplet height profile h(x, t) corresponds to the free energy of a DP in a random medium going from the point (0, 0) to (x, t) (see Fig. 2). ‡ Note that if one is interested only in the large time limit, any finite w 0 > 0 leads to the same result, up to a constant shift.  Figure 1. Plot of the Tracy-Widom distribution F 2 (σ) (red line) compared with its tail for positive σ given by F 2 (σ) (black line). Inset: difference between the Tracy-Widom distribution F 2 (σ) and its tail F At large time the KPZ field grows linearly in time with O(t 1/3 ) fluctuations. At the level of the fluctuations in one space point (choosing here x = 0) these are governed by the GUE Tracy Widom of PDF, f 2 (σ) = F 2 (σ), i.e. one has at large t [2,14] h where F 2 (σ) is given by a Fredholm determinant involving the Airy Kernel K Ai : and where P σ (v) = θ(v − σ) is the projector on [σ, +∞[. To get rid of the part linear in time we will, from now on, redefine the KPZ field, and the DP partition sum, at all times, as and for notational simplicity, we will omit the tilde in what follow.

Tail approximation
Since we obtain below some exact results for the tail of the two-time JPDF, it is useful to introduce here the tail of the single time GUE-TW distribution. We define the tail approximations of the CDF of the GUE-TW distribution as the function 3 σ 3/2 ) for σ → +∞ and the corrections are of higher (stretched) exponential order As can be seen on Fig. 1, this approximation is reasonably good (with error less than 10 −3 ) for any σ > −1. More generally below we will always use the superscript (1) to indicate the tail approximation (13). It is important to note that keeping the full trace in (12) is a considerably better approximation that its leading large σ limit. As discussed below, this approximation also consists in keeping only single string states in the RBA method, leading to huge simplifications in the application of the method.
3. Two-time problem: summary of the results

Notations
Schematic representation of the two-time setting: we consider the joint probability distribution JPDF of the free energy of a DP in a random medium starting at (0, 0) and ending at (X 1 , t 1 ) and another polymer in the same random medium starting at the same position and ending at (X 2 , t 2 ). In the following, we focus on the case X 1 = 0 and X 2 = X.
Here we study the KPZ height field h(x, t) with droplet initial condition at two different times t = t 1 and t = t 2 > t 1 , and two different space points x = 0 and x = X, and we denote In particular we will be also interested in the difference of the two height functions built over the time difference t 21 ≡ t 2 − t 1 . In the large time limit, namely when both t 1 and t 2 are sent to +∞, the relevant parameter characterizing the JPDF of the two heights will be the time difference rescaled by the earlier time, denoted as Let us start with X = 0, the general case being discussed below. From the previous Section, we know that the two heights grow in time as H 1 ∼ t Here h 1 , h 2 , h are random variables whose joint distributions we want to determine. From the previous Section we know that at large time each individual height is distributed according to the GUE-TW distribution We define the joint probability distribution function (JPDF) at two arbitrary finite times 0 < t 1 < t 2 of the rescaled KPZ field h 1 and height difference h as where σ 1 , σ are two real numbers and λ parametrizes the time t 1 as λ = (t 1 /4) 1/3 . We will be particularly interested in the limit of all times being large keeping ∆ finite, such that t 1 ∆ ∼ t 1 , which we denote as

The two-time joint PDF
As mentioned in the introduction the full JPDF, P ∆ (σ 1 , σ), is very difficult to compute exactly, and here we obtain this function only in the region of large positive σ 1 , and for arbitrary σ. More precisely we find that where P (1) We are able to derive a relatively simple expression for the tail approximation P (1) ∆ (σ 1 , σ). We report the result as follows: expressed in terms of the Airy kernel, as well as of a novel kernel where we recall that F 2 (σ) and F 2 (σ 1 ) are respectively the GUE-TW CDF and its tail approximation, given respectively by (9) and (12). Below, the function f ∆ (σ 1 , σ) in (23) is traded for the function calledĝ Our exact result for the JPDF satisfies two important properties. In the limit of infinite time difference t 2 /t 1 → +∞ (which corresponds to the limit ∆ → +∞), it converges to the product of two GUE-TW distributions (25) and in the limit of small (scaled) time separation (t 2 − t 1 )/t 1 1 (which translates into ∆ → 0 + ) it also decouples as follows where F 0 (σ) is the Baik-Rains (BR) probability distribution [1,24] which governs the stationary growth profile in the infinite time limit, and whose explicit expression is recalled in (157). In fact, reexamining the problem in terms of Airy processes in Section 7, we argue that in these two limits our result is valid for any value of σ 1 , i.e. we predict the exact property of the full JPDF hence our exact result is able to reproduce this property. Note that the second property has not been found previously. It is non-trivial since here the large time limits t 1 , t 2 → ∞ have been taken. A systematic expansion in the large ∆ limit is performed in Section 6.6, and the results are obtained in an expansion in 1/∆ 1/3 . The small ∆ limit is analyzed in Section 6.5 where (26) is shown, and the leading correction (in ∆ 1/3 ) are obtained.
All the above results can be extended to an arbitrary endpoint X, in terms of the scaling variableX In such a case the definition of the scaled height difference variable h becomes where we recall ∆t 1 = t 2 − t 1 . Then definition (20) and formula (23) and (23) still hold, with a more general expression for the kernel K ∆ σ 1 , given in (140). The property (25) also holds. In the property (26), F 0 (σ) is replaced by F 0 (σ −X 2 ;X) ≡ H(σ;X 2 , −X 2 ), the extended Baik-Rains distribution (i.e. its CDF). The function H is given in [1]. It describes the one point PDF of the (scaled) height of the KPZ equation with tilted stationary initial conditions, in the large time limit, as recalled in (208) and in (288).
3.3. Conditional probability and moments: numerical evaluations, and series expansions 3.3.1. Conditional probability and its moments. In addition to the JPDF, of great interest is the conditional PDF (CPDF) for h = σ given that h 1 = σ 1 , together with its moments From what we discussed above as ∆ increases from zero to infinity it is predicted to crossover from the Baik-Rains distribution F 0 to the GUE-TW distribution F 2 (and their moments). It is useful to recall the known values of the lowest cumulants, as well as the skewness and kurtosis of the two limiting cases as well as [58] h BR = 0 , h 2 c BR = 1.15039 At the two limit points ∆ = 0 and ∆ = +∞ the CPDF and the conditional moments (cumulants) are independent of σ 1 since the JPDF factorizes. For intermediate values of ∆ these depend on σ 1 and the result that we obtain here is their value in the limit of large positive σ 1 . More precisely we calculate which deviate from their exact values with corrections that are exponentially suppressed in σ 1 at large positive σ 1 In particular we are able to obtain analytical expressions for the moments of the conditional probability in the limit ∆ → 0 (here given forX = 0) where the first three coefficients {A p } 3 p=1 are calculated as {−0.5751, −0.2214, −2.1755} and the general formula is given in (160).
In the limit ∆ → ∞ they converge to the moments of the GUE-TW distribution, and we have obtained the corrections as with where the next two orders in the ∆ −1/3 expansion are reported and analyzed, together with the cumulants, in (171) and the equations below it. Note that it is shown there that these coefficients are independent ofX. Indeed the dependence of the cumulants on the endpointX vanishes exponentially in ∆, see Fig. 6, therefore the coefficients of the large ∆ expansion of the first cumulants are all independent of the final endpoint. As discussed below this is a manifestation of the persistence of correlations and ergodicity breaking during the time evolution with droplet initial condition. In order to have an easier comparison with numerics or experiments it may be more advantageous to define a cumulative conditioning, i.e. (42) for which more statistics can be accumulated, making it easier to compare with numerical or experimental results, see Fig. 5. On the other hand the integrated cumulants have properties very similar to the conditional cumulants.
Note that an important sum rule satisfied by the exact JPDF is Similarly, our result for the tail function of the JPDF, P ∆ , satisfies the sum rule (shown in I) Therefore we can use the difference F 2 (σ c ) − F 2 (σ c ) as a good measure of the error on the cumulative distribution (42) due the approximating the JPDF with its tail for large σ 1 . Encouraged by Fig. 1, we surmise that the predictions remain quite precise for a quite broad range of values of σ c and especially in the domain σ c > 0.

Unconditioned mean value
While here we obtain results for the conditional first moment of h, it is useful to point out that the unconditionned first moment can be calculated exactly for arbitrary fixed ∆ and t 1 → +∞, from pure scaling. Indeed, using that at large time one immediately obtains from the definitions (18) that where the second form is adapted to obtaining readily the series expansion at large ∆. Note that this behavior in 1/∆ holds equally well for a larger class of initial conditions, e.g. for flat or stationary initial condition (replacing σ GUE → σ GOE and σ GUE → σ BR respectively). Formula (47) immediately leads to the following asymptotic behaviors: Appendix 7 for Airy processes arguments). Interestingly, similar behavior for the unconditioned cumulants at small and large ∆ was found numerically and experimentally, not for the droplet geometry, but for the flat initial condition in [36], although no determination of the prefactors (predicted by the present arguments) was discussed there.
3.3.3. Persistence of correlations and conditional two-time correlations The first order correction O(∆ −1/3 ) of the conditional moments to their respective GUE values, i.e. the term ∆ −1/3 R 1/3 (σ 1 ) in (40), is at the origin of the property of "persistence of correlations" for droplet initial conditions. This property, also called breaking of ergodicity in the temporal evolution, states that the two-time connected correlation of the KPZ heights H 1 , H 2 at two large and well separated times t 1 , t 2 do not decay to zero in the limit t 2 /t 1 → +∞ [34,37,45,52], i.e. the following dimensionless connected correlation has a finite limit In the present work we define a conditional version of the above function involving correlations restricted to realizations such that h 1 > σ c . Our main result for the JPDF (23) allows to numerically evaluate this function for general ∆, see Fig. 8. In Section 6.7 we furthermore obtain an analytic expression for its large ∆ limit, c σc , which is exact for large positive σ c and expected to be accurate up to σ c ≈ −1.5 (see Fig. 12 there). An extrapolation to σ c → −∞ provides an estimate for c = lim σc→−∞ c σc ≈ 0.58 ± 0.05 close to the observation in Fig. 12 of [34]. Note that our result provides a proof by explicit calculation that the large ∆ limit of this ratio is finite. A complementary study is made using Airy processes in Section 7 with satisfactory agreement. Both methods show that there is no dependence of this ratio in the endpoint positionX. Note that since this property does not occur for flat initial conditions [34,37,45,52] we surmise that in that case h σ 1 flat = h GOE + o(∆ −1/3 ). Fig. 3 we report the results for the conditional probability distribution function obtained by numerical evaluations of equation (22) as a function of ∆ and at fixed σ 1 . In Fig. 4) we plot the lowest conditional cumulants of the scaled height difference h = σ, i.e. the mean value, standard deviation, skewness and kurtosis, for a fixed value of h 1 = σ 1 (in particular σ 1 = 0, 1) where the approximation (13) is believed to be correct (see Fig. 4) up to deviation of 10 −4 (an estimate based on the approximation of the GUE probability distribution with its tail, see Fig.  1). We note the genuine feature that the variance is slightly non-monotonous, even if the deviations from a monotonic behavior are very small (at least for the values of σ 1 considered here). On the other hand the other cumulants are monotonous decreasing functions of ∆. In Fig. 6 we study the effect of the endpoint displacementX on the 1 Figure 4. Plot of the lowest cumulants (mean, variance, skewness and kurtosis) of the conditional probability distribution P (1) (σ|σ 1 ), defined in (37), at fixed σ 1 = 0 (red lines) and σ 1 = 1 (black lines) as a function of ∆ 1/3 . The small horizontal lines on the right in each plot indicate the asymptotic values for infinite ∆, namely the GUE expectation values reported in (32). Note the weak non-monotonicity of the variance which is a genuine effect.

Numerical evaluations In
first two cumulants as a function of ∆. Since at ∆ = 0 we recover the BR distribution with an arbitrary endpointX (see Appendix E) we also plot here for completeness its first two cumulants as a function ofX in Fig. 7.
The numerical absolute error on the cumulants is estimated (for both cases of conditional and integrated probability) to be around 10 −6 for the mean value, 10 −5 for the variance, 10 −4 for the skewness and 10 −3 for the kurtosis. This decrease of precision for higher cumulants is due to the finite size effects induced by the discretization of the kernels and which are larger along the tails of the distribution. For more details on the numerical evaluations see Appendix J.  Figure 5. Plot of the lowest cumulants (mean, variance, skewness and kurtosis) of the integrated conditional probability distribution P (1) (σ|h 1 > σ c ), defined in (42), at fixed σ c = −1 (red lines) and σ c = 0 (black lines) as a function of ∆ 1/3 . The small horizontal lines on the right in each plot indicate the asymptotic values for infinite ∆, namely the GUE expectation values reported in (32). The small deviation from its asymptotic value of the kurtosis at σ c = −1 is due to numerical error.

Replica
Bethe ansatz approach to the two-time correlations

Joint integer moments from the δ−Bose gas
As for the calculation of the one-time observables, the RBA method starts with the calculation of the integer moments of the DP partition sums, from which one aims to extract the joint PDF of KPZ height fields. Here we consider the joint moments with n 1 , n 2 ≥ 0, as defined in (14). Here we have allowed for arbitrary endpoints X 1 , X 2 , however thanks to the so-called statistical tilt symmetry, one can choose with no loss of generality, X 1 = 0 and X 2 = X, which we will do later on § § As compared to the result for (X 1 , X 2 ) = (0, X), a general X 1 only results in a simultaneous deterministic shift of H 1 by −X 2 1 /(4t 1 ) and of H 2 by (X 2 − (X − X 1 ) 2 )/(4(t 2 − t 1 )).  Figure 6. Plot of the mean and variance of the integrated conditional probability distribution P (1) (σ|h 1 > σ c ) at fixed σ c = 0 and endpoint position (29)X = 0 (black lines),X = 0.5 (red lines),X = 1 (blue lines) as a function of ∆ 1/3 . At ∆ = 0 the cumulants converge to the mean and variance of the Baik-Rains distribution for with a displacementX (see Fig. 7). In the other limit, for large ∆, the dependence onX is exponentially (as ∼ e −const ∆ 1/3 ) vanishing.  For calculational convenience, we will consider a slightly more general partition sum Z 2 defined as with w > 0. From the definition of the partition sums as path integrals (4) it is easy to see that in the limit w → 0 one recovers  (40). For any value of σ c the data approach their asymptotic values with corrections of order ∆ −2/3 (see Section 6.7 for details).
hence one obtains (52). We will perform the calculation with w > 0 (implicitly, for notational simplicity) and perform the limit w → 0 at the end. It turns out that all the n−th joint moments of the DP can be mapped to quantum expectation values in a model of n attractive bosons on the line [53,54]. This is known as the attractive Lieb-Liniger model [49] for n bosons with Hamiltonian on the line x ∈ R. In our chosen unitc = 1, which we will use from now on. For definiteness one fixes periodic boundary conditions with period L, where L → +∞ is considered afterwards. The joints moments of the DP can then be expressed as the matrix elements of the quantum evolution operator with imaginary time t Using these relations, and the statistical independence of the partition sum (i.e. of the noise) on the time intervals [0, t 1 ] and [t 1 , t 2 ], we can rewrite the joint two-time moments as with t 2 − t 1 = ∆t 1 (see Fig. 9). Here we have assumed n 1 ≥ 1 and n 2 ≥ 1. Clearly Z n 1 ,0 (t 1 , t 2 ) = Z n 1 1 = Z n 1 (t 1 ) and Z 0,n 2 (t 1 , t 2 ) = Z n 2 2 = Z n 2 (t 2 ) reduce to a one-time problem Z n (t), which was studied in [14] (for w = 0) and which we do not recall in details here. Let us now use the decomposition of the evolution operator into the eigenstates of H n , namely e −tHn = µ∈Λn |µ e −tEµ µ|. Here Λ n denotes a complete orthogonal basis of eigenstates |µ of H n , E µ their eigenenergy and ||µ|| their norms. The wave-functions in the coordinate basis are denoted as ψ µ (x 1 , ..x n ) ≡ x 1 , ..x n |µ . We obtain where F n 2 ;n 1 +n 2 µ;γ is a "form factor" defined as For w = 0, it is proportional to the form factor of the operator [Ψ(X)] n 1 in the secondquantized formulation (which we will not use in our work). Until now the formula are very general, they hold for arbitrary L, and even apply to a broader class of time uncorrelated noise η(x, t) upon a suitable generalization of H n .

Bethe ansatz in large size limit and strings
We now use the known properties of the LL model [49], which allows to treat the continuum KPZ equation when the noise η(x, t) is a Gaussian white noise. The eigenstates are parametrized by a set of distinct (in general complex) quasi-momenta or rapidities µ ≡ {λ 1 , ..λ n }. These are determined by the coupled Bethe equations, which are the quantization rules for the rapidities and that enforce the periodicity of the eigenstate |{λ 1 , ..λ n } . For |µ ∈ Λ n , the (un-normalized) eigenfunctions are totally symmetric in the x α , and are such that, on the sector x 1 ≤ x 2 ≤ .. ≤ x n , they are equal to (forc = 1) where the sum runs over all n! permutations P of the rapidities λ j . The eigenenergies associated to (60) are given by a sum over all the single particle energies E µ = n α=1 λ 2 α . Note that the eigenfunction (60) computed with all the particles sitting at the same point is given by ψ µ (x, ..x) = n!e ix α λα . Focusing on the limit L → +∞, the general eigenstate is built by partitioning the n particles into a set of n s ≤ n bound states called strings [55] formed by m j ≥ 1 particles with n = ns j=1 m j . Their associated rapidities are where m j k j ∈ R is the total momentum of string j. Here, a j = 1, ..., m j labels the rapidities within the string j = 1, . . . n s . From now on a string state is denoted as |µ = |k, m and labeled by the set of (k j , m j ) j=1,..ns . Inserting (61) in (60) leads to the Bethe eigenstates of the infinite system with corresponding eigenvalues: Note that these are the eigenvalues of the shifted Hamiltonian H n − 1 12 n, which is the proper Hamiltonian to consider, in accordance with the redefinition (11). The (inverse squared) norm of these eigenstates are given by [56] 1 In the large size limit L → ∞ the string momenta k j follow the same quantization rules as the ones for free particles k j = 2πI j L (where I j are integer numbers), therefore the sums over the eigenstates in (58) can be recast into integrals over the whole real line, one for each string momenta. In particular we can write lim L→+∞ µ∈Λn We are now finally in position to write the final formula for the joint moments (58) in the large L limit. The sums in (58) are now involving string states |γ = |q, m γ , for the time interval [0, t 1 ], and |µ = |p, m µ for the time interval [t 1 , t 2 ]. Using (64) for each sum, it leads to where F n 2 ;n 1 +n 2 p,m µ ;q,m γ is the form factor given in (59) specialized to the case X 1 = 0 and evaluated on generic strings states. The explicit dependence on X 1 , X 2 has been extracted, using that for any Bethe state ψ µ (x 1 + X, .., x n + X) = e iX n α=1 λα . Apart from the form factors all the other factors are explicit functions.

Two-time generating function: definitions
In order to obtain the probability distribution from the moments, we now introduce a two-time generating function. To avoid proliferation of factors 2 1/3 in the calculation, we remind that we parametrize the time t 1 as and, from now on, we use interchangeably the two variables σ and s (ors, in the case X 2 = X = 0, X 1 = 0), related via so that λs 1 = t 1/3 1 σ 1 and λs = t 1/3 1 σ. We define the generating function of the two-time problem as It can be rewritten in terms of the (minus) the scaled DP free energies, i.e. the scaled KPZ heights h 1 and h 21 ≡ h defined in (14),(213), (18), and when X = 0 in (30), leading to where we have definedh 1 = 2 2/3 h 1 andh = 2 2/3 (h −X 2 ). If this function is known, then one can retrieve in principle the JPDF of h 1 and h, following e.g. [14]. For now let us indicate how to do it only in the limit of both times being large, i.e. λ large, at fixed ∆, which is our main focus here.
In that case the two double exponentials in (69) converge to the Heaviside step function θ(x) in the limit λ → ∞ leading to It is then easy to see that the JPDF of h 1 and h, i.e. P ∆ (h 1 , h) defined in (20)- (21), is obtained by taking derivatives respect to σ 1 , σ as follows

Generating function: expression as a series in the number of strings
We can now expand the exponential in (68) and write the following expansion of the generating function in terms of the two-time joint moments: where we have separated the pieces which can be expressed in terms of g λ (s; X), the generating function of the one-time problem (i.e. point to point DP problem, droplet IC for KPZ), namely for n 1 = 0 or n 2 = 0, and the piece which genuinely depends on the two timeŝ Note that we have written the second term in (72) directly in the limit w = 0. The analysis of this term for w > 0, and a generalization thereof, is slighlty more involved and is performed in the Appendix. Note that in the limit of both times being large one thus obtains (setting X 1 = X 2 = 0 in that equation for simplicity) where, we recall F 2 (s) is the CDF of the GUE-TW distribution andĝ ∆ (s 1 , s) = lim λ→∞ĝ∆,λ (s 1 , s). We note that the one-time terms, namely the ones that depend only on one single rescaled height at either t 1 , i.e. only on σ 1 (from H 1 ), or at t 2 , i.e. only on σ 1 + ∆ 1/3 σ (from H 2 ), do not produce any contribution to the JPDF, hence (71) leads to We can now use (65) to express the generating function for arbitrary finite times as a series of terms with different values of the numbers of strings n µ s , n γ s aŝ 4.5. Generating function: tail approximation at large h 1 and single string contribution In principle the formula (77), (78), if properly evaluated, should lead to the full result for the two-time JPDF via (76). This requires however to obtain closed expressions for all the form factors. While such a program is in progress, it is not yet available. In particular, the recent attempt in [46], is, as we will show, incorrect. Our aim will thus be more modest. We will focus on the tail approximation of the JPDF for large h 1 , but arbitrary h. To this aim, instead of calculating the complete sum (77), we will obtain an exact expression for the partial sum i.e. restrict to states such that the initial string state (in the interval [0, t 1 ]) is made of only one string, n γ s = 1, of length m γ 1 = n 1 + n 2 , as shown in Fig. 10. The main reason of this restriction is that this is a case where the form factor can be explicitly calculated, as shown in the next Section.
We will then be able to compute explicitly in the limit of both times large the partial sumĝ (1) ∆,λ (s 1 , s) and we will argue that it contains the complete leading large s 1 correction, witĥ 1 ) uniformly in s. Hence, via (76), we will obtain the leading tail of the JPDF, P ∆,λ (s 1 , s), as described in (22) , with It is interesting to note that there is also a double tail approximation where both h 1 and h are large (i.e. both σ and σ 1 are large): such that is the leading term in the tail of the JPDF in both variables. The calculation of this particular term will be performed independently in Appendix C. Figure 10. Graphical representation of the one-string approximation. We restrict the sum over all the eigenstates |γ in the sector t ≤ t 1 to a reduced sum over all the states such that all the n 1 + n 2 bosons are in the same bound state (string). On the other hand we give no restrictions for the eigenstates in the time interval t > t 1 , where we sum over all the possible configurations with n 2 bosons for the state |µ . The one-string to multiple-string form factor F n2;n1+n2 µ;γ (85) connects the two states at t = t 1 in a non-trivial manner. The summations can now be performed exactly.

A one-string to multiple-string form factor in the δ−Bose Gas
Consider now the form factor F n 2 ;n 1 +n 2 µ;γ = F n 2 ;n 1 +n 2 p,m µ ;q,m γ introduced in (65), and evaluated on two string states such that |γ contains a single string, i.e n γ s = 1 and |γ = |q, n 1 +n 2 while the second string state |µ = |p, m is an arbitrary string state with n µ s ≥ 1. Clearly the strings in the state |µ must satisfy the constraint n µ s j=1 m j = n 2 . It is shown in Appendix B that this form factor admits a relatively simple expression, namely F n 2 ;n 1 +n 2 p,m µ ;q,n 1 +n 2 = n 2 !(n 1 + n 2 )!(2w + n 1 ) n 2 ns j=1 Its calculation turns out to be very similar to the one of the overlap in the case of KPZ with a stationary initial condition for KPZ. Interestingly enough, as discussed later, the stationary KPZ will also play an important role in the physics of the present problem.
6. Calculation of the tail of the two-time JPDF

Aim and starting formula
Here our aim is to calculate the functionĝ ∆,λ (s 1 , s) defined in (79). As discussed in Section 4.5 it is a partial sum over states, which contains only contribution of singlestring states |γ = |q, n 1 + n 2 , i.e. n γ s = 1, with arbitrary string states |µ and we will denote here for simplicity n µ s ≡ n s and m µ j ≡ m j . We first obtain expressions valid for arbitrary λ, then we study large λ, i.e. calculateĝ (1) ∆ (s 1 , s), which contains the desired result about the tail of the JPDF via (81) . We recall that it is given as a series in increasing number of strings n s aŝ where from now on we use the shorthand notationẐ n γ s =1,n µ s (s 1 , s) ≡Ẑ 1,ns (s 1 , s). Using (78) and the expression for the form factor (85) we obtain our starting formula for the partition sums at fixed n s aŝ where we have set X 1 = 0 and X 2 = X with no loss of generality. We have defined the scaled positionX replaced t 1 = 4λ 3 everywhere, and then performed the rescalings p j → p j /(2∆ 1/3 λ), q → q/(2λ), convenient for later use. The challenge is now to simplifyẐ 1,ns (s 1 , s) so as to be able to perform the sum. This will be done below in the limit of large time λ → +∞. Before doing so let us point out that the single string termẐ 1,1 (s 1 , s) is somewhat easier to evaluate, and its detailed calculation is given in C. It will be checked that it can be recovered from the full calculation in the limit of large and positive s. In addition, as discussed above and below, it contains information about the double tail of the JPDF.

Calculation ofẐ 1ns (s 1 , s)
One of the main difficulty in performing the sums in (87) is that, due to the Kronecker delta, δ n 2 , ns j=1 m j , the sum over the m j is constrained and not independent from the sum over n 2 . Alternatively replacing n 2 by ns j=1 m j in the expression for the form factor lead to a complicated coupling of the different strings, making then impossible to sum over all the possible string lengths. To perform the decoupling we use the identity encircling the zero in the positive sense. We can now invoke the expression for the Cauchy determinant and the consequent identities used in previous works [14] to decouple the strings in the inverse norm factor where a set of n s real (positive) auxiliary variables has been introduced. Using properties of determinants upon line or column multiplication all summations over m j and integrations over p j can be brought inside the determinant a, using the identity This leads tô where the Kernel is given by where we have used the representation of the Pochhammer symbols (x) n in terms of Gamma functions, i.e. (x) n = Γ(x + n)/Γ(x). Note that this representation allows for an analytic continuation to complex n, which is used later. It turns out that this kernel has a very similar structure as the kernel obtained in the study of stationary KPZ in [24]. We now make this more apparent by providing an equivalent representation involving the deformed Airy function Ai Γ Γ (a, b, c, d) introduced in [24]. This function is defined as For convenience, we use an integration variable z = −iz where z denotes the one used in [24].
with Γc,d b goes from −i∞ to +i∞ passing on the left of the poles of the Gamma functions, e.g. a possible choice being a straight line from −i∞+α to +i∞+α with 0 ≤ α < Re[c/b] and 0 ≤ α < Re[d/b]. Note that the condition α ≥ 0 is necessary in order for the integral to be convergent ¶.
Our kernel can be equivalently written in terms of this deformed Airy function using the generalization of the Airy trick introduced in [24], which is valid for p ∈ R, m, ∆ 1/3 λ ≥ 0 and a, b > n/2, namely Inserting in (93), our kernel becomes where we have defined the following coefficients This expression allows to sum over the string length m ≥ 1, after performing the shift y → y + p 2 + u i + u j + s − iXp + , leading to An equivalent form for the kernel can be given, which involves now two (simpler) deformed Airy functions. These are defined as follows with similar contours Γd b as in (94) that goes from −i∞ to +i∞ passing on the left of the poles of the Gamma functions. We can now use the doubling formula given ¶ Note also that b > 0 has been assumed (so that all poles of the Gamma function are then to the right of the contour), a condition which will always be satisfied here. + forX = 0 one can check that the complex shift is legitimate, as in Eq. (68) of [67] in [24], which generalizes a similar standard identity for Airy functions. It states that, for a, b, x ∈ R and w > 0, one has where Ai Γ Γ is defined above. We have slightly generalized this identity using a change of variable p → p + iX/2 followed by a shift of contour back to the real axis, thereby generalizing the Eq.(83) of [67] for standard Airy functions. We can now use the relation (100) to integrate over p in (98) and obtain with, using (97), Note that the kernel in (102) should be multiplied by the factor e − 1 2 (u i −u j )X . However it is easy to check that this factor is immaterial when calculating the Fredholm determinant, hence we omit it and we reintroduce it only at the very end of the calculation. This kernel is thus very similar to the one for the stationary case, e.g.K in (4.22) in [24]. The parameter γ t ∼ t 1/3 there, is replaced by 2 2/3 ∆ 1/3 λ = (t 2 − t 1 ) 1/3 here. The drifts v, w there are replaced byc + ,c − . There is however an important difference. In the present case, there are additional degrees of freedom, iq and n 1 + n 2 (inc + ,c − ) that need to be summed over. To recapitulate, we have obtained an expression forẐ 1ns (s 1 , s), Eq. (92), in terms of a n s -uple integral over auxiliary variables u j , of the determinant of a kernel, K n 1 ,n 2 ,q,z , and also involving "external" integration (or summation) variables n 1 , n 2 , q, z. We obtained three equivalent expressions for this kernel. It is thus now possible to formally sum over n s to obtain a representation forĝ × e −λs 1 (n 1 +n 2 ) e 1 3 λ 3 (n 1 +n 2 ) 3 −λ(n 1 +n 2 )q 2 n 2 !(n 1 + n 2 )!(2w + n 1 ) n 2 Det[I − K n 1 ,n 2 ,q,z ]

Calculation ofĝ
It is now possible to choose (or move) the integration contour of z to encircle the negative real axis to obtain an equivalent formula. The details of the procedure are given in the Appendix F, where more general considerations about similar constrained sums are also discussed. As a result we obtain with where we recallc ± = w + n 1 +n 2 2 ± iq 2λ ±X 4∆ 1/3 λ ands = s +X 2 4 . These expressions are valid for arbitrary times t 1 , t 2 , i.e. arbitrary λ, ∆. The remaining summations over n 1 , n 2 still involves a cubic exponential which we can now eliminate using the standard Airy trick e 1 3 λ 3 (n 1 +n 2 ) 3 = dyAi(y 1 )e λ(n 1 +n 2 )y 1 For convenience we also eliminate one denominator using the identity 1/(λ(n 1 + n 2 )) = +∞ 0 dve −vλ(n 1 +n 2 ) . Upon shifting y 1 → y 1 + s 1 + q 2 + v we obtain This is our final result for arbitrary times, i.e. arbitrary λ and ∆. In the next Section we will give an equivalent expression where summations over n 1 , n 2 are replaced by contour integrals. Further progress on these formula will then be achieved, but only in the limit of large times, i.e. λ → ∞.

Large times limit ofĝ
(1) ∆,λ (s 1 , s) We now use the so-called Mellin-Barnes transformation of the sums over n 1 and n 2 . The basic identity is where each contour encircles all the solutions of sin(πz) = 0 with Re(z) > 0, i.e. all strictly positive integers. As such it is strictly true for a function f (z 1 , z 2 ) which is analytic near each of these points. The usual use of this formula is then to move these contours parallel to the imaginary axis as C j = a + iR with 0 < a < 1, and we will assume in the following that this can be done. It is useful for the following to rescale the new variables z 1 , z 2 by 1/λ, and to define a rescaled slopew as We then obtain Ai Γ Γ 2 −2/3 (y +s + 2u j ), Ai Γ Γ 2 −2/3 (y +s + 2u i ), Ai Γ Γ 2 −2/3 (y +s + 2u j ), So until now, the formula are still valid for arbitrary times, i.e. arbitrary λ and ∆. We will now take the limit of large times t 1 , t 2 → ∞, at fixed ∆, i.e. the limit of large λ in the above formula. A first simplification which occurs in that limit is that: where we recall that a choice for the contour Γ d/b is a straight line from −i∞ In addition in the denominators λ sin(πz i /λ) → πz i and the contours of integration become C 1,2 → C = + iR with > 0 infinitesimal, i.e. they pass to the right of zero. We thus obtain (using again the notation σ = 2 −2/3s and σ 1 = 2 −2/3 s 1 ) Here we have performed manipulations as in the Appendix, going there from (308) to (310) where we have changed κ → −κ and then the integration variable inside the kernel y → y − σ. For convenience we definê The new kernel is given by (see the analogous formula (309) of the Appendix) Note that we have performed a change of variable y → 2 2/3 y and a similarity transformation in the Fredholm determinant u i → 2 −1/3 u i (which absorbs also a factor 2 1/3 in front of the kernel from the integration measure).
We can now change variables to z 1 + z 2 → z and z 2 → z 2 in the complex integral, as well as κ → κ − σ, which allows to rewrite the expression aŝ We now take the limitw → 0 + . The case of a finitew is studied in the Appendix D. Here we obtain where now We note that the last arguments of the kernel are given by b =X + (4∆) 1/3 z 2 + i q 2 and c = −X + (4∆) 1/3 z 2 − i q 2 , therefore is possible to choose the contours Γ b , Γ c as X + + iR and −X + iR such that we can integrate over y using that [ and by rearranging 1 we can rewrite the kernel as in terms of the Airy kernel K Ai defined in (10). We have also used the integral representation of the Airy function (with α ≥ 0) and defined the functions B a (x) as follows The second expression is valid for Re(a) > 0, which is the domain of interest here, however the first expression is valid for all a. This function can also be determined as the solution of (a + ∂ x )B a (x) = Ai(x) which vanishes at x = +∞. Note that (128) is formula (6.22) in [24] but here it is extended to the needed complex values of the argument.
Using (125) we can rewrite the above kernel (122) as follows (128) involves a rank one projector, one can rewrite the Fredholm determinant as This expansion generates then two terms which, once inserted in equation (121), can be integrated over z and z 2 . For the first term we have simply We now use that to integrate the second piece as follows the integration over z gives a δ y + (4∆) 1/3 2 (y 1 + y 2 ) and the one over z 2 produces a θ(κ − σ), leading to Plugging this second piece in (283) we can integrate over y and move the derivative on the θ function (by integration by parts) We can then finally write an expression for the two-time generating function in the large time limit where have used the identity [57] We have used everywhere the first definition of the Airy kernel in (10). Now, recalling the definition (12) of the tail of the TW-GUE cumulative distribution we can simplify the formula, to obtain our final result for the generating function where we have introduced a new kernel K ∆ σ 1 defined as Using our starting formula (81) we then arrive to an expression for the JPDF Note that to the precision which we use here, i.e. up to higher order corrections in σ 1 of order O(e − 8 3 σ 3/2 1 ), the last term in the expression (139) can be rewritten as a difference between two Fredholm determinantŝ g (1) suitable for numerical evaluations. The kernel K ∆ σ 1 is trace-class (i.e its trace is finite) for ∆ > 0 (see below), and it is interesting to note that it can also be written as a square K ∆ Another observation, used below, is that one can write formally where the function B a (u) has been defined in (281) and (143) is taken in the sense of operators.
We will now study this result in various limiting cases, namely large and small ∆ and large σ 1 .

Small ∆ limit
In this Section we study the expansion for small ∆ of our final expression (139). For simplicity we first set X = 0, the general case is sketched at the end.
If one performs the limit ∆ → 0 directly on the expression (140) of the new kernel K ∆ σ 1 one finds, keeping also the leading correction where we recall B a (u) = 0 −∞ dye ay Ai(y + u) = e a 3 /3−ua − +∞ 0 dye ay Ai(y + u). In particular we use later ∂ a B a (u)| a=0 = −uB 0 (u) + Ai (u). The expression (144) is correct at fixed u, v, but formal as an operator since the resulting limit is not trace class (the integral of B 0 (u) diverges at large u). To make the manipulations more rigorous in the limit we now put the kernel K ∆ σ 1 in an alternative form, more suitable to take the ∆ → 0 in formula (139). First we rewrite We then use the following integration formula for the Airy functions obtained using the integral representation (126) Expanding the factors, this allows to rewrite the kernel as follows Using this formula one can easily recovers (144) at fixed u, v. Next, using this new representation we can explicity calculate the building blocks in (139) in an expansion at small ∆ and we find with where the second, third and fourth terms come from the first, second and third line, respectively in (147). The next order involves the function + ∞ σ du +∞ 0 dy 1 dy 2 y 1 Ai(y 1 + u)Ai(y 2 + u) We also find, for any integer q ≥ 1, Here we can equivalently use (144), since the convergence of the integrals involved in the trace is guaranteed by the properties of the Airy kernel. Here and below we denote interchangeably the rank one projector We can now substitute in (142) where we represent (I − P σ K Ai P σ ) −1 = I + q≥1 (P σ K Ai P σ ) q , and we obtain the expansion where a cancellation of two leading terms O(∆ 0 ) has occured. The trace in the second line can be rewritten using determinants as follows (supressing temporarily the O( We have used that P σ K Ai P σ B 0 B T 0 = |P σ K Ai P σ B 0 B 0 | is also a rank one projector, hence one can apply formula (320).
We can now calculate the JPDF for small ∆ using the formula (81) Let us first discuss ∆ = 0. Clearly in the limit ∆ → 0 the first term vanishes and the only non-zero term is given by the second term. It leads to a probability distribution factorized in a term dependent on σ 1 and a one (the non-trivial one) depending on σ, which takes the form of a second derivative, and we finally find where F 2 (σ 1 ) = K Ai (σ 1 , σ 1 ) is the tail of the PDF of the TW-GUE distribution, and F 0 (σ) is precisely the PDF of the Baik-Rains distribution, of associated CDF in a form similar to the one given in [24] in the study of the KPZ equation with Brownian initial conditions (see Appendix 7, D and E for an alternative derivation of this result, a more detailed comparison with [24] (the reader should be aware of the presence of misprints in Eqs (2.37-2.38) in [24]). The above result (156) corresponds to a product of a GUE random variable and a Baik-Rains random variable. Thus, we find here that in the limit where 1 t 2 − t 1 t 1 , the height increment h(t 2 , 0) − h(t 1 , 0) becomes statistically independent of the height at the earlier time, and furthermore, that it is distributed as in large time stationary KPZ. This result will be confirmed by a different approach in Appendix 7.
The above results about the limit ∆ = 0 extend for a finiteX, and even in presence of a slopeŵ as is shown in Appendix D and E . The BR distribution F 0 (σ) in (156) is then replaced by the extended BR distribution, ∂ σ F 0 (σ −X 2 ;X,ŵ).
We now obtain the first correction to the GUE-TW times BR distribution at small but finite ∆. Applying the differentiation operator (155) to (153) and keeping the next order correction we find Quite remarkably, a comparison with (289) shows that where we recall that F This means that the term in the first order in ∆ 1/3 is related to a shift in a finite slopeŵ of the BR distribution. This is perfectly consistent with the formula (143). Furthermore, since F 2 (σ 1 ), in the double limit of large σ 1 and small ∆ such thatŵ σ 1 ,∆ = ∆ 1/3 σ 1/2 1 is constant, the JPDF in (158)-(159) remains decoupled into the product of (the tail of) the GUE-TW distribution and the extended BR distribution with a finite slope given byŵ σ 1 ,∆ . While this property is shown here for a fixed but small value ofŵ σ 1 ,∆ 1 it is shown in Section 6.8 that it holds for arbitrary value of w σ 1 ,∆ .
From the result (159) one can calculate conditional cumulants of σ at fixed σ 1 . One finds with 6.6. Large ∆ limit Let consider again the generating function in (139). To perform the large ∆ limit it is convenient to make a change of variable y j → y j ∆ −1/3 in the definition of K ∆ σ 1 and expand as with φ pq (σ 1 ) = φ qp (σ 1 ), and we have defined what appears to be the relevant scale along x in the large ∆ limit. Note that we are using here the symmetrized version K (4),sym σ 1 of the operator K ∆ σ 1 , where eX (y 2 −y 1 ) has been replaced by cosh(X(y 2 − y 1 )). This is allowed since one can use in (139) the symmetry of the operator (I − P σ K Ai P σ ) −1 to symmetrize in u, v, hence in y 1 , y 2 and one can use there either form (symmetrized or unsymmetrized) of K ∆ σ 1 . Each of the terms in (162) is a rank one projector, hence using (320) we can expand (using that the Airy kernel is symmetric) where we have used similar notation as in (152) to denote rank one projectors. We can now expand the generating function including all terms up to O(1/∆) and obtain where we have used formula (324) and (325) to simplify the first two subleading orders in terms of the derivatives of the TW-GUE distribution.
The first consequence is that in the infinite time separation limit the joint PDF is simply the product of two PDF of GUE distributions as expected. Note that it is obtained here in the tail region for σ 1 but it is expected to hold for any σ 1 . In addition we obtain the corrections to the JPDF as where we have defined and the functions C 1 (σ 1 , σ) and C 4/3 (σ 1 , σ) have more complicated expressions, given in the Appendix H. Dividing by F 2 (σ 1 ) = K Ai (σ 1 , σ 1 ), this leads to a rather simple exact formula for the three leading orders in the large ∆ expansion of the conditional moments h p h 1 defined in Section. 3.3.1. For convenience here we use the notation σ p σ 1 ≡ h p h 1 . Using the above, and integration by parts we obtain for any integer p ≥ 1 where we recall that the values of the cumulants of the GUE-TW and BR distribution are given in (32) and (35). We have defined the functions R 11 , R 12 , R 13 and the coefficients a p , b p being given in the Appendix H. The functions R α contain all the dependence inX. Note that (171) is exact at large σ 1 , with corrections expected to be only of the order O(e − 4 3 σ 3/2 1 ). Let us now discuss some properties of the functions R 1/3 , R 2/3 and their consequences. Let us focus on the caseX = 0. Then one has and one finds that at large σ 1 → +∞ and that keeping this three term approximation gives 1 percent accuracy around σ 1 = 4.5 and 10 percent around σ 1 = 2.5. Specific values are R 1/3 (1) = 0.984861, R 1/3 (0) = 1.20144, R 1/3 (−1) = 1.552 and R( σ GUE ) = 1.87372. One also finds that at large σ 1 → +∞ and that keeping this two term approximation gives 1 percent accuracy around σ 1 = 10. and 10 percent around σ 1 = 4.5. Specific values are R 2/3 (1) = 0.40778, R 2/3 (0) = 0.59372 and R 2/3 (−1) = 0.971346. From the above formula for the moments, we obtain the expansion for the conditional cumulants Note that the term ∆ −1/3 cancels in the variance. The correction to the variance is thus O(∆ −2/3 ) with a negative coefficient (equal to −0.466009, −0.256029 and −0.154392 for σ 1 = −1, 0, 1 respectively and − 1 4σ 1 for large σ 1 ). Since the variance of the BR distribution is larger than the GUE one, this implies non-monotonicity in ∆ for the variance. For the higher cumulants we find and the O( 1 ∆ ) term is calculated explicitly in the Appendix H. Since the variance is already corrected to order O(∆ −2/3 ) we easily obtain the leading contribution to the skewness and kurtosis The case X = 0 can be analyzed similarly. One finds that as long asX < 1, i.e.X < ∆ 1/3 , there is essentially no change in the JPDF. WhenX increases beyond unity, the function R α keeps increasing, which presumably means that the present tail approximation is valid only for larger values of σ 1 asX 1.
Finally, we can study the integrated conditional moments h p h 1 >σc as defined in Section. 3.3.1. The large ∆ expansion formula for the moments and cumulants, Eqs (171),(177),(178),(181), are still valid but with the functions R α replaced as follows In the case of the mean value, it is possible to compare the above result for the conditional mean (177) and the integrated conditional mean (182), with the exact result (47) and its large ∆ expansion in Eq. (49). The comparison of the O(∆ −1/3 ) coefficient is plotted in Fig. 11. At eyesight it gives a rough estimate of the value of σ c where our tail approximation in the RBA calculation must break down (around σ c ≈ −1.5). Interestingly, the conditioned and unconditioned mean share the absence of a O(∆ −2/3 ) correction term.

Two-time conditional correlations at large ∆: persistence of correlations
Here we introduce and study the correlation functions for the two height profiles, conditioned to realizations for which H 1 is larger than a certain value. To be more precise we recall the meaning of the conditional averages as: which has been studied in experiments (see Fig. 12 in [34]) and found there to reach a non-zero limit, c = lim ∆→+∞,σc→−∞ at large ∆ as also argued from theory (see Eq. (229) in Appendix 7 here, and [45,52]). This phenomenon of breaking of ergodicity can be probed in a much finer way here, by conditioning it to the value of h 1 > σ c . Indeed, the larger σ c the more the two polymer paths are expected to overlap, and the higher the ratio should be. We now calculate this ratio as a function of σ c : our result is exact for large positive σ c and extrapolates nicely to negative values of σ c , leading to a quite reasonable estimate for c (see discussion below). Using the result of the previous Section we can rewrite the factors in (184) respectively as We can now substitute h σ 1 ≡ σ σ 1 from Eq. (177). The leading term σ GUE cancels in the numerator (186) and we obtain the limit y). The first order corrections in 1/∆ are given by and the coefficients a 1 , R 11 , b 1 , R 12 defined in Appendix H. Note that here the existence of this finite limit, and breaking of ergodicity, is demonstrated from an explicit calculation based on the RBA method. We also note that it exhibits no dependence inX, hence the ratio c σc is independent of the endpoint position in the scaling regionX = O(1). The ratio c σc is plotted as a function of σ c in Fig. 12 and it approaches 1 with corrections of order ∼ σ −1 c for large and positive σ c . As σ c decreases below ≈ −1.5 the tail approximation becomes unreliable and indeed the curve c σc in Fig. 12 exhibits a minimum at σ c ≈ −2 while the exact minimum is expected at σ c ≈ −∞. Considering the value at the minimum, and in view of Fig. 11 (where the maximum of the curve was found not far from the exact value) we can estimate the value c = c σc=−∞ = 0.58 ± 0.05 for the unconditionned ratio.

Limit of σ 1 ∼ ∆ −2/3 large, and extended Baik-Rains distribution with a slope
Here we compute the leading tail of JPDF when σ 1 is large and positive, at fixed σ. We use the asymptotics of the Airy kernel for σ 1 → +∞ at fixed y 1,2 ∆ 1/3 as follows Inserting in the definition (140) this leads to where we recall that for a > 0 Thus it is clear that the most interesting scaling regime is in the double limit where σ 1 → +∞ while ∆ → 0 so that the combination remains fixed, as well as fixedX. If instead one fixes ∆ and takes σ 1 → +∞ the expression degenerates to the one studied previously in the large ∆ expansion (see below). We thus now consider only this scaling limit. Let us setX = 0 for notational simplicity, and restore it at the end. We obtain We are now in position to compute the JPDF in this scaling limit, that we denote asPŵ we obtain which can be rewritten as (we reintroduce now the dependence onX) in terms of the kernel Hence we find that in this scaling regime the (leading order in the tail) JPDF "simplifies" into the product of the of the GUE in the variable σ 1 , and, in the variable σ the generalized Baik-Rains probability distribution (with a fixed slopeŵ σ 1 ,∆ = σ 1/2 1 ∆ 1/3 ), ∂ σ F 0 (σ −X 2 ;X,ŵ σ 1 ,∆ ) as recalled in formula (279) (and studied there). In that regime the effect of a large σ 1 is thus to bias the longer polymer path to remain closer to x = 0 at time t = t 1 . Note also that balancingX andŵ σ 1 ,∆ produces a scale, X ∼ ∆t Note that whenŵ σ 1 ,∆ 1 (e.g. large σ 1 at fixed ∆ or large ∆ at fixed σ 1 ) we obtain simply a product of two GUE-TW probability distributionŝ which can be seen by noticing that the leading behavior of the tail of a GUE-TW distribution in σ 1 is given by

Arguments using the Airy processes
In this Section we analyze the two-time problem in terms of Airy processes [2,32,59].
In the two-time problem this approach was pioneered by Corwin and Hammond [69], and recently in [45]. There the focus was on the second moment of the two-time height correlations. Here we push the approach a bit forward, and use this to discuss some properties of the full JPDF, as well as to confirm and explain some of the results obtained in this paper. We also derive new results for the large ∆ correction to the mean, and for the dependence on the endpoint positionX. Everywhere in this Section we consider the KPZ height function in the large time limit.

Airy 2 process and droplet initial condition
Let us recall a few facts about Airy processes. At large time the droplet solution of the KPZ equation takes the form Here A 2 (x) is the Airy process, which is stationary, parity invariant, with power law correlations over a scale of order unity. Its marginal PDF at one point is the GUE-TW distribution, and explicit Fredholm determinant formula exist for its multi-point correlations (in terms of the extended Airy kernel). Equation (203) is true as a process, either in x at fixed y, or y at fixed x. This means that it allows to obtain multi-point correlations of h(x, t|y, 0) in x, for a fixed value of y, or the reverse. It is not true if h(x, t|y, 0) is seen as a process in both x, y * (that process is called the Airy sheet and is much less is known about it [62,70]).

Airy stat process and stationary initial condition
Consider now the solution h w L ,w R stat (x, t) of the KPZ equation with Brownian initial condition plus drifts, h(y, t = 0) = B 0 (y) + y(w L θ(−y) − w R θ(y)), where B 0 (y) is a two-sided unit Brownian with B 0 (0) = 0 and w ± are the drifts. At large time where, denotingŵ L,R = w L,R t 1/3 , using (203) )) (205) * This implies that if information about varying y is used, e.g. as in (209), only the one-point PDF for one value of x can be obtained

Its one point CDF is given by the (extended) Baik-Rains distribution
where the function H(x; w + , w − ) is defined in definition 3 in [1] and is symmetric in w + , w − . Here we have slightly extended the Theorem 1.5 in [64], using the fact that the STS symmetry fixes some of the dependence in the variables (σ,ŵ L,R ,x). We have also introduced an alternative notation for the extended BR distribution. The RBA solution of the KPZ equation with these initial conditions, obtained in [24], shows that remembering the shift ζ = σ −x 2 in the definition (2.5) in [24].
Let us now consider now the true stationary limit w L,R = 0. At large time as a process in x. Here A stat (x) = A 0,0 stat (x) is the so-called Airy process with stationary initial data [59,63], which has been much studied. Since it arises from the Brownian initial condition with no drifts, one has [59] A stat (x) = max where the second equality is obtained by performing the changeŷ →x −ŷ in the first one. In the third one we have definedB(ŷ) ≡ B(x −ŷ) − B(x) which is also a two-sided Brownian motion (as can be checked by computing its two point correlator). Thus one can write in law that A stat (x) = √ 2B(x) + A stat (0), i.e. a sum of a two-sided unit Brownian motion and a uniform global shift, A stat (0), distributed according to the Baik-Rains distribution [1]. Since B(x) andB(x) are correlated, B(x) and A stat (0) are also correlated. Indeed A stat (x) is a non-trivial process whose multi-point correlations are given by an explicit Fredholm determinant formula (see [63]). The equalities (210)-(212) are valid a priori only for the one point PDF of A stat (x), not as a process in x (the same is true for the process A w L ,w R stat (x) above, less studied). The one point PDF of A stat (x) (see e.g. Section 2.4 and formula (63) in [65]) is the extended BR distribution (with no drifts), a specialization of (206) Note also the notation Fx(ζ) introduced in [66] (see Appendix A and formula (1.20) there), to which we prefer the notation F 0 (ζ;x). This distribution has zero mean for all x [1], and is even inx, a consequence of the statistical parity invariance of the process A stat (x).

Two-time problem
From (203) we have, for the earlier time For the later time t 2 , we write where the two terms are statistically independent. In the large time limit it becomes whereÃ 2 is a second Airy process, independent of A 2 . Here, and from now on we definê In other words, the scaled height variables introduced in (18) and (30) are expressed in terms of Airy processes This correspondence allows to make some statements about the large and small ∆ limit, following, and then extending, the analysis of [45]. In both cases one uses the property that the Airy process is locally a Brownian, i.e. one has, for a 1 and as a process inŷ The precise mathematical statements can be found in [60,61]. A crucial point, whose consequences were not emphasized before, is that B(ŷ) and A 2 (0) are uncorrelated . Note however that the regime of scales aŷ where the (uncorrelated) Brownian approximation holds presumably becomes more narrow when conditioning to a larger positive value of A 2 (0) (as studied below). Much less is known about the next order process C(ŷ;ẑ).

Large ∆ limit
We will consider the limit of large ∆. In a first stage we will keepX fixed, i.e. X ∼ (∆t 1 ) 2/3 which is the standard KPZ scaling. In that limit we can thus replace in (218), as a process inŷ a simple qualitative argument for that was pointed to us by I. Corwin. The Airy process is the limit of whereB is a unit Brownian, uncorrelated with bothÃ 2 (X) and with with A 2 (ŷ). Denoting for convenience and from now on σ 1 = h 1 and σ = h the two random variables whose JPDF is studied in this work, we arrive at the equalities in law Thus we see that to retain a non-trivial dependence in the space variable of the leading correction in the large ∆ limit we must indeed introduceX = ∆ −1/3X as was found in the RBA calculation in Section 6.6, see Eq. (164). From (222) we obtain the following expansion in powers of ∆ −1/3 Since the pointŷ * where the maximum is attained is O(1), the quadratic term ∆ −1ŷ2 in (222) can be neglected, being subdominant. One word of caution is needed here. Since we have little control on the neglected term O(∆ −1/3 ) in (222) (originating from the neglected term in (220)) the dependence onX given here and below is tentative. The only certain result is that the correction C 1 is independent ofX at fixedX (i.e. the considerations below are safest at the pointX = 0). The first and important conclusion of (223) is that for ∆ → +∞ the JPDF of σ, σ 1 becomes the product of two independent GUE-TW distributions. Let us examine the corrections. First let us compare (224) and (210). One can rewrite, as an equality in law at fixedX where we have used the stationarity of the Airy process, i.e. that A 2 (ŷ) has the same statistics as A 2 (X −ŷ), and (210). Hence we see that the first term in (224), equal to C 1 + A 2 (0), is distributed, up to a shift, as A stat (X), i.e. according to the (extended) Baik-Rains distribution. It is however correlated with A 2 (0) in a non-trivial way. Our RBA result for the conditional first moment (177) gives some information about this correlation: it predicts that at large σ 1 In fact one can also check the reverse, i.e. calculate from (223) the leading O(∆ −1/3 ) corrections to P (σ, σ 1 ). One precisely finds the decoupled form −∆ −1/3 F 2 (σ)R 1/3 (σ 1 ) the top curve of independent Brownians conditioned not to intersect. Since at any typical point there is a finite gap between the two top curves, the process is locally Brownian and uncorrelated to its height.
obtained in (168), if one assumes (226) as well as statistical independence of C 1 and A 2 (X). Thus, to that order, there is perfect agreement with our RBA result. Note that in the limit σ 1 → +∞ one can argue that the position of the maximum in (226) approaches zero asŷ * ∼ 1/σ 1 , correctly accounting for the large σ 1 asymptotics of our RBA result R 1/3 (σ 1 ) ∼ σ −1/2 1 † † Note that, since the BR distribution has zero mean, (224) implies that the leading correction O(∆ −1/3 ) to the (unconditionned) first moment is consistent with (49).
Our result for the variance (178) shows the absence of correction of O(∆ −1/3 ), which implies that the following conditional covariance vanishes This is, again, in perfect agreement with fact that the BrownianB(ŷ) andÃ 2 (X) are not correlated, as mentionned above, hence that C 1 andÃ 2 (X) are also uncorrelated. The leading correction to the variance is thus of the form and since the total coefficient of the O(∆ −2/3 ) was found in (178) to be negative, it means that there must be a non-zero (negative) covariance ofÃ 2 (X) with the next order correction C 2 . Note from (224) that the unconditionned variance of C 1 can be bound between the sum and (absolute value of) the difference of the BR and GUE variances, recalled in (32)-(35).

Persistence of correlations
The correction term ∆ −1/3 C 1 in (224) is at the origin of the property of "persistence of correlations" in the case of the droplet initial condition. The latter is usually expressed as the property that the two-time connected correlation of the KPZ heights H 1 , H 2 at two large and well separated times t 1 , t 2 , and same space point X 1 = X 2 = 0, do not decay to zero in the limit t 2 /t 1 → +∞ [34,37,45,52]. For the cross-second moment one can write where the constant c is given as a universal strictly positive constant of O(1). More generally, the non-trivial JPDF of A 2 (0) and A stat (0) = maxŷ[A 2 (ŷ)−ŷ 2 + √ 2B(ŷ)] characterizes the property of persistence † † Indeed one can argue that for events such that A 2 (0) = σ 1 1 the macroscopic profile of the Airy process minus the parabola, A 2 (ŷ) −ŷ 2 , is changed away, in a region of extension ∼ σ 1 , from minus the parabola −ŷ 2 into a triangular shape ∼ −σ 1/2 1 |ŷ|. On each side near zero the problem thus becomes similar to finding the maximum of a unit one sided Brownian with a negative drift ∼ −σ 1/2 1 , which is attained at a distance y * ∼ σ −1/2 1 near the origin. We thank I. Corwin for help in setting up this argument. of correlations. In particular the asymptotic value of the conditional correlation ratio (184) can be expressed in terms of the Airy process as i.e. the ratio of two covariances conditioned on events such that A 2 (0) > σ c . We have performed a preliminary numerical simulation to measure this quantity directly from (231) using a lattice DP model, which led to a value of c ≈ 0.6 in the range of values obtained using the RBA arguments (see previous Section) and leave a more complete and precise determination of the function c ∆=+∞,σc to the future. In the DP language, it can be viewed as a measure of the tendency of the two DP of lengths t 1 and t 2 to share a finite fraction 0 < q < 1 of the shorter DP path, since both must come through space coordinate 0 at time t = 0. This is not the case for the flat initial condition, where correlation does decay to zero. In terms of optimal paths this is particularly clear (since those tend to coalesce, see e.g. [70]) and both behaviors (for droplet and flat) were observed in a numerical simulation in [52].
What we further find here is that the constant c in (229) (and the full JPDF mentioned above) is also independent of the spatial separation of the points in the KPZ regime X ∼ (∆t 1 ) 2/3 (fixedX,X = 0). One needs a much larger spatial separation X ∼ ∆ 1/3 (∆t 1 ) 2/3 to "unbind" the two DP paths from each others, an effect which is visible in the numerical data of Fig. 6.

Small ∆ limit
To study the limit ∆ → 0, we now write the equations (218) in terms of the variableŝ X andŶ , defined in (217), as In the small ∆ limit we can use (219) and expand the second line as an equality valid in law. Hence as ∆ → 0, the one point PDF of σ converges to the (extended) Baik-Rains distribution ∂ σ F 0 (σ −X 2 ;X). In addition, since A 2 (0) and B(ỹ) are uncorrelated, the full JPDF of σ 1 , σ decouples in that limit into the product of the GUE-TW distribution and of the BR distributions. This is perfectly consistent with what is found in this paper from the RBA method in (156) and in Appendix D. There this property was shown in the tail, but the above argument shows that it holds for all couples σ 1 , σ.
Note that using the RBA we obtain that in the double limit of large σ 1 and small ∆ with fixedŵ σ 1 ,∆ = ∆ 1/3 σ 1/2 1 , the JPDF of σ 1 and σ remains a product of a GUE-TW distribution times a BR distribution in presence of a wedge of slopeŵ σ 1 ,∆ . This means that conditioning in (232) to a large positive value of A 2 (0) leads to an effective wedge, i.e. (for sayX = 0) an additional term −ŵ σ 1 ,∆ |Ŷ | in (233). It would be desirable to confirm this result using finer properties of the Airy process.
The two point connected correlation is easily estimated at small ∆ [36,45] from where the only property used is that h approaches the BR distribution for small ∆ (which must be replaced by the extended BR in the case whereX = 0). Other such identities can be obtained, e.g. from expanding H 3 21 c one obtains the time reversal antisymmetric correlation Note however that the statistical independence of H 1 and H 21 in the limit ∆ → 0 is a stronger prediction, leading to the vanishing of all lim ∆→0 H p , which could be checked in numerics and experiments.

Conclusion
In summary we have obtained some exact results for the correlation functions at two different times in the KPZ equation with droplet initial conditions, when both times are large, but their ratio t 2 /t 1 is fixed. We have used the replica Bethe ansatz method, which has been very successful for the one-time problem, but until now, failed to produce correct and testable results for the two-time problem, the main difficulty being the calculation of the so-called form factors and the summation of the resulting expressions over all eigenstates.
Our method could be described as "catching the devil by the tail". We have restricted the summation over a simpler set of states, enough to be able to perform the calculation with no further approximations, but still rich enough to capture the exact tail of the two-time joint probability distribution function (JPDF). More precisely, calling h 1 the scaled height at time t 1 and h the difference between the height at time t 2 and the one at time t 1 , we have obtained the exact behavior of the JPDF, P ∆ (h 1 , h), for large and positive h 1 and for any h. Furthermore it is expected to be quite good even for moderate values of h 1 , not necessarily very large, in particular we expect very precise results in the whole region h 1 > −1.
Our result for the tail is, therefore, rich enough to capture a notable physical aspect, namely that the JPDF factorizes into two GUE-TW distributions when the ratio t 2 /t 2 diverges, and into the product of a GUE-TW and a Baik-Rains distribution when t 2 /t 1 → 1 + . Hence the feature of stationarity at two close times, which has been conjectured and observed (to some extent) in numerics and experiments, is quite easily recovered. To our knowledge this is the first time that this feature is shown through an analytic calculation, and here we obtain it from the RBA method. Our formula therefore provides the full crossover, in the large h 1 region, from BR stationarity to GUE-TW as the time separation increases. It produces a number of results for the cumulants of h, conditioned to a fixed value, or interval of values, of h 1 that can be numerically evaluated. All these predictions being universal, we hope that they will be tested soon in numerics and experiments.
Although this still leaves the formidable challenge of obtaining the complete JPDF P ∆ (h 1 , h) as an open problem, the partial summation carried here, and the consequent results, will certainly prove themselves as a useful guide for further progress. In particular the constraint that stationarity must be recovered at close times emerges quite naturally within our approximation, and it remains to understand how it emerges more generally within the RBA. On the other hand, extensions to other classes of initial conditions seem possible, and work in that direction is in progress.

A. Comment on previous works on the two-time problem
We now discuss relations of the present work with the results reported in references [46,47]. In formula (73) of reference [46] is presented an expression for the infinite time cumulative distribution W (f 1 , f 2 , ∆) of the two-time problem with X = 0. Using the following change of notation (16) there, corresponds to the JPDF computed in the present paper, defined in (68), is while ∆ is a common notation here and there. The cumulative distribution W (f 1 , f 2 , ∆) is computed in principle by summing over all the number of replicas, i.e. all the possible string configurations for the state |γ and the state |µ (in our notations). However in Ref. [46] a substantial simplification is performed (not mentioned in the text), namely the author imposes on the starting formula (37) the following condition on the string configurations of the two states n γ s ≥ n µ s Ref [46] (238) where in the notation of reference [46] M 1 = n µ s and M = M 1 + M 2 = n γ s , with M 1 , M 2 ≥ 0. Therefore we claim that formula (73) in [46] cannot be considered a priori as the exact expression of the cumulative two-time distribution Ref [46] = P ∆ (σ 1 , σ)dσdσ 1 (239) since all the terms with different types of string contents as n γ s < n µ s must also be taken into account. In this work we focus on the cases and we show indeed that these are non-trivial terms which lead to a finite contribution to the probability distribution in the limit of large times, and for any ∆. We emphasize that the only contribution to the cumulative distribution function that is present in both works is the first diagonal contribution 1 = n γ s = n µ s which corresponds to the term g (1,1) ∆,λ (s 1 , s) in this paper, and which indeed is equal (after implementing the change of variable (236)) to the contribution M 1 = 1, M 2 = 0 in formula (73) in [46] W (f 1 , f 2 , ∆) In reference [47] another cumulative distribution is considered, however also in this case a similar approximation is done by imposing n γ s = n µ s Ref [47] (242) where in principle all the contribution with different n γ s , n µ s ≥ 0 can contribute to the final formula (20) in [47].
A non-trivial check of our claim (239) can be performed in the limit ∆ → 0. The limit ∆ → 0 of formula (73) in [46] was obtained in Ref. [48], see equation (19) there. However, this formula is not the product of a GUE-TW distribution and a Baik-Rains distributions as predicted in the present paper (see equations (26) and 28). The occurence of stationarity, hence of the Baik-Rains distribution, in the limit ∆ → 0, is an important physics requirement (also discussed in [36,45]) which is missed by the neglect of important terms in [46][47][48].
B. Calculation of one-string to multiple-string form factor

B.1. A magic identity for stationary KPZ
In the RBA study of the KPZ equation with an initial condition given by a two-sided Brownian [24], with drifts w L and w R respectively, the following integrals are needed such that G R w [λ 1 , .., λ p ] = G L w [−λ p , .., −λ 1 ] (by convention for p = 0, G R = G L = 1). In each line, the last identity is valid in the domain of parameters where the integral converge, to which we will restrict from now on. It was shown there [24] that the following "magic" identity holds P ∈Sn where the sum is over all permutations P of the n rapidities, and the factors A P are given by (60). It is the existence of this identity which allows for an explicit solution of the KPZ equation with stationary initial conditions, as obtained in [24].
Consider now the case where |γ is a single string |γ = |q, N = n 1 + n 2 . Recall that the wavefunction of a single string is given in the sector z 1 ≤ .. ≤ z N , by Hence in the sector y 1 < .. < y p < x 1 = ..x n 1 = 0 < y p+1 < ..y n 2 it becomes Putting all together we need to calculate Hence we can identify For instance, we can identify λ a ≡ −µ * a and w L ≡ w L + iq + N/2 and w R = w R − iq − N/2 + n 1 + n 2 = −iq + N/2, and rewrite F n 2 ;n 1 +n 2 where we have used that A * P = A P | µa→−µ * a . We can now use the magic identity (244) and, remembering that N = n 1 + n 2 , we obtain our final result for the form factor with n γ s = 1 and arbitrary eigenstate |µ as F n 2 ;n 1 +n 2 µ;γ = n 2 !(n 1 + n 2 )!2 2n 2 (256) If we now specify that |µ is also a string state with |µ = |p, m µ we can express the product in terms of Pocchammer symbols (x) n = x(x + 1)..(x + n − 1) = Γ(x + n)/Γ(x) leading to F n 2 ;n 1 +n 2 p,m µ ;q,n 1 +n 2 = n 2 !(n 1 + n 2 )!(w L + w R + n 1 ) n 2 which, for w L = w R = w, reduces to the formula (85) in the text. Finally note that since Re(w L,R ) ≥ Re(w L,R ) in the above identifications, there is no additional convergence problem with respect to the calculations in [24] (convergence issues were discussed and resolved there).
C. Calculation ofẐ 1,1 (s 1 , s) In this Appendix we present an independent calculation of the contribution of a single string both for the µ and γ eigenstates. It is related to the "double tail" of the two times JPDF, where both h 1 and h 2 − h 1 are large. (−1) n 1 +n 2 n 1 !n 2 ! dp 2πn 2 dq 2π(n 1 + n 2 ) × e 1 12 n 3 2 ∆t 1 −n 2 p 2 ∆t 1 − 1 12 (n 1 +n 2 ) 3 t 1 −(n 1 +n 2 )q 2 t 1 −λs 1 (n 1 +n 2 )−λ∆ 1/3 n 2 s+2iλ 2 ∆ 2/3X pn 2 We now apply Mellin-Barnes in order to compute the sums over n 1 , n 2 (259) where we assumed that there are no extra poles of f (z 1 , z 2 ) on the plane z 1,2 > 0. (Note that > 0 since we exclude the pole in z 1,2 = 0). We can rescale z 1 , z 2 → z 1 λ , z 2 λ and take the infinite λ limit lim λ→∞ n 1 ≥1,n 2 ≥1 (−1) n 1 +n 2 f (n 1 , n 2 ) = +i∞+ −i∞+ We also perform the change of variable p → p/(2λ), q → q/(2λ) and used the rescaled slopew = wλ∆ 1/3 . We obtain then in the limit λ → +∞: where in the last line we have used the Airy trick () and a shift of the variables y 1 and y 2 . The limitw → 0 + can be taken by just settingw = 0 in the above complex integral over z 1 , z 2 , since they have a small positive real part. This leads to some simplifications, and the integral is then easy to calculate, by expanding the product in the last line of (262). This leads tô In both terms we change p → p∆ −1/3 , v 2 → v 2 ∆ 1/3 , and in the first one we also change y 2 → y 2 ∆ 1/3 , while in the second one y 2 = −y 1 because of the delta function. The sum can then be rewritten aŝ We have used that We can now check that this result can also be obtained from our more general result (139) and "expanding in the number of Airy functions" containing σ, i.e. writinĝ g should identify withẐ 11 (s 1 , s). Using (137) is clear that the first terms in both formula (264 ) and (266) are identical. To check that the second terms are also the same, one first introduces the representation sin(ay) (264 ), then uses (136) for both the integrals over p and q. Integration over v 1 and v 2 yields two Airy kernels, and the change of variable y = 2 −1/3 ∆ 1/3 (y 1 + y 2 ), u = 2 −1/3 ∆ 1/3 (y 1 − y 2 ) allows rewrite the second term as

D. Small ∆ limit and extended Baik-Rains distribution with finite drifts
In this Section we first obtain a general formula for the generating function for both times large, and for arbitrary ∆ and arbitrary rescaled slopew (also called drift). Then we study its small ∆ limit and show that the JPDF converges to the product of the GUE-TW distribution and the extended Baik-Rains distribution with a finite drift. Let us start again from (120) and keep the rescaled slopew finite. We can then rewrite the factors in the complex integrals as follows We can then integrate by parts and move the derivative on the determinant and rewrite formula (120) aŝ which is a general formula for finite slopew, where the kernel is given in (119) Let us consider now the small ∆ limit. In that limit the kernel in (119) loses its dependence on z and q and becomes Then we can integrate over z and z 2 and obtain g (1) This generating function is such that when we apply the differential operator where we have first shifted y → y − 2 2/3 σ 1 , and then used that as well as, using (136), Hence we obtain lim ∆→0 P (1) where is the one-time, one-point CDF of the (scaled) height for the KPZ equation with Brownian initial conditions with (scaled) slopeŵ and (scaled) positionX [24]. It is the so-called extended BR distribution, also discussed in (208). In the limitX = 0 and w → 0, it becomes the CDF of the Baik-Rains distribution.
The derivation of the explicit distribution in this limit is recalled in the Appendix E.

E. Baik-Rains distribution as a limit
In this Appendix we expand the formula (277) for smallŵ. We recover results of Ref. [24] forŵ = 0 and obtain the next order. We want to obtain explicit formula for the first two terms in the expansion where we recall Using the identity (125) we can rewrite where we recall that Since the second term in (281) is a rank one projector, using (320), we can rewrite the Fredholm determinant in (277) as We now consider the limitŵ → 0. Note that the first part ew 3 /3−uw leads to divergences in the integrals over u in the second term in (283) in the limitŵ → 0. From (282) one obtains the following behavior asŵ → 0 The next order correction is given by  (149) and (150). By contrast, the first term in (283) behaves well in that limit, as the Airy kernel makes the integrals convergent. We thus obtain where a cancellation of the O(1) terms has occured. Let us apply now the operator 1 + ∂ σ 1 2ŵ and first focus on the order zero term, i.e. the limitŵ → 0. Then, using again (320) one obtains the (extended) Baik-Rains distribution in terms of the Airy kernel together with the two functions BX(u), defined in (282), and LX(σ) is defined in (285). It generalizes Eq. (157) to non-zeroX, in agreement with [24] and provides the one-point PDF of the A stat (X) process (see Appendix 7). One can push to the next order inŵ. One finds which is used in the text in (159) forX = 0 to compare to the effect of a small ∆.

F. Constrained summations
In this Appendix we will study constrained summations of the type for kernels K with some properties discussed below. In the absence of constraint, i.e. without the factor δ n, ns j=1 m j , the sum is simply a Fredholm determinant Det[I +K] with K(u i , u j ) = m≥1 K(u i , u j ; m). For illustration, we will consider kernels of the form so that the summation over the m variable is of geometric type and but the method is more general. Such constrained summations naturally occur in RBA calculations, where m j are string length variables and n the total number of replica, i.e. of particles in the associated LL problem. In the DP/KPZ problem n is the order of some moment of the partition sum/exponential of height field, averaged over the disorder/noise, respectively. Two examples which we detail are: (i) the single time DP/KPZ problem with point to point/droplet initial condition, studied e.g. in [14]. There (−1) n F n is directly proportional to the moment of the partition sum Z(t) = e h(0,t) = Z η (0, t|0, 0) and it has the form (290) and (291) with where λ = (t/4) 1/3 . (ii) for the present two-time problem, in Section (6) we have encountered a sum over string lengths {m j } ns j=1 constrained to have a fixed total sum j m j = n 2 . There n 2 is the number of particles/replica in the time interval [t 1 , t 2 ] and is associated to a moment of order n 2 of Z 2 . The associated F n has the form (290), although it now involves additional external integration/sum factors (e.g. over q, n 1 , n 2 ), and the kernel has the form (291), see e.g. Eq.(102) (K n 1 ,n 2 ,q,z (u i , u j ) there is the kernelK z defined below in (297), with λ → λ∆ 1/3 ).
The general method is to rewrite the constrain as a countour integral around the pole z = 0 using the identity This allows to decouple the sums and one can then perform the sum over m j leading to the countour integral of a Fredholm determinant where the KernelK z is defined as a sum over m so thatK z=1 =K. While the method is more general, for kernels of the form (291) explicit formula can be obtained. Performing the geometric sum one obtains K z (u i , u j ) = −θ(u i )θ(u j ) dy 1 1 + ze −λy A(u i , u j ; y) so that the z-dependence is explicit. In particular the factor 1 1+ze −λy has a pole for z = −e λy on the negative real axis, which, after integration over y results in a branch cut forK z along the real negative axis. Using 1 X−i0 + = P V 1 X + iπδ(X) we have, for x ∈ R 1 1 + ze −λy | z=x±i0 + = P V 1 x + e −λy ∓ iπδ(1 + xe −λy ) hence we can write explicitly the Kernel just above and below the branch cut as in terms of two new real kernels B x (u i , u j ) = −θ(u i )θ(u j )P V dy 1 1 + xe −λy A(u i , u j ; y) (300) P x (u i , u j ) = λθ(u i )θ(u j ) dy δ(1 + xe −λy )A(u i , u j ; y) where we have inserted the factor λ so thatP x has a well defined large λ limit (see below).
To perform the z contour integral in (295) we will choose the contour for the z integration to enclose the negative real axis, infinitesimally close from above and below. ♣ Hence we write, for small > 0 where f (z) = z n−1 Det[1 +K z ]. Hence the final result involves the difference of two FD Now we further specialize to the case where A(u i , u j ; y) = A(u i , y)Ā(u j , y) which is the case for both examples (i) and (ii) above, A being Airy, or deformed Airy functions, respectively. It is also convenient to choose the following parameterization for the real variable x and denote B f =B x and P f =P x . In that case we have B f (u i , u j ) = P V dyA(u i , y)Ā(u j , y) 1 e −λ(y+f ) − 1 θ(u i )θ(u j ) (306) Since P f = |U V | is a rank one projector, we can use the matrix determinant lemma (320) which implies that the dependence in P f is only linear, hence we obtain the final result It is interesting to note that: lim λ→+∞ B f (u i , u j ) = B f ,∞ (u i , u j ) = − dyA(u i , y)Ā(u j , y)θ(y + f ) (309) ♣ we can either choose it from the start, or start from a large contour encircling zero, and assume that and we see that P f = −∂ f B f ,∞ . Hence in the large time limit Let us now discuss the applications (i) and (ii). For (i) one can read off from (308) the one time PDF of the scaled free energy/height field for the point to point DP/droplet KPZ at any finite time. Indeed since Z n = n!e λns (−1) n F n for integer n, (308) provides an analytic continuation of the moments Z n for arbitrary complex n. Furthermore it is equivalent to state that: where G is an independent unit Gumbel variable (i.e. of PDF p(G) = e −G−e −G ), and f is a random variable of PDF where B f and P f are given by equations (306) and (307) with A(u i , y)Ā(u j , y) = A(u i , u j ; y) given by (293). One recovers exactly the result obtained in [14] (u there being −f here), by a completely different method. In the large time limit one easily recovers the GUE-TW distribution The application (ii) to Section (6) is now very similar, with A(u i , y) andĀ(u j , y) given by deformed Airy functions. This leads to transform formula (104 ) into (105) in the text, where n = n 2 and we use the notation κ = −f , with additional indices for the kernels.
On the other hand let us now consider the sum written as a constrained sum (−1) n 1 e λyn 1 e −λs 1 n 1 ∞ m=1 (−1) m δ n 2 ,m e −λs 2 m+yλm 1 λ(n 1 + m) Although in this simple case we know the explicit result, we can pretend to not be able to resolve the constraint δ n 2 ,m and therefore to choose a different trick to perform the sum. Let us now use the contour integral representation of the Kronecker delta. This gives which gives the correct result (315). Hence we have verified that the procedure to compute a constrained summation as the one in (316) leads to a correct result, at least in the case where one single constrained sum has to be performed.

G.1. Determinant identities
We recall two identities which are known for matrices but we assume extend readily to operators. The so-called matrix determinant lemma for any invertible operator A, and rank one projector |U V | (in quantum mechanics notation) and the so-called Sherman-Morrison formula the only singularities of Det[1 +K z ] arise from those of the factor 1 1+ze −λy , i.e. is a branch cut along the real negative axis.

I. A sum rule for the JPDF
We here show that We use the definition of the JPDF P whereĝ (1) ∆ is given in equation (139) g ∆ (s 1 , s) = 1 − F Therefore we have which shows (343).

J. Numerical details
For the numerical evaluation of equation (23) we chose a Gaussian quadrature discretization method to evaluate the inverse of the matrix I − P σ K Ai P σ as well as the final trace and we numerically perform the derivatives necessary to compute the probability density function. This procedure is very slowly converging, as compared to the calculation of the usual GUE distribution (see [68]), due to the definition of the kernel K ∆ σ 1 in equation (24) which involves integration of the Airy functions on the negative axis, where they oscillate. On the other hand, by using two different representation of the kernel (24) (equation (147)), we obtain good convergence of the first cumulants already when the number of points n in the discretization is around n ∼ 50 (See Fig.  13).
∆ (σ|σ 1 ), defined in (42), at fixed σ 1 = 1 and ∆ = 1 as a function of the number n of points used for the discretization of the kernels of equation (23).