format: text

text:
# Thermodynamic stability of Fe/O solid solution at inner-core conditions

\[

## Abstract

We present a new technique which allows the fully ab initio calculation of the chemical potential of a substitutional impurity in a high-temperature crystal, including harmonic and anharmonic lattice vibrations. The technique uses the combination of thermodynamic integration and reference models developed recently for the ab initio calculation of the free energy of liquids and anharmonic solids. We apply the technique to the case of the substitutional oxygen impurity in h.c.p. iron under Earth’s core conditions, which earlier static ab initio calculations indicated to be thermodynamically very unstable. Our results show that entropic effects arising from the large vibrational amplitude of the oxygen impurity give a major reduction of the oxygen chemical potential, so that oxygen dissolved in h.c.p. iron may be stabilised at concentrations up a few mol % under core conditions.

\]

The thermodynamic stability of oxygen dissolved in iron is a key factor in considering the physics and chemistry of the Earth’s core. We present here a new technique which allows the ab initio calculation of the chemical potential of an impurity in a high-temperature solid solution, including harmonic and anharmonic lattice vibrations. We report the application of the technique to substitutional oxygen dissolved in hexagonal close-packed (h.c.p.) iron at Earth’s core conditions, and we show that the Fe/O solid solution is thermodynamically far more stable than expected from earlier work. The new technique should find wide application to a range of other earth-science problems.

It has long been recognised that the liquid outer core must contain a substantial fraction of light impurities, since its density is $`610`$ % less than that estimated for pure liquid Fe ; similar arguments suggest that the inner core contains a smaller, but still appreciable impurity fraction . The leading impurity candidates are S, Si and O, and arguments have been presented for and against each of them . Ringwood  and others  have argued strongly on grounds of geochemistry that oxygen must account for a large part of the impurity content. However, it has proved difficult to assess these ideas, because the Fe/O phase diagram is so poorly known at Earth’s core conditions. (For reference, we note that the pressures at the core-mantle boundary, the inner-core boundary (ICB) and the centre of the Earth are 136, 330 and 364 GPa respectively; the temperatures at the core-mantle boundary and the ICB are poorly established, but are believed to be in the region of 4000 and 6000 K respectively.)

The thermodynamic stability of dissolved oxygen is governed by the free energy change in the reaction

$$(n1)\mathrm{Fe}(\mathrm{solid})+\mathrm{FeO}(\mathrm{solid})\mathrm{Fe}_n\mathrm{O}(\mathrm{solid}\mathrm{solution})$$
(1)
Let $`\mathrm{\Delta }G`$ be the increase of Gibbs free energy as this reaction goes from left to right, excluding the configurational contribution associated with the randomness of the lattice sites occupied by dissolved O. Then the maximum concentration (number of O atoms per crystal lattice site) at which dissolved O is thermodynamically stable with respect to precipitation of FeO is $`c_{\mathrm{max}}=\mathrm{exp}(\mathrm{\Delta }G/k_\mathrm{B}T)`$. Several years ago, Sherman  used ab initio calculations based on density functional theory (DFT) to calculate the zero-temperature limit of $`\mathrm{\Delta }G`$, i.e. the enthalpy $`\mathrm{\Delta }H`$ of reaction (1). He found that $`\mathrm{\Delta }H`$ is very large ($`5`$ eV at the ICB pressure of 330 GPa), and concluded that the concentration of dissolved O in the inner core must be completely negligible. His argument has been widely cited. However, these were static, zero-temperature calculations, which entirely ignored entropic effects. We shall show here that the high-temperature entropy of dissolved O produces such a large reduction of free energy that Sherman’s argument should be treated with caution when considering core temperatures.

Our ab initio calculations are based on the well established DFT methods used in virtually all ab initio investigations of solid and liquid Fe , including Sherman’s . We employ the generalised gradient approximation for exchange-correlation energy, as formulated by Perdew et al. , which is known to give very accurate results for the low-pressure elastic, vibrational and magnetic properties of body-centred cubic (b.c.c.) Fe, the b.c.c. $``$ h.c.p. transition pressure, and the pressure-volume relation for h.c.p. Fe up to over 300 GPa . We use the ultra-soft pseudopotential implementation  of DFT with plane-wave basis sets, an approach which has been demonstrated to give results for solid Fe that are virtually identical to those of all-electron DFT methods . Our calculations are performed using the VASP code , which is exceptionally stable and efficient for metals. The technical details of pseudopotentials, plane-wave cut-offs, etc. are the same as in our previous work on the Fe/O system .

We first report static zero-temperature results for the enthalpy $`\mathrm{\Delta }H`$ of the reaction (1). Sherman’s results , later confirmed by the present authors , were obtained for the high O concentration of 25 mol %, corresponding to $`n=3`$, but here we wish to focus on the dilute limit, and we take $`n=63`$, which gives a mole fraction of 1.5 %. To do this, we treat a 64-atom supercell with the h.c.p. structure containing a single O substitutional, and we calculate the total ground-state energy and pressure for a sequence of atomic volumes, with all atoms relaxed to their equilibrium positions at each volume. The enthalpy of the pure iron system is obtained from total-energy and pressure calculations for a single unit cell of the h.c.p. crystal. For the FeO crystal, we obtain the enthalpy from total-energy and pressure calculations on a unit cell of the NiAs structure. (The high-pressure stable structure of FeO is believed to be either NiAs or inverse-NiAs; the relative stability of the two structures has been controversial , but our own ab initio calculations indicate that the NiAs structure is slightly more stable at pressures above ca. 145 GPa.) The enthalpy $`\mathrm{\Delta }H`$ is reported as a function of pressure in Fig. 1, where we also show Sherman’s results and our own for the 25 mol % case. We see that $`\mathrm{\Delta }H`$ at 1.5 mol % is between 0.5 and 1.0 eV lower than at 25 mol %, having a value of ca. 4.7 eV at 330 GPa, but this is still very large and Sherman’s arguments would remain valid if this represented a good estimate of the free energy of reaction (1) at Earth’s core temperatures.

To get an idea of the freedom of movement of the substitutional O atom, and hence its vibrational entropy, we now perform a series of calculations in which the O atom is displaced by different amounts from its equilibrium site, with all other atoms held fixed at the equilibrium positions they have when the O atom is at its own equilibrium position. Results for the energy variation with displacement along the $`c`$-axis and in the basal plane for the crystal volume of 6.97 Å3/atom are shown in Fig. 2. Since we are interested in core temperatures, we also mark the energy $`k_\mathrm{B}T`$ for $`T=6000`$ K. We see that in the energy range set by this $`T`$ the energy surface is extremely anharmonic, with almost vanishing curvature at the equilibrium site and large curvatures for large displacements. Clearly, the vibrational root-mean-square (r.m.s.) O displacement will be $`0.45`$ Å (our direct ab initio molecular-dynamics calculations confirm this value). For comparison, we estimate the r.m.s. O displacement in FeO at the same $`P`$ and $`T`$ to be $`0.23`$ Å. This implies a substantial vibrational entropy for the O impurity, because it fits so loosely into the Fe lattice. For the same reason, the vibrations of the 12 Fe neighbours of the O impurity will be softened, and this will also increase the entropy.

To make quantitative statements about high-temperature behaviour, we need to calculate the ab initio Gibbs free energies, rather than zero-temperature enthalpies. The ab initio calculation of the Helmholtz free energies $`F`$ of the Fe and FeO perfect crystals is straightforward in the harmonic approximation, since this requires only the energy of the static lattice and the ab initio lattice vibrational frequencies, which we calculate by the small-displacement method, as discussed in several previous papers . From $`F`$, we then directly obtain the Gibbs free energy $`G=F+PV`$, by calculating the pressure $`P`$ as $`(F/V)_T`$, with $`V`$ the volume. The only difficult part of the present problem is therefore the calculation of $`F`$ (and hence $`G`$) for the $`\mathrm{Fe}_n\mathrm{O}`$ crystal containing the substitutional O atom. This free energy must include the vibrations of many shells of neighbours of the O impurity. The harmonic approximation will clearly not suffice. We meet this challenge by drawing on recently developed ab initio methods for calculating the free energies of liquids and anharmonic solids . These methods rely on two things: empirical reference models, parameterised to accurately mimic the ab initio energies; and the technique of ‘thermodynamic integration’, used to determine free energy differences. Our overall strategy will be to obtain the ab initio free energy $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}`$ of the O-substitutional system by starting from the ab initio free energy $`F_{\mathrm{Fe}}^{\mathrm{AI}}`$ of pure Fe and using thermodynamic integration to compute the free energy change $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}}^{\mathrm{AI}}`$ that results from converting a single Fe atom into an O atom.

We recall briefly that thermodynamic integration  is a general technique for calculating the difference of free energies $`F_1F_0`$ of two systems containing the same number $`N`$ of atoms but having different total-energy functions $`U_0(𝐫_1,𝐫_2,\mathrm{}𝐫_N)`$ and $`U_1(𝐫_1,𝐫_2,\mathrm{}𝐫_N)`$, with $`𝐫_i(i=1,2\mathrm{}N)`$ the atomic positions. The technique relies on the equivalence between the free energy difference and the reversible work done on switching the total energy function continuously from $`U_0`$ to $`U_1`$. The work done is:

$$F_1F_0=_0^1𝑑\lambda U_1U_0_\lambda ,$$
(2)
where the thermal average $`_\lambda `$ is evaluated in the canonical ensemble generated by the switched energy function $`U_\lambda `$ defined by:

$$U_\lambda =(1\lambda )U_0+\lambda U_1.$$
(3)
To apply this in practice, we use molecular dynamics simulation to evaluate the average $`U_1U_0_\lambda `$ at a sequence of $`\lambda `$ values and we perform the integration over $`\lambda `$ numerically.

In principle, we could calculate $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}}^{\mathrm{AI}}`$ by identifying $`U_0`$ and $`U_1`$ as the ab initio total energy functions $`U_{\mathrm{Fe}}^{\mathrm{AI}}`$ and $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}`$ of the pure-Fe and O-substitutional systems, but this brute-force approach is computationally prohibitive at present. It is also unnecessary, since exactly the same result can be achieved much more cheaply by using empirical reference models. In our recent work on liquid Fe , we have found that a simple inverse-power pair potential $`\varphi (r)=A/r^\alpha `$ reproduces the ab initio total energy very accurately; for the anharmonic high-temperature Fe crystal, a linear combination of this pair-potential model with an ab initio harmonic description has been very effective . We denote the total energy of this latter anharmonic model by $`U_{\mathrm{Fe}}^{\mathrm{ref}}(𝐫_1,\mathrm{}𝐫_N)`$, where $`𝐫_i`$ are the atomic positions.

