text:
# Acceleration of energetic particles by large-scale compressible magnetohydrodynamic turbulence

## 1 Introduction

Particle acceleration by magnetohydrodynamic (MHD) turbulence has been investigated by a number of authors (e.g., Fermi 1949, Hall & Sturrock 1967, Wentzel 1968, Kulsrud & Ferrari 1971, Eichler 1979, Achterberg 1981, Schlickeiser & Miller 1998), and has been studied in a number of astrophysical settings, including the interstellar medium (Ptuskin 1988, Schlickeiser & Miller 1998), solar flares (LaRosa et al 1994, Miller et al 1996, Miller et al 1997, Miller 1998, Schlickeiser & Miller 1998, Miller 1998), active galactic nuclei (Blackman 1999, Gruzinov & Quataert 1999), diffuse intracluster plasma in galaxy clusters (Eilek & Weatherall 1999, Weatherall & Eilek 1999, Dogiel 1999, Brunetti et al 1999, Ohno et al 2002, Fujita et al 2003), and the lobes of extragalactic radio sources (Eilek & Henriksen 1984, Manolakou et al 1999).

The way in which magnetohydrodynamic (MHD) turbulence accelerates energetic particles depends on the relative magnitudes of $`\lambda _{}`$ and $`\lambda _{\mathrm{mfp}}`$, where $`\lambda _{}`$ is the length of a turbulent eddy (or the wave length of a weakly damped wave) measured along the magnetic field, and $`\lambda _{\mathrm{mfp}}`$ is the energetic-particle scattering mean free path. In the weak-scattering limit, $`\lambda _{}\lambda _{\mathrm{mfp}}`$, a particle’s velocity $`v_{}`$ along the background magnetic field is approximately constant as the particle travels many wave lengths along the magnetic field (assuming that the amplitudes of the fluctuations are sufficiently small), and the particle interacts resonantly with a weakly-damped wave when the wave frequency in the frame of reference moving along the magnetic field at speed $`v_{}`$ is an integer multiple of the particle’s gyrofrequency. In the strong-scattering limit, $`\lambda _{}\lambda _{\mathrm{mfp}}`$, a particle travels diffusively across an eddy of length $`\lambda _{}`$, and acceleration is non-resonant.

Particle acceleration in the strong-scattering limit is the focus of this paper. Ptuskin (1988) studied particle acceleration in the strong-scattering limit by a power-law spectrum of acoustic or magnetosonic waves extending from a large scale $`l`$ to a smaller scale $`l_1`$, with rms fluctuation velocity $`u_{\mathrm{rms}}`$ dominated by the waves at scale $`l`$. He assumed that particle diffusion in space is isotropic with diffusion coefficient $`D`$, as opposed to primarily along field lines, and that some type of small-scale wave distinct from the acoustic/magnetosonic waves scatters the particles, with $`\lambda _{\mathrm{mfp}}D/vl_1`$, where $`v`$ is the particle speed. He found that for $`Dv_\mathrm{w}l`$, where $`v_\mathrm{w}`$ is the wave phase velocity, particles diffuse in momentum space according to the equation

$$\frac{f}{t}=\frac{1}{p^2}\frac{}{p}\left(p^2D_p\frac{f}{p}\right),$$
(1)
where $`f`$ is the particle distribution function averaged over the waves, and the momentum diffusion coefficient $`D_p`$ is given by

$$D_p=\frac{2p^2u_{\mathrm{rms}}^2}{9D},$$
(2)
where $`p`$ is the particle momentum. The condition $`Dv_\mathrm{w}l`$ implies that particles diffuse a distance $`l`$ in a time much shorter than the wave period at scale $`l`$. In this limit, the largest-scale waves make the dominant contribution to $`D_p`$.

This paper extends Ptuskin’s analysis to strong compressible MHD turbulence taking into account the reversible motion of particles along magnetic field lines. It is assumed that the rms turbulent velocity $`u_{\mathrm{rms}}`$ equals or exceeds the Alfvén speed $`v_\mathrm{A}`$, and that there is a single dominant length scale (outer scale) denoted $`l`$ for both velocity and magnetic fluctuations. Only the fluctuations at scale $`l`$ are considered. As in Ptuskin’s (1988) study, it is assumed that some type of small-scale wave scatters particles, and that $`\lambda _{\mathrm{mfp}}l`$. It is assumed that eddies are randomized in a time $`\tau _{\mathrm{rand}}l/u_{\mathrm{rms}}`$, and that $`D_{}u_{\mathrm{rms}}l`$, where $`D_{}`$ is the diffusion coefficient for particle motion along field lines. Thus, a particle diffuses a distance greater than $`l`$ during a time $`\tau _{\mathrm{rand}}`$, and can diffuse through and then return to an eddy multiple times before the eddy is randomized in the turbulent flow. It is shown in section 2 that this increases the coherence of stochastic particle acceleration, effectively increasing $`D_p`$ above the value in equation (2) by a factor of $`N`$, where $`N\sqrt{D_{}/ul}`$ is the typical number of times a particle diffuses through and then returns to an eddy before the eddy is randomized, assuming a particle is tied to a single field line. The discussion of section 2 consists of estimates based on physical arguments; these estimates are tested and validated by Monte Carlo particle simulations in section 3. In section 4 the results of the paper are applied to particle acceleration in the interstellar medium and the lobes of radio galaxies.

It is known that small cross-field motions and the divergence of neighboring magnetic field lines in turbulent plasmas allow a particle to escape from its initial field line by traveling sufficiently far along the magnetic field. This places an upper limit on $`N`$, which is explained in the Appendix. It is shown that this upper limit is not important in the astrophysical examples considered in section 4 given the assumed values of $`D_{}`$. However, for sufficiently large values of $`D_{}`$, the upper limit on $`N`$ becomes important.

## 2 Phenomenological description of momentum diffusion

Since the scattering mean free path is assumed to be $`l`$ and only those eddies at scale $`l`$ are considered, a particle’s momentum $`p`$ evolves according to the equation

$$\frac{dy}{dt}=\frac{𝐮}{3},$$
(3)
where

$$y=\mathrm{ln}(p/p_0),$$
(4)
$`p_0`$ is a particle’s initial momentum, and $`𝐮`$ is the turbulent velocity associated with the large-scale eddies (Skilling 1975, Ptuskin 1988). As a particle diffuses in space, it encounters a series of eddies that cause random increments in $`y`$, leading to diffusion in $`y`$-space with diffusion coefficient

$$D_y=\frac{(\delta y)_{\mathrm{rms}}^2}{2\delta t},$$
(5)
where $`\delta t`$ is the time during which increments to $`y`$ remain correlated, and $`(\delta y)_{\mathrm{rms}}`$ is the rms change in $`y`$ during a time $`\delta t`$. It is assumed that the turbulent velocities are randomized in a time

$$\tau _{\mathrm{rand}}=\frac{bl}{u_{\mathrm{rms}}},$$
(6)
where $`b`$ is a constant of order unity. The characteristic time for a particle to diffuse through an eddy is

$$\tau _{\mathrm{diff}}=\frac{l^2}{2D_{}}.$$
(7)
It is assumed that $`D_{}u_{\mathrm{rms}}l`$, and thus

$$\tau _{\mathrm{diff}}\tau _{\mathrm{rand}}.$$
(8)
In this section, it is assumed that a particle remains tied to the same magnetic field line. Thus, given equation (8), a particle typically diffuses through an eddy and later returns to the same eddy before the eddy is randomized. During such successive encounters, the values of $`𝐮`$ within the eddy are correlated, and the corresponding increments to $`y`$ are correlated. The correlation time $`\delta t`$ for random steps in $`y`$-space is thus $`\tau _{\mathrm{rand}}`$, and not $`\tau _{\mathrm{diff}}`$.

The increment to $`y`$ during a time $`\tau _{\mathrm{rand}}`$ is

$$\delta y=\frac{1}{3}_0^{\tau _{\mathrm{rand}}}𝐮𝑑t.$$
(9)
It is assumed that there is no average compression or decompression in the turbulent flow, so that

$$\delta y=0,$$
(10)
where $`\mathrm{}`$ indicates an average over an arbitrarily large number of random steps. During a time $`\tau _{\mathrm{rand}}`$, a particle’s position is distributed over an interval of length $`Nl`$, where

$$N=\frac{\sqrt{D_{}\tau _{\mathrm{rand}}}}{l}.$$
(11)
Since the values of $`𝐮`$ in different eddies are uncorrelated, the value of $`𝐮`$ averaged over a distance $`Nl`$ is $`N^{1/2}(𝐮)_{\mathrm{rms}}`$, where $`(𝐮)_{\mathrm{rms}}`$ is the rms value of $`𝐮`$. Thus, the rms value of $`\delta y`$ is approximately

$$(\delta y)_{\mathrm{rms}}\frac{N^{1/2}(𝐮)_{\mathrm{rms}}\tau _{\mathrm{rand}}}{3}.$$
(12)
Since only eddies of scale $`l`$ are considered, and since $`u_{\mathrm{rms}}v_\mathrm{A}`$,

$$(𝐮)_{\mathrm{rms}}=\frac{au_{\mathrm{rms}}}{l},$$
(13)
where $`a`$ is a constant of order unity. Thus, if $`N1`$, then $`(\delta y)_{\mathrm{rms}}1`$, and the probability $`P(y,t)dy`$ that $`\mathrm{ln}(p/p_0)(y,y+dy)`$ at time $`t`$ satisfies the Fokker-Planck equation

$$\frac{P}{t}=\frac{^2}{y^2}(D_yP),$$
(14)
where $`D_y`$, which is defined in equation (5), is given by

$$D_y\frac{Na^2u_{\mathrm{rms}}^2}{18D_{}}.$$
(15)
The log-momentum diffusion coefficient $`D_y`$ is analogous to Ptuskin’s (1988) momentum diffusion coefficient $`D_{p,\mathrm{Ptuskin}}`$ divided by $`p^2`$. If Ptuskin’s coefficient $`D`$ for isotropic spatial diffusion is set equal to $`D_{}`$, then in the large-$`N`$ limit

