Quantum R\'enyi entropy by optimal thermodynamic integration paths

Despite being a well-established operational approach to quantify entanglement, R\'enyi entropy calculations have been plagued by their computational complexity. We introduce here a theoretical framework based on an optimal thermodynamic integration scheme, where the R\'enyi entropy can be efficiently evaluated using regularizing paths. This approach avoids slowly convergent fluctuating contributions and leads to low-variance estimates. In this way, large system sizes and high levels of entanglement in model or first-principles Hamiltonians are within our reach. We demonstrate it in the one-dimensional quantum Ising model and perform the evaluation of entanglement entropy in the formic acid dimer, by discovering that its two shared protons are entangled even above room temperature.

Introduction.-Measuring the entanglement of a quantum state or the entropy of a quantum system at thermal equilibrium has always been a challenge. With the aim of achieving this goal, several methods have been proposed so far [1][2][3][4][5][6][7][8][9]. One of the most promising approaches is based upon the evaluation of the quantum Rényi entropy. For the subsystem A of a quantum system, itis defined as [1,[10][11][12] where ρ and ρ A are density matrices of the full system and of its subsystem A, respectively, with α ∈ R >0 \ {1}.
Despite being a fundamental proxy to understand the thermodynamics of quantum systems, the Rényi entropy is not measurable in experimental setups, apart from rare successes [22,26]. The same holds for the Rényi entropy evaluation via analytical treatments or numerical methods, such as density matrix renormalization group (DMRG) and stochastic sampling frameworks [2-5, 9, 18, 27, 28]. In the former case, it is limited mostly to integrable models [14,[29][30][31][32], while in the latter, it is limited to low-dimensional systems with sufficiently low * miha.srdinsek@upmc.fr entanglement [33][34][35]. Quantifying the Rényi entropy in more general complex systems has remained so far an unattainable task, hampered by the exponentially high energy barriers in stochastic sampling frameworks.
Of particular interest are hydrogen-rich materials, like liquid water [36], where quantum effects arise due to the light mass of hydrogen. It has been shown that nuclear quantum effects pilot phase transitions in these systems, leading, for instance, to phase VIII of water ice and to the superconducting phases of hydrides, such as LaH 10 [37], YH n [38], and H 3 S [39]. For the same reasons, entanglement is supposed to play a role also in biochemical systems, like the formic acid dimer [40][41][42][43] and base pairs in DNA [44][45][46].
In this Letter, we present an alternative approach that overcomes previous limitations and allows one to compute the Rényi entropy S α A for complex quantum systems, described by either model or ab initio Hamiltonians. The latter is demonstrated by producing the first evidence of entanglement in the formic acid dimer. Our approach is based upon the combination of the path integral (PI) formalism and thermodynamic integration along appropriately defined paths.
Rényi entropy with path integrals.-In the PI formulations of statistical mechanics, the quantum density matrix ρ at time t is mapped to a classical counterpart by Wick rotating (t → −iβ) the quantum action and discretizing it into a classical Hamiltonian H [47,48]. In practical implementations, the imaginary time interval [0, β = 1/(k B T )], T being the temperature, is divided into a finite number of time steps, which are individual snapshots (or beads) that interact particle-wise in the imaginary time direction. One can then use importance sampling algorithms to evaluate the Rényi entropy of α ∈ N \ {1} by calculating the free energy difference between two statistical ensembles, namely S α A = log(Z A /Z ∅ )/(1 − α) [14]. Z ∅ = [Trρ] α is the partition function of the ∅ ensemble, consisting of α independent Z ∅ and ZA for a quantum particle with α = 2 are shown in the upper boxes. S 2 can then be computed by averaging the "swap" operator or by introducing intermediate steps (grey arrows) and averaging the spring interactions highlighted in colors (red, blue) over an appropriately defined path. The springs strength is proportional to color intensity. copies of the system, and Z A = Trρ α A is the one for the joint ensemble, where each particle belonging to the subsystem A is replaced by one particle living in all the α copies. Thus, Rényi entropy depends only on the free energy cost of changing the boundary conditions through the "swap" operator (see Fig. 1) in the imaginary time direction.
This free energy cost can be estimated by running a PI Monte Carlo simulation in one ensemble, say Z ∅ , and averaging the exponent of the energy difference between the two boundary conditions at given β [2,8,49], in order to evaluate H A,∅ are classical Hamiltonians arising from the discretization of their corresponding quantum actions Z A,∅ . Though Eq. 2 is in principle exact, it implies that for increasing energy differences one has to wait exponentially longer times to sample high-energy configurations. A possible remedy for this slowing down is to split the calculation in shorter increments of smaller free-energy differences [2], and to sample the ratio of the Metropolis-Hastings transition probabilities min(1, exp(−β∆H)) Z ∅ / min(1, exp(β∆H)) Z A [50,51], rather than Eq. 2. Following these procedures, reliable and groundbreaking entanglement entropy estimations were obtained in one-dimensional (1D) and twodimensional (2D) spin chains [3,9], and cold atoms [27]. Nevertheless, in larger or more strongly entangled systems, the consequently higher energy barriers can cause the aforementioned approach to fail. Thermodynamic integration.-An alternative scheme is the thermodynamic integration based on a new partition function Z[λ], a differentiable function of the parameter λ ∈ [0, 1], that connects the two ensembles [5,8].
Then, by using the relation:  Fig. 2c). Thus, to get reliable results high precision and a large number of integration steps are needed. The thermodynamic integration based on the simplest integration path is therefore also doomed to failure when faced with large systems [3]. In order to understand the origin of these drawbacks, we focus on the harmonic oscillator and its Rényi entropy of second order (i.e. with α = 2), the exactly solvable minimal model, where the behaviour is reproduced. In the PI formulation of Z ∅ , each copy of the quantum harmonic oscillator of mass m and frequency ω is described by a ring polymer with beads {x k } k=1,...,P connected by harmonic springs [48], such that In the following, we will solve this model for m = 1 and ω = 1. In Eq. 4, ζ = β/P is the imaginary time step which controls the discretization error and represents the inverse temperature of the classical analogue. Switching from the split ensemble to the joint one amounts to neglecting the harmonic interaction in H ∅ between the beads x P and x 1 in both copies, and adding the interaction between x P of the first copy and x 1 of the second one, and its crossed term, in H A (see Fig. 1). The average energy of the joint interaction computed over Z ∅ is greater than the split one, because copies do not interact with each other in Z ∅ and, thus, they can be at relatively large distances. This contribution grows with the square of the inter-bead distance and with the phasespace size. Moreover, as the temperature increases, the interaction energy grows as 1/β and de Broglie wavelength shrinks. Copies collapse to almost point-like particles and this further contributes to the diverging cost of joining them. This shows up in the large value of the integrand at λ = 0 in the simplest integration scheme. Similarly, the variance of this contribution is large, since Z ∅ minimises the action without H A − H ∅ . An analogous contribution but of opposite sign appears at λ = 1, further increasing the overall uncertainty of the integral. Nevertheless, we found that if a more general path is considered, the pathological behavior of the integrand at the edges and its large variance can be avoided all together.
Path regularization.-In search of enhanced paths, we focus on the 2D parameter space (g, h) of Hamiltonians H(g, h) = H ∅ + (g − 1)K ∅ + hK A . The K ∅,A operators correspond to the interactions that enforce the boundary conditions: K ∅ (K A ) drives the intra (inter) -copy closure of the particle rings. These terms are depicted in red (blue) in Fig. 1. Different paths connecting (1, 0) (H = H ∅ ) to (0, 1) (H = H A ) can be compared by inspecting the gradient field of the energy in the (g, h) plane, shown in Figs. 2(a) and 2(b), respectively. At each (g, h) point, the direction of the field indicates the path p ≡ (g(λ), h(λ)) that would yield the largest possible increment or the largest possible uncertainty to the line integral in Eq. 3, with magnitude represented by the color of the stream plot. For the 1D harmonic oscillator (Fig. 2c), we can clearly see that the path connecting linearly the two endpoints, i.e. the one that has commonly been employed so far, produces two spikes in ∂ λ H(λ) Z[λ] at (1, 0) and (0, 1), exactly where the scalar product of the path direction (p ≡ ∂ λ p) with the gradient field is the largest. This is understood once the integral in Eq. 3 is recast in −β 1 0 p · Kdλ. The origin of spikes can be traced down to the K A growth along (g, 0) and K ∅ growth along (0, h), as the temperature increases. These drawbacks are present in any system where intra-bead interaction is large. In order to remove the spikes, one should therefore avoid moving towards directions where previously non-existing interactions are switched on in H(g, h).
This analysis suggests that the optimal integration path is the one which minimizes over the entire path. Another criterion can be derived by minimizing the full variance [53], which implies the line minimization of Either choice makes the integrand flat as a function of λ. This scheme acts therefore as a path regularization. Not only the spikes at the endpoints are cut off, but the line  (2)) computed with thermodynamic integration via path regularization (red crosses) and transition probability sampling based on the swap operator (blue crosses) compared with analytical results (dashed-black line) for L = 64 and β = 3. Inset: Comparison of system-size scaling for both methods at different r.
integral can also be computed on a much coarser grid, speeding up the calculation and reducing both deterministic and stochastic errors.
It is clear that the full path optimization would not be affordable in the most general case. However, the simple 1D quantum oscillator, where the path search can be carried out systematically, provides a shape that is transferable to complex quantum many-body systems. p can then be parametrized as differentiable curve: If rescaled by the proton mass, the reported behavior of the 1D harmonic oscillator spans a physically relevant range of parameters with temperatures going from 4.3 K (β = 40) to 344 K (β = 0.5) and with a vibrational frequency of 5122 cm −1 . In this realistic regime, it turns out that l = 2 optimizes the path based on |p ·K| ( Fig. 3(a)), while l = 3 is the optimal power law based on var[K] ( Fig. 3(b)). Nevertheless, the latter is the best choice in ab initio systems over a large range of temperatures [53]. Indeed, these systems share a similar intra-bead interaction to the one of the harmonic oscillator studied above, provided by the quantum kinetic term of the ab initio action.
1D Ising model in transverse magnetic field.-The harmonic oscillator is however a too simple model to possess any entanglement. In a system with a larger number of interacting particles, the coefficients (g, h) change the boundary condition of the full subsystem. To explicitly check the generality of the path regularization approach, we now study the 1D Ising model in transverse magnetic field with periodic boundary conditions. Its Hamiltonian reads where σ x,z i are Pauli matrices acting on i-th site, and r the strength of the magnetic field in the x direction. Due to its integrability (it can be analytically solved using Jordan-Wigner transformation [54]), it represents an ideal benchmark for our path regularization in an extended system. It undergoes a quantum phase transition at r = 1.
Within the PI framework, Eq. 8 is mapped into the 2D classical anisotropic Ising model [53]. The interaction in the imaginary time direction of the classical counterpart diverges as r → 0, causing effects analogous to the temperature increase in the harmonic oscillator. At small r, the path in Eq. 7 successfully kills large fluctuations appearing at the endpoints if l = 3 is used. In this case, there is still room for improving the thermodynamic integration path. Indeed, in this model system, the magnitude of the endpoints still grows with the interaction strength. This can be cured by rescaling the λ parameter [53], which provides an additional freedom to optimize the path. Rescaling λ leaves the shape of the path unchanged but it varies the "speed" (i.e. the integration points density) along the thermodynamic trajectory. To compare our method with the sampling algorithms based on the swap operator (Eq. 2), like the one of Ref. 3, we computed [55] the Rényi entropy of second order S 2 full with the whole system as a subsystem. It is the hardest quantity to evaluate. The full entropies obtained with the swap scheme and the path regularization both agree well with analytical results over a wide range of magnetic field strengths (Fig. 4). S 2 full nicely captures the quantum phase transition at r = 1, by showing a clear peak. In order to estimate the maximum system size that the path regularization procedure can afford, we pushed the entropy calculation to very long spin chains, where the thermodynamic limit is reached. In the limit of large L, the full entropy scales as L. From the inset of Fig. 4, it can be seen that this limit is reached at relatively small system sizes. It is also seen that our procedure outperforms the swap-based one, since the former can still be applied to systems larger than 600 sites where S 2 full exceeds the value of 50, while the latter is broken already before 400 sites [3]. Indeed, one of the strengths of the path regularization method is that the number of integration steps remains constant with the subsystem size. By increasing the level of entanglement, the time cost grows linearly, while in methods based on Eq. 2, it grows exponentially.
A realistic system: the formic acid dimer.-By means of the path regularization scheme, the calculation of the Rényi entropy for realistic systems becomes feasible. As an illustrative example, we take here the formic acid dimer, a minimal model of biochemical physics. This is a system of two molecules that form a dimer via a double hydrogen bond (Fig. 5c) Due to the 180 • rotation symmetry around the axis connecting the carbon atoms, the potential energy surface (PES) has two minima, which correspond to the two hydrogen configurations shown in Fig. 5. They are separated by a barrier that grows with the inter-dimer distance d between the oxygen atoms (d = 2.7Å in equilibrium [41,56]), which determines also the hydrogen bond stretch. During the double proton transfer the two molecules get closer, up to d = 2.4Å, and the barrier dwindles. Due to the light hydrogen mass, at intermediate distances quantum effect become prominent [40][41][42][43] and the barrier is low enough that the two configurations are expected to be entangled, leading to a proton concerted motion.
In our simulations, we restricted the protons to move along the hydrogen bond and we evaluated the PES as a function of d by means of the coupled cluster single double perturbative triple (CCSD(T)) method [53]. At each distance, we ran a PI Monte Carlo simulation of the two protons in a projected PES with our thermodynamic integration scheme [55]. As the potential is simple enough, the outcome of these simulations can be compared against results obtained with the exact diagonalization in a discretized space. This comparison further assesses the robustness of the path regularization in Eq.7, with the optimal power of l = 3 for ab initio systems. Since the Rényi entropy of the 2-proton subspace ("full" entropy) is nonzero at distances where the system is in a mixed state, we can extract a lower bound for quantum correlations [57] by computing the difference between the full and the Rényi entropy of the single proton subspace (entanglement entropy if in a pure state). The results show that the entanglement is present along the full range of distances, explored during the double proton transfer (d = 2.4 − 2.7Å [41]). However, the temperature plays a major role. At room temperature the entanglement is considerably lower compared to the one at 100K. Remarkably, it persists up to temperatures as large as 500K for some intermolecular distance. We note that at high temperatures the thermal motion of the full molecular complex must be taken into account. Nevertheless, the effect is so strong that the entanglement should still be relevant well above room temperature even in this case.
Conclusions.-In this Letter, we introduced a path regularization scheme that allows an efficient and stable calculation of the Quantum Rényi entropy via a thermodynamic integration performed within the path integral framework. The proposed regularization defines an optimal thermodynamic path that smoothly changes the imaginary time boundary conditions of the quantum partition function, by avoiding slowly convergent contributions and by yielding a low-variance estimate of the Rényi entropy. The method has been shown to be efficient in the 1D Ising model with transverse magnetic field, where we reached very large subsystem sizes. It also allowed us to perform a first ab initio evaluation of the Rényi entropy for the concerted hydrogen motion in the formic acid dimer. The path regularization makes the evaluation of the Rényi entropy feasible for model and real systems, comprising large sizes and/or complex interactions, for which the Rényi entropy analysis was previously inaccessible.

