text:
# Magnus Expansion and Three-Neutrino Oscillations in Matter 1footnote 11footnote 1Presented at Mexican School of Astrophysics (EMA), Guanajuato, Mexico, July 31 - August 7, 2002. Final version to appear in the Proceedings of IX Mexican Workshop on Particles and Fields “Physics Beyond the Standard Model”, Colima Col. México, November 17-22, 2003.

## I Introduction

The need to introduce the three generations of neutrinos into a computation of the transition probabilities of these particles while traversing a medium has been recognised long ago need  to try to accomodate the observations of different experiments studying neutrino oscillations. Analytical treatments of three neutrino oscillations in matter with varying density in the three generation scheme have been studied in the past tres  aiming to deduce expressions for the oscillation probabilities of one type of neutrino into another. It has been shown, first in the case of two generations petcov\_1 ; toshev , and later in the case of three generations osland  that the differential equation describing the evolution of the neutrino state in an exponentialy varying density profile could be solved analytically in terms of confluent hypergeometric functions. Corrections to the mixing parameters in matter calculated as series expansions have been performed by Freund freund , and a different perturbative analysis has been done by Narayan in narayan . Global analyses of the recent experimental dathave been extensively studied both, in the two and three generations cases concha . In a previous paper by D’Olivo and Oteo dolivo\_oteo  an approximate expression to the evolution operator using the Magnus expansion was found, only for the non-adiabatic regime in the exponentially varying density profile. In this paper we present a complete semi-analytical computation of the evolution operator for neutrinos in the same density profile, using the Magnus expansion approximation which will work propperly in the case of adiabatic and non-adiabatic evolution of the neutrino state. This paper describes some important results of the M. C. Thesis work presented by Luis G. Cabral-Rosetti in tesis . Let $`H`$ denote the Hamiltonian of a quantum systen and $`U=U(t,t_0)`$ the time evolution operator satisfying the Schrödinger equation

$$i\mathrm{}\frac{}{t}U=HU,U(t_0,t_0)=I.$$
(1)
When $`H`$ is independent of time, or more generally when $`[_{t_0}^t𝑑t^{}H(t^{}),H(t)]=0`$, the solution of Eq.(1) is formally $`U=\mathrm{Exp}[i/\mathrm{}_{t_0}^t𝑑t^{}H(t^{})]`$. Then it is natural to ask wether a solution of the form $`U=\mathrm{Exp}\mathrm{\Omega }`$ would always be possible. A method for finding such a true exponential solution (without time ordering) is supplied by the Magnus Expansion (ME) magnus . The magnus operator $`\mathrm{\Omega }=\mathrm{ln}U`$ satisfies a differential equation which in turn is solved through a series expansion: $`\mathrm{\Omega }=_{n=1}^{\mathrm{}}\mathrm{\Omega }_n`$, where each term $`\mathrm{\Omega }_n`$ is of order $`\mathrm{}^n`$. The first two contributions are explicitly given by

$$\mathrm{\Omega }_1=\frac{i}{\mathrm{}}_{t_0}^t𝑑t_1H(t_1),$$

$$\mathrm{\Omega }_2=\frac{i}{2\mathrm{}^2}_{t_0}^t𝑑t_1_{t_0}^{t_1}𝑑t_2[H(t_1),H(t_2)].$$
(2)
Recursive methods to obtain the succesive terms have been extensively worked out in the literature burum . Because of the anti-Hermitian character of every $`\mathrm{\Omega }_n`$, each approximate time-evolution operator obtained as $`UU_k=\mathrm{Exp}(_{n=1}^k\mathrm{\Omega }_n)`$ will be unitary. Here we use the first Magnus approximant to obtain (approximate) analytical solutions to the problem of 3 neutrinos oscillations in a medium with varying density as the Sun.

## II Survival Probability for electron Neutrinos

From the dependency of $`\mathrm{sin}\mathrm{\hspace{0.33em}2}\theta _{12}^m(t)`$, and $`\mathrm{sin}\mathrm{\hspace{0.33em}2}\theta _{13}^m(t)`$ on $`V(t)`$ (see Ref. tesis  and nosotros  for all details), we can distinguish five interesting regions where the oscillations behave differently.

(1) Region of Extreemely Low Density
If $`0V(t)\mathrm{\Delta }_{21}`$, both mixing angles, $`\theta _{12}^m(t)`$, and $`\theta _{13}^m(t)`$, are close to their values in vacuum (see Fig.(1)), this is $`\theta _{12}^m(t)\theta _{12}`$ ; $`\theta _{13}^m(t)\theta _{13}`$, giving

$$𝒫_\alpha (t)c_{12}^2c_{13}^2,\mathrm{and}𝒫_\gamma (t)s_{13}^2,$$
(3)
leading us to the well known result of vacuum oscillations of three neutrinos

$$<𝒫(\nu _e\nu _e)>=c_{12}^4c_{13}^4+s_{12}^4c_{13}^4+s_{13}^4,$$
(4)
(2) Low Density Resonance Region
For $`V(t)\mathrm{\Delta }_{21}`$, $`\theta _{13}^m(t)`$ is still close to its vacuum value ($`\theta _{13}^m(t)\theta _{13}`$), while $`\theta _{12}^m(t)`$ is at the low density resonance, i.e. $`\theta _{12}^m(t)\frac{\pi }{4}`$, making

$`𝒫_\alpha (t)`$ $``$ $`\mathrm{cos}^2\mathrm{\Pi }(t)\mathrm{cos}^2\theta _{12}^m(t_0)\mathrm{cos}^2\theta _{13}^m(t_0)`$ (5)
$`+\mathrm{sin}^2\mathrm{\Pi }(t)\mathrm{sin}^2\theta _{12}^m(t_0)\mathrm{cos}^2\theta _{13}^m(t_0),`$
and

$$𝒫_\gamma (t)\mathrm{sin}^2\theta _{13}^m(t_0),$$
(6)
which leads to the expression

$$<𝒫_{(\nu _e\nu _e)}>=\frac{1}{2}+\frac{1}{2}\left(12𝒫_c^l\right)\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{12}^m(t_0)\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{12},$$
(7)
where

$$\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{12}^m(t_0)=\frac{V_lV(t_0)}{\sqrt{(V_lV(t_0))^2+B_l^2}},$$
(8)
and

$$𝒫_c^l=\mathrm{sin}^2\left[\left(\theta _{12}^m(t_0)\theta _{12}^m(T)\right)\text{Exp}(\kappa _l)\right].$$
(9)
The quantity $`𝒫_c^l`$ represents the tramsition probability between the states $`\nu _{2m}`$ y $`\nu _{1m}`$, and has been derived by D’Olivo D'Olivo3  in the case of oscillations between two neutrino species. The adiabatic result can be recovered by seting $`𝒫_c^l=0`$. For $`\kappa _l1`$, $`𝒫_c^l`$ is exponentially supressed as expected in the asimptotic regime. On the other hand, for $`\kappa _l<1`$ there are significant corrections to the adiabatic approximation which reduce the effect of the resonant transition. Eq.(7) was first obtained by Parke Parke4  for two species using the Landau-Zener approximation for the crossing probability: $`𝒫_c^l=Exp\left(\frac{\pi }{2}\kappa _l\right)`$. In the extreeme non adiabatic case, $`\mathrm{\Delta }_{21}0`$, making $`\kappa _l0`$, and from Eq. (8) it follows that $`\theta _{12}^m(t_0)\frac{\pi }{2}`$, making $`𝒫_c^l=\mathrm{cos}^2\theta _{12}`$. Introducing this value of $`𝒫_c^l`$ in Eq.(7) we recover the vacuum result for two neutrinos:

$$<𝒫_{(\nu _e\nu _e)}>=1\frac{1}{2}\mathrm{sin}^2\mathrm{\hspace{0.33em}2}\theta _{12},$$
(10)
This should be contrasted with the result

$$<𝒫(\nu _e\nu _e)>=\mathrm{cos}^2\theta _{12},$$
(11)
predicted by the Landau-Zener formula (and, in general, by any result derived from the Dykhne’s formula), which deviates from the correct limit given by Eq.(10), when $`\theta _{12}`$ is large. The correct value for the extreeme nonadiabatic case has been derived before, under the assumption that the transition between the adiabatic eigenstates occurs instantaneously at the time $`t=t_l`$ Kim2 .

However, in this case the corresponding result for $`𝒫_c^l`$ approaches $`\mathrm{cos}^2\theta _{12}`$ as $`(\frac{\delta m^2}{E})^2`$ instead of linearly, as in Eq. (9). Non adiabatic effects start to become important in this region when $`\kappa _l`$ is comparable to 1, whenever the neutrinos cross the low density resonance. If $`E<\frac{\delta m_{21}^2}{2V(t_0)}\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{12}`$, $`𝒫_c^l=0`$, and the propagation of neutrinos will be adiabatic for $`\kappa _l<1`$. For this reason, the asimptotic exponential expression ($`\mathrm{}0`$) for $`𝒫_c^l`$, must be modified by hand to consider this situation. An effective way of implementing such modification is to multiply $`𝒫_c^l`$ by a step function $`\mathrm{\Theta }(V(t_0)\mathrm{\Delta }_{21}cos\mathrm{\hspace{0.33em}2}\theta _{12})`$, in such a way that the transition probability vanishes if a neutrino is produced after the resonance. It is worthnoting that such modification is not necessary with the Magnus result Eq. (9) given that, as a function of $`\frac{\delta m_{21}^2}{2E}`$, the difference $`\theta _{12}^m(t_0)\theta _{12}^m(T)`$ behaves as a continuous step.

(3) Intermediate Density Region
For $`\mathrm{\Delta }_{21}V(t)\mathrm{\Delta }_{31}`$, the mixing angles in matter are $`\theta _{12}^m(t)\frac{\pi }{2}`$, and $`\theta _{13}^m(t)\theta _{13}`$, making

$$𝒫_\alpha (t)c_{12}^2c_{13}^2,$$
(12)
and

$$𝒫_\gamma (t)s_{13}^2,$$
(13)
giving the result shown in Eq.(4), provided there exists a clear separation between the resonance regions: $`\mathrm{\Delta }_{31}\mathrm{\Delta }_{21}`$.

