Fidelity susceptibility made simple: A unified quantum Monte Carlo approach

The fidelity susceptibility is a general purpose probe of phase transitions. With its origin in quantum information and in the differential geometry perspective of quantum states, the fidelity susceptibility can indicate the presence of a phase transition without prior knowledge of the local order parameter, as well as reveal the universal properties of a critical point. The wide applicability of the fidelity susceptibility to quantum many-body systems is, however, hindered by the limited computational tools to evaluate it. We present a generic, efficient, and elegant approach to compute the fidelity susceptibility of correlated fermions, bosons, and quantum spin systems in a broad range of quantum Monte Carlo methods. It can be applied both to the ground-state and non-zero temperature cases. The Monte Carlo estimator has a simple yet universal form, which can be efficiently evaluated in simulations. We demonstrate the power of this approach with applications to the Bose-Hubbard model, the spin-$1/2$ XXZ model, and use it to examine the hypothetical intermediate spin-liquid phase in the Hubbard model on the honeycomb lattice.


I. INTRODUCTION
Phase transitions highlight the beauty of universality, despite the great diversity of nature. For example, one finds a unified description for systems ranging from ultracold bosons [1,2] to magnetic insulators [3][4][5] on the verge of a phase transition. Phase transitions originate from the competition between different tendencies when a macroscopic system tries to organize itself. Thermal fluctuations can drive classical phase transitions at nonzero temperatures, while quantum phase transitions can occur even at zero temperature because of the competition between noncommuting terms in the quantum-mechanical Hamiltonian [6,7]. At the phase transition point, physical observables often exhibit singular behavior. In this respect, phase transitions are the most dramatic manifestation of the laws of statistical and quantum mechanics.
Traditional descriptions of phase transitions are based on low-energy effective theories of local order parameters, which have had enormous success in explaining various phase transitions of superfluids, superconductors [8], and quantum magnets [9,10]. However, in recent years, exceptions to this Ginzburg-Landau-Wilson paradigm have emerged [11]. In particular, topological phase transitions [12][13][14] do not have a local order parameter on either side of the phase transition. Therefore, new theoretical tools are needed to search for and characterize these new quantum phases and phase transitions. Many concepts in quantum information science [15], such as quantum fidelity and quantum entanglement, have proven to be useful [16,17]. Having a point of view that is totally different from the traditional condensed matter approach, they do not assume the presence of a local order parameter and thus offer new perspectives of the phase transitions and their universalities.
Specifically, we consider the following one-parameter family of Hamiltonians with a driving parameter λ: As λ changes, the system may go through one or several phase transition(s) because of the competition betweenĤ 0 andĤ 1 . The quantum fidelity measures the distance on the manifold of λ, which is defined as the overlap between the ground-state wave functions at two different values of the driving parameter, Since in general the quantum fidelity vanishes exponentially with the system size for a many-body system, it is more convenient to study the change of its logarithm with respect to the driving parameter, called the fidelity susceptibility [21]: : ð3Þ The first-order derivative vanishes because F is at its maximum when ϵ ¼ 0. In general, the fidelity susceptibility is an extensive quantity away from the critical point, but it exhibits a maximum or even diverges at the critical point, thus indicating a quantum phase transition [22,23]. Similar to conventional thermodynamic quantities, it also follows a scaling law close to the critical point [22][23][24][25], which can be used to extract universal information about the phase transition. An important feature of the fidelity susceptibility is that it can reveal a phase transition without prior knowledge of the local order parameter. This makes it suitable for detection of topological phase transitions [26][27][28][29] and Berezinsky-Kosterlitz-Thouless-type transitions [30][31][32] as well as for tackling challenging cases where an in-depth understanding of the underlying physics is still lacking [33,34]. Interestingly, the fidelity susceptibility may also be accessible to experiments [35][36][37]. Despite its appealing features, the difficulty in calculating the fidelity susceptibility has hindered its use in numerical simulations. Many previous studies were thus limited to the cases where the ground-state wave function overlap could be calculated from the analytical solution, exact diagonalization, or density-matrix renormalizationgroup (DMRG) methods [16].
There are several equivalent formulations of the fidelity susceptibility Eq. (3), which reveal different aspects of the quantity. From a computational point of view, they offer direct ways to calculate the fidelity susceptibility without the need to perform numerical derivatives of the fidelity as in Eq. (2). (a) Expanding jΨ 0 ðλ þ ϵÞi for small ϵ, one can cast the definition Eq. (3) into an explicit form [22,38]: The above form does not assume properly normalized wave functions jΨ 0 i. Equation (4) reveals the geometric content of the fidelity susceptibility [22,38], since this expression is the real part of the quantum geometric tensor [39]. (b) Alternatively, one can calculate the first-order perturbation for jΨ 0 ðλ þ ϵÞi and get Compared to Eq. (4), Eq. (5) does not contain derivatives but involves all eigenstates and the full spectrum. It explicitly shows that χ F ≥ 0 and suggests the divergence of χ F when the energy gap of the system closes. (c) Reference [21] views Eq. (5) as the zero-frequency component of a spectral representation; thus, a Fourier transform is performed to obtain an alternative expression, whereĤ 1 ðτÞ ¼ eĤ τĤ 1 e −Ĥτ . Equation (6) has the form of a linear-response formula and is computationally more friendly than Eq. (4) or Eq. (5). References [24,25] generalize it to nonzero temperature by replacing the integration limit with β=2: where hÁ Á Ái denotes the thermal average at inverse temperature β. Besides reducing to Eq. (6) as β → ∞, Eq. (7) is nevertheless a well-defined quantity at nonzero temperatures. It bounds the divergence of an alternative "mixed state" fidelity susceptibility [40,41], which is based on the Uhlmann fidelity [42,43], and both quantities follow the same scaling law close to a quantum critical point [24,25]. In general, the evaluation of Eq. (7) is still a formidable computational task, which requires ad hoc implementation depending on the details of the Hamiltonian. For example, the fidelity susceptibility for twodimensional quantum spin systems is calculated in Refs. [24,25] using a quantum Monte Carlo (QMC) method, while for a one-dimensional quantum spin system, Ref. [44] computed it using the transfermatrix DMRG method. The nontrivial implementation of these specific approaches and the overhead in the calculation still limits the wide applicability of the fidelity susceptibility approach to a broad range of quantum many-body systems. In this paper, we present a simple yet generic approach to compute the fidelity susceptibility in a large variety of modern quantum Monte Carlo methods, including the continuoustime worldline [45][46][47][48] and stochastic series expansion (SSE) [49] methods for bosons and quantum spins, and the diagrammatic determinantal methods for quantum impurity [50][51][52][53] and fermion lattice models [54][55][56].
In all cases, the Monte Carlo estimator is generic and the implementations are straightforward. As long as the quantum Monte Carlo simulation is feasible (not hindered by the sign problem), the fidelity susceptibility can be easily calculated. Our finding can boost the investigation of quantum phase transitions from a quantum information perspective and becomes especially advantageous for the exploration of exotic phases beyond the Ginzburg-Landau-Wilson paradigm.
The organization of the paper is as follows. In Sec. II, we present our estimator for the fidelity susceptibility and discuss its implementations in various quantum Monte Carlo methods. Section III presents derivations of the estimator. In Sec. IV, we demonstrate the power of the fidelity susceptibility approach with applications to various models, including correlated bosons, fermions, and quantum spins, using a variety of quantum Monte Carlo methods. Section V discusses the relation between the zerotemperature and nonzero-temperature estimators for the fidelity susceptibility and compares them to the previous approaches [24,25]. We conclude with future prospects in Sec. VI.

