Bose-Fermi Anderson Model with SU(2) Symmetry: Continuous-Time Quantum Monte Carlo Study
Ang Cai Qimiao Si





Department of Physics and Astronomy, Rice Center for Quantum Materials, Rice University,
Houston, Texas, 77005, USA


Department of Physics and Astronomy, Rice Center for Quantum Materials, Rice University,
Houston, Texas, 77005, USA


In quantum critical heavy fermion systems, local moments are coupled to both collective spin fluctuations 
and conduction electrons. As such, the Bose-Fermi Kondo model, describing the coupling of a local moment
to both a bosonic and a fermionic bath, has been of extensive interest. For the model in the presence 
of SU(2) spin rotational symmetry, questions have been raised about its phase diagram.
Here we develop a version of continuous-time Quantum Monte Carlo (CT-QMC) 
method suitable for addressing this issue; this procedure
 can reach sufficiently low temperatures
while preserving the SU(2) symmetry.
Using this method for the Bose-Fermi Anderson model,
we clarify the renormalization-group fixed points and the phase diagram for the case 
with a constant fermionic-bath density of states and
  a power-law bosonic-bath spectral function ρ_b(ω) ∝ω^s (0<s<1).
We find two types of Kondo destruction QCP, depending on the power-law exponent s in the bosonic bath spectrum.
For s^*<s<1, both types of QCPs exist and, in the parameter regime accessible by an analytical
ϵ-expansion renormalization-group calculation (here ϵ=1-s), the  CT-QMC
result is fully consistent 
with prior predictions by the latter method. 
For s<s^*, there is only one type of QCP.
At both type of Kondo destruction QCPs, 
we find that the exponent of the local spin susceptibility 
η obeys the relation η=ϵ, which
 has important implications
for Kondo destruction QCP in the Kondo lattice problem.


Bose-Fermi Anderson Model with SU(2) Symmetry: Continuous-Time Quantum Monte Carlo Study
    Qimiao Si
    December 30, 2023
========================================================================================



§ INTRODUCTION

Heavy fermion systems serve as a prototype system to study quantum criticality 
<cit.>.
 Experimental discoveries in various heavy fermion compounds 
open up the opportunity to explore
beyond-Landu type quantum critical points (QCP) 
in the context of antiferromagnetic Kondo lattice systems. One prominent example is the Kondo destruction QCP 
<cit.>,
where the phase transition 
at zero temperature not only involves the 
magnetic order parameter, but also the localization to delocalization transition of the 4f electrons constituting  the local moments. Some of the hallmarks of Kondo destruction type QCP involves ω/T scaling of the dynamical spin susceptibility as seen from inelastic neutron scattering, jump of the fermi surface volume from magnetotransport and quantum oscillation measurement <cit.>.
Such properties 
 are inconsistent with predictions from the traditional spin-density-wave type QCP <cit.>.

One of the simplest models that contain a Kondo destruction type QCP is the Bose-Fermi Kondo model (BFKM)
<cit.>. It arises in the context of understanding the competition between Kondo effect and 
magnetic fluctuations in Kondo lattice model using extended dynamical mean field theory (EDMFT) <cit.>. It describes a local moment coupled to both itinerant electrons 
as well as free bosons, which are usually referred to as fermionic bath and bosonic bath. 
Typically the fermionic bath will assume a constant density of states, and the bosonic bath has a sub-ohmic spectrum:
its density of states at low frequencies (ω) have a power-law form,
ρ_b(ω) ∝ω^s with s<1. It characterized the softened spectrum of the magnons 
near the magnetic QCP,  which competes with the conduction electrons in their couplings to the local moment
and causes the suppression of the Kondo effect.

This model is first treated with ϵ-expansion renormalization group (RG) method, 
using ϵ=1-s as a small parameter <cit.>. 
It turns out the fixed point structure will depend on the symmetry of the spin boson coupling: for the SU(2) and XY symmetric cases, it has a Kondo screened stable fixed point (K) at strong coupling, a bosonic bath dominated stable fixed point (L) at intermediate coupling (so called critical phase), and an unstable critical point (C) describing the quantum phase transition. Both L and C can be accessed by the ϵ-expansion; for the Ising anisotropic case, on the other hand, the critical phase controlled by L is unstable and is replaced by the local moment fixed point (L^') at strong coupling. In all three cases, it is predicted that at the critical point where the Kondo effect is critically destroyed, the local spin correlation function will behave as χ_spin(τ) ∼ (1/τ)^η, with an exact relation η=ϵ <cit.>. This has important implications for the EDMFT calculation of the Kondo lattice problem. For two dimensional magnetic fluctuations, it predicts a Kondo destruction QCP solution, provided that the relation η=ϵ will remain valid at ϵ→ 1^-.

The numerical calculations of the Bose-Fermi Kondo model and the closely related Bose-Fermi Anderson model (BFAM)
include treating it either as a standalone model using numerical renormalization group (NRG) <cit.> and continuous-time quantum Monte Carlo (CT-QMC) <cit.>, or as an effective model under EDMFT <cit.>.
Our focus in this work is on the CT-QMC method,
from which a seeming controversy existed for the 
SU(2) symmetric BFAM <cit.>:
for s=0.2, it was shown that 
the Kondo-destruction phase has the local-moment character instead of being critical;
 in the temperature dependence of the local spin susceptibility in this Kondo-destruction phase,
 it was  found χ^spin∼ 1/T instead of the χ^spin∼ 1/T^s behavior 
 predicted by ϵ-expansion RG<cit.> 
 for fixed point L.






To resolve this seeming inconsistency,
we start with the observation that, if s is close to 1,
 the CT-QMC result must be consistent with that of the 
 ϵ-expansion RG in the range of coupling constants 
accessed by this expansion (again ϵ=1-s).
 To make progress,
in this article we 
develop the CT-QMC procedure for the BFAM such that it can reach sufficiently low temperatures
while preserving the SU(2) symmetry.
Using this procedure,
 we  carry out
 a comprehensive study of the SU(2) BFAM
 for s ranging from close to 0 to close to 1.
 We 
 study 
a variety of
 observables in order to identify all the QCPs between different phases, 
 combined with detailed finite size scaling analysis to extract critical exponents. 
 
 Our analysis shows that the ϵ-expansion <cit.>
  and CT-QMC results are fully compatible with each other.
 Our results are summarized by the RG-flow diagrams of figure <ref>. 
 For the s>s^* regime, we identify i) the critical point C separating the Kondo screened phase and critical phase, as predicted from ϵ-expansion RG for the coupling constants accessible by the latter method; and ii)
 a separate critical point C^' and stable fixed point L^', which occurs for larger values 
 of the bosonic-Kondo coupling g.
 For s<s^*, there exists only a type ii) quantum phase transition <cit.>.