In order to make a reference system for the Fe/O system containing a single substitutional O atom, whose total energy is $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}(𝐫_1,\mathrm{}𝐫_N)`$ ($`𝐫_1`$ is the position of the O atom), we simply delete all terms in $`U_{\mathrm{Fe}}^{\mathrm{ref}}`$ involving atom 1 and replace them with a pair interaction potential $`\chi (r)=B/r^\beta `$. All parts of $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$ involving the $`N1`$ Fe atoms $`2,3,\mathrm{}N`$ remain precisely as they are in $`U_{\mathrm{Fe}}^{\mathrm{ref}}`$. Our procedure for determining $`B`$ and $`\beta `$ starts by requiring that the pair potential $`\chi (r)`$ should reproduce as well as possible the dependence of the total energy on position of the O atom with all Fe atoms held fixed, i.e. the curves shown in Fig. 2. The values for $`B`$ and $`\beta `$ thus obtained give us an initial form for $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$. This initial $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$ is then used in a classical MD simulation to generate a long trajectory at the temperature of interest, from which we take 100 statistically independent configurations. The full ab initio energies are calculated for these configurations, and the $`B`$ and $`\beta `$ parameters are readjusted to give a least squares fit to these energies. Finally, the $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$ obtained from the new $`B`$, $`\beta `$ is used to generate a further 100 statistically independent configurations, and $`B`$ and $`\beta `$ are adjusted once more to fit the ab initio energies of these configurations. The $`B`$ and $`\beta `$ produced by this final step are found to be essentially identical to those in the previous step, and we accept them as the optimal values for the assumed form of $`\chi (r)`$.

The free energy difference $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}}^{\mathrm{AI}}`$ between the Fe/O and pure Fe systems is now expressed as the sum of three contributions: the difference $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}F_{\mathrm{Fe}}^{\mathrm{ref}}`$ between the reference systems, and the two differences $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$ and $`F_{\mathrm{Fe}}^{\mathrm{ref}}F_{\mathrm{Fe}}^{\mathrm{AI}}`$ between the ab initio and reference systems. All three of these differences are calculated by thermodynamic integration. We emphasise that, although the reference systems play a vital role, the final result for $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}}^{\mathrm{AI}}`$ does not depend on how they are chosen.

