Calculation of R´ enyi entropy in realistic quantum systems Miha Srdinˇ sek 1,2,3 , ∗ Michele Casula 2 , and Rodolphe Vuilleumier 3 1 Institut des sciences du calcul et des donn´ ees (ISCD), Sorbonne Universit´ e, 4 Place Jussieu, 75005 Paris. 2 Sorbonne Universit´ e, Institut de min´ eralogie, de physique des mat´ eriaux et de cosmochimie (IMPMC), CNRS UMR 7590, MNHM, 4 Place Jussieu, 75005 Paris. 3 Processus d’Activation S´ electif par Transfert d’Energie Uni-´ electronique ou Radiative (PASTEUR), CNRS UMR 8640, D´ epartement de Chimie, ´ Ecole Normale Superieure, 24 rue Lhomond, 75005 Paris. (Dated: December 30, 2021) R´ enyi entanglement entropy can be expressed as free energy cost of changing imaginary time boundary conditions in the path integral formulation of quantum mechanics. Despite being a well- established operational approach to quantify entanglement, R´ enyi entropy calculations have been plagued by their computational complexity. We show that R´ enyi entropy can be efficiently evaluated using thermodynamic integration if a regularizing path, which avoids slowly convergent fluctuating contributions, is used. 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 first evaluation of entanglement entropy in the formic acid dimer, by discovering that its two shared protons are entangled even above room temperature. Introduction.— Entanglement is a unique feature of quantum mechanics, which naturally arises in interact- ing many-body systems. To quantify it, one of the most robust ways is based on the R´ enyi entanglement entropy [1–3], defined as S α A = 1 1 − α log Trρ α A [Trρ] α , (1) where ρ and ρ A are density matrices of the full system and of its subsystem A, respectively, with α ∈ R >0 \{1}. Often studied in spin systems, entanglement has been found to play a major role in quantum phase transi- tions [4–12] and thermalization under unitary time evo- lution [13–15]. It has been used to classify different ground states, topological phases[6, 16], and to extract critical exponents. Quantum mechanics spans a wide range of length and energy scales, and - due to the light mass of hydrogen - quantum effects arise also in many hydrogen-rich materials, like liquid water[17]. In par- ticular, it has been shown that nuclear quantum effects pilot phase transitions in these systems, leading, for in- stance, to phase VIII of water ice and to the supercon- ducting phases of hydrides, such as LaH 10 [18], YH n [19], and H 3 S[20]. For the same reasons, entanglement is sup- posed to play a role also in biochemical systems, like the formic acid dimer [21–24] and base pairs in DNA [25–27]. Despite being a fundamental proxy to understand the thermodynamics of quantum systems, the entanglement entropy is largely not measurable in experimental setups, apart from rare successes[28, 29]. The same holds for the R´ enyi entropy evaluation via analytical treatments or nu- merical methods, such as density matrix renormalization group (DMRG). In the former case, it is limited mostly to integrable models[5, 30–33], while in the latter, it is * miha.srdinsek@upmc.fr limited to low-dimensional systems with sufficiently low entanglement [34–36]. Quantifying the R´ enyi entropy in more general or realistic systems has remained so far an unattainable task, hampered by the exponentially high energy barriers in stochastic sampling frameworks. In this Letter, we present an alternative approach that over- comes previous limitations and allows one to compute the R´ enyi entropy S α A for complex model or ab initio Hamil- tonians. Our approach is based upon the combination of the path integral (PI) formalism and thermodynamic integration along appropriately defined paths. R´ enyi entropy with path integrals.— In the PI formu- lations of statistical mechanics, the quantum density ma- trix ρ at time t is mapped to a classical counterpart by Wick rotating (t →−iβ) the quantum action and dis- cretizing it into a classical Hamiltonian H[37, 38]. 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 individ- ual snapshots (or beads) that interact particle-wise in the imaginary time direction. One can then use impor- tance sampling algorithms to evaluate the R´ enyi entropy of α ∈ N \{1} by calculating the free energy differ- ence between two statistical ensembles, namely S α A = log(Z A /Z ∅ )/(1 − α)[5]. Z ∅ = [Trρ] α is the partition func- tion of the ∅ ensemble, consisting of α independent copies of the system, and Z A = Trρ α A is the one for the joint en- semble, where each particle belonging to the subsystem A is replaced by one particle living in all the α copies. Thus, R´ enyi entropy just depends on the freee energy cost of changing the boundary conditions through the “swap” operator (see Fig. 1) in the imaginary time direction. That 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 β[8, 39, 40], in arXiv:2112.14199v1 [cond-mat.stat-mech] 28 Dec 2021