text:
# 1 Introduction

## 1 Introduction

The radiative decays of low mass vector mesons into two neutral pseudoscalars, $`VP^0P^0\gamma `$, are known to be a useful tool to investigate the complicated dynamics governing meson physics around 1 GeV. Particularly interesting are those decays proceeding mainly by the exchange of scalar resonances because of the unconventional  and controversial nature of these states, as exemplified, for instance, in the very recent analysis by Close and Törnqvist . Experimental data on the $`\varphi \pi ^0\pi ^0\gamma `$ decays where, after the emission of a photon, a $`K^+K^{}`$ pair rescatters into $`\pi ^0\pi ^0`$ in a process dominated by the exchange of isoscalar scalar resonances in the $`s`$-channel have thus been watched for with great interest .

The first measurements of this $`\varphi \pi ^0\pi ^0\gamma `$ decay have been reported by the SND and CMD-2 Collaborations. For the branching ratio, they obtain

$$B(\varphi \pi ^0\pi ^0\gamma )=\{\begin{array}{cc}(1.22\pm 0.10\pm 0.06)\times 10^4\hfill & \text{SND }\text{[5]}\hfill \\ (0.92\pm 0.08\pm 0.06)\times 10^4\hfill & \text{CMD-2 }\text{[6]}\hfill \end{array}$$
(1)
for $`m_{\pi ^0\pi ^0}>700`$ MeV in the latter case. More recently, the KLOE Collaboration has measured

$$B(\varphi \pi ^0\pi ^0\gamma )=(1.09\pm 0.03\pm 0.05)\times 10^4,$$
(2)
in agreement with (1) but with a considerably smaller error. In all the cases, the spectrum is clearly peaked at $`m_{\pi ^0\pi ^0}970`$ MeV, as expected from an important $`f_0(980)`$ contribution. This and other radiative vector meson decays are now being further investigated at the Frascati $`\varphi `$-factory DA$`\mathrm{\Phi }`$NE  with higher accuracy (see for example Ref. ).

On the theoretical side, the $`\varphi \pi ^0\pi ^0\gamma `$ decay has been considered by a number of authors . Early calculations of the vector meson dominance (VMD) amplitude for these processes, i.e. the contributions proceeding through the decay chain $`\varphi \pi ^0\rho ^0\pi ^0\pi ^0\gamma `$, were summarized in Ref. . The $`\varphi `$-decay vertex is forbidden by the Zweig rule and the width and branching ratios predicted by VMD, $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}}=51`$ eV and $`B(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}}=0.12\times 10^4`$ , are found to be substantially smaller than the experimental results quoted in Eqs. (1) and (2). Moreover the VMD spectrum peaks at low $`m_{\pi \pi }`$ values well below the $`f_0(980)`$ resonance mass.

The possibility of an enhancement in this branching ratio through the $`\varphi K^+K^{}(\gamma )\pi ^0\pi ^0\gamma `$ mechanism  was pointed out in Ref.  and further discussed in Ref.  in a Chiral Perturbation Theory (ChPT) context enlarged to include on-shell vector mesons. This formalism gives well defined predictions for the various $`VP^0P^0\gamma `$ decays in terms of $`P^+P^{}P^0P^0`$ rescattering amplitudes, which are easily calculated in strict ChPT, and a loop integral over the intermediate $`P^+P^{}`$ pair. In this approach, the $`\varphi \pi ^0\pi ^0\gamma `$ decay is dominated by (Zweig rule allowed) kaon loops leading to $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )_\chi =224`$ eV, while (Zweig rule forbidden) pion loop contributions are strongly suppressed. The corresponding spectrum peaks at moderate values of $`m_{\pi \pi }`$ but the interference between this kaon loop contribution and the previous VMD amplitude turns out to be quite small leading globally to $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}+\chi }=269`$ eV and $`B(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}+\chi }=0.61\times 10^4`$ , which are still below the experimental results quoted in Eqs. (1) and (2). Additional contributions are thus certainly required. The most natural candidates for closing this gap between theory and experiment are the contributions coming from the exchange of scalar resonances, such as the well established $`f_0(980)`$ and the more controversial $`\sigma (500)`$ (or $`f_0`$(400–1200)) mesons .

The purpose of this paper is to study the effects of the low mass scalar states in the $`\varphi \pi ^0\pi ^0\gamma `$ decays following the ChPT inspired context introduced in Refs.  to account similarly for the scalar exchange contributions to $`\varphi \pi ^0\eta \gamma `$ and $`\rho ,\omega \pi ^0\pi ^0\gamma `$. In this context, one takes advantage of the common origin of ChPT and the L$`\sigma `$M to improve the chiral loop predictions for $`VP^0P^0\gamma `$ exploiting the complementarity of both approaches for these specific processes. As a result, a simple analytic amplitude, $`𝒜(\varphi \pi ^0\pi ^0\gamma )_{\text{L}\sigma \text{M}}`$, will be obtained which includes the effects of the scalar meson poles and also shows the appropriate behaviour expected from ChPT at low dipion invariant masses. There also exists the contribution to $`\varphi \pi ^0\pi ^0\gamma `$ coming from the previously mentioned vector meson exchange. This VMD amplitude, $`𝒜(\varphi \pi ^0\pi ^0\gamma )_{\text{VMD}}`$, is well known and scarcely interesting but has to be added to $`𝒜(\varphi \pi ^0\pi ^0\gamma )_{\text{L}\sigma \text{M}}`$, i.e. to the relevant amplitudes containing the scalar meson effects, in order to compare with available and forthcoming data. We will conclude that data on the $`\varphi \pi ^0\pi ^0\gamma `$ channel can contribute decisively to improve our knowledge on the $`f_0(980)`$ scalar resonance, as well as on some specific features of its elusive and controversial partner $`\sigma (500)`$.

## 2 Chiral loop contributions to $`\varphi \pi ^0\pi ^0\gamma `$

The vector meson initiated $`VP^0P^0\gamma `$ decays cannot be treated in strict Chiral Perturbation Theory (ChPT). This theory has to be extended to incorporate on-shell vector meson fields. At lowest order, this may be easily achieved by means of the $`𝒪(p^2)`$ ChPT Lagrangian

$$_2=\frac{f^2}{4}D_\mu U^{}D^\mu U+M(U+U^{}),$$
(3)
where $`U=\mathrm{exp}(i\sqrt{2}P/f)`$, $`P`$ is the usual pseudoscalar octet matrix and, at this order, $`M=\text{diag}(m_\pi ^2,m_\pi ^2,2m_K^2m_\pi ^2)`$ and $`f=f_\pi =92.4`$ MeV. The covariant derivative, now enlarged to include vector mesons, is defined as $`D_\mu U=_\mu UieA_\mu [Q,U]ig[V_\mu ,U]`$ with $`Q=\text{diag}(2/3,1/3,1/3)`$ being the quark charge matrix and $`V_\mu `$ the additional matrix containing the nonet of vector meson fields. We follow the conventional normalization for the vector nonet matrix such that the diagonal elements are $`(\rho ^0+\omega _0)/\sqrt{2},(\rho ^0+\omega _0)/\sqrt{2}`$ and $`\varphi _0`$, where $`\omega _0=(u\overline{u}+d\overline{d})/\sqrt{2}`$ and $`\varphi _0=s\overline{s}`$ stand for the ideally mixed states. The physical $`\varphi `$-field is approximately $`\varphi _0+ϵ\omega _0`$, where $`ϵ+0.059\pm 0.004`$ accounts for the $`\omega `$-$`\varphi `$ mixing angle in the flavour basis .

One easily observes that there is no tree-level contribution from the Lagrangian (3) to the $`\varphi \pi ^0\pi ^0\gamma `$ amplitude and that, at the one-loop level, one needs to compute the set of diagrams shown in Ref. . We do not take into account pion loop contributions here since they proceed only through the small $`ϵ`$ piece of the $`\varphi `$-field (because of the Zweig rule) which is further suppressed by $`G`$-parity, as in $`\omega \pi ^+\pi ^{}`$. These pion loop contributions are known to be negligible as compared to those from kaon loops proceeding through the Zweig rule allowed piece, $`\varphi _0=s\overline{s}`$, not restricted by $`G`$-parity arguments .

A straightforward calculation leads to the following finite amplitude for $`\varphi (q^{},ϵ^{})\pi ^0(p)\pi ^0(p^{})\gamma (q,ϵ)`$ (see Ref.  for further details):

$$𝒜(\varphi \pi ^0\pi ^0\gamma )_\chi =\frac{eg_s}{2\pi ^2m_{K^+}^2}\{a\}L(m_{\pi ^0\pi ^0}^2)\times 𝒜(K^+K^{}\pi ^0\pi ^0)_\chi ,$$
(4)
where $`\{a\}=(ϵ^{}ϵ)(q^{}q)(ϵ^{}q)(ϵq^{})`$, $`m_{\pi ^0\pi ^0}^2s(p+p^{})^2=(q^{}q)^2`$ is the invariant mass of the final dipion system and $`L(m_{\pi ^0\pi ^0}^2)`$ is the loop integral function defined as

$$\begin{array}{cc}\hfill L(m_{\pi ^0\pi ^0}^2)& =\frac{1}{2(ab)}\frac{2}{(ab)^2}\left[f\left(\frac{1}{b}\right)f\left(\frac{1}{a}\right)\right]\hfill \\ & +\frac{a}{(ab)^2}\left[g\left(\frac{1}{b}\right)g\left(\frac{1}{a}\right)\right]\hfill \end{array}$$
(5)
with

$$\begin{array}{c}f(z)=\{\begin{array}{cc}\left[\mathrm{arcsin}\left(\frac{1}{2\sqrt{z}}\right)\right]^2\hfill & z>\frac{1}{4}\hfill \\ \frac{1}{4}\left(\mathrm{log}\frac{\eta _+}{\eta _{}}i\pi \right)^2\hfill & z<\frac{1}{4}\hfill \end{array}\hfill \\ g(z)=\{\begin{array}{cc}\sqrt{4z1}\mathrm{arcsin}\left(\frac{1}{2\sqrt{z}}\right)\hfill & z>\frac{1}{4}\hfill \\ \frac{1}{2}\sqrt{14z}\left(\mathrm{log}\frac{\eta _+}{\eta _{}}i\pi \right)\hfill & z<\frac{1}{4}\hfill \end{array}\hfill \end{array}$$
(6)
and $`\eta _\pm =\frac{1}{2}(1\pm \sqrt{14z})`$, $`a=m_\varphi ^2/m_{K^+}^2`$ and $`b=m_{\pi ^0\pi ^0}^2/m_{K^+}^2`$. The coupling constant $`g_s`$ comes from the strong decay amplitude $`𝒜(\varphi K^+K^{})=g_sϵ^{}(p_+p_{})`$ and takes the value $`|g_s|4.6`$ to agree with $`\mathrm{\Gamma }(\varphi K^+K^{})_{\mathrm{exp}}=2.19`$ MeV . This coupling is the part beyond standard ChPT which we have fixed phenomenologically.

The four-pseudoscalar amplitude is instead a standard ChPT amplitude which is found to depend linearly on the variable $`s=m_{\pi ^0\pi ^0}^2`$ only:

$$𝒜(K^+K^{}\pi ^0\pi ^0)_\chi =\frac{s}{4f_\pi f_K},$$
(7)
with $`f_K=1.22f_\pi `$. It is important to notice that in deriving this ChPT amplitude no use has been made of the relation $`s+t+u=2m_K^2+2m_\pi ^2`$ valid only for on-shell particles but not in our case with off-shell kaons in the chiral loop. Thanks to this genuine $`t`$ and $`u`$ independence, the chiral amplitude (7) factorizes in Eq. (4).

Integrating the invariant mass distribution for the $`\varphi \pi ^0\pi ^0\gamma `$ decay over the whole physical region one obtains $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )_\chi =219`$ eV and

$$B(\varphi \pi ^0\pi ^0\gamma )_\chi =0.49\times 10^4.$$
(8)
These results improve the prediction for this process given in Ref.  where $`SU(3)`$-breaking effects were ignored111 In the good $`SU(3)`$ limit, i.e. with the chiral amplitude $`𝒜(K^+K^{}\pi ^0\pi ^0)_\chi =s/4f_\pi f_Ks/4f_\pi ^2`$ and $`g_sg4.2`$ for the $`\rho \pi \pi `$ strong coupling constant used in Ref. , one obtains $`B(\varphi \pi ^0\pi ^0\gamma )_\chi =0.54\times 10^4`$..

## 3 Scalar meson exchange in $`\varphi \pi ^0\pi ^0\gamma `$