ACKNOWLEDGMENTS
We thank the support of the HPCaVe computational platform of Sorbonne University, where the main calculations have been performed. We are grateful for the environment provided by ISCD and its MAESTRO junior team. This work was partially supported by the European Centre of Excellence in Exascale Computing TREX-Targeting Real Chemical Accuracy at the Exascale, funded by the European Union's Horizon 2020 Research and Innovation program under Grant Agreement No. 952165.

QUANTUM RÉNYI ENTROPY BY OPTIMAL THERMODYNAMIC INTEGRATION PATHS SUPPLEMENTAL MATERIAL
The supplemental material is organized as follows. In Sec. I we introduce the computational details for the quantum Harmonic oscillator, in Sec. II we explain how we simulated the one-dimensional Ising model, while in Sec. III we detail some numerical aspects of the ab initio calculations of the formic acid dimer.

I. HARMONIC OSCILLATOR
The one-dimensional (1D) quantum harmonic oscillator was used as the minimal model that retains the pathological behavior of the integrand at the edges of the thermodynamic integration in Eq. 3 of the main text. The expression in Eq. 3, proportional to the Rényi entropy, can be recast in the following form: where the components of K ≡ ( K ∅ λ , K A λ ) are the energies of the springs connecting the beads at the boundaries of ∅ and A ensembles, respectively, the quantum averages in K are taken over the the distribution generated by H(g, h), and the integral is performed over the path p : To find the optimal shape of the path that best regularizes the integrand in this system, we evaluated three cost functionals: In the above expressions, |p |dλ is the length of the full path, p 2 ≡ ((∂g/∂λ) 2 , (∂h/∂λ) 2 ), and var is the variance of both K components.
Since the integration in Eq. S1 depends only on the initial (g, h)(0) = (1, 0) and final (g, h)(1) = (0, 1) points of the path p, the behavior of the integrand can be tested with the functional F Abs in Eq. S2 and Fig. S1a. This functional will disclose paths that lead to cancellation between large contributions of opposite sign. Indeed, it will reach its minimal (and optimal) value when its positive and negative contributions are minimized in absolute value. With the aim at performing an efficient integration, not only the value but also the shape of the integrand is relevant. Indeed, an oscillating integrand requires more integration steps. The functional F Slope in Eq. S3 and Fig. S1b sums up the absolute value of the integrand derivative along the path and reaches the minimum when the integrand is linear. Finally, the functional F Variance in Eq. S4 and Fig. S1c is the variance of the integral, assuming no covariance between K ∅ and K A , for simplicity. It originates from the fluctuations in the stochastic evaluation of K, which accumulate along the path. It is an estimate of the stochastic rather than deterministic uncertainty in the computation of the integral.
The smallest variance at high temperature (β = 0.5) is reached when l = 3 (Fig. S1c), but all the paths with l > 1 are substantially better than the linear one (with l = 1). Very similar observations can be made for the slope functional (Fig. S1b). This outcome is expected, given the behavior of the integrands shown in Figs. 2c-2e of the Letter. Indeed, at l = 3, the integrand is almost linear at low temperatures. While F Slope and F Variance agree in yielding l = 3 as best path, F Abs favors l = 2 (see Fig. S1a), because negative contributions to the integral almost vanish in that case. In F Abs , the value of l = 3 is, however, not much worse compared to l = 2, particularly at high temperature. Thus, l = 3 is the best compromise between the different cost functionals presented here and, therefore, it is the optimal power law to be used in the regularizing path, not only for the harmonic oscillator but also for more challenging realistic systems.