(4) High Density Resonance Region
When $`V(t)\mathrm{\Delta }_{31}`$, the mixing angle $`\theta _{12}^m(t)`$ is approximately equal to $`\frac{\pi }{2}`$, while $`\theta _{13}^m(t)`$ is at its high density resonance value of $`\theta _{13}^m(t)\frac{\pi }{4}`$. In this case

$$𝒫_\alpha (t)\mathrm{cos}^2\theta _{13}^m(t_0)\mathrm{cos}^2\theta _{12}^m(t_0),$$
(14)
and

$$𝒫_\gamma (t)cos^2\mathrm{\Pi }(t)sin^2\theta _{13}^m(t_0),$$
(15)
leading to

$$<𝒫_{(\nu _e\nu _e)}>=\frac{1}{2}+\frac{1}{2}\left(12𝒫_c^h\right)\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{13}^m(t_0)\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{13}$$
(16)
where

$$\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{13}^m(t_0)=\frac{V_hV(t_0)}{\sqrt{(V(t_0)V_h)^2+B_h^2}},$$
(17)
and

$$𝒫_c^h=\mathrm{sin}^2\left[\left(\theta _{13}^m(t_0)\theta _{13}^m(T)\right)\text{Exp}(\kappa _h)\right],$$
(18)
which again implies the two neutrino oscillations result. The crossing probability between the instantaneous eigenstates $`\nu _{3m}`$ and $`\nu _{2m}`$ given by $`𝒫_c^h`$ is perfectly analogue to that studied in the low density resonance region. We can recover the adiabatic case if we set $`𝒫_c^h=0`$, and for $`\kappa _h1`$, $`𝒫_c^h`$ is exponentially supressed. The extreeme nonadiabatic case requires $`\mathrm{\Delta }_{31}0`$, and from Eq. (8) we have $`\theta _{13}^m(t_0)\frac{\pi }{2}`$ leading to $`𝒫_c^h=\mathrm{cos}^2\theta _{13}`$. All these observations will lead us to an expression of the form of Eq.(10) with $`\theta _{12}`$ replaced by $`\theta _{13}`$. Non adiabatic effects become important when $`\kappa _h1`$ provided the neutrino crosses the High density resonance. If $`E<\frac{\delta m_{31}^2}{2V(t_0)}(\mathrm{cos}\mathrm{\hspace{0.33em}2}\theta _{13}R^{\mathrm{\hspace{0.33em}1}}\mathrm{sin}^2\theta _{12}\mathrm{cos}\mathrm{\hspace{0.17em}2}\theta _{13})`$, with $`R=\mathrm{\Delta }_{31}/\mathrm{\Delta }_{21}`$, no transitions between the instantaneous eigenstates can occur ($`𝒫_c^h=0`$), and the propagation is adiabatic for $`\kappa _h<1`$. Again, the difference $`\theta _{13}^m(t_0)\theta _{13}^m(T)`$ behaves as a continuous step giving $`𝒫_c^h`$ the appropiate behaviour.

(5) Extreemely High Density Region
For $`V(t)\mathrm{\Delta }_{31}`$ the $`(\nu _e)`$ oscillations are strongly supressed due to the fact that $`\theta _{12}^m(t)\frac{\pi }{2}`$, and $`\theta _{13}^m(t)\frac{\pi }{2}`$, making

$$𝒫_\alpha (t)c_{12}^2c_{13}^2$$
(19)
and

$$𝒫_\gamma (t)s_{13}^2,$$
(20)
giving the result of Eq.(4)

$$<𝒫(\nu _e\nu _e)>=c_{12}^4c_{13}^4+s_{12}^4c_{13}^4+s_{13}^4,$$
(21)
Plots of the $`\nu _e`$ survival probability $`<P_{(\nu _e\nu _e)}>`$ as a function of $`\frac{\delta m^2}{2E}`$ are shown in Figs.2–3, for different values of $`\mathrm{sin}^22\theta _{12}`$, $`\mathrm{sin}^22\theta _{13}`$, and $`R`$. We use the exponential profile bahcall10  $`N_e(r)=245Exp(10.54r/R_{})N_{Avo}cm^{3}`$, where $`N_{Avo}`$ is the Avogadro’s number, r is the radial distance measured from the center of the Sun, and $`R_{}`$ is the solar radius ($`R_{}=6.96\times 10^5Km`$). Except for those regions close to the center or the surface, this is a good approximation of the electron density in the Sun. We further assume that the $`\nu _e`$ are produced at $`r_0=0.08635R_{}`$, where $`N_e`$ is the central electron density predicted by the Standard Solar Model (SSM). In Fig.4 we show the allowed region in the $`\mathrm{\Delta }_{21}`$ vs $`\mathrm{sin}^22\theta _{12}`$ plane obtained by simply plotting the points in the plane wit survival probabilities lying within the $`1\sigma `$ range of values extracted from recent joint analyses of the results of SNO and Super Kamiokande implications1 . To produce this region we used the estimated value $`P_e=(0.34\pm 0.05)`$ for the survival probability, and an average neutrino energy of 5 MeV. We computed similar regions using the survival probabilities estimated by the experiments SAGE, GALEX and HOMESTAKE, taking the average energies for this experiments as in ohlsson , and looked for an overlap of these regions in the parameter plane (see Fig. 5). No intention to give statistical significance to this region exists, but only it is shown that our result it is consistent with those achieved by rigurous analyses.

###### Acknowledgements.

This work has been supported in part by DGAPA-UNAM under Grant PAPIIT proyect No. IN109001 and in part by CoNaCyT (México) under Grant No. I30307-E.


format: text

text:
# UV Continuum, Physical Conditions and Filling Factor in Active Galactic Nuclei

## 1 Introduction

The intrinsic spectral energy distribution (SED) of active galactic nuclei (AGN) extends from radio up to $`\gamma `$ rays with observations available in a wide wavelength range. However, due to Galactic and intrinsic absorption, the SED above the Lyman limit cannot be observed. Since the main mechanism powering the emission-line spectrum of AGN is photoionization, it is important to know the shape of the UV radiation. As suggested by the observed continuum, it is usually assumed that the SED from ultraviolet to soft X-rays is basically a power law with at least two other components corresponding to the so-called “big blue bump”, usually represented by a blackbody radiation with temperature between 1.5 to 3.0 $`\times `$ 105 K (see, for instance Oliva et al. 1994, Prieto & Viegas 2000), and to the soft X-ray excess, frequently described as a thermal component (see, for instance, Terashima et al. 2002). In addition, the SED is also important to test the emission by supermassive black holes and to discriminate among different models (Narayan, Kato & Homna 1997), as well as to help finding the signature of shocks in the narrow emission-line region (Contini, Viegas e Prieto 2002, and references therein).

Many observations of Seyfert galaxies point to the existence of three kinematically distinct emission regions: (a) an unresolved broad line region (BLR), with typical line widths of $``$ 5000 km s-1, densities above 108 cm-3 and ionization parameter U $``$ 0.03-1; (b) a narrow line region (NLR), with lower densities ($``$ 107 cm-3) and lower velocities ($``$ 500 km s-1) gas clouds with U in the range 0.03 to 0.1 embedded in a more diluted gas; (c) an extended emission-line region (EELR) with typical line widths of less than 50 km s-1, densities below 103 cm-3, and U smaller than 0.005 (Schulz and Komossa 1993).

The emission-line spectra of AGN show a rich variety of optical and UV lines, including those from highly ionized species (coronal lines). The observed line ratios point to photoionization as the dominant mechanism powering these objects. Recently, however, near infrared observations of narrow-line Sy 1 galaxies tend to favor shocks as the mechanism producing coronal lines (Rodriguez-Ardila et al. 2002). Nevertheless, photoionization being the main mechanism contributing to the emission-line spectra of AGN, photoionization models can be used to constrain the SED by fitting the observed lines.

Photoionization codes presently available were developed in the 70’s. They are one-dimensional (1-D) codes that usually assume a homogeneous gas distribution with a well-defined symmetry (plan-parallel or spherically symmetric). The standard input parameters are the shape and intensity of the ionizing radiation spectrum, the density of the gas, and the elemental abundances. Therefore, 1-D models are extremely limited to simulate clumpy regions like HII regions or the narrow line region of AGN.

In order to simulate an inhomogeneous region, two types of models are usually proposed: (a) those assuming a homogeneous distribution of gas and a filling factor less than unity (see, for instance Osterbrock 1989); (b) a multi-cloud model generated by a weighted average emission-line spectrum from single-cloud simulations (with $`ϵ`$ = 1), as firstly proposed by Viegas-Aldrovandi & Gruenwald (1988), and more recently used by other authors (Baldwin et al. 1995, Binette, Wilson & Storchi-Bergmann 1996, Fergunson et al. 1997, Komossa & Shulz 1997, Kraemer & Crenshaw 2000). Multi-cloud models of the NLR, based on sigle-cloud models which take into account the coupled effect of photoionization and shocks, have also been used to analyze the NLR of AGN (Contini, Prieto & Viegas 1998a, 1998b, Contini & Viegas 2001, Contini et al. 2002).

The idea of a filling factor was first introduced by Osterbrock & Flatter (1959) in order to explain the radio brightness temperature of HII regions lower than expected, indicating that the gas has a clumpy distribution. The filling factor $`ϵ`$ is the fraction of the total volume occupied by the clumps or condensations. Since 1-D photoionization codes are unable to properly deal with an inhomogeneous gas distribution, the optical depth of each shell across the nebula is multiplied by the chosen filling factor (see, for instance, Osterbrock 1989). The net effect corresponds to a decrease of the average density, increasing the size of the ionized region, particularly the low ionization one, as it will be discussed below.

In the following, the expression ‘model with a filling factor’ refers to a model with a filling factor less than unity. Because of the increase of the low ionization zone relative to the high ionization zone, this kind of model tends to overestimate the low-ionization lines. Therefore, when simulating the NLR we must keep in mind this systematic effect and how it can affect the conclusions reached

A multi-cloud model offers a more realistic view of the NLR, as long as a proper choice of single-clouds (densities and ionization parameters) is achieved. Assuming that the chemical composition is the same, each group of single-clouds contributing to the emission-line spectrum is characterized by the density and ionization parameter (ratio of ionizing photons to gas density). The simulated emission-line spectrum corresponds to the weighted average of the line intensities coming from each group. The quality of the fit is related to the number of observed emission-lines. A better fit to the observed spectrum is obtained when a larger number of observed lines is available, allowing a fine tuning of the physical conditions producing the lines, leading to a better choice of the groups of single clouds.

