3 The cluster evolution

 We follow the cluster evolution starting from a clean stationary state (no dissipative perturbation), which corresponds to a spheroidal stellar system gravitationally coupled with the central black hole (Bahcall & Wolf 1976). Hence, number density n₀(a) of the satellites is assumed proportional to a^1/4, while the initial distribution in eccentricities and inclinations conforms with the condition of a randomly generated sample. Common parameters of the computations are the mass $M=10^8~M_{\odot}$ of the central body, and the orbiters' column density $\Sigma _{\ast }=\Sigma _{\odot }$. Now, the perturbation effects are switched on, depending on the prescribed disc type and accretion rate (typically $\mu_{_{\rm {}E}}=1$, but we also examined other situations and less conventional values of the above parameters; cf. Subr 2001).

As a consequence of the dissipative action, the system departs from the initial configuration. We concentrate our main attention on the range of radii , well below the outer edge of the disc $r_{\rm {d}}\approx10^4 r_{\rm {}g}$ (where $\Sigma_{{\rm d}}$ becomes negligible). We find a new quasi-stationary, modified cluster distribution developing in that region. On the outer boundary, fresh satellites could be inserted in the system from an external reservoir in order to maintain the steady flow towards the center, but in this work we alternatively consider a closed system with a large but finite number of its members, N₀, which are not replenished. The actual value of N₀ is not very important in the present work because we neglect gravitational interaction among the satellites themselves (N₀ will play a role when two-body collisions are taken into account).

3.1 A role of the disc model

Several different models of the disc were adopted. We specify them in terms of their corresponding density profiles, $\Sigma_{{\rm d}}(r)$, and vertical thicknesses, h(r):

(i)

Radiation pressure dominated, geometrically thin disc (Shakura-Syunaev $\alpha$-viscosity model; Frank et al. 1992);

(ii)

Gas pressure dominated standard disc with opacity due to electron scattering;

(iii)

The same as in (ii) but with opacity dominated by free-free scattering;

(iv)

Gravitationally unstable discs in the region of a self-gravitating disc beyond the Toomre radius; see Sect. (2.3) of Collin & Huré (1999) for analysis of the marginally unstable ($\zeta=5$) self-gravitating disc with solar metallicity and large optical thickness;

(v)

The same as in (iv) but for optically thin medium with zero metallicity.

The choice of different models is dictated by their applicability for the description of accretion regimes relevant under different conditions, namely, the cases (iv)-(v) are appropriate for black-hole discs at large ( $10^3r_{\rm {}g}$) radii. One could consider other models relevant in various situations but the adopted examples cover qualitatively different evolutions which are frequently encountered. In particular, the above models comprise a range of $\Sigma_{{\rm d}}/\Sigma_{\ast}$ dependencies (dimensionless factor determining the evolutionary time-scales) and different $(h/r)(\Sigma_{{\rm d}}/\Sigma_{\ast})=\epsilon/\tau$ whose radial dependence determines area in the disc where the gap is formed.

First, we carried out simulations while sticking to one of these models in the whole range of radii from the centre up to the outer edge, even if it extends somewhat out of the zone of validity of the particular model. This helps us to avoid additional complexities which could stem from a complicated prescription for the disc properties. For example, one expects a different rate of radial drift across boundaries where the disc properties are switched, but we refrain from this complication for the moment. This effect will be observed later on: to achieve a more realistic description, a further step will be carried out by joining different regions of the disc in a unified model. In that case, the radiation pressure supported standard disc might be limited to the innermost part of the system while gravitationally unstable part takes over at distances of the order of a few hundreds  $r_{\rm {g}}$.

Parameters defining radial migration in the models (i)-(v) were given in Table 1. Here we only remark that the factor B determines the pace in absolute time units with which the evolution proceeds, while power-law indices q_iinfluence the form of the satellites number density distribution n(a).

Figure 3: Number densities nt(a) of the cluster. Left panel: quasi-stationary states for $M_{\ast}=1~M_{\odot}$ (solid line), $3~M_{\odot}$ (dashed), and $10~M_{\odot}$ (dotted). Right panel: snapshots of nt(a) are shown for $M_{\ast}=1~M_{\odot}$ at t=106(solid), 108 (dashed), 1010 (dotted), and 1012 (dash-dotted). In all cases, $\Sigma _{\ast }=\Sigma _{\odot }$ (solar-type compactness of the satellite). On the vertical axis, units are arbitrary and scaled with respect to the initial profile $n_0(a)\,\propto \,a^{1/4}$. Units of GM/c2 are used for a on the horizontal axis.

