text:
# The stellar velocity distribution in the solar neighbourhood

## 1 Introduction

The velocity distribution function (df) of stars in the solar neighbourhood provides unique insights into the Galactic potential field, the dynamical history of the Galactic disk, and the relationships between kinematics, age, and metallicity for disk stars.

Let us define the Local Standard of Rest (LSR) to be a fictitious point that coincides with the Sun at the present instant and travels in a circular orbit around the Galactic centre. We introduce a rotating Cartesian coordinate system with origin at the LSR, $`x`$-axis pointing radially outward, $`y`$-axis pointing in the direction of Galactic rotation, and $`z`$-axis pointing to the south Galactic pole. The df $`f(𝐯)`$ is defined so that $`f(𝐯)d𝐯`$ is the number of stars per unit volume with velocities in the range $`[𝐯,𝐯+d𝐯]`$. The standard empirical model for the past century has been the Schwarzschild (1907) df, which is a triaxial Gaussian of the form

$$f(𝐯)\mathrm{exp}\left[\frac{1}{2}\underset{i,j=1}{\overset{3}{}}\alpha _{ij}(v_i\overline{v}_i)(v_j\overline{v}_j)\right],$$
(1)
where $`𝐯(v_x,v_y,v_z)`$. The two lowest moments of the df are the mean velocity, $`\overline{𝐯}=(\overline{v}_x,\overline{v}_y,\overline{v}_z)`$, and the velocity-dispersion tensor,

$$\sigma _{ij}^2\overline{(v_i\overline{v}_i)(v_j\overline{v}_j)};$$
(2)
where $`\overline{X}f(𝐯)X(𝐯)𝑑𝐯/f(𝐯)𝑑𝐯`$, and the matrix $`\sigma _{ij}^2`$ is the inverse of the matrix $`\alpha _{ij}`$.

In a steady-state, axisymmetric galaxy (i) the mean velocity relative to the LSR is tangential, so $`\overline{v}_x=\overline{v}_z=0`$; (ii) the tensors $`\alpha _{ij}`$ and $`\sigma _{ij}^2`$ are both diagonal; (iii) the ratio $`\sigma _{xx}^2/\sigma _{yy}^2`$ is determined by the local gravitational force and its radial gradient (e.g. Chandrasekhar 1942, Binney & Tremaine 1987; see also eq. 14 below).

The velocity dispersion of stars in the solar neighbourhood increases with age, probably because the disk is “heated” by one or more dynamical mechanisms. However, the interpretation of the observed age-velocity dispersion relation (avr) is uncertain: (i) One school models the avr as a smooth power law, $`\sigma _{xx}(t)t^p`$, and interprets this behavior as evidence for a continuous heating mechanism. Estimates of the exponent $`p`$ in the literature span the range 0.2–0.5. (ii) Some authors argue that the dispersion rises steeply with age for stars $`5\mathrm{Gyr}`$ old, and thereafter is relatively flat (Carlberg et al., 1985; Gómez et al., 1997), perhaps because the continuous heating mechanism saturates once the dispersion is large enough. (iii) A third model is that the dispersion does not increase smoothly with age. For example, Freeman (1991) and Quillen & Garnett (2001) argue that there are three discrete age groups ($`t3\mathrm{Gyr}`$, $`3\mathrm{Gyr}t10\mathrm{Gyr}`$, $`t10\mathrm{Gyr}`$) with different dispersions, but within each group there is no evidence for a correlation between dispersion and age. Such groups might arise because the continuous heating saturates after only 3 Gyr, and the higher dispersion of the oldest stars is due to a discrete event such as a merger.

A wide variety of mechanisms for disk heating has been discussed (see Lacey 1991 for a review): (i) Spitzer & Schwarzschild (1951, 1953) suggested that massive gas clouds could gravitationally scatter stars, leading to a steady increase in velocity dispersion with age, and thereby predicted the existence of giant molecular clouds (gmcs). However, heating by gmcs alone fails to explain several observations (Lacey, 1984; Villumsen, 1985; Jenkins & Binney, 1990; Lacey, 1991; Jenkins, 1992): the predicted ratio $`\sigma _{zz}/\sigma _{xx}`$ of vertical to radial dispersion may be too high, roughly 0.72 compared to the observed value of 0.5 for old stars (but see Ida, Kokubo & Makino 1993 for an opposing view); the predicted exponent in the avr is somewhat too low, $`p0.25`$; and the masses and number density of gmcs determined from CO observations are too low to explain the observed dispersion, probably by a factor of five or so but with substantial uncertainty. (ii) Transient spiral waves lead to potential fluctuations in the disk that excite the random motions of disk stars (Barbanis & Woltjer, 1967; Sellwood & Carlberg, 1984; Carlberg & Sellwood, 1985). However, spiral waves only excite the horizontal ($`x`$ and $`y`$) velocity components effectively, since their characteristic spatial and temporal scales are much larger than the amplitude or period of vertical oscillations for young stellar populations. (iii) These considerations lead naturally to a hybrid model, in which spiral waves excite non-circular velocities in the plane, and the velocities are then redistributed between horizontal and vertical motion through gmc scattering (Carlberg, 1987). The hybrid model has been investigated thoroughly by Jenkins & Binney (1990) and Jenkins (1992) using the Fokker-Planck equation. They find that they can successfully reproduce the observed axis ratio $`\sigma _{zz}/\sigma _{xx}`$, the exponent $`p`$ in the avr and the radial dispersion of old stars. (iv) Other possible heating mechanisms, all of which rely to some extent on hypothetical or poorly understood components or features of the Galaxy, include scattering by massive compact halo objects or halo substructure (Lacey & Ostriker, 1985), mergers with dwarf galaxies (Tóth & Ostriker, 1992; Walker, Mihos & Hernquist, 1996; Huang & Carlberg, 1997), or the outer Lindblad resonance from the Galactic bar (Kalnajs, 1991; Dehnen, 1999, 2000; Fux, 2001; Quillen, 2003).

The Schwarzschild df (1) only approximates the velocity distribution on the largest scales in velocity space. On smaller scales, there is substructure, which is most prominent in the youngest stars but present in stars of all ages. Discussion of substructure in the velocity df dates back to Kapteyn’s (1905) model of “two star streams”, and for decades Eggen advocated the case for substructure in the form of “moving groups” in the solar neighbourhood (Eggen 1996 and references therein).

Eggen and others have usually explained moving groups as the result of inhomogeneous star formation in the disk: in this model, stars in a moving group are born at a common place and time, and then disperse into a stream that happens to intersect the solar neighbourhood. This model predicts that stars in a moving group should have the same age, metallicity, and azimuthal velocity (i.e. the same angular momentum, since this determines their mean angular velocity). The existence and membership of these groups was controversial until the Hipparcos satellite measured reliable distances and proper motions for a large, homogeneous database of nearby stars, and verified the presence of rich substructure in the velocity df of both young and old stars, including a number of features that coincide with moving groups already identified by Eggen (e.g. Dehnen 1998, Chereul et al. 1998,1999)

In this paper we explore a quite different explanation for substructure in the velocity df. We suggest that substructure arises naturally from the same spiral gravitational fluctuations that excite the growth of the velocity dispersion. In other words, substructure is caused by homogeneous star formation in an irregular potential, as opposed to inhomogeneous star formation in a regular potential in the traditional model.

We investigate this hypothesis by simulating the evolution of the velocity df induced by transient spiral structure in a simple two-dimensional model of the local Galaxy. We restrict ourselves to two dimensions because spiral structure does not excite vertical motions efficiently, and because the velocity df appears to be well-mixed in the vertical direction (Dehnen, 1998).

Following Eggen, we shall use the term “moving group” to denote substructure in the velocity df of old stars ($`1`$ Gyr) at a given position. Unfortunately, the same term is sometimes also applied to OB associations, which are spatially localized concentrations of much younger stars (e.g., de Zeeuw et al. 1999).

Section 3 describes our simplified dynamical model. The results of our simulations are analyzed and compared to observed data in §4. Section 5 examines briefly the traditional hypothesis that substructure arises from inhomogeneous star formation. Section 6 contains a brief discussion of the closely related process of radial migration of stars, and §7 contains concluding remarks.

## 2 Properties of spiral structure

Our model depends on several properties of the Galaxy’s spiral structure, such as the number of arms, the arm strength, and the pitch angle. We are concerned with spiral structure in the disk surface density (rather than, say, in the distribution of young stars or gas). The properties of this structure in external galaxies are revealed by near-infrared images, which are dominated by the stars that contribute most of the mass (Rix & Rieke 1993, but see also Rhoads 1998 and James & Seigar 1999 for qualifications).

Rix & Zaritsky (1995) examined the K-band spiral structure of 18 face-on galaxies. They found that almost half had strong two-arm ($`m=2`$) spirals, with arm-interarm contrasts $`I_{\mathrm{max}}/I_{\mathrm{min}}2`$. If the arms are sinusoidal, with fractional amplitude $`ϵ`$ relative to the axisymmetric background, the arm-interarm contrast is $`I_{\mathrm{max}}/I_{\mathrm{min}}=(1+ϵ)/(1ϵ)`$, so a contrast of 2 corresponds to $`ϵ0.3`$.

Block & Puerari (1999) examined K-band spiral structure in 19 spirals. They found pitch angles ranging from $`8\mathrm{°}`$ to $`49\mathrm{°}`$ with a median value of $`22\mathrm{°}`$, and $`m=2`$ amplitudes ranging from $`ϵ=0.03`$ to $`ϵ=0.5`$, with a median $`ϵ=0.1`$.

Seigar & James (1998) conducted a similar survey of 45 galaxies. They found that the dominant Fourier mode usually has $`m=1`$ (36%) or $`m=2`$ (31%), and the median pitch angle of the spiral structure was $`8\mathrm{°}`$, with little or no correlation with Hubble type. They measure the strength of the spiral arms in terms of the “equivalent angle”; their median equivalent angle of $`14\mathrm{°}`$ corresponds to an amplitude $`ϵ0.07`$ for a sinusoidal $`m=2`$ spiral.

Elmegreen et al. (1999) stress that the amplitude of $`m=2`$ spiral structure in the near-infrared depends on whether the optical spiral arms are classified as flocculent or grand-design (Elmegreen, 1998). Grand-design spirals have arm-interarm contrasts of 1.5–6, corresponding to amplitude $`ϵ=0.2`$–0.7, while flocculent galaxies have contrast $`1.7`$, corresponding to $`ϵ0.25`$.

At optical wavebands, Ma et al. (1999) find that the mean pitch angle for 51 Sbc galaxies (the same Hubble type as the Galaxy) is $`15\mathrm{°}`$.

The measurement of spiral structure in our own Galaxy is more difficult than in face-on external galaxies. Drimmel (2000) uses K-band photometry of the Galactic plane to conclude that the Galaxy contains a two-arm spiral with pitch angle $`18\mathrm{°}`$. Vallée (2002) reviews a number of studies of the Galaxy’s spiral structure, mostly based on young stars, gas and dust; these yield a range of pitch angles from $`6\mathrm{°}`$$`20\mathrm{°}`$, but Vallée concludes that the best overall fit is provided by an $`m=4`$ spiral with pitch angle of $`12\mathrm{°}`$. The arm structure in the Galaxy appears to be intermediate between grand-design and flocculent (Elmegreen, 1998).

## 3 A numerical model of disk heating

We model disk heating in the sheared sheet, which approximates the local dynamics of a differentially rotating disk (Spitzer & Schwarzschild, 1953; Goldreich & Lynden-Bell, 1965; Julian & Toomre, 1966). The LSR is assumed to travel around the Galactic centre in a circular orbit of radius $`R_0`$ and angular frequency $`\mathrm{\Omega }>0`$. We use the same Cartesian coordinate system $`(x,y,z)`$ introduced in the last Section, restrict ourselves to the $`z=0`$ plane, and denote position and velocity by $`𝐱=(x,y)`$ and $`𝐯=(\dot{x},\dot{y})`$. For $`|x|,|y|R_0`$ the equations of motion of a test particle are