We now turn to the contributions coming from scalar resonance exchange. From a ChPT perspective their effects are encoded in the low energy constants of the higher order pieces of the ChPT Lagrangian. But the effects of the $`f_0(980)`$ —and its lower mass partner $`\sigma (500)`$ meson, if confirmed— should manifest in the $`\varphi \pi ^0\pi ^0\gamma `$ decays not as a constant term but rather through more complex resonant amplitudes. In this section we propose an amplitude which not only obeys the ChPT dictates in the lowest part of the $`\pi ^0\pi ^0`$ spectrum, as it must, but also generates the scalar meson effects for the higher part of the spectrum where the resonant poles should dominate.

The Linear Sigma Model (L$`\sigma `$M) , which we take as a first guidance probably requiring future refinements, will be shown to be particularly appropriate for our purposes. It is a well-defined $`U(3)\times U(3)`$ chiral model which incorporates ab initio both the nonet of pseudoscalar mesons together with its chiral partner, the scalar mesons nonet. In this context, the $`VP^0P^0\gamma `$ decays proceed through a loop of charged pseudoscalar mesons emitted by the initial vector. Because of the additional emission of a photon, charged pseudoscalar pairs with the initial $`J^{PC}=1^{}`$ quantum numbers can rescatter into $`J^{PC}=0^{++}`$ pairs of charged or neutral pseudoscalars. For the $`\varphi \pi ^0\pi ^0\gamma `$ decay the contributions from charged pion loops are again negligible compared to those from kaon loops and do not need to be considered for the same reasons discussed in the previous section. The elusive $`\sigma (500)`$ and the clearly established $`f_0(980)`$ scalar resonances are then expected, at least in principle, to play the central role in this $`K^+K^{}\pi ^0\pi ^0`$ rescattering process. Their contributions will be conveniently parametrized in terms of L$`\sigma `$M amplitudes compatible with ChPT for low dipion invariant masses.

The $`K^+K^{}\pi ^0\pi ^0`$ amplitude in the L$`\sigma `$M turns out to be

$$\begin{array}{ccc}\hfill 𝒜(K^+K^{}& & \pi ^0\pi ^0)_{\text{L}\sigma \text{M}}=g_{K^+K^{}\pi ^0\pi ^0}\hfill \\ & & \frac{g_{\sigma K^+K^{}}g_{\sigma \pi ^0\pi ^0}}{sm_\sigma ^2}\frac{g_{f_0K^+K^{}}g_{f_0\pi ^0\pi ^0}}{sm_{f_0}^2}\frac{g_{\kappa ^{}K^\pm \pi ^0}^2}{tm_\kappa ^2}\frac{g_{\kappa ^{}K^\pm \pi ^0}^2}{um_\kappa ^2},\hfill \end{array}$$
(9)
where the various coupling constants are fixed within the model and can be expressed in terms of $`f_\pi `$, $`f_K`$, the masses of the pseudoscalar and scalar mesons involved in the process, and the scalar meson mixing angle in the flavour basis $`\varphi _S`$ . In particular, the $`g_{K^+K^{}\pi ^0\pi ^0}`$ coupling accounting for the constant four-pseudoscalar amplitude can be expressed in a more convenient form by imposing that the amplitude $`𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}`$ vanishes in the soft-pion limit (either $`p0`$ or $`p^{}0`$). Then, the amplitude (9) can be rewritten as the sum of three terms each one depending only on $`s`$, $`t`$ and $`u`$:

$$\begin{array}{c}𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}𝒜_{\text{L}\sigma \text{M}}^s(s)+𝒜_{\text{L}\sigma \text{M}}^t(t)+𝒜_{\text{L}\sigma \text{M}}^u(u)\hfill \\ =\frac{sm_\pi ^2}{2f_\pi f_K}\left[\frac{m_K^2m_\sigma ^2}{D_\sigma (s)}\mathrm{c}\varphi _S(\mathrm{c}\varphi _S\sqrt{2}\mathrm{s}\varphi _S)+\frac{m_K^2m_{f_0}^2}{D_{f_0}(s)}\mathrm{s}\varphi _S(\mathrm{s}\varphi _S+\sqrt{2}\mathrm{c}\varphi _S)\right]\hfill \\ +\frac{tm_K^2}{4f_\pi f_K}\frac{m_\pi ^2m_\kappa ^2}{D_\kappa (t)}+\frac{um_K^2}{4f_\pi f_K}\frac{m_\pi ^2m_\kappa ^2}{D_\kappa (u)},\hfill \end{array}$$
(10)
where $`D_S(s)=sm_S^2+im_S\mathrm{\Gamma }_S`$ are the $`S=\sigma ,f_0`$ propagators —similar expressions hold for $`D_\kappa (t,u)`$— and $`(\mathrm{c}\varphi _S,\mathrm{s}\varphi _S)(\mathrm{cos}\varphi _S,\mathrm{sin}\varphi _S)`$. Notice that we start using simple Breit-Wigner expressions for the propagators in spite of the well known difficulties associated to threshold effects for the $`f_0`$. Moreover, the amplitude (10) is not strictly unitary. However, we will not discuss this issue in detail but refer the reader to Ref.  where a careful study of full unitarity in $`K\pi `$ scattering in the framework of the L$`\sigma `$M can be found. We should stress that the main concern of our analysis is to explicitly incorporate the scalar resonances and that, for this specific purpose, the tree level (and not fully unitary) amplitude (10) with simple expressions for the scalar propagators will be shown to be a reasonable starting point. A more sophisticated treatment of the $`f_0`$ propagator  accounting for the $`K\overline{K}`$ threshold effects ($`m_{f_0}2m_K`$) will be discussed later on.

A few remarks on the four-pseudoscalar amplitudes in Eqs. (9,10) and on their comparison with the ChPT amplitude in Eq. (7) are of interest:

* for $`m_S\mathrm{}`$ ($`S=\sigma ,f_0,\kappa `$), the L$`\sigma `$M amplitude (10) reduces to

$$\frac{sm_\pi ^2}{2f_\pi f_K}+\frac{t+u2m_K^2}{4f_\pi f_K},$$

which would coincide with the ChPT amplitude (7) if the on-shell condition $`s+t+u=2m_\pi ^2+2m_K^2`$ could be invoked. As shown in Eq. (9), $`𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}`$ consists of a constant four-pseudoscalar vertex plus three terms whose $`s,t,u`$ dependence is generated by the scalar propagators $`D_S(s,t,u)`$. Combining each one of these three terms with its corresponding part from the first, four-pseudoscalar vertex leads immediately to Eq. (10). For $`m_S\mathrm{}`$, this amplitude is linear in the $`s,t,u`$ variables, mimics perfectly the effects of the derivative and massive terms in the Lagrangian (3) and leads to the ChPT amplitude (7) once $`s+t+u=2m_K^2+2m_\pi ^2`$ is used. This corresponds to the aforementioned complementarity between ChPT and the L$`\sigma `$M, thus making the whole analysis quite reliable.
* however, the relation $`s+t+u=2m_K^2+2m_\pi ^2`$ cannot be invoked when plugging our L$`\sigma `$M amplitude (10) into the loop of virtual kaons leading to the $`\varphi \pi ^0\pi ^0\gamma `$ decay amplitude. This is in sharp contrast with the processes studied in Refs. , governed by $`t,u`$-independent L$`\sigma `$M amplitudes, and thus requires a new treatment for our case. To this aim, one has to distinguish between the $`𝒜_{\text{L}\sigma \text{M}}^s`$ term in Eq. (10) and the $`𝒜_{\text{L}\sigma \text{M}}^t+𝒜_{\text{L}\sigma \text{M}}^u`$ terms of the same amplitude. The former, being only $`s`$-dependent, can be directly introduced in Eq. (4) instead of the $`𝒜(K^+K^{}\pi ^0\pi ^0)_\chi `$ amplitude. Fortunately, this term contains most of the relevant dynamics due to the scalar resonances —the $`f_0`$(980) and, eventually, the $`\sigma `$(500) poles— for the process $`\varphi \pi ^0\pi ^0\gamma `$ where the dipion mass spectrum covers the range $`4m_\pi ^2m_{\pi \pi }^2m_\varphi ^2`$. The remaining two terms $`𝒜_{\text{L}\sigma \text{M}}^t+𝒜_{\text{L}\sigma \text{M}}^u`$, being $`t,u`$-dependent, cannot be simply plugged and evaluated into the loop integral222A more involved four-propagator loop integral would be required to compute the contributions of the $`t,u`$-channel in a precise way. but their contribution to the $`\varphi \pi ^0\pi ^0\gamma `$ amplitude can be identified and reasonably estimated. Indeed, the $`\kappa `$ contribution to the $`𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}`$ amplitude can be fixed by subtracting from the chiral-loop amplitude in Eq. (7) the contributions from the $`\sigma `$ and $`f_0`$ terms in Eq. (10) taking $`m_{\sigma ,f_0}\mathrm{}`$. The resulting expression corresponds to the desired $`\kappa `$ contribution in the $`m_\kappa \mathrm{}`$ limit. The difference between this $`\kappa `$ contribution for $`m_\kappa \mathrm{}`$ and the same contribution for a $`\kappa `$ of finite mass can be considered as negligible due to the high mass of the $`\kappa `$ and the lack of scalar poles in the $`t`$\- and $`u`$-channel of the $`\varphi \pi ^0\pi ^0\gamma `$ decay process333 The $`\kappa `$ has been identified in the literature as a broad resonance with a mass around 900 MeV (see for example Ref. ) or with the $`K_0^{}(1430)`$ scalar resonance as in Ref. . A model-independent analysis based on the analytic continuation of $`\pi K`$ scattering finds that there is room for a $`K_0^{}(1430)`$ scalar state but not for a $`\kappa (900)`$ one .. In other words, while the $`K^+K^{}\pi ^0\pi ^0`$ rescattering amplitude evaluated at lowest order in ChPT is necessarily a poor approximation to the (pole dominated) $`s`$-channel dynamics in $`\varphi \pi ^0\pi ^0\gamma `$ decays, the same lowest order amplitude can be taken as a reasonable estimate for the $`t`$\- and $`u`$-channel contributions where the $`\kappa `$ pole position cannot be kinematically approached. Therefore, an improved expression for the chiral-loop amplitude in Eq. (7) is the following:

$$\begin{array}{c}𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}=\frac{m_\pi ^2s/2}{2f_\pi f_K}+\frac{sm_\pi ^2}{2f_\pi f_K}\hfill \\ \times \left[\frac{m_K^2m_\sigma ^2}{D_\sigma (s)}\mathrm{c}\varphi _S(\mathrm{c}\varphi _S\sqrt{2}\mathrm{s}\varphi _S)+\frac{m_K^2m_{f_0}^2}{D_{f_0}(s)}\mathrm{s}\varphi _S(\mathrm{s}\varphi _S+\sqrt{2}\mathrm{c}\varphi _S)\right].\hfill \end{array}$$
(11)
As desired, it has no $`t,u`$ dependence and can be thus plugged in Eq. (4) instead of $`𝒜(K^+K^{}\pi ^0\pi ^0)_\chi `$. One then obtains the $`s`$-dependent amplitude

$$𝒜(\varphi \pi ^0\pi ^0\gamma )_{\text{L}\sigma \text{M}}=\frac{eg_s}{2\pi ^2m_{K^+}^2}\{a\}L(s)\times 𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}},$$
(12)
which will be used from now on. Notice that for large scalar masses one recovers Eq. (4), that the $`f_0`$(980) (and, eventually, $`\sigma `$(500)) $`s`$-channel poles do now appear, and that the remaining term, $`(m_\pi ^2s/2)/2f_\pi f_K`$, accounts for $`\kappa `$ exchange effects in the $`t,u`$-channel where, as stated before, no poles are approached and has thus been approximated by the corresponding part of the simple chiral-loop amplitude.
* the large widths of the scalar resonances break chiral symmetry if they are naively introduced in Eq. (9), an effect already noticed in Ref. . Accordingly, we introduce the $`\sigma (500)`$ and $`f_0(980)`$ widths in the propagators only after chiral cancellation of constant terms in the amplitude. In this way the pseudo-Goldstone nature of pions is preserved.
* the $`\pi ^0\pi ^0`$ invariant mass spectrum for the $`\varphi \pi ^0\pi ^0\gamma `$ decay covers the region where the presence of the $`f_0(980)`$ (and the $`\sigma (500)`$) meson(s) should manifest. Because of the presence of the corresponding propagators in Eq. (11) —closely linked to the ChPT amplitude and thus expected to account for the lowest part of the $`\pi ^0\pi ^0`$ spectrum— one should also be able to reproduce the effects of the $`f_0`$ (and the $`\sigma (500)`$) pole(s) at higher $`\pi ^0\pi ^0`$ invariant mass values.

