Footprints of impurity quantum phase transitions in quantum Monte Carlo statistics

Interacting single-level quantum dot connected to BCS superconducting leads represents a well-controllable system to study the interplay between the correlation effects and the electron pairing that can result in a $0-\pi$ (singlet-doublet) quantum phase transition. The physics of this system can be well described by the impurity Anderson model. We present an analysis of the statistics of the perturbation expansion order of the continuous-time hybridization expansion quantum Monte Carlo algorithm in the vicinity of such a first-order impurity quantum phase transition. By calculating the moments of the histograms of the expansion order which deviate from the ideal Gaussian shape, we provide an insight into the thermodynamic mixing of the two phases at low but finite temperatures which is reflected in the bimodal nature of the histograms.


I. INTRODUCTION
The advances in the fabrication of nanodevices enabled to connect quantum dots with superconducting leads, forming superconducting quantum dot nanostructures. Many experimental realizations utilizing carbon nanotubes, semiconducting nanowires or even single molecules as central active elements demonstrate the large versatility of the concept which attracted a lot of attention over the last 20 years (for the overview see Refs. [1][2][3]).
In many cases the system can be well described by the single-impurity Anderson model (SIAM) coupled to superconducting BCS leads [4]. This model may exhibit, depending on the parameters, a so-called 0 − π transition signalled by the sign reversal of the supercurrent. This transition is governed by the underlying quantum critical point (QCP) and can be induced by changing the experimentally controllable model parameters as local energy level (gate voltage) or superconducting phase difference (magnetic flux in SQUID geometry) as well as by changing parameters that are hard to control in experiment, e.g. the local on-dot Coulomb interaction or the tunneling rate between a lead and the dot.
The usual method of choice for solving the superconducting SIAM is the numerical renormalization group (NRG) method [5]. It works at zero and low temperatures and provides direct access to the one-particle spectral function. The drawback of this method is the inability to solve this model using reasonable computational resources for more than two lead channels, e.g., in a case with an additional metallic lead. Despite the recent progress in this area [6] it still hinders the usability of NRG in multi-terminal and multi-dot setups.
Another class of numerically exact methods to solve SIAM are the quantum Monte Carlo (QMC) techniques. Many flavours of QMC were already employed to solve the superconducting SIAM including the Hirsch-Fye [7] and the continuous-time implementations in either interaction expansion (CT-INT) [4,8] or hybridization expan-sion (CT-HYB) [9,10] formalism. All these algorithms work in imaginary time domain and are therefore bound to finite temperatures with unpleasant scaling properties with decreasing temperature. Also, obtaining the spectral function requires to perform an analytic continuation to real frequency domain from stochastic imaginary time data which is a notorious ill-defined problem [11]. Nevertheless, QMC can provide unbiased finite-temperature results and invaluable insight into the thermal effects.
The most interesting effects in superconducting quantum dots are realized in the vicinity of the 0 − π impurity quantum phase transition (QPT). This transition is of the first order and is induced by a crossing of the two lowest-lying many-body eigenstates -a spin-singlet and a spin-doublet [12]. A convenient tool to study a QPT is the so-called fidelity which is a measure of the overlap of ground states of quantum systems sharing the same Hamiltonian but associated with different Hamiltonian parameters [13][14][15][16]. This quantity is by definition very sensitive to the change of the ground state induced by a variation of a control parameter. A method how to extract fidelity and its second derivative w.r.t the change of the control parameter -the fidelity susceptibility, from various QMC methods was used to study QPT in bulk (in Hubbard and Heisenberg models [17,18]) as well as impurity (in the two-orbital Anderson model [19]) systems. Unfortunately the calculation of the fidelity susceptibility is not implemented by default in most of the publicly available QMC solvers (except the iQIST package [20]).
In this paper we show that a different and valuable insight into the properties of a system close to a QPT can be obtained by studying the statistics of the perturbation expansion order from a CT-HYB calculation. We utilize the CT-HYB method to study the 0 − π transition in superconducting impurity Anderson model (SCIAM). This scenario is rather different from QPT in bulk models (Hubbard or Heisenberg) as we are dealing with an impurity QPT in a zero-dimensional system where the usual concepts like correlation length and finite-size scaling are not defined.
We quantify the unusual non-Gaussian shape of the histograms of the expansion order using the first four moments: mean, variance, skewness and excess kurtosis and discuss their behavior by changing the control parameters and the temperature. We show how the position of the QPT at zero temperatures can be extracted from the average expansion order and how is the thermodynamic mixing of the two phases reflected in the higher moments.
The paper is organized as follows. In Sec. II we introduce the SCIAM and the transformation that allows us to use CT-HYB to solve it. Sec. III introduces the moments of the expansion order and explains how they describe the statistics of a random variable. In Sec. IV we present and discuss results on the behavior of the system using the local energy as a control parameter and their connection to previous results from Ref. [21]. We explain the unusual bimodal shape of the histogram of the expansion order and how is it reflected in the moments. The main points are then summarized in Sec. V. Appendix. A presents more data obtained by changing other control parameters, namely the interaction strength and the superconducting phase difference, showing the universality of the results from Sec. IV.