$$\frac{p^2D_y}{D_{p,\mathrm{Ptuskin}}}N.$$
(16)
This ratio can be understood as follows. Ptuskin’s result is approximately recovered by taking $`\delta t\tau _{\mathrm{diff}}`$. Such a relation is valid if spatial diffusion is isotropic, since there is only a small chance that a particle returns to the same eddy multiple times before the eddy is randomized when $`Dlu_{\mathrm{rms}}`$. However, when a particle diffuses along a single field line, the particle’s location during a time $`\tau _{\mathrm{rand}}=N^2\tau _{\mathrm{diff}}`$ is distributed among roughly $`N`$ eddies, and thus a particle spends a time $`N\tau _{\mathrm{diff}}`$ in an eddy before the eddy is randomized. The particle’s random walk in $`y`$-space can thus be thought of consisting of individual steps of duration $`N\tau _{\mathrm{diff}}`$ and magnitude $`\delta y=N\tau _{\mathrm{diff}}(𝐮)_{\mathrm{rms}}/3`$. The resulting value of $`D_y`$ is a factor of $`N`$ larger than in the isotropic-spatial-diffusion case in which a particle takes individual random steps in $`y`$-space of duration $`\tau _{\mathrm{diff}}`$ and magnitude $`\delta y=\tau _{\mathrm{diff}}(𝐮)_{\mathrm{rms}}/3`$. Thus, the reversible motion of a particle along a field line makes the particle’s random walk in $`y`$-space more coherent, and its diffusion in $`y`$-space more rapid, than in the isotropic-spatial-diffusion case.

One consequence of equation (15) is that $`D_y\mathrm{}`$ as $`\tau _{\mathrm{rand}}\mathrm{}`$. The reason is that if $`𝐮`$ were constant in time along a field line, increments to $`y`$ would remain correlated forever if a particle were tied to that field line, and the particle would undergo superdiffusion in $`y`$-space. This phenomenon can be described approximately as follows. During a time $`t`$, a particle’s position is distributed over $`M`$ eddies, where

$$M=\frac{\sqrt{D_{}t}}{l}.$$
(17)
Since the values of $`𝐮`$ within different eddies are uncorrelated, the typical value of $`𝐮`$ averaged over $`M`$ eddies is $`M^{1/2}(𝐮)_{\mathrm{rms}}`$. For large $`M`$, the rms value of $`y`$ for a particle obeys the equation

$$\frac{dy_{\mathrm{rms}}}{dt}\frac{M^{1/2}(𝐮)_{\mathrm{rms}}}{3}t^{1/4}.$$
(18)
Integrating and squaring the result, one obtains

$$y_{\mathrm{rms}}^2\frac{16l(𝐮)_{\mathrm{rms}}^2t^{3/2}}{81\sqrt{D_{}}}.$$
(19)
## 3 Numerical simulations

In this section, particle acceleration by compressible turbulence is investigated with the use of Monte Carlo particle simulations. Particles are assumed to follow a single field line. Distance along the field line is denoted by the coordinate $`x`$. The values of $`𝐮`$ at space-time grid points $`(x,t)=(ml,n\tau _{\mathrm{rand}})`$ are set to either $`+\tau _{\mathrm{rand}}^1`$ or $`\tau _{\mathrm{rand}}^1`$ with equal probability, where $`m`$ and $`n`$ are integers. The value of $`𝐮`$ between grid points is obtained by linear interpolation in space and time, which gives $`(𝐮)_{\mathrm{rms}}\tau _{\mathrm{rand}}=2/3`$, or, equivalently, $`ab=2/3`$, where $`a`$ and $`b`$ are defined in equations (6) and (13). During each time step of duration $`\delta t`$, a particle’s $`y`$-value is incremented by an amount $`𝐮\delta t/3`$, and there is an equal probability that its $`x`$-coordinate is either increased or decreased by an amount $`\sqrt{2D_{}\delta t}`$.

In figure 1, $`D_y`$ is plotted for different values of $`D_{}`$. To express the results in terms of $`u_{\mathrm{rms}}`$ and $`l`$, it is assumed that $`a=1`$. The values of $`D_y`$ in the numerical simulations, depicted with solid triangles in figure 1, are obtained by dividing the average of $`y^2`$ for $`10^4`$ simulated particles by $`2t`$ at $`t=130\tau _{\mathrm{rand}}`$. The solid line represents the value of $`D_y`$ in equation (15). The dashed line represents the scaling $`D_yD_{}^1`$, analogous to Ptuskin’s (1988) scaling $`D_pD^1`$ for isotropic spatial diffusion of particles.

The increase of $`y_{\mathrm{rms}}^2`$ with time when $`𝐮`$ is constant in time along a field line is shown in figure 2. The solid line is the value of $`y_{\mathrm{rms}}^2`$ in equation (19). The averages of $`y^2`$ for $`10^4`$ simulated particles are given by the solid triangles and show superdiffusive behavior in approximate agreement with equation (19).

## 4 Application to Galactic cosmic rays and relativistic electrons in radio galaxies

From equations (14) and (15) the characteristic time for particle acceleration by large-scale compressible turbulence is

$$\tau _{\mathrm{acc}}=D_y^1\frac{18D_{}^{1/2}l^{1/2}}{u_{\mathrm{rms}}^{3/2}},$$
(20)
where it is assumed that $`a=1`$ in equation (13) and $`b=1`$ in equation (6). The values of $`D_{}`$ for Galactic cosmic rays and relativistic electrons in the lobes of radio galaxies are not well known. There is disagreement as to the correct way to model the waves that cause energetic particle scattering (see, e.g., Bieber et al 1994, Schlickeiser & Miller 1998, Chandran 2000, Felice & Kulsrud 2001, Yan & Lazarian 2002). Some studies indicate that $`D`$ and $`D_{}`$ have different values in different phases of the ISM (e.g., Felice & Kulsrud 2001), and the values of $`D`$ and $`D_{}`$ may be different in the disk and halo of the Galaxy (Berezinskii et al 1990). In the hot coronal plasma of the ISM, $`u_{\mathrm{rms}}10^7`$ cm/s (Shelton et al 2001). This value is comparable to the sound speed, and thus the turbulence is compressible. It is assumed that $`l3\times 10^{20}`$ cm. If $`D_{}=10^{28}\text{ cm}^2/\text{s}`$, then $`\tau _{\mathrm{acc}}=3\times 10^7`$ yr. If $`D_{}=10^{29}\text{ cm}^2/\text{s}`$, then $`\tau _{\mathrm{acc}}=10^8`$ yr. These values of $`\tau _{\mathrm{acc}}`$ are comparable to the Galactic confinement time of $`1`$-GeV cosmic rays, and are the types of time scales required if reacceleration is to explain the peaking of secondary/primary ratios at a few GeV (Ptuskin 2001). For relativistic electrons in radio lobes, it is simply assumed that $`D_{}=v_\mathrm{A}l_{\mathrm{lobe}}`$, where $`l_{\mathrm{lobe}}=100`$ kpc is the typical size of a radio lobe. For this value of $`D_{}`$, the time for electrons to diffuse out of a radio lobe is of order the Alfvén crossing time of the lobe. The thermal plasma density is taken to be $`5\times 10^5\text{ cm}^3`$, which is an an upper limit based upon Faraday depolarization, and $`B`$ is taken to be 20 $`\mu `$G, a typical equipartition field (Spangler & Pogge 1984, Spangler & Sakurai 1985). These values yield $`v_\mathrm{A}=6\times 10^8`$ cm/s, and it is assumed that $`u_{\mathrm{rms}}=v_\mathrm{A}`$. The value of $`l`$ for radio-lobe turbulence is taken to be $`l_{\mathrm{lobe}}/30`$ based on the standard deviation of the linear Stokes parameters in 3C 166 (Spangler 1983). For these parameters, $`\tau _{\mathrm{acc}}5\times 10^7`$ yr. This type of time scale is comparable to the ages of relativistic-electron populations obtained from synchrotron aging estimates, which suggests that reacceleration may have to be taken into account in models of radio-lobe electrons.

## 5 Summary

It is shown that the reversible motion of energetic particles along the magnetic field enhances the rate of particle acceleration by large-scale compressible turbulence relative to the isotropic-spatial-diffusion case. When this effect is taken into account, the time scale for acceleration of relativistic electrons by large-scale compressible turbulence in the lobes of radio galaxies is $`5\times 10^7`$ yr, and the time scale for acceleration of 1-GeV protons by large-scale compressible turbulence in the coronal plasma of the ISM is $`310\times 10^7`$ yr. The main source of uncertainty in these estimates is the uncertainty in $`D_{}`$; other sources include the uncertainty in $`l`$, $`u_{\mathrm{rms}}`$, and the dimensionless constants $`a`$ and $`b`$ that relate $`\tau _{\mathrm{rand}}`$ and $`(𝐮)_{\mathrm{rms}}`$ to $`u_{\mathrm{rms}}`$ and $`l`$.

I thank Eliot Quataert, Steve Spangler, Torsten Ensslin, and Eric Blackman for helpful comments. This work was supported by NSF grant AST-0098086 and DOE grants DE-FG02-01ER54658 and DE-FC02-01ER54651.

## Appendix A Particle escape from field lines through cross-field motion and field-line divergence

Particles with gyroradii $`\rho _g`$ much less than the dominant magnetic field scale length $`l`$ drift slowly across magnetic field lines because of gradients in the magnetic-field strength, curvature of field lines, and wave pitch-angle scattering. Cross-field motions enhanced by the divergence of neighboring field lines allow particles to “escape” from their initial field lines, that is, to separate from their initial field lines by a distance $`l`$ (Rechester & Rosenbluth 1978, Chuvilgin & Ptuskin 1996, Chandran & Cowley 1998). This can be seen with the aid of figure 3.

Suppose a particle starts out at point P on field line F1, and that its diffusive motion along the field initially moves it towards point Q. After moving a short distance, field gradients and collisions cause the particle to take a step of length $`\rho _g`$ across the magnetic field, from F1 to a new field line F2. Although the particle continuously drifts across the field to new field lines, let us assume for the moment that it remains attached to F2. As the particle diffuses along F2, F2 diverges from F1. Let $`z_\mathrm{s}`$ be the distance that F2 must be followed before F2 separates from F1 by a distance $`l`$. (Because the particle continuously drifts across the field, it typically separates from F1 after traveling a distance somewhat less than $`z_\mathrm{s}`$ along the field; this effect is ignored in this discussion.) After the particle moves a distance $`z_\mathrm{s}`$ along F2, its subsequent motion is not correlated with $`F_1`$. The particle proceeds to point R, and then its diffusive motion along the magnetic field changes direction, bringing it back towards point Q. Instead of following F2 back to point Q, however, the particle drifts across the field and ends up on a new field line F3. After following F3 for a distance $`z_\mathrm{s}`$, the particle separates from F2 by a distance $`l`$ and proceeds to point S. In this example, after the particle moves significantly farther than $`z_\mathrm{s}`$ along the magnetic field, it has a negligible probability of returning to within a distance $`l`$ of its initial location. The maximum time a particle can remain correlated with an eddy is thus of order

