text:
# Galactic Structure and Evolution: a decade of surveys

## 1. Introduction

The mass assembly history of galaxies remains one of the critical issues in observational cosmology: did galaxies reach their present stellar mass only recently (say, at $`z1`$)? Or were most (massive) galaxies already in place by $`z1`$? Is late/current accretion significant or merely a perturbation? These questions may best be addressed by a combination of detailed surveys of the Local Group, where quantitative resolved information may be obtained, together with spectroscopic surveys of galaxies as a function of redshift. Substantial progress on both fronts is being made through the wealth of surveys underway and recently reported.

Surveys are indeed the currently favoured means to progress in astrophysics: a search through the journals or astro-ph for ‘survey’ immediately overloads any reading list. Equally noticeable is the current fashion to add a speculation about current mergers to any paper studying any object in the Local Group, irrespective of its direct relevance to the specific data at hand. More fundamental is to consider if the rate and amplitude of the (clearly still happening) current accretion into the Milky Way is signal or noise, and is or is not consistent with model expectations.

That is, we should not aspire to show that current formation models correctly predict the Local Group, and its current merger activity: we know, from direct observation, that current models require many better constraints to allow their improvement, to correspond to real galaxies. Our goal is to consider the observations objectively, to quantify those constraints.

Spectroscopic surveys of faint galaxies, especially those selected in the $`K`$-band, currently offer an excellent opportunity to address these questions for the distant universe (Broadhurst et al. 1992; Cimatti etal 2003: see especially http://www.arcetri.astro.it/$``$k20/). The results of these surveys, together with current simulations of galaxy formation, provide the context in which Local Group studies can be interpreted, by highlighting those aspects of current models where most progress is required.

One specific result of current deep surveys, which has been long established in the Galaxy, is the apparently long time since the galaxies were assembled. This result, long known locally, and the subject of many continuing local investigations (eg Sandage, Lubin and VandenBerg 2003; see also Nördstrom etal 2004) identifies one of the most promising ways in which our understanding of galaxy formation can be improved and extended. In essence, galaxies seem to have formed earlier than current models expect, and have remained (relatively) undisturbed for a longer time than predicted. Extending and explaining this contradiction promises to improve significantly our understanding not only of galaxy evolution, but also of the coupled question, the distribution of Dark Matter on small scales, and its corresponding temperature.

In one of the most impressive of recent studies, Cimatti etal (2002) present the redshift distribution of a complete sample of 480 galaxies with $`K_s<20`$ distributed over two independent fields covering a total area of 52 arcmin2. A “blind” comparison was made with the predictions of a set of the most recent $`\mathrm{\Lambda }`$CDM hierarchical merging (HMM) and a set of pure luminosity evolution (PLE) models. The hierarchical merging models overpredict and underpredict the number of galaxies at low-$`z`$ and high-$`z`$ respectively, whereas the PLE models match the median redshift and the low-$`z`$ distribution, while still being able to follow the high-$`z`$ tail of $`N(z)`$. These results are summarised in Figure 1.

These observational surveys illustrate the areas where current paradigmatic models can be most improved. In particular, $`\mathrm{\Lambda }`$CDM hierarchical merging models as they are currently implemented overpredict the total number of galaxies with $`K_s<20`$ in the K20 survey area by factors of 30-45%, and are inconsistent with the observed galaxy redshift distribution by substantial factors, especially at high redshifts. A Kolmogorov-Smirnov test showed that all the hierarchical merging models disagree with the observations at $`>99`$% level.

Nonetheless, the hierarchical merging $`\mathrm{\Lambda }`$CDM scenario has spectacular success in reproducing large scale structure and CMB results, and is certainly the paradigm of choice in which we should analyse survey and galactic structure data today. Our ambition is to improve the ability of this framework to become consistent with the real physics of small-scale structure and galaxy formation, by identifying areas in which the input physics can be improved. The observed discrepancies derived from the high-redshift studies are related to what is essentially a purely heuristic set of algorithms adopted to describe the star formation processes and their feedback, both within individual galaxies and in their environment. The high-redshift results suggest that the galaxy formation recipes should allow galaxy formation in a CDM dominated universe to occur in a manner described by the old-fashioned pre-CDM monolithic collapse scenario. This requires the recipes to allow enhanced merging and star formation in massive haloes at high redshift (say, $`z3`$), while at the same time suppressing star formation in low-mass haloes.

We now ask if this improved recipe would be consistent with the Milky Way, an L galaxy whose properties should not be too rare in any plausible formation scenario.

## 2.  Mergers, phase-space structures, fashion accessories

Search for and interpretation of phase-space structure in the Galaxy has always been important, but has recently become a major industry. Significant structure is unexpected in a context where galaxy photometrists provid convincing evidence for symmetric and smooth surface brightness profiles (eg Kormendy & Djorgovski 1989), where density profiles for special tracers (RR Lyraes being popular) show smooth profiles over very many scale lengths, and where local stellar kinematics looks closely Gaussian (cf Gilmore, Wyse & Kuijken 1989): clearly whatever process makes galaxies makes smooth-looking relaxed-looking galaxies. Nonetheless, testing the paradigm is important. The current fashion is perhaps to start from the other extreme (cf Freeman & Bland-hawthorn 2002) where some aspire to “associate components of the Galaxy to elements of the protocloud”, an interesting alternative to the statistical distribution function analyses on which the science of stellar dynamics has been built.

Earlier searches for structure in the stellar distribution used techniques such as correlation functions (eg Gilmore, G.; Reid, N.; Hewett, P. 1985) or presumed local evolutionary processes were the default explanation of identified field halo clustering (eg Arnold & Gilmore 1992). The many studies of Eggen on moving groups were perhaps not given due recognition, since most were of thin disk stars, disk structure is inevitably lumpy in phase space, and tells little of global or long-term galaxy evolution. His study of the Arcturus Moving Group was however long understood as significant for thick-disk studies (cf Fuhrmann 2004, sect 2.4) and appears as such in encyclopeadias, though is still being rediscovered. \[Eggen however (private communication) was deeply sceptical of the distinction between thin and thick disk moving groups.\] gould (2003) has recently reminded us that an apparently smooth distribution is either really smooth - in which case its history is lost into statistical distribution functions – or the central limit theorem applies – in which case its history is lost into statistical distribution functions.

The field changed rapidly after discovery of the Sgr dwarf galaxy (Ibata, Gilmore & Irwin 1994, 1995), an object genuinely discovered in phase-space in real-time at the AAT where Ibata and Gilmore were plotting colour-magnitude-velocity diagrams during their study of the Galactic Bulge (Ibata & Gilmore 1995a, 1995b). While this discovery certainly activated the field (though the forthcoming 2MASS and SDSS surveys would have made the same discovery through traditional photometric, rather than phase-space, techniques) the reson it was possible seems not to have been fully appreciated by those attempting to emulate the discovery: Sgr was found because its stellar population is very different from that of the Galactic bulge or the Galactic halo. It is an intermediate age and intermediate metallicity population, and is unique in that respect.

One (should have) immediately deduced that Sgr was a rare event, not a paradigm for the average (see Unavane, Wyse & Gilmore 1996). The outcome of a truly impressive and truly vast amount of work since its discovery mapping the distribution of Sgr-related debris trails around the Galaxy has reinforced this early conclusion: while much of Sgr is strewn around the Galaxy, nothing else like it has left any trace of its demise in the last Hubble time. Except, of course, the Thick Disk (Gilmore & Reid 1983), whenever that got into place.

## 3. The chemical elements as archaeological tracers

The recent surveys which have delivered the most quantitative information on Galactic history are those which provide high-precision elemental abundances for many stars. There are many such recent surveys since the arrival of 8m telescopes with excellent echelle spectrographs, with many more due to be published soon. These range from studies of the most extreme metal-poor stars, through the ‘halo’, ‘thick disk’ and ‘thin disk’ metallicity ranges, with only the Bulge and inner disk still awaiting useful statistics. The bulge and inner disk studies are especially needed, since the common (modellers’) presumption that ‘the average star forms in a disk, and then merges into a spheroid’ merits a check.

The results from these studies are outstanding science, and provide information which must be considerde fully in any discussion of Galaxy formation and evolution.

The common feature to all studies is the extremely well-defined pattern of element ratios, spanning huge dynamic ranges. One example is shown in Figure 2, taken from a superb study by Gratton etal 2003. This figure shows the common result: different stellar population show smooth systematic element ratio patterns across very wide dynamic ranges.

The systematic chemical abundance patterns evident here, and the different systematics between different stellar populations, provide undisputable evidence that the formation of the field stars which are today in the Galactic halo happened in a single environment, or a range of environments, in which the rate of star formation, the rate of mass loss and the rate of chemical enrichment had an extremely small range. Noone has yet suggested how this is possible if indeed the halo is built up over time from a large number of systems with uncorrelated histories.

Other recent studies extend this conclusion to the stars of the thick disk and the old thin disk.

## 4. Merger rates

Mergers happen. The Hubble-Toomre sequence is built that way.

The Sagittarius dwarf is direct evidence. In M32 the debris from the stripped elliptical M32 has recently been observed directly, decades after its prediction. The Magellanic Stream exists. It is possible that a further debris remnant has been discovered at low latitudes, although at the time of writing the ‘normal’ complex structure of an outer warped disk is a more likely explanation of the observations. Similarly, speculations that structures, such as clusters, could be trapped into the Galactic disk from a !merger parent’ need to treated with some caution. The depth of a potential scales as the square of its velocity dispersion, Thus the Galactic halo potential (dispersion $``$100km/s) is an order of magnitude deeper than that of the disk (dispersion $``$30km/s). Trapping something into the outer disk is not an obvious dynamical process without very special geometries, and adroit use of dynamical friction.

Is Sgr all there is?

Phase space structures are not an indication of mergers. The most dramatic local phase-space structures are spiral arms, Gould’s belt, and star-formation regions. Even in the halo, very considerable phase-space structure must exist even if there has never been a single merger in the Hubble time.

Globular clusters evolve, and evaporate eventually. It is often speculated that many hundreds of dissolved globular clusters exist in the halo. A globular cluster, Pal5, has recently been observed in the final stages of dissolution (Odenkirchen etal 2003) and is shown in Figure 3.

The same limits on accretion histories which can be deduced from the difference of the stellar populations of Sgr, and are quantified for field stars by Unavane, Wyse and Gilmore (1996), can be applied to other tracers. The severe difficulty of hiding stars with different enrichment histories in the element ratio data is noted above. One may do the same thing for globular clusters. Bellazzini and collaborators are extending the census of former members of Sgr utilising that fact that Sgr stellar populations are different from field stars. One may do the same slightly more generally.

Figure 4 shows the size distribution (core radii, from of the different groups of Galactic globular clusters, and those old clsuters in Galactic satellites. Most clusters in Galactic satellites are young, thereby trivially excluding any (non-Sgr) mergers in the last $``$6-8Gyr with anything resembling anything which exists in the Local group now, or its precursors. Even older mergers are interesting in this context: the only Galactic globulars which are consistent, on size grounds, with having arrived in a merger are the few younger ones, including those associated already with Sgr.

Once again, we conclude that mergers with sufficient mass to perturb the history of the Galaxy in an observable way must be rare in recent times.

## 5. Field star kinematic surveys: much still to learn

Kinematic surveys of field stars in the Galaxy remain as rare as is evidence for multiple mergers in the Milky Way. Each of the few attempted has made significant discoveries (eg SDSS) which should be an encouragement. One for which the first results have recently been reported indicates the discovery potential. Gilmore, Wyse and Norris (2002) report first results from a survey using 2dF on the AAT. their first results are summarised in Figure 5, and show significant disagreement with expectation.

The mean star a few kpc from the Plane has angular momentum intermediate between that of the halo and the thin disk. This argues strongly for an independent origin, and is consistent with - but does not require - a merger origin.

## 6. Conclusion

Mergers happen. The Hubble-Toomre sequence is built that way (Figure 6).

This is an age of surveys. We are learning to process and distribute vast data volumes. We are learning to extract astrophysics. While fashions still reign in interpretation, the story in the data is there to be read: it tells the story of Galaxy formation.

## References

Arnold, R,. Gilmore, G 1992 MNRAS 257 225

Cimatti etal 2002 A+A 391 L1 2002

Cimatti etal ESO Messenger, March 2003, No. 111, p. 29.

Freeman, K., Bland-Hawthorn, J., 2002 ARAA 40 487

Fuhrmann, K. 2004 Astr Nachr 325 1

Gilmore,G., Reid, I.N 1983 MNRAS 202 1025

Gilmore, G.; Reid, N.; Hewett, P. 1985 MNRAS 213 257

Gilmore, G., Wyse, R., & Norris, J. 2002 ApJ 574 L39

Gilmore, G., Wyse, R., & Kuijken, K. 1989 ARAA 27 555

Gould, A. 2003 ApJ 529 L63

Gratton etal 2003 2003A&A 406 131

Ibata, R., Gilmore G 1995a MNRAS 275 591

Ibata, R., Gilmore G 1995b MNRAS 275 601

Ibata, R., Gilmore, G., Irwin, M. 1994 Nature 370 194

Ibata, R., Gilmore, G., & Irwin, M. 1995 MNRAS 277 781

Kormendy, J., & Djorgovski, G., 1989 ARAA 27 235

Mackey, D., & Gilmore, G. 2003a MNRAS 338 85

Mackey, D., & Gilmore, G. 2003b MNRAS 338 120

Mackey, D., & Gilmore, G. 2003c MNRAS 340 175

Nordstrom, B etal 2004 ‘ the Geneva-Copenhagen survey of the Solar neighbourhood, A+A 2004 in press

Odenkirchen. M., etal 2003 Astron.J. 126 2385

Sandage, A., Lubin, L., VandenBerg, D., PASP 115 1187 2003

Unavane, M., Wyse, R., & Gilmore, G 1996 MNRAS 278 727


format: text

text:
# Angular Momentum Conservative Algorithm of Collisional Process in DSMC Method

## 1 Introduction

The traditional algorithm of collisional process in DSMC method, used for calculation of velocities of molecules after collision, is based on the conservation laws for linear momentum and total energy . This algorithm does not provide conservation of angular momentum with respect to some axis, but for non-rotating flow this circumstance is insignificant, since the mean value of angular momentum of such kind of flow is equal to zero.

The situation changes for axisymmetric flow with rotation. In this case the traditional procedure leads to some internal source or sink of angular momentum that may distort the flowfield.

In the paper the discussed effect is recognized and investigated. The angular momentum conservative (AMC) algorithm of collisional process is proposed. The draft version of the paper is published in .

## 2 Problem Formulation

The investigations were made for one-dimensional axisymmetric flow of monoatomic gas in the tube with specular wall. The evolution of initially swirling flow during the time was studied. The gas in the tube with uniform initial density $`n_0`$ and temperature $`T_0`$ and equilibrium distribution function was assumed to start rotating as a solid body, the initial tangential velocity on the wall being equal to $`v_\phi `$. For given molecular model the problem contains two governing parameters: the Knudsen number $`Kn`$, defined by the ratio of the mean free path $`l_0`$ in the gas at $`t=0`$ to the radius of the tube $`r_t`$

$$Kn=\frac{l_0}{r_t}$$

and the speed ratio

$$S_\phi =v_\phi /\sqrt{2RT_0}$$

The simulations were made on uniform grid with $`10^4`$ simulated molecules. Standard NTC DSMC procedure  was applied.

The initial values of total mass $`M_0`$, angular momentum $`I_0`$ and energy $`E_0`$ of the gas may be defined by the relations

$$M_0=mN_0=mn_0\pi r_t^2$$
(1)
$$I_0=m\underset{i=1}{\overset{N_0}{}}v_ir_i$$
(2)
$$E_0=m\underset{i=1}{\overset{N_0}{}}\frac{𝑽_i^2}{2}$$
(3)
where $`N_0`$ is the total number of simulated molecules, $`m`$ \- the mass of molecule, $`𝑽`$ \- the velocity vector, $`v`$ \- the tangential component of this vector. By the same relations may be determined the total energy $`E`$ and total angular momentum $`I`$ of the gas during the time.

For any values of $`Kn`$ and $`S_\phi `$ there should be no temporal dependence of $`E`$ and $`I`$. The computations show that this condition is satisfied only for $`E`$.

Fig.1 illustrates the dependencies of $`I/I_0`$ on the time for some of considered variants for $`S_\phi =\mathrm{\hspace{0.17em}1}`$. The time is normalized by the value of $`t_0=\mathrm{\hspace{0.17em}2}\pi r_t/v_\phi `$, so $`t/t_0`$ represents the number of revolutions of solid body with $`v(r_t)=v_\phi `$.

As can be seen from Fig.1, for all of the variants the decreasing of $`I/I_0`$ with increasing of $`t/t_0`$ is observed. For given $`Kn`$ the effect is more pronounced for solid sphere molecules $`(SS)`$ in comparison with $`VSS`$ molecular model for Maxwell molecules $`(MM)`$. The effect depends on the Knudsen number: the less $`Kn`$ the more the effect. But the highest sensitivity of the effect is observed to the grid size: changing the number $`N`$ of cells from $`50`$ to $`25`$ causes the same effect as the decrease of the Knudsen number approximately by $`4`$ times.

Due to the dependence of $`I(t)`$ it is difficult to obtain the steady solution of the problem. The dissipation of rotational contribution of energy into the heat takes place during the time. The only steady solution that may be obtained with traditional collisional procedure is terminal one, when at $`t\mathrm{}`$ the gas will stop rotating, having uniform density $`n_0`$ and getting warm to the temperature

$$T=\mathrm{\hspace{0.17em}2}E_0/\mathrm{\hspace{0.17em}3}RM_0$$

## 3 AMC Algorithm

The reason of temporal dependence of total angular momentum is quite clear: the colliding molecules have different radial coordinates $`r_1,r_2`$, that leads to the difference in angular momentum

$$m_1v_1r_1+m_2v_2r_2$$

of these molecules before and after collision ($`v`$ denotes the tangential component of velocity of the molecule). To provide the conservation of angular momentum the algorithm of collisional process should be changed. One possible variant of AMC algorithm is as follows.

Let us denote by $`𝒖`$ the vector containing axial and radial components of the velocity of molecule. The post-collision velocities calculated by traditional procedure will be denoted by symbol $``$, while the symbol  will be used for these velocities in the proposed algorithm.

The algorithm is based on conservation law of angular momentum $`mvr`$ for tangential component of velocity instead of linear momentum $`mv`$. That is why the post-collision tangential velocities are calculated by usual relations , with the value of $`m_i`$ being replaced by the product $`m_ir_i`$. The obtained velocities $`v_1^{},v_2^{}`$ provide precise conservation of angular momentum, but violate the energy conservation. To provide energy conservation the correction of velocities $`𝒖_1^{}`$ and $`𝒖_2^{}`$ is needed. To do this correction the analysis of energy defects

$$dE_i=m_i(v_i^2v_i^2)$$

should be made. Three possible relations between energy defects $`dE_i`$ and velocities $`𝒖_i^{}`$ are as follows.

1. If

$$m_1u_1^2>dE_1$$

and

$$m_2u_2^2>dE_2$$

both components of vectors $`𝒖_i^{}`$ should be corrected by the factors $`c_i`$

$$c_i=\sqrt{1dE_i/m_iu_i^2}$$

2. If

$$m_1u_1^2+m_2u_2^2>dE_1+dE_2$$

the correction factor $`c`$ for components of both velocities $`𝒖_i^{}`$ should be the same and equal to

$$c=\sqrt{1(dE_1+dE_2)/(m_1u_1^2+m_2u_2^2)}$$

3. If none of the above conditions is satisfied, the collision is considered to be ”bad” and is not performed, i. e. the molecules conserve their pre-collision velocities.

The computations show that the relative number of collisions of types 1 and 2 is about $`99\%`$ and $`1\%`$, respectively. The relative number of ”bad” collisions never exceeds $`10^4`$, so the effect caused by the neglecting of these collisions is small.

## 4 Steady Solution

The described algorithm provides precise conservation of angular momentum and energy. All the curves, presented in Fig.1 transform into straight line $`I/I_0=\mathrm{\hspace{0.17em}1}`$, if this algorithm is applied.

AMC algorithm enables to obtain the steady solution of the considered problem. This solution proves to be the same, as the prediction of Navier-Stokes theory and has the following features.

1. The flow is isothermal

$$T=const$$
(4)
2. The gas rotates as a solid body

$$v=v_wr/r_t$$
(5)
3. The radial distribution of density is described by the relation

$$n(r)=n(0)\mathrm{exp}(v^2/2RT)$$
(6)
4. The solution is completely determined by the values of $`M_0`$, $`I_0`$, $`E_0`$ and does not depend on the way of initial swirling.

5. The solution does not depend either on the molecular model or on the Knudsen number.

These features of the considered flow are important for clarifying the nature of the Ranque effect .

Three unknown values $`v_w,T,n(0)`$ determining the radial distribution of parameters of steady solution may be defined based on the values of $`M_0`$, $`I_0`$ and $`E_0`$ from the relations

$$\pi r_t^2mn(0)\frac{(\mathrm{exp}S_w^21)}{S_w^2}=M_0$$
(7)
$$M_0v_wr_t\frac{\mathrm{exp}S_w^2(S_w^21)+1}{S_w^2(\mathrm{exp}S_w^21)}=I_0$$
(8)
$$\frac{3}{2}RTM_0+\frac{v_wI_0}{2r_t}=E_0$$
(9)
where speed ratio $`S_w`$ of the flow at $`r=r_t`$ is defined by

$$S_w=v_w/\sqrt{2RT}$$
(10)
The left hand sides of the relations 7 \- 9 representing the total mass, angular momentum and energy of the flow may be obtained based on the relations 4 \- 6 by simple integration.

The comparison between numerical and analytical radial distribution of parameters of steady flow for $`VSS`$ molecules, $`Kn=\mathrm{\hspace{0.17em}0.1}`$, $`S_\phi =\mathrm{\hspace{0.17em}1}`$, $`N=\mathrm{\hspace{0.17em}50}`$ is made on Figs. 2 \- 4 for density, tangential velocity and temperature, respectively. The numerical results are shown by solid circles while solid lines represent the results of analytical solution 7 \- 9. As can be seen from these Figures, numerical and analytical results are in excellent agreement. For $`Kn=\mathrm{\hspace{0.17em}10}`$ the steady solution was found to be the same as for $`Kn=\mathrm{\hspace{0.17em}0.1}`$.

It should be noted that the steady flow in the considered problem is characterized by locally Maxwellian distribution function for any molecular models and $`Kn<\mathrm{}`$.

AMC algorithm of collisional process may be applied for simulation of non-rotating flows too. The experience of such application to the flow in gasdynamic window  shows that the results obtained with proposed and traditional algorithms are the same within the expected statistical scatter.

To stay the effect of non-conservative nature of traditional collisional procedure in the real conditions, some test computations of considered flow were made by AMC algorithm for adiabatic tube wall with diffuse-specular scattering of molecules for different values of accommodation coefficient $`\sigma `$. Fig.5 illustrates the dependencies of $`I/I_0`$ on the time for $`VSS`$ molecules, $`Kn=\mathrm{\hspace{0.17em}0.1}`$, $`S_\phi =\mathrm{\hspace{0.17em}1}`$, $`N=\mathrm{\hspace{0.17em}50}`$ and $`4`$ values of $`\sigma `$ in the range $`\sigma =\mathrm{\hspace{0.17em}2}10^3\mathrm{\hspace{0.17em}10}^4`$. The comparison of these dependencies with those presented in Fig.1 reveals that the traditional collisional procedure causes the same change in angular momentum as the diffuse-specular wall with the accommodation coefficient $`\sigma `$ in the range $`\sigma =\mathrm{\hspace{0.17em}10}^3\mathrm{\hspace{0.17em}10}^4`$. Therefore, for flows with real surface $`(\sigma \mathrm{\hspace{0.17em}1})`$ in the flowfield the studied effect is small, its relative value is about $`10^3\mathrm{\hspace{0.17em}10}^4`$.

## 5 Conclusion

The proposed algorithm of collisional process provides precise conservation of angular momentum and energy and may be used for direct simulation of flows with rotation.

It should be noted, however, that the improvement of the results, obtained by this algorithm may be important only for the problems without external source or sink of angular momentum. If solid surface with diffuse reflection is present in the flowfield, the changes in angular momentum, caused by this surface will be much greater than those, caused by non-conservative nature of traditional algorithm of collisional process. For such kind of flows there is no need to employ AMC algorithm, especially taking into account that it takes $`10\%\mathrm{\hspace{0.17em}15}\%`$ more CPU time.

Nevertheless, for some problems the application of AMC algorithm is advisable. An example of such kind of problem is the Ranque effect , which nature is determined by conservation of angular momentum of the gas flow in the vortex tube .


text:
# Estimating mixing properties of local Hamiltonian dynamics and continuous quantum random walks is PSPACE-hard

## 1 Introduction

Ergodic theory studies the long-term average behavior of dynamical systems. The continuous time-evolution is described by a one-parameter family $`\{T_t:t\}`$ of maps of the phase space $`X`$ into itself. A major topic of classical ergodic theory is to examine qualitatively how the orbits $`\{T_tx\}`$ (where $`x`$ is the initial state) penetrate the phase space $`X`$. If the orbits “spread uniformly” over $`X`$, then the system is called ergodic. Ergodicity of dynamical systems plays an important role in statistical mechanics and thermodynamics because in the ergodic case the time averages and ensemble averages of physical variables are equal (see  for a rigorous study of these questions).

In this paper we consider the case that $`X`$ is the Hilbert space of many qubits, and $`\{T_t:t\}`$ is the unitary one-parameter group $`\{e^{iHt}:t\}`$ generated by the system Hamiltonian $`H`$. We study the computational complexity of estimating the von Neumann entropy generated by Hamiltonian time evolutions if we average over the time. The motivation for considering the entropy of the time average is that it quantifies in a certain sense what portion of the state space and how uniformly it is filled out by the trajectory according to the time evolution. A second interesting property of the entropy of the time-average is that it could be used to quantify how far a pure state $`|\mathrm{\Phi }`$ is from all eigenvectors of a Hamiltonian (this question was raised in ). The reason is that the entropy of the time-average of the initial pure state $`|\mathrm{\Phi }\mathrm{\Phi }|`$ can only be large if the initial state is a superposition of many eigenvectors (corresponding to different eigenvalues) of $`H`$ with equal amplitudes. We explain this later in more detail.

The time-average is defined as follows:

###### Definition 1 (Time average)

Let $`H`$ be the Hamiltonian of a quantum system $``$. For an arbitrary state $`\rho `$ on $``$ its time average is the state $`\overline{\rho }`$ defined by

$$\overline{\rho }:=\underset{T\mathrm{}}{lim}\frac{1}{T}_{t=0}^Te^{iHt}\rho e^{iHt}dt.$$
(1)
It is well-known that the time average depends on the spectral decomposition of the Hamiltonian (a proof of the following lemma is included in the appendix for completeness).

###### Lemma 1

The time-average $`\overline{\rho }`$ of $`\rho `$ is given by

$$\overline{\rho }=\underset{\lambda }{}P_\lambda \rho P_\lambda ,$$

where $`P_\lambda `$ are the orthogonal projections onto the eigenspaces of $`H`$ for different eigenvalues $`\lambda `$.

Before we can formulate question on the complexity of estimating certain properties of time averages we have to restrict the class of considered Hamiltonians. This is because physically realistic Hamiltonians must satisfy some locality condition. One way to formalize this is the following definition.

###### Definition 2 (local operator)

Let $`:=_1\mathrm{}_n`$ be the tensor product Hilbert space of $`n`$ Hilbert spaces. We call an operator $`k`$-local if it is a sum of operators that act on at most $`k`$ tensor components non-trivially.

If each tensor component is the state space of a physical particle, then the natural interactions are usually pair-interactions, that is, $`2`$-local operators in the sense of the above definition.

Nevertheless, for the Hilbert space $`:=(^2)^n`$ of $`n`$ interacting qubits it may also be physically reasonable to allow $`k`$-local interaction for $`k>2`$. One reason is that they may describe effective Hamiltonians. A second reason is that there is not necessarily a one-to-one correspondence between qubits and physical particles (for instance, several qubits could describe the state of one particle).

We consider time-evolutions according to $`4`$-local Hamiltonians and continuous quantum random walks. The time evolution $`e^{iHt}`$ is called a continuous quantum random walk if $`H`$ is the adjacency matrix of some graph . Recall that a matrix $`A`$ is called an adjacency matrix if it is symmetric, has only $`0`$s and $`1`$s as entries, and all its entries on the diagonal are $`0`$. Furthermore, we say that a matrix is a $`k`$-local adjacency matrix if it is a sum of $`k`$-local operators that are adjacency matrices themselves. We will consider $`5`$-local adjacency matrices.

We show that the problem estimating the von-Neumann entropy of time-averages of computational basis states, where the time-evolution is described by $`4`$-local Hamiltonians and $`5`$-local quantum random walks, is PSPACE hard. The problem of estimating the density matrix of the time-average reduced to a qubit is also PSPACE-hard.

Due to Lemma 1 it is clear that the entropy generation depends on the whole spectrum of the Hamiltonians. Therefore, these results should be compared to the results on the complexity of estimating the minimal eigenvalue (ground state energy) of local Hamiltonians. The problem of determining the lowest energy value of a (classical) spin-spin interaction of Ising type is known to be NP-complete . For interacting qubits the problem of determining the lowest energy value is even QMA-complete (“Quantum-NP”) if one allows $`3`$-local interactions only . Recall that in these NP and Quantum-NP problems the task is not to determine the lowest eigenvalues with high precision. The demanded accuracy is only inverse polynomially in the number $`n`$ of interacting qubits.

Note that the entropy of the time-average of a pure state $`\rho :=|\mathrm{\Phi }\mathrm{\Phi }|`$ could be used to measure how “far away” $`|\mathrm{\Phi }`$ is from all eigenstates of $`H`$. This question was studied in ; the authors constructed quantum states that are far away from any eigenstate of any non-trivial local Hamiltonian, in the sense that $`|\mathrm{\Phi }|\mathrm{\Psi }`$ is greater than some constant lower bound for all eigenstates $`|\mathrm{\Psi }`$ of all $`k`$-local (for a fixed $`k`$) Hamiltonians, independent of the form of the Hamiltonians. The measure for the distance of a pure quantum state $`|\mathrm{\Phi }`$ from the spectrum of a Hamiltonian used in  is

$$d(|\mathrm{\Phi },H):=\underset{\{|\mathrm{\Psi }_\lambda \}}{\mathrm{min}}\underset{\lambda }{\mathrm{min}}|\mathrm{\Phi }|\mathrm{\Psi }_\lambda ,$$

where $`\{|\mathrm{\Psi }_\lambda \}`$ runs over all sets consisting of eigenvectors of $`H`$ corresponding to different eigenvalues $`\lambda `$. One easily verifies that the distance is equivalent to

$$\underset{\{|\mathrm{\Psi }_\lambda \}}{\mathrm{min}}\underset{\lambda }{\mathrm{min}}2(1|\alpha _\lambda |^2),$$

where $`|\mathrm{\Phi }=_\lambda \alpha _\lambda |\mathrm{\Psi }_\lambda `$. Therefore, the bound $`d(|\mathrm{\Phi },H)`$ is maximal if and only if $`|\mathrm{\Phi }`$ is an equally weighted superposition of eigenvectors corresponding to different eigenvalues of $`H`$. Using Lemma 1 one see that this is equivalent that the entropy of the time average of $`|\mathrm{\Phi }`$ is maximal; the maximally possible entropy is the logarithm of the number of different eigenvalues of $`H`$. Based on these observations, our results may be also interpreted as a proof of the computational difficulty of deciding for a given quantum state how far it is from all eigenstates of a local Hamiltonian.

## 2 Stating the problem

The techniques we use here to prove PSPACE-hardness of estimating the entropy of time-averages are closed related to those in  where we studied the complexity of measuring local observables. The general idea is to describe the computational steps carried out by a (reversible) space-bounded Turing machine by a quantum circuit consisting of elementary quantum gates or Toffoli gates and then to encode the circuit into a local Hamiltonian (this will be explained in Sections $`3`$ and $`4`$).

###### Theorem 1 (PSPACE-hardness)

Let $`(A_l)`$ be a uniformly generated family of $`4`$-local Hamiltonians (or $`5`$-local adjacency matrices) on the Hilbert spaces $`H_l`$ consisting of $`poly(l)`$ qubits. Then the problems of (1) estimating the von-Neumann entropy of the time-average $`S(\overline{\rho })`$ and (2) estimating the reduced state of $`\overline{\rho }`$ to single qubits are PSPACE-hard.

More precisely: (1) Deciding if the entropy of the time-average of a computational basis state of $`_l`$ according the dynamics given by $`A_l`$ is smaller than $`a_l`$ or greater than $`b_l`$, where $`b_la_l>poly(l)`$, is PSPACE-hard.

Let $`\overline{\rho }_i`$ denote the reduced state of the time-average $`\overline{\rho }`$ to the $`i`$th qubit of $`_l`$. (2) Deciding if $`1|\overline{\rho }_i|1<1/d_l`$ or $`1|\overline{\rho }_i|1>1/21/d_l`$, where $`d_l=poly(l)`$, is PSPACE-hard.

We prove this theorem in Section 4. The proof is based on a characterization of the complexity class PSPACE with quantum circuits that is explained in the next section.

## 3 Characterizing PSPACE by circuits

The complexity class PSPACE is usually defined with respect to the Turing machine model . PSPACE is the class of all languages recognizable by polynomial space bounded deterministic Turing machines that halt on all inputs .

In the following we work with a characterization of PSPACE with respect to circuits consisting of elementary quantum gates (Toffoli gates only). This characterization is proved in Theorem 1 in .

###### Theorem 2 (PSPACE)

For every language $`L`$ in PSPACE there is a polynomial-time uniformly generated family of quantum circuits $`(V_l)_l`$. Here $`l`$ denotes the length of the input $`x`$. Each $`V_l`$ consists of $`s_l=poly(l)`$ elementary quantum gates (or Toffoli gates111In  we have used elementary quantum gates instead of Toffoli gates. But it is clear that the transitions performed reversible Turing machines can be simulated by a circuit consisting of Toffoli gates only. This is because the Toffoli gate is universal for reversible computation. The reason why we use the Toffoli gates is that it is possible to construct from such a circuit a Hamiltonian is an adjacency matrix.) and acts on $`m_l=poly(l)`$ many qubits. The circuit $`V_l`$ decides whether an input string $`x`$ is an element of $`L`$ in the following sense.

There is a polynomial-time computable natural number $`r_l`$ such that the $`r_l`$-fold concatenation of $`V_l`$ solves the corresponding PSPACE problem, i.e.

$$V_l^{r_l}(|x|y|00\mathrm{}0)=|x|yf(x)|00\mathrm{}0,$$

where $`f`$ is the characteristic function of $`L`$. That is $`f(x)=1`$ if $`xL`$ and $`f(x)=0`$ otherwise. The vector $`|x`$ is the basis state given by the binary word $`x\{0,1\}^l`$, the vector $`|y`$ is the state of the output qubit and $`|00\mathrm{}0`$ is the initial state of $`m_ll1`$ ancilla qubits.

The intuition behind this construction is as follows. Recall that irreversible TMs (used to define PSPACE) can be simulated space-efficiently by reversible TMs . The computational steps of a reversible TM can be simulated efficiently by a quantum circuit $`U`$. This circuit acts on registers representing the state of the head of TM, the current tape position, and a polynomial portion of the tape (see  for an explicit construction of the quantum circuit $`U`$) and realizes permutations of basis states corresponding to the transitions of the TM. By augmenting $`U`$ with some control logic, we can achieve that once the TM moves into a final state, that the new circuit performs some idle cycles, flips the output qubit if and only if the answer is yes after exactly $`r_l/2`$ steps, then performs some idle cycles, and reverses the computation. This is schematically illustrated in Figure 1. It is necessary to perform idle cycles because we sometime do not know exactly after how many steps the TM moves into a final step, but only know a lower bound. By including the idle steps we guarantee the computational time depends only on the input length.

###### Remark 1

We can consider an output register instead of an output qubit. It is also possible to construct a quantum circuit such that

$$V_l^{r_l}(|x|y|00\mathrm{}0)=|x|y+f(x)mod2^c|00\mathrm{}0,$$

for any $`c`$ that is polynomial in $`l`$.

## 4 Constructing the Hamiltonians /
the quantum random walks

Starting from the family $`(V_l)`$ of circuits we construct a uniformly generated family of $`4`$-local Hamiltonians ($`5`$-local adjacency matrices) $`(A_l)_l`$ such that the time-average of the state

$$|x|0|00\mathrm{}0|100\mathrm{}0$$

and its entropy encodes the solution the instance $`x`$ of a PSPACE problem; the first three tensor components are as in Theorem 2 and the last one corresponds to the new register clock that is needed in our construction. The construction is based on Feynman’s construction of a computer whose dynamics is an autonomous time-evolution .

To explain the construction of the Hamiltonian (adjacency matrix $`A`$) let $`V`$ be a circuit as in Theorem 2 and $`s`$ be its size $`s`$, that is, the number of elementary gates (Toffoli gates). We need a register clock indicating which gate is applied. It consists of $`s`$ qubits. The allowed states of the register clock are of the form $`|0\mathrm{}010\mathrm{}0`$ indicating which gate of $`V`$ is applied currently. We denote by $`T_k`$ the elementary quantum gates (Toffoli gates) of $`V`$.

We first define the forward-time operator

$`F`$ $`=`$ $`T_0|1_20|_1|0_11|_0+`$
$`T_1|1_30|_2|0_21|_1+`$
$`\mathrm{}`$
$`T_{s1}|1_10|_0|0_s1|_{s1},`$
where the operators $`|0_k1|_k`$ and $`|1_k0|_k`$ are the annihilation $`a`$ and creation operators $`a^{}`$, respectively, acting on the $`i`$th qubit of the clock. The backward-time operator is defined as the adjoint of $`F`$. The Hamiltonian is defined as

$$A:=F+F^{}.$$

Now we show that if $`V`$ consists of Toffoli gates only, then $`A`$ is a $`5`$-local adjacency matrix of some graph in the computational basis, that is, $`A`$ is symmetric, has only $`0`$s and $`1`$s as entries, and its entries on the diagonal are all $`0`$s. Since the matrix describing the action of the Toffoli gate (in the computational basis) contains only $`0`$s and $`1`$s, the sums

$$T_k|1_k0|_k|0_{k+1}1|_{k+1}+T_k|0_k1|_k|1_{k+1}0|_{k+1}$$

are adjacency matrices that are $`5`$-local operators. $`A`$ is the sum of the above matrices. Therefore, it is symmetric and has only $`0`$s on the diagonal. It remains to show that no $`1`$s of these matrices can meet when summing. This is done by checking that the operators

$$a_k^{}a_{k+1mods},a_ka_{k+1mods}^{},$$

lead to orthogonal states for all $`k=0,\mathrm{},s1`$ when applied to any computational basis state. We see this by observing that if $`klmods2`$ then we have

$$b|a_k^{}a_{k+1}a_la_{l+1}^{}|b=0\text{ and }b|a_ka_{k+1}^{}a_la_{l+1}^{}|b=0$$

because the operators $`a_k^{}a_{k+1}a_la_{l+1}^{}`$ and $`a_ka_{k+1}^{}a_la_{l+1}^{}`$ either map $`|b`$ onto the zero vector or change the bits of $`|b`$ at the non-overlapping positions $`(k,k+1)`$ and $`(l,l+1)`$. Similarly, if $`l=k+1`$ then we have

$$b|a_k^{}a_{k+1}a_{k+1}a_{k+2}^{}|b=0\text{ and }b|a_ka_{k+1}^{}a_{k+1}a_{k+2}^{}|b=0$$

because $`a_k^{}a_{k+1}a_{k+1}a_{k+2}^{}`$ is the zero operator and $`a_ka_{k+1}^{}a_{k+1}a_{k+2}^{}`$ either maps $`|b`$ onto the zero vector or changes the bits of $`|b`$ at the non-overlapping positions $`k`$ and $`k+2`$.

Therefore, only computational basis states appear in $`A|b`$ with coefficients all equal to $`1`$. Consequently, we have either $`b|A|b=0`$ or $`b|A|b=1`$. This show that $`A`$ is a $`5`$-local adjacency matrix if $`V`$ contains Toffoli gates only.

We denote the linear span of the vectors

$$F^j|\mathrm{\Psi }_0\text{for }j$$

where $`|\mathrm{\Psi }_0:=|x\mathrm{\hspace{0.17em}00}\mathrm{}0|100\mathrm{}0`$ as $`𝒪`$. Here the first part of the tensor product denotes the register where $`V`$ acts on (input and ancilla registers) and the second component is the register clock. All states of this orbit are orthogonal until one has a recurrence to the initial state $`|\mathrm{\Psi }_0`$. This can be seen as follows: If the register clock is in an allowed state there is only one summand of $`F`$ that is relevant. Its action on the clock is simple since it moves the symbol $`1`$ to the next qubit (this dynamics of the clock may be interpreted as a propagation of a spin-wave). Therefore it is clear that the first $`s1`$ states are orthogonal. The whole circuit $`V`$ is a classical logical operation which permutes basis states. Therefore the state $`F^s|\mathrm{\Psi }_0`$ is either orthogonal to $`|\mathrm{\Psi }_0`$ or both states coincide. Along the same line we argue that all states of the orbit are orthogonal until a state coincides with the initial state. Hence $`F`$ acts as a cyclic shift on $`𝒪`$.

For our construction is is essential that the dimension of the orbit depends on the solution of the PSPACE problem. It is $`2sr`$ if $`f(x)=1`$ and $`sr`$ if $`f(x)=0`$, where $`f(x)`$ is the answer for the instance $`x`$ of the considered PSPACE problem. We denote the dimension of the orbit $`𝒪`$ by $`d`$. Let $`\omega `$ be a primitive complex $`d`$-th root of unity. The eigenvalues of $`F`$ restricted to $`𝒪`$ are

$$\omega ^0,\omega ^1,\omega _2,\mathrm{},\omega ^{d1}.$$

Furthermore, $`|\mathrm{\Psi }_0`$ is a superposition of all eigenvectors of $`F`$ restricted to $`𝒪`$ with equal amplitudes. All this follows from properties of the cyclic shift operator. Since $`F`$ and $`F^{}`$ commute on $`𝒪`$ the eigenvalues of $`A`$ restricted to $`𝒪`$ are $`\omega ^k+\overline{\omega }^k=2\mathrm{cos}(2\pi j/d)`$ for $`k=0,\mathrm{},d1`$.

Set $`\rho :=|\mathrm{\Psi }_0\mathrm{\Psi }_0|`$ and $`\lambda _k:=2\mathrm{cos}(2\pi k/d)`$. Note that $`d`$ is always even; this is because $`r`$ in Theorem 2 is even. Therefore, there are two eigenvalues with multiplicity $`1`$ (corresponding to the real eigenvalues $`+1`$ and $`1`$ of $`F`$) and $`d2`$ eigenvalues with multiplicity $`2`$ (corresponding to the complex eigenvalues of $`F`$). We have

$$\rho =\frac{1}{d}\underset{k,l=0}{\overset{d1}{}}|kl|,$$

where $`|k`$ denotes the $`k`$th eigenvalue of $`F`$ restricted to $`𝒪`$. For the time-average we obtain with Lemma 1

$$\overline{\rho }=\frac{1}{d}\underset{k,l=0}{\overset{d1}{}}[\lambda _k=\lambda _l]|kl|,$$

where $`[\lambda _k=\lambda _l]`$ is $`1`$ if both eigenvalues are equal and $`0`$ otherwise. If we permute the order of the eigenvectors $`|k`$ such that eigenvectors corresponding to the same eigenvalues are adjacent, then $`\rho `$ looks in this basis as follows:

$$\frac{1}{d}\left(\begin{array}{ccccccc}1& & & & & & \\ & 1& & & & & \\ & & 1& 1& & & \\ & & 1& 1& & & \\ & & & & \mathrm{}& & \\ & & & & & 1& 1\\ & & & & & 1& 1\end{array}\right)$$
(2)
Using the form of the time-average in (2), its entropy is easily computed

$$S(\overline{\rho })=\mathrm{log}_2d(d2)/d\mathrm{log}_2d1.$$

Therefore, the difference in entropy of the time-average of the computational basis state $`|\mathrm{\Psi }_0`$ encoding a “yes”-instance and a “no”-instance is at least $`c/2`$ bits, where $`c`$ is the size of output register (see Remark 1). This proves the PSPACE-hardness of estimating the entropy of the time-average as stated in Theorem 2.

Let the output register consist of only one qubit. Obviously, the reduced state of time-average to the output qubit is $`|00|`$ for a “no”-instance. For a “yes”-instance we know that the output qubit of the states $`|\mathrm{\Psi }_j:=F^j|\mathrm{\Psi }_0`$ is $`1`$ for $`j=d/4,\mathrm{},3d/41`$ and $`0`$ for all other $`j`$. To see this, recall how $`V`$ works: at step $`r_l/2`$ the circuit $`V`$ changes the output qubit from $`0`$ to $`1`$, at step $`r_l`$ we obtain the initial state with the output flipped, at step $`r_l+r_l/2`$ the output is flipped from $`1`$ to $`0`$, and finally in step $`2r_l`$ we obtain the original input state.

Therefore, the occupation probability of the level $`|1`$ of $`\overline{\rho }_{\mathrm{out}}`$ ($`\overline{\rho }`$ restricted to the output qubit) is exactly

$`1|\overline{\rho }_{\mathrm{out}}|1`$ $`=`$ $`{\displaystyle \underset{j=d/4}{\overset{3d/41}{}}}\mathrm{\Psi }_j|\overline{\rho }|\mathrm{\Psi }_j.`$
Let us denote by $`P(j)`$ the probability to observe the state $`|\mathrm{\Psi }_j`$ if we measure the time-average in the computational basis. The probability $`P(j)`$ is given by

$`P(j)`$ $`=`$ $`\mathrm{\Psi }_j|\overline{\rho }|\mathrm{\Psi }_j`$
$`=`$ $`{\displaystyle \frac{1}{d^2}}{\displaystyle \underset{k,l=0}{\overset{d1}{}}}\omega ^{(kl)j}[\lambda _k=\lambda _l]`$
$`=`$ $`{\displaystyle \frac{1}{d}}+{\displaystyle \frac{1}{d^2}}{\displaystyle \underset{kl}{\overset{d1}{}}}\omega ^{(kl)j}[\lambda _k=\lambda _l]`$
This seen by observing that $`|\mathrm{\Psi }_j=\frac{1}{\sqrt{d}}_{k=0}^{d1}\omega ^{jk}|k`$. Now we make use of the Diaconis-Shahshahani bound (see Appendix) applied to the cyclic group $`_d`$. The characters of $`_d`$ are given by $`\chi _a(j)=\omega ^{aj}`$. The coefficients of the Fourier transform of the probability density $`P=(P(0),\mathrm{},P(d1))`$ are

$`\widehat{P}(a)`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{d}{}}}P(j)\chi _a(j)`$
$`=`$ $`{\displaystyle \frac{1}{d^2}}{\displaystyle \underset{kl:\lambda _k=\lambda _l}{}}{\displaystyle \underset{j=0}{\overset{d1}{}}}\omega ^{(kl+a)j}`$
Note that the possible values for $`kl`$ are only the even numbers $`2,4,\mathrm{},d/2`$ since we sum over all pairs $`(k,l)`$ satisfying $`kl`$ and the condition $`\mathrm{cos}(2\pi k/d)=\mathrm{cos}(2\pi l/d)`$ being equivalent to $`k=j`$.

If $`a`$ is odd, then $`kl+a0modd`$. But if $`a`$ is even, then there is exactly one pair of $`(k,l)`$ satisfying above properties and $`kl+a0modd`$. With these observations we obtain

$$\widehat{P}(a)=\{\begin{array}{cc}0& \text{ if }a\text{ odd}\\ \frac{1}{d}& \text{ if }a\text{ even}\end{array}$$
(3)
Now we obtain with Diaconis-Shahshani bound a lower bound on the total variation distance between $`P`$ and the uniform distribution $`U`$ on the orbit states $`F^j|\mathrm{\Psi }_0`$

$$PU_{\mathrm{TV}}^2\frac{1}{4}(d/21)\frac{1}{d^2}\frac{1}{8d}$$

Taking the square root and multiplying both by $`2`$ we obtain

$$\underset{j=0}{\overset{d1}{}}|P(j)\frac{1}{d}|\frac{1}{\sqrt{2d}}.$$

Using the triangle inequality the above bound implies that

$$1|\overline{\rho }_{\mathrm{out}}|1=1|\overline{\rho }|1=\underset{j=d/4}{\overset{3d/41}{}}P(j)[\frac{1}{2}\frac{1}{\sqrt{2d}},\frac{1}{2}+\frac{1}{\sqrt{2d}}].$$

(Note that this is also true if we sum over any set containing half of all possible $`j`$.)

In summary, we have proved that for a “yes”-instance the occupation probability $`1|\overline{\rho }_{\mathrm{out}}|1`$ is almost $`1/2`$, whereas it is $`0`$ for a “no”-instance.

## 5 Appendix

Proof. Lemma 1: Let $`\{|k\}`$ be a basis consisting of eigenvectors of $`H`$ with corresponding eigenvalues $`\lambda _k`$. Express the state $`\rho `$ in this basis

$$\rho =\underset{k,l}{}\gamma _{kl}|kl|.$$

For the time-average we obtain

$`\overline{\rho }`$ $`=`$ $`\underset{T\mathrm{}}{lim}{\displaystyle \frac{1}{T}}{\displaystyle _{t=0}^T}{\displaystyle \underset{k,l=0}{\overset{d1}{}}}e^{i(\lambda _k\lambda _l)t}\gamma _{kl}|kł|`$
$`=`$ $`{\displaystyle \underset{k,l=0}{\overset{d1}{}}}\left(\underset{T\mathrm{}}{lim}{\displaystyle \frac{1}{T}}{\displaystyle _{t=0}^T}e^{i(\lambda _k\lambda _l)t}\right)\gamma _{kl}|kł|`$
$`=`$ $`{\displaystyle \frac{1}{d}}{\displaystyle \underset{k,l=0}{\overset{d1}{}}}[\lambda _k=\lambda _l]\gamma _{kl}|kl|`$
$`=`$ $`{\displaystyle \underset{\lambda }{}}P_\lambda \rho P_\lambda ,`$
where $`[\lambda _k=\lambda _l]`$ is $`1`$ if both eigenvalues are equal and $`0`$ otherwise, $`\lambda `$ are the different eigenvalues of $`H`$, and $`P_\lambda `$ are the orthogonal projection onto the corresponding eigenspaces. $`\mathrm{}`$

The proof of the lemma can be found in .

###### Lemma 2 (Diaconis-Shahshahani bound)

Let $`A`$ be an arbitrary abelian group, $`P`$ an arbitrary probability distribution on $`A`$ and $`U`$ the uniform distribution on $`A`$. The Diaconis-Shahshahani bound on the total variation distance

$$PU_{\mathrm{TV}}:=\frac{1}{2}\underset{gA}{}|P(g)U(g)|$$

is

$$PU_{\mathrm{TV}}^2\frac{1}{4}\underset{\chi \chi _0}{}|\widehat{P}(\chi )|^2,$$
(4)
where the sum runs over all non-trivial irreducible characters of $`A`$ and $`\widehat{P}`$ denotes the Fourier transform of $`P`$, that is, $`\widehat{P}(\chi )=_{gA}\chi (a)P(g)`$.

### Acknowledgements

We would like to thank Dominik Janzing for proposing us to consider the problem of estimating the entropy of time-averages and many interesting discussions.

This work has been supported by grants of the Landesstiftung Baden-Württemberg (project “Kontinuierliche Modelle der Quanteninformationsverarbeitung”).


text:
# X-Ray Timing of the Young Pulsar in 3C 58

## 1 Introduction

PSR J0205$`+`$6449 was discovered in the heart of 3C 58 in Chandra HRC data and immediately confirmed in archival RXTE data by Murray et al. Murray et al. (2002). Camilo et al. (Camilo et al., 2002) soon thereafter found extremely faint radio pulsations from the pulsar using the new Green Bank Telescope (GBT). Chandra ACIS-S observations taken by Slane, Helfand, & Murray (Slane et al., 2002) showed that the temperature of the neutron star (NS) is significantly below that predicted by standard NS cooling models given the very young age of the pulsar ($`820`$ yrs). Measurement of the spin period (65.686 ms) and spin-down rate ($`\dot{P}=1.94\times 10^{13}s/s`$) of the pulsar have allowed the estimation of several key parameters using the canonical expressions from pulsar theory:

* Spin-down energy

$$\dot{E}=4\pi ^2I\dot{P}/P^3=2.7\times 10^{37}\mathrm{erg}/\mathrm{s}$$
(1)
* Surface magnetic field strength

$$B=\sqrt{\frac{3c^3IP\dot{P}}{8\pi ^2R^6}}=3.6\times 10^{12}\mathrm{Gauss}$$
(2)
* Characteristic age

$$\tau =P/(2\dot{P})=5400\mathrm{yr}$$
(3)

Each of these parameters is typical for young pulsars. One unusual point is that if the pulsar really is 822 yrs old and has been spinning down with the canonical braking index (due to magnetic dipole braking) of $`n=3`$, then it must have been born spinning relatively slowly with $`P_060`$ ms. Interestingly, the 65-ms PSR J1811$``$1925 in G11.2$``$0.5 seems to be an identical twin in many respects (see Roberts et al., these proceedings).

## 2 Observations and Data Preparation

We analyzed all available RXTE PCA data on 3C 58, which totaled more than $``$410 ks of on-source time. The data came from OBSIDs P20259 ($``$20 ks), P60130 ($``$100 ks), and P70089 ($``$300 ks), and were taken in GoodXenon mode which provided events with time resolutions of $``$122 $`\mu `$s and 256 energy channels in the energy range $``$2$``$60 keV. Most observations were taken using 3 PCUs.

We processed the RXTE data using standard FTOOLs111http://heasarc.gsfc.nasa.gov/docs/software/ftools/ftools\_menu.html in order to filter out non-Level 1 events or those recorded outside of Good Time Intervals and to convert their arrival times to the Solar System barycenter (using the DE200 planetary ephemeris (Standish, 1982)). We determined the X-ray pulse times-of-arrival (TOAs) with a maximum-likelihood technique that we developed which assumed a two-Gaussian model of the X-ray pulse profile. A signal-to-noise ratio optimization procedure allowed us to determine that the best energy range to use for the timing analysis is $``$2$``$16 keV (energy channels 4$``$39 inclusive), and so we filtered out events outside of this energy range (see the profile in Figure 1). A single time-of-arrival (TOA) typically comprised one or two RXTE orbits depending on the background levels and the number of PCUs on. Since the observations consisted of between 4$``$6 orbits, we measured 3$``$4 TOAs each month with precisions of 100$``$200 $`\mu `$s each.

The GBT data were taken with the Berkeley-Caltech Pulsar Machine (BCPM)222http://www.gb.nrao.edu/~dbacker backend at either 820 MHz or 1400 MHz. The BCPM is an analog/digital filterbank which samples each of 2$`\times `$96 channels using 4-bits at flexible sampling rates and channel bandwidths (Backer et al., 1997). The data were recorded using 134 MHz of bandwidth and 50 $`\mu `$s samples for the 1400 MHz observations or 50 MHz of bandwidth and 72 $`\mu `$s samples for the 820 MHz observations. Typical integrations lasted between 4$``$8 hrs. We de-dispersed and folded all of the data using the pulsar analysis package PRESTO (Ransom, 2001). Each observation yielded 2$``$3 TOAs which we measured by correlating (in the frequency domain) the folded pulse profile with a Gaussian of fractional width 0.04 in phase. The typical precision of the radio TOAs was 200$``$400 $`\mu `$s.

## 3 Timing

The timing behavior and therefore the spin-down rate of PSR J0205$`+`$6449 is extremely noisy. Determining a phase-connected timing solution would have been impossible using the RXTE or GBT data alone. Thankfully, the cadence of the monthly observations in the radio and X-ray bands was sufficiently different such that the observations complemented each other and usually allowed unambiguous pulse numbering between observations.

We applied time corrections to transform all the TOAs to UTC at the Solar System Barycenter using faxbary and TEMPO333http://pulsar.princeton.edu/tempo for the RXTE and GBT data respectively. This absolute timing shows that the X-ray main pulse is delayed relative to the radio pulse (Camilo et al., in prep). We performed all timing model fits with TEMPO.

Our data show that a “giant” glitch ($`\mathrm{\Delta }\nu /\nu `$1$`\times `$10-6) occurred near MJD 52555 (see Figure 2). The timing solution prior to the glitch consists of a 6 parameter fit (spin frequency $`\nu `$ plus 5 frequency derivatives) resulting in RMS timing residuals of $``$0.6 ms. The post-glitch timing solution consists of the spin frequency plus 2 frequency derivatives resulting in RMS timing residuals of $``$1.2 ms. Additional frequency derivatives in the post-glitch solution improve the residuals but are not required for phase connection444While we believe that the post-glitch data is phase-connected, the sparse sampling and high levels of timing noise make it possible that we have miscounted pulses between one or more of the observations.. The residuals from the GBT and RXTE for both timing solutions are shown in Figure 3.

## 4 Phase-Resolved Spectroscopy

Figure 4 shows the results of absorbed power-law fits to the phase-resolved spectra of the pulsar’s main pulse and interpulse as determined using the FTOOL fasebin. We analyzed events from each PCU independently for each OBSID, combined these for each AO, and then combined the three AOs together. We fit the spectra with xspec after subtracting the average off-pulse background as shown in the hatched region in Figure 1 and freezing the value of $`N_H=3.62\times 10^{21}`$ cm-2 based on the recent Chandra ACIS observation of 3C 58 (Slane et al., in prep). Fits made to the summed PCU data and simultaneously to the 5 PCUs independently were consistent within the statistical errors to one another.

The measured photon indices $`\mathrm{\Gamma }=0.84_{0.15}^{+0.06}`$ for the main pulse and $`\mathrm{\Gamma }=1.0_{0.3}^{+0.4}`$ for the interpulse are surprisingly hard. These values are difficult to reconcile with the $`\mathrm{\Gamma }=1.73\pm 0.07`$ photon index measured with Chandra for the central point source (Slane et al., 2002), although significant nebular emission within the selection region could explain some of the spectral softening. Recent results from the 350 ks Chandra observation of 3C 58 place an approximate upper limit of $`1.5\times 10^{14}`$ erg cm-2 s-1 for any hard component in the range 0.5$``$10 keV. But our fits to the main pulse alone predict fluxes in the same energy band at least one hundred times larger. We note that our measurements are consistent within the errors with the phase-resolved spectral fits of Gotthelf (Gotthelf, 2003) ($`\mathrm{\Gamma }=1.11\pm 0.34`$ for the pulsed emission) based on the RXTE data from P60130.

## 5 Future Prospects

With its young age, low temperature, very noisy spin-down, high $`\dot{E}`$, hard pulsed spectrum, and large glitches, the pulsar at the center of 3C 58 is certainly an object worth watching in the future. Is the extremely high level of timing noise from the pulsar somehow related to its “cool” temperature? Are both of these strange properties related to a more fundamental characteristic of the neutron star such as its mass? We are optimistic that continued monitoring of this unusual source will tell us something about the interior structure and internal dynamics of neutron stars. In addition, we note that the narrow pulse profile and hard spectrum of the pulsed emission from J0205$`+`$6449 make it an ideal and possibly unique target for INTEGRAL and other upcoming hard X-ray missions.


text:
# Theoretical Fe i/ii/iii Emission-Line Strengths from Active Galactic Nuclei with Broad-Line Regions

## 1 Introduction

The ultraviolet spectra of active galactic nuclei (AGN) with broad-line regions (BLRs) exhibit a quasi-continuum of thousands of blended iron emission lines, dominated by Fe II (Wills et al.  1980, Wills, Netzer, & Wills 1980); numerous Fe II transitions are also present in the optical (Boroson & Green 1992, Véron-Cetty, Joly & Véron 2004) and near-infrared (Rudy et al. 2000, Rudy et al. 2001, Rodríguez-Ardila et al. 2002). Since their recognition as a significant component of the BLR emission spectrum, these Fe II transitions have been the subject of intense observational scrutiny and theoretical modelling (Netzer 1980; Kwan & Krolik 1981; Netzer & Wills 1983; Wills, Netzer & Wills 1985; Elitzur & Netzer 1985; Collin-Souffrin et al. 1986; Penston 1987; Collin-Souffrin, Hameury & Joly 1988; Dumont & Collin-Souffrin 1990; Netzer 1990; Sigut & Pradhan 1998; Verner et al. 1999; Sigut & Pradhan 2003). Reproducing the strength of the Fe II emission has been a considerable challenge for traditional photoionized models for the BLR clouds; typically the Fe II/H$`\beta `$ flux ratio is observed to be $`\mathrm{\hspace{0.17em}10}`$, ranging from $`\mathrm{\hspace{0.17em}2}`$ to near $`\mathrm{\hspace{0.17em}30}`$ in the case of super-strong Fe II emitters (Joly 1993; Graham, Clowes & Campusano 1996).

However Fe II is not the only component of the low-ionization emission spectrum of the BLR. Ultraviolet transitions of Fe III are also well established (Baldwin et al. 1996, Laor et al. 1997, Vestergaard & Wilkes 2001). Laor et al. (1997) identify a strong feature near $`\lambda \mathrm{\hspace{0.17em}2418}`$Å in the spectrum of the prototypical narrow-line Seyfert 1 (NLS1) galaxy I Zw 1 as Fe III multiplet UV47. Most recently, Vestergaard & Wilkes (2001) provide a detailed analysis of the UV Fe III emission from the I Zw 1 and use their observations to empirically derive Fe II and Fe III flux templates.

In our previous work (Sigut & Pradhan 2003, hereafter SP03), we constructed a detailed non-LTE radiative transfer model for Fe II, including 829 atomic levels, and predicted a large number ($`23,000`$) of Fe II emission line strengths for conditions typical of BLRs of AGN. In this paper, we extend our earlier work on Fe II to compute detailed non-LTE models for Fe i-iv. We compare these predictions to the observed UV I Zw 1 template.

## 2 Atomic Data

The four lowest ionization stages of iron, Fe i-iv, were explicitly included in the calculations (although Fe IV was included as only a one-level atom). A summary of the complete iron model atom is given in Table 1. The final non-LTE atomic model includes 944 energy levels and 14,962 radiative bound-bound (rbb) transitions. With the exception of level energies, the majority of radiative and collisional atomic data was computed using methods developed under the Iron Project111Complete references for the Iron Project are given at www.astronomy.ohio-state.edu/$``$pradhan (Hummer et al. 1993, IP.I). The R-matrix method employed in the Iron Project (IP) ensures relatively uniform accuracy for the computed data. In the following sections, the IP atomic data will also be identified according to the number of the paper in the ongoing IP series published in Astronomy and Astrophysics (e.g. IP.I for the first paper in this series).

### 2.1 Fe I

A minimal Fe I atomic model was included in the calculation, mainly to predict the ionization fraction of Fe I in the various BLR models and provide preliminary flux estimates. While rare, Fe I fluxes have been tentatively identified for a few Fe II-strong quasars (Boroson & Green 1992, Kwan et al. 1995). The Fe I atomic model consisted of 77 fine-structure levels (see Figure 1) comprising 19 low-energy LS levels from the triplet, quintet, and septet spin symmetries. Energy levels and Einstein spontaneous transition probabilities ($`A_{ji}`$ values) for the 185 transitions were taken from the NIST compilation222http://physics.nist.gov/PhysRefData/contents-atomic.html. Photoionization cross sections were adopted from the R-matrix calculations of Bautista & Pradhan (1997, IP.XX). Low-temperature, IP effective collision strengths for electron impact excitation of the first 10 meta-stable Fe I levels have been computed by Pelan & Berrington (1997, IP.XXI). All remaining collision strengths were estimated with the effective Gaunt factor approximation of van Regemorter (1962).

### 2.2 Fe II

Atomic data for Fe II has been extensively discussed by Sigut & Pradhan (1998, 2003). The IP data for collisions strengths for electron impact excitation were computed by Zhang & Pradhan (1995, IP.VI), for radiative transition probabilities, by Nahar (1995, IP.VII), and for photoionization cross sections, by Nahar & Pradhan (1994). We have adopted the smaller, 256 fine-structure level atom of Sigut & Pradhan (1998) for use in this calculation. The smaller number of Fe II levels allows a maximal number of Fe III energy levels to be simultaneously included. However, this smaller Fe II atom is still sufficiently large for realistic estimates of the flux and the inclusion of $`\mathrm{Ly}\alpha `$ fluorescent excitation.

Among the atomic data for Fe II, we call particular attention to the $`\mathrm{Fe}^{2+}+\mathrm{H}\mathrm{Fe}^++\mathrm{H}^+`$ charge-exchange reaction rates of Neufeld & Dalgarno (1987). These large rates, computed in the Landau-Zener approximation, play an important role in the Fe II-Fe III ionization balance. However, their accuracy and efficacy in determining emission line strengths can be a source of uncertainly in the present calculations, as discussed later.

### 2.3 Fe III

Fe III was represented by 581 fine-structure energy levels (see Figure 2) which includes all Fe III levels with known energies (see the NIST compilation; also Nahar & Pradhan 1996). We have used four principal sources to compile the required $`A_{ji}`$ values for Fe III transitions: the NIST critically-evaluated database; the R-matrix results of Nahar & Pradhan (1996, IP.XVII) for dipole-allowed ($`\mathrm{\Delta }S=0`$) transitions; the Hartree-Fock calculations of Ekberg (1993), which also treats inter-system ($`\mathrm{\Delta }S\mathrm{\hspace{0.17em}0}`$) transitions, and the extensive semi-empirical calculations of Kurucz (1992).

Photoionization cross sections were adopted from Nahar & Pradhan (1996, IP.XVII). Collision strengths for electron impact excitation were adopted from Zhang (1996, IP.XVIII). All remaining collision strengths were estimated with the Gaunt factor approximation. Charge exchange reactions with Fe IV were also included following Kingdon & Ferland (1996).

## 3 Calculations

The iron line fluxes were calculated using the procedure of SP03. Briefly, a background temperature and density structure for a single BLR cloud of a given ionization parameter and total particle density was computed with cloudy (Ferland, 1992). The shape of the photoionizing continuum was taken to be that of Mathews & Ferland (1987). Given this fixed, one-dimensional model, the coupled equations of radiative transfer and statistical equilibrium were solved for a self-consistent set of iron level populations and line fluxes. A complete discussion of the system of equations solved, including the treatment of fluorescent excitation by Lyman-$`\alpha `$ and Lyman-$`\beta `$, and the implementation details, are given in SP03. We note that our treatment of fluorescent excitation involves first the construction the full frequency-dependent, PRD (partial-redistribution) source functions for Lyman-$`\alpha `$ and Lyman-$`\beta `$, and then their inclusion in the total monochromatic source functions used in the radiative transfer solutions for all overlapping iron rbb transitions.

In the current work, each set of iron fluxes is specified by four model parameters representing the BLR cloud: the ionization parameter, the total particle density, the (internal) cloud turbulent velocity333Turbulent velocity is used to refer to what is commonly denoted microturbulence by stellar astronomers; this velocity represents the width of a Gaussian distribution of turbulent motions small in scale compared to unit optical depth. These motions act to broaden the atomic absorption profile and thus can alter the radiative transfer solution in each line., and the iron abundance relative to hydrogen. Table 2 identifies the different BLR models constructed in this work and the associated ionization parameters and particle densities.

## 4 Iron Ionization Balance

The ionization balance among the four lowest ionization stages of iron was self-consistently calculated simultaneously with the level populations. Expression for the total and level-specific photo and collisional ionization-recombination rates are given in Paper I. A small correction was applied at each depth for the fraction of iron more highly ionized than Fe IV, with this correction being estimated by cloudy.

Our self-consistent procedure is different from Verner et al. (1999) who first solve for the iron ionization balance in isolation, balancing photoionization from the ground states against effective recombination rates, and then solving for the individual (Fe II) level populations using a fixed ionization structure. While this approximation is typical of nebular work where the population in the excited states is very small, densities in the BLR may approach $`10^{+9}10^{+11}`$ particles per cm3, high enough to warrant a detailed treatment of photoionization from excited states in the calculation of the ionization balance. Another advantage of our approach is that no a priori assumptions are made as to the optical depths in the radiative transitions forming the recombination cascade (assumptions usually expressed as “Case A” or “Case” B, following Osterbrock, 1989, for example); a self-consistent solution is automatically obtained.

Figure 3 shows the Fe I-Fe IV ionization balance for all the models of Table 2 assuming a solar iron abundance and turbulent velocity of 10 $`\mathrm{km}\mathrm{s}^1`$(at $`10^4`$K, the thermal widths of the iron lines are $`2`$$`\mathrm{km}\mathrm{s}^1`$). At large optical depths in the Lyman continuum, Fe II is always the dominant iron species. However, closer to the illuminated face, the ionization parameter plays a strong role; for the lowest ionization parameter considered, Fe III dominates in the outer layers, while at the higher ionization parameters, Fe IV is the dominant species, with Fe III confined to a narrow zone at intermediate Lyman continuum optical depths. Also shown in Figure 3 is the hydrogen ionization balance as predicted by cloudy. The ionization of hydrogen has a strong influence on the iron ionization balance not only through the obvious effect of its opacity on the local ionizing radiation field, but also through strong charge-exchange reactions, $`\mathrm{Fe}^{2+}+\mathrm{H}\mathrm{Fe}^++\mathrm{H}^+`$.

Figure 4 shows the predicted ionization fractions in the same BLR models as Figure 3 except neglecting Fe-H charge-transfer reactions. There are significant differences, most notably for the higher ionization parameters at the lower particle density considered, models u20h96 and u13h96. The absence of charge exchange reactions in these models leads to Fe III, and not Fe II, as the dominant iron ionization stage throughout the cloud model. The influence of the charge-exchange reactions on the predicted iron fluxes will be discussed in a later section; such models may be partially supported by observation of I Zw 1 as they are better able to account for an intense feature near $`\lambda \mathrm{\hspace{0.17em}2418}`$Å identified by Laor et al. (1997) as Fe III UV47. We note that the charge exchange ionization and recombination rates for $`\mathrm{Fe}^{2+}+\mathrm{H}\mathrm{Fe}^++\mathrm{H}^+`$ predicted by Neufeld & Dalgarno (1987) were obtained in the rather crude Landau-Zener approximation (see Flower 1990, also Kingdon & Ferland 1996) and thus it is not inconceivable that this rate may be significantly overestimated.

Finally we note that the predicted fraction of Fe I in all cases, including the lowest ionization parameter considered, is very small, $`10^3`$. The question of weather or not such a small neutral fraction can account for the tentatively identified Fe I features in AGN spectra will be addressed in a later section.

## 5 Predicted Fe III Fluxes

Table 3 gives the total iron fluxes in the wavelength interval of 1500-5000 Å for the single-zone BLR cloud models of Table 2. The effect of increasing the cloud turbulent velocity from 10 to 50 $`\mathrm{km}\mathrm{s}^1`$ and of tripling the iron abundance relative to hydrogen are also given. The total flux is broken down into the percentage contributions of each of the iron ions. The higher density models with the largest ionization parameters, models u20h11 and u13h11, predict the largest contribution of the flux arising from Fe III, approaching $`35`$% in model u13h11. As also seen from the table, the percentage contribution of Fe III tends to decrease with increasing turbulent velocity but increases with an enhanced iron abundance; the u13h11 model with a threefold enhancement of the iron abundance comes close to having an equal flux split between Fe II and Fe III. It should be borne in mind that these comparisons are with the significantly smaller atom Fe II of Sigut & Pradhan (1998) to allow the simultaneous treatment of the current large Fe III atom.

Figures 5 and 6 show the wavelength distributions of the emitted iron line fluxes for the u20h96 and u13h11 models, respectively, in the case of a solar iron abundance and 10 km/s of turbulent velocity. The lower wavelength limit has been extended down from the 1500 Å of Table 3 to 500 Å in order to show the complete Fe III spectrum. The separate panels in these plots show the total flux as well as the contributions from the individual iron ions. There is a striking difference in the Fe III fluxes predicted for these two models. For u20h96, the strongest Fe III emission occurs at wavelengths near 1900 Å. These transitions represent decays from the lowest odd-parity Fe III level, $`\mathrm{z}^7\mathrm{P}^\mathrm{o}`$, to $`\mathrm{a}^7\mathrm{S}^\mathrm{e}`$ and are the only possible $`\mathrm{\Delta }S=0`$ transitions (see Figure 2). Inter-system ($`\mathrm{\Delta }S0`$) transitions from $`\mathrm{z}^7\mathrm{P}^\mathrm{o}`$ to $`\mathrm{a}^5\mathrm{D}^\mathrm{e}`$, the Fe III ground state, giving lines near $`\lambda \mathrm{\hspace{0.17em}1220}`$Å, and to $`\mathrm{a}^\mathrm{S}5^\mathrm{e}`$, giving lines near $`\lambda \mathrm{\hspace{0.17em}2400}`$Å (multiplet UV47), are also present but weak. The multiplet UV47 transitions identified as the strongest Fe III features in the spectrum of I Zw 1 by Laor et al. (1997) will be further discussed in Section 8.1. Decays from the only other low-lying odd-parity level, $`\mathrm{z}^5\mathrm{P}^\mathrm{o}`$, about $`0.8`$eV above $`\mathrm{z}^7\mathrm{P}^\mathrm{o}`$, are not prominent in the predicted spectrum for u20h96. A number of $`\mathrm{\Delta }S=0`$ decay paths from this level are possible, including the Fe III ground state $`\mathrm{a}^5\mathrm{D}^\mathrm{e}`$ resulting in transitions near $`\lambda \mathrm{\hspace{0.17em}1120}`$Å, as well as decays to $`\mathrm{a}^5\mathrm{P}^\mathrm{e}`$ and $`\mathrm{b}^5\mathrm{D}^\mathrm{e}`$ giving rise to transitions with $`\lambda >4000`$Å.

The higher ionization parameter, higher density model u13h11 presents a radically different picture. Now, the region around 1000 Å dominates the Fe III emission. This region is filled with decays from the numerous odd-parity levels near $`15`$eV to the lowest even parity levels ($`<5`$eV) of each spin system. The spectral region around $`\lambda \mathrm{\hspace{0.17em}2000}`$Å is now filled with strong Fe III emission arising principally from the decay of odd parity levels to the numerous even parity levels near $`10`$eV. As noted in Table 3, a significant fraction of the total iron flux comes out in the form of Fe III transitions.

In both of these models, fluorescent excitation by Lyman-$`\alpha `$ is unimportant due to the absence of significant line-overlap with any Fe III transitions. However, fluorescent excitation by Lyman-$`\beta `$ offers more potential. There are transitions from several low-lying, even-parity triplet states within a few Ångstroms of Lyman-$`\beta `$ (at $`1025.72`$Å): $`\mathrm{a}^3\mathrm{P}^\mathrm{e}`$-$`\mathrm{z}^5\mathrm{S}^\mathrm{o}`$ at $`1025.71`$Å ($`f_{ij}=0.0022`$), $`\mathrm{a}^3\mathrm{D}^\mathrm{e}`$-$`\mathrm{y}^3\mathrm{P}^\mathrm{o}`$ at $`1024.11`$Å ($`f_{ij}=0.023`$), and $`\mathrm{a}^3\mathrm{D}^\mathrm{e}`$-$`\mathrm{y}^3\mathrm{P}^\mathrm{o}`$ at $`1026.79`$, $`1026.88`$Å ($`f_{ij}=0.026`$,$`f_{ij}=0.0011`$). Re-performing the u20h96 model without Lyman-$`\beta `$ fluorescent excitation resulted in a reduction in the total predicted Fe III flux of only about 12%; the reduction in the case of u13h11 from the absence of Lyman-$`\beta `$ fluorescent excitation was negligible.

## 6 Predicted Fe I Fluxes

According to Table 3, $`0.15`$% of the iron line flux comes out in transitions of Fe I. In Figure 6, representing model u13h11, the predicted emission due to Fe I in the $`2700`$-$`4500`$Å region is, although quite weak, not entirely negligible. This is illustrated in Figure 7; the strongest transitions present in this figure are just those expected from the Fe I Grotrian diagram of Figure 1, namely decays from low-lying odd parity levels between $`3.5`$ to $`4.5`$eV to the three lowest even parity states, the $`\mathrm{a}^5\mathrm{D}^\mathrm{e}`$ ground state, $`\mathrm{a}^5\mathrm{F}^\mathrm{e}`$, and $`\mathrm{a}^3\mathrm{F}^\mathrm{e}`$. The strongest predicted multiplet is 23, $`\mathrm{z}^5\mathrm{G}^\mathrm{o}`$ \- $`\mathrm{a}^5\mathrm{F}^\mathrm{e}`$, giving lines near $`\lambda \mathrm{\hspace{0.17em}3600}`$Å. Such multiplets of Fe I have been tentatively identified in PHL 1092 by Bergeron & Knuth (1980) and Kwan et al. (1995) and in IRAS 07598+6508 by Kwan et al.  Both objects are Fe II-strong quasars. The current calculations tend to support these identifications and the occurrence of Fe I transitions in some AGN spectra. Comparing our predictions to the IRAS 07598+6508 spectrum given by Kwan et al. (their Figure 1), we note that the observed Fe I flux relative to the nearby Fe II flux between $`\lambda \mathrm{\hspace{0.17em}3600}3800`$Å is much larger than predicted by our models, even when considering the effect of the tentative Ti I blends identified by Kwan et al.  However, it is likely that our models underestimate the Fe I flux; the very limited Fe I atomic model employed artificially suppresses recombination by omitting energy levels within $`3`$eV of the continuum. This acts both to reduce the recombination contribution the line fluxes and to lower the Fe I ionization fraction which also leads to weaker lines. We plan to extend our treatment of the Fe I atomic model in the future to provide more realistic flux estimates to compare with these observations.

## 7 Dependence of Iron Flux on BLR Cloud Parameters

Figure 8 summarizes the dependence of the Fe I, Fe II, and Fe III fluxes on the various model parameters. The trend of increasing iron line fluxes with both ionization parameter and particle density is evident, as well as the maximization of the Fe III fluxes for the higher ionization parameter models. The figure also shows increased flux in Fe II and Fe III for either an increased cloud turbulent velocity or an increased iron abundance. Interestingly, however, the predicted Fe I flux does not depend on the internal cloud turbulent velocity. The lines of Fe I are weak and unsaturated so are insensitive to the turbulent velocity, contrary to the lines of the other ionization stages. From this figure it might seem that the internal cloud turbulent velocity might be deducible from the spectrum of an object exhibiting lines of both Fe I and Fe II. In fact, this is the classic technique in stellar atmospheres to determine the turbulent velocity dispersion (often dubbed microturbulence)–forcing weak and strong lines to be reproduced in strength by the same abundance. It might be worthwhile to attempt this in AGN; however, the situation is much more complex as the entire BLR spectrum is likely not formed within a single cloud or ensemble of identical clouds (Baldwin et al. 1995); the question would arise as to whether clouds of different ionization parameters and particle densities would have similar turbulent velocities. Nevertheless it would be still worthwhile to attempt this investigation. Knowledge of the intrinsic width of an emission line for individual clouds in the BLR has important implications on the uncertain nature of the BLR clouds themselves. Narrow, near thermal widths of $`10`$$`\mathrm{km}\mathrm{s}^1`$ imply large numbers of BLR clouds in oder to reproduce the smooth observed line profiles. Much larger intrinsic widths would require fewer clouds to produce smooth profiles, perhaps allowing stellar models for the BLR clouds to become viable (Peterson, Pogge & Wanders, 1999).

## 8 Comparison to I Zw 1 Fe II-Fe III UV Template

Vestergaard & Wilkes (2001) provide an empirical Fe II-Fe III UV template ($`1250\lambda \mathrm{\hspace{0.17em}3080}`$Å) derived from I Zw 1, the prototypical narrow-line Seyfert 1 galaxy. One of the main results of this work is the careful documentation of extensive Fe III emission from I Zw 1, building on the identifications of Laor et al. (1997). Figure 9 compares four of our basic models with a solar iron abundance and minimal turbulent velocity (10 $`\mathrm{km}\mathrm{s}^1`$) with the I Zw 1 template. In these comparisons, the calculations have been broadened by convolution with a Gaussian to $`500`$$`\mathrm{km}\mathrm{s}^1`$(FWHM), and the template has been normalized to the median model flux in the region of strongest Fe III flux included in the template wavelength region, $`1800`$-$`2000`$Å. Only our highest-density, highest ionization parameter model, u13h11, correctly reproduces the overall level of the Fe II UV emission line strength in the $`2000`$-$`2500`$Å region. Thus we find that there is, in principle, no problem in accounting for the basic level of Fe III flux from I Zw 1, even with a crude, single-zone model.

While inspection of Figure 9 clearly reveals a strong correlation between the model and the empirical template, there are several disagreements in detail. In the principle Fe III wavelength region, $`1800`$-$`2000`$Å, the relative strengths of individual features are not correctly reproduced. Most notable, however, is the failure of the basic model to correctly reproduce the strength of the Fe III feature near $`\lambda \mathrm{\hspace{0.17em}2418}`$Å identified by Laor et al. (1995). This strong feature is associated with Fe III multiplet UV47, the inter-system transition $`\mathrm{z}^7\mathrm{P}^\mathrm{o}`$-$`\mathrm{a}^5\mathrm{S}^\mathrm{e}`$, giving (identified) features near $`\lambda \mathrm{\hspace{0.17em}2418.58}`$Å ($`J=32`$, $`f_{ij}=0.0027`$) and $`\lambda \mathrm{\hspace{0.17em}2438.18}`$Å ($`J=22`$, $`f_{ij}=0.0011`$). Finally, some of the disagreement in the $`2200`$-$`2600`$Å region can be traced to our use of a limited Fe II model atom.

To assess the seriousness of these discrepancies, we have performed additional calculations along two fronts: first, we considered predictions of models in which the highly-uncertain Fe-H charge-exchange reactions were omitted. As noted in section 4, omission of these rates can substantially alter the iron ionization balance for certain models; secondly, we have performed a series of Monte Carlo simulations to assess the effect of errors in the basic atomic data have in the predicted Fe III line fluxes. We shall first discuss models omitted charge-exchange reactions.

### 8.1 Charge-Exchange and Fe III UV47 in I Zw 1

The important role played by the rather uncertain Fe-H charge exchange reaction rates motivated us to compute a series of models in which these rates were omitted. The predicted iron ionization balances have already been discussed in section 4 with the results shown in Figure 4. The predicted flux of model u20h96 without charge-exchange reactions is shown in Figure 10 and these predictions may be compared with Figure 5. As expected by the change in the ionization balance shown in Figures 3 and 4, the model without charge-exchange reactions has a lower Fe II flux and a higher Fe III flux. It is still the case that normalizing the predicted iron fluxes to the observed template in the interval $`1800`$-$`2000`$Å results in an over-prediction of the Fe II flux, but the over-prediction is now less. Now conspicuous in the model omitting charge-exchange reactions is the presence of strong lines of Fe III multiplet UV47. However, this is not translated into a significantly better fit to the Vestergaard & Wilkes template as shown by the upper panel of Figure 10. While there is now a strong predicted feature at $`\lambda \mathrm{\hspace{0.17em}2418}`$Å corresponding to Fe III, the fit to the Fe II flux retains strong features for $`\lambda 2400`$Å, centred around Fe II $`\lambda \mathrm{\hspace{0.17em}2382.04}`$Å ($`\mathrm{z}^6\mathrm{F}^\mathrm{o}`$-$`\mathrm{a}^6\mathrm{D}^\mathrm{e}`$), which have no counterpart (in terms of observed strength) in the empirical template. We have also explored the possibility that the Fe III $`\lambda \mathrm{\hspace{0.17em}2400}`$Å feature is strengthened by numerous, blended Fe II transitions by merging our predicted Fe III fluxes with the more extensive Fe II fluxes of SP03 based on an 829-level Fe II model atom. However, such merged line-lists did not help to improve the fit by much.

The overall fit to the Fe II flux distribution is generally better for the higher density models ($`\mathrm{log}N=11.6`$) but unfortunately these models do not predict strong Fe III features for multiplet UV47, either with or without charge-exchange reactions.

### 8.2 Monte Carlo Simulation Estimates for Error Bounds

Given that all of the models presented have shortcomings when compared to the Vestergaard & Wilkes template, it is important to assess the accuracy to which the iron line fluxes can be computed. This is a complex question, encompassing everything from the accuracy of the underlying BLR cloud model(s) and the approximations used in the numerical methods to the accuracy of the underlying atomic data. In this work, we will address only the latter issue by asking how accurately the iron line fluxes can be predicted given the sometimes large uncertainties in the basic atomic data. Even this question is not completely straightforward. For example, a single $`A_{ji}`$ value clearly affects the flux in the $`ji`$ radiative transition. However, through the coupling of the atomic level populations by the statistical equilibrium equations, and the global coupling of the emitting volume through the transfer of radiation, a single $`A_{ji}`$ value can possibly affect the fluxes in many lines. And similar observations can be made, in principle, for all of the remaining atomic data, including the photoionization and recombination rates and the collisional excitation and de-excitation rates. To address all of these interconnections in a consistent manner, we have turned to the Monte Carlo simulation technique used by Sigut (1996). Each fundamental atomic parameter ($`A_{ji}`$ value, effective collision strength, photoionization cross section, collisional ionization cross section, and charge-transfer reaction rate) is assigned an uncertainty. Given these uncertainties, a set of atomic data is randomly realized and used to solve the radiative transfer-statistical equilibrium problem for the iron fluxes. Then, a new set of atomic data is randomly realized and the predicted fluxes for this new data found. This sequence is repeated, and the width of the distribution of each line flux can be taken as a measure of its uncertainly.

Table 4 lists the basic atomic data uncertainties adopted for our collection of iron atomic data. The error assignment is kept deliberately simple; in principle, uncertainties could be assigned on a transition by transition basis. The R-matrix collision strengths (cbb rates) for low ionization states computed under the Iron Project are the most accurate ones available. However, they also display extremely complicated structure due to autoionizing resonances in the near-threshold region which dominates the rate coefficient at $`T10^4`$K. Therefore, uncertainties as large as a factor of $`23`$ for individual transitions can not be ruled out. The $`A_{ji}`$ values (rbb rates) are expected to be more accurate; however, significant errors may still be present, particularly in several calculations where relativistic effects have not been considered in an ab initio manner. Such calculations are now in progress, but the computations are about an order of magnitude more difficult that those in LS coupling (see Nahar, 2003).

The distribution of the random scalings is also kept very simple: a set of uniform random deviates, $`r`$, is computed with the ran2 algorithm of Press et al. (1986). Given the uncertainty assignment, $`p`$, from Table 4, a uniform set of deviates for the logarithm of the scaling is found from the linear relation $`l=a+r(ba)`$ where $`a=\mathrm{log}_{10}(p)`$ and $`b=\mathrm{log}_{10}(p)`$. The actual set of scalings, $`d`$, is taken as $`d=10^l`$. This procedure ensures the following common-sense property of the scalings: for example, if a charge-exchange rate is assumed accurate to within a factor of 100, then one might expect the a scaling from $`0.11.0`$ to be as likely as one from $`110`$. However, this procedure also results in the mean scaling being different from one, namely $`\overline{d}=(pp^1)/(ln\mathrm{\hspace{0.17em}10}(ba))`$; for example, the mean scaling for $`p=100`$ is about $`10.8`$. Thus the most probable value for the rate is not that adopted in the default (unscaled) atomic model. However, given the uncertainty (a factor of 100 in this case), we do not consider this deviation significant.

Figure 11 summarizes the uncertainties due to the atomic data uncertainties in Table 4 on the predicted Fe II-Fe III flux for model u13h11. This model BLR was chosen as it gave the best fit to the overall Fe II-Fe III flux levels in the Vestergaard & Wilkes template. Charge-exchange reactions between iron and hydrogen were included in this calculation, although assumed uncertain to within a factor of 10 (Table 4). The top panel displays the minimum and maximum predicted flux at each wavelength and the middle panel shows the standard deviation of each predicted flux expressed as a percentage of the average flux at that wavelength. The largest uncertainties approach 30% for the Fe III features near 2000 Å. This reflects mostly the scaling of the charge-exchange reaction rates which strongly affects the Fe II-Fe III ionization balance. Much smaller uncertainties are predicted for the bulk of the Fe II emission; this smaller uncertainty reflects, to some extent, the large number of individual Fe II transitions contributing to each wavelength in the broadened spectrum (see the bottom panel of Figure 11), each with an individually scaled uncertainty following Table 4. In light of this figure, the detailed discrepancies between the I Zw 1 template of Vestergaard & Wilkes and the computed model model shown in Figure 9 are mostly explained as uncertainties in the underlying atomic data.

## 9 Discussion

The primary aim of this continuing project is to incorporate advanced methods of non-LTE radiative transfer from stellar astrophysics into emission line analysis of AGN. In the current work, we have provided predictions for the entire low-ionization spectrum of iron, focusing on Fe III. Given the preliminary nature of the models, and the use of only single-zone BLR cloud models, a reasonable fit to the Fe II-Fe III UV template of Vestergaard & Wilkes for I Zw 1 was obtained.

An aspect in which our treatment is incomplete is in the detailed treatment of the Fe IV atom, the dominant ionization state in the fully ionized optically thin region of BLR models. The atomic data for Fe IV are in hand, computed from the Iron Project, and we plan to incorporate those into a further extension of the present work on par with Fe II and Fe III. The Fe IV lines may also contribute to UV spectra of BLR; for example, see observations of the Orion nebula using HST (Rubin et al. 1997).

In addition, the $`\lambda \mathrm{\hspace{0.17em}2418}`$Å feature in the computed and observed Fe III spectra offers a diagnostics of atomic processes and physical conditions and of uncertainties in atomic data. The combined role of these factors in determining the intensity of this emission feature may be illustrated by considering level-specific $`\mathrm{e}^{}+\mathrm{Fe}^{3+}\mathrm{Fe}^{2+}`$ recombination. An increase in the recombination-cascade rates might lead to: (a) less of a role for charge-exchange, as inferred by the non-charge-exchange models, (b) enhancement of the $`\lambda \mathrm{\hspace{0.17em}2418}`$ feature due to recombination contribution to the Fe III line, and (c) less flux in Fe II short-ward of $`2400`$Å, consistent with observations. The present recombination rates were derived from level-specific photoionization cross sections of Fe III computed in LS coupling (Nahar & Pradhan 1994). Later work has shown (Pradhan, Nahar & Zhang 2001, Nahar & Pradhan 2003) that near-threshold resonance complexes may affect the recombination rates by up to several factors. It is necessary to consider relativistic fine structure and resolve the resonance complexes in great detail (excited metastable state cross sections are affected much more than than the ground state). Nahar & Pradhan (2003) have developed a scheme for unified calculations of electron-ion recombination, including both the radiative and the dielectronic recombination processes. Such level-specific calculations for total recombination into the Fe III levels may possibly resolve the discrepancy in the observed and computed intensities of the $`\lambda \mathrm{\hspace{0.17em}2418}`$Å feature.

It is therefore not an exaggeration to say that in spite of the present scale of this effort, both in terms of atomic physics and radiative transfer, considerable work is still needed to improve the models.

We would like to thank Marianne Vestergaard for providing the UV I Zw 1 iron template. This work was supported by the Natural Sciences and Engineering Research Council of Canada (TAAS) and by the U.S. National Science Foundation and NASA (AKP & SN).