text:
# The Stellar Initial Mass Function and Beyond

## 1. Introduction

It is almost 50 years since Salpeter (1955) published the first and still most famous paper on the stellar Initial Mass Function, or distribution of masses with which stars are formed, and showed that for masses between 0.4 and 10 solar masses it can be approximated by a declining power law. Since then, this ‘Salpeter law,’ or variants or extensions of it, have been used in countless efforts to model the properties and evolution of stellar systems. Recently, for example, it has been used to model the evolution of galaxies at high redshifts. Because of its importance in many areas of astronomy, the stellar IMF has been a subject of intense study for almost five decades, and much effort has been invested in trying to determine how valid and how universal the Salpeter law or other approximations to the IMF might be, and also whether the IMF is variable, depending perhaps on parameters such as metallicity or location.

The most important result of all this work can be stated simply: Salpeter was basically right, and the power-law approximation that he proposed to the IMF between 0.4 and 10 solar masses has yet to be conclusively falsified or superseded. However, we have also learned that this power law does not apply at all masses: the IMF may not continue to follow the Salpeter law at the very largest masses, and it definitely falls below this law at the smallest masses, becoming nearly flat (in number of stars per unit log mass) below about 0.5 solar masses. This flattening of the IMF at low masses implies that most of the mass in the IMF is contained in stars whose masses are within an order of magnitude of one solar mass. Variability of the IMF has frequently been suggested, and it probably does occur at some level, but it has yet to be conclusively established. Extensive reviews of the literature on the stellar IMF have been presented by Miller & Scalo (1979), Scalo (1986, 1998), and Kroupa (2001, 2002); a briefer review has been given by Larson (1999), and a review with emphasis on young clusters has been given by Meyer et al. (2000).

## 2. Summary of The Evidence

Evidence concerning the IMF of the field stars in the solar vicinity has been collected and reviewed by Scalo (1986, 1998) and Kroupa (2001, 2002). Several proposed approximations to the IMF are all consistent with the data for masses between about 0.5 and 10 solar masses, including the original Salpeter law and the lognormal approximation of Miller & Scalo (1979). However, for masses above 10 M, the lognormal approximation now seems excluded because it falls off too steeply; although some results suggest a high-mass slope that is steeper than the Salpeter law (Scalo 1986, 1998), other results remain consistent with it (Kroupa 2001, 2002), and the Salpeter law therefore continues to be widely used as an approximation to the IMF at all masses above a solar mass. Below about 0.5 solar masses, however, the field IMF becomes almost flat in logarithmic units, and it may even decline in the brown dwarf regime. If we define the slope $`x`$ of the IMF by $`dN/d\mathrm{log}mm^x`$, the Salpeter slope is $`x=1.35`$ (where the second decimal place has never been significant and the uncertainty is at least $`\pm 0.3`$), while the slope at masses below about 0.5 M is approximately $`x0`$. Kroupa (2002) suggests an approximation to the field IMF consisting of three power laws: $`x=1.3`$ (essentially the Salpeter law) for $`m>0.5`$ M, $`x=0.3`$ (a nearly flat IMF) for $`0.08<m<0.5`$ M, and $`x=0.7`$ (an IMF declining toward small masses) for $`m<0.08`$ M. Many often inconsistent definitions and notations have been used for the IMF slope, but some alternative notations that have been widely used are the parameter $`\mathrm{\Gamma }=x`$ defined by Scalo and the parameter $`\alpha =1+x`$ used by Kroupa.

Scalo (1986) noted a possible tendency for the IMFs of star clusters to be flatter than the field IMF and more similar to the Salpeter law for masses above a solar mass, although the data show a large scatter. The scatter in the IMF slopes of clusters was illustrated graphically by Scalo (1998), who concluded that either there are real variations in the IMF or the measurements have very large uncertainties. Since the slopes found at masses above 1 M scatter roughly equally above and below the Salpeter value, most researchers have preferred the conservative interpretation that the data are consistent, within the large uncertainties, with a universal Salpeter-like IMF for masses above 1 M. However the cluster IMFs, like the field IMF, clearly flatten at smaller masses, and Kroupa’s (2002) plot of IMF slope versus average stellar mass suggests that the slope varies continuously below 1 M, and that the IMF eventually even declines in the brown dwarf regime. Recent studies of young clusters have generally yielded results consistent with a universal IMF (e.g., Meyer et al. 2000; Luhman et al. 2000; Muench et al. 2002), but some studies continue to suggest variability of the IMF; for example, the Taurus region may have a deficiency of brown dwarfs (Briceño et al. 2002), while the massive young clusters in the Galactic Center region may have a relatively flat upper IMF (Figer et al. 1999; Figer 2001; Stolte et al. 2003). Some studies of field populations have also suggested that some old populations like the Galactic Thick Disk may have a relatively flat IMF (Reylé & Robin 2001; Robin & Reylé 2003). However, many past claims of variability have not withstood the test of time, and therefore the default assumption of a universal IMF is still widely adopted.

The observed IMF is of interest also because of what it may tell us about the physics of star formation. The implications of the IMF for star formation processes are best illustrated by considering the amount of mass that goes into stars per unit logarithmic mass interval. This function is flatter by one power of mass than the number of stars per unit log mass discussed above, and it is broadly peaked around one solar mass, with a slow decline toward larger masses and a much faster decline below 0.1 M. This implies that the processes of star formation turn most of the star-forming gas into stars whose masses are within an order of magnitude of one solar mass; although significant mass also goes into the more massive stars, very little goes into low-mass stars or brown dwarfs. Therefore, in terms of where most of the mass goes, there are two main features of the IMF that need to be understood: (1) there is a characteristic stellar mass of the order of one solar mass, and (2) at larger masses the IMF has a power-law tail similar to the original Salpeter law that also contains significant mass.

## 3. The Origin of the Characteristic Mass

What processes might account for a characteristic stellar mass of the order of one solar mass? Two basic types of hypotheses have been proposed to explain this characteristic mass: (1) One possibility is that stars form from small molecular gas clumps whose masses are similar to those of the stars that from them; the stars then derive their typical masses from those of the clumps, whose properties are presumably determined by cloud fragmentation processes (Larson 1985, 1996, 1999). (2) Another possibility is that forming stars may continue to accrete mass indefinitely from an extended medium, and they may ‘determine their own masses’ by producing outflows that terminate the accretion process at some stage (Shu, Adams, & Lizano 1987; Adams & Fatuzzo 1996). Neither of these possibilities can yet be definitely excluded, and both may well play some role in determining stellar masses. However, the hypothesis that typical stellar masses are determined by cloud fragmentation processes is more amenable to direct observational test, and it appears to receive support from some recent observations.

Motte, André, & Neri (1998) have made high-resolution millimeter continuum maps of the rho Ophiuchus star-forming cloud that show many small, dense, apparently pre-stellar clumps with masses extending down to well below a solar mass (see also André et al. 1999; André, Ward-Thompson, & Barsony 2000). They found that the mass spectrum of these clumps closely resembles the stellar IMF, and they suggested on this basis that the clumps are direct stellar progenitors. Luhman & Rieke (1999) and Motte & André (2001) have noted in addition that the mass spectrum of these clumps resembles not only the standard stellar IMF described above, but also the IMF of the young stars in the $`\rho `$ Oph cloud itself. It may be premature to claim that there is quantitative agreement between the clump mass spectrum and the stellar IMF, since a similar study of the $`\rho `$ Oph clumps by Johnstone et al. (2000) finds somewhat larger clump masses and questions whether the smallest clumps found by Motte et al. (1998) are gravitationally bound. Studies of the dense clumps in the Orion clouds by Motte & André (2001) and Johnstone et al. (2001) have also yielded somewhat larger typical clump masses. However, the results of all of these studies are qualitatively very similar, and they all imply the existence of a typical clump mass of the order of one solar mass. Thus they all appear to be consistent with the view that typical stellar masses derive from the masses of small pre-stellar clumps created by cloud fragmentation processes. These results may pose a problem for the hypothesis that forming stars continue to accrete mass indefinitely from an extended medium until the accretion is shut off by outflows, since the observed clump masses are not much larger than those of the stars that form from them; however, it cannot be ruled out that outflows still play a significant role in limiting the efficiency with which stars can form in these clumps.

If the small apparently pre-stellar clumps observed in molecular clouds are created by cloud fragmentation processes, and if they are gravitationally bound, then some version of the Jeans criterion, which balances gravity against thermal pressure, almost certainly plays some role in determining their typical sizes and masses. The original Jeans (1929) analysis of the stability of an infinite uniform medium predicted a minimum mass for a collapsing clump that depends on the temperature and density of the medium, but the relevance of this analysis has often been questioned because it was not mathematically self-consistent, neglecting the overall collapse of the fragmenting medium (Spitzer 1978). However, rigorous analyses of the stability of equilibrium configurations that do not undergo overall collapse, such as sheets, filaments, and disks, yield results that are dimensionally equivalent to the original Jeans criterion (Larson 1985). Another type of rigorous analysis considers the stability of an equilibrium isothermal sphere with a given temperature and boundary pressure, and derives the size and mass of a marginally stable or ‘Bonnor-Ebert’ isothermal sphere (Spitzer 1968). The result is again dimensionally equivalent to the Jeans criterion, so this type of criterion seems likely to be of quite general relevance, at least in an approximate dimensional sense.

If star-forming clumps are produced by turbulent compression in clouds that have supersonic internal turbulent motions, as suggested by Larson (1981), then the clumps may initially have boundary pressures that are comparable to the typical turbulent ram pressure in these clouds. Assuming a temperature of 10 K and a typical turbulent cloud pressure of $`3\times 10^5`$ cm<sup>-3</sup> K, the mass and radius of a marginally stable ‘Bonnor-Ebert’ sphere are predicted to be about 0.7 M and 0.03 pc, similar to the typical masses and sizes of the pre-stellar clumps actually observed in star-forming clouds (Larson 1996, 1998, 1999). Simulations of supersonic turbulence provide at least qualitative support for the idea that star-forming clumps are created by turbulent compression, perhaps typically in clumpy filaments such as are found quite generally in many types of simulations (e.g., Ostriker, Gammie, & Stone 1999; Vázquez-Semadeni et al. 2000; Klessen 2001), regardless of the possible role of magnetic fields, but much more work is needed to clarify the quantitative details of this process and the initial properties of the collapsing clumps.