3.2 Modified cluster distributions

 The main results are shown in Fig. 3, in which the disc model (iii) is adopted. This will serve as a representative test case (more examples are deferred to Appendix). It is particularly interesting to observe transition regions where the mode of the mutual satellite-disc interaction is changed. At the transition, the efficiency of radial transfer of the satellites is changed, and this effect produces sudden changes (either dips or concentrations) in n_t(a)profiles. If the average inclination of the satellites is small (the case of a flattened configuration), a ring is created in the place where the orbital decay is stalled.

The computations were launched from the initial distribution $n_0(a)\,\propto \,a^{1/4}$. In the left panel, Fig. 3a, new n_t(a) profiles are shown at a sufficiently late time, t=10¹², when a quasi-stationary state has been established. Three cases are plotted with increasing satellite masses: $M_{\ast}=1~M_{\odot}$, $3~M_{\odot}$, and $10~M_{\odot}$. After being inclined to the disc, low mass satellites proceed via density waves (e.g., $1~M_{\ast}$ case with $n_t(a)\,\propto\,a^{-1/2}$ at small a). On the other hand, very massive satellites develop a gap already very far from the centre and they proceed in this mode the whole way down ( $10~M_{\ast}$ case with $n_t(a)\,\propto\,a^{1/4}$, which is incidentally parallel with the initial n₀ distribution). Intermediate satellites show an evident dip in n_t(a): its formation is connected with the gap in the disc which arises at sufficiently small a. Exact location of the dip depends on the disc model too; for $3~M_{\odot}$ satellites it occurs at  .

The right panel, Fig. 3b, is complementary to the left one, showing temporal evolution of n_t(a) for one solar mass satellites. One notices two areas with different n_t(a) slopes in the graph: at small a the overall distribution is eventually dominated by satellites residing in the disc (corresponding to negative slope), while a whole mixture of satellite inclinations persists farther out (a positive slope). Transition between the two regions occurs at $a\approx10^2$ (its location moves gradually towards larger a as time proceeds; see the figure). A terminal dip is seen at large a (on the right of the graphs); it is caused by depletion of the sample in the simulation.

A realistic situation must involve further changes of the n_t(a) slope caused by varying the accretion mode (see Appendix and Fig. A.1 for details). In this respect, the form of quasi-stationary distributions is of particular interest (the persisting quasi-stationary profile of n_t(a) is given more detailed discussion in the next section, Sect. 3.3). We remark that the long-term evolution of the orbits is influenced by collisions with the disc especially if the satellite's cross-sectional area is large, e.g. a giant star or a solar-size body, but quite the opposite conclusion can be drawn for compact objects such as neutron stars and stellar-mass black holes for which the cross-section is much smaller, diminishing their collisional interaction with the disc. The whole cluster is modified by this kind of dissipative process in a selective manner, depending on $M_{\ast }$.

In Fig. 4a we present the resulting form of the distribution n_t=10¹²(a) for different masses of compact bodies. On the other hand, panel 4b concerns a cluster consisting of one solar-mass satellites, $M_{\ast}=1~M_{\odot}$, shown at different times (several consecutive moments are examined with a logarithmic spacing in t). Now, there is no need to consider separately that part of the distribution of trajectories inclined into the disc; indeed, no compact bodies reach the disc plane, even at late phases of the evolution. This implies that the orbits evolve predominantly by gravitational radiation and the resulting graphs come out almost identical for all the disc models under consideration. We are thus led to restrict ourselves to non-compact satellites which suffer relatively intense drag from the disc.

Figure 4: Left panel: the shape of the distribution nt(a) of a cluster containing compact orbiters of the masses $M_{\ast}=1.5~M_{\odot}$(solid line), $3~M_{\odot}$ (dashed) and $10~M_{\odot}$ (dotted) at t=1012. Right panel: subsequent phases of evolution for a cluster of neutron stars with $M_{\ast}=1.5~M_{\odot}$. The power-law is represented by a line in the log-log plane of nt (normalized by N0, the initial number of cluster members) vs. a (units of  $r_{\rm {g}}$). Again, a thin dotted line stands for the initial power-law distribution, n0. Other curves correspond, from top to bottom, to t=109, 1011, 1013, 1015, and 1017.

