text:
# The wave equation on the Schwarzschild metric II: Local Decay for the spin 2 Regge Wheeler equation

## 1 Co-ordinates and equations

The Schwarzschild manifold describes a static black hole solution to the Einstein equation. The exterior of the black hole is most easily described by $`(t,r,\theta ,\varphi )𝐑\times (2M,\mathrm{})\times S^2`$ with the metric

$$ds^2=(1\frac{2M}{r})dt^2(1\frac{2M}{r})^1dr^2r^2ds_{S^2}^2$$
(1.1)
To simplify the problem, Regge and Wheeler () introduced a new radial co-ordinate, $`r_{}`$, satisfying

$$\frac{r}{r_{}}=(1\frac{2M}{r})$$
(1.2)
This allows the definition of a space like manifold

$$𝔐=𝐑\times S^2$$
(1.3)
The old co-ordinate $`r`$ is now treated as a function of $`r_{}`$.

In these new co-ordinates, the Regge-Wheeler equation for a scalar field $`u:𝐑\times 𝔐𝐑`$ which determines the behavior of the odd-type tensor fields is

$$\ddot{u}+Hu=0$$
(1.4)
where

$`H=`$ $`{\displaystyle \underset{j=1}{\overset{3}{}}}H_j`$ (1.5)
$`H_1=`$ $`{\displaystyle \frac{^2}{r_{}^2}}`$ (1.6)
$`H_2=`$ $`(1s^2)V`$ (1.7)
$`V=`$ $`{\displaystyle \frac{2M}{r^3}}(1{\displaystyle \frac{2M}{r}})`$ (1.8)
$`H_3=`$ $`V_L(\mathrm{\Delta }_{S^2})=V_L{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}l(l+1)P_l`$ (1.9)
$`V_L=`$ $`{\displaystyle \frac{1}{r^2}}(1{\displaystyle \frac{2M}{r}})`$ (1.10)
and where $`s=2`$ for the case of the tensor field and $`P_l`$ is projection onto spherical harmonics with total angular momentum $`l`$. The case $`s=0`$ is the scalar field previously considered and $`s=1`$ is for the odd-type vector (Maxwell) fields.

Because of the way the scalar field $`u`$ is defined it is not possible for it to have any component with spherical harmonic component $`l=0`$. It has also been shown that the $`l=1`$ component corresponds to changing the non-rotating Schwarzschild background to a rotating Kerr solution and to gauge transformations . For this reason, we only consider $`u`$ with no $`l=0`$ or $`l=1`$ spherical harmonic component. This provides a lower bound on the spherical Laplace-Beltrami operator.

$$\mathrm{\Delta }_{S^2}2(2+1)=6$$
(1.11)
For the scalar wave equation, Bachelot and Nicolas have proven global existence  in both an energy space and in $`C^{\mathrm{}}`$. The assumption of global existence in $`C^{\mathrm{}}`$ greatly simplifies all the following arguments and will be assumed although we are not yet aware of a published proof. However, the method of Bachelot and Nicolas should extend to the Regge-Wheeler equation without difficulty. The assumption of global existence in $`C^{\mathrm{}}`$ means that all solutions are assumed to be $`C^{\mathrm{}}(𝔐)H^1(𝔐,dr_{}d^2\omega _{S^2})`$, are infinitely differentiable in $`t`$ and have time derivative in $`C^{\mathrm{}}(𝔐)L^2(𝔐,dr_{}d^2\omega _{S^2})`$. The notation $`u(t)`$ denotes the function from $`𝔐𝐑`$ corresponding to $`u`$ evaluated at time $`t`$. The measure $`dr_{}d^2\omega _{S^2}`$ is used for all norms and inner products unless otherwise specified.

## 2 The Heisenberg-type relation and preliminary estimates

For the Schrödinger equation, the Heisenberg relation describes the time evolution of expectation values for an operator and gives conserved quantities from symmetries of the Hamiltonian. A similar relation exists for the wave equation .

###### Theorem 2.1 (Heisenberg-like relation).

For a time independent operator $`A`$ and a solution to the linear wave equation $`\ddot{u}+Hu=0`$ such that $`u`$ and $`Hu`$ are in the domain of $`A`$, and $`u`$ and $`Au`$ are in the domain of $`H`$,

$$\frac{d}{dt}(u,A\dot{u}\dot{u},Au)=u,[H,A]u$$
(2.12)
Proof: in Theorem 1 of

The first and most important application of theorem 2.1 is conservation of energy. As usual it is generated by time translation symmetry. This result is already well known .

###### Theorem 2.2 (Energy Conservation).

The Regge-Wheeler equation, equation 1.4, has a conserved quantity $`u_{}^2`$ which we call the energy.

$$u_{}^2=\dot{u},\dot{u}+u^{},u^{}+u,3Vu+\underset{l=2}{\overset{\mathrm{}}{}}P_lu,l(l+1)V_LP_lu$$
(2.13)
This acts as a metric on the space $`=\{uL^2:u_{}<\mathrm{}\}`$.

Proof: The conservation of energy follows from the Heisenberg-like relation with the multiplier $`A=\frac{d}{dt}`$. This acts as a metric because the $`u,(H_2+H_3)u`$ term is positive. The positivity of this potential was known to Regge and Wheeler  and is verified here.

Since only functions orthogonal to the spherical harmonics with $`l=0`$ and $`l=1`$ are considered

$`u,3Vu+{\displaystyle \underset{l=2}{\overset{\mathrm{}}{}}}P_lu,l(l+1)V_lP_lu`$ $`u,3Vu+u,6V_lu`$
$``$ $`u,({\displaystyle \frac{6M}{r^3}}+{\displaystyle \frac{6}{r^2}})(1{\displaystyle \frac{2M}{r}})u`$
$``$ $`u,{\displaystyle \frac{6}{r^2}}(1{\displaystyle \frac{M}{r}})(1{\displaystyle \frac{2M}{r}})u`$
Since $`\frac{6}{r^2}(1\frac{M}{r})(1\frac{2M}{r})`$ is always positive $`u,(H_2+H_3)u`$ is positive definite. Since $`\dot{u},\dot{u}`$, $`u^{},u^{}`$, and $`u,(H_2+H_3)u`$ are all positive each is defined if $`u_{}`$ is finite, and $`u_{}`$ is a metric.

$`\mathrm{}`$

As stated in the proof of energy conservation, the energy controls certain derivative norms and this can be used to control the growth of the $`L^2`$ norm.

###### Theorem 2.3.

If $`u`$ is a real valued solution of the Regge-Wheeler equation (equation 1.4) and $`u(t)_{L^2}`$ is the norm of $`u`$ at time $`t`$, then for $`t0`$

$`\dot{u}_{L^2}`$ $`u_{}`$
$`u(t)`$ $`tu_{}+u(0)_{L^2}`$
Proof: Since $`u^{},u^{}`$ and $`u,(H_2+H_3)u`$ are strictly positive $`\dot{u}(t)`$ is controlled by the energy. This is used to control the growth rate of $`u(t)_{L^2}`$

$`{\displaystyle \frac{d}{dt}}u(t)_{L^2}^2=`$ $`{\displaystyle \frac{d}{dt}}u,u=2u,\dot{u}`$
$`2u_{L^2}{\displaystyle \frac{d}{dt}}u_{L^2}`$ $`2u_{L^2}\dot{u}_{L^2}`$
$`{\displaystyle \frac{d}{dt}}u_{L^2}`$ $`\dot{u}_{L^2}u_{}`$
$`\mathrm{}`$

## 3 Local Decay

For the scalar wave equation, a radial differential operator $`\gamma `$ was introduced to prove the local decay estimate . It is used here for the same purpose. This multiplier is centered at the peak of the angular potential $`V_L`$. To simplify calculations the standard $`r_{}`$ co-ordinate is translated to a new one, $`\rho _{}=r_{}\alpha _{}`$ so that $`\rho _{}=0`$ at the peak of $`V_L`$. This centrifugal tortoise co-ordinate satisfies the same differential definition as $`r_{}`$, equation 1.2.

###### Definition 3.1.

The centrifugal origin $`\alpha `$, centrifugal tortoise origin $`\alpha _{}`$, and centrifugal tortoise radius $`\rho _{}`$ are defined by

$`\alpha `$ $`3M`$ (3.14)
$`\alpha _{}`$ $`r_{}|_{r=\alpha =3M}`$ (3.15)
$`\rho _{}`$ $`r_{}\alpha _{}`$ (3.16)
###### Definition 3.2.

Given $`\sigma (\frac{1}{2},1]`$, the Morawetz type multiplier $`\gamma _\sigma `$ is defined by

$`g_\sigma (r_{})`$ $`{\displaystyle _0^{\frac{r_{}\alpha _{}}{2M}}}(1+\tau ^2)^\sigma 𝑑\tau `$ (3.17)
$`\gamma _\sigma `$ $`{\displaystyle \frac{i}{2}}(g_\sigma {\displaystyle \frac{}{r_{}}}+{\displaystyle \frac{}{r_{}}}g_\sigma )`$ (3.18)
As before, $`C^{\mathrm{}}`$ solutions are assumed so that there are no domain issues. In all cases. the value of $`\sigma `$ will be fixed and the notation $`g=g_\sigma `$ and $`\gamma =\gamma _\sigma `$ will be used.

###### Theorem 3.3.

If $`uH^1(𝔐)`$ and $`\sigma (\frac{1}{2},1]`$, then

$$u,\gamma _\sigma u=0$$
(3.19)
and there is a constant $`C_\sigma =lim_r_{}\mathrm{}g(r_{})`$ such that

$$\gamma uC_\sigma u_{}+\frac{1}{2}(1+(\frac{r_{}\alpha _{}}{2M})^2)^\sigma u_{L^2}$$
(3.20)
Proof: Equation 3.19 is proven in theorem 16 of  (the statement of which includes the additional, but unnecessary, assumption that $`u`$ satisfy the scalar wave equation). For equation 3.20, theorem 17 of  does not directly apply since the space $``$ defined there involves different potentials. However, the same argument applies. It is first noted that since $`\sigma >\frac{1}{2}`$ and the integrand in the definition of $`g`$ is positive and even, $`|g|`$ is bounded by $`C_\sigma =lim_r_{}\mathrm{}g(r_{})`$. Now, by direct computation,

$`\gamma u=`$ $`gu^{}+{\displaystyle \frac{1}{2}}g^{}u`$
$``$ $`gu^{}+{\displaystyle \frac{1}{2}}g^{}u`$
$``$ $`C_\sigma u_{}+{\displaystyle \frac{1}{2}}(1+({\displaystyle \frac{r_{}\alpha _{}}{2M}})^2)^\sigma u_{L^2}`$
$`\mathrm{}`$

The Heisenberg-like relation will be applied to the multiplier $`\gamma `$. To do this it is necessary to estimate the commutator $`[_{j=1}^3H_j,\gamma ]`$.

###### Lemma 3.4.

For $`\sigma (\frac{1}{2},1]`$

$$i[\frac{1}{r^2}(1\frac{2M}{r})\mathrm{\Delta }_{S^2},\gamma ]=g(r_{})[\frac{3M}{r}1]\frac{2}{r^3}(1\frac{2M}{r})\mathrm{\Delta }_{S^2}0$$
(3.21)
Proof: in Lemma 18 of

###### Lemma 3.5.

For $`\sigma (\frac{1}{2},1]`$ and $`u`$ in the domain of $`\gamma `$ and $`H`$

$$u,i[\frac{^2}{r_{}^2},\gamma ]uu,\frac{\sigma }{(1+(\frac{\rho _{}\alpha _{}}{2M})^2)^{\sigma +2}}\frac{1}{(2M)^3}[5+(32\sigma )(\frac{\rho _{}\alpha _{}}{2M})^2]u$$
(3.22)
Proof: in Lemma 20 of

###### Lemma 3.6.

For $`\sigma (\frac{1}{2},1]`$ and $`u`$ in the domain of $`H`$ and $`\gamma `$, there is a constant $`c_\sigma `$ so that

