Stationary measure for the open KPZ equation
Ivan Corwin

[Ivan Corwin] Department of Mathematics, Columbia University,
New York, NY, USA
E-mail<EMAIL_ADDRESS>
Alisa Knizel

[Alisa Knizel] Department of Statistics,
the University of Chicago,
Chicago, IL, USA.
E-mail<EMAIL_ADDRESS>
We provide the first construction of stationary measures for the open KPZ equation on the spatial interval $[0,1]$ with general inhomogeneous Neumann boundary conditions at $0$ and $1$ depending on real parameters $u$ and $v$, respectively. When $u+v\geq 0$ we uniquely characterize the constructed stationary measures through their multipoint Laplace transform which we prove is given in terms of a stochastic process that we call the continuous dual Hahn process.

§ INTRODUCTION

The open KPZ equation models stochastic interface growth on an interval $[0,1]$ subject to mixed inhomogeneous Neumann boundary conditions at $0$ and $1$. It is expected that for each choice of boundary parameters $u$ and $v$, there exists a unique stationary measure (i.e. the height increment process remains invariant in law with respect to the time evolution).

Theorem <ref> (<ref>)-(<ref>), shows that the open WASEP stationary measures are tight and all subsequential limits are stationary measures for the open KPZ equation. This relies on stochastic PDE methods and ASEP couplings. Theorem <ref> (<ref>), shows that if the boundary parameters satisfy $u+v\geq 0$ then there is a unique limit point for the corresponding WASEP stationary measures and that limit point can be completely characterized in terms of an expect formula for its multi-point Laplace transform which is written in terms of the continuous dual Hahn process.

To illustrate the type of formulas we find, if $H_{u,v}(X)$ has the law of the above described limit point open KPZ stationary measure with $u,v>0$ we show that for $c\in (0,2\vv)$,

\begin{equation}
\E\left[e^{-cH_{u,v}(1)}\right] = e^{c^2/4}\,\cdot\,\dfrac{\displaystyle\int_0^{\infty} e^{-\frac{1}{4}r}\cdot
\frac{\Big|\Gamma\Big(\frac{c}{2}+\uu+\iu\frac{\sqrt{r}}{2}\Big)\cdot\Gamma\Big(-\frac{c}{2}+\vv+\iu\frac{\sqrt{r}}{2}\Big)\Big|^2}{\sqrt{r}
\cdot \big|\Gamma(\iu \sqrt{r})\big|^2}dr}{\displaystyle\int_0^{\infty} e^{-\frac{1}{4}r}\cdot
\frac{\Big|\Gamma\Big(\uu+\iu\frac{\sqrt{r}}{2}\Big)\cdot\Gamma\Big(\vv+\iu\frac{\sqrt{r}}{2}\Big)\Big|^2}{\sqrt{r}
\cdot \big|\Gamma(\iu \sqrt{r})\big|^2}dr}
\end{equation}

Since we fix $H_{u,v}(0)=0$, the random variable $H_{u,v}(1)$ records the net height change for the open KPZ equation stationary measure across the full interval $[0,1]$. The above formula uniquely characterizes the law of this random variable. Our general result allows for $u+v>0$ without assuming both $u,v>0$, and provides formulas for the multi-point Laplace transform. If either $u$ or $v$ are negative, then the integral formulas become considerably more complicated and involve a mixture of continuous parts and discrete atomic parts. All of our formulas also involve restrictions on the $c$ variables in the Laplace transform. While it should be possible to analytically continue to formulas valid for all $c$ variables, we do not pursue this here and we expect that these analytic continuations will also involve considerably more atomic terms.

We now briefly recall the open KPZ equation and ASEP, and introduce the simplest case of the  continuous dual Hahn process in which there are no atoms. This special case is sufficient for the case where $u,v>0$. However, for $u+v>0$ with either $u$ or $v$ negative, we must deal with a more general version of this process with atoms in its support. This version is introduced later in Section <ref>. Armed with these processes (KPZ, ASEP and continuous dual Hahn) we then state our main result, Theorem <ref>.

§.§ The open KPZ equation

The stochastic heat equation (SHE) with inhomogeneous Robin boundary conditions parameterized by $u,v\in \R$ is formally written as

\begin{equation}\label{eq_SHE_informal}
\partial_T Z(T,X) = \tfrac{1}{2}\partial_X^2 Z(T,X)+  \xi(T,X) Z(T,X)
\end{equation}

for $T\geq 0$ and $X\in [0,1]$ with boundary conditions that for all $T>0$,

\begin{equation}\label{eq_SHE_bcs}
\partial_X Z(T,X)\big\vert_{X=0} =(u-\tfrac{1}{2})  Z(T,0),\qquad\textrm{and}\qquad  \partial_X Z(T,X)\big\vert_{X=1} = -(v-\tfrac{1}{2}) Z(T,0).
\end{equation}

Here $\xi(T,X)$ is space-time Gaussian white noise.
A proper definition and construction of (<ref>) was given in [15] and is recalled in Section <ref>. This require working with an integrated version of (<ref>) and (<ref>) since $Z(T,X)$ is not differentiable.

The Hopf-Cole solution to the open Kardar-Parisi-Zhang (KPZ) with inhomogeneous Neumann boundary conditions parameterized by  $u,v\in \R$ is defined as the logarithm of the SHE

\begin{equation}\label{eq_HC_KPZ}
H(T,X):=\log Z(T,X).
\end{equation}

We will generally drop the term “Hopf-Cole” and just refer to $H$ defined as above as the open KPZ equation.

At a formal level (which will not be used at all in this paper) $H(T,X)$ can be thought of as a solution to the following stochastic PDE (KPZ equation):

\begin{equation}\label{eq_KPZ_informal}
\partial_T H(T,X) = \tfrac{1}{2}\partial_X^2 H(T,X) + \tfrac{1}{2}\big(\partial_X H(T,X)\big)^2 +  \xi(T,X)
\end{equation}

where, for all $T>0$, we impose the boundary conditions

\begin{equation}\label{eq_KPZ_bcs}
\partial_X H(T,X)\big\vert_{X=0} = u,\qquad\textrm{and}\qquad  \partial_X H(T,X)\big\vert_{X=1} = -v.
\end{equation}

Note that making sense of this equation directly is significantly harder than making sense of the SHE.
We have removing a factor of $1/2$ in going from the SHE to the KPZ equation boundary condition. This is purely a matter of convention to match the convention for KPZ boundary conditions present in [30] and [33] which studies the open KPZ equation directly and not in terms of the Hopf-Cole transform. Again, since we only work with the Hopf-Cole transform, this is inconsequential to this work. Finally, note that if we let (again formally) $U(T,X) = \partial_X H(T,X)$, then $U$ will solve the stochastic Burgers equation with $U(T,0)=u$ and $U(T,1)=-v$. As we will later explain, when $u+v=0$, the open KPZ equation has Brownian motion stationary measure. The convention of dropping the $1/2$ factors is implemented so that in this case, formally $U(T,0)=U(T,1)$.

A stationary measure for the open KPZ equation is the law on a random function $H_0:[0,1]\to \R$ with $H_0(0)=0$ and with the property that if we run the open KPZ equation with initial data $H(0,X) = H_0(X)$ for all $X\in [0,1]$, the law of $X\mapsto H(T,X)-H(T,0)$, as a process in $X$, is $T$-independent for all $T\geq 0$. For the stochastic Burgers equation $U_0(X)=\partial_X H_0(X)$ will be an invariant measure in the usual sense and for the SHE started with $Z_0(X)=e^{H_0(X)}$, we will have that the law of  $X\mapsto Z(T,X)/Z(T,0)$ is $T$-independent.

This paper constructs stationary measures for the open KPZ equation with general parameters $u$ and $v$ as a limit of the stationary measure of the open ASEP height function process. When $u+v=0$ the stationary measure constructed in this manner is quite simple—a Brownian motion with drift $u=-v$. When $u+v\geq 0$ (the so called fan region), we are able to provide a complete description of the stationary measure in terms of Laplace transform formulas for its finite dimensional distributions. In fact, we prove a sort of duality formula between the stationary measure and a process that we call the continuous dual Hahn process. The fan region, $u+v\geq 0$, covers the entire limit of the ASEP maximal current region (when $u,v>0$), as well as portions of the low and high density regions (the lower density regime is when $u<v$ and $u<0$; the high density regime is when $u>v$ and $v<0$). We also find a remarkable stochastic ordering of the stationary measure increments with respect to the boundary parameters $u,v$. For the shock region, $u+v<0$, we do not have explicit formulas.

§.§ The open ASEP

We consider the open ASEP on $N$ sites $\{1,\ldots, N\}$ with left jump rate $q<1$ right jump rate $1$, and jump rates $\alpha,\beta,\gamma,\delta$ for particles entering at $1$, exiting at $N$, exiting at $1$ and entering at $N$, respectively. To this particle system we associate a height function process. See Figure <ref> for an illustration of the process and the various rates, or Section <ref> for a more detailed definition of the model. Importantly,  note that there is a unique stationary measure for this height process (or equivalently, a unique invariant measure for the open ASEP occupation process).

Open ASEP with system size $N=10$ and its height function  $h_N$ (the $\pm1$ sloped function drawn above the particle configuration). The red arrows indicate some possible moves with rates labeled.

In place of  $(q,\alpha,\beta,\gamma,\delta)$ we can parameterize open ASEP by $(q,A,B,C,D)$ where

\begin{equation}\label{eq_ABCD}
A=\kappa^+(q,\beta,\delta), \quad B=\kappa^-(q,\beta,\delta),\quad C=\kappa^+(q,\alpha,\gamma) ,\quad D=\kappa^-(q,\alpha,\gamma)
\end{equation}

\begin{equation}\label{eqnkappa}
\kappa^{\pm}(q,x,y):=\dfrac{1}{2
x}\left(1-q-x+y\pm\sqrt{(1-q-x+y)^2+4 xy} \right).
\end{equation}

It is easy to check that given $q$, (<ref>) provides a bijection between $\big\{(\alpha,\beta,\gamma,\delta):  \alpha, \beta>0,\, \gamma, \delta\geq 0\big\}$ and $\big\{(A,B,C,D): A, C>0,\, B,D\in (-1,0]\big\}$.

We make the following assumptions on the rates and how they vary with $N$:

\begin{equation}\label{eq:assumerates}
q=\exp\Big(-\frac{2}{\sqrt{N}}\Big), \qquad A=q^{\uu},\qquad B=-q,\qquad C=q^{\vv},\qquad D=-q.
\end{equation}

We may solve for $\alpha,\beta,\gamma$ and $\delta$ in terms of $q,u$ and $v$ as

\begin{equation}\label{eq:explicitparam}
\alpha = \frac{1}{1+q^u},\qquad \beta = \frac{1}{1+q^v},\qquad \gamma = \frac{q^{u+1}}{1+q^u},\qquad \delta = \frac{q^{v+1}}{1+q^v},
\end{equation}

thus under (<ref>), $q,u,v$ parameterize a family of open ASEPs. We define the following diffusive scaling of the open ASEP stationary measure height function correspond to $u$ and $v$:

\begin{equation}\label{eqn_hscale}
\hscale(X) := N^{-1/2} h_N(NX),
\end{equation}

where $h_N$ is the open ASEP height function (see (<ref>) for its definition or Figure <ref> for an illustration).
To keep track of the dependence on $u$ and $v$ we will sometimes write $H^{(N)}_{u,v}$.

The scaling of $q$ in (<ref>) in conjunction with the height function scaling in (<ref>) is called weak asymmetry scaling. The conditions on $A,B,C$ and $D$ in (<ref>) correspond to $\alpha,\beta,\gamma$ and $\delta$ which satisfy Liggett's condition that $\alpha +\gamma/q =1$ and $\beta+\delta/q=1$. Moreover, from (<ref>) we see that $\alpha,\beta,\gamma$ and $\delta$ satisfy triple point scaling, which means that as $N\to \infty$,

\begin{equation}
\label{eq_trip_exp}
\begin{aligned}
&\alpha = \frac{1}{2} + \frac{u}{2} N^{-1/2} + o(N^{-1/2}), &\beta = \frac{1}{2} + \frac{v}{2} N^{-1/2} +o(N^{-1/2}),\\
&\gamma = \frac{1}{2} - \frac{u}{2} N^{-1/2} + o(N^{-1/2}), &\delta = \frac{1}{2} - \frac{v}{2} N^{-1/2} + o(N^{-1/2}).
\end{aligned}
\end{equation}

§.§ Continuous dual Hahn process

The open KPZ stationary measures that we construct in this work will be characterized via a form of duality with another stochastic process which we call the continuous dual Hahn process (denoted below by $\T_s$). This is a special limit of the Askey-Wilson processes constructed by Bryc and Wesołowski [13]; see Section <ref> for more on this. The continuous dual Hahn process depends on two parameters $u,v\in\R$ which are assumed throughout to satisfy the relation $u+v>0$.

The definition of this process is simplest when $u,v>0$ and thus for the sake of this introduction we will only define it in that case here. Under this restriction on parameters, we will be able to work entirely with absolutely continuous distributions. If we allow for either $u$ or $v$ to be negative, we must introduce atomic parts to our measures and the complexity of formulas becomes considerable. We do address these complicated cases in Section <ref>. When $u,v>0$, the result in Theorem <ref> (<ref>) can be made sense of based on the definitions below (which also appears in [10]); when either $u$ or $v$ become negative one needs to turn to the definitions in Section <ref> which involves atoms (and is not present in [10]).

We will only define the continuous dual Hahn process $\T_s$ for $s\in [0,\Cuv)$ where

\begin{equation}\label{eq:Cuv}
\Cuv :=
\begin{cases}
2&\textrm{if }\vv\leq 0\textrm{ or }\vv\geq 1,\\
2\vv&\textrm{if }\vv\in(0,1).
\end{cases}
\end{equation}

The reason for this additional restriction also stems from the fact that formulas become considerably more involved outside of this window. This domain of definition will be sufficient for our purposes.

Assuming that $u,v>0$ and $s\in [0,\Cuv)$ define a measure $\p_s$ with density given by

\begin{equation}
\p_s(r):=\dfrac{(\uu+\vv)(\uu+\vv+1)}{8\pi}\cdot
\frac{\Big|\Gamma\Big(\dfrac{s}{2}+\uu+\iu\dfrac{\sqrt{r}}{2}\Big)\cdot\Gamma\Big(-\dfrac{s}{2}+\vv+\iu\dfrac{\sqrt{r}}{2}\Big)\Big|^2}{\sqrt{r}
\cdot \big|\Gamma(\iu \sqrt{r})\big|^2}{\bf 1}_{r>0}.
\end{equation}

This infinite mass measure will be the marginal distribution of our Markov process.

Following [51], we define the orthogonality probability measure for the continuous dual Hahn orthogonal polynomials as follows: For $a\in \R$ and $b=\bar c\in \C\setminus \R$ with $\textup{Re}(b)=\textup{Re}(c)>0$ let

\begin{equation}
\CDHtransition(x; a, b, c):=\dfrac{1}{8\pi}\,\cdot\,\dfrac{\bigg\rvert \Gamma\left(a+\iu  \dfrac{\sqrt x}{2}\right)\cdot\Gamma\left(b+\iu
\dfrac{\sqrt x}{2}\right)\cdot\Gamma\left(c+\iu  \dfrac{\sqrt x}{2}\right)\bigg \rvert^2}{\Gamma(a+b)\cdot \Gamma(a+c)\cdot\Gamma (b+c)\cdot \sqrt x\cdot\left\rvert \Gamma\left(\iu \sqrt x\right)\right\rvert^2}{\bf 1}_{x> 0}
\end{equation}

denote the density in $x$ of the measure $\CDHtransition(\cdot; a, b, c)$.

For $s,t\in [0,\Cuv)$ with $s<t$ and $m,r\in (0,\infty)$, define a measure $\p_{s,t}(m,\cdot)$ with density  given by

\begin{equation}
\p_{s, t}(m, r):=\CDHtransition\left(r; \vv-\dfrac{t}{2}, \dfrac{t-s}{2}+\iu\dfrac{\sqrt m}{2}, \dfrac{t-s}{2}-\iu\dfrac{\sqrt m}{2}\right).
\end{equation}

The continuous dual Hahn process is the Markov process $\{\T_s\}_{s\in [0,\Cuv)}$ with distribution at time $0$ given by $\p_0$ and transition probability given by $\p_{s,t}$. Lemma <ref> verifies that the $\p_{s,t}$ satisfy the Chapman-Kolmogorov equation, and that the marginal law of $\T_s$ is $\p_s$ for all $s\in [0,\Cuv)$.

§.§ Statement of the main result

For $d\in \Z_{\geq 1}$ we will assume throughout the paper that

\begin{align}
\vec{X} &= (X_0,\ldots, X_{d+1})\textrm{ where }0=X_0<X_1<\cdots<X_d\leq X_{d+1}=1,\\
\vec{c} &= (c_1,\ldots, c_d)\textrm{ where }c_1,\ldots, c_d>0,\label{eq:Xcs}\\
\vec{s} &= (s_1> \cdots>s_{d+1})\textrm{ where }s_k = c_k+\cdots +c_d \textrm{ and }s_{d+1} = 0.
\end{align}

We are now prepared to state our main theorem. Part (<ref>) shows the existence of open KPZ stationary measures as limits of WASEP stationary measures. Part (<ref>) records a coupling between these stationary measures in which height differences are stochastically ordered relative to a certain ordering of the boundary parameters. Part (<ref>) records the simple Brownian case for the stationary measure which occurs when $u+v=0$. Part (<ref>) records a duality which comes from switching the roles of $u$ and $v$. Part (<ref>) provides a unique characterization of the WASEP-stationary measure for the open KPZ equation provided that $u+v>0$ (the $u+v=0$ case was already addressed in part (<ref>)). This characterization is given by a remarkable duality formula which relates the stationary measure to the tangent Askey-Wilson process introduced in Section <ref>. This shows that provided $u+v>0$, there is a unique limit point in part (<ref>).

