Introduction — Long-range correlations observed in high-multiplicity p+p [1], [2] and p+Pb collisions [3]–[5] resemble correlation structures associated with collective dynamics in heavy-ion collisions [6], but appear in systems that are much smaller. This raises the question whether they reflect collective dynamics or arise from other long-range sources [7], [8]. Interpreting these measurements requires separating collective dynamics from nonflow and finite-multiplicity effects. Experimentally, azimuthal and transverse-momentum correlations are often measured with pseudorapidity gaps or subevent methods to suppress short-range nonflow contributions [9]–[12]. These procedures reduce local correlations from jets, resonance decays, and other short-range sources, but they do not remove correlations imposed by exact global conservation laws [13]–[17]. Conservation-induced correlations are intrinsically long range because particles in separated rapidity intervals must share the same globally conserved energy, momentum, and charges. Since conservation effects are less suppressed in small systems and peripheral events, they constitute an unavoidable contribution to measured long-range correlations [18], [19]. Existing treatments often compute exact-conservation effects as corrections to observables [14], [18], [19]. Here we derive leading conservation-induced correlations for a general set of additively conserved quantities.

The derivation follows from formulating equilibrium statistics as conditional probability under additive conservation laws, rather than through ensembles, reservoirs, or entropy maximization [20], [21]. This viewpoint is closely related to conditional limit theorems, the Gibbs conditioning principle, and large-deviation formulations of statistical mechanics [22]–[29]. Here we formulate a many-mode system as an exact conditional probability law under additive conservation constraints. This provides a unified framework for equilibrium occupation statistics and conservation-induced correlations. The exact one-mode conditional law gives relativistic Maxwell–Boltzmann (MB), Bose–Einstein (BE), and Fermi–Dirac (FD) statistics at leading saddle order. The conserved quantities fix the common exponential tilt, while the allowed occupations and microscopic occupation weights determine the statistics. To best of our knowledge, this is the first derivation in which the relativistic Maxwell–Jüttner exponential form [30], [31] arises from conditional probability, rather than by postulating a canonical ensemble or maximizing entropy.

The same conditional law also gives the correlations imposed by exact conservation. Expanding the exact two-mode conditional law to Gaussian order gives the leading conservation-induced covariance between modes. For observables linear in the mode occupations, this covariance gives the leading exact-conservation contribution to long-range correlations. It depends on each observable only through its overlap with the selected conserved quantities. This makes it possible to construct projected observables whose leading exact-conservation contribution vanishes by construction. We use this construction to define experimentally accessible long-range observables and demonstrate it with PYTHIA8/Angantyr-generated [32], [33] p+Pb events at \(\sqrt{s_{\mathrm{NN}}}=5.02~\mathrm{TeV}\). The projection removes the conservation-aligned part of the covariance while leaving the conservation-orthogonal covariance unchanged.

Conditional framework — We consider a closed system with particles distributed among single-particle modes \(i\). Let \(n_i\) be the occupation number of mode \(i\). We denote the exactly conserved quantities by an \(r\)-component vector \(Q\), with components \[Q^A=\sum_i q_i^A n_i \,, \qquad A=1,\ldots,r \,,\] where \(q_i^A\) is the \(A\)-th component of the conserved-quantity vector carried by mode \(i\). A configuration \(\{n_i\}\) is counted with microscopic weight \(\mathbb{W}(\{n_i\})\), and the constrained density of configurations is \[\Omega(Q) = \sum_{\{n_i\}} \mathbb{W}(\{n_i\}) \delta^{(r)}\!\left(\sum_i q_i n_i-Q\right). \label{eq:Omega}\tag{1}\] Here \(\delta^{(r)}\) denotes the appropriate product of Dirac or Kronecker delta functions. We take the unconstrained reference measure to be independent across modes, \[\mathbb{W}(\{n_i\})=\prod_i W_i(n_i) \,. \label{eq:wf}\tag{2}\] Here \(n_i\in \mathcal{A}_i\), where \(\mathcal{A}_i\) is the set of allowed occupations of mode \(i\). This factorization is appropriate for an ideal-gas reference measure, where the unconstrained occupation weights contain no intermode correlations. The allowed occupations \(\mathcal{A}_i\) and the single-mode weights \(W_i(n_i)\) distinguish the MB, BE, and FD cases.

One-mode law and equilibrium statistics — Fixing the occupation \(n_k\) of one mode forces the remaining modes to carry \(Q-n_k q_k\). Therefore the exact one-mode conditional probability is \[\mathbb{P}_k(n_k\mid Q) = W_k(n_k) \frac{\Omega_{\neq k}(Q - n_k q_k)}{\Omega(Q)} \,. \label{eq:exact95one95mode}\tag{3}\] This exact identity is the starting point for the leading occupation laws. To evaluate the constrained densities, we introduce variables \(\chi_A\) conjugate to the conserved quantities \(Q^A\). The Laplace transform of \(\Omega(Q)\) is \[\widehat{\Omega}(\chi) = \int d^r Q\, e^{-\chi_{_A} Q^A}\Omega(Q). \label{eq:laplace95def}\tag{4}\] For discrete conserved quantities, this notation denotes the corresponding generating function, and the inverse transform below should be understood as coefficient extraction. Using the definition of \(\Omega(Q)\) and the factorized reference measure, Eqs. 1 and 2 , one obtains \[\widehat{\Omega}(\chi) = \prod_i z_i(\chi), \quad z_i(\chi) = \sum_{n_i \in \mathcal{A}_i} W_i(n_i)e^{-n_i \chi_{_A} q_i^A} \,. \label{eq:sm95z95general}\tag{5}\] The constrained density is recovered from the inverse transform \[\Omega(Q) = \int_\Gamma \frac{d^r\chi}{(2\pi i)^r} \exp\!\left[\Phi(\chi;Q)\right], \label{eq:inverse95laplace}\tag{6}\] where the contour \(\Gamma\) is chosen in the convergence domain of the transform, and \(\Phi(\chi;Q) = \chi_{_A} Q^A + \sum_i \log z_i(\chi)\). The constrained density is evaluated by expanding Eq. 6 around its stationary point \(\chi_{_A}^*\). The stationarity condition fixes \(\chi_{_A}^*\) through \[Q^A=\sum_i q_i^A \langle n_i\rangle_* \,, \label{eq:saddle95constraints}\tag{7}\] where \[\langle n_i\rangle_* = \frac{1}{z_i(\chi^*)} \sum_{n_i \in \mathcal{A}_i} n_i W_i(n_i) e^{-n_i \chi_{_A}^* q_i^A} \,. \label{eq:nbar95saddle}\tag{8}\] Thus the saddle variables are fixed by the imposed conserved quantities.

At leading saddle order, the numerator in Eq. 3 is expanded around the saddle fixed by the total conserved quantities \(Q\). The \(n_k\)-dependent part of the density-of-states ratio is \[\Omega_{\neq k}(Q - n_k q_k) \propto \frac{e^{-n_k \chi_{_A}^* q_k^A}}{z_k(\chi^*)} \,,\] where the proportionality constant is independent of \(n_k\). Therefore \[\mathbb{P}_k^{(0)}(n_k\mid Q) = \frac{1}{z_k(\chi^*)} W_k(n_k) e^{-n_k \chi_{_A}^* q_k^A} \,. \label{eq:leading95one95mode}\tag{9}\] This is the conditional equilibrium law. The conserved quantities determine the exponential tilt through \(\chi_{_A}^* q_k^A\), while the occupations \(\mathcal{A}_k\) and weights \(W_k(n_k)\) determine whether the resulting law is Maxwell–Boltzmann, Bose–Einstein, or Fermi–Dirac.

The stability of the saddle and the leading Gaussian response to conserved-quantity fluctuations are controlled by the Hessian of \(\Phi\), \[H^{AB} = \left. \frac{\partial^2\Phi}{\partial \chi_{_A} \partial \chi_{_B}} \right|_{\chi^*} = \sum_i q_i^A q_i^B \sigma_{i,*}^2 \,. \label{eq:H}\tag{10}\] Here \(\sigma_{i,*}^2 = \langle n_i^2 \rangle_* - \langle n_i\rangle_*^2\), evaluated with Eq. 9 . Thus \(H^{AB}\) is the susceptibility matrix of the conserved quantities in the leading product measure. A regular saddle requires \(H^{AB}\) to be nonsingular on the constrained subspace. Redundant conserved directions, or directions with zero fluctuations, should be removed before using the inverse susceptibility in the covariance projection below.