$$u,i[\underset{i=j}{\overset{3}{}}H_j,\gamma ]u,\frac{c_\sigma }{(1+(\frac{r_{}\alpha _{}}{2M})^2)^{\sigma +1}}u$$
(3.23)
Proof: Since $`\sigma 1`$

$$\frac{\sigma }{(1+(\frac{\rho _{}\alpha _{}}{2M})^2)^{\sigma +2}}\frac{1}{(2M)^3}[5+(32\sigma )(\frac{\rho _{}\alpha _{}}{2M})^2]<\frac{\sigma }{(1+(\frac{\rho _{}}{2M})^2)^{\sigma +1}}$$

In the proof of lemma 21 of  it is shown that

$$i[V,\gamma ]=g(3\frac{8M}{r})\frac{2M}{r^4}(1\frac{2M}{r})$$
(3.24)
Since $`g<0`$ for $`r<3M`$ and $`g>0`$ for $`r>3M`$, $`i[H_2,\gamma ]=3i[V,\gamma ]`$ is negative for $`r<\frac{8M}{3}`$, positive for $`\frac{8M}{3}<r<3M`$, and negative for $`3M<r`$. In the region $`\frac{8M}{3}<r<3M`$, all the terms of the form $`i[H_j,\gamma ]`$ are positive so an estimate of the form 3.23 holds. The other $`r`$ values are now treated.

It is useful to note that a term relating $`H_2`$ to $`H_3`$ is decreasing since

$$\frac{d}{dr}\frac{3\frac{8M}{r}}{1\frac{3M}{r}}=\frac{M}{(r3M)^2}<0$$

At $`r=2M`$, $`3\frac{8M}{r}=1=2(1\frac{3M}{r})`$. Therefore in $`2M<r\frac{8M}{3}`$ and for $`l2`$

$`|3{\displaystyle \frac{8M}{r}}|`$ $`2|1{\displaystyle \frac{3M}{r}}|`$
$`|r{\displaystyle \frac{1}{r^3}}(1{\displaystyle \frac{2M}{r}})({\displaystyle \frac{2M}{r}})(3{\displaystyle \frac{8M}{r}})|<`$ $`|g{\displaystyle \frac{1}{r^3}}(1{\displaystyle \frac{2M}{r}})2(1{\displaystyle \frac{3M}{r}})|`$
$`|i[H_2,\gamma ]|<`$ $`{\displaystyle \frac{1}{6}}i[H_3,\gamma ]`$
The $`\frac{1}{6}`$ factor is present due to the restriction that $`l2`$ and hence $`\mathrm{\Delta }_{S^2}6`$.

At $`r=3.3M`$

$$\left|\frac{2(1\frac{3M}{r})6}{3(3\frac{8M}{r})}\right|=\frac{1\frac{1}{11}}{1\frac{8}{11}}=\frac{12}{19}>\frac{2M}{3.3M}$$

Therefore for $`r>3.3M`$ and $`l2`$

$`|3(3{\displaystyle \frac{8M}{r}}){\displaystyle \frac{2M}{3.3M}}|<`$ $`2(1{\displaystyle \frac{3M}{r}})6`$
$`|{\displaystyle \frac{g}{r^3}}(1{\displaystyle \frac{2M}{r}})(3)(3{\displaystyle \frac{8M}{r}}){\displaystyle \frac{2M}{r}}|<`$ $`{\displaystyle \frac{g}{r^3}}(1{\displaystyle \frac{2M}{r}})(1{\displaystyle \frac{3M}{r}})(\mathrm{\Delta }_{S^2})`$
$`|i[H_2,\gamma ]|<`$ $`i[H_3,\gamma ]`$
Finally for $`3Mr3.3M`$, since $`i[H_3,\gamma ]`$ vanishes quadratically in $`(r3M)`$ where as $`i[H_2,\gamma ]`$ vanishes only linearly it is necessary to bound $`i[H_2,\gamma ]`$ by $`i[H_1,\gamma ]`$. On this interval $`F^1<3`$, $`\rho _{}<.9M`$, and $`g<.9M`$. Again, assuming $`l2`$

$`|i[H_2,\gamma ]|=`$ $`3g(3{\displaystyle \frac{8M}{r}})({\displaystyle \frac{2M}{r}})^4{\displaystyle \frac{1}{(2M)^3}}(1{\displaystyle \frac{2M}{r}})`$
$`<(3)(.9)({\displaystyle \frac{19}{33}})({\displaystyle \frac{2}{3}})^4({\displaystyle \frac{13}{33}}){\displaystyle \frac{1}{(2M)^3}}<.121{\displaystyle \frac{1}{(2M)^3}}`$
$`i[H_3,\gamma ]>`$ $`{\displaystyle \frac{5\sigma }{(1+(\frac{\rho _{}}{2M})^2)^{\sigma +2}}}{\displaystyle \frac{1}{(2M)^3}}>{\displaystyle \frac{5\sigma }{(1.2025)^{\sigma +2}}}{\displaystyle \frac{1}{(2M)^3}}>1.43{\displaystyle \frac{1}{(2M)^3}}`$
In summary, for $`r<3M`$ and $`r>3.3M`$, $`i[H_2+H_3,\gamma ]>0`$ and for $`3Mr3.3M`$, $`i[H_1+H_2,\gamma ]`$ is strictly positive. Since for $`r<3M`$ and $`r>3.3M`$, $`i[H,\gamma ]>i[H_1,\gamma ]>C(1+(\frac{\rho _{}}{2M})^2)^{\sigma 1}`$ and for $`3Mr3.3M`$, $`i[H,\gamma ]`$ is strictly positive, there is a constant $`C`$ so that

$$u,i[H,\gamma ]uu,\frac{C}{(1+(\frac{r_{}\alpha _{}}{2M})^2)^{\sigma +1}}u$$

$`\mathrm{}`$

It is now possible to apply the Heisenberg type relation to $`\gamma `$ and integrate the result to prove local decay.

###### Theorem 3.7 (Local Decay).

If $`u`$ is a solution to the Regge-Wheeler equation (equation 1.4), $`u_{}=E`$, $`u(0)=f`$, and $`\beta >\frac{3}{2}`$, then there is a constant $`D_\sigma `$ such that

$$_0^{\mathrm{}}(1+(\frac{r_{}}{2M})^2)^{\frac{\beta }{2}}u^2𝑑tD_\sigma E^{\frac{1}{2}}(E^{\frac{1}{2}}+f_{L^2})$$
(3.25)
Proof: Initially the result will be proven with $`\beta =\sigma +1`$ and $`\sigma (\frac{1}{2},1]`$ and $`\rho _{}`$ in place of $`r_{}`$. By integrating lemma 3.6 and applying the Heisenberg like relation, Theorem 2.1, it is possible to bound the time integral of the local decay term by an inner product evaluated at time $`T`$. Despite the explicit factors of $`i`$ appearing in the following, all terms are real valued.

$`{\displaystyle _0^{\mathrm{}}}(1+({\displaystyle \frac{r_{}\alpha _{}}{2M}})^2)^{\frac{\sigma +1}{2}}u^2𝑑t`$
$`{\displaystyle _0^T}i[H,\gamma ]𝑑t`$
$`{\displaystyle _0^T}{\displaystyle \frac{d}{dt}}(u,i\gamma \dot{u}\dot{u},i\gamma u)𝑑t`$
$`{\displaystyle _0^T}i{\displaystyle \frac{d}{dt}}({\displaystyle \frac{d}{dt}}u,\gamma u2\dot{u},\gamma u)𝑑t`$
$`{\displaystyle _0^T}2i{\displaystyle \frac{d}{dt}}\dot{u},\gamma udt`$
$`2(\dot{u}\gamma u)|_{t=T}+2(\dot{u}\gamma u)|_{t=0}`$
$`E^{\frac{1}{2}}(4C_\sigma E^{\frac{1}{2}}+(1+({\displaystyle \frac{r_{}\alpha _{}}{2M}})^2)^\sigma f)+E^{\frac{1}{2}}(1+({\displaystyle \frac{r_{}\alpha _{}}{2M}})^2)^\sigma u(T)`$ (3.26)
Since $`\sigma >\frac{1}{2}`$, $`q`$ can be chosen so that $`\frac{1}{2\sigma }+\frac{1}{2}<q<\frac{3}{2}`$. If $`p`$ is the conjugate exponent to $`q`$ and $`\kappa \frac{2}{p}`$, then

$$\frac{1}{p}>1\frac{2}{3}=\frac{1}{3},\frac{2\kappa }{2}q=1,q\sigma >\frac{\sigma +1}{2}$$

Hölder’s inequality can now be applied to the last norm in line 3.26.

$`(1+({\displaystyle \frac{r_{}\alpha _{}}{2M}})^2)^\sigma u^2=`$ $`{\displaystyle _𝔐}{\displaystyle \frac{|u|^\kappa |u|^{2\kappa }}{(1+(\frac{\rho _{}}{2M})^2)^{2\sigma }}}𝑑r_{}d^2\omega _{S^2}`$
$``$ $`\left({\displaystyle _𝔐}|u|^{p\kappa }𝑑r_{}d^2\omega _{S^2}\right)^{\frac{1}{p}}\left({\displaystyle _𝔐}{\displaystyle \frac{|u|^{(2\kappa )q}}{(1+(\frac{\rho _{}}{2M})^2)^{2\sigma q}}}𝑑r_{}d^2\omega _{S^2}\right)^{\frac{1}{q}}`$
$`(1+({\displaystyle \frac{r_{}\alpha _{}}{2M}})^2)^\sigma u`$ $`u^{\frac{1}{p}}{\displaystyle \frac{|u|}{(1+(\frac{\rho _{}}{2M})^2)^{\sigma q}}}^{1\frac{1}{p}}`$
$``$ $`(E^{\frac{1}{2}}T+f)^{\frac{1}{p}}{\displaystyle \frac{u}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}}}^{1\frac{1}{p}}`$
For sufficiently large $`T`$, there is a constant $`F`$ so that

$$(1+(\frac{r_{}\alpha _{}}{2M})^2)^\sigma uFT^{\frac{1}{p}}\frac{u}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}}^{1\frac{1}{p}}$$

$`{\displaystyle _0^T}(1+({\displaystyle \frac{\rho _{}}{2M}})^2)^{\frac{\sigma +1}{2}}u^2𝑑t`$ $`E^{\frac{1}{2}}(4C_\sigma E^{\frac{1}{2}}+f)+FT^{\frac{1}{p}}{\displaystyle \frac{u}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}}}^{1\frac{1}{p}}`$ (3.27)
This establishes an integral relation between the local decay norm and its square integral. The local decay norm has bounded derivative since

$$\frac{d}{dt}\frac{u}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}}^2=2\frac{u}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}},\frac{\dot{u}}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}}2\frac{u}{(1+(\frac{\rho _{}}{2M})^2)^{\frac{\sigma +1}{2}}}E^{\frac{1}{2}}$$
(3.28)
These two conditions are sufficient to apply lemma 25 of . That lemma states that for $`\theta :𝐑𝐑^+`$ with uniformly bounded derivative, $`ϵ(0,\frac{1}{3})`$, if $`_0^t\theta (\tau )^2𝑑\tau C_1+C_2t^ϵ\theta ^{1ϵ}`$ then $`t^ϵ\theta (t)^{1ϵ}`$ goes to zero sequentially and hence $`_0^t\theta (\tau )^2𝑑\tau C_1`$. The lemma can be applied with $`\theta `$ as the local decay norm, $`\frac{1}{p}=ϵ`$, and $`C_1`$ and $`C_2`$ as in 3.27.

This proves the result for $`\beta (\frac{3}{2},2]`$ and for $`\rho _{}`$ instead of $`r_{}`$. Since $`(1+(\frac{\rho _{}}{2M})^2)^\beta `$ is a decreasing function of $`\beta `$, the result holds for all $`\beta >\frac{3}{2}`$. Finally since for any $`\beta `$ there is a constant so that for all $`r_{}`$, $`(1+(\frac{r_{}\alpha _{}}{2M})^2)^\beta C(1+(\frac{r_{}}{2M})^2)^\beta `$ the statement of the theorem holds. $`\mathrm{}`$

