S Florens

A Freyn

D Venturelli

R Narayanan

Dissipative spin dynamics near a quantum critical point: Numerical Renormalization Group and Majorana diagrammatics

We provide an extensive study of the sub-ohmic spin-boson model with power law density of states J(ω) ∝ ω s (with 0 < s ≤ 1), focusing on the equilibrium dynamics of the three possible spin components, from very weak dissipation to the quantum critical regime. Two complementary methods, the bosonic Numerical Renormalization Group (NRG) and Majorana diagrammatics, are used to explore the physical properties in a wide range of parameters. We show that the bosonic self-energy is the crucial ingredient for the description of critical fluctuations, but that many-body vertex corrections need to be incorporated as well in order to obtain quantitative agreement of the diagrammatics with the numerical simulations. Our results suggest that the out-of-equilibrium dynamics in dissipative models beyond the Bloch-Redfield regime should be reconsidered in the long-time limit. Regarding also the spin-boson Hamiltonian as a toy model of quantum criticality, some of the insights gained here may be relevant for field theories of electrons coupled to bosons in higher dimensions.

I. INTRODUCTION

Quantum dissipative models, introduced on a phenomenological basis in many areas of physics, have proved to be an invaluable tool for studying the quantum mechanical friction of an external bosonic environment on a generic two-state system. [START_REF] Leggett | [END_REF][START_REF] Weiss | Quantum Dissipative Systems[END_REF] While most studies on quantum dissipative models investigated questions pertaining to decoherence at short time scales, their longtime behavior, in particular near possible quantum critical points, is less well understood, see Refs. 3-5 for recent reviews. The motivation to investigate such regimes comes from the recent interest in zero temperature (quantum) phase transitions, [START_REF] Sachdev | Quantum Phase Transitions[END_REF] where anomalous low energy properties emerge due to the presence of quantum critical modes. Despite their apparent simplicity, dissipative impurity models are deeply concerned by this phenomenon, displaying often non-trivial long-time dynamics. They can be used to describe quantum criticality not only at the level of a single magnetic impurity diluted in a correlated system (such as Mott insulators [START_REF] Clarke | [END_REF][8][9] or magnetic metals 10,11 ), but also in bulk materials themselves (e.g. in quantum spin glasses, 12 or heavy fermion compounds 13,14 ) via the framework of the Dynamical Mean Field Theory. 15 Although we restrict our discussion here to a two-level impurity model coupled to a bosonic continuum, it is hard not to mention the host of fascinating physical phenomena in fermionic models, with the Kondo problem [START_REF] Hewson | The Kondo Problem to Heavy Fermions[END_REF] and its variations being another major example for dissipation (via magnetic screening) and quantum critical phenomena (see Ref. 3 for a recent theoretical review and Ref. 17 for a discussion of some experimental realizations in quantum dots).

Despite their quantum nature, the phase transitions in dissipative quantum impurity models are generally believed to be well described by an associated long-range classical one-dimensional Ising model. 18 Indeed the subohmic spin-boson model, in which a single spin S = 1/2 is subjected to a magnetic field ∆ perpendicular to the quantization axis, and coupled (along the quantization axis) to a bosonic continuum with vanishing density of states J(ω) = 2παω 1-s c ω s θ(ω)θ(ω cω), obeys such quantum-to-classical mapping in the range 1/2 < s < 1, as confirmed by recent direct simulations of the quantum model. 19 In this case, the quantum critical behavior follows non-trivial classical exponents of an Ising chain with long-range couplings, decaying as 1/(ττ ) 1+s in the imaginary time domain (τ variables), and is described by an interacting fixed point. 20,21 However, the behavior in the strongly sub-ohmic case, 0 < s < 1/2, has been recently debated. According to the quantum-toclassical mapping, the model should fall above its upper critical dimension and show mean-field exponents with a violation of hyperscaling laws. Initial evidence for a different behavior, namely non-mean-field exponents above the upper critical dimension for the strongly sub-ohmic spin-boson model, came with surprise from Numerical Renormalization Group (NRG) calculations. 22 Ref. 23 also provided similar puzzling results from Monte Carlo simulations based on a truncation of the long-range in-teraction (which may or may not be in the same universality class as the spin-boson problem). At this juncture it is to be noted that these results have appeared recently to be controversial, [24][25][26][27][28][29] and at present no final consensus has appeared on the definite nature of this continuous quantum phase transition (QPT). Still, there are hints that the NRG displays truncation errors in the long-range ordered phase 25,26,29 , so that the quantumto-classical mapping may hold after all. We note that a discontinuous first order transition, as obtained from variational methods, 30 is likely an artifact resulting from a breakdown of this approximation in the quantum critical regime, and that the localization transition in the spin-boson model is currently believed to be second order.

In the present paper we will not attempt to answer all these unresolved issues, yet we hope to make progress on the quantitative understanding of the spin-boson model near its quantum critical point, focusing on the so-called delocalized regime, which corresponds to the classically disordered phase. Two complementary tools will be developed, namely NRG calculations 22,[31][32][33] and Majoranafermion diagrammatic theory, 34,35 that will be benchmarked against each other, both far from and close to the quantum critical point. Although the Majorana modes appear here as a technical device to reproduce the spin dynamics of a standard two-level system, the methodology developed here may be useful too in the emerging field of Majorana qubits. 36 Our work will also aim at an exhaustive study of the physical properties of the model, and we will examine in detail the zero-temperature equilibrium spin dynamics in all three possible spin directions, which was previously not achieved to our knowledge. In this context, it is also interesting to use the NRG as a testbed for diagrammatic methods in fermionboson models, and we will discuss how various GW-like schemes compare with the numerics. Interestingly, we will show that the phase transition is driven by a Bosecondensation of the mode mediating the interaction between Majorana fermions, and that ladder resummation within the bosonic self-energy is the crucial ingredient to obtain the precise location of the quantum phase transition. This allows us to establish (at two-loop order) a phase boundary given by the following critical dissipation α c = (s/2 + s 2 /4)|∆/ω c | 1-s + O(s 3 ), that describes the NRG phase diagram quite accurately for all values of ∆/ω c < 1 as long as s < ∼ 0.8. We believe that the precise understanding achieved here will be important to make further progress in the elucidation of quantum critical properties of the model. It is also worth stressing that the non-equilibrium dynamics, 1 more often considered in the context of spin-qubits, [START_REF] Shnirman | Quantum Noise in Mesoscopic Physics[END_REF] should also be deeply affected at long-times in case of proximity to a quantum critical points, see e.g. Ref. 38 for a related study. Standard master equation methods at the level of Bloch-Redfield approximations (analogous to lowest order perturbation theory in our equilibrium computations) should fail in this regime, and ought to be reconsidered seriously, pos-sibly at the light of further developments of both perturbative techniques 39,40 and out-of-equilibrium NRG. 41 The plan of the paper is as follows. Section II will present the general properties of the spin-boson model, and its solution at zero-temperature by the NRG. An improved broadening method 42 will be used to extract the dynamical spin susceptibilities (the longitudinal as well as the transverse ones). Section III will present the Majorana-fermion diagrammatic method, based on perturbation theory in the dissipation strength, focusing on the regime of very weak dissipation far from the quantum critical point. Then, we will extend the diagrammatics up to the vicinity of the quantum phase transition in Section IV. Leading logarithmic corrections in the form of ladder diagrams will be obtained in a two-loop calculation of the bosonic self-energy, leading to the accurate analytical formula for the phase boundary quoted above. Such renormalization effects, taking into account quantitative vertex corrections, will allow to extend the Majorana diagrammatics non-perturbatively, and match the numerical data both at low (critical modes) and high energy (dissipative features). Further directions of research will be proposed as a conclusion, and several appendices will contain some technical details, such as a derivation of Shiba's relation for the sub-ohmic spin-boson model, and an heuristic discussion of the quantum-to-classical mapping via the effective bosonic action.

II. SPIN-BOSON MODEL AND ITS SOLUTION BY THE NRG

A. Hamiltonian and physical aspects

The (sub)-ohmic spin-boson Hamiltonian involves a three-component half-integer spin S = h 2 σ (we set h = 1 in what follows), with σ the three Pauli matrices, and a continuous bath of bosonic oscillators a † i :

H = ∆ 2 σ x + 2 σ z + λ 2 σ z i (a † i + a i ) + i ω i a † i a i . (1)
∆ and are magnetic fields applied to the spin (respectively orthogonal and parallel to the quantization axis), and the coupling constant λ controls the strength of the dissipation to the environment. The bosonic spectrum (which includes by standard definition the dissipative coupling λ) is assumed to be a power-law up to a sharp high-energy cutoff ω c :

