\section{Introduction} 
The symmetric and asymmetric simple exclusion processes have received a lot of attention in the last decades. One of the reasons is that they are exactly solvable and thus exhibit a profound underlying mathematical structure. A particularly interesting setup that preserves the integrability is when the SSEP/ASEP is in contact with two reservoirs at the left and right end of the chain. In this setup the model is still integrable but more difficult to study, e.g. the steady state does not factorise any longer. The processes are known to be related to the XXX and the XXZ Heisenberg spin chain and thus can be formulated in the framework of the quantum inverse scattering method \cite{Faddeev:1996iy} that was generalised  in \cite{Sklyanin:1988yz} to the case of open boundaries which are relevant to model the reservoirs of the processes. 

More recently an integrable generalisation of the ASEP was proposed that is known as the q-Hahn process. In a certain regime of parameters this model has already appeared two decades ago \cite{Sasamoto} but the full jump probabilities  were only recently obtained by Povolotsky \cite{Povolotsky}. The rates and the extension to right and left particle jumps was then given in \cite{barraquand}.
  Similar as in the case of the ordinary ASEP, the generator of the q-Hahn process is related to the  Hamiltonian of the $U_q(sl_2)$ invariant non-compact  XXZ spin chain with a certain choice of parameters by a similarity transformation as discussed in \cite{Frassek:2019isa}. The relevant representations on the infinite-dimensional vector space $V$ whose tensor product forms the  quantum space of the spin chain is defined by the following action of the generators
\begin{equation}\label{eq:sl2act}
 S_+|m\rangle=[m+2s]|m+1\rangle\,,\quad S_-|m\rangle=[m]|m-1\rangle\,,\quad S_0|m\rangle=(m+s)|m\rangle\,.
\end{equation}  
Here the spin label is taken to be $s>0$. The  generators obey the standard $U_q(sl_2)$ commutation relations
\begin{equation}\label{eq:sl2}
 [S_+,S_-]=-[2S_0]\,,\quad\quad [S_0,S_\pm]=\pm S_\pm\,,
\end{equation} 
with the q-number defined via
\begin{equation}
 [x]=\frac{q^{x}-q^{-x}}{q-q^{-1}}\,.
\end{equation}  
The algebraic formulation is advantageous in particular when studying stochastic dualities of the processes see e.g. \cite{Schtz1994NonAbelianSO,giardina2009duality} but not to describe the actual process which for a given representation arises in a specific choice of basis.

The main goal of this note is to introduce boundary reservoirs for the continuous time q-Hahn process that preserve the integrable structure. This has been done in \cite{Frassek:2019vjt} for the rational $q\to1$ limit following the  quantum inverse scattering method. The authors derived a suitable K-matrix written in terms of $sl_2$ generators (similar to the evaluation map in the case of the ordinary Yangian) and were able to show that it yields the desired boundary terms following Sklyanin's prescription \cite{Sklyanin:1988yz} by computing the relevant matrix elements of the resulting boundary Hamiltonian. It should be possible to apply the same strategy for the q-deformed case but we are currently lacking the algebraic expression of the most general K-matrix, see \cite{Tsuboi:2019qac} for some partial results. To circumvent the construction of the K-matrix we follow an idea initially proposed by Sutherland  \cite{sutherland1970two} for the Hamiltonian density and generalise it to our setup with boundaries. 
For this purpose  it is convenient to fix the basis and work right from the beginning in components. We thus follow closely the conventions of  \cite{Mangazeev:2019rzf}  to avoid algebraic expressions.
 
  
  The note is organised as follows: Section~\ref{sec:qhahn} summarises our results and presents  the Markov generator of  the open q-Hahn process. Section~\ref{sec:manga} reviews the formalism and results of \cite{Mangazeev:2019rzf} that are needed in order to argue that the proposed boundary conditions are integrable. This latter point is discussed in Section~\ref{sec:ham}. The argument is based on solving the boundary Yang-Baxter equation and its first derivative with respect to the spectral parameter at special points that are relevant in the construction of the spin chain Hamiltonian. Section~\ref{sec:otasep} briefly discusses the open q-TASEP limit and Section~\ref{sec:conc} contains our conclusions. Further details of computations can be found in the appendix.
 
 


\section{The open q-Hahn process}\label{sec:qhahn}
In this section we propose integrable boundary conditions for the q-Hahn process.  The full stochastic Hamiltonian of the open non-compact XXZ spin chain with $N$ sites and reservoirs at site $1$ and $N$ is of the form
\begin{equation}\label{eq:ham}
H=B_L+\sum_{i=1}^{N-1}\mathcal{H}_{i,i+1}+B_R
\end{equation} 
where 
\begin{equation}\label{eq:qspace} 
 H: V\otimes \ldots \otimes V \rightarrow V\otimes \ldots \otimes V\,.
\end{equation} 
It relates to the continuous-time Markov generator $M$ with standard conventions via the simple relation
\begin{equation}
 M=-H^{t}\,.
