Exact non-equilibrium steady state of the boundary-driven multispecies harmonic process: factorised form and mixture representation

The boundary-driven multispecies harmonic process is a continuous-time Markov chain on a chain of $N$ sites, each hosting an unbounded number of particles of $M$ species, whose Markov generator coincides with the stochastic Hamiltonian of an open rational $gl(M+1)$ Heisenberg spin chain in non-compact (infinite-dimensional) symmetric-tensor representations. Casini, Frassek and Giardinà (2026) constructed the underlying factorised $R$-matrix, new non-diagonal stochastic $K$-matrices solving the boundary Yang-Baxter equation, and the Sklyanin double-row transfer matrix that generates the process, together with absorbing, hidden-parameter and heat-conduction dual processes; however, they left the steady state of the process undetermined, explicitly deferring its study to a follow-up article. The open problem is to determine the explicit non-equilibrium steady state, including its factorised form and the proof of the mixture representation in terms of products of negative multinomials. The monospecies case is completely solved in the literature; the steady state of the multispecies process remains open.

Background

Integrable stochastic particle systems are Markov chains whose generators coincide with the Hamiltonians of integrable quantum spin chains, which makes the tools of the quantum inverse scattering method — transfer matrices, Bethe ansatz, $Q$-operators — available for exact non-equilibrium computations. The harmonic process of Frassek, Giardinà and Kurchan (2020) is such a chain on a one-dimensional lattice of $N$ sites in which each site can hold an unbounded number of particles: a site occupied by $m$ particles sheds $k$ of them, chosen uniformly, to a given neighbour at rate $\varphi_s(k,m)=\Gamma(k)\,\Gamma(2s+m-k)\,\Gamma(m+1)/\big(\Gamma(2s+m)\,\Gamma(k+1)\,\Gamma(m-k+1)\big)$ with spin parameter $2s>0$; in the canonical harmonic case $2s=1$ this rate is $1/k$, and the diagonal rate is $\psi(2s+m)-\psi(2s)$ with $\psi$ the digamma function, which becomes the harmonic number $H_m$. Casini, Frassek and Giardinà (2026) introduced the multispecies ($M$-colour) version, in which a site holding $m_a$ particles of species $a$ moves a vector $(k_1,\dots,k_M)$ of particles with rate $\varphi_s(k,m)=\frac{\Gamma(|k|)\,\Gamma(2s+|m|-|k|)}{\Gamma(2s+|m|)}\prod_{a=1}^{M}\binom{m_a}{k_a}$ (the integrable deformation of an improper Dirichlet-multinomial sampling), and put the chain in contact with boundary reservoirs, encoded by parameters $\beta_{l,a}$ and $\beta_{r,a}$, that inject and remove particles. They proved the integrability of the resulting process: the Markov generator $\mathcal H$ is the stochastic Hamiltonian of an open rational $gl(M+1)$ Heisenberg spin chain acting on non-compact, that is infinite-dimensional, symmetric-tensor representations of $gl(M+1)$ (the infinite dimensionality encodes the unbounded local occupation numbers). The construction rests on three ingredients. First, an $R$-matrix, an operator solving the Yang-Baxter equation on a tensor product of two infinite-dimensional spaces, which factorises into two $R$-operators associated with right- and left-moving particles and is realised through $M$ conjugate Heisenberg (oscillator) pairs. Second, $K$-matrices, solutions of the boundary Yang-Baxter (reflection) equations that encode the reservoirs; these are non-diagonal, allowing injection and extraction of particles at the chain ends, and become triangular when the right-reservoir strength $\rho_r$ vanishes, in which case they match the rational limit of the triangular $U_q(\widehat{sl_n})$ reflection-equation solutions of Kolyaskin and Mangazeev. Third, Sklyanin's double-row transfer matrix $T(x)=\mathrm{tr}_a\big(\widetilde{K}(x)\,\mathcal R_{a,1}(x)\cdots\mathcal R_{a,N}(x)\,\mathcal K_a(x)\,\mathcal R_{N,a}(x)\cdots\mathcal R_{1,a}(x)\big)$, whose logarithmic derivative at $x=0$ yields $\mathcal H$, so the process possesses a commuting family of conserved quantities; a companion transfer matrix $T_f(x)$ with the fundamental (finite-dimensional) $gl(M+1)$ representation in the auxiliary space commutes with $T(x)$ and is the object that plays the central role in the algebraic Bethe ansatz. What is already known for a single species ($M=1$) sets the target for the multispecies answer. The non-equilibrium steady state (NESS), that is the unique stationary distribution of the boundary-driven chain, was found in closed factorised form through a similarity transformation to an isospectral Hamiltonian $\mathcal H''$ with triangular boundary terms whose ground state factorises (Frassek and Giardinà 2022); equivalent nested-integral and unifying matrix-product representations are known (Carinci et al. 2023; Frassek 2025). The NESS is a mixture, i.e. a convex combination, of products of geometric distributions whose local parameters are random and distributed as uniform order statistics, equivalently as an ordered Dirichlet distribution; the parameters themselves evolve as a Markov process on the parameter space, the 'hidden parameter model', and this mixture structure propagates in time (Giardinà, Redig and van Tol 2025). The ordered Dirichlet law is the unique mixing measure with a spatial Markov property under natural scaling assumptions (Redig and van Tol 2024), and the resulting NESS supports full law-of-large-numbers, central-limit and large-deviation results (Redig and van Tol 2026; Carinci et al. 2023). On the spectral side, Bethe-ansatz eigenstates of triangularisable (conjugate, by a similarity transformation, to a triangular matrix) open XXX (isotropic Heisenberg) spin chains underlie the closed-form steady states of the open simple exclusion process (SSEP) and of the boundary-driven multiparticle asymmetric diffusion model (Frassek 2020; Frassek and Szécsényi 2023), and Baxter $Q$-operators — auxiliary transfer matrices whose eigenvalues encode the Bethe roots — are known for the open ASEP (Lazarescu and Pasquier 2014), for open XXX chains of arbitrary spin (Frassek and Szécsényi 2022), for non-compact closed chains of arbitrary rank (Derkachov and Manashov 2006), and very recently for open non-compact $SL(2,\mathbb C)$ chains with boundary $K$-matrix (Antonenko, Derkachov and Valinevich 2025); none of these covers open non-compact $gl(M+1)$ chains with the stochastic $K$-matrices relevant here. For $M\ge 2$ species nothing beyond the integrability structure and duality relations is known. The source paper derives three dual processes — an absorbing dual, a hidden-parameter model, and a heat-conduction model — where the absorbing duality expresses moments of order $n$ of the occupation numbers through $n$ dual particles that are eventually absorbed, and it observes that when all left and right reservoir parameters coincide the process is reversible with a product of negative multinomial distributions (the multivariate generalisation of the geometric distribution) as reversible measure. The same paper states that the broader classification of stochastic integrable boundary conditions for non-compact vertex models in terms of $K$-matrices is an open problem with only partial results (diagonal solutions, triangular solutions, $q$-Hahn boundaries); that classification program is adjacent context for this problem rather than part of its completion boundary, since the $K$-matrices needed for the harmonic process are already constructed there.

