Entanglement scaling at first order phase transitions

First order quantum phase transitions (1QPTs) are signaled, in the thermodynamic limit, by discontinuous changes in the ground state properties. These discontinuities affect expectation values of observables, including spatial correlations. When a 1QPT is crossed in the vicinity of a second order one (2QPT), due to the correlation length divergence of the latter, the corresponding ground state is modified and it becomes increasingly difficult to determine the order of the transition when the size of the system is finite. Here we show that, in such situations, it is possible to apply finite size scaling to entanglement measures, as it has recently been done for the order parameters and the energy gap, in order to recover the correct thermodynamic limit. Such a finite size scaling can unambigously discriminate between first and second order phase transitions in the vicinity of multricritical points even when the singularities displayed by entanglement measures lead to controversial results.


I. INTRODUCTION
Understanding how many-body interacting systems order into different quantum phases as well as the transitions between them remains one of the most challenging open problems in modern physics. Quantum phase transitions (QPTs) are associated with the non-analytical behavior of some observable and/or correlator which can be either local or non local. With the discovery of topological and new exotic phases [2], which fall outside Landau's symmetry breaking paradigm, the use of entanglement to describe quantum matter seems to be, by all means, necessary [3]. Here, we restrict ourselves to a large class of QPTs which, in analogy to their classical counterparts, are signalled by singularities in the derivatives of the (free) energy [4]. In such cases, phase transitions are classified according to the minimum order of the derivative of the ground state energy which is not continuous. Correspondingly, 1QPTs (2QPTs) show singular behavior on the first (second) derivative of the ground state energy. For 1QPTs, the singular behavior translates into an abrupt discontinuity of some local observables while for 2QPTs the order parameters change continuously with a power law scaling. This scaling is characterized by the critical exponents allowing the classification of apparently different 2QPTs into the same universality class. With such definitions at hand it looks straightforward to distinguish if a given QPT is of first or second order. But this is not always the case when finite size effects are present. This question becomes especially relevant in the vicinity of multicritical points where several QPTs of different order coexist in a narrow range of Hamiltonian parameters.
First studies of the entanglement behavior close to a QPT were performed by analyzing bipartite entangle- * Electronic address: abel.yuste@uab.cat ment in simple spin models like the quantum Ising spin-1/2 chain which exhibits a 2QPT. In [5,6], the authors show that, similarly to the ground state energy behavior, the bipartite entanglement between two adjecent spins, as measured by the concurrence (C), is a continuous function of the parameter driving the transition but its first derivative is singular at the critical point in the thermodynamical limit. Moreover, they showed that, for finite size systems, a finite size scaling (FSS) can be performed on the concurrence allowing the correct extraction of the critical exponents corresponding to the Ising transition. These seminal contributions opened a new path in the analysis of QPTs using entanglement. Since then, a large amount of work has been devoted to deepen the connections between quantum information and QPTs, see for instance [7][8][9][10][11][12][13][14][15][16][17]. All these previous results were brought to a more general form in [18] using the Kohn theorem which links the ground state energy of a local Hamiltonian,Ĥ(λ) = kĤk (λ), where eacĥ H k involves as a maximum k-body interactions, with the corresponding k-reduced density matrices. In the usual case of local Hamiltonians involving just two-body interactions of two particles i and j,Ĥ i j , the energy of the ground state |Ψ 0 can be written as where ρ i j is the corresponding reduced density matrix. If the local Hamiltonians,Ĥ i j , depend smoothly on the parameter λ, then a one-to-one correspondence can be made between the singularities of E 0 (λ) arising in a QPT and the singularities of (the matrix elements) of ρ i j . In turn, this fact translates into singularities of entanglement measures, like the concurrence, which depends exclusively on ρ i j . Therefore, ifĤ i j are smooth functions of the λ parameter that drives the phase transition, it generally follows that ∂ λ E 0 ∼ ρ i j . In other words, for a first order quantum phase transition, the singularity on the ground state energy should arise from a singu-larity in ρ i j . As a consequence, the singularity should be evident in any bipartite entanglement measures depending on the reduced density matrix, ρ i j , for example in the concurrence. By the same reasoning, since a 2QPT has a singularity in the second derivative of the ground state energy and ∂ 2 λ E 0 ∼ ∂ λ ρ i j , it follows that a 2QPT should display a singularity in ∂ λ C. This was nicely illustrated in [5,6], demonstrating the fact that quantum phase transitions are related to entanglement. The above results were cast into a theorem in [18] stating that, unless there exist accidental divergences, the order and properties (e.g. critical exponents, critical point) of first and second order quantum phase transitions are signaled by singularities on entanglement measures depending on the appropriate reduced density matrix of the corresponding ground state. The theorem works in both directions, i.e., a discontinuity in e.g. the concurrence signals a 1QPT while a discontinuity/divergence in its derivative signals a 2QPT. The above theorem has some known caveats. For instance, entanglement measures often involve optimization procedures which may lead to accidental nonanalytic behavior. In the spin-1/2 XXZ model, for example, there is a continuous C and a discontinuous derivative along the 1QPT between the ferromagnetic and critical phase [19]. For the same transition, it has been shown that a symmetry breaking in the ferromagnetic phase also modifies the origin of the non-analytic behavior of the concurrence [20]. Also, a non-analytic C in the absence of a QPT for a model with three body interactions has been reported in [21]. Here, we analyze scaling properties of bipartite entanglement measures near multicritical points. Although finite size scaling is a tool to obtain the thermodynamical properties of the system for continuous phase transitions, here we show that such a tool can be employed also in entanglement measures for 1QPTs. Further we demonstrate that when finite size effects are important, it is precisely the scaling of the entanglement measure and not the measure itself which determines the correct order of the transition. This fact is especially relevant for 1QPT crossed in the vicinity of 2QPT and it is in accordance with the recent results reported by Campostrini and coauthors [1] showing that the order parameter of a 1QPT can be continuous for finite systems and admits an appropriate finite size scaling. The paper is organized as follows. In Sec. II we report how the concurrence scales in the spin-1/2 Ising chain in transverse and longitudinal field. In Sec. III we analyze the negativity in a spin-1 model with several 1QPTs. Finally, in Sec. IV we discuss the results and conclude.