\end{equation} 
where $t$ denotes the transposition in the quantum space $V\otimes \ldots\otimes V$.
The Hamiltonian \eqref{eq:ham} is stochastic since all matrix elements outside
the diagonal are negative or zero. The sum over the columns vanishes. 
The bulk part has been introduced in \cite{Sasamoto,Povolotsky,barraquand}, see also \cite{Frassek:2019isa}, and is discussed in Section~\ref{sec:bulkH}. The boundary part is given in Section~\ref{sec:bnd}. Each boundary term introduced  at site $1$ and $N$ is governed by one additional  parameter $\rho_L$ and $\rho_R$ respectively.
\subsection{Bulk}\label{sec:bulkH}
The  q-Hahn process in the bulk is obtained from the stochastic Hamiltonian density defined 
by  the action on two neighboring sites with  $m$ particles on one site and $m'$ particles on the other via
\begin{equation}\label{eq:hact}
\begin{split}
\mathcal{H}|m\rangle \otimes |m'\rangle=\left(\alpha_-(m)+\alpha_+(m')\right)|m\rangle \otimes |m'\rangle&-\sum_{k=1}^m \beta_-(m,k) |m-k\rangle \otimes |m'+k\rangle\\&-\sum_{k=1}^{m'} \beta_+(m',k) |m+k\rangle \otimes |m'-k\rangle\,.
 \end{split}
\end{equation} 
The rates for particles moving to the left are of the form
\begin{equation}\label{eq:bp}
\beta_+(m,k)= 
\frac{\mu ^{k} }{1-\gamma^k}\frac{(\gamma;\gamma)_m (\mu ;\gamma)_{m-k}}{ (\gamma;\gamma)_{m-k} (\mu ;\gamma)_m}\,,
\end{equation} 
while the rates for particles that move to the right read
\begin{equation}\label{eq:bm}
\beta_-(m,k)= 
\frac{1}{1-\gamma^k}\frac{(\gamma;\gamma)_m (\mu ;\gamma)_{m-k}}{ (\gamma;\gamma)_{m-k} (\mu ;\gamma)_m}\,.
\end{equation}  
Here the parameters take values $0<\gamma,\mu<1$ and  $(a;\gamma)_m=\prod_{j=0}^{m-1}(1-a\gamma^j)$ denotes
the q-Pochhammer symbol.
The diagonal terms are expressed in terms of finite sums as
\begin{equation}\label{eq:alphas}
\alpha_+(m)
=\sum _{k=0}^{m-1} \frac{\gamma^{k}}{\mu^{-1}-\gamma^{k}}\,,\qquad 
 \alpha_-(m)
 =\sum _{k=0}^{m-1} \frac{1}{1-\gamma^{k}\mu}\,,
\end{equation} 
see also  Appendix~\ref{app:sumrates}.
For later purposes we set 
\begin{equation}
 \gamma=q^2\,,\qquad \mu=q^{4s}\,.
\end{equation}  
\subsection{Boundary}\label{sec:bnd}
Taking inspiration from \cite{Frassek:2019vjt}, we introduce the boundary terms 
\begin{equation}\label{eq:bd1}
\begin{split}
B_L|m\rangle =\left(\alpha_+(m) +\sum_{k=1}^\infty\frac{\rho_L^k}{1-\gamma^{k}}\right)|m\rangle&
-\sum_{k=1}^m \beta_+(m,k)|m-k\rangle -\sum_{k=1}^\infty\frac{\rho_L^k}{1-\gamma^{k}}|m+k\rangle \,,
\end{split}
\end{equation} 
and
\begin{equation}\label{eq:bd2}
\begin{split}
B_R|m\rangle =\left(\alpha_-(m) +\sum_{k=1}^\infty\frac{\rho_R^k}{1-\gamma^{k}}\right)|m\rangle&
-\sum_{k=1}^m\beta_-(m,k)|m-k\rangle -\sum_{k=1}^\infty\frac{\rho_R^k}{1-\gamma^{k}}|m+k\rangle \,.
\end{split}
\end{equation} 
Particles are removed from the first and last site with the  bulk rates into the reservoirs. The insertion rates are governed by one real parameter  at each boundary $\rho_L$ and $\rho_R$ with $0<\rho_L,\rho_R<1$ as in the rational case.
We are free to shift $\rho_R\to\rho_R\mu$ in order to achieve an expression more similar to the bulk  rates in \eqref{eq:bp} and \eqref{eq:bm}.



\section{Quantum inverse scattering method}\label{sec:manga}

In the framework of the quantum inverse scattering method the Hamiltonian is derived from the fundamental transfer matrix. More precisely it is obtained from the logarithmic derivative at the permutation point $x_0$, i.e.
\begin{equation}\label{eq:logderiv1}
H=-\frac{1}{4}\frac{\partial}{\partial x}\ln T(x)\big|_{x=x_0} +c_0\,I \,,
\end{equation} 
with a constant $c_0$ that  is proportional to the identity and which we will fix later on.
 

The transfer matrix is given by Sklyanin's double-row construction \cite{Sklyanin:1988yz}   and is expressed in terms of R- and K-matrices  via
\begin{equation}\label{eq:transm}
 T(x)=\tr_a \bar{\mathcal{K}}_a(x)\mathcal{R}_{a,1}(x)\cdots \mathcal{R}_{a,N}(x)\mathcal{K}_a(x)\mathcal{R}_{N,a}(x^{-1})\cdots \mathcal{R}_{1,a}(x^{-1})\,.
\end{equation} 
The R-matrices act non-trivially on the tensor product of two vector spaces $V\otimes V$ as introduced above   in \eqref{eq:sl2}  while the K-matrices act non-trivially on a single copy of the auxiliary space $V_a\sim V$ isomorphic to a single site of the quantum space. This   transfer matrix  contains the Hamiltonian. 