The usual equilibrium distributions follow directly from Eq. 9 . For four-momentum and particle-number constraints, \(Q=(P^\mu,N)\) and \(q_k=(p_k^\mu,1)\). Writing \(\chi_{_A}^*=(\beta_\mu^*,-\alpha^*)\), the saddle variable entering the one-mode law is \[x_k^* \equiv \chi_{_A}^* q_k^A = \beta_\mu^* p_k^\mu - \alpha^* \,.\] Here \(\beta_\mu^*\) and \(\alpha^*\) are the saddle values of the Laplace variables conjugate to \(P^\mu\) and \(N\). For an unbounded relativistic spectrum, convergence of the transform requires \(\beta_\mu p^\mu>0\) for all future-directed momenta, so \(\beta^\mu\) lies in the future-timelike domain. The leading conditional law becomes \[\mathbb{P}_k^{(0)}(n_k\mid P,N) = \frac{1}{z_k(x_k^*)} W_k(n_k) e^{-n_k x_k^*} \,. \label{eq:Pi095xk}\tag{11}\] Above equation shows for the first time that the relativistic Maxwell–Jüttner exponential form [30], [31] arises as a conditional tilt fixed by the imposed four-momentum and particle-number constraints, rather than by entropy maximization postulate.

For Maxwell–Boltzmann statistics [34], [35], any number of particles may occupy the same mode, \(\mathcal{A}_k=\{0,1,2,\ldots\}\), and the Gibbs factor removes overcounting by permutations, \(W_k(n_k)=1/n_k!\) [36]. Equation 11 then gives \[\langle n_k\rangle^{\rm MB}_* = \sum_{n_k=0}^{\infty} n_k\,\mathbb{P}_{k,{\rm MB}}^{(0)}(n_k \mid P,N) = e^{-x_k^*} \,. \label{eq:MB95result}\tag{12}\]

For Bose–Einstein statistics [37], [38], any number of identical bosons may occupy the same mode, \(\mathcal{A}_k=\{0,1,2,\ldots\}\), and each occupation number has unit weight, \(W_k(n_k)=1\). Equation 11 gives \[\langle n_k\rangle^{\rm BE}_* = \sum_{n_k=0}^{\infty} n_k\,\mathbb{P}_{k,{\rm BE}}^{(0)}(n_k \mid P,N) = \frac{1}{e^{x_k^*}-1} \,. \label{eq:BE95result}\tag{13}\]

For Fermi–Dirac statistics [39], [40], the Pauli exclusion rule restricts each mode to occupations \(\mathcal{A}_k=\{0,1\}\), with unit weight \(W_k(n_k)=1\). Equation 11 gives \[\langle n_k\rangle^{\rm FD}_* = \sum_{n_k=0}^{1} n_k\,\mathbb{P}_{k,{\rm FD}}^{(0)}(n_k \mid P,N) = \frac{1}{e^{x_k^*}+1} \,. \label{eq:FD95result}\tag{14}\]

For bosons, convergence requires \(x_k^*>0\) for each mode. In the local rest frame, \(\beta^{*\mu}=u^\mu/T\) and \(\alpha^*=\mu/T\), so \(x_k^*=(E_k-\mu)/T\). The boundary \(x_0\to0^+\) for the lowest mode corresponds to Bose accumulation [37], [38], [41]. In contrast, the FD occupation remains bounded and approaches the saturated Fermi surface as \(T\to 0\) [42]. Thus accumulation and saturation appear here as boundary behavior of the same conditional saddle.

Two-mode law and conservation-induced correlations — We next apply the same conditional construction to two modes. For two distinct modes \(k \neq \ell\), fixing \(n_k\) and \(n_\ell\) forces the remaining modes to carry \(Q - n_k q_k - n_\ell q_\ell\). The exact joint conditional probability is therefore \[\mathbb{P}_{k\ell}(n_k, n_\ell \mid Q) = W_k(n_k) W_\ell(n_\ell) \frac{\Omega_{\neq k,\ell} (Q-n_k q_k - n_\ell q_\ell)}{\Omega(Q)} \,. \label{eq:exact95two95mode}\tag{15}\] This is the two-mode analogue of Eq. 3 . The leading one-mode law factorizes, but the exact constraint couples the two modes because their conserved-quantity fluctuations must be compensated by the remaining system.

A saddle expansion of the density-of-states ratio in Eq. 15 gives a Gaussian cost for the conserved quantity carried by the deviations \(\Delta_i=n_i-\langle n_i\rangle_*\). The mixed part of this cost is \[-\Delta_k \Delta_\ell\, q_k^A (H^{-1})_{AB} \,q_\ell^B \,,\] where \(H^{AB}\) is the susceptibility matrix in Eq. 10 . Since the leading law is factorized, this mixed term gives the leading connected covariance for \(k \neq \ell\), \[\mathrm{Cov}_Q(n_k,n_\ell) = -\sigma_{k,*}^2 \sigma_{\ell,*}^2 q_k^A (H^{-1})_{AB}\, q_\ell^B +\cdots . \label{eq:offdiag95cov95main}\tag{16}\] The omitted terms include local one-mode corrections and higher-order terms in the saddle expansion. They do not change the leading off-diagonal covariance.

Combining this off-diagonal result with leading local variances gives the conditional-Gaussian projection form \[C_{ij}^{Q} = \sigma_{i,*}^2 \delta_{ij} - \sigma_{i,*}^2 \sigma_{j,*}^2 q_i^A (H^{-1})_{AB} \,q_j^B \,. \label{eq:C95projection}\tag{17}\] This form makes exact conservation manifest as \(\sum_i q_i^A C_{ij}^{Q}=0\). Thus fluctuations along the exactly conserved directions are projected out. Equivalently, the conditional covariance has no component along the conserved directions. The detailed derivation is given in the Supplemental Material.

The same covariance form can be contracted with observables that are linear in the mode occupations. Let \[X=\sum_i f_i n_i,\qquad Y=\sum_i g_i n_i \,,\] where \(f_i\) and \(g_i\) are fixed coefficients that specify the selected modes and their relative weights in each observable. Contracting Eq. 17 with these weights gives the conditional-Gaussian covariance \[\mathrm{Cov}^{\rm cg}_Q(X,Y) = \sum_i f_i g_i \sigma_{i,*}^2 - U_X^A (H^{-1})_{AB} U_Y^B , \label{eq:linear95proj95main}\tag{18}\] where \[U_X^A=\sum_i f_i\sigma_{i,*}^2 q_i^A,\qquad U_Y^A=\sum_i g_i\sigma_{i,*}^2 q_i^A \,. \label{eq:U95def}\tag{19}\] The vectors \(U_X^A\) and \(U_Y^A\) measure the overlap of the observables with the conserved quantities. If \(X\) and \(Y\) are built from disjoint sets of modes, the local diagonal term in Eq. 18 is absent, and the conservation-induced covariance reduces to \[\mathrm{Cov}^{\rm cons}_Q(X,Y) = -U_X^A(H^{-1})_{AB}U_Y^B \,. \label{eq:linear95cons95main}\tag{20}\] Thus exact conservation produces a finite-rank covariance determined by the conserved-quantity overlaps \(U_X^A\) and \(U_Y^A\), and the inverse susceptibility matrix.

The finite-rank form of Eq. 20 gives a direct way to remove this leading conservation component. Suppose the projection is performed on the same set of modes used to define \(H^{AB}\). For the observable \(X=\sum_i f_i n_i\), we define projected coefficients \[f_i^\perp = f_i - q_i^A (H^{-1})_{AB} U_X^B \,. \label{eq:fperp}\tag{21}\] The corresponding projected observable is \[X^\perp = \sum_i f_i^\perp n_i \,. \label{eq:Xperp}\tag{22}\] By construction, \(U_{X^\perp}^A=\sum_i f_i^\perp \sigma_{i,*}^2 q_i^A=0\). Equation 20 gives \(\mathrm{Cov}^{\rm cons}_Q(X^\perp,Y)=0\) for any disjoint linear observable \(Y\), at this leading order and for the selected conserved quantities. Thus the projection removes only the weight component aligned with the selected conserved quantities. Consequently, the leading covariance generated by the selected global conservation laws vanishes, while local statistical and genuine dynamical correlations can remain.

