text:
# Supermassive Stars in Quasar Disks

## 1 Introduction

Accretion disks massive enough to fuel bright quasars are expected to be self-gravitating beyond a few hundred to a few thousand Schwarzschild radii ($`R_\mathrm{S}`$), (*e.g.*  Shlosman & Begelman, 1987) Even under extreme assumptions, it is difficult to see how such disks could extend beyond $`0.11\mathrm{pc}`$ without fragmenting completely into gravitationally bound objects (Goodman, 2003; Sirko & Goodman, 2003). Self-gravity is less problematic beyond $`10\mathrm{pc}`$ because the stellar bulge dominates the rotation curve, and because the disk eventually becomes optically thin.

Quasar accretion disks are not resolved and there are few direct constraints on their structure. The principal reason for believing in them at all is that one does not know of any other plausible mode of accretion that converts mass to radiation with high efficiency. The overall size of the accretion flow is however constrained by the need for sufficient surface area to radiate the observed continuum. At wavelength $`\lambda `$, this area is

$$10^{34}\left(\frac{\lambda L_\lambda }{10^{46}\mathrm{erg}\mathrm{s}^1}\right)\left(\frac{T_B}{2\times 10^4\mathrm{K}}\right)^1\left(\frac{\lambda }{1\mu \mathrm{m}}\right)^3\mathrm{cm}^2,$$

assuming incoherent emission by gas at temperature $`T_B`$. Thus the optically luminous part of the disk should extend to at least $`10^3R_\mathrm{S}`$, where self-gravity becomes important.

We shall assume that a conventional thin disk does exist and extends to $`10^3R_\mathrm{S}`$ with an accretion rate sufficient to support a typical bright QSO of mass $`10^8M_{}`$ and luminosity $`10^{46}\mathrm{erg}\mathrm{s}^1`$. We shall argue that a likely consequence of the incipient self-gravity at the outer edge of such a disk is formation of supermassive stars.

The term “supermassive” is truly justified here. While fragments may start with only $`10^2M_{}`$, they seem likely to grow to a significant fraction of the disk mass, perhaps $`10^5M_{}`$. Yet it is quite possible that even such extreme objects would escape notice until appropriate gravitational-wave detectors are built (§5).

Supermassive stars of even larger masses ($`10^8`$) were once proposed as models for quasars themselves (Hoyle & Fowler, 1963; Zel’dovich & Novikov, 1971). This notion has long been abandoned on grounds of stability (§2), and more conclusively because statistical arguments demonstrate that quasars convert mass to energy much more efficiently than nuclear fusion (Soltan, 1982; Chokshi & Turner, 1992; Yu & Tremaine, 2002). It remains possible that supermassive stars may be seeds or precursors to quasars (Rees, 1984).

In this paper, we accept the conventional view that the energy source for quasars and AGN is accretion onto massive black holes. But the disks required by this interpretation of the facts are so massive and dynamically cool that they are likely to form substructures. We expect this to occur *via* the standard local dynamical instability of self-gravitating disks, *i.e.* the gaseous version of the Toomre instability (Toomre, 1964), although growth of very large masses may perhaps be seeded by other means, for example capture of stars (§4.3).

In steady accretion, gravitational instability becomes progressively more severe towards larger radii (*e.g.* Goodman, 2003, henceforth Paper I). Thus one might chose to consider a strongly unstable region where the disk is likely to fragment completely: in other words, where most of the mass would be found in clumps rather than a smooth layer. Others have investigated this regime (Shlosman & Begelman, 1989; Shlosman et al., 1990; Kumar, 1999). However, we have chosen to focus on a transitional region where self-gravity is mild and bound structures, where they exist, are likely to be well separated. There are several reasons for our choice. For one thing, we doubt whether very strongly self-gravitating material can form an accretion flow at all, except as stars (Paper I). For another, there is little direct evidence that steady disk accretion extends beyond $`10^3R_\mathrm{S}`$, where self-gravity becomes strong. If it does, then a marginally self-gravitating region surely exists by continuity. Finally, marginal self-gravity offers the methodological advantage that one can study it as a perturbation to non-self-gravitating, thin-disk accretion, which is somewhat understood.

The plan of the paper is as follows. In §2, we review the structure of supermassive stars and what is known about their stability. §3 describes the physical conditions of a QSO accretion disk in the region where it is likely to be marginally self-gravitating, assuming an $`\alpha `$ disk in steady state (Shakura & Sunyaev, 1973). Scalings with black hole mass, accretion rate, and viscosity parameter are given. §4 builds upon these standard results to discuss the formation, growth, and fate of bound fragments (“stars”) that form in marginally self-gravitating regions, with particular emphasis on their maximum mass. §5 sums up.

## 2 Very massive stars

This section collects some elementary but important structural formulae for very massive stars; somewhat more accurate forms were given by Bond et al. (1984, henceforth BAC). We also summarize what is known about the stability of such objects, both before and after they reach the main sequence.

Very massive stars tend to be radiation-pressure dominated. To the extent that they are chemically homogeneous and the opacity ($`\kappa `$) is constant, in hydrostatic and radiative (but not necessarily nuclear) equilibrium, they are approximated by Eddington models:

$`p_{\mathrm{gas}}`$ $`=`$ $`\beta _{}p+\text{constant},`$
$`L_{}`$ $`=`$ $`(1\beta _{})L_{\mathrm{Edd},}=(1\beta _{}){\displaystyle \frac{4\pi GM_{}c}{\kappa }},`$ (1)
in which $`p_{\mathrm{gas}}`$ and $`p_{\mathrm{rad}}=pp_{\mathrm{gas}}`$ are the partial pressures of gas and radiation, and the asterisk is used to distinguish parameters of the star from those of the surrounding disk. The dimensionless parameter $`\beta _{}<1`$ is nearly constant within a given star. The additive constant on the first line above is unimportant throughout most of the interior, so that we may regard $`\beta _{}`$ approximately as $`p_{\mathrm{gas}}/p`$.

Except near the surface where the additive constant must be taken into account, and in the core where luminosity may not be proportional to mass, it follows that the stars are $`n=3`$ polytropes,

$$pK\rho ^{4/3},K=\left[\frac{3}{a}\left(\frac{k_\mathrm{B}}{\mu m_\mathrm{H}}\right)^4\frac{1\beta _{}}{\beta _{}^4}\right]^{1/3},$$
(2)
where $`\mu `$ is the molecular weight relative to hydrogen, so that $`\mu 0.62`$ for a fully ionized gas of solar metallicity. From standard results for such polytropes,

$$\frac{M_{}}{M_{}}47\frac{\sqrt{1\beta _{}}}{\beta _{}^2}\left(\frac{\mu _{}}{\mu _{}}\right)^2.$$
(3)
It is important that this relation does not depend upon the stellar radius ($`R_{}`$), so that it holds even if the star has not reached the main sequence. Thus for extreme masses $`M100M_{}`$, one has $`\beta _{}M^{1/2}`$. The binding energy

$$E_{}=\frac{3}{4}\beta _{}\frac{GM_{}^2}{R_{}}$$
(4)
is small because the gravitational and internal energy nearly cancel, and the Kelvin-Helmholtz time is correspondingly reduced:

$$t_{\mathrm{KH}}\frac{|E_{}|}{L_{}}=\frac{3\kappa M_{}}{16\pi cR_{}}\frac{\beta _{}}{1\beta _{}}.$$
(5)
The entropy per unit mass is

$$S=\frac{k_\mathrm{B}}{2\mu m_\mathrm{H}}\left[\frac{8(1\beta _{})}{\beta _{}}+\mathrm{ln}\left(\frac{1\beta _{}}{\beta _{}\mu ^2\overline{\rho }}\right)\right]+\text{constant},$$
(6)
so that the effective adiabatic index for perturbations is

$$\mathrm{\Gamma }_1=\frac{3224\beta _{}3\beta _{}^2}{3(87\beta _{})}\frac{4}{3}+\frac{1}{6}\beta _{}+O(\beta _{}^2).$$
(7)
On the main sequence, the central temperature is nearly constant. Since $`GM_{}^2/R_{}^4p_\mathrm{c}T_\mathrm{c}^4/(1\beta _{})`$, one has $`R_{}M_{}^{1/2}`$. More accurately (BAC)

$$\frac{R_{,\mathrm{ms}}}{R_{}}=10\frac{(1\beta _{})^{0.39}}{\beta _{}^{0.95}}\left(\frac{X_{CN}}{0.01}\right)^{0.05}1.6\left(\frac{M_{}}{M_{}}\right)^{0.47}\left(\frac{\mu _{}}{\mu _{}}\right)^{0.95}(M_{}10^3M_{}),$$
(8)
where $`X_{CN}`$ is the mass fraction in CNO elements. The Kelvin-Helmholtz time on the main sequence is therefore

$$t_{\mathrm{KH},\mathrm{ms}}3300\left(\frac{\kappa }{0.4\mathrm{g}\mathrm{cm}^2}\right)\left(\frac{\mu _{}}{\mu _{}}\right)^2\beta _{}^{0.053}(1\beta _{})^{0.89}\mathrm{yr}.$$
(9)
The dependence on $`M_{}`$, which enters through $`\beta _{}`$ and eq. (3), is so weak that we shall generally ignore it.

The main-sequence lifetime is also nearly constant (BAC),

$$t_{\mathrm{ms}}(23)\times 10^6\mathrm{yr}\frac{\epsilon _{\mathrm{fusion}}M_{}c^2}{L_{\mathrm{Edd},}},$$
(10)
unless cut short by catastrophic instabilities.

### 2.1 Stability of massive stars

Instabilities can be categorized according to if they occur before, during, or after the main-sequence phase.

The most famous pre-main-sequence instability is a relativistic one. Since the adiabatic index (eq. 7) is close to $`4/3`$, changes in internal and gravitational energy are almost equal and opposite if the star is homologously and adiabatically compressed. Small corrections to the energy can have a large effect on stability. One such correction is general relativity, which increases (in absolute magnitude) the potential energy by a fraction $`(GM_{}/c^2R_{})`$. A very crude estimate of the minimum unstable mass can be obtained by evaluating $`GM_{}/R_{}`$ on the main sequence and equating this to $`\beta _{}/6`$. Using equations (8) and (3), this yields $`M_{,GR}6\times 10^5M_{}`$. More detailed calculations predict a rather smaller mass (Chandrasekhar, 1964), $`M_{,GR}5\times 10^4M_{}`$, for nonrotating stars.

