text:
# The Galactic Center as a Dark Matter Gamma-Ray Source

## 1 Introduction

Experimental cosmology has been steadily progressing over the latest years. The emerging picture has been recently reinforced by the data from the Wilkinson Microwave Anisotropy Probe (WMAP) which have pinned down several fundamental parameters to a remarkable level of precision. In particular, in the latest global fit, the contribution to the critical density of non-relativistic matter has been found in the range $`\mathrm{\Omega }_mh^2=0.135_{0.008}^{+0.009}`$ (here $`h`$ is the Hubble constant in units of 100 km s-1 Mpc-1; $`h=0.71_{0.03}^{+0.04}`$), much larger than the baryonic term, $`\mathrm{\Omega }_bh^2=0.0224\pm 0.0009`$.

Unveiling the nature of non-baryonic cold dark matter (CDM) is one of the major challenges in science today. Weakly interacting massive particles (WIMPs) are among the leading dark matter candidates. They would naturally appear as another of the thermal leftovers from the early Universe, and, at the same time, their existence is predicted in several classes of extensions of the Standard Model of particle physics. The most popular of such candidates is the lightest neutralino in R-parity conserving supersymmetric models. Considerable effort has been put in the search for dark matter WIMPs in the last decade, with several complementary techniques applied (for a recent review, see, e.g., ). Among them, indirect detection through the identification of the yields of WIMP pair annihilations in dark matter structures seems to be a very promising method. In particular, we will focus here, as a signature to identify dark matter, on the possible distortion of the spectrum of the diffuse $`\gamma `$-ray flux in the Galaxy due to a WIMP-induced component, extending up to an energy equal to the WIMP mass (a list of other recent analysis on this topic includes Refs. ).

The EGRET telescope on board of the Compton Gamma-Ray Observatory has mapped the $`\gamma `$-ray sky up to an energy of about 20 GeV over a period of 5 years. The data collected by EGRET toward the Galactic center (GC) region show  high statistical evidence for a gamma-ray source, possibly diffuse rather than point-like, located within $`1.5^{}`$ of the GC ($`l=b=0^{}`$). The detected flux largely exceeds the diffuse $`\gamma `$-ray component expected in the GC direction with a standard modeling of the interaction of primary cosmic rays with the interstellar medium (see, e.g., ); the latter fails also to reproduce the spectral shape of the GC source. Although other plausible explanations have been formulated, it is very intriguing that the EGRET GeV excess shows, as basic features, the kind of distortion of the diffuse $`\gamma `$-ray spectrum one would expect from a WIMP-induced component, assuming that the dark matter halo profile is peaked toward the GC. We will identify for which classes of WIMP compositions and masses, fair agreement with the measured flux can be obtained.

No firm conclusion about the nature of the GC excess can be driven from data available at present; on the other hand, the picture is going to become much clearer in the near future. The Gamma-ray Large Area Space Telescope (GLAST) has been selected by NASA as the next major $`\gamma `$-ray mission,111A list of people and institutions involved in the collaboration together with the on-line status of the project is available at http://www-glast.stanford.edu. For a detailed description of the apparatus see ; a discussion of the main scientific goals can be found in . and is scheduled for launch in the first half of 2006. Compared to EGRET, GLAST will have a much larger effective area, better energy and angular resolutions, as well as it will cover a much wider energy range. GLAST will perform an all-sky survey of $`\gamma `$-ray sources, with scientific objectives including the study of blazars, $`\gamma `$-ray bursts, supernova remnants, pulsars, the diffuse radiation in the Galaxy, and unidentified high-energy sources. The identification of dark matter sources has been indicated as one of its main scientific goals. We illustrate here the conditions under which it may be feasible that GLAST will single out a dark matter source located at the GC.

The paper is organized as follows. In Section 2 we show that the GC EGRET $`\gamma `$-ray excess can be modeled with a component due to WIMP annihilations. In Section 3 we discuss what kind of information GLAST could provide to confirm such hypothesis. In Section 4 we illustrate the perspectives of dark matter detection with GLAST in the case in which the EGRET excess will be found to be due to another kind of source. In the first part of the paper we discuss features for the dark matter signal in the contest of a generic toy-model; in Section 5 we will focus on neutralino dark matter candidates in the minimal supergravity framework, applying to this specific case the tools developed in the other Sections. Conclusions are in Section 6.

## 2 A Dark Matter Source at the Galactic Center?

In Table 1 we report the flux per energy bin for the GC gamma-ray source as measured by EGRET, together with the expected flux from cosmic ray interactions in a standard scenario  (see also Table 2 and Fig. 4 in ). As already mentioned there is a significant mismatch between the two: we entertain here the possibility that the bulk of the high energy flux is due to pair annihilations of dark matter WIMP’s in the GC region. Hence, we assume that the total flux measured by EGRET can be described as the superposition of two contributions:

* i) a component due the interaction of primary cosmic rays with the interstellar medium, with spectral shape defined by the function $`S_b(E_\gamma )`$ (background contribution)
* ii) a component due to WIMP annihilation in the dark matter halo, whose energy spectrum is defined by $`S_\chi (E_\gamma )`$ (signal contribution).

We write the total $`\gamma `$-ray flux as:

$$\varphi _\gamma =\varphi _b+\varphi _\chi =N_bS_b+N_\chi S_\chi ,$$
(1)
where $`N_b`$ and $`N_\chi `$ are dimensionless normalization parameters, respectively, for the standard and exotic flux, which we define below.

### 2.1 The background component

There are three mechanisms which give rise to diffuse $`\gamma `$-ray radiation in the Galaxy: production and decay of $`\pi ^0`$s, inverse Compton scattering and bremsstrahlung (see, e.g., ). According to standard scenarios, in the energy range we will consider, $`E_\gamma >1\mathrm{Gev}`$, the dominant background source is given by $`\pi ^0`$ decays. Photons are generated in the interaction of primary cosmic rays with the interstellar medium via:

$$p+X\mathrm{}\pi ^02\gamma $$

$$He+X\mathrm{}\pi ^02\gamma ,$$

where $`X`$ stands for an interstellar atom, mainly $`H`$ and $`He`$. We have simulated the induced $`\gamma `$ray yield according to standard treatments (see, e.g.,) and as implemented in the Galprop software package . We assume that the $`p`$ and $`He`$ cosmic ray fluxes in the Galaxy have the same energy spectra and relative normalization as those measured in the local neighborhood , and that the $`He`$ fraction of the interstellar gas is 7%. We write the background flux, splitting it into two factors:

$$S_b(E_\gamma )=\frac{1}{\left(1\mathrm{cm}^2\mathrm{sr}\right)}Em(E_\gamma )$$
(2)
and

$$N_b=\frac{1}{\left(1\mathrm{cm}^2\mathrm{sr}^1\right)}_{l.o.s.}𝑑l\frac{n_H(l)}{4\pi }\frac{\varphi _p^{prim}(l)}{\varphi _p^{prim}(l=0)}.$$
(3)
Here $`Em(E_\gamma )`$ \[GeV-1 s-1\] is the local emissivity per hydrogen atom, i.e. the number of secondary photons with energy in the range ($`E_\gamma `$, $`E_\gamma +dE_\gamma `$) emitted per unit time per target hydrogen atom, for an incident flux of protons and helium nuclei equal to the locally measured primary proton and helium fluxes. The factor $`N_b`$ is instead associated to the interstellar hydrogen column density $`n_H(l)`$, integrated along the line of sight and weighted over the proton primary flux at the location $`l`$, $`\varphi _p^{prim}(l)`$, normalized to the local value $`\varphi _p^{prim}(l=0)`$.

Above an energy of about 1 GeV the background spectrum (and therefore the function $`\varphi _b`$) can be described by a power law of the form $`\varphi _{back}E_\gamma ^\alpha `$ with the same spectral index as the dominant primary component, i.e. the proton spectral index $`\alpha =2.72`$.

The relative normalization of the primary components in different places in the Galaxy can be estimated once a radial distribution of primary sources is chosen (following, for instance, the radial distribution of supernovae) and then by propagating the injected fluxes with an appropriate transport equation (this is what is done, e.g., in the Galprop code ). On the other hand, the hydrogen column density toward the Galactic center is very uncertain; we choose therefore to define the spectral shape of the background through the function $`S_b`$, and to keep $`N_b`$ as a free normalization parameter.

### 2.2 The signal component

The production of $`\gamma `$-rays in a dark matter halos made of WIMP’s follows essentially by the definition of WIMP, regardless of any specific scenario one has in mind. The signal scales linearly with the pair annihilation rate in the limit of non-relativistic particles.

We consider a generic framework in which the dark matter in the Galactic halo is made of particles $`\chi `$, WIMP dark matter candidates with mass $`M_\chi `$ and total pair annihilation rate into lighter Standard Model particles $`\sigma v`$ (in the limit of zero relative velocity). Among the kinematically-allowed tree-level final states, the leading channels are often $`b\overline{b},t\overline{t},\tau ^+\tau ^{},W^+W^{},Z^0Z^0`$. This is the case, e.g., for neutralinos and, more generically, for any Majorana fermion WIMP, as for such particles the S-wave annihilation rate into the light fermion species is suppressed by the factor $`m_f^2/M_\chi ^2`$, where $`m_f`$ is the mass of the fermion in the final state. The fragmentation and/or the decay of the tree-level annihilation states gives rise to photons. Again the dominant intermediate step is the generation of neutral pions and their decay into $`2\gamma `$. The simulation of the photon yield is standard; we take advantage of a simulation performed with the Lund Monte Carlo program Pythia 6.202  implemented in the DarkSUSY package .

Suppose that the dark matter halo is roughly spherical and consider the induced $`\gamma `$-ray flux in the direction that forms an angle $`\psi `$ with the direction of the Galactic center; the WIMP induced photon flux is the sum of the contributions along the line of sight (l.o.s):

$$\varphi _\chi (E,\psi )=\frac{\sigma v}{4\pi }\underset{f}{}\frac{dN_f}{dE}B_f_{l.o.s}𝑑l(\psi )\frac{1}{2}\frac{\rho (l)^2}{M_\chi ^2}$$
(4)
where $`B_f`$ is the branching ratio into the tree-level annihilation final state $`f`$, while $`dN_f/dE`$ is the relative differential photon yield. The WIMP mass density along the line of sight, $`\rho (l)`$, enters critically in the prediction for the flux, as the number of WIMP pairs is equal to $`1/2\rho (l)^2/M_\chi ^2`$. It is then useful to factorize the flux in Eq.(4) into two pieces, one depending only on the the particle physics setup, i.e. on the cross section, the branching ratios and the WIMP mass, and the other depending on the WIMP distribution in the galactic halo. We rewrite Eq.(4) as :

$`\varphi _\chi (E,\psi )`$ $`=`$ $`3.7410^{10}\left({\displaystyle \frac{\sigma v}{10^{26}\mathrm{cm}^3\mathrm{s}^1}}\right)\left({\displaystyle \frac{50\mathrm{GeV}}{M_\chi }}\right)^2{\displaystyle \underset{f}{}}{\displaystyle \frac{dN_f}{dE}}B_f`$ (5)
$`J(\psi )\mathrm{cm}^2\mathrm{s}^1\mathrm{GeV}^1\mathrm{sr}^1`$
where we introduced the dimensionless function $`J`$, containing the dependence on the halo density profile,

$$J(\psi )=\frac{1}{8.5\mathrm{kpc}}\left(\frac{1}{0.3\mathrm{GeVcm}^3}\right)\rho ^2(l)𝑑l(\psi )$$
(6)
More precisely, given a detector with angular acceptance $`\mathrm{\Delta }\mathrm{\Omega }`$, we have to consider the average of $`J(\psi )`$ over the solid angle $`\mathrm{\Delta }\mathrm{\Omega }`$ around the direction $`\psi `$:

$$J(\psi )_{\mathrm{\Delta }\mathrm{\Omega }}=\frac{1}{\mathrm{\Delta }\mathrm{\Omega }}J(\psi )𝑑\mathrm{\Omega }$$
(7)
To compare with the GC EGRET gamma-ray source, we will consider $`\mathrm{\Delta }\mathrm{\Omega }10^3`$ sr, i.e. the same magnitude as the angular region probed by the EGRET telescope.

As for the background component, we split the signal into a term which fixes the spectral shape of the flux, plus a normalization factor. In the notation introduced in Eq.(1), we label $`N_\chi J(\psi )_{\mathrm{\Delta }\mathrm{\Omega }}`$ and define $`S_\chi \varphi _\chi /N_\chi `$. The dependence on $`\rho (l)`$ has been included in the term we treat as a free normalization parameter $`N_\chi `$, as $`\rho (l)`$ is very uncertain both from the theoretical and the observational points of view. Although there is a some spread in the predictions for the $`\gamma `$-ray flux when coming to specific WIMP models, its spectral features are rather generic. As most photons are produced in the hadronization and decay of $`\pi ^0`$s, the shape of the photon spectrum is always peaked at half the mass of the pion, about 70 MeV, and it is symmetric around it on a logarithmic scale (sometimes this feature is called the “$`\pi ^0`$ bump”, see, e.g., ). The same is true for the background, but still it may be possible to tell signal from background: the signal arises in processes which have all the same energy scale, i.e. $`2M_\chi `$, therefore the WIMP induced flux, contrary to the background, is spectral index free and shows a sharp cutoff when $`E_\gamma `$ approaches the WIMP mass. This is shown in the left panel of Fig. 1, where we plot the differential photon yield per annihilation times the inverse of WIMP mass squared, for a few values of the WIMP mass, and assuming that the $`\chi `$’s have a single dominant annihilation channel ($`b\overline{b}`$ in the case displayed). In the same figure, for comparison, the spectral shape of the background is shown: as it can be seen, one may hope to identify the WIMP induced component as a distortion of the background spectrum at relatively high energies. For a given WIMP mass, the photon yields in the different annihilation channels are analogous, as shown in the right panel of Fig. 1: solid curves indicate the total photon yield, while dashed curves indicate the photon yield in radiative processes, i.e. in all processes rather than $`\pi ^0`$ decays. The spectrum for the $`t\overline{t}`$ and $`W^+W^{}`$ channels are very close to the one for $`b\overline{b}`$ (differences are mainly given by prompt decays before hadronization); only in the $`\tau ^+\tau ^{}`$ case, radiative photon emission is dominant, still with a large bump due to the hadronic decay modes of $`\tau `$ leptons.

### 2.3 EGRET data fit

EGRET has performed measures in the energy range $`30\mathrm{MeV}÷10\mathrm{GeV}`$ with few bins in the high energy region. Given the paucity of the data in the highest end of the energy region in a first approximation it is not sensible to discriminate different annihilation channels leading to photons in the final state.

It is then convenient to keep the discussion as general as possible and consider a simplified context (a *toy-model*) with a thermal relic $`\chi `$ for which only one intermediate annihilation channel is open ($`B_f=1`$ in that channel). Furthermore we assume the total annihilation cross section to scale with the inverse of the relic abundance $`\mathrm{\Omega }_\chi `$ :

$$\sigma v\sigma v\frac{310^{27}\mathrm{cm}^3\mathrm{s}^1}{\mathrm{\Omega }_\chi h^2}310^{26}\mathrm{cm}^3\mathrm{s}^1,$$
(8)
where $`\sigma v`$ is the thermally averaged annihilation cross section.

Eq.(8) is not valid in presence of resonances or thresholds near the kinematically released energy in the annihilation $`2M_\chi `$ and of coannihilation effects. In the presence of these conditions we can have large deviations from this approximate scaling.

Moreover there are cases in which the inverse proportionality between $`\mathrm{\Omega }_\chi `$ and the annihilation rate only gives a lower bound for the latter since a non-thermal relic component can provide the relic density needed to account for $`\mathrm{\Omega }_{CDM}`$.

In our toy-model the WIMP mass is kept as a free parameter , as we have shown that the photon spectrum is rather sensitive to it. The results we discuss below are then dependent on a mass scale and an overall normalization parameter and can be easily rescaled to fit an explicit model for which $`M_\chi `$ and $`\sigma v`$ are specified.

We then fit the EGRET data on the GC excess taken from Table 1 using our toy-model. We are not using the two lowest energy bins as for $`E1`$ GeV the background is most probably dominated by inverse Compton and bremsstrahlung rather than by $`\pi ^0`$ production as we assumed. We repeat this fit for different values of the WIMP mass $`M_\chi `$ and for a few tree-level annihilation channels. $`N_b`$ and $`N_\chi `$ in Eq.(1) are treated as free parameters which are varied in order to minimize the statistical $`\chi ^2`$ variable. The allowed range of variation for the background normalization, $`N_b`$, is between $`3.210^{20}`$ and $`1.810^{21}`$. The two extrema of this variation interval are taken to correspond to a best fit of the background in agreement with the standard scenario (column 2 in Table 1) and to the best fit for the EGRET data from the GC (column 3 in Table 1) with $`N_\chi =0`$.

In Fig.2 we show our best fits for two sample values of $`M_\chi `$ and two intermediate channels, $`b\overline{b}`$ and $`W^+W^{}`$. On the qualitative side, the agreement with the data is rather good, even if the reduced $`\chi ^2`$ (for 6 degrees of freedom) in the two examples displayed is still rather large (of the order of 5). This may depend on an underestimate of error bars or also on the fact that we are neglecting uncertainties in the theoretical predictions for the spectral shapes. It is clear, on the other hand, that adding a component due to WIMP annihilations on top of the background component greatly improves the agreement between the expected spectral form and the one found in the EGRET measurement. This is shown in Fig. 3 where the reduced $`\chi ^2`$ for our best fits is shown as a function of WIMP mass. These values should be compared with those obtained when only the background contribution is included ($`N_\chi =0`$). This case is represented by the horizontal line in the figure and gives a reduced $`\chi ^2150`$. Fig. 3 also indicates that light WIMP masses are marginally favored over heavier WIMP’s. Results for other tree-level annihilation states are analogous and show the same trends.

The best fit curves displayed correspond to rather large values of the normalization parameter $`N_\chi `$, a few times $`10^4`$ for the two cases shown in Fig. 3. $`N_\chi `$ tends to increase going to heavier WIMP’s as a function inversely proportional to the WIMP mass raised to a certain power. This power is a little smaller than the value 2 one would have inferred from the fact that the WIMP pair density decreases as $`M_\chi ^2`$. An investigation of the dependence of the statistical variable $`\chi ^2`$ from the WIMP mass and $`N_\chi `$ is shown in Fig. 4, for the $`b\overline{b}`$ annihilation channel. Given that $`\chi ^2=\chi ^2(N_b,N_\chi ,M_\chi )`$, we minimize it with respect to the parameter $`N_b`$ for fixed values of $`N_\chi `$ and $`M_\chi `$. The isolevel curves for the reduced $`\chi ^2`$ are plotted in Fig. 4. Fixing a value for the WIMP mass and going from the left-hand side to the right-hand side of the figure, we move from the case in which the WIMP signal is marginal with respect of the background flux to the case where the sum of the two reproduces most closely the data. Further increasing the value of $`N_\chi `$ we reach a region where the WIMP signal exceeds the flux detected by EGRET.

We recall that in our toy-model $`N_\chi `$ is identified with the halo model dependent function $`J(\psi )_{\mathrm{\Delta }\mathrm{\Omega }}`$ for the Galactic center direction $`\psi =0^{}`$ and EGRET angular acceptance $`\mathrm{\Delta }\mathrm{\Omega }=10^3`$ sr. As already mentioned, the distribution of dark matter in the inner part of the Galaxy is still a controversial issue. Dynamical measurements show that dark matter is, in the Galactic center region, just a subdominant component with respect to baryonic matter, but lack the resolution we need for an estimate of $`J`$. On the other hand, N-body simulation of structure formation in a CDM Universe, find dark matter density profiles which are singular toward the center of the Galaxy, possibly scaling as $`1/r`$ (profile of Navarro, Frenk & White  (NFW)) or $`1/r^{1.5}`$ (Moore et al. profile ) as the galactocentric distance $`r0`$. The corresponding values of $`J`$, and hence of $`N_\chi `$ are listed in the second column of Table 2 assuming, as commonly done, that the dark matter density at the Sun’s galactocentric distance is equal to $`0.3\mathrm{GeV}\mathrm{cm}^3`$. One should however keep in mind that it has been questioned whether the NFW and the Moore et al. profiles can be used to describe inner dark matter halos (especially for smaller galaxies, see, e.g.,  for a review). These halo models give a snapshot of the Galaxy before the baryon infall. The appearance of a massive black hole at the Galactic center and of the stellar components may sensibly modify such pictures with further enhancements (but a depletion is possible as well) of the central dark matter density . For comparison, in Table 2 we give the value of $`J`$ for the modified isothermal sphere profile, which is non-singular toward the Galactic center and, as well known, give a normalization for the dark matter induced fluxes well below the background and the sensitivity of even next-generation detectors. All three halo profiles listed in the Table are consistent with available dynamical constraints on the Galaxy. We conclude then that $`N_\chi `$ can be at the level needed in our toy-model to reproduce the EGRET excess, but at the same time that both larger or smaller values are feasible as well.

## 3 EGRET Excess as Mapped by GLAST

Much more information on the nature of the EGRET excess at the Galactic center will be provided by the GLAST telescope. With respect to EGRET, GLAST will cover a wider energy range, with an increased effective area and better energy and angular resolution. Besides pinning down features in the WIMP-induced flux mediated by $`\pi ^0`$ decays, GLAST will have the power to search for the monochromatic gamma-ray flux eventually arising from the pair annihilation, at 1-loop level, of non-relativistic WIMP’s into a two-body final state containing a photon (for neutralinos as WIMP dark matter candidates, two such states are allowed: $`\chi \chi 2\gamma `$  and $`\chi \chi \gamma Z`$ , producing photons with energy equal to, respectively, $`M_\chi `$ and $`M_\chi (1M_z^2/4M_\chi ^2)`$). The discussion of GLAST potential to detect the monochromatic components is postponed to a future analysis; we focus here on the term with continuum energy spectrum.

We start by supposing that the GC excess as mapped by EGRET is indeed due to WIMP annihilations and extrapolate what kind of data GLAST would collect about it. For the performance of the detector, we rely on a simplified picture emerging from the latest simulations . We assume that, on average in the energy interval of interest to us, the instrument has an energy resolution of 10%, angular resolution of 10-5 sr ($``$0.1), and a peak effective area of 11000 cm2. GLAST will perform an all sky survey, rather than operating in the pointing mode. We assume a data acquisition time of 2 years and derive the fraction of time the GC center is visible by the instrument. We simulate a sky survey with $`\pm `$ 35 rocking and take into account the loss of exposure due to the South Atlantic Anomaly (SAA) passages, which corresponds to about 14.2% of the orbiting time. We find that the fraction of time that the GC is within the GLAST field of view is equal to 0.592, and that the net fraction of time that the source can be observed is 0.508; the reduction in effective area for sources which are not located at the instrument zenith gives a mean effective area equal to 60% of the peak effective area.

Fig. 5 shows a projection for the GC flux which could be measured by GLAST, assuming the spectrum and normalization for the dark matter source and the normalization for the background as derived from the fit of the EGRET data in Fig. 2 (lower panel). The simulated data points are derived by assuming as angular acceptance the EGRET angular resolution ($`\mathrm{\Delta }\mathrm{\Omega }10^3`$ sr), and choosing the energy bin widths to be of the order of $`10\%`$ of their central values to take advantage of the GLAST energy resolution. The error bars displayed are associated to the statistical error only. If we try to fit the simulated data with the spectral shape of the background and a free normalization, the reduced $`\chi ^2`$ we get is higher than $`210^3`$, much larger than the reduced $`\chi ^2`$ we obtained in the corresponding fit for the EGRET data set.