A. Quantum-to-classical isomorphism
The 1D quantum Ising model, studied in the Letter, was analysed using Path Integral Monte Carlo simulations with action expanded up to the second order in the Suzuki-Trotter breakup. At this order of the breakup, the model can be formally interpreted as a twodimensional (2D) classical anisotropic Ising model [58] with Hamiltonian of the form where ζ = β/P is the inverse temperature of the classical analogue, P is the number of beads and N is the number of spins in the 1D quantum model. Only when P → ∞ the equivalence is formally exact. For numerical investigations, we restrict ourselves to finite P with error proportional to O(β 3 /P 3 ). The first index refers to the spin position in the quantum chain and the second one to the position in the imaginary-time axis. Each σ i,j is allowed to take values ∈ {−1, 1}. In the original 1D spin model we used periodic boundary conditions, which translate into σ z N +1,j = σ z 1,j in the classical model. r is the strength of the magnetic field in the x direction, and causes the interaction term to diverge when r → 0.

B. Sampling schemes
We evaluated the Rényi entanglement entropy by measuring the free energy cost of modifying the topology on a N × αP grid, by changing the boundary conditions in the imaginary-time direction. In the Z ∅ ensemble, they lead to imaginary-time indices defined as with f : N → N a function such that f (i) = i for i ∈ {1, . . . , P } and f (P + 1) = 1, which is a consequence of the trace. In the joint Z A ensemble, where a subset A of particles is chosen, the boundary conditions stay the same for the complement i ∈Ā, but change for i A ∈ A such that This is a consequence of first raising the reduced density matrix to the power α and only then taking the trace over the subsystem. When the thermodynamic integration is performed, both boundary conditions are used, but the interaction is weighted with parameters g and h, defining the path (see Eq. 5 of the main text).
The most efficient sampling scheme for this system is not the standard Metropolis-Hastings algorithm, but the cluster algorithm adapted to the anisotropic 2D classical Ising model [59], where all proposed moves are accepted. Instead of flipping one spin or the whole row, a cluster of spins is created and each spin added with a probability that fulfills the detailed balance condition. After the cluster is constructed, all its spins are flipped. The algorithm has to be slightly modified for the purpose of thermodynamic integration, where the spins lying on the boundary have three neighbours in the imaginary-time direction. The cluster can grow in both directions there, while the probability of growing in each direction is controlled by the corresponding weights.
Evaluation of full entropy with the "swap" operator is inefficient when entropy is large, but this can be partially cured, by splitting the calculation in parts for L equal to the subsystem size. Every ratio can be computed in a separate simulation and only at the end the product is performed.
In order to make a fair comparison between the thermodynamic integration with the path regulatization and the direct evaluation of the partition function ratios in Eq. S8, the "swap" operator was split in a number of steps equal to the one in the thermodynamic integration, and the same sampling scheme was used, with the same number of cluster algorithm steps. We ran 16 integration steps and 8 intermediate points with "swap" operator, based on Metropolis-Hastings transition probabilities min(1, exp(−β∆H)) Z ∅ / min(1, exp(β∆H)) Z A . The results of this comparison are shown in Fig. 3 of the Letter.
Although it is not really necessary for the simulation of this system, a multi-walker algorithm was used, where each walker follows its own Markov chain, and contributes to the total average. The final error is evaluated as the standard deviation of the averages given by each walker.

