text:
# MASS SEGREGATION IN STAR CLUSTERS: ANALYTIC ESTIMATION OF THE TIMESCALE

## 1 INTRODUCTION

Within a star cluster, massive objects segregate into the cluster core as they lose kinetic energies to the less massive objects during approach toward energy equipartition. The result is an increase of the number density of the massive objects in the core. This phenomenon of mass segregation has long been known (Spitzer 1969; Binney & Tremaine 1987, p. 531) but has attracted new interests in recent years. The large number density induces successive mergers of those massive objects, which could lead to the formation of an exotic object such as a very massive black hole.

The important parameter of mass segregation is its timescale. In order to explain an observation in terms of mass segregation, the timescale has to be less than the age of the star cluster.

However, mass segregation has been studied only for specific cases using numerical methods, i.e., Fokker-Planck, Monte Carlo, and $`N`$-body simulations. The exceptional analytical works that are available are those of Spitzer (1969) and Tremaine, Ostriker, & Spitzer (1975). Spitzer (1969) considered an isothermal cluster of uniform density and estimated the timescale of mass segregation as the timescale of energy equipartition between stars of different masses. This model is too idealized for a star cluster, which is not isothermal or of uniform density. Tremaine et al. (1975) studied a decaying circular motion of an object in the singular isothermal sphere. This is originally a model for a motion of a globular cluster around the center of a galaxy, and it is not a model for random motions of stars in a cluster. Therefore, an analytical study of a more appropriate situation is desirable.

We analytically study a two-component cluster, which consists of two masses $`m_0`$ and $`m`$ ($`m_0<m`$). The major component is the less massive stars. Their total mass $`M_0`$ dominates over the total mass $`M`$ of the more massive stars, $`M_0M`$. Here the subscript 0 is used to indicate that the quantity is relevant to the major component. The cluster is spherically symmetric in all its properties and completely isotropic in the velocity space.

Since the timescale of mass segregation is expected to be less than the relaxation timescale (Spitzer 1969), the star cluster is regarded as almost collisionless. The individual stars move almost freely in the mean potential. The gravitational encounters occur only in a stochastic manner. Our study is accordingly based on the overall behaviors of the probability distribution functions (PDFs) of the major and minor components. They are approximated by PDFs of completely collisionless systems. The minor-component PDF is from a parameterized family of steady-state PDFs for which the degree of mass segregation is allowed to change. The major-component PDF is kept the same (§2). These PDFs are used to obtain the energy loss per unit time from the minor component due to encounters with the major component (§3). The energy loss rate yields the rate of change of the minor-component PDF, which in turn yields the timescale of mass segregation (§4, eq. ).

We compare our result with those for relevant timescales, including the pioneering work of Spitzer (1969). We also use our result to discuss the origin of a black hole with mass of $`10^3M_{\mathrm{}}`$ found in a young star cluster of the starburst galaxy M82 (§5).

## 2 PROBABILITY DISTRIBUTION FUNCTIONS

The major component, i.e., the less massive stars, is assumed to follow Plummer’s model (polytrope with $`n=5`$; Binney & Tremaine 1987, p. 223). This steady-state solution for a collisionless spherical system has a simple analytic form but well reproduces observations of some globular clusters. The PDF in the phase space is determined by two parameters, e.g., the total mass $`M_0`$ and the core radius $`r_\mathrm{c}`$:

$$f_0=\frac{32\sqrt{2}}{7\pi ^2}\frac{3M_0}{4\pi r_\mathrm{c}^3}\left(\frac{GM_0}{r_\mathrm{c}}\right)^5\left(\mathrm{\Psi }_0\frac{v^2}{2}\right)^{7/2},$$
(1)
with $`f_0=0`$ at $`v(2\mathrm{\Psi }_0)^{1/2}`$. Here $`G`$ is the gravitational constant and $`v`$ is the velocity. The relative potential $`\mathrm{\Psi }_0`$ is a function of the radius $`r`$ and defined as the negative of the gravitational potential $`\mathrm{\Phi }_00`$:

$$\mathrm{\Psi }_0(r)=\mathrm{\Phi }_0(r)=\frac{GM_0}{r_\mathrm{c}}\frac{1}{\left(1+r^2/r_\mathrm{c}^2\right)^{1/2}}.$$
(2)
The corresponding mass density $`\rho _0`$ is

$$\rho _0(r)=4\pi _0^{\mathrm{}}v^2f_0𝑑v=\frac{3M_0}{4\pi r_\mathrm{c}^3}\frac{1}{\left(1+r^2/r_\mathrm{c}^2\right)^{5/2}}.$$
(3)
Here the integral $`4\pi _0^{\mathrm{}}r^2\rho _0𝑑r`$ gives the total mass $`M_0`$. The squared velocity dispersion $`\sigma _0^2`$ is

$$\sigma _0^2(r)=\frac{4\pi _0^{\mathrm{}}v^4f_0𝑑v}{4\pi _0^{\mathrm{}}v^2f_0𝑑v}=\frac{GM_0}{2r_\mathrm{c}}\frac{1}{\left(1+r^2/r_\mathrm{c}^2\right)^{1/2}}.$$
(4)
The minor component, i.e., the more massive stars, is assumed to have no contribution to the gravitational field of the star cluster. Then the gravitational field is determined solely by the relative potential $`\mathrm{\Psi }_0`$ of the major component. According to the Jeans theorem for a steady-state collisionless spherical system (Binney & Tremaine 1987, p. 221), any function of $`\mathrm{\Psi }_0v^2/2`$ serves as the minor-component PDF. The PDF assumed in our study is from a function family that has a parameter $`\beta 1`$:

$$f_\beta =\frac{32\sqrt{2}}{7\pi ^2}F_1F_2\frac{3M}{4\pi r_\mathrm{c}^3}\left(\frac{GM_0}{r_\mathrm{c}}\right)^{(7\beta +3)/2}\left(\mathrm{\Psi }_0\frac{v^2}{2}\right)^{7\beta /2},$$
(5)
with $`f_\beta =0`$ at $`v(2\mathrm{\Psi }_0)^{1/2}`$. The factors $`F_1`$ and $`F_2`$ are functions of $`\beta `$ defined as

$$_0^{\mathrm{}}\frac{x^2dx}{\left(1+x^2\right)^{(7\beta +3)/4}}=\frac{1}{3}\frac{1}{F_1(\beta )},$$
(6)
and

$$_0^1x^2\left(1x^2\right)^{7\beta /2}𝑑x=\frac{7\pi }{512}\frac{1}{F_2(\beta )},$$
(7)
where $`F_1=F_2=1`$ at $`\beta =1`$. Their numerical values are shown in Table 1 and Figure 1. The corresponding mass density $`\rho _\beta `$ is

$$\rho _\beta (r)=F_1\frac{3M}{4\pi r_\mathrm{c}^3}\frac{1}{\left(1+r^2/r_\mathrm{c}^2\right)^{(7\beta +3)/4}},$$
(8)
with the total mass $`M`$. The squared velocity dispersion $`\sigma _\beta ^2`$ is

$$\sigma _\beta ^2(r)=\frac{12}{7\beta +5}\frac{GM_0}{2r_\mathrm{c}}\frac{1}{\left(1+r^2/r_\mathrm{c}^2\right)^{1/2}}\sigma _0^2(r).$$
(9)
Figure 2 shows the radial profiles of the mass density and squared velocity dispersion for $`\beta =1`$, 10, and 100. The total energy $`E_\beta `$ is obtained by integrating $`\rho _\beta (\sigma _\beta ^2+\mathrm{\Phi }_0)=\rho _\beta (\sigma _\beta ^2\mathrm{\Psi }_0)`$ over the entire volume:

$$E_\beta =\frac{3\pi }{32}\frac{2(7\beta 1)}{7\beta +5}F_1F_3\frac{GM_0M}{r_\mathrm{c}}.$$
(10)
The factor $`F_3`$ is a function of $`\beta `$ defined as

$$_0^{\mathrm{}}\frac{x^2dx}{\left(1+x^2\right)^{(7\beta +5)/4}}=\frac{\pi }{16}F_3(\beta ),$$
(11)
where $`F_3=1`$ at $`\beta =1`$. We define the initial state as $`f_{\beta =1}`$ so that the radial profile of mass density and the velocity dispersion are the same as those of the major component. With an increase of the $`\beta `$ value, the mass density of the minor component becomes more centrally concentrated than that of the major component. Simultaneously, the velocity dispersion of the minor component becomes smaller (Fig. 2). These behaviors represent mass segregation. The radial profile of velocity dispersion of the minor component remains the same as that of the major component. This property is essential to establishing the final state of mass segregation where no net energy exchange exists between the two components (§3). Thus, with a single parameter $`\beta 1`$, the function family $`f_\beta `$ provides the simplest model to trace the evolution of the minor component from the initial to final states of mass segregation. At an intermediate state, the minor-component PDF is set to be just one of the functions, $`f=f_\beta `$, or a linear combination of them, $`f=(M_\beta /M)f_\beta `$ with $`M=M_\beta `$ (§4).

## 3 ENERGY CHANGE RATE

Consider a single object in the minor component, which loses kinetic energy due to gravitational encounters with objects in the major component. If the velocity of the single object is $`v`$, the rate of change of kinetic energy per unit mass is obtained as

$$\frac{}{t}\left(\frac{v^2}{2}\right)=16\pi ^2G^2\mathrm{ln}\mathrm{\Lambda }\left(m_0_v^{\mathrm{}}v_0f_0𝑑v_0m_0^v\frac{v_0^2}{v}f_0𝑑v_0\right),$$
(12)
with $`f_0=f_0(\mathrm{\Psi }_0v_0^2/2)`$ (Binney & Tremaine 1987, eq. \[8-66\]). Here $``$ denotes an ensemble average. The cutoff factor $`\mathrm{\Lambda }`$ is the ratio of the maximum and minimum impact parameters, $`b_{\mathrm{max}}`$ and $`b_{\mathrm{min}}`$:

$$\mathrm{\Lambda }=\frac{b_{\mathrm{max}}}{b_{\mathrm{min}}}=\frac{b_{\mathrm{max}}v_{\mathrm{rel}}^2}{G(m+m_0)}\frac{m_0}{m+m_0}N_0.$$
(13)
Here $`v_{\mathrm{rel}}`$ is a typical relative velocity between objects in the cluster (Binney & Tremaine 1987, p. 511), and $`N_0=M_0/m_0`$ is the total number of objects in the major component. We have assumed $`v_{\mathrm{rel}}^2GM_0/b_{\mathrm{max}}=Gm_0N_0/b_{\mathrm{max}}`$.

By multiplying equation (12) by the minor-component PDF $`f`$ and integrating over the $`v`$ space, we obtain the rate of change of kinetic energy per unit volume of the minor component. For the PDF $`f=f_\beta (\mathrm{\Psi }_0v^2/2)`$, the result is

$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\rho \sigma ^2}{2}}\right)`$ $`=`$ $`{\displaystyle \frac{2^{19}\sqrt{2}}{49\pi }}{\displaystyle \frac{16!!}{19!!}}F_1F_2F_4{\displaystyle \frac{1}{\left(1+r^2/r_\mathrm{c}^2\right)^{(7\beta +12)/4}}}`$ (14)
$`\times G^2\mathrm{ln}\mathrm{\Lambda }\left({\displaystyle \frac{m}{7\beta +2}}{\displaystyle \frac{m_0}{9}}\right){\displaystyle \frac{3M_0}{4\pi r_\mathrm{c}^3}}{\displaystyle \frac{3M}{4\pi r_\mathrm{c}^3}}\left({\displaystyle \frac{GM_0}{r_\mathrm{c}}}\right)^{1/2}.`$
The factor $`F_4`$ is a function of $`\beta `$ defined as

$$_0^1x^2\left(1x^2\right)^{(7\beta +9)/2}𝑑x=\frac{16!!}{19!!}F_4(\beta ),$$
(15)
where $`F_4=1`$ at $`\beta =1`$. By integrating equation (14) over the entire volume, we obtain the rate of change of total energy $`E`$ of the minor component:

$`{\displaystyle \frac{dE}{dt}}`$ $`=`$ $`{\displaystyle \frac{2^{21}\sqrt{2\pi }}{3\times 5\times 7^2}}{\displaystyle \frac{16!!}{19!!}}{\displaystyle \frac{\mathrm{\Gamma }(13/4)}{\mathrm{\Gamma }(15/4)}}F_1F_2F_4F_5r_\mathrm{c}^3`$ (16)
$`\times G^2\mathrm{ln}\mathrm{\Lambda }\left({\displaystyle \frac{m}{7\beta +2}}{\displaystyle \frac{m_0}{9}}\right){\displaystyle \frac{3M_0}{4\pi r_\mathrm{c}^3}}{\displaystyle \frac{3M}{4\pi r_\mathrm{c}^3}}\left({\displaystyle \frac{GM_0}{r_\mathrm{c}}}\right)^{1/2}.`$
The factor $`F_5`$ is a function of $`\beta `$ defined as

$$_0^{\mathrm{}}\frac{x^2dx}{\left(1+x^2\right)^{(7\beta +12)/4}}=\frac{\sqrt{\pi }}{15}\frac{\mathrm{\Gamma }(13/4)}{\mathrm{\Gamma }(15/4)}F_5(\beta ),$$
(17)
where $`F_5=1`$ at $`\beta =1`$.

Recall that gravitational encounters in a star cluster are stochastic. The energy change rates (12) and (14) have been obtained as ensemble averages. They are not locally defined quantities. The energy change rate for the entire minor component (16) has been obtained by assuming $`dE/dt=dE/dt`$ because the cluster contains many stars.

The first and second integrals in equation (12) represent heating and cooling by the major component, respectively. They yield the terms $`m_0/9`$ and $`m/(7\beta +2)`$, respectively, in equations (14) and (16). When the two terms balance, the total energy of the minor component settles to an equilibrium, $`dE/dt=0`$. Hence the $`\beta `$ value corresponding to the final state of mass segregation is

$$\beta _{\mathrm{fin}}=\frac{9m2m_0}{7m_0}.$$
(18)
This relation does not imply equipartition of kinetic energies, $`m\sigma ^2=m_0\sigma _0^2`$ (see eqs.  and ), because the major component is kept to follow Plummer’s nonthermal model. The condition for energy equipartition in an equilibrium state is that each of the two components has a thermal velocity distribution (Chapman & Cowling 1970, p. 80). A qualitatively consistent result for very large values of $`m/m_0`$ was obtained by Chatterjee, Hernquist, & Loeb (2002).111 Chatterjee et al. obtained the squared velocity dispersion $`\sigma ^2`$ at $`r0`$ for a Brownian motion of a single massive object that is in equilibrium with the background stars. Their study is mathematically equivalent to ours because we have ignored encounters among objects in the minor component. The background stars were assumed to follow Plummer’s model as in our present study.

## 4 TIMESCALE OF MASS SEGREGATION

By using the energy change rate $`dE/dt`$ in equation (16), we estimate the timescale of mass segregation $`\tau _{\mathrm{ms}}`$. This timescale is defined in terms of the secular increase of the mass of the minor component, i.e., the more massive stars, in the cluster core $`M_\mathrm{c}`$:

$$\tau _{\mathrm{ms}}=\frac{M_\mathrm{c}}{dM_\mathrm{c}/dt}|_{t=0},$$
(19)
with

$$M_\mathrm{c}=4\pi _0^{r_\mathrm{c}}r^2\rho 𝑑r=\frac{F_1F_6}{2\sqrt{2}}M,$$
(20)
for $`f=f_\beta `$. The factor $`F_6`$ is a function of $`\beta `$ defined as

$$_0^1\frac{x^2dx}{\left(1+x^2\right)^{(7\beta +3)/4}}=\frac{1}{6\sqrt{2}}F_6(\beta ),$$
(21)
where $`F_6=1`$ at $`\beta =1`$.

Since it is not analytically determined how the minor-component PDF changes in the course of mass segregation, we study two extreme cases. The one is that the minor-component PDF is characterized by a single $`\beta `$ value, which increases from the initial value of 1 toward the final value $`\beta _{\mathrm{fin}}`$, i.e., $`f=f_\beta `$ with $`\beta =1`$ at $`t=0`$ ($`1\beta \beta _{\mathrm{fin}}`$). Equations (6), (20), and (21) yield $`dM_\mathrm{c}/d\beta `$, while equations (6), (10), and (11) yield $`dE/d\beta `$. Then $`dM_\mathrm{c}/dt`$ is obtained as $`(dM_\mathrm{c}/d\beta )(d\beta /dt)=(dM_\mathrm{c}/d\beta )(dE/dt)/(dE/d\beta )`$. Using their initial values, we obtain

$`\tau _{\mathrm{ms}}`$ $`=`$ $`{\displaystyle \frac{3^3\times 5\times 7^2\pi \sqrt{\pi }}{2^{23}\sqrt{2}}}{\displaystyle \frac{19!!}{16!!}}{\displaystyle \frac{\mathrm{\Gamma }(15/4)}{\mathrm{\Gamma }(13/4)}}{\displaystyle \frac{78\mathrm{ln}2}{8\mathrm{ln}(1+\sqrt{2})3\sqrt{2}\mathrm{ln}22\sqrt{2}}}`$ (22)
$`\times {\displaystyle \frac{1}{G^2\mathrm{ln}\mathrm{\Lambda }}}{\displaystyle \frac{1}{mm_0}}\left({\displaystyle \frac{3M_0}{4\pi r_\mathrm{c}^3}}\right)^1\left({\displaystyle \frac{GM_0}{2r_\mathrm{c}}}\right)^{3/2}.`$
The numerical factor is about 0.3878.

The other extreme case is that the state $`\beta =1`$ evolves directly into the state $`\beta =\beta _{\mathrm{fin}}`$. These two states serve as the bases for the minor-component PDF, i.e., $`f=(M_{\beta =1}/M)f_{\beta =1}+(M_{\beta =\beta _{\mathrm{fin}}}/M)f_{\beta =\beta _{\mathrm{fin}}}=(1\gamma )f_{\beta =1}+\gamma f_{\beta =\beta _{\mathrm{fin}}}`$ with $`\gamma =0`$ at $`t=0`$ ($`0\gamma 1`$). The core mass $`M_\mathrm{c}`$ and the total energy $`E`$ of the minor component are

$$M_\mathrm{c}=(1\gamma )\frac{1}{2\sqrt{2}}M+\gamma M,$$
(23)
and

$$E=(1\gamma )\frac{3\pi }{32}\frac{GM_0M}{r_\mathrm{c}}\gamma \frac{GM_0M}{r_\mathrm{c}}.$$
(24)
Here we have made approximations for the state $`\beta =\beta _{\mathrm{fin}}`$ taking account of the fact that the minor component is well concentrated in the core and has a negligible velocity dispersion. Then $`dM_\mathrm{c}/dt`$ is obtained as $`(dM_\mathrm{c}/d\gamma )(d\gamma /dt)=(dM_\mathrm{c}/d\gamma )(dE/dt)/(dE/d\gamma )`$. Using their initial values, we obtain

$`\tau _{\mathrm{ms}}`$ $`=`$ $`{\displaystyle \frac{3^2\times 5\times 7^2\sqrt{\pi }}{2^{23}}}{\displaystyle \frac{19!!}{16!!}}{\displaystyle \frac{\mathrm{\Gamma }(15/4)}{\mathrm{\Gamma }(13/4)}}{\displaystyle \frac{323\pi }{2\sqrt{2}1}}`$ (25)
$`\times {\displaystyle \frac{1}{G^2\mathrm{ln}\mathrm{\Lambda }}}{\displaystyle \frac{1}{mm_0}}\left({\displaystyle \frac{3M_0}{4\pi r_\mathrm{c}^3}}\right)^1\left({\displaystyle \frac{GM_0}{2r_\mathrm{c}}}\right)^{3/2}.`$
The numerical factor is about 0.6331. The present case has a somewhat larger timescale than the previous case because the effective energy difference is somewhat larger, i.e., $`dE/d\gamma <dE/d\beta <0`$ while $`dM_\mathrm{c}/d\gamma dM_\mathrm{c}/d\beta `$ in the initial state.

Mass segregation in a real star cluster is expected to be intermediate between the above two cases. According to Fokker-Planck and Monte Carlo simulations of early stages of mass segregation (Inagaki & Wiyanto 1984; Fregeau et al. 2002), the ratio of velocity dispersions $`\sigma /\sigma _0`$ averaged over the whole cluster does not change appreciably, while that averaged over the core alone approaches steadily to the final value. The reason is that gravitational encounters occur preferentially among stars on orbits that pass through the dense core. Although this fact has been ignored in our present calculation, the behavior of $`\sigma /\sigma _0`$ would be reproduced if the minor-component PDF were a superposition of various states between those of $`\beta =1`$ and $`\beta =\beta _{\mathrm{fin}}`$ with the average $`\beta `$ value increasing with the time.

Equations (22) and (25) have the same dependence on physical quantities. The values of the numerical factors are also similar to each other and well approximated to be $`0.5\pm 0.1`$ or $`2^1`$. Thus the timescale of mass segregation is

$$\tau _{\mathrm{ms}}=\frac{1}{2}\frac{1}{G^2\mathrm{ln}\mathrm{\Lambda }}\frac{1}{mm_0}\left(\frac{3M_0}{4\pi r_\mathrm{c}^3}\right)^1\left(\frac{GM_0}{2r_\mathrm{c}}\right)^{3/2},$$
(26)
where $`3M_0/4\pi r_\mathrm{c}^3`$ and $`(GM_0/2r_\mathrm{c})^{1/2}`$ correspond to the central mass density $`\rho _0(0)`$ and central velocity dispersion $`\sigma _0(0)`$ of the major component, respectively (eqs.  and ).

## 5 DISCUSSION

### 5.1 Comparison with Relevant Timescales

The timescale of mass segregation is greater than the core crossing timescale. Since the core radius is $`r_\mathrm{c}`$ and the central velocity dispersion is $`(GM_0/2r_\mathrm{c})^{1/2}`$, the core crossing timescale $`\tau _{\mathrm{cc}}`$ is

$$\tau _{\mathrm{cc}}=\left(\frac{2r_\mathrm{c}^3}{GM_0}\right)^{1/2}.$$
(27)
Then equation (26) for the timescale of mass segregation $`\tau _{\mathrm{ms}}`$ is rewritten as

$$\tau _{\mathrm{ms}}=\frac{\pi }{6}\frac{m_0N_0/(mm_0)}{\mathrm{ln}\left[m_0N_0/(m+m_0)\right]}\tau _{\mathrm{cc}}\tau _{\mathrm{cc}},$$
(28)
where we have assumed $`m_0/m10^1`$$`10^2`$ and $`N_010^5`$$`10^6`$ as a practical situation. On the other hand, the timescale of mass segregation is less than the relaxation timescale $`\tau _{\mathrm{rlx}}`$ of the major component, which is about $`N_0/8\mathrm{ln}N_0`$ times the crossing timescale $`\tau _{\mathrm{cc}}`$ (Binney & Tremaine 1987, p. 190):

$$\tau _{\mathrm{rlx}}\frac{1}{8}\frac{N_0}{\mathrm{ln}N_0}\tau _{\mathrm{cc}}>\tau _{\mathrm{ms}}.$$
(29)
These results justifies our assumption that the star cluster is almost collisionless and the major-component PDF remains the same during mass segregation.

The timescale of mass segregation is close to the timescale of energy loss from the minor component, $`\tau _{\mathrm{el}}=E/(dE/dt)`$ at $`t=0`$. Using equations (10) and (16), we obtain

$`\tau _{\mathrm{el}}`$ $`=`$ $`{\displaystyle \frac{3^3\times 5\times 7^2\pi \sqrt{\pi }}{2^{23}}}{\displaystyle \frac{19!!}{16!!}}{\displaystyle \frac{\mathrm{\Gamma }(15/4)}{\mathrm{\Gamma }(13/4)}}{\displaystyle \frac{1}{G^2\mathrm{ln}\mathrm{\Lambda }}}{\displaystyle \frac{1}{mm_0}}\left({\displaystyle \frac{3M_0}{4\pi r_\mathrm{c}^3}}\right)^1\left({\displaystyle \frac{GM_0}{2r_\mathrm{c}}}\right)^{3/2}.`$ (30)
The dependence on physical quantities is the same as that of equation (26). This is because the energy loss from the minor component and the increase of its central concentration are of the same phenomenon. The numerical factor of equation (30) is about 0.4832, which happens to be in agreement with the numerical factor 2-1 of equation (26). Since this energy loss timescale is obtained without any assumption on the evolution of the minor-component PDF, it is a robust measure.

### 5.2 Comparison with Spitzer (1969)

Spitzer (1969) considered an isothermal cluster of uniform density and estimated the timescale of mass segregation as the timescale of approach toward energy equipartition between the major and minor components,222 Defined as the time it takes for the difference between $`m\sigma ^2`$ and $`m_0\sigma _0^2`$ to decrease by a factor of $`e^1`$. which was about $`m_0/m`$ times the relaxation timescale $`\tau _{\mathrm{rlx}}`$ of the major component. Using equation (29), we obtain the corresponding timescale in our model cluster:

$$\tau _{\mathrm{ee}(\mathrm{Spitzer})}\frac{1}{8}\frac{m_0N_0/m}{\mathrm{ln}N_0}\tau _{\mathrm{cc}}.$$
(31)
If $`mm_0`$, this energy equipartition timescale is by definition equivalent to our energy loss timescale $`\tau _{\mathrm{el}}`$ in equation (30). Since the numerical factor is $`2^1`$, the equation is rewritten as

$$\tau _{\mathrm{el}}=\frac{\pi }{6}\frac{m_0N_0/m}{\mathrm{ln}\left(m_0N_0/m\right)}\tau _{\mathrm{cc}},$$
(32)
where we have replaced $`m\pm m_0`$ with $`m`$. Except for an unimportant difference in the argument of the logarithm, the timescales (31) and (32) have the same dependence on physical quantities. On the other hand, the numerical factors are different. This difference is not serious because there is ambiguity in the definition of the relaxation timescale. Our relaxation timescale (29) is based on the core crossing timescale $`\tau _{\mathrm{cc}}`$ and might not be sufficiently large to yield the energy equipartition timescale for the whole cluster. Overall, although Spitzer (1969) considered a very idealized cluster, his timescale is expected to be valid for energy equipartition in a real star cluster at least as an order-of-magnitude estimation. This is also the case for mass segregation, which occurs during approach toward energy equipartition.

Spitzer (1969) also predicted that, if the total mass of the minor component is much less than that of the major component, $`MM_0`$, the two components eventually achieve energy equipartition, $`m\sigma ^2=m_0\sigma _0^2`$. His prediction has been confirmed at least for the cluster core by Fokker-Planck and Monte Carlo simulations (Inagaki & Wiyanto 1984; Watters, Joshi, & Rasio 2000: Fregeau et al. 2002). Since it is required that both the major and minor components evolve to have thermal velocity distributions (§3), the energy equipartition is achieved at $`t\tau _{\mathrm{rlx}}`$.333 This energy equipartition in the cluster core is not exact. The difference between $`m\sigma ^2`$ and $`m_0\sigma _0^2`$ is very small but still existent ($`1`$%; Inagaki & Wiyanto 1984; Watters et al. 2000). The reason is that the velocity distributions do not become exactly thermal even in the cluster core. This is beyond the limit of our model ($`t<\tau _{\mathrm{rlx}})`$, where the cluster is assumed to be almost collisionless. The major component is kept to follow Plummer’s nonthermal model and hence does not achieve energy equipartition with the minor component, regardless of the total mass ratio $`M/M_0`$. This fact is not serious to our estimation of the timescale of mass segregation, which is less than the relaxation timescale $`\tau _{\mathrm{rlx}}`$ and has been defined at the initial time of the evolution $`t=0`$ in equation (19).

### 5.3 Application to Observations

If the central mass density $`\rho _0(0)`$ and central velocity dispersion $`\sigma _0(0)`$ have values typical of globular clusters (Binney & Tremaine 1987, p. 26, see also p. 423 for the cutoff factor $`\mathrm{\Lambda }`$), equation (26) is written as

$$\tau _{\mathrm{ms}}=3\times 10^7\mathrm{yr}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{10}\right)^1\left(\frac{m}{10M_{\mathrm{}}}\right)^1\left(\frac{\rho _0(0)}{10^4M_{\mathrm{}}\mathrm{pc}^3}\right)^1\left(\frac{\sigma _0(0)}{10\mathrm{km}\mathrm{s}^1}\right)^3.$$
(33)
Here we have assumed $`mm_0`$ and replaced $`mm_0`$ with $`m`$. If we instead use the typical values of the core radius $`r_\mathrm{c}`$ and central velocity dispersion $`\sigma _0(0)`$, both of which are observable quantities, equation (26) is rewritten as

$$\tau _{\mathrm{ms}}=2\times 10^7\mathrm{yr}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{10}\right)^1\left(\frac{m}{10M_{\mathrm{}}}\right)^1\left(\frac{r_\mathrm{c}}{1\mathrm{pc}}\right)^2\left(\frac{\sigma _0(0)}{10\mathrm{km}\mathrm{s}^1}\right).$$
(34)
Thus mass segregation in a typical star cluster has the timescale $`\tau _{\mathrm{ms}}`$ of order $`10^7`$ yr.

The above timescale of mass segregation is used to discuss the origin of a very massive black hole found in the starburst galaxy M82 (Kaaret et al. 2001; Matsumoto et al. 2001; see also Matsushita et al. 2000). At the 2 µm secondary peak, i.e., an active site of star formation, there is a source of compact X-ray emission. The observed strong variability implies that the source is an accreting black hole. The observed luminosity of $`10^{41}`$ ergs s-1 implies that the mass is greater than $`10^3`$$`M_{\mathrm{}}`$ if the emission is isotropic and its luminosity is below the Eddington limit. This black hole should have been formed via successive mergers of massive stars and black holes that had segregated into the core of a star cluster (Taniguchi et al. 2000). Since the 2 µm secondary peak has a starburst age of about $`1\times 10^7`$ yr (Satyapal et al. 1997), the mass segregation should have occurred well within this short duration. The timescales given in equations (33) and (34) are too long. We suspect that the star cluster is exceptionally dense and compact. For example, the young star cluster R136 in the Large Magellanic Cloud has $`\rho _0(0)10^6M_{\mathrm{}}`$ pc-3, $`\sigma _0(0)10^1`$ km s-1, and $`r_\mathrm{c}10^2`$$`10^1`$ pc (Portegies Zwart et al. 1999), which yields $`\tau _{\mathrm{ms}}10^5`$ yr. In fact, only one black hole heavier than $`10^3M_{\mathrm{}}`$ has been found among more than 100 young star clusters of M82. The objects that had segregated into the cluster core should have merged in a runaway manner (Quinlan & Shapiro 1989; Portegies Zwart et al. 1999; Mouri & Taniguchi 2002) because the merger probability is higher for more massive objects.

This work has been supported in part by the Japanese Ministry of Education, Science, and Culture under grants 10044052 and 10304013. We are grateful to N. Kanaeda and the referee for helpful comments.

## FIGURE CAPTIONS


format: text

text:
# Deep spectroscopy of 𝐳∼𝟏 6C radio galaxies - II. Breaking the redshift–radio power degeneracy

## 1 Introduction

Powerful radio sources are typically associated with massive ellipticals, and are often observed to have extended emission line regions aligned along the radio source axis. The luminosity and kinematic properties of the emission line regions are generally observed to be more extreme for radio sources at high rather than at low redshifts. It is important to establish exactly how the properties of the extended emission line regions (ionization mechanism, kinematics and physical extent) are functions of redshift, radio power and radio size. This will add to our understanding of these complex systems.

Baum & McCarthy (2000) carried out a spectroscopic study of a sample of 52 radio galaxies covering a large range of redshifts, most of which are selected from the 3CR sample. Their analysis of the kinematic and morphological properties of these objects demonstrated a number of important correlations. The kinematic properties of the emission line gas in these sources were found to vary strongly with redshift and/or radio power, the higher redshift galaxies generally displaying greater line widths and velocity amplitudes. The inferred mass of ionized gas and the enclosed dynamical gas were also seen to increase with redshift. Although jet-cloud interactions are not excluded by the data, Baum & McCarthy preferred a gravitational origin for the observed kinematics.

Best, Röttgering & Longair (2000a,b) carried out a study of the emission line regions of high redshift 3CR galaxies at $`z1`$. Small radio sources (those with a radio size $`D_{rad}<120\mathrm{kpc}`$) were observed to exist in a lower ionization state, and to possess more distorted velocity profiles and boosted low-ionization (e.g. \[Oii\]) emission line luminosities in comparision with larger radio galaxies ($`D_{rad}>120\mathrm{kpc}`$). Their spectra were consistent with small radio sources being predominantly shock ionized, and large radio sources photoionized by the AGN. Emission line regions were generally larger in spatial extent in the smaller radio sources. The kinematic and ionization properties of these emission line regions were also found to be strongly correlated.

Best et al also compared the properties of the emission line regions of high redshift 3CR galaxies at $`z1`$ with a sample of low redshift 3CR galaxies with z $`\mathrm{}<0.2`$ (from Baum et al 1992). Similar radio size trends were observed in the low redshift sample, although some trends with either radio luminosity and/or redshift were also observed. Line luminosities and equivalent widths were generally less extreme in the low redshift sources. Tadhunter et al (1998) showed that a similar correlation exists between line luminosity and radio luminosity in a subsample of 2 Jy sources from the sample of Wall & Peacock (1985). Assuming a simple quasar illumination model, the observed ionization state of the emission line gas implied that this correlation could not be explained by a simple increase in the flux of ionizing photons associated with more powerful AGN at higher redshifts. The observations could however be reconciled with the simple photoionization model if emission line region cloud density is enhanced at higher redshifts, either due to changes in the host galaxy environment or radio source shocks increasing in importance.

The flux–limited 3CR sample suffers, however, from Malmquist bias: the strong correlation between redshift and radio power in a flux limited sample prevents the effects of changes in these two parameters from being disentangled from each other. This degeneracy needs to be broken, in order to determine to what extent changes in redshift and/or the host galaxy environment, in particular any possible variations in the size and structure of the emission line regions and the density of the IGM, influence the ionization and kinematics of the extended emission line regions of these galaxies. Higher radio power galaxies generally have intrinsically more powerful jets; the bulk kinetic power of the radio jets has been shown to be strongly correlated with narrow line region luminosities (Rawlings & Saunders 1991). It is important that we fully understand the effects of changes in radio power on the ionization state and kinematic properties of the extended emission line regions (EELRs).

Over recent years we have been involved in an ongoing programme of multiwavelength observations investigating a subsample of 11 6C radio galaxies, (Eales 1985; Best et al 1999; Inskip et al (hereafter Paper 1)). These sources were selected in order to be well matched to a previously well studied subsample of 28 3CR galaxies with redshift $`z1`$, (Longair, Best & Röttgering 1995; Best, Longair & Röttgering 1996, 1997; Best et al 1998, 2000a,b). In Paper 1 the results of deep spectroscopic observations of eight of these sources using the William Herschel Telescope (WHT) were presented. We refer the reader to Paper 1 for details of sample selection, observations and data reduction, reduced 1- and 2-dimensional spectra, tabulated line fluxes and composite spectra. Also included in Paper 1 is a discussion of the kinematic and ionization properties of the 6C galaxies, contrasting them with those of the spectroscopic sample of 3CR galaxies at $`z1`$.

The major results of Paper 1 are briefly summarized below:

1. The kinematical properties of the EELRs of 6C galaxies are similar to those of the more powerful 3CR sources studied at the same redshift. Small radio sources generally possess more extensive emission line regions with a more distorted velocity profile and more extreme kinematics than those of larger radio sources.
2. \[Oii\] emission line luminosity is anticorrelated with the size (or age) of the radio source.
3. Ionization state varies similarly with radio size for both subsamples, despite the decrease in radio power of the 6C sources. This is interpreted as being due to a changing contribution of ionizing photons from the shock front associated with the expanding radio source. The optical/UV spectra of the EELRs associated with large radio sources are dominated by photoionization by the AGN and those of small sources by shocks.

There is a great deal of evidence that shocks associated with the expanding radio source play a major role in creating the extreme gas kinematics observed in the spectra of radio galaxies. In addition to the results outlined above, the extreme line widths observed in the extended emission line regions are coincident with the radio source structures. This argues against the gravitational origin for the observed kinematics of the gas preferred by Baum & McCarthy, at least for the more extreme sources.

In this paper we compare the results of the 3CR and 6C $`1`$ subsamples with lower redshift 3CR sources, selected to have similar radio powers to those in the 6C subsample. By considering samples of galaxies covering a larger region of the $`P`$$`z`$ plane (where $`P`$ is the radio luminosity of the galaxies) we are able to break the degeneracy between redshift and radio power.

The structure of the paper is as follows. In section 2, we discuss the selection of low redshift galaxies for comparison with the $`z1`$ 3CR and 6C subsamples. The results of our comparison between low and high redshift sources are presented in section 3, and discussed in section 4. Conclusions are drawn in section 5. Values for the cosmological parameters $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`H_0=65\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ are assumed.