A. Model formulation
The Hamiltonian of a spinful, single-level quantum dot connected to left (L) and right (R) phase-biased superconducting leads reads where the quantum dot is described by where d † σ creates an on-dot electron with spin σ ∈ {↑, ↓}, ε σ = ε − σB is the local energy level, B is the magnetic field, U is the on-site Coulomb interaction and the local energy is shifted in a way that ε = 0 always represents a half-filled dot. The non-interacting electrons in the leads are described by BCS Hamiltonians where c † αkσ creates an electron in lead α ∈ {L, R} with spin σ and energy ε αkσ , and ∆e iϕα is the complex superconducting order parameter for lead α with amplitude ∆ and phase ϕ α . Here we assume that the amplitude is the same in both leads (i.e. they are made from the same material) but the phases can differ. We note that physical observables can only depend on the phase difference ϕ = ϕ L − ϕ R and not on the individual values, i.e. the system must be invariant w.r.t. a global phase shift ϕ α → ϕ α + ϕ sh [3]. We also denotedk = −k to save space in equations. The coupling between the dot and the leads is governed by where V αkσ is the tunnel matrix element. As Hamiltonian (1) does not conserve charge, we cannot use the standard CT-HYB algorithm right away to solve it. Therefore we perform a canonical particle-hole transformation in the spin-down segment of the Hilbert space. Following Ref. [8] we define which maps the SCIAM (1) on the SIAM with nonsuperconducting bath but with attractive on-site interaction strength −U . We can define Nambu-like spinors in this basis as From now on we drop the spin index in the dispersion ε αk and the tunnel coupling element V αk . Let us define matrices Hamiltonian (1) can be rewritten in the form where σ z is a Pauli matrix. The occupation numbers transform as where n σ = d † σ d σ is the spin-resolved electron density, n is the total charge, m is the magnetization and ν = d ↓ d ↑ is the on-dot induced pairing.
From now on we assume the tunneling density of states to be constant in the relevant energy window, where W is the half-bandwidth of the non-interacting band. As the relevant energy scale for this model is the superconducting gap ∆ which is usually of order of 100µeV, the assumption of constant DoS over such a narrow energy window is justifiable in most experimentally relevant situations. Equation (10) also defines the tunneling rates Γ α which, for k-independent tunnel matrix elements V can be expressed as Γ α = π|V α | 2 /(2W ). We will focus only on the situations with symmetric coupling Γ L = Γ R = Γ/2 as any asymmetric setup with Γ L = Γ R can be mapped on the symmetric one using a simple formula [22].
The total effect of a lead α on the impurity can be summed into a hybridization function which describes the hopping from the impurity to the lead α, the propagation through the lead and the hopping back to the impurity. It can be represented in the imaginary-frequency (Matsubara) domain as where I 2 is the 2 × 2 unit matrix, ω n = (2n + 1)π/β is the n-th Matsubara frequency and β = 1/(k B T ) is the inverse temperature. Under the assumption of a constant tunneling density of states (10), the hybridization function of the lead α ∈ {L, R} can be expressed as where A(iω n ) = (2/π) arctan(W/ ω 2 n + ∆ 2 ) is the correction due to finite bandwidth that approaches unity for W → ∞.
The so-called 0 − π quantum phase transition corresponds to the change of the ground state from a nonmagnetic singlet (0-phase) to a spin-degenerate doublet (π-phase). It is signalled at zero temperatures by the abrupt change of the sign of both the Josephson current and the on-dot induced pairing from positive in 0-phase to negative in π-phase [23][24][25], and the crossing of the subgap, Andreev bound states at the Fermi energy [26][27][28]. As this transition is of the first order, the finitetemperature state is a simple thermodynamic mixture of the two phases and the transition is smeared out by increasing temperature into a crossover.
The method of extraction of the position of the QPT at zero temperature from finite temperature (experimental or QMC) data was already discussed in Ref. [21]. For low temperatures, only the lowest-lying many-body states are significant. Due to the presence of a superconducting pairing, these states are discrete and separated from the continuous part by a gap of the order of ∆. In the absence of magnetic field, we are dealing with a one spin singlet and one spin doublet. The second discrete singlet state and the continuum of states lie higher in energy and can be neglected at low temperatures. The canonical average of an observable F then reads where Z is the partition function, F s (F d ) are zerotemperature values of F in singlet (doublet) state and x is a control parameter (e.g. local energy ε), changing which we cross the QPT at critical value x c . At QPT, the two discrete many-body states cross, i.e. E s (x c ) = E d (x c ) and the average of an observable F reads This value is independent of temperature, which implies that all curves F (x) for low enough temperatures cross at one point and this point marks the position of the QPT at zero temperatures.