A. Universal covariance estimator
Many modern QMC methods [45][46][47][48][49][50][51][52][53][54][55][56] share a unified conceptual framework, namely that the partition function is calculated as a perturbative series expansion for the λĤ 1 term, where the second summation runs over all the Monte Carlo configurations of a given expansion order k. The detailed meaning of the configuration depends on the specific QMC algorithm and will be explained in the next subsection. Figure 1(a) depicts a generic configuration, where the k objects residing on the periodic imaginary-time axis represent the vertices λĤ 1 in the expansion, with a Monte Carlo weight λ k wðC k Þ for this configuration. QMC simulations [45][46][47][48][49][50][51][52][53][54][55][56] sample the summation over k and C k on an equal footing. Specific algorithms differ by the detailed form of wðC k Þ and by the sampling schemes. Nevertheless, these QMC methods share a unified framework provided by Eq. (8), which is the only requirement for the estimator of the fidelity susceptibility Eq. (7) to possess an appealing universal form in nonzero-temperature QMC simulations: where k L and k R are the number of vertices residing in the range ½β=2; βÞ and ½0; β=2Þ of the imaginary-time axis, respectively, shown in Fig. 1(a). The Monte Carlo average of an observable is defined as hOi ¼ 1 where OðC k Þ denotes the value measured for the configuration C k . In practice, because of the periodic boundary condition on the imaginary-time axis, the division of the time axis to halves may be done at an arbitrary location. Moreover, it is even possible to perform multiple measurements on the same configuration by generating several random divisions.
QMC methods [45][46][47][48][49][50][51][52][53][54][55][56] can also be utilized at zero temperature, where the unnormalized ground-state wave function is obtained from an imaginary-time projection: Here, β is a projection parameter and the trial wave funciton jΨ T i shall not be orthogonal to the true ground state. A similar framework as Eq. (8) applies, except that one now samples from the overlap hΨ T je −βĤ jΨ T i instead of the partition function. In the projection scheme, the fidelity susceptibility has the estimator where k L and k R are the number of vertices for the bra and ket states, which reside in the range ½β=2; βÞ and ½0; β=2Þ of the imaginary-time axis respectively, shown in Fig. 1(b).
FIG. 1. Measurement of the fidelity susceptibility in the (a) nonzero-temperature formalism and in the (b) ground-state projection scheme. Each red object represents a term in the driving Hamiltonian λĤ 1 , denoted as a vertex. To measure the fidelity susceptibility Eqs. (9) and (11), we divide the imaginarytime axis into two halves and count the number of vertices k L and k R , respectively. The nonzero-temperature formalism allows an arbitrary division because of the periodic boundary condition in the imaginary-time axis, while in the ground-state projection scheme the division has to be at β=2.
Since the fidelity susceptibility is non-negative, the covariance formula Eq. (11) reveals positive correlation of k L and k R in a Monte Carlo simulation. Equations (9) and (11) are the central results of this paper. As is obvious from the discussions in this section, neither the details of the Hamiltonian nor the statistics of the system need to be specified. These estimators are thus general and can be readily implemented in a variety of QMC methods for correlated fermionic, bosonic, or quantum spin systems [45][46][47][48][49][50][51][52][53][54][55][56].