II. SPIN-1/2 ISING MODEL WITH LONGITUDINAL FIELD
The first spin model we explore is the spin-1/2 Ising model with a longitudinal field, whereσ α i are the Pauli matrices for spin i and we set J = 1 and B z ≥ 0. In Fig. 1, we provide the phase diagram of the model. For B x = 0, where the system reduces to the integrable Ising model, there is a 2QPT at B z = 1 between the ferromagnetic (B z < 1) and paramagnetic phases (B z > 1). This 2QPT was the first one studied by means of bipartite entanglement, [5]. When B x 0, the system is no longer integrable and we use the Density Matrix Renormalization Group (DMRG) with open boundary conditions (OBC) [22][23][24] and exact diagonalization (ED) calculations. When the system is in the ferromagnetic phase, a 1QPT takes place at B x = 0 between the two ferromagnetic ground states, ferromagnetic ↑ and ferromagnetic ↓. This transition can be detected with a discontinuity in the magnetization, M x = i σ x i , which passes from positive to negative values. For finite systems, numerical calculations in a region sufficiently close to B x = 0 lead to a smooth slope in M x instead of a discontinuity. To deal with this effect, in Ref. [1] a FSS is proposed for transitions driven by a magnetic field, h, with the argument that the relevant scaling variable, κ, should depend on the ratio between the energy contribution of h, and the gap at the critical point, ∆ L , This quantity can be interpreted as a zoom in the h range where the change in the system energy is comparable to ∆ L . When hL ∼ ∆ L , quantum fluctuations mix both ground states and avoid the energy crossing. As a result, the ground state is continuous along the transition and so is the order parameter associated to the transition. This feature is obviously stronger as we get closer to a 2QPT where the correlation length of the system (ξ) diverges, enhancing the nearby finite size effects and thus increasing ∆ L . As we get away from the 2QPT, ξ decreases, ∆ L → 0, and the relevant range of h shrinks, effectively inducing a discontinuity even for small systems. In [1,25], it is shown that both, the first energy gap and the magnetization obey the scaling ansatz, where f ∆ (κ) and f M (κ) are continuous universal functions for all L and B z . For the model in (2), which is integrable for B x = 0, the scaling variable can be defined as [1], where ∆ L ≈ 2(1 − B 2 z )B L z is the first gap for OBC at the critical point, B x = 0, and Here we apply the concepts developed in Ref. [1] to study how entanglement behaves across the same transition. We use, as a measure of entanglement, the concurrence [26], where λ i are the eigenvalues, in decreasing order, of Since C is supposed to be discontinuous at the critical point for 1QPTs, naively we might expect that it will follow a similar scaling behavior as M x , Eq. (5). In Fig. 2, we show C for the two central spins which, if border effects are neglected, will also hold for the rest of neighboring spin pairs. In the left panel, we observe that C is a continuous function with a spike at the critical point which signals a singularity in the first derivative, ∂ λ C. The spike is more pronounced as the transition is closer to the 2QPT at B z = 1.
It is worth pointing out that such behavior bears strong similarities with the geometric entanglement, a collective measure of entanglement indicating how much the ground state differs from a separable state [14]. Notice that for B z > 1, when the system is in the paramagnetic phase independently of the sign of the longitudinal field B x , (as indicated in Fig. 2 for B z = 1.5) , the concurrence becomes an analytical function of B x . In the central panel of Fig. 2, we plot the scaling of the concurrence normalized by its maximum value, C = C/C max , as a function of the scaling variable κ 1 to determine whether C fulfills a similar scaling relation as in Eq. (5). It is enough to investigate what happens for κ 1 > 0 because C(κ 1 ) has even parity, i.e., C(κ 1 ) = C(−κ 1 ). Finally, in panel c) we display the scaling of the derivative of the concurrence normalized by its minimum value ∂ Bx C = ∂ Bx C/[∂ Bx C] min . Interestingly enough, the concurrence does not scale with the fitting parameter κ 1 , but its derivative does. In the central panel, we can see that the data for different B z and L does not collapse in a universal function, while, in the right panel, there is a good data collapse for different values of B z and L. Thus, ∂ Bx C fulfills the scaling anstatz, where g(κ 1 ) is a universal function for any L and B z . We further discuss all these results in Sec. IV.