In any case, many simulations of cloud collapse and fragmentation, both with and without turbulence, have demonstrated the formation of bound clumps whose masses are typically comparable to the initial Jeans mass, and in some cases the distribution of clump masses around this value resembles qualitatively the lower stellar IMF (Larson 1978; Klessen 2001; Bate, Bonnell, & Bromm 2002a,b, 2003). For example, the bound clumps that appear in Klessen’s (2001) simulations of the collapse and fragmentation of turbulent clouds have a mass spectrum that is roughly centered on the Jeans mass and that declines at larger masses but is flatter at lower masses, as observed; similar results have also been found in the higher-resolution simulations of Bate et al. (2002a,b, 2003). In the latter simulations, significant numbers of objects are formed with masses all the way down to the opacity limit for fragmentation, which is about 0.01 solar masses; these very low-mass objects or ‘proto-brown dwarfs’ are formed by the fragmentation of the thin dense spiral filaments that appear frequently in these simulations. Thus ordinary gravitational or Jeans fragmentation, when simulated with sufficiently high spatial resolution, can yield a considerable range of masses and may be able to account for the basic observed features of the lower IMF, including the observed numbers of brown dwarfs.

Evidence for the existence of a preferred size and mass scale in the cloud fragmentation process also comes from studies of the clustering of newly formed stars (Gomez et al. 1993; Larson 1995; Simon 1997), and from studies of the scaling properties of the kinematics (Goodman et al. 1998) and the spatial structure (Williams, Blitz, & McKee 2000) of the gas in star-forming molecular clouds. As discussed by Goodman et al. (1998), all of this evidence is consistent with a picture in which star-forming clouds are dominated on large scales by turbulence and chaotic dynamics, possibly self-similar, while on scales smaller than about 0.1 parsecs the motions become more orderly, and ‘coherent cores’ can form and collapse to form individual stars or binary systems. These ‘coherent cores’ are observed to have sizes comparable to the Jeans length, as might be expected if they form on scales that are small enough for turbulence to have become unimportant but large enough that the Jeans criterion is still satisfied. While many of the details remain to be clarified, such a picture appears to offer at least a framework for understanding many features of the formation of low-mass stars, including the origin of the characteristic stellar mass.

## 4. The Origin of the Power-Law Upper IMF

The most important effects of star formation on the environment are produced by the most massive stars, so it is important also to understand the origin of the power-law upper IMF, which may have a universal Salpeter-like slope. The fact that a significant amount of mass goes into this power law means that the processes that form massive stars must be fairly efficient, and the possible universality of the Salpeter slope also suggests that these processes may be universal in nature and not dependent on local conditions. Numerous efforts have been made over the years to explain the Salpeter law, and some of the types of models that might naturally produce such a power law IMF have been reviewed by Larson (1991, 1999). These include (1) clump coagulation models, which tend under a variety of assumptions to produce a power-law clump mass spectrum; (2) models involving the continuing accretional growth of forming stars, which in simple cases can build up a power-law tail on the IMF; and (3) models postulating self-similar clustering or cloud structure, which may also naturally produce a power-law mass spectrum.

The older literature on the subject contains many models of clump coagulation that can explain the relatively flat mass spectrum ($`x0.7`$) of the CO clumps in molecular clouds, but they do not directly explain the steeper Salpeter slope ($`x1.35`$) without additional assumptions. These models do not now seem very relevant to the stellar IMF, because the CO clumps whose mass spectra they explain are of relatively low density and are not closely associated with star formation, unlike the smaller and denser continuum clumps discussed above. On the other hand, models postulating the continuing accretional growth of the most massive forming stars, such as the model suggested by Zinnecker (1982), now seem more promising when implemented in the context of realistic simulations of cluster formation. Zinnecker (1982) noted that if accreting protostars moving at a typical speed $`V`$ through a medium of density $`\rho `$ gain mass at the Bondi-Hoyle rate $`dM/dtG\rho M^2/V^3`$, then an initially peaked mass function gets stretched out toward larger masses and develops a power-law tail with a slope $`x=1`$ that is not greatly different from the Salpeter slope $`x=1.35`$. A steeper slope would result if the accretion rate were more strongly dependent on stellar mass, and this might in fact be expected since the more massive stars tend to form preferentially in the central regions of clusters where the ambient density $`\rho `$ is higher (Larson 1982; Zinnecker, McCaughrean, & Wilking 1993; Hillenbrand & Hartmann 1998; Bonnell & Davies 1998). The simulations by Bonnell et al. (1997, 2001; see also Bonnell 2000) of the accretional growth of stars in a forming cluster confirm this expectation and predict the development of a power-law tail on the IMF that has a slope similar to the Salpeter slope. Thus, continuing accretional growth of the more massive stars in a forming cluster seems likely to play an important role in the development of the upper IMF.

The tendency of stars to form in a hierarchy of groupings, and the tendency of the more massive stars to form in the larger groupings (Larson 1982; Zinnecker et al. 1993; Elmegreen et al. 2000), must also play some role in the origin of the upper IMF. A power-law IMF could result if, for example, stars form in a self-similar clustering hierarchy, and if the mass of the most massive star that forms in each subgroup of the hierarchy increases as some power $`n`$ of the mass of the subgroup; then a power-law IMF results whose slope is $`x=1/n`$ (Larson 1991). An example of this type of model would be the formation of stars in a cloud with a fractal mass distribution, such as a cloud consisting of a fractal network of filaments; for example, if stars form by the accumulation of matter at the nodes of such a network, and if they acquire masses proportional to the lengths of the filaments intersecting there, then the resulting IMF is a power law whose slope $`x`$ is equal to the fractal dimension $`D`$ of the network (Larson 1992). If this dimension were similar to the fractal dimension $`D1.4`$ characterizing the clustering of the T Tauri stars in the Taurus-Auriga clouds (Larson 1995), then an IMF with a slope $`x1.4`$ would be predicted by this model. Other types of fractal models to explain the Salpeter law, based on different assumptions, have been developed by Elmegreen (1997, 1999, 2000a).

A model based only on geometrical concepts is, however, likely to be too simple since it neglects the effects of the dynamical interactions that occur in systems of forming stars. For example, if star formation involves the accumulation of matter at the nodes of a filamentary network, then stars as well as gas will accumulate at these nodes and will interact strongly there. Interactions among the forming stars can play an important role in the star formation process itself by helping to redistribute angular momentum and enabling continuing accretion to occur (Larson 1990, 2002). Interactions may play an even more important role in the formation of massive stars than they do for low-mass stars because more interactions are likely to be involved (Larson 2002). An example of a simple model whereby a sequence of accretion events in a hierarchy of interacting protostars can lead to a Salpeter-like IMF was suggested by Larson (1999). However, the ease with which it is possible to devise explanations of the Salpeter law shows that this power law does not by itself place strong constraints on the physics of massive star formation; numerous theories have been shown to be equally capable of explaining the Salpeter law.

Important progress has recently come from high-resolution numerical simulations of cloud collapse and star formation that include the effects of continuing gas accretion as well as the effects of interactions and mergers in dense forming systems of stars (Bonnell 2002; Bonnell & Bate 2002). In the detailed simulation of cluster formation by Bonnell & Bate (2002), all of these processes occur and play important roles in the development of the upper stellar IMF. A striking result is the occurrence of a runaway increase in the density of a central core region containing some of the most massive forming stars, similar to what had been suggested by Larson (1990); this is accompanied by a runaway growth in the masses of these stars that is due to continuing gas accretion and mergers, as had been suggested by Bonnell, Bate, & Zinnecker (1998). The collapsing cluster-forming cloud also develops a filamentary structure, and groups of stars form at the nodes of the resulting network of filaments, similar to what happens in simulations of galaxy formation in cosmological models; these results may thus reflect quite general features of the evolution of systems of gas and stars under the influence of gravity and dissipation.

Like the earlier simulations of cluster formation by Bonnell et al. (1997, 2001), the more detailed simulation by Bonnell & Bate (2002) yields a power-law upper IMF that is very similar to the Salpeter law. Because many processes occur in this simulation and the dynamics becomes highly complex and chaotic, the reason for this result is not immediately apparent, but almost certainly a combination of processes is involved, including all of those discussed above and perhaps others as well. It is nevertheless worth noting, as these authors point out, that the basic physics included in this simulation is very simple – only isothermal gas dynamics, Newtonian gravity, and the merging of protostars at close distances, but magnetic fields are not included, and neither are any effects of stellar evolution or radiative feedback. Because isothermal gas dynamics and gravity introduce no new mass scale larger than the Jeans mass, the dynamics of the system becomes nearly scale-free on large scales, so it is not surprising that a power-law upper IMF should result. However, it is not obvious why the slope should be very similar to the Salpeter law, and more work will again be needed to clarify the processes involved.

## 5. The Most Massive Objects

Clearly, it is particularly important to understand better how the most massive stars are formed. It is also relevant at a meeting on “Star Formation Across the Stellar Mass Spectrum” to ask how far the stellar mass spectrum extends, and whether there is an upper limit on stellar masses. If a Salpeter IMF were to extend indefinitely to very large masses, this would imply that there is no mass limit and that stars of larger and larger mass form in systems of larger and larger size. However, the Salpeter law cannot continue indefinitely, because this would imply the existence of some stars in our Galaxy with masses of 1000 solar masses or more, which are not observed (Elmegreen 2000b). The formation of massive stars remains, unfortunately, a very poorly understood subject because massive stars form in very complex environments that are difficult to study observationally and model theoretically (Evans 1999; Garay & Lizano 1999; Stahler, Palla, & Ho 2000). Thus it may be of interest to see what can be learned just from the phenomenology of massive star formation, and from observed trends such as the dependence of maximum stellar mass on system mass.