Recently, NLR photoionization models have been used to constrain the shape of the spectral energy distribution of active nuclei (Oliva et al. 1994, Alexander et al. 1999 and 2000, Prieto & Viegas 2000). Here, our goal is to discuss the dependence of the SED on the type of model used to describe the NLR. Before analyzing the emission-line spectrum of two Seyfert galaxies (NGC 4151 and NGC 1068), we discuss the filling factor problem related to photoionization models for extragalactic HII regions. For these objects the ionizing radiation spectrum is better known and the chemical composition obtained by empirical methods is reliable. Photoionization simulations of the NW knot of I Zw 18 provide a good example of the effect of a filling factor on the gas temperature and how models with a filling factor may induce to misleading results.

In §2, models for an extragalactic HII region are used to illustrate the effect of the filling factor, and a model for the NW knot of I Zw 18 is discussed. In §3, the observational data for NGC 1068 and NGC 4151 are presented, whereas photoionization models appear in §4. The concluding remarks are presented in §5.

## 2 Filling Factor Models

The effect of $`ϵ`$ $`<`$ 1 in 1-D models is to increase the size of the low-ionization zone, consequently increasing the ionized region and reducing the average temperature of the ionized gas. Models with $`ϵ`$ = 1 should provide a higher average gas temperature, increasing the calculated \[O III\] line ratio. In order to illustrate the effect of a filling factor less than unit, the results of HII region models with $`ϵ`$ $``$ 1, obtained by the 1-D code Aangaba (Gruenwald & Viegas 1992), are discussed.

We adopt the ionizing radiation of a zero age stellar cluster from Cid-Fernandes et al. (1992) shown in Figure 1 (top panel), with a number of ionizing photons QH = 1051 s-1, a uniform distribution of the gas with hydrogen density equal to 100 cm-3, and chemical abundances of I Zw 18 (Izotov & Thuan 1998). The results obtained for the outer radius of the HII region, as well as the average temperature TOIII, as a function of $`ϵ`$, are shown in Figure 1 (middle panel), while the ionic fractions of O0, O+ and O++ appear in Figure 1 (bottom panel). As expected, the outer radius, O0 and O+ increase with decreasing $`ϵ`$, while TOIII and O++ decrease, corresponding to the expansion of the low-ionization zone for smaller values of $`ϵ`$.

The shape of the ionizing radiation is another factor defining the relative size of the high-ionization to the low-ionization zone. To exemplify this point the zero age stellar cluster continuum spectrum was modified between 1 Ryd and 4 Ryd (Fig. 1, top panel) and the results are also shown in Figure 1 (middle and bottom panels). The filling factor effect is similar for the two types of ionizing spectra. However, the ionic fractions of O+ and O++ are larger while the O0 fraction is lower when the spectrum has less photons near the Lyman limit. Notice that, in this case, O++/O increases by a larger amount. Ions of other elements follow the same behavior. Therefore, the increase of the low-ionization zone due to a filling factor less than unity, leading to stronger low-ionization emission-lines, can be compensated by adopting an ionizing radiation with less photons near the Lyman limit.

### 2.1 The Star Forming Galaxy I Zw 18

Recently, 1-D photoionization models for the NW knot of I Zw 18, a metal poor extragalactic HII region, have been used by Stasinska & Schaerer (1999, hereafter SS99) to analyze its emission-line spectrum. They concluded that in addition to photoionization, another heating mechanism is necessary in order to explain all observed features, i. e., the size of the ionized region, the H$`\beta `$ luminosity and the emission-line ratios of the more significant lines.

In order to model the NW knot, SS99 adopted a homogeneous gas distribution with density equal to 100 cm-3, derived from the observed \[S II\] line ratio, as well as the chemical abundances obtained by Izotov & Thuan (1998).The best fit model should reproduce the significant emission-line ratios, the H$`\beta `$ luminosity, as well as the size of the ionized region. Adjusting the total initial mass of the stellar cluster in order to reproduce the luminosity, the authors show that a homogeneous model may explain some of the line ratios, but the calculated size of the ionized region is too small unless a filling factor of the order of 0.01 is adopted. Their best fit model corresponds to an homogeneous model with this filling factor, combined with condensations simulated by higher density models. All the observed features are explained except for the observed \[O III\] $`\lambda `$4363/5007 line ratio, the calculated ratio being too low. Because this ratio is a temperature indicator, the authors concluded that an additional heating source must be present in the NW knot and suggested the possible presence of shocks.

Following Viegas, Contini & Contini (1999), shocks must be contributing to the observed emission-lines of starburst galaxies depending on the evolutionary phase of the stellar cluster and could be identified as an additional heating mechanism. However, I Zw 18 is a metal poor extragalactic HII galaxy where the shock contribution to the physical conditions is negligible. Such effects could only be significant in evolved star forming regions. In order to analyze the apparent discrepancy between the need for an additional heating source and the characteristics of I Zw 18, the SS99 model assumptions have been reviewed (see Viegas 2002 for more details), using the code Aangaba.

Following SS99, $`ϵ`$ $`<`$ 1 is necessary in order to match the size of the ionized region. However, a model with $`ϵ`$ = 1 and a lower density could reproduce it. Notice that a density of about 30 cm-3 is also in agreement with the observed \[S II\] ratio. In fact, Viegas (2002) showed that a uniform density model for a spherical region with n = 30 cm-3 and filling factor equal to unit, photoionized by a 5.4 Myr stellar cluster, reproduces the size of the NW knot and the H$`\beta `$ luminosity, although some of the observed line ratios are not explained. The same kind of discrepancy between observational and theoretical line ratios was reported by SS99 concerning the homogenous model. These authors suggested the presence of condensations to improve the fitting. In fact, an inhomogeneous gas distribution is suggested by the HST image of I Zw 18. Two other models with a higher density (104 cm-3), mimicking possible condensations, are also presented by Viegas (2002), showing that the emission-line intensities, the H$`\beta `$ luminosity, as well as the size of the emission region, are well reproduced with no need of an additional heating mechanism.

No specific combination of dense and dilute gas is provided by SS99, but they show that the presence of condensations improves the results for the \[O I\] emission line, which was too low in the homogeneous model. The results obtained by Viegas (2002) with dense gas (104 cm-3) at different distances from the stellar cluster indicated that the best choice is a set of condensations located in the outer zone of the ionized region. Notice that a large fraction of H$`\beta `$ is emitted in the inner zone. Taking these facts into account, an example of a composite model for the NW knot of I Zw 18 is presented in Figure 2, where the results for a dilute uniform gas (30 cm-3, filled circle) and a denser gas (104 cm-3, filled square) are combined (circle) assuming that 60% of H$`\beta `$ comes from the low density gas and the remaining 40% from condensations. This model provides a very good fit of the line ratios, as well as of the size of the ionized region and of the average TOIII. The results presented in Fig. 2 make clear the potential of multi-cloud models to reproduce inhomogeneous regions without the problems associated with models with a filling factor less than unity. It is not our intention to claim that this is the best model for the NW knot of I Zw 18, but just to exemplify how an artificial problem created by models with filling factor can be solved by composite models. Notice, however, that a reliable composite model requires the fitting of large number of observed lines, not only the most significant ones.

In brief, the adoption of a filling factor led to a “heating problem” and, consequently, to the suggestion of the presence of shocks, whereas a composite model, built from single-cloud models with different densities and $`ϵ`$ = 1, explains all the observable features without requiring an additional heating mechanism, in agreement with the conclusions from Viegas et al. (1999). In this case, the misleading effect of a filling factor less than unity – artificially increasing the low-ionization zone of the cloud and decreasing the average temperature of the gas – could be easily identified because the ionizing radiation spectrum is well known. However, for the NLR of AGN the “filling factor effect” could influence the choice of the shape of the unknown UV ionizing radiation spectrum inferred from the models, as discussed below for two Seyfert galaxies, NGC 1068 and NGC 4151.

## 3 NGC 4151 and NGC 1068 Observational Data

NGC 4151 is a spiral galaxy, almost face-on, with z = 0.0033 and magnitude mV = 11.5. Its continuum spectrum is dominated by a non-stellar component (Kaspi et al. 1996) which is extremely variable, as are the broad emission lines. This galaxy is usually given as an example of a classic Seyfert galaxy, but it shows a variety of physical properties of different AGN classes. It could then be more complex than the average Seyfert galaxies. Originally classified as an intermediate Sy 1.5 galaxy (Osterbrock & Koski 1976), it presented the characteristics of a Sy 2 during its extremely low state in 1984 (Penson & Perez 1984). Currently NGC 4151 is classified as a Sy 1 galaxy.

NGC 1068 is one of the brightest and closest known Sy 2 galaxies. It is also the archetypal “hidden Seyfert 1” galaxy, following the unified model for active galactic nuclei (Antonucci & Miller 1985, Antonucci 1993). The narrow-line region is approximately co-spatial with a linear radio source with two lobes (Wilson & Ulvestad 1983). Most of the NLR and extended emission line region are likely powered by photoionization by a central source (Marconi et al. 1996). Moreover, high-resolution observations suggest kinematic disturbance and possibly shock excitation of regions close to the radio outflow (e.g. Axon et al. 1998).

Both objects have already been modeled by Alexander et al. (1999) and Alexander et al. (2000) (hereafter A99 and A00, respectively), using photoionization models with a filling factor. These authors used the observed emission-lines to constrain the ionizing continuum spectral shape, in particular to discuss the presence of a blue bump. On the other hand, two-component photoionization models, with a broken power law, have been used by Kraemer & Crenshaw (2000) to reproduce the emission-line spectra of NGC 1068. Multi-cloud models, which take into account the coupled effect of photoionization by a power-law ionizing radiation and shocks, are also discussed by Contini et al. (2002) for NGC 4151.

In order to discuss the filling factor effect on the inferred SED of the ionizing radiation, NLR of these two AGN are modeled in this work. We will adopt the list of observed emission-lines of NGC 4151 and NGC 1068 from A99 and A00, respectively. They provide a list of emission-lines carefully selected according to the following criteria: (a) only lines that are primarily formed by photoionization due to an UV – soft X-ray continuum from a central source are included; (b) all lines emitted from ions with ionization potential less than 13.6 eV are excluded, since these ions can be easily ionized by other processes; (c) only lines with reliable flux measurements are included, requiring that the line measurements are not too inconsistent with the large SWS and LWS (ISO Short- and Long-Wavelength Spectrometer, respectively) apertures; and finally, (d) they excluded the \[FeX\] 6734 $`\AA `$ and \[FeXI\] 7892 $`\AA `$ coronal lines, because the collision strength values are highly uncertain. Whenever more than one measurement for these lines exists, they quote the average flux and use rms scatter as an error estimate.