J(ω) ≡ i πλ 2 δ(ω -ω i ) = πλ 2 ω s ω 1+s c θ(ω)θ(ω c -ω) J(ω) = 2παω 1-s c ω s θ(ω)θ(ω c -ω) (2)
thus defining the dimensionless dissipation strength α = λ 2 /(2ω 2 c ), an important parameter of the model. In all what follows, we will consider a vanishing magnetic field parallel to the bosonic bath = 0, focusing only on the interplay of dissipation and spin precession, which are driven respectively by the parameters α and ∆. Temperature will also be set to zero, so that the model (1) displays two distinct kinds of ground state: i) a delocalized phase, with zero magnetization along the z axis, S z = 0, in which the ground state is non-degenerate and is adiabatically connected to a coherent superposition of the two spin up/down states as the dissipation strength α is increased from zero (keeping a finite value of the magnetic field ∆); ii) a localized phase, with a doubly degenerate ground state, where dissipation is strong enough to polarize the spin, leading to a non-zero magnetization S z = 0, associated to a spontaneous symmetry breaking. We emphasize that in the present quantum problem, long-range correlations associated to the slow decay in time of the bath correlations can be sufficient to generate long-range order at zero-temperature. When 0 < s ≤ 1, a quantum phase transition therefore occurs between the two phases at a critical value α c of the dissipation strength. Having in mind a perturbative expansion in the dissipation, we will focus in this work on the delocalized phase only up to the quantum critical point, i.e. α < α c . We also note that in both phases, a finite magnetization along the applied magnetic field ∆ is always present, S x = 0.

Several exact results are known for the model (1), see the reviews in Refs. 1,3. For s = 1 (ohmic case), the quantum phase transition is of Kosterlitz-Thouless type, and is intimately connected to the ferromagneticantiferromagnetic transition in the Kondo model. This connection, that can be precisely shown with the use of bosonization, [START_REF] Leggett | [END_REF]43 makes it possible to demonstrate that the critical dissipation strength is exactly α c = 1 in the limit ∆ ω c . A particularity of the ohmic case (s = 1) is that the crossover from coherent to incoherent spin dynamics (corresponding respectively to weakly damped and overdamped Rabi oscillations in the transverse spin autocorrelation functions) occurs before the onset of the QPT, so that the localized phase has always an incoherent character, see Fig. 1. In contrast, the sub-ohmic model (0 < s < 1) can undergo localization before or after the disappearance of the coherent oscillations, depending on the precise values of the microscopic parameters, [START_REF] Nalbach | [END_REF]41 see Fig. 1. In particular, in the limit of small magnetic field ∆/ω c 1, the quantum phase transition occurs at small dissipation α 1, so that low-energy dynamics (related to quantum criticality) and high energy physics (the damping of the Rabi oscillations) are governed by different energy scales, and both can be treated in principle by perturbative methods. We will focus in this section on a Numerical renormalization group (NRG) solution that encompasses both weak (α 1) and nonperturbative dissipation (α 0.1) regimes, before making comparison with analytical methods in Sections III and IV for the weak dissipation case. cases respectively, as a function of the dissipation parameter α and the dimensionless magnetic field ∆/ωc. The superohmic (s > 1) regime shows only the crossover from coherent to incoherent dynamics, without a localization transition. In all cases, the coherent/incoherent crossover happens at comparable coupling strength, namely α 0.5.

B. Bosonic NRG algorithm

The implementation of the NRG method follows the standard procedure initially introduced for fermionic models, 31,32 and recently extended to bosonic Hamiltonians. 33 First, the spin-boson model ( 1) is rewritten in a continuous form:

H = ∆ 2 σ x + 1 2 σ z ωc 0 d h( ) [a † + a ] + ωc 0 d a † a (3) with J( )/π = [h( )] 2 .
The bosonic bath J( ) is then logarithmically discretized using the Wilson parameter Λ > 1, first on the highest energy interval near the cutoff [Λ -z ω c , ω c ], and then iteratively on successive decreasing energy windows [ω n+1 , ω n ] with ω n = Λ -n-z+1 ω c (for n strictly positive integer). This choice of discretization introduces the so-called z-parameter (with 0 < z ≤ 1) that is used to average over N z different realizations of the Wilson chain 44 (this amounts to gradual changes of boundary conditions, allowing to obtain better statistics from the numerical simulations). Taking N z different z values, usually uniformly distributed in the ]0, 1] interval, allows to obtain an order N z improvement of the spectral resolution at finite energy. For very narrow spectral structures in the correlation functions, this method becomes too expensive and and one lays recourse to improved broadening techniques, see Ref. 42 and Sec. II C.

The bosonic fields are then decomposed in Fourier modes (with p a positive or negative integer) on each interval ω n+1 < < ω n of width d n = (1 -Λ -1 )Λ -n-z for n > 0 (note that d 0 = 1 -Λ -1 ):

a † = n,p e i2πp /dn √ d n a † n,p . (4)
This step is of course exact, and the first NRG approximation consists in neglecting all p = 0 modes, keeping only the operator a † n ≡ a † n,0 (the full Hamiltonian is thus only recovered in the Λ → 1 limit. 19 ) This leads to the "star"-Hamiltonian:

H star = ∆ 2 σ x + 1 2 σ z +∞ n=0 γ n √ π [a † n + a n ] + +∞ n=0 ξ n a † n a n (5)
with the "impurity" coupling strength (the expressions below apply for n > 0 only, the n = 0 value can be obtained analogously)

γ 2 n = ωn ωn+1 dω J(ω) = 2πα 1 -Λ -(s+1) s + 1 ω 2 c Λ -(n+z)(s+1) (6
) and the typical energy ξ n in each Wilson shell

ξ n = 1 γ 2 n ωn ωn+1 dω ω J(ω) = s + 1 s + 2 1 -Λ -(s+2) 1 -Λ -(s+1) ω c Λ -n-z . (7)
One then performs an exact mapping onto the so-called chain Hamiltonian:

H chain = ∆ 2 σ x + 1 2 σ z η 0 π [b † 0 + b 0 ] (8)
+ +∞ n=0 [ n b † n b n + t n (b † n+1 b n + b † n b n+1 )]
where η 0 = d J( ) and the on-site energies n and intersite hoppings t n are determined numerically by a tridiagonalization of the bosonic part of the star Hamiltonian (5). The remaining step of the NRG follows Ref. 33, adding successively sites into H chain starting with n = 0. The site of order n roughly describes the physics for energies of the order ω c Λ -n . This iterative procedure allows to reach exponentially small scales in a linear effort, an important achievement for the study of Kondo physics 31 or impurity quantum critical behavior. [START_REF] Vojta | [END_REF] It is important to stress that the success of the NRG relies on the exponential decrease of the chain parameters n and t n (related to the analogous behavior of the star parameters ξ n and γ n ), which makes the iterative diagonalization reliable and stable (see however Ref. 45 for the extension of NRG to other classes of models). Because the size of the Hilbert space needs to remain finite, truncation constraints need to be implemented. Typical calculations are usually performed with N b = 10 bosonic states on each successive bosonic site of the Wilson chain that is added during the renormalization procedure, 33,42 and N = 160 kept NRG states (truncation). Matrix sizes thus do not exceed 2000 × 2000 during the whole computation, and we proceed to N = 30 NRG iterations with Λ = 2 , so that a typical low energy scale ω min = ω c Λ -N 10 -9 can be reached. The resulting discrete energy spectra at successive NRG iterations are combined using the interpolation scheme proposed in Ref. 46 (see however the more rigorous implementation given by full density-matrix NRG calculations 47,48 ), leading to a set of z-dependent many-body energy levels k,z labelled by quantum number k. The spin-spin correlation function along the i-axis (i = x, y, z) is defined from the equilibrium spin autocorrelation functions in real time and their Fourier transform to real frequency:

C i (t) = 1 2 [σ i (t), σ i (0)] (9)
C i (ω) = dω 2π C i (t) cos(ωt). ( 10
)
An alternative way to introduce the spin dynamics is from the imaginary time spin correlation functions and their Fourier transform onto Matsubara frequency:

χ(τ ) = S i (τ )S i (0) (11)
χ(iν) = β 0 dτ e iντ χ(τ ) (12)
χ(ω + i0 + ) = χ (ω) + iχ (ω) (13)
where analytic continuation was performed in the last equation in order to obtain the retarded spin susceptibility. These susceptibilities are simply related to the spin autocorrelation functions (10) by the fluctuationdissipation theorem, which reads at equilibrium and for zero temperature (we consider the T = 0 limit in all what follows):

C i (ω) = Sign(ω)χ i (ω). ( 14
)
We will therefore use equivalently both nomenclatures in the rest of the paper. There exists an important equation, known as Shiba's relation, 33,43 which connects real and imaginary parts of the spin susceptibility at low frequency. This exact relation is usually derived for the ohmic (s = 1) spin-boson model (or equivalently for the anisotropic Kondo model) and reads at small frequency:

C z (ω) = 2πα|ω|[χ z (0)] 2 . ( 15
)
We will prove in appendix A that this result can be generalized to the sub-ohmic case (0 < s < 1) in the following form

C z (ω) = J(|ω|)[χ z (0)] 2 (16)
valid in the small ω limit. An useful byproduct of Shiba's relation is that, at a magnetic ordering of the spin (corresponding to the quantum critical point), the susceptibility χ z (0) diverges, so that the low-energy power law behavior C z (ω) ∝ J(|ω|) ∝ |ω| s obeyed in the whole delocalized phase should turn into a different (and diverging) power law.

A last important technical step is that all these correlation functions can be computed at zero temperature from the raw NRG data using Lehmann's decomposition rule:

C i (ω) = 1 2N z k,z | 0, z|σ i |k, z | 2 δ (|ω| + 0,z -k,z ) (17)
where 0,z is the ground state energy. This leads to a superposition of sharp δ-peaks which need to be broadened. This sensitive issue is examined now.

