Half-Space Stationary Kardar–Parisi–Zhang Equation

We study the solution of the Kardar–Parisi–Zhang (KPZ) equation for the stochastic growth of an interface of height h(x, t) on the positive half line, equivalently the free energy of the continuum directed polymer in a half space with a wall at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x=0$$\end{document}x=0. The boundary condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial _x h(x,t)|_{x=0}=A$$\end{document}∂xh(x,t)|x=0=A corresponds to an attractive wall for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A<0$$\end{document}A<0, and leads to the binding of the polymer to the wall below the critical value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=-1/2$$\end{document}A=-1/2. Here we choose the initial condition h(x, 0) to be a Brownian motion in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x>0$$\end{document}x>0 with drift \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-(B+1/2)$$\end{document}-(B+1/2). When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A+B \rightarrow -1$$\end{document}A+B→-1, the solution is stationary, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h(\cdot ,t)$$\end{document}h(·,t) remains at all times a Brownian motion with the same drift, up to a global height shift h(0, t). We show that the distribution of this height shift is invariant under the exchange of parameters A and B. For any \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A,B > - 1/2$$\end{document}A,B>-1/2, we provide an exact formula characterizing the distribution of h(0, t) at any time t, using two methods: the replica Bethe ansatz and a discretization called the log-gamma polymer, for which moment formulae were obtained. We analyze its large time asymptotics for various ranges of parameters A, B. In particular, when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A, B) \rightarrow (-1/2, -1/2)$$\end{document}(A,B)→(-1/2,-1/2), the critical stationary case, the fluctuations of the interface are governed by a universal distribution akin to the Baik–Rains distribution arising in stationary growth on the full-line. It can be expressed in terms of a simple Fredholm determinant, or equivalently in terms of the Painlevé II transcendent. This provides an analog for the KPZ equation, of some of the results recently obtained by Betea–Ferrari–Occelli in the context of stationary half-space last-passage-percolation. From universality, we expect that limiting distributions found in both models can be shown to coincide.


Introduction
The Kardar-Parisi-Zhang equation [1] describes the stochastic dynamics of the height field, h(x, t), of a growing interface in the continuum, as a function of time, or equivalently, of the free energy of a continuum directed polymer in a random potential as a function of its length. In one space dimension, x ∈ R, it is at the center of a vast universality class, the KPZ class, which contains numerous solvable discrete models with the same large scale behavior. Examples of such solvable models include random permutations and the associated PNG growth model, interacting particle systems (TASEP, ASEP, and variants), the stochastic sixvertex model and other stochastic vertex models, and several discrete models of directed polymers (DP). There has been many papers on the subject and we refer the reader to the reviews and lecture notes [2][3][4][5][6][7][8][9][10][11][12][13][14]. Exact results have also been obtained for the KPZ equation itself, and its equivalent system, the continuous directed polymer [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. These results have been obtained on the full line, x ∈ R, and have allowed proving the existing conjectures for the scaling exponents of the height fluctuations, δh ∼ t 1/3 ∼ x 1/2 , and to predict and classify the universal probability distributions which arise for various initial conditions. Most notably, the droplet and flat initial conditions were shown to lead, at large time, to the Tracy-Widom distributions [33,34] for the height at one point (centered and scaled by t 1/3 ) associated respectively to the Gaussian unitary and orthogonal ensembles, GUE and GOE, of random matrix theory. Some of these predictions have been successfully tested in experiments [35][36][37][38][39][40]. It is also interesting for applications [41] to study models in the KPZ class restricted to the half line x ∈ R + , for which fewer results are available at present. For some specific coupling to the boundary (at x = 0) the solvability properties can sometimes be preserved. One can usually define a parameter, that we will call A, and which will be defined below in Eq. (2.3) for the KPZ equation, characterizing this coupling. For the KPZ equation, A = 0 and A = +∞ correspond respectively to Neumann and to Dirichlet boundary conditions. The half-line KPZ equation is equivalent to a continuum directed polymer on a half-space with a wall at x = 0, which is repulsive for A > 0, and attractive for A < 0. A remarkable feature of the half-space problem is the existence of a critical value A c of the parameter A at which a phase transition occurs. In the polymer language it corresponds to a binding of the polymer to the wall if the attraction is strong enough A < A c . The existence of this transition was predicted by Kardar in [42] for the continuum directed polymer, using the replica Bethe ansatz. No prediction for the height distribution was obtained however. In mathematics, in a pioneering paper in 1999, Baik and Rains [2] proved the existence of a similar transition in the context of the longest increasing sub-sequences (LIS) of symmetrized random permutations. There it was shown that, in the unbound phase, for the droplet initial condition, and near x = 0, the (scaled) height fluctuations obey the Tracy-Widom distribution associated to the Gaussian symplectic ensemble (GSE). In the bound phase A < A c , the fluctuations are simply Gaussian. Exactly at the transition, A = A c , fluctuations also obey the Tracy-Widom distribution but the one associated to the Gaussian orthogonal ensemble (GOE). These results were extended in other models, e.g. a GSE-GUE crossover was shown as the endpoint position is varied towards the bulk in the PNG model [43]. For the TASEP in a half-space, equivalent to last passage percolation in a half-quadrant [3], similar results were obtained in [44,45] using Pfaffian-Schur processes [46], in particular concerning the crossover as the parameter controlling the boundary is varied simultaneously with the distance to the boundary. All these half-space models discussed above can be studied via the framework of Pfaffian point processes and random matrix theoretic techniques (these mod-els are free-fermionic). However these models do not converge to the KPZ equation. The models converging to the KPZ equation such as the asymmetric simple exclusion process (ASEP) or directed polymers are not directly related to Pfaffian point processes (they are nonfree-fermionic). Among those, ASEP was studied in [47][48][49] and the half-space log-gamma polymer was studied in [50][51][52]. For the KPZ equation in the half-space, a solution for the one-point height distribution near the wall valid at any time, was obtained for A = +∞, for the droplet initial condition using the replica Bethe ansatz in [53] (see also [54]). Another solution (also non-rigorous) was obtained for A = 0 in [55], with related but different methods using nested contour integral representations. In both cases the large time limit is found to be GSE Tracy-Widom , consistent with the general picture obtained by Baik and Rains as discussed above. Next, taking the limit from a half-space ASEP, a rigorous solution for the critical case A = − 1 2 was obtained [49], leading to GOE Tracy-Widom fluctuations. More recently, in [56], a solution valid for any time was found using the replica Bethe ansatz for any A > −1/2, and which leads to the GOE Tracy-Widom distribution at the critical point A = −1/2. In [57] a solution from the RBA, taking into account bound states, was obtained, in agreement with the results of [56]. A remarkable property of the KPZ equation on the full line is that the stationary measure is the Brownian motion in the sense that if the initial condition h(x, 0) is a two sided Brownian motion (with the appropriate amplitude) the PDF of the height difference between two space points is time independent. This leaves a uniform shift h(0, t) whose fluctuations, scaled by t 1/3 , was shown to follow the so-called Baik-Rains distribution, which is universal over the KPZ class. For the KPZ equation, the solution for all time with Brownian initial condition was found using the RBA in [30] and proved rigorously in [31]. Early investigations on stationary models in the KPZ universality class started with [58] in the context of the polynuclear growth model, which introduced the Baik-Rains distribution as the limiting distribution of height fluctuations. For a very similar model (TASEP), the spatial correlations were investigated in [59]. Outside the class of free fermionic models, besides [30,31] that we have already discussed, let us also mention [60] which proved the one point convergence of ASEP height function towards the Baik-Rains distribution. The aim of the present paper is to address the same problem, but for the half-line. We study the KPZ equation on the half-line with a boundary parameter A and an initial condition chosen as a unit one-sided Brownian motion with a drift, which we denote for later convenience as −(B + 1 2 ). The problem is thus determined by two parameters, A and B and we will study the phase diagram in the (A, B) plane. As we show, on the line A + B + 1 = 0 the initial condition is stationary, in the same sense as above, i.e. the PDF of the height difference between two space points remains at all time the one of the same Brownian motion. We show furthermore that the distribution of the height shift, h(0, t), is invariant under the exchange of parameters A and B. We will use two methods to obtain the exact generating function which characterizes the distribution of h(0, t), at any time t. The first one is the replica Bethe ansatz and is a generalization of the calculation presented in [56]. The second one starts from a known moment formula for the so-called log-gamma polymer [61], and takes the continuum limit to the KPZ equation. The obtained formula is valid for any point in the quadrant A, B > −1/2. We then study its large time limit, leading to the phase diagram of Fig. 1. In the (A, B) plane, the point (A, B) = (−1/2, −1/2) plays a special role, as the system is critical with respect to the boundary, and at the same time stationary, and will be called critical stationary. At this point, we show that in the large time limit the fluctuations of the interface are governed by a universal distribution, that we express in terms of a simple Fredholm determinant, equivalently in terms of the Painlevé II transcendent. In a sense, it is the analog for the half-line problem, of the Baik-Rains distribution for the full line. Inside the quadrant A, B > −1/2 the distribution is obtained as the GSE Tracy-Widom distribution and on scales A + 1 2 , B + 1 2 ∼ t −1/3 , there is a universal two-parameter crossover distribution that we obtain in the quadrant. In the so-called bound phase A < −1/2 or B < −1/2 away from the critical stationary point (including along the line A + B + 1 = 0), the fluctuations are expected to be Gaussian, except when A = B < −1/2 where we expect the fluctuations to be distributed as the maximum eigenvalue of a 2 × 2 GUE matrix. The crossover to this behavior is however beyond the scope of this paper. Finally, note that our results will be consistent in the limit B → +∞ with all previous results for the droplet initial condition, in particular with the ones in [56] for A ≥ −1/2. The level of mathematical rigor of all these results is discussed in Sect. 2.4. It is important to mention that very recently Betea et al. [62] studied stationary half-space KPZ growth for a discrete model, the last-passage-percolation with exponential weights (i.e. a zero-temperature polymer). They obtained a formula for the asymptotic height distribution, depending on several parameters controlling the distance to the boundary and the position on the line A + B + 1 near the critical stationary point. We expect, from the universality within the KPZ class, that our present result and theirs should match.  (74)] . Note that our formula allows for a very easy numerical evaluation of the CDF, given below in Fig. 3, and of the first moments, and allows us to determine tail estimates.
It would be interesting to study the distribution of h(x, t) when x is at a distance of order t 2/3 from 0, and A, B are scaled close to −1/2. This would correspond to varying the parameter η in [62]. However, while we can obtain some integral formula for the moments of Z (x, t) in the case x > 0, see (4.19) below, we do not expect that it can be rewritten as a Pfaffian formula and the asymptotic analysis would require to develop other methods. We leave this for future consideration.