## Appendix A Numerical verification of the positivity of the commutator

The key step in proving the local decay estimate is the lower bound for the commutator $`i[H,\gamma ]`$ proven in lemma 3.6. From the asymptotics of $`i[H_2,\gamma ]r^4`$ and $`i[H_3,\gamma ]r^3`$ it is clear that the negative contributions from $`i[H_2,\gamma ]`$ will be dominated eventually and it is sufficient to show $`i[H,\gamma ]`$ is positive in some finite domain. To verify positivity of the commutator, the sum of the exact form for $`i[H_2,\gamma ]`$ from equation 3.24, the lower bounds for $`i[H_1,\gamma ]`$ from equation 3.22, and the lower bound for $`i[H_3,\gamma ]`$ from equation 3.21 and $`l2`$ is plotted for $`M=1`$ and $`\sigma =1`$. From the graph it is clear that the total commutator is positive. The graph decays because all the terms involved decay.

Acknowledgments: This work was completed when the authors were visiting the Fields Institute in Toronto; we thank the organizers of the special PDE program for their hospitality. We would like to thank the organizers of the Centre de Recherches Mathematiques’s workshop on the interaction of gravity with classical fields where this project was started. We thank B.F. Whiting and S. T. Yau for valuable discussions.

This work was partially supported by NSF grant # DMS-0100490.


format: text

text:
# On magnetic field generation in Kolmogorov turbulence

## Abstract

We analyze the initial, kinematic stage of magnetic field evolution in an isotropic and homogeneous turbulent conducting fluid with a rough velocity field, $`v(l)l^\alpha `$, $`\alpha <1`$. We propose that in the limit of small magnetic Prandtl number, i.e. when ohmic resistivity is much larger than viscosity, the smaller the roughness exponent $`\alpha `$, the larger the magnetic Reynolds number that is needed to excite magnetic fluctuations. This implies that numerical or experimental investigations of magnetohydrodynamic turbulence with small Prandtl numbers need to achieve extremely high resolution in order to describe magnetic phenomena adequately.

1. Introduction. In a turbulent highly conducting fluid, magnetic fields may be amplified since the field lines are generally stretched by randomly moving fluid elements in which these lines are frozen . Such a mechanism of turbulent dynamo is expected to work in a variety of astrophysical systems (galaxy clusters, interstellar medium, stars, planets), is confirmed numerically, and is consistent with simple analytical models.

Valuable insight can be gained from considering the so-called kinematic stage of the dynamo, when magnetic field is amplified from an initially weak ‘seed’ field. As the simplest example, assume that the resistive scale is smaller than or of the same order as the viscous scale of the fluid. In the turbulent Kolmogorov velocity field, the smallest eddies have the highest shearing rate given by $`v(l)/ll^{2/3}`$, where $`l`$ is the size of the eddy, and $`v(l)l^{1/3}`$ is its turn-over velocity. Therefore, at this stage the magnetic field grows predominantly on small scales. Fig. 1 shows schematically the initial stage of the evolution of the magnetic spectrum, see, e.g., .

The magnetic energy collapses toward small, resistive scale during this initial evolution, until the field is strong enough to affect the dynamics of fluid through the Lorentz force. Such behavior is the evidence of small-scale turbulent dynamo; it is firmly established in numerical experiments, and can be derived analytically . Since in this example, the resistive scale is smaller than the viscous one, the dynamo is essentially governed by a smooth, viscous-scale velocity field.

To study the growth of magnetic energy on larger scales we have to understand how the magnetic field is generated in the inertial interval of the turbulence (the interval of scales much smaller than the external scale where the turbulence is excited, and much larger than the viscous scale where the turbulent energy is dissipated). This question is non-trivial since in this interval the velocity field is not smooth, i.e., $`v(l)l^\alpha `$, with $`\alpha <1`$. This situation is especially relevant for the case of small magnetic Prandtl number ($`\text{Pm}=\nu /\eta `$, where $`\nu `$ is fluid viscosity and $`\eta `$ is resistivity), where the magnetic energy is concentrated mostly in the inertial interval of the velocity field.

The question was first addressed by Batchelor  who used the analogy between magnetic field lines and fluid vorticity lines to conclude that when $`\nu \eta `$, magnetic energy is not amplified. This analogy was criticized in  because the magnetic field can have arbitrary initial conditions, while the vorticity field is related to the velocity field. It was further argued in  that magnetic field line stretching generally dominates resistive reconnection, thus making dynamo possible for $`\nu \eta `$. Direct numerical simulations of MHD turbulence with $`\text{Pm}=1`$ confirm that week magnetic fluctuations are generally amplified by the Navier-Stokes velocity field except for a special case when the initial magnetic field is close to the vorticity field .

Recently, there appeared the number of numerical simulations of MHD turbulence with small magnetic Prandtl numbers, where magnetic fluctuations were not amplified , which revived the claims that dynamo does not exist in the Kolmogorov turbulence with $`\text{Pm}1`$, see e.g., . However, astrophysical observations show that magnetic fields are generated by turbulent motion rather effectively in planets and stars where magnetic Prandtl numbers are small (e.g., in the geo-dynamo, $`\text{Pm}10^510^6`$, in the solar photosphere, $`\text{Pm}10^7`$). This apparent contradiction motivated our interest to the problem.

In this paper we argue that dynamo always exists in a rough velocity field, such as the Kolmogorov field. We however find that the magnetic Reynolds number, $`\text{Rm}=v(L_0)L_0/\eta `$, and the numerical resolution required to generate magnetic field strongly depend on the roughness exponent of the velocity field, $`\alpha `$; the rougher the velocity field, the larger the required resolution. This result explains why dynamo action is hard to achieve in experiments with small magnetic Prandtl numbers, while it is easily achieved when the magnetic Prandtl number is large, $`\text{Pm}1`$. In the latter case the velocity field is effectively smooth at the resistive scales, where the magnetic energy is concentrated.

Our analysis is based on the so-called Kazantsev model of kinematic dynamo, which is exactly solvable and which allows us to change the velocity roughness exponent, $`\alpha `$.

2. The Kazantsev model for a rough velocity field. Consider the induction equation for the magnetic field:

$`B^i/t+v^jB^i/x^jB^jv^i/x^j=\eta \mathrm{\Delta }B^i,`$ (1)
where $`B^i(𝐱,t)`$ is the magnetic field, $`v^i(𝐱,t)`$ is the fluid velocity, $`\eta `$ is small ohmic resistivity, and we sum over the repeated vector indices. At the initial stage of evolution, the weak magnetic field is passively advected by the fluid. Therefore, one can use simplifying assumptions about the fluid velocity to make the problem manageable. Kazantsev  and Kraichnan  introduced the model based on the Gaussian, short-time correlated velocity field, with zero mean and the covariance

$`v^i(𝐱,t)v^j(𝐱^{},t^{})=\kappa ^{ij}(𝐱𝐱^{})\delta (tt^{}).`$ (2)
This model is a valuable tool for the analytical investigation of kinematic dynamo; direct numerical simulations reveal that a purely Gaussian velocity field amplifies small magnetic fluctuations in a similar manner as the true Navier-Stokes field .

Assuming isotropy and homogeneity, the velocity correlation function has the form

$`\kappa ^{ij}(𝐫)=\kappa _N(r)\left(\delta ^{ij}{\displaystyle \frac{r^ir^j}{r^2}}\right)+\kappa _L(r){\displaystyle \frac{r^ir^j}{r^2}},`$ (3)
where $`𝐫=𝐱𝐱^{}`$. If we further assume that the velocity field is incompressible, we have $`\kappa _N=\kappa _L+(r\kappa _L^{})/2`$, and velocity statistics can be characterized by only one scalar function, $`\kappa _L(r)`$.

The model defined by (1), (2), and (3) allows one to write a closed equation for the correlation function of the magnetic field,

$`B^i(𝐱,t)B^j(𝐱^{},t)=H^{ij}(𝐱𝐱^{},t),`$ (4)
where, analogously to (3), the $`H^{ij}`$ function can be represented as

$`H^{ij}=H_N(r,t)\left(\delta ^{ij}{\displaystyle \frac{r^ir^j}{r^2}}\right)+H_L(r,t){\displaystyle \frac{r^ir^j}{r^2}},`$ (5)
furthermore, the condition $`𝐁=0`$ gives $`H_N=H_L+(rH_L^{})/2`$. We will characterize the magnetic field correlator by the function $`H_L(r,t)`$. The equation for this function can be found by differentiating (4) with respect to time, and by using Eqs. (1)-(3). A rather tedious but essentially straightforward calculation gives:

$`_tH_L=\kappa H_L^{\prime \prime }+\left({\displaystyle \frac{4}{r}}\kappa +\kappa ^{}\right)H_L^{}+\left(\kappa ^{\prime \prime }+{\displaystyle \frac{4}{r}}\kappa ^{}\right)H_L,`$ (6)
where primes denote the derivatives with respect to $`r`$ and we have introduced the ‘renormalized’ velocity correlation function $`\kappa (r)=2\eta +\kappa _L(0)\kappa _L(r)`$. Equation (6) was originally derived by Kazantsev , and can be rewritten in a different form that formally coincides with the Shrödinger equation in the imaginary time. Effecting the change of variable, $`H_L=\psi (r,t)r^2\kappa (r)^{1/2}`$, one obtains:

$`_t\psi =\kappa (r)\psi ^{\prime \prime }V(r)\psi ,`$ (7)
which describes the wave function of a quantum particle with the variable mass, $`m(r)=1/[2\kappa (r)]`$, in a one-dimensional potential ($`r>0`$):

$`V(r)={\displaystyle \frac{2}{r^2}}\kappa (r){\displaystyle \frac{1}{2}}\kappa ^{\prime \prime }(r){\displaystyle \frac{2}{r}}\kappa ^{}(r){\displaystyle \frac{(\kappa ^{}(r))^2}{4\kappa (r)}}.`$ (8)
This equation can be investigated for different choices of $`\kappa (r)`$, however, we restrict ourselves to the inertial interval of the turbulence, where the velocity correlator has power-law asymptotics, $`\kappa (r)r^\beta `$. The exponent $`\beta `$ can be found from the scaling of turbulent diffusivity, $`Dv(r)rr^{1+\alpha }`$. Indeed, in the derivation of Eq. (6) we used the integral $`D=_0^{\mathrm{}}[v(x,t)v(x^{},t)][v(x,t^{})v(x^{},t^{})]d(tt^{})=\kappa (r)`$, which is the turbulent diffusivity . Comparing the two expressions we obtain $`\beta =1+\alpha `$.

The Shrödinger equation (7) has the effective potential $`U_{eff}(r)=V(r)/\kappa (r)=A(\beta )/r^2`$, where $`A(\beta )=23\beta /23\beta ^2/4`$. When $`A(\beta )<1/4`$, the quantum particle falls toward the origin  and its wave function is therefore concentrated at the smallest, resistive scale. This behavior is the manifestation of the dynamo mechanism that we discussed in the introduction. We obtain that $`A(1+\alpha )<1/4`$ for any roughness exponent of the velocity field, $`0<\alpha <1`$, therefore, dynamo is always possible; the same conclusion was made in .

3. Magnetic field correlator and dynamo growth rates. To find the wave function (the magnetic field correlator) we need to specify the boundary conditions. The small-scale regularization is naturally given by ohmic resistivity. For scales much smaller than the correlation scale of the velocity field, the $`\kappa `$ function can be expanded as $`\kappa (r)\kappa _L(0)(2\eta +r^{1+\alpha })`$, which corresponds to the limit of infinitely small Prandtl number. In this formula and below we use the dimensionless variables: $`\eta `$ is measured in the units of the large-scale turbulent diffusivity, $`\kappa _L(0)`$, while $`r`$ is in the units of the integral scale, $`L_0`$.