C. Different integration paths in 1D quantum Ising model
The Ising model substantially differs from the ab initio Hamiltonian already by the fact that spins are allowed to take only two values. This bounds the interaction in the imaginary-time direction and the divergence like the one observed in the harmonic oscillator does not appear. However, the interaction strength still diverges with r → 0, and a move connecting two split ensembles is too costly to be accepted. The regularizing path introduced in the Letter is therefore relevant again, because it reduces significantly the variance and allows one to produce accurate results up to very large system sizes. Although the path regularization with l = 3 leads to good results with low variance in the 1D Ising model (Fig. S2b), one is still left with the freedom of adjusting the density of points along the path, based on some quadrature rule. Thus, more points can be added close to the difficult endpoints, resulting in a slower pace along the path and, hence, reducing the weight of integration steps with high-density points. This reshapes the integrand and can further reduce the integration error [49]. The drawback of this approach is that a denser integration grid could lead to an increased computational cost, if there is no particular gain in making the integrand smoother.
In our analysis of the Ising model, λ was rescaled according to  from the unfavorable shape of the integrand (Fig. S2a). We note that the path regularization proposed in the Letter, based on a flexible path direction in the (g, h) Hamiltonian space, is already sufficient to reduce the variance.
Rescaling the parameter λ according to Eq. S9 on top of the path regularization can merely be used to further enhance the integration efficiency by reducing the number of integration steps required, as shown in Fig. S2a. The need of path regularisation is demonstrated in Fig. S3, where linear interpolation (l = 1) is compared with the regularised path. It is also apparent that rescaling λ according to Eq. S9 is less effective than changing the shape, i.e. tuning l, in order to regularize the thermodynamic integration.