$$\tau _{\mathrm{esc}}=\frac{z_\mathrm{s}^2}{D_{}},$$
(A1)
where $`z_\mathrm{s}`$ is the average value of $`z_\mathrm{s}`$ over a large sample of field-line pairs.

Several studies have calculated $`z_\mathrm{s}`$. Jokipii (1973) and Skilling, McIvor, & Holmes (1974) found that $`z_\mathrm{s}l`$ for isotropic MHD turbulence with $`\delta BB_0`$, where $`\delta B`$ is the fluctuating magnetic field and $`B_0`$ is the strength of any background field in the system that is coherent over distances $`l`$. Narayan & Medvedev (2001) found $`z_\mathrm{s}l`$ for locally anisotropic MHD turbulence with $`\delta BB_0`$. A related quantity, the Kolmogorov entropy, which describes the exponential divergence of neighboring field lines while their separation is smaller than the dissipation scale of the turbulence, has been calculated by several authors (e.g., Rechester & Rosenbluth 1979, Zimbardo et al 1995, Barghouty & Jokipii 1996, Casse, Lemoine, & Pelletier 2001). Another related quantity, the diffusion coefficient for the wandering of single magnetic field line, has been studied by a number of authors (e.g., Matthaeus et al 1995, Michalek & Ostrowski 1998).

Maron, Chandran, & Blackman (2003) and Chandran & Maron (2003) calculated $`z_\mathrm{s}`$ by tracking field lines in direct numerical simulations of incompressible MHD turbulence with an inertial range extending from a dominant length scale $`l`$ to a much smaller dissipation scale $`l_d`$. In simulations with $`B_0=0`$, they found that $`z_\mathrm{s}`$ asymptotes to a value of order several $`l`$ as $`\rho _g`$ is decreased towards the dissipation scale $`l_d`$ in the large-$`l/l_d`$ limit. From figures 5 and 7 of Chandran & Maron (2003), this value is $`57l`$. We assume that $`z_s`$ has a similar value in turbulence with a mean magnetic field provided $`\delta BB_0`$.

The escape of particles from field lines modifies the estimate of $`D_y`$ presented in section 2 when $`\sqrt{D_{}\tau _{\mathrm{rand}}}>z_\mathrm{s}`$. In this case, a particle remains correlated with a magnetic field line for a time $`\tau _{\mathrm{esc}}`$ that is less than $`\tau _{\mathrm{rand}}`$. During the time $`\tau _{\mathrm{esc}}`$ the particle travels a distance $`z_\mathrm{s}`$, and its location is spread out over roughly

$$N^{}=\frac{z_\mathrm{s}}{l}$$
(A2)
eddies. The arguments leading to equations (14) and (15) can be repeated, replacing $`\tau _{\mathrm{rand}}`$ with $`\tau _{\mathrm{esc}}`$, to obtain

$$D_y\frac{N^{}a^2u_{\mathrm{rms}}^2}{18D_{}}.$$
(A3)
To summarize, if $`\tau _{\mathrm{esc}}<\tau _{\mathrm{rand}}`$, then equation (A3) is valid. If $`\tau _{\mathrm{esc}}>\tau _{\mathrm{rand}}`$, then equation (15) is valid.

For the astrophysical examples considered in section 4, $`\tau _{\mathrm{esc}}\tau _{\mathrm{rand}}`$, and thus equation (15) can be applied. However, for larger values of $`D_{}`$, $`\tau _{\mathrm{esc}}`$ can be reduced below $`\tau _{\mathrm{rand}}`$, in which case cross-field motion becomes important and equation (A3) applies.


format: text

text:
# Velocity and hierarchical spread of epidemic outbreaks in scale-free networks

## Abstract

We study the effect of the connectivity pattern of complex networks on the propagation dynamics of epidemics. The growth time scale of outbreaks is inversely proportional to the network degree fluctuations, signaling that epidemics spread almost instantaneously in networks with scale-free degree distributions. This feature is associated to an epidemic propagation that follows a precise hierarchical dynamics. Once the highly connected hubs are reached, the infection pervades the network in a progressive cascade across smaller degree classes. The present results are relevant for the development of adaptive containment strategies.

The connectivity pattern of the network of individual contacts has been acknowledged since long as a relevant factor in determining the properties of epidemic spreading phenomena Hethcote:1984 ; May:1992 . This issue assumes a particular relevance in the case of networks characterized by complex topologies and very heterogeneous structures barabasi:2000 ; mdbook ; psvbook ; lil ; may01  that in many cases present us with new epidemic propagation scenarios Moore:2000 ; Abramson:2001 ; Pastor:2001 ; lloyd01 ; Moreno:2002 ; newman02 . A striking example of this situation is provided by scale-free networks characterized by large fluctuations in the number of connections (degree) $`k`$ of each vertex. This feature usually finds its signature in a heavy-tailed degree distribution with power-law behavior of the form $`P(k)k^\gamma `$, with $`2\gamma 3`$, that implies a non-vanishing probability of finding vertices with very large degrees barabasi:2000 ; mdbook ; psvbook ; Amaral:2000 . The latter are the “hubs” or “superspreaders” Hethcote:1984  responsible for the proliferation of infected individuals whatever the rate of infection characterizing the epidemic, eventually leading to the absence of any epidemic threshold below which the infection cannot initiate a major outbreak Pastor:2001 . This new scenario is of practical interest in computer virus diffusion and the spreading of diseases in heterogeneous populations mdbook ; psvbook ; lil ; may01 . It also raises new questions on how to protect the network and find optimal strategies for the deployment of immunization resources Pastor:2002 ; dezso02 ; havlin03 . So far, however, studies of epidemic models in complex networks have been focused on the stationary properties of endemic states or the final prevalence (number of infected individuals) of epidemics. The dynamical evolution of the outbreaks has been instead far less investigated and a detailed inspection of the temporal behavior in complex topologies is still missing.

In this paper we intend to fill this gap by providing a first analysis of the time evolution of epidemic outbreaks in complex networks with highly heterogeneous connectivity patterns. We consider the time behavior of epidemic outbreaks and find that the growth of infected individuals is governed by a time scale $`\tau `$ proportional to the ratio between the first and second moment of the network’s degree distribution, $`\tau k/k^2`$. This implies an instantaneous rise of the prevalence in very heterogeneous networks where $`k^2\mathrm{}`$ in the infinite size limit. In particular, this result shows that scale-free networks with $`2\gamma 3`$ exhibit, along with the lack of an intrinsic epidemic threshold, a virtually infinite propagation velocity of the infection. Furthermore, we study the detailed propagation in time of the infection through the different degree classes in the population. We find a striking hierarchical dynamics in which the infection propagates via a cascade that progresses from higher to lower degree classes. This infection hierarchy might be used to develop dynamical ad-hoc strategies for network protection.

In order to study the dynamical evolution of epidemic outbreaks we shall focus on the susceptible-infected (SI) model in which individuals can be in two discrete states, either susceptible or infected May:1992 ; Murray:1993 . Each individual is represented by a vertex of the network and the edges are the connections between individuals along which the infection may spread. The total population (the size of the network) $`N`$ is assumed to be constant and if $`S(t)`$ and $`I(t)`$ are the number of susceptible and infected individuals at time $`t`$, respectively, then $`N=S(t)+I(t)`$. In the SI model the infection transmission is defined by the spreading rate $`\lambda `$ at which susceptible individuals acquire the infection from an infected neighbor. In this model infected individuals remain always infective, an approximation that is useful to describe early epidemic stages in which no control measures are deployed. We shall see in the following that the results obtained for the SI model can be readily generalized to more complex schemes.

A first analytical description of the SI model can be undertaken within the homogeneous mixing hypothesis May:1992 ; Murray:1993 , consisting in a mean-field description of the system in which all vertices are considered as being equivalent. In this case the reaction rate equation for the average density of infected individuals $`i(t)=I(t)/N`$ (the prevalence) reads as

$$\frac{di(t)}{dt}=\lambda ki(t)\left[1i(t)\right].$$
(1)
The above equation states that the growth rate of infected individuals is proportional to the spreading rate $`\lambda `$, the density of susceptible vertices that may become infected, $`s(t)=1i(t)`$, and the number of infected individuals in contact with any susceptible vertex. The homogeneous mixing hypothesis considers that this last term is simply the product of the number of neighbors $`k`$ and the average density $`i(t)`$. Obviously, this approximation neglects correlations among individuals and considers that all vertices have the same number of neighbors $`k`$, i.e. it assumes a perfectly homogeneous network. At small times, when the density of infected vertices is very small, we can neglect terms of order $`𝒪(i^2)`$ and obtain the leading behavior $`i(t)i_0e^{t/\tau _H}`$, where $`i_0`$ is the initial density of infected individuals and $`\tau _H=(\lambda k)^1`$ is the time scale governing the growth of the infection in a homogeneous network.

The above calculations recover that the outbreak’s time scale $`\tau _H`$ is inversely proportional to the epidemic reproduction rate; i.e. the mean number of infections generated by each infected individual. However, in populations with a heterogeneous connectivity pattern, it is known that the reproduction rate depends upon the contacts’ fluctuations May:1992 , and it is natural to expect that also the outbreaks’ time scale is analogously affected. Indeed, in heterogeneous networks the degree $`k`$ of vertices is highly fluctuating and the average degree is not anymore a meaningful characterization of the network properties. In order to take fully into account the degree heterogeneity in the dynamical evolution, it is possible to write down the whole set of reaction rate equations for the average densities of infected vertices of degree $`k`$, $`i_k(t)=I_k(t)/N_k`$, where $`N_k`$ and $`I_k(t)`$ are the number of vertices and infected vertices within each degree class $`k`$, respectively Pastor:2001 . For the SI model the evolution equations read as

$$\frac{di_k(t)}{dt}=\lambda k\left[1i_k(t)\right]\mathrm{\Theta }_k(t),$$
(2)
where the creation term is proportional to the spreading rate $`\lambda `$, the degree $`k`$, the probability $`1i_k`$ that a vertex with degree $`k`$ is not infected, and the density $`\mathrm{\Theta }_k`$ of infected neighbors of a vertex of degree $`k`$. The latter term is thus the average probability that any given neighbor of a vertex of degree $`k`$ is infected. In the case of uncorrelated networks corr , $`\mathrm{\Theta }_k\mathrm{\Theta }`$ is independent of the degree of the emanating edge. The probability that each edge of a susceptible is pointing to an infected vertex of degree $`k^{}`$ is proportional to the fraction of edges emanated from these vertices. By considering that at least one of the edge of each infected vertex is pointing to another infected vertex, from which the infection has been transmitted, one obtains

