Finite density lattice QCD without extrapolation: Bulk thermodynamics with physical quark masses from the canonical ensemble

This paper presents the first lattice QCD results for bulk thermodynamics at finite baryon density with physical quark masses using a canonical ensemble approach that avoids extrapolation and the sign problem up to $\mu_B \approx 500$ MeV.

The Gist

Imagine trying to understand the behavior of a crowd at a massive concert. You want to know how the crowd moves, how loud it gets, and what happens when the music changes. In the world of physics, this "crowd" is made of quarks and gluons (the building blocks of matter), and the "concert" is the extreme environment inside a neutron star or the moments right after the Big Bang.

This paper is about a new, clever way to study this crowd without getting lost in a mathematical nightmare.

The Problem: The "Ghost" in the Machine

For decades, physicists have used a supercomputer method called Lattice QCD to simulate this crowd. Usually, they simulate the crowd in a "Grand Canonical" way. Think of this like a party where the number of guests can change freely, but the host (the computer) has to keep a running tally of the "chemical potential" (a fancy term for how much the guests want to be there).

The problem is that when the party gets crowded (high density), the math turns into complex numbers (numbers with imaginary parts). In the world of probability and computer simulations, you can't have a "negative" or "imaginary" chance of a guest showing up. It's like trying to bake a cake with a recipe that calls for "-3 eggs." The computer gets confused, and the results become unreliable. This is known as the Sign Problem.

To get around this, scientists usually tried to simulate the empty party (zero density) and then use math tricks (extrapolation) to guess what happens when the party gets crowded. But guessing is risky; it's like trying to predict a hurricane by looking at a gentle breeze.

The Solution: The "Canonical" Detour

This paper introduces a new strategy: The Canonical Ensemble.

Instead of letting the guest count fluctuate, the researchers decide to simulate the party with exactly 1 guest, then exactly 2 guests, then exactly 3, and so on.

• The Analogy: Imagine you are studying a dance floor. Instead of watching a chaotic, shifting crowd, you lock the door and say, "Okay, today we only have exactly 5 dancers." You study their moves. Then you let 6 dancers in and study them.
• Why it works: When you fix the number of guests, the "ghostly" imaginary numbers disappear. The math becomes clean and real. You can simulate these specific scenarios directly on the computer without the sign problem ruining the data.

The Magic Trick: The "Replica" Bridge

Here is the tricky part: The real world (and the Big Bang) doesn't have a fixed number of guests. The number fluctuates. So, how do we get from "Exactly 5 guests" back to "A fluctuating crowd"?

The authors developed a brilliant mathematical bridge.

1. Step 1: They simulated the party with imaginary chemical potentials (a safe, mathematical version of the problem) to get a baseline.
2. Step 2: They used a Fourier Transform (think of this as a complex prism that splits light into colors) to convert that baseline data into the "Exactly N guests" scenarios.
3. Step 3 (The Innovation): They realized that if you take two copies of your "5-guest" simulation and glue them together, you get a "10-guest" system. If you glue them together mathematically in a specific way, you can simulate a system that is twice as big, or three times as big.
   • The Analogy: Imagine you have a small, perfect model of a city with 5 people. You can't just build a bigger city easily. But if you take your model, make a copy, and mathematically "stretch" the space between them, you can simulate a city with 10 people, then 20, then 100.
   • By doing this "stretching" (which they call the replica method) and then shrinking the scale back down to infinity, they can reconstruct the behavior of the real, fluctuating crowd without ever having to guess or extrapolate.

Why This Matters

1. No More Guessing: Previous methods were like looking at a map and guessing the terrain ahead. This method is like driving the car and seeing the road directly. They get the results for high-density matter without relying on mathematical approximations that might break down.
2. Physical Masses: They did this with "real" quark masses (the actual weight of the particles in nature), not just theoretical approximations. This makes the results much more relevant to real-world physics.
3. The Phase Diagram: They successfully mapped out the "Phase Diagram" of matter up to a certain density (about 500 MeV). This is like drawing a map that shows exactly when matter turns from a solid (like a neutron star) into a hot soup (quark-gluon plasma).

The Bottom Line