## 4 Vector meson exchange in $`\varphi \pi ^0\pi ^0\gamma `$

In addition to the L$`\sigma `$M contributions, which can be viewed as an improved version of the chiral-loop predictions, the analysis should be extended to include vector meson exchange in the $`t`$\- and $`u`$-channel. These VMD contributions were already considered in Ref. . In this framework, $`\varphi \pi ^0\pi ^0\gamma `$ proceeds through $`\omega `$-$`\varphi `$ mixing followed by the exchange of intermediate $`\rho ^0`$ mesons in the direct or crossed channel of the $`\varphi \pi ^0\rho ^0\pi ^0\pi ^0\gamma `$ decay chain.

In order to describe these vector meson contributions we use the $`SU(3)`$-symmetric Lagrangians

$$\begin{array}{ccc}\hfill _{\mathrm{VVP}}& =& \frac{G}{\sqrt{2}}ϵ^{\mu \nu \alpha \beta }_\mu V_\nu _\alpha V_\beta P,\hfill \\ \hfill _{\mathrm{V}\gamma }& =& 4f^2egA_\mu QV^\mu ,\hfill \end{array}$$
(13)
where $`G=\frac{3g^2}{4\pi ^2f}`$ is the $`\omega \rho \pi `$ coupling constant and $`|g|4.0`$ as follows from various $`\rho `$ and $`\omega `$ decay data . The VMD amplitude for $`\varphi (q^{},ϵ^{})\pi ^0(p)\pi ^0(p^{})\gamma (q,ϵ)`$ is then found to be

$$𝒜(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}}=\frac{ϵ}{3}\frac{G^2e}{\sqrt{2}g}\left(\frac{P^2\{a\}+\{b(P)\}}{M_\rho ^2P^2iM_\rho \mathrm{\Gamma }_\rho }+\frac{P_{}^{}{}_{}{}^{2}\{a\}+\{b(P^{})\}}{M_\rho ^2P_{}^{}{}_{}{}^{2}iM_\rho \mathrm{\Gamma }_\rho }\right),$$
(14)
with $`\{a\}`$ as in Eq. (4) and a new amplitude

$$\begin{array}{ccc}\hfill \{b(P)\}& =& (ϵ^{}ϵ)(q^{}P)(qP)(ϵ^{}P)(ϵP)(q^{}q)\hfill \\ & & +(ϵ^{}q)(ϵP)(q^{}P)+(ϵq^{})(ϵ^{}P)(qP),\hfill \end{array}$$
(15)
where $`P=p+q`$ and $`P^{}=p^{}+q`$ are the momenta of the intermediate $`\rho ^0`$ meson in the $`t`$\- and $`u`$-channel, respectively. From this VMD amplitude one easily obtains $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}}=37`$ eV and

$$B(\varphi \pi ^0\pi ^0\gamma )_{\mathrm{VMD}}=8.3\times 10^6,$$
(16)
in agreement with the results in Ref.  once the same numerical inputs are used. Further details on these contributions are given in the following section.

## 5 Final results

Our final results for $`𝒜(\varphi \pi ^0\pi ^0\gamma )`$ are thus the sum of the VMD contribution in Eq. (14) plus the L$`\sigma `$M contribution containing the scalar resonance effects in Eq. (12). The final $`\pi ^0\pi ^0`$ invariant mass distribution for the process is444 In terms of the photon energy, $`E_\gamma =(m_\varphi ^2m_{\pi ^0\pi ^0}^2)/(2m_\varphi )`$, the photonic spectrum is written as $`d\mathrm{\Gamma }/dE_\gamma =(m_\varphi /m_{\pi ^0\pi ^0})\times d\mathrm{\Gamma }/dm_{\pi ^0\pi ^0}`$.

$$\frac{d\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )}{dm_{\pi ^0\pi ^0}}=\frac{d\mathrm{\Gamma }_{\text{L}\sigma \text{M}}}{dm_{\pi ^0\pi ^0}}+\frac{d\mathrm{\Gamma }_{\text{VMD}}}{dm_{\pi ^0\pi ^0}}+\frac{d\mathrm{\Gamma }_{\text{int}}}{dm_{\pi ^0\pi ^0}},$$
(17)
where the L$`\sigma `$M term (the scalar or *signal* contribution) is

$$\begin{array}{cc}\hfill \frac{d\mathrm{\Gamma }_{\text{L}\sigma \text{M}}}{dm_{\pi ^0\pi ^0}}& =\frac{1}{2}\frac{\alpha }{192\pi ^5}\frac{g_s^2}{4\pi }\frac{m_\varphi ^4}{m_{K^+}^4}\frac{m_{\pi ^0\pi ^0}}{m_\varphi }\left(1\frac{m_{\pi ^0\pi ^0}^2}{m_\varphi ^2}\right)^3\sqrt{1\frac{4m_{\pi ^0}^2}{m_{\pi ^0\pi ^0}^2}}\hfill \\ & \times |L(m_{\pi ^0\pi ^0}^2)|^2|𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}|^2,\hfill \end{array}$$
(18)
with the four-pseudoscalar amplitude taken from Eq. (11). The VMD term (the vector or *background* contribution) and the contribution resulting from the interference of both amplitudes are given by

$$\begin{array}{cc}\hfill \frac{d\mathrm{\Gamma }_{\text{[VMD,int]}}}{dm_{\pi ^0\pi ^0}}& =\frac{1}{2}\frac{1}{256\pi ^3}\frac{m_{\pi ^0\pi ^0}}{m_\varphi }\left(1\frac{m_{\pi ^0\pi ^0}^2}{m_\varphi ^2}\right)\sqrt{1\frac{4m_{\pi ^0}^2}{m_{\pi ^0\pi ^0}^2}}\hfill \\ & \times _1^1dxA_{[\mathrm{VMD},\mathrm{int}]}(m_{\pi ^0\pi ^0},x),\hfill \end{array}$$
(19)
where one explicitly has

$$\begin{array}{c}A_{\mathrm{VMD}}(m_{\pi ^0\pi ^0},x)\frac{1}{3}_{\mathrm{pol}}|𝒜_{\mathrm{VMD}}|^2=\frac{1}{3}\left(\frac{ϵ}{3}\frac{G^2e}{\sqrt{2}g}\right)^2\hfill \\ \times [\frac{1}{8}(m_{\pi ^0}^82m_\varphi ^2m_{\pi ^0}^6+4\stackrel{~}{m}_\rho ^4m_{\pi ^0}^4+m_\varphi ^4m_{\pi ^0}^42\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^2m_{\pi ^0}^44\stackrel{~}{m}_\rho ^6m_{\pi ^0}^2\hfill \\ 2\stackrel{~}{m}_\rho ^2m_\varphi ^4m_{\pi ^0}^2+4\stackrel{~}{m}_\rho ^4m_\varphi ^2m_{\pi ^0}^24\stackrel{~}{m}_\rho ^4\stackrel{~}{m}_\rho ^2m_{\pi ^0}^2+2\stackrel{~}{m}_\rho ^2m_\varphi ^2\stackrel{~}{m}_\rho ^2m_{\pi ^0}^2+\stackrel{~}{m}_\rho ^8\hfill \\ +\stackrel{~}{m}_\rho ^4m_\varphi ^4+2\stackrel{~}{m}_\rho ^4\stackrel{~}{m}_\rho ^42\stackrel{~}{m}_\rho ^6m_\varphi ^2+2\stackrel{~}{m}_\rho ^6\stackrel{~}{m}_\rho ^22\stackrel{~}{m}_\rho ^4m_\varphi ^2\stackrel{~}{m}_\rho ^2)\times \frac{1}{|D_\rho (\stackrel{~}{m}_\rho ^2)|^2}\hfill \\ +\frac{1}{16}(m_{\pi ^0}^8\stackrel{~}{m}_\rho ^2m_\varphi ^2m_{\pi ^0}^4+2\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^2m_{\pi ^0}^4m_\varphi ^2\stackrel{~}{m}_\rho ^2m_{\pi ^0}^44\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^4m_{\pi ^0}^2\hfill \\ 4\stackrel{~}{m}_\rho ^4\stackrel{~}{m}_\rho ^2m_{\pi ^0}^2+4\stackrel{~}{m}_\rho ^2m_\varphi ^2\stackrel{~}{m}_\rho ^2m_{\pi ^0}^2+\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^6+3\stackrel{~}{m}_\rho ^4\stackrel{~}{m}_\rho ^4\stackrel{~}{m}_\rho ^2m_\varphi ^2\stackrel{~}{m}_\rho ^4\hfill \\ +\stackrel{~}{m}_\rho ^6\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^4m_\varphi ^2\stackrel{~}{m}_\rho ^2)\times 2\text{Re}\left(\frac{1}{D_\rho (\stackrel{~}{m}_\rho ^2)D_\rho ^{}(\stackrel{~}{m}_\rho ^2)}\right)+(\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^2)],\hfill \end{array}$$
(20)
and

$$\begin{array}{c}A_{\mathrm{int}}(m_{\pi ^0\pi ^0},x)\frac{2}{3}\text{Re}_{\mathrm{pol}}𝒜_{\text{L}\sigma \text{M}}𝒜_{\mathrm{VMD}}^{}=\frac{1}{3}\left(\frac{eg_s}{2\pi ^2m_{K^+}^2}\right)\left(\frac{ϵ}{3}\frac{G^2e}{\sqrt{2}g}\right)\hfill \\ \times \mathrm{\hspace{0.17em}2}\text{Re}\{L(m_{\pi ^0\pi ^0}^2)𝒜(K^+K^{}\pi ^0\pi ^0)_{\text{L}\sigma \text{M}}\hfill \\ \times \frac{1}{4}[\frac{\stackrel{~}{m}_\rho ^2(\stackrel{~}{m}_\rho ^2+\stackrel{~}{m}_\rho ^22m_{\pi ^0}^2)^2m_\varphi ^2(\stackrel{~}{m}_\rho ^2m_{\pi ^0}^2)^2}{D_\rho ^{}(\stackrel{~}{m}_\rho ^2)}+(\stackrel{~}{m}_\rho ^2\stackrel{~}{m}_\rho ^2)]\},\hfill \end{array}$$
(21)
with

$$\begin{array}{ccc}\hfill P^2\stackrel{~}{m}_\rho ^2& =& m_{\pi ^0}^2+\frac{m_\varphi ^2m_{\pi ^0\pi ^0}^2}{2}\left(1x\sqrt{1\frac{4m_{\pi ^0}^2}{m_{\pi ^0\pi ^0}^2}}\right),\hfill \\ \hfill P_{}^{}{}_{}{}^{2}\stackrel{~}{m}_\rho ^2& =& m_{\pi ^0}^2+\frac{m_\varphi ^2m_{\pi ^0\pi ^0}^2}{2}\left(1+x\sqrt{1\frac{4m_{\pi ^0}^2}{m_{\pi ^0\pi ^0}^2}}\right)\hfill \\ & =& 2m_{\pi ^0}^2+m_\varphi ^2m_{\pi ^0\pi ^0}^2\stackrel{~}{m}_\rho ^2.\hfill \end{array}$$
(22)
The separate contributions from L$`\sigma `$M, VMD and their interference, as well as the total result are shown in Fig. 1. In this figure, whose main purpose is to illustrate the possibilities of our approach, we use simple Breit-Wigner formulae for the $`\sigma `$ and the $`f_0`$ propagators. We take $`m_\sigma =478`$ MeV and $`\mathrm{\Gamma }_\sigma =324`$ MeV, which are the central values of the measurements quoted in , $`m_{f_0}=980`$ MeV and $`\mathrm{\Gamma }_{f_0}=70`$ MeV, from Ref. , and $`\varphi _S=6^{}`$ for the scalar mixing angle. The global agreement with the data is rather good although the model will require further improvements. As expected, $`f_0(980)`$ scalar meson exchange contributes decisively to achieve this agreement. Indeed, for the integrated decay width one now obtains $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )_{\text{L}\sigma \text{M+VMD}}=530`$ eV and for the branching ratio

$$B(\varphi \pi ^0\pi ^0\gamma )_{\text{L}\sigma \text{M+VMD}}=1.19\times 10^4,$$
(23)
quite in line with the experimental results quoted in Eqs. (1) and (2). Without the VMD contribution, our result decreases only by some 10%, $`B(\varphi \pi ^0\pi ^0\gamma )_{\text{L}\sigma \text{M}}=1.07\times 10^4`$, but it still remains well above the value quoted in Eq. (8), $`B(\varphi \pi ^0\pi ^0\gamma )_\chi =0.49\times 10^4`$, which did not contain the crucial scalar contributions.