Outline
First in Sect. 2 we define the models and make a summary of the main results and formulae obtained in this paper. In the two following sections we compute the moments of the polymer partition sum, i.e. the exponential moments of the KPZ field, by two methods. In Sect. 3 we present the derivation using the Bethe ansatz. In Sect. 4 we obtain the moments starting from the log-gamma polymer, and check that the two moment formulae coincide. In Sect. 5 we obtain, a Pfaffian formula for the Laplace transform generating function by summing up the moments. This leads to our first result, valid for all times and any A, B > −1/2, for the generating function as a Fredholm Pfaffian in terms of a matrix kernel. The large time limit of this formula, and of the matrix kernel, is studied in Sect. 6. In Sect. 7 we extend the method described in [54] to obtain a formula for the Laplace transform generating function in terms of a scalar kernel, valid for all times and A, B > −1/2. In Sect. 7.2 we perform the large time limit on this scalar kernel, which leads to a two parameter family of interpolating kernels near the point (A, B) = (−1/2, −1/2). From it we obtain various limits, including our formula for the critical stationary distribution.

Model
In this paper we study the KPZ equation, which reads, in dimensionless units understood here with the Ito prescription. Equation (2.2) means that Z (x, t) can be seen as a partition sum over continuum directed paths in the random potential − √ 2 ξ(x, t), with the endpoint at time t fixed at position x. Definition 2. 1 We consider the SHE on the half-line x ≥ 0 with boundary parameter A and (x, t) → Z (x, t) to be the solution to (2.2) (it can be shown that the solution is unique, see [63,64] and references therein) with the boundary condition and with the Brownian initial data, in presence of a drift −1/2 − B We have shifted by 1/2 the drift parameter to make more explicit a remarkable symmetry between parameters A and B. Indeed, we show (see Claim 4.9) that for any A, B ∈ R, we have the equality in distribution for any t > 0. When B goes to +∞, we recover a result recently proved in [65,Theorem 1.1].
When A + B + 1 = 0, the model defined in Definition 2.1 is stationary in the sense that for any fixed time t, the spatial process {Z (x, t)/Z (0, t)} x≥0 has the same distribution as {Z (x, 0)} x≥0 , that is the exponential of a standard Brownian motion with drift −1/2 − B. Equivalently, the distribution of the slope field ∂ x h(x, t) is time-stationary. Let us explain where this condition A + B + 1 = 0 comes from. In Sect. 4, we consider a discretization of the KPZ equation, the log-gamma directed polymer. We identify in Sect. 4.2 initial and boundary conditions for the log-gamma polymer which make increments of the partition function stationary in time, using a result from [61] which deals with the full-space case. The log-gamma polymer partition function converges weakly to the stochastic heat equation from Definition 2.1 at high temperature, see details in Sect. 4.3 (note that we provide only a sketch of proof of this result based on the combination of results from [66] to [65]). Hence we may pass to the limit, and taking into account the precise scalings, we obtain that the spatial process {Z (x, t)/Z (0, t)} x≥0 is stationary when A + B + 1 = 0.
There may exist other initial conditions for the half-space KPZ equation (or equivalently the stochastic heat equation) such that the slope field ∂ x h(x, t) is time-stationary. Indeed, the KPZ equation arises as a scaling limit of the height function of particle systems such as ASEP for which other stationary distributions exist (see Refs. [67][68][69][70][71]). Before presenting our main results, let us clarify the meaning of the boundary condition (2.3). As a process in x, Z (x, t) has the same regularity as a Brownian motion, hence ∂ x Z (x, t) cannot be associated to a real value. To make sense of (2.3), we say [63] where the last integral is an Ito integral and p A t (x, y) is the heat kernel on the positive half-line (i.e. it solves the equation ∂ t u = u with initial data δ y ) that satisfies the boundary condition The main consequence that we will use below is that which can be obtained by replacing Z (x i , t) using (2.5) inside the expectation and differentiating with respect to x i . In terms of directed polymers, Z (x, t) can be represented as a partition sum over directed paths (2.8) where η is a space-time white noise and D denotes the "measure on paths"; more precisely the path integral is defined as an expectation value over reflected Brownian bridges x(τ ) ∈ R + (reflected at x = 0) for a given realization of the Brownian initial condition B(y), followed by an expectation over the Brownian B. The extra δ interaction ensures the proper boundary condition at x = 0 for Z (x, t), see [55, Sect. 3.2].

Presentation of the Main Results
Our main results concern the height at x = 0, h(0, t). In general its large time behavior is expected to be where χ is an O(1) random variable, and β the growth fluctuation exponent. In the quadrant A, B ≥ −1/2, to which our exact results are restricted, one has β = 1/3 and v A,B ∞ = − 1 12 . Hence to present these results, we define everywhere the shifted variable Note however that in the so-called bound phase, which will not be studied in great detail here, we expect a different value of v A,B ∞ with β = 1/2 and different distributions for χ (see Sect. 4.6). In the limit B → +∞, i.e. for the droplet initial condition, it was found [57] that v A,+∞

Finite Time: Fredholm Pfaffian of Matrix Kernel
Our main result valid for all time t ≥ 0 and all A, B > − 1 2 is that the following generating function defined for ς > 0 can be written as a Fredholm Pfaffian Here W is a random variable, with an inverse gamma distribution of parameter A + B +1, see (3.22), independent from H (t), which enters in the construction of the generating function. The kernel K is matrix valued and represented by a 2 × 2 block matrix with elements where the dependence in parameters A, B only appears in the function (2.13) and the contour C is an upwardly oriented vertical line parallel to the imaginary axis with real part between 0 and min{A + 1 2 , B + 1 2 , 1}. The series in (2.11) can also be interpreted as a Fredholm Pfaffian, see Eq. (5.17) and Appendix B. The kernel (2.12) has a similar structure as the kernel defining the GSE Tracy-Widom distribution [72]. It is not entirely obvious that the integrals over the r i in (2.11) are well-defined, but this is the case. Indeed, one can show that (i) all the entries of K have exponential decay as r , r go too +∞, using a standard contour shift argument, see e.g. [51,Lemma 6.4] (ii) all entries of K grow at most polynomially with |r |, |r |, which can be shown using a variant of [51,Lemma 7.11].

Finite Time: Result in Terms of a Scalar Kernel
The matrix kernel in (2.12) has the structure of a Schur Pfaffian. Following [54] and Appendix B, we are able to express the generating function in terms of the Fredholm determinant of a scalar kernel (2.14) were the kernelK t,ς defined for all (x, y) ∈ R 2 + as

