\begin{equation}\label{eq:HN4}
\begin{cases}
\lambda g(z)+\mathcal{H}(\phi(z)+\lambda  \overline{g(z)})=\lambda \psi(z), \\
\overline{\eta(z)}+\mathcal{H}(\eta(z))=\overline{\phi(z)}+\lambda \psi(z).
\end{cases}
\end{equation}

We now solve for $\eta$ and $g$ in terms of $\psi$ and $\phi$. Using the Lewy transform

\begin{align*}
L_\sigma(z)=\overline{z+\nabla \sigma(z)}
\end{align*}

we have by <cit.>

\begin{align*}
L_\sigma^{-1}(w)=\frac{1}{2}\bigg(\overline{w}+\mathcal{H}(w)\bigg)\in N^\circ, \quad \overline{w}\in \text{dom}\, (\mathcal{H}).
\end{align*}

and so the equation $L_\sigma^{-1}(w)=\zeta$ has solution $w=L_{\sigma}(\zeta)$, which gives

\begin{align*}
\eta(z)&=L_\sigma\bigg(\frac{1}{2}(\overline{\phi(z)}+\lambda \psi(z))\bigg)\\
g(z)&=\frac{\overline{\eta(z)}-\overline{\phi(z)}}{\lambda}=\frac{1}{\lambda}(\overline{L_{\sigma}((\overline{\phi(z)}+\lambda \psi(z)))/2}-\overline{\phi(z)}).
\end{align*}

It is tempting to use the parametrisation of $\eta$ and $g$ in terms of the pair of holomorphic functions $(\phi,\psi)$ to deduce regularity of the height function $h$ and the boundary of the liquid domain $\dv \LL$. Unfortunately, it seems very difficult to determine which holomorphic functions $\psi$ and $\phi$ that gives rise to homeomorphic $g$ and proper $\eta$. Furthermore it seems difficult to connect to the case when $\lambda=0$ and $g$ is teleomorphic. What we want is to view $g$ in the case $\lambda\neq 0$ as a natural deformation of a teleomorphic $g$ as in [3]. Luckily it seems that a certain combination of $\phi$ and $\psi$ can be determined. We now define

\begin{align*}
\omega(z):=\frac{1}{2}(\overline{\phi(z)}+\lambda \psi(z)).
\end{align*}

Then $\Delta \omega=0$ so $\omega$ is a harmonic map. Furthermore, since $L_\sigma^{-1}$ is a homeomorphism and $u$ is proper and $\omega=L_\sigma(u)$ it follows that $\omega: \Di \to N$ is a proper harmonic map. Hence $\omega$ is $k:1$ for some $k$ and an orientation preserving map. Furthermore, since $\mathcal{H}'$ is a proper holomorphic map it follows that $\mathcal{H}'(\eta)$ is proper as well. The case of harmonic homeomorphisms $\omega:\Di \to N$, where $N$ is a convex compact polygon in the plane is well-understod, see <cit.>, although the authors are unaware of a complete classification of $k:1$ harmonic mappings $\omega: \Di \to N$. Assuming one has a good understanding of this class of harmonic maps, instead of studying (<ref>) one could study the linear equation

\begin{align*}
\end{align*}

where $\omega$ is now any proper harmonic map $\omega: \Di \to N$ and determine the boundary regularity for $g$. A key idea in [3] for studying the boundary regularity in the case when $\lambda=0$ is to consider the map

\begin{align*}
\gamma(z)=g(z)+\mathcal{H}'(\eta(z))\overline{g(z)},
\end{align*}

which is holomorphic in this case. This is no longer the case when $\lambda\neq 0$. In this case one can show that $\gamma$ solves the equation

\begin{align*}
\gamma_{\dbb}(z)&=\alpha(z)\gamma(z)+\beta(z)\overline{\gamma(z)},
\end{align*}

\begin{align*}
\alpha(z)&=\frac{\mathcal{H}''(\eta(z))\eta_{\dbb}(z)}{1-\vert \mathcal{H}'(\eta(z))\vert^2}\overline{\mathcal{H}'(\eta(z))},\\
\beta(z)&=\frac{\mathcal{H}''(\eta(z))\eta_{\dbb}(z)}{1-\vert \mathcal{H}'(\eta(z))\vert^2}.
\end{align*}

Thus $\gamma$ is a generalized analytic function in the sense of Vekua. Then Vekua's factorization theorem implies that

\begin{align*}
\gamma(z)=\Phi(z)e^{\kappa(z)}.
\end{align*}

for some holomorphic $\Phi$ and Hölder continuous $\kappa$. Solving for $g$ gives

\begin{align*}
g(z)&=\frac{1}{1-\vert \mathcal{H}'\eta(z))\vert^2}(\gamma(z)-\mathcal{H}'(\eta(z))\overline{\gamma(z)}).
\end{align*}

To understand the geometry of the frozen boundary one must then study the singular limit

\begin{align*}
\lim_{z\in \Di\to z_0\in \dv \Di}\frac{1}{1-\vert \mathcal{H}'(\eta(z))\vert^2}(\gamma(z)-\mathcal{H}'(\eta(z))\overline{\gamma(z)}).
\end{align*}

If one can fully understand the case when the weight function $w$ is constant, it is conceivable that the case for general non-singular $w$ can be reduced to that for a constant $w$ by approximation and localization arguments. The case with singular weights $w$ however would require completely different methods.

§.§.§ Complex structure and complex Burger's equation

We next show that the complex structure  $\mathcal H'(f)$  is directly related to the complex slope which solves an inhomogeneous complex Burgers equation (as expected). To this end, we first  define the complex valued functions $\Upsilon,\Om$ on $\LL$ through

\begin{align}\label{eq:UO}
\log \Upsilon&=\pi \sigma_1(\nabla h)-i\pi h_y\\
\log \Om&= \pi \sigma_2(\nabla h)+i\pi h_x.\label{eq:UOp}
\end{align}

Note that our $\sigma$ is normalized differently than in [34]. As in [34] one verifies that

\begin{align} \label{eq:UO2}
\Upsilon+\Om\equiv 1,
\end{align}

and that $(\Upsilon,\Om)$ solves the complex Burger equation

\begin{align}\label{eq:CB}
\frac{\Upsilon_x(z)}{\Upsilon(z)}+\frac{\Om_y(z)}{\Om(z)}=w(z),\quad z\in \LL
\end{align}

whenever $h$ solves the Euler-Lagrange equation (<ref>) in $\LL$.

From (<ref>), (<ref>) and (<ref>) we find that the limiting density for the lozenges can be read off from  the triangle  formed by $\Omega$, $0$ and $1$. We  relate $\Upsilon$ and $\Om$ to the complex structure $\mathcal{H}'(f)$.

Let $(\Upsilon,\Om)$ be given by (<ref>) and (<ref>). Then

\begin{align}\label{eq:ComplexS}
\mathcal{H}'(f)=-\frac{i+(1-i)\Upsilon}{i-(1+i)\Upsilon}=-\frac{1+(i-1)\Om}{-1+(1+i)\Om}.
\end{align}

Furthermore, $\Upsilon$ and $\Om$ are related to $f$ in Theorem <ref> through

\begin{align*}
\Upsilon&=e^{\frac{\pi}{2}(f-\mathcal{H}(f))},\\
\Om&=e^{i\frac{\pi}{2}(f+\mathcal{H}(f))}.
\end{align*}

Identifying $\nabla \sigma(z)$ with $\sigma_x+i\sigma_y$, we have

\begin{align*}
\pi \sigma_x(z)-i\pi y&=\frac{\pi}{2}(\nabla \sigma(z)+\overline{\nabla \sigma(z)})-\frac{\pi}{2}(z-\overline{z})\\
&=-\frac{\pi}{2}(z-\nabla \sigma(z))+\frac{\pi}{2}\overline{(z+\nabla \sigma(z))}\\
\pi \sigma_y(z)+i\pi x&=i\frac{\pi}{2}(\overline{\nabla \sigma(z)}-\nabla \sigma(z))+i\frac{\pi}{2}(z+\overline{z})\\
&=i\frac{\pi}{2}(z-\nabla \sigma(z))+i\frac{\pi}{2}\overline{(z+\nabla \sigma(z))}
\end{align*}

We recall that

\begin{align*}
\mathcal{H}(z)=(I-\nabla \sigma)(I+\nabla \sigma)^{-1}(\overline{z}),
\end{align*}

and that by Theorem <ref>,

\begin{align*}
\nabla h(z)=U(\overline{f(z)}),
\end{align*}

where $U=((I+\nabla \sigma)^{-1})$ is the inverse of the Lewy transform $L(z)=z+\nabla \sigma(z)$. Set $z=U(\overline{w})$. Then

\begin{align*}
\log \gamma(w)&:=-\frac{\pi}{2}((U(\overline{w}))-\nabla \sigma(U(\overline{w})))+\frac{\pi}{2}\overline{L(U(\overline{w}))}\\
\log \omega(w)&:=i\frac{\pi}{2}(U(\overline{w})-\nabla \sigma(U(\overline{w})))+i\frac{\pi}{2}\overline{L(U(\overline{w}))}\\
\end{align*}

and $\Upsilon(z)=\gamma(f(z))$ and $\Om(z)=\omega(f(z))$. Since $\gamma(w)+\omega(w)\equiv 1$ differentiating this relation with respect to $\frac{\dv}{\dv w}$ gives

\begin{align*}
\frac{\pi}{2}(1-\mathcal{H}'(w))e^{\frac{\pi}{2}(w-\mathcal{H}(w))}+i\frac{\pi}{2}(1+\mathcal{H}'(w))e^{i\frac{\pi}{2}(w+\mathcal{H}(w))}=0
\end{align*}

or equivalently

\begin{align*}
\end{align*}

Solving for $\mathcal{H}'(w)$ gives using that $\omega=1-\zeta$

\begin{align*}
\mathcal{H}'(w)=-\frac{\zeta(w)+i\omega(w)}{-\zeta(w)+i\omega(w)}=-\frac{i+(1-i)\zeta(w)}{i-(1+i)\zeta(w)}=-\frac{1+(i-1)\omega(w)}{-1+(1+i)\omega(w)}
\end{align*}

This completes the proof.

§ LOZENGE TILINGS OF A HEXAGON: ASYMPTOTIC ANALYSIS

In this section we return to the lozenge tilings of the hexagon with weight (<ref>).  As discussed in Introduction, the position of the white points in Figure <ref> from the lozenge tiling with the $q$-Racah weight can be described as a dynamic $q$-Racah ensemble. This is shown in Theorem 7.5. in [7] and we will briefly recall their definition. We then use Theorems <ref>, <ref> and <ref> to study the asymptotic behavior of this point process.

§.§ Extended $q$-Racah ensemble

Let $\mathbf a, \mathbf b, \mathbf c\in \mathbb N$ and consider the probability measure on all possible tilings of the hegagon defined by the weight (<ref>).
Let  $\{(\unscaledtime,y_j(\unscaledtime))\}_{\unscaledtime=1,j=1}^{\mathbf{b}+\mathbf{c}-1,\mathbf{a}}$ be the point process for the lozenge tiling model, which describes the positions of the white points in Figure <ref> shifted by $\frac 12$.   As mentioned in the Introduction, for fixed $1\leq s \leq {\bf b}+{\bf c}-1$, the points $\{y_j(s)\}_{j=1}^{{\mathbf a}}$ form a $q$-Racah ensemble. Following [7], the parameters for the ensemble, $\alpha$, $\beta$, $\gamma$ and $\delta$ can be expressed in terms of the parameters $\kappa, {\bf a}, {\bf b}$ and ${\bf c}$ as indicated in Table <ref>.  There are four cases in Table <ref> depending on $s$ and the dimensions of the hexagon ${\bf a},{\bf b}$ and ${\bf c}$.  If one assumes $\mathbf b<\mathbf c$, only  cases i) iii) and iv) are involved, whereas if $\mathbf b>\mathbf c$, only  cases iii) and iv) are involved. If $\mathbf b=\mathbf c$ then only cases i) and 4) occurs.   Important to note is that we always count from the bottom of the hexagon. This means that for $s>{\bf b}$, we shift the $q$-Racah ensemble by $s-{\bf b}$, which happens for cases iii) and iv).

In [7] the probability distribution for the extended process $\{(\unscaledtime,y_j(\unscaledtime))\}_{\unscaledtime=1,j=1}^{\mathbf{b}+\mathbf{c}-1,\mathbf{a}}$ was also computed and we will recall their definitions.
The starting point are transition matrices, for $s=0,\ldots,{\mathbf b}+{\mathbf c}-1$,

\begin{align*}
T_{\unscaledtime}(y_i,y_j)  \coloneqq\sum_{k=0}^{{\bf a}-1}\frac{c_{k,\unscaledtime+1}}{c_{k,\unscaledtime}}r_{k,\unscaledtime}(y_i)r_{k,\unscaledtime+1}(y_j),
\end{align*}

where, for $s=0,\ldots, {\mathbf b}+{\mathbf c}$,  the function $r_{k,\unscaledtime}(y) $ is the $q$-Racah orthonormal polynomial with parameters specified in Table <ref>.  So  the $r_{k,s}$ are  polynomials in $\nu_s(y)$ orthogonal with respect to the $q$-Racah weight $w_{n,\unscaledtime}(y)$ as defined in (<ref>). Important to note is that $s > {\bf b}$ we shifted the  $q$-Racah polynomial and its weight $w$ by $s-{\bf b}$. To avoid cumbersome notation we do not indicate this shift explicitly and we trust this does not lead to confusion. The coefficients $c_{k,s}$ are determined by

\begin{align} \label{eq:expressionsck}
c_{k,0}  \coloneqq \sqrt{\prod_{r=0}^k\frac{1-\mathbf{q}^{-(\mathbf{a}+\mathbf{b}+\mathbf{c}-r)}}{1-\mathbf{q}^{-(\mathbf{a}-r)}}},\qquad
c_{k,\unscaledtime} \coloneqq c_{k,0}\sqrt{\prod_{r=0}^{\unscaledtime-1}(1-\mathbf{q}^{-\mathbf{a}-r+k})(\mathbf{q}^{(\mathbf{a}+\mathbf{b}+\mathbf{c}-r-k-1)}-1)} .
\end{align}

The probability measure on the point process becomes the so called dynamic $q$-Racah polynomial ensemble, which is proportional to