GLAST will also allow to search for an eventual angular signature of the dark matter source. A 10-5 sr angular resolution implies that the telescope will map the GC resolving regions with a precision of 7 pc (assuming the sun galactocentric distance is 8.5 kpc).

Most models for the distribution of dark matter in the Galaxy, such as the NFW and the Moore et al. models, predict an enhancement in the dark matter distribution toward the Galactic center on a scale which exceeds this size. There is then the chance that WIMP annihilations in the GC region may produce a detectable flux on an angular scale exceeding the angular resolution of GLAST. It is then sensible to investigate whether GLAST will detect the flux in the example shown in the Fig. 5 as coming from a point source located at the GC or from a diffuse source with degrading intensity increasing the angle between the line of sight and the GC direction. To do that we need to focus on a specific model for the dark matter density profile in the inner Galaxy. Inspired by the NFW and the Moore et al. profiles, we assume that the WIMP density toward the GC scales like:

$$\rho (r)=\{\begin{array}{cc}\rho _0\left(\frac{r}{r_0}\right)^\gamma ,r>r_{min}\hfill & \\ \rho _0\left(\frac{r_{min}}{r_0}\right)^\gamma ,rr_{min}\hfill & \end{array}.$$
(9)
We normalize $`\rho (r)`$ by fixing the local WIMP density, i.e. the density at the galactocentric distance $`r_0=8.5`$ kpc, to be equal to $`\rho _0=0.3`$ GeV/cm3. $`\gamma `$ is kept as a free parameter. To avoid the singularity in $`r=0`$ we have introduced a lower cut-off $`r_{min}=10^5`$ kpc, corresponding to a distance from the GC below which we assume that the power law behavior cannot be trusted.

For the sample toy-model shown in Fig. 2 (lower panel) and in Fig. 5, $`N_\chi =J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^3\mathrm{sr})=8.510^4`$ from which we infer that $`\gamma =1.54`$. Values of $`J(\psi )_{\mathrm{\Delta }\mathrm{\Omega }}`$, for such $`\gamma `$ and for the dark matter density profile specified in Eq.(9), are shown in Fig. 6 as a function of $`\psi `$ and for a few values of $`\mathrm{\Delta }\mathrm{\Omega }`$. We can then calculate the expected flux for the model obtained in the fit of EGRET data, fixing the direction of observation and the angular acceptance. Finally we compute the expected flux that, with the above provisions, GLAST will collect. We have also included a background component independent from $`\psi `$ which is equal to that we have estimated from the fit of the data by EGRET.

In Fig. 7, as a function of the direction $`\psi `$ and for a few values of $`\mathrm{\Delta }\mathrm{\Omega }`$, we plot the value of the reduced $`\chi ^2`$ we get by fitting in each case the expected flux with a component that has the spectral shape of the background and free normalization. As $`\chi ^21`$ on angular scales much larger than the angular resolution, this indicates that GLAST will resolve the angular structure of the dark matter source at the GC.

¿From the figure we see also that, for $`\psi =0`$, the largest $`\chi ^2`$ are obtained for the minimum $`\mathrm{\Delta }\mathrm{\Omega }`$ considered, i.e. an angular acceptance equal to the angular resolution of the detector. This indicates that, for the specific halo model considered, the ratio of the signal to square root of the background increases going to smaller and smaller angular acceptances.

Values of the reduced $`\chi ^2`$ have been obtained so far by supposing that the spectral shape of the background is known. Actually, what we have done is to fix it according to one of the currently favored scenarios. Slight discrepancies with respect to this model are plausible. On the other hand, GLAST will perform an all sky survey in which the background component will be accurately measured at all longitudes and latitudes. It might still be problematic to choose the background normalization in the Galactic center direction, if an excess is indeed found in that direction. But it will be possible to relax the assumption that the spectral shape of the background is known from theory, as it will be possible to extrapolate it from the data at higher latitudes and longitudes. A (measured) spectral shape which is different from what we assumed would slightly change the numerical predictions we derived so far; on the other hand the general features would remain exactly the same. Keeping this in mind, one should not take a reduced $`\chi ^2`$ of 1 as a strict discriminator, e.g. in Fig. 7, to decide whether the signal could be resolved from the background.

## 4 GLAST Performance for a Weaker Source

Although we have shown that the flux measured by EGRET is compatible with being due to WIMP annihilations, until more accurate data are available it will not be possible to discriminate this solution from other plausible scenarios. In particular, given the rather poor angular resolution of the EGRET instrument, there is even the possibility that the EGRET excess is not actually associated to a source located at the Galactic center. This is the case, e.g., if the flux is identified with inverse Compton emission from the electrons responsible for the synchrotron emission in the radio arc around the Galactic center . Should the next generation of gamma-ray observations solve the puzzle in this direction, there still would be room for finding a component due to WIMP annihilations in the Galactic center region, maybe associated to a dark matter source weaker than the one we postulated so far. Hence we extend the analysis we performed to investigate the potential of GLAST to single out such a source.

As we have shown, for a halo profile of the type given in Eq.(9), it is advantageous to focus on a region which is as small as possible around the GC. Hence we consider a survey with angular acceptance equal to the angular resolution of GLAST, $`10^5`$ sr. We assume that the background component is still due to diffuse emission from primary cosmic rays interacting with the interstellar medium. Hence we keep the spectral form implemented so far for $`S_b`$, with a normalization such that it matches at least the higher latitude measured flux, i.e. $`N_b>3.210^{20}`$ for a background at least at the level of the flux reported in the second column in Table 1. For a given WIMP model, i.e. fixing $`S_\chi `$, we search then for the minimum ratio between the two normalization factors $`N_\chi /N_b`$ that is needed to eventually discriminate with the GLAST telescope the WIMP annihilation signal from the background. $`N_\chi `$ is now equal to $`J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^5\mathrm{sr})`$: sample values for this quantity are listed in the third column of Table 2 for the three halo profiles introduced in Section 2.3. Again, the huge spread in the predictions (possibly even further amplified by the redistribution of dark matter particles during the formation of the central black hole ) reflects our lack of knowledge about the dark matter density in the Galactic center region.

For each pair $`N_\chi `$ and $`N_b`$, assuming the same instrumental performance, exposure time and energy binning specified in Section 3, we simulate the corresponding data set that GLAST would obtain, i.e. the expected flux measurements with the associated statistical error for the chosen energy binning, analogously to what is shown in Fig. 5 for the EGRET GC source. The criterion to discriminate whether the WIMP component would be singled out is based on the usual $`\chi ^2`$ test statistic. We have computed the reduced $`\chi ^2`$ between the number of counts expected in each energy bin for the two hypothesis: WIMP signal plus background and background only. Taking into account the number of degrees of freedom, which in our case is equal to the number of energy bins, the signal plus background curve is distinguishable from the background only curve, for a reduced $`\chi ^2>constant`$. This constant in uniquely determined by the number of degrees of freedom and by the confidence level we want to reach. We have also checked our results against those obtained with the likelihood ratio method , obtaining no discrepancies.222We thank G.Ganis for providing the software package we used to perform this analysis. This latter method is especially suited for the case we have at hand: to decide whether a certain event belongs to the background only hypothesis ($`H_0`$) or to signal plus background hypothesis ($`H_1`$), one starts by constructing two probability distributions, $`P_0`$ and $`P_1`$, for an estimator $`F=L(H_1)/L(H_0)`$, which is the ratio between the likelihoods $`L`$ of the two hypotheses. In our case, since we are interested in counting, we can choose the Poisson distribution to obtain the likelihood. Comparing the two distributions one can decide, at a certain confidence level, if they will result distinguishable or not, once it is fixed the accuracy of the experimental data that will be used for the discrimination. The likelihood ratio method is in general more powerful than the $`\chi ^2`$ one, since, in addition to giving the probability of a certain set of data to belong to the signal plus background probability distribution, it allows to compute the probability to be wrong when accepting such hypothesis, the so called power of the test, considering the background only hypothesis as the *true* one.

As a sample test case, we consider again the toy-model introduced in Section 2.3 for a single WIMP annihilation channel, e.g., $`b\overline{b}`$. In Fig. 8 we plot, as a function of the WIMP mass, the minimum value of the ratio $`N_\chi /N_b`$ needed for a discrimination of the WIMP signal from the background. The corresponding values of $`N_\chi =J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^5\mathrm{sr})`$ are of the order of $`10^3`$ for a WIMP mass of about 50 GeV. They increase approximately as $`M_\chi ^{1.8}`$ for heavier particles. These values are larger than the ones one would obtain from a smooth profile, see Table 2. A local enhancement in the WIMP dark matter density at the GC is then required to match our values of $`N_\chi `$. Such enhancement seems to be smaller than the one needed to fit the EGRET excess in the previous section ($`N_\chi `$ of the order of $`10^4÷10^5`$). As a word of caution we remark that the values of $`N_\chi `$ we obtained in the two analysis should not be directly compared in principle as the angular acceptance in the two cases is different: $`\mathrm{\Delta }\mathrm{\Omega }=10^5\mathrm{sr}`$ in the current case, while we had $`\mathrm{\Delta }\mathrm{\Omega }=10^3\mathrm{sr}`$ when we fitted the EGRET data set. A halo model has to be specified to translate one into the other.

## 5 A Specific WIMP: the Lightest Neutralino in the mSUGRA Framework

As a sample application of the generic tool we discussed so far, we focus now on the most widely studied WIMP dark matter candidate, the lightest neutralino, in the most restrictive supersymmetric extension of the Standard Model, the minimal supergravity (mSUGRA) framework . This setup has been considered extensively in the contest of dark matter detection (a list of recent references includes, e.g., ) and therefore the comparison of our results with previous work and other complementary techniques should be transparent in this case.

In the general framework of the minimal supersymmetric extension of the Standard Model (MSSM), the lightest neutralino is the lightest mass eigenstate obtained from the superposition of four interaction eigenstates, the supersymmetric partners of the neutral gauge bosons (the bino and the wino) and Higgs bosons (two Higgsinos). Its mass, composition and couplings with Standard Model particles and other superpartners are a function of the several free parameters one needs to introduce to define such supersymmetric extension. In the mSUGRA model, universality at the grand unification scale is imposed. With this assumption the number of free parameters is limited to five:

$$m_{1/2},m_0,sign(\mu ),A_0\mathrm{and}\mathrm{tan}\beta ,$$

where $`m_0`$ is the common scalar mass, $`m_{1/2}`$ is the common gaugino mass and $`A_0`$ is the proportionality factor between the supersymmetry breaking trilinear couplings and the Yukawa couplings. $`\mathrm{tan}\beta `$ denotes the ratio of the vacuum expectation values of the two neutral components of the SU(2) Higgs doublet, while the Higgs mixing $`\mu `$ is determined (up to a sign) by imposing the Electro-Weak Symmetry Breaking (EWSB) conditions at the weak scale. In this context the MSSM can be regarded as an effective low energy theory. The parameters at the weak energy scale are determined by the evolution of those at the unification scale, according to the renormalization group equations (RGEs). For this purpose, we have made use of the ISASUGRA RGE package in the ISAJET 7.64 software .

After fixing the five mSUGRA parameters at the unification scale, we extract from the ISASUGRA output the weak-scale supersymmetric mass spectrum and the relative mixings. Cases in which the lightest neutralino is not the lightest supersymmetric particle or there is no radiative EWSB are disregarded. The ISASUGRA output is then used as an input in the DarkSUSY package. The latter is exploited to:

* reject models which violate limits recommended by the Particle Data Group 2002 (PDG) ;
* compute the neutralino relic abundance, with full numerical solution of the density evolution equation including resonances, threshold effects and all possible coannihilation processes ;
* compute the neutralino annihilation rate at zero temperature in all kinematically allowed tree-level final states (including fermions, gauge bosons and Higgs bosons;
* estimate the induced gamma-ray yield by linking to the results of the simulations performed with the Lund Monte Carlo program Pythia  as implemented in the DarkSUSYpackage.

Note that none of the approximations implemented in the toy-model introduced above, regarding the estimate of the relic density or the annihilation cross section, are applied here.

We are ready then to exploit the procedure outlined in Section 4 to study in what region of the mSUGRA parameter space, for a given dark matter halo profile, the induced continuum $`\gamma `$ray flux would be detectable by GLAST. Fixing $`\mathrm{tan}\beta `$, $`A_0`$ and $`sgn(\mu )`$, we have performed a scan in the $`(m_0,m_{1/2})`$ plane searching for the minimum dark matter density, in the GC region, needed to be able to single out the neutralino annihilation signal with GLAST. To do this we have followed the same discrimination criteria described in Section 4 that we recapitulate here for the reader’s benefit. First we estimate the statistical error (1$`\sigma `$) on GLAST data to be the square root of the number of events. To compute the latter we multiply the flux by the effective area of the detector, by the total observational time and the angular resolution $`\mathrm{\Delta }\mathrm{\Omega }=10^5\mathrm{sr}`$. Then for each value of the pair $`(m_0,m_{1/2})`$ we compute the difference between the fluxes $`\varphi _\gamma =\varphi _b+\varphi _\chi =N_bS_b+N_\chi S_\chi `$ and $`\varphi _\gamma ^{}=\varphi _b=N_bS_b`$. If $`\varphi _\gamma \varphi _\gamma ^{}>3\sigma `$ we consider the SUSY model with those values of $`(m_0,m_{1/2})`$ to be detectable by GLAST. In Figs. 9 and 10 we show the isolevel curves for the minimum allowed value of $`N_\chi =J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^5\mathrm{sr})`$ for the signal detection, in the $`(m_0,m_{1/2})`$ plane and for five sample sets of the other parameters.

The colored regions in the figures represent portions of the parameter space where the neutralino has the right cosmological relic density to constitute CDM: in the green region $`0.13\mathrm{\Omega }_\chi h^20.3`$, while the red region corresponds to the WMAP range $`0.09\mathrm{\Omega }_\chi h^20.13`$ . It is interesting to remark that in the cosmologically favored regions the values of $`J`$ which allow the detection of the WIMP signal by GLAST are among the smallest in the interval of variability of the variable $`J`$ itself. To some extent this was anticipated in Eq.(8) that was one of the approximations of our toy-model: low values of the relic abundance lead to high values of the annihilation cross section.

The cosmologically favored regions correspond to two regimes which have been extensively studied in the literature. The first regime is given by models where the neutralino tends to be a very pure bino. It comprises the region where the neutralino is light (lower parts of each panel, compare with Fig. 12) and the region where the stau coannihilation is active (the low $`m_0`$ region on the left of the first three panels, see, e.g., ). The second regime, sometimes dubbed as “focus-point” region , is the region on the right-hand side of the panels with $`\mathrm{tan}\beta `$ equal to 10, 50 and 55, close to the region in parameter space where there is no EWSB, in which the neutralino has a relevant Higgsino component.

Comparing the values of $`J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^5\mathrm{sr})`$ obtained with those listed in Table 2, we see that in the case of a moderate enhancement of the dark matter density toward the GC as predicted by the NFW profile or an even slightly milder singular profile, there is a fair portion of the mSUGRA parameter space which gives fluxes detectable with GLAST. It is worth to observe that high $`\mathrm{tan}\beta `$ models (*i.e.* $`\mathrm{tan}\beta =50`$ and $`\mathrm{tan}\beta =55`$) possess the largest cosmologically favored regions. Furthermore they give the highest chances to single out the $`\gamma `$-ray signal from neutralino annihilations with GLAST. We can compare our results with those of . There it is assumed that GLAST has a certain sensitivity to the integrated continuum $`\gamma `$ray flux from a region around the GC, of an extension $`10^2`$ times wider than the GLAST angular resolution. The neutralino-induced signal is then supposed to be detectable if its integrated flux exceeds such sensitivity. In Figs. 18 and 19 of  are given the regions which can be probed by GLAST in case of $`J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^3\mathrm{sr})=500`$. They are in qualitative agreement with our corresponding predictions (first and fourth panels of Fig. 9) for $`\mathrm{\Delta }\mathrm{\Omega }=10^5`$ sr.

In addition to this study of the GLAST sensitivity, we have tried to single out the regions of the mSUGRA parameter space ($`(m_0,m_{1/2})`$ for fixed $`\mathrm{tan}(\beta )`$, $`A_0`$ and $`sign(\mu )`$), which are already experimentally excluded, due to a neutralino-induced $`\gamma `$-ray flux exceeding the GC EGRET data of table 1. In Fig. 11 we show isolevel curves for the maximal $`J(0)_{\mathrm{\Delta }\mathrm{\Omega }}(\mathrm{\Delta }\mathrm{\Omega }=10^3\mathrm{sr})`$ marginally consistent with the data at the $`5\sigma `$ level. We can observe that the values we obtain are generally significantly larger than the corresponding sample values in Table 9. This result implies that rather weak constraints on the mSUGRA parameter space can actually be imposed, on the basis of the GC $`\gamma `$ray flux measured by EGRET.

One could finally wonder how the GC EGRET data could be fitted in the context of the mSUGRA models. In order to answer to this question, we rely on the general analysis that we have performed in the toy-model scheme (see section 2). In particular we could extract from Fig. 4 the values of the neutralino mass that allow the best fit of the GC EGRET data333Recall that the possibility to fit EGRET data is quite insensitive to the dominant annihilation channel.. Looking at Fig. 12 we can then single out the $`(m_0,m_{1/2})`$ parameter regions that corresponds to such values.

## 6 Conclusions

We have found that the excess in the $`\gamma `$-ray flux detected by the EGRET telescope toward the Galactic center shows spectral features which are compatible with an exotic component due to WIMP annihilations, especially for WIMP masses in the lower end of the mass range currently considered for WIMP dark matter candidates. For the WIMP-induced flux to be at the level of the measured flux, a fairly large dark matter density is needed in the Galactic center region; indeed, such density enhancements are found in N-body simulations of halo profiles in cold dark matter cosmologies.

Although it is not possible with present data on the Galactic center excess to discriminate between the interpretation we propose here and other viable explanations, we have shown that, with the data that will be collected by the GLAST, the next major $`\gamma `$-ray telescope in space, it will be possible to identify both spectral and angular signature expected for a WIMP-induced component. If on the other hand the data will point to an alternative explanation, there will still be the chance for the GLAST telescope to single out a (weaker) dark matter source. The potentials of GLAST have been explored both in the contest of a generic simplified toy-model for WIMP dark matter and for one of the most widely studied WIMP dark matter candidate, the lightest neutralino, in the minimal supergravity framework. We find that even in case of moderately singular dark matter profiles, there are regions in the parameter space which will be probed by GLAST, especially in the high $`\mathrm{tan}\beta `$ case. We find, on the contrary, that limits from current EGRET data are rather weak.

## Acknowledgments

We would like to thank all the component of the GLAST Dark Matter working group for lots of discussion on the subject, in particular Elliot Bloom and Eduardo do Couto e Silva. We also thank Hans Mayer-Hasselwander for providing a table of the EGRET data toward the Galactic center and G.Ganis and C.Pittori for many useful discussions.


format: text

text:
# A first-principles study of oxygen vacancy pinning of domain walls in PbTiO3

## I INTRODUCTION

Ferroelectric materials are of intense interest for use in nonvolatile memory applications, in which the electric polarization of an array element is used to store a bit of information. However, the switchable polarization tends to decrease after many cycles of polarization reversal during device operation, a problem that is known as polarization fatigue. In the last decade, polarization fatigue in ferroelectrics has been under intensive study. Although several models have been proposed to explain this phenomenon, there is still no consensus, and many details of the fatigue process remain unclear. Nevertheless, it is generally believed that defects, especially charged defects, play an important role. For example, a series of experiments has provided some understanding of how such defects may pin the domain walls. In particular, attention has been drawn to oxygen vacancies, which are often the most common and mobile defects in perovskite ferroelectrics. Moreover, it has been found that fatigue resistance can be greatly improved by replacing the Pt electrodes with conductive-oxide electrodes, which can be explained in terms of the ability of oxide electrodes to control the concentration of oxygen vacancies in the sample. Furthermore, the fatigue rate in Pb(Zr,Ti)O3 films was found to be very sensitive to the oxygen partial pressure above the sample, suggesting that oxygen vacancies strongly affect the fatigue process.

Mechanisms for polarization fatigue based on pinning of domain walls by oxygen vacancies have been discussed phenomenologically by several authors. However, in order to put such phenomenologically theories on a firm atomistic basis, it is important to have detailed first-principles calculations that can provide information about the structure and energetics of the ferroelectric domain walls, of the oxygen vacancies, and of the interactions between the two.

While domain walls have been studied using Landau-type continuum theories in earlier works, first-principles calculations are essential for an accurate microscopic description of the domain walls. Most recently, Meyer and Vanderbilt studied 180 and 90 domain walls in PbTiO3 using first-principles methods,  establishing the geometry of the domain walls at the atomic level and calculating the creation energy of the domain walls. The domain wall was found to be extremely narrow, with a width of the order of the lattice constant $`a`$; the positions of the atoms change rapidly inside the domain wall and converge to their bulk value very quickly outside. As for oxygen vacancies, recent calculations have provided very useful information about the structure of oxygen-vacancy defects in this class of materials. However, to our knowledge, direct first-principles studies of the interactions between vacancies and ferroelectric domain walls have not yet appeared.

Thus, the goal of this work is to use first-principles calculations to examine how an oxygen vacancy would interact with a ferroelectric domain wall, and thus to shed some light on how oxygen vacancies might affect the switching process and cause polarization fatigue. Indeed, as the atomistic structure and polarization profile are very different in a domain wall than in the bulk of a ferroelectric domain, one may expect that the oxygen vacancies would also behave differently in such very different environments. To explore these issues, we adopt PbTiO3 as a model system for this study, and calculate the formation energies for neutral oxygen vacancies in the bulk and in 180 domain walls of tetragonal PbTiO3. Our calculations indicate that the vacancies do have a lower formation energy in the domain wall than in the bulk, thereby confirming the tendency of these defects to migrate to, and pin, the domain walls. The order of magnitude of the computed pinning energy is $`10^1`$eV, with substantial variations depending on geometrical configuration.

The rest of the paper is organized as follows. In Sec. II, we describe the technical details of our computational method and the supercells we used to model the domain walls. We present the results on the vacancy pinning energies of each of the three possible orientations of the original Ti–O–Ti bonds from first-principles calculations in Sec. III. A simple continuum model is introduced and used to help understand these results in Sec. IV. In Sec. V we discuss briefly the anticipated role of domain-wall pinning in fatigue, and discuss the case of charged vs. neutral defects. Finally, we summarize in Sec. VI.

## II METHODOLOGY

Our calculations are based on standard density-functional theory (DFT) within the local-density approximation (LDA). We use a planewave pseudopotential code implemented by the Vienna ab-initio Simulation Package (VASP). Vanderbilt ultrasoft pseudopotentials  are used, with Pb $`5d`$ and Ti $`3s`$, $`3p`$ electrons included explicitly. The 29 Ry cut-off used here was well tested and determined to be adequate for this material and these pseudopotentials. The structure is considered to be relaxed when the forces are less than 0.05 eV/Å; the change of total energy at this time is typically less than 1 meV.

We first carried out reference calculations for domain walls without vacancies, following closely the approach of Ref. . We describe the 180 domain wall using an 8$`\times `$1$`\times `$1 supercell (long direction along $`\widehat{x}`$) with 2$`\times `$4$`\times `$4 k-point sampling; the tetragonal axis and polarization are along $`\pm \widehat{z}`$, and two $`yz`$-oriented domain walls divide the supercell into up and down domains of equal width.

We also carried out reference calculations for vacancies in the bulk. There are three kinds of oxygen vacancies that can be formed, depending on whether the removed oxygen atom had its Ti–O bonds along $`\widehat{x}`$, $`\widehat{y}`$, or $`\widehat{z}`$, which we denote as V($`x`$), V($`y`$), or V($`z`$), respectively. These bulk vacancies were studied using a 2$`\times `$2$`\times `$2 supercell with the tetragonal axis (and polarization) chosen along $`\widehat{z}`$, and with 2$`\times `$2$`\times `$2 k-point sampling. (Thus, V($`x`$) and V($`y`$) are equivalent by symmetry in the bulk.) To obtain the vacancy energy, we first calculate the total energy of the pure tetragonal 2$`\times `$2$`\times `$2 supercell as the reference energy. One oxygen is then removed from the supercell in such a way that the supercell remains net neutral, and the resulting structure is relaxed. We keep the volume and shape of the supercell fixed and allow only the ion positions to relax. The vacancy energy is then calculated by comparing the total energy difference before and after removing the oxygen atoms from the supercell.

In order to study the interaction of the vacancies with the 180 domain wall, we constructed an 8$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell by doubling the 8$`\times `$1$`\times `$1 supercell in the $`y`$-$`z`$ plane in a c$`2`$$`\times `$$`2`$ or $`\sqrt{2}`$$`\times `$$`\sqrt{2}`$R(45) arrangement, and then removing one oxygen atom (and its periodic images) from the interior of the domain-wall structure. This is illustrated schematically in Fig. 1, referring here to vacancies labeled “Vdw” (i.e., located in the domain wall) in Fig. 1(a). The vacancies form sheets in the $`y`$-$`z`$ plane, with individual vacancies separated by a distance of about $`a\sqrt{2}`$, as indicated in Fig. 1(b). Even sparser arrangements would be desirable, but the supercell already contains 80 atoms, and computational limitations make it difficult to increase the separation further. A 1$`\times `$3$`\times `$3 k-point mesh is used for these supercells, and once again the ion positions are allowed to relax while keeping the supercell volume and shape fixed. (A single isolated domain wall decorated by vacancies could not expand or contract in the $`\widehat{y}`$ or $`\widehat{z}`$ directions because of the epitaxy constraint to the bulk, and our experience, consistent with Ref. , indicates that the lattice would expand very little in the $`\widehat{x}`$ direction if allowed to do so.)

Finally, in order to calculate the energy difference for an oxygen vacancy to be inside the domain wall, relative to being in bulk, we find that it is advantageous to carry out corresponding 8$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell calculations in which the oxygen vacancy has been removed from a bulk-like region of the supercell. In this way, we reduce the systematic errors associated with supercell shape, k-point sampling, etc. Our results for the binding energies of oxygen vacancies in the domain walls will normally be based on calculations of this kind, except where unavailable as noted below.

## III RESULTS

We first report our calculations on isolated vacancies in bulk ferroelectric PbTiO3 in the 2$`\times `$2$`\times `$2 supercell. Using the theoretical values of the lattice constants ($`a=3.86`$ Å, $`c/a=1.0466`$) as obtained in Ref. , our calculations show that the oxygen vacancies of type V($`z`$) are more stable than V($`x`$) and V($`y`$) by 0.3 eV, similar to what was found by Park and Chadi for somewhat different conditions. As indicated in Sec. II, however, we prefer to calculate pinning energies for oxygen vacancies in domain walls by using reference bulk vacancy calculations in an 8$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell that can be used with and without domain walls, in order to provide maximum cancellation of systematic errors. Such calculations will be the basis for the the results given in this section. However, we will return to the use of the 2$`\times `$2$`\times `$2 supercell in Sec. IV for some calculations relevant to the modeling of vacancies in a perturbed environment.

We then confirm the structure of the relaxed 180 domain wall. Our results for the 180 domain wall are very close to those of the previous calculation. The domain wall is located on the Pb-O plane. The atomic displacements converge rapidly to their bulk values, with most of the displacements ocurring within one unit cell of the domain wall center. Consequently, the polarization reverses sharply within approximately one unit cell of the domain wall, and quickly saturates to a bulk value further away. Thus, the domain wall is extremely narrow, only about two lattice constant wide, and we find that it displays a significant $`x`$-$`z`$ shear in addition to the reduced polarization. With this justification, we can heuristically divide the supercell into two different regions, a bulk region and a domain-wall region, as sketched in oversimplified form in Fig. 1(a).

To minimize the interactions between neighboring vacancies, the domain-wall supercell is doubled as described in Sec. II. The shape of the supercell in the $`y`$-$`z`$ plane is indicated by the dashed lines in Fig. 1(b). The total energy E0 of this domain wall structure (80-atom $`8`$$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell) is calculated for use as the reference energy.

To calculate the oxygen vacancy energy, we remove one oxygen from the supercell and relax the resulting structure. Recall that there are three types of oxygen vacancies, V($`x`$), V($`y`$), and V($`z`$), according to the orientation of the Ti–O–Ti bond from which the oxygen atom has been removed. For each type, we choose one vacancy as close as possible to the center of the bulk region, and another as close as possible to the central domain-wall plane. We label these as ‘b’ (bulk) and ‘dw’ (domain wall), respectively. Thus, Vdw($`x`$) is an $`x`$-oriented vacancy at the domain wall, etc. In all calculations, we keep the supercell charge-neutral.

All unrelaxed vacancies have an $`M_y`$ mirror symmetry through the defect site. In addition, the unrelaxed Vdw($`x`$) and Vb($`x`$) defects have extra $`C_2^{(y)}`$ and $`M_x`$ symmetries respectively, because they lie precisely in, or else half-way between, the domain-wall planes. When relaxing the structure, we observe that the two Ti ions that neighbor the vacancy site relax most, as expected. In the bulk region, the relaxation of Vb($`x`$) and Vb($`z`$) does not lead to any symmetry breaking, while for Vb($`y`$) there are Pb ion displacements that break the $`M_y`$ symmetry and contribute about 10 meV to the relaxation energy. In the domain-wall region, Vdw($`y`$) also has similar distortions, but the energy is only lowered by about 5 meV. As for Vdw($`x`$), which starts with a relatively high unrelaxed symmetry (four-element point group C2h), the relaxation breaks the $`M_y`$ and $`C_2^{(y)}`$ symmetries so that the only remaining symmetry is inversion (point group C1h), and the total energy is lowered by about 20 meV.

The total energy of the oxygen-vacancy supercell is then calculated for each configuration, and a vacancy formation energy is computed using an appropriate reference. The results are listed in Table I. Here $`\mathrm{\Delta }E_{\mathrm{dw}}=E_{\mathrm{v}\mathrm{dw}}+E_{\mathrm{oxy}}E_{\mathrm{dw}}`$, where $`E_{\mathrm{v}\mathrm{dw}}`$ is the energy of the supercell with the vacancy in the domain wall, $`E_{\mathrm{oxy}}`$ is the energy of a free oxygen atom, and $`E_{\mathrm{dw}}`$ is the energy of a vacancy-free supercell with domain walls; and either $`\mathrm{\Delta }E_\mathrm{b}=E_{\mathrm{v}+\mathrm{dw}}+E_{\mathrm{oxy}}E_{\mathrm{dw}}`$ (“Ref. with DW”) or $`\mathrm{\Delta }E_\mathrm{b}=E_\mathrm{v}+E_{\mathrm{oxy}}E_0`$ (“Ref. without DW”), where $`E_{\mathrm{v}+\mathrm{dw}}`$ is the energy of a vacancy in the bulk-like region of a supercell containing domain walls, $`E_\mathrm{v}`$ is the energy of a vacancy in a domain-wall-free supercell, and $`E_0`$ is the energy of the corresponding bulk supercell containing neither domain walls nor vacancy.

As shown in Table I, all the vacancies in the domain wall have lower formation energies than their counterparts in the bulk. This indicates that there is an attractive interaction between the vacancy and the domain wall and implies a positive pinning energy $`E_\mathrm{p}=\mathrm{\Delta }E_\mathrm{b}\mathrm{\Delta }E_{\mathrm{dw}}`$.

Where available, the values for $`E_\mathrm{p}`$ obtained using the reference supercell with domain walls is to be preferred (last two columns of Table I), since one expects a more systematic cancellation of errors in this case. However, in some situations this turns out not to be possible. For the vacancy V($`z`$), in particular, a problem arises. We find that when we attempt to compute the relaxed energy $`E_{\mathrm{v}+\mathrm{dw}}`$ of the supercell in which the vacancy has been placed as far from the domain walls as possible, the domain wall actually shifts its position during the relaxation in order to coincide with the vacancy, thus spontaneously converting the supercell from $`E_{\mathrm{v}+\mathrm{dw}}`$ to $`E_{\mathrm{v}\mathrm{dw}}`$. This will be discussed further in Sec. III A. For this case, our best value for $`E_\mathrm{p}`$ of 261 meV is obtained by falling back to the use of a reference supercell without domain walls (middle columns of Table I). By looking at other cases (i.e., V($`x`$), V($`y`$), and unrelaxed cases), it can be seen that an uncertainty of approximately 25 meV is introduced by the use of the less preferable reference supercell.

It may be noted in Table I that the formation energies $`\mathrm{\Delta }E_\mathrm{b}`$ are significantly different for Vb($`x`$) and Vb($`y`$) ($``$200 meV), although by symmetry they should be equal in a true bulk environment. Here the differences arise from a supercell size effect connected with the arrangement of vacancies into sheets of fairly high density in the $`y`$$`z`$ plane. For a sufficiently large supercell we would expect these energies to become equal, because the local environments of Vb($`x`$) and Vb($`y`$) would be almost identical and the interactions between them would be negligible. However, in our case the vacancies are only separated by a distance of about $`a\sqrt{2}`$, so the interactions are not negligible. On the other hand, the vacancies in the domain walls should have similar interactions, and we can expect some cancellation of errors when arriving at the pinning energy. Thus, we have more confidence in the $`E_\mathrm{p}`$ values than in the relative formation energies of V($`x`$), V($`y`$) and V($`z`$).

### A Domain-wall shift

As indicated in the previous subsection, when we attempt to calculate the energy $`E_{\mathrm{v}+\mathrm{dw}}`$ of the Vb($`z`$) vacancy in the bulk-like region of an $`8`$$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell containing domain walls, the domain wall spontaneously shifts toward the vacancy during the relaxation process. We carried out further tests using a $`10`$$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell and observed the same phenomenon, as shown in Fig. 2. First, we put the vacancy on a TiO2 plane inside the domain wall ($`x=4.5a`$) and allow relaxation. We observe that the Pb-centered domain wall shifts towards the vacancy and becomes a Ti-centered domain wall. We then put the vacancy on a TiO2 plane near the domain wall ($`x=3.5a`$), and observe that the domain wall center (originally at $`x=5a`$) again shifts to the left and becomes centered roughly at the TiO2 plane at $`x=3.5a`$, ending up with almost same total energy as before. When we attempt to put the vacancy even farther from the domain wall at $`x=2.5a`$, we find that a new domain wall forms at the vacancy position. Clearly, the domain wall is trying to shift its position in each case so as to minimize the polarization at the position of the vacancy, thereby demonstrating directly the pinning effect of oxygen vacancies and ferroelectric domain walls.

This effect is most pronounced for the case of Vb($`z`$) because it has the strongest pinning energy, as can be seen in Table I. In the case of V$`(y)`$, we find a similar but weaker effect. That is, if the V$`(y)`$ vacancy is placed close enough to the domain wall (e.g., at $`x=4.5a`$ in Fig. 2), a similar shift can occur; but no shift occurs if the defect is placed farther away.

To pursue the calculation of $`E_{\mathrm{v}+\mathrm{dw}}`$ in order to obtain $`\mathrm{\Delta }E_\mathrm{b}`$ for Vb($`z`$), it would be necessary for us to use a supercell larger than $`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ in the $`y`$$`z`$ plane. However, as this would be computationally prohibitive, we have instead chosen to recalculate the bulk vacancy energies using an $`8`$$`\times `$$`\sqrt{2}`$$`\times `$$`\sqrt{2}`$ supercell without domain walls, as indicated in the previous subsection. This provides an alternative reference energy which, though less accurate, is available in all cases.

## IV MODEL CALCULATION

Our first-principles calculations give us a rough picture of the interactions between oxygen vacancies and domain walls. These calculations show that the domain walls can indeed be pinned by oxygen vacancies. We obtain estimates for the pinning energies, and find that the pinning energy for V($`z`$) is much larger than for V($`x`$) or V($`y`$). We would like to understand better the physics underlying these results, and to appreciate which results might generalize to other situations (e.g., other ferroelectrics materials, or other domain-wall structures such as 90 boundaries).

With this motivation, we consider a model in which the vacancy formation energy depends on the immediate vacancy environment as characterized by local polarizations and strains. In particular, we consider a continuum description of the polarization and strain fields in the domain wall, as would occur in a Landau-type model, and assume that the vacancy energy can be expressed as a function of the local strain and polarization only. (A more sophisticated model might involve also a dependence on the local gradients of these fields, but we have not pursued this here.) Thus, in general we would write the vacancy formation energy as $`E_\mathrm{v}(\eta ,𝐏)`$, where $`\eta `$ and $`𝐏`$ are the strain tensor and polarization vector describing the state of the local environment before removal of the oxygen atom. However, we specialize here to the case of interest, a 180 wall lying in a $`y`$$`z`$ plane separating tetragonal phases with polarizations along $`\pm \widehat{z}`$ on either side. Thus, we focus on only the $`z`$ component of polarization, and for convenience we define a dimensionless reduced polarization $`p=P_z/P_{\mathrm{bulk}}`$. As for the strain tensor, we have $`\eta _{yy}=\eta _{yz}=\eta _{zz}=0`$ by lattice continuity and $`\eta _{xy}=0`$ by symmetry. Moreover, we find that $`\eta _{xx}`$ remains quite small in the domain-wall region. Thus, we focus only on the $`xz`$ shear strain component and let $`\eta _s=(a/c)\eta _{xz}`$. The vacancy formation energy $`E_\mathrm{v}(\eta _s,p)`$ is thus considered as a function of the local $`\eta _s`$ and $`p`$.

Expanding $`E(\eta _s,p)`$ in powers of $`p`$,

$$E_\mathrm{v}(\eta _s,p)=A(\eta _s)+B(\eta _s)p^2+𝒪(p^4),$$
(1)
where odd powers in $`p`$ have been dropped by symmetry, and only terms up to $`𝒪(p^2)`$ are retained henceforth. We then expand $`A(\eta _s)`$ and $`B(\eta _s)`$ in powers of $`\eta _s`$ as

$$A(\eta _s)=a_0+a_1\eta _s+a_2\eta _s^2+𝒪(\eta _s^3)$$
(2)
and

$$B(\eta _s)=b_0+b_1\eta _s+b_2\eta _s^2+𝒪(\eta _s^3).$$
(3)
Dropping terms beyond quadratic order in $`\eta _s`$, we take the coefficients $`a_1`$, $`a_2`$, $`a_3`$, $`b_1`$, $`b_2`$ and $`b_3`$, in Eqs. (1-3) to constitute the parameters of our model.

To obtain these six parameters, we do a series of calculations at a set of different $`\eta _s`$ values for both $`p=0`$ and $`p=1`$. To do this, we construct $`2`$$`\times `$$`2`$$`\times `$$`2`$ supercells at different fixed values of $`\eta _s`$ and calculate the vacancy energies as in the last section, allowing relaxation of ionic positions but not strains. A $`2`$$`\times `$$`2`$$`\times `$$`2`$ k-mesh is used in all these calculations. We do calculations at $`p=0`$ by enforcing inversion symmetry of the lattice. If there is no constraint of symmetry imposed, the ions relax freely, resulting in $`p=1`$. The results of these calculations are shown in Fig. 3.

It might be thought that the coefficients $`a_1`$ and $`b_1`$ in Eqs. (2-3) should vanish by symmetry, but if the oxygen vacancy induces a distortion which lowers the lattice symmetry as discussed in Sec. III, this may not always be true. Consider, for example, the case of V($`x`$) at $`\eta _s=0`$ and $`p=0`$ (inversion symmetry imposed). While the unrelaxed defect has D2h symmetry, the relaxation lowers the symmetry to C2h ($`E`$, $`I`$, $`M_y`$ and $`C_2^y`$) and reduces the energy by about 65 meV relative to the case where no symmetry breaking is allowed. Actually, there are two equivalent defects, related to each other via the broken symmetry $`M_x`$, having degenerate energies at $`\eta _s=0`$ but having $`a_1`$ coefficients of opposite sign (i.e., opposite response to an applied $`xz`$ strain $`\eta _s`$.) In Fig. (3) only the energy of the more stable of these two defects is plotted as a function of $`\eta _s>0`$, and the dashed line is a fit to Eq. (2). The resulting (negative) value of $`a_1`$ is given in Table II.

Considering the other vacancies at $`p=0`$, vacancy V($`y`$) has a similar symmetry breaking but its orientation is such that the degeneracy would be split by a $`\eta _{yz}`$ strain, which is absent here. For V($`z`$) we observe no symmetry breaking at $`p=0`$. Thus $`a_1`$ vanishes for these cases.

Turning now to $`p=1`$, we find a reversed situation: V($`x`$) and V($`y`$) show no breaking of their C2v symmetry at $`\eta _s=0`$, while V($`z`$) breaks from C4v to C1h after relaxation. Thus, $`b_1=0`$ for V($`y`$) but not for V($`z`$). For non-zero $`\eta _s`$ in Fig. (3), no further symmetry breaking is observed beyond what was already present at $`\eta _s=0`$.

The parameters resulting from all the fits of Eqs. (1-3) are listed in Table II. From the fact that $`a_1`$ and $`a_2`$ are negative, we see that the vacancies will prefer an environment of high shear strain. Similarly, since $`b_0`$ is positive, they will prefer an environment of low polarization. Thus, the parameters are suggestive of a tendency for the vacancies to pin the domain walls.

In order to model quantitatively the vacancy formation energy in the domain wall, we now have to estimate the values of $`p`$ and $`\eta _s`$ that occur in a vacancy-free domain wall at the location where the vacancy would occur. Since Vdw($`x`$) lies in the Pb-O symmetry plane, $`p=0`$ there. From the first-principles calculations of Ref. , we can estimate that $`p0.8`$ at the neighboring TiO2 plane where Vdw($`y`$) and Vdw($`z`$) are located.

The estimation of $`\eta _s`$ is more subtle. The problem is that the shear strain is not well defined in the domain walls, since it depends strongly on which of the sublattice we follow. For example, if we define the shear strain to be $`\eta _s=\delta _z/c`$, where $`\delta _z`$, is the displacement of the adjacent atoms of same type in the $`\widehat{z}`$ direction, we estimate that $`\eta _s(\mathrm{Pb})=0.068`$, $`\eta _s(\mathrm{Ti})=0.054`$, $`\eta _s(\mathrm{O}(\mathrm{x}))=0.029`$, $`\eta _s(\mathrm{O}(\mathrm{y}))=0.047`$, and $`\eta _s(\mathrm{O}(\mathrm{z}))=0.027`$ in the center of the domain wall. This variation reflects the reversal of the polarization-related displacements along $`\widehat{z}`$ as one passes through the domain wall along $`\widehat{x}`$. We could define a mean shear as $`\widehat{\eta }_s=(1/5)_i\eta _s(i)0.005`$, which as expected is about half of the “geometrical offset” defined in Ref.  (the offset occurs over approximately two unit cells). An alternative choice is the root-mean-square shear strain $`\overline{\eta }_s=[(1/5)_i\eta _s^2(i)]^{1/2}=0.048`$. We expect that the most reasonable choice of an effective $`\eta _s^{\mathrm{eff}}`$ should lie somewhere between these two limits. For Vdw($`x`$), the pinning energy we get from first-principles calculation is 97 meV. Comparing with Eqs. (1-3), we find that $`\eta _s^{\mathrm{eff}}=0.03`$ gives a reasonable agreement with the first-principles result for this case, and we thus adopt this value. The shear strain should be slightly smaller at the location of Vdw($`y`$) or Vdw($`z`$), half a lattice constant away from the domain-wall center, but for simplicity we retain the same value of $`\eta _s=0.03`$ for all three defects.

Using the parameters from Table II and the values of $`\eta _s`$ and $`p`$ discussed above, we may estimate the pinning energy $`E_\mathrm{p}=E_\mathrm{v}(0,1)E_\mathrm{v}(\eta _s,p)`$ via

$$E_\mathrm{p}=(1p^2)b_0(a_1+b_1p^2)\eta _s(a_2+b_2p^2)\eta _s^2.$$
(4)
The resulting pinning energies are reported, and compared with the direct first-principles calculations, in Table III.

Recall that the good agreement for V($`x`$) occurs by construction. The agreement for V($`z`$) is fair and that for V($`y`$) is somewhat poor, although at least we have the correct sign of the pinning energy even in the worst case of V($`y`$). Thus, we find that the present system is sufficiently complex that our simple model description is only partially successful.

There are several reasons why this may be. As discussed in Sec. III A, the domain walls may shift if the density of oxygen vacancies is high, and the shift increases the pinning energy. This may help explain the underestimation of the pinning energies for V($`y`$) and V($`z`$). Also, the use of a 2$`\times `$2$`\times `$2 supercell for the calculations of Fig. 3, from which the values of Table II were obtained, means that the vacancies were much closer to the dilute limit than was the case for the vacancy-in-domain-wall calculations of Table I. This is most likely the dominant source of the discrepancy for the case of Vdw($`y`$), which is less sensitive to $`\eta _s`$. Indeed, the fact that a higher density of oxygen vacancies in the domain wall leads to a larger pinning energy is consistent with a picture in which oxygen vacancies would tend to make planar clusters in the domain wall, thus acting to increase the pinning energy.  Future tests on larger supercells in the $`y`$$`z`$ directions might help clarify these issues, although these still remain intractable for the time being.

Another obvious source of the discrepancies may be the limitations of the model. It is unsatisfying that the choice of $`\eta _s`$ is so ambiguous, and it is unclear whether two variables ($`\eta _s`$ and $`p`$) should suffice to describe the local environment. After all, the structural distortions change rapidly as one passes through the domain wall, so that it is not clear whether a Landau-type continuum model should be expected to capture the details of the energetics. It might be interesting to see whether a model more like the effective-Hamiltonian description of Bellaiche et al., which includes compositional disorder in order to treat alloys, could be successfully used here.

Nevertheless, we believe that our model captures some of the essential physics of the pinning mechanism of oxygen vacancies in the 180 domain walls. It gives the correct sign and overall order of magnitude for the pinning energy, and correctly reflects that the pinning is strongest for Vdw($`z`$), intermediate for Vdw($`x`$), and weakest for Vdw($`y`$). It also helps clarify the relative roles of strain and polarization effects in the pinning mechanism. We thus expect that it may be of some use for understanding other ferroelectric materials as well.

## V DISCUSSION

We now briefly discuss how oxygen vacancies may affect the ferroelectric switching process. As is well known, switching in ferroelectrics occurs not through a homogeneous concerted reversal of the polarization in the bulk, but through the motion of domain walls separating regions of different polarization. Thus, insofar as these domain walls become pinned, the switching will be suppressed.

In a pristine ferroelectric material that has a robust hysteresis loop and a large remanant polarization, the ferroelectric domain walls are presumably only weakly pinned by some pre-existing defects. Our calculations indicate that oxygen vacancies will tend to migrate into these domain walls over time, since they experience a binding energy to the domain wall of between 100 and 250 meV. In fact, once inside the domain wall, we would expect the vacancy to hop into the V($`z`$) configuration, since this is lower in energy than the V($`x`$) or V($`y`$) configurations. Thus, the effective binding energy is $``$250 meV, the value associated with the V($`z`$) vacancy.

If a significant number of vacancies accumulate in the domain wall, they in turn can act to pin the domain wall. When the density of such vacancies is low, they will not pin the domain walls strongly, and switching will still be able to occur. As the areal density of vacancies increases, however, an increasing fraction of domain walls (or increasingly large portions of individual domain walls) may become immobile, resulting in the decay of the switchable polarization.

Of course, there are many limitations of our theoretical analysis, and the real experimental situation could be much more complicated. We have studied neutral oxygen vacancies, which correspond to vacancies of charge +2$`e`$ neutralized by electrons residing in nearby states of mainly Ti $`3d`$ character. In the absence of neutralizing electrons, the vacancies may pin the domain walls even more strongly. (In this case, the situation becomes more complicated, since we can expect that the domain walls may acquire a tilt in order to compensate the vacancy charges. That is, if the 180 domain wall is not exactly parallel to the polarization, there is a bound charge $`\mathrm{\Delta }𝐏\widehat{n}`$ that can help neutralize the vacancy charges, an effect which may contribute to the strength of the pinning effect.) Alternatively, the oxygen vacancies may tend to aggregate into clusters or to form defect complexes of various kinds. Finally, the polarization switching process is rather a complicated dynamical process that we have not attempted to model in detail.

Nevertheless, we believe that our first-principles results serve as a first step towards understanding the possible role of oxygen vacancies in the pinning of ferroelectric domain walls. While we have investigated only one class of defects that may be involved in pinning, at least our calculations provide a lower limit for the strength of the pinning effect, which may be stronger if other defects or defect complexes play the dominant role. Our results may also serve as input for more complex modeling and simulation. For example, it could be used to extend a model such as that of Ahluwalia and Cao, who have done simulations of domain-wall formation in a 2D model simulation, by the inclusion of vacancies into the simulation. Our results might also be useful in the formulation of an effective-Hamiltonian approach that could be used to carry out finite-temperature simulations of the domain-wall behavior. Such studies might help quantify the effective strength of the pinning effect under more realistic conditions.

## VI SUMMARY

In summary, we have used first-principles density-functional calculations to investigate the interaction of oxygen vacancies and 180 domain walls in tetragonal PbTiO3. Our calculations indicate that the vacancies do have a lower formation energy in the domain wall than in the bulk, thereby confirming the tendency of these defects to migrate to, and pin, the domain walls. The pinning energies are calculated for each of the three possible orientations of the original Ti–O–Ti bonds, and are found to be 97 meV, 76 meV and 261 meV for V($`x`$), V($`y`$) and V($`z`$) respectively. We also introduce a simple continuum model with only two parameters ($`p`$, $`\eta _s`$) to model the results. This simple model gives pinning energies that agree qualitatively with the first-principles calculations, and we expect that it may prove useful for other ferroelectric systems as well.

This work was supported by ONR grant N0014-97-1-0048.


text:
# General pairing interactions and pair truncation approximations for fermions in a single-𝑗 shell

## I Introduction

Since pairing has proven to be important in atomic, nuclear, and condensed matter physics, pair truncation approximations to many body wavefunctions have been extensively studied. The first example is the seniority scheme, introduced by Racah , for the classification of states in atomic spectra and later applied extensively in nuclear physics, where $`S`$ pairs with angular momentum zero are related with a strong and attractive monopole pairing interaction. The second example is the interacting boson model (IBM) introduced by Arima and Iachello , where the low-lying excitations of complex even-even nuclei are described successfully by $`s`$ bosons which correspond to correlated $`S`$ nucleon pairs with angular momentum zero and $`d`$ bosons which correspond to $`D`$ nucleon pairs with angular momentum two. Again, the success of the IBM in nuclear physics partly comes from the validity of the pairing plus quadrupole-quadrupole force for effective interactions between valence nucleons. Monopole and quadrupole pairing are important as well in low and high temperature superconductiviy in materials .

In this paper we investigate the general pair truncation approximation for fermions in a single-$`j`$ shell. The examples explored may provide a clue as how to classify the states which come from the diagonalization of an attractive pairing interaction for which two fermions are coupled to an angular momentum $`J`$,

$`H_J={\displaystyle \underset{M=J}{\overset{J}{}}}A_M^JA_M^J`$ (1)
$`A_M^J={\displaystyle \frac{1}{\sqrt{2}}}\left[a_j^{}\times a_j^{}\right]^J,A_M^J=(1)^M{\displaystyle \frac{1}{\sqrt{2}}}\left[\stackrel{~}{a}_j\times \stackrel{~}{a}_j\right]_M^J.`$ (2)
where $`[]_M^J`$ means coupled to angular momentum $`J`$ and projection $`M`$. Most of examples pursued in this paper have $`n`$ = 4, where $`n`$ is the nucleon number.

We first point out in this paper that the low-lying eigenstates of (2) can be approximated by wavefunctions with pairs with angular momentum $`J`$ only. We shall next show that a large array of eigenvalues of four nucleons in a single-$`j`$ shell are asymptotic integers when $`JJ_{max}=2j1`$ and the total angular momentum, $`I`$, is not very close to $`I_{max}=4j6`$, and that this phenomenon originates from validity of the pair truncation scheme and special features of coupling coefficients. We shall finally prove that the pair Hamiltonian (2) has exactly one and only one non-zero eigenvalue for four fermion eigenstates with total angular momentum zero, $`I=0`$. This sheds light on the problem of angular momentum zero ground state dominance in many-body systems interacting by random interactions .

## II Comparison of the Pair Approximation to the Shell Model

Fig. 1 compares the exact ground state angular momentum $`I`$ for four nucleons in a single-$`j`$ shell interacting by the attractive pair Hamiltonian $`H_J`$ with the angular momentum $`I`$ of the ground state in the truncated space of pairs coupled to angular momenutm $`J`$ only. We have examined all the cases up to $`J=20`$ and $`j31/2`$ but here we show only two typical examples with $`J=6`$ and 14. In the case of $`J=0`$, the seniority scheme ($`S`$-pair approximation) produces the exact ground state. In case of $`J=2`$ and $`n=4`$, a $`D`$-pair approximation is found to be always very good. When $`J>2`$, the $`J`$-pair approximation of low-lying states is not perfect but always very reasonable. In Fig. 1, most of ground state angular momenta are correctly reproduced by the $`J`$-pair truncation. In Fig. 1b), there are two exceptions in which the ground state is not correctly given by two $`J=14`$ pairs: $`2j=19`$ and 25.

Even in those cases where the ground state angular momenta are not correctly predicted by the $`J`$-pair truncation, the low-lying state energies are reasonably reproduced (including the binding energies). As a “bad” example in which the ground state angular momentum is not reproduced, we show in Fig. 2 the case of four nucleons in a $`j=25/2`$ shell with $`J=14`$. The calculated levels using two $`J=14`$ pairs are shown in the first column. The next two columns are the shell model states obtained by diagonalizing the Hamiltonian in the full space. The states in the second column are the shell model states correponding to the pair truncation states in the first column. All the levels below $`0^+`$ are included. One sees that the lowest four states $`2_1^+`$, $`6_1^+`$, $`12_1^+`$, and $`10_1^+`$ are well approximated by two $`J=14`$ pairs, although the ground state angular momentum is not correctly given by this $`J=14`$ pair approximation. The $`0_1^+`$ state is always precisely reproduced (as we shall see there is only one non-zero eigenvalue for $`I=0`$ states and that eigenstate is constructed by $`J`$ pairs). The states with odd $`I`$ are always outside the pair truncation space but their energies are quite high in all cases that we have checked. In Fig. 2, the $`9_1^+`$ is the lowest state with odd value of $`I`$. Below the $`9_1^+`$ there are 10 states with even values of $`I`$ and most of which are nicely described by the $`J=14`$ pair approximation. The angular momentum of those states for which excited energies are not well reproduced by two $`J`$ pairs are labled by italic font in Fig. 2.

For $`J=J_{max}=2j1`$, the $`I=I_{max}=4j6`$ or $`I=I_{max}2`$ states are always the lowest. These two states may be constructed by pairs with angular momentum either $`J_{max}`$ or $`J_{max}2`$. However, pairs with angular momentum $`J_{max}2`$ do not present a good classification for other $`I`$ states while those with angular momentum $`J_{max}`$ do.

For n=3, the $`J`$-pair truncation describes the low-lying states precisely. For $`n`$=5 and 6 cases that we have examined, the low-lying states are reasonably approximated by the $`J`$-pair truncation. We note without details that bosons with spin $`l`$ exhibit a similar situation. It would be interesting to know the situation in more complicated systems.

## III Integer Eigenvalues

We next report a very interesting regularity in the spectrum of the Hamiltonian (2) with $`J=J_{max}=2j1`$. The eigenvalues of most states with $`I2j3`$ are found to be very close to integers corresponding to the number of pairs with angular momentum $`J_{max}`$ except for a very few eigenvalues. Taking four nucleons in a $`j=31/2`$ shell as an example, the diagonalization of $`H_J`$ ($`J=30`$) gives “integer” eigenvalues for low $`I`$ states$``$ all eigenvalues are 0, $`1`$, and $`2`$ to within a precision of 0.01 for all states with $`I<22`$. For states with $`22I52`$, these three “integer” eigenvalues continue to be valid except that 7 eigenvalues which are not close to 0, $`1`$ or $`2`$ come in. These “non-integer” eigenvalues are very stable in magnitude for states with $`22I52`$. The states with $`I53`$ are one dimensional, so the corresponding eigenvalues (saturate quickly with $`j`$) may be analytically derived . These “integer” eigenvalues are best seen in case of $`J=J_{max}`$ and becomes less dominant for smaller $`J`$ and the same single-$`j`$ shell.

To understand the validity of the pair approximation and the occurrence of these “integer” eigenvalues we consider the pair basis of four nucleons

$$|j^4[J_1J_2]I,M=\frac{1}{\sqrt{N_{J_1J_2;J_1J_2}^{(I)}}}\left(A^{J_1}\times A^{J_2}\right)_M^{(I)}|0,$$
(3)
where $`N_{J_1J_2;J_1J_2}^{(I)}`$ is the diagonal matrtix element of the normalization matrix

$`N_{J_1^{}J_2^{};J_1J_2}^{(I)}=0|\left(A^{J_1^{}}\times A^{J_2^{}}\right)_M^{(I)}\left(A^{J_1}\times A^{J_2}\right)_M^{(I)}|0`$ (4)
$`=\delta _{J_1,J_1^{}}\delta _{J_2,J_2^{}}+()^I\delta _{J_1,J_2^{}}\delta _{J_2,J_1^{}}4\widehat{J}_1\widehat{J}_2\widehat{J_1^{}}\widehat{J_2^{}}\left\{\begin{array}{ccc}j& j& J_1\\ j& j& J_2\\ J_1^{}& J_2^{}& I.\end{array}\right\}.`$ (8)
In general this basis is overcomplete and the normalization matrix may have zero eigenvalues for a given $`I`$.

The matrix elements of $`H_J`$ are

$$j^4[J_1^{}J_2^{}]I,M|H_J|j^4[J_1J_2]I,M=\frac{1}{\sqrt{N_{J_1J_2;J_1J_2}^{(I)}N_{J_1^{}J_2^{};J_1^{}J_2^{}}^{(I)}}}\underset{J^{}=\mathrm{even}}{}N_{J_1J_2;JJ^{}}^{(I)}N_{J_1^{}J_2^{};JJ^{}}^{(I)}$$
(9)
where $`\widehat{J}_1`$ is a short hand notation of $`\sqrt{2J_1+1}`$. There are two terms in $`N_{J_1J_2;JJ^{}}^{(I)}`$: the second term is a nine-$`j`$ coefficient which is usually much less than unity in magnitude, in particular when $`J`$ is large and $`I`$ not large (refer to Appendix A). Neglecting this nine-$`j`$ symbol, the allowed states are $`|j^4[J_1J_2]I,M\left(A^{J_1}\times A^{J_2}\right)_M^{(I)}|0`$ for $`J_2<J_1`$, and $`|j^4[J_1J_1]I,M\frac{1}{\sqrt{2}}\left(A^{J_1}\times A^{J_1}\right)_M^{(I)}|0,I`$ even only, and the Hamiltonian matrix becomes

$$j^4[J_1^{}J_2^{}]I,M|H_J|j^4[J_1J_2]I,M\frac{(\delta _{J_1,J_1^{}}\delta _{J_2,J_2^{}}+(1)^I\delta _{J_1,J_2^{}}\delta _{J_2,J_1^{}})}{\sqrt{1+(1)^I\delta _{J_1,J_2}}\sqrt{1+(1)^I\delta _{J_1^{},J_2^{}}}}(\delta _{J_1,J}+\delta _{J_2,J}).$$
(10)
First of all we see that the matrix is diagonal which validates the pair approximation. Secondly the eigenvalues are either 0, $`1`$, or $`2`$ with their corresponding wavefunctions having 0, 1 or 2 pairs with angular momentum $`J`$, respectively. Therefore, the integer eigenvalues of many $`I`$ states originate from both the special properties of these nine-$`j`$ symbols and the validity of $`J`$-pair truncation.

From the $`J`$-pair coupling scheme, the number of $`|j^4[J_1J_2]I,M`$ with $`J_1=J_2=2j1`$ is $`\frac{1+()^I}{2}`$, and the number of $`|j^4[J_1J_2]I,M`$ with $`J_1=2j1`$, $`J_2<J_1`$ and $`I<2j1`$ is $`[\frac{I}{2}]`$ (the largest integer not exceeding $`\frac{I}{2}`$). According to the above discussion, the number of states with eigenvalues $`2`$ is $`\frac{1+()^I}{2}`$ and the number of states with eigenvalues $`1`$ is $`[\frac{I}{2}]`$. This is confirmed in all cases with $`I<2j8`$. Eigenvalues not close to integers arise in states with $`2j8I4j12`$. These “non-integer” eigenvalues are found to be almost the same for the $`2j8I4j12`$ states; an understanding of this regularity is in progress.

From a more general expression of Eq. (9), say, Eq. (5.8) of , we expect that the “integer” eigenvalues appear not only in even systems but also appear in odd-A systems. According to our numerical results, the pattern of “integer” eigenvalues also appear in the states with small $`I`$ for $`j11/2`$ and $`n=3`$, and for $`j23/2`$ and $`n=5`$ etc. For cases with $`n=3`$, a similar proof is readily obtained in terms of 6-$`j`$ symbols. An explicit proof for more nucleons will be quite complicated.

## IV Exact Results for Angular Momentum Zero States

We now come to the last point of this paper by pointing out that there is only one non-zero eigenvalue for $`I=0`$ and $`n=4`$. We define a new basis for the $`I=0`$ states for the pair Hamiltonian $`H_J`$

$$|j^4J_1=|j^4[J_1J_1]I=M=0\frac{N_{J_1J_1;JJ}^{(0)}}{\sqrt{N_{JJ;JJ}^{(0)}N_{J_1J_1;J_1J_1}^{(0)}}}|j^4[JJ]I=M=0,(J_1J),$$
(11)
$$|j^4J=|j^4[JJ]I=M=0.$$
(12)
The unnormalized states $`|j^4J_1,J_1J`$, are orthogonal with respect to $`\mathrm{𝑜𝑛𝑙𝑦}`$ $`|j^4[JJ]I=M=0`$. From (9) one has

$$j^4[J_1^{}J_1^{}]0,0|H_J|j^4[J_1J_1]0,0=\frac{N_{J_1J_1;JJ}^{(0)}N_{J_1J_1;JJ}^{(0)}}{\sqrt{N_{J_1^{}J_1^{};JJ}^{(0)}N_{J_1J_1;JJ}^{(0)}}},$$
(13)
where $`J_1^{},J_1=0,2,\mathrm{},2j1`$. Using this formula, one easily confirms that all matrix elements of the Hamiltonian in the basis (11), $`j^4J_1^{}|H_J|j^4J_1,J_1^{}J`$, are zero. Therefore all the eigenvalues for $`n=4`$ and $`I=0`$ are zero except for the single state with both pairs having angular momentum $`J`$ and its eigenvalue is $`E_0^{J(j)}=N_{JJ;JJ}^{(0)}`$.

Since $`H_J`$ is a negative definite operator, its eigenvalues will be negative or zero. From above we see that, for $`I`$ = 0, there is only one state with a non-zero eigenvalue, $`E_0^{J(j)}`$. Therefore one expects this eigenvalue to be the lowest in the spectrum because the eigenvalues of $`I0`$ states are more or less scattered in many states, generally speaking. Thus the probability that the $`I=0^+`$ is the lowest state of four fermions in a single-$`j`$ shell in the presence of random two-body interactions is expected to be larger than the probability for all other angular momentum $`I`$, according to the empirical rule of Ref. . For $`j31/2`$ there are only two exceptions , $`j=7/2`$ and $`13/2`$.

The sum rule of diagonal matrix elements  gives $`_JE_0^{J(j)}=\frac{1}{2}n(n1)D_0^{(j)}=6D_0(j)`$, where $`D_0^{(j)}`$ is the number of $`I=0`$ states and here $`n=4`$. For $`n`$ = 4 the number of states is $`D_0^{(j)}=[(2j+3)/6]`$  which gives 1, 1, 1, 2, 2, 2, 3, 3, 3, $`\mathrm{}`$ for $`2j`$=3, 5, 7, 9, 11, 13, 15, $`\mathrm{}`$, etc., regularly. Thus the staggering in the number of states is expected to be reflected in the staggering of the energy which was pointed out in  but without an explanation.

## V Summary

In this paper we have shown that an attractive $`J`$-th pairing interactions favors pairs with angular momentum $`J`$ in low-lying states of fermions in a single-$`j`$ shell. Therefore, one may use pairs with angular momentum $`J`$ as building blocks of wavefunctions of low-lying states. This is in contrast to repulsive pair interactions, used, for example, in the fractional quantum Hall effect , for which the pair truncation approximation is not valid.

In addtion we discovered that the eigenvalues of states with low angular momentum $`I`$ for pair Hamiltonians with $`J`$ larger than $`j`$ are approximately integers. We explain the origin of this fact for four nucleons in a single-$`j`$ shell from the validity of pair truncation and special properties of nine-$`j`$ symbols. We point out without details that the same holds for $`n=3`$ and $`j11/2`$ and for $`n=5`$ ($`j23/2`$).

Finally we pointed out that there is only one non-zero eigenvalue for $`I=0`$ and $`n=4`$. This result, together with the empirical rule of ref. , provides a simple argument for the large probability of angular momentum zero states to be the lowest in energy for $`n=4`$. Also for the first time, we demonstrated that the staggering of this probability for $`n=4`$ vs. $`j`$ is directly related to the staggering of number of $`I=0`$ states vs. $`j`$.

## VI Acknowledgements

One of the authors (ZYM) is grateful to the Japan Society for the Promotion of Science (contract No. P01021) for supporting his work, and he also thanks Los Alamos National Laboratory for the warm hospitality extended to him. This work was partly supported by U.S. Department of Energy under contract W-7405-ENG-36.

## VII Appendix A: Some nine-$`j`$ symbol formulae

These formulae are obtained in the following two steps: First, rewrite the nine-$`j`$ in terms of six-$`j`$ symbols, i.e.,

$`\left\{\begin{array}{ccc}j& j& r_1\\ j& j& r_2\\ s_1& s_2& I\end{array}\right\}={\displaystyle \underset{t}{}}()^{2t}(2t+1)\left\{\begin{array}{ccc}j& j& r_1\\ r_1& I& t\end{array}\right\}\left\{\begin{array}{ccc}j& j& r_2\\ j& t& s_2\end{array}\right\}\left\{\begin{array}{ccc}s_1& s_2& I\\ t& j& j\end{array}\right\},`$ (23)
and second, make use of the analytical formulas of six-$`j`$. Through these examples (though we are unable to get a universal formulae), one sees that the nine-$`j`$ symbols in Eq. (3) are much less than unity and may be neglected in (9) when $`I`$ is not very large.

$`\left\{\begin{array}{ccc}j& j& 2j1\\ j& j& 2j1\\ 2j1& 2j1& 0\end{array}\right\}={\displaystyle \frac{j(4j3)\left[(2j1)!\right]^2}{(4j1)(4j1)!}};`$ (27)
$`\left\{\begin{array}{ccc}j& j& 2j1\\ j& j& 2j1\\ 2j1& 2j1& 2\end{array}\right\}={\displaystyle \frac{j(8j^26j3)\left[(2j1)!\right]^2}{(4j1)(4j3)(4j1)!}};`$ (31)
$`\left\{\begin{array}{ccc}j& j& 2j1\\ j& j& 2j1\\ 2j1& 2j1& 4\end{array}\right\}={\displaystyle \frac{3j(2j+1)(4j^23j5)\left[(2j1)!\right]^2}{(4j1)(4j3)(4j5)(4j1)!}};`$ (35)
$`\left\{\begin{array}{ccc}j& j& 2j1\\ j& j& 2j1\\ 2j1& 2j1& 6\end{array}\right\}={\displaystyle \frac{5j(j+1)(2j+1)(8j^26j21)\left[(2j1)!\right]^2}{(4j1)(4j3)(4j5)(4j7)(4j1)!}};`$ (39)
$`\left\{\begin{array}{ccc}j& j& 2j1\\ j& j& 2j3\\ 2j1& 2j3& 2\end{array}\right\}={\displaystyle \frac{j\left(36+j(4j9)\left(19+2j(9+4j)\right)\right)\left[(2j1)!\right]^2}{3(4j3)(4j5)(4j1)!}};`$ (43)
$`\left\{\begin{array}{ccc}j& j& 2j1\\ j& j& 2j3\\ 2j1& 2j3& 3\end{array}\right\}={\displaystyle \frac{j(4j3)\left(1+2j(4j9)\right)\left[(2j1)!\right]^2}{6(4j5)(4j1)!}}.`$ (47)
From the above formulae one sees that these nine-$`j`$ symbols are proportional to $`\frac{\left[(2j1)!\right]^2}{(4j1)!}`$, and are very close to zero when $`j`$ becomes considerably large. For example, the absolute values of these nine-$`j`$ symbols are less than $`10^{15}`$ when $`j=31/2`$.

## VIII Appendix B: A New Sum Rule for a six-$`j`$ symbol

From Eq. (3) one obtains

$`N_{JJ;JJ}^{(0)}=2+4(2J+1)\left\{\begin{array}{ccc}j& j& J\\ j& j& J\end{array}\right\}.`$ (50)
Since $`N_{JJ;JJ}^{(0)}`$ is also the unique eigenvalue of the $`I=0`$ eigenstate of $`H_J`$, one has a sum rule that

$`{\displaystyle \underset{\mathrm{even}J}{}}N_{JJ;JJ}^{(0)}={\displaystyle \frac{1}{2}}n(n1)D_0^{(j)},`$ (51)
where $`n=4`$, $`D_0^{(j)}`$ = $`\left[\frac{2j+3}{6}\right]`$ . One finally obtains that

$`{\displaystyle \underset{\mathrm{even}J}{}}(2J+1)\left\{\begin{array}{ccc}j& j& J\\ j& j& J\end{array}\right\}=\{\begin{array}{cc}\frac{1}{2},& \mathrm{if}2j=3k,\hfill \\ 0,& \mathrm{if}2j+2=3k,\hfill \\ \frac{1}{2},& \mathrm{if}2j2=3k.\hfill \end{array}`$ (57)