Large Time Limit and Phase Diagram
The phase diagram in the (A, B) plane in the large time limit is shown in Fig. 1. Qualitatively there are three regions. In the region A < −1/2 with A < B, we expect, from the results of [57] for B = +∞, that the polymer is "bound to the wall" and that the (scaled) height distribution at large time is Gaussian (see also analogous results for other models in [74], [44,Sect. 6], [51, Sect. 8.1]). By "bound to the wall", we mean that the polymer path spends most of its time at the boundary and does not significantly venture into the bulk, this phenomenon was predicted by Kardar [42] who studied the depinning of the polymer by the random environment. By symmetry the same can be expected for the region B < −1/2 with B < A, which corresponds to a polymer "bound to the Brownian". In the special case where A = B < −1/2, the nature of fluctuations is different, there is a competition between the boundary and the initial condition and we expect that the fluctuations have the same distribution as the largest eigenvalue of a 2 × 2 GUE matrix, based on heuristic arguments presented in Sect. 4.6. We have, however, no exact formula for the region A < −1/2 or B < −1/2 called the bound phase. Our exact results concern the third region, the quadrant A, B ≥ −1/2.  By symmetry the same holds for B = −1/2 and any A > −1/2. These results are natural to expect. Indeed, when B > −1/2, the initial condition has a drift so negative that the asymptotics of the height function should be the same as for the narrow wedge initial data. The limiting distribution then depends on the value of the boundary parameter A according to the Baik-Rains transition discussed in the Introduction, hence the GSE and GOE Tracy-Widom asymptotics.
To study the "critical stationary" point (A, B) = (−1/2, −1/2) we write 1 and consider the large time limit at fixed values of a, b. This corresponds to zooming around the critical stationary point as represented in Fig. 2. It is natural to expect -and we indeed show -that in the large time limit there is a two parameter family of CDFs F (a,b) (s), indexed by a, b such that Here, we first obtain a Fredholm Pfaffian formula for the CDF F (a,b) (s) with a, b > 0, by taking the large time limit of the matrix kernel formula (2.11). It reads where the large time matrix kernel K (a,b) is given in (6.5). Equivalently, a more convenient formula is obtained by taking the large time limit of the scalar kernel formula (2.14). We then obtain the CDF of the one-point KPZ height H (t) in the critical region, in terms of a Fredholm determinant, for a, b > 0 where the scalar transition kernelK (a,b) takes the form (2.22) Here the function A (a,b) (x) is defined by the integral representation where the contour is a upwardly oriented vertical line with real part between 0 and min{a, b}. Finally, introducing the operatorÂ s with kernelÂ s (x, y) = A (a,b) (x + y + s), the final and simplest expression for the cross-over CDF obtained from an algebraic manipulations of (2.21) is It is clear on this formula that if a, b → +∞ simultaneously, then A (a,b) (x) converges to the standard Airy function, andK (a,b) to the kernel associated to the GSE Tracy-Widom distribution (in the form found in [53]). This thus matches smoothly with the result (2.16) valid for any fixed A, B > −1/2. Another interesting limit, that we call F (a) (s) = lim b→+∞ F (a,b) (s), is the limit b → +∞ at fixed a, which corresponds to the droplet initial condition in the critical region for the wall parameter. A formula for that CDF was obtained, for a ≥ 0, using the RBA in [56]. It was conjectured to coincide with the GSE-GOE-Gaussian crossover introduced by Baik and Rains [74], see also [44,45] in the context of last passage percolation. This crossover was also studied in the context of spiked models of random matrices from the GSE [75]. By the A ↔ B symmetry, the case A → +∞ is similar, and we obtain the same distribution F (b) (s) = lim a→+∞ F (a,b) (s), which corresponds to the model with Brownian initial data in presence of an infinite repulsive wall (see [64] for a more mathematical interpretation). This is consistent with Tracy-Widom GOE fluctuations for any fixed A = −1/2 and B > −1/2 (a = 0, b = ∞) or fixed B = −1/2 and A > −1/2 (b = 0, a = ∞).

Stationary Critical Distribution
At the point (A, B) = (−1/2, −1/2) we have found a remarkable universal distribution, corresponding to the CDF F(s) := F (0,0) (s). Taking the limit a, b → 0 is delicate (as it is also to obtain the Baik-Rains distribution) and we found that the representation of the kernel K (a,b) as in (2.22) was crucial. The limit is performed in Sect. 7.3. We have shown that the limit is well defined, i.e. independent of the ratio r = b/a. We have obtained the result in several equivalent forms. We recall that we are characterizing the CDF F(s) such that The first form is in terms of the sum of two Fredholm determinants. Defining the following two kernels acting on functions in L 2 (R + ), The second form is expressed in terms of the CDF's of the GOE and GUE Tracy-Widom distributions F 1 and F 2 respectively, as which is very reminiscent of the formula for the Baik-Rains distribution for the full space stationary problem [recalled in (C. 12)]. The third form is expressed in terms of the Hastings-McLeod solution q(s) to the Painlevé II equation as (2.29) The first moments and cumulants are given in the Table 1 and we plot in Fig. 3, the CDF F together with its derivative, the PDF. Finally, we computed and plotted in Appendix D the asymptotics of the CDF F(s) • for large positive s  (2.31) As we mentioned above, it remains to be shown that our formula is equivalent to the Fredholm Pfaffian formula obtained in [62, Theorem 2.7] (setting there δ = 0 and u = 0) as expected from universality.

Mathematical Aspects
The results presented in this article rely on a combination of physics and mathematics methods, but we focus in this article on physics results and make clear here that most of our results are not proved according to the standards of rigor of the mathematics literature (in particular the results stated as Claims below). It remains a challenge to turn the arguments that we present here into mathematical theorems. Let us comment further on these aspects for the mathematically inclined reader.
The first difficulty from the mathematical point of view is that we cannot rigorously characterize the distribution of the KPZ equation through its moments, because they grow too fast to uniquely determine the distribution. This is why the moment generating series that we consider in Sect. 5 are actually divergent series, but the formal power series become convergent after certain manipulations and exchanges of series/integrals. For the full-space KPZ equation, it has been proved (see e.g. [55]) that these manipulations lead to the correct answer. It could be possible to overcome this issue in our case by working on a model for which the moment problem is well-defined, and take a scaling limit to the KPZ equation. Such a strategy has been implemented for instance in [22,31,78,79] for the full-space KPZ equation and in [49] for the half-space KPZ equation with A = −1/2 and droplet initial data. Another possible approach is provided by the framework of half-space Macdonald processes [51] which allows to prove Laplace transform formulae despite the divergence of moments.
The second obstacle is that in order to prove the results of Sect. 3 below, one would need to prove the completeness of the Bethe ansatz eigenfunctions. Actually, we present in Sect. 4 another approach to obtain the same moment formulae. It relies on rigorous formulae for the log-gamma polymer from [51], and we take a scaling limit to the KPZ equation. We obtain a nested contour integral formulae for the moments of Z B A (x, t) in Claim 4.7. Note that this formula allows to take x > 0. Then, for x = 0, we may move the contours together appealing to a combinatorial conjecture from Borodin, Bufetov and Corwin [55,Conjecture 5.2], and a Pfaffian structure appears. Hence we see that assuming completeness of Bethe ansatz eigenfunctions or using this conjecture leads to the same moment formula. The results from [80,81] suggest that the two problems are indeed related.
Finally, the asymptotic analysis of Fredholm Pfaffians such as (2.11) is delicate, especially in the critical stationary regime. In particular, we have first performed a large time limit for positive a, b, and then let the parameters a, b go to 0. This allowed us to benefit from the structure of the kernel (2.22), which eventually lead to a very simple formula for the distribution F(s). It would be interesting to find a generalization of the form (2.22) at finite time, and prove that the limits commute, i.e. one can take first A, B → −1/2 and then study the large time limit.