Conservation-orthogonal observables in high-energy collisions — We now specialize the projection to long-range correlation measurements in high-energy collisions. The aim is to construct modified analysis weights so that their overlap with a chosen set of exact conservation laws vanishes, and then compare the projected covariance with the ordinary covariance.

For an experimental analysis, a mode \(i\) can be represented by an analysis bin. Such a bin may be specified by particle species, \(p_T\), pseudorapidity, and azimuth, and \(n_i\) denotes the event-by-event number of particles in that bin. Particles in bin \(i\) are assigned conserved quantities \(q_i^A\), such as momentum components \(p_i^x, p_i^y\), energy \(E_i\), electric charge, and other conserved charges. In a fixed multiplicity or centrality class \(\mathcal{C}\), we write \[\bar{n}_i = \langle n_i\rangle_\mathcal{C},\qquad \delta n_i=n_i - \bar{n}_i \,,\] where \(\langle\cdots\rangle_\mathcal{C}\) denotes an event average over that class.

Let \(L\) and \(R\) be two non-overlapping rapidity-separated measured subevents in \(\mathcal{C}\). The mean multiplicity in left subevent is \(\overline{N}_L=\sum_{i\in L} \bar{n}_i\). We define the left-subevent fluctuation, \[X_L[f]=\frac{1}{\overline{N}_L}\sum_{i\in L} f_i\,\delta n_i \,. \label{eq:XL95HI}\tag{23}\] The coefficients \(f_i\) define the observable. For example, \(f_i=p_{T,i} - \bar{p}_{T,L}\), with \(\bar{p}_{T,L} = \sum_{i\in L}p_{T,i} \bar{n}_i / \sum_{i\in L} \bar{n}_i\), gives a mean-\(p_T\)-type fluctuation, while \(f_i=\cos(n\phi_i)\) or \(\sin(n\phi_i)\) gives a harmonic weight. The overlap of this left-subevent observable with the conserved quantities is \[U_L^A=\sum_{i\in L} f_i d_i q_i^A \,. \label{eq:UL95HI}\tag{24}\] Here \(d_i\) is the reference variance assigned to bin \(i\). It plays the role of \(\sigma_{i,*}^2\) in the formal projection 19 and determines how strongly bin \(i\) enters the conserved-quantity overlap. The independent-particle choice \(d_i=\bar n_i\) corresponds to a Poisson baseline for the bin occupation. In an experimental analysis, one may instead use measured or efficiency-corrected bin variances. The dependence on this choice can be used to estimate a systematic uncertainty of the projection. The right-subevent observable \(X_R[g]\) and overlap \(U_R^A\) are defined analogously.

Since the measured particles are only a subset of the full event, the leading covariance induced by exact conservation is governed by the full-event susceptibility \[H_{\rm full}^{AB}=\sum_{a\in{\rm full\;event}} d_a\, q_a^A q_a^B \,, \label{eq:Hfull95HI}\tag{25}\] which is generally inaccessible from finite-rapidity data. The leading global-conservation contribution between the two measured subevents has the form 20 \[C_{\rm cons}^{LR}[f,g]=-\frac{1}{\overline{N}_L \overline{N}_R}\,U_L^A (H_{\rm full}^{-1})_{AB} U_R^B \,. \label{eq:LR95cons95HI}\tag{26}\] Thus the unknown full-event susceptibility is required to predict the magnitude of the conservation-induced covariance, but not to construct a projected observable for which the leading conservation contribution vanishes. If the measured overlaps \(U_L^A\) and \(U_R^A\) vanish for the conserved quantities being projected out, the leading global-conservation contribution vanishes independently of the unknown \(H_{\rm full}^{AB}\).

We use only the measured subevent to define the weights \[S_L^{AB} = \sum_{i\in L} d_i q_i^A q_i^B \,. \label{eq:SL95HI}\tag{27}\] This matrix is not a full-event susceptibility and is used only to remove from \(f_i\) the component aligned with the selected conserved quantities in the measured left-subevent. If some selected conserved direction is redundant in this measured set, the inverse is taken after removing that direction. The projected left-subevent coefficients are \[f_i^\perp = f_i - q_i^A (S_L^{-1})_{AB} U_L^B \,, \qquad i\in L \,. \label{eq:fperp95HI}\tag{28}\] With this choice, \[U_{L,\perp}^A \equiv \sum_{i\in L} f_i^\perp d_i q_i^A=0 \,. \label{eq:Uperp95HI}\tag{29}\] Thus the projected observable has no overlap with the selected conserved quantities. Similarly, we define the right-subevent coefficients \(g_j^\perp\). The projected subevent fluctuations are \[X_L^\perp=\frac{1}{\overline{N}_L}\sum_{i\in L}f_i^\perp\,\delta n_i \,, \qquad X_R^\perp=\frac{1}{\overline{N}_R}\sum_{j\in R}g_j^\perp\,\delta n_j \,. \label{eq:Xperp95HI}\tag{30}\] The projected long-range covariance is \[C_\perp^{LR}[f,g] =\! \left\langle X_L^\perp[f]\,X_R^\perp[g] \right\rangle_{\mathcal{C}} \!\!= \frac{1}{\overline{N}_L\overline{N}_R} \sum_{i\in L}\sum_{j\in R} f_i^\perp g_j^\perp \!\left\langle \delta n_i\delta n_j\right\rangle_{\mathcal{C}}. \label{eq:Cperp95HI}\tag{31}\] Since \(U_{L,\perp}^A = U_{R,\perp}^A = 0\), the leading exact-conservation contribution from the selected conserved quantities vanishes in \(C_\perp^{LR}\). Note that the ordinary covariance, \[C^{LR}[f,g] = \mathrm{Cov}_\mathcal{C}\!\left(X_L[f], X_R[g]\right)\] contains the full measured long-range correlation, including any component aligned with the selected exact-conservation constraints. If exact conservation gives a significant contribution, the difference \(C^{LR} - C_\perp^{LR}\) should increase toward lower multiplicity, where the same conserved quantities are shared among fewer particles. A signal surviving the projection is free of leading exact-conservation backgrounds from the selected conserved quantities, while retaining dynamical correlations such as jets, resonance decays, and local charge conservation.

We show the implementation in Fig. 1, using \(10^5\) PYTHIA8/Angantyr [32], [33] p+Pb events. This analysis provides a closure test in a controlled event sample where the projection can be checked against transverse-momentum recoil and baryon-number conservation. The projection is applied separately in two rapidity-separated subevents, with \(d_i=\bar{n}_i\) and bins in \(p_T\) and \(\phi\). Projection against \(q_i^A=(p_{x,i},p_{y,i})\) removes the first-harmonic covariance, as expected from transverse-momentum recoil, while leaving the second harmonic and mean-\(p_T\) covariance unchanged. For the baryon-number observable \(f_i=B_i\), projecting only \((p_x,p_y)\) leaves the covariance unchanged, while including \(B_i\) among the projected quantities removes it. Thus projection removes correlations aligned with selected conservation laws, while leaving conservation-orthogonal long-range components unchanged.

Discussion — The conditional framework presented here separates two aspects in equilibrium statistics. Conditional-limit reasoning explains why additive constraints produce an exponential tilt [25], [28]. The microscopic occupation measure then determines the statistics. Poisson weights give MB statistics, unrestricted unit weights give BE statistics, and the exclusion constraint gives FD statistics. The Gaussian expansion around the same conditional saddle gives the susceptibility matrix that controls finite-size fluctuations and conservation-induced correlations.

The same framework turns exact conservation laws into projection variables. In high-energy collisions, the projection separates the component of a long-range correlation aligned with selected exact conservation laws from the conservation-orthogonal component. This is useful in small systems and peripheral collisions, where conservation effects are less suppressed by multiplicity. The construction can also be applied to conserved-charge fluctuation measurements near the QCD critical point, where global charge conservation must be separated from fluctuations generated by critical dynamics [43], [44]. Such a separation is crucial for the reliable identification of critical-point signatures. Few-body ultracold-atom expansions provide another setting where collective-like behavior has been reported in systems with very small particle number [45], [46]. In such systems, correlations implied by fixed particle number, conserved energy–momentum, and fixed internal-state populations can be comparable to dynamical correlations. Applying the projection would test which part of the observed collective-like response survives after the leading exact-conservation component is removed. The framework developed here opens new avenues for identifying genuine collective dynamics in finite many-body systems.