We also determine the correlation length exponent ν.
Additionally,
we find another unstable fixed point LC that controls the transition from fixed point L and fixed point L^'. 
Finally, we quantitatively estimate s^* and conclude that the result at s=0.2 falls outside the 
regime that is controlled by 
the ϵ-expansion.

The remainder of the paper is organized as follows. In Sec. <ref> we introduce the SU(2) 
Bose-Fermi Anderson model, and give an overview of the CT-QMC method as well as the physical 
quantities we will investigate in this work. We will present the numerical results in Sec. <ref>. 
We will start with a detailed study for the s=0.6 case in Sec. <ref>, 
followed by the s=0.2 case in Sec. <ref>,
before carrying through the analysis that leads to an estimate for the value of s^* in Sec. <ref>. 
We will discus the implication of our results in Sec. <ref> and conclude the article 
 in Sec. <ref>.



§ MODEL AND METHOD


The Hamiltonian for the SU(2) symmetric BFAM reads,

    H=H_c+H_b+H_d+H_g+H_V,

where H_c and H_b describes the bosonic and fermionic bath part, respectively,

    H_c=∑_k,σϵ_k c_k,σ^† c_k,σ , H_b=∑_α H_b^α=∑_p,αω_qϕ_p^α^†ϕ_p^α.

H_d contains the local electron part,

    H_d=∑_σϵ_d d_σ^† d_σ + U d^†_↑ d_↑ d^†_↓ d_↓ .

H_V and H_g couples the local orbital to the bosonic and fermionic bath,

    H_V=∑_k,σ  V d_σ^† c_k,σ + h.c. , 
    H_g=∑_p,α g S_α ( ϕ^α_p^† + ϕ_-p^α )

where the summation over α runs through x,y,z, S_α=d^†_στ^α_σσ^' d_σ^', and τ^α_σσ^' being the three components Pauli matrices.


The properties of the fermionic and bosonic bath are specified by their density of states. For the fermionic bath, we choose a constant density of states,

    ρ_F(ϵ)=∑_kδ(ϵ-ϵ_k) = ρ_0Θ (|D-ϵ|),

which leads to a hybridization function Γ(ϵ)= Γ_0Θ(|D-ϵ|), with Γ_0=πρ_0V^2. 

Unless specified otherwise, the density of states for the sub-Ohmic bosonic bath has an exponential cutoff, given by the following,

    ρ_b(ω)=∑_qδ(ω-ω_q ) = K_0ω^s e^-ω/ΛΘ(ω).


Throughout the text we fix D=1, Λ=1, and stays at the particle-hole symmetric point U=-2ϵ_d=0.1. 
The prefactor ρ_0 and K_0 in the density of states of the fermionic bath and bosonic bath are determined 
from the normalization condition ∫_-D^Dρ_F(ϵ) d ϵ = 1 and ∫_0^∞ρ_b(ω) 
dω =1. We will use either the amplitude of the hybridization function Γ_0 or the spin-boson 
coupling g as our tuning parameter.



 §.§ Monte-Carlo procedure

We will employ the CT-QMC algorithm, first introduced in reference <cit.> and then generalized 
to treat the BFAM in references <cit.>.
We start with removing the z component of the spin-boson coupling 
by employing a Firsov-Lang transformation H̃=e^S 
H e^-S with S=g S_z∑_p1/ω_p ( ϕ_p^z^† +ϕ_-p^z ) 
(similar to Ref. werner2007) and work with the transformed Hamiltonian 
 H̃,

    H̃   =    H_c+H_b+H̃_d+H̃_V+H̃_g
    H̃_d   =   ∑_σϵ̃d_σ^† d_σ + Ũd^†_↑d_↑d^†_↓d_↓
    H̃_V   =    V ∑_k,σ(d_σ^† c_k,σ 
    e^∑_pgs_σ/ω_p(ϕ^z_p^†-ϕ^z_p) +h.c.) 
    H̃_g   =   ∑_p (g / √(2) )
    (
    S_+ϕ_p^- e^∑_pg /ω_p(ϕ^z_p^†-ϕ^z_p) .
    
       +   .
    S_-ϕ_p^+ e^-∑_pg /ω_p(ϕ^z_p^†-ϕ^z_p) ),
 
where we have defined the renormalized parameters ϵ̃_d=ϵ_d-(g^2/4) ∑_q (1/ω_q)^2, Ũ=U+(g^2/2) ∑_q (1/ω_q)^2, s_σ=± 1/2 for σ=↑/↓.
and recombined the x and y components of S_α and ϕ_α into
S_+=d^†_↑ d_↓, S_-=d^†_↓ d_↑, ϕ_p^±=(1/√(2))((ϕ_p^x^†+ϕ_p^x)± i (ϕ_p^y^†+ϕ_p^y)). 

The partition function is constructed by expanding in the non-diagonal terms
<cit.>,
 H̃_V and H̃_g 