Quantum Mechanics and Bethe Ansatz
In this Section we use the replica Bethe ansatz method to calculate the integer moments of the partition sum. The equal time multi-point moments of the solution of the SHE, Z (x, t), over the KPZ noise can be expressed [82] as a matrix element of the quantum mechanical evolution operator in imaginary time of the Lieb-Liniger model [83] Here H n is the Hamiltonian of the Lieb Liniger model [83] for n quantum particles with attractive delta function interactions of strength c = −c < 0 with here an below, in our unitsc = 1. The initial state | (t = 0) is such that Since here we are considering the Brownian initial condition and interested in averages both over the Brownian and the KPZ noise we must take the initial state | (t = 0) as A simple calculation shows that 0 (x 1 , . . . , x n ) is the fully symmetric function which in the sector 0 ≤ x 1 < · · · ≤ x n takes the form We can now rewrite (3.1) at coinciding points using the decomposition of the evolution operator e −t H n in terms of the eigenstates of H n as Here the un-normalized eigenfunctions of H n are denoted μ (of norm denoted ||μ||) with eigenenergies E μ . Here we used the fact that only symmetric (i.e. bosonic) eigenstates contribute since the initial and final states are fully symmetric in the x i . Hence the μ denotes a sum over all bosonic eigenstates of the Lieb-Liniger model, also called delta Bose gas, and μ | 0 denotes the overlap, i.e. the Hermitian scalar product of the initial state (3.5) with the eigenstate μ .
We should remember now that H n is defined on the half-line x ≥ 0. The boundary condition at the wall with parameter A translates into the same boundary condition for the wavefunctions (in each of their coordinate). This half-line quantum mechanical problem can be solved by the Bethe ansatz for A = +∞, i.e. for Dirichlet boundary condition [47,84,85] (see also [86,Sect. 5.1]) and this fact was used in [87]. It can also be solved for arbitrary A, [55,86,[88][89][90][91][92][93] which led to the moment formula in [94], [56] and [57].
From the Bethe ansatz the eigenstates μ are thus Bethe states, i.e. superpositions of plane waves over all permutations P of the n rapidities λ j for j ∈ [1, n] with an additional summation over opposite pairs ±λ j due to the infinite hard wall. The bosonic (fully symmetric) eigenstates can be obtained everywhere from their expression in the sector 0 ≤ x 1 ≤ · · · ≤ x n , which reads This wavefunction automatically satisfies both 1. The matching condition arising from the δ( The allowed values for the rapidities λ i , which parametrize the true physical eigenstates are determined by the Bethe equations arising from the boundary conditions at x = L as discussed below. One will find that the normalized eigenstates ψ μ = μ /||μ|| vanish as (λ i − λ j ) or (λ i + λ j ) when two rapidities become equal or opposite: hence the rapidities obey an exclusion principle. The detailed Bethe equations, which determine the allowed values for the set of rapidities {λ j }, depend on the choice of boundary condition at x = L. However, in the L → +∞ limit, these details do not matter. For simplicity we choose a hardwall at x = L. The Bethe equations then read In the case of the infinite hardwall, these equations are also given in Ref. [85] and their solutions in the large L limit were studied in Ref. [95]. The structure of the states for infinite L is found similar to the standard case, i.e. the general eigenstates are built by partitioning the n particles into a set of n s bound-states formed by m j ≥ 1 particles with n = n s j=1 m j . Each bound state μ is indexed by a set of {k j , m j } j=1...n s where the k j 's are real numbers. These states are perfect strings [96] , i.e. a set of rapidities where a = 1, . . . , m j labels the rapidities within the string. Such eigenstates have momentum and energy The ground-state corresponds to a single n-string with k 1 = 0. The difference with the standard case is that the states are now invariant by a sign change of any of the momenta λ j → −λ j , i.e. k j → −k j . From now on, we will denote the wavefunctions of the string states as {k ,m } . It is important to note that although for A = +∞ the strings are the only solutions of the Bethe equations at large L, for finite A there are other solutions which correspond to so-called boundary bound states. These solutions have been obtained and studied in details in [57]. As we will see below we will not need them in this work.

Moment Formula
To calculate the n-th moments of Z (x, t) from formula (3.6), we need to perform a summation over the eigenstates. For A < +∞ these eigenstates contain both the string states and the boundary bound states mentioned above. Our strategy here will be similar to the one in [56], i.e. we will calculate the moments for n < 2 A + 1 which turn out to be sufficient to perform the analytic continuation in n and obtain the generating function for any A > −1/2, using a method similar to the one in [30]. The nice feature is that when n < 2 A + 1 there are no boundary bound states. To see that, consider the Table 1 in [57] which contains the classification of the boundary bound states for this problem. For A > −1/2, a m particle boundary bound state must obey m ≥ 2 A +2 ≥ 2 A+1. On the other hand from n < 2 A+1, one must have m ≤ n < 2 A + 1, which excludes the bound state. Hence we need to consider only the string states.
A formula for the inverse of the squared norm of an arbitrary string state was obtained for A < +∞ in [94] and [57], consistent with the results of [56], as with S k,1 = 1. Note that we have only kept the leading term in L as L → +∞. Inserting the norm formula (3.12) into (3.6), we obtain the starting formula for the integer moments of the partition sum with Brownian weight on the endpoint in the limit L → +∞ Here the Kronecker delta enforces the constraint n s j=1 m j = n with m j ≥ 1 and in the summation over states we used k j → m j L R dk 2π which holds also here in the large L limit: the momenta sums become continuous and one can use that the string momenta m j k j correspond to free particles as in Refs. [15,24,25,57,87].
We can simplify the factor {k ,m } (x, . . . , x) in (3.13). For the general Bethe state (3.7) (before insertion of the string solution), the x = 0 limit then reads reads Inserting the string solution we see that we can replace in (3.13) at the wall To obtain the n-th moment in (3.13) we still need to calculate the overlap {k ,m } | 0 where 0 is given in (3.5). In general it involves sums over permutations and leads to complicated expressions but in our case, a simple structure emerges akin to the one known in full-space for a few initial conditions (droplet, half-flat, Brownian). Here, as we find in the Appendix A, the result in the half-space for Brownian initial conditions is quite simple This holds under the condition that the integral converge, that is n 2 < B + 1 2 , which we will also assume from now on. Inserting the rapidities λ j of the string state one see that the denominator in the product in the overlap (3.17) read while the numerator was already calculated in (3.15). We can thus define C k, j = A 2 k, j E k, j and putting all together we obtain the starting expression for the integer moments, denoting here and below Z ( where we recall the constraint n s j=1 m j = n. Let us usec = 1 from now on. Denoting The starting formula for the moments is then where B k,m is given in (3.20) and D k i ,m i ,k j ,m j is given in (3.12) and where we recall the constraint n s j=1 m j = n.

Decorated Moments
As in Refs. [30,97], it is useful to eliminate the Gamma factor (A+B+1) (A+B−n+1) in (3.21). To this aim we introduce a random variable W ∼ Gamma −1 (A + B + 1), independent of the KPZ height, which is inverse gamma distributed with parameter A + B + 1, in this case The n-th moment of W is given by As we will see E W n Z (t) n will serve as the basis to form a Fredhom Pfaffian.

Moments in Terms of a Pfaffian
An important identity, which makes the problem solvable in the end, is that the inverse norms of the states can be expressed as a Schur Pfaffian. Introducing the reduced variables where we recall that the Pfaffian of an anti-symmetric matrix A of size N × N is defined by (3.25) and that the Schur Pfaffian is given by (see Ref. [98]) Hence the starting formula for the moments now becomes:

Moments from the Log-Gamma Polymer
In this section we compute again the moments of the solution of the SHE, Z (x, t), taking a limit of a known formula for the moments of the partition function of the log-gamma polymer on the half-quadrant square lattice. This method uses the convergence of the log-gamma  Fig. 4), where the probability of an admissible path π between (1, 1) and (n, m) is given by

Moment Formula for the Log-Gamma Polymer
and where w i, j i≥ j is a family of independent random variables such that for i > j, The notation Gamma −1 (θ ) denotes the inverse of a Gamma distributed random variable with shape parameter θ . The partition function Z(n, m) is given by The moments of the partition function were computed using half-space Macdonald processes in [51].

Remark 4.3
One may also compute mixed moments of the partition functions at several points along a down-right path.

Stationary Structure for the Log-Gamma Polymer
In this paragraph, we will need to assume α • + α 1 = 0. In order to do so, we need to consider a modified partition function where we have removed the weight w 1,1 , i.e. we define Following [61], we define horizontal and vertical increments of the partition function as The partition function satisfies the recurrence (4.5) We will need the following lemma from [61] where the stationary structure for the full-space log-gamma polymer was introduced.

Lemma 4.4 ( [61, Lemma 3.2]) Let U , V , w be independent random variables. Let
If for some α > 0 and θ ∈ (−α, α), Coming back to our model Z stat (n, m), when α • + α 1 = 0 the model is stationary in the following sense.   Fig. 5). We associate to this down-right path increments I j 1≤ j≤k−1 where Then the increments I j 1≤ j≤k−1 are all independent and distributed as I j ∼ Gamma −1 (α 1 + α) when I j is a horizontal U increment, and I j ∼ Gamma −1 (α • + α) when I j is a vertical V increment. In particular, for any m, the increments U n,m n≥m+1 are independent and distributed as U n,m ∼ Gamma −1 (α 1 + α).
Proof The distribution of increments along the first row is completely constrained by the definition of the model. Indeed, we have that Z stat (n, 1) = n i=2 w i,1 , so that the increments along the first row are given by U n,1 = w n,1 and the definition of the model implies that weights w n,1 ∼ Gamma −1 (α 1 + α) are independent. Hence, for m = 1, the increments U n,m n≥2 are independent and distributed as Gamma −1 (α 1 + α) as claimed.
In other terms, we have seen that the statement of the Proposition is true for the infinite path going through the points (n, 1) for all n ≥ 1. We will show that the property is preserved under two types of local transformation of paths, depicted in It is clear that any infinite down-right path can be obtained by iteration of these local transformations, starting from the path going through the points (n, 1) for all n ≥ 1. Furthermore, if the statement of the Proposition is true for any infinite down-right path, it is true as well for any subpath of the form {(n i , m i )} 1≤i≤k as in the statement of the Proposition. Hence we only need to show that the distribution of increments is preserved under the two local moves.
The distribution of increments on the boundary is constrained by the definition of the model. We have that Z stat (n, n) = w n,n Z stat (n, n − 1), so that V n,n = w n,n and we recall that w n,n ∼ Gamma −1 (α • + α n ) is independent from all other weights. Hence, for any n = m, V n,m is distributed as Gamma −1 (α • + α n ) and is independent from the increments {U n ,m , V n ,m } for m < m (since those increments are independent from w n,n ). Hence, after a boundary update, the distribution of increments is preserved.
In order to show that the property is preserved under bulk update, we use Lemma 4.4. After a bulk update, the increments are updated according to (4.5), where w n,m is independent from the increments on the earlier path and distributed as Gamma −1 (2α) (recall that we have assumed that α 2 = α 3 = · · · = α). If α • +α 1 = 0, we may set θ = α 1 = −α • and Lemma 4.4 implies that increments along the new path will be distributed as U n,m ∼ Gamma −1 (α + α 1 ), V n,m ∼ Gamma −1 (α +α • ). This shows that the distribution of increments is preserved under bulk update. Because before the bulk update, the variables (U , V , w) are independent from the rest of the increments I j by induction, and the new random variables (U , V ) are just measurable functions of (U , V , w), the new variables are also independent of the other increments I j . This concludes the proof.
One consequence of the stationary structure is that we may compute the expectation of log Z stat (n, m). We assume that parameters α i are chosen as in Proposition 4.5. Observe that log Z stat (n, m) is equal to the sum of the logarithms of increments of the partition function along any path from (1, 1) to (n, m). These increments are not independent, so that their sum is a highly non trivial random variable, but we know the expectation of each increment. Since the vertical increments are distributed as Gamma −1 (α • + α) and the horizontal increments are distributed as Gamma −1 (α 1 + α), we have that (for α • + α 1 = 0) where we have used that E[log(Gamma −1 (θ ))] = −ψ(θ) and ψ is the digamma function.

Convergence to the Half-Space KPZ Equation
At high temperature (when the parameters of inverse gamma random variables go to infinity and space-time coordinates are rescaled appropriately), the partition function Z(n, m) converges to the multiplicative noise stochastic heat equation on R ≥0 with Robin type boundary condition [66].
Although the convergence of discrete directed polymers to half-space KPZ equation was proved rigorously in [66] (based on the full-space analogous result in [99]), we will rederive (heuristically) this convergence in order to adapt it to our units and initial condition (which is not covered in [66]). Let us change coordinates and use more natural time and space coordinates τ = n + m − 2 and κ = n − m. The partition function Z d (κ, τ) := Z(n, m) satisfies the discrete version of the stochastic heat equation where w κ,τ ∼ Gamma −1 (γ κ,τ ) with parameter γ κ,τ = α n + β m (independent for each κ, τ). The boundary condition at κ = 0 is given by Let us renormalize Z d and define Z r (κ, τ) = C −τ Z d (κ, τ). The correct factor to use is such that C τ behaves asymptotically as the point to line partition function where weights would be replaced by their average. Hence we set C = 2E[w κ,τ ]. Note that in the following, we will choose parameters so that E[w κ,τ ] does not depend on τ, κ (except along the lines τ = κ or κ = 0). We may rewrite (4.8) as where b κ,τ = 2w κ,τ C − 1, ∇ τ is the discrete time derivative and κ is the discrete Laplacian. Let us fix α • ∈ R, α 1 ∈ R and set α i = 1/2 + √ n/2 for all i ≥ 2 and use the scalings In this case, we may choose C = 2/ √ n and the family of random variables w κ,τ rescales to a white noise in the sense that n b κ,τ ⇒ √ 2ξ(x, t). At τ = κ, we have that for large κ, where the inverse Gamma random variable (coming from w 1,1 ) and the Brownian motion B(x) are independent. It is then natural to define the continuous limit Under the scalings that we consider, the boundary condition (4.9) becomes (4.14) Let us take the average on both sides of (4.14). We use that E[ w C ] = 1 1+ 2α•−1 √ n , and consider that the weight w is independent from Z ∞ ( 2 √ n , t), as this is true in (4.9). We obtain which, by Taylor approximation, leads to Note that one may also obtain the more general boundary condition (2.7) for mixed moments by multiplying both sides of (4.14) by Z ∞ (x 2 , t) . . . Z ∞ (x n , t) before taking the average. Finally, multiplying (4.10) by n we obtain, when taking formally the n → ∞ limit, that Z ∞ (x, t) should satisfy the SHE (2.2). Thus, we have arrived at the following.   (A + B + 1)).
Note that the derivation that we presented above is only heuristic. We will not provide a complete proof of this result, though we indicate where the needed arguments can be found. The convergence of the polymer partition function for the half-space log-gamma polymer to the multiplicative noise stochastic heat equation is proved in [66] using a chaos series representation of the polymer partition function, see in particular Section 5 therein. However, the setting of [66] restricts to delta initial data B = +∞ (and Robin type boundary with arbitrary parameter A). The convergence for Brownian initial data with arbitrary parameter B was proven in [65, Theorem 2.2], though [65] works only in the case of Dirichlet boundary condition, that is in the case A = +∞. Hence, one needs to combine the arguments from [66] and [65] to deduce this result.
Using the stationary structure from Sect. 4.2 together with Claim 4.6, we obtain the following. Let Z (x, t) be as in Claim 4.6 and assume that A + B + 1 = 0. Then, for any time t > 0, Z (x, t)/Z (0, t) is the exponential of a Brownian motion with drift −B − 1/2.
We may also compute the expectation of h(x, t) = log Z (x, t) in the stationary case when A + B + 1 = 0. Using (4.7) and plugging there the scalings of Claim 4.6, we obtain that In particular, when A = B = −1/2, we have that E [h(0, t)] = −t/12.

