text:
# Non-factorizable effects in the 𝐵-𝐵̄ mixing

## I Introduction

The study of the oscillations in the system of neutral $`B`$ mesons provides important information on the pattern of the CP violation in the Standard model and its extentions (a detailed discussion can be found in ). The main uncertainty in the theoretical description of the process arises from the nonperturbative long-distance contributions to the mixing amplitude. In the $`B\overline{B}`$ mixing there are two kinds of such contributions: first, effects related to the presence of the $`B`$-mesons in the initial and final states, and, second, corrections to the weak 4-quark amplitude due to the soft gluon exchanges between quarks of the initial $`B`$ and the final $`\overline{B}`$ meson. By neglecting the second contribution the amplitude factorizes. In this case all nonperturbative effects connected with the meson formation are reduced to only one quantity - the leptonic decay constant $`f_B`$.

Soft-gluon corrections to the weak 4-quark amplitude give rise to non-factorizable contributions, such that the influence of the $`B`$-meson structure is no longer described by the decay constant only. The understanding of the actual size of the non-factorizable effects and their adequate description in the heavy quark systems is an important and challenging task.

The theoretical analysis of $`B\overline{B}`$ mixing was first performed within the QCD sum rules  and later using lattice QCD, see  and refs therein. The results from the sum rules are in all cases well compatible with factorization (vacuum saturation), within rather large errors. Recent lattice analyses definitely reported negative non-factorizable effects, around 5-10% for the central values at the scale $`\mu m_b`$. Still, the errors of the calculations have a similar size. Results from these approaches are shown in Table I.

In this letter we analyse corrections to factorization in the $`B\overline{B}`$ mixing due to soft gluon exchanges in the leading $`\alpha _s`$-order, assuming that the main effect of such exchanges can be described by the local gluon condensate . We show that within this approximation the parameter $`\mathrm{\Delta }B_B`$ which describes correction to factorization, is determined by the meson transition form factors at zero momentum transfer. It is important that $`\mathrm{\Delta }B_B`$ turns out to be negative, independently of the particular values of these form factors.

As a next step we calculate the relevant form factors. To this end we make use of the relativistic dispersion approach based on the constituent quark picture . We obtain spectral representations for the form factors in terms of the $`B`$-meson wave function and the propagator of the color-octet $`b\overline{q}`$-system.

Parameters of the model, such as the quark masses and the wave functions, are known from the analysis of the weak $`B`$-meson form factors  within the same dispersion approach. For numerical estimates we also need the propagator of the color-octet $`q\overline{b}`$-system in the nonperturbative region, so we discuss the possible form of this propagator. The lack of reliable information about its details in the nonperturbative region leads to the main uncertainties in our numerical estimates.

The paper is organised as follows: The next Section presents basic formulas for the $`B\overline{B}`$ mixing. In Section III we calculate the non-factorizable contribution to the mixing amplitude induced by the local gluon condensate and show that the correction is negative. In Section IV we calculate the $`B`$-meson transition form factors within the dispersion approach. Section V gives numerical estimates.

## II Effective Hamiltonian and the structure of the amplitude

The effective Hamiltonian which summarizes the perturbative QCD corrections for the $`B^0\overline{B}^0`$ mixing has the following form

$$H_{eff}^{\mathrm{\Delta }B=2}=\frac{G_F^2M_W^2}{\sqrt{2}}\left(V_{tb}^{}V_{td}\right)^2C(\mu )\left(\overline{b}O_\sigma d\right)\left(\overline{b}O_\sigma d\right)+h.c.$$
(1)
where $`G_F`$ is the Fermi constant, $`O_\sigma =\gamma _\sigma (1\gamma _5)`$, $`\gamma ^5=i\gamma ^0\gamma ^1\gamma ^2\gamma ^3`$. The parameter $`\mu `$ stands for the renormalization scale which separates the hard region taken into account perturbatively and the soft region. The Wilson coefficient $`C(\mu )`$ includes perturbative corrections above the scale $`\mu `$. The explicit expression for $`C(\mu )`$ can be taken from . The soft contributions are contained in the $`B`$-meson matrix elements of the operators in the effective Hamiltonian

$`A=\overline{B}^0\left|\left(\overline{b}O_\sigma d\right)\left(\overline{b}O_\sigma d\right)\right|B^0.`$ (2)
The diagrammatic representation for the amplitude $`A`$ of Eq. (2) is shown in Fig 1.

Using the language of the hadronic intermediate states one finds that the contribution of the hadronic vacuum to this amplitude gives

$$\overline{B}^0|H_{eff}^{\mathrm{\Delta }B=2}|B^0=\frac{8}{3}\frac{G_F^2M_W^2M_B^2}{\sqrt{2}}(V_{tb}^{}V_{td})^2C(\mu )f_B^2,$$
(3)
with $`f_B`$ the leptonic decay constant of the $`B`$ meson defined according to the relation $`0|\overline{b}\gamma _\mu \gamma _5d|B(p)=if_Bp_\mu `$. It is convenient to parametrize the full amplitude as follows

$$\overline{B}^0|H_{eff}^{\mathrm{\Delta }B=2}|B^0=\frac{8}{3}\frac{G_F^2M_W^2M_B^2}{\sqrt{2}}(V_{tb}^{}V_{td})^2C(\mu )f_B^2B_B,$$
(4)
such that the quantity $`B_B1\mathrm{\Delta }B_B`$ measures contributions of the non-vacuum intermediate hadronic states.

Using the language of quarks and gluons, the corrections to the amplitude of Fig. 1 due to soft gluon exchanges are obtained by inserting the (soft) gluons between the quark lines in Fig 1. Soft gluon exchanges between the quarks of the same loop lead to the $`\alpha _s`$-corrections either to the meson vertices or to the quark propagators. They only contribute to the leptonic decay constant $`f_B`$ but do not lead to non-factorizable effects.

The non-factorizable effects originate from the soft gluon exchanges between quarks of different loops. To lowest $`\alpha _s`$-order these effects are described by the 4 diagrams shown in Fig. 2.

We assume that the main effect of the soft gluon exchange can be described by the local gluon condensate. In the next section we demonstrate that in this approximation the $`\alpha _sGG`$-correction to the factorization is negative.

## III $`\mathrm{\Delta }B`$ in terms of the local gluon condensate

Let us consider the exchange of a soft gluon between different quark loops assuming the dominance of the local gluon condensate . In this case a typical graph of Fig 2 describing the soft-gluon contribution is reduced to the product of the three-point diagrams as shown in Fig 3.

The corresponding amplitude has the form

$`A^{(1)}=g^2{\displaystyle 𝑑q𝑑q^{}T_{\sigma \alpha }^{ab(l)}(p,q)T_{\sigma \alpha ^{}}^{ab^{}(r)}(p,q)\frac{dx}{(2\pi )^4}\frac{dx^{}}{(2\pi )^4}\mathrm{e}^{iqx+iq^{}x^{}}A_\alpha ^b(x)A_\alpha ^{}^b^{}(x^{})},`$ (5)
where $`T_{\sigma \alpha }^{ab(l)}(p,q)`$ with the superscript $`l`$ stands for the amplitude of Fig 3a with the $`B`$ in the initial state, and $`T_{\sigma \alpha ^{}}^{ab^{}(r)}(p,q^{})`$ with the superscript $`r`$ for the amplitude of Fig 3b with the $`\overline{B}`$ in the final state.

It is convenient to use the fixed-point gauge $`x_\mu A_\mu ^b(x)=0`$ . In this case the gluon potential reads

$`A_\alpha ^b(x)={\displaystyle \frac{1}{2}}x_\beta G_{\alpha \beta }^b+\mathrm{},`$ (6)
where dots stand for terms with higher order derivatives. The average over the hadronic vacuum is performed according to the relation

$`G_{\alpha \beta }^bG_{\alpha ^{}\beta ^{}}^b^{}={\displaystyle \frac{\delta ^{bb^{}}}{96}}(g_{\alpha \alpha ^{}}g_{\beta \beta ^{}}g_{\alpha \beta ^{}}g_{\alpha ^{}\beta })GG,`$ (7)
with $`\delta ^{aa}=N_c^21`$. Here $`GG`$ is the positive-valued gluon condensate $`\frac{\alpha _s}{\pi }GG=0.012\mathrm{GeV}^4`$ .

The amplitude of Fig 3 then takes the form

$`A^{(1)}`$ $`=`$ $`{\displaystyle \frac{g^2}{4}}{\displaystyle 𝑑q𝑑q^{}T_{\sigma \alpha }^{ab(l)}(p,q)\frac{}{q_\beta }\delta (q)T_{\sigma \alpha ^{}}^{ab^{}(r)}(p,q^{})\frac{}{q_\beta ^{}^{}}\delta (q^{})G_{\alpha \beta }^bG_{\alpha ^{}\beta ^{}}^b^{}}`$ (8)
$`=`$ $`{\displaystyle \frac{g^2}{4}}G_{\alpha \beta }^bG_{\alpha ^{}\beta ^{}}^b^{}{\displaystyle \frac{\delta ^{ab}}{2\sqrt{N_c}}}T_{\sigma \alpha ,\beta }^{(l)}(p){\displaystyle \frac{\delta ^{ab^{}}}{2\sqrt{N_c}}}T_{\sigma \alpha ^{},\beta ^{}}^{(r)}(p),`$ (9)
where we have denoted

$`{\displaystyle \frac{}{q_\beta }}T_{\sigma \alpha }^{ab(l,r)}(p,q)|_{q=0}{\displaystyle \frac{\delta ^{ab}}{2\sqrt{N_c}}}T_{\sigma \alpha ,\beta }^{(l,r)}(p).`$ (10)
Notice that only the parts of the amplitudes $`T^{(l,r)}`$ antisymmetric in indices $`\alpha `$ and $`\beta `$ give a nonvanishing contribution. These antisymmetric parts have the following general Lorentz structure:

$`T_{\sigma \alpha ,\beta }^{A(l,q)}`$ $`=`$ $`4p^\nu \left[f_1^qϵ_{\alpha \beta \sigma \nu }+if_2^q(g_{\alpha \nu }g_{\beta \sigma }g_{\beta \nu }g_{\alpha \sigma })\right],`$ (11)
$`T_{\sigma \alpha ,\beta }^{A(l,\overline{q})}`$ $`=`$ $`4p^\nu \left[f_1^{\overline{q}}ϵ_{\alpha \beta \sigma \nu }+if_2^{\overline{q}}(g_{\alpha \nu }g_{\beta \sigma }g_{\beta \nu }g_{\alpha \sigma })\right],`$ (12)
$`T_{\sigma \alpha ,\beta }^{A(r,q)}`$ $`=`$ $`T_{\sigma \alpha ,\beta }^{A(l,\overline{q})},`$ (13)
$`T_{\sigma \alpha ,\beta }^{A(r,\overline{q})}`$ $`=`$ $`T_{\sigma \alpha ,\beta }^{A(l,q)},`$ (15)
$`f_1^q=f_1^{\overline{q}},f_2^q=f_2^{\overline{q}}.`$
The superscript $`q(\overline{q})`$ in these expressions denotes the flavour of the quark (antiquark) to which the soft external gluon is attached in diagrams of Fig 3. The quantities $`f_{1,2}^{q,\overline{q}}`$ are real constants. The last three lines in Eq. (11) are the consequence of the $`C`$-invariance of the strong interaction $`S`$-matrix. As we discuss in the next section, the constants $`f_{1,2}^{q,\overline{q}}`$ can be represented as specific $`B`$-meson transition form factors at zero momentum transfer.

We have to collect now contributions of the four diagrams in Fig 2. For instance, the contribution to the amplitude of the subprocess in which the soft gluons are attached to the $`b`$ quark in the left loop and the $`\overline{b}`$ quark in the right loop reads

$`T_{\sigma \alpha ;\beta }^{A(l,b)}(p)T_{\sigma \alpha ^{};\beta ^{}}^{A(r,\overline{b})}(p)(g_{\alpha \alpha ^{}}g_{\beta \beta ^{}}g_{\alpha \beta ^{}}g_{\alpha ^{}\beta })=192M_B^2\left((f_1^b)^2+(f_2^b)^2\right).`$ (16)
Similar expressions are easily obtained for other diagrams using the relations (11). Taking into account the relevant symmetry factors as shown in Fig 1, we finally arrive at the following expressions (we also list the factorizable amplitude in the same normalization)