Continuous-time worldline and diagrammatic determinantal approaches
Continuous-time worldline methods [45][46][47][48] are widely used to simulate boson and quantum spin systems, while the diagrammatic determinantal approaches are the stateof-the-art methods for solving quantum impurity [50][51][52][53] and fermion lattice models [54][55][56]. A common feature of these methods is to split the Hamiltonian in the form of Eq. (1) and perform a time-dependent expansion in λĤ 1 , which obviously fits in the general framework of Eq. (8).
In the continuous-time worldline approach [45][46][47][48], thê H 1 term corresponds to hoppings of bosons or spin flips, depicted as kinks of the worldlines in Fig. 2(a). In continuous-time diagrammatic determinantal approaches [50][51][52][53][54][55][56], λĤ 1 contains the fermion interactions, drawn as interaction vertices in Fig. 2(b). Equation (12) has the form of a grand canonical partition function for a classical gas, where λ plays the role of fugacity and k is the number of certain classical objects (kinks or vertices) residing on the imaginary-time axis. Typical updates of continuous-time diagrammatic determinantal approaches [50][51][52][53][54][55][56] consist of randomly inserting or removing vertices, which are identical to the updates of the grand canonical Monte Carlo method for molecular simulations [57,58]. For bosons and quantum spins there are more effective nonlocal updates, such as the worm and directed loop updates [45][46][47][48][49]. In any case, the Monte Carlo estimators Eqs. (9) and (11) are independent to the detailed sampling procedures. It suffices to count k L and k R of Monte Carlo configurations to calculate the fidelity susceptibility. Examples are presented in Secs. IVA and IV C.