Problem Statement

For the boundary-driven multispecies harmonic process of arXiv:2607.26262 — the continuous-time Markov chain on $N$ sites with $M\ge 2$ particle species and unbounded local occupation numbers, with bulk jump rates $\varphi_s(k,m)=\frac{\Gamma(|k|)\,\Gamma(2s+|m|-|k|)}{\Gamma(2s+|m|)}\prod_{a=1}^{M}\binom{m_a}{k_a}$ for jumps $0<|k|\le|m|$ at admissible spin $2s>0$, and with reservoirs injecting and removing particles at the two ends according to the boundary parameters $\beta_{l,a},\beta_{r,a}\in(0,1)$ satisfying $\sum_{a=1}^{M}\beta_{l,a}<1$ and $\sum_{a=1}^{M}\beta_{r,a}<1$ — whose Markov generator $\mathcal H$ is obtained, up to the additive identity shift and overall normalisation fixed in that paper, from the logarithmic derivative at $x=0$ of the Sklyanin double-row transfer matrix built from the factorised non-compact $gl(M+1)$ $R$-matrix and the stochastic $K$-matrices constructed in that paper, determine the explicit non-equilibrium steady state of the process.

Non-equilibrium steady state. Determine the stationary vector $|\Psi\rangle$, $\mathcal H|\Psi\rangle=0$, explicitly. This comprises: (i) its factorised form — the source reduces this problem, through the similarity transformations given there, to the ground state of the isospectral Hamiltonian $\mathcal H''$ with triangular and diagonal boundary terms, conjugated back to the stochastic Hamiltonian; and (ii) the proof of the mixture representation in analogy with the monospecies case: show that the stationary measure is a mixture of products of negative multinomial distributions whose parameter vectors are distributed according to the stationary law of the multispecies hidden-parameter model, and identify that stationary mixing law (the $M$-species analogue of the ordered Dirichlet distribution of uniform order statistics that appears at $M=1$).