Additionally we introduce the transfer matrix with the two-dimensional defining representation in the auxiliary space 
\begin{equation}\label{eq:bat}
 T_\square(x)=\tr_\square  \bar{K}_\square (x)\mathcal{L}_{\square,1}(x)\cdots \mathcal{L}_{\square ,N}(x){K}_\square (x)\mathcal{L}_{N,\square }(x^{-1})\cdots \mathcal{L}_{1,\square }(x^{-1})\,,
\end{equation} 
here indicated with $\square$.
The Lax matrices $\mathcal{L}$ act non trivially on $\mathbb{C}^2$ in the auxiliary space and the quantum space at a single site $V$. The K-matrices act non-trivially on the auxiliary space and are matrices of the size $2\times2$. This transfer matrix is usually used in the context of the algebraic Bethe ansatz. 

To ensure commutativity of the transfer matrix the boundary Yang-Baxter relations 
\begin{equation}\label{eq:bybe1}
 \mathcal{L}_{a,\square }(x/y)\mathcal{K}_a(x) \mathcal{L}_{\square,a}(xy)K_\square(y)=K_\square(y) \mathcal{L}_{a,\square}(xy)\mathcal{K}_a(x)  \mathcal{L}_{\square,a}(x/y)
\end{equation} 
and
\begin{equation}\label{eq:bybe2}
  \mathcal{L}_{a,\square}(y/x)\bar{\mathcal{ K}}_a(x) \bar{ \mathcal{L}}_{\square,a}(1/(xy))\bar K_\square(y)=\bar K_\square(y)\bar{ \mathcal{L}}_{a,\square}(1/(xy))\bar {\mathcal{K}}_a(x) { \mathcal{L}}_{\square,a}(y/x)
\end{equation} 
for the left and right boundary have to be obeyed. For a detailed discussion that also takes into account the bulk part we refer the reader to \cite{Mangazeev:2019rzf}. 

The Lax matrices $\bar{\mathcal{L}}$ in the second boundary Yang-Baxter relation \eqref{eq:bybe2}  is defined via
\begin{equation}\label{eq:Lbars}
 \bar{ \mathcal{L}}_{\square, a}(x)=\left(\left(\mathcal{L}_{a,\square}^{t_\square}(1/x)\right)^{-1}\right)^{t_\square}\,,\qquad \bar{ \mathcal{L}}_{a,\square}(x)=\left(\left(\mathcal{L}_{\square ,a}^{t_a}(1/x)\right)^{-1}\right)^{t_a} \,.