III. MOMENTS OF THE EXPANSION ORDER
A complete overview of the CT-HYB method can be found, e.g., in Ref [29]. It is based on an expansion of the partition function in powers of the hybridization term H hyb (4) around H dot +H lead . The order of the expansion k (the number of hybridization events in the given random configuration) depends rather strongly on the model parameters. It scales linearly with inverse temperature β and in general decreases with increasing interaction strength U . The statistics of k can be accumulated during the Monte Carlo run and represented in a form of a histogram H(k). The moments of k can be then calculated and their analysis provides a deeper insight into the performance of the solver and also the properties of the solved model.
We focus on four moments of the expansion order: mean, variance, skewness and excess kurtosis. The behavior of the variance and the excess kurtosis of k around a phase transition point has been briefly discussed by Huang et al. in an appendix of Ref. [18] for the case of the Mott transition in the single-band Hubbard model within the dynamical mean-field theory (DMFT). As the 0 − π impurity QPT in the superconducting Anderson model is a way simpler scenario than a Mott transition in the Hubbard model, it provides an opportunity to examine their behavior in a well-controllable setup and shed more light on the statistics of the expansion order.
The n-th central moment of a discrete random variable k ∈ {0, 1, 2 . . .} is defined as where p(k) is the number of occurrences of given k in the sample, c n is the n-th central value and N is the number of measurements. The mean µ = µ 1 is the first moment about zero (c 1 = 0), and variance as the second moment about the mean (c 2 = µ), For higher moments it is useful to define standardized moments (normalized moments about the mean)μ n = µ n /σ n . Skewness and kurtosis are then defined as the third and fourth standardized moments, γ =μ 3 and κ 0 =μ 4 , respectively. As the kurtosis of the Gaussian distribution equals 3, it is often convenient to define the excess kurtosis κ = κ 0 − 3. To summarize, The first two moments are standard tools in mathematical statistics and their meaning is clear, but skewness and kurtosis are more abstract quantities. Therefore we illustrate them in Fig. 1. Skewness measures the asymmetry of the tails of a distribution. In panel (a) we plotted distributions with skewness γ > 0 (left-leaning, right tail is longer), γ = 0 (no skew) and γ < 0 (rightleaning, left tail is longer). The excess kurtosis measures how fast the distribution vanishes at higher values compared to the Gaussian. Distributions with κ > 0 (with "fat tails") are usually referred to as leptokurtic, with κ = 0 as mesokurtic and with κ < 0 ("thin tails") as platykurtic. Examples of such distributions are plotted in panel (b).
It is worth noting that an alternative way how to describe a distribution, complementary to the moments, are the statistical cumulants. The first four cumulants K n of a random variable k are connected to the standardized moments as defined in Eq. (17) as K 1 = µ, K 2 = σ 2 , K 3 = γσ 3 and K 4 = κσ 4 and the cumulant analysis would be analogous to the moment analysis presented below.
The averaged expansion order µ from CT-HYB also has a physical meaning. In the case of the impurity Anderson model it is connected to the hybridization energy as k = βE hyb where [30] This is one of the differences between a calculation of the impurity Anderson model and a DMFT calculation of the Hubbard model, for which the averaged expansion order from CT-HYB can be identified with the kinetic energy.  This formula can be evaluated numerically and used as a benchmark of the accuracy of the measurement of the local Green function w.r.t. the averaged expansion order.

