Main

Condensed 4He undergoes a transition from a normal liquid to a superfluid phase at a critical temperature Tc 2.17 K, at its saturated vapour pressure9,10. Superfluid 4He was the first experimentally realized, and remains the most extensively studied, quantum phase of matter. Anomalous phenomena such as dissipationless flow, non-classical rotational inertia, quantized vortices and the Josephson effect have been thoroughly experimentally characterized11. Early theoretical work demonstrated the quantum mechanical origin of these phenomena12,13,14 where the Hamiltonian of liquid 4He is that of interacting spinless, non-relativistic bosons. Continuous space quantum Monte Carlo methods enable the precise computation of a wide range of its microscopic and thermodynamic properties, confirming theoretical predictions and reproducing experimental observations15. Moving beyond conventional simulations, recent algorithmic advances have opened up the possibility of measuring entanglement, the non-classical information shared between parts of a quantum state, in numerical experiments16,17. We combine these two technologies to measure the entanglement entropy in the superfluid phase of bulk 4He at zero temperature. Its ground state, | Ψ〉, in a cubic volume can be bipartitioned into a spherical subregion A and its complement as shown in Fig. 1. The standard measure of entanglement between A and is the Rényi entropy, Sα(A) ≡ log(TrρAα)/(1 − α), where ρA is the reduced density matrix of the subsystem: . The α = 1 case is most commonly known as the von Neumann entropy. The integer α ≥ 2 entropies have special physical significance, since they are related to the expectation value of an operator18. This allows for their evaluation using conventional measurement techniques, without resorting to full-state tomography. Since, from a many-body physics perspective, relevant features of entanglement (such as the area law) are quantifiable for any Rényi entropy, the natural choice, measured both in numerical simulations as well as recent experiments19, is α = 2. Here, we investigate its dependence on the radius R of the spherical subregion over a range of densities in the superfluid and find a dominant area law scaling: S2 R2.

Figure 1: Entanglement across a spherical boundary.
figure 1

A container of superfluid 4He is bipartitioned into a spherical subregion A of radius R and its complement at fixed density n ≡ 1/r03. The entanglement between A and is dominated by an area law, scaling with area of the bounding surface.

While there is no proof of the area law outside of a restrictive case in one spatial dimension20, it is the leading contribution to a scaling form that can be argued to arise from a few fundamental physical principles6,21,22,23,24. These are: S2(A) arises from correlations local to the entangling surface; and it has contributions at all length scales ranging from the microscopic scale of the interactions, r0, up to the characteristic size of the system, rf = min[R, ξ], where ξ is a correlation length. From these, a simple phenomenological scaling theory can be inferred for a spherical boundary of radius R in three dimensions (Fig. 1). For each infinitesimal region of the bounding surface dΣ, there is a local contribution to S2 from each length scale . For a given , the lowest order dimensionless quantity that can contribute to S2 is dΣ/2; when integrated over the surface this provides a contribution R2/2. To account for contributions at all length scales, defined in the renormalization group sense, we integrate using a logarithmic measure23. This measure can be checked, for example, in the case of one dimension for a critical system with entangled region of size L, where it reproduces the appropriate scaling S2 log(L) (ref. 18). In three dimensions the resulting integral is S 2 ~ r 0 r f ( R / l ) 2 d ( log l ) ~ R 2 : an area law. While higher order corrections lead to a power series in R, the symmetry of the entanglement entropy between complementary regions of pure states, , limits this expansion to even powers (in odd dimensions). This leads to the generic scaling form,

where a, b and c are dimensionless numbers. a is non-universal and depends on the microscopic details of the system, while b and c potentially encode universal information that is independent of the short-distance physics. Note that the leading-order area (as opposed to a volume) scaling and the absence of a subleading linear term are both features of the underlying physical principles in three dimensions.

We perform numerical tests of the general scaling form of equation (1) via high-performance simulations of superfluid 4He. Measuring S2 is significantly more computationally complex than for conventional estimators such as the energy and required the development of a new algorithm described in the Methods. The combined results of S2(R) for RL/2 are shown in Fig. 2 for the ground state of 4He at equilibrium density. We find that the entanglement entropy for different numbers of particles N collapses to a nearly universal curve. Before finite-size effects dominate near R L/2, and for R r0, we can fit the data to the scaling form in equation (1) with the dominant behaviour captured by a two-parameter fit (a, c) shown as a solid line. The extracted value of a is robust within 5% for fits including b ≠ 0 (details are provided in the Supplementary Information).