The most massive stars form near the centers of the densest and most massive clusters (Zinnecker et al. 1993; Hillenbrand & Hartmann 1998; Elmegreen et al. 2000; Clarke, Bonnell, & Hillenbrand 2000), and the mass of the most massive star present also tends to increase systematically with the mass of the cluster (Larson 1982; Testi, Palla, & Natta 1999; Testi 2003). This is illustrated by comparing the $`\rho `$ Oph cluster, which has a mass of about 100 M and whose most massive star has a mass of about 8 M (Wilking, Lada, & Young 1989), with the Orion Nebula cluster, which has a mass of about 2000 M and whose most massive star, located in the Trapezium at its center, has a mass of about 40 M (Hillenbrand & Hartmann 1998). Even the Orion Nebula cluster is only a modest example of massive star formation, however; for example, the Arches and Quintuplet clusters near the Galactic Center have masses of at least $`10^4`$ solar masses and contain stars with masses well over 100 M (Figer et al. 1999). The most luminous and massive star known may be the ‘Pistol star’ in the Quintuplet cluster; the mass of this star is estimated to be at least 150 M, and the mass of the associated cluster is about $`2\times 10^4`$ M (Figer et al. 1998, 1999; Figer & Kim 2002). A similarly luminous star, LBV 1860-20, has been found in a similarly massive cluster that is not near the Galactic Center (Eikenberry et al. 2001). The R136 cluster in the Large Magellanic Cloud, which has a mass of about $`5\times 10^4`$ M, also contains stars with masses up to about 150 M (Hunter et al. 1995; Massey & Hunter 1998).

The examples mentioned above, namely the $`\rho `$ Oph cluster, the Orion Nebula cluster, and the Quintuplet and R136 clusters, together with a number of other young systems for which data were collected by Larson (1982), define a rough power-law correlation between the mass of the most massive star and the mass of the associated cluster which is approximately $`M_{\mathrm{star},\mathrm{max}}1.2M_{\mathrm{cluster}}^{0.45}`$. From these very limited data it is not clear that any limiting stellar mass has yet been reached, since the mass of the most massive star increases more or less continually with cluster mass up to the most massive young clusters known in our Galaxy and the LMC. However, the dependence of maximum stellar mass on cluster mass found here is not consistent with a simple extension of the Salpeter law to very large masses, since this would predict that the maximum stellar mass should vary as $`M_{\mathrm{cluster}}^{0.74}`$, rather than only as $`M_{\mathrm{cluster}}^{0.45}`$. If the data used here are reliable and representative, they suggest that the upper IMF eventually becomes steeper than the Salpeter law and approaches a slope that is closer to $`x=1/0.45=2.2`$ at very large masses. Thus the IMF may not be truly self-similar for the most massive stars but may have a slope that increases with mass, possibly reflecting an increasing difficulty in forming progressively more massive stars that is due to the effects of radiation pressure and winds.

Do stars with even larger masses form in more massive clusters, such as young globular clusters? It has long been suspected that the massive centrally condensed globular cluster M15 has a central black hole with a mass of the order of 1000 M, and this possibility is supported by recent radial velocity measurements of stars only a few tenths of an arcsecond from the center (Gebhardt et al. 2000; Gerssen et al. 2001). According to Gebhardt et al. (2000), this evidence is consistent with a black hole with a mass of $`2000`$ M, although other interpretations are not excluded, and the STIS measurements of Gerssen et al. (2001) provide further support for an increase in stellar velocity dispersion toward the center. If M15 does indeed contain such a central black hole, one explanation could be that a similarly massive star formed at the center of this cluster and then collapsed to a black hole. Stars with masses larger than 250 M are predicted to collapse completely to black holes at the end of their lifetimes, provided that they do not first lose most of their mass in winds (Fryer, Woosley, & Heger 2001), and it may be relevant in this case that mass loss is expected to be relatively unimportant for metal-poor stars (Baraffe, Heger, & Woosley 2001) and that M15, with $`[\mathrm{Fe}/\mathrm{H}]=2.0`$, is one of the most metal-poor clusters known. M15 might thus have provided a particularly favorable environment for the formation of a very massive star that could have collapsed to a massive black hole.

Other centrally condensed systems are also known that contain central black holes: supermassive black holes are found at the centers of most, if not all, galaxy bulges, and the masses of these nuclear black holes increase systematically with the mass of the bulge (Kormendy & Richstone 1995; Kormendy & Gebhardt 2001). If M15 contains a black hole of mass $`2000`$ M, it falls on the same relation between black hole mass and system mass that was found by Kormendy & Richstone (1995) for galaxy bulges, since these authors found that the mass of the central black hole is typically about 0.002 times the mass of the system, which is also approximately true for M15. The updated reviews of Kormendy & Gebhardt (2001) and McLure & Dunlop (2002) find, for a much larger sample of galaxies, a very similar relation between black hole mass and bulge mass, with a slightly smaller value of 0.0013 for the typical ratio of black hole mass to bulge mass.

Could there be a connection between the formation of massive stars at the centers of clusters and the formation of supermassive black holes at the centers of galaxies? In both cases, the mass of the central object increases systematically with the mass of the system, and this suggests that similar processes may be involved in building up the central object (Larson 2002). For example, runaway accretion and coalescence processes like those discussed above could play a role not only in forming massive stars in clusters but in building supermassive black holes in galactic nuclei. The relation between central object mass and system mass is not identical in the two cases, however, since the mass of the most massive star in a cluster increases with system mass only as $`1.2M_{\mathrm{system}}^{0.45}`$, while the mass of the central black hole in a galaxy bulge increases more nearly linearly with system mass and is approximately equal to $`0.0015M_{\mathrm{system}}`$. These two relations intersect at a system mass of about $`2\times 10^5`$ M and a central object mass of about 300 M, suggesting that a transition may occur at about this point from a regime of normal massive star formation in systems of smaller mass to a regime of black hole formation and growth in systems of larger mass. It is consistent with this possibility that the transition appears to occur at a maximum object mass of around 300 M, since this is approximately the same as the mass of 250 M above which stars are predicted to collapse completely to black holes; this suggests that systems more massive than a few times $`10^5`$ M may sometimes make stars massive enough to collapse completely to black holes.

Once a massive black hole has formed, the possibilities for the further growth in mass of this object are increased because such a black hole, unlike a massive star, is a permanent feature of the system that can only grow in mass by any continuing interaction with its surroundings. For example, a central black hole might continue to accrete matter from neighboring massive stars in a central subcluster whose members might otherwise have merged. Such processes might have contributed to the mass of the possible black hole in M15, since the maximum stellar mass predicted by the above relation for a cluster with the mass of M15 is only about 600 M, and this suggests that some further growth in mass might have been required. In the case of galaxy bulges, which are built up at least partly by mergers of smaller systems, the mergers may lead to the continuing growth of the nuclear black holes through the accretion of gas, stars, or other black holes that become concentrated at the center during each merger, and such processes might account for the existence of a relation between black hole mass and bulge mass (Cattaneo, Haehnelt, & Rees 1999; Kauffmann & Haehnelt 2000; Larson 2000).

Even without gas dynamics, runaway stellar mergers can occur in massive and dense young star clusters when they experience a stellar dynamical ‘core collapse’, and this can lead to the building up of very massive stars and black holes at their centers. Portegies Zwart et al. (1999, 2002) and Figer & Kim (2002) suggest that stellar mergers might account for the existence of such very massive stars as the Pistol star, and Ebisuzaki et al. (2001) and Portegies Zwart & McMillan (2002) suggest that black holes with masses of hundreds to thousands of solar masses may be built up by similar processes at the centers of massive young clusters. X-ray observations provide evidence that such ‘intermediate mass black holes’ exist in some very luminous young clusters in starburst galaxies (Ebisuzaki et al. 2001; Portegies Zwart & McMillan 2002), and this evidence supports the possibility that massive black holes can indeed form at the centers of massive young clusters. Portegies Zwart & McMillan (2002) speculate further that the accumulation of such clusters at the centers of galaxies, accompanied by the merging of their black holes, could lead to the building up of supermassive black holes in galactic nuclei, and they argue that central objects with masses of the order of 0.1 percent of the system mass should result quite generally from such processes.

## 6. Summary

Much evidence indicates that there are two basic general features of the stellar IMF that need to be understood: there is a characteristic stellar mass of the order of one solar mass, and the IMF also has a power-law tail toward larger masses that is similar to the original Salpeter law. Variability of the IMF has often been suggested, but has not been conclusively established. The characteristic mass may derive from the typical masses of the small, dense, apparently pre-stellar clumps that have been observed in a number of star-forming clouds. In many simulations of cloud collapse and fragmentation, the typical masses of the clumps formed are similar to the Jeans mass, regardless of whether turbulence plays an important role; magnetic fields also do not appear to make a major difference. Gravitational fragmentation produces in addition significant numbers of proto-brown dwarfs by the fragmentation of thin dense spiral filaments, and it may thus be able to account for the basic observed features of the lower stellar IMF. The Salpeter power-law that approximates the IMF at higher masses is less well understood, but it probably results from a combination of processes including the effects of continuing gas accretion by the more massive forming stars and the effects of interactions and mergers in dense forming systems of stars. Simulations of the formation of star clusters show that the dynamics of these accretional growth processes becomes almost scale-free at large masses, and that these processes can produce a power-law upper IMF that is very similar to the Salpeter law.

Observations also indicate that the mass of the most massive star that forms in a cluster increases systematically with the mass of the cluster, although not as rapidly as is predicted by the Salpeter law, and there is no clear indication of an upper stellar mass limit. Some of the most massive known clusters, including the globular cluster M15 and some luminous young clusters in starburst galaxies, may contain massive black holes that originated from the collapse of very massive stars that formed near the centers of these systems. The observed dependence of maximum object mass on system mass suggests that stellar systems more massive than a few times $`10^5`$ solar masses may sometimes make stars massive enough to collapse completely to black holes.

## References

Adams, F. C., & Fatuzzo, M. 1996, ApJ, 464, 256

André, P., Motte, F., Bacmann, A., & Belloche, A. 1999, in Star Formation 1999, ed. T. Nakamoto (Nobeyama Radio Observatory, Nobeyama), p. 145

André, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 59

Baraffe, I., Heger, A., & Woosley, S. E. 2001, ApJ, 550, 890

Bate, M. R., Bonnell, I. A., & Bromm, V. 2002a, MNRAS, 332, L65

Bate, M. R., Bonnell, I. A., & Bromm, V. 2002b, in The Origins of Stars and Planets: The VLT View, ed. J. F. Alves & M. J. McCaughrean (Springer-Verlag, Berlin), p. 139

Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, in Galactic Star Formation Across the Stellar Mass Spectrum, ASP Conf. Ser. Vol. 287, ed. J. M. De Buizer & N. S. van der Bliek (ASP, San Francisco), p. 427 (this volume)

Bonnell, I. A. 2000, in Star Formation from the Small to the Large Scale, 33rd ESLAB Symposium, ed. F. Favata, A. A. Kaas, & A. Wilson (ESA SP-445, ESA, Noordwijk), p. 273

