text:
# The OGLE-II event sc5_2859 – An example of disk-disk microlensing

## 1 Introduction

It is nearly ten years since the first gravitational microlensing event was detected toward the Galactic bulge (Udalski et al. 1993). Since then, microlensing has proved to be a useful tool for many astrophysical applications (for a review of Galactic microlensing, see Paczyński 1996). In the past, however, studies of microlensing statistics have been limited by the small number of detected events. This situation is improving as new microlensing projects begin operation, e.g., the OGLE-III111http://www.astrouw.edu.pl/~ogle/ogle3/ews/ews.html project, which is hoping to detect up to 1000 per year.

Two ways that microlensing statistics can be used to investigate Galactic models are through the study of the observed microlensing optical depth and the distribution of event time-scales. The optical depth toward the Galactic bulge has been studied by many collaborations (for example, Alcock et al. 2000 and Sumi et al. 2002) and the resulting values are significantly larger than predicted estimates (see, for example, Binney, Bissantz & Gerhard 2000; Evans & Belokurov 2002; Klypin, Zhao & Somerville 2002), although recent work has questioned the significance of this discrepancy (Popowski 2002; Afonso et al. 2003).

The distribution of event time-scales, which can be used to investigate the mass spectrum of lensing objects (e.g. Han & Gould 1996; Peale 1998), also seems to disagree with predicted estimates; it has long been suspected that microlensing searches appear to be identifying an unexpectedly large proportion of long-duration events. An important early study into the distribution of microlensing time-scales was carried out by Han & Gould (1996; see also Zhao, Rich & Spergel 1996). They found that in one-year of microlensing observations from the OGLE and MACHO collaborations 8% of events had time-scale greater than 70 days, whereas the largest fraction predicted by their theoretical models was 2%. Bennett et al. (2002a) analysed a more recent distribution of events timescales from MACHO data and found a similar excess of long duration events. Although no systematic analysis of event time-scales has been performed for the full 4-year OGLE-II catalogue, a preliminary investigation using the first 3-years of data also indicates the existence of an excess (Udalski et al. 2000).

To date, two microlensing events have been identified with time-scales greater than one year (OGLE-1999-BUL-32/MACHO-99-BLG-22, Mao et al. 2002, Bennett et al. 2002b; OGLE-1999-BUL-19, Smith et al. 2002). We will introduce a third event, sc5\_2859, in this paper. It has been proposed that some of these extreme long-duration events could be due to lensing by massive stellar remnants, e.g. black holes (Agol et al. 2002; Bennett et al. 2002a), although it is unlikely that this hypothesis could account for every long-duration event, for example OGLE-1999-BUL-19 (Smith et al. 2002).

In this paper we present the analysis of event sc5\_2859. We begin by describing the observational data (Section 2), before proceeding to fit the event with both the standard and parallax microlensing models (Section 3). We then investigate whether any constraints can be placed on the lens mass; we first attempt to utilise finite source size effects (Section 4.1), before considering more general arguments based on the relative transverse velocity of the lens and limits on the lens brightness (Section 4.2). In Section 4.1 we also discuss the possible location and spectral type of the source star. Section 5 contains a brief discussion regarding the nature of event sc5\_2859 and the possibility that there are an excess of long duration microlensing events in the OGLE-II catalogue. This section also mentions possible scientific returns from such long duration events. We conclude with a summary (Section 6).

## 2 Observational data

The event sc5\_2859 was identified during the second phase of the OGLE project (in which it was named BUL\_SC5 244353; Udalski et al. 2000) toward the Galactic bulge. The observations were carried out with the 1.3 m Warsaw telescope at the Las Campanas Observatory, Chile, which is operated by the Carnegie Institution of Washington. The instrumentation of the telescope and CCD camera are described in detail by Udalski, Kubiak, & Szymański(1997). The position of the source is RA=17:50:36.09 and Dec=$``$30:01:46.6 (J2000), corresponding to Galactic coordinates $`l=359.6235^{}`$ and $`b=1.4930^{}`$. This event is also included in the catalogue of Woźniak et al. (2001), which employs the Difference Image Analysis method of data reduction (Alard & Lupton 1998). This method generally appears to result in greater accuracy compared to the classical point spread function approach (Schechter, Mateo & Saha 1993), and so we have chosen to use the Difference Image Analysis data in the following analysis222 This analysis uses the full 4-year Difference Image Analysis data, generously provided for our use by the OGLE collaboration. The partial 3-year Difference Image Analysis data, along with the basic calibration data for this event, is publicly available online from http://astro.princeton.edu/~wozniak/dia/lens.. Unfortunately, this event was only detected after the peak in magnification had already occurred, which means that less than half of the full light curve is available for analysis.

The baseline magnitude for the source is given by $`I=18.53\pm 0.02`$ mag333 It should be noted that there is a discrepancy in the baseline magnitude of the source between the two photometric methods: the full 4-year Difference Image Analysis data have a baseline of $`I=18.53\pm 0.02`$ mag, but the 3-year classical point spread function data have a baseline of $`I=18.86\pm 0.02`$ mag. However, the conclusions regarding the nature of event sc5\_2859 remain unchanged if the classical point spread function data is used., and the colour is $`(VI)=1.97\pm 0.02`$ mag. In Fig. 1 we present the colour-magnitude diagram for the field around sc5\_2859. From this figure it can be seen that the field around sc5\_2859 is subject to a large amount of interstellar extinction, owing to the fact that the source is located very close to the Galactic plane. It appears that the lensed star is located in the main sequence branch of the colour-magnitude diagram. We will return to the issue of the possible source location and spectral type in Section 4.1.

## 3 Model fitting

We initially fit the event with the standard microlensing model, which assumes that the observer, lens and (point) source all move with constant velocities. The magnification is given by (see, for example, Paczyński 1986),

$$A(t)=\frac{u^2+2}{u\sqrt{u^2+4}},u(t)\sqrt{u_0^2+\tau (t)^2},$$
(1)
where $`u_0`$ is the impact parameter (in units of the Einstein radius) and,

$$\tau (t)=\frac{tt_0}{t_\mathrm{E}},$$
(2)
with $`t_0`$ being the time of the closest approach (i.e., maximum, magnification), and $`t_\mathrm{E}`$ the event time-scale. The time-scale is defined such that,

$$t_\mathrm{E}=\frac{r_\mathrm{E}}{v}=\frac{\stackrel{~}{r}_\mathrm{E}}{\stackrel{~}{v}},$$
(3)
where $`r_\mathrm{E}`$ is the lens’ Einstein radius, $`v`$ is the lens’ transverse velocity relative to the observer-source line of sight, and $`\stackrel{~}{r}_\mathrm{E}`$ and $`\stackrel{~}{v}`$ are the values of these two quantities projected onto the observer plane. Therefore, $`t_\mathrm{E}`$ corresponds to the time it takes for the lens to move a distance equal to its Einstein radius444An alternate definition for the time-scale is sometimes employed (particularly by the MACHO collaboration), which uses $`\widehat{t}`$, the Einstein diameter crossing time, i.e. $`\widehat{t}=2t_\mathrm{E}`$..

The Einstein radius projected onto the observer plane is related to the mass according to the following equation,

$$\stackrel{~}{r}_\mathrm{E}=\sqrt{\frac{4GMD_\mathrm{S}x}{c^2(1x)}},$$
(4)
where $`M`$ is the lens mass, $`D_\mathrm{S}`$ the distance to the source and $`x=D_\mathrm{L}/D_\mathrm{S}`$ is the ratio of the distance to the lens and the distance to the source. This shows the well-known degeneracy inherent in the standard microlensing formalism; the quantities $`\stackrel{~}{v}`$, $`M`$ and $`x`$ cannot be determined uniquely from a given microlensing light curve, even if the source distance is known.

The Difference Image Analysis flux is given by,

$$F(t)=F_0\left(f_\mathrm{S}\left[A(t)1\right]+1\right)F_{\mathrm{Ref}},$$
(5)
where $`F_0`$ is the total baseline flux from the source plus any blended star(s), if present, $`f_\mathrm{S}`$ is the ratio of the baseline source flux to the total baseline flux (i.e. a measure of the blending), and $`F_{\mathrm{Ref}}=416.15`$ is the flux of the reference image. All the fluxes here are in units of 10 ADU and can be converted into $`I`$-band magnitudes using the following transformation (given in Woźniak et al. 2001),

$$I(t)=I_{\mathrm{Ref}}2.5\mathrm{log}_{10}\frac{F(t)}{F_{\mathrm{Ref}}},$$
(6)
where $`I_{\mathrm{Ref}}=16.76`$ is the $`I`$-band magnitude of the reference image. Note that the reference image is brighter than the baseline magnitude for sc5\_2859.

We fit the light curve with the above five parameters for the standard model, i.e., $`t_0`$, $`t_\mathrm{E}`$, $`u_0`$, $`F_0`$, $`f_\mathrm{S}`$. The best-fit parameters are given in Table 1, and the corresponding light curve is shown in Fig. 2. Clearly this model is unable to provide a suitable fit for the data, with the best standard $`\chi ^2`$ per degree of freedom greater than 5.

Since the duration of this event is longer than one year, the next logical step is to fit the light curve with a model that incorporates the parallax effect (Gould 1992). This effect, which arises when the Earth’s motion around the Sun is considered, is described in detail in Soszyński et al. (2001); see also Alcock et al. (1995), Dominik (1998). It requires two additional parameters, the Einstein radius projected onto the observer plane, $`\stackrel{~}{r}_\mathrm{E}`$, and an angle in the ecliptic plane, $`\psi `$, describing the orientation of the lens trajectory (given by the angle between the heliocentric ecliptic $`x`$-axis and the normal to the trajectory).

The parallax model produces a drastic reduction in $`\chi ^2`$. The best-fit parameters are given in Table 1, and the corresponding light curve is shown in Fig. 2. The $`\chi ^2`$ per degree of freedom is 1.23, which indicates that the parallax model provides a reasonable fit to the data.

We also fit this event with a slight variation on the above parallax model. In the above model we describe the geometrical properties in the ecliptic plane and then project these quantities into the lens plane. However, since the ecliptic plane intersects the Galactic bulge, this can lead to the projection being almost singular. Therefore, to avoid this potential singularity, one can instead take the more conventional approach and describe the geometrical properties in the plane perpendicular to the line-of-sight (see, for example, Dominik 1998, Alcock et al. 1995); to describe the lens plane coordinate system we form a right-handed set with the $`x`$-axis chosen to correspond to the North Ecliptic Pole projected onto the lens plane and the $`z`$-axis chosen to be the observer-source line-of-sight, which implies that the $`y`$-axis corresponds to the intersection of the lens plane and the Ecliptic plane. We find that the best-fit parameters are practically identical and therefore we do not present them here. However, for this coordinate system, the angle describing the relative lens-source trajectory is found to be $`\theta =35.96_{1.20}^{+0.66}`$ degrees, where $`\theta `$ is the angle between the trajectory and the $`x`$-axis (measured from the positive $`x`$-axis towards the positive $`y`$-axis).

From Table 1, the parallax parameters that provide information regarding the lens properties are,

$$\stackrel{~}{r}_\mathrm{E}=10.84_{0.41}^{+0.44}\mathrm{au},t_\mathrm{E}=547.6_{7.8}^{+22.6}\mathrm{days},$$
(7)
which implies that the transverse velocity of the lens relative to the source, projected onto the observer plane is,