The boundary condition at the origin follows from finiteness of magnetic energy, $`\underset{r0}{lim}H_L(r,t)=H_0(t)<\mathrm{}`$. The boundary condition at large scales follows from the absence of the mean magnetic field, $`\underset{r\mathrm{}}{lim}H_L(r,t)=0`$. We will see in a moment that at the kinematic stage magnetic energy is concentrated at the resistive scale, and the corresponding wave function decays exponentially fast in the inertial interval, $`rr_\eta `$, so it is indeed independent of the large-scale properties of the velocity field, as we expected.

Problem (6) can be cast into the Sturm–Liouville form by changing $`H_L(r,t)=h(r,t)/r^2`$, which teaches us that its solution can be expanded in the eigenfunctions of the corresponding Sturm–Liouville operator. The maximum growth rate is therefore given by the largest eigenvalue of the operator on the right-hand side of Eqs. (6) or (7). The corresponding wave function is the ground state of the potential $`U_{eff}(r)`$. Thus, following Kazantsev , we look for the solution of (7) in the form $`\psi (r,t)=\phi (r)\mathrm{exp}(\lambda t)`$.

In the inertial interval, $`2\eta r^\beta `$ , the equation for the $`\phi `$ function reads

$`\phi ^{^{\prime \prime }}+\left[{\displaystyle \frac{3\beta ^26\beta +8}{4r^2}}+{\displaystyle \frac{\lambda }{r^\beta }}\right]\phi =0,`$ (9)
where $`\lambda `$ is dimensionless and is measured in the units of $`\kappa _L(0)/L_0^2`$. For small $`\lambda `$, the potential is dominated by the first term in the square brackets, and one can easily check that the corresponding wave function oscillates in the inertial interval. However, the ground-state wave function cannot oscillate, therefore the ground-state growth rate, $`\lambda `$, will be such that the second term in the square brackets of (9) dominate the first one in the inertial interval. At the resistive scale, $`r_\eta =(2\eta )^{1/\beta }`$, both terms should be of the same order, therefore, $`\lambda \eta ^{(\beta 2)/\beta }`$.

The wave function corresponding to the growing solution, $`\lambda >0`$, decays exponentially fast for $`r\eta ^{1/\beta }`$, $`\phi \mathrm{exp}\left(\frac{2\sqrt{\lambda }r^{(2\beta )/2}}{2\beta }\right)`$. This function is concentrated at the resistive scale and its growth rate is of the order of the eddy turn-over time at this scale; in agreement with our qualitative discussion in the introduction. When one changes the magnetic Prandtl number (by changing viscosity, for instance) the effective roughness of the velocity field at the resistive scale changes. We therefore suggest that the effect of different Prandtl numbers can be studied in terms of Eq. (9) with different roughness exponents. By this analogy, the limit of the smooth velocity field, $`\beta =2`$, corresponds to large Pm, while the other limit of the Kolmogorov-scaled velocity field, $`\beta =4/3`$, corresponds to small Pm.

We observe that the rougher the velocity field, the broader the wave function compared to the resistive scale, and, therefore, the larger the magnetic Reynolds number necessary to generate magnetic fluctuations. To estimate the critical resolution for the dynamo onset, we solved Eq. (6) numerically with the large-scale boundary condition, $`H_L(\mathrm{})=0`$. For given $`\beta `$ and $`\eta `$, we increased the ‘system size’, $`\mathrm{}`$, until $`H_L`$ started to grow and we thus found the critical value $`\mathrm{}=\mathrm{}_c`$. An analogous calculation could be done by fixing $`\mathrm{}`$ and changing $`\eta `$.

To characterize the inertial range, we introduce the resolution parameter, $`R=\mathrm{}_c/r_\eta =\mathrm{}_c/(2\eta )^{1/\beta }`$; this parameter is universal, i.e., it is independent of large and small-scale properties of the velocity field. We obtain that when the velocity roughness increases from $`\beta =2.0`$ to $`\beta =1.3`$, the corresponding resolution parameter increases from $`R3.8`$ to $`R29`$, see Fig. (2). This means that if we simulated the random equation (1) directly, the required numerical resolution should increase by about an order of magnitude in each spatial direction as we go from the smooth to the Kolmogorov-scaled velocity field.

The subresistive scales, $`r\eta ^{1/\beta }`$, can be easily described using Eq. (6). By direct substitution, one can check that the magnetic field correlator can be expanded as $`H_L(r,t)=H_0[1r^\beta /(2\eta )+\mathrm{}]\mathrm{exp}(\lambda t)`$, which implies that the spectrum of the magnetic field decays as $`H_kk^2|B_k|^2k^{1\beta }`$ in the subresistive region, $`\eta ^{1/\beta }k\nu ^{1/\beta }`$. With such a spectrum, the rate of magnetic energy dissipation at the viscous scale is bigger than the rate at the resistive scale by the factor $`k_\eta ^5|B_{k_\eta }|^2/k_\nu ^5|B_{k_\nu }|^2\text{Pm}^{(2\beta )/\beta }`$. Thus, the effective resistivity may be ‘anomalously large’ in turbulent flows with small Prandtl numbers.

Strictly speaking, the last result requires short-time correlated small-scale eddies and is not applicable in the Kolmogorov turbulence where velocity is correlated at the eddy turn-over time, $`\tau _{corr}(l)l/v(l)`$. At the subresistive scales, $`l<l_\eta `$, this correlation time is larger than the resistive relaxation time, $`\tau _\eta l^2/\eta `$, while in the Kazantsev model the velocity correlation time is smaller than the resistive relaxation time. However, model (2) is physically self-consistent. The large magnetic energy dissipation is balanced by the large energy transfer from the small-scale velocity eddies to the small scale magnetic field; this is possible since the fluid energy is formally infinite in (2). In a practical situation where the velocity correlation time is not infinitely small, the asymptotic tail $`H_kk^{1\beta }`$ will hold up to the scales where this correlation time is comparable to the resistive time, $`k1/\sqrt{\eta \tau _{corr}}`$. One can expect that for smaller scales, the magnetic energy spectrum will have a steeper decay, see, e.g, .

4. Conclusions. We proposed that magnetic fluctuations are always generated in isotropic and homogeneous three-dimensional turbulence if the magnetic Reynolds number is large enough. The required critical magnetic Reynolds number sharply increases with velocity roughness, which explains the ‘no-dynamo’ outcomes reported in numerical simulations with small Prandtl numbers, see e.g.,. We do not see any physical reason for the absence of dynamo in these numerical simulations rather than lack of resolution. Our results also suggest that obtaining small-scale dynamo will be a serious challenge for laboratory experiments, where the magnetic Prandtl number is small, $`\text{Pm}10^5`$, while the magnetic Reynolds number is rather moderate, $`\text{Rm}100`$ .

As we mentioned in the introduction, the Batchelor analogy of magnetic field lines and vorticity lines would hold for $`\nu =\eta `$ if the initial magnetic field were proportional to the vorticity field. In this special case, magnetic energy would not be amplified. In the model we investigate there are infinitely many such special initial conditions; any magnetic correlator $`H_L(r)`$ lacking the growing eigenfunctions in its expansion will not be amplified. This, of course, does not mean that dynamo does not exist.

We also note that the Schrödinger-type equation of the form (7) was considered for the magnetic field correlator in the case of a general velocity field in . In this sense the applicability of the equation may not be restricted to the Gaussian, short-time correlated velocity field only. For example, we also investigated the modified Kazantsev model, with a more realistic, finite-time correlated velocity field. The velocity correlation time at a given wave number, $`k`$, was chosen to be of the order of the eddy turnover time at the corresponding scale. This model could not be solved analytically, however, its numerical integration gave qualitatively the same results as (1,2). We will report these results elsewhere.

We would like to thank S. C. Cowley, R. Rosner, A. Schekochihin, and S. I. Vainshtein for many useful discussions. This work was supported by the NSF Center for magnetic self-organization in astrophysical and laboratory plasmas at the University of Chicago. S.B. would like to thank the Aspen Center for Physics, where this work was initiated.


text:
# Evidence for the Black Hole Event Horizon1footnote 11footnote 1George Darwin Lecture presented at the Royal Astronomical Society, London, December 13, 2002.

## 1 Introduction

The story of astrophysical black holes begins with a simple question: How does a star, or any other gravitating object, hold itself up against its own self-gravity? In the case of the Sun, the answer is easy. The Sun is hot because of thermonuclear reactions in its interior, and the resulting thermal pressure counteracts the compressive action of gravity. The same is true of all the stars we see shining in the night sky. When a star runs out of nuclear fuel and dies, it must find other ways to fight gravity. Dead stars with masses up to the Chandrasekhar limit, $`M_{\mathrm{Ch}}=1.4M_{}`$, become white dwarfs, where electron degeneracy supplies the necessary pressure. Above the Chandrasekhar limit, and up to a second mass limit, $`M_{\mathrm{NS},\mathrm{max}}23M_{}`$, dead stars become neutron stars, where neutron degeneracy pressure holds them up. But that, according to conventional physics, is the end of the road. If a dead star has a mass $`M>M_{\mathrm{NS},\mathrm{max}}`$, there is no known force that can hold the star up. What we have then is a black hole, one of the most extraordinary concepts in physics.

A black hole represents the ultimate victory of gravity, where all the mass in the object collapses down to a “singularity”, a true geometrical point (at least within classical physics). The object has no material surface. Instead, surrounding the singularity is an “event horizon”, which plays the role of a virtual surface. The event horizon is a one-way membrane through which matter and energy can fall in from the outside, but nothing, not even light, can escape from within. The region inside the horizon is thus causally cut off from the outside world. In a real sense, the horizon serves as an effective surface, even though there is no actual material there. For a non-spinning black hole of mass $`M`$, the radius of the horizon is given by Schwarzschild’s result, $`R_S=2.95(M/M_{})`$ km. Spinning black holes have a somewhat smaller radius for the same mass.

Since Nature almost certainly makes dead stars with $`M>M_{\mathrm{NS},\mathrm{max}}`$, black holes must exist in our Galaxy and in other galaxies. Finding these black holes has been a central goal of high energy astrophysicists for the last few decades, and indeed dozens of excellent candidates have been discovered, many in the last ten years or so. These discoveries are briefly summarized in §2. We then discuss in §§3-5 the main theme of this Lecture: How can we confirm that the black hole candidates discovered by astrophysics are truly black holes? Specifically, do the objects possess the defining characteristic of a black hole, namely an event horizon?

## 2 Astrophysical Black Holes

Many excellent black hole candidates have been discovered in a class of objects called X-ray binaries. These are double stars in which a compact primary star, either a neutron star or a black hole, accretes mass from a normal secondary companion star (see Fig. 1), and radiates the accretion luminosity in X-ray, ultraviolet and optical radiation.

A particular class of X-ray binaries, called soft X-ray transients or SXTs for short, has turned out to be especially helpful in the hunt for black holes. In these binaries, the mass accretion rate $`\dot{M}`$, and consequently the accretion luminosity $`L_{\mathrm{acc}}`$, vary with time. For most of the time, an SXT is in a very low luminosity state with $`L_{\mathrm{acc}}10^610^8L_{\mathrm{Edd}}`$, where $`L_{\mathrm{Edd}}=1.25\times 10^{38}(M/M_{})\mathrm{erg}\mathrm{s}^1`$ is the Eddington luminosity and $`M`$ is the mass of the accreting star. Every once in a while, however, the system goes into an accretion outburst and becomes very bright, with a luminosity almost equal to $`L_{\mathrm{Edd}}`$. It is during the accretion outburst that the binary is usually discovered. After the outburst, the luminosity slowly declines over a period of several months, and the binary reverts to quiescence until the next accretion outburst.

Whereas the optical light from an SXT during outburst comes mostly from the accretion disk, the emission during quiescence is dominated by the secondary star. Therefore, in quiescence, it is possible to observe the absorption lines in the stellar spectrum and use the Doppler shifts on these to measure the orbital motion of the secondary. From the orbital period $`P_{\mathrm{orb}}`$ of the binary and the semi-amplitude $`K_s`$ of the secondary’s line-of-sight velocity, one can then obtain the mass function