Acknowledgments — We thank Richard Furnstahl, Chun Shen, Soumangsu Chakraborty, Trakshu Sharma, and Chandrodoy Chattopadhyay for engaging discussions. S.J. gratefully acknowledges support from the CSSI program under Award No. OAC-2004601 (BAND Collaboration [47]). S.J. acknowledges the warm hospitality of NISER Bhubaneswar where this work was initiated. A.J. acknowledges Department of Atomic Energy (DAE), India for financial support.

We present the derivation of the equations appearing in the main text.

1 Notation↩︎

We denote the exactly conserved quantities by an \(r\)-component vector \(Q\), with components \[Q^A=\sum_i q_i^A n_i, \qquad A=1,\ldots,r .\] Here \(n_i\) is the occupation number in mode \(i\), and \(q_i^A\) is the \(A\)th component of the corresponding single-particle conserved-charge vector \(q_i\). Capital indices \(A,B,\ldots\) label the space of conserved quantities and are summed over when repeated: \(\chi_{_A} q_i^A = \sum_{A=1}^r \chi_{_A} q_i^A\). Mode labels \(i,j,\ldots\) are summed or multiplied over only when an explicit \(\sum_i\) or \(\prod_i\) is written. They should not be confused with Lorentz indices, which are denoted by Greek letters \(\mu,\nu,\ldots\).

For the four-momentum and particle-number case used in the main text, we may take \[Q=(P^\mu, N), \qquad q_i=(p_i^\mu, 1).\] The last component of \(q_i\) is unity because each particle occupying mode \(i\) contributes one unit to the total particle number \(N\). The corresponding conjugate variables are \(\chi_{_A}=(\beta_\mu, -\alpha)\), so that \[\chi_{_A} Q^A = \beta_\mu P^\mu - \alpha N, \qquad \chi_{_A} q_i^A = \beta_\mu p_i^\mu - \alpha .\] For a hadron gas with exact four-momentum, baryon number, strangeness, and other conserved charges, one may instead take \(Q=(P^\mu, B, S, \ldots)\) with \(q_i=(p_i^\mu, B_i, S_i, \ldots)\). The signs of the corresponding chemical potentials are absorbed into the definition of \(\chi_{_A}\).

2 Equilibrium statistics as conditional limit laws↩︎

We consider occupation configurations \(\{n_i\}\), with allowed occupations \(n_i \in \mathcal{A}_i\). A configuration \(\{n_i\}\) is counted with a weight \(W(\{n_i\})\), which encodes the microscopic counting. The constrained density of configurations at fixed \(Q\) is \[\Omega(Q) = \sum_{\{n_i\}} \mathbb{W}(\{n_i\})\, \delta^{(r)} \left( \sum_i q_i n_i-Q \right). \label{eqA:Omega95Q}\tag{32}\] The notation \(\delta^{(r)}\) denotes the appropriate product of Dirac or Kronecker delta functions, depending on whether the corresponding conserved quantities are continuous or discrete. For example, if \(Q=(P^\mu,N)\) and \(q_i=(p_i^\mu,1)\), then \[\delta^{(r)}\left(\sum_i q_i n_i-Q\right) = \delta^{(4)}\left(\sum_i p_i^\mu n_i-P^\mu\right) \delta_{\sum_i n_i,N}.\] For the Maxwell–Boltzmann (MB), Bose–Einstein (BE), and Fermi–Dirac (FD) counting rules considered in this work, the unconstrained microscopic counting assigns a separate weight to each mode. Hence \[\mathbb{W}(\{n_i\}) = \prod_i W_i(n_i). \label{eqA:weight95factorized}\tag{33}\]

2.1 Exact single-mode conditional distribution↩︎

The factorization in Eq. 33 allows us to isolate the exact conditional distribution of a single mode. If mode \(k\) has occupation \(n_k \in \mathcal{A}_k\), the remaining modes must carry the conserved quantity \(Q - n_k q_k\). Therefore the exact conditional probability for the occupation of mode \(k\) is \[\mathbb{P}_k(n_k \mid Q) = W_k(n_k) \frac{\Omega_{\neq k}(Q-n_k q_k)}{\Omega(Q)} . \label{eqA:exact95marginal95Q}\tag{34}\] Here \(\Omega_{\neq k}\) denotes the constrained density of configurations of all modes except \(k\). For the four-momentum and particle-number case, \(Q=(P^\mu,N)\) and \(q_k=(p_k^\mu,1)\), Eq. 34 becomes \[\mathbb{P}_k(n_k\mid P,N) = W_k(n_k) \frac{\Omega_{\neq k}(P-n_k p_k,N-n_k)}{\Omega(P,N)} .\]

2.2 Laplace transform and conditional saddle↩︎

To evaluate the constrained densities in the large-system limit, we introduce variables \(\chi_{_A}\) conjugate to the conserved quantities \(Q^A\). We define \[\widehat{\Omega}(\chi) = \int d^r Q\, e^{-\chi_{_A} Q^A}\Omega(Q). \label{eqA:Laplace95def}\tag{35}\] For conserved quantities with discrete components, this notation should be read as the corresponding generating function. The inverse transform then amounts to coefficient extraction in those discrete variables.

Using Eqs. 32 and 33 , we obtain \[\begin{align} \widehat\Omega(\chi) &= \sum_{\{n_i\}} \left[ \prod_i W_i(n_i)\right] \int d^r Q\, e^{-\chi_{_A} Q^A} \delta^{(r)} \left( \sum_i q_i n_i-Q \right) \nonumber\\ &= \sum_{\{n_i\}} \left[ \prod_i W_i(n_i)\right] \exp \left( -\chi_{_A}\sum_i q_i^A n_i \right) \end{align}\] Since \[\exp \left( -\chi_{_A}\sum_i q_i^A n_i \right) = \prod_i \exp\left(-n_i\chi_{_A} q_i^A\right),\] we have \[\begin{align} \widehat\Omega(\chi) = \sum_{\{n_i\}} \prod_i \left[W_i(n_i)\, \exp\left(-n_i\chi_{_A} q_i^A\right)\right] = \prod_i \sum_{n_i\in\mathcal{A}_i} W_i(n_i) \, \exp\left(-n_i\chi_{_A} q_i^A\right) = \prod_i z_i(\chi) \,. \label{eqA:Laplace95factorization} \end{align}\tag{36}\] Here we have defined the single-mode generating function \[z_i(\chi) = \sum_{n_i\in\mathcal{A}_i} W_i(n_i)\, e^{-n_i\chi_{_A} q_i^A} \,. \label{eqA:single95mode95z}\tag{37}\] The constrained density is recovered by the inverse transform \[\Omega(Q) = \int_{\Gamma} \frac{d^r\chi}{(2\pi i)^r} \exp\left[\Phi(\chi;Q)\right], \label{eqA:inverse95transform}\tag{38}\] with \[\Phi(\chi;Q) = \chi_{_A} Q^A + \sum_i \log z_i(\chi). \label{eqA:Phi95def}\tag{39}\] The contour \(\Gamma\) is a product of Bromwich contours chosen within the convergence domain of the transform.

The constrained density is obtained by a saddle expansion of the inverse-transform integral. The dominant contribution comes from the stationary point \(\chi_{_A}^*\) of the exponent \(\Phi(\chi;Q)\), whose saddle condition is \[0= \left. \frac{\partial\Phi}{\partial\chi_{_A}} \right|_{\chi^*} = Q^A + \sum_i \left. \frac{\partial\log z_i}{\partial\chi_{_A}} \right|_{\chi^*}. \label{eqA:saddle95cond}\tag{40}\] Using Eq. 37 , \[\frac{\partial\log z_i}{\partial\chi_{_A}} = -q_i^A \frac{\sum_{n_i\in\mathcal{A}_i} n_i W_i(n_i)e^{-n_i \chi_{_C} q_i^C}}{\sum_{n_i\in\mathcal{A}_i} W_i(n_i) e^{-n_i \chi_{_C} q_i^C}} \,. \label{eqA:dlogz}\tag{41}\] Thus the saddle equation becomes \[Q^A=\sum_i q_i^A \langle n_i\rangle_* , \label{eqA:saddle95Q}\tag{42}\] where \[\langle n_i\rangle_* = \frac{1}{z_i(\chi^*)} \sum_{n_i\in\mathcal{A}_i} n_i W_i(n_i)e^{-n_i \chi_{_A}^* q_i^A} \,.\]