In order to show the sensitivity of our treatment on the $`f_0`$ decay width and on the scalar mixing angle we have plotted in Fig. 2 our predictions for the $`\pi ^0\pi ^0`$ invariant mass distribution in $`\varphi \pi ^0\pi ^0\gamma `$ for various values of $`\mathrm{\Gamma }_{f_0}`$ and $`\varphi _S`$. The prediction including chiral-loops (without scalar poles) and VMD corrections is also shown for comparison. The invariant mass spectrum and the branching ratio are very sensitive to $`\mathrm{\Gamma }_{f_0}`$ and, even more, to $`\varphi _S`$. Present data on $`\varphi \pi ^0\pi ^0\gamma `$  favour $`\varphi _S6^{}`$ —which is in the expected range $`14^{}\varphi _S3^{}`$ — and $`\mathrm{\Gamma }_{f_0}70`$ MeV for the total width of the $`f_0`$ Breit-Wigner. This latter value is compatible with the PDG estimates, $`\mathrm{\Gamma }_{f_0}40`$$`100`$ MeV, and somewhat above other recent values obtained with usual Breit-Wigner expressions: $`\mathrm{\Gamma }_{f_0}=56\pm 20`$ MeV  or $`44\pm 2\pm 2`$ MeV . Notice that in our approach $`\varphi _S6^{}`$ implies $`\mathrm{\Gamma }(f_0\pi \pi )35`$ MeV accounting for some 50% of the total $`f_0`$-width. Other decay channels, such as $`f_0K\overline{K}`$ for which our treatment predicts a large coupling, should account for the remaining 50%. But this requires a two-channel analysis around the $`f_0K\overline{K}`$ threshold to which we now turn.

The original two-channel treatment of scalar resonances near the $`K\overline{K}`$ threshold due to Flatté  has been improved in Refs.  and is now widely applied in experimental analyses of $`VPP\gamma `$ decays . It amounts to the substitution of the simple expression $`D_{f_0}(s)=sm_{f_0}^2+im_{f_0}\mathrm{\Gamma }_{f_0}`$ by a complete one-loop $`f_0`$ propagator $`D_{f_0}(s)=sm_{f_0}^2\text{Re}\mathrm{\Pi }(m_{f_0}^2)+\mathrm{\Pi }(s)`$, where the term $`\mathrm{\Pi }(s)\text{Re}\mathrm{\Pi }(m_{f_0}^2)`$ takes into account the finite-width corrections due to the $`\pi \pi `$ and $`K\overline{K}`$ channels. Within our L$`\sigma `$M approach, the $`f_0\pi \pi `$ and $`f_0K\overline{K}`$ partial widths can be expressed in terms of $`m_{f_0}`$, $`\varphi _S`$ and accurately measured pseudoscalar masses and decay constants. For $`m_{f_0}=985`$ MeV and $`\varphi _S=9^{}`$ we obtain the various curves shown in Fig. 3 which explain quite reasonably the available data. Integrating over the whole kinematical region one obtains $`\mathrm{\Gamma }(\varphi \pi ^0\pi ^0\gamma )=518`$ eV and

$$B(\varphi \pi ^0\pi ^0\gamma )=1.16\times 10^4.$$
(24)
The shape of the $`\pi \pi `$ mass spectrum depends strongly on the value for $`\varphi _S`$, as shown in Fig. 4, and follows the same pattern previously shown in Fig. 2.

## 6 Comments and conclusions

The recent publication by experimental groups working at VEPP-2M  and DA$`\mathrm{\Phi }`$NE  of the first data on the $`\varphi \pi ^0\pi ^0\gamma `$ decay opens the possibility to improve our understanding of the controversial $`I=0`$, scalar resonances below and around the 1 GeV mass region. In our present approach we deal with these scalar resonances by combining ChPT requirements with L$`\sigma `$M dynamics. This treatment was first applied with some success to $`\varphi \pi ^0\eta \gamma `$  and later to $`\rho ,\omega \pi ^0\pi ^0\gamma `$ . The main feature in our approach is the use of an amplitude which agrees with ChPT in the low part of the two-pseudoscalar invariant-mass spectrum but which also includes the scalar poles dominating at higher invariant-mass values. According to our final amplitude (12), the $`\varphi \pi ^0\pi ^0\gamma `$ decay proceeds through a loop of charged kaons (as discussed, among others, in Refs. ) strongly coupled to the initial $`\varphi `$ and to the final dipion system through scalar resonance formation in the $`s`$-channel. Additional vector meson contributions, acting as a kind of background to the previous and more interesting scalar exchange, have been recalculated in the framework of VMD thus confirming older results .

The measured values for the branching ratio $`B(\varphi \pi ^0\pi ^0\gamma )`$ in Eqs. (1,2) and the shape of the $`\pi \pi `$ mass spectrum  cannot be explained without a sizable contribution from the $`f_0`$, which depends on the expression of the $`f_0`$ propagator. If a simple Breit-Wigner is used, a good description of the data is achieved for $`\mathrm{\Gamma }_{f_0}70`$ MeV and $`\varphi _S6^{}`$. Although these two values are only marginally consistent within our model, the whole situation looks rather satisfactory if one takes into account the simplicity of our approach and the well known difficulties one has to face in scalar meson physics. If a complete one-loop expression is used for the $`f_0`$ propagator , a fully consistent description of the data is achieved with $`\varphi _S9^{}`$. The range of values for $`\varphi _S`$ is thus considerably restricted regardless the expression adopted for the $`f_0`$ propagator.

In turn, such small values for $`|\varphi _S|`$ imply a large $`\sigma \pi \pi `$ coupling which, in principle, should clearly manifest in the low part of the dipion mass spectrum. But this is not the case in a L$`\sigma `$M approach like ours where the $`\sigma (500)`$ contribution contains a $`\sigma K\overline{K}`$ coupling which is proportional to $`m_\sigma ^2m_K^2`$ and thus almost vanishing. This depletion effect of the $`\sigma (500)`$ is somehow visible in the DA$`\mathrm{\Phi }`$NE data which show very small values for $`dB(\varphi \pi ^0\pi ^0\gamma )/dm_{\pi \pi }`$ in the $`m_{\pi \pi }m_\sigma 500`$ MeV region. Indeed, the chiral-loop amplitude (without scalar poles) plus the VMD background are known to predict values  which are two or three times larger than the data in this range of masses. This feature, namely, the prediction of small $`\sigma K\overline{K}`$ and large $`\sigma \pi \pi `$ couplings (because of their proportionality to $`m_\sigma ^2m_K^2`$ and to $`m_\sigma ^2m_\pi ^2`$, respectively) is specific of our L$`\sigma `$M approach. It admits a crucial and simple test. Contrary to what happens for $`m_{\pi \pi }`$ around 500 MeV in $`\varphi \pi ^0\pi ^0\gamma `$ decays, data on $`\rho \pi ^0\pi ^0\gamma `$ —dominated in this case by pion rather than kaon loops— should show a non-suppressed $`\sigma (500)`$ contribution around $`m_{\pi \pi }500`$ MeV .

In summary, the recently published data on $`\varphi \pi ^0\pi ^0\gamma `$ decays can be reasonably described in a rather simple context based on the complementary between ChPT and the L$`\sigma `$M. Future refinements of our approach and comparison with forthcoming and more accurate data should contribute considerably to clarify one of the most challenging aspects of present day hadron physics, namely, the structure of the lowest lying scalar states.

## Acknowledgements

Work partly supported by the EEC, TMR-CT98-0169, EURODAPHNE network. M. N. was supported by CONCYTEG-Mexico under project 00-16-CONCYTEG/CONACYT-075. J. L. L. M. acknowledges partial financial support from CONACyT and CONCyTEG. We thank N. N. Achasov for useful comments.


format: text

text:
# Eddington limit and radiative transfer in highly inhomogeneous atmospheres

## 1 Introduction

Radiation-dominated atmospheres of accretion disks and massive stars permeated by a moderately strong magnetic field are susceptible to the “photon bubble” instability (Arons 1992; Gammie 1998; Begelman 2001). In this process, high density regions tend to be pulled downward along the field lines and low density regions are pushed upward by radiation forces. The gas elements accelerated by radiation forces enter density maxima, where radiation forces decrease, and then progress downward again completing the cycle. Because subsequent acceleration episodes are increasingly large and magnetic tension prevents high density regions from spreading sideways, the density contrast increases. This leads to density inhomogeneities on scales shorter than the characteristic scale height of the accretion disk or stellar atmosphere. Under such conditions, radiation tends to bypass high density regions and travel more freely through tenuous ones. If the low- and high-density regions are dynamically coupled and most of the mass is in the high density phase, then the flux necessary to support the atmosphere against gravity can exceed Eddington limit. Photon bubble instability may be applicable to objects suspected of having super-Eddington luminosities, such as “ultraluminous X-ray sources” (Begelman 2002) and narrow-line Seyfert 1 galaxies (King and Puchnarewicz 2002). It has also been argued that the instability may occur in simple non-magnetized Thomson atmospheres as they approach the Eddington limit (Shaviv 2001b). The origin of this instability and the statistical properties of the inhomogeneities, such as their size relative to the size of the atmosphere, are different than that predicted by the magnetic photon bubble instability. This mechanism has been invoked to explain the discrepancy between relatively high outburst luminosities and relatively low outflow velocities in novae (Shaviv 1998, 2000, 2001a).
Radiative transfer in a medium where parts of the gas are optically thin, such as the upper regions of atmospheres in particular, cannot be treated using the standard diffusion approximation. The application of this approximation in the optically thin regime may lead to the radiation propagation rate exceeding the free-streaming rate $`|𝐅|=cu`$, where $`𝐅`$ and $`u`$ are the radiation flux and energy density, respectively. Previous work on radiative transfer under such conditions focused on the development of flux-limited diffusion approximations (e.g., Levermore and Pomraning (1981); Melia and Zylstra (1991); Anile and Romano (1992)). A serious limitation of these methods is that they are applicable only in cases where the angular distribution of the specific intensity is a slowly varying function of space and time.
In this paper, we focus on radiative transfer deep within atmospheres (i.e., at high optical depth), where the effects of flux-limited diffusion are less important. But we consider the case where large density gradients – on scales of order the photon mean free path – lead to rapid fluctuations in the angular distribution function. It is easy to see why the standard diffusion approximations (with or without flux-limiting) fail under these circumstances. The standard diffusion approach predicts that the flux responds instantaneously to local changes in density. But in reality the flux can only respond to inhomogeneities provided that the density changes on scales much larger than the photon mean free path. When density inhomogeneities are optically thin, radiation does not “see” any density contrasts. Using a simple ansatz, we have derived a modified diffusion equation and show that the effects of “photon viscosity” (i.e., off-diagonal elements of the radiation stress tensor) play an important role in radiative transfer in this case. Using analytical and Monte Carlo techniques we find that our analytical approach is much more accurate than standard diffusion or multi-stream approximations. Our results can be applied to the study of dynamical coupling between radiation forces and density inhomogeneities in radiation-dominated accretion disks.
The structure of this paper is as follows. In section 2 we discuss three-dimensional radiative transfer and show that the “photon viscosity” terms are important. In section 3 we constrain ourselves to the case of periodic, planar density waves embedded deep within an atmosphere and derive the corresponding Eddington enhancement factors under the assumption of global dynamical equilibrium. In section 4 we compare analytical formulas with Monte Carlo calculations and discuss our results. We propose that our method could be incorporated into existing radiation hydrodynamics codes such as the RHD module for ZEUS developed by Turner and Stone (2001).

## 2 Radiative transfer in an inhomogeneous atmosphere

We consider radiative transfer within infinite highly inhomogeneous atmospheres where the effects of flux-limited diffusion are very much reduced. Moreover, in a porous atmosphere, the radiation diffusion timescale is most likely much shorter than the characteristic gas dynamical timescale, i.e., radiation is probably not trapped by the motion of gaseous inhomogeneities (Begelman 2001). Therefore, for the purpose of calculating radiative transfer, we can neglect the time-dependence and motion of the gas distribution. We approximate the intensity distribution as

