text:
# The WIMP annual modulation signal and non-standard halo models

## I Introduction

The rotation curves of spiral galaxies are typically flat out to about $`30`$ kpc. This implies that the mass enclosed increases linearly with radius, with a halo of dark matter extending beyond the luminous matter . The nature of the dark matter is unknown , with possible candidates including MAssive Compact Halo Objects (MACHOs), such as brown dwarves, Jupiters or black holes and elementary particles, known as Weakly Interacting Massive Particles (WIMPs), such as axions and the neutralino, the lightest supersymmetric particle.

WIMPs can be directly detected via their elastic scattering off target nuclei. In the long term the directional dependence of the detector recoil will provide the best means of direct WIMP detection . Currently the best prospect for distinguishing a WIMP signal from the detector background is via the annual modulation of the signal, which occurs due to the Earth’s rotation around the sun . The DAMA collaboration , using a detector consisting of radiopure NaI crystal scintillators at the Gran Sasso Laboratory, have recently reported a $`4\sigma `$ annual modulation signal consistent with the detection of a WIMP with mass $`m_\chi =52_8^{+10}`$ GeV , assuming a Maxwellian halo with velocity dispersion roughly equal to the local rotation velocity, $`v_0=220\mathrm{k}\mathrm{m}\mathrm{s}^1`$ .

The values of the WIMP mass and interaction cross section found from the annual modulation signal are known to depend somewhat on the assumed values of poorly known astrophysical parameters. The effects of bulk rotation  and varying the velocity dispersion , within the observationally allowed range $`v_0=220\pm 40\mathrm{k}\mathrm{m}\mathrm{s}^1`$ , have been examined for a Maxwellian halo. The DAMA collaboration have subsequently included these uncertainties in the analysis of their latest data .

The standard Maxwellian halo model has a number of deficiencies (for example see Ref.  and references therein). The halo may not be spherical; N body simulations of gravitational collapse produce axisymmetric or triaxial halos, and indeed other spiral galaxies appear to have flattened halos . The power-law halo models of Evans  provide an analytically tractable framework for investigating the effect of varying halo properties such as the flattening, and have previously been used to investigate variations in the mean total  and directional  WIMP detection rates. Another possible modification to the standard halo model is an asymmetric velocity distribution  and the annual modulation of the total  and directional  WIMP signals have been calculated for various asymmetric velocity distributions.

In this paper we compare the variation of the annual modulation signal for power-law halo models, asymmetric velocity distributions and Maxwellian halos with bulk rotation. We then examine how the range of WIMP masses consistent with the latest DAMA data depends on the halo model assumed.

## II Annual Modulation signal

The WIMP detection rate depends on the speed distribution of the WIMPs in the rest frame of the detector, $`f_v`$. This is found from the halo velocity distribution, $`f(𝐯)`$ by making a Galilean transformation $`𝐯\stackrel{~}{𝐯}=𝐯+𝐯_\mathrm{e}`$, where $`𝐯_\mathrm{e}`$ is the Earth’s velocity relative to the galactic rest frame, and then integrating over the angular distribution ***The effect of the Sun’s gravity on the WIMP distribution can be neglected . In galactic co-ordinates the axis of the ecliptic lies very close to the $`\varphi z`$ plane and is inclined at an angle $`\gamma 29.80^{}`$ to the $`\varphi r`$ plane . Including all components of the Earth’s motion, not just that parallel to the galactic rotation:

$`𝐯_\mathrm{e}`$ $`=`$ $`v_1\mathrm{sin}\alpha \widehat{r}+`$ (2)
$`(v_0+v_1\mathrm{cos}\alpha \mathrm{sin}\gamma )\widehat{\varphi }v_1\mathrm{cos}\alpha \mathrm{cos}\gamma \widehat{z},`$
where $`v_0232\mathrm{k}\mathrm{m}\mathrm{s}^1`$ is the speed of the sun with respect to the galactic rest frame, $`v_130\mathrm{k}\mathrm{m}\mathrm{s}^1`$ is the orbital speed of the Earth around the Sun and $`\alpha =2\pi (tt_0)/T`$, with $`T=1`$ year and $`t_0153`$ days (June 2nd), when the component of the Earth’s velocity parallel to the Sun’s motion is largest.

In the range of masses and interaction cross sections accessible to current direct detection experiments the best motivated WIMP candidate is the neutralino, for which the event rate is dominated by the scalar contribution. The differential event rate can then be written in terms of the WIMP cross section on the proton,

$$\sigma _\mathrm{p}=\frac{4m_\mathrm{p}^2m_\chi ^2}{\pi (m_\mathrm{p}+m_\chi )^2}f_\mathrm{p}^2,$$
(3)
where $`m_\mathrm{p}`$ is the proton mass and $`f_\mathrm{p}`$ the effective WIMP cross section on the proton. The differential event rate simplifies to (see the Appendix for more details):

$$\frac{\mathrm{d}R}{\mathrm{d}E}=\xi \sigma _\mathrm{p}\left[\frac{\rho _{0.3}}{\sqrt{\pi }v_0}\frac{(m_\mathrm{p}+m_\chi )^2}{m_\mathrm{p}^2m_\chi ^3}A^2T(E)F^2(q)\right],$$
(4)
where the local WIMP density, $`\rho _\chi `$ has been normalised to a fiducial value $`\rho _{0.3}=0.3\mathrm{GeV}/\mathrm{cm}^3`$, such that $`\xi =\rho _\chi /\rho _{0.3}`$, $`E`$ is the recoil energy of the detector nucleus and $`T(E)`$ is defined as

$$T(E)=\frac{\sqrt{\pi }v_0}{2}_{v_{\mathrm{min}}}^{\mathrm{}}\frac{f_v}{v}dv,$$
(5)
where $`v_{\mathrm{min}}`$ is the minimum detectable WIMP velocity

$$v_{\mathrm{min}}=\left(\frac{E(m_\chi +m_A)^2}{2m_\chi ^2m_A}\right)^{1/2},$$
(6)
and $`m_A`$ the atomic mass of the target nuclei.

In order to compare the theoretical signal with that observed we need to take into account the response of the detector. The electron equivalent energy, $`E_{\mathrm{ee}}`$, which is actually measured is a fixed fraction of the recoil energy: $`E_{\mathrm{ee}}=q_\mathrm{A}E`$. The quenching factors for I and Na are $`q_\mathrm{I}=0.09`$ and $`q_{\mathrm{Na}}=0.30`$ respectively . The energy resolution of the detector  is already taken into account in the data released by the DAMA collaboration.

The expected experimental spectrum per energy bin for the DAMA collaboration set-up is then given by

$`{\displaystyle \frac{\mathrm{\Delta }R}{\mathrm{\Delta }E}}(E)`$ $`=`$ $`r_{\mathrm{Na}}{\displaystyle _{E/q_{\mathrm{na}}}^{(E+\mathrm{\Delta }E)/q_{\mathrm{na}}}}{\displaystyle \frac{\mathrm{d}R_{\mathrm{Na}}}{\mathrm{d}E_{\mathrm{ee}}}}(E_{\mathrm{ee}}){\displaystyle \frac{\mathrm{d}E_{\mathrm{ee}}}{\mathrm{\Delta }E}}`$ (8)
$`+r_\mathrm{I}{\displaystyle _{E/q_\mathrm{I}}^{(E+\mathrm{\Delta }E)/q_\mathrm{I}}}{\displaystyle \frac{\mathrm{d}R_\mathrm{I}}{\mathrm{d}E_{\mathrm{ee}}}}(E_{\mathrm{ee}}){\displaystyle \frac{\mathrm{d}E_{\mathrm{ee}}}{\mathrm{\Delta }E}},`$
where $`r_{\mathrm{Na}}=0.153`$ and $`r_\mathrm{I}=0.847`$ are the mass fractions of Na and I respectively.

Since $`v_0v_1`$ the differential event rate in the k-th energy bin can be expanded in a Taylor series in $`\mathrm{cos}\alpha `$ :

$$\frac{\mathrm{\Delta }R}{\mathrm{\Delta }E}(E_\mathrm{k})S_{0,\mathrm{k}}+S_{\mathrm{m},\mathrm{k}}\mathrm{cos}\alpha .$$
(9)
The DAMA collaboration use a maximum likelihood method to separate the time-independent background from the WIMP signal and have released the resulting values of $`S_{0,\mathrm{k}}`$ and $`S_{\mathrm{m},\mathrm{k}}`$, and the errors on them, for each energy bin  (see Table 1).

## III Halo models

Since the resolution of numerical simulations is not yet large enough to allow the numerical calculation of the annual modulation signal , we have to use analytic models for the velocity distribution of the dark matter halo. If the galactic halo contains large amounts of substructure, as N-body simulations appear to indicate , then the WIMP detection rate may be dramatically altered if we are currently passing through a clump of substructure. The extent to which substructure is actually present in the halo is currently a matter of debate however.

The effect of the halo model on the annual modulation signal can be assessed most simply, with least recourse to the detector properties, in terms of the dimensionless function $`T(E)`$ The shape of the energy spectrum measured also depends on the nuclear form factor $`F(q)`$ however., as defined in Eq. (5). The yearly averaged value of $`T(E)`$ depends weakly on the halo model  decreasing with increasing recoil energy, with the decrease being most rapid for small WIMP masses. For each model we plot the annual variation in $`T(E)`$:

$$\mathrm{\Delta }T(E)=\frac{(T(E)_{\mathrm{max}}T(E)_{\mathrm{av}})}{T(E)_{\mathrm{av}}},$$
(10)
defined such that $`\mathrm{\Delta }T(E)`$ is taken to be positive if $`T(E)`$ is largest in June ($`\alpha =0`$) as found by the DAMA collaboration, as a function of recoil energy for a $`\mathrm{Ge}^{76}`$ detector and four values of the WIMP mass: $`m_\chi =30,50,100,200`$ GeV. Values for other monatomic detectors can be found by rescaling the x-axis by $`m_\mathrm{A}/(m_\mathrm{A}+m_\chi )^2`$.

### A Asymmetric velocity distribution

An extension of the Michie model can be used as a simple model of velocity anisotropy

$`f(𝐯)`$ $`=`$ $`N\left[\mathrm{exp}\left({\displaystyle \frac{v^2}{\sigma ^2}}\right)\mathrm{exp}\left({\displaystyle \frac{v_{\mathrm{esc}}}{\sigma ^2}}\right)\right]`$ (12)
$`\times \mathrm{exp}\left(\lambda {\displaystyle \frac{v_\varphi ^2+v_\mathrm{z}^2}{\sigma ^2}}\right),`$
where $`v_{\mathrm{esc}}`$ is the escape velocity and $`N`$ is a normalisation constant. The deviation of the velocity distribution from isotropic is parameterised by $`\lambda `$, the standard Maxwellian halo, cut-off at the escape velocity, is recovered for $`\lambda =0`$. Since the halo is formed by gravitational collapse the deviation is most likely to be towards radial orbits, with $`0<\lambda <1`$ . We take $`v_{\mathrm{esc}}=650\mathrm{k}\mathrm{m}\mathrm{s}^1`$, although the effect of variations in the value of the escape velocity is small .