$$f(M)\frac{P_{\mathrm{orb}}K_s^3}{2\pi G}=\frac{\mathrm{sin}^3i}{(1+M_s/M_p)^2}M_p,$$
(1)
where $`M_s`$ and $`M_p`$ are the masses of the secondary and the compact primary, respectively, $`i`$ is the inclination angle of the binary orbit, and the final relation on the right is derived by applying Newton’s laws of motion to a binary system bound by gravity in a circular orbit. Regardless of the values of $`M_s`$ and $`i`$, one sees that $`f(M)`$ is a strict lower limit on $`M_p`$.

McClintock & Remillard (1986) showed that the SXT A0620–00 has a mass function close to $`3M_{}`$, above the likely maximum mass $`M_{\mathrm{NS},\mathrm{max}}`$ of a neutron star. A0620–00 was the first truly bona fide black hole candidate, and the sensational discovery of this object opened up a new avenue in the hunt for black holes. In the years since, a dozen or more additional black hole candidates have been discovered in SXTs, all identified as having mass functions $`f(M)M_{\mathrm{NS},\mathrm{max}}`$. Special mention should be made of the SXT V404 Cyg (Casares, Charles & Naylor 1992) which for many years was the champion among black hole candidates, with the largest known mass function of $`6.1M_{}`$. Clearly, the compact stars in A0620–00, V404 Cyg, and the other SXTs with large mass functions, are all excellent black hole candidates. In a few SXTs and other kinds of X-ray binaries, $`f(M)`$ itself is not very large, but there are independent estimates of $`i`$ and $`M_s`$ which, combined with the measured $`f(M)`$, indicate that $`M_p>M_{\mathrm{NS},\mathrm{max}}`$. These objects are also very good black hole candidates, though clearly not as good as the systems with large mass functions. Narayan, Garcia & McClintock (2002) list the best black hole SXTs as of 2001.

The typical black hole mass in SXTs is $`515M_{}`$. In addition to these stellar-mass black holes, supermassive black holes have been identified in the nuclei of galaxies. At the center of our own Milky Way Galaxy, for instance, the object Sagittarius A\* has been very convincingly shown to be dark and compact, and to have a mass $`4\times 10^6M_{}`$ (see the remarkable observations described in Schödel et al. 2002 and Ghez et al. 2003). Another excellent candidate is in the nucleus of the galaxy NGC 4258 (Miyoshi et al. 1995). A couple of dozen additional candidates have been confimed in the nuclei of other nearby galaxies (Gebhardt et al. 2003, and references therein), and today it is generally agreed that virtually every galaxy in the universe has a supermassive dark compact mass in its nucleus.

## 3 Is This Enough?

As the above summary indicates, astrophysicists have certainly discovered many compact stars that are too massive to be neutron stars. Can we therefore claim victory in the search for black holes? In the opinion of the speaker, it would be premature to do so.

It is true that physics, specifically General Relativity combined with our knowledge of the properties of matter up to nuclear density, tells us that a compact star with $`M>M_{\mathrm{NS},\mathrm{max}}`$ cannot be a neutron star and must therefore be a black hole. However, what makes a black hole unique — indeed the reason why astronomers spend enormous resources searching for it, and the public is so fascinated by it — is not that it is “not a neutron star” but that it possesses an event horizon. It is the event horizon that makes the black hole so special. Therefore, before claiming victory, it would be prudent to look for independent evidence that the black hole candidates identified so far on the basis of their masses actually do have event horizons. The speaker and his collaborators, Dr. Michael Garcia, Dr. Jeremy Heyl and Dr. Jeffrey McClintock, have searched for such evidence during the last several years.

What would constitute plausible evidence? In brief, we need an observed phenomenon — or lack of it — that is a unique signature of an event horizon. Ideally, we should compare black hole candidates to a control sample of objects, say neutron stars, that are known to have surfaces. We should show that some observable characteristic is distinctly different in the two classes of objects, and that the difference is consistent with the notion that one class (black hole candidates) has event horizons and the other class (neutron stars) has surfaces. Furthermore, the difference should not have any other plausible explanation. Of course, for this line of argument to work, we need a working definition of what is plausible and what is not. This tends to be somewhat subjective, though scientists by and large agree on the matter.

X-ray binaries are particularly good for such investigations since some X-ray binaries contain black hole candidates and some contain neutron stars. It is thus possible to find two well-matched sets of binaries, one with black hole candidates that we think may have event horizons, and the other with neutron stars that we know for sure have surfaces. The speaker and his collaborators pioneered such comparisons of matched black hole and neutron star X-ray binareis, and showed in an early study that quiescent black hole SXTs are very much dimmer than quiescent neutron star SXTs (Narayan, Garcia & McClintock 1997). The latest observations are summarized in Fig. 2, taken from McClintock et al. (2003). As explained in Garcia et al. (2001) and Narayan et al. (2002), the large luminosity difference — a factor of 100 to 1000 in Eddington units — is natural if black hole candidates have event horizons. Various alternative explanations have been advanced in an effort to get around this argument, but many of these explanations appear to be ruled out by subsequent observations (Narayan et al. 2002).

Following our initial work, other signatures have been identified that distinguish black hole and neutron star binaries. Sunyaev & Revnivtsev (2000) showed that neutron star systems have strong variability up to frequencies $`1`$ kHz, presumably from boundary layer radiation close to the stellar surface, whereas black hole systems have a significant decline in flux variations above $`1050`$ Hz, presumably because they have no surface and hence do not possess boundary layers. In another study, Done & Gierlinsky (2003) considered the luminosities and spectra of black hole and neutron star binaries and showed that there are large differences between the two classes. They suggested that the differences arise because black hole candidates have event horizons while neutron stars do not.

We have recently begun a detailed study of Type I X-ray bursts in X-ray binaries with a view to using the bursts to verify the presence of event horizons in black hole candidates. This is the topic of the rest of the Lecture.

## 4 Type I X-ray Bursts

When gas accretes on the surface of a neutron star, it is compressed by the strong surface gravity of the star. As the gas sinks under the weight of continued accretion, it becomes denser and hotter, until the conditions are right for igniting thermonuclear reactions. If the gas consists largely of hydrogen and helium as it usually does, the nuclear burning tends to be unstable (Hansen & van Horn 1975, see Lewin, Taam & van Paradijs 1993, Bildsten 1998, for reviews). The instability causes the accreted layer of gas to burn its nuclear fuel explosively within a very short time, leading to a burst of X-ray emission. After the fuel is consumed, the star reverts to its accretion phase, in which fuel accumulates on the surface, until the next thermonuclear instability is triggered. The star thus undergoes a semi-regular series of thermonuclear explosions.

Bursts of X-ray emission from X-ray binaries were first discovered by Grindlay et al. (1976) and were immediately identified with thermonuclear explosions on the surface of a neutron star. These explosions are known as Type I bursts (to distinguish them from a different kind of burst, called Type II bursts, which is not thermonuclear in origin). In a typical Type I burst, the luminosity of the neutron star increases to nearly the Eddington limit within less than a second, and the flux then declines over a period of seconds to tens of seconds. The time interval between bursts is usually several hours to perhaps a day or two. Note that Type I bursts are very distinct from the accretion outbursts described earlier. The latter are associated with changes in mass accretion rate and have durations of months, not seconds, and recurrence times of many years, not hours or days.

Type I bursts are very common and have been seen in many neutron star X-ray binaries. However, it is a remarkable fact that no black hole candidate in any X-ray binary has ever had a Type I burst. In some sense, it is obvious that these objects should not experience bursts. A Type I burst requires a surface where matter is compressed and heated until a thermonuclear instability is triggered. A black hole has no surface; matter simply falls in through the event horizon and disappears. Therefore, a black hole cannot have Type I bursts. Since black hole candidates are indeed observed not to have bursts, does it then prove that they have event horizons? The answer is, unfortunately, “No!”

An object that has a Type I burst must have a surface and therefore cannot be a black hole. This statement is uncontroversial. However, an object that does not have Type I bursts does not necessarily lack a surface and therefore is not necessarily a black hole. For instance, there are bona fide neutron star systems, such as most X-ray pulsars, that do not exhibit bursts even though they certainly have surfaces. In order to use the lack of bursts as evidence for the presence of event horizons, we need to carry out the following steps:

(i) Develop a detailed theory of thermonuclear stability of accreting gas on a compact star.

(ii) Use the theory to explain observations of Type I bursts in neutron star systems and, in particular, explain why some neutron stars have bursts and some do not.

(iii) Apply the theory to black hole candidates and demonstrate that these systems would certainly exhibit bursts if they had surfaces.

(iv) Show that the only plausible explanation for the lack of bursts in black hole candidates is that the objects have event horizons, and that all other explanations either are ruled out or are very implausible.

Until we successfully accomplish these four stages, we cannot claim that the lack of bursts in black hole binaries is strong evidence for the event horizon. Over the last two years, we have embarked on a systematic attack on this problem. As part of this work, we have developed a theoretical framework for understanding the stability of nuclear burning on compact stars (Narayan & Heyl 2002, 2003).

We follow a star as it accretes gas of specified nuclear composition (typically solar, corresponding to 70% by mass of hydrogen, 28% by mass of helium, and 2% by mass of heavier elements). As the accretion column grows, at each instant we solve for the quasi-equilibrium configuration of the gas layer and obtain the density, temperature and composition of the gas as a function of depth. We then carry out a linear stability analysis to determine if the layer is stable to small perturbations. If for some particular column depth the layer is unstable, we say that the system undergoes a Type I burst. From the amount of unburnt fuel available in the layer at the moment of instability, we estimate the net energy emitted in the burst, and from the column density and the mass accretion rate, we estimate the recurrence time between the bursts. If the accreted layer is stable for all column depths up to some very large value (say $`10^{13}\mathrm{g}\mathrm{cm}^2`$), then we say that the system is able to burn nuclear fuel stably and does not undergo bursts. The main innovation in our work is the application of a formal linear stability analysis. The approach is more rigorous than previous analyses, which employed heuristic prescriptions to determine the stability of the accreted gas.

### 4.1 Application to Neutron Star X-ray Binaries

Figure 3 summarizes the results obtained with this model for a neutron star of mass $`1.4M_{}`$ and various radii from $`10^{0.2}R_S=6.5`$ km to $`10^{0.6}R_S=16.4`$ km. The calculations correspond to a core temperature of $`T_{\mathrm{core}}=10^8`$ K, which is approximately the temperature expected for a core that cools via the modified URCA process of neutrino emission. This temperature is also close to the upper limit to $`T_{\mathrm{core}}`$ in neutron star SXTs (see Narayan & Heyl 2002, 2003). Figure 3 shows that, almost independent of radius, accreting neutron stars are unstable to bursts at low accretion rates $`L_{\mathrm{acc}}0.25L_{\mathrm{Edd}}`$ and stable at higher accretion rates. Observations indicate that bursts occur up to a luminosity of about $`0.250.3L_{\mathrm{Edd}}`$ (van Paradijs et al. 1988; Cornelisse et al. 2003; Tournear et al. 2003), and that systems brighter than this limit have a much reduced bursting frequency. The model is thus consistent with the observations, with a tendency perhaps to underestimate slightly the prevalence of burst behavior.

X-ray pulsars are generally found not to burst. Although these systems are usually not very luminous, nevertheless, because the accretion is channeled onto the magnetic poles by a strong magnetic field, the local mass accretion rate $`\dot{\mathrm{\Sigma }}`$ ($`\mathrm{g}\mathrm{cm}^2\mathrm{s}^1`$) is generally above the Eddington rate. Therefore, the lack of bursts is consistent with the predictions of the model. In the case of millisecond X-ray pulsars, however, the field is weak and the accretion flow is only mildly channeled. We would expect these sources to burst, and indeed they do.

One could compare model predictions with more detailed observations of burst recurrence times, burst durations, and other observables, as discussed in Narayan & Heyl (2003). We do not describe the results here, but merely note that the agreement overall is fairly good. In particular, some puzzling observations for bright systems with $`L_{\mathrm{acc}}>0.1L_{\mathrm{Edd}}`$, that previously had not been explained, appear to be predicted satisfactorily by the new model. For faint sources with luminosities below about $`0.03L_{\mathrm{Edd}}`$ there appear to be some discrepancies in the detailed model predictions, and for $`L_{\mathrm{acc}}<0.01L_{\mathrm{Edd}}`$ there are very few observations presently available to check predictions. This could be improved in the future. In any case, we believe we have successfully accomplished stages (i) and (ii) as listed in §4, and we feel that the model is sufficiently reliable to embark on stage (iii).

