text:
# Observation of the scissors mode in the quasicontinuum

## I Introduction

The phenomena ’pygmy resonance’ (PY) and ’scissors mode’ (SC) have been discussed separately in the literature. In this work, we investigate for the first time quantitatively, if those two experimentally observed phenomena have the same physical origin.

The PY was introduced as a name for non-statistical features in quasicontinuous $`\gamma `$-ray spectra at energies below the giant dipole resonance. It is e.g. broadly discussed in a review of radiative strength functions as early as in 1973 . These non-statistical features observed in different nuclear mass regions and at different energies have different physical origins. A PY in rare earth nuclei at an energy around 3.5 MeV was first reported in the $`\gamma `$-ray spectrum of 170Tm following a neutron capture reaction . A systematic investigation of the PY parameters in deformed rare earth nuclei by radiative neutron capture has been carried out in , and very recently new investigations on the PY have appeared . An important application of the PY is the calculation of astrophysical (n,$`\gamma `$) reaction rates, where the existence of a soft dipole mode can greatly enhance the calculated reactions rates in e.g. the astrophysical r-process .

The SC was discovered in 1984 in inelastic electron scattering experiments on 156Gd . The bulk of information on the SC comes, however, from $`(\gamma ,\gamma ^{})`$ experiments, also called nuclear resonance fluorescence (NRF). These data are compiled in Ref. , which also gives systematics of average excitation energies of SC states and summed orbital M1 strengths around 3 MeV.

It exists no satisfactory theoretical explanation for the PY, and it is an interesting thought that the PY and the SC might emerge due to the same physical origin, i.e. if the PY can be described in terms of the observed orbital M1 strength from NRF experiments. In order to investigate this idea, we will compare the energy, the spreading width, and the resonance strength of the two phenomena.

## II Scissors mode

In NRF experiments, the energy-integrated cross-section $`I_f`$ of resonant scattered $`\gamma `$-rays on the ground state, populating a state at excitation energy $`E`$ and then decaying down to a low-lying final state $`f`$, is given by