The annual variation in $`T(E)`$, $`\mathrm{\Delta }T(E)`$, is plotted in Fig. 1 for $`\lambda =0,0.5,1`$. The mean value of $`T(E)`$ is slightly higher for the asymmetric models than for the standard Maxwellian halo model.

### B Power-law halos

Evans’ family of axisymmetric distribution functions  lead to velocity distribution functions, in the rest frame of the galaxy, of the form :

$`f(𝐯)`$ $`=`$ $`(AR^2v_\varphi ^2+B){\displaystyle \frac{\mathrm{exp}\left[2(v/v_0)^2\right]}{(R^2+R_\mathrm{c}^2)^2}}`$ (14)
$`+{\displaystyle \frac{\mathrm{exp}\left[2(v/v_0)^2\right]}{R^2+R_\mathrm{c}^2}},`$
with

$`A`$ $`=`$ $`\left({\displaystyle \frac{2}{\pi }}\right)^{5/2}{\displaystyle \frac{(1q^2)}{Gq^2v_0^3}},B=\left({\displaystyle \frac{2}{\pi ^5}}\right)^{1/2}{\displaystyle \frac{R_\mathrm{c}}{Gq^2v_0}},`$ (15)
$`C`$ $`=`$ $`{\displaystyle \frac{2q^21}{4\pi Gq^2v_0}},`$ (16)
where $`R_\mathrm{c}`$ is the core radius, $`R_0`$ is the solar radius and $`q`$ is a flattening parameter, which varies between 1, for spherical halo, and $`1/\sqrt{2}=0.707`$. Following Ref.  we take $`R_\mathrm{c}=8.5`$kpc and $`R_0=7`$kpc and explore the effect of varying $`q`$.

The annual variation in $`T(E)`$, $`\mathrm{\Delta }T(E)`$, is plotted in Fig. 2 for $`q=1,0.85,0.707`$. The mean value of $`T(E)`$ is slightly lower for the flattened halos than for the Maxwellian halo. The change in the annual variation due to flattening is larger than that due to an asymmetric velocity distribution, particularly so for large recoil energies and small WIMP masses.

### C Bulk rotation

Halo models with bulk rotation can be constructed by taking linear combinations of the velocity distribution function :

$$f_{\mathrm{rot}}(𝐯)=a_{\mathrm{rot}}f_+(𝐯)+(1a_{\mathrm{rot}})f_{}(𝐯),$$
(17)
where