### 4.2 Application to Black Hole X-ray Binaries

Let us assume that black hole candidates in SXTs are not black holes, but are some unusual kind of compact objects with surfaces. Would these objects have Type I X-ray bursts? We can answer this question using our burst model.

First, we need to choose the core temperature and the stellar radius. From the very low luminosities of black hole SXTs in quiescence, we know that their core temperatures are no larger than $`10^{7.5}`$ K (Narayan & Heyl 2002). We present results corresponding to this core temperature, noting that the conclusions do not change very much for lower temperatures. We consider a range of likely radii for the compact mass, from $`(9/8)R_S`$ (the smallest allowed radius, see Shapiro & Teukolsky 1983) up to about $`2.8R_S`$.

Figure 4 shows the results for a $`10M_{}`$ object with a hypothetical surface. We see that such an object would have bursts over a wide range of mass accretion rate. That is, if SXT black hole candidates have surfaces, they ought to produce bursts. Consequently, the absence of bursts is inconsistent with the presence of a surface. This takes care of stage (iii) in §4 and goes a long way towards proving that black hole candidates must have event horizons. What is left to accomplish is the difficult stage (iv).

## 5 Why Do Black Hole Candidates Not Have Type I Bursts?

The lack of bursts in black hole candidates is certainly consistent with the presence of event horizons in these objects. This is not in dispute. But, is the presence of event horizons the only plausible reason for the lack of bursts? This is the critical question, and it is where the argument becomes most difficult and subjective. The following discussion covers all the possible explanations for the lack of bursts that the speaker considers plausible.

### 5.1 Obvious Explanations

Could black hole candidates be accreting the wrong kind of nuclear fuel, e.g., something other than hydrogen and helium? This is ruled out because the companion stars that supply the gas are very normal, with standard nuclear abundances. Moreover, a strong hydrogen line is seen in the optical spectra of all black hole SXT accretion disks, confirming the presence of hydrogen in the accreting gas, and helium lines are seen as well in many systems (Dr. J. E. McClintock, private communication).

Could the core temperature be such that bursts are quenched? Our model calculations indicate that bursts are quenched only if the core temperature is above several times $`10^8`$ K (Narayan & Heyl 2002, 2003). The very low X-ray luminosities of quiescent black hole SXTs indicate that their cores are much cooler, $`T_{\mathrm{core}}10^{7.5}`$ K. So this is not a promising explanation.

Could the burst recurrence times be extremely long, thereby making it very difficult to detect bursts? Figure 5 shows the dependence of the recurrence time on the mass accretion rate for a typical black hole candidate with a radius equal to $`2R_S`$. The recurrence times are a factor of a few longer than for neutron stars, as one might expect since the accreting objects have larger masses, but the bursts are correspondingly also more powerful by a factor of a few. There is no reason why these bursts should be missed by observations. Figure 5 indicates that the bursts are most easily observed for $`L_{\mathrm{acc}}0.1L_{\mathrm{Edd}}`$. For lower accretion rates, the systems are still burst-unstable (see Fig. 4), but the recurrence times are very long and it is possible to miss these bursts.

### 5.2 Wrong Accretion Rate?

Figures 4 and 5 show that for certain accretion rates very close to the Eddington rate bursts are either absent altogether or are very rare. Could it be that our black hole candidate systems spend all their time in these particular regions of parameter space and that this is why they do not burst? In other words, is it just plain bad luck that we have not seen bursts from these objects? This explanation can be eliminated easily.

Recall from §2 that all our black hole candidates are found in transient binary systems, so that each of them has experienced one or more accretion outbursts during which its mass accretion rate varied over a wide range. At the start of an accretion outburst, a typical system has a very low accretion rate, off beyond the left edge of Fig. 4. Then, as the outburst builds up, the accretion rate increases and the system moves to the right almost to the right edge of the figure, near the limit of Eddington accretion. This journey typically takes about a week. Then, the accretion rate slowly declines over a period of six months or so, during which the system moves slowly back towards the left edge of the figure. Finally, the system plunges back into quiescence, beyond the left edge of the figure.

Figure 5 indicates that there is a favorable range of mass accretion rate, between about 0.1 and 0.6 Eddington, for which bursts are expected to occur quite frequently. Although the particular calculation shown corresponds to a specific choice of the radius, this favorable region of bursting activity is present for all radii that we have tried. Furthermore, most black hole candidates have crossed this region during their accretion outbursts, and each system has spent about a month or so in the region. During this time it should have had several Type I bursts according to the calculations. Therefore, the fact that not a single source has experienced any Type I bursts is very significant, and cannot be “just plain bad luck.”

### 5.3 Rotation and Magnetic Fields

Could rapid rotation somehow eliminate bursts? Rotation has the effect of introducing a variation in the effective surface gravitational acceleration $`g`$ as a function of latitude. However, even a maximally rotating compact star has only a factor of $`2`$ variation in $`g`$ between the equator and the pole. Such a modest variation does not have a serious effect on the burst instability.

Rotation might influence the propagation of the burning front once the instability has been triggered at a single point, and might thereby have an effect on the development of the burst. This is a complex problem since the propagation of deflagration fronts is poorly understood. Although current calculations (e.g., Spitkowsky, Levin & Ushomirsky 2002), when applied to black hole candidates, seem to indicate that bursts ought to develop readily in these systems, it must still be considered an open question.

Neutron stars with strong magnetic fields do not burst. The accretion flows in these systems are channeled by the field, thereby increasing the effective local mass accretion rate $`\dot{\mathrm{\Sigma }}`$ to above the Eddington rate. Could black hole candidates have similarly strong fields and thereby avoid bursts? The possibility can be ruled out immediately. When accretion is channeled, the radiation comes out in beams from the magnetic poles, causing the observed X-ray flux to be modulated at the rotation period of the star. This is, in fact, the explanation for the modulations seen in X-ray pulsars. If black hole candidates have strong magnetic fields, they ought to exhibit X-ray pulsations. Since none of these objects has ever shown coherent pulsations, we can rule out the magnetically channeled accretion scenario. In principle, if the magnetic axis is aligned perfectly with the rotation axis there would be no observable pulsation, but it requires such an unusual conspiracy in so many objects that we consider it implausible.

### 5.4 Exotic Stars

Conventional physics tells us that black hole candidates cannot be normal compact stars such as white dwarfs or neutron stars. It is therefore reasonable to suppose that, if they are not black holes, then they must be exotic stars of some kind. Could they be exotic in such a manner as to prevent Type I bursts when they accrete gas?

The nature of matter deep in the core of a neutron star is uncertain (Shapiro & Teukolsky 1983; Glendenning 2000). It might be normal baryonic matter (mostly neutrons plus some protons and electrons), or it might be a pion condensate, or it might be something more unusual like Q-matter (Bahcall, Lynn & Selipsky 1990), quark matter (Glendenning 2000) or a color-locked superconducting state (Alford, Rajagopal & Wilczek 1998; Rapp et al. 1998). Could black hole candidates have cores made of some very exotic material, and could this prevent bursts? No, because bursts are very much a surface phenomenon, restricted to low densities $`10^6\mathrm{g}\mathrm{cm}^3`$. The unusual kinds of matter that are invoked for neutron star interiors typically occur at very high pressure, when the density is above nuclear density ($`10^{15}\mathrm{g}\mathrm{cm}^3`$) and certainly above the neutron drip point ($`10^{11}\mathrm{g}\mathrm{cm}^3`$, despite the claim to the contrary by Abramowicz, Kluzniak & Lasota 2002 who seem to state that the color-locked superconducting state can exist even at low densities). Changes at such high densities and such large depths have no effect on bursts at the surface.

Could black hole candidates be strange stars or some other kind of exotic stars where exotic matter extends all the way to zero pressure at the surface? Models of such stars have indeed been proposed (Bodmer 1971; Witten 1984; see Glendenning 2000). However, when the stars accrete, it is believed that the newly added gas is supported electrostatically as a separate surface layer of normal matter, getting converted to exotic matter only at high density. The low-density layer would behave exactly as in our models and would produce bursts just as in the calculations.

Is there any form of exotic matter that instantly converts infalling baryonic matter into exotic matter even at densities below $`10^6\mathrm{g}\mathrm{cm}^3`$? If a star were to be made of such material, there would be no bursts since there would be no nuclei on the surface to undergo thermonuclear burning. To the knowledge of the speaker, no such matter is presently known.

Finally, the black hole candidate may consist of some kind of dark matter that does not interact with baryonic gas other than via gravity (just like the dark matter in the universe). We could then have a $`10M_{}`$ compact object that is completely “porous” to ordinary matter. Infalling gas would fall freely through the dark matter and collect at the center. (The speaker thanks Professor M. J. Rees for suggesting this possibility.) Would such an object have Type I bursts? To answer this question, we need to calculate the radius, surface gravity and redshift of the baryonic gas ball that forms at the center and then check whether or not the object would have Type I bursts. We have carried out some preliminary calculations which suggest that the central gas sphere has characteristics not very different from a neutron star and that the accreting gas does undergo thermonuclear instability and Type I bursts.

## 6 Conclusion

The lack of Type I bursts in black hole SXTs is an important clue to the nature of these compact stars. Based on detailed models we conclude that, if these objects possess surfaces, they should exhibit widespread burst activity. Why then do we not see Type I bursts in black hole SXTs?

We have considered in §5 a variety of explanations for the lack of bursts and find that most can either be ruled out or can be ignored on grounds of implausibility. A couple of explanations are still viable. One is that, once a thermonuclear instability is triggered, the burning front, for some mysterious reason, is unable to propagate quickly on the surface of a black hole candidate (§5.3). This explanation, however, involves a considerable degree of special pleading. Burning fronts move quite rapidly on the surfaces of neutron stars and white dwarfs, so why should they have trouble on black hole candidates whose surface gravities lie in between? The other possibility is that black hole candidates are made of non-interacting dark matter so that there is no hard surface where the accreting gas can accumulate and undergo a thermonuclear instability (§5.4). However, in this case, the gas would still accumulate at the center and form a baryonic ball with a surface. We have carried out preliminary calculations of the bursting behavior of such objects and our results suggest that bursts should be common. Detailed results are awaited.

Leaving aside the possibilities mentioned above, by far the most plausible explanation for the lack of bursts in black hole candidates is that the objects simply have no surfaces. They must then have event horizons. Though not yet strong enough to qualify as proof, this line of reasoning must surely be considered compelling evidence for the reality of the event horizon.

This work was supported in part by NASA grant NAG 5-10780 and NSF grant AST 0307433.


text:
# Heavy quark masses from finite volume effectsTalk given at Lattice 2003

## 1 Introduction

Quark masses are fundamental parameters of the QCD Lagrangian. Their accurate knowledge is required in order to give quantitative predictions of fundamental processes. A direct experimental measurement of quark masses is not possible because of confinement, and their determination can only be inferred from a theoretical understanding of the hadron phenomenology. The calculations can be numerically performed through different strategies, depending upon the quark flavor and the available computational facilities. Accurate non–perturbative measurements of the $`u`$, $`d`$, $`s`$ and $`c`$ quark masses have been obtained from a straightforward comparison of hadron spectroscopy and lattice QCD predictions.

The situation is different for the $`b`$ quark, with present computers capabilities. In principle, the $`b`$ quark mass could be extracted, similarly to lighter flavors, by looking at the heavy–light and heavy–heavy meson spectrum. However, heavy–light mesons are characterized by the presence of two different scales, i.e. $`\mathrm{\Lambda }_{QCD}`$, that sets the wavelengths of the light quark, and the heavy $`b`$–quark mass. Managing these two scales in a naive way would require a very large lattice ($`O(100^4)`$ points).