## 2 selection of the low redshift sample

The comparison of the 6C spectroscopic data with those of a more powerful sample of 3CR galaxies at the same redshift gives us an insight into the variation of the ionization states and kinematical properties of radio galaxies with radio power. By selecting another subsample of galaxies at lower redshift, matched in radio power to the 6C subsample at $`z1`$, we can investigate the dependence of these galaxy properties with redshift alone, as the combination of three such samples breaks the $`P`$$`z`$ degeneracy. Unfortunately there are no data available in the literature exactly matching our criteria of full ionization and kinematic data for a large number of galaxies matched in radio power to the 6C subsample. The low redshift data used in this paper has therefore been taken from two different samples: ionization data from the sample of Tadhunter et al (1998), and the low redshift kinematic data from the sample of Baum & McCarthy (2000).

For ionization studies we have selected narrow line radio galaxies from the steep spectrum selected subsample of Tadhunter et al (1993), which is itself taken from the 2 Jy sample of Wall & Peacock (1985). The Tadhunter sample includes all sources from the Wall & Peacock sample with 2.7 GHz flux density $`>2`$Jy, declinations $`\delta <10^{}`$ and redshifts $`z<0.7`$. Whilst the full 2-Jy sample is not spectroscopically complete, the remaining objects without identified redshifts are very faint, and unlikely to lie within the $`z<0.7`$ subsample of Tadhunter et al (1993). This upper redshift limit is also the lower redshift limit of the spectroscopic subsample of 3CR galaxies of Best et al (2000a), which spanned the redshift range $`0.7<z<1.25`$. The eighteen sources selected have both \[Oiii\]/\[Oii\] and \[Oiii\]/H$`\beta `$ line ratios available, as presented by Tadhunter et al (1998). Eight of these sources lie in a similar radio luminosity range to the 6C $`z1`$ subsample, and thus can be used in a statistical comparison with the high redshift sources. Two sources were excluded on the basis of incomplete emission line information, but neither of these would have fallen in the luminosity range selected for comparison to the 6C galaxies. The sources included in the statistical comparison subsample are 0023-26, 0039-44, 0859-25, 1306-09, 1934-63, 2250-41, 2314+03 (3C459) and 2356-61. The full subsample of eighteen galaxies is 90% complete, and unbiased. Although this sample provides ionization data suitable for comparison with the $`z1`$ subsamples, we need to look elsewhere for the kinematics of low redshift radio sources.