$`\ddot{x}2\mathrm{\Omega }\dot{y}4\mathrm{\Omega }Ax`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }}{x}},`$
$`\ddot{y}+2\mathrm{\Omega }\dot{x}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }}{y}},`$ (3)
where $`\mathrm{\Phi }(𝐱,t)`$ is the gravitational potential due to sources other than the axisymmetric Galactic disk, and $`A>0`$ is the Oort constant

$$A=\frac{1}{2}\left(R\frac{d\mathrm{\Omega }}{dR}\right)_{R_0},$$
(4)
where $`R\mathrm{\Omega }(R)`$ is the circular speed at radius $`R`$ and $`\mathrm{\Omega }(R_0)=\mathrm{\Omega }`$.

If $`\mathrm{\Phi }`$ = 0 the trajectories governed by equations (3) have two integrals of motion related to energy and angular momentum,

$$E\frac{1}{2}(\dot{x}^2+\dot{y}^24\mathrm{\Omega }Ax^2),H\dot{y}+2\mathrm{\Omega }x.$$
(5)
The solutions of the equations of motion are

$`x`$ $`=`$ $`x_g+a\mathrm{cos}(\kappa t+\varphi ),`$
$`y`$ $`=`$ $`y_g(t){\displaystyle \frac{2\mathrm{\Omega }}{\kappa }}a\mathrm{sin}(\kappa t+\varphi ),`$
$`y_g(t)`$ $`=`$ $`y_{g0}2Ax_gt,`$ (6)
where $`[x_g,y_g(t)]`$ are the coordinates of the guiding centre, $`a`$ is the epicycle amplitude, $`\varphi `$ is a phase constant, and $`\kappa `$ is the radial or epicycle frequency,

$$\kappa =\mathrm{\Omega }\left(4+2\frac{d\mathrm{ln}\mathrm{\Omega }}{d\mathrm{ln}R}\right)_{R_0}^{1/2}\text{or}\kappa ^2=4\mathrm{\Omega }(\mathrm{\Omega }A).$$
(7)
The energy and angular momentum are related to the guiding-centre radius and epicycle amplitude by

$$E=2A(A\mathrm{\Omega })x_g^2+\frac{1}{2}\kappa ^2a^2,H=2(\mathrm{\Omega }A)x_g,$$
(8)
and the guiding-centre radius is related to the phase-space coordinates by

$$x_g=\frac{\dot{y}+2\mathrm{\Omega }x}{2(\mathrm{\Omega }A)}.$$
(9)
The epicycle energy is defined as

$`E_x`$ $``$ $`\frac{1}{2}\left[\dot{x}^2+\kappa ^2(xx_g)^2\right]`$ (10)
$`=`$ $`E+{\displaystyle \frac{AH^2}{2(\mathrm{\Omega }A)}}`$
$`=`$ $`\frac{1}{2}\kappa ^2a^2`$
$`=`$ $`\frac{1}{2}\dot{x}^2+{\displaystyle \frac{2\mathrm{\Omega }^2}{\kappa ^2}}(\dot{y}+2Ax)^2.`$
The epicycle energy is closely related to the radial action

$$I=\frac{1}{2}\kappa a^2=E_x/\kappa .$$
(11)
For a particle in a circular orbit,

$$(x,y,\dot{x},\dot{y})=(x_g,y_{g0}2Ax_gt,0,2Ax_g),$$
(12)
thus $`E_x=0`$.

The approximations used in deriving the linearized equations of motion (3) are only marginally valid in the solar neighbourhood, since the epicycle amplitudes $`a`$ can be a significant fraction of $`R_0`$ (for a population with radial velocity dispersion $`\sigma _{xx}`$ in a galaxy with a flat rotation curve,

$$\overline{a^2}=\sigma _{xx}^2/\mathrm{\Omega }^2.$$
(13)
Thus old disk stars in the solar neighbourhood, with $`\sigma _{xx}40\mathrm{km}\mathrm{s}^1`$, have $`(\overline{a^2})^{1/2}0.18R_0`$). Nevertheless, we believe that the sheared sheet accurately captures the most important features of the evolution of disk-star kinematics.

The kinematics of a population of stars is described by the df $`f(𝐱,𝐯,t)`$, where $`fd𝐱d𝐯`$ is the number of stars in the interval $`(𝐱,𝐯)(𝐱+d𝐱,𝐯+d𝐯)`$. According to Jeans’s theorem (Binney & Tremaine, 1987), a stationary df can only depend on the integrals of motion $`E`$ and $`H`$ (or $`E_x`$, $`x_g`$, $`a`$, etc.). In the solar neighbourhood ($`x=y=0`$), the integrals are $`H=\dot{y}`$ and $`E_x=\frac{1}{2}\dot{x}^2+2\mathrm{\Omega }^2\dot{y}^2/\kappa ^2`$. Thus in a steady state the velocity distribution must be an even function of $`\dot{x}`$.

It is straightforward to show that the mean velocity and velocity-dispersion tensor of any stationary, spatially homogeneous df in the sheared sheet must satisfy the relations

$$\overline{v}_x=0,\overline{v}_y=2Ax,\sigma _{xy}^2=\sigma _{yx}^2=0,\frac{\sigma _{xx}^2}{\sigma _{yy}^2}=\frac{\mathrm{\Omega }}{\mathrm{\Omega }A}.$$
(14)
A useful model df for the sheared sheet is

$$f(𝐱,𝐯)\mathrm{exp}\left(E_x/\sigma _0^2\right),$$
(15)
where $`E_x`$ is defined in equation (10). At a given location, the df (15) leads to a two-dimensional version of the triaxial Gaussian velocity distribution (1), in which the mean velocity and velocity-dispersion tensor satisfy the relations (14). Thus the df (15) is also sometimes (confusingly) called the Schwarzschild df.

In our simulations the relations (14) are not satisfied exactly, so it is useful to introduce the principal-axis system $`(x_1,x_2)`$ in which the velocity-dispersion tensor is diagonal; the $`x_1`$-axis is chosen to lie within $`45\mathrm{°}`$ of the $`x`$-axis. The corresponding velocity dispersions are $`\sigma _1`$ and $`\sigma _2`$, which we will normally plot instead of $`\sigma _{xx}`$ and $`\sigma _{yy}`$; usually in our simulations and in the data we shall find that $`\sigma _1>\sigma _2`$. The vertex deviation $`l_v`$ is the Galactic longitude of the $`x_1`$-axis, and is given by

$$l_v\frac{1}{2}\mathrm{arctan}\left(\frac{2\sigma _{xy}^2}{\sigma _{xx}^2\sigma _{yy}^2}\right),$$
(16)
where $`|l_v|<45\mathrm{°}`$. Stellar populations in the solar neighbourhood have vertex deviations in the range $`0l_v30\mathrm{°}`$ (Binney & Merrifield, 1998).

Throughout the paper, we assume that the rotational curve of the underlying axisymmetric galaxy is flat, $`R\mathrm{\Omega }(R)=\text{constant}`$, so that $`A=\frac{1}{2}\mathrm{\Omega }`$ and $`\kappa =\sqrt{2}\mathrm{\Omega }`$. We shall assume that the surface density of the disk in the solar neighbourhood is $`\mathrm{\Sigma }_d=50M_{}\mathrm{pc}^2`$; recent observational estimates are $`\mathrm{\Sigma }_d=40M_{}\mathrm{pc}^2`$ (Crézé et al., 1998), $`48M_{}\mathrm{pc}^2`$ (Holmberg & Flynn, 2000), and $`42\pm 6M_{}\mathrm{pc}^2`$ (Korchagin et al., 2003). The assumed surface density only enters our analysis through the definition of the fractional spiral amplitude $`ϵ`$ below (eq. 21).

### 3.1 Spiral waves

We approximate spiral structure as a superposition of waves with surface density and potential

$`\mathrm{\Sigma }(x,y,t)`$ $`=`$ $`\mathrm{\Sigma }_s(t)\mathrm{exp}i(k_xx+k_yy),`$
$`\mathrm{\Phi }(x,y,t)`$ $`=`$ $`\mathrm{\Phi }_s(t)\mathrm{exp}i(k_xx+k_yy),`$ (17)
where $`𝐤=(k_x,k_y)`$ is the wavenumber, and only the real part of $`\mathrm{\Sigma }`$ or $`\mathrm{\Phi }`$ is physical. Without loss of generality we can assume that $`k_y0`$; then the spiral is trailing if $`k_x>0`$ and leading if $`k_x<0`$. In this paper, we only consider trailing spiral waves. The number of arms $`m`$ and the pitch angle $`\alpha `$ are determined by $`𝐤`$ through the relations

$$k_y=\frac{m}{R_0},\left|\frac{k_y}{k_x}\right|=\mathrm{tan}\alpha .$$
(18)
The relation between potential and surface density for spiral waves is given by (e.g. Binney & Tremaine 1987)

$$\mathrm{\Phi }_s(t)=\frac{2\pi G\mathrm{\Sigma }_s(t)}{k},$$
(19)
where $`k=|𝐤|`$.

Steady spirals, in which $`\mathrm{\Sigma }_s(t),\mathrm{\Phi }_s(t)\mathrm{exp}(i\omega t)`$, heat stars on nearly circular orbits only at the Lindblad resonances (Lynden-Bell & Kalnajs, 1972), which occur at the guiding centre radii given by

$$\omega =2Ak_yx_g\pm \kappa .$$
(20)
We focus instead on transient spiral structure, which can heat stars over a range of radii. We model each transient spiral using a Gaussian amplitude dependence centred at time $`t_s`$ with standard deviation $`\sigma _s`$:

$$\mathrm{\Phi }_s(t)=\frac{2\pi ϵG\mathrm{\Sigma }_d}{k}\mathrm{exp}\left[\frac{(tt_s)^2}{2\sigma _s^2}+i(\theta +2Ak_yx_ct)\right].$$
(21)
Here $`\theta `$ is a phase constant, and $`x_c`$ is the corotation radius relative to the LSR (the phase of the spiral is constant for an observer on a circular orbit at $`x_c`$). The parameter $`ϵ`$ measures the amplitude of the spiral, and the normalizing constants are chosen so that the maximum surface density in the spiral is a fraction $`ϵ`$ of the surface density of the underlying disk (eq. 19).

In each simulation the trajectories of the stars were followed between time $`t=0`$ and $`t=t_0`$. During this interval the stars were perturbed by $`N_s`$ transient spiral waves. Each wave had phase constant $`\theta `$ chosen randomly from $`[0,2\pi ]`$, corotation radius $`x_c`$ chosen randomly from a Gaussian distribution centred on the LSR with standard deviation $`\sigma _c`$, and central time $`t_s`$ chosen randomly from $`[2\sigma _s,t_0+2\sigma _s]`$ (this is slightly longer than the integration interval, to include the effects of transients whose wings are inside the integration interval although their central times are not).

The fractional amplitude of the surface-density perturbation due to a single transient spiral is $`ϵ`$ at the wave peak. However, in some ways a better quantity to compare with the observational data on spiral amplitudes is the root-mean-square (rms) time average of the fractional surface-density amplitude,

$$ϵ_{\mathrm{rms}}=ϵ\left(\frac{\pi ^{1/2}N_s\sigma _s}{t_0}\right)^{1/2}.$$
(22)
A closely related quantity is the rms potential perturbation,

$$\mathrm{\Phi }_{\mathrm{rms}}=2\pi Gϵ_{\mathrm{rms}}\mathrm{\Sigma }_dR_0\frac{|\mathrm{sin}\alpha |}{\sqrt{2}m};$$
(23)
(the factor $`\sqrt{2}`$ arises because $`\mathrm{\Phi }_{\mathrm{rms}}`$ is the rms fluctuation rather than the rms amplitude, which is smaller by $`\sqrt{2}`$). Jenkins & Binney (1990) use Fokker-Planck calculations of heating by transient spiral structure to estimate that $`\mathrm{\Phi }_{\mathrm{rms}}=(9\text{}13\mathrm{km}\mathrm{s}^1)^2`$.

The power spectrum of each transient is a Gaussian with standard deviation $`(2\sigma _s^2)^{1/2}`$. The central frequencies of the power spectra of the transients follow a Gaussian distribution with standard deviation $`2Ak_y\sigma _c`$. Thus the power spectrum of the ensemble of $`N_s`$ transients is smooth if the overlap factor

$$CN_s\frac{(2\sigma _s^2)^{1/2}}{2Ak_y\sigma _c}=\frac{N_s}{2^{3/2}Ak_y\sigma _c\sigma _s}$$
(24)
is large compared to unity; on the other hand if $`C1`$ the power spectrum and hence the heating is localized at narrow resonances. All of our simulations have $`C1`$ (see Table 1).

The dispersion of the power spectrum of all the transients is $`\sigma _t`$, given by

$$\sigma _t^2=\frac{1}{2\sigma _s^2}+(2Ak_y\sigma _c)^2.$$
(25)
Note that in this model problem different azimuthal wavenumbers $`m`$ and $`m^{}=fm`$ are equivalent if we rescale the other variables by: $`𝐤^{}=f𝐤`$, $`\alpha ^{}=\alpha `$, $`t^{}=t`$, $`(x^{},y^{})=(x/f,y/f)`$, $`\mathrm{\Phi }_s^{}=\mathrm{\Phi }_s/f^2`$, $`\mathrm{\Sigma }_s^{}=\mathrm{\Sigma }_s/f`$, $`ϵ^{}=ϵ/f`$, $`\sigma _{ij}^{}=\sigma _{ij}/f`$ $`\sigma _c^{}=\sigma _c/f`$, $`\sigma _0^{}=\sigma _0/f`$, $`\sigma _s^{}=\sigma _s`$.

### 3.2 Determining the solar neighbourhood velocity distribution

Assume that the disk stars are formed on circular orbits at time $`0`$. We wish to determine the velocity distribution at the present time $`t_0`$ in the solar neighbourhood, $`(x,y)=(0,0)`$. This is a two-point boundary-value problem rather than an initial-value problem: the boundary conditions on the phase-space coordinates are $`(x,y,\dot{x},\dot{y})=(x_g,y_{g0},0,2Ax_g)`$ at $`t=0`$ (on circular orbits) and $`(x,y)=(0,0)`$ at $`t=t_0`$ (in the solar neighbourhood). The solutions to the boundary-value problem will be a set of distinct initial positions $`[x_i(t=0),y_i(t=0)]`$ or final velocities $`[\dot{x}_i(t=t_0),\dot{y}_i(t=t_0)](u_i,v_i)`$, $`i=1,2,\mathrm{}`$ (the sign of $`u`$ is chosen to agree with the usual convention that stars moving towards the Galactic centre have $`u>0`$). Because the motion of stars in the fluctuating gravitational field of the transient spirals is complicated, we expect that there will be many solutions. In this idealized model the velocity-space df in the solar neighbourhood will consist of a set of delta-functions at the velocities $`(u_i,v_i)`$, but in practice these will be smeared into a continuous distribution by observational errors, the non-zero volume of the “solar neighbourhood”, the small initial velocity dispersion of the stars when they are formed, etc.

Normally, boundary-value problems are solved most efficiently by iterative methods. In this case, however, most of the solutions have only a tiny domain of attraction in the $`(u,v)`$ plane; thus iterative methods are inefficient. We have therefore adopted a different approach, based on Monte Carlo sampling (Dehnen, 2000).

Stars are not born on precisely circular orbits, in part because star-forming clouds are not on circular orbits. We may therefore assume that the stars initially have a Schwarzschild df (eq. 15), with a small but non-zero initial velocity dispersion $`\sigma _0`$.

Following Dehnen, we integrate orbits backward in time, starting at $`t=t_0`$ (which we call the initial integration time or “present epoch”) and ending at $`t=t_0`$ (the final integration time or “formation epoch”). The initial conditions are chosen at random from $`(x,y,\dot{x},\dot{y})=(0,0,u,v)`$, with $`|u|,|v|<v_{\mathrm{max}}`$. We choose $`v_{\mathrm{max}}=110\mathrm{km}\mathrm{s}^1`$, large enough to include almost all disk stars in the solar neighbourhood. We then integrate equations (3) backwards in time to the formation epoch. The collisionless Boltzmann equation (Binney & Tremaine, 1987) states that the phase-space density around a trajectory is time-invariant. Therefore the phase-space density $`F_i`$ at $`(u_i,v_i)`$ at the present epoch is given by (15), where the epicycle energy $`E_x`$ is measured at the formation epoch. We assume that the spatial volume of our survey of the solar neighbourhood is independent of velocity, so the density in velocity space at the present epoch is proportional to $`F_i`$. Dehnen’s procedure therefore provides a Monte Carlo sampling of the velocity-space df at the present epoch.

We convert our Monte Carlo realization to a smooth df by convolution with observational errors. We replace the point solutions $`𝐯_i=(u_i,v_i)`$ by Gaussians, that is,

$$f(𝐯)=\underset{i}{}F_i\delta (𝐯𝐯_i)$$
(26)
is replaced by

$$f(𝐯)=\frac{1}{2\pi \sigma _{\mathrm{ob}}^2}\underset{i}{}F_i\mathrm{exp}\left[\frac{(𝐯𝐯_i)^2}{2\sigma _{\mathrm{ob}}^2}\right],$$
(27)
where $`\sigma _{\mathrm{ob}}`$ is the assumed observational error.

We have described the procedure for estimating the df of stars of age $`t_0`$, using equation (15) and the epicycle energies at $`t=0`$. However, we can also determine the df for stars having an intermediate age $`t_0t_m`$ in the course of the same orbit integration, using the epicycle energies of the stars at the intermediate time $`t_m`$. This approach measures the dependence of the df on stellar age, and thus determines the avr for a given realization of the set of transient spiral arms, using only a single set of orbit integrations.

### 3.3 Units and model parameters

The time unit is chosen to be $`\mathrm{\Omega }^1`$ and the distance unit is chosen to be $`R_0`$. The azimuthal orbital period is then $`2\pi /\mathrm{\Omega }=2\pi `$. For conversion to physical units we assume that $`R_0=8\mathrm{kpc}`$ and $`\mathrm{\Omega }R_0=220\mathrm{km}\mathrm{s}^1`$, so that the unit of time is 35.6 Myr. Most of our numerical integrations ran for an interval $`t_0=281`$, corresponding to 40.8 rotation periods of the LSR or 10.0 Gyr.

We consider only trailing spirals with $`m=2`$ or $`m=4`$. The standard deviation of the distribution of corotation radii for the transient spirals was $`\sigma _c=0.25`$ corresponding to 2 kpc. With these choices the overlap factor is $`C=1.41N_s(m/2)/(\mathrm{\Omega }\sigma _s)`$.

The radial velocity dispersion of stars at the time of their birth is taken to be $`\sigma _0=3\mathrm{km}\mathrm{s}^1`$. The observational error in the velocities is taken to be $`\sigma _{\mathrm{ob}}=3\mathrm{km}\mathrm{s}^1`$; for the typical Hipparcos parallax error of 1 mas, $`\sigma _{\mathrm{ob}}`$ corresponds to the velocity error arising from the distance uncertainty for a star with transverse velocity $`30\mathrm{km}\mathrm{s}^1`$ at a distance of 100 pc.

In each simulation we launch either $`1\times 10^7`$ or $`7\times 10^7`$ stars from the solar neighbourhood and integrate back for 10 Gyr (the smaller simulations are used to estimate moments, such as those in Figure 1, and the larger simulations are used to estimate the df, as in Figure 7). In the larger simulations, typically 10000–15000 stars contribute significantly to the solution, in the sense that their epicycle energy at the formation epoch is $`<3\sigma _0^2`$. This number is larger than the number of stars in the Hipparcos samples to which we compare (e.g. 4600 in the largest color-selected sample analysed by Dehnen 1998). We have verified that the conclusions below are unaffected if we double the number of stars in the simulation or change the random-number seed used to select the sample of stars (keeping the same seed to select the spiral transients).

We carried out a number of simulations with different values of the spiral-wave parameters $`m`$ (number of arms), $`N_s`$ (number of spirals), $`ϵ`$ (fractional surface-density perturbation), $`\sigma _s`$ (duration of the transient), $`\sigma _c`$ (dispersion in corotation radius), and $`\alpha `$ (pitch angle). The parameters of the simulations are shown in Table 1. Each simulation is defined by these six parameters and the random number seeds for both the stars and the spiral arms.

Our runs have pitch angles of $`5\mathrm{°}`$, $`10\mathrm{°}`$, $`20\mathrm{°}`$, and $`40\mathrm{°}`$, which span the range of observed pitch angles in spiral galaxies. We then choose the amplitude $`ϵ`$ and the number of transients $`N_s`$ to approximately reproduce the observed radial velocity dispersion of old main-sequence stars in the solar neighbourhood<sup>1</sup><sup>1</sup>1More precisely, main-sequence stars redward of the Parenago discontinuity at $`BV=0.6`$ are all believed to have the same mean age, and have a color-independent radial dispersion $`38\mathrm{km}\mathrm{s}^1`$ (Dehnen & Binney, 1998). We match this to the radial dispersion in our simulation by assuming a constant star-formation rate, which is consistent with solar neighbourhood data within the large uncertainties (Binney, Dehnen & Bertelli, 2000)..

All of our runs have two-armed spirals ($`m=2`$), except for run R<sub>m4</sub>, which has the same parameters as run R<sub>20</sub> except that $`m=4`$.

We also have run R<sub>sim</sub>, in which we have used a larger number of weaker transients than in run R<sub>20</sub>, and have also adjusted some of the central times of the spiral transients to better reproduce the features of the observed solar neighbourhood df (see §4.4).

The rms potential fluctuation $`\mathrm{\Phi }_{\mathrm{rms}}=(13\text{}20\mathrm{km}\mathrm{s}^1)^2`$ in our four runs, higher than the velue $`\mathrm{\Phi }_{\mathrm{rms}}=(9\text{}13\mathrm{km}\mathrm{s}^1)^2`$ derived by Jenkins & Binney (1990) from Fokker-Planck calculations by about a factor of 2–2.4. A minor component of this difference arises because Jenkins & Binney (1990) include gmcs, which contribute about 20% of their total heating rate. A more important effect, pointed out to us by A. J. Kalnajs, is that in our model the heating is spatially inhomogeneous: our assumptions place the Sun fairly close to corotation for most of the transients, so that heating at the Lindblad resonances is more effective inside and outside the solar circle than it is at the solar circle.

The fractional surface-density fluctuation $`ϵ`$ required to reproduce the avr in run R<sub>5</sub> exceeds unity, and the rms surface-density fluctuation is close to unity; these high values are unrealistic and suggest that a larger number of weaker transients is required if the pich angle is as small as $`5\mathrm{°}`$.

Runs R<sub>10</sub>, R<sub>20</sub>, R<sub>40</sub>, and R<sub>m4</sub> have peak amplitudes $`ϵ=0.7`$$`0.8`$, similar to the strongest grand-design spirals, but these large amplitudes are only present for a small fraction of the time ($`N_s\sigma _s/t_0=0.03`$–0.06). The relatively small number of high-amplitude transients in these runs is consistent with a model in which the Galaxy has experienced several short-lived grand-design spiral phases, possibly caused by minor mergers. In contrast, run R<sub>sim</sub> has a much larger number of weaker transients ($`ϵ=0.3`$, $`N_s\sigma _s/t_0=0.17`$).

Each run was repeated seven times (realizations a-g), using different random-number seeds to generate both the stars and the spiral transients, in order to distinguish systematic from stochastic variations. To study extremely long-term behavior, we also ran a realization h, in which both the integration time and the number of transients are a factor of three larger than in a-g.

## 4 Results

### 4.1 Age-velocity dispersion relation

In §3.2, we have shown that a single set of orbit integrations can be used to derive the df for stars of all ages between 0 and $`t_0`$. Thus we may combine dfs of different ages to produce the df that corresponds to a uniform distribution of ages between 0 and 10 Gyr (hereafter the “combined df”), as would result from a uniform star-formation rate. The average moments of this df for our runs are shown in Table 2, along with the observed values taken from Dehnen (1998). Of course, the close agreement of the radial velocity dispersion $`\sigma _{xx}`$ in each run with the observed value of $`37\mathrm{km}\mathrm{s}^1`$ arises because the strength and number of spiral transients were chosen to produce the correct radial dispersion.

The evolution of the principal axes $`\sigma _1`$ and $`\sigma _2`$ of the velocity-dispersion tensor in runs R<sub>5</sub>, R<sub>10</sub>, R<sub>20</sub> and R<sub>40</sub> is shown in Figure 1. The three rows of panels show the results for three different realizations, i.e., three different sets of random-number seeds used to generate both the stars and the spiral transients. In general, the results in this and other figures are independent of the realization of the stellar distribution, and all of the variations seen are due to changes in the realization of the transient spiral structure—a reflection of the fact that the Monte Carlo simulation of the stellar distribution is a numerical method while the simulation of the distribution of spiral transients is a model of a stochastic physical process.

The avrs illustrate the following:

* Open spirals (large pitch angle) heat the disk stars more effectively than tightly wound spirals. For example, both the number $`N_s`$ and the surface-density amplitude $`ϵ`$ of the spiral transients are larger in R<sub>5</sub> than the other runs shown (see Table 1), even though all runs produce the same velocity dispersion by design.
* The heating occurs in discrete steps separated by periods of zero growth. This feature arises because we have only a small number of spiral transients over the course of the simulation, each of which has the same amplitude and lasts only for a characteristic time $`\sigma _s=1`$ out of the total integration time of $`281`$. During the intervals between the spiral transients, the df is stationary so the velocity dispersions are constant.
* In run R<sub>5</sub> the growth of the dispersion slows sharply after the radial dispersion reaches $`30\mathrm{km}\mathrm{s}^1`$. This effect probably arises because the epicycle amplitude becomes comparable to the radial spacing between the spiral arms, so that the forces from the spirals tend to average to zero over an epicycle oscillation. More specifically, the radial spacing of the spiral arms is $`\mathrm{\Delta }x=2\pi R_0\mathrm{tan}\alpha /m`$; the rms epicycle amplitude $`(\overline{a^2})^{1/2}`$ is given by equation (13), and we expect the spiral heating to become ineffective when the rms epicycle amplitude exceeds the peak-to-trough distance between arms, that is, when $`(\overline{a^2})^{1/2}/\mathrm{\Delta }x0.5`$. For small pitch angle $`\alpha `$,

$$\frac{(\overline{a^2})^{1/2}}{\mathrm{\Delta }x}=\frac{m\sigma _{xx}}{2\pi \mathrm{\Omega }R_0\mathrm{tan}\alpha }0.50\frac{m}{2}\frac{\sigma _1}{30\mathrm{km}\mathrm{s}^1}\frac{5\mathrm{°}}{\alpha }.$$
(28)
* There are substantial differences among the realizations of a given run; in other words there are significant stochastic effects in spiral-wave heating (see the discussion of the exponent $`p`$ below for further detail). The stochastic variations in the heating rate become smaller as the velocity dispersion increases (see also Figure 2, which follows the avrs over an even longer time interval).

For comparisons with theoretical models, it is useful to fit the simulated avrs with a parametrized function. A common parametrization (e.g. Lacey 1991) is

$$\sigma _1(t)=(\sigma _0^{1/p}+Ct)^p,$$
(29)
where $`\sigma _0`$ is the initial dispersion along the major principal axis, fixed at $`3\mathrm{km}\mathrm{s}^1`$ as in §3.3. The motivation for this form is the assumption that the evolution of the dispersion is governed by a diffusion equation with velocity-dependent diffusion coefficient,

$$\frac{d\sigma _1^2(t)}{dt}=C_1\sigma _1^q(t),$$
(30)
where $`C_1=2pC`$ and $`q=p^12`$, with $`C`$, $`p`$ and $`q`$ constants. For $`\sigma _1\sigma _0`$, equation (29) implies $`\sigma _1(t)t^p`$.

The dashed lines in Figure 1 show the results of least-squares fits for the parameters $`C`$ and $`p`$ in equation (29). In each case, the fitting function successfully preserves the overall shape of the avr but smooths over fluctuations due to the finite number of transients. We fit each of the seven realizations of each run and obtained the values of the avr exponent $`p`$ given in Table 3. The standard deviation in $`p`$ can be as large as 35% (in run R<sub>20</sub>) among realizations of a single parameter set, and the variation due to changes in the spiral-structure parameters is just as large. Thus, the exponent of the avr is not likely to discriminate well between heating by transient spiral structure and heating by other mechanisms.

Figure 2 shows the avrs of realization h, which was run for $`t_0=30`$ Gyr, three times as long as the other runs. The purpose of this long run was to see how well the parametrization (29) worked as the heating process continued. The dashed lines show that equation (29) still provides a reasonable fit; however, the best-fit exponent $`p`$ (see Table 3) is systematically smaller than in the shorter realizations a–g, except in run R<sub>5</sub>. Therefore, equation (29) should be considered as a fitting function valid over a limited time interval, rather than a formula with predictive power.

All of our runs have the same value of the parameter $`\sigma _s`$ (Table 1), which represents the duration of the spiral transients. We have experimented with a range of values of $`\sigma _s`$. As shown in Figure 3, the dispersion for the combined df with a uniform distribution of ages between 0 and 10 Gyr is almost independent of $`\sigma _s`$ so long as $`\sigma _s0.5`$. A probable explanation is that the width of the power spectrum of the transient perturbations, given by equation (25), is dominated by the dispersion in pattern speeds rather than the duration of a single transient once $`\sigma _s1/(2\sqrt{2}Ak_y\sigma _c)`$, which occurs when $`\sigma _s1.4`$ for our parameters. Of course, when $`\sigma _s`$ becomes very large, the overlap factor of equation (24), will decrease below unity and the heating will be localized at discrete resonances.

### 4.2 The velocity ellipsoid

The time evolution of the vertex deviation is shown in Figure 4. Not surprisingly, the stochastic variations in the vertex deviation are larger for runs R<sub>10</sub>, R<sub>20</sub> and R<sub>40</sub>, which have larger pitch angle and fewer transients, and the fluctuations in the vertex deviation decline with time as the growth rate of the velocity dispersion declines. The small but non-zero vertex deviations seen near the end of the runs are roughly consistent with the vertex deviation $`10\mathrm{°}`$ seen in the old stellar population in the solar neighbourhood. There is no evidence of a systematic preference for positive or negative vertex deviations for old stars.

Figure 5 shows the evolution of the ratio of the major and minor principal axes of the velocity ellipsoid. As in the case of the vertex deviation, the fluctuations are much larger in the runs with larger pitch angle. In a steady-state axisymmetric galaxy with a flat rotation curve, the axis ratio should be $`\sigma _1/\sigma _2=[\mathrm{\Omega }/(\mathrm{\Omega }A)]^{1/2}=2^{1/2}=1.414`$ (eq. 14), and the observed axis ratio gradually approaches this value as the runs progress.

Figure 6 shows the evolution of the mean radial and azimuthal velocity in the solar neighbourhood. There are significant fluctuations, up to $`10\mathrm{km}\mathrm{s}^1`$ for the larger pitch angles, even in the oldest stellar populations. These fluctuations limit the accuracy with which we can construct axisymmetric models of local Galactic kinematics.

We have also plotted the axis ratio of the velocity ellipsoid against the vertex deviation but have found no systematic correlation between these two quantities.

The parameters of the df for a population of stars with a uniform distribution of ages between 0 and 10 Gyr are shown in Table 2. They are all consistent with the observations, except for the mean azimuthal velocity $`\overline{v}`$. The discrepancy in $`\overline{v}`$ reflects the fact that our simulations do not show asymmetric drift because of the approximations inherent in the sheared sheet \[eqs. (3) with the perturbing potential $`\mathrm{\Phi }=0`$ are symmetric under reflection through the origin\].

### 4.3 Shape of the velocity distribution

Figure 7 shows the df for stars with age 5 Gyr, from runs R<sub>5</sub>a, R<sub>10</sub>a, R<sub>20</sub>a and R<sub>40</sub>a. We present the df in two forms: (i) we display all of the points from our sampling of velocity space at the present epoch that yield nearly circular orbits at the formation epoch (i.e. those with epicycle energy $`3\sigma _0^2`$, $`\sigma _0=3\mathrm{km}\mathrm{s}^1`$). (ii) We show contour plots of the smooth df that is expected with Hipparcos observational errors at $`100\mathrm{pc}`$ (see eq. 27), for two different realizations of the stellar distribution and the same realization of the spiral transients. The difference between these two realizations (the second and third rows) is an indication of the numerical noise in our simulations, which is small.

The simulated dfs in the second and third row of panels do not resemble the Gaussian or Schwarzschild df (eq. 1). They are much lumpier, and the lumps are distributed more-or-less uniformly over a large area in velocity space, rather than being concentrated at zero velocity. Qualitatively, the Schwarzschild dfs resemble a contour map of a single large mountain, whereas the simulated dfs resemble an entire mountain range. The “lumpiness” or degree of substructure in the df is a strong function of the pitch angle of the spiral arms, and in the range of pitch angles $`10\mathrm{°}`$$`20\mathrm{°}`$ that observations suggest for our Galaxy the dfs are highly irregular.

The dfs in Figure 7 exhibit prominent streaks in the direction $`v=v_y=`$constant, particularly at the larger pitch angles. The explanation of these streaks is straightforward, and dates back to Woolley (1961). It is straightforward to show from equations (6) that for stars in the solar neighbourhood ($`x=y=0`$) the radial and tangential velocities $`(u,v)`$ are related to the guiding-center coordinates $`(x_g,y_g)`$ by

$$(x_g,y_g)=\frac{1}{2(\mathrm{\Omega }A)}(v,u).$$
(31)
Suppose that scattering from a spiral transient leads to an enhancement in the density of stars over the range of guiding-centre radii $`x_gx_g+\mathrm{\Delta }x_g`$. Since the guiding-centre velocity is $`\dot{y}_g=2Ax_g`$ (eq. 6), the stars will slowly spread out into an extended streamer in physical space; after time $`t`$ the streamer will have width of order the epicycle amplitude and length $`\mathrm{\Delta }y2A\mathrm{\Delta }x_gt`$. The stars found at a given location along the streamer (for example the solar neighbourhood) will thus have very similar values of $`x_g`$ and thus of $`v`$. In other words, a tidal streamer will present itself in the local velocity distribution as a streak with constant $`v`$.

The 1:1 relation between $`(u,v)`$ and $`(y_g,x_g)`$ implied by equation (31) implies that the structure seen in Figure 7 is closely related to the guiding center correlations discussed by Toomre & Kalnajs (1991).

Skuljan, Hearnshaw & Cottrell (1999) find that the observed df in the solar neighbourhood exhibits parallel “branches” of enhanced density in velocity space. It is tempting to identify these with the streaks predicted by Woolley (1961) and seen in Figure 7. However, the branches found by Skuljan, Hearnshaw & Cottrell (1999) are tilted with respect to the expected orientation $`v=`$constant; they find instead that $`v=\text{constant}0.47u`$. Thus the correct interpretation of these branches is unclear: they may arise from some physical effect we have not modeled, or they may be artifacts arising from our predisposition to find order in random patterns. Note, for example, that parallel tilted “branches” appear to be present in the contour plots of the df in run R<sub>20</sub>a shown in Figure 7, with a slope about half of that claimed by Skuljan et al.

The principal conclusion from Figure 7 is that the lumps in the simulated dfs, at least in the simulations with pitch angles $`10\mathrm{°}`$, bear a striking resemblance to the moving groups in the solar neighbourhood (see for example Dehnen 1998, Skuljan, Hearnshaw & Cottrell 1999, or Chereul et al. 1998, 1999), which strongly suggests that at least the older moving groups are a consequence of the same transient spiral structure that heats the stellar distribution. Thus we suggest that old moving groups arise from fluctuations in the gravitational potential rather than fluctuations in the star-formation rate, as has usually been assumed. This result was anticipated to some extent by Kalnajs (1991) and Dehnen (1999, 2000), who argued that some of the most prominent substructure in the solar neighbourhood df arose from the outer Lindblad resonance with the Galactic bar. The distinction is that they focused on a single steady bar-like potential whereas we have studied stochastic spiral transients.

### 4.4 Evolution of moving groups

One of the striking features of the solar neighbourhood df is that the same moving groups appear to be present in stellar populations of different ages (see Table 2 of Dehnen 1998 or the top row of Figure 9). This result is difficult to understand in the traditional view that moving groups arise from inhomogeneous star formation. Thus it is worthwhile to explore the temporal evolution of the df and the survival of moving groups in our simulations.

Our standard four simulations, R<sub>5</sub>, R<sub>10</sub>, R<sub>20</sub>, and R<sub>40</sub>, are not ideal for this task. As we have seen, run R<sub>5</sub> has relatively weak moving groups and unrealistic parameters. The other runs show a strong moving group at the origin in the combined df which is not present in the observed df. This moving group arises because the typical interval between spiral transients in these simulations is 0.5–1 Gyr (Table 1). Thus in most of these runs there has been no spiral transient in the last $`0.5`$ Gyr, and all of the stars formed since the last transient still have their birth epicycle energy, which is nearly zero. Moreover, these runs feature a relatively small number of strong spirals, which hit the disk very infrequently but hard. They may therefore overestimate the persistence of moving groups.

To address these concerns, we have run another simulation R<sub>sim</sub>. In this run, we increase the number of spirals to $`N_s=48`$ and set $`ϵ=0.31`$. We have also adjusted the central times of the spiral transients as follows: (1) We have set the central time $`t_s`$ of the most recent transient to be the present time $`t_0`$; (2) Because Dehnen’s subsample B1, in which the oldest stars have age 0.4 Gyr, already shows two prominent moving groups (top left panel of Figure 9), we have set the central time of another transient to $`t_00.2`$ Gyr; (3) The other central times are chosen randomly between $`t_s=0`$ and $`t_s=t_00.2`$ Gyr. This strategy suppresses the strong moving group at the origin and makes the other features in the df more prominent.

Figure 8 shows the velocity distribution in run R<sub>sim</sub>a as a function of stellar age. “Moving groups” (i.e. the same local concentration of stars in velocity space) can be seen among stars of any age, even the oldest stars. Generally with increasing age there are more but weaker moving groups, however, a strong “moving group” can be present even in a very old sample (see panel 17, 8.5 Gyr). Notice that the same “moving group” can be present in stars with a wide range of ages; for example, the moving group “B” in panel 13 can be traced from panels 11 to 14, a range of at least 1.5 Gyr in age. This result is consistent with the observation that the same moving group is seen in observational sub-samples of different age (e.g., Dehnen 1998, Table 2).

In Figure 9 we show the observed and simulated dfs, divided into age bins. The observed df, shown in the top row, is taken from Dehnen (1998). The panels represent subsamples determined by color cuts, which should correspond to age cuts. The maximum ages in the three subsamples B1, B2, and B3 (the first, second, and fourth panels) are 0.4 Gyr, 2 Gyr, and 8 Gyr (Table 1 of Dehnen 1998). The sample B4GI (fifth panel) is a combination of red ($`BV>0.6`$) main-sequence stars and giant stars, which should have a uniform distribution of ages if the star-formation rate has been uniform. The third panel, labeled B23, is obtained by subtracting the scaled df of B2 from the df of B3, in order to represent approximately the df of stars with ages 2–8 Gyr.

In the bottom two rows of Figure 9 we show the simulated dfs from two runs: run R<sub>sim</sub>, and a modified version of R<sub>sim</sub> in which the width $`\sigma _0`$ of the epicycle energy at formation (eq. 15) is increased from $`3\mathrm{km}\mathrm{s}^1`$ to $`8\mathrm{km}\mathrm{s}^1`$. The different panels show the dfs corresponding to the same age ranges as in the subsamples of the data. The simulated dfs reproduce most of the prominent features of the observed dfs: (1) The stars clump together in moving groups. (2) The number of moving groups increases as the mean stellar age increases, from left to right in the figure. (3) The same moving groups are present in stars with a wide range of ages. (4) The moving groups appear to be elongated along lines of constant $`v`$, although this effect is significantly stronger in the simulations than the observations.

Increasing the initial epicycle energy (compare the second and third row of Figure 9) spreads out the moving groups somewhat, and enhances their elongation along lines of constant $`v`$, but they remain prominent; the moments of the df given in Table 2 are almost unchanged.

In Figure 10 we plot the birthplaces of groups “A” and “B” shown in panels 2 and 12 of Figure 8. The stars in the younger group, A (age 1 Gyr), come from a relatively small arc, while the stars in the older moving group B (age 6 Gyr) were born with a wide range of Galactic azimuths. In contrast, Eggen’s model for moving groups predicts that the stars in a group should be born in a small region and thus should have the same age and metallicity.

We also ran a modified version of R<sub>20</sub>, which has the same parameters as R<sub>20</sub> but adjusts the central times of the transients as described above for R<sub>sim</sub>. This run gives more prominent moving groups than R<sub>sim</sub>, and the groups last for a longer time (up to 3.5 Gyr). This result suggests that a small number of strong spiral transients tend to produce strong moving groups and a lumpy df, while a large number of weak transients tend to produce weak or undetectable moving groups and a smooth df. Therefore, the observed moving groups suggest that at most a few tens of spiral transients have affected the stars in the solar neighbourhood.

We also looked for moving groups in run R<sub>m4</sub>, which has the same parameters (including the random number seeds) as the modified run R<sub>20</sub>, except that the spirals have four arms<sup>3</sup><sup>3</sup>3It is striking that run R<sub>m4</sub> yields almost the same velocity dispersion as run R<sub>20</sub> (Table 2). The reason is as follows: consider a simulation with azimuthal wavenumber $`m=2`$ that yields a power-law age-velocity dispersion relation of the form (29), $`\sigma (t)=(Ct)^p`$. We expect that the constant $`C`$ is proportional to the diffusion coefficient, and that this in turn is proportional to the mean-square fluctuating force $`F^2(k\mathrm{\Phi }_s)^2`$; thus $`C=cF^2`$ where $`c`$ is a constant and the age-velocity dispersion relation can be written as $`\sigma (t)=F^{2p}(ct)^p`$. According to the scaling relations at the end of §3.1, a simulation with azimuthal wavenumber $`m`$ is equivalent to one with $`m^{}=fm`$ (here $`f=2`$) if $`\sigma ^{}=\sigma /f`$ and $`F^{}=F/f`$. Thus $`\sigma ^{}(t)=F^{2p}(ct)^p/f=f^{2p1}(F^{})^{2p}(ct)^p`$. Since the exponent $`p0.5`$ for run R<sub>20</sub>, $`f^{2p1}1`$ and we expect the age-velocity dispersion relations for runs R<sub>20</sub> and R<sub>m4</sub> to be nearly the same.. We found that substructure in the df in run R<sub>m4</sub> was much weaker than in the modified run R<sub>20</sub>, presumably because the wavelength $`2\pi /k`$ of the arms is smaller so the df is more thoroughly mixed for a given amount of heating. The absence of substructure suggests that transient spirals in the Galaxy may be predominantly two-armed rather than four-armed. Further exploration of the effects of the spiral-structure parameters on the structure of the local df is obviously worthwhile.

We have also looked for substructure in the velocity df in runs in which the duration $`\sigma _s`$ of the spiral transients was varied (Figure 3). We found that the strength of the substructure declines as $`\sigma _s`$ increases; in other words, spiral transients with duration much longer than $`\mathrm{\Omega }^1`$ tend to heat the disk without generating moving groups.

We have also investigated heating by gmcs rather than spiral structure, and find that this process does not yield significant substructure in the df (De Simone, 2000). Thus the observed strong substructure suggests that spiral structure rather than gmcs is the dominant heating mechanism in the solar neighbourhood, a conclusion reached by others from independent arguments (see §1).

## 5 Inhomogeneous star formation

So far in this paper we have assumed that the star-formation rate is uniform in both space and time, and that both the avr and the structure of the velocity distribution in the solar neighbourhood are due to scattering by the large-scale potential fluctuations associated with transient spiral arms. The traditional view, which we explore briefly in this Section, is that the structure in the local velocity distribution results from spatial and temporal variations in the star-formation rate. Of course, some elements of both pictures must be correct, since we know both that spiral arms are present and that star formation is concentrated in gmcs.

The observed avr requires heating by some mechanism—whether spiral arms, gmcs, or one of the other mechanisms discussed in §1. However, in this Section we ignore the details of this heating, to focus on the properties of the irregularities in the velocity distribution that are due to inhomogeneous star formation alone. We consider the following crude model:

1. All stars form in clusters of $`N`$ stars, of limited spatial extent.
2. The cluster dissociates after the stars form. When the cluster dissociates, the df of its stars is Gaussian in both position and velocity, with dispersions $`\sigma _q`$ and $`\sigma _p`$ respectively. The stars do not interact with one another once the cluster has dissolved; thus, their motion is governed by equations (6). Typical values of $`\sigma _q`$ and $`\sigma _p`$ are a few pc or $`\mathrm{km}\mathrm{s}^1`$.
3. The heating mechanism does not disrupt or spread the phase-space stream of stars from the dissolved cluster (any such effects will make many of the conclusions below even stronger). We ignore the details of the heating mechanism, and account for heating by simply assuming that the dissociated cluster is born at position $`(x_0,y_0)`$, on a non-circular orbit with velocities $`u_0`$ and $`v_0`$. Thus, the initial df of the dissociated cluster is

$$f_0(𝐳)=\frac{N}{(2\pi \sigma _p\sigma _q)^2}\mathrm{exp}\left[\frac{1}{2}(𝐳𝐳_0)𝐐(𝐳^T𝐳_0^T)\right].$$
(32)
where $`𝐳=(x,y,u,v)`$, $`T`$ denotes the transpose, and

$$𝐐=\left(\begin{array}{cccc}\sigma _q^2& 0& 0& 0\\ 0& \sigma _q^2& 0& 0\\ 0& 0& \sigma _p^2& 0\\ 0& 0& 0& \sigma _p^2\end{array}\right).$$
(33)
The covariance matrix is

$$(𝐳^T𝐳_0^T)(𝐳𝐳_0)=𝐐^1,$$
(34)
where the average $``$ is over the df $`f_0(𝐳)`$.

The trajectories (6) can be rewritten in the form (Asiain, Figueras & Torra, 1999)

$$𝐳(t)=𝐀(t)𝐳(0),$$
(35)
where $`t`$ is the time elapsed since the dissociation of the cluster, and (recall that $`u=\dot{x}`$)

$$𝐀(t)=\left(\begin{array}{cccc}\frac{\mathrm{\Omega }A\mathrm{cos}\tau }{\mathrm{\Omega }A}& 0& \frac{\mathrm{sin}\tau }{\kappa }& \frac{1\mathrm{cos}\tau }{2(\mathrm{\Omega }A)}\\ \frac{2A\mathrm{\Omega }(\mathrm{sin}\tau \tau )}{\kappa (\mathrm{\Omega }A)}& 1& \frac{2\mathrm{\Omega }(1\mathrm{cos}\tau )}{\kappa ^2}& \frac{\mathrm{\Omega }\mathrm{sin}\tau A\tau }{\kappa (\mathrm{\Omega }A)}\\ \frac{A\kappa \mathrm{sin}\tau }{\mathrm{\Omega }A}& 0& \mathrm{cos}\tau & \frac{\kappa \mathrm{sin}\tau }{2(\mathrm{\Omega }A)}\\ \frac{2A\mathrm{\Omega }(\mathrm{cos}\tau 1)}{\mathrm{\Omega }A}& 0& \frac{2\mathrm{\Omega }\mathrm{sin}\tau }{\kappa }& \frac{\mathrm{\Omega }\mathrm{cos}\tau A}{\mathrm{\Omega }A}\end{array}\right).$$
(36)
Here $`\tau \kappa t`$. The matrix $`𝐀(t)`$ has the following properties: $`𝐀(0)=𝐈`$ (by definition), $`|𝐀(t)|=1`$ (by Liouville’s theorem), and $`𝐀^1(t)=𝐀(t)`$ (by time-reversibility).

Since phase-space density is conserved along the particle trajectories, the initial df is Gaussian (eq. 32), and $`𝐳(t)`$ is a linear function of $`𝐳(0)`$, the df of the dissociated cluster will remain Gaussian for all time. The centroid of the Gaussian will be $`𝐳_c(t)=𝐳_0𝐀^T(t)`$ and the covariance matrix will be

$$𝐂(t)[𝐳^T𝐳_c^T(t)][𝐳𝐳_c(t)]=𝐀(t)𝐐^1𝐀^T(t).$$
(37)
The disrupted cluster is greatly elongated in the azimuthal direction. Therefore to estimate the maximum number of stars from a cluster that could be visible in the solar neighbourhood, it is useful to compute the linear number density at time $`t`$, $`n(y,t)=f(𝐳,t)𝑑x𝑑𝐯`$, where $`𝐳=(x,y,u,v)=(𝐱,𝐯)`$. Since the df is Gaussian, the linear density must also be Gaussian,

$$n(y,t)=\frac{N}{(2\pi \sigma _y^2)^{1/2}}\mathrm{exp}\left[(yy_c)^2/2\sigma _y^2\right],$$
(38)
where

$$\sigma _y^2=[yy_c(t)]^2=C_{22}(t);$$
(39)
at large times, $`\kappa t1`$, we have (Asiain, Figueras & Torra, 1999)

$$\sigma _y\frac{4A\mathrm{\Omega }t}{\kappa ^2}(\sigma _p^2+4\mathrm{\Omega }^2\sigma _q^2)^{1/2}.$$
(40)
For a typical star-forming cluster, $`\mathrm{\Omega }\sigma _q\sigma _p`$ ($`\mathrm{\Omega }\sigma _q=0.03\mathrm{km}\mathrm{s}^1(\sigma _q/1\mathrm{pc})`$ compared to typical values $`\sigma _p=1\mathrm{km}\mathrm{s}^1`$). Thus we may drop the term involving $`\sigma _q`$, so the maximum linear density from a disrupted cluster of age $`t`$ is

$$n_{\mathrm{max}}(t)=\frac{N}{(2\pi )^{1/2}\sigma _pt}\left(\frac{\kappa ^2}{4A\mathrm{\Omega }}\right),$$
(41)
where the factor in brackets is unity for a flat rotation curve. In the solar neighbourhood,

$$n_{\mathrm{max}}(t)=4\times 10^5\mathrm{pc}^1N\left(\frac{1\mathrm{km}\mathrm{s}^1}{\sigma _p}\right)\left(\frac{10\mathrm{Gyr}}{t}\right).$$
(42)
With the same assumptions, the rms spatial extent of the cluster in the radial direction is

$$\sigma _x=\frac{\sigma _p}{\mathrm{\Omega }}\left[\frac{1}{2}\mathrm{sin}^2\kappa t+(1\mathrm{cos}\kappa t)^2\right]^{1/2};$$
(43)
averaging over the epicycle period we have $`\sigma _x=(\frac{7}{4})^{1/2}\sigma _p/\mathrm{\Omega }=50\mathrm{pc}(\sigma _p/1\mathrm{km}\mathrm{s}^1)`$.

Now suppose that our survey includes all stars in the solar neighbourhood to a distance $`L`$. Then if $`L\sigma _x`$ we can detect up to $`2n_{\mathrm{max}}L`$ stars in the moving group, while if $`L\sigma _x`$ we will detect up to $`(\pi /2)^{1/2}n_{\mathrm{max}}L^2/\sigma _x`$ stars in the group (these upper limits are achieved if the Sun lies exactly on the centroid $`x_c(t)`$). Extrapolating these formulas until they meet, we find that the upper limit to the number of detected stars in the moving group is given by

$$n_{\mathrm{mg}}<8\times 10^3N\frac{L}{100\mathrm{pc}}\frac{1\mathrm{km}\mathrm{s}^1}{\sigma _p}\frac{10\mathrm{Gyr}}{t}$$
(44)
for $`L>80\mathrm{pc}(\sigma _p/1\mathrm{km}\mathrm{s}^1)`$, and

$$n_{\mathrm{mg}}<1.0\times 10^2N\left(\frac{L}{100\mathrm{pc}}\right)^2\left(\frac{1\mathrm{km}\mathrm{s}^1}{\sigma _p}\right)^2\frac{10\mathrm{Gyr}}{t}$$
(45)
for $`L<80\mathrm{pc}(\sigma _p/1\mathrm{km}\mathrm{s}^1)`$. The sample of Dehnen (1998) was limited to distances $`L100\mathrm{pc}`$ (the median parallax error in the Hipparcos catalog is $`\sigma _\pi =1`$ mas and Dehnen’s sample is restricted to stars with $`\sigma _\pi /\pi 0.1`$). If we assume that at least $`n_{\mathrm{mg}}=50`$ stars are needed to define a moving group, then any such group identified in Dehnen’s catalog must have come from a cluster with initial population

$$N6\times 10^3\frac{\sigma _p}{1\mathrm{km}\mathrm{s}^1}\frac{t}{10\mathrm{Gyr}}\text{if }\sigma _p<1.3\mathrm{km}\mathrm{s}^1,$$
(46)
or

$$N5\times 10^3\left(\frac{\sigma _p}{1\mathrm{km}\mathrm{s}^1}\right)^2\frac{t}{10\mathrm{Gyr}}\text{if }\sigma _p>1.3\mathrm{km}\mathrm{s}^1.$$
(47)
These are conservative lower limits, since (i) we have neglected the spreading of the stream due to interactions with small-scale fluctuations in the gravitational potential, such as gmcs; (ii) we have assumed that the radial coordinate of the stream of stars is centred precisely on the Sun.

Infrared observations show that most star formation occurs in embedded star clusters within gmcs. These clusters have masses in the range $`50`$$`10^3M_{}`$ (Lada & Lada, 2003) and virial dispersions $`\sigma _p1\mathrm{km}\mathrm{s}^1`$ and thus are too small to contribute detectable moving groups with an age exceeding $`1\mathrm{Gyr}`$. An alternative is to suppose that the dissociated cluster is the gmc itself. In this case $`\sigma _p5\mathrm{km}\mathrm{s}^1`$ which requires $`N1.3\times 10^5`$; given that the largest gmcs have gas masses $`10^5`$$`10^6M_{}`$ we would require extremely high star-formation efficiency to produce a sufficient number of stars. In other words it is difficult to find a plausible astrophysical source for the star-forming regions which are supposed to be the progenitors of the old moving groups.

## 6 Migration

Sellwood & Binney (2002) discuss the dynamics of radial migration due to transient spiral arms, and the relationship between radial migration and disk heating. They estimate that old stars formed in the solar neighbourhood should now be scattered nearly uniformly throughout the annulus from $`R_0\pm 4\mathrm{kpc}`$.

We have examined the radial migration or diffusion of stars that is induced by spiral-wave scattering, using the spiral-arm parameters for our runs R<sub>5</sub>, R<sub>10</sub>, R<sub>20</sub> and R<sub>40</sub>. In each case we started 100,000 test stars in circular orbits near the solar neighbourhood ($`\pm 0.2\mathrm{kpc}`$) and followed them for 10 Gyr. In Figure 11 we plot the rms change in guiding-centre radius and the distribution of stars as functions of final epicycle energy. The present in-plane epicycle energy of the Sun is only $`E_x=64(\mathrm{km}\mathrm{s}^1)^2`$, which is nearly zero on the scale of these plots. If these simulations correctly represent the properties of the local spiral structure, we expect that the radial migration of the Sun since its formation 4.5 Gyr ago is in the range $`0.45\pm 0.25\mathrm{kpc}`$ (pitch angle $`10\mathrm{°}`$) to $`1.4\pm 0.9\mathrm{kpc}`$ (pitch angle $`20\mathrm{°}`$). The sheared sheet is symmetric in $`\pm x`$ so our simulations cannot address the question of whether the migration is likely to be inward or outward; Sellwood & Binney (2002) argue that migration is generally outward because of the surface-density gradient in the disk. These estimates of the migration are consistent with the estimate made by Wielen, Fuchs & Dettbarn (1996), on the basis of the Sun’s metallicity and the radial metallicity gradient in the Galactic disk, that the Sun has migrated outward by $`1.9\pm 0.9\mathrm{kpc}`$ in the past 4.5 Gyr.

## 7 Discussion

We have explored the hypothesis that transient spiral arms are the dominant mechanism that drives the evolution of the velocity distribution in the solar neighbourhood (Barbanis & Woltjer, 1967; Sellwood & Carlberg, 1984; Carlberg & Sellwood, 1985). The most thorough investigations of this mechanism so far have been based on the Fokker-Planck equation (Jenkins & Binney, 1990; Jenkins, 1992). Our investigation is based on direct numerical integration of test-particle orbits in the sheared sheet (eqs. 3). We impose a randomly distributed sequence of trailing $`m=2,4`$ spiral waves with pitch angle $`\alpha `$ and a strength that is Gaussian in time (eq. 21), and apply the boundary conditions that the test particles must be on circular orbits at their formation time and in the solar neighbourhood at the present time.

We confirm that transient spiral waves can heat the galactic disk in velocity space. The configurations of the spiral waves in our simulations, such as the number ($`N_s=9`$–48), duration ($`\mathrm{\Omega }\sigma _s=0.5`$–2.5), rms fluctuation amplitude ($`ϵ_{\mathrm{rms}}=0.16`$–0.86), and pitch angle ($`\alpha =5\mathrm{°}`$$`40\mathrm{°}`$) are all plausible. Spiral arms can lead to a wide range in the observed exponent of the age-velocity dispersion relation (eq. 29), $`p=0.2`$–0.76, so this exponent is not a strong discriminator between different mechanisms of disk heating. The variations of the vertex deviation, mean radial velocity, and axis ratios of the velocity ellipsoid in the solar neighbourhood are consistent with observations.

The small-scale structure of the local stellar velocity distribution cannot be interpreted in terms of axisymmetric, steady-state models. Our simulated distribution functions do not have the Schwarzschild form (1) and in addition exhibit rich small-scale structure, similar to the “moving groups” and “branches” that observers describe in the solar neighbourhood distribution function. This result suggests that moving groups arise primarily from smooth star formation in a lumpy potential, in contrast to the traditional assumption that they reflect lumpy star formation in a smooth potential—although of course elements of both pictures must be present to some extent.

All these findings strongly support the conclusion that transient spiral waves are the dominant heating mechanism of the disk in the solar neighbourhood.

Strong substructure in the local velocity distribution is most easily produced when the Galaxy experiences a small number of short-lived but strong (fractional surface-density amplitude $`ϵ1`$) spiral transients, as might be caused by minor mergers. A model in which spiral transients are continuously present at a lower amplitude could reproduce the age-velocity dispersion relation without significant substructure.

We also find some other results:

1. The age-velocity dispersion relation depends strongly on the stochastic properties of the spiral waves, especially in cases where the heating is due mainly to a relatively small number of spiral transients (e.g. $`10`$ in runs R<sub>20</sub> and R<sub>40</sub>). Therefore, it is difficult, if not impossible, to infer the details of the heating mechanism from observations of the age-velocity dispersion relation, no matter how accurate.
2. Spiral arms lead to radial migration of stars (Sellwood & Binney 2002 and references therein). For the spiral-wave parameters we used, the rms radial migration of a star like the Sun is a strong function of pitch angle $`\alpha `$, ranging from $`0.45\pm 0.25\mathrm{kpc}`$ for $`\alpha =10\mathrm{°}`$ to $`1.4\pm 0.9\mathrm{kpc}`$ for $`\alpha =20\mathrm{°}`$. In other words the relation between the amount of heating of the solar neighborhood and the amount of radial migration is not unique, a conclusion already reached by Sellwood & Binney (2002) using other arguments.
3. Stars in old moving groups did not form at a common place and time, and therefore do not necessarily share a common age and metallicity.
4. The traditional model, in which structure in the velocity distribution function is the result of inhomogeneous star formation, is difficult to reconcile with at least two observations: 1. In this model, all of the stars in a moving group should have the same age and metallicity, while in fact the stars in prominent moving groups show a wide range of ages. 2. It is difficult to find a suitable astrophysical source that is small enough, cold enough, and produces enough stars to be the progenitor of the old moving groups.
5. We have also investigated heating by giant molecular clouds alone, and find that this process produces much weaker substructure in the distribution function than is observed (De Simone, 2000); thus giant molecular clouds are unlikely to be the dominant heating mechanism in the solar neighbourhood.

Although heating by spiral waves explains most of the properties of the solar neighbourhood velocity distribution, there is still room for other heating mechanisms, such as the Galactic bar, giant molecular clouds, and halo substructure or minor mergers.

## 8 Acknowledgments

We thank James Binney, Walter Dehnen, and Jerry Sellwood for discussions and advice, and the referee, Agris Kalnajs, for many insightful comments. This research was supported in part by NSF Grant AST-9900316.


format: text

text:
# Far Infrared mapping of NGC 891Based on observations with the Infrared Space Observatory (ISO), an ESA project with instruments funded by ESA member states (especially the PI countries: France, Germany, the Netherlands, and the United Kingdom) and with the participation of ISAS and NASA.

## 1 Introduction

NGC 891 is one of the most extensively observed and studied edge-on spiral galaxy in the nearby universe. At a distance of 9.5 Mpc (van der Kruit & Searle 1981), it has been classified as an Sb galaxy by Sandage (1961) and is often quoted to be very similar to our own Galaxy. However it has also some unusual characteristics. For example NGC 891 contains one of the most spectacular layers of extraplanar diffuse ionised gas (DIG) (Dettmar 1990; Rand et al. 1990; Keppel et al. 1991; Pildis et al. 1994; Rand 1997, 1998; Hoopes et al. 1999; Howk & Savage 2000, Otte et al. 2001). The galaxy also has a radio continuum halo (Allen et al. 1978; Hummel et al. 1991), an HI halo (Swaters et al. 1997) and an X-ray halo (Bregman & Pildis 1994; Bregman & Houck 1997).

The disk of NGC 891 is a strong source of both CO and 13CO emission (Sofue et al. 1987, García-Burillo et al. 1992) and the global distribution and kinematics of molecular gas have been investigated by Handa et al. (1992), Scoville et al. (1993), García-Burillo & Guélin (1995), Sakamoto et al. (1997). The first extragalactic direct detection of large-scale molecular hydrogen was established in the disk of NGC 891 (Valentijn & van der Werf 1999), based on observations with the Short-Wavelength Spectrometer (SWS) instrument aboard the Infrared Space Observatory (ISO). The spectrum of the unidentified infrared (UIR) emission bands between 5.9 and 11.7 $`\mu `$m has been also observed for the first time in the disk of an external galaxy in NGC 891 (Mattila et al. 1999), using the low-resolution spectrometer of the ISOPHOT instrument aboard ISO. The mapping of the Unidentified Infrared Bands emitted by NGC 891 was also done with the ISOCAM instrument on board ISO (Le Coupanec et al. 1999).

The distribution of the cold dust in NGC 891 has been observed at submillimeter(submm)/millimeter(mm) wavelengths using the IRAM 30 m telescope (Guélin et al. 1993) and the Submillimeter Common-User Bolometer Array (SCUBA) at the James Clerk Maxwell Telescope (JCMT)(Alton et al. 1998; Israel et al. 1999). The deep SCUBA maps revealed dust emission over 2/3 of the optical disk, but did not shed light on the existence of an extraplanar dust emission (Alton et al. 2000). Large amount of cold dust ($`15`$ K) was found in the disk (Alton et al. 1998) - an order of magnitude more than the amount of warm dust detected by Wainscoat et al. (1987) using the Infrared Astronomical Satellite (IRAS).

However, the peak of the spectral energy distribution (SED) from cold dust, lying at the long far-infrared (FIR) wavelengths, is outside the wavelength coverage of IRAS, and thus has not been readily accessible until recently. With a wavelength coverage extending to 240 $`\mu `$m and a superior intrinsic sensitivity as compared to IRAS, the ISOPHOT instrument (Lemke et al. 1996) on board ISO (Kessler et al. 1996) was the first to directly measure the peak of the FIR SED for a number of nearby galaxies (Tuffs et al. 1996 for NGC 6946, Haas et al. 1998 for M 31, Wilke et al. 2003 for the SMC, Hippelein et al. 2003 for M 33), for smaller statistical samples (Krügel et al. 1998; Siebenmorgen et al. 1999; Contursi et al. 2001), for the ISOPHOT Virgo Cluster Deep Sample (Tuffs et al. 2002a,b; Popescu et al. 2002; Popescu & Tuffs 2002a) and for the sample of Bright Revised Shapley Ames galaxies (Bendo et al. 2002, 2003).

Here we used the dedicated mapping mode P32 (Tuffs & Gabriel 2003) of ISOPHOT to obtain deep maps of NGC 891 at 170 and 200 $`\mu `$m wavelengths. These data, taken at the peak of the SED from cold dust, provide a precise measurement of the stellar light re-radiated by grains in NGC 891. They also allow cold dust to be probed in regions of low surface brightness, such as the disk periphery and halo of NGC 891, regions which are currently inaccessible to submm facilities. When combined with the higher column density presented by the edge-on orientation, the ISOPHOT observations of NGC 891 have the capability of tracing cold dust to higher galactic radii. Thus, the goals of this study are: 1) to derive flux densities at longer FIR wavelengths, and thereby to directly measure the peak of the FIR SED in NGC 891; 2) to measure the physical extent of the dust disk and in particular to search deep for a cold dust counterpart to the extended HI disk; 3) to search for cold extraplanar FIR emission; 4) to compare the brightness and colour profile of the dust disk with predictions from the three dimensional model of stellar and dust distributions in NGC 891 proposed by Popescu et al. (2000a). This paper concentrates on the 1st, 2nd, and 4th goals, while the search for a dust counterpart to the extended HI disk is presented in Popescu & Tuffs (2003). In Sect. 2 we present the observations and data reduction, including the derivation of the integrated flux densities from the maps. In Sect. 3 we give a detailed comparison between the data and the model predictions for NGC 891. Some implications of this comparison are discussed in Sect. 4. A summary is given in Sect. 5.

## 2 Observations and data reduction

The observations were made using the ISOPHOT-C200 2$`\times 2`$ pixel array in the C160 and C200 filters111We warn the reader that the name of the C160 filter does not correspond to the central wavelength, unlike the case of the C200 filter., which respectively cover passbands of $`130\mathrm{\hspace{0.17em}218}`$ and $`170\mathrm{\hspace{0.17em}239}\mu \mathrm{m}`$ and have central wavelengths of 170 and 200 $`\mu `$m. The “P32” mapping mode was used to provide near Nyquist sampling over a large area encompassing the optical galaxy and the extended HI disk, as well as the surrounding background. The HI disk extends $`10`$ arcmin from the nucleus in the southern half of the galaxy and $`7`$ arcmin in the northern half (Swaters et al. 1997). In order to cover the whole southern field, and also to have a symmetrical map, a field of radius $`\pm 13.5`$ arcmin ($`\pm 40`$ kpc) was mapped along the major axis of the galaxy.

A total of 19,000 seconds were spent in three overlapping fields: north, south and central. Due to scheduling constraints, the coverage of the northern field was made in two shorter overlapping observations at each wavelength. Thus, 4 maps were made at each wavelength. The parameters of each of these observations are summarised in Table 1. Part of these data sets have been presented by Popescu & Tuffs (2002b), Tuffs & Popescu (2003) and by Dupac et al. (2003).

The data for each observation were separately processed using the latest P32 reduction package (Tuffs & Gabriel 2003), which corrects for the transient response of the detector pixels. This allowed high dynamic range maps to be constructed to levels of 1 percent of the peak disk brightness. A time-dependent flat-field correction was made for each map, by fitting a cubic function to the response of the detector pixels to the background. Calibration was made using V8.1 of the PHOT Interactive Analysis (PIA) Package (Gabriel et al. 1997). We emphasise that although the maps are oversampled, independent data contribute to each map pixel.

Finally, the maps were combined into a mosaic spanning half a degree parallel to the plane of the galaxy. To achieve a common photometric scale over all the fields, the northern and southern fields were first scaled to the background brightness found from colour-corrected COBE/DIRBE maps within a 1.5 degree radius circle of NGC 891, taking into account the temporal variation of the zodiacal light component of the background. Values of the observed ISO backgrounds and the COBE/DIRBE backgrounds to which they were scaled are given in Table 1. Scaling factors between the central field and the southern and northern fields were then calculated from source structure in the overlapping regions of the maps, and applied to the data of the central field as a linear interpolation in position between the overlapping regions. Finally, backgrounds were subtracted from the southern and northern fields, and a constant background from the central field such that the absolute rms discrepancy in the overlap regions was minimised.

Overall, the self consistency of the pixel responses to the backgrounds suggests a systematic uncertainty of 5 and 10 percent in the C160 and C200 filters, respectively. The random noise variations in the southern and northern regions of the mosaic are 0.1 and 0.2 MJy/sterad, respectively, some 1000 times fainter than the peak emission. Thus the maps are dynamic-range limited by the sidelobes of the beam response to the bright disk emission.

The derived mosaic maps at 170 and 200 $`\mu `$m are shown in Fig. 1. The central field is dominated by the FIR emission from the optical disk of the galaxy. Fainter emission is seen out to at least the edge of the optical disk (360′′ from the nucleus). A detailed investigation of the FIR emission beyond the optical disk is presented in Popescu & Tuffs (2003). Another feature seen in the maps is a pair of wings (highlighted by the outer contours in Fig 1) separated from the nucleus by $`\pm 320`$ arcsec in the scan direction. The feature is more prominent at 170 $`\mu `$m than at 200 $`\mu `$m. The alignment of this feature with the scan direction points to an instrumental origin. However a specific cause could not be identified and it cannot be completely ruled out that extraplanar emission has not been detected.

The integrated fluxes are given in Table 2, together with their uncertainties ($`ϵ(F)`$), peak brightness ($`B_{\mathrm{peak}}`$) and the corresponding uncertainties ($`ϵ(B_{\mathrm{peak}})`$) for both C160 and C200 maps. The last column represents an estimate of the systematic uncertainty in the detector responsivity derived from scatter of the response of the individual detector pixels to the background.

## 3 Comparison between the data and the model for the FIR emission of NGC 891

Especially because it is so well observed, but also because of its edge-on geometry, which allows the study of the vertical structure of the different stellar, dust and gas components, NGC 891 was chosen as a prototype since the study of van der Kruit & Searle (1981), Bahcall & Kylafis (1985) and Kylafis & Bahcall (1987). More recently, Xilouris et al. (1998, 1999) fitted the optical and NIR images of NGC 891 with simulated maps produced by radiative transfer calculations, in an attempt to derive the scale lengths and heights of the stellar and dust distributions. All these studies concentrated on modelling the optical appearance of NGC 891.

Popescu et al. (2000a) developed a modelling technique for the whole spectral energy distribution, from the UV to FIR and sub-mm, which was tested and also first applied to NGC 891. This model successfully reproduced the observed IRAS and sub-mm flux densities of NGC 891. Furthermore, it is the only model which makes direct predictions for the spatial distribution of the FIR emission. Before comparing the predictions of this model with the new ISOPHOT data, we review the main characteristics of the dust and stellar distributions in NGC 891, as derived from this model.

Both a diffuse component and a clumpy component associated with the star forming regions are calculated. The model for the diffuse component includes a consistent treatment of grain heating and emission, solves the radiation transfer problem for a finite disk and bulge, and self-consistently calculates the stochastic heating of grains placed in the resulting radiation field. Compared with the other self-consistent treatments of the UV/optical–FIR/submm SEDs our model is unique in that it analyses the surface brightness distribution of optical and NIR images to constrain the intrinsic distributions of the old stellar populations and associated dust (see also Misiriotis et al. 2001). The appearance of the optical/NIR images was produced in an iterative optimisation procedure using the technique for solving the radiation transfer equation for direct and multiply scattered light for arbitrary geometries by Kylafis & Bachall (1987). The calculation was done independently for a number of optical and NIR images, thus determining the extinction law for diffuse dust empirically.

The “young” stellar population was then determined primarily from the FIR/submm data. The emissivity of the “young stellar disk” was parameterised in terms of the current star formation rate ($`SFR`$). A second dust disk of grain mass $`M_{\mathrm{dust}}`$ was associated with the young stellar population. For the clumpy component a third primary parameter, $`F`$, was introduced to denote the fraction of non-ionising UV which is locally absorbed in HII regions around the massive stars. These three parameters -$`SFR`$, $`M_{\mathrm{dust}}`$ and $`F`$ \- are sufficient to fully determine the FIR-submm SED, due to the precise constraints on the distribution of stellar emissivity in the optical-NIR and on the distribution and opacity of dust in the “old dust disk” yielded by the radiation transfer analysis of the highly resolved optical-NIR images, coupled with the simple treatment of the young stellar population and associated dust.

Comparison with IRAS and submm data yielded a best solution for NGC 891 with $`SFR=3.8M_{}`$/yr, $`F=0.22`$ and $`M_{\mathrm{dust}}=7\times 10^7M_{}`$. This corresponds to a disk central face-on opacity in the V band of $`\tau _V^f=3.1`$ and a non-ionising UV luminosity of $`8.2\times 10^{36}`$ W. The luminosity of the diffuse dust emission is $`4.07\times 10^{36}`$ W, which accounts for 69$`\%`$ of the total FIR luminosity, and the luminosity of dust in the clumpy component is $`1.82\times 10^{36}`$ W, making up the remaining 31$`\%`$ of the FIR luminosity.

### 3.1 The integrated FIR flux densities

The predicted FIR SED of NGC 891 is presented in Fig. 2, together with the IRAS, SCUBA (Alton et al. 1998), IRAM (Guélin et al. 1993) and our ISOPHOT data. Both the model and the data were integrated only within $`\pm 225^{\prime \prime }`$, in order to match the SCUBA and the IRAS flux densities derived by Alton et al. (1998). At 170 and 200 $`\mu `$m the integrated flux densities within $`\pm 225^{\prime \prime }`$ are 166 and 144 Jy, respectively. Our model SED, derived and checked on IRAS and submm data by Popescu et al. (2000a), is found to be in excellent agreement with the ISOPHOT flux densities as well.

### 3.2 The FIR maps

A more stringent test of the model is to compare its predictions for the morphology of the dust emission with the observed morphology. For this purpose simulated FIR maps were produced using the actual pointing data to scan the diffuse disk model. The model map was then convolved with empirical PSFs derived from point source measurements. The comparison between the observed maps at 170 and 200 $`\mu `$m (Fig. 3a and 3b, respectively) and the simulated maps at the same wavelengths (Fig. 3c and 3d, respectively) show a remarkable agreement. To search for small differences between the model and the observations, not detectable in the maps due to the high dynamical range of the displayed data, we present in Fig. 3e and 3f the residual maps of the difference between the observed and the simulated maps at 170 and 200 $`\mu `$m, respectively.

At 170 $`\mu `$m, the main feature in the residuals is a localised, unresolved source in the northern side of the disk, with a peak of 52.3 MJy/sr. This localised source is probably a giant molecular cloud complex - associated with one of the spiral arms, and not considered in the simulated map, which only includes the diffuse component of the model. At this FIR wavelength our model predicts an $`11\%`$ contribution from the star-forming complexes. The integration of the unresolved source gives a flux density of 13.6 Jy, which is $`7\%`$ of the total flux density. Furthermore a faint extended source is seen in the southern side of the galaxy, of 9.5 Jy integrated flux density. This makes another $`5\%`$ of the total integrated emission. Thus the faint localised sources seen in the residual maps sum up to $`12\%`$ of the total FIR emission. This is in reasonable agreement with the prediction of our model, which reassures us that the template used in our model and scaled according to our model parameters ($`SFR`$ and the $`F`$ factor), is indeed a good representation for the galaxy. Apart from the two sources, a faint extended halo (at $`1\%`$ brightness level) is seen extending at large heights perpendicular to the disk.

At 200 $`\mu `$m, the residual map (Fig. 3f) shows the same localised sources present at 170 $`\mu `$. This time the northern localised source has a flux density of 8.4 Jy (5.1$`\%`$ of the total emission) and the southern localised source has a flux density of 6.9 Jy (4.2$`\%`$). Overall the percentage contribution of both localised sources to the total flux density is smaller at 200 $`\mu `$m than it is at 170 $`\mu `$m, which means that the localised sources are warmer than the diffuse component. This is to be expected if we consider that the localised sources are associated with the HII regions, and thus consistent with the model of Popescu et al. (2000a).

### 3.3 The FIR averaged profiles

Another way to compare the predictions of the model with the data is to look at the averaged profiles, both in the radial and vertical direction. Firstly, profiles were produced from the observed maps and from the maps of the predicted diffuse emission (Fig. 3a,b,c,d), by integrating the flux parallel to the minor axis of the galaxy for each bin along the major axis. The resulting profiles, which we refer to as “radial profiles”, are given in Fig. 4a and 4b in the 170 and 200 $`\mu `$m bands, respectively. Again, a good agreement between the data points and the model prediction is apparent at both wavelengths. Centred at +100 arcsec there is an excess emission in the observed profile with respect to the predicted one, which is due to the localised source in the northern half of the galaxy, already identified in the residual maps (Fig. 3e,f). This asymmetry is also seen in the 850 $`\mu `$m map of Alton et al. (1998) and in the 1300 $`\mu `$m map of Guélin et al. (1993). Likewise, the fainter localised source in the southern half of the galaxy manifests itself through the excess observed emission seen around -150 arcsec radius. At both wavelengths there is some faint emission apparent in the observed profile extending beyond the edge of the optical disk (360 arcsec). This emission is considered in Popescu & Tuffs (2003).

Similar to the construction of the radial profiles, profiles perpendicular to the major axis of the galaxy were calculated by integrating the flux parallel to the major axis of the galaxy for each bin along the minor axis. The resulting profiles, which we refer to as “vertical profiles”, are given in Fig. 4c,d. As expected, the predicted disk emission is unresolved by ISOPHOT, and therefore the predicted profile coincides with the beam profile. The observed profile is also unresolved up to $`\pm 100`$ arcsec, which corresponds to a factor of 10 in dynamic range. At larger heights above the disk there is some faint diffuse emission (at a level of a few percent of the peak brightness) in excess of the model predictions for the disk. We have already discussed in Sect. 2 that this emission may be of an instrumental nature, although a detection of extraplanar emission cannot be discounted completely. Such an emission may be expected to arise from dust grains carried out by a gentle wind, as described in Popescu et al. (2000b).

The measured and predicted radial colour profiles 200/170 are plotted in Fig. 4e. The predicted profile has a smooth progression towards colder emission with increasing radial distance. The measured profile broadly follows this prediction, although there may be evidence for a tendency for a flatter profile within $`\pm 200`$. The larger error bars on the data points beyond +200 arcsec are due to the lower brightness of the 170 $`\mu `$m outer disk on the northern side.

## 4 Discussion

At the wavelength of the ISOPHOT measurements presented here the model for NGC 891 predicts that the bulk of the FIR dust emission is from the diffuse component. The close agreement between the data and the model predictions, both in integrated flux densities, but especially in terms of the spatial distribution, constitutes a strong evidence that the large scale distribution of stellar emissivity and dust predicted by the model is in fact a good representation of NGC 891. In turn, this supports the prediction of the model that the dust emission in NGC 891 is predominantly powered by UV photons.

Depending on the FIR/submm wavelength, the UV powered dust emission arises in different proportions from within the clumpy component and from the diffuse component. For example at 60 $`\mu `$m, 61$`\%`$ of the FIR emission is powered by UV photons locally absorbed in star forming complexes, 19$`\%`$ by diffuse UV photons in the weak radiation fields in the outer disk (where stochastic emission predominates), and 20$`\%`$ by diffuse optical photons in high energy densities in the inner part of the disk and bulge. At 100 $`\mu `$m there are approximately equal contributions from the diffuse UV, diffuse optical and locally absorbed UV photons. At 170, 200 $`\mu `$m and submm wavelengths, most of the dust emission in NGC 891 is powered by the diffuse UV photons. So our analysis does not support the preconception that the weakly heated cold dust (including the dust emitting near the peak of the SED sampled by the ISOPHOT measurements presented here) should be predominantly powered by optical rather than UV photons. The reason is as follows: the coldest grains are those which are in weaker radiation fields, either in the outer optically thin regions of the disk, or because they are shielded from radiation by optical depth effects. In the first situation the absorption probabilities of photons are controlled by the optical properties of the grains, so the UV photons will dominate the heating. The second situation arises for dust associated with the young stellar population, where the UV emissivity far exceeds the optical emissivity.

## 5 Summary

Here we have presented deep FIR maps of the edge-on spiral galaxy NGC 891, obtained with the ISOPHOT instrument on board ISO. The observations were done near the peak of the FIR SED, namely at 170 and 200 $`\mu `$m. Both the integrated flux densities and the surface brightness distributions obtained from our ISOPHOT maps were found to be in excellent agreement with the predictions of the model for the optical/FIR/submm SED of NGC 891 of Popescu et al. (2000a). Furthermore, the model distribution of emitters and absorbers in NGC 891 is able to reproduce both the observed optical/NIR images and the FIR maps. NGC 891 is the first galaxy for which an intrinsic distribution of dust and stars was derived under such strong observational constraints.

###### Acknowledgements.

We would like to thank the anonymous referee for his useful comments and suggestions.