Notice that A99 compare the emission-line absolute fluxes of NGC 4151 to the model results, instead of the line intensity relative to a Balmer line. Their argument is that the narrow H$`\beta `$ flux measurement is very uncertain, due to the difficulty of decomposing the narrow and broad components. This is true indeed. Notice however that some of the references regarding the observational data present the emission lines relative to H$`\beta `$, without providing H$`\beta `$ absolute value (for example, Osterbrock & Koski 1976). In this case, A99 uses the absolute value from another data set. Furthermore, some of the references are extremely old, and clearly show serious problems with the measurements, mainly with the decomposition of the narrow and broad line components (for example, Oke & Sargent 1968 where H$`\alpha `$/H$`\beta `$ = 1.61). Because of these problems, we decided to exclude data coming from references where H$`\alpha `$/H$`\beta `$ is less than 2.8, since this clearly indicates a problem with the decomposition. In addition, as pointed out by A99, the narrow component of the UV lines (Kriss et al. 1992) show a FWHM which is 1.5 to 3 times larger than that of the NLR forbidden lines, which could also indicate a problem with the decomposition.

For NGC 1068, the absolute fluxes of the lines are provided by the authors in all references listed by A00. However, some measurements seemed to be problematic as well. Most of the observations were made in the 70’s, with photographic plates. Some of the measurements are very discordant, which can be seen by the rms scatter in the A00 list. Taking this into account, we decided to use only the very strong optical and infrared lines, excluding problematic data from the average.

Therefore, averaged line intensities in our sample are not exactly the same as in A99 or A00. We choose to compare the observed and calculated line ratios, and not the absolute fluxes, adopting the line intensities relative to H$`\beta `$. The line ratios corrected for reddening are presented in Table 1a and 1b, respectively, assuming E(B-V) = 0.03 for NGC 4151 and E(B-V) = 0.30 for NGC 1068, in agreement with those adopted by A99 and A00. In Tables 1a and 1b, the doublets are identified by the wavelength of one of the lines followed by a plus sign.

## 4 The Models

The photoionization code Aangaba (Gruenwald & Viegas 1992) is used to model the NLR of NGC 4151 and NGC 1068, assuming a plan-parallel symmetry. For each AGN, we first reproduce the results of the best fit given by A99 and A00, using the same input parameters, as listed in Table 2, where U is the ionization parameter, n the gas density, $`ϵ`$ the filling factor, and Ri the distance of the gas to the central source. Solar chemical abundances are adopted for modeling NGC 4151. On the other hand, for NGC 1068 a low oxygen abundance is assumed, as suggested by Netzer (1997) and adopted by A00 for their best fit model. From their models with a filling factor, A99 and A00 proposed that the continuum of these objects should be represented by a composition of power laws, as shown in Figures 3a and 3b, for NGC 4151 and NGC 1068, respectively. The continuum spectral shapes above 103 eV (X-rays) and below 10 eV (optical-UV) are obtained from observations, while between 10 and 104 eV several different shapes have been tested.

After reproducing the results of the best fit models proposed by A99 and A00 for NGC 4151 and NGC 1068, our next step is to explain the emission-line spectra using a multi-cloud model (based on single-cloud models with $`ϵ`$ = 1), in order to test the spectral shapes suggested by A99 and A00, who excluded the presence of a blue bump.

### 4.1 NGC 4151 models

As claimed by A99, the spectrum showed by a solid line in Figure 3a (spectrum A) provides the best fit model for the observed narrow emission-lines of NGC 4151. Those authors did not find indications of a big blue bump in the ionizing continuum, which should appear beyond the Lyman limit. On the contrary, their best fit model suggests a trough in a power law.

As pointed out in §2, the trough in the ionizing spectrum could be an artifact which compensates for the systematic effect of a filling factor less than unity in the 1-D models. Thus, we look for multi-cloud models with another ionizing spectrum shape in order to verify if they are also able to reproduce the observed emission-lines.

For sake of comparison to a multi-cloud model, we reproduce the A99 best fit model using the code Aangaba with the input parameters listed in Table 2. Notice that our data set is slightly different from that used by A99, as explained above (§3), in particular the H$`\alpha `$/H$`\beta `$ ratio, which is equal to 2 in their data set. In spite of the differences, the fitting to the data is as good as that presented by A99, as seen in Figure 4 where observed and calculated line ratios are confronted (to be compared with Fig. 8 of A99).

Models using other spectral shapes were also obtained. Among them, the model with spectrum B (Fig. 3a), showing a small blue bump, also provides a good fit to the observed lines, except for the high ionization lines, as expected from models with such type of continuum and assuming a filling factor less than unit. However, as in the case of the extragalactic HII region I Zw 18, we cannot exclude this kind of spectral shape, since the underestimation of the high ionization lines can also be due to $`ϵ`$ less than unit. In order to test this issue, multi-cloud models are also obtained. For this a grid of single-cloud models is built varying the density and the ionization parameter U.

The simplest composite model includes a low- and a high-density cloud, as adopted by Kraemer & Crenshaw (2000) in order to reproduce the the optical-UV emission-line spectrum of NGC 1068.

A grid of models for single-clouds characterized by the density (n) and the ionization parameter (U) was obtained. The models were generated with these parameters in the range 102cm-3 $``$ n $``$ 105cm-3 and 0.005 $``$ U $``$ 4.0. To restrict the physical characteristics of the models we used the strong lines \[OIII\]$`\lambda `$5007+ and \[MgVIII\]$`\lambda `$ 3.03 $`\mathrm{\mu m}`$. For a given value of U and varying the density, we look for the best two-density models which can explain these emission-lines. Once the densities are derived, we used all the observed emission-line intensities (Table 1a) to obtain the values of U for the low and high density clouds. The emission-line spectrum is reproduced by a four-cloud model. Notice that line ratios are not method dependent, i.e., the same cloud composition is found if we first define the values of U and then the cloud densities. The line ratios are obtained by a weighted average from the contributions from each of the four types of cloud. The ionization parameter, the gas density, and the H$`\beta `$ flux, as well as the fraction of H$`\beta `$ (gi) coming from each type of cloud, are listed in Table 3. This fraction is related to the area of the corresponding group of clouds photoionized by the central source. The largest contribution to the emission-line spectrum comes from the low-ionization, high density clouds.

A comparison between the observed and calculated emission-line ratios appears in Figure 5, where the results of A99 best fit are also plotted. It is easily seen that both the A99 model, using an ionizing SED with a trough and $`ϵ`$ $`<`$ 1, and our four-cloud model, using an ionizing radiation SED with a bump (spectrum B) and $`ϵ`$ = 1, provide a similar fit to the observations. This can be verified by the $`\chi ^2`$ value, which are 2.27 and 2.28 respectively for the A99 best fit model and our multi-cloud model.

The multi-cloud model presented above is an example of the method used to fit the emission-lines and shows that the ionizing spectrum may show a bump instead of a trough as claimed by A99. In order to have a better determination of the ionizing continuum, additional constraints, coming from the observed continuum spectrum from radio to X-rays, could be added to the model, as well as the contribution of an additional energy source, as discussed by Contini et al. (2002).

### 4.2 NGC 1068 models

In the case of NGC 1068, A00 do not find a single-component model that satisfactorily reproduces the observed emission lines. For this object, their best fit is a two-component model; one of the components is calculated with a radially varying filling factor. It is out of the scope of this paper to reproduce such a complex model for this object, since our goal is just to test the two approaches used to reproduce the emission-line spectrum of an inhomogeneous region with 1-D photoionization simulations. Therefore, we used input parameters of the best single-component model of A00 to discuss these two approaches, using the two shapes of the ionizing continuum plotted in Fig. 3b.

Spectrum A (Fig. 3b), with a trough, was assumed in the A00 single-cloud model which gives the best fit to the observations with the input parameters listed in Table 2. Spectrum B is an example of a bumpy continuum.

Following the procedure used for NGC 4151, the emission-line spectrum of NGC 1068 is first modeled using an uniform density distribution and a filling factor $`ϵ`$ $`<`$ 1. The input parameters are listed in Table 2. The general conclusion is similar to the one obtained for NGC 4151, i.e., regardless of the spectral shape of the ionizing radiation, single-cloud models with $`ϵ`$ $`<`$ 1 overestimate the low ionization lines with respect to the high ionization lines, as already pointed out by A00. The comparison between observed and calculated emission-line intensities are shown in Figure 6 for both types of ionizing spectra.

We now look for the best fit multi-cloud model using the bumpy spectrum (spectrum B, Fig.3 ) in order to compare the results to those of A00 best fit model with ionizing spectrum A and $`ϵ`$ $`<`$ 1. Following the same procedure adopted for NGC 4151 multi-cloud model, four types of clouds are defined. The characteristics of the single clouds, as well as the fraction of H$`\beta `$ coming from each type of cloud, are presented in Table 4. The larger contribution to the total H$`\beta `$ comes from the low-ionization low-density cloud.

The comparison between observed and calculated line ratios (Figure 7) shows that, except for the infrared lines \[NeIII\] 15.56 and \[NeIII\] 36.01, our four-cloud fit, photoionized by a bumpy spectrum, is better than that provided by the best A00 single-cloud model. The corresponding $`\chi ^2`$ value for these models are 2.88 and 7.89, respectively.

Notice that A00 find that this simple filling factor model was not good enough to reproduce the observed spectrum. They assumed a more complex two-component model, obtained from models with a filling factor. The low-density component has a constant filling factor, while the denser one a filling factor dependent on the distance. It is clear that our simpler multi-cloud model provides a very good fit to the observed line spectrum without the problems associated to and complexity required by their models with a filling factor.

## 5 Concluding Remarks

The emission-line spectra of the two well-known Seyfert galaxies NGC 4151 and NGC 1068 are used to test 1-D photoionization models for the NLR and to illustrate the misunderstanding introduced by the use of models with a filling factor less than unity to simulate inhomogeneous regions.