\begin{equation} \label{eq: tiling dpp}
\sim \prod_{\unscaledtime=0}^{\mathbf{b}+\mathbf{c}-1}\det(T_{\unscaledtime}(y_{i}(\unscaledtime),y_{j}(\unscaledtime+1)))_{i,j=1}^{\mathbf{a}} \prod_{\unscaledtime=1}^{\textbf{b}+\mathbf{c}-1}\prod_{j=1}^{\mathbf{a}} w_{\unscaledtime}(y_j(\unscaledtime)),
\end{equation}

with $y_j(0)=j-1$ and $y_j({\mathbf b}+{\mathbf c})={\mathbf c}+j-1$ fixed for $j=1,\ldots,{\mathbf a}$.

$\alpha$     $\beta$     $\gamma$     $\delta$     $y$
i)    $\unscaledtime<\min\{\mathbf{b},\mathbf{c}\}$,
$0\leq y+\min\{0,\mathbf{b}-\unscaledtime\}\leq \unscaledtime+\mathbf{a}-1$      $\mathbf{q}^{-\mathbf{c}-\mathbf{a}}$     $\mathbf{q}^{-\mathbf{b}-\mathbf{a}}$     $\mathbf{q}^{-\unscaledtime-\mathbf{a}}$      $\kappa^2\mathbf{q}^{-\mathbf{c}+\mathbf{a}}$     $y$

ii)    if $\mathbf{b}>\mathbf{c}$,   $\mathbf{c}-1<\unscaledtime<\mathbf{b}-1$,
$0\leq y+\min\{0,\mathbf{b}-\unscaledtime\}\leq \mathbf{c}+\mathbf{a}-1$      $\mathbf{q}^{-\unscaledtime-\mathbf{a}}$     $\mathbf{q}^{\unscaledtime-\mathbf{a}-\mathbf{b}-\mathbf{c}}$     $\mathbf{q}^{-\mathbf{a}-\mathbf{c}}$      $\kappa^2\mathbf{q}^{-\unscaledtime+\mathbf{a}}$     $y$

iii)    if $\mathbf{b}\leq \mathbf{c}$, $\mathbf{b}-1<\unscaledtime<\mathbf{c}$,
$0\leq y+\min\{0,\mathbf{b}-\unscaledtime\}\leq \mathbf{b}+\mathbf{a}-1$      $\mathbf{q}^{\unscaledtime-\mathbf{a}-\mathbf{b}-\mathbf{c}}$     $\mathbf{q}^{-\unscaledtime-\mathbf{a}}$     $\mathbf{q}^{-\mathbf{b}-\mathbf{a}}$      $\kappa^2\mathbf{q}^{\unscaledtime-\mathbf{b}-\mathbf{c}+\mathbf{a}}$     $\mathbf{b}-\unscaledtime+y$
iv)     $\max\{\mathbf{b},\mathbf{c}\}<\unscaledtime+1$,
$0\leq y+\min\{0,\mathbf{b}-\unscaledtime\}\leq \mathbf{b}+\mathbf{c}-\unscaledtime+\mathbf{a}-1$      $\mathbf{q}^{-\mathbf{b}-\mathbf{a}}$     $\mathbf{q}^{-\mathbf{c}-\mathbf{a}}$     $\mathbf{q}^{\unscaledtime-\mathbf{a}-\mathbf{b}-\mathbf{c}}$      $\kappa^2\mathbf{q}^{-\mathbf{b}+\mathbf{a}}$     $\mathbf{b}-\unscaledtime+y$

Parameters for Lozenge Tiling with $q$-Racah weight

To establish a probability measure one also needs a configuration of $({\mathbf q}, \kappa)$. Here we consider

* Real case: ${\mathbf q}>0, \kappa\in \mathbb{R}$ and $\kappa$ can not lie in $[{\mathbf q}^{-{\mathbf a}+\frac{1}{2}},{\mathbf q}^{\frac{{\mathbf b}+{\mathbf c}}{2}-\frac{1}{2}}]$ if $\mathbf q>1$ or $[{\mathbf q}^{\frac{\mathbf b+{\mathbf c}}{2}-\frac{1}{2}}, {\mathbf q}^{-\mathbf a+\frac{1}{2}}]$ if ${\mathbf q}<1$.

* Imaginary case : ${\mathbf q}>0, \kappa\in i \mathbb{R}$.

* Trigonometric case : ${\mathbf q}=e^{i \lambda}, \kappa=e^{i \eta}$. To establish a probability measure,  it is enough to let $\lambda,\eta$ be such that $\eta-\frac{\lambda({\mathbf b}+ {\mathbf c }-1)}{2}, \eta+\lambda(\mathbf{a}-\frac{1}{2}) \in(k\pi,(k+1)\pi)$ for some $k$ to be an integer.

* Singular trigonometric case: similarly in the trigonometric case with $k\pi \leq \eta-\frac{\lambda({\mathbf b}+ {\mathbf c }-1)}{2}, \eta+\lambda(\mathbf{a}-\frac{1}{2})\leq(k+1)\pi$ and at least one of the inequality being an equality.

* Singular real case:   $\mathbf q>0$ and $\kappa = {\mathbf q}^{-\mathbf{a}+\frac{1}{2}}$ or ${\mathbf q}^{\frac{{\mathbf b}+{\mathbf c}-1}{2}}$.

§.§ Large hexagon

To study the asymptotic behavior we need to scale the parameters. And the limit regime we are considering are as such

\begin{align}\label{eq: scale parameters}
\mathbf{a}=n, \quad \mathbf{b}=nb, \quad\mathbf{c}=nc, \quad\mathbf{q}=q^{\frac{1}{n}}
\end{align}

for some $b,c>0$ and $q>0$ or $q$ on the unite circle. We further scale the lattice by $\frac{1}{n}$ via

\begin{align}
y =nx,  &\qquad x\in\frac{1}{n}\mathbb{Z}\\
\unscaledtime=\lfloor nt \rfloor, &\qquad t\in[0,b+c].
\end{align}

Instead of looking at $(s,y)$ plane, we are looking at the $(t,x)$ plane. Set
We are going to study the process $\{(t,x_j(t))\}_{t,j}$.

Following the analysis of Section <ref> the limiting behavior of a random tiling is determined by the asymptotic behavior of the recurrence coefficients for the $q$-Racah polynomials and the coefficients$c_{k,s}$. We will summarize the relevant results first and apply them to find the limiting behavior of the height function and its fluctuations in later sections.  We will only work out the case of real and imaginary cases of parameter listed above, and we leave the trigonometric case to the reader

The recurrence coefficients were are already studied in Section <ref>. We will recall their limits and express them in the parameters of the tiling model. To this end, define

\begin{align} \label{eq:defaballcases}
a(\xi;t) = \sqrt{A(\xi)C(\xi)} \qquad, \qquad b(\xi;t) = -A(\xi)-C(\xi)+1+\gamma\delta,
\end{align}

where $A(\xi)$ and $C(\xi)$ are defined by

\begin{align}
A(\xi) & =\frac{\left(1-q^{-c+\xi-1}\right) \left(1-q^{\xi-t-1}\right) \left(1-q^{-b-c+\xi-2}\right) \left(1-\kappa ^2 q^{-b-c+\xi}\right)}{\left(1-q^{-b-c+2 \xi-2}\right)^2} \label{eq:recurrance A}\\
C(\xi) & = \frac{\left(1-q^\xi\right) \left(1-q^{-b+\xi-1}\right) \left(\kappa ^2 q^{1-c}-q^{-c+\xi-1}\right) \left(q^{-t-1}-q^{-b-c+\xi-2}\right)}{\left(1-q^{-b-c+2 \xi-2}\right)^2}.\label{eq:recurrance CC}
\end{align}

These expressions are the limits of the coefficients $A_{k,n}$ and $C_{k,n}$, but only in case i) and ii). For case iii) and iv) we find different limits, but somewhat remarkably, the expressions above still give the right limiting behavior for $a(\xi;t)$ and $b(\xi;t)$ (up to a constant) as stated in the following lemma.

Let $a(\xi;t)$ and $b(\xi;t)$ be as in (<ref>). Then

\begin{equation}
\lim_{\frac{j}{n}\to\xi} a_{k,n}^{(s)} =\begin{cases}
a(\xi;t), & 0< t\leq b,\\
q^{t-b} a(\xi;t), &b<t<b+c,
\end{cases}
\end{equation}

\begin{equation}
\lim_{\frac{j}{n}\to\xi} b_{k,n}^{(s)} =\begin{cases}
b(\xi;t), & 0< t\leq b,\\
q^{t-b} b(\xi;t), &b<t<b+c.
\end{cases}
\end{equation}

Moreover, with $\mu(t,x)=q^{-x}+\kappa^2q^{-t-c-x}$  we have

\begin{equation}
\mu_n(x) =\begin{cases}
\mu(x;t), & 0< t\leq b,\\
q^{t-b} \mu(x;t), &b<t<b+c.
\end{cases}
\end{equation}

This follows by a straightforward computation using (<ref>), (<ref>) and Table (<ref>). The details are left to the reader.

The following lemma deals with the limiting behavior of the coefficients $c_{k,s}$.

Let $t \in (0,b+c)$ and $k,l \in \mathbb N$. Then, with $s=\lfloor tn\rfloor $ we have

\begin{equation}
\lim_{n\to \infty}\frac{c_{n+k,t}}{c_{n+l,t}}=e^{\tau(t)(l-k)}
\end{equation}

\begin{equation}
\tau(t)  = \frac{1}{2}\log\left|\frac{q^{t}-1}{1-q^{t-b-c}}\right| \label{eq: tau}.
\end{equation}

For $k<l$ and $\unscaledtime>l-k+1$ we obtain, from (<ref>)

\begin{align}
% \frac{c_{k,\unscaledtime+1}}{c_{k,\unscaledtime}} &  =\sqrt{-(1-\mathbf{q}^{-\mathbf{a}-\unscaledtime+k})(1-\mathbf{q}^{(\mathbf{a}+\mathbf{b}+\mathbf{c}-\unscaledtime-k-1)})}\\
\frac{c_{k,\unscaledtime}}{c_{l,\unscaledtime}} & =\frac{c_{k,0}}{c_{l,0}} \sqrt{\prod_{r=k+1}^{l}\frac{(1-\mathbf{q}^{-(\mathbf{a}+\unscaledtime-r)})(\mathbf{q}^{\mathbf{a}+\mathbf{b}+\mathbf{c}-r}-1)}{(1-\mathbf{q}^{-(\mathbf{a}-r)})(\mathbf{q}^{\mathbf{a}+\mathbf{b}+\mathbf{c}-\unscaledtime-r}-1)}} = \sqrt{\mathbf{q}^{\unscaledtime(l-k)}\prod_{r=k+1}^{l}\frac{1-\mathbf{q}^{-(\mathbf{a}+\unscaledtime-r)}}{1-\mathbf{q}^{-(\mathbf{a}+\mathbf{b}+\mathbf{c}-s-r)}}}
\end{align}

In the first equality we used the second expression in (<ref>) and the fact that we have a telescoping product. In the second, we used the explicit expression for $c_{k,0}$ and $c_{l,0}$.
Therefore, with $k<l$, $t\in[0,b+c]$ and $\unscaledtime \coloneqq \lfloor nt\rfloor$,  we have that

\begin{align}
\frac{c_{n+k,\unscaledtime}}{c_{n+l,\unscaledtime}} & = \sqrt{q^{\unscaledtime (l-k)/n}\prod_{r=k+1}^{l}\frac{1-q^{-(\unscaledtime/n-r/n)}}{1-q^{-(b+c-\unscaledtime/n-r/n)}}} \to \sqrt{\left|\frac{q^{t}-1}{1-q^{-(b+c-t)}}\right|^{l-k}},
\end{align}

as $n \to \infty$, and this proves the statement.

Real Case $q=0.05, \kappa=0.01$

Imaginary Case $q=0.5, \kappa=29i$

Singular Trigonometric Case
$q=e^{i(\frac{\pi}{2}-\frac{1}{n})}, \kappa=e^{i\frac{\pi}{2}}$

Singular Real Case $q=e, \kappa=e$

Different Cases of $q$-Racah Weights with size $n=100$

§.§ Limit Shape

We will now compute the limiting height function. We recall the definition of the height function $h_n$ in (<ref>). Then for any fixed $(t,x)$, we can rewrite the height function as

\begin{align*}
h_n(t,x)=\sum_{j=1}^{\mathbf{a}}\chi_{(-\infty, x]}(x_j(t)) =X_n(\chi_{(-\infty, x]}),
\end{align*}

where  $\chi_{A}$ stands for the characteristic function of the set $A$. Hence $h_n$ is the linear statistic for the step function with the position of the step at $x$.

We recall from Section <ref> that the limiting density of points can be computed by solving

\begin{equation}\label{eq:centraleqall}
\end{equation}

where $a(\xi;t)$ and $b(\xi;t)$ are defined as in (<ref>). We recall that these expressions are only valid for $0<t\leq b$, but by Lemma <ref> we see that for $t>b$ equation (<ref>) has an overall factor $q^{t-b}$ which we can divide out. Thus, (<ref>) is the correct equation for all $t$.

Let $h_n$ be the height function for the rescaled $q$-Racah lozenge tiling in the hexagon defined in (<ref>). For any given point $(t,x)$ we have the almost sure convergence of the height function

\begin{align}
\frac{1}{n}h_n(t,x)\to \int_{\max\{0,t-b\}}^x\rho(t,u)du
\end{align}

as $n\to\infty$
where $\rho(t,x)$ is the density defined by

\begin{align}\label{eq:limitshapetiling}
\rho(x)= \begin{cases} 1 & \quad \xi_-(x)<\xi_+(x)\leq 1\\ \frac{1}{\pi}\arccos\left(\frac{2y(1)-y(\xi_-(x))-y(\xi_+(x))}{\left|y(\xi_-(x))-y(\xi_+(x))\right|}\right)& \quad \xi_-(x)<1<\xi_+(x) \\ 0 &\quad 1\leq \xi_-(x) <\xi_+(x),
\end{cases}
\end{align}

where $\xi_-(x)\leq \xi_+(x)$ are the solutions to (<ref>) and $y(\xi)=q^{-\xi}+q^{\xi-b-c-2}$.

Since the hexagon is compact domain the moment condition (<ref>) is satisfied.
Hence, the conditions in Theorem <ref> are satisfied. Note that $a(\xi;t)$ has is analytic in $\xi$ and (<ref>) is satisfied. Note also that a step function is Riemann integrable. From here we can apply our Theorem <ref> and we use Proposition <ref> to further compute the limiting density $\rho(t,x)$.