Rotation lends stability because its energy scales with $`p_{\mathrm{gas}}`$ rather than $`p_{\mathrm{rad}}`$ under homologous contractions, that is, $`E_{\mathrm{rot}}R_{}^2\rho _{}^{2/3}`$. The maximum rotational energy of a uniformly rotating $`n=3`$ polytrope is small, $`E_{\mathrm{rot}}/|E|0.007`$, at which value centrifugal and gravitational forces balance at the equator (Zel’dovich & Novikov, 1971). General-relativistic effects are also small, however, so rotation is important, especially since stars forming out of a disk will probably rotate strongly. We therefore take the minimum mass for general-relativistic instability to be the value appropriate to maximal uniform rotation (Baumgarte & Shapiro, 1999)

$$M_{,GR}10^6M_{}.$$
(11)
Stars with $`M_{}<M_{,GR}`$ reach the main sequence on a timescale no longer than the Kelvin time (eq. 9). Thermonuclear reactions commence on the CNO cycle, since quasar accretion disks are at least as metal rich as the Sun (Dietrich et al., 2003, and references therein). At $`M_{}100M_{}`$, the sensitivity of the thermonuclear reaction rate to central temperature drives a linear overstability to pulsations (Eddington, 1926; Ledoux, 1941; Schwarzschild & Härm, 1959). Non-linear calculations suggest that this instability drives mass loss (Appenzeller, 1970b, a; Papaloizou, 1973). The rate of mass loss is uncertain, but is probably no greater than $`\dot{m}_{}=5\times 10^5M_{}\mathrm{yr}^1`$ for a $`130M_{}`$ star and $`\dot{m}_{}=5\times 10^4M_{}\mathrm{yr}^1`$ for a $`270M_{}`$ star (Appenzeller, 1970b), and could be as little as $`10^6M_{}\mathrm{yr}^1`$ for similar masses (Papaloizou, 1973). The scaling of $`\dot{m}_{}`$ to much larger stellar masses is also uncertain.

The above mass-loss rates may be compared to those resulting from radiation pressure driven winds. Consider a $`120M_{}`$ zero age main sequence (ZAMS) star with $`T_{\mathrm{eff}}=5.3\times 10^4\mathrm{K}`$ and luminosity $`L_{}=1.8\times 10^6L_{}`$ (Schaller et al., 1992). The analytic relation of Vink et al. (2000) for radiation pressure driven winds including multiple scattering, predicts a mass-loss rate of $`1.65\times 10^5M_{}\mathrm{yr}^1`$. This mass-loss rate scales approximately as $`L_{}^{2.2}`$ and $`m_{}^{1.3}`$, with an additional weaker increase resulting from hotter surface temperatures. It can be considered as a lower limit to the actual mass-loss rates from massive stars, since these typically increase during the later periods of stellar evolution, particularly during the Wolf-Rayet phase. For example, at the above rate the stellar mass would decrease to about $`70M_{}`$ in the lifetime of $`3\mathrm{Myr}`$. However, more detailed evolutionary models including mass loss rates calculated consistently for each stage of the evolution give a final stellar mass of only $`8M_{}`$ (Schaller et al., 1992); half of the total mass loss occurs in the last 0.5 Myr of evolution. It should be noted that these results are for non-rotating stars, evolving with no additional mass accretion.

If the star manages to survive the main-sequence, then in the post-main-sequence phase an instability arises by the creation of $`e^+e^{}`$ pairs when the central temperature $`10^8\mathrm{K}`$ (Zel’dovich & Novikov, 1971). This is believed to result in an explosion if $`M_{}200M_{}`$ (BAC) but complete collapse to a black hole if $`M_{}300M_{}`$ (Fryer et al., 2001).

## 3 Initial conditions for massive star formation

Several physical properties differentiate the star-forming environment of a quasar disk from that of a typical region of Galactic massive star formation (Plume et al., 1997; McKee & Tan, 2003). Firstly the densities and pressures are many orders of magnitude greater. Secondly the environment is subject to strong shear, due to the Keplerian orbits around the massive black hole. The dominant source of heating of the gas, at least in the inner regions of the disk, is due to its accretion in the QSO’s potential. At the typical temperatures in the QSO disk we expect dust grains and molecules to have been destroyed, and finally, radiation pressure can lend significant support to the gas.

### 3.1 Physical conditions in the disk

From Paper I, we have the following properties for geometrically thin, optically thick, Keplerian disks that are heated by viscous dissipation only. The radial dependence of the surface density, $`\mathrm{\Sigma }`$, and the midplane temperature, $`T`$, are given by:

$`T`$ $`=`$ $`\left({\displaystyle \frac{\kappa \mu m_H}{16\pi ^2\alpha \beta ^{b1}k_\mathrm{B}\sigma }}\right)^{1/5}\dot{M}^{2/5}\mathrm{\Omega }^{3/5}`$
$``$ $`5.27\times 10^4\left({\displaystyle \frac{l_\mathrm{E}^2\widehat{\kappa }\mu }{\epsilon _{0.1}^2\alpha _{0.3}\beta ^{b1}}}\right)^{1/5}M_8^{1/5}r_3^{9/10}\mathrm{K},`$
$`\mathrm{\Sigma }`$ $`=`$ $`{\displaystyle \frac{2^{4/5}}{3\pi ^{3/5}}}\left({\displaystyle \frac{\mu ^4m_H^4\sigma }{k_\mathrm{B}^4}}\right)^{1/5}(\alpha \beta ^{b1})^{4/5}\kappa ^{1/5}\dot{M}^{3/5}\mathrm{\Omega }^{2/5}`$
$``$ $`2.56\times 10^5(\alpha _{0.3}\beta ^{b1})^{4/5}l_\mathrm{E}^{3/5}\epsilon _{0.1}^{3/5}\widehat{\kappa }^{1/5}\mu ^{4/5}M_8^{1/5}r_3^{3/5}\text{g cm}^2,`$
The meaning of the various symbols is as follows: $`\mathrm{\Omega }=(GM/r^3)^{1/2}`$, $`\widehat{\kappa }`$ is the opacity relative to the electron-scattering value ($`0.4\mathrm{cm}^2\mathrm{g}^1`$); $`\mu `$ is the mean particle mass in units of the Hydrogen atom mass, $`m_H`$; $`0.3\alpha _{0.3}`$ is the Shakura-Sunyaev viscosity parameter \[scaled to the largest likely value in a selfgravitating, but not fragmented, thin disk (Gammie, 2001)\]; $`\beta p_{\mathrm{gas}}/p`$; $`0.1\epsilon _{0.1}`$ is the radiative efficiency of the disk ($`=L/\dot{M}c^2`$); $`l_\mathrm{E}`$ is the disk luminosity relative to the Eddington limit; $`M=10^8M_8M_{}`$ is the black hole mass; and $`R_\mathrm{S}=2GM/c^210^5M_8\mathrm{pc}`$ is the Schwarzschild radius, with $`r_3`$ the orbital radius in units of $`10^3R_\mathrm{S}`$. Also $`b`$ takes a value of either 0 or 1, depending on whether the viscosity is proportional to the total or gas pressure, respectively.

By including the opacity as an explicit parameter, we avoid its sometimes complicated functional dependence on density and temperature, thereby obtaining simple algebraic formulae for the disk properties. The numerical disk models of Sirko & Goodman (2003, henceforth Paper II), which took realistic opacities, show that $`\kappa 1\mathrm{g}\mathrm{cm}^3`$ throughout the region of the disk of interest to us, namely $`r10^3R_\mathrm{S}`$. (Indeed, $`\widehat{\kappa }1`$ out to $`r10^4R_\mathrm{S}`$ if nonviscous auxiliary heating maintains $`Q1`$ to that radius.) Furthermore, as will be seen, many quantities of interest depend rather weakly on $`\widehat{\kappa }`$.

If viscosity scales with gas pressure ($`b=1`$) then eqs. (3.1) & (3.1) do not depend on $`\beta `$, which in any case is a known function of density and temperature:

$$\frac{\beta }{1\beta }=\frac{p_{\mathrm{gas}}}{p_{\mathrm{rad}}}=\frac{3ck_\mathrm{B}}{4\sigma m}\frac{\rho }{T^3}=\frac{3c}{8\sigma }\left(\frac{k_\mathrm{B}}{m}\right)^{1/2}\beta ^{1/2}\frac{\mathrm{\Sigma }\mathrm{\Omega }}{T^{7/2}};$$

using eqs. (3.1) & (3.1),

$`{\displaystyle \frac{\beta ^{(1/2)+(b1)/10}}{1\beta }}`$ $`=`$ $`(2^3\pi ^4)^{1/5}\alpha ^{1/10}c(k_\mathrm{B}/m)^{2/5}\sigma ^{1/10}\kappa ^{9/10}\mathrm{\Omega }^{7/10}\dot{M}^{4/5}`$ (14)
$``$ $`0.311\alpha _{0.3}^{1/10}\widehat{\kappa }^{9/10}\mu ^{2/5}(\epsilon _{0.1}/l_\mathrm{E})^{4/5}M_8^{1/10}r_3^{21/20}.`$
Using eq. (3.1) to eliminate the isothermal sound speed, $`c_s=\sqrt{p/\rho }=kT/(\mu m_H)`$ from Toomre (1964)’s gravitational stability parameter yields

$`Q`$ $`=`$ $`3(4\pi )^{3/5}\alpha ^{7/10}\beta ^{(7b12)/10}\left({\displaystyle \frac{k_\mathrm{B}}{\mu m_H}}\right)^{6/5}\sigma ^{3/10}G^1\kappa ^{3/10}\dot{M}^{2/5}\mathrm{\Omega }^{9/10}`$ (15)
$``$ $`1.0\alpha _{0.3}^{7/10}\beta ^{(7b12)/10}\epsilon _{0.1}^{2/5}l_\mathrm{E}^{2/5}\widehat{\kappa }^{3/10}\mu ^{6/5}M_8^{13/10}r_3^{27/20}.`$
Consider quasar accretion disks at the inner radius of gravitational instability where $`Q=1`$:

$$r_{\mathrm{s}.\mathrm{g}.}915\alpha _{0.3}^{14/27}\widehat{\kappa }^{2/9}\mu ^{8/9}(\epsilon _{0.1}/l_\mathrm{E})^{8/27}M_8^{26/27}\beta ^{(14b24)/27}R_\mathrm{S}.$$
(16)
This equation differs slightly from eq. (10) of Goodman (2003) because some approximations assuming $`\beta 1`$ were used there.

Substituting for $`r`$ from eq. (16) into eq. (14) yields