An answer must be consistent with the exactly known limits: at $M=1$ it must reduce to the known monospecies steady state in its closed, matrix-product and mixture forms with ordered-Dirichlet mixing, and when all left and right reservoir parameters coincide, $\beta_{l,a}=\beta_{r,a}$ for all $a$, it must reproduce the reversible product of negative multinomials. The broader classification of stochastic integrable $K$-matrices for non-compact vertex models, also stated as open in the source, is not part of this problem's completion boundary: the $K$-matrices required above are already available. Likewise, a spectral solution via Baxter $Q$-operators or the algebraic Bethe ansatz for the commuting transfer matrices $T(x)$ and $T_f(x)$, although a natural direction suggested by the integrability structure constructed in the source, is not posed there as an open problem and is not part of this problem's completion boundary.

References

1. F. Casini, R. Frassek, C. Giardinà, 'The boundary-driven multispecies harmonic process', arXiv:2607.26262 (2026). https://arxiv.org/abs/2607.26262
2. R. Frassek, C. Giardinà, J. Kurchan, 'Non-compact quantum spin chains as integrable stochastic particle processes', Journal of Statistical Physics 180, 135–171 (2020), arXiv:1904.01048, DOI 10.1007/s10955-019-02375-4. https://arxiv.org/abs/1904.01048
3. R. Frassek, C. Giardinà, 'Exact solution of an integrable non-equilibrium particle system', Journal of Mathematical Physics 63, 103301 (2022), arXiv:2107.01720, DOI 10.1063/5.0086715. https://arxiv.org/abs/2107.01720
4. R. Frassek, 'Steady state representations for the harmonic process', Journal of Physics A: Mathematical and Theoretical 59, 085202 (2026), arXiv:2507.20436, DOI 10.1088/1751-8121/ae4799. https://arxiv.org/abs/2507.20436
5. G. Carinci, C. Franceschini, R. Frassek, C. Giardinà, F. Redig, 'Large deviations and additivity principle for the open harmonic process', Communications in Mathematical Physics 406, 103 (2025), arXiv:2307.14975, DOI 10.1007/s00220-025-05271-z. https://arxiv.org/abs/2307.14975
6. G. Carinci, C. Franceschini, D. Gabrielli, C. Giardinà, D. Tsagkarogiannis, 'Solvable stationary non equilibrium states', Journal of Statistical Physics 191, 10 (2024), arXiv:2307.02793, DOI 10.1007/s10955-023-03226-z. https://arxiv.org/abs/2307.02793
7. C. Giardinà, F. Redig, B. van Tol, 'Intertwining and propagation of mixtures for generalized KMP models and harmonic models', Journal of Statistical Physics 192, 21 (2025), arXiv:2406.01160, DOI 10.1007/s10955-025-03393-1. https://arxiv.org/abs/2406.01160
8. F. Redig, B. van Tol, 'Non-equilibrium steady states with a spatial Markov structure', Journal of Statistical Physics 192, 89 (2025), arXiv:2411.11425, DOI 10.1007/s10955-025-03471-4. https://arxiv.org/abs/2411.11425
9. F. Redig, B. van Tol, 'Ergodic properties of the harmonic process', arXiv:2604.17469 (2026). https://arxiv.org/abs/2604.17469
10. R. Frassek, 'Eigenstates of triangularisable open XXX spin chains and closed-form solutions for the steady state of the open SSEP', Journal of Statistical Mechanics (2020) 053104, arXiv:1910.13163, DOI 10.1088/1742-5468/ab7af3. https://arxiv.org/abs/1910.13163
11. R. Frassek, I. M. Szécsényi, 'The steady state of the boundary-driven multiparticle asymmetric diffusion model', Journal of Physics A: Mathematical and Theoretical 57, 095205 (2024), arXiv:2311.03603, DOI 10.1088/1751-8121/ad2725. https://arxiv.org/abs/2311.03603
12. F. Casini, R. Frassek, C. Giardinà, 'Duality for the multispecies stirring process with open boundaries', Journal of Physics A: Mathematical and Theoretical 57, 295001 (2024), arXiv:2312.15532, DOI 10.1088/1751-8121/ad5ac9. https://arxiv.org/abs/2312.15532
13. A. Lazarescu, V. Pasquier, 'Bethe Ansatz and Q-operator for the open ASEP', Journal of Physics A 47, 295202 (2014), arXiv:1403.6963, DOI 10.1088/1751-8113-47-29/295202. https://arxiv.org/abs/1403.6963
14. R. Frassek, I. M. Szécsényi, 'Algebraic Bethe ansatz for Q-operators of the open XXX Heisenberg chain with arbitrary spin', Journal of Physics A: Mathematical and Theoretical 55, 505201 (2022), arXiv:2208.02197, DOI 10.1088/1751-8121/aca5d3. https://arxiv.org/abs/2208.02197
15. S. E. Derkachov, A. N. Manashov, 'R-matrix and Baxter Q-operators for the noncompact SL(N,C) invariant spin chain', SIGMA 2, 084 (2006), arXiv:nlin/0612003, DOI 10.3842/SIGMA.2006.084. https://arxiv.org/abs/nlin/0612003
16. P. Antonenko, S. Derkachov, P. Valinevich, 'BC-type open SL(2,C) spin chain', Annals of Henri Poincaré (2026), arXiv:2508.04972, DOI 10.1007/s00023-025-01653-0. https://arxiv.org/abs/2508.04972
17. P. Antonenko, N. Belousov, S. Derkachov, S. Khoroshkin, 'Reflection operator and hypergeometry I: SL(2,R) spin chain', Zapiski Nauchnykh Seminarov POMI 532, 5–46 (2024); English translation in Journal of Mathematical Sciences 299, 691–719 (2026), arXiv:2406.19862, DOI 10.1007/s10958-026-08467-3. https://arxiv.org/abs/2406.19862
18. V. V. Mangazeev, X. Lu, 'Boundary matrices for the higher spin six vertex model', Nuclear Physics B 945, 114665 (2019), arXiv:1903.00274, DOI 10.1016/j.nuclphysb.2019.114665. https://arxiv.org/abs/1903.00274
19. D. Kolyaskin, V. V. Mangazeev, 'Triangular solutions to the reflection equation for U_q(sl_n)', Journal of Physics A 57, 245201 (2024), arXiv:2402.05442, DOI 10.1088/1751-8121/ad4d2f. https://arxiv.org/abs/2402.05442
20. R. Frassek, 'Integrable boundaries for the q-Hahn process', Journal of Physics A 55, 404008 (2022), arXiv:2205.10512, DOI 10.1088/1751-8121/ac901b. https://arxiv.org/abs/2205.10512
21. Z. Tsuboi, 'On diagonal solutions of the reflection equation', Journal of Physics A 52, 155201 (2019), arXiv:1811.10407, DOI 10.1088/1751-8121/ab0b6d. https://arxiv.org/abs/1811.10407
22. E. K. Sklyanin, 'Boundary conditions for integrable quantum systems', Journal of Physics A 21, 2375–2389 (1988), DOI 10.1088/0305-4470/21/10/015. https://doi.org/10.1088/0305-4470/21/10/015
23. C. Giardinà, M. Möhle, 'On a Λ-mutation model and the harmonic model', arXiv:2606.00781 (2026). https://arxiv.org/abs/2606.00781