$$\stackrel{~}{v}=\frac{\stackrel{~}{r}_\mathrm{E}}{t_\mathrm{E}}=34.2_{1.9}^{+1.5}\mathrm{kms}^1.$$
(8)
If we convert the direction of $`\stackrel{~}{v}`$ from the lens plane ($`\theta =36.0`$ degrees) into the plane perpendicular to the line-of-sight to the Galactic centre, we find that $`\stackrel{~}{v}`$ is directed almost parallel to the Galactic plane in the direction of rotation ($`\theta _\mathrm{G}=6.1_{1.2}^{+0.7}`$ degrees, where $`\theta _\mathrm{G}`$ is measured from the Galactic plane towards the North Galactic Pole). In this aspect sc5\_2859 is similar to the strong parallax events presented in Bennett et al. (2002a), all of which had $`90^{}<\theta _\mathrm{G}<0^{}`$.

This timescale of 548 days is the second longest time-scale ever identified (after the event OGLE-1999-BUL-32/MACHO-99-BLG-22; see Mao et al. 2001, Bennett et al. 2002b). In addition, the value of $`\stackrel{~}{r}_\mathrm{E}`$ is also unusually large, which implies that this event could be another black-hole microlens candidate (c.f. the three current black-hole microlens candidates, which have $`\stackrel{~}{r}_\mathrm{E}=29.3,\mathrm{\hspace{0.17em}11.2}\mathrm{and}8.7\mathrm{au}`$; see Agol et al. 2002). Using equation (4), the mass of the lens for this event is given by,

$$M=7.21M_{}\left(\frac{\stackrel{~}{r}_\mathrm{E}}{10.84\mathrm{au}}\right)^2\left(\frac{\pi _{\mathrm{rel}}}{0.5\mathrm{mas}}\right),\pi _{\mathrm{rel}}\frac{\mathrm{AU}}{D_\mathrm{L}}\frac{\mathrm{AU}}{D_\mathrm{S}}.$$
(9)
A value of $`\pi _{\mathrm{rel}}=0.5`$ mas corresponds to a disk source with $`D_\mathrm{S}=2\mathrm{kpc}`$ and $`D_\mathrm{L}/D_\mathrm{S}=0.5`$. There are various approaches that can be employed to constrain the location of the source and lens, and hence the lens mass; these are considered in Section 4.

To verify the viability of our parallax fit, we proceed to fit the event with a model that incorporates a constant acceleration term, instead of the Earth’s centripetal acceleration (see Smith, Mao & Paczyński 2003). This model is unable to provide a feasible fit for sc5\_2859, meaning that the parallactic nature of the deviations appears to be secure.

## 4 Constraints on the lens mass

### 4.1 Finite-source effects

Since the peak magnification is predicted to be greater than 40, one may suspect that this event could be affected by finite-source effects (Gould 1994; Witt & Mao 1994; Nemiroff & Wickramasinghe 1994). Such effects become apparent when the lens passes sufficiently close to the source, resulting in an invalidation of the assumption that the source is point-like. However, there is one significant drawback for sc5\_2859, namely that there is no coverage for the peak of the light-curve, i.e. the point at which the lens and source are in closest alignment and therefore where the finite-source effects should be most prominent.

To implement the finite-source model requires an additional parameter, $`\rho _{}`$, which denotes the source radius in units of the lens’ angular Einstein radius. Despite the lack of coverage around the peak, constraints can be placed on $`\rho _{}`$555 Formally, it is possible to place a $`2\sigma `$ lower-limit on $`\rho _{}`$ ($`\rho _{}>0.024`$), although we consider this to be unphysical since it would result in a highly unlikely value of $`x=D_\mathrm{L}/D_\mathrm{S}>0.99`$.,

$$\rho _{}<0.0418\text{at the 3}\sigma \text{ confidence level}.$$
(10)
This constraint on $`\rho _{}`$ can be used to place a lower-limit on the mass of the lensing object. To do this, we first require an estimate of the angular size of the source star. From the colour-magnitude diagram presented in Fig. 1, it appears that the source lies on the main-sequence branch. The extinction and reddening for typical bulge stars in this field can be estimated from the position of the red clump region on the observed colour-magnitude diagram (see, for example, Albrow et al. 2000). The location of the centre of the red clump region for this field is given by $`(VI)_{\mathrm{cl},\mathrm{obs}}3.82`$ mag. Since the intrinsic dereddened colour of the red-clump region is $`(VI)_{\mathrm{cl},0}1.00`$ mag (Popowski 2000), this implies that the clump is reddened by $`E(VI)_{\mathrm{cl}}2.82`$ mag. The slope of the reddening line for this field ($`A_I/E[VI]=0.95`$ for $`l=0.1^{}`$, $`b=1.8^{}`$, Udalski 2002), yields the extinction for the red clump region, $`A_{I,\mathrm{cl}}2.68`$ mag. Therefore a star located in the red clump region, i.e. in the Galactic bulge, will undergo extinction and reddening of,

$$A_{I,\mathrm{cl}}2.68\mathrm{mag},E(VI)_{\mathrm{cl}}2.82\mathrm{mag}.$$
(11)
In Section 2 we stated that the observed magnitude and colour of sc5\_2859 is $`I_{\mathrm{obs}}18.5`$ mag and $`(VI)_{\mathrm{obs}}1.97`$ mag. Since the best-fit parallax model presented in Section 3 predicts that there is no blending ($`f_\mathrm{S}=1.000_{0.065}^{+0}`$, i.e. all of the observed flux comes from the source), this implies that the observed magnitude and colour of the source is $`I_{\mathrm{S},\mathrm{obs}}18.5`$ mag and $`(VI)_{\mathrm{S},\mathrm{obs}}1.97`$ mag. Therefore, if the source were located in the bulge (with $`D_\mathrm{S}8`$ kpc) and underwent the same amount of reddening and extinction as given in equation (11), the absolute magnitude and intrinsic colour would be $`M_I1.3`$ mag and $`(VI)_{\mathrm{S},0}0.85`$ mag. However, this is clearly incompatible with spectral types known to be in the bulge. Therefore, we conclude that the source is unlikely to be a Galactic bulge star, i.e., sc5\_2859 is more-likely an example of disk-disk lensing. By taking a simple model for the extinction we can obtain the absolute magnitude and colour of the source as a function of $`D_\mathrm{S}`$. We model the extinction using an exponential dust sheet of scale height 130pc (e.g. Drimmel & Spergel 2001), with the Sun located 20pc above the Galactic plane (Humphreys & Larsen 1995). This simple analysis suggests that the source is consistent with a K- or G-type dwarf at a distance of approximately $`1.52.5`$ kpc. For example, for $`D_\mathrm{S}2.0`$ kpc, this gives $`M_{I,\mathrm{S}}6.0`$, $`I_{\mathrm{S},0}17.5`$ and $`(VI)_{\mathrm{S},0}0.86`$, with $`A_I1.1`$ mag.

We can check this conclusion by utilising a different approach. If we assume that the source is a typical G5 main-sequence star with absolute magnitude $`M_I=4.41`$ mag and $`(VI)_\mathrm{C}=0.69`$ mag (from Cox 2000, converted into the Cousins system using Bessell 1979), then, given the slope of the reddening line for this field $`A_I/E[VI]=0.95`$ (Udalski 2002), we can calculate the predicted source distance. This method suggests that $`D_\mathrm{S}3.8`$ kpc and $`A_I1.2`$ mag. Applying the same approach to a fainter K5 main-sequence star with $`M_I=6.09`$ mag and $`(VI)_\mathrm{C}=1.26`$ mag (from Cox 2000) suggests that $`D_\mathrm{S}2.3`$ kpc and $`A_I0.7`$ mag.

Clearly, the exact brightness and colour of the source will vary depending on the assumptions; for the following analysis we shall proceed with the values calculated above using the simple extinction model and taking $`D_\mathrm{S}2.0`$ kpc. Spectroscopic observations of the source would be very useful to determine its spectral-type and hence constrain the distance to the source.

Our prediction for the intrinsic brightness and colour of the source can be used to estimate its angular size through an available empirical relationship. For example, van Belle (1999) provides the following relationship for B- to G-type main-sequence stars with $`0.4<(VK)_0<+1.5`$ mag,

$$\mathrm{log}(2\theta _{})+\frac{\mathrm{V}_{\mathrm{S},0}}{5}=0.500\pm 0.023+(0.264\pm 0.012)\times [(\mathrm{V}\mathrm{K})_{\mathrm{S},0}].$$
(12)
The source star’s intrinsic colour of $`(VI)_{\mathrm{S},0}0.86`$ can be converted into $`(VK)_{\mathrm{S},0}1.91`$ (by using, for example, Table II of Bessell & Brett 1988), which implies the star’s angular radius is,

$$\theta _{}1.08\mu \mathrm{as},$$
(13)
If the distance to the source is $`2.0`$ kpc, this corresponds to a physical source radius of $`0.5R_{}`$, i.e. slightly less than typical K- or G-type dwarfs (e.g. Table 15.8 of Cox  gives the physical radius of K- and G-type dwarfs to be between $`0.72R_{}`$ and $`1.1R_{}`$).

This value for $`\theta _{}`$ can be used to estimate the $`3\sigma `$ lower limit on the lens mass through the formula,

$$M=c^2(4G)^1\stackrel{~}{r}_\mathrm{E}\theta _\mathrm{E},$$
(14)
since the lens’ angular Einstein radius is given by $`\theta _\mathrm{E}=\theta _{}/\rho _{}`$. Using the $`\rho _{}`$ constraint provided in equation (10) gives,

$$M>0.034M_{}\text{at the 3}\sigma \text{ confidence level},$$
(15)
which corresponds to $`x=D_\mathrm{L}/D_\mathrm{S}<0.995`$.

Therefore we are only able to place a weak lower-limit on the lens mass from finite-source considerations. This outcome can be understood when one considers the likely values of the parameter $`\rho _{}`$. From equation (14), we can obtain the following relation for $`\rho _{}`$,

$$\rho _{}=0.003\left(\frac{\stackrel{~}{r}_\mathrm{E}}{10.84\mathrm{au}}\right)\left(\frac{R_{}}{R_{}}\right)\left(\frac{D_\mathrm{S}}{2.0\mathrm{kpc}}\right)^1\left(\frac{M}{M_{}}\right)^1,$$
(16)
where $`R_{}`$ is the physical source radius, which we expect to be $`1R_{}`$ if the source for sc5\_2859 is a main-sequence star. Gould (1994) showed that the detection of finite-source signatures is only possible provided $`\rho _{}>1/A_{\mathrm{max}}`$, where $`A_{\mathrm{max}}`$ is the peak magnification. For sc5\_2859 we have $`1/A_{\mathrm{max}}0.023`$, which implies that finite-source signatures are extremely unlikely to be detected.

### 4.2 Additional considerations regarding the lens mass

In addition to the above weak constraint placed on the lens mass from the finite-source analysis, there is other information that we can use from the properties of the best-fit parallax model, i.e. independent of the above finite-source fit. First, we shall consider the parallax velocity parameter ($`\stackrel{~}{v}`$), and then incorporate our knowledge of the limits on the blended flux.