IV. RESULTS
All CT-HYB calculations were performed using the TRIQS/CTHYB 2.2 solver [32]. We set B = 0, W = 100∆ and the cutoff in Matsubara frequencies ω max n ≥ 300∆. We encountered no fermionic sign problem during the calculations. The total charge n and the induced pairing ν were evaluated by tracing the measured impurity density matrix and using Eqs. (9).
In Fig. 2 we plotted the on-dot induced pairing ν as a function of the local energy ε for U = 8∆, Γ = ∆ and ϕ = 0 for three values of inverse temperature β∆ = 40, 20 and 10. Here ε = 0 represents the half-filled dot. This plot illustrates the typical behavior of SCIAM around the QPT. The induced paring is positive in the 0-phase and negative in the π-phase and all lines cross at ε c ≈ 2.67∆. The position of such a crossing is consistent with the position of the QPT at zero-temperature as explained in Ref. [21]. The inset shows the total on-dot electron density n that decreases as we move away from the halffilled case. Again, all lines cross at the same value of ε, proving the applicability of the formula (14).
The histograms of the expansion order k normalized to the number of QMC measurements (4.8 × 10 8 in this case) for β∆ = 40 are plotted in Fig. 3(a) We plotted few histograms around the critical value ε c (solid lines) together with histograms for values further from the transition point, ε = ∆ and ε = 4∆ (dashed lines) for comparison. As the averaged expansion order is rather large, the histograms resemble more a continuous curve than a discrete set of values. It is also the reason why the histograms are not distorted much by the fact that the expansion order is bounded from below by zero. As the averaged expansion order increases with decreasing temperature, it is always possible to run the simulation at low enough temperatures to minimize the distortion.
We see that histograms for ε = ∆ and 4∆ are of the Gaussian shape. This shape does not change until we approach the critical point. As we move closer, the histogram starts to deviate strongly from a Gaussian, first developing a shoulder and later turning into a broad two-peaked structure. Similar double-peaked histograms were also reported in Ref. [18] for the Hubbard model in the vicinity of the Mott transition. The histograms for values of ε close to the critical value also intersect at one point at k ≈ 99. As the two-peak structure can be well fitted by a pair of Gaussians, we can assume these histograms can be also seen as thermal averages over the singlet and doublet ground states, We denote the intersection point of the two histograms as k 0 , so H s (k 0 , ε c ) = H d (k 0 , ε c ) where we replaced the local energy by its critical value, assuming that the constituent histograms H s and H d don't depend much on ε around the critical point. As E d (ε c ) = E s (ε c ), we obtain in the vicinity of the critical point (ε ≈ ε c ), for all values of ε close enough to the critical point that we can replace H s,d (k, ε) ≈ H s,d (k, ε c ). This feature represents itself as the histograms crossing at the same point at k = k 0 for ε ≈ ε c . In panels (b)-(e) we plotted selected histograms in the vicinity of the critical point fitted by a pair of Gaussians. This illustrates the feature that the constituent histograms H s and H d don't depend much on ε, only their weights are changing as we cross the transition point. The first four standardized moments (mean, variance, skewness and excess kurtosis) as defined by Eq. (17) calculated from the histograms for the same parameters as in Figs. 2,3 are plotted in Fig. 4. These quantities were calculated using Eq. (15) where we approximated p(k) by the histogram H(k). Panel (a) shows the behavior of the average expansion order µ scaled to inverse temperature β for three values β∆ = 40, 20 and 10. All curves cross at ε c , similarly as curves for the induced pairing. This is not a surprise since we identified the physical meaning of µ/β as the hybridization energy and at low temperatures Eq. (18) also reduces to a relation in a form of Eq. (14). It proves that the position of the QCP is en- coded already in the statistics of the expansion order and can be obtained without actually measuring any physical observable or the one-particle Green function.
We can obtain the same result that lines µ(ε)/β for different low enough temperatures cross at the same point using the assumption that the histogram in the vicinity of the critical point can be decomposed into parts from the singlet and from the doublet, H(k, ε ≈ ε c ) = [H s (k, ε) + 2H d (k, ε)]/3. Then Eq. (15) for n = 1 and c 1 = 0 reduces to a weighted sum of two contributions, µ ≡ k = ( k s + 2 k d )/3 where k s,d is an average of k w.r.t. H s,d .
The broadening of the histogram can be quantified by the scaled variance σ 2 /β plotted in Fig. 4(b) that shows a peak above ε c . For β∆ = 40 the peak lies at ε ≈ 2.72∆ (between the blue and green curves in Fig. 3(a)). The curves for different temperatures do not cross at the same point. Using the same reasoning as for the averaged expansion order we obtain The first term scales as β with temperature while the second scales as β 2 . As a result, curves of scaled variance σ 2 /β at the transition point don't cross as there is a linear offset due to the last term in Eq. (22). As the width of the crossover region between the two phases is increasing with increasing temperature, the peak in variance becomes wider, less pronounced and the maximum moves away from ε c . As this maximum lies away from the transition point at any finite temperature, variance is not a good estimator of the position of the QPT. We would like to point out that the fact that the scaled average expansion order µ/β can be identified with the hybridization energy E hyb does not automatically imply that the variance of the expansion order is identical to the variance of E hyb . The question whether the variance has a physical meaning is beyond the scope of this article.
The overall shape of the histogram can be well quantified by the skewness γ which is plotted in Fig. 4(c). Histograms for low temperatures have almost always a slight positive skewness (leaning to the left). This is probably due to the fact that the expansion order k is bound from below but not from above. As we approach the QPT by moving away from half-filling (from the πphase to the 0-phase), skewness shows a prominent positive peak connected with the development of the shoulder on the histogram indicating increasing presence of the 0phase in the thermal mixture. The skewness is zero at the same point where variance has a maximum, above the zero-temperature QPT. This is the point where the histogram is symmetric due to equal admixtures of the zero and the π phases at given temperature. As we move across this point, the scenario is reversed with a negative skewness indicating decreasing presence of the π phase. All the features are quickly smeared out by increasing temperature and for β∆ = 10 the skewness is positive for all values of ε.
The curves for different temperatures seem to cross at one point very close to ε c within the QMC error bars. As skewness in our definition is normalized to σ 3 , it is not possible to perform a simple analysis as in the case of the variance and the reason of this behavior remains unknown.
The excess kurtosis κ provides additional insight into the deviations from a Gaussian shape. It is, in contrast to the skewness, insensitive to which tail of the distribution is deviating from the Gaussian, therefore the curve is rather symmetric around ε c . Excess kurtosis is almost zero at low temperatures and far from the transition point indicating that the tails of the histogram are of Gaussian shape. As we move closer to QPT from any side, histogram becomes leptokurtic (fat-tailed), again because of the presence of the shoulder, as evident from Figs. 3(b),(e). This indicates that configurations with large or small number of hybridization events are encountered here more often during the QMC sampling. Then it shows a well developed minimum at the same point where skewness crosses zero, i.e. where the histogram is symmetric.
In summary, the abovementioned statistical moments quantify well the shape of the histogram of the expansion order. It changes its shape from nearly ideal Gaussian to bimodal which can be fitted by a pair of Gaussians as a result of the thermodynamic mixing of the two phases at low but finite temperatures. The position of the zerotemperature QCP can be extracted from the crossing of the curves of the first moment for two low-enough temperatures. The extremal points of higher moments (maximum of σ 2 , zero of γ, minimum of κ) lie always away from the transition point in the direction of the 0-phase and approach it only slowly with decreasing temperature, therefore they are not good estimators of the position of the QCP, although they provide a way to quantify the mixing of the phases in the crossover region.