$`A^{(0)}`$ $``$ $`{\displaystyle \frac{8}{3}}M_B^2f_B^2`$ (17)
$`A^{(1)}`$ $``$ $`8C_FM_B^2{\displaystyle \frac{\alpha _s}{\pi }}GG\pi ^2\left[(f_1^b+f_1^d)^2+(f_2^bf_2^d)^2\right],`$ (18)
with the color factor $`C_F=\frac{N_c^21}{4N_c}`$. Finally, for $`\mathrm{\Delta }B_B`$ we find the expression

$`\mathrm{\Delta }B_B(\mu )={\displaystyle \frac{\frac{\alpha _s}{\pi }GG}{f_B^2}}2\pi ^2\left[\left(f_1^b+f_1^d\right)^2+\left(f_2^bf_2^d\right)^2\right].`$ (19)
Obviously, the correction to the factorizable amplitude due to the local gluon condensate is negative. This agrees with all numerical results for $`\mathrm{\Delta }B`$. Notice that the constants $`f_1`$ and $`f_2`$ depend on the renormalization scale $`\mu `$.

## IV Form factors in the dispersion approach

We now proceed to the calcuation of the form factors $`f_1`$ and $`f_2`$. To this end we make use of the relativistic dispersion approach based on the constituent quark picture. Within this approach all observables are given by the spectral representations in terms of the soft wave function of the participating mesons.

We start with the amplitude $`T_{\sigma \alpha }^{(l)ab}(p,q)`$. It is diaginal in color indices, so we write $`T_{\sigma \alpha }^{(l)ab}(p,q)=\frac{\delta ^{ab}}{\sqrt{N_c}}T_{\sigma \alpha }^{(l)}(p,q)`$. The triangle diagram for $`T_{\sigma \alpha }^{(l)}(p,q)`$ is shown in Fig. 4a. Throughout this section we shall omit the superscript $`l`$. The quark structure of the $`B`$-meson is described by the vertex $`i\overline{q}_1\gamma _5q_2G(s)/\sqrt{N_c}`$, where $`s`$ is the invariant mass of the $`q_1\overline{q}_2`$ pair . One should set $`m_1=m_b`$, $`m_2=m_d`$ for the calculation of $`f^{(b)}`$, and $`m_1=m_d`$, $`m_2=m_b`$ for $`f^{(d)}`$.

We need to calculate $`T_{\sigma \alpha ,\beta }(p)=\frac{}{q^\beta }T_{\sigma \alpha }(p,q)|_{q=0}`$ (10). For the $`T_{\sigma \alpha }(p,q)`$ shown in Fig 4a, its derivative $`T_{\sigma \alpha ,\beta }(p)`$ is given by the diagram in Fig 4b. The trace corresponding to the Feynman diagram of Fig 4b reads:

$`S_{\sigma \alpha ;\beta }=\mathrm{Sp}\left[i\gamma _5(m_2\widehat{k}_2)O_\sigma (m_1+\widehat{k}_1)\gamma _\alpha (m_1+\widehat{k}_1)\gamma _\beta (m_1+\widehat{k}_1)\right],`$ (20)
where $`O_\sigma =\gamma _\sigma (1\gamma _5)`$. Since this expression is further multiplied by $`G_{\alpha \beta }`$, only its antisymmetric part in $`\alpha `$ and $`\beta `$ is necessary. The latter reads

$`S_{\sigma \alpha ;\beta }^A`$ $`=`$ $`(m_1^2k_1^2)\mathrm{Sp}\left(i\gamma _5(m_2\widehat{k}_2)\gamma _\sigma (1\gamma _5)\gamma _\beta (m_1\widehat{k}_1)\right)`$ (21)
$`=`$ $`(m_1^2k_1^2)\left[4ϵ_{\beta \alpha \sigma \nu }(m_2k_1^\nu +m_1k_2^\nu )+4im_1k_2^\nu (g_{\alpha \sigma }g_{\beta \nu }g_{\beta \sigma }g_{\alpha \nu })\right].`$ (22)
The corresponding expression for the Feynman amplitude takes the form

$`T_{\sigma \alpha ;\beta }^A(p)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^4i}}{\displaystyle 𝑑k_1𝑑k_2\delta (pk_1k_2)\frac{1}{(m_1^2k_1^2)^2(m_2^2k_2^2)}}`$ (24)
$`\times \left[4ϵ_{\beta \alpha \sigma \nu }(m_2k_1^\nu +m_1k_2^\nu )+4im_1k_2^\nu (g_{\alpha \sigma }g_{\beta \nu }g_{\beta \sigma }g_{\alpha \nu })\right].`$
We define $`f_v`$ and $`f_s`$ by the relations

$`{\displaystyle \frac{1}{(2\pi )^4i}}{\displaystyle 𝑑k_1𝑑k_2\delta (pk_1k_2)\frac{k_1^\nu }{(m_1^2k_1^2)^2(m_2^2k_2^2)}}`$ $`=`$ $`p^\nu f_v(p^2),`$ (25)
$`{\displaystyle \frac{1}{(2\pi )^4i}}{\displaystyle 𝑑k_1𝑑k_2\delta (pk_1k_2)\frac{1}{(m_1^2k_1^2)^2(m_2^2k_2^2)}}`$ $`=`$ $`f_s(p^2).`$ (26)
After the $`k`$-integration in Eq. (24) one recovers the structure of the amplitude from Eq. (11)

$`T_{\sigma \alpha ;\beta }^A(p)=4p^\nu \left[f_1ϵ_{\beta \alpha \sigma \nu }+if_2(g_{\alpha \sigma }g_{\beta \nu }g_{\beta \sigma }g_{\alpha \nu })\right],`$ (27)
with

$`f_1`$ $`=`$ $`(m_2m_1)f_v+m_1f_s,`$ (28)
$`f_2`$ $`=`$ $`m_1(f_sf_v).`$ (29)
Let us obtain now spectral representations for the form factors $`f_v^F`$ and $`f_s^F`$. We first notice that the expressions (25) correspond to the triangle Feynman diagram at zero momentum transfer $`q=0`$. This condition means that both $`q^2=0`$ and $`p^2=p^2`$, where $`p^{}=pq`$, so we can write $`f_i^F(p^2)=f_i^F(q^2=0,p^2,p^2=p^2)`$, with $`i=v,s`$.

It is convenient to start with the case $`q^20`$ and $`p^2p^2`$. Then the form factors $`f_i`$ can be written in the form of the double spectral representation :

$`f_i^F(q^2,p^2,p^2)={\displaystyle \frac{ds}{sp^2}\frac{ds^{}}{s^{}p^2}\mathrm{\Delta }_i(s,s^{},q^2)},`$ (30)
where the double spectral density $`\mathrm{\Delta }_i`$ can be calculated for any of the form factors from the Feynman integral.

As the next step, we set $`q^2=0`$, but still treat $`p^2`$ and $`p^2`$ as independent variables. In this case the double spectral representations of the form (30) simplify to the following single spectral representations

$`f_s^F(p^2,p^2)`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2}}{\displaystyle 𝑑s\frac{1}{sp^2}\frac{1}{sp^2}\mathrm{log}\left(\frac{s+m_1^2m_2^2+\lambda ^{1/2}(s,m_1^2,m_2^2)}{s+m_1^2m_2^2\lambda ^{1/2}(s,m_1^2,m_2^2)}\right)},`$ (31)
$`f_v^F(p^2,p^2)`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2}}{\displaystyle 𝑑s\frac{1}{sp^2}\frac{1}{sp^2}\frac{\lambda ^{1/2}(s,m_1^2,m_2^2)}{s}}.`$ (32)
The representations for the form factors in the form (32) are valid however only in the region of $`p^2`$ and $`p^2`$ (far) below the threshold. The reason for that is the following: the representations (30) and (32) are based on the Feynman form of the quark propagators. This form is however valid only for highly virtual particles, while in the soft region it is strongly distorted by nonperturbative effects. In particular, the pole at $`k^2=m^2`$ in the propagator of a color object, like quark and gluon, is absent. Recently, this was confirmed in a lattice study of the gluon propagator .

The modification of the quark propagator in the soft region leads to the change of the spectral representation for the form factors in the region of $`p^2`$ and $`p^2`$ near and above the $`b\overline{q}`$ threshold.

Notice that the quantities $`1/(sp^2)`$ and $`1/(s^{}p^2)`$ in Eq. (30) are the propagators of the initial and final $`b\overline{q}`$ states with virtualities $`s`$ and $`s^{}`$, and the squared masses $`p^2`$ and $`p^2`$, respectively. The nonperturbative effects which modify the quark and gluon propagators in the soft region, modify the propagators of the $`b\overline{q}`$ states as well.

We know the character of these changes in the $`p^2`$-channel corresponding to the $`B`$-meson: As discussed in , soft interaction of quarks effectively replace the factor $`1/(sM_B^2)`$ with a regular $`B`$-meson soft wave function $`\varphi _B(s)`$ in the integrand of (32). The $`B`$-meson soft wave function is normalized as follows

$`{\displaystyle \frac{1}{8\pi ^2}}{\displaystyle 𝑑s\varphi _B^2(s)\frac{\lambda ^{1/2}(s,m_b^2,m_q^2)}{2s}(s(m_bm_q)^2)}=1.`$ (33)
The leptonic decay constant is given by the expression

$`f_B={\displaystyle \frac{\sqrt{N_c}(m_b+m_q)}{8\pi ^2}}{\displaystyle 𝑑s\varphi _B(s)\frac{\lambda ^{1/2}(s,m_b^2,m_q^2)}{s}\frac{s(m_bm_q)^2}{s}}.`$ (34)
In the color-octet $`p^2`$-channel the Feynman propagator $`1/(sp^2)`$ in Eq (32) is replaced by the propagator $`D(s,p^2)`$ which involves proper modifications in the soft region. We therefore find the following representation for the form factors $`f_i(p^2)f_i(q^2=0,p^2=M_B^2,p^2)`$

$`f_v(p^2)`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2}}{\displaystyle 𝑑s\varphi _B(s)D(s,p^2)\frac{\lambda ^{1/2}(s,m_1^2,m_2^2)}{s}},`$ (35)
$`f_s(p^2)`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2}}{\displaystyle 𝑑s\varphi _B(s)D(s,p^2)\mathrm{log}\left(\frac{s+m_1^2m_2^2+\lambda ^{1/2}(s,m_1^2,m_2^2)}{s+m_1^2m_2^2\lambda ^{1/2}(s,m_1^2,m_2^2)}\right)}.`$ (36)
We know that $`D(s,p^2)1/(sp^2)`$ at large values of $`sp^2`$, and also that $`D(s,p^2)`$ is finite at $`s=p^2`$. However, we do not know details of $`D(s,p^2)`$ in the soft region. Motivated by the discussion in , we assume here that the nonperturbative effects can be described by the following modification of the propagator

$`D(s,p^2)={\displaystyle \frac{1}{sp^2+M_0^2}}`$ (37)
where $`M_0`$ is a mass parameter. In order to guarantee the absence of the pole in $`D(s,p^2)`$ in the heavy quark limit, $`M_0`$ should scale with the heavy quark mass as follows: $`M_0^2=O(\mathrm{\Lambda }_{QCD}m_Q)`$. So it is convenient to write $`M_0`$ in the form

$`M_0^2=wm_dm_b,`$ (38)
where $`m_d\mathrm{\Lambda }_{QCD}`$ is the constituent mass of the light quark, and $`w`$ is a parameter of order unity.*** Assuming that for the color-octet light $`q\overline{q}`$ system Eq (37) remains valid with $`M_0^2=wm_d^2`$, we find $`D(k^2=0)15\mathrm{GeV}^2`$ for $`w=1`$. This is close to the value of the gluon propagator $`D(k^2=0)18\mathrm{GeV}^2`$ as found in . This agreement seems to be reasonable since one can expect the propagator of the light $`q\overline{q}`$ color-octet system to have a structure in the nonperturbative region similar to the structure of the gluon propagator.