Baum & McCarthy (2000) provide an analysis of the kinematical and morphological properties of a sample of 52 radio galaxies. Their sample is a compilation of two other samples: low redshift galaxies are taken from the samples of Baum, Heckman & van Breugel (1990, 1992) and the intermediate–high redshift galaxies from the sample of McCarthy et al (1996). Both of these were drawn in turn from the emission line imaging survey of Baum et al (1988) and imaging by McCarthy et al (1995). The Baum & McCarthy sample is slightly biased against sources with small emission line regions, particularly at higher redshifts; at redshifts greater than $`0.2`$, it includes most 3CR sources with emission line imaging sizes greater than 5<sup>′′</sup>. For comparison with the kinematical properties of $`z1`$ radio galaxies, we have selected twenty-two FRII galaxies from the sample of Baum & McCarthy (2000). These are all FRII sources with redshifts of $`z<0.7`$, the same upper redshift limit as used in the Tadhunter et al subsample. Of the twenty-two galaxies selected, twenty-one have a measured FWHM, but only eight are used in statistical comparisons with the 6C galaxies at $`z1`$, as these have similar radio luminosities. These are the sources 3C79, 3C169.1, 3C299, 3C300, 3C306.1, 3C337, 3C435 and 3C458.

Fig. 1 shows the radio–power vs. redshift distribution of the four samples considered in this paper, illustrating which of the low redshift galaxies are suitable for comparison with the $`z1`$ 6C and 3CR subsamples.