V. CONCLUSIONS
Superconducting quantum dots provide a simple testbed for studying impurity QPT and this physics can be captured very well by the SCIAM. Numerically exact CT-HYB QMC simulation can provide an unbiased insight into finite-temperature behavior in the vicinity of the critical point. The change of the ground state from non-magnetic singlet to spin-degenerate doublet leaves a unique footprint on the statistics of the hybridization expansion order which is accumulated during the simulation and can be quantified by the moments of its histogram.
The position of the zero-temperature QCP can be extracted from finite-temperature results by measuring two sets of data for any physical observable at two different low-enough temperatures as already discussed in Ref. [21]. As the average expansion order in CT-HYB can be identified with the hybridization energy, it is also a physical observable and therefore position of the QCP can be extracted from its statistics without performing any QMC measurement at all (except a trivial measurement of the average sign). In the vicinity of the critical point the histogram of the expansion order deviates from the ideal Gaussian shape and this deviation can be quantified by the higher moments (variance, skewness, excess kurtosis). We show that the histogram can be separated into Gaussian contributions from the two ground states. The weights of these contributions equalize at the point that linearly approaches the QCP with decreasing temperature and can be identified from the maximum of variance, zero of the skewness or the minimum of the excess kurtosis.
We believe this analysis can provide additional insight into the physics around an impurity QPT beyond the scope of standard tools like the fidelity susceptibility for systems that can be mapped onto an impurity Anderson model. Furthermore, this analysis can be performed using any available CT-HYB solver without any modification or implementation of new features. The results also hold for any system exhibiting a first-order impurity QPT as well as bulk systems that are treated using the DMFT technique of mapping a lattice model on an effective impurity problem.
Appendix A: Dependence on the interaction strength and phase difference A superconducting quantum dot can be driven through a QPT by changing several control parameters. Here we provide additional data on the behavior of the histograms and moments of the expansion order as functions of the interaction strength U (Fig. 5) and phase difference ϕ (Fig. 6). The behavior is analogous to the case discussed above, except that here we move in both cases from zero to the π-phase by increasing the parameter value. In both cases, panel (a) shows the induced pairing for three different temperatures. The curves cross at the same point proving the applicability of the two-level model for the selected temperatures. We added the available zerotemperature NRG result to Fig. 5 to further illustrate this feature. Panel (b) shows histograms for β∆ = 40 in the vicinity of the QPT (solid lines) that also cross at one point, supplemented by Gaussian-like histograms from calculations away from the QPT (dashed lines). Panels (c-e) show the behavior of the first four standardized moments as defined in Eq. (17). The analysis is again analogous to the ε-dependent results from the main text. skewness (e) and excess kurtosis (f). Lines are splines of QMC data and serve only as a guide for the eye. The critical phase difference ϕc ≈ 0.48π matches the available zero-temperature NRG data (not plotted).