$$\frac{\beta _{\mathrm{s}.\mathrm{g}.}^{1b/3}}{(1\beta _{\mathrm{s}.\mathrm{g}.})^{3/4}}=0.389\alpha _{0.3}^{1/3}\widehat{\kappa }^{1/2}\mu ^1(\epsilon _{0.1}/l_\mathrm{E})^{5/6}M_8^{5/6}.$$
(17)
Thus gas and radiation pressure are comparable at $`r_{\mathrm{s}.\mathrm{g}.}`$. Equations (17) & (16) imply $`r_{\mathrm{s}.\mathrm{g}.}`$ is not very sensitive to the choice of viscosity law, *e.g.* $`r_{\mathrm{s}.\mathrm{g}}2700R_\mathrm{S}`$ for $`b=0`$ *vs.* $`1700R_\mathrm{S}`$ for $`b=1`$ at the fiducial values of the other parameters. On the other hand, radiation pressure is at least moderately dominant at all $`rr_{\mathrm{s}.\mathrm{g}.}`$, so that we may replace $`1\beta `$ with unity in the quartic (14); this makes $`\beta `$ a simple power law in all other parameters. To further simplify the discussion, we shall also take the viscosity to be proportional to gas pressure except where stated otherwise.

With these assumptions ($`\beta 1`$, $`b=1`$), the properties of the disk at $`rr_{\mathrm{s}.\mathrm{g}.}`$ are as follows. The midplane density is

$$\rho =\frac{\mathrm{\Sigma }\mathrm{\Omega }}{2c_\mathrm{s}}4.31\times 10^{10}\alpha _{0.3}^{4/5}\widehat{\kappa }^{6/5}\mu ^{4/5}(\epsilon _{0.1}/l_\mathrm{E})^{2/5}M_8^{4/5}r_3^{3/5}\mathrm{g}\mathrm{cm}^3$$
(18)
(corresponding to number densities of hydrogen nuclei $`10^{14}\mathrm{cm}^3`$), and the total pressure is

$$p/k_\mathrm{B}1.41\times 10^{20}\alpha _{0.3}^{4/5}\widehat{\kappa }^{4/5}\mu ^{4/5}(l_\mathrm{E}/\epsilon _{0.1})^{8/5}M_8^{4/5}r_3^{18/5}\mathrm{cm}^3\mathrm{K}.$$
(19)
The midplane temperature is given by eq. (3.1) with $`b=1`$ (note that at this temperature dust and molecules will be destroyed), which corresponds to an isothermal sound speed, $`c_s=\sqrt{p/\rho }`$, of

$$c_\mathrm{s}67.2\widehat{\kappa }^{1/5}l_\mathrm{E}\epsilon _{0.1}^1r_3^{9/10}\mathrm{km}\mathrm{s}^1.$$
(20)
The disk scale height is

$$h=\frac{c_\mathrm{s}}{\mathrm{\Omega }}10\widehat{\kappa }l_\mathrm{E}\epsilon _{0.1}^1R_\mathrm{S},$$
(21)
the surface density is given by eq. (3.1), and the total disk mass inside $`r`$ is

$$M_{\mathrm{disk}}(<r)5.03\times 10^5\alpha _{0.3}^{4/5}\widehat{\kappa }^{1/5}\mu ^{4/5}(l_\mathrm{E}/\epsilon _{0.1})^{3/5}M_8^{11/5}r_3^{7/5}M_{}.$$
(22)
The orbital timescale is

$$t_{\mathrm{orb}}8.79M_8r_3^{3/2}\mathrm{yr},$$
(23)
whereas the accretion time is

$$t_{\mathrm{visc}}=\frac{M_{\mathrm{disk}}(<r)}{\dot{M}}2.28\times 10^5\alpha _{0.3}^{4/5}\widehat{\kappa }^{1/5}\mu ^{4/5}(\epsilon _{0.1}/l_\mathrm{E})^{2/5}M_8^{6/5}r_3^{7/5}\mathrm{yr},$$
(24)
and the circular velocity is

$$v_{\mathrm{circ}}6710r_3^{1/2}\mathrm{km}\mathrm{s}^1.$$
(25)
Finally, we may eliminate $`\beta `$ from eq. (16) to express the radius of marginal self-gravity itself as

$$r_{\mathrm{s}.\mathrm{g}.}1550(\alpha _{0.3}/\mu )^{1/3}(l_\mathrm{E}/\epsilon _{0.1})^{1/6}\widehat{\kappa }^{1/2}M_8^{1/2}R_\mathrm{S}.$$
(26)
It is worth emphasizing once again that eqs. (18)-(26) assume $`\nu p_{\mathrm{gas}}`$ and $`\beta 1`$ at $`rr_{\mathrm{s}.\mathrm{g}.}`$; in view of eq. (17), the latter condition is expected to be satisfied at high accretion rates, $`\dot{M}1M_{}\mathrm{yr}^1`$.

### 3.2 Fragmentation

We assume that the initial sizes and masses of protostars forming in the disk are determined by Toomre (1964)’s dynamical instability. The shear stabilizes long wavelength perturbations from collapse, while the pressure stabilizes short wavelength modes. The most unstable mode has radial wave number $`k_{\mathrm{m}.\mathrm{u}.}=\pi G\mathrm{\Sigma }/c_s^2=(Qh)^1`$ so that

$$\lambda _{\mathrm{m}.\mathrm{u}.}(r_{\mathrm{s}.\mathrm{g}.})=\frac{2\pi }{k_{\mathrm{m}.\mathrm{u}.}}63\widehat{\kappa }(l_\mathrm{E}/\epsilon _{0.1})R_\mathrm{S}.$$
(27)
The numerical expression assumes $`\beta _{\mathrm{s}.\mathrm{g}.}1`$ so that $`h`$ is given by eq. (21), which is probably a good approximation for $`M_81`$ and $`l_\mathrm{E}1`$. The instability is axisymmetric, but fragments having comparable azimuthal and radial dimensions probably result. The corresponding mass at $`r=r_{\mathrm{s}.\mathrm{g}.}`$ is then of order

$$M_{\mathrm{m}.\mathrm{u}.}\mathrm{\Sigma }\lambda _{\mathrm{m}.\mathrm{u}.}^2300\alpha _{0.3}^1\widehat{\kappa }^{3/2}\mu (l_\mathrm{E}M_8/\epsilon _{0.1})^{5/2}M_{}.$$
(28)
A second way to estimate the fragment mass is to note that $`\beta `$ is approximately conserved under an adiabatic contraction, since as eq. (6) shows, the entropy depends only logarithmically on density at fixed $`\beta `$. Therefore one may set $`\beta _{}\beta (r_{\mathrm{s}.\mathrm{g}.})`$ in the relation eq. (3) for Eddington stellar models. This yields an initial stellar mass about four times larger than eq. (28) but with the same scalings. The star is likely to gain a great deal more mass by accretion, as we shall explain next.

## 4 Ultimate masses

In a strongly unstable case where $`Q<1`$, the disk might break up completely into closely spaced fragments. The subsequent evolution would then be dominated by collisional agglomeration. Let us consider instead a less unstable situation ($`Q1`$) so that the initial fragments are well separated and contain only a small fraction of the disk mass. Then a protostar gains mass mainly by accretion from the ambient gaseous disk.

We first consider whether the fragment (with mass given by eq. 28) will itself fragment into smaller masses. The usual condition for this to happen is that the cooling time be less than the the dynamical time. The initial density ($`\rho _{}`$) and optical depth ($`\tau _{}\kappa \rho _{}R_{}`$) of the fragment will be comparable to those of the ambient disk, and its radius $`R_{}\lambda _{\mathrm{m}.\mathrm{u}.}2h`$, the disk half-thickness. The initial dynamical time is $`(G\rho _{})^{1/2}(G\rho _{\mathrm{disk}})^{1/2}=Q\mathrm{\Omega }`$, and the cooling time of the fragment is of order the local thermal time of the disk, $`t_{\mathrm{th}}=\alpha ^1\mathrm{\Omega }`$ (Pringle, 1981). Therefore, provided $`\alpha Q<1`$ then $`t_{\mathrm{cool}}>t_{\mathrm{dyn}}`$. Furthermore, if the fragment contracts at constant mass, then $`t_{\mathrm{dyn}}R_{}^{3/2}`$, whereas $`t_{\mathrm{cool}}R_{}^1`$ if $`\kappa `$ is constant, so that $`t_{\mathrm{cool}}/t_{\mathrm{dyn}}R_{}^{5/2}`$, making further fragmentation increasingly unlikely as the radius shrinks. Even if subfragmentation were to begin, rapid accretion from the surrounding disk (see below) would probably drive the fragments together again by absorbing energy and angular momentum from the relative motion of the fragments. In short, the mass of the protostar is more likely to grow by accretion than to decrease by fragmentation. Nevertheless, it would be useful to confirm this by numerical simulations in two or three dimensions with allowance for cooling.

### 4.1 Accretion and the isolation mass

The accretion rate will be estimated by analogy with growth of terrestrial planets in a planetesimal disk (Lissauer, 1987). An important lengthscale in the latter problem is the Hill radius $`R_H(M_\mathrm{p}/3M_{})^{1/3}r`$, which is essentially the size of the Roche lobe surrounding the planet (mass $`M_\mathrm{p}`$) at an orbital distance $`r`$ from the central star (mass $`M_{}`$). At distances $`R_H`$ from the planet, motions of planetesimals are dominated by the gravitational field of the planet rather than the star. In a quasar disk, the role of the planet is played by the protostar (mass $`M_{}10^{25}M_{}`$), and that of the Sun, by the central black hole (mass $`M10^8M_{}`$). Thus

$$R_H\left(\frac{M_{}}{3M}\right)^{1/3}r.$$
(29)
Tidal torques exerted by the growing protostar will eventually open a gap in the disk, shutting off or at least reducing the rate of accretion. These torques, however, act only on gas that has encountered the protostar at least once. Material on an orbit whose semimajor axis differs from that of the protostar by $`R_H`$ stands a good chance of being accreted on its first passage. Let us therefore assume that the protostar accretes the entire annulus $`|rr_{}|f_HR_H`$, where $`f_HO(1)`$. (Lissauer infers $`f_H34`$ from simulations of planetesimal growth.) The mass of this annulus is larger than the initial mass (eq. 28) of the protostar by a large factor $`r/h`$ if $`Q1`$. Therefore, the initial width of the gap that is cleared is proportional to Hill radius of the mass accreted. This condition defines the isolation mass (Lissauer, 1987),

$$M_{\mathrm{iso}}=\frac{(4\pi r^2\mathrm{\Sigma }f_H)^{3/2}}{(3M)^{1/2}}0.96\times 10^5f_H^{3/2}\alpha _{0.3}^{6/5}\beta ^{6(1b)/5}\widehat{\kappa }^{3/10}\mu ^{6/5}(l_\mathrm{E}/\epsilon _{0.1})^{9/10}M_8^{14/5}r_3^{21/10}M_{}.$$
(30)
Numerically, if $`f_H1`$, this is comparable to the disk mass (eq. 22), but formally $`M_{\mathrm{iso}}M_{\mathrm{disk}}\times (M_{\mathrm{disk}}/M)^{1/2}`$. The dependence of $`M_{\mathrm{iso}}`$ with $`r`$ and $`M`$ is shown in Figures 1 and 2. The isolation mass is startlingly large and one wants to examine very critically the assumptions that have led to it.

