text:
# Using peak distribution of the cosmic microwave background for WMAP and Planck data analysis: formalism and simulations

## 1 Introduction.

The sub-degree structure of the CMB probes linear scales that were inside the horizon at the last scattering epoch. The CMB fluctuations on these scales carry a signature of causal processes during the last scattering and thereby provide a very important constraint on the physics of the early Universe and the models for structure formation. The most popular of these models is the cold dark matter (CDM) set of models based on the inflationary paradigm for the evolution of the early Universe. The models are very appealing, not only because of their relative simplicity, but also because they provide a clear-cut set of predictions that can be verified by observations. One (and perhaps the most critical) of these predictions is the sub-degree structure of the CMB anisotropies. In the framework of CDM models with adiabatic density fluctuations the structure of the CMB power spectrum reflects the linear physics of sound waves and initial density perturbations. While on scales outside the horizon at de-coupling the CMB field preserves the initial power spectrum, on sub-degree scales the interaction between the photon fluid and matter leads to a series of acoustic peaks. The relative height, width and spacing of these peaks depend on a final set of the cosmological parameters ($`\mathrm{\Omega }_{\mathrm{total}},\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_\mathrm{b},h`$) and also serve to validate the cosmological CDM paradigm (see Hu & Dodelson 2002 for a recent review). It is then important to measure the sub-degree structure of the CMB with high accuracy.

After its first year of operation, the WMAP experiment has measured the Cosmic Microwave Background (hereafter CMB) at five different frequencies with the maximal angular resolution of $`0.21\mathrm{°}`$ and sensitivity close to $`175\mu \mathrm{K}`$ per 7’ pixel, (Bennet et al. 2003a). With these sensitivity and angular resolution levels, the cosmological parameters as $`\mathrm{\Omega }_0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }_{\mathrm{baryon}}`$, $`H_0`$, $`n`$, $`\tau _{\mathrm{reio}}`$ were determined with high accuracy (Spergel et al 2003). These measurements should be further improved with the planned ESA’s mission Planck, where maps will contain up to $`10^7`$ pixels and extend to higher angular resolution. Angular resolution and noise are both important for determining the cosmological parameters from the CMB.

A major challenge to understanding current and future CMB measurements is to find an efficient algorithm that can reduce these enormous datasets: $`N10^5`$ pixels in balloon experiments, $`3\times 10^6`$ for WMAP, to $`5\times 10^7`$ for the Planck HFI data. (For comparison the COBE DMR data analysis was based on only 4,144 pixels). Traditional methods require inverting the covariance matrix and need $`N^3`$ operations making them unfeasible for the current generation of computers. Thus alternatives have been developed for estimating the CMB multipoles from Gaussian sky maps: Tegmark (1997) proposed a least variance method that yields $`C_l`$’s directly from the temperature map in O($`N^2`$) operations, while Oh, Spergel & Hinshaw (1999, hereafter OSH) perform an iterative maximization of the likelihood of the temperature map, also in O($`N^2`$) operations. The first year WMAP results were analyzed with the OSH method (Hinshaw et al 2003). Bond, Jaffe & Knox (2000) and Wandelt, Hivon & Górski (2001) concentrate on the statistics of $`C_l`$’s, once the temperature map has been Fourier transformed. Hivon, et al. (2002) studied the likelihood of the power spectrum as obtained by direct Fast Fourier Transform of the available portion of the sky. Though it requires $`O(N^{3/2}\mathrm{ln}N)`$ operations, the accuracy of their results depends on the fraction and geometry of the sky covered. A different approach consists in computing the correlation function directly from the data in $`O(N^2)`$ operations. Smoot et al. (1992) used this type of analysis for the first year release of the COBE data. More recently, Szapudi et al. (2001), Szapudi, Prunet & Colombi (2001) have developed it for mega-pixel CMB data sets. The last contribution in this field comes from Penn (2003), who applies existing (Padmanabhan, Seljak & Penn 2002) iterative algorithms on CMB data analysis, and manages to estimate to power spectra in $`𝒪`$(N log N) operations.

A different method to compute the CMB power spectrum in a fast and accurate manner was proposed by us (Kashlinsky, Hernández-Monteagudo & Atrio-Barandela, 2001, hereafter KHA). The method exploits Gaussian properties of the CMB and noise fields and uses high peaks (and troughs) of the CMB field whose abundance is much smaller than the total number of pixels, $`N`$, and whose correlation properties are strongly amplified in a way that depends on the underlying power spectrum. This simultaneously achieves two important goals: reducing the number of computational steps for analysis and the good accuracy of the measured parameters. KHA have shown that the tiny fraction of the remaining pixels (peaks and troughs) contains the required information on the CMB power spectrum in the small scales. The peaks also trace, by default, the pixels with high signal-to-noise ratio and keep most of the information about the power spectrum of the signal. Although this method may not provide as accurate results as other more time consumings methods, it provides a new, non-standard and independent tool to unveil the sub-degree structure of the CMB.

In this paper we present detailed numerical simulations with the application of the KHA method to WMAP and Planck datasets in the presence of realistic components, such as inhomogeneous noise and non-Gaussian features expected from the Galactic foreground residuals. We show that with this method we can determine the CMB power spectrum outside the beam ($`l640`$ for WMAP and higher for Planck) in only $`<10^3N^2`$ operations. We show that for the projected two-year WMAP noise levels our error bars at each $`l`$ are, on average, comparable to OSH, but have larger correlations which may be a small price to pay for the significantly shorter computational time. For both WMAP and Planck parameters the KHA method is also faster than the direct computation of the CMB power spectrum in $`O(N^{3/2}\mathrm{log}N)N^{1.61}`$ steps.

The structure of this paper is as follows: In Sect. 2 we briefly review the KHA formalism and in Sect. 3 we discuss its numerical implementation. Sect. 4 deals with application of the method to both idealized and realistic WMAP data. We show there that the method is immune to noise inhomogeneities and non-Gaussian features of Galactic foregrounds. Sect. 5 follows with results for application of the KHA method to Planck and in Sect. 6 we present our main conclusions.

## 2 Mathematical formalism: an overview

For completeness we give a brief overview of the KHA method; more details are in KHA.