$$f_+(𝐯)=\{\begin{array}{cc}f(𝐯)v_\varphi >0\hfill & \\ 0v_\varphi <0\hfill & \end{array}$$
(18)
$$f_{}(𝐯)=\{\begin{array}{cc}0v_\varphi >0\hfill & \\ f(𝐯)v_\varphi <0\hfill & \end{array}$$
(19)
and $`a_{\mathrm{rot}}`$ is related to the dimensionless galactic angular momentum, $`\lambda _{\mathrm{rot}}`$: $`\lambda _{\mathrm{rot}}=0.36|a_{\mathrm{rot}}0.5|`$. Numerical studies of galaxy formation find that $`|\lambda _{\mathrm{rot}}|<0.05`$ , corresponding to $`0.36<a_{\mathrm{rot}}<0.64`$. A non-rotating halo has $`a_{\mathrm{rot}}=0.5`$, whilst a counter-rotating (co-rotating) has $`a_{\mathrm{rot}}<0.5`$ ($`a_{\mathrm{rot}}>0.5`$).

The annual variation in $`T(E)`$, $`\mathrm{\Delta }T(E)`$, is plotted in Fig. 3 for Maxwellian halos with $`a_{\mathrm{rot}}=0.36,0.5,0.64`$. The mean value for counter(co)-rotation is lower (higher) than for the non-rotating Maxwellian halo at small recoil energies, and higher (lower) at large recoil energies. The change in the annual variation, which is suppressed (enlarged) for counter(co)-rotation, is largest for small recoil energies.

For each of the halo models studied the annual variation is negative for small recoil energies ($`E4`$ keV for $`M_\chi =30`$GeV, $`E30`$ keV for $`M_\chi =200`$ GeV) increasing as the recoil energy increases. For fixed recoil energy the annual variation is largest for small WIMP masses.

## IV Analysis of the DAMA data

Brhlik and Roszkowski  have devised a technique for comparing the experimental results released by DAMA with theoretical predictions for the annual modulation signal, in the absence of detailed information about the experimental set-up, such as the efficiency of each NaI crystal. They define a function $`\kappa `$:

$$\kappa =\mathrm{\Sigma }_k\frac{(S_{0,\mathrm{k}}^{\mathrm{th}}S_{0,\mathrm{k}}^{\mathrm{exp}})^2}{\sigma _{0,\mathrm{k}}^2}+\mathrm{\Sigma }_k\frac{(S_{\mathrm{m},\mathrm{k}}^{\mathrm{th}}S_{\mathrm{m},\mathrm{k}}^{\mathrm{exp}})^2}{\sigma _{\mathrm{m},\mathrm{k}}^2},$$
(20)
where the experimental errors on the time dependent and independent parts of the signal, $`\sigma _{\mathrm{m},\mathrm{k}}`$ and $`\sigma _{0,\mathrm{k}}`$ respectively, serve as weights. The contour, in the $`m_\chi \xi \sigma _\mathrm{p}`$ plane, $`\kappa =35`$ agrees reasonably well with the DAMA collaborations $`3\sigma `$ contour. Whilst this approach does not give accurate confidence limits on $`m_\chi `$ and $`\xi \sigma _\mathrm{p}`$ it does illustrate the qualitative effect of varying the properties of the halo model on the values of the WIMP parameters obtained from the data. A full likelihood analysis of the DAMA data has been carried out for models with bulk rotation in Ref. , allowing us to check the reliability of Brhlik and Roszkowski’s technique.

Given the experimental difficulties of extracting a small annual variation from, possibly time dependent, backgrounds, concerns have been expressed about the interpretation of the earlier (1 and 2 year) DAMA data as evidence for a WIMP signal . Furthermore the Cryogenic Dark Matter Search (CDMS) collaboration have recently released limits on the WIMP-Nucleon cross section which exclude, at 85$`\%`$ confidence, the entire DAMA $`3\sigma `$ allowed region . The DAMA collaboration have, however, performed a thorough analysis of the various sources of possible systematic errors  and, since the DAMA and CDMS experiments use different target nuclei, assumptions are required to perform a direct comparison. In any case even if the interpretation of the DAMA annual modulation signal as WIMP scattering is eventually found to be erroneous, our results will still indicate the qualitative effect of these non-standard halo models on the analysis of a WIMP annual modulation signal.

In Figs. 46 we plot contours of $`\kappa =35`$, delineating the region of $`m_\chi `$-$`\xi \sigma `$ parameter space compatible with the DAMA annual modulation signal, for asymmetric, flattened and rotating halo models respectively. In Fig. 4 we see that as the asymmetry of the velocity distribution is increased the allowed region is enlarged and moves to larger masses and slightly smaller interaction cross sections. We saw in Fig. 2 that the change in the annual modulation signal for the flattened halo model is larger than that for the asymmetry velocity distribution. Consequentially the change in the allowed region for the flattened halo model is smaller, extending to larger masses and slightly larger interaction cross sections. In Fig. 6 we can see that counter-rotation ($`a_{\mathrm{rot}}<0.5`$) contracts the allowed region and shifts it to smaller cross sections whilst co-rotation ($`a_{\mathrm{rot}}>0.5`$) expands it and shifts it to larger cross sections. This is in good agreement with the results of Ref.  on the effect of bulk rotation, which were found via a full likelihood analysis of the DAMA data.

## V Conclusions

In this paper we have examined the effect of an asymmetric halo velocity distribution, halo flattening and bulk halo rotation on the WIMP annual modulation signal. Whilst for each of the halo models the change in the mean signal is small, $`<10\%`$ for the range of recoil energies probed by the DAMA experiment, the amplitude of the annual modulation can change significantly. The magnitude of the change depends on the WIMP mass and recoil energy, as well the nature of the deviation of the halo model from the standard Maxwellian halo. With bulk rotation the change in the annual variation, relative to the Maxwellian halo model, is largest for small recoil energies, whilst halo flattening produces the largest change for larger recoil energies.

The range of WIMP masses consistent with the DAMA annual modulation signal is enlarged significantly, to roughly $`30<m_\chi <150`$ GeV at $`3\sigma `$, for each of the models considered. The shape of the allowed region in the $`m_\chi \xi \sigma _\mathrm{p}`$ plane is different for each model however. This indicates that the uncertainties in halo modelling have a significant effect on the WIMP mass determined from the annual modulation signal. Therefore more sophisticated halo models (such as those in Ref. ) need to be developed and used in the analysis of annual modulation data.

###### Acknowledgements.

The author is supported by PPARC and acknowledges the use of the Starlink facilities at QMW. The author would also like to thank J. D. Vergados for useful discussions.

## A Interaction Cross Section

The differential cross section for neutralino scattering off a target nucleus is dominated by the scalar and axial terms :

$$\frac{\mathrm{d}\sigma }{\mathrm{d}𝐪^2}=\frac{\mathrm{d}\sigma ^{\mathrm{scalar}}}{\mathrm{d}𝐪^2}+\frac{\mathrm{d}\sigma ^{\mathrm{axial}}}{\mathrm{d}𝐪^2},$$
(A1)
where $`𝐪=(m_\mathrm{A}m_\chi )/(m_\mathrm{A}+m_\chi )𝐯`$ is the momentum transfered.

The scalar differential cross section, which arises due to Higgs boson and squark exchange, is given by

$$\frac{\mathrm{d}\sigma ^{\mathrm{scalar}}}{\mathrm{d}𝐪^2}=\frac{1}{\pi v^2}\left[Zf_\mathrm{p}+(AZ)f_\mathrm{n}\right]^2F^2(q),$$
(A2)
where $`f_\mathrm{p}`$ and $`f_\mathrm{n}`$ are the effective neutralino couplings to the proton and the neutron respectively and $`F(q)`$ is the scalar nuclear form factor. The form factor for Na is usually taken to be equal to unity, whilst for I the Saxon-Woods form factor

$$F(q)=\frac{3j_1(qR_1)}{qR_1}\mathrm{exp}\left[(qs)^2/2\right],$$
(A3)
where $`R_1=\sqrt{R_\mathrm{A}^25s^2}`$, $`R_\mathrm{A}=A^{1/3}\times 1.2`$fm and $`s=1`$fm, is used.

The axial differential cross section, which arises due to $`Z_0`$ and squark exchange, is given by

$$\frac{\mathrm{d}\sigma ^{\mathrm{axial}}}{\mathrm{d}𝐪^2}=\frac{8}{\pi v^2}\mathrm{\Lambda }^2J(J+1)S(q),$$
(A4)
where $`J`$ is the total angular momentum, $`S(q)`$ is the spin form factor and $`\mathrm{\Lambda }`$ depends on the axial couplings of the neutralino to the quarks (see Ref.  for more details and explicit expressions).

The differential event rate for a given detector can then be expressed as :

$`{\displaystyle \frac{\mathrm{d}R}{\mathrm{d}E}}`$ $`=`$ $`{\displaystyle \frac{4}{\pi ^{3/2}}}{\displaystyle \frac{\rho _\chi }{m_\chi }}T(E)`$ (A6)
$`\times \{\left[Zf_\mathrm{p}+(AZ)f_\mathrm{n}\right]^2F^2(q)+8\mathrm{\Lambda }^2J(J+1)S(q)\}.`$


format: text

text:
# Stellar Mass-to-Light Ratios and the Tully-Fisher relation

## 1. Galaxy Evolution Models and Mass-to-Light Ratios

We have used the galaxy evolution models described by Bell & Bower (2000) to investigate stellar $`M/L`$ ratios of galaxies. These models were tuned to fit the observed trends between the colors and the structural parameters of spiral galaxies (Bell & de Jong 2000). Using a local gas density dependent star formation law, the photometric evolution is calculated, taking chemical evolution into account. As well as a closed box model we have models with gas infall and outflow, mass dependent formation epochs and star bursts. All models show large variations in $`M/L`$, amounting from a factor 8 in $`B`$ to 2 in $`K`$, but in all models we find a strong correlation between $`M/L`$ and optical color (e.g. Fig. 1a). The slope of the color-$`M/L`$ relation is very robust against the particular stellar population synthesis model and against the exact details of the galaxy evolution model. The main uncertainty in the correlation is the zero-point, which is determined by the assumed IMF; most notably by the relative amount of low mass stars, which contribute to the mass but not to the luminosity and color.

A constraint on the color-$`M/L`$ correlation zero-point can be obtained from galaxy rotation curves. The stellar disk in a galaxy cannot be more massive than allowed by its rotation curve, resulting in a maximum-disk $`M/L`$. In Fig. 1b we show the maximum-disk $`M/L`$ values versus the extinction corrected $`B`$$`R`$ color for the Verheijen (1998) galaxy sample. These $`M/L`$ values are truly upper limits: any mass not accounted for in the rotation curve decompositions will push the stellar $`M/L`$ even lower. The solid line in Fig. 1b shows the fit to the color-$`M/L`$ relation for our best model using a standard Salpeter IMF. Clearly this model over-predicts the maximum allowed mass for many galaxies, as many galaxies lie below the line. Using a Salpeter IMF with a flat slope below 0.6 $`M_{}`$ as suggested by recent observations results in the dotted line which is consistent with the observations (for $`H_0=70`$ km s-1 Mpc-1 or $`D_{\mathrm{UrsaMajor}}=20`$ Mpc).

## 2. Tully-Fisher Relations

Observed TF-relations are known to have a passband dependence, both in slope and in zero-point. Applying Tully et al. (1998) extinction corrections and our color-$`M/L`$ correlations, we can calculate stellar mass TF-relations from the observed TF-relations. We find that the stellar mass TF-relations derived from the different passbands are equal to within the uncertainties. By adding in the HI gas mass we can calculate the baryonic TF-relations (Fig. 1c). We find that the slope of the baryonic TF-relation must be less then $`3.5\pm 0.2`$, significantly lower than found by McGaugh et al. (2000), mainly due to our use of stellar M/Ls consistent with maximum disk constraints and our exclusion of low luminosity dwarfs with poorly determined inclinations and rotation velocities.

### Acknowledgments.

Support for RSdJ was provided by NASA through Hubble Fellowship grant #HF-01106.01-A from the Space Telescope Science Institute, which is operated by the AURA, Inc., under NASA contract NAS5-26555. Support for EFB was provided by NASA LTSA grant NAG5-8426.

## References

Bell, E. F., & Bower, R. G. 2000, MNRAS in press (astro-ph/0005294)

Bell, E. F., & de Jong, R. S. 2000, MNRAS, 312, 497

McGaugh, S. S. et al. 2000, ApJ, 533, L99

Tully, R. B. et al. 1998, AJ, 115, 2264

Verheijen, M. A. W. 1997, PhD thesis, Univ. of Groningen


text:
# HST Optical-NIR Colors of Nearby 𝑅^{1/4} and Exponential Bulges1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by Association of Universities for Research in Astronomy, Inc. (AURA), under NASA contract NAS5-26555

## 1 Introduction

The size of the central bulge compared to the size of the disk is one of the classification criteria of the entire Hubble sequence (Sandage & Bedke 1994). Furthermore, bulges host a large fraction of the baryons which have been converted into stars during the Hubble time (e.g., Fukugita, Hogan & Peebles 1998). It is also likely that the formation of bulges plays a role in polluting the intergalactic medium with processed matter and radiation, and thus affects the formation of luminous structure on larger scales. Yet, to date, it is not understood not only ‘how’ bulges form, but even whether they form before, contemporaneously with, or after the disks (e.g., see the detailed review by Wyse, Gilmore & Franx 1997, and also Bouwens, Cayon & Silk 1999). In the one extreme, bulges may arise from the collapse of a primordial gas cloud into clumps which then merge together, and the disk be a later accretion by gas infall (e.g., Eggen et al. 1962; Larson 1975; Carlberg 1984). In the opposite extreme, bulges may form by vertical dynamical instabilities of the disks (e.g., Combes et al. 1990; Pfenniger 1993).

Historically, ‘bulges’ have been recognized as such because they appear as high surface brightness, concentrated objects in the considerably larger and fainter disks. These ‘bulges by definition’ often also have quite a distinct color, and, as the name suggests, are considerably fatter than the disks. Superficially, they appear to share many properties with small elliptical galaxies: The relatively few and typically rather massive disk-embedded spheroids which have been studied in detail show an $`R^{1/4}`$ radial light profile, relatively old and metal-enriched stellar populations, and are rotationally-supported systems (Wyse et al. 1997 and references therein). Despite the temptation of historical inertia, clearly wanting us to restrict the use of the name ‘bulge’ to this sort-of-well-defined class of stellar systems only, it is however unwise, in the quest for clarifying the origin of these fundamental components of one-Hubble-time-old galaxies, to ignore that different kinds of stellar structures are also often found in the centers of the local-universe disks (Wyse et al. 1997).

Significant differences between the ‘bulges’ of intermediate- to late-type disks and the more massive spheroids have in fact been found. First, preliminary studies suggest that the stellar populations of small bulges in general may be younger and more metal-poor than those of the large bulges (Trager et al. 1999). Furthermore, several studies have revealed ‘pseudo-bulges’ with e.g., cold kinematics (Kormendy 1993), a peanut-shape morphology associated with bar-like kinematics (e.g., Kuijken & Merrifield 1995; Bureau & Freeman 1999) or an exponential – rather than $`R^{1/4}`$ – fall-off of the light distribution (hereafter exponential bulges; e.g., Kormendy & Bruzual 1978; Shaw & Gilmore 1989; Courteau et al. 1996 and references therein).

The pseudo-bulges also appear as clearly distinct structures from the disks, and replace the ‘classical’ bulges in disks that, at least superficially, look rather similar to others that host instead the canonical $`R^{1/4}`$ structures. It is unclear how these strange central structures relate to their $`R^{1/4}`$ relatives. However, the fact that they ‘substitute for’ the $`R^{1/4}`$ bulges in otherwise normal intermediate-type disks provides an argument for including them in any thorough exploration which is aimed at clarifying the origin and evolution of the spheroidal components of the local disk galaxies. These pseudo-bulges may lead us to understand by comparison which are the required circumstances for forming a dense elliptical-like spheroid (failing which, only a pseudo-bulge is formed), or may indicate the occurrence of fundamental changes with cosmic time in the evolutionary paths of galaxies, or may turn out to be consistent with being evolutionary related to the ‘classical’ bulges for which we seek the origin and progenitors.

In order to help clarify the origin of the structural diversity between the largest and the smallest disk-embedded spheroids, we performed a Hubble Space Telescope (HST) snapshot survey of $`80`$ spiral galaxies randomly-selected out of a complete sample of 134 targets. The central regions of the selected objects were imaged in the visual (with WFPC2 and the F606W filter) and in the near-infrared (NIR; with NICMOS and mostly the F160W and occasionally the F110W filter). The observations and the data reduction are fully described in a series of papers, separately for the WFPC2 survey (Carollo et al. 1997; Carollo et al. 1998; Carollo & Stiavelli 1998; Carollo 1999) and the NICMOS survey (Carollo et al. 2000, paper I; Seigar et al. 2000, paper II). The high angular resolution of the HST data allowed us to mask out from several of the images patches of dust and knots of star formation, and thus derive for those targets ‘clean’ isophotal models for the underlying galaxian light. Two-component ‘bulge plus disk’ analytical fits were applied to the so-derived radial light profiles. Echoing previous voices, we considered as a ‘bulge’ that distinct central structure which ‘contains all the light in excess of the inward extrapolation of a constant scalelength exponential disk’ (Wyse et al. 1997; Gilmore 1999). Within this more general definition of bulge, we found that the bulge light profile was typically best modeled with an $`R^{1/4}`$-law in the early-type disks and with an exponential law in the later types, in agreement with previous studies (e.g., Courteau et al. 1996).

In this paper we combine the NIR measurements presented in papers I and II with the WFPC2 measurements so as to investigate the optical-NIR colors of the exponential versus $`R^{1/4}`$-law bulges found in our HST survey, as well as the colors of the photometrically-distinct nuclei which we found embedded in the dynamical centers of all the exponential bulges. The objects included in this study are listed in Table 1, together with the information on whether they have an $`R^{1/4}`$-law or an exponential bulge, on the $`VH`$ color for the host bulges, and on the $`VH`$ color for the embedded nuclei. The intrinsic dust-age-metallicity degeneracy of broad-band colors for intermediate-to-old age stellar populations is a well-known problem. Nonetheless, although the calibrations of absolute ages and metallicities do remain elusive, photometric indicators can still be a powerful bench-mark for ‘ranking’ the stellar populations among diverse kinds of bulges, and among bulges and other galactic sub-components such as the photometrically-distinct nuclei. We briefly discuss plausible implications of our results for the formation of bulges during the course of cosmic history.

## 2 Results

### 2.1 The Colors of Exponential and $`R^{1/4}`$-law Bulges

The majority of the $`R^{1/4}`$-law bulges found in our sample are bright (massive) systems embedded in early-type disks, in contrast with the fainter exponential bulges which are found mostly in Sb to Sc hosts, in agreement with previous studies (see Figure 1, where the absolute $`V`$ magnitude of exponential and $`R^{1/4}`$-law bulges is plotted versus their $`VH`$ colors). Figure 2 shows the $`VH`$ color distribution for the $`R^{1/4}`$-law (solid line) and exponential (dashed line) bulges detected in our HST survey. The measurements are expressed in $`AB`$ magnitudes. The colors are obtained by integrating the smooth analytical fits to the bulge light profiles inside $`1/2R_e`$ (or within $`6^{\prime \prime }`$ if that radius exceeds the radial extension of the NICMOS data); therefore they exclude (or at least minimize) any contribution from e.g., the nuclei, strong dust lanes/patches and knots of recent star formation.

The $`VH`$ color distribution of the $`R^{1/4}`$ bulges peaks around $`<VH>1.3`$, with a sigma $`\mathrm{\Delta }(VH)0.1`$ magnitudes. Bruzual & Charlot models, calibrated into the F606W and F160W filters with Synphot/IRAF, indicate that the average metallicity of these systems must be at least half-solar for their age to be smaller than that of the universe. Assuming a solar metallicity, the average $`<VH>`$ color corresponds to an average age for the population of $`6`$ Gyrs, with a spread of about 3 Gyrs around this average. The presence of a small fraction of younger stars cannot be clearly excluded on the basis of the available data alone. In the absence of dust, the bluest colors that we observe in the R1/4 bulges could be obtained by contaminating, e.g., a 9 Gyrs old underlying population with a 10%, 1% or 0.1% in mass of stars 1 Gyr, 100 Myrs or less-than-5 Myrs old, respectively. Possible ‘solutions’ clearly exist also in the presence of dust: Using the ‘screen-approximation’ and a Cardelli extinction law (Cardelli, Clayton & Mathis 1989) as a benchmark, with e.g., an extinction in the $`V`$ band of $`A_V=3`$, about 10% (in mass) of the population should be less than 5 Myrs old in order to reproduce the bluest end of the $`R^{1/4}`$ bulges $`(VH)`$ distribution. A more modest $`A_V=1.5`$ could accomodate a fraction as high as $`70\%`$ of young stars, if these were about 100-Myrs-old. On the other hand, the smooth morphology of this class of objects and their color maps, virtually featureless at all radii, make it rather implausible that these systems contain large amounts of dust. In fact, the only features detected in the color maps for the $`R^{1/4}`$ bulges are some mostly-nuclear patches of dust. The effects of these dust patches are removed by measuring the colors by integrating over the light profiles (since the latter are obtained by means of isophotal fits in which the patches of dust are masked out). The most plausible solution for the $`R^{1/4}`$ bulges is therefore that at least more than 90% of their stars are older than several Gyrs. This finding echos previous results for the $`R^{1/4}`$-law bulges (e.g., Peletier et al. 1999 and references therein) by supporting the idea that many of these dense spheroids are, at least population-wise, similar to small ellipticals.

The $`VH`$ color distribution of the exponential bulges, in contrast, presents some surprises. The average $`<VH>`$ of these systems is in fact $`0.4`$ magnitudes bluer than the corresponding average for, and the sigma $`\mathrm{\Delta }(VH)0.4`$ magnitudes is significantly larger than that of, the $`R^{1/4}`$ bulges. Furthermore, although a few of the exponential bulges are actually as red as the reddest of the $`R^{1/4}`$-law bulges, the former extend to significantly bluer colors than the latter, down to a $`(VH)0`$. In principle, adding a small fraction of young stars to an underlying old stellar population could again explain the observed colors. For example, in the absence of dust, a 10% of 100 Myr old stars mixed to the remaining 9-Gyr-old population would be sufficient to obtain a $`(VH)0.5`$ magnitudes; the required percentage of young stars would clearly decrease for even younger stellar populations. In the presence of substantial amounts of dust, say $`A_V`$ up to 2 magnitudes, the age of the population should decrease down to $`<\mathrm{\hspace{0.25em}5}`$ Myr in order to obtain a $`(VH)0.5`$ with a 10% or smaller contamination of young stars. Still, thanks to the high angular resolution of the HST images which allow to perform an accurate masking of any localized feature in the images, the effects of dust patches and possible knots of star formation are also in this case negligible in the colors derived from the light profiles (from isophotal fits). Therefore, while a diffuse younger component well-mixed with the underlying old population cannot clearly be ruled out by the $`(VH)`$ color alone, the most appealing solution for the exponential bulges is that they are on average significantly younger than their $`R^{1/4}`$ relatives.

### 2.2 Colors of the Nuclei

About a third of our sample was imaged in all three filters $`J`$, $`H`$ and $`V`$. In Figure 3 we plot the $`JH`$ color versus the $`VH`$ color for the compact nuclei found embedded in intermediate- to late-type galaxies, including the nuclei which are found in the centers of the exponential bulges. Actually, each exponential bulge in our sample hosts one such a nucleus (Carollo 1999). In the Figure, squares are the nuclei embedded in the exponential bulges. The triangles represent instead the nuclei embedded in systems with no isophotal fit, i.e., systems whose central regions are heavily obscured by dust or complicated by recent star formation, so that for them no reliable isophotal fit could be performed (see previous papers for the corresponding list of measurements). Dust extinction affects points on this diagram by shifting them along the arrow shown in the upper-left corner of the Figure. The arrow is drawn for a 1 mag extinction in the $`V`$ band. Theoretical tracks derived from Bruzual & Charlot models for metallicities $`Z=0.02Z_{}`$ (dotted line) and $`Z=Z_{}`$ (solid line) are also shown in Figure 3. Stellar ages increase upward along the tracks from $``$100 Myr to 18 Gyr.

Many central compact nuclei (typically in hosts without an analytical fit to the bulge) have $`VH`$ and $`JH`$ colors which are incompatible with arising from unobscured stellar populations. These nuclei are embedded in very complex circum-nuclear structures, e.g., in strong dust lanes/arms or rings/arms of recent star formation. For a few of these systems information on the central spectral properties is available at ground-based resolution, i.e., at a resolution covering an area $`100`$ times larger than that covered by the nuclei, an area which includes this complex circum-nuclear structure as well (making the spectra not suitable for studying the physical properties of the central compact nuclei, for which we have to rely on the HST photometry only). On ground-based scales of $``$1kpc, these galaxies show an active-type central spectrum, either of AGN- or of HII-type. It is therefore not surprising that the optical/NIR light from their central resolved nuclei can also be polluted by several magnitudes of extinction, and/or can contain non-thermal emission from a central AGN.

In contrast, most of the central nuclei which are embedded in the exponential bulges have $`VH`$ and $`JH`$ colors which are compatible with arising from nearly unobscured stellar populations. Spectral information is not available for most of their host galaxies, not even from the ground. In principle, large amounts of dust extinction cannot be ruled out even for these nuclei. However, in contrast with the complexity of structure generally underlying the central nuclei with very red $`VH`$ and $`JH`$ colors, the exponential bulges are generally rather smooth, amorphous systems, showing not nuch evidence for the presence of large amounts of dust. This supports the idea that the optical/NIR light of the nuclei of exponential bulges is mostly stellar, i.e., it is not heavily polluted by a non-thermal component or by large amounts of dust. Under this assumption, relatively old ages – of the order of $``$ 1Gy and above – are suggested by the combined $`VH`$ and $`JH`$ colors for several of these nuclei embedded in the exponential bulges.

If spectroscopically confirmed for at least some of these nuclei, these relatively old ages would have interesting consequences for the masses of these systems. In Figure 4 we plot the absolute magnitude $`V`$ against the $`VH`$ color for all the nuclei for which both a $`V`$ and an $`H`$ image are available (for several of these a $`J`$ image is not available; this is the reason why there are more points here compared to Figure 3). Age increases from top to bottom along the plotted Bruzual & Charlot tracks; these refer to a $`2.5\times 10^6M_{}`$ star cluster of metallicities $`Z=0.02Z_{}`$ (dotted line) and $`Z=Z_{}`$ (solid line). Dust extinction shifts points along the direction indicated by the upper-left arrow. The mass of a stellar cluster and the age of its stellar population have nearly fully degenerate effects on a color-magnitude diagram: The tracks can be shifted upward (downward) by increasing (decreasing) the mass of the star cluster (a 2.5 magnitudes shift implies a factor 10 variation in mass); correspondingly, a larger (lower) stellar age will be associated with any fixed point encompassed by the tracks on the $`V`$ vs $`VH`$ diagram. In this enlarged sample, again most of the central compact nuclei which are embedded in the complex circum-nuclear structure of those hosts without an analytical fit have $`VH`$ colors which are incompatible with arising from unobscured/unpolluted stellar populations. And again most of the nuclei embedded in the exponential bulges are instead consistent with being relatively unobscured stellar structures. They are also, at face value, more enriched than e.g., the Milky Way globular clusters, since they typically lie rightward of the metal-poor stellar track. If the large ages estimated from the $`(JH)`$$`(VH)`$ diagram are assumed for (some of) these nuclei, their location on the $`V`$$`(VH)`$ plane implies masses in the thereabouts of a few 106 to a few 107 M, depending on the exact dating of the stellar structure between $`1`$ and many Gyrs.

### 2.3 Comparison between Nuclei and Host Exponential Bulges

In Figure 5 we show the comparison between the $`VH`$ color of the exponential bulges and the $`VH`$ color of their hosted nuclei. Within the error bars there is a positive trend between the colors of the nuclei and those of their host bulges (Spearman’s and Kendall’s rank order tests give a probability of correlation greater than 94% and 95%, respectively). The slope of the correlation is compatible with being unity, but the error bars are too large to determine it reliably. Different bulges-nuclei pairs cover a large range of $`VH`$ colors (consistently with what discussed in Figure 2), but most bulge-nucleus pairs seem to have possibly rather similar $`VH`$ colors. Conspiracies may be at play; on the other hand, this result suggests that the nuclei embedded in the exponential bulges may have stellar populations similar to those of their host bulges.

## 3 Discussion

It is by now a decade ago that e.g., Hasan & Norman (1990), Norman & Hasan (1990) and Pfenniger & Norman (1990), proposed some basic dynamical mechanisms that produce an amplification of the accretion rate of gas clouds into the central regions of barred galaxies, and discussed their effects on the evolution of disk galaxies and their implications for bulge formation. Only a few years later John Kormendy presented the first evidence for cold kinematics in intermediate-type bulges at the IAU Symposium 153, thereby providing extremely compelling observational “evidence that some bulges are really disks” (Kormendy 1993). Since then, several other studies have been performed that further support a close connection between bulges and disks in Sb-Sc galaxies (e.g., Kuijken & Merrifield 1995; Bureau & Freeman 1999; Courteau et al. 1996 and references therein), and thus possibly the growth of bulges inside the disks due to the secular evolution of the latter.

The theoretical suggestions for how the ‘disk-driven’ bulge formation may take place invoke the presence of non-axisymmetric perturbations, i.e. ‘bars’, generated by some kind of disk instabilities. First, the stars initially in a bar can be kicked off the plane of the disk by ‘buckling’ instabilities, and form a three-dimensional bulge (e.g., Combes et al. 1990). Furthermore, numerical experiments indicate that a very efficient way to induce significant evolution along the Hubble sequence is to accrete mass in the centers of galaxies hosting a tumbling triaxial structure such as a bar. The addition of a central mass concentration of $`1\%`$ of the total mass in fact disrupts orbits essential to its existence (e.g., Pfenniger & Norman 1990; Hasan, Pfenniger & Norman 1993; Norman, Sellwood & Hasan 1996). Bars do greatly enhance the efficiency of dissipative mass transfer to the centers of the disks (Shlossman 1994; Shlosman & Robinson 1995; Shlosman, Begelman & Frank 1990). Thus, provided there is an efficient mechanism for funneling the dissipative disk material down to the very center, the bar itself is dissolved, and its stars are deflected out of the plane of the disk, giving birth to a three-dimensional bulge-like structure.

The broad-band photometric analysis presented here is clearly afflicted by the well-known full-degeneracy of mass, age, metallicity, dust, emission lines. The arguments of plausibility which lead to the age and mass estimates for the exponential bulges and their central nuclei discussed above are indeed to be taken with care: Spectroscopic confirmation is needed. Still, within the caveats, those assumptions are a plausible solution with several implications for the formation of bulges as a function of cosmic time.

First, if dust effects are not at play, the fact that some exponential bulges appear to be as red as the elliptical-like bulges suggests that this ‘mode’ of bulge formation may have been active even at the early stages in our universe. However, the fact that the exponential bulges appear to be bluer on average than the elliptical-like bulges does suggest that, as a class, they are younger (and likely less metal-enriched) than the latter. Therefore, on the one hand these new observations support the idea that the growth of central structure within the disks is still going on in our universe. On the other hand, this growth seems to occur in the form of ‘exponential-type’ structure only, i.e., a structure quite dissimilar from the dense $`R^{1/4}`$-law relics of the early times.

The mechanisms that grow the central exponential bulgelike structures remain clearly indetermined. However, at face value, the ‘photometrically-estimated’ central masses of the central nuclei embedded in the exponential bulges are well-matched to those which would have been required, in the bar formation/dissolution scenario, to disrupt progenitors bars with masses of about a few $`10^7`$-$`10^9M_{}`$. These masses are consistent with the masses inferred for the exponential bulges from their total $`V`$ luminosities (assuming similar $`M/L_V`$ ratios for the exponential bulges and their nuclei; see also Carollo 1999).

The above arguments clearly do not prove that the bar formation/dissolution mechanism operates: For example, the central nuclei may have grown inside the exponential bulges due to dynamical friction of globular clusters, as suggested by Tremaine, Ostriker & Spitzer (1975) for the nucleus of M31. At a time $`t`$, any globular cluster within a radius $`r(Gmt/\sigma )^{1/2}`$ will fall into the center dragged by dynamical friction. For plausible values of $`m`$, the mass of the infalling globular cluster, and $`\sigma `$, the velocity dispersion of the exponential bulge,444No kinematic data are yet available for the exponential bulges, so that ‘plausible’ in this context means ‘as derived from the Faber-Jackson (1976) relation and allowing for generous ($`50\%`$) variations around the Faber-Jackson value. all the globular clusters within a few to several half-light radii of the exponential bulges should have converged to their centers in about 1 to a few Gyrs. The range of masses inferred from the broad-band colors imply that about 10–100 globular clusters should have formed within such a radial range, a value which is not ruled out given that many more globular clusters will have formed than those which would be currently observed (and given that the low central stellar densities inferred for the exponential bulges, would allow the newly-borned clusters to sink in the galaxy centers without being tidally-destroyed during their infall, see Carollo & Stiavelli 1998).

On the other hand, the fact that the nuclei appear to be typically more enriched than the Milky Way globular clusters, and the similar $`VH`$ colors of the nuclei and their host exponential bulges, does argue in favour of similar stellar populations and thus of a synchronous formation of these two galactic sub-components. A natural way to synchronously form the exponential bulge and its central nucleus would indeed be a scenario where the stars in the bulge and the central nucleus are grown together by dissipative material infalling into the galaxy center.

In summary, the high-resolution photometric analysis of intermediate-type bulges, while far from being ‘conclusive’, does indicate that ‘non-classical’ bulge formation is taking place in the universe around us, and that this bulge formation is consistent with arising from secular evolution processes within the disks. In contrast, there is growing evidence that the early-type, massive bulges not only are ‘old’, but may even be as old as the Coma cluster ellipticals, with an internal age-spread of only $`2`$ Gyr (Peletier et al. 1999). Together, these results suggest that the processes that have formed the dense elliptical-like spheroids in the centers of the early-type disks are quite dissimilar from the processes which are growing the disk-like central structures within the later-type disks today.

This research has been partially funded by Grants GO-06359.01-95A and GO-0731.02-96A awarded by STScI, and has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with NASA.


text:
# Origin of Magnetoresistance Anomalies in Antiferromagnetic 𝑌𝐵𝑎₂𝐶𝑢₃𝑂_{6+𝑥}

## 1 Introduction

The transport properties of a single hole in a strongly correlated antiferromagnetically ordered quasi-2D cuprate have been the topic of much debate, both theoretically and experimentally. Unusual magnetoresistance anomalies in the heavily underdoped antiferromagnetic YBa2Cu3O6+x ($`x=0.30;0.32`$) crystals were reported recently by Y. Ando et al. . The in-plane resistivity $`\rho _{ab}`$ exhibits unconventional metal-dielectric duality with the high-temperature ($`T>50`$K) metal-like behavior in contrast with the low-$`T`$ insulating one which is not compatible both with that for a simple band insulator and for an Anderson insulator.

The crystals demonstrate an unusual behavior of the in-plane magnetoresistance, $`\mathrm{\Delta }\rho _{ab}/\rho _{ab}`$ when the magnetic field $`\stackrel{}{H}`$ is applied along the CuO2 planes. These are a striking d-wave shaped ($`\mathrm{cos}2\varphi `$) angular dependence, anomalous low-field behavior with saturation above a well-defined threshold field, and hysteretic effects at low temperatures. As the temperature decreases below $`2025`$K, the magnetic field dependence of $`\rho _{ab}`$ becomes essentially irreversible, and the system acquires a memory. The application of the field results in a persistent change in the resistivity. The authors consider qualitatively these features to be a manifestation of the ’charge stripe’ ferromagnetic structure in this system, which could be easily rotated by a rather small external magnetic field. As the temperature is lowered, it is expected that the stripe dynamics slows down and the appropriate texture in the CuO2 layers is frozen, forming somewhat like a cluster spin glass.

The antiferromagnetic domain structure in insulating YBa2Cu3O6+x ($`x<\mathrm{\hspace{0.17em}0.15}`$) with 1% Gd3+ as a ESR probe has been recently studied in detail by A. Janossy et al. . According to the ESR data the easy axes of the antiferromagnetic order are along , and  are hard axes in the ($`a,b`$) plane. In zero magnetic field the single crystal consists of the equal amounts of domains oriented along two possible easy axes. External magnetic field in the ($`a,b`$) plane induces a spin-flop reorientation transition, and at $`h>h_{sf}=\mathrm{\hspace{0.17em}5}`$T ($`T\mathrm{\hspace{0.17em}20}`$K) practically all domains arrange perpendicular to the field. The authors  consider the two models of domain texture: 1) magnetic domains separated by charge domain walls within the ($`a,b`$) plane; and 2) perfect antiferromagnetic plane domains separated by defective ($`a,b`$) planes, and they conclude their ESR data qualitatively evidence in favor of the second model.