$$\mathrm{\Theta }(t)=\frac{_k^{}(k^{}1)P(k^{})i_k^{}(t)}{k},$$
(3)
where $`k=_k^{}k^{}P(k^{})`$ is the proper normalization factor dictated by the total number of edges. A reaction rate equation for $`\mathrm{\Theta }(t)`$ can be obtained from Eqs. (2) and (3). In the initial epidemic stages, we neglect terms of order $`𝒪(i^2)`$, obtaining the following set of equations

$`{\displaystyle \frac{di_k(t)}{dt}}`$ $`=`$ $`\lambda k\mathrm{\Theta }(t),`$ (4)
$`{\displaystyle \frac{d\mathrm{\Theta }(t)}{dt}}`$ $`=`$ $`\lambda \left({\displaystyle \frac{k^2}{k}}1\right)\mathrm{\Theta }(t).`$ (5)
These equations can be solved with the uniform initial condition $`i_k(t=0)=i_0`$ yielding for the total average prevalence $`i(t)=_kP(k)i_k(t)`$

$$i(t)=i_0\left[1+\frac{k^2k}{k^2k}(e^{t/\tau }1)\right],$$
(6)
where

$$\tau =\frac{k}{\lambda (k^2k)}.$$
(7)
This readily implies that the growth time scale of an epidemic outbreak is related to the graph heterogeneity. Indeed, the ratio $`\kappa =k^2/k`$ is the parameter defining the level of heterogeneity of the network, since the normalized degree variance can be expressed as $`\kappa /k1`$ and therefore high levels of fluctuations correspond to $`\kappa k`$. In networks with a Poisson degree distribution in which $`\kappa =k+1`$, we recover the result $`\tau (\lambda k)^1`$. Instead, in networks with very heterogeneous connectivity patterns $`\kappa `$ is very large and the outbreak time-scale $`\tau `$ is very small, signaling a very fast diffusion of the infection. In particular, in scale-free networks characterized by a degree exponent $`2\gamma 3`$ we have that $`\kappa k^2\mathrm{}`$ with the network size $`N\mathrm{}`$. Therefore, while $`\tau `$ is a function of the finite size $`N`$, for large scale-free networks we face a virtually instantaneous rise of the epidemic incidence.

It is worth stressing that these results can be easily extended to the susceptible-infected-susceptible (SIS) and the susceptible-infected-removed models (SIR) May:1992 , in which the Eq. (2) contains an extra term $`\mu i_k(t)`$ defining the rate at which infected individuals of degree $`k`$ recover and become again susceptible or removed from the population, respectively. In addition, in the SIR model the normalization imposes that $`s_k(t)=1i_k(t)r_k(t)`$, where $`r_k(t)`$ is the density of removed individuals of degree $`k`$. The inclusion of the decaying term $`\mu i_k`$, does not change the picture obtained in the SI model. By using the same approximations, the time scale is found to behave as $`\tau k/(\lambda k^2(\mu +1)k)`$. In the case of diverging fluctuations the time-scale behavior is therefore still dominated by $`k^2`$ and $`\tau `$ is always positive whatever the spreading rate $`\lambda `$. This allows to recover the absence of an epidemic threshold, i.e. the lack of a decreasing prevalence region in the parameter space.

In order to check these analytical results we have performed numerical simulations of the SI model on networks generated with the Barabási-Albert (BA) algorithm barabasi99 . These networks have a scale-free degree distribution $`P(k)k^3`$ and we use different network sizes $`N`$ and minimum degree values $`m`$ in order to change the level of heterogeneity, that in this case is given by $`\kappa m\mathrm{ln}N`$. Simulations use an agent-based modeling strategy in which at each time step the SI dynamics is applied to each vertex by considering the actual state of the vertex and its neighbors. It is then possible to measure the evolution of the average number of infected individuals and other quantities, by averaging over $`10^3`$ realizations of the dynamics. In addition, given the stochastic nature of the model, different initial conditions and networks realizations can be used to obtain averaged quantities. In Fig. 1 we report the early time behavior of outbreaks in networks with different heterogeneity levels and the behavior of the measured $`\tau `$ with respect to $`k/\lambda (k^2k)`$. The numerical results recover the analytical prediction with great accuracy. Indeed, the BA network is a good example of heterogeneous network in which the approximations used in the calculations are satisfied corr . In networks with correlations we expect to find different quantitative results but a qualitatively similar framework as it happens in the case of the epidemic threshold evaluation corr .

The previous results show that the heterogeneity of scale-free connectivity patterns favors epidemic spreading not only by suppressing the epidemic threshold, but also by accelerating the virus propagation in the population. The velocity of the spreading leaves us with very short response times in the deployment of control measures and a detailed knowledge of the way epidemics propagate through the network could be very valuable in the definition of adaptive strategies. Indeed, the epidemic diffusion is far from homogeneous. The simple formal integration of Eq. (2) written for $`s_k=1i_k`$ yields $`s_k(t)=s_k^0e^{\lambda k\mathrm{\Phi }(t)}`$ where $`\mathrm{\Phi }(t)=_0^t𝑑t^{}\mathrm{\Theta }(t^{})`$. This result is valid for any value of the degree $`k`$ and the function $`\mathrm{\Phi }(t)`$ is positive and monotonically increasing. This last fact implies that $`s_k`$ is decreasing monotonically towards zero when time grows. For any two values $`k>k^{}`$, and whatever the initial conditions $`s_k^0`$ and $`s_k^{}^0`$ are, there exists a time $`t^{}`$ after which $`s_k(t)<s_k^{}(t)`$. This means that vertices belonging to higher degree classes are generally infected more quickly.

A more precise characterization of the epidemic diffusion through the network can be achieved by studying some convenient quantities in numerical spreading experiments in BA networks. First, we measure the average degree of the newly infected nodes at time $`t`$, defined as

$$\overline{k}_{\mathrm{inf}}(t)=\frac{_kk[I_k(t)I_k(t1)]}{I(t)I(t1)}.$$
(8)
In Fig. 2a) we plot this quantity for BA networks as a function of the rescaled time $`t/\tau `$. The curves show an initial plateau that can be easily understood by considering that at very low density of infected individuals $`i`$, each vertex will infect a fraction of its neighbors without correlations with the spreading from other vertices. In this case each edge points to a vertex with degree $`k`$ with probability $`kP(k)/k`$ and the average degree of newly infected vertices is given by $`\overline{k}_{\mathrm{inf}}(t)=k^2/k`$. After this initial regime, $`\overline{k}_{\mathrm{inf}}(t)`$ decreases smoothly when time increases. The dynamical spreading process is therefore clear; after the hubs are very quickly infected the spread is going always towards smaller values of $`k`$. This is confirmed by the large time regime that settles in a plateau $`\overline{k}_{\mathrm{inf}}(t)=m`$; i.e. the vertices with the lowest degree are the last to be infected.

Further information on the infection propagation is provided by the inverse participation ratio $`Y_2(t)`$ Derrida:1987 . We first define the weight of infected individuals in each class degree $`k`$ by $`w_k(t)=I_k(t)/I(t)`$. The quantity $`Y_2`$ is then defined as

$$Y_2(t)=\underset{k}{}w_k^2(t).$$
(9)
If $`Y_21/k_{\mathrm{max}}`$, infected vertices are homogeneously distributed among all degree classes. In contrast, if $`Y_2`$ is not small (of order unity) then the infection is localized on some specific degree classes that dominate the sum of Eq. (9). In Fig. 2b) we report the behavior of $`Y_2`$ versus time for BA networks with different minimum degree. The function $`Y_2`$ has a maximum at the early time stage, indicating that the infection is localized on the large $`k`$ classes, as we infer from the plot of $`k^{\mathrm{inf}}(t)`$, see Fig. 2a). Afterwards $`Y_2`$ decreases, with the infection progressively invading lower degree classes, and providing a more homogeneous diffusion of infected vertices in the various $`k`$ classes. Finally, the last stage of the process corresponds to the capillary invasion of the lowest degree classes which have a larger number of vertices and thus provide a larger weight. In the very large time limit, when the whole network is infected, $`Y_2(t=\mathrm{})=_kP(k)^2`$. Noticeably, curves for different levels of heterogeneity have the same time profile in the rescaled variable $`t/\tau `$. This implies that, despite the various approximations used in the calculations, the whole spreading process is dominated by the time-scale defined in the early exponential regime of the outbreak.

The presented results provide a clear picture of the infection propagation in heterogeneous networks. First the infection takes control of the large degree vertices in the network. Then it rapidly invades the network via a cascade through progressively smaller degree classes. The dynamical structure of the spreading is therefore characterized by a hierarchical cascade from hubs to intermediate $`k`$ and finally to small $`k`$ classes. This result, along with the very fast growth rate of epidemic outbreaks, could be of practical importance in the set-up of dynamic control strategies in populations with heterogeneous connectivity patterns. In particular, targeted immunization strategies that evolve with time might be particularly effective in epidemics control.

###### Acknowledgements.

A.B, A. V, and R.P.-S. . are partially funded by the EC-Fet Open project COSIN IST-2001-33555. R.P.-S. acknowledges support from the MCyT (Spain), and from the DURSI, Generalitat de Catalunya (Spain).


text:
# Proper motion surveys of the young open clusters Alpha Persei and the Pleiades

## 1 Introduction

The young open cluster Alpha Per lies at a distance of 183 pc ($`(Mm)_o=6.31`$, van Leeuwen 1999) with an age of 90 Myr (Stauffer, Barrado & Bouvier 1999). We assume in this paper that $`A_R=0.23`$(O’Dell, Hendry & Collier Cameron 1994). It has been reasonably well studied in the past using various techniques. Prosser’s catalogue (Prosser 1993) contains a comprehensive list of candidate member stars identified by photometry, spectroscopy and radial velocity measurements. The faintest of these has $`V21.8`$($`M_V15.5`$). These candidate stars have been supplemented by those identified in similar studies by Zapatero Osorio et al. (1996), Stauffer et al. (1999) and Rebolo et al. (1992). In addition Prosser, Randich & Simon (1998) identified candidate stars which were optical counterparts to X-ray sources within the cluster.
An early proper motion study by Heckmann, Dieckvoss & Kox (1956) identified stars with proper motions consistent with cluster membership. However this study only went down to $`V12`$ and did not calculate formal membership probabilities. Prosser (1992) also used stellar proper motions as an initial selection technique. Although this study was deeper than the earlier one ($`V=18.8`$) it also did not calculate formal membership probabilities.
Most recently Barrado y Navascués et al (2001) produced a deep CCD survey of the cluster. This identified several candidate low mass stars and Brown Dwarfs down to $`I_c=22`$. They then used this data to produce a mass function. The mass function was fitted with a power law whose slope was found to be $`\alpha =0.59`$.
The Pleiades is the best studied open cluster, here we assume $`(Mm)_o=5.4`$ (van Leeuwen 1999) and $`A_R=0.15`$ (O’Dell, Hendry & Collier Cameron 1994). In recent years several deep CCD surveys of the cluster such as Dobbie et al (2002) and Moraux et al (2003) have produced candidates down to $`0.03M_{}`$. These have been complimented by proper motion studies such as those of Hambly et al (1999) and Moraux, Bouvier and Stauffer (2002). However the only wide field, high precision proper motion survey is that of Hambly, Hawkins and Jameson (1991). Adams et al (2001) used 2MASS data along with proper motions to produce a wide field study. This contained many candidate stars outside the tidal radius of the cluster showing background star contamination is a problem for many techniques.

