Calculating n-Point Charge Correlations in Evolving Systems

TL;DR

This paper introduces a diagrammatic formalism to compute the time evolution of non-local charge correlations in dynamic systems, aiding the interpretation of charge fluctuation measurements in heavy-ion collisions.

Contribution

It presents a novel diagrammatic approach to calculate evolving non-local charge correlations assuming local features are equilibrated.

Findings

Provides a method to compute $n$-point charge correlations at different positions.

Accounts for correlations involving up, down, and strange charges.

Aims to interpret charge fluctuation data from the Relativistic Heavy-Ion Collider.

Abstract

In dynamic systems, charge susceptibilities and local charge correlations change with time. These changes are accompanied by non-local correlations which spread diffusively with time and are constrained by local charge conservation. Assuming the local features of the correlation, which for a gas would be the correlation of charges within the same particle, are equilibrated, a diagrammatic formalism is presented for calculating the evolution of the associated non-local correlations. These provide correlations of $n$ density operators at different positions for arbitrary $n$. The techniques were developed with an eye towards relativistic heavy-ion collisions, and can account for correlations indexed by up, down and strange charges. Understanding the evolution of such correlations is crucial if one is to interpret measurements of charge fluctuations from the Relativistic Heavy-Ion Collider.

Equations

$${\mathcal{C}}_{ab}^{(\mathrm{tot})}({\mathbf{r}}_1,{\mathbf{r}}_2,t)$$

$${\mathbf{r}}_{ij}$$

$${\chi }_{ab}^{(2)}$$

$$D_1{C}_{ab}^{(\mathrm{tot})}({\mathbf{r}}_1,t_1,{\mathbf{r}}_2,t_2)$$

$$D_i$$

$$C_{a;b}^{(1;1)}({\mathbf{r}}_1+{\mathbf{v}}_{1}dt,t_1+dt,{\mathbf{r}}_2+{\mathbf{v}}_{2}dt,t_2+dt)$$

$$D_t$$

$$D_t{C}_{ab}^{(\mathrm{tot})}({\mathbf{r}}_1,{\mathbf{r}}_2,t)$$

$${\mathbf{r}}_{12}$$

$$D_t{C}_{a;b}^{(1;1)}({\mathbf{r}}_1,{\mathbf{r}}_2,t)$$

$$S_{ab}^{(2)}({\mathbf{r}}_{12},t)$$

$$\int d^3{r}_1{d}^3{r}_2{C}_{a;b}^{(1;1)}({\mathbf{r}}_1,{\mathbf{r}}_2,t)$$

$$\left[{\partial }_t+\mathbf{v}\cdot \nabla \right]\delta V$$

$$D_t{\chi }_{ab}^{(2)}(\mathbf{r},t)$$