under the interaction representation of H_0≡ H_b+H_c+H̃_d. It has the following form
<cit.>:

    Z   =   Z_0∑_m∫∏_i=1^m dτ^s_i dτ^s^'_i∏_σ=↑,↓( ∫∏_i=1^n_σdτ^dσ_i
    dτ^d^'σ_i)
    
          
    w_d({τ^tot}_n_tot )
    ∏_σ=↑,↓
    w^σ_c({τ^dσ}_n_σ,{τ^d^'σ}_n_σ)
    
           
    w_z({τ^tot}_n_tot )
    w_p({τ^s}_m,{τ^s^'}_m),

where Z_0=Tr[e^-β H_c] Tr[e^-β H^z_B] Tr[e^-β (H^x_B+H^y_B)] is the partition function of the bath, β being the inverse temperature: β=1/T.
∫∏_i=1^m dτ^α_i dτ^α^'_i=∫_0^β dτ^α_1⋯∫_τ^α_N-1^β dτ^α_N∫_0^β dτ^α^'_1⋯∫_τ^α^'_N-1^β dτ^α^'_N.
{τ^α}_n denotes the set of imaginary time of all the operators of a given type α in the expansion: {τ^α}_n={τ^α_1,τ^α_2…,τ^α_n}. α∈{s,s^',dσ, dσ^'} represents S_+, S_-, d^†_σ, or d_σ. n=m or n_σ denotes the number of pairs of S_+, S_- or d^†_σ, d_σ, also labeling the expansion order.
{τ^tot}_n_tot refers to all the {τ^α}_n combined,



with n_tot=2(∑_σ n_σ + m). The integrand, or so-called weight, factorizes into multiple components. In the following we will present the form of each part explicitly.

w_d({τ^tot}_n_tot ) is the contribution from the local d electron part. It describes valence and spin fluctuations of the local orbitals,

    w_d   =   Tr[e^-βH̃_dT_τ  S_-(τ^s^'_m) S_+(τ^s_m) ⋯ S_-(τ^s^'_1) S_+(τ^s_1)  
       ×   ∏_σ d_σ(τ^d^'σ_n_σ) d_σ^†(τ^dσ_n_σ) ⋯ d_σ(τ^d^'σ_1) d_σ^†(τ^dσ_1)].

Here for a given operator O, O(τ) denotes the corresponding operator in the interaction representation O(τ)=e^τ H_0 O e^-τ H_0.

w^σ_c({τ^dσ}_n_σ,{τ^d^'σ}_n_σ) is the contribution from the conduction electron with spin index σ,

    w^σ_c   =    V^2 n_σ( ∏_i=1^n_σ∑_k_i,k^'_i) Tr[T_τ e^-β H_c  c^†_k_n_σ,σ(τ^d^'σ_n_σ)
    
       ×    c_k^'_n_σ,σ(τ^dσ_n_σ) ⋯ c^†_k_1,σ(τ^d^'σ_1)
    c_k^'_1,σ(τ^dσ_1)]/Tr[e^-β H_c]
    
       =   (F^σ).

It can be expressed as a determinant of matrix F^σ, whose matrix element is given by 

    F^σ_ij=-∑_kV^2 Tr[e^-β H_c T_τ c_k,σ(τ^dσ_j) c^†_k,σ(τ^d^'σ_i)]/Tr[e^-β H_c].


w_z({τ^tot}_n_tot ) comes from the z component bosonic bath part <cit.>,

    w_z   =   Tr[e^-β H_B^z∏_i=1^n_tot  e^s_i∑_p (g^z/ω_p) (ϕ_p^z^†(τ^tot_i)-ϕ_p^z(τ^tot_i ))  
      ]/Tr[e^-β H_B^z ]
       =   exp( -g^2∑_1<i<j<n_tot s_is_j(B(τ_i-τ_j)-B(0))   )
      ,

where s_i=± s_σ or ± 1 when the operator O(τ_i^tot) at τ_i^tot corresponds to d_σ^† / d_σ or S^±, and

    B(τ_j-τ_i)   =   ∑_p Tr[T_τ e^-β H_B^zϕ^z_p(τ_i) ϕ^z_p^†(τ_j) ]/ω_p^2 Tr[e^-β H_B^z ]
       +    (τ_i↔τ_j).


Finally,
w_p({τ^s}_m,{τ^s^'}_m) involves the bosonic bath in the transverse direction <cit.>,
forming a permanent,

    w_p   =   (g/√(2))^2m( ∏_i=1^m ∑_p_i,p^'_i) Tr[e^-β(H_B^x+H_B^y)T_τϕ_p_m^+(τ^s^'_m)
    
       ×   ϕ^-_p^'_m(τ^s_m)⋯ϕ^+_p_1(τ^s^'_1)ϕ^-_p^'_1(τ^s_1) ]/Tr[e^-β(H_B^x+H_B^y) ]
    
       =   ∑_p ∈ S_m∏_i=1^m P_i,p(i).

The summation extends over S_m, representing all permutations of 1,2,⋯,m. 
The matrix element of P is the following,

    P_ij   =   (g^2/2)∑_pTr[ e^-β(H_B^x+H_B^y) T_τϕ^-_p(τ^s_j)  ϕ^+_p(τ^s^'_i)]
    /Tr[e^-β(H_B^x+H_B^y)]
       ≡   (g^2/2) J(τ_j^s-τ_i^s^').


Now the partition function can be interpreted as integrating a probability distribution function over some configuration space. Here, each configuration is specified by all sets of different {τ^α}_n and a particular permutation p ∈ S_m, which is then sampled through a Metropolis algorithm with a probability proportional to w_d× w_z× w_c^↑× w_c^↓×∏_i=1^m P_i,p(i). 






We now describe the Monte Carlo updates. We inherit the updates from the Ising BFAM
<cit.>, namely
 the insertion, removal and shift of  d_σ^† c_k,σ/c^†_k,σ d_σ pair,
 and also  adopt the insertion/removal of S^+ϕ^-/S^-ϕ^+ and the sampling of the permutation 
 S_m introduced in reference <cit.> (named updates (a)-(c) there).
 In addition we introduce 
 a swap update that swaps
  S_+ (S_-) with a pair of d_↑ ^† and d_↓ ( d_↓ ^† and d_↑  ). For example consider the S_+ case. We first randomly pick a pair of S_+(τ^s_i), S_-(τ^s^'_j) from the m pairs of S_+ and S_- that is connected by one of J(τ). Then we choose a d_↑^†(τ^d↑_k) with a probability P_k=J (τ^d↑_k-τ^s^'_j) / ( ∑_n=1,n_↑ J(τ^d↑_n - τ^s^'_j ) ) from the n_↑ of d^†_↑ operators. We then swap the position of  S_+(τ^s_i) and d_↑^†(τ^d↑_k). Finally, we find the d_↓ (τ^d^'↓_l) that is closest to d_↑^†(τ^d↑_k) before the swap, and move it to d_↓ (τ^d^'↓_new). τ^d^'↓_new is randomly selected within an interval of length l_max, which is the distance between two creation operators in the σ=↓ orbital next to S_+ before the swap. The corresponding proposal probability is given by

    P_prop= 1/l_max m × J (τ^d↑_k-τ^s^'_j) /∑_n=1^n_↑ J(τ^d_n - τ^s^'_j ) .
 
Likewise we can find the proposal probability for the inverse update,

    P_prop^inv=  1/ l^'_max m× J (τ^s_i-τ^s^'_j ) /∑_ n=1, 
     n≠ k ^n_↑ J(τ^d_n - τ^s^'_j ) + J(τ^s_i- τ^s^'_j ) .
 

The weight ratio between the proposed configuration and the current configuration is given by

    w_new/w_old   =     w^↑_c( {τ^d↑}^new_n_↑ , {τ^d^'↑}_n_↑  ) w^↓_c( {τ^d↓}_n_↓ , {τ^d^'↓}^new_n_↓  )  /  w^↑_c( {τ^d↑}_n_↑ , {τ^d^'↑}_n_↑  ) w^↓_c( {τ^d↓}_n_↓ , {τ^d^'↓}_n_↓  )   
       ×   w_d( {τ^tot}^new_n_tot ) w_z( {τ^tot}^new_n_tot ) J (τ^d↑_k-τ^s^'_j) /w_d( {τ^tot}_n_tot )  w_z( {τ^tot}_n_tot ) J (τ^s_i -τ^s^'_j) .

where {τ^d↑}^new_n_↑ is {τ^d↑}_n_↑ with τ_k^d↑ replaced by τ_i^s, {τ^d↓}^new_n_↓ is {τ^d↓}_n_↓ with τ_l^d^'↓ replaced by τ_lnew^d^'↓, and {τ^tot}^new_n_tot is {τ^tot}_n_tot with the above two substitutions, plus τ_i^s replaced by τ_k^d↑.

The detailed balance condition is satisfied by the adopting the acceptance ratio max[R,1], with R given by

    R=w_new/w_old×P_prop^inv/P_prop.



The reason that we choose the proposal probability to be the form in equation (<ref>) and equation (<ref>) is to cancel out the J (τ^d↑_k-τ^s^'_j)/J (τ^s_i -τ^s^'_j) factor in the weight ratio in equation (<ref>), such that the acceptance ratio R is of order 1. Otherwise if we select d^†_↑ ( τ_k^d↑) using a uniform distribution from 0 to β, since J(τ) ∼ 1/τ^1+s, on average J(τ^d↑_k-τ^s^'_j) ∼ 1/ β^s, while the average value of J (τ^s_i -τ^s^'_j) is β independent, as a result R will be suppressed by a factor of 1/β^s. Similar ideas have been introduced in reference <cit.>.

In practice we have tested that the swap update introduced here can replace the role of update (d) in reference <cit.>, which breaks up one S^+ (S^-) into a pair of d^†_↑ and d_↓ ( d^†_↓ and d_↑ ) at two different time. Both of these updates are introducing shortcuts between configurations that are connected by a large number of other updates. But unlike update (d) whose acceptance ratio decreases with β as a power-law, the swap update has an acceptance ratio that does not depend on β. This facilitates the task of reaching low enough temperatures and access the scaling 
regime. We have verified that our procedure preserves the SU(2) symmetry.

We now make a few remarks on how to evaluate J(τ) and B(τ) in the numerical calculation. This is important because in the current expansion scheme the weight contribution from the bosonic bath in the transverse direction ϕ^± and in the z direction ϕ^z enters differently. Thereby the SU(2) symmetry of the model has to be recovered dynamically in the sampling process. In actual calculation we find that in order to maintain the SU(2) symmetry, it is crucial to evaluate B(τ) and J(τ) to sufficiently high accuracy.

Starting with the Fourier components of J(τ) in the matsubara frequency domain,

    J(iν_n)=∑_p 2 ω_p/ω_p^2-ν_n^2,

where ν_n=2π n/β, n ∈ℤ is the matsubara frequencies. There are two ways to calculate J(τ). We can either perform the integration over the density of states first,

    J(iν_n ) = ∫_0^∞2 ω/ω^2-ν_n^2ρ_b(ω) d ω,

followed by the matsubara summation,

    J(τ)=1/β J(i ν_n =0 ) + 2/β∑_ν_n>0 J(iν_n) cos (ν_nτ) .

Or we can first do the matsubara summation, then integrate over the density of states,

    J (τ) =∫_0^∞ e^(β-τ)ω + e^τω/ e^βω -1 ρ_b(ω) d ω .

In practice we find the summation in equation (<ref>) converges too slow when β is large. So using equation (<ref>) is recommended.

On the other hand, J(τ) is related to B(τ) by being its second derivative: J(τ)=d^2 B(τ) /d τ^2. B(τ) is most easily evaluated using the following formula,

    B(τ)-B(0) 
       =    J(iν_n=0) τ(τ-β)/2β 
    + ∑_n ≠ 0 J (iν_n) 1-cos(ν_nτ)/βν^2_n   .

Because of the extra 1/ν_n^2 factor here, the summation actually converges very quickly.






 §.§ Observables

In this subsection we introduce all the quantities we will calculate using CT-QMC.

We start with the local magnetization,

    ⟨ m_α⟩ = ⟨1/β∫_0^β S _ α (τ) dτ⟩ ,  α=x,y,z,

which is related to most of the quantities we discussed below.

Because the sampling will preserve spin rotation symmetry, the actual measured ⟨ m_α⟩ is always 0. Instead we measure its root mean square,

    σ_α= √(⟨ m_α ^2⟩),

which is also related to the static spin susceptibility 
χ^spin_α (T) = ∫_0^βχ^s_α(τ)  dτ=  ∫_0^β⟨ T_τ S _ α (τ)S _ α⟩  dτ
by,

    χ^spin_α=βσ_α^2,

where we have also defined the dynamical spin correlation function χ^s_α(τ). From χ^s_α(τ) we can also extract the spin correlation length ξ_α along the imaginary-time axis,

    ξ_α=1/ν_1√(χ^s_α(ν_0)/χ^s_α(ν_1)  -1).

Here χ^s_α(ν_n) is the Fourier transform of χ^s_α(τ). This is in close analogy with extracting the
spatial
correlation length from the momentum dependence of the structure factor <cit.>. This can be understood by considering an ansatz χ^s_α(ν_n)∝ (ν_n^2 + E^*^2 )^-x/2. At criticality, the crossover energy scale E^* vanishes and χ^s_α(ν_n) will diverge as χ^s_α(ν_n)∝ν_n^-x. Away from criticality E^* is finite and will contribute a factor e^-E^*τ to χ_α^s(τ) when transformed back to imaginary-time domain. Then using equation (<ref>) the crossover scale E^* is inversely proportional to the correlation length ξ_α∝ 1/E^*.

As we will always preserve spin SU(2) symmetry, in the following we will drop the subscript α labeling different spin components in any vector quantity.

We will also look at the Binder cumulant <cit.>, generalized to a n-components order parameter <cit.>,

    U_2=n+2/2( 1-n/n+2⟨ (  m· m )^2⟩/⟨ m· m⟩^2),

which is defined such that U_2 approaches 1 in the ordered state and 0 in the disordered state. Quantities like ⟨ (  m· m )^2⟩ will involve 4-point correlation functions of different components of S_α which would require implementing worm type algorithm <cit.>. 
In the presence of spin SU(2) symmetry, we can utilize the relation ⟨ (  m· m )^2⟩=5 ⟨ m_z^4⟩ and ⟨ m· m⟩ = 3⟨ m_z^2⟩ to simplify the expression (for n=3),

    U_2=5/2( 1- 1/3⟨ m_z^4⟩/⟨ m_z^2⟩^2).


Another interesting quantity is the fidelity susceptibility χ^λ_f. Suppose the Hamiltonian is composed of two parts H=H_λ=0+λ H_λ, with λ being some tuning parameter. Then χ^λ_f is defined as <cit.>

    χ^λ_f(T)=∫_0^β/2( ⟨ T_τ H_λ(τ) H_λ⟩ - ⟨ H_λ⟩^2) τ   dτ.


At a second order quantum phase transition,
⟨ :H_λ(τ) : : H_λ : ⟩∼( 1/τ) ^2 Dim[H_λ].
Here :H_λ: denotes normal ordering :H_λ:=H_λ-⟨ H_λ⟩ and Dim[H_λ] denotes scaling dimension of H_λ. As we require ∫_0^β dτλ H_λ is scale invariant, we have Dim[H_λ]=1-Dim[λ], so 
⟨ :H_λ(τ) : : H_λ : ⟩τ∼( 1/τ) ^1-2 Dim[λ]. We see that if λ is relevant at the critical point, in which case it is usually identified as the correlation length exponent ν^-1, Dim[λ]=ν^-1>0, then χ^λ_f(T) will diverge,

    χ^λ_f(T) ∝β^2/ν.


Therefore χ^λ_f can be used to detect the location of a QPT, without knowing the actual order parameter. It turns out for hybridization expansion CT-QMC, if we choose λ to be the hybridization strength V, then the corresponding fidelity susceptibility, which we denoted by χ^V_f, can be calculate by a very simple formula <cit.>,

    χ^V_f=⟨ k_L k_R⟩ - ⟨ k_L⟩⟨ k_R⟩/ 2 V^2,

where we have considered dividing the imaginary-time axis into two pieces, and k_L and k_R are the number of H_λ between [0,β/2] and [β/2,β] at each Monte Carlo step, respectively. 





§ QUANTUM PHASE TRANSITIONS AND PHASE DIAGRAM



We now present the CT-QMC results. We describe the details of our analysis in the representative cases
of s=0.6 in section <ref> and s=0.2 in section <ref>. We then consider the dependence
on s in the range 0<s<1 appropriate for sub-ohmic bosonic bath in section <ref>.



 §.§ s=0.6


We start by presenting our analysis at s=0.6, which belongs to the case of RG flow specified in figure <ref> (a).
Alongside with C^' that controls the transition from local moment phase to Kondo phase, 
due to the appearance of a stable fixed point L representing the critical phase, 
we have two additional unstable fixed points C and LC, 
each describing the transition from critical phase to Kondo phase and critical phase to local moment phase, 
respectively. In the following, we will present numerical evidence for each of the three QCPs.



  §.§.§ Critical phase-Kondo transition


First we stay at g=0.5, and gradually increase Γ_0.
In figure <ref>(a) we plot ξ/β versus Γ_0 from β=200 all the way to β=6400. For Γ_0≲ 0.08, we find ξ/β is almost independent of β (system size), suggesting the system being scale invariant for a range of Γ_0. This is the signature of the critical phase. At larger Γ_0, ξ grows slower than the system size β, signifying short time correlation between the impurity spin as the impurity is Kondo screened. At some critical value of Γ_0 we expect a quantum phase transition separating the two phases. But the exact location is hard to pin-point, as we do not see any crossing in ξ/β. In section <ref> we will show that ξ/β does have a crossing at the local moment to Kondo QCP.






One observable we can utilize is the root mean square magnetization σ defined in equation (<ref>). We expect a scaling form as follows should hold,

    σ(Γ_0,β) = β^-(1-x)/2σ̃(β^1/ν(Γ_0-Γ_c)/Γ_c + A/β^ϕ/ν).

where σ̃ is the universal function and A/ β^ϕ/ν is the sub-leading terms.

In the universal function σ̃ the dependence of the tuning parameter only comes in through the combination of β^1/ν (Γ_0-Γ_c) (ignoring sub-leading corrections). This can be justified from RG or understood phenomenologically bases on the consideration that at QCP the system only depend on the ratio β/ξ and ξ diverges with ξ∝ |Γ_0-Γ_c|^-ν. One subtlety here is that the correlation length diverges in the entire critical phase. So one could question whether such a scaling form still apply in the region of Γ_0<Γ_c.  The prefactor β^-(1-x)/2 comes from equation (<ref>) and that at the QCP we expect χ^spin∝β^x with the exact relation x=s based on ϵ-expansion RG result. Here instead of imposing this relation we allow x to adjust freely. As shown in figure <ref>(b), the quality of the scaling collapse suggests that equation (<ref>) is the correct scaling hypothesis. In addition the correspondingly determined Γ_c=0.08(1) and ν^-1=0.26(3) are consistent with what we obtained from χ^V_f. We also find x=0.63(4), consistent with the prediction x=s. 

From ϵ-expansion calculation to second order <cit.>, we obtain ν^-1= ϵ/2 + ϵ^2 /6 ≃ 0.23, in reasonably good agreement with the numerical value.
 
Unlike the χ^spin(T) ∼  1/T local moment behavior in the s=0.2 case previously found in reference <cit.>, here the temperature dependence of the spin susceptibility obeys a nontrivial power-law, as shown in figure <ref>. We find x=0.66,0.67,0.66,0.65 for Γ_0=0.04,0.05,0.06,0.07 respectively. We interpret this as all the Γ_0 < Γ_c points under RG will flow towards the critical phase fixed point L with χ^spin(T) ∼ A_1/T^s. Notice that according to ϵ-expansion the leading irrelevant operator has a very small scaling dimension y_i=-ϵ/2+O(ϵ^2), so the deviation from the exact relation x=s is most likely due to corrections to scaling. At Γ_0=Γ_c, we have x=0.61. This is also consistent with the predicted critical behavior χ^spin(T) ∼ A_2/T^s at fixed point C from ϵ-expansion RG.







  §.§.§ Critical phase-local moment transition

So far we have considered the regime accessible by the ϵ-expansion of the SU(2) model, 
namely when both the fermionic and bosonic couplings are small.
Unlike the Coulomb-gas expansion of the Ising case <cit.>,
the ϵ-expansion here does not reach the regime of large g.
In order to simplify the calculation we set Γ_0=0 in this section. 
We have also performed calculation at small but
nonzero Γ_0 and the conclusion remains the same.

First let us look at the behavior of the correlation length as a function of g, plotted in figure <ref>. The low temperature behavior of ξ/β for g≲ 0.5 resembles the critical phase behavior in figure <ref>(a), both converging to a value around 0.3. For g≳ 0.8, on the other hand, ξ/β rises as temperature decreases, which suggests local moment phase behavior. 




A more quantitative way of studying the transition between these two phases is by looking at the temperature dependence of the mean square magnetization σ^2. 
Following reference <cit.>, the low temperature behavior of χ^spin(T) can be described by the following ansatz,

    χ^spin(T)=M_0/T + 1/T^x T_B^1-x.

Here M_0 is the Curie constant, T_B the crossover temperature scale above which the critical fluctuation part T^-x will dominate. This together with equation (<ref>) leads to

    σ^2(T)=M_0 + ( T/T_B)^1-x.


Our result for σ^2(T) is plotted in figure <ref>. 
For g≤ 0.5, the data can be described by equation (<ref>) with M_0=0 and x=0.68,0.67,0.66 for g=0.3,0.4,0.5. This is the critical phase and the exponent is very close to what we obtained at Sec.<ref>. For g ≥ 0.8, fitting σ^2(T) using the same equation gives a finite M_0. This indicates we are entering the local moment phase. While we have obtained x=0.60 for g>1, we have x=0.67,0.65,0.64 for g=0.8, 0.9, 1.0, reflecting corrections to scaling not captured by equation (<ref>).

The extracted M_0 is plotted in figure <ref>(b). Close to the transition point at g=g_c, we expect M_0 will vanish as M_0∝ (g-g_c)^β_1. We attempt to use this relation to find the value of g_c by fitting over the M_0 versus g data. Bearing in mind that for 0.8≤ g ≤ 1 the value of x obtained from equation (<ref>) is larger than s, it is likely that we will be overestimating M_0 in this region, so we only use M_0 down to g≥ g_min, and vary g_min from 0.8 to 1. Depending on the cutoff g_min, the obtained g_c lands within the range g_c∈ [0.74,0.88]. Notice that the fitting with different g_min all describe the g ≥ 1 part of the data quite well. We thus take our final estimate of g_c to be g_c = 0.8 ± 0.1.







  §.§.§ Local moment-Kondo transition 


Now that we have established that the system resides in the local moment phase for g>g_c≃ 0.8 at Γ_0=0, we consider a path to the Kondo screened phase by turning on the hybridization while fixing g=1. As expected, we observe a crossing in ξ/β, and a divergence in χ_f^V, both around Γ_0=0.4 (cf. figure <ref>). 





Similar to what we have done for σ(Γ_0,β) in equation (<ref>), we consider the following finite size scaling hypothesis for ξ and χ_f^V,

    ξ(Γ_0,β)   =   βξ̃(β^1/ν(Γ_0-Γ_c)/Γ_c+A/β^ϕ/ν) , 
    χ^V_f(Γ_0,β)    =   β^2/νχ̃(β^1/ν(Γ_0-Γ_c)/Γ_c+A/β^ϕ/ν).
    .



As seen in figure <ref>, close to the critical point the data fall nicely under a single universal curve. We obtain Γ_c=0.35(2), ν^-1=0.39(6) from ξ and Γ_c=0.34(2), ν^-1=0.37(5) from χ_f^V. Our final estimated value are Γ_c=0.35(2) and ν^-1=0.38(5). The value of ν^-1 obtained here for critical point C^' is in sharp contrast with that for critical point C with ν=0.25(4). This further established that C and C^' are two distinct critical points.








We now turn to the critical behavior of spin susceptibility. We plot χ^spin vs. T at different Γ_0. At the critical coupling Γ_0=Γ_c≃ 0.35, χ^spin(T) can be fitted with a power law χ^spin(T) ∝ T^-x with x=0.65. Inside the local moment phase at Γ_0=0.1, it can be described by equation (<ref>) with a finite M_0=0.10 for the M_0/T term and a sub-leading 1/T^xT_B^1-x term with x=0.62. These are consistent with the critical spin fluctuations being dominated by a T^-s behavior. Thus we think the local spin susceptibility at C^' should also diverge as χ^spin∼ 1/T^s.







 §.§ s=0.2


We now turn to the s=0.2 case. 
This is also the case investigated in reference <cit.> at the U=∞ limit. We will fix g=0.5 and gradually increase Γ_0 to find the QCP from local moment phase to Kondo screened phase.

We first plot the dependence on Γ_0 of the Binder cumulant U_2 and the reduced correlation length ξ/β in figure <ref>(a) and figure <ref>(b), where we have identified crossing points for both quantities. 
This suggests a transition from a local moment phase at small Γ_0 to a Kondo screened phase at large Γ_0. 
The crossing points have a sizable drift as we lower the temperature, which can be seen more clearly by plotting the crossing points between curves at β and 2β in figure <ref>. 
We see that the crossing points obtained from U_2 and ξ/β are approaching to the same critical value Γ_c in the T=0 limit from opposite directions. By extrapolating the crossing points Γ_cross to T=0 using a simple power-law relation Γ_cross=Γ_c+aT^b, we find Γ_c=0.48(1).









We can then repeat the analysis done in Sec.<ref> for the same type of transition at s=0.2 by considering scaling collapse of the form in equation (<ref>) for the correlation length ξ and similarly for the Binder cumulant U_2, 

    U_2(Γ_0,β)= Ũ_2(β^1/ν(Γ_0-Γ_c)/Γ_c+A/β^ϕ/ν)

where the presence of the sub-leading term A/β^ϕ/ν can take into account the finite temperature shift of the crossing point.

It turns out these ansatzes describe the data very well. The collapsed data using equation (<ref>) and equation (<ref>) is plotted in figure <ref>(a), and they give consistent estimates for the value of the critical coupling Γ_c and correlation length exponent ν. We obtain Γ_c=0.49(1), ν^-1=0.42(3) from U_2 and Γ_c=0.48(1), ν^-1=0.43(3) from ξ. 

We further test the applicability of the fidelity susceptibility in this case, which serves as another independent tool to detect the QCP. As shown in figure <ref>(c) the measured χ^V_f appears to diverge near our estimated Γ_c. A finite size scaling analysis can be performed as well. For consistency we consider the same type of scaling form of χ^V_f as appeared in equation (<ref>),





The result, plotted in figure <ref>(b), gives Γ_c=0.46(2) and ν^-1=0.48(3), in fairly good agreement with what we have obtained from U_2 and ξ. Our final estimates are Γ_c=0.48(1) and ν^-1=0.44(5).

After identified the location of the QCP, we now look at the temperature dependence of the spin susceptibility χ^spin across the QCP, shown in figure <ref>(a).  It turns out the critical behavior of χ^spin is much harder to study for the s=0.2 case compared to the s=0.6 case. For Γ_0<Γ_c, the dominant behavior of χ^spin(T) is Curie-Weiss like, reflecting the localized nature of the impurity spin. For Γ_0>Γ_c, χ^spin(T) will saturate at low T, corresponding to Kondo singlet formation. In between, we can see some indication of quantum critical behavior χ^spin(T) ∝ T^-s at Γ_0=0.50, slightly away from our estimated Γ_c. We think this is due to the fact that χ^spin(T) at Γ_0=Γ_c is still in the initial cross-over regime. To see this, we may define a transient power law exponent by α(T)=-d log(χ^spin(T) )/ d log(T). For Γ_0≤ 0.46 we find α(T) is increasing as T is lowered while for Γ_0≥ 0.48 it is decreasing. 

We note that the calculation in reference <cit.> has assumed the relation χ^spin(T) ∝ T^-s and use it as a tool to locate the QCP by looking for the crossing point of T^sχ^spin(T) at different T. But there the crossing point has significant drift versus temperature, which is consistent with an evolving α(T) in our calculation. Here we determine the critical coupling Γ_c via a variety of independent methods and obtained unambiguous results for the presence and the location of the QCP. Then we attempt to verify the critical behavior of χ^spin(T) directly. Unfortunately from figure <ref>(b) it seems that, in contrast to the case of s=0.6, accessing the asymptotic
critical regime requires even lower temperatures for s=0.2:  it appears we 
would need at least two more decades below the lowest temperature we can obtain T=1/6400 
in the temperature dependence of χ^spin(T) to access the true critical behavior for s=0.2.
We interpret this as the s=0.2 case has a extremely low entry point T<10^-5 to reach the asymptotic 
critical regime for χ^spin(T), even though the bare Kondo temperature T_K^0 is about 
T_K^0=2.5 at Γ=Γ_c in the absence of bosonic coupling. 
We have seen earlier that in the s=0.6 case it is much easier to access the asymptotic critical behavior 
of χ^spin(T).









 §.§ Phase diagram upon varying the power of the sub-ohmic spectrum


The phase diagram, as specified by the two types of RG flows given in figure <ref>, can be
determined for any given 0<s<1 once we have estimated s^*. For this purpose,
we can turn to the pure bosonic problem by setting Γ_0=0, and vary both the bosonic coupling g 
as well as the bosonic bath exponent s. 
As s → 0^+, the procedure to obtain J(τ) defined in equation (<ref>) 
and equation (<ref>) will encounter convergence issue. 
As the critical property only depends on the long time asymptotic behavior of J(τ), we directly adopt 
a J(τ) that has the correct 1/τ^1+s dependence as our input without specifying the actual form 
of ρ_b(ω). To be specific, we choose J(τ) to be the following,

    J(τ)=[  π/β/sin(πτ/β)  (1+e^-β-e^-τ-e^-(β-τ)) ]^1+s.


The exponential factor will make J(τ) finite at the end points: lim_τ→ 0 J(τ) = lim_τ→β J(τ) = 1. Also J(τ) is even under reflection about τ=β/2.

We can then integrate J(τ) twice to get B(τ),

    B(τ)-B(0)=∫_0^τ∫_0^τ^' J(τ^'') d τ^' d τ^'' + aτ.

with a=-∫_0^β/2 J(τ^'') d τ^'' determined from the condition d B(τ) /d τ |_τ=β/2 =0.





Using equation (<ref>) and equation (<ref>) as input we have obtained the dynamical spin correlation function χ^s(τ) for different value of g and s. In figure <ref>(a) we present the result of χ^s(τ) vs. τ at several different value of β for the specific case of g=0.4, s=0.2. At each β, χ^s(τ) drops from 1/4 at τ=0 and reaches its minimum at τ=β/2. As β is increased, χ^s(τ=β/2) converge to a finite value. 







We then plot the evolution of χ^s(τ=β/2) obtained at low temperature, as a function of s for four different choices of g in figure <ref>(b), up to the smallest value of  χ^s(τ=β/2) that we can reach convergence. We can identify χ^s(τ=β/2) obtained here as an effective Curie constant, and use it as the order parameter for the local moment phase. We see that for fixed g, χ^s(τ=β/2) decreases smoothly as a function of s. Furthermore, we can extrapolate each curve to larger value of s until χ^s(τ=β/2) vanishes at some critical value s=s_c(g). This gives the value of s where the corresponding g is the critical value between the local moment phase and the critical phase.

The dependence of s_c( g ) on g  maps out the phase boundary between the local moment phase and the critical phase, which is shown in figure <ref>. Note that the shape of the phase boundary will depend on the specific form of J(τ) that is employed.
 As we can see the dependence of s_c (g) on the value of s is fairly weak and it reaches the g^2=0 axis at around s=s^*≃ 0.47.  Although we do need to admit the simple extrapolation scheme performed in figure <ref>(a) could introduce some error. 



§ DISCUSSION


Our result has important implications for the Kondo lattice model. In the EDMFT solution of the Kondo lattice model, the Kondo destruction QCP of the lattice problem is embedded in the impurity QCP of an effective BFKM. The self-consistency condition is satisfied at s→ 0^+, or ϵ→ 1^-, provided the relation η=ϵ holds, which initially is an statement made at critical point C from ϵ-expansion perspective. Our calculation implies that C should disappear before ϵ reaches 1, and that the actual impurity QCP encountered in the EDFMT calculation should be C^' instead. Nonetheless, the relation η=ϵ is still true at C^', even though C and C^' have different correlation length exponent, thus belonging to different universality classes. This is quite surprising until we realize that the argument that leads to η=ϵ only relies on the condition η=ϵ + 2 β(g) /g |_g=g^*,J=J^*, which is shown to be valid to all orders in ϵ in reference<cit.>. The relation η=ϵ then follows at any intermediate coupling fixed point g=g^*, J=J^*, where β(g) /g |_g=g^*,J=J^*=0, regardless of whether g^* and J=J^* is of the order ϵ. Thereby this argument can be extended to C^' as well.



§ CONCLUSIONS


We have studied the SU(2) Bose-Fermi Anderson model using CT-QMC, focusing on the Kondo destruction type QCP. We find two type of such QCPs: one from Kondo screened phase to a local moment phase, the other to a critical phase. The second type QCP only exists when s>s^*, in which case the critical properties we have calculated 
agree with those from an ϵ-expansion RG. At both types of QCP, our results suggest the 
spin correlation function obeys the power law χ^spin(τ) ∼ (1/τ)^η with η=1-s.



§ ACKNOWLEDGMENT

We thank H. Hu, J. Pixley and K. Ingersent for useful discussions.
The work was in part supported by the NSF (DMR-1611392), the Robert A. Welch Foundation (C-1411), 
Computing time was supported in part by the Data Analysis and Visualization Cyberinfrastructure funded by NSF under grant OCI-0959097 and an IBM Shared University Research (SUR) Award at Rice University, and by the Extreme Science and Engineering
Discovery Environment (XSEDE) by NSF under Grants No. DMR170109.
Q.S. acknowledges the hospitality  and support by a Ulam Scholarship from the Center for Nonlinear Studies at Los Alamos National Laboratory, and the hospitality of the Aspen Center for Physics, which is supported by NSF grant No. PHY-1607611.