III. SPIN-1 XXZ MODEL WITH ON-SITE ANISOTROPY
In this section, we extend the study of entanglement along 1QPTs to a spin-1 system. Since C cannot be used to measure the entanglement of two spin-1 particles, we use the negativity, a lower bound of the actual entanglement [27,28]: where ||ρ T B || 1 is the sum of the absolute value of all singular values of the partially transposed density matrix. whereŜ α l are the spin-1 matrices for spin l and D is the uniaxial single-ion anisotropy which we take as positive. We choose this model because of the richness of its phase diagram [29], schematically shown in Fig. 3, with several 1QPTs which we depict with dashed lines. The behaviour of N along the different 1QPTs present in the model can be very different depending on their closeness to a multicritical point which also involves 2QPTs. We start by examining the negativity as a function of the J z for a constant uniaxial field D = 3.5 and L = 8. Two 1QPTs are crossed at such value of D as indicated in the phase diagram by a grey horizontal line (see Fig. 3). The first one corresponds to the transition from ferromagnetic order to the large D phase, which is crossed approximately at J z = −4.2. Another 1QPT appears between large D/Néel at approx J z = 3.8. As clearly shown in Fig. 4, the former phase transition is clearly signalled by a discontinuity in N, whereas the latter shows a smooth slope along the transition.
In Fig. 5, we show in more detail the behavior of N along the ferromagnetic/large D phase transition. This transition is depicted in Fig. 3 by number (1). The results, in this case, are for D = 2 and show that the behavior is the same along all the transition line, even if we are close to the XY phase. In the left panel, we show N very close to the phase transition critical point, J zc , with a very small step in J z of δ = 10 −13 . We can observe that N is still discontinuous even with such a high precision around J zc . In the right panel, we show that there is a neat level crossing between the two-fold degenerated fer- We observe a discontinuity in the 1QPT from ferromagnetic to large D phases and a smooth slope for the 1QPT between the large D and Néel phases. romagnetic ground state, for J z < J zc , and the large D ground state, for J z > J zc . This means that, at J zc , there is a sudden change in the ground state which we can detect with a discontinuity in N. Therefore, in Fig. 5, we observe the expected discontinuous behavior for N along a 1QPT even for a system of just 8 spins without any finite size effects. In Fig. 4, though, we show that a smooth slope in N rather than a discontinuity is observed for the the large D/Néel transition. In Fig. 6, we study with more detail this transition fixing J z = 3.8 and varying D around the critical point D c,L depicted by (2) in Fig. 3. We add the subindex L to indicate that this quantity now depends on the system size. In order to show the behavior for larger systems, we have used DMRG calculations for L = 32, 64, 128 and 150. In the two top panels, we observe that both, N and the staggered magnetization, , change smoothly around the transition point. As L is increased, the slope is more pronounced getting closer to a discontinuity and we need values of D closer to the critical point to observe the continuous slope. For instance, the necessary step in D to observe a continuous behavior is δ ∼ 10 −4 and δ ∼ 10 −6 for L = 64 and L = 150, respectively. Note that in this In the first row, we plot the quantities as a function of D. We observe a smooth slope of both the negativity and the staggered magnetization where the 1QPT is expected. Bottom row, quantites plotted as a function of κ 2 showing the tendency to converge towards a universal function case, as in Sec. II, we are very close to a 2QPT. As we get further from it, for J z 4, this effect becomes less important and both N and M st z are effectively discontinuous. Since the transition is known to be of first order, we propose a similar FSS as in the previous section, defining a relevant scaling variable, κ 2 , as the ratio between the energy contribution of D along the transition and the gap at the critical point, where now, ∆ L is obtained numerically and (D−D c,L )L is a bare estimation for the energy contribution of the parameter D. In the bottom panels, we plot N and M st z as a function of this scaling variable κ 2 . As we can observe, both quantities seem to converge, though not perfectly, towards a universal scaling, as described by (5). It is worth mentioning that (12) is an approximation whereas, in the previous section we had analytic expressions for ∆ L and m 0 in (6). Actually, in Ref. [25], a similar FSS, with nonanalytic expressions, is proposed for the Potts chain with a similar convergence. It seems reasonable, thus, to state that for this 1QPT, when we are close to the 2QPT, N is continuous due to finite size effects and that it obeys the scaling ansatz for 1QPT.
Finally, we want to study N in a similar case as in the previous section for the spin-1/2 Ising model with transverse field, i.e., a 1QPT which is a consequence of a discrete Z 2 symmetry breaking of a two fold-degenerated ground state. In the spin-1 model in (11), this can be done by adding a magnetic field, B z L i=1Ŝ z i , in the ferromagnetic phase leading to a ferromagnetic ↑/ferromagnetic ↓ transition, or a staggered field, B st z L i=1 (−1) iŜ z i , in the Néel phase leading to a Néel/Anti-Néel transition. These regions are highlighted by the ovals along D = 0 in Fig. 3. The order parameter of both transitions, M z and M st z , respectively, are discontinuous in the thermodynamic limit. For the ferromagnetic phase the entanglement remains always constant and zero. More interesting features appear in the Néel/Anti-Néel transition whose results are summarized in Fig. 7. In panel a) and b) we show, respectively, M st z and N as a function of B st z . While M st z has the expected discontinuous behavior, N has a dip at the critical point, signalling a continuous N but a discontinuous derivative. In order to apply a similar FSS ansatz that was used in Sec. II, we define as the relevant scaling variable, Since the 2QPT between Haldane/Néel, which is depicted by dotted lines in Fig. 3, is known to be of the same universality class as the spin-1/2 Ising, we can expect that the analytic function used in (7) is valid now is the 2QPT critical point for D = 0, [30]. In panel c) and d) of Fig. 7, we show, respectively, M st z and N as a function of κ 3 . The scaling ansatz, (5), works fine for M st z , but, as it happened for C in the spin-1/2 Ising model, it fails for N. In panels e) and f) we show how, ∆ L fulfills the scaling ansatz, (4), and that the derivative of the negativity, ∂ B st z N, also verifies the scaling ansatz, (9), which was fulfilled by ∂ B x C in the previous section. To check the correctness of our results, in Fig. 7 we compare the numerical data for M st (κ 3 ) and ∆ L (κ 3 ) the analytic functions derived in Ref. [1] as a result of a two level theory.