In our opinion, the experimental data  are unlikely to be consistent with a conventional model of the well isolated spin and orbital Zhang-Rice (ZR) singlet $`{}_{}{}^{1}A_{1g}^{}`$  believed to be a ground state of the hole-doped CuO4 center in the CuO2 layers. Here, it should be noted that when speaking of a Zhang-Rice singlet as being ’well isolated’, one implies that the $`{}_{}{}^{1}A_{1g}^{}`$ ground state for the CuO4 plaquette with the two holes of the $`b_{1g}(d_{x^2y^2})`$ symmetry is well separated from any other excited two-hole states. Indeed, a simple s-like spin and orbital symmetry implies the tetragonally isotropic ’s-wave’ transport properties of the well isolated ZR singlet in the CuO2 layers. The unconventional d-wave magnetoresistance anisotropy displayed by insulating cuprates, as well as many other experimental data and theoretical model considerations evidence in favor of the more complicated structure of the valence multiplet for the hole-doped CuO4 center rather than simple ZR singlet albeit namely the latter is a guideline in the overwhelming majority of current model approaches.

The nature of the valent hole states in doped cuprates is considered as being of great importance for the high-$`T_c`$ problem. Having solved the problem we could justify the choice of the relevant effective Hamiltonian together with the opportunities of a mapping to the single band $`tJ`$ or Hubbard model.