The ionizing radiation of AGN is not well known, but its shape, combined with the filling factor as it is introduced in 1- D photoionization codes, has a direct effect on the calculated emission lines. In the case of the extragalactic HII regions, for which the ionizing continuum is better known, we show that a multi-cloud model can explain all the observed features (size of the emitting region, H$`\beta `$ luminosity and the significant emission-line ratios) without requiring an additional heating mechanism as inferred by SS99 who used models with $`ϵ`$ $`<`$ 1.

Regarding the two Seyfert galaxies analyzed in this paper, we show that multi-cloud models are more appropriate to simulate the NLR. They provide a good fit to the observed emission-line spectra and do not exclude a bump in the SED of the ionizing radiation, as suggested by A99 and A00.

Two-component models have been used by Kraemer & Crenshaw (2000) to reproduce the optical-UV emission-line spectrum from several resolved zones of the NLR of NGC 1068. The hydrogen density of the dense components are in the range 2.0 $`\times `$ 104 to 2.5 $`\times `$ 105 cm-3 and the ionization parameter 10-3.8 $``$ U $``$ 10-2.5 , while for the tenuous component the density varies from 8.0 $`\times `$ 102 to 2.3 $`\times `$ 104 cm-3 and 10-2 $``$ U $``$ 10-1.6. The values of the density and ionization parameters used in our multi-cloud model are in the range of those obtained for the two-cloud models proposed by Kraemer & Crenshaw for several resolved zones of the NLR of NGC 1068. Notice that our model simulates the whole NLR and also explains the infrared lines not included in Kraemer & Crhenshaw discussion.

A more comprehensive multi-cloud model of the emission line spectrum of NGC 4151 is presented by Contini et al. (2002). The code used in the simulations includes the effect of shocks in the physical conditions of the emitting gas. The observed continuum data from radio to X-rays are used to constrain the model and to find the contribution from the different clouds to the observed continuum and emission-line spectra. Noticed that the presence of shocks is also suggested by Kraemer & Crenshaw (2000).

In brief, in the last years the improvement of observational techniques, allowing observations with high spatial and spectral resolution, and providing a large amount of data in a wide range of wavelengths, has not been followed by a similar effort to improve the interpretation of the emitting regions with more realistic models. The standard 1-D photoionization models are becoming inappropriate to simulate emission-line regions like HII regions, planetary nebulae and active galactic nuclei. The new observations require more complex models regarding both additional mechanisms powering the emission (as shocks) and a better simulation of the geometry of the region with a more realistic gas distribution. It is clear that 3-D photoionization simulations would be more appropriate to analyze the NLR of AGN, as shown in the case of planetary nebulae (Gruenwald, Viegas & Broguière 1997, Monteiro et al. 2000). Such 3-D models will certainly contribute to a better discrimination of the possible nature of the ionizing radiation.

We are thankful to T. Alexander and H. Netzer for clarifying some of their results and data. We thank also Claus Leitherer for the careful reading of this paper. This paper is partially supported by FAPESP (00/06695-0, 00/122117-4) and CNPq (304077/77-1 and 306122/88-0), PRONEX/CNPq.


text:
# Superconductivity by Kinetic Energy Saving ?

## I Model independent properties of the superconducting state

### I.1 Internal energy

When we cool down a superconductor below the critical temperature, the material enters a qualitatively different state of matter, manifested by quantum coherence over macroscopic distances. Because the critical temperature represents a special point in the evolution of the internal energy versus temperature, the internal energy departs from the temperature dependence seen in the normal state when superconductivity occurs. Because for $`T<T_c`$ the superconducting state is the stable equilibrium state, the internal energy in equilibrium at $`T=0`$ is an absolute minimum. Hence cooling down from above the phase transition one would expect a drop in the internal energy when at the critical temperature. This drop of internal energy stabilizes the superconducting phase among all alternative states of matter. An experimental example of this well known behavior is displayed in Fig. 1, where the internal energy was calculated from the electronic specific heat dataloram01  according to the relation $`E_{int}(T)=_0^Tc(T^{})𝑑T^{}`$.

The broadened appearance of the phase transition suggests that the superconducting correlations disappear rather gradually when the temperature is increased above the phase transition, an effect which remains noticeable in the specific heat graph up to 125 K.
Understanding the mechanism of superconductivity means to understand what stabilizes the internal energy of the superconducting state. In BCS theory superconductivity arises as a result of a net attractive interaction between the quasi-particles of the normal state. Note, that implicitly this approach is firmly rooted in the paradigm of a Fermi-liquid type normal state. On the other hand, the school based on Anderson’s original workanderson87  asserts, that a strong on-site repulsive interaction can also give rise to high $`T_c`$ superconductivity. The latter models typically require that the material is not a Fermi liquid when superconductivity is muted, either by raising the temperature or by other means.

### I.2 Pair-correlations

Without loss of generality, i.e. independent of the details of the mechanism which leads to superconductivity, it is possible to provide a microscopic definition of the superconducting state. For this purpose let us consider the correlation function defined as

$$G(r,R_1,R_2)=\psi _{}^{}(R_1+r)\psi _{}^{}(R_1)\psi _{}(R_2)\psi _{}(R_2+r)$$
(1)
where $`\psi _\sigma ^{}(R_j)`$ is a single electron creation operator. This function defines two types of correlation: (i) the electron pair-correlation as a function of the relative coordinate $`\stackrel{}{r}`$, and (ii)v the correlation between two pairs located at center of mass coordinates $`\stackrel{}{R}_1`$ and $`\stackrel{}{R}_2`$. In the normal state there is no correlation of the phase of $`G(r,R_1,R_2)`$ over long distances $`|R_1R_2|`$ due to the finite mean free path of the electrons. As a result the integral over the center of mass coordinates of the correlation function

$$g(r)=\frac{1}{V^2}d^3R_1d^3R_1G(r,R_1,R_2)$$
(2)
averages to zero in the normal state. In contrast the superconducting state is characterized by long range phase coherence of the center of mass coordinates, implying (among other things) that the correlation function averaged over all center of mass coordinates, $`g(r)`$, is a finite number.

## II Internal energy and its decomposition using the BCS model

### II.1 BCS: Correlation function, and correlation energy for d-wave pairing

In the weak coupling scenario of BCS theory the electrons have an effective attractive interaction, as a result of which they tend to form pairs. For the purpose of the present discussion we will assume that the interaction is of the form

$$H^i=V(rr^{})\widehat{n}(r)\widehat{n}(r^{})$$
(3)
where $`\widehat{n}(r)`$ is the electron density operator. The interaction energy in the superconducting state becomes lower than in the normal state, due to the fact that the effective attractive interaction favors a state with enhanced pair-correlations. The value of the interaction energy of the superconducting state, relative to the normal state is

$$H^i_sH^i_n=d^3rg(r)V(r)=\underset{k}{}g_kV_k$$
(4)
where $`g_k`$ and $`V_k`$ are the Fourier transforms of $`g(r)`$ and $`V(r)`$ respectively. Using the Bogoliubov transformation the correlation function can be expressedparticle-hole  in terms of the gap-function $`\mathrm{\Delta }_k`$ and the single particle energies $`E_k=\{(ϵ_k\mu )^2+\mathrm{\Delta }_k^2\}^{1/2}`$.

$$g_k=\underset{q}{}\frac{\mathrm{\Delta }_{q+k}\mathrm{\Delta }_q^{}}{4E_{q+k}E_q}$$
(5)
This corresponds to the conversion of a pair $`(q,q)`$ to a pair with quantum numbers $`(q+k,qk)`$. In the expression for the correlation energy, Eq. 4, the transferred momentum $`k`$ is carried by interaction kernel $`V_k`$.
Starting from a model expression for the single electron energy-momentum dispersion $`ϵ_k`$, and the gap-function $`\mathrm{\Delta }_k`$, it is a straightforward numerical exercise to calculate the summations in Eq. 5. Adopting the nearest neighbor tight-binding model with a d-wave gap, and adopting the ratio $`\mathrm{\Delta }(\pi ,0)/W=0.2`$, where $`W`$ is the bandwidth, the correlation function $`g_k`$ can be easily calculated, and the result is shown in Fig. 2.

We see from this graph that a negative value of $`H^i_sH^i_n`$ requires either (i) $`V_k<0`$ for $`k`$ in the neighborhood of the origin, or (ii) $`V_k>0`$ for $`k`$ in the vicinity of $`(\pi ,\pi )`$. The corresponding representation in real space, $`g(r)`$, shown in Fig. 2, illustrates that the dominant correlation of the d-wave superconducting state is of pairs where the two electrons occupy a nearest neighboring site, while the on-site amplitude is zero.
Combining the information of Fig. 2 with Eq. 2, it is clear that the strongest saving of correlation energy is expected if the electrons interact with an interaction of the form $`V(r_1,r_2)=V_0_a\delta (\stackrel{}{r}_1\stackrel{}{r}_2+\stackrel{}{a})`$ where the vector $`a`$ runs over nearest neighbor sites, and $`V_0`$ is a negative (meaning attractive) interaction between electrons on nearest neighbor sites.

### II.2 BCS: The gap-equation, specific heat and internal energy

To illustrate the predictions of BCS theory for the temperature dependence of the correlation energy and the kinetic energy, we start by solving the gap equation

$$\mathrm{\Delta }_k=\underset{q}{}\frac{V_{kq}\mathrm{\Delta }_q}{2E_q}\mathrm{tanh}\left(\frac{E_k}{2k_BT}\right)$$
(6)
which must solved together with the constraint that the average number of particles is temperature independent. This requires that the chemical potential must be adjusted to keep the thermal average of $`_k\widehat{n}_k=N`$ at a constant value. In the numerical examples of this paper we have adopted $`N=0.85`$ corresponding to a hole doping of 0.15. This solution of the gap equation is shown in Fig. 3.

Notice, that the temperature dependence of the chemical potential follows closely the experimental observations reported for Y123 in Ref. rietveld92
To check that the band-parameters used here are reasonable, we display in Fig. 4 the corresponding prediction for the specific heat and the internal energy, using the same parameters as for Fig. 3. If we compare this to Fig. 1, we see that the band-parameters adopted here quantitatively reproduce the observed specific heat. Hence the present set of band-parameters ($`t`$,$`t^{}`$) and the coupling parameter $`V_0`$ represent the best phenomenological choice for quantitative testing of the BCS-model.

### II.3 BCS: Temperature dependence of the correlation energy and the kinetic energy