Moments of the Half-Space KPZ Equation with Brownian Initial Data
Using the moment formula from Proposition 4.2 and the convergence result from Claim 4.6, we obtain the following moment formula for the half-space stochastic heat equation Z (x, t).
Note that the formula is valid for any x ≥ 0.

Claim 4.7 Let Z (x, t) be the solution to the half-space stochastic heat equation (Definition 2.1) with Brownian initial data with drift −1/2 − B and boundary parameter A.
Assume that B > 0, and A + B > k − 1. Then, we have where the contours are chosen so that B > r 1 >

Remark 4.8 One may also compute mixed moments of
Proof We start from the moment formula given in Proposition 4.2. Under the scalings considered in Claim 4.6, the second line of (4.1) becomes Using dominated convergence, one readily obtains that Finally, we have used the change of variablesz i = −z k−i+1 to obtain the statement of the Proposition.

Symmetry Between Drift and Boundary Parameters
We now exploit one of the symmetries of the log-gamma polymer, arising from more general symmetries of so-called half-space Macdonald processes [51], which, in the KPZ limit, extends a result of Parekh [65].

Claim 4.9 Let us denote here by Z B A (x, t) (with explicit dependence in the parameters A, B) the solution to the multiplicative noise SHE with boundary parameter A and initial drift
is defined as in Sect. 4.2 and w 1,1 has the same distribution in both cases. This implies that we have also the equality in distribution (4.23) The distribution of the random variable in (4.23) depends on parametersα • ,α 1 , α 2 , . . . , α n . The equality in distribution can be analytically extended to all parameters such that α i +α j > 0 for any 2 ≤ i < j ≤ n,α 1 + α i > 0 for any 2 ≤ i ≤ n, andα • + α i > 0 for any 2 ≤ i ≤ n.
In particular, we do not require anymore thatα • +α 1 > 0. Passing to the limit in (4.23) using Claim 4.6, we obtain the desired result.