Thermodynamic integration to get $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}F_{\mathrm{Fe}}^{\mathrm{ref}}`$ is easy and rapid, since only simple potential models are involved. We have tested size effects by using simulated systems containing up to 768 atoms, but we find that with only 288 atoms the size-effect errors are less than the statistical errors of ca. 30 meV. For the ab-initio/reference differences $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$ and $`F_{\mathrm{Fe}}^{\mathrm{AI}}F_{\mathrm{Fe}}^{\mathrm{ref}}`$, the fluctuations of $`U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}U_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{ref}}`$ and $`U_{\mathrm{Fe}}^{\mathrm{AI}}U_{\mathrm{Fe}}^{\mathrm{ref}}`$ are so small that explicit thermodynamic integration over $`\lambda `$ is unnecessary, and we can use instead the small-$`\lambda `$ approximation explained elsewhere . We have studied the size errors for these ab initio/reference differences, and we find that results obtained with 36, 64 and 96 atoms are identical within statistical errors. The overall statistical error on the ab initio difference $`F_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}F_{\mathrm{Fe}}^{\mathrm{AI}}`$ is ca. 90 meV. We have repeated all the above calculations at the four volumes 6.86, 6.97, 7.20 and 7.40 Å3/atom, and from the dependence on volume we obtain the pressure change on replacing Fe by O and hence the Gibbs free energy difference $`G_{\mathrm{Fe}/\mathrm{O}}^{\mathrm{AI}}G_{\mathrm{Fe}}^{\mathrm{AI}}`$.

Our calculated Gibbs free energies $`\mathrm{\Delta }G`$ for reaction (1) are displayed in Fig. 1 for the two temperatures 5000 and 7000 K. We note the very large entropic lowering of $`\mathrm{\Delta }G`$, which, at $`P=330`$ GPa comes down from 4.7 to ca. 1.7 eV at the temperature $`T6000`$ K expected at the ICB. This is still a substantial positive value, but implies that the stability-limit concentration $`c_{\mathrm{max}}=\mathrm{exp}(\mathrm{\Delta }G/k_\mathrm{B}T)`$ is ca. 3 mol %, which is far from negligible.

In assessing our $`c_{\mathrm{max}}`$ value, one should note the remaining uncertainties in our calculations. First, we have ignored anharmonicity in the pure Fe and FeO crystals. Our recent work on the effect of anharmonicity in pure Fe  showed that at the melting point, anharmonicity can contribute as much as 70 meV/atom to the free energy; the same might be true of FeO. These effects could shift $`\mathrm{\Delta }G`$ by perhaps 0.15 eV. Second, there is the question of strong electronic correlation in FeO, which is a prototypical Mott insulator at low pressures. Such correlation effects will be much weakened at Earth’s core pressures, but could still shift $`\mathrm{\Delta }G`$ by a few tenths of an eV. This means that our prediction for $`c_{\mathrm{max}}`$ at a given temperature is probably not reliable to better than a factor of 3. We are therefore cautious about the detailed numerical value of $`c_{\mathrm{max}}`$, and claim only that it could be a few mol % at ICB pressure and temperature.

In summary, we conclude that, because substitutional oxygen atom fits so loosely into the Fe lattice and has so much freedom of movement, it undergoes a very large entropic lowering of free energy at high temperatures, this lowering being as much as 3 eV at 6000 K and 330 GPa. The consequence is that substitutional O dissolved in h.c.p. Fe may be thermodynamically stabilised at concentrations up to a few mol %. Earlier ab initio calculations  which ignored entropic effects should therefore not be taken at face value. Finally, we point out that a wide range of geological problems depend on an understanding of chemical potentials – for example, all problems concerned with the partitioning of elements between coexisting phases. The ab initio techniques for calculating chemical potentials outlined here should therefore be of wide interest.

Acknowledgments. The work of D.A. was supported by NERC Grant No. GR3/12083. Allocations of time on the Cray T3E machines at the Manchester CSAR service and at Edinburgh Parallel Computer Centre were provided by the UK Car-Parrinello Consortium (EPSRC grant GR/M01753) and the Minerals Physics Consortium (NERC grant GST/02/1002). Some of the calculations were performed at the UCL HiPerSPACE Centre, partially funded by the Joint Research Equipment Initiative.


text:
# 1 Introduction

## 1 Introduction

Recently there have been proposals to reconstruct the quantum state of electromagnetic fields inside cavities . The reconstruction of non- classical states is a central topic in quantum optics and related fields, and there have been a number of proposals to achieve it (see for instance ). In fact, the full reconstruction of nonclassical field states  as well as of (motional) states of an ion  have been experimentally accomplished. The reconstruction is normally achieved through a finite set of either field homodyne measurements, or selective measurements of atomic states .

However, the presence of noise and dissipation has normally destructive effects. In fact, the reconstruction schemes themeselves also indicate loss of coherence in quantum systems . Schemes for compensation of losses have already been proposed  and the relation between losses and $`s`$-parametrized quasiprobabilities has been pointed out in ref . A method of reconstruction of the Wigner function that takes into account losses (at $`T=0`$) has also been presented .

We consider here a single mode high-$`Q`$ cavity where we suppose that a nonclassical field is prepared. The first step of our method consists in driving the generated state by a coherent pulse. The reconstruction of the field is done after turning-off the driving field, i.e. at a time when the cavity field has interacted with its environment at finite temperature. We show that by measuring the density matrix diagonal elements and properly weighting them, we can obtain directly the Wigner function even at $`T0`$. We should remark that to know a state, one has to have information about all the density matrix elements (diagonal and off-diagonal), however, with the method presented here (see also ), it is only necessary to have information about diagonal matrix elements.

## 2 Master equation and its solution

The master equation in the interaction picture for the reduced density operator $`\widehat{\rho }`$ relative to a driven cavity mode, taking into account cavity losses at non-zero temperature and under the Born-Markov approximation is given by  (in a frame rotating at the field frequency $`\omega `$)

$$\frac{\widehat{\rho }}{t}=(\widehat{}+\widehat{})\widehat{\rho },$$
(1)
where

$$\widehat{}\widehat{\rho }=(\widehat{}_1+\widehat{}_2)\widehat{\rho }$$
(2)
with

$$\widehat{}_1\widehat{\rho }=\frac{\gamma (\overline{n}+1)}{2}\left(2\widehat{a}\widehat{\rho }\widehat{a}^{}\widehat{a}^{}\widehat{a}\widehat{\rho }\widehat{\rho }\widehat{a}^{}\widehat{a}\right),\widehat{}_2\widehat{\rho }=\frac{\gamma \overline{n}}{2}\left(2\widehat{a}^{}\widehat{\rho }\widehat{a}\widehat{a}\widehat{a}^{}\widehat{\rho }\widehat{\rho }\widehat{a}\widehat{a}^{}\right),$$
(3)
and

$$\widehat{}\widehat{\rho }=\frac{i}{\mathrm{}}[\widehat{H},\widehat{\rho }],$$
(4)
where

$$\widehat{H}=i\mathrm{}\left(\alpha ^{}\widehat{a}\alpha \widehat{a}^{}\right).$$
(5)
$`\widehat{a}`$ and $`\widehat{a}^{}`$ are the annihilation and creation operators, $`\gamma `$ the (cavity) decay constant, $`\overline{n}`$ is the mean number of thermal photons and $`\alpha `$ the amplitude of the driving field.

### 2.1 Displacing the field

The formal solution to Eq. (1) is given by (see for instance )

$$\widehat{\rho }(t)=\mathrm{exp}\left[(\widehat{}+\widehat{})t\right]\widehat{\rho }(0).$$
(6)
It is not difficult to show that Eq. (6) can be factorized in the product of two exponentials, one containing the reservoir (super) operators and the other the interaction (5), the latter one yielding an effective displacement on the initial field. In order to show this we calculate the commutator

$$[\widehat{},\widehat{}]\widehat{\rho }=\frac{\gamma }{2}\widehat{}\widehat{\rho },$$
(7)
which allows the factorization

$$\widehat{\rho }(t)=\mathrm{exp}(\widehat{}t)\mathrm{exp}\left[\frac{2\widehat{}}{\gamma }(1e^{\gamma t/2})\right]\widehat{\rho }(0).$$
(8)
After driving the initial field during a time $`t`$, the resulting field density operator will read

$$\widehat{\rho }(t)=e^{\widehat{}t}\widehat{\rho }_\beta (0),$$
(9)
where

$$\widehat{\rho }_\beta (0)=\widehat{D}^{}(\beta )\widehat{\rho }(0)\widehat{D}(\beta ),$$
(10)
with the effective (displacing) amplitude

$$\beta =2\alpha \frac{1e^{\gamma t/2}}{\gamma }.$$
(11)
### 2.2 Solution to the Master Equation at finite temperature

We now obtain the density matrix (4), by defining

$$\widehat{J}_{}\widehat{\rho }=\widehat{a}\widehat{\rho }\widehat{a}^{},\widehat{J}_+\widehat{\rho }=\widehat{a}^{}\widehat{\rho }\widehat{a},\widehat{J}_3\widehat{\rho }=\widehat{a}^{}\widehat{a}\widehat{\rho }+\widehat{\rho }\widehat{a}^{}\widehat{a}+\widehat{\rho },$$
(12)
where the superoperators $`\widehat{J}_{}`$, $`\widehat{J}_+`$ and $`\widehat{J}_3`$ obey the commutation relations $`[\widehat{J}_{},\widehat{J}_+]\widehat{\rho }=\widehat{J}_3\widehat{\rho }`$ and $`[\widehat{J}_3,\widehat{J}_\pm ]\widehat{\rho }=\pm 2\widehat{J}_\pm \widehat{\rho }`$, which can be written as

$$\widehat{\rho }(t)=e^{\frac{\gamma t}{2}}e^{\mathrm{\Gamma }_{\overline{n}}(t)\widehat{J}_+}\left[\frac{e^{\gamma t/2}}{1+N_t}\right]^{\widehat{J}_3}e^{\mathrm{\Gamma }_{\overline{n}+1}(t)\widehat{J}_{}}\widehat{\rho }_\beta (0),$$
(13)
and we have defined

$$\mathrm{\Gamma }_{\overline{n}}(t)=\frac{\overline{n}(1e^{\gamma t})}{1+N_t},\mathrm{\Gamma }_{\overline{n}+1}(t)=\frac{(\overline{n}+1)(1e^{\gamma t})}{1+N_t},$$
(14)
with $`N_t=\overline{n}(1e^{\gamma t})`$.

## 3 The reconstruction method

We now calculate the diagonal density matrix elements $`<m|\widehat{\rho }(t)|m>`$ from (13)

$$m|\widehat{\rho }_\beta (t)|m=\frac{1}{1+N_t}\underset{k=0}{\overset{\mathrm{}}{}}\underset{n=0}{\overset{\mathrm{}}{}}\left(\begin{array}{c}k\\ n\end{array}\right)\left(\begin{array}{c}m\\ n\end{array}\right)[\mathrm{\Gamma }_{\overline{n}}(t)]^{mn}[\mathrm{\Gamma }_{\overline{n}+1}(t)]^{kn}\frac{e^{n\gamma t}}{[1+N_t]^{2n}}k|\widehat{\rho }_\beta (0)|k.$$
(15)
Multipliying the expression above by powers of the weight function $`\chi _s`$, where

$$\chi _s=\frac{\frac{s+1}{s1}\mathrm{\Gamma }_{\overline{n}+1}(t)}{\frac{e^{\gamma t}}{[1+N_t]^2}+\mathrm{\Gamma }_{\overline{n}}(t)\left(\frac{s+1}{s1}\mathrm{\Gamma }_{\overline{n}+1}(t)\right)}$$
(16)
and adding over $`m`$ we obtain

$$F(\beta ;s)=\underset{m=0}{\overset{\mathrm{}}{}}\chi _s^mm|\widehat{\rho }_\beta (t)|m=\frac{1}{[1+N_t][1\chi _s\mathrm{\Gamma }_{\overline{n}}(t)]}\underset{k=0}{\overset{\mathrm{}}{}}\left(\frac{s+1}{s1}\right)^kk|\widehat{\rho }_\beta (0)|k,$$
(17)
where we have used the fact that

$$\underset{m=0}{\overset{\mathrm{}}{}}\left(\begin{array}{c}m\\ n\end{array}\right)x^m=\frac{x^n}{[1x]^{n+1}}.$$
(18)
If we now multyply $`F(\beta ;s)`$ by the quantity

$$\frac{2[1+N_t][1\chi _s\mathrm{\Gamma }_{\overline{n}}(t)]}{\pi (s1)},$$
(19)
we finally obtain

$$W(\beta ;s)=\frac{2}{\pi (s1)}\underset{k=0}{\overset{\mathrm{}}{}}\left(\frac{s+1}{s1}\right)^kk|\widehat{\rho }_\beta (0)|k,$$
(20)
which is the $`s`$-parametrized quasiprobability distribution .

Therefore, by measuring the diagonal elements of the evolved field density matrix, Eq. (15), we may obtain complete information on the initial state.

We should remark that in this case (thermal environment), more information on the $`P_m(t)=<m|\widehat{\rho }(t)|m>`$ is required than in the zero temperature case. Nevertheless, it is possible to find a weight function that allows the reconstruction of the initial field.

### 3.1 Measuring the photon distribution

For completness, we suggest a way to measure the photon number distribution $`P_m(\beta ,t)`$ of Eq. (17). It is not difficult to show that the atomic inversion for the case of a three-level atom in a cascade configuration with the upper and the lower levels having the same parity and satisfaying the two-photon resonance condition is given by (see for instance )

$$W(\beta ;t+\tau )\underset{m=0}{\overset{\mathrm{}}{}}P_m(\beta ,t)\mathrm{cos}([2m+3]\lambda \tau ),$$
(21)
where $`\lambda `$ is the atom-field coupling constant. In order to obtain $`P_m`$ from a family of measured population inversions, we invert the Fourier series in Eq. (21), or

$$P_m(\beta ,t)=\frac{2\lambda }{\pi }_0^{\tau _{max}}𝑑\tau W(t+\tau )\mathrm{cos}([2m+3]\lambda \tau ).$$
(22)
We need a maximum interaction time $`\tau _{max}=\pi /\lambda `$ much shorter than the cavity decay time, which implies we must be in the strong-coupling regime, i.e. $`\lambda >>\gamma `$.

## 4 Conclusions

We have presented a method to reconstruct the Wigner function (and in general any quasiprobability distribution) of an initial nonclassical state at times when the field would have normally lost its quantum coherence because of the interaction with an environment at finite temperature. This is an extension of our previous work  where we considered the interaction with the environment at zero temperature. The crucial point of our approach is the apropriate weighting of the evolved (driven and decayed) photon number distribution. Driving the initial field immediately after preparation, is not only useful for covering a region in phase space but also makes possible, together with the weight function, $`\chi _s`$, to store quantum coherences in the diagonal elements of the time evolved density matrix. We have shown here that it is possible to find a weight function that allows to reconstruct the initial field even at finite temperature, and this is the main result of our paper. Similar conclusions may be reached by employing the method of generating functions . The possibility of reconstructing quantum states when the system interacts with an environment maybe relevant for applications in quantum computing. Loss of coherence for such interactions is likely to occur in such devices, and our method could be used, for instance, as a scheme to refresh the state of a quantum computer  in order to minimize the destructive action of a hot environment.

Acknowledgments

This work was partially supported by CONACYT (Mexico), FAPESP and CNPq (Brazil).


text:
# Fast Quantum Search Algorithms in Protein Sequence Comparison - Quantum Biocomputing

## Abstract

Quantum search algorithms are considered in the context of protein sequence comparison in biocomputing. Given a sample protein sequence of length $`m`$ (i.e $`m`$ residues), the problem considered is to find an optimal match in a large database containing $`N`$ residues. Initially, Grover’s quantum search algorithm is applied to a simple illustrative case - namely where the database forms a complete set of states over the $`2^m`$ basis states of a $`m`$ qubit register, and thus is known to contain the exact sequence of interest. This example demonstrates explicitly the typical $`O(\sqrt{N})`$ speedup on the classical $`O(N)`$ requirements. An algorithm is then presented for the (more realistic) case where the database may contain repeat sequences, and may not necessarily contain an exact match to the sample sequence. In terms of minimizing the Hamming distance between the sample sequence and the database subsequences the algorithm finds an optimal alignment, in $`O(\sqrt{N})`$ steps, by employing an extension of Grover’s algorithm, due to Boyer, Brassard, H$`ø`$yer and Tapp for the case when the number of matches is not a priori known.

The fantastic possibilities of quantum parallelism in computing, suggested by the convergence of quantum mechanics and information theory in the past two decades, are fast being enumerated in the guise of quantum algorithms. First, and foremost, among these is the factoring algorithm of Shor, which provided great impetus to the field of quantum computing. Shor’s algorithm applied to a given number $`N`$ requires $`O((\mathrm{log}N)^3)`$ steps, and represents an exponential speed-up over the best classical algorithms. Another important result, due to Grover, was the discovery of a quantum search algorithm for finding a particular element in an unordered set of $`N`$ elements in only $`O(\sqrt{N})`$ steps - a significant improvement over the classical cost $`O(N)`$.

In this paper the application of quantum search algorithms to an important problem at the heart of biocomputing (or bioinformatics), that of protein sequence comparison and alignment, is considered. As the mapping and sequencing of the human genome (some $`3\times 10^9`$ base pairs) nears completion, the relatively new field of biocomputing has become obvious in its importance to the quantitative analysis of this vast amount of data. Some fundamental tasks in biocomputing involving sequence analysis include: searching databases in order to compare a new sub-sequence to existing sequences, inferring protein sequence from DNA sequence, and calculation of sequence alignment in the analysis of protein structure and function. A tremendous amount of computing is required, much of which is devoted to search-type problems, either directly in large databases, or in configuration space of alignment possibilities. While it is possible that all of these problems may be amenable to quantum algorithmic speed-up, it is explicitly demonstrated in this work how the fundamental task of sequence alignment can be approached using a quantum computer. Indeed, this problem is a very natural application of the quantum search algorithm (perhaps a strange reflection of the possibility that the machinery of DNA itself may actually function using quantum search algorithms).

In general terms Grover’s search algorithm relies on the existence of a quantum computer $`Q`$ operating using an oracle function, F. The set of search possibilities is represented by states in the Hilbert space of $`Q`$. The oracle function simply tests whether a given state is the actual target state. Grover found a unitary operator $`U`$ (involving the oracle function test) which evolves the quantum computer in such a way that the amplitude of the target state in the wave function of $`Q`$ is amplified. Furthermore, Grover showed that there exists a number $`k<\sqrt{N}`$, such that after $`k`$ applications of $`U`$, the probability of finding the target state is at least 1/2. Subsequently, Boyer, Brassard, H$`ø`$yer and Tapp (BBHT) proved a tighter bound: one must iterate the algorithm on average at least $`\left(\mathrm{sin}\frac{\pi }{8}\right)\sqrt{N}`$ times to achieve a probability of 1/2 for finding the target.

To begin the application of quantum search algorithms to protein sequence analysis, the problem of sequence alignment to a large database of sequence domains is considered. That is, given a sample sequence the task is to find out the location in the database of an exact or closest match (with respect to some defined measure). Application of the Grover algorithm directly to this search task would cause trouble immediately because, by definition, it is not known if the target exists in the database, or if it actually exists multiple times. If there are actually $`N_t`$ solutions, the number of iterations required to find a solution with probability 1/2 is $`\left(\mathrm{sin}\frac{\pi }{8}\right)\sqrt{N/N_t}`$. Thus, if one does not know the number of solutions at the outset, the computer may inadvertently be halted when the amplitude of the target states is very small. This happens because the process of amplitude amplification is not monotonic, but rather oscillates with the number of iterations. Fortunately, this difficult impasse has been solved by BBHT and they provide an algorithm, based on Grover’s algorithm as a subroutine, for finding a solution in the case where the number of solutions is unknown. This result allows for the application of quantum search algorithms to the field of biocomputing.

In terms of protein sequences, the human genome is composed of about 150,000 domains, each containing on average 300 residues (amino acids). An interesting feature of approaching the sequence analysis problem using a quantum computer is that the entire database could in principle be stored in a single wave function superposition, and then be presented simultaneously for inspection. To illustrate the basic idea, a very simple case of sequence comparison is initially considered, followed by a more realistic problem later. Consider a database, $`D`$, constructed from the domains of the human genome placed end-to-end, so that a continuous list of $`N`$ residues, $`D=\{R_0,R_1,\mathrm{},R_{N1}\}`$, is created. Independently, a sample sequence is given $`s=\{r_0,r_1,\mathrm{},r_{m1}\}`$ composed of $`m`$ residues; the task is to compare this with the database. Each residue is labeled by a letter of the 20-letter amino acid alphabet, so in order to encode the database 5 bits per residue is needed. Thus, the residues $`R_i`$ and $`r_i`$ are represented by bit strings, $`_{\alpha =0}^4B_{i\alpha }`$ and $`_{\alpha =0}^4b_{i\alpha }`$, respectively.

The quantum computer to analyze this system is composed of two registers, with number of qubits $`Q_1`$ and $`Q_2`$, respectively. The bit-wise representation of the protein sequences will be encoded into the qubits of this system. Leaving issues of data transfer aside, the entire database is represented by a quantum superposition over the two registers:

$$|\mathrm{\Psi }_D\frac{1}{\sqrt{Nm+1}}\underset{i=0}{\overset{Nm}{}}|\varphi _i|i,$$
(1)
where all the consecutive sub-sequences in the database of length $`m`$ are encoded in the first register with $`Q_1=5m`$ as

$$|\varphi _i=\underset{\alpha =i}{\overset{i+m1}{}}\underset{\beta =0}{\overset{4}{}}|B_{\alpha \beta }\underset{\alpha =0}{\overset{5m1}{}}|q_{i\alpha }.$$
(2)
That is, from from the database of length $`N`$ residues, $`Nm+1`$ sub-sequences of length $`m`$ are constructed by moving along from the first position (allowing domain crossing). Position information of the sub-sequences is meanwhile tagged explicitly by binary numbers, $`|i`$, in the second register, and is accessed by an operator, $`\widehat{X}`$, acting in the Hilbert space of the second register, which gives the position as $`\widehat{X}|i=i|i`$ ($`0iNm`$). In order that this register can encode all positions $`Q_2`$ must satisfy $`2^{Q_2}>Nm`$. The number of qubits required in this register is relatively small: taking the database size to be that for the number of residues in the human genome implies $`Q_2=26`$ suffices. In the first register, typical sequence comparison problems require $`mO(300)`$.

The next step in the initialization process is the coding of a table, $`T[0\mathrm{}Nm]`$, into the quantum state, which measures the difference between the database states $`|\varphi _i`$ and the sample sequence state in terms on the total number of bit flips required to transform any database state into the sample sequence. In other words, $`T[0\mathrm{}Nm]`$ is the set of Hamming distances. Remarkably, the set of Hamming distances for the entire database can be created by simply acting on each qubit of the computer with a CNOT operation with respect to the sample sequence state:

$$|\mathrm{\Psi }_H=U_{\mathrm{CNOT}}(s)|\mathrm{\Psi }_D\frac{1}{\sqrt{Nm+1}}\underset{i=0}{\overset{Nm}{}}|\overline{\varphi _i}|i.$$
(3)
Denoting the individual qubits of the “Hamming states” $`|\overline{\varphi _i}`$ by

$$|\overline{\varphi _i}=\underset{\alpha =0}{\overset{5m1}{}}|\overline{q}_{i\alpha },$$
(4)
an operator, $`\widehat{T}`$, is introduced which, acting on a state $`|\overline{\varphi _i}`$, gives the Hamming distance table value $`T[i]`$ as:

$$\widehat{T}:\widehat{T}|\overline{\varphi _i}=T[i]|\overline{\varphi _i},T[i]=\underset{\alpha =0}{\overset{5m1}{}}\overline{q}_{i\alpha }.$$
(5)
With the computer design completed and initialized, a simple search problem can be defined in order to demonstrate how the computer works. First, the database is taken to be of length $`N=2^m+m1`$ so that there are exactly $`2^m`$ states in the superposition, and furthermore demand that all these states are distinct. The problem is to search the database for the sub-sequence $`s`$, which occurs exactly once, but at an unknown location. Classically, this would require $`O(N)`$ steps. However, by using Grover’s search algorithm, the match can be found in $`O(\sqrt{N})`$ steps. In this example, the database decomposition has been artificially arranged to be over a complete set of states of the first register, which means that Grover’s search algorithm can be applied directly.

The problem defined by Grover has been modified slightly, but the applicability of the search algorithm remains. The original problem was defined in terms of an oracle function, $`F(x)`$, over a set of values $`x\{0,\mathrm{},N1\}`$, which is zero everywhere except at some value $`t`$, the target of the search, where $`F(t)=1`$. The sequence comparison problem here has been re-structured so that a value of $`x`$ represents a subsequence of the database, and the oracle function is just a direct comparison with the sample sequence. In a sense, the black box nature of the oracle function has been simplified, at the cost of increasing the complexity of the initial wave function with position information. It remains to be seen whether this is a feasible way of coding a sequence database. Of course, an alternative is to sweep all details of the database look-up and comparison into the oracle function. The difference is subtle, and perhaps non-trivial in practice. The advantage of the latter approach might be in the initialization of the quantum computer state. The algorithms presented here would still apply in this case.

In the computer design defined here, Grover’s search algorithm is applied to the first register containing the sub-sequence state superposition. The problem is to find the state $`|\overline{s}=U_{\mathrm{CNOT}}|s=|0\mathrm{}0`$ (zeros in all $`m`$ qubits of the first register) with table value $`T[i_s]=0`$, occurring at position $`i_s`$ (as yet unknown). Once the state is found, the location of the sequence in the database can be determined by making a measurement of $`\widehat{X}`$ on the second register.

To illustrate the working of the algorithm the geometrical picture, which is particularly transparent, is applied to this framework. The search algorithm is initiated by decomposing the state $`|\mathrm{\Psi }_H`$ into orthogonal components with respect to $`|\overline{s}`$ as

$$|\mathrm{\Psi }_H=\sqrt{\frac{Nm}{Nm+1}}|R+\frac{1}{\sqrt{Nm+1}}|S,$$
(6)
where

$`|S`$ $`=`$ $`|\overline{s}|i_s`$ (7)
$`|R`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{Nm}}}{\displaystyle \underset{ii_s}{}}|\overline{\varphi _i}|i.`$ (8)
The evolution of the quantum computer representing the search algorithm occurs in the first register, the second register lying dormant, yet through quantum entanglement carrying the position information required at the end. The operator $`U`$ is constructed from reflection operators in the Hilbert space of the first register,