C. Optimized broadening method

A single NRG calculation (for a given value of z) usually provides a set of energy levels which come by packets located around each Wilson shell. The energy resolution at scale ω n = Λ -n-z is thus usually thought to be of the order

d n = ω n -ω n+1 = (1 -Λ -1
)ω n , and degrades at higher energy. In order to generate smooth NRG spectra, the delta-peaks in the Lehmann formula (17) are therefore usually broadened 32 at energy ω n on a scale of the same order (the broadening scale is typically bω n , with b 0.7) using the substitution:

δ (|ω| -ω n ) → e -b 2 /4 ω n b √ π e -[log(|ω|/ωn) 2 /b] 2 . ( 18
)
As a result, spectral features at a given energy ω that are sharper than ω itself cannot be well resolved, and will come out overbroadened. Some relative degree of improvement can be achieved by combining N z NRG runs with different values of the z parameter (z-averaging procedure, or z-trick 44 ). Although the mathematical justification for this procedure is still unclear, in the optimal case the resulting energy packets turn out to be uniformly distributed if the set of z-values is also chosen uniformly, and this allows a decrease of the broadening parameter to b 0.7/N z . In the case of very sharp features, this method requires a large number of NRG runs, and will be prohibitive. It may also become problematic in regions where the z-averaging does not provide a uniform distribution of states (this may occur from the accidental disappearance of some NRG states), in which case uncontrolled oscillations of period Λ can be generated (situation of underbroadening). These problems have led authors to speculate that quantitative NRG spectral functions can be extracted only in the continuum limit, [49][50][51] either with Λ → 1 or N z → +∞ (only the former is mathematically sound, but the NRG algorithm cannot be managed anymore as scale separation breaks down). Both limitations of the z-averaging can however be lifted thanks to a simple remark made by two of us in a previous publication: 42 for a single NRG run, the levels within a given logarithmic energy window are not uniformly distributed, but rather tend to bunch together close to sharp resonances. This effect becomes quite obvious when the width of such resonance becomes extremely small, see Fig. 2, but is always present in all NRG calculations. Figure 2 the spin susceptibilities of the spin-boson model, computed with N z = 20 interleaved z-averaged NRG runs (dots) and smoothed using two standard NRG broadenings (dash-dotted curves) with respective parameters b = 0.7 and b = 0.7/N z = 0.035. The largest broadening clearly fails to produce the expected peak (as a testbed, an accurate analytical expression is provided by the full line, see Sec. III for further details), while the smallest broadening signals the peak, but overestimates its width by more than one order of magnitude! Large scale z-averaging with N z 10 4 would probably allow to resolve the correct spectral structure, but is unreasonably expensive. The key observation 42 is that the NRG eigenvalues do not always come in packets uniformly distributed within the Wilson shell, but on the contrary tend to cluster close to resonances, as clearly seen on the raw data (dots) in Fig. 2. This provides the missing information which is needed in order to achieve a good spectral resolution with a limited numerical cost. The main strategy to be used is that the broadening parameter b must become energy-dependent, which we previously called the "b-trick". 42 However, the generic implementation of this simple idea is not fully clarified: we present here a simple scheme that seems to be quite efficient, but we emphasize that a more generic and robust method is still to be found. The required algorithm must indeed adapt the broadening parameter b(ω) in Eq. ( 18) to the frequency dependence of the density of highest weight NRG peaks. Our approach, which draws inspiration from the NRG logarithmic discretization, is to extract b(ω) from the logarithmic derivative of the integrated spectrum up to frequency ω:

b(ω) = b 0 2   q + d log ω 0 C d log ω -1 + q + d log ω +∞ C d log ω -1   (19)
where q 1 is a regularization parameter whose precise value does not matter much to the final result, and b 0 provides the typical broadening at low and high frequencies (far from the atomic resonances). Note that we have used here two different frequency sweeps, one from ω = 0 and one from ω = +∞, in order to treat on an equal basis low and high frequency tails (this procedure is not satisfactory in cases where several resonances are present). Because the actual NRG data is fully discrete, see Eq. ( 17), we compute b(ω) recursively using Eq. ( 19) on the broadened NRG spectra. This procedure converges after few iterations to the results displayed on Fig. 2. It is clearly remarkable that such a good agreement can be reached with so little numerical effort, demonstrating that the raw NRG data encode much more information than previously thought (see Ref. 42 for further details). One important advantage of the b-trick is that b(ω) takes very small values near resonances, greatly enhancing the resolution, while keeping large values away from the peak, thus avoiding the usual underbroadened NRG oscillations. The main tuning parameter in our method is the low-frequency broadening b 0 in Eq. ( 19), which is determined by a simple convergence method. We finally emphasize that some degree of z-averaging is always needed in order to generate smooth spectra with the b-trick, but in our experience z-averaging alone is never sufficient to obtain accurate results near sharp finite-frequency resonances.

D. NRG results for the dynamical susceptibilities in the three spin directions

We present here the NRG data for a fixed value of the magnetic field ∆/ω c = 0.1, varying the dissipation from values α α c (far from the quantum critical point) up to α c . We choose three different values of the bath exponent s = 0.1, 0.5, 1.0, which span the complete generality of the model. We remark that the case s = 0 is special, as the spin is localized for all α > 0 at zero temperature 41 , and will not be considered here, as we aim at understanding the delocalized phase of the model only. We will also not examine the situation of s > 1, which shows for weak dissipation similar features as the s = 1 case, and does not present a phase transition at larger coupling. We stress already that depending whether α c 1 or not, the transition will or won't take place at weak coupling, and this can lead to different physical pictures.

Let us consider first the case s = 0.1 and the spin susceptibilities along each spin direction (i = x, y, z) shown in Fig. 3 (the plot also includes a comparison to lowest order pertubation theory, introduced in Section III B). The 0.01 0.1 10 100 3: Spin susceptibilities Cz(ω), Cy(ω) and Cx(ω) (top to bottom panels) for ∆/ωc = 0.1, s = 0.1 and increasing dissipation strength α = 0.001, 0.004, 0.007 (bottom to top). The critical point is located at αc 0.007. The NRG plots were generated with Nz = 20 interleaved z-averages and used the b-trick discussed in Sec. II C, while the perturbative expressions (G 0 W 0 scheme) are given in Eqs. (36-37). Note that the upper panel has a different range in ω than the two other panels.

C z (ω) 10 -8 10 -6 10 -2 10 0 ω NRG G 0 W 0 scheme 10 -6 10 -4 10 0 10 2 C y (ω) 0.001 0.01 0.1 1 ω NRG G 0 W 0 scheme 10 -4 10 -3 10 -1 C x (ω) 0.001 0.01 0.1 1 ω NRG G 0 W 0 scheme FIG.
transverse (orthogonal to the field and parallel to the bath) spin susceptibility C z (ω), already considered, 41,42 shows two interesting features. First, the resonance at ω = ∆ broadened by the bath is always sharply defined, because here α c 0.007 is so small that one remains in the perturbative regime for all α < α c values regarding this high energy features. This is consistent with the near perfect agreement with the perturbative result in this range of frequencies, and weakly damped oscillation should occur in the real-time spin dynamics. The low-energy part of this spectral function is yet more interesting. For α α c , a slow decay with C z (ω) ∝ |ω| s is observed, also consistent with the perturbative result. Increasing α, a crossover at a scale ω is seen between the |ω| s behavior at ω ω and a different power law |ω| -s at ω ω ∆. At the quantum critical point α = α c , ω vanishes and a complete |ω| -s divergence is realized throughout. Clearly, the perturbative expression fails to reproduce this behavior. We stress here that even values relatively far from the critical point, such as α = 0.004, show large quantitative deviation from lowest order perturbation theory at low frequency, despite the good agreement on the coherent peak. This shows that the range of validity of perturbative methods (such as Bloch-Redfield), although possibly accurate at short time scales, is very limited in the long time limit. Clearly the high energy dissipation mechanism does not care for the complex low-energy behavior of the spin dynamics, even at the quantum critical point. 26,33 Turning to the second transverse susceptibility C y (ω) (orthogonal to both the magnetic field and the spin-bath term), we observe similar behavior as in C z (ω) for the resonant peak, and good agreement with perturbation theory. The low energy part of the spectrum is however much less dramatic, because critical modes, related to fluctuations driven by the bath, pertain mainly to the z-component of the spin. What is actually going on is a very mild crossover from ω 2+s at ω ω to ω 2-s at ω ω ∆ (the latter behavior includes the whole range of low frequencies at the critical point), due to the small value of s = 0.1. Agreement with perturbation theory appears thus much better, although small deviations, due to the tiny changes of the power law exponent, can be seen at low frequency. These precise (and exact) values for the various power laws will be demonstrated in Section III. Finally, the longitudinal (parallel to the bath) spin susceptibility C x (ω) presents quite different features. Clearly in the absence of dissipation (α = 0), a δ(ω) peak at zero frequency occurs, in contrast to the δ(ω -∆) peak at the magnetic field frequency displayed by both C y (ω) and C z (ω). Thus no sharp resonance is generated at small dissipation, and rather a broad shoulder emerges, quite reminiscent of the T-matrix of the Kondo model in a magnetic field. 52 This structure is well accounted for by the perturbative calculation for ∆ < ω, however, stark discrepancy is seen at low frequency for |ω| < ∆: the NRG data present a tail of low energy modes, while perturbation theory provides a gap in this range. This new feature given by the NRG will be elucidated in Section III C, in connection with multiparticle effects in the diagrammatics, that are neglected in the current lowest order calculation.

