Harmonic Process in Nonequilibrium Systems

• The harmonic process is a boundary-driven continuous-time Markov model with unbounded site occupancy and explicit algebraic steady-state representations.
• It offers three exact formulations—closed gamma-product, nested-integral, and matrix-product ansatz—that enable precise analysis of asymptotic behavior and fluctuations.
• Its explicit stationary measure and mixture-of-products structure underpin rigorous derivations of local equilibrium, central limit theorems, and large deviation principles in non-reversible systems.

The harmonic process is an exactly solvable boundary-driven continuous-time Markov model on a one-dimensional chain with unbounded site occupancy, developed in nonequilibrium statistical mechanics and integrable probability as a tractable non-equilibrium steady-state (NESS) system (Frassek, 27 Jul 2025). In its algebraic formulation, the dynamics are written in “Hamiltonian” form and depend on a spin label $2s\in\mathbb N$; in the boundary-driven setting, the stationary law admits several exact representations, including a closed gamma-product formula, a nested-integral or mixed-measure form, a matrix-product ansatz, and a mixture of product geometric measures (Frassek, 27 Jul 2025, Redig et al., 19 Apr 2026). These representations support a detailed theory of local equilibrium, laws of large numbers, central limit theorems, and large deviations for local observables in a non-reversible interacting system (Redig et al., 19 Apr 2026).

1. Definition of the model and its dynamics

In the formulation emphasized in the steady-state analysis, the state space consists of configurations

$|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$

on a chain of length $N$, with unbounded occupancy at each site. The process is a continuous-time Markov chain whose generator is written as

$$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$

acting on

$$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$

and the steady state satisfies $H|\mu\rangle=0$ (Frassek, 27 Jul 2025).

The bulk interaction between neighboring sites is governed by the harmonic weights

$h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$

and the jump kernel

$$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$

The local bulk operator acts by redistributing $k$ particles from one site to its neighbor, with rates determined by $\phi_s$. Boundary reservoirs inject according to parameters $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$0 and extract with matching harmonic rates. The parameter $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$1 controls both the bulk-jump spectrum and the harmonic weights (Frassek, 27 Jul 2025).

A boundary-driven formulation used for ergodic analysis is written on

$|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$2

with configuration $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$3. There the bulk dynamics replace a pair $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$4 by a uniformly chosen split of the conserved total $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$5, while the left and right reservoirs replace the boundary occupation by geometric laws $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$6 and $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$7, respectively. The bulk conserves total mass, whereas the reservoirs exchange mass with the environment and create a genuine NESS when $|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$8 (Redig et al., 19 Apr 2026).

2. Exact stationary measure

The stationary state is explicitly computable. In the quantum-inverse-scattering treatment, the unique normalized zero-mode can be represented as two nonlocal rotations acting on a telescopic product vector, with reservoir variables

$|m_1,\dots,m_N\rangle,\qquad m_i\in\mathbb N,$9

The resulting stationary weight can be written in fully explicit form as

$N$0

where

$N$1

and

$N$2

This is the closed-form formula derived by Frassek et al. (Frassek et al., 2021).

For the boundary-driven harmonic model treated probabilistically, the NESS can also be expressed as a mixture of product geometric laws. Writing

$N$3

with mean $N$4, one has

$N$5

In the harmonic NESS, the random vector $N$6 is distributed as the order statistics of $N$7 independent Uniform$N$8 variables, equivalently

$N$9

where $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$0 are the usual order statistics. Conditional on $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$1, the sites are independent with laws $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$2 (Redig et al., 19 Apr 2026).

3. Three exact steady-state representations

A central structural result is that the harmonic-process NESS admits three equivalent exact representations: a closed gamma-product formula, a nested-integral representation, and a matrix-product ansatz (Frassek, 27 Jul 2025, Carinci et al., 2023, Carinci et al., 2023).

Carinci–Giardinà et al. rewrote the stationary weight as an $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$3-fold iterated integral over auxiliary variables $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$4: $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$5

$$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$6

with $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$7 and $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$8. This representation gives the stationary state a mixed-measure interpretation in terms of auxiliary reservoir variables (Carinci et al., 2023, Carinci et al., 2023).

The new contribution of the matrix-product treatment is an explicit oscillator realization. One seeks operators $$H=\mathcal B_L+\sum_{i=1}^{N-1}\mathcal H_{i,i+1}+\mathcal B_R,$$9 and boundary vectors $$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$0, $$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$1 such that

$$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$2

Using oscillators $$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$3 with $$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$4 and Fock basis $$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$5, an explicit choice is

$$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$6

with boundary vectors

$$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$7

For this choice, the matrix product reproduces the closed gamma-product weight exactly (Frassek, 27 Jul 2025).

The relative roles of the three representations are distinct.

Representation	Core form	Typical use
Closed form	Telescopic gamma-product	Asymptotic analysis, large-charge limits
Nested integrals	$$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$8-fold mixed-measure integral	Functional identities, correlation integral formulas
Matrix-product ansatz	Infinite-dimensional oscillator MPA	Correlation functions, algebraic steady-state computations

The closed form is the most explicit. The nested-integral form is flexible for functional relations in the reservoir parameters. The MPA is the most algebraic, naturally connected with the integrability framework and direct expectation-value computations (Frassek, 27 Jul 2025).