$$I_f=\frac{\pi ^2\mathrm{}^2c^2g}{E^2}\frac{\mathrm{\Gamma }_0\mathrm{\Gamma }_f}{\mathrm{\Gamma }}.$$
(1)
Here $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_f`$ are the partial radiative decay widths for transitions from the excited state to the ground state and the final state, respectively, $`\mathrm{\Gamma }`$ is the total radiative width of the excited state and $`g`$ is a spin-factor. In NRF experiments on even-even nuclei, the excited states are usually observed to decay into two states, the ground state and the first excited state, denoted by ’1’ (in deformed nuclei the first rotational $`2^+`$ state). The experimental information from NRF experiments on even-even nuclei consists therefore generally of the energy and the values $`I_0`$ and $`I_1`$ for every observed state. Often, instead of listing $`I_1`$ values, the branching ratio $`R`$ with

$$\frac{I_1}{I_0}=\frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }_0}=R\left(\frac{E_1}{E}\right)^3$$
(2)
is given , where $`E_1`$ is the energy of the $`\gamma `$-ray transition, populating the first excited state. Summing over all final states in Eq. (1), the total energy-integrated photon-absorption cross-section $`I_{\mathrm{M1}/\mathrm{E1},\mathrm{t}}=I_0+I_1+\mathrm{}`$ is obtained. Experimentally, one observes the partial energy-integrated photon-absorption cross-sections of excited states $`I_{\mathrm{M1}/\mathrm{E1},\mathrm{p}}^{\mathrm{NRF}}=I_0+I_1`$ where the word ’partial’ refers to the fact, that only the two decay branches to the ground state and to the first excited state are taken into account in the sum. As can be seen from Eq. (1),

$$\frac{I_{\mathrm{M1}/\mathrm{E1},\mathrm{p}}^{\mathrm{NRF}}}{I_{\mathrm{M1}/\mathrm{E1},\mathrm{t}}}=\frac{\mathrm{\Gamma }_0+\mathrm{\Gamma }_1}{\mathrm{\Gamma }}.$$
(3)
The decay pattern in odd-nuclei is much more complex, since one can populate by dipole radiation excited levels with, in general, three different spins. In NRF experiments on odd gadolinium and dysprosium isotopes , branching ratios have been reported to several different low-lying states. However, very few individual excited states show $`\gamma `$ decay to more than one of these states at the same time. In addition, the observed branching ratios scatter strongly. For many excited states, only the ground state transition has been detected. A probable explanation for this observation is that the dipole strength in odd nuclei is highly fragmented, thus ground state transitions in odd nuclei are in general much weaker than in even nuclei. The even weaker branches to low-lying states might therefore easily fall below the experimental detection limit. On the other hand, due to the high fragmentation, the chances will increase that accidentally two excited states are separated in energy by the excitation energy of some low-lying state, and that the two respective peaks in the $`\gamma `$-ray spectrum are therefore interpreted as two branches from one excited state. This problem has e.g. already been observed in the even nucleus 172Yb .

In general, data on odd nuclei are very scarce (see e.g.  and references therein), and the results are not conclusive. The summed observed M1 strength in odd nuclei is by far weaker than in even-even nuclei. This has been explained in terms of unresolved strength hiding in the experimental background. A fluctuation analysis on the background has been performed , and indeed the resulting summed M1 strength in 157Gd has been found to be approximately the same value as in the neighboring even nuclei. However, it has not yet been shown that similar results can be achieved for the much less fragmented M1 strength in 161Dy, which has been investigated with a quite similar detection threshold as 157Gd , still, a new analysis is in progress and might prove the opposite . In addition, the observed irregular branching properties of excited states in odd nuclei show that the NRF technique is taken to its very limit in the investigation of odd nuclei, and that global SC observables, like summed M1 strengths and average excitation energies of SC states from such experiments, have to be taken with great care. As a last and very important point one should mention that the M1 character of the observed radiation in NRF experiments on odd nuclei has never been measured experimentally. In fact, E1 and even E2 radiation can not be excluded . In the further discussion, we will therefore mainly concentrate on even nuclei, where the deduced data are based on a more transparent experimental situation.

## III Pygmy resonance

The PY photon-absorption cross-section is usually parameterized by a Lorentzian function $`\sigma _{\mathrm{py}}(E)`$

$$\sigma _{\mathrm{py}}(E)=\sigma _{\mathrm{py}}\left(1+\frac{(E^2E_{\mathrm{py}}^2)^2}{E^2\mathrm{\Gamma }_{\mathrm{py}}^2}\right)^1,$$
(4)
where $`\sigma _{\mathrm{py}}`$, $`E_{\mathrm{py}}`$ and $`\mathrm{\Gamma }_{\mathrm{py}}`$ are the strength, centroid and width of the PY, respectively. These parameters have been obtained by fits to experimental data , which are shown in Fig. 1. Previously investigated M1 strength arises due to the tail of the spin-flip resonance (GMDR) . This resonance can be described by a Lorentzian function $`\sigma _{\mathrm{sf}}(E)`$ as well

$$\sigma _{\mathrm{sf}}(E)=\sigma _{\mathrm{sf}}\left(1+\frac{(E^2E_{\mathrm{sf}}^2)^2}{E^2\mathrm{\Gamma }_{\mathrm{sf}}^2}\right)^1,$$
(5)
where $`\sigma _{\mathrm{sf}}`$, $`E_{\mathrm{sf}}`$ and $`\mathrm{\Gamma }_{\mathrm{sf}}`$ are the strength, centroid and width of the GMDR, respectively.

Assuming that the PY is due to magnetic dipole strength, one can add the two cross-sections incoherently and obtain the total M1 photon-absorption cross-section $`\sigma _{\mathrm{M1},\mathrm{t}}(E)`$

$$\sigma _{\mathrm{M1},\mathrm{t}}(E)=\sigma _{\mathrm{py}}(E)+\sigma _{\mathrm{sf}}(E).$$
(6)
One can, at this point, discuss if deviations from the Lorentzian description of the PY and interference between the PY (presumably orbital M1 strength) and the spin-flip resonance might occur. Neither of the two effects is, however, manifested in the experimental data in Refs.  (see e.g. Fig. 1), and therefore they must be very limited in size and are neglected in the further calculation. In the upper panels of Fig. 2, the two contributions to the total M1 photon-absorption cross-section and their sums for the nuclei 162Dy and 172Yb are shown.

When comparing to the results of the NRF experiments, one should, however, in analogy with Eq. (3) define the partial M1 photon-absorption cross-section $`\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E)`$, where branching only to the ground state and the first excited state is taken into account, i.e.

$$\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E)=\sigma _{\mathrm{M1},\mathrm{t}}(E)\frac{\mathrm{\Gamma }_0+\mathrm{\Gamma }_1}{\mathrm{\Gamma }}.$$
(7)
Here, the index ’QC’ refers to the fact that the semi-experimental quantity $`\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E)`$ is calculated by means of strength functions obtained in the quasicontinuum. The total radiative width in Eq. (7) is given by

$`\mathrm{\Gamma }`$ $`=`$ $`D_{1^+}{\displaystyle \underset{\mathrm{X}L}{}}{\displaystyle \underset{I_f,\mathrm{\Pi }_f}{}}{\displaystyle _0^E}dE_\gamma `$ (9)
$`E_\gamma ^{2L+1}f_{\mathrm{X}L}(E_\gamma )\rho (EE_\gamma ,I_f,\mathrm{\Pi }_f),`$
where the second sum is going over all final levels with spin and parity $`I_f,\mathrm{\Pi }_f`$ which are accessible by multipole radiation of the type X$`L`$ from a $`1^+`$ state, $`f_{\mathrm{X}L}`$ are the respective radiative strength functions, $`\rho `$ is the level density and $`D_{1^+}`$ is the average spacing of $`1^+`$ states. The sum of the partial radiative widths $`\mathrm{\Gamma }_0+\mathrm{\Gamma }_1`$ can be calculated by taking into account in Eq. (9) only the two lowest states in the level density, i.e.

$`\rho (EE_\gamma ,I_f,\mathrm{\Pi }_f)`$ $`=`$ $`\delta (EE_\gamma )\delta _{I_f\mathrm{\Pi }_f,0^+}`$ (11)
$`+\delta (E_1E_\gamma )\delta _{I_f\mathrm{\Pi }_f,2^+},`$
and performing the integral and the two sums, yielding finally

$$\mathrm{\Gamma }_0+\mathrm{\Gamma }_1=D_{1^+}\left[f_{\mathrm{M1}}(E)E^3+f_{\mathrm{M1}}(E_1)E_1^3\right].$$
(12)
In Eq. (9), electric quadrupole and higher-order multipole contributions can usually be neglected as has been done in Eq. (12).

In order to be able to evaluate Eqs. (7), (9) and (12), the level density and the radiative strength functions $`f_{\mathrm{E1}}`$ and $`f_{\mathrm{M1}}`$ have to be known. Recently, these data have been obtained experimentally  (see Fig. 1) for the four nuclei 161,162Dy and 171,172Yb from the (3He,$`\alpha `$) reaction by using an iteration technique. The parameters of the PY have also been extracted from these experiments. The quality of the data allows us to calculate the partial photon-absorption cross-section and make comparison to NRF data.

A problem encountered here is that radiative strength functions, obtained for $`\gamma `$ transitions in the quasicontinuum, have to be used at low excitation energies where nuclear structure effects can be predominant. For instance, the observed branching ratios in NRF experiments  for even nuclei with $`\mathrm{\Gamma }_1/\mathrm{\Gamma }_00.5`$ and 2 for $`K=1`$ and 0 states, respectively, can not be reproduced accurately within this approach, where one obtains branching ratios of approximately 1 for both cases. The branching ratios are, on the other hand, well estimated by the Alaga rules, describing the nucleus in the rotational limit and thereby assuming a good $`K`$ quantum number. Unfortunately, a nuclear structure effect represented by e.g. a good $`K`$-quantum number can not be easily incorporated in the common definition of the radiative strength function. Therefore, we use here the statistical approach of Eq. (12). Since the sum in Eq. (12) is over two low-lying states only, we can not expect that all possible nuclear structure effects are averaged out by our treatment.

In Fig. 2, the total and partial M1 photon-absorption cross-sections, calculated according to Eqs. (6) and (7), and their ratio are shown for the nuclei 162Dy and 172Yb. Branching ratios further further discussed in Sect. V. The partial energy-integrated M1 photon-absorption cross-section $`I_{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}`$ which should be compared to the experimental SC value $`I_{\mathrm{M1},\mathrm{p}}^{\mathrm{NRF}}`$, can now be calculated by

$$I_{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}=dE\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E),$$
(13)
where the integral covers the appropriate energy interval of interest.

## IV Comparison

In the following, experimental data of the PY will be compared to the global observables of SC data. First, a comparison of the centroids and widths is performed. For the SC, one obtains an average excitation energy of 2.87 MeV for 162Dy . This value changes by less than 10 keV when taking into account possible additional strength where no parity assignment was possible. The NRF experiment on 172Yb was performed without polarimeter, thus, no model-independent parity assignment could be performed. If one, however, assumes that all states with $`K=1`$ are populated by M1 radiationHere, and in the following, we will make use of the hypothesis two in Ref.  in the disentanglement of overlapping elastic and inelastic peaks in the experimental NRF $`\gamma `$-ray spectrum. It seems that also the authors of Ref.  are tending to this hypothesis (see e.g. their Fig. 5). Using hypothesis one would not change the conclusion drawn in this article. the average excitation energy of the SC states is 3.09 MeV . Including states with uncertain $`K`$ assignment yields a higher value of 3.26 MeV, mainly due to the strong transition at 3.863 MeV. These values compare well with the centroids of the PY being 2.73(5) MeV and 3.48(7) MeV for 162Dy and 172Yb, respectively .

The experimentally observed spreading of SC states in 162Dy and 172Yb, which is estimated by the energy difference between the highest and the lowest observed SC state in the NRF experiment, are 0.67 MeV and 1.0 MeV, respectively . The latter value increases to 1.3 MeV when taking the transition at 3.863 MeV into account. These values are in good agreement with the widths of the PY, which were determined to be 1.35 MeV and 1.30 MeV for 162Dy and 172Yb, respectively .

Figure 3 shows the distribution of SC states for both nuclei. In the same figure, also $`\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E)`$ is displayed. The distribution of observed SC states coincides very well with the maximum and the width of the $`\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E)`$ curves.

In the second step, the strength of the PY is compared to the summed cross-section of SC states. For 162Dy, one obtains $`I_{\mathrm{M1},\mathrm{p}}^{\mathrm{NRF}}=0.41(4)\mathrm{MeV}\mathrm{mb}`$ for all transitions with M1 assignment and $`0.44(5)\mathrm{MeV}\mathrm{mb}`$, when also taking into account transitions to states with undetermined parity. For 172Yb, these quantities read $`0.25(10)\mathrm{MeV}\mathrm{mb}`$ and $`0.34(15)\mathrm{MeV}\mathrm{mb}`$, respectively. The difference of some $`0.09\mathrm{MeV}\mathrm{mb}`$ between the latter values is mostly due to the state at 3.863 MeV with uncertain $`K`$ assignment and hence uncertain parity. Here, one should recall once more that the M1 character of transitions in 172Yb is not determined model independently by a polarimeter. For the PY, one has to integrate $`\sigma _{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}(E)`$ over an appropriate energy interval. We choose, for both nuclei, a 2 MeV wide energy region between 2 MeV and 4 MeV. The integration intervals are marked by vertical lines in Fig. 3. The choice of the intervals is motivated by the fact that in NRF experiments the experimental conditions allow to observe possible SC states only in this energy region. The chosen integration region coincides also quite well with $`E_{\mathrm{py}}\pm \mathrm{\Gamma }_{\mathrm{py}}`$, and is therefore covering most of the area of the PY. We obtain for 162Dy and 172Yb values of $`I_{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}=0.46\mathrm{MeV}\mathrm{mb}`$ and $`0.42\mathrm{MeV}\mathrm{mb}`$, respectively. For convenience, all data are also collected in Table I. It is now difficult to judge the errors of these values. Certainly, the Lorentzian shapes of the PY and the GMDR are somehow idealized with respect to nature. In addition, the strength functions used for evaluating Eq. (7) are all obtained experimentally in the quasicontinuum, but they have to be applied to transitions to low-lying states in e.g. Eq. (12). Some shortcomings of the statistical approach to the region of discrete levels have already been discussed in Sect. III. Probably, a considerable systematical error is therefore introduced in $`I_{\mathrm{M1},\mathrm{p}}^{\mathrm{QC}}`$. Nevertheless, the difference of 20–40% between the calculated and the experimental values for the 172Yb nucleus and the agreement within 10% for the 162Dy nucleus allow to hope that the systematic error of the calculation is not too high and within an approximate error of 30%, we find good agreement with the SC data.

Other reasons for discrepancy between estimated and experimental SC strength may exist. Firstly, some strength of the PY could in fact be E1 radiation. Such a resonance has been suggested by P. Van Isacker et al.  in a schematic calculation. Within their model, the valence neutrons oscillate with respect to the core of neutrons and protons. However, this excitation mode has not yet been observed in NRF experiments. Another possibility is that in NRF experiments SC states are populated from the ground state where pairing correlations are important. The PY strength functions are, on the other hand, obtained in the quasicontinuum, where depairing  has been observed. It has been shown, that depairing yields in general a higher summed M1 strength and also higher excitation energies of the SC .

Unfortunately, the M1 character of the PY in rare earth nuclei could not yet be determined experimentally. In fact, the PY in iron and lead nuclei have recently been measured to be of E1 character . Since the PY in these nuclei, which have other mass numbers and deformations than rare earth nuclei, occur at higher excitation energies and consequently may have different physical origins, no strong conclusions to rare earth nuclei should be drawn from these results. An indirect proof for the M1 character of the PY in rare earth nuclei comes from (n,2$`\gamma `$) experiments on 162Dy with thermal neutrons, using the sum-coincidence method. Here, F. Bečvář et al.  have found good agreement with experiment when including a resonant component around 3 MeV in their M1 strength-function model. Since they only use schematic level density functions in their analysis, this result can, however, only serve as a strong hint of the M1 character of the PY in rare earth nuclei. Still, one can conclude that different experimental observations commonly labeled as ’pygmy resonances’ might have different physical origins.

## V Applications

### A Numerical values of $`\mathrm{\Gamma }_0`$ and $`B(\mathrm{M1})`$ from NRF experiments

In NRF experiments, often ground-state decay widths $`\mathrm{\Gamma }_0`$ or reduced transition probabilities $`B(\mathrm{M1}/\mathrm{E1})`$ are given as results. These quantities can be calculated from the total energy-integrated photon-absorption cross-section $`I_{\mathrm{M1}/\mathrm{E1},\mathrm{t}}`$ by

$$I_{\mathrm{M1}/\mathrm{E1},\mathrm{t}}=\frac{\pi ^2\mathrm{}^2c^2g}{E^2}\mathrm{\Gamma }_0$$
(14)
and

$$\mathrm{\Gamma }_0=\frac{16\pi }{9g}\left(\frac{E}{\mathrm{}c}\right)^3B(\mathrm{M1}),$$
(15)
where the latter Equation is only valid for M1 transitions. However, only the partial energy-integrated photon-absorption cross-section $`I_{\mathrm{M1}/\mathrm{E1},\mathrm{p}}^{\mathrm{NRF}}`$ is measured in NRF experiments, as discussed in Sect. II. Usually, one therefore assumes that $`I_{\mathrm{M1}/\mathrm{E1},\mathrm{t}}=I_{\mathrm{M1}/\mathrm{E1},\mathrm{p}}^{\mathrm{NRF}}`$ , i.e.

$$\mathrm{\Gamma }_0+\mathrm{\Gamma }_1=\mathrm{\Gamma }$$
(16)
over the whole energy region where excited states are observed. On the other hand, Fig. 2 shows that this assumption may become incorrect for energies above $``$3 MeV. Of course, the estimates in the lower panels of Fig. 2 are quite rough, since experimental radiative strength functions obtained in the quasicontinuum are applied to the energy region of discrete levels. Still, the large deviation from the common assumption (16) can not be neglected because we have obtained good agreement between estimated and experimental partial energy-integrated photon-absorption cross-sections, see Table I. As a consequence, too low $`B(\mathrm{M1})`$ values may have been reported in NRF works. We will therefore in the following estimate the summed $`B(\mathrm{M1})`$ value of the SC, based on the assumption that the PY is entirely due to orbital M1 strength.

Combining Eqs. (14) and (15) and summing over all individual SC states, one obtains

$$I_{\mathrm{M1},\mathrm{t}}=\frac{16\pi ^3}{9\mathrm{}c}EB(\mathrm{M1})=\frac{16\pi ^3}{9\mathrm{}c}\overline{E}B(\mathrm{M1}),$$
(17)
where $`\overline{E}`$ is the average excitation energy of the SC states. Now, we estimate $`I_{\mathrm{M1},\mathrm{t}}`$ by integrating over the Lorentzian PY model of Eq. (4) from zero to infinity, yielding

$$I_{\mathrm{M1},\mathrm{t}}=\frac{\pi }{2}\sigma _{\mathrm{py}}\mathrm{\Gamma }_{\mathrm{py}}.$$
(18)
Assuming that the average excitation energy of the SC states $`\overline{E}`$ is equal to the centroid $`E_{\mathrm{py}}`$ of the PY, we finally obtain an estimate for the summed $`B(\mathrm{M1})`$ value of SC states

$$B(\mathrm{M1})=\frac{9\mathrm{}c}{32\pi ^2}\frac{\sigma _{\mathrm{py}}\mathrm{\Gamma }_{\mathrm{py}}}{E_{\mathrm{py}}}.$$
(19)
Using experimental PY parameters from Ref. , summed $`B(\mathrm{M1})`$ strengths of $`7\pm 2\mu _\mathrm{N}^2`$ and $`6\pm 2\mu _\mathrm{N}^2`$ for 162Dy and 172Yb, respectively, are found. Here, one should stress that these values follow rather directly from experiment, since the PY parameters are obtained by a fit to normalized, purely experimental strength function data. Neither branching ratios are needed, nor the assumption of a Lorentzian shape of the PY does influence this result much. Although these data are rather high compared to recent estimates  using effective relative $`g`$-factors of $`0.7`$, theoretical values around e.g. 6–7$`\mu _\mathrm{N}^2`$ have indeed been discussed in the literature for dysprosium and ytterbium nuclei in a sum-rule approach , assuming free values of the $`g`$-factors. If our assumption for calculating the summed SC strength holds, the reason for observing summed strengths around $`3\mu _\mathrm{N}^2`$ in NRF experiments would not be the manifestation of effective $`g`$ factors, but rather partly the lack of detected branching from SC states to higher lying states than the first excited state and partly the fact that some strength lies outside the investigated energy region of 2-4 MeV.

### B E1 strength in NRF experiments on even-even nuclei

The assumption of a dominant M1 character of the PY can also be confirmed by a reproduction of the observed E1 strength in NRF experiments using the strength functions from the quasicontinuum. For this reason, we calculate the partial E1 photon-absorption cross-section $`\sigma _{\mathrm{E1},\mathrm{p}}^{\mathrm{QC}}(E)`$ according to the model of S.G. Kadmenskiĭ et al.

$$\sigma _{\mathrm{E1},\mathrm{p}}^{\mathrm{QC}}(E)=\underset{i=1,2}{}\frac{0.7\sigma _{i,\mathrm{E1}}\mathrm{\Gamma }_{i,\mathrm{E1}}^2(E^2+4\pi ^2T^2)}{EE_{i,\mathrm{E1}}(E^2E_{i,\mathrm{E1}}^2)^2}\frac{\mathrm{\Gamma }_0+\mathrm{\Gamma }_1}{\mathrm{\Gamma }}$$
(20)
which describes well the tail of the giant electric dipole resonance and reproduces the non-vanishing E1 strength at low energies, as observed by Yu.P. Popov  in (n,$`\gamma \alpha `$) reactions. In Eq. (20), $`\sigma _{i,\mathrm{E1}}`$, $`\mathrm{\Gamma }_{i,\mathrm{E1}}`$ and $`E_{i,\mathrm{E1}}`$ are the cross-sections, the widths and the centroids of the two parts of the giant electric dipole resonance, respectively, and $`T`$ is the nuclear temperature. In the following, we use a constant value for the nuclear temperature, which has been obtained by a fit to the experimental data of Fig. 1 . Similar approaches have also been very successful in describing $`\gamma `$-ray spectra of radiative neutron-capture experiments and in the calculation of isomeric cross-sections . Integrating Eq. (20) over an appropriate energy interval should now yield a value approximately equal to the sum of all partial energy-integrated E1 photon-absorption cross-sections $`I_{\mathrm{E1},\mathrm{p}}^{\mathrm{NRF}}`$ observed in NRF experiments and defined in Eq. (3).

In Fig. 4 the distribution of partial energy-integrated E1 photon-absorption cross-sections is displayed. From the NRF experiments, we obtain $`I_{\mathrm{E1},\mathrm{p}}^{\mathrm{NRF}}=0.07(1)\mathrm{MeV}\mathrm{mb}`$ and $`0.6(3)\mathrm{MeV}\mathrm{mb}`$ for 162Dy and 172Yb, respectively, when taking into account only levels with certain parity assignments in the case of 162Dy and certain $`K`$ assignments in the case of 172Yb. These values read $`0.10(2)\mathrm{MeV}\mathrm{mb}`$ and $`0.7(3)\mathrm{MeV}\mathrm{mb}`$ for the two nuclei when also taking into account levels with uncertain assignments. Integrating Eq. (20) from 2 MeV to 4 MeV, as it is done in the case of M1 radiation in Sect. IV, we obtain from the quasicontinuum $`0.5\mathrm{MeV}\mathrm{mb}`$ and $`0.7\mathrm{MeV}\mathrm{mb}`$ for 162Dy and 172Yb, respectively. Also here, we judge the errors of these quantities to be around 30% due to systematical errors, see Sect. IV. For easier comparison all values can also be found in Table I.

Clearly, the observed strength in NRF experiments on 172Yb is very well described by the strength functions from the quasicontinuum, whereas this is not at all true for 162Dy. We have no good explanation for the discrepancy in the latter nucleus. One can of course argue that the E1 strength function from the quasicontinuum might not describe well the observed E1 strength in the region of discrete levels. We find, however, that our approach works well for 172Yb. Certainly, the situation would be even worse when assuming a substantial part of the PY to be of E1 character. The good agreement in terms of the M1 strength in both nuclei and the E1 strength in 172Yb would be degraded, while the overestimation of observed E1 strength in NRF experiments on 162Dy by the strength function from the quasicontinuum would even increase. Another possibility is that a substantial part of the E1 strength in the NRF experiment is hidden in the experimental background. If we, however, accept a detection limit of some $`0.005\mathrm{MeV}\mathrm{mb}`$, also this seems quite unrealistic. At last, one can imagine some distinct nuclear structure effect responsible for hindering the population of $`1^{}`$ states at around 3 MeV from the ground state. The physical nature of this possible effect is, however, not clear at all.

### C Missing strength in NRF experiments on odd nuclei

An important issue of the SC is the question of missing strength in odd nuclei observed in NRF experiments. The solution of this problem is that in the case of odd nuclei, indeed a large portion of M1 strength is hidden in the experimental background . As a consequence, the summed M1 strength in odd nuclei might be comparable to the observed strength in neighboring even nuclei. The method for recovering the full strength yields, however, large errors and gives at the present unsatisfactory results for 161Dy . The latter statement might be disproved be a new experimental analysis . The PY parameters, on the other hand, show no odd-even effect at all  and the summed $`B(\mathrm{M1})`$ strengths in odd nuclei, calculated according to Eq. (19) are $`9\pm 2\mu _\mathrm{N}^2`$ and $`7\pm 2\mu _\mathrm{N}^2`$ for 161Dy and 171Yb, respectively, which are equal within the errors to those of neighboring even nuclei. In the following, a possible explanation of the observations in NRF experiments on odd nuclei is given by means of strength functions from the quasicontinuum.

The calculated partial M1 and E1 photon-absorption cross-sections using Eqs. (7) and (20) are displayed in Fig. 5. Also, the sum of all energy-integrated partial photon-absorption cross-sections from Ref.  is calculated, thereby taking into account all observed branching. Integrating the E1 and M1 curves of Fig. 5 as usual from 2 MeV to 4 MeV yields $`0.062\mathrm{MeV}\mathrm{mb}`$ and $`0.083\mathrm{MeV}\mathrm{mb}`$, respectively, and again a systematical error of 30% is assumed. The sum of both contributions is $`0.145\mathrm{MeV}\mathrm{mb}`$ and agrees well with the value from the NRF experiment $`0.11(2)\mathrm{MeV}\mathrm{mb}`$. Unfortunately, in NRF experiments on odd nuclei, neither parity nor $`K`$ assignments can be performed reliably . Thus one can not distinguish between M1 and E1 strength. Again, the numerical values are collected in Table I for convenience.

Figure 5 also shows the branching ratio $`\mathrm{\Gamma }_0+\mathrm{\Gamma }_1/\mathrm{\Gamma }`$ in 161Dy. This branching ratio is generally much smaller in odd nuclei than in even nuclei 162Dy due to the higher level density in odd nuclei entering Eq. (9). Here, the physical reason for the missing strength in odd nuclei becomes obvious. We find that the missing strength in NRF experiments might not only be due to unobserved elastic photon scattering below the experimental threshold, as suggested in Refs. , but also due to unobserved branching to higher lying states than the first excited states; a possibility which was not mentioned in earlier works and which can contribute much more to the summed strength than in even nuclei. These branches could, in principle, be detected in a similar type of fluctuation analysis on the experimental background as in Ref. ; still, one has to expect a general tendency of these branches to appear at lower $`\gamma `$-ray energies in the $`\gamma `$-ray spectrum.

Additional sources of uncertainty in the calculation of partial photon-absorption cross-sections in odd nuclei comes from the use of the branching ratio $`\mathrm{\Gamma }_0+\mathrm{\Gamma }_1/\mathrm{\Gamma }`$ and the fact that only 7/2 states are taken into account as intermediate statesSince the ground state spin of 161Dy is 5/2, one can also populate 3/2 and 5/2 states by absorption of dipole radiation.. As discussed in Sect. II, observed branching for odd nuclei in NRF experiments is much more complex than for even nuclei since it is dependent of the spin of the intermediate state and it does not only involve decay to the first excited state. The use of $`\mathrm{\Gamma }_0+\mathrm{\Gamma }_1/\mathrm{\Gamma }`$ calculated for intermediate 7/2 states should therefore be regarded as an approximation to the average branching properties of all observed states in odd nuclei.

## VI Conclusion

In conclusion, global parameters of the SC have been compared to PY parameters. A good agreement has been found for the centroids, the widths and the partial energy-integrated M1 photon-absorption cross-sections. Therefore, although the M1 character of the PY has not been verified experimentally, we conclude, that the PY in rare earth nuclei with high probability originates from orbital M1 strength. Still, direct measurements of the M1 character of the PY by e.g. the sum-coincidence method, using experimental level densities, is very desirable.

Further, we have shown that branching to highly excited states from SC states might be important in the estimation of B(M1) strengths. Our approach, which does not depend on the knowledge of branching ratios, gives for even and odd nuclei summed $`B(\mathrm{M1})`$ strengths of 6-9$`\mu _\mathrm{N}^2`$ for the SC in the quasicontinuum, which follow rather directly from experimental data and are in good agreement with a sum-rule approach using free $`g`$-factors. The E1 strength observed in NRF experiments could also be reproduced by the strength functions from the quasicontinuum in the case of 172Yb, whereas in the case of 162Dy a discrepancy has been observed.

Also the missing strength in NRF experiments on odd nuclei compared to even nuclei has been interpreted in terms of additional branching to higher lying states than the first excited states. Since branching might play a crucial role in the interpretation of NRF data, we would like to confront our branching estimates to experimental measurements of branching ratios from the sum-coincidence method. Here, one should especially keep in mind that the radiative strength functions used for estimating branching ratios in this work are only based on experimental data down to energies of $`1.5`$ MeV, whereas below, extrapolations of the experimental data by theoretical models have to be used.

The investigation of the SC in the quasicontinuum is complementary to the conventional NRF method. It might yield valuable data on branching ratios, summed M1 strengths and possible odd-even effects, as well as the influence of depairing on the summed M1 strength.

## ACKNOWLEDGMENTS

The authors are grateful to J. Kopecky, E. Lipparini, P. von Neumann-Cosel and J. Enders for interesting discussions. This work is supported by the Norwegian Research Council (NFR).


format: text

text:
# 1 Introdution

## 1 Introdution

Many extensions of the Standard Model (SM) contain exotic fermions. Strongly interacting exotics, such as heavy quarks, can be produced in abundance at the Tevatron or the LHC. However, particles which are not strongly interacting, such as heavy charged leptons, can best be produced at an electron-positron collider. In general, studies of heavy charged leptons at such colliders focus on s-channel production, through a $`\gamma `$, $`Z`$, $`Z^{}`$, etc. The phenomenology of exotic particles has been considered widely \[1-8\]. A good report can be found in Ref. .

In this paper, we note that a model which arises from superstring-inspired $`E_6`$ grand unification models will allow pair production of heavy charged leptons in the t-channel. We discuss this model, and study the forward-backward and left-right asymmetries at linear colliders. For simplicity, we neglect mixing between extra particles (bosons or fermions) and the normal particles of the SM, since such mixing angles are generally small.

## 2 The model

There are many phenomenologically acceptable low energy models which arise from $`E_6`$.

$`(a)E_6`$ $``$ $`SU(3)_C\times SU(2)_L\times U(1)_Y\times U(1)_Y^{}`$
$`(b)E_6`$ $``$ $`SO(10)\times U(1)_\psi SU(5)\times U(1)_\chi \times U(1)_\psi `$
$`(c)E_6`$ $``$ $`SU(3)_C\times SU(2)_L\times SU(2)_R\times U(1)_L\times U(1)_R`$
$`(c^{})E_6`$ $``$ $`SU(3)_C\times SU(2)_L\times SU(2)_I\times U(1)_Y\times U(1)_Y^{}`$
where $`U(1)_\psi `$ and $`U(1)_\chi `$ can be combined into $`U(1)_\theta `$ in model (b), reducing it to the effective rank-5 model $`SU(3)_C\times SU(2)_L\times U(1)_Y\times U(1)_\theta `$, which is most often considered in the literature. Models (c) and $`(c^{})`$ come from the subgroup $`SU(3)_C\times SU(3)_L\times SU(3)_R`$. The 27-dimensional fundamental representation has the branching rule

$$\mathrm{𝟐𝟕}=\underset{q}{\underset{}{(\mathrm{𝟑}^c,\mathrm{𝟑},\mathrm{𝟏})}}+\underset{\overline{q}}{\underset{}{(\overline{\mathrm{𝟑}}^c,\mathrm{𝟏},\overline{\mathrm{𝟑}})}}+\underset{l}{\underset{}{(\mathrm{𝟏}^c,\overline{\mathrm{𝟑}},\mathrm{𝟑})}}$$
(1)
and the particles of the first family are assigned as

$$\left(\begin{array}{c}u\\ d\\ h\end{array}\right)+\left(\begin{array}{ccc}u^c\hfill & d^c& \hfill h^c\end{array}\right)+\left(\begin{array}{ccc}E^c& \nu & N\\ N^c& e& E\\ e^c& \nu ^c& S^c\end{array}\right)$$

where $`SU(3)_L`$ operates vertically and $`SU(3)_R`$ operates horizontally. (Different symbols for these particles may be used in the literature.)

The most common method of breaking the $`SU(3)_R`$ factor is to break the 3 of $`SU(3)_R`$ into 2+1, so that $`(u^c,d^c)`$ forms an $`SU(2)_R`$ doublet with $`h^c`$ as a $`SU(2)_R`$ singlet. This gives Model (c), the familiar left-right symmetric model . Model (c) can be reduced further to an effective rank-5 model with $`U(1)_{V=L+R}`$. Another possibility, resulting in Model $`(c^{})`$, is to break the 3 of the $`SU(3)_R`$ into $`\mathrm{𝟏}+\mathrm{𝟐}`$ so that $`(d^c,h^c)`$ forms an $`SU(2)`$ doublet with $`u^c`$ as a singlet. In this option, the $`SU(2)`$ doesn’t contribute to the electromagnetic charge operator and it is called $`SU(2)_I`$ (I stands for Inert). Then the vector gauge bosons corresponding to $`SU(2)_I`$ are neutral. Model $`(c^{})`$ can be reduced to an effective rank-5 model $`SU(3)_C\times SU(2)_L\times U(1)_Y\times SU(2)_I`$. Both of them will be considered in this paper.

At the $`SU(2)_L\times SU(2)_I\times U(1)_Y\times U(1)_Y^{}`$ level, a single generation of fermions can be represented as

$$\left(\begin{array}{cc}\nu & N\\ e^{}& E^{}\end{array}\right)_L,\left(\begin{array}{c}u\\ d\end{array}\right)_L,\left(\begin{array}{cc}d^c\hfill & \hfill h^c\end{array}\right)_L,\left(\begin{array}{c}E^c\\ N^c\end{array}\right)_L,\left(\begin{array}{cc}\nu ^c\hfill & \hfill S^c\end{array}\right)_L,h_L,e_L^c,u_L^c$$

where $`SU(2)_{L(I)}`$ acts vertically (horizontally). Note that additional heavy leptons $`\left(\genfrac{}{}{0pt}{}{N}{E}\right)`$ and its conjugate $`\left(\genfrac{}{}{0pt}{}{E^c}{N^c}\right)`$ form two new isodoublets under $`SU(2)_L`$.

## 3 Cross Section Production and Asymmetries

The relevant interactions for the process $`e^+e^{}E^+E^{}`$ are

$``$ $`=`$ $`{\displaystyle \underset{f=e,E}{}}Q_f\overline{f}_\alpha \gamma ^\mu f_\alpha A_\mu +{\displaystyle \frac{g}{\mathrm{cos}\theta _W}}\overline{e}_\alpha \gamma ^\mu (T_{e_\alpha }^3Q_e\mathrm{sin}^2\theta _W)e_\alpha Z_\mu `$ (2)
$`+{\displaystyle \frac{g}{2\mathrm{cos}\theta _W}}\overline{E}_\alpha \gamma ^\mu (12\mathrm{sin}^2\theta _W)E_\alpha Z_\mu `$
$`+{\displaystyle \frac{g_I}{2\sqrt{2}}}\overline{e}\gamma ^\mu (1\gamma _5)EW_{I\mu }+H.c.`$
$`+{\displaystyle \frac{g_I}{4}}(\overline{E}\gamma ^\mu (1\gamma _5)E\overline{e}\gamma ^\mu (1\gamma _5)e)Z_{I\mu }`$
$`+{\displaystyle \underset{f=e,E}{}}g_Y^{}{\displaystyle \frac{Y_{f_\alpha }^{}}{2}}\overline{f}_\alpha \gamma ^\mu f_\alpha Z_\mu ^{}`$
where $`\alpha =L`$ or $`R`$. $`g`$, $`g_I`$ and $`g_Y^{}`$ are coupling constants and $`\theta _W`$ is the electroweak mixing angle. For simplicity, we will assume that $`g_I=g`$ and $`g_Y^{}=g_Y`$ in our numerical results, it is straightforward to relax this assumption. The first two lines are couplings between fermions and standard $`\gamma `$ and Z. The rest are couplings with extra neutral gauge bosons. The $`e^+e^{}E^+E^{}`$ process can proceed via s-channel exchange of a $`\gamma ,Z,Z^{}`$ or $`Z_I`$, and can also proceed via t-channel exchange of a $`W_I`$. Each amplitude can be written as the form of

$$C_i\overline{v}_e\gamma ^\mu (1a_i\gamma _5)u_e\overline{u}_E\gamma _\mu (1b_i\gamma _5)v_E.$$
(3)
Note that the $`W_I`$ leads to a t-channel process unsuppressed by small mixing angles. This is unique to this model. Note that if one considered production of the heavy charged leptons which form an $`SU(2)_I`$ doublet with the muon or the tau, then the processes would be identical except that the t-channel process would be absent.

The differential cross section for this process is given by

$$\frac{d\sigma }{d\mathrm{cos}\theta }=\frac{1}{8\pi s}\sqrt{\frac{1}{4}\frac{m_E^2}{s}}\{D_1(m_E^2u)^2+D_2(m_E^2t)^2+2D_3m_E^2s\}$$
(4)
s, t and u are the Mandelstam variables, and with

$`D_1`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{5}{}}}C_iC_j\{(1+a_ia_j)(1+b_ib_j)+(a_i+a_j)(b_i+b_j)\}`$
$`D_2`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{5}{}}}C_iC_j\{(1+a_ia_j)(1+b_ib_j)(a_i+a_j)(b_i+b_j)\}`$
$`D_3`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{5}{}}}C_iC_j\{(1+a_ia_j)(1b_ib_j)\}`$ (5)
where the $`C_i`$, $`a_i`$ and $`b_i`$ are given in Table 1.