2.3 Leading occupation law↩︎

The leading one-mode law follows from the saddle expansion of the exact density-of-states ratio in Eq. 34 . For fixed \(n_k\in\mathcal{A}_k\), the numerator is evaluated at the shifted conserved quantity \(Q-n_k q_k\). At leading saddle order, we evaluate the numerator and denominator at the same stationary point \(\chi^*\).

The numerator is then approximated by \[\Omega_{\neq k}(Q-n_k q_k) \simeq \exp\left[ \Phi_{\neq k}(\chi^*; Q-n_k q_k) \right],\] where \[\begin{align} \Phi_{\neq k}(\chi^*; Q-n_k q_k) &= \chi_{_A}^* Q^A + \sum_{i\neq k}\log z_i(\chi^*) - n_k \chi_{_A}^* q_k^A = \chi_{_A}^* Q^A + \sum_{i}\log z_i(\chi^*) - \log z_k(\chi^*) - n_k \chi_{_A}^* q_k^A \nonumber\\ &= \Phi(\chi^*;Q) - \log z_k(\chi^*) - n_k \chi_{_A}^* q_k^A . \end{align}\] Similarly, the denominator is approximated by \[\Omega(Q) \simeq \exp\left[\Phi(\chi^*;Q)\right].\] Substituting these into Eq. 34 , we obtain \[\mathbb{P}_k^{(0)}(n_k \mid Q) = W_k(n_k) \left. \frac{\Omega_{\neq k}(Q-n_k q_k)}{\Omega(Q)} \right|_{\chi^*} = \frac{1}{z_k(\chi^*)} W_k(n_k) e^{-n_k\chi_{_A}^* q_k^A}. \label{eqA:Pk095def}\tag{43}\] The exponential factor is determined by the saddle values \(\chi_{_A}^*\), which are fixed by the imposed conserved quantities through Eq. 42 . For the four-momentum and particle-number case, \(Q=(P^\mu,N)\) and \(q_k=(p_k^\mu,1)\), we write \(\chi_{_A}^*=(\beta_\mu^*,-\alpha^*)\). Here \(\beta_\mu\) and \(\alpha\) are Laplace variables conjugate to \(P^\mu\) and \(N\). For an unbounded relativistic spectrum, convergence requires \(\beta_\mu p^\mu>0\) for all future-directed on-shell momenta \(p^\mu\), which places \(\beta^\mu\) in the future-timelike domain.

2.4 Hessian and saddle stability↩︎

The stability of the conditional saddle is determined by the quadratic variation of the exponent around the stationary point. We write \[\chi_{_A} = \chi_{_A}^* + \delta \chi_{_A} .\] Expanding the exponent about the stationary point gives \[\Phi(\chi;Q) = \Phi(\chi^*;Q) + \frac{1}{2}\, \delta\chi_{_A} H^{AB}\delta\chi_{_B} +\cdots , \label{eqA:Phi95gaussian95general}\tag{44}\] where the linear term is absent by the saddle equation 40 , and \[H^{AB} \equiv \left. \frac{\partial^2\Phi}{\partial\chi_{_A} \partial\chi_{_B}} \right|_{\chi^*} \,, \label{eqA:Hessian95def95general}\tag{45}\] is the Hessian of the exponent at the saddle. Since \(\Phi(\chi;Q)=\chi_{_A}Q^A+\sum_i\log z_i(\chi)\), \[\begin{align} H^{AB} &= \left.\frac{\partial}{\partial \chi_{_B}} \left(\frac{\partial \Phi}{\partial\chi_{_A}} \right) \right|_{\chi^*} = \left.\frac{\partial}{\partial \chi_{_B}} \left(Q^A + \sum_i \frac{\partial \log z_i}{\partial\chi_{_A}} \right) \right|_{\chi^*} = \left. \sum_i \frac{\partial^2 \log z_i}{\partial\chi_{_A} \partial\chi_{_B}} \right|_{\chi^*} \end{align}\] Using Eq. 37 , the first derivative is \[\frac{\partial \log z_i}{\partial \chi_{_A}} = -q_i^A \langle n_i\rangle ,\] where \[\langle n_i\rangle = \frac{1}{z_i(\chi)} \sum_{n_i\in\mathcal{A}_i} n_i W_i(n_i)e^{-n_i\chi_{_C} q_i^C}.\] Taking one more derivative gives \[\begin{align} \frac{\partial^2 \log z_i}{\partial \chi_{_A} \partial\chi_{_B}} &= -q_i^A \frac{\partial \langle n_i\rangle}{\partial\chi_{_B}} = -q_i^A \frac{\partial}{\partial\chi_{_B}} \left[ \frac{ \sum_{n_i\in\mathcal{A}_i} n_i W_i(n_i)e^{-n_i\chi_{_C} q_i^C} }{ z_i(\chi) } \right] \nonumber\\ &= -q_i^A \left[ \frac{1}{z_i(\chi)} \frac{\partial}{\partial\chi_{_B}} \sum_{n_i\in\mathcal{A}_i} n_i W_i(n_i)e^{-n_i\chi_{_C} q_i^C} - \frac{ \sum_{n_i\in\mathcal{A}_i} n_i W_i(n_i)e^{-n_i\chi_{_C} q_i^C} }{ z_i(\chi)^2 } \frac{\partial z_i(\chi)}{\partial\chi_{_B}} \right] \nonumber\\ &= -q_i^A \left[ -q_i^B \langle n_i^2\rangle + q_i^B \langle n_i\rangle^2 \right] = q_i^A q_i^B \left( \langle n_i^2\rangle-\langle n_i\rangle^2 \right). \end{align}\] Therefore, at the saddle, \[H^{AB} = \left. \frac{\partial^2\Phi}{\partial\chi_{_A}\partial\chi_{_B}} \right|_{\chi^*} = \sum_i q_i^A q_i^B \sigma^2_{i,*}, \label{eqA:Hessian95general}\tag{46}\] where \[\sigma^2_{i,*} = \langle n_i^2\rangle_* - \langle n_i\rangle_*^2\] is the variance of the leading one-mode law in Eq. 43 .

Equivalently, we can define the fluctuation of the conserved quantities in the leading product measure by \[\delta Q^A = \sum_i q_i^A \left( n_i-\langle n_i\rangle_* \right),\] Then \[H^{AB} = \left\langle \delta Q^A \delta Q^B \right\rangle_* . \label{eqA:Hessian95as95covariance}\tag{47}\] Thus \(H^{AB}\) is the covariance matrix of the conserved quantities before the exact constraint is imposed. The inverse of \(H^{AB}\) exists only if the imposed conserved quantities are independent and have nonzero fluctuations. In other words, for every nonzero vector \(v_A\) in the conserved-charge space, \[v_AH^{AB}v_B = \sum_i \sigma_{i,*}^2 \left( v_A q_i^A \right)^2 > 0 . \label{eqA:Hessian95positive95condition}\tag{48}\] This condition is the saddle-stability condition in the constrained directions. It fails if one of the constraints is redundant, or if some imposed conserved combination does not fluctuate in the leading product measure. In that case the corresponding direction must be removed before using the projection formula.

When Eq. 48 holds, we define the inverse Hessian as \[H^{AC}(H^{-1})_{CB} = \delta^A_{\;B}. \label{eqA:H95inverse95def}\tag{49}\] The position of capital indices is not tensorial. \(H^{AB}\) and \((H^{-1})_{AB}\) denote a matrix and its inverse.

2.5 Equilibrium statistics↩︎

For the four-momentum and particle-number constraint, we take \[Q=(P^\mu,N), \qquad q_k=(p_k^\mu,1), \qquad \chi_{_A}^*=(\beta_\mu^*,-\alpha^*) .\] The leading one-mode law given in Eq. 43 becomes \[\mathbb{P}_k^{(0)}(n_k\mid P,N) = \frac{1}{z_k(\beta^*,\alpha^*)} W_k(n_k) \exp\left[-n_k(\beta_\mu^*p_k^\mu-\alpha^*)\right]. \label{eqA:Pk095PN}\tag{50}\] Defining \[x_k \equiv \beta_\mu^*p_k^\mu-\alpha^* ,\] the single-mode generating function is \[z_k(\beta^*,\alpha^*) = \sum_{n_k\in\mathcal{A}_k} W_k(n_k)e^{-n_k x_k}.\]