In the next section we make use of Eqs. (35) and (37) to analyse $`f_{v,s}(p^2)`$ for $`B_d`$ and $`B_s`$ mesons.

## V Numerical results.

The parametrization for the $`B`$-meson wave function and values of the quark masses have been determined from the analysis of the exclusive processes in . The $`B`$-meson soft wave function can be written in the form

$$\phi (s)=\frac{\pi }{\sqrt{2}}\frac{\sqrt{s^2(m_q^2m_{\overline{q}}^2)^2}}{\sqrt{s(m_qm_{\overline{q}})^2}}\frac{w(k^2)}{s^{3/4}}$$
(39)
with $`k^2=\lambda (s,m_q^2,m_{\overline{q}}^2)/4s`$. As found in  a simple Gaussian parametrization $`w(k^2)\mathrm{exp}(k^2/2\beta ^2)`$ gives a good description of the $`B`$-meson properties. The numerical parameters of the model from  are listed in Table II.

The procedure of Ref.  determines the wave function by fitting the lattice data on the weak transition form factors for large momentum transfers at the normalization point $`\mu 5`$ GeV. Therefore, the $`B`$-meson soft wave function and the form factors $`f_{s,v}(p^2)`$ at this scale are determined.

Whereas the $`B`$-meson wave function is known quite well, good information about the details of $`D(s,p^2)`$ is lacking. We therefore use the simple Ansatz (37) for $`D(s,p^2)`$ in the full range of $`s`$ and $`p^2`$ and treat $`w`$ as a free parameter of order unity. We expect that the variation of $`w`$ in the interval $`w=0.5÷2`$ provides reasonable error estimates for the form factors and $`\mathrm{\Delta }B`$. Notice that the value of $`\mathrm{\Delta }B`$ corresponding to $`w=1`$ agrees favourably with the recent lattice estimates, see Table I. Fig 5 shows the form factor $`f_s^{(d)}`$ which gives the main contribution to $`\mathrm{\Delta }B_B`$.

The calculated form factors are shown in Table III. The relative magnitudes of these form factors can be easily understood taking into account their scaling propertiers in the heavy quark limit: $`f_s^{(d)}m_b^{1/2}`$, $`f_s^{(b)}m_b^{3/2}`$, $`f_v^{(b)}/f_s^{(b)}1`$.

The value of $`\mathrm{\Delta }B`$ is obtained from Eq. (19) taking $`\frac{\alpha _s}{\pi }GG=0.012\mathrm{GeV}^4`$ from . Using the Wilson coefficent at the scale $`\mu =m_b`$ we obtain the renorm-invariant factors $`\widehat{B}_{B_d}`$ and $`\widehat{B}_{B_s}`$ (for definition see ) listed in Table I.

Our results are in good agreement with the lattice results, except for a slightly different estimate of the SU(3) violating effects in the $`B_d\overline{B}_d`$ and $`B_s\overline{B}_s`$ cases. This nice agreement gives support to our Ansatz for the description of the nonperturbative effects in the propagator of the color-octet $`q\overline{q}`$ system.

We would like to notice that taking the Ansatz (37) and (38) for the propagator of the $`s\overline{q}`$ system and setting $`M_0^2=wm_sm_{\overline{q}}`$, with $`w1`$ allows a parameter-free estimate of the $`K^0\overline{K}^0`$ mixing. In this case the result is much more stable with respect to the particular value of $`w`$ and therefore to the details of the propagator in the nonperturbative region. We obtained this way $`\mathrm{\Delta }B_K(1\mathrm{GeV})=0.21\pm 0.04`$ .

## VI Conclusion.

We considered corrections to factorization in the $`B\overline{B}`$ mixing amplitude due to soft-gluon exchanges, assuming that the main effect of such exchanges can be described in terms of the local gluon condensate.

1. It was demonstrated that within this approximation correction of the order $`\alpha _sGG`$ to the factorization is negative. It can be expressed through the specific $`B`$-meson transition form factors at zero momentum transfer.

2. A relativistic dispersion approach based on the constituent quark picture has been used to calculate these form factors. Spectral representations for the form factors in terms of the $`B`$-meson soft wave function and the propagator of the color-octet $`q\overline{b}`$ system were obtained. The behaviour of this propagator in the nonperturbative region was discussed. The obtained numerical estimates for $`\mathrm{\Delta }B_B`$ favourably compare with the recent lattice calculations.

Let us point out that the proposed approach can be extended to the description of more complicated problems. In particular, it can be applied to the analysis of the non-factorizable effects in non-leptonic $`B`$ decays.

###### Acknowledgements.

The authors are grateful to H. G. Dosch, O. Nachtmann, and B. Stech for discussions and interest in this work. D.M. was supported by the BMBF under project 05 HT 9 HVA3.


format: text

text:
# Numerical hydrodynamics: SPH vs. AMR

## 1. From Henyey to SPH and AMR

A history of numerical hydrodynamics in application to a problem of protostellar collapse begins more than three decades ago. It is remarkable to note that already at that time numerical simulations followed two distinct paths. Bodenheimer & Sweigart (1968) and Bodenheimer (1968) used a Lagrangian implicit method closely following Henyey’s approach (Henyey, Forbes & Gould 1964) commonly used for stellar evolution computations. On the other hand, Larson (1969) developed a variant of an Eulerian scheme for radiation hydrodynamics in the optically thick regime. Despite the use of very different numerical schemes, both studies yielded much similar results.

However, problems like the massive removal of angular momentum from the clump of molecular gas collapsing into stellar densities, could not be addressed within the one-dimensional framework. Among the simplest ways of removing excess angular momentum is the deposition of angular momentum in the orbital motion of planets formed around a central star, or by the formation of a multiple stellar system. Neither of the two ways of removal of angular momentum can be studied in detail in one dimension. To overcome this problem Lucy (1977) proposed a method in which the distribution of hydrodynamical quantities is approximated with a discrete set of diffuse clouds of matter, Smoothed Particle Hydrodynamics (Monaghan 1992). The SPH method is purely Lagrangian, can be easily implemented on a computer, and allows for performing 3-dimensional simulations even on relatively small systems.

On the other hand, the solution of the hydrodynamical equations in multidimensions poses severe problems for grid based codes. In the Lagrangian approach the grid cells follow individual fluid elements and whenever a significant amount of shear is present in the flow the numerical grid becomes heavily distorted. This problem is nicely illustrated by the work of Woodward (1976) in his study of a shock-cloud interaction problem. Since a strong distortion of the computational grid leads to a loss of accuracy of the discretization, most of the later research on Lagrangian schemes was devoted to the development of sophisticated rezoning modules which would allow to adequately follow individual fluid elements with an accuracy independent of their shapes.

An overall quality of the solution in the Eulerian model is determined by both the order of the advection scheme and the resolution provided by the numerical grid. Woodward and Colella (1984) demonstrated that substantial savings in terms of both memory and CPU time can be achieved when advection schemes of high-order are used instead of their low-order counterparts. These modern advection algorithms, shock-capturing schemes, gained much popularity during the last decade and follow the seminal works of Godunov (1959), van Leer (1976) and Harten (1983). The shock-capturing methods are characterized by at least a second order accuracy in the smooth part of the flow, with flow discontinuities being resolved in one or two zones, and no artificial viscosity is required to obtain a physically correct solution.

It took almost 20 years for Woodward’s pioneering study of a shock-cloud interaction to be successfully addressed with an Eulerian code (Klein, McKee & Colella 1994). However, even the use of a shock-capturing scheme was still not enough to adequately resolve the flow structure and to obtain a convergent solution: a novel technique of local adaptive mesh refinement (AMR; Berger & Colella 1989) was employed. In essence, AMR is a way of efficient discretization which allows to concentrate the computational effort in regions in which the errors of the solution are large. From a conceptual point of view AMR is similar to multigrid techniques with local refinements (Brandt 1984), in which the computational volume is covered with a hierarchy of completely nested grid patches and the resolution increases between parent and child grid levels by some predefined integer factor. However, due to the conservative character of the time-dependent hydrodynamical equations, the AMR method has to fulfill additional requirements making its actual computer implementation, to say the least, difficult.

## 2. Far side of numerical modelling

The planning of numerical simulations is much similar to the concept of project design known from industry rather than science. However, in both cases, given that objectives of the project have been clearly defined, a proper choice of the method of solution is necessary. Once an adequate tool is available one needs to find a way for using it. In practice this demand translates into having access to a computer installation which has to provide enough resources for obtaining the solution to the problem at hand. At first this last issue may appear to be a far side of numerical modelling, but to those involved in computationally intensive calculations this is that side which ultimately makes numerical simulation possible to accomplish.

Large simulations represent a separate class of computational problems which are characterized by a high value of the product of computer memory and processor time. This practical definition follows the observation that simulations done on grids with just a few thousand zones (or particles) might be demanding once the modelled system has to be evolved for a long time. A typical example here are simulations of accretion disks which suffer from short evolutionary timescales present near the inner edge of the disk. On the other hand, large memory is required for virtually all three-dimensional problems, even with memory savings indirectly offered by modern advection methods, or whenever very high, even if only locally, resolution is needed. The case of planet formation in protoplanetary disks serves here as an example.

Furthermore, the large size of the simulation demands a certain minimum speed of the processor and the main memory: program data has to be accessed and subsequently processed within acceptable time. Most of the large simulations put much more strain on the memory subsystem, including the system bus, than on the processor and, therefore, are “bandwidth” rather than CPU limited. Certainly, not only some of those available but even the largest existing computer installation may not satisfy the requirement imposed by a “memory times processing speed” limit: numerical experiments are limited by the availability of computer resources. Finally, also storage and analysis of large ($`10^2`$ times problem memory) data sets requires specialized hardware and software.

The variety of existing computer architectures does not necessarily make numerical simulations easier to plan and conduct. The performance of a computer code may differ vastly between scalar and vector architectures. Distributed and shared memory machines follow different programming paradigms, may ask for deep modifications to the existing codes. Such “technicalities” should be considered at the very beginning of any project which involves intensive computations as they might be difficult or impossible to overcome later. Keeping an eye wide open on current trends in computer architecture design, usually reflected in the policy of supercomputer centers, seems to be nowadays as important for numerical modelers as the development and application of new algorithms.

## 3. A few words of critics

In my introductory remarks I sketched a path that astrophysical hydrodynamics followed from relatively simple one-dimensional Lagrangian models to full three-dimensional simulations done with help of the SPH and AMR methods. I find this development remarkable even though, as we could learn during this conference, after three decades the problem of one-dimensional stellar collapse may still offer enough of material for debates. As for the comparison between SPH and AMR, the two so extremely distinct approaches to hydrodynamics, there appears to exist a sharp boundary dividing the groups of users of these two methods. However, both methods have their own limits and aspects which ask for special attention. I will list them briefly.

The Lagrangian nature of the SPH method makes spatial resolution impossible to maintain comparable in all regions of the computational volume, a worry for those who aim at studying properties of voids in cosmological simulations. Resolution near the flow discontinuities is usually poor, equivalent to a few smoothing lengths ($``$ particle “radii”), in most implementations of the SPH method shocks are handled with help of an artificial viscosity. Also the resolution and stability (in multidimensions) of contact discontinuities is problematic. SPH particles may “penetrate” each other increasing the diffusivity of the scheme and making studies which involve mixing and shear flows difficult.

Contrary to SPH, the resolution of an Eulerian model is defined by the smallest spatial scale and, would AMR not have been invented, no Eulerian calculation could achieve a resolution in mass comparable to that of large SPH simulations. However, the nonuniformity of the spatial resolution makes the AMR method unsuitable for studying turbulent flows which demand a smooth sampling of the spatial scales. Also, the roughness of the AMR grids adds a certain amount of vorticity to the model since each boundary between the levels in the AMR hierarchy appears as an obstacle to the flow and acts as a source of purely spurious vorticity.