For parallax microlensing events the velocity parameter $`\stackrel{~}{v}`$ (i.e. the transverse velocity of the lens relative to the source, projected onto the observer plane), must be sufficiently small so that the Earth’s orbital motion is able to affect the light curve. For sc5\_2859, this transverse velocity is $`\stackrel{~}{v}=34.2_{1.9}^{+1.5}\mathrm{kms}^1`$ directed almost parallel to the Galactic plane in the direction of rotation ($`\theta _\mathrm{G}=6.1`$ degrees; see Section 3). The velocity vector $`\stackrel{~}{𝐯}`$ is related to the transverse velocities of the observer (i.e. the Sun, $`𝐯_{}`$), source ($`𝐯_\mathrm{S}`$) and lens ($`𝐯_\mathrm{L}`$), where $`𝐯_{}`$, $`𝐯_\mathrm{S}`$ and $`𝐯_\mathrm{L}`$ are the 2-dimensional velocities perpendicular to the line-of-sight (i.e. in the lens plane), through the equation,

$$\stackrel{~}{𝐯}=\frac{𝐯_\mathrm{L}x𝐯_\mathrm{S}}{1x}𝐯_{}.$$
(17)
Previous studies of parallax microlensing events (e.g. Alcock et al. 1995) have used the parameter $`\stackrel{~}{𝐯}`$ and equation (17) to obtain a likelihood function for $`x`$ (and hence, from equation , a likelihood function for the lens mass),

$$L(x;\stackrel{~}{𝐯})\sqrt{x(1x)}\rho _\mathrm{L}(x)|\stackrel{~}{𝐯}|(1x)^3f_{\mathrm{v},\mathrm{S}}(𝐯_\mathrm{S})f_{\mathrm{v},\mathrm{L}}([1x][𝐯_{}+\stackrel{~}{𝐯}]+x𝐯_\mathrm{S})𝑑𝐯_\mathrm{S},$$
(18)
where $`f_{\mathrm{v},\mathrm{L}}`$ and $`f_{\mathrm{v},\mathrm{S}}`$ are the lens and source velocity distributions, $`\rho _\mathrm{L}(x)`$ is the density of lenses at a distance $`x`$, and all vectors are described in the plane perpendicular to the line-of-sight to the Galactic centre (i.e. two dimensional). A more thorough analysis may be obtained by including a mass function prior (see, for example, Agol et al. 2002). However, we do not attempt such an analysis here; for the purposes of this work we follow the approach of Bennett et al. (2002b) and assert that the likelihood function represents all of our knowledge about the lens mass and location (i.e. we select a uniform prior), meaning that the above likelihood function can be interpreted as the probability distribution for $`x`$.

We evaluate this function by assuming that both the lens and source reside in the disk, using density and mass distributions from Belokurov & Evans (2002) and taking the Sun’s two dimensional peculiar velocity to be 8.89 km/s in a direction $`\theta _\mathrm{G}=53.8^{}`$ (Dehnen & Binney 1998). The results of this calculation are presented in Figure 3. From this figure, we can see that the distribution for $`x`$ is not very narrow, implying that the lens can take a wide range of masses. For example, for $`D_\mathrm{S}2.0`$ kpc, this analysis suggests that $`x=0.32_{0.17}^{+0.20}`$ and hence the lens mass $`M=15.6_{8.8}^{+27.4}M_{}`$; moreover, the $`2\sigma `$ lower limit on the lens mass is $`3.1M_{}`$, suggesting that a low mass lens is strongly disfavoured.

Additional constraints can be placed on the lens nature by considering the blended flux parameter, $`f_\mathrm{S}`$. The best-fit parallax model given in Section 3 predicts that there is no blending (i.e., $`f_\mathrm{S}=1`$), meaning that all of the flux is being emitted by the lensed source. This is important because it suggests that the lens may not be a main sequence star, since if this was true one would expect some blending due to the light from the lens, i.e. one would expect $`f_\mathrm{S}<1`$.

This possibility can be investigated by considering the limits on the $`f_\mathrm{S}`$ parameter. From the best-fit parallax model, we can say that the $`2\sigma `$ and $`3\sigma `$ limits on the lens brightness are $`I_{\mathrm{L},\mathrm{obs}}>20.14`$ and $`I_{\mathrm{L},\mathrm{obs}}>19.62`$, respectively. If we model the lens with a main-sequence mass-luminosity relation $`LM^4`$, we can use this to place a lower limit on $`x`$ through equation (4). However, to do this we need to know the source distance, $`D_\mathrm{S}`$, and also how the extinction varies with distance. If we again take $`D_\mathrm{S}2.0`$ kpc and assume the previous simple extinction model (see Section 4.1), we obtain $`2\sigma `$ and $`3\sigma `$ limits of $`x>0.944`$ and $`x>0.937`$, respectively. This conclusion depends only weakly on the assumptions regarding the extinction and mass-luminosity relation. These constraints on $`x`$ appear to contradict the above conclusions from the $`\stackrel{~}{v}`$ likelihood analysis, which suggested that values of $`x`$ close to 1 are strongly disfavoured due to the low value of $`\stackrel{~}{v}`$. Figure 3 shows how the $`2\sigma `$ limits on the lens brightness compare to the likelihood distribution described in equation (18) for three values of $`D_\mathrm{S}`$. These constraints appear to be incompatible with the above likelihood analysis (for example, for $`D_\mathrm{S}=2.0`$ kpc, the likelihood analysis gives the probability that $`x`$ satisfies the $`2\sigma `$ constraints on the lens brightness to be $`<10^4`$).

Therefore, we tentatively conclude that the lens for sc5\_2859 is unlikely to be a main-sequence star, which implies that this event could be a white dwarf (although this seems unlikely, given that typical white dwarf masses of $`1M_{}`$ \[Bergeron, Ruiz, & Leggett 1997\] are disfavoured from our likelihood analysis), a neutron star, or possibly an example of lensing by a stellar mass black hole (Agol et al. 2002; Bennett et al. 2002a).

## 5 Discussion

### 5.1 The nature of event sc5\_2859

One approach that can be utilised to probe the nature of the lens for sc5\_2859 is through proper motion analysis. If the proper motion of the source can be measured, then this can be combined with the parallax model’s prediction for the transverse velocity of the lens relative to the source (i.e. $`\stackrel{~}{v}`$) to determine the proper motion of the lens. The importance of this method is that it can be applied even in the case where the lens is not luminous. By analysing one OGLE-II field, Sumi, Eyer & Woźniak (2003) have shown that proper motion measurements can be determined for objects in the OGLE-II catalogue. Preliminary results from the analysis of all OGLE-II fields suggests that the source star for sc5\_2859 may have a large proper motion (Sumi, private communication), although this requires verification as the source star is relatively faint.

Despite the fact that the light curve for sc5\_2859 appears to be reasonably well fit by the parallax microlensing model (see Section 3), it is not inconceivable that the variation in brightness could be due to some phenomena other than microlensing. Whilst this ambiguity can affect any microlensing event, for sc5\_2859 it is particularly severe since data only exist for the declining branch of the light curve. Indeed, recently Cieslinski et al. (2003) carried out a search of OGLE-II light curves in an attempt to identify new cataclysmic variables. In this work they suggested that the shape of the light curve for sc5\_2859 indicates that this event could be a candidate nova. However, it is known that microlensing events, unlike many intrinsically variable stars, are expected to be achromatic (this may not be true for heavily blended events, but for sc5\_2859 the best-fit model predicts no blended flux). For sc5\_2859 the OGLE data only provide colour information for two epochs, namely $`\mathrm{JD}2450000608`$ days, and $`\mathrm{JD}24500001020`$ days. Reassuringly, the $`(VI)`$ colour appears to stay roughly constant for these two epochs, with $`(VI)=1.97\pm 0.02`$ and $`(VI)=2.06\pm 0.07`$, respectively, which supports our belief that sc5\_2859 is a microlensing event.

### 5.2 Scientific returns from long duration events such as sc5\_2859