## 3 Comparison with galaxies at low redshift

### 3.1 The ionization state of the emission line gas

#### 3.1.1 \[Oii\]/H$`\beta `$ as a diagnostic test

For the Tadhunter et al subsample, the line ratios \[Oiii\]/\[Oii\] and \[Oiii\]/H$`\beta `$ are available. These can be used to determine \[Oii\]/H$`\beta `$, a line ratio which can be inferred for the $`z1`$ data based on our measurements of \[Oii\] and H$`\gamma `$ or other Balmer lines, assuming the Balmer line ratios for case B recombination. The composite spectra produced for the $`z1`$ samples (Paper 1) suggest that this line ratio is quite different for large and small radio sources, and thus may be a useful diagnostic of the ionization mechanism. To investigate this, we used the mappingsii (Sutherland, Bicknell & Dopita 1993) results for \[Oii\]/H$`\beta `$ and Ciii\]1909/Cii\]2326 to create a new ionization mechanism diagnostic diagram (Fig. 2).

Marked on the diagram are the predictions of shock (Dopita & Sutherland 1996) and photoionization models (Binette, Wilson & Storchi-Bergmann 1996; Allen, Dopita & Tsvetanov 1998; digitised data kindly provided by Mark Allen). The 6C and 3CR data points are also included. Where H$`\gamma `$, H$`\delta `$ or H$`\zeta `$ could not be measured in our spectra, an upper limit was determined for the H$`\beta `$ flux for these sources. In these cases, the \[Oii\]/H$`\beta `$ line ratio is given as a lower limit. No correction has been made for internal extinction. Assuming similar levels of internal extinction to that measured for Cygnus A (Osterbrock, 1989), this is likely to have only a minimal impact on our results.

The diagnostic plot in Paper 1 utilised the predictions of the line ratios Ciii\]1909/Cii\]2326 and \[Neiii\]3869/\[Nev\]3426 for these models. Large sources were generally well matched to pure photoionization by the AGN, with an ionization parameter $`2\mathrm{}<\mathrm{log}_{10}\mathrm{U}\mathrm{}<1`$ and spectral index $`\alpha =1`$. Small sources fitted well the predictions of the shock models incorporating precursor ionization regions. Alternatively, the photoionization models combining matter bounded and ionization bounded clouds of Binette et al also seemed to provide a good fit to the data, with the smaller radio sources typically having a smaller fraction of matter bounded clouds. The positions of the data points in relation to the different models can be compared between Fig. 2 and the diagnostic plot from Paper 1, allowing the consistency of the models in fitting the observed data to be gauged.

When comparing two or more line ratio diagnostic diagrams, it quickly becomes apparent that not all the models explain the spectra of the observed sources consistently from one plot to another. However, some of the modelled tracks appear to fit the data reasonably well on both plots. The five ‘photoionized’ sources on the diagnostic plot in Paper 1 appear close to the predicted line ratios of the same model in Fig. 2. Most of the points which were previously well matched by the predictions of the shock plus precursor region models remain so in this diagram as well. As a general rule, small sources show lower ionization parameters than larger sources, and their positions on the diagnostic plot are consistent with the predictions of shock ionization. However, the shock ionization and AGN photoionization models cannot be completely separated in Fig. 2. These results clearly illustrate that the power of line ratio plots as a diagnostic of the ionization mechanism is greatly enhanced when several different plots are considered together, as ambiguities on one or more plots may be resolved by additional modelled line ratios.

One clear result from this plot is that the preferred region of the photionization tracks predicts lower values of \[Oii\]/H$`\beta `$ than the predictions of the shock models. Thus, the \[Oii\]/H$`\beta `$ line ratio can provide a suitable indication of any changes of ionization state which may be observed in the low redshift sources, just as the Ciii\]1909/Cii\]2326 line ratio was used to do so in Paper 1.

#### 3.1.2 Balmer lines vs. \[OII\]3727 - investigating the variation in ionization state within the three samples.

In Fig. 3a the \[Oii\]/H$`\beta `$ line ratio data for the Tadhunter et al subsample and the two $`z1`$ samples are plotted against radio size. Fig. 3b displays the variation of this line ratio with radio power for the same samples; Fig. 3c plots this parameter against redshift. It is interesting to note that there is a much larger spread in the \[Oii\]/H$`\beta `$ line ratio for the lower radio power sources.

Spearman rank correlation tests have been carried out for these data and the results tabulated in Table 1. For three 6C sources, the \[Oii\]/H$`\beta `$ line ratio could only be calculated as a lower limit; these values have been treated as measured values for all statistical analyses unless stated otherwise. However, as the \[Oii\]/H$`\beta `$ line ratios for these sources are amongst the highest in the samples, the statistics are unlikely to suffer. The line ratios of the galaxies in all three subsamples combined are anticorrelated with radio size at the 99.8% significance level, suggesting that this line ratio is a suitable alternative to Ciii\]1909/Cii\]2326 as an indicator of changes in ionization state with radio size. The data are slightly less strongly correlated with radio power and redshift.

We have also carried out partial rank correlations for these galaxies, using the method described by Macklin (1982), which enables the true dependencies of the different properties to be determined, by removing the effects of correlations with one or more of any other parameters.

Our results show strong anticorrelations of the \[Oii\]/H$`\beta `$ line ratio with both radio size and luminosity (see Table 1). Although a Spearman rank test initially found a significant correlation with redshift, partial rank analysis shows that redshift is an unimportant parameter and the correlation was induced by the strong correlation between $`z`$ and radio power. We can conclude that there is no strong dependence of ionization state on any changes with redshift in the properties of the host galaxy or its environment.

Another way of determining the most important correlations within a set of data is a principal component analysis (also known as a Karhunen-Loeve transform). This is described by Kendall (1980) and Efstathiou & Fall (1984). For a data set with p variables for each of N parameters, a new set of uncorrelated parameters can be determined (the principal components of the data) which are linear functions of the original variables. If the original variables are standardized to zero mean and unit variance, the principal components will also be independent of each other. The first principal component of the data has the largest possible variance of any linear function, and can be thought of as defining the line of closest fit to the p–dimensional set of N data points, in that the sum of squares of the distances of the data points in any direction from such a line will be minimised. The magnitude of this eigenvector indicates what fraction of the dataset can be explained by these correlations alone. The second principal component has the largest variance subject to being uncorrelated (i.e. orthogonal) to the first. Further principal components are uncorrelated with all the preceding principal components.

The results of a principal component analysis of our data are displayed in Table 2, and confirm the results of the partial rank coefficients. This analysis has been carried out twice for the three galaxy subsamples; first using only the low-redshift data matched in radio power to the 6C subsample, and secondly using the full low redshift subsample. Although there are subtle differences between the two sets of results, the overall correlations determined are the same in both cases.

In both cases, the eigenvector of the first principal component accounts for $`45`$% of the variance within the sample, that of the second component $`30\%`$, the third component $`20\%`$ and the fourth $`<10\%`$. A positive correlation between $`P`$ and $`z`$ is clearly the strongest correlation within the data, being inverted significantly only in the final eigenvector, in order to remove the remaining scatter in the data. The advantage of using the chosen combination of three samples is that the effects of this correlation can be accounted for, and the redshift–radio power degeneracy successfully broken. Anti-correlations of the \[Oii\]/H$`\beta `$ line ratio with both redshift and radio power are implied by the first principal component. However the second component reverses and weakens these correlations, particularly for redshift. In the first eigenvector, ionization state and radio size ($`D_{rad}`$) are seen to be anticorrelated for the matched samples, and very weakly correlated for the full low redshift sample. The second eigenvector shows a strong anticorrelation between the ionization state and size of the radio source in both cases, suggesting that this is indeed a significant correlation. This difference is perhaps due to an underlying anticorrelation between radio size and either $`P`$ or $`z`$, which weakens the relationship between radio size and the \[Oii\]/H$`\beta `$ line ratio, although it is unclear from these data whether anticorrelations between either $`P`$ and $`D_{rad}`$ or between $`z`$ and $`D_{rad}`$ are real. An anticorrelation between $`D_{rad}`$ and $`z`$ is observed by Neeser et al (1995) and Blundell et al (1999), and is likely to be due in part to the selection effects inherent within flux limited samples, which are biased towards smaller sources at high redshift due to Malmquist bias and the decrease in flux density with radio size for individual sources (Kaiser, Dennett–Thorpe & Alexander 1997).