$$I(𝐱,\widehat{\mathrm{\Omega }})=I_0(𝐱,\widehat{\mathrm{\Omega }})+\frac{3}{4\pi }\widehat{\mathrm{\Omega }}𝐅(𝐱),$$
(1)
where $`𝐅`$ is the local flux vector, $`\widehat{\mathrm{\Omega }}`$ is the directional unit vector, and $`𝐱`$ is the position vector. This approximation means that the first nontrivial correction to intensity $`I_0(𝐱,\widehat{\mathrm{\Omega }})`$ is symmetric with respect to the direction defined by the local flux. This has to be contrasted with a familiar case of a plane-parallel atmosphere. In such a case, it is customary to assume that the directional distribution of intensity $`I(𝐱,\widehat{\mathrm{\Omega }})`$ is symmetric relative to the vertical direction which coincides with the direction of the flux vector. However, in the case of highly inhomogeneous atmosphere, radiation bypasses denser regions and the flux vector can rapidly change its orientation and only the volume-averaged flux is vertical. Thus, in our approach, the local “symmetry axis” is changing direction throughout the atmosphere and no global symmetry is required. We also assume that all the odd moments of $`I_0`$ vanish, i.e., $`\widehat{\mathrm{\Omega }}_iI_0𝑑\mathrm{\Omega }=\widehat{\mathrm{\Omega }}_i\widehat{\mathrm{\Omega }}_j\widehat{\mathrm{\Omega }}_kI_0𝑑\mathrm{\Omega }=0`$, but otherwise make no assumptions about the directional dependence of $`I_0`$. The equation of radiative transfer for a 3D scattering atmosphere reads

$$\frac{1}{\sigma }\widehat{\mathrm{\Omega }}I=I+\frac{1}{4\pi }I𝑑\mathrm{\Omega },$$
(2)
where $`\sigma =\rho \kappa `$ is the scattering coefficient. Taking zeroth and first moments of eq. (2) we obtain

$$𝐅=0$$
(3)
and

$$𝐅=\frac{c}{\sigma }𝐓,$$
(4)
where $`𝐓`$ is the radiation stress tensor, the components of which are

$$T_{ij}=\frac{1}{c}\widehat{\mathrm{\Omega }}_i\widehat{\mathrm{\Omega }}_jI_0𝑑\mathrm{\Omega }.$$
(5)
The closure relation for $`T_{ij}`$ in terms of $`F_i`$ and $`J\frac{1}{4\pi }I𝑑\mathrm{\Omega }`$ can be obtained by calculating the second moment of eq. (2), assuming the form of the intensity given by eq. (1). This leads to the equation for the radiation stress tensor:

$$T_{ij}=\{\begin{array}{cc}\frac{u}{3}\frac{1}{5\sigma c}\left(\frac{F_i}{x^j}+\frac{F_j}{x^i}\right)\hfill & \text{for}i=j\hfill \\ \frac{1}{10\sigma c}\left(\frac{F_i}{x^j}+\frac{F_j}{x^i}\right)\hfill & \text{for}ij\hfill \end{array}$$
(6)
where $`u=4\pi J/c`$ is the energy density and $`J`$ is the mean intensity. Note that (i) the diagonal terms $`T_{ii}`$ may be different from one another, and (ii) the off-diagonal elements of the radiation stress tensor do not vanish. The first point implies that this approach incorporates variable Eddington factors. The off-diagonal elements are responsible for “photon viscosity” and are non-zero even though bulk gas motions in the atmosphere were assumed to be negligible. This is because the “photon fluid” is moving through the gas, and may exert shear stresses. Substituting eq. (6) into eq. (4) and using eq. (3), we obtain:

$`F_i`$ $`=`$ $`{\displaystyle \frac{c}{3\sigma }}{\displaystyle \frac{u}{x^i}}+{\displaystyle \frac{1}{10\sigma ^2}}[{\displaystyle \frac{^2F_i}{x_jx^j}}+2{\displaystyle \frac{^2F_i}{x_i^2}}`$ (7)
$`{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{\sigma }{x_j}}({\displaystyle \frac{F_i}{x^j}}+{\displaystyle \frac{F_j}{x^i}}){\displaystyle \frac{2}{\sigma }}{\displaystyle \frac{\sigma }{x^i}}{\displaystyle \frac{F_i}{x_i}}].`$
In the above equation, summation is imposed only over repeated $`j`$ indices. Equation (7) and equation (3) are the governing equations of radiative transfer in our approximation.

## 3 2D density inhomogeneities

In order to illustrate our method, we now consider the simplified case of a plane-parallel wave density pattern (see Fig. 1). The scattering coefficient is assumed to be only a function of the distance $`\xi `$ perpendicular to the slabs, i.e., $`\sigma =\rho \kappa =\sigma (\xi )`$ with $`\xi =\mu x+(1\mu ^2)^{1/2}z`$, where $`\mu =\mathrm{cos}\psi `$. We also assume that $`/y=0`$ and that the components of flux depend only on $`\xi `$. We consider an atmosphere which is in global (i.e., volume-averaged) hydrostatic equilibrium, where radiation pressure balances gravity $`g\widehat{z}`$. This implies that the gradient of energy density must be of the form

$$u=u^{}(\xi )\xi \frac{3\sigma }{c}F_{\mathrm{Edd}}\widehat{z},$$
(8)
where a prime denotes differentiation with respect to $`\xi `$, $`\sigma `$ is the volume average of the scattering coefficient, and $`F_{\mathrm{Edd}}=gc/\kappa `$ is the Eddington flux. From equation (7) and equation (3) we obtain:

$$\mu F_z(1\mu ^2)^{1/2}F_x=\mu \frac{\sigma }{\sigma }F_{\mathrm{Edd}}+\frac{1}{\alpha ^2\sigma \mu }\left(\frac{F_z^{}}{\sigma }\right)^{},$$
(9)
where $`\alpha ^2(\mu )=10/(1+4\mu ^24\mu ^4)`$. Integrating eq. (3) over $`\xi `$ we obtain

$$\mu F_x+(1\mu ^2)^{1/2}F_z=const(1\mu ^2)^{1/2}F_0,$$
(10)
where $`F_0`$ is the integration constant. Using eq. (10) to eliminate $`F_x`$ from eq. (9) we get

$$F_z=(1\mu ^2)F_0+\mu ^2\frac{\sigma }{\sigma }F_{\mathrm{Edd}}+\frac{1}{\alpha ^2\sigma }\left(\frac{F_z^{}}{\sigma }\right)^{}.$$
(11)
Multiplying equation (11) by $`\sigma `$ and then volume-averaging it and demanding that the solution be bounded, we have

$$\sigma F_{\mathrm{Edd}}=\sigma F_z=(1\mu ^2)\sigma F_0+\mu ^2\sigma F_{\mathrm{Edd}},$$
(12)
where the first equality in equation (12) comes from the requirement that the atmosphere be in global hydrostatic equilibrium. Therefore, $`F_0=F_{\mathrm{Edd}}`$ and the final equation for the vertical flux is

$$F_z=\left(1\mu ^2+\mu ^2\frac{\sigma }{\sigma }\right)F_{\mathrm{Edd}}+\frac{1}{\alpha ^2\sigma }\left(\frac{F_z^{}}{\sigma }\right)^{}.$$
(13)
### 3.1 Eddington enhancement factor

Using equation (13) we can now calculate the Eddington enhancement factor $`lF_z/F_{\mathrm{Edd}}`$. We can simplify and non-dimensionalize equation (13) by defining

$$f\frac{F_z}{F_{\mathrm{Edd}}}(1\mu ^2).$$
(14)
Thus, the Eddington factor is

$$l(1\mu ^2)+f,$$
(15)
where $`f`$ is the volume average of $`f`$. Defining optical depth $`d\tau =\sigma d\xi `$, letting the prime now denote differentiation with respect to $`\tau `$, and normalizing the scattering coefficient to its mean value (i.e. $`\sigma \sigma /\sigma `$), we obtain a simplified form of equation (13)

$$f^{\prime \prime }\alpha ^2f=\frac{\alpha ^2\mu ^2}{\sigma }.$$
(16)
We now consider a periodic slab model in which the scattering coefficient is given by:

$$\sigma =\{\begin{array}{cc}\sigma _1\hfill & \text{if }\tau _1<\tau <0\text{, region 1}\hfill \\ \sigma _2\hfill & \text{if }\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<\tau <\tau _2\text{, region 2 .}\hfill \end{array}$$
(17)
The solutions to equation (16) in regions 1 and 2 are

$$f_{1,2}(\tau )=a_{1,2}\mathrm{sinh}(\alpha \tau )+b_{1,2}\mathrm{cosh}(\alpha \tau )+\frac{\mu ^2}{\sigma _{1,2}},$$
(18)
where $`a_{1,2}`$ and $`b_{1,2}`$ are the integration constants. Functions $`f`$ and $`f^{}`$ have to be continuous across each slab boundary. Therefore the matching conditions are:

$`f_1(0)=f_2(0),`$ $`f_1^{}(0)=f_2^{}(0)`$ (19)
$`f_1(\tau _1)=f_2(\tau _2),`$ $`f_1^{}(\tau _1)=f_2^{}(\tau _2).`$ (20)
Using the above matching conditions to derive the integration constants $`a_{1,2}`$ and $`b_{1,2}`$ and then volume averaging the solution for $`f`$, we obtain the expression for $`f`$

$$f=\frac{\mu ^2\sigma }{\xi _1+\xi _2}\left[\frac{\xi _1}{\sigma _1}+\frac{\xi _2}{\sigma _2}\frac{2}{\alpha }\frac{(\mathrm{\Delta }\sigma )^2}{\sigma _1^2\sigma _2^2}\frac{\mathrm{sinh}x\mathrm{sinh}y}{\mathrm{sinh}(x+y)}\right],$$
(21)
where $`x=\alpha \tau _1/2`$, $`y=\alpha \tau _2/2`$, $`\mathrm{\Delta }\sigma =\sigma _2\sigma _1`$, and $`\xi _{1,2}=\tau _{1,2}/\sigma _{1,2}`$. The Eddington enhancement factor $`l`$ follows directly from equations (21) and (15).

## 4 Comparison with Monte Carlo simulations

We now compare our analytical results with Monte Carlo simulations. The setup of the numerical experiment was as follows. We instantaneously injected a large number of photons in the equatorial plane (i.e., $`z=0`$) of a very flat, three-dimensional box (i.e., $`z[z_o,z_o]`$; $`x,y[w_o,w_o]`$, where $`z_ow_o`$). Densities and height of the computational box were chosen in such a way as to assure that the optical thickness in the vertical direction would always be large throughout the box. The initial photon angular distribution was uniform. Although we included the anisotropy due to the Thomson scattering cross section, we found this effect to have no influence on our final results. We followed the trajectories of all photons and calculated photon travel times between scatterings. Momentum transfer for every scattering was calculated using the method of weights (Pozdnyakov at al. 1983). We then computed the force exerted on the atmosphere as a function of the time delay following the instantaneous photon injection. Of course, as photons diffuse out of the atmosphere, the force exerted on the gas gradually declines. Therefore, the total force, corresponding to a continuous photon flux, was calculated by superposing many such time-dependent force distributions due to groups of photons injected (instantaneously) at uniform time intervals. The total force exerted on the atmosphere was characterized by a gradual increase with time followed by a flat maximum. The Eddington enhancement factor is then given by the ratio of the “saturated” total force (i.e., total, constant force at late times) exerted on a homogeneous atmosphere, to the total force acting on an inhomogeneous atmosphere characterized by the same mean density. Note that the Eddington enhancement factor, defined in this way, can also be interpreted as the ratio of fluxes necessary to exert the same amount of force on an inhomogeneous atmosphere as on the corresponding homogeneous atmosphere. This ratio is the same even if the homogeneous atmosphere is sub-Eddington, i.e., is only partially supported against gravity by radiation.

Figure 2 shows the Eddington enhancement factor for variable inclination of slabs (upper panels) and for changing density contrast of vertical slabs (lower panels). In all cases, the Thomson depth of the high density slabs $`\tau _h`$ is constant but the optical depth across low density regions increases from values $`\tau _l<1`$ to $`\tau _l>1`$ from left to right (see caption of Fig. 2 for more details). As expected, the Eddington factor increases as the slabs rotate toward the vertical direction because the atmosphere effectively becomes more porous (see upper panels). When the slabs are vertical, the flux enhancement factor increases as the density contrast $`\sigma _h/\sigma _l`$ becomes larger for constant mean density. This is due to the fact that the volume filling factor of the high-density gas decreases while that of the low-density gas increases, but the respective masses of the two density phases remain the same. Therefore, the mean, volume-weighted, flux is

$$F=(1f_v)F_l+f_vF_hF_l,$$
(22)
where $`f_v`$ is the volume filling factor of the dense gas and $`F_l`$ and $`F_h`$ are the fluxes propagating through tenuous and dense regions, respectively. As the density contrast increases and $`f_v`$ decreases, radiation tends to “flow” primarily through the low density channels and, therefore, more flux is necessary to exert the same total force as in the homogeneous case because radiation interacts less efficiently with tenuous gas. Quantitatively, in the diffusion limit, we have (Shaviv 1998)

$$l=\rho \frac{1}{\rho }=\left(\frac{\xi _l}{\sigma _l}+\frac{\xi _h}{\sigma _h}\right)\frac{\sigma }{\xi _l+\xi _h}.$$
(23)
When most volume is in the low-density phase but most mass is in the high-density phase (i.e., $`f_v0`$ or $`\xi _h/\xi _l0`$), then $`lf_v\rho _h/\rho _l`$ if $`f_v\rho _l/\rho _h`$. This qualitatively explains why $`l`$ decreases with $`\tau _l1`$ at constant density contrast $`\sigma _h/\sigma _l`$ (cf. third and fourth columns on Fig. 2). At small optical depth $`\tau _l`$, equation (23) would lead to very inaccurate answers. For example, equation (23) predicts $`l23`$ for vertical slabs with $`\tau _l=0.1`$ and $`\sigma _h/\sigma _l=100`$, compared to the actual value $`l5`$ and our analytic result $`4`$ from eq. (21) (cf. lower left panel). This discrepancy is due largely to neglect of the anisotropy of the radiation field, whereas our approach gives much more accurate results even in such an extreme case. Moreover, note that the “anisotropy term” in our expression for the flux enhancement factor, which is proportional to $`(\mathrm{\Delta }\sigma )^2=(\sigma _h\sigma _l)^2`$, vanishes for large Thomson depths and thus equations (21) and (15) reduce to equation (23) in the diffusion limit. We also considered “multi-stream” approximation schemes in order to account for the radiation anisotropy, but found the “intensity moment” approach developed here to be in significantly better agreement with Monte Carlo simulations.

## 5 Summary

We have considered radiative transfer deep within extremely inhomogeneous atmospheres, and have demonstrated that, under such conditions, radiation viscosity – i.e., the off-diagonal elements of the radiation stress tensor – plays an important role. Our approach is significantly more accurate than approaches based on the diffusion equation and multi-stream approximation. The technique developed here can be applied to the nonlinear evolution of radiation-driven instabilities in accretion disks. In particular, it can be used to study the dynamical coupling of matter and radiation in order to determine the characteristic length scales and density contrasts arising from “photon bubble” instability. This, in turn, will permit a self-consistent determination of the magnitude of the Eddington enhancement factor in radiation-dominated accretion disks. We also suggest that our method could be incorporated into radiation hydrodynamics codes such as the RHD module for ZEUS (Turner and Stone 2001).

This work was supported in part by NSF grant AST–9876887. We thank Nir Shaviv and Neal Turner for comments on the manuscript.


text:
# Could supermassive black holes be quintessential primordial black holes?

## I Introduction

There is growing evidence for the presence of supermassive black holes (SMBHs) in the center of most galaxies  including our own  (but see  for a more skeptical view). The origin of these black holes (and their relation to the host galaxy) is far from certain but several theories have been advanced. Of relevance is the the observational fact that there is a proportionality relation between the mass of the nuclear black hole, $`M_{bh}`$, and that of the bulge. The bulge mass appears to be about a thousand times larger the black hole mass, a relation which holds over 3 to 4 orders of magnitude.

Nevertheless this close relationship may be interpreted variously, and ultimately one is confronted with a chicken and egg problem - what came first: the host galaxies or the SMBHs? It is not inconceivable that SMBHs are purely the result of the internal galactic dynamics and their merging history; yet no one has proposed a concrete mechanism for converting stellar mass objects into objects 6 to 10 orders of magnitude larger (with the possible exception of ). But it could also be that central black holes preceded any luminous activity and that black holes led the formation of the first galactic bulges and Quasar (QSO) activity . If the latter is true one must then find an explanation for the origin of the primordial population of black holes.

In both scenarios it is inescapable that black hole growth has taken place in the recent cosmic history. Even if the budge matter is well virialized, and whether or not it fuels the SMBH, every time galaxies merge and their nuclear black holes coalesce part of the bulge matter ends up in the central black holes . Thus, as galactic merging proceeds, the comoving density of SMBHs decreases and their masses increase. In Section II we spell out the uncertainties of this chicken and egg process, identifying the conditions under which a primordial population of black holes is necessary. We then compute the density and mass profile required of this population, in order to explain the observed SMBHs.

The rest of our paper is devoted to proposing a cosmological mechanics for producing the required pre-galactic black holes. According to our theory SMBHs are descendants of primordial black holes (PBHs) produced in the very early universe. PBHs are produced, for example, at the end of inflation, in double inflation scenarios , or in first order phase transitions . To fix ideas, and although this is not essential for our paper, we shall follow  and assume that PBHs are produced at the end of a period of inflation with a blue spectrum of fluctuations (with the possibility of a running spectral index $`n_s`$). Whatever their origin, all PBHs previously considered in the literature are much lighter than SMBHs, with masses of the order $`MM_{}(T/1\mathrm{G}\mathrm{e}\mathrm{v})^2`$ for PBHs formed at temperature $`T`$. Hence only a very unrealistic phase transition, at $`T1\mathrm{M}\mathrm{e}\mathrm{v}`$ could produce SMBHs with masses of the right order of magnitude.

However the standard argument assumes that for all relevant cases PBH masses remain constant once they’re formed, since evaporation and accretion can be neglected . We show that this is not necessarily true and focus on scenarios in which significant accretion occurs during the lifetime of PBHs. Specifically we assume that the universe is pervaded by a quintessence field . Black holes cannot support static scalar fields in their vicinity and will try to “eat” them; quintessence is no exception. In the process their mass increases, so that the seeds which led to the SMBHs we observe today could be PBHs which have eaten too much quintessence.

We examine this possibility in Sections III to V. In Section III we estimate the effects of evaporation and accretion in the presence of a quintessence field. In Section IV we compute the comoving density of PBHs in our model. Finally in Sec. V we compute the PBH mass spectrum, and constrain the free parameters in our model in order to fit the requirement derived in Section II.

## II Galactic black hole accretion and merger histories

As outlined above, the close relationship between the mass of the central galactic black hole $`M_{bh}`$ and that of the galactic bulge may be interpreted variously. It could signal that the central black hole was formed by the inflow of bulge matter (stars and gas), but it could also be that the central black hole was there initially  and led the formation of galactic bulges. Most probably there was a combination of the two scenarios, and the central black hole could indeed be primordial, but at the same time it was also fed by outside matter.

In this section we show that such a combination is indeed highly feasible. We devolve the merger history of the galactic halos , and use a simple merging/ accretion prescription for the behaviour of the central black holes . We thereby show that the majority of galactic black holes today could have originated prior to the formation of significant halos, whilst still being in agreement with the observed mass evolution seen in QSO’s .

To devolve the galactic black holes, we use a merger tree approach  to establish the histories of the halos and the supermassive black holes within them. The method involves prescribing a number density distribution of halos nowadays using the Press-Schecter formula ,

$`N(M,z)dM=`$ $`\left({\displaystyle \frac{1}{2\pi }}\right)^{\frac{1}{2}}\times `$ (1)
$`\left({\displaystyle \frac{\rho }{M}}\right)\left({\displaystyle \frac{\omega }{\sigma }}\right)\left|{\displaystyle \frac{dln\sigma }{dlnM}}\right|`$ $`\mathrm{exp}\left\{{\displaystyle \frac{\omega ^2}{2\sigma ^2}}\right\}dM.`$ (2)
An adaptation of the basic equation provides the conditional probability that a halo of mass $`M_0`$ at redshift $`z_0`$ evolved from a progenitor of mass $`M_1=M_0\mathrm{\Delta }M`$ at redshift $`z_1=z_0+\mathrm{\Delta }z`$ ,

$`P(M_0,z_0|M_1,z_1)=`$ $`\left({\displaystyle \frac{1}{2\pi }}\right)^{\frac{1}{2}}\times `$ (3)
$`{\displaystyle \frac{\omega (z_0)\omega (z_1)}{[\sigma ^2(M_1)\sigma ^2(M_0)]^{\frac{3}{2}}}}`$ $`\mathrm{exp}\left\{{\displaystyle \frac{[\omega (z_0)\omega (z_1)]^2}{2[\sigma ^2(M_1)\sigma ^2(M_0)]}}\right\}.`$ (4)
We use the notation of Lacey and Cole where in equations (2) and (4), $`\rho (z)`$ is the energy density, $`\sigma ^2(M)`$ is the variance of density fluctuations on a spherical scale enclosing a mass $`M`$, and $`\omega (z)=\delta _c(1+z)`$ where $`\delta _c1.686`$ is the over density threshold for density fluctuations to collapse. We take $`\omega (z)`$ from a modelled matter power spectrum using $`H_0=75kms^1Mpc^1`$ and $`\mathrm{\Omega }_c=0.25`$,$`\mathrm{\Omega }_b=0.05`$,and $`\mathrm{\Omega }_\varphi =0.7`$. The redshift step $`\mathrm{\Delta }z`$ is mass dependent and represents a realistic merging timescale for the chosen halo, we follow  and take $`\mathrm{\Delta }z=\mathrm{\Delta }\omega (M)/\delta _c1`$ where

$`\mathrm{\Delta }\omega (M)\left({\displaystyle \frac{\sigma ^2(M)}{M}}\left|{\displaystyle \frac{dln\sigma ^2(M)}{dlnM}}\right|\right)^{1/2}.`$ (5)
As a halo is deconstructed only those with a mass greater than a limiting mass $`M_l`$ are traced, any progenitor with $`M<M_l`$ is treated as unbound matter accreted onto the halo. We consider a range of limiting mass scales $`M_l=10^910^{10}M_{}`$. This is akin to putting in a lower threshold on the halo’s velocity dispersion $`\sigma _{}v_{vir}=(GM/r_{vir})^{1/2}`$ of $`4050kms^1`$ where $`r_{vir}`$ is the virial radius and where we assume a spherical bulge for the halo so that $`M=4\pi r_{vir}^3\rho _{halo}/3`$. We initially restrict our analysis to typical galaxy mass scales today $`M10^{10}10^{12}M_{}`$.

We assume that a central black hole is present in all halos $`M>M_l`$. To assign black hole masses in each halo today, we use the powerful correlation recently found between the black hole mass, $`M_{bh}`$, at the center of galaxies and the line of sight velocity dispersion $`\sigma _{}`$ . Specifically we use the best fit relation, $`M_{bh}=1.2\times 10^8M_{}(\sigma _{}/200kms^1)^{3.75}`$.

The evolution of black holes through merging events and accretion of halo matter is a complex one, the timescales over which merging occurs will be intricately dependent on, amongst others, the size of halos and the ferocity of the merging event. While the mechanism for accretion will be dependent upon the halos properties (for example redshift and halo velocity dispersion). We do not endeavour to model these complex processes here and instead take a simplistic approach focused on current observational constraints. We assume halo mergers are violent events allowing black holes to merge on timescales significantly less than the time between halo mergers. We then model accretion using a relation proposed by  whereby a fraction of the halo gas is accreted onto the black hole,

$`M_{acc}=ϵ_{acc}(1+z)^2M_{halo}\mathrm{exp}[(\sigma _{}/300kms^1)^4]`$ (6)
where $`ϵ_{acc}`$ is the accretion efficiency factor, and the velocity dispersion dependence is introduced to account for the reduction in the halo gas’s ability to cool in the gravitational well around the centre at late times.

We are interested in the prospect that primordial black holes are progenitors for the SMBHs nowadays. We assume that only halos with $`M>M_l`$ contribute significantly to accretion onto the central black hole. Subsequently, if a halo only has progenitors of $`M<M_l`$ and retains a black hole then the black hole is assumed to have not undergone any gas accretion at higher redshifts and is treated as primordial.

Figure 1 demonstrates, for the case of a single halo, how the accretion efficiency, $`ϵ_{acc}`$, in (6) has a strong influence on how much of the black hole mass could be present at higher redshifts. In figure 2 we show the evolution of the total comoving energy and number densities for black holes with redshift. There is an intuitive play off, shown in figure 2, between the accretion efficiency and the fraction of black holes originating at early times. The higher the accretion rate, the higher the chance is of black holes being formed at late times, during the halo merging, as opposed to being primordial. As discussed in , constraints can be placed on the efficiency coefficient using inferred accretion rates from QSO luminosity functions . One can see that the QSO evolution data places a tight constraint on the accretion efficiency parameter: for $`M_l=10^{10}M_{}`$ one requires $`ϵ_{acc}1.85\times 10^6`$. In this scenario we see that the vast majority of black holes, using this accretion prescription, are present prior to halo merging.

An alternative scenario extends the galaxies able to support black holes down to dwarf scales. This can be modeled by lowering the mass threshold to $`10^9M_{}`$. Inevitably this will alter the evolutionary profile; accretion in smaller galaxies will be comparatively more efficient under this model as their gravitational well is not as prohibitive. As is shown in figure 3 this leads to a slight reduction in the predicted comoving energy density of black holes at early times. However, the general conclusions are the same in that the PBHs still contribute a significant fraction of the energy density today and comprise the majority of the total number of black holes involved in the evolution.

In figure 4 we provide the main output of this section. We show the expected mass distribution for the primordial black hole for $`M_l=10^{10}M_{}`$ and $`M_l=10^9M_{}`$. Their masses, typically $`M_{bh}10^310^6M_{}`$, are still many orders of magnitude larger than stellar mass scales, even considering the potential merging of early massive stars. We propose that such black hole masses could be generated through early time accretion of a quintessential scalar field onto PBHs. In the following sections we investigate such a scenario and show that mass distributions such as those is figure 4 can be obtained.

## III The effects of evaporation and accretion on primordial black holes

There has been intense discourse regarding whether or not PBHs are capable of accreting radiation. Carr and Hawking gave arguments for negligible accretion , but these were later disputed (see  and references therein). In any case these arguments only apply to perfect fluids and a scalar field $`\varphi `$ is not a perfect fluid (even though an isotropic, homogeneous scalar field does behaves like a perfect fluid). In Appendix 2 we show that even for the simplest potentials $`V(\varphi )`$ the field is indeed absorbed by the black hole leading to a mass increase with rate:

$$\frac{dM}{dt}=4\pi (2M)^2\dot{\varphi }_c^2$$
(7)
where $`\varphi _c`$ is the cosmological solution for $`\varphi `$. For more general potentials it may well be the case that the proportionality constant in the equation receives corrections of order 1, but these shall not matter for the rest of the argument in this paper.

From eqn. 7 we note that kinetic (but not potential) scalar energy leads to black hole growth. This is consistent with the result that the presence of a cosmological constant does not lead to equivalent growth.

Assuming now a potential of the form $`V=V_0e^{\lambda \varphi }`$, we have

$$\varphi =\frac{2}{\lambda \sqrt{8\pi }}\mathrm{log}\frac{t}{t_0}$$
(8)
(using $`G=1`$, not $`8\pi G=1`$), leading to:

$`{\displaystyle \frac{dM}{dt}}`$ $`=`$ $`\kappa {\displaystyle \frac{M^2}{t^2}}`$ (9)
$`\kappa `$ $`=`$ $`{\displaystyle \frac{8}{\lambda ^2}}`$ (10)
Eqn. 9 integrates to

$$\frac{1}{M}=\frac{1}{M_0}+\kappa \left(\frac{1}{t}\frac{1}{t_0}\right)$$
(11)
leading to the asymptotic mass:

$$M_{\mathrm{}}=\frac{M_0}{1\frac{\kappa M_0}{t_0}}$$
(12)
for black holes with a mass smaller than a critical mass $`M_{crit}=t/\kappa `$. These black holes eventually stop accreting and therefore are subdominant with respect to all others. For $`M=M_{crit}`$ BHs grow like $`t`$. Above this value the black holes would seem to grow faster that $`t`$, however clearly the approximations used must break down at this stage. In  it was shown that causality constrains these to grow like $`t`$ as well, a result we shall use for the rest of our calculations.

Hence there is a critical mass at time of formation when PBHs may grow proportionally to the horizon mass. This critical mass separates those black holes which will be of relevance for out scenario and those which won’t.

In addition PBHs may experience significant evaporation, via Hawking radiation. This leads to a decrease in their mass, at rate:

$`{\displaystyle \frac{dM}{dt}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{M^2}}`$ (13)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{15360\pi }}`$ (14)
This equation can be integrated into

$$M=(M_03\alpha (tt_0))^{1/3}$$
(15)
implying an evaporation time of

$$\tau =\frac{M_0^3}{3\alpha }$$
(16)
or

$$\frac{\tau }{10^{17}sec}\left(\frac{M}{10^{15}g}\right)^3$$
(17)
For black holes formed at temperatures much smaller than the Planck temperature this effect can be ignored.

Considering now both accretion and evaporation the black hole mass rate equation becomes

$$\frac{dM}{dt}=\frac{\alpha }{M^2}+\kappa \frac{M^2}{t^2}$$
(18)
which does not have analytical solutions. However, if the temperature at which the black holes are produced is not close to the Planck temperature the interplay of accretion and evaporation is very simple. Black holes with $`M_0>M_{crit}`$ will grow with $`Mt`$, and since their mass was never too small they never experience significant evaporation. Black holes with $`M_0<M_{crit}`$ will stop growing at some $`M_{\mathrm{}}`$, following (11). If this is smaller than $`10^{15}g`$ they will evaporate before today following (15). Accretion and evaporation happen at very distinct times, so although we do not have an exact solution to (18) it is an excellent approximation simply to glue together back to back (11) and (15).

For reheating temperatures of the order of the Planck temperature the situation is more complicated as evaporation may be significant while black holes are accreting. The result of a numerical integration is plotted in Fig. 6 (to be compared with Fig. 5). We see that the effect of evaporation is then to shift upwards the critical mass $`M_{crit}`$ above which the black hole mass scales like $`t`$. In addition, black holes with masses lower than $`M_{crit}`$ are quickly removed by the effects of evaporation even while they are accreting.

## IV PBH formation and their comoving density

Having identified the conditions required from PBHs in terms of early-time accretion and late-time merging, we now proceed to construct a cosmological scenario in which they are satisfied. Although not compulsory, for definiteness we consider an inflationary scenario with a reheating temperature $`T_r`$ and a tilted spectrum of scalar perturbations with a spectral index $`n_s>1`$. The index $`n_s`$ need not be a constant and indeed many inflationary models predict a running tilt, varying from scale to scale. Bearing this in mind we stress that the constraints upon $`n_s`$ discussed in this section refer to very small scales (the horizon size at $`T_r`$), widely different from the scales probed by CMB fluctuations or large scale structure surveys.

As shown in  for suitable $`n_s`$ one may have production of PBHs in a short time window immediately after reheating. Typically a black hole is formed if the density contrast on the horizon scale exceeds a given critical value, $`\delta >\delta _c`$. Its mass is given by:

$$M_{bh}=kM_H(\delta \delta _c)^{1/\gamma }.$$
(19)
Numerical studies with PBHs formed in a pure radiation background have identified $`\delta _c`$, $`k`$ and $`\gamma `$ ; it may be that quintessence modifies these values slightly however for simplicity we take the calculated values $`k=3.3,\gamma =0.34`$. But given the uncertainty we consider two values of $`\delta _c`$ in this section $`\delta _c=0.25,1`$. The correct value should be somewhere in between. In section V we assume the value $`\delta _c=0.67`$ to outline a precise example.

Note that the horizon scale is important because it is the relevant Jeans scale for radiation but also for the quintessence field. In Appendix 1 we present a simple model for the formation of PBHs with quintessence: some peculiarities are found, but they should not affect the rest of the argument.

The mass variance on the horizon scale at temperature $`T`$ is :

$$\sigma _H=10^4\left(10^{17}\left(\frac{T}{\mathrm{Mev}}\right)^2\right)^{\frac{1n_s}{4}}$$
(20)
To a good approximation, as shown in figure 9, we may assume that all PBHs are formed immediately after reheating, since $`\sigma _H`$ then decreases sharply making black holes rarer.

We have seen that we expect no black hole coalescence for redshifts higher than $`z10`$. Hence the comoving density of PBHs can be directly related to the probability of an accreting black hole being formed, $`p_{acc}`$. This is the probability of $`\delta >\delta _{crit}`$, i.e. the probability for black hole formation and accretion in each horizon.

$$p_{acc}=_{\delta _{crit}}^{\mathrm{}}\frac{e^{\frac{\delta ^2}{2\sigma _H^2}}}{\sqrt{2}\pi \sigma _H}𝑑\delta \frac{e^{\frac{\delta _{crit}^2}{2\sigma _H^2}}}{2\pi }$$
(21)
(where the last approximation comes from ) and where $`\delta _{crit}`$ is the value of $`\delta `$, in equation (19), that would create a black hole of critical mass, $`M_{crit}`$. The comoving density of accreting black holes is $`np_{acc}/(\frac{4}{3}\pi r_H^3)`$ where $`r_H`$ is the comoving horizon radius at $`T_r`$. The latter is given by $`r_H=2ta_0/a80(T/\mathrm{Mev})^1Pc`$ so that

$$n=4.6\times 10^{11}p_{acc}\left(\frac{T}{\mathrm{Mev}}\right)^3\mathrm{Mpc}^3$$
(22)
Combining (20), (21) and (22), and to be consistent with Fig. 2, requiring that $`n9\times 10^3\mathrm{Mpc}^3`$ we obtain a value for $`n_s`$ for each value of $`T_r`$. In Fig. 7 we plot the required $`n_s(T)`$ for the two values $`\delta _c=0.25,1`$, using a fixed value of $`M_{crit}(\lambda )`$. As mentioned above, the correct value of $`\delta _c`$ should lie somewhere in this region. In Fig. 7 we plot the required tilt for three different critical masses with fixed $`\delta _c=0.67`$.

The tilts required fit within the constraints of , and can accommodate the recent reionisation tilt constraint of  with $`n_s<1.27`$ with a high reheat temperature $`>10^{15}Gev)`$.

Note that the primordial black hole comoving density imposed in our considerations is much smaller than that considered in . Indeed the black holes studied by  do not grow (as opposed to ours ). The idea in  is to produce a larger density of much lighter black holes, suitable to promote them to candidates for dark matter. Our purpose is to generate a much lower density of much heavier PBHs, so that they could supply the primordial eggs for the merging history of SMBHs.

## V The PBH mass spectrum

As an example we calculate the mass spectrum following the methods of  to evaluate the initial mass spectrum. This involves using a Monte-Carlo technique to obtain black hole masses which, unfortunately, is computationally unfeasible for low $`p_{acc}`$ $`(<10^7)`$ required to give the correct comoving densities for $`\lambda >2`$. With this in mind, we consider two scenarios which provide the correct comoving number density of PBHs ($`9\times 10^3Mpc^3`$ assuming a limiting mass of $`10^{10}M_{}`$ in section II), one with a reheat temperature of the order $`T_r=10^{10}Gev`$, tilt $`n_s=1.33`$ and a second with reheat temperature $`T_r=10^4Gev`$, $`n_s=1.57`$ with $`\lambda =2`$ in each case. As Fig. 9 shows most black holes, and in particular those that are able to accrete, are then formed immediately after reheating, with a slight spread in their masses, which we further display in Fig. 10.

We then evolve this initial mass distribution considering accretion and evaporation in a quintessence model as described in Section III. We assume, however, that the quintessence field goes off scaling and kinates at a temperature of the order $`T=1Mev`$, since if the field continues to scale after this time one ends up with too large masses. Nevertheless if the field goes off scaling the black holes stop accreting (cf. Eqn 7), and even if the field starts scaling subsequently, BHs will not grow significantly again, since their growth pause has rendered them subcritical. Hence in what follows what we need is simply a brief pause in scaling at a temperature of around $`T=1Mev`$. We shall consider two examples here, one with $`T_{off}4Mev`$, and another with $`T_{off}=1Mev`$.

In Figs. 11 and 12 we find that the projected mass spectra at $`z10`$ (or indeed at any time after $`T_{off}`$ but before the galactic merging history started). In Fig. 11 we plot results for PBHs formed at a reheating temprature of $`T_r=10^4`$ Gev as described above. In Fig. 12 we plot the corresponding distribution for the scenarios in which $`T_r=10^{10}`$ Gev. In both cases note that the existence of a critical mass for accretion implies that the final distribution mimicks the initial one clipped at the critical mass. For the first scenario considered the cut off is at the peak of the distribution; hence the final distribution is very skewed. For the second scenario the cut off is to the left of the peak - so the final distribution is more symmetric. As expected if the field goes of scaling later, the final PBH masses are much larger: for $`T_{off}=4`$ Mev we find final PBH masses of the order of $`10^410^5M_{}`$; for $`T_{off}=1`$ Mev these masses grow to $`10^510^6M_{}`$.

In either case these plots are entirely consistent with those obtained from the merging history in section II (Fig. 4). Quintessence could have therefore provided the primordial seeds which then turn into the SMBHs we observe today.

## VI Conclusions

We have demonstrated that a scenario in which primordial black holes attain super massive size through the accretion of a cosmological scalar field is wholly consistent with current observational constraints. Such a model can generate the correct comoving number density and mass distribution for SMBH’s, given a standard prescription for late-time merging and matter accretion and with reasonable choices of tilt.

Existing schemes explaining SMBHs require very contrived choices of parameters (e.g.). The root of all the evil lies in the huge mass discrepancy between cosmic scale masses (such as those considered in previous studies of PBHs) and those of their supermassive cousins - the cosmological counterpart of the discrepancy between solar mass objects and SMBHs found in astrophysical schemes. We attempt to bridge the mass scale gap in our model by allowing primordial holes to grow by accreting quintessence. Still, we find that we have to switch off this process at a carefully tuned time, so that there is not too much growth.

The fact that this requirement is equivalent to requesting that quintessence goes off scaling just before nucleosynthesis is of some consolation: such a feature is already present in non-minimally coupled quintessence models such as those studied by . In these scenarios the present acceleration of the universe results from a coupling between quintessence and dark matter. It switches on close to the radiation to matter transition, but is affected by a long transient. This explains why the universe did not start accelerating until nowadays. The fact that this transient is symmetric around equality, and that equality is an equal number of expansion times from us and from nucleosynthesis, then makes the field kinate away at nucleosynthesis time.

In addition, in order to obtain the correct comoving density we have to tune carefully the value of the scalar tilt, $`n_s`$, on the scale of the horizon at the time of primordial black hole formation. However this fine tuning is a problem with any theory employing primordial black holes for astrophysical purposes, such as as candidates for dark matter , and is no better or worse in our theory.

In spite of these fine tuning problems we believe that this is an interesting scenario which deserves further work. The effects of angular momentum upon the whole picture are perhaps the most important next issue to consider.

Acknowledgements We would like to thank Andrew Jaffe, Lev Kofman, Andrew Liddle, David Spergel and Paul Steinhardt for their helpful thoughts and comments. RB is supported by PPARC.

## VII Appendix I - PBH creation and quintessence

In this appendix we the examine the effects of quintessence on PBH formation using the spherical model. The idea is to follow non-linear collapse of a super-Jeans size spherical region by modelling the overdense region as a portion of a Friedmann closed model pasted onto a flat model. Setting up this model entails studying the dynamics of quintessence in closed models. For completeness we shall also consider open models, which may be of relevance for modelling the void structure of out universe.

### A Quintessence in open and closed models

The relevant equations are:

$`\left({\displaystyle \frac{\dot{a}}{a}}\right)^2`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left(\rho +{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )\right){\displaystyle \frac{K}{a^2}}`$ (23)
$`\rho +3{\displaystyle \frac{\dot{a}}{a}}(\rho +p)`$ $`=`$ $`0`$ (24)
$`\ddot{\varphi }+3{\displaystyle \frac{\dot{a}}{a}}\dot{\varphi }`$ $`=`$ $`\lambda e^{\lambda \varphi }`$ (25)
where dots represent derivatives with respect to proper time. It is well know that if $`K=0`$ with $`p=w\rho `$ we have:

$`a`$ $``$ $`t^{\frac{2}{3(1+w)}}`$ (26)
$`\rho `$ $``$ $`{\displaystyle \frac{1}{a^{3(1+w)}}}`$ (27)
$`w_\varphi `$ $`=`$ $`w`$ (28)
$`\mathrm{\Omega }_\varphi `$ $`=`$ $`{\displaystyle \frac{\rho _\varphi }{\rho +\rho _\varphi }}={\displaystyle \frac{3}{\lambda ^2}}(1+w)`$ (29)
$`\varphi `$ $`=`$ $`{\displaystyle \frac{2}{\lambda }}\mathrm{log}{\displaystyle \frac{t}{t_0}}`$ (30)
We call this solution the scaling solution. With $`K0`$ this is also the solution at early times, when curvature is subdominant.

Without quintessence open models ($`K=1`$) at late times become vacuum dominated - the so called Milne Universe for which $`at`$. However we find that in the presence of quintessence the onset of negative curvature domination leads to another scaling solution: one in which curvature and quintessence remain proportional. We find that:

$`a`$ $``$ $`t`$ (31)
$`\rho _\varphi `$ $``$ $`{\displaystyle \frac{1}{a^2}}`$ (32)
$`w_\varphi `$ $`=`$ $`{\displaystyle \frac{1}{3}}`$ (33)
$`\mathrm{\Omega }_\varphi `$ $`=`$ $`{\displaystyle \frac{\rho _\varphi }{\rho _c}}={\displaystyle \frac{2}{\lambda ^2}}`$ (34)
$`\varphi `$ $`=`$ $`{\displaystyle \frac{2}{\lambda }}\mathrm{log}{\displaystyle \frac{t}{t_0}}`$ (35)
which implies that open universes at late times are devoid of normal matter, but not of quintessence. This type of behavior can be understood from the $`K=0`$ scaling solution, considering that open curvature behaves like a fluid with $`w=1/3`$. In the same way that quintessence locks on to matter or radiation in a scaling solution, it locks on to open curvature.