The dynamics of a pressure-supported gas are somewhat different from those of a nearly collisionless swarm of planetesimals. Rafikov (2003) has recently demonstrated that the growth rate of a planetary embryo slows significantly at a mass $`M_{\mathrm{iso}}`$ because gravitational encounters with the embryo tend to “heat” the epicyclic velocity dispersion of the planetesimals and thereby decrease their rate of collisions with the embryo. In this respect a quasar disk is more favorable for rapid growth because the gas can cool radiatively. In fact, the cooling time is $`(\alpha \mathrm{\Omega })^1`$, which is probably short compared to the time between successive encounters with the protostar, $`2\pi (\mathrm{\Omega }f_HR_H/r)^1`$. Also, the accreting protostar probably nearly fills its Hill sphere during the growth phase (see below) so that it presents a relatively large cross section for “collisions” compared to a solid planetary embryo whose density is determined by chemical bonds.

On the other hand, whereas impacts of planetesimals on planetary embryos are perhaps fully inelastic, it is not immediately obvious that the accreting protostar can accept gas at the rate that it arrives *via* differential rotation. To address this question, we compare the rate at which mass enters the Hill sphere, the rate at which this material sheds its “spin” angular momentum, and the rate at which it sheds energy. Since masses much larger than (28) are concerned, we estimate these rates in the limit $`R_HH`$.

We begin by estimating the rate of growth of the Hill sphere. Material on an orbit at radial separation $`\mathrm{\Delta }r`$ from the center of the Hill sphere approaches it at azimuthal speed

$$v_{rel}=r\frac{d\mathrm{\Omega }}{dr}\mathrm{\Delta }r\frac{3}{2}\mathrm{\Omega }\mathrm{\Delta }r.$$

The surface density on such a streamline is $`\mathrm{\Sigma }`$, the background surface density of the disk, if the material is approaching the star for the first time, that is to say,

$$t<\frac{2\pi r}{|v_{rel}|}=\frac{4\pi }{3}\frac{r}{f_HR_H}\mathrm{\Omega }^1,$$
(31)
where $`t`$ is the time since accretion begins. We have set $`\mathrm{\Delta }r=f_HR_H`$ because if the inequality is satisfied at this separation, then it is satisfied at smaller separations, and more distant material is not likely to be much perturbed. Hence mass enters the Hill sphere at the rate

$$\frac{dM_H}{dt}\mathrm{\Sigma }\underset{rf_HR_H}{\overset{r+f_HR_H}{}}|v_{rel}|𝑑r\frac{3}{2}\mathrm{\Sigma }\mathrm{\Omega }f_H^2R_H^2.$$

Since $`R_Hr(M_H/3M)^{1/3}`$, this can be rewritten as

$$\frac{d}{dt}R_H=\frac{\mathrm{\Omega }\mathrm{\Sigma }r^3}{6M}=\mathrm{constant}.$$
(32)
(We use $`M_H`$ rather than $`M_{}`$ to denote the mass within the Hill sphere so as not to prejudge the question as to the timescale on which this mass concentrates into a quasi-spherical star rather than a centrifugally supported disk: see below.) Thus $`R_H`$ increases linearly with time until the two sides of the inequality (31) become equal, at which point accretion presumably slows or stops. This defines the asymptote of $`R_H`$ and hence the isolation mass:

$`R_H(M_{\mathrm{iso}})`$ $`=`$ $`\left({\displaystyle \frac{2\pi \mathrm{\Sigma }r^2}{9M}}\right)^{1/2}r,`$
$`M_{\mathrm{iso}}`$ $`=`$ $`{\displaystyle \frac{(2\pi f_H\mathrm{\Sigma }r^2)^{3/2}}{9M^{1/2}}}.`$
The mass is smaller than our previous estimate (30) by a factor $`6^{3/2}`$, which is an indication of the crudeness of both arguments. The approximate time at which this mass is achieved is found by interpreting eq. (31) as an equality with $`R_H`$ as given above:

$`t_{\mathrm{iso}}`$ $``$ $`{\displaystyle \frac{2\pi }{\mathrm{\Omega }}}\left({\displaystyle \frac{2M}{f_H^3\pi \mathrm{\Sigma }r^2}}\right)^{1/2}`$ (33)
$``$ $`210f_H^{3/2}\alpha _{0.3}^{2/5}\widehat{\kappa }^{1/10}\mu ^{2/5}(\epsilon _{0.1}/l_\mathrm{E})^{3/10}M_8^{2/5}r_3^{4/5}\mathrm{yr}.`$
Since $`M/\mathrm{\Sigma }r^2r/H`$ at $`Q=1`$, $`t_{\mathrm{iso}}`$ is longer than the local orbital period by $`(r/H)^{1/2}`$.

Having estimated the timescale on which the isolation mass accumulates within its Hill sphere, we now consider the timescale on which it contracts to radii $`R_H`$. Unless the contraction time is smaller than the accumulation time, it is questionable whether the Hill sphere can accept the full mass $`M_{\mathrm{iso}}`$ in the first place. On the other hand, a contraction time much shorter than $`t_{\mathrm{iso}}`$ would suggest fragmentation.

In order to contract, the mass $`M_H`$ must shed both energy and the angular momentum associated with rotation around its own axis *i.e.* spin. We will now show that the timescale for loss of angular momentum is likely to be shorter than that for loss of energy. Suppose the contrary: then the material within the Hill sphere settles into a centrifugally supported disk, which will be called the “stellar disk” to distinguish it from the larger disk that orbits the black hole. Let the mass-weighted average thickness of this disk be $`\overline{h}`$, presumed to be $`R_H`$ but not necessarily equal to $`H`$, the thickness of the ambient black-hole disk; let its median (half-mass) radius be $`\overline{R}`$; and let its median angular velocity around its own axis be $`\overline{\omega }`$. Then the viscous time in the stellar disk is $`(\alpha \overline{\omega })^1(\overline{R}/\overline{h})^2`$. Suppose provisionally that this is long compared to $`t_{\mathrm{iso}}`$, as it would be if $`\overline{R}R_H`$ and $`\overline{h}H`$. Then $`\overline{R}R_H`$ and $`\overline{\omega }\mathrm{\Omega }`$ as $`R_H`$ expands, since the specific angular momentum of the material falling into the Hill sphere is $`\mathrm{\Omega }R_H^2`$. \[In fact, the stellar disk should be no larger than the streamline-crossing radius, which is $`R_{\mathrm{sc}}0.4R_H`$ at these extreme mass ratios $`M_H:M10^3`$ (Paczyński, 1977). At and beyond $`R_{\mathrm{sc}}`$, angular momentum is efficiently removed by the tidal potential.\] Thus $`\overline{h}\overline{c}_s/\mathrm{\Omega }`$, where $`\overline{c}_s`$ is the median sound speed in the stellar disk. The cooling rate of the disk by vertical radiative diffusion is therefore

$$\overline{t}_{\mathrm{cool}}^1=f_1\frac{c}{\kappa \overline{\mathrm{\Sigma }}\overline{h}},$$
(34)
where $`f_1`$ is a positive dimensionless constant, $`c`$ the speed of light, $`\kappa `$ the opacity, and $`\overline{\mathrm{\Sigma }}`$ the median surface density. We have assumed that the stellar disk is radiation-pressure dominated, like the ambient black-hole disk, but if not, then the cooling time would be longer by a factor $`(1\overline{\beta })^1`$.

New material joining the disk does so through a shock that endows it with internal energy $`GM_H/\overline{R}`$ per unit mass. Thus the median internal energy per unit area of the circumstellar disk evolves as

$$\frac{d}{dt}\overline{\mathrm{\Sigma }}\overline{c}_s^2=f_1\frac{c}{\kappa \overline{h}}\overline{c}_s^2+f_2\frac{GM_H\dot{M}_H}{\overline{R}^3},$$
(35)
in which $`f_2`$ is another positive constant. From the discussion above, $`\overline{R}R_Ht`$, $`M_Ht^3`$, and therefore $`\dot{M}_Ht^2`$ and $`\overline{\mathrm{\Sigma }}M_H/R_H^2t`$. So if $`\overline{h}`$ were constant in time, then $`\overline{c}_s`$ would also be constant, and the three terms in eq. (35) would scale with time at different rates: the first two would be constant, but the third $`t^2`$. Thus the equation could not hold for all $`t`$. The solution is to let $`\overline{c}_st`$ so that the first and third terms balance, in which case $`\overline{h}R_H`$.

In short, a centrifugally supported stellar disk would not remain thin but would evolve toward a quasi-spherical, pressure-supported configuration as long as it is rapidly accreting from the ambient black-hole disk. Once $`\overline{h}R_H`$, the viscous time becomes $`(\alpha \mathrm{\Omega })^1`$ at most. In fact, if the “star” fills its Hill sphere, then the angular momentum of its outer parts is not conserved, while its inner parts may transfer angular momentum outward by global bar instabilities on their own dynamical timescale, which is $`\mathrm{\Omega }^1`$.

At this point, it appears that energy is more problematic than angular momentum. The estimate (34) suggests that the cooling time of a quasispherical configuration should scale as $`t^2`$, whereas $`R_H/\dot{R}_Ht`$. This would seem to imply that cooling cannot keep pace with accretion, and therefore that the star should fill or overflow its Hill sphere. But the matter is more subtle than this because of the peculiarities of $`\mathrm{\Gamma }4/3`$ polytropes.

For a completely nonrotating star of fixed mass, the contraction timescale is the Kelvin-Helmholtz time obtained from eqs. (5) and (3); putting $`R_{}=R_H`$,

$$t_{\mathrm{KH},H}53.0\alpha _{0.3}^{1/5}\widehat{\kappa }^{19/20}\mu ^{4/5}(l_\mathrm{E}/\epsilon _{0.1})^{3/20}M_8^{1/5}r_3^{13/20}\left(\frac{M_{}}{M_{\mathrm{iso}}}\right)^{1/6}\mathrm{yr},$$
(36)
which is somewhat shorter than the accumulation time $`t_{\mathrm{iso}}`$ \[eq. (33)\]. (Henceforth we write $`M_{}`$ rather than $`M_H`$ for the mass within the Hill sphere.) The critical difference between the two timescales is the factor $`\beta _{}M_{}^{1/2}`$ in eqs. (4) & (5), which reflects the nearly zero binding energy of a spherical high-mass star, so that only a small fraction of the internal energy need be radiated in order to cause the star to contract substantially. This factor does not occur in the vertical contraction of a disk, of course.

