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

TL;DR

This study develops a specialized continuous-time Quantum Monte Carlo method to analyze the phase diagram and critical points of the Bose-Fermi Anderson model with SU(2) symmetry, revealing two types of Kondo destruction quantum critical points depending on the bosonic bath spectrum.

Contribution

The paper introduces a low-temperature, SU(2)-symmetric CT-QMC method for the Bose-Fermi Anderson model, clarifying its phase diagram and critical points with new insights into Kondo destruction.

Findings

Identified two types of Kondo destruction quantum critical points depending on the bosonic bath exponent.

Confirmed the relation η=ε at the quantum critical points, consistent with analytical predictions.

Validated the CT-QMC results against renormalization-group calculations for certain parameter regimes.

Abstract

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 $\rho_{b}(\omega) \propto \omega^{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 $\epsilon$-expansion renormalization-group calculation (here $\epsilon=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 $\eta$ obeys the relation $\eta=\epsilon$, which has important implications for Kondo destruction QCP in the Kondo lattice problem.

Equations

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

$$H_c=\sum _{k,\sigma }{ϵ}_k{c}_{k,\sigma }^{†}{c}_{k,\sigma },H_b=\sum _{\alpha }{H}_b^{\alpha }=\sum _{p,\alpha }{\omega }_q{{\varphi }_p^{\alpha }}^{†}{\varphi }_p^{\alpha }.$$

$$H_d=\sum _{\sigma }{ϵ}_d{d}_{\sigma }^{†}{d}_{\sigma }+U{d}_{\uparrow }^{†}{d}_{\uparrow }{d}_{\downarrow }^{†}{d}_{\downarrow }.$$

$$H_V=\sum _{k,\sigma }V{d}_{\sigma }^{†}{c}_{k,\sigma }+h.c.,H_g=\sum _{p,\alpha }g{S}_{\alpha }({{\varphi }_p^{\alpha }}^{†}+{\varphi }_{-p}^{\alpha })$$

$${\rho }_F(ϵ)=\sum _k\delta (ϵ-{ϵ}_k)={\rho }_0\Theta (|D-ϵ|),$$

$${\rho }_b(\omega )=\sum _q\delta (\omega -{\omega }_q)=K_0{\omega }^s{e}^{-\omega /\Lambda }\Theta (\omega ).$$

$$\tilde{H}$$

$${\tilde{H}}_d$$

$${\tilde{H}}_V$$

$${\tilde{H}}_g$$

$$Z$$

$$w_d$$

$$w_c^{\sigma }$$

$$F_{ij}^{\sigma }=\frac{-{\sum }_k{V}^{2}Tr[e^{-\beta H_c}{T}_{\tau }{c}_{k,\sigma }({\tau }_j^{d\sigma })c_{k,\sigma }^{†}({\tau }_i^{d^{\prime }\sigma })]}{Tr[e^{-\beta H_c}]}.$$

$$w_z$$

$$B({\tau }_j-{\tau }_i)$$

$$w_p$$

$$P_{ij}$$

$$P_{prop}=\frac{1}{l_{max}m}\times \frac{J({\tau }_k^{d\uparrow }-{\tau }_j^{s^{\prime }})}{{\sum }_{n=1}^{n_{\uparrow }}J({\tau }_n^d-{\tau }_j^{s^{\prime }})}.$$

$$P_{prop}^{inv}=\frac{1}{l_{max}^{\prime }m}\times \frac{J({\tau }_i^s-{\tau }_j^{s^{\prime }})}{{\sum }_{\begin{array}{c}n=1,\\ n\ne k\end{array}}^{n_{\uparrow }}J({\tau }_n^d-{\tau }_j^{s^{\prime }})+J({\tau }_i^s-{\tau }_j^{s^{\prime }})}.$$

$$\frac{w_{new}}{w_{old}}$$

$$R=\frac{w_{new}}{w_{old}}\times \frac{P_{prop}^{inv}}{P_{prop}}.$$

$$J(i{\nu }_n)=\sum _p\frac{2{\omega }_p}{{{\omega }_p}^2-{\nu }_n^2},$$

$$J(i{\nu }_n)={\int }_0^{\infty }\frac{2\omega }{{\omega }^2-{\nu }_n^2}{\rho }_b(\omega )d\omega ,$$

$$J(\tau )=\frac{1}{\beta }J(i{\nu }_n=0)+\frac{2}{\beta }\sum _{{\nu }_n>0}J(i{\nu }_n)\cos ({\nu }_n\tau ).$$

$$J(\tau )={\int }_0^{\infty }\frac{e^{(\beta -\tau )\omega }+e^{\tau \omega }}{e^{\beta \omega }-1}{\rho }_b(\omega )d\omega .$$

$$B(\tau )-B(0)$$

$$⟨m_{\alpha }⟩=⟨\frac{1}{\beta }{\int }_0^{\beta }{S}_{\alpha }(\tau )d\tau ⟩,\alpha =x,y,z,$$

$${\sigma }_{\alpha }=\sqrt{⟨m_{\alpha }^2⟩},$$

$${\chi }_{\alpha }^{spin}=\beta {\sigma }_{\alpha }^2,$$

$${\xi }_{\alpha }=\frac{1}{{\nu }_1}\sqrt{\frac{{\chi }_{\alpha }^s({\nu }_0)}{{\chi }_{\alpha }^s({\nu }_1)}-1}.$$

$$U_2=\frac{n+2}{2}\left(1-\frac{n}{n+2}\frac{⟨(\mathbf{m}\cdot \mathbf{m}{)}^2⟩}{⟨\mathbf{m}\cdot \mathbf{m}{⟩}^2}\right),$$

Full text

I Introduction

Heavy fermion systems serve as a prototype system to study quantum criticality Si and Steglich (2010); Coleman and Schofield (2005). 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 Si et al. (2001); Coleman et al. (2001); Senthil et al. (2004), 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 $\omega/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 Si and Paschen (2013). Such properties are inconsistent with predictions from the traditional spin-density-wave type QCP Hertz (1976); Millis (1993); Moriya (2012).

One of the simplest models that contain a Kondo destruction type QCP is the Bose-Fermi Kondo model (BFKM) Si et al. (2014). 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) Si et al. (2001, 2003). 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 ($\omega$) have a power-law form, $\rho_{b}(\omega)\propto\omega^{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 $\epsilon$-expansion renormalization group (RG) method, using $\epsilon=1-s$ as a small parameter Si and Smith (1996); Smith and Si (1999); Sengupta (2000); Si et al. (2001, 2003); Zhu and Si (2002); Zaránd and Demler (2002). 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 $\epsilon$-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 $\chi_{spin}(\tau)\sim(1/\tau)^{\eta}$, with an exact relation $\eta=\epsilon$ Zhu and Si (2002); Zaránd and Demler (2002). 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 $\eta=\epsilon$ will remain valid at $\epsilon\rightarrow 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) Glossop and Ingersent (2005); Glossop and Ingersent (2007a) and continuous-time quantum Monte Carlo (CT-QMC) Pixley et al. (2011, 2013); Otsuki (2013), or as an effective model under EDMFT Grempel and Si (2003); Zhu et al. (2003); Glossop and Ingersent (2007b); Zhu et al. (2007). Our focus in this work is on the CT-QMC method, from which a seeming controversy existed for the SU(2) symmetric BFAM Otsuki (2013): 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 $\chi^{spin}\sim 1/T$ instead of the $\chi^{spin}\sim 1/T^{s}$ behavior predicted by $\epsilon$-expansion RGZhu and Si (2002); Zaránd and Demler (2002) 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 $\epsilon$-expansion RG in the range of coupling constants accessed by this expansion (again $\epsilon=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 [math] 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 $\epsilon$-expansion Zhu and Si (2002); Zaránd and Demler (2002) and CT-QMC results are fully compatible with each other. Our results are summarized by the RG-flow diagrams of figure 1. For the $s>s^{*}$ regime, we identify i) the critical point C separating the Kondo screened phase and critical phase, as predicted from $\epsilon$-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 Otsuki (2013). We also determine the correlation length exponent $\nu$. 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 $\epsilon$-expansion.

The remainder of the paper is organized as follows. In Sec. II 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. III. We will start with a detailed study for the $s=0.6$ case in Sec. III.1, followed by the $s=0.2$ case in Sec. III.2, before carrying through the analysis that leads to an estimate for the value of $s^{*}$ in Sec. III.3. We will discus the implication of our results in Sec. IV and conclude the article in Sec. V.

II Model and Method

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

[TABLE]

where $H_{c}$ and $H_{b}$ describes the bosonic and fermionic bath part, respectively,

[TABLE]

$H_{d}$ contains the local electron part,

[TABLE]

$H_{V}$ and $H_{g}$ couples the local orbital to the bosonic and fermionic bath,

[TABLE]

where the summation over $\alpha$ runs through $x$,$y$,$z$, $S_{\alpha}=d^{\dagger}_{\sigma}\tau^{\alpha}_{\sigma\sigma^{\prime}}d_{\sigma^{\prime}}$, and $\tau^{\alpha}_{\sigma\sigma^{\prime}}$ 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,

[TABLE]

which leads to a hybridization function $\Gamma(\epsilon)=\Gamma_{0}\Theta(|D-\epsilon|)$, with $\Gamma_{0}=\pi\rho_{0}V^{2}$.

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

[TABLE]

Throughout the text we fix $D=1$, $\Lambda=1$, and stays at the particle-hole symmetric point $U=-2\epsilon_{d}=0.1$. The prefactor $\rho_{0}$ and $K_{0}$ in the density of states of the fermionic bath and bosonic bath are determined from the normalization condition $\int_{-D}^{D}\rho_{F}(\epsilon)d\epsilon=1$ and $\int_{0}^{\infty}\rho_{b}(\omega)d\omega=1$. We will use either the amplitude of the hybridization function $\Gamma_{0}$ or the spin-boson coupling $g$ as our tuning parameter.

II.1 Monte-Carlo procedure

We will employ the CT-QMC algorithm, first introduced in reference Werner et al. (2006); Werner and Millis (2006) and then generalized to treat the BFAM in references Pixley et al. (2011, 2013); Otsuki (2013). We start with removing the $z$ component of the spin-boson coupling by employing a Firsov-Lang transformation $\tilde{H}=e^{S}He^{-S}$ with $S=gS_{z}\sum_{p}\frac{1}{\omega_{p}}({\phi_{p}^{z}}^{\dagger}+\phi_{-p}^{z})$ (similar to Ref. Werner and Millis, 2007) and work with the transformed Hamiltonian $\tilde{H}$,

[TABLE]

where we have defined the renormalized parameters $\tilde{\epsilon}_{d}=\epsilon_{d}-(g^{2}/4)\sum_{q}(1/\omega_{q})^{2}$, $\tilde{U}=U+(g^{2}/2)\sum_{q}(1/\omega_{q})^{2}$, $s_{\sigma}=\pm 1/2$ for $\sigma=\uparrow/\downarrow$. and recombined the $x$ and $y$ components of $S_{\alpha}$ and $\phi_{\alpha}$ into $S_{+}=d^{\dagger}_{\uparrow}d_{\downarrow}$, $S_{-}=d^{\dagger}_{\downarrow}d_{\uparrow}$, $\phi_{p}^{\pm}=(1/\sqrt{2})\left(({\phi_{p}^{x}}^{\dagger}+\phi_{p}^{x})\pm i({\phi_{p}^{y}}^{\dagger}+\phi_{p}^{y})\right)$. The partition function is constructed by expanding in the non-diagonal terms Werner et al. (2006); Werner and Millis (2006); Pixley et al. (2011, 2013); Otsuki (2013), $\tilde{H}_{V}$ and $\tilde{H}_{g}$ under the interaction representation of $H_{0}\equiv H_{b}+H_{c}+\tilde{H}_{d}$. It has the following form Pixley et al. (2011, 2013); Otsuki (2013):

[TABLE]

where $Z_{0}=Tr[e^{-\beta H_{c}}]Tr[e^{-\beta H^{z}_{B}}]Tr[e^{-\beta(H^{x}_{B}+H^{y}_{B})}]$ is the partition function of the bath, $\beta$ being the inverse temperature: $\beta=1/T$. $\int\prod_{i=1}^{m}d\tau^{\alpha}_{i}d\tau^{\alpha^{\prime}}_{i}=\int_{0}^{\beta}d\tau^{\alpha}_{1}\cdots\int_{\tau^{\alpha}_{N-1}}^{\beta}d\tau^{\alpha}_{N}\int_{0}^{\beta}d\tau^{\alpha^{\prime}}_{1}\cdots\int_{\tau^{\alpha^{\prime}}_{N-1}}^{\beta}d\tau^{\alpha^{\prime}}_{N}$. $\{\tau^{\alpha}\}_{n}$ denotes the set of imaginary time of all the operators of a given type $\alpha$ in the expansion: $\{\tau^{\alpha}\}_{n}=\{\tau^{\alpha}_{1},\tau^{\alpha}_{2}\dots,\tau^{\alpha}_{n}\}$. $\alpha\in\{s,s^{\prime},d\sigma,d\sigma^{\prime}\}$ represents $S_{+}$, $S_{-}$, $d^{\dagger}_{\sigma}$, or $d_{\sigma}$. $n=m$ or $n_{\sigma}$ denotes the number of pairs of $S_{+},S_{-}$ or $d^{\dagger}_{\sigma},d_{\sigma}$, also labeling the expansion order. $\{\tau^{tot}\}_{n_{tot}}$ refers to all the $\{\tau^{\alpha}\}_{n}$ combined, with $n_{tot}=2(\sum_{\sigma}n_{\sigma}+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}(\{\tau^{tot}\}_{n_{tot}})$ is the contribution from the local d electron part. It describes valence and spin fluctuations of the local orbitals,

[TABLE]

Here for a given operator $O$, $O(\tau)$ denotes the corresponding operator in the interaction representation $O(\tau)=e^{\tau H_{0}}Oe^{-\tau H_{0}}$.

$w^{\sigma}_{c}(\{\tau^{d\sigma}\}_{n_{\sigma}},\{\tau^{d^{\prime}\sigma}\}_{n_{\sigma}})$ is the contribution from the conduction electron with spin index $\sigma$,

[TABLE]

It can be expressed as a determinant of matrix $F^{\sigma}$, whose matrix element is given by

[TABLE]

$w_{z}(\{\tau^{tot}\}_{n_{tot}})$ comes from the $z$ component bosonic bath part Pixley et al. (2011, 2013),

[TABLE]

where $s_{i}=\pm s_{\sigma}$ or $\pm 1$ when the operator $O(\tau_{i}^{tot})$ at $\tau_{i}^{tot}$ corresponds to $d_{\sigma}^{\dagger}/d_{\sigma}$ or $S^{\pm}$, and

[TABLE]

Finally, $w_{p}(\{\tau^{s}\}_{m},\{\tau^{s^{\prime}}\}_{m})$ involves the bosonic bath in the transverse direction Otsuki (2013), forming a permanent,

[TABLE]

The summation extends over $S_{m}$, representing all permutations of $1,2,\cdots,m$. The matrix element of $P$ is the following,

[TABLE]

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 $\{\tau^{\alpha}\}_{n}$ and a particular permutation $p\in S_{m}$, which is then sampled through a Metropolis algorithm with a probability proportional to $w_{d}\times w_{z}\times w_{c}^{\uparrow}\times w_{c}^{\downarrow}\times\prod_{i=1}^{m}P_{i,p(i)}$.

We now describe the Monte Carlo updates. We inherit the updates from the Ising BFAM Pixley et al. (2011, 2013), namely the insertion, removal and shift of $d_{\sigma}^{\dagger}c_{k,\sigma}$/$c^{\dagger}_{k,\sigma}d_{\sigma}$ pair, and also adopt the insertion/removal of $S^{+}\phi^{-}$/$S^{-}\phi^{+}$ and the sampling of the permutation $S_{m}$ introduced in reference Otsuki (2013) (named updates (a)-(c) there). In addition we introduce a swap update that swaps $S_{+}(S_{-})$ with a pair of $d_{\uparrow}^{\dagger}$ and $d_{\downarrow}$ ( $d_{\downarrow}^{\dagger}$ and $d_{\uparrow}$ ). For example consider the $S_{+}$ case. We first randomly pick a pair of $S_{+}(\tau^{s}_{i})$, $S_{-}(\tau^{s^{\prime}}_{j})$ from the $m$ pairs of $S_{+}$ and $S_{-}$ that is connected by one of $J(\tau)$. Then we choose a $d_{\uparrow}^{\dagger}(\tau^{d\uparrow}_{k})$ with a probability $P_{k}=J(\tau^{d\uparrow}_{k}-\tau^{s^{\prime}}_{j})/\left(\sum_{n=1,n_{\uparrow}}J(\tau^{d\uparrow}_{n}-\tau^{s^{\prime}}_{j})\right)$ from the $n_{\uparrow}$ of $d^{\dagger}_{\uparrow}$ operators. We then swap the position of $S_{+}(\tau^{s}_{i})$ and $d_{\uparrow}^{\dagger}(\tau^{d\uparrow}_{k})$. Finally, we find the $d_{\downarrow}(\tau^{d^{\prime}\downarrow}_{l})$ that is closest to $d_{\uparrow}^{\dagger}(\tau^{d\uparrow}_{k})$ before the swap, and move it to $d_{\downarrow}(\tau^{d^{\prime}\downarrow}_{new})$. $\tau^{d^{\prime}\downarrow}_{new}$ is randomly selected within an interval of length $l_{max}$, which is the distance between two creation operators in the $\sigma=\downarrow$ orbital next to $S_{+}$ before the swap. The corresponding proposal probability is given by

[TABLE]

Likewise we can find the proposal probability for the inverse update,

[TABLE]

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

[TABLE]

where $\{\tau^{d\uparrow}\}^{new}_{n_{\uparrow}}$ is $\{\tau^{d\uparrow}\}_{n_{\uparrow}}$ with $\tau_{k}^{d\uparrow}$ replaced by $\tau_{i}^{s}$, $\{\tau^{d\downarrow}\}^{new}_{n_{\downarrow}}$ is $\{\tau^{d\downarrow}\}_{n_{\downarrow}}$ with $\tau_{l}^{d^{\prime}\downarrow}$ replaced by $\tau_{lnew}^{d^{\prime}\downarrow}$, and $\{\tau^{tot}\}^{new}_{n_{tot}}$ is $\{\tau^{tot}\}_{n_{tot}}$ with the above two substitutions, plus $\tau_{i}^{s}$ replaced by $\tau_{k}^{d\uparrow}$.

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

[TABLE]

The reason that we choose the proposal probability to be the form in equation (14) and equation (15) is to cancel out the $J(\tau^{d\uparrow}_{k}-\tau^{s^{\prime}}_{j})/J(\tau^{s}_{i}-\tau^{s^{\prime}}_{j})$ factor in the weight ratio in equation (16), such that the acceptance ratio $R$ is of order 1. Otherwise if we select $d^{\dagger}_{\uparrow}(\tau_{k}^{d\uparrow})$ using a uniform distribution from [math] to $\beta$, since $J(\tau)\sim 1/{\tau^{1+s}}$, on average $J(\tau^{d\uparrow}_{k}-\tau^{s^{\prime}}_{j})\sim 1/\beta^{s}$, while the average value of $J(\tau^{s}_{i}-\tau^{s^{\prime}}_{j})$ is $\beta$ independent, as a result $R$ will be suppressed by a factor of $1/\beta^{s}$. Similar ideas have been introduced in reference Steiner et al. (2015).

In practice we have tested that the swap update introduced here can replace the role of update (d) in reference Otsuki (2013), which breaks up one $S^{+}$ ($S^{-}$) into a pair of $d^{\dagger}_{\uparrow}$ and $d_{\downarrow}$ ( $d^{\dagger}_{\downarrow}$ and $d_{\uparrow}$ ) 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 $\beta$ as a power-law, the swap update has an acceptance ratio that does not depend on $\beta$. 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(\tau)$ and $B(\tau)$ in the numerical calculation. This is important because in the current expansion scheme the weight contribution from the bosonic bath in the transverse direction $\phi^{\pm}$ and in the $z$ direction $\phi^{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(\tau)$ and $J(\tau)$ to sufficiently high accuracy.

Starting with the Fourier components of $J(\tau)$ in the matsubara frequency domain,

[TABLE]

where $\nu_{n}=2\pi n/\beta,\ n\in\mathbb{Z}$ is the matsubara frequencies. There are two ways to calculate $J(\tau)$. We can either perform the integration over the density of states first,

[TABLE]

followed by the matsubara summation,

[TABLE]

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

[TABLE]

In practice we find the summation in equation (20) converges too slow when $\beta$ is large. So using equation (21) is recommended.

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

[TABLE]

Because of the extra $1/\nu_{n}^{2}$ factor here, the summation actually converges very quickly.

II.2 Observables

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

We start with the local magnetization,

[TABLE]

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

Because the sampling will preserve spin rotation symmetry, the actual measured $\langle m_{\alpha}\rangle$ is always 0. Instead we measure its root mean square,

[TABLE]

which is also related to the static spin susceptibility $\chi^{spin}_{\alpha}(T)=\int_{0}^{\beta}\chi^{s}_{\alpha}(\tau)d\tau=\int_{0}^{\beta}\langle T_{\tau}{S}_{\alpha}(\tau)S_{\alpha}\rangle d\tau$ by,

[TABLE]

where we have also defined the dynamical spin correlation function $\chi^{s}_{\alpha}(\tau)$. From $\chi^{s}_{\alpha}(\tau)$ we can also extract the spin correlation length $\xi_{\alpha}$ along the imaginary-time axis,

[TABLE]

Here $\chi^{s}_{\alpha}(\nu_{n})$ is the Fourier transform of $\chi^{s}_{\alpha}(\tau)$. This is in close analogy with extracting the spatial correlation length from the momentum dependence of the structure factor Sandvik et al. (2010). This can be understood by considering an ansatz $\chi^{s}_{\alpha}(\nu_{n})\propto(\nu_{n}^{2}+{E^{*}}^{2})^{-x/2}$. At criticality, the crossover energy scale $E^{*}$ vanishes and $\chi^{s}_{\alpha}(\nu_{n})$ will diverge as $\chi^{s}_{\alpha}(\nu_{n})\propto\nu_{n}^{-x}$. Away from criticality $E^{*}$ is finite and will contribute a factor $e^{-E^{*}\tau}$ to $\chi_{\alpha}^{s}(\tau)$ when transformed back to imaginary-time domain. Then using equation (26) the crossover scale $E^{*}$ is inversely proportional to the correlation length $\xi_{\alpha}\propto 1/E^{*}$.

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

We will also look at the Binder cumulant Binder (1981), generalized to a n-components order parameter Sandvik et al. (2010),

[TABLE]

which is defined such that $U_{2}$ approaches $1$ in the ordered state and [math] in the disordered state. Quantities like $\langle({\bf m}\cdot{\bf m})^{2}\rangle$ will involve 4-point correlation functions of different components of $S_{\alpha}$ which would require implementing worm type algorithm Gunacker et al. (2015, 2016). In the presence of spin SU(2) symmetry, we can utilize the relation $\langle({\bf m}\cdot{\bf m})^{2}\rangle=5\langle m_{z}^{4}\rangle$ and $\langle{\bf m}\cdot{\bf m}\rangle=3\langle m_{z}^{2}\rangle$ to simplify the expression (for n=3),

[TABLE]

Another interesting quantity is the fidelity susceptibility $\chi^{\lambda}_{f}$. Suppose the Hamiltonian is composed of two parts $H=H_{\lambda=0}+\lambda H_{\lambda}$, with $\lambda$ being some tuning parameter. Then $\chi^{\lambda}_{f}$ is defined as Albuquerque et al. (2010)

[TABLE]

At a second order quantum phase transition, $\langle:H_{\lambda}(\tau)::H_{\lambda}:\rangle\sim\left({1}/{\tau}\right)^{2Dim[H_{\lambda}]}$. Here $:H_{\lambda}:$ denotes normal ordering $:H_{\lambda}:=H_{\lambda}-\langle H_{\lambda}\rangle$ and $Dim[H_{\lambda}]$ denotes scaling dimension of $H_{\lambda}$. As we require $\int_{0}^{\beta}d\tau\lambda H_{\lambda}$ is scale invariant, we have $Dim[H_{\lambda}]=1-Dim[\lambda]$, so $\langle:H_{\lambda}(\tau)::H_{\lambda}:\rangle\tau\sim\left({1}/{\tau}\right)^{1-2Dim[\lambda]}$. We see that if $\lambda$ is relevant at the critical point, in which case it is usually identified as the correlation length exponent $\nu^{-1}$, $Dim[\lambda]=\nu^{-1}>0$, then $\chi^{\lambda}_{f}(T)$ will diverge,

[TABLE]

Therefore $\chi^{\lambda}_{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 $\lambda$ to be the hybridization strength $V$, then the corresponding fidelity susceptibility, which we denoted by $\chi^{V}_{f}$, can be calculate by a very simple formula Wang et al. (2015a, b),

[TABLE]

where we have considered dividing the imaginary-time axis into two pieces, and $k_{L}$ and $k_{R}$ are the number of $H_{\lambda}$ between $[0,\beta/2]$ and $[\beta/2,\beta]$ at each Monte Carlo step, respectively.

III 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 III.1 and $s=0.2$ in section III.2. We then consider the dependence on $s$ in the range $0<s<1$ appropriate for sub-ohmic bosonic bath in section III.3.

III.1 s=0.6

We start by presenting our analysis at $s=0.6$, which belongs to the case of RG flow specified in figure 1 (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.

III.1.1 Critical phase-Kondo transition

First we stay at $g=0.5$, and gradually increase $\Gamma_{0}$. In figure 3(a) we plot $\xi/\beta$ versus $\Gamma_{0}$ from $\beta=200$ all the way to $\beta=6400$. For $\Gamma_{0}\lesssim 0.08$, we find $\xi/\beta$ is almost independent of $\beta$ (system size), suggesting the system being scale invariant for a range of $\Gamma_{0}$. This is the signature of the critical phase. At larger $\Gamma_{0}$, $\xi$ grows slower than the system size $\beta$, signifying short time correlation between the impurity spin as the impurity is Kondo screened. At some critical value of $\Gamma_{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 $\xi/\beta$. In section III.1.3 we will show that $\xi/\beta$ does have a crossing at the local moment to Kondo QCP.

One observable we can utilize is the root mean square magnetization $\sigma$ defined in equation (24). We expect a scaling form as follows should hold,

[TABLE]

where $\tilde{\sigma}$ is the universal function and $A/\beta^{\phi/\nu}$ is the sub-leading terms.

In the universal function $\tilde{\sigma}$ the dependence of the tuning parameter only comes in through the combination of $\beta^{1/\nu}(\Gamma_{0}-\Gamma_{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 $\beta/\xi$ and $\xi$ diverges with $\xi\propto|\Gamma_{0}-\Gamma_{c}|^{-\nu}$. 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 $\Gamma_{0}<\Gamma_{c}$. The prefactor $\beta^{-(1-x)/2}$ comes from equation (25) and that at the QCP we expect $\chi^{spin}\propto\beta^{x}$ with the exact relation $x=s$ based on $\epsilon$-expansion RG result. Here instead of imposing this relation we allow $x$ to adjust freely. As shown in figure 3(b), the quality of the scaling collapse suggests that equation (32) is the correct scaling hypothesis. In addition the correspondingly determined $\Gamma_{c}=0.08(1)$ and $\nu^{-1}=0.26(3)$ are consistent with what we obtained from $\chi^{V}_{f}$. We also find $x=0.63(4)$, consistent with the prediction $x=s$.

From $\epsilon$-expansion calculation to second order Zhu and Si (2002); Zaránd and Demler (2002), we obtain $\nu^{-1}=\epsilon/2+\epsilon^{2}/6\simeq 0.23$, in reasonably good agreement with the numerical value.

Unlike the $\chi^{spin}(T)\sim 1/T$ local moment behavior in the $s=0.2$ case previously found in reference Otsuki (2013), here the temperature dependence of the spin susceptibility obeys a nontrivial power-law, as shown in figure 4. We find $x=0.66,0.67,0.66,0.65$ for $\Gamma_{0}=0.04,0.05,0.06,0.07$ respectively. We interpret this as all the $\Gamma_{0}<\Gamma_{c}$ points under RG will flow towards the critical phase fixed point L with $\chi^{spin}(T)\sim A_{1}/T^{s}$. Notice that according to $\epsilon$-expansion the leading irrelevant operator has a very small scaling dimension $y_{i}=-\epsilon/2+O(\epsilon^{2})$, so the deviation from the exact relation $x=s$ is most likely due to corrections to scaling. At $\Gamma_{0}=\Gamma_{c}$, we have $x=0.61$. This is also consistent with the predicted critical behavior $\chi^{spin}(T)\sim A_{2}/T^{s}$ at fixed point C from $\epsilon$-expansion RG.

III.1.2 Critical phase-local moment transition

So far we have considered the regime accessible by the $\epsilon$-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 Si and Smith (1996); Smith and Si (1999); Zhu and Si (2002); Zaránd and Demler (2002), the $\epsilon$-expansion here does not reach the regime of large $g$. In order to simplify the calculation we set $\Gamma_{0}=0$ in this section. We have also performed calculation at small but nonzero $\Gamma_{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 5. The low temperature behavior of $\xi/\beta$ for $g\lesssim 0.5$ resembles the critical phase behavior in figure 3(a), both converging to a value around $0.3$. For $g\gtrsim 0.8$, on the other hand, $\xi/\beta$ 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 $\sigma^{2}$. Following reference Otsuki (2013), the low temperature behavior of $\chi^{spin}(T)$ can be described by the following ansatz,

[TABLE]

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 (25) leads to

[TABLE]

Our result for $\sigma^{2}(T)$ is plotted in figure 6. For $g\leq 0.5$, the data can be described by equation (34) 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.III.1.1. For $g\geq 0.8$, fitting $\sigma^{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 (34).

The extracted $M_{0}$ is plotted in figure 5(b). Close to the transition point at $g=g_{c}$, we expect $M_{0}$ will vanish as $M_{0}\propto(g-g_{c})^{\beta_{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\leq g\leq 1$ the value of $x$ obtained from equation (34) 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\geq 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}\in[0.74,0.88]$. Notice that the fitting with different $g_{min}$ all describe the $g\geq 1$ part of the data quite well. We thus take our final estimate of $g_{c}$ to be $g_{c}=0.8\pm 0.1$.

III.1.3 Local moment-Kondo transition

Now that we have established that the system resides in the local moment phase for $g>g_{c}\simeq 0.8$ at $\Gamma_{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 $\xi/\beta$, and a divergence in $\chi_{f}^{V}$, both around $\Gamma_{0}=0.4$ (cf. figure 7).

Similar to what we have done for $\sigma(\Gamma_{0},\beta)$ in equation (32), we consider the following finite size scaling hypothesis for $\xi$ and $\chi_{f}^{V}$,

[TABLE]

As seen in figure 8, close to the critical point the data fall nicely under a single universal curve. We obtain $\Gamma_{c}=0.35(2),\nu^{-1}=0.39(6)$ from $\xi$ and $\Gamma_{c}=0.34(2),\nu^{-1}=0.37(5)$ from $\chi_{f}^{V}$. Our final estimated value are $\Gamma_{c}=0.35(2)$ and $\nu^{-1}=0.38(5)$. The value of $\nu^{-1}$ obtained here for critical point C*′* is in sharp contrast with that for critical point C with $\nu=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 $\chi^{spin}$ vs. $T$ at different $\Gamma_{0}$. At the critical coupling $\Gamma_{0}=\Gamma_{c}\simeq 0.35$, $\chi^{spin}(T)$ can be fitted with a power law $\chi^{spin}(T)\propto T^{-x}$ with $x=0.65$. Inside the local moment phase at $\Gamma_{0}=0.1$, it can be described by equation (33) with a finite $M_{0}=0.10$ for the $M_{0}/T$ term and a sub-leading $1/T^{x}T_{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 $\chi^{spin}\sim 1/T^{s}$.

III.2 s=0.2

We now turn to the $s=0.2$ case. This is also the case investigated in reference Otsuki (2013) at the $U=\infty$ limit. We will fix $g=0.5$ and gradually increase $\Gamma_{0}$ to find the QCP from local moment phase to Kondo screened phase.

We first plot the dependence on $\Gamma_{0}$ of the Binder cumulant $U_{2}$ and the reduced correlation length $\xi/\beta$ in figure 10(a) and figure 10(b), where we have identified crossing points for both quantities. This suggests a transition from a local moment phase at small $\Gamma_{0}$ to a Kondo screened phase at large $\Gamma_{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 $\beta$ and $2\beta$ in figure 11. We see that the crossing points obtained from $U_{2}$ and $\xi/\beta$ are approaching to the same critical value $\Gamma_{c}$ in the $T=0$ limit from opposite directions. By extrapolating the crossing points $\Gamma_{cross}$ to $T=0$ using a simple power-law relation $\Gamma_{cross}=\Gamma_{c}+aT^{b}$, we find $\Gamma_{c}=0.48(1)$.

We can then repeat the analysis done in Sec.III.1.1 for the same type of transition at $s=0.2$ by considering scaling collapse of the form in equation (35) for the correlation length $\xi$ and similarly for the Binder cumulant $U_{2}$,

[TABLE]

where the presence of the sub-leading term $A/\beta^{\phi/\nu}$ 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 (37) and equation (35) is plotted in figure 12(a), and they give consistent estimates for the value of the critical coupling $\Gamma_{c}$ and correlation length exponent $\nu$. We obtain $\Gamma_{c}=0.49(1)$, $\nu^{-1}=0.42(3)$ from $U_{2}$ and $\Gamma_{c}=0.48(1)$, $\nu^{-1}=0.43(3)$ from $\xi$.

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 10(c) the measured $\chi^{V}_{f}$ appears to diverge near our estimated $\Gamma_{c}$. A finite size scaling analysis can be performed as well. For consistency we consider the same type of scaling form of $\chi^{V}_{f}$ as appeared in equation (36),

The result, plotted in figure 12(b), gives $\Gamma_{c}=0.46(2)$ and $\nu^{-1}=0.48(3)$, in fairly good agreement with what we have obtained from $U_{2}$ and $\xi$. Our final estimates are $\Gamma_{c}=0.48(1)$ and $\nu^{-1}=0.44(5)$.

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

We note that the calculation in reference Otsuki (2013) has assumed the relation $\chi^{spin}(T)\propto T^{-s}$ and use it as a tool to locate the QCP by looking for the crossing point of $T^{s}\chi^{spin}(T)$ at different $T$. But there the crossing point has significant drift versus temperature, which is consistent with an evolving $\alpha(T)$ in our calculation. Here we determine the critical coupling $\Gamma_{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 $\chi^{spin}(T)$ directly. Unfortunately from figure 13(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 $\chi^{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 $\chi^{spin}(T)$, even though the bare Kondo temperature $T_{K}^{0}$ is about $T_{K}^{0}=2.5$ at $\Gamma=\Gamma_{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 $\chi^{spin}(T)$.

III.3 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 1, 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 $\Gamma_{0}=0$, and vary both the bosonic coupling $g$ as well as the bosonic bath exponent $s$. As $s\rightarrow 0^{+}$, the procedure to obtain $J(\tau)$ defined in equation (19) and equation (21) will encounter convergence issue. As the critical property only depends on the long time asymptotic behavior of $J(\tau)$, we directly adopt a $J(\tau)$ that has the correct $1/\tau^{1+s}$ dependence as our input without specifying the actual form of $\rho_{b}(\omega)$. To be specific, we choose $J(\tau)$ to be the following,

[TABLE]

The exponential factor will make $J(\tau)$ finite at the end points: $\lim_{\tau\rightarrow 0}J(\tau)=\lim_{\tau\rightarrow\beta}J(\tau)=1$. Also $J(\tau)$ is even under reflection about $\tau=\beta/2$.

We can then integrate $J(\tau)$ twice to get $B(\tau)$,

[TABLE]

with $a=-\int_{0}^{\beta/2}J(\tau^{\prime\prime})d\tau^{\prime\prime}$ determined from the condition $dB(\tau)/d\tau|_{\tau=\beta/2}=0$.

Using equation (38) and equation (39) as input we have obtained the dynamical spin correlation function $\chi^{s}(\tau)$ for different value of $g$ and $s$. In figure 14(a) we present the result of $\chi^{s}(\tau)$ vs. $\tau$ at several different value of $\beta$ for the specific case of $g=0.4,s=0.2$. At each $\beta$, $\chi^{s}(\tau)$ drops from $1/4$ at $\tau=0$ and reaches its minimum at $\tau=\beta/2$. As $\beta$ is increased, $\chi^{s}(\tau=\beta/2)$ converge to a finite value.

We then plot the evolution of $\chi^{s}(\tau=\beta/2)$ obtained at low temperature, as a function of $s$ for four different choices of $g$ in figure 14(b), up to the smallest value of $\chi^{s}(\tau=\beta/2)$ that we can reach convergence. We can identify $\chi^{s}(\tau=\beta/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$, $\chi^{s}(\tau=\beta/2)$ decreases smoothly as a function of $s$. Furthermore, we can extrapolate each curve to larger value of $s$ until $\chi^{s}(\tau=\beta/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 15. Note that the shape of the phase boundary will depend on the specific form of $J(\tau)$ 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^{*}\simeq 0.47$. Although we do need to admit the simple extrapolation scheme performed in figure 15(a) could introduce some error.

IV 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\rightarrow 0^{+}$, or $\epsilon\rightarrow 1^{-}$, provided the relation $\eta=\epsilon$ holds, which initially is an statement made at critical point C from $\epsilon$-expansion perspective. Our calculation implies that C should disappear before $\epsilon$ reaches $1$, and that the actual impurity QCP encountered in the EDFMT calculation should be C*′* instead. Nonetheless, the relation $\eta=\epsilon$ 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 $\eta=\epsilon$ only relies on the condition $\eta=\epsilon+2\beta(g)/g|_{g=g^{*},J=J^{*}}$, which is shown to be valid to all orders in $\epsilon$ in referenceZhu and Si (2002). The relation $\eta=\epsilon$ then follows at any intermediate coupling fixed point $g=g^{*},J=J^{*}$, where $\beta(g)/g|_{g=g^{*},J=J^{*}}=0$, regardless of whether $g^{*}$ and $J=J^{*}$ is of the order $\epsilon$. Thereby this argument can be extended to C*′* as well.

V 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 $\epsilon$-expansion RG. At both types of QCP, our results suggest the spin correlation function obeys the power law $\chi^{spin}(\tau)\sim(1/\tau)^{\eta}$ with $\eta=1-s$.

VI 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.