As is well known, without quintessence closed models ($`K=1`$) expand and eventually turn around and collapse in a Big Crunch. This type of behavior is not changed by the presence of quintessence; however the scaling behavior of quintessence itself is drastically changed. As the universe comes to an halt the friction term in the $`\varphi `$ equations (that is the term $`3\frac{\dot{a}}{a}\dot{\varphi }`$) is withdrawn and the field becomes kinetic energy dominated (i.e. it kinates). Hence as the universe turns around quintessence becomes subdominant, as it scales like $`\rho _\varphi 1/a^6`$ in contrast with $`\rho 1/a^3`$.

However as the universe enters the contracting phase what used to be a friction term starts to drive the field, since now $`\frac{\dot{a}}{a}<0`$. This leading to runaway kination, since it is precisely the balance of braking and the slope of $`V(\varphi )`$ what usually moderates the balance of kinetic and potential energy. As before, kinetic energy domination implies $`w_\varphi =1`$ and a faster decay rate with expansion $`\rho 1/a^6`$. However contraction reverses the argument, and whatever has stronger dilution rate during expansion, will have higher compression rate during contraction. Therefore, at late stages of collapse quintessence dominates. Since curvature and background matter can be ignored we have the solution:

$`a`$ $``$ $`(t_ct)^{1/3}`$ (36)
$`\rho _\varphi `$ $``$ $`{\displaystyle \frac{1}{a^6}}`$ (37)
$`w_\varphi `$ $`=`$ $`1`$ (38)
$`\mathrm{\Omega }_\varphi `$ $``$ $`1`$ (39)
$`\varphi `$ $``$ $`a^{1/\sqrt{6}}`$ (40)
in which $`t_c`$ is the crunch time.

### B Implications for structure formation and PBH formation

Qualitatively these results indicate that quintessence has a leading role in the strongly non-linear stages of structure formation. Voids should be filled with quintessence, judging from what happens to the $`K=1`$ case. Also it would appear that black hole formation would be led by quintessence and not by matter (as implied by the $`K=1`$ case). Interestingly, quintessence appears to behave like a stiff fluid ($`p_\varphi =\rho _\varphi `$) during collapse, so we should be able to simplify the calculations. The black hole mass seems to be dominated by the amount of quintessence accreted. Also because $`M_\varphi \rho _\varphi a^31/a^3`$, there must be a mass enhancement during collapse (due to gravity acting against quintessence’s pressure).

All of these effects may at most introduce factors of order 1 in the calculations in the main body of the text, and so interesting as these results might be we have relegated them to Appendix.

## VIII Appendix II - Quintessence around a black hole

Given the large Jeans mass of quintessence, its fluctuations should be very small even in the interior of non-linear objects, such as the Solar system. Indeed the relativistic restoration force induced by the $`\varphi `$ gradients should ensure that quintessence fluctuations remain linear, even under highly non-linear gravitational forces. One obvious exception is the vicinity of a BHs horizon, where $`\varphi `$ has to change drastically.

The interaction of BHs and quintessence then becomes a problem of boundary matching: on one side the homogeneous cosmological solution; on the other the infinite radius limit of the Schwarzchild solution. This matter has been examined in the literature in the context of the gravitational memory problem in Brans-Dicke theories . Our analysis will closed mimic that of Jacobson .

We assume that the cosmological scalar field generates a much weaker gravitational field than the black hole, so that we can impose a Schwarzchild metric, with a quasi-stationary mass parameter. Hence the equation for the scalar field is the vicinity of the BH is:

$$\mathrm{}\varphi =\frac{\ddot{\varphi }}{A}+\frac{1}{r^2}(r^2A\varphi ^{})^{}=\frac{V}{\varphi }$$
(41)
where $`A(r)=12M/r`$, and dots and dashes are derivatives with respect to time and $`r`$ respectively. For free waves ($`V=0`$) the system separates giving:

$$\varphi =e^{i\omega t}R(r)$$
(42)
with:

$$\frac{1}{r^2}(r^2AR^{})^{}+\frac{\omega ^2}{A}R=0$$
(43)
or introducing a Kruskal coordinate $`r^{}`$, such that $`Adr^{}=dr`$ (or $`r^{}=r+2M\mathrm{log}(\frac{r}{2M}1)`$):

$$\frac{d^2R}{dr^2}+\frac{2A}{r}\frac{dR}{dr^{}}+\omega ^2R=0$$
(44)
It is clear that far away the solutions are:

$$\varphi =\frac{e^{i\omega (t\pm r)}}{r}$$
(45)
and that near the horizon they become:

$$\varphi =e^{i\omega (t\pm r^{})}$$
(46)
Focusing on solutions regular on the horizon $`^+`$, we introduce the advanced time coordinate $`v=t+r^{}`$, so that the relevant oscillatory solutions take the form $`e^{i\omega v}`$.

In addition there are non-oscillatory solutions regular at $`r=2M`$, such as those studied by Jacobson. These can be obtained noting that (41) can be rewritten as:

$$2\varphi _{,rv}+\frac{2}{r}\varphi _{,v}+\frac{1}{r^2}(r^2A\varphi _{,r})_{,r}=V^{}(\varphi )$$
(47)
For a general $`V(\varphi )`$ separability is lost and numerical work is necessary. In order to further the analytical approach we model the rolling potential by its gradient at a point $`V=\mu \varphi `$ (we could add a constant here without loss of generality). Equation (47) then becomes:

$$2\varphi _{,rv}+\frac{2}{r}\varphi _{,v}+\frac{1}{r^2}(r^2A\varphi _{,r})_{,r}=\mu $$
(48)
Setting $`\varphi =f(v)+g(r)`$ leads to $`f=Bv+D`$, and

$$g_{,r}=\frac{\mu r^2}{3(r2M)}\frac{Br}{r2M}+\frac{C}{r(r2M)}$$
(49)
For the solution to be regular at the horizon we must have:

$$C=4M^2\left(B+\frac{2}{3}\mu M\right)$$
(50)
Integrating finally leads to:

$`\varphi `$ $`=`$ $`B[vr2M\mathrm{log}r]`$ (52)
$`{\displaystyle \frac{\mu }{3}}\left({\displaystyle \frac{r^2}{2}}+2Mr+4M^2\mathrm{log}r\right)+D`$
which generalized Jacobson’s solution. The general solution is a superposition of free oscillatory solutions and this solution.

We must now impose a boundary condition of the form:

$$\varphi =\varphi _c+\dot{\varphi }_ct$$
(53)
where $`\varphi _c`$ refers to the quintessence cosmological solution. If we fix $`B`$ and $`D`$ so that:

$`\varphi `$ $`=`$ $`\varphi _c+\dot{\varphi }_c\left[vr2M\mathrm{log}{\displaystyle \frac{r}{2M}}\right]`$ (55)
$`{\displaystyle \frac{\mu }{3}}\left({\displaystyle \frac{r^2}{2}}+2Mr+4M^2\mathrm{log}{\displaystyle \frac{r}{2M}}\right)`$
we have asymptotically:

$$\varphi =\varphi _c+\dot{\varphi }_ct\frac{\mu }{6}r^2$$
(56)
A small error is made in matching the two conditions, but in a quintessence scenario with

$$\mu \frac{V}{\varphi }=\lambda V=\lambda \mathrm{\Omega }_\varphi \rho $$
(57)
(where $`\lambda `$ has units of $`1/\varphi `$ and is expressed in Planck units and $`\rho `$ is the cosmological density) we find thatthe scale on which the error becomes significant is the cosmological horizon scale. Hence this small error can be neglected.

Computed the flux of energy through the BH horizon associated with this solution:

$$=T_{vv}=\dot{\varphi }_c^2$$
(58)
we finally find an equation for the BH mass:

$$\frac{dM}{dt}=4\pi (2M)^2\varphi _c^2$$
(59)
which is the same formula derived by Jacobson. Hence we conclude that the potential $`V(\varphi )`$ only has indirect impact upon the BH mass, via its effect upon the time evolution of the cosmological solution $`\varphi _c`$.