Table 1 Coefficients appearing in Eq. (5)

| i | $`C_i`$ | $`a_i`$ | $`b_i`$ |
| --- | --- | --- | --- |
| 1 | $`\frac{e^2}{s}`$ | 0 | 0 |
| 2 | $`\frac{g^2(14\mathrm{sin}^2\theta _W)(12\mathrm{sin}^2\theta _W)}{8\mathrm{cos}^2\theta _W(sm_Z^2)}`$ | $`\frac{1}{14\mathrm{sin}^2\theta _W}`$ | 0 |
| 3 | $`\frac{g_I^2}{16(sm_{Z_I}^2)}`$ | 1 | 1 |
| 4 | $`\frac{9g_Y^{}^2}{144(sm_Z^{}^2)}`$ | $`\frac{1}{3}`$ | $`\frac{1}{5}`$ |
| 5 | $`\frac{g_I^2}{8(tm_{W_I}^2)}`$ | 1 | 1 |

The forward-back asymmetry is defined by

$$A_{FB}=\frac{_0^1\frac{d\sigma }{d\mathrm{cos}\theta }d\mathrm{cos}\theta _1^0\frac{d\sigma }{d\mathrm{cos}\theta }d\mathrm{cos}\theta }{_1^1\frac{d\sigma }{d\mathrm{cos}\theta }d\mathrm{cos}\theta }$$
(6)
and the left-right asymmetry is defined by