Another Conjectural Identity in Law
The symmetry between parameters A and B stated in Claim 4.9 relies on a similar property for the log-gamma polymer (symmetry between α • and α 1 ) based on the theory of half-space Macdonald processes ans stated as [51,Proposition 8.1]. This result was stated in [51] for the log-gamma polymer model where α 2 = · · · = α n . However, the same property actually holds for general half-space Macdonald process and for any choice of parameters, see [51, Proposition 2.6], as long as we restrict to the partition function on the boundary. Furthermore, the law of Z(n, n) is symmetric with respect to permutation of the parameters α i . Thus, we claim that one can also exchange the roles of the parameter α • and the parameter α 2 .
When scaling parameters to the (multiplicative noise) stochastic heat equation in Claim 4.6, we set α • = A − 1 2 , α 1 = B − 1 2 and α i = √ n 2 . If we exchange parameters α • and α 2 , we expect that we will obtain the stochastic heat equation with boundary parameter equal to +∞ (i.e. with Dirichlet boundary condition Z (0, t) = 0) and initial condition given by where B 1 and B 2 are independent Brownian motions with respective drifts −(B + 1/2) and −(A + 1/2). Let us call Z A,B Dir the solution to the heat equation with Dirichlet boundary condition and the initial condition given above in (4.24). We refer to [64] regarding the exact meaning of the Dirichlet boundary condition in this context. Following similar arguments as in the proof of [64, Theorem 1.1], we conjecture that for any t > 0, we have the identity in distribution The identity in law (4.25) allows to predict the law of large numbers for h(0, t) = log Z B A (0, t) in the bound phase. Recall that we expect that h(0, t) follows the asymptotics where χ is an O(1) random variable, and β the growth fluctuation exponent. Using (4.25), should also follow the same asymptotics. When A < −1/2 or B < −1/2, we see that Z A,B Dir (x, 0) grows as e max{|A+ 1 2 |,|B+ 1 2 |}x . Hence, the polymer partition function will be dominated by paths from (x, 0) to (0, t) where x is of order t. More precisely, since the point to point free energy from (x, 0) to (0, t) behaves asymptotically as − t 12 − x 2 4t , the partition function will be dominated by paths leaving from We can even predict the fluctuation exponent β and the nature of fluctuations. When A < −1/2 or B < −1/2 with A = B (say A < B for simplicity), the initial condition for Z A,B Dir (x, t) in (4.24) will essentially be the Brownian motion, i.e. B 2 (x) with drift −(A+1/2). The fluctuations of the initial condition at the optimal point x * t are thus ≈ B 2 (x * t ) i.e. Gaussian on the scale t 1/2 , and they will dominate the fluctuations in the partition function, hence we find that β = 1/2 and χ is Gaussian. The situation is completely similar when B < A.
When A = B < −1/2, the situation is a bit more delicate since we cannot approximate the initial condition (4.24) by a Brownian motion. Instead, we notice that (4.24) can be interpreted as the partition function in the O'Connell-Yor directed polymer model [101]. For large values of x, it behaves as It was proved [102,103] that the latter quantity behaves asymptotically as the largest eigenvalue of a 2 × 2 GUE matrix in the scale x 1/2 . We must now evaluate this quantity at the optimal point x = x * t = |2B + 1|t, hence we find that β = 1/2 and that χ has the same distribution as the top eigenvalue of a 2 × 2 GUE matrix, discussed for instance in [101,Sect. 6].

Residue Expansion
In this section, we restrict to x = 0 and denote Z (0, t) = Z (t) as in the previous sections. The moment formula (4.19) is not convenient for asymptotic analysis because the contours are different for each variable, so that the complexity of the formula significantly increases with k. To overcome this issue, one has to deform the contours to all lie on a fixed vertical line and take into account the residues encountered during this contour deformation. This procedure was implemented in [55], but the computation of residues is very involved and the result relies on a conjectural combinatorial simplification. Applying this result [55,Conjecture 5.2] to the moment formula (4.19), we conjecture that for A > k − 1 and B > k − 1, we have where we use the Pochhammer notation for rising factorials (4.30) and where (4.31) It turns out that the symmetrization can be performed using [55,Equation (54)], relying on the theory of BC-symmetric Hall-Littlewood polynomials [104], and one finds (4.32)

Remark 4.10
It is now apparent that the moment formulae are invariant with respect to the transformation (A, B) → (B, A).
Comparing with the formula (3.27) obtained from the replica Bethe ansatz, we see that (4.35) and (3.27) agree after the substitutions Indeed, under this change of variables, we have where X j = m p + 2ik p for j = 2 p − 1, X j = m p − 2ik p for j = 2 p. Thus we have that One may also check that Using the reflection formula for the Gamma function, we may write the hyperbolic sine in the definition of B k,m in (3.20) as a product of Gamma functions so that we have (dropping indices of variables k j , m j , λ j , w j ) (4.39) Using the duplication formula for the Gamma function and the reflection formula a few times, we arrive at Thus, we have shown that (4.35) and (3.27) agree as claimed.

Generating Function in Terms of a Fredholm Pfaffian
We will now write the moment generating function of Z (t). We define, for ς > 0, Ignoring the fact that the summation over n cannot be exchanged with the expectation due to the divergence of moments, we will consider the following formal power series that we will again denote by g(ς).
Equation (5.3) is a positive temperature analogue of [105,Eq. (4.5)] which we will use in the next section. This type of arguments can be traced back to the work of Baik and Rains [58].

Remark 5.2
Alternatively, one may use the density of the inverse Gamma distribution and write We may rewrite this equation using the following representation of the Bessel function for x > 0, where K ν denotes the modified Bessel K function. Exchanging the expectation with the integral over u in (5.4), and using this integral representation, we arrive at where K ν (z) denotes the modified Bessel K function. Note that similar integral transforms involving the modified Bessel function K appear in [31,Theorem 2.9]. It is plausible that this expression can be inverted to compute the distribution of Z (t). For A + B + 1 = 0, an inversion formula is provided in [31, Appendix E]. We will see in Sect. 6 that we will actually not need to perform this inversion.
The constraint n s i=1 m i = n in (3.27) can then be relaxed by reorganizing the series according to the number of strings: where Z (n s , ς) is the partition sum at fixed number of strings n s , calculated below. We now show that one can write the generating function as a Fredholm Pfaffian. It will be possible thanks to the Schur Pfaffian identity, (3.24), given above. The partition sum at fixed number of strings, expressed in terms of the reduced variables X 2 p−1 = m p +2ik p and X 2 p = m p −2ik p for p ∈ [1, n s ], reads where B k,m was given in (3.20). The summation over the variables m p can be done using the Mellin-Barnes summation trick similarly to Refs. [30,97]. The barrier A > (n − 1)/2 is overcome exactly as in Ref. [97] (see Lemma. 6 and the discussion therein) from an analytic continuation of Gamma functions included in the B k,m factor, the introduction of a particular contour C 0 and a final requirement for the drift A + 1/2 > 0. Indeed, define the contour C 0 = a + iR with a ∈ (0, min{2B + 1, 2 A + 1, 1}), then for any holomorphic function f having sufficient decay at infinity and in particular denoting the summand of Eq. (5.8) by the function f (m p ), we have For each m p we therefore introduce two variables r p and w p and we redefine the reduced variables X 2 p and X 2 p−1 under the minimal replacement m p → w p imposed by the Mellin-Barnes formula, which we will apply despite the presence of poles on the right of the contour C 0 . This is an a priori illegal step, but it will exactly turn the diverging moment generating series into a well-defined and converging series equal to the Laplace transform. We refer to [55,Sect. 6] where this procedure and its degree of rigor is discussed in great details. This leads to the following rewriting of the coefficient Z (n s , ς) as (see section 5 in [54] for similar manipulations)