It should be noted that our result does not imply that there is a unique stationary measure for the open KPZ equation—rather it just implies that if such an alternative stationary measure exists, it does not arise as the limit of the open ASEP stationary measure. This is good reason to believe that the stationary measure that we construct is, in fact, unique—Conjecture <ref> contains such a statement and the paragraphs following it identify the evidence for this form of uniqueness.

Assume that the open ASEP parameters $q,\alpha,\beta,\gamma,\delta$ and $N$ satisfy Assumption <ref>.

*  Tightness and construction of WASEP-stationary measures:  For any $u,v\in \R$, the $N$-indexed sequence of laws of $H^{(N)}_{u,v}(\cdot)$ are tight in the space of measures on continuous functions on $[0,1]$. All subsequential limits $H_{u,v}$ are stationary measures for the open KPZ equation and are Hölder $1/2-$ almost surely. We call any such subsequential limit a WASEP-stationary measure for the open KPZ equation.

*  Coupling: For any $u,v\in \R$ the increments of the WASEP-stationary measures are stochastically ordered as follows. For any $M\in \Z_{\geq 2}$, $u_1\leq \cdots \leq u_M$ and $v_1\geq \cdots \geq v_M$, let $\left\{H_{u_i,v_i}\right\}_{i=1}^{M}$ be a collection of WASEP-stationary measures for the open KPZ equation with boundary parameters $u_i$ and $v_i$. Assume that these measures all arise in part (<ref>) along the same subsequence of values of $N$. Then there exists a coupling (i.e. a common probability space supporting all of these random functions) with the property that for all $0\leq X\leq X'\leq 1$ and $1\leq i\leq j\leq M$,

\begin{equation}\label{eqnHcoupling}
H_{u_i,v_i}(X')-H_{u_i,v_i}(X)\leq H_{u_j,v_j}(X')-H_{u_j,v_j}(X).
\end{equation}

*  Brownian case: For $u+v=0$, there is a unique WASEP-stationary measure  $H_{u,v}$ for the open KPZ equation whose law is that of a standard Brownian motion of drift $u=-v$.

*  Duality: For any $u,v\in \R$ let $H_{u,v}$ and $H_{v,u}$ be a pair of WASEP-stationary measures for the open KPZ equation which arise in part (<ref>) along the same subsequence of values of $N$. Then as processes in $X\in [0,1]$, $H_{u,v}(X)$ and $H_{v,u}(1-X)-H_{v,u}(1)$ have the same law.

* Explicit characterization via Laplace transform formula: For $u,v\in \R$ with $u+v>0$, there is a unique WASEP-stationary measure $H_{u,v}$ for the open KPZ equation (recall from part (<ref>)) whose law is determined by the following formula for its multi-point  Laplace transform. For any $d\in \Z_{\geq 1}$, and $\vec{X},\vec{c}$ and $\vec{s}$ as in (<ref>), provided that $c_1+\cdots +c_d =s_1<\Cuv$, see (<ref>),

\begin{align}\label{eq_mainthmformphi}
\displaystyle\E\left[\exp\left(-\sum_{k=1}^{d} c_k H_{u,v}(X_k)\right)\right] = \frac{\displaystyle\E\left[\exp\left(\frac{1}{4}\sum_{k=1}^{d+1} (s_k^2-\T_{s_k})(X_k-X_{k-1})\right)\right]}{\displaystyle\E\left[\exp\left(-\frac{1}{4}\T_0\right)\right]} =:\phi_{u,v}(\vec{c},\vec{X})\qquad
\end{align}

where $\T_s$ is the continuous dual Hahn process (see Section <ref> or Section <ref>).

§.§.§ Outline

Section <ref> contains a discussion of Theorem <ref> and highlights key ingredients in its proof.

The proofs of Theorem <ref> (<ref>)-(<ref>) are given in Section <ref> and rely on Section <ref> (which introduces the open ASEP and gives some background about its phase diagram and behavior), Section <ref> (which describes the weak asymmetry scaling under which open ASEP converges to the open KPZ equation) and Section <ref> (which describes certain coupling results for open ASEP).

The proof of Theorem <ref> (<ref>) is given in Section <ref>. The starting point in this proof (given in Section <ref>) is Corollary <ref> (see also Proposition <ref>  and <cit.>) which relates the generating function for the open ASEP invariant measure to the Askey-Wilson process. Section <ref> defines the continuous dual Hahn process which arises as a special limit of the Askey-Wilson process. The main calculation in the proof of Theorem <ref> (<ref>) is Proposition <ref> which computes the limit of the open ASEP invariant measure generating function. Combining this with some of the results in Theorem <ref> (<ref>)-(<ref>) we show weak convergence and that the limit of the open ASEP formula gives the open KPZ Laplace transform formula claimed in the theorem. The proof of this limit is given in Section <ref> and relies heavily on precise asymptotics for $q$-Pochhammer symbols. These asymptotics are stated as Proposition <ref> and proven in Section <ref>.

§.§.§ Notation

We adopt the notational device that lower case variables will always refer to discrete objects while upper case variables will always refer to continuum objects, or (when a superscript $(N)$ is included) discrete objects rescaled in such a way as to converge to continuum objects. For instance $h_{N}(t,x)$ denotes the height function of the open ASEP with $N$ sites at time $t\geq 0$ and position $x\in [0,N]$, while $H^{(N)}(T,X)$ denotes the scaled and centered height function given in (<ref>) where $X\in [0,1]$ corresponds to $x=NX$ and $T\geq 0$ corresponds to $t=N^2T$ (there is also scaling of $h_N$ by $N^{-1/2}$). For integers $a\leq b$ we write $\llbracket a,b\rrbracket$ for the set $\{a,\ldots, b\}$. We use the notation $\Z_{\geq a}$ to denote all integers greater than or equal to $a$ (and likewise for $\Z$ replaced by $\R$ and $\geq$ replaced by $>$, $\leq $ or $<$). In a nod to the common notation in the physics literature, we will use the notation $\langle \cdot \rangle_N$ to denote the expectation of a function of the open ASEP invariant measure. For the open KPZ equation we will utilize the more traditional $\E[\cdot]$ notation, and likewise for the Askey-Wilson process. Various measures that we will encounter will involve absolutely continuous parts and discrete atomic parts. We will use the convention of a superscript $c$ for the density of the continuous part and superscript $d$ for the mass function for the discrete part.

We will use the standard notation for the Pochhammer and $q$-Pochhhamer symbols: For $j\in \Z_{\geq 0}$ and $x\in \R$, define $[x]_j:=(x)(x+1)\cdots(x+j-1)$ with the convention that $[x]_0:=1$. For multiple arguments $x_1,\ldots x_n\in \R$, define $[x_1, \ldots, x_n]_j:=[x_1]_j\cdots [x_n]_j$. For $a, q\in \C$, with $|q|<1$, and $j\in \Z_{\geq 0}\cup\{\infty\}$, define $(a; q)_j:=(1-a)(1-a q)\cdots(1-a q^{j-1})$ and $(a_1, \dots, a_n; q)_j=(a_1;q)_j\cdots (a_j;q)_j$. We will often omit the dependence on $q$ and write $(a)_j$ or $(a_1,\dots,a_n)_j$. We will denote $\iu := \sqrt{-1}$.

§.§.§ Acknowledgements

I. Corwin was supported by the NSF through DMS:1811143, DMS:1664650 and DMS:1937254, a Packard Fellowship in Science and Engineering and a W. M. Keck Foundation Science and Engineering Grant.  A. Knizel was supported by the NSF through DMS:1704186.
The authors wish to thank Wlodek Bryc and Yizao Wang for discussions surrounding Askey-Wilson processes and their tangent processes, and in particular to Wlodek Bryc for visiting Columbia and inspiring us to work on this subject and Yizao Wang for illuminating discussions around the atoms which occur in our work (in particular, see Remark <ref>).
The authors also wish to thank Alex Dunlap, Alexandra Florea, Hao Shen and Kevin Yang for helpful discussions and comments related to this work.

§ DISCUSSION AND KEY IDEAS FOR THEOREM <REF>

We collect remarks about the relation between Theorem <ref> and existing literature and possible directions of further interest. Some related discussion on open ASEP is contained in Section <ref>.

§.§ Discussion

§.§.§ Uniqueness

Our tightness result implies existence of open KPZ stationary measures. Uniqueness is an open problem. Even for parameters where we prove that ASEP-stationary measures are unique (i.e. uniqueness of the limit points of the scaled open ASEP stationary measures), we do not rule out the existence of other stationary measures for the open KPZ equation with the same boundary parameters. However, based on related results in the literature, we conjecture that the stationary measures we have constructed are unique for all choices of $u$ and $v$.

Consider any two random functions $H_0,\tilde H_0\in C([0,1])$, supported on the same probability space. On this probability space, define a space-time white noise $\xi(T,X)$ with the time variable $T\in \R$ and let $H(T,X;-T_0)$ and $\tilde H(T,X;-T_0)$ denote the solution to the open KPZ equation started at time $-T_0$ with (respectively) initial data $H_0$ and $\tilde H_0$. Then the following one force one solution principle holds: Almost surely as $T_0\to \infty$, the random functions $(T,X)\mapsto H(T,X)-H(T,0)$ and $(T,X)\mapsto \tilde H(T,X)-\tilde H(T,0)$ converge to the same limit on any compact region of $(T,X)$. In particular, there exists a unique stationary measure for the open KPZ equation.

In the context of the KPZ equation with periodic boundary conditions on the torus, [35] showed uniqueness of  Brownian bridge  as the stationary measure while [32] constructed the infinitesimal generator for this version of the KPZ equation and estimated its spectral gap, establishing $L^2$ exponential ergodicity. Further, [48] demonstrated the one force one solution principle. In much earlier work, [28] demonstrated the one force one solution principle for the Burgers equation with compact stochastic forcing, and [6] extended that to the non-compact full line setting. Let us also mention recent work of [23] which deals with the stochastic Burgers equation on the full line with colored noise (i.e., space-time white noise mollified in space) and proves existence and uniqueness of stationary solutions for each choice of specified mean.

The conjecture relates to the mixing properties of the open KPZ equation. For open ASEP this problem was studied in terms of the relaxation time (inverse spectral gap) in the physics literature via Bethe ansatz in work of [22] which (non-rigorously) showed that in the maximal current phase (see Section <ref>), the relaxing time is of order $N^{3/2}$, whereas in the other phases, it is of order $N$. Under our rate assumptions we are around the triple-point between the open ASEP phases and hence the same time scales as in the maximal current phase should hold. With the weak asymmetry of order $N^{-1/2}$, it is reasonable to then expect that mixing should occur on time of order $N^{-2}$ (we need to compensate for the weak asymmetry). This is precisely the scaling of time under which the open ASEP converges to the open KPZ equation (Proposition <ref>) and hence agrees with the conjecture. The only rigorous progress towards mixing for the open ASEP is in recent work of [31]. While they were able to prove the order $N$ mixing time in the high and low density phases, in the maximal current phase they were able to upper-bound the mixing time by $N^3 \log(N)$, not the conjectured $N^{3/2}$ scale. Note $\alpha=\beta=\gamma=\delta=0$ (i.e. reflecting boundary conditions), the order $N^{3/2}$ mixing time and exact cutoff profile for ASEP on a segment has been established in [39] and, respectively, [9].

§.§.§ Coupling and Brownian comparison

Theorem <ref> (<ref>) and (<ref>) can be combined to show that increments of the $H_{u,v}$ stationary measures are stochastically sandwiched between Brownian motions of different drifts. Take $M=3$ and let $u_1=-v, u_2=u, u_3=u$ and $v_1=-u, v_2=v, v_3=v$. Then Theorem <ref> (<ref>) implies that along any subsequence $\hscale_{u_1,v_1}$ converges in distribution to  $H_{u_1,v_1}$ which has the law of a drift $-v$ Brownian motion, denoted by $B_{-v}$; similarly $\hscale_{u_3,v_3}$ converges in distribution to  $H_{u_3,v_3}$ which has the law of a drift $u$ Brownian motion, denoted by $B_u$. There is a subsequence along which $\hscale_{u_2,v_2}$ converges to $H_{u_2,v_2}$. Thus,  Theorem <ref> (<ref>)  implies that for all $0\leq X\leq X'\leq 1$,