Stochastic series expansion
SSE is based on a Taylor expansion of the partition function [49], which may seem to be different from the framework of Eq. (8). However, as is shown in Ref. [59], one can formally treat the SSE as the time-dependent expansion Eq. (12) with respect to the full HamiltonianĤ ¼Ĥ 0 þ λĤ 1 .
In the implementation of SSE, one truncates the sum to a large number M and pads M − n identity operators in the square bracket of Eq. (13). SSE then samples operators in the fixed-length operator string. To map to a Monte Carlo configuration in the continuous-time formalism, one can assign an imaginary time to each operator, as shown in the bottom of Fig. 3. As long as the mapping keeps the relative order in the original operator string, the Monte Carlo weight remains unchanged [60]. In particular, the configuration is sampled with a weight proportional to λ k if there are k of λĤ 1 operators in the operator string. In this way, although the sampling of SSE is carried out differently from Eq. (12), the general framework of Eq. (8) still applies. The fidelity susceptibility is then measured easily by counting the numbers k L and k R of operators associated with λĤ 1 in the two halves of the imaginary-time axis after the mapping.
From Fig. 3 it is clear that even though one performs an equal bipartition in the imaginary-time axis, the corresponding location of division is not always in the center of the operator string. In fact, it is easier to directly sample the location of division in the operator string, as shown in the upper part of Fig. 3. A division at the lth position (l ¼ 0; 1; …; M) means that there are l slots being mapped to one half of the imaginary-time axis and M − l slots to the other half. Therefore, the division l itself follows a binomial distribution pðlÞ ¼ 1 2 M ð M l Þ that can be sampled directly. In this way, the fidelity susceptibility can be efficiently calculated in SSE similar to the continuous-time QMC approaches discussed in Sec. II B 1. As M → ∞, the binomial distribution approaches a delta function peaked at the center of the operator string. Only in this limit the position in the operator string can be directly interpreted as the imaginary time and a bipartition of the operator string in the center will yield the correct result for the fidelity susceptibility.

III. DERIVATIONS
In this section, we derive the estimators for the fidelity susceptibility at nonzero temperature Eq. (9) and for ground-state projector formalism Eq. (11). Readers may skip this section and continue reading with the following sections.

A. Ground state
Using the diagrammatic expansion for the projection operator, the unnormalized ground-state wave function Eq. (10) has the following form: Substituting it into Eq. (4) and noticing that the partial derivatives affect only the prefactor λ k , one obtains Eq. (11). The estimator also holds for the continuous-time auxiliary field expansion methods (CT-AUX [52] and LCT-AUX [56]), because one can cast the ground-state wave function to a similar form as Eq. (14), assuming the shift parameter used in these methods [54] to be proportional to λ.

B. Nonzero temperature
We present two derivations of the nonzero temperature estimator Eq. (9).

Derivation based on the definition of nonzero-temperature fidelity
The nonzero-temperature estimator Eq. (9) can be obtained directly from the definition of the nonzerotemperature fidelity [61]: This is a nonzero-temperature generalization of Eq. (2) and leads to Eq. (7) by using the definition of the fidelity susceptibility Eq. (3) [44]. We expand the traces of the density matrices around the partition function Z ¼ Trðe −βĤðλÞ Þ to Oðϵ 2 Þ, where the notation∂ λ indicates that the partial derivative acts only on operators in the imaginary-time interval 0 ≤ τ < β=2. Substituting the above two expansions into Eq. (15) and keeping terms up to Oðϵ 2 Þ, one obtains To obtain the second line, we use the partition function expansion Eq. (8) and the fact that the partial derivatives act only on the prefactor λ k ¼ λ k L þk R . For the third line, we use k ¼ k L þ k R and hk L i ¼ hk R i ¼ hki=2. This derivation is abstract and is independent of the details of a QMC scheme.
Carrying out a similar procedure starting from Eq. (2), one can also prove the ground-state estimator Eq. (11).
FIG. 3. Division of the operator string in SSE to measure the fidelity susceptibility. The slots represent the fixed-length operator string where empty slots hold identity operators, the red circles (blue squares) correspond to the operators inĤ 1 (Ĥ 0 ). These operators can be mapped to a continuous-time configuration indicated by the arrows. A division can then be made on the imaginary-time axis, for example, at β=2. An equivalent approach without explicit mapping to continuous time is to divide the operator string at the location indicated by the vertical dashed line, where the integer l is drawn from a binomial distribution. For the estimators Eqs. (9) and (11), one counts the number of red circles (operators in λĤ 1 ) in both sides for k L and k R . In this example, M ¼ 12; n ¼ 6; l ¼ 7, and k L ¼ k R ¼ 2.
2. Derivation based on the imaginary-time correlator Eq. (7) This derivation starts from the definition of fidelity susceptibility based on the imaginary-time correlator Eq. (7). We utilize its connection to the Monte Carlo weight that appears in Eq. (12) to derive the nonzerotemperature estimator Eq. (9). First of all, the second term in the square bracket of Eq. (7) can be measured directly from the average expansion order [49,50,53,54]: Integrating over the imaginary-time and using hk L i ¼ We then consider the QMC estimator for the first term of Eq. (7), where T is the time-ordering operator and τ i and τ j are the imaginary times of two vertices in the Monte Carlo configuration. Integrating both sides of Eq. (19) and using the fact that Gðτ 1 − τ 2 Þ depends only on jτ 1 − τ 2 j, one has where kðΛÞ is the number of vertices in the range of 0 ≤ τ < Λ. For example, kðβÞ ¼ k and kðβ=2Þ ¼ k R . When choosing Λ ¼ β and using GðτÞ ¼ Gðβ − τÞ, Eq. (20) becomes When setting Λ ¼ β=2, Eq. (20) reads Together with Eq. (21), it leads to In combination with Eqs. (7) and (18), we obtain Eq. (9).