Remark 5.3
Note that the contour C 0 passes to the left of the poles of the Gamma function at X = 2 A + 1, 2B + 1.
We observe that the integrals are almost separable in X 2 p and X 2 p−1 except for the sine function which couples them. By anti-symmetrization and similarly to [54,Sect. 5], we can proceed to the replacement 2 The last manipulations consist in rescaling all variables X by a factor 2 and replacing the contours of integration by C = a 2 + iR. Hence we have There are a few last steps before we introduce the Fredholm Pfaffian. First define the functions Using a known property of Pfaffians (see De Bruijn [106]), we can rewrite the partition sum at fixed number of strings itself as a Pfaffian, i.e. we use that 14) This leads to the definition of a 2n s × 2n s matrix M such that Since a variable r p will be shared between four elements of this matrix, it is more convenient to view M as composed of 2 × 2 blocks which we denote K , whose elements are presented in Eqs. (2.12). Finally, the string-replicated partition function is given by an infinite series of Pfaffians This series is a Fredholm Pfaffian, where K is given in (2.12), the function σ ς is given by σ ς (r ) = ς ς +e −r and the 2 × 2 kernel J is given by J (r , r ) = 0 1 −1 0 1 r =r . For the precise definition and properties of Fredholm Pfaffians see Section 8 in [107], as well as e.g. Section 2.2. in [44], Appendix B in [108] and Appendix G in [24,25].

Large Time Limit of the Fredholm Pfaffian and the Distribution of the KPZ Height: Crossover Kernel
We will now study the large time limit of our kernel. To understand the scaling required at large time, let us recall the expression of the partition sum at fixed number of strings At large time, we want to eliminate the time factor in the exponential, hence we perform the change of variables In the large time limit, the Gamma, cosine and sine functions simplify using that for small positive argument Under these simplifications, the partition sum at fixed number of strings reads in the limit t → +∞ (dropping all tildes) The contours C have now to be understood as C =ã + iR, whereã ∈ (0, min{a, b}). We emphasize that the contours all lie at the left of the poles at X = a and X = b. We may now use (5.14) and (5.16) to write the Laplace transform g(ζ ) under the scalings in (6.2) as the Pfaffian Pf[J − σ ζ K (a,b) ] L 2 (R) . The matrix valued kernel K (a,b) reads in this limit

Remark 6.1
The kernel K (a,b) has a particular structure, indeed its elements are related through derivative identities: K K (a,b) can be obtained equivalently from the kernel (2.12) by rescaling, as in [56].

Remark 6.2 The kernel
where is the Theta Heaviside function. Note that since W is an inverse Gamma random variable with parameter A + B + 1 = t −1/3 (a + b), the random variable log(W )/t 1/3 weakly converges to an exponential random variable with parameter a + b.
At this point, we obtain that for any a, b > 0, where E is an exponential random variable with parameter a + b independent from H (t), and the matrix kernel K (a,b) is given in (6.5). Using the density of the exponential distribution, and denoting F (a,b) (s) = lim t→∞ P H (t) t 1/3 ≤ s , we may rewrite (6.8) as Following [105, Eq. (4. 3)] (see also Remark 5.1), we differentiate in s in (6.9) and use integration by parts in the left hand side. We obtain Finally, we can write K (a,b) ) L 2 (s,+∞) := F (a,b) (s). (6.10) In the next sections, we show that the distribution F (a,b) interpolates between various known distribution as b or a goes to infinity. The most interesting case, corresponding to stationary growth, is when b, a go to zero and will be studied in details in Sect. 7.

Remark 6.3
Performing the large time limit at any fixed A > −1/2 and B > −1/2 corresponds -modulo an exchange of limits -to the scaling considered above with a, b = +∞. One can show that in the limit a, b → +∞, the kernel (6.5) converges to the GSE matrix kernel, by the same manipulations as in [56,Sect. 4.1]. Indeed, as the contours of kernel K (a,b) are parallel to the imaginary axis and cross the real axis between 0 and min{a, b}, we can take the limit a, b → ∞ in the integrand without affecting the contours. All rational functions involving the parameter a, b in the large time limit of the kernel K (a,b) in (6.5) converge to the value −1. Hence in this limit we obtain the kernel K ∞ given in [56,Eq. (113)], which is precisely the kernel associated to the Gaussian Symplectic Ensemble (GSE) of random matrices as also given in Lemma 2.7. of [44]. Hence this shows that the distribution of the height at x = 0 converges at large time for boundary conditions such that a, b → ∞ (e.g. for any fixed A, B > −1/2) to the GSE Tracy-Widom distribution, as we will also show below.

Remark 6.4
In the limit b → +∞, the kernel K (a,b) in (6.5) converges to the kernel K with = a obtained in [56,Eq. (64)] in the study of the droplet initial condition. The Fredholm Pfaffian F (a) (s) := Pf(J − K (a,+∞) ) L 2 (s,+∞) interpolates between the CDF of the GSE Tracy-Widom distribution, F 4 (s), (at a → +∞) and CDF of the GOE Tracy-Widom distribution function, F 1 (s), (at a = 0). Note that the GOE kernel obtained in [56,Eq. (83)] was found to provide a new representation of the GOE Tracy-Widom distribution. If instead we take a → +∞ for fixed b, we obtain again the same kernel K b because of the symmetry a ↔ b. Physically, it describes the CDF of the distribution of the rescaled height of the KPZ equation with Brownian initial data and Dirichlet boundary condition.

Solution for the KPZ Generating Function at All Times for Generic A, B in Terms of a Scalar Kernel
The general kernel we have obtained in ( and the kernelK t,ς such that for all (x, y) ∈ R 2 sin(π(z − w)) sin(π(z + w)) ς w+z e −xz−yw+t w 3 +z 3 3 (7.2) Then, the Laplace transform of the one-point distribution of the exponential of the KPZ height admits the following representation:

Large Time Limit of the Scalar Kernel
To deduce the large time asymptotics of H (t) from the determinantal formula (7.3), one performs the same rescaling as in Sec. 6, namely one chooses ς = e −t 1/3 s and one rescales (w, z) → t −1/3 (w, z), r → t 1/3 r . The kernelK t,ς becomes where the contour C is an upwardly oriented vertical line with real part between 0 and min{a, b} as previously. Then, using the defintion of F (a,b) from (6.10), the determinantal formula (7.3) implies the following: for any a, b > 0, where the contour is a vertical line with real part between 0 and min{a, b}. Note that the function A (a,b) has exponential decay at +∞, that is for any c ∈ (0, min{a, b}), there exist C ∈ R such that A (a,b) (x) ≤ Ce −cx . Let us introduce an operatorÂ s acting on L 2 (0, +∞) with kernelÂ and an operatorK where all operators act on L 2 (0, +∞).
Proof Equation 7.5 implies that, as operators acting on L 2 (0, +∞), we havē At this point, we recognize that the operatorK (a,b) s has the same structure as in [87] and we may apply the same steps as Equations (32) A (a,b) for fixed a, b > 0).

CDF in Terms of a Fredholm Determinant
Moving the contour to the right in the definition of A (a,b) in (7.6), we obtain 13) where in the integral over z, the contour passes to the right of b, a. We introduce the operator A s acting on L 2 (0, +∞) with kernel A s (x, y) =Ã (a,b) (s + x + y).

We thus writê
14) with f b (x) = √ be b 3 /6−bs/2−bx . Using the matrix determinant lemma, we have We now consider the limit b, a → 0 with a fixed arbitrary ratio r = a/b. We use the exact expressions for the scalar products Plugging these asymptotics in (7.15) we obtain  (7.27) which is given in the Introduction in (2.26).

CDF in Terms of the Solution of the Painlevé II Equation
In this Section, we show that the distribution F(s) can also be written in terms of the Tracy-Widom distributions F 1 (s) and F 2 (s) only. We also provide formulae in terms of the . (7.30) where the last identity is from [112]. The next result allows to rewrite the resolvant R +,0 s appearing in (7.24) and (7.25) in terms of q.

1|
Ai Proof Before proving this proposition, we need the following known lemma.

Lemma 7.5 ([109, Lemma 3]) We have the following relation
where D is the derivative operator.
Since D |1 = 0 we have  Since R +,0 s → 0 in the limit s → +∞, and using (C.2) and the asymptotics (C.5) and (C.6), we obtain that κ = 0, which concludes the proof. Using Proposition 7.4 we arrive at the following equivalent expressions for F(s) (defined in (7.24)). The first one reads where the first line was given in (2.28). The second formula is expressed in terms of the Hastings-McLeod solution of the Painlevé II equation (see Appendix C) and reads 2 ) . (7.38) which was given in (2.29). A third formula, useful for the asymptotics, is obtained applying the derivative in front of the bracket and using Eqs. (C.2) and (C.3)

