Particlization with local event-by-event conservation laws

TL;DR

This paper introduces a novel sampling method for the transition from hydrodynamics to particle transport that ensures local conservation laws are preserved event-by-event, enabling more accurate studies of fluctuations and small systems.

Contribution

The paper presents a new Metropolis-based particlization method that maintains local conservation laws during the transition from hydrodynamics to particle transport.

Findings

Preserves local conservation laws event-by-event.

Enables accurate fluctuation and correlation studies.

Applicable to small systems and stochastic hydrodynamics.

Abstract

We present a sampling method for the transition from relativistic hydrodynamics to particle transport, commonly referred to as particlization, which preserves the local event-by-event conservation of energy, momentum, baryon number, strangeness, and electric charge. The proposed method is essential for studying fluctuations and correlations by means of stochastic hydrodynamics. It is also useful for studying small systems. The method is based on Metropolis sampling applied to particles within distinct patches of the hypersurface, where hydrodynamics and kinetic evolutions are matched.

Equations

$$\left(\begin{array}{c}{P}_{tot}^{\mu }\\ B_{tot}\\ S_{tot}\\ Q_{tot}\end{array}\right)=\sum _{\begin{array}{c}\mathrm{cells}\\ i\end{array}}\int \left(\begin{array}{c}{p}_i^{\mu }\\ B_i\\ S_i\\ Q_i\end{array}\right)\frac{p^{\nu }d{\sigma }_{\nu }}{p^0}{f}_i(p^{\alpha }{u}_{\alpha },T,{\mu }_i)\frac{g_i{d}^{3}p}{(2\pi \hslash {)}^3}$$

$$P(N,\{N_s{\}}^{\mathrm{species}},\{x_i{\}}_{i=1}^N,\{p_i{\}}_{i=1}^N)=\mathcal{N}$$

$$\left(\prod _s\frac{1}{N_s!}\right)\prod _{i=1}^N\frac{g_i}{(2\pi \hslash {)}^3}\frac{d^3{p}_i}{p_i^0}{p}_i^{\mu }d{\sigma }_{\mu }{f}_i(p_i^{\nu }{u}_{\nu },T,{\mu }_i)\times$$

$${\delta }^{(4)}(\sum _i{p}^{\mu }-P_{tot}^{\mu }){\delta }_{{\sum }_i{B}_i}^{B_{tot}}{\delta }_{{\sum }_i{S}_i}^{S_{tot}}{\delta }_{{\sum }_i{Q}_i}^{Q_{tot}}$$

$$P^{t+1}(\xi )-P^t(\xi )=\sum _{{\xi }^{\prime }}w({\xi }^{\prime }\to \xi )P^t({\xi }^{\prime })-$$

$$w(\xi \to {\xi }^{\prime })P^t(\xi ).$$

$$\frac{P({\xi }^{\prime })}{P(\xi )}=\frac{w(\xi \to {\xi }^{\prime })}{w({\xi }^{\prime }\to \xi )}=\frac{T(\xi \to {\xi }^{\prime })A(\xi \to {\xi }^{\prime })}{T({\xi }^{\prime }\to \xi )A(\xi \to {\xi }^{\prime })}.$$

$$a\equiv A(\xi \to {\xi }^{\prime })=\text{min}\left(1,\frac{P({\xi }^{\prime })T({\xi }^{\prime }\to \xi )}{P(\xi )T(\xi \to {\xi }^{\prime })}\right).$$

$$d{R}_n(\sqrt{s},m_1,m_2,\dots ,m_n)=\frac{(2\pi {)}^4}{(2\pi {)}^{3n}}$$

$$\frac{d^3{p}_1}{2E_1}\frac{d^3{p}_2}{2E_2}\dots \frac{d^3{p}_n}{2E_n}{\delta }^{(4)}(P_{tot}^{\mu }-\sum P_i^{\mu }),$$

$$T(2\to 3)=\frac{1}{2}\frac{G_2^{ch}}{G_2}\frac{1}{N_3^{ch}}\frac{d{R}_3^{ch}}{R_3^{ch}}\frac{1}{N_{cells}^3}$$

$$T(3\to 2)=\frac{1}{2}\frac{G_3^{ch}}{G_3}\frac{1}{N_2^{ch}}\frac{d{R}_2^{ch}}{R_2^{ch}}\frac{1}{N_{cells}^2},$$

$$a_{n\to m}=\frac{N_m^{ch}{R}_m}{N_n^{ch}{R}_n}\frac{N!}{(N+m-n)!}\frac{m!}{n!}\frac{k_m^{id}!}{k_n^{id}!}\times$$