$$A_{LR}=\frac{\sigma _L\sigma _R}{\sigma _L+\sigma _R},$$
(7)
Note that the $`C_i`$, $`a_i`$ and $`b_i`$ will be somewhat different for $`\sigma _L`$ and $`\sigma _R`$ due to the insertion of the projection operator in Eq. (3). Both $`A_{FB}`$ and $`A_{LR}`$ at $`e^+e^{}`$ colliders were studied in Ref. , but only s-channel contributions were considered.

## 4 Results

The electroweak part $`SU(2)_L\times U(1)_Y`$ has been measured precisely. Let us first consider the rank 5 case. Setting $`g_Y^{}=0`$, we have two gauge boson mass parameters $`m_{W_I}`$ and $`m_{Z_I}`$. We will assume that these masses are equal and thus there is only one mass parameter remaining, which we choose to be near the experimental lower bound for direct production, $`m_{Z_I}=650`$ GeV. This is basically the same as assuming that the gauge bosons do not substantially mix with each other. The numerical results for cross section, forward-backward and left-right asymmetries are shown in Figs. 1-3. We have plotted the results for $`E^+E^{}`$ and $`M^+M^{}`$ production, where $`M`$ is the $`SU(2)_I`$ partner of the muon or tau (the only difference will be due to the t-channel process). For comparison, we also include the standard model results for both a vectorlike heavy lepton and a chiral heavy lepton. Although we have assumed that the $`Z_I`$ mass ($`E,M`$ mass) is $`650`$ GeV ($`200`$ GeV), it is easy to see how the figures will be qualitatively modified if these assumptions are relaxed.

In the rank 6 model, one has an additional mass scale and additional coupling. If we assume that the $`g_Y^{}`$ coupling is the same as $`g_Y`$, and that the mass of the $`Z^{}`$ is $`\frac{5g_Y}{3g_I}M_{Z_I}`$, then one can recalculate the cross section, forward-backward and left-right asymmetries. We find that there is not a substantial difference from the rank 5 case, except in the immediate vicinity of the $`Z^{}`$ mass.

## 5 Conclusions

How does one detect these leptons? The main decay modes depend sensitively on the masses and mixing angles. Since the $`E`$ and its isodoublet partner $`N`$ are degenerate in the limit of no mixing, one expects the $`ENW^{}`$ to be into a virtual $`W`$, leading to a three-body decay. Since the allowed three-body phase space is very small, this decay will be negligible unless the mixing with the lighter generations is extremely small. In the more natural case, in which such mixing is not very small, the two-body decays $`E\nu _eW`$ and $`EeZ`$ would dominate. A detailed analysis of the lifetimes and the decay modes can be found in Ref. . There, it was shown that the ratio of $`\mathrm{\Gamma }(EeZ)`$ to $`\mathrm{\Gamma }(E\nu _eW)`$ is given by the ratio of $`|U_{Ee}|^2`$ to $`|U_{E\nu _e}|^2`$. This is very model-dependent.