3.3 A quasi-stationary state

 Time-scales of the perturbed cluster evolution depend on the disc model and, as shown above, a number of other factors characterizing the cluster members - their typical mass, compactness and efficiency of the dissipative interaction with the disc medium. With solar-type (and less compact) satellites, the disc plays a dominant role in determining the form of n_t(a). The starting canonical distribution $n_0(a)\,\propto \,a^{1/4}$ (corresponding to $\rho\,\propto\,r^{-7/4}$in the notation of Bahcall & Wolf 1976) changes to n_t(a), where it stands almost frozen in a quasi-stationary state for a prolonged interval of time. This is true also in the case of initial profiles reasonably different from the canonical one. In particular, we verified that $n_0(a)={\rm {}const.}$converges to identical slopes as those shown in Fig. A.1 at $t\approx10^{12}$ (cf. the bottom row of plots in the Appendix).

Clear-cut results emerge provided the distribution is influenced by a single mode of star-disc interaction. The evolution time-scales are shorter in the region where substantial fraction of the satellites are already inclined into the disc plane, so that n_t(a) forgets about the exact initial distribution much earlier, at about -10⁹. In that case the final slope  $q_{{\rm f}}$ ranges from the initial value up to the limiting value q₄ of the corresponding model and the mode of migration (the values are listed in Table 1). The actual final value of  $q_{{\rm f}}$ reaches q₄ if two conditions are satisfied in the region where the power-law distribution is established: first, the flow of the satellites in this region continues in the disc plane (all bodies are inclined to the disc); second, the satellites inflow is faster than the changes of their number density at the outer boundary. These conditions are satisfied in the late stages until the sample is exhausted. In other words, the modified cluster distribution reflects the adopted disc model while it is rather insensitive to the initial form of n₀(a).

Figure 5: Left panel: a quasi-stationary nt(a) profile is plotted (by crosses) for the model (i), t=1013. A power-law approximation, $n_t(a)\,\propto \,a^{q_{{\rm f}}}$, is constructed in appropriate ranges of radius (solid lines). Suppression of the profile occurs at small a in consequence of quite efficient gravitational radiation losses near the center. Right panel: the slope evolves from the initial value, q=1/4, until the quasi-stationary state of $q\rightarrow {q_{{\rm f}}}$ is reached. Two snapshots are shown for the model (iii): t=1010 (upper profile), and t=1012 (lower profile). See the text for corresponding values of $q_{{\rm f}}$.

This behaviour is illustrated in Fig. 5 in terms of a power-law approximation to the computed distributions. Given the disc model, different power-law relations may be established in separate regions depending on the mode of star-disc interaction applicable at that distance. In the case (i) and for solar-type satellites, we find that the slope $q_{{\rm f}}=-4.9$ arises in the region (the evolution is driven by density waves, accumulating the satellites in an evident peak), while $q_{{\rm f}}=2.3$ suits for (i.e. beyond the transition point where the star-disc interaction is determined by the gap; here, as an example, $q_{{\rm f}}=2.3$ is to be compared with the corresponding $q_4=\mbox{$\,^{5}/_{2}$ }$ in the sixth row in Table 1). Furthermore, the slope is influenced by gravitational radiation which dominates the evolution on very small radii: $q_{{\rm f}}=2.8$ for $r\leq20$. The three regions are distinguished in the graph: gw (dominated by gravitational waves), dw (density waves), gap (gap formation). The large span of values acquired by $q_{{\rm f}}$ reflects the presence of ripples in n_t(a) which develop in certain areas of the disc.

The right panel of Fig. 5 shows how the n_t(a) profile is gradually adjusted as time proceeds. Here we plot the case (iii) with for illustration, and we find q=-0.71(t=10¹⁰, ), and q=-0.51 (t=10¹², ), respectively. On the vertical axis, one can read the fraction of satellites with semi-major axes falling in interval $\langle{a,a+{\rm {}d}a}\rangle$. One observes different slopes for more massive satellites, , because of gap formation (the slope reaches $q\,\dot=\,0.25$); cp. with Fig. 3 to see the effect of $M_{\ast }$. On the other hand, the situation gets simpler with the other models where such a transition does not occur, and a single power-law relation prevails at late stages for all . The emerging slopes are: (ii)  $q_{{\rm f}}=-0.79$, and (iv) 0.30, for which we found that the quasi-stationary state is established at , t=10¹². The case (v) develops more slowly and reaches the expected value $q_{{\rm f}}=-0.67$somewhat later, at time t=10¹³.