IV. APPLICATIONS
We first demonstrate the power of the new approach by identifying quantum and thermal phase transitions in the Bose-Hubbard model and in the spin-1=2 XXZ model. Then, we use the fidelity susceptibility to address the presence of the intermediate quantum spin-liquid state in the Hubbard model on the honeycomb lattice. In all cases, this requires only minimal modifications to existing codes. We purposely chose a variety of QMC methods in the following to demonstrate the wide applicability of the covariance estimators. Because of the flexibility of the nonzero-temperature estimator, we use Eq. (9). In Sec. VA, we compare it to the zero-temperature scheme.

A. Quantum phase transition in the Bose-Hubbard model
First, we use the fidelity susceptibility to probe the quantum phase transition in the Bose-Hubbard model, where U is the on-site interaction and μ is the chemical potential. The driving parameter λ has the physical meaning of a tunneling amplitude. The Bose-Hubbard model has a well-known quantum phase transition between the Mott insulating state and the superfluid state as λ=U increases [1,2]. In particular, for integer fillings, the system has an emergent Lorentz invariance at the critical point and the dynamical critical exponent is z ¼ 1 [7]. The fidelity susceptibility has previously been calculated using DMRG methods for the one-dimensional Bose-Hubbard model [62][63][64]. We now calculate the fidelity susceptibility on a square lattice with N ¼ L 2 sites at unit filling by tuning μ. In accordance with the dynamical critical exponent z ¼ 1, we scale the inverse temperature proportionally to the system length βU ¼ 4L. The simulation employs the directed worm algorithm [46,65,66]. We utilize Eq. (9) to sample the fidelity susceptibility by counting the number of kinks in the worldline configuration, as illustrated in Fig. 2(a). Figure 4 shows that, as the system size increases, the peak in the fidelity susceptibility (as a function of the driving parameter λ) is becoming more pronounced around the previously determined critical point ðλ=UÞ c ¼ 0.05974ð3Þ [67]. The inset of Fig. 4 shows the scaled fidelity susceptibility χ F L −2=ν versus ½ðλ=UÞ − ðλ=UÞ c L 1=ν according to the scaling law of Refs. [22,24,25]. The data collapses well under the critical exponent ν ¼ 0.6715 of the 3D XY universality class [67]. The ability to calculate the fidelity susceptibility using the state-of-the-art directed worm algorithm [46,65,66] will greatly advance the study of quantum phase transitions of ultracold bosons. It is worth pointing out that the fidelity susceptibility is related to the quantity (kineticenergy correlator) previously calculated in the study of Higgs mode in a two-dimensional superfluid [68].