Bonnell, I. A. 2002, in Hot Star Workshop III: The Earliest Stages of Massive Star Birth, ASP Conf. Ser. Vol. 267, ed. P. A. Crowther (ASP, San Francisco), p. 163

Bonnell, I. A., & Bate, M. R. 2002, MNRAS, 336, 659

Bonnell, I. A., & Davies, M. B. 1998, MNRAS, 295, 691

Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS, 285, 201

Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS, 298, 93

Bonnell, I. A., Clarke, C. J., Bate, M. R., & Pringle, J. E. 2001, MNRAS, 324, 573

Briceño, C., Luhman, K. L., Hartmann, L., Stauffer, J. R., & Kirkpatrick, J. D. 2002, ApJ, 580, 317

Cattaneo, A., Haehnelt, M. G., & Rees, M. J. 1999, MNRAS, 308, 77

Clarke, C. J., Bonnell, I. A., & Hillenbrand, L. A. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 151

Ebisuzaki, T., Makino, J., Tsuru, T., Funato, Y., Portegies Zwart, S., Hut, P., McMillan, S., Matsushita, S., Matsumoto, H., & Kawabe, R. 2001, ApJ, 562, L19

Eikenberry, S. S., Matthews, K., Garske, M. A., Hu, D., Jackson, M. A., Patel, S. G., Barry, D. J., Colonno, M. R., Houck, J. R., & Wilson, J. C. 2001, BAAS, 33, 1448

Elmegreen, B. G. 1997, ApJ, 486, 944

Elmegreen, B. G. 1999, ApJ, 515. 323

Elmegreen, B. G. 2000a, MNRAS, 311, L5

Elmegreen, B. G. 2000b, ApJ, 539, 342

Elmegreen, B. G., Efremov, Y., Pudritz, R. E., & Zinnecker, H. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 179

Evans, N. J. 1999, ARA&A, 37, 311

Figer, D. F. 2001, in Starburst Galaxies: Near and Far, ed. L. Tacconi & D. Lutz (Springer-Verlag, Berlin), p. 13

Figer, D. F., & Kim, S. S. 2002, in Stellar Collisions, Mergers, and Their Consequences, ASP Conf. Ser. Vol. 263, ed. M. M. Shara (ASP, San Francisco), p. 287

Figer, D. F., Najarro, F., Morris, M., McLean, I. S., Geballe, T. R., Ghez, A. M., & Langer, N. 1998, ApJ, 506, 384

Figer, D. F., Kim, S. S., Morris, M., Serabyn, E., Rich, R. M., & McLean, I. S. 1999, ApJ, 525, 750

Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372

Garay, G., & Lizano, S. 1999, PASP, 111, 1049

Gebhardt, K., Pryor, C., O’Connell, R. D., Williams, T. B., & Hesser, J. E. 2000, AJ, 119, 1268

Gerssen, J., van der Marel, R. P., Dubath, P., Gebhardt, K., Guhathakurta, P., Peterson, R., & Pryor, C. 2001, BAAS, 33, 1385

Gomez, M., Hartmann, L., Kenyon, S. J., & Hewett, R. 1993, AJ, 105, 1927

Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223

Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540

Hunter, D. A., Shaya, E. J., Holtzman, J. A., Light, R. M., O’Neil, E. J., & Lynds, R. 1995, ApJ, 448, 179

Jeans, J. H. 1929, Astronomy and Cosmogony (Cambridge Univ. Press, Cambridge; reprinted by Dover, New York, 1961)

Johnstone, D., Wilson, C. D., Moriarty-Schieven, G., Joncas, G., Smith, G., Gregersen, E., & Fich, M. 2000, ApJ, 545, 327

Johnstone, D., Fich, M., Mitchell, G. F., & Moriarty-Schieven, G. 2001, ApJ, 559, 307

Kauffmann, G., & Haehnelt, M. 2000, MNRAS, 311, 576

Klessen, R. S. 2001b, ApJ, 556, 837

Kormendy, J., & Gebhardt, K. 2001, in The 20th Texas Symposium on Relativistic Astrophysics, ed. J. C. Wheeler & H. Martel (AIP Conference Proceedings Vol. 586, AIP Press, New York), p. 363

Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581

Kroupa, P. 2001, MNRAS, 322, 231

Kroupa, P. 2002, Science, 295, 82

Larson, R. B. 1978, MNRAS, 184, 69

Larson, R. B. 1981, MNRAS, 194, 809

Larson, R. B. 1982, MNRAS, 200, 159

Larson, R. B. 1985, MNRAS, 214, 379

Larson, R. B. 1990, in Physical Processes in Fragmentation and Star Formation, ed. R. Capuzzo-Dolcetta, C. Chiosi, & A. Di Fazio (Kluwer, Dordrecht), p. 389

Larson, R. B. 1991, in Fragmentation of Molecular Clouds and Star Formation, IAU Symposium No. 147, ed. E. Falgarone, F. Boulanger & G. Duvert (Kluwer, Dordrecht), p. 261

Larson, R. B. 1992, MNRAS, 256, 641

Larson, R. B. 1995, MNRAS, 272, 213

Larson, R. B. 1996, in The Interplay Between Massive Star Formation, the ISM and Galaxy Evolution, ed. D. Kunth, B. Guiderdoni, M. Heydari-Malayeri, & T. X. Thuan (Editions Frontières, Gif sur Yvette), p. 3

Larson, R. B. 1998, unpublished, in The Orion Complex Revisited, ed. M. J. McCaughrean & A. Burkert (ASP Conference Series, in press), available at www.astro.yale.edu/larson/pubs.html

Larson, R. B. 1999, in Star Formation 1999, ed. T. Nakamoto (Nobeyama Radio Observatory, Nobeyama), p. 336

Larson, R. B. 2000, in Star Formation from the Small to the Large Scale, 33rd ESLAB Symposium, ed. F. Favata, A. A. Kaas, & A. Wilson (ESA SP-445, ESA Publications, Noordwijk), p. 13

Larson, R. B. 2002, MNRAS, 332, 155

Luhman, K. L., & Rieke, G. H. 1999, ApJ, 525, 440

Luhman, K. L., Rieke, G. H., Young. E. T., Cotera, A. S., Chen, H., Rieke, M. J., Schneider, G., & Thompson, R. I. 2000, ApJ, 540, 1016

Massey, P., & Hunter, D. A. 1998, ApJ, 493, 180

McLure, R. J., & Dunlop, J. S. 2002, MNRAS, 331, 795

Meyer, M. R., Adams, F. C., Hillenbrand, L. A., Carpenter, J. M., & Larson, R. B. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 121

Miller, G. E., & Scalo, J. M. 1979, ApJS, 41, 513

Motte, F., André, P., & Neri, R. 1998, A&A, 336, 150

Motte, F., & André, P., 2001, in From Darkness to Light: Origin and Evolution of Young Stellar Clusters, ASP Conf. Ser. Vol. 243, ed. T. Montmerle & P. André (ASP, San Francisco), p. 301

Muench, A. A., Lada, E. A., Lada, C. J., & Alves, J. 2002, ApJ, 573, 366

Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 513, 259

Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899

Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 1999, A&A, 348, 117

Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 2002, in Stellar Collisions, Mergers, and Their Consequences, ASP Conf. Ser. Vol. 263, ed. M. M. Shara (ASP, San Francisco), p. 95

Reylé, C., & Robin, A. C. 2001, A&A, 373, 886

Robin, A. C., & Reylé, C. 2003, in Galactic Star Formation Across the Stellar Mass Spectrum, ASP Conf. Ser. Vol. 287, ed. J. M. De Buizer & N. S. van der Bliek (ASP, San Francisco), p. 104 (this volume)

Salpeter, E. E. 1955, ApJ, 121, 161

Scalo, J. 1986, Fundam. Cosmic Phys., 11, 1

Scalo, J. 1998, in The Stellar Initial Mass Function, ASP Conf. Ser. Vol. 142, ed. G. Gilmore & D. Howell (ASP, San Francisco), p. 201

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

Simon, M. 1997, ApJ, 482, L81

Spitzer, L. 1968, in Nebulae and Interstellar Matter, ed. B. M. Middlehurst & L. H. Aller (Univ. of Chicago Press, Chicago), p. 1

Spitzer, L. 1978, Physical Processes in the Interstellar Medium (Wiley-Interscience, New York)

Stahler, S. W., Palla, F., & Ho, P. T. P. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 327

Stolte, A., Grebel, E. K., Brandner, W., & Figer, D. F. 2003, in Galactic Star Formation Across the Stellar Mass Spectrum, ASP Conf. Ser. Vol. 287, ed. J. M. De Buizer & N. S. van der Bliek (ASP, San Francisco), p. 433 (this volume),

Testi, L., Palla. F., & Natta, A. 1999, A&A, 342, 515

Testi, L. 2003, in Galactic Star Formation Across the Stellar Mass Spectrum, ASP Conf. Ser. Vol. 287, ed. J. M. De Buizer & N. S. van der Bliek (ASP, San Francisco), p. 163 (this volume)

Vázquez-Semadeni, E., Ostriker, E. C., Passot, T., Gammie, C. F., & Stone, J. M. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 3

Wilking, B. A., Lada, C. J., & Young, E. T. 1989, ApJ, 340, 823

Williams, J. P., Blitz, L., & McKee, C. F. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona Press, Tucson), p. 97

Zinnecker, H. 1982, in Symposium on the Orion Nebula to Honor Henry Draper, ed. A. E. Glassgold, P. J. Huggins, & E. L. Schucking (Ann. New York Academy of Sciences, Vol. 395, New York), p. 226

Zinnecker, H., McCaughrean, M. J., & Wilking, B. A. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine (Univ. of Arizona Press, Tucson), p. 429


format: text

text:
# Correlated X-ray and optical variability in KV UMa

## 1 Introduction

The transient X-ray source XTE J1118+480 (=KV UMa) appeared in January 2000 (Remillard et al. 2000) and (with an interruption in February) was active as an X-ray source till end of July. The source appears to be a relatively nearby (2 kpc) binary with an orbital period of 4 hrs, and a K-star companion. Its mass function of 6 M makes the primary a solid black hole candidate, (Mc Clintock et al. 2001a, Wagner et al. 2001). The X-ray outburst was somewhat unusual compared with ‘standard’ transients such as X-ray Nova Muscae. The X-ray flux was relatively constant during the main outburst episode (March–July). The X-ray luminosity (1-100 keV) was only about $`10^{36}`$ erg/s, and the source was always found in a X-ray hard state (most of the X-ray energy flux appearing around 100keV). Since hard state spectra are the rule during the decay phase of X-ray transients, this suggests that the accretion rate of the system was significantly below the Eddington value.