Examples of $\rho(t,x)$ are graphed in Figure <ref>.  As expected, the corners are frozen in all cases and the density there is either zero or one.    The arctic curve is characterised by

\begin{align}
\{(t,x): 4a(1;t)^2-(\mu(t,x)-b(1;t))^2=0\}.
\end{align}

Note that this is not an algebraic curve in the variables $t$ and $x$, but it is a degree two curve in the variable variables $\mu(t,x)$ and $q^{-t}$.  Furthermore, the arctic curve has 6 points tangent to the hexagon in both real and imaginary cases. In the singular trigonometric case there are 2 tangential points and two non-differential points arctic curve, while in the singular real case there are 5 tangential points and one non-differential point on the arctic curve. The tangential points $(t,x )$ can be classified in two ways.
Intersection points of the arctic curve and the boundary of hexagon are given by

\begin{align}
t=0 \text{ and }&  x= \mu^{-1,(0)}(b(1;0)) \label{eq: intersection 1} \\
t=b+c \text{ and }&  x= \mu^{-1,(b+c)}(b(1;b+c)) \label{eq: intersection 2}
\end{align}

\begin{align}
t= & \log_q \left(\frac{\kappa ^2 q^{b+c+1}-\kappa ^2 q^{b+c+2}-q^{b+c}+q^{b+2 c+1}-\kappa ^2 q^{c+1}+\kappa ^2 q}{-q^{b+c+1}+q^{b+2 c+1}-\kappa ^2 q^{c+2}+q^{c+1}-q^c+\kappa ^2 q}\right) \text{ and }  x= 0 \\
t= & \log_q \left(-\frac{q^c \left(-q^{b+c}+q^{b+c+1}-\kappa ^2 q^{b+2}+q^b+\kappa ^2 q-1\right)}{\kappa ^2 q^{b+c+2}-q^{b+c+1}-\kappa ^2 q^{c+2}+q^c+\kappa ^2 q^2-\kappa ^2 q}\right) \text{ and }  x= 1+t \\
t= & \log_q \left(\frac{q^{b+c} \left(-q^{b+c+1}+q^{2 b+c+1}-\kappa ^2 q^{b+2}+q^{b+1}-q^b+\kappa ^2 q\right)}{\kappa ^2 q^{b+c+1}-\kappa ^2 q^{b+c+2}-q^{b+c}+q^{2 b+c+1}-\kappa ^2 q^{b+1}+\kappa ^2 q}\right) \text{ and }  x= t-b \\
t= & \log_q \left(-\frac{q^c \left(\kappa ^2 q^{b+c+2}-q^{b+c+1}-\kappa ^2 q^{b+2}+q^b+\kappa ^2 q^2-\kappa ^2 q\right)}{-q^{b+c}+q^{b+c+1}-\kappa ^2 q^{c+2}+q^c+\kappa ^2 q-1}\right) \text{ and }  x= 1+c.
\end{align}

The desired points are the intersection of the arctic curve with the boundary of the hexagon. For $t=0$  we have $A(1)=0$. For  $t=b+c$ we have $C(1)=0$. Hence in both cases we have $a(1;0)=a(1;b+c)=0$. Thus two points as (<ref>) and (<ref>) are found and given by $b(1;0)$ and $b(1;b+c)$. For the other points, we note that  on the boundary of the hexagon the discriminant defined in (<ref>) is always zero. Recall $\xi_-(x)$ and $\xi_+(x)$ defined in Proposition <ref>. This means at these points $\xi_-(x)=\xi_+(x)$. Moreover, the point $(t,x)$ lying on the arctic curve also means either $y(\xi_-(x))=y(1)$ or $y(\xi_+(x))=y(1)$. Then solve the equations $y(\xi_+(x_0))=y(1)$ with respect to $t$ for $x_0=0,1+t,t-b, 1+c$ respectively and we will find the rest of the four points.

Real Case with size $q=0.05$, $\kappa=0.01$

Imaginary Case with $q=0.5$, $\kappa=29i$

Singular Trigonometric Case with $q=e^{i\pi/2}$, $\kappa=e^{i\pi/2}$

Singular Real Case $q=e$, $\kappa=e$

Limiting densities  (<ref>) for $q$-Racah weight with $b=c=1$

§.§ Height Fluctuation and Gaussian Free Field

In this subsection we will show how the fluctuation of the height function in the liquid region  can be described by  Gaussian Free Field on the strip $S= \mathbb{R}\times (0,\pi)$, with Dirichlet boundary conditions.
Recall that on the Sobolev space $W^{1,2}_0(S)$ we have a bilinear paring

\begin{align*}
(f,g)_\nabla=\pi \int_S\langle \nabla f(x),\nabla g(x)\rangle dx.
\end{align*}

This bilinear paring can be extended to a well-defined bilinear pairing between the space of distributions $W^{-1,2}(S)\subset \mathcal{D}'(S)$ (where $\mathcal{D}'(S)$ denotes the space of distributions on $S$ in the sense of L. Schwartz) and $W^{1,2}_0(S)$ which we still denote by $\langle \cdot,\cdot \rangle_\nabla$, (see [1] for definitions on negative order Sobolev spaces).
The Gaussian Free Field (we refer to [41] for more details) is a random variable $\Phi$ with values in $W^{-1,2}(S)$ such that for every $\phi\in W^{1,2}_0(S)$ the random variable

\begin{equation} \label{eq:dirichlet}
\langle \Phi,\phi\rangle_\nabla \sim N\bigg(0,\pi\int |\nabla \phi|^2\bigg).
\end{equation}

For our statement, we need to introduce new coordinates for the liquid region

\begin{align*}
\mathcal{E} \coloneqq \left\{(t,x):4a(1; t)^2-\left(\mu(t,x)-b(1; t)\right)^2>0\right\}.
\end{align*}

To this end, with $S= \mathbb{R}\times (0,\pi)$, we define

\begin{equation}
F  :\mathcal{E}\to
S : (t,x) \mapsto(\tau,\theta) =\left(\tau(t),\arccos\left(\frac{\mu(t, x)-b(1; t)}{2a(1; t)}\right)\right), \label{pullback_qracah}
\end{equation}

where $\tau$ is as defined in (<ref>). One readily verifies that $F$ is a diffeomorphism  with

\begin{equation}
F^{-1}  :S\to\mathcal{E}:
(\tau,\theta) \mapsto (t,x)=(t(\tau),\mu^{-1, (t(\tau))}\left(b\left(1; t(\tau)\right)+2a\left(1; t(\tau))\cos\theta\right)\right),
\end{equation}

where $\mu^{-1,(t)}(x)$ is the inverse of $x\mapsto \mu(t,x)$.

We will show that $h_n\circ F^{-1} -\mathbb E[ h_n \circ F^{-1}]$ converges to the Gaussian Free Field on $S$ (with Dirichlet boundary conditions). More precisely, we will prove that
$$\langle h_n\circ F^{-1},\phi \rangle_\nabla -\mathbb E [\langle h_n \circ F^{-1}, \phi \rangle_\nabla] \to \langle \Phi, \phi\rangle_\nabla,$$

for all $C^2$ functions $\phi$ on $S$ such that $\Delta \phi$ has compact support and is piecewise constant in the vertical direction. That is, we assume there exists

0<t_0 <t_1 < \ldots < t_N<b+c,

such that

\begin{equation}\label{eq:piece}
\Delta  \phi(\tau, \theta)=
\begin{cases}
\Delta \phi(\tau(t_m),\theta), & \tau \in (\tau(t_{m-1}),\tau(t_m)], \qquad m=1,\dots,N\\
0,& \text{ otherwise.}
\end{cases}
\end{equation}

For such $\phi$, the pairing (<ref>) can be written as

\begin{align*}
\langle \Phi,\phi\rangle_\nabla  \coloneqq - \pi\sum_{m=1}^N(\tau(t_{m})-\tau(t_{m-1}))\int_0^{\pi}  \Phi(\tau(t_m),\theta) \Delta\phi(\tau(t_m),\theta)d\theta.
\end{align*}

Functions $\phi$ that satisfy (<ref>) are dense in the space of test functions with Dirichlet norm.

Consider a non singular $q$-Racah tiling model. That is the partial derivative $\mu_x(t,x)\neq 0$ for any $(t,x)$ inside the hexagon. Let $\phi$ be any twice differentiable function with compact support in $\mathbb{R}\times (0,\pi)$ such that $\Delta \phi$ is piecewise constant (as in (<ref>)).

\begin{align*}
\langle h_n\circ F^{-1}, \phi \rangle_\nabla-\mathbb{E}[\langle h_n\circ F^{-1}, \phi\rangle] \to \mathcal{N}\left(0,\pi\int |\nabla \phi|^2 \right)
\end{align*}

as $n\to\infty$ in distribution.

The proof is based on Theorem <ref>.

To apply Theorem <ref> we first show that $\langle h_n \circ F^{-1}, \Delta \phi\rangle$ is a linear statistic. To this end, note that $h_n(t,x)=\sum_{j=0}^{\mathbf{a}-1}\chi_{\{x_j(t)\leq x\}}$ where we used the notation $\chi_A$ for the characteristic function of the set $A$. Then,

\begin{align*}
\langle h_n\circ F^{-1}, \phi \rangle =-\pi\sum_{m=1}^N\sum_{j=1}^{n}(\tau(t_m)-\tau(t_{m-1}))\int_{\arccos \left(\frac{\mu(t_m,x_j(t_m))-b(1;t)}{2 a(1;t)}\right)}^\infty  \Delta  \phi(\tau(t_m),\theta) d \theta.
\end{align*}

If we thus define
$$f(t_m,y)=-\pi(\tau(t_m)-\tau(t_{m-1}))\int_{\arccos \left(\frac{\mu(t_m,x_j(t_m))-b(1;t)}{2 a(1;t)}\right)}^\infty  \Delta  \phi(\tau(t_m),\theta) d \theta,$$
then  the above can be rewritten as a linear statistic for the  point process $\{(t_m,x_j(t_m))\}_{m=1,j=0}^{N,\mathbf{a}-1}$. Note that this is a marginal of (<ref>) and the marginal density can be computed using repeated application of the Cauchy-Binet identity and using orthogonality of the polynomials. The result is exactly of the form (<ref>)
and (<ref>) with the $q$-Racah polynomials $r_{k,t_m}$ and coefficients $c_{k,t_m}$ as in (<ref>).

Now we are to check the conditions of Theorem <ref>. First, the hexagon forms a compact domain, hence (<ref>) is satisfied. Recall that $\mu_n(t,x)=q^{-x}+ \kappa^2 q^{-c-t+x+\frac{1}{n}}$, which converges uniformly to $\mu(t,x)=q^{-x}+ \kappa^2 q^{-c-t+x}$. Hence (<ref>) and (<ref>) is satisfied. Moreover

\begin{align}
\mu(t,\mu_n^{-1,(t)}(y))-y = \frac{1}{2} q^{-1/n} \left(q^{1/n}-1\right) \left(\sqrt{y^2-4 \kappa^2q^{-c-t+1/n}}-y\right),
\end{align}

and for any $y_1,y_2\in\{\mu_n(t,x): (t,x)\in\left(\frac{1}{n}\mathbb{Z}\times(\frac{1}{n}\mathbb{Z}+\frac{1}{2n})\right) \cap Hexagon\}$, we have $|y_1+y_2|$ is uniformly bounded. Moreover, for non singular $q$-Racah, we have $\mu_x(t,x)\neq 0$ for any $(t,x)$ inside the hexagon. Since $\mu(t,x)^2-4\kappa^2q^{-c-t}=(\log q)^{-2} \mu_x (t,x)^2$ and compactness, this implies that

\begin{align} \label{eq:lip bound y1 y2}
\left|\frac{\sqrt{y_1^2-4\kappa^2q^{-c-t-1/n}}-\sqrt{y_2^2-4\kappa^2q^{-c-t-1/n}}}{y_1-y_2} \right| &=  \left|\frac{y_1+y_2}{\sqrt{y_1^2-4\kappa^2q^{-c-t-1/n}}+\sqrt{y_2^2-4\kappa^2q^{-c-t-1/n}}} \right|\\
&\leq \frac{\sup_{(t,x)}|\mu(t,x)|}{\inf_{(t,x)}\sqrt{\mu(t,x)^2-4\kappa^2q^{-t-c}}}+O(n^{-1}), \nonumber
\end{align}

as $n\to \infty$,  which is a uniform bound. In the last step we used that $y_j=\mu(t,x_j) $ for some $x_j$ and $j=1,2$. Thus $L(\mu(t,\mu_n^{-1,(t)}))-id)=O(n^{-1})$, which shows the condition (<ref>). Moreover, (<ref>) and (<ref>) hold by Lemmas <ref> and <ref>.

Concluding, all conditions of Theorem <ref> are satisfied and we get the following convergence in distribution,

\begin{align*}
X_n(f)- \mathbb{E}[X_n(f)] \to \mathcal{N}\left(0,\sigma^2(f\circ\mu^{-1})\right),
\end{align*}

as $n\to \infty,$
where $\sigma^2(f\circ\mu^{-1})$ is the variance in the limit as defined as (<ref>) and (<ref>).  It remains to show $\sigma^2(f\circ\mu^{-1})$ matches the desired variance.

To this end,  we first recall a standard Fourier transform:

\begin{align*}
e^{-k|\tau|}=\frac{1}{\pi}\int_{-\infty}^\infty \frac{k}{k^2+\omega^2}e^{-i\omega \tau} d\omega.
\end{align*}

This allows us to rewrite the variance  as

\begin{align}\label{eq: variance Fourier}
\sigma^2(f\circ\mu^{-1}) = \frac{1}{\pi}\sum_{k=1}^\infty\int_{\mathbb{R}}\frac{k^2}{\omega^2+k^2}\left|\sum_{m=1}^Ne^{-i\tau(t_m)\omega}\widehat{f\circ \mu^{-1}}_k^{(m)}\right|^2d\omega
\end{align}

Next note that by construction we have

\begin{align*}
\theta\left(t,\mu^{-1,(t)}(b(1; t)+2a(1; t)\cos(\theta))\right)=\theta,
\end{align*}

which implies