4. Transformations between representations

The equivalence of the three steady-state formulas is itself an exact structural result. From the closed gamma-product formula, one obtains the nested-integral representation by repeatedly applying the Beta-function identity

$$|P\rangle=\sum_{\vec m}P(\vec m)|m_1,\dots,m_N\rangle,$$9

Each ratio of gamma functions is converted into a one-fold integral; after $H|\mu\rangle=0$0 iterations, the telescopic product becomes an $H|\mu\rangle=0$1-fold nested integral (Frassek, 27 Jul 2025).

The passage from nested integrals to the MPA is mediated by an integral operator

$H|\mu\rangle=0$2

Interpreting $H|\mu\rangle=0$3 as an infinite matrix yields entries coinciding with $H|\mu\rangle=0$4, and the full nested integral becomes

$H|\mu\rangle=0$5

Conversely, the MPA recovers the closed form either by diagonalizing the oscillators or by computing $H|\mu\rangle=0$6 through successive action on the Fock basis, which reproduces the gamma-product expression (Frassek, 27 Jul 2025).

At the algebraic level, stationarity is encoded by a bulk relation and matching boundary relations involving $H|\mu\rangle=0$7 and an auxiliary family $H|\mu\rangle=0$8. In the summary form,

$H|\mu\rangle=0$9

together with left and right boundary identities for $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$0 and $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$1. These relations guarantee $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$2 and place the harmonic process within the same matrix-product and Zamolodchikov-algebra paradigm that appears in other integrable boundary-driven systems (Frassek, 27 Jul 2025).

5. Local equilibrium, fluctuations, and large deviations

Because the stationary measure is an explicit mixture of product states, the harmonic process admits a detailed ergodic theory. The local-equilibrium profile is linear: $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$3 For every fixed multi-index $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$4 and $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$5,

$h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$6

almost surely. A quantitative concentration bound for the order-statistic parameters is

$h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$7

For self-duality polynomials $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$8, the deviation from local equilibrium is controlled at order $h_s(m)=\sum_{k=1}^m \frac{1}{2s+k-1},$9 uniformly in $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$0 (Redig et al., 19 Apr 2026).

For a local observable $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$1 of range $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$2 and a test function $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$3, the empirical field

$$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$4

satisfies a strong law of large numbers: $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$5 The proof decomposes $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$6 into its conditional expectation given $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$7 and a centered remainder; the first part is handled as an $$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$8-statistic in order statistics, the second as a triangular array with uniform moment bounds (Redig et al., 19 Apr 2026).

The centered fluctuation field

$$\phi_s(m,k)=\frac1k\cdot \frac{\Gamma(m+1)\Gamma(m-k+2s)}{\Gamma(m-k+1)\Gamma(m+2s)}.$$9

converges in law to a Gaussian with variance

$k$0

Here

$k$1

is the equilibrium white-noise contribution, while

$k$2

comes from fluctuations of the random environment $k$3, which converge to a Brownian bridge. For the density field $k$4, the limit splits into a Brownian-bridge term plus centered white noise (Redig et al., 19 Apr 2026).

Large deviations are also explicit. For

$k$5

one obtains an LDP at speed $k$6. The random parameter profile $k$7 has rate function

$k$8

when $k$9 is absolutely continuous, strictly increasing, and satisfies $\phi_s$0, $\phi_s$1, and $\phi_s$2 otherwise. The final rate function for the mixed measure is the variational formula

$\phi_s$3

This yields LLN, CLT, and LDP within a single exact framework (Redig et al., 19 Apr 2026).

6. Analytical significance and methodological role

The harmonic process occupies a distinctive position because its NESS is simultaneously accessible by integrable, probabilistic, and operator-algebraic methods. On the integrability side, the closed-form stationary state was obtained by the quantum-inverse-scattering method, the nested-integral representation exposes a mixed-measure structure, and the matrix-product ansatz furnishes a new oscillator-based auxiliary-space realization that had not previously been available for the model (Frassek et al., 2021, Carinci et al., 2023, Carinci et al., 2023, Frassek, 27 Jul 2025).

On the probabilistic side, the mixture-of-products structure reduces fluctuation questions to two components: the product-measure part and the random profile formed by order statistics. The proofs then combine empirical-process theory for order statistics, strong laws and central limit theorems for triangular arrays, a positive maximum-principle perspective, self-duality, and Varadhan’s lemma (Redig et al., 19 Apr 2026). In this sense, the harmonic process provides a complete microscopic derivation of macroscopic fluctuation results for a non-reversible, boundary-driven interacting system (Redig et al., 19 Apr 2026).

The three exact stationary representations are complementary rather than redundant. The closed form is particularly suited to asymptotic analysis as $\phi_s$4 or in large-charge regimes. The nested-integral form is adapted to functional identities and correlation-function integral formulas. The matrix-product ansatz is adapted to algebraic computations of observables, including 1- and 2-point functions, without combinatorial explosion, even though its auxiliary space is infinite-dimensional and realized on oscillator Fock space (Frassek, 27 Jul 2025).

Taken together, these results define the harmonic process as a rare exactly solvable nonequilibrium lattice model for which stationary structure, algebraic representation, and asymptotic fluctuation theory can all be developed explicitly and in mutually compatible forms.