Figure 2: Entanglement entropy of superfluid 4He.
figure 2

The ground state of N superfluid helium-4 atoms in a cubic cell with linear dimension L and periodic boundary conditions is bipartitioned into a sphere of radius R and its complement at fixed equilibrium density n0 = N/L3 0.02186 Å−3. The Rényi entanglement entropy S2 in the limit R r0 is driven by particle fluctuations into the spherical subregion and is well described by free bosons (dotted line). For r0 R L/2 we perform the simplest fit of the numerical data to the leading-order behaviour of the scaling form in equation (1), using a two-parameter fit (a, c) with b = 0, shown as a solid line. Deviations from the universal curve for each system size occur as R approaches L/2, due to finite-size effects.

The efficacy of this fit and its confirmation of the leading-order scaling behaviour of the entanglement is investigated by computing the residuals between the simulation data and two scaling forms as shown in Fig. 3. We explore the area law predicted by equation (1), and a volume law that would be expected for an extensive entropy of thermodynamic origin. The residuals for the area law are consistent with zero, while strong deviations for the volume law exclude this as a candidate for the leading-order scaling. We observe no evidence of a subleading linear correction to the area law as predicted by equation (1). The systematic investigation of further subleading (non-constant) terms commensurate with equation (1) would require simulations of larger system sizes, providing a wider range of length scales.

Figure 3: Area or volume law?
figure 3

The residuals of the Rényi entanglement entropy computed via quantum Monte Carlo (QMC) and two types of two-parameter fit corresponding to an area law: S2fit = 4πa(R/r0)2 + c and a volume law: S2fit = (4πa/3)(R/r0)3 + c for N = 64 4He atoms at the equilibrium density n0. Residuals are shown for values of the spherical subregion radius R over which the fits were performed. The data are poorly described by the volume law and support the area law scaling predicted in equation (1).

To understand the physical origin of the area law scaling coefficient a, we define an entanglement length scale . From a fit to equation (1), e 1.3r0 5 Å at n = n0. This strongly suggests that the short-distance physics of the potential hard core and adjacent attractive minima dominate the area law scaling behaviour.

To confirm, we study the effect of the density on the entanglement by computing S2 for superfluid 4He over a range of densities near n0 corresponding to positive (n > n0) and negative (n < n0) pressures. Performing a two-parameter fit to the area law scaling for each density, we plot an ‘entanglement equation of state’ in Fig. 4. a is an increasing function of density, and thus a monotonic function of pressure (see inset). We find that e depends both on the nature of short-distance interaction as well as the interparticle separation. We can contrast this behaviour to the non-interacting Bose gas, where S2(R) is a pure function of the aspect ratio R/L.

Figure 4: Entanglement equation of state.
figure 4

The area law coefficient a versus density n in the ground state of 4He in the superfluid regime, as computed by two-parameter fits to quantum Monte Carlo data for N = 64 (symbols); the line is a guide for the eye. Inset: the pressure P of the ground state of superfluid 4He as a function of density, from ref. 30. This suggests that a is a monotonic function of pressure. The dashed vertical lines correspond to the equilibrium density n0 0.02186 Å−3.

In conclusion, we have demonstrated that the prototypical quantum fluid, superfluid 4He, displays area law scaling of its entanglement entropy. Using large-scale, exact microscopic simulations, we have extracted the numerical coefficient of the area law term and find that it is a monotonically increasing function of density. This confirms that fluctuations and interactions local to the entangling boundary drive the physics of the area law. These fluctuations also play an important role in constraining the subleading scaling of the entanglement entropy, which contains new universal physics. For example, it is predicted that logarithmic corrections should arise due to the existence of a spontaneously broken continuous symmetry in the thermodynamic limit, contributing a universal coefficient due to the presence of a low-energy ‘tower of states’ spectrum and a Goldstone boson25,26,27. For superfluid 4He with a spherical entangling surface, this b coefficient will combine with another universal number arising from the vacuum theory governing the bosonic fluctuations. This latter quantity encodes one of the two central charges that characterize a three-dimensional conformal field theory28, believed to be the fundamental constant that quantifies how entropymonotonically decreases under renormalization group flow29. A curved bounding surface such as a sphere without defects is only possible in the spatial continuum. It is thus possible that fundamental physical quantities that arise for smooth geometries are inaccessible in simple lattice models and can be probed only in Galilean-invariant quantum liquids.