For Maxwell–Boltzmann counting, \(\mathcal{A}_k=\{0,1,2,\ldots\}\) and \(W_k(n_k)=1/n_k!\). Hence \[\begin{align} \langle n_k\rangle^{\rm MB}_* &= \sum_{n_k=0}^{\infty} n_k\,\mathbb{P}_{k,{\rm MB}}^{(0)}(n_k \mid P,N) = \frac{\sum_{n_k=0}^{\infty} n_k\,e^{-n_k x_k}/n_k!}{\sum_{n_k=0}^{\infty} e^{-n_k x_k}/n_k!} = e^{-x_k} \,. \end{align}\]

For Bose–Einstein counting, \(\mathcal{A}_k=\{0,1,2,\ldots\}\) and \(W_k(n_k)=1\). Therefore \[\begin{align} \langle n_k\rangle^{\rm BE}_* &= \sum_{n_k=0}^{\infty} n_k\,\mathbb{P}_{k,{\rm BE}}^{(0)}(n_k \mid P,N) = \frac{\sum_{n_k=0}^{\infty} n_k\,e^{-n_k x_k}}{\sum_{n_k=0}^{\infty} e^{-n_k x_k}} = \frac{1}{e^{x_k}-1} \,. \end{align}\] Note that the sums converge only for \(x_k>0\).

For Fermi–Dirac counting, \(\mathcal{A}_k=\{0,1\}\) and \(W_k(n_k)=1\). Thus \[\begin{align} \langle n_k\rangle^{\rm FD}_* &= \sum_{n_k=0}^{1} n_k\,\mathbb{P}_{k,{\rm FD}}^{(0)}(n_k \mid P,N) = \frac{\sum_{n_k=0}^{1} n_k\,e^{-n_k x_k}}{\sum_{n_k=0}^{1} e^{-n_k x_k}} = \frac{1}{e^{x_k}+1} \,. \end{align}\]

3 Two-mode conditional law and induced covariance↩︎

We now derive the leading correlation between two modes induced solely by the exact conservation constraint. Consider two distinct modes \(k\neq \ell\), with allowed occupations \(n_k \in \mathcal{A}_k\) and \(n_\ell \in \mathcal{A}_\ell\). If these occupations are fixed, the remaining modes must carry the conserved quantity \(Q-n_k q_k-n_\ell q_\ell\). Using the factorized microscopic weight, the exact joint conditional probability is \[\mathbb{P}_{k\ell}(n_k,n_\ell\mid Q) = W_k(n_k) W_\ell(n_\ell) \frac{\Omega_{\neq k,\ell}(Q-n_k q_k-n_\ell q_\ell)}{\Omega(Q) } . \label{eqA:two95mode95exact}\tag{51}\] Here \(\Omega_{\neq k,\ell}\) denotes the constrained density of configurations of all modes except \(k\) and \(\ell\).

3.1 Saddle expansion of the two-mode law↩︎

We now expand both the numerator and denominator in Eq. 51 around the same full-system saddle \(\chi^*\). The denominator is written as the inverse transform \[\Omega(Q) = \int_{\Gamma}\frac{d^r\chi}{(2\pi i)^r} \exp\left[\Phi(\chi;Q)\right], \label{eqA:TM95denominator95inverse}\tag{52}\] with \[\Phi(\chi;Q) = \chi_{_A}Q^A+\sum_i\log z_i(\chi). \label{eqA:TM95Phi95denominator}\tag{53}\] Introducing fluctuations around the saddle, \[\chi_{_A}=\chi_{_A}^* + \delta\chi_{_A} ,\] the denominator exponent expands as \[\Phi(\chi;Q) = \Phi(\chi^*;Q) + \frac{1}{2}\delta\chi_{_A} H^{AB}\delta\chi_{_B} +\cdots , \label{eqA:TM95denominator95Taylor}\tag{54}\] where the linear term is absent by the saddle equation, and \(H^{AB}\) is the Hessian defined in Eq. 46 . Near the saddle, the inverse-transform contour is locally parametrized by \[\delta\chi_{_A}=i y_A ,\] with real \(y_A\). The denominator exponent then becomes \[\Phi(\chi;Q) = \Phi(\chi^*;Q) - \frac{1}{2}y_AH^{AB}y_B+\cdots . \label{eqA:TM95denominator95gaussian95y}\tag{55}\] Thus, to Gaussian order, \[\Omega(Q) \simeq e^{\Phi(\chi^*;Q)} \int \frac{d^r y}{(2\pi)^r} \exp\left[ -\frac{1}{2}y_AH^{AB}y_B \right] \left[ 1 + \cdots \right] = e^{\Phi(\chi^*;Q)} \frac{\mathcal{C}}{\sqrt{\det H}} \left[ 1 + \cdots \right]. \label{eqA:TM95denominator95gaussian95integral}\tag{56}\] Here \(\mathcal{C}\) denotes a source-independent normalization factor determined by the Gaussian integration measure and contour convention. Here and below, the ellipsis denotes non-Gaussian saddle corrections from cubic and higher derivatives of the exponent.

For the numerator, we define the shifted conserved quantity \[Q'^A = Q^A - n_k q_k^A - n_\ell q_\ell^A . \label{eqA:TM95Qprime95def}\tag{57}\] Then \[\Omega_{\neq k,\ell}(Q') = \int_{\Gamma}\frac{d^r\chi}{(2\pi i)^r} \exp\left[\Phi_{\neq k,\ell}(\chi;Q')\right], \label{eqA:TM95numerator95inverse}\tag{58}\] where \[\Phi_{\neq k,\ell}(\chi;Q') = \chi_{_A} Q'^A + \sum_{i\neq k,\ell}\log z_i(\chi). \label{eqA:TM95Phi95numerator}\tag{59}\] At the full-system saddle, \[\begin{align} \Phi_{\neq k,\ell}(\chi^*;Q') &= \chi_{_A}^* \left( Q^A-n_k q_k^A-n_\ell q_\ell^A \right) + \sum_{i\neq k,\ell}\log z_i(\chi^*) \nonumber\\ &= \Phi(\chi^*;Q) - \log z_k(\chi^*) - \log z_\ell(\chi^*) - n_k\chi_{_A}^*q_k^A - n_\ell\chi_{_A}^*q_\ell^A . \label{eqA:TM95reduced95Phi95at95saddle} \end{align}\tag{60}\] Unlike the denominator exponent, the reduced exponent is not stationary at \(\chi^*\) because the occupations of modes \(k\) and \(\ell\) have been fixed. Its first derivative at \(\chi^*\) is \[\begin{align} \left. \frac{\partial\Phi_{\neq k,\ell}}{\partial\chi_{_A}} \right|_{\chi^*} &= Q'^A + \sum_{i\neq k,\ell} \left. \frac{\partial\log z_i}{\partial\chi_{_A}} \right|_{\chi^*} \nonumber\\ &= Q^A - n_k q_k^A - n_\ell q_\ell^A - \sum_{i\neq k,\ell} q_i^A\langle n_i\rangle_* \nonumber\\ &= Q^A - \sum_i q_i^A\langle n_i\rangle_* - (n_k-\langle n_k\rangle_*)q_k^A - (n_\ell-\langle n_\ell\rangle_*)q_\ell^A . \end{align}\] Using the saddle equation 42 , the first two terms cancel. Therefore \[\left. \frac{\partial\Phi_{\neq k,\ell}}{\partial\chi_{_A}} \right|_{\chi^*} = - (n_k-\langle n_k\rangle_*)q_k^A - (n_\ell-\langle n_\ell\rangle_*)q_\ell^A . \label{eqA:TM95reduced95linear}\tag{61}\] It is useful to define \[R_{k\ell}^A = (n_k-\langle n_k\rangle_*)q_k^A + (n_\ell-\langle n_\ell\rangle_*)q_\ell^A . \label{eqA:Rkl95def}\tag{62}\] Then the numerator exponent expands as \[\Phi_{\neq k,\ell}(\chi;Q') = \Phi_{\neq k,\ell}(\chi^*;Q') - \delta\chi_{_A} R_{k\ell}^A + \frac{1}{2}\delta\chi_{_A} H^{AB}\delta\chi_{_B} +\cdots . \label{eqA:TM95numerator95Taylor}\tag{63}\] Using the same local contour parametrization \(\delta\chi_{_A}=i y_A\), this becomes \[\Phi_{\neq k,\ell}(\chi;Q') = \Phi_{\neq k,\ell}(\chi^*;Q') - i y_A R_{k\ell}^A - \frac{1}{2}y_AH^{AB}y_B +\cdots . \label{eqA:TM95numerator95gaussian95y}\tag{64}\] Thus, to Gaussian order, \[\begin{align} \Omega_{\neq k,\ell}(Q') &\simeq e^{\Phi_{\neq k,\ell}(\chi^*;Q')} \int \frac{d^r y}{(2\pi)^r} \exp\left[ -\frac{1}{2}y_AH^{AB}y_B - i y_A R_{k\ell}^A \right] \left[ 1 + \cdots \right] \nonumber\\ &= e^{\Phi_{\neq k,\ell}(\chi^*;Q')} \frac{\mathcal{C}}{\sqrt{\det H}} \exp\left[ -\frac{1}{2} R_{k\ell}^A(H^{-1})_{AB}R_{k\ell}^B \right] \left[ 1 + \cdots \right]. \label{eqA:TM95numerator95gaussian95integral} \end{align}\tag{65}\] In obtaining the second line, we used the standard Gaussian identity \[\int d^r u\, \exp\left[ -\frac{1}{2}u_AH^{AB}u_B - iJ^A u_A \right] = \frac{\mathcal{C}}{\sqrt{\det H}} \exp\left[ -\frac{1}{2}J^A(H^{-1})_{AB}J^B \right], \label{eqA:gaussian95identity}\tag{66}\] with \(J^A=R_{k\ell}^A\).