Last but not least, the numerically “correct” results are those which correspond to the convergent solution. The Eulerian routine of doubling the resolution from one model to another and comparing low- and high-resolution models does not translate into doubling the number of particles in a SPH simulation. A concept of the smoothing length, usually mentioned when the spatial resolution of the SPH model is presented, is in this respect misleading since a twofold decrease of the smoothing length requires an increase in the number of neighbouring particles by a factor of $`2^D`$, where $`D`$ is the number of spatial dimensions. Yet, in most cases convergence studies (if present at all!) are done with the number of particles being varied by a factor of just a few.

I decided not to mention problems related to radiation transport or the incorporation of magnetic fields into hydrodynamical simulations. I did it on purpose since their complexity will likely add to the confusion; it is better to leave out their discussion for the years to come and closely follow the development of more efficient and correct algorithms. At this very moment it is my conviction that even pure hydrodynamics still poses challenging problems and obtaining physically meaningful results is, and as it has ever been, a difficult task.

### Acknowledgments.

I thank the organizing committee for financial support and hospitality during the conference. This work was partly supported by grant 2.P03D.014.19 from the Polish Committee for Scientific Research.

## References

Berger, M. J., & Colella, P. 1989, J. Comput. Phys., 82, 64

Bodenheimer, P. 1968, ApJ, 153, 483

Bodenheimer, P., & Sweigart, A. 1968, ApJ, 152, 515

Brandt, A., 1984, Multigrid Techniques: 1984 Guide, with Applications to Fluid Dynamics, GMD-Studien Nr. 85, Sankt Augustin, GMD

Godunov, S. K. 1959, Mat. Sb., 47, 217

Harten, A. 1983, J. Comput. Phys., 49, 357

Henyey, L. G., Forbes, J. E., & Gould, N. L., 1964, ApJ, 139, 306

Klein, R. I., McKee, C. F., & Colella, P. 1994, ApJ, 420, 213

Larson, R. B. 1969, MNRAS, 145, 271

Monaghan, J. J. 1992, ARA&A, 30, 543

van Leer, B. 1977, J. Comput. Phys., 23, 276

Woodward, P. R. 1976, ApJ, 207, 484

Woodward, P. R., Colella, P. 1984, J. Comput. Phys., 54, 115


text:
# Experimental Evidence of Time Delay Induced Death in Coupled Limit Cycle Oscillators

## Abstract

Experimental observations of time delay induced amplitude death in a pair of coupled nonlinear electronic circuits that are individually capable of exhibiting limit cycle oscillations are described. In particular, the existence of multiply connected death islands in the parameter space of the coupling strength and the time delay parameter for coupled identical oscillators is established. The existence of such regions was predicted earlier on theoretical grounds in \[Phys. Rev. Lett. 80, 5109 (1998); Physica 129D, 15 (1999)\]. The experiments also reveal the occurrence of multiple frequency states, frequency suppression of oscillations with increased time delay and the onset of both in-phase and anti-phase collective oscillations.

PACS numbers: 05.45.Xt, 87.10.+e

Coupled limit cycle oscillator models have been extensively studied in recent years because of the useful insights they provide into the collective behaviour of many physical, chemical and biological systems . One of the simplest such models, the so-called Kuramoto model , which retains only the phase information of each oscillator and is valid in the limit of weak mutual couplings, displays a spontaneous transition to a synchronized collective state of a single frequency when the coupling strength exceeds a critical value. Similar collective behaviour is observed in many natural systems, such as the synchronous flashing of fireflies, the phase locking of cardiac pacemaker cells, and the collective chirping of crickets . Phase locking has also been demonstrated experimentally in arrays of coupled nonlinear electronic circuits . When the coupling becomes stronger, amplitude effects become important and give rise to other interesting collective states, such as that of amplitude death in which the various oscillators pull each other off their periodic states and collapse to a state of zero amplitude. The condition for such a state to occur is for the oscillators to have a broad dispersion in their natural frequencies and for the coupling strength to exceed a threshold value. Thus, as has been pointed out in a number of theoretical studies , a collection of identical limit cycle oscillators cannot display amplitude death. However, more recent investigations  indicate that the presence of finite propagation time delays in the coupling removes this restriction and predict the possibility of inducing the death state even in a system of two coupled identical limit cycle oscillators. Time delay is ubiquitous in most physical systems due to finite propagation speeds of signals, finite chemical reaction times, finite response times of synapses, etc., and its influence on the collective dynamics of coupled systems can have wide-ranging implications . It is important therefore to establish the experimental feasibility of such a death phenomenon.

In this Letter we present experimental observations on time delay induced death in two coupled nonlinear circuits that are individually capable of exhibiting limit cycle oscillations. A specially designed digital delay line gives precise control over the delay time in the coupling and permits us to explore a large area of parameter space. Observations on the phenomenon of amplitude death have recently been reported  for a pair of optothermal oscillators that are thermally coupled and for which the occurrence of death for strong couplings and its relation to Hopf bifurcations of the uncoupled and coupled systems have been experimentally verified. However, as the authors themselves point out, their experimental results are not conclusive about the role of delay in the death phenomenon and in particular they have not investigated the phenomenon of death islands or, indeed, that of the multiplicity of death islands predicted by the theory of time delay induced death . In our experiments we provide clear evidence of time delay induced death islands and their multiple connectedness in the parameter space defined by the coupling strength, time delay and frequency, for the case of two coupled identical oscillators. We also find the existence of multiple frequency states, frequency suppression of oscillations with increased time delay and the onset of both in-phase and anti-phase collective oscillations.

Our experimental system, schematically shown in Fig. 1, consists of two nonlinear $`LCR`$ circuits coupled through a digital delay line (DDL). The individual oscillator circuits are a variant of the so-called Chua circuit, which has been widely employed in a number of nonlinear dynamical studies . The nonlinear resistive element $`R_N`$ has the typical $`VI`$ characteristic, with negative resistance around the origin and positive slopes away from the origin as shown in Fig. 2. This enables each individual oscillator circuit in the uncoupled state to sustain limit cycle oscillations with a characteristic frequency of $`1/\sqrt{LC}`$. The two channel DDL taps the input signal at a certain rate using an analog to digital converter (in this case an ADC0809 with a conversion rate of $`100\mu `$s) and stores it in a random access memory (RAM) bank. The stored bits are read using a logical circuit and converted back to an analog signal to be fed back into the circuit. The coupling between the oscillators is linear, resistive and proportional to the difference in the signal strengths of the two oscillators with a time delay. The coupling strength is varied by changing the resistances ($`R_K`$) that join the two oscillators. By using the current summation rules at nodes $`V_1`$ and $`V_2`$, a theoretical description of the circuit can be given as,

$$\ddot{V}_i+g(V_i)\dot{V}_i+\omega _i^2V_i=K_i[\dot{V}_j(t\tau )\dot{V}_i(t)],$$
(1)
where $`i,j=1,2,ij`$ and $`g(V_i)=(1/C_i)(0.95+0.66V_i^2+322.0V_i^4)\times 10^5`$ Amp/F, the nonlinear damping factor, is obtained from an empirical fit of the $`VI`$ characteristics of the nonlinear resistive element as shown in Fig. 2. The frequencies are $`\omega _j=1/\sqrt{L_jC_j}`$ and the coupling strengths $`K_j=1/(C_jR_{K_j})`$. For our experiments on identical oscillators we have fixed $`C_1=C_2=0.1\mu F`$ and varied $`L=L_1=L_2`$ and $`R=R_{K_1}=R_{K_2}`$ such that $`\omega =\omega _1=\omega _2`$ spans the range of $`100s^1`$ to $`1000s^1`$.

The coupling strength varies in the range of $`10s^1`$ to $`10^6s^1`$. Examples of our experimental results for the time evolution of the oscillator voltages as a function of the delay parameter are shown in Fig. 3 for $`K=1000s^1`$ and $`\omega =837s^1`$. The identical oscillators acquire an in-phase locked state (for $`\tau =0.514`$ ms), death (for $`\tau =2`$ ms) and an anti-phase locked state (for $`\tau =4.428`$ ms). A series of such observations are recorded for a fixed frequency and at different fixed values of $`K`$ by varying the time delay parameter. The output voltage $`V`$ of the oscillator is monitored on the oscilloscope up to the point where it shows a sudden drop in value by several orders of magnitude, i.e., nearly to zero. The value of $`\tau `$ is further increased up to the point where the signal from the oscillator again suddenly revives to a large value. An entire parameter space in $`K\tau `$ for amplitude death is thereby explored. In Fig. 4 the solid squares are the experimental data points marking the loci of the death regions. They clearly reveal closed death island regions in the plane of the coupling strength, $`K`$, and the time delay, $`\tau `$, for the common frequency of $`\omega =837s^1`$.

The solid and dashed lines are analytic curves obtained from a linear stability analysis of the origin for the model equation (1). These curves are derived in a standard manner  from the characteristic equation for the eigenvalue $`\lambda `$, namely, $`\lambda ^2(\lambda \frac{a_1}{C}+K)^2+2\omega ^2(\lambda \frac{a_1}{C}+K)\lambda +\omega ^4K^2\lambda ^2e^{2\lambda \tau }=0`$, and are given by,

$`\tau _1(n,K)`$ $`=`$ $`{\displaystyle \frac{(n1)\pi +\mathrm{cos}^1(1\frac{a_1}{CK})}{\frac{A}{2}+\sqrt{(\frac{A}{2})^2+\omega ^2}}},`$ (2)
$`\tau _2(n,K)`$ $`=`$ $`{\displaystyle \frac{n\pi \mathrm{cos}^1(1\frac{a_1}{CK})}{\frac{A}{2}+\sqrt{(\frac{A}{2})^2+\omega ^2}}},`$ (3)
where $`A=\sqrt{\frac{2a_1}{C}K(\frac{a_1}{C})^2}`$, $`a_1=0.95\times 10^5`$ and $`n=1,2,\mathrm{}`$. The region enclosed by the intersection of $`\tau _1(1,K)`$ and $`\tau _2(1,K)`$ is the primary death island and that enclosed between $`\tau _1(2,K)`$ and $`\tau _2(2,K)`$ is the secondary (higher order) death island. We observe that the experimental observations not only display topological similarity to the analytic predictions but are also in good quantitative agreement with them. The deviations observed are primarily due to two reasons. First, the analytic domains are sensitive to the values of the coefficients used in fitting the $`VI`$ characteristic of the nonlinear resistive element. The relative downward shift and larger width of the experimental islands compared to the analytic region are due to these estimation errors. Second, the deviation at very high values of $`K`$ (and hence low $`R_K`$) is related to limitations on the maximum current that can flow through $`R_{K_j}`$ for the rated biased voltages of our experimental circuit. The existence of the higher order island is particularly significant since it confirms the theoretical finding of multiple connectivity of the stability region of the origin .

In Fig. 5, the experimental death region is shown at a constant coupling strength $`K=500s^1`$ for different values of the frequency. The death region ceases to exist below a certain threshold frequency, in accordance with theoretical predictions . On either side of this region, there are phase locked states exhibiting in-phase locked oscillations for small time delays and anti-phase locked oscillations for longer time delays.

Another important characteristic of time delay systems in general , and delay coupled limit cycle oscillators in particular , is the existence of multiple phase locked states. When the oscillators are identical, these phase locked states coincide with in-phase (phase difference, $`\varphi =0`$) and anti-phase ($`\varphi =\pi `$) solutions. In Fig. 6(a), the domains of existence of these solutions in $`\tau `$ space are indicated. Here $`I_n`$ represent the in-phase solutions and $`A_n`$ represent the anti-phase solutions. These domains have been determined experimentally. For a fixed $`\tau `$ the various existing multiple states can be accessed with slight perturbations to the system, i.e., by slightly changing the initial conditions. We have also observed hysteretic behaviour between the various branches. In Fig. 6(b), the frequencies of the first few multiple states are plotted for increasing values of $`\tau `$. Note that the collective states can be much higher than the individual oscillator frequencies and that

their magnitude decreases with increasing values of time delay. Such a phenomenon of time delay induced frequency suppression has been discussed earlier in the context of phase only oscillators .