Due to the low galactic absorption at the position of the source, its spectral energy distribution (SED) could be determined from the IR to X-rays, including the crucial EUV part of the spectrum (Hynes et al. 2000, McClintock et al. 2000b), a unique result for X-ray transients. These observations showed that the source had an unusually high optical-to-X-ray flux ratio compared with other transients. The SED appears to peak at 100 eV, and the inferred UV-to-EUV flux of $`10^{36}`$ erg/s was about the same as the X-ray flux. The optical-to-X-ray flux ratio was about 0.1. High optical brightness during a hard X-ray state has been seen before in the unusual long-term variable source GX 339-4 (Motch et al. 1983).

The long duration of this transient, its location above the galactic plane, and its relatively close distance made it an ideal target for simultaneous X-ray and optical observations. Though such coordinated observations have been attempted regularly, successful observations at high (sub-second) time resolution have been reported previously in only two cases. Motch et al. (1981) reported correlated fast optical and X-ray variability from 96s of observations of GX 339-4, when it was in one of its rapidly variable and optically bright states. The X-ray-optical cross-correlation of these data was reported to show a negative correlation preceding the X-rays by a few seconds. Hynes et al. (2000) reported simultaneous HST and XTE observations of the black hole transient GRO J1655-40. The cross-correlation showed, in contrast with the case of 339-4, a positive correlation peak following the X-rays, interpreted by the authors as evidence of reprocessing of X-rays by the accretion disk.

The optical emission of black hole candidates in outburst thus appears puzzling, both in its unpredictable amplitude, and its contradictory correlations with the X-rays. Theoretically, the optical light is interesting because it may provide a way to probe the structure of the accretion flow, if the light is caused by reprocessing of X-rays from the central object (e.g. Arons & King 1993). Evidence for such reprocessing of X-rays by an accretion disk has been found in type-1 X-ray bursts from accreting neutron stars (Matsuoka et al. 1984, Turner et al. 1985, van Paradijs et al. 1990, Kong et al. 2000).

On the other hand, lower limits to the optical brightness temperatures inferred from the short optical variability in GX 339-4 led Fabian et al. (1982) to conclude that a more likely mechanism for the optical radiation in the Motch et al. data is thermal cyclosynchrotron emission in a strong magnetic field. We reached a similar conclusion in a first description and interpretation of the results reported here (Kanbach et al. 2001). If this interpretation is correct, the regions contributing most to the optical-UV light would be much closer to the black hole than in a reprocessing model, at least in objects like KV UMa and GX 339-4, and optical light of black hole transients would potentially be a powerful new diagnostic of the inner accretion flow.

In Kanbach et al. (2001) we found that the average X-ray/optical cross correlation function consists of a peak starting at zero lag, preceded by a ‘precognition dip’. It thus combines both of the effects seen previously by Motch et al. and Hynes et al. This raises the question how these two components of the cross correlation are related, in particular since the dip at negative lag is not very easily explained in models proposed so far. In this paper we address this question by a more detailed statistical analysis of our observations of KV UMa. The theoretical interpretation will be developed further in a future paper.

## 2 Observations