B. Thermal phase transition in the XXZ model
Next, we consider the spin-1=2 antiferromagnetic XXZ model on a square lattice with N ¼ L 2 sites, where the driving parameter λ plays the role of the coupling strength in the XY plane. When λ dominates, the Hamiltonian favors Néel order in the XY plane, while if J z dominates, the system has an antiferromagnetic Ising ground state. The Heisenberg point λ ¼ J z is a quantum critical point, which separates the XY order and the Ising order. This quantum critical point can be easily located from the peak of the fidelity susceptibility (not shown). Our approach makes it possible to obtain the fidelity susceptibility in much larger systems compared to the previous exact diagonalization study [69], and thus can enable a more accurate scaling analysis. At nonzero temperature, thermal fluctuations will destroy the antiferromagnetic Ising phase at a second-order phase transition. Since one can cross the phase boundary by changing either λ or the temperature T, we see that the fidelity susceptibility can also indicate thermal phase transitions. As a demonstration, we fix λ ¼ 1, J z ¼ 1.5 and scan the temperature T to drive a phase transition from the low-temperature antiferromagnetic Ising phase to the high-temperature disordered phase. Figure 5 shows the fidelity susceptibility calculated using Eq. (9) via the SSE method [49,71]. The peak in the fidelity susceptibility correctly singles out the previously determined critical temperature ðT=λÞ c ≈ 0.75 [70].

C. Intermediate phase in the Hubbard model on the honeycomb lattice
Finally, we apply the fidelity susceptibility estimator to a more challenging and controversial example-the Hubbard model on the honeycomb lattice, where λ has the meaning of on-site Hubbard interaction strength. The simulation employs the recently developed efficient continuous-time QMC method for lattice fermions (LCT-INT) [56,72]. We consider lattices with N ¼ 2L 2 sites, with L ¼ 6; 9; 12, and scale the inverse temperature βt ¼ L.
The ground-state phase diagram of the Hubbard model on the honeycomb lattice [74] has been controversial. It was suggested to possess an intermediate nonmagnetic spin-liquid phase for λ=t ∈ ½3.5; 4.3 [75]. However, more recent QMC studies on larger systems [76] and with improved observables [77,78] suggest a single continuous phase transition at λ=t ≈ 3.8 belonging to the Gross-Neveu universality class [79]. Other less unbiased methods, such as quantum cluster approaches, give conflicting results on the presence of the intermediate phase [80][81][82][83][84][85], depending on implementation details.
The fidelity susceptibility offers a new perspective on the debate about the phase diagram. In the scenario with an intermediate phase, there shall be two features in χ F when λ=t approaches the two phase boundaries. This consideration is independent of the presence of a local-orderparameter description of the possible intermediate phase. Figure 6 shows the fidelity susceptibility per site for various system sizes obtained using Eq. (9). The calculation of the fidelity susceptibility for fermions is more challenging than that for bosons (Sec. IVA) and quantum spins (Sec. IV B), because of the unsatisfactory local Monte Carlo updates and the cubic (instead of linear) scaling with respect to the system size N. Given the limited system sizes, in Fig. 6 we can identify only a single broad peak for small systems (and high temperature). The peak becomes sharper and shifts towards smaller interaction strength as the system size increases; however, the examined system sizes are insufficient for a reliable determination of the phase diagram. Further algorithmic improvements (nonlocal updates and better scaling) on the fermionic QMC methods are needed to access accurate fidelity susceptibility data for larger systems and lower temperatures to better locate phase transitions without knowing the local order parameter. Moreover, by using the histogram reweighting [86,87] or the quantum Wang-Landau approach [88], it is possible to obtain the fidelity susceptibility in a continuous range of λ. Such data may be used to precisely determine the critical point and even the critical exponent of the quantum phase transition.

V. DISCUSSION
To help the reader gain a better understanding of the estimators Eqs. (9) and (11), we first discuss their relationship then compare them with the previous approach adopted in SSE calculations [24,25]. Finally, we compare the fidelity susceptibility approach with other generic approaches for detecting phase transitions.