We now turn to the intermediate value of s = 0.5, giving rise to a somewhat larger critical dissipation strength α c 0.105, see Fig. 4. In that case, the smallest chosen value of α = 0.02 still lies in the perturbative regime, so that again near perfect agreement with perturbation theory can be checked. However, the approach to the quantum critical point shows larger quantitative deviations, both at high and low energy. Clearly, the coherent peak is more strongly smeared in the NRG for the The NRG plots were generated with Nz = 10 interleaved z-averages and used the btrick discussed in Sec. II C, while the perturbative expressions (G 0 W 0 scheme) are given in Eqs. (36-37). Note that the upper panel has a different range in ω than the two other panels. larger shown values α = 0.06, 0.1 compared to perturbation theory, due to the now relatively important magnitude of α. Nevertheless, the peak structure survives up to the critical point, so that moderately damped oscillations should show up in the time domain, even close to α c . Regarding the low energy part of the spectrum, the same crossover from ω s to ω -s behavior (for ω ω and ω ω ∆ respectively) is seen in C z (ω), now more clearly defined due to the larger value of s = 0.5. Again perturbation theory fails in this regime, and this is now more easily seen on C y (ω) as well, which presents respectively a crossover from ω 2+s to ω 2-s in the frequency dependency. We note also a similar enhancement of the low energy tails for the last component C x (ω) below the shoulder at ω < ∆.

We finally consider the so-called ohmic case s = 1, for which α c 1, clearly out of the perturbative range, see Fig. 10 -4 10 -2 10 0 ω NRG G 0 W 0 scheme FIG. 5: Spin susceptibilities Cz(ω), Cy(ω) and Cx(ω) (top to bottom panels) for ∆/ωc = 0.1, s = 1 and increasing dissipation strength α = 0.1, 0.5, 0.9 (bottom to top). The critical point is located at αc 1. The NRG plots were generated with Nz = 3 interleaved z-averages and used the b-trick discussed in Sec. II C, while the perturbative expressions (G 0 W 0 scheme) are given in Eqs. (36-37).

defined resonance and reasonable agreement with perturbation theory, although some deviations are visible. However, all values of α > ∼ 0.5 show a complete disappearance of the coherent peak at ω = ∆, leaving only the crossover scale ω , which separates regimes where C z (ω) ∝ ω for ω ω to C z (ω) ∝ ω -1 (up to logarithmic corrections) for ω ω ω c . In that case, oscillations are completely overdamped, and coherence is lost already well before the quantum critical point α c = 1 where the spin becomes logarithmically free (this point corresponds precisely to the ferromagnetic scaling limit in the associated Kondo model). Again, C y (ω) shows a similar behavior as observed for smaller s values, with a crossover from ω 2+s = ω 3 to ω 2-s = ω power laws for ω ω and ω ω ∆ respectively. Finally C z (ω) presents also a similar crossover in the low-energy tails, and a broad shoulder at the magnetic field value (for weak dissipation).

We now turn to the analytical description of this physics, and show that both qualitative and quantitative aspects of the numerical data can be well understood on the basis of the Majorana-fermion diagrammatic method.

III. MAJORANA DIAGRAMMATICS

A. Spin representation and correlation functions

The idea we henceforth develop is to use a Majorana description 34,35 for the impurity quantum spin S = 1 2 σ, in order to capture quantitatively the spin fluctuations in the whole delocalized phase up to the quantum critical point using a weak dissipation expansion (understanding the quantum phase transition coming in from the localized phase is a more delicate issue, that we will not attempt here). Technically, our methodology involves the introduction of a triplet of real fermions η j (with j = x, y, z) that satisfy the anticommutation relations {η j , η k } = δ jk , leading to the faithful representation S j = -(i/2) kl jkl η k η l . We emphasize that Majorana qubits 36 involve a doublet of real fermions, while the mathematical decomposition of a two-level qubits requires a triplet of real fermions. The main motivation behind considering the Majorana representation instead of more standard ones (Abrikosov or Popov fermions, Holstein-Primakoff or Schwinger bosons) is that one can alternatively write S = χ η, where χ = -2iη x η y η z is a fermionic operator which commutes with the Hamiltonian (1). This implies the crucial property that the spin susceptibilities Eq. ( 10) are obtained as simple identities involving standard one-particle Majorana Green's function (here at zero temperature):

C i (ω) = - 1 π Im G i (ω), (20)
where G i (ω) is the retarded real-frequency Majoranafermion Green's function, obtained from the Fourier transformation of the time-ordered Green's function G i (t) = -iθ(t) {η i (t), η i (0)} . This key relation between C i and G i is in great contrast to all other spin representations where the spin correlation functions are given in terms of four-fermion correlators, making the diagrammatics much more cumbersome. Therefore, thanks to Eq. ( 20), both the dissipative features at high energy and the low-energy quantum dynamics close to the quantum critical point can be simply encoded via Dyson equation by a Majorana-fermion self-energy. To obtain these important relations, we need to reexpress the spin-boson model Eq. ( 1) in the Majorana language:

H = -i∆η y η z -iλη x η y j (a † j + a j ) + j ω i a † j a j . (21)
The "diagonal" part given by the magnetic field is readily diagonalized into "free" (or bare) Majorana Green's functions:

G 0 x (ω) = 1 ω + i0 + (22) G 0 y (ω) = ω (ω + i0 + ) 2 -∆ 2 (23)
G 0 z (ω) = ω (ω + i0 + ) 2 -∆ 2 (24)
which are readily interpreted by a spin precession in the yz plane around the magnetic field pointing in the x direction. Interaction between the Majorana fermions is generated by the coupling to the bosonic bath, which can be expressed as a bare vertex involving the η x and η y field, see Fig. 6. This coupling results in three Majorana self-energies Σ i (ω), but clearly the z-component vanishes altogether, Σ z (ω) = 0, due to the form of the interaction.

Inverting the Majorana propagator (written in matrix form):

Ĝ =   ω -Σ x (ω) 0 0 0 ω -Σ y (ω) i∆ 0 -i∆ ω   -1 (25)
provides Dyson's equation for the full retarded Majorana Green's functions that encode the complete spin dynamics:

G x (ω) = 1 ω -Σ x (ω) (26)
G y (ω) = ω ω[ω -Σ y (ω)] -∆ 2 (27) G z (ω) = ω -Σ y (ω) ω[ω -Σ y (ω)] -∆ 2 (28)
The bare Green's functions (22-24) are correctly reproduced in the zero-dissipation limit, replacing Σ i (ω) by a vanishing imaginary part -i0 + . Because the interaction term in the Hamiltonian Eq. ( 21) involves a particular combination of the bosonic modes, it is useful to define the full bosonic propagator G Φ (t) = -iθ(t) [Φ(t), Φ(0)] of the "local" field Φ ≡ -i λ √ α j (a † j + a j ). This allows us to introduce a bosonic self-energy Σ Φ (ω) for the full Green's function

G Φ (ω) = G 0 Φ (ω) 1 -G 0 Φ (ω)Σ Φ (ω) (29)
in terms of the bare bosonic propagator

G 0 Φ = 1 α j λ 2 ω + i0 + -ω j - λ 2 ω + i0 + + ω j . (30)
Using the convention chosen to normalize Φ, the bare vertex takes simply the value √ α (see Fig. 6), and the bare bosonic Green's function can be expressed through the bosonic spectral density J(ω) given by Eq. ( 2):

G 0 Φ (ω) = - 1 α d π J(| |) Sign( ) ω + i0 + - = -2ω 1-s c ωc -ωc d | | s Sign( ) ω + i0 + - (31)
a quantity which is indeed independent of the dissipation strength α, now encoded (by convention) only in the bare vertex. We thus find an explicit expression for the imaginary part of the bosonic Green's function:

ImG 0 Φ (ω) = 2πω c Sign(ω) ω ω c s θ(ω 2 c -ω 2 ). (32)
The real part of the integral (31) can be computed analytically from Kramers-Kronig's relation only in the lowenergy limit |ω| ω c :

ReG 0 Φ (ω) = 4ω c s -2πω c cotan πs 2 ω ω c s ( 33
)
although the complete formula (31) will be used for later numerical computations. The low-energy expression (33) breaks down for s = 1, where additional logarithmic corrections arise (hinting at the known connection to the Kondo model 1 ). The exact resummation of perturbation theory requires the accurate knowledge of the full fermionic and bosonic self-energies, as well as the full three-particle vertex Γ (a quantity that depends on three independent external frequencies). The general diagrammatics is given formally in Fig. 7. We note that replacement of the bare propagators by the full ones without considering together the vertex corrections must be done in general with care, and can lead to serious pathologies, see Ref. 53. An exception to this rule is given in Sec. III C, and the precise role of the interaction vertex will be elucidated in Sec. IV E.

B. Lowest order perturbative self-energies in the Majorana diagrammatics

