Entanglement and classical fluctuations at finite-temperature critical points

We investigate several entanglement-related quantities at finite-temperature criticality in the three-dimensional quantum spherical model, both as a function of temperature $T$ and of the quantum parameter $g$, which measures the strength of quantum fluctuations. While the von Neumann and the R\'enyi entropies exhibit the volume-law for any $g$ and $T$, the mutual information obeys an area law. The prefactors of the volume-law and of the area-law are regular across the transition, reflecting that universal singular terms vanish at the transition. This implies that the mutual information is dominated by nonuniversal contributions. This hampers its use as a witness of criticality, at least in the spherical model. We also study the logarithmic negativity. For any value of $g,T$, the negativity exhibits an area-law. The negativity vanishes deep in the paramagnetic phase, it is larger at small temperature, and it decreases upon increasing the temperature. For any $g$, it exhibits the so-called sudden death, i.e., it is exactly zero for large enough $T$. The vanishing of the negativity defines a"death line", which we characterise analytically at small $g$. Importantly, for any finite $T$ the area-law prefactor is regular across the transition, whereas it develops a cusp-like singularity in the limit $T\to 0$. Finally, we consider the single-particle entanglement and negativity spectra. The vast majority of the levels are regular across the transition. Only the larger ones exhibit singularities. These are related to the presence of a zero mode, which reflects the symmetry breaking. This implies the presence of sub-leading singular terms in the entanglement entropies. Interestingly, since the larger levels do not contribute to the negativity, sub-leading singular corrections are expected to be suppressed for the negativity.