III. FORMIC ACID DIMER
A. Potential energy surface evaluation In the formic acid dimer simulations, we explored the possibility of evaluating quantum entanglement in a truly ab initio system. It is made of of two molecules that form a dimer via a double hydrogen bond. The Rényi entanglement entropy calculation was based upon the determination of the potential energy surface (PES) for the two protons involved in the bond. We constrained them to move only along the hydrogen bond direction, and we fixed the positions of all the other atoms, thus reducing the PES to two dimensions.
The position of the other atoms was obtained by optimizing the geometry of the system in a symmetric configuration, i.e. at the transition state, where the two protons are equally shared. This could affect the results at larger distances, but we expect these effects to be small close to the symmetric transition state. After fixing the H-C<OO configuration for each monomer, only the distance between them is varied and, at each distance, the PES is evaluated. The PES constructed in this way has a global minimum close to the one of the transition state (see Fig. S4).
The PES was computed with the coupled cluster single double perturbative triple (CCSD(T)) method, using the standard library implemented in the Python-based Simulations of Chemistry Framework (PySCF) package [60,61] with the 'STO-3G' basis set. Because the method is time consuming, the PES is first determined on a 30×30 points grid, and then Bspline-interpolated on a denser grid of 10, 000 × 10, 000 points.
The PES dependence on the distance d between faceto-face oxygen atoms (a proxy for the distance between the two monomers) is shown in Fig. S4. With increasing d, the PES minimum splits in two, and entanglement between the protons drastically increases. The configurations with maximally entangled protons are not the configurations with the lowest energy.