The CMB sky is expected to be highly Gaussian (Komatsu et al 2003) and this property, Eq. (1) below, is also widely used in standard maximum likelihood methods. For a Gaussian ensemble of $`N`$ data points (e.g. pixels) describing the CMB data $`\delta TT`$ one will find a fraction $`f(\nu )=`$erfc($`\nu /\sqrt{2})`$ with $`|\delta |\nu \sigma `$, where $`\sigma ^2=\delta ^2`$ is the variance of the field and erfc is the complementary error function. The fraction of peaks, $`f(\nu )`$, is a rapidly decreasing function for $`\nu \genfrac{}{}{0pt}{}{_>}{^{}}1`$ and is e.g. $`f(\nu )=(4.5,1,0.1)\times 10^2`$ for $`\nu =(2,2.5,3)`$ respectively. The joint probability density of finding two pixels within $`d\delta _{1,2}`$ of $`\delta _{1,2}`$ and separated by the angular distance $`\theta `$ is given by the bivariate Gaussian:

$$p(\delta _1,\delta _2)=\frac{1}{2\pi \sqrt{𝑪}}\mathrm{exp}(\frac{1}{2}𝜹𝑪^1𝜹)$$
(1)
where $`C`$ is the covariance matrix of the temperature field. We model the temperature correlation function (or matrix, if defined on a set of pixels) as

$$C(\theta _{ij})=C^{\mathrm{CMB}}(\theta _{ij})+𝒩_i𝒩_j$$
(2)
where $`C^{\mathrm{CMB}}(\theta _{ij})`$ and $`𝒩_i𝒩_j`$ are the CMB and the noise correlation matrices, respectively. We further define the total variance of the map as $`C_0C^{\mathrm{CMB}}(0)+𝒩^2`$. The case of inhomogeneous noise, i.e., noise whose variance varies across the sky, will be discussed in Sect. 4.2. Note that the power spectrum of the map is nothing but the Legendre transform of $`C(\theta )`$, (see Eq. (5) below).

The distribution of peaks of a Gaussian field is strongly clustered (Rice 1954, Kaiser 1984, Jensen & Szalay 1986, Bardeen et al. 1986, Kashlinsky 1987). Their angular clustering can be characterized by the 2-point correlation function, $`\xi `$, describing the excess probability of finding two events at the given separation. The correlation function of such regions is:

$$\xi _\nu (\theta )=\frac{2_{\nu \sigma }^{\mathrm{}}_{\nu \sigma }^{\mathrm{}}[p(\delta _1,\delta _2)+p(\delta _1,\delta _2)]𝑑\delta _1𝑑\delta _2}{[2_{\nu \sigma }^{\mathrm{}}p(\delta )𝑑\delta ]^2}1=A_\nu (\frac{C}{C_0})$$
(3)
where $`A_\nu `$ is evaluated in detail in KHA to be:

$$A_\nu (x)=\frac{1}{H_1^2(\frac{\nu }{\sqrt{2}})}\underset{k=1}{\overset{\mathrm{}}{}}\frac{x^{2k}}{2^{2k}(2k)!}H_{2k1}^2(\frac{\nu }{\sqrt{2}})$$
(4)
Here $`H_n(x)=()^n\mathrm{exp}(x^2)d^n/dx^n\mathrm{exp}(x^2)`$ is the Hermite polynomial; $`H_1(x)\frac{\sqrt{\pi }}{2}\mathrm{exp}(x^2)`$erfc($`x)`$.

At each angular scale the value of $`\xi _\nu `$ for every $`\nu `$ is determined uniquely by $`C`$ at the same angular scale. In the limit of the entire map ($`\nu `$=0) our statistic is $`\xi _\nu `$=0 and our method becomes meaningless; the new statistic has meaning only for sufficiently high $`\nu `$. One should distinguish between the 2-point correlation function, $`\xi `$, we directly determine from the maps, and the commonly used statistics in CMB studies, the temperature correlation function, $`C`$.

KHA thus suggested the following procedure to determine the power spectrum of CMB in only $`f^2(\nu )N^2N^2`$ operations:

1. Determine the variance of the CMB temperature, $`C_0`$, from the data in $`N`$ operations;

2. Choose sufficiently high value of $`\nu `$ when $`f(\nu )`$ is small but at the same time enough pixels are left in the map for robust measurement of $`\xi _\nu (\theta )`$;

3. Determine $`\xi _\nu (\theta )`$ in $`[f(\nu )]^2N^2`$ operations.

4. Finally, solve the equation $`A_\nu (C/C_0)=\xi _\nu (\theta )`$ to obtain $`C^{\mathrm{CMB}}(\theta )`$ and from it $`C_l`$. This last step assumes that the experimental noise correlation function $`𝒩_i𝒩_j`$ has been properly determined, at least in the set of pixels selected in the analysis, so that it can be subtracted from the total correlation function $`C(\theta _{ij})`$.

Formally speaking this procedure would require $`\frac{1}{2}[f(\nu )]^2N^2`$ operations because of the symmetry in counting pairs. For simplicity we will be referring to this as O$`(f^2N^2)`$ method implying a gain factor of $`[f(\nu )]^2`$ compared to other O$`(N^2)`$ methods.

## 3 Numerical implementation

### 3.1 Modeling the CMB sky

Maps of the CMB sky were generated using the hierarchical equal area isolatitude pixelization of the celestial sphere implemented in HEALPix111HEALPix’s URL site: (http://www.eso.org/science/healpix/). Pixels have equal area and are arranged in “constant latitude rings”. Maps can be constructed with varying resolution, the number of pixels given by $`N=12\times N_{\mathrm{side}}^2`$, being $`N_{\mathrm{side}}`$ the number of times in which each side of a pixel will be divided in two, starting from a given initial configuration. For a value of $`N_{\mathrm{side}}`$ of 512, one generates a map of 3,145,728 pixels of size of seven arcminutes. In the case of PLANCK, we used $`N_{\mathrm{side}}=1024`$ or 12,582,912 pixels of 3.5 arcminutes size. The CMB spectrum was obtained using the CMBFAST code (Seljak & Zaldarriaga, 1996).

For WMAP and PLANCK, we performed two types of simulations including: (1) cosmological signal plus white noise, and (2) cosmological signal, foregrounds and inhomogeneous white noise. We chose three different thresholds: $`\nu =2.0,2.27,2.525`$ for WMAP and $`\nu =2.55,2.78,3.0`$ for Planck. At each threshold, the pixels selected are (on average) $`f(\nu )N`$ and their number doubles with decreasing threshold. The thresholds for WMAP and Planck were chosen such that they select the same number of pixels in both experiments. In this way, we can analyze the effect of pixel number on the accuracy of the method, i.e., we can estimate how the multipole error bars scale with threshold.

From the peak spatial distribution, $`\xi _\nu (\theta )`$ was computed on a grid of 31,415 equally spaced bins. The power spectrum is computed using a Gauss-Legendre integration:

$$C_l=2\pi _0^\pi 𝑑\theta \mathrm{sin}\theta C(\theta )P_l(\mathrm{cos}\theta ).$$
(5)
Accurate integration requires the evaluation of $`C(\theta )`$ at the roots of the Legendre polynomial out to a maximum order $`l_{\mathrm{max}}`$ (Press et al 1992) and this property was used by Szapudi et al (2001), Szapudi, Prunet, & Colombi (2001) and KHA. The value of $`l_{\mathrm{max}}=800`$ was chosen as a compromise between the beam scale and the pixel angular scale, in order to prevent the first root of the Legendre polynomial from being within the pixel size where there is no information on $`\xi _\nu (\theta )`$. Prior to inverting Eq. (3) $`\xi _\nu (\theta )`$ was smoothed using a Gaussian of width 4’ centered on the roots of the Legendre polynomials. We verified that other average schemes lead to essentially identical results.

Fig. (1a) shows $`\xi _{\nu =2}(\theta )`$ for a SCDM model, with $`\mathrm{\Omega }_{\mathrm{total}}=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, $`\mathrm{\Omega }_\mathrm{b}=0.04`$, $`H_0=55\mathrm{km}\mathrm{s}^1\mathrm{Mpc}`$ and scalar perturbation spectral index $`n_S=1`$. The solid line corresponds to the theoretical prediction and the dashed line is a smooth average obtained from the data. The agreement is very good out to $`\theta 10^{}`$, allowing an accurate reconstruction of the power spectrum in the range of interest: $`l\genfrac{}{}{0pt}{}{_>}{^{}}30`$.

### 3.2 Accuracy of $`\xi `$ and correlations

There are numerous ways to estimate $`\xi _\nu (\theta )`$ from the data, each way coming with its own uncertainty and systematics. Ideally, the uncertainty in $`\xi _\nu (\theta )`$ should reach $`1/\sqrt{N_{\mathrm{pairs}}}`$ but in practice systematic and other effects become dominant (Peebles 1980). Because the KHA method involves a complex and non-linear algorithm for inverting $`C_l`$’s from $`\xi _\nu (\theta )`$, the accuracy of the resultant CMB power spectrum may depend sensitively on the details, uncertainties and correlations of a particular estimator of $`\xi _\nu `$.

The ultimate accuracy with which the CMB power spectrum can be determined at each $`l`$ is given by the so called optimal variance containing both the cosmic variance and the instrument noise (Knox 1995):

$$\sigma _{\mathrm{OV}}(l)=\sqrt{\frac{2}{(2l+1)f_{\mathrm{sky}}}}[C_l+\frac{4\pi 𝒩^2B_l^2}{N}].$$
(6)
In this equation, $`N`$ is the number of pixels, $`f_{\mathrm{sky}}`$ is the fraction of the sky covered by the experiment and $`B_l`$ is the window function due to the finite beam resolution. This uncertainty also limits the accuracy with which one can determine the correlation functions $`C`$ and $`\xi `$. The latter are given by:

$$\sigma _{C(\theta )}^2=\underset{l}{}|\frac{C(\theta )}{C_l}|^2\sigma _{\mathrm{OV}}^2(l),\sigma _{\xi _\nu (\theta )}^2=|\frac{\xi _\nu (\theta )}{C(\theta )}|^2\sigma _{C(\theta )}^2.$$
(7)
The shaded region in Fig. (1) shows an example of the optimal variance uncertainty in $`\xi _{\nu =2}(\theta )`$. (Note that it is a function of $`\nu `$).

In KHA, we evaluated $`\xi _\nu `$ as the ratio of number of peak pairs at a given angular scale with respect to a poissonian catalog:

$$\widehat{\xi _\nu }(\theta )=DD/RR1.$$
(8)
At each angular separation, $`DD`$ is the number of pairs of peaks and $`RR`$ the number of pairs of $`f(\nu )N`$ points randomly located on the sphere Landy & Szalay (1993) argue that this estimator is neither optimal nor unbiased. They proposed a more accurate estimator defined by:

$$\stackrel{~}{\xi _\nu }(\theta )=(DD+RR2DR)/RR,$$
(9)
with $`DR`$ is the number of pairs given a cross-correlation between the $`f(\nu )N`$ peaks of the CMB map and the same number of points Poisson distributed on the sky. We performed simulations for WMAP as described in Sect. 4 and computed $`\xi _\nu `$ using both estimators. We found that, for the range of interest ($`\theta <10^{}`$), both estimators show scatter equally close to Eq. (7). Hence, in what follows we shall use Eq. (8) to be consistent with KHA.

In Fig. (1b) we plot the 1-$`\sigma `$ errors in $`\xi _{\nu =2}`$ for the 94 GHz WMAP channel given by the noise and cosmic variances. The smooth average obtained from the data (dashed line) is in good agreement with theoretical value. Therefore, we will expect that the radiation power spectrum obtained from $`\xi _{\nu =2}`$ as described below will be very close to that of the input model.

The sampling variance associated with the low number of pairs - compared with methods based on the correlation function - is much smaller that the cosmic variance on $`\xi _\nu `$. Fig. (2) shows $`\sqrt{N_{\mathrm{pairs}}}\sigma _{\xi ,\mathrm{OV}}`$ which is the ratio of the optimal variance uncertainty on $`\xi `$ to the statistical uncertainty, $`1/\sqrt{N_{\mathrm{pairs}}}`$. The lines correspond to various $`\nu `$=1.5 (dashes), 2 (dots) and 3 (solid line). The ratio was computed for WMAP number of pixels. The plot shows that for any mega-pixel CMB map the method can determine the new statistic, $`\xi _\nu `$, out to the angles of interest optimally even for threshold as low as $`\nu =1.5`$. This is the reason why our results are not sensitive to the particular estimator used in calculating $`\xi _\nu `$.

At $`\theta \genfrac{}{}{0pt}{}{_>}{^{}}10{}_{}{}^{}20^{}`$ the value of $`\xi _\nu (\theta )`$ is small and its value becomes dominated by shot noise. Hence, we restrict the analysis to $`\theta <10{}_{}{}^{}20^{}`$ by introducing a taper function that cuts out the contribution of the correlation function for $`\theta \genfrac{}{}{0pt}{}{_>}{^{}}10^{}`$. As a result (a) we will not recover very accurately multipoles below $`l30`$ and (b) tapering will introduce additional correlations among the different $`C_l`$’s. The first limitation is irrelevant, since one can always degrade the map to smaller resolution and apply standard techniques to recover the power spectrum at $`l2040`$. The second limitation is common to all methods that compute the power spectrum by means of the correlation function or in the presence of Galactic (and other) cut. The intrinsic multipoles need to be de-convolved from the tapering function, i.e., if $`C_l^{\mathrm{intrinsic}}`$ are the multipoles of the sky radiation power spectrum, and $`C_l`$ are obtained using Eq. (5), then

$$C_l=\underset{l^{}}{}C_l^{}^{\mathrm{intrinsic}}_{ll^{}}$$
(10)
where $`_l`$ is the Legendre transform of the tapering function $`F`$.

Fig. (3) shows the Legendre transform of the tapering function as function of $`\mathrm{\Delta }lll^{}`$ for Gaussian and top-hat tapering. The advantage of Gaussian tapering is quite obvious as it leads to significantly less prominent side-lobes. Hence, we adopted it in our computations. As the figure shows, it would lead to FWHM correlation width of only $`\mathrm{\Delta }l`$ (a few) for tapering angles $`\genfrac{}{}{0pt}{}{_<}{^{}}10^{}`$. Correlations on the $`C_l`$’s would be dominated by tapering and correlated errors on $`\xi _\nu `$.

## 4 Application to WMAP

WMAP is observing at 5 frequency bands: 23, 33, 41, 61, and 94 GHz, (Jarosik et al. 2003, Bennett et al. 2003a). The beam response at each band are given by a gain pattern, G, which is asymmetric and non-gaussian. From the solid angle beam, one can always define a FWHM, which, for increasing frequency channels, are equal to 0.82, 0.62, 0.49, 0.33 & 0.21 degrees, (Page et al. 2003). The sky maps based on the full two-years of data are expected to have an $`rms`$ noise of $``$35 $`\mu `$K per $`0.3{}_{}{}^{}\times 0.3^{}`$ pixel. By design, the noise will be essentially uncorrelated from pixel to pixel and it is expected to be highly Gaussian (Hinshaw 2000, Hinshaw et al. 2003). Due to the sky scanning strategy, the noise is reasonably uniform across the sky, except at the small regions near the ecliptic poles and at ecliptic latitude $`45^{}`$ where the sensitivity will be somewhat higher. The WMAP radiometers produce raw temperature measurements that are the differences between two points on the sky separated $`140^{}`$ . Since a given pixel $`i`$ is observed with up to 1000 different pixels $`j`$, the covariance between any given pair of pixels $`(i,j)`$ is much less than 1% of the variance of pixel $`i`$. The noise covariance in the final sky maps will, by design, be very nearly diagonal.

### 4.1 Homogeneous noise results

First, we performed a set of 500 simulations for $`\nu =(2,2.27,2.525)`$ including only cosmological signal and white noise. For each simulation, we computed the correlation function and the power spectrum. Fig. (4) summarizes the results for the threshold $`\nu =2`$. In panel (a) the dotted line shows the power spectrum recovered from a single simulation, whereas the solid line gives the input power spectrum. Tapering angle was $`\theta _t=12^{}`$. In (b), we show the average results for 500 realizations: the solid and dashed lines represent the input model and the average power spectrum of the 500 simulations, respectively. The shaded area shows the 1-$`\sigma `$ error bar region for each $`C_l`$ and the dashed lines show the optimal variance. Our method traces the acoustic peaks up to the beam scale, $`l_{\mathrm{beam}}640`$. Its accuracy is roughly similar to the optimal variance $`\sigma _{\mathrm{OV}}`$ for $`l300`$. The spectrum of Fig. (4a) has been obtained in $`10^3N^2`$ operations, using a taper window of $`12^{}`$ FWHM. The accuracy can be improved by lowering the threshold $`\nu `$ and hence increasing the computational time.

To see the effect of tapering, in Fig. (4c) we plot the power spectra obtained for two different tapering angles for the same simulation as in Fig. (4a). The solid line corresponds to $`\theta _t=8^{}`$ while the dashed line to $`\theta _t=18^{}`$. Increasing the tapering angle leads to larger oscillations, but also decreases the correlations between different $`l`$. In Fig. (4d) we plot the band power average of the previous spectra: diamonds correspond to the tapering angle of $`8^{}`$ and triangles to $`18^{}`$. The solid line is the input model and the dashed lines correspond to the optimal variance error bars as in (b). As expected, both results are almost identical since all the information at high $`l`$’s is contained at angular scales smaller than $`1^{}`$.

The correlation matrix, $`𝒞_{l,l^{}}`$, among different multipoles is given by:

$$𝒞_{l,l^{}}=\frac{\delta C_l\delta C_l^{}}{\sqrt{\delta C_l^2\delta C_l^{}^2}},$$
(11)
with $`\delta C_l=C_lC_l`$. The correlation coefficient matrix has been computed after using a Gaussian taper function $`F(\theta )=\mathrm{exp}(\theta /\theta _t)^2`$. It is highly diagonal with FWHM width $`\mathrm{\Delta }l10`$ for $`\theta _t=8^{}`$, independently  of the value of the multipole $`l`$. At $`\mathrm{\Delta }l20`$ the correlation drops to the level of $`5\%`$, independently of $`\nu `$. Outside the central diagonal strip the residual correlations are due to shot noise from the finite number of simulations. For consistency, we repeated the analysis with only 125 simulations. The width of the diagonal remained the same, but the off-diagonal terms grew by a factor of two, consistent with Poisson statistics.

For larger/smaller tapering angles, the correlation scale $`\mathrm{\Delta }l`$ will be smaller/larger but, as demonstrated in Fig. (4d), it will not change our estimate of the power spectrum. We also checked that our results did not change if we use the correlation function from the coarser (and computationally tractable) map from $`\theta =10^{}`$ out to $`180^{}`$ instead of tapering. The correlation on the final $`C_l`$’s are dominated by the correlated errors in $`\xi _\nu `$ at small angles.

In Fig. (5) we compare the accuracy of our method for WMAP and $`\nu =2`$ with that of OSH. We compute the ratio of the uncertainties in the multipoles recovered by both methods. Our method gives a comparable precision to that of OSH but requires much fewer operations. On the other hand, it gives fewer independent data points. In direct methods, such as OSH, the necessary Galactic cut leads to $`C_l`$’s that are correlated on scales of $`\mathrm{\Delta }l23`$, while our method gives a correlation length of about $`\mathrm{\Delta }l10`$.

In Fig. (6) we show the variation of the error bar with bandwidth for the multipole at $`l=200`$. Solid, dashed and dotted lines correspond to $`\nu =2.5,2.27`$ and $`2`$, respectively. Due to the presence of this correlation among the multipole estimates, this $`\sigma _{C_l}(\mathrm{\Delta }l)^{1/2}`$ behavior is obtained for $`\mathrm{\Delta }l\genfrac{}{}{0pt}{}{_>}{^{}}10`$, the correlation scale introduced by our method.

In Fig. (7) we show $`\sigma _{C_l}`$ at $`l=200,300,400,500`$ for the three thresholds. Smaller error bars correspond to smaller thresholds, i.e., to larger number of peaks. A power-law fit using all multipoles gives a power law behavior of the form $`\sigma _{C_l}N_{\mathrm{peaks}}^\beta `$, with $`\beta =0.405\pm 0.016`$. These relations can be used to estimate the amplitude of the error bar attached to each multipole for a wide range of values of $`\nu `$ and $`\mathrm{\Delta }l`$. In particular, it can be used to find what values of $`\nu `$ and $`\mathrm{\Delta }l`$ are necessary to achieve a given degree of accuracy in the power spectra and the time required in the computation.

We have shown that for a given experiment the accuracy of the $`C_l`$’s computed spectrum depends mainly on one parameter: the peak threshold allowing to choose between larger accuracy or speed.

### 4.2 Including foregrounds and realistic noise

We now generalize the method to include realistic foreground emission and inhomogeneous noise and will demonstrate that also in this case the KHA method gives similar estimates of the CMB power spectrum.

#### 4.2.1 Foregrounds

After the first year data release, the WMAP team made public foreground maps built after combining the five different frequency maps provided by the instrument, (Bennett et al. 2003b). These foreground maps were constructed after applying a Maximum Entropy Method using existing templates as priors, and the resulting model for galactic emission matches the observed emission to $`<1`$%. Using this model for Galactic emission, we get that with the fiducial value of 20% for the amplitude of the foreground residuals above $`b=10^{}`$, they are much smaller than the intrinsic CMB temperature fluctuations associated with peaks ($`\delta T\genfrac{}{}{0pt}{}{_>}{^{}}200\mu `$K for $`\nu =2`$).

To model the effect of foregrounds, we used simulated maps provided to us earlier by G. Hinshaw of the WMAP science team (2002, private communication). These were produced by combining the Haslam 408 MHz map and the Schlegel, Finkbeiner & Davis (1988, hereafter SFD) IRAS/DIRBE 100 $`\mu `$m dust map. The Haslam map was used as the template for synchrotron emission and the SFD map as the template for dust and free-free. These maps were scaled to microwave frequencies using the COBE DMR-based fits of these templates (with 7 degree resolution). The frequency by frequency fit results are in Table 1 of Kogut et al. (1996). In detail, the Haslam map was scaled using a power law index $`\alpha _{\mathrm{syn}}=3`$. The SFD free-free map was scaled using $`\alpha _{\mathrm{ff}}=2.15`$ and the SFD dust model was scaled using an index $`\alpha _{\mathrm{dust}}=2.0`$. This model is known to over-predict the plane emission at DMR resolution, thus it is likely to be conservative. In the foreground maps, we did not include point sources. Sokasian, Gawiser & Smoot (2001) have compiled and analyzed the available extra-galactic point source data and have concluded that such sources will contribute negligibly to the angular power spectrum at 94 GHz for $`l<800`$.

Fig. (8) shows the histogram of foreground contributions to the WMAP 94GHz channel for 3 values of Galactic cut: the solid, dotted and dashed lines correspond to $`|b|_{\mathrm{cut}}=5^{}`$, $`10^{}`$ and $`20^{}`$, respectively. For comparison at $`\nu =2`$ the value CMB contribution to the remaining pixels will be $`200\mu `$K. Clearly, foregrounds are not likely to significantly affect the method.

#### 4.2.2 Noise model

Our statistic $`\xi _\nu `$ has been worked out for homogeneous Gaussian fields. The inversion of $`\xi _\nu (C(\theta )/C_0)`$ into the correlation function $`C(\theta )`$ is accurate provided that signal and noise correlation functions are uniform across the sky. For most experiments this is not the case as the sky coverage is not homogeneous and the noise variance changes with location. For this more realistic case we adopt the following strategy: if the noise variance $`𝒩_i𝒩_i`$ changes according to the number of times a pixel has been observed ($`𝒩_i𝒩_it_{\mathrm{observ}}^{1/2}`$), pixel $`i`$ will be selected as peak above a threshold $`\nu `$ if it verifies $`|\delta T_i|\genfrac{}{}{0pt}{}{_>}{^=}\nu (C^{\mathrm{CMB}}(0)+𝒩_i𝒩_i)^{1/2}`$, i.e., we take into account the local noise variation. In the expression, $`C^{\mathrm{CMB}}(0)`$ is the variance of the cosmological signal. We apply this formalism to the noise model of the WMAP two-year scanning strategy. This model gives the number of observations ($`N_{\mathrm{obs}}`$) in each pixel at the end of the second mission year. Fig. (9) shows the distribution of the r.m.s. noise distribution for WMAP 94 GHz channel, for $`7`$ arcmin-sized pixels. In the figure, the noise average is $`\sqrt{𝒩_i𝒩_i_{\mathrm{full}\mathrm{sky}}}97\mu `$K per pixel.

Fig. (10) shows that this gives similar results as for homogeneous noise. The thick solid line corresponds to the theoretical estimate of $`\xi _{\nu =2}`$ according to eq (3). For a fixed cosmological signal, if we add homogeneous noise to a CMB map, the dot-dashed line will be the numerical estimate, like in Fig. (1). The thin solid line gives the estimate of $`\xi _{\nu =2}`$ assuming the noise is homogeneous when it is not, while the long-dashed line gives the estimate when the pixels have been selected according to the criterion defined in the previous paragraph. The latter result is comparable to the case of homogeneous noise and, in this sense, this procedure gives an optimal estimator.

#### 4.2.3 Results

We performed 150 simulations of WMAP data using the SCDM model, with inhomogeneous noise and foreground residuals. We choose three different galactic cuts and two different amplitudes for the foreground residuals: 10% and 20% amplitude of the original foreground data at 94GHz. In the first case, we performed simulations with a Galactic cut at $`|b|=5^{}`$ and $`10^{}`$ and for the 20% residuals we imposed a cut at $`|b|=20^{}`$. We have checked that these models give a foreground residual level comparable to the actual contribution found in the data by Bennett et al. (2003c).

In Fig. (11) we show the results of 150 simulations. As in Fig. (4), the dotted lines shows the mean power spectrum of all the simulations, while the solid line shows the input model. Shaded areas give the 1-$`\sigma `$ confidence level and dashed lines the optimal variance 1-$`\sigma `$ error bars (see Eq. (6)). In the whole range, the error bars were almost identical (less than 5% increase) to the case with no foregrounds and homogeneous noise. The correlations in $`l`$-space are practically the same as before ($`\mathrm{\Delta }l10`$). The deviations of the mean from the input model are from the smaller number of simulations.

To quantify the amplitude of the non-Gaussian components, we computed the skewness and kurtosis of our simulations. The skewness and kurtosis are defined (cf. Press et al. 1992) as $`\sigma _3=(\delta T/\sigma _2)^3`$, $`\sigma _4=(\delta T/\sigma _2)^43`$, respectively, with $`\sigma _2`$ being the variance of the data. In our case, the non-gaussian signal is dominated by the foreground residuals. We computed the skewness and kurtosis for a randomly selected map of the previous set of simulations. For a Galactic cut at $`10^{}`$ and 10% foreground residuals amplitude, the values found were $`\sigma _3=0.12`$ and $`\sigma _4=1.2`$. For comparison, for the same CMB map re-simulated with a realistic noise component, galactic cut and no foregrounds, the values were $`\sigma _3=0.02`$ and $`\sigma _4=0.03`$.

To conclude, the addition of a non-Gaussian signal (the foreground residuals) and inhomogeneous noise does not affect significantly the performance of the method. The reason that the introduction of Galactic foregrounds produces only small variations in the results is at the very core of our procedure: with the new KHA statistic, Eq. 3 we select regions with high S/N ratio, the peaks, and small departures from gaussianity neither degrades the quality of our subset of data, nor introduces additional correlations.

The results presented in Fig. (11) would improve for lower thresholds. As the method provides correlated $`C_l`$’s, it is necessary to bin the estimates into band powers. Other methods share, to some extent, this limitation; direct computation of the power spectrum from a map which necessarily has a cut due to bright foreground regions, is possible only with a finite band-width determined by size of the Galactic cut. This was avoided in the COBE/DMR data with the Gramm-Schmidt orthogonalization of the base functions (Górski 1994), but the procedure is impractical for the mega-pixel CMB datasets. For the WMAP 94 GHz channel, our method provides, under reasonable accuracy, around 12-14 independent data points for the first two Doppler peaks of the CMB power spectrum, and with a gain factor of $`10^310^4`$ in CPU time compared to standard methods.

We have also demonstrated that our method retains its accuracy in the presence of foregrounds and for the realistic/inhomogeneous WMAP noise. The treatment in the presence of the inhomogeneous noise can also be improved by removing the pixels with significantly fewer observations. Furthermore, because the peaks method computes a correlation function, $`\xi _\nu `$, it is immune to masking. The latter would allow us to reduce the foregrounds signal by removing from the CMB maps more isolated regions with higher foreground contribution.

## 5 Application to Planck

The Planck satellite222(http://astro.estec.esa.nl/Planck/) is due to be launched in early 2007 to map the all sky CMB distribution with an even finer resolution than WMAP. It will have two instruments: the low-frequency instrument (LFI), that operates at three frequency channels between 30 and 70 GHz, and the high-frequency instrument (HFI) in six frequency channels between 100 and 857 GHz. The highest Planck resolution will be 5 arcmin and the instrument noise after one year of operations is expected to be around 6-10 $`\mu `$K.

For this case we simulated the 217 GHz channel, with a beam of $`5.5^{}`$ FWHM and an average noise level of $`11.4\mu `$K per beam area. We tested this channel in two different cosmological models: for consistency, we performed 125 simulations with the same SCDM model used for WMAP, (top row of Fig. (12), but we also tested our method with the $`\mathrm{\Lambda }`$CDM concordance model, given by $`\mathrm{\Omega }_\mathrm{m}=0.292`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.708`$, $`\mathrm{\Omega }_\mathrm{b}=0.044`$ and $`h=0.72`$, (bottom row of the same figure). We assumed that the instrumental noise was the sum of two different components with white noise and $`1/f`$ contributions (see Maino et al. 1999 for a detailed account of systematic effect on the Planck LFI instrument).

The $`1/f`$ component is characterized by a knee frequency $`f_k`$, for which the power spectra of both white and $`1/f`$ noise components are equal. If the spin frequency $`f_s`$ is not much smaller than the knee frequency $`f_k`$, then we can use the fact that the telescope spends sixty minutes in each ring to simply consider the average noise pattern for every ring, (Maino et al. 2002). In other words, we can neglect the noise high frequency components present in the same ring, as these will be averaged out. We model the $`1/f`$ component as a different baseline present in every scanned ring, giving a pattern of stripes in the noise map (Janssen et al. 1996). We used the Planck 217 GHz channel noise model, provided to us by the Max Planck Institute für Astrophysik (MPA) at Garching, that includes all systematic effects. For each simulation, we added a realization of this $`1/f`$ component to the cosmological and white noise components. For simplicity we did not include foreground residuals as in the previous section we have shown to give negligible contribution. For similar reasons we also did not include point sources.

The maximum multipole to which a beam of $`5.5^{}`$ FWHM is sensitive is $`l_{\mathrm{beam}}=1472`$. The HEALPix configuration chosen to pixelize the CMB sky as seen by this channel was $`N_{\mathrm{side}}=1024`$, which yields pixels of 3.5 arcmins. In Fig. (12) we show the results of 500 simulations. As indicated, we have chosen three threshold levels that have the same number of pixels that the three $`\nu `$’s considered for WMAP. We plot the results for those thresholds following the conventions of Fig. (4). Note that, while we obtain the radiation power spectrum from the same number of pixels as in the case of WMAP, the method recovers it up to $`l=1472`$, the largest possible multipole determined by the beam size. This could be expected, since the clustering of peaks and troughs is a statistic particularly sensitive to small angular scales. In the range $`\theta [\theta _{\mathrm{beam}},10^0]`$, or equivalently $`l[40,l_{\mathrm{beam}}]`$, our method reconstructs the power spectrum with an accuracy depending only on the threshold. The range in $`l`$-space probed by the method depends on the resolution of each experiment but is independent of the total number of peaks used in the analysis.

Using the data from the three set of simulations we can compute how the band power error bars scale with the threshold $`\nu `$. In this case, a fit of the form $`\sigma _{C_l}N_{\mathrm{peaks}}^\beta `$ to the error bars of all multipoles gives: $`\beta =0.65\pm 0.011`$, close to the expected behavior $`\sigma _{C_l}1/\sqrt{N_{\mathrm{peaks}}}`$. As mentioned before, this scaling can be used to decide what threshold to choose to achieve a preselected accuracy. The correlations on the $`C_l`$’s were similar to those of WMAP, of the order of $`\mathrm{\Delta }l10`$.

## 6 Discussion and conclusions

We have presented the detailed implementation of the peaks KHA method for computing the CMB power spectrum. The method is based on the correlation properties of peaks on the CMB and assumes that temperature fluctuations are Gaussian. We have shown, that the method is robust in the presence of Galactic (non-Gaussian) foregrounds and realistic WMAP and Planck noise inhomogeneities.

This method provides a new and independent way to characterize the sub-degree structure of the CMB. Its accuracy is generally lower than standard methods, but it is remarkably faster for the range of values of $`\nu `$ required by WMAP and Planck. This is demonstrated in Fig. (13), where we compare speed of a standard $`𝒪`$(N log N) $`N^{1.61}`$ method with ours.

For WMAP the KHA method should work most optimally for peak thresholds of $`\nu 1.82.5`$ or in only $`(2.5\times 10^38\times 10^5)N^2`$ operations. Applying the method to the simulated PLANCK sky maps has shown that we could probe, with the same number of points, the power spectrum out to much higher multipoles and without significant loss of accuracy. For Planck mission parameters, the KHA method should work most optimally for $`\nu 23`$, reducing the number of steps down to $`4\times 10^6N^2`$ operations for thresholds as high as $`\nu 3`$. This can potentially lead to a gain of up to $`10^5`$ in speed compared to the existing O$`(N^2)`$ methods.

A significant advantage is that because we use a $`\xi _\nu `$ statistic, the computed CMB spectrum is immune to masking, allowing to remove high noise and foreground emission parts of the sky. This also makes straightforward the analysis of particular isolated regions in the sky, enabling its application to future CMB experiments of very high angular resolution scanning small celestial patches.

With the scalings of $`\sigma _{C_l}`$ with $`\nu `$ for both WMAP and Planck one can estimate the uncertainties of the different multipoles at still lower peak threshold values. In Fig. (14) we show what would be the spectrum recovered from the 2 year WMAP data release using a threshold $`\nu =1.6`$. The whole calculation would require only $`6\times 10^3N^2`$ operations.

###### Acknowledgements.

We are grateful to Gary Hinshaw for providing us with the Galactic foregrounds model data and the WMAP observations model used in the simulations. We also thank the Planck simulation center at the Max Planck Institute für Astrophysik for allowing us to use their 217 GHz noise model for Planck. Some of the results of this paper have been derived using the HEALPix package, (Górski, Hivon & Wandelt 1999). C.H.M. and F.A.B. acknowledge support of Junta de Castilla y León (project SA002/03) and Ministerio de Educación y Cultura (project BFM2000-1322). C.H.M acknowledges the financial support provided through the European Community’s Human Potential Programme under contract HPRN-CT-2002-00124, CMBNET. C.H.M. also thanks the Astrophysikalisches Institut Potsdam for allowing the use of their computer resources and J.A.Rubiño-Martín for encouraging the final submission of the manuscript.

REFERENCES
Bardeen, J.M., Bond, J.R., Kaiser, N. & Szalay, A.S. 1986, ApJ, 304, 15
Bennett, C., Halpern, M., Hinshaw, G. et al. 2003a, ApJ, accepted
Bennett, C., Hill, R.S., Hinshaw, G. et al. 2003b, ApJ, accepted
Bond, J.R., Jaffe, A.H. & Knox, L. 2000, ApJ, 533, 19.
Finkbeiner, D.P., Davis, M., Schlegel 1999, ApJ, 524, 867
Górski, K.M. 1994, ApJ 430, L85
Górski, K.M, Hivon, E. & Wandelt, B.D., 1999, Proceedings of the MPA/ESO Cosmology Conference ”Evolution of the Large Scale Structure”, eds. A.J.Banday, R.S.Sheth and L.Da Costa, PrintPartners Ipskamp, NL, pp.37-42.
Hinshaw, G. 2000, Proceedings of the MPA/ESO/MPE Workshop “Mining the Sky”, held at Garching, Germany, 31st July - 4th August 2000. Edited by A.J.Banday, S.Zaroubi and M.Bartelmann. Heidelberg: Springer-Verlag, 2001, p.447
Hivon, E., Górski, K.M, Netterfield, C.B. et al. 2002 ApJ, 567, 2
Hu, W. & Dodelson, S. 2002, A.R.A.A., 40, 171
Janssen, M., Scott, D., White, M. et al. 1996, astro-ph/9602009
Jarosik, N., Bennett, C.L., Halpern, M. et al., 2003, ApJS, 145, 413
Jensen, L.G. & Szalay, A.S. 1986 ApJ, 305, L5
Kaiser, N. 1984 ApJ, 282, L9
Kashlinsky, A. 1987, ApJ 317, 19
Kashlinsky, A., Hernández-Monteagudo, C. & Atrio-Barandela, F. 2001 ApJ, 557, L1 (KHA)
Knox, L., 1995, Phys. Rev. D, 52, 4307
Kogut, A., Banday, A.J., Bennett, C.L. et al. 1996, ApJ, 464, L5
Landy, S.D & Szalay, A.D. 1993 ApJ 412, 64
Maino, D., Burinaga, C., Maltoni, M. et al. 1999, A&AS, 140, 383
Maino, D., Burinaga, C., Górski, K.M, Mandolesi, N. & Bersanelli, M. 2002, A&A, 387, 356
Oh, S.P., Spergel, D.N. & Hinshaw, G. 1999 ApJ, 510, 551
Padmanabhan, N., Seljak, U. & Penn, U.L. 2003 NewA 8 581
Page, L., Nolta, N.R, Barnes, C. et al., 2003, ApJ, accepted
Peebles, P.J.E. 1980 “The Large Scale Structure of the Universe”, Princeton University Press, Princeton
Penn, U., 2003, astro-ph/0304513
Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P. 1992 “Numerical Recipes”, 2nd edition, Cambridge University Press, Cambridge
Rice, S.O. 1954, in “Noise and Stochastic Processes”, ed. Wax, N., p.133 Dover (NY)
Seljak, U. & Zaldarriaga, M. 1996, Ap. J., 469, 437 Smoot, G.F., Bennett, C.L., Kogut, A. et al. 1992, ApJ, 396, L1
Sokasian, A., Gawiser, E. & Smoot G.F. 2001, ApJ 562 885
Szapudi, I., Prunet, S., Pogosyan, D. et al 2001, ApJ, 548, L115
Szapudi, I., Prunet, S. & Colombi, S. 2001, ApJ, 561, L11
Schlegel, D.J., Finkbeiner, D.P. & Davis, M. 1998, ApJ, 500, 525
Tegmark, M. 1997, Phys.Rev.D., 55, 5898
Wandelt, B.D., Hivon, E. & Górski, K. 2001, Phys.Rev.D, 64, 083003


format: text

text:
# Critical local moment fluctuations in the Bose-Fermi Kondo model

## I Introduction

The interest in quantum criticality in strongly correlated metals arises primarily because it provides a mechanism for non-Fermi liquid behavior. While the problem appears to be important for a number of correlated electron systems – including high temperature superconductors – the issues are particularly well-defined in heavy fermion metals. Here, many materials have been shown to explicitly display a magnetic quantum critical point (QCP), including a growing list of stoichiometric or nearly stoichiometric ones . A particularly striking puzzle has emerged from inelastic neutron-scattering experiments and magnetization measurements in some of these (nearly) stoichiometric materials. The dynamical spin susceptibility displays an $`\omega /T`$ scaling with a fractional energy/temperature exponent. In addition, the same exponent is seen not only at the antiferromagnetic ordering wavevector but also essentially everywhere else in the Brillouin zone.

In a recent work, a “locally-critical point” is found in Kondo lattice systems. This picture appears to explain the salient features of the aforementioned experiments. At a locally critical point, spatially extended critical fluctuations co-exist with spatially local ones. More microscopically, the criticality of the local Kondo physics is “embedded” into that associated with the magnetic ordering of the lattice: the divergence of the spatial correlation length is accompanied by the destruction of the Kondo resonance. The microscopic analysis was carried out within an extended dynamical mean field theory (EDMFT), in which the Kondo lattice system is treated in terms of a single-impurity Bose-Fermi Kondo model supplemented by a self-consistency condition. The impurity model describes a local moment coupled at once to a conduction electron band and a dissipative bath of vector bosons; the bosonic bath describes a fluctuating magnetic field generated by the neighboring local moments. The locally critical point is identified as a self-consistent solution when the Bose-Fermi Kondo model is treated to the first order in $`ϵ`$ within an $`ϵ`$expansion renormalization group (RG) procedure. Here $`ϵ1\gamma `$, where $`\gamma `$ is the power-law exponent of the spectrum of the dissipative bosonic bath (defined in Eq. (4)).

There are several important questions that remain. First, does a self-consistent locally-critical solution arise at higher orders in $`ϵ`$? Second, what happens to the locally critical point in spin-anisotropic situations? Such anisotropy occurs in heavy fermions since the spin-orbit coupling is usually strong in these systems.

These issues are addressed in the present paper. We carry out a detailed analysis of the Bose-Fermi Kondo model and show that an unstable fixed point exists not only in the spin-isotropic case but also in the presence of spin-anisotropy (both xy and Ising cases). In each case, this critical point describes a continuous transition from a Kondo phase, in which the local moment is quenched by the spins of the conduction electrons, to a “local-moment” phase where there is no Kondo resonance. For the isotropic and xy cases, we calculate the local spin susceptibility to the order $`ϵ^2`$ and, in addition, we determine the associated critical exponent, $`\eta `$ (defined in Eq. (38)), to all orders in $`ϵ`$. In the Ising case, a slightly different RG approach turns out to be useful. In all three cases, we find that the exponent $`\eta `$ at the unstable fixed point is

$`\eta =ϵ`$ (1)
This result turns out to guarantee the existence of a self-consistent solution of the locally-critical type in the Kondo lattice systems.

We note in passing that the Bose-Fermi Kondo model is also of interest in other contexts. Historically, the model, in its Ising version, was first studied in the context of an EDMFT treatment of a spinless model in ref. , using an $`ϵ`$expansion. The spinful version was subsequently considered, also through an $`ϵ`$expansion, in refs. . The Bose-only Kondo model was later extensively analyzed within a similar $`ϵ`$expansion in ref. , which went to higher orders in $`ϵ`$ and established Eq. (1) for a stable fixed point of that model to all orders. The Bose-only Kondo model appears implicitly in a single-impurity spin model first studied (using a large N approach) in ref. ; see also ref. . The Fermi-only Kondo model is of course a standard textbook problem.

The remainder of the paper is organized as follows. In Section II, we define the model and introduce the formalism. Section III is devoted to the RG analysis, to order $`ϵ^2`$, in the isotropic case. Section IV presents the explicit calculation of the local spin susceptibility to order $`ϵ^2`$, as well as the analysis of the associated critical exponent $`\eta `$ to all orders. In Section V, we consider the effect of spin-anisotropy. Section VI discusses the consequence of our results for the locally critical solution in a Kondo lattice and Section VII provides a brief summary of our results. Some of the technical details are relegated to Appendices A, B C, and D.

## II The model and formalism

The Bose-Fermi Kondo model is defined as follows

$`_{\text{BFK}}=`$ $`J𝐒𝐬_c+{\displaystyle \underset{p,\sigma }{}}E_pc_{p\sigma }^{}c_{p\sigma }`$ (3)
$`+g{\displaystyle \underset{p}{}}𝐒\left(\stackrel{}{\varphi }_p+\stackrel{}{\varphi }_p^{}\right)+{\displaystyle \underset{p}{}}w_p\stackrel{}{\varphi }_p^{}\stackrel{}{\varphi }_p.`$
A spin-$`\frac{1}{2}`$ local moment, $`𝐒`$, is coupled to both a fermionic bath ($`c_{p\sigma }`$), through the Kondo interaction $`J`$, and a dissipative vector-bosonic bath ($`\stackrel{}{\varphi }_p`$) with a coupling constant $`g`$. $`J`$ is positive, i.e., antiferromagnetic, but $`g`$ can be either positive or negative. (A sign change in $`g`$ can be absorbed by a corresponding sign change in $`\varphi `$.) The spectral function of the bosonic bath is taken to have a sublinear power-law dependence on energy, at sufficiently low energies:

$`{\displaystyle \underset{p}{}}[\delta (\omega w_p)\delta (\omega +w_p)]=(K_0^2/\pi )|\omega |^\gamma \mathrm{sgn}\omega `$ (4)
for $`|\omega |<\mathrm{\Lambda }`$. Here, the power-law exponent

$`0<\gamma 1ϵ<1.`$ (5)
The density of states of the conduction electron band near the Fermi energy is taken to be a constant:

$$\underset{p}{}\delta (\omega E_p)=N_0.$$
(6)
The fluctuating magnetic field competes against the Kondo-singlet formation. When the fluctuations are sufficiently slow such that $`ϵ>0`$, $`g`$ is a relevant coupling in the RG sense and it can compete with the marginally-relevant $`J`$ coupling. The physics of this competition is amenable to an $`ϵ`$expansion. This is in contrast to the Kondo fixed point, which occurs at an infinite coupling and can only be addressed by strong-coupling methods.

Following Smith and Si, we adopt the Abrikosov representation of the spin in terms of pseudo-$`f`$-electrons,

$$𝐒=\underset{\sigma \sigma ^{}}{}f_\sigma ^{}\frac{\stackrel{}{\tau }_{\sigma \sigma ^{}}}{2}f_\sigma ^{},$$
(7)
where $`\tau ^{x,y,z}`$ are the Pauli matrices. In this representation, the Bose-Fermi Kondo Hamiltonian takes the following form,

$`\stackrel{~}{}`$ $`=_0+_J+_g,`$ (8)
$`_0`$ $`={\displaystyle \underset{\sigma }{}}\lambda f_\sigma ^{}f_\sigma +{\displaystyle \underset{p,\sigma }{}}E_pc_{p\sigma }^{}c_{p\sigma }+{\displaystyle \underset{p}{}}w_p\stackrel{}{\varphi }_p^{}\stackrel{}{\varphi }_p,`$ (9)
$`_J`$ $`={\displaystyle \frac{J}{4}}{\displaystyle \underset{\sigma \sigma ^{}}{}}\sigma \sigma ^{}f_\sigma ^{}f_\sigma c_\sigma ^{}^{}c_\sigma ^{}+{\displaystyle \frac{J}{2}}(f_{}^{}f_{}c_{}^{}c_{}+\text{H.c.}),`$ (11)
$`_g`$ $`={\displaystyle \frac{g}{2}}{\displaystyle \underset{\sigma }{}}\sigma f_\sigma ^{}f_\sigma \varphi ^z+{\displaystyle \frac{g}{\sqrt{2}}}(f_{}^{}f_{}\varphi ^{}+\text{H.c.}),`$ (12)
where $`\lambda `$ is an energy level for the $`f`$electron that will be set to $`\mathrm{}`$ at the end, $`\sigma =\pm 1`$, $`\stackrel{}{\varphi }_p(\stackrel{}{\varphi }_p+\stackrel{}{\varphi }_p^{})`$, and $`\varphi ^\pm =(\varphi ^x\pm \varphi ^y)/\sqrt{2}`$.

To analyze the critical behavior of the Bose-Fermi Kondo model, we carry out an RG procedure based on a dimensional regularization with a minimal subtraction (MS) of poles. We define a renormalized field $`f`$ and a dimensionless coupling constant $`g`$ by

$`f_B`$ $`=`$ $`Z_f^{1/2}f,`$ (13)
$`g_B`$ $`=`$ $`gZ_f^1Z_g\mu ^{ϵ/2},`$ (14)
where $`Z_f`$ is the wave-function renormalization factor for $`f`$ electrons, $`g_B`$ is the bare coupling constant, $`Z_g`$ is a coupling constant renormalization of $`g`$, and $`\mu `$ is a renormalization energy scale.

The “dimensional” regularization for the conduction electron density of states is done by introducing an $`ϵ^{}`$ through

$$\underset{p}{}\delta (\omega E_p)=N_0|\omega |^ϵ^{},$$
(15)
Accordingly, the dimensionless coupling constant $`J`$ is defined by

$$J_B=JZ_f^1Z_J\mu ^ϵ^{}.$$
(16)
Note that $`ϵ^{}`$ is introduced strictly for the purpose of facilitating the MS analysis; it is set to zero at the final stage of the calculation.

## III RG analysis

In this section, we carry out the RG analysis to order $`ϵ^2`$.

### A RG flow equations

It is known that, to the linear order of $`ϵ`$, there is a critical point at $`(K_0g^{})^2=ϵ/2`$, $`(N_0J^{})=ϵ/2`$. To obtain the corrections to the order $`ϵ^2`$, we need to calculate the self-energy and vertex corrections beyond the orders $`J`$ and $`g^2`$ and include also terms to the orders $`g^4`$, $`J^2`$, and $`g^2J`$. The details of our calculation are given in Appendix A.

The renormalization factors are obtained as

$`Z_f=`$ $`1{\displaystyle \frac{3}{4ϵ}}(K_0g)^2{\displaystyle \frac{3}{16ϵ^{}}}(N_0J)^2`$ (18)
$`{\displaystyle \frac{15}{32ϵ^2}}(K_0g)^4+{\displaystyle \frac{3}{8ϵ}}(K_0g)^4`$
$`Z_g=`$ $`1+{\displaystyle \frac{1}{4ϵ}}(K_0g)^2+{\displaystyle \frac{1}{16ϵ^{}}}(N_0J)^2`$ (20)
$`+{\displaystyle \frac{9}{32ϵ^2}}(K_0g)^4{\displaystyle \frac{1}{8ϵ}}(K_0g)^4,`$
$`Z_J=`$ $`1{\displaystyle \frac{1}{ϵ^{}}}N_0J`$ (24)
$`+{\displaystyle \frac{1}{4ϵ}}(K_0g)^2+{\displaystyle \frac{1}{16ϵ^{}}}(N_0J)^2+{\displaystyle \frac{1}{ϵ^2}}(N_0J)^2`$
$`\left({\displaystyle \frac{1}{ϵ^{}(ϵ+ϵ^{})}}+{\displaystyle \frac{1}{4ϵϵ^{}}}\right)(K_0g)^2(N_0J)`$
$`+{\displaystyle \frac{9}{32ϵ^2}}(K_0g)^4{\displaystyle \frac{1}{8ϵ}}(K_0g)^4.`$
They can then be used to determine the beta functions for the coupling constants $`g`$ and $`J`$, defined as

$`\beta (g)=`$ $`\mu {\displaystyle \frac{dg}{d\mu }}|_{g_B,J_B},`$ (25)
$`\beta (J)=`$ $`\mu {\displaystyle \frac{dJ}{d\mu }}|_{g_B,J_B}.`$ (26)
The results are given as follows

$`\beta (g)=`$ $`g\left({\displaystyle \frac{ϵ}{2}}(K_0g)^2+(K_0g)^4{\displaystyle \frac{(N_0J)^2}{2}}\right),`$ (27)
$`\beta (J)=`$ $`J\left((N_0J){\displaystyle \frac{(N_0J)^2}{2}}\right)`$ (29)
$`J\left((K_0g)^2+(K_0g)^4\right).`$
### B The phase diagram

The RG equations (29) yield several fixed points. There are two stable ones. The fixed point located at $`g^{}=0`$ and large $`J^{}`$ (“K” in Fig. 1), as usual, specifies the Kondo phase. Here, the local moment is quenched by the spins of the conduction electrons, leading to the development of a Kondo resonance.

Another stable fixed point (“L” in Fig. 1) is located at

$`(K_0g^{})^2`$ $`=`$ $`{\displaystyle \frac{ϵ}{2}}+{\displaystyle \frac{1}{4}}ϵ^2+O(ϵ^3)`$ (30)
$`(N_0J^{})`$ $`=`$ $`0`$ (31)
We will call it a “local-moment” fixed point, emphasizing the fact that it describes a phase in which no Kondo resonance arises. The dynamics in the local-moment phase is controlled by the coupling of the local moment to the dissipative bosonic bath alone.

A separatrix specifies the boundary of the domains of attractions for these two phases in the J-g parameter space. Lying on this separatrix is an unstable fixed point, or a critical point (“C” in Fig. 1). It is located at

$`(K_0g^{})^2`$ $`=`$ $`{\displaystyle \frac{ϵ}{2}}+{\displaystyle \frac{1}{8}}ϵ^2+O(ϵ^3)`$ (32)
$`(N_0J^{})`$ $`=`$ $`{\displaystyle \frac{ϵ}{2}}+O(ϵ^3)`$ (33)
The critical point marks the point where the spectral weight of the Kondo resonance goes to zero. It captures the competition between the Kondo coupling of the local moment to the conduction electrons on the one hand, and its coupling to the fluctuating magnetic field on the other.

The RG flow can be solved numerically. In Fig. 1, we give the result for $`ϵ=0.1`$.

## IV The local dynamical spin susceptibility

We now turn to the calculation of the local susceptibility.

### A Critical susceptibility to order $`ϵ^2`$

Consider the local spin susceptibility to order $`ϵ^2`$. We first calculate the bare spin susceptibility, i.e., the autocorrelation function of the unrenormalized local spin operator specified by Eq. (7),

$`\chi _B(\tau )`$ $`{\displaystyle \frac{1}{2}}T_\tau S^{}(\tau )S^+(0)=\underset{\lambda \mathrm{}}{lim}{\displaystyle \frac{1}{2}}e^{\beta \lambda }\stackrel{~}{\chi }_B(\tau ),`$ (34)
$`\stackrel{~}{\chi }_B(\tau )`$ $`{\displaystyle \frac{1}{2}}T_\tau f_{}^{}(\tau )f_{}(\tau )f_{}^{}(0)f_{}(0),`$ (35)
where $`S^\pm (S^x\pm iS^y)`$. (For notational simplicity, we’ve dropped the subscript $`B`$ in the operators used in Eq. (35).)

Correspondingly, we define a spin renormalization factor $`Z`$ which relates the renormalized spin susceptibility, $`\chi `$, to the bare susceptibility,

$$\chi (\tau )=\frac{1}{Z}\chi _B(\tau ),$$
(36)
The calculational details are given in Appendix B. From Eq. (B12), $`Z`$ is given by

$$Z=1\frac{2(K_0g)^2}{ϵ}+\frac{(K_0g)^4}{ϵ}\frac{(N_0J)^2}{2ϵ^{}}.$$
(37)
At the critical point, the local spin susceptibility is expected to have a power-law form:

$`\chi (\tau )A_ϵ\left({\displaystyle \frac{\tau _0}{|\tau |}}\right)^\eta ,`$ (38)
for $`|\tau |\tau _0`$, where $`\tau _0=1/\mathrm{\Lambda }`$ is a cut-off scale. Here, $`\eta `$ is the anomalous dimension. It can be calculated from the renormalization factor $`Z`$,

$`\eta =`$ $`\mu {\displaystyle \frac{d\mathrm{ln}Z}{d\mu }}|_{g=g^{},J=J^{}}`$ (39)
$`=`$ $`\left(\beta (g){\displaystyle \frac{\mathrm{ln}Z}{g}}+\beta (J){\displaystyle \frac{\mathrm{ln}Z}{J}}\right)|_{g=g^{},J=J^{}}`$ (40)
$`=`$ $`2(K_0g^{})^22(K_0g^{})^4+(N_0J^{})^2`$ (41)
$`=`$ $`ϵ.`$ (42)
It is striking that, although the order $`ϵ^2`$ corrections are present both in the beta functions and in the critical coupling constants, they cancel out in the critical exponent $`\eta `$, leaving only a non-zero linear term.

The critical amplitude is given in Appendix B.

For completeness, we also briefly discuss the local spin susceptibility at the local-moment fixed point. We assume that the susceptibility here also has a power-law form:

$$\chi (\tau )A_ϵ^L\left(\frac{\tau _0}{|\tau |}\right)^{\eta _L},$$
(43)
with an anomalous dimension $`\eta _L`$. $`\eta _L`$ has the same dependence on $`g^{},J^{}`$ as in Eq. (42). Even though the coupling constants $`g^{},J^{}`$ now take the values specified by Eq. (31), it turns out that the $`ϵ^2`$ terms again cancel with each other leaving

$$\eta _L=ϵ.$$
(44)
The critical amplitude is also given in Appendix B.

### B Critical exponent to all orders in $`ϵ`$

In this section, we show that $`\eta =ϵ`$ given in Eq. (42) is in fact exact to all orders.

First, we explore the reason of the cancellation of the $`ϵ^2`$ terms in $`\eta `$ seen in Eq. (42). Using the equations (24), (37), we can easily verify that

$$Z^1=(Z_f^1Z_g)^2$$
(45)
both to linear and second orders in $`ϵ`$. Combining Eq. (45) with Eq. (14) lead to

$$g_B=gZ^{1/2}\mu ^{ϵ/2}.$$
(46)
Differentiating the logarithm of both sides with respect to $`\mathrm{ln}\mu `$, keeping $`g_B`$ and $`J_B`$ fixed, we have

$$\eta =ϵ+\frac{2\beta (g)}{g}|_{g=g^{},J=J^{}}.$$
(47)
At the fixed point, $`\beta (g)=0`$, which gives

$$\eta =ϵ.$$
(48)
There is also a consistency check for our result. We can start from Eq. (16). Combining this equation with Eq. (45) lead to

$`J_B=JZ^{1/2}(Z_J/Z_g).`$ (49)
(We have set $`ϵ^{}=0`$.) Again, differentiating the logarithm of both sides with respect to $`\mathrm{ln}\mu `$, and again keeping $`g_B`$ and $`J_B`$ fixed, we end up with

$$\eta =2\frac{d\mathrm{ln}(Z_J/Z_g)}{d\mathrm{ln}\mu }|_{g=g^{},J=J^{}}.$$
(50)
To be consistent, the RHS should be equal to $`ϵ`$. This is indeed verified to both linear and quadratic orders in $`ϵ`$.

We now proceed to higher order contributions. We show in Appendix B that Eq. (45) is valid to all orders. Eq. (47) is then valid to all orders in $`ϵ`$. By extension, $`\eta `$ equals to $`ϵ`$ to all orders in $`ϵ`$. Our reasoning basically parallels that of ref.  for the stable fixed point in the Bose-only Kondo problem.

## V Anisotropic cases (xy and Ising)

So far, we have discussed the spin-isotropic case. We now turn to the effect of anisotropy in spin-space.

### A RG equations

We introduce separate parameters for the transverse and longitudinal spin couplings as follows:

$`_J`$ $`={\displaystyle \frac{J_z}{4}}{\displaystyle \underset{\sigma \sigma ^{}}{}}\sigma \sigma ^{}f_\sigma ^{}f_\sigma c_\sigma ^{}^{}c_\sigma ^{}+{\displaystyle \frac{J_{}}{2}}(f_{}^{}f_{}c_{}^{}c_{}+\text{H.c.}),`$ (51)
$`_g`$ $`={\displaystyle \frac{g_z}{2}}{\displaystyle \underset{\sigma }{}}\sigma f_\sigma ^{}f_\sigma \varphi ^z+{\displaystyle \frac{g_{}}{\sqrt{2}}}(f_{}^{}f_{}\varphi ^{}+\text{H.c.}),`$ (53)
The RG procedure parallels that for the isotropic case, except that we need to differentiate between the longitudinal and transverse couplings in all the contributions, and associate different renormalization factors for the different types of vertices:

$`g_i^B`$ $`=`$ $`g_iZ_f^1Z_{g_i}\mu ^{ϵ/2}`$ (54)
$`J_i^B`$ $`=`$ $`J_iZ_f^1Z_{J_i}\mu ^ϵ^{},`$ (55)
where $`i=,z`$. Up to the order of our interest, the results for the renormalization factors are given in Appendix C. The resulting RG equations are given as follows,

$`\beta (g_{})`$ $`=`$ $`{\displaystyle \frac{ϵ}{2}}g_{}+{\displaystyle \frac{(K_0g_{})^2+(K_0g_z)^2}{2}}g_{}`$ (56)
$``$ $`{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2+(K_0g_{})^4}{2}}g_{}`$ (57)
$`+`$ $`{\displaystyle \frac{(N_0J_{})^2+(N_0J_z)^2}{4}}g_{}`$ (58)
$`\beta (g_z)`$ $`=`$ $`{\displaystyle \frac{ϵ}{2}}g_z+(K_0g_{})^2g_z`$ (59)
$``$ $`(K_0g_{})^2(K_0g_z)^2g_z`$ (60)
$`+`$ $`{\displaystyle \frac{(N_0J_{})^2}{2}}g_z`$ (61)
$`\beta (J_{})`$ $`=`$ $`{\displaystyle \frac{(K_0g_{})^2+(K_0g_z)^2}{2}}J_{}`$ (62)
$``$ $`{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2+(K_0g_{})^4}{2}}J_{}`$ (63)
$``$ $`(N_0J_z)J_{}+{\displaystyle \frac{(N_0J_{})^2+(N_0J_z)^2}{4}}J_{}`$ (64)
$`\beta (J_z)`$ $`=`$ $`(K_0g_{})^2J_z`$ (65)
$``$ $`(K_0g_{})^2(K_0g_z)^2J_z`$ (66)
$``$ $`(N_0J_{})J_{}+{\displaystyle \frac{(N_0J_{})^2}{2}}J_z`$ (67)
Consider first the $`J=0`$ case, i.e. the Bose-only Kondo model. The RG flow is given in Fig. 2. There are two more fixed points in addition to the isotropic one. The xy fixed point is located at

$`(K_0g_{}^{})^2`$ $`=ϵ+ϵ^2+O(ϵ^3)`$ (68)
$`(K_0g_z^{})`$ $`=0.`$ (69)
While the Ising fixed point is nominally located at $`g_{}^{}=0`$ and $`g_z^{}=\mathrm{}`$. Both the isotropic fixed point and the xy fixed point are accessible by the $`ϵ`$expansion. The Ising fixed point, on the other hand, is beyond the reach of the perturbative RG scheme.

The $`J`$coupling can lead away from the “local-moment” fixed point not only in the isotropic case, but also in the xy and Ising cases. The isotropic case has already been discussed in the previous sections. We now turn to the xy and Ising cases, respectively.

### B Critical behavior in the xy case

We will now set $`g_z=0`$ in Eqs. (58)- (67). \[Eq. (61) implies that $`g_z`$ will stay at zero under the RG transformation when its initial value vanishes.\] An unstable fixed point occurs at

$`(K_0g_{}^{})^2`$ $`=ϵ+{\displaystyle \frac{5ϵ^2}{8}}+O(ϵ^3)`$ (70)
$`N_0J_{}^{}`$ $`={\displaystyle \frac{ϵ}{\sqrt{2}}}+{\displaystyle \frac{7}{16\sqrt{2}}}ϵ^2+O(ϵ^3)`$ (71)
$`N_0J_z^{}`$ $`={\displaystyle \frac{ϵ}{2}}+O(ϵ^3)`$ (72)
Note that $`g_z`$ remains irrelevant near this fixed point, establishing the consistency of our analysis. The RG flows are shown in Fig. 3.

The perturbative correction to the local spin susceptibility can be calculated as in the isotropic case. The most general expression for the susceptibility as well as its renormalization factor $`Z`$ are given in Appendix C. The corresponding expression for $`\eta `$ is also given there, in Eq. (C34). Using the values of the coupling constants at this fixed point, Eq. (72), we find that $`\eta `$ remains to be $`ϵ`$ to order $`ϵ^2`$.

The same argument made in Appendix B for the isotropic case carries through here for the xy-case as well, resulting in

$$Z^1=(Z_f^1Z_g_{})^2,$$
(73)
to infinite orders in $`ϵ`$. So $`\eta =ϵ`$ is again valid to all orders in $`ϵ`$.

We now again briefly examine the local spin susceptibility at the local-moment fixed point in this case, where $`J_{}^{}=J_z^{}=0`$, and

$$(K_0g_{}^{})^2=ϵ+ϵ^2+O(ϵ^3).$$
(74)
The anomalous dimension $`\eta _L`$ has the same expression as in Eq. (C41). Substituting the values of the stable fixed point, we find that $`\eta _L`$ is equal to $`ϵ`$, just as in the isotropic case.

### C Critical behavior in the Ising case

None of the non-trivial fixed points in the Ising case is within the reach of the perturbative RG scheme. For the unstable fixed point, for instance, setting $`g_{}=0`$ in Eqs. (58-67) will yield a $`(N_0J_{}^{})^2ϵ`$ that is still small, but $`J_z^{}`$ and $`g_z^{}`$ that are of order unity. The latter violates our starting assumption.

This problem, however, has already been studied in ref. . A finite $`J_z`$ and $`g_z`$ can be handled by introducing a so-called “kink-gas” representation. In this representation, the unstable fixed point is still accessible by an $`ϵ`$expansion. The calculation is outlined in Appendix D. The RG equations are no longer constructed in terms of the bare couplings. Instead, they are given in terms of the stiffness constants $`\kappa _j`$, $`\kappa _g`$, and the fugacity $`y_j`$, whose initial values are specified by the bare parameters as follows,

$`\kappa _j^B`$ $`=[1{\displaystyle \frac{1}{\pi }}\mathrm{tan}^1({\displaystyle \frac{\pi }{4}}N_0J_z^B)]^2`$ (75)
$`\kappa _g^B`$ $`={\displaystyle \frac{\mathrm{\Gamma }(\gamma )}{4}}\tau _0^{1\gamma }(K_0g_z^B)^2`$ (76)
$`y_j^B`$ $`={\displaystyle \frac{N_0J_{}^B}{2}}`$ (77)
The RG flow projected onto the $`y_j\kappa _g`$ plane is shown in Fig. 4. In the vicinity of the unstable fixed point, $`\kappa _j`$ is irrelevant.

The exponent $`\eta `$ at the unstable fixed point is calculated in Appendix D to the linear order in $`ϵ`$. Again, it is found that

$`\eta =ϵ`$ (78)
In the same Appendix, it is also shown that $`\eta `$ is in fact the same as that for the critical point of a classical ferromagnetic Ising chain with a long-range interaction that decays in distance in terms of a power-law exponent $`2ϵ`$. For this problem, it has long been held that Eq. (78) is in fact exact.

In passing, we note that this analogy also allows us to state what happens to the local susceptibility at the local-moment fixed point shown in Fig. (4). Here, $`\chi (\tau )`$ picks up a constant ($`\tau `$independent) piece. At the same time, the “connected” susceptibility (defined such that the constant piece does not appear) decays with an exponent $`\eta _L=2ϵ`$.

## VI Locally critical point of the Kondo lattice

We now discuss the consequences of our results for the locally critical point of the Kondo lattices. A self-consistent microscopic treatment of the Kondo lattice model has been presented in detail elsewhere. This was carried out within the EDMFT approach developed in. (This approach is a generalization of the dynamical mean field theory such that magnetic fluctuations are also taken into account.) The self-consistent equation reads as follows,

$`<\chi (𝐪,\omega )>_𝐪=\chi _{\mathrm{loc}}(\omega ),`$ (79)
Here, the lattice dynamical spin susceptibility has the form

$`\chi (𝐪,\omega )={\displaystyle \frac{1}{M(\omega )+I_𝐪}},`$ (80)
where $`M(\omega )`$ is the “spin self-energy” and $`I_𝐪`$ describes the exchange (RKKY) interaction between the local moments. Eq. (79) is simply a statement of translational invariance, as it equates the on-site spin susceptibility of the lattice system (LHS) to the susceptibility of any local moment (RHS). The locally critical solution arises when $`\chi _{\mathrm{loc}}(\omega )`$ is singular – signaling the destruction of the Kondo resonance as discussed in previous sections – at the point where the peak susceptibility $`\chi (𝐐_{AF},\omega )`$ just diverges. It is shown in refs.  that one requirement for a locally critical point is that the magnetic fluctuations are two-dimensional. This can be easily seen from Eqs. (79,80). The fact that the peak susceptibility $`\chi (𝐐_{AF},\omega )`$ diverges means that $`M(0)=I_{𝐐_{AF}}`$. Recognizing that for $`𝐪𝐐_{AF}`$, the interaction has the generic form, $`I_𝐪=I_{𝐐_{AF}}+a(𝐪𝐐_{AF})^2`$, we have

$`\chi (𝐪𝐐_{AF},\omega 0)={\displaystyle \frac{1}{a(𝐪𝐐_{AF})^2}}`$ (81)
Now, the fact that the local susceptibility, $`\chi _{\mathrm{loc}}(\omega 0)`$, is singular demands that the $`𝐪`$averaging of the LHS of Eq. (79) should diverge. In other words, integrating the RHS of Eq. (81) over $`𝐪`$ should yield a singular result. Two-dimensionality provides the phase space for this.

This condition, however, is not sufficient. It turns out that a second condition must be met for the existence of such a locally-critical solution. The condition is precisely Eq. (1). The origin of this second condition is somewhat subtle, but has to do with the fact that $`M(\omega )`$ and $`\chi _{\mathrm{loc}}`$, in addition to being related through Eqs. (79,80), must also satisfy a Dyson-like equation of the (self-consistent) Bose-Fermi Kondo model. The details can be found in ref. . By proving that Eq. (1) is valid to all orders in $`ϵ`$, and also in spin-anisotropic cases, we have then established that the locally critical point is indeed a robust result within these microscopic considerations. (The arguments for the robustness of the locally critical point beyond the microscopic approaches can be found in refs. .)

We stress that a key element of the locally critical point of the Kondo lattice is that, the criticality of the local degrees of freedom is embedded in the criticality associated with the long-range ordering in the lattice. In other words, the point where a local energy scale turns to zero coincides with where the spatial correlation length just diverges. This is very different from the self-consistent spin liquid solutions discussed in other contexts, which correspond to a phase instead of a critical point. This difference is also reflected in the fact that stabilizing the locally critical solution in Kondo lattices requires two-dimensional magnetic fluctuations.

## VII Summary

In short, we have carried out a detailed analysis of the Bose-Fermi Kondo model both when the spin-rotational invariance is satisfied (isotropic) and when it is broken (xy and Ising). In each case, we have identified a critical point that separates a Kondo phase, where the local moment is quenched by the spins of the conduction electrons, and a local moment phase where there exists no Kondo resonance. This unstable fixed point marks the point where the spectral weight associated with the Kondo resonance has just gone down to zero.

In all three cases, we find that the exponent for the local spin susceptibility at the unstable fixed point is equal to $`ϵ`$. We note that the three cases correspond to three different one-dimensional statistical mechanical problems. It is remarkable that, the susceptibility exponent at the unstable fixed point is insensitive to these differences. However, the same cannot be said about the stable fixed points: here, the susceptibility exponent for the Ising case is very different from the isotropic and xy cases.

Our results have important consequences for the locally critical behavior in Kondo lattices. In particular, it guarantees that a quantum phase transition of a locally-critical type is a robust microscopic solution to the Kondo lattices: it arises to all orders in $`ϵ`$ and, in addition, not only when the system is spin-rotationally invariant but also in spin-anisotropic situations.

###### Acknowledgements.

This work has been supported in part by TCSUH, the Robert A. Welch Foundation, and NSF Grant No. DMR-0090071. We would also like to acknowledge the support of Argonne National Laboratory, the University of Chicago, and the University of Illinois at Urbana-Champaign during our leave of absence from Rice University. This work was briefly reported in ref. . We have recently learned that G. Zarand and E. Demler have also analyzed this problem, with results which are similar to ours.
## A RG procedure to order $`ϵ^2`$

In this appendix, we give the details of the RG analysis. We adopt the dimensional regularization and minimal subtraction scheme .

We first carry out a renormalized perturbation calculation. The quantities of interest are the $`J`$ and $`g`$ vertices:

$`\mathrm{\Gamma }_J(\omega +\lambda )`$ $`J\mu ^ϵ^{}\gamma _J(\omega +\lambda )`$ (A1)
$`\mathrm{\Gamma }_g(\omega +\lambda )`$ $`g\mu ^{ϵ/2}\gamma _g(\omega +\lambda )`$ (A2)
as well as the $`f`$electron self-energy, $`\mathrm{\Sigma }_f(\omega +\lambda )`$.

$$G^1(\omega +\lambda )=\omega \mathrm{\Sigma }_f(\omega +\lambda )=Z_fG_B^1(\omega +\lambda )$$
(A3)
Within the renormalized perturbation calculation, the bare parameters in Eq. (LABEL:H-BFK-tilde) are replaced with the corresponding renormalized ones. We end up with a Hamiltonian containing one part that has the form of Eq. (3) (the tree level Hamiltonian) and another part which describes the counterterms:

$`_{ct}=`$ $`{\displaystyle \underset{\sigma }{}}(Z_f1)\lambda f_\sigma ^{}f_\sigma +(Z_J1)J\mu ^ϵ^{}𝐒𝐬_c`$ (A5)
$`+(Z_g1)g\mu ^{ϵ/2}𝐒\stackrel{}{\varphi },`$
There are then two types of perturbative diagrams, one coming from the tree level Hamiltonian, which we will call direct perturbative contributions below, and the other from the counterterms.

The singular contributions are kept track of through poles as a function of $`ϵ`$ and $`ϵ^{}`$. By demanding that such poles are “minimally” removed so that the renormalized quantities $`\mathrm{\Gamma }_g`$, $`\mathrm{\Gamma }_J`$, and $`G_f`$ are regular (as a function of $`ϵ`$ and $`ϵ^{}`$), we fix the wavefunction and vertex renormalization factors $`Z_f`$, $`Z_J`$, and $`Z_g`$ order by order in the coupling constants. From these renormalization factors, we can determine the RG beta functions as well as the anomalous dimensions. Note that what we are carrying through is a double expansion in $`ϵ`$ and $`ϵ^{}`$.

In the following, we separate our discussions according to the order in $`J`$ and $`g`$.

### 1 First order result

To the first order of perturbation, there are only corrections to the $`J`$coupling from processes that involve the $`J`$coupling alone. The diagrams are specified by Figs. 5(a)(b) in Ref. . The result is given by

$`\gamma _J^{(1)}(\omega +\lambda )=`$ $`{\displaystyle \frac{N_0J}{ϵ^{}}}({\displaystyle \frac{\mu }{\omega }})^ϵ^{}`$ (A6)
$`=`$ $`{\displaystyle \frac{N_0J}{ϵ^{}}}+N_0J\mathrm{ln}({\displaystyle \frac{\mu }{\omega }})+\mathrm{}`$ (A7)
Adding the counterterm $`(Z_J1)`$ and demanding that the pole terms cancel, we find that

$$(Z_J1)^{(1)}=\frac{1}{ϵ^{}}(N_0J).$$
(A8)
### 2 Second order result

To the order $`g^2`$ and $`J^2`$, the contributions to the $`f`$electron self-energy are shown in Figs. 5a and 5b and they yield the following expression:

$`\mathrm{\Sigma }_f^{(2)}(\omega +\lambda )=`$ $`\omega {\displaystyle \frac{3}{4ϵ}}(K_0g)^2({\displaystyle \frac{\mu }{\omega }})^ϵ`$ (A10)
$`\omega {\displaystyle \frac{3}{16ϵ^{}}}(N_0J)^2({\displaystyle \frac{\mu }{\omega }})^{2ϵ^{}},`$
By demanding that these poles of the self-energy be cancelled by the counterterm contribution $`(Z_f1)\omega `$, we fix

$$(Z_f1)^{(2)}=\frac{3}{4ϵ}(K_0g)^2\frac{3}{16ϵ^{}}(N_0J)^2.$$
(A11)
To the same order, the corrections to the coupling constant $`g`$ are specified by Fig. 6 of Ref.  and they yield

$`\gamma _g^{(2)}(\omega +\lambda )=`$ $`{\displaystyle \frac{1}{4ϵ}}(1ϵ)(K_0g)^2({\displaystyle \frac{\mu }{\omega }})^ϵ`$ (A13)
$`{\displaystyle \frac{1}{16ϵ^{}}}(N_0J)^2({\displaystyle \frac{\mu }{\omega }})^{2ϵ^{}},`$
This fixes the second-order contributionn to $`Z_g`$,

$$(Z_g1)^{(2)}=\frac{1}{4ϵ}(K_0g)^2+\frac{1}{16ϵ^{}}(N_0J)^2.$$
(A14)
The corrections to the coupling constant $`J`$ include not only the direct perturbation diagrams, but also contributions from the first order counterterm $`(Z_J1)^{(1)}J`$. Figs. 5(c)(d) in Ref.  give

$`\gamma _J^{(2)I}(\omega +\lambda )=`$ $`{\displaystyle \frac{1}{4ϵ}}(1ϵ)(K_0g)^2({\displaystyle \frac{\mu }{\omega }})^ϵ`$ (A16)
$`{\displaystyle \frac{1}{16ϵ^{}}}(N_0J)^2({\displaystyle \frac{\mu }{\omega }})^{2ϵ^{}}.`$
The so-called “parquet” diagrams as in Kondo problems, shown in Figs. 6(a)-(f), give

$$\gamma _J^{(2)II}(\omega +\lambda )=\frac{(N_0J)^2}{ϵ^2}(\frac{\mu }{\omega })^{2ϵ^{}}.$$
(A17)
The diagrams with the counterterm $`(Z_J1)^{(1)}J`$, shown in Figs. 6(g)(h), give

$`\gamma _J^{(2)III}(\omega +\lambda )=`$ $`2(Z_J1)^{(1)}\gamma _J^{(1)}(\omega +\lambda )`$ (A18)
$`=`$ $`{\displaystyle \frac{2(N_0J)^2}{ϵ^2}}({\displaystyle \frac{\mu }{\omega }})^ϵ^{},`$ (A19)
where the prefactor 2 is a symmetry factor. From these contributions, we obtain the second order correction to $`Z_J`$:

$$(Z_J1)^{(2)}=\frac{1}{4ϵ}(K_0g)^2+\frac{1}{16ϵ^{}}(N_0J)^2+\frac{1}{ϵ^2}(N_0J)^2.$$
(A20)
### 3 Third order result

To the third order of perturbation, we need only to consider the $`g^2J`$ order corrections to the $`J`$ vertices. (There are also $`J^3`$ order corrections to the $`J`$ vertices, but these contributions are higher than $`ϵ^2`$ order.) The direct perturbative diagrams are shown in Fig. 7. Summing up the contributions of these diagrams, we get

$`\gamma _J^{(3)I}(\omega +\lambda )=`$ $`\left({\displaystyle \frac{1}{4ϵ^{}}}{\displaystyle \frac{1}{4ϵϵ^{}}}{\displaystyle \frac{1}{ϵ(ϵ+ϵ^{})}}\right)`$ (A22)
$`\times (K_0g)^2(N_0J)({\displaystyle \frac{\mu }{\omega }})^{ϵ+ϵ^{}}.`$
Fig. 8 shows the diagrams with counterterms to $`g^2J`$ order. The contributions from Figs. 8(a)(b) give rise to

$`\gamma _J^{(3)II}(\omega +\lambda )=`$ $`(Z_f1)^{(2)}\gamma _J^{(1)}(\omega +\lambda )`$ (A23)
$`=`$ $`{\displaystyle \frac{3(K_0g)^2(N_0J)}{4ϵϵ^{}}}({\displaystyle \frac{\mu }{\omega }})^ϵ^{};`$ (A24)
those from Figs. 8(c)(d) yield

$`\gamma _J^{(3)III}(\omega +\lambda )=`$ $`2(Z_J1)^{(2)}\gamma _J^{(1)}(\omega +\lambda )`$ (A25)
$`=`$ $`{\displaystyle \frac{(K_0g)^2(N_0J)}{2ϵϵ^{}}}({\displaystyle \frac{\mu }{\omega }})^ϵ^{};`$ (A26)
and finally, Fig. 8(e) leads to

$`\gamma _J^{(3)IV}(\omega +\lambda )=`$ $`(Z_J1)^{(1)}\gamma _J^{(2)}(\omega +\lambda )`$ (A27)
$`=`$ $`\left({\displaystyle \frac{1}{4ϵϵ^{}}}{\displaystyle \frac{1}{4ϵ^{}}}\right)(K_0g)^2(N_0J)({\displaystyle \frac{\mu }{\omega }})^ϵ.`$ (A28)
Summing up $`\gamma _J^{(3)I,II,III,IV}`$ and subtracting poles, we have

$$(Z_J1)^{(3)}=\left(\frac{1}{ϵ^{}(ϵ+ϵ^{})}+\frac{1}{4ϵϵ^{}}\right)(K_0g)^2(N_0J).$$
(A29)
Note that the simple poles in Eqs. (A22,A28) cancel with each other, leaving only double poles in the net result, Eq. (A29). Therefore, no $`g^2J`$ terms will appear in the beta functions.

### 4 Fourth order result

We now turn to contributions of order $`g^4`$, the only terms we are interested in to this order. First, we consider the $`f`$electron self-energy. Figs. 5(c), (d) are the direct perturbative contributions to $`f`$electron self-energy of this order. Figs. 5(e), (f), on the other hand, are the diagrams with counter terms of order $`g^2`$ and $`J^2`$. Keeping only the $`g^4`$ order terms, we have:

$`\mathrm{\Sigma }_f^{(4)(\text{5}c+\text{5}d)}(\omega +\lambda )=`$ $`\omega \left({\displaystyle \frac{15}{32ϵ^2}}+{\displaystyle \frac{3}{8ϵ}}\right)(K_0g)^4({\displaystyle \frac{\mu }{\omega }})^{2ϵ},`$ (A30)
$`\mathrm{\Sigma }_f^{(4)(\text{5}e)}(\omega +\lambda )=`$ $`(Z_f1)^{(2)}\mathrm{\Sigma }_f^{(2)}(\omega +\lambda )`$ (A31)
$`=`$ $`\omega {\displaystyle \frac{9}{16ϵ^2}}(K_0g)^4({\displaystyle \frac{\mu }{\omega }})^ϵ`$ (A32)
$`\mathrm{\Sigma }_f^{(4)(\text{5}f)}(\omega +\lambda )=`$ $`2(Z_g1)^{(2)}\mathrm{\Sigma }_f^{(2)}(\omega +\lambda )`$ (A33)
$`=`$ $`\omega {\displaystyle \frac{3}{8ϵ^2}}(K_0g)^4({\displaystyle \frac{\mu }{\omega }})^ϵ`$ (A34)
The resulting corrections to $`Z_f`$ in this order is given by

$$(Z_f1)^{(4)}=\frac{15}{32ϵ^2}(K_0g)^4+\frac{3}{8ϵ}(K_0g)^4.$$
(A35)
Next, we consider the $`g`$vertex corrections. They are specified in Fig. 9. The contributions from the direct perturbative diagrams, Figs. 9(a)-(f), sum up to

$$\gamma _g^{(4)af}(\omega +\lambda )=\left(\frac{9}{32ϵ^2}\frac{7}{16ϵ}\right)(K_0g)^4(\frac{\mu }{\omega })^{2ϵ};$$
(A36)
while those from the counterterms, Figs. 9(g)(h), yield,

$`\gamma _g^{(4)g,h}(\omega +\lambda )=`$ $`3(Z_g1)^{(2)}\gamma _g^{(2)}(\omega +\lambda )`$ (A38)
$`2(Z_f1)^{(2)}\gamma _g^{(2)}(\omega +\lambda )`$
$`=`$ $`\left({\displaystyle \frac{9}{16ϵ^2}}+{\displaystyle \frac{9}{16ϵ}}\right)(K_0g)^4({\displaystyle \frac{\mu }{\omega }})^ϵ.`$ (A39)
Therefore, the correction to $`Z_g`$ is given by

$$(Z_g1)^{(4)}=\frac{9}{32ϵ^2}(K_0g)^4\frac{1}{8ϵ}(K_0g)^4.$$
(A40)
Similarly, the $`J`$coupling corrections contain both direct perturbative diagrams (Fig. 10) and diagrams with counterterms (Fig. 11) which are given by

$$\gamma _J^{(4)I}(\omega +\lambda )=\left(\frac{9}{32ϵ^2}\frac{7}{16ϵ}\right)(K_0g)^4(\frac{\mu }{\omega })^{2ϵ},$$
(A41)
and

$`\gamma _J^{(4)II}(\omega +\lambda )=`$ $`\left(2(Z_g1)^{(2)}2(Z_f1)^{(2)}\right)\gamma _J^{(2)}(\omega +\lambda )`$ (A43)
$`+(Z_J1)^{(2)}\gamma _J^{(2)}(\omega +\lambda )`$
$`=`$ $`\left({\displaystyle \frac{9}{16ϵ^2}}+{\displaystyle \frac{9}{16ϵ}}\right)(K_0g)^4({\displaystyle \frac{\mu }{\omega }})^ϵ,`$ (A44)
respectively. The resulting correction to $`Z_J`$ is

$$(Z_J1)^{(4)}=\frac{9}{32ϵ^2}(K_0g)^4\frac{1}{8ϵ}(K_0g)^4.$$
(A45)
We can now collect the renormalization factors $`Z_f`$, $`Z_g`$, and $`Z_J`$ to the orders of our interest. The results are given in Eq. (24).

### 5 The beta functions

After obtaining the renormalization factors, we can calculate the beta functions, defined in Eq. (26).

Taking the $`\mu `$ derivative of Eqs. (14,16) at the fixed bare couplings, we have

$`0=`$ $`{\displaystyle \frac{\beta (g)}{g}}+\beta (g){\displaystyle \frac{\mathrm{ln}(Z_f^1Z_g)}{g}}`$ (A47)
$`+\beta (J){\displaystyle \frac{\mathrm{ln}(Z_f^1Z_g)}{J}}+{\displaystyle \frac{ϵ}{2}},`$
$`0=`$ $`{\displaystyle \frac{\beta (J)}{J}}+\beta (g){\displaystyle \frac{\mathrm{ln}(Z_f^1Z_J)}{g}}`$ (A49)
$`+\beta (J){\displaystyle \frac{\mathrm{ln}(Z_f^1Z_J)}{J}}+ϵ^{}.`$
After some algebra, and setting $`ϵ^{}=0`$, we obtain the beta functions of $`g`$ and $`J`$ as in Eq. (29).

### 6 Renormalizability

We close this appendix by commenting on two issues regarding renormalizability. Our perturbative results for the self-energy and vertex corrections as well as the renormalization factors contain double poles such as $`1/ϵ^2`$, $`1/(ϵ^{})^2`$ and $`1/ϵϵ^{}`$. However, we find that terms of the form $`\frac{1}{ϵ}\mathrm{ln}\frac{\mu }{\omega }`$ all cancel out after we carry out expansions such as

$`({\displaystyle \frac{\mu }{\omega }})^ϵ=1+ϵ\mathrm{ln}{\displaystyle \frac{\mu }{\omega }}+O(ϵ^2).`$ (A50)
This is in accordance with the requirement of renormalizability . In addition, we find that the double poles in the renormalization factors cancel out in the beta functions. This is also a consistency check for the minimal subtraction procedure, and is in accordance with the requirement that beta functions are analytic.

## B Calculations of the local spin susceptibility

### 1 The local susceptibility to order $`ϵ^2`$

The local spin correlation function is defined in Eq. (35). We find it convenient to carry out our calculation in $`\tau `$ space. (The self-energy and vertex correction diagrams in Appendix A were calculated in $`\omega `$ space.) So we take the expressions for the boson and conduction electron propagators as

$`G_\varphi ^0(\tau )=`$ $`\stackrel{~}{K}_{0}^{}{}_{}{}^{2}\left({\displaystyle \frac{\pi /\beta }{\mathrm{sin}\frac{\pi \tau }{\beta }}}\right)^{2ϵ}`$ (B1)
$`G_c^0(\tau )=`$ $`\stackrel{~}{N}_0\left({\displaystyle \frac{\pi /\beta }{\mathrm{sin}\frac{\pi \tau }{\beta }}}\right)^{1ϵ^{}}`$ (B2)
where $`\stackrel{~}{K}_0`$ and $`\stackrel{~}{N}_0`$ can be determined from $`K_0`$ and $`N_0`$ by

$`\stackrel{~}{K}_{0}^{}{}_{}{}^{2}=K_0^2\mathrm{\Gamma }(2ϵ)`$ (B3)
$`\stackrel{~}{N}_{0}^{}{}_{}{}^{2}=N_0^2{\displaystyle \frac{\mathrm{\Gamma }(22ϵ^{})}{\mathrm{\Gamma }^2(1ϵ^{})}},`$ (B4)
where $`\mathrm{\Gamma }(x)`$ is the Euler Gamma function.

We now calculate the bare spin susceptibility in terms of the bare Hamiltonian (instead of the renormalized Hamiltonian used in Appendix A). The zero-th order result is given by

$`\stackrel{~}{\chi }_B^{(0)}(\tau )`$ $`={\displaystyle \frac{1}{2}}{\displaystyle _0^\beta }𝑑\tau G_f^0(\tau )G_f^0(\tau )={\displaystyle \frac{1}{2}}e^{\beta \lambda }`$ (B5)
$`\chi _B^{(0)}(\tau )`$ $`=\underset{\lambda \mathrm{}}{lim}{\displaystyle \frac{1}{2}}e^{\beta \lambda }\stackrel{~}{\chi }_B(\tau )={\displaystyle \frac{1}{4}},`$ (B6)
where $`G_f^0(\tau )`$ is the bare $`f`$electron propagator. This result is simply $`S(S+1)/3`$ with $`S=1/2`$.

To the order $`g^2`$ and $`J^2`$, the diagrams have been shown in Fig. 8 of Ref. . Using the dimensional regularization adopted here, we find

$`\chi _B^{(2)}(\tau )=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \frac{2(K_0g_B)^2\tau ^ϵ}{ϵ}}\left(1+Cϵ+{\displaystyle \frac{1}{2}}(C^2+{\displaystyle \frac{\pi ^2}{6}})ϵ^2\right)`$ (B8)
$`{\displaystyle \frac{1}{4}}{\displaystyle \frac{(N_0J_B)^2\tau ^{2ϵ^{}}}{2ϵ^{}}}(1+{\displaystyle \frac{\pi ^2}{6}}ϵ^2),`$
where $`C=0.577215\mathrm{}`$ is the Euler Gamma constant.

To $`g^4`$ order, the diagrams are shown in Fig. 12. (Structurally, they are similar to those considered in ref. .) The result is

$`\chi _B^{(4)}(\tau )=`$ $`{\displaystyle \frac{1}{4}}({\displaystyle \frac{4}{ϵ^2}}+{\displaystyle \frac{1+8C}{ϵ}}+2+2C+8C^2+{\displaystyle \frac{\pi ^2}{3}})\times `$ (B10)
$`\times (K_0g_B)^4\tau ^{2ϵ}`$
We now express $`\chi _B`$ in terms of the renormalized coupling constants, $`g`$ and $`J`$, by using Eqs. (14,16) in combination with our results for the renormalization factors, Eq. (24). We can then determine the renormalization factor $`Z`$, defined in Eq. (36), by demanding that all the pole terms in $`\chi _B`$ are removed by those of $`Z^1`$ so that the renormalized susceptibility, $`\chi =Z^1\chi _B`$, is regular. The result is

$`Z^1=`$ $`1+{\displaystyle \frac{2(K_0g)^2}{ϵ}}+{\displaystyle \frac{(N_0J)^2}{2ϵ^{}}}.`$ (B12)
$`+{\displaystyle \frac{4(K_0g)^4}{ϵ^2}}{\displaystyle \frac{(K_0g)^4}{ϵ}}.`$
Eq. (37) then follows.

The renormalized susceptibility, $`\chi =Z^1\chi _B`$, is equal to

$$\chi (\tau )=\chi ^{(0)}\left[A(g,J,ϵ)B(g,J,ϵ)\mathrm{ln}(\mu \tau )+O(\mathrm{ln}^2(\mu \tau ))\right],$$
(B13)
where, to the $`ϵ^2`$ order,

$`A(g,J,ϵ)=`$ $`12(K_0g)^2(C+{\displaystyle \frac{1}{2}}(C^2+{\displaystyle \frac{\pi ^2}{6}})ϵ)`$ (B15)
$`+(K_0g)^4(2+2C+4C^2{\displaystyle \frac{\pi ^2}{3}})`$
$`B(g,J,ϵ)=`$ $`2(K_0g)^2(1+Cϵ)+(N_0J)^22(K_0g)^4(1+4C).`$ (B16)
At the critical point(“C” in Fig. 1), the spin correlation has the form

$$\chi (\tau )=\chi ^{(0)}A(g^{},J^{},ϵ)\frac{1}{(\mu \tau )^\eta },$$
(B17)
where the exponent $`\eta `$ can either be calculated from the renormalization factor $`Z`$ in Eq. (42), or by $`B(g^{},J^{},ϵ)/A(g^{},J^{},ϵ)`$, which also yields $`\eta =ϵ`$. The renormalized amplitude factor $`A(g^{},J^{},ϵ)`$ is given by

$$A(g^{},J^{},ϵ)=1Cϵ+\left(\frac{1}{2}+\frac{1}{4}C+\frac{1}{2}C^2\frac{1}{6}\pi ^2\right)ϵ^2.$$
(B18)
For large values of $`ϵ`$, it would require an appropriate resummation such as Padé approximation.

For the local moment fixed point, the local spin correlation has the form

$`\chi (\tau )=`$ $`\chi ^{(0)}\left[A_L(ϵ)B_L(ϵ)\mathrm{ln}(\mu \tau )+O(\mathrm{ln}^2(\mu \tau ))\right],`$ (B19)
$`A_L(ϵ)=`$ $`A(g,J,ϵ)|_{g=g_L^{},J=J_L^{}=0}`$ (B20)
$`=`$ $`1Cϵ+\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}C^2{\displaystyle \frac{1}{6}}\pi ^2\right)ϵ^2,`$ (B21)
$`B_L(ϵ)=`$ $`B(g,J,ϵ)|_{g=g_L^{},J=J_L^{}=0}`$ (B22)
$`=`$ $`ϵCϵ^2.`$ (B23)
Again, $`B_L(ϵ)/A_L(ϵ)=ϵ`$ to the $`ϵ^2`$ order.

### 2 The local susceptibility to all orders in $`ϵ`$

We now discuss the contributions to $`\chi `$, to all orders of the perturbation theory. We can group all the diagrams in a manner illustrated in Fig. 13. Here, each double line represents the full $`f`$electron propagator, while a shaded area denotes the full vertex. The renormalization factor for the spin, $`Z`$, is the product of the renormalization factors for the full f-electron propagators and those for the vertices. The renormalization factor for a full f-electron propagator is equal to $`Z_f`$. By inspecting the diagrams for each vertex, it is straightforward to see that to each order, this vertex is identical to that of the $`g`$vertex. As a result, each vertex contributes a factor $`Z_g^1`$. Taking these two things together, we end up with Eq. (45), to all orders. As discussed in Section IV B, Eq. (45) then leads to Eqs. (47,48).

## C RG analysis for anisotropic cases

In the anisotropic case, the diagrams are topologically the same as those in the isotropic case, except for different coupling vertices. In the following, we will simply list the results.

The renormalization factors are given as

$`Z_f`$ $`=`$ $`1{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^2}{2}}+{\displaystyle \frac{(K_0g_z)^2}{4}}\right)`$ (C1)
$``$ $`{\displaystyle \frac{1}{ϵ^2}}\left({\displaystyle \frac{(K_0g_{})^4}{8}}+{\displaystyle \frac{3(K_0g_{})^2(K_0g_z)^2}{8}}{\displaystyle \frac{(K_0g_z)^4}{32}}\right)`$ (C2)
$`+`$ $`{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^4}{8}}+{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{4}}\right)`$ (C3)
$``$ $`{\displaystyle \frac{1}{ϵ^{}}}\left({\displaystyle \frac{(N_0J_{})^2}{8}}+{\displaystyle \frac{(N_0J_z)^2}{16}}\right),`$ (C4)
$`Z_g_{}`$ $`=`$ $`1+{\displaystyle \frac{(K_0g_z)^2}{4ϵ}}`$ (C5)
$`+`$ $`{\displaystyle \frac{1}{ϵ^2}}\left({\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{4}}+{\displaystyle \frac{(K_0g_z)^4}{32}}\right)`$ (C6)
$``$ $`{\displaystyle \frac{(K_0g_{})^4}{8ϵ}}+{\displaystyle \frac{(N_0J_z)^2}{16ϵ^{}}},`$ (C7)
$`Z_{g_z}`$ $`=`$ $`1+{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^2}{2}}{\displaystyle \frac{(K_0g_z)^2}{4}}\right)`$ (C11)
$`+{\displaystyle \frac{1}{ϵ^2}}\left({\displaystyle \frac{3(K_0g_{})^4}{8}}{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{8}}+{\displaystyle \frac{(K_0g_z)^4}{32}}\right)`$
$`+{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^4}{8}}{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{4}}\right)`$
$`+{\displaystyle \frac{1}{ϵ^{}}}\left({\displaystyle \frac{(N_0J_{})^2}{8}}{\displaystyle \frac{(N_0J_z)^2}{16}}\right),`$
$`Z_J_{}`$ $`=`$ $`1{\displaystyle \frac{N_0J_z}{ϵ^{}}}`$ (C16)
$`+{\displaystyle \frac{(K_0g_z)^2}{4ϵ}}+{\displaystyle \frac{(N_0J_z)^2}{16ϵ^{}}}+{\displaystyle \frac{(N_0J_{})^2+(N_0J_z)^2}{2ϵ^2}}`$
$`{\displaystyle \frac{(K_0g_{})^2(N_0J_z)}{ϵ^{}(ϵ+ϵ^{})}}{\displaystyle \frac{(K_0g_z)^2(N_0J_z)}{4ϵϵ^{}}}`$
$`+{\displaystyle \frac{1}{ϵ^2}}\left({\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{4}}+{\displaystyle \frac{(K_0g_z)^4}{32}}\right)`$
$`{\displaystyle \frac{(K_0g_{})^4}{8ϵ}}`$
$`Z_{J_z}=`$ $`1{\displaystyle \frac{N_0(J_{}^2/J_z)}{ϵ^{}}}+{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^2}{2}}{\displaystyle \frac{(K_0g_z)^2}{4}}\right)`$ (C22)
$`+{\displaystyle \frac{1}{ϵ^{}}}\left({\displaystyle \frac{(N_0J_{})^2}{8}}{\displaystyle \frac{(N_0J_z)^2}{16}}\right)+{\displaystyle \frac{(N_0J_{})^2}{ϵ^2}}`$
$`{\displaystyle \frac{(K_0g_z)^2(N_0J_z)}{ϵ^{}(ϵ+ϵ^{})}}`$
$`{\displaystyle \frac{\left(2(K_0g_{})^2(K_0g_z)^2\right)(N_0J_{}^2/J_z)}{4ϵϵ^{}}}`$
$`+{\displaystyle \frac{1}{ϵ^2}}\left({\displaystyle \frac{3(K_0g_{})^4}{8}}{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{8}}+{\displaystyle \frac{(K_0g_z)^4}{32}}\right)`$
$`+{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^4}{8}}{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{4}}\right)`$
Eqs. (C4)- (C22) reduce to the isotropic results, Eq. (24), when we set $`g_z=g_{}`$, $`J_z=J_{}`$. Similar to the isotropic case, we obtain the beta functions for each coupling constant from these renormalization factors, which are given in Eqs. (58)- (67). We can then set $`g_z=0`$ ($`g_{}=0`$) to discuss the xy (Ising) case.

We now turn to the calculation of the local spin susceptibility, focusing on the xy case only. Like in the isotropic case, we use the bare Hamiltonian. The contributions from the perturbative diagrams to the $`g^2`$, $`J^2`$, and $`g^4`$ orders are

$`\chi _B^{+(2)}(\tau )=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \frac{(K_0g_{}^B)^2+(K_0g_z^B)^2}{ϵ}}\tau ^ϵ`$ (C25)
$`\times \left(1+Cϵ+{\displaystyle \frac{1}{2}}(C^2+{\displaystyle \frac{\pi ^2}{6}})ϵ^2\right)`$
$`{\displaystyle \frac{(N_0J_{}^B)^2+(N_0J_z^B)^2}{16ϵ^{}}}\tau ^{2ϵ^{}}(1+{\displaystyle \frac{\pi ^2}{6}}ϵ^2),`$
$`\chi _B^{+(4)}(\tau )=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \frac{(K_0g_{}^B)^4+(K_0g_{}^B)^2(K_0g_z^B)^2}{2}}\tau ^{2ϵ}\times `$ (C27)
$`\times \left({\displaystyle \frac{4}{ϵ^2}}+{\displaystyle \frac{1+8C}{ϵ}}+2+2C+8C^2+{\displaystyle \frac{\pi ^2}{3}}\right).`$
(There are no contributions at the $`J`$ or $`g^2J`$ order.) Following the same RG procedure as in Appendix B, we can calculate the renormalization factor $`Z`$:

$`Z`$ $`=`$ $`1{\displaystyle \frac{1}{ϵ}}\left((K_0g_{})^2+(K_0g_z)^2\right)`$ (C28)
$`+`$ $`{\displaystyle \frac{1}{ϵ}}\left({\displaystyle \frac{(K_0g_{})^4}{2}}+{\displaystyle \frac{(K_0g_{})^2(K_0g_z)^2}{2}}\right)`$ (C29)
$``$ $`{\displaystyle \frac{1}{ϵ^{}}}\left({\displaystyle \frac{(N_0J_{})^2}{4}}+{\displaystyle \frac{(N_0J_z)^2}{4}}\right),`$ (C30)
and the anomalous dimension $`\eta `$ at the critical point:

$`\eta `$ $`=`$ $`{\displaystyle \frac{d\mathrm{log}Z}{d\mathrm{log}\mu }}|_{g_i=g_i^{},J_i=J_i^{}}`$ (C31)
$`=`$ $`(K_0g_{}^{})^2+(K_0g_z^{})^2(K_0g_{}^{})^4`$ (C32)
$``$ $`(K_0g_{}^{})^2(K_0g_z^{})^2+{\displaystyle \frac{(N_0J_{}^{})^2}{2}}+{\displaystyle \frac{(N_0J_z^{})^2}{2}}`$ (C33)
Setting $`g_z=0`$ (for the xy case), the exponent becomes

$$\eta =(K_0g_{}^{})^2(K_0g_{}^{})^4+\frac{(N_0J_{}^{})^2+(N_0J_z^{})^2}{2}.$$
(C34)
The renormalized amplitude $`A(g,J,ϵ)`$ is given by

$`A(g,J,ϵ)=`$ $`1(K_0g_{})^2(C+{\displaystyle \frac{1}{2}}(C^2+{\displaystyle \frac{\pi ^2}{6}})ϵ)`$ (C36)
$`+(K_0g_{})^4(1+C+2C^2{\displaystyle \frac{\pi ^2}{6}})`$
At the critical point, taking account of Eq. (72), we find

$`\eta `$ $`=ϵ,`$ (C37)
$`A_C(ϵ)`$ $`=1Cϵ+\left(1+{\displaystyle \frac{3}{8}}C+{\displaystyle \frac{3}{2}}C^2{\displaystyle \frac{1}{4}}\pi ^2\right)ϵ^2.`$ (C38)
Near the bosonic stable fixed point where $`J_{}^{}=J_z^{}=0`$, which is given by Eq. (74), the results are

$`\eta _L=`$ $`(K_0g_{}^{})^2(K_0g_{}^{})^4`$ (C39)
$`=`$ $`ϵ,`$ (C40)
$`A_L(ϵ)=`$ $`1Cϵ+\left(1+{\displaystyle \frac{3}{2}}C^2{\displaystyle \frac{1}{4}}\pi ^2\right)ϵ^2.`$ (C41)
## D Kink-gas analysis of the Ising case

We now turn to the Ising case, by setting $`g_{}=0`$. The mapping to the kink-gas action is similar to that given in ref. . The resulting action, describing a one-dimensional long-ranged statistical-mechanical model, is similar to Eq. (7) of ref.  and, in our notation, takes the following form:

$`S(\tau _{2n},`$ $`\mathrm{},\tau _1)=2n\mathrm{ln}(y_j)+{\displaystyle }_l(1)^lh(\tau _{l+1}\tau _l)/\tau _0`$ (D2)
$`+{\displaystyle \underset{l<m}{}}(1)^{l+m}[2\kappa _j\mathrm{ln}(\tau _m\tau _l)/\tau _0+K(\tau _m\tau _l)]`$
where $`[\tau _{2n},\mathrm{},\tau _1]`$, for $`n=1,2,\mathrm{}`$, labels a sequence of spin flips (kinks) along the imaginary time axis, and $`y_j`$ and $`\kappa _j`$ are defined in Eq. (77). The last term, $`K(\tau )`$, originates from the $`g_z`$coupling and takes the form:

$`K(\tau )=\kappa _g[(\tau /\tau _0)^ϵ1]/ϵ`$ (D3)
where $`\kappa _g`$ is also specified by Eq. (77). The RG equations for this kink-gas problem has already been derived in ref.  and are given as follows:

$`\beta (y_j)`$ $`=y_j(1\kappa _j\kappa _g/2)`$ (D4)
$`\beta (\kappa _j)`$ $`=4\kappa _jy^2`$ (D5)
$`\beta (\kappa _g)`$ $`=\kappa _g(ϵ4y_j^2)`$ (D6)
$`\beta (h)`$ $`=h(12y_j^2)`$ (D7)
This procedure is valid for arbitrary values of the stiffness constants $`\kappa _j`$ and $`\kappa _g`$, provided the fugacity $`y_j`$ is small.

An unstable fixed point still exists. To the linear order in $`ϵ`$, it is located at

$`(y_j^{})^2`$ $`=ϵ/4`$ (D8)
$`\kappa _g^{}`$ $`=2`$ (D9)
$`\kappa _j^{}`$ $`=0`$ (D10)
From Eq. (D7), we can easily see that $`\kappa _g`$ is in fact irrelevant around this fixed point. We can then determine the critical exponents by staying within the $`y_j\kappa _g`$ plane: the projection of the RG flow on the $`y_j\kappa _g`$ plane is given in Fig. 4. Within the $`ϵ`$expansion, then, the critical exponents of the unstable fixed point become identical to their counterparts for a similar fixed point of the classical ferromagnetic Ising chain with an long-range interaction that decays in distance in terms of a power-law exponent $`2ϵ`$. The critical exponent $`\eta `$, for instance, can be straightforwardly calculated by combing Eq. (D7) and Eq. (D10). The result is $`\eta =ϵ`$, the same value as it takes at the critical point of the corresponding long-ranged Ising chain.