Below we show that the model of the valence $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1}E_{u}^{}`$ multiplet, developed in Refs., provides a consistent explanation of the d-wave magnetoresistance anisotropy in the CuO2 layers of the doped cuprates believed to be the most important result of Ref..

## 2 The model of the valence $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1}E_{u}^{}`$ multiplet

The model implies a quasi-degeneracy in the ground state of the two-hole CuO$`{}_{}{}^{5}{}_{4}{}^{}`$ center with the two close in energy $`{}_{}{}^{1}A_{1g}^{}`$ and $`{}_{}{}^{1}E_{u}^{}`$ terms of $`b_{1g}^2`$ and $`b_{1g}e_u`$ configurations, respectively. In other words, one implies two near equivalent locations for the additional hole, either to the Cu3dO2p hybrid $`b_{1g}(d_{x^2y^2})`$ state to form ZR singlet $`{}_{}{}^{1}A_{1g}^{}`$, or to purely oxygen nonbonding doublet $`e_ux,y`$ state with peculiar Cu2+-Cu3+ valence resonance. The electron density distribution for two valent hole states $`b_{1g}`$ and $`e_u`$, respectively, together with the quantitative picture of the energy spectrum for valence multiplet is shown in Fig.1. It should be noted that the symmetry of the O2pe$`{}_{u}{}^{}x,y`$ states coincides with that of for Cu4p$`x,y`$ states.

In a sense, the valence $`(b_{1g}^2){}_{}{}^{1}A_{1g}^{}(b_{1g}e_u){}_{}{}^{1}E_{u}^{}`$ multiplet for the hole CuO$`{}_{}{}^{5}{}_{4}{}^{}`$ center implies an unconventional state with $`Cu`$ valence resonating between Cu3+ and Cu2+, or ’ionic-covalent’ bonding . In other words, the CuO4 center with the valence $`(b_{1g}^2){}_{}{}^{1}A_{1g}^{}(b_{1g}e_u){}_{}{}^{1}E_{u}^{}`$ multiplet represents a specific version of the ’correlation’ polaron, introduced by Goodenough and Zhou .

The model is supported both by local-density-functional calculations , ab initio unrestricted Hartree-Fock self-consistent field MO method (UHF-SCF) for copper-oxygen clusters , and a large variety of experimental data. To the best of our knowledge the one of the first quantitative conclusions on a competitive role of the hybrid copper-oxygen $`b_{1g}(d_{x^2y^2})`$ orbital and purely oxygen O2pπ orbitals in the formation of valent states near the Fermi level in the CuO2 planes has been made by A.K. McMahan et al.  and J. Tanaka et al. . Namely these orbitals, as they state, define the low-energy physics of copper oxides.