Certainly, the signature for $`EeZ`$ would be quite dramatic. Even if the $`Z`$ decays hadronically or invisibly, the monochromatic electron, plus the invariant mass of the $`Z`$ decay products, would allow for virtually background-free detection. The signature for $`E\nu _eW`$ is less dramatic, but would lead to $`W^+W^{}`$ plus missing transverse momentum. As discussed in Ref. , requiring that the $`W`$’s decay leptonically gives a signal of $`l^+l^{}`$, where $`l=(e,\mu )`$. The backgrounds, due to $`e^+e^{}\tau ^+\tau ^{},W^+W^{}`$ and $`ZZ`$, can be eliminated by calculating the invariant mass of the charged fermion pair. The signal would be striking since it would consists of a pair of $`l^+l^{}`$ with approximately the same invariant mass.

Suppose these leptons are found. One would first learn the cross section. Unless one is in the vicinity of the $`Z_I`$ resonance, the cross section in this model would be somewhat higher than the standard model. For example, at an NLC of $`\sqrt{s}=500`$ GeV and luminosity of $`6\times 10^4`$ pb<sup>-1</sup>/yr and for a heavy lepton of $`200`$ GeV, one expects approximately $`2\times 10^4`$ SM vectorlike fermion pairs produced per year, whereas one has $`3\times 10^4`$ $`E^+E^{}`$ pairs and $`5\times 10^4`$ $`M^+M^{}`$ pairs (note that the t-channel process destructively interferes). In the vicinity of the resonance, of course, the cross section can be much larger. As discussed in the previous paragraph, if the main decay is into $`\nu W`$, then a very clear signature arises if both $`W`$’s decay into $`e\nu _e`$ or $`\mu \nu _\mu `$. This will occur approximately $`5\%`$ of the time, giving a few thousand such events per year. Necessary cuts on the transverse missing energy will reduce the number of usable events, but it should still be several hundred per year, with very low background. If the main decay is into $`eZ`$ or $`\mu Z`$, then the signature is even more dramatic.

There is no forward-backward asymmetry for the pair production of SM vectorlike fermions, while the polarization asymmetry for heavy SM chiral fermions is very small. Therefore, combining $`A_{FB}`$ with $`A_{LR}`$ would make it very straightforward to distinguish $`E^+E^{}`$ and $`M^+M^{}`$ pairs from SM fermions. The behavior of the asymmetries for each of these is very different at high $`\sqrt{s}`$.

An important point is to note that the statistical uncertainty, $`\left(\frac{1A^2}{N}\right)^{1/2}`$, is very small for this model. With the approximate number of reconstructed events being between several hundred and several thousand, this gives a statistical uncertainty of between $`1`$ and $`10`$ percent. This will be even smaller in the vicinity of the resonance. From the figures, it is clear that this uncertainty is small enough that the various models can be distinguished, even off-resonance.

We thank JoAnne Hewett for a useful conversation. This work was supported by the National Science Foundation grant NSF-PHY-9900657


text:
# The Red-Sequence Cluster Survey:The Value of Ω_𝑚 and 𝜎₈

## 1. The RCS and Cosmology

###

The RCS will provide a large sample of optically selected galaxy clusters over the redshift range $`0<z<1.4`$. Clusters are selected using a technique based on locating the red-sequence of early-type galaxies in the cluster core (Gladders & Yee 2000a). Simulations show that the RCS will be complete for clusters of line-of-sight velocity dispersions of 750 km sec<sup>-1</sup> to $`z1.1`$ (Gladders & Yee 2000b), regardless of the details of the cluster galaxy properties.

###

The large and well-understood cluster sample from the RCS will be well suited to measuring the evolution of the cluster mass spectrum, $`N(M,z)`$, and hence $`\mathrm{\Omega }_m`$ and $`\sigma _8`$. Individual cluster redshifts can be estimated photometrically from the survey data to an accuracy of $`7\%`$ (Gladders & Yee 2000a). Masses can also be estimated photometrically, using the $`B_{gc}`$ richness estimator, which has been shown to be well correlated with velocity dispersion (Yee & López-Cruz 1999) and hence mass. The $`B_{gc}`$ \- mass correlation is currently calibrated to $`z0.6`$; once calibration is complete over the entire RCS redshift range the whole survey should yield estimates of $`\mathrm{\Omega }_m`$ and $`\sigma _8`$ accurate to $`5\%`$ \- more accurate and independent of current results based on relatively small X-ray selected samples (e.g., see Henry 2000 and references therein).

###

As a qualitative illustration of the cosmology indicated by the RCS, we have used the RCS selection functions (Gladders & Yee 2000a,b) to integrate out the mass dependence in model predictions of $`N(M,z)`$ and compared the RCS measure of $`N(z)`$ to predictions from two typical cosmologies. Predictions of $`N(M,z)`$ are made from the standard Press-Schechter formalism, and multiplied by the RCS selection functions (expressed in mass and redshift, assuming a cluster with a typical luminosity function shape, galaxy mixture and cluster shape, concentration and size). These selection functions tail off to zero probability at a lower mass limit which becomes progressively more massive at higher redshift, limiting the contribution to $`N(M,z)`$ from lower-mass clusters and groups. The result of these computations, as well as the actual counts from a portion of the RCS, are shown in Figure 1. Clearly, and not surprisingly, the low $`\mathrm{\Omega }_m`$ and high $`\sigma _8`$ model is preferred. Regardless of the precise details of the shapes of the modeled curves (which are dependent on knowing the selection functions), the simple abundance of candidates at high redshift in the real data is a strong argument for a low $`\mathrm{\Omega }_m`$ universe.

## References

Gladders, M.D., & Yee, H.K.C. 2000a, AJ, 120, 2148

Gladders, M.D., & Yee, H.K.C. 2000b, to be submitted to AJ

Henry, J.P. 2000, ApJ, 534, 565

Yee, H.K.C, & López-Cruz, O. 1999, AJ, 117, 1985


text:
# Band-by-band decompositions of the Born effective charges

## I Introduction

The Born effective charge  ($`Z^{}`$) is the equivalent, for crystalline insulating solid, of the atomic polar tensor introduced for isolated molecules . It is a dynamical charge related to the macroscopic polarization induced by the collective displacements of nuclei belonging to a given sublattice. In the study of the lattice dynamics of insulating crystals, it is considered as a fundamental quantity, because it governs the amplitude of the long-range Coulomb interaction between nuclei, and the splitting between longitudinal (LO) and transverse (TO) optic phonon modes.

In simple materials, like $`A^NB^{8N}`$ binary crystals , the phonon eigenvectors are imposed by symmetry. Infra-red measurements of the splitting between LO and TO modes allows an accurate estimation of $`|Z^{}|^2/ϵ_{\mathrm{}}`$ and offers therefore an unambiguous way to extract the amplitude of $`Z^{}`$ from the experiment. However, in more complex materials like ABO3 compounds, LO and TO mode eigenvectors are not necessarily equivalent. The determination of $`Z^{}`$ from the experimental data is consequently not straightforward and requires the use of some approximations . For such compounds, the development of theoretical methods giving direct access to $`Z^{}`$ acquires therefore a specific interest.

Conventionally, the Born effective charge tensor $`Z_{\kappa ,\alpha \beta }^{}`$ of nuclei belonging to the sublattice $`\kappa `$ is defined as the coefficient of proportionality relating, under the condition of zero macroscopic electric field, the change in macroscopic polarization $`𝒫_\beta `$ along the direction $`\beta `$ and the collective nuclear displacements of atoms $`\kappa `$ along direction $`\alpha `$, times the unit cell volume $`\mathrm{\Omega }_0`$ :

$$Z_{\kappa ,\alpha \beta }^{}=\mathrm{\Omega }_0\frac{𝒫_\beta }{\tau _{\kappa \alpha }}|_{=0}.$$
(1)
However, a thermodynamical equality relates the macroscopic polarization to a derivative of the electric enthalpy $`\stackrel{~}{E}`$ and another relationship connects the forces $`F_\kappa `$ on the nuclei $`\kappa `$ to a derivative of the electric enthalpy so that $`Z^{}`$ can be alternatively defined as follows:

$$Z_{\kappa ,\alpha \beta }^{}=\frac{^2\stackrel{~}{E}}{_\beta \tau _{\kappa \alpha }}=\frac{F_{\kappa ,\alpha }}{_\beta }|_{\tau _{\kappa \alpha }=0}.$$
(2)
From these relationships, Z can therefore be thought either as (i) the change of polarization induced by the collective displacements of atoms $`\kappa `$, under zero field (ii) a mixed second derivative of the electric enthalpy or (iii) the derivative of the force induced on a nucleus $`\kappa `$ by an homogeneous effective electric field $`_\beta `$, at zero atomic displacements.

The three definitions are formally equivalent. However, each of them can lead to different algorithms for the computation of $`Z^{}`$ from first-principles. Among the approaches which are the most widely used, a first powerful and systematic procedure was introduced by Baroni, Giannozzi and Testa , who suggested to determine $`Z^{}`$ from a linear response formalism grounded on a Sternheimer equation. A different algorithm, based on a variational principle, was later reported by Gonze, Allan and Teter , yielding a new alternative expression for $`Z^{}`$. Thanks to progress in the theory of the macroscopic polarization, $`Z^{}`$ is also directly accessible from finite difference of polarization . The first two algorithms were exclusively implemented within the density functional formalism (DFT) while the last one also allowed calculations of changes in polarization within different other one-electron schemes (Hartree-Fock method , model GW approximations to many-body theory , Harrison tight-binding model ) and the Hubbard tight-binding model .

Accurate predictions of the Born effective charges have been reported for a large variety of materials. In some of these studies , the decomposition of $`Z^{}`$ in individual contributions from separate groups of occupied bands appeared as a powerful tool to identify the microscopic mechanisms monitoring its amplitude. However, the physical interpretation of these decompositions was never explicitely discussed. Moreover, unlike the total $`Z^{}`$, contributions from individual groups of bands are not uniquely defined.

In the present paper, we aim at presenting the links between the theoretical frameworks used nowadays to compute $`Z^{}`$, and at deducing from this comparison the correct way to develop a band-by-band analysis. We demonstrate that the natural decompositions arising from Eq. (1) and Eq. (2) differs, although independent approches (linear response, Berry phase approach, electronic Wannier functions) to Eq. (2) are strictly equivalent provided the phase of the wavefunctions are correctly chosen.

The paper is organized as follows. In Section II, we describe the mathematical links existing between the different expressions that can be used to determine the global value of $`Z_{\kappa ,\alpha \beta }`$ within the density functional formalism without yet referring to band-by-band decompositions. We adopt the notations of Ref. . In Section III, we discuss how contributions from isolated sets of bands can be separated from each others. We identify different expressions and discuss their meaning in terms of electronic Wannier functions. In Section IV, we illustrate our results on a numerical example, emphasizing that independent decompositions yield in practice radically different values. Finally, in Section V, we conclude on the physically correct way of performing band-by-band decomposition of $`Z^{}`$.

## II Different formulations of Z\*

For practical purposes, in what follows, the Born effective charge will be decomposed into two contributions :

$$Z_{\kappa ,\alpha \beta }^{}=Z_\kappa \delta _{\alpha \beta }+Z_{\kappa ,\alpha \beta }^{el}.$$
(3)
The first term, $`Z_\kappa `$, is the charge of the nuclei (or pseudo-ion, in case of pseudopotential calculations), and can be trivially assigned. The second, $`Z_{\kappa ,\alpha \beta }^{el}`$, is the contribution due to the electrons.

### A First derivative of the polarization

A first straightforward approach for the determination of $`Z_\kappa ^{}`$ consists in computing the difference of macroscopic polarization between a reference state, and a state where the atoms belonging to the sublattice $`\kappa `$ have been displaced by a small but finite distance $`\mathrm{\Delta }\tau _{\kappa ,\alpha }`$. The electronic contribution to $`Z^{}`$ can be obtained as :