$${\left(\frac{2N_{cells}}{{\hslash }^3}\right)}^{m-n}\frac{\prod _{i=1}^m{g}_i{f}_i({\mu }_i-p_i^{\alpha }{u}_{\alpha },T)p_i^{\mu }d{\sigma }_{\mu }}{\prod _{j=1}^n{g}_j{f}_j({\mu }_j-p_j^{\alpha }{u}_{\alpha },T)p_j^{\mu }d{\sigma }_{\mu }}$$

Full text

One of the key goals of modern heavy ion collision experiments is to search for a phase transition between a hadron gas and a quark-gluon plasma (QGP), and to locate the corresponding critical point. The vicinity of the critical point is characterized by enhanced event-by-event fluctuations LL_T5_par146 ; Stephanov:2008qz . Therefore, considerable attention is devoted to correlation and fluctuation observables, such as proton, net-proton, net-charge, and kaon cumulants Luo:2015ewa ; Adamczyk:2017wsl , fluctuations of various particle ratios Abdelwahab:2014yha , transverse momentum correlations Adam:2019rsf , and charge balance functions Adamczyk:2015yga . Since a heavy ion reaction is a dynamical process, it is essential to study these observables within a dynamical framework. A very successful dynamical treatment is a hybrid approach Hirano:2012kj which combines the relativistic hydrodynamic evolution of the high (energy) density QGP phase with the kinetic transport for a more dilute hadronic phase. This approach successfully reproduces bulk observables such as particle spectra and flow (see e.g. Ryu:2017qzn ). Switching from the continuous relativistic hydrodynamics to the discrete particle transport, often referred to as “particlization”, is usually performed on a hypersurface characterized by constant energy density, temperature, or Knudsen number Ahmad:2016ods . In the existing relativistic models the switching only occurs in one direction: from hydrodynamics to particles, but not vice versa. This is in contrast to non-relativistic hybrid approaches, where dynamical domain decomposition methods are routinely applied and the switching is performed in both directions (see e.g. Donev:2009 ). To study correlations and fluctuations within the hybrid approach, hydrodynamics has to be either extended by stochastic terms directly Kapusta:2011gt ; Murase:2013tma ; Kumar:2013twa ; Nahrgang:2018afz or coupled to a non-equilibrium field with a stochastic noise Nahrgang:2011mg ; Herold:2014zoa . In both cases, it is essential that the particlization preserves the fluctuations generated by such models. This is a non-trivial task, which so far has not been done in the context of relativistic hydrodynamics. In the non-relativistic case, this problem is addressed in several ways Donev:2009 . One of them is to exactly match the fluxes at the interface, which in the relativistic case corresponds to local event-by-event conservation laws, or in other words, microcanonical sampling. The standard Cooper-Frye sampling used in relativistic models (for a detailed description see e.g. Huovinen:2012is ), on the other hand, is grand-canonical. It combines the Cooper-Frye formula for the momentum distribution in a hypersurface cell Cooper:1974mv with Poissonian sampling of the multiplicity distributions. Together they result in total energy, momentum, and charges fluctuating around correct means, as it is illustrated in Fig. 5 of Schwarz:2017bdg . Additionally, in the standard procedure particles in different cells are sampled independently, and thus are uncorrelated, although their velocities may still be correlated via a common flow velocity profile. The scaled variances of multiplicitites are $\omega\equiv\frac{\langle N^{2}\rangle-\langle N\rangle^{2}}{\langle N\rangle}=1$ by construction. In contrast, in the microcanonical case $\omega<1$ and particles should be correlated due to the conservation laws. Moreover, event-by-event conservation laws are important not only for correlations and fluctuations. For example, in small systems they can also affect mean values. Therefore, constructing and realizing a microcanonical sampling algorithm is the purpose of this work.

Attempts to introduce (micro)canonical particlization have been undertaken previously, but they rely on intuition-based ad hoc modifications to the standard sampling algorithm such as: introducing local charge conservation by sampling particle-antiparticle pairs Bozek:2012en , trying to satisfy conservation laws one by one in a “mode sampling” algorithm Huovinen:2012is , and rejecting particles that increase the deviation from the desired conserved charges in the SPREW algorithm Schwarz:2017bdg . The multiplicity distribution sampled by these algorithms is not known precisely and does not correspond to a canonical or microcanonical ensemble. The SER algorithm Schwarz:2017bdg ; Oliinychenko:2016vkg samples the correct canonical distribution, but extending it to the microcanonical ensemble is impossible. Note that energy and momentum conservation, which distinguish the microcanonical from the canonical ensemble, influence not only $p_{T}$ fluctuations, but also the fluctuations of multiplicities. This is why we propose a method to conserve all charges, energy, and momentum simultaneously.