\begin{equation}\label{eq:browniancomparison}
B_{-v}(X')-B_{-v}(X) \leq H_{u,v}(X')-H_{u,v}(X) \leq B_{u}(X')-B_{u}(X),
\end{equation}

hence moments of increments of $H_{u,v}$ are bounded by moments of increments of Brownian motions.

The Brownian case in Theorem <ref> follows easily from the known fact that Bernoulli product measure is the invariant measure for open ASEP. This is the only case when the invariant measure for open ASEP is simple and of product form. For general parameters, the invariant measure is quite complicated, hence our use of Askey-Wilson processes. Note that for ASEP on the full line or torus, the invariant measure is product Bernoulli (in fact, on the torus, since open ASEP conserves the number of particles, the uniform law on particle configurations given a total number of particles is also invariant). For full line [29] and period [35] KPZ, the stationary measure is Brownian (two-sided Brownian motion with general drift, or Brownian bridge with any fixed height shift).

Our proof of the coupling in Theorem <ref> (<ref>) proceeds from a similar result for our WASEP approximation. However, we expect that it should be possible to prove such a statement directly in the continuum. On the line, such a proof is contained in <cit.>.

§.§.§ Range of the Laplace transform formula

The Laplace transform formula Theorem <ref> (<ref>) requires that $c_1,\ldots, c_d>0$ satisfies $c_1+\cdots +c_d=s_1 <\Cuv$ (a strictly positive constant depending only on $v$ as in (<ref>)). We expect that it should be possible to compute the Laplace transform for all $c_1,\ldots, c_d>0$ without this restriction on their sum. This would involve extending the definition of the continuous dual Hahn process $\T$ as well the convergence result of Proposition <ref> (see Remark <ref>). The Laplace transform formula we prove holds for an open set of parameters which by uniqueness of analytic continuations is sufficient to characterize  the stationary measure.

As indicated by Corollary <ref> and Remark <ref>, it should be possible to compute a slightly more general Laplace transform formula than in (<ref>), namely
$\E\left[\exp\left(-\sum_{k=1}^{d} c_k H_{u,v}(X_k) -\tilde{c} H_{u,v}(1)\right)\right]$ with $\tilde{c}\in \R$ and the $c_1,\ldots,c_d$ are above. We do not, however, know how to prove a similar formula when any of the $c_1,\ldots, c_d$ are negative. This is despite the fact that the Laplace transform is finite and analytic in those cases (this follows from the Brownian sandwiching). Therefore, in principle, one could access these more general formulas by analytic continuation of (<ref>). We do not pursue this here, though remark that analytically continuing a formula defined in terms of complex integrals is not always straight-forward. A simple example is to consider $f(z) = \int_{-\infty}^{\infty} e^{-t^2}(z-t)^{-1}dt$ for $\textup{Im}(z)>0$. Its analytic continuation to $\textup{Im}(z)\leq 0$ is given by $f(z) = \int_{-\infty}^{\infty} e^{-t^2}(z-t)^{-1}dt + 2\pi \iu e^{-z^2}$.

§.§.§ Probabilistic description

Theorem <ref> (<ref>) uniquely characterizes the finite dimensional distributions of the WASEP-stationary measure for the open KPZ equation when $u+v>0$. However, unlike when $u+v=0$, we do not have a simple probabilistic description of this height function process. As explained in Proposition <ref>, for ASEP with $q$ is fixed (not scaled to $1$ as we consider), the limit of the stationary height function measure can be described in terms of Brownian motion plus independent excursions or meanders (depending on the $u,v$ parameters). This is shown in [14].

In our case, we can show that the law of $X\mapsto H_{u,v}(X)$, as a process in $C([0,1],\R)$, is equal to that of $X\mapsto 2^{-1/2} B(X) + \tilde{H}_{u,v}(X)$ where $B$ is a standard Brownian motion and $\tilde{H}_{u,v}$ is independent of $B$, has  $\tilde{H}_{u,v}(0)=0$, and has finite dimensional distributions uniquely characterized by (assume the same conditions on $d, \vec{X}, \vec{s}$ as in (<ref>))

\begin{equation}\label{eq_mainthmform}
\E\left[e^{-\sum_{k=1}^{d}c_k \tilde{H}_{u,v}(X_k)}\right] =
\frac{\E\left[e^{-\frac{1}{4}\sum_{k=1}^{d+1}\T_{s_k}(X_k-X_{k-1})}\right]}{\E\left[e^{-\frac{1}{4}\T_{0}}\right]}.
\end{equation}

This equation can be reformulated in the following compelling form:

\begin{equation}
\E\left[e^{-\int_0^1 \tilde{H}_{u,v}(\chi(r))\sigma(dr)}\right] =\frac{\E\left[e^{-\frac{1}{4}\int_0^1 \T_{\sigma(r)}\chi(dr)}\right]}{\E\left[e^{-\frac{1}{4}\T_{0}}\right]}
\end{equation}

where we assume that $\chi(r)$ is a piecewise constant monotone increasing right-continuous function with $\chi(0)=0$ and $\chi(1)=1$, and that $\sigma$ is a piecewise constant monotone decreasing  right-continuous function with $\sigma(0)<\Cuv$ and $\sigma(1)=0$.
It is natural to expect that the piecewise constant condition can be relaxed. Proving this would require controlling the regularity of both $\tilde{H}_{u,v}$ and $\T$—the first of which is controlled by the Brownian comparison in Theorem <ref> and the second of which may be accessible through analysis of the exact transition probability formulas.

Another natural question about our Laplace transform formula is whether we can show directly using it that the associated stochastic process $H_{u,v}$ is stationary for the open KPZ equation. For many reasons this is very challenging. We would really need to work at the level of the derivative, i.e., the stochastic Burgers equation. The generator for the periodic version of that equation was recently constructed in [32], though it is unclear if those techniques would work well with the open boundary conditions. If that were doable, then a probabilistic description of $H_{u,v}$ could be quite useful in checking that the generator applied to that law gives zero.

§.§.§ Shock behavior

Our Laplace transform formula characterization in part  (<ref>) of Theorem <ref> is only shown for $u+v\geq 0$. This is known as the fan region in open ASEP (see Section <ref>) and it is in this region of parameters that Askey-Wilson processes can be used to characterize the open ASEP invariant measure. The complementary shock region is when $u+v<0$. We do not have exact formulas in this region (see also Remark <ref>); our results are limited to tightness and coupling.

§.§.§ Half space

With our methods, it should be possible to access the stationary measure for the open KPZ equation on an interval $[0,M]$ for any $M>0$. Taking the $M\to \infty$ limit should make contact with the KPZ equation in a half-space. The half-space model should have a rich family of stationary measures indexed by the boundary parameters at 0 and the slope at infinity. It would be interested to construct and explicitly characterize them. Just as with open KPZ, there is a special family of boundary conditions for half-space KPZ in which the stationary measure is Brownian. This has been studied in [8]. The study of the half ASEP invariant measure goes back to  [37]. The matrix product ansatz was developed for this model in [34].

§.§.§ Continuous dual Hahn process

Our formula in Theorem <ref> (<ref>) is given in terms of the continuous dual Hahn process. In Section <ref> we provide a definition of this process provide its time parameter is limited to lie in the non-empty interval $[0,\Cuv)$. The process arises as a tangent process to the Askey-Wilson process through rather involved asymptotic analysis. We restricted the time interval on which we perform these asymptotics since it reduces the amount of bookkeeping needed to keep track of all of the possible atoms in the marginal and transition measures. This restriction does not impede our proof Theorem <ref> (<ref>), though it does limit the range of the Laplace variables for which we have an explicit formula. We do not touch herein the compelling problem of defining the continuous dual Hahn process on all of $\R$ and showing that it can be used to extend the Laplace transform formula to a more general range of parameters. Another interesting question is to construct nice versions of the continuous dual Hahn process, e.g. those with piece-wise continuous trajectories which only jump between atoms and are otherwise continuous.

§.§ Key ideas in proving Theorem <ref> (<ref>)-(<ref>)

All of the key ideas rely on stochastic analysis:

* When $u+v=0$, the open ASEP invariant measure is product Bernoulli. This implies tightness of the height function $\hscale_{u,v}$ and Hölder bounds for it, and yields Theorem <ref> (<ref>).

* There exists an attractive coupling between versions of open ASEP with different boundary rates. This coupling (see Section <ref>) tells us that the finite $N$ version of the coupling in Theorem <ref> (<ref>) holds. We then bound the increments of the height function by random walk increments using the $u+v=0$ result above and deduce similar Hölder bounds. The attractive coupling that we provide is proved by use of a multi-species version of open ASEP. A simpler particle/hole duality for open ASEP implies the prelimiting version of Theorem <ref> (<ref>).

* The open ASEP converges to the open KPZ equation under suitable weak asymmetry and triple point scaling. This result is basically contained in [15, 45] (see Section <ref>) though to apply it to the open ASEP height function stationary measure, we need to verify Hölder bounds on the exponential of $\hscale_{u,v}(X)$. These can be deduced from the coupling and random walk bounds in increments discussed above. Thus, we show that the open ASEP height function process started from the stationary measure is tight as a space time process and any subsequential limit will be a Hopf-Cole solution to the open KPZ equation. This implies Theorem <ref> (<ref>)—that subsequential limits of the open ASEP stationary measure are stationary measures for the open KPZ equation.

§.§ Sketch of the ideas in proving Theorem <ref> (<ref>)

There are three (other) key ideas here, the first two of which rely on certain integrable structure and exact formulas:

* The generating function for the open ASEP invariant measure can be written in terms expectations over the Askey-Wilson process.
In fact, the true starting point is the matrix product ansatz (see Section <ref>) which describes the invariant measure in terms of certain non-commuting matrix products. This expression has a very algebraic flavor and for asymptotic problems it is natural to look for more analytic formulas. This was accomplished by Bryc and Wesołowski [13] who rephrased the matrix product ansatz in terms of the Askey-Wilson processes. We record their main result as Proposition <ref> in Section <ref> (see Section <ref> for definitions). Corollary <ref>, shows how this proposition yields a Laplace transform formula for the finite dimensional distributions of the open ASEP height function stationary measure. The major technical challenge in our proof is to extract the asymptotic behavior of this formula under the scaling which yields the open KPZ equation—this is accomplished in Proposition <ref>. The starting point for its proof is the fact that the Askey-Wilson process marginal distributions and transition probabilities are written in terms of the Askey-Wilson orthogonality measure, which in turn is written in terms of $q$-gamma functions (i.e., certain $q$-Pochhammer symbols).

* The $q$-gamma function has an explicit asymptotic expansion as $q\to 1$ with precise error bounds which can be controlled uniformly over all arguments. Point-wise $q$-gamma function asymptotics are well-known in the literature. However, our formulas involve taking integrals against these measures and hence we must show that we can apply the dominated convergence theorem. To do this, we need to control the $q$-gamma asymptotics as the argument varies in imaginary strips (i.e., $z$ with arbitrary imaginary part, and real part bounded in some interval). Proposition <ref> achieves this control (in fact we only need to use the simplest case when $m=1$, though we include an all order expansion as we believe it may be of independent interest). The proof of this result relies on methods from analytic number theory and along the way we also derive new uniform bounds on the Hurwitz zeta function (see Proposition <ref>). Even with the $q$-gamma convergence result from Proposition <ref>, there is considerable work in extracting well-controlled asymptotics of the Askey-Wilson process marginal and transition probabilities—this is accomplished in Section <ref>.

* Our Laplace transform convergence on any open interval implies weak convergence. Proposition <ref> proves a limit for the Laplace transform of the finite dimensional distributions for the open ASEP height function stationary measure. However, due to limitations on the parameters in Proposition <ref>, the Laplace transform variables $c_1,\ldots,c_d$ are restricted to lie in an interval $(0,\Cduv)$ for some $\Cduv>0$ (see Remark <ref>). This Laplace transform convergence along with some inputs from  Theorem <ref> (<ref>)-(<ref>) turns out to be sufficient to show weak convergence and that the Laplace transform of the limit is the limit of the Laplace transforms.

§.§.§ Matrix product ansatz

Since open ASEP is a finite state Markov chain, for any $N$ it has a unique invariant measure which can be found from its transition matrix. While for small $N$ this may be feasible, it becomes useless when it comes to taking a large $N$ limit—that is, unless there is a representation for the invariant measure which behaves in a controllable manner as $N$ grows. The first step in this direction is the matrix product ansatz which goes back to Derrida-Evans-Hakim-Pasquier [20]. They showed that the invariant measure $\pi^{\mathrm{ASEP}}_N$ (see Section <ref> for definitions) can be written as (here $\tau_1,\ldots, \tau_N$ are the occupation variables taking values in $\{0,1\}$)

\pi^{\mathrm{ASEP}}_N(\tau_1,\ldots, \tau_N) = \frac{\big\langle W\big\vert \big(\tau_1 D + (1-\tau_1)E\big)\cdots \big(\tau_N D + (1-\tau_N)E\big)\big\vert V\big\rangle}{\big\langle W\big\vert \big(D + E\big)^N\big\vert V\big\rangle}

provided that $\langle W|, |V\rangle, D$ and $E$ satisfy the relations

DE-qED=D+E,\qquad (\beta D-\delta E)|V\rangle = |V\rangle,\qquad \langle W|(\alpha E-\gamma D)= \langle W|.

This formula is verified by appealing directly to the Kolmogorov forward equation. It is not a priori clear how useful this formula is. For instance, do there even exist representations of the $\langle W|, |V\rangle, D$ and $E$ which satisfy the desired relation? Over time various (all infinite dimensional) representations for $\langle W|, |V\rangle, D$ and $E$ have been discovered. Some of these actually allow for calculation of various asymptotics. For $q=0$ such representations and calculations were developed in [20]. Later, Uchiyama-Sasamoto-Wadati [50] determined useful representations for the general $q>0$ case. From these representations, it was possible to compute asymptotics needed to demonstrate the phase diagram from Section <ref>. The $q=0$ version of Proposition <ref> was demonstrated by Derrida-Enaud-Lebowitz [21] in this manner as well. Based on a different matrix representation (a solution to a slightly more general version of the matrix product ansatz), Corteel-Williams [17] developed a combinatorial description for the open ASEP invariant measure in terms of tableaux.

The Askey-Wilson process formulas that we recall now are another implication of these representations, though of a more analytic flavor. This turns out to be very useful for asymptotics.

§.§.§ Askey-Wilson processes

Askey-Wilson polynomials form a remarkable family of $q$-orthogonal polynomials specified by four free parameters $a,b,c,d\in \C$ in addition to $q$. They resides on the top of the hierarchy of one-variable $q$-orthogonal polynomials in the Askey scheme [36]. These polynomials (in fact, their Jacobi matrix) were used to construct the infinite dimensional representation for $\langle W|, |V\rangle, D$ and $E$ in  [50].

In an a priori different direction,  Bryc and Wesołowski [12] used the orthogonality measure for these polynomials to construct a class of Markov processes called Askey-Wilson processes which enjoy the characterizing property of being quadratic harnesses. The Askey-Wilson polynomials are orthogonal martingale polynomials for these processes (just like how Hermite polynomials serve this role for Brownian motion, a special limit of the Askey-Wilson processes). We will discuss the construction in detail in Section <ref>.

These two perspectives were united by Bryc and Wesołowski [13] who showed the following remarkable identity that serves as the starting point for our analysis.

Let $ \big \langle \cdot  \big\rangle_N$
denote the expectation with respect to the invariant measure $\pi^{\mathrm{ASEP}}_N$ of open ASEP parameterized by $q$ and $(A,B,C,D)$ as in (<ref>).
Assume that $AC<1$. Then for $0 < t_1 \leq t_2 \leq\cdots\leq t_n$, the joint generating function of the invariant measure for open ASEP can be expressed as

\begin{equation}\label{eq:key_identity}
\bigg \langle\prod_{j=1}^N t_j^{\tau_j}  \bigg
\rangle _N=\dfrac{\E\left[\prod_{j=1}^N\left(1+t_j+2\sqrt t_j
\Y_{t_j}\right) \right]}{2^N \E
\left [\left (1+\Y_1 \right )^N\right ]},
\end{equation}

where $\{\Y_t\}_{t\geq 0}$ is the Askey-Wilson process with parameters $(A, B, C, D, q)$, see Section <ref>.

§.§.§ Asymptotics

The Askey-Wilson orthogonality measure is written in terms of $q$-Pochhammer symbols. Recall the q-Pochhammer symbol notation for $a, q\in \C,$ $|q|<1$

$$(a; q)_j:=(1-a)(1-a q)\dots\left(1-a q^{j-1}\right), (a_1, \dots, a_n; q)_j=(a_1;q)_j\dots (a_j;q)_j, \quad 1\leq j\leq \infty, n \in \Z_{\geq 1}.$$

We will generally omit the dependence on $q$  and write $(a)_j$ or $(a_1,\dots,a_n)_j$ if it is clear from the context.
One of the key technical ingredients in our study is the asymptotics of the $q$-Pochhammer symbols as $q\rightarrow 1$ which we present next. After the statement of the result, we will briefly discuss how it relates to some existing literature, in particular asymptotics of the $q$-Gamma function.

Let us define for $z\in \C$ and $\kappa >0$ the following functions:

\begin{align}
\AqP^+[\kappa, z]&= -\dfrac{\pi^2}{6 \kappa}-\left(z-\dfrac{1}{2}\right)\log \kappa- \log \left[\dfrac{\Gamma (z)}{\sqrt{2\pi}}\right ],\label{eq:AqPp}\\
\AqP^-[\kappa, z]&= \dfrac{\pi^2}{12\kappa}-\left(z-\dfrac{1}{2}\right)\log 2.\label{eq:AqPm}
\end{align}

Let $q, \kappa\in (0,1)$ and $q=e^{-\kappa}$. Then for $z\in\C$ and $m\in\Z_{\geq 1}$

\begin{align}
\log (q^z; q)_\infty&=\AqP^+[\kappa, z]-\sum\limits_{n=1}^{m-1}\dfrac{B_{n+1}(z)B_{n}}{n(n+1)!}\kappa^{n}+\Error_m^+[\kappa, z],\label{eq:factorialspos}\\
\log (-q^z; q)_\infty&=\AqP^-[\kappa, z]-\sum\limits_{n=1}^{m-1}\left(2^{n}-1\right)\dfrac{B_{n+1}(z)B_{n}}{n(n+1)!}\kappa^{n}+\Error_m^-[\kappa, z].\label{eq:factorialsneg}
\end{align}

where $B_{k}(z)$ and $B_{k}$ denote Bernoulli polynomials and Bernoulli numbers (see Section <ref>).
Moreover, for any $\e\in (0,1/2)$ and any $b\in (m-1,m)$, there exists $C>0$ and $\kappa_0>0$  so that for all $\kappa\in(0,\kappa_0)$ and all $z\in\C$ with $|\textup{Im}(z)|<\tfrac{5}{\kappa}$,

\begin{equation}\label{eq:errorbounds}
\Big \lvert \Error^{\pm}_m[\kappa, z]\Big \rvert \leq  C\left(\kappa(1+|z|)^2 + \kappa^b(1+|z|)^{1+2b+\e}\right)
% \kappa^{1+m-\varepsilon}\left(|z|^{2+m}+|z|^{m}\right).
\end{equation}

In this paper we will only make use of the above asymptotic result when $m=1$. In that case, the sum involving the Bernoulli polynomials and numbers is not present. (In that sum terms where $n$ is odd and greater than $1$ are all zero (since $B_n=0$ for odd $n>1$.)

In the course of the our proof of this proposition, we actually prove a more precise asymptotic which depends on whether $\textup{Re}(z)\geq 1/2$ or $\textup{Re}(z)< 1/2$ ($1/2$ could be replaced by any real strictly between $0$ and $1$). Namely, when $\textup{Re}(z)\geq 1/2$ the error bound in (<ref>) can be replaced by $ C\kappa^b(1+|z|)^{1+2b+\e}$. However, when $\textup{Re}(z)< 1/2$, we must make use of a functional identity involving the Jacobi theta functions, see Section <ref>. This introduces some additional terms in the asymptotics:

\Error^{\pm}_m[\kappa, z] = \sum\limits_{n=1}^{m-1}\dfrac{B_{n+1}(1)B_{n}}{n(n+1)!}\kappa^{n}+\dfrac{\kappa}{8}(1-2z)^2+ \widetilde{\Error}^+_m[\kappa, z]\quad \textrm{for}\quad
\Big\lvert\widetilde{\Error}^+_m[\kappa, z]\Big\rvert \leq  C \kappa^b(1+|z|)^{1+2b+\e}.

The function $(q^z;q)_{\infty}$ is closely related to the $q$-gamma function which is defined as

\Gamma_q(z) = (1-q)^{1-z} \frac{(q;q)_{\infty}}{(q^z;q)_{\infty}},

thus our result can be seen as an asymptotic result for the $q$-gamma function as well. Asymptotics of $\Gamma_q$ have been studied in a number of contexts previously, e.g. [41, 40, 18, 52]. In all of those works (and others) the error bounds are either for $z$ fixed as $\kappa$ goes to zero, or $\kappa$ fixed as $z$ goes to infinity in some direction. To our knowledge, there has been no analysis of how these two limits balance. This balance, however, is incredibly important for us since we will deal with measures that are defined with respect to these $q$-Pochhammer symbols and certain key asymptotics that we perform in Section <ref> will involve probing $|z|$ of order $\kappa^{-1}$, with $\kappa$ going to zero.

The proof of Proposition <ref> that we provide in Section <ref> relies heavily on complex analytic methods often used in analytic number theory [43, 47] such as the Melin transform and the use of gamma, zeta, Hurwitz zeta and Jacobi theta functions, as well as functional identities for them. The formula (<ref>) can already be found in <cit.>, though the error bound stated there involves fixed $z$ with $\kappa$ tending to zero. In fact, the proof of <cit.> relies on an incorrect result, <cit.>, which claims that $\zeta(s,z) = \mathcal{O}(|t|^{2N+1})$ as $|t|\to \infty$ uniformly for $z$ in any compact subset of $\textup{Re}(z)>0$. Here $\zeta(s,z)$ is the Hurwitz zeta function and $s=\sigma+\iu t$ with real $\sigma >-2N$ for $N\in \Z_{\geq 1}$ and $t\in\R$. Proposition <ref> provides a correct bound on the Hurwitz zeta function with exponential growth $e^{|\textup{Arg(z)}\cdot t|}$ and important polynomial factors which ultimately translate into our error bound above.

The analysis of $(-q^z;q)_{\infty}$ is similar, though it involves the Dirichlet eta function as well.

§ OPEN ASEP

Though the titular focus of this work is on the open Kardar-Parisi-Zhang (KPZ) equation, we will almost entirely focus on the open ASEP. Everything that we have to say about the open KPZ equation is derived by taking appropriate asymptotics of results pertaining to the height function for the open ASEP. For a reader who is new to the study of open ASEP, we have included in Section <ref> an explanation for the general physical motivation for its study. Section <ref> contains a precise definition of open ASEP. Section <ref> returns to the physical study and describes the phase diagram for the current in stationarity, and Section <ref> records what is known of the fluctuations of the stationary measure. A reader who purely wishes to follow the proof of Theorem <ref> need only read Section <ref>. The other discussion provides useful context to understand our result.

§.§ Physical motivation

Consider the exclusion process in Figure <ref>. Particles (red dots) occupy vertices of a graph and jump along edges according to exponential rates subject to the exclusion rule (jumps to occupied sites are suppressed). There are reservoirs (grey squares) which (according to exponential rates) insert particles into unoccupied neighboring sites, or remove particles from occupied neighboring sites. This models transport through a network with sources and sinks.

An exclusion process on a network with reservoirs. Dotted red lines show the possible jumps of the left-most particle and the dotted red line with an $\times$ through it, shows an excluded jump. Other possible jumps and jump rates are not shown.

The presence of reservoirs leads to remarkable physical behavior. Since the total number of particles is not conserved, even when in its steady state (i.e., started its unique invariant measure) there will typically be a net current of particles flowing through this system—like how a flowing stream may have a stationary density profile while still moving water from source to sink. In statistical physics such models are said to have a non-equilibrium steady state (see the review [7]). Because the Markovian dynamics which produce these steady states are not reversible, non-equilibrium steady states do not take the Boltzmann weight form common in equilibrium statistical mechanics.  This renders the study of thermodynamic limits of non-equilibrium steady states quite challenging.

The presence of reservoirs can also produce phase transitions as the number of nodes ($N$ in our case) in the network grows. We are unaware of any precise statements to this effect when dealing with general networks. We will now, and for the rest of this work, focus entirely on the one-dimensional exclusion process with reservoirs where such a phase transition is understood. We will simply refer to this one-dimensional model as the open ASEP. Besides serving as transport model, the height function for open ASEP is also connected with stochastic interface growth and, through exponentiation, the stochastic heat equation—both important interpretations in our study.

Open ASEP additional enjoys exactly solvable structure. In particular, its invariant measure can be characterized via the matrix product ansatz and in certain cases, by Askey-Wilson processes.

§.§ Defining open ASEP

In one-dimension, there are four well-studied variants of ASEP: the full line ASEP, the period ASEP, the half-line ASEP and the open ASEP. We will focus on the last of these models which involves the most parameters and, in particular, has two sources/sinks.

Fix six parameters $q\in [0,1)$, $\alpha, \beta>0$, $\gamma, \delta\geq 0$, and $N\in \Z_{\geq 1}$.  Open ASEP is a continuous-time Markov process on the state space $\{0,1\}^{\llbracket 1,N\rrbracket}$. The state at time $t$ is denoted by $\tau(t)=\big(\tau_1(t),\ldots, \tau_N(t)\big)$;  sites $x\in \{1,\ldots N\}$ where $\tau_x(t)=1$ are said to be occupied by a particle, and those where $\tau_x(t)=0$ are unoccupied. The process $\tau(t)$ is defined by its infinitesimal generator $\mathcal{L}$ which acts on functions $f(\tau)$ of the state space as

\begin{align*}
\big(\mathcal{L}f\big)(\tau) &= \Big(\alpha(1-\tau_1)+\gamma\tau_1\Big)\cdot \Big(f(\tau^1)-f(\tau)\Big)\\
&\quad+\sum_{x=1}^{N-1} \Big( \tau_x(1-\tau_{x+1}) + q \big(1-\tau_x)\tau_{x+1} \Big)\cdot \Big(f(\tau^{x,x+1})-f(\tau)\Big)\\
&\quad+\Big(\beta\tau_N+ \delta(1-\tau_N)\Big) \cdot\Big(f(\tau^N)-f(\tau)\Big).
\end{align*}

where $\tau^1$ means to replace $\tau_1$ by $1-\tau_1$, $\tau^{x,x+1}$ means to switch the values of $\tau_x$ and $\tau_{x+1}$, and $\tau^N$ means to replace $\tau_N$ by $1-\tau_N$.
In words, particles jump left or right from occupied sites to unoccupied sites within $\llbracket 1,N\rrbracket$ at rate $q$ or $p=1$, respective; at the left boundary, sites become occupied (if presently unoccupied) at site $1$ at rate $\alpha$ and become unoccupied (if presently occupied) at rate $\gamma$; at the right boundary, particles become occupied (if unoccupied) at site $N$  at rate $\delta$ and become unoccupied (if occupied) at rate $\beta$. All moves are from independent exponential clocks.

From the occupation variable process $\tau(t)$ we define the height function process $h_N(t,x)$. The subscript indicates the lattice size $N$, and the time $t$ and spatial location $x$ are both arguments. The dependence of $h_N(t,x)$ on the other parameters $q,\alpha,\beta,\gamma,\delta$ will be generally suppressed from the notation. There will be two exceptions to this notation. The first will occur in Section <ref> where we will write $h_{N}(t,x;p,q,\alpha,\beta,\gamma,\delta)$ for the height function corresponding to the specified parameters. The other exception will occur in  Section <ref> where we will write $h_{N;u,v}$ for the height function which is uniquely labeled by $u$ and $v$.
The height function is defined for $t\geq 0$ and $x\in \llbracket 0, N\rrbracket$ as

\begin{equation}\label{eq_hN}
h_N(t,x) := h_N(t,0) + \sum_{i=1}^{x} (2\tau_i(t)-1),\qquad \textrm{where}\quad h_N(t,0):=-2\mathcal{N}_N(t)
\end{equation}

where the net current $\mathcal{N}_N(t)$ equals the number of particles to enter into site $1$ from the left reservoir minus the number of particles to exit from site $1$ into the left reservoir, up to time $t$. The height function definition is extended to $x\in [0,N]$ by linear interpolation—see  Figure <ref> for an illustration.

The open ASEP has a unique invariant probability measure. In the physics literature this is often denote by $\pi_{N}(\tau)$, with the dependence on the other parameters $q, \alpha, \beta, \gamma, \delta$ implicit. We will denote the expectation of a function $f:\{0,1\}^{\llbracket 1,N\rrbracket}\to \R$ under the invariant measure by

\big\langle f\big\rangle_N:= \!\!\!\!\sum_{\tau\in \{0,1\}^{\llbracket 1,N\rrbracket}} f(\tau)\cdot \pi^{\mathrm{ASEP}}_N(\tau).

For $\tau$ distributed according to $\pi^{\mathrm{ASEP}}_N$ we write $h_N(x)$ for the associated random height function and call this the stationary measure for the open ASEP height function process.
We use the term invariant measure when referring to the open ASEP particle process $\tau(t)$ and stationary measure when referring to the open ASEP height function process. The open ASEP height function process does not have an invariant measure since $\mathcal{N}_N(t)$ has a drift. However, if we start $h_N(0,x) = h_N(x)$ then at a later time $t>0$, $h_N(t,x)-h_N(t,0)$ will still have the law of $h_N(x)$ as a process in $x$.

§.§ Phase diagram

The boundary reservoirs play a key role in determining the limiting behavior of open ASEP as $N$ grows to infinity.
Define the stationary current to be

J_N:= \frac{\big\langle \alpha(1-\tau_1)-\gamma\tau_1 \big\rangle_N}{1-q}.

If we start open ASEP from $\pi^{\mathrm{ASEP}}_N$, then it follows from stationarity that for all $t$

J_N= \frac{\big\langle \mathcal{N}_N(t)\big\rangle_N}{t(1-q)}.

Here we have slightly abused notation since we are using the $\langle\cdot\rangle_N$ to represent the expectation of open ASEP started from in its stationary state, not just the expectation of a function of the state space.
$J_N$ represents the average (normalized by the drift $1-q$) net current of particles entering the system. By the internal conservation of particles (they are only created/removed at the boundaries), $J_N$ also represents the average net current moving through any other bond.

The phase diagram for open ASEP describes how the large $N$ limit of $J_N$ depends on the boundary rates. The existence of a limit is non-trivial, let alone computing its value as a function of $q,\alpha,\beta,\gamma,\delta$. The following phase diagram was developed using the matrix product ansatz in [20] for $q=0$ and [50] for general $q$; see also related work on the phase diagram in [19, 25, 26, 49]. A version of this phase diagram arose in a related context in work of [38].

The limit $J:=\lim_{N\to \infty} J_N$ exists and is a function of two parameters $\rho_{\ell}$ and $\rho_r$ which have a physical interpretation as effective boundary densities at the left and right boundaries. We will first describe heuristically how $\rho_{\ell}$ and $\rho_r$ are determined from $q,\alpha,\beta,\gamma,\delta$, and then how they determine $J$. The proof of the phase diagram is more complicated and proceeds through asymptotics for exact formulas for $J_N$ afforded by the matrix product ansatz introduced in Section <ref>.

Consider $N$ large and focus on the invariant measure near the left-boundary near $1$. Let us postulate that locally, asymptotically in $N\to \infty$, the invariant measure becomes product Bernoulli with density $\rho_\ell$ for some $\rho_\ell\in (0,1)$. Provided this is the case it is possible to determine the value of $\rho_\ell$ through a simple consideration. Since particles are conserved within the bulk of open ASEP, we know that in stationarity, the net number of particles moving from the left reservoir into site $1$ must equal the net number moving from site $1$ to $2$. This implies a simple conservation equation

\alpha (1-\rho_\ell) -\gamma \rho_\ell = (1-q) \rho_\ell(1-\rho_\ell).

On the left, the first term accounts for the rate $\alpha$ at which particles enter from the boundary to site $1$ provided it is unoccupied (which happens in stationarity with probability $1-\rho_\ell$) and the second term accounts for the rate $\gamma$ at which particles exit to the boundary from site $1$ provided it is occupied (which happens in stationarity with probability $\rho_\ell$). On the right, the factor $\rho_\ell(1-\rho_\ell)$ is the probability (under the Bernoulli product measure assumption) of having a particle and a hole next to each other in any prescribed order. The factor $(1-q)$ comes from the rate of particles jumping right minus the rate jumping left. Solving this quadratic equation gives us $\rho_\ell = \rho_\ell(q,\alpha,\gamma)$. A similar consideration allows us to solve a quadratic equation for $\rho_r=\rho_r(q,\beta,\delta)$.  In both cases, there is a single positive solution which gives the effective density.

We will write the solutions to the above quadratic equations explicitly in the following form:

\begin{equation}\label{eq_rho}
\rho_\ell = \frac{1}{1+C},\qquad \rho_r = \frac{A}{1+A}
\end{equation}

where $A,B,C$ and $D$ are given in (<ref>). Since  $A, C>0$, the densities $\rho_\ell$ and $\rho_r$ are in $(0,1)$.

Open ASEP phase diagram for the current $J$.

The phase diagram for the value of $J$ is determined entirely by $\rho_\ell$ and $\rho_r$ (or equivalently $A$ and $C$) as follows (see also Figure <ref> for an illustration):

* Maximal current: Here $\rho_\ell>1/2$ and $1-\rho_r>1/2$, or equivalently $C<1$ and $A<1$. The left boundary creates particles at a fast enough rate (and does not  remove them too quickly) and the right boundary removes particles at a fast enough rate (and does not create them too quickly) so that the system is able to transport particles from left to right in the bulk at its maximal efficiency. The maximal rate of transport for particles is $1/4$ and is achieved when the bulk density is $1/2$. To see this, note that if we have local product measure of density $\rho$, then the rate at which particles move (normalized by $1-q$) will be $\rho(1-\rho)$. This is maximized at $\rho=1/2$ and takes value $1/4$ there. Thus, in this phase, $J=1/4$ and the density will look locally in the bulk like product Bernoulli with density $1/2$.

* Low density: Here $\rho_\ell+\rho_r<1$,  or equivalently $C<A$ and $A>1$. The left boundary creates particles relatively slowly and the right boundary removes them fast enough so that they do not build up there. As such, the density in the bulk is determined by the effective density on the left, $\rho_\ell$ and the current $J=\rho_\ell(1-\rho_\ell)$.

* High density: Here $\rho_\ell+\rho_r>1$, or equivalently $C>A$ and $C>1$. By applying a particle-hole transform this is equivalent to the low density phase for the holes. In particular, now the right boundary removes particles so slowly and there is enough input from the left boundary that the density builds up and becomes $\rho_r$. Thus, the current $J=\rho_r(1-\rho_r)$.

Though there is interesting behavior at the line $\rho_\ell+\rho_r=1$, we will not address it here. Another point of interest is the triple point when $\rho_\ell=\rho_r=1/2$. Our KPZ equation results will ultimately come from zooming into that point.

Besides the phases addressed above, there are two other ways of categorizing the phase diagram:

* Fan: Here $\rho_\ell>\rho_r$, or equivalently $AC<1$. Very close (in a scaling going to zero relative to $N$) to the left boundary, the density exceeds that of the bulk. This is called the fan region since going from high to low density in the Burgers equation produces a rarefaction fan.

* Shock: Here  $\rho_\ell<\rho_r$, or equivalently $AC>1$.  Very close (in a scaling going to zero relative to $N$) to the left boundary, the density is lower than that of the bulk. Thus, in the same spirit of the Burgers equation, one sees shock-type behavior here.

The boundary between the fan and shock regions, when $\rho_\ell=\rho_r$ or equivalently $AC=1$, is special since it is the only part of the phase diagram where the invariant measure is simple. Along this line, in the finite $N$ open ASEP, the invariant measure is Bernoulli product with density $\rho=\rho_\ell=\rho_r$.

§.§ Fluctuations

The above discussion has pertained to the law of large number (or hydrodynamic) behavior of the invariant measure for open ASEP as $N\to \infty$. It is natural to consider fluctuations around that.
We recall the following result which provides a limit theorem for the fluctuations of the invariant measure in the maximal current phase around its law of large number limit density. This result was proved in the special case of TASEP when $q=0$ by Derrida-Lebowitz-Enaud [21] and recently extended to all $q\in (0,1)$ by Bryc and Wang [11].

Recall $\hscale(X)$ from (<ref>), the diffusively rescaled stationary measure for the open ASEP height function. Assume that the parameters $\alpha,\beta,\gamma,\delta$ a fixed and correspond to the interior of the maximal current phase, then

\begin{equation}\label{eq_WB_fluct}
\hscale(X) \Longrightarrow \frac{B(X)+B^{ex}(X)}{\sqrt{2}}.
\end{equation}

Here, the convergence occurs as $N\to \infty$ and is in finite dimensional distributions. The limit on the right-hand side is the sum of a standard Brownian motion $B(X)$ and an independent standard Brownian excursion $B^{ex}(X)$ (a Brownian bridge on the interval $[0,1]$, started at $X=0$ and ended at $X=1$ at zero, and conditioned to stay positive for all $X\in(0,1)$).

Limiting fluctuations of the invariant measure height function around the hydrodynamic limit, as proved in <cit.>. Here $B^{me}$ and $\hat{B}^{me}(X)$ denote versions of Brownian meanders.

Bryc and Wang <cit.> actually provides a characterization of the fluctuations the entire fan region, not just the maximal current region (which is included in the fan region). We do not reproduce all of their results here though Figure <ref> provides a graphical depiction of the various limits based on the region of the phase diagram. The open KPZ stationary measures that we construct in this current paper arise as the fluctuations of the open ASEP stationary measure when $q$ goes to $1$ and the boundary parameters are scaled so as to focus in on the triple point in the phase diagram where the maximal current and high / low density phases meet.

The above result is proved by performing asymptotic analysis on the Askey-Wilson process formulas for the open ASEP stationary measure generating function. This procedure produces formulas for the Laplace transform of the finite dimensional distributions of the limit. Identifying the Laplace transform formulas with the sum of a Brownian motion and excursion is achieved in [14] and amount a duality between the limit of the Askey-Wilson process and the Brownian processes. The proof of Theorem <ref> (<ref>) follows the same road map, though the asymptotic analysis is considerably more involved and we currently do not have a nice probabilistic description for the limiting process (which is only defined in terms of its Laplace transform).

§ WEAK ASYMMETRY LIMIT TO THE OPEN KPZ EQUATION

Section <ref> provides a precise definition of what it means to solve the open KPZ equation as introduced in [15]. Section <ref> formulates the key result of this section, Proposition <ref>, which shows how the open KPZ equation arises as a limit of the open ASEP under weak asymmetry scaling, near the triple point of the phase diagram. This result is proved in Section <ref> by appealing to the work of [15, 45] and suitably matching parameters.

§.§ Defining the open KPZ equation

We now define the Hopf-Cole solution to the open KPZ equation (<ref>) with inhomogeneous Neumann boundary parameters $u$ and $v$ (<ref>). There are two impediments to defining a notion of solution here, both of which stem from the fact that $H(T,X)$ should have regularity in $X$ like a Brownian motion. Thus, the non-linearity (the square of the derivative in $X$) will not be classically defined. Similarly, the boundary conditions, which involve the derivative in $X$, will not make sense. Instead of trying to make direct sense of this equation, we proceed through the Hopf-Cole transform and make sense of the resulting stochastic heat equation through its mild formulation, incorporating the boundary conditions into the heat kernel used in that formulation.  The open KPZ equation is then defined as the logarithm of the stochastic heat equation solution. This approach was first taken for the open KPZ equation in [15]. A direct approach to making sense of the open KPZ equation was developed via regularity structures in [30] which results in an equivalent notion of solution as we provide here.

We say that a stochastic process $H(T,X)$, with $T\geq 0$ and $X\in [0,1]$, solves the open KPZ equation with inhomogeneous Neumann boundary parameters $u$ and $v$ if $H(T,X) = \log Z(T,X)$ where
$Z$ is a mild solution to the stochastic heat equation  (<ref>) with inhomogeneous Robin boundary conditions (<ref>). We say that a stochastic process $Z(T,X)$, with $T\geq 0$ and $X\in [0,1]$, is a mild solution to the stochastic heat equation with boundary parameters $u$ and $v$ and initial data $Z_0$ if (1) $Z(T,\cdot)$ is adapted to the filtration generated by $Z_0$ and the space-time white noise $\xi$ up to time $T$, and if (2) $Z$ satisfies

\begin{equation}\label{eq:mildSHE}
Z(T,X)= \int_0^1 P_{u,v}(T,X,Y) Z_0(Y)dY + \int_0^1 \int_0^\infty P_{u,v}(T-S,X,Y)Z(S,Y) \xi(dS,dY)
\end{equation}

for all $T>0$ and $X\in [0,1]$. Here, the integral against $\xi$ is in the sense of Itô, and $P_{u,v}(T,X,Y)$ is the Gaussian heat kernel on $[0,1]$ with inhomogeneous Robin boundary conditions. Specifically, for $T\geq 0$ and $X,Y\in [0,1]$, $P_{u,v}(T,X,Y)$ satisfies

\partial_T P_{u,v}(T,X,Y) = \partial_X^2 P_{u,v}(T,X,Y)\qquad P_{u,v}(0,X,Y) = \delta_{X=Y}

with Dirac delta initial data
and boundary conditions for all $T>0$ and $Y\in [0,1]$,

\partial_X P_{u,v}(T,X,Y) \big\vert_{X=0} = (u-\tfrac{1}{2}) P_{u,v}(T,0,Y),\quad \textrm{and}\quad \partial_X P_{u,v}(T,X,Y) \big\vert_{X=1} = -(v-\tfrac{1}{2}) P_{u,v}(T,1,Y).

The following proves the existence, uniqueness and strict positivity for the solution of the SHE. It is <cit.> for $u, v\geq \tfrac{1}{2}$) and <cit.> for general $u,v\in \R$.

For a function-valued initial data $Z_0$ with $\sup\big\{\E\big[Z_0(X)^2\big]:X\in[0,1]\big\}<\infty$ there exists a mild solution to the SHE (see (<ref>)) which is unique in the class of adapted processes $Z$ which satisfy $\sup\big\{\E\big[Z(T,X)^2\big]:X\in[0,1],T\in[0,T_0]\big\}<\infty$ for all $T_0>0$. Furthermore, if $\P\big(Z_0(X)>0\, \forall X\in[0,1]\big)=1$ (i.e. $Z_0$ is almost surely strictly positive) then $\P\big(Z(T,X)>0 \,\forall X\in[0,1], T\geq 0\big)=1$ (i.e., the mild solution $Z$ is also almost surely strictly positive).

§.§ Formulating the weak asymmetry limit result

The open ASEP height function process converges to the Hopf-Cole solution to the open KPZ equation under the following assumptions.

* Weak asymmetry scaling:

\begin{equation}\label{eq_qN}
q= \exp\left(-\frac{2}{\sqrt{N}}\right).
\end{equation}

* Liggett's condition:

\begin{equation}\label{eq_Liggett}
\alpha +\gamma/q=1\qquad \textrm{and}\qquad \beta+\delta/q=1.
\end{equation}

This condition first arose in [37, 38]. Substituting (<ref>) into (<ref>) we find that
$B=D=-q$. On the other hand, $A$ and $C$, or equivalently $\rho_\ell$ and $\rho_r$, remain free parameters. This means that even when enforcing Liggett's condition we are able to access the entire phase diagram.

* Triple point scaling: The boundary densities $\rho_\ell$ and $\rho_r$ are scaled into an $N^{-1/2}$ window of the triple point. Specifically, we assume that for some $u, v\in \R$, as $N\to \infty$

\begin{equation}\label{eq_triple_point_scaling}
\rho_\ell = \frac{1}{2} +\frac{u}{2} N^{-1/2} + o(N^{-1/2}),\qquad \textrm{and}\qquad \rho_r = \frac{1}{2} -\frac{v}{2} N^{-1/2} + o(N^{-1/2}).
\end{equation}

The term $o(N^{-1/2})$ is any quantity which, when divided by $N^{-1/2}$, goes to zero as $N\to \infty$.

* $4:2:1$ height function scaling and Hopf-Cole transform: For $T\geq 0$ and $X\in [0,1]$ define

\begin{align}\label{eq_hNtx_scaled}
\hscale(T,X) &:= N^{-1/2} h_N(\tfrac{1}{2} e^{N^{-1/2}} N^2 T,N X) + (\tfrac{1}{2} N^{-1} +\tfrac{1}{24}) T,\\
\zscale(T,X) &:= e^{\hscale(T,X)},
\end{align}

where $h_N$ is the ASEP height function process.

* Hölder bounds on the initial data: We have an $N$-indexed sequence of open ASEP initial data $h_N(x)$ which satisfies the following Hölder bounds which are stated in terms of $\zscale(X):= e^{\hscale(X)}$ where $\hscale(X) := N^{-1/2} h_N(NX)$. For all $\theta\in (0,\tfrac{1}{2}\big)$  and every $n\in \Z_{\geq 1}$, there exist finite constants $C(n), C(\theta,n)>0$ such that for every $X,X'\in[0,1]$ and every $N\in \Z_{\geq 1}$

\begin{equation}\label{eq_Holder}
\|\zscale(X)\|_n\leq C(n),\qquad \textrm{and}\qquad \|\zscale(X)-\zscale(X')\|_n\leq C(\theta,n)|X-X'|^{\theta}.
\end{equation}

Here we have employed the notation $\|\cdot\|_n:= \E\big[|\cdot|^n\big]^{1/n}$ for $n\in \Z_{\geq 1}$ where $\E$ represents the expectation over the possibly random ASEP height function initial data $h_N(x)$.

The following proposition is proved in Section <ref>.

Consider any $N$-indexed sequence of open ASEP processes with parameters and initial datum satisfying all of the conditions of Assumption <ref>. Then the law of $\big\{\zscale(\cdot,\cdot)\big\}_{N}$ is tight in the Skorokhod space $D\big([0,T_0],C([0,1])\big)$ for any fixed time horizon $T_0>0$ and any limit point is in $C\big([0,T_0],C([0,1])\big)$.
If we further assume that there exists a (possibly random) non-negative-valued function $Z\in C([0,1])$ such that, as $N\to \infty$, $\zscale(X)\Longrightarrow Z(X)$ in the space of continuous processes of $X\in [0,1]$), then, as $N\to \infty$,
\zscale(T,X) \Longrightarrow Z(T,X)
as a space-time process in $D\big([0,T_0],C([0,1])\big)$ for any fixed time horizon $T_0>0$.  Here the limit process $Z(T,X)$ is the unique mild solution to the stochastic heat equation in Definition <ref> with boundary parameters $u$ and $v$, and initial data $Z(X)$.

§.§ Proof of Proposition <ref>

This result is essentially contained in [15] for $u,v\geq 1/2$ and [45] for general $u,v\in \R$. The tightness is from <cit.> and <cit.> while the convergence result is from
<cit.> and <cit.>. The only difference from those works is that we have used a different parametrization. For the tightness, the boundary parameters play no role and hence our result follows immediately from that of [15, 45]. For the convergence we explain below how to relate our parametrization to that used in [15, 45].
We relate $\AA$ and $\BB$ from (<ref>) to our $u$ and $v$ via (<ref>) (note, we are using san-serif $\AA$ and $\BB$ here to distinguish these variables from those in equation (<ref>)). The results of <cit.> and <cit.> are stated for $\AA$ and $\BB$ fixed. However, as indicated by (<ref>), for our matching $\AA=\AA_N$ and $\BB=\BB_N$ depend on $N$ and tend to some fixed limits (namely $u-\frac{1}{2}$ and $v-\frac{1}{2}$) as $N\to \infty$. The proofs of <cit.> and <cit.> readily generalize to allow for this diminishing $N$ dependence (we do not reproduce them here).

In order to explain the matching with [15, 45], we will introduce a second version of open ASEP in which the jump rates and the height function will be labeled with a tilde. We will describe how this process converges to the open KPZ equation, and then explain how to relate the desired $h_N$ process convergence to the established $\tilde{h}_N$ convergence, thus proving the proposition.

Let $\tilde{h}_N(t,x)$ denote the height function for $N$-site open ASEP with left and right jump rates $\tilde{q}$ and $\tilde{p}$, and boundary rates $\tilde{\alpha}, \tilde{\beta}, \tilde{\gamma}, \tilde{\delta}$. The bulk jump rates are specified in terms of $N$ as

\begin{equation}\label{eq:qptildescale}
\tilde{q} := \frac{1}{2}\exp\left(-\frac{1}{\sqrt{N}}\right),\qquad\textrm{and} \qquad\tilde{p} := \frac{1}{2}\exp\left(\frac{1}{\sqrt{N}}\right).
\end{equation}

We assume that the boundary rates satisfy Liggett's condition—in [15, 45] notation:

\frac{\tilde\alpha}{\tilde p} + \frac{\tilde\gamma}{\tilde q}=1,\qquad\textrm{and}\qquad\frac{\tilde\beta}{\tilde p} + \frac{\tilde\delta}{\tilde q}=1.

Given these relations, the boundary rates can be parameterized by $\tilde{\mu}_{\AA}, \tilde{\mu}_{\BB}\in \big[\sqrt{\tilde{q}/\tilde{p}},\sqrt{\tilde{p}/\tilde{q}}\big]$ as

\begin{align}\label{eq:tildealpha}
&\tilde\alpha := \frac{\tilde{p}^{3/2}(\tilde{p}^{1/2} - \tilde{\mu}_{\AA} \tilde{q}^{1/2})}{\tilde p-\tilde q},
& \tilde\beta := \frac{\tilde{p}^{3/2}(\tilde{p}^{1/2} - \tilde{\mu}_{\BB} \tilde{q}^{1/2})}{\tilde p-\tilde q},\\
&\tilde\gamma := \frac{\tilde{q}^{3/2}(\tilde{q}^{1/2} - \tilde{\mu}_{\AA} \tilde{p}^{1/2})}{\tilde q-\tilde p},
&\tilde\delta := \frac{\tilde{q}^{3/2}(\tilde{q}^{1/2} - \tilde{\mu}_{\BB} \tilde{p}^{1/2})}{\tilde q-\tilde p}.
\end{align}

We scale the parameters $\tilde{\mu}_{\AA}$ and $\tilde{\mu}_{\BB}$ with $N$ as

\begin{equation}\label{eq:tildemu}
\tilde{\mu}_{\AA} = 1- \AA/N,\qquad \textrm{and}\qquad  \tilde{\mu}_{\BB} = 1- \BB/N
\end{equation}

for some $\AA,\BB\in \R$ fixed. With these scalings, [15, 45] then defines (we are using slightly different notation here than used there)

\begin{equation}\label{eq:tildezn}
\tilde{Z}_N(t,x) :=\exp\Big(-\tilde{\lambda} \tilde{h}_N(t,x) + \tilde{\nu}t\Big)
\end{equation}

where $t\geq 0$, $x\in [0,N]$ and

\begin{equation}\label{eq_lambda_nu}
\tilde{\lambda} := \frac{1}{2} \log\frac{\tilde{q}}{\tilde{p}} = -N^{-1/2},\qquad
\tilde{\nu} := \tilde{p}+\tilde{q} - 2\sqrt{\tilde{p}\tilde{q}}= \frac{1}{2}N^{-1} + \frac{1}{24} N^{-2} + o(N^{-2}).
\end{equation}

The scaled version of $\tilde{Z}_N(t,x)$ is given by

\tilde{Z}^{(N)}(T,X) = \tilde{Z}_N(N^{2}T,N X)

where $T\geq 0$ and $X\in [0,1]$.

The following is shown as <cit.> ($\AA,\BB\geq 0$) and <cit.> ($\AA,\BB,\in \R$).

Consider any $N$-indexed sequence of open ASEP initial datum $\tilde{h}_N(x)$ that satisfies the following Hölder bounds, formulated in terms of $\tzscale(X):= \exp\big(N^{-1/2}\tilde{h}_N(N^2T,NX)\big)$: For all $\theta\in (0,\tfrac{1}{2}\big)$  and every $n\in \Z_{\geq 1}$, there exists a finite constant $C(\theta,n)>0$ such that for every $X,X'\in[0,1]$ and every $N\in \Z_{\geq 1}$

\begin{equation}\label{eq_Holderagain}
\|\tzscale(X)\|_n\leq C(\theta,n),\quad \textrm{and}\quad \|\tzscale(X)-\tzscale(X')\|_n\leq C(\theta,n)|X-X'|^{\theta}.
\end{equation}

Further, assume that there exists a (possibly random) non-negative-valued function $\tilde{Z}\in C([0,1])$ such that $\tzscale(X)\longrightarrow \tilde{Z}(X)$ in the space of continuous process of $X\in [0,1]$. Then, as $N\to \infty$

\tzscale(T,X)\longrightarrow\tilde{Z}(T,X)

as a space-time process in the Skorokhod space $D([0,T_0],C([0,1]))$ for any fixed time horizon $T_0>0$. Here the limit $\tilde{Z}(T,X)$ is the unique mild solution to the open stochastic heat equation in Definition <ref> with boundary parameters $\AA+\tfrac{1}{2}$ and $\BB+\tfrac{1}{2}$, and initial data $\tilde{Z}(X)$.

In order to apply this result and Proposition <ref> we need to match the parametrization of open ASEP in Section <ref> with that given above. To ease with this, let us augment the notation for the height function $h_N$ and $\tilde{h}_N$ so as to contain the rate parameters, i.e. we write $h_{N}(t,x;p,q,\alpha,\beta,\gamma,\delta)$.

We claim that if we specify

\begin{equation}\label{eq:notildetildequality}
q = \tilde{q} \tilde{p}^{-1},\quad \alpha = \tilde{\alpha} \tilde{p}^{-1},\quad   \beta = \tilde{\beta} \tilde{p}^{-1},\quad   \gamma = \tilde{\gamma} \tilde{p}^{-1},\quad   \delta = \tilde{\delta} \tilde{p}^{-1}.
\end{equation}

then there exists a coupling where, as processes in $t\geq 0$ and $x\in [0,N]$

\begin{equation}\label{eq:hnmatch}
h_N(t,x;1,q,\alpha,\beta,\gamma,\delta) = \tilde{h}_{\tilde{p}^{-1} t}(x;\tilde{p},\tilde{q},\tilde{\alpha},\tilde{\beta},\tilde{\gamma},\tilde{\delta})
\end{equation}

provided that the equality holds for $t=0$ (i.e., the initial data). This matching follows by coupling the driving Poisson point processes used in the graphical construction of open ASEP and noting that rescaling the time of a Poisson point process by a constant results (in law) in another Poisson point process with rate rescaled by that constant.

With our scaling of $q$ in (<ref>) of Assumption <ref>, and the scaling of  $\tilde{q}$ and $\tilde{p}$ from (<ref>), we already have the equality $q = \tilde{q} \tilde{p}^{-1}$. Equation (<ref>) of Assumption <ref> parameterizes $\alpha,\beta,\gamma,\delta$ in terms of $u$ and $v$ in such a way that

\begin{align}
&\alpha=\frac{1}{2}+\frac{u}{2} N^{-1/2} + O(N^{-1}), &\beta=\frac{1}{2}+\frac{v}{2} N^{-1/2} + o(N^{-1/2}),\\
&\gamma=\frac{1}{2}-\frac{u}{2} N^{-1/2} + O(N^{-1}), &\delta=\frac{1}{2}-\frac{v}{2} N^{-1/2} + o(N^{-1/2}).
\end{align}

Similarly the parametrization in (<ref>) and (<ref>) of $\tilde{\alpha}, \tilde{\beta}, \tilde{\gamma}$ and $\tilde{\delta}$ in terms of $\AA$ and $\BB$ implies that

\begin{align}
&\tilde\alpha=\frac{1}{4}+\Big(\frac{3}{8}+\frac{1}{4}\AA\Big) N^{-1/2} + O(N^{-1}), &\tilde\beta=\frac{1}{4}+\Big(\frac{3}{8}+\frac{1}{4}\BB\Big) N^{-1/2} + O(N^{-1}),\\
&\tilde\gamma=\frac{1}{4}-\Big(\frac{3}{8}+\frac{1}{4}\AA\Big) N^{-1/2} + O(N^{-1}), &\tilde\delta=\frac{1}{4}-\Big(\frac{3}{8}+\frac{1}{4}\BB\Big) N^{-1/2} + O(N^{-1}).
\end{align}

Thus, using the expansion $\tilde{p}^{-1} = 2+2N^{-1/2}+O(N^{-1})$, we find that in order for (<ref>) to hold, we must take $\AA=\AA_N$ and $\BB=\BB_N$ dependent on $N$ in such a way that

\begin{equation}\label{eq:ANBN}
\AA_N = u- \tfrac{1}{2} +o(1),\qquad \BB_N = v- \tfrac{1}{2} +o(1).
\end{equation}

Proposition <ref> is stated for $\AA$ and $\BB$ fixed. However, the proof of <cit.> and <cit.> permits a straightforward generalization (the details of which we forego repeating here) to allow $\AA_N$ and $\BB_N$ to tend to limiting values $A$ and $B$ as $N\to \infty$ (as in (<ref>)), with no change at all to the result. With this, we may now complete the proof of Proposition <ref>.

In light of the matching (<ref>) and the expansion in (<ref>) of $\tilde{\nu}$, provided that (<ref>) holds, we may now identify (as processes in $T$ and $X$)

\begin{equation}\label{eq:zmatching}
\tzscale(T,X) = \zscale(T,X)\big(1+o(1)\big)
\end{equation}

where the $o(1)$ term is deterministic, coming from the expansion of $\tilde\nu$, and goes to $0$ as $N\to\infty$. By the hypotheses in Proposition <ref>, the initial data $\zscale(X)$ satisfies Hölder bounds and converges (as a spatial process) to some $Z(X)$. Thus, the same holds for $\tzscale(X)$, which means that we can apply Proposition <ref> with $\AA=\AA_N$ and $\BB=\BB_N$ tending, as in (<ref>), to respective limits $\AA=u-\tfrac{1}{2}$ and $\BB=v-\tfrac{1}{2}$ as $N\to \infty$. This, along with the relation (<ref>) implies that $\zscale(T,X)$ converges to $Z(T,X)$, the mild solution to the open stochastic heat equation in Definition <ref> with boundary parameters $\AA+\tfrac{1}{2}=u$ and $\BB+\tfrac{1}{2}=v$, and initial data $Z(X)$.

§ ATTRACTIVE COUPLING OF OPEN ASEPS WITH DIFFERENT BOUNDARY PARAMETERS

The goal of this section is to prepare for the proof of Theorem <ref> (<ref>)-(<ref>) in Section <ref>. Section <ref> contains our main open ASEP coupling result which is proved by appealing to a multi-species open ASEP. Some implications (e.g. coupling and Hölder bounds) of this result in terms of the open ASEP height function and its Hopf-Cole transform are given in Section <ref>.

§.§ Coupling via multi-species open ASEP

We give an attractive coupling of open ASEPs with different boundary conditions. Attractive means that if the occupation variables start ordered between different ASEPs, then they will remain ordered. As is standard in proving attractive couplings (e.g. [3]), we appeal to a multi-species version of the model.

Fix $q\geq 0$, any $M\geq 2$ and any collection of $M$-tuples of non-negative real numbers $\{\alpha^i\big\}_{i=1}^{M}$,  $\{\beta^i\big\}_{i=1}^{M}$,  $\{\gamma^i\big\}_{i=1}^{M}$ and  $\{\delta^i\big\}_{i=1}^{M}$ such that for all $i<j$,

\alpha^i\leq \alpha^j,\qquad \beta^i\geq \beta^j,\qquad \gamma^i\geq \gamma^j,\qquad \delta^i\leq \delta^j.

For each $i\in \{1,\ldots, M\}$ fix any initial data $\tau^i = (\tau^i_x)_{x\in \llbracket 1,N\rrbracket}\in \{0,1\}^{N}$ such that for all $1\leq i<j\leq M$,  $\tau^i\leq \tau^j$ (i.e., $\tau^i_x\leq \tau^j_x$ for all $x\in \llbracket 1, N\rrbracket$). Let $\tau^i(\cdot)$ denote the $N$ site open ASEP with parameter $(q,\alpha^i,\beta^i,\gamma^i,\delta^i)$ started with $\tau^i(0)=\tau^i$. Then:

* There exists an attractive coupling of all $M$ open ASEP processes $\tau^1(\cdot),\ldots, \tau^M(\cdot)$. In other words, there exists a common probability space which supports the $M$ processes and has the property that for all $t\geq 0$ and  $1\leq i<j\leq M$,  $\tau^i(t)\leq \tau^j(t)$.

* Let $\tilde{\tau}^i$ denote an occupation vector whose law is that of the invariant measure for the $N$ site open ASEP with parameters $(q,\alpha^i,\beta^i,\gamma^i,\delta^i)$. Then there exists a coupling of all $M$ invariant measures such that for all  $1\leq i<j\leq M$,  $\tilde\tau^i\leq \tilde\tau^j$.

The second claim in the lemma—the coupling of the invariant measures—follows easily from the first. Start with any initial datum such that $\tau^i(0)\leq \tau^j(0)$ for all $1\leq i<j\leq M$. By the first result of the lemma, the evolution can be coupled so as to maintain those inequalities for all $t$. However, for each $i$, the law of $\tau^i(t)$ converges to the law of the invariant measure $\tilde\tau^i$ as $t$ goes to infinity. This constructs the desired coupling of $\big(\tilde\tau^i\big)_{i=1}^{M}$  as $t$ goes to infinity limits of

The rest of this proof constructs the attractive coupling claimed in the first part of the lemma. The key idea is to introduce a suitable multi-class particles. We will introduce an $M+1$ species version of open ASEP. Sites $x\in \llbracket 1,N\rrbracket$ can be occupied by a single particle of type $1$ through $M+1$. The higher the type the lower the priority. In particular, this process has the following transition rates. Below we use the conventions that $\alpha^0=\beta^{M+1}=\gamma^{M+1}=\delta^0=0$.

* In the bulk (for $x\in\llbracket 1,N-1\rrbracket$) if the configuration at sites $x$ and $x+1$ is given by $AB$ (i.e., $A$ at site $x$ and $B$ at site $x+1$), then this becomes $BA$ with rate $1$ if $A<B$ and rate $q$ if $A>B$.

* At the left boundary, if the configuration at site $1$ is $A$, then this becomes $B$ at rate $\alpha^B-\alpha^{B-1}$ if $A>B$ and at rate $\gamma^B-\gamma^{B+1}$ of $A<B$.

* At the right boundary, if the configuration at site $N$ is $A$, then this becomes $B$ at rate $\beta^B-\beta^{B+1}$ if $A<B$ and at rate $\delta^B-\delta^{B-1}$ if $A>B$.

Denote the occupation variables for this process by $\eta^i_x(t)$ where $i\in \llbracket 1, M+1\rrbracket$, $x\in \llbracket 1,N\rrbracket$ and $t\geq 0$. In words, $\eta^i_x(t)$ is equal to $1$ if there is a type $i$ particle at position $x$ at time $t$, and equal to $0$ otherwise. From these multi-species occupation variables we define

\begin{equation}\label{eq_tautx}
\tau^i_x(t) = \sum_{j=1}^{i} \eta^j_x(t).
\end{equation}

From the $\tau^i_x$ in the statement of the lemma, we define initial data for the multi-species ASEP by

\begin{equation}\label{eq_etatau}
\eta^j_x(0) = \tau^i_x-\tau^{i-1}_x.
\end{equation}

Clearly, combining (<ref>)  with (<ref>) implies that $\tau^i_x(0)=\tau^i_x$. We claim that marginally, for each $i\in \llbracket 1,M\rrbracket$, $\tau^i(t)$ evolves as a process in $t$ precisely as open ASEP with parameters $(q,\alpha^i,\beta^i,\gamma^i,\delta^i)$. Provided this is true, the desired attractive coupling follows since for any $1\leq i<j\leq N$, the difference $\tau^j_x(t)-\tau^i_x(t) = \sum_{k=i+1}^{j} \eta^k_x(t)$ is positive, hence $\tau^i_x(t)\leq \tau^j_x(t)$ as desired.

The claim that $\tau^i(t)$ is marginally distributed as open ASEP with parameters $(q,\alpha^i,\beta^i,\gamma^i,\delta^i)$ is easily verified. In the bulk this is the standard multi-species coupling and well known. Let us consider the left boundary (the right one is similar). If $\tau^i_1(t)=1$ there is must be a multi-species particle $A$ at site $1$ with $A\leq i$. The rate at which $\tau^i_1(t)$ changes from $1$ to $0$ is then equal to the rate at which $A$ changes to any $B>i$. Since these are each exponential times, the rate of the minimum of them is equal to the sum of the rates. That is to say, the rate of $A$ changing to $B>i$ is equal to $\sum_{j=i}^{M}(\gamma^{j}-\gamma^{j+1})$ which (using the assumed convention $\gamma^{M+1}=0$) telescopes to equal $\gamma^i$. This is precisely the desired rate at which $\tau^i_1(t)$ changes from $1$ to $0$. Similarly, if $\tau^i_1(t)=0$ this means that there is a multi-species particle $A$ at site $1$ and $A> i$. The rate at which $\tau^i_1(t)$ changes from $0$ to $1$ is equal to the rate at which $A$ changes to any $B\leq i$. By the same argument as before, the rate of this change is equal to $\sum_{j=1}^{i}(\alpha^{j}-\alpha^{j-1})$ which telescopes to equal $\alpha^i$, the desired rate at which $\tau^i_1(t)$ changes from $0$ to $1$. This proves claim and completes the proof of the lemma.

§.§ Height function coupling and its implications

Armed with the attractive coupling of Lemma  <ref> we may now prove the following results.

Fix $u,v\in \R$ and consider the invariant measure for $N$ site open ASEP parameterized as in (<ref>) of Assumption <ref> by $u$ and $v$. Denote the associated height function by $h_{N;u,v}(x)$ for $x\in [0,N]$, centered so that $h_{N;u,v}(0)=0$. As in (<ref>), define $\hscale_{u,v}(X):=N^{-1/2} h_{N;u,v}(NX)$ for $X\in [0,1]$, and $\zscale_{u,v}(X) := \exp(\hscale_{u,v}(X))$.
Then the following holds:

*  For any $u\in \R$, $h_{N;u,-u}(x)$ has the law of a discrete time simple random walk with $x\in [0,N]$ and iid jumps in which the height increases by 1 with probability $\rho(u) = (1+q^u)^{-1}$ and decreases by 1 with probability $1-\rho(u)$. Consequently, $\hscale_{u,-u}(X)\Longrightarrow B_u(X)$ as a process in $X$, where $B_u$ is a drift $u$ Brownian motion with $X\in [0,1]$.

*  For any $u,v\in \R$, $h_{N;u,v}(x)$ has the same law as a process in $x\in [0,N]$ as $h_{N;v,u}(N-x)-h_{N;v,u}(N)$.

*  For all $u$ and $v$ such that $u+v\geq 0$ there exists a coupling of $\hscale_{u,v}$  with $\hscale_{u,-u}$ and $\hscale_{-v,v}$ such that for all $X,X'\in [0,1]$ with $X\leq X'$,

\hscale_{-v,v}(X')-\hscale_{-v,v}(X)\leq \hscale_{u,v}(X')-\hscale_{u,v}(X)\leq \hscale_{u,-u}(X')-\hscale_{u,-u}(X).

For all $u$ and $v$ such that $u+v\leq 0$ there exists a coupling of $\hscale_{u,v}$  with $\hscale_{u,-u}$ and $\hscale_{-v,v}$ such that for all $X,X'\in [0,1]$ with $X\leq X'$,

\hscale_{-v,v}(X')-\hscale_{-v,v}(X)\geq \hscale_{u,v}(X')-\hscale_{u,v}(X)\geq \hscale_{u,-u}(X')-\hscale_{u,-u}(X).

*  For all $u,v\in \R$ we have the following Hölder bound. For all $\theta\in (0,\tfrac{1}{2})$ and every $n\in \Z_{\geq 1}$, there exists a constant $C(\theta,n,u,v)>0$ such that for every $X,X'\in [0,1]$ and every $N\in \Z_{\geq 1}$,

\left\| \hscale_{u,v}(X)-\hscale_{u,v}(X')\right\|_n \leq C(\theta,n,u,v)\left| X-X'\right|^{\theta}.

*   For all $u,v\in \R$ we have the following Hölder bounds. For all $n\in \Z$ there exists $C(n,u,v)>0$ such that for all $N\in\Z_{\geq 1}$ and all  $X\in [0,1]$

\begin{equation}\label{eq:LNbound}
\| \zscale_{u,v}(X)\|_n \leq C(n,u,v),
\end{equation}

and for all $\theta\in (0,\tfrac{1}{2})$ and every $n\in \Z_{\geq 1}$, there exists a constant $C(\theta,n,u,v)>0$ such that for every $X,X'\in [0,1]$ and every $N\in \Z_{\geq 1}$,

\begin{equation}\label{eq:LNbound2}
\left\| \zscale_{u,v}(X)-\zscale_{u,v}(X')\right\|_n \leq C(\theta,n,u,v)\left|X-X'\right|^{\theta}.
\end{equation}

There are five parts to this proof. We prove each one sequentially.

Part (<ref>): When $v=-u$, the product $AC=1$ (recall $A$ and $C$ from (<ref>)) or equivalently, $\rho_\ell=\rho_r$ (recall $\rho_\ell$ and $\rho_r$ from (<ref>)). As explained in <cit.> or [27] (see also [33]), in this case the open ASEP has invariant measure which is Bernoulli product form with particle density $\rho=(1+A)^{-1}$. This implies that random walk description for $h_{N;u,-u}(x)$. The Brownian motion limit follows from the convergence of simple random walks to Brownian motion. Notice that $\rho(u) = \tfrac{1}{2}\big(1-u N^{-1/2} + O(N^{-1})\big)$, thus the mean of the jump distribution for $h_{N;u,-u}(x)$ is $uN^{-1/2} + O(N^{-1})$ and the variance is $1+O(N^{-1/2})$. This implies that the diffusively scale $\hscale_{u,-u}(X)$ converges to the desired Brownian limit.

Part (<ref>): This immediately follows from the particle/hole duality for open ASEP. Namely, for $i\in \llbracket 1,N\rrbracket$ let $\eta_{i}(t) = 1-\tau_{N-i+1}(t)$ denote the holes of open ASEP, after a space reversal. It is evident that if $\tau(t)$ evolves according to open ASEP with parameters $q,\alpha,\beta,\gamma,\delta$, then $\e(t)$ evolves according to open ASEP with parameters $q,\beta,\alpha,\delta,\gamma$. In terms of our $u,v$ parametrization, this amounts to interchanging $u$ and $v$. This implies that the invariant measures for open ASEP with parameters $u,v$ and $v,u$ are related by a particle/hole duality and space reversal. Translating this into the height function implies the claimed equality.

Part (<ref>): We combine (<ref>) with the parametrization in (<ref>) of $(A,B,C,D)$ in terms of $(\tilde u,\tilde v)$. (We are using $(\tilde u,\tilde v)$ here in place of $(u,v)$ since below we will want to assign $(\tilde u,\tilde v)$ to various values such as $(u,-u), (u,v)$ and $(-v,v)$). Using the above parametrization, we see that $(\alpha,\beta,\gamma,\delta)$ is expressed in terms of $(\tilde u,\tilde v)$ in such a way that when $\tilde u$ increases, $\alpha$ increases and $\gamma$ decrease, and when $\tilde v$ increases, $\beta$ increases and $\delta$ decreases (this agrees with the Taylor expansion in (<ref>).

Let us first consider when $u+v\geq 0$, or in other words when $u\geq -v$. We will use the coupling from  Lemma <ref> with $M=3$. Using the parametrization of $(\alpha,\beta,\gamma,\delta)$ in terms of $(\tilde u,\tilde v)$ from the previous paragraph, let us define three sets of parameters: $(\alpha^1,\beta^1,\gamma^1,\delta^1)$ is defined with respect to $(\tilde u,\tilde v) =(-v,v))$,  $(\alpha^2,\beta^2,\gamma^2,\delta^2)$ is defined with respect to $(\tilde u,\tilde v) = (u,v)$, and $(\alpha^3,\beta^3,\gamma^3,\delta^3)$ is defined with respect to $(\tilde  u,\tilde v) =(u,-u)$. From the properties of how the parametrization varies with changing $\tilde u$ and $\tilde v$, it follows that

\begin{equation}\label{eq_ordering}
\alpha^1< \alpha^2=\alpha^3,\quad \beta^1=\beta^2> \beta^3, \quad \gamma^1> \gamma^2=\gamma^3, \quad\delta^1=\delta^2<\delta^3.
\end{equation}

Thus, the conditions of Lemma <ref> holds with $M=3$. The lemma implies that the invariant measure of the three open ASEPs are ordered so that $\tilde{\tau}^1\leq \tilde{\tau}^2\leq \tilde{\tau}^3$. Since the associated height functions $h_{N;-v,v}$, $h_{N;u,v}$ and $h_{N;u,-u}$ are defined by taking sums of these respective occupation variables, the claimed coupling of the height function differences follows immediately (after diffusive scaling).

The case when $u+v\leq 0$ follows similarly. Since $u\leq -v$, the ordering of the $(\alpha^i,\beta^i,\gamma^i,\delta^i)$ is reversed from (<ref>) which is why the inequalities between the height function differences switch.

Part (<ref>): This makes use of results we have already demonstrated in parts (<ref>) and (<ref>) above. By part (<ref>), $\hscale_{-v,v}$ and $\hscale_{u,-u}$ are both diffusively rescaled random walks. By part (<ref>), the increments of $\hscale_{u,v}$ are bounded below and above in our coupled probability space by the increments  $\hscale_{-v,v}$ and $\hscale_{u,-u}$ (depending on the sign of $u+v$, the ordering is switched). Taking absolute values and the $n$-th power, this inequality yields

\begin{equation}\label{eq_hscale_hold}
\Big\vert \hscale_{u,v}(X')-\hscale_{u,v}(X)\Big\vert^n\leq \max\Big(\Big\vert \hscale_{-v,v}(X')-\hscale_{-v,v}(X)\Big\vert^n, \Big\vert \hscale_{u,-u}(X')-\hscale_{u,-u}(X)\Big\vert^n\Big).
\end{equation}

Given two non-negative numbers $a$ and $b$, note that $\max(a,b) \leq a+b\leq 2\max(a,b)$. Now, consider two non-negative random variables $A$ and $B$. Using the first inequality we find that  $\E\big[\max(A,B)\big] \leq \E[A]+\E[B]$ and using the second inequality we find that $\E[A]+\E[B]\leq 2\max\big(\E[A],\E[B]\big)$. Applying this to compute the expectation of the right-hand side of (<ref>) we find that

\begin{equation}
\left\| \hscale_{u,v}(X')-\hscale_{u,v}(X)\right\|^n_n\leq 2 \max\left(\left\| \hscale_{-v,v}(X')-\hscale_{-v,v}(X)\right\|^n_n, \left\| \hscale_{u,-u}(X')-\hscale_{u,-u}(X)\right\|^n_n\right).
\end{equation}

So, in order to prove our Hölder bound $\hscale_{u,v}$, it suffices to show it for both $\hscale_{-v,v}$ and $\hscale_{u,-u}$. These, however, are diffusively rescaled simple random walks which converge to drifted Brownian motions. Therefore, the Hölder bounds are quite standard (involving moments of Binomial distributions) and we do not reproduce them here.

Part (<ref>): We start by showing (<ref>) by, as in part (<ref>), appealing to the coupling from part (<ref>) and the random walk reduction from part (<ref>). The coupling in part (<ref>) with $X'=0$, implies that $\hscale_{u,v}(X)$ is bounded between $\hscale_{-v,v}(X)$ and $\hscale_{u,-u}(X)$ (with the ordering depending on the sign of $u+v$). Exponentiating we see that

\zscale_{u,v}(X) \leq \max\left(\zscale_{-v,v}(X),\zscale_{u,-u}(X)\right).

As argued in part (<ref>), this implies that

\begin{equation}\label{eq_zscalemax}
\left\|\zscale_{u,v}(X)\right\|_n^n \leq 2  \max\left(\left\|\zscale_{-v,v}(X)\right\|,\left\|\zscale_{u,-u}(X)\right\|\right).
\end{equation}

So, it is sufficient to show (<ref>)  with $\zscale_{u,v}(X)$ replaced by both $\zscale_{-v,v}(X)$ and $\zscale_{u,-u}(X)$. Since each of these expressions involves the exponential of a diffusively rescaled simple random walk, this is a standard calculation. As we did not include the analogous part (<ref>) calculation, we include this one.

Let $y_1,\ldots, y_N$ be i.i.d. random variables with probability $\P(y_1=1)=\rho:=(1+q^u)^{-1}$ and $\P(y_1=-1)=1-\rho$. Then we claim that for all $n\in \Z_{\geq 1}$ there exists $C(n,u)>0$ such that

\begin{equation}\label{eq:iidbound}
\E\Big[\Big(e^{N^{-1/2} (y_1+\cdots+y_M)}\Big)^n\Big] \leq C(n,u)
\end{equation}

for all $M\leq N\in\Z_{\geq 1}$.
It is clear that (<ref>) immediately follows by combining this with (<ref>) and the random walk description of both $\zscale_{-v,v}(X)$ and $\zscale_{u,-u}(X)$. So, we now focus on proving (<ref>).

By independence of the $y_i$, we can rewrite the left-hand side of (<ref>) as

\begin{equation}\label{eq:iidbound2}
\E\Big[\Big(e^{N^{-1/2} (y_1+\cdots+y_M)}\Big)^n\Big] = \E\Big[ \big(e^{N^{-1/2} y}\big)^n\Big]^{M} = \left(\rho e^{nN^{-1/2}}+(1-\rho)e^{-nN^{-1/2}}\right)^M
\end{equation}

where $y$ has the law of each of the $y_i$. We wish to bound the right-hand side above by $C(n,u)$ for all $M\leq N\in\Z_{\geq 1}$. Taylor expanding $\rho e^{nN^{-1/2}}+(1-\rho)e^{-nN^{-1/2}}=1+(n^2/2+un)N^{-1} +o(N^{-1})$. Raising this to the $M^{th}$ power yields a bound on the right-hand side of (<ref>) like $e^{(n^2/2+un)X}$ where $X=M/N\in [0,1]$. The maximum occurs at either $X=0$ or $X=1$ depending on the sign of $n^2/2+un$, and thus letting $C(n,u) = \max(1, e^{n^2/2+un})$ seems to work.

To make the above argument rigorous we just need to control the Taylor expansion error. If $ \rho e^{nN^{-1/2}}+(1-\rho)e^{-nN^{-1/2}}\leq 1$ then the right-hand side in (<ref>) is maximized when $M=0$, otherwise it is maximized when $M=N$. In the first case when $M=0$, the right-hand side of (<ref>) is obviously bounded by $1$. So, it suffices to bound the right-hand side of (<ref>) when $M=N$.

We recall a few basic inequalities which control the Taylor expansion and provide bounds which hold for all $n,N\in \Z_{\geq 1}$ and $u\in \R$. Below, $C$ will denote a positive constant which is independent of $n,N$ and $u$. Recalling that $q=e^{-2N^{-1/2}}$. First, we can bound $\rho=(1+q^u)^{-1}$ as

\left|\rho - (\tfrac{1}{2} - \tfrac{u}{2} N^{-1/2})\right| \leq C|u|^3N^{-3/2}.

With this we can show that

\left|\rho e^{nN^{-1/2}}+(1-\rho)e^{-nN^{-1/2}} - \left(\cosh(nN^{-1/2}) + \sinh(nN^{-1/2}) uN^{-1/2}\right) \right|\leq C\cosh(nN^{-1/2})|u|^3N^{-3/2}.

Now we can use the hyperbolic cosine and sine bounds

\left|\cosh(nN^{-1/2}) - \left(1 + \frac{n^2}{2} N^{-1}\right)\right| \leq C n^4 N^{-2},\qquad \left|\sinh(nN^{-1/2}) - nN^{-1/2} \right| \leq C n^3 N^{-3/2}.

Noting that $\cosh(nN^{-1/2})\leq \cosh(n)$ for $N\in \Z_{\geq 1}$, and writing $C(u,n)$ for a positive real number that now depends on $u$ and $n$ (and will vary between lines), we deduce from the above bounds that

\left|\rho e^{nN^{-1/2}}+(1-\rho)e^{-nN^{-1/2}} - \left(1 + \left(\frac{n^2}{2}+un\right) N^{-1}\right)\right|\leq C(u,n) N^{-3/2}.

From this it readily follows that

\left(\rho e^{nN^{-1/2}}+(1-\rho)e^{-nN^{-1/2}}\right)^N \leq C(u,n)

which is precisely the bound we sought to show in order to conclude (<ref>).

We now turn show (<ref>). This bound follows by combining the result of part (<ref>) with the already showed inequality (<ref>) in part (<ref>). First note that for all $a,b\in \R$,

\begin{equation}\label{eq_abineq}
\left| e^a-e^b\right| \leq \max\left(e^a,e^b\right) \left| a-b\right|.
\end{equation}

To see this, recall that for $c\geq 0$, $|1-e^{-c}|\leq c$. Thus, for $a<b$,
$| e^a-e^b| = e^a| 1-e^{b-a}|\leq e^a | a-b|$.
Similarly, for $b<a$,$\left| e^a-e^b\right| \leq e^b \left| a-b\right|$. Putting these together yields (<ref>).

Now, setting $a= \hscale_{u,v}(X)$ and $b=\hscale_{u,v}(X')$ in (<ref>) we find that

\left|\zscale_{u,v}(X)-\zscale_{u,v}(X')\right| \leq \max\left( \zscale_{u,v}(X),\zscale_{u,v}(X')\right) \left| \hscale_{u,v}(X)-\hscale_{u,v}(X')\right|.

Taking the expectation of the $n$-th power and applying Cauchy-Schwartz yields

\left\| \zscale_{u,v}(X)-\zscale_{u,v}(X')\right\|_n \leq \left\| \max\left(\zscale_{u,v}(X),\zscale_{u,v}(X')\right)\right\|_{2n} \cdot \left\|\hscale_{u,v}(X)-\hscale_{u,v}(X')\right\|_{2n}.

The result of part (<ref>) provides us with control over the second term on the right-hand side, yielding

\left\| \zscale_{u,v}(X)-\zscale_{u,v}(X')\right\|_n \leq \left\| \max\left(\zscale_{u,v}(X),\zscale_{u,v}(X')\right)\right\|_{2n} \cdot C(\theta,n,u,v)\left| X-X'\right|^{\theta}.

Thus it remains to prove that there exists $C'(n,u,v)>0$ such that for all $X,X'\in [0,1]$

\begin{equation}\label{eqmaxsza}
\left\| \max\left(\zscale_{u,v}(X),\zscale_{u,v}(x')\right)\right\|_{2n} \leq C'(n,u,v).
\end{equation}

We already know how to control the $n$-norms of the individual terms inside of the $\max$ in the left-hand side of (<ref>) by the already proved bound (<ref>). Thus, we need to extend this control to the $n$-norm of the $\max$.
For positive random variables $A$ and $B$, have that for any $m\in \Z_{\geq 1}$,

\begin{equation}\label{eq_ABmineq}
\E\Big[\max(A,B)^m\Big] \leq \E\Big[(A+B)^m\Big] \leq 2^m \max_{k\in \llbracket 0,m\rrbracket } \E\left[A^kB^{m-k}\right]\leq 2^m \max_{k\in \llbracket 0,m\rrbracket } \E\left[A^{2k}\right]^{1/2} \E\left[B^{2(m-k)}\right]^{1/2}.
\end{equation}

The first inequality uses $\max(A,B)\leq A+B$, the second bounds the $2^m$ terms in the binomial expansion by the maximal term times $2^m$ and the third inequality follows from Cauchy-Schwartz.

Now substitute $m=2n$, $A=\zscale_{u,v}(X)$ and $B=\zscale_{u,v}(X')$. By the bound already proved in (<ref>) it follows that $\E\left[A^{2k}\right]^{1/2}\leq C(2k,u,v)^k$ and likewise $\E\left[B^{2(m-k)}\right]^{1/2} \leq C(2(m-k),u,v)^{m-k}$. Thus,

\left\| \max\left(\zscale(x),\zscale(x')\right)\right\|_{2n} \leq 2 \max_{k\in \{0,\ldots,2n\}} C(2k,u,v)^{\frac{k}{2n}}C(2(2n-k),u,v)^{\frac{2n-k}{2n}}=:C'(n,u,v).

This is the desired bound (<ref>) and thus completes the proof of (<ref>).

§ PROOF OF THEOREM <REF> (<REF>)-(<REF>)

Proof of Theorem <ref>  (<ref>): The first claim is that the law of $\hscale_{u,v}(\cdot)$ is tight as an $N$-indexed sequence of measures in the space of continuous functions on $[0,1]$. This follows from Proposition <ref> (<ref>) and the fact that $\hscale_{u,v}(0)=0$ by applying the Kolmogorov continuity theorem. That theorem, in fact, further implies that all subsequential limits $H_{u,v}$ are almost surely Hölder $1/2-$.

The next claim is that any subsequential limits $H_{u,v}$ is a stationary measure for the open KPZ equation. Fix a subsequence $\{N_k\}_{k=1}^{\infty}$ along which $H^{(N_k)}_{u,v}$ converges to the limit point $H_{u,v}$. We claim that all of the conditions of Assumption <ref> are satisfied. The weak asymmetry, Liggitt's condition and the triple point scaling assumptions all follow from the choice of parameters we have made in Assumption <ref>. The only thing to verify is the Hölder bounds on the initial data, but that is precisely what is shown in Proposition <ref> (<ref>).
We may therefore apply Proposition <ref>. By our choice of the subsequence $\{N_k\}_{k=1}^{\infty}$, we have that  $Z^{(N_k)}_{u,v}(X)\Rightarrow Z_{u,v}(X)$ (as a process in $X$). Thus, the proposition implies that $Z^{(N_k)}_{u,v}(T,X)\Rightarrow Z_{u,v}(T,X)$ (as a process in $T$ and $X$), the unique mild solution to the stochastic heat equation with boundary parameters $u$ and $v$, and initial data $Z(X)$.

By the stationarity of the open ASEP initial data, the law of $Z^{(N_k)}_{u,v}(T,X)/Z^{(N_k)}_{u,v}(T,0)$ (as a process in $X$)  is $T$-independent. By the convergence to the stochastic heat equation, the same holds true for  $Z_{u,v}(T,X)/Z_{u,v}(T,0)$. Taking logarithms, this implies that the law of $H_{u,v}(T,X)-H_{u,v}(T,0)$ is likewise $T$-independent, which implies that the limit point $H_{u,v}(X)$ that we have constructed is a stationary measure for the open KPZ equation, as desired.

Proof of Theorem <ref>  (<ref>):  Proposition <ref> (<ref>) implies that the desired coupling holds along any subsequence $\{N_k\}_{k=1}^{\infty}$ such that all of the $\{H^{(N_k)}_{u_i,v_i}\}_{i=1}^M$ converges to the limit points $\{H_{u_i,v_i}\}_{i=1}^{M}$. Thus, the coupling passes to the limit points as desired.

Proof of Theorem <ref>  (<ref>): This follows from Proposition <ref> (<ref>) immediately in conjunction with the result we already proved in part (<ref>) of this theorem.

Proof of Theorem <ref>  (<ref>): Proposition <ref> (<ref>) implies that the desired equality in law holds along any subsequence $\{N_k\}_{k=1}^{\infty}$ such that $H^{(N_k)}_{u,v}$ and $H^{(N_k)}_{v,u}$ converges to their respective limit points $H_{u,v}$ and $H_{v,u}$. Thus, the coupling passes to the limit points as desired

§ ASKEY-WILSON PROCESS AND THE OPEN ASEP INVARIANT MEASURE

The focus of this paper now pivots to the proof of Theorem <ref> (<ref>). In Section <ref> we start by defining the Askey-Wilson process (Section <ref>) used in the statement of Proposition <ref>; in Section <ref> we record Corollary <ref> of that proposition. Finally, in Section <ref> we define the continuous dual Hahn process which is a special limit of the Askey-Wilson process that arises in our scaling limit.

§.§ Definition of the Askey-Wilson process

Fix $q\in (-1, 1)$ and $a, b, c, d\in \C$ be such that

\begin{equation}\label{eq:assumptabcd}
ac,ad,bc,bd,qac,qad,qbc,qbd,abcd,qabcd \in \C\setminus [1, \infty),
\end{equation}

and such that either all $a, b, c, d$ are real, or two of them are real and two form a complex conjugate pair (e.g. $a,b\in \R$ and $c=\bar{d}$), or they form two complex conjugate pairs (e.g. $a=\bar{c}$ and $b=\bar{d}$).

§.§.§ Askey-Wilson orthogonal polynomials

Following we are interested in the family of Askey-Wilson polynomials $\{w_n\}_{n\in \Z_{\geq 0}}$ defined through the three-term recurrence relation

$$2 x w_n(x)=A_nw_{n+1}+B_n w_n(x)+C_nw_{n-1}(x), \quad n\geq 0, \quad w_{-1}(x)=0, w_0(x)=1,$$

where the coefficients $A_n,B_n$ and $C_n$ are given by

\begin{equation*}
\begin{split}
A_n=&\dfrac{\left(1-a b c d q^{n-1}\right)\left(1-a b q^n \right)}{\left(1-a b c d q^{2n-1}\right)\left(1-a b c d q^{2n}\right)},\\
B_n=&a+\dfrac{1}{a}-\dfrac{\left(1-a b c d q^{n-1}\right)\left(1-a b q^n \right) \left(1-a c q^n \right) \left(1-a d q^n \right)}{a \left(1-a b c d q^{2n-1}\right) \left(1-a b c d q^{2n}\right)}\\
&-a\dfrac{\left(1-q^n \right) \left(1-b c q^{n-1} \right) \left(1-b d q^{n-1}\right) \left(1-c d q^{n-1} \right)}{ \left(1-a b c d q^{2n-2}\right) \left(1-a b c d q^{2n-1}\right)},\\
C_n=&\dfrac{\left(1-q^n \right) \left(1-a c q^{n-1} \right) \left(1-a d q^{n-1} \right) \left(1-b c q^{n-1} \right)  \left(1-b d q^{n-1} \right) \left(1- c d q^{n-1} \right)}{ \left(1-a b c d q^{2n-2}\right) \left(1-a b c d q^{2n-1}\right)}.
\end{split}
\end{equation*}

These remarkable polynomials were introduced in [4] (we use a variant on them from [12]).

§.§.§ Askey-Wilson distribution

Under our assumption on parameters in (<ref>), $\prod_{k=0}^n A_k C_k\geq 0$ for all $n$. In this case the Askey-Wilson polynomials are orthogonal with respect to a unique compactly supported probability measure. For any Borel set $V\subset \R$ that write the measure as $AW(V; a, b, c, d, q)$ and call it the the Askey-Wilson distribution. As shown in  <cit.> or <cit.>, the Askey-Wilson distribution can be decomposed into an absolutely continuous part and a discrete finitely support atomic part.
The density, relative to the Lebesgue measure, of the absolutely continuous part will be denoted by $AW^c$ and the discrete part, which is a sum of weighted Dirac delta functions, will be denoted by $AW^d$. Thus,

\begin{equation}
AW(V; a, b, c, d, q)=\int_{V}AW^c(x; a,b,c,d,q)dx+AW^d(V;a,b,c,d,q).
\end{equation}

The density $AW^c$ is supported on $x\in \Supp^c:=[-1,1]$ and given by

\begin{equation}\label{eq:AW}
AW^c(x;a,b,c,d,q)={\bf 1}_{|x|<1}\frac{(q,ab,ac,ad,bc,bd,cd)_\infty}{2\pi(abcd)_\infty \sqrt{1-x^2}}\bigg\lvert\frac{(e^{2\iu \theta_x})_\infty}{(a
e^{\iu\theta_x}, be^{\iu\theta_x}, c e^{\iu\theta_x}, d e^{\iu\theta_x}
\end{equation}

with $x=\cos\theta_x$. (Recall that we are using the convention of dropping the $q$ in $q$-Pochhammer symbols.)
The discrete part $AW^d$ is given by

\begin{equation}\label{eq:AWddef}
AW^d(V;a,b,c,d,q) = \sum\limits_{y\in V\cap \Supp^d(a,b,c,d,q)}AW^d(y;a,b,c,d,q),
\end{equation}

i.e., it is a sum of Dirac delta functions at locations $y\in \Supp^d(a,b,c,d,q)$ with masses $AW^d(y;a,b,c,d,q)$.
The set $\Supp^d(a,b,c,d,q)$ and the mass $AW^d$ are given as follows. If $|a|, |b|, |c|, |d|<1$ the set $\Supp^d(a,b,c,d,q)$ is empty. By (<ref>), if $\chi \in\{a,b,c,d\}$ has $|\chi|>1$ then it must be real. Each  $\chi \in\{a,b,c,d\}$ with $|\chi|>1$ generate its own set of atoms, the union of which constitutes $\Supp^d(a,b,c,d,q)$. By (<ref>), all elements $\chi\in\{a,b,c,d\}$ with $|\chi|>1$ must be pairwise unique. There are a finite number of atoms generated by such $\chi$ and they are given by

\begin{equation}\label{eq:atomchi}
y^\chi_j=\dfrac{1}{2}\left(\chi q^j+\dfrac{1}{\chi q^j}\right), \text{ for } j\in\Z_{\geq 0} \text{ such that }|\chi q^j|\geq 1.
\end{equation}

Each atom has a different mass. When $\chi=a$, these are given as

\begin{align}\label{eq:AWdmasses}
AW^d(y^a_0;a,b,c,d,q)&=\dfrac{\left(a^{-2}, b c, b d, c d\right)_\infty}{\left(b/a,c/a, d/a, a b c d\right)_\infty},\\ AW^d(y^a_j;a,b,c,d,q)&=AW^d(y^a_0;a,b,c,d,q)\dfrac{\left(a^{2}, a b, a c, a d\right)_j\left(1-a^2 q^{2j}\right)}
{\left(q, q a/b, qa/c, q a/d \right)_j \left(1-a^2\right)} \left(\dfrac{q}{a b c d}\right)^j,
\end{align}

where, in the second line, we assume that $j\in\Z_{\geq 1}$ such that $|\chi q^j|\geq 1$.
For other values of $\chi$, the masses are given as above except with $a$ and $\chi$ swapped.

§.§.§ Askey-Wilson process

Bryc and Wesołowski [12] used the orthogonality of the Askey-Wilson polynomials to construct a class of Markov processes that they called Askey-Wilson processes. These enjoy the characterizing property of being standard quadratic harnesses, and the Askey-Wilson polynomials form a family of martingale orthogonal polynomials for these processes. Bryc and Wesołowski [13] showed that the generating function for the invariant measure of open ASEP could be encoded in terms of this process.

We define the Askey-Wilson process which will relate to open ASEP. This is a Markov process that depends on parameters $A,  B, C, D, q$.
We assume that $A,B,C,D$ correspond to the ASEP parameters $\alpha,\beta,\gamma,\delta$ via (<ref>), in which case $A,C>0$ and $B,D\in (-1,0]$. Additionally, we will assume that $AC<1$, as is necessary for the existence of the Askey-Wilson processes with which we work. This puts us in the fan region of the open ASEP phase diagram (see Section <ref>).

Under the above assumptions on $A,B,C,D$, the Askey-Wilson process is defined to be the Markov process $\{\Y_s\}_{s\in[0, \infty)}$ whose  marginal and transitional distributions are given as follows.
Define the support of the continuous and discrete atomic part at time $s$ as

\begin{equation}\label{eq:Fs}
\Supp^c_s=\Supp^c:=[-1,1], \qquad \Supp^d_s:=\Supp^d\left(A\sqrt{s}, B\sqrt{s}, \dfrac{C}{\sqrt{s}}, \dfrac{D}{\sqrt{s}},q\right)
\end{equation}

where $\Supp^d(a,b,c,d,q)$ is defined as in Section <ref>. The marginal distribution of $\Y_s$ is supported on $\Supp_s:= \Supp^c_s\cup \Supp^d_s$. For any Borel $V\subset\R$, the measure is given by

\begin{equation}\label{eq:def_pi}
\pi_s(V)=AW\left(V; A\sqrt{s}, B\sqrt{s}, \dfrac{C}{\sqrt{s}}, \dfrac{D}{\sqrt{s}}\right).
\end{equation}

The transition distribution from $x\in \Supp_s$ to any Borel $V\subset\R$ is given by

\begin{equation}\label{eq:def_p}
\pi_{s,t}( x, V)=AW \left (V; A\sqrt{t}, B\sqrt{t},\sqrt{ \dfrac{s}{t} }\left(x+\sqrt{x^2-1}\right), \sqrt{ \dfrac{s}{t} }\left(x-\sqrt{x^2-1}\right)\right),
\end{equation}

where, for $x\in\Supp^c$, we interpret $x\pm\sqrt{x^2-1}$ as $e^{\pm \iu \theta_x}$ with $\theta_x$ defined through $x=\cos \theta_x$.

Both $\pi_s$ and $\pi_{s,t}$ have absolutely continuous and discrete atomic parts. For $x\in \Supp^c$ we denote the density of $\pi_s$  by $\pi^c_s(x)$ and for $x\in \Supp^d_s$ we denote the mass $\pi_s(\cdot)$ assigns to $x$ by $\pi^d_s(x)$.

For the transition distribution, there are four cases:

* If $x\in \Supp^c$ and $y\in \Supp^c$ then we write $\pi^{c,c}_{s,t}(x,y)$ for the density  in $y$.

* If $x\in \Supp^c$ and $y\in \Supp^d_t$ then we write $\pi^{c,d}_{s,t}(x,y)$ for the mass assigned to $y$.

* If $x\in \Supp^d_s$  and $y\in \Supp^c$ then we write $\pi^{d,c}_{s,t}(x,y)$ for the density in $y$

* If $x\in \Supp^d_s$ and $y\in \Supp^d_t$ then we write $\pi^{d,d}_{s,t}(x,y)$ for the mass assigned to $y$.

For all other values of $x$ and $y$, we declare that these functions are zero.

The existence and uniqueness of a Markov process with these marginals and transition probabilities is shown in <cit.>. In particular, <cit.> shows that for all Borel sets $V\subset \R$ (and in the second equation, for all $x\in \Supp_s$)

\begin{align}
\int_{\R}\pi_s(dx)\pi_{s,t}(x,V) = \pi_t(V),\qquad \textrm{and} \qquad
\int_{\R}\pi_{s,t}(x,dy)\pi_{t,u}(y,V) = \pi_{s,u}(x,V).
\end{align}

The first identity implies that the Askey-Wilson process started under $\pi_s$ at time $s$ has marginal distribution $\pi_t$ at time $t$, while the second is the Chapman-Kolmogorov identity needed to define the Markov process. We will only be dealing with the Askey-Wilson process on a finite time interval and we assume that at the starting time $s$ the Markov chain is initialized according to the marginal distribution $\pi_s$. Then, for all later times $t$, the marginal distribution will necessarily be $\pi_t$.

We remark that the Askey-Wilson process is defined in terms of its consistent family of finite-dimensional distribution. We only rely on these finite dimensional distributions and not on the existence of a continuous version of the process. In fact, the continuity of sample paths is complicated by the discrete support of the measures and should involve a mixture of continuous trajectories within the continuous support, and jump processes in the discrete support.

§.§ Generating function for open ASEP invariant measure

Proposition <ref> contains a concise formula for the open ASEP invariant measure generating function in terms of the above defined Askey-Wilson process. As a corollary of the result, we record here a formula for the Laplace transform of the scaled stationary measure for the open ASEP height function $\hscale(x)$, see (<ref>).

As in (<ref>), let $ \big \langle \cdot  \big\rangle_N$ denote the expectation with respect to the invariant measure $\pi^{\mathrm{ASEP}}_N$ of open ASEP parameterized by $q$ and $(A,B,C,D)$. Assume that $AC<1$.
For any $d\in \Z_{\geq 1}$, let $\vec{X},\vec{c}$ and $\vec{s}$ be as in (<ref>), $\tilde{c}\in\R$. Then

\begin{equation}\label{eq:pihn}
\phi^{(N)}(\vec{c},\tilde{c},\vec{X}):= \left\langle e^{-\sum_{k=1}^d c_k\hscale(X^{(N)}_k) -\tilde{c}\hscale(1)} \right\rangle_N = \dfrac{\displaystyle\E\left[\,\,\prod_{k=1}^{d+1}\left(\cosh\left(\tilde{s}^{(N)}_k\right)+
\Y_{e^{-2 \tilde{s}^{(N)}_k}}\right)^{n_k-n_{k-1}} \right]}{\E  \left [\left (1+\Y_1 \right )^N\right ]}
\end{equation}

where $\tilde{s}^{(N)}_k = N^{-1/2} (s_k + \tilde{c})$ for $k\in \llbracket 1,d+1\rrbracket$, $X_k^{(N)}:= N^{-1} \lfloor NX_k\rfloor$, $n_k=\lfloor NX\rfloor$  for $k\in \llbracket 1,d\rrbracket$, and $\Y$ is the Askey-Wilson process from Definition <ref> with parameters $A,B,C,D,q$ matching those of the ASEP we are considering. In the special case when $\tilde{c}=0$ we will write $\phi^{(N)}(\vec{c},\vec{X})$ instead of $\phi^{(N)}(\vec{c},0,\vec{X})$ and use $s_k$ instead of $\tilde{s}_k$ on the right-hand side of (<ref>).

Recall that $\hscale(X)$ is written in terms of the occupation variables $\tau_1,\ldots, \tau_N$ (distributed according to $\pi^{\mathrm{ASEP}}_N$) in such a way that for $X\in [0,1]\cap \Z/N$, $\hscale(X) = N^{-1/2} \sum_{i=1}^{NX} (2\tau_i -1)$. Assume that $X_k\in \Z/N$ for $k\in \llbracket 1,d\rrbracket$ and recall that $s_{d+1}= 0$ hence $\tilde{s}_{d+1}=\tilde{c}$.  Then,

\begin{align}\label{eq:lalace_height}
\phi^{(N)}(\vec{c},\tilde{c},\vec{X})&:=\Bigg\langle \exp\bigg(-\sum_{k=1}^{d+1}\,\,\,
\sum_{j=n_{k-1}+1}^{n_k}(2\tau_j-1)\tilde{s}^{(N)}_k\bigg)\Bigg \rangle_N\\
&=\exp\left(\,\,\sum_{k=1}^{d+1} \tilde{s}^{(N)}_k(n_k-n_{k-1}) \right) \Bigg\langle
\prod_{k=1}^{d+1}\,\,\prod_{j=n_{k-1}+1}^{n_k}\left(e^{-2 \tilde{s}^{(N)}_k}\right)^{\tau_j}
\Bigg \rangle_N\\
&=\exp\left(\,\,\sum_{k=1}^{d+1} \tilde{s}^{(N)}_k(n_k-n_{k-1}) \right) \dfrac{\displaystyle\E\left[\,\,\prod_{k=1}^{d+1}\left(1+e^{-2\tilde{s}^{(N)}_k}+2 e^{-\tilde{s}^{(N)}_k}
\Y_{e^{-2 \tilde{s}^{(N)}_k}}\right)^{n_k-n_{k-1}} \right]}{2^N \E
\left [\left (1+\Y_1 \right )^N\right ]}\\
\dfrac{\displaystyle\E\left[\,\,\prod_{k=1}^{d+1}\left(\cosh\left(\tilde{s}^{(N)}_k\right)+
\Y_{e^{-2 \tilde{s}^{(N)}_k}}\right)^{n_k-n_{k-1}} \right]}{\E  \left [\left (1+\Y_1 \right )^N\right ]},
\end{align}

where we used Proposition <ref> between the second and third line. The extension to general $\vec{X}$ amounts to just rounding each $X_k$ down to the closest multiple of $N^{-1}$.

For the rest of the paper we will only work with the $\tilde{c}=0$ case of Corollary <ref> since this is sufficient for our purposes. It should be possible to use the general $\tilde{c}$ result and the connection to Askey-Wilson processes to extract asymptotics involving the joint distribution of the increments of the open KPZ invariant measure along with its value at $X=1$. We do not pursue that here.

Notice also that in Corollary (<ref>) we require that the Laplace variables $c_k$ are strictly positive, recall (<ref>). The formula trivially extends when some $c_k=0$ since that just decreases the value of $d$. However, we do not know how to prove a formula of this type when any of the $c_k$ are strictly negative. This limitation on parameters comes from the fact that, in Proposition <ref>, the $t$'s need to be increasing. That restriction, in turn, comes from the proof of this result from [13] which relies on using conditional expectations for the Askey-Wilson process, and the orthogonal martingale property for Askey-Wilson polynomials with respect to this process. We attempted to generalize this argument to allow for time reversal, but had no luck. It is an interesting question as to whether a similar formula as in Proposition <ref> holds when the $c_k$ are allowed to be negative. We do know that the Laplace transform is defined for all $c_k$ and is analytic in those variables. The question is whether the Askey-Wilson formulas (which are written in terms of certain integrals) have a nice analytic continuation which can be similarly explicit.
Finally, we note that despite the limitation on our variables $c_k$ needing to be positive, this will not prevent us from deducing convergence in distribution of $\hscale$ as $N\to \infty$.

§.§ The continuous dual Hahn process

The open KPZ stationary measures that we construct in Theorem <ref> (<ref>) is characterized by a duality with the continuous dual Hahn process that we define here (the special case of this when $u,v>0$ was already briefly defined in Section <ref>). This process is a certain limit of the Askey-Wilson processes [13] (see Section <ref>).

We will use the standard notation for the Pochhhamer symbol: For $j\in \Z_{\geq 0}$ and $x\in \R$, define $[x]_j:=(x)(x+1)\cdots(x+j-1)$ with the convention that $[x]_0:=1$. For multiple arguments $x_1,\ldots x_n\in \R$, define $[x_1, \ldots, x_n]_j:=[x_1]_j\cdots [x_n]_j$.

To define the continuous dual Hahn process, we will give its marginal and transition probability distributions, and then show that they are consistent. These distributions all of related to generalized beta integrals in [2, 51]. We start by defining the orthogonality probability measure for the continuous dual Hanh orthogonal polynomials. Its definition depends on the nature of its parameters, so we must consider a few cases. In particular, there may be an absolutely continuous part $\CDHtransition^c$ and a discrete atomic part $\CDHtransition^d$ to the measure.

Assume that $a, b, c\in \R$ with $a+b, a+c>0 $ or  $a\in \R$ and $b=\bar c\in
\C\setminus \R$ with $\textup{Re}(b)=\textup{Re}(c)>0$. For such $a,b,c$, define

\begin{equation}
\CDHtransition^c (x; a, b, c):=\dfrac{1}{8\pi}\,\cdot\,\dfrac{\bigg\rvert \Gamma\left(a+\iu  \dfrac{\sqrt x}{2}\right)\cdot\Gamma\left(b+\iu
\dfrac{\sqrt x}{2}\right)\cdot\Gamma\left(c+\iu  \dfrac{\sqrt x}{2}\right)\bigg \rvert^2}{\Gamma(a+b)\cdot \Gamma(a+c)\cdot\Gamma (b+c)\cdot \sqrt x\cdot\left\rvert \Gamma\left(\iu \sqrt x\right)\right\rvert^2}{\bf 1}_{x> 0}.
\end{equation}

Assume that $a<0$ and that $\SuppCDH$ is a finite subset of $\R$
with size $|\SuppCDH| = \lfloor -a\rfloor+1$ and elements $x_0<
x_1<\cdots<x_{\lfloor -a\rfloor}$. Assume that $b,c$ and $j$ satisfy one of the three conditions:

* $b=\bar c\in \C\setminus \R$ with
$\textup{Re}(b),\textup{Re}(c)>0$, and $j\in \llbracket 0,\lfloor -a\rfloor\rrbracket$;

* $b, c\in \R$ with
$a+b, a+c>0 $, and $j\in \llbracket 0,\lfloor -a\rfloor\rrbracket$;

* $b, c\in \R$ and $b, b+c, c-a>0,$ with $a+b=-k$ for $k\in \Z_{\geq 0}$, and $j\in\llbracket 0,k\rrbracket$.

For $a,b,c$,  $\SuppCDH$, and $j$ as above, define

\begin{equation}
\CDHtransition^d (x_j; a, b, c;\SuppCDH):=\dfrac{[2 a,a+b,a+c]_j\cdot(a+ j)\cdot
\Gamma(b-a)\cdot\Gamma(c-a)}{[1,a-b+1,a-c+1]_j \cdot a \cdot\Gamma(-2
\end{equation}

For all other first arguments, define $\CDHtransition^d (\cdot; a, b, c;\SuppCDH)$ as zero.

We will define the measure $\CDHtransition(V;a,b,c;\SuppCDH)$ (for Borel subsets $V$ of $\R$) under three possible sets of conditions on parameters $a,b,c$ and sets $\SuppCDH$:

P: For $a\geq 0$; $\SuppCDH=\varnothing$; and either $b=\bar c\in \C\setminus \R$ with
$\textup{Re}(b),\textup{Re}(c)>0$ or $b, c\in \R_{>0}$, define

\begin{equation}
\CDHtransition(V; a,b,c;\SuppCDH):=\int_V \CDHtransition^c (x; a, b, c)dx.
\end{equation}

N1: For $a<0$; $\SuppCDH$ a finite subset of $\R$
with size $|\SuppCDH| = \lfloor -a\rfloor+1$ and elements $x_0<x_1<\cdots<x_{\lfloor -a\rfloor}$; and either $b=\bar c\in \C\setminus \R$ with
$\textup{Re}(b),\textup{Re}(c)>0$  or $b, c\in \R$ with $a+b, a+c>0$, define

\begin{equation}
\CDHtransition(V;a,b,c;\SuppCDH):=\int_V \CDHtransition^c (x; a,
b, c)dx+\sum\limits_{x\in \SuppCDH\cap V}\CDHtransition^d(x; a,b,c;\SuppCDH)
\end{equation}

N2: For $a<0$;  $\SuppCDH$ a finite subset of $\R$
with size $|\SuppCDH| = \lfloor -a\rfloor+1$ and elements $x_0<x_1<\cdots<x_{\lfloor -a\rfloor}$; and $b, c\in \R$ with $b, b+c, c-a>0$ and $a+b=-k$ for $k\in \Z_{\geq 0}$, define

\begin{equation}
\CDHtransition(V;a,b,c;\SuppCDH):=\sum\limits_{x\in \SuppCDH\cap V}\CDHtransition^d(x; a,b,c;\SuppCDH).
\end{equation}

In all three cases  of Definition <ref>, $\CDHtransition(\cdot; a,b,c;\SuppCDH)$ is a probability measure on $\R$.

This follows from the identities (3.1), (3.3) and (3.4) in
[51] as a limit when Wilson polynomials
degenerate to dual continuous Hahn polynomials. In particular, this limit corresponds to taking one of the parameters of the Wilson polynomials to infinity. The case when $a=0$ does not seem to be covered therein, but can be recovered by taking the limit as $a\to 0$.

The $\CDHtransition$ distribution will be used to define the transition probability distribution for the continuous dual Hahn process. Before doing that, we will define what will be its marginal distributions. We will only define the continuous dual Hahn process for times in a certain interval $[0,\Cuv)$ where we recall from (<ref>) that $\Cuv=2$ if $\vv\leq 0$ or $\vv\geq1$ and $\Cuv=2\vv$ if $\vv\in (0,1)$.
For all $\vv\in \R$, $[0,\Cuv)$ has a non-empty interior.
It should be possible to extend the definition to all times in $\R$, though it would require introducing more complicated distributions. For the purpose of stating Theorem <ref> (<ref>), we only need to have the process defined on some non-empty interval.