B. Path integral Monte Carlo
Rényi entanglement entropy was evaluated using two methods: Monte Carlo sampling, used to get the main results reported in the Letter, and exact diagonalization, used for benchmarking when accessible.
For the simulation of the ab initio Hamiltonian only a simple algorithm was used in path integral Monte Carlo, consisting of two types of moves. The first move consists in displacing a randomly chosen pair of beads with a random vector drawn from the uniform distribution of width ζ 2 /(2m). The second possible move is a mirror transformation, obtained by changing the sign in front of the particles position of one copy. The move is always accepted for particles that are not connected, while for the connected particles the only energy difference involved comes from the boundary. This second type of move was necessary for sampling the double well potential, where particles get stuck in one of the wells if the barrier is too high. In this case, ergodicity would not be preserved. This additional Monte Carlo move was enough to pre-serve ergodicity in our simulations.

C. Exact diagonalization
In order to evaluate the entanglement entropy via the exact diagonalization, we used the Hamiltonian truncation method. With this method the infinite Hilbert space is truncated in such a way that the space becomes finite and the low-energy properties stay the same. One can test the correctness of the results by gradually increasing the space and by looking at the eigenvalues. Since the Hilbert space of the system is a product of two oneparticle states, V 1 ⊗ V 2 , the space can be truncated by introducing only a finite number of possible positions for each particle. In the position space, the kinetic term is represented by finite differences while the potential is diagonal in the basis. We performed calculations on grids spanning from 30 × 30 to 100 × 100 points. For the potentials considered, small grids were already extremely accurate.
Once the low-lying eigenvalues are known (at low temperature the high-energy contributions can be neglected), the evaluation of the full entropy is straightforward. However, for the entanglement entropy the partial trace has to be performed. Because we know how the space is constructed, the trace can be expressed as x ,x y c i,x ,y c * i,x ,y |x x |, (S10) where i runs over the eigenstates, c i,x ,y are the coefficients that determine the eigenstates expansion over the position basis set: |v i = x y c i,x,y |x ⊗ |y . The partial trace requires just the sum over y. After this is computed for each pair x , x , a reduced density matrix is obtained. The remaining object to calculate is the trace of ρ A raised to the power α which can be done without any diagonalization.