X-ray and optical observations were collected on the nights of 4,5,6 and 7 July 2001. X-ray observations were made by the Rossi X-ray Timing Explorer, for about one orbit on each of these nights. The total duration of useful simultaneous observations was 2.5 hours. The optical observations were done with the OPTIMA photometer (Straubmeier et al. 2001) attached to the 1.3 m telescope on Mt. Skinakas, Crete (http://observatory.physics.uoc.gr/), which is operated jointly by the Foundation for Research and Technology Hellas (FORTH) and the Max-Planck-Institut für Extraterrestrische Physik. The photometer consists of a cluster of fiber-coupled avalanche diodes. The diodes are sensitive from 450 to 950 nm at a mean quantum efficiency of 50%. Individual photon arrival times are recorded to an accuracy of 2$`\mu `$s with a GPS-based clock, at a maximum count rate of 200 000/s. On several of the nights, the telescope was facing into a fairly strong wind, causing modulations of the signal at frequencies around 9 and 40 Hz due to telescope vibration. Most of the results reported here are based on the correlation between optical and X-ray signals, for which this artefact has little consequence. In the cross-correlation functions it appears as a source of noise but has no systematic effects.

The average total X-ray count rate was 500 s-1 (3 detectors), typical for the brightness of the source during most of its outburst. The optical brightness was around V=13, resulting in a count rate of about 24000 $`s^1`$. While the X-ray flux was highly variable, with amplitudes typical of the hard states of black hole transients or Cyg X-1, the rms variability of the optical flux was only about 10% on time scales of 100ms–30s.

## 3 Analysis

An example of the X-ray and optical light curves is shown in Fig. 1. The data show the original count rates, without background or baseline subtractions. Peaks in the X-rays with durations of a few seconds as seen in this example were fairly typical. Visual inspection of such events suggests that the corresponding optical signal consists of a dip, preceding the X-rays by several seconds, followed by a positive response peaking a few seconds after the X-rays.

In order to quantify this impression in a more objective way, an automatic peak-finding algorithm was devised to collect and superimpose these events. The X-ray time series was compared with a box-car smoothed version of itself (with a width of 30s). Time intervals where the actual signal was more than twice the smoothed signal were identified as peaks. The 100 events found in this way were aligned in time on the maxima of the X-ray peaks.

The average X-ray and optical light curves of the peak events found in this way is shown in Fig. 2. This average confirms the subjective impression given by Fig. 1. The dip in the optical preceding the X-rays appears to be a systematic effect, though its amplitude varies between events.

### 3.1 Cross-correlation

The superposition of events shown in Fig. 2 is limited in its time resolution by the width of the X-ray peaks, of the order of a second. A view at much higher time resolution, at the expense of a somewhat less direct interpretation, is given by the cross correlation between the data. The cross-correlation of the X-ray time series $`S_\mathrm{X}`$ and the optical series $`S_\mathrm{O}`$, $`C_{XO}(t)=S_\mathrm{X}S_\mathrm{O}`$, where a $``$ denotes the convolution, was computed for 297 independent segments of data of 25s length. The average of these cross-correlations is shown in Fig. 3.

The cross correlation contains a residual of the photon noise in the original time series. In Figs. 3, 4 this shows up as the weak high-frequency noise component. Its amplitude is much lower, however, than the noise due to the intrinsic variability of the cross correlation, as Fig. 4 shows.

The cross-correlation appears to contain significant signal for lags (optical later than X-rays) from about -5 to +5s. On these time scales, the correlation shows the same properties as the direct comparison of the signals shown in Fig. 2: a dip followed by a peak of about the same amplitude.

As the right panel of Fig. 3 shows, the change from the dip to the peak takes place quite close to zero lag. The slope of the correlation increases by a factor 10 within about 30 ms at $`t=0`$, and 50% of the peak is reached at $`t=100`$ ms.

### 3.2 Dip and peak

A peak starting at $`t=0`$ is understandable in a reprocessing model for the optical light, and its duration is also of the order expected for reprocessing of X-rays by an accretion disk. Evidence for such reprocessing has been seen in X-ray bursters (Matsuoka et al. 1984, van Paradijs et al. 1990, Kong et al. 2000). It has also been advanced to explain the optical response in some black hole transients (O’Brien & Horne 2000, Hynes et al. 1998). A delayed optical response is also a natural result in the outflow model of Kanbach et al. (2001). In this model, the unusually large optical flux of KV UMa does not come from the accretion disk but instead is produced by synchrotron emission in a magnetic outflow. The dip in the optical light preceding the X-rays is not very natural in either of these explanations, however. This raises the question whether the dip should be regarded as an independent phenomenon, or as a part of the correlation function associated directly with the peak.

The question whether the dip represents an independent physical process from that producing the peak can be answered to some extent by investigating the variability of the dip and the peak between different samples of the cross correlation. If their amplitudes turn out to be statistically uncorrelated, it is unlikely that they are due to the same process. Significant variability is in fact present in both the dip and the peak. This is shown in Fig. 4, which presents a random subset of 16 of the 300 correlations. From this figure it is also clear, however, that a large part of the variability of the cross correlation appears to be independent of both the dip and the peak. In many cases, this additional component masks the dip and/or the peak. From simple inspection of figures like these, it is hard to determine if this additional variability is just some form of intrinsic random variability in the source, or whether it contains additional systematic variations.

## 4 Principal components of the cross correlation

A more quantitative way of investigating the variability of a signal like the cross correlations in Fig. 4 is by a principal component analysis (PCA, e.g. Kendall 1980, not to be confused here with the PCA detectors on board of XTE). If the signal is variable, and a sufficiently large set of samples of it available, this method can decompose the samples into components whose amplitudes vary statistically independently through the set. The method has been used in astronomy, among other applications, for the classification of galaxy spectra from large surveys (e.g. Folkes et al. 1999), AGN spectra (Boroson 2002), Chepheid light curves (Kanbur et al. 2002), and velocity distributions in molecular clouds (Brunt & Kerton 2002). The components found in this way can often be interpreted more easily in terms of variations in physically quantities than the original (large) datasets.

In practice, significant care has to be excercised in interpretating the principal components. In the following we do this by comparison with forward modeling. First we apply PCA to our set of 297 samples of the cross correlation and interpret a few of the most significant components. We then verify the stability and significance of these components by constructing sets of synthetic cross correlation functions, and processing them in the same way. These synthetic samples consist of superpositions of a few input components whose shapes are guessed, with statistical properties similar to those of the real set of cross correlations. The input components and the statistics of their amplitudes are adjusted to find the minimal number of ‘true’ components, and their shapes, that are needed to reproduce the pattern of components found by the PCA of the real data set.

### 4.1 Method

For our application, the PCA is used as follows. Let the X-O cross correlation function of sample $`j`$ be denoted by $`𝐚^j(t)`$. Let $`m`$ be the number of available samples of the correlation function ($`m=297`$ in the present application), discretize the time-delay axis $`t`$ in $`n`$ bins $`t_i`$ ($`12.5<t_i<12.5`$s, $`n=150`$ in our case), and write $`a_i^j=𝐚^j(t_i)`$. We treat the $`n`$ time bins as the variables whose statistical properties are to be found by PCA. Since neighboring time bins are highly correlated, the number of degrees of freedom is much less than $`n`$. The degree of correlation of the time bins with each other contains information on the actual number of degrees of freedom in the signal. We measure the correlation between time bins $`i`$ and $`k`$ by the covariance matrix $`C_{ik}`$:

$$C_{ik}=\underset{j}{}(a_i^j\overline{a}_i)(a_k^j\overline{a}_k),$$
(1)
where $`\overline{𝐚}=_j𝐚^j/m`$ is the average of $`𝐚`$ in the set. The principal components of $`𝐚`$ are then conventionally defined as those linear combinations of the columns of $`a_i^j`$ that diagonalize the covariance matrix. \[In many PCA applications, the correlation matrix is used instead of the covariance matrix. This is important if the variables measured are different physical quantities. For our case, this makes practically no difference, since the variables $`a_i`$ are the discretized values of a single continuous function.\] If the coefficients of these linear combinations are denoted by the matrix $`c_{ik}`$, the principal components $`p_k`$ are given by

$$p_k^j=\underset{i}{}c_{ki}a_i^j,$$
(2)
and the samples themselves are given, in terms of the principal components, by the inverse of $`c`$,

$$a_i^j=\underset{k}{}c_{ik}^1p_k^j.$$
(3)
These components are somewhat abstract, since a linear combination of different time bins is not a physically very meaningful quantity in our application. A more useful representation of the principal components is found by noting that the quantity

$$a_{ki}^j=c_{ik}^1p_k^j$$
(4)
describes the amplitude of principal component $`k`$ in the $`i`$-th time bin of the $`j`$-th sample. The average of this amplitude over the set of samples,

$$P_{ki}\underset{j}{}a_{ki}^j/m,$$
(5)
describes the principal components in the ‘time domain’, i.e. as a decomposition of our average cross correlation function. We use this representation for the rest of our analysis. The average of the cross correlation function over the set, $`\overline{𝐚}`$, is the sum of these principal components:

$$\overline{a}_i=\underset{k}{}P_{ki}.$$
(6)
### 4.2 Results

The covariance matrix is singular if the number of time bins $`n`$ is larger than the number of samples $`m`$. This limits the time resolution that can be used for the discretization of our cross correlation functions to 300 points111This limitation can be circumvented by data-compression of the samples. A convenient compression is obtained by applying a wavelet-transform to the samples. PCA is then applied on the first $`n_0m`$ wavelet coefficients, and the resulting principal components transformed back to the time domain by the inverse wavelet transform. For $`n_0=n`$ this process is loss-free and yields identical results. The results depend somewhat on the range in $`t`$ used for the analysis. Most of the signal of interest appears between $`t=\pm 5`$s. Using this range, and bins of 30ms yields the principal components shown in Fig. 5.

The amplitude $`A`$ of components varies both in absolute value and sign. For the first component we find that the average is larger than the variation, but for the higher components the average is generally small compared to the variation. Instead of showing the components at their average amplitude $`\overline{A}`$, we have therefore plotted them at a scale that takes into account both the average and the variation. The value used for this, somewhat arbitrarily, is $`\stackrel{~}{A}=(\overline{A}^2+\sigma _A^2)^{1/2}`$, where $`\sigma _A`$ is the rms deviation of $`A`$.

The first two components together account for 67% of the variance between the samples of the cross correlation function, if a time-lag interval $`(\mathrm{\Delta }t<t<\mathrm{\Delta }t)`$ with $`\mathrm{\Delta }t=5`$ s is used. The first component appears quite similar to the average cross correlation. The second component peaks around the ‘precognition dip’, though it has the opposite sign and extends somewhat to positive lag as well. If the time interval used is shorter, the first two components account for a larger fraction of the variation, and less for a longer interval. For $`\mathrm{\Delta }t=2s`$ the first two components account for 86% of the variance, for $`\mathrm{\Delta }t=12`$ s only 44%. This is understandable if the cross correlation consists of a signal of relatively stable shape concentrated around $`t=0`$, plus an unsystematically varying component that is spread out over a wider range in the time delay coordinate $`t`$.

The nature of the remaining components is not obvious at this point. Some of them, such as $`P_3`$ and $`P_5`$, appear to contain signals related to the sharp rise of the cross correlation at $`t=0`$. The PCA also interprets pure noise in the data in terms of additional components, and it is not a priori obvious which components are systematic and which may be just artefacts. A useful test to see which of the components is stable with respect to noise is to divide the set of samples in two halves and comparing the components derived from each of these subsets. This is shown in Fig. 6. The four curves shown result from two different ways of dividing the set in halves. In the first, the two subsets consist of the even numbered and the odd numbered samples of the (chronologically ordered) set. The difference between these two subsets measures the intrinsic variability of the cross correlation on a time scale of 25s (the duration of a sample). In the second division, the first half of the set is compared with the second half. This measures the variability on time scales of days.

As Fig. 6 shows, the components, especially the higher ones, are remarkably stable on both the short and the long time scales. From this one would be tempted to conclude that all six components shown are real. The process by which the components are derived in PCA, however, leaves the possibility open that some of the components are artefacts. This is because the components are constructed as orthogonal directions in the abstract $`n`$dimensional space of variables. The directions of insignificant higher components (constructed out of noise) are thus biased in a systematic way by the directions of the significant lower components.

The only way to test which of the components is real in this sense is by comparison of the results with a synthetic data set.

#### 4.2.1 Tests with synthetic data

To test the significance of the components, we have constructed synthetic data and applied PCA on this data in the same way as was done with the real data. By varying the input data, the resulting PCs can be adjusted for agreement. Since many of the components found are suspected to be spurious, the aim of the excercise was to find the minimum number of independently varying signals needed to reproduce the components found by PCA of the real data. It is not neccessary to start the synthetic process with time series of the count rates, since the questions to be answered concern the behavior of principal component analysis as applied to our cross-correlation functions. Thus, the synthetic data constructed consist of samples of the cross correlation function. Each of the samples consists of a linear combination of a number of adjustable functions of the time delay coordinate $`t`$, here called ‘signals’. The shape of these signals is the same in all samples, but their amplitudes differ. In addition to these signals, we add a random noise component that is different in each sample.

The amplitudes of the signals are assigned randomly between the samples, normally distributed with specified means and deviations. The noise component consists of white noise, low-pass filtered by boxcar averaging with a width $`\delta t`$. The parameters of the synthetic data are thus: the shapes of the signals, the means $`A_1`$ and $`A_2`$ and deviations $`\sigma _1`$, $`\sigma _2`$ of their amplitudes, the amplitude $`A_\mathrm{n}`$ of the noise component, and $`\delta t`$.

PCA on a set of 297 of these synthetic cross correlations yields the components shown in Fig. 8, for a case with two input signals plus a noise component. The shapes of the two input signals were adjusted to fit the shapes of the output components to those of the real data (Fig. 5). The final result for these input signals is shown in fig. 7

The amplitude of first two PCs of the simulated data are related directly to the average amplitudes of the two input signals. The amplitudes of PCs 3 and higher, however, turn out to be proportional to the amplitude of the boxcar-filtered noise component instead, showing that they are artefacts of the PCA. PCA somehow constructs these higher components out of the noise component. The shapes of these higher components, however, depend sensitively on the shape of the two input signals. The input signals can be determined rather accurately by the fitting process. Two signals plus a noise component turned out to be the minimum number needed. Tests with only one signal plus noise have also been done. These are able to reproduce the first two PCA components, but leave most of the detail in the higher components unexplained.

Except for the first, none of the components found by the PCA resembles the input signals very well. It is thus clear that PCA without comparison with a simulation using synthetic data can not be used to identify the statistically independent components reliably in the present problem. Unfortunately, the forward modeling process was quite cumbersome, owing to the large number of parameters to be adjusted, and is not easily automated.

From this test we conclude that just two independent signal components are needed to explain the observed variability of the cross correlation, apart from an aditional unsystematic noise component. The interpretation of these components presents some puzzles. One might have hoped that the actual independent signals would be the ‘spike’ and the ‘precognition dip’, but this turned out not to be the case. Extensive attempts were made to reproduce the PCA components with input signals resembling a pure dip before $`t=0`$ and a spike after $`t=0`$, but no satisfactory fit was found. When the shape of the signals was relaxed away from this prejudice, iteration converged onto the signals shown in Fig. 7. Their shapes could be determined rather accurately by this fitting process.

Both signal components consist of a dip followed by a peak. The second signal component resembles a compressed (in time) version of the first, but with the peak being somewhat less pronounced relative to the dip. A possible interpretation is thus that the intrinsic shape of the optical response is a combination of a dip followed by a sharp peak, while the time scale of this pattern varies, by a factor of 3 or so.

## 5 Power spectra

Information complementary to the cross correlations is given by the power spectra. The average X-ray and optical power spectra are shown in Fig. 9. The broad peak near 0.1 Hz in the X-ray power spectrum is associated with a QPO around this frequency. This QPO has been found in most of the X-ray data from KV UMa (Revnivtsev et al. 2000, Wood et al. 2000), and with high amplitude also in the optical light curve (Haswell et al. 2000, Pavlenko et al. 2001). In our data its amplitude was significantly smaller than in these earlier observations, and undetectable in the optical. Its frequency varied a bit during our observations, with some preference for frequencies of 0.08 and 0.2 Hz.

The variation of the power spectrum can be studied with dynamical power spectra; such studies have been reported for many sources (cf. van der Klis 2000). In the spirit of our analysis of the cross correlation, the variability in the power spectra can also by studied by PCA. Treating the 297 individual power spectra in the same way as the cross correlations yields components shown in Fig. 10.

The first component of the X-ray power spectrum is rather similar to the average, with a broad peak around 0.1 Hz. The second component shows a dip at 0.1 Hz and peaks at 0.08 Hz and 0.2 Hz. This could be explained by the presence of QPOs at frequencies of 0.08, 0.1 and 0.2 Hz, with the amplitude of the 0.1 Hz peak decreasing when the others are increasing. Tests with synthetic data show, however, that PCA also produces such a pattern if there is just a single a peak at a somewhat variable frequency. This seems a more likely interpretation in view of what is known about QPOs in similar sources.

The optical power spectra are affected by telescope vibrations at frequencies near 10 and 40 Hz. As explained above, these vibrations were caused by the rather strong wind loading of the telescope on several of the nights. The optical variability is also contaminated with seeing effects, which have frequencies in the range range of interest here. Since the X-O cross correlation reaches maximum values around 40% it is unlikely that seeing effects dominate the power spectrum, but its precise shape can not be reliably deduced from our observations.

With this caveat, the optical power spectrum shown in Fig. 9 is much broader than the X-ray power spectrum, and there is no clear evidence for the QPO seen in the X-rays. It thus seems unlikely that a QPO contributes strongly to the X-O cross correlations. The effect of the X-ray QPO’s on our cross correlations may well be just a a random noise contribution.

## 6 Conclusions and discussion

We have studied the correlation of X-ray and optical properties of the black hole transient KV UMa. Thanks to the length of the dataset, the statistical properties of this correlation could be studied in detail. By analysis of the time series themselves, we find that the optical response to an X-ray peak consists of a broad dip followed by a sharper peak. By cross-correlating the two time series the same pattern is found, but with a much higher signal-to-noise ratio. From the cross correlation we find that the onset of the peaks in the optical response coincides with the X-ray maximum to within 30 ms, and their initial rise is very fast. The minimum of the dip in the optical precedes the X-rays by about 2 seconds.

The agreement of the pattern of the cross-correlation with that found directly in the time series establishes, on the one hand, that the dip and the delayed peak are both properties of the optical light curve rather than the X-rays, and on the other hand that they are not due to artefacts of background subtraction or the cross-correlation process.

The shape of the cross correlation was found to be highly variable on time scales from minutes to days. A principal component analysis (PCA) of this variablity was found to be of some value in studying its statistical properties. It turned out that the principal components themselves, except for the first, are not very useful since their shape is highly affected by artefacts arising from a strong, randomly variable additive noise component in the cross correlation.

We find that the PCA components can be used, however, as an intermediate representation of the data, which can then be interpreted by a forward modeling process. The result of such a simulation shows that there are two systematic components which vary statistically independently. Both have the shape of a broad dip followed by a sharper peak, but their duration differs by a factor of about 3.

The most puzzling aspect of the optical variations is of course the ‘precognition dips’. If the optical light is due to reprocessing of X-rays, then an optical signal preceding the X-rays would have to be due to something in the accretion flow preceding an X-ray event. In principle, this could happen if the increase in mass flux which somewhat later is observed as an increase of the X-ray flux, starts at some distance from the hole, such that the accretion time scale from that distance is of the order of a few seconds. For typical accretion disk models, this would be a distance intermediate between the X-ray emitting region and the optical reprocessing regions in the outer disk. If this event is associated with a local thickening of the accretion disk, and the X-ray emitting region is small enough, such a ‘bump’ on the disk surface might be sufficient to shield the outer disk temporarily from the X-rays, leading to a dip in the optical emission.

While this possibility may be attractive in the reprocessing model, reprocessing itself has difficulties explaining the shape of the optical response. The sudden change in slope at zero lag turns out to be impossible to reproduce with any positive definite disk reprocessing function. The main reason for this is that the sudden change in slope requires much more variability on very short time scales ($`<100`$ ms) than the X-ray light curve actually has. A second difficulty is that the optical-to-X-ray flux ratio in KV UMa was quite high, or the order 10%. If, as the spectral energy distribution of McClintock et al. suggests, the UV-EUV emission is part of the same process producing the visible light, the reprocessed flux would have to be as large as the X-rays themselves.

From these difficulties it appears that alternatives may have to be explored to explain the visible light in black hole accreters like KV UMa and GX 339-4. In a future paper we will explore the possibility that the visible light is produced as thermal cyclo-synchrotron emission in the inner 20000 km of the accretion disk, as first suggested by Fabian et al. (1982) for GX 339-4 (see also Merloni et al. 2001 for an interpretation of the spectrum of KV UMa). If this interpretation is correct, it suggests a new and potentially very powerful means of probing the inner accretion flow by means of optical, and especiallly coordinated X-ray and optical observations.

###### Acknowledgements.

We thank the XTE time allocation team and Skinakas observatory for their support and flexibility in making these observations possible, Dr. C. Straubmeier and F. Schrey, who carried out the optical observations and Dr. T. Belloni, who provided the XTE count rate data. We thank Dr. Christian Motch for his comments on an earlier version of this text. This work was done with support from the European Commission under grant no ERB-FMRX-CT98-0195.


text:
# Edge-on boxy profiles in non-barred disc galaxies

## 1 introduction

A significant number of edge-on disc galaxies display in their inner regions a ‘boxy’ or ‘peanut’-shaped (hereafter b/p) structure. It is believed (Lütticke, Dettmar & Pohlen 2000), that more than 45% of disc galaxies have this kind of edge-on profile. Typical examples are NGC 2424, NGC 6771, NGC 5746, IC 4767, Hickson 87a and the Milky Way. Several authors (Combes & Sanders 1981; Pfenniger 1984, 1985; Combes, Debbasch, Friedli et al. 1990, Pfenniger & Friedli 1991, Raha, Sellwood J., James et al. 1991, Kuijken & Merrifield 1995, Bureau & Freeman 1999) have related these profiles to the presence of a strong bar. The tangential force in these bars is typically of the order of 25% of the axisymmetric one (Combes & Sanders 1981).

Patsis & Grosbøl (1996) have shown that b/p orbital profiles appear also in cases with a spiral instead of a bar perturbation. Recently Athanassoula (unpublished) made a large number of $`N`$-body simulations to follow the formation and evolution of bars in isolated bar-unstable discs. Some of the models show in their edge-on views a conspicuous b/p morphology, while their face-on views show clearly that they are non-barred, and even in some cases almost axisymmetric. It is clear that, at least in these cases, the b/p morphology is due to the internal dynamics of the self-consistent model and not to merging phenomena.

In this paper we will first describe a particularly illustrative simulation (section 2). We then use orbital theory to understand the dynamics of b/p structures. We use a purely axisymmetric potential consisting of a disc and a halo component (section 3), in order to identify the orbits that constitute the backbone of the b/p structure (section 4). The edge-on profiles are discussed in section 4. We discuss our results in section 6 and we stress the fact that it is the existence of vertical resonances per se and not the kind of perturbation that gives rise to the appearance of b/p morphologies.

## 2 The $`N`$-body model

Athanassoula’s $`N`$-body stellar model starts with initial conditions created by the method of Hernquist (1993) and consists of an exponential disc with a $`sech^2`$ vertical dependence, and a halo profile proposed by Hernquist (1993). The detailed description of the disc and halo initial density distributions can be found in Athanassoula & Misiriotis (2002; equations (1) and (3)).

In computer units the mass of the disc is taken $`M_d=1`$, the halo mass $`M_h=5`$, the disc scale length $`h=1`$, the scale height of the disc $`z_0=0.2`$ and the halo scale length $`\gamma =5`$. The rest of the parameters of the model in equations (1) and (3) of Athanassoula & Misiriotis (2002) are as described in Section 2 of the above mentioned paper. In the simulation we present here we can assume the unit of mass to be $`5\times 10^{10}M_{\mathrm{}}`$, the length unit 3.5 kpc, the unit of velocity 248 km s-1 and the time unit $`1.4\times 10^7`$ yr. The simulation has 200000 particles in the disc and the halo is live and composed of 931206 particles. The Toomre parameter Q is 1.8 (Toomre 1964). We underline the lack of an explicit bulge component. The simulation was carried out on a Marseille Observatory GRAPE-5 system using a tree-code similar to the one described in Athanassoula, Bosma, Lambert et al. (1998). For more details we again refer the reader to Athanassoula & Misiriotis (2002).

A characteristic snapshot of the simulation after 8.4 Gyr from the start is given in Fig. 1. The upper panel gives the circular velocity curve, which is rising and reaches approximately 0.68 (i.e. $``$ 170 km s-1) at 5 unit lengths (i.e. 17.5 kpc). The dashed and dotted lines give the disc and the halo contribution, respectively. The second, third and fourth panel from the top give the isodensities of the disc particles projected face-on, side-on and end-on respectively. The lower panel gives the dot-plot of the particles in the ($`x,y`$) plane. The side of the box for the face-on views is 10 units (i.e. 35 kpc), so the height of the boxes for the edge-on views 3.33 units (i.e. $``$ 11.65 kpc). From the face-on view it is evident that the model essentially does not have a bar. One can speak about a weak overall oval distortion. Such a weak perturbation can be detected in a large fraction of disc galaxies as an m=2 component of low amplitude. Nevertheless, both side-on and end-on views are obviously boxy. The fact that the relative extent of the b/p feature is about the same when projected on the horizontal axes reflects the fact that the isocontours of the density in the face-on view are nearly round.

## 3 The orbital model

For our orbital calculations we will use a general axisymmetric potential consisting of a disk and a halo component. We adopt a Miyamoto and Nagai (1975) disc potential, $`\mathrm{\Phi }_D`$, which in cylindrical coordinates has the form:

$$\mathrm{\Phi }_D(r,z)=\frac{GM_D}{\sqrt{r^2+[a+\sqrt{z^2+b^2}]^2}}$$
(1)
In the above the parameter $`M_D`$ refers to the disc mass, $`a`$ and $`b`$ are the horizontal and vertical scale lengths respectively, and $`r`$ and $`z`$ are the cylindrical coordinates. The halo potential, $`\mathrm{\Phi }_H`$, is given by

$$\mathrm{\Phi }_H(r,z)=\frac{v_H^2}{2}\mathrm{ln}\left(1+\frac{1}{r_c^2}(r^2+z^2)\right)$$
(2)
where $`v_H`$ is the limiting circular velocity as $`r\mathrm{}`$, and $`r_c`$ is the core radius of the halo. The total potential used for the orbital calculations is of the form: $`\mathrm{\Phi }=\mathrm{\Phi }_D+\mathrm{\Phi }_H`$, and the adopted values of the parameters are $`M_D=6\times 10^{10}M_{\mathrm{}}`$, $`a=3`$kpc, $`b=1.5`$kpc, $`r_c=18`$kpc and $`v_H=176.8`$km s-1. The rotation curve (Fig. 2) reproduces fairly well that of the $`N`$-body model presented in the previous section.

One can clearly see that we have a model with a maximum disc, as in the case of the $`N`$-body model. We remind that 5 length units in the $`N`$-body model correspond to 17.5 kpc, and its velocity unit is 248 km s-1.

## 4 Orbits and orbital behavior

The calculations have been done in cartesian coordinates in a frame of reference rotating around the $`z`$-axis. Thus, the Hamiltonian governing the motion of the test particles is:

$$H=1/2(\dot{x}^2+\dot{y}^2+\dot{z}^2)+\mathrm{\Phi }(x,y,z)1/2\mathrm{\Omega }_p^2(x^2+y^2),$$
(3)
where $`\mathrm{\Omega }_p`$ is the pattern speed. In the following we will denote by $`E_j`$ the numerical value of $`H`$ and adopt for the pattern speed the value $`\mathrm{\Omega }_p=11.9`$ km s-1 kpc-1, which places corotation approximately at 14 kpc.

On the equatorial plane we can define the radial resonances between the epicyclic frequency and the angular velocity in the rotating frame. In a 3D model we can also define vertical resonances, which involve the vertical instead of the epicyclic frequency (for definitions see e.g. Binney & Tremaine (1987). The first vertical resonance in our model is the 3:1.

The main family in our model is the family of direct circular periodic orbits on the $`z`$=0 plane. These orbits, in the presence of a spiral or barred perturbation become ellipses, the well known x1 orbits (see e.g. Contopoulos and Grosbøl 1986, 1989). By analogy, we will call this family x1 also in our axisymmetric model. In 2D models the orbits of this family support the spiral or bar structure. The stability of a periodic orbit is characterized by the behaviour of two indices $`b_1`$ and $`b_2`$. In this study, the stability index $`b_1`$ is associated to the motion perpendicular to the equatorial plane, while $`b_2`$ is associated to radial perturbations. A family is stable if both stability indices $`b_i`$ are $`2<b_i<2`$ (Hadjidemetriou 1975). For more details on the stability of families of periodic orbits in 3D systems the reader may refer to Contopoulos and Magnenat (1985).

The most important families of periodic orbits for the dynamics of a 3D disc galaxy are the central family and those bifurcated from it at the vertical $`n:1`$ resonances, where $`n`$ is a small integer. At these resonances, in the axisymmetric case, index $`b_1`$ of family x1 becomes tangent to the $`b=2`$ axis. The bifurcated families come actually in pairs111The two families, which bifurcate at a tangency of a stability index with the $`2`$ axis in the axisymmetric case, have a stable and an unstable counterpart in the non-axisymmetric case and are in this case marginally stable because they always have one of their stability indices on the $`b=2`$ axis, while the other remains always between $`2`$ and $`2`$. The variation of the stability indices $`b_1`$ and $`b_2`$ of the x1 family, as well as those of the families bifurcated at the vertical resonances 3:1, 4:1, 5:1, 6:1 and 7:1 are given in Fig. 3 as a function of the Jacobian $`E_j`$. The vertical black arrows indicate the points where index $`b_1`$ of x1 becomes tangent to the $`b=2`$ axis at the vertical resonances. Horizontal arrows point to the index of each bifurcating family which oscillates between $`2`$ and 2. They are given always close to the value of the Hamiltonian at which the family bifurcates from the x1. The introduction of a barred or spiral perturbation in the model brings in the system a stable and an unstable family. For the stable one the index which in the axisymmetric case was lying on the $`b=2`$ axis now becomes absolutely smaller than 2, so the family remains stable over a large radial region (Patsis & Grosbøl 1996; Skokos, Patsis, Athanassoula 2002a,b). The vertical white arrows indicate the tangencies of the index $`b_2`$ of the x1 family with the $`b=2`$ axis at the radial resonances.

## 5 The edge-on profile

Most stars in a galaxy will move along non-periodic orbits trapped around stable periodic orbits (Poincaré 1892). Thus, the topology of the main families of periodic orbits will determine the basic features in the galaxy. In the present case the main families are the central family and the families bifurcated at the vertical resonances, which exist as two branches symmetric with respect to the equatorial ($`z=0`$) plane. The response density profile is a weighted average of the individual orbit contributions. As a weight, we use the disc density, $`\rho _D(r,z)`$, where $`\rho _D`$ is the density corresponding to the Miyamoto disc calculated at the position $`(r,z)`$, where $`r`$ is the radius of the circular equatorial-plane orbit which has the same $`E_j`$ as the orbit to be weighted and $`z`$ is the mean vertical distance of the orbit from the equatorial plane. We have adopted the density of the luminous (disc) component because the orbital profiles will have to reproduce the light distribution, when we compare to real galaxies, or the disc component of the $`N`$-body simulations.

Since our model is axisymmetric we can use arbitrarily any two orthogonal axes on the equatorial plane in order to calculate periodic orbits. This means that a periodic orbit rotated around the axis of symmetry ($`z`$-axis) is also a periodic orbit of the system. The projections of these periodic orbits on a given axis will always be confined within certain limits determined by a maximum length and a maximum height. We also note that the orbits of the two families that bifurcate from the central family at the tangencies of the $`b_1`$ index with the $`b=2`$ axis (Fig. 3), at a given $`E_j`$ value, are topologically similar.

We have calculated at each energy ($`E_j`$) only one periodic orbit per family along our $`x`$-axis. We started with the orbit which has the same radius with the circular orbit at the $`E_j`$ where the family is bifurcated, and we followed the evolution of the family by finding the orbits along the specific $`x`$-axis. These orbits have been used for constructing the profiles in Fig. 4. Nevertheless, at each energy, one can find, by rotation, an infinite number of representatives of the same family, due to the axisymmetric nature of the potential. It is as viewing a given orbit from all possible view angles by rotating our point of view around the $`z`$-axis, while staying always on the equatorial plane.

In the $`(y,z)`$ profile, created by (marginally) stable orbits of the families considered in Fig. 3 (Fig. 4a), we observe a ‘stair-type’ structure with a central boxy region. This central region is formed by the orbits of the two branches (symmetric with respect to the equatorial plane) of the 3D family bifurcated at the vertical 3:1 resonance. By applying a smoothing filter on the image with the weighted orbits we obtain a blurred profile, which, to a first approximation, can be considered as the profile for the density distribution of the model. This blurred profile is shown in Fig. 4b. By considering all possible orbits at each energy, the area inside the rectangle drawn with a dashed line in Fig. 4b will be filled. The $`(x,z)`$ profile, when considering all orbits, will be identical to the $`(y,z)`$ one, as expected.

In Fig. 4b we illustrate the resulting morphology by using isocontours. Clearly the profiles of the orbital model (Fig. 4) and the profiles of the $`N`$-body simulation (Fig. 1) have many similarities. First of all we have a boxy structure confined in radii less than 2 length units of the $`N`$-body simulation, which correspond to 7 kpc. This is practically the radius in Fig. 4 inside which we find the boxy structure of the orbital model. Also the ‘stair-type’ edge-on profiles of the weighted orbits is in good agreement with the outer parts of the edge-on profiles of the $`N`$-body model.

## 6 Discussion

In the present study we propose a mechanism, which can account for boxy structures in the edge-on profiles of non-barred disc galaxies. In the example we present here, it is the presence of the vertical 3:1 resonance which introduces this structure in the system. It is, however, not necessary to have a particular vertical resonance in order to form a b/p profile. Furthermore, since radial and vertical resonances can be defined even in an axisymmetric case, one can find even in the axisymmetric model the bifurcating orbits that give rise to these features. In fact, the only two necessary ingredients are:

\- the presence of a vertical $`n:1`$ resonance, where $`n`$ a small integer, so that a new family, bifurcating at this point, is introduced in the system

\- the bifurcated family should be stable over a sufficiently large $`E_j`$ interval and should trap a sufficient number of regular orbits around it.

If these two conditions are fulfilled, the ‘stair-type’ orbital profile follows naturally, because the successive 3D families lower their mean heights as the energy at which they are bifurcated from x1 increases. In addition, as energy increases, the successive orbits of a bifurcating 3D family increase their mean spherical radius by growing more in $`z`$ than in their cylindrical radius and, beyond a critical energy, by growing practically only in the vertical direction. Since their cylindrical radius – or extent along the equatorial plane – is thus limited, while their vertical extent is large, at least for large values of the Jacobi energy $`E_J`$, they contribute to a boxy profile. As we have shown here, it is not necessary to have a strong perturbation in order for the vertical extent to be important. Thus the boxy feature can be strong and clearly defined, even in a purely axisymmetric model.

If we introduce a perturbation, the topology of the relevant orbital families remains the same and they will be stable rather than marginally stable (Patsis & Grosbøl 1996; Skokos et al. 2002a,b). There will, however, be one morphological difference, due to the fact that in the non-axisymmetric case the rotational symmetry is broken. In this case it will not be possible to have orbits of any desired azimuthal orientation, and as a result we will have a ‘peanut-shaped’ profile, at least for a range of viewing angles, instead of a boxy one.

Lack of a boxy structure in a model indicates one of the following three possibilities: Either vertical resonances with small $`n`$:1 do not exist in the system. Or the family bifurcated at the vertical resonance with the lowest $`E_j`$ value has too large unstable parts. Or, for reasons that could be linked to the formation history of the galaxy, too few stars are on orbits trapped around the stable periodic $`n`$:1 orbits. The conditions needed for building a b/p profile may be favored by the presence of the bar, but the bar per se is not the reason that edge-on disc galaxies have boxy profiles. Of course our mechanism relies on the existence of a pattern speed which does not vary much with time. In other words it implicitly assumes the existence, or past existence, of some non-axisymmetric feature, albeit of perhaps infinitesimal amplitude.

Bureau & Athanassoula (1999) and Athanassoula & Bureau (1999) developed diagnostics to detect the presence and orientation of a bar in edge-on disc galaxies. They detected in most of the peanut shaped edge-on galaxies in Bureau & Freeman (1999) the signature of a ‘x2-flow’, in the position-velocity diagrams. This was taken as the manifestation of the presence of a bar. However, in a few cases, this feature was absent. Athanassoula & Bureau accounted the lack of such a feature either to a lack of an ILR or to a lack of emitting gas around the ILR region. The present study adds a third possibility, namely that the galaxy is not barred.

## Acknowledgments

We acknowledge fruitful discussions and very useful comments by Prof. G. Contopoulos. We also thank the anonymous referee for her/his comments which improved the paper. This work has been supported by the Research Committee of the Academy of Athens. E. Athanassoula would also like to thank the IGRAP, the Region PACA, the INSU/CNRS and the University of Aix-Marseille I for funds to develop the GRAPE computing facilities used for the simulations discussed in this paper. P.A. Patsis and Ch. Skokos thank the Laboratoire d’Astrophysique de Marseille, for an invitation. Ch. Skokos was supported by the “Karatheodory” fellowship No 2794 of the University of Patras.