First of all, we define regions over which conservation laws should be applied; we call these regions “patches”. There are two requirements for the patch size $b$. First, it should be comparable to the hydrodynamics scale discussed in the context of fluctuating hydrodynamics An:2019osr , therefore it should be much larger then the mean free path in a weakly-coupled system or a thermal length in a strongly coupled system: $b\gg 1/T$, where $T$ is temperature. Second, one patch should contain many particles, therefore $b^{3}n\gg 1$, where $n$ is particle density. A patch should not necessarily be much smaller than the system size. If a system is small, then conservation laws should be applied over the whole system and by definition the whole system is one patch. However, if a patch is comparable to the system size, then the applicability of hydrodynamics may be questionable. In the usual applications of hydrodynamics it is reasonable to use every computational grid cell as a patch. Indeed, even in simulations of micro- and nanofluids the number of particles per computational cell is of the order of 100 Donev:2009 . However, in typical simulations of relativistic ion collisions the average number of particles per computational cell is of order $10^{-3}-10^{-1}$, given a typical cell size between (0.2 fm)3 and (0.5 fm)3. While the strategy of sampling “fractional” particles is possible Steinheimer:2017dpb , it leads to complications in the treatment of the fluctuations in the subsequent transport evolution. Our strategy is to split the hypersurface into independent patches of similar size, and to require conservation laws in every patch. In practice this implies that a patch consists of roughly $50-1000$ computational cells. A typical hypersurface consists of the order of $10^{6}$ cells, resulting in about $10^{3}-10^{4}$ patches. The exact number of patches, and therefore the “localness” of the conservation laws, can be treated as a parameter. For a given patch, the quantities to be conserved are

[TABLE]

where the sum runs over all hadron species $i$ including resonances, and over all cells of the patch. Here $d\sigma_{\mu}$ is a normal four-vector of the cell (see definition in Huovinen:2012is ), $u^{\mu}$ is the collective velocity of the cell, $T$ is temperature of the cell, and $\mu_{B,S,Q}(x)$ are the chemical potentials of the cell, responsible for the conservation of baryon number $B$, strangeness $S$, and electric charge $Q$. The chemical potential of the species $i$ is defined as $\mu_{i}=\mu_{B}B_{i}+\mu_{S}S_{i}+\mu_{Q}Q_{i}$, while $g_{i}$ is the degeneracy of the species. It may seem surprising to consider a local temperature and chemical potentials in a microcanonical sampling. However, there is no contradiction here. Conservation laws are imposed only over the whole patch. Variations in energy density, quantum number densities, and collective velocities from cell to cell within a patch are allowed and characterized by local values of $T$ and $\mu$. Preserving these local variations is important to ensure a faithful description of higher order azimuthal anisotropies Gardim:2011xv , which otherwise would be smeared. The probability $P$ of a given particle configuration in a patch is a product of the usual Cooper-Frye formulas and global delta-functions which guarantee conservation laws over the patch:

[TABLE]

Note that here the number of particles of each hadron species $N_{s}$ is not fixed, and neither is the total number of particles $N=\sum_{s}N_{s}$. Instead, both are distributed according to Eq. (Microcanonical particlization with local conservation laws). The quantities $d\sigma_{\mu}$, $u^{\mu}$, $T$, and $\mu_{B,S,Q}$ depend on the spatial position of a particle $x_{i}$. The charges $B_{tot}$, $S_{tot}$ and $Q_{tot}$ are computed using Eq. (1). In practice these charges are real numbers, not integers. To address this problem, we suggest to either round $B_{tot}$, $S_{tot}$, and $Q_{tot}$ to nearest integers or distribute the non-integer parts according to a multinomial distribution, which is guaranteed not to obfuscate total charges on the hypersurface 111A detailed discussion of how to generate the patches, what their optimal size are and how to address the rounding problem will reported on in a separate detailed publication. .