## 2 Observational data and reduction

### 2.1 Obseravtions

The SuperCOSMOS facility at the Royal Observatory Edinburgh has (using plates from the United Kingdom Schmidt Telescope) produced complete southern sky surveys in $`B_J`$, $`R`$ and $`I`$ with an additional second epoch $`R`$ survey. These surveys are now publicly available (Hambly et al 2001). The scanning program has now moved on to the northern hemisphere, using film and glass copies of plates taken by the Oschin Schmidt Telescope on Mount Palomar, California. These data will soon be publicly available. The declination of Alpha Per ($`\delta 49^{}`$) and the Pleiades ($`\delta 24^{}`$) means this prerelease northern hemisphere data must be used. Details of these plates are given in Table 1 for Alpha Per and Table 2 for the Pleiades. Unfortunately due to saturation caused by the bright core stars of the cluster and a satellite track, 3.1% (2.3 sq deg) of the Pleiades survey area was unusable. The region of the CO cloud near Merope was used in the proper motion survey, it is assumed that the increased reddening in this region will not affect the results. The areas covered by both surveys are shown in figures 1 and 2.

### 2.2 Sample selection

The sample of stars used was selected to maximise completeness and reliability (e.g. by minimising astrometric errors). Highly elliptical objects, non-stellar objects and poor quality images near bright stars were all removed. Deblended images were included. This is because Alpha Per (lying at a galactic latitude of $`7^{}`$) is in a crowded field meaning many stellar images will be merged with others on the plate. Unfortunately the deblending algorithm is rather crude and this can lead to increased astrometric and photometric errors. It is necessary to calculate the astrometric errors caused by the deblending alogorithm before deciding whether to include deblended objects.

### 2.3 Errors in the proper motions

In the Alpha Per data set there were several overlap regions between the plates used in this study. Hence objects in this region would have two different measurements of their position. These were used to calculate the RMS error on the proper motions. The errors were calculated in both axes ($`\xi `$ and $`\eta `$) and for objects that had been deblended and those that had not. Figure 3 shows how these errors vary with magnitude. Notice the higher errors at the bright end and the faint end. The overlap regions in the Pleiades were too small to be used for a similar calculation. Hence a small selection of faint blue stars was selected. These are expected to have negligible proper motions, so any measured proper motion would be the result of random errors. This yielded an error estimates of $`\sigma _x=6.2`$ mas/yr and $`\sigma _y=6.3`$ mas/yr . It should be noted that these errors are high because the stars used are faint, Figure 3 shows that for fainter stars errors are higher. When a study was done using candidate stars from other studies for Alpha Per it was found that 172 out of 323 were merged with other images. As it has been shown that there is no significant increase in the errors for deblended stars it was decided to include these to maximise completeness.

### 2.4 Photometric calibration

The photometric calibration used by the SuperCOSMOS survey is subject to large systematic errors especially at the bright end. This situation can be improved when large numbers of photoelectric measurements are availible, as is typically the case in open clusters, using a process of photometric recalibration. To do this the photographic magnitudes from the sample were compared with their photoelectric magnitudes. The photoelectric data for recalibrating the Alpha Per data set was provided by John Stauffer and with others being taken from Barrado y Navascués et al (2003). Figure 4 shows a plot of photoelectric magnitude vs. photographic, the line shown is a the fourth degree polynomial fit to the data. The best fit line shown in the graph was then used to recalibrate the photographic magnitudes. Unfortunately the data from Barrado y Navascués’ only includes stars which fall in one field (200). Hence the recalibration was carried out on this plate and propagated through to the other plates. The scatter from the fitted polynomial was used to estimate the photometric errors to be, $`\sigma _R=0.20`$ and $`\sigma _I=0.20`$. A similar process was carried out for the Pleiades using data from Jameson & Skillen (1989) and Stauffer (1984). Figure 5 shows a plot of photoelectric magnitude vs. photographic, again the line shown is a the fourth degree polynomial fit to the data. The photometric errors were found to be, $`\sigma _R=0.17`$ and $`\sigma _I=0.21`$.

### 2.5 Colour selection

After photometric recalibration, data cuts could be made on the basis of colour and brightness. Firstly it was decided that as stars brighter than $`R=10`$ produce high centroiding errors these would be excluded. In the Alpha Per study the low separation between the cluster and the field in the Proper Motion Vector Point Diagram produced a large amount of background star contamination at the faint end. Hence stars fainter than $`R=18`$ were excluded. With it’s better separation between the cluster and field such a cut was not necessary for the Pleiades data set. To reduce the amount of background star contamination colour selections were made. Figures 6 and 7 show colour magnitude diagrams for each sample. The dots represent stars in the sample while the star symbols are records from the sample which have been identified as candidates from both our study and previous studies, the triangles are star thought to be candidiates by previous studies but not by our study. It can clearly be seen that the candidate stars from other surveys form a main sequence-like line distinct from the large mass of background stars. Only stars above the line shown were used in this study. In the case of the Pleiades figure 7 shows the colour magnitude diagram. The large bulk of background stars are easily identified and removed. The line shown is again the colour cut used. Both lines used in the colour cuts were chosen by inspection As astrometric errors vary with magnitude the data was divided into the following magnitude ranges, for Alpha Per $`10<R<12`$, $`12<R<14`$, $`14<R<16`$ and $`16<R<18`$ and for the Pleiades, $`10<R<12`$, $`12<R<14`$, $`14<R<16`$, $`16<R<17.5`$ and $`R>17.5`$ . Each of these was treated independently in the fitting process.

### 2.6 The fitting process

The membership probabilities were calculated by a process similar to that outlined by Sanders (1971). The full mathematical detail is outlined in Appendix A. To ensure that the data fitting process is robust a series of twenty sets of simulated data were produced and run through the fitting program. The data sets were generated by taking uniform random number distributions. These were converted into gaussian and exponential distributions with the correct scale lengths, means and variances to model the expected distributions for the field and cluster stars. Each data set included 350 cluster stars and 350 field stars, both typical numbers for the real data. Table 3 shows the results obtained. It is apparent that there is no significant offset in any of the calculated parameter values.

## 3 Results

The proper motion vector point diagrams for each of the magnitude ranges for both clusters are shown in Figure 8 (the Pleiades) and Figure 9 (Alpha Per) along with probability histograms. In the case of Alpha Per the cluster is especially apparent in the vector point diagram $`16<R<18`$ and consequently it also has the highest number of high probability cluster members. As the Pleiades study goes far deeper we can see the cluster becoming more apparent with increasing magnitude before becoming engulfed by the mass of field stars in the lowest magnitude interval. Accordingly the number of high probability members increases with increasing magnitude, the high number of field stars in the lowest magnitude interval means there are no stars with membership probabilities $`>90\%`$. A full list of stars with calculated membership probabilities greater than $`60\%`$ is available for both clusters from CDS in Strasbourg. Example tables are shown in Appendix B.

## 4 Analysis

### 4.1 Luminosity Functions

As Schmidt plates have poor detection rates close to the plate limit it was necessary to carry out a completeness estimate. This was done by measuring how the number of stars in the full stellar sample increased with $`R`$ magnitude. For a uniform distribution of stars the number in each sample should rise exponentially with increasing magnitude with any dropoff being caused by incompleteness. An estimate of this incompleteness can be found by taking the logarithm of the number of stars in each magnitude interval, fitting a best fit line up to the point where the dropoff begins and then using the deficit to estimate the incompleteness. Figure 10 shows such a fit for the Alpha Per data set and Figure 11 that for the Pleiades sample.
The luminosity functions were found by summing the membership probabilities of stars in each one magnitude wide interval. Each interval was then corrected for incompleteness. Figure 12 shows the Luminosity Function produced for Alpha Per. Binarity was not taken into account when producing this. The luminosity function is smooth and rises with increasing magnitude. The Pleiades luminostiy function (Figure 13) by contrast has distinct features, the most obvious of which is the peak at $`16<R<17`$. The errors shown are Poisson errors which also take into account errors in the number of backround stars.

### 4.2 Mass Function

To convert the luminosty function to a mass function a mass luminosity relationship is required. This is provided by the models of Barraffe et al (1998). Figure 14 shows the mass luminosity relation taken from these models for a 90Myr old population such as Alpha Per. The line shown is a sixth degree polynomial least squares fit to the data which is used to find the mass of a star of a particular luminosity. The mass function is found from the luminosity function, $`\mathrm{\Phi }(M)`$. In some particular magnitude interval of width $`dM`$ there will be a number of stars $`dN=\mathrm{\Phi }(M)dM`$. In the corresponding mass interval the same number $`dN`$ is given by $`dN=\xi (m)dm`$, where $`dm`$ is the width of the interval and $`\xi (m)`$ is the mass function. This can be calculated by taking the number of stars in a particular luminosity interval and dividing by the width of that interval in mass. Again binarity was not taken into account. The mass function for Alpha Per is shown in Figure 17. The mass function was then fitted by a power law fit of the form,

$$\xi (m)m^\alpha $$
(1)
It was found that $`\alpha =0.86_{0.19}^{+0.14}`$. As mentioned earlier Barrado y Navascués’ et al (2002) derived a mass function for Alpha Per at lower masses (down to $`M=0.035M_{}`$). They also fitted a power law function of the form, finding that $`\alpha =0.59`$. It therefore appears that the mass function flattens off significantly towards lower masses.
A mass function was also produced for the Pleiades. This time the mass luminosity relation was that for a 120 Myr old population and is shown in Figure 15. Again this was taken from Barraffe et al (1998). The mass function is shown in Figure 16, the data points from this survey are supplemented by those from Hambly et al (1999) and references therein. It has been suggested by Adams and Fatuzzo (1996) that the initial stellar mass function can be approximated by a log normal function as shown below,