The authors found a way to bypass the "Sign Problem" by changing the rules of the game. Instead of trying to simulate a chaotic, fluctuating crowd directly (which breaks the math), they simulated fixed-size groups, used a clever mathematical "stretching" trick to connect the groups, and then reconstructed the chaotic crowd from those pieces.

It's a bit like solving a jigsaw puzzle by first building perfect, small sections of the picture, and then using a special glue to snap them together into the final, massive image, rather than trying to force all the pieces into place at once.

This opens the door to understanding the densest matter in the universe with a level of precision that was previously thought impossible.

Technical Summary

1. Problem Statement

Lattice Quantum Chromodynamics (QCD) is the primary tool for studying the phase diagram of strongly interacting matter. However, simulating QCD at finite baryon density (real baryo-chemical potential, $\mu_B > 0$) is hindered by the sign problem.

• The Sign Problem: At real $\mu_B$, the fermion determinant becomes complex, rendering standard importance sampling (Monte Carlo) impossible because the probability weight is not real and positive.
• Current Limitations: Existing methods rely on indirect approaches:
  • Taylor Expansion: Expanding observables around $\mu_B = 0$. This suffers from a cancellation problem (sign problem) that grows exponentially with the order of expansion and volume. Furthermore, it requires truncation, introducing systematic errors.
  • Imaginary Chemical Potential: Simulating at imaginary $\mu_B$ and analytically continuing to real values. This is limited by the radius of convergence and potential phase transitions (e.g., Roberge-Weiss transition).
  • Reweighting: Directly reweighting from $\mu_B = 0$ to real $\mu_B$. This introduces severe overlap problems and, crucially, is ambiguous for rooted staggered fermions (the standard formulation for physical quark masses), leading to unphysical lattice artifacts (fake transitions).

The authors aim to provide direct, truncation-free results for finite density QCD with physical quark masses without relying on extrapolation in $\mu_B$, specifically addressing the rooting ambiguity of staggered fermions.

2. Methodology

The authors propose a three-step hybrid approach combining the canonical ensemble with a novel limit procedure to recover grand canonical results.

A. Canonical Ensemble Formulation

Instead of simulating at real $\mu_B$, the authors work in the canonical ensemble where the net baryon number $N_B$ is fixed to an integer.

• Input: They start with high-statistics grand canonical ensembles generated at $\mu_B = 0$ on a $16^3 \times 8$ lattice with physical quark masses ($N_f = 2+1$) using the 4HEX staggered action.
• Fourier Transform: The canonical partition function $Z_C(N_B)$ is obtained by Fourier transforming the grand canonical ratio $Z_{GCR}(\phi)$ (where $\phi$ is the imaginary chemical potential) over the interval $[0, 2\pi]$:
  $$Z_C(N_B) = \int_0^{2\pi} \frac{d\phi}{2\pi} e^{-i\phi N_B} Z_{GCR}(\phi)$$
• Rooting Solution: By performing the Fourier transform after the ensemble average (rather than on individual configurations), the rooting procedure (taking the 1/4 power of the determinant) is applied to real, positive numbers (the ensemble average of determinants at imaginary $\mu$). This circumvents the ambiguity of rooting complex numbers, a major flaw in previous reweighting attempts.

B. Algorithmic Optimization: Center Sector Restoration

A critical technical challenge is that dynamical simulations at $\mu_B=0$ break the $Z_3$ center symmetry, leading to poor sampling of non-trivial center sectors. This causes large statistical errors in the Fourier transform.

• Solution: The authors implement a sector partitioning algorithm. They classify configurations into three center sectors ($S_0, S_1, S_2$) based on which imaginary chemical potential ($0, \pm i2\pi/3$) yields the largest determinant weight.
• Restoration: They reconstruct the full expectation value by combining results from the dominant sector ($S_0$) with shifted versions of the data, effectively restoring the $Z_3$ symmetry and ensuring all sectors are sampled correctly without requiring separate simulations.

C. Recovering the Grand Canonical Limit (The $\alpha$-Replica Method)

Canonical results are volume-dependent and defined only for integer $N_B$. To obtain physical, volume-independent grand canonical observables (like pressure and chemical potential at a specific density $n_B$), the authors introduce a replica parameter $\alpha$.