I. INTRODUCTION
Understanding the interplay between entanglement and classical fluctuations in quantum many-body systems is a challenging task. Several entanglement-related quantities have been investigated [1][2][3][4]. Arguably, the most popular ones are the Rényi entropies and the von Neumann entropy. Given a partition of the system in subsystems A and its complement B ≡Ā (see Fig. 1 (a) for a three-dimensional setup), the Rényi entropies S n are defined as with ρ A the reduced density matrix of part A. In the limit n → 1 one obtains the so-called von Neumann entropy S = −Trρ A ln ρ A . If the total system A ∪Ā is in a pure state, S n are good measures of the entanglement between A and the rest. At zero temperature the famous area-law behaviour S n ∝ L d−1 holds true, with L being the linear size of the system, and d the dimensionality. Logarithmic corrections are present for gapless fermionic systems [5][6][7]. This is strikingly different at finite temperature, where S n cannot distinguish genuine quantum correlations from thermal ones. At finite temperature the volume-law behaviour S n ∝ L d holds, reflecting, for instance, the transformation between entanglement and standard thermodynamic entropy. In this situation, the mutual information I n,A∶Ā (see section II for its definition) can be used to measure the total correlation, i.e., both classical and quantum, between A andĀ. The mutual information exhibits the area-law behaviour I n ∝ L d−1 . Still, I n is not a proper measure of the mutual entanglement between A andĀ. This reflects A ∪Ā being in a mixed state. The so-called logarithmic negativity [8][9][10] is the most popular entanglement measure for mixed states. For systems described by Conformal Field Theory, its behaviour has been fully characterized, both at zero temperature [11], and at finite temperature [12]. Similar to the mutual information, for any temperature the negativity exhibits the area-law scaling, which has been checked in several systems [13][14][15][16][17]. An important question, which is key in this paper, is how quantum and classical fluctuations are intertwined at a finite temperature critical point. An important remark is that typical local one dimensional models do not exhibit finite-temperature or mid-spectrum criticality, with the notable exception of models that exhibit the many-body localisation transition [18]. It has been argued that entanglement-related quantities can be used to detect this transition [4,[19][20][21], although no conclusion has been reached, due the severe finite-size effects. On the other hand, although standard finite-temperature phase transitions are possible in higher dimensions, the interplay between classical and quantum fluctuations at finite-temperature criticality has been scarcely studied, despite the growing interest [22][23][24].
Much efforts focused on the behaviour of the mutual information. For instance, it has been suggested [25][26][27][28][29][30] that the mutual information exhibits a crossing for different sizes at a finite-temperature critical point, similar to more traditional tools in critical phenomena, such as the Binder cumulant [31].
Importantly, there is growing interest to understand the behaviour of the logarithmic negativity at finite-temperature criticality. An important question is whether the logarithmic negativity exhibits signatures of criticality. Since finitetemperature phase transitions are driven by classical fluctuations, and the negativity is only sensitive to true quantum correlations, one should expect the negativity to be smooth across a finite-temperature critical point. Surprisingly, it has been observed that the negativity can exhibit "weak" singularities, such as cusp-like features, at the critical point [16,17,32]. On the other hand, similar singularities have been also observed in the area-law prefactor of the Rényi entropies at quantum phase transitions [33,34].
In this paper we investigate the behaviour of several entanglement-related observables at finite-temperature criticality. Specifically, we consider the Rényi entropies (and the von Neumann entropy), the mutual information, and the logarithmic negativity. We also discuss the so-called single-particle entanglement spectrum and the negativity spectrum. To be specific, here we focus on the paradigmatic quantum spherical model (QSM) in three dimensions. Spherical models have a venerable history, as they served as test ground for the theory of critical phenomena and the theory of finite-size scaling [35,36]. The phase diagram of the model (see Fig. 2) is spanned by the temperature T and the so-called quantum coupling g. The latter measures the strength of the quantum fluctuations. In three dimensions, the model exhibits a line of finite-temperature second-order phase transitions between a paramagnetic phase and a ferromagnetically ordered phase. The universality class of the transition has been characterised analytically, since the model is exactly solvable, and it is that of the N -vector model [31] in the limit N → ∞ limit [37]. Surprisingly, the entanglement properties of spherical models have not been explored much.
First, we verify that the entanglement entropy and the Rényi entropies exhibit the expected volume-law in the whole phase diagram. In the thermodynamic limit, the prefactor of the volume-law and its first derivative with respect to the model parameters and the temperature are continuous across the transition. We refer to this behaviour as regular. Conversely, we refer to non-analytic behaviour as singular. This means that, although the entanglement entropies contain singular contributions, these vanish at the transition as g − g c κ (or as T − T c κ ) with κ > 1, and T c and g c marking the critical point. The mutual information between two regions exhibits the expected area-law behaviour for any value of g and T . Moreover, although the prefactor of the area law contains both regular and singular contributions at criticality, again, we observe regular behaviour across the transition, suggesting that the singular terms vanish. Interestingly, for finite-size systems we observe that the mutual information does not exhibit a universal crossing at the transition. This reflects that the singular terms, which contain universal information, are vanishing at the critical point, and the mutual information is dominated by non-universal contributions. We find a similar behaviour for the logarithmic negativity. The negativity exhibits the expected area-law in the whole phase diagram. The area-law prefactor has a maximum at the phase transition. In the paramagnetic phase the negativity vanishes as 1 g in the limit g → ∞. In the ordered phase it depends mildly on the system parameters for low enough temperature and large enough values of g. For any fixed value of g, the negativity exhibits a "sudden death" upon increasing the temperature. Specifically, there is a negativity "death line", above which the negativity is exactly zero. Interestingly, in the limit g → 0, i.e., when the model becomes classical, the death line exhibits the behaviour ∝ g 1 2 . Importantly, at finite temperature the prefactor of the area-law is regular across the para-ferro transition, whereas it develops a cusp-like feature upon lowering the temperature.
Finally, we discuss the single-particle entanglement spectrum [4] and the negativity spectrum [38]. For two nearestneighbour sites in an infinite system, the spectrum contains only two levels. One of the levels is regular across the transition, whereas the other exhibits a cusp singularity. This is different for the negativity between two extended subsystems. Specifically, we observe that the majority of the single-particle entanglement spectrum levels show regular behaviour at the transition. Only the larger levels exhibit a cusp-like behaviour. Interestingly, the corresponding eigenvectors show a flat-in-space structure, suggesting that the singular levels originate from the zero mode, which is associated with the symmetry breaking. As it is well known, this zero mode leads to sub-leading logarithmic terms [39,40] in the entropies. This also means that, although the prefactor of the volume law behaviour of the entropies is regular across the transition, the sub-leading corrections exhibit signatures of the criticality. We observe a similar structure for the negativity spectrum. However, since the larger negativity spectrum levels do not contribute to the logarithmic negativity, this suggests that the sub-leading singular terms are suppressed.
The manuscript is organised as follows. In section II we introduce the entanglement-related quantities that we consider in the paper. In section III we review the quantum spherical model. Specifically, in section III A we provide the results for the two-point correlation functions, discussing their behaviour in the paramagnetic phase, and in the ordered phase (in section III B and section III C, respectively). In section IV we address the critical behaviour of the model. Criticality is encoded in the so-called spherical parameter that we discuss in section IV A. In section IV B we focus on the correlators near the para-ferro transition. In order to compare entanglement-related quantities with standard thermodynamic ones, in section IV C and section IV D we discuss the singular behaviour of the free energy, Periodic boundary conditions are used in all directions. In (b) the total system is divided into two parts as A ∪ B, and B is traced over. Here x × y × z is the volume of A. In (c) subsystem A is further divided as A = A1 ∪ A2. Here we consider the entanglement between the two adjacent cubic blocks A1 and A2. We consider both the case of finite B, as well as the limit of B infinite.
and of the universal ratio R ξ . In section V we focus on the interplay between entanglement and finite-temperature criticality. In section V B and section V C we study the entanglement entropy, and the mutual information. In section VI we present our results for the logarithmic negativity. In section VI A we focus on the negativity between two sites embedded in an infinite system, whereas in section VI B we discuss the negativity between two extended regions. In section VII we investigate the single-particle entanglement spectrum and the negativity spectrum. Finally, we conclude in section VIII.

II. ENTANGLEMENT ENTROPIES, MUTUAL INFORMATION & LOGARITHMIC NEGATIVITY: DEFINITIONS
In this section we introduce the definitions of the entanglement and Rényi entropies, the mutual information, and the so-called logarithmic negativity. Here, we always consider a system on a three-dimensional cubic lattice of linear size L (see Fig. 1 (a)), prepared in a thermal state at temperature T ≡ 1 β.
In order to define the entanglement entropy (von Neumann entropy) and the Rényi entropies, we consider a bipartition of the system into two parts A and B (see Fig. 1 (b)). From the reduced density matrix for A, ρ A ≡ Tr B ρ, with ρ being the density matrix of the full system, we define the Rényi entropies S n,A as [1][2][3][4] In the limit n → 1 one obtains the von Neumann entropy as Both, the von Neumann and the Rényi entropies, are good entanglement measures for pure systems, for instance bipartite systems at zero temperature. This is not the case anymore if A ∪Ā, withĀ the complement of A, is in a mixed state, for instance a thermal one. In a thermal state, S n,A exhibits a volume-law behaviour S n,A ∝ L d . For S vN , the prefactor of the volume-law is the same as that of the thermal entropy, i.e., the von Neumann entanglement entropy becomes the thermodynamic entropy at finite temperature [41][42][43]. The mutual information between two subsystems is defined for the typical geometry depicted in Fig. 1 (c). The subsystem A is divided into two complementary parts viz A 1 ∪ A 2 = A, withĀ in general not empty. The mutual information I n,A1∶A2 is defined as Here S n,X is the Rényi entropy of subsystem X. It is important to note that the mutual information is only a measure of the total correlation shared between the two subsystems, and not of the shared entanglement. This is because I n,A1∶A2 is sensitive to both, quantum and classical correlations, which is reflected in the fact that A is in a mixed state. Importantly, this is also the case at T = 0 if A 1 ∪ A 2 is not the full system, i.e., if A 1 and A 2 are not complementary intervals. In this case the classical correlation originates from the trace over the environment, i.e.,Ā in (4). It can be shown, however, that the mutual information provides a bound to the mutual entanglement [9]. Both the entanglement entropy and the mutual information have been studied intensely as possible witnesses of critical behaviour in several systems (see Ref. 4 for a review), both at finite-temperature phase transitions [25,26,28,29,[44][45][46], or at T = 0 phase transitions [27,30,33,34,47]. Specifically, in Ref. 26 (see also 25) it has been shown that the mutual information densities I n,A1∶A2 L d−1 exhibit a universal crossing at T = T c and T = nT c . On the other hand, at zero temperature, it has been shown that the area-law coefficient of the entanglement entropy displays a cusp-like singularity across a second order phase transition (see for instance Ref. 33).
Here we also consider the so-called logarithmic negativity [8][9][10][48][49][50]. Unlike the mutual information, the negativity is a good entanglement measure for mixed states. The negativity allows, for instance, to quantify the entanglement in a bipartite system at finite temperature (see Fig. 1 (b)), or the entanglement between two non-complementary subsystems (as in Fig. 1 (c)) at zero temperature. The negativity is defined from the so-called partial transpose. Given the partition of A as A = A 1 ∪ A 2 (see Fig. 1 (c)), the matrix elements of the partial transpose ρ T2 A with respect to the degrees of freedom of A 2 are defined as Here {ϕ 1 } and {ϕ 2 } are two orthonormal bases for A 1 and A 2 , respectively. Interestingly, the eigenvalues ζ i of ρ T2 A can be both positive and negative. This is in contrast with the eigenvalues of ρ A , which can be only positive. The logarithmic negativity E A1∶A2 is defined as Recently, the negativity has become a useful tool to characterise universal aspects of quantum many-body systems [11,12,32,38,, also because it can be computed with Matrix Product States (MPS) methods [51,55,57]. The negativity (6) can be obtained for free-bosonic models in arbitrary dimension by correlation matrix techniques [13,14,79]. This is not possible for free-fermion models [80][81][82][83][84][85][86][87]. An alternative entanglement measure, which is effectively calculable using free-fermion techniques, has been introduced [15,85,86,[88][89][90], and it is also an upper bound for the negativity [87]. Very recently, much attention has been focused to study the behaviour of the negativity at a finite-temperature phase transition. It has been suggested in Ref. 16 and Ref. 17 that in some cases the logarithmic negativity exhibits a cusp-like singularity, i.e., it is sensitive to the classical criticality, although it is a measure of quantum entanglement.

III. QUANTUM SPHERICAL MODEL
The quantum spherical model [91,92] (QSM) is defined on a 3D cubic lattice of volume V = L 3 , with L the lattice linear size (see Fig. 1). The Hamiltonian reads In Eq. (7), n = (n x , n y , n z ) denotes a generic lattice site, and ⟨n, m⟩ a lattice bond joining two nearest-neighbour sites.
Here s i and p i are canonically conjugated variables satisfying the commutation relations In Eq. (7), µ and g are real parameters. The first term in Eq. (7) is a kinetic term, which makes the model quantum.
The parameter g is the quantum coupling. For g = 0 the model becomes classical and it reduces to the famous spherical model [93,94]. The spherical parameter µ is fixed by imposing the spherical constraint as where ⟨⋅⟩ denotes the average over the thermal ensemble. The additive shift by 3 in the definition of the spherical parameter in (7) is solely for later convenience and routinely not performed in literature. Critical properties (for instance, critical exponents) of the QSM are determined by the behaviour of µ. In order to diagonalise H, one exploits the translational invariance of the system by defining the Fourier transformed operators π k and q k as Here the sum over k ≡ (k x , k y , k z ) runs in the first Brillouin zone k i ≡ π Lj, with j ∈ [−L, L] integer. In terms of q k , π k the Fourier representation of the Hamiltonian becomes The single-particle dispersion Λ k reads To completely diagonalise (11) we introduce ladder operators b k and b † k as with α 2 k = g 2Λ −1 k . The operators b k obey standard bosonic commutation relations. By using (13), the Hamiltonian (11) becomes diagonal, and it is given as In an equilibrium thermal ensemble, the modes k are occupied according to the Bose-Einstein distribution where β = 1 T is the inverse temperature. The quantum spherical model has been studied extensively in any dimension. In three dimensions, the model exhibits a finite-temperature second-order phase transition between a high-temperature paramagnetic phase and a low-temperature ferromagnetic (ordered) phase. The universality class of the transition is the same as that of the classical spherical model, and it has been fully characterised [35] (see also [36]). The universality class is the same as that of the N -vector model at N → ∞. At T = 0 the model undergoes a second-order quantum phase transition at a critical g c . The universality class of the transition is the same as that of the classical spherical model in 3 + 1 dimensions [95], as expected from renormalisation group arguments. Critical properties of the model are determined by the behaviour of µ. For any temperature, µ ≥ 0 for g > g c . On the other hand, one has µ = 0 for g < g c , i.e., in the ordered phase, which signals non-analytic behaviour. The phase diagram of the model is reported in Fig. 2. The continuous line is the critical line marking the second-order phase transition between the ferromagnetic phase at small g and low temperature, and the standard paramagnetic phase.

A. Two-point correlation functions
The key ingredients to study entanglement-related quantities in the quantum spherical model are the two-point correlators of the operators s n and p n (cf. (7)) at equilibrium. They can be readily obtained by first expressing s n , p n in terms of b k , b † k that diagonalise the model, and using (15). One obtains [96] ⟨s n s m ⟩ = 1 2V k e i(n−m)⋅k α 2 k coth(βE k 2) (16) Here we defined k ≡ (k x , k y , k z ). Note that in the thermodynamic limit, at the critical point one has that µ = 0, and the contribution of the zero mode with k x = k y = k z = 0 diverges. This is not the case at finite L, because for a finite system µ is nonzero for any g and β (see section IV A). Clearly, one can rewrite (16) as where F −1 d (x) denotes the inverse Fourier transform of x in d dimensions. Eq. (19) is more suitable than (16) for numerical computations because there are very efficient methods for evaluating the Fourier transform.
From (16), the constraint (9) for the spherical parameter µ reads The free energy F of the QSM reads as The expressions for (16)(17)(18) and (21) in the thermodynamic limit L → ∞ are obtained, as usual, by replacing In the following sections we discuss the behaviour of the correlators (16)(17)(18) in the paramagnetic phase at large g, and in the low-temperature ordered phase.

B. Paramagnetic phase: Large g expansion
Here we discuss the large g expansion of the correlators in the paramagnetic phase, in the thermodynamic limit. To this purpose, one has to first determine the behaviour of the spherical parameter µ, which is obtained by solving Eq. (30), in the large g limit. We numerically verified that µ → ∞ for g → ∞. Moreover, by taking the limit g → ∞ in (30), and using the definition of E k (cf. (11)), one obtains that the leading behaviour of µ is µ + 3 = g 8 + o(g). To derive the higher order corrections, it is natural to conjecture that µ has the following expansion The coefficients C 1 and C 3 can be determined by substituting the ansatz (24) in (30), expanding for g → ∞, and equating the terms with the same power of g. After doing that, and after neglecting exponentially suppressed terms, one obtains Higher order terms can be obtained in a similar way. The behaviour of the correlation functions (16) and (17) at large g is easily obtained after substituting (24), and expanding at large g. This yields Here δ n−m ,1 = ∏ j=x,y,z δ nj −mj ,1 . Clearly, longer-range correlations are suppressed with higher powers of 1 g, as expected because the model is paramagnetic, and correlation functions decay exponentially at large distances.

C. Ferromagnetic phase: Low-temperature expansion
We now discuss the behaviour of the model at low temperature in the thermodynamic limit. In the ordered phase one has µ = 0. In the limit β → ∞ one can replace coth(βx) → 1 in (16) and (17). The next-to-the-leading behaviour at large β is obtained by using the standard saddle point method. One obtains the expansions for the correlators (16) and (17) as As it is clear from (28) and (29), the leading behaviour of ⟨s n s m ⟩ is determined by the integral of 1 √ ω k . The first ). Note that the sub-leading corrections do not depend on n, m. A similar behaviour occurs at criticality (see section IV), and it has important consequences for the singularity structure of entanglement-related quantities.

IV. CRITICAL BEHAVIOUR OF THE QUANTUM SPHERICAL MODEL
Here we are interested in the critical behaviour of the quantum spherical model at finite temperature. The main goal of this section is to derive the behaviour of the two-point correlators in the vicinity of the para-ferro transition (see Fig. 2). To do that we first derive the expansion for the spherical parameter near the critical point (see section IV A). In section IV B we present our results for the correlators. In order to compare the behaviour of entanglement-related quantities at the finite-temperature phase transition with that of standard quantities, in section IV C and section IV D we discuss the scaling of the free energy and the universal ratio R ξ at criticality.

A. Spherical parameter
The critical behaviour of both the classical and the quantum spherical model is determined by the spherical parameter µ (cf. (9)). It is somewhat easier to work in the thermodynamic limit, although the finite-size behaviour of µ can be derived using standard techniques [97]. In the thermodynamic limit, Eq. (9) becomes (see section III A) To proceed let us first rewrite (30) as with ω k as defined in (12). Near criticality, on the paramagnetic side, one has µ → 0, whereas µ = 0 everywhere in the ordered phase. This reflects the presence of singular terms in the expansion of µ near the transition. To derive these Spherical parameter µ near the para-ferro transition on the paramagnetic side. In (a) and (b) we approach the transition at fixed temperature and fixed g, respectively. The symbols are exact numerical results in the thermodynamic limit obtained by solving (30). The lines are the analytical results near the critical point (cf. (38) and (40)).
terms at the leading order in g − g c , it is convenient to expand (31) for small µ + ω k . The reason is that the singular terms are determined by the singularity at small k of the integrand in (31). One has Here a(µ, β, g) is an analytic function of its arguments. The first few terms of its series expansion read where the dots denote higher order terms. Importantly, at the leading order in µ one has a k = O(µ). The first term in (33) encodes the critical behaviour of the model. The integral is the celebrated Watson integral [98]. The same integral appears in the classical spherical model, reflecting that for nonzero g the universality class of the transition is the same as that of the classical model. It is interesting to observe that the integral can be expressed explicitly in terms of hypergeometric functions [99]. For instance for µ = 0 one has To proceed, one subtract from (33) the expansion of the equation for the spherical parameter (30) at g = g c . One The leading behaviour in µ of the integral in (35) is obtained by expanding ω k for k → 0. One now has In the last steps we changed variables to spherical coordinates with y ≡ k , performing the integration on a sphere of radius π √ µ, instead of the cube [−π √ µ, π √ µ]. This introduces a O(µ) term. Since the leading order is O( √ µ), this is negligible in the limit g → g c . Note that the term √ µ is singular at g c .
Since we are interested in the leading O( √ µ) behaviour, and a k (µ, β, g) = O(µ) we can set µ = 0 in a k . We also observe that At the leading order in g − g c , from (33)(35)(37) we obtain It is also convenient to expand the integrands in (38) at small k, to obtain In a similar way, one can derive the expression for µ if the transition is approached at fixed g by varying the temperature. One obtains Again, after expanding for small k, one has Note that the expected behaviours [95] µ ∝ (g − g c ) 2 and µ ∝ (β − β c ) 2 hold. Finally, a similar calculation [97] should allow, in principle, to extract the finite-size behaviour of µ.
The correctness of (38) and (40) is verified in Fig. 3. The symbols in the figure are the values of µ obtained by solving (30). The dashed-dotted lines are the analytical results (38) and (40). Note that (38) and (40)  It is also useful to check how the thermodynamic limit is approached. In Fig. 4 we report numerical results for µ for finite-size systems. These are obtained by solving numerically (20). In the figure we plot 3 + µ versus g. The data are for fixed temperature T = 1, although similar results are obtained for different temperatures. The vertical dotted line denotes the critical coupling g c . The dashed-dotted line is the result for µ in the thermodynamic limit, i.e., obtained by solving (30). Clearly, µ = 0 for g ≤ g c , whereas µ → ∞ in the limit g → ∞.
The different symbols are the finite-size results for µ. One obtains a nonzero value for µ for any g. As it is clear from the figure, both in the paramagnetic and in the ferromagnetic regions the data quickly converge to the thermodynamic limit (dashed-dotted line). Around the critical point, large finite-size effects are present. The perturbative result for µ (cf. (24)) is reported in Fig. 4 as dashed-dotted line. Although the agreement is not perfect for the values of g reported in the figure, we checked that the result (24) is recovered upon increasing g.

B. Two-point correlators
By using the expansion for µ near the critical line that we derived in the previous section, it is now straightforward to obtain the two-point correlation functions. Again, the idea is to expand ω k around k = 0 in (16) and (17) in the thermodynamic limit. For the case in which the transition is approached at fixed T , one obtains ω k e ik(n−m) If the critical point is approached at fixed g, one has For the correlators ⟨p n p m ⟩, a similar calculation yields If one approaches the transition along the temperature direction, one obtains Clearly, in all cases (cf. (42)(43)(44) (45)) the correlation functions are finite at the critical point, and exhibit a O(g − g c ) and O(β − β c ) behaviour. Crucially, in both (42) and (43) one has the singular contribution ∝ √ µ. On the other hand, this is not present in the expansion of ⟨p n p m ⟩. This implies that the correlator ⟨s n s m ⟩ has a cusp-like singularity across the critical point, whereas ⟨p n p m ⟩ is regular. Finally, one should observe that the singular term in (42) and (43) does not depend on the position. This means that the spatial structure of the correlators appears only in sub-leading contributions in g − g c that we are neglecting.

C. Scaling of the free energy
It is instructive to investigate the singular contributions to the free energy. A similar qualitative behaviour will be observed for entanglement-related quantities (see section V B, section V C and section VI). Let us focus on the situation in which T is kept fixed, and let us study the behaviour around the critical point as a function of g. At a second order phase transition the free energy density contains both analytic and singular terms. In general it can be written as [31] f ≡ F L 3 = f an (g) + f sing (g).
The term f an is the smooth part of the free energy, which depends analytically on the system parameters. The second part f sing contains the singularities and the universal properties of the model. f sing obeys the scaling form (see Ref. [31] for a review) Here b > 0 is an arbitrary number, and u n are the scaling fields, which are analytic functions of the system's parameters. y n are the scaling dimensions associated with the fields u n . Here we assume that there are only two relevant scaling fields, u 1 and u L . u 1 is associated with changing the coupling g, and u L with the finite-size scaling. u L has scaling dimension y L = 1. We define y 1 = 1 ν. In (47), y n < 0 for n > 1, i.e., u n are irrelevant, which give non-analytic scaling 18  corrections to the free energy. Several scaling laws can be derived from (47). For instance, the finite-size scaling form of the free energy is obtained by choosing b = 1 u L . We obtain Using that u L ≈ 1 L, and that close to the phase transition u 1 ≈ g − g c , we find the standard result For a second-order phase transition, the free energy density is finite everywhere in the phase diagram. The singular part f sing vanishes at the critical point. Indeed, for g = g c one obtains from (49) the contribution f sing (0)L −3 , which is vanishing because f sing is finite. By choosing u 1 b y1 = 1 in (47) in the thermodynamic limit one obtains the scaling behaviour Note the dependence on the sign of g − g c . In (50) we used the hyperscaling relation dν = 2 − α, with α the critical exponent of the specific heat. In the limit L → ∞ one has to recover (50) from (49). This implies that f sing ((g − g c )L 1 ν ) ≈ (g − g c )L dν f sing (±∞) for L → ∞. We should stress that in the derivations above we neglected all the scaling corrections. It is interesting to derive the singular behaviour of the free energy in the quantum spherical model. We work in the thermodynamic limit. The density of free energy of the QSM reads (cf. (21)) As for the spherical constraint and for the correlators (see section IV A and section IV B), the idea is to expand (51) for small k. After expanding (51), one obtains where b k denotes an analytic function of its arguments. At the leading order it is given as The leading singular behaviour is encoded in the first term in (52). To extract the singularity one can use the trivial identity Now, the integration over dk in (52) can be performed exactly to give 1 2β Here I 0 is the modified Bessel function of the first kind. Now one can split the integration range as [0, ∞) = [0, t 0 ] ∪ [t 0 , ∞). For t 0 large enough, by using the asymptotic behaviour of the Bessel function I 0 (t) ∝ e t √ 2πt, one obtains The singularity of the free energy is extracted from the second term in (56). In the limit µ → 0 this gives Here Γ(0, t 0 ) is the incomplete Gamma function. Note that both analytic and nonanalytic terms are present in (57). The term µ 3 2 gives the well-known singularity of the free energy in the quantum spherical model [95] f sing = − From (58) we obtain ν = 1. The prefactor of (g − g c ) 3 can be easily extracted from (57) and (38) (40). Eq. (58) can be also derived by observing that the spherical constraint (see section IV A) is obtained as ∂f ∂µ = V , and by integrating with respect to µ the singular contribution in (36). From (58), one has that the specific heat exponent α = 2 − dν = −1 is negative. Eq. (58) implies that the free energy and its first derivative with respect to β and g are continuous functions at the critical point.
We illustrate the behaviour of the free energy across the finite temperature transition in Fig. 5. In the two panels (a) and (b) we consider two different limits. Specifically, in (a) we calculate the free energy density from the finitesize expression (21), but using the value of µ calculated in the thermodynamic limit, i.e., by solving (30). This is convenient because it is straightforward to numerically solve (30), whereas solving (20) is a nontrivial task since the sum over k cannot be performed explicitly. On the other hand, since in the thermodynamic limit one has µ = 0 below the transition, the contribution of the zero mode in (21) is divergent and it has to be regularized by hand. The strategy that we use is to introduce a small mass putting µ = 10 −6 . While this ensures that one has the correct thermodynamic limit results, it affects the sub-leading contributions. This is clear from Fig. 5 (a). In the figure we show data for T = 0.6 and several system sizes. The data for L = 4 show a jump at g = g c (vertical line). This is an artifact of the regularization. We checked that the jump becomes sharper and sharper as the mass term is sent to zero. At fixed mass, upon increasing L, the contribution of the zero-mode becomes negligible and the data approach the thermodynamic limit result (dashed-dotted line in the figure). Clearly, in the thermodynamic limit the free energy and its first derivative are continuous, as expected from (58). In Fig 5 (b) we show the free energy calculated by using the finite-size value of µ, obtained by solving numerically (20). Now the contribution of the zero mode is finite because of the finite L.

D. Universal ratio R ξ
One important goal of this paper is to investigate the effectiveness of entanglement-related observables to detect finite-temperature criticality. In this section we briefly review the behaviour of the universal ratio ξ 2nd L, with ξ 2nd the second-moment correlation length. This is a standard tool used in numerical simulations to analyse criticality [31]. For instance, it allows to detect second-order phase transitions via the so-called crossing method, and it can be used to extract the critical exponent ν by the usual data collapse analysis. The ratio R ξ is defined as The so-called second-moment correlation length ξ 2nd is extracted from the long-distance behaviour of the correlation function. Its definition reads  withG(k) the Fourier transform of the spin-spin correlation function (cf. (16) for the spherical model). In (60), q min is the minimum nonzero lattice momentum q min = (2π L, 0, 0). Note thatG(0) is the spin susceptibility χ defined as Alternatively, instead of ξ 2nd , one can define R ξ by using ξ gap , which is defined from the gap between the ground state and the first excited state in the energy spectrum. For the quantum spherical model R ξ can be expressed analytically as a function of µ. By using (16), one obtains that the correlation length is given as We now discuss the finite-size scaling properties of R ξ . The finite-size scaling ansatz for R ξ reads as [31] R ξ = K((g − g c )L 1 ν ).
Here we neglect both analytic and nonanalytic scaling corrections, that, however, one can include. The function K(x) is universal apart from a renormalisation of its argument. Assuming analyticity of K(x) one can expand (63) as The value R * ξ depends only on the universality class of the transition and on the geometry (for instance, the boundary conditions). Crucially, Eq. (64) implies that the curves for R ξ at different finite sizes L exhibit a crossing at the critical point. Moreover, when plotted against the scaling variable X ≡ (g − g c )L 1 ν , the data for different Ls collapse on the same curve, at least in the limit L → ∞, when corrections to scaling can be neglected.
It is important to understand the behaviour of R ξ , especially with the moderately small system sizes that are available in our simulations. In Fig. 6 we show numerical data for R ξ plotted as a function of g. The data are for T = 1. The curves for different L exhibit a crossing at g ≈ g c , as expected. The value R * ξ ≈ 0.5 at the crossing is universal, and it could be calculated by using the results of Ref. 97. To our knowledge R * ξ is not known exactly. Note that scaling corrections are present. Indeed, for the smaller Ls the crossing is not at g c (vertical line). Upon increasing L, the crossing point exhibits a systematic drifts towards g c . In Fig. 7 we show the same data as in Fig. 6, now plotted versus the scaling variable X = (g − g c )L 1 ν . As expected, due to the scaling corrections, the smaller lattice sizes do not show data collapse. However, upon increasing L the quality of the collapse improves significantly. The scaling is satisfactory for the larger sizes. In the inset we provide data also for T = 2.2. We observe that scaling corrections are smaller as compared with T = 1, which is reflected in a better data collapse. Note that the scaling function K (cf. (63)) depends on the temperature, although the value at X = 0 is the same for both temperatures, as expected. However, since the finite-temperature universality class is the same on the whole para-ferro transition line (see Fig. 2), the two scaling functions should coincide after an analytic redefinition of the scaling field u 1 as u 1 ≈ c T (g − g c ), where c T depends on the temperature. We numerically verified that c T =2.2 c T =1 ≈ 1 2.

V. ENTANGLEMENT SCALING IN THE CRITICAL SPHERICAL MODEL
We now discuss the behaviour of entanglement-motivated quantities in the finite-temperature critical QSM. Specifically, in section V A we briefly review how to calculate entanglement-related observables. In section V B we focus on the von Neumann entropy. In section V C we discuss the von Neumann and Rényi mutual information. In section VI we investigate the the logarithmic negativity. We consider both the negativity between two spins in an infinite system (in section VI A), as well as the negativity between two extended adjacent blocks (in section VI B). Finally, in section VII we focus on the single-particle entanglement spectra and negativity spectra.

A. Computation of entanglement-related observables in the QSM
The quantum spherical model is mappable to a system of free bosons (see section III). This implies that entanglement-related quantities are obtained from the two-point correlation functions (see Ref. [100] for a review).
The key ingredients are the matrices Q nm ≡ ⟨s n s m ⟩ and P nm ≡ ⟨p n p m ⟩ constructed from the two-point correlation functions (cf. (16) and (17)).
Let us define the matrices Q[A] and P[A], where n, m are now restricted to subsystem A (see Fig. 1 (b) and (c)). The Rényi entropies of A are constructed from the eigenvalues λ 2 j of the matrix Q[A] ⋅ P[A]. As it is common in the literature, we refer to the λ j as the single-particle entanglement spectrum levels. In terms of λ j , the Rényi entropies are given as The von Neumann entropy is obtained by performing the analytic continuation n → 1, which yields The mutual information (see section II) is obtained by using (65) (66) and (4). We now discuss the logarithmic negativity for a partition of A as A = A 1 ∪ A 2 (see Fig. 1 (c)). One first defines the transposed matrix P[A T2 ] as Here the matrix R[A T2 ] acts as the identity matrix I A1 on A 1 and as −I A2 on A 2 . The eigenvalues ] form the single-particle negativity spectrum. In terms of ν 2 i , the negativity is given as Note that while ν 2 i > 0, − ln(2ν i ) can be both positive and negative.

B. Von Neumann entropy
Let us start discussing the behaviour of the von Neumann entropy. Numerical data for the half-system entropy are reported in Fig. 8, for fixed temperature T = 2.2. Similar to the free energy (see Fig. 5) and to the thermal entropy, the von Neumann entropy exhibits a volume law at any g. The figure shows the density of entropy S vN L 3 , plotted as a function of g. Different symbols are for different system sizes. Finite-size effects decay very quickly with L. The data for the larger sizes L = 12−16 collapse on the same curve. The entropy density exhibits regular behaviour around the transition. This is expected because at finite temperature the density of von Neumann entanglement entropy becomes the same as the thermal entropy, which is not singular at the transition. We should stress, however, that sub-leading contributions can be singular, for instance due to the presence of the zero mode. A natural scaling ansatz for the entropy density reads as Here the functions s sing and s an encode the singular and regular terms. In (69) we neglect scaling corrections. Eq. (69) implies that the singular term vanishes at the transition when increasing L, and only the regular term survives. We should stress that in (69) we also neglect logarithmic contributions that can arise because of the presence of the zero mode [39], which reflects the symmetry breaking. We will investigate these contributions in section VII discussing the single-particle entanglement spectrum. Moreover, in principle, there can be extra logarithmic corrections if the bipartition has corners (see, for instance Ref. 30). These are not present in our case because the boundary between A and its complement is smooth (see Fig. 1).

C. Mutual information
We now turn to the mutual information (see (4) for its definition). Here we consider a bipartite system at finite temperature (as in Fig. 1 (b)). The system is divided into two equal parts A and B =Ā. We consider the mutual information between A and its complement. As it is clear from (4), the volume-law contribution of the entropies  cancels out. This cancellation happens for any value of g and T . Thus, the mutual information exhibits an area law. A natural scaling ansatz for the density of the mutual information is Eq. (70) is compatible with the scaling ansatz proposed in Ref. 25. Here ν is the critical exponent as in (69). Similar to (69), here we are neglecting scaling corrections. In (70), q an is the analytic contribution. We present our results for the von Neumann mutual information I vN in Fig. 9. We plot the density of mutual information I vN L 2 versus g. Data are at fixed T = 2.2. The data exhibit a crossing point around g ≈ 13.5, which is incompatible with the critical point at g c ≈ 14.77. Moreover, the position of the crossing point does not change upon increasing system size, in contrast with the behaviour of R ξ (see Fig. 6). Importantly, the curves for different Ls do not "fan out" as L increases, in contrast with the behaviour for the ratio R ξ (see 6). As for the von Neumann entropy, the mutual information is dominated by the regular part (cf. (70)). More precisely, although the singular term in (70) gives a universal crossing, this is not visible because it is suppressed as L −2 in the limit L → ∞.
Similar behaviour is observed for the Rényi mutual information. This is discussed in Fig. 10 focusing on I 2 . As for I vN , the data in Fig. 10 collapse on the same curve upon increasing L. Specifically, in the paramagnetic phase for g > g c and close to the critical point, finite-size effects decay dramatically with L, and the data with L ≳ 10 are Survey of the behaviour of the logarithmic negativity in the quantum spherical model. On the x-axis g is the quantum coupling. On the y-axis T is the temperature. The figure shows the density plot of the half-system negativity. The system is a cube of length L = 2. The dashed-dotted line is the critical line dividing the paramagnetic phase from the ordered phase at low temperature. The dotted line is the "death line" of the negativity. Above the death line the negativity is exactly zero. The continuous line is the death-line calculated from the negativity between two adjacent spins embedded in an infinite system. The behaviour as ∝ g 1 2 at small g is also reported (dashed line). already indistinguishable from the thermodynamic limit. Also no universal crossing is visible within the system sizes presented in the Figure. It is interesting to compare our results with Ref. 25. There is large evidence, based on quantum Monte Carlo simulations, that at a finite-temperature phase transition the ratio I 2 L d−1 exhibits two crossing at T c and 2T c , with T c the critical temperature. The scaling ansatz for the mutual information presented in Ref. 25 is compatible with (70). However, as stressed in Ref. 25 the presence of the crossing relies on q (n) an changing sign across the phase transition. Our results suggest that this does not happen in the quantum spherical model.

VI. LOGARITHMIC NEGATIVITY
We now discuss the logarithmic negativity between two complementary subsystems. This allows us to study the interplay between genuine quantum fluctuations and thermal fluctuations. One of the goals of this section is to map out the role of entanglement in the different regions of the phase diagram of the quantum spherical model.
The generic behaviour of the negativity is illustrated in Fig. 11. The figure shows a density plot of E as a function of g and T . The dashed-dotted line is the critical line dividing the paramagnetic phase from the ordered phase (see Fig. 2). In the figure we show the half-system negativity for a cube of linear size L = 2. The data are obtained by using the value of µ in the thermodynamic limit and the finite-size formulas for the correlators (cf. (16)(17)). The negativity is large at the quantum critical point and in the ordered phase, and it quickly decays upon increasing the temperature. The dotted line is the negativity "death line". Above the death line the half-system negativity is exactly zero. The death line that we report in the figure is obtained by considering the half-system negativity in the limit L → ∞. In the figure we also report the death line calculated from the negativity between two adjacent sites embedded in an infinite system (continuous line), which provides only a bound. Interestingly, the death line exhibits the behaviour ∝ g 1 2 at small g. For the case of two sites embedded in the infinite system the precise behaviour can be calculated analytically and it is reported in the figure (dashed line). Finally, it is interesting to observe that the death line increases with g, although the negativity decreases as 1 g (see section VI A 1). Logarithmic negativity between two nearest-neighbour spins embedded in an infinite system. We show E along the critical T = Tc(g) line of the para-ferro transition. Note the sudden death of the negativity at g ≈ 13.6.

A. Two-site negativity
Several of the generic features of the negativity discussed in Fig. 11 can be extracted by considering two adjacent spins embedded in an infinite system. In this section we derive analytically the behaviour of the logarithmic negativity in this situation.

Large g expansion
It is straightforward to derive the behaviour of the negativity between the two sites in the large g regime (see Fig. 2). The starting point of the analysis are the correlators ⟨s n s m ⟩ and ⟨p n p m ⟩ in the large g limit. These are reported in (26) and (27). Now the matrices Q[A] and P[A T2 ] (see section II for their definitions) are two-by-two matrices. The negativity spectrum contains only two levels. These are the eigenvalues ν 2 ± of Q[A] ⋅ P[A T2 ], which in the large g limit are given as Since ν + > 1 2, only ν − contributes to the negativity (see (68)), which is given as By expanding in the large g limit, one obtains the behaviour E ∝ 1 g.

Low-temperature expansion
It is also interesting to discuss the limit of low-temperature in the ordered phase of the quantum spherical model. To do that, we exploit the expansion of the correlators ⟨s n s m ⟩ and ⟨p n p m ⟩ (cf. (28) and (29)). Let us first consider the negativity between site n ≡ (n x , n y , n z ) and m ≡ (m x , m y , m z ). The eigenvalues ν 2 i of Q[A] ⋅ P[A T2 ] are given as Logarithmic negativity between two adjacent spins in the quantum spherical model: Negativity E plotted as a function of g (in (a)) and of the temperature T (in (b)). The figure shows results in the thermodynamic limit. The vertical lines denote the critical point. Data in (a) are for fixed T = 1. In (b) we fixed g = g(Tc = 1) (see Fig. 2). Note in (a) the slow decay at g → ∞ and the "death" of E at g → 0. Note that in (b) the negativity is identically zero for T ≳ 3.
Here we defined the Watson-type integrals W ′ n as Here ω k is defined in (12). Interestingly, the integral W ′ n can be calculated analytically in terms of hypergeometric functions [98]. We now restrict ourselves to the negativity between two nearest-neighbour spins, i.e., with n − m = 1. Specifically, we choose n = (0, 0, 0) and m = (1, 0, 0, ). One can numerically check that ν 2 1 < 1 4, implying that ν 1 does not contribute to the negativity (see (68)). On the other hand, for a given β, one has that ν 2 contributes to the negativity only for sufficiently large g. Specifically, we observe that the condition g > g * (β) has to hold, with the "critical" value g * given as where W ′ 0 ≡ W ′ (0,0,0) , W ′ 1 ≡ W ′ (1,0,0) , and similarly for W ′′ . For each value of temperature, the negativity is exactly zero for g < g * (T ). Note in (76) the behaviour g * ∝ 1 β 2 . This behaviour persists when considering the negativity between two extended systems, although the prefactor in (76) is different. In particular, for two extended systems one obtains a larger value of g * in (76), as it is clear from Fig. 11.
Note that for large enough temperature the negativity vanishes even on the critical line (see Fig. 2). This is shown explicitly in Fig. 12. The figure shows E calculated on the critical line T c (g), with T c the critical temperature at the given value of g. The negativity is exactly zero for g ≲ 14, whereas it is finite nonzero up to g ≈ 19, which corresponds to zero temperature.
The behaviour of the negativity between two nearest-neighbour sites as a function of both temperature and quantum coupling g is summarized in Fig. 13. Panel (a) shows E versus g, at fixed T = 1. In the figure we plot the negativity and not its density. The vanishing behaviour below g * ≈ 2, as predicted by Eq. (76), is clearly visible, as well as the slow decay as 1 g in the paramagnetic phase. Panel (b) shows E at fixed g = g c (T = 1) as a function of T . The negativity exhibits a weak dependence on temperature in the ordered phase, whereas it suddenly drops to zero in the paramagnetic phase (sudden death).

Critical region
An important feature in Fig. 13 (a) is that the negativity exhibits regular behaviour across the para-ferro transition. On the other hand, recently it has been observed that the negativity can exhibit a cusp-like singularity at a finite temperature phase transition [16,17]. Moreover, it has been suggested that the singular part of the negativity obeys the scaling form as with α the specific heat critical exponent. Our results are consistent with Ref. 16  FIG. 14. Single-particle negativity spectrum for two adjacent spins the embedded in an infinite system. We plot − ln(2νn), with ν 2 n the eigenvalues of the correlation matrix, versus g. Data are for fixed T ≈ 1.486. The spectrum consists of two eigenvalues shown in (a) and (b). Note that in (b) − ln(2νn) < 0, implying that the level does not contribute to the negativity. Note the kink at the critical point (vertical line). In contrast, the eigenvalue in (a) is regular, implying that the negativity is not singular at the critical point.
singularity. On the other hand, for α = 0, which is the case investigated in Ref. 17, one has the cusp-like behaviour as To clarify the behaviour of E for two nearest-neighbour sites, it is useful to consider the single-particle negativity spectrum. Using the expansions for the correlators near the critical point (cf. (43)(42)(44), and (45)), it is straightforward to obtain the spectrum levels. These are given as Here we defined Note that Z sing contains the singular term √ µ that we derived in section IV A. Interestingly, this affects only ν 2 , whereas ν 1 is regular.
In Fig. 14 we show the levels of the single-particle negativity spectrum (cf. (80) and (81)). In the figure we plot − ln(2ν n ) versus g. In (a) we show the regular eigenvalue ν 1 , whereas ν 2 is reported in (b). The singularity as g −g c in (b) is clearly visible. Importantly, one has that − ln(2ν 2 ) < 0, implying that ν 2 does not contribute to the logarithmic negativity. This suggests that in the case of two extended subsystems only a subset of the single-particle negativity spectrum levels will exhibit singular behaviour. These levels, however, do not contribute to the negativity, which is regular at the finite-temperature transition. This point will be better clarified when we will discuss the negativity spectrum of two extended regions in section VII.

B. Half-system negativity
We now discuss the half-system logarithmic negativity. Our results are reported in Fig. 15. The figure provides an overview of E at fixed temperature T = 1 and T = 0.2 ((a) and (b), respectively) as a function of g. In both panels we plot the negativity density E L 2 of half system. The different lines are for different sizes L. Similar to section IV C, we calculate the negativity by using the finite-size result for the correlators (see (16) and (17)), and the value of The effect of the system size L is small, especially for T = 1. The qualitative behaviour is similar to Fig. 13 (compare (a) and (b)). For T = 1 the negativity exhibits regular behaviour at gc, whereas in the limit T → 0 a cusp-like feature appears at gc. In (b) the inset shows the data for L = 10, zooming around the critical point.
the spherical parameter µ in the thermodynamic limit. To regularize the divergent contribution of the zero mode we fix µ = 10 −6 for g ≤ g c . We also checked that the results for the negativity do not depend on the choice of the regularization. Clearly, in Fig. 15 finite-size effects are small, especially for T = 1 (see (a)). In (a) the data for L = 8 cannot be distinguished from the result in the thermodynamic limit. In both (a) and (b), the negativity has its maximum value at the critical point. In the paramagnetic phase the negativity decays as 1 g, similar to the case of two spins (see section VI A). Within the ordered phase E exhibits a mild dependence on g, except at small g, where it drops dramatically. Specifically, for g < g * ≈ 2 the negativity is exactly zero. Note that g * decreases with decreasing the temperature (compare (a) and (b) in the figure). This behaviour is similar to what observed for the two-site negativity (see Fig. 13). Again, this implies that for any value of the temperature there is a "critical" value of the coupling g below which the quantum fluctuations are not strong enough to give a finite entanglement negativity. An important remark, however, is that since the negativity gives only a bound on the entanglement, the vanishing of E below g * does not imply the absence of entanglement. From Fig. 15 (a) it is clear that at T = 1, the negativity exhibits regular behaviour at g c . On the other hand, upon decreasing the temperature the negativity develops a cusp. This is clear from Fig. 15 (b). Similar (weak) singular behaviour was observed in the area-law prefactor of the von Neumann entropy in ground-state quantum phase transitions [27,33,34].
It is interesting to investigate the behaviour of the logarithmic negativity at fixed g, i.e., approaching criticality by changing the temperature. This is illustrated in Fig. 16. We plot E L 2 versus T for fixed g = g c (T = 1). Now the negativity exhibits a fast decay with T in the paramagnetic region, and already for T ≈ 3 it is exactly zero. This is in contrast with the result at fixed T (see Fig. 15), where the behaviour as 1 g is observed. Finite-size effects are larger upon decreasing the temperature in the ordered phase. Note that for the smaller system sizes one has a quite large value of E in the limit T → 0, although the area-law behaviour E ∝ L 2 is expected to hold also at zero temperature. Finally, the negativity is not singular at the critical point, as in Fig. 15.
Finally, we investigate the behaviour of the logarithmic negativity in finite-size systems. To address this point, we provide results for the half-system negativity obtained by using the finite-size value of the spherical parameter µ, i.e., by numerically solving (20). Our results are reported in Fig. 17. The figure plots the negativity density E L 2 versus g, for several sizes (symbols in the figure). Similar to the mutual informations (cf. Fig. 9 and 10), the data show small finite-size effects both in the paramagnetic phase and in the ordered phase. At the critical point (vertical line) finite-size effects are larger. However, the data appear to converge to the thermodynamic limit result (dashed-dotted line). The latter is obtained from Fig. 15. Note that the data for different Ls do not exhibit any crossing at the critical point.

VII. ENTANGLEMENT SPECTRA, NEGATIVITY SPECTRA, AND THE ZERO MODE
In section VI A 3 we showed that the negativity spectrum for two nearest-neighbor sites contains one regular and one singular levels. The latter exhibits a cusp-like singularity across the finite temperature transition. The singular level, however, does not contribute to the logarithmic negativity. Moreover, we observed in section VI B that the negativity between two extended blocks shows regular behaviour across the transition, similar to the entanglement plotted as a function of temperature T . Data are obtained by using the value of the spherical parameter µ in the thermodynamic limit. The vertical line is the critical temperature Tc = 1 at gc (see Fig. 2). Note that E = 0 for T ≳≈ 3.5. In the limit T → 0 quantum fluctuations are enhanced and E L 2 increases, although area-law E ∝ L 2 behaviour should persist even at zero temperature. entropy (see Fig. 8). It is interesting to investigate how this is reflected in the negativity spectrum of two extended regions. It is also interesting to compare the singularity structure of the negativity spectrum and the entanglement spectrum.

A. Entanglement spectra
Let us start discussing the entanglement spectrum. Let us consider the two matrices Q nm = ⟨s n s m ⟩ and P nm = ⟨p n p m ⟩ (see section II). Here we work in the thermodynamic limit. Since we are interested in the cusp-like singularity of the entanglement spectrum, it is convenient to expand Q and P as n ) ′ ± the left and right derivative of λ 2 n with respect to g calculated at gc. Data are for = 6. Note that δn ≠ 0 only for a small subset of levels. FIG. 19. Single-particle negativity spectrum between two adjacent blocks embedded in an infinite system. The figure shows the single-particle negativity spectrum levels ν 2 n . Subsystem A is a cube of size . A1 and A2 are the two halves of the cube of linear sizes x = 2, y = , z = (see Fig. 1 (c)). The results are at the critical point at gc ≈ 12.77. Note that only the levels ν 2 n < 1 4 (horizontal dashed line) contribute to the logarithmic negativity. Inset: δn ≡ (ν 2 n ) ′ + − (ν 2 n ) ′ − , with (ν 2 n ) ′ ± the left and right derivative of ν 2 n with respect to g at the critical point gc. Note that δn = 0 for large n, i.e., for the levels that contribute to the negativity.
where the dots denote higher order terms in powers of g − g c , which we neglect. The matrices Q 0 and P 0 are obtained from S (0) nm and P (0) nm (cf. (80) and (82)), whereas Q 1 and P 1 are easily derived from S (1) nm and P (1) nm (see (81) and (83)). The ± in Q ± 1 is to stress that the matrix is different on the two sides of the transition, because of the term Z sing (see (84)), which is only present for g > g c . We now have at the leading order O(g − g c ) Q ⋅ P = Q 0 ⋅ P 0 + (g − g c )(Q 0 ⋅ P 1 + Q ± 1 ⋅ P 0 ) + . . .
By using standard perturbation theory, one obtains that the corrections to the eigenvalues λ 2 n of the matrix Q 0 ⋅ P 0 are given as δλ 2 n = (g − g c )⟨φ n Q 0 ⋅ P 1 + Q ± 1 ⋅ P 0 φ n ⟩, where φ n ⟩ is the eigenvector of Q 0 ⋅ P 0 corresponding to eigenvalue λ 2 n . Clearly, the singularity in the negativity is determined by the second term in (88), which depends on Z sing . Now we observe that Z sing does not depend on the positions n, m. This has striking consequences for the single-particle entanglement spectrum. First, the singular contribution in (88) is given as [δλ 2 n ] sing = (g − g c )Z sing ⟨φ n (1, 1, . . . ) ⊗ (1, 1, . . . ) ⋅ P φ n ⟩. (89) Note that the vector (1, 1, . . . ) ⊗ (1, 1, . . . ) ⋅ P φ n ⟩ is flat, i.e, all the elements are equal. Moreover, we numerically observed that the eigenvector φ 0 ⟩ is approximately flat, i.e., all its components are 1 3 2 , which reflects the presence of a zero mode. This implies that for n = 0 the expectation value in (89) is nonzero. On the other hand, the components of φ n ⟩ with n > 0 are real and are orthogonal to φ 0 ⟩. This implies that the expectation value in (89) approximately vanishes for n > 0. This allows us to conclude that the only few spectrum levels are singular across the transition, i.e., the ones related to the zero mode. We summarize our results in Fig. 18. The figure shows the single-particle entanglement spectrum levels λ 2 n . Here subsystem A is the cube of size (see Fig. 1 (c)) embedded in an infinite system. One has that λ 2 n > 1 4 ∀n. The eigenvalues λ 2 n quickly decay with their index n. For most of the eigenvalues, one has λ 2 n ≈ 1 4. In the inset of Fig. 18 we plot δ n defined as Here (λ 2 n ) ′ ± are the right and left derivatives with respect to g of λ 2 n , calculated at g c . A nonzero value of δ n signals that the level λ n is singular across the transition. The results in the figure are obtained by using (89). Clearly, one has that δ n ≠ 0 only for small n.

B. Negativity spectra
We now discuss the negativity spectrum. Subsystem A is now divided into two parts A 1 , A 2 of sizes x = 2, y = , z = (see Fig. 1 (c)). The partial transposition is performed with respect to A 2 . Our results for the negativity spectrum are shown in Fig. 19. In the main figure we show the eigenvalues ν 2 n of the matrix Q 0 ⋅ P 0 [A T2 ] versus n 3 . For most of the eigenvalues one has that ν 2 n > 1 4, implying that they do not contribute to the negativity (see section V A). As for the entanglement spectrum, the flat vector is an approximate eigenvector of Q 0 ⋅ P 0 [A T2 ], with eigenvalue ν 2 0 . For instance, we numerically checked that for = 4, 6 the eigenvector associated with ν 2 0 has components are in the interval [−0.15, −0.10], clustering around the expected value 1 3 2 = 0.125. This behaviour is reflected in that of δ n . The definition of δ n is the same as in (90) after replacing λ n → ν n . Clearly, the expansions (85) also hold upon replacing P 0 → P 0 [A T2 ] and P 1 → P 1 [A T2 ]. For each negativity spectrum level, δ n is plotted in the inset in Fig. 19. In contrast with the entanglement spectrum (compare with Fig. 18), one has δ n ≠ 0 for a large subset of levels at small n. For instance, for = 6 one has δ n = 0 only for n ≳ 10. As it is clear from Fig. 18 all the levels of the negativity spectrum with δ n ≠ 0 do not contribute to the negativity because they correspond to λ 2 n > 1 4. This suggests a suppression of the logarithmic correction to the boundary-law scaling of the logarithmic negativity, in contrast with the entanglement entropy (see section VII A).

VIII. CONCLUSIONS
We investigated the interplay between entanglement and classical fluctuations at finite-temperature critical points. Specifically, we focused on the three dimensional quantum spherical model, which has a finite-temperature transition between a paramagnetic phase and a ferromagnetically ordered phase. We considered several entanglement-related observables, such as the von Neumann and Rényi entropies, the mutual information, and the logarithmic negativity. In particular, we characterized the behaviour of the logarithmic negativity in all the different phases and at the transition. We also investigated how the behaviour of the entropies and of the negativity is reflected in the singleparticle entanglement spectrum and on the negativity spectrum.
We now mention several important directions for future research. First, it would be important to explore the behaviour of entanglement-related observables in different dimensions. For instance, in d = 4 the para-ferro transition becomes mean field. It would be interesting to investigate how this is reflected in the singularity structure of entanglement. An interesting direction is to investigate the logarithmic negativity at the quantum phase transition in d = 3. Our results suggest that the logarithmic negativity exhibits a cusp-like singularity. It would be useful to understand how this is reflected in the structure of the negativity spectrum. Another interesting direction is to study the crossover from quantum to classical criticality. Moreover, it would be enlightening to generalize the perturbative analysis developed in section VII including higher-order corrections. This would allow, in principle, to characterise the full singularity structure of entanglement-related quantities.
It would be also useful to extend our analysis to other entanglement-related quantities, such as the quantum Fisher information, especially in the light of the results of Ref. 22. Interestingly, the quantum spherical model remains exactly solvable even in the presence of long-range interactions. This opens the possibility of studying entanglement-related quantities in long-range models.
Another important direction is to study the out-of-equilibrium dynamics of entanglement-related quantities after a quantum quench. In recent years, it has been shown [101][102][103][104] that for integrable models it is possible to describe quantitatively the entanglement dynamics by combining integrability with a quasiparticle picture. It would be interesting to investigate the validity of this picture for d > 1. In d = 3 this could also allow to investigate whether the presence of finite-temperature criticality affects the entanglement dynamics. Finally, it has been suggested in Ref. [105] that in the O(N ) model with N → ∞ the steady-state arising after a quantum quench is not described by the so-called Generalized Gibbs Ensemble. The quantum spherical model provides an ideal framework to clarify this issue.