\section{Comparison to previous works}\label{Section9} In this section we give a formal comparison between the arguments in the present paper and those in \cite{Dot13,PSpohn,ISS13} and \cite{NZ}. The purpose of this section is to help readers familiar with some of these works relate them to this paper, but also to present some of our arguments in an informal and more accessible way.


\subsection{Prelimit formulas via difference operators}\label{Section9.1} In Section \ref{Section9.1.1} we informally apply the method of the Macdonald difference operators and derive a prelimit formula given in (\ref{GenFunExpand3}), comparing the result with formula (\ref{GenFunForm}), which was obtained for the KPZ equation in \cite{Dot13,PSpohn,ISS13}. The discussion in this section is more conceptual and we do not discuss issues of convergence. In Section \ref{Section9.1.2} we explain how to rewrite (\ref{GenFunExpand3}) in terms of contour integrals and match our prelimit formula from Section \ref{Section3}.


\subsubsection{KPZ equation}\label{Section9.1.1}
The two-point large time distribution of the KPZ equation \ref{KPZEq} started from narrow wedge initial data was investigated independently in \cite{Dot13} and \cite{PSpohn}, and the authors obtained different formulas, which were ultimately reconciled in \cite{ISS13}. In this section we discuss the approach in these papers and compare it to the one in our paper. The exposition below follows \cite{ISS13}.

Let $\mathcal{H}(X,T)$ denote the solution to the KPZ equation \ref{KPZEq} with narrow wedge initial data. The Cole-Hopf transform $\mathcal{Z} = e^{\mathcal{H}}$ becomes a solution to the {\em stochastic heat equation}
\begin{equation}\label{SHEEq}
\partial_T\mathcal{Z}(X,T) = \frac{1}{2} \partial^2_X \mathcal{Z}(X,T)+  \xi (X,T)\mathcal{Z}(X,T),
\end{equation}
with initial condition $\mathcal{Z}(x,0) = \delta(x)$. Let us set $\gamma_T = (T/2)^{1/3}$ and define $\tilde{\mathcal{H}}$ through
$$\mathcal{H}(X,T) = - \frac{T}{24} + \gamma_T \tilde{\mathcal{H}}(X,T).$$
We also introduce the generating function 
\begin{equation}\label{GenFun}
\begin{split}
&G_{T,X_1, X_2}(s_1, s_2) = \mathbb{E} \left[ \exp \left( - e^{\gamma_T \tilde{\mathcal{H}}( X_1,T) - s_1 }- e^{\gamma_T \tilde{\mathcal{H}}( X_2,T) - s_2 }\right)\right] = \\
& \mathbb{E} \left[ \exp \left( - e^{(T/24) - s_1} \mathcal{Z}(X_1,T)- e^{(T/24) - s_2} \mathcal{Z}(X_2, T) \right) \right].
\end{split}
\end{equation}
Provided the limit exists, we have the following asymptotic equality
\begin{equation}\label{S1lim}
\begin{split}
&\lim_{T \rightarrow \infty} G_{T,2\gamma_T^2 X_1, 2\gamma_T^2 X_2}(\gamma_T (s_1 - X_1^2), \gamma_T(s_2 - X_2^2) )= \\
&\lim_{T \rightarrow \infty} \mathbb{P} \left( \tilde{\mathcal{H}}(X_1, T) + X_1^2 \leq s_1, \tilde{\mathcal{H}}(X_2, T) + X_2^2 \leq s_2 \right).
\end{split}
\end{equation}
Equation (\ref{S1lim}) reduces the question of understanding the two-point large time distribution of the KPZ equation to understanding the limit of the generating function $G_{T,X_1, X_2}(s_1, s_2) $. In \cite{Dot13} and \cite{PSpohn} the authors found different ways to rewrite this generating series so that the limit can be studied. 

In the present paper we consider the stochastic six vertex model, which can be thought of as a discrete analogue of the KPZ equation. Indeed, the convergence of the height function $h$ of the six vertex model to the KPZ equation under a weak asymmetric scaling ($b_2 \in (0,1)$-fixed and $b_1 \rightarrow b_2^-$) was recently proved in \cite{CGST} for several, but in particular step, initial conditions. We remark that earlier in \cite[Theorem 12.3]{BO17} the one-point marginals of the height function $h$ were shown to converge to the KPZ equation in the limit $b_1, b_2 \rightarrow 1^-$. For the stochastic six-vertex model $\mathcal{P}(b_1, b_2)$ from Section \ref{Section1.2} we have the following formal correspondence with the above notation:
\begin{equation}\label{Corr}
\begin{split}
&M \leftrightarrow T,\hspace{5mm} M - h(\cdot, M) \leftrightarrow \mathcal{H}(\cdot, T), \hspace{5mm} t^{h(\cdot , M) - M} \leftrightarrow \mathcal{Z}(\cdot, T)\\
&\mathbb{E}\left[ \frac{1}{(u_1t^{h(n_1 + 1 , M) - M};t)_{\infty}} \frac{1}{(u_2t^{h(n_2 + 1 , M) - M};t)_{\infty}} \right] \leftrightarrow G_{T,X_1, X_2}(s_1, s_2) .
\end{split}
\end{equation}
We recall that $t = b_1/b_2$ and the weak asymmetric limit taking the six-vertex model to the KPZ equation corresponds to taking $t \rightarrow 1^-$. In view of the above correspondence, we see that in our paper we are in a sense studying the asymptotics of a discrete analogue of the KPZ equation, by taking the limit of a discrete analogue of the joint Laplace transform.\\

The first step in rewriting $G_{T,X_1, X_2}(s_1, s_2)$ in a form suitable for asymptotics is to form the moment expansion
\begin{equation}\label{GenFunExpand}
\begin{split}
&G_{T,X_1, X_2}(s_1, s_2) = \sum_{L = 0}^\infty \sum_{R =0}^\infty \frac{(-e^{-s_1 +T/24} )^L}{L!} \frac{(-e^{-s_2 + T/24})^R}{R!}  \mathbb{E} \left[ \mathcal{Z}(X_1,T)^L \mathcal{Z}(X_2,T)^R \right].
\end{split}
\end{equation}
Equation (\ref{GenFunExpand}) is already mathematically ill-posed since $\log \mathbb{E} \left[\mathcal{Z}(X,T)^n \right] \sim n^3$, which means that the radius of convergence of the above expansion in the variables $-e^{-s_1 + T/24}$ and  $-e^{-s_2 + T/24}$ is zero. The issue of formulas involving divergent sums is ubiquitous in the physics replica approach and we will ignore any convergence issues from here on out when we discuss it.

The analogue of (\ref{GenFunExpand}) in our setting is given by
\begin{equation}\label{GenFunExpand2}
\begin{split}
&\mathbb{E}\left[ \frac{1}{(u_1t^{h(n_1 + 1 , M) - M};t)_{\infty}} \frac{1}{(u_2t^{h(n_2 + 1 , M) - M};t)_{\infty}} \right]= \\
&\sum_{L = 0}^\infty \sum_{R =0}^\infty \frac{u_1^L(1-t)^{-L}}{L_t!}\frac{u_2^R(1-t)^{-R}}{R_t!} \mathbb{E}\left[ t^{L (h(n_1 + 1 , M) - M)} t^{R (h(n_2 + 1 , M) - M)}\right],
\end{split}
\end{equation}
where we used the power series expansion of the $t$-exponential function (\ref{S1tExpon}). Note that equation (\ref{GenFunExpand2}) as a power series in $u_1, u_2$ has a positive radius of convergence, since $h \in [0,M]$ with probability $1$, making the above moments bounded by $t^{-M(R+L)}$. In particular, the moment expansion in our discrete setting is rigorous and one can view the expansion in (\ref{GenFunExpand}) as a formal $t \rightarrow 1^-$ limit of the one in (\ref{GenFunExpand2}).\\