$$Z_{\kappa ,\alpha \beta }^{el}=\mathrm{\Omega }_o\underset{\mathrm{\Delta }\tau _{\kappa ,\alpha }0}{lim}\frac{\mathrm{\Delta }𝒫_\beta ^{el}}{\mathrm{\Delta }\tau _{\kappa ,\alpha }}$$
(4)
In periodic systems, the change in electronic polarization in zero field can be computed from the King-Smith and Vanderbilt formula  :

$$𝒫_\beta ^{el}=\frac{1}{(2\pi )^3}i\underset{n}{\overset{occ}{}}s_{BZ}u_{n𝐤}|\frac{}{k_\beta }|u_{n𝐤}𝑑𝐤$$
(5)
where $`s`$ is the occupation number of states in the valence bands ($`s=2`$ in spin-degenerate system) and $`u_{n𝐤}`$ is the periodic part of the Bloch functions. Taken independently, the matrix elements of the previous equation are ill-defined because the phase of the wavefunctions at a given wavevector of the Brillouin zone is arbitrary, and thus unrelated with the phases at neighbouring k points. However, the integral of the right-hand side is a well-defined quantity, which takes the form of a Berry phase of band $`n`$, as discussed by Zak .

The King-Smith and Vanderbilt definition is valid only under the constraint that the wavefunctions fulfill the periodic gauge condition. This means that the periodic part of Bloch functions must satisfy

$$u_{n𝐤}(𝐫)=e^{i𝐆.𝐫}u_{n𝐤+𝐆}(𝐫),$$
(6)
This condition does not fix unambiguously the phase of the wavefunctions at a given k-point (even not at neighbouring k-points) but it imposes a constraint for wavefunctions at distant wavevectors. It defines a topology in k-space, within which the polarization takes the convenient form of a Berry phase.

When working within one-electron schemes (DFT, Hartree-Fock, …), a second choice of phase is present at another level. For the ground-state, the Lagrange multiplier method applied to the minimization of the Hohenberg and Kohn fonctional under orthonormalization conditions on the wavefunctions , gives the following equations :

$$H_𝐤|u_{m𝐤}=\underset{n}{\overset{occ}{}}\mathrm{\Lambda }_{mn,𝐤}|u_{n𝐤}$$
(7)
This condition, associated with the minimisation of the Hohenberg and Kohn energy functional, means that the application of the Hamiltonian to a given wavefunction generates a vector which must stay within the Hilbert space defined by the set of $`u_{n𝐤}`$ wavefunctions. We observe that a unitary transformation between the wavefunctions will leave that Hilbert space invariant, and Eq. (7) will remain satisfied provided the matrix of Lagrange multiplier $`\mathrm{\Lambda }_{mn,𝐤}`$ is transformed accordingly. In order to build Kohn-Sham band structures, the unitary transform is implicitely chosen such as to guarantee

$$\mathrm{\Lambda }_{mn,𝐤}=\delta _{mn}ϵ_{m,𝐤}$$
(8)
in which case $`ϵ_{m,𝐤}`$ correspond to the eigenvalues of the Kohn-Sham Hamiltonian and the associated functions $`u_{d,m𝐤}`$ are the Kohn-Sham orbitals. This choice is called the diagonal gauge condition. All along this work, it will be emphasized by a “$`d`$” subscript.

We note that the periodic gauge condition connect wavefunctions at different k-points, while the diagonal gauge condition fixes wavefunctions at a given k-point. The choice defined by Eq. (8) is not mandatory, and the computation of the total energy, the density, or the Berry phase (Eq. (5)) will give the same value independently of the fulfillment of Eq. (8). The diagonal gauge is the natural choice for the ground-state wavefunctions while, as it will be discussed later, another choice is usually preferred for the change in wavefunctions in linear-response calculations.

Instead of approximating Eq. (1) from finite differences, it can be computed directly. The combination of Eqs. (1), (3) and (5) gives :

$`Z_{\kappa ,\alpha \beta }^{el}={\displaystyle \frac{\mathrm{\Omega }_o}{(2\pi )^3}}i{\displaystyle \underset{m}{\overset{occ}{}}}s{\displaystyle _{BZ}}[{\displaystyle \frac{u_{n𝐤}}{\tau _{\kappa ,\alpha }}}|{\displaystyle \frac{u_{n𝐤}}{k_\beta }}+u_{n𝐤}|{\displaystyle \frac{}{k_\beta }}|{\displaystyle \frac{u_{n𝐤}}{\tau _{\kappa ,\alpha }}}]𝑑𝐤`$ (9)
where the second expectation value can be worked out :

$`{\displaystyle _{BZ}}u_{n𝐤}|{\displaystyle \frac{}{k_\beta }}|{\displaystyle \frac{u_{n𝐤}}{\tau _{\kappa ,\alpha }}}𝑑𝐤={\displaystyle _{BZ}}[{\displaystyle \frac{}{k_\beta }}u_{n𝐤}|{\displaystyle \frac{u_{n𝐤}}{\tau _{\kappa ,\alpha }}}{\displaystyle \frac{u_{n𝐤}}{k_\beta }}|{\displaystyle \frac{u_{n𝐤}}{\tau _{\kappa ,\alpha }}}]𝑑𝐤.`$ (10)
In the last expression, the first term of the right-hand side is the gradient of a periodic quantity integrated over the Brillouin zone. Within any periodic gauge, its contribution will be zero. Using the time-reversal symmetry, we arrive therefore at the final expression :

$`Z_{\kappa ,\alpha \beta }^{el}=2{\displaystyle \frac{\mathrm{\Omega }_o}{(2\pi )^3}}i{\displaystyle \underset{n}{\overset{occ}{}}}s{\displaystyle _{BZ}}{\displaystyle \frac{u_{n𝐤}}{\tau _{\kappa ,\alpha }}}|{\displaystyle \frac{u_{n𝐤}}{k_\beta }}𝑑𝐤`$ (11)
The first-derivatives of the wave functions, $`u_{n𝐤}/\tau _{\kappa ,\alpha }`$ and $`u_{n𝐤}/k_\beta `$, appearing in this expression, can be computed by linear-response techniques either by solving a first-order Sternheimer equation  or by the direct minimization of a variational expression as described in Ref. .

We note that the choice of gauge will influence the value of the first-derivative of $`u_{n𝐤}`$, although the integrated quantity $`Z_{\kappa ,\alpha \beta }^{el}`$ must remain independent of this choice (in any periodic gauge). Usually, the following choice is preferred in linear-response calculations :

$`{\displaystyle \frac{u_{n𝐤}}{\lambda }}|_p|u_{m𝐤}=0`$ (12)
for $`m`$ and $`n`$ labelling occupied states, and $`\lambda `$ representing either the derivative with respect to the wavevector or to atomic displacements. As emphasized by the “$`p`$” subscript, this condition defines what is called the parallel gauge and insures that the changes in the occupied wavefunctions are orthogonal to the space of the ground-state occupied wavefunctions. This projection on the conduction bands is not reproduced within the diagonal gauge defined by the generalization of Eq. (7-8) at the first order of perturbation, as elaborated in Ref. .

### B Mixed second derivatrive of the electric enthalpy

The Born effective charge also appears as a mixed second derivative of the electric enthalpy. Therefore, as reported in Eq. (41) of Ref. , $`Z_{\kappa ,\alpha \beta }^{el}`$ can be alternatively formulated in terms of a stationary expression, involving the first-order derivative of the wavefunctions with respect to a collective displacement of atoms of the sublattice $`\kappa `$ and the first-order derivatives of the wavefunctions with respect to an electric field and to their wavevector  .

The mathematical equivalence between Eq. (11) and Eq. (41) of Ref.  is a consequence of the interchange theorem but can be directly highlighted from the stationary character of the latter. Indeed, as the error on $`Z^{}`$ is proportional to the product of the errors on the first-order change in wavefunctions, if $`u_{n𝐤}/\tau _{\kappa ,\alpha }`$ was known perfectly, a correct estimation of $`Z_{\kappa ,\alpha \beta }^{el}`$ should be obtained independently of the knowledge of $`u_{n𝐤}/_\beta `$. Putting therefore to zero $`u_{n𝐤}/_\beta `$ and the corresponding density changes in Eq. (41) of Ref. , most of the terms cancel out and we recover Eq. (11), which evaluated for the exact $`u_{n𝐤}/\tau _{\kappa ,\alpha }`$ must still correspond to a valid expression for $`Z^{}`$.

### C First derivative of the atomic force

By the same token as above, we can choose alternatively for $`u_{n𝐤}/\tau _{\kappa ,\alpha }`$ and the associated density derivative to vanish in Eq. (41) of Ref. , and we still obtain a valid expression for $`Z^{}`$ :

$`Z_{\kappa ,\alpha \beta }^{el}=`$ $`2\left[{\displaystyle \frac{\mathrm{\Omega }_0}{(2\pi )^3}}{\displaystyle _{BZ}}{\displaystyle \underset{n}{\overset{occ}{}}}su_{n𝐤}|{\displaystyle \frac{v_{\mathrm{ext},𝐤}^{}}{\tau _{\kappa \alpha }}}|{\displaystyle \frac{u_{n𝐤}}{_\beta }}d𝐤+{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{\Omega }_0}}[{\displaystyle \frac{v_{\mathrm{xc0}}}{\tau _{\kappa \alpha }}}(𝐫)][{\displaystyle \frac{n}{_\beta }}(𝐫)]^{}𝑑𝐫\right]`$ (13)
This equation corresponds to the third formulation of $`Z^{}`$ in which it appears as the first derivative of the force on the atoms $`\kappa `$ with respect to an electric field (Eq. (2)). Indeed, it is directly connected to the following expression of the force, deduced from the Hellmann-Feynman theorem :

$`F_{\kappa ,\alpha }^{el}=`$ $`{\displaystyle \frac{\mathrm{\Omega }_0}{(2\pi )^3}}{\displaystyle _{BZ}}{\displaystyle \underset{n}{\overset{occ}{}}}su_{n𝐤}|{\displaystyle \frac{v_{\mathrm{ext},𝐤}^{}}{\tau _{\kappa \alpha }}}|u_{n𝐤}d𝐤+{\displaystyle _{\mathrm{\Omega }_0}}[{\displaystyle \frac{v_{\mathrm{xc0}}}{\tau _{\kappa \alpha }}}(𝐫)][n(𝐫)]^{}𝑑𝐫`$ (14)
Compared to Eq. (11) and Eq. (41) of Ref. , Eq. (13) has the advantage that the computation of the first-order wavefunction derivative with respect to the electric field perturbation is the only computationally intensive step needed to deduce the full set of effective charges. We note however that the implementation of Eq. (13), rather easy within a plane wave – pseudopotential approach, is not so straightforward when the basis set is dependent on the atomic positions, as in LAPW methods (additional Pulay terms must be introduced).

## III Band-by-band decompositions

### A Displacement of the center of gravity of Wannier functions

Inspired by a previous discussion by Zak , Vanderbilt and King-Smith  emphasized that the macroscopic electronic polarization acquires a particular meaning when expressed in terms of localized Wannier functions. The periodic part of Bloch functions $`u_{n𝐤}(𝐫)`$ are related to the Wannier functions $`W_n(𝐫)`$ through the following transformations :