In conclusion, we have experimentally demonstrated for the first time the existence of delay induced amplitude death in a pair of coupled identical nonlinear oscillators . We have also demonstrated the coexistence of in-phase and anti-phase states with multiple frequencies and the suppression of the common frequency with time delay. The observation of multiply connected death islands (for identical oscillators) in the parameter space of the coupling strength and time delay is our most important result. This result not only bears topological resemblance to earlier model calculations , but also shows good quantitative agreement with analytic estimates obtained from a theoretical model of the experimental circuits. It should be noted that our experimental system of coupled nonlinear electronic circuits (and its theoretical model) is capable of richer nonlinear behaviour and goes beyond the simple normal form model used in previous theoretical studies . In this sense the results on delay induced amplitude death can be said to have a broader validity and therefore wider applicability to physically realistic systems.

We now briefly discuss a few interesting directions for future experiments that are relevant to our present results as well as to practical applications. A natural next step is to couple more than two oscillators to explore some of the other theoretical predictions discussed in . The principal technical challenge and the major expensive elements for such an endeavor are the construction of additional delay lines with enhanced number of channels. Such an effort is presently underway. Another interesting possibility is to move away from the limit cycle regime and investigate the effect of time delay on the synchronization of coupled chaotic oscillators - an area of much current scientific interest . In principle this should be possible with our present experimental system since the basic Chua circuits are capable of displaying chaotic oscillations. It would be necessary however to employ a faster delay line. Time delay could facilitate phase synchronization in view of its significant influence on phenomena like phase slips . A final area of practical interest is the stability of these limit cycle collective states to external noise. Past numerical studies indicate that time delay induced higher frequency states can exhibit metastable behaviour and decay to the lowest periodic state in the presence of sufficient external noise . We do not observe this behaviour since by design our system has very low intrinsic (thermal) noise in order to facilitate access to the higher frequency states. It would be interesting to introduce a suitably designed noise source in the individual oscillator circuits to experimentally test these theoretical predictions including the influence of noise on the threshold for onset of synchronization .

We acknowledge the important contribution of H. S. Mazumdar in the design, development and construction of the digital delay line.


text:
# Natural Mass Generation for the Sterile Neutrino

## Acknowledgments