$$log_{10}\xi (m)=a_0+a_1log_{10}m+a_2(log_{10}m)^2$$
(2)
It is clear that a power law fit is not appropriate for the mass function derived here; a log normal was fitted to the Mass function. The parameters are $`a_0=2.213`$, $`a_1=2.069`$ and $`a_2=0.745`$. Hambly et al (1999) produced a similar log normal fit to the mass function with parameters, $`a_0=2.341`$,$`a_1=2.313`$ and $`a_2=1.191`$. These are in generally good agreement with the discrepancy in $`a_2`$ probably due to the poorer constaint of the Hambly study at the faint end. The referee has pointed out the possibilty of systematic errors in the mass-luminosity relationship from the Baraffe et al (1998) models. We checked for such errors by taking an empirical $`BC_R`$ vs $`T_{eff}`$ relation derived using the same data and techniques as described in Hambly et al. (1999) and using it to estimate the $`R`$ magnitude from $`M_{bol}`$ for each of the points taken from Baraffe’s model. The mass of these points was then compared to that given by the mass-luminosity relation for the recalculated $`R`$. The results of this are shown in Figure 18.

### 4.3 Comparison with other studies

As mention earlier there are several other methods to determine cluster membership. Examining the membership probabilities of stars from other studies can help both to confirm individual star’s cluster membership and to examine the relevance of both study’s results. It would be worrying if candidates from previous studies did not have, by and large, high membership probabilities. The candidate stars were paired with the appropriate records in the catalogue. The maximum pairing radius was set to 6 arcseconds. A table of those candidate stars for both clusters with calculated membership probabilities is availible from CDS in Strasbourg, examples are provided in Appendix B. Figure 19 shows a histogram of the membership probabilities for these candidate stars for Alpha Per while Figure 20 shows the same graph for the Pleiades study. It is encouraging that both show that most stars thought to be members in previous studies have high membership probabilities and that the large number of probable nonmembers shown in the membership probability histograms (see Figure 8 and Figure 9) do not appear here. Both Figure 19 and Figure 20 show a sharp increase in the number of candidate stars from previous studies at p=60%. For this reason it was decided to include only stars with a membership probability greater than this in the catalogues of probable members for both clusters (available on CDS). It should be noted that several HHJ objects (see Hambly, Jameson & Hawkins 1991) do not have calculated membership probabilities. This is because these objects appeared highly elliptical on the first epoch plates SuperCOSMOS measurements and hence were not included in the present proper motion survey.

## 5 Conclusion

We have presented the first proper motion survey of the cluster Alpha Persei to calculate formal membership probabilities. A Luminosity Function for the range $`R=10`$ to $`R=18`$ has been produced and a Mass Fuction derived from this. The Mass Function was fitted with a power law distribution with $`\alpha =0.86_{0.19}^{+0.14}`$. A catalogue of 339 high probability members has been created. This has been cross-correlated with previous studies of the cluster showing that a large number of stars previously thought to be cluster members have high calculated membership probabilities. A similar study has also been undertaken on the Pleiades and while not the first of it’s type it has produced several low mass star candidates and yielded a better constrained Mass Function at the faint end.

###### Acknowledgements.

The authors would like to thank John Stauffer for supplying data on candidate member stars from other surveys and for a prompt and thorough referee’s report.

## Appendix A Mathematics

In order to calculate membership probabilities, distributions were fitted to both the cluster stars and the field stars in the vector point diagram. The cluster stars were fitted with a circularly symmetric gaussian distribution as in Sanders (1971) however the field stars were fitted with a decaying exponential in the direction of cluster proper motion and a gaussian perpendicular to this direction as in Hambly et al (1995). In order to fit a distribution for the field stars the vector point diagram was rotated so that the cluster lay on the y axis. The total distribution function $`\mathrm{\Phi }`$ is given by equation 3 where $`\mathrm{\Phi }_f`$ is the field star distribution, $`\mathrm{\Phi }_c`$ is the cluster star distribution and $`f`$ is the fraction of stars which are field stars. $`\mathrm{\Phi }_f`$ is given in (4) and $`\mathrm{\Phi }_c`$ is given in (5). $`c_o`$ is the normalisation for the exponential between the limits $`\mu _1`$ and $`\mu _2`$and is given by (7).

$$\mathrm{\Phi }=f\mathrm{\Phi }_f+(1f)\mathrm{\Phi }_c$$
(3)
$$\mathrm{\Phi }_f=\frac{c_o}{\sqrt{2\pi }\mathrm{\Sigma }_x}exp\left(\frac{(\mu _x\mu _{xf})^2}{2\mathrm{\Sigma }_x^2}\frac{\mu _y}{\tau }\right)$$
(4)
$$\mathrm{\Phi }_c=\frac{1}{2\pi \sigma ^2}exp\left(\frac{(\mu _x\mu _{xc})^2+(\mu _y\mu _{yc})^2}{2\sigma ^2}\right)$$
(5)
$$c_o=\frac{1}{\tau (e^{\frac{\mu _2}{\tau }}e^{\frac{\mu _1}{\tau }})}$$
(6)
$$c_o_{\mu _1}^{\mu _2}e^{\mu _y/\tau }𝑑\mu _y=1$$
(7)
Using the maximum likelihood method

$$\underset{i}{}\frac{\delta ln\mathrm{\Phi }_i}{\delta \mathrm{\Theta }}=0$$
(8)
(where $`\mathrm{\Theta }`$ is some parameter) the following set of nonlinear equations are found.

$$f:\underset{i}{}\frac{\mathrm{\Phi }_f\mathrm{\Phi }_c}{\mathrm{\Phi }}=0$$
(9)
$$\sigma :\underset{i}{}\frac{\mathrm{\Phi }_c}{\mathrm{\Phi }}\left(\frac{(\mu _x\mu _{xc})^2+(\mu _y\mu _{yc})^2}{\sigma ^2}2\right)=0$$
(10)
$$\mathrm{\Sigma }_x:\underset{i}{}\frac{\mathrm{\Phi }_f}{\mathrm{\Phi }}\left(\frac{(\mu _x\mu _{xf})^2}{\mathrm{\Sigma }_x^2}1\right)=0$$
(11)
$$\mu _{xf}:\underset{i}{}\frac{\mathrm{\Phi }_f}{\mathrm{\Phi }}(\mu _x\mu _{xf})=0$$
(12)
$$\mu _{xc}:\underset{i}{}\frac{\mathrm{\Phi }_c}{\mathrm{\Phi }}(\mu _x\mu _{xc})=0$$
(13)
$$\mu _{yc}:\underset{i}{}\frac{\mathrm{\Phi }_c}{\mathrm{\Phi }}(\mu _y\mu _{yc})=0$$
(14)
$$\tau :\underset{i}{}\frac{\mathrm{\Phi }_f}{\mathrm{\Phi }}\left(\frac{\mu _y}{\tau }1c_o(\mu _1e^{\frac{\mu _1}{\tau }}\mu _2e^{\frac{\mu _2}{\tau }})\right)=0$$
(15)
These are solved by a simple bisection algorithm. Applying the calculated values of the parameters, the probability of the ith star being a cluster member is,

$$p_i=\frac{(1f)\mathrm{\Phi }_{ci}}{\mathrm{\Phi }_i}$$
(16)
The fitted values for the parameters for each magnitude interval are shown in table 4 for Alpha Per and table 5 for the Pleiades.

## Appendix B Example Tables

Example Tables of the catalogue produced. All coordinates are J2000 and all magnitudes are recalibrated photographic magnitudes in the natural photographic ($`R_{59F}`$, $`I_{IVN}`$) system. Some membership probabilities may be underestimated at the bright end due to the low separation of the main sequence from the background stars at such magnitudes in the colour magnitude diagram.


text:
# Electronic Structure of Transition-Metal Dicyanamides Me[N(CN)2]2 (Me = Mn, Fe, Co, Ni, Cu)

## I Introduction

Molecule-based magnets have been a subject of interest in the past decade, after the discovery of room-temperature magnetic ordering in some of these materials.Manriquez ; Ferlay  Interest is driven not only by the potential for use in applications like information storage and processing, but also by fundamental questions about the origin of magnetic ordering in these materials. The transition-metal dicyanamides, Me\[N(CN)2\]2 where Me is a $`3d`$ transition metal, provide an example where, even within an isostructural series of compounds, the magnetic ordering can differ substantially depending on the metal species and details of the crystal structure. Experimental studies of the structural and magnetic properties of the dicyanamides with Me = Mn, Fe, Co, Ni, and Cu have been performed by several groups.Kmety3 ; Batten ; Kurmoo ; Kmety1 ; Kmety2 ; thesis  Mn\[N(CN)2\]2 and Fe\[N(CN)2\]2 have been found to be canted antiferromagnets, while Co\[N(CN)2\]2 and Ni\[N(CN)2\]2 have been found to be ferromagnets. Cu\[N(CN)2\]2 has been reported to be paramagnetic, at least down to a few Kelvin, and a recent experiment has found that it becomes a ferromagnet below 1.8 K. thesis  Thus the compounds with transition-metal ions having six or fewer $`3d`$ electrons order as canted antiferromagnets, while those with seven or more $`3d`$ electrons order as ferromagnets. The Curie (or Néel) temperatures of these materials range from 1.8 K for Cu\[N(CN)2\]2 to 22.7 K for Ni\[N(CN)2\]2. Ni\[N(CN)2\]2 and Fe\[N(CN)2\]2 have been found to have the largest coercive fields of all known metal-organic magnets – 7975 and 17800 Oe, respectively.

The molecular building block of the bulk dicyanamides consists of a divalent metal ion that is two-fold coordinated by two N(CN)2 groups. When the molecular units pack to form the bulk crystal, the negatively charged CN groups at the end of neighboring molecules interact with the positive central metal ion in such a way that the metal ion becomes surrounded by an octahedron of six N atoms. The local symmetry around the metal site is tetragonal, with the original molecular axis corresponding to the long axis. The primitive cell contains two formula units with their long axes rotated in opposite directions in the $`xy`$ plane, leading to an orthorhombic cell. The metal ions are the source of magnetic moments while the organic species provide superexchange pathways between the magnetic centers.