Furthermore, the star could contract even without radiative losses because it accretes material with negative total energy. It can be shown that when $`R_HH`$, so that fluid streamlines can be approximated by ballistic test-particle trajectories, material accretes onto the Hill sphere with a negative Jacobi constant per unit mass,

$$H_J\frac{1}{2}𝐯^2\frac{GM_{}}{|𝐫𝐫_{}|}\frac{3}{2}(\mathrm{\Omega }\mathrm{\Delta }r)^2,$$
(37)
where $`𝐫_{}`$ is the center of the star, $`\mathrm{\Omega }`$ its orbital angular velocity, and $`𝐯`$ is the test-particle velocity relative to the star in a frame rotating at $`\mathrm{\Omega }`$. (The Jacobi constant of the accreting streamlines cannot be arbitrarily negative, since $`H_J(9/2)(\mathrm{\Omega }R_H)^2`$ in order to cross the inner or outer Lagrange point.) Very close to the star, at $`|𝐫𝐫_{}|R_H`$, $`H_J`$ reduces to the usual kinetic-plus-potential energy in the potential of an isolated mass $`M_{}`$. Thus, the putatively spherical star accretes material that has negative total energy even without radiative losses. This would cause contraction on a timescale even shorter than (36) if rotation could be ignored altogether. However, we expect that the star will in fact rotate with a median angular velocity scaling in proportion to $`\mathrm{\Omega }`$ as it gains mass. Consequently, the star’s rotational energy cannot be entirely neglected unless it is smaller than the binding energy given by eq. (4), which is marginally unlikely since $`\beta _{}10^2`$ at $`M_{\mathrm{iso}}`$.

To recap, the growing star cannot cool fast enough to become a disk but will be more nearly spherical and therefore able to shed angular momentum quickly; on the other hand, the residual rotational energy, though small compared to its gravitational energy, seems likely to stabilize it against rapid contraction. The interplay between energy loss, angular-momentum loss, and radiation-pressure dominance (so that $`\mathrm{\Gamma }4/3`$) is so intricate that it will be difficult to decide whether the star overflows or detaches from its Hill sphere without elaborate calculations that are beyond the scope of this paper.

Again, it would be useful to have recourse to numerical simulation. We are not aware of any calculations for disks supported largely by radiation pressure, but Bate et al. (2003) and Lufkin et al. (2003) have recently simulated three-dimensional accretion from a protostellar disk onto embedded planets ranging in mass from that of earth ($`M_{}`$) to that of Jupiter ($`M_J`$). Both studies place their planets at $`r_J=5.2\text{au}`$ and follow the evolution for $`10^2`$ orbits in a disk with $`h/r=0.05`$. Bate *et al.* use an eulerian finite-difference method, whereas Lufkin *et al.* use smooth particle hydrodynamics. Neither includes cooling, but this is not an impediment to accretion because the former study removes removes the gas at each time step from the grid cells closest to the planet (well within the Hill sphere), while the latter uses an isothermal equation of state and includes the selfgravity of the gas, which allows indefinite compression of the planet’s atmosphere. The accretion rates measured in these simulations are indeed comparable, up to factors $`2`$, to what is predicted by eq. (33) with $`f_H=1`$, even though the most massive planets may exceed the isolation mass (eq. 30), $`M_{\mathrm{iso}}=(0.089,0.25)f_H^{3/2}M_J`$, corresponding to the surface densities $`\mathrm{\Sigma }(r_J)(0.23,0.46)M_J/r_J^2`$ in the two studies. (Bate *et al.* impose a partially cleared gap around the planet in their initial conditions; we have quoted $`\mathrm{\Sigma }`$ as it would be if there were no gap.)

Unfortunately for the purposes of this comparison, the Hill radii of the most massive planets in these simulations is scarcely larger than the disk thickness, so that the accretion is not very far into the two-dimensional regime contemplated by equation (33). To express this another way, the planetary isolation mass quoted above is comparable to the mass $`M_h`$ defined by $`R_H(M_h)=h`$, *viz.* $`M_h=3(h/r)^3M_{}0.38M_J`$. The ratio $`M_h/M_{\mathrm{iso}}`$ scales as $`(Qh/r)^{3/2}`$, which is of order unity for gaseous protoplanetary disks but very small at $`r10^3R_\mathrm{S}`$ in quasar disks because both $`h/r`$ and $`Q`$ are smaller there. Satellites $`M_h`$ probably cannot open a gap regardless of the viscosity of the disk (Lin & Papaloizou, 1993); they will undergo Bondi accretion provided that the accumulated gas cools rapidly enough. Thus, the applicability of the concept of isolation mass in quasar disks remains to be tested.

Once the isolation mass is achieved, or possibly even earlier, the star will detach from its Hill sphere and contract until it reaches the main sequence after a time of order the Kelvin time (eq. 9). It is important to note that, at least for our nominal parameters, the time required to reach the main sequence is longer than the time to accrete the isolation mass: therefore we do not expect nuclear burning to impede accretion. The formation time is also short compared to the viscous evolution time of the disk (eq. 24).

### 4.2 Orbital migration and tidal disruption

When fusion begins, the star may disintegrate *via* pulsational instability (§2.1). It would certainly be interesting to understand this stage: will the star disrupt completely or leave behind a remnant? How much will the disk be enriched? Since the answers to these questions are not available at present, we assume provisionally that the star survives the instability with much of its mass intact. If so, then it can be expected to undergo Type II radial migration: that is, while maintaining an annular gap in the disk by tidal torques, the star will drift inwards on the viscous timescale (eq. 24) by analogy with jovian planets in protostellar disks (Ward, 1997). For our nominal parameters, the latter timescale is shorter than the main-sequence lifetime (eq. 10). Thus the migration may be completed before the star leaves the main sequence.

The star will be tidally disrupted at an orbital radius where $`R_H`$ equals the main-sequence radius (eq. 8),

$$r_{\mathrm{tid}}12\left(\frac{\mu _{}}{\mu _{}}\right)^{0.95}\left(\frac{M_{}}{10^5M_{}}\right)^{0.14}M_8^{2/3}R_\mathrm{S},$$
(38)
which is comparable to the radius of the marginally stable orbit for a Schwarzschild black hole with $`M=10^9M_{}`$.

### 4.3 Capture

The concept of isolation mass is independent of the initial $`Q`$ of the disk, although it is true that $`M_{\mathrm{iso}}`$ increases with the overall mass of the disk \[eq. (30)\]. Thus mechanisms other than gravitational instability may provide the initial seed from which the isolation mass develops by accretion. In the primordial solar nebula, for example, where most likely $`Q10^2`$ at the present position of Jupiter, that planet is believed to have been initiated by collisional agglomeration of solid bodies. We do not expect solids in the non-self-gravitating parts of quasar disks, but there may be other ways to form seeds that might grow to the isolation mass. One such process is capture of pre-existing stars in collisions with the disk.

The capture process has been examined at length by Syer et al. (1991) and Artymowicz et al. (1993), and we have only a few remarks to add. In the following discussion, all disk and black-hole parameters are assumed equal to their nominal values in §3 except the orbital radius, for which we assume $`r_3<1`$ so that self-gravity is slight and radiation pressure dominant.

The first remark is that in these regions, stars may perhaps more often be destroyed by their encounters with the disk than captured. The energy that must be dissipated to circularize a typical stellar orbit exceeds the binding energy of a solar-type star by a factor $`10^3r_3^2`$. Most of the dissipation should occur in a bow shock driven into the tenuous disk gas rather than within the star. From §3.1, the peak postshock pressure is $`\rho v^210^8r_3^{8/5}\text{dyne cm}^2`$; at just $`10^9`$ of the star’s central pressure, this is only a gentle squeeze. However, except near the stagnation point, the postshock flow passes the star at $`10^4r_3^{1/3}\mathrm{km}\mathrm{s}^1`$ and at a temperature between $`10^6r_3^{2/3}\mathrm{K}`$ and $`10^9r_3^{8/3}\mathrm{K}`$ (depending whether the radiation field comes to equilibrium with the plasma; the mean free path for creation of soft photons by bremsstrahlung is comparable to the standoff distance of the shock), so that it may erode the stellar surface, which is not rigid but simply a contact discontinuity. Of course, in view of the very low density of the disk gas ($`10^9r_3^{3/5}\text{g cm}^3`$), the mass lost from the star in a single disk passage is bound to be very small, but so is the orbital energy; some $`\mathrm{\Sigma }R_{}^2/M_{}10^6r_3^{3/5}`$ passages and $`10^7r_3^{21/10}\mathrm{yr}`$ elapse before circularization. A small minority of stars on orbits of low inclination and low eccentricity may encounter the disk at less than their own surface escape speed. Clearly the capture process presents further opportunities for interesting hydrodynamical, or radiation-hydro, simulations; Armitage et al. (1996) have simulated the passage of giant stars through a disk and find that much of the envelope can be stripped in a single encounter, but we are not aware of any hydrodynamic simulations for main-sequence stars.

The second remark concerns the rate of accretion by the star once it settles into the disk. Syer et al. (1991) and Artymowicz et al. (1993) estimate this at the Bondi rate,

$$\dot{M}_{,\mathrm{B}}=\frac{4\pi (GM_{})^2}{c_\mathrm{s}^3}5\times 10^3r_3^{21/10}M_{}\mathrm{yr}^1,$$

provided only that the Hill radius is smaller than the disk thickness, *i.e.* $`M_{}300r_3^3M_{}`$. But accretion rate will be smaller if the gas is unable to cool sufficiently quickly; in other words, the rate should not exceed Eddington rate,

$$\dot{M}_{,\mathrm{E}}=\frac{4\pi c}{\kappa }R_{}1.0\times 10^3\frac{M_{}/R_{}}{M_{}/r_{}}M_{}\mathrm{yr}^1.$$

Each of the first several doublings of the stellar mass is likely to require $`10^3\mathrm{yr}`$, which may be gradual enough so that a gap will have formed by the time $`M_{}M_h`$.

## 5 Observable consequences

Exotic and luminous as stars of $`10^5M_{}`$ would be, they could easily have escaped attention. Such stars radiate at their Eddington luminosity, but that is smaller than the Eddington limit of the quasar itself by the mass ratio $`M_{}/M10^3`$. Their effective temperature on the main sequence is $`7\times 10^4\mathrm{K}`$, based on the radius (eq. 8), which is determined essentially by virial considerations; a detailed atmospheric model might give a somewhat different photospheric radius and temperature. This is about twice the inferred $`T_{\mathrm{eff}}`$ of bright QSOs (Malkan, 1983).