For very weak dissipation α 1 and far from the quantum critical point α α c , the spin dynamics should in 7: Exact formal diagrammatic expression for the fermionic and bosonic self-energies in terms of the full (exact) fermionic and bosonic Green's functions (Gi and GΦ respectively) and the full (exact) three-particle vertex Γ.

Σ exact y = G x G Φ Γ Σ exact x = G y G Φ Γ Σ exact Φ = G x G y Γ FIG.
Σ G 0 W 0 y = G 0 x G 0 Φ Σ G 0 W 0 x = G 0 y G 0 Φ FIG. 8:
Lowest order perturbative Majorana self-energies, obtained from bare fermion and bare bosonic propagators (G 0 W 0 scheme).

principle be well described by the lowest order perturbative Majorana self-energies (we will see in Sec. III C that this naive statement is not completely correct in some specific frequency range). These self-energies are obtained at order α by a single boson exchange, see Fig. 8, where both fermionic and bosonic lines are given by the bare propagators Eqs. (22-24) and Eq. ( 31), respectively. In terms of resummation of diagrams, this level of approximation corresponds to the so-called G 0 W 0 scheme for the generic fermion-boson models arising in the fields of heavy fermions [START_REF] Hertz | [END_REF] or strongly interacting electron liquids [START_REF] Giuliani | Quantum Theory of the Electron Liquid Cambridge[END_REF] (G 0 denotes the bare fermionic propagator and W 0 ≡ G 0 Φ the bare bosonic one). It is alternatively called the (non-self-consistent) Born approximation in the context of disordered systems, [START_REF] Castro | The Electron Liquid Paradigm in Condensed Matter Physics[END_REF] and it roughly corresponds to the weak-coupling Bloch-Redfield approximation for the real-time non-equilibrium dynamics of the spin-boson model. [START_REF] Shnirman | Quantum Noise in Mesoscopic Physics[END_REF] Note that the present Majorana diagrammatics is formulated for the equilibrium dynamics only. Extension to temporal dynamics under a sudden switching (of the bath or transverse field) could be incorporated in our formalism by using two-times Keldysh Majorana Green's functions, which lies however beyond the scope of the present work. Standard evaluation of the diagrams pro-vides the following imaginary part of the Majorana selfenergies at zero temperature:

ImΣ G 0 W 0 x (ω) = α ω 0 d π ImG 0 y ( )ImG 0 Φ (ω -)(34)
ImΣ G 0 W 0 y (ω) = α ω 0 d π ImG 0 x ( )ImG 0 Φ (ω -)(35)
Using the bare propagators (22-24) and ( 30), one readily obtains:

ImΣ G 0 W 0 x (ω) = -παω c θ[ω 2 c -(|ω| -∆) 2 ]θ(ω 2 -∆ 2 ) × ||ω| -∆| s ω s c ( 36
)
ImΣ G 0 W 0 y (ω) = -παω c θ(ω 2 c -ω 2 ) |ω| s ω s c ( 37
)
We are now equipped to start understanding the NRG results discussed in Section II D. Owing to the general relation Eq. ( 20) between spin autocorrelation functions and Majorana propagators, we restrict the discussion to the latter in the rest of the paper. For α 1 one can safely neglect the small real part of the self-energy (corresponding to a tiny Knight shift of the spin precession frequency ∆), so that we get from Dyson's equations (26-28) for |ω| < ω c :

ImG G 0 W 0 x (ω) - παω 1-s c ω 2 ||ω| -∆| s θ(ω 2 -∆ 2 ) ( 38
)
ImG G 0 W 0 y (ω) - παω 1-s c |ω| 2+s (ω 2 -∆ 2 ) 2 + ω 2 παω 1-s c |ω| s (39)
ImG G 0 W 0 z (ω) - παω 1-s c ∆ 2 |ω| s (ω 2 -∆ 2 ) 2 + ω 2 παω 1-s c |ω| s (40)
Following the discussion given in Sec. II D, we start by interpreting the transverse susceptibilities ImG G 0 W 0 y and ImG G 0 W 0 z . Clearly the δ(ω -∆) peak is replaced by a sharp resonance, due to the small but finite lifetime given by ImΣ y (ω ∆) in both expressions (39-40). The low-energy behavior is however different in those two quantities, as we find ImG G 0 W 0 y ∝ |ω| 2+s and ImG G 0 W 0 z ∝ |ω| s , in agreement with the NRG results, see the two topmost panels in Figs. 345and the related discussion in Sec. II D. A higher density of lowenergy modes is indeed obtained in C z than in C y , due to the coupling of the bosonic continuum to the spin along the z-axis. This leads therefore to the following decay of the z-component of the spin autocorrelation function C z (t) ∝ 1/t 1+s in the long-time limit. In the perturbative regime α 1 and α α c , it is not surprising to find very accurate agreement between the NRG and our analytical expressions (note that the real parts of the self-energies have been included in all plots). Turning to the longitudinal susceptibility ImG G 0 W 0 x , we observe in Eq. ( 38) a broad shoulder for |ω| > ∆, which is also seen in the numerics, but we obtain instead a hard gap in the range 9: Selfconsistent fermionic self-energies (analogous to the GW 0 approximation), allowing to improve the calculation in frequency ranges where lowest order perturbation incorrectly predicts the absence of spin excitations. This scheme recovers multiparticle excitations that fill in the spurious gaps of bare perturbation theory.

Σ GW 0 y = G x G 0 Φ Σ GW 0 x = G y G 0 Φ FIG.
|ω| < ∆, in disagreement with the NRG data. This surprising failure of the perturbative method for G G 0 W 0 x at weak dissipation comes from multiparticle effects beyond leading order, which do play a qualitative role whenever the spectral density is identically zero at lowest order in perturbation theory (otherwise higher order self-energy corrections are always small at weak coupling). We now examine this question in greater detail.

C. Multiparticle effects

We still consider here the spin dynamics at very weak dissipation α 1 and α α c , with the goal to improve the diagrammatics for the the lowest order longitudinal susceptibility ImG G 0 W 0 x , whose spectrum is spuriously gapped at lowest order in perturbation theory. This inconsistency can be easily tracked to expression (22), where the bare transverse propagator G 0 y resumes to a δ(|ω| -∆) peak and misses the ω 2+s low-energy modes obtained from the coupling to the bosonic bath, see equation (39). The discrepancy is thus resolved by reinjecting the full Majorana-Green's functions G i into the Majorana self-energies (this is the so-called GW 0 scheme, see Fig. 9):

ImΣ GW 0 x (ω) = α ω 0 d π ImG y ( )ImG 0 Φ (ω -)( 41
)
ImΣ GW 0 y (ω) = α ω 0 d π ImG x ( )ImG 0 Φ (ω -)(42)
The low-energy tails of the transverse susceptibility will therefore fill in the gap of the longitudinal susceptibility, so that in principle a single iteration of the above selfconsistent equation should be sufficient required. The results are shown in Fig. 10, demonstrating our correct physical interpretation of the missing low energy tails in G x . We note that the two transverse susceptibilities are very mildly affected by the self-consistency at weak coupling, and were therefore not shown. A second instance where multiparticle effects are relevant occurs in the frequency range above the high energy cutoff, now for all three spin correlations. Indeed, lowest order perturbation theory provides again a gap for excitations with |ω| > ω c , see Eqs. (36-37) in the case of one of the transverse spin correlation function (similar results are obtained for the other spin susceptibilities). In 1 and α = 0.002, comparing the raw NRG data to the self-consistent GW 0 scheme. In contrast to bare perturbation theory, which leads to a gap for ω > ωc, the self-consistent resummation correcly describes the tail arising from multiparticle excitations above the threshold ωc (note the linear frequency scale).

contrast, the NRG data show a tail with small yet nonzero spectral weight above the cutoff ω c , see Fig. 11. This continuum of magnetic excitations decays very quickly at increasing energy, and displays furthermore well-defined thresholds, pointing again to multiparticle effects, where more and more bosons are exchanged with the spin degrees of freedom, in agreement with recent observation. [START_REF] Nalbach | [END_REF] Using the fermionic self-consistent scheme proposed in Eqs. (41-42), we indeed reproduce numerically both the rapidly decaying tail and the multiparticle thresholds at quantized values of the cutoff ω c (we note that the multi-particle thresholds seen in the raw NRG data shift do not precisely match the quantized values, possibly an artifact of the discretization procedure). This final comparison concludes our study of the weak coupling regime α α c of the spin-boson model.

IV. DYNAMICS NEAR THE QUANTUM CRITICAL POINT

A. Role of the bosonic mass

We provide here general considerations on the structure of perturbation theory in terms of the coupled field theory of fermions and bosons, with the emphasis on the regime near the quantum critical point. The main motivation is to understand precisely where the singularities associated to critical modes will occur in the diagrammatics. From the knowledge gained in the previous weak coupling analysis, we understand that the fermionic sector is immune to singularities as long as the bosonic propagator remains regular. Clearly self-consistency for the fermionic self-energies (GW 0 scheme) provides only small perturbative correction and cannot account for critical behavior. One can convince oneself that the frequencydependent vertex Γ (see Fig. 7) is also regular as dissipation is increased. The only possibility left is that only the bosonic propagator encapsulates the singular dynamics via its self-energy (see general diagrammatic expression on the lower panel in Fig. 7). Indeed using the exact Dyson's equation (29) and the asymptotic bare propagator (32-33), we get the low frequency retarded bosonic Green's function:

G Φ (ω) = 1 m Φ + [a + ibSign(ω)]|ω| s -Σ Φ (ω) + Σ Φ (0) (43)
m Φ = s 4ω c -Σ Φ (0) (44)
with m Φ the bosonic mass, a and b real coefficients that can be read from Eqs. (32-33). This simple expression allows to understand the development of critical fluctuations in the spin-boson model: the quantum critical point is just associated to a vanishing of the mass m Φ in the denominator of the bosonic propagator, due to the progressive building up of the bosonic self-energy. At the quantum critical point, Σ Φ (0) = s/(4ω c ), and the correlation function now diverges as G Φ (ω) ∝ |ω| -s at low frequency, owing to expression (43). One can check that the frequency dependence of the bosonic self-energy Σ Φ (ω) is subdominant with respect to ω s in all orders in perturbation theory, so that the -s exponent of the critical fluctuations turns out to be exact. In contrast, the bosonic Green's function has a non-divergent power law behavior at low frequency G Φ (ω) -G Φ (0) ∝ |ω| s in the whole delocalized phase. This change of exponent near the quantum critical point can be used to explain the behavior of the transverse spin susceptibility, calculated

Σ G 0 W y = G 0 x G Φ Σ G 0 W x = G 0 y G Φ FIG.
12: Fermionic self-energies (without self-consistency) including a dressed bosonic propagator (this is analogous to the G 0 W approximation). This scheme is necessary to recover the critical modes near the quantum critical point, but asks for a careful evaluation of the bosonic self-energy (see Fig. 13 and the related discussion in Sec. IV B and Sec. IV D).

previously with the NRG. When using the full bosonic propagator within the fermionic self-energy (see Fig. 12), we get:

ImΣ G 0 W y (ω) = - 1 2 Sign(ω)ImG Φ (ω) (45)
∝ |ω| -s at α = α c (46)
From the exact Dyson equation for the η z Majorana Green's function Eq. ( 28), we then get at low frequency:

ImG G 0 W z (ω) 1 ∆ 2 ImΣ G 0 W y (ω) (47)
∝ |ω| -s at α = α c (48)
The divergence of the transverse spin susceptibility at the quantum critical point is thus explained (again, we stress that the above exponent is exact). Considering also the second transverse susceptibility Eq. ( 27), we get at low frequency:

ImG G 0 W y (ω) ω 2 ∆ 4 ImΣ G 0 W y (ω) (49)
∝ |ω| 2-s at α = α c (50)
as also observed with the NRG calculations of Sec. II D. This discussion shows that even in the weak dissipation regime α 1, the presence of such critical modes for α < ∼ α c will invalidate bare perturbation theory. It is thus necessary to investigate in more detail the bosonic correlations, in order to improve the agreement between numerics and diagrammatic calculation. We start by an analytic computation of the bosonic self-energy up to two-loop order.

B. Critical line at two-loop order

The lowest order (one-loop) contribution to the bosonic self-energy (first diagram in Fig. 13) is readily computed:

Σ 1loop φ (ω) = α 2∆ ∆ 2 ∆ 2 -(ω + i0 + ) 2 . ( 51
) Σ 2loop φ = G 0 x G 0 y + G 0 x G 0 y G 0 x G 0 y G 0 φ + G 0 x G 0 x G 0 x G 0 φ G 0 y + G 0 y G 0 x G 0 x G 0 x G 0 φ FIG.
13: Diagrammatic representation for the bosonic selfenergy computed up to two-loop order, leading to Eq. ( 53).

This level of approximation in the diagrammatics is equivalent to the standard RPA for the interacting electron-gas. [START_REF] Giuliani | Quantum Theory of the Electron Liquid Cambridge[END_REF] From vanishing of the mass in Eq. ( 44), the above result gives the critical dissipation at one-loop order

α 1loop c = s 2 ∆ ω c + O(s 2 ), (52)
recovering previous results. 22 We recognize here the multi-scale nature of the perturbative expansion in the dissipation strength. While the high-energy part of the spectrum is well described as long as α 1, the low energy sector will be immune to quantum critical fluctuations as long as α α c . While α c vanishes for s → 0, showing that the quantum critical point is perturbatively accessible, the above result is only accurate for very small s. As we will demonstrate now, an anomalous power-law also shows up in the dependence of α c with respect to the dimensionless parameter ∆/ω c . To establish this result, we need to push the weak-coupling calculation of the bosonic mass to next to leading order (two-loops), see Fig. 13. After straigthforward but lengthy calculations, we find at zero frequency:

Σ 2loop φ (0) = α 2∆ + α 2 4∆ 2 +∞ 0 d π ImG 0 Φ ( ) 3∆ + 2 (∆ + ) 2 α 2∆ + α 2 s ω c ∆ 2 1 -1 - s 2 ∆ ω c s (53)
where the last equation applies for ∆ < ∼ ω c . We now solve the equation s/(4ω c ) = Σ 1loop φ (0)+Σ 2loop φ (0) determining the quantum critical point, which gives in the small s limit:

α 2loop c (ω) = s 2 ∆ ω c 1 + s 2 -s log ∆ ω c + O(s 3 ). ( 54
)
This precise form confirms that the weak-coupling diagrammatics is controlled as a systematic small-s expansion. The presence of a logarithmic correction hints for the following anomalous power-law behavior:

α 2loop c (ω) s 2 1 + s 2 ∆ ω c 1-s . ( 55
)
We will first examine this analytical result with respect to the NRG phase diagram, before turning to a more rigorous derivation.

C. Phase diagram of the sub-ohmic spin-boson model

According to the above result, the formula (55) for the critical dissipation α c remains perturbative either for small s (at arbitrary ∆/ω c < ∼ 1) or for small ∆/ω c . In the latter limit, s can take relatively large values, but should not be too close to 1, otherwise the small coefficient (∆/ω c ) 1-s becomes of order 1, and perturbation theory breaks down. In this case, an alternative approach, based on a development at small h ≡ ∆/ω c but finite α can be put forward. 19,21 This results in the well-known one-loop scaling equations (which were initially obtained from the long-range unidimensional Ising model):

dα dl = (1 -s)α -αh 2 (56)
dh dl = (1 -α)h (57)
where the parameter l characterizes the change of scale during the renormalization procedure. We stress that these flow equations are valid for all values of α, as long as h 1, and are thus complementary to the regime of validity of our small α expansion at fixed h. Combining both equations, one obtains:

d dl log(α) -α + h 2 2 -(1 -s) log(h) = 0 (58)
From Eqs. ( 56) and ( 57) the renormalization fixed point occurs at α * = 1 and h * = √ 1s for s ≤ 1, so that the critical line α c (h) is obtained by the condition:

log(α c )+(1-α c )+ h 2 -(1 -s) 2 -(1-s) log h √ 1 -s = 0
(59) leading to the behavior α c ∝ h 1-s at h 1 for all s. We also recover from this analysis the exact value α c = 1 for s = 1 in the h → 0 limit. A joint comparison of the small α expansion and of the Kosterlitz flow to the NRG data is provided in Fig. 14 for the case h = 0.1. While the two expansions work where they are expected to, the twoloop small α-development is quantitatively accurate on the wide range 0 < s < 0.8, while the one-loop small hdevelopment is restricted to the small range 0.9 < s < 1. Although pushing the small h expansion to higher order would enhance the range of validity of the method, it is unlikely that this would recover the expected vanishing of α c (s/2)h at small s. In contrast, the small αexpansion converges very rapidly even for the largest s values: the extreme case s = 1 leads from Eq. ( 55 59) (giving reliable results for 0.9 < s < 1), and to the one-loop and two-loop formulas ( 52) and ( 55) (the latter is accurate on the wide range 0 < s < 0.8).

D. Ladder resummation in the bosonic self-energy

We have just checked that our analytical expression for the critical dissipation (55) is in excellent agreement with the NRG data. However, the derivation of this result was performed by a strict expansion at two-loops, and a re-exponentiation of the logarithmically singular terms. In the spirit of a numerical scheme based only on Green's functions, such as developed in Sec. III, the correct phase boundary will not be recovered if strict fourth order bosonic self-energies are used, and a clever resummation scheme should be found, that would be a diagrammatic equivalent to the re-exponentiation of the analytical result. Insight on this issue can be gained from the analytical derivation, because the logarithmic term in Eq. ( 54) can be traced back to a single diagram at two-loop order, namely the one with a single-rung ladder exchange in Fig. 13 (the last two diagrams with Majorana self-energy corrections are regular, and provide the perturbative corrections of order s 2 in Eq. ( 55)). This hints that the correct bosonic self-energy, and hence a Approximation:

≈ 1 ω ω + 1 ω + 1 + 2 1 2 ω ω 1 ω + 2 Σ ladder φ = + + + . .
FIG. 16: Upper panel: ladder series for the bosonic self-energy that is needed to obtain the correct phase boundary within the diagrammatic calculation. Lower panel: single frequency approximation used at the two-rung level (and similarly for higher order contributions), see text.