The task of determining the $`b`$–quark mass has been faced in literature by resorting to some approximations of the full theory, e.g. HQET on the lattice , lattice NRQCD , or QCD sum rules . A novel approach recently introduced is based on the non–perturbative renormalization of the static theory and its matching to QCD .

An alternative approach to the bottom quark physics, based on finite size scaling, has been proposed in a previous paper , where it has been applied to the heavy–light meson decay constants (see also ). The main advantages of this step scaling method (SSM) are that the entire computation is performed with the relativistic QCD Lagrangian and that the continuum limit can be taken, avoiding the unfeasible direct calculation. In order to implement the SSM, a finite size scheme is required, and we adopt the Schrödinger Functional (SF) as the most useful framework.

## 2 Step scaling functions and HQET

The SSM has been designed in order to deal with two scale problems in lattice QCD  and it has been discussed in detail in . The main assumption of the method is that the finite size effects affecting the meson masses, $`M_{\{P,V\}}`$ have a mild dependence upon variations of the heavy quark mass and are controllable from a numerical point of view. Finite size effects can be obtained from

$`\sigma _P(L,m_1,m_2)={\displaystyle \frac{M_P(m_1,m_2)|_{2L}}{M_P(m_1,m_2)|_L}}`$
$`\sigma _V(L,m_1,m_2)={\displaystyle \frac{M_V(m_1,m_2)|_{2L}}{M_V(m_1,m_2)|_L}}`$ (1)
for the pseudoscalar and vector meson masses respectively.
To validate the hypothesis of low sensitivity upon the high energy scale, we can make use of the HQET predictions on the heavy–light meson masses. In the infinite volume the pseudoscalar and vector meson masses have the following expansion in terms of the heavy–quark mass:

$$M_X(m_h,m_l)=m_h+\overline{\mathrm{\Lambda }}(m_l)+\frac{\alpha _X(m_l)}{m_h}+\mathrm{}$$
(2)
where $`X\{P,V\}`$. Assuming the contribution of the $`1/m_h^2`$ corrections to be negligible, at finite volume one has

$$M_X(m_h,m_l,L)=m_h+\overline{\mathrm{\Lambda }}_X(m_l,L)+\frac{\alpha _X(m_l,L)}{m_h}$$
(3)
where $`\overline{\mathrm{\Lambda }}_X(m_l,L)`$ depends upon the spin of the meson state because of the contamination of the excited states to the finite volume correlations. Using eqs. (1) and (3) we obtain the HQET predictions for the step scaling functions of the heavy-light meson masses

$$\sigma _X(L,m_h,m_l)=1+\frac{\sigma _X^{(0)}(m_l,L)}{m_h}+\frac{\sigma _X^{(1)}(m_l,L)}{m_h^2}$$
(4)
This result requires some considerations. First we want to stress that in the infinite heavy–quark mass limit the step scaling functions have to be exactly equal to one, $`\sigma _X(L,m_l,m_h\mathrm{})=1`$111We thank A. Kronfeld and R. Sommer for having pointed out this property of $`\sigma _X`$.. This represents a strong constraint for the fits of the heavy–quark mass dependence of the step scaling functions.
The second observation concerns the number of terms to be considered in eq. (4). At order $`O(1/m_h)`$ one has

$$\sigma _X^{(0)}(m_l,L)=\overline{\mathrm{\Lambda }}_X(m_l,2L)\overline{\mathrm{\Lambda }}_X(m_l,L)$$
(5)
corresponding to the static approximation in eq. (3). By increasing the physical volume $`L`$, the difference between $`\overline{\mathrm{\Lambda }}_X(m_l,2L)`$ and $`\overline{\mathrm{\Lambda }}_X(m_l,L)`$ decreases because the two quantities have to be equal in the infinite volume limit, making the heavy–quark mass expansion of the finite volume effects rapidly convergent. The same arguments apply to the coefficient $`\sigma _X^{(1)}(m_l,L)`$ that has to be considered when in the expansion of the meson masses, eq. (3), the order $`O(1/m_h)`$ is taken into account. In our calculation we have performed the fits of the step scaling functions considering the $`O(1/m_h^2)`$ term, $`\sigma _X^{(1)}`$, beyond the so called static evolution (SE) that retains $`\sigma _X^{(0)}`$ only. We will also report for comparison the fits for the SE that are anyway compatible.

## 3 Results and Discussion

The physical numbers for the heavy meson spectrum have been obtained combining the results of a small volume ($`L_0=0.4`$ fm) calculation, where the simulations have been performed at the physical values of the heavy quark masses, with the results of the step scaling functions according to the identity:

$$M_X(L_{\mathrm{}})=M_X(L_0)\sigma _X(L_0)\sigma _X(2L_0)\mathrm{}$$
(6)
The step scaling functions at the values of the heavy quark masses simulated on the small volume have been obtained by interpolation in the heavy–light case using eq. (4) while, in the heavy–heavy case, the results are linearly extrapolated. The plots of $`\sigma _P(L_0)`$ and $`\sigma _P(2L_0)`$ as functions of the inverse quark mass are shown in figs. 1 and 2 respectively.

Comparing the numerical heavy meson spectrum with the experimental results  we have obtained different determinations of the $`b`$–quark mass, depending upon the physical state used as experimental input.
Within the quenched approximation, the determinations of the quark masses coming from the heavy–heavy or from the heavy–light spectrum in principle differ because the theory does not account for the fermion loops. The numbers we quote as final results are obtained by averaging the results for the heavy-heavy and the heavy-strange vector and pseudoscalar mesons and by keeping the typical error of a single case:

$`m_b^{RGI}`$ $`=`$ $`6.73(16)\text{ GeV}`$
$`m_b^{\overline{MS}}(m_b^{\overline{MS}})`$ $`=`$ $`4.33(10)\text{ GeV}`$ (7)
for the $`b`$–quark and

$`m_c^{RGI}`$ $`=`$ $`1.681(36)\text{ GeV}`$
$`m_c^{\overline{MS}}(m_c^{\overline{MS}})`$ $`=`$ $`1.319(28)\text{ GeV}`$ (8)
for the charm. The latter results compare favorably with the results of the direct computations .

Our error estimate includes both the *statistical* error from the Monte Carlo simulation as well as the *systematic* error coming from the uncertainty on the lattice spacing corresponding at a given $`\beta `$ value (we have used $`r_0=0.5`$ fm ) and from the uncertainty on the renormalization constants. The final errors on the continuum quantities, of the order of $`2\%`$ percent for the renormalization constants and of about $`1\%`$ percent for the scale, are added in quadrature and then linearly added to the statistical errors. The evolution to the $`\overline{MS}`$ scheme has been done using four–loop renormalization group equations.


text:
# Dynamic density functional study of driven colloidal particles: the effect of the system dimension

The dynamical properties of classical fluids at pico-liter scale attract experimentally and theoretically much attention in the soft-matter and biophysics communities, due to the appearance of the microfluidics ref.mic , also called ’lab-on-a-chip’ ref.labchip , as one of the most relevance $`21`$th century technologies. Theses devices try to put a full chemical laboratory onto your hand thought a sophisticated dispositive that includes particles pump ref.pp  and gates, chemical reaction ref.Yeomans , etc.. In this work we focus our attention on one of the key factors for the design and construction of this type of devices: the system dimension. We consider a generic system formed by a dilute solution of colloidal particles dragged ( at a constant rate c in the $`\widehat{𝐳}`$ direction) by a moving potential barrier, modeled by a time dependent external potential ( $`V_{ext}(𝐫,t)=V_{ext}(𝐫ct\widehat{𝐳})`$) acting on the colloidal particles but with no effect on the solvent. We base our results on a new technique named Dynamic density functional (DDF) theory ref.ddf ; ddf2  which is a generalization of Density functional (DF) theory to out of equilibrium, for systems with relaxative molecule dynamics (those in which hydrodynamic modes are not relevant), under the hypothesis that the dynamic correlations can be approximated by those in a system with the same density distribution at equilibrium; the DDF is a promising way of describe classical fluids dynamics at a microscopic level including our knowledge of the very well developed DF theory for the equilibrium.

The central deterministic DDF equation is,

$`{\displaystyle \frac{\rho (𝐫,t)}{t}}=\mathrm{\Gamma }_o\left[\rho (𝐫,t)\left({\displaystyle \frac{\delta [\rho ]}{\delta \rho (𝐫,t)}}+V_{\mathrm{ext}}(𝐫,t)\right)\right];`$ (1)
where $``$ is the Helmholtz free energy in DF theory, usually split onto the ideal gas contribution and the excess due to the interactions ($`\mathrm{\Delta }`$), and $`\mathrm{\Gamma }_o`$ is the mobility of an isolated colloidal particle in the solvent.
In the case of the shifting external potential we are dealing with, $`V_{\mathrm{ext}}(𝐫ct\widehat{𝐳})V_{\mathrm{ext}}(𝐫^{})`$ (where $`𝐫^{}`$ is the coordinate in the reference framework of the external potential), the most relevant feature comes from the stationary state (achieved after a transient period from any initial state ) in which $`\rho (𝐫,t)=\rho (𝐫ct\widehat{𝐳})\rho (𝐫^{})`$ i.e. when the density distribution is shifted at the same rate c as the potential barrier, so that, eq.(1) reduces to,

$`\left[\rho (𝐫^{})\left({\displaystyle \frac{\delta [\rho ]}{\delta \rho (𝐫^{})}}+V_{\mathrm{ext}}(𝐫^{})+{\displaystyle \frac{cz^{}}{\mathrm{\Gamma }_0}}\right)\right]=0.`$ (2)
For non interacting particles, such that $`\mathrm{\Delta }=0`$, eq.(2) becomes a linear Fokker-Planck equation $`^2\rho +(\rho \beta V_k)=0`$, with a $`\mathrm{"}`$kinetic potential$`\mathrm{"}`$ given by $`\beta V_k(𝐫^{})=\beta V_{\mathrm{ext}}(𝐫^{})+\overline{c}z^{}`$, as a function of the reduced shifting rate $`\overline{c}\beta c/\mathrm{\Gamma }_o`$, with inverse length units and $`\beta =\frac{1}{k_bT}`$. In this case, the effective dimension behaviour is completely determined by the difference between the system dimension (container) and the dimensions over which the external potential is homogeneous, so that an effective $`1D`$ system got thought a $`2D`$ plane barrier dragging particles in a $`3D`$ container, is equivalent to a point-like barrier dragging particles in a narrow channel ref.PRLcha ; similarly, an effective $`2D`$ system may be observed for a cylindrical barrier in a $`3D`$ container, or with a point-like barrier in a $`2D`$ container. For $`\mathrm{\Delta }0`$ (interacting particles) ref.1d ; joe2 ; prlfrus  this is not longer true, because the excess free energy depends explicitly on the dimension of the colloidal particles container.

Since the system dynamic behaviour is determined mainly through the particles interaction and the effective system dimension, and we are interested here in the characteristics induce by this last factor, we consider in this work only non interacting particles which would represent diluted solutions of colloidal particles. We set the notation before we start solving eq.(2) for each dimension. In the following, we refer to the front (wake) region as that away from the external potential, well beyond (behind) the z’ positions with $`V_{ext}0`$; if we assume that this external potential barrier is located around the origin of the shifting reference framework, and restricted to values $`|z^{}|\sigma `$ with $`\sigma `$ the molecular diameter, then the points with $`z^{}\sigma `$ are in the front while those with $`z^{}\sigma `$ are in the wake. In both regions $`V_{ext}=0`$, eq.(2) reduces to $`^2\rho +\overline{c}\frac{\rho }{z^{}}=0`$, and the analytic form of $`\rho (𝐫^{})`$ may be obtained as a series with the appropriate asymptotic limit in each dimension. Notice that only for effective $`1D`$ systems, the front and the wake are fully separated by the external potential. In higher dimensions, they are connected through those regions with $`|z^{}|\sigma `$ but $`|𝐫^{}|\sigma `$, for which $`V_{ext}=0`$. Our result are presented as follow: first we solve analytically eq.(2) (where $`V_{ext}=0`$), them we find numerically the full solution, and finally we compare both.