If the star has opened up a gap in the accretion disk, then this will affect the quasar’s spectrum, particularly when the gap is in the innermost regions. In fact, because of the short viscous time of the inner disk, we expect the gas inside the gap to be accreted relatively quickly, leaving an evacuated region between the central black hole and supermassive star, which is surrounded by what is effectively a circumbinary disk. The quasar will be relatively dim and red just before disruption or accretion of the star: indeed the source may not have been previously identified as a quasar, given the migration time of $`10^5\mathrm{yr}`$. This is the time when any supermassive star would be most easily observed. For favorable inclinations, the star’s spectral features would exhibit very fast orbital velocities of order several thousand $`\mathrm{km}\mathrm{s}^1`$ with periods greater than several days. The emission from the QSO disk, would gradually brighten and harden. Figure 3 shows an approximate comparison between the spectrum of the supermassive star and QSO disk as the inner gas radius decreases.

The merger of a supermassive star (or supermassive remnant) with the central black hole would be a strong source of gravitational waves. We estimate the “modified characteristic amplitude” for quadrupole gravitational waves as defined by Finn & Thorne (2000) from the inspiral of a compact source with a mass equal to the isolation mass formed at $`10^3R_s`$ from the central black hole (Figure 4). The trajectories in the frequency-amplitude diagram depend on the angular momentum of the central black hole, since the innermost stable orbit is much tighter for maximally rotating Kerr black holes: in the extreme case shown in Figure 4 with rotation parameter a=0.999, this radius is 1.18 $`GM/c^2`$, while it is $`6GM/c^2`$ for the Schwarzschild black hole. For the systems shown in Figure 4, a main sequence supermassive star would be disrupted before the differences due to angular momentum show themselves. Nevertheless, at and shortly after disruption, the waveform is bound to differ from that expected for an inspiraling black-hole companion. It would be interesting to try to predict these differences in some quantitative detail, as they might be the most direct way to confirm the existence of supermassive main-sequence stars in these disks.

The dependence of the isolation mass on the mass of the black hole means that we expect the supermassive stellar companions of more massive quasars to be more readily detected, although the signal is shifted towards lower frequencies. The cosmological redshift of the sources is important, as illustrated in Figure 4. The evolution of the gravitational wave signal from inspiral, should also be accompanied by a brightening and hardening of the quasar’s electromagnetic radiation, as described above, and may be marked by an outburst if the star is disrupted.

We can make a crude estimate of the upper limit to the event rate of supermassive star merger events from the known quasar population that may be seen by LISA. Using the quasar luminosity function of Boyle et al. (2000), we estimate that the space density of “typical” quasars (within a couple of magnitudes of the knee of the luminosity function) is $`10^6\mathrm{Mpc}^3`$ in a flat $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ cosmology. The comoving volume from $`0.5<z<1.5`$ is $`320\mathrm{Gpc}^3`$ and from $`1.5<z<2.5`$ is $`480\mathrm{Gpc}^3`$ so there are $`2.3,3.5\times 10^5`$ typical QSOs in these intervals, respectively, that can be seen over the whole sky. The supermassive star formation rate cannot be much greater than one per viscous time at $`r_{\mathrm{s}.\mathrm{g}.}`$ (and could be substantially less), which is about $`1\times 10^5\mathrm{yr}^1`$ per quasar. Thus the maximum event rate for LISA is of order a few supermassive star - black hole mergers per year. Of course the observed quasar luminosity function may underestimate the true number because of beaming effects.

## 6 Conclusions

We have argued that conditions in QSO accretion disks are likely to lead to the formation of massive stars. If the final masses of these stars are limited by the opening of a gap in the disk near the innermost radius where the disk is self-gravitating, then they would have $`10^5M_{}`$. This is much more massive than any stellar object currently known. A key point is that the main sequence lifetime of these stars is much longer than their formation time, and modestly longer than the time they take to migrate to the black hole. Thus the end point of stellar evolution does not prevent the existence of these supermassive stars, and indeed the evolution may be interrupted only by tidal disruption and merger near the central black hole. A second important point is that once a star has formed in the disk with initial mass of order the Toomre mass, $`100M_{}`$, then the subsequent accretion rates of surrounding disk material are much larger ($`M_{}\mathrm{yr}^1`$) than the mass loss rates caused by thermonuclear instabilities ($`<10^3M_{}\mathrm{yr}^1`$).

Radiative feedback from the forming massive star (from luminosity generated either by nuclear burning or gravitational contraction) will have little influence while conditions in the disk and accretion flow are optically thick. Gas streamlines will only be significantly perturbed as the stellar surroundings become optically thin, probably close to the time of gap formation, when the isolation mass has already been reached.

Though they are very luminous, supermassive stars would easily be lost in the glare of the quasar. Periodic modulations of the light curve due to such stars would have periods ranging from a few years to a few days but with millimagnitude amplitudes. Even if the quasar was in a quiet state, perhaps because of the disruption of its disk by the supermassive star, then the light from the central regions of the host galaxy might dominate. The best candidates for periodicity searches would be quasars that appear to lack a blue bump, as this might be due to a gap opened in the disk by such a star. Space-based gravitational wave detectors are promising tools in the search for supermassive stars, particularly since the signal should also be accompanied by a rapid brightening of the quasar as the star is disrupted.

Very recently, ideas similar to those presented here have been presented by Levin (2003). Although we have discussed our respective drafts, our two efforts have been independent. There is much we agree on, but there are some differences in assumptions and emphasis. Whereas Levin considers a disk that fragments into many stars of up to a few hundred solar masses, we consider the formation of a single dominant mass that accretes a significant fraction of the entire disk. Naturally we prefer our own scenario, but much more work will be required to decide which (if either) is correct. It is also possible that the outcome depends upon the initial conditions of the disk before fragmentation begins. Some insight into the early stages of fragmentation could be gained by numerical simulations of marginally self-gravitating disks; there have been many such studies, but none that we are aware of that include radiation-pressure dominance and radiative diffusion. Another difference between our studies is that and Levin considers parameters appropriate to the Galactic Center, whereas we focus on bright and massive QSOs. Collin & Zahn (1999a) and Collin & Zahn (1999b) have also discussed star formation in AGN disks recently, but they too consider stars which, though massive by ordinary standards, are much less than the isolation mass considered here.

Another area of theoretical uncertainty, which is perhaps even more challenging, is the nonlinear stability and evolution of extremely massive stars. In recent years, apart from some important contributions cited above, work on this topic has languished for lack of a clear astrophysical motivation. Convincing observational evidence for stars above $`10^2M_{}`$ does not exist. But, as we have already stressed, the physical conditions in a quasar accretion disk are very much more extreme than those of a conventional star-forming region, and we have argued that the self-gravity and short timescales of these disks favor the rapid assembly of truly supermassive protostars. It would be very interesting to revisit the question of whether such objects can survive on the main sequence.

We thank Yuri Levin and two anonymous referees for helpful discussions. JCT is supported by a Spitzer-Cotsen fellowship from Princeton University and by NASA grant NAG 5-10811. This work was supported in part by NSF grant AST-0307558 to JG.


format: text

text:
# Phase Motion of the Scalar 𝜋𝜋 Amplitudes in 𝐷⁺, 𝐷_𝑠⁺→𝜋⁻𝜋⁺𝜋⁺ Decays 11footnote 1Talk presented at the Scalar Meson Workshop, May 2003, SUNYIT, Utica, NY

## 1 Introduction

Fermilab experiment E791, with a full Dalitz plot analysis, showed strong evidence for the existence of light and broad scalar resonances in charm $`D^+`$ meson decay prl ; kappa . The $`\pi ^+\pi ^{}`$ resonance is compatible with the isoscalar meson $`\sigma (500)`$, and was observed in the Cabbibo-suppressed decay $`D^+\pi ^{}\pi ^+\pi ^+`$ . To get a good fit quality in this analysis, it was necessary to include an extra scalar particle, other than the well established dipion resonances pdg . For the new state, modeled by a Breit-Wigner amplitude, it was measured a mass and width of $`478_{23}^{+24}\pm 17`$ MeV/$`c^2`$ and $`324_{40}^{+42}\pm 21`$ MeV/$`c^2`$respectively . The $`D^+\sigma (500)\pi ^+`$ decay contribution is dominant, accounting for approximately half of this particular $`D^+\pi ^{}\pi ^+\pi ^+`$ decay. We found also evidence for a scalar $`K^{}\pi ^+`$ resonance, or $`\kappa `$, in the Cabibbo-allowed decay $`D^+K^{}\pi ^+\pi ^+`$ kappa . Further studies about $`\kappa `$ are discussed in this proceeding gobel .

In full Dalitz plot analyses, each possible resonance amplitude is represented by a Breit-Wigner function multiplied by angular distributions associated with the spin of the resonance. The various contributions are combined in a coherent sum with complex coefficients that are extracted from maximum likelihood fits to the data. The absolute value of the coefficients are related to the relative fraction of each contribution and the phases take into account the final state interaction (FSI) between the resonance and the third pion.

Due to the importance of this scalar meson in many areas of particle and nuclear physics, it is desirable to be able to confirm the amplitude’s phase motion in a direct observation, without having to assume, a priori, the Breit-Wigner phase approximation for low-mass and broad resonances ochs ; oller ; polosa . Recently, a method was proposed to extract the phase motion of a complex amplitude in three body heavy meson decays ad . The phase variation of a complex amplitude can be directly revealed through the interference in the Dalitz-plot region where it crosses with a well established resonant state, represented by a Breit-Wigner.

Here we begin with a simple example, showing that the AD method can be applied to extract the resonant phase motion of the scalar amplitude due to the resonance $`f_0(980)`$, using the same $`f_0(980)`$ resonance in the crossing channel in the Dalitz plot of the decay $`D_s^+\pi ^{}\pi ^+\pi ^+`$ using E791 data prlds . This example shows the ability of this method to extract the phase motion of an amplitude. Then we apply the AD method using the well known $`f_2(1270)`$ tensor meson in the crossing channel, as the base resonance, to extract the phase motion of the scalar low mass $`\pi \pi `$ amplitude in $`D^+\pi ^{}\pi ^+\pi ^+`$ , confirming the $`\sigma (500)`$ suggested by the E791 full Dalitz plot analysis prl .

## 2 Extracting $`f_0(980)`$ phase motion with the AD method.