$`u_{n𝐤}(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{𝐑}{}}e^{i𝐤.(𝐫𝐑)}W_n(𝐫𝐑)`$ (15)
$`W_n(𝐫)`$ $`=`$ $`{\displaystyle \frac{\sqrt{N}\mathrm{\Omega }_o}{(2\pi )^3}}{\displaystyle _{BZ}}e^{i𝐤.𝐫}u_{n𝐤}(𝐫)𝑑𝐤`$ (16)
From this definition, we deduce that :

$`{\displaystyle \frac{}{𝐤_\beta }}u_{n𝐤}(𝐫)={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \underset{R}{}}[i(r_\beta R_\beta )]e^{i𝐤.(𝐫𝐑)}W_n(𝐫𝐑)`$ (17)
where R runs over all real space lattice vectors. Introducing this result in Eq. (5), we obtain :

$`𝒫_\beta ^{el}={\displaystyle \frac{s}{\mathrm{\Omega }_o}}{\displaystyle \underset{n}{\overset{occ}{}}}{\displaystyle r_\beta }.|W_n(𝐫)|^2d𝐫`$ (18)
From this equation, the electronic part of the polarization is simply deduced from the position of the center of gravity of the electronic charge distribution, as expressed in terms of localized Wannier functions. In other words, for the purpose of determining the polarization, “the true quantum mechanical electronic system can be considered as an effective classical system of quantized point charges, located at the centers of gravity associated with the occupied Wannier functions in each unit cell” .

We observe that Eqs. (15) and (16) establish a one-to-one correspondence between $`u_{n𝐤}`$ and $`W_n`$. As previously emphasized in Section III, when working within the diagonal gauge, $`u_{d,n𝐤}`$ becomes identified with the Kohn-Sham orbitals so that the associated $`W_{d,n}`$ will correspond to a single band Wannier function. Within this specific gauge, we can therefore isolate $`P_{m,\beta }`$, the contribution of band $`m`$ to the $`\beta `$ component of the polarization, by separating the different term in the sum appearing in Eq. (18) :

$`𝒫_{m,\beta }^{el}={\displaystyle \frac{s}{\mathrm{\Omega }_o}}{\displaystyle r_\beta }.|W_{d,m}(𝐫)|^2d𝐫`$ (19)
If we take the derivative of the polarization with respect to a collective atomic displacement, $`Z_{\kappa ,\alpha \beta }^{el}`$ can be written in terms of Wannier functions as :

$`Z_{\kappa ,\alpha \beta }^{el}={\displaystyle \underset{n}{\overset{occ}{}}}s{\displaystyle r_\beta [(\frac{W_n(𝐫)}{\tau _{\kappa ,\alpha }}|_d)^{}W_{d,n}(𝐫)+(W_{d,n}(𝐫))^{}\frac{W_n(𝐫)}{\tau _{\kappa ,\alpha }}|_d]𝑑𝐫}`$ (20)
As for the polarization, this equation has a simple physical meaning. In response to an atomic displacement, the electronic distribution is modified and the electronic contribution to $`Z^{}`$ can be identified from the displacement of the center of gravity of the occupied Wannier functions. Working within the diagonal gauge at any order of perturbation, we will be able to follow the change of single band Wannier functions all along the path of atomic displacements. In the previous expression, the contribution of band $`m`$ to $`Z_{\kappa ,\alpha \beta }^{el}`$ can be isolated :

$`[Z_{\kappa ,\alpha \beta }^{el}]_m=s{\displaystyle r_\beta [(\frac{W_m(𝐫)}{\tau _{\kappa ,\alpha }}|_d)^{}W_{d,m}(𝐫)+(W_{d,m}(𝐫))^{}\frac{W_m(𝐫)}{\tau _{\kappa ,\alpha }}|_d]𝑑𝐫}`$ (21)
This equation identifies the contribution from band $`m`$ to the Born effective charge as $`\mathrm{\Omega }_0`$ times the change of polarization corresponding to the displacement of a point charge $`s`$ on a distance equal to the displacement of the Wannier center of this band. Eq. (21) can also be estimated from finite difference by combining Eq. (4) and (19), providing an easy way to decompose $`Z_{\kappa ,\alpha \beta }^{el}`$ as soon as the Wannier functions of the system are known .

Alternatively, Eq. (21) can also easily be evaluated in reciprocal space :

$`[Z_{\kappa ,\alpha \beta }^{el}]_m=2{\displaystyle \frac{\mathrm{\Omega }_o}{(2\pi )^3}}is{\displaystyle _{BZ}}{\displaystyle \frac{u_{m𝐤}}{\tau _{\kappa ,\alpha }}}|_d|{\displaystyle \frac{u_{m𝐤}}{k_\beta }}|_d𝑑𝐤`$ (22)
As Bloch and Wannier functions are related through a band-by-band transformation, the contribution from band $`m`$ to $`Z_{\kappa ,\alpha \beta }^{}`$ in Eq. (22) keeps the same clear physical meaning as in Eq. (21) :

$`[Z_{\kappa ,\alpha \beta }^{el}]_m=\mathrm{\Omega }_0.\mathrm{\Delta }𝒫_{m,\beta }^{el}=\mathrm{\Omega }_0.s.\mathrm{\Delta }d_\beta `$ (23)
where $`\mathrm{\Delta }d_\beta `$ is the displacement in direction $`\beta `$ of the Wannier center of band $`m`$ induced by the unitary displacement of the sublattice of atoms $`\kappa `$ in direction $`\alpha `$. This decomposition is strictly equivalent to what is obtained when computing $`\mathrm{\Delta }𝒫_{m,\beta }^{el}`$ from finite differences either in real space, using Wannier functions and Eq. (19), as reported by Marzari et al., or within the Berry phase approach when separating band by band contribution to Eq. (5), the reciprocal space equivalent of Eq. (19).

In practical calculations, where each band can be thought as a combinaison of well-known orbitals, the displacement of the Wannier center is associated to the admixture of a new orbital character to the band and must be attributed to dynamical changes of orbital hybridizations. As illustrated in some recent studies , the decomposition of $`Z^{}`$ appears therefore as a powerful tool for the microscopic characterisation of the bonding in solids.

Let us emphasize again that the previous decomposition in terms of a single band is valid only if the diagonal gauge was used to define the Kohn-Sham wavefunctions, hence the “$`d`$” subscript in Eq. (19), (21) and (22). The ground-state wavefunctions are conventionally computed within the diagonal gauge. However, in most calculations, the first-derivatives of these wavefunctions are computed within the parallel gauge. Within this choice, the change in each Bloch function will be a mixing of different Kohn-Sham orbitals when the perturbation is applied so that the associated change in functions $`W_n`$, defined from Eq. (16), will correspond to the change of a multi-band Wannier function. Evaluating Eq. (21) or (22) within such a gauge, we will identify the displacement of a complex of bands rather than that of a single band. In practice, the first-order derivative of wavefunctions in the diagonal gauge $`\frac{du_{n𝐤}}{d\lambda }|_d`$ can be deduced from those in the parallel gauge $`\frac{du_{n𝐤}}{d\lambda }|_p`$ and the ground-state wavefunctions in the diagonal gauge $`u_{d,n𝐤}`$, by adding contributions from the subspace of the occupied bands :

$$\frac{du_{m𝐤}}{d\lambda }|_d=\frac{du_{m𝐤}}{d\lambda }|_p\underset{nm}{\overset{occ}{}}\frac{u_{d,n𝐤}|\frac{H}{\lambda }|u_{d,m𝐤}}{(ϵ_{n𝐤}ϵ_{m𝐤})}u_{d,n𝐤}$$
(24)
We note that this transformation (Eq. (24)) can present some problems when the denominator vanishes : this happens when the valence energies are degenerated. The problem can be partly bypassed by keeping a parallel transport gauge within the space of degenerated wavefunctions. Practically, this means that we will only be able to separate the contributions of disconnected set of bands.

### B Other band-by-band decompositions

After focusing on Eq. (11), we now investigate the possible band-by-band decompositions of Eq. (41) of Ref.  and Eq. (13). These expressions, unlike Eq. (11), are not written as simple sums of matrix elements, each related with a single band. However, individual contributions to Eq. (13) can be identified using the following decomposition of the density :

$`n(𝐫)=`$ $`{\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle _{BZ}}{\displaystyle \underset{n}{\overset{occ}{}}}su_{n𝐤}^{}(𝐫)u_{n𝐤}(𝐫)d𝐤.`$ (25)
It gives :

$`Z_{\kappa ,\alpha \beta }^{el}=`$ $`2{\displaystyle \frac{\mathrm{\Omega }_0}{(2\pi )^3}}{\displaystyle _{BZ}}{\displaystyle \underset{n}{\overset{occ}{}}}su_{n𝐤}|{\displaystyle \frac{v_{\mathrm{ext},𝐤}^{}}{\tau _{\kappa \alpha }}}+{\displaystyle \frac{v_{\mathrm{xc0}}}{\tau _{\kappa \alpha }}}|{\displaystyle \frac{u_{n𝐤}}{_\beta }}d𝐤`$ (26)
for which the following decomposition is obtained, using the diagonal gauge wavefunctions :

$`[\stackrel{~}{Z}_{\kappa ,\alpha \beta }^{el}]_m=`$ $`2{\displaystyle \frac{\mathrm{\Omega }_0}{(2\pi )^3}}{\displaystyle _{BZ}}su_{d,m𝐤}|{\displaystyle \frac{v_{\mathrm{ext},𝐤}^{}}{\tau _{\kappa \alpha }}}+{\displaystyle \frac{v_{\mathrm{xc0}}}{\tau _{\kappa \alpha }}}|{\displaystyle \frac{u_{m𝐤}}{_\beta }}|_d𝑑𝐤`$ (27)
This expression corresponds to the contribution of the electrons of band $`m`$ to the force induced on atom $`\kappa `$ by a macroscopic field $`_\beta `$. However, it is not equivalent to Eqs. (21) or (22). Indeed, for a particular band $`m`$, the difference between matrix elements present in Eq. (22) and (27) is (within a given gauge) :

$`\left[u_{m𝐤}|{\displaystyle \frac{v_{\mathrm{ext},𝐤}^{}}{\tau _{\kappa \alpha }}}+{\displaystyle \frac{v_{\mathrm{xc0}}}{\tau _{\kappa \alpha }}}|{\displaystyle \frac{u_{m𝐤}}{_\beta }}\right]\left[{\displaystyle \frac{u_{m𝐤}}{\tau _{\kappa ,\alpha }}}|i{\displaystyle \frac{u_{m𝐤}}{k_\beta }}\right]=`$ (28)
$`{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{\Omega }_0}}K_{\mathrm{xc}}(𝐫,𝐫^{})[{\displaystyle \frac{n}{\tau _{\kappa \alpha }}}(𝐫)]^{}{\displaystyle \frac{n_{m𝐤}}{_\beta }}(𝐫^{})𝑑𝐫𝑑𝐫^{}`$ (29)
$`+{\displaystyle \frac{1}{2}}{\displaystyle _{\mathrm{\Omega }_0}}K_{\mathrm{xc}}(𝐫,𝐫^{})[{\displaystyle \frac{n_{m𝐤}}{\tau _{\kappa \alpha }}}(𝐫)]^{}{\displaystyle \frac{n}{_\beta }}(𝐫^{})𝑑𝐫𝑑𝐫^{}`$ (30)
where $`n_{m𝐤}(𝐫)`$ is a short notation for $`u_{m𝐤}^{}(𝐫)u_{m𝐤}(𝐫)`$. The summation of these differences on all the bands and integration on the Brillouin zone gives zero, as expected. However, the band-by-band difference, Eq. (28), does not vanish. This demonstrates that the quantity defined from Eq. (27) is independent from that of Eq. (21) and has therefore no specific meaning in terms of Wannier functions. Unlike Eq. (13), Eq. (41) of Ref.  is not naturally convertible in a sum of independent band contributions.