\begin{align}
f(t,\mu^{-1,(t_m)}(b(1; t_m)+2a(1; t_m)\cos\theta)) & =- \pi \left(\tau(t_m)-\tau(t_{m-1})\right)\int_\theta^{\pi}\Delta \phi(\tau(t_m),\theta')d\theta',\nonumber
\end{align}

and thus

\begin{align*}
\widehat{f\circ \mu^{-1}}_k^{(m)} & = -\left(\tau(t_m)-\tau(t_{m-1})\right)\int_0^{\pi}\int_{\theta}^\pi \Delta \phi(\tau(t_m),\theta_q)d\theta_q\cos k\theta \ d\theta.
\end{align*}

Using integration by parts we get

\begin{align} \label{eq: fourier 1}
\widehat{f\circ \mu^{-1}}_k^{(m)}& = -\frac{\left(\tau(t_m)-\tau(t_{m-1})\right)}{ k}\int_0^\pi \Delta \phi(\tau(t_m),\theta)\sin k\theta \ d\theta.
\end{align}

By inserting (<ref>) into (<ref>), we find

\begin{multline} \label{eq:variance fourier last}
\sigma^2(f\circ\mu^{-1}) = \frac{1}{\pi}\sum_{k=1}^\infty\int_{\mathbb{R}}\frac{1}{\omega^2+k^2}\left|\sum_{m=1}^N \left(\tau(t_m)-\tau(t_{m-1})\right) e^{-i\tau(t_m)\omega}\int_0^\pi \Delta \phi(\tau(t_m),\theta)\sin k\theta \ d\theta\right|^2d\omega\\
= \frac{1}{\pi}\sum_{k=1}^\infty\int_{\mathbb{R}}\frac{1}{\omega^2+k^2}\left| \int_{\mathbb R} \int_0^\pi  e^{-i\tau\omega} \Delta \phi(\tau ,\theta)\sin k\theta \ d\theta d\tau\right|^2d\omega,
\end{multline}

where in the last step we used that $\Delta \phi$ is piecewise constant in $\tau$ and thus the sum over $m$ can be replaced by an integral. Finally, the innermost double integral at the right-hand side of (<ref>) is the Fourier transform on the strip S for functions that vanish on the boundary. Moreover, multiplication by  $(\omega^2+k^2)^{-1}$
on the Fourier transform is equivalent to applying $-\Delta ^{-1}$, and thus one one Laplace operator cancels. By Plancherel  we thus find that

$$\sigma^2(f\circ\mu^{-1})=-\pi \int \phi \Delta \phi$$
and this proves the statement.

For the singular case, the CLT theorem fails since (<ref>) is no longer bounded. The singular case allows the situation where there exist two different points $y_1, y_2$ such that $y_i = 2\sqrt{4\kappa^2q^{-c-t-1/n}}+O(n^{-1})$ for $i=1,2$. Hence

\begin{align*}
\frac{\sqrt{y_1^2-4\kappa^2q^{-c-t-1/n}}-\sqrt{y_2^2-4\kappa^2q^{-c-t-1/n}}}{y_1-y_2} =\frac{y_1+y_2} {\sqrt{y_1^2-4\kappa^2q^{-c-t-1/n}}+\sqrt{y_2^2-4\kappa^2q^{-c-t-1/n}}}
\end{align*}

diverges linearly with $n$.  Hence $L(\mu(t,\mu_n^{-1,(t)}))-id)$ is bounded below away from zero.

§.§ Complex Structure

In the remaing part of the paper we return to the complex structure of the liquid region. With the limiting height function at hand from Theorem <ref>, we mention  that, for $(t,x)$ in liquid region, the densities of the lozenges can be found by

\begin{align}
\rho_{I}(t,x) = & -\partial_t h(t,x)\pi, \label{eq:def_rhoi}\\
\rho_{II}(t,x) = &\partial_t h(t,x)\pi +\partial_x h(t,x)\pi, \label{eq:def_rhoii}\\
\rho_{III}(t,x) = & \pi-\partial_x h(t,x)\pi, \label{eq:def_rhoiii}
\end{align}

where $\rho_{I,II,III}$ stand for the limiting densities (up to factor $\pi$) of the lozenge of type $I$, $II$ and $III$ (this is based on fairly general principle and, for example, explained in  [14]). We included the factor $\pi$ for  convenience in the upcoming analysis. Now let $\Omega \in \mathbb C$ such that the angles of triangle $0,1$ and $\Omega$ has angle $\rho_I$ at $0$, $\rho_{II}$ at $\Omega$ and $\rho_{III}$ at $1$.  Note that we can write $\Omega$ explicitly as

\begin{align} \label{complexslope}
\Omega \coloneqq  \frac{\tan(\rho_{III})}{\tan(\rho_{III})+\tan(\rho_{I})}+i\frac{\tan(\rho_{III})\tan(\rho_{I})}{\tan(\rho_{III})+\tan(\rho_{I})}
\end{align}

We recall that $\Omega$ already appear in (<ref>) (in which $h=h_{tiles}$ and not the $h$ above, but they are related (<ref>)).
In the following we will show that this complex slope can be derived directly from the diffeomorphism $F$ in (<ref>). More precisely, we will show that the relation between $F$ and $\Omega$ as predicted from the combination of (<ref>) and (<ref>) indeed holds, by explicit computations on the recurrence coefficients.

Let $F$ be the diffeomorphism as (<ref>) and $\Omega$ as (<ref>). Then we have

\begin{align}\label{eq:complex structure}
\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}=-\frac{1+(i-1)\Om}{-1+(1+i)\Om},
\end{align}

where $\partial_z=(\partial_t-i\partial_x)/2$ and $\partial_{\bar{z}}=(\partial_t+i\partial_x)/2$ are the Wirtinger derivatives.

We will also show that the $\Omega$ satisfies the complex Burgers equation.
The complex Burgers' equation holds for $\Omega$,

\begin{align}\label{eq:Burgers}
\frac{\Omega_x}{\Omega}-\frac{\Omega_t}{1-\Omega}=-\log (q)\frac{ q^{c+t}+\kappa ^2 q^{2 x}}{q^{c+t}-\kappa ^2 q^{2 x}}.
\end{align}

In the remaining part we will prove Theorems <ref> and <ref>. Important is that in our proof we will only use the limiting behavior  recurrence coefficients of the $q$-Racah polynomials. With a little abuse of notation let

\begin{align} \label{eq:At}
A(t)  \equiv & A(\xi)|_{\xi=1} = \frac{\left(1-q^{-c}\right) \left(1-q^{-t}\right) \left(1-q^{-b-c-1}\right) \left(1-\kappa ^2 q^{-b-c+1}\right)}{\left(1-q^{-b-c}\right)^2} \\
C(t)   \equiv & C(\xi)|_{\xi=1} = \frac{\left(1-\frac{1}{q}\right) \left(1-q^{-b}\right) \left(1-\kappa ^2 q\right) q^{-c-t} \left(1-q^{-b-c+t}\right)}{\left(1-q^{-b-c}\right)^2}.  \label{eq:Ct}
\end{align}

We start the proofs by rewriting the diffeomorphism $F$ in terms of $A$ and $C$.

In the lozenge tiling with $q$-Racah weight, the diffeomorphism $F$ defined in (<ref>) can be rewritten as

\begin{align}
F(t,x)=\frac{1}{2}\log\left(\frac{A(t)}{C(t)}\right)+i\arccos{\frac{\mu(t,x)-1-\kappa^2 q^{-t-c}+A(t)+C(t)}{2\sqrt{A(t)C(t)}}} + \text{constant} \label{eq: working F}.
\end{align}

And $\frac{\partial_tF(t,x)}{\partial_xF(t,x)}$ lies in the lower upper plane for $(t,x)$ inside the liquid region.

By (<ref>), (<ref>) and (<ref>) we find

\begin{align*}
\tau(t) -\frac{1}{2}\log\left(\frac{A(t)}{C(t)}\right) = -\frac{1}{2}\log\left|\frac{(q^c-1)(1-q^{-b-c-1})(1-\kappa^2q^{-b-c+1})}{(1-q^{-1})(1-q^{-b})(1-\kappa^2q)}\right|,
\end{align*}

which we set to be the "$constant$". Recall that we have $a(1;t)=\sqrt{A(t)C(t)}$ and $b(1;t)=1+\gamma\delta-A(t)-C(t)$ for the recurrence coefficients. Then (<ref>) holds. Strictly speaking, this only only holds for $0<t\leq b$, but by Lemma <ref> we see that the extra factors in the argument of  the arccosine  cancel out.

To show the second statement we note that the imaginary part of $F$ only depends on $x$ and thus  $\Im\left(\frac{\partial_tF(t,x)}{\partial_xF(t,x)}\right)=-\frac{\partial_t\Re F(t,x)}{\partial_x\Im F(t,x)}$. Now we do the following computation

\begin{align*}
\partial_x \Im F(t,x)= & \frac{\partial_x\mu(t,x)2\sqrt{A(t)C(t)}}{\sqrt{4A(t)C(t)-(\mu(t,x)-1-\kappa^2 q^{-t-c}+A(t)+C(t))^2} }\\
\partial_t \Re F(t,x)= & \frac{\log(q)}{2}\left(\frac{q^t}{|q^t-1|}+\frac{q^{-(b+c-t)}}{|1-q^{-(b+c-t)}|}\right).
\end{align*}

Hence the imaginary part can be computed as

\begin{align*}
\Im\left(\frac{\partial_tF(t,x)}{\partial_xF(t,x)}\right) = -\frac{\log(q)}{\partial_x\mu(t,x)}\frac{\sqrt{4A(t)C(t)-(\mu(t,x)-1-\kappa^2 q^{-t-c}+A(t)+C(t))^2} }{4\sqrt{A(t)C(t)}}\left(\frac{q^t}{|q^t-1|}+\frac{q^{-(b+c-t)}}{|1-q^{-(b+c-t)}|}\right)
\end{align*}

Recall that $\mu(t,x)=q^{x}+\kappa^2 q^{-t-c} q^{-x}$. By the choice of the parameters, inside the hexagon we have $\frac{\partial_x\mu(t,x)}{\log(q)}$ is non-negative.

In the rest of the section we will use the subscript to denote the partial derivatives with respect to $x$ or $t$ and $'$ to denote the derivatives of single variable functions i.e.,

\begin{align*}
\mu_x \coloneqq \partial_x\mu, \qquad \mu_{xx} \coloneqq \partial_x^2\mu, \qquad \mu_{tx} \coloneqq \partial_x\partial_t\mu, \\
A'\coloneqq\frac{d}{dt}A, \qquad  A''\coloneqq\frac{d^2}{dt^2}A, \qquad C'\coloneqq\frac{d}{dt}C , \qquad  C''\coloneqq\frac{d^2}{dt^2}C.
\end{align*}

The following lemma is straightforward but will be important in our proofs.

Let $Y(t,x)=\mu(t,x)+1+\kappa^2 q^{-t-c}-A(t)-C(t)$.  $A(t), C(t)$ and $Y(t,x)$
satisfy the following differential equations

\begin{align}
Y_{xx}+2Y_{tx}+	Y_{tt}& =-\log(q) (Y_x+Y_t)\label{diff:mutx}\\
Y_{tt} & =-\log(q) Y_t\label{diff:mutt}\\
A''& =-\log(q)	\ A' \label{diff:recurrence1}\\
C''& =-\log(q)\ C'.\label{diff:recurrence2}
\end{align}

We compute that

\begin{align*}
A'(t)=\frac{\log (q)}{q^{t}-1}A(t) ,\quad C'(t)=-\frac{\log (q) q^{b+c}}{q^{b+c}-q^{t}}C(t),
\end{align*}

\begin{align*}
A''(t)=-\log (q) A'(t), \quad  C''(t)=-\log (q) C'(t).
\end{align*}

Recall that $\mu(t,x)=q^{-x}+\kappa ^2 q^{-c-t+x}$. Let $\tilde{\mu}(t,x)=\mu(t,x)+1+\kappa^2 q^{-t-c}$. Then

\begin{align*}
\partial_t^2\tilde{\mu}(t,x)=-\log(q)\partial_t\tilde{\mu}(t,x)\\
\partial_x^2\tilde{\mu}(t,x)+2\partial_x\partial_t\tilde{\mu}(t,x)=-\log(q)\partial_x\tilde{\mu}(t,x).
\end{align*}

This gives
$Y_{tt}=-\log(q) Y_t$ and $Y_{xx}+2Y_{tx}+Y_{tt}=-\log(q)(Y_x+Y_t)$.

We will also need the following.

\begin{align}
\frac{1}{\Omega(t,x) }=	\frac{\tan(\rho_{II})}{\tan(\rho_{II})+\tan(\rho_{I})}-i\frac{\tan(\rho_{II})\tan(\rho_{I})}{\tan(\rho_{II})+\tan(\rho_{I})}
\end{align}

Recall the definition of $\Omega$ in (<ref>). Then the statement holds by applying the following sum formula for tangent

\begin{align}\label{eq: sum tangent}
\tan(\rho_{I}+\rho_{II})=\frac{\tan (\rho_{I})+\tan( \rho_{II})}{1-\tan (\rho_{I})\tan (\rho_{II})},
\end{align}

and $\tan(\rho_{I}+\rho_{II})=-\tan(\rho_{III})$.

Now we are ready for the proof of Theorem. <ref>.

Consider the Möbius transformation from the complement of unit disc to the upper half plane given by

\begin{align}\label{mobius}
L(z)=\frac{(1-i) z+(1+i)}{i-i z}.
\end{align}

Instead of proving (<ref>) directly we show the equivalent statement

\begin{align}\label{mobiustoslope}
L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)=	\frac{1}{\Omega}=\frac{\tan(\rho_{II})}{\tan(\rho_{II})+\tan(\rho_{I})}-i\frac{\tan(\rho_{II})\tan(\rho_{I})}{\tan(\rho_{II})+\tan(\rho_{I})},
\end{align}

where  we used Lemma <ref> in the second equality. We denote the angles of the triangle formed by $0$,  $L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)$ and $1$, by $\hat \rho_I$ at 0, $\hat \rho_{II}$ at $L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)$ and $\hat \rho_{III}$ at $1$. The statement of the theorem is saying that $\rho_j=\hat \rho_j$ for and to show that this hold it is sufficient proof that $\hat \rho_1=\rho_1$ and $\hat \rho_{III}=\rho_{III}$.

my angle/.style = draw,
angle radius=5mm,
angle eccentricity=1.7,
right, inner sep=1pt,
(0,0) coordinate[label=above:$0$] (a) –
(5,0) coordinate[label=above:$1$] (c) –
(4,-3) coordinate[label=below:$	L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)$] (b) – cycle;
[my angle, "$\hat\rho_{I}$"] angle = b–a–c;
[my angle, "$\hat\rho_{II}$"] angle = a–c–b;
[my angle, "$\hat\rho_{III}$"] angle = c–b–a;
(7,0) coordinate[label=above:$0$] (a) –
(12,0) coordinate[label=above:$1$] (c) –
(11,-3) coordinate[label=below:$ \frac{1}{\Omega(t,x)}$] (b) – cycle;
[my angle, "$\rho_{I}$"] angle = b–a–c;
[my angle, "$\rho_{II}$"] angle = a–c–b;
[my angle, "$\rho_{III}$"] angle = c–b–a;