A. Relation of the ground-state and nonzerotemperature estimators
The factor of 2 difference in Eqs. (9) and (11) is due to the different boundary conditions of the imaginary-time axis in the ground-state projection and nonzero-temperature QMC formalisms; see Fig. 1. We use a four-site Hubbard model [Eq. (26)] as an illustrative example. Consider the integrand of Eq. (7), the correlator GðτÞ ¼ hĤ 1 ðτÞĤ 1 i is related to the distribution of the vertices on the imaginary-time axis. For a given configuration, the probability of finding two vertices with a time difference τ is proportional to λ 2 GðτÞ. If we equally divide the imaginary-time axis into two halves and impose the additional constraint that the two vertices reside in different halves (denoted as a separable vertex pair), the joint probability changes to λ 2 GðτÞ minfτ; β − τg. Figure 7(a) shows the histogram of separable vertex pairs accumulated in the imaginary time, which indeed agrees with the exact curve. Summing up the histogram gives the total number of separable vertex pairs, which equals the following integration: Since GðτÞ is symmetric around τ ¼ β=2 in the nonzerotemperature simulation, the two terms of Eq. (27) are equal. Thus, Eq. (27) reduces to Eq. (23). Furthermore, Fig. 7(b) shows the correlator GðτÞ sampled by accumulating the histograms of distances between vertices [51,53] together with the exact results (solid black line). The correlation between vertices decays rapidly with imaginary-time distance and approaches the uncorrelated value hĤ 1 i 2 (dashed blue line). However, in the zero-temperature limit, the correlator GðτÞ decays monotonically with τ and two vertices will decorrelate for τ ≥ β=2, where β → ∞ in the projection scheme. Therefore, the second term of Eq. (27) reduces to ðλ 2 β 2 hĤ 1 i 2 Þ=8 ¼ hki 2 =8 and cancels half of the second term in the estimator Eq. (11), resolving the apparent difference by the factor of 2. In practical calculations, it is, however, crucial to adopt the correct formula to obtain consistent results, as illustrated in Fig. 8. The fidelity susceptibility calculated using Eq. (9) in a nonzerotemperature LCT-INT [56] simulation agrees perfectly with exact diagonalization results. The blue square shows the value obtained using Eq. (11) in a projector LCT-INT calculation [89], which correctly reproduces the exact value of the ground-state fidelity susceptibility. Figure 7(b) also reveals the difficulty of computing the fidelity susceptibility. If the decay of GðτÞ is faster than 1=τ, the integrand of Eq. (7) has vanishing contributions at large τ. However, as the two terms in Eq. (7) are sampled independently in the actual QMC simulations, uncorrelated vertices at large imaginary-time distance will cause noises in the fidelity susceptibility signal. For the applications in Sec. IV, we thus perform the calculations at nonzero temperature as it provides a natural cutoff.

B. Comparison to previous approaches
The present approach to sample the fidelity susceptibility is more generic and efficient than those developed in Refs. [24,25] specifically for the SSE method. It is nevertheless instructive to compare them in detail. The key difference lies in the sampling of the first term of Eq. (7). References [24,25] employ the SSE estimator [90,91] where GðnÞ is the number of occurrences of two operators fromĤ 1 that are separated by n positions in the fixedlength operator string (n ¼ 0 if they are next to each other).
Multiplying both sides with maxfτ; β − τg and integrating over the imaginary time, one finds where the weight function WðnÞ is written in terms of the regularized incomplete beta function I x ða; bÞ [92], References [24,25] explicitly go through kðk − 1Þ=2 pairs of vertices to accumulate hGðnÞi and multiply it with the weight function WðnÞ. In Eq. (27), however, the multiplication by the imaginary time τ is taken into account implicitly by the sampling procedure (which requires separable vertices). Besides being more generic, our approach reduces the computational cost from Oðk 2 Þ to OðkÞ, which is crucial for the simulation of bosonic and quantum spin systems. In this sense, the specification of our general result Eq. (9) for the SSE method can be regarded as an improved estimator of Eq. (28), which by itself already improves the approach of Refs. [24,25] in several aspects [93]. The improved estimator Eq. (9) not only unifies the SSE approach in a broader context of continuous-time diagrammatic QMC methods, it also gives better statistics with less computational cost compared to Eq. (28). Figure 9 shows the weight function WðnÞ for various truncation lengths. As M increases, it approaches two straight lines, and a division in the center of the operator string would yield increasingly accurate results for the fidelity susceptibility, consistent with the discussion in Sec. II B 2 concerning the large M limit.