quantitative phase diagram, will be obtained upon performing a ladder resummation to all orders, similar to the diffuson mode in disordered systems, [START_REF] Castro | The Electron Liquid Paradigm in Condensed Matter Physics[END_REF] see Fig. 16.

Performing such calculation, even numerically, turns out to be difficult, because this requires the self-consistent determination of a four-electron vertex (depending of three independent external frequencies). A standard approximation in the field of disordered systems, [START_REF] Castro | The Electron Liquid Paradigm in Condensed Matter Physics[END_REF] but also in the context of interacting electron liquids [START_REF] Giuliani | Quantum Theory of the Electron Liquid Cambridge[END_REF] (called the Ng-Singwi scheme, or GW eff approximation) amounts to keeping a single running frequency within the vertex, see lower panel of Fig. 16 in the case of a two-rung ladder. In our case, this approximation is vindicated by the fact that the low frequency part of the bosonic Green's function provides the most important contribution to the Feynman diagram, so that strict energy conservation can be omitted from one rung to the next. We stress again that the inclusion of such vertex corrections in the fermionic self-energies are regular, and can be neglected at weak coupling, contrary to the singular bosonic propagator that we treat here. The ladder series including the full three particle vertex can now be explicitely carried out for the zero-temperature Matsubara bosonic self-energy:

Σ ladder Φ (iν) α dω 2π G 0 y (iω + iν)G 0 x (-iω) 1 1 -Λ 0 (iω, iν) (60) where
Λ 0 (iω, iν) α d 2π G 0 Φ (i )G 0 y (-iω -i )G 0 x (iω + iν + i ). ( 61
) Evaluation of the zero-frequency self-energy can be carried out analytically in the small ∆/ω c limit, leading to:

Σ ladder Φ (0) = α 2∆ 1 1 + 2αωc s∆ 1 -∆ ωc s . (62)
The phase boundary can now be obtained by solving Σ ladder Φ (0) = 1/G 0 Φ (0), leading to the critical value α ladder c = (s/2)(∆/ω c ) 1-s , in agreement with Eq. ( 55). Note that the small corrections at order s 2 are missing in the above derivation, because fermionic self-energy insertions were discarded in the present calculation (these were computed explicitely in Sec. IV B). This correct result for the critical dissipation confirms that the ladder resummation is the good strategy to reproduce the correct phase boundary in a fully diagrammatic calculation. We finish the paper by completing this study of the phase diagram, which will require a careful discussion of vertex corrections in the fermionic self-energies as well.

E. Full Majorana diagrammatics and the three-particle vertex

We are now equipped to collect all the needed ingredients for a successful implementation of the Majorana diagrammatics up to the quantum critical point (as long as α 1, i.e. for values of s not too close to the ohmic regime). The next step beyond the previous levels of approximation (G 0 W 0 and GW 0 schemes) is the inclusion of the vertex corrections (at ladder level) in the bosonic self-energy Eq. ( 62), which is required to obtain the correct phase boundary and to describe quantitatively the critical low-energy fluctuations. We note that the GW scheme with fully self-consistent propagators but without vertex corrections does not recover the correct critical point, and will not be further considered. We therefore propose a GWΓ approximation, which should incorporate the vertex corrections both for the bosonic and fermionic self-energies (see Fig. 7). Indeed, in order to emphasize the importance of the vertex corrections in the fermionic sector, let us first rewrite the bosonic selfenergy as Σ ladder in Fig. 17 for two values of s as a function of α < α c , which shows that the full vertex Γ remains finite up to the critical point, in agreement with our general expectation that this quantity remains regular at all orders in perturbation theory. Interestingly, the renormalization of the coupling constant α can be relatively large, up to 20% and 75% for the two cases s = 0.1 and s = 0.5 considered here. Such a sizeable effect will impact as well the fermionic self-energies, so that our consistent GWΓ scheme will finally rely on the following set of equations:

ImΣ GW Γ x (ω) = α R ω 0 d π ImG y ( )ImG Φ (ω -) (64) ImΣ GW Γ y (ω) = α R ω 0 d π ImG x ( )ImG Φ (ω -) (65) Σ GW Γ Φ (ω) = α R 2∆ ( 66
)
with α R given in Eq. ( 63). The numerical implementation of the above equation, with self-consistency imposed by Dyson's equations (26-29), is readily implemented, and gives the results displayed in Figs. 18 and19 (we note again that fermionic self-consistency is only used in order to fill in the spurious gap of C x , but it is not too important here). In the case s = 0.1, the outcome of the Majorana diagrammatics is at par with the NRG results for all three spin susceptibilities, both at low and high energies, and including the quantum critical regime. For the larger s = 0.5 value, agreement is still quite good, and the visible deviations in the quantum critical regime (for α = 0.1) arise from slight discrepancies in the critical α c between the diagrammatics and the NRG, so that the crossover scale is not accurately captured. We emphasize that this quantitative description of the spin dynamics relies of the weak coupling nature of the spin boson model for small s values. The key step here was to understand that the most important contributions to the diagrammatics are the renormalization (beyond leading order) of the bosonic mass shift and the vertex corrections. It is possible that analytical renormalization group techniques (such as developed in Ref. 57 for the quantum phase transition in the pseudogap Kondo problem at weak coupling) may be applied to the spin-boson Hamiltonian based on the present understanding, but this would likely require more complex machinery to obtain similar results. This successful comparaison seals our study of the spin dynamics in delocalized phase of the spin-boson model.

V. CONCLUSION

We have investigated the equilibrium spin dynamics of the sub-ohmic spin-boson model, which generically describes a two-level quantum system coupled to a dissipative environment of harmonic oscillators (characterized by a power-law spectrum J(ω) = 2παω 1-s c ω s θ(ω)θ(ω cω) with s ≤ 1) and submitted to a perpendicular magnetic field ∆. Two complementary techniques were employed, that were quantitatively benchmarked against each other. First, the bosonic Numerical Renormaliza- tion Group (NRG) was used with an optimized broadening method ("b-trick") to obtain reliably the dynamical spin susceptibilities from weak (α 1) to strong (α 1) dissipation. Here, we considered not only the usual transverse spin susceptibility (orthogonal to the field and parallel to the bosonic bath), but also for the first time the longitudinal (parallel to the field) and the second transverse (orthogonal to both bath and field) ones. Second, we established a general Majorana-fermion diagrammatic perturbation theory, controlled in the weak dissipation limit α 1. Because the spin-localization quantum phase transition takes place at a critical value α c which becomes perturbatively accessible (α c 1) either in the small field limit (∆ ω c for all s < 1) or in the strong sub-ohmic regime (s 1 for all ∆/ω c ), it was possible to extend the diagrammatic method up to the quantum critical point. We then examined precisely how various flavors of GW-like approximations perform against the numerically exact results, in two different regimes. For α α c , bare perturbation theory (analogous to the standard Bloch-Redfield approximation for the non-equilibrium dynamics) was sufficient to describe well the finite energy dissipative features near the applied magnetic field ∆, but fermionic self-consistency was needed in order to improve both on the low and high energy parts of the spectrum. Self-consistency was required because some of the bare self-energies can show spurious gaps that are in reality filled by multiparticle excitations. For α < ∼ α c , huge deviations appeared in the low-energy limit between lowest order perturbation theory and the NRG data, due to the proximity to the quantum critical point. We showed that the effective interaction mediated by the bosonic field was the crucial ingredient to correctly describe the physics close to criticality, and that it should include vertex corrections (at the ladder level) in the bosonic and fermionic selfenergies for quantitative agreement with the numerics. Such diagrammatic contributions are difficult to evaluate for higher dimensional field theories of coupled fermions and bosons [START_REF] Giuliani | Quantum Theory of the Electron Liquid Cambridge[END_REF][START_REF] Hertz | [END_REF] and are rarely considered (see however Ref. 58). A two-loop analysis allowed us also to obtain a simple equation for the critical line, that describes the NRG phase diagram very accurately in a wide range of parameters. We also presented a proof of Shiba's relation generalized to the sub-ohmic case and the derivation of an effective bosonic action that may support the validity of the quantum-to-classical mapping.

We finally discuss possible outlooks for the methodology and the ideas developed here. First, the addition of a magnetic field component S z parallel to the quantization axis may be considered with special interest 59 . Indeed, understanding the magnetization process is the next step towards a good description of the localized phase of the spin-boson model, which was non considered in the present work, and which has led to recent debate. [22][23][24][25][26][27][28][29] Second, further extensions of the model to multiple bosonic environments present very interesting physics in the ohmic case, 8 which remains to be investigated in greater generality. Third, an important consequence of our results is that the long time equilibrium dynamics is never accurately captured by lowest order perturbation theory in the coupling between the spin and the environment (unless dissipation is vanishingly small). This remark should have deep implications for the out-of-equilibrium long-time dynamics in dissipative models, which ought to be re-examined beyond the Bloch-Redfield regime. Extension of the Majorana diagrammatic method to two-times Keldysh Green's function could be also considered in the context of sudden quenches in the spin boson model near the quantum critical point. Lastly, considering the spin-boson Hamiltonian as a toy model for quantum criticality, the quantitative understanding achieved here may be relevant for field theories of electrons coupled to bosons in higher dimensions.