Methods

The Hamiltonian of bulk liquid 4He is that of spinless, non-relativistic bosons interacting with a two-body interatomic potential:

where 2 is the reduced Planck’s constant and m is the mass. Accurate microscopic interatomic potentials V for 4He have been developed31 that include a repulsive hard core of radius σ 2.6 Å, an attractive power-law tail, and a minimum of depth ε 11 K at radius rm 3.0 Å. At the equilibrium number density n0 0.02186 Å−3, the mean interparticle separation r0 3.6 Å is slightly larger than rm. The use of V in conjunction with path integral quantum Monte Carlo methods has allowed for precise calculations of a wide range of ground state and finite temperature properties of liquid helium15 for N of order 104 atoms32.

In recent work we have extended these Monte Carlo methods to compute Rényi entropies in systems of itinerant particles in the spatial continuum16,17,33,34. Due to the increased computational difficulty of computing Rényi entropies, the data presented here are limited to N ≤ 64 4He atoms. However, systems of this size have been demonstrated to be sufficiently large to display fundamental macroscopic features of superfluid 4He (ref. 15).

While our previously published algorithm was focused on one spatial dimension, for this work we have developed a new variant that allows for the computation of Rényi entropies in the ground state of systems of interacting bosons in three-dimensional continuous space. We use a path integral method15,35 that gives access to group state properties via imaginary-time projection on a trial state |ΨT〉:

We label the classical configuration space of N bosons by R, which is a vector of three-dimensional particle coordinates. Considering Hamiltonians of the form equation (2), we use a standard approximation to the imaginary-time propagator

which is accurate to fourth order in the short time τ (refs 36,37). Because equation (4) is non-negative for bosonic systems, we can Monte Carlo sample discrete imaginary-time world-line configurations, where we use P discrete time steps with 2β = .

For the second Rényi entropy, we define a replicated Hilbert space of two non-interacting copies of the system, . We may compute S2 under a bipartition of the system into A and its complement from the expectation value of a ‘swap operator’ that swaps the configuration of A between the two replicas:

where such that RA () is a vector of the coordinates of the particles in A (). The estimator for S2 is then simply related to the expectation value of the swap operator16:

To compute equation (5) with path-integral ground-state Monte Carlo, we use equation (3), and consider imaginary-time paths of length 2β, capped by ΨT on either end of the path. We Monte Carlo sample an extended configuration space of imaginary-time world lines that includes configurations where world lines that pass through A at imaginary time β may swap between replicas. That is, world lines that pass through the spatial subregion at time β are always propagated in imaginary time along the same replica, but particles that pass through the A subsystem at time β may be connected via ρτ to a world line in or RA at time β + τ. Such swapped world-line configurations have weight of the form

where ρτA () are the reduced propagators for the A () subsystem17,34 and Rj is the configuration of N particles at imaginary time . By including updates that allow for the world-line connectivity to interchange between swapped (s) and unswapped (u) configurations, we may measure S2 from equation (5) from the ratio of the swapped and unswapped generalized partition functions:

We find this estimator to be more efficient than previous variants for systems above one spatial dimension. To further improve its performance we used a ‘ratio method’ to build up A from smaller increments33,38. The systematic errors due to finite τ and β can be made arbitrarily small by increasing P at a computational cost that is polynomial in P and N with details shown in the Supplementary Information. All results shown were computed using β = 0.48 K−1, τ = 0.005 K−1, and a constant trial wavefunction.

Data availability.

All quantum Monte Carlo data that were used to generate the plots within this paper and other findings of this study are available from the corresponding author on reasonable request.

Additional Information

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.