\end{equation} 
We will express them in terms of the Lax matrices $\mathcal{L}$ in Section~\ref{sec:cros} using the appropriate crossing relations.
The relevant R- and Lax matrices  can be found in the component form in \cite{Mangazeev:2019rzf}.  
The Lax matrices read
\begin{equation}
\langle m|\mathcal{L}_{\square ,a}(x)|m'\rangle=\frac{1}{ \left[x q^{\frac{1}{2}-s}\right]}\left(\begin{array}{cc}
                      \delta_{m,m'}\left[x q^{\frac{1}{2}-s-m}\right]q^{m}&\delta_{m,m'+1}x\left[q^{1-2s-m}\right] q^{m-\frac{1}{2}-s}\\\delta_{m,m'-1}
                      x^{-1}\left[q^{m+1}\right] q^{\frac{1}{2}+s+m}&\delta_{m,m'}\left[x q^{m+s+\frac{1}{2}}\right]q^{m+2s}
                     \end{array}
\right)
\end{equation} 
and
\begin{equation}
\langle m|\mathcal{L}_{a,\square}(x)|m'\rangle=\frac{1}{ \left[x q^{\frac{1}{2}-s}\right]}\left(\begin{array}{cc}
                      \delta_{m,m'}\left[x q^{\frac{1}{2}-s-m}\right]q^{-m}&\delta_{m,m'+1}x^{-1}\left[q^{1-2s-m}\right] q^{-m-s+\frac{1}{2}}\\\delta_{m,m'-1}
                      x\left[q^{m+1}\right] q^{-\frac{1}{2}-s-m}&\delta_{m,m'}\left[x q^{m+s+\frac{1}{2}}\right]q^{-m-2s}
                     \end{array}
\right)\,.
\end{equation} 
They obey the unitarity relations $\mathcal{L}_{a,\square}(x)
 \mathcal{L}_{\square,a}(x^{-1})=I$. 
The R-matrix for the fundamental transfer matrix \eqref{eq:transm} is expressed through its components  in the same article \cite{Mangazeev:2019rzf}. It is written in terms of a generalised hypergeometric function as
\begin{equation}\label{eq:Rmatrix1}
\begin{split}
 \langle m,n|\mathcal{R}(x)|m',n'\rangle =\delta_{m+n,m'+n'}q^{4ms} &
 \left[\begin{array}{c}
 m+n                                                                      \\
m                                                                           \end{array}
\right]_{q^2}
\frac{\left(x^{-2};q^{2}\right)_{n'}\left(x^{-2};q^2\right)_{m}\left(q^{4s};q^2\right)_{n}}{\left(x^{-2}q^{4s};q^2\right)_{m+n}\left(q^{4s};q^2\right)_{n'}}\\
&\hspace{-1.15cm}\times\,\,\,\pFq{4}{3}{ q^{-2m},q^{-2n'},x^2 q^{4s},x^2 q^{2-4s-2m-2n}}{q^{-2m-2n},x q^{2-2m},x q^{2-2n'}}{q^2,q^2}\,.
\end{split}
\end{equation} 
where the q-binomial is defined by
\begin{equation}
  \left[\begin{array}{c}
 n                                                                      \\
m                                                                           \end{array}
\right]_{q}=\frac{(q;q)_n}{(q;q)_{n-m}(q;q)_m}\,.
\end{equation} 
The  R-matrix is not symmetric under the exchange $m\leftrightarrow n$ and $m'\leftrightarrow n'$. In the definition of the transfer matrix \eqref{eq:transm} we  use the notation
\begin{equation}
 \mathcal{R}(x)= \mathcal{R}_{a,i}(x)\,,
\end{equation} 
and 
\begin{equation}\label{eq:symm}
 \mathcal{R}_{i,a}(x)=\mathcal{P}\mathcal{R}_{a,i}(x)\mathcal{P}=\mathcal{R}^{-1}(x^{-1})\,,
\end{equation} 
where $i=1,\ldots, N$ denote the spin chain sites and $\mathcal{P}$ the permutation operator. 
We further remark that there is a diagonal Drinfeld twist involved, cf.~\cite{borodin2016stochastic}. This is not immediately relevant for us but convenient as the logarithmic derivative will immediately yield a Hamiltonian on the stochastic line. The twist can also be introduced  at a later stage and appears in the form of a similarity transformation on the level of the Hamiltonian, cf.~\cite{deGier:2005zz,Frassek:2019isa}.


For later purposes we also introduce an  alternative expression  for the R-matrix. 
In  \cite{Mangazeev:2019rzf} it has been written as
\begin{equation}\label{eq:Rmatrix2}
\begin{split}
 \langle m,n|\mathcal{R}_{a,i}(x)|m',n'\rangle =\delta_{m+n,m'+n'} \sum_{k=0}^{m}\Phi_{q^2}(k|k+n;x^{-2},x^{-2}q^{4s})\Phi_{q^2}(m-k|n';x^2q^{4s},q^{4s})
\end{split}
\end{equation} 
with 
\begin{equation}
 \Phi_{q}(m|n;x,y)=\left(\frac{y}{x}\right)^m
\frac{\left(x;q\right)_{m}\left(y/x;q\right)_{n-m} }{\left(y;q\right)_{n }}
 \left[\begin{array}{c}
 n                                                              \\
m                                                                    \end{array}
\right]_{q}\,.
\end{equation} 
This factorisation reminds of Derkachov's factorisation in relation to Q-operators \cite{Derkachov:1999pz}. We will see later in Section~\ref{sec:bulkT} that the product of the two terms  turns into a sum at the level of the Hamiltonian density. One term yields the left moving particles and the other the right moving particles which combine into the q-Hahn process.  

Finally we also introduce the K-matrices in the two-dimensional defining representation that appears in the definition of the transfer matrix \eqref{eq:bat}. It depends on the four parameters $t^\pm_R,\nu_R,\kappa_R$ that will be adjusted later and reads
\begin{equation}
 K_\square(x)=\left(\begin{array}{cc}
                     \frac{t^-_R}{q\nu_R}-q\nu_R t^+_R+x^2(t^-_R-t^+_R)&\kappa_R^{-1}t_R^-(x^2-x^{-2})\\
                     \kappa_R t_R^+(x^2-x^{-2})& \frac{t_R^-}{q\nu_R}-q\nu_R t_R^++x^{-2}(t_R^--t_R^+)
                    \end{array}
\right)\,.
\end{equation} 
As discussed in   \cite{Mangazeev:2019rzf}  we can express the other K-matrix via
\begin{equation}\label{eq:crossK}
 \bar K_\square(x)=\left.D_\square^{-1}K_\square(1/(qx))\right|_{R\to L}\,,
\end{equation} 
where we introduced another four variables $t^\pm_L,\nu_L,\kappa_L$ and the matrix 
\begin{equation}
 D_\square=\left(\begin{array}{cc}
                                                         1&0\\
                                                         0&q^2
                                                        \end{array}
\right) \,.
\end{equation} 
To our knowledge, the K-matrices $\mathcal{K}$ and $\bar{\mathcal{K}}$ in \eqref{eq:bybe1} and \eqref{eq:bybe2} seem currently not available in the literature in the infinite-dimensional representations of interest.    
 

\subsection{Crossing relations and the other boundary}\label{sec:cros}
In this subsection we discuss the crossing relations for the Lax matrices introduced above and argue that the boundary Yang-Baxter equation \eqref{eq:bybe2} can be brought to the form \eqref{eq:bybe1}. From the latter statement it then follows that the left and right K-matrices are determined by the same equation \eqref{eq:bybe1}. 
The crossing relations have been given in \cite{Mangazeev:2019rzf}. The first one reads
\begin{equation} 
\left(
 \mathcal{L}_{s,\square}^{t_\square}(1/ x)\right)^{-1}=g(x) D_\square  \mathcal{L}_{\square,s}^{t_\square}(x q^{-2})D_\square^{-1}
\end{equation} 
with 
\begin{equation}
 g(x)=\frac{\left(q^{s+3}-x^2\right) \left(x^2 q^s-q\right)}{\left(q^{s+1}-x^2\right) \left(x^2 q^s-q^3\right)}\,.
\end{equation}  
It can be used  to express  the Lax matrix appearing in \eqref{eq:bybe2} as
\begin{equation}
 \bar{\mathcal{L}}_{\square,s}(x)=g(x)D_\square^{-1}\mathcal{L}_{\square,s} (x q^{-2})D_\square\,.
\end{equation}  
The second crossing relation relevant is of the form
\begin{equation} 
\left(
 \mathcal{L}_{\square,s}^{t_s}(1/x)\right)^{-1}=g(x) D_s\mathcal{L}_{s,\square}^{t_s}(x q^{-2})D_s^{-1}
\end{equation} 
with the same function $g(x)$. The similarity transformation is diagonal and  of the form
\begin{equation}\label{eq:dmat}
 \langle m|D_s|m'\rangle =q^{2m}\delta_{m,m'} \,.
\end{equation} 
Similar as before we obtain
\begin{equation}
\bar{ \mathcal{L}}_{s,\square}(x)=g(x)D_s^{-1}S_{s,\square}(x q^{-2})D_s\,,
\end{equation} 
cf.~\eqref{eq:bybe2}. 
We insert the obtained relations into the boundary Yang-Baxter equation \eqref{eq:bybe2}. Further,  inverting the K-matrix $K_\square$ and $D_\square$ and shifting the spectral parameters $y\to y q^{-1}$ and $x\to (xq)^{-1}$ we obtain
\begin{equation}
   \bar K_\square^{-1}(y/q)D_\square^{-1} S_{s,\square}(y x)D_s\bar K_s(1/(xq)) {S}_{\square,s}(x/y)={S}_{s,\square}(x/y)D_s\bar K_s(1/(xq)) S_{\square,s}(xy)\bar K_\square^{-1}(y/q)D_\square^{-1}\,.
\end{equation} 
Here we also used the invariance condition $
 [D_s D_\square,S_{s,\square}(y x)]=0$.
As  the K-matrix obeys the relation \eqref{eq:crossK}, it follows that the inverse of $\bar K_\square$ can be expressed as
\begin{equation}
 \bar K_\square^{-1}(y/q)D_\square^{-1} =\left.K_\square^{-1} (1/y)\right|_{R\to L}\,.
\end{equation} 
Further it is easy to verify that 
\begin{equation}
 K_\square^{-1} (1/y)=\frac{q^2 \nu_R^2 y^4}{\left(q \nu_R+y^2\right) \left(q \nu_R y^2+1\right) \left(q t^+_R\nu_R-t^-_R y^2\right) \left(q t^+_R\nu_R y^2-t^-_R\right)} K_\square (y)\,.
\end{equation} 
Inserting also this relation in the second boundary Yang-Baxter equation we find that solutions to the two equations can be related via
\begin{equation}
 \bar{\mathcal{K}}_s(x)=\left.D_s^{-1}\mathcal{K}_s(1/(qx))\right|_{R\to L}\,.
\end{equation} 
This relation will be of importance in Section~\ref{sec:ham}.

\subsection{Component form of BYBE}
We end this  section with the boundary Yang-Baxter equation \eqref{eq:bybe1} in component form as given in \cite{Mangazeev:2019rzf}. The equations below can be obtained taking coefficients in spectral parameter $x$:
\begin{equation}\label{manga1}
 \begin{split}
 \kappa_R t^+_R q^{2-4s}\left(1-q^{2(j+2s)}\right)\mathcal{K}_{j,l+1}(y)
 +\kappa^{-1}_Rt^-_Rq^{-4s}\left(1-q^{2+2l}\right)\mathcal{K}_{j+1,l}(y)\\
 +\nu^{-1}_Ry^2q^{2-2s}\left(q^{2j}-q^{2l}\right)\left(t^-_R-\nu^2_Rq^2t^+_R\right)\mathcal{K}_{j+1,l+1}(y)
 -\kappa_R t^+_R y^4q^{2-4s}\left(1-q^{2(1+l+2s)}\right)\mathcal{K}_{j+1,l+2}(y)\\
 -\kappa^{-1}_Rt^-_Ry^4 q^{-4s}\left(1-q^{4+2j}\right)\mathcal{K}_{j+2,l+1}(y)=0\,,
 \end{split}
\end{equation} 
and
\begin{equation}\label{manga2}
 \begin{split}
 \kappa_R t^+_R q^{2(2+l-2s)}\left(1-q^{2(j+2s)}\right)\mathcal{K}_{j,l+1}(y)
 +\kappa^{-1}_Rt^-_Rq^{2(1+j-2s)}\left(1-q^{2+2l}\right)\mathcal{K}_{j+1,l}(y)\\
 +q^{2-4s}\left(q^{2j}-q^{2l}\right)\left(t^+_R-t^-_R\right)\mathcal{K}_{j+1,l+1}(y)
 -\kappa_R t^+_R q^{2(1+j-2s)}\left(1-q^{2(1+l+2s)}\right)\mathcal{K}_{j+1,l+2}(y)\\
 -\kappa^{-1}_Rt^-_R q^{-4s+2l}\left(1-q^{4+2j}\right)\mathcal{K}_{j+2,l+1}(y)=0\,.
 \end{split}
\end{equation} 
Here we use the notation $\mathcal{K}_{i,j}(x)=\langle i|\mathcal{K}(x)|j\rangle $. The equations above yield conditions on five neighboring entries of the K-matrix.  This is discussed in further detail in Section~\ref{sec:rib}.



\section{Stochastic Hamiltonian from the transfer matrix}\label{sec:ham}



Let us compute the logarithmic derivative of the fundamental transfer matrix \eqref{eq:transm}.   
Assuming that the K-matrix at $x=1$ is proportional to the identity $\mathcal{K}(1)\proptoI$ we find 
\begin{equation}\label{eq:logderiv2}
-\frac{1}{4}\frac{\partial}{\partial x}\ln T(x)\big|_{x=1}=-\frac{1}{4}\frac{\tr_a \bar{\mathcal{ K}}_a'(1)}{\tr_a \bar{\mathcal{K}}_a(1)}-\frac{1}{2}\frac{\tr_a \bar{\mathcal{K}}_a(1)\mathcal{R}'_{a,1}(1)\mathcal{P}_{a,1}}{\tr_a \bar{\mathcal{K}}_a(1)}-\frac{1}{2}\sum_{k=1}^{N-1}\mathcal{R}'_{k-1,k}(1)\mathcal{P}_{k-1,k}-\frac{1}{4}\frac{{{\mathcal{K}}}_N'(1)}{{\mathcal{K}}_N(1)}\,,
\end{equation} 
see~\cite{Sklyanin:1988yz}. Here we used the symmetry of the R-matrix \eqref{eq:symm}. The first term in \eqref{eq:logderiv2} is proportional to the identity and can be subtracted without spoiling commutativity by setting 
\begin{equation}
 c_0=\frac{1}{4}\frac{\tr_a \bar{\mathcal{ K}}_a'(1)}{\tr_a \bar{\mathcal{K}}_a(1)}
\end{equation} 
in the Hamiltonian \eqref{eq:logderiv1}. The interesting points of the K-matrix that enter the Hamiltonian  are thus $\bar{ \mathcal{K}}(1)=M^{-1} \mathcal{K}(q^{-1})|_{R\to L}$ , $\mathcal{K}(1)$ and $\mathcal{K}'(1)$. 
As we are currently not aware of a generic expression of the K-matrix, we follow a procedure similar to the one proposed by Sutherland to determine the Hamiltonian density \cite{sutherland1970two}. 
More precisely we will study the boundary Yang-Baxter equations in the form given in \eqref{manga1} and \eqref{manga2} for the K-matrices at the points $y=q^{-1}$ and $y=1$.  A similar equation can be obtained for after taking the derivative of these equations also for $\mathcal{K}'(1)$ by differentiating the relations \eqref{manga1} and \eqref{manga2}. We will not attempt to compute the scalar $c_0$ proportional to the identity  as it is not relevant for our argument. 



\subsection{Bulk}\label{sec:bulkT}



At the special point $x=1$ the R-matrix above reduces to the permutation operator such that
\begin{equation} 
 \langle m,n|\mathcal{R}(1)|m',n'\rangle =\delta_{m,n'}\delta_{n,m'}
\end{equation} 
or $\mathcal{R}(1)=\mathcal{P}$. 
This follows from the property $
 \Phi_{q}(i|j;y;y)=\delta_{i,j}$ discussed in \cite{Mangazeev:2019rzf}. 
Using $
 \Phi_{q}(i|j;1;y)=\delta_{i,0}$
 which can also be found in \cite{Mangazeev:2019rzf}  we  obtain  
\begin{equation}
\begin{split}
 \langle m,n|\mathcal{R}'(1)\mathcal{P}|m',n'\rangle =\delta_{m+n,m'+n'} \bigg(&\partial_x\Phi_{q^2}(n'-n|n';x^{-2},x^{-2}q^{4s}) +\partial_x\Phi_{q^2}(m|m';x^2q^{4s},q^{4s})\bigg)\bigg|_{x=1}\,.