Properties of F(s): First Moments
We check in Appendix C using the formulae (7.39) that the function F(s) has the behaviour at s → ±∞ that is required for a CDF, i.e. its limit at s = −∞ is 0 and its limit at s = +∞ This formula can be used to obtain the moments and the cumulants through a numerical evaluation of F(s). One notices that the mean vanishes, M 1 = 0, which indeed must be the case since E[H (t)] = 0 in the critical stationary case (see Sect. 4.3 where we have computed the expectation of h(x, t) using the stationary structure of the log-gamma polymer). We used two numerical methods to evaluate F(s). The first one uses Eq. (7.27) where the Fredholm determinants are calculated using the method described in Ref. [76,113]. The second method uses the formula (7.37) and uses the Mathematica routines for F 1,2 (s), and is in agreement with the first one. The CDF F(s) and its derivative F (s) are plotted in Fig. 3. The mean, variance, skewness and excess kurtosis are given in Table 1.

Limit a, b → +∞: Convergence to the GSE
As we already discussed in Remark 6.3 the limit a, b → +∞ can be performed on the Fredholm pfaffian formula and leads to GSE Tracy-Widom fluctuations. This limit can also be performed on the formula (7.11) involving the scalar kernelÂ s defined in (7.7) in terms of the function A (a,b) (x) defined in (7.6

Limit (a, b) → (0, +∞): Convergence to the GOE
Another interesting limit, that we call F (a) (s) = lim b→+∞ F (a,b) (s), is the limit b → +∞ at fixed a and by the A ↔ B symmetry, the case A → +∞ at fixed b is similar. In particular we consider now the limit a → 0. The manipulations follow closely the ones of Section. 7.3. We start by moving the contour to the right in the definition of A (a,∞) in (7.6) already taking into account the b → +∞ limit. We obtain A (a,∞) (x) =Ȃ (a,∞) (x) + 2h a (y), with h a (x) = ae −xa+a 3 /3 and (7.41) where in the integral over z, the contour passes to the right of a. We introduce the operatoȓ A s acting on L 2 (0, +∞) with kernelȂ s (x, y) =Ȃ (a,∞) (s + x + y). We thus writê Looking at formula (7.11) in the limit b → +∞ we see that we need only the following estimates from ( These integrals can be explicitly evaluated G n,w (λ 1 , . . . , λ n ) = n j=1 −1 − j(B + 1/2) + iλ n + · · · + iλ n+1− j + j 2 /2 (A.5) Now in (A.3) for each permutation P we can relabel all the ε p → ε P( p) and denoting by ε={±1} n the operation of summation over all the variables ε i (an operation independent of their labeling) we can rewrite (A.3) as where we have used the fact that the products are independent of the permutation P. Now we use the following symmetrization identity, given in [114] Applying it to the set {ε k λ k } we obtain It remains to perform the symmetrization over ε.
Remark A. 2 In the limit A → +∞, this formula yields back the overlap for Dirichlet boundary condition, see [54].
Proof This identity was essentially known in the context of Hall-Littlewood polynomials of type BC. In particular, a generalization of it was proved in [104,Theorem 2.6]. We refer to [55,Eq. (54)] for more details about how to degenerate Venkateswaran's symmetrization identity to the one we need (trigonometric to rational limit). For any parameters A, B, C ∈ C and complex variables (z i ) 1≤i≤n , we have ε∈{±1} n P∈S n k<l ε P(k) z P(k) + ε P(l) z P(l) − C ε P(k) z P(k) + ε P(l) z P(l) Let us perform first the symmetrization over P ∈ S n . Since the first product in the left hand side of (A.12) is S n -invariant, and using the symmetrization identity ( [117, Chap. III, Eq.
Applying Lemma A.1 in (A.10), we obtain the formula for the overlap (3.17) given in the text.
We note that the overlap formula it is a priori valid before the insertion of the solution of the Bethe equations, i.e. it is valid for any set of complex λ j such that the overlap integral converges. The condition for that to be true can be read from (A.9) as Re(iλ j ) < B (A.14) for all j ∈ [1, n]. Inserting a string state, labeled by {k j , m j } j=1,...,n s , the condition becomes max j 1 2 (2m j − 1) ≤ 1 2 (2n − 1) < B, which leads to the condition n 2 < B + 1 2 .

B From Two-Dimensional Kernels to Scalar Kernels
We present in this section an equivalent representation of a class of Fredholm Pfaffians with 2 × 2 block kernels in terms of a Fredholm determinant with a scalar valued kernel. Consider a measure dμ on a contour C in the complex plane and another measure dν ς on the real line R, depending on a real parameter ς. Consider the quantity g(ς) defined by the series and Then, following Ref. [54], we have the following equivalent representations for g(ς).

Lemma B.1 (Fredholm Schur Pfaffian) g(ς) is equal to a Fredholm
Pfaffian with a 2 × 2 matrix valued skew-symmetric kernel For (r , r ) ∈ R 2 the matrix kernel K is given by which are assumed to be in L 2 (R + ), the scalar kernelK is given, for (x, y) ∈ R 2 + , bȳ K (x, y) = 2∂ x R dν ς (r ) [f even (x + r )f odd (r + y) − f odd (x + r )f even (r + y)] (B.7) and the scalar kernel I is the identity kernel I (x, y) = 1 x=y .

C.1 The Hastings-McLeod Solution of the Painlevé II Equation
Define q(s) to be the solution of the Painlevé II equation for s ∈ R, q (s) = sq(s) + 2q(s) 3 (C.1) which satisfies the asymptotic condition q(s) ∼ s→+∞ Ai(s). The unique smooth solution of the Painlevé II equation with that asymptotic condition is called the Hastings-McLeod solution [110]. Let us indicate the two following formula which we use in the manuscript. The first one is given in Ref. [

C.2 Relations Between the Tracy-Widom, Baik-Rains Distributions and the Hastings-McLeod Solution of the Painlevé II Equation
The Tracy-Widom distributions for β = 2, 1, 4 and the Baik-Rains distribution (denoted BR) are given by [109,118] -For β = 2 It allows to obtain the following relation between the distributions: and The next sections provide tail asymptotics for F 1 and F 2 .

D.1 Right Tail Using the Determinantal Formula
Let us perform the trace expansion for s → +∞ on the form

D.2 Right Tail Using the Asymptotics of q(s) and the Tracy-Widom Distributions
We can now compare with the formula (7.39) which we recast into Since q(s) behaves asymptotically for large positive s as the Airy function, we obtain at first order in q(s) the expansion Using the asymptotics of the Tracy-Widom distributions (C.14) and (C.15) and the asymptotics of the Airy function, we obtain

E Extended Kernel and the Kadomtsev-Petviashvili Equation
Recently it was shown that the height CDF at the (large time) KPZ fixed point for the full space problem is related to scale-invariant solutions of the Kadomtsev-Petviashvili (KP) equation [121]. A related observation was made for the periodic KPZ fixed point [122]. This connection to the KP equation extends to the generating function at arbitrary time for the KPZ equation in full space, for some particular initial conditions, droplet, half-Brownian [121] and Brownian [123]. A similar relation was also obtained for a class of linear statistics associated to the Airy process [123]. However at present no such relation is known for the half-space problem.
Here we provide an extended version of our kernel which can be related to the KP equation. It involves however an additional variable which plays the role of a "fictitious space". Although at this stage we could not see a physical interpretation for this variable, we believe this fact is curious enough to be reported.
We recall that for the half-space problem at all times, the generating function of the exponential of the KPZ height at the origin, x = 0, reads (7.3) E exp(−ς W e H (t) ) = Det(I −K t,ς ) L 2 (R + ) . (E.1) Denoting ς = e −r , one rewrites the kernelK as sin(π(z − w)) sin(π(z + w)) e −(x+r )z−(y+r )w+t w 3 +z We can now extend the kernel by introducing a fictitious variable u in the following way sin(π(z − w)) sin(π(z + w)) e −(x+r )z−(y+r )w+ u 2 (w 2 −z 2 )+t w 3 +z 3 3 (E.5) so that K t,e −r ,u=0 (x, y) =K t,e −r (x, y) (E.6) The new kernel K t,e −r ,u verifies the same set of differential equations (E.4) and in addition verifies a third one ∂ u K t,e −r ,u (x, y) = 1 2 [∂ 2 y − ∂ 2 x ]K t,e −r ,u (x, y) (E.7) It was shown in [121], and known earlier in [124] (see discussion in [123,Appendix E]), that the three conditions (E.4) and (E.7) imply that the Fredholm determinant associated to K is a τ -function of the KP equation. Thus, denoting F(r , t, u) = Det(I − K t,e −r ,u ) L 2 (R + ) , the function φ(r , t, u) := 2∂ 2 r log(F(r , t, u)) solves the KP equation for (r , t, u) Note that the knowledge of φ(r , t, u = 0) is equivalent to the knowledge of the half-space generating function for the KPZ equation. For the full-space KPZ equation, the variable u was interpreted as a spatial variable whereas in our case, it is a fictitious variable with no obvious direct interpretation.