One of the most exciting experimental evidences in favor of the model with the valence $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1}E_{u}^{}`$ multiplet is associated with the observation of the midinfrared (MIR) absorption bands which polarization features are compatible with those for $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1}E_{u}^{}`$ intra-multiplet dipole transitions . The corresponding transition energies ($`\mathrm{\Delta }_{AE}`$) observed for various cuprates are of the order of a few tenths of eV, that yields a typical energy scale for the valence multiplet.

The $`e_u`$ hole can be coupled with the $`b_{1g}`$ hole both antiferro- and ferromagnetically. This simple consideration indicates clearly a necessity to incorporate in the valence multiplet both the spin singlet $`(b_{1g}e_u){}_{}{}^{1}E_{u}^{}`$ and the spin triplet $`(b_{1g}e_u){}_{}{}^{3}E_{u}^{}`$, which energy could be even lower due to ferromagnetic $`b_{1g}e_u`$ exchange. Indeed, the low-lying spin triplet state for the two-hole CuO$`{}_{}{}^{5}{}_{4}{}^{}`$ center was detected by 63,65Cu NQR in La2Cu0.5Li0.5O4 with a singlet-triplet separation $`\mathrm{\Delta }_{ST}=0.13`$ eV . The indirect manifestations of O2p$`\pi `$, or $`e_u`$ valent states were detected in the Knight shift measurements by NMR for 123-YBaCuO system . In connection with the valence $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1,3}E_{u}^{}`$ multiplet model for copper oxides one should note and comment the results of paper by Tjeng et al. , where the authors state that they ”are able to unravel the different spin states in the single-particle excitation spectrum of antiferromagnetic CuO and show that the top of the valence band is of pure singlet character, which provides strong support for the existence and stability of Zhang-Rice singlets in high-$`T_c`$ cuprates”. However, in their photoemission studies they made use of the Cu2p3/2($`L_3`$) resonance condition that allows to detect unambiguously only copper photo-hole states, hence they cannot see the purely oxygen photo-hole $`e_u`$ states.

It should be noted that the complicated $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1,3}E_{u}^{}`$ structure of the valence multiplet for the two-hole CuO$`{}_{}{}^{5}{}_{4}{}^{}`$ center could be revealed in the photoemission spectra, all the more that the odd $`{}_{}{}^{1}E_{u}^{}`$ terms play here a principal role, namely these yield a nonzero contribution to the ARPES for $`\stackrel{}{k}=0`$, or in other words at $`\mathrm{\Gamma }`$ point. In this connection one should note the experimental measurements of the photoemission spectra in Sr2CuO2Cl2  and Ca2CuO2Cl2 . All these clearly detect a nonzero photocurrent intensity in the BZ center, thus supporting the $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1,3}E_{u}^{}`$ structure of the ground state valence multiplet. Overall, the model of extended $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1,3}E_{u}^{}`$ valence multiplet allows to explain consistently many puzzling properties both of insulating and superconducting cuprates: midinfrared absorption bands , (pseudo)Jahn-Teller (JT) effect and related phenomena , spin properties .

## 3 Anomalous magnetotransport in $`b_{1g}e_u:{}_{}{}^{3}E_{u}^{}`$ hole state