Direct sampling of the $N$-particle probability distribution expressed by Eq. (Microcanonical particlization with local conservation laws) is difficult due to the unknown normalization factor $\mathcal{N}$ and the $\delta$-functions. To sample it we apply a Metropolis algorithm, a Markov chain Monte Carlo method, which in our case is closely related to solving the Boltzmann equation with the stochastic rate method Seifert:2017oyb . The state of our Markov chain $\xi$ depends on multiplicities, coordinates and momenta of all particles: $\xi=\xi(N,\{N_{s}\}^{\mathrm{species}},\{x_{i}\}_{i=1}^{N},\{p_{i}\}_{i=1}^{N})$. The initial state is an arbitrary set of particles that satisfy the required conservation laws (Eq. 1). Charge conservation for initial state is fulfilled by an ad hoc heuristic algorithm picking lightest particles of necessary charges, while the energy-momentum conservation is achieved by rescaling momenta as in Schwarz:2017bdg . Given a state $\xi$ we propose a state $\xi^{\prime}$ with probability $T(\xi\to\xi^{\prime})$ and then decide, if this state should be accepted, with probability $A(\xi\to\xi^{\prime})$. Therefore, the probability to obtain a state $\xi^{\prime}$ from $\xi$ is $w(\xi\to\xi^{\prime})=T(\xi\to\xi^{\prime})A(\xi\to\xi^{\prime})$. The master equation, connecting the probability to obtain the state $\xi$ at steps $t$ and $t+1$ is

[TABLE]

After many steps the probability $P^{t\to\infty}(\xi)$ should converge to $P(\xi)$ given by Eq. (Microcanonical particlization with local conservation laws). A sufficient condition for this is known as the detailed balance condition:

[TABLE]

This condition is satisfied if

[TABLE]

There is some freedom to select the proposal matrix $T(\xi\to\xi^{\prime})$. We choose it such that it conserves energy, momentum, and quantum numbers. Consequently, our Markov chain never leaves the desired subspace where conservation laws are fulfilled. Our proposal matrix may be viewed as $2\to{3}$ and $3\to{2}$ stochastic “collisions” Seifert:2017oyb on the hypersurface. However, we note, that there is no real time involved and “collisions” are not related to any physical process. They are simply a mathematical method to sample the distribution of Eq. (Microcanonical particlization with local conservation laws). The proposal procedure is the following:

With 50% probability choose a $2\to 3$ or $3\to 2$ transition. 2. 2.

Select the “incoming” particles by uniformly picking one of all possible pairs or triples. 3. 3.

Select the outgoing channel democratically with probability $1/N^{ch}$, where $N^{ch}$ is the number of possible channels, satisfying both quantum number and energy-momentum conservation. 4. 4.

For the selected channel sample the “collision” kinematics uniformly from the available phase space with probability $\frac{dR_{n}}{R_{n}}$, $n=2$ or $3$. 5. 5.

Choose a cell for each of the outgoing particles uniformly from all cells in the patch. Note that this choice matters for the acceptance probability, because the corresponding temperatures, chemical potentials, velocities $u^{\mu}$, and normal 4-vectors $d\sigma_{\mu}$ in the Eq. (Microcanonical particlization with local conservation laws) will be taken at the cells, where the outgoing particles are thrown.

Here $R_{n}$ is a phase-space integral for outgoing particles defined as the integral over $dR_{n}$:

[TABLE]

where $\sqrt{s}=(P_{tot}^{\mu}P^{tot}_{\mu})^{1/2}$. The integration of $dR_{2}$ and $dR_{3}$ is possible analytically Bauberger:1994nk ; Seifert:2017oyb . Our proposal procedure generates the following probabilities for $2\to 3$ and $3\to 2$ proposals:

[TABLE]

where $G_{2}=\frac{N(N-1)}{2!}$ and $G_{3}=\frac{N(N-1)(N-2)}{3!}$ denote total numbers of incoming pairs and triplets of any species, while $G_{2}^{ch}$ and $G_{3}^{ch}$ are the numbers of ways to select a given incoming particle species. Consequently, $\frac{G_{2}^{ch}}{G_{2}}$ and $\frac{G_{3}^{ch}}{G_{3}}$ represent the probabilities to obtain pairs and triplets of a given incoming species. The number of possible triplets and pairs of outgoing species with appropriate quantum numbers are denoted by $N^{ch}_{3}$ and $N^{ch}_{2}$. Inserting the proposal probabilities, Eqs. (7) and (8), as well as the desired probability distribution, Eq. (Microcanonical particlization with local conservation laws), into the expression for the acceptance probability, Eq. (5), we arrive, after some algebra, at

[TABLE]