## IV Numerical comparison

The previous theoretical results can now be illustrated on a numerical example. In what follows, we will consider the case of barium titanate (BaTiO3), a well-known ferroelectric material which is stable at high temperature in a cubic perovskite structure and exhibits non-trivial values of $`Z^{}`$ .

Our calculations have been performed within the density functional theory and the local density approximation . For the exchange-correlation energy, we used a polynomial parametrization  of Ceperley-Alder  homogeneous electron gas data. We adopted a planewave-pseudopotential approach. We choose highly transferable extended norm-conserving pseudopotentials as described in Ref. . The Ba 5s, Ba 5p, Ba 6s, Ti 3s, Ti 3p, Ti 3d, Ti 4p, O 2s, and O2p levels have been treated has valence states. The electronic wavefunction has been expanded in plane-waves up to a kinetic energy cutoff of 35 hartrees. Integrals over the Brillouin zone have been replaced by sums on a $`6\times 6\times 6`$ mesh of special k-points. The Born effective charges have been computed in the cubic phase at the optimized lattice parameter of 3.94 Å. They have been obtained by linear response following the scheme described in Ref. .

In Table I, we summarize the results obtained from independent formulations for the titanium charge ($`Z_{Ti}^{}`$). The decomposition of total $`Z^{}`$ is provided, according to Eq. (22) and Eq. (27) in the diagonal gauge and to Eq. (22) in the parallel gauge. We also compare our results to those reported independently by Marzari et al.  from a direct computation of the displacement of the center of gravity of the electronic Wannier functions.

As expected in this class of compounds , the total charge on the Ti atom is anomalously large ($`+7.25`$) and comparable in amplitude to the value of $`+7.16`$ reported independently using the Berry phase approach . The main anomalous contribution (deviation from the nominal value of the second column) is located in the O 2p bands. Similarly, the oxygen charge along the Ti–O bond is anomalously large and equal to $`5.71`$. Both Ti and O anomalous charge contributions are related to each other and can be assigned to dynamical changes of hybridization between O 2p and Ti 3d orbitals . This was explicitely demonstrated for a parent compound (KNbO3) by Posternak et al. .

We observe that the global charge is equivalent independently of the approach while it is not the case for partial contributions coming from different isolated sets of bands. First, the band-by-band decompositions obtained within the diagonal and parallel gauges are not similar. This means that the unitary transform performed when changing the gauge strongly mix the different bands. Second, the results deduced from Eq. (22) and (27) within the diagonal gauge are significantly different, demonstrating that the amplitude of the quantity defined in Eq. (28) is not negligible. Third, the results obtained from Eq. (22) within the diagonal gauge are comparable to those of Marzari et al. who explicitely computed the electronic Wannier functions and estimated Eq. (21) using a finite difference technique combining Eq. (4) and Eq. (19). This illustrates the physical interpretation of Eq. (22) in terms of localized Wannier functions : the contributions describe the displacement of the Wannier center of each given set of bands, induced in response to the displacement of the Ti atom.

## V Conclusions

In conclusion, the Born effective charges can be computed from first-principles using different techniques and algorithms. The global charge is a gauge invariant quantity and is obtained independently of the approach while a special care is needed to separate individual contributions from separate groups of occupied bands. When using linear response techniques, the identification of band-by-band contributions, equivalent to those obtained within the Berry phase approach, requires the use of Eq. (22), when working within the diagonal gauge. The contribution $`[Z_{\kappa ,\alpha \beta }^{el}]_m`$ is then directly related to the displacement in direction $`\beta `$ of the Wannier center of band $`m`$, when displacing the sublattice of atoms $`\kappa `$ in direction $`\alpha `$. The diagonal gauge condition is mandatory to identify single band contributions. The results obtained are conceptually and numerically different from those computed when using Eq. (27), independently of the gauge choice.

## VI Acknowledgments

We thank N. Marzari for interesting discussions. X.G. acknowledges financial support from the F.N.R.S (Belgium), from the program “Pole d’attraction interuniversitaire” P4/10 of the Belgian SSTC, and from the FNRS-FRFC project 2.4556.99.


text:
# Timescales of Disk Evolution and Planet Formation

## 1. Introduction

Planetary systems form out of circumstellar disks that are the remnants of star formation (Shu, Adams, & Lizano 1987). The young T Tauri stars with ages of $``$1 Myr, which frequently have optically-thick, actively-accreting disks (Strom, Edwards & Skrutskie 1993), are thought to represent a stage prior to the main epoch of planet formation. In contrast, the dust disks of the so-called Vega-like objects, main sequence stars whose ages could be as large as 1 Gyr, appear to be more evolved; most of the dust has presumably coagulated into planets or planetesimals, and the remaining dust in the optically-thin disk is continually replenished by collisions between larger bodies (Lagrange, Backman & Artymowicz 2000). However, the age at which this transition occurs is not well constrained. One reason is a paucity of intermediate-age ($`t`$10 Myr) nearby stars in current samples. Unfortunately, the age estimates of early-type isolated main sequence stars are highly uncertain.

The two competing theories for making gas giant planets in disks predict very different timescales. Core accretion mechanism, in which a solid core builds up from icy and rocky planetesimals and later acquires an envelope of nebular gas, is expected to take $``$10 Myr (e.g., Pollack et al. 1996). On the other hand, the disk instability mechanism, where a gravitationally unstable disk fragments directly into self-gravitating protoplanetary clumps, is a much faster process, forming giant planets in $``$0.1 Myr (e.g., Boss 2000). Therefore, better observational constraints on disk evolution timescales may allow us to distinguish between the two scenarios.

## 2. Insights from Nearby Stellar Groups

The ROSAT all-sky X-ray survey, combined with follow-up spectroscopy and proper motion data, has led to the recent identification of several groups of young stars within 100 pc of the sun. The proximity and age differences of these groups make them suitable laboratories for detailed studies of disk evolution and planet formation.

The TW Hydrae Association (TWA), at a distance of $``$55 pc and not associated with any obvious parent cloud, is of particular interest because of its $``$10-Myr age, as estimated using a variety of techniques. Our mid-infrared observations, carried out using the OSCIR instrument on Keck and the Cerro Tololo 4-meter telescope over the past two years, show that many of the TWA stars have little or no disk emission at 10$`\mu `$m (Jayawardhana et al. 1999b). Even among the five stellar systems with 10$`\mu `$m excesses, most show some evidence of inner disk evolution. The disk around the A0V star HR 4796A has an $`r50`$ AU central hole in mid-infrared images (Jayawardhana et al. 1998; Koerner et al. 1998). The SEDs of HD 98800 and Hen 3-600A also suggest possible inner disk holes (Jayawardhana et al. 1999a). The modest excess we detected from CD -337795 could well be due to a faint companion. Only TW Hya itself appears to harbor an optically thick, actively accreting disk of the kind observed in $``$1-Myr-old classical T Tauri stars; it is the only one with a large H$`\alpha `$ equivalent width (-220 Å).

If the TWA stars are indeed 10 million years old, our results suggest that their inner disks have already depleted either through coagulation of dust or accretion on to the central star. The fact that only one (TW Hya) out of a sample of 16 TWA members shows classical T Tauri characteristics –compared to $``$50% –90% of $``$1-Myr-old stars in star-forming regions– argues for relatively rapid evolution of inner disks in pre-main-sequence (PMS) stars. Observations at far-infrared and sub-millimeter wavelengths may reveal whether most TWA stars still retain their outer disks.

We have also obtained mid-infrared photometry and high-resolution optical spectra for eight confirmed members of the MBM 12 high-latitude cloud at $``$65 pc. These stars are estimated to be 1-3 Myr in age. Interestingly, in contrast to the TWA sample, the majority of MBM 12 stars do show signatures of optically thick, actively accreting disks (Jayawardhana et al. 2000).

## 3. Comparison to Other Studies

Our results for TWA and MBM 12 are consistent with those of Skrutskie et al. (1990) who studied infrared excess in a sample of “solar-type” PMS stars in Taurus-Auriga. Of stars estimated to be younger than 3 Myr, Skrutskie et al. found that roughly half have optically thick disks extending close to the stellar surface whereas fewer than 10% of somewhat older stars showed evidence for similar disks. Therefore, the authors suggested that few disks survive beyond 10 Myr.

More recently, Brandner et al. (2000) have compared disk properties between Chamaeleon and Sco-Cen star-forming regions to derive disk evolution timescales. They find that the 5-15 Myr-old stars in Sco-Cen have a 15$`\mu `$m/7$`\mu `$m spectral index intermediate between that of younger (1-5 Myr-old) Chamaeleon objects and that of pure stellar photospheres. Brandner et al. suggest that the difference is due to a deficiency of small dust grains relative to large dust grains in the Sco-Cen disks. If it is an evolutionary effect, dust depletion of disks must start at an age between 5 and 15 Myr, again in agreement with our findings for the TWA and MBM 12 groups. Meteoritic evidence also suggests that planetary material was made over a period of $``$ 10 Myr in the solar nebula (Podosek & Cassen 1994; Wadhwa & Russell 2000).

On the other hand, in high-mass star-forming regions such as Orion, disks may dissipate much faster (in $`<`$ 1 Myr) due to irradiation from neighboring luminous stars (e.g., Bally et al. 1998). It is not clear whether planets can form in such rapidly photoevaporating disks.

## 4. Future Prospects

It is important to study a sufficiently large sample of stars with a range of ages and environments to place strong constraints on disk evolution and planet formation timescales. The identification of additional groups of young stars in the solar neighborhood will help, as will the unprecedented sensitivity of SIRTF and SOFIA to detect even optically thin disk emission from clusters at much greater distances than hitherto possible. There may not be a universal timescale for disk evolution, especially when the effects of binary companions and external irradiation are taken into account. More studies are also needed to better understand the dissipation of the gaseous component of disks and to determine whether it persists long enough to allow for slow accumulation of gas onto giant protoplanets, as predicted by the core accretion scenario.

## References

Bally, J., Testi, L., Sargent, A., & Carlstrom, J. 1998, AJ, 116, 854

Boss, A.P. 2000, ApJ, 536, L101

Brandner, W., et al. 2000, AJ, 120, 950

Jayawardhana, R., et al. 1998, ApJ, 503, L79

Jayawardhana, R., et al. 1999a, ApJ, 520, L41

Jayawardhana, R., et al. 1999b, ApJ, 521, L129

Jayawardhana, R., et al. 2000, in preparation

Koerner, D., et al. 1998, ApJ, 503, L83

Lagrange, A.-M., Backman, D.E., & Artymowicz, P. 2000, in Protostars and Planets IV, ed. V. Mannings, A.P. Boss, & S.S. Russell (Tucson: University of Arizona Press), 639

Podosek, F.A., & Cassen, P. 1994 , Meteoritics, 29, 6

Pollack, J. B., et al. 1996, Icarus, 124, 62

Shu, F.H., Adams, F.C., & Lizano, S. 1987, ARA&A, 25, 23

Skrutskie, M.F., et al. 1990, AJ, 99, 1187

Strom, S.E., Edwards, S., & Skrutskie, M.F. 1993, in Protostars and Planets III, ed. E.H. Levy & J.I. Lunine (Tucson: Univ. Arizona Press), 837

Wadhwa, M., & Russell, S.S. 2000 in Protostars and Planets IV, ed. V. Mannings, A.P. Boss, & S.S. Russell (Tucson: University of Arizona Press), 995