• Concept: They define a generalized free energy $F_\alpha$ for a system with volume $\alpha V$ and particle number $\alpha N$.
• Asymptotic Limit: Using the relation derived from the Legendre transform, they show that the grand canonical free energy is recovered in the limit $\alpha \to \infty$:
  $$F(T, V, N) \approx \lim_{\alpha \to \infty} \frac{1}{\alpha} F(T, \alpha V, \alpha N)$$
• Implementation: They compute observables for various integer values of $\alpha$ (effectively scaling the density $n_B = N/(\alpha V)$) and perform a numerical extrapolation to $\alpha^{-1} \to 0$. This allows them to map discrete canonical data points to a continuous grand canonical equation of state.

3. Key Contributions

1. First Physical Mass Results: This is the first time canonical ensemble results have been presented for QCD with physical quark masses (pion mass $\approx 135$ MeV).
2. Rooting-Free Finite Density: The method successfully bypasses the "rooting problem" at finite density. By rooting the ensemble average of real determinants (via imaginary $\mu$) rather than complex ones, they avoid the fake transitions and lattice artifacts associated with rooted staggered fermions at real $\mu_B$.
3. No Extrapolation in $\mu_B$: Unlike Taylor expansion, this method does not rely on a truncated series expansion in $\mu_B$. It provides direct results for specific densities, limited only by the sign problem's severity at high densities, not by the order of a polynomial fit.
4. Center Symmetry Restoration: The introduction of the optimal sector partitioning algorithm significantly reduces statistical errors in the Fourier transform, making the method feasible for physical volumes.

4. Results

The authors mapped the QCD phase diagram and calculated bulk thermodynamic observables up to a baryo-chemical potential of $\mu_B \approx 500$ MeV.

• Phase Diagram: They generated contours of constant net-baryon density ($n_B$) on the $T-\mu_B$ plane. The results are consistent with other theoretical approaches (Dyson-Schwinger, Functional Renormalization Group) in the crossover region.
• Pressure and Chemical Potential:
  • They calculated the pressure ($p$) and chemical potential ($\mu_B$) as functions of density.
  • The results show excellent agreement with the Hadron Resonance Gas (HRG) model at low temperatures and densities.
  • At higher densities, deviations from HRG appear, consistent with repulsive baryon interactions.
• Comparison with Taylor Expansion:
  • The canonical results are compared to Taylor expansions (up to order $N^{10}$).
  • The canonical method yields smaller statistical errors and more stable results, particularly at higher temperatures and densities where the Taylor expansion suffers from large cancellations and truncation errors.
  • The Taylor expansion exhibits non-physical non-monotonic behavior in some regions, whereas the canonical results remain smooth.
• Susceptibility: The baryon susceptibility ($\chi_2$) was extracted via numerical differentiation of the $\mu_B(n_B)$ curve, showing the expected sigmoid shape at low $\mu_B$ and the development of a peak at higher $\mu_B$.

5. Significance and Outlook

• Direct Access to Finite Density: This work demonstrates that lattice QCD can provide direct, non-extrapolated results for finite density matter, moving beyond the role of merely providing expansion coefficients.
• Feasibility: The sign problem is manageable with current computing power up to $\mu_B \approx 500$ MeV for the studied volume ($16^3 \times 8$). This covers a significant portion of the parameter space relevant to the RHIC Beam Energy Scan program.
• Future Directions:
  • Volume Scaling: The current volume is small ($L \approx 2$ fm). Increasing the volume is necessary to access chiral observables and reduce finite-size effects, though this exponentially increases the sign problem severity.
  • Strangeness: The formalism can be extended to include strangeness neutrality ($\mu_S \neq 0$) to better match heavy-ion collision conditions.
  • Continuum Limit: While the current results are at finite lattice spacing ($N_\tau = 8$), the method is applicable to finer lattices, though eigenvalue extraction becomes computationally more expensive.

In conclusion, the paper presents a robust, technically optimized framework for studying finite density QCD with physical parameters, offering a viable alternative to Taylor expansions and reweighting, and providing the first direct lattice constraints on the QCD equation of state at $\mu_B \approx 500$ MeV.