The pseudo-JT polaronic nature of the spin-singlet $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1}E_{u}^{}`$ ground state  favors their localization. In addition, one should account for the antiferromagnetic background which leads to the crucial enhancement of the effective mass for the moving spin singlets. So, a spin-singlet small pseudo-JT polaron as a hole ground state is likely to be immobile. In such a situation the most effective channel for the hole transport could be related to the low-lying excited spin-triplet $`b_{1g}e_u:{}_{}{}^{3}E_{u}^{}`$ term. This gives rise to a thermo-activated hole conductivity actually observed in most of slightly doped cuprates.

In order to obtain magnetoresistivity effect we’ll consider spin and spin-orbital interactions for the spin-triplet $`b_{1g}e_u:{}_{}{}^{3}E_{u}^{}`$ state. The $`e_u`$ hole is strongly exchange-coupled both with the $`b_{1g}`$ hole on the same CuO4 center and with the nearest neighboring CuO4 centers. The spin state of the isolated CuO4 center with the $`b_{1g}e_u`$ hole configuration is described by two spin operators

$$\stackrel{}{S}=\stackrel{}{s}_{b_{1g}}+\stackrel{}{s}_{e_u},\stackrel{}{V}=\stackrel{}{s}_{b_{1g}}\stackrel{}{s}_{e_u}$$
(1)
($`\stackrel{}{S}^2+\stackrel{}{V}^2=\mathrm{\hspace{0.17em}3}/2,(\stackrel{}{S}\stackrel{}{V})=0`$), and corresponding order parameters . For the description of the orbital $`E_u`$ doublet one could use the pseudospin $`s=1/2`$ formalism with the Pauli matrices $`\sigma _{x,y,z}`$ having simple transformation properties within the $`E_ux,y`$ doublet:111We make use of standard notations for the $`D_{4h}`$ point group representations.

$$\sigma _x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)b_{2g},\sigma _y=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)a_{2g},\sigma _z=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)b_{1g}.$$

Then the effective spin-Hamiltonian for the spin-triplet $`{}_{}{}^{3}E_{u}^{}`$ state of the hole center can be represented as follows

$$\widehat{H}_S=\lambda S_z\sigma _y+DS_z^2+a(S_x^2S_y^2)\sigma _z+b\stackrel{~}{S_xS_y}\sigma _x\mu _B\stackrel{}{h}\widehat{g}_S\stackrel{}{S}2\mu _B\stackrel{}{H}_{ex}\stackrel{}{S},$$
(2)
where $`\stackrel{~}{S_xS_y}=1/2(S_xS_y+S_yS_x)`$, $`\lambda `$ is an effective spin-orbital coupling parameter for the $`{}_{}{}^{3}E_{u}^{}`$ term, $`D,a,b`$ the spin anisotropy constants, $`\widehat{g}_S`$ the effective $`g`$-tensor, $`\stackrel{}{H}_{ex}`$ the internal effective spin exchange field, $`\stackrel{}{h}`$ the external magnetic field. Hereafter, we will assume for simplicity the ideal planar CuO2 layers, isotropic spin $`g`$-factor and the planar direction of magnetic fields. In the frames of the strong molecular field approximation the spin operators can be replaced by appropriate averages

$$S_x^2S_y^2=\mathrm{cos}2\mathrm{\Phi },\stackrel{~}{S_xS_y}=\mathrm{sin}2\mathrm{\Phi },$$

and the spin-Hamiltonian (2) transforms into an effective Hamiltonian of the spin-induced low-symmetry crystalline field

$$\widehat{H}_S=a\mathrm{cos}2\mathrm{\Phi }\sigma _z+b\mathrm{sin}2\mathrm{\Phi }\sigma _x,$$
(3)
where $`\mathrm{\Phi }`$ is the azimuthal angle of the orientation for the total magnetic field $`\stackrel{}{H}=\stackrel{}{H}_{ex}(\stackrel{}{h})+\stackrel{}{h}`$. Eigenvectors and eigenvalues for such a simple Hamiltonian are

$`\mathrm{\Psi }_+=\mathrm{cos}\alpha |x+\mathrm{sin}\alpha |y,\mathrm{\Psi }_{}=\mathrm{sin}\alpha |x\mathrm{cos}\alpha |y`$ (4)
$`\mathrm{tan}2\alpha ={\displaystyle \frac{b}{a}}\mathrm{tan}2\mathrm{\Phi },E_\pm =\pm \mathrm{\Delta }(\mathrm{\Phi }),\mathrm{\Delta }(\mathrm{\Phi })=[b^2+(a^2b^2)\mathrm{cos}^22\mathrm{\Phi }]^{\frac{1}{2}},`$
respectively, where $`|x,y|E_ux,y`$. The quantum-mechanic and thermodynamic averages for $`\sigma _{z,x}`$, describing the orbital polarization (quadrupole ordering) effect, are

$`\sigma _z_\pm =\mathrm{\Psi }_\pm |\sigma _z|\mathrm{\Psi }_\pm =\pm \mathrm{cos}2\alpha ,`$
$`\sigma _z=\rho _+\sigma _z_++\rho _{}\sigma _z_{}={\displaystyle \frac{a\mathrm{cos}2\mathrm{\Phi }}{\mathrm{\Delta }(\mathrm{\Phi })}}\mathrm{tanh}\beta \mathrm{\Delta }(\mathrm{\Phi }),`$ (5)
$`\sigma _x_\pm =\pm \mathrm{sin}2\alpha ,\sigma _x={\displaystyle \frac{b\mathrm{sin}2\mathrm{\Phi }}{\mathrm{\Delta }(\mathrm{\Phi })}}\mathrm{tanh}\beta \mathrm{\Delta }(\mathrm{\Phi }),`$
respectively. Here, $`\rho _\pm `$ is the statistical weight of the $`\mathrm{\Psi }_\pm `$ states, $`\beta =1/kT`$. One should note a specific $`d`$-wave ($`d_{x^2y^2}`$ and $`d_{xy}`$) angular $`\mathrm{\Phi }`$-dependence for the thermodynamic averages $`\sigma _z`$ and $`\sigma _x`$, respectively (see Fig.2).

So, in frames of our approximation, the orbital state of the spin triplet $`{}_{}{}^{3}E_{u}^{}`$ is easily governed by the external magnetic field, thus providing an effective “magneto-orbital” transformation. Obviously, this spin-induced orbital polarization has considerable magnitude only at rather low temperatures. In the high-temperature limit ($`|a|,|b|kT`$):

$$\sigma _z\frac{a}{kT}\mathrm{cos}2\mathrm{\Phi }.$$
(6)
It should be noted that Hamiltonian (3) possess an overall tetragonal symmetry, however, this could describe the effects of spontaneous breaking of this symmetry separately in spin and orbital subspaces, or the effects of the spin-induced breaking of tetragonal symmetry in the space of orbital $`E_u\pm `$ states.

Assuming the only nonzero (Cu3d-O2p)σ-bonding contribution to the $`e_u`$ hole transfer between the $`E_u`$ states on the nearest neighboring CuO4 plaquettes, we can represent the matrix of the $`e_ue_u`$ hole transfer integral $`\widehat{t}(e_ue_u)`$ on the $`|x,y`$ basis set as follows

$$\widehat{t}_a(e_ue_u)=t_\sigma \left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right),\widehat{t}_b(e_ue_u)=t_\sigma \left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)$$

for the transfer in the  and  directions, respectively. Here, $`t_\sigma =t_a(e_uxe_ux)=t_b(e_uye_uy)`$. In frames of the tetragonal anzatz the $`e_u`$ hole energy spectrum in the tight binding approximation with the $`nn`$ transfer consists of two bands $`E^{x,y}`$ formed by the transfer $`e_uxe_ux`$ and $`e_uye_uy`$, respectively

$$E_𝐤^x=2t_\sigma \mathrm{cos}(k_xa),E_𝐤^y=2t_\sigma \mathrm{cos}(k_ya).$$
(7)
These 1D bands are characterized by maximally anisotropic effective mass, and describe the 1D hole motion in $`a`$\- and $`b`$\- directions, respectively. Nevertheless, with account for rigorous $`xy`$ ($`ab`$) symmetry it is easy to see that the transport properties appear to be tetragonally isotropic. The spin-induced low-symmetry crystalline field (3) results in a mixing of the two bands. Indeed, instead of simple expressions for the bare hole transfer integrals we obtain matrices for the renormalized hole transfer integrals on the $`|\pm `$ basis set

$`\widehat{t}_a(e_ue_u)`$ $`=`$ $`t_\sigma \left(\begin{array}{cc}\mathrm{cos}^2\alpha & \frac{1}{2}\mathrm{sin}2\alpha \\ \frac{1}{2}\mathrm{sin}2\alpha & \mathrm{sin}^2\alpha \end{array}\right),`$
$`\widehat{t}_b(e_ue_u)`$ $`=`$ $`t_\sigma \left(\begin{array}{cc}\mathrm{sin}^2\alpha & \frac{1}{2}\mathrm{sin}2\alpha \\ \frac{1}{2}\mathrm{sin}2\alpha & \mathrm{cos}^2\alpha \end{array}\right),`$ (8)
where we have simply made use of Exps.(3) for the $`\mathrm{\Psi }_\pm `$ functions.

Thus, after simple algebra we can obtain for the renormalized bands

$`E_\stackrel{}{k}^\pm =t_+(\stackrel{}{k})\left[t_{}^2(\stackrel{}{k})+\mathrm{\Delta }^2(\mathrm{\Phi })+2\mathrm{\Delta }(\mathrm{\Phi })t_{}(\stackrel{}{k})\mathrm{cos}2\alpha \right]^{\frac{1}{2}},`$
$`t_\pm (\stackrel{}{k})=t_\sigma \left(\mathrm{cos}(k_xa)\pm \mathrm{cos}(k_ya)\right).`$ (9)
The inverse in-plane effective-mass tensors for the two bands can be written as follows

$$\widehat{\left(\frac{1}{𝐦^{}}\right)}_\pm =\frac{1}{m_0^{}}\left[\mathbf{\hspace{0.17em}1}\pm \mathrm{cos}2\alpha \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\right],$$
(10)
where $`\frac{1}{m_0^{}}=\frac{4a^2}{\mathrm{}^2}t_\sigma `$.

Thus, we obtain the spin-induced shift/splitting of the bare bands, which results in the breaking of the tetragonal $`xy`$ symmetry with appearance of the spin-dependent anisotropy in effective mass.

Taking into account the proportionality relation between the single-band conductivity and the inverse effective-mass tensor, one might obtain a surprisingly simple expression for the in-plane magnetoresistance

$$\frac{\delta \rho _{a,b}}{\rho _{a,b}}=\mathrm{cos}2\alpha =\sigma _z,$$
(11)
where the upper (lower) sign corresponds to the  () directions, respectively. In other words, the in-plane magnetoresistance appears to be straightforwardly linked to the spin-induced orbital (quadrupole) polarization of the $`e_u`$ states in the CuO4 centers. Rigorously speaking, the expression (11) is valid in the high-temperature region $`\mathrm{\Delta }(\mathrm{\Phi })kT`$, when this easily (see (6)) reduces to

$$\frac{\delta \rho _{a,b}}{\rho _{a,b}}\frac{a}{kT}\mathrm{cos}2\varphi ,$$
(12)
where, in contrast to (6), $`\varphi `$ is an in-plane azimuthal angle for the external field $`\stackrel{}{h}`$, and we take into account the near orthogonality of $`\stackrel{}{H}`$ and $`\stackrel{}{h}`$. This extremely simple expression describes all the essential features of the magnetoresistance anisotropy in the CuO2 layers and represents a main result of our model theory. Indeed, when comparing with experimental data  for $`YBa_2Cu_3O_{6+x}`$ ($`x0.3`$) one might see that the experimental angular dependence of the magnetoresistance for a rather strong external field exceeding the $`h_{flop}5T`$ obeys the ($`\varphi ,T`$) dependence (12) rather well given the relatively small albeit reasonable value of the parameter $`a`$ of the spin anisotropy: $`a+\mathrm{\hspace{0.17em}0.1}K`$.

Observed $`T`$-dependent deviations from a simple $`\mathrm{cos}2\varphi `$ law can be related to the JT effect within the $`E_u`$ doublet  which results in strong spin-vibronic effects providing the spin-induced distortions of the CuO4 center. Moreover, the JT effect results in a two-well adiabatic potential with the stabilization of the hybrid structural-orbital modes of the $`b_{1g}`$, or $`b_{2g}`$ symmetry, providing the $`\mathrm{cos}2\varphi `$, or $`\mathrm{sin}2\varphi `$ dependence of the magnetoresistance, respectively. These problems will be considered in details elsewhere.

We did not address important issues related to the origin of magnetic anisotropy and AF domain textures in doped insulating cuprates. In our opinion, the decisive role here is played by the doping-induced nucleation of the domains of a novel phase, which percolation at $`x0.4`$ results in superconductivity . The effective nucleation of these domains in the CuO2 layers is promoted by a strong in-plane charge inhomogeneity generated at least by three linearly arranged $`nn`$ chain oxygen atoms. Namely this leads to a quasi-1D stripe-like structure of the in-plane domains. The $`[100]`$ and $`[010]`$ oriented stripe-like domains induce the orthorhombic crystal structure distortions in the surrounding antiferromagnetic tetragonal matrix and generate an appropriate in-plane magnetic anisotropy. Thus, the CuO2 layers of the underdoped tetragonal insulating YBa2Cu3O6+x could be considered as antiferromagnets with the stripe-induced fluctuating in-plane magnetic anisotropy with the competition of two easy axes $`[100]`$ and $`[010]`$, respectively.

For the well developed phase separation regime with the comparable volume fractions of both phases the in-plane resistivity could be qualitatively represented as a sum of three contributions: the semiconducting one, the stripe contribution, and the contact resistivity of the interface transport. One should note that the stripe domains in the hole doped cuprates could be considered as a source of the hole carriers (donors) for the semiconducting matrix.

The observed effects such as the low-temperature magnetic hysteresis and the memory could be associated with the quenching of the stripe texture below $`T20`$ K due to the sharp slowing down of the inter-phase boundary relaxation revealed by the 63,65Cu NQR measurements . It is of interest to note, that the low-temperature ($`20÷25`$ K) spin-glass like transition has been reported for heavily underdoped antiferromagnetic Y1-yCayBa2Cu3O6 .

## 4 Conclusion

In conclusion, we propose a model where specific $`d`$-wave angular dependence of the in-plane magnetoresistance in antiferromagnetic tetragonal YBa2Cu3O6+x ($`x0.3`$) on the orientation of the external magnetic field in ($`a,b`$) plane is assigned to the effective hole transport through the low-lying excited purely oxygen doublet O2peu state, rather than the ground $`b_{1g}(d_{x^2y^2})`$ state. Such an approach implies that the ground state of the doped CuO4 plaquette in cuprates could have substantially more complicated singlet-triplet $`{}_{}{}^{1}A_{1g}^{}{}_{}{}^{1,3}E_{u}^{}`$ structure instead of the well-isolated simple $`{}_{}{}^{1}A_{1g}^{}`$ ZR singlet. External magnetic field determines the orientation of the strong exchange field for the spin-triplet $`b_{1g}e_u:{}_{}{}^{3}E_{u}^{}`$ state of the hole CuO4 center and due to the spin-orbital coupling results in a specific orbital polarization of the $`E_u`$ doublet giving rise to the spatial anisotropy of the hole transport. The observed $`d`$-wave angular dependence of the magnetoresistance on the orientation of the external magnetic field is associated with the shift and splitting of two 1D bands, and the anisotropy of the effective mass which in turn results from the spin-induced effective low-symmetry crystalline field for the $`e_u`$ hole.

## 5 Acknowledgements

The research described in this publication was made possible in part by CRDF grant No.REC-005. The authors acknowledge a partial support from the Russian Ministry of Education, grant E00-3.4-280, and Russian Foundation for Basic Researches, grant 01-02-96404. One of us (A.S.M.) would like to thank for hospitality Max-Planck Institute für Physik komplexer Systeme and Institut für Festkörper- und Werkstofforschung Dresden where part of this work was made.


text:
# Effect of resonances on the transport properties of two-dimensional disordered systems

## Abstract

We study both analytically and numerically how the electronic structure and the transport properties of a two-dimensional disordered system are modified in the presence of resonances. The energy dependence of the density of states and the localization length at different resonance energies and strengths of coupling between resonances and random states are determined. The results show, that at energy equals to the resonance energy there is an enhancement in the density of states. In contrast, the localization length remains unaffected from the presence of the resonances and is similar to the one of the standard Anderson model. Finally, we calculate the diffusion constant as a function of energy and we reveal interesting analogies with experimental results on light scattering in the presence of Mie resonances.

PACS: 71.55.Jv, 72.15.Rn

The study of waves propagating in disordered lattices, during the past years has attracted much attention and many interesting phenomena have been, more or less, well understood. However, until recently the interest of the physics community was primarily focusing on quantum, i.e. electronic waves. It was the observation of the coherent backscattering effect in classical wave systems , the analogous effect to weak localization in the electronic case, which triggered a burst of interest in further studies of scattering disordered classical wave systems. The question of localization of classical waves has attracted attention for two reasons. First, the properties of classical waves such as light waves, microwaves, and acoustic waves in random media are of fundamental interest for their own sake. Second, classical waves can serve as a model system for testing the theory of Anderson localization of electrons experimentally in a clean way, without the complication of strong inelastic scattering and other effects of electron- electron and electron-phonon interaction. Existing theories predict the localization of classical waves under certain circumstances  leading thus to the conclusion that the analogy between quantum and classical waves works reasonably well. However, there is no conclusive experimental evidence, and thus the full correspondence is not yet established beyond any doubt (a thorough discussion may be found in Refs. ).

The present paper has been motivated by an experimental work, performed several years ago  on the scattering of light by a disordered medium in the presence of Mie resonances. A low concentration of Mie resonances leads to a strong change in the transport properties of the system, reflected in a strong reduction of the transport velocity, or equivalently, of the diffusion constant $`D(E)`$, near the resonance energies. With increasing concentration of resonant scatterers, the dip in $`D(E)`$ widens and becomes less deep. Qualitatively speaking, the transport velocity is reduced because a lot of energy is temporarily stored inside the resonance or equivalently the wave spends a lot of time (dwell time) inside the Mie resonances.

To understand this surprising result, different theoretical approaches based on different considerations were developed , , . However, there are controversies as to whether this can be applied to the quantum case, i.e., electronic waves. In fact, at low concentration of Mie resonances , in the spirit of the coherent–potential approximation (CPA) , it was shown that, unlike electronic systems, the diffusion constant $`D(E)`$ for classical waves decreases sharply close to the resonance energies. This result was viewed as due to the different Ward identity caused by an energy–dependent scattering potential in the classical case. An extension of the CPA to the strong concentration limit , the so called coated CPA, shows that the effect of the Mie resonances decreases at this regime. On the other hand, in a recent work , it was shown that an extension of the Random Matrix Theory (RMT), capable of accounting for the presence of resonance scattering, is able to explain the obtained experimental result and generalize it to both quantum and classical waves, although the corresponding Ward identities are different .

In the present paper, we extend our previous analysis on the effect of resonances on electrons in a one-dimensional disordered medium , to higher dimensions. The quantities that will monopolize our interest, will be the density of states (DOS) $`\rho (E)`$ and the localization length $`l_{\mathrm{}}(E)`$. In the one-dimensional case , at low concentration and weak coupling both these quantities are affected strongly by the presence of resonances. The DOS exhibits a Lorentzian peak at the resonance energy and the localization length is drastically decreased except at the resonance energy which remains essentially the same with the one corresponding to the standard Anderson model. We will show that in higher dimensions, the density of states exhibits the same resonant enhancement while the localization length remains unaffected in contrast to the one-dimensional case. Finally, by making use of the well known results for quasi-one-dimensional systems we will compute the diffusion constant $`D(E)`$ and identify similarities with the corresponding problem of light scattering in the presence of Mie resonances.

The mathematical model we consider is the tight-binding Anderson Hamiltonian on a two-dimensional lattice,

$$H=\underset{l}{}ϵ_l|l><l|+\underset{l,l^{}}{}t_{ll^{}}|l><l^{}|$$
(1)
where $`(l)`$ denote the sites of a two-dimensional lattice $`L\times M`$. The local site energies $`ϵ_l`$ are taken independently at random within the interval $`[W/2,W/2]`$. At some randomly chosen sites $`k`$, ($`m`$ in total), the energies $`ϵ_k`$ are taken to be fixed i.e. $`ϵ_k=E_r`$. These energies correspond to the resonance energies. To simplify the problem, we assume that all resonances have the same energy $`E_r`$. The hopping matrix elements $`t_{ll^{}}`$ are restricted to nearest neighbors. When the coupling is between random states ($`N`$ in total), $`t_{ll^{}}`$ are constant taken to be the unit of energy while the coupling between a resonant state and another state $`t_{ll^{}}`$ has a different constant value $`v`$.

To evaluate the density of states $`\rho (E)`$ of this system, we start by writing the original $`(N+m)\times (N+m)`$ Hamiltonian matrix $`H`$ in the following form

$$H=\left(\begin{array}{cc}\begin{array}{c}H_0\end{array}& V\\ V^T& E_r\times I_m\end{array}\right).$$
(2)
This form can be obtained easily by separating the resonant states from random states. $`V^T`$ denotes the transpose of $`V`$ and $`H_0`$ is a $`N`$– dimensional matrix that describes the part of the Hamiltonian $`H`$ without resonances. Here $`I_m`$ is the $`m`$–dimensional unit matrix, and $`E_r`$ is the resonance energy of each of the $`m`$ resonant scatterers. The rectangular matrix $`V`$ couples the $`m`$ resonances to $`H_0`$. All the matrix elements of each column of $`V`$ are zero except four with value $`v`$.
The Green’s function corresponding to the Hamiltonian (2), can be written as

$`G`$ $`=`$ $`\left(\begin{array}{cc}\begin{array}{c}G_0\end{array}& G_V\\ G_V^T& G_r\end{array}\right)`$ (6)
$`=`$ $`\left(\begin{array}{cc}\begin{array}{c}E+i\eta H_0\end{array}& V\\ V^T& (EE_r)\times I_m+i\eta \end{array}\right)^1.`$ (10)
To determine the different matrices in $`G`$, we use the previous relation to obtain the following set of equations:

$`(EH_0+i\eta )G_VVG_r=0,`$ (11)
$`V^TG_V+[(EE_r)I+i\eta ]G_r=1.`$ (12)
Solving (12) with respect to $`G_r`$ we get the expression

$$G_r=\frac{1}{(EE_r)I+i\eta V^TG_aV}$$
(13)
where $`G_a=\frac{1}{EH_0+i\eta }`$. In the limit of weak coupling $`v1`$ and low concentration of resonances one can make the approximation $`G_aG_0`$. Then the density of states is given by

$$\rho (E)=\frac{1}{i\pi }\left(TrG_0+Tr\frac{1}{(EE_r)I+i\eta V^TG_0V}\right)$$
(14)
which, in the limit of $`v1`$ and low concentration of resonances, can be written as

$$\rho (E)=\rho _0(E)+m\frac{\mathrm{\Gamma }/2\pi }{(EE_r)^2+(\frac{\mathrm{\Gamma }}{2})^2},$$
(15)
where $`\rho _0(E)=\frac{1}{i\pi }TrG_0`$. To a good approximation, $`\rho _0(E)`$ can be identified as the density of states of the standard two–dimensional Anderson model without resonances. The effect of the resonances is given by the second term on the right hand side of Eq. (15). It exhibits a resonance enhancement near the resonance energy $`E_r`$ described by a Lorentzian peak. The width, $`\mathrm{\Gamma }2\pi \rho _0v^2`$, is given by the decay width of individual resonances. Based on physical grounds, we expect that in the limit of weak coupling $`(v1)`$ there will be a strong degeneracy due to the resonance states. This degeneracy is responsible for the resonance enhancement of the DOS at $`E=E_r`$. Gradually, as the coupling with the other states is increased, the degeneracy and consequently the peak in $`\rho (E_r)`$ decrease. These expectations are in full agreement with Eq. (15). The same result has been found in the extended RMT, Ref. , with the use of the super-symmetry formalism. Here however Eq. (15) appears natural under the frame of the tight-binding picture.

To verify the validity of the approximations used in the derivation of Eq. (15) and to go beyond the weak coupling limit, we have performed numerical simulations of the 2d–Anderson model with $`m`$ randomly distributed resonant states. Our main concern in this study is to investigate how the DOS and the localization length behave as a function of energy for various resonance concentration $`f=m/(N+m)`$ and coupling strengths $`v`$. To evaluate these quantities, we used the iterative procedure developed in . This method is particularly suitable for the calculation of various quantities for a macroscopically large system described by a microscopic Hamiltonian. The system we considered is a cylinder of width $`M`$ in $`y`$ direction and length $`L`$ in $`x`$ direction. For $`L`$ large compared to $`M`$ the sample is essentially a quasi–one–dimensional system. For the above calculations, we have used samples of length $`L=10^5`$ and width up to $`M=100`$.

We first examine the behavior of $`\rho (E)`$ as a function of the coupling constant $`v`$ and the resonance concentration $`f`$. Our numerical results for low resonance concentration and various values of the coupling $`v`$ are reported in Fig. 1. As we can see, at $`E=E_r`$ the DOS shows a resonance enhancement which is very well reproduced by the analytical formula (15). As the resonance coupling $`v`$ increases we observe a gradual decrease of the hight of the Lorentzian peak as well as an increase of its width. In the limit $`v=1`$ we recover the Anderson case. We also investigated the dependence of the DOS on the concentration $`f`$ of resonances in the disordered sample. In the weak coupling limit, we found that by increasing the concentration $`f`$ the Lorentzian peak at the DOS increases in high (see Fig. 2) in agreement with the analytical formula (15). Our numerical calculations, for small resonance coupling $`v=0.2`$, are presented in Fig. 2. for two particular values of concentration. We remind that exactly the same behavior of $`\rho (E)`$ was observed for the one-dimensional case .

Next, we calculate the localization length $`l_{\mathrm{}}(W,E)`$, needed for the determination of diffusion constant $`D(E)`$. Using the same numerical approach  we define a localization length $`l_M(W,E)`$ depending on $`M`$, $`W`$ and $`E`$. The asymptotic localization length $`l_{\mathrm{}}(W,E)`$ corresponding to the infinite system, is then defined in the limit of $`M\mathrm{}`$ i.e.

$$l_{\mathrm{}}(W,E)l_M\mathrm{}(W,E)$$
(16)
The numerical results, performed for different values of resonance concentration, show that $`l_M(E)`$ can be considered approximately independent from the concentration $`f`$ and from the resonance coupling $`v`$ (as long as the concentration is not extremely large). Our results remain practically the same for all values of $`M`$ we have considered in our calculations (up to $`M=100`$) and $`f`$. In Fig. 3 we present the behavior of the localization length $`l_M(E)`$ , for one particular value of $`f`$, vs. the corresponding one of the standard Anderson model. This behavior of the localization length is totally different from the one found in the one-dimensional case where $`l_{\mathrm{}}(E)`$ decreases drastically everywhere except close to the resonance $`E_r`$ where it is approximately equal to the corresponding one obtained for the standard Anderson model . The underline physical reason is that the electron in the two-dimensional geometry can bypass the resonances following a different path. We expect that the same behavior will appear also in three-dimensions. In contrast, in the one-dimensional case the electron has to pass through the resonance and thus to suffer resonance tunneling. Thus, the one dimensional structure acts as a filter allowing only electrons with energy $`E_r`$ to transmit.

We now turn our attention to the calculation of the diffusion constant which will be the linking observable with the experimental results . In the quasi-one dimensional multichannel case the diffusive constant $`D(E)`$ is related with $`l_{\mathrm{}}(E)`$ and $`\rho (E)`$ through the relation

$$D(E)=\frac{l_{\mathrm{}}(E)}{4\pi \rho (E)}.$$
(17)
Then, it becomes clear from the previous analysis that $`D(E)`$ will show a dip near the resonance energy $`E_r`$ (see also ), due to the resonance enhancement of the density of states $`\rho (E)`$. However, in order to compare the behavior of $`D(E)`$ with the experimental results  we have in addition, to take into account the fact that the strength of the resonance coupling $`v`$ depend on the concentration of resonances $`f`$ as $`vf^{5/6}`$ . Fig. 4 shows $`D(E)`$ vs. energy for different concentrations. In agreement with the experimental results , we find that an increase of $`f`$ causes a decrease of $`D`$ far away from the resonance while the dip near the resonance energy $`E_r`$ widens and eventually disappears for $`v1`$. This result identifies the resonant behavior of the DOS as the main explanation for the dip in the diffusion constant. Moreover, it establishes the equivalence of transport properties in the presence of resonances, between classical and quantum waves.

In conclusion, we have studied, both analytically and numerically, the effect of resonances on the transport properties of electrons in two–dimensional tight–binding Anderson model. At low concentration and weak coupling the DOS is affected strongly by exhibiting a Lorentzian peak at resonance energy. As the coupling strength increases, the width $`\mathrm{\Gamma }`$ of the Lorentzian increases and the resonance structure in DOS is smeared out. Contrary, the localization length is unaffected from the resonances. An application of the formula, $`l_{\mathrm{}}(E)=4\pi \rho (E)D(E)`$ , identifies the behavior of $`\rho (E)`$ as the explanation for the observed experimental results of the diffusion constant. This establishes the analogy between classical and quantum waves. Our results are similar to those obtained within the ERMT .

Acknowledgements We appreciate discussion with Y. Imry and C. Tejedor. B. E gratefully acknowledge the support of the European Union through the Training and Mobility of Researchers Ultrafast Network .