From the original $`2\times 10^{10}`$ event data collected in 1991/92 by Fermilab experiment E791 from $`500GeV/c`$ $`\pi ^{}nucleon`$ interactions ref791 , and after reconstruction and selection criteria, we obtained the $`\pi ^{}\pi ^+\pi ^+`$ sample shown in Figure 1. To study the resonant structure of these three-body decays we consider the 1686 events with invariant mass between 1.85 and 1.89 GeV/c2, for the $`D^+`$ analysis and the 937 events with invariant mass between 1.95 and 1.99 GeV/c2 for the $`D_s^+`$. Figure 2(a) shows the Dalitz-plot for the $`D_s^+\pi ^{}\pi ^+\pi ^+`$ selected events and Figure 2(b) the Dalitz-plot for $`D^+\pi ^{}\pi ^+\pi ^+`$ events. The two axes are the squared invariant-mass combinations for $`\pi ^{}\pi ^+`$, and the plot is symmetrical with respect to the two identical pions.

We can see in Figure 2(a) the scalar $`f_0(980)`$ in $`s_{12}`$, the square invariant mass, crossing the $`f_0(980)`$ in $`s_{13}`$, forming an interference region around $`s_{13}=s_{12}=`$ 0.95GeV2. The AD method uses the interference region, between two crossing resonances, to extract the phase motion of one of them, and Final State Interaction (FSI) phase, provided that the second is represented by a Breit-Wigner ad . In fact we are using a Bootstrap approach; that is, using a well established resonance $`f_0(980)`$ in $`s_{12}`$ to extract its phase motion in $`s_{13}`$. It is a nice and didactic example to show the ability of this method to extract the phase motion of an amplitude and the FSI phase within the E791 data sample.

The coherent amplitude to describe the crossing between a well known scalar resonance, represented by a Breit-Wigner in $`s_{12}`$, and a complex amplitude under study in $`s_{13}`$ in a limited region of the phase space, where we can neglect any other contributions, is given by:

$$𝒜(s_{12},s_{13})=a_R𝒲(s_{12})+a_s/(p^{}/\sqrt{s}_{13})sin\delta (s_{13})e^{i(\delta (s_{13})+\gamma )}$$
(1)
$`p^{}/\sqrt{s}_{13}`$ is a phase space factor to make this description compatible with $`\pi \pi `$ scattering, $`\gamma `$ is the final state interaction (FSI) phase difference between the two amplitudes, $`a_R`$ and $`a_s`$ are respectively the real magnitudes of the resonance and the under-study complex amplitude. Finally $`sin\delta (s_{13})e^{i\delta (s_{13})}`$ represents the most general amplitude for a two-body hadronic interaction.

The Breit Wigner distribution is given by:

$$𝒲=\frac{m_0\mathrm{\Gamma }_0}{m_0^2sim_0\mathrm{\Gamma }(m)}$$

Taking the amplitude square of Equation 1 we get:

$`𝒜(s_{12},s_{13})^2=a_R^2𝒲_{f_0(980)}(s_{12})^2+a_s^2/p^2/s_{13}sin^2\delta (s_{13})`$
$`+{\displaystyle \frac{2a_Ra_sm_0\mathrm{\Gamma }_0sin\delta (s_{13})/(p^{}/\sqrt{s}_{13})}{(m_0^2s_{12})^2+m_0^2\mathrm{\Gamma }^2(s_{12})}}\times [(m_0^2s_{12})cos(\delta (s_{13})+\gamma )+m_0\mathrm{\Gamma }_0sin(\delta (s_{13})+\gamma )]`$ (2)
Since the Breit-Wigner is approximately symmetrical around $`m_0`$ as seen in Figure 3 (the asymmetries would come from $`\mathrm{\Gamma }(m)`$, and is negligible for the narrow $`f_0(980)`$). We can divide our $`f_0(980)`$ mass distribution in two pieces, one for $`m_0+ϵ`$ and the other with $`m_0ϵ`$. From Equation 2 and noticing that the pure Breit-Wigner term will cancel we can write:

$`𝒜(m_0^2+ϵ,s_{13})^2𝒜(m_0^2ϵ,s_{13})^2={\displaystyle \frac{4a_sa_R/(p^{}/\sqrt{s}_{13})ϵm_0\mathrm{\Gamma }_0}{ϵ^2+m_0^2\mathrm{\Gamma }_0^2}}(sin(2\delta (s_{13})+\gamma )sin\gamma )`$ (3)
Only the real part of the interference term in Equation 2 remains.

To extract the phase motion of the scalar amplitude in $`s_{13}`$ through the $`f_0(980)`$ in $`s_{12}`$, represented by a Breit-Wigner, we took the events in $`s_{12}`$ between 0.7 and 1.2 GeV2 and divided them into two bins, as presented in Figure 3. The $`s_{13}`$ distribution for the events of the $`s_{12}`$ region integrated between 0.95 and 1.2 GeV2, is shown in Fig. 4a and the same in Figure 4b for events integrated between 0.7 and 0.95 GeV2.

We can see that the peaks in these two plots are in different $`s_{13}`$ positions. The subtraction of these distributions, corresponds to the integration of Equation 3, that we can write as:

$`\mathrm{\Delta }{\displaystyle 𝒜^2}={\displaystyle _{m_{eff}^2}^{m_{eff}^2+ϵ}}𝒜(s_{12},s_{13})^2𝑑s_{12}{\displaystyle _{m_{eff}^2ϵ}^{m_{eff}^2}}𝒜(s_{12},s_{13})^2𝑑s_{12}`$
$`𝒞(sin(2\delta (s_{13})+\gamma )sin\gamma )`$ (4)
where $`𝒞`$ is a constant factor coming from the constant and integrated factors of Equation 3, to be determined from data. The variation of the phase space in the integral was considered negligible for the $`f_0(980)`$ resonance. $`\mathrm{\Delta }𝒜^2`$ directly reflects the behavior of $`\delta (s_{13})`$. A constant $`\mathrm{\Delta }𝒜^2`$ would imply constant $`\delta (s_{13})`$. This would be the case for a non-resonant contribution. The same way a slow phase motion will produce a slowly varying $`\mathrm{\Delta }𝒜^2`$ and a full resonance phase motion produces a clear signature in $`\mathrm{\Delta }𝒜^2`$ with the presence of zero, maximum and minimum values.

The subtracted distribution, corresponding to Equation 4, is shown in Figure 5. There is a significant difference between the minimum (bin3) and maximum (bin4) of $`\mathrm{\Delta }𝒜^2`$.

We can see in Equation 4 that the zeros occur when $`\delta (s_{13})`$ = $`0^0`$, $`180^0`$ or $`\pi /2\gamma `$. In Figure 5 we can see a zero at $`s_{13}`$ near 0.5GeV2, another one at $`s_{13}`$ = 0.95GeV2 and a third zero near 1.4GeV2. Assuming $`\delta (s_{13})`$ is an analytical function of $`s_{13}`$, Equation 4 allow us to set the two following conditions at the maximum and minimum values of $`\mathrm{\Delta }𝒜^2`$ respectively:

$$\mathrm{\Delta }𝒜_{max}^2sin(2\delta (s_{13})+\gamma )=1$$
(5)
$$\mathrm{\Delta }𝒜_{min}^2sin(2\delta (s_{13})+\gamma )=1$$
(6)
with these two conditions we get $`𝒞`$ and $`\gamma `$, calculated from the maximum and minimum values of the $`\mathrm{\Delta }𝒜^2`$ distribution in Figure 5:

$$𝒞=(\mathrm{\Delta }𝒜_{max}^2\mathrm{\Delta }𝒜_{min}^2)/2$$
(7)
$$\gamma =sin^1(\frac{\mathrm{\Delta }𝒜_{max}^2+\mathrm{\Delta }𝒜_{min}^2}{\mathrm{\Delta }𝒜_{min}^2\mathrm{\Delta }𝒜_{max}^2})$$
(8)
From Figure 5 and using the equations above, we measure $`\gamma =0.15\pm 0.31`$, that is compatible with zero, as should be since we are crossing the same resonances with, of course, the same final state interaction phase.

With the above conditions we solve Equation 4 for $`\delta (s_{13})`$:

$$\delta (s_{13})=\frac{1}{2}(sin^1(\frac{1}{𝒞}\mathrm{\Delta }𝒜(s_{13})^2+sin(\gamma ))\gamma )$$
(9)
Assuming that $`\delta (s_{13})`$ is an increasing function of $`s_{13}`$, we can extract directly the $`\delta (s_{13})`$ value from each bin of Figure 5, creating the $`f_0(980)`$ phase motion shown in Figure 6. The errors in the plot were produced by generating statistically compatible experiments, allowing each bin of $`_{m_0^2}^{m_0^2+ϵ}𝒜(s_{12},s_{13})^2`$ (Figure4a) and $`_{m_0^2ϵ}^{m_0^2}𝒜(s_{12},s_{13})^2`$ (Figure4b) to fluctuate randomly following a Poisson law. We then solve the problem for each ”experiment”. The error in each bin of $`\delta (s_{13})`$ will be the RMS of the distributions generated by the ”experiments”.

From Figure 6 we can see what one could expect, that is the scalar amplitude near 970GeV with a phase motion of about $`180^0`$ degrees. This example demonstrates the ability of AD method to extract the phase motion of an amplitude with E791 statistics.

## 3 Extracting the scalar low mass $`\pi \pi `$ amplitude phase motion with the AD method.

In the preceding section, we showed how to apply the AD method to extract the phase motion of an amplitude, from nonleptonic charm-meson three-body decay. Here we apply the same method to extract the phase motion of the scalar low-mass $`\pi \pi `$ amplitude in $`D^+\pi ^{}\pi ^+\pi ^+`$ decay, where we previously found strong experimental evidence for the existence of a light and broad isoscalar resonance prl . To start this analysis, we have to decide what is the best well-known resonance to be used for crossing the low mass amplitude under study. Taking a look at Figure 2b we can see the signature of three resonances that in principle could be used, the $`\rho (770)`$, $`f_0(980)`$ and $`f_2(1270)`$. In fact, the E791 analysis of this Dalitz plot found a significant contribution from these three resonances in $`D^+\pi ^{}\pi ^+\pi ^+`$ decay prl . Since this $`D^+`$ decay is symmetric for the exchange of the $`\pi ^+`$ meson, each resonance in $`s_{12}`$ is present also in $`s_{13}`$. Then if we use $`\rho (770)`$ as the base resonance in $`s_{12}`$, we have also the presence of the $`\rho (770)`$ in same mass square distribution of the $`\sigma (500)`$ in $`s_{13}`$. The proximity of the $`\rho (770)`$ with the $`\sigma (500)`$, both broad resonances, creates an overlap between them such that we are not able to separate the phase motion of one from the other. We could use the $`f_0(980)`$ as a base resonance, but again the presence of the $`\rho (770)`$ overlapping with the $`\sigma (500)`$ creates the same problem.