$`I_S`$ $`=`$ $`12|SS|`$ (9)
$`I_H`$ $`=`$ $`12|\mathrm{\Psi }_H\mathrm{\Psi }_H|.`$ (10)
The operator $`I_S`$ contains the query to the oracle function, $`F(i)`$, and acts on the Hamming states $`|\overline{\varphi _i}`$ with a phase shift dependent on the search criteria $`T[i_s]=0`$:

$$I_S|\overline{\varphi _i}=(1)^{F(i)}|\overline{\varphi _i}=\{\begin{array}{ccc}|\overline{\varphi _i}\hfill & & \mathrm{if}T[i]=0\hfill \\ & & \\ |\overline{\varphi _i}\hfill & & \mathrm{otherwise}.\hfill \end{array}$$
(11)
In terms of these reflection operators, the unitary operator evolving the system through one step of the search algorithm is given by $`U=I_HI_S`$. The evolution of the computer proceeds through application of the operator, $`U`$, a number of times on the initial state, $`|\mathrm{\Psi }_H`$. The effect of this evolution is to amplify the component of the target state, $`|S`$ in the superposition. It is important to understand the nature of this process in order to appreciate how the quantum computer functions. To see this point it is convenient to express $`U`$ in the representation of the subspace $`\{|S,|R\}`$:

$`U=\left[\begin{array}{ccc}\frac{Nm1}{Nm+1}& & \frac{2\sqrt{Nm}}{Nm+1}\\ & & \\ \frac{2\sqrt{Nm}}{Nm+1}& & \frac{Nm1}{Nm+1}\end{array}\right]=\left[\begin{array}{ccc}\mathrm{cos}\theta & & \mathrm{sin}\theta \\ & & \\ \mathrm{sin}\theta & & \mathrm{cos}\theta \end{array}\right].`$ (18)
where $`\mathrm{sin}\theta 2\sqrt{Nm}/(Nm+1)`$.

After $`k`$ steps of the algorithm the state of the computer is given by

$$|\mathrm{\Psi }_k=U^k|\mathrm{\Psi }_H=\underset{i=0}{\overset{Nm}{}}c_i^{(k)}|\overline{\varphi _i}|i.$$
(19)
The amplitude of the target state, $`c_{i_s}^{(k)}`$, can be easily calculated using the matrix representation for $`U`$. One obtains:

$$c_{i_s}^{(k)}=\mathrm{cos}\left(k\theta \alpha \right),\mathrm{cos}\alpha \frac{1}{\sqrt{Nm+1}}.$$
(20)
The component along $`|S`$ is amplified to near unity at $`k_{\mathrm{max}}\frac{\pi }{4}\sqrt{N}`$ (for $`N>>m`$). A measurement of $`\widehat{T}`$ on the first register will give a result $`T[i]`$ with probability $`|c_i^{(k)}|^2`$. If $`T[i]=0`$ then the algorithm has succeeded - i.e the sample sequence has been found - and a subsequent measurement of $`\widehat{X}`$ in the second register will give the position, $`i_s`$, of the sequence in the database. A crucial point is that one has to be careful interpreting the number of steps required to obtain a successful outcome - merely increasing the number of steps beyond $`k_{\mathrm{max}}`$ does not improve the chances of success because the amplification is not monotonic. Indeed, the probability of success actually decreases when $`k_{\mathrm{max}}`$ is exceeded. The search may therefore have to be run several times, however, for large $`N`$ the savings in computer time compared to a classical computer are clear, even if the search is repeated several times.

While the above example serves to display the potential of quantum search algorithms in the context of sequence matching to a large database, it does not contain an important concept in bioinformatics - optimal alignment. Generally, the sample sequence may not be contained exactly in the database, and so one is interested in how close is the best match (or matches) with respect to a well defined distance measure. Often this measure involves editing of strings by insertion of gaps in order to minimize the distance; in practice this process is very complicated. In the first instance, the problem is extended to that of finding an optimal alignment with respect to the Hamming distance, without editing of sequences (which can be incorporated at a later stage).

Let us first define the problem using, as far as possible, the same notation as previously. The database is taken to be of size $`N>>m`$, but the restriction that the set of database sub-sequence states is equal to $`2^m`$ is relaxed, and the possibility is allowed that the set of sub-sequences may contain repeats, and, more importantly, may or may not contain the sample sequence. The problem then is to find an optimal alignment of the sample sequence to a sub-sequence in the database. An optimal alignment here is defined in the sense of finding the smallest Hamming distance $`T[i]`$ with respect to the sample sequence state.

In terms of our quantum computer, the database state in this case is also described by the state $`|\mathrm{\Psi }_D`$. An important point is that the state is still normalised by the factor $`1/\sqrt{Nm+1}`$ because the repeats occur at different locations, and thus each state in the product space of the two registers is distinct. The introduction of the position register $`Q_2`$ has ensured this. Using the CNOT operation on $`|\mathrm{\Psi }_D`$ the superposition, $`|\mathrm{\Psi }_H`$, of Hamming states is once again obtained. The algorithm strategy is to search for alignments of increasing Hamming distance. At the start of each search it is not known how many solutions exist, or if there exist matches at all, and so Grover’s algorithm cannot be used directly. However, we use now the extension of Grover’s algorithm due to Boyer, Brassard, H$`ø`$yer and Tapp, which performs a search with an a priori unknown number of solutions $`N_t`$, and finds a match (if it exists) in $`O(\sqrt{N/N_t})`$ steps. During the course of the algorithm the computer’s evolution must be tailored to accommodate the fact that the search is now based on all the target states that satisfy $`T[i]=n`$ where $`n`$ is some pre-defined Hamming distance determined by the algorithm. In order to apply the search algorithm in this case the operator $`I_SI_S(n)`$ is modified such that:

$$I_S(n)|\overline{\varphi _i}=\{\begin{array}{ccc}|\overline{\varphi _i}\hfill & & \mathrm{if}T[i]=n\hfill \\ & & \\ |\overline{\varphi _i}\hfill & & \mathrm{otherwise}.\hfill \end{array}$$
(21)
At each iteration the BBHT algorithm is employed, with a repeat index $`r`$ as a pre-determined measure of the search confidence level.

The optimal alignment algorithm is as follows:

1. 0th iteration: search for an occurrence of the state with zero Hamming distance, $`T[i]=0`$. If successful measure position and exit, if unsuccessful after $`r`$ repeats of the BBHT search algorithm go to the next iteration.
2. $`n`$th iteration: search for a state with $`T[i]=n`$ using $`U=I_HI_S(n)`$. If successful locate position and exit, if unsuccessful after $`r`$ repeats of the BBHT search algorithm go to the next iteration, by setting $`nn+1`$.
3. Upon exit at some iteration $`n=k`$, one optimal alignment $`T[i_k]=k`$, and its position $`i_k`$ has been found.

The total number of steps required is $`O(rk\sqrt{N})`$, discounting the effect of sequence repeats (which reduces the required number of iterations). At more cost a sub-loop may be introduced to search for the other optimally aligned sequences. In practice, the number of iterations required is $`k<<m`$, as one would determine a maximum Hamming distance on biological grounds, beyond which searching for an aligned state is pointless.

While the focus of this paper has been on protein sequence comparison, the framework can be easily translated into that for nucleotide sequence comparison in DNA. In this case representing the four letter nucleotide alphabet requires only two qubits.

Although only the algorithmic aspect of the application of quantum computing to sequence analysis has been dealt with here, an obvious point to raise is the feasibility of building such a device. With the ever increasing ability to manipulate systems at the quantum level there has been great progress in the demonstration of quantum computation at the two qubit level. Quantum logic gates were demonstrated using ion traps in 1995, and two years later in nuclear magnetic resonance (NMR) systems. In 1998 the actual experimental realization of a quantum computer solving Deutsch’s problem was reported by two groups using NMR. This was closely followed by NMR implementations of the quantum search algorithm. Of course, a realistic quantum computer needs to be scaled up significantly on these two qubit configurations. Perhaps the most promising prospect for a scalable quantum computer capable of running the algorithms presented here is based on the solid state design of Kane. The creation of a superposition representing the human genome database would be another considerable challenge.

To conclude, in this work the application of quantum search algorithms in the context of biocomputing has been studied, at least at the rather simple level of sequence alignment with respect to the Hamming distance. Actual alignment problems would include alignment through editing of sequences - i.e. insertion of gaps. It is quite possible that this procedure can be achieved using a multi-qubit representation (which includes gap characters) within the quantum search algorithm process by suitable choice of qubit evolution operators. Work in this direction is in progress.

The author wishes to thank S. O’Donoghue for many helpful discussions. Thanks also to G. Akemann and L. Benet Fernandez for kindly reading through the manuscript, and the theory group at MPIK for their hospitality. This work was supported by the Alexander von Humboldt Foundation.


text:
# The ages of quasar host galaxies.

## 1 Introduction

Determining the nature of the host galaxies of powerful active galactic nuclei (AGN) is of importance not only for improving our understanding of different manifestations of AGN activity, but also for exploring possible relationships between nuclear activity and the evolution of massive galaxies. The recent affirmation that black-hole mass appears approximately proportional to spheroid mass in nearby inactive glaxies (Maggorian et al. 1998) has further strengthened the motivation for exploring the link between active AGN and the dynamical and spectral properties of their hosts.

Over the last few years, improvements in imaging resolution offered by both space-based and ground-based optical/infrared telescopes has stimulated a great deal of research activity aimed at determining the basic structural parameters (i.e. luminosity, size and morphological type) of the hosts of radio-loud quasars, radio-quiet quasars, and lower-luminosity X-ray-selected and optically-selected AGN (e.g. Disney et al. 1995; Hutchings & Morris 1995; Bahcall et al. 1997; Hooper et al. 1997; McLure et al. 1999, Schade et al. 2000). However, relatively little corresponding effort has been invested in spectroscopic investigations of AGN hosts, despite the fact that this offers an independent way of classifying these galaxies, as well as a means of estimating the age of their stellar populations.

Our own work in this field to date has focussed on the investigation of the hosts of matched samples of radio-quiet quasars (RQQs), radio-loud quasars (RLQs) and radio galaxies (RGs) at relatively modest redshift ($`z=0.2`$). Details of these samples can be found in Dunlop et al. (1993). In brief, the sub-samples of quasars (i.e. RQQs and RLQs) have been selected to be indistinguishable in terms of their two-dimensional distribution on the $`Vz`$ plane, while the sub-samples of radio-loud AGN (i.e. RLQs and RGs) have been selected to be indistinguishable in terms of their two-dimensional distribution on the $`P_{5GHz}z`$ plane (as well as having indistinguishable spectral-index distributions). Deep infrared imaging of these samples (Dunlop et al. 1993; Taylor et al. 1996) has recently been complemented by deep WFPC2 HST optical imaging (McLure et al. 1999), the final results of which are reported by Dunlop et al. (2000). As well as demonstrating that, dynamically, the hosts of all three types of luminous AGN appear indistinguishable from normal ellipticals, this work has enabled us to deduce crude spectral information on the host galaxies in the form of optical-infrared colours. However, broad-baseline colour information can clearly be most powerfully exploited if combined with detailed optical spectroscopy. Over the past few years we have therefore attempted to complement our imaging studies with a programme of deep optical off-nuclear spectroscopy of this same sample of AGN.

Details of the observed samples, spectroscopic observations, and the basic properties of the observed off-nuclear spectra are given in a companion paper (Hughes et al. 2000). As discussed by Hughes et al. (2000), the key feature of this study (other than its size, depth, and sample control) is that we have endeavoured to obtain spectra from positions further off-nucleus ($`5`$ arcsec) than previous workers, in an attempt to better minimize the need for accurate subtraction of contaminating nuclear light. This approach was made possible by our deep infrared imaging data, which allowed us to select slit positions $`5`$ arcsec off-nucleus which still intercepted the brighter isophotes of the host galaxies (slit positions are shown, superimposed on the infrared images, in Hughes et al. 2000).

In this paper we present the results of attempting to fit the resulting off-nuclear spectra with evolutionary synthesis models of galaxy evolution. Our primary aim was to determine whether, in each host galaxy, the optical spectrum could be explained by the same model as the optical-infrared colour, and, if so, to derive an estimate of the age of the dynamically dominant stellar population. However, this study also offered the prospect of determining whether the hosts of different classes of AGN differ in terms of their more recent star-formation activity.

It is worth noting that our off-nuclear spectroscopy obviously does not enable us to say anything about the level of star-formation activity in the nucleus of a given host galaxy. Rather, any derived estimates of the level of star-formation activity refer to the region probed by our observations $`5`$ arcsec off-nucleus. However, given the large scalelengths of the hosts, their relatively modest redshift, and the fact that our spectra are derived from long-slit observations, it is reasonable to regard our conclusions as applying to fairly typical regions, still located well within the bulk of the host galaxies under investigation.

The layout of this paper is as follows. In section 2 we provide details of the adopted models, and how they have been fitted to the data. The results are presented in section 3, along with detailed notes on the fitting of individual spectra. The implications of the model fits are then discussed in section 4, focussing on a comparison of the typical derived host-galaxy ages in the three AGN subsamples. Finally our conclusions are summarized in section 5. The detailed model fits, along with corresponding chi-squared plots are presented in Appendix A.

## 2 Spectral fitting.

The stellar population synthesis models adopted for age-dating the AGN host galaxy stellar populations are the solar metallicity, instantaneous starburst models of Jimenez et al. (2000). We have endeavoured to fit each off-nuclear host galaxy optical spectrum by a combination of two single starburst components. For the first component, age is fitted as a free parameter, while for the second component the age is fixed at 0.1 Gyr, and the normalization relative to the first component is the only free parameter. This dual-component approach was adopted because single-age models are not able to adequately represent the data, and because it allows the age of the dynamically dominant component to be determined in a way which is not overly reliant on the level of ultraviolet flux which might be contributed either by a recent burst of star-formation, or by contamination of the slit by scattered light from the quasar nucleus. After experimentation it was found that the spectral shape of the blue light was better represented by the 0.1 Gyr-old (solar metallicity) model of Jimenez than by models of greater (intermediate) age (e.g. 1 Gyr). Moreover, the further addition of a third intermediate age component, $``$ 1 Gyr, did not significantly improve the quality of the fits achieved.

In general, our data are of insufficient quality to tell whether the component which dominates at $`\lambda 3000`$Å really is due to young stars, or is produced by direct or scattered quasar light. However, while the origin of the blue light is obviously of some interest, it has little impact on the main results presented in this paper, which refer to the age of the dynamically dominant stellar population which dominates the spectrum from longward of the 4000Å break through to the near infrared. The robustness of the age determination of the dominant population is demonstrated by comparing the results of the two stellar-population model fit with the results of a fit allowi ng for a nuclear contribution as well as the two stellar populations.

The model parameters determined were therefore the age of the dominant stellar population, which we can reasonably call the age of the galaxy, and the fraction, by (visible) baryonic mass, of the 0.1 Gyr component. For the fit including a nuclear contribution, the fraction of the total flux contributed by quasar light was a third parameter. The red end of each SED was further constrained by fitting $`RK`$ colour simultaneously with the optical spectral energy distribution. The fitting process is described below.

First, the observed off-nuclear spectra were corrected for redshift and transformed to the rest frame. They were then rebinned to the spectral resolution of the model spectra. The rebinned flux is then the mean flux per unit wavelength in the new bin and the statistical error on each new bin is the standard error in this mean. The data were normalised to a mean flux per unit wavelength of unity across the wavelength range 5020 $``$ 5500 $`\mathrm{\AA }`$.

The two-component model was built from the instantaneous-burst stellar-population synthesis model SEDs, so that

$$F_{\lambda ,age,\alpha }=const(\alpha f_{\lambda ,0.1}+(1\alpha )f_{\lambda ,age})$$

where $`f_{\lambda ,age}`$ is the mean flux per unit wavelength in the bin centred on wavelength $`\lambda `$ for a single burst model of $`age`$ Gyr, $`\alpha `$ is the fraction by mass of the young (0.1 Gyr) component, and $`F_{\lambda ,age,\alpha }`$ is then the new, mean twin stellar-population flux per unit wavelength in the bin centred on wavelength $`\lambda `$ for a model of $`age`$ Gyr. This composite spectrum was then normalised in the same way as the observed spectra.

A $`\chi ^2`$ fit was used to determine the age of the older stellar population and the mass-fraction, $`\alpha `$, of the younger population, for each host galaxy in the sample. The whole parameter space was searched, with the best-fit values quoted being those parameter values at the point on the grid with the minimum calculated $`\chi ^2`$. The normalization of the model spectra was allowed to float during the fitting process, to allow the best-fitting continuum shape to be determined in an unbiased way.

The models were fitted across the observed (rest-frame) spectral range, within the wavelength ranges listed in Table 1. As a result of the optimisation of the instruments with which the objects were observed, some of the galaxy spectra contain a ‘splice’ region where the red and blue end of the spectrum have been observed separately and then joined together (see Hughes et al. (2000) for details). These splice regions, defined in Table 1, were masked out of the fit in order to guard against the fitting procedure being dominated by data-points whose flux calibration was potentially less robust. The main emission lines, present due to nuclear light contamination or nebular emission from within the host, were also masked out, over the wavelength ranges given in Table 2.

The fit including a nuclear component was carried out in the same way, with the model flux in this case being

$$F_{\lambda ,age,\alpha ,\eta }=const(\alpha f_{\lambda ,0.1}+(1\alpha )f_{\lambda ,age}+\eta f0054_\lambda )$$

where $`\eta `$ is the fraction contributed to the total model flux by the nucleus, and $`f0054_\lambda `$ is the observed flux of the nucleus of the radio quiet quasar 0054+144. $`\alpha `$ is the fraction by mass of the total stellar population contributed by the 0.1 Gyr population. The wavelength range of the observed nuclear flux of 0054+144 is 3890$``$6950 $`\mathrm{\AA }`$. This was extended to 3500$``$8500 $`\mathrm{\AA }`$ by smooth extrapolation over the wavelength ranges 3500$``$3890 $`\mathrm{\AA }`$ and 6950$``$8500 $`\mathrm{\AA }`$ in order to carry out the $`\chi ^2`$ fit across the full wavelength range of the observed off-nuclear spectr a.

$`RK`$ colour was fitted in both cases with a typical error of a few percent. The observed $`RK`$ colours for the host galaxies are obtained from ukirt  and hst images (McLure et al. , 1999, Dunlop et al. , 2000), and define the basic shape of the host galaxy SED out to $`\lambda 2\mu m`$. The composite model spectra were appropriately red-shifted before calculating the colour, so that they could be compared to the observed colours without introducing uncertainties in k-correction. The $`R`$ band was simulated using the filter function, including system response and CCD quantum efficiency, for the hst wfpc2 F675W filter, and the K band was reproduced using the filter data for the ircam3 K Ocli filter at 77K combined with Mauna Kea atmosphere.

## 3 Results

The plots showing fits to individual spectra and $`\chi ^2`$as a function of fitted age are given in Appendix A. The plots for the two-component fit are presented in Figure A1, and those for the two stellar-component plus nuclear contribution are in Figure A2. The results for each object are summarized below, under their IAU names, with alternative names given in parentheses. Objects are listed in order of increasing right ascension, within each AGN class (radio loud quasars, radio quiet quasars and radio galaxies). The telescopes with which the spectra were obtained are also noted; M4M denotes the Mayall 4m Telescope at Kitt Peak, and WHT denotes the 4.2m William Herschel Telescope on La Palma. Where there are two spectra, the first plot is for the spectrum observed with the Mayall 4m Telescope, and the second is for the spectrum taken with the William Herschel Telescope.

### 3.1 Notes on individual objects

#### 3.1.1 Radio loud quasars

0137+012 (L1093) M4M
The models give a good fit at 13 Gyr, which is clearly improved by the inclusion of a small percentage (0.25%) of young stars. There is no significant nuclear contribution to the spectrum. HST imaging has shown that this host galaxy is a large elliptical, with a half-light radius $`r_e=13`$ kpc (McLure et al. 1999).
0736+017 (S0736+01, OI061) M4M,WHT
0736+017 has been observed with both telescopes, and fits to the two observed spectra are in good agreement. Both indicate an age of 12 Gyr. The M4M spectrum requires a somewhat larger young blue population (0.75%) than the WHT spectrum (0.125%). This may be due to poorer seeing at Kitt Peak leading to slightly more nuclear contamination of the slit, or to the use of slightly different slit positions at the two telescopes. However, this difference between the observed spectra shortward of 4000Å leaves the basic form of $`\chi ^2`$ versus age, and the best-fitting age of 12 Gyr unaffected. Inclusion of a nuclear component gives a much better fit to the blue end of the M4M spectrum, without changing the age estimation. The size of the fitted young populations are in much better agreement in this case. HST imaging has shown that morphologically this host galaxy is a large elliptical, with a half-light radius $`r_e=13`$ kpc (McLure et al. 1999).
1004+130 (S1004+13, OL107.7, 4C13.41) WHT
The spectrum of this luminous quasar certainly appears to display significant nuclear contamination below the 4000$`\mathrm{\AA }`$ break. As a result a relatively large young population is required to attempt (not completely successfully) to reproduce the blue end of the spectrum . However, the models predict that the underlying stellar population is old (12 Gyr). Allowing a nuclear component to be fit reproduces the blue end of the spectrum much more successfully, without changing the best-fit age estimation. HST $`R`$-band imaging indicates the morphology of the host galaxy is dominated by a large ($`r_e=8`$ kpc) spheroidal component, but subtraction of this best-fit model reveals two spiral-arm-type features on either side of the nucleus (McLure et al. 1999), which may be associated with the young stellar component required to explain the spectrum.
1020$``$103 (S1020$``$103, OL133) M4M
This object has the second bluest $`RK`$ of this sample, which leads to a much younger inferred age than the majority of the rest of the sample (5 Gyr), despite the presence of a rather clear 4000Å break in the optical spectrum. Ages greater than $``$10 Gyr are rejected by Jimenez’ models, primarily on the basis of $`RK`$ colour. HST imaging has shown that this host has an elliptical morphology, and a half-light radius of $`r_e=7`$ kpc. (Dunlop et al. 2000).
1217+023 (S1217+02, UM492) WHT
Nuclear contamination can again be seen bluewards of 4000$`\mathrm{\AA }`$, with a correspondingly large young population prediction for the purely stellar population model, which still fails to account for the very steep rise towards 3000Å. Hence, a large nuclear contribution is required to reproduce the blue end of the spectrum. The fit achieved by the models suggests that the dominant population is old, with a best-fit age of 12 Gyr. HST imaging has shown that this host has an elliptical morphology, and a half-light radius of $`r_e=11`$ kpc. (Dunlop et al. 2000).
2135+147 (S2135$``$14, PHL1657) WHT
2135+147 has a very noisy spectrum, but a constrained fit has still been achieved, and an old population is preferred. 2135+147 requires a large $`\alpha `$, even when a nuclear contribution is fitted. HST imaging has shown that this host has an elliptical morphology, and a half-light radius of $`r_e=12`$ kpc. (Dunlop et al. 2000).
2141+175 (OX169) WHT
This is another noisy spectrum, which has a relatively large quasar light contribution. An old population is again indicated by the model fits. ¿From optical and infrared imaging this object is known to be complex, but HST images indicate that it is dominated by a moderate sized ($`r_e=4`$ kpc) elliptical component (see McLure et al. (1999) for further details).
2247+140 (PKS2247+14, 4C14.82) M4M,WHT
2247+140 has been observed with both telescopes. The model fitting indicates an old population is required by both spectra - although the two observations do not agree precisely, the general level of agreement is very good, the two $`\chi ^2`$ plots have a very similar form, and the difference in $`\chi ^2`$ between the alternative best-fitting ages of 8 Gyr and 12 Gyr is very small. No significant nuclear contribution to the flux is present. HST imaging has shown that this host has an elliptical morphology, and a half-light radius of $`r_e=14`$ kpc (Dunlop et al. 2000).
2349$``$014 (PKS2349$``$01, PB5564) WHT
This is a very good fit to a good-quality spectrum, showing an obvious improvement when the low-level young population is added. Jimenez’ models clearly predict that the dominant population is old, with a well-constrained age of 12 Gyr. A very small nuclear contribution ($`\eta `$ = 0.050) does not significantly change the results. HST imaging of this object strongly suggests that it is involved in a major interaction, with a massive tidal tail extending to the north of the galaxy. However, the dominant morphological component is a spheroid with a half-light radius of $`r_e=18`$ kpc (McLure et al. 1999).

#### 3.1.2 Radio quiet quasars

0054+144 (PHL909) M4M,WHT
There is evidence of relatively large contamination from nuclear emission in the spectrum of this luminous quasar taken on both telescopes, and the age is not well-constrained, although the fit to the WHT spectrum derived from the models again suggests an old age. The $`\chi ^2`$ plots serve to emphasize how similar the two spectra of this object actually are (as also discussed by Hughes et al. 2000). Inclusion of a nuclear component in the model better constrains the age and improves the goodness of the fit. HST imaging of this object has shown that, morphologically, it is undoubtedly an elliptical galaxy, with a half-light radius $`r_e=8`$ kpc (McLure et al. 1999).
0157+001 (Mrk 1014) M4M,WHT
The age inferred from both the M4M and WHT spectrum of 0157+001 is again 12 Gyr. The apparently more nuclear-contaminated WHT spectrum does not give such a good fit, but 0157+001 is a complex object known to have extended regions of nebular emission, and the slit positions used for the two observations were not identical (Hughes et al. 2000). The age is much better constrained from the more passive M4M spectrum, to which Jimenez’ models provide a very good fit. Again, it seems that the nuclear contamination does not have a great influence on the predicted age of the old population, although the fit to the WHT spectrum is greatly improved by including a nuclear component in the model. Despite its apparent complexity in both ground-based and HST images, this host galaxy does again seem to be dominated by a large spheroidal component, of half-light radius $`r_e=8`$ kpc (McLure et al. 1999).
0204+292 (3C59) WHT
Jimenez’ models fit the spectrum of this object well, indicating an old underlying stellar population ($`>6`$Gyr), with a best-fit age of 13 Gyr for the stellar population plus nuclear component model, and a very small young population. HST imaging has shown this galaxy to be an elliptical, with half-light radius $`r_e=9`$ kpc (Dunlop et al. 2000).
0244+194 (MS 02448+19) WHT
The colour derived from the optical and infrared imaging of this host galaxy is rather blue ($`(RK)_{obs}`$ = 2.34) and this in part leads to a fairly young (5 Gyr) age prediction. However, the spectrum is very noisy, and, as indicated by the very flat $`\chi ^2`$ plot, the age is not strongly constrained. No nuclear flux contamination is fitted. HST imaging has shown this galaxy to have an elliptical morphology, with half-light radius $`r_e=9`$ kpc (McLure et al. 1999).
0923+201 WHT
The spectrum of 0923+201 is noisy, and it also appears to have some nuclear contamination. The fit is therefore improved by inclusion of a nuclear component. An old age is strongly preferred by the form of the $`\chi ^2`$ / age plot, with a best-fit value of 12 Gyr. HST imaging has shown this galaxy to have an elliptical morphology, with half-light radius $`r_e=8`$ kpc (McLure et al. 1999).
1549+203 (1E15498+203, LB906, MS 15498+20) WHT
This is a good fit, which is clearly improved by the addition of the younger population. The slope of the $`\chi ^2`$ / age plot strongly indicates an old dominant population, with a best-fit age of 12 Gyr. There is very little evidence of nuclear contamination. HST imaging has shown this galaxy to be a moderate-sized elliptical $`r_e=5`$ kpc (Dunlop et al. 2000).
1635+119 (MC1635+119, MC2) WHT
This is another very successful fit. An old age (12 Gyr) is inferred. Again, HST imaging has shown this galaxy to be a moderate-sized elliptical $`r_e=6`$ kpc (McLure et al. 1999).
2215$``$037 (MS 22152$``$03, EX2215$``$037) WHT
2215$``$037 has a noisy spectrum, to which an acceptable fit has nevertheless been possible. Jimenez’ models suggest an old population, with a best-fitting age of 14 Gyr, but this age is not well constrained. This is the only object where the inclusion of a nuclear contribution to the fitted spectrum substantially changes the age predicted of the dominant stellar population. However, the $`\chi ^2`$ plots are very flat after an age of 5 Gyr, and the age is not strongly constrained in either case. HST imaging has shown this galaxy to have an elliptical morphology, with $`r_e=7`$ kpc (Dunlop et al. 2000).
2344+184 (E2344+184) M4M,WHT
2344+184 has been observed with both telescopes, and the fits to both spectra are in good agreement $``$ although, formally, two different ages are predicted. This is because the $`\chi ^2`$ / age plots are fairly flat after about 8 Gyr. The fits to both observations suggest an old dominant population, with a small young blue population improving the fit. There is no significant change in the predictions when a nuclear flux component is included in the model. HST imaging has shown this to be one of the few host galaxies in the current sample to be disc-dominated. However, the nuclear component is in fact sufficiently weak that this object should really be classified as a Seyfert galaxy rather than an RQQ (McLure et al. 1999).

#### 3.1.3 Radio galaxies

0230$``$027 (PKS0230$``$027) WHT
0230$``$027 has a very noisy spectrum, and the colour of the host as derived from optical and infrared imaging is very blue ($`(RK)_{obs}`$ = 2.09). Consequently the best-fitting age derived using the models is 1 Gyr, with no younger component, but it is clear that little reliance can be placed on the accuracy of this result. Allowing for a contribution from quasar light does not improve the fit. HST imaging has shown this galaxy to have an elliptical morphology, with $`r_e=8`$ kpc (Dunlop et al. 2000).
0345+337 (3C93.1) WHT
0345+337 requires no young component or nuclear flux contribution at all, and an old age, of 12 Gyr, is clearly indicated by the models. HST imaging has shown this galaxy to be a large elliptical, with $`r_e=13`$ kpc (McLure et al. 1999).
0917+459 (3C219, 3C219.0) WHT
This is an excellent fit, with an old age produced by the models (12 Gyr), together with a very small young population. HST imaging has shown this galaxy to be a large elliptical, with $`r_e=11`$ kpc (McLure et al. 1999).
1215$``$033 WHT
Jimenez’ models suggest that the population of 1215$``$033 is universally old (best-fit age, 13 Gyr), with no young component or nuclear contribution required. HST imaging has shown this galaxy to be a large elliptical, with $`r_e=9`$ kpc (Dunlop et al. 2000).
1330+022 (3C287.1) M4M
An excellent fit to the data is produced by Jimenez’ models, indicating that the dominant population is old, with a best-fit age of 8 Gyr. HST imaging has shown this galaxy to be a large elliptical, with $`r_e=16`$ kpc (Dunlop et al. 2000).
2141+279 (3C436) M4M,WHT
This is another very successful, and well-constrained fit, with an inferred age of 12 Gyr. HST imaging has shown this galaxy to be a large elliptical, with $`r_e=21`$ kpc (McLure et al. 1999).

### 3.2 Sample overview

The results illustrated in Appendix A are summarised in Table 3. It should be noted that the 4000$`\mathrm{\AA }`$ break typical of evolved stellar populations is present in the majority of the observed spectra (see Appendix A and Hughes et al. 2000), so we can be confident in fitting stellar population models to the data. The plots clearly show that the addition of even a very small amount of secondary star formation to the simple, near-instantaneous star-burst models reproduces the blue end of the observed host galaxy spectra much more successfully (and in most cases very well) than does a single stellar population. Including a nuclear component further improves the fit to the blue end, especially for those spectra not well fit by purely stellar light. At the same time, the red end of the spectra, plus the observed $`RK`$ colours generally require that the underlying stellar populations are old.

The most meaningful output from the model fitting is a constraint on the minimum age of the host galaxies; the $`\chi ^2`$ plots show that often the best-fit age is not strongly constrained but the trend is clearly towards old ($``$ 8 Gyr) stellar populations. In general, young populations are strongly excluded.

The peaks and troughs in $`\chi ^2`$ as a function of model age appear to be the result of real features of the population evolution synthesis, rather than being due to, for example, poor sampling (Jimenez, private communication).

## 4 Discussion

Figure 1 shows the distribution of best-fit ages estimated using Jimenez’ models. The left-hand panel shows the ages of the dominant stellar component which result from fitting the two-component (stars only) model, while the right-hand panel shows the ages of the dominant stellar component in the case of the 3 component model (i.e. two stellar components + a contribution from scattered nuclear light). The host galaxies of the AGN in each sample are predominantly old, and this result is unaffected by whether or not one chooses to include some nuclear contribution to assist in fitting the blue end of the spectrum. This general result also appears to be relatively unaffected by the precise choice of stellar population model; the models of Bruzual and Charlot (2000), fitted to the host galaxies with the same process, reproduce the results for Jimenez’ models to within a typical accuracy of 1-3 Gyr. We have not attempted to fit the stellar population models of Yi et al. (2000), because of problems we have found in their MS rate of evolution (Nolan et al., 2000).

Inclusion of a nuclear contribution in the 3-component models obviously raises the question of whether a young stellar population component is really necessary at all. We thus also explored the results of fitting a two-component model consisting of a single stellar population plus nuclear contribution. Such a model adequately reproduces the two-population-plus-nuclear-component results for the redder galaxies, but in fact the spectra of the bluest galaxies cannot be adequately reproduced by these models; the resulting serious increase in minimum $`\chi ^2`$ demonstrates that inclusion of a young population component is necessary to achieve a statistically acceptable fit to these data.

Figure 2 shows the distribution of percentage contribution (by mass), $`\alpha `$, of the young (0.1 Gyr) component to the spectra of the hosts. The left-hand panel shows the values of $`\alpha `$ as derived from fitting the two-component (stars only) model, while the right-hand panel shows the values of $`\alpha `$ produced by the 3 component model (i.e. two stellar components + a contribution from scattered nuclear light).

Where two spectra of the same object have been obtained with alternative telescopes/instruments, the derived ages of the dominant stellar componenets are reproduced reassuringly well. There are, however, small discrepancies in the estimated percentage of young population present (see Table 3). This effect is suggestive that at least some of the blue light might be due to a scattered nuclear contribution, the strength of which would be highly dependent on the seeing at the time of observation, and on the precise repeatability of slit placement relative to the galaxy core. Interestingly, when a nuclear component is included with the stellar flux model, an even smaller percentage of 0.1 Gyr stellar population is required to fit the blue end of the spectra, and (more importantly) the difference in $`\alpha `$ between two spectra of the same object is generally reduced. This provides further support for the suggestion that some of the bluest quasar host spectra remain contaminated by quasar light at the shortest wavelengths, and indicates that the right-hand panel of Figure 2 provides a more realistic estimate of the level of on-going star-formation in the host galaxies. While this figure still appears to suggest that at least some quasar hosts display higher levels of ongoing star-fomration activity than do radio galaxies, statistically this ‘result’ is not significant.

There are three galaxies which have very low age estimates, namely 0230$``$027, 0244+1944 and 1020$``$103. These objects have the bluest observed $`RK`$ colours, so it may be expected that the fitted ages would be younger than the rest of the sample, and that these populations are genuinely young. As discussed above, it may be that these objects are bluer because of scattered nuclear light contaminating the host galaxy spectrum. However, it seems unlikely that their fitted ages are low simply because of this, because elsewhere in our sample, where two spectra have been obtained of the same object, the amount of nuclear contamination present does not significantly affect the age estimation (e.g. 0736+017 and 0157+001). Moreover, the inclusion of a nuclear component to the fit does not change the estimated age distribution of the host galaxies. If these ages are in error, then a more likely explanation, supported by the relative compactness of these particular host galaxies, is that nuclear and host contributions have been imperfectly separated in the $`K`$-band images, leading to an under-estimate of the near-infrared luminosity of the host.

The result of the three-component fitting process which also allows a contribution from scattered nuclear light is that there are in fact only 3 host galaxies in the sample for which there is evidence that $`\alpha >0.5`$. One of these is the host of a radio-loud quasar (2135$``$147) but, as can be seen from Fig A2, this spectrum is one of the poorest (along with 0230$``$027, the only apparently young radio galaxy) in the entire dataset. The 2 convincing cases are both the hosts of radio-quiet quasars, namely 0157$`+`$001 and 2344$`+`$184.

Within the somewhat larger sample of 13 RQQs imaged with the HST by McLure et al. (1999) and Dunlop et al. (2000), 4 objects showed evidence for a disk component in addition to a bulge, namely 0052+251, 0157+001, 0257+024, and 2344+184. Since we do not possess spectra of 0052+251 and 0257+024 this means that there is a 1:1 correspondence between the objects which we have identified on the basis of this spectroscopic study as having recent star-formation activity, and those which would be highlighted on the basis of HST imaging as possessing a significant disk component. This straightforward match clearly provides us with considerable confidence that the spectral decomposition attempted here has been effective and robust. Finally we note that it is almost certainly significant that 0157+001, which has the largest starburst component ($`\alpha =1.1`$) based on this spectroscopic analysis, is also the only IRAS source in the sample.

## 5 Conclusions

We conclude that the hosts of all three major classes of AGN contain predominantly old stellar populations ($`11`$ Gyr) by $`z0.2`$. This agrees well with the results of McLure et al. (1999), and Dunlop et al. (2000) who compare host galaxy morphologies, luminosities, scale lengths and colours in the same sample, and conclude that the hosts are, to first order, indistinguishable from ‘normal’ quiescent giant elliptical galaxies.

The best-fitting age of the dominant stellar population is not a function of AGN class. For the purely stellar models, the fitted percentage contribution of the blue component is, however, greater in the quasar hosts than in the radio galaxies; the median values are 0.6% for the 9 radio-loud quasars, 0.6% for the 9 radio-quiet quasars, and 0.05% for the 6 radio galaxies. However, when a nuclear component is included, the median values are 0.3% for the radio-loud quasars, 0.3% for the radio-quiet quasars, and 0.00% for the radio galaxies. Performing a Kolmogorov-Smirnov test on these results yields a probability greater than 0.2 that the percentage of young stellar population in host galaxies is in fact also not a function of AGN class.

These results strongly support the conclusion that the host galaxies of all three major classes of AGN are massive ellipticals, dominated by old stellar populations.
ACKNOWLEDGEMENTS
Louisa Nolan acknowledges the support provided by the award of a PPARC Studentship. Marek Kukula and David Hughes acknowledge the support provided by the award PPARC PDRAs, while Raul Jimenez acknowledges the award of a PPARC Advanced Fellowship. We thank an anonymous referee for perceptive comments which helped to clarify the robustness of our results and improved the clarity of the paper.
REFERENCES
Bahcall J.N., Kirhakos S., Saxe D.H., Schneider D.P., 1997, ApJ, 479, 642
Bruzual A.G., Charlot S., 2000, in preparation
Disney M.J. et al., 1995, Nat, 376, 150
Dunlop J.S., 2000, In: ‘The Hy-Redshift Universe: Galaxy Formation and Evolution at High Redshift’, ASP Conf. Ser., Vol 193, eds. A.J. Bunker & W.J.M. van Breugel, in press, (astro-ph/9912380)
Dunlop J.S., McLure R.J., Kukula M.J., Baum S.A., O’Dea C.P., Hughes D.H., 2000, MNRAS, in press
Dunlop J.S., Taylor G.L., Hughes D.H., Robson E.I., 1993, MNRAS, 264, 455
Hooper E.J., Impey C.D., Foltz C.B., 1997, ApJ, 480, L95
Hughes D. H., Kukula M.J., Dunlop J.S., Boroson T., 2000, MNRAS, in press
Hutchings J.B., Morris S.C., 1995, AJ, 109, 1541
Jimenez R., Dunlop J.S., Peacock J.A., Padoan P., MacDonald J., J$`ø`$rgensen U.G., 2000, MNRAS, in press
Jimenez R., Padoan P., Matteucci F., Heavens A.F., 1998, MNRAS, 299, 123
Magorrian J., et al., 1998, AJ, 115, 2285
McLure R.J., Kukula M.J., Dunlop J.S., Baum S.A., O’Dea C.P., Hughes D.H., 1999, MNRAS, in press, (astro-ph/9809030)
Nolan L.A., Dunlop J.S., Jimenez R., 2000, MNRAS, submitted, (astro-ph/0004325) Schade, D.J., Boyle, B.J., Letawsky, M., 2000, MNRAS, in press
Taylor G.T., Dunlop J.S., Hughes D.H., Robson E.I., 1996, MNRAS, 283, 930
Yi S., Brown T.M., Heap S., Hubeny I., Landsman W., Lanz T., Sweigart A., 2000, ApJ, submitted

## Appendix A Spectra and $`\chi ^2`$ plots

The fits for all the off-nuclear spectra are given. In Fig A1, the rest frame spectra are in the first column (black), with the best-fitting two-component model spectra (Jimenez et al., 2000) superimposed (green). The spectra of the single-aged old popula tion (red) is given for comparison. In Fig A2, the fits allowing for an additional contribution to the flux from the nucleus are presented. The key is the same as in A1, with the additional blue line representing the best-fitting two-component model flux plus the nuclear flux contribution.

The second column of plots shows the $`\chi ^2`$evolution with age for the dominant older population. The third column shows the best-fit $`\chi ^2`$as a function of percentage young population, $`\alpha `$, for fixed ages of the dominant component. The models have solar metallicity. The subscript $`\eta `$ denotes results obtained by including the nuclear contribution.

Where there are two spectra of the same object, the spectrum given first is the one observed on the Mayall 4m Telescope, and the second is that observed using the William Herschel Telescope.The data for the following objects have been smoothed using a Han ning function: 2135+147 (RLQ), 2141+175 (RLQ), 0244+194 (RQQ), 0923+201 (RQQ), 1549+203 (RQQ), 2215-037 (RQQ), 0230-027 (RG) and 0345+337 (RG).