The BCS prediction for the temperature dependence of the average interaction energy follows from Eq. 4. We then use the BCS variational wave-function for the statistical average of Eq. 5, resulting in

$$H^i_sH^i_n=\underset{k}{}\frac{|\mathrm{\Delta }_k|^2}{2E_k}\mathrm{tanh}\left(\frac{E_k}{2k_BT}\right)$$
(7)
Simultaneously there is an increase of the ’kinetic’ energy

$$H^{kin}=\underset{k}{}ϵ_k\left\{1\frac{ϵ_k\mu }{E_k}\mathrm{tanh}\left(\frac{E_k}{2k_BT}\right)\right\}$$
(8)
In Fig. 5 this is displayed, using the same parameters as for Fig. 3.

## III Relationship between intra-band spectral weight and kinetic energy

A measure of the kinetic energy is provided by the following relationnorman02 ; kubo57 ; drew97

$$_\mathrm{\Omega }^\mathrm{\Omega }𝑑\omega \text{Re}\sigma (\omega )𝑑\omega =\pi \frac{e^2}{\mathrm{}^2V}\underset{k,\sigma }{}\widehat{n}_{k,\sigma }\frac{^2ϵ(\stackrel{}{k})}{k^2}$$
(9)
where the high frequency limit indicates that the integral should include only the intra-valence band transitions, and the condensate peak at $`\omega =0`$ if the material is a superconductor. The integral over negative and positive frequencies (note that $`\sigma (\omega )=\sigma ^{}(\omega )`$) avoids ambiguity about the way the spectral weight in the condensate peak should be counted. If the band structure is described by a nearest neighbor tight-binding model, Eq. 9 leads to the simple relation

$$\begin{array}{c}\rho _L\frac{\mathrm{}^2}{a^2\pi e^2}_\mathrm{\Omega }^\mathrm{\Omega }\text{Re}\sigma (\omega )𝑑\omega =H_{kin}\hfill \end{array}$$
(10)
Hence in the nearest neighbor tight-binding limit the partial f-sum provides the kinetic energy contribution, which depends both on the number of particles and the hopping parameter $`t`$maldague77 ; carmelo88 ; basov99 ; chakravarty99 .

However, if the band-structure has both nearest neighbor hopping and next nearest neighbor hopping, Eq. 10 is not an exact relation, and instead Eq. 9 should be compared directly to the experiments. In Fig. 6 we compare the spectral weight, calculated directly using Eq. 9 to the result of Eq. 10, using the same parameters as for Fig. 3.

Note that the kinetic energy has to be divided by a factor two, as we are interested in the projection along one of the two axes in the ab-plane, which can be compared directly to the experimental value of $`\rho _L`$. From Fig. 6 we can conclude, that the effect of including $`t^{}`$ in the calculation is rather small, and it is still OK to identify $`\rho _L`$ with the kinetic energy apart from a minus sign.

## IV Experimental determination of the intra-band spectral weight

In two recent experimental papers measurements of $`\rho _L`$ in Bi2212 have been reported molegraaf02 ; santander02 . The values of the kinetic energy change in the superconducting state were in quantitative agreement with each other, and both papers arrived at the same conclusion: Contrary to the BCS prediction, the kinetic energy of the superconducting state is lower than in the normal state (taking into account a correction for the temperature trends of the normal state). In Fig. 7 the data of Ref. molegraaf02  have been reproduced. Comparing this with the BCS prediction clearly demonstrates the large qualitative discrepancy between theory and experiment. Clearly the type of mechanism assumed in BCS theory is not at work here.

## V Implications of the experimental data

The trend seen in the experimental data has been predicted by Hirsch in 1992hirsch92a ; hirsch92b ; hirsch98 . The model assumption made by Hirsch was, that the hopping probability of a single hole between two sites becomes larger if one of the two sites is already occupied by a hole. Although this model provides good qualitative agreement with the optical experiments, it has one serious deficiency: It also predicts s-wave symmetry for the order parameter, in sharp contrast to a large body of experimental data which show that the superconducting gap in the cuprates has d-wave symmetry.
In a recent set of calculations based on the Hubbard model, Jarrell et aljarrell03  obtained a similar effect as seen in our experiments, both for underdoped and optimally doped samples. Crudely speaking the mechanism is believed due to the frustrated motion of single carriers in a background with short-range (RVB-type) spin-correlations, which is released once pairs are formed.

## VI conclusions

We have made a quantitative comparison between the BCS-prediction for d-wave pairing in a band with nearest neighbor and next-nearest neighbor hoppping and various experiments, in particular specific heat and measurements of the optical ab-plane sumrule. We have shown that the BCS-prediction produces the wrong sign for the kink at $`T_c`$ of the ab-plane intraband spectral weight, while the model calculation is in good agreement with the experimental specific heat data.

###### Acknowledgements.

DvdM gratefully acknowledges M. Norman, N. Bontemps, J. Hirsch, and M. Jarrell for stimulating discussions, and J. W. Loram for making his data files of the specific heat of Bi2212 available.


text:
# Ab initio study of Curie temperatures of diluted III-V magnetic semiconductors

## I INTRODUCTION

Because of their ferromagnetic properties, diluted magnetic III-V semiconductors (DMS) represent a new class of materials with potential technological applications in spintronics . The physics of DMS is in particular interesting because of the yet unclear origin of the occurring ferromagnetism and the even less understood role of disorder. It is now generally accepted that the ferromagnetic coupling in III-V DMS containing Mn impurities is mediated via valence band holes in the host semiconductor. Mn atoms that substitute three-valent cations in the host act as acceptors which in turn create holes in the valence band. For finite concentration of Mn atoms the Fermi energy is then pinned within the valence band. The ferromagnetic coupling of Mn atoms is often qualitatively explained by the fact that because of the small size of the corresponding hole Fermi surface the period of the Ruderman-Kittel-Kasua-Yosida (RKKY) oscillations exceeds the typical distance between Mn atoms.

Theoretical approaches to the DMS can be divided into two groups, namely model studies using the effective-mass approximation, and studies based on the detailed knowledge of the electronic structure. The model studies mostly employ the kinetic-exchange model (KEM) in the connection with a continuum approximation for the distribution of Mn atoms and other defects , a model that yields a disorder-free problem, although it was recently refined in order to include disorder via a supercell method by means of Monte Carlo simulations . A simpler model assuming that holes hop only between Mn-sites, where they interact with the Mn moments via phenomenological exchange interactions, was also developed .

Electronic structure studies based on the spin-density functional theory (DFT) in the framework of the local spin-density approximation (LSDA) represent a further step towards a more detailed, parameter-free description of the properties of DMS. One possibility of describing such system is a supercell approach , in which large cells are used to simulate semiconductor crystals with magnetic atoms and other impurities. The supercell approach has limited applicability because it leads to excessively large supercells for low concentrations of impurities characteristic for DMS. Alternatively, one can employ Green function methods in which the configurational averaging is performed within the coherent potential approximation (CPA)  and the resulting averaged Green function can be used to determine any physical quantity of interest. The Green function approach can be implemented, e.g., in the framework of the Korringa-Kohn-Rostoker (KKR) method  or the tight-binding linear muffin-tin orbital method (TB-LMTO)  as in the present case. DMS with Curie temperatures ($`T_c`$) in the range of the room temperature are needed for practical applications; the study of realistic models for evaluation of $`T_c`$ is thus of great importance. The evaluation of the Curie temperature on ab initio level is still a challenging task, in particular for complex alloy systems like the DMS. Conventional approaches known from model studies  cannot be directly employed for evaluation of $`T_c`$. Therefore, a different approach has to adopted for the LSDA Hamiltonian. It is well-known  that transversal fluctuations of the local magnetization have to be taken into account to describe temperature-dependent properties of magnetic itinerant systems. The direct evaluation of the corresponding enhanced frequency-dependent magnetic susceptibility  is very complicated, in particular for such complex random systems like DMS. A tractable approach is based on the adiabatic treatment of atomic moments in which noncollinear configurations of magnetic moments are taken into account. The adiabatic approximation is well justified for magnetic atoms with large exchange splittings, like, e.g., Mn-atoms. This method can be implemented either in the real space  or in the reciprocal space (the frozen-magnon approach) . The real-space approach can be naturally applied to random system in the framework of the CPA while in the reciprocal-space approach the supercell method has to be adopted .

In the real-space approach employed here the total energy is mapped onto a classical Heisenberg model whose statistical treatment represents a next step. The mapping is further simplified by using the magnetic force-theorem  which allows to use the band energy of the calculated ground state as a good estimate also for total-energy differences of the perturbed ground state. The separated treatment of statistical degrees of freedom is another advantage of this approach as it can be done on different levels of sophistication. It is one of the goals of the present paper to describe evaluation of the Curie temperature of DMS in the framework of the above mentioned two-step model.

DMS containing Mn atoms, such as (Ga,Mn)As mixed crystals, are highly compensated systems. The experimentally observed number of holes in the valence band is much lower than the concentration of Mn-impurities and indicates the presence of other lattice defects acting as donors. The As-antisites, i.e., As atoms on the Ga-sublattice, add two electrons to the valence band for each defect compensating thus two holes. Interstitial Mn atoms acting as double donors have the same effect  as the As-antisites. Recent theoretical studies indicate that magnetic structures with partial disorder of local magnetic moments can lower the total energy of system as compared to the ferromagnetic state. The stability of the ferromagnetic state was studied in  within the KEM model concluding that non-collinear ferromagnetism is common in DMS. A new magnetic state stabilized by the As-antisites and characterized by a partial disorder in the orientations of the local Mn-magnetic moments was found recently also in first-principles studies of Ga0.96Mn0.04As alloys . All this indicates the great importance of various defects for the electronic structure, magnetic properties, and the Curie temperature of the DMS. A detailed understanding of such influence calls for a parameter-free approach. A detailed study of the influence of As-antisites is another goal of the present paper.

## II ELECTRONIC STRUCTURE

In this section we present results of electronic structure calculations, in particular the magnetic phase diagram, magnetic moments, and densities of states for the whole range of physically interesting concentrations of Mn-impurities and As-antisites.

### A Model and computational details