These results are fully consistent with the results of the partial rank correlation tests, and suggest that the following are real: a strong positive correlation between $`P`$ and $`z`$, and weaker anticorrelations between $`D_{rad}`$ and $`z`$ or $`P`$, due to the selection effects in defining the sample. Anticorrelations between $`I`$ and $`P`$, and between $`I`$ and $`D_{rad}`$ are also real, where ‘$`I`$’ represents the \[Oii\]/H$`\beta `$ line ratio as a measure of the ionization state within the emission line gas, a large value for $`I`$ implying a low ionization state. Of the implied correlations with ionization state, the \[Oii\]/H$`\beta `$ line ratio is most strongly anticorrelated with radio size, followed by radio power.

### 3.2 The Baum & McCarthy sample: Kinematics comparison

The first parameter we consider in our comparison of the kinematics of radio galaxies at high and low redshift is the maximum value of the velocity FWHM, which varies with spatial position. This parameter is very robust, easily determined and provides a straightforward measure of the kinematic effects of any shocks associated with jet-cloud interactions. The emission line gas FWHM is plotted against radio size, radio power and redshift in Fig. 4.

In general, the lower redshift sources have smaller FWHM than the two high redshift samples. The FWHM of both 6C and 3CR galaxies are anticorrelated with radio size, whilst the low redshift galaxies selected from the Baum & McCarthy sample are not correlated very significantly with this parameter. Baum & McCarthy claim that there is little evidence for an anticorrelation between radio size and the velocity FWHM of the emission line gas in their complete sample, such as that seen in both $`z1`$ subsamples. It is not surprising that this result is obtained, as by considering galaxies over their full range of redshifts, any positive correlation between redshift and FWHM such as that observed for our three samples will act to diminish any anticorrelation between radio size and FWHM. Best et al (2001) separated the full Baum & McCarthy sample into subsamples of sources with $`z<0.6`$ and $`z>0.6`$; the high and low redshift subsamples each displayed a strong anticorrelation between velocity FWHM and radio size. In the redshift range we consider, the sample is nearly complete for sources with large emission line regions and an anticorrelation between FWHM and redshift is observed. Considering all three samples together, FWHM is clearly anticorrelated with radio size, and correlated with redshift and radio power. The strengths of these correlations are listed in Table 3.

As before, some of these correlations may not indicate an intrinsic connection between FWHM and redshift, radio power or radio size, but instead be due to secondary correlations between these other parameters. This can be disentangled by carrying out partial rank correlation tests on the data, removing the influence of one or more of the other parameters. The partial rank analysis shows that the strongest relationship between the parameters is an anticorrelation with radio size ($`>99.9`$% significant), followed by a positive correlation with redshift ($`>98`$% significant) and a slightly weaker correlation with radio power ($`>94`$% significant). The relative importance of the trends with redshift and radio power can also be observed in Fig. 4. Fig. 4b shows that the kinematic properties of the 6C subsample are clearly more extreme than those of the low redshift 3CR subsample. These samples are matched in radio power, but not in redshift. The differences between the $`z1`$ 6C and 3CR subsamples, matched in redshift but not radio power, are much less obvious, suggesting that the gas kinematics are indeed more strongly dependent on redshift than on radio power.

We have also carried out a principle component analysis for these data sets, the results of which are tabulated in Table 4. The first principal component accounts for over 60% of the variance in the kinematic data. Practically all the scatter in the data can be accounted for with the first two eigenvectors. The reduction in the scatter of the data points compared with the ionization state data provides the first indication that any correlation between $`P`$, $`D`$ or $`z`$ with the emission region gas FWHM is likely to be much more significant than the correlations observed between these parameters and the ionization state of the gas. The first eigenvector for both data sets gives strong positive correlations between $`P`$, $`z`$ and FWHM, and a strong anticorrelation between each of these parameters and $`D_{rad}`$. In the second eigenvector the correlation between $`P`$ and $`z`$, and the anticorrelation between FWHM and $`D_{rad}`$ remain strong, but the other correlations implied by the first eigenvector are reversed and hence greatly weakened (particularly the FWHM/radio power correlation). These results fully confirm those of the partial rank analysis: the most important kinematic correlations are FWHM anticorrelated with radio size, and positively correlated with redshift and radio power, in decreasing order of significance. As for the ionization data we also observe a correlation of redshift with radio power, as well as an anticorrelation of radio size with redshift.

The variation of velocity range (defined as the difference between the most positive and negative velocity components of the \[Oii\] line) with radio size is plotted in Fig. 5 for all three samples. The distribution of all three samples appears very similar, suggesting that velocity range has no dependence on either redshift or radio power. A visual examination of Fig. 5 suggests that the greatest range of velocities is observed in the smaller radio galaxies. However, a Kolmogorov-Smirnov test, which is more suitable than a Spearman rank test for an analysis of the significance of the suspected trend, shows that the decrease in velocity range with radio size is not significant for any or all of the subsamples combined; nor is the decrease in maximum velocity range from the 3C to 6C sources statistically significant. It can be safely assumed that the velocity range data for any sample is drawn from the same population, regardless of redshift, radio size or radio power.

Another parameter we have studied previously is the luminosity of the \[Oii\] line scaled by the ratio of the radio power of a source to the mean value for the sample, and its variation with radio power (see Paper 1). Due to different observational techniques (long-slit spectroscopy versus narrow band imaging), a strict comparison between the \[Oii\] luminosities of the low and high redshift galaxies is not feasible. Nonetheless there is little variation in scaled $`L_{[\mathrm{OII}]}`$ with radio size for the low redshift sources. A decrease in scaled $`L_{[\mathrm{OII}]}`$ with radio size as seen for the $`z1`$ data is *not* observed, suggesting that any boosting of line luminosity by shocks in the smaller sources is much less important at lower redshifts. This would be consistent with the less extreme kinematics observed in the low redshift subsample. If the kinematic results are due in part to radio power/shock strength effects, as found from the strong independent correlation of FWHM with radio power, a weaker or non-existent decrease in $`L_{[\mathrm{OII}]}`$ with radio size would also be expected.

## 4 Discussion

In this paper, we have compared the ionization properties of both the 6C and 3CR subsamples at $`z1`$ with the ionization properties of the low–intermediate redshift galaxies ($`z<0.7`$) in the sample of Tadhunter et al (1998). We have shown that, whilst not as effective as the Ciii\]1909/Cii\]2326 line ratio, \[Oii\]/H$`\beta `$ can be used as a useful diagnostic of the emission line ionization mechanism. The \[Oii\]/H$`\beta `$ ratios of the Tadhunter sample galaxies at low redshift clearly display the same variation between the degree of ionization and the radio source size, or alternatively, source age and ionization mechanism. We have contrasted the kinematics of our $`z1`$ samples with galaxies with $`z<0.7`$ selected from the sample of Baum & McCarthy (2000). The data are strongly correlated with both redshift and radio luminosity, as well as being strongly anticorrelated with radio size. Even considered individually, the three subsamples display similar trends, although the kinematics of the low redshift galaxies are clearly less extreme than their high redshift counterparts. In the low redshift objects, the lack of any anticorrelation between radio size and \[Oii\] luminosity, scaled by the radio power of the host galaxy, suggests a much weaker (or perhaps non-existent) relation between emission line strength and radio source age. It also implies that shock induced boosting of emission line luminosities is much less important in the low redshift sources than for similar radio power sources (i.e. 6C galaxies) at $`z1`$.

The variations with radio size are apparent across the full range of redshifts, and have been well documented in previous studies (e.g. Best et al (2000b) and Paper 1). In general, ionization state, EELR size, velocity profile and FWHM all vary with radio size. Kinematic and ionization properties of the EELRs are found to be intrinsically correlated, despite the selection biases. It is clear that photoionization by the AGN is the only important mechanism for the more evolved, larger radio sources, which also display relatively undisturbed kinematics. The influence of shocks on the kinematics of the EELR gas clouds is certainly of importance, providing an explanation of the very distorted kinematics seen in the smaller radio sources. The clear importance of radio source shocks in these sources suggests that a gravitational origin for the observed extreme kinematics is unlikely. However, in the case of large, photoionized sources, where radio source shocks are less effective, the increased FWHM observed at high redshifts could in fact be due to a larger gravitational potentials, for example if the radio galaxies are associated with clusters (e.g. Best 2000) rather than just the potential of an individual galaxy.

We will now consider the dependence of the ionization and kinematic properties of the EELRs on redshift and radio power independently. The inclusion of 6C data in the combined sample weakens the inherent $`P`$$`z`$ correlation rather than removing it entirely. However, using the combined sample to investigate the relationships between the EELR properties and these parameters does allow the redshift–radio power degeneracy to be successfully broken. With this achieved, we can create a fuller picture of the processes occurring within the emission line regions of all three subsamples of radio galaxies. The ionization state of the gas shows no significant variation with redshift, and depends only on radio power, in addition to the size of the radio source. The interplay between the strength of any shocks involved and the ionizing UV flux of the AGN is still unclear; nevertheless, the emission line regions of the powerful $`z1`$ 3CR radio sources generally exist in a higher ionization state than the sources in the low redshift or 6C subsamples, as shown by the variation of the observed \[Oii\]/H$`\beta `$, Ciii\]1909/Cii\]2326 and \[Neiii\]3869/\[Nev\]3426 line ratios. The ionization state is clearly dependent on AGN properties rather than other host galaxy properties or the environment of the radio source.

The ionization parameter $`U`$, which is proportional to the ratio of ionizing photons impinging on the surface of a cloud to the gas density, is approximately twice as great for the large 3CR radio sources as for the large 6C sources, although the distributions of sources from both subsamples cover a range of values for $`U`$ (Paper 1). Although the difference in ionization parameter is less than the difference in radio power between the two samples (a factor of $`5`$ at $`z1`$), this anomaly could perhaps be explained by less dense clouds populating the EELRs of 6C galaxies compared with their more powerful 3CR counterparts. The sensitivity of the observations could also be an important factor. For the more powerful 3CR sources at $`z1`$, photoionized gas would be observable out to a greater distance from the central AGN than for the less powerful 6C sources at the same redshift. As the flux of ionizing photons impinging on a cloud per unit surface area decreases with distance from the AGN, the ionization parameter for more distant clouds will be lower than that for clouds closer to the AGN. The increase in the observed extent of the EELRs of more powerful radio sources will therefore lead to the average value of the ionization parameter being lower than that for a less powerful radio source observed to the same sensitivity. However, if the clouds are in pressure equilibrium, their density is also likely to decrease with distance from the AGN, reducing the variation in $`U`$ with distance. Alternatively, the ionizing photon flux may scale non-linearly with radio luminosity, as suggested by the work of Serjeant et al (1998) and Willott et al (1998). They find that the optical luminosity of quasars scales roughly as $`L_{\mathrm{rad}}^{0.5}`$. This result suggests that the difference in ionization parameter for the two $`z1`$ subsamples can be expected to be less than the factor of five difference in their radio powers.