An important aspect of sc5\_2859 that is worth noting is that the above analysis is possible even though data exist for less than half of the microlensing light curve. Despite the fact that observations commenced after the peak, it is still possible to identify the parallactic signatures and gain well-constrained measurements of $`\stackrel{~}{r}_\mathrm{E}`$ and $`\theta `$. Obviously, this is only possible for sc5\_2859 due to the high-quality data and the long duration of the event. In future, even with initially sparse sampling, long-duration events can provide useful information provided high-quality data are obtained for the latter part of the light curve. This highlights the importance of real-time microlensing alert systems (such as the OGLE-III Early Warning System666http://www.astrouw.edu.pl/~ogle/ogle3/ews/ews.html or the MOA Transient Alert Page777http://www.roe.ac.uk/~iab/alert/alert.html) and also the importance of high-quality follow-up observations (such as those performed by the PLANET collaboration888http://mplanet.anu.edu.au).

On the other hand, if real-time alert systems can detect events at a suitably early stage, then it could be possible to make direct determinations of the lens mass. For example, once the ESO Very Large Telescope Interferometer becomes fully operational, its high sensitivity may enable measurements to be made of the angular separation of the two microlensed images (Delplancke, Górski & Richichi 2001). Coupled with a measurement of the parallax effect, this would completely break the lens degeneracy and unambiguously determine the mass of the lensing object. In future this approach could prove to be very important for long-duration microlensing events.

Another approach that would benefit from timely follow-up observations is the method of combining the parallax model with finite-source effects. As was shown in Section 4.1, if firm measurements can be made of both parallax and finite-source effects, is it possible to completely break the lens degeneracy and determine the lens mass. Although this was not possible for sc5\_2859, it has recently been shown to be feasible (for example, this approach was applied to the OGLE-II event sc26\_2218, resulting in a tentative lens mass determination \[Smith, Mao & Woźniak 2003\]). If real-time alert systems are able to identify high-magnification long-duration events suitably early, it may be possible to obtain good quality data for the peak of the light curve and hence determine the lens mass.

### 5.3 An excess of long-duration events in OGLE-II?

With the addition of sc5\_2859, this brings the total number of known events in the OGLE-II catalogue with $`t_\mathrm{E}>1`$ year to three: OGLE-1999-BUL-32999 Event OGLE-1999-BUL-32 was independently detected by the MACHO collaboration, by which it was named MACHO-99-BLG-22. After combining both MACHO and OGLE datasets, along with additional follow-up data from GMAN (Becker 2000) and MPS (Rhie et al. 1999), they find $`t_\mathrm{E}=560\pm 45`$ days (Bennett et al. 2002b). ($`t_\mathrm{E}=640_{54}^{+68}`$ days; Mao et al. 2002), sc5\_2859 ($`t_\mathrm{E}=547.6_{7.8}^{+22.6}`$ days; see Section 3 above), and OGLE-1999-BUL-19 ($`t_\mathrm{E}=372.0\pm 3.3`$ days; Smith et al. 2002). Information regarding these three events is presented in Table 2.

Table 2 shows that these three OGLE-II events are all highly magnified, with peak magnifications ranging from 10.2 to greater than 40. Figure 4 illustrates how these magnifications compare to a high-quality subset of the OGLE-II catalogue. By plotting the total baseline magnitude vs the change in magnitude at peak magnification, it can clearly be seen that these three events do not appear to lie within the main cluster. This could suggest that some, if not all, of these three long duration events are artifacts, i.e. the variability may not be due to microlensing. However, we believe that this is more-likely due to a selection effect, which implies that there could be additional long-duration events residing in the main lower-magnification cluster, i.e., events that have been omitted from existing OGLE-II microlensing catalogues due to insufficient magnification (low signal-to-noise ratio). In addition, from this table it is interesting to note that that these events do not appear to be clustered in any one direction (c.f. Popowski 2002, where it is claimed that an abundance of long duration MACHO events appear to be clustered in one particular direction).

Previous work has indicated the existence of an excess of long duration microlensing events toward the Galactic bulge (e.g., Han & Gould 1996; Zhao, Rich & Spergel 1996; Bennett et al. 2002a). However, is it possible to quantify the extent to which these three longest-duration OGLE-II events constitute an excess? If one naïvely extrapolates the behaviour of the timescale distribution of moderately long duration OGLE-II events to long timescales using a $`t_\mathrm{E}^3`$ power-law (Mao & Paczyński 1996), this predicts far fewer than three events should have $`t_\mathrm{E}>1`$ year. However, to perform this calculation thoroughly requires the detection efficiency of the OGLE-II catalogue, which is currently undetermined. We aim to investigate this important issue and report our findings in a future paper.

One implication of this potential excess is that many long-duration events could easily be overlooked by microlensing search projects. Microlensing detection algorithms usually rely on the presence of a constant baseline to differentiate true events from variable stars. However, such long-duration events may exhibit noticeable magnification over the course of many years, which can mean that a constant baseline is not observed. For example, the detection algorithm of Woźniak et al. (2001) initially failed to identify the three longest-duration events (i.e., those presented in Table 2) owing to the fact that they did not pass the constant-baseline criteria. These three events were only found when this OGLE-II catalogue was compared with a previous, independently constructed, catalogue of OGLE-II events. Also, as was mentioned above, it is conceivable that there are additional long-duration events in the OGLE-II catalogue that have been omitted due to having insufficient magnification (i.e., low signal-to-noise ratio). However, it should be noted that if some long-duration events are being overlooked in microlensing catalogues then this would increase the observed optical depth, which is already significantly larger than current theoretical predictions (see, for example, Binney, Bissantz & Gerhard 2000; Evans & Belokurov 2002; Klypin, Zhao & Somerville 2002).

## 6 Conclusion

This paper has presented a new long-duration parallax event from the OGLE-II database, sc5\_2859, which has $`t_\mathrm{E}=547.6_{7.8}^{+22.6}`$ days. We believe that both the lens and source reside in the Galactic disk, making event sc5\_2859 one of the first confirmed examples of so-called disk-disk lensing. In this aspect it differs from the two other longest-duration events, which are both believed to have source stars located in the Galactic bulge (Mao et al. 2002, Bennett et al. 2002b; Smith et al. 2002).

The source star for sc5\_2859 is likely to be located at a distance of $`D_\mathrm{S}2.0`$ kpc, although we are not able to provide a definitive value (see Section 4.1). This is the major source of uncertainty in the analysis presented in Section 4. We strongly recommend that spectroscopic observations of the source are carried out in order to determine the spectral-type of the source and hence improve the estimate of $`D_\mathrm{S}`$.

In Section 4.2 we showed that the transverse velocity of the lens relative to the source ($`\stackrel{~}{v}`$) can be used to construct a likelihood function for $`x=D_\mathrm{L}/D_\mathrm{S}`$ and hence the lens mass. This analysis suggests that lens is unlikely to be a main-sequence star, since the resulting luminosity of the lens would exceed the limits imposed on the blended flux from the parallax model. Therefore sc5\_2859 could be another candidate for microlensing by a stellar mass black hole (Agol et al. 2002; Bennett et al. 2002a).

## Acknowledgements

This work was supported by a PPARC grant. I would like to thank Shude Mao for many insightful comments and for detailed criticisms of draft versions of this manuscript. I am also grateful to Bohdan Paczyński for helpful suggestions and comments, particularly concerning the importance of event sc5\_2859 and Fig. 4, and to Andy Gould for his advice and guidance. The full 4-year data for event sc5\_2859 was generously provided by the OGLE collaboration and I gratefully acknowledge the help of Igor Soszyński and Przemysław Woźniak for their assistance. In addition, I would like to thank Takahiro Sumi for providing preliminary data from his proper motion analysis.

## Note added in proof

After this paper was accepted for publication, it was discovered that the EROS microlensing collaboration also observed the field of sc5\_2859 during this event. Preliminary analysis of these data suggest that they are in good agreement with the OGLE-II data during the period when the observations overlapped. However, there are a number of EROS data-points prior to the beginning of the OGLE-II observations (i.e. for the period $`\mathrm{JD}245000<550`$ d), and these data appear to be incompatible with the microlensing model presented in this paper. It is hoped that the definitive analysis of the EROS data, when combined with relevant follow-up observations, will help to clarify the nature of this event. The results of these findings will be presented at a later date.


format: text

text:
# 1 Initial parameters of the numerical models

Mechanisms of the Vertical Secular Heating of a Stellar Disk
N. Ya. Sotnikova (nsot@astro.spbu.ru)
S. A. Rodionov (seger@astro.spbu.ru)
Astronomical Institute, St. Petersburg State University,
Universitetskii pr. 28, Petrodvorets, St. Petersburg, 198904 Russia

Abstract

We investigate the nonlinear growth stages of bending instability in stellar disks with exponential radial density profiles.We found that the unstable modes are global (the wavelengths are larger than the disk scale lengths) and that the instability saturation level is much higher than that following from a linear criterion. The instability saturation time scales are of the order of one billion years or more. For this reason, the bending instability can play an important role in the secular heating of a stellar disk in the $`z`$ direction. In an extensive series of numerical $`N`$-body simulations with a high spatial resolution, we were able to scan in detail the space of key parameters (the initial disk thickness $`z_0`$, the Toomre parameter $`Q`$, and the ratio of dark halo mass to disk mass $`M_\mathrm{h}/M_\mathrm{d}`$). We revealed three distinct mechanisms of disk heating in the $`z`$ direction: bending instability of the entire disk, bending instability of the bar, and heating on vertical inhomogeneities in the distribution of stellar matter.

STELLAR RELAXATION IN THE DISKS OF SPIRAL GALAXIES

According to the theory of stellar evolution, the correlation between the spectral types of stars in the solar neighborhood and their kinematic parameters found more than half a century ago (Parenago 1950) is the reflection of another correlation between the stellar age and the random velocity dispersion (for a review of the currently available data, see, e.g., Fuchs et al. 2000). The latter correlation has the following pattern: the velocity dispersion for old stars is, on average, larger than that for young stars. All stars are born with a small spread in random velocities (about 5–10 km s-1. In this case, the increase in the velocity dispersion with time (dynamic heating) is the reflection of a random scattering of stars (continuous or occasional) by massive objects or large-scale density inhomogeneities.

Two pioneering studies (Spitzer and Schwarzshild 1951, 1953), in which large gaseous structures discovered twenty years later as giant molecular clouds (GMCs) from CO observations were postulated as scattering objects, were based on this idea. Subsequently, other relaxation mechanisms, for example, interaction with close satellites, were also suggested. The latter mechanism relies on observational data, according to which the disks of galaxies in interacting systems are a factor of 1.5 or 2 thicker than the disks of isolated galaxies (Reshetnikov and Combes 1997). Its efficiency is confirmed by numerical simulations (Walker et al. 1999; Velazquez and White 1999).

For isolated galaxies, apart from scattering by GMCs, the heating on inhomogeneities in the distribution of stellar matter that arise during the development of internal instabilities in the disk itself is commonly considered. The heating can be produced by transient spiral arms or a growing bar and by bending instability.

In suggesting a particular dynamic heating mechanism, the following two well-established observational facts should be borne in mind.

(1) The stellar velocity dispersion depends on time as $`\sigma _{\mathrm{tot}}t^\alpha `$, where $`\sigma _{\mathrm{tot}}`$ is the three-dimensional velocity dispersion and $`\alpha 0.330.5`$ (Fuchs et al. 2000; Binney et al. 2000).

(2) The ratio of the vertical and radial stellar velocity dispersions is $`\sigma _z/\sigma _R<1`$.

The second fact requires an explanation. It is well known that in an equilibrium axisymmetric galaxy in the absence of the third integral of motion, $`\sigma _z`$ and $`\sigma _R`$ must be equal at a given point of the disk (see, e.g., Saslaw 1989). Observational data indicate that this is not the case. The $`\sigma _z/\sigma _R`$ ratio for stars in the solar neighborhood does not depend on their spectral type (and, hence, on their age) and is $`\sigma _z/\sigma _R=0.53\pm 0.07`$ (see, e.g., Dehnen and Binney 1998). Until recently, no similar data have been available for external galaxies. At present, this ratio is known for two more nearby galaxies111In both cases, we give galaxy-averaged data.: $`\sigma _z/\sigma _R=0.70\pm 0.19`$ for NGC 488 (morphological type Sb) (Gerssen et al. 1997) and $`\sigma _z/\sigma _R=0.85\pm 0.1`$ for NGC 2985 (morphological type Sab) (Gerssen et al. 2000); i.e., it increases as one passes to earlier-type galaxies.

Spitzer and Schwarzshild (1951, 1953) showed that when stars are scattered by GMCs, the stellar velocity dispersion increases as $`\sigma _{\mathrm{tot}}t^{1/3}`$; i.e., the exponent corresponds to the lower limit of this quantity, which is still in agreement with observational data. If, however, we take into account the fact that the GMCs are confined in a thin layer and the stars lie outside this layer most of the time, then the efficiency of the relaxation mechanism associated with GMCs becomes even lower: $`\sigma _{\mathrm{tot}}t^{0.25}`$ (Lacey 1984). At the same time, GMCs are capable of effectively converting the energy of random motions in the disk plane into the energy of random motions in the perpendicular direction. Lacey (1984) found that when stars are scattered by GMCs, a ratio $`\sigma _z/\sigma _R0.8`$ is rapidly established. This ratio is much larger than the value observed for the solar neighborhood, although it is close to that obtained for the galaxy NGC 2985. However, no CO-line emission was detected for this galaxy (see, e.g., Merrifield et al. 2000).

As regards the spiral arms, numerical simulations (the first numerical results were obtained by Sellwood and Carlberg 1984) and analytical calculations (Binney and Lacey 1988; Jenkins and Binney 1990; Jenkins 1992) show good agreement between the theoretical and observational age – velocity dispersion relations. However, the following should be remembered: in the above theoretical studies, the spiral arms were assumed to effectively scatter stars only in the disk plane. Additional mechanisms are commonly invoked to convert part of the energy in the disk plane into the energy of vertical motions. The combined effect from spiral arms and GMCs yields a result that is in good agreement with observational data (Binney and Lacey 1988), but this does not rule out other vertical disk heating mechanisms. For isolated galaxies, bending instability can serve as such a mechanism. In addition, as we show below, spiral arms are capable of producing density inhomogeneities in the $`z`$ direction. The interaction of stars with these density enhancements also result in an increase of $`\sigma _\mathrm{z}`$.

Here, we study in detail the nonlinear growth stages of bending instability and show that it can be responsible for the secular increase of the stellar velocity dispersion in the vertical direction.

BENDING INSTABILITY AS A MECHANISM OF THE SECULAR HEATING OF STELLAR DISKS
The Linear Theory

Stellar disks are known to be unstable against the formation of spiral arms and bars. In addition, conditions for the growth of bending instability exist in thin equilibrium disks. This instability arises in systems with highly anisotropic particle motions and is similar to the fire-hose instability in plasma.

The bending instability of an infinitely thin gravitating layer with a nonzero velocity dispersion was first investigated by Toomre (1966). Toomre is believed to have been sceptical about his results. Therefore, they were published in a barely accessible paper and only several years ago did they become widely known after their detailed presentation by Merrit and Sellwood (1994). The bending instability of flat stellar systems has been rediscovered several times. Kulsrud et al. (1971) independently obtained a result similar to that of Toomre (1966). As regards equilibrium cold disks (with a zero velocity dispersion), Hunter and Toomre (1969) showed that they (in contrast to hot disks) are stable against the growth of bending perturbations. This fact had long been misused as an argument against any astrophysical applications related to bending instability. Only in 1977 did Polyachenko and Shukhman construct an exact linear theory for a homogeneous thin layer with a nonzero stellar velocity dispersion. A more realistic model of a thin layer with the vertical density profile described by the function $`\mathrm{sech}^2(z/z_0)`$ was analyzed by Araki (1985)222The presentation of the results of this analysis by Sellwood and Merritt (1994) made it accessible for a wide circle of researchers..

Toomre (1966) was the first to derive the dispersion relation for long-wavelength ($`\lambda =2\pi /k>>h`$, where $`2h`$ is the layer thickness) bending perturbations

$$\omega ^2=2\pi G\mathrm{\Sigma }|k|\sigma _x^2k^2,$$
(1)
where $`\mathrm{\Sigma }`$ is the surface density of the layer stars and $`\sigma _x^2`$ is the velocity dispersion along a particular coordinate in the layer plane. It follows from Eq. (1) that perturbations with wavelengths $`\lambda >\lambda _\mathrm{J}\sigma _x^2/G\mathrm{\Sigma }`$ are stable, because $`\omega ^2>0`$ in this range.

Short-wavelength perturbations

$$\lambda <\lambda _2h\frac{\sigma _x}{\sigma _z}=h\alpha ,$$

where

$$\alpha \frac{\sigma _x}{\sigma _z},$$
(2)
is the control parameter of the bending instability, must also be stable. In this case, the time of one vertical stellar oscillation, $`t_{}=h/\sigma _z`$, is longer than the time it takes for the star to traverse one wavelength, $`t_{}=\lambda /\sigma _x`$. Therefore, the star traverses a distance of several wavelengths in the time $`t_{}`$. As a result, a mismatch in the coherent particle motion arises, causing the perturbation to decay. Intermediate-wavelength ($`\lambda _2<\lambda <\lambda _\mathrm{J}`$) perturbations are unstable; the bending only grows.

When $`\lambda _2=\lambda _\mathrm{J}`$, the instability region disappears and the disk becomes stable against bending perturbations of any wavelengths. The following analytical estimate is valid for the parameter $`\alpha `$ defined by formula (2) in the linear approximation (Toomre 1966; Kulsrus et al. 1971; Polyachenko and Shukhman 1977; Araki 1985):

$$\lambda _2=\lambda _\mathrm{J}\text{for}\alpha =\alpha _{\mathrm{cr}}3.0,$$

or

$$\left(\frac{\sigma _z}{\sigma _x}\right)_{\mathrm{cr}}0.3.$$
(3)
The instability is completely suppressed for $`\sigma _z/\sigma _x>0.3`$ and grows for $`\sigma _z/\sigma _x<0.3`$.

Stars are born out of a gaseous medium and initially have low random velocities. The increase of the velocity dispersion in the radial and azimuthal directions through scattering by spiral density waves can lead to an anisotropy of the particle motions in the plane and in the vertical direction. As a result, bending instability can develop in the disk; this instability causes an increase in $`\sigma _z`$ to the level corresponding to instability saturation and an increase in the stellardisk thickness. The linear criterion gives a low instability saturation level (3). However, as was mentioned above, $`\sigma _z/\sigma _R0.5`$ in the solar neighborhood. That is why Toomre was sceptical about his discovery. However, the following two things should be remembered. First, $`(\sigma _z/\sigma _R)_{\mathrm{cr}}0.3`$ was obtained from a linear analysis. Numerical studies of the nonlinear growth stages of bending instability (Raha et al. 1991; Sellwood and Merritt 1994; Merrit and Sellwood 1994; and our numerical results) show that the instability is saturated at much larger $`\sigma _z/\sigma _R`$. Second, the growth rate of the instability is very low. The saturation time scales are several billion years.

For these two reasons, bending instability can play an important role in the secular disk heating in the $`z`$ direction.

The Numerical Model

A high spatial resolution in the $`z`$ direction is required to properly simulate the growth of bending instability in numerical simulations. In the $`N`$-body problem, this is possible only for a large number of gravitationally interacting particles. In the first numerical models, $`N`$ was taken to be 100 000 (Raha et al. 1991; Sellwood and Merritt 1994; Merrit and Sellwood 1994). We used a considerably larger number of particles, $`N=\mathrm{300\hspace{0.17em}000}\mathrm{500\hspace{0.17em}000}`$, which made it possible to reach a higher resolution and to trace the evolution of the stellar disk on time scales of about 5 Gyr. In addition, our model of a disk galaxy is much more realistic: we considered a rotating disk with an exponential density profile and assumed the existence of an additional spherical component (a dark halo). Finally, we scanned in detail the space of control parameters (see the next section) and found many of the growth features of bending instability that were overlooked both in the linear analysis and in the numerical simulations.

In our numerical simulations, we simulated the evolution of an isolated disk galaxy by using an algorithm that is based on the data structuring in the form of an hierarchical tree (Barnes and Hut 1986). Its implementation was taken from the NEMO package (http://astro.udm.edu/nemo; Teuben 1995). The package was adapted to personal computers and significantly expanded by including original data analysis and visualization programs.

In specifying a galaxy model, we separated two subsystems: the stellar disk and the spherically symmetric component (a dark halo). Star formation was disregarded. The disk was represented as a system of gravitating bodies with $`R`$ and $`z`$ density profiles, which corresponds to the observed brightness profiles for spiral galaxies:

$$\rho _\mathrm{d}(R,z)=\frac{M_\mathrm{d}}{4\pi h^2z_0}\mathrm{exp}\left(\frac{R}{h}\right)\mathrm{sech}^2\left(\frac{z}{z_0}\right),$$
(4)
where $`h`$ is the exponential disk scale length, $`z_0`$ is the typical scale of density variation in the $`z`$ direction, and $`M_\mathrm{d}`$ is the total disk mass.

The spherical component was described in terms of the external static potential

$$\mathrm{\Phi }_\mathrm{h}(r)=\frac{v_{\mathrm{}}^2}{2}\mathrm{ln}(r^2+a_\mathrm{h}^2),$$
(5)
where $`a_\mathrm{h}`$ is the typical scale and $`v_{\mathrm{}}`$ is the velocity of a particle on a circular orbit in this potential for $`r\mathrm{}`$ (($`v_{\mathrm{}}`$ can be expressed in terms of the halo mass within a sphere of a given radius and the parameter $`a_\mathrm{h}`$).

We specified the initial particle velocities (the rotation velocity and the random component) on the basis of equilibrium Jeans equations using a standard technique with a specified dependence $`\sigma _R^2\mathrm{\Sigma }(R)`$, where $`\mathrm{\Sigma }(R)`$ is the disk star surface density (see, e.g., Hernquist 1993).

All of the parameters specified in a numerical simulation can be divided into three groups: input parameters of the algorithm for solving the $`N`$-body problem, parameters of the initial model, and control parameters of the problem (the parameters that significantly affect the processes under consideration).

(1) Parameters of the algorithm: $`\delta t=0.5\times 10^6`$ yr (occasionally, $`0.25\times 10^6`$ yr) is the integration step; $`T_{\mathrm{end}}=30005000`$ Myr is the total integration time, $`\theta =0.7`$ is the parameter responsible for the accuracy of calculating the force (see, e.g., Hernquist 1987); in addition, in our computations, we took not one but the first two terms in the Laplace expansion of the potential: the monopole and quadrupole terms (Hernquist 1987); $`eps=0.02`$ kpc is the potential smoothing parameter.

(2) Parameters of the initial model: $`h=3.5`$ kpc is the exponential disk scale length; $`a_\mathrm{h}=2`$ kpc is the dark-halo scale parameter; and $`M_\mathrm{d}=(48)\times 10^{10}M_{}`$ is the disk mass.

(3) Control parameters of the problem. Let us consider them in more detail.

(A) $`z_0=0.10.5`$ kpc is the initial disk halfthickness (this quantity can be related to the initial velocity dispersion $`\sigma _\mathrm{z}`$ by assuming that the system is vertically isothermal, $`\sigma _\mathrm{z}^2=\pi Gz_0\mathrm{\Sigma }(R)`$). Since we are concerned with the increase of the stellar velocity dispersion in the $`z`$ direction, we constructed an initially thin equilibrium galaxy with a small value of $`z_0`$. This, in turn, implies a small velocity dispersion in the $`z`$ direction. Next, we observed the development of bending instability, which leads to stellar relaxation. By varying $`z_0`$, we were able to determine the instability saturation level independent of the initial conditions.

The assumption that the system is isothermal and the dependence $`\sigma _\mathrm{R}^2\mathrm{\Sigma }(R)`$ taken from empirical considerations yields the relation $`\sigma _\mathrm{z}/\sigma _\mathrm{R}=const`$ for the initial time. The ratio $`\sigma _\mathrm{z}/\sigma _\mathrm{R}`$ has always been lower than the value that follows from the linear criterion for bending instability. As a result, we were able to subsequently analyze the instability saturation level for various $`R`$.

(B) $`Q_{8.5}=1.52.2`$ is the Toomre parameter (Toomre 1964) at $`R_{\mathrm{ref}}=8.5`$ kpc (this quantity characterizes the initial radial velocity dispersion at this radius $`R_{\mathrm{ref}}`$). It follows from the dependence $`\sigma _\mathrm{R}^2\mathrm{\Sigma }(R)`$ that

$$Q=\sigma _\mathrm{R}/\sigma _\mathrm{R}^{\mathrm{cr}}\mathrm{\Sigma }^{1/2}\frac{\kappa }{\mathrm{\Sigma }}\kappa (R)\mathrm{exp}(R/2h),$$

where $`\kappa `$ is the epicyclic frequency; i.e., in contrast to $`\sigma _\mathrm{z}/\sigma _\mathrm{R}`$, $`Q`$ is not constant for the entire disk. The function $`Q(R)`$ for an exponential disk has a broad maximum in the range $`h<R<3h`$ (see Fig. 1 in Hernquist 1993). Specifying $`R_{\mathrm{ref}}2.5h8.5`$ kpc gives the condition $`Q(R)Q_{8.5}`$. This condition, in turn, ensures a level of stability against perturbations in the disk plane that is not lower than the level specified at $`R_{\mathrm{ref}}`$.

Here, we were concerned with two cases:

— a large value of $`Q_{8.5}`$ ($`2.0`$) or, in other words, a large initial radial velocity dispersion at which the formation of a bar in the disk is suppressed; this allows the development of bending instability to be traced in pure form (without the influence of a bar).

— a moderate value of $`Q_{8.5}`$ ($`1.5`$); in this case, the disk is unstable against the growth of a bar mode and we were able to investigate the influence of a bar on the relaxation in the $`z`$ direction.

(C) $`M_\mathrm{h}(4h)/M_\mathrm{d}(4h)=0.63.0`$ is the ratio of dark-halo mass to disk mass within a sphere of radius $`4h`$. The spherical component is a stabilizing factor during the development of bending instability; we consider its influence on the secular galaxy heating in the section entitled “Stellar Relaxation in Models…”.

For our computations, we used PC-compatible computers of the Astronomical Institute of the St. Petersburg State University. Some of the numerical simulations were carried out in cluster mode. The table gives the numbers of all of the models analyzed below333The models with equal control parameters have di.erent random realizations of the initial conditions — particle positions and random velocities..

Results of the Numerical Simulations

We carried out a large series of stellar-dynamics simulations and studied the growth of bending instability in thin stellar disks as a function of the control parameters of the problem. Analysis of the results of our computations revealed three distinct mechanisms of stellar relaxation in the vertical direction. The following system of units is used in all of the figures that illustrate our conclusions: the unit of time is 1 Myr; the unit of velocity is 978 km s-1; and the unit of length is 1 kpc.

Large-Scale Bending Instability of the Disk. The large-scale bending instability of the entire disk is the first galaxy heating mechanism. It is most typical of galaxies with low-mass spherical components that are initially hot in the plane, i.e., of galaxies in which the formation of a bar was suppressed (for bar-mode suppression mechanisms, see, e.g., the monograph of Binney and Tremaine 1987).

Below, we describe the scenario for the development of bending instability using Model 26\_1 as an example. All of the key evolutionary features of the stellar disk that were shown by this model were also observed in other hot models.

If we decompose $`\overline{z(R)}`$ (the mean particle deviation from the $`z=0`$ plane) into Fourier harmonics, $`\overline{z}_m(R)=A_m\mathrm{exp}(im\phi )`$, then we can calculate the amplitudes of the first three harmonics ($`m=0`$ is the axisymmetric bending or, alternatively, the bell mode; $`m=1`$ is the bending mode; and $`m=2`$ is the saddle-type mode). The change of $`A_m`$ with time describes the various bending formation stages shown in Figs. 1 and 2 (left panels). These stages are more clearly distinguished in the color two-dimensional histograms that can be found at http://www.astro.spbu.ru/staff/seger/articles/warps\_2002/fig1\_web.html. In these histograms, the mean disk particle deviation from the $`z=0`$ plane is represented by different colors (the shades of yellow and blue indicate upward and downward deviations, respectively). We see in Figs. 1 and 2 (as in the two-dimensional histograms) that a large-scale bending is produced in the galaxy at a time $`t200`$ Myr. At $`t400`$ Myr, the bending amplitude reaches its local maximum. The bending perturbation wavelength (the radial extent) is comparable to the disk scale length. The bending is not axisymmetric; the amplitude of the zero mode, i.e., the mode with the azimuthal number $`m=0`$, is small, at least less than the amplitudes of other harmonics. This relationship between the harmonic amplitudes is preserved until $`t800`$ Myr.

Next, a steadily growing axisymmetric bending mode ($`m=0`$) is clearly revealed, while all of the still large (in amplitude) nonaxisymmetric modes ($`m=1`$ and $`m=2`$) are displaced to the galactic periphery. Subsequently ($`t8001200`$ Myr), the nonaxisymmetric bending modes are gradually damped, while the axisymmetric bending reaches its maximum ($`t1000`$ Myr).

For the subsequent 1500 Myr of its evolution, the disk freezes in a shape similar to the shape of circles on water. A similar vertical disk structure, with box-shaped iosphotes in central regions, was also observed in the numerical simulations of Sellwood and Merritt (1994) for initially hot Kuzmin – Toomre disk models and in other simulations in which the bar formation was suppressed through a large Toomre parameter $`Q`$ (Patsis et al. 2002). In our computations, the amplitude of the axisymmetric bending slowly decreased and, to all appearances, the bending must subsequently disappear altogether.

Two evolutionary stages of the bending can be separated in this model: the initial stage with a nonaxisymmetric bending and the main stage with an axisymmetric bending.

Let us turn to the plot of the vertical ($`\sigma _z`$) and radial ($`\sigma _R`$) velocity dispersions against time (Fig. 3a). The value of $`\sigma _z`$ approximately doubled in time $`t=3000`$ Myr. Two times of rapid increase in the velocity dispersion $`\sigma _z`$ can be distinguished: one at $`t500`$ Myr (the time the amplitude of the initial nonaxisymmetric bending reaches its maximum) and the other at $`t10001200`$ Myr (the time the amplitude of the main axisymmetric bending reaches its maximum). As we see, the two times of change in $`\sigma _z`$ are closely related to the characteristic growth stages of bending instability. Such a relationship exists in all of the models without bars in which bending instability developed. Thus, we can conclude that in this model, the bending instability is responsible for the increase in $`\sigma _z`$; when the bell ($`m=0`$) mode appears, the energy of the random stellar velocity in the disk plane is converted most effectively into the energy of random vertical motions.

Note that the integrated quantities $`\sigma _z`$ and $`\sigma _R`$ calculated for the entire galaxy were used to construct the plots in Fig. 3. These quantities allow us to judge only the general run of the processes in the disk. Figure 4 shows the radial profiles of the azimuthally averaged $`\sigma _z`$ for various times. We see that the relaxation effect is most pronounced for central regions.

In Fig. 5, the $`\sigma _z/\sigma _R`$ ratio is plotted against $`R`$ for several times (here, we also averaged $`\sigma _z/\sigma _R`$ in concentric rings). We can conclude from our data that the saturation level of the bending instability in a region of about two exponential disk scale lengths in size (about 7 – 7.5 kpc) ismuch higher than its level predicted by the linear criterion (3). For these regions, $`\sigma _z/\sigma _R0.40.8`$.

In Figs. 1 and 2 (right panels), the edge-on image of a model galaxy is represented by isophotes for several times. The image was vertically magnified by a factor of 4. As a result, the disk bending is clearly seen. Note that the galactic evolution at late stages gives rise to a family of X-shaped orbits in central regions of the disk ($`t=1200,1400,3000`$; the left frames; the X-shaped structures are more clearly seen on the pseudoimage of the model galaxy, http://www.astro.spbu.ru/staff/seger/articles/warps\_2002/fig1\_web.html). The main axisymmetric bending may be saturated at a galaxy thickness at which a resonance arises between the frequencies of the stellar oscillations across the disk and in the rotation plane. This resonance can produce X-shaped orbits, which have repeatedly been observed in numerical simulations (Combes et al. 1990; Pfenninger and Friedli 1991; Patsis et al. 2002).

In other models, the initial and main bends manifested themselves differently that those in Model 26\_1. The behavior of all models can be assigned to one of the following cases.

(1) Both the initial and main bends arose almost simultaneously. In the long run, only the axisymmetric bending survived. Abrupt disk heating took place at the time the bending amplitude was increasing. This scenario was observed for initially thin galaxies, i.e., for the models that started from initial conditions far from the linear saturation level. Almost all of the models with $`Q_{8.5}=2.0`$ and $`z_0<0.3`$ kpc show such a pattern of evolution.

(2) The occurrences of the initial and main bends were well separated in time, as for Model 26\_1 described in detail.

(3) There was no initial bending at all, while the main bending could be observed at very late evolutionary stages. This scenario is characteristic of initially thick disks, i.e., of the models that started from conditions close to the instability saturation limit. This result closely corresponds to the following theoretical result: at a given ratio of the masses of the spherical and disk components, the mode with the azimuthal number $`m=1`$ becomes stable earlier than the other modes as the disk thickness increases (Fridman and Polyachenko 1984). As the disk thickness increases further, the m = 2 mode is ultimately saturated and the $`m=0`$ mode remains at the fore. Note that the results of our numerical simulations and the theoretical calculations of Fridman and Polyachenko (1984) are inconsistent with the conclusions of Merrit and Sellwood (1994). These authors argue that the $`m=0`$ mode is primarily stabilized in thick (radially hot) disks and the $`m=1`$ mode must remain at the fore. In our numerical simulations, as the disk thickened, the growth rate of the $`m=0`$ mode decreased sharply. A low-amplitude $`m=1`$ mode was initially observed against the background of a slowly growing $`m=0`$ mode; after the $`m=1`$ mode displaced to the periphery, an axisymmetric bending clearly showed up. The vertical velocity dispersion in thick models began to increase only after the axisymmetric bending developed. Thus, for example, the disk bending in Model 32 was observed only 3000 Myr after the beginning of the evolution.

(4) In some of the models in which a bar was formed, only a small increase in the vertical velocity dispersion was associated with the initial bending of the entire disk. This applies to all of the models with $`Q_{8.5}=1.5`$ and $`z_0<0.3`$ kpc. For thicker galaxies, no initial bending of the entire disk whatsoever was observed. In all of the models with $`Q_{8.5}=1.5`$, the bending instability of the bar itself, which is investigated in the next subsection, was the leading heating mechanism.

The Bending Instability of Bars. The second vertical disk heating mechanism revealed by our numerical simulations is related to the bending instability of bars. A bar bending was first detected by Raha et al. in three-dimensional numerical simulations. Its appearance was explained by the development of firehose instability in the bar. In our simulations, this instability was responsible for the secular disk heating almost in all of the models with bars444In some of the models that started from $`Q_{8.5}=2.0`$, despite the large velocity dispersion in the disk plane, a distinct bar was formed at late stages. However, by this time, the vertical velocity dispersion was so large (because of the relaxation associated with the bending instability of the entire disk) that no bending was formed in the bar..

Consider the disk evolution scenario using Model 8\_1 as an example. The main evolutionary stages of the bending in this initially moderately hot model are shown in the color two-dimensional histograms (Fig. 6 and 7). At early stages ($`t200400`$ Myr), the initial bending of the entire disk comes to the fore (the group of frames in the middle column). Gradually, the initial bending perturbation reaches its saturation level, is drifted to the galactic periphery, and decays ($`t600`$ Myr). A distinct bar has already been formed in the galaxy by this time (Figs. 6 and 7; the group of left frames; all images were oriented in such a way that the bar major axis was horizontal). At $`t800`$ Myr, the bending is generated already in the bar. At $`t1000`$ Myr, its amplitude reaches a maximum, after which the bending perturbation rapidly decays ($`t12001600`$ Myr).

The bar bending effect becomes understandable if we look at the group of right frames in Figs. 6 and 7. The shades of gray in these figures indicate the disk thickness in different disk regions. The thickness was calculated as $`\sqrt{\overline{z^2(R,\phi )}\overline{z(R,\phi )}^2}`$. We see that the bar was much thinner than the rest of the galaxy, which can be explained as follows. Since most of the disk stars were captured into the bar, it has a high surface density. Self-gravity causes the disk to become thinner in the bar region, which gives rise to a bending.

If we look at the plot of the vertical velocity dispersion against time (Fig. 3b), then we will see the curve rises twice: first at $`t200`$ Myr, which is clearly associated with the appearance of the initial bending of the entire disk, and, second, at $`t10001200`$ Myr, which coincides with the time the amplitude of the bar bending reaches its maximum.

In Fig. 5, $`\sigma _z/\sigma _R`$ is plotted against $`R`$ for several times. We see that for relatively cold disks (as in the case of hot disks; Fig. 5, Model 26\_1), the general saturation level of the bending instability is higher than the linear level ((3). Merrit and Sellwood (1994) pointed out that moderately hot models ($`Q1`$) behave virtually as predicted by the linear theory. Large deviations take place only for hot models. In our simulations, we observed a similar picture. The saturation level of the bending instability associated with the entire disk in the models with $`Q_{8.5}=1.5`$ is close to the linear level (see Fig. 5, $`t=800`$, Model 8\_1). However, when the bending instability of the forming bar came into effect, the final vertical disk heating was almost the same as that for the models with $`Q_{8.5}=2.0`$ (see Fig. 5, $`t=3000`$, Models 26\_1 and 8\_1).

For the remaining models with $`Q_{8.5}=1.5`$ and a low dark-halo mass ($`M_\mathrm{h}(4h)/M_\mathrm{d}(4h)<1`$), we observed the same growth stages of bending perturbations as those in Model 8\_1. If a bar was formed in the disk, then bending instability sooner or later began to develop in it. The bending amplitude rapidly increased and the instability was ultimately saturated (all of these processes took place on time scales of the order of one billion years). Note that even the bending shape was in each case similar to that observed in Model 8\_1. The models differed only by the bar formation time and by the duration of the stage that preceded the onset of bending formation in the bar itself. The larger was the dark-halo mass, the later the bar was formed and the farther was the time of bending generation in the bar from this time. The larger was the initial disk thickness, the later the bar was formed.

Stellar Relaxation in Models with a Massive Halo. Consider the stellar relaxation for models with a massive halo. A massive spherical component is known to effectively suppress the growth of bar-like instability (Ostriker and Peebles 1973) and to have a stabilizing effect on the growth of bending perturbations (Zasov et al. 1991; Sellwood 1996). The stellar relaxation in such models (if it takes place) must be produced by some other factors. Let us trace the evolution of the vertical structure of a stellar disk embedded in a massive halo using Model 9\_1 as an example.

Because of the presence of a massive spheroidal component in the galaxy, the bending instability of its disk was suppressed. No bar was formed either (at least on a time scale of 5 Gyr). Since a massive halo primarily stabilizes the disk against the growth of perturbations in the disk plane with azimuthal numbers $`m2`$, higher-order modes come to the fore. They manifest themselves in the form of multiple shortlived spiral waves and persist in the disk for several disk rotations (see, e.g., the left frame in Fig. 8). Their amplitude is initially large but the spiral pattern is blurred almost completely by the time $`t=3000`$ Myr.

Transient spirals are produced by collective processes and they are responsible for the heating of the disk in its plane. This relationship was first pointed out by Sellwood and Carlberg (1984). In Fig. 3c, the radial ($`\sigma _R`$) and vertical ($`\sigma _z`$) velocity dispersions are plotted against time. As would be expected, the rate of increase in the radial velocity dispersion slowly decreases with decreasing amplitude of the spiral arms and, after $`t=1500`$ Myr, the stellar relaxation in the disk plane related to transient spiral perturbations becomes ineffective. At the same time, for all five billion years of evolution, $`\sigma _z`$ slowly increases, although the rate of increase gradually decreases. The observed vertical secular disk heating is not the result of numerical pair relaxation, because the levels and the patterns of increase of $`\sigma _z`$ and $`\sigma _R`$ differ greatly ($`\sigma _R`$ is almost constant at late evolutionary stages). The cause of the increase of $`\sigma _z`$ in our simulations is completely different.

Let us analyze the vertical disk structure in more detail. Figure 8 shows three frames for one of the times: the galaxy seen face-on, the two-dimensional disk thickness map, and the two-dimensional vertical velocity dispersion map. We clearly see from these frames that the disk thickness in the regions where the spiral arms are located is smaller than the disk thickness in the interarm space. In other words, inhomogeneities in the distribution of stars in the plane produce inhomogeneities in the vertical distribution of stars. The effect is similar in nature to the bar thinning described in the preceding section: a self-gravitating disk is thinner where the surface density is higher.

We assume that the observed increase of the vertical velocity dispersion is related to the scattering of stars by inhomogeneities in the distribution of matter in the $`z`$ direction. Apart from the above explanation, there are two more additional facts that are indicative of this.

If we look at the two-dimensional $`\sigma _z`$ distribution (Fig. 8), then we can see that $`\sigma _z`$ in the spiral arms is higher than that in the interarm space. In addition, it follows from our computations that the decrease in the rate of increase of $`\sigma _z`$ well correlates with the decrease in the intensity of transient spirals.

Thus, we believe that the slow increase of the vertical velocity dispersion in our described models is attributable to the scattering of stars by transient spiral perturbations. The latter produce inhomogeneities in the vertical distribution of matter due to self-gravity.

Vertical Secular Disk Heating Mechanisms. During our numerical stellar-dynamics simulations, we scanned the space of control parameters (the initial disk half-thickness, the degree of disk heating in the disk plane, and the relative mass of the spheroidal subsystem). A detailed analysis of our numerical results revealed three distinct mechanisms of secular disk heating in the z direction.

(1) The large-scale bending instability of the entire disk.

(2) The bending instability of bars.

(3) The heating due to inhomogeneities in the vertical distribution of stars produced by matter inhomogeneities in the plane.

The action of a particular heating mechanism depends on the control parameters. If the bar mode is suppressed in the galaxy, then the heating is attributable to the large-scale bending instability associated with the entire disk. This is the case of a hot disk ($`Q_{8.5}>1.5`$) and (or) a moderately massive halo. In this case, the heating of the central regions is particularly strong. The saturation level of the bending instability in these regions is almost a factor of 2.5 higher than the level that follows from the linear theory. The bending mode with the azimuthal number $`m=0`$ plays a key role in the secular heating. The time of its occurrence depends on the initial disk thickness: the thicker is the disk, the farther in time is the secondary rise in disk temperature related to the growth of the bell mode.

If a bar mode develops in the disk, then the bending instability of the bar gives the largest contribution to the heating. Until a bar is formed, the maximum value of $`\sigma _z/\sigma _R`$ is 0.37 in the central regions and 0.3 on the periphery. This is the level that the bending instability associated with the entire disk reaches. It is in good agreement with the linear criterion. After the formation of a bar, at the time of its bending in the central regions, $`\sigma _z/\sigma _R`$ rises to 0.7 – 0.8. The heating on vertical bar bending perturbations dominates in initially moderately hot disks ($`Q_{8.5}=1.5`$) with a low-mass halo.

The heating on vertical inhomogeneities take place when the bending modes and the bar mode were suppressed. This was observed for moderately hot models ($`Q_{8.5}=1.5`$) with a massive halo. Note that the secular disk heating in all directions in similar simulations has been observed more than once. However, we are probably the first to associate the increase of $`\sigma _z`$ with the scattering of stars by vertical inhomogeneities.

CONCLUSIONS

We have numerically analyzed the nonlinear growth stages of bending instability in stellar disks with exponential radial density profiles and found significant deviations from the linear theory.

(1) All of the observed modes are global; i.e., the scale of unstable perturbations is larger than the typical scale of density variations in the disk. Our conclusion agrees with the conclusions of Sellwood (1996). It implies that, although the dispersion relation (1) derived for a homogeneous layer is also locally valid for inhomogeneous disks (in particular, exponential disks, as in our numerical simulations), it would be inappropriate to use this relation to analyze the saturation level of long-wavelength perturbations. This also suggests that bending instability will develop in inhomogeneous disks differently in different parts of the galaxy.

(2) The value of $`(\sigma _z/\sigma _R)_{\mathrm{cr}}0.3`$ was obtained from a linear analysis. As our numerical computations show, the saturation level of large-scale bending perturbations is a factor of 2 or 3 higher than the linear level. The differences are best seen in central regions of the stellar disk. The lower the dark-halo mass, the higher the instability saturation level and the larger the ratio of the vertical and radial velocity dispersions averaged within two exponential disk scale lengths. A similar dependence was pointed out by Mikhailova et al. (2001), although without discussing the underlying mechanisms. Early-type spiral galaxies have, on average, a smaller dark-to-luminous mass ratio and, as follows from the observational data given above, a larger value of $`\sigma _z/\sigma _R`$, in agreement with our computations.

(3) The instability saturation time scales are several billion years.

(4) Our numerical simulations revealed three distinct mechanisms of the secular heating of stellar disks in the $`z`$ direction. We confirmed the existence of bar bending instability that was first detected by Raha et al. (1991). For a large series of models, we showed that the bar bending is an inevitable stage of its evolution.

Thus, we can conclude that bending instability can play an important role in the secular disk heating in the $`z`$ direction.

ACKNOWLEDGMENTS

This work was supported in part by the Program “Leading Scientific Schools” (project no. 00-15-96607), the Federal Program “Astronomy” (project no. 40.022.1.1.1101), and the UR grant no. 02.01.006.

REFERENCES

Translated by V. Astakhov


text:
# Small World Graphs by the iterated “My Friends are Your Friends” Principle

## 1 Introduction

In 1998 Watts and Strogatz introduced in a seminal paper Small World Networks and provided a simple model of such graphs . The key property of Small World Graphs is the simultaneous presence of small diameter (at most logarithmic in the vertex size) and high clustering coefficient. Extensive investigations over the last decade of real networks like the WorldWideWeb, citation networks, friendship networks and many others have shown that these graphs are really of Small World type. But there is a second remarkable property shared by most real networks namely that they have a power law like fat tail distribution (see  for an excellent recent survey about the whole subject).

There are many models which produce graphs with small diameter and power law distribution. If additionally local search rules are incorporated one easily obtains high clustering as well. But up to now it was common believe that for the small diameter property a certain amount of essentially independent random edge formations like in classical random graphs where necessary. To obtain a scale free distribution for the degree the basic conviction was that the random edge formation should be biased according to the “the richer you are the richer you get- principle”  where richness is just measured by the degree of a vertex. In the present paper we will present a model which shows that even by entirely local edge formation rules without preferential choices graphs are obtainable which admit all the three above mentioned properties. The results are entirely based on numerical simulations since an analytic treatment of the model seems to be rather difficult. Nevertheless we hope to present in a forthcoming paper a more theoretical study of the main properties.

## 2 The iterated “My Friends are Your Friends” Principle

Looking at the edge formation process in social networks there are essential two ways to get a new contact. First by the random event of meeting somebody for instance in a train or airplane. Humans refer to this kind of relation creation often as fate. Second there are contacts which where created by the local rule “let me introduce you to one of my friends”, certainly a very common process in real life. In the language of graphs this translates to the formation of an edge between two vertices say $`x`$ and $`y`$ which had distance 2 with respect to the underlying graph metric and we refer to this

kind of edge formation as the “My Friends are your Friends” Principle (FFP). Graphs which where partially build by using the FF-Principle are already described and analyzed in . In the following we want to describe a model based on the iteration of this principle - the IFF-model. A typical example of an essential FF-graph (although a directed one and with increasing number of vertices) is the citation network.

## 3 The model description and results

In detail, we want to discuss the following model of a time dependent random graph space $`𝒢_t(N,t_0)`$ with $`N`$ vertices and an explicit integer parameter $`t_0`$ tuning the total number of edges. We start with a circular chain of $`N`$ vertices with the usual nearest neighbor edges, i.e. vertex $`1`$ is connected to vertex $`N`$ and $`2`$, vertex $`2`$ is connected to vertex $`1`$ and $`3`$ and so on. We call this configuration a circle and refer to the corresponding edges as fixed base edges. The random graph space $`𝒢_t(N,t_0)`$ will be defined through the following algorithm: In each of the first $`t_0`$ time-steps, a vertex say $`x`$ is chosen randomly. From the set of neighbors $`N_1(x,t)`$ of this vertex (at the given time $`t`$) one element say $`y`$ is chosen randomly and then randomly from $`N_1(y,t)\left\{x\right\}`$ – the neighbor set of $`y`$ at time $`t`$ with the exclusion of $`x`$ – another element say $`z`$. If there was no edge between $`z`$ and $`x`$ a new one is created between the two vertices otherwise nothing happens. In the first case we say that $`x`$ has created an FF-edge to $`z`$ respectively $`z`$ was chosen by $`x`$. After repeating $`t_0`$ times the above procedure one continues for $`t>t_0`$ the same way with the additional rule that whenever a vertex say $`x`$ has created a new FF-edge randomly an FF-edge containing $`x`$ is deleted. Note that the total number of edges stays therefore constant for $`t>t_0.`$

### 3.1 The diameter

The first surprising observation is the collapse of the diameter from $`constN`$ to $`const\mathrm{log}N`$ for sufficiently large $`t`$, $`t_0`$ and typical elements of $`𝒢_t(N,t_0)`$. This happens despite the fact that only local edge formations where used. Fig. 1 shows the diameter versus $`t`$ for different $`N`$ and $`t_0`$. Our simulations indicate that in the case $`t_0/N2`$ one gets an asymptotic regime in which the diameter is very small and shows no strong dependence on $`N`$. Tab. 1 shows for three choices of $`t_0/N`$ the mean diameter $`\overline{d}`$ in the asymptotic regime. The diameter $`\overline{d}`$ is approximately proportional to $`\mathrm{log}N`$ which is one of the characteristics of small world graphs. In the case $`t_0/N=2`$ an asymptotic regime with still small diameter but large fluctuations starts approximately at $`t/N10^4`$. Note that this value is much bigger than the collapse time in the cases $`t/N=4,10`$. It seems that $`t_0/N`$ $`2`$ is just the borderline where a collapse of the diameter appears. This corresponds to the case where the mean number of non-base edges equals $`N`$. For smaller values of $`t_0`$ we found the diameter to stay of order $`f\left(t_0\right)N`$ with an almost linear function $`f`$ (see Fig. 2). This indicates an interesting phase transition in the total number of edges reminiscent to the famous phase transition in the size of the largest component in classical random graph spaces like the $`G(N,M)`$ – the space of all graphs with $`N`$ vertices and $`M`$ edges equipped with the uniform probability measure.

Fig. 3 shows one realization of the time evolution of a typical graph from $`𝒢_t(N,t_0)`$ with $`t_0/N=4`$ and $`N=10^3`$. It illustrates how the diameter of the graph becomes very small only by using the IFF principle. Before the asymptotic state is reached an interesting symmetry breaking can be observed: the appearance of a few components with long range edges (taken the circle as a reference system) and hence small diameter. The components itself are only connected via the circle skeleton. It takes a relative long time till the components finally merge and the circle is filled uniformly with edges. The same phenomenon can be seen for larger values of $`N`$ and different $`t_0/N`$.

It is interesting to compare our model with a version of the random graph space $`G(N,M)`$ where an underlying circle is added to keep the graph connected. To be precise let $`G_c(N,M)`$ be the random graph space where $`M`$ edges are randomly added to a circle graph with $`N`$ vertices. Fig. 4 shows the diameters $`\overline{d}_{random}`$ of the average of a few samples of $`G_c(N,M)`$ as a function of $`M/N`$. Here the diameter clearly decays like $`\frac{constN}{M}`$ and reaches the $`\mathrm{log}N`$ regime already for very small values of $`M/N`$ ($`0.05`$). Note that this value is much smaller than the threshold value $`M/N`$ $`0.5`$ for the emergence of a giant component in $`G(N,M)`$. Let us remark that the $`G_c(N,M)`$ model is very close in spirit to the original Small World model of Watts and Strogatz. The main difference is with respect to the clustering coefficient since the underlying skeleton – the circle – has zero clustering coefficient. Replacing the standard circle with next neighbor connections by one where also next-next neighbor connections are edges would give an essentially equivalent model to the Watts-Strogatz one.

### 3.2 Fat tail of the degree distribution

Another unexpected property of our model is the fat tail of the degree distribution. Fig. 5a/b show the degree distribution for $`N=10^5`$ and $`t_0/N=4\left(10\right)`$ at various times.

Let us compare the degree distribution in the asymptotic regime with those of the random graph $`G_c(N,M\left(t_0\right))`$, where $`M\left(t_0\right)`$ is the expectation of the number of edges in the IFF random graph space $`𝒢_t(N,t_0)`$ (see Fig. 2 for a plot of $`M\left(t_0\right)/N`$). The resulting degree distributions are shown in Fig. 5a and Fig. 5b by big dotted lines. At the beginning (that is for $`tt_0`$) the IFF-model has essential the degree distribution of $`G_c(N,M\left(t_0\right))`$ but in the asymptotic regime the distributions differ drastically since the IFF-model has gotten a fat tail. At the moment the limited numerical data still don’t allow to check whether it is really a power law.

We want to close this section with a remark about the degree preferences if a vertex $`x`$ was chosen in the formation of a new edge. Although never explicitly included in the model there is a strong numerical evidence (see Fig. 6) that the probability of a vertex to be chosen is proportional to its degree. That is one of the basic assumptions in the Albert–Barabási model. A heuristic explanation for this property in our model is the following argument based on two independence assumptions. Let $`N_1`$ and $`N_2`$ be the expected values for the first and second neighborhood sizes. Assume that the conditional expectation of $`N_2\left(x\right)`$ for $`x`$ having degree $`k`$ equals $`kN_1`$ and assume further that for $`zN_2\left(xd\left(x\right)=k\right)`$ the expected value of $`N_2\left(z\right)`$ equals $`N_2`$. Clearly with these conditions one can easily compute the probability of a vertex $`x`$ with degree $`k`$ to be chosen within the process of an edge formation to be equal to $`\frac{kN_1}{N}\frac{1}{N_2}`$ .

### 3.3 Clustering properties

A basic quantity to measure the local clustering around a vertex $`x`$ is the number of triangles $`C_3\left(x\right)`$ containing $`x`$ as a vertex. Note that for classical random graph spaces the expectation of this number is is of order $`N^1`$. We are interested in the averaged number of triangles in the asymptotic regime of the IFF-model. Fig. 7 shows the mean number of triangles $`C_3`$ for all vertices with a fixed degree for $`t/N=600`$ and $`t/N=4000`$. The mean number of triangles per vertex (independent of its degree) is shown by the big symbols on the figure frame. The plot for the two time-values coincides practically. The number of triangles increases nearly linear with the degree $`k`$. The reason for the fluctuations for high $`k`$ is the low number of vertices with such high degrees. Fig. 7 shows also the mean clustering coefficient $`C_c`$ which is directly connected to the mean number of triangles. Namely the clustering coefficient $`C_c\left(x\right)`$ of a vertex $`x`$ with degree $`k`$ is given by the normalized triangle coefficient $`\frac{2C_3\left(x\right)}{k\left(k1\right)}`$ .

There is another interesting quantity, defined via the second shortest distance between two vertices, which characterizes additional clustering properties. For pairs of vertices with a common edge let $`d_2(x,y)`$ be the distance between $`x`$ and $`y`$ after removal of the connecting edge. Define for a given graph $`G`$ the quantity $`\phi ^{(n)}\left(G\right)`$ as the fraction of pairs of vertices with distance $`1`$ whose $`d_2`$-distance is larger $`n`$ (for a given random graph space $`𝒢`$ let $`\phi ^{(n)}\left(𝒢\right)`$ the expectation of $`\phi ^{(n)}`$). In case when the diameter of a random graph space is small due to the presence of independently generated edges the $`\phi ^{(2)}`$ value is usually very large since no short second-shortest paths between two vertices connected by an independently generated edge exist. This remains also true for the Watts-Strogatz Small World graph. In our model we get by the very construction process a small value of $`\phi `$ (3) similar to the situation met in real networks (see Tab. 2).

### 3.4 Replacing the circle by a torus

Up to now, the system of base edges forms a circle. It is natural to ask if one gets qualitatively the same results with a two-dimensional system of base edges like a torus lattice. This means that at the beginning the vertices are connected in form of a torus or a lattice with periodic boundaries. The algorithm to generate (and define) $`𝒢_t(N,t_0)`$ remains the same.

For the diameter one gets again an asymptotic regime with collapse up to logarithmic size but the asymptotics starts earlier as in the “circle”-case. Even for $`t_0/N=2`$ we get an asymptotic regime starting at $`t/N200`$. Tab. 3 shows the mean diameter $`\overline{d}`$ in the case of a torus base for various values of $`t_0`$ and $`N`$.

There is another difference to the one-dimensional case which can be seen in the realization shown in Fig. 8: the graph evolution is more uniformly with a two-dimensional base edge system than with a one-dimensional (compare to Fig. 3).

We obtain qualitatively the same results as in the one-dimensional case with respect to the degree distribution and the clustering properties for $`N=300^2`$ and $`t_0/N=4`$, $`t_0/N=10`$. In the asymptotic regime one sees again a clear fat tail distribution. The distribution of the number of triangles and the clustering coefficient for $`N=300^2`$ and $`t_0/N=4`$ gives qualitatively the same results as in the “circle”-case, i.e. the distributions are nearly equal for $`t/N=600`$ and $`t/N=4000`$ and we have a linear dependence of the number of triangles on the degree $`k.`$

### 3.5 Removal of the base edge system

Finally we want to discuss the stability of the graphs in the asymptotic regime under removal of the circle edges or torus edges, i.e. at $`t/N=1000`$ we remove the base edges and examine in which way the graph breaks into connected components and how these components evolve by continuing using our algorithm. Fig. 9 shows the number of connected components with different size $`s`$ (number of vertices) and different diameter $`d`$ for $`t_0/N=4`$ directly after deleting the base edges at $`t/N=1000`$ for the “circle” and the “torus” case. There is still one big component. The time evolution of the biggest component can be seen in Fig. 10 which shows its size $`s`$ and diameter $`d`$ versus time. An interesting result is that the diameter of the biggest component is nearly constant for $`1500<t/N<4000`$, nearly independent of $`t_0/N`$ and also nearly independent of using a circle or a torus as base edge system. Let us note that without the skeleton there can be no merging of different components by the very nature of the IFF principle. Therefore, as can easily be seen, the only stable asymptotic configurations are components which are totally connected (i.e. every vertex has an edge to all other vertices in the component) since no further changes in such components can happen. But it is very likely that the time scale till this phenomenon can be seen is huge compared to the time scales we where studying.

To compare with the classical random graph situation observe that the removal of the circle in $`G_c(N,M)`$ gives just the model $`G(N,M).`$ The component distribution for this random graph space is well known. Especially for $`M>N/2`$ there is always a giant component of size $`constN.`$ But in contrast to our model, the diameter $`d`$ of the biggest component depends via the constant on $`M\left(t_0\right)`$.

## 4 Summary and outlook

We have presented a model where the iteration of the ”My Friends are Your Friends (FF)” principle has produced random graphs with small world properties (logarithmic diameter, high clustering) and a fat tail distribution for the degree. The model used a fixed regular graph of large diameter as a skeleton, reminiscent to a pre-given geographical structure with pre-given neighborhood relations for vertex pairs. But in contrast to the small world graphs by Watts and Strogatz no random, global edge formation processes are involved to obtain small diameter. The high clustering is an immediate consequence of the FF -principle. The fat tail distribution is surprising since no preferential choice mechanisms are contained in the graph generation algorithm. So far our investigations are entirely numerical and clearly a more theoretical explanation of the observed phenomena is the desired next step.

There are several natural variants of the model which we would like to mention shortly.

First, in our model the total number of edges is kept fixed after the build up phase. Instead of that one could use probabilistic rules which keep the number of edges only fixed in mean. This would match better real situations where the FF-principle is of relevance.

Second, like in the Albert & Barabási network , a growing number of vertices could be considered. Growth should happen here in form of offsprings of already existing vertices to be able to apply the FF-principle to the ”newcomers”. In this situation it could well be the case that even without the iteration of the FF-principle – every vertex when entering the network forms just once a number of edges according to the FF-rule – Small World graphs plus fat tail for the degree distribution are obtainable. This model variant is actually the typical situation for the growth and formation of the network of citations or collaboration.