There remains only the tensor meson $`f_2(1270)`$ candidate at $`m_0^2=1.61GeV^2`$, which is placed where the $`\rho (770)`$ contribution reaches a minimum due the angular distribution in the middle of the Dalitz plot, as we can see from the $`D^+\rho (770)\pi ^+`$ decay Monte Carlo simulation shown in Figure 7.

The amplitude for the crossing of the $`f_2(1270)`$ in $`s_{12}`$ and the complex amplitude under study in $`s_{13}`$ is given in the same way as in Equation 1:

$`𝒜(s_{12},s_{13})=a_R𝒲_{f_2(1270)}(s_{12}){}_{}{}^{j=2}_{f_2(1270)}^{}(s_{12},s_{13})+`$ (10)
$`+a_s/(p^{}/\sqrt{s}_{13})sin\delta (s_{13})e^{i(\delta (s_{13})+\gamma )}`$
where $`{}_{}{}^{j=2}_{f_2(1270)}^{}(s_{12},s_{13})`$ is the angular function for the $`f_2(1270)`$ tensor resonance. The amplitude under study represents the scalar low mass $`\pi \pi `$ amplitude in a limited region of the phase space, where we can neglect the other amplitude contributions.

Both the width $`\mathrm{\Gamma }(s_{12})`$ and the angular function $`{}_{}{}^{j=2}_{f_2(1270)}^{}`$ from this resonance produce asymmetries in $`s_{12}`$ and consequently we can not use the nominal $`f_2(1270)`$ mass to divide our sample into two slices, as we did for the $`f_0(980)`$ example. So we performed a Monte Carlo study to determine the effective mass we must use. The $`s_{12}`$ Monte Carlo projection of the $`f_2(1270)`$ in $`D^+\pi ^{}\pi ^+\pi ^+`$ decay is shown in Fig. 8. We can see the asymmetry created around the nominal $`f_2(1270)`$ mass value due to $`\mathrm{\Gamma }(s_{12})`$ and $`{}_{}{}^{j=2}_{f_2(1270)}^{}`$ contributions to the amplitude.

Here we require an effective mass squared ($`m_{eff}`$), such that the number of events integrated between $`m_{eff}^2`$ and $`m_{eff}^2+ϵ`$ is equal, by construction, to the number of events integrated between $`m_{eff}^2`$ and $`m_{eff}^2ϵ`$. We choose, using the $`f_2(1270)`$ Monte Carlo distribution, a mass of $`m_{eff}^2=1.535GeV^2`$, within $`\pm 0.26GeV^2`$ 222Within this mass region, the amount of $`\rho (770)`$ events was estimate to be around 5%, in such way that we can write:

$$_{m_{eff}^2}^{m_{eff}^2+ϵ}𝒲_{f_2(1270)}(s_{12}){}_{}{}^{j=2}_{f_2(1270)}^{}^2𝑑s_{12}=_{m_{eff}^2ϵ}^{m_{eff}^2}𝒲_{f_2(1270)}(s_{12}){}_{}{}^{j=2}_{f_2(1270)}^{}^2𝑑s_{12}$$
(11)
The effective mass squared $`m_{eff}`$ and the separation between $`m_{eff}^2+ϵ`$ (red) and $`m_{eff}^2ϵ`$ (blue) are shown in Figure 8.

The $`{}_{}{}^{j=2}_{f_2(1270)}^{}`$ function in $`s_{13}`$ is presented in Figure 9333 Since we divided our data sample by this function, we represent this function in a histogram with the same binning of data.. The distribution between $`m_{eff}^2`$ and $`m_{eff}^2+ϵ`$ is shown in Fig. 9a, for events between $`m_{eff}^2`$ and $`m_{eff}^2ϵ`$ in Fig. 9b. We can see that these two plots are slightly different. However we considered the approximation $`{}_{}{}^{j=2}_{f_2(1270)}^{+}(s_{13})^{j=2}_{f_2(1270)}^{}(s_{13})`$ and take an average function $`{}_{}{}^{j=2}\overline{}_{f_2(1270)}^{}(s_{13})`$. Another important effect, that we had to take into account, is the zero of this function at $`s_{13}=0.48GeV^2`$. Below we discuss the consequences of that in the AD method.

With the above considerations about the $`f_2(1270)`$ in $`s_{12}`$ and $`s_{13}`$ we can write the integrated amplitude-square difference as:

$`\mathrm{\Delta }{\displaystyle 𝒜^2}={\displaystyle _{m_{eff}^2}^{m_{eff}^2+ϵ}}𝒜(s_{12},s_{13})^2𝑑s_{12}{\displaystyle _{m_{eff}^2ϵ}^{m_{eff}^2}}𝒜(s_{12},s_{13})^2𝑑s_{12}`$
$`𝒞(sin(2\delta (s_{13})+\gamma )sin\gamma ){}_{}{}^{j=2}\overline{}_{f_2(1270)}^{}(s_{13})/(p^{}/\sqrt{s}_{13})`$ (12)
This Equation is similar to Equation 4, with an extra angular function term $`{}_{}{}^{j=2}\overline{}_{f_2(1270)}^{}(s_{13})`$ 444 For short we shall use, from here on $`{}_{}{}^{j=2}\overline{}_{f_2(1270)}^{}(s_{13})=\overline{}`$ and $`p^{}/\sqrt{s}_{13}=p^{}`$..

The background and the acceptance are similar between $`m_{eff}^2`$ and $`m_{eff}^2+ϵ`$ and $`m_{eff}^2`$ and $`m_{eff}^2ϵ`$. Since we are subtracting these two distributions, we do not take into account these effects in our analysis.

The $`𝒜^2`$ in $`s_{13}`$, for events integrated in $`s_{12}`$ between $`m_{eff}^2=1.535GeV^2`$ and $`m_{eff}^2+ϵ`$ and $`m_{eff}^2`$ and $`m_{eff}^2ϵ`$, with $`ϵ=0.26GeV^2`$ are presented in Figure 10a and b respectively.

Subtracting these two histograms, in the same way we did for the $`f_0(980)`$ example, gives the $`\mathrm{\Delta }𝒜^2`$ of the Equation 12. The result is shown in Figure 11.

Here we can not extract directly the phase motion from Figure 11, as we did for the $`f_0(980)`$ example using the conditions 5 and 6. We have to divide the $`\mathrm{\Delta }𝒜^2`$ by $`\overline{}`$ (average of the distributions in Figures 9a and b) and multiply by $`p^{}`$, since phase space here is an important effect. By doing this the only $`s_{13}`$ dependence of the right hand side of Equation 12 is in $`\delta (s_{13})`$. However, as we could see in Figure 9, there is a zero about 0.48GeV2 in the angular function, which means a singularity around this value in $`\mathrm{\Delta }𝒜^2/\overline{}`$. To avoid this singularity, we first produced a binning in such a way that the singularity is placed in the middle of one bin. In Figure 12 we show the $`\mathrm{\Delta }𝒜^2`$ by $`\overline{}`$ distribution. We can see that the 6th bin (around 0.48GeV2), has a huge error, that corresponds to the bin of the singularity. Due to the singularity we decided not to use this region (bin 6) further in this analysis. The consequences of this choice are going to be taken care of in systematic error studies. In any case, the singularity could only affect the position of the minimum of Figure 12. It does not interfere with the general feature of starting at zero, having statistically significant maximum and minimum values, and coming back to zero, indicating a strong phase variation. Bins 2 and 5 are respectively the maximum and minimum value of $`\mathrm{\Delta }𝒜^2p^{}/\overline{}`$ of Figure 12 where we use the Equation 5 and 6 conditions.

With the same assumptions used for the $`f_0(980)`$, that is $`\delta (s_{13})`$ is an analytical and increasing function of $`s_{13}`$, and using Equation 7, 8 and 9 (multiplied by $`p^{}`$ and divided by $`\overline{}`$), we can extract $`\gamma `$ and $`\delta (s_{13})`$ from Figure 12. For the FSI phase we found $`\gamma =3.26\pm 0.33`$, that is somewhat bigger than found by the E791 full Dalitz-plot analysis prl  ($`\gamma _{Dalitz}=2.59\pm 0.19`$). The fact that we used the effective mass for the $`f_2(1270)`$ = 1.535 GeV2 instead of the nominal mass is responsable for the shift observed in the relative phase. To verify this statement we generated 1000 samples of fast-MC, with only two amplitudes, $`f_2(1270)`$ and $`\sigma (500)`$. For both we used Briet-Wigner functions and the E791 parameters. We generated the phase difference of 2.59 rad, measured by the E791. For these 1000 samples, we measure $`\gamma `$ using the method presented here. The result has a mean value of 3.07. We can say that the difference between the generated and measured $`\gamma `$ value is a correction factor due to the use of an effective mass. Using this offset factor ( 2.59 - 3.07 = -0.48) we correct the measurement $`\gamma =3.26\pm 0.33`$ to $`\gamma _{corr}=2.78\pm 0.33`$. So the observed $`\gamma `$ difference between Dalitz analysis and the $`\gamma _{corr}`$ is in good agreement, with a difference below one standard deviation.

The $`\delta (s_{13})`$ was extracted bin by bin, with the same approach for the errors used in the $`f_0(980)`$ example, and we got the distribution shown in Figure 13 555 In Figure 13 there is no the 6th bin because of the singularity we mentioned above.. We can see a strong phase variation of about 1800 around the mass for the $`\sigma (500)`$, showing a phase motion compatible with a resonance.

## 4 Conclusions

We showed that the AD method can be applied to E791 data to extract the phase motion of the resonance $`f_0(980)`$ in the Dalitz plot of the decay $`D_s^+\pi ^{}\pi ^+\pi ^+`$ . This example demonstrates the ability of this method to extract the phase motion of a resonance amplitude.

Preliminary E791 results present a direct and model-independent approach, obtained with the AD method, and confirms our previous result of the evidence of an important contribution of the isoscalar $`\sigma (500)`$ meson in $`D^+\pi ^{}\pi ^+\pi ^+`$ decay prl . We use the well known $`f_2(1270)`$ tensor meson in the crossing channel, as the base resonance, to extract the phase motion of the low mass $`\pi \pi `$ scalar amplitude. We obtain a $`\delta (s_{13})`$ variation of about 1800 consistent with a resonant $`\sigma (500)`$ contribution. We also obtain good agreement between the FSI $`\gamma _{corr}`$ observed with AD method and the $`\gamma `$ observed in the full Dalitz plot analysis.

We wish to thank the Scalar Meson Workshop organizers, specially A. Fariborz, for the invitation and to offer us a pleasent workshop in Utica. We would like to thank also Wolfgang Ochs for suggestions and important comments.