We also observe the effects of radio power on the kinematics of the emission line regions, which are independently correlated with radio power once the correlations with other parameters (redshift and radio size) are removed. If the disturbed kinematics of smaller radio sources are due in the main to the influence of shocks, parameters such as the gas velocity FWHM are expected to be lower for the weaker shocks associated with less powerful radio galaxies. In addition, the more powerful radio sources may involve an increased number of jet–cloud interactions.

Finally, the observed decrease in gas velocity FWHM with redshift needs to be explained. For the kinematics, redshift is a much more important parameter than radio power, as shown by the results of the partial rank correlations and principal component analysis. This result seems to point to variations in either the host galaxy properties, e.g. mass of gas available in the emission line regions, or changes in, say, the IGM density with redshift. Changes in cloud sizes and/or densities would clearly affect any acceleration of the clouds due to shocks associated with the expanding radio source. An increase in jet–cloud interactions may be due to a higher gas density in the ICM for the higher redshift galaxies. High redshift 3CR galaxies are generally found in clusters, whereas low redshift 3CR sources are more likely to be located in smaller groups (e.g. Hill & Lilly 1991, Best 2000). A change in radio source environment with cosmic epoch may have an effect on the observed kinematics. The more extreme kinematics of the EELRs of higher redshift radio galaxies can plausibly be attributed in part to their richer environments as compared with low redshift radio galaxies.

## 5 Conclusions

Our conclusions are as follows:

1. A comparison between the emission line diagnostic plots in Paper 1 and the current paper suggests that the shock model with a precursor model best explains the spectra of small double radio sources. The emission line regions of larger radio sources are well explained by photoionization by the UV flux from the AGN. An accurate description of the spectra of small radio sources is likely to require a combination of both shock ionization and photoionization, even where shock ionization is likely to be the most important ionization mechanism.
2. In addition to the known variation in EELR gas kinematics with radio size, the kinematic properties are strongly correlated with redshift and radio power independently. The more extreme gas kinematics observed at earlier cosmic epochs suggest that the structure and/or composition of the EELR gas clouds may vary with redshift. Evolution in the environment of the radio source (e.g. changes in IGM/ICM density profile) may also be occuring. The dependence of the EELR kinematics on radio power is likely to be due to the influence of shocks.
3. The ionization state of the gas is strongly independently correlated with radio size, as well as less strongly correlated with radio power. The data show no dependence on the redshift of the source once radio power effects have been removed. The ionization state, and hence the dominant source of ionizing photons, is therefore only dependent on the properties of the AGN (i.e. radio size and radio power) and not on cosmic epoch.
4. Whilst there is a fairly tight distribution of scaled $`L_{[\mathrm{OII}]}`$ with radio size in the low redshift subsample, the weak positive correlation of this data suggests that unlike the strong anticorrelation of emission line luminosity with source age observed in the redshift $`z1`$ sources, this correlation is much weaker or non-existent for the lower redshift radio galaxies. This indicates that shocks are less important for the radio galaxies at low redshifts, and is consistent with our other results.

## Acknowledgements

This work was supported in part by the Formation and Evolution of Galaxies network set up by the European Commission under contract ERB FMRX–CT96–086 of its TMR programme. We would like to thank Mark Allen for providing digitised mappingsii line ratios. KJI acknowledges the support of a PPARC research studentship. PNB is grateful for the generous support offered by a Royal Society Research Fellowship. The William Herschel Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. We would like to thank the referee for useful comments.


text:
# 1. Introduction

## 1. Introduction

Dirac’s Large Numbers Hypothesis (LNH) is the origin of many theoretical explorations of time-varying $`G`$. According to the LNH, the value of $`\dot{G}/G`$ should approximately coincide with the Hubble rate. Although it has become clear in the recent decades that the Hubble rate is too high to be compatible with experiment, the enduring legacy of Dirac’s bold stroke is the acceptance by modern theories of non-zero values of $`\dot{G}/G`$ as being potentially consistent with physical reality.

There are three problems related to $`G`$, whose origin lies mainly in unified model predictions: 1) absolute $`G`$ measurements, 2) possible time variations of $`G`$, 3) possible range variations of $`G`$, i.e., non-Newtonian, or new interactions. For 1) and 3) see .

After the original Dirac hypothesis some new concepts appeared and also some generalized theories of gravitation admitting variations of the effective gravitational coupling. We can single out three stages in the development of this field:

1. Study of theories and hypotheses with variations of fundamental physical constants, their predictions and confrontation with experiments (1937-1977).
2. Creation of theories admitting variations of an effective gravitational constant in a particular system of units, analyses of experimental and observational data within these theories  (1977-present).
3. Analyses of variations of fundamental physical constants within unified models  (present).

Different theoretical schemes lead to temporal variations of the effective gravitational constant:

1. Empirical models and theories of Dirac type, where $`G`$ is simply replaced with $`G(t)`$.
2. Numerous scalar-tensor theories (STT) of Jordan-Brans-Dicke type, with $`G`$ depending on the scalar field $`\varphi (t)`$ or a number of scalar fields.
3. Gravitational theories with a nonminimally coupled (in particular, conformal) scalar field arising in different approaches  (they can actually be treated as special cases of STT).
4. Multidimensional unified theories in which there are dilatonic fields and effective scalar fields appearing in our 4-dimensional spacetime from extra dimensions . They may also help one in solving the problem of a variable cosmological constant from Planckian to present values and the cosmic coincidence problem.

A striking feature of the present status of theoretical physics is that there is no satisfactory theory unifying all four known interactions; most modern unification theories do not admit unique and universal constant values of physical constants and of the Newtonian gravitational coupling constant $`G`$ in particular. In this paper we discuss the bounds that may be suggested by a general class of STT. One can mention that STT are among the viable alternatives to general relativity; on the one hand, they are widely used for comparison with observations and, on the other, their different versions emerge in the field limits of the candidate “theories of everything”.

Although the bounds on $`\dot{G}`$ and $`G(r)`$ are in some classes of theories rather wide on purely theoretical grounds, since any theoretical model contains a number of adjustable parameters, we note that observational data concerning other phenomena, in particular, cosmological data, place limits on the possible ranges of these adjustable parameters.

Here we restrict ourselves to the problem of $`\dot{G}`$ (for $`G(r)`$ see \[1–4\]). We show that various theories predict the value of $`\dot{G}/G`$ to be $`10^{12}/`$yr or less. The significance of this fact for experimental and observational determinations of the value of or upper bound on $`\dot{G}`$ is the following: any determination with error bounds significantly better than $`10^{12}/`$yr (combined with experimental bounds on other parameters) will typically be compatible with only a small portion of existing theoretical models and will therefore cast serious doubt on the viability of all other models. In short, a tight bound on $`\dot{G}`$, in conjunction with other astrophysical observations, will be a very effective “theory killer” and/or significantly reduce the class of viable theories. Any step forward in this direction will be of utmost significance.

Some estimations of $`\dot{G}`$ had been done long ago in the framework of general scalar-tensor and multidimensional theories using the values of cosmological parameters ($`\mathrm{\Omega }`$, $`H`$, $`q`$ etc) known at that time . It is easy to show that for modern values they predict $`\dot{G}/G`$ at the level of $`10^{12}/`$yr and less (see also recent estimations of A. Miyazaki , predicting time variations of $`G`$ at the level of $`10^{13}\mathrm{yr}^1`$ for a Machian-type cosmological solution in the Brans-Dicke theory).

The most reliable, by now, experimental bounds on $`\dot{G}/G`$ (spacecraft radar ranging ) and lunar laser ranging ) give a limit of $`10^{12}/`$yr), so any results at this level or better will be very important for solving the fundamental problem of variations of constants and for discriminating between viable unified theories. So, realization of such multipurpose new generation type space experiments like Satellite Energy Exchange (SEE) for measuring $`\dot{G}`$ and also absolute value of $`G`$ and Yukawa type forces at the ranges of metres and the Earth radius  become extremely topical.

## 2. Scalar-tensor cosmology and variations of $`G`$

We are going to estimate the order of magnitude of variations of the gravitational constant $`G`$ due to cosmological expansion in the framework of scalar-tensor theories (STT) of gravity, using modern data on the cosmological parameters.

Consider the general (Bermann-Wagoner-Nordtvedt) class of STT where gravity is characterized by the metric $`g_{\mu \nu }`$ and the scalar field $`\varphi `$; the action is

$`S={\displaystyle }d^4x\sqrt{g}\{f(\varphi )[g]`$
$`+h(\varphi )g^{\mu \nu }\varphi _{,\mu }\varphi _{,\nu }2U(\varphi )+L_m\}.`$ (1)
Here $`[g]`$ is the scalar curvature, $`g=|det(g_{\mu \nu })|`$; $`f,h`$ and $`U`$ are certain functions of $`\varphi `$, varying from theory to theory, $`L_m`$ is the matter Lagrangian.

This formulation of the theory corresponds to the Jordan conformal frame, in which matter particles move along geodesics and hence the weak equivalence principle is valid, and non-gravitational fundamental constants do not change. In other words, this is the frame well describing the existing laboratory, geophysical and cosmological observations.