Both in our paper and in \cite{ISS13} the essential goal is to rewrite the $(L,R)$-th summand in (\ref{GenFunExpand2}) and (\ref{GenFunExpand}) respectively, and rearrange the sum in a form that is suitable for asymptotic analysis. The way this is achieved in \cite{ISS13} is by utilizing the Feynman-Kac formula, which implies that the $n$-point correlation function $\mathbb{E} \left[ \mathcal{Z}(X_1, t) \cdots \mathcal{Z}(X_n,T)\right]$ satisfies the {\em imaginary-time Schr{\"o}dinger} equation
$$ -\partial_T \mathbb{E} \left[ \mathcal{Z}(X_1, t) \cdots \mathcal{Z}(X_n,T)\right] = H_n \mathbb{E} \left[ \mathcal{Z}(X_1, T) \cdots \mathcal{Z}(X_n,T)\right],$$
with initial condition $\mathcal{Z}(X_i,0) = \delta(X_i)$, where $H_n$ denotes the Lieb-Liniger quantum Hamiltonian of $n$ particles on the line with attractive $\delta$-interaction \cite{Lieb, McGuire}:
$$H_n = -\frac{1}{2} \sum_{i = 1}^n \partial_{X_i}^2 - \frac{1}{2} \sum_{i \neq j = 1}^n \delta (X_i - X_j).$$
By using the Bethe ansatz \cite{Dot10} derived a set of orthonormal eigenfunctions and corresponding eigenvalues for the Lieb-Liniger quantum Hamiltonian $H_n$. We mention that the orthonormal set found in \cite{Dot10} for the attractive case is not proved therein to be complete, see the last paragraph in \cite[Section B.3]{Dot10}; however, completeness was established in \cite{HO97, Oxford79, PS11}. Once the eigenbasis is obtained, one can express the $(L,R)$-th summand in (\ref{GenFunExpand}) in terms of it and rearrange the resulting sum to obtain
\begin{equation}\label{GenFunForm}
\begin{split}
&G_{T,X_1, X_2}(s_1, s_2) = \sum_{M = 0}^\infty \frac{(-1)^M}{M!} \prod_{\alpha = 1}^M \int_{\mathbb{R}} \frac{dq_\alpha}{2\pi} \sum_{m_\alpha + l_\alpha > 0} (-1)^{m_\alpha + l_\alpha - 1} \binom{m_\alpha + l_\alpha}{m_\alpha} \times \\
&\det \left[ \frac{1}{\frac{1}{2} (m_\alpha + l_\alpha) + \frac{1}{2} (m_\beta + l_\beta) + \iota (q_\alpha - q_\beta)}\right]_{\alpha, \beta = 1}^M \times \\
& e^{(m_\alpha + l_\alpha)^3 T/24 - \frac{1}{2}(m_\alpha + l_\alpha)q_\alpha^2T - m_\alpha s_1 - l_\alpha s_2 + \iota X_1 m_\alpha q_\alpha + \iota X_2 l_\alpha q_\alpha - m_{\alpha} l_\alpha(X_2 - X_1)/2} \times CT_{KPZ},
\end{split}
\end{equation}
where 
\begin{equation}\label{S1CT2}
\begin{split}
&CT_{KPZ} = \prod_{\alpha \neq \beta} \frac{\Gamma\left[1 + \frac{1}{2} (m_\alpha + A_{\alpha, \beta} + l_\beta) +B_{\alpha, \beta}  \right] \Gamma\left[1 + \frac{1}{2} (-m_\alpha + A_{\alpha, \beta} - l_\beta) + B_{\alpha, \beta}\right]}{\Gamma\left[1 + \frac{1}{2} (-m_\alpha + A_{\alpha, \beta} + l_\beta) +B_{\alpha, \beta}  \right] \Gamma\left[1 + \frac{1}{2} (m_\alpha + A_{\alpha, \beta} - l_\beta) + B_{\alpha, \beta}\right]},
\end{split}
\end{equation}
and we have set $A_{\alpha, \beta} = m_\beta + l_\alpha$, $B_{\alpha, \beta} = \iota (q_\alpha - q_\beta)$. 

We will not have much to say about equation (\ref{GenFunForm}) apart from the fact that the eigenbasis for the Lieb-Liniger quantum Hamiltonian $H_n$ is labeled by complex momenta, which have continuous real parts (labeled by $q_\alpha$'s) and discrete imaginary parts (labeled by $n_{\alpha}$'s) and the latter upon re-groupings and symmetrizations become the discrete labels $m_\alpha$'s and $l_\alpha$'s in (\ref{GenFunForm}). Correspondingly, the summations and integrals in (\ref{GenFunForm}) represent a decomposition of (\ref{GenFunExpand}) onto the eigenbasis. \\

The approach in the present paper formally corresponds to rewriting the $(L,R)$-th summand in (\ref{GenFunExpand2}) using the relation
\begin{equation}\label{S1MomentFun}
\begin{split}
\mathbb{E}\left[ t^{L (h(n_1 + 1 , M) - M)} t^{R (h(n_2 + 1 , M) - M)}\right] = \frac{\mathcal{D}^R_{n_2}\mathcal{D}^L_{n_1} \Pi(X;Y)}{ \Pi(X;Y)} \Big{\vert}_{x_i =a = y_j},
\end{split}
\end{equation}
where $a^2 = \frac{1- b_2}{1-b_1}$, $\mathcal{D}_{n}$ is an affine shifted degeneration of the Macdonald difference operator acting on the variables $(x_1, \dots, x_n)$, and $\Pi(X;Y) =\prod_{i = 1}^N \prod_{j = 1}^M \frac{1 - t x_i y_j}{1 - x_iy_j}$. The reason we use the term ``formally'' here is that one needs to assume for example that $u_1, u_2$ are suffuciently small in (\ref{GenFunExpand2}), and then do certain analytic continuation arguments. For the purposes of this section, we will ignore these convergence issues and continue with the discussion.

The right side of (\ref{S1MomentFun}) can be written as a $(L+R)$-fold nested contour integral using Proposition \ref{milestone1}. The result is 
\begin{equation}\label{S1MomentFun2}
\begin{split}
&\mathbb{E}\left[ t^{L (h(n_1 + 1 , M) - M)} t^{R (h(n_2 + 1 , M) - M)}\right] = \frac{1}{(2\pi \iota)^{R+L}} \int_{C^1_{0,1}} \cdots \int_{C^1_{0,L}}  \int_{C^2_{0,1}} \cdots \int_{C^2_{0,R}}\\
&   \prod_{1 \leq \alpha < \beta \leq L}  \frac{\omega_\alpha - \omega_\beta}{\omega_\alpha - \omega_\beta t^{-1}} \prod_{1 \leq \alpha < \beta \leq R}  \frac{\zeta_\alpha - \zeta_\beta}{\zeta_\alpha - \zeta_\beta t^{-1}}  \prod_{i =1}^L \left( \frac{1- a\omega_i^{-1}t^{-1}}{1 - a\omega_i^{-1}} \right)^{n_1}\left(\frac{1 - a\omega_i}{1 - at\omega_i} \right)^M \\
& \prod_{i =1}^R \left( \frac{1- a\zeta_i^{-1}t^{-1}}{1 - a\zeta_i^{-1}} \right)^{n_2} \left( \frac{1 - a\zeta_i}{1 - at\zeta_i} \right)^M \prod_{i = 1}^L\prod_{j =1}^R    \frac{\omega_i - \zeta_j}{\omega_i- \zeta_jt^{-1}} \prod_{i =1}^R\frac{d\zeta_i}{\zeta_i}\prod_{i =1}^L\frac{d\omega_i}{\omega_i}.
\end{split}
\end{equation}
In (\ref{S1MomentFun2}) the contours are nested zero-centered positively oriented circles with $C_{0,r}^i$ containing $t^{-1} \cdot C_{0,{r+1}}^i$ and $C_{0,1}^2$ contained in $C^1_{0,L}$. In addition, one has that $a$ needs to be inside the smallest contour $C^2_{0,R}$ and $a^{-1}$ needs to be outside the biggest contour $C_{0,1}^1$. 

We next proceed to deform the $C_{0,i}^1$ contours to $C_{0,L}^1$, which can be taken to be a circle of radius $1+\varepsilon$ and also we deform the $C_{0,i}^2$ contours to $C_{0,1}^2$, which can be taken to be a circle of radius $1 - \varepsilon$. In plain words we are deforming the inner circles $C_{0,i}^2$ in the outward direction and the outer circles $C_{0,i}^1$ in the inward direction and the two almost meet in the middle. We explain in Section \ref{Section9.1.2} why we make this particular choice of deformation of the contours. As we deform the contours, we pick up poles from $\zeta_i = \zeta_jt^{-1}$ and $\omega_i = \omega_jt^{-1}$ in (\ref{S1MomentFun2}). The result is that the dimension of the contour integrals goes down and we end up integrating over residue subspaces, which for the $\omega$-contours are labelled by partitions $\lambda \vdash L$ and for the $\zeta$-contours by partitions $\mu \vdash R$. This contour deformation procedure is an instance of the contour integral ansatz, see Lemma \ref{nestedContours} in the main text (its origin is \cite[Proposition 7.2]{BBC}). The result of applying this ansatz is
\begin{equation}\label{S1MomentFun3}
\begin{split}
&\mathbb{E}\left[ t^{L (h(n_1 + 1 , M) - M)} t^{R (h(n_2 + 1 , M) - M)}\right] = \sum_{\lambda \vdash L} \sum_{\mu \vdash R} \frac{(t^{-1} - 1)^{L}(1-t)^RR_t! L_t!}{\prod_{i = 1}^\infty m_i(\lambda)! m_i(\mu)!}  \int_{C_{1 + \varepsilon}^{\ell(\lambda)}} \int_{C_{1-\varepsilon}^{\ell(\mu)}} \\
&  \det \left[ \frac{1}{w_it^{-\lambda_i} - w_j}\right]_{i,j = 1}^{\ell(\lambda)} \prod_{i = 1}^{\ell(\lambda)} \left( \frac{1 - aw_it^{- \lambda_i}}{1 - aw_i } \right)^M  \left(\frac{1 - a w_i^{-1} }{1 - a w_i^{-1}  t^{\lambda_i}}\right)^{n_1} \\
&  \det \left[ \frac{1}{\hat{z}_i - \hat{z}_j t^{\mu_j}}\right]_{i,j = 1}^{\ell(\mu)} \prod_{i = 1}^{\ell(\mu)} \left( \frac{1 - a\hat{z}_i }{1 - a \hat{z}_it^{\mu_i}} \right)^M    \left( \frac{1 - a\hat{z}^{-1}_it^{-\mu_i}}{1 - a\hat{z}^{-1}_i } \right)^{n_2} \\
& \prod_{i = 1}^{\ell(\lambda)}  \prod_{j = 1}^{\ell(\mu)}   \frac{(\hat{z}_j w_i^{-1} t^{\lambda_i}; t)_\infty }{(\hat{z}_j w_i^{-1} t^{ \lambda_i + \mu_j}; t)_\infty }\frac{(\hat{z}_j w_i^{-1} t^{\mu_j}; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty } \prod_{i = 1}^{\ell(\mu)} \frac{d \hat{z}_i}{2\pi \iota}\prod_{i = 1}^{\ell(\lambda)} \frac{d w_i}{2\pi \iota},
\end{split}
\end{equation}
where $k_t! = \frac{(1-t)(1-t^2) \cdots (1- t^k)}{(1-t)^k}$ and $m_i(\lambda)$ denotes the multiplicity of $i$ in the partition $\lambda$ (see Section \ref{Section2.1} for more notation about parititons). The contours $C_{1-\varepsilon}$ and $C_{1+ \varepsilon}$ are positively oriented, zero-centered circles of radius $1-\varepsilon$ and $1 + \varepsilon$ respectively. A priori, equation (\ref{S1MomentFun3}) holds only if $a$ is sufficiently close to zero (depending on $L$ and $R$), because only then is equation (\ref{S1MomentFun2}) valid. However, the immediate benefit of performing the contour deformation is that the right side of (\ref{S1MomentFun3}) is analytic in $a$ and then the equality in (\ref{S1MomentFun3}) extends to all $a \in \mathbb{C}$ such that $|a|^{-1} > 1 + \varepsilon$ and $|a| < 1 - \varepsilon$ by analyticity (the pole $w_i = a^{-1}$ needs to be outside the contour $C_{1 +\varepsilon}$ and the pole $\hat{z}_i = a$ needs to be inside the contour $C_{1 - \varepsilon}$). 

We may now symmetrize (\ref{S1MomentFun3}) in the parts $(\lambda_1, \dots, \lambda_{\ell(\lambda)})$ and $(\mu_1, \dots, \mu_{\ell(\mu)})$ -- observe that there are precisely $\frac{\ell(\lambda)! \ell(\mu)!}{\prod_{i = 1}^\infty m_i(\lambda)! m_i(\mu)!}$ ways to do this. Afterwards we perform the summation in (\ref{GenFunExpand2}) and rearrange the resulting sum according to the values $N_1 = \ell(\lambda)$ and $N_2 = \ell(\mu)$. The result is
\begin{equation}\label{GenFunExpand3}
\begin{split}
&\mathbb{E}\left[ \frac{1}{(u_1t^{h(n_1 + 1 , M) - M};t)_{\infty}} \frac{1}{(u_2t^{h(n_2 + 1 , M) - M};t)_{\infty}} \right]=\sum_{N_1 = 0}^\infty \sum_{N_2 = 0}^{\infty} \frac{1}{N_1! N_2!} \sum_{\lambda_1, \dots, \lambda_{N_1} = 1}^\infty  \sum_{\mu_1, \dots, \mu_{N_2} = 1}^\infty\\
&\int_{C_{1 + \varepsilon}^{N_1}} \int_{C_{1 - \varepsilon}^{N_2}}   \det \left[ \frac{1}{w_it^{-\lambda_i} - w_j}\right]_{i,j = 1}^{N_1} \prod_{i = 1}^{N_1} \left( \frac{1 - aw_it^{- \lambda_i}}{1 - aw_i } \right)^M  \left(\frac{1 - a w_i^{-1} }{1 - a w_i^{-1}  t^{\lambda_i}}\right)^{n_1} [u_1t^{-1}]^{\lambda_i} \\
&  \det \left[ \frac{1}{\hat{z}_i- \hat{z}_j t^{\mu_j}}\right]_{i,j = 1}^{N_2} \prod_{i = 1}^{N_2} \left( \frac{1 - a\hat{z}_i }{1 - a \hat{z}_it^{\mu_i}} \right)^M    \left( \frac{1 - a\hat{z}^{-1}_it^{-\mu_i}}{1 - a\hat{z}^{-1}_i } \right)^{n_2}  [u_2]^{\mu_i}\\
& \prod_{i = 1}^{N_1}  \prod_{j = 1}^{N_2}   \frac{(\hat{z}_j w_i^{-1} t^{\lambda_i}; t)_\infty }{(\hat{z}_j w_i^{-1} t^{ \lambda_i + \mu_j}; t)_\infty }\frac{(\hat{z}_j w_i^{-1} t^{\mu_j}; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty } \prod_{i = 1}^{N_2} \frac{d \hat{z}_i}{2\pi \iota}\prod_{i = 1}^{N_1} \frac{d w_i}{2\pi \iota},
\end{split}
\end{equation}
where the $N_1 = N_2 = 0$ term in (\ref{GenFunExpand3}) is just $1$. 

At this time, we can do a formal comparison between (\ref{GenFunExpand3}) and (\ref{GenFunForm}). Firstly, the integration over $w_i$ and $\hat{z}_i$ roughly corresponds to the integration over $q_\alpha$, while the summation over $\mu_i$ and $\lambda_i$ roughly corresponds to the one over $m_\alpha$ and $l_\alpha$. The cross term that appears in the last line and the two Cauchy determinants in (\ref{GenFunExpand3}) can be rewritten as
\begin{equation}\label{CrossV1}
\begin{split}
&\det \left[ \frac{1}{w_it^{-\lambda_i} - w_j}\right]_{i,j = 1}^{N_1} \prod_{i = 1}^{N_1} \det \left[ \frac{1}{\hat{z}_i - \hat{z}_j t^{\mu_j}}\right]_{i,j = 1}^{N_2} \prod_{j = 1}^{N_2}   \frac{(\hat{z}_j w_i^{-1} t^{\lambda_i}; t)_\infty }{(\hat{z}_j w_i^{-1} t^{ \lambda_i + \mu_j}; t)_\infty }\frac{(\hat{z}_j w_i^{-1} t^{\mu_j}; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty } = \\
& \det \left[ \frac{1}{Z_i - W_j}\right]_{i,j = 1}^{N_1 + N_2} \cdot \frac{ \Gamma_t(1 +\log_t \hat{z}_j - \log_t w_i) \Gamma_t (1 + \log_t \hat{w}_j - \log_t z_i)}{\Gamma_t(1 + \log_t \hat{z}_j - \log_tz_i) \Gamma_t(1 + \log_t \hat{w}_j - \log_t w_i )},
\end{split}
\end{equation}
where $Z = (w_1 t^{-\lambda_1}, \dots, w_{N_1} t^{-\lambda_{N_1}}, \hat{z}_1, \dots, \hat{z}_{N_2})$ and $W = (w_1, \dots, w_{N_1}, \hat{z}_1 t^{\mu_1}, \dots, \hat{z}_{N_2}t^{\mu_{N_2}})$ and $\Gamma_t$ denotes the $t$-Gamma function, see \cite[(10.3.3)]{Andrews},
$$\Gamma_t(x) = \frac{(t;t)_\infty}{(t^x;t)_\infty}(1-t)^{1-x}.$$
In deriving (\ref{CrossV1}) we used the identity $(a;t)_\infty = (1-a) (at;t)_\infty$ and the Cauchy determinant formula, see e.g. \cite[1.3]{Prasolov}. In veiw of (\ref{CrossV1}) we can now recognize in (\ref{GenFunExpand3}) both an analogue of the Cauchy determinant in (\ref{GenFunForm}) as well as an analogue of the cross term $CT_{KPZ}$ with Gamma functions being replaced by their discrete analogue.

As we can see, there are many similarities between formulas (\ref{GenFunExpand3}) and (\ref{GenFunForm}), although we do not have a good conceptual understanding of why that is the case. The starting point for these two formulas (i.e. equations (\ref{GenFunExpand}) and (\ref{GenFunExpand2}) ) are clearly analogues of each other; however, the frameworks that get us from those formulas to (\ref{GenFunForm}) and (\ref{GenFunExpand3}) respectively are different. In simple words, we do not have a good way to translate the imaginary time Schr{\"o}dinger equation + Bethe ansatz approach to the one of Macdonald difference operators + contour integral ansatz. The question is further obscured by the different roles played by the joint moments in the two frameworks. Specifically, in the context of the KPZ equation, the moments of the solutions of the SHE are {\em eigenfunctions} for the  imaginary time Schr{\"o}dinger operator $- \partial_T - H_n$, while the $t$-moments in our setup are {\em eigenvalues} for the Macdonald difference operators (the eigenfunctions being the Macdonald symmetric functions). We believe that searching for deeper connections between the two frameworks is warranted, but we will leave this task for the future.


\subsubsection{Contour integral formulas}\label{Section9.1.2}
Equation (\ref{GenFunForm}) can be found as \cite[(3.7)]{ISS13}, and upon a change of variables matches \cite[(25)]{Dot13}. As explained after equation (3.7) in \cite{ISS13}, the authors of \cite{PSpohn} derived an analogous formula (see \cite[(3.8)]{ISS13})  but only under a certain factorization assumption, which corresponds to setting the cross term $CT_{KPZ} = 1$ in (\ref{GenFunForm}). The way the argument proceeds after \cite[(25)]{Dot13} and after \cite[(3.7)]{ISS13} is to rewrite the discrete sums in (\ref{GenFunForm}) as suitable contour integrals and then perform a formal asymptotic analysis demonstrating the emergence of the Fredholm determinant formula for the two-point cdf of the Airy process. We mention here that in rewriting the discrete sums in terms of integrals, \cite{Dot13} ignores the contribution of the infinitely many poles coming from the Gamma functions in the cross term $CT_{KPZ}$. The formal limit shows that the cross term $CT_{KPZ}$ converges to $1$, which {\em a posteriori} gives some credence to why one can ignore the contribution of these poles, and why one can make the factorization assumption in \cite{PSpohn} that corresponds to setting $CT_{KPZ} = 1$ before taking the limit.

In the remainder of this section, we explain what guided our choice of contour deformation in going from (\ref{S1MomentFun2}) to (\ref{S1MomentFun3}) and explain how starting from (\ref{GenFunExpand3}) we can obtain the formula in Theorem \ref{PrelimitT}, which is the one we use in the paper to do asymptotics. In the course of our discussion, we will see that when rewriting (\ref{GenFunExpand3}) as a contour integral we do not encounter the infinite cross-term pole problem from \cite{Dot13}.\\

Starting from (\ref{S1MomentFun2}) there are a few natural ways one can attempt to deform the contours using the contour integral ansatz. The first is to simply deform all $L+R$ contours down to the same contour $C^2_{0,R}$ (the innermost circle). If one does this, the residue subspaces become labeled by partitions $\nu \vdash (L +R)$. The structure that is lost in performing this ansatz is that the resulting integral is no longer (visibly) analytic in the parts $(\nu_1, \dots, \nu_{\ell(\nu)})$. Since eventually, we will want to rewrite any discrete sum over partitions as contour integrals, it is desirable for the resulting expression after the contour integral ansatz to be an analytic function of the parts of these partitions. As one can see from (\ref{S1MomentFun3}) this is clearly the case for the deformation we performed. Another possibility is to deform all contours $C^1_{0,i}$ down to $C^{1}_{0,L}$ and all contours $C^{2}_{0,i}$ down to $C^2_{0,R}$. After doing this contour integral ansatz, one obtains an expression that is very similar to (\ref{S1MomentFun3}). The problem that arises in this case is that the requirement for $a$ to be inside $C^2_{0,R}$ and $a^{-1}$ outside $C^{1}_{0,L}$ makes the analogue of (\ref{S1MomentFun3}) valid only for small enough (depending on $R+L$) values of $a$. In particular, we cannot find contours that work for all the terms in (\ref{GenFunExpand2}) simultaneously unless $a = 0$, which trivializes the model. There are a few other possible ways of deforming the contours, but they all run into one of the above two problems and the only working solution we have is through equation (\ref{S1MomentFun3}).

We finally, discuss how to rewrite the sums over $\mu_i$ and $\lambda_i$ in (\ref{GenFunExpand3}) as contour integrals. The keen reader may have noticed that there is an extra $t^{-1}$ next to $u_1$ in the second line of (\ref{GenFunExpand3}) compared to the third line. This extra term is not a typo and after re-expressing the discrete sums in (\ref{GenFunExpand3}) with integrals below, we will obtain a symmetric expression. The key identity we use is as follows:
$$\sum_{n_1, \dots, n_k = 1}^\infty f(t^{n_1}, \dots, t^{n_k}) y_1^{n_1} \cdots y_k^{n_k} = \frac{1}{(2 \pi \iota)^k} \int_{\ell_\delta^k} \prod_{i = 1}^k \frac{\pi [-y_i ]^{s_i} }{\sin (-\pi s_i)} f(t^{s_1}, \dots, t^{s_k}) ds_1 \cdots ds_k,$$
where $\ell_\delta$ is a vertical contour passing through $\delta \in (0,1)$, with the orientation of increasing imaginary part. The latter equation essentially follows from $Res_{s = n} \frac{\pi}{\sin(-\pi s)} = (-1)^{n+1}$. Fixing a small $\delta > 0$ and applying the above formula to (\ref{GenFunExpand3}) we obtain
\begin{equation}\label{GenFunExpand4}
\begin{split}
&\mathbb{E}\left[ \frac{1}{(u_1t^{h(n_1 + 1 , M) - M};t)_{\infty}} \frac{1}{(u_2t^{h(n_2 + 1 , M) - M};t)_{\infty}} \right]=\sum_{N_1 = 0}^\infty \sum_{N_2 = 0}^{\infty} \frac{1}{N_1! N_2!}\int_{C_{1 + \varepsilon}^{N_1}} \int_{C_{1 - \varepsilon}^{N_2}} \int_{\ell^{N_1}_\delta}\int_{\ell^{N_2}_\delta} \\
&  \det \left[ \frac{1}{w_it^{-s_i} - w_j}\right]_{i,j = 1}^{N_1} \prod_{i = 1}^{N_1} \left( \frac{1 - aw_it^{- s_i}}{1 - aw_i } \right)^M  \left(\frac{1 - a w_i^{-1} }{1 - a w_i^{-1}  t^{s_i}}\right)^{n_1}\frac{\pi [-u_1t^{-1}]^{s_i}}{\sin ( -\pi s_i)} \\
&  \det \left[ \frac{1}{\hat{z}_i - \hat{z}_j t^{v_j}}\right]_{i,j = 1}^{N_2} \prod_{i = 1}^{N_2} \left( \frac{1 - a\hat{z}_i }{1 - a \hat{z}_it^{v_i}} \right)^M    \left( \frac{1 - a\hat{z}^{-1}_it^{-v_i}}{1 - a\hat{z}^{-1}_i } \right)^{n_2} \frac{ [-u_2]^{v_i}}{\sin(-\pi v_i)}\\
& \prod_{i = 1}^{N_1}  \prod_{j = 1}^{N_2}   \frac{(\hat{z}_j w_i^{-1} t^{s_i}; t)_\infty }{(\hat{z}_j w_i^{-1} t^{ s_i + v_j}; t)_\infty }\frac{(\hat{z}_j w_i^{-1} t^{v_j}; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty } \prod_{i = 1}^{N_2} \frac{d v_i}{2\pi \iota}\prod_{i = 1}^{N_1} \frac{d s_i}{2\pi \iota} \prod_{i = 1}^{N_2} \frac{d \hat{z}_i}{2\pi \iota}\prod_{i = 1}^{N_1} \frac{d w_i}{2\pi \iota}.
\end{split}
\end{equation}
We point out that the integrand in (\ref{GenFunExpand4}) has no poles in the right half-plane in either of the $s_i$ or $v_i$ variables. In particular, unlike \cite{Dot13}, our contour representation does not include any poles from the cross term. Finally, we may perform the change of variables $z_i = w_i t^{-s_i}$ for $i = 1, \dots, N_1$ and $\hat{w}_i = \hat{z}_i t^{v_i}$ for $i = 1, \dots, N_2$ to rewrite (\ref{GenFunExpand4}) as 
\begin{equation}\label{GenFunExpand5}
\begin{split}
&\mathbb{E}\left[ \frac{1}{(u_1t^{h(n_1 + 1 , M) - M};t)_{\infty}} \frac{1}{(u_2t^{h(n_2 + 1 , M) - M};t)_{\infty}} \right]=\sum_{N_1 = 0}^\infty \sum_{N_2 = 0}^{\infty} \frac{1}{N_1! N_2!}\int_{C_{1 + \varepsilon}^{N_1}}\int_{C_{1 - \varepsilon}^{N_2}}  \\
& \int_{C_{(1+\varepsilon)t^{-\delta}}^{N_1}}\int_{C_{(1-\varepsilon)t^{\delta}}^{N_2}}  \det \left[ \frac{1}{z_i - w_j}\right]_{i,j = 1}^{N_1} \prod_{i = 1}^{N_1} \left( \frac{1 - az_i}{1 - aw_i } \right)^M  \left(\frac{1 - a w_i^{-1} }{1 - a z_i^{-1}  }\right)^{n_1}\frac{S(w_i, z_i; u_1, t)}{[-\log t] w_i} \\
&  \det \left[ \frac{1}{\hat{z}_i - \hat{w}_j }\right]_{i,j = 1}^{N_2} \prod_{i = 1}^{N_2} \left( \frac{1 - a\hat{z}_i }{1 - a \hat{w}_i} \right)^M    \left( \frac{1 - a\hat{w}^{-1}_i}{1 - a\hat{z}^{-1}_i } \right)^{n_2}\frac{S(\hat{w}_i, \hat{z}_i; u_2, t)}{[-\log t] \hat{w}_i}\\
& \prod_{i = 1}^{N_1}  \prod_{j = 1}^{N_2}   \frac{(\hat{z}_j z_i^{-1} ; t)_\infty }{(\hat{w}_j z_i^{-1} ; t)_\infty }\frac{(\hat{w}_j w_i^{-1} ; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty } \prod_{i = 1}^{N_2} \frac{d \hat{w}_i}{2\pi \iota}\prod_{i = 1}^{N_1} \frac{d z_i}{2\pi \iota} \prod_{i = 1}^{N_2} \frac{d \hat{z}_i}{2\pi \iota}\prod_{i = 1}^{N_1} \frac{d w_i}{2\pi \iota},
\end{split}
\end{equation}
where $C_r$ is a positively oriented circle of radius $r$, centered at the origin, and $S(w,z;u,t)$ is the function from Definition \ref{DefFunS}:
\begin{equation}\label{S1SpiralDef}
S(w, z; u,t) =   \sum_{m \in \mathbb{Z}} \frac{\pi \cdot [ - u ]^{[\log w - \log z] [\log t]^{-1} -  2m \pi  \iota [\log t ]^{-1}}}{\sin(-\pi [[\log w - \log z] [\log t]^{-1} -  2m \pi  \iota [\log t ]^{-1}])}.
\end{equation}
Notice that (\ref{GenFunExpand5}) is now symmetric in the $u_1, u_2$ variables (compared to (\ref{GenFunExpand4})) and matches equation (\ref{PrelimitEq}) from Theorem \ref{PrelimitT} for the case when $x_i = a = y_j$. 


\subsection{Convergence and asymptotics}\label{Section9.2}
In Section \ref{Section9.2.1} we explain how our prelimit formulas compare to the ones for the log-gamma polymer in \cite{NZ}. In Section \ref{Section9.2.2} we address the questions of convergence of our formulas, and give an overview of our asymptotic analysis.


\subsubsection{Log-gamma polymer}\label{Section9.2.1}\label{Section9.2.1} In this section we summarize the approach of studying the two-point asymptotics of the log-gamma polymer model from \cite{NZ} and compare it to the present paper. The log-gamma polymer depends on two sequences of real parameters $\{ \alpha_i \}_{i \geq 1}$ and $\{\hat{\alpha}_i\}_{i \geq 1}$ such that $\alpha_i + \hat{\alpha_j} > 0$ for all $i, j \in \mathbb{N}$. One takes a countable collection of independent random variables $w_{i,j}$ indexed by $(i,j) \in \mathbb{Z}^2_{\geq 1}$ such that $w_{i,j}$ has the {\em inverse-Gamma distribution with parameter} $\theta_{i,j} = \alpha_i + \hat{\alpha}_i$, meaning that 
$$\mathbb{P}(w_{i,j} \in dx) ={\bf 1}\{ x > 0 \} \cdot \frac{1}{\Gamma(\theta_{i,j})} x^{-\theta_{i,j} -1} \exp (-1/x) dx.$$
From this data one defines the (random) {\em log-gamma partition function} as
\begin{equation}\label{S1LGPP}
Z_{m,n}: = \sum_{\pi \in \Pi_{m,n}} \prod_{(i,j) \in \pi} w_{i,j},
\end{equation}
where $(m,n) \in \mathbb{Z}^2_{\geq 1}$, the sum is over all up-right paths from $(1,1)$ to $(m,n)$ and the product is over the vertices that belong to $\pi$. We remark that in \cite{NZ} the convention is for the numbers to increase in the downward direction so that the paths are down-right. This should cause no confusion.

\cite[Conjecture 4.1]{NZ} states that if $(m_1, n_1) = (N - t_1 N^{2/3}, N + t_1 N^{2/3})$ and $(m_2, n_2) = (N + t_2N^{2/3}, N - t_2N^{2/3})$ with $t_1, t_2 > 0$ then the vectors $(\log Z_{m_1, n_1}, \log Z_{m_2, n_2})$ under suitable shifts and scales should weakly converge to the two-point distribution of the Airy process. The starting point of the analysis is the following formula for the joint Laplace transform for $Z_{m_1, n_1}$ and $Z_{m_2, n_2}$
\begin{equation}\label{S1LGJL}
\begin{split}
&\mathbb{E}\hspace{-1mm} \left[ e^{-u_1 Z_{m_1, n_1} - u_2 Z_{m_2, n_2}}\right] \hspace{-1mm}= \hspace{-1mm}\int_{\ell_{\delta}^{m_1}}\hspace{-1mm}\int_{\ell_{\delta+ \gamma}^{n_2}}  s_{m_1}(\vec{\lambda}) s_{n_2}(\vec{\mu})  F^{LG}(\vec{\lambda},\vec{\mu};u_1, u_2)  \prod_{i = 1}^{m_1} \prod_{j = 1}^{n_2} \frac{\Gamma(\lambda_i + \mu_j)}{\Gamma(\alpha_i + \hat{\alpha}_j)}d\vec{\mu}  d\vec{\lambda},
\end{split}
\end{equation}
which is obtained via the geometric Robinson-Schensted-Knuth correspondence. In (\ref{S1LGJL}) the contour $\ell_x$ is a vertically oriented line passing through $x$ in the complex plane, and $s_n(\vec{\lambda})$ denotes the {\em Sklyanin measure}
$$s_n(\lambda) = \frac{1}{(2\pi \iota)^n n!} \prod_{i \neq j } \Gamma(\lambda_i - \lambda_j)^{-1}.$$
The precise formula (i.e. the definition of $\delta, \gamma$ and $F^{LG}$ in (\ref{S1LGJL})) can be found as \cite[(2.21)]{NZ}. Equation (\ref{S1LGJL}) is not suitable for asymptotic analysis since the number of contours goes to infinity, and in order to re-express it \cite{NZ} uses \cite[Theorem 2]{BCR}, which allows one to rewrite the $d\vec{\mu}$ and $d\vec{\lambda}$ integrals in (\ref{S1LGJL}) as Fredholm determinants. 

By applying \cite[Theorem 2]{BCR} twice to (\ref{S1LGJL}), expanding the corresponding Fredholm determinants and formally rearranging the resulting contours (see \cite[(4.13)]{NZ} and the discussion surrounding that equation for more information on why the rearrangement is formal) one obtains
\begin{equation}\label{S1LGJLV2}
\begin{split}
&\mathbb{E}\left[ e^{-u_1 Z_{m_1, n_1} - u_2 Z_{m_2, n_2}}\right] = \sum_{N_1 = 0}^{m_1} \sum_{N_2 = 0}^{n_2} \frac{1}{N_1! N_2!} \frac{1}{ (2\pi \iota)^{2N_1 + 2N_2}} \int_{\ell_{\delta_1}^{N_1}} \int_{C_{\delta_1}^{N_1}}\int_{\ell_{\delta_1}^{N_2}} \int_{C_{\delta_1}^{N_2}}   \\
&D( \vec{w}, \vec{z})  {G}^{LG}_{n_1,m_1}(\vec{w},\vec{z} , u_1) \cdot D( \vec{\hat{w}}, \vec{\hat{z}})  {G}^{LG}_{n_2, m_2}( \vec{\hat{w}},\vec{\hat{z}}, u_2)  \cdot CT_{LG}(\vec{w}, \vec{z}; \vec{\hat{w}}, \vec{\hat{z}}) d\vec{\hat{w}} d\vec{\hat{z}}d\vec{w} d\vec{z},
\end{split}
\end{equation}
where $C_{\delta}$ is a positively oriented, zero-centered circle of radius $\delta$, $D(\vec{a}, \vec{b}) = \det \left[\frac{1}{b_i - a_j} \right]_{i,j = 1}^{n}$ is the Cauchy determinant and the cross term $CT_{LG}$ is given by
\begin{equation}\label{S1LGCT}
\begin{split}
&CT_{LG}(\vec{w}, \vec{z}; \vec{\hat{w}}, \vec{\hat{z}}) = \prod_{i = 1}^{N_1} \prod_{j = 1}^{N_2} \frac{\Gamma(\gamma - w_i - \hat{w}_j) \Gamma(\gamma -z_i - \hat{z}_j)}{\Gamma(\gamma - w_i - \hat{z}_j) \Gamma(\gamma -z_i - \hat{w}_j)}
\end{split}
\end{equation}
In equation (\ref{S1LGJLV2}) all the $\alpha$ parameters of the model are set to $0$, and all the $\hat{\alpha}$ parameters are set to $\gamma > 0$, and $\delta \in (0, \gamma/2)$, while $\delta_1 \in (0, \min (\delta, 1- \delta))$. We will not write down the formulas for ${G}^{LG}_{n_i,m_i}$ but refer the interested reader to \cite[Proposition 4.5]{NZ}. \\

The analogue of (\ref{S1LGJLV2}) we derive in this paper formulated in terms of the stochastic six-vertex model is as follows
\begin{equation}\label{GenFunExpand6}
\begin{split}
&\mathbb{E}\left[ \frac{1}{(u_1t^{h(n_1 + 1 , M) - M};t)_{\infty}} \frac{1}{(u_2t^{h(n_2 + 1 , M) - M};t)_{\infty}} \right]=\sum_{N_1 = 0}^\infty \sum_{N_2 = 0}^{\infty}  I_{M}(N_1,N_2), \mbox{ where }\\
&I_M(N_1, N_2) = \frac{1}{N_1! N_2!}\int_{\gamma_1^{N_1}} \int_{\gamma_2^{N_1}}  \int_{\gamma_3^{N_2}}\int_{\gamma_4^{N_2}} D(\vec{w}, \vec{z})G(\vec{w}, \vec{z}, n_1, u_1)      \\
&  D(\vec{\hat{w}}, \vec{\hat{z}})G(\vec{\hat{w}}, \vec{\hat{z}},n_2, u_2) \cdot  CT(\vec{w}, \vec{z}; \vec{\hat{w}}, \vec{\hat{z}}) \prod_{i = 1}^{N_2}\frac{d\hat{w}_i}{2\pi \iota}\prod_{i = 1}^{N_2}\frac{d\hat{z}_i}{2\pi \iota}\prod_{i = 1}^{N_1}\frac{dw_i}{2\pi \iota}\prod_{i = 1}^{N_1}\frac{dz_i}{2\pi \iota}.
\end{split}
\end{equation}
In (\ref{GenFunExpand6}) the contours $\gamma_i$ are positively oriented zero-centered circles of radius $r_i \in (a, a^{-1})$ such that $r_1 > r_2 > r_3 > r_4 > tr_1$, where we recall that $t = b_1/b_2$ and $a^2 = \frac{1 - b_1}{1-b_2}$. In addition, we have
\begin{equation}\label{S1MPLTGv2}
\begin{split}
& D(\vec{w}, \vec{z}) =   \det \left[ \frac{1}{z_i - w_j}\right]_{i,j = 1}^{N_1} \hspace{5mm} CT(\vec{w}, \vec{z}; \vec{\hat{w}}, \vec{\hat{z}}) =  \prod_{i = 1}^{N_1}\prod_{j = 1}^{N_2}  \frac{(\hat{z}_j z_i^{-1}; t)_\infty }{(\hat{w}_j z_i^{-1}; t)_\infty }\frac{(\hat{w}_j w_i^{-1} ; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty },  \\
&G(\vec{w}, \vec{z}, n, u) =  \prod_{i = 1}^{N_1} \left( \frac{1 - az_i}{1 - aw_i} \right)^M \cdot \left( \frac{1 - a/ w_i}{1 - a/ z_i}\right)^n \cdot \frac{ S(w_i, z_i; u,t )}{ -\log t \cdot w_i} , \\
\end{split}
\end{equation}
where $S(w,z;u,t)$ is the function from Definition \ref{DefFunS} (alternatively see equation (\ref{S1SpiralDef}). If $N_1 = N_2 = 0$ then we adopt the convention $I_M(0,0) = 1$. Equation (\ref{GenFunExpand6}) is formulated for the ascending Hall-Littlewood process in Theorem \ref{PrelimitT}, and one deduces the analogous result from the distributional equality of the latter and the stochastic six-vertex model, see the end of Section \ref{Section2.2}.


Comparing (\ref{S1LGJLV2}) with (\ref{GenFunExpand6}) we observe striking structural similarities. On the left side of (\ref{GenFunExpand6})  we have a joint $t$-Laplace transform, which is a discrete analogue of the joint Laplace transform in (\ref{S1LGJLV2}). On the right side of both (\ref{S1LGJLV2}) and (\ref{GenFunExpand6})  we have an analogous double sum of $(2N_1 + 2N_2)$-fold contour integrals, where the integrand is a product of similar three types of terms. The third term in both formulas is a cross term, which in the log-gamma case is a double product of Gamma functions, while in our case is a double product of $t$-Gamma functions, cf. (\ref{CrossV1}). We believe that the similarities between (\ref{S1LGJLV2}) and (\ref{GenFunExpand6}) are not coincidental but come from the known connections between the Macdonald processes and the log-gamma polymer, cf. \cite[Section 4]{BCFV}. Specifically, we expect that one can carry out a similar program to the one in the present paper and derive joint $q$-Laplace transform formulas for the $q$-Whittaker process (a certain dual to the Hall-Littlewood process considered in our paper). Since the $q$-Whittaker process converges to the log-gamma polymer model, as shown in \cite[Section 4]{BCFV}, we expect such joint $q$-Laplace transform formulas to asymptotically give the formulas in \cite{NZ}. At this time this approach is purely conjectural, but if correct it would explain the strong similarities we see between our formulas and those of \cite{NZ}.


\subsubsection{Handling the cross term}\label{Section9.2.2}
We first discuss the convergence issues in equation (\ref{GenFunExpand6}). There are four types of terms that influence the growth of $I_M(N_1,N_2)$: (1) the functions $G$, (2) the cross term $CT$, (3) the Cauchy determinants $D$ and (4) the factorials $N_1!N_2!$ in the front. Unlike (\ref{S1LGJLV2}) all of our contours are compact and so obtaining point-wise estimates for the integrand essentially suffices for controlling $I_M(N_1,N_2)$. Out of the four terms, the functions $G$ are benign, because they grow as $e^{C(N_1 + N_2)}$ and are thus controlled by the factorials that decay like $e^{-N_1 \log N_1 - N_2 \log N_2}$. In the past, the known way to control the Cauchy determinant is via Hadamard's inequality, which states that 
\begin{equation}\label{S1Had}
|\det A| \leq \prod_{i = 1}^N \|v_i\|
\end{equation}
for an $N \times N$ matrix $A$ with column vectors $v_1, \dots ,v_N$. Consequently, if we ignore the cross term in (\ref{GenFunExpand6}) all other terms can be bounded by
$$\exp\left( O(N_1 + N_2) - \frac{N_1 \log N_1}{2} - \frac{N_2 \log N_2}{2} \right),$$
which would be summable over $(N_1, N_2) \in \mathbb{Z}_{\geq 0}^2.$ The problem we run into is that the cross term $CT$ is pointwise of order $e^{cN_1N_2}$ and so the factorials are no longer enough to control it. 

It is worth pointing out that the prelimit formula (\ref{S1LGJLV2}) from \cite{NZ} does not encounter this cross term problem, because the sum is {\em finite}. The reason \cite{NZ} obtained a finite sum is a consequence of the fact that in applying \cite[Theorem 2]{BCR} the kernels involved are finite rank, and the resulting Fredholm determinant epansion series terminate after finitely many terms. The $G^{LG}$ functions in (\ref{S1LGJLV2}) all decay sufficiently fast near infinity and so even though the contours in that formula are infinite each summand is easily seen to be finite, and so the whole sum is finite. As one takes the $N \rightarrow \infty$ limit in (\ref{S1LGJLV2})  the cross terms will start being problematic, but at least at a finite $N$ level they do not cause problems unlike our setup.

The way we deal with $CT$ in our formula is by utilizing some {\em hidden built-in} decay in the Cauchy determinants $D$, which to our knowledge has not been previously recognized. In particular, we have the following estimate for any $r, R \in (0,\infty)$ with $R > r$
\begin{equation}\label{S1CDetGood}
\left|\det \left[ \frac{1}{z_i - w_j}\right]_{i,j = 1}^N \right| \leq R^{-N} \cdot \frac{N^N \cdot (r/R)^{\binom{N}{2}}}{(1-r/R)^{N^2}},
\end{equation}
where $z_i, w_i \in \mathbb{C}$ are such that $|z_i| = R$ and $|w_i|  \leq r$ for $i = 1, \dots, N$. This results appears as Lemma \ref{DetBounds} in the main text.

The way we use (\ref{S1CDetGood}) to control the summands in (\ref{GenFunExpand6}) is by deforming the $\gamma_i$ contours so that $r_2/r_1$ and $r_4/r_3$ are both very small and then one obtains the following bound on the Cauchy determinants 
\begin{equation}\label{S1Had2}
|D(\vec{w}, \vec{z}) D(\vec{\hat{w}}, \vec{\hat{z}})| \leq  \exp( - c (N_1^2  +N_2^2)),
\end{equation}
where the constant $c$ increases as the ratios $r_2/r_1$ and $r_4/r_3$ become smaller. The restrictions $r_4 > tr_1$ and $r_i \in (a, a^{-1})$ prevent us from freely deforming the contours $\gamma_i$ so that these ratios are arbitrarily small, but at least if $t$ and $a$ are small enough, this is possible and one is able to control the cross term $CT$ using the Cauchy determinants $D(\vec{w}, \vec{z})$, $D(\vec{\hat{w}}, \vec{\hat{z}})$. This is one of the sources of the parameter restriction in Theorem \ref{thmMain}, and we can only prove (\ref{GenFunExpand6}) in Theorem \ref{PrelimitT} for a special range of parameters. See Remark \ref{RemarkRest} for more details.\\

We conclude this section with a discussion about the asymptotic analysis of equations (\ref{S1LGJLV2}) and (\ref{GenFunExpand6}). In both equations  (\ref{S1LGJLV2}) and  (\ref{GenFunExpand6}) it is relatively easy to show that the right sides term-wise converge to a suitable Fredholm determinant expansion for the Airy process (as $N \rightarrow \infty$ in (\ref{S1LGJLV2}) and $M \rightarrow \infty$ in (\ref{GenFunExpand6})). The arguments in both cases involve deforming contours to descent ones for the $G$ functions and applying a careful steepest descent argument. In the present paper, the term-wise convergence is stated as Proposition \ref{PropTermConv} and proved in Section \ref{Section5}.  The identification of the resulting sum with the Fredholm determinant expansion for the Airy process is formulated as Proposition \ref{PropTermLimit} and proved in Section \ref{Section7.2}.


The essential ingredient missing and making the proof in \cite{NZ} incomplete is a uniform in $N$ estimate on the growth of the terms (in terms of $N_1$ and $N_2$) in the series (\ref{S1LGJLV2}) that would allow one to exchange the order of the sum and the limit $N \rightarrow \infty$. The difficulty in obtaining such growth estimates comes from the presence of the cross term $CT_{LG}$, whose behavior as $N_1, N_2$ become large is very complicated. Part of the progress made in our paper is the ability to control the cross terms $CT$ in our formulas (at least for some range of parameters) and obtain the necessary bounds that would allow one to exchange the sums and the $M \rightarrow \infty$ limit in (\ref{GenFunExpand6}). The precise bounds we can obtain are given in Proposition \ref{PropTermBound} and its proof can be found in Section \ref{Section6}. 

While the proof of Proposition \ref{PropTermBound} is fairly technical, we try to give a rough account of the ideas contained in it below. The way the proof of Proposition \ref{PropTermBound} goes is by splitting the problem into three cases: (1) $\Delta \gtrsim M$, (2) $M^{1/2} \lesssim \Delta \lesssim M$ and (3) $\Delta \lesssim M^{1/2}$, where $\Delta = N_1 + N_2$. The goal is then to deform the contours $\gamma_1, \dots, \gamma_4$ in a suitable way so that the different terms ($D$, $G$ and $CT$) can be balanced in a favorable way so that a rapid enough decay in $\Delta$ is achieved. The Cauchy determinants $D$ are always helpful (i.e. they are sources of decay), the cross term $CT$ is always harmful (it is a source of growth) and the $G$ terms can be helpful if the contours $\gamma_i$ are descent contours or harmful if they are not. 

In the case $\Delta \gtrsim  M$ the idea is simply to spread apart the $\gamma_i$ so that the ratios of radii $r_2/r_1$ and $r_4/r_3$ are very small. By doing this we obtain
$$ \left| D(\vec{w}, \vec{z}) D(\vec{\hat{w}}, \vec{\hat{z}}) \right| \leq \exp( -c_1 \Delta^2), \left|G(\vec{w}, \vec{z}, n_1, u_1) G(\vec{\hat{w}}, \vec{\hat{z}},n_2, u_2) \right| \leq \exp ( c_2 \Delta M), \left| CT \right|\leq  \exp ( c_3 \Delta^2).$$
Notice that by spreading apart the contours $\gamma_i$, we move them away from the descent contours for the $G$ functions and so they grow; however, if $a$ is sufficiently small and $\Delta \gtrsim M$ the Cauchy determinants are decaying fast enough to offset both the $G$ functions and the cross term $CT$, and one obtains 
$$|I_M(N_1,N_2)| \leq \exp \left( - \epsilon \Delta^2 \right).$$

In the case $\Delta \lesssim M^{1/2}$ the idea is simply to deform $\gamma_i$ to descent contours. In this case, the $G$ functions decay fast enough to offset the cross term, and one can bound the Cauchy determinants by Hadamard's inequality (\ref{S1Had}) and use the factorials $N_1! N_2!$ in $I_M(N_1, N_2)$ to prove
$$|I_M(N_1, N_2)| \leq \exp \left( O(\Delta) - \frac{\Delta \log \Delta}{4}\right).$$

The case $M^{1/2} \lesssim \Delta \lesssim M$ is the more involved. In this case we spread the $\gamma_i$ contours a little bit, so that {\em some} decay is extracted from the Cauchy determinants $D$, but at the same time the growth of the $G$-functions, which are big away from their descent contours, is still manageable. In this case, we need to Taylor expand $\log CT$ upto second order around the critical point of $\log G$, and use the $G$ functions to control $CT$ when enough of the $z$ and $w$ variables are away from the critical point. If many of these variables are close to the critical point, then that implies that they are close to each other and one obtains {\em improved} decay estimates for the Cauchy determinants $D$. This extra bit of decay can be traced to the presence of the Vandermonde determinants in the numerator of the Cauchy determinant formula, which makes the $D$ terms even smaller if many of the $z$, $w$ variables start being close to each other. Ultimately, there is a delicate balance of all three types of terms $D, G$ and $CT$ that is favorably resolved leading to the estimate
$$|I_M(N_1,N_2)| \leq \exp \left( - \epsilon \Delta^2 \right),$$
provided the $t$ and $a$ parameters are sufficiently small. 

The above exposition, aimed to illustrate on a very basic level how the proof of Proposition \ref{PropTermBound} is structured. The interested reader is referred to Section \ref{Section6} for the technical details.
{"config": "arxiv", "file": "2006.15934/Section9.tex"}