FIG. 1 :

 1 FIG.1: Sketch of the generic phase diagram of the spin-boson model, in the ohmic (s = 1) and sub-ohmic (0 < s < 1) cases respectively, as a function of the dissipation parameter α and the dimensionless magnetic field ∆/ωc. The superohmic (s > 1) regime shows only the crossover from coherent to incoherent dynamics, without a localization transition. In all cases, the coherent/incoherent crossover happens at comparable coupling strength, namely α 0.5.

  FIG. 2: (Color online) The transverse (along the bath) spin susceptibility Cz(ω) of the sub-ohmic spin-boson model at s = 0.1, ∆ = 0.1ωc and α = 0.000125. Raw data | 0, z|σz|k, z | 2 /[2Nz( k,z -0,z )] are given as circles for Nz = 20 combined NRG calculations, solid line is the perturbative result, and the three dashed lines are the various NRG broadenings discussed in the text. All frequency-dependent plots are given in units of ωc = 1.

FIG. 4 :

 4 FIG.4: Spin susceptibilities Cz(ω), Cy(ω) and Cx(ω) (top to bottom panels) for ∆/ωc = 0.1, s = 0.5 and increasing dissipation strength α = 0.02, 0.06, 0.1 (bottom to top). The critical point is located at αc 0.105. The NRG plots were generated with Nz = 10 interleaved z-averages and used the btrick discussed in Sec. II C, while the perturbative expressions (G 0 W 0 scheme) are given in Eqs.(36-37). Note that the upper panel has a different range in ω than the two other panels.

  5. Again, small values of α = 0.1 present well-

FIG. 6 :

 6 FIG. 6: Bare vertex of the spin-boson model in the Majorana language.

FIG. 10 :

 10 FIG.10: Longitudinal spin susceptibility Cx(ω) for ∆/ωc = 0.1, s = 0.1 and α = 0.001, 0.004, 0.007 (bottom to top), comparing the NRG result to the self-consistent GW 0 scheme, allowing to recover the low-energy tail below the threshold ω = ∆ (results are quantitatively correct only for the smallest α value, see Sec. IV E for the improvements that are required near the quantum critical point).

FIG. 11 :

 11 FIG.11: Transverse spin susceptibility Cz(ω) for ∆/ωc = 0.1, s = 0.1 and α = 0.002, comparing the raw NRG data to the self-consistent GW 0 scheme. In contrast to bare perturbation theory, which leads to a gap for ω > ωc, the self-consistent resummation correcly describes the tail arising from multiparticle excitations above the threshold ωc (note the linear frequency scale).

FIG. 14 :

 14 FIG.14:Phase diagram in the (s,α) plane for h ≡ ∆/ωc = 0.1, comparing the NRG data to the Kosterlitz flow obtained from the resolution of Eq. (59) (giving reliable results for 0.9 < s < 1), and to the one-loop and two-loop formulas (52) and (55) (the latter is accurate on the wide range 0 < s < 0.8).

FIG. 15 :

 15 FIG. 15: Phase diagram in the (s,α) plane for h ≡ ∆/ωc = 0.00001, 0.0001, 0.001, 0.01, 0.1 (bottom to top), comparing the NRG simulations versus the two-loop weak dissipation result (valid for s 1).

Φ 1 FIG. 17 :

 117 FIG. 17: Renormalization factor α R /α of the dissipation strength from Eq. (63) computed for ∆/ωc = 0.1 as a function of α up to the critical point (αc = 0.007, 0.105 respectively for the two values s = 0.1, 0.5).

FIG. 18 :

 18 FIG.18: Spin susceptibilities Cz(ω), Cy(ω) and Cx(ω) (top to bottom panels) for ∆/ωc = 0.1, s = 0.1 and increasing dissipation strength α = 0.001, 0.004, 0.007 (bottom to top), compared to the full Majorana diagrammatics (GWΓ scheme).

FIG. 19 :

 19 FIG.19: Spin susceptibilities Cz(ω), Cy(ω) and Cx(ω) (top to bottom panels) for ∆/ωc = 0.1, s = 0.5 and increasing dissipation strength α = 0.02, 0.06, 0.1 (bottom to top), compared to the full Majorana diagrammatics (GWΓ scheme).

Acknowledgments

We wish to thank R. Bulla, F. Evers, L. Fritz, A. Shnirman, N.-H. Tong, T. Vojta, and R. Whitney for useful discussions, and especially M. Vojta for many stimulating exchanges. We acknowledge financial support from the ERC Advanced Grant Mol-NanoSpin n226558, from MIUR under the FIRB IDEAS Project No. RBID08B3FM, and by the Universita Italo Francese/Université Franco Italienne (UIF/UFI) under the Program VINCI 2008 (Chapter II).

Appendix A: Shiba's relations for the sub-ohmic spin-boson model

The goal of this appendix is to provide a simple derivation of Shiba's relation Eq. ( 16) generalized to the subohmic case. The first step of the proof is to note that the linear coupling between the z component of the spin and the bosonic continuum, see Hamiltonian Eq. (1), implies an exact relation between the longitudinal spin susceptibility and the full retarded bosonic Green's function:

where the bare bosonic Green's function G 0 Φ (ω) was introduced in the previous equation (30). This connection can be easily established from the use of equations of motions, or equivalently by using standard manipulation of a path-integral representation of the problem. The second step lies in a low-energy analysis of the bosonic propagator G Φ (ω) from the diagrammatic perspective.

The key point is that the bosonic self-energy Σ Φ (ω), which was defined in Eq. ( 29) and expressed diagrammatically in Fig. 7, has a regular low-frequency behavior, owing to the role of the transverse magnetic field ∆ in the bare Majorana Green's function G 0 y (ω), see Eq. (23). As a consequence, all bosonic self-energy diagrams are cut in the long-time limit due to this energy gap, so that Σ Φ (ω) = Σ Φ (0) + O(ω) at low energy. Although the gap will be in reality filled by low energy excitations if one considers the full Majorana propagator G y (ω), we have seen in Sec. III B that a quickly vanishing tail of excitations is obtained, so that the above argument persists non-perturbatively (provided that perturbation theory is convergent). As a concrete illustration, equation (51) computed at one-loop order demonstrates the regular character of the bosonic self-energy in the small frequency limit. Accounting for the static part of the self-energy within the explicit Dyson equation (43), and neglecting the term Σ Φ (ω) -Σ Φ (0) at low energy, one obtains:

where the full mass m Φ was given in Eq. (44). Similarly, the bare correlation function obeys a similar behavior:

now involving the bare mass m 0 Φ . We finally use the exact relation (A1) in the low energy limit, which imposes two conditions by taking the real and imaginary parts:

A5) Solving this system to eliminate the mass m Φ leads to the simple expression:

A6)

Appendix B: Quantum-to-classical mapping

We can reformulate the usual quantum-to-classical mapping 18 in terms of an effective quantum action in imaginary time calculated perturbatively in the coupling constant α. Indeed, in the → 0 limit, we first notice that the linear form of the coupling in the Hamiltonian Eq. ( 1) leads to the exact relation:

which can be derived by the same strategy as discussed in Appendix A. This result implies that the critical behavior of the impurity spin magnetization S z can be extracted from the knowledge of the bosonic static average. The static critical properties can thus be alternatively understood from the point of view of an effective action for the boson Φ only, which is obtained by integrating out the spin degrees of freedom perturbatively in α. This can be done diagrammatically using the Majorana fermions, leading after some standard manipulations to the following expression in imaginary frequency:

. The first important parameter here is the bosonic mass m Φ = s/(4ω c )-α/(2∆)+O(α 2 ), which agrees with the one-loop result of Eq. ( 51), and controls the distance to the quantum critical point at leading order. The second crucial term is the Φ 4 interaction with coefficient u = α 2 /(16∆ 3 ) + O(α 3 ).

Within the delocalized phase, i.e. as long as the renormalized mass m R Φ does not vanish, the bosonic spectral function obeys the low-frequency behavior

e. the slow decay in imaginary time G φ (τ ) ∝ 1/τ 1+s . Owing to the exact relation (A1), we recover the correct small frequency behavior of the longitudinal spin susceptibility, see Eq. (40). The low-energy sector of the spin-boson model is thus well captured by the effective bosonic action (B2), which turns out to be equivalent to an Ising model in imaginary time with long-range spin exchange decaying as 1/(ττ ) 1+s , exactly as expected from the quantum/classical equivalence 18,19 .

Let us finally discuss the low energy behavior at the quantum critical point. Scaling analysis relative to the massless theory (m Φ = 0) leads to the bare scaling dimension [u] = 2s -1. Thus, u is relevant when 1/2 < s and irrelevant otherwise. Therefore, in the weakly subohmic regime 1/2 < s < 1, the phase transition is described by an interacting fixed point 20 and the magnetization S z obeys non-trivial classical exponents, as verified in the numerical study of Ref. 19. On the contrary, in the strongly sub-ohmic range 0 < s < 1/2, u is irrelevant and mean-field theory should apply, a result that was recently debated. [22][23][24][25][26][27][28][29] The effective action (B2) would therefore tend to strengthen the idea that the quantumto-classical mapping is robust, unless the Landau parameter u picks up non-analyticities (such as discussed in Ref. 60,61 for field theories involving coupled fermions and bosons in higher dimensions). However, the regularity of the diagrammatics would lead us to believe that this possibility is unlikely for the spin-boson model.