The triangles in the proof of Theorem <ref>. In the proof we show $\hat \rho_j=\rho_j$ and thus $L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)=1/\Omega$.

Before we come these proofs we  will first rewrite the left-hand side of (<ref>). To suppress the notation, we denote

\begin{align*}
Y(t,x) \coloneqq & \mu(x)-1-\kappa^2q^{-t-c} -A(t)-C(t).
\end{align*}

Note that $Y_x=\mu_x$ and $Y_{xx}=\mu_{xx}(x)$, which are the partial derivatives with respect to $x$ of the first and second order respectively. One can then compute the partial derivatives

\begin{align}
F_t(t,x) = & \frac{1}{2} \left(\frac{A'(t)}{A(t)}-\frac{C'(t)}{C(t)}\right)+i\frac{ \left(A'(t) C(t)+A(t)C'(t)\right)Y(t,x)-2 A(t)C(t) Y_t(t,x)}{2 (A(t) C(t))\sqrt{4A(t) C(t)-Y(t,x)^2}}, \\
F_x(t,x) = & -\frac{i Y_x(t,x)}{ \sqrt{4A(t) C(t)-Y(t,x)^2}}.
\end{align}

By the definition of $\partial_z, \partial_{\bar{z}}$ and $L$  we have

\begin{align}
L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)= \frac{F_t(t,x)}{F_x(t,x)}+1 .
\end{align}

Define two differential operators

\begin{align}\label{eq:def: d1 d2}
\mathcal{D}_1=\partial_t+\partial_x \quad \mathcal{D}_2=\partial_t .
\end{align}

Since $A(t)$ and $C(t)$ do not depend on $x$, we have $\mathcal{D}_1(AC)=\mathcal{D}_2(AC)$.  Then one can rewrite  $L$ as

\begin{align}
L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) & =\frac{2AC\mathcal{D}_1(Y)-\mathcal{D}_1(AC)Y}{2AC Y_x}+i\frac{(A'C-C'A)\sqrt{4AC-Y^2}}{2ACY_x} \label{eq:L structure 1}\\
L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) -1& =\frac{2AC\mathcal{D}_2(Y)-\mathcal{D}_2(AC)Y}{2AC Y_x}+i\frac{(A'C-C'A)\sqrt{4AC-Y^2}}{2ACY_x}.
\end{align}

Now we turn to the proofs of $\rho_I=\hat \rho_I$ and $\rho_{III}= \hat \rho_{III}$ starting with the latter.

STEP 1 : $\hat \rho_{III}=\rho_{III}$

Recall that we have that the density of the white points in (<ref>) and a location either  white point or lozenge of type III (in other words, they are complementary).  Thus the density of the type III lozenge is

\begin{align*}
\rho_{III}(t,x)= & \pi-\arccos\left(\frac{2y(1)-y(\xi_-(x))-y(\xi_+(x))}{\left|y(\xi_-(x))-y(\xi_+(x))\right|}\right).
\end{align*}

More explicitly, by (<ref>)  we can rewrite above to be

\begin{align*}
\rho_{III}(t,x)= & \arccos\left(\frac{\alpha_1}{\alpha_2}\right), \\
\alpha_1 \coloneqq & -\left(2\left( \frac{\mu'(x)}{\log(q)}\right)^2(q^{-1}+\alpha\beta q)-h(x)\right), \\
\alpha_2 \coloneqq & \sqrt{h^2(x)+4\left( \frac{\mu'(x)}{\log(q)}\right)^2p(x)},
\end{align*}

where $h(x)$ and $p(x)$ are defined by (<ref>) and (<ref>). By the trigonometric identity we have

\begin{align*}
\tan \rho_{III}(t,x)=\frac{\sqrt{\alpha_2^2-\alpha_1^2}}{\alpha_1} .
\end{align*}

By the relation (<ref>) we have $\alpha_2^2-\alpha_1^2=4\left(\frac{\partial_x\mu(t,x)}{\log(q)}\right)^2(q^{-1}-\alpha\beta q)^2(4 A(t) C(t)-Y(t,x)^2)$, where $(q^{-1}-\alpha\beta q)=(q^{-1}-q^{-b-c-1})$ and

\begin{align*}
\alpha_1= & 2 q^{-b-c-t-2 x-2} \left(-q^{b+c+t+1}+q^{b+c+x+1}+q^{b+t+x+1}+q^{b+x}-2 q^{b+2 x}+q^{t+x}-q^{t+1}\right) \\
& + \kappa^2 2 q^{-b-2 c-2 t-1} \left(\left(q^{b+1}+q^{t+1}+1+q^{t-c}\right) q^{c+t-x}+2 q^t \left(q^{b+c}+1\right)+q^{x-1} \left(q^{b+t+1}+q^b+q^t+q^{b+c+1}\right)\right) \\
& + \kappa^2 2 q^{-b-2 c-2 t-1} \left( -2 q^{t-1} \left(q^{b+1}+q^{t+1}+1\right)-2 q^c \left(q^{b+t+1}+q^b+q^t\right)\right) \\
& + \kappa^4 2 q^{-b-3 c-2 t-1} \left(q^{b+c+x+1}-q^{b+c+2 x}+q^{c+t+x+1}-2 q^{c+t+1}+q^{c+x}+q^{t+x}-q^{2 x}\right).
\end{align*}

Next we turn to $\hat \rho_{III}$. A simple computation shows that

\begin{align} \label{eq:rhoiii tri upsidedown}
\tan \hat \rho_{III}   = & \frac{-\Im\left(L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)\right)}{\Im\left(L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)\right)^2+\Re\left(L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)\right)(\Re\left(L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)\right)-1)}.
\end{align}

Then  by
evaluating $A(t)$, $C(t)$ and $Y(t,x)$ and inserting  it into (<ref>) we get

\begin{align}			L\left(\frac{\partial_{\bar{z}}F(t,x)}{\partial_{z}F(t,x)}\right)= & \frac{q^{c+2 t+1} \left(-q^{b+c}+2 q^t-1\right)+\left(q^{b+c+1}-q^{b+t+1}+q^b-q^t\right) q^{c+t+x}}{2 q \left(q^t-1\right) \left(q^t-q^{b+c}\right) \left(q^{c+t}-\kappa ^2 q^{2 x}\right)} \label{eq:Lexplicit}\\
& +\frac{\kappa ^2 q^{x+1} \left(q^{b+c+t+1}+\left(q^t-2\right) q^{b+c+x}+q^t \left(-q^{c+t+1}+q^c+q^x-q^t\right)\right)}{2 q \left(q^t-1\right) \left(q^t-q^{b+c}\right) \left(q^{c+t}-\kappa ^2 q^{2 x}\right)}\nonumber\\
& +i \frac{ q^{c+2 t+x}\left(q^{b+c}-1\right)\sqrt{(4 A(t) C(t)-Y(t,x)^2)}}{2 \left(q^t-1\right) \left(q^t-q^{b+c}\right) \left(q^{c+t}-\kappa ^2 q^{2 x}\right)} \nonumber
\end{align}

Then by  a cumbersome but straightforward computation we find that $\tan \rho_{III}= \tan \hat \rho_{III}$.

STEP 2 : $\hat \rho_I=\rho_I$

It is straightforward from (<ref>) that $\rho_I$ satisfies the following initial value problem

\begin{cases}
\partial_x \rho_{I}(t,x)=\mathcal D_2\rho_{III}\\
\tan \rho_{I}(t,x^*(t))=0.
\end{cases},$$

where $x^*(t)=\mu^{(-1),t}(b(1;t)-2a(1;t))$ is the boundary of the liquid region. In fact, this initial value problem defines $\rho_{I}$ uniquely. We will complete the proof by showing that $\hat\rho_I$ also solves this equation.

First let us define

\begin{align}
R_j\coloneqq & 2AC\mathcal{D}_j(Y)-\mathcal{D}_j(AC)Y \text{for }j=1,2\\
I\coloneqq & (A'C-C'A)\sqrt{4AC-Y^2},
\end{align}

and write

\begin{align*}
L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) =\frac{R_1}{2ACY_x}+i\frac{I}{2ACY_x} ,\quad
L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)-1 =\frac{R_2}{2ACY_x}+i\frac{I}{2ACY_x}.
\end{align*}

Then $ \mathcal{D}_1\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) \right) =  \mathcal{D}_1\arctan \frac{I}{R_1}$ and  $
\mathcal{D}_2\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) -1\right) =  \mathcal{D}_2\arctan \frac{I}{R_2}$. Using the for differential equations given by Lemma <ref> we have the following handy relations

\begin{align}
\mathcal{D}_j(R_j)  = &c_1 R_j -\mathcal{D}_j(AC)\mathcal{D}_j(Y)+2A'C'Y \label{eq: handy1}\\
\mathcal{D}_j(I) = & I\left(c_1  +\frac{2\mathcal{D}_j(AC)-Y\mathcal{D}_j(Y)}{4AC-Y^2}\right)\label{eq: handy2}
\end{align}

Surprisingly, these differential equations are the only ingredients of the rest of the proof.  We start with

\begin{align*}
\mathcal{D}_j\arctan \frac{I}{R_j}  =  \frac{I}{4AC-Y^2}\frac{\left(2\mathcal{D}_j(AC)-Y\mathcal{D}_j(Y)\right)R_j-(4AC-Y^2)\left(-\mathcal{D}_j(AC)\mathcal{D}_j(Y)+2A'C'Y\right)}{I^2+R_j^2}
\end{align*}

and mention that

\begin{align}
I^2+R_j^2=4AC\left(A'C'Y^2-\mathcal{D}_j(AC)\mathcal{D}_j(Y)Y+AC\mathcal{D}_j(Y)^2+\mathcal{D}_j(AC)^2-4AA'CC'\right). \label{eq: handy3}
\end{align}

Then we have for $j=1,2$

\begin{align}
\mathcal{D}_j\arctan \frac{I}{R_j} = \frac{YI}{2AC(4AC-Y^2)}
\end{align}

which shows

\begin{align}
\mathcal{D}_2\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) -1\right) =  \mathcal{D}_1\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) \right).\label{eq:Burgers0 imaginary}
\end{align}

Since $\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) -1\right)=\pi-\hat \rho_{III}+\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) \right)$. Hence $\hat \rho_1=-\arg\left(L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) \right)$ we have

\begin{align}		\mathcal{D}_1\hat{\rho_I}+\mathcal{D}_2(\pi-\hat {\rho_I}-\hat \rho_{III}) = 0.
\end{align}

And since we already proved that $\hat \rho_{III}= \rho_{III}$  and $\mathcal{D}_1-\mathcal{D}_2=\partial_x$, we find that

$$ \partial_x \hat \rho_{I}(t,x)=\mathcal D_1\rho_{III}.$$

Hence $\hat \rho_I$ and $\rho_I$ satisfy the same differential equation. Moreover, by (<ref>) we easily verify that they share the boundary value at $x^*(t)$. Since the initial value problem is obviously unique we conclude that $\hat \rho_I=\rho_I$ and this finishes the proof.

We end the paper with the proof of Theorem <ref>.

Theorem <ref> shows that $L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)=\frac{1}{\Omega}$. Then, by the chain rule, to show (<ref>) is equivalent to shows that  $	L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)$ is a solution to the following differential equation,

\begin{align}\label{Burgers0}
-\frac{\mathcal{D}_1\mathcal{L}}{\mathcal{L}}+\frac{\mathcal{D}_2\mathcal{L}}{\mathcal{L}-1}=  \frac{\mu_{xx}}{\mu_x}.
\end{align}

Theorem <ref> also shows that

\begin{align}
\end{align}

which is the imaginary part of (<ref>). To show (<ref>), it remains to show the real part. Recall that we define $R_j=2AC\mathcal{D}_j(Y)-\mathcal{D}_j(AC)Y$ and $I=(A'C-C'A)\sqrt{4AC-Y^2}$ and rewrite
\begin{align*}
L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right) =\frac{R_1}{2ACY_x}+i\frac{I}{2ACY_x} ,\quad
L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)-1 =\frac{R_2}{2ACY_x}+i\frac{I}{2ACY_x}.
\end{align*}

Now note that by (<ref>), (<ref>), one get

\begin{align*}
\frac{I\mathcal{D}_j(I)+R_j\mathcal{D}_j(R_j)}{I^2+R_j^2} = -\log(q) + \frac{I^2\left(2\mathcal{D}_j(AC)-Y\mathcal{D}_j(Y)\right)+R_j(4AC-Y^2)\left(-\mathcal{D}_j(AC)\mathcal{D}_j(Y)+2A'C'Y\right)}{(I^2+R_j^2)(4AC-Y^2)}.
\end{align*}

Further using the definition of $I, R_j$ to replace the numerator and the formula $I^2+R_j^2$ computed above as (<ref>) to replace the denominator, we have

\begin{align*}
\frac{I\mathcal{D}_j(I)+R_j\mathcal{D}_j(R_j)}{I^2+R_j^2} = -\log(q) + \frac{\mathcal{D}_j(AC)}{2AC}.
\end{align*}

\begin{align*}
\frac{1}{2}\mathcal{D}_j\log \left|\frac{I^2+R_j^2}{(2ACY_x)^2}\right| =  & \frac{I\mathcal{D}_j(I)+R_j\mathcal{D}_j(R_j)}{I^2+R_j^2}-\mathcal{D}_j\log(2AC)-\frac{\mathcal{D}_j(Y_x)}{Y_x}\\
= & -\log(q) - \frac{\mathcal{D}_j(AC)}{2AC}-\frac{\mathcal{D}_j(Y_x)}{Y_x}.
\end{align*}

Since $A $ and $C$ are only functions in $t$, $\mathcal{D}_1\log(AC)=\mathcal{D}_2\log(AC)$. Additionally, $\mathcal{D}_1(Y_x)-\mathcal{D}_2(Y_x)=Y_{xx}$ and hence we have

\begin{align}
-\mathcal{D}_1\log \left|L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)\right|+\mathcal{D}_2\log \left|L\left(\frac{\partial_{\bar{z}}F}{\partial_{z}F}\right)-1\right| = & \frac{\mu_{xx}}{\mu_x}, \label{eq:Burgers0 real}
\end{align}