Among the three functions of $`\varphi `$ entering into (2.) only two are independent since there is a freedom of transformations $`\varphi =\varphi (\varphi _{\mathrm{new}})`$. We use this arbitrariness, choosing $`h(\varphi )1`$, as is done, e.g., in Ref. . Another standard parametrization is to put $`f(\varphi )=\varphi `$ and $`h(\varphi )=\omega (\varphi )/\varphi `$ (the Brans-Dicke parametrization of the general theory (2.)). In our parametrization $`h1`$, the Brans-Dicke function $`\omega (\varphi )`$ is $`\omega (\varphi )=f/f_\varphi ^2`$; here and henceforth, the subscript $`\varphi `$ denotes a derivative with respect to $`\varphi `$. The Brans-Dicke STT is the particular case $`\omega =\mathrm{const}`$, so that in (2.)

$$f(\varphi )=\varphi ^2/(4\omega ),h1.$$
(2)
The field equations that follow from (2.) read

$`\varphi \frac{1}{2}f_\varphi +U_\varphi =0,`$ (3)
$`f(\varphi )\left(_\mu ^\nu \frac{1}{2}\delta _\mu ^\nu \right)=\varphi _{,\mu }\varphi ^{,\nu }+\frac{1}{2}\delta _\mu ^\nu \varphi ^{,\alpha }\varphi _{,\alpha }`$
$`\delta _\mu ^\nu U(\varphi )+(_\mu ^\nu \delta _\mu ^\nu )fT_\mu ^\nu {}_{(\mathrm{m})}{}^{},`$ (4)
where $``$ is the D’Alembert operator, and the last term in (2.) is the energy-momentum tensor of matter.

Consider now isotropic cosmological models with the standard FRW metric

$`ds^2=dt^2a^2(t)\left[{\displaystyle \frac{dr^2}{1kr^2}}+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2)\right],`$
where $`a(t)`$ is the scale factor of the Universe, and $`k=1,0,1`$ for closed, spatially flat and hyperbolic models, respectively. Accordingly, we assume $`\varphi =\varphi (t)`$ and the energy-momentum tensor of matter in the perfect fluid form $`T_\mu ^\nu {}_{(\mathrm{m})}{}^{}=diag(\rho ,p,p,p)`$ ($`\rho `$ is the density and $`p`$ is the presuure).

The field equations in this case can be written as follows (the dot denotes $`d/dt`$:

$`\ddot{\varphi }+3{\displaystyle \frac{\dot{a}}{a}}\dot{\varphi }{\displaystyle \frac{3}{a^2}}(a\ddot{a}+\dot{a}{}_{}{}^{2}+k)+U_\varphi =0,`$ (6)
$`{\displaystyle \frac{3f}{a^2}}(\dot{a}{}_{}{}^{2}+k)=\frac{1}{2}\dot{\varphi }{}_{}{}^{2}+U3{\displaystyle \frac{\dot{a}}{a}}\dot{f}+\rho ,`$ (7)
$`{\displaystyle \frac{f}{a^2}}(2a\ddot{a}+\dot{a}{}_{}{}^{2}+k)=\frac{1}{2}\dot{\varphi }{}_{}{}^{2}+U\ddot{f}2{\displaystyle \frac{\dot{a}}{a}}\dot{f}p.`$ (8)
To connect these equations with observations, let us fix the cosmic time $`t`$ at the present epoch (i.e., consider the instantaneous values of all quantities) and introduce the standard observables:

$`H=\dot{a}/a`$ (the Hubble parameter),

$`q=a\ddot{a}/\dot{a}^2`$ (the deceleration parameter),

$`\mathrm{\Omega }_m=\rho /\rho _{\mathrm{cr}}`$ (the matter density parameter),

where $`\rho _{\mathrm{cr}}`$ is the critical density, or, in our model, the r.h.s. of Eq. (7) in case $`k=0`$: $`\rho _{\mathrm{cr}}=3fH^2`$. This is slightly different from the usual definition $`\rho _{\mathrm{cr}}=3H^2/8\pi G`$ where $`G`$ is the Newtonian gravitational constant. The point is that the locally measured Newtonian constant in STT differs from $`1/(8\pi f)`$; provided the derivatives $`U_{\varphi \varphi }`$ and $`f_{\varphi \varphi }`$ are sufficiently small, one has

$$8\pi G_{\mathrm{eff}}=\frac{1}{f}\frac{2\omega +4}{2\omega +3}.$$
(9)
(more details can be found in Refs.  where the connection between $`G_{\mathrm{eff}}`$ and $`\omega `$ was studied on the basis of cosmological solutions with local inhomogeneities and the equations of particle motion.)

Since, according to the solar-system experiments, $`\omega 2500`$, for our order-of-magnitude reasoning we can safely put $`8\pi G=1/f`$, and, in particular, our definition of $`\rho _{\mathrm{cr}}`$ now coincides with the standard one.

The time variation of $`G`$, to a good approxiamtion, is

$$\dot{G}/G\dot{f}/f=gH,$$
(10)
where, for convenience, we have introduced the coefficient $`g`$ expressing $`\dot{G}/G`$ in terms of the Hubble parameter $`H`$.

Eqs. (6)–(8) contain too many arbitrary parameters for making a good estimate of $`g`$. Let us now introduce some restrictions according to the current state of observational cosmology:

(i) $`k=0`$ (a spatially flat cosmological model, so that the total density of matter equals $`\rho _{\mathrm{cr}}`$);

(ii) $`p=0`$ (the pressure of ordinary matter is negligible compared to the energy density);

(iii) $`\rho =0.3\rho _{\mathrm{cr}}`$ (the ordinary matter, including its dark component, contributes to only 0.3 of the critical density; unusual matter, which is here represented by the scalar field, comprises the remaining 70 per cent).

Then Eqs. (7) and (8) can be rewritten in the form

$`\frac{1}{2}\dot{\varphi }{}_{}{}^{2}+U3H\dot{f}`$ $`=`$ $`2.1H^2f,`$ (11)
$`\frac{1}{2}\dot{\varphi }{}_{}{}^{2}+U2H\dot{f}\ddot{f}`$ $`=`$ $`(12q)H^2f.`$ (12)
Subtracting (12) from (11), we exclude the “cosmological constant” $`U`$, which can be quite large but whose precise value is hard to estimate. We obtain

$$\dot{\varphi }{}_{}{}^{2}H\dot{f}+\ddot{f}=(1.1+2q)H^2f.$$
(13)
The first term in Eq. (13) can be represented in the form

$$\dot{\varphi }{}_{}{}^{2}=\dot{f}{}_{}{}^{2}(df/d\varphi )_{}^{2}=\dot{f}{}_{}{}^{2}\omega /f,$$

and $`\dot{f}/f`$ can be replaced with $`gH`$. The term $`\ddot{f}`$ can be neglected for our estimation purposes. To see this, let us use as an example the Brans-Dicke theory, in which $`f=\varphi ^2/(4\omega )`$. We then have

$$\ddot{f}=(\dot{\varphi }{}_{}{}^{2}+\varphi \ddot{\varphi })/(2\omega );$$

here the first term is the same as the first term in Eq. (13), times the small parameter $`1/(2\omega )`$. Assuming that $`\varphi \ddot{\varphi }`$ is of the same order of magnitude as $`\dot{\varphi }^2`$ (or only slightly greater), we see that, generically, $`|\ddot{f}|\dot{\varphi }^2`$. Note that our consideration is not restricted to the Brans-Dicke theory and concerns the model (2.) with an arbitrary function $`f(\varphi )`$ and an arbitrary potential $`U(\varphi )`$.

Neglecting $`\ddot{f}`$, we see that (13), divided by $`H^2f`$, leads to a quadratic equation with respect to $`g`$:

$$\omega g^2+gq^{}=0,$$
(14)
where $`q^{}=1.1+2q`$.

According to modern observations, the Universe is expanding with an acceleration, so that the parameter $`q`$ is, roughly, $`0.5\pm 0.2`$, hence we can take $`|q^{}|0.4`$. (Note that this condition is only plausible rather than certain.)

In case $`q^{}=0`$ we simply obtain $`g=1/\omega `$. Assuming

$$H=h_{100}100\mathrm{km}/(\mathrm{s}.\mathrm{Mpc})h_{100}10^{10}\mathrm{yr}^1$$

and $`\omega 2500`$, we come to the estimate

$$|\dot{G}/G|410^{14}h_{100}\mathrm{yr}^1,$$
(15)
where $`h_{100}`$ is, by modern views, close to 0.7. So (15) becomes

$$|\dot{G}/G|310^{14}\mathrm{yr}^1.$$
(16)
For nonzero values of $`q^{}`$, solving the quadratic equation (14) and assuming $`q^{}\omega 1`$, we arrive at the estimate $`|g|\sqrt{q^{}/\omega }`$, so that, taking $`q^{}=0.4`$ and again $`\omega 2500`$, we have instead of (15)

$$|\dot{G}/G|1.310^{12}h_{100}\mathrm{yr}^10.910^{12}\mathrm{yr}^1,$$
(17)
where we have again put $`h_{100}=0.7`$.

We conclude that, in the framework of the general STT, modern cosmological observations, taking into account the solar-system data, restrict the possible variation of $`G`$ to values within $`10^{12}`$/yr. This estimate may be considerably tightened if the matter density parameter $`\mathrm{\Omega }_m`$ and the (negative) deceleration parameter $`q`$ will be determined more precisely.

Our estimates are rather universal since they do not use special solutions to the field equations, but actually rest on the well-justified assumption that the expansion of the Universe occurs without abrupt changes in its parameters during a fairly long period before now.

## 3. Discussion

Summarizing the above considerations, we can conclude that restrictions on possible nonzero values of $`\dot{G}`$ give no bound on the possible class of generalized gravitation theories, but in the framework of some fixed theory any restriction on $`\dot{G}`$ restricts the possible class of models.

We note that similar estimations of $`\dot{G}`$ can be made for different multicomponent multidimensional models , giving a result on the level of $`10^{12}`$/yr and less for, e.g., dust and p-brane matter sources.

We can also mention that the behaviour of the gravitational constant can actually be much more complex and intriguing than a simple time (or even range) dependence: very recently, there appeared two papers, which may open a new series of theoretical and experimental studies related to possible anisotropy in the absolute value of $`G`$  and/or its possible dependence on the latitude and longitude of the laboratory where $`G`$ was measured .