Dividing Eq. 65 by Eq. 56 , the leading Gaussian normalization cancels in the ratio. The residual determinant ratio from \(H_{\neq k,\ell}\) is independent of \(n_k,n_\ell\) at the order relevant for the mixed covariance and is absorbed into the normalization of the two mode distribution. We obtain \[\frac{\Omega_{\neq k,\ell}(Q')}{\Omega(Q)} \simeq \exp\left[ \Phi_{\neq k,\ell}(\chi^*;Q') - \Phi(\chi^*;Q) \right] \exp\left[ -\frac{1}{2} R_{k\ell}^A(H^{-1})_{AB}R_{k\ell}^B \right] \left[ 1 + \cdots \right].\] Using Eq. 60 , this becomes \[\frac{\Omega_{\neq k,\ell}(Q')}{\Omega(Q)} \simeq \frac{e^{-n_k\chi_{_A}^*q_k^A}}{ z_k(\chi^*)} \frac{e^{-n_\ell\chi_{_A}^*q_\ell^A}}{z_\ell(\chi^*)} \exp\left[ -\frac{1}{2} R_{k\ell}^A(H^{-1})_{AB}R_{k\ell}^B \right] \left[ 1 + \cdots \right]. \label{eqA:TM95density95ratio95final}\tag{67}\] Substituting this into the exact two-mode law, Eq. 51 , gives \[\begin{align} \mathbb{P}_{k\ell}(n_k,n_\ell\mid Q) &= \frac{1}{\mathcal{Z}_{k\ell}} \mathbb{P}_k^{(0)}(n_k\mid Q) \mathbb{P}_\ell^{(0)}(n_\ell\mid Q) \exp\left[ -\frac{1}{2} R_{k\ell}^A(H^{-1})_{AB}R_{k\ell}^B \right] \left[ 1 + \cdots \right]. \label{eqA:TM95joint95gaussian} \end{align}\tag{68}\] where \(\mathcal{Z}_{k\ell}\) normalizes the two mode distribution.

3.2 Conservation-induced covariance projection↩︎

Equation 68 is the conditional-Gaussian form of the two-mode law. Corrections beyond this form can modify local one-mode terms and higher cumulants. The leading connected correlation between two distinct modes is obtained from the mixed part of the quadratic compensation term.

Defining the occupation deviations \[\Delta_k = n_k-\langle n_k\rangle_*, \qquad \Delta_\ell = n_\ell-\langle n_\ell\rangle_* \,, \label{eqA:Delta95defs}\tag{69}\] Eq. 62 becomes \[R_{k\ell}^A = q_k^A\Delta_k+q_\ell^A\Delta_\ell . \label{eqA:Rkl95Delta}\tag{70}\] We also define the shorthand \[\Lambda_{ij} = q_i^A(H^{-1})_{AB}q_j^B \,. \label{eqA:Lambda95def}\tag{71}\] Then the quadratic form in Eq. 68 is \[\begin{align} \frac{1}{2} R_{k\ell}^A(H^{-1})_{AB}R_{k\ell}^B &= \frac{1}{2}\Lambda_{kk}\Delta_k^2 + \frac{1}{2}\Lambda_{\ell\ell}\Delta_\ell^2 + \Lambda_{k\ell}\Delta_k\Delta_\ell . \label{eqA:compensation95expanded} \end{align}\tag{72}\] The first two terms depend only on a single occupation variable and can be absorbed into local one-mode corrections and the normalization \(\mathcal{Z}_{k\ell}\). Edgeworth-type corrections from cubic and higher derivatives of the saddle exponent can also modify local terms and higher cumulants. They do not change the leading mixed Gaussian compensation term proportional to \(\Lambda_{k\ell}\Delta_k\Delta_\ell\).

Expanding Eq. 68 to the order needed for the leading off-diagonal covariance, and retaining only the mixed term between \(\Delta_k\) and \(\Delta_\ell\), gives \[\begin{align} \mathbb{P}_{k\ell}(n_k,n_\ell\mid Q) &\simeq \mathbb{P}_k^{(0)}(n_k\mid Q) \mathbb{P}_\ell^{(0)}(n_\ell\mid Q) \left[ 1-\Lambda_{k\ell}\Delta_k\Delta_\ell+\cdots \right]. \label{eqA:two95mode95mixed95expanded} \end{align}\tag{73}\] The omitted terms are either local in one of the two modes or beyond the conditional-Gaussian contribution relevant for the leading off-diagonal correlation.

We denote expectation with respect to the conditional two-mode law \(\mathbb{P}_{k\ell}(n_k,n_\ell\mid Q)\) by \(\langle\cdots\rangle_Q\). The conditional covariance between two distinct modes is \[\mathrm{Cov}_Q(n_k,n_\ell) = \langle n_k n_\ell\rangle_Q - \langle n_k\rangle_Q\langle n_\ell\rangle_Q, \qquad k\neq \ell . \label{eqA:cov95def}\tag{74}\] Since \(\Delta_k\) and \(\Delta_\ell\) differ from \(n_k\) and \(n_\ell\) only by constants, the same covariance can be written as \[\mathrm{Cov}_Q(n_k,n_\ell) = \langle \Delta_k\Delta_\ell\rangle_Q - \langle \Delta_k\rangle_Q \langle \Delta_\ell\rangle_Q . \label{eqA:cov95def95Delta}\tag{75}\]

Let \(\langle\cdots\rangle_0\) denote expectation with respect to the factorized leading law \[\mathbb{P}_k^{(0)}(n_k\mid Q)\, \mathbb{P}_\ell^{(0)}(n_\ell\mid Q).\] In particular, \[\langle \Delta_k\rangle_0=0, \qquad \langle \Delta_\ell\rangle_0=0, \qquad \langle \Delta_k^2\rangle_0=\sigma_{k,*}^2, \qquad \langle \Delta_\ell^2\rangle_0=\sigma_{\ell,*}^2 . \label{eqA:leading95Delta95moments}\tag{76}\] Using Eq. 73 , the one-mode shift of mode \(k\) is \[\begin{align} \langle \Delta_k\rangle_Q &= \sum_{n_k\in\mathcal{A}_k} \sum_{n_\ell\in\mathcal{A}_\ell} \Delta_k\, \mathbb{P}_k^{(0)}(n_k\mid Q) \mathbb{P}_\ell^{(0)}(n_\ell\mid Q) \left[ 1-\Lambda_{k\ell}\Delta_k\Delta_\ell+\cdots \right] \nonumber\\ &= \langle \Delta_k\rangle_0 - \Lambda_{k\ell} \langle \Delta_k^2\rangle_0 \langle \Delta_\ell\rangle_0 +\cdots = 0+\cdots . \label{eqA:Delta95k95shift} \end{align}\tag{77}\] Similarly, \(\langle \Delta_\ell\rangle_Q=0+\cdots\). Thus the product \(\langle \Delta_k\rangle_Q\langle \Delta_\ell\rangle_Q\) does not contribute to the leading off-diagonal covariance in Eq. 75 .

The only term in Eq. 73 that couples the two modes is the mixed term proportional to \(\Lambda_{k\ell}\Delta_k\Delta_\ell\). Its contribution to the remaining factor in Eq. 75 is \[\begin{align} \langle \Delta_k\Delta_\ell\rangle_Q &= \sum_{n_k\in\mathcal{A}_k} \sum_{n_\ell\in\mathcal{A}_\ell} \Delta_k\Delta_\ell\, \mathbb{P}_k^{(0)}(n_k\mid Q) \mathbb{P}_\ell^{(0)}(n_\ell\mid Q) \left[ 1-\Lambda_{k\ell}\Delta_k\Delta_\ell+\cdots \right] \nonumber\\ &= \langle \Delta_k\rangle_0 \langle \Delta_\ell\rangle_0 - \Lambda_{k\ell} \langle \Delta_k^2\rangle_0 \langle \Delta_\ell^2\rangle_0 +\cdots = - \Lambda_{k\ell}\sigma_{k,*}^2\sigma_{\ell,*}^2 +\cdots . \label{eqA:Delta95mixed95second95moment} \end{align}\tag{78}\] Combining Eqs. 75 , 77 , and 78 , the leading off-diagonal covariance is \[\mathrm{Cov}_Q(n_k,n_\ell) = - \Lambda_{k\ell}\sigma_{k,*}^2\sigma_{\ell,*}^2 +\cdots \qquad k\neq \ell . \label{eqA:offdiag95cov95Lambda}\tag{79}\] Using the definition \(\Lambda_{k\ell} = q_k^A(H^{-1})_{AB}q_\ell^B\) we obtain \[\mathrm{Cov}_Q(n_k,n_\ell) = - \sigma_{k,*}^2\sigma_{\ell,*}^2 q_k^A(H^{-1})_{AB}q_\ell^B +\cdots , \qquad k\neq \ell . \label{eqA:offdiag95cov95final}\tag{80}\] This is the leading conservation-induced intermode covariance. Including the intrinsic one-mode variance on the diagonal, the same result can be written in the compact conditional-Gaussian projection form \[C_{ij}^{Q} = \sigma_{i,*}^2\delta_{ij} - \sigma_{i,*}^2\sigma_{j,*}^2 q_i^A(H^{-1})_{AB}q_j^B . \label{eqA:projection95components}\tag{81}\] For \(i\neq j\), Eq. 81 reduces to Eq. 80 . For \(i=j\), it gives the conditional-Gaussian projection of the local variance.

This form makes exact conservation manifest. Contracting with a conserved-charge vector gives \[\begin{align} \sum_i q_i^A C_{ij}^{Q} &= \sigma_{j,*}^2 q_j^A - \sigma_{j,*}^2 \sum_i \sigma_{i,*}^2 q_i^A q_i^B (H^{-1})_{BC}q_j^C \nonumber\\ &= \sigma_{j,*}^2 q_j^A - \sigma_{j,*}^2 H^{AB}(H^{-1})_{BC}q_j^C = 0 . \label{eqA:projection95null95left} \end{align}\tag{82}\] We used the definition 46 : \(H^{AB} =\sum_i q_i^A q_i^B \sigma^2_{i,*}\) to go to the second equality. Similarly \(\sum_j C_{ij}^{Q}q_j^A=0\). Thus, the projected covariance has no component along any exactly conserved direction.

3.3 Linear observables and conservation-orthogonal projections↩︎

The covariance projection can be applied directly to measured observables that are linear in the mode occupations. Recall that \(n_i\) is the occupation number of mode \(i\). We define two such observables by \[X=\sum_i f_i n_i, \qquad Y=\sum_i g_i n_i . \label{eqA:linear95observables}\tag{83}\] Here \(f_i\) and \(g_i\) are analysis coefficients assigned to mode \(i\); they are not the microscopic counting weights \(W_i(n_i)\). For example, in heavy-ion collisions, \(f_i=1\) gives a multiplicity-type observable in the selected set of modes, while \(f_i=p_{T,i}\) gives a transverse-momentum-weighted observable.

Using Eq. 81 , the projected covariance of \(X\) and \(Y\) is \[\begin{align} \mathrm{Cov}^{\rm cg}_Q(X,Y) &= \sum_{ij} f_i g_j C_{ij}^{Q} = \sum_i f_i g_i \sigma_{i,*}^2 - U_X^A(H^{-1})_{AB}U_Y^B , \label{eqA:linear95cov95projection} \end{align}\tag{84}\] where \[U_X^A=\sum_i f_i\sigma_{i,*}^2 q_i^A, \qquad U_Y^A=\sum_i g_i\sigma_{i,*}^2 q_i^A . \label{eqA:U95defs}\tag{85}\] The vectors \(U_X^A\) and \(U_Y^A\) measure how strongly the observables overlap with the exactly conserved quantities. If \(X\) and \(Y\) are built from disjoint sets of modes, the diagonal term in Eq. 84 is absent. In that case, the long-range conservation-induced covariance is \[\mathrm{Cov}^{\rm cons}_Q(X,Y) = - U_X^A(H^{-1})_{AB}U_Y^B . \label{eqA:linear95cons95cov}\tag{86}\]

This form gives a direct prescription for constructing observables with no leading overlap with the conserved quantities. We seek modified coefficients \(f_i^\perp\) satisfying \[\sum_i \sigma_{i,*}^2 f_i^\perp q_i^C =0 \qquad \text{for all } C . \label{eqA:orthogonality95condition}\tag{87}\] Starting from a given set of analysis coefficients \(f_i\), subtract a linear combination of conserved-charge vectors, \[f_i^\perp = f_i-q_i^A \lambda_A . \label{eqA:orthogonal95ansatz}\tag{88}\] Imposing Eq. 87 gives \[\begin{align} 0 &= \sum_i \sigma_{i,*}^2 f_i q_i^C - \lambda_A \sum_i \sigma_{i,*}^2 q_i^A q_i^C = U_X^C-H^{CA}\lambda_A . \end{align}\] Hence \[\lambda_A=(H^{-1})_{AB}U_X^B ,\] and therefore \[f_i^\perp = f_i-q_i^A(H^{-1})_{AB}U_X^B . \label{eqA:orthogonal95coeff}\tag{89}\] Thus the projected observable is \[X^\perp=\sum_i f_i^\perp n_i \,.\] Using Eqs. 85 , its conserved-charge overlap vanishes, \[\begin{align} U_{X^\perp}^C = \sum_i \sigma_{i,*}^2 f_i^\perp q_i^C &= \sum_i \sigma_{i,*}^2 f_i q_i^C - \sum_i \sigma_{i,*}^2 q_i^Cq_i^A (H^{-1})_{AB}U_X^B \nonumber\\ &= U_X^C - H^{CA}(H^{-1})_{AB}U_X^B = 0 . \label{eqA:orthogonality95check} \end{align}\tag{90}\] Therefore, for any linear observable \(Y\), \[\begin{align} \mathrm{Cov}^{\rm cg}_Q(X^\perp,Y) &= \sum_i f_i^\perp g_i\sigma_{i,*}^2 - U_{X^\perp}^A(H^{-1})_{AB}U_Y^B = \sum_i f_i^\perp g_i\sigma_{i,*}^2 . \label{eqA:projected95covariance} \end{align}\tag{91}\] Therefore, the finite-rank conservation contribution vanishes once one observable is projected. If the two observables are built from disjoint sets of modes, the diagonal term is absent, and hence \[\mathrm{Cov}^{\rm cons}_Q(X^\perp,Y)=0 \,.\] Thus the projected observable has no leading conservation-induced long-range covariance with any disjoint linear observable. Local statistical correlations and dynamical correlations, if present, are not removed by this construction. This provides a practical way to design weighted observables that are orthogonal to exact conservation laws before comparing with dynamical correlations.

References↩︎