i.e., the real part of the equation (<ref>). Then using the chain rule one can get (<ref>).

§.§ Acknowledgements

We thank Kari Astala for allowing us to present Theorem <ref> in this work. We also thank Christophe Charlier for allowing us to modify his program to generate pictures of random tiling in this paper. We thank Istvan Prause for pointing out reference [40].

[1]
R.A. Adams Sobolev spaces.
Academic Press [A subsidiary
of Harcourt Brace Jovanovich, Publishers], New York-London,
1975. Pure and Applied Mathematics, Vol. 65.

[2]
A. Aggarwal,
Universality for lozenge tiling local statistics,
arXiv preprint arXiv:1907.09991, 2019.

[3]
K. Astala, E. Duse, I. Prause, and X. Zhong,
Dimer models and conformal structures,
arXiv preprint arXiv:2004.02599, 2020.

[4]
K. Astala, T. Iwaniec and G. Martin,
Elliptic Partial Differential Equations and Quasiconformal Mappings in the Plane
Princeton University Press (2009) NJ.

[5]
J. Baik, T. Kriecherbauer, K. D. T-R. McLaughlin, and P. D. Miller,
Discrete Orthogonal Polynomials.(AM-164): Asymptotics and Applications (AM-164),
Princeton University Press, 2007.

[6]
D. Betea,
Elliptically distributed lozenge tilings of a hexagon,
SIGMA. Symmetry, Integrability and Geometry: Methods and Applications, Vol. 14, 2018, p. 032.

[7]
A. Borodin, V. Gorin, and E. M. Rains,
q-Distributions on boxed plane partitions,
Selecta Mathematica, Vol. 16, No. 4, 2010, pp. 731–789.

[8]
A. Borodin, I. Corwin, and T. Sasamoto,
Spectral theory for the q-Boson particle system,
Compositio Mathematica, Vol. 152, No. 10, 2016, pp. 2329–2382.

[9]
C. Boutillier and Z. Li,
Limit shape and height fluctuations of random perfect matchings on square-hexagon lattices,
Ann. Inst. Fourier (Grenoble) 71 (2021), 2305-2386.

[10]
J. Breuer and M. Duits,
The Nevai condition and a local law of large numbers for orthogonal polynomial ensembles,
Advances in Mathematics, Vol. 265, 2014, pp. 441–484.

[11]
J. Breuer and M. Duits,
Universality of mesoscopic fluctuations for orthogonal polynomial ensembles,
Communications in Mathematical Physics, Vol. 342, No. 2, 2016, pp. 491–531.

[12]
J. Breuer and M. Duits,
Central limit theorems for biorthogonal ensembles and asymptotics of recurrence coefficients,
Journal of the American Mathematical Society, Vol. 30, No. 1, 2017, pp. 27–66.

[13]
A. Bufetov and A. Knizel,
Asymptotics of random domino tilings of rectangular Aztec diamonds, Annales de l'Institut Henri Poincaré, 54 (2018), 1250-1290.

[14]
C. Charlier, M. Duits, A. B. Kuijlaars, and J. Lenells,
A periodic hexagon tiling model and non-Hermitian orthogonal polynomials,
Communications in Mathematical Physics, Vol. 378 (2020), 401–466.

[15]
S. Chhita and B. Young,
Coupling functions for domino tilings of Aztec diamonds,
Advances in mathematics, Vol. 259, 2014, pp. 173–251.

[16]
H. Cohn, R. Kenyon and J. Propp,
A variational principle for domino tilings,
J. Amer. Math. Soc.,
14 (2001), 297-346.

[17]
H. Cohn, M. Larsen and J. Propp,
The shape of a typical boxed plane partition,
New York J. Math. 4 (1998) 137-165.

[18]
D. De Silva and  O. Savin,
Minimizers of convex functionals arising in random surfaces
Duke Mathematical Journal, Vol. 151, No. 3 (2010), 487$-$532

[19]
P. Deift, T. Kriecherbauer, K. T-R. McLaughlin, S. Venakides, and X. Zhou,
Strong asymptotics of orthogonal polynomials with respect to exponential weights,
Comm. Pure Appl. Math., Vol. 52, No. 12, 1999, pp. 1491–1552.

[20]
P. Deift and X. Zhou,
A universality theorem for the spacing of zeros of orthogonal polynomials,
Internat. Math. Res. Notices, Vol. 5, 1998, pp. 285–299.
[21]
J.F. van Diejen and J.V. Stokman,
Multivariable $ q $-Racah polynomials,
Duke Math. J., Vol. 95, No. 1 (1998), 89–136.

[22]
E. Dimitrov and A. Knizel,
Log-gases on quadratic lattices via discrete loop equations and q-boxed plane partitions,
Journal of Functional Analysis, Vol. 276, No. 10, 2019, pp. 3067–3169.

[23]
J. Duchon, R. Robert, and V. Vargas,
Forecasting volatility with the multifractal random walk model,
Mathematical Finance: An International Journal of Mathematics, Statistics, and Financial Economics, Vol. 22.

[24]
P. Duren,
Harmonic Mappings in the Plane.
Cambridge University Press, 2004.

[25]
M. Duits,
On global fluctuations for non-colliding processes
The Annals of Probability,
Vol. 46, No. 3 (2018), 1279–1350.

[26]
V. Gorin,
Lectures on Random Lozenge Tilings
Cambridge Studies in Advanced Mathematics, Series Number 193

[27]
V. Gorin,
Nonintersecting paths and the Hahn orthogonal polynomial ensemble. Funct. Anal. Appl. 42 (2008), no. 3, 180–197.

[28]
V. Gorin and G. Olshanski
Determinantal measures related to big q-Jacobi polynomials
Functional Analysis and Its Applications,
Vol. 49 (2015), pp. 214-217.

[29]
V. Gorin and J. Huang
Dynamical Loop Equation
arXiv preprint arXiv:2205.15785,(2022)

[30]
A. Hardy,
Polynomial ensembles and recurrence coefficients,
Constructive Approximation,
Vol. 48, No. 1 (2018), 137–162.

[31]
K. Johansson,
Non-intersecting, simple, symmetric random walks and the extended Hahn kernel, Ann. Inst. Fourier 55 (2005), no. 6, 2129–2145.
[32]
K. Johansson, The arctic circle boundary and the Airy process,
Ann. Probab. 33 (2005), no. 1, 130.

[33]
K. Johansson, Random matrices and determinantal processes, Mathematical Statistical Physics, Elsevier B.V. Amsterdam (2006) 1–55.

[34]
R. Kenyon and  A. Okounkov,
Limit shapes and the complex Burgers equation. Acta Math., 199 (2007), 263-302.

[35]
R. Kenyon, A. Okounkov and S. Sheffield,
Dimers and amoebae.
Ann. of Math.,
163 (2006), 1019–1056.

[36]
R. Koekoek, P. A. Lesky, and R. Swarttouw,
Hypergeometric orthogonal polynomials and their q-analogues,
Springer Science & Business Media (2010).

[37]
N. Kuchumov,
A variational principle for domino tilings of multiply-connected domains
arXiv: 2110.06896

[38]
A.B. Kuijlaars and W. Van Assche,
The asymptotic zero distribution of orthogonal polynomials with varying recurrence coefficients,
Journal of Approximation Theory,
Vol. 99, No. 1 (1999), 167–197.

[39]

W. König, Orthogonal polynomial ensembles in probability theory, Probab. Surveys 2 (2005), 385–447.

[40]
A. Morales, I. Pak and M. Tassy,  Asymptotics for the number of standard tableaux of skew shape and for weighted lozenge tilings. Combinatorics, Probability and Computing,  31(4), 550-573. (2022).

[41]
E. Powell and W. Werner
Lecture notes on the Gaussian Free Field

[42]
L. Petrov
Asymptotics of uniformly random lozenge tilings of polygons. Gaussian free field
Annals of Probability,
Vol. 43, No. 1, pp. 1-43, 2015.

[43]
G. Schwarz
Hodge Decompositions - A method for solving boundary value problems
Springer Lecture Notes in Mathematics (1995) 1607

§ ANOTHER EXAMPLE: STATIONARY BIRTH-DEATH PROCESS

Let $X=\{0,1,\dots, M\}$ be a state space. Let

\begin{align*}
B,D:X\mapsto [0,+\infty) \\
\end{align*}

Define an $M+1$ dimensional tridiagonal matrix

\begin{align*}
-B(0)-D(0) & B(0) & 0 & \dots\\
D(1) & -B(1)-D(1) & B(1) & 0 & \dots\\
0 & D(2) & -B(2)-D(2) & B(2) & 0 & \dots \\
\dots\\
& & & 0 & D(M) & -B(M)-D(M)
\end{pmatrix}
\end{align*}

and a self-adjoint negative tridiagonal matrix

\begin{align}
\tilde{L}=\begin{pmatrix}
-B(0)-D(0) & \sqrt{B(0)D(1)} & 0 & \dots\\
\sqrt{B(0)D(1)} & -B(1)-D(1) & \sqrt{B(1)D(2)} & 0 & \dots\\
0 & \sqrt{B(1)D(2)} & -B(2)-D(2) & \sqrt{B(2)D(3)} & 0 & \dots \\
\dots
\end{pmatrix}
\end{align}

Let $\lambda_0,\dots,\lambda_N$ be eigenvalues for the non-negative operator $-\tilde{L}$ and the corresponding eigenvectors are $\varphi_0,\varphi_1,\dots, \varphi_N$. Then $L=D^{-1}\tilde{L}D$ where $D=\text{diag}(\varphi_0(0),\varphi_0(1),\dots,\varphi_0(M))$.  Do $\varphi_0(k)$ denote the entries of the eigenvector? I could not follow the above relation between $L$ and $\tilde L$ It is not hard to verify that $L$ is singular with this configuration and successive column operations.  What does this mean? Without loss of generality we let $\lambda_0=0$ and $\varphi_0(x)=\sqrt{\prod_{y=0}^{x-1}\frac{B(y)}{D(y+1)}}$, $\varphi_0(0)=1$. Then we define the transition matrix
\begin{align*}
P_t(x,y) = \frac{\varphi_0(y)}{\varphi_0(x)}\sum_{j=0}^M e^{-t\lambda_j}\varphi_j(x)\varphi_j(y)
\end{align*}

that is the probability of particles located at $y$ after time $t$ starting from position $x$. Then the process generated by $P_t$ has the Markov property, i.e.,

\begin{align*}
\sum_{y=0}^MP_s(x,y)P_t(y,z) = P_{s+t}(x,z)
\end{align*}

since $\tilde{L}$ is self adjoint and thus different eigenvectors are orthogonal. Moreover the stationary distribution can be found as,

\begin{align*}
\lim_{t\to\infty}P_{t}(x,y)=\varphi_0(y)^2
\end{align*}

The process described above is solvable if there exist a polynomials $p_j$ and some real function $\mu(x)$ (called sinusoidal coordinate) for all $j=0,\dots, N$ such that

\begin{align*}
\varphi_j(x) = \varphi_0(x)p_j(\mu(x)).
\end{align*}

Note that by construction , $p_j(\mu(x))$ is an orthonormal polynomial with weight $\varphi_0$. The solvability depends on $B, D$. For solvable models one can refer to <cit.>.  incomplete references
agraphStationary Birth Death Process on the Weyl Chamber

Let the Weyl Chamber be $W=\{x_1<x_2<\dots <x_n: x_j\in X\}$ for some $n<N$. Define an $n$ particles birth-death process on the Weyl chamber as

\begin{align*}
P_t(\mathbf{x},\mathbf{y}) = \frac{1}{Z_t}\frac{\prod_{1\leq i<j\leq n} \mu(y_i)-\mu(y_j)}{\prod_{1\leq i<j\leq n} \mu(x_i)-\mu(x_j)}\det (P_t(x_i,y_j))_{i,j=1}^n
\end{align*}

where by Cauchy-Binet formula and the orthonormality of $p_j$