In this paper, we report a combined experimental and theoretical study of the electronic structure of Me\[N(CN)2\]2 with Me = Mn, Fe, Co, Ni, Cu. We have measured x-ray emission spectra (XES) and x-ray photoemission spectra (XPS) across the series. Generally, XPS probes the total occupied density of states (DOS) whereas XES probes the site-selective partial DOS due to the dipole selection rules. In the case of carbon and nitrogen $`K_\alpha `$ XES, which correspond to $`2p1s`$ transitions, the occupied $`2p`$ states can be studied. Similarly, Me $`L_{2,3}`$ XES, related to the $`3d4s2p`$ transition, gives information about occupied $`3d`$ bands. The comparison of XPS and XES in the binding-energy scale yields information about the distribution of the partial and total DOS in multicomponent systems. The measured spectra are compared with results of density-functional calculations. In a recent paper focusing on the electronic structure of Mn\[N(CN)2\]2, reasonable agreement was found between theory and experiment for the electronic structure, magnetic ordering, and magnetic anisotropy.Pederson  Here we find that the direct comparison of density-functional calculations with experimental data is not as straightforward for some of the other compounds in the series.

## II Experimental and Calculation Details

The x-ray fluorescence measurements were performed at Beamline 8.0 at the Advanced Light Source (ALS) of Lawrence Berkeley National Laboratory (LBNL). Carbon and nitrogen $`K_\alpha `$ ($`2p1s`$ transition) and Me $`L_{2,3}`$ ($`3d4s2p`$ transition) x-ray emission spectra were taken, employing the University of Tennessee at Knoxville soft x-ray fluorescence endstation.Jia  Photons with energy of 300 eV above the carbon $`K`$ edge were delivered to the endstation via the beamline’s 89-period, 5 centimeter undulator insertion device and spherical monochromator. For nitrogen $`K_\alpha `$ and Me $`L_{2,3}`$ XES, an excitation-energy dependence was measured near the N $`1s`$ and Me $`2p`$ thresholds. The carbon and nitrogen $`K_\alpha `$ spectra were obtained with an energy resolution of 0.3-0.4 eV. The Me $`L_{2,3}`$ XES were measured with an energy resolution of 0.7-0.8 eV.

The XPS measurements were performed with an ESCA spectrometer of Physical Electronics (PHI 5600 ci) with monochromatized Al $`K_\alpha `$ radiation of a 0.3 eV FWHM. The energy resolution of the analyzer was 1.5% of the pass energy. The pressure in the vacuum chamber during the measurements was below $`5\times 10^9`$ mbar. Prior to XPS measurements, the samples were fractured in ultra-high vacuum. All photoemission studies were performed at room temperature on freshly cleaved surfaces. Charging of the insulating Me\[N(CN)2\]2 samples was compensated with an electron gun. Because there is no saturated carbon in these molecules, the binding energies can not be corrected using the C $`1s`$ line. Therefore we have used an alternative wayMoulder  of fixing the N $`1s`$ line of the cyanide group at 399.6 eV.

Polycrystalline pressed pellets of Me\[N(CN)2\]2 compounds were used for XES and XPS measurements. Details about the sample preparation can be found in Refs. Kmety1,  and Kmety2, .

Diffuse reflectance spectroscopy on powder samples was used to measure the lowest electronic transition observable in the visible region, which can be interpreted as possible band gap energies. The reflectance spectra were taken on powdered samples at room temperature in a 1-mm quartz cuvette which was illuminated by light from a 75-watt Xe high-pressure arc lamp. The reflected light was admitted into an Acton Spectropro 300i 300 mm spectrometer with a Hammamatsu CCD detector, which gave light intensities versus wavelength from 300 to 800 nm. The spectrum of a “white” standard, consisting of powdered NaCl, was used for determination of the absorbance, $`\mathrm{log}(I_0/I_s)`$, as a function of wavelength for each sample, where $`I_0`$ is the intensity of the “white” standard and $`I_s`$ is the intensity of the sample compound.

For comparison with the experimental results, electronic-structure calculations were performed for ferromagnetic (FM), antiferromagnetic (AF), and nonmagnetic (NM) phases using the crystal structures determined from x-ray and neutron diffraction experiments.Kurmoo ; Kmety1 ; Kmety2 ; Kmety3 ; thesis  Spins were restricted to be collinear. Calculations were carried out using the VASP program,vasp  a density-functional-based code employing planewave basis sets and ultrasoft pseudopotentials. The electron-electron interaction was treated in the generalized gradient approximation (GGA),GGA1 ; GGA2  and electronic wave functions were expanded in planewaves up to a kinetic energy cutoff of 435 eV. Monkhorst-Pack meshesMonkhorst-Pack  of up to $`8\times 8\times 8`$ k-points were used to sample the Brillouin zone, and the linear tetrahedron method was used to integrate over the Brillouin zone.

## III Magnetic Properties

Table 1 summarizes the results of the total-energy calculations. Density-functional-theory (DFT) calculations for terminated cluster models yield similar results.Pederson1  Experimentally, Mn\[N(CN)2\]2 and Fe\[N(CN)2\]2 are both canted antiferromagnets. In the Mn compound, the canting angle is reported to be at most a few degrees,Kmety1  while in the Fe compound, the canting angle is estimated to be 20-23.thesis  Our calculations find Mn\[N(CN)2\]2 to be antiferromagnetic, but neglect of both the canting and spin-orbit interaction in Fe\[N(CN)2\]2 results in the FM phase being favored. Similar results were recently obtained in Ref. Ruiz, . We expect that when spins are not restricted to be collinear, density functional theory would yield an AF ground state for Fe\[N(CN)2\]2. For Ni\[N(CN)2\]2 and Co\[N(CN)2\]2, we find the ground state to be ferromagnetic, in agreement with experiments. In the case of Cu\[N(CN)2\]2, the FM and AF phases are calculated to be very close in energy and favored over the NM phase. Magnetization studies show no long-range ordering above 2 K,Kurmoo  but a recent study thesis  finds that Cu\[N(CN)2\]2 becomes ferromagnetic at lower temperatures with $`T_c1.8`$ K. Hence we report results of our calculations for the FM phase.

In superexchange systems, the magnitude and sign of the one-electron contribution to the exchange interaction depend on the filling of those orbitals on neighboring magnetic centers that interact via a bridging ligand. Anderson ; Goodenough ; Gudel  The relevant orbitals, i.e., those with non-zero interaction matrix elements, depend in turn on the geometry of the superexchange pathway. The magnetic ordering is thus determined by an interplay between the occupation of the $`d`$ shell and the superexchange angle. In Me\[N(CN)2\]2, neighboring Me sites are connected through a Me-\[N-C-N\]-Me pathway, where the two N sites along the path are inequivalent. The Me-\[N-C-N\]-Me angle increases monotonically in going from Mn\[N(CN)2\]2 to Ni\[N(CN)2\]2, with a total change of about 2.5 degrees. Kmety1  It has been suggested that this structural change in the pathway is responsible for the different magnetic ground states in the dicyanamides, with the angle passing through a critical value in going across the $`3d`$ series. Kmety1  Our calculations find that such small changes in structure do not change the energetic ordering of the magnetic phases. The occupation of the $`3d`$ shell is evidently a more significant factor in determining the sign of the superexchange interaction in these systems.

Density functional theory typically underestimates magnetic moments, but, as shown in Table 1, discrepancies for Me = Fe and Ni appear unusually large. This is likely because the measured magnetic moments in the table, determined from neutron powder diffraction studies,Kmety1 ; thesis  include both spin and orbital contributions. For octahedrally coordinated Me2+ cations without spin-orbit coupling, the orbital moment is expected to be completely quenched for Me = Mn, Ni, and Cu, and unquenched or partially quenched for Me = Fe and Co.Figgis  With spin-orbit coupling, the quenching due to the crystal field is incomplete, resulting in a Landé g-factor that can deviate from 2.thesis ; Kmety2  With no mechanism for unquenching the orbital moment within our calculations, a direct comparison between the theoretical and experimental moments in Table 1 is not appropriate.

It is well known that the local spin magnetic moment is related to spectral splitting of the $`3s`$ core-level x-ray photoemission spectra in transition metals and their compounds. The splitting originates from the exchange coupling between the $`3s`$ hole and the $`3d`$ electrons and the magnitude of the splitting is proportional to $`(2S+1)`$, where $`S`$ is the local spin of the $`3d`$ electrons in the ground state. VanVleck  The $`3s`$ spectra of the Me\[N(CN)2\]2 materials together with those of simple oxides MnO, FeO, CoO, NiO, and CuO are shown in Fig. 1. Peaks A and B denote $`3s^13d^{n+1}L^1`$ final-state configurations, and C and D denote $`3s^13d^n`$ final-state configurations, where $`L^1`$ denotes a ligand hole due to the charge-transfer process. For Mn\[N(CN)2\]2, fitting curves and their component peaks are shown. The Mn $`3s`$ spectrum shows two peaks, C and D, which can be attributed to the $`3s^13d^5`$ final state. The magnitude of the $`3s`$ exchange splitting in Mn\[N(CN)2\]2, 6.5 eV, is very similar to that in MnO, 6.2 eV, Galakhov2  which indicates the Mn $`3d^5`$ ground-state configuration for Mn\[N(CN)2\]2. This correlates well with calculated and measured spin moments of Mn\[N(CN)2\]2 presented in Table 1. The Me $`3s`$ spectra of Ni\[N(CN)2\]2, Co\[N(CN)2\]2, and Fe\[N(CN)2\]2 show a complex structure formed by both $`3s^13d^n`$ and $`3s^13d^{n+1}L^1`$ final-state configurations. Coexistence of both exchange and charge-transfer effects makes the comparison of the Me $`3s`$ splitting with the local magnetic moments difficult. Similarities with the spectra of the monoxides suggest that the Me $`3s`$ energy splitting is systematically reduced with increasing Me atomic number, as expected for transition metal ions with nominal $`2^+`$ valence (Table 1). For Cu\[N(CN)2\]2, the exchange splitting is expected for the $`3s^13d^9`$ configuration only, which is localized in the energy interval from 130 to 135 eV. Low energy resolution limits our ability to separate this splitting in the present data, but it has been observed in copper oxides (see Fig. 1 and Ref. Galakhov, ).

## IV Electronic Structure

Figure 2 shows the measured carbon and nitrogen $`K_\alpha `$ XES of the Me\[N(CN)2\]2 compounds. As seen, the carbon $`K_\alpha `$ XES (Fig. 2a) are very similar across the series because of the same local atomic and electronic structure of carbon atoms in cyano groups. On the other hand, the nitrogen atoms occupy two non-equivalent sites in the crystal structure of Me\[N(CN)2\]2 compounds, where each Me atom has two N(1) neighbors along the long axis and four N(2) neighbors from CN anions along the other axes.Kmety1  The nitrogen $`K_\alpha `$ XES of Me\[N(CN)2\]2 (Fig. 2b) has contributions from both sites. The N $`K_\alpha `$ spectrum of the reference compound Na4Fe(CN)6(OH2O) (Fig. 2b) shows N contributions from cyano groups and is quite different from the non-resonant nitrogen spectra of Me\[N(CN)2\]2. In our XPS measurements, the N $`1s`$ spectrum is split into two main lines separated by approximately 1.6 eV, corresponding to the difference in the N $`1s`$ binding energies of the two non-equivalent N sites (Table 2).