IV. DISCUSSION AND CONCLUSIONS
In this work, we have analyzed the behavior of bipartite entanglement measures in diverse 1QPTs. We have focused in 1QPTs which are influenced by a nearby 2QPT. In the first model we have considered, the spin-1/2 Ising model with longitudinal field, the concurrence, C, shows always a continuous behavior along the 1QPT, with a spike signalling a discontinuous first derivative. Following [18], a singularity in the first derivative of C should signal a 2QPT given that this singularity comes strictly from the elements of ρ i j , and not from the computation of C itself. In fact, in Ref. [10], the authors compute C along a 2QPT and obtain a similar behavior to our results in Fig. 2. To determine the origin of this unusual behavior, in Fig. 8, we look into the elements of ρ A,B of the two central spins. In the first column, we plot two of these elements as a function of B x . The range of B x values used, is such that |κ 1 | 1 avoiding, thus, the finite size effects described in Sec. II. The first plotted element, ρ(1, 2) = ↑↑| ρ A,B |↑↓ , is discontinuous along the transition while the second, ρ(1, 1) = ↑↑| ρ A,B |↑↑ , has a spike signalling a singularity in the first derivative. In the second column we show that both singularities become smooth if we look inside the range κ 1 ∼ 1. We have added data for different L and B z to show that these singularities fulfill the scaling for 1QPT described in [1]. When computing C, the discontinuities present in elements such as ρ(1, 2) cancel out, but the spikes present in elements such as ρ(1, 1) do not. As a result, C has a singularity in the first derivative as would happen in a 2QPT and it is precisely ∂ Bx C the quantity which fulfills the scaling ansatz and not C itself. Therefore, in this case, a non-analyticity in the derivative of the concurrence/negativity (given that the concurrence/negativity are continuous analytic functions), is not a sufficient condition to signal a 2QPT as it was conjectured, [18]. The FSS shown in panel c) of Fig. 2 enables us to link the spike in C to a 1QPT despite that, at first sight, this behavior is the one expected for 2QPTs. We have then applied our analysis to the spin-1 XXZ chain with on-site anisotropy. This system has several 1QPTs and we have found 4 patterns for N along these transitions. These patterns are: a discontinuity, a smooth slope, a dip (with a singularity in the first derivative) and a constant N. The discontinuous behavior, which is the one expected for 1QPTs, is observed in the ferromagnetic/large D transition where no finite size effects are present, even for small systems. For the large D/Néel transition, when crossing close enough to the tricritical point, a smooth slope in N, instead of a discontinuity, appears due to finite size effects. Hence, one has to be careful when making a one-to-one connection between 1QPTs and discontinuous bipartite entanglement for finite systems, especially when a 2QPT is close. The FSS for 1QPTs can link this smooth slope to a discontinuity in the thermodynamic limit. For the Néel/Anti-Néel, a dip in N is observed signalling a singularity in the first derivative. This behavior has the same origin as in the spin-1/2 case and, therefore, it is the derivative of N the one which fulfils the FSS for 1QPTs. Finally, we mention the ferromagnetic ↑/ferromagnetic ↓ transition where bipartite entanglement is constant and zero along the transition. In this case the discontinuities of the elements of ρ A,B trivially cancel out giving a constant and zero N and, thus, bipartite entanglement does not signal the QPT.
In conclusion, finite size effects become relevant when studying a 1QPT which is close to a 2QPT. Our main result is to have shown that C and N scale in a 1QPT as described in Ref. [1]. We have observed that, close to a 1QPT, the discontinuities present in some elements of ρ A,B cancel out when computing the measures of bipartite entanglement but that the singularities present in the first derivative of some other elements do not. In this controversial example, where C has a spike at the critical point showing the apparent behavior of a 2QPT, the finite size scaling can unequivocally determine the order of the transition.