We have determined the electronic structure of DMS in the framework of the ab initio all-electron tight-binding linear muffin-tin orbital (TB-LMTO) method within the atomic-sphere approximation . The valence basis consists of $`s`$-, $`p`$-, and $`d`$-orbitals, we include scalar-relativistic corrections but neglect the spin-orbit effects. The substitutional disorder due to Mn-atoms and As-antisites is treated in the CPA. Although the CPA neglects local environment effects and possible lattice relaxations it correctly reproduces concentration dependent trends and can also treat systems with small but finite concentrations of defects typical for DMS. We introduce empty spheres  into the interstitial positions of the zinc-blende GaAs semiconductor for a better space filling. The lattice constant of the pure GaAs ($`a=5.653\mathrm{\AA }`$) was used in all calculations but we have verified that we can neglect a weak dependence of the sample volume on defect concentrations. We used equal Wigner-Seitz radii for all atoms and empty spheres. Charge selfconsistency is achieved within the framework of the local spin density approximation using the Vosko-Wilk-Nusair parametrization for the exchange-correlation potential . The details of the method can be found in Ref. .

In general, in addition to the chemical disorder, the DMS are characterized also by some degree of magnetic disorder. We treat this disorder in the framework of the disordered local moment (DLM) method  which is the simplest way of including disorder in spin orientations. The DLM can be included in the framework of the CPA: the Mn atoms have collinear but random positive (Mn+) and negative (Mn-) orientations. The corresponding concentrations, $`x^+`$ and $`x^{}`$, fulfill the condition $`x=x^++x^{}`$, where $`x`$ is the total Mn-concentration. The degree of magnetic disorder can be characterized by the order parameter $`r=(x^+x^{})/x`$ which is directly related to the macroscopic magnetization. Clearly, $`x^\pm =(1\pm r)x/2`$. The GaAs mixed crystal with prescribed concentrations of $`x`$ Mn- and $`y`$ As-atoms on the Ga-sublattice can thus be treated formally as a multicomponent (Ga1-x-yMn$`{}_{(1+r)x/2}{}^{}{}_{}{}^{+}`$Mn$`{}_{(1r)x/2}{}^{}{}_{}{}^{}`$Asy)As alloy with a varying parameter $`r`$, $`0r1`$. In the saturated ferromagnetic (FM) state, $`r=1`$, all magnetic moments are aligned in the direction of a global magnetization. The paramagnetic (PM) state, $`r=0`$, is characterized by a complete disorder of spin directions with vanishing total magnetization. The magnetic state with $`0<r<1`$ , the partial ferromagnetic (pFM) state, is characterized by some degree of disorder in the spin-orientations . It should be noted that the DLM method is used here only to determine the corresponding ground magnetic state in the presence of impurities. The Curie temperature is, however, determined from the Heisenberg Hamiltonian corresponding to a given magnetic ground state where fluctuations of local magnetization are taken into account.

### B Magnetic phase diagram

We have determined the order parameter $`r`$ that corresponds to the ground state at $`T=0`$ for the whole range of physically relevant concentrations $`(x,y)`$ in steps of $`\mathrm{\Delta }x=0.01`$ and $`\mathrm{\Delta }y=0.0025`$. In order to find the order parameter $`r`$ we varied for each composition $`(x,y)`$ the ratio $`r`$ in steps $`\mathrm{\Delta }r=0.05`$ and evaluated the corresponding total energy. The magnetic phase diagram determined in this way is presented in Fig. 1a. It should be noted that in addition to the FM and PM states there exist the pFM state with a certain degree of disorder in the spin-orientations. The pFM state, which separates the FM phase from the PM phase, appears in the concentration region where the number of holes is reduced by As-antisites ($`y>0`$). Without As-antisites the ground state is the FM state. For $`y>x/2`$, i.e., in the $`n`$-type materials, we did not find any indication for an electron-induced ferromagnetism even at the highest studied concentrations of donors. Further details can be read off from Fig. 1b, in which for a fixed Mn concentration of $`x=0.05`$ the total energy with respect to the ground state is displayed as a function of the order parameter $`r`$. For less than 1% of As-antisites the ground state is a saturated ferromagnet, for more than 2.35% of As-antisites a paramagnet, with a partial ferromagnet inbetween. Empty circles for $`y`$=0.015 in Fig. 1b denote results obtained for the experimental lattice constant of (Ga0.95Mn0.05)As alloy ($`a=5.672\mathrm{\AA }`$) and demonstrates that the ground state is not changed by a variation of the lattice constant.

### C Magnetic moments

The magnetic state of disordered magnetic semiconductors can be characterized by the local magnetic moment at the sites occupied by Mn atoms and by the total magnetic magnetization which involves both localized moments at magnetic impurities and the induced moments at non-magnetic atoms. The calculations show that the total magnetic moment per unit cell is proportional to the concentration of Mn atoms. This is the reason why we consider the total magnetic moment normalized to one Mn atom.

The calculated total and the local Mn-magnetic moments for (Ga0.95-yMn0.05Asy)As are presented in Fig. 2 as a function of the concentration $`y`$ of As-antisites. The size of the local moments on Mn atoms is almost independent of the chemical composition and depends only weakly on the direction of the moments. This justifies, a posteriori, the validity of the DLM model in the present case .

The total magnetic moment per Mn-atom equals to 4 $`\mu _B`$ in the system without As-antisites. The integer value of the moment is a consequence of the half-metallic character of the system. It is in agreement with the expectation that the localized moment of 5 $`\mu _B`$ combines with the opposite moment of one hole in the majority-spin valence band. The difference between total and local moments shows the degree of the magnetic polarization of the host material. Although the induced moments on non-magnetic atoms are generally very small, they add to a non-negligible contribution to the total magnetic moment.

In systems with compensating impurities, the total magnetic moment per Mn-atom increases linearly with $`y`$ and extrapolates to the value of 5 $`\mu _B`$ in the completely compensated ferromagnetic state. In reality, however, this situation is not reached because the transition from the FM to the pFM state takes place at $`y0.01`$. The transition is marked by a drop of the magnetization in the total magnetic moment for $`y>0.01`$. The reason for the decrease of the magnetization in the pFM state is obvious: the averaged local moment on Mn atoms, $`m^{\mathrm{Mn}}=x^+m^{\mathrm{Mn}^+}+x^{}m^{\mathrm{Mn}^{}}`$, is strongly reduced in the presence of Mn-sites with reversed moments. Clearly, the induced part of the magnetization also decreases proportionally to $`m^{\mathrm{Mn}}`$.

### D Densities of states

The total density of states (DOS) per unit cell and the local Mn-density of states for (Ga0.95Mn0.05)As in the FM state are presented in Fig. 3. We find a good agreement with existing results (see, e.g., Ref. ). In addition, we present the Mn-local DOS on the Ga-sublattice resolved with respect to orbital symmetries. The different behavior of $`t_2`$\- and $`e`$-states which are split due to the tetragonal environment is clearly seen. The $`t_2`$-states just above the top of the GaAs host valence band were experimentally observed . The influence of As-antisites is illustrated in Fig. 4. We have chosen $`y=0.015`$, i.e., the ground state is the pFM state with a partial disorder in local moments (see Fig. 1a). The local Mn+-DOS is quite similar to that in Fig. 3 because $`x^+0.044`$. On the contrary, the local Mn--DOS exhibits a sharp peak which is characteristic for a single-impurity limit ($`x^{}0.006`$). Such behavior also qualitatively explains the stability of the FM phase for systems without As-antisites: the Fermi energy in the pFM state without antisites, i.e., with larger number of holes in the valence band, would be shifted downwards and would be situated within a sharp peak of the minority Mn- states which is energetically unfavorable.

## III  CURIE TEMPERATURE

In this section we present two different ways of evaluating the mean-field Curie temperature using the selfconsistent DFT calculations. They are both implicitly based on the classical Heisenberg model, however, without the need for an explicit evaluation of the corresponding pair exchange interactions.

### A  Simple estimate of the Curie temperature

A simple estimate of $`T_c`$ can be obtained from DFT calculations by relating $`T_c`$ to the total energy difference per unit cell between the paramagnetic and the ferromagnetic state as obtained from selfconsistent calculations, i.e., $`\mathrm{\Delta }=E_{\mathrm{PM}}E_{\mathrm{FM}}`$. The DFT describes only the ground state properties, therefore an additional procedure is needed to estimate $`T_c`$. The classical Heisenberg model $`H=_{ij}J_{ij}𝐞_i𝐞_j`$, where $`𝐞_i`$ and $`𝐞_j`$ are unit vectors characterizing the orientations of the local magnetic moments at sites $`i`$ and $`j`$, and the $`J_{ij}`$ denote the effective exchange interactions between a pair of magnetic atoms, provides a suitable tool to describe magnetic properties at finite temperatures. We can identify $`\mathrm{\Delta }`$ with the energy resulting from a classical Heisenberg model since the corresponding energies in the FM and PM states are given by $`E_{\mathrm{FM}}=x^2_{i0}J_{0i}^{\mathrm{Mn},\mathrm{Mn}}`$ and $`E_{\mathrm{PM}}=0`$, respectively. $`J_{0i}^{\mathrm{Mn},\mathrm{Mn}}`$ denote the exchange interactions between Mn-atoms and all other interactions are neglected. On the other hand, the mean-field theory of the classical Heisenberg model yields $`k_\mathrm{B}\stackrel{~}{T}_c=2x_{i0}J_{0i}^{\mathrm{Mn},\mathrm{Mn}}/3`$, so that we have the following simple estimate,

$$k_\mathrm{B}\stackrel{~}{T}_c=2\mathrm{\Delta }/3x.$$
(1)
It should be noted that non-magnetic As-antisites influence the Curie temperature only indirectly via $`\mathrm{\Delta }`$, which depends on the concentrations $`x`$ and $`y`$ of Mn-atoms and As-antisites, respectively.

### B  Effective on-site exchange parameters

A more systematic approach to the evaluation of $`T_c`$ was developed by Liechtenstein and coworkers  and it is based on the evaluation of magnetic excitation energies from the ground state which are mapped onto the effective Heisenberg Hamiltonian. The Curie temperature is then obtained from its consequent statistical study. A simpler approach, which avoids the direct evaluation of exchange interactions, consists in the determination of the effective exchange field acting on the magnetic moment at this site. From this field the Curie temperature is estimated using the generalized molecular field theory . In the framework of the mean-field approximation adopted here are both approaches equivalent , but the latter one is much simpler when applied to complex systems with few sublattices. Recently, we have applied this approach to evaluate Curie temperatures of bulk ferromagnets  as well as low-dimensional systems such as ultrathin films  and achieved a good agreement with available experimental data.