C. Relationship to other quantities
The fidelity susceptibility is related to the second-order derivative of the free energy A ¼ −ð1=βÞ ln Z [25,94]. FIG. 8. QMC results for the fidelity susceptibility compared with exact results (solid line). The nonzero-temperature QMC data (red dots) are obtained from Eq. (9), while the ground-state data (blue square at β −1 ¼ 0) are obtained from Eq. (11) in a projector LCT-INT calculation [89]. The system is the same as in Fig. 7.
Because of the Hellmann-Feynman theorem [95,96], hĤ 1 i equals to the first-order derivative of the free energy with respect to λ. A further derivative following the Kubo formula gives The third equality follows from Eqs. (17) and (21) [97]. The quantity resembles the widely used SSE estimator for the specific heat [49,98], but can be used to probe quantum phase transitions [25]. At zero temperature ð∂hĤ 1 iÞ=∂λ ¼ −ð1=λÞð∂hĤ 0 iÞ=∂λ, and the latter quantity is computed using numerical differentiation of the kinetic energy, so as to address the quantum phase transition in the Hubbard model on the honeycomb lattice [75,85]. As is pointed out in Refs. [25,94], the fidelity susceptibility has a stronger singularity compared to the second-order derivative of the free energy and is thus a better indicator of quantum phase transitions. In particular, notice that the integrand of Eq. (30) differs from Eq. (7) by a multiplicative factor τ. One can show that in the critical region the singular part of χ F diverges faster than Eq. (30) by a factor of L z , where z is the dynamical critical exponent [22,25]. There are concrete examples in a class of topological phase transitions, which do not exhibit singularity in the secondorder derivative of the ground-state energy [99], but can still be detected using the fidelity susceptibility [100].
The covariance that appears in the estimators Eqs. (9) and (11) can also be written as where VarðxÞ ¼ hx 2 i − hxi 2 is the variance of x. This expression has an appealing meaning, i.e., the distributions of the vertices residing on the whole and on the halves of the imaginary-time axis have different widths, and the difference in these widths gives the estimate. This form resembles the bipartite fluctuation [101], which was proposed to be a diagnostic tool for phase transitions [102] because of its relation to the entanglement entropy. However, there are important differences. First, the fidelity susceptibility estimator requires a division in the imaginary-time axis for vertices, not in the real space for the physical particles. Second, the total number of vertices is fluctuating in the QMC simulations as opposed to being conserved in the case of bipartition fluctuations. Third, it is easier to locate the critical point using the fidelity susceptibility. As is shown in this paper and in many previous studies [16], the fidelity susceptibility exhibits an increasingly sharp peak at a phase transition as the system size enlarges. On the other hand, to utilize bipartite fluctuations and entanglement entropy for phase transition, one typically needs to resolve the scaling or subleading behavior with the system size, which is often difficult in finite size simulations.

VI. OUTLOOK
We present a general approach to compute the fidelity susceptibility of correlated fermions, bosons, and quantum spin systems in a broad class of quantum Monte Carlo methods [45][46][47][48][49][50][51][52][53][54][55][56]. The calculation of the fidelity susceptibility is surprisingly simple yet generic. It provides a general purpose indicator of quantum phase transitions without the need for a prior knowledge of the local order parameter.
Conceptually, our work shows it is rewarding to view the modern QMC methods [45][46][47][48][49][50][51][52][53][54][55][56] in a unified framework provided by Eq. (8), which deals with the same type of classical statistical problem irrespective of microscopic details of the original quantum system. In the QMC simulations, a quantum phase transition manifests itself as a particle condensation transition driven by changing of the fugacity of the corresponding classical model. This connection suggests generic ways to detect and characterize quantum phase transition through studying classical particle condensations. For example, Eq. (30) actually relates the second-order derivative of free energy of a quantum system to the particle compressibility of a virtual classical system. In this respect, the significance of the covariance estimators Eqs. (9) and (11) is evident because they capture the key critical fluctuation upon a particle condensation transition.
It is straightforward to generalize our Eqs. (9) and (11) to cases with multiple driving parameters, where one needs to count the vertices of different types (as is already done in the SSE calculations in Secs. II B 2 and IV B). It is interesting to find out whether this can lead to a general approach to measure the Berry curvature (the imaginary part of the quantum geometric tensor) in quantum Monte Carlo simulations. Related to these efforts, the nonequilibrium QMC method has been developed in recent years [103,104] to study the nonadiabatic response of quantum systems in imaginary time. In particular, it also allows the extraction of the fidelity susceptibility and the Berry curvature [103,105]. It would be interesting to compare the nonequilibrium QMC approach [103,104] to the equilibrium one presented in this paper.