Exploring quantum phase transitions by the cross derivative of the ground state energy

In this work, the cross derivative of the Gibbs free energy, initially proposed for phase transitions in classical spin models [Phys. Rev. B 101, 165123 (2020)], is extended for quantum systems. We take the spin-1 XXZ chain with anisotropies as an example to demonstrate its effectiveness and convenience for the Gaussian-type quantum phase transitions therein. These higher-order transitions are very challenging to determine by conventional methods. From the cross derivative with respect to the two anisotropic strengths, a single valley structure is observed clearly in each system size. The finite-size extrapolation of the valley depth shows a perfect logarithmic divergence, signaling the onset of a phase transition. Meanwhile, the critical point and the critical exponent for the correlation length are obtained by a power-law fitting of the valley location in each size. The results are well consistent with the best estimations in the literature. Its application to other quantum systems with continuous phase transitions is also discussed briefly.


I. INTRODUCTION
Exploring novel phases of matter and phase transitions has always been one of the central topics in statistical and condensed matter physics. It has greatly enriched our understanding of matter phases and attracted much attention and effort since the discovery of the topological phases and phase transitions 1-5 beyond Landau's symmetry-breaking theory. Previously, we have proposed and demonstrated that the cross derivative of the Gibbs free energy is efficient and convenient for detecting various phase transitions in classical spin models 6 , no matter whether a transition is conventional or exotic with topological excitations. Its success inspires us to extend its applicability to quantum cases, wherein the driving forces of a phase transition are much more diversified. Accordinglly, the Gibbs free energy is reduced to the ground state energy, and the temperature fluctuation is replaced by the frustration effect between multiple competing interactions.
Usually, the introduction of competition into a quantum system raises complexity and difficulty in identifying phase transitions and the phase diagram. In particular, the transitions higher than 2nd-order are tough to determine precisely. One typical example is the spin-1 XXZ chain with the single-ion anisotropies, where the 3rd-and 5th-order Gaussian-type quantum phase transitions have been suggested for different parameters [7][8][9][10] . For these transitions, the conventional 2nd-order differential of the ground state energy doesn't work as usual, and the fidelity susceptibility is not necessarily applicable. 11 Because of its rich phase diagram, this model has been widely utilized to study effective one-dimensional spin-1 magnetic materials, as reviewed in Ref. [12].
It also serves as a testing ground for sophisticated numerical methods. Recently, the tunable single-ion anisotropic effects have been realized experimentally in spin-1 models with ultra-cold atoms 13 , and in the compound [Ni(HF 2 )(3-Clpyradine) 4 ]BF 4 with inelastic neutron scattering by pressure variance 14 . There are more details about the anisotropic effects and the phase diagram of this system in Refs. [7,8,[15][16][17][18][19]. Due to the inefficiency of the conventional differential tool, it is usually studied with unique tactics such as entanglement entropy, fidelity susceptibility, etc. [7][8][9]15,20,21 .
In this work, we use this model as an example to demonstrate that, by making use of competing interactions, the cross derivative captures the essential characteristics of the phase transitions. By tuning the uniaxial anisotropy strength D and the rhombic one E, we first investigate the 3rd-order Gaussian-type transitions in this model when J z = 1. In the (D, E)-plane, the critical point between the Haldane phase and the Large-D phase is determined at (0.9687, 0), and the Haldane-Large-E critical point locates at (−0.4862, 0.4862). Both are well consistent with previous predictions 7, 8,19 . Moreover, the critical exponent for the correlation length is obtained simultaneously. We also study the more difficult 5thorder Gaussian-type Haldane-Large-D transition when J z = 0.5. The critical point is estimated precisely at (0.6197, 0), which agrees well with best estimations in the literature 8,9,19 . Thus briefly, our cross derivative method provides a convenient, efficient, and universal tool to detect phase transitions not only in classical spin models 6 but also in quantum systems, whether a transition is conventional or exotic with higher order.
The rest of the paper is organized as follows. Sec. II introduces the model and the method employed in our arXiv:2210.08250v2 [cond-mat.str-el] 11 Apr 2023 work. The results are presented in Sec. III. Finally, Sec. IV gives a summary and discussion.

II. MODEL AND METHOD
In this work, we study the spin-1 XXZ chain with anisotropies, whose Hamiltonian reads as where S α i (α = x, y, z) are spin-1 operators on the i-th site, and L is the length of the spin chain. The strength of the Heisenberg exchange interaction in xy-plane is set to 1 for convenience. J z is the anisotropy strength in z-direction with respect to xy-plane. In this work, we only consider two special cases, i.e. J z = 1 and J z = 0.5. D and E are the strengths of the uniaxial and rhombic single-ion anisotropies, respectively.
To calculate the ground state energy, we employed the density matrix renormalization group (DMRG) method [22][23][24] , by using ladder scheme and encoding parity symmetry 25 with a periodic boundary condition. Regarding the accuracy with respect to the bond-dimension m, we follow the strategy of Ref. [25] to choose big enough m to make sure the truncation error is smaller than 10 −9 for each system size, thus to guarantee the obtained energy is precise enough for further differential computation. Practically, m also increases gradually with system size, i.e., m = 700 for L = 40, and m = 1200 for L = 80. The obtained result is also a good reflection of the computing precision. For illustration, we usually adopt an intermediate system size with L = 40. Afterward, a finite size extrapolation to the thermodynamic limit is performed. The quantity proposed to detect the quantum phase transitions, is the cross derivative of the ground state energy density ( = H /L) with respect to the anisotropic strengths D and E, i.e., ∂ 2 /∂D∂E. Mathematically, to get complete information of a 3-dimensional surface f (x, y, z) = 0 at a given point (x, y, z), one needs in principle not only the curvatures in x-and y-directions but also the twist, namely, where ∂/∂t ≡ α∂/∂x + β∂/∂y is the slope operator in an arbitrary direction. An equal weight is chosen for simplicity, i.e., α = β. As demonstrated below, this cross derivative contains the contributions from both orthogonal directions, and then is able to detect quantum phase transitions, especially those higher-order ones, for which the curvature in either principal direction is inadequate.
To obtain this quantity at any point in the (D, E)plane, we apply a normal central differential formula as where h is set to 10 −3 . The error induced by this formula is of order O(h 2 ). Since the Hamiltonian is symmetric between S x and S y , D,E and D,−E are equal, and thus the above cross derivative on the line E = 0 is zero. In this situation, we instead use a rotated one as where the parameters X = (D + E)/ √ 2, and Y = (D − E)/ √ 2. It is equivalent to rotating the coordinate frame clockwise by π/4 , as schematically shown in Fig. 1.

III. RESULTS
In the following two subsections, we will study the phase transitions in this model with J z = 1 by using the above cross derivative. A rough phase diagram is shown in Fig. 2 below to stress the focus, as labelled in different colors.
A. Jz = 1, 2nd-order phase transition When J z = 1 and E = 0, the Eq. (1) model reduces to the spin-1 Haldane chain with uniaxial anisotropy and has a 3rd-order Gaussian-type phase transition between the Haldane and the Large-D phases at D c1 ∼ 0.9684 7,8 . As shown in Fig. 2, in the (D, E)-plane, the y-Néel (x-Néel) phase lies above (below) this point (D c1 , 0), and between them, there is a 2nd-order quantum phase transition cross the point (D c1 , 0) along the vertical D = D c1 blue line.
Here, we first demonstrate the effectiveness of the cross derivative for this 2nd-order phase transition with the system size L = 40. Fig. 3 shows the results, where the transition point D = 0.893730 is utilized (see Fig.  4(a)). In Fig. 3(a), one can clearly see that the 1storder differential of the ground state energy ∂ /∂E is continuous, while the 2nd-order one ∂ 2 /∂E 2 is divergent at E = 0, with D fixed at D . This is the smoking evidence of a 2nd-order phase transition. In Fig. 3(b), both the normal and the rotated cross derivatives show clear divergence at the same point (D , E = 0), rendering identical information as in Fig. 3(a). So, it shows clearly that in the study of 2nd-order phase transitions, the cross derivative can work equally well with the conventional derivative method.
B. Jz = 1, 3rd-order Gaussian-type transitions In this subsection, we will focus on the 3rd-order Gaussian-type quantum phase transitions in this model with J z = 1 and further illustrate the versatility of the proposed novel function. As is well known, the 2nd-order differential of the ground state energy with respect to D, namely ∂ 2 /∂D 2 , is incapable of detecting the transition between the Haldane and Large-D phases, as shown in Fig. 7 of Sec. VI, and in Fig. 3(c) of Ref. 8 and Fig. 4(b) of Ref. 9. More explicitly, Fig. 7(b) has no power-law fitting property as the cross derivative does. Fig. 3(c) of Ref. 8 shows ∂ 2 /∂D 2 near the critical point, but neither any peak structure nor discontinuity shows up. In Fig. 4(b) of Ref. 9, although d 2 (L)/dD 2 exhibits a broad peak-like structure, its location is somehow invariant as L grows and is far away from the expected critical point. The peak height even decreases slightly instead of diverging, increasing L from 100 to 220. None of these phenomena matches with a phase transition scenario.
According to Refs. [7-9, and 15], when J z = 1, there are indeed three 3rd-order phase transitions from the Haldane phase to the Large-D phase at (D c1 , 0), and two Large-E phases at (D c2 , ∓D c2 ) respectively, as shown in Fig. 2. This also explains why the 2nd-order differential ∂ 2 /∂D 2 is inadequate to identify these transitions. Here, we calculate the cross derivative ∂ 2 /∂D∂E (or ∂ 2 /∂X∂Y ) to investigate these transitions. Since the Hamiltonian is invariant for opposite Es, the two Large-E phases are symmetric about the line E = 0, as well as the two Haldane-Large-E critical points. So, we pick up one of these two. Figure 4 shows the transition between the Haldane phase and the Large-D phase on the line E = 0, where the rotated cross derivative is adopted, as explained earlier. The result of ∂ 2 /∂X∂Y for L = 40 is presented as an illustration in Fig. 4(a), where a single valley shows up clearly. In the inset, a quadratic fitting is performed around the valley minimum as, and the valley depth is obtained as H p (L = 40) = −1.039945, with the location at D p (L = 40) = 0.893730. We repeat this process for different system sizes from 24 to 80, and the valley depths for each size H p (L) are drawn in Fig. 4(b), which matches a logarithmic fitting perfectly as with a = −0.92(2), b = 48(2), and c = 4.9(2). The logarithmic divergence of the valley depth with increasing the system size signals a phase transition therein. It then verifies the validity of the cross derivative for a phase transition. Furthermore, the position of the valley minimum for each size D p (L) is collected in Fig. 4(c). Then, a power-law fitting 26 is performed as By the finite-size extrapolation, the critical point is determined at D c1 = 0.9687(4), and the critical exponent for the correlation length is ν = 0.817 (5). The estimated critical point is consistent with the most precise predictions to date 7,8,20 , as listed in Tab. I below. The critical   (5) stiffness(QMC) 17 0.96845 (8) entropy(DMRG) 8 0.9684713(1) level spectroscopy(DMRG) 7 0.9685 (2) tangential finite size scaling 20 0.9687 (4) cross derivative(this work) ing the valley depth H p (L) for different sizes from 24 to 64, we can observe a clear logarithmic divergence, as fitted in Fig. 5(b). Again, this logarithmic divergence of the valley depth with the increase of the system size indicates a phase transition. At the same time, with the valley position of each size D p (L) shown in Fig. 5(c), we obtain the critical point at (D c2 , −D c2 ) with D c2 = −0.4862 (4) and the critical exponent for the correlation length as ν = 0.886 (7), by the finite size extrapolation. As expected, the transition point agrees well with the previous prediction 7 . An interesting and worth mentioning thing is that the logarithmic fitting function here is precisely half of that in Fig. 4(b). The difference is consistent with the definitions, as sketched in Fig. 1. This fact that these two logarithmic fitting functions coincide well with each other verifies the two separate calculations are consistent and indicates the two transitions are the same type.
So far, we have shown that the proposed cross derivative of the Gibbs free energy can detect and locate the 3rd-order quantum Gaussian-type transitions in this model. Furthermore, convenience, precision, and efficiency are explicitly displayed. Combining with Figs. 3 and 7, we should note that the cross derivative contains contributions from both orthogonal directions, where the divergence behavior comes from the nature of the 2ndorder transition in E-direction, as expressed in Eq. (4). This may explain the reason why the cross derivative works well.
C. Jz = 0.5, 5th-order Gaussian-type transition In this model with J z = 0.5, a 5th-order Gaussian-type quantum phase transition between the Haldane phase and the Large-D phase has also been reported along E = 0 at D c ≈ 0.63 [8][9][10]19 . However, this transition and the critical point are even more challenging to detect and locate. Recent researches have adopted the entanglement entropy, the fidelity susceptibility, or other complex quantities.
We follow the same logic for this case, and compute the rotated cross derivative ∂ 2 /∂X∂Y by fixing E = 0. The results are presented in Fig. 6, and the valley lies at D p = 0.3500 with depth H p = −2.862 for L = 40. A logarithmic fitting is performed from the collected data of H p (L) with L ranging from 24 to 72. They match each other perfectly, and the logarithmic divergence denotes the phase transition. By collecting D p (L) with different sizes, a finite size extrapolation is carried out to obtain the critical point at D c = 0.6197(3) and the critical exponent for the correlation length as ν = 1.138 (1) simultaneously. Both are pretty close to the estimations of D c = 0.63 and ν = 1.51 in Ref. [9], and D c = 0.635 in Refs. [8 and 19].
Once again, the validity and efficiency of the cross derivative are illustrated for this difficult 5th-order Gaussian-type quantum phase transition.

IV. SUMMARY AND DISCUSSION
In brief, we extend the scope of the cross derivative of the Gibbs free energy, initially proposed for phase transitions in classical spin models 6 , to the study of quantum cases. Its validity and efficiency have been demonstrated by the typical and challenging higher-order Gaussian-type phase transitions in spin-1 XXZ chain with anisotropies.
When J z = 1, the 3rd-order Gaussian-type phase transition is precisely located at (0.9687, 0) for the Haldane-Large-D transition and at (-0.4862, 0.4862) between the Haldane phase and the Large-E phase. The obtained critical points agree well with the most accurate estimations to date in the literature. As for J z = 0.5, the 5th-order Gaussian-type transition is determined at (0.6194, 0), also consistent with the previous predictions.
In both transitions, the critical exponent for the correlation length is a little smaller than the predictions in the literature. We should note that, for higher-order continuous phase transition, the critical exponent for the correlation length is more difficult to determine precisely than the location of the critical point. When approaching the critical point, the correlation length grows rapidly and becomes much larger than the size used in this work. To accurately estimate the critical exponent for the correlation length, one may need larger system sizes for the finite size scaling to eliminate the small size effect, as mentioned in Ref. [8], wherein the biggest system size utilized is 10 4 . To further improve the accuracy of the estimated critical exponent with a larger system size, one may try other numerical algorithms, like the recently developed variational corner transfer matrix renormalization group method 28 or the (infinite) time-evolving block decimation method [29][30][31] to even deal with an infinite onedimensional system directly. The other minor possibility is the equal weight (i.e., α = β) adopted in Eq. (2) for simplicity, which may not be optimal and requires further investigations. These will be presented in subsequent studies. Also, the order of the quantum critical point can be altered by changing the direction in the phase diagram across the critical point, which may also affect the critical exponent.
Given its simplicity and convenience, the cross derivative is efficient and universal to investigate the phase transitions in quantum spin systems, whether it is the conventional 2nd-order one or the complicated Gaussiantype one. Moreover, the predictions will be more accurate if the system size and the precision of the ground state energy can be improved further. The method is also readily applied to other complex systems, like the argued 3rd-order phase transition in Bose-Einstein condensation [32][33][34] . Requiring only the precise ground state energy, this quantity is much easier to deal with than wave function, correlation functions or order parameters used in other methods. (a) Illustration of the 2nd order differential ∂ 2 /∂D 2 along D direction with L = 40, and a quadratic fitting near the peak to give the peak position; (b) Power fitting of the peak position of the differential ∂ 2 /∂D 2 , and the fitting still fails even we neglect two values for small size systems. For comparison, the rotated cross derivative, i.e. Fig. 4(c) is also included in blue.
According to the statistical physics, an n-th order phase transition is detected from the divergence of the nth order differential of the Gibbs free energy, while whose (n − 1)-th order differential is continuous. So naturally, a 2nd-order differential is not able to identify a higherorder phase transition. Here, we show below in Fig. 7 the failure of the 2nd order differential of the Gibbs free energy with respect to D, to detect the 3rd order Gaussian transition in this model along D direction when J z = 1. For comparison, the rotated cross derivative (Fig. 4(c)) is also included.