\end{split}
\end{equation} 
The two terms exactly lead  the rates of particles moving to the left and to the right. We find 
\begin{equation}\label{eq:hleft}
 \left.\partial_x\Phi_{q^2}(n'-n|n';x^{-2},x^{-2}q^{4s})\right|_{x=1}\delta_{m+n,m'+n'}=-2\alpha_+(n)\delta_{n,n'}\delta_{m,m'}+2 \beta_{+}(n',n'-n)\delta_{m+n,m'+n'}\delta_{n<n'}
\end{equation} 
\begin{equation}\label{eq:hright}
 \left.\partial_x\Phi_{q^2}(m|m';x^2q^{4s},q^{4s})\right|_{x=1}\delta_{m+n,m'+n'}=-2\alpha_-(m)\delta_{n,n'}\delta_{m,m'}+2 \beta_{-}(m',m'-m)\delta_{m+n,m'+n'}\delta_{m<m'}
\end{equation}  
and thus
\begin{equation}
\mathcal{H}=-\frac{1}{2}\mathcal{R}'(1)\mathcal{P}\,.
\end{equation} 
The   two relations \eqref{eq:hleft} and \eqref{eq:hright} are shown in Appendix~\ref{app:hopper}.
 


We remark that in the construction of the stochastic R-matrix \cite{Kuniba:2016fpi} the rates  arise from special points of the transfer operators leading to either left or right moving particles. Here we consider the full transfer matrix to obtain the natural generalisation to the q-Hahn process with hoppings in both directions from the logarithmic derivative. For the closed chain the q-Hahn Hamiltonian  has been constructed as a superposition   in \cite{Kuniba:2016fpi}.
 


 
 

\subsection{Left boundary}
The goal of this subsection is to obtain $\bar{\mathcal{K}}(1)$ such that the boundary term $B_L$  in \eqref{eq:bd1} arises from  \eqref{eq:logderiv2}. More precisely
\begin{equation}\label{eq:b1term}
B_L= 
 \frac{\tr_a \bar{\mathcal{K}}_a(1)\mathcal{H}_{a,1}}{\tr_a \bar{\mathcal{K}}_a(1)}\,,
\end{equation} 
where 
\begin{equation}
 \bar{\mathcal{K}}(1)=\left.D^{-1}\mathcal{K}(q^{-1})\right|_{R\to L}\,,
\end{equation}
with $D$ given in \eqref{eq:dmat}.
Here we take inspiration from the rational case where it was found that the K-matrix at the special point becomes a rank $1$ matrix.  

Inserting the equations \eqref{manga1} and \eqref{manga2} into a computer algebra system at $y=q^{-1}$ we find that for
\begin{equation}
 t^+_R=-\rho_R\,,\qquad  t^-_R=1\,,\qquad \kappa_R=1\,,\qquad \nu_R=-q^{2s-2}\,,
\end{equation} 
there exists such degenerate solution
\begin{equation}
 {\mathcal{K}}_{ij}\left(q^{-1}\right)=c_1\left(q^{2} \rho_R \right)^i \frac{  \left(q^{4 s};q^2\right){}_i}{\left(q^2;q^2\right){}_i}\,,
\end{equation} 
where $c_1$ denotes an overall constant that drops from the final result for the boundary term. It follows  that the other K-matrix can be expressed as
\begin{equation}
 \bar{\mathcal{K}}(1)=c_1\,  \rho_L  ^i \frac{  \left(q^{4 s};q^2\right){}_i}{\left(q^2;q^2\right){}_i}\,.
\end{equation} 
The trace immediately  yields
\begin{equation}
 \tr  \bar{\mathcal{K}}\left(1\right)=c_1\,  \frac{  \left(q^{4 s}\rho_L ;q^2\right){}_\infty}{\left(\rho_L;q^2\right){}_\infty}\,.
\end{equation} 
We can now compute the expression in \eqref{eq:b1term} using the explicit form of the Hamiltonian density obtained in the previous Section~\ref{sec:bulkT}. One finds
\begin{equation}\label{eq:trkh}
\begin{split}
 \frac{\tr_a \left( \bar{\mathcal{K}}_a\left(1\right)\mathcal{H}_{a1}\right)}{ \tr  \bar {\mathcal{K}}\left(1\right)}|m_1\rangle=\left(\alpha_+(m_1) +\sum_{k=1}^\infty\frac{\rho_L^k}{1-\gamma^{k}}\right)|m_1\rangle&
-\sum_{k=1}^{m_1} \beta_+(m_1,k)|m_1-k\rangle\\& -\sum_{k=1}^\infty\frac{\rho_L^k}{1-\gamma^{k}}|m_1+k\rangle \,.
\end{split}
\end{equation} 
The details of the computation are given in Appendix~\ref{app:leftbnd}.


\subsection{Right boundary}\label{sec:rib}

To show that the boundary term $B_R$ arises from the transfer matrix via
\begin{equation}\label{eq:BRR}
 B_R=-\frac{1}{4}\frac{{{\mathcal{K}}}_N'(1)}{{\mathcal{K}}_N(1)}\,,
\end{equation} 
 we differentiate the two relations \eqref{manga1} and \eqref{manga2}.
We impose that the K-matrix is proportional to the identity at the special point $y=1$
\begin{equation}
 \mathcal{K}(1)=\frac{1}{4}I\,,
\end{equation} 
which is consistent with  \eqref{manga1} and \eqref{manga2}
and set
\begin{equation}\label{eq:identification}
  t^+_R=-\rho_R\,,\qquad  t^-_R=1\,,\qquad \kappa_R=1\,,\qquad \nu_R=-q^{-2s}\,,\,.
\end{equation}  
Then the equations arising from the boundary Yang-Baxter equation above turn into
\begin{equation}\label{manga1h}
 \begin{split}
  -\rho_R q^{2-4s}\left(1-q^{2(j+2s)}\right)\mathcal{K}'_{j,l+1}(1)
 +q^{-4s}\left(1-q^{2+2l}\right)\mathcal{K}'_{j+1,l}(1)\\
 -q^{2}\left(q^{2j}-q^{2l}\right)\left(1+q^{2-4s}\rho_R\right)\left(\mathcal{K}'_{j+1,l+1}(1)+\frac{1}{2}\delta_{j+1,l+1}\right)\\
 +\rho_R q^{2-4s}\left(1-q^{2(1+l+2s)}\right)\left(\mathcal{K}'_{j+1,l+2}(1)+\delta_{j+1,l+2}\right)\\
 - q^{-4s}\left(1-q^{4+2j}\right)\left(\mathcal{K}'_{j+2,l+1}(1)+\delta_{j+2,l+1}\right)=0\,,
 \end{split}
\end{equation} 
and
\begin{equation}\label{manga2h}
 \begin{split}
 -\rho_R q^{2(2+l)}\left(1-q^{2(j+2s)}\right)\mathcal{K}'_{j,l+1}(1)
 +q^{2(1+j)}\left(1-q^{2+2l}\right)\mathcal{K}'_{j+1,l}(1)\\
 -q^{2}\left(q^{2j}-q^{2l}\right)\left(\rho_R+1\right)\mathcal{K}'_{j+1,l+1}(1)
 +\rho_R q^{2(1+j)}\left(1-q^{2(1+l+2s)}\right)\mathcal{K}'_{j+1,l+2}(1)\\
 - q^{2l}\left(1-q^{4+2j}\right)\mathcal{K}'_{j+2,l+1}(1)=0\,.
 \end{split}
\end{equation} 
 Using a computer algebra system we can verify that 
\begin{equation}\label{eq:trk}
 \mathcal{K}'_{i,j}(1)=-\left(\alpha _-(i)+\sum_{k=1}^\infty\frac{\rho_R^k}{q^{-2k}-1}\right)\delta_{i,j}+\beta_-(j,j-i)\delta_{i< j}+\frac{\rho_R^{i-j}}{q^{-2(i-j)}-1}\delta_{i>j}\,,
\end{equation}  
solves both of the equations \eqref{manga1h} and \eqref{manga2h}. To keep track of the differerent contributions on the diagonal, upper and lower triangular part  in \eqref{eq:trk} we distinguish five cases that are exemplified in Figure~\ref{fig:mat} for which the equations are verified separately.
We thus recover the right boundary in \eqref{eq:bd2} from \eqref{eq:BRR}.
 


\begin{figure*}[t!]
    \centering
    \begin{minipage}{.19\columnwidth}
        \pfive
        \caption*{$j-l>1\qquad$}
    \end{minipage}
    \begin{minipage}{.19\columnwidth}
      \pfour
        \caption*{$j-l=1\qquad$}
    \end{minipage}
    \begin{minipage}{.19\columnwidth}
    \pthree
        \caption*{$j=l\qquad$}
    \end{minipage}
    \begin{minipage}{.19\columnwidth}
        \ptwo
        \caption*{$l-j=1\qquad$}
    \end{minipage}
    \begin{minipage}{.19\columnwidth}
        \pone
        \caption*{$l-j>1\qquad$}
    \end{minipage}
    \caption{Case distinctions of the relations \eqref{manga1h} and \eqref{manga2h}. The boxes denote the matrix entries of the K-matrix with the diagonal  colored in  grey. The entries affected by the relations \eqref{manga1h} and \eqref{manga2h} are indicated in red color.}
    \label{fig:mat}
\end{figure*}

 

\section{Limiting case of open q-Hahn TASEP}\label{sec:otasep}
The open q-Hahn process \eqref{eq:ham} naturally yields a totally asymmetric process in the limit when $\mu$ tends to zero. The bulk term turns into
\begin{equation}\label{eq:hactb}
\begin{split}
 \mathcal{H}|m\rangle \otimes |m'\rangle=\alpha_-(m) |m\rangle \otimes |m'\rangle-\sum_{k=1}^m \beta_-(m,k) |m-k\rangle \otimes |m'+k\rangle \,,
 \end{split}
\end{equation} 
while the boundaries become 
\begin{equation}\label{eq:bd1b}
\begin{split}
B_L|m\rangle = \sum_{k=1}^\infty\frac{\rho_L^k}{1-\gamma^{k}}|m\rangle
  -\sum_{k=1}^\infty\frac{\rho_L^k}{1-\gamma^{k}}|m+k\rangle \,,
\end{split}
\end{equation} 
and
\begin{equation}\label{eq:bd2b}
\begin{split}
B_R|m\rangle = \alpha_-(m)  |m\rangle
-\sum_{k=1}^m\beta_-(m,k)|m-k\rangle  \,.
\end{split}
\end{equation}  
This Hamiltonian describes a process where particles are inserted at the first site with a rate governed by the parameter $\rho_L$ and are then transported to site $N$ at rate $\beta_-$. At the last site they are extracted at the same rate $\beta_-$.

\section{Conclusion}\label{sec:conc}

The definition of the q-Hahn process with open boundary conditions allows for the study of its microscopic properties along the lines of its prominent finite-dimensional cousin known as the Asymmetric Simple Exclusion Process (ASEP). In the $q\to 1$ limit we recover the symmetric hopping model with boundaries introduced in \cite{Frassek:2019vjt} while $\mu\to0$ leads to an analogy of the open TASEP. It would be interesting to generalise the results obtained so far in the symmetric setting to the asymmetric one. In particular, for the symmetric model the steady state has been computed exactly in \cite{Frassek:2021yxb} and a dual Levy process was obtained in \cite{Franceschini}. 

Further, it would be desirable to derive the K-matrix that yields the introduced boundary conditions and understand the relation to the K-matrices introduced in \cite{Lazarescu,Mangazeev:2019rzf}. As in the rational case that was studied in \cite{Frassek:2019vjt} the K-matrices should arise from a universal K-operator first introduced in the same article. In case of the q-deformed model,  such K-operator has only been found for the triangular case \cite{Tsuboi:2019qac}. 
One may also ask whether the boundary conditions we introduced are the most general stochastic boundaries. This question can likely be answered once the K-matrix is known.

Finally  it would be very interesting to study whether our prescription to obtain boundary terms of the integrable Hamiltonian  can be formalised leading to a method that allows bootstrapping integrable Hamiltonians for open chains in analogy to Sutherland's approach.  In particular we are not aware of a reference about the special point where the K-matrix reduces to rank $1$ in the literature.
 


\paragraph{Acknowledgments}
I like to thank Cristian Giardin\`{a}, Istv\'{a}n M. Sz\'{e}cs\'{e}nyi and  especially 
Rodrigo A. Pimenta for inspiration and fruitful discussions. Further I thank the organisers and participants of the GGI workshop \textit{``Randomness, Integrability and Universality''}, where part of the work has been carried out, for creating a stimulating environment.