The basic quantity needed for a determination of $`T_c`$ is the effective on-site exchange parameter for Mn-atoms, $`\overline{J}_i^{\mathrm{Mn}}`$, which represents the Weiss field acting on site $`i`$. In the framework of the TB-LMTO approach we obtain

$$\overline{J}_i^{\mathrm{Mn}}=\frac{1}{4\pi }\mathrm{Im}_C\mathrm{tr}_L\left[\delta _i^{\mathrm{Mn}}(z)(\overline{g}_{}^{\mathrm{Mn},}{}_{ii}{}^{}(z)\overline{g}_{}^{\mathrm{Mn},}{}_{ii}{}^{}(z))+\delta _i^{\mathrm{Mn}}(z)\overline{g}_{}^{\mathrm{Mn},}{}_{ii}{}^{}(z)\delta _i^{\mathrm{Mn}}(z)\overline{g}_{}^{Mn,}{}_{ii}{}^{}(z)\right]dz.$$
(2)
Here $`\mathrm{tr}_L`$ denotes the trace over angular momenta $`L=(\mathrm{}m)`$, the energy integration is performed in the upper half of the complex energy plane over a contour $`C`$ starting below the bottom of the valence band and ending at the Fermi energy, and $`\delta _i^{\mathrm{Mn}}(z)=P_i^{\mathrm{Mn},}(z)P_i^{\mathrm{Mn},}(z)`$, where $`P_i^{\mathrm{Mn},\sigma }(z)`$ are the $`L`$-diagonal potential functions of Mn atoms corresponding to the spin-index $`\sigma `$ ($`\sigma =,`$). The quantity $`\delta _i^{\mathrm{Mn}}(z)`$ is proportional to the corresponding exchange splitting. Finally, the quantities $`\overline{g}_{}^{Mn,\sigma }{}_{ii}{}^{}(z)`$ refer to site-diagonal blocks of the conditionally averaged Green function, namely, the average of the Green function over all configurations with a Mn-atom fixed at the site $`i`$ . The quantity $`\overline{J}_{i}^{}{}_{}{}^{\mathrm{Mn}}`$ reflects the details of the electronic and lattice structures of the system, its carrier concentration, as well as the effect of various kinds of impurities on the resulting electronic structure. The mean-field value of the Curie temperature is given by

$$k_\mathrm{B}T_c=\frac{2}{3}\overline{J}_{i}^{}{}_{}{}^{\mathrm{Mn}}.$$
(3)
The effective on-site exchange parameters corresponding to non-magnetic atoms are at least two-orders of magnitude smaller than $`\overline{J}_{i}^{}{}_{}{}^{\mathrm{Mn}}`$. In Table I we compare $`\overline{T}_c`$ calculated from expressions (1) and (3): both approaches give similar values.

A general trend of Curie temperatures in the FM phase as a function of the chemical composition $`(x,y)`$ is illustrated in Fig. 5. $`T_c`$ increases with $`x`$ in systems without As-antisites and reaches the room temperature at approximately $`x=0.05`$ and then decreases. The same trend and similar values of $`T_c`$ were obtained recently in the framework of the frozen-magnon approach applied to large ordered supercells . The Curie temperature for a fixed Mn-concentration is strongly reduced with increasing concentration of As-antisites. This result clearly indicates a strong correlation between $`T_c`$ and the number of carriers known also from model studies . Detailed results for (Ga1-x-yMnxAsy)As alloys without ($`y=0.0`$) and with ($`y=0.01`$) As-antisites are presented in Fig. 6a. In both cases $`T_c`$ saturates for higher Mn-concentrations which in turn indicates that the effective exchange interactions $`\overline{J}_{0i}^{\mathrm{Mn},\mathrm{Mn}}`$ should decrease with increasing Mn doping. For $`y=0.01`$ the Curie temperature is reduced by approximately 100 K$``$150 K over the whole range of Mn-concentrations as compared to the case of $`y=0.0`$; the results are in a reasonably good agreement with experimental data for ferromagnetic metallic samples ($`0.035<x<0.055`$) .

The concentration of As-antisites $`y`$ in Mn-enriched samples is not well known from the experiment . The detailed knowledge of $`T_c`$ as a function of $`x`$ and $`y`$, Fig. 5, allows to make an estimate of the relation between $`x`$ and $`y`$ in such systems. We have inserted the experimental points into Fig. 5: for each experimental $`T_c`$ corresponding to a given $`x`$ we have determined the concentration $`y`$ such that the calculated and experimental  Curie temperatures coincide. Points obtained in this way follow approximately a straight line. This allows to make an important conclusion, namely that the number of As-antisites increases with the concentration of Mn-atoms. A recent evaluation of the formation energy of the As-antisite defect in GaMnAs  confirms this conclusion.

The dependence of $`T_c`$ on the concentration of As-antisites for a fixed Mn-concentration is presented in Fig. 6b. Full circles correspond to the FM state. There is a monotonic decrease of $`T_c`$ with increasing $`y`$ due to the decreasing number of holes in the valence band which mediate the ferromagnetic coupling between Mn atoms. This result is in a qualitative agreement with predictions of the KEM model . In the present case, however, the concentration of holes is not a free parameter but it is given by the chemical composition ($`x,y`$). We have also determined $`T_c`$ for the pFM state which is the ground state for $`y>0.011`$ (see Fig. 1b). We have now two magnetic atoms Mn+ and Mn-, but the Curie temperature $`T_c`$ can be still estimated from Eq. (3). It should be noted that now the larger of the two values $`\overline{J}_i^{\mathrm{Mn}^\pm }`$ has to be inserted in (3) since the magnetism of such alloy is dominated by this constituent  (empty circles). The transition between the FM and pFM states is manifested by the change of the slope of the curve (see also Fig. 2). The regions of parameters for which the ferromagnetic, partial ferromagnetic, and paramagnetic phases exist are also shown. It should be noted that in the pFM case the Curie temperature is shifted to higher values of $`y`$ than in the FM phase.

### C  Effect of disorder on Curie temperature

Further we investigate the effect of disorder, neglected in the framework of the KEM, on the Curie temperature. For (Ga0.95-yMn0.05Asy)As alloys with constant Mn content is the effect of impurities in the framework of the KEM reduced simply to the change of the number of carriers. The present theory includes the effect of disorder on the electronic structure and hence also on $`T_c`$ in the framework of the CPA. The Curie temperature is influenced by disorder via quantities $`\overline{g}_{}^{\mathrm{Mn},\sigma }{}_{ii}{}^{}(z)`$ in Eq. (2) which depend on the varying concentration of Mn-atoms and As-antisites. We will demonstrate that the effect of disorder on $`T_c`$ is non-negligible.

To this end we have evaluated the Curie temperature as a function of the number of holes in the valence band while keeping the same reference density of states, i.e., by employing the rigid-band model. All quantities in (2) correspond to the reference system but the Fermi energy is an adjustable parameter corresponding to a given alloy composition. The results are summarized in Table II assuming (Ga0.95Mn0.05)As as a reference system and are compared with selfconsistent calculations in which the effect of impurities is fully taken into account. We will first consider the effect of As-antisites assuming (Ga0.95-yMn0.05Asy)As alloys for which the number of holes in the valence band varies as $`n_h=0.052y`$ while the concentration of Mn atoms is kept constant. The selfconsistent calculations and the rigid-band values follow the same trend but the latter ones are systematically smaller and quantitative differences increase with $`y`$: for $`y=0.005`$ and $`y=0.01`$ are the rigid-band values 88% and 62% of the selfconsistent ones, respectively.

It is also possible to use the same approach for estimating $`T_c`$ for (Ga1-xMnx)As alloys without As-antisites. In this case, however, we have to scale the calculated rigid-band values in the spirit of the KEM model by a factor $`x_{\mathrm{ref}}/x`$, where the reference Mn-concentration $`x_{\mathrm{ref}}=0.05`$. The results are similar to the previous case: we observe the same trends in both cases but the rigid-band values are lower than in the case in which the effect of disorder is fully taken into the account. We can conclude that the neglect of disorder does not change the qualitative trends of the calculated $`T_c`$ as a function of Mn- and As-concentrations. On the other hand, the disorder leads in general to an enhancement of $`T_c`$ as compared to the model in which the effect of disorder is treated approximately in the framework of the rigid-band model. This is not surprising as the rigid-band model is not able to describe properly the effect of As- and Mn-impurities which both represent a strong perturbation of the electronic structure.

## IV CONCLUSIONS

Based on first-principles calculations we have investigated the effect of As-antisites on the electronic and magnetic properties of (Ga,Mn)As alloys. In addition to chemical disorder magnetic disorder was described using the disordered local moment model. Both the chemical and the magnetic disorder are treated within the framework of the coherent potential approximation. The mean-field Curie temperatures were estimated from the effective exchange fields acting on a given magnetic atom as well as from total energy differences between the ferromagnetic and paramagnetic state and showed a reasonable agreement between both approaches.

The main results can be summarized as follows: (i) there exists an additional phase with partially ordered local magnetic moments which separates the FM and PM states in the magnetic phase diagram. This phase is stabilized by As-antisites and has a reduced magnetization; (ii) our detailed study of the local magnetic moments justifies a posteriori the validity of the disordered local moment model for determination of magnetic properties of DMS; (iii) a strong reduction of the Curie temperature occurs with increasing concentration of the As-antisites; (iv) with increasing concentration of Mn-atoms the calculations indicate a reduction of the exchange interactions between Mn atoms; (v) a comparison of the calculated and the measured concentration dependences of the Curie temperature indicates a correlation between the concentrations of Mn-impurities and As-antisites, namely an increase of the donor concentration with an increase of the Mn-content; and (vi) a proper account of the influence of impurities on the electronic structure is needed for a quantitative estimate of the Curie temperature, in particular, the Curie temperature is reduced if the effect of disorder on the electronic structure is neglected.

## Acknowledgements

The financial support provided by the Grant Agency of the Czech Republic (No. 202/00/0122), the Grant Agency of the Academy of Sciences of the Czech Republic (No. A1010203, A1010214), the Czech Ministry of Education, Youth, and Sports (COST P5.30, MSM 113200002), the CMS Vienna (GZ 45.504), and the RTN Computational Magnetoelectronics (HPRN-CT-2000-00143) is gratefully acknowledged.