For a shifting repulsive barrier in effective $`1D`$ systems ref.1d , the solution of eq.(2) is given by $`\rho (z^{})=\rho _0+A\mathrm{exp}(\overline{c}z^{})`$ at the front, where $`A/\rho _0`$ is a constant, fixed by the nature of the barrier and the boundary conditions, and by $`\rho (z^{})=\rho _0`$, without any kind of structure, in the wake region. It may be proved by a formal reduction of the stationary DDF equation with $`\mathrm{\Delta }=0`$ in any dimension, that this exponential decay behaviour at the front, and this absence of excess molecules in the wake, are generic characteristics of the total excess density distribution, under constant bulk boundary conditions ($`lim_{|𝐫^{}|\mathrm{}}\rho =\rho _o`$), for the integrated excess density distribution across the transverse directions, i.e. $`\frac{\mathrm{\Delta }n(z^{})}{\rho _o}=\frac{1}{\sigma }𝑑x(\frac{\rho (𝐫^{})}{\rho _o}1)`$ in effective $`2D`$ systems, and $`\frac{\mathrm{\Delta }N(z^{})}{\rho _o}=\frac{2\pi }{\sigma ^2}DrR(\frac{\rho (𝐫^{})}{\rho _o}1)`$ in effective $`3D`$ systems.
In effective $`2D`$ systems, for a shifting cylindric external potential (e.g. the dielectric potential created by a laser beam through a bulk solution), eq.(2), for the region with $`V_{\mathrm{ext}}(𝐫^{})=0`$, reduces to,

$$\frac{^2\rho }{x^2}+\frac{^2\rho }{z^2}+\overline{c}\frac{\rho }{z^{}}=0.$$
(3)
Through a Fourier transform, the solution is given by,

$$\rho (x,z^{})=\rho _0+_0^{\mathrm{}}𝑑\alpha \mathrm{cos}[\alpha x]e^{\beta z^{}}g(\alpha );$$
(4)
where,

$$\beta _\pm =\frac{\overline{c}}{2}\pm \sqrt{\left(\frac{\overline{c}}{2}\right)^2+\alpha ^2}.$$
(5)
Far away from the external potential, a fast convergence of the Fourier components $`g(\alpha )`$ is observed, and the relevant features come from their behavior for $`\alpha \overline{c}`$. Hence we may use $`\beta _+\overline{c}+\frac{\alpha ^2}{\overline{c}}`$ and $`\beta _{}\frac{\alpha ^2}{\overline{c}}`$, for the front and wake regions respectively. The expansion of $`g_\pm (\alpha )`$ as an even polynomial function for small $`\alpha `$, and the zero wake requirement, lead in the wake to $`g_{}(\alpha )A_1\alpha ^2+A_2\alpha ^4+\mathrm{}`$, while at the front we expect $`g_+(\alpha )B_0+B_1\alpha ^2+B_2\alpha ^4+\mathrm{}`$ Thus,

$`\rho (x,z^{})\rho _0+\sqrt{\pi }e^{\overline{c}z^{}}e^{\frac{a^2}{4}}\left[B_0{\displaystyle \frac{w}{2}}+B_1{\displaystyle \frac{w^3}{8}}\left(a^22\right)+B_2{\displaystyle \frac{w^5}{32}}\left(a^412a^2+12\right)+\mathrm{}\right],`$ (6)
for $`z^{}\sigma `$, and

$`\rho (x,z^{})\rho _0+\sqrt{\pi }e^{\frac{a^2}{4}}\left[A_1{\displaystyle \frac{w^3}{8}}\left(a^22\right)+A_2{\displaystyle \frac{w^5}{32}}\left(a^412a^2+12\right)+\mathrm{}\right];`$ (7)
for $`z^{}\sigma `$; where $`w=\sqrt{\overline{c}/z^{}}`$ and$`a=wx`$.

In effective 3D system, for a shifting spherical external potential (e.g. a driven colloid particle perturbing a bulk solution of others colloids), we use cylindrical coordinates $`𝐫^{}=(R,\varphi ,z^{})`$, and eq.(2), for the region with $`V_{\mathrm{ext}}(𝐫^{})=0`$, reduces to,

$`{\displaystyle \frac{1}{R}}{\displaystyle \frac{}{R}}\left(R{\displaystyle \frac{\rho (𝐫^{})}{R}}\right)+{\displaystyle \frac{^2\rho (𝐫^{})}{z^2}}+\overline{c}{\displaystyle \frac{\rho (𝐫^{})}{z^{}}}=0.`$ (8)
Through a Hankel transform, the solution is given by ,

$`\rho (R,z^{})=\rho _0+{\displaystyle _0^{\mathrm{}}}𝑑\alpha \alpha f(\alpha )J_0[\alpha R]e^{\beta z^{}};`$ (9)
where $`J_0`$ is the zeroth order Bessel function. Proceeding with the Hankel components as in the effective $`2D`$ with the Fourier ones, we might take $`f_{}(\alpha )C_1\alpha ^2+C_2\alpha ^4+\mathrm{}`$ in the wake, and at the front $`f_+(\alpha )D_0+D_1\alpha ^2+D_2\alpha ^4+\mathrm{}`$. Thus,

$`\rho (R,z^{})\rho _0+e^{\overline{c}z^{}}e^{\frac{h^2}{4}}\left[D_0{\displaystyle \frac{w^2}{2}}+D_1{\displaystyle \frac{w^4}{8}}\left(h^24\right)+D_2{\displaystyle \frac{w^6}{32}}\left(h^416h^2+32\right)+\mathrm{}\right],`$ (10)
for $`z^{}\sigma `$, and

$`\rho (R,z^{})\rho _0+e^{\frac{h^2}{4}}\left[C_1{\displaystyle \frac{w^4}{8}}\left(h^24\right)+C_2{\displaystyle \frac{w^6}{32}}\left(h^416h^2+32\right)+\mathrm{}\right];`$ (11)
for $`z^{}\sigma `$; where $`h=wR`$.

The coefficients (A,B,C and D) are determined by the particular choice of the external potential but, for each dimension, the asymptotic decay forms are generic. In the transverse plane (for fixed $`z^{}`$) the excess stationary density distribution ( $`\rho (𝐫^{})\rho _o`$) is given by a gaussian times a polynomial function, and along the direction of the movement (z’ axis) for a fixed distance to the origin ($`a=const`$ in effective $`2D`$ systems and $`h=const`$ in the $`3D`$ ) is given, at the front, by an exponential decay ($`\mathrm{exp}(\overline{c}z^{})`$) and an inverse power in the wake. Over the transverse directions, at the front, the maximum density is at the $`z^{}`$ axis, while in the wake, this axis is a minimum between two symmetric maxima which produce a cap-like structure with parabolic shape. The distance between the maximum and the z’ axis is $`2\sqrt{2z^{}/\overline{c}}\sigma `$ in effective $`3D`$ systems and $`\sqrt{6z^{}/\overline{c}}\sigma `$ in the $`2D`$. In fig.1 we present the numerical solution for the stationary density distribution of colloidal particles under the influence of an external potential $`V_{ext}=V_oe^{|𝐫^{}/\sigma |^6}`$, with $`\beta V_o=10`$, shifting at the velocity $`\overline{c}\sigma =10`$, in a$`)`$ an effective $`3D`$ and b$`)`$ $`2D`$ system. The clear cap-like structure, extended further away from the external potential, is clearly observed.

We have mentioned before that only in effective $`1D`$ systems the front and wake are fully separated by the external potential, which implicate that as the solution of eq.(2), the front and wake should share an unique analytic expression in effective dimensions beyond $`1D`$. This unique expression produces, at the front (wake), an exponential growth of the $`\beta _+(\beta _{})`$ contribution, which only by the appropriate behaviour of $`g(\alpha )(f(\alpha ))`$ for large $`\alpha `$ be canceled. We have compared our analytical prediction with the numerical solution of eq.(2), the concordance reveal a local asymptotic convergence which makes useful those asymptotic expressions (eqs.(6),(7),(10),(11)), to represent separately (at the front and in the wake regions) the stationary state. The use of polar (spherical) coordinates, instead of the rectangular (cylindric) ones, would not solve the problem, because the parabolic structure of the wake implies the entanglement of the radial and the angular coordinates.
Because this work may be usefully in the construction of particles pumps, we want to point out some effects to be considered in the design of an efficient device to drag particles. For the same radial form of the external potential it is clear that an effective $`2D`$ system is more efficient dragging particles than a $`3D`$, because it produces a higher density at the front peak. In fig.2 we present the numerical stationary density distribution of an effective $`3D`$ (dashed and dotted lines) and $`2D`$ system (solid and dotted-dashed lines), for different fixed values of x and R. The higher front density peak in the effective $`2D`$ system is clearly observed. For the same reason, an effective $`1D`$ system is more efficient than the $`2D`$.

For an effective $`1D`$ system we obtained in our previous work ref.1d  a non-linear dependence of the total current ($`j=\overline{c}(\rho (𝐫^{})\rho _o)𝑑𝐫^{}`$) with the drift velocity. If we denoted by $`\overline{c}_{max}`$, the velocity at which the total current is maximum for a given external potential, the efficiency of a device as a particles pump increases with $`\overline{c}`$ until $`\overline{c}_{max}`$ is reached, and then diminishes for higher $`\overline{c}`$. As we have mentioned before, the integrated excess of the stationary density distribution across the transverse directions, $`\mathrm{\Delta }n`$ in effective $`2D`$ systems and $`\mathrm{\Delta }N`$ in the $`3D`$, satisfies the same stationary DDF equation as the $`1D`$ density distribution, so we may expect for the total current this non-linear behaviour too. To support numerically our analytical result in fig.3 we present the integrated excess of the stationary density distribution, a$`)`$ in an effective $`2D`$ and b$`)`$ $`3D`$ system, the exponential decay behaviour at the front and the absence of structure in the wake as for the $`1D`$ stationary density distribution is clearly observed. We have estimated that the maximum total currents are obtain for $`\sigma \overline{c}_{max}`$ $`0.5`$ in $`1D`$, $`12.41`$ in $`2D`$, and $`14.76`$ in $`3D`$, which implicate that, apart from the dependence of $`\overline{c}_{max}`$ with the external potential shape, it is a function of the effective system dimension. $`\overline{c}_{max}`$ determines when for a given external potential the system is saturated by the drift velocity. Increasing the dimension, we increase the paths for the molecules to avoid the external potential so a higher velocity is needed to make the molecules go thought the external potential.

We have considered the geometrical effects of the effective dimension in systems of non-interacting colloidal particles driven by external barriers with constant drift rates. The main difference between the 2D and 3D cases analyzed here and the previously reported 1D case is the appearance of long wake structures behind the moving barrier. In contrast with the common exponential decay $`\rho (𝐫^{})\rho _oexp(\overline{c}z^{})`$ for the structure at barrier front, the structure of the wake decays only as an inverse power, $`\rho (𝐫^{})\rho _o1/|z^{}|^2`$ for 3D and $`\rho (𝐫^{})\rho _o1/|z^{}|^{3/2}`$ for 2D. However, the transverse integral of the wake, over the directions other than $`z^{}`$, gives no excess of particles, with respect to the unperturbed density $`\rho _o`$ far from the potential barrier. This is in agreement with the complete lack of wake structure behind a 1D barrier, and it implies that the global current induced by the shifting potentials, calculated as the integral of $`\overline{c}(\rho (𝐫^{})\rho _o)`$ comes from the balance between the positive contribution from the density peak at the front, and the void ($`\rho (𝐫^{})<\rho _o`$) within the barrier. This total current follows a non-monotonic dependence with the shifting rates, as reported in 1D 1d , with an optimum shifting rate $`\overline{c}_{max}`$ which increases with the effective dimensions. The colloid-colloid interactions effects, which have being leaved out of this work, may be included within the general framework of the DDF 1d ; joe2 ; prlfrus  and may help to the design of efficient particle pump devices.

F. P gratefully acknowledges discussions with J.D. This work has been supported by the Dirección General de Investigación Científica (MCyT) under grant BMF2001-1679-C03-02, and a FPU grant AP2001-0074 from the MECD of Spain.