where we made use of the relation $\prod\frac{d^{3}p_{i}}{(2\pi\hbar)^{3}p^{0}_{i}}\delta^{(4)}(P_{tot}^{\mu}-\sum P_{i}^{\mu})=2^{n}\frac{dR_{n}}{(2\pi)^{4}}$. Here $n=2,3$ and $m=3,2$ are the numbers of incoming and outgoing particles, and $N$ is the total number of particles before proposing the Markov chain step. The product in the numerator is taken over the outgoing particles and the one in the denominator is taken over the incoming particles. The quantities $d\sigma$, $u$, $T$, $\mu$ should be evaluated in the cell where the particles are proposed to be, or coming from. The total number of particles in the entire patch is given by $N$, and $k^{id}_{m}$ and $k^{id}_{n}$ are the numbers of outgoing and incoming identical species in the reaction. Note that the sampling accounts for the variations in temperature and chemical potential within the patch. Also, and equally important, the distribution function $f$ may contain viscous corrections. To summarize, the algorithm consists of multiple Markov chain steps, where the step is proposed according to Eqs. (7) and (8) and accepted with probability given by Eq. (Microcanonical particlization with local conservation laws).

We have tested the above sampling algorithm both in a simplified and a somewhat more realistic setup. First, we consider a patch consisting of one cell with $u^{\mu}=(1,0,0,0)$, $d\sigma_{\mu}=(V,0,0,0)$, and $f(p)=\exp(-p^{0}/T)$. In this case Eq. (Microcanonical particlization with local conservation laws) represents a well-known microcanonical distribution, which has been sampled before Werner:1995mx ; Becattini:2004rq . Additionally, analytical expectations of scaled variances of hadron multiplicities in the large volume limit are available in this case Hauer:2007ju . Our algorithm reproduces the analytical results for both means and scaled variances well. The resulting momentum distributions are very close to Boltzmann, as expected. This is a non-trivial result, because multiplicity and momentum distributions are a not a direct input to the sampling; the only input is volume $V$, total energy, momentum and conserved charges.

Next, we demonstrate our sampling for a more realistic scenario, where we consider a patch consisting of three cells with non-trivial values for $u^{\mu}$, $d\sigma_{\mu}$, and $T$, which also vary from cell to cell 222We have also successfully tested the above algorithm including patch selection on realistic hyper-surfaces, which will be reported on in detail in a subsequent publication. . Conservation laws are imposed over the entire patch, while the local energy density and charge densities vary from cell to cell. As shown in Fig. 1(a), we obtain the expected means in each cell, agreeing with the grand-canonical standard Cooper-Frye sampling. This is a demonstration that we correctly reproduce the temperature variation from cell to cell. The scaled variances, $\omega$, of multiplicities in the patch, shown in Fig. 1(b), are already drastically different from the standard grand-canonical Cooper-Frye sampling, where $\omega=1$ by construction. In our case the variances agree with the microcanonical analytical expectation from Hauer:2007ju . Finally, in Fig. 1(c) we demonstrate the non-trivial correlations emerging from conservation laws. Unlike for variances, to our knowledge, there is no analytical calculation of correlations in a microcanonical ensemble, although in principle such calculations are possible using techniques developed in Hauer:2007ju . Beside testing the sampler implementation, Fig.1 1(b,c) also shows the expected effect of conservation laws on fluctuations and correlations in heavy ion collisions.

In summary, we have proposed, implemented and tested a particlization method which takes into account local event-by-event conservation laws in a systematic fashion. Localness is achieved by splitting the hypersurface into patches and enforcing conservation laws in every patch. Event-by-event conservation of total energy, momentum, baryon number, strangeness, and electric charge over the patch is guaranteed by the algorithm. At the same time local cell-by-cell variations of energy and charge densities within a patch are preserved, ensuring that observables sensitive to these variations, such as higher order azimuthal asymmetries Gardim:2011xv , are not smeared out. The proposed method is essential for studies of correlations and fluctuations, especially in combination with stochastic hydrodynamics, since it will not obfuscate its correlations and fluctuations. It may also be applied to exploring small systems, where the impact of event-by-event conservation laws is large, as well as charge dependent correlations relevant for the the chiral magnetic effect Voloshin:2004vk . The code used for sampling is publicly available Oliinychenko:code . We have checked that its execution time for realistic hypersurfaces is not impractical 333The overall runtime to generate $10^{5}$ samples for our test scenario constituted around 2 minutes. For a realistic hypersurface expected from a heavy ion collision of $\sqrt{s}\simeq 10\,\rm GeV$. where we have about 150 patches of $\sim 10\,\rm GeV$, the present implementation takes about 40 minutes to generate 10000 samples As a next step, we will apply it to search for observable effects of critical fluctuations in heavy ion collisions. Through accounting for local conservation laws we might be able to detect critical fluctuations in coordinate space via correlations in momentum space which were previously not visible when using standard particlization.