One of the authors (T.Y) is grateful to F. Borzumati, H. Murayama and Y. Nomura for useful discussions. K.B acknowledges the Theory Group at Tokyo University for its warm hospitality. The work of T.Y is supported in part by the Grant-in-Aid, Priority Area “Supersymmetry and Unified Theory of Elementary Particles”(#707), K.B is supported in part by Department of Energy Grant No. DE-FG03-98ER41076 and by a grant from the Research Corporation.


text:
# The gluon condensate from gauge invariant vortex vacuum texture

## 1 Introduction

A revival of the vortex picture of quark confinement arose with the construction of the $`p`$-vortices which are defined after adopting the so-called center gauge  by projecting the gauge fixed link variables onto center elements . In this case, it was observed that the vortices are sensible degrees of freedom in the continuum limit : the (area) density of the $`p`$-vortices as well as their binary interactions extrapolate to the continuum . The $`p`$-vortex picture of the Yang-Mills ground state also provides an appealing explanation of the deconfinement phase transition at finite temperatures . Recently, it was pointed out that the center gauge fixing which is prior to identify the physical vortex structure might be plagued by a so-called practical Gribov problem . For this reason, as well as for rating the phenomenological importance of the vortices, a gauge invariant definition of the vortex vacuum texture is highly desired.

In this talk, I propose a method for obtaining gauge invariant vortex structures and investigate whether there is a genuine relation between these structures and the gluon condensate in SU(2) Yang-Mills theory. Extensive results and an adequate referencing can be found in a subsequent publication .

## 2 Vortices, gluons and vortex induced condensates

In order to exhibit the degrees of freedom responsible for confinement in SU(2) lattice Yang-Mills according to the $`Z(2)`$ vortex mechanism, fractionizing $`SU(2)\widehat{=}Z_2\times SO(3)`$ is instrumental. The corresponding ”coordinates” are center vortices and coset fields. Since the $`SO(3)`$ coset fields are isomorphic to algebra valued fields, these degrees of freedom are identified with the gluonic ones.

A new self–restricted cooling algorithm which reduces the $`SO(3)`$ action of the coset fields facilitates the gradual removal of the gluon fields from the lattice configurations while preserving the center degrees of freedom. For these purposes, I define a gluonic action density per link by

$$s_{x,\mu }^{gl}=\underset{\overline{\nu }\pm \mu }{}\left\{1\frac{1}{3}\text{tr}_AO_{x,\mu \overline{\nu }}\right\}=\frac{1}{3}\underset{\overline{\nu }\pm \mu }{}F_{\mu \overline{\nu }}^a[A]F_{\mu \overline{\nu }}^a[A]a^4+𝒪(a^6),$$
(1)
where $`O_{x,\mu \nu }`$ is the plaquette calculated in terms of the $`SO(3)`$ part of the $`SU(2)`$ link elements. The sum over $`\overline{\nu }`$ runs from $`4\mathrm{}4`$. $`F_{\mu \nu }^a[A]`$ is the (continuum) field strength functional of the (continuum) gluon fields $`A_\mu (x)`$ and $`a`$ is the lattice spacing. A local cooling step amounts for minimizing the action density (1) with respect to the coset fields at a given link. A self-restriction is imposed by rejecting the cooling of the adjoint link iff the gluonic action is smaller than some threshold value $`s_{x,\mu }^{gl}<\mathrm{\hspace{0.33em}8}\kappa ^4a^4.`$ Thereby $`\kappa `$ is a gauge invariant cooling scale of mass dimension one. For $`\kappa =0`$, the cooling procedure completely removes the gluon fields from the SU(2) lattice configurations leaving only center degrees of freedom. For $`\kappa `$ smaller than the string tension, a clear signal of the vortex vacuum texture is noticed by a large SU(2) action density which is accumulated at the two dimensional vortex world sheets of the 4-dimensional space time. Since the new cooling procedure is gauge covariant and since the vortex structure manifests itself in the gauge invariant $`SU(2)`$ action density, the vortex texture is gauge invariant .

In order to study the relevance of this vortex texture, I calculated the force between a static quark anti-quark pair from $`SO(3)`$ cooled configurations.

I find that the cooling procedure has strongly affected the force at short distances. This is expected since the behavior at small $`r`$ is dominated by the exchange of gluons, which are already partially eliminated by cooling. Most important is, however, that the string tension is unchanged by the cooling procedure.

Hadronic correlators at high momentum transfer are accessible by experiments. On the other, these correlators are sufficiently under the control of the operator product expansion (OPE), where non–perturbative properties of the Yang–Mills (or QCD) vacuum are parameterized by so–called condensates . The mass dimension four operator, $`GC`$, is provided by the trace of the energy momentum tensor, and, hence proportional to the $`SU(2)`$ action density in the present case. Here, I will use the new cooling method to uncover the contribution of the vortices to the mass dimension four condensate of $`SU(2)`$ Yang–Mills theory. It turns out, that the vortex contribution to $`GC`$ shows appropriate scaling towards the continuum limit . The mass dimension four condensate $`GC`$ in physical units is also shown in figure 1 as function of $`\kappa ^4`$. At large values of $`\kappa `$, the contribution of the coset, say gluon, fields to the condensate is large and $`GC`$ linearly rises with $`\kappa ^4`$. However, the most striking feature is that $`GC`$ approaches a finite value in the limit $`\kappa 0`$, thus showing a non-trivial dimension four condensate which is of pure vortex origin. The limiting value of $`0.19\text{GeV}^4`$, is in rough agreement with recent values for the gluon condensate quoted in the literature .

## Acknowledgments

It is a pleasure to thank my collaborators E.-M. Ilgenfritz and H. Reinhardt.

## References


text:
# 1 Introduction

## 1 Introduction

In the standard afterglow shock model (for a review see ), a gamma-ray burst (GRB) afterglow is usually believed to be produced by synchrotron radiation or inverse Compton scattering in an ultrarelativistic shock wave expanding in a homogeneous medium. As more and more ambient matter is swept up, the shock gradually decelerates while the emission from such a shock fades down, dominating at the beginning in X-rays and progressively at optical to radio energy band. The standard model is based on four basic assumptions: (1) the total energy of the shock is released impulsively before its formation; (2) the medium swept up by the shock is homogeneous and its density ($`n`$) is the one of the interstellar medium $`1\mathrm{cm}^3`$; (3) the electron and magnetic field energy fractions of the shocked medium and the index ($`p`$) in the accelerated electrons’ power-law distribution are constant during the whole evolution stage; and (4) the shock is spherical.

Each of these assumptions has been varied to discuss why some observed afterglows deviate from that expected by the standard afterglow model. For example, the R-band light curve of GRB 970508 afterglow peaks around two days after the burst, and there is a rather rapid rise before the peak which is followed by a long power-law decay. There are two models explaining this special feature: (i) It was envisioned  that a postburst fireball may contain shells with a continuous distribution of Lorentz factors. As the external forward shock sweeps up ambient matter and decelerates, internal shells will catch up with the shock and supply energy into it. A detailed calculation shows that this model can explain well this special feature . (ii) We considered continuous energy injection from a strongly magnetized millisecond pulsar into the shock through magnetic dipole radiation . This model can also account for well the observations. It is very clear that these models don’t use basic assumption (1).

There are several models in the literature that discuss the effect of inhomogeneous media on afterglows , dropping the second assumption. Generally, an $`nr^k`$ ($`k>0`$) medium is expected to steepen an afterglow’s temporal decay. We studied, for the first time, the wind model for afterglows based on the Blandford-McKee self-similar solution  of a relativistic adiabatic shock, and suggested that GRB 970616 is an interactor with a stellar wind of $`nr^2`$ . It was found  that a Wolf-Rayet star wind likely leads to an $`nr^2`$ medium, and thus if GRB 980519 resulted from the explosion of such a massive star, subsequent evolution of a relativistic shock in this medium is consistent with the steep decay in the R-band light curve of the afterglow from this burst. Another way of dropping the second assumption is that the density of an ambient medium is invoked to be as high as $`n10^6\mathrm{cm}^3`$. The temporal decay of the R-band afterglow of GRB 990123 has been detected to steepen about 2.5 days after this burst . We proposed a plausible model in which a shock expanding in a dense medium has evolved from a relativistic phase to a nonrelativistic phase . We found that this model fits well the observational data if the medium density is about $`3\times 10^6\mathrm{cm}^3`$. We further suggested that such a medium could be a supernova or supranova or hypernova ejecta.

In basic assumption (3), the electron and magnetic field energy fractions of the shocked medium may not be varied during whole evolution, as argued in , where all the observational data including both the prompt optical flash and the afterglow of GRB 990123 were analyzed.

The steepening in the light curves of the afterglows of some bursts may also be due to lateral spreading of a jet, as analyzed in  when the jet expands in a homogeneous interstellar medium (ISM). This in fact drops basic assumption (4). However, numerical studies of  show that the break of the light curve is weaker and smoother than the one analytically predicted when the light travel effects related to the lateral size of the jet and a realistic expression of the lateral expansion speed are taken into account. In the case of a jet expanding in a wind, the calculated light curve is even much weaker and smoother than the ananlytical one . We recently calculated light curves for GRB afterglows when anisotropic jets ($`dE/d\mathrm{\Omega }\theta ^k`$) expand both in the interstellar medium and in the wind medium . We found that in each type of medium, one break appears in the late-time afterglow light curve for small $`k`$ but becomes weaker and smoother as $`k`$ increases. When $`k2`$, the break seems to disappear but the afterglow decays rapidly. Thus, we expect that the emission from expanding, highly anisotropic jets provides a plausible explanation for some rapidly fading afteglows whose light curves seem to have no break.

We discussed the dense medium model in more details , by taking into account both the synchrotron self-absorption effect in the shocked medium and the energy injection effect of . Recently, we combined the dense medium model with the pulsar energy injection effect to provide a good fit to the optical afterglow data of GRB 000301C . Here we want to give a brief review of some of our studies on GRB afterglows.

## 2 Shock Evolution

It is well known that the evolution of a partially radiative shock depends on both the efficiency with which the shock transfers its bulk kinetic energy to electrons and magnetic fields and on the efficiency with which the electrons radiate their energy. In 1998, we proposed, for the first time, a unified model for dynamical evolution of a partially radiative shock . This model is not only valid during the whole evolution stage including the Sedov phase for an adiabatic shock, but also can describe well an adiabatic shock as well as a highly radiatve shock. This model was later re-investigated and referred to as a generic one in . For simplicity, we here assume that a relativistic shock expanding in a dense medium is adiabatic. This assumption is correct particularly for a low electron energy density fraction in the shocked medium . The Blandford-McKee self-similar solution  gives the Lorentz factor of an adiabatic relativistic shock,

$$\gamma =\frac{1}{4}\left[\frac{17E_0(1+z)^3}{\pi nm_pc^5t_{}^3}\right]^{1/8}=1E_{52}^{1/8}n_5^{1/8}t_{}^{3/8}[(1+z)/2]^{3/8},$$
(1)
where $`E_0=E_{52}\times 10^{52}\mathrm{ergs}`$ is the total isotropic energy, $`n_5=n/10^5\mathrm{cm}^3`$, $`t_{}`$ is the observer’s time since the gamma-ray trigger in units of 1 day, $`z`$ is the the redshift of the source generating this shock, and $`m_p`$ is the proton mass. We assume $`\gamma =1`$ when $`t_{}=t_b`$. This implies

$$n_5=E_{52}t_b^3[(1+z)/2]^3.$$
(2)
For $`t_{}>t_b`$, the shock will be in a nonrelativistic phase. In the following we will discuss the spectrum and light curve during the non-relativistic phase.

As usual, only synchrotron radiation from the shock is considered. To analyze the spectrum and light curve, one needs to know three crucial frequencies: the synchrotron peak frequency ($`\nu _m`$), the cooling frequency ($`\nu _c`$), and the self-absorption frequency ($`\nu _a`$). We assume a power law distribution of the electrons accelerated by the shock: $`dn_e^{}/d\gamma _e\gamma _e^p`$ for $`\gamma _e\gamma _{em}`$, where $`\gamma _e`$ is the electron Lorentz factor and $`\gamma _{em}=610ϵ_e(\gamma 1)`$ is the minimum Lorentz factor. We further assume that $`ϵ_e`$ and $`ϵ_B`$ are the electron and magnetic energy density fractions of the shocked medium respectively. The $`\nu _m`$ is the characteristic synchrotron frequency of an electron with Lorentz factor of $`\gamma _{em}`$, while the $`\nu _c`$ is the characteristic synchrotron frequency of an electron which cools on the dynamical age of the shock. According to Sari et al. , we have derived the synchrotron peak frequency, the cooling frequency and the synchotron self-absorption frequency, measured in the observer’s frame . They are correct for the whole evolution stage.

Now we give the spectrum and light curve of the afterglow during the non-relativistic phase. First, for the case without energy injection, the shock velocity decays as $`t_{}^{3/5}`$ and thus we have

$$F_\nu =\{\begin{array}{ccc}(\nu _a/\nu _m)^{(p1)/2}(\nu /\nu _a)^{5/2}F_{\nu _m}\nu ^{5/2}t_{}^{11/10}\hfill & \mathrm{if}\nu <\nu _a\hfill & \\ (\nu /\nu _m)^{(p1)/2}F_{\nu _m}\nu ^{(p1)/2}t_{}^{(2115p)/10}\hfill & \mathrm{if}\nu _a<\nu <\nu _c\hfill & \\ (\nu _c/\nu _m)^{(p1)/2}(\nu /\nu _c)^{p/2}F_{\nu _m}\nu ^{p/2}t_{}^{(43p)/2}\hfill & \mathrm{if}\nu >\nu _c.\hfill & \end{array}$$
(3)
We easily see that for high-frequency radiation the temporal decay index $`\alpha =(2115p)/10`$ for emission from slow-cooling electrons or $`\alpha =(43p)/2`$ for emission from fast-cooling electrons. If $`p2.8`$, then $`\alpha 2.1`$ or $`2.2`$. Comparing this with the relativistic result, we conclude that the afterglow decay steepens at the nonrelativistic stage.

Some models for GRB energy sources (for a brief review see ) predict that during the formation of an ultrarelativistic fireball required by GRB, a strongly magnetized millisecond pulsar will be born. If so, the pulsar will continuously input its rotational energy into the forward shock of the postburst fireball through magnetic dipole radiation because electromagnetic waves radiated by the pulsar will be absorbed in the shocked medium . Since an initially ultrarelativistic shock discussed in  rapidly becomes nonrelativistic in a dense medium, we next investigate the evolution of a nonrelativistic adiabatic shock with energy injection from a pulsar. The total energy of the shock is the sum of the initial energy and the energy which the shock has obtained from the pulsar:

$$E_0+_0^t_{}L𝑑t_{}=E_{\mathrm{tot}}v^2r^3,$$
(4)
where $`L`$ is the stellar spindown power $`(1+t_{}/T)^2`$ ($`T`$ is the initial spindown time scale). The term on the right-hand side is consistent with the Sedov solution. Please note that $`L`$ can be thought of as a constant for $`t_{}<T`$, while $`L`$ decays as $`t_{}^2`$ for $`t_{}T`$. Because of this feature, we easily integrate the second term on the left-hand side of equation (17). We now define a time at which the shock has obtained energy $`E_0`$ from the pulsar, $`t_c=E_0/L`$, and assume $`t_cT`$. The evolution of the afterglow from such a shock can be divided into three stages.

Stage (i): $`t_{}t_c`$, viz., the second term on the left-hand side of equation (4) can be neglected. The evolution of the afterglow is the same as in the above case without any energy injection.

Stage (ii): for $`T>t_{}t_c`$, the term $`E_0`$ in equation (4) can be neglected. At this stage, the shock’s velocity $`vt_{}^{2/5}`$, we have derived the spectrum and light curve of the afterglow

$$F_\nu =\{\begin{array}{ccc}(\nu _a/\nu _m)^{(p1)/2}(\nu /\nu _a)^{5/2}F_{\nu _m}\nu ^{5/2}t_{}^{7/5}\hfill & \mathrm{if}\nu <\nu _a\hfill & \\ (\nu /\nu _m)^{(p1)/2}F_{\nu _m}\nu ^{(p1)/2}t_{}^{(125p)/5}\hfill & \mathrm{if}\nu _a<\nu <\nu _c\hfill & \\ (\nu _c/\nu _m)^{(p1)/2}(\nu /\nu _c)^{p/2}F_{\nu _m}\nu ^{p/2}t_{}^{2p}\hfill & \mathrm{if}\nu >\nu _c.\hfill & \end{array}$$
(5)
It can be seen that for high-frequency radiation the temporal decay index $`\alpha =(125p)/50.4`$ for emission from slow-cooling electrons or $`\alpha =2p0.8`$ for emission from fast-cooling electrons if $`p2.8`$. This shows that the afterglow decay may significantly flatten due to the effect of the pulsar.

Stage (iii): for $`t_{}T`$, the power of the pulsar due to magnetic dipole radiation rapidly decreases as $`Lt_{}^2`$, and the evolution of the shock is hardly affected by the stellar radiation. Thus, the evolution of the afterglow at this stage will be the same as in the above case without any energy injection.

In summary, as an adiabatic shock expands in a dense medium from an ultrarelativistic phase to a nonrelativistic phase, the decay of radiation from such a shock will steepen, subsequently may flatten if a strongly magnetic millisecond pulsar continuously inputs its rotational energy into the shock through magnetic dipole radiation, and finally the decay will steepen again due to disappearance of the stellar effect. In the next section, we will see how to explain some unusual afterglows based on the above conclusion.

## 3 Some Unusual Afterglows

### 3.1 GRB 980519

The optical afterglow $`8.5`$ hours after GRB 980519 decayed as $`t_{}^{2.05\pm 0.04}`$ in BVRI , while the power-law decay index of the X-ray afterglow $`\alpha _X=2.07\pm 0.11`$ , in agreement with the optical. The spectrum in optical band alone is well fitted by a power low $`\nu ^{1.20\pm 0.25}`$, while the optical and X-ray spectra together can also be fitted by a single power law of the form $`\nu ^{1.05\pm 0.10}`$. In addition, the radio afterglow of this burst was observed by the VLA at 8.3 GHz, and its temporal evolution $`t_{}^{0.9\pm 0.3}`$ between 1998 May 19.8UT and 22.3UT .

We now analyze the observed afterglow data of GRB 980519 based on our model. We assume that for this burst, the forward shock evolved from an ultrarelativistic phase to a nonrelativistic phase in a dense medium at $`8`$ hr after the burst. So, the detected afterglow, in fact, was the radiation from a nonrelativistic shock. This implies $`\gamma 1`$ at $`t_b1/3`$ days. From equation (2), therefore, we find

$$n_527E_{52}[(1+z)/2]^3.$$
(6)
If $`p2.8`$, and if the observed optical afterglow was emitted by slow-cooling electrons and the X-ray afterglow from fast-cooling electrons, then according to equation (3), the decay index $`\alpha _R=(2115p)/102.1`$ and $`\alpha _X=(43p)/22.2`$, in excellent agreement with observations. Furthermore, the model spectral index at the optical to X-ray band and the decay index at the radio band, $`\beta =(p1)/20.9`$ and $`\alpha =1.1`$, are quite consistent with the observed ones, $`1.05\pm 0.10`$ and $`0.9\pm 0.3`$, respectively.

We  took into account three observed data which correpond to the radio, R-band and X-ray frequencies respectively, and inferred intrinsic parameters of the shock and the redshift of the burst, $`ϵ_e0.16`$, $`ϵ_B2.8\times 10^4`$, $`E_{52}0.27`$, $`n_53.4`$, and $`z0.55`$. After considering these reasonable parameters, we numerically studied the trans-relativistic evolution of the shock  and found that our dense medium model can provide an excellent fit to all the observational data of the radio afterglow from GRB 980519 shown in .

### 3.2 GRB 000301C

The optical afterglow data of GRB 000301c were presented in . In addition, the spectral index $`\beta =1.1\pm 0.1`$. We  combined the dense medium model with the pulsar energy injection effect to explain the unusual optical afterglow. For stage (ii), if $`p=3.4`$, then $`\alpha =(125p)/5=1.0`$ and $`\beta =(p1)/2=1.2`$ are consistent with the GRB 000301c R-band afterglow data in initial 7.5 days after the burst. These data indicate $`\alpha _11.1`$, which implies $`\alpha _{\mathrm{obs}}\beta _{\mathrm{obs}}`$ at early times. If the afterglow were radiated by fast-cooling electrons in the shocked medium, we would find $`\alpha =2(1\beta )`$, which is clearly inconsistent with the observational result. Therefore, the GRB 000301c R-band afterglow arose from those slow-cooling electrons in the shocked medium. For stage (iii), in the case of $`p=3.4`$, the model’s time index $`\alpha =(2115p)/10=3.0`$ is quite consistent with the observational data of the GRB 000301c R-band afterglow at late times, $`\alpha _2=3.01\pm 0.53`$ . We also carried out simulations of the evolution of a shock with energy injection from a pulsar and the resulting emission . Our numerical results indeed show one sharp break in the late-time afterglow light curve and give a good fit to the R-band afterglow data of GRB 000301c.

## 4 Conclusions

We discussed the evolution of an adiabatic shock expanding in a dense medium from an ultrarelativistic phase to a nonrelativistic phase in more details in this paper. In particular, we discussed the effects of synchrotron self absorption and energy injection on the afterglow from this shock. In a dense medium, the shock becomes nonrelativistic rapidly after a short relativistic phase. This transition time varies from several hours to a few days when the medium density is from $`10^5`$ to a few $`\times 10^6\mathrm{cm}^3`$, and the shock energy from $`10^{51}`$ to $`10^{54}`$ ergs. The afterglow from the shock at the nonrelativistic stage decays more rapidly than at the relativistic stage, while the decay index varies from $`1.35`$ to $`2.1`$ if the spectral index of the accelerated electron distribution, $`p=2.8`$, and the radiation comes from those slow-cooling electrons. Since some models mentioned above predict that a strongly magnetic millisecond pulsar may be born during the formation of GRB, we also discuss the effect of such a pulsar on the evolution of a nonrelativistic shock through magnetic dipole radiation, in contrast to the case discussed in . We found that after the energy which the shock obtains from the pulsar is much more than the initial energy of the shock, the afterglow decay will flatten significantly and the decay index will become $`0.4`$. When the pulsar energy input effect disappears, the index will still be $`2.1`$. These features are in excellent agreement with the afterglow of GRB 980519. Furthermore, our model fits very well all the observational data of GRB 980519 including all the radio data. Our model also provides a good fit to the R-band afterglow data of GRB 000301c.

### Acknowledgments

We would like to thank T. Lu, D. M. Wei and Y. F. Huang for helpful discussions. This work was supported by the National Natural Science Foundation of China (grant 19825109) and by the National 973 Project.


text:
# 1 The time-dependent 𝛾₋₁ index for three classes of fluids we are interested in: vacuum, radiation and matter. The plot corresponds to an open universe of 𝜆=0.2 along the comoving time axis. There is a small region associated to each 𝛾₋₁ where it is negative (Chaplygin-like fluid).

DARBOUX COSMOLOGICAL FLUIDS IN COMOVING TIME

RENATO KLIPPERT1,3, HARET C. ROSU2,3, REMO RUFFINI3

1 Brazilian Center for Research in Physics, R. Dr. Xavier Sigaud 150 Urca, 22290-180 Rio de Janeiro RJ, Brazil

2 Instituto de Física, Universidad de Guanajuato, Apdo Postal E-143, León, Gto, Mexico

3 University of Rome “La Sapienza” and ICRA, I-00185 Rome, Italy

The barotropic indices and the corresponding FRW scale factors of the so-called Darboux cosmological fluids are presented in the comoving time axis, which is the natural one for the phenomenology related to the cosmological data. Some useful comments on the features of the plots are included.

In a previous work, an interesting class of Darboux cosmological fluids in closed and open FRW models of the $`\mathrm{\Lambda }=0`$ cosmology has been introduced by one of the authors.1 However, the results have been displayed in the conformal time axis, whereas for comparison with definite cosmological data the comoving time axis is usually needed. The purpose of this letter is to present the cosmological comoving evolution of that class of fluids and add several heuristic comments on the obtained results. We plot here the comoving time dependent barotropic indices and the corresponding scale factors for open FRW universes,2 respectively. Our general conclusions are presented in the following.

(i) For the open cases, we find that $`\gamma _1(t,\lambda )`$ displays a single deep structure where an accelerating universe can exist, $`\gamma _1(t,\lambda )<2/3`$. Moreover, depending on the value of the $`\lambda `$ parameter,1 a negative region of limited extent may occur indicating the presence of a Chaplygin-like fluid whose origin may be due to the contribution of quantum field extended objects, such as d-branes.3

(ii) For the closed, matter- and radiation-dominated cases, we find a more complicated damped periodic $`\gamma _{+1}(t,\lambda )`$, for which the accelerating region is present only in the first period. No Chaplygin-like region occurs in these cases, which points to more stable cosmological conditions. The vacuum-dominated closed case has a Chaplygin-like (unstable) oscillatory behavior forever.

A very definite prediction of the model is that there is a single localized accelerating region (SLAR) for all types of non flat FRW universes, except for the closed vacuum case. Concerning a possible application to cosmological data,4 this means that only the sample of supernovae with redshifts corresponding to the SLAR can show a clear-cut accelerating effect. Such behavior can be tested by means of a precise calibration of the supernovae cosmological data, which should take into account the determination of the time scale $`T`$ introduced in the plots.

Acknowledgements

RK would like to acknowledge Brazilian CAPES foundation for a grant.

References


text:
# Dynamics of Black-Hole Nuclei

## 1. Introduction

Supermassive black holes (BHs) are generic components of galaxies, and they appear to be linked in fundamental ways to the dynamics of the stellar component, both on large and small scales. An astonishingly strong correlation exists between BH mass and stellar velocity dispersion, $`M_{}\sigma ^\alpha `$, $`\alpha 5`$ (Ferrarese & Merritt 2000; Gebhardt et al. 2000b). The $`M_{}\sigma `$ relation as defined by the best-determined BH masses is so tight that it surpasses in predictive accuracy what can be achieved from detailed dynamical modelling of high-quality stellar-kinematical data in most galaxies. On small scales, the approximately power-law variation of stellar density with radius (Ferrarese et al. 1994) would be difficult to understand in the absence of BHs. Low luminosity ellipticals and bulges have steep nuclear density profiles, $`\rho r^\gamma `$, $`\gamma 2`$; such steep cusps form naturally in stellar systems where the BHs grow on time scales long compared to crossing times (Peebles 1972; Quinlan, Hernquist & Sigurdsson 1995). Brighter galaxies typically have weaker cusps, $`\rho r^1`$ (Merritt & Fridman 1995; Gebhardt et al. 1996). While no universally accepted model has yet been proposed for the origin of the weak cusps, we note here one feature that suggests a link to BHs. Two galaxies, NGC 3379 and M87, have weak cusps with well-determined structural parameters and also have BHs with accurately determined masses. Table 1 gives for each galaxy the “break” radius $`r_b`$ at which the central power law turns over to a steeper outer profile; and also the stellar mass $`M_{}`$ contained within $`r_b`$. In both galaxies, $`M_{}(r_b)`$ is identical within the uncertainties with $`M_{}`$, the mass of the BH; this is in spite of a factor of $`30`$ difference in $`M_{}`$ and $`6`$ in $`r_b`$.

The rough equality between $`M_{}`$ and $`M_{}(r_b)`$, and the exclusive association of weak cusps with bright galaxies, is consistent with the suggestion that weak cusps are relics of galaxy mergers (Begelman et al. 1980; Ebisuzaki et al. 1991). In this model, the two BHs from the merging galaxies fall to the center and heat the stars, reducing the stellar density and forming a tight binary. The “damage” done by the binary, and hence the extent of the shallow cusp, would be roughly proportional to the total BH mass.

New instruments and observational techniques, both ground-based (e.g. SAURON; Miller 2000) and space-based (HST; e.g. Joseph et al. 2000), are giving us increasingly detailed glimpses of the stellar kinematics very near the centers of galaxies. These data are useful both for constraining BH masses, but also as fossil signatures of the stellar dynamical processes that shaped nuclei. Interpreting these new data in the context of models for galaxy formation will be a challenging computational task. Stellar densities in a steep stellar cusp are $`>10^7M_{}`$ pc<sup>-3</sup> implying dynamical times of $`<10^4`$ yr, and even less in the vicinity of the BH; thus an $`N`$-body code that simulates both the large- and small-scale evolution of merging galaxies must be able to handle a range of $`10^5`$ in time scales as well as a sufficiently large number of particles to resolve the stellar cusps. Until recently, $`N`$-body codes satisfying these requirements were almost non-existent, but happily that situation has now changed. Tree codes that efficiently implement variable time steps on parallel architectures now exist (e.g. GADGET, Springel et al. 2000); such codes are ideal for following the early stages of galaxy mergers before the BHs form a bound pair. In the later stages, when the BH binary begins to interact “collisionally” with individual stars, a code like NBODY6++ (Spurzem & Baumgardt 1999) is more appropriate: this code implements Aarseth’s efficient hierarchical algorithm using block time steps on parallel computers and can treat particle-particle interactions on arbitrarily small scales.

Merger simulations will need to reproduce a number of observed regularities in the properties of the stellar cusps. One is the tight correlation between $`M_{}`$ and $`\sigma `$ mentioned above. Another observational result is the dependence of cusp slope $`\gamma `$ on galaxy luminosity (Gebhardt et al. 1996): faint ellipticals, $`|M_v|>20`$, have $`\gamma 2`$ while $`\gamma `$ for bright ellipticals varies from $`2`$ to $`0`$. There are at least three interesting questions here: (1) What determines the characteristic luminosity separating galaxies with strong and weak cusps? (2) Assuming that all galaxies begin life with strong cusps, how do the weak cusps form? (3) How are the weak cusps in bright galaxies maintained in spite of mergers with dense, low-luminosity galaxies? Another empirical result is the rough proportionality between break radius $`r_b`$ and galaxy luminosity (Faber et al. 1997). A direct link between cusp slope and BH mass has also been suggested (van der Marel 1999).

In addition to these well-established parameter relations, the new, high-resolution kinematical data will provide detailed stellar velocity fields in galactic nuclei which can be compared to $`N`$-body simulations.

In this article, we summarize the existing $`N`$-body work on interactions between galaxies with high central densities and BHs. We then present preliminary results from new simulations which, for the first time, follow the evolution of BH nuclei through a realistic merger to the formation and decay of a hard BH binary.

## 2. $`N`$-Body Studies of Mergers with Stellar Cusps and Black Holes

In the merger of two galaxies each containing a supermassive nuclear BH, the dynamical evolution can usefully be divided into two stages, before and after the formation of a hard BH binary. Quinlan (1996) defines a “hard” BH binary as one in which the semimajor axis $`a`$ satisfies

$$a<a_h=\frac{Gm_2}{4\sigma ^2}=2.7\mathrm{pc}\left(\frac{m_2}{10^8M_{}}\right)\left(\frac{200\mathrm{k}\mathrm{m}\mathrm{s}^1}{\sigma }\right)^2$$
(1)
with $`m_2`$ the mass of the smaller of the two BHs and $`\sigma `$ the stellar velocity dispersion. Before the binary becomes hard, each BH acts as an independent point mass and its orbit decays in roughly the way described by Chandrasekhar’s dynamical friction formula (except that the relevant mass is that of the BH plus any stars from the original cusp that remain bound to it). The evolution in this regime is essentially collisionless, i.e., independent of the masses of the background stars. After the formation of the hard binary, the two BHs move about the galactic center like a single particle of mass $`M_{12}=m_1+m_2`$ and their interaction with passing stars is dominated by three-body encounters. The latter almost always result in the transfer of energy from the BH binary to the stars; as a consequence, stars are ejected from the nucleus and the binary shrinks (Hills 1992; Quinlan 1996).

Simulating the merger of two galaxies, each containing a dense stellar cusp and a central point mass, is computationally demanding, and essentially no such simulations have appeared in the literature to date. Barnes (1999) described mergers between identical spherical galaxies with Dehnen’s (1993) density law, having central power-law cusps with indices $`\gamma =(1,1.5,2)`$, but no BHs. He used a tree code with fixed time step. Barnes found that the $`\gamma =1`$ and $`\gamma =1.5`$ cusps became slightly shallower following the merger but that the $`\gamma =2`$ cusps were nearly unchanged. Holley-Bockelmann & Richstone (1999) obtained similar results in simulations of unequal-mass mergers, also without BHs; they fixed the potential corresponding to the larger galaxy and assumed that the smaller galaxy remained spherical as its orbit decayed. They found that the secondary remained intact as long as its initial density was higher than that of the primary, as expected in mergers between real galaxies. Neither of these studies was able to produce weak cusps in galaxies that did not contain them initially.

Merger simulations including nuclear BHs have been carried out but only from rather unrealistic initial conditions, typically consisting of galaxies with large, constant-density cores (Ebisuzaki et al. 1991; Makino et al. 1993; Governato et al. 1994; Makino & Ebisuzaki 1996). The BHs were found to generate still larger cores in the merger products due to heating of the stars; the core tended to grow in such a way as to maintain rough proportionality between core radius and half-mass radius, and between core mass and (total) BH mass. This result is consistent with observed parameter relations (Faber et al. 1997) if core radius is identified with break radius $`r_b`$. The central density profiles in some of these simulations showed weak, power-law cusps. Nakano & Makino (1999a,b) argued that this was due simply to the removal of low-energy stars by heating from one or both of the BHs. Holley-Bockelmann & Richstone (2000) extended their earlier, BH-free simulations to the case of a fixed primary containing a central point mass and an evolving secondary with no BH. They found that the secondary was usually disrupted by tidal forces from the primary’s BH; however this result might be modified if a BH were added to the secondary as well.

A collisional $`N`$-body code is required to follow the evolution once the BH binary becomes hard. This statement is at first sight surprising, since the rate of evolution of a hard binary in a background of low-mass perturbers depends only on their total density, not their individual masses (Hills 1992; Quinlan 1996). However the rate of binary decay depends also on how widely the binary wanders about the center of the potential. A fixed binary rapidly ejects all of the stars that pass near it; continued evolution requires that the binary wander to regions of higher density. It is able to do this because the low stellar density implies a small gravitational restoring force, and because the binary receives kicks from ejected stars. The resultant random walk depends on the size of the kicks, hence on the mass ratio $`m_{}/M_{}`$ and the number of particles $`N`$. The decay rates observed in the simulations with largest $`N`$, $`N3\times 10^5`$ (Makino 1997; Quinlan & Hernquist 1997; Hemsendorf 1999), imply a binary lifetime (until gravitational radiation coalescence) of order a Hubble time (Merritt 2000); if one extrapolates the decay rate to still larger $`N`$ using the $`N`$-dependence estimated by Makino (1997), $`N^{1/2}`$, one concludes that BH binaries in real galaxies might be expected to often survive long enough for a third BH to find its way to the center (Valtonen 1996). However Quinlan & Hernquist (1997) saw evidence that the decay rate may saturate for $`N>10^5`$. Another possibility is that gas dynamical processes accelerate the decay (Gould & Rix 2000).

Although the timescale for evolution of massive BH binaries to complete coalescence is still unclear, most of the damage that the binary inflicts on the stellar cusp takes place shortly after it becomes hard. The mass ejected is

$$M_{ej}M_{12}\mathrm{ln}(a_h/a)$$
(2)
(Quinlan 1996), hence the binary ejects of order its own mass after shrinking by a factor of only two. This is probably sufficient to convert a strong cusp into a weak cusp (Table 1). Quinlan & Hernquist (1997) presented density profiles that look qualitatively similar to those in bright galaxies, with a clear break radius; some of these profiles even exhibit central minima. The decaying binary also has a strong influence on the nuclear kinematics, since ejected stars are predominantly on eccentric orbits; the remaining stars exhibit a tangentially-biassed velocity ellipsoid. Detecting this velocity polarization in real galaxies would be difficult due to the low surface brightnesses of weak-cusp galaxies.

Taken together, these simulations probably reproduce much of the interesting dynamics that takes place during the merger of two galaxies containing nuclear BHs. Among the interesting questions yet to be answered are: (1) Can a shallow power-law cusp be produced from the merger of two galaxies with steep cusps and nuclear BHs? (2) Can a shallow cusp survive the accretion of a small dense galaxy with a nuclear BH? (3) What is the precise $`N`$-dependence of the wandering amplitude and decay rate of a BH binary in a galactic nucleus? (4) Is the structure of the resultant stellar core also $`N`$-dependent?

In addressing these questions, one would ideally like to simulate mergers in a continuous way from the largest, kiloparsec scales down to the sub-parsec scales that characterize the BH binary. We describe in the next section a first attempt to carry out this ambitious program.

## 3. A Merger Simulation

We report the simulation of a merger between two spherical, equal-mass galaxies with steep central density cusps surrounding nuclear BHs. We first integrated the entire merger using the tree code GADGET (Springel et al. 2000), a parallel algorithm with continuously variable time steps. GADGET has a fixed spatial resolution determined by a softening length $`h`$, which was chosen to be slightly smaller than $`a_h`$, the expected separation of the hard BH binary. To follow the evolution of the binary at late times when $`a<h`$, we continued the integration using NBODY6++, a parallel, multiple-time-step version of Aarseth’s direct-summation code. The particle coordinates as computed with GADGET were extracted at a time before the formation of the hard binary and used to generate initial conditions for NBODY6++. The early evolution as computed by NBODY6++ was compared to that of GADGET over the same time interval to ensure that the tree code had accurately integrated the equations of motion prior to extraction of the coordinates.

The initial galaxies were generated from Dehnen’s (1993) law with central density dependence $`\rho r^\gamma `$, $`\gamma =2`$. Initial velocities were assigned based on an isotropic distribution function that accounts for the presence of the BH (Tremaine et al. 1994). Each galaxy was given $`N=1.31\times 10^5`$ equal-mass particles; the two particles representing the BHs were given masses $`1\%`$ that of their host galaxies. Two, randomly-selected subsets of particles were extracted for the later integration with NBODY6++, with $`N=32,768`$ and $`N=65,536`$. The initial separation between the galaxies was four times the half-mass radius of each, where $`r_{1/2}=1.0`$, and the initial orbital velocity was $`1/2`$ that of a circular orbit. The internal orbital period at the half-mass radius of either galaxy was 8.8 in $`N`$-body units; the GADGET integration was carried out for $`2`$ of these periods and required approximately 10 days using 8 processors on the Rutgers Sun E10000 supercomputer. At the time of this writing, the NBODY6++ integration with $`N=32768`$ had advanced until a time of $`20`$ and the integration with $`N=65,536`$ until a time of $`12`$.

The separation between the two BH particles as a function of time is shown in Fig. 1a. As computed by GADGET, the binary separation stalls at $`ah`$, the softening length, soon after the BHs have fallen to the center and formed a hard binary. With NBODY6++, the binary continues to slowly decay as it exchanges energy with the stars. However most of the damage done by the binary to the cusp takes place very shortly after the binary forms. The central density (Fig. 1b) drops by a large factor in just a few orbital periods of the binary. This behavior is consistent with a simple model (Merritt 2000) which predicts a mass ejected by the binary of

$$M_{ej}=M_{12}\mathrm{ln}\left[\frac{4(tt_h)}{t_0}\right]$$
(3)
with $`t_h`$ the time of formation of the hard binary and $`t_01.76GM_{12}/\sigma _{}^30.02`$ in $`N`$-body units. The transformation of two density cusps into one is shown in Fig. 2; here again, the very sudden drop in density at $`tt_h`$ is apparent.

Fig. 3 shows the change in the stellar density profile during the time of its most rapid evolution, $`t>t_h`$. The cusp evolves quickly from its initial $`r^2`$ dependence to a shallower power law, $`\rho r^1`$, and a break appears at a radius $`0.032.4\times a_h`$. Thereafter the profile continues to evolve toward shallower slopes but very slowly, as the binary continues to eject stars. This figure confirms for the first time that the merger of two galaxies with steep cusps and BHs can produce a weak, power-law cusp. Bright ellipticals and bulges sometimes have even shallower central profiles (Gebhardt et al. 1996). Increasing the masses of the BHs would probably produce weaker cusps (e.g. Quinlan & Hernquist 1997), although the BHs in our simulation are already uncomfortably large. The continued, slow evolution in the central density as the binary ejects stars (Fig. 1b) might also produce a shallower profile; repeated mergers would probably have have a similar effect. Future work should be directed toward understanding which of these explanations accounts for the observed trend of cusp slope with galaxy luminosity.

The sharply lowered density at $`t>t_h`$ causes another change in the behavior of the binary, as shown in Fig. 4. The binary begins suddenly to wander over a region much larger than $`a_h`$. Wandering was inferred by Quinlan & Hernquist (1997) and Makino (1997) in their $`N`$-body studies but little is known about its properties. The wandering is enabled by the reduced force gradients in the low-density core; it is driven by the ejection of single stars that interact strongly with the binary and hence its amplitude depends on the mass ratio $`m_{}/M_{12}`$. Fig. 4 suggests that the wandering amplitude is constant or even increasing with time. The $`N`$-dependence of the wandering amplitude is another important topic for future study.

The kinematical structure of the merged galaxy at late times is shown in Fig. 5. The rotation curve reflects primarily the initial orbital angular momentum of the galaxies and not the presence of the BHs; it peaks at a radius of $``$ a few times $`r_b`$, similar to what is seen in merger simulations without BHs (e.g. Bendo & Barnes 2000) and in real, high-luminosity ellipticals (e.g. Statler et al. 1999). It is very different from the rotation curves seen in fainter galaxies (e.g. M32, Joseph et al. 2000) which remain flat or rising into the sphere of influence of the BH. The velocity dispersions do continue to rise toward the center, roughly as $`r^{1/2}`$, as they did in the initial galaxy models.

This work was supported by NASA through grant NAG5-6037 and by the NSF through grant AST 96-17088. We thank S. Aarseth and R. Spurzem for their advice and help with NBODY6++.

## References

Barnes, J. E. 1999, in Galaxy Interactions at Low and High Redshifts, IAU Symp. 186, ed. J. E. Barnes & D. B. Sanders (Dordrecht: Kluwer), 137

Begelman, M. C., Blandford, R. D. & Rees, M. J. 1980, Nature, 287, 307

Bendo, G. J. & Barnes, J. 2000, astro-ph/0003248

Dehnen, W. 1993, MNRAS, 265, 250

Ebisuzaki, T., Makino, J. & Okumura, S. K. 1991, Nature, 354, 212

Faber, S. M. et al. 1997, AJ, 114, 1771

Ferrarese, L., van den Bosch, F. C., Ford, H. C. & O’Connell, R. W. 1994, AJ, 108, 1598

Ferrarese, L. & Merritt, D. 2000, ApJ, in press

Gebhardt, K. et al. 1996, AJ, 112, 105

Gebhardt, K. et al. 2000a, AJ, 119, 1157

Gebhardt, K. et al. 2000b, ApJ, in press

Gould, A. & Rix, H.-W. 2000, ApJ, 532, L29

Governato, F., Colpi, M. & Maraschi, L. 1994, MNRAS, 271, 317

Hemsendorf, M. 1999, PhD Thesis, Universität Heidelberg

Hills, J. G. 1992, AJ, 103, 1955

Holley-Bockelmann & Richstone 1999, ApJ, 517, 92

Holley-Bockelmann & Richstone 2000, ApJ, 531, 232

Joseph, C. et al. 2000, astro-ph/0005530

Macchetto, F. et al. 1997, ApJ, 489, 579

Makino, J. 1997, ApJ, 478, 58

Makino, J. & Ebisuzaki, T. 1996, ApJ, 465, 527

Makino, J., Fukushige, T., Okumura, S. K. & Ebisuzaki, T. 1993, PASP, 45, 303

Merritt, D. 2000, in Dynamics of Galaxies: From the Early Universe to the Present, ASP Conf. Ser. Vol. 197, ed. F. Combes, G. A. Mamon & V. Charmandaris, 221

Merritt, D. & Fridman, T. 1995, in Fresh Views of Elliptical Galaxies, A.S.P. Conf. Ser. Vol. 86 (ed. A. Buzzoni, A. Renzini & A. Serrano), 22

Miller, B. W. et al. 2000, in Imaging the Universe in Three Dimensions, ASP Conf. Ser. Vol. 195 (ed. W. van Breugel & J. Bland-Hawthorn), 158

Nakano, T. & Makino, J. 1999a, ApJ, 510, 155

Nakano, T. & Makino, J. 1999b, astro-ph/9906131

Peebles, P. J. E. 1972, Gen. Rel. Grav. 3, 63

Quinlan, G. D., 1996, NewA, 1, 35

Quinlan, G. D. & Hernquist, L. 1997, New A, 2, 533

Quinlan, G. D., Hernquist, L. & Sigurdsson, S. 1995, ApJ, 440, 554

Springel, V., Yoshida, N. & White, S. D. M. 2000, astro-ph/0003162

Spurzem, R. & Baumgardt, H. 1999, “A parallel implementation of an Aarseth $`N`$-body integrator on general and special purpose supercomputers,” preprint

Statler, T. S. & Smecker-Hane, T. 1999, AJ, 117, 839

Tremaine, S. et al. 1994, AJ, 107, 634

Valtonen, M. 1996, Comm. Ap., 18, 191

van der Marel, R. 1994, MNRAS, 270, 271

van der Marel, R. 1999, AJ, 117, 744


text:
# Reply to “Comment on ‘Critical behavior of a two-species reaction-diffusion problem’ ”

## Acknowledgments

The authors acknowledge discussion with F. van Wijland. This work is part of the French-Brazilian scientific cooperation project CAPES-COFECUB 229/97. The authors also thank CNPq and Projeto Nordeste de Pesquisa for support.