Additional measurements were made for resonantly excited nitrogen $`K_\alpha `$ XES in Mn\[N(CN)2\]2. The energy dependence of the spectra when the excitation energy was tuned near the N $`1s`$ threshold is shown in Fig. 3a. The intensities of the peaks at 389.9, 392.7 and 394.8 eV change with excitation energy. The non-resonant spectrum $`d`$ excited far from the N $`1s`$ threshold can be described as a superposition of spectra $`b`$ and $`c`$ (Fig. 3b), corresponding to contributions from selectively excited non-equivalent nitrogen atoms. Based on our measurements of the non-resonant N $`K_\alpha `$ XES of Na4Fe(CN)6(OH2O) (Fig. 2b), we can attribute the origin of the resonant spectrum of Mn\[N(CN)2\]2 excited near peaks $`a`$ and $`b`$ of the N $`1s`$ TEY to N(2) atoms located in cyano groups, whereas the spectrum excited near peak $`c`$ of the N $`1s`$ TEY can be explained by contributions from both N(1) and N(2) atoms.

The valence band (VB) XPS and the XES of the constituents of Me\[N(CN)2\]2 are compared in Figs. 4-8. We have selected for this comparison the non-resonant carbon and nitrogen $`K_\alpha `$ XES shown in Fig. 2 and the Me $`L_3`$ XES ($`3d4s2p_{3/2}`$ transition) which were excited just above the $`L_3`$ threshold to exclude the excitation of the Me $`L_2`$ XES ($`3d4s2p_{1/2}`$ transition) that overlaps the Me $`L_3`$ XES due to spin-orbit splitting. To convert the XES of the constituents to the binding-energy scale, we have subtracted emission energies from the XPS binding energies of core levels (Table 2). According to this comparison, the top of the valence band (a) is derived mainly from Me $`3d`$ and N $`2p`$ states. Towards the middle of the valence band the Me $`3d`$, N $`2p`$ and C $`2p`$ states are strongly hybridized, forming a broad structure (b-d). At the bottom of the valence band (e) are mixed carbon and nitrogen $`2p`$ states. Atomic-like carbon and nitrogen $`2s`$ states are located in region (f), with binding energy around 20-25 eV.

For comparison, we present theoretical results for Ni\[N(CN)2\]2 as an example. The calculated electronic density of states (DOS) of Ni\[N(CN)2\]2 is shown in Fig. 9. The total DOS is in the lowest panel, and site- and $`l`$-projected DOS are plotted in the other panels. In the lowest panel, a Gaussian broadening of 1.0 eV has been introduced to simulate the effect of finite instrumental resolution. The partial DOS in the upper three panels have been broadened by 0.2 eV to help in identifying features seen in the experimental XES spectra. Features in the DOS are labeled (a)-(f) as in the experimental XPS spectra. The nature of the states at different energies is generally consistent with what is found in the XPS/XES comparison. Regions (a) and (b) near the top of the valence band are dominated by Ni $`3d`$ states, with some hybridization with C and N $`2p`$ states, and they correspond to the peak in the Ni $`L_3`$ XES spectrum. Regions (c) and (d) contain strong contributions from C $`2p`$ and N $`2p`$ states, as seen in the C and N XES spectra, as well as some Ni $`3d`$ states. Peak (e), approximately 12 eV below the top of the valence band, is primarily of C $`2p`$ and C $`2s`$ character, and N $`2p`$ contributions from N(1) sites. Region (f) has significant contributions from the N $`2s`$, C $`2s`$, and C $`2p`$ states and is observed in the C $`K_\alpha `$ XES spectrum.

For the other compounds, the electronic DOS is very similar to that of Ni\[N(CN)2\]2 from about 25 eV below the top of the valence band up to region (d) where $`3d`$ states of the transition metal start playing a significant role. The hybridization between N $`2p`$ and Me $`3d`$ states increases with the Me atomic number. For Me = Mn, Fe, and Co, states near the top of the valence band are overwhelmingly of Me $`d`$ character, while for Me = Ni and Cu, the relative contribution of N $`p`$ states is increased. In both the calculations and the experiments, the weight of the $`3d`$-derived bands shifts with respect to the top of the valence band. Based on the XPS valence bands, the maximum of the $`3d`$ bands varies from about 4 eV for Mn\[N(CN)2\]2 to about 6 eV for Cu\[N(CN)2\]2 (see Figs. 4-8 and Table 2). This result is confirmed by the Me $`L_3`$ XES converted to the binding-energy scale. In addition, in the calculations, we find, for example, about a 1.5 eV difference in the position of the main maxima of the occupied $`d`$-DOS weights in Cu\[N(CN)2\]2 compared to Ni\[N(CN)2\]2 (Fig. 10), which is similar to what is observed experimentally (Table 2). The deepening of $`3d`$ states found for Cu\[N(CN)2\]2 is due to the localization $`d`$ electrons at the end of the $`3d`$ series. This is observed in molecular crystals because the $`dd`$ and $`dsp`$ interactions are weaker than in conventional divalent (and even monovalent) copper compounds. The decreasing exchange splitting is likely an additional factor affecting energy deepening of $`3d`$ states.

The deepening of $`d`$ states across the series is accompanied by a shift towards higher binding energies of the maxima of the Me $`2p`$ and Me $`3s`$ XPS spectra of Me\[N(CN)2\]2 in comparison with those for $`3d`$ monoxides.Galakhov2 ; Galakhov  This shift increases from about 0 eV for Mn\[N(CN)2\]2 to about 3 eV for Cu\[N(CN)2\]2. Due to weak metal-metal interactions in molecular crystals, especially in Cu\[N(CN)2\]2 where $`3d`$ states are strongly localized, the screening of the final states is weaker than in the other cases, resulting in larger binding energies of the valence bands and the core levels.

At the top of the valence band the detailed electronic structure is determined by the filling of the transition metal $`3d`$ states. The local environment around Me sites is a distorted octahedron of N sites, and the ligand field splits the $`3d`$ states by approximately 1.5 eV into $`t_{2g}`$ and $`e_g`$ manifolds. In Mn\[N(CN)2\]2, with an exchange splitting of over 3 eV, the majority $`d`$ bands are fully occupied and there is a significant gap to the unoccupied minority $`d`$ bands.Pederson  As the number of $`d`$ electrons increases, the ligand-field splitting stays about the same while the exchange splitting decreases. The majority bands nevertheless remain full. The minority $`t_{2g}`$ bands are narrow and split by the tetragonal distortion of the local octahedral environment and by the longer-range crystal field. In the FM configuration, intersublattice $`dd`$ interactions broaden the majority $`t_{2g}`$ bands, resulting in a two-fold subband of $`d_{xz,yz}`$ states with a gap to the $`d_{xy}`$ subband. (Axes are referenced to the local octahedral environment.) In the AF configuration, intersublattice $`dd`$ interactions are suppressed, which narrows the $`t_{2g}`$ bands. For Me = Fe, the ordering of the one-fold and two-fold subbands reverses in the FM phase. Thus, with one minority $`t_{2g}`$ electron per Fe site, AF Fe\[N(CN)2\]2 is metallic while FM Fe\[N(CN)2\]2 is insulating, as shown in Fig. 11. The lack of a gap in AF Fe\[N(CN)2\]2 with collinear spins contributes to its instability with respect to the insulating FM phase. In FM Co\[N(CN)2\]2, the two-fold subband fills, and in FM Ni\[N(CN)2\]2, the $`t_{2g}`$ manifold is completely occupied. The Jahn-Teller effect significantly enhances the tetragonal distortion of the octahedral environment around Cu sites in Cu\[N(CN)2\]2,Kurmoo ; Batten  which increases the splitting between $`e_g`$ bands. The calculated GGA band gaps are listed in Table 1, along with estimated gaps from diffuse reflectance spectroscopy. The agreement is better for the compounds where $`t_{2g}`$ band is either completely empty or completely full.

The narrow width of some of the $`d`$ bands (Fig. 11) suggests that electron correlations could be important in these materials. Additional evidence for this is the GGA’s large underestimation of the band gaps in some of the compounds, particularly Co\[N(CN)2\]2. Furthermore, the fact that the measured magnetic moments in Fe\[N(CN)2\]2 and Ni\[N(CN)2\]2 are larger than the maximum moments that could be created by spin alone indicates the importance of the orbital moment, and therefore methods incorporating orbital dependence of the potential could significantly improve the results. The situation is reminiscent of the monoxides of the $`3d`$ transition metals. As with the oxides, it is likely that methods beyond DFT-GGA are needed for a better description of the electronic and magnetic structure of the dicyanamides in the solid state.

## V Conclusions

We have carried out a theoretical and experimental investigation of the electronic structure of the Me\[N(CN)2\]2 (Me = Mn, Fe, Co, Ni, Cu) molecular magnets. The general features of the valence bands are similar across the series, with the primary exceptions being the location of the $`3d`$ bands with respect to the Fermi level and the degree of hybridization between Me $`3d`$ states and N $`2p`$ states. Going across the row from Mn to Cu, the $`d`$ bands lie deeper in energy, and the degree of hybridization with N $`p`$ states increases.

Total-energy calculations indicate that the small differences in structural parameters (for example, the angle of the superexchange pathway) found in these isostructural materials is not sufficient to account for differences in their magnetic ground states. Rather, the occupation of the $`3d`$ shell is likely a more significant factor.

For most of the materials in the series, our calculations underestimate either the band gap or the local magnetic moment. Together with the existence of narrow d-bands in these materials, this suggests that a better description of the magnetic and electronic structure of these materials will require a more accurate treatment of electron-electron interactions.

###### Acknowledgements.

The authors thank M. R. Pederson for helpful discussions and for communicating unpublished results. This work was supported by the Research Council of the President of the Russian Federation under Grant NSH-1026.2003.2, the Russian Foundation for Basic Research under Project 02-02-16674, the President’s NSERC fund of the University of Saskatchewan, Deutsche Forschungsgemeinschaft (DFG) under the priority programme SPP 1137 Molekularer Magnetismus, the US National Science Foundation under Grant DMR-0210717, and the US Office of Naval Research under Grant N00014-02-1-1046. Work at the Advanced Light Source at Lawrence Berkeley National Laboratory was supported by the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. Work at the Advanced Photon Source is supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Contract No. W-31-109-ENG-38.