\begin{align*}
Z_t & = \sum_{\mathbf{y}\in W} \frac{\det(p_i(\mu(y_j))_{i,j=1}^n \det(P_t(x_j,y_i))_{i,j=1}^n}{\det(p_i(\mu(x_j))_{i,j=1}^n} \\
& = n!\frac{\det(\sum_{y=0}^Mp_i(\mu(y))P_t(x_j,y))_{i,j=1}^n}{\det(p_i(\mu(x_j))_{i,j=1}^n}\\
& = n!e^{\sum_{j=0}^M-\lambda_jt}
\end{align*}

Then the stationary process can thus be defined as the one governed by $P_t$ above with the initial value

\begin{align*}
P_0(\mathbf{x}) & =\frac{1}{n!}\prod_{1\leq i<j\leq n}(\mu(x_i)-\mu(x_j))^2\prod_{j=1}^N\varphi_0(x_j)^2
\end{align*}

This is stationary since $\sum_{x\in W} P_0(\mathbf{x}) P_t(\mathbf{x},\mathbf{y}) = P_0(\mathbf{y})$ by Cauchy-Binet. Now consider a sequence $t_0<t_1<t_2<\dots<t_N\in\mathbb{R}$ with stationary initial condition at $t_0$. Then the probability measure for this Markov process becomes

\begin{align*}
\end{align*}

which, is (<ref>) with $w_m(x)=\varphi_0(x)$ and $c_{j,m}=e^{-\lambda_jt_m}$.
A concrete example is as the follows
agraphq-Charlier Stationary Birth-Death Process  incomplete reference

\begin{align*}
B=q^x \quad D=1-q^x
\end{align*}

Then the eigenvector of $-\tilde{L}$ consists of the q-Charlier orthonormal polynomials,  $p_k(\mu(x))$ whose recurrence relation is

\begin{align*}
\mu(x)p_k(\mu(x))=\sqrt{A_{k+1}C_k}p_{k+1}(\mu(x))-(A_k+C_k)p_k(\mu(x))+\sqrt{A_kC_{k-1}}p_{k-1}(\mu(x))\\
\mu(x)=q^{-x}-1, \quad A(k)=-q^{-2 k-1}  \quad C(k)=-q^{-2 k} \left(1-q^k\right) \left(q^k+1\right)
\end{align*}

In order to get a non-trivial process in the limit we need to scale $q$ to $q^{\frac{1}{n}}$, $x$ to $nx$ and $t$ to $nt$. Let $\xi=\lim_{n\to\infty}\frac{k}{n}$.  Then the limiting recurrence coefficients after scaling are as

\begin{align*}
A(\xi)=-q^{-2 \xi } \quad C(\xi)=-q^{-2 \xi } \left(1-q^{\xi }\right) \left(q^{\xi }+1\right)
\end{align*}

Theorem <ref> gives the limit shape density (see Figure <ref>)

\begin{align*}
\rho(t, x)= & \frac{-q^{-x}\log q}{\pi}\int_{\min(\xi_0(x),1)}^{1}\frac{1}{\sqrt{q^{-2 (\xi +x)} \left(-q^{2 \xi }-4 q^{2 x}+4 q^x\right)}}d\xi\\
= &  \begin{cases}
\frac{1}{\pi}\arccos \left(\frac{q}{\sqrt{-4 q^{2 x}+4 q^x}}\right) \quad &,  -q^{2 }-4 q^{2 x}+4 q^x>0\\
0 \quad &, -q^{2 }-4 q^{2 x}+4 q^x\leq0
\end{cases}
\end{align*}

where $\xi_0(x)$ is such that $-q^{2 \xi_0(x) }-4 q^{2 x}+4 q^x=0$.

Limiting Density of Paths of q-Charlier Stationary BD Process $q=0.6$

Moreover, the eigenvalues are $\lambda_{k}=1-q^{\frac{k}{n}}$ and we have $\tau(t)=t(-q\log q)$. The linear statistics has the Gaussian fluctuation  in the limit by Theorem <ref>,

\begin{align*}
X_n(f)-\mathbb{E}[X_n(f)]\to^d \mathcal{N}(0,\sum_{m_1,m_2=0}^N\sum_kke^{-|(t_{m_1}-t_{m_2})q\log q|k}\left(\widehat{f\circ \mu^{-1}}_k^{(t_{m_1})}\widehat{f\circ \mu^{-1}}_{-k}^{(t_{m_2})}\right))
\end{align*}

Then the diffeomorphism towards the Gaussian free field is given by

\begin{align*}
F(t,x) = \log (-tq\log q)+i \arccos\left(\frac{q^{2-x}-2}{2\sqrt{1-q^2}}\right)
\end{align*}

$\mathcal{L}(t,x)\equiv\frac{\partial_tF(t,x)}{\partial_x F(t,x)}+1=1+i\frac{ \sqrt{4 q^{x}-q^{2}-4q^{2x}}}{ -t q \log (q)}$.

§ INTRODUCTION TO Q-RACAH POLYNOMIAL

q-Shifted Factorial
Let $q\neq 0$ nor $1$.

* q-extended Pochhammer-symbol $(a;q)_k= (1-a)(1-aq)\dots (1-aq^{k-1})$. Whereas, $\lim_{q\to1}\frac{(q^\alpha,q)_k}{(1-q)^k}=\alpha(\alpha+1).\dots.(\alpha+k-1)\equiv(\alpha)_k$

* q-shifted factorials $(a;q)_{-k}= \frac{1}{\prod_{i=1}^k(1-aq^{-i})}$

Noting that the usual algebra is recovered by taking $q\to 1$.
agraphSome useful results:

\begin{equation}\label{q-shifted_factorial_q^-n}
(q^{-n};q)_k=\frac{(q;q)_n}{(q;q)_{n-k}}(-1)^kq^{n(n-1)/2-nk} \text{ , k=0,1,2\dots,n}
\end{equation}

Basic hypergeometric (or q-hypergemometric function) $_r\phi_s$ is defined as

\begin{equation}
_r\phi_s\bigg(\begin{matrix}a_1,\dots, a_r\\ b_1,\dots, b_s\end{matrix} ;q,z\bigg)= \sum_{k=0}^\infty \frac{(a_1,\dots, a_r;q)_k}{(b_1,\dots, b_s;q)_k}(-1)^{(1+s-r)k}q^{(1+s-r)\binom{k}{2}}\frac{z^k}{(q;q)_k}
\end{equation}

where $(a_1,\dots,a_r;q)_k=(a_1;q)_k.\dots.(a_r;q)_k$

§.§ q-Racah

Racah is among the top hierarchy in the Askey scheme with classic polynomials followed by, with orders, Hahn, Jacobi, Laguerre and Hermite polynomials. Also in the hexagon tiling, q-Racah provides the most general weights only depending on the vertical location of the lozenge giving an orthogonal polynomial ensemble. This gives an interpretation of "the height of the boxes". For biorthogonal one should look into more general elliptic weights.

\begin{align}
R_n(\nu(x)) & = _4\phi_3\bigg(\begin{matrix}q^{-n}, \alpha\beta q^{n+1},q^{-x},\delta \gamma q^{x+1}\\\alpha q,\delta \beta q,\gamma q\end{matrix}; q,q \bigg) \textbf{   ,  } n=0,1,2,\dots, -\log_q(\gamma q) \\
& = \sum_{k\geq 0}\frac{(q^{-n},\alpha\beta q^{n+1},q^{-x},\delta \gamma q^{x+1};q)_k}{(\alpha q, \delta \beta q,\gamma q;q)_k}\frac{q^k}{(q;q)_k} \\
& = \sum_{k\geq 0}\frac{(q^{-n};q)_k(\alpha\beta q^{n+1};q)_k(q^{-x};q)_k(\delta \gamma q^{x+1};q)_k}{(\alpha q;q)_k(\delta \beta q;q)_k(\gamma q;q)_k}\frac{q^k}{(q;q)_k} \nonumber  \\
\nu(x) & = q^{-x}+\delta \gamma q^{x+1}
\end{align}

* Here we consider the highest order of the polynomial to be $-\log_q(\gamma q) $.

* $R_n$ is a polynomial of order n with respect to $\nu(x)$ ($n=0,\dots, \gamma+1$)

Denote $p_n(\nu(x))$ by the monic q-Racah polynomials i.e. $p_n(\nu(x))=R_n(\nu(x))\frac{(q^{\alpha+1},\delta q^{\beta+1}, q^{-N};q)_n}{( q^{\alpha+\beta}q^{n+1};q)_n}$. And $r_n(x)$ the q-Racah orthonormal polynomials. Then the three term recurrence relation can be written as

\begin{align}
\nu(x)p_n(\nu(x)) & =p_{n+1}(\nu(x))+(1+\delta\gamma q-A_n-C_n)p_n(\nu(x))+A_{n-1}C_np_{n-1}(\nu(x)),
\end{align}

\begin{align*}
A_n & =\frac{(1-\gamma q^{n+1})(1-\alpha q^{n+1})(1- \alpha\beta q^{n+1})(1-\delta \beta q^{n+1})}{(1-\alpha\beta q^{2n+1})(1-\alpha\beta q^{2n+2})}\\
C_n & =\frac{q\left(1-q^n\right) \left(1-\beta  q^n\right) \left(\delta -\alpha  q^n\right) \left(\gamma -\alpha  \beta  q^n\right)}{(1-\alpha\beta q^{2n})(1- \alpha\beta q^{2n+1})}.
\end{align*}

Then by normalization we have

\begin{align}
\nu(x)r_n(x) & =\sqrt{A_nC_{n+1}}r_{n+1}(x)+(1+\delta\gamma q-A_n-C_n)r_n(x)+\sqrt{A_{n-1}C_n}r_{n-1}(x).
\end{align}

We can also rewrite it into

\begin{align}
\frac{\nu(x)-1-\gamma\delta q}{q^{-1}-1}r_n(x) & =\frac{\sqrt{A_nC_{n+1}}}{q^{-1}-1}r_{n+1}(x)+\frac{(-A_n-C_n)}{q^{-1}-1}r_n(x)+\frac{\sqrt{A_{n-1}C_n}}{q^{-1}-1}r_{n-1}(x) \label{normalised_recurrance}
\end{align}

From the last version of the recurrence relation we one can easily recover the classic Hahn polynomial recurrence relation by taking $q\to 1$ and $\delta=0$.

We further define the $\nu^{-1}$ function by

\begin{align}
x & \equiv -\frac{(1-q^{-\nu^{-1}(x)})(1-\delta \gamma q^{\nu^{-1}(x)+1})}{1-q^{-1}} \\
\nu^{-1}(x) & =-\frac{1}{\log(q)}\log\left(\frac{1+\delta \gamma q+(1-q^{-1})x+\sqrt{(1+\delta \gamma q+(1-q^{-1})x)^2-4\delta \gamma q}}{2}\right) \label{g_q_racah}\\
\nu^{-1}(x) & \to^{q\to 1}\frac{x}{1-\delta}
\end{align}

* There should be two solutions, and we pick one that works even for $\delta<0$. And such case is the purely imaginary q-Racah case.

* $\nu^{-1}$ is increasing in $x$ for all $q>0$.

* Let $\delta \to 0$, then $g(x)=\frac{\log(1+(1-q^{-1})x)}{\log(q^{-1})}$, which fits in the q-Hahn case.

§ VARIANCE STRUCTURE FOR THE GAUSSIAN FREE FIELD

Note that by construction we have

\begin{align*}
\theta_q\left(t,\mu^{-1,(t)}(b(1; t)+2a(1; t)cos(\theta))\right)=\theta
\end{align*}

By either inverse Fourier transformation of $e^{-k|\tau|}$ or residual theorem we have the identity

\begin{align*}
e^{-k|\tau|}=\frac{1}{\pi}\int_{-\infty}^\infty \frac{k}{k^2+\omega^2}e^{-i\omega \tau} d\omega
\end{align*}

Then the variance can be rewritten as

\begin{align}
\sigma^2(f\circ\mu^{-1}) & =\frac{1}{\pi}\sum_{k=1}^\infty\int_{\mathbb{R}}\frac{k^2}{\omega^2+k^2}|\int_Ie^{-i\tau(t)\omega}\widehat{f\circ \mu^{-1}}_k^{(t)}dt|^2d\omega \\
f(t,y) & =-\pi\int_y^{g(b(1; t)+2a(1; t))}\Delta\phi\left(\tau(t),\theta_q(t,x)\right)\frac{d(\tau,\theta_q)}{d(t,x)}dx \\\nonumber
& = \pi \tau'(t)\int_0^{\theta_q(t,y)}\Delta\phi(\tau(t),\theta_q)d\theta_q \\\nonumber
f(t,g(b(1; t(\tau))+2a(1; t(\tau))cos\theta)) & = \pi \tau'(t)\int_0^{\theta}\Delta\phi(\tau(t),\theta_q)d\theta_q\nonumber
\end{align}

So we can rewrite

\begin{align*}
\widehat{f\circ \mu^{-1}}_k^{(t)} & = \frac{1}{\pi}\int_0^\pi \pi \tau'(t)\int_0^{\theta}\Delta\phi(\tau(t),\theta_q)d\theta_qcos(k\theta)d\theta \\
& = \frac{\tau'(t)}{k}\int_0^\pi \Delta\phi(\tau(t),\theta)\sin(k\theta)d\theta\\
\int_Ie^{-i\tau(t)\omega}\widehat{f\circ \mu^{-1}}_k^{(t)}dt & = \frac{1}{k} \int_\mathbb{R}\int_0^\pi \Delta\phi(\tau,\theta)\sin(k\theta)e^{-i\tau\omega}d\theta d\tau
\end{align*}

Define an unitary transform $\mathcal{G}:L^2(\mathbb{R}\times (0,\pi))\to L^2(\mathbb{R})\times l^2(\mathbb{N})$ by $\mathcal{G}\phi(\omega,k)=\frac{1}{\pi}\int_\mathbb{R}\int_0^\pi \phi(\tau,\theta)e^{-i\omega\tau}\sin(k\theta)d\theta d\tau$.
Then we have $\mathcal{G}\Delta\phi(\omega,k)=-(\omega^2+k^2)\mathcal{G}\phi(\omega,k)$. By Plancherel's identity we conclude

\begin{align*}
\sigma^2(f\circ\mu^{-1}) &  =\pi\sum_{k=1}^\infty\int_{\mathbb{R}}\frac{1}{\omega^2+k^2}|\mathcal{G}\Delta\phi(\omega,k)|^2d\omega \\
& = -\pi \sum_{k=1}^\infty\int_{\mathbb{R}}\mathcal{G}\Delta\phi(\omega,k)\overline{\mathcal{G}\phi(\omega,k)}d\omega \\
& =\|\phi\|_{\nabla}^2
\end{align*}

From here we can conclude there is Gaussian free field in the pull-back of the liquid region.

§ SYMMETRY FOR THE COMPLEX STRUCTURE OF THE Q-HAHN LOZENGE TILING

There is a simpler argument to the complex structure if there is some symmetry in the model. It is the case for the regular equilateral hexagon partition with the q-Hahn weights. Intuitively, q-Hahn weights on type III lozenge on the $(x,t)$ plane is equivalent to counting the number of unite boxes in the $(x,y,z)$ 3-D space. And the counting of boxes is rotational invariant. Then the partition becomes an equilateral hexagon on the $(\tilde{t},\tilde{x})$ plane, which is an affine transformation given by $(t,x)=(\frac{2}{\sqrt{3}}\tilde{t}+1,\tilde{x}+1+\frac{1}{\sqrt{3}}\tilde{t})$ see figure(<ref>). Define $\tilde{\rho}_{I,II,III}(\tilde{t},\tilde{x})  \equiv \rho_{I,II,III}(\frac{2}{\sqrt{3}}\tilde{t}+1,\tilde{x}+1+\frac{1}{\sqrt{3}}\tilde{t})$. In this new coordinate, densities of different types at the same point are the same up to $3\pi/2$ rotation, by the rotation invariance. Therefore, the following identity holds,

\begin{align*}
\tilde{\rho}_{II}(\tilde{t},\tilde{x})   =  \tilde{\rho}_{III}(-\frac{1}{2}\tilde{t}+\frac{\sqrt{3}}{2}\tilde{x},-\frac{\sqrt{3}}{2}\tilde{t}-\frac{1}{2}\tilde{x})
\end{align*}

Change of Coordinates

For $t,x$ in the liquid region, we have

\begin{align}
\rho_{III}(t,x)& = \frac{1}{\pi} \arccos \frac{ q^{t+x}+q^{t+x+2}-q^{t+1}-q^{t+3}+q^{x+1}+q^{x+3}-2 q^{2 x+1}}{2 \sqrt{q \left(q^x-1\right) \left(q^x-q^2\right) \left(-q^{t+x}-q^{t+x+2}+q^{2 t+1}+q^{2 x+1}\right)}}\\
\rho_{II}(t,x) & = \frac{1}{\pi} \arccos \frac{-q^{t+x}-q^{t+x+2}+q^{t+1}+q^{t+3}+q^{x+1}+q^{x+3}-2 q^3}{2 \sqrt{q^2\left(q^t-1\right) \left(q^t-q^2\right) \left(q^x-1\right) \left(q^x-q^2\right)}}\\
\rho_{I}(t,x) & = \frac{1}{\pi} \arccos \frac{q^{t+x}+q^{t+x+2}+q^{t+1}+q^{t+3}-2 q^{2 t+1}-q^{x+1}-q^{x+3}}{2\sqrt{q \left(q^t-1\right) \left(q^t-q^2\right) \left(-q^{t+x}-q^{t+x+2}+q^{2 t+1}+q^{2 x+1}\right)}}
\end{align}

And since $\tan(\theta)=\frac{\sqrt{1-cos(\theta)^2}}{cos(\theta)}$, we have

\begin{align*}
\tan(\rho_{III}(t,x))& =  \frac{\sqrt{C}}{q^{t+x}+q^{t+x+2}-q^{t+1}-q^{t+3}+q^{x+1}+q^{x+3}-2 q^{2 x+1}}\\
\tan(\rho_{II}(t,x)) & = \frac{\sqrt{C}}{-q^{t+x}-q^{t+x+2}+q^{t+1}+q^{t+3}+q^{x+1}+q^{x+3}-2 q^3}\\
\tan(\rho_{I}(t,x)) & = \frac{\sqrt{C}}{q^{t+x}+q^{t+x+2}+q^{t+1}+q^{t+3}-2 q^{2 t+1}-q^{x+1}-q^{x+3}}\\
\text{where }  C & =-\left(q^2-1\right)^2 q^{2 t+2}+2 (q-1)^2 \left(q^3+q\right) \left(q^t+q\right) q^{t+x}\\
& +\left(-\left(q-q^3\right)^2-\left(q^2-1\right)^2 q^{2 t}+2 (q-1)^2 \left(q^3+q\right) q^t\right) q^{2 x}
\end{align*}

Then $\Omega(t,x)$ is the solution to

\begin{align}
\label{eq:complexstructure}
\left (\Omega-\frac{q^{t+x}+q^{t+x+2}+q^{t+1}+q^{t+3}-2 q^{2 t+1}-q^{x+1}-q^{x+3}}{-2 \left(q^x-q^{t+1}\right) \left(q^{x+1}-q^t\right)}\right)^2   =\frac{-C}{(-2 \left(q^x-q^{t+1}\right) \left(q^{x+1}-q^t\right))^2}
\end{align}

which is an algebraic equation to $q^x$ and $q^t$.  Also note that $\frac{q^{-2 (t+x+1)}}{\left(q^2-1\right)^2}*C=4A(t)C(t)-(\mu(x)+A(t)+C(t))^2$, where $A(t),C(t), v$ are defined as in (<ref>). Here we are to verify the complex Burgers equation. That is  for equation (<ref>) the following holds

\begin{align}
\frac{\partial_x\Omega}{\Omega}-\frac{\partial_t\Omega}{1-\Omega}=-\log q
\end{align}

This is equivalent to show that

\begin{align*}
\partial_x\arg \Omega+\partial_t\arg (1-\Omega) & =0\\
\partial_x\log| \Omega|+\partial_t\log|1-\Omega| & =-\log q
\end{align*}

And then it is not hard to verify the following with the formulas given above,

\begin{align*}
&\partial_x  \rho_I(t,x)-\partial_t\rho_{III}(t,x) =0\\
&\partial_x\log\left(\frac{Denominator(\tan (\rho_{I}))^2+C}{(Denominator(\tan (\rho_{I}))+Denominator(\tan (\rho_{III})))^2}\right)\\
& +\partial_t\log\left(\frac{Denominator(\tan (\rho_{III}))^2+C}{(Denominator(\tan (\rho_{I})^2+Denominator(\tan (\rho_{III})))^2}\right)  =-2\log q
\end{align*}

And this justifies the Burgers equation for $\kappa=0$.

§.§ q-Racah Weighted Tilings and Non-Intersecting Paths

A tiling means a covering of the partition with lozenges inside the partition adjacent to each other. And a typical configuration of such is the following picture

Example of a tiling configuration and non-intersecting paths

As the left of the picture indicates each of the tiling can also be viewed as boxes piled in a corner of a room. We will consider a particular probability measure on the set of all possible tilings of the hexagon associated with the "volumes of boxes". To introduce this probability measure, we need to recall that each tiling can be represented by a collection of non-intersecting paths on a weighted directed graph.

Path assignment

Paths are assigned according to types of lozenges. For Type I lozenges we draw a north-east path through it and  through Type II an east moving path. There are no paths in Type III lozenges, which are usually called holes (see Figure (<ref>) ). It is well known that for each tiling, such arrangement gives rise to a collection of north-east moving non-intersecting paths starting at $(0,x(0)+\frac{1}{2})$ where $x(0)=0,1,\dots, \mathbf{a}-1$ and terminating at $(\mathbf{b}+\mathbf{c},x(\mathbf{b}+\mathbf{c})+\frac{1}{2})$ where  $x(\mathbf{b}+\mathbf{c})=\mathbf{c},\mathbf{c}+1,\dots, \mathbf{c}+\mathbf{a}-1$. Shift all paths $\frac{1}{2}$ downwards in the $x$ direction and then we have the following picture (Figure <ref>). The weighted directed graph $G=(V,E)$ has vertices and edges

\begin{align*}
V & = \{(t,x)\in\mathbb{Z}^2: 0\leq t\leq \mathbf{b}+\mathbf{c}, 0\leq x\leq \mathbf{a}+\mathbf{c}-1, -\mathbf{b}\leq x-t\leq \mathbf{a}-1 \} \\
E & = \{(t,x)\to (t+1,x), (t,x)\to (t+1,x+1): (t,x)\in V\}
\end{align*}

Let $\mathbf{q}$ be either positive real or on the unit circle. The weights on the edges $(t,x)\to (t+1,x+1)$ are defined as

\begin{align}
w_1(t,x) & = \mathbf{q}^{\frac{\mathbf{a}+\mathbf{b}+\mathbf{c}-t-1}{2}}\frac{1-\kappa^2\mathbf{q}^{x-\mathbf{b}-\mathbf{c}+1}}{1-\kappa^2\mathbf{q}^{2x-t-\mathbf{c}+1}} (1-\mathbf{q}^{-\mathbf{a}-\mathbf{c}+x+1}) \label{eq:weight_w1}
\end{align}

The weights on the edges $(t,x)\to (t+1,x)$ are defined as

\begin{align}
w_0(t,x) & = \mathbf{q}^{\frac{t+1-\mathbf{a}-\mathbf{b}-\mathbf{c}}{2}} \frac{1-\kappa^2\mathbf{q}^{x+\mathbf{a}-t}}{1-\kappa^2\mathbf{q}^{2x-t-\mathbf{c}+1}} (\mathbf{q}^{x-t+\mathbf{b}}-1). \label{eq:weight_w2}
\end{align}

These weights have three factors each that can be interpreted as the q-volume part, Racah part and the boundary restriction part. The last of which means that the paths are restricted into the hexagon drawn above without any specification in the vertices in the north-east phase of the $(t,x)$ plane.

We further define the weight of a path by the product of all weights of the edges and the weight of a configuration of paths as the product of all weights of each path

\begin{align*}
w(\pi_0,\dots, \pi_{\mathbf{a}-1})=\prod_{j=0}^{\mathbf{a}-1}\prod_{e_j\in\pi_j}w(e_j)
\end{align*}

agraphThe Probability Measure of a collection of paths $\{\pi_j: j=0,\dots, \mathbf{a}-1\}$ is proportional to its weight. And we consider the conditional measure conditioned on the fact that no path collides to any other in the partition. According to the Lindström–Gessel–Viennot theorem (LGV), this implies that the probability of the non-intersecting paths $\{\pi_j\}_{j=0}^{\mathbf{a}-1}$ passing through $\{(t,x_j(t))\}_{j=0,t=0}^{\mathbf{a}-1,\mathbf{b}+\mathbf{c}}$ is

\begin{align}
\mathbb{P}(\pi_j(t)=x_j(t),j=0,\dots,\mathbf{a}-1, t=0,\dots,\mathbf{b}+\mathbf{c}| \{\pi_j\} \text{ non -intersecting})
& = \frac{1}{Z}\prod_{t=0}^{\mathbf{b}+\mathbf{c}}\det(T_t(x_{i}(t),x_{j}(t+1)))_{i,j=0}^{\mathbf{a}-1}
\end{align}

where $Z$ is the normalizing factor and the transition weight $T_t$ is defined as

\begin{align*}
T_t(x,y)& =\begin{cases}w_1(t,x) \quad  y-x=1 \\ w_0(t,x) \quad y=x \\0 \quad otherwise\end{cases}
\end{align*}

[7] lemma 7.4. using the q-difference relation and the orthogonality of q-Racah polynomials showed that

\begin{align*}
U_t(x,y) & = \frac{1}{w^{(t+1)}(y)}T_{t}(x,y)\\ & = \frac{w_1(t,x)}{w^{(t+1)}(x+1)}\delta_{y=x+1}+\frac{w_0(t,x)}{w^{(t+1)}(x)}\delta_{y=x} \\
& =\sum_{k=0}^M\frac{c_{k,t+1}}{c_{k,t}}r_{k,t}(x)r_{k,t+1}(y) \\
\end{align*}

where $ r_{k,t}(x) =const_4\phi_3\bigg(\begin{matrix}\mathbf{q}^{-k}, \alpha\beta \mathbf{q}^{k+1},\mathbf{q}^{-x},\delta \gamma \mathbf{q}^{x+1}\\\alpha \mathbf{q},\delta \beta \mathbf{q},\gamma \mathbf{q}\end{matrix}; \mathbf{q},\mathbf{q} \bigg)$ is the q-Racah orthonormal polynomial with respect to the q-Racah weight $w^{(t)}(x)$ defined in (<ref>). Note that here we denote $w^{(t)}\equiv w$ to emphasise that the parameters are depending on $t$ as shown in Table <ref>. $M=-log_\mathbf{q}\gamma$ denotes the highest order, and $c_{k,t}$ is defined as

\begin{align*}
c_{k,0} &  = \sqrt{\prod_{s=0}^k\frac{1-\mathbf{q}^{-(\mathbf{a}+\mathbf{b}+\mathbf{c}-s)}}{1-\mathbf{q}^{-(\mathbf{a}-s)}}}\\
c_{k,t} & = c_{k,0}\sqrt{\prod_{s=0}^{t-1}(1-\mathbf{q}^{-\mathbf{a}-s+k})(\mathbf{q}^{(\mathbf{a}+\mathbf{b}+\mathbf{c}-s-k-1)}-1)} .
\end{align*}

One can also think of $r_{k,t}(x)$ as a polynomial in $\nu(x)$ of order $k$,

\begin{align}
\nu(x)= \mathbf{q}^{-x}+\delta\gamma \mathbf{q}^{x+1}
\end{align}

According to [7] there are four cases to consider (see Table <ref>). If one assumes $b<c$, only i,iii, iv cases are involved, whereas if $b>c$, only i ii iv cases.  However, there is some symmetry of the parameter specification, i.e., $\alpha\beta=q^{-2a-b-c}$ and $\delta=\kappa^2q^{2a}\alpha$ for all four cases, $\delta\gamma$ equals to $\kappa^2q^{-t/n-c}$ in case i and ii (or $\kappa^2q^{t/n-b-c}$ in case iii and iv), and $\alpha$ in case i (or case iii) equals to $\gamma$ in case ii (or case iv).  Then the case i and ii give the same q-Racah polynomial ensemble and so do case iii and iv. Indeed, the weights of the lozenges only depend on their vertical position. And the hexagon domain can be translated to the one with the lower boundary on the line $x=0$ and one can observe a symmetry along $t=\frac{\mathbf{b}+\mathbf{c}}{2}$. Actually, one can verify that the in case i and ii $\sqrt{4A(\xi)C(\xi)-(\nu(x)-1-\gamma\delta \mathbf{q} +A+(\xi)C(\xi))^2}$ is $\mathbf{q}^{\mathbf{b}-t}$ times of that in case iii and iv. Therefore, we can reduce four cases into one and focus on case i where $(t,x)$ can be extended to the entire hexagon domain  in the computation.

To establish a probability measure one also needs a configuration of $(q, \kappa)$ so as to $w(\pi_0,\dots, \pi_{\mathbf{a}-1})>0$. Typically there are three cases one can consider according to [7]

* Real case: $q>0, \kappa\in \mathbb{R}$ and $\kappa$ can not lie in $[q^{-a+\frac{1}{2n}},q^{\frac{b+c}{2}-\frac{1}{2n}}]$ if $q>1$ or $[q^{\frac{b+c}{2}-\frac{1}{2n}}, q^{-a+\frac{1}{2n}}]$ if $q<1$.

* Imaginary case : $q>0, \kappa\in i \mathbb{R}$.

* Trigonometric case : $q=e^{i \lambda}, \kappa=e^{i \eta}$. The normalising factor in (<ref>) and (<ref>) is to make sure, after the common factor $2i$ being removed, the weights are real numbers as,

\begin{align}
w_1(t,x) = & \frac{\sin\left(\eta+\frac{\lambda}{2}(\frac{x+1}{n}-b-c)\right)\sin\left(\frac{a+c-\frac{x}{n}}{2}\right)}{\sin \left(\eta +\frac{\lambda}{2}(2x-c-\frac{t-1}{n})\right)}\\
w_0(t,x) = & \frac{\sin\left(\eta+\frac{\lambda}{2}(\frac{x-t}{n}+a)\right)\sin\left(\frac{x-t}{2n}+\frac{b}{2}\right)}{\sin \left(\eta +\frac{\lambda}{2}(\frac{2x-t+1}{n}-c)\right)}&
\end{align}

To make all weights above have the same sign for all $(t,x)$ in the $abc$ hexagon partition,  it is enough to let $\lambda,\eta$ be such that $k\pi \leq \eta-\lambda(b+c)/2, \eta+\lambda a\leq (k+1)\pi$ for some $k$ to be an integer. We call it a singular trigonometric case, whenever one of the inequalities binds,

* Singular real case:   $q>0$ and $\kappa = q^{-a}$ or $q^{(b+c)/2}$.

The relation of q-Racah weights of paths and the weight of horizontal lozenges are connected as the followings

\begin{align*}
\frac{w_1(t,x)w_0(t+1,x+1)}{w_1(t+1,x)w_0(t,x)}=\frac{\partial_x\nu(x-1; t-1)}{\partial_x\nu(x; t-1)}
\end{align*}