$${\mathcal{C}}_{abc}^{(\mathrm{tot})}({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3,t)$$

$$\delta Q_d$$

$$\delta N_s$$

$$\delta Q_d$$

$$\delta {\mu }_a$$

$$\delta N_s$$

$${\int }_{\in \delta V}{d}^3{r}_{12}{d}^3({\mathbf{r}}_1-{\mathbf{r}}_2)\delta {\rho }_a({\mathbf{r}}_1)\delta {\rho }_b({\mathbf{r}}_2)\delta ({\mathbf{r}}_1-{\mathbf{r}}_2)$$

$$C_{ab;c}^{(2;1)}({\mathbf{r}}_{12},{\mathbf{r}}_3,t)$$

$$L_{ab,e}^{(2)}({\mathbf{r}}_{12},t)$$

$$C_{a\cdots c;\cdots }^{(m;\cdots )}(\mathbf{r},\cdots ,t)$$

$$L_{a\cdots c,d}^{(m)}(\mathbf{r},t)$$

$$⟨\delta Q_a\delta Q_b\delta Q_c⟩$$

$$\frac{d}{dt}⟨\delta Q_a\delta Q_b\delta Q_c{⟩}_V$$

$$D_t{\mathcal{C}}_{abc}^{(\mathrm{tot})}({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3,t)$$

$$D_t{C}_{ab;c}^{(2;1)}({\mathbf{r}}_{12},{\mathbf{r}}_3,t)$$

$$[d_t{L}_{ab,d}^{(2)}(\mathbf{r},t)]⟨\delta {\rho }_d(\mathbf{r},t)X⟩$$

$$D_t{C}_{a;b;c}^{(1;1;1)}({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3,t)$$

$$S_{abc}^{(3)}(\mathbf{r},t)$$

$$S_{ab;c}^{(2;1)}(\mathbf{r},{\mathbf{r}}^{\prime },t)$$

$$⟨\delta {\rho }_b({\mathbf{r}}_2,t_2)⟩$$

$$\int d^3{r}_2{G}_{ab}({\mathbf{r}}_1,t_1,{\mathbf{r}}_2,t_2)$$

Full text

I Introduction

Correlations and fluctuations of conserved charges play a central role in heavy-ion collisions. Charge fluctuations represent defining property of any bulk system, particularly in characterizing phase transitions. As temperatures rise above $\approx 160$ MeV, matter undergoes a transition from a hadronic gas to a strongly interacting plasma of quarks and gluons, the quark-gluon plasma (QGP). For neutral matter, equal numbers of particles and anti-particles, lattice gauge theory has shown that the transition is a smooth crossover occurring within a rather narrow window of temperatures, $150\lesssim T\lesssim 180$ MeV. For temperatures below 150 MeV charge susceptibilities from lattice calculations are consistent with expectations for a weakly interacting gas of hadrons, and for temperatures above 180 MeV, they become consistent with a weakly interacting gas of light up, down and strange quarks Koch:2008ia ; Bzdak:2012an ; Borsanyi:2011sw ; Bazavov:2012jq ; Bazavov:2014xya . This consistency is rather surprising given the experimental evidence that the system behaves like a nearly ideal liquid, with mean free paths on the order of the thermal wavelength Romatschke:2007mq ; Song:2010mg ; Pratt:2015zsa ; Auvinen:2017fjw . Fluctuations between baryon charge and strangeness Koch:2005vg , and baryon fluctuations of third or fourth order, e.g. $\langle\delta Q^{3}\rangle$ or $\langle\delta Q^{4}\rangle$, are especially illuminating Athanasiou:2010kw . They suggest that charges fluctuate in units of one third baryon number rather than in units of baryon number once temperatures rise above the aforementioned window Bazavov:2014xya . At finite baryon density, where lattice calculations struggle due to a sign problem, the possibility exists for a first-order phase transition, which would culminate at a critical point Stephanov:2008qz ; Stephanov:2011pb ; Vovchenko:2015pya . The baryon density of the critical point might be several times normal nuclear density and the critical temperature would likely be moderately less than the temperature window quoted above for the smooth transition at zero baryon density. Baryon fluctuations should represent identifying properties of such a phase transition, particularly near the critical point. Fluctuations of baryon number in high energy heavy ion collisions have been analyzed as a function of beam energy at the Relativistic Heavy Ion Collider (RHIC) Aggarwal:2010wy ; Adamczyk:2013dal ; Adamczyk:2014fia ; Adamczyk:2017wsl ; McDonald:2013aoa . At the highest RHIC energies, the incoming beams are insufficiently stopped to contribute large numbers of baryons to the mid-rapidity region. Combined with rampant particle production at such energies, net baryon densities are much lower than entropy densities and experiments are able to investigate the properties of matter with nearly zero baryon chemical potential, $\mu_{B}\approx 0$. Baryon densities increase for lower beam energies, which provides the opportunity to study the properties of matter as a function of baryon density at high temperature.

Unfortunately, interpreting charge fluctuations is greatly complicated by the dynamic nature of the collision. A system’s baryon density and temperature traverse a swath through the density-temperature plane. Thus, any measurement reflects on the bulk properties over a range of density and temperature. Further, the short lifetime of these environments restricts charge fluctuations within any volume from attaining equilibrium values because charge is locally conserved and requires significant time to diffuse across the volume. Any quantitative model of charge fluctuations must therefore describe the evolution of charge correlations, $\langle\delta\rho(\bm{r}_{1},t)\cdots\delta\rho(\bm{r}_{n},t)\rangle$, in order to understand how the equilibrium properties one would study in a static system would become manifest in the finite-size and finite-volume systems created in heavy-ion collisions, where measurement of the conserved charges are confined to the final state.

Models have addressed the challenges outlined in the previous paragraph, but mainly for two-point correlations Stephanov:2017ghc ; Stephanov:2017wlw ; Pratt:2017oyf ; Pratt:2018ebf . As reviewed in the next section, two-point correlations can be split into two pieces. The first piece is the short-range contribution. This is when the two density operators in the correlation, $\langle\delta\rho(\bm{r}_{1},t)\delta\rho(\bm{r}_{2},t)\rangle$, refer to charges within the same particle. Throughout this paper, the word “particle” can be extended to any short-range feature that might locally equilibrate. For example, in a hadron gas this would be the case where the two density operators referred to the same hadron, and for a QGP gas, this would be when the two density operators address the same quark. Assuming knowledge of the local part of the correlation part, e.g. assuming a chemically equilibrated gas, the remainder of the correlation function is constrained by the fact that the net correlation function must integrate to zero due to local charge conservation. If one assumes that charges move diffusively, the evolution of the non-local part of the correlation function can be modeled by the diffusion equation with a source term given by the rate of change of the local part. Such an approach was superimposed onto a hydrodynamic description of a heavy-ion collision in Pratt:2017oyf and then extended to include a hybrid hydrodynamic model interfaced to a hadronic simulation that simulated the break-up stage of a heavy-ion collision Pratt:2018ebf ; Pratt:2019pnd . These approaches roughly reproduced several experimental measures of charge correlations from the STAR collaboration at RHIC, which should translate into reproducing experimental measures of charge fluctuations to order $\langle\delta Q^{2}\rangle$, because fluctuations are determined by integrating over the correlations.

The aim of this work is to provide a theoretical foundation to extend the treatment of two-point correlations to three-point, four-point and $n$-point correlations. This is more difficult than the case for two-point correlations. In that case the correlation was divided into a local part, where the two density operators referred to the same particle, and a non-local part. For three-point correlations, $\langle\delta\rho_{a}(\bm{r}_{1},t)\delta\rho_{b}(\bm{r}_{2},t)\delta\rho(\bm{r}_{3},t)\rangle$, one must consider three cases: where all three density operators refer to the same particle, where two of the three refer to the same particle, and where all three refer to different particles. Even if one makes an assumption about the local part for the three-point function, i.e. it reflects chemical equilibrium, one must understand how to split the remaining correlation over the other two possibilities. Four-point and higher correlations offer even more possible splittings. Section III shows how three-point correlations can be modeled. The question of how to spread correlations amongst the various splittings is answered by assuming that a differential charge results in differential changes in the various species according to equilibrium. This then allows correlations of order $\delta\rho^{3}$ due to a two-point function to be determined by the two-point function contribution to correlations of order $\delta\rho^{2}$. The algebra in Sec. III is rather lengthy, and a similar exposition for $n$-point correlations with $n>3$ would be much more so. Fortunately, the expressions can be represented diagrammatically. A diagrammatic description, which is extendible to higher $n$, is presented in Sec. IV.

In Pratt:2017oyf the diffusion equation for the two-point correlation function was addressed by noting the equivalence with a random walk. In Sec. VI the benefits of a random-walk algorithm vs. a mesh-based description of the correlation function is discussed. In a quark gas, up, down and strange quarks represent good quasi-particles and the three-by-three diffusivity tensor is diagonal. In a hadron gas, hadrons carry multiple quarks and the diffusivity tensor, just like the susceptibility, is no longer diagonal in the $u,d,s$ basis. A strategy for applying a random walk algorithm in a situation where the diffusivity tensor is not diagonal is also provided in Sec. VI.

The final section, Sec. VII, presents a discussion of the applicability of the relations from Sec.s III and IV. The role of assuming chemical equilibrium is emphasized. Finally, strategies are presented for handling both local and non-local contributions to the susceptibility. Critical phenomena involves correlation on longer length scales and seems well suited for hydrodynamic treatments Stephanov:2017ghc ; Stephanov:2017wlw ; Stephanov:2012ki ; Pratt:2017lce ; Paech:2003fe ; Paech:2005cx ; Nahrgang:2011mg ; Kapusta:2014dja ; Kapusta:2012sd ; Young:2014pka ; Ling:2013ksb . Phase separation dynamics might also be addressed with such an approach Steinheimer:2013xxa ; Steinheimer:2013gla ; Steinheimer:2012gc ; Randrup:2010ax ; Heiselberg:AnnPhys ; Heiselberg:1988oha ; Napolitani:2014ima ; Chomaz:2003dz ; Borderie:2001jg ; Colonna:2002ti ; Guarnera:1996svb . However, hydrodynamics, noisy or not, is a clumsy means by which to model the correlation of a particle with itself because hydrodynamics is based on gradients, which implies that correlations have a length scale greater than the inter-particle separation. Section VII describes the possibility of combining a hydrodynamics approach to account for the non-local contribution to the susceptibility and the formalism presented here to account for the local part. This study considers only evolving the correlations in coordinate space, whereas measurements are restricted to the asymptotic momenta. Techniques for translating correlations to momentum space has been described and implemented in Pratt:2011bc ; Pratt:2015jsa , and it would be straight-forward to extend these methods to project $n-$point correlations into momentum space. Implementations of the formalism presented here will be deferred for another study.

II Two-point functions

Before launching into a formalism for three- and four-point functions, that for two-point functions is reviewed here. This has been appeared in Pratt:2017oyf and applied to a hydrodynamic evolution of a heavy-ion collision in Pratt:2017oyf ; Pratt:2018ebf ; Pratt:2019pnd .

First, the definitions,

[TABLE]

The subscript $a$ denotes the various charges, perhaps the $u,d,s$ charges on quarks. Here, $\delta\rho_{a}=\rho_{a}-\langle\delta\rho_{a}\rangle$ so that $\langle\delta\rho_{a}\rangle=0$. The subscripts on the correlations indicate whether the charges are on the same particle, as $\chi^{(2)}(\bm{r},t)$ describes the contributions where both charges come from the same point, or from the same particle, whereas $C^{(1;1)}$ encapsulates the contribution to the correlation when the charges are on different particles. The semicolon in $C^{(1;1)}_{a;b}$ emphasizes that the two charges $a$ and $b$ are not on the same particle. For the considerations of this paper, it will be assumed that the local part is understood, i.e. it could be the equilibrated susceptibility if the particles are well defined and are in chemical equilibrium. In the gaseous limit, $\chi^{(2)}_{ab}$ is the correlation of the charges within a particle,

[TABLE]

where $n_{s}$ is the number density of species $s$, and $q_{s,a}$ is the charge of type $a$ on a particle of species $s$. For the example of a hadron gas, the contribution from $\pi^{+}$ mesons to $\chi^{(2)}_{ud}$ is $-n_{\pi+}$, where $n_{\pi+}$ is the density of $\pi^{+}$ mesons. The negative sign ensues because the $\pi^{+}$ meson has an up quark and an anti-down quark. Even for individual quarks, one finds a contribution to $\chi^{(2)}$ from the correlations of quarks with themselves. For a gas of quarks, $\chi^{(2)}_{ss}=n_{s}+n_{\bar{s}}$, the density of strange plus that of anti-strange quarks. In this paper, the quasi-particles that carry charge will be referred to as particles. Particles could refer to point charges, hadrons, atoms, molecules, or could even include a local polarization cloud. The non-local part, $C^{(1;1)}$, will diffuse and spread over large relative coordinates. Providing the theoretical structure for calculating the evolution of $C^{(1;1)}$, for the case of two-particle correlations, and $C^{(1;1;1)}$ or $C^{(1;1;1;1)}$ for three- or four-particle correlations, is the principal goal of this paper. The local correlation, whose strength is $\chi^{(n)}$, will be assumed to be given, by assuming local chemical equilibrium.

The evolution of the correlation is guided by the equation,

[TABLE]

Here, $\bm{v}$ is the local velocity of the fluid, and $\bm{j}_{a}$ is the current measured in the fluid frame, i.e. it neglects the part of the current from $\delta\rho_{a}\bm{v}$. The definition of $D_{i}$ differs from the usual definition of a co-moving derivative because of the presence of the term $\nabla\cdot\bm{v}$. That term accounts for the current $\bm{j}$ being measured relative to the local frame of the fluid. If one were to include the term $\delta\rho\bm{v}$ to the current, this additional contribution to $D_{i}$ would not be necessary. With this definition, $\bm{j}$ can be considered as the diffusive contribution to the current, i.e. it ignores the part from simple fluid movement. Because the r.h.s of Eq. (3) is a divergence, this represents local charge conservation. A corresponding equation is also true for $D_{2}$. To propagate the equal-time correlation forward,

[TABLE]

or in terms of $C^{(1;1)}$,

[TABLE]

The last term, with $S^{(2)}_{ab}(\bm{r}_{12},t)$, behaves like a source function for $C^{(1;1)}$,

[TABLE]

The source term contributes to the strength of the correlation at $\bm{r}_{1}-\bm{r}_{2}=0$. For a small fluid element of volume $\delta V$ that expands with the fluid, the source contributes when the product $\chi\delta V$ changes with time. This is a consequence of the definition of $D_{t}$ including the $\nabla\cdot\bm{v}$ term, because the usual comoving derivative, $\partial_{t}+\bm{v}\cdot\nabla$, acting on $\delta V$ gives

[TABLE]

For ideal hydrodynamics, the entropy within $\delta V$, which equals $s\delta V$, would remain constant. In that case,

[TABLE]

and one can see that the source term is principally a function of whether the ratio $\chi/s$ rises or falls as one moves with the fluid. If entropy is not conserved, the source term differs somewhat. In the treatments of Pratt:2011bc ; Pratt:2019pnd ; Pratt:2011bc the hydrodynamic evolution was viscous, and the source term was calculated with the full expression given in Eq. (II). Nonetheless, the approximate form in Eq. (8) is insightful, as plotting $\chi_{ab}/s$ as a function of temperature describes at what points in the trajectory the source term becomes significant. Further, this ratio can be calculated in lattice gauge theory Pratt:2015jsa .

Although the expressions involve two powers of the density, the evolution of $C^{(1;1)}$ is described by a linear equation including a source term. For each differential contribution to the source function, $-D_{t}\chi^{(2)}_{ab}d^{3}r~{}dt$, one can solve for its contribution of $C^{1,1}$. Finally, one can sum each contribution by integrating over the source function. If the evolution is diffusive, $\bm{j}_{a}(\bm{r})=-\mathcal{D}_{ab}(\bm{r})\nabla\delta\rho_{b}(\bm{r})$, where $\mathcal{D}$ is the diffusivity tensor. If the diffusivity tensor is diagonal each of these contributions can be represented by two sample charges $a$ and $b$ undergoing a random walk with the parameters of the random walk set by the diffusivity tensor. The positions of the two charges can then be used to construct the correlation function in coordinate space. It is then straight-forward to design a Monte Carlo procedure to generate pairs for each contribution from the source function at some point $\bm{r}_{s}$ and time $t_{s}$. Because only the charges originating from the same source point are correlated with one another, there is no combinatoric noise to overcome. This approach was applied in Pratt:2017oyf ; Pratt:2018ebf ; Pratt:2019pnd . A method for handling non-diagonal diffusivity tensors is provided in Sec. VI.

III Three-Point Correlators

Without loss of generality, the three-point correlator can be written as

[TABLE]

Here, $\bm{r}_{123}\equiv(\bm{r}_{1}+\bm{r}_{2}+\bm{r}_{3})/3$. The correlation $C^{(1;1;1)}$ describes correlations when all three positions are different, i.e. the density operators refer to different particles, and $C^{(2;1)}_{ab;c}$ describes the correlations when two positions are the same, i.e. the two charges $a$ and $b$ are on the same particle and $c$ is on a separate particle. The correlation when all three points are the same, or all three charges are on the same particle, is described by $\chi^{(3)}_{abc}$. Just as with $\chi^{(2)}_{ab}$, this will be identified as the equilibrium susceptibility here. If the particles are molecules, assigning $\chi^{(n)}$ as the equilibrium susceptibility represents an assumption of chemical equilibrium.

The correlator, $C^{(2;1)}_{ab;c}(\bm{r}_{12},\bm{r}_{3},t)$, describes the correlation between a charge of type $c$ at $\bm{r}_{3}$ and a particle at position $\bm{r}_{12}$ carrying a product of charges $Q_{a}Q_{b}$. Here, we show that it is directly determined by $C^{(1;1)}_{d;c}(\bm{r}_{12},\bm{r}_{3},t)$ and the susceptibilities. To demonstrate this relation, one can consider a particle of species $s$. The charge $\delta Q_{d}$ due to the increased probability of having a particle $\delta N_{s}$, is

[TABLE]

where $q_{sd}$ is the charge of type $d$ on the particle of type $s$. If the particle probability is equilibrated in response to the small charge,

[TABLE]

where $\delta\mu_{a}$ is the chemical potential inspired by the small charges, divided by the temperature. One can insert Eq. (11) into (10),

[TABLE]

where $[\chi^{(2)}]^{-1}_{ab}$ is the inverse two-point susceptibility matrix. Inserting this into Eq. (11),

[TABLE]

This expresses how many extra particles of type $s$, $\delta N_{s}$, one would generate in a volume when a small charge, $\delta Q_{a}$, is added to the volume.

One can now calculate the additional product of charges $\delta(Q_{a}Q_{b})$ due to $\delta Q_{c}$. To that end, one can consider a small volume $\delta V$ restricting the position $\bm{r}_{12}$. The delta function, should not be of zero extent, but should have a range large enough to fit in a quasi-particle.

[TABLE]

Here, $\delta(Q_{a}Q_{b})$ refers to the charges inside the volume $\delta V$. It is indeed this product of charges in a single particle at $\bm{r}_{12}$ that is described by $C^{(2;1)}_{ab;c}(\bm{r}_{12},\bm{r}_{3},t)$. Thus,

[TABLE]

By assuming that $\chi^{(3)}$ is consistent with chemical equilibrium, all mention of the individual particles and their charges has disappeared, and $C^{(2;1)}_{ab;c}(\bm{r}_{12},\bm{r}_{3},t)$ is determined by the correlation and the susceptibilities evaluated at $\bm{r}_{12}$. For future reference, one can readily show that for any product of $m$ charge densities $\rho_{a}\cdots\rho_{c}$,

[TABLE]

Our principal goal is to determine the evolution of $C^{(1;1;1)}_{a;b;c}(\bm{r}_{1},\bm{r}_{2},\bm{r}_{3},t)$. Assuming that the two-point correlation $C^{(1;1)}_{a;b}$ was already determined using the methods of Sec. II, all terms from the r.h.s. of Eq. (III) involving two-point functions can be calculated from Eq. (15). Here, we first solve for $D_{t}=D_{1}+D_{2}+D_{3}$ of the l.h.s. of the equation, i.e. $D_{t}$ acting on the total correlation. Then applying $D_{t}$ to the r.h.s. will provide an expression for $D_{t}C^{(1;1;1)}$.

Before applying $D_{t}$ to the l.h.s. of Eq. (III), one can surround the points $\bm{r}_{1}$, $\bm{r}_{2}$ and $\bm{r}_{3}$ with surfaces and consider the net correlation of the product of charges within the enclosing volumes $V_{1}$ $V_{2}$ and $V_{3}$,

[TABLE]

The rate of change of $C_{abc}$ is determined by the rate at which charge flows out of the small encircling volumes,

[TABLE]

Here, the prime on the averages $\langle\cdots\rangle^{\prime}$ restricts the integrals to not include charges in the same particle. The last several terms used Eq. (15) to relate how the product of charges carried by a single particle is determined by the single-charge charge density. If the volumes are moving and expanding with the fluid, and if the currents $\bm{j}_{a}$ are defined relative to the fluid, one can use the divergence theorem to rewrite Eq. (III) in differential form with $d/dt$ replaced by $D_{t}$,

[TABLE]

Putting all these terms together gives the result for applying $D_{t}$ to the l.h.s. of Eq. (III),

[TABLE]

Next, one applies $D_{t}$ to the r.h.s. of Eq. (III). First, a sample term is considered where two of the charges are carried by the same particle,

[TABLE]

Here, the definition of $d_{t}$ includes the statistical average to its right, $\langle\delta\rho(\bm{r},t)X\rangle$,

[TABLE]

Here, $X$ could refer to any operator away from the position $\bm{r}$. The quantitiy $\langle{\bm{j}}(\bm{r},t)_{d}X\rangle$ is reexpressed as a ratio over $\langle\delta\rho_{d}(\bm{r},t)X\rangle$ multiplied the same quantity. This is motivated so that one can see that $d_{t}$ is effectively the co-moving derivative, but co-moving in the frame of the current, which is not necessarily the same as the frame of the fluid. Thus, if $\delta\rho$ is represented by Monte Carlo sampling, the derivative $d_{t}$ would refer to the rate of change according to an observer moving with the sampling particles.

Comparing Eq. (III) to Eq. (III) one can see that many of the terms cancel. The resulting equation expresses the evolution of $C^{(1;1;1)}$,

[TABLE]

The first three terms in Eq. (III) describe how the correlations evolve when all three coordinates differ. If the current is diffusive, $\bm{j}_{a}=-\mathcal{D}_{ab}\nabla\delta\rho_{b}$, the correlations spread with time. The remaining terms represent source terms for $C^{(1;1;1)}$. In the absence of the source terms $C^{(1;1;1)}$ would integrate to a constant. The last four terms describe the sourcing of $C^{(1;1;1)}$ for instances when at least two of the coordinates are equal. The term proportional to $D_{t}\chi^{(3)}$ was expected because $\chi^{(3)}$ describes the correlation when all three charges are on the same particle. The three terms proportional to $L^{(2)}S^{(2)}$ describe how some of the correlation of $\delta Q_{a}\delta Q_{b}\delta Q_{b}$ is absorbed by the change of the two-point function, i.e. two of the charges are on one particle and the third on a second particle. The sources $S^{(2;1)}$ describe how the three point function can be seeded with two points on the same particle, and one on a separated particle. The factor $d_{t}L^{(2)}$ describes how the correlation of two charges carried by one particle split onto two particles if $L^{(2)}=\chi^{(3)}[\chi^{(2)}]^{-1}$ would change with time. In the next section, a graphical scheme is presented, which provides some visual delineation of the various terms above, while providing the means to write down the corresponding terms for four-point or $n-$point correlations.

IV Graphical representations and higher-order correlations

In the previous section, equations of motion were found for three-point correlation functions in Eq. (III). Combined with the expressions for two-point functions in Eq. (II), and using Eq. (15), one can find all correlations of order $\delta\rho^{3}$. The evolution of the two and three-point functions, described in Eq.s (II) and (III), can also be expressed graphically. The elements of the graphs are lines connected by vertices, with the vertices having either zero or one incoming lines and $n$ outgoing lines. The lines will connect space time points $\bm{r}_{1},t_{1}$ and $\bm{r}_{2},t_{2}$ and are Green’s functions describing how charge a charge $\delta Q_{a}$, placed at $\bm{r}_{1},t_{1}$ would affect the density, $\delta\rho_{b}$, at a point $\bm{r}_{2},t_{2}$, where $t_{2}>t_{1}$.

[TABLE]

The Green function is normalized,

[TABLE]

and obeys the boundary condition at $t_{1}=t_{2}$,

[TABLE]

For a diffusive equation, $\bm{j}_{a}=-\mathcal{D}_{ab}\nabla\delta\rho_{b}$, the Green’s function can be calculated by solving the differential equation,

[TABLE]

For any realistic dynamic system, it is unlikely $G$ can be found analytically. The choices are either to solve the differential equation numerically on a three-dimensional mesh, or to sample the diffusive spread as a random walk. The sources of the Green’s function for two- and three-point functions are listed in Eq.s (II) and (III) respectively. Sources for the Green’s functions can be represented diagramatically, with vertices representing sources at points in space time, and lines between the vertices representing Green’s functions. These graphical elements are illustrated in Fig. 1 and the vertices are defined below,

[TABLE]

Figure 2 shows the diagrams for calculating two-point, three-point and four-point functions. Many of the diagrams are topologically identical and are related by permuting the final-state labels. In those cases the similar diagrams are noted by the number of permutations for that topology.

As an example, the contribution to the three-point diagram from the second three-point diagram in Fig. 2 is the integral

[TABLE]

Each vertex in the diagram is assigned a space-time point, in this case $y_{1}$ and $y_{2}$. Integrations are performed over those coordinates. Each internal line is assigned two charge indices which then determine the charge indices for the vertices. All diagrams begin with a vertex $V^{(0\rightarrow n)}$, and end with open Green’s functions denoted by the desired measurement.

V Relation to Charge Fluctuations

Within some large volume $V$, charge fluctuations are defined

[TABLE]

Each charge $Q_{a}$ can expressed as an integral over the charge density $\delta\rho_{a}$. For the order $Q^{n}$ fluctuation, one obtains contributions from the two-point, three-point, up to $n-$point functions. The contribution from the $n-$point function is simply the integral over all the external coordinates in the diagrams from Fig. 2. The contributions from the $(n-1)-$point functions with final-state charge indices $a$ and $b$ can be found by attaching an operator $L^{(2)}_{ab,a^{\prime}}(x)$ to any external Green’s function $G_{d^{\prime}a^{\prime}}(y,x)$ where $x$ is a final-state coordinate and $a^{\prime}$ denotes the measured charge. Thus, each 3-point diagram from Fig. 2 contributes to $F^{(4)}$. The contributions to $F^{(4)}$ from two-point functions come from either attaching $L^{(2)}$ to both of the external lines, or by attaching $L^{(3)}$ to either external line. Finally, $F^{(4)}$ has a contribution from all four charges being on the same particle, which would be represented by $\chi^{(4)}$.

Experimentally, the contributions to $F^{(4)}$ from four-point functions come from summing over all combinations of four final-state particles, never using the same particle twice in the same term. The contribution to $F^{(4)}$ from three-point functions would be found by summing over all sets of three final-state particles then requiring one particle to provide two powers of the charge. The contributions from two-point functions describes the case where the sum extends over all pairs, with each particle contributing an order $Q^{2}$ contribution or for one particle to provide an order $Q$ and the second providing an order $Q^{3}$ contribution. Finally, summing over the particles individually, one would add the contribution of $Q_{a}Q_{b}Q_{c}Q_{d}$ for that particle. Aside from the contribution to $F^{(n)}$ from the $n-$point function, all other contributions are determined by correlations of fewer coordinates, and thus do not represent additional information beyond what would have been gathered by $(n-1)-$point functions.

VI Algorithms

In principle, one could solve the differential equations for the correlation functions. The differential equation would involve solving for all points on a grid with three spatial dimensions and one time dimension. If the three-dimensional space-time grid was represented by $N\times N\times N$ grid points calculated for $N_{t}$ values of the time, an $n$-point correlation function would involve of the order $N^{3n}N_{t}$ grid points. This would be likely be prohibitively expensive.

Another possibility for calculating $n$-point correlations would be to solve $n$ separate one-point diffusion equations on $n$ meshes. For each source point, one would increment the correlations on each of the meshes. For each source point, $S^{(n)}_{ab\cdots c}d^{4}x$, one could increment the charges on the corresponding mesh points by amounts $\delta Q_{a},\cdots\delta Q_{c}$, such that the product of the charges reproduced $S^{(n)}_{ab\cdots c}d^{4}x$, but so that the incremented charge on any of the individual mesh had a random sign. At the final time, one would construct the $n-$point correlation function by using the charges from each of the $n$ meshes. Unless the contributions came from the same source point, they would, on average, cancel. Unfortunately, the cancellation would require repeating the procedure many times to combat combinatoric noise. This would be especially true for $n>2$. For $n>2$ it would be more efficient to evolve the contribution from a single source-point, then construct the correlation. One would perform a Monte Carlo sampling over the many source points using $|S^{(n)}_{ab\cdots c}d^{4}x|$ as the probability to choose the sampling points, then use $S^{(n)}_{ab\cdots c}/|S^{(n)}_{ab\cdots c}|$ as a weight to increment the correlations functions.

A second approach, built on the assumption that the dynamics is diffusive, is to represent the correlations with clusters of sample particles undergoing random walks. For $n-$point correlations, the clusters would involve $n$ charges. Sample charges move with some velocity $\bm{v}$ and then have their directions reoriented randomly. The probability that a particle is thus scattered during a time interval $dt$ is $dt/\tau$, with $\tau=6\mathcal{D}/v^{2}$. In the limit $v\rightarrow\infty$ the random walk approaches the diffusion equation. By setting $v$ to the speed of light, it is causal and approaches the diffusion equation after several scatterings. If the diffusivity, $\mathcal{D}_{ab}$, is not diagonal, a more sophisticated representation would need to be invoked, and is described below. Each cluster of particles would be created via a Monte Carlo procedure weighted by the source function, and would evolve as a random walk to mimic the diffusion equation. Each particle carries a unit charge, and each group of particles would carry a weight, which could be either positive or negative, describing the contribution from the original vertex after accounting for the Monte Carlo weight. When calculating the correlations, only those particles within the same cluster need to be combined, which results in low combinatoric noise. Such an approach was applied in Pratt:2019pnd ; Pratt:2018ebf ; Pratt:2017oyf . In those instances, only two-point functions were considered. Calculations for three-and four-point functions would involve accounting for a larger number of diagrams. One advantage of this approach is that vertices of the form $V^{(1\rightarrow n)}$ would be rather straight forward to calculate. These vertices behave as $d_{t}L^{(n)}_{a\cdots c,b}(x)$. The derivative $d_{t}$, defined in Eq. (III), is the time-derivative that co-moves with the current, or in this case is co-moving with the sampling particles. Thus, thus the $(1\rightarrow n)$ vertices involve calculating how $L_{a\cdots c,d}$ changes according to an observer moving with the sampling particles. This simplifies sampling the secondary vertices with Monte Carlo.

For either approach, the $n-$point functions must be addressed in order. The two-point evolution can be used to calculate the three-point evolution, and the two- and three-point correlations serve as a basis for the four-point function. For the random-walk representation, one stores the correlated clusters of correlated particles. When evolving a pair of particles to represent the two-point function, one could bifurcate one of the particles carrying charge $d$ into two with charges $a$ and $b$ during a time interval $dt$ with probability $V^{(1\rightarrow 2)}_{d,ab}(\bm{r},t)dt$. One would continue to simulate the non-split trajectory for calculation of the two-point function, and would add the bifurcated trajectory into a list of samplings for the three-point function. Such a three-point trajectory, generated from an initially two-point trajectory, would represent the second diagram for 3-point correlations in Fig. 2. This would be added to the purely three-point trajectory described by the first three-point diagram in Fig. 2. Similarly, one can calculate four-point functions. Given the lack of combinatoric noise, such a calculation would be tenable and require only modest computational resources. However, the simple random walk approach needs to be altered if the diffusivity is not diagonal, as described below.

VI.1 Non-diagonal diffusivity matrix

It is straight-forward to model the diffusive evolution to the $n-$point contribution of the $\delta\rho^{n}$ correlation, $C^{(1;1;\cdots 1)}_{a;b;\cdots c}(\bm{r}_{1},\bm{r}_{2},\cdots\bm{r}_{n},t)$, for the case where the diffusivity matrix is diagonal, if one is given the source function $S^{(n)}_{ab\cdots c}(\bm{r}_{1},\bm{r}_{2}\cdots\bm{r}_{n})$. Because diffusion represents a random walk, one simply creates a set of unit charges, $a\cdots c$ at the points $\bm{r}_{1}\cdots\bm{r}_{n}$, and assigns a weight to the group. The weight is given by $S^{(n)}_{ab\cdots c}d^{4}x/P_{MC}$. The probability $P_{MC}$ accounts for the fact that in each four-volume element $d^{4}x$, one may choose whether or not to create the sampling charges. If the Monte Carlo probability, $P_{MC}$, is chosen as $|S^{(n)}_{ab\cdots c}d^{4}x|$, then the weights are $\pm 1$, depending on whether the sources are positive or negative. Each charge is then propagated as a random walk. The sample charges move with velocity $v$ relative to the medium, and with random directions. The charges then reorient randomly according to a lifetime, $\tau_{\rm coll}=6\mathcal{D}/v^{2}$. I.e., in each time step $dt$ the particle reorients with probability $dt/\tau_{\rm coll}$.

There are significant advantages to using a random walk representation of the diffusion equation. First, such implementations tend to be simpler to implement than the solving the differential equation on a four-dimensional space-time grid. Second, one can make the evolution causal by setting $v$ to the speed of light. In the limit of $v\rightarrow\infty$ the random walk exactly reproduces the diffusion equation, but by lowering $v$ the method prunes the acausal tail of the correlation function. The difference between causal and acausal treatments matters only for short diffusion times. For long times, the random walk approaches the solution to the diffusion equation as long as the evolutions involve many reorientations for each test charge. For the calculations in Pratt:2017oyf ; Pratt:2018ebf ; Pratt:2019pnd , the number of such reorientations was on the order of a half dozen. Finally, the random walk makes it easy to label the contribution to the correlation from the same source point. Contributions from different source points should cancel, so by only incrementing contributions from the same source point the combinatoric noise is greatly reduced. This becomes increasingly important as one considers correlations of increasing order.

Treating diffusion as a random walk is more complicated once the diffusivity tensor becomes non-diagonal,

[TABLE]

The method will be based on considering sample charges in a basis where $\mathcal{D}$ is diagonal. Here, the eigenvectors of $\mathcal{D}$ are labeled $u^{(i)}$. The source functions and susceptibilities can be expressed in this new basis, and the sample charges are labeled by eigenvectors. If the basis were constant throughout the evolution, the algorithm would then be unchanged from what was described above. If one applies Eq. (13) to translate a sample charge into particles of a specific species, the unit sample charge $\delta Q_{a}$ is $u^{(i)}_{a}$, where $i$ refers to the specific eigenvector representing the sample charge.

The non-diagonal elements are thus rather straight-forward to accommodate if the eigenvectors of the diffusivity tensor do not change as the sample charge traverses the medium. However, when the eigenvectors transform they must be reformulated in terms of the new eigenvectors. Let’s assume the original normalized eigenvectors were $\hat{a}$, $\hat{b}$ and $\hat{c}$. The new eigenvectors will be $\hat{a}^{\prime}$, $\hat{b}^{\prime}$ and $\hat{c}^{\prime}$. Also, one can assume the charge is originally in the state $u=\hat{a}$. Using completeness,

[TABLE]

For the Monte Carlo treatment one can probabilistically choose which new eigenvector to use along with an adjustment of the weight so that on average the charge is still in $\hat{a}$. Here, $w$ designates the original weight assigned to the group of sampling charges, and $w^{\prime}$ will be the new weight after the charge has been re-designated in the new basis. One can generate a random number $r$ such that $0<r<1$. Using $r$, the following algorithm should maintain the continuity of the charge,

[TABLE]

If one averages over values of $r$, the result averages to $u\rightarrow(\hat{a}\cdot\hat{a}^{\prime})\hat{a}^{\prime}+(\hat{a}\cdot\hat{b}^{\prime})\hat{b}^{\prime}+(\hat{a}\cdot\hat{c}^{\prime})\hat{c}^{\prime}$, which indeed equals $u=\hat{a}$. One need only check whether to reassign the basis with time steps sufficiently small so that the change in the diffusive movement is small during that time step.

VII Applicability

The approach here was inspired by understanding how three- and four-point correlations measured in heavy-ion collisions could be modeled. In particular, the goal was to understand the role of local charge conservation. During a heavy-ion collision, one changes phases from a quark-gluon plasma to a hadronic gas. At zero baryon chemical potential, calculations of the susceptibilities from lattice gauge theory suggest that for temperatures below 150 MeV the system is reasonably represented as a hadronic gas, whereas for temperatures above 200 MeV quarks are reasonable quasi-particles. For intermediate temperatures, the transition appears smooth. During a central collision of heavy ions at LHC energies or at RHIC energies, the system traverses a range of temperatures from well above 200 MeV to approximately 100 MeV, and undergoes a radical change in chemistry during that time. The changing number of up, down and strange charges, and the combination of such charges into hadrons induces a rich evolution of charge correlations. If such correlations are short range, and if the expansion is not too fast, it seems a reasonable approximation to assume chemical equilibrium, at least until the temperatures fall below 150 MeV, at which point chemical rates fall below the expansion rate. If the chemistry is equilibrated, one would expect the short range, $\lesssim 1$ fm, correlations to match that of an equilibrated gas. However, beyond one Fermi correlations due to local charge conservation persist. The treatment presented here would thus seem a reasonable approximation to reality, and comparing predictions to measurements would provide a stringent test of the assumption of local chemical equilibrium. Indeed, for two-point correlations this approach has matched a range of experimental measurements.

If chemical equilibrium is not attained the approach can become invalid. For two-point correlations one can assume the local correlation is some function $\chi_{ab}\sim\delta(\bm{r}-\bm{r}^{\prime})$, and if the function $\chi_{ab}$ can be modeled the approach can still be applied as chemical equilibrium was not an essential approximation. As long as the correlation is local, compared to the size of the system, the same approach, but with a different model for the local correlation, would remain valid. However, the derivations for three- and four-body correlations did rest on the assumptions of chemical equilibrium. Even if one had a model of the local correlations, $\chi^{(3)},\chi^{(4)},\cdots$, that would not be sufficient to understand how a charge $\langle\delta Q_{a}\rangle$ on a particle would translate into knowing $\delta Q_{b}\delta Q_{c}$ on the same particle. Local chemical equilibrium was critical in deriving Eq. (15).

Aside from chemical equilibrium, the second assumption is that the correlation is sufficiently local to separate it from the balancing correlation. This should be true in most cases, but would fail for correlations associated with phase transitions. In the critical region correlations fall as power laws, effectively with infinite extent. When inside the coexistence region, bubbles and drops represent macroscopic structures that should not be described by an expansion of $\langle\delta\rho^{n}\rangle$. For high-energy heavy-ion physics, there remains the possibility that a phase transition exists at finite baryon density, and might be accessible at the lower range of beam energies at RHIC. If the correlations from bulk structure are sufficiently long range, it is possible that short-range correlations, e.g. those from charge conservation, might be superimposed onto a model where a one-body description Steinheimer:2013xxa ; Steinheimer:2013gla ; Steinheimer:2012gc ; Randrup:2010ax ; Heiselberg:AnnPhys ; Chomaz:2003dz ; Borderie:2001jg ; Colonna:2002ti ; Guarnera:1996svb ; Paech:2005cx ; Nahrgang:2011mg , including some with implementations of noise Napolitani:2014ima ; Gavin:2016hmv ; Kapusta:2012sd ; Young:2014pka ; Pratt:2017lce . The strategy would be then to first treat the bulk correlations, including the critical correlations and those related to phase separation, using some form of hydrodynamics. Highly local correlations, including their contribution to the susceptibilities, would be ignored for this first pass. The methods presented here could then be applied to account for the remainder of the correlation, i.e. those from short-range correlations and the associated balancing charge.

Even if one’s main motivation for analyzing multi-charge correlations and fluctuations is to search for evidence of phenomena related to phase transitions, it is crucial to estimate the degree to which the short-range correlations and the associated charge balance affect the result. The methods presented here provide a means to calculate that background. If one’s goal is to investigate the chemical evolution of a heavy-ion collision through $n-$point correlations in a system where there are only short-range correlations, the methods here make it possible to extend such studies to $n>2$.

Acknowledgements.

This work was supported by the Department of Energy Office of Science through grant number DE-FG02-03ER41259, and benefited from conversations within the Beam Energy Scan Theory (BEST) Topical Collaboration, also supported by the Department of Energy.