format: text

text:
# Electron Mass Operator in a Strong Magnetic Field and Dynamical Chiral Symmetry Breaking

(April 2, 2002)

## Abstract

The electron mass operator in a strong magnetic field is calculated. The contribution of higher Landau levels of virtual electrons, along with the ground Landau level, is shown to be essential in the leading log approximation. The effect of the electron dynamical mass generation by a magnetic field is investigated. In a model with $`N`$ charged fermions, it is shown that some critical number $`N_{cr}`$ exists for any value of the electromagnetic coupling constant $`\alpha `$, such that the fermion dynamical mass is generated with a doublet splitting for $`N<N_{cr}`$, and the dynamical mass does not arise at all for $`N>N_{cr}`$, thus leaving the chiral symmetry unbroken.

preprint: YARU-HE-02/03
hep-ph/0204201

Asymptotic properties of the QED diagrams and operators in strong magnetic fields $`BB_e,B_e=m_e^2/e4.4110^{13}\text{G}`$ ($`e`$ is the elementary charge) are of conceptual interest both from the standpoint of the searches of the perturbation theory applicability borders, and also in view of possible applications in astrophysics and in cosmology of the early Universe. The investigations of this type are being performed by many authors during a rather long time. For example, a history of calculations of the electron mass operator in a strong magnetic field lasts more than 30 years already. However, as we show in this Letter, it is too early to put the final point in the problem.

One-loop contribution into the electron mass operator in a strong magnetic field was obtained for the first time by Jancovici Jancovici (1969) in the leading (double) log approximation. Later on, in the papers by Loskutov and Skobelev Loskutov and Skobelev (1979, 1982) the attempts were performed to calculate the two-loop contribution and to summarize all the many-loop contributions in the same approximation. A correct formula for the electron mass operator in this approximation was obtained recently in the paper Gusynin and Smilga (1999). However, the double log approximation becomes invalid in asymptotically strong magnetic fields Loskutov and Skobelev (1983), because of the crucial influence of the strong magnetic field on the virtual photon polarization operator. This influence provides an appearance of the effective photon mass, $`m_\gamma ^2=(2\alpha /\pi )eB`$, which replaces the electron mass in one of the two logarithms. A correct expression for the electron mass operator in the leading (single) log approximation was obtained recently in our paper Kuznetsov et al. (2002) by the summation of the rainbow Feynman diagrams, in the form:

$$M=\frac{m_0}{1(\alpha _R/2\pi )[\mathrm{ln}(\pi /\alpha _R)\gamma _\text{E}]\mathrm{ln}(eB/m_0^2)},$$
(1)
where $`m_0`$ is the electron mass without field, $`\gamma _\text{E}=0.577\mathrm{}`$ is the Euler constant, $`\alpha _R`$ is the electromagnetic coupling constant renormalized by the field

$$\alpha _R=\frac{\alpha }{1(\alpha /3\pi )\mathrm{ln}(eB/m_0^2)}.$$
(2)
However, the formula (1) has a restricted area of applicability with respect to the parameter $`\eta _0=\frac{\alpha }{2\pi }\mathrm{ln}\frac{eB}{m_0^2}`$, because the mass operator (1) tends to infinity at some finite value of this parameter. On the one hand, as it was mentioned in Gusynin and Smilga (1999), this could be the signal of a new physics, the spontaneous breaking of chiral symmetry. On the other hand, it indicates that the formula (1) is not self-consistent at the field values for which $`\eta _01`$. It means that the transitions to the nonperturbative regime over the parameter $`\eta _0`$ is necessary. In terms of the diagram technique it leads to the necessity to consider, along with the rainbow diagrams giving the leading log contributions, also the diagrams of a next-to-leading order. As the analysis shows, this extension reduces to a substitution of the electron mass operator instead of mass into the diagrams. It is more convenient, however, to use the Schwinger-Dyson equations as it was performed in Refs. Gusynin et al. (1999a, b); Alexandre et al. (2000, 2001).

In this Letter, we realize this program by calculating the electron mass operator to be valid in any asymptotically strong magnetic fields. To obtain a full system of the Schwinger-Dyson equations in the presence of external magnetic field, which corresponds to the irreducible diagrams depicted in Figs. 1 and 2, one should know the exact vertex $`\mathrm{\Gamma }_\mu `$ which contains, as is well-known, an infinite number of irreducible diagrams.

Fortunately, the problem is simplified essentially in the strong field limit, as it was shown for the first time in Ref. Loskutov and Skobelev (1984), see also Gusynin et al. (1999a, b). Namely, the exact vertex $`\mathrm{\Gamma }_\mu `$ can be reduced to the bare vertex $`\gamma _\mu `$ in an appropriate gauge.

In the strong field limit, the main contribution into the electron mass operator originates from virtual electrons occupying the lowest Landau level Loskutov and Skobelev (1976). In this case the Fourier transform of the translational invariant part of the exact electron propagator in a magnetic field can be presented in the form

$$G_0^{(e)}(k)=i\mathrm{exp}\left(\frac{k_{}^2}{2eB}\right)\left[\widehat{k}_{}\widehat{M}(k_{})\right]^1O^{()},$$
(3)
where $`\widehat{k}_{}=k^0\gamma ^0k^3\gamma ^3`$ (the field is directed along the 3d axis), $`\widehat{M}(k_{})`$ is the sought mass operator, $`O^{()}=(1i\gamma ^1\gamma ^2)/2`$ is the projecting operator corresponding to the electron state with the spin antiparallel to the external field direction.

On the other hand, photons of the only one mode (’transversal’ in the Adler’s notation Adler (1971)) interact with electrons occupying the lowest Landau level. In this case, in the gauge where $`\mathrm{\Gamma }_\mu \gamma _\mu `$ the exact photon propagator can be presented in the form Loskutov and Skobelev (1983)

$`G_{\mu \nu }^{(\gamma )}(q)=i𝒟(q_{}^2,q_{}^2)\stackrel{~}{\mathrm{\Lambda }}_{\mu \nu },`$ (4)
$`𝒟(q_{}^2,q_{}^2)={\displaystyle \frac{1}{q^2𝒫(q_{}^2,q_{}^2)}}.`$ (5)
Here $`q_{}^2=q_\mu \mathrm{\Lambda }_{\mu \nu }q_\nu `$, $`q_{}^2=q_\mu \stackrel{~}{\mathrm{\Lambda }}_{\mu \nu }q_\nu `$, $`q^2=q_{}^2q_{}^2`$, $`\mathrm{\Lambda }_{\mu \nu }=\phi _{\mu \rho }\phi _{\rho \nu }`$, $`\stackrel{~}{\mathrm{\Lambda }}_{\mu \nu }=\stackrel{~}{\phi }_{\mu \rho }\stackrel{~}{\phi }_{\rho \nu }`$, $`\phi _{\alpha \beta }`$ is the dimensionless tensor of external magnetic field, $`\phi _{\alpha \beta }=F_{\alpha \beta }/\sqrt{F_{\mu \nu }^2/2}`$, $`\stackrel{~}{\phi }_{\alpha \beta }=\frac{1}{2}\epsilon _{\alpha \beta \rho \sigma }\phi _{\rho \sigma }`$ is the dual tensor. The function $`𝒫(q_{}^2,q_{}^2)`$ is the eigenvalue of the photon polarization operator $`𝒫_{\mu \nu }(q)`$, which is depicted by the Feynman diagram Fig. 2. With the reduction $`\mathrm{\Gamma }_\mu =\gamma _\mu `$, the operator has the form

$`𝒫_{\mu \nu }(q)`$ $`=`$ $`i{\displaystyle \frac{\alpha }{4\pi ^3}}{\displaystyle d^4kTr\left[\gamma _\mu G^{(e)}(k)\gamma _\nu G^{(e)}(kq)\right]}`$ (6)
$`=`$ $`\left(\stackrel{~}{\mathrm{\Lambda }}_{\mu \nu }{\displaystyle \frac{q_\mu q_\nu }{q_{}^2}}\right)𝒫(q_{}^2,q_{}^2)+\mathrm{},`$
where dots denote the contribution of the other photon modes. Thus the polarization operator (6) is reduced in fact to the one-loop operator.

The mass operator $`\widehat{M}(p_{})`$ corresponding to the irreducible diagram depicted in Fig. 1, in the gauge $`\mathrm{\Gamma }_\mu \gamma _\mu `$ and in terms of Eqs. (3) and (4) is reduced to a function $`M(p_{}^2)`$, see Refs. Gusynin et al. (1999a, b). The following integral equation arises for this function

$$M(p_{}^2)=m_0i\frac{\alpha }{2\pi ^3}d^4k\mathrm{exp}\left(\frac{k_{}^2}{2eB}\right)\frac{M(k_{}^2)}{k_{}^2[M(k_{}^2)]^2}\frac{1}{(kp)_{}^2k_{}^2𝒫(k_{}^2,(kp)_{}^2)}.$$
(7)
Integral in the expression (7) does not contain an ultraviolet divergency, because the integration over the momenta $`k_{}`$ transversal to the field direction has the cutoff $`k_{}\sqrt{eB}`$.

On the other hand, the photon polarization operator (6), in general, does contain the ultraviolet divergency. As a result, virtual electrons occupying both the lowest and higher Landau levels contribute to the integral. This fact was not taken into account in all previous publications in the field. However, as will be shown below, it leads to very interesting physical consequences.

The function $`𝒫(q_{}^2,q_{}^2)`$ in a strong magnetic field in the one-loop approximation can be extracted, for example, from Ref. Tsai (1974), see also Shabad (1992), where the sought function $`M(q_{}^2)`$ should be inserted instead of the field-free mass $`m_0`$:

$`𝒫(q_{}^2,q_{}^2)`$ $`=`$ $`{\displaystyle \frac{2\alpha }{\pi }}eB\mathrm{exp}\left({\displaystyle \frac{q_{}^2}{2eB}}\right)H\left({\displaystyle \frac{q_{}^2}{4[M(q_{}^2)]^2}}\right)`$ (8)
$`+`$ $`{\displaystyle \frac{\alpha }{3\pi }}q^2\mathrm{ln}{\displaystyle \frac{eB}{[M(q_{}^2)]^2}}.`$
Here the first term is caused by virtual electrons occupying the lowest Landau level, while the second term contains the contribution from higher Landau levels. The function $`H(z)`$ has the form

$$H(z)=\frac{1}{2\sqrt{z(1z)}}\mathrm{ln}\frac{\sqrt{1z}+\sqrt{z}}{\sqrt{1z}\sqrt{z}}1.$$
(9)
In the analysis of Eq. (7), we are interested in the region of parameters $`q_{}^2<0`$, $`|q_{}^2|[M(q_{}^2)]^2`$. For large negative values of the argument, the function $`H(z)`$ is simplified, $`H(z)1`$. The first term in Eq. (8) acquires in this case the meaning of the photon mass squared, $`m_\gamma ^2=(2\alpha /\pi )eB`$, induced by a magnetic field. As for the second term in Eq. (8), containing the contribution from higher Landau levels into the photon polarization operator, its role is in renormalization of the electromagnetic constant $`\alpha `$ in a magnetic field, $`\alpha \alpha _R`$. The expression for the renormalized constant can be obtained from Eq. (2) by the replacement $`m_0M(q_{}^2)`$.

Turning back to Eq. (7), it should be mentioned that the integral in the right-hand side exactly corresponds to the one-loop field-induced correction to the electron mass, with the replacement of the field-free mass $`m_0`$ by the sought mass operator $`M(k_{}^2)`$ in the integrand. In the superstrong field limit, the main contribution into the integral (7) in the form of a big logarithm $`\mathrm{ln}(eB/M^2)`$ originates from the region $`[M(k_{}^2)]^2|k_{}^2|m_\gamma ^2\alpha eB`$, $`m_\gamma ^2k_{}^2eB`$. In view of this, the calculation of the integral in Eq. (7) with logarithmic accuracy gives the following result for the mass operator $`M(p_{}^2)`$ in the Euclidean region $`p_{}^2<0,|p_{}^2|m_\gamma ^2`$:

$$M(p_{}^2)M(0)\left[1+\frac{p_{}^2}{4eB}\mathrm{ln}\frac{2eB}{[M(0)]^2}\right].$$
(10)
The expression (10) shows that the electron physical mass which is defined, strictly speaking, as the solution of the dispersion equation $`m=M(m^2)`$, can be taken with a great accuracy in the zero point, $`mM(0)`$.

In the leading log approximation, we have obtained the following transcendental equation for the electron physical mass from the integral equation (7):

$$m=m_0+m\frac{\alpha _R}{2\pi }\left(\mathrm{ln}\frac{\pi }{\alpha _R}\gamma _\text{E}\right)\mathrm{ln}\frac{eB}{m^2},$$
(11)
where $`\alpha _R`$ is taken in the point $`q_{}^2=0`$, i.e. obtained from Eq. (2) by the replacement $`m_0M(0)=m`$. It is interesting to note that the formula (11) reproduces exactly our result (1), with the substitution of the electron physical mass $`m`$ instead of the field-free mass $`m_0`$ under the logarithms. The equation (11) solves the problem of finding the electron physical mass for any large values of the magnetic field. It is free of a singularity, unlike the Eq. (1).

An analysis of Eq. (11) shows that its solution in asymptotically strong fields when $`mm_0`$, becomes independent on $`m_0`$ and is reduced in fact to the solution at $`m_0=0`$. This would mean the generation of the dynamical mass of the initially massless electron in a magnetic field. This effect which is also called the dynamical chiral symmetry breaking, was studied in refs. Gusynin and Smilga (1999); Gusynin et al. (1999a, b); Alexandre et al. (2000, 2001), however, the contribution from higher Landau levels into the photon polarization operator was not considered there.

Let us show, that this contribution changes the behaviour of the dynamical mass essentially. Let us extend our analysis to a model with $`N`$ initially massless fermions ($`m_0=0`$) with equal charges $`e`$ (in a case of different fermion charges $`Q_fe`$, the parameter $`N`$ has the meaning of $`_fQ_f^2`$). In this case the photon polarization operator (8) is the sum over all fermion loops, i.e. it should be multiplied by $`N`$. The transcendental equation for the fermion dynamical mass ($`m0`$) can be obtained from Eq. (11) as follows

$$\frac{\alpha _R}{2\pi }\left(\mathrm{ln}\frac{\pi }{N\alpha _R}\gamma _\text{E}\right)\mathrm{ln}\frac{eB}{m^2}=1,$$
(12)
where

$$\alpha _R=\frac{\alpha }{1(N\alpha /3\pi )\mathrm{ln}(eB/m^2)}.$$
(13)
The expression (12) allows to reproduce the result of Refs. Gusynin et al. (1999a, b), if the actual dependence of the coupling constant $`\alpha _R`$ on the ratio $`eB/m^2`$ is formally ignored and $`\alpha _R=\alpha `$ is taken in it. It is remarkable that the constant $`C_1`$ obtained there by a numerical calculation as $`C_11.82\pm 0.06`$, appears to be $`C_1=\pi \mathrm{exp}(\gamma _\text{E})=1.763877\mathrm{}`$.

In Fig. 3 the behaviour of the fermion dynamical mass $`m`$ divided by $`\sqrt{eB}`$ is shown versus the number of fermions $`N`$ and the field-free coupling constant $`\alpha `$, considered as free parameters of the model.

The dependence is seen to differ essentially from the results of Refs. Gusynin et al. (1999a, b). Namely, for any value of the coupling constant $`\alpha `$, such a critical number of fermions $`N_{cr}`$ exists that for $`N<N_{cr}`$ two values of the fermion dynamical mass are generated. For $`N>N_{cr}`$ the equation (12) does not have a solution at all, thus the chiral symmetry is kept unbroken.

The dependence of the critical number $`N_{cr}`$ on the value of the coupling constant $`\alpha `$ takes the form

$`N_{cr}(\alpha )`$ $`=`$ $`\sqrt{{\displaystyle \frac{3\pi }{2\mathrm{exp}(\gamma _\text{E}+1)}}}{\displaystyle \frac{1}{\sqrt{\alpha }}}{\displaystyle \frac{3}{4}}`$ (14)
$`+`$ $`{\displaystyle \frac{9}{16}}\sqrt{{\displaystyle \frac{3\mathrm{exp}(\gamma _\text{E}+1)}{2\pi }}}\sqrt{\alpha }+O(\alpha ).`$
It is quite remarkable that the doublet splitting of the fermion dynamical mass can be rather large. For example, for $`\alpha =0.1`$ and $`N=1`$ the mass difference is of 15 orders of magnitude. If one considers for the purposes of illustration the magnetic field value $`10^{33}\text{G}`$ Ambjørn and Olesen (1993), the two values of the fermion dynamical mass are $`10^2\text{GeV}`$ and $`10^4\text{eV}`$.

We believe that the doublet splitting of the fermion dynamical mass and the conservation of the chiral symmetry at $`N>N_{cr}`$ as well are new interesting physical phenomena in QED in strong external magnetic field.

###### Acknowledgements.

We are grateful to V.A. Rubakov and M. V. Chistyakov for helpful discussions. This work was supported in part by the Russian Foundation for Basic Research under the Grant No. 01-02-17334 and by the Ministry of Education of Russian Federation under the Grant No. E00-11.0-5.


text:
# Multiple lepton pair production in relativistic ion collisions

## 1 Introduction

The multiplicity and the distribution of lepton pairs produced in the Coulomb fields of two colliding relativistic heavy ions are closely connected to the problem of unitarity. When heavy ions collide at relativistic velocities their Lorentz contracted electromagnetic fields are sufficiently intense to produce a large numbers of such pairs. Usually the process of lepton pairs production is considered as pair creation in the classic Coulomb potential of a charge moving along a straight line. Such an approach allows one  to investigate the impact parameter dependent total probability of the pair creation $`P(b)`$, which by definition is connected with the total cross section $`\sigma =P(b)d^2b`$. As was noticed in  the probability of single pair production calculated to lowest order in fine structure constant $`\alpha =\frac{e^2}{4\pi }=\frac{1}{137}`$ (throughout we use units for which $`c=\mathrm{}=1`$) at small impact parameters exceeds one, thus violating unitarity. This excess begins at impact parameters smaller than the Compton wavelength of the electron $`\lambda _c=\frac{1}{m}=386`$fm and at energies of practical interest (RHIC & LHC).

Allowance for the finite size of colliding nuclei doesn’t remedy the situation, because that would affect only the impact parameters comparable to t0he nuclei radii, which are much smaller than the Compton wavelength of the electron $`R\lambda _c`$.

As was shown in  this problem can be solved by taking into account the possibility of multiple pair production, whose relative contribution grows with energy and dominates at small impact parameters. Since this early publication there has been much work has in this direction (see e.g.  and references therein) with the common for all statement: the probability to produce n lepton pairs in the Coulomb field of heavy ions colliding at fixed impact parameter $`b`$ can be approximately represented as a Poisson distribution i.e. $`P(b,n)=\frac{W^n(b)}{n!}e^{W(b)}`$, where W(b) is the average multiplicity of pairs at a fixed impact parameter.

Because of the somewhat controversial situation in the subject and its importance for the operation of relativistic heavy ion colliders (RHIC & LHC), in the present communication we revisit the multiple pair production based on the powerful Sudakov technique of calculation of high energy Feynman diagrams. Recently it has been applied  to the calculation of Coulomb corrections to single lepton pair production in relativistic heavy ion collisions. Here we extend the Sudakov approach to get the probability of multiple pair ($`n2`$) production in relativistic heavy ions collisions. For heavy ions with charge the numbers satisfy the conditions $`Z_1\alpha Z_2\alpha 1`$, $`Z_1Z_2\alpha 1`$, one needs full allowance of multiple Coulomb interaction of colliding nuclei, whereas the secondary interaction of produced pairs (real or virtual) with the Coulomb fields of colliding ions can be neglected.

The paper is organized as follows. In section 2 we introduce the necessary definitions and discussed briefly the Sudakov’derivation of the Born amplitude for single pair creation. Then We show how one can sum the perturbative series to get the amplitude for the process of $`n`$ pairs production in a compact analytical form.

In Section 3 we consider the Coulomb interaction between the colliding ions and show how it affects the amplitude and cross section of the process under consideration. Apart from the conventional multiple photon exchange between ions, this interaction include also the so called ”screening corrections” which are the result of the ions interaction through the virtual lepton pairs and are responsible for the probability of vacuum–to–vacuum transition. In the language of the Abramovski–Gribov–Kancheli unitarity rules , the diagrams with lepton pair loops can be associated with the QED model for the pomeron, and the n-pair production amplitudes can be associated with the unitarity cuts through n exchanged pomerons. In section 4 we report the multiple production amplitudes in the impact representation in which the summation of the perturbation series for the n-pair production can be done in a compact closed form. Finally, using this amplitude we obtain the probability of $`n`$ pair production which can be cast in the Poisson form, and the total cross section of all pair production process.

## 2 The amplitude of the $`n`$ pairs production process.

The typical Feynman diagram (FD) describing the $`n`$ lepton pairs production in the collision of relativistic nuclei with atomic numbers $`A_1,A_2`$, with $`n_e`$ exchanged photons between colliding nuclei as well as screening effects e.g. the insertions of $`n_s`$ light by light (LBL) scattering blocks is drawn in Fig. 1.

Upper and lower blocks in Fig. 1 describes many virtual photons interaction amplitudes with nuclei. They contain the complete set of $`(n_e+n+2n_s)!`$ Feynman diagrams. To avoid the multiple counting in what follows, we will multiply the relevant amplitude by the factor $`1/(n_e!n_s!(2!)^{2n_s})`$.

As was mentioned above we restrict ourselves to ions with charge numbers such that

$$Z_{1,2}1,Z_1\alpha Z_2\alpha 1,Z_1Z_2\alpha >1,$$
(1)
which permits us to omit the multiphoton exchanges between the produced pairs and colliding ions. The dominant mechanism is a production of a single lepton pair per collision of equivalent photons, i.e., one lepton loop per two-photon ladder. The alternative mechanism of multiple pair production per collision of equivalent photons, i.e., multiple lepton loops per two-photon ladder, is suppressed by inverse powers of $`Z_1Z_2`$.

For the description of a peripheral process of n lepton pairs creation i.e. the process

$$\begin{array}{c}A_1(Z_1,p_1)+A_2(Z_2,p_2)A_1(Z_1,p_1^{})+A_2(Z_2,p_2^{})+e_+e_{}(r_1)+\mathrm{}+e_+e_{}(r_n),\hfill \\ \hfill r_i=q_+^i+q_{}^i\end{array}$$
(2)
It is convenient to use the Sudakov parameterization for the four-momentum of all exchanged photons (for details see )

$$k_i=\alpha _i\stackrel{~}{p}_2+\beta _i\stackrel{~}{p}_1+k_i,d^4k_i=\frac{s}{2}d\alpha _id\beta _id^2k_i,$$
(3)
$$s=(p_1+p_2)^2,sp_i^2=M_i^2m^2,\stackrel{~}{p}_1=p_1p_2\frac{p_1^2}{s},\stackrel{~}{p}_2=p_2p_1\frac{p_2^2}{s},$$

$$\stackrel{~}{p}_1^2=\stackrel{~}{p}_2^2=O\left(\frac{m^6}{s^2}\right),\stackrel{~}{p}_1k_i=\stackrel{~}{p}_2k_i=0,s=2p_1p_2=2\stackrel{~}{p}_1\stackrel{~}{p}_2.$$

Here $`\stackrel{~}{p}_i`$ are light-like four-vectors build from $`p_i`$, $`M_i`$ are the masses of colliding nuclei, $`m`$ and $`s`$ are the electron mass and the total center mass energy.

The denominators of intermediate states of the nucleon Green functions for upper and lower blocks are the same and have the following form

$$s\underset{i}{}\alpha _i\left(\underset{i}{}𝐤_𝐢\right)^2+i0,s\underset{i}{}\beta _i\left(\underset{i}{}𝐤_𝐢\right)^2+i0.$$
(4)
Peripheral process is characterized by small values of longitudinal Sudakov parameters $`\alpha _i`$, $`\beta _i`$ and the transverse momenta of the order of electron mass

$$|\alpha _i||\beta _i|1,k_i^2m^2.$$
(5)
Further simplification follows from the form of the nominators of the exchanged photons Green function (we work in the Feynman gauge). Using the Gribov’s representation for the metric tensors

$$g_{\mu \nu }=g_{\mu \nu }^{}+\frac{2}{s}\left(\stackrel{~}{p_1}^\mu \stackrel{~}{p_2}^\nu +\stackrel{~}{p_1}^\nu \stackrel{~}{p_2}^\mu \right).$$
(6)
it is easy to show that for the typical conversion of the nuclei currents $`J_\mu (p_1)J^\mu (p_2)`$ only one term, which contains the scalar products of a nucleus current with a four–momentum of another nucleus, becomes relevant (with power accuracy)

$$J_\mu (p_1)J^\mu (p_2)\frac{2}{s}J_\lambda (p_1)p_2^\lambda J_\sigma (p_2)p_1^\sigma \left(1+O\left(\frac{m^2}{s}\right)\right).$$
(7)
It can be seen that the quantity $`J_\mu (p_1)J^\mu (p_2)/s`$ remains finite with the large values of $`s`$. This fact provides great simplification of the spinor structure of the amplitude

$$\overline{u}(p_1^{})\stackrel{~}{p}_2(p_1+\chi _1+M_1)\stackrel{~}{p}_2\mathrm{}(p_2+\chi _N+M_1)\stackrel{~}{p}_2u(p_1)s^{N+1}N_1,$$
(8)
$$\overline{u}(p_2^{})\stackrel{~}{p}_1(p_2+\eta _1+M_2)\stackrel{~}{p}_1\mathrm{}(p_2+\eta _N+M_2)\stackrel{~}{p}_1u(p_2)s^{N+1}N_2,$$

$$N_1=\frac{1}{s}\overline{u}(p_1^{})\widehat{p}_2u(p_1),N_2=\frac{1}{s}\overline{u}(p_2^{})\widehat{p}_1u(p_2).$$

Besides we have $`|N_1|^2=|N_2|^2=2`$ for the nuclei with the spin 1/2 and $`|N_1|^2=|N_2|^2=1`$ for the scalar one. Using the identity

$$\underset{perm}{}\frac{1}{\alpha _{i_1}}\frac{1}{\alpha _{i_1}+\alpha _{i_2}}\mathrm{}\frac{1}{{\displaystyle \underset{j=1}{\overset{N}{}}}\alpha _{i_j}}=\underset{i=1}{\overset{N}{}}\frac{1}{\alpha _i}$$
(9)
one can be convinced that the amplitude describing the upper and lower blocks in Fig. 1 can be put in the form

$`I_1=N_1{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{s}{s\alpha _i+i0}}+{\displaystyle \frac{s}{s\alpha _i+i0}}\right),I_2=N_2{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{s}{s\beta _i+i0}}+{\displaystyle \frac{s}{s\beta _i+i0}}\right)`$ (10)
with $`N=n_e+2n_s+n1`$.

This expressions contain all dependence on Sudakov parameters $`\alpha _i`$, $`\beta _i`$ (the 4-momenta of exchanged photons in the peripheral kinematics in denominators of their Green functions can be considered as Euclidean two-vectors $`k_i^2=s\alpha _i\beta _i+k_i^2k_i^2=𝐤_{𝐢}^{}{}_{}{}^{2}`$).

At this stage the integration over Sudakov parameters can be done, because the dependence of the amplitude on $`\alpha _i`$, $`\beta _i`$ provides the convergence of the relevant integrals

$$I_1\underset{i=1}{\overset{N}{}}d\alpha _i=(2\pi i)^NN_1,I_2\underset{i=1}{\overset{N}{}}d\beta _i=(2\pi i)^NN_2.$$
(11)
Let us now consider the single pair production. The amplitude of the process (1) in its lowest order (Born approximation) reads

$$M_{(0)}^{(1)}=is(8\pi \alpha )^2Z_1Z_2N_1N_2\frac{B_{\alpha \beta }p_1^\alpha p_2^\beta }{s𝐪_\mathrm{𝟏}{}_{}{}^{2}𝐪_{\mathrm{𝟐}}^{}^2},B_{\alpha \beta }=\overline{v}(q_+)O_{\alpha \beta }u(q_{}),$$
(12)
$`B_{\alpha \beta }`$ is the Compton tensor  for pair creation by two virtual photons with polarization vectors $`e_1(q_1)`$, $`e_2(q_2)`$. $`q_{1(2)}`$ are the 4–momenta of exchange photons and $`r=q_++q_{}`$. Strictly speaking the squares of these 4–vectors $`q_i^2`$ do not vanish in the limit $`𝐪_𝐢{}_{}{}^{2}0`$. This fact becomes essential when one calculates the total cross section of a single pair production process. For the case of two or more pairs production (which is our case) the replacement $`q_i^2=𝐪_{𝐢}^{}{}_{}{}^{2}`$ can safely be done.

Using the gauge invariance

$$q_1^\alpha B_{\alpha \beta }=q_2^\beta B_{\alpha \beta }=0$$
(13)
one can perform the replacement

$$\frac{B_{\alpha \beta }p_1^\alpha p_2^\beta }{s}=\frac{B_{\alpha \beta }e_1^\alpha e_2^\beta }{\stackrel{~}{s}_1}|𝐪_\mathrm{𝟏}||𝐪_\mathrm{𝟐}|,e_i^\alpha =\frac{q_i^\alpha }{|𝐪_𝐢|},\stackrel{~}{s}_1=s\alpha _2\beta _1,$$
(14)
The quantity $`\stackrel{~}{s}_1`$ is related to the square of the invariant mass of a pair

$$s_1=(q_1+q_2)^2=\stackrel{~}{s}_1(𝐪_\mathrm{𝟏}+𝐪_\mathrm{𝟐})^2=(q_++q_{})^2.$$
(15)
Two-dimensional vectors $`e_i`$ can be interpreted as a polarization vectors of exchanged virtual photons.

Using (13-15) one can rewrite the Born amplitude (12) in a form

$$M_{(0)}^{(1)}=isN_1N_2B(q_1,q_2),$$
(16a)
$$B(q_1,q_2)=(8\pi \alpha )^2Z_1Z_2\frac{B_{\alpha \beta }e_1^\alpha e_2^\beta }{\stackrel{~}{s}_1|𝐪_\mathrm{𝟏}||𝐪_\mathrm{𝟐}|}$$
(16b)
Now we are able to construct the amplitude for the process of $`n`$ pairs production. Bearing in mind the expressions (11) the matrix element of two pairs production can be represented as the convolution of two Born terms from expression (16b)

$$M_{(0)}^{(2)}=i^2sN_1N_2B(k,kr_1)B(qk,qr_2k)\frac{d^2𝐤}{8\pi ^2}.$$
(17)
A straightforward generalization to the matrix element in the case of $`n`$ pairs production reads

$$\begin{array}{c}M_{(0)}^{(n)}=i^nsN_1N_2\underset{i=1}{\overset{n1}{}}\left(B(k_i,k_ir_i)\frac{d^2k_i}{8\pi ^2}\right)B(h,hr_n),\hfill \\ \hfill h=q\underset{i=1}{\overset{n1}{}}k_i.\end{array}$$
(18)
Thus one can see that the amplitude for multiple pair production is solely determined by the convolution of the amplitudes corresponding to single pair production.

In the language of the AGK unitarity rules, this result can be interpreted as unitarity cut through all exchanged pomerons.

## 3 The Coulomb exchanges between ions

Let us now consider the effect of $`m`$ photon exchanges between nuclei $`A_1`$, $`A_2`$. The arguments given above leads to the following matrix element for the process of $`n`$ pairs production with m photons exchanges among the colliding ions

$`M_{(m)}^{(n)}`$ $`={\displaystyle \frac{i^nsN_1N_2}{m!}}{\displaystyle \underset{j=1}{\overset{m}{}}\left(\frac{i\alpha Z_1Z_2}{𝝌_j^2+\lambda ^2}\frac{d^2𝝌_j}{\pi }\right)}`$ (19)
$`{\displaystyle }{\displaystyle \underset{i=1}{\overset{n1}{}}}\left(B(k_i,k_ir_i){\displaystyle \frac{d^2𝐤_𝐢}{8\pi ^2}}\right)B(k_n,k_nr_n),k_n=q{\displaystyle \underset{i=1}{\overset{n1}{}}}k_i{\displaystyle \underset{i=1}{\overset{m}{}}}\chi _i.`$
Another effect which we take into account is the possibility of the ion–ion interaction through the LBL blocks (screening effect in Fig. 1), which in the AGK language  is equivalent to the exchange by additional un–cut pomerons. It is associated with the iteration of a typical kernel

$$LY(𝐥)d^2𝐥,$$
(20a)
$$Y(𝐥)=\frac{(\alpha ^2Z_1Z_2)^2}{32\pi ^4}\frac{P}{|𝐥_\mathrm{𝟏}||𝐥_\mathrm{𝟐}||𝐥_\mathrm{𝟏}𝐥||𝐥_\mathrm{𝟐}+𝐥|}d^2𝐥_\mathrm{𝟏}d^2𝐥_\mathrm{𝟐}\frac{d\stackrel{~}{s_1}}{\stackrel{~}{s_1}^2},$$
(20b)
$$P=\mathrm{\Pi }^{\alpha \beta \gamma \delta }e_1^\alpha (l_1)e_2^\beta (l_1l)e_3^\gamma (l_2)e_4^\delta (l_2+l),$$
(20c)
where we rearrange the ”extra” phase volume of longitudinal Sudakov parameters $`\alpha _2`$, $`\beta _1`$ in terms of invariant mass square of LBL block and extract explicitly the boost degree of freedom of LBL block

$$\frac{d\alpha _2d\beta _1}{s(\alpha _2\beta _1)^2}=\frac{d\beta _1}{\beta _1}\frac{d\stackrel{~}{s}_1}{\stackrel{~}{s}_1^2}=L\frac{d\stackrel{~}{s}_1}{\stackrel{~}{s}_1^2},$$
(21)
$$L=\mathrm{ln}\gamma _1\gamma _2,\stackrel{~}{s}_1=(l_1+l_2)^2>4m^2.$$

The structure (20) entered the matrix element (19) in the compact form

$$\frac{1}{n_s!}\underset{i=1}{\overset{n_s}{}}\left(LY(𝐥_𝐢)d^2𝐥_𝐢\right)$$
(22)
Real part of (20c) (equal to one half of its s–channel discontinuity) is related with (16b)

$$\begin{array}{c}\mathrm{}\mathrm{\Pi }^{\alpha \beta \gamma \delta }e_1^\alpha (l_1)e_2^\beta (l_1l)e_3^\gamma (l_2)e_4^\delta (l_2+l)=\frac{1}{2}B_{\alpha \beta }(l_1,rl_1)e_1^\alpha (l_1)e_2^\beta (rl_1)\hfill \\ \hfill B_{\gamma \delta }(l_1l,r+ll_1)e_1^\gamma (l_1l)e_2^\delta (r+ll_1)d\mathrm{\Phi }_r,\end{array}$$
(23)
where $`d\mathrm{\Phi }_r`$ is the phase volume of the intermediate pair

$$d\mathrm{\Phi }_r=\frac{\delta ^4(rq_+q_{})d^3q_+d^3q_{}}{(2\pi )^22\epsilon _+2\epsilon _{}}.$$
(24)
We will show later that the imaginary part of $`P`$ is irrelevant either for the total cross section or for the probability of $`n`$ pairs production distribution.

## 4 The impact parameter representation for the amplitude

The last step in building the matrix element for the process of $`n`$ real pairs creation consist in transformation of the obtained above expressions in the impact–parameter representation and summation over all eikonal photons and LBL blocks. For this we introduce the identity

$$\begin{array}{c}\delta ^2\left(𝐤_𝐧𝐪+\underset{i=1}{\overset{n1}{}}𝐤_𝐢+\underset{i=1}{\overset{n_s}{}}\chi _𝐢+\underset{i=1}{\overset{n_v}{}}𝐥_𝐦\right)d^2𝐤_𝐧=\hfill \\ \hfill \frac{1}{4}e^{i𝐪𝝆}\mathrm{exp}\left[i𝝆(𝐤_𝐧+\underset{i=1}{\overset{n1}{}}𝐤_𝐢+\underset{i=1}{\overset{n_s}{}}\chi _𝐢+\underset{i=1}{\overset{n_v}{}}𝐥_𝐦)\right]\frac{d^2𝐤_𝐧}{\pi }\frac{d^2𝝆}{\pi }=1\end{array}$$
(25)
Using this expression the summation in $`n_e`$ and $`n_s`$ can be easily done with the result

$$M^{(n)}=\frac{i^n\pi s}{2}N_1N_2e^{i𝐪𝝆}e^{i\mathrm{\Psi }(\rho \lambda )}e^{L[A(\rho )/2+i\phi (\rho )]}\underset{i=1}{\overset{n}{}}\stackrel{~}{B}(𝝆,𝐫_𝐢)\frac{d^2𝝆}{\pi },n2,$$
(26)
with

$$\frac{A(\rho )}{2}+i\phi (\rho )=Y(𝐥)e^{i𝐥𝝆}\frac{d^2𝐥}{\pi },\stackrel{~}{B}(𝝆,𝐫_𝐢)=B(𝐤,𝐤𝐫_𝐢)e^{i𝐤𝝆}\frac{d^2𝐤}{8\pi ^2},$$
(27)
$$\mathrm{\Psi }(\rho \lambda )=\alpha Z_1Z_2\frac{e^{i𝝌𝝆}}{𝝌^2+\lambda ^2}\frac{d^2𝝌}{\pi }=2\alpha Z_1Z_2K_0(\lambda \rho )$$
(28)
where $`K_0(\lambda \rho )`$ is modified Bessel function (Mac-Donald function).

The phase volume of final state which consists from the scattered nuclei and $`n`$ pairs can be written in the following form

$`d\mathrm{\Gamma }_{n+2}`$ $`={\displaystyle \underset{i=1}{\overset{n}{}}}\left({\displaystyle \frac{d^3q_+d^3q_{}}{(2\pi )^62\epsilon _+2\epsilon _{}}}\right){\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle \frac{d^3p_1^{}d^3p_2^{}}{2\epsilon _1^{}2\epsilon _2^{}}}\delta ^4\left(p_1+p_2p_1^{}p_2^{}{\displaystyle \underset{i=1}{\overset{n}{}}}r_i\right)`$
$`={\displaystyle \underset{i=1}{\overset{n}{}}}\left({\displaystyle \frac{L}{(2\pi )^4}}{\displaystyle \frac{d^2𝐫_𝐢}{2}}ds_id\mathrm{\Phi }_i\right){\displaystyle \frac{d^2𝐪}{2s(2\pi )^2}}`$ (29)
with $`q=p_1^{}`$.

Cross section of $`n`$ pairs production has the form

$$d\sigma _n=\frac{1}{8s}\frac{|M_n|^2}{n!}d\mathrm{\Gamma }_{n+2}.$$
(30)
Statistical factor $`1/n!`$ is included to take into account the identity of pairs. Using the expressions (25)-(4) we get

$$\frac{d\sigma _n}{d^2\rho }=P_n(\rho ),P_n(\rho )=\frac{(LA_1(\rho ))^n}{n!}e^{LA_1(\rho )},n2$$
(31)
with

$$A_1(\rho )=\frac{1}{2^5\pi ^4}|B(\rho ,r)|^2𝑑s_1d^2𝐫𝑑\mathrm{\Phi }_r.$$
(32)
It can be easily recognized that for $`A(\rho )`$ from (27)

$$A(\rho )=A_1(\rho ),$$
(33)
thus confirming the Poisson character of probability distribution in impact-parameter representation.

Let us mention that the effect of eikonal photons as well as the imaginary part of the amplitude corresponding to LBL blocks do not modify the total cross section as well as differential cross section integrated over phase volume of final nuclei. Really, integrating the square of the amplitude (26) over the phase volume one immediately obtains

$$\begin{array}{c}e^{i𝐪(𝝆_1𝝆_2)}e^{i\{[\psi (\rho _1)\psi (\rho _2)]L[\phi (\rho _1)\phi (\rho _2)]\}}f(\rho _1)f(\rho _2)\frac{d^2𝐪}{\pi }\frac{d^2𝝆_1}{\pi }\frac{d^2𝝆_2}{\pi }=\hfill \\ \hfill f^2(\rho )\frac{d^2𝝆}{\pi }\end{array}$$
(34)
Nevertheless, the exclusive cross section is sensitive to both these factors.

Expression (32) can be simplified if one neglects the dependence of Compton tensor (14) on external photons virtualities

$$B_{\alpha \beta }(k,r)e^\alpha e^\beta B_{\alpha \beta }(0,r)e^\alpha e^\beta $$
(35)
and use the well known relation

$$\underset{4m^2}{\overset{\mathrm{}}{}}|B_{\alpha \beta }(0,r)e^\alpha e^\beta |^2\frac{ds_1d\mathrm{\Phi }_r}{s_1^2}=\frac{7}{36\pi m^2}$$
(36)
As a result the expression (32) can be cast in the form

$$A(\rho )=\frac{7}{18\pi ^2m^2}(\alpha ^2Z_1Z_2)^2I(\rho ),$$
(37)
$`I(\rho )`$ $`={\displaystyle \stackrel{m}{}}{\displaystyle \frac{e^{i𝝆(𝐤_\mathrm{𝟏}𝐤_\mathrm{𝟐})}}{|𝐤_\mathrm{𝟏}||𝐤_\mathrm{𝟏}𝐫||𝐤_\mathrm{𝟐}||𝐤_\mathrm{𝟐}𝐫|}}{\displaystyle \frac{d^2𝐫}{\pi }}{\displaystyle \frac{d^2𝐤_\mathrm{𝟏}}{\pi }}{\displaystyle \frac{d^2𝐤_\mathrm{𝟐}}{\pi }}`$
$`={\displaystyle \stackrel{m}{}}{\displaystyle \frac{e^{i𝝌𝝆}}{|𝐤||𝐤𝝌||𝐤^{}||𝐤^{}+𝝌|}}{\displaystyle \frac{d^2𝝌}{\pi }}{\displaystyle \frac{d^2𝐤}{\pi }}{\displaystyle \frac{d^2𝐤^{}}{\pi }}`$
$`={\displaystyle e^{i𝝌𝝆}\mathrm{ln}^2\left(\frac{m^2}{𝝌^2}\right)\frac{d^2𝝌}{\pi }}{\displaystyle \frac{16}{\rho ^2}}\left(\mathrm{ln}(\rho m)+O(1)\right),`$
where we introduce the cut–off parameter $`|k|<m`$ as a result of the fast decreasing of matrix element of pair production by two photons. For the case of heavy leptons production ($`\mu `$ or $`\tau `$) the upper limit must be replaced by quantity $`Q`$ which can be associated with maximal momentum transferred to nucleus without its disintegration. For the case $`\rho m1`$ one has

$$A(\rho )\frac{56}{9}\frac{(\alpha ^2Z_1Z_2)^2}{\pi ^2(\rho m)^2}\left(\mathrm{ln}(\rho m)+O(1)\right)$$
(38)
which is in the agreement with .

Our formula (31) can be applied for the case $`n=0`$ (the probability of elastic nuclei scattering). The case $`n=1`$ needs a bit accurate consideration. The expression for $`\sigma _1`$ can be cast in the following form

$$\sigma _1=\sigma _B+LA(\rho )\left(e^{LA(\rho )}1\right)d^2𝝆.$$
(39)
First term $`\sigma _B`$ corresponds to the leading order of the Racah formula (see for instance ) for the cross section in Born approximation. Inferring (39), one has to take into account the longitudinal components of momenta of exchanged photons which create the pair. The second term takes into account the unitarity corrections to the total cross section.

## Conclusions

In conclusions let us enumerate very shortly the main results of the present work. In our paper we obtained the general form for the amplitude of $`n`$ lepton pairs production, accounting for the mutual Coulomb interaction of relativistic ions. We get the probability of vacuum–vacuum transition in the closed analytic form. We confirmed that the probability of $`n`$ pair production can be approximated by Poisson distribution. Although this result has been claimed before in a number of publications, the mutual interaction of ions has not been taken into account in these works. At last we have shown that the Coulomb interaction between the colliding nuclei, although important in the exclusive kinematics, does not affect the integral yield of lepton pairs.

## Acknowledgements

We are grateful to Valery Serbo, Alexander Tarasov and Lev Lipatov for their interest to the subject and useful comments. The work is supported by INTAS 97-30494, INTAS-00366 and SR-2000 grants. E. K. is grateful to NCPI, JINR for support.


text:
# Advantages of modified ADM formulation: constraint propagation analysis of Baumgarte-Shapiro-Shibata-Nakamura system

## I Introduction

One of the current most important topics in the field of numerical relativity is to find a formulation of the Einstein equations which gives us stable and accurate longterm evolution. We all know that simulating space-time and matter based on general relativity is the essential research direction to go in the future, but we do not have a definite recipe for controlling numerical blow-ups. We concentrate our discussion on free evolution of the Einstein equations based on the 3+1 (space + time) decomposition of space-time, which requires solving the constraints only on the initial hypersurface and monitors the violation (error) of the calculation by checking constraints during the evolution.

Over the decades, the Arnowitt-Deser-Misner (ADM) adm  formulation has been treated as the default for numerical relativists. (More precisely, the version introduced by Smarr and York ADM-York  was taken as the default, which we denote the standard ADM formulation hereafter.) Although the ADM formulation mostly works for gravitational collapses or cosmological models in numerical treatments, it does not satisfy the requirement for longterm evolution e.g. the studies of gravitational wave sources.

As we mentioned in our previous paper adjADMsch , we think we can classify the current efforts of formulating equations for numerical relativity in the following three ways: (1) apply a modified ADM (BSSN) formulation SN ; BS , (2) apply a first-order hyperbolic formulation (see the references e.g. in ronbun1 ; reviewLehner ; KST ), or (3) apply an asymptotically constrained system BFHR ; SiebelHuebner ; SY-asympAsh ; ronbun2 .

The first refers to using a modified ADM formulation, originally proposed by Nakamura in late 80s, and subsequently modified by Nakamura-Oohara and Shibata-Nakamura SN . This introduces conformal decomposition of the ADM variables, a new variable for calculating Ricci curvature, and adjusts the equations of motion using constraints. The advantage of this formulation was re-introduced by Baumgarte and Shapiro BS , and therefore this is often cited as the BSSN formulation, which we follow also. The BSSN equations are now widely used in the large scale numerical computations, including coalescence of binary neutron stars binaryNS  and binary black holes binaryBH .

The second and third efforts use similar modifications such as introductions of new variables and/or adjustments of the equations, but differ in their purposes: to construct a hyperbolic formulation or to construct a formulation which constraints will decay or propagate away. The latter is intended to control numerical evolution so as the constrained manifold is its attractor. While the hyperbolic formulations have been extensively studied in this direction, we think the worrisome point in the discussion is the treatment of the non-principal part which is ignored in the hyperbolic formulation. As Kidder, Scheel and Teukolsky KST  reported recently, unless we reduce the effect of the non-principal part of the equations we may not get an advantages of hyperbolic formulation for numerical results ronbun1 ; HernPHD .

Through the series of studies adjADMsch ; ronbun1 ; ronbun2 ; adjADM , we propose a systematic treatment for constructing a robust evolution system against perturbative error. We call it an asymptotically constrained (or asymptotically stable) system if the error decays itself. The idea is to adjust evolution equations using constraints (we term it an adjusted system), and to decide the coefficients (multipliers) by analyzing constraint propagation equations. We propose to apply an eigenvalue analysis of the propagation equations of the constraints, especially in its Fourier components so as to include the non-principal part in the analysis. The characters of eigenvalues will be changed according to the adjustments to the original evolution equations. We conjectured that the constraint violation which was occurred during the evolution will decay (if negative real eigenvalues) or propagate away (if pure imaginary eigenvalues).

This conjecture was affirmatively confirmed to explain the numerical behaviors: wave propagation in the Maxwell equations ronbun2 , in the Ashtekar version of the Einstein equations ronbun2 , and in the ADM formulation (flat space-time background) adjADM . The advantage of this construction scheme is that it can be applied to a formulation which is not a first-order hyperbolic form, such as to the ADM formulation adjADM ; adjADMsch . We think therefore our proposal is an alternative way to control/predict the violation of constraints. (We believe that the idea of the constraint propagation analysis first appeared in Frittelli Fri-con , where she made hyperbolicity classification for the standard ADM formulation).

The purpose of this article is to apply this constraint propagation analysis to the BSSN formulation, and understand how each improvement contributes to more stable numerical evolution. Together with numerical comparisons with the standard ADM case potsdam0003 ; LHG , this topic has been studied by many groups with different approaches. Using numerical test evolutios, Alcubierre et al. potsdam9908  found that the essential improvement is in the process of replacing terms by constraints, and that the eigenvalues of the BSSN evolution equations has fewer “zero eigenvalues” than those of ADM, and they conjectured that the instability can be caused by “zero eigenvalues” that violate “gauge mode”. Miller Miller  applied von Neumann’s stability analysis to the plane wave propagation, and reported that BSSN has a wider range of parameters that give us stable evolution. These studies provide some supports regarding the advantage of BSSN, while it was also shown an example of an ill-posed solution in BSSN (as well in ADM) FrittelliGomez . (Inspired by BSSN’s conformal decomposition, several related hyperbolic formulations have also been proposed ArBona ; FR99 ; ABMS .)

We think our analysis will offer a new vantage point on the topic, and contribute an alternative understanding of its background. Consequently, we propose more effective improvement of the BSSN system which has not yet been tried in numerical simulations.

The construction of this paper is as follows. We review the BSSN system in §II, and there also we discuss where the adjustments are applied. In §III we apply our constraint propagation analysis to show how each improvement works in the BSSN equations, and in §IV we extend our study to seek a better formulation which might be obtained by small steps. We only consider the vacuum space-time thoughout the article, but the inclusion of matter is straightforward.

## II BSSN equations and their constraint propagation equations

### II.1 BSSN equations

We start presenting the standard ADM formulation, which expresses the space-time with a pair of 3-metric $`\gamma _{ij}`$ and extrinsic curvature $`K_{ij}`$. The evolution equations become

$`_t^A\gamma _{ij}`$ $`=`$ $`2\alpha K_{ij}+D_i\beta _j+D_j\beta _i,`$ (1)
$`_t^AK_{ij}`$ $`=`$ $`\alpha R_{ij}^{ADM}+\alpha KK_{ij}2\alpha K_{ik}K_{}^{k}{}_{j}{}^{}D_iD_j\alpha `$ (2)
$`+(D_i\beta ^k)K_{kj}+(D_j\beta ^k)K_{ki}+\beta ^kD_kK_{ij}`$
where $`\alpha ,\beta _i`$ are the lapse and shift function and $`D_i`$ is the covariant derivative on 3-space. The symbol $`_t^A`$ means the time derivative defined by these equations, and we distinguish them from those of the BSSN equations $`_t^B`$, which will be defined in (15)-(19). The associated constraints are the Hamiltonian constraint $``$ and the momentum constraints $`_i`$:

$`^{ADM}`$ $`=`$ $`R^{ADM}+K^2K_{ij}K^{ij},`$ (3)
$`_i^{ADM}`$ $`=`$ $`D_jK^j{}_{i}{}^{}D_iK.`$ (4)
The widely used notationSN ; BS  is to introduce the variables ($`\phi ,\stackrel{~}{\gamma }_{ij}`$,$`K`$,$`\stackrel{~}{A}_{ij}`$,$`\stackrel{~}{\mathrm{\Gamma }}^i`$) instead of ($`\gamma _{ij}`$,$`K_{ij}`$), where

$`\phi `$ $`=`$ $`(1/12)\mathrm{log}(\mathrm{det}\gamma _{ij}),`$ (5)
$`\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`e^{4\phi }\gamma _{ij},`$ (6)
$`K`$ $`=`$ $`\gamma ^{ij}K_{ij},`$ (7)
$`\stackrel{~}{A}_{ij}`$ $`=`$ $`e^{4\phi }(K_{ij}(1/3)\gamma _{ij}K),`$ (8)
$`\stackrel{~}{\mathrm{\Gamma }}^i`$ $`=`$ $`\stackrel{~}{\mathrm{\Gamma }}_{jk}^i\stackrel{~}{\gamma }^{jk}.`$ (9)
The new variable $`\stackrel{~}{\mathrm{\Gamma }}^i`$ was introduced in order to calculate Ricci curvature more accurately. $`\stackrel{~}{\mathrm{\Gamma }}^i`$ also contributes to make the system re-produce wave equations in its linear limit. In the BSSN formulation, Ricci curvature is not calculated as

$$R_{ij}^{ADM}=_k\mathrm{\Gamma }_{ij}^k_i\mathrm{\Gamma }_{kj}^k+\mathrm{\Gamma }_{ij}^l\mathrm{\Gamma }_{lk}^k\mathrm{\Gamma }_{kj}^l\mathrm{\Gamma }_{li}^k,$$
(10)
but

$`R_{ij}^{BSSN}`$ $`=`$ $`R_{ij}^\phi +\stackrel{~}{R}_{ij},`$ (11)
$`R_{ij}^\phi `$ $`=`$ $`2\stackrel{~}{D}_i\stackrel{~}{D}_j\phi 2\stackrel{~}{\gamma }_{ij}\stackrel{~}{D}^k\stackrel{~}{D}_k\phi `$ (12)
$`+4(\stackrel{~}{D}_i\phi )(\stackrel{~}{D}_j\phi )4\stackrel{~}{\gamma }_{ij}(\stackrel{~}{D}^k\phi )(\stackrel{~}{D}_k\phi ),`$
$`\stackrel{~}{R}_{ij}`$ $`=`$ $`(1/2)\stackrel{~}{\gamma }^{lk}_l_k\stackrel{~}{\gamma }_{ij}+\stackrel{~}{\gamma }_{k(i}_{j)}\stackrel{~}{\mathrm{\Gamma }}^k+\stackrel{~}{\mathrm{\Gamma }}^k\stackrel{~}{\mathrm{\Gamma }}_{(ij)k}`$ (13)
$`+2\stackrel{~}{\gamma }^{lm}\stackrel{~}{\mathrm{\Gamma }}_{l(i}^k\stackrel{~}{\mathrm{\Gamma }}_{j)km}+\stackrel{~}{\gamma }^{lm}\stackrel{~}{\mathrm{\Gamma }}_{im}^k\stackrel{~}{\mathrm{\Gamma }}_{klj},`$
where $`\stackrel{~}{D}_i`$ is covariant derivative associated with $`\stackrel{~}{\gamma }_{ij}`$. These are weakly equivalent, but $`R_{ij}^{BSSN}`$ does have wave operator apparently in the flat background limit, so that we can expect more natural wave propagation behavior.

Additionally, BSSN requires us to impose the conformal factor as

$$\stackrel{~}{\gamma }(:=\mathrm{det}\stackrel{~}{\gamma }_{ij})=1,$$
(14)
during the evolution. This is a kind of definition, but can also be thought of as a constraint. We will return to this point shortly.

BSSN’s improvements are not only the introductions of new variables, but also the replacement of terms in the evolution equations using the constraints. The purpose of this article is to understand and to identify which improvement works for the stability. Before doing that we first show the standard set of the BSSN evolution equations:

$`_t^B\phi `$ $`=`$ $`(1/6)\alpha K+(1/6)\beta ^i(_i\phi )+(_i\beta ^i),`$ (15)
$`_t^B\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`2\alpha \stackrel{~}{A}_{ij}+\stackrel{~}{\gamma }_{ik}(_j\beta ^k)+\stackrel{~}{\gamma }_{jk}(_i\beta ^k)(2/3)\stackrel{~}{\gamma }_{ij}(_k\beta ^k)+\beta ^k(_k\stackrel{~}{\gamma }_{ij}),`$ (16)
$`_t^BK`$ $`=`$ $`D^iD_i\alpha +\alpha \stackrel{~}{A}_{ij}\stackrel{~}{A}^{ij}+(1/3)\alpha K^2+\beta ^i(_iK),`$ (17)
$`_t^B\stackrel{~}{A}_{ij}`$ $`=`$ $`e^{4\phi }(D_iD_j\alpha )^{TF}+e^{4\phi }\alpha (R_{ij}^{BSSN})^{TF}+\alpha K\stackrel{~}{A}_{ij}2\alpha \stackrel{~}{A}_{ik}\stackrel{~}{A}^k{}_{j}{}^{}+(_i\beta ^k)\stackrel{~}{A}_{kj}+(_j\beta ^k)\stackrel{~}{A}_{ki}`$ (18)
$`(2/3)(_k\beta ^k)\stackrel{~}{A}_{ij}+\beta ^k(_k\stackrel{~}{A}_{ij}),`$
$`_t^B\stackrel{~}{\mathrm{\Gamma }}^i`$ $`=`$ $`2(_j\alpha )\stackrel{~}{A}^{ij}+2\alpha (\stackrel{~}{\mathrm{\Gamma }}_{jk}^i\stackrel{~}{A}^{kj}(2/3)\stackrel{~}{\gamma }^{ij}(_jK)+6\stackrel{~}{A}^{ij}(_j\phi ))_j(\beta ^k(_k\stackrel{~}{\gamma }^{ij})\stackrel{~}{\gamma }^{kj}(_k\beta ^i)`$ (19)
$`\stackrel{~}{\gamma }^{ki}(_k\beta ^j)+(2/3)\stackrel{~}{\gamma }^{ij}(_k\beta ^k)).`$
We next summarize the constraints in this system. The normal Hamiltonian and momentum constraints (the “kinematic” constraints) are naturally written as

$`^{BSSN}`$ $`=`$ $`R^{BSSN}+K^2K_{ij}K^{ij},`$ (20)
$`_i^{BSSN}`$ $`=`$ $`_i^{ADM},`$ (21)
where we use Ricci scalar defined by (11). Additionally, we regard the following three as the constraints (the “algebraic” constraints):

$`𝒢^i`$ $`=`$ $`\stackrel{~}{\mathrm{\Gamma }}^i\stackrel{~}{\gamma }^{jk}\stackrel{~}{\mathrm{\Gamma }}_{jk}^i,`$ (22)
$`𝒜`$ $`=`$ $`\stackrel{~}{A}_{ij}\stackrel{~}{\gamma }^{ij},`$ (23)
$`𝒮`$ $`=`$ $`\stackrel{~}{\gamma }1,`$ (24)
where the first two are from the algebraic definition of the variables (8) and (9), and the (24) is from the requirement of (14). Hereafter we write $`^{BSSN}`$ and $`^{BSSN}`$ simply as $``$ and $``$ respectively.

Taking careful account of these constraints, (20) and (21) can be expressed directly as

$``$ $`=`$ $`e^{4\phi }\stackrel{~}{R}8e^{4\phi }\stackrel{~}{D}^j\stackrel{~}{D}_j\phi 8e^{4\phi }(\stackrel{~}{D}^j\phi )(\stackrel{~}{D}_j\phi )+(2/3)K^2\stackrel{~}{A}_{ij}\stackrel{~}{A}^{ij}(2/3)𝒜K,`$ (25)
$`_i`$ $`=`$ $`6\stackrel{~}{A}^j{}_{i}{}^{}(\stackrel{~}{D}_j\phi )2𝒜(\stackrel{~}{D}_i\phi )(2/3)(\stackrel{~}{D}_iK)+\stackrel{~}{\gamma }^{kj}(\stackrel{~}{D}_j\stackrel{~}{A}_{ki}).`$ (26)
In summary, the fundamental dynamical variables in BSSN are ($`\phi ,\stackrel{~}{\gamma }_{ij}`$, $`K`$,$`\stackrel{~}{A}_{ij}`$,$`\stackrel{~}{\mathrm{\Gamma }}^i`$), total 17. The gauge quantities are ($`\alpha ,\beta ^i`$) which is 4, and the constraints are $`(,_i,𝒢^i,𝒜,𝒮)`$, i.e. 9 components. As a result, 4 (2 by 2) components are left which correspond to two gravitational polarization modes.

### II.2 Adjustments in evolution equations

Next, we show the BSSN evolution equations (15)-(19) again, identifying where the terms are replaced using the constraints, (20)-(24).

By a straightforward calculation, we get:

$`_t^B\phi `$ $`=`$ $`_t^A\phi +(1/6)\alpha 𝒜(1/12)\stackrel{~}{\gamma }^1(_j𝒮)\beta ^j,`$ (27)
$`_t^B\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`_t^A\stackrel{~}{\gamma }_{ij}(2/3)\alpha \stackrel{~}{\gamma }_{ij}𝒜+(1/3)\stackrel{~}{\gamma }^1(_k𝒮)\beta ^k\stackrel{~}{\gamma }_{ij},`$ (28)
$`_t^BK`$ $`=`$ $`_t^AK(2/3)\alpha K𝒜\alpha +\alpha e^{4\phi }(\stackrel{~}{D}_j𝒢^j),`$ (29)
$`_t^B\stackrel{~}{A}_{ij}`$ $`=`$ $`_t^A\stackrel{~}{A}_{ij}+\left((1/3)\alpha \stackrel{~}{\gamma }_{ij}K(2/3)\alpha \stackrel{~}{A}_{ij}\right)𝒜+\left((1/2)\alpha e^{4\phi }(_k\stackrel{~}{\gamma }_{ij})(1/6)\alpha e^{4\phi }\stackrel{~}{\gamma }_{ij}\stackrel{~}{\gamma }^1(_k𝒮)\right)𝒢^k`$ (30)
$`+\alpha e^{4\phi }\stackrel{~}{\gamma }_{k(i}(_{j)}𝒢^k)(1/3)\alpha e^{4\phi }\stackrel{~}{\gamma }_{ij}(_k𝒢^k),`$
$`_t^B\stackrel{~}{\mathrm{\Gamma }}^i`$ $`=`$ $`_t^A\stackrel{~}{\mathrm{\Gamma }}^i+\left((2/3)(_j\alpha )\stackrel{~}{\gamma }^{ji}(2/3)\alpha (_j\stackrel{~}{\gamma }^{ji})(1/3)\alpha \stackrel{~}{\gamma }^{ji}\stackrel{~}{\gamma }^1(_j𝒮)+4\alpha \stackrel{~}{\gamma }^{ij}(_j\phi )\right)𝒜(2/3)\alpha \stackrel{~}{\gamma }^{ji}(_j𝒜)`$ (31)
$`+2\alpha \stackrel{~}{\gamma }^{ij}_j(1/2)(_k\beta ^i)\stackrel{~}{\gamma }^{kj}\stackrel{~}{\gamma }^1(_j𝒮)+(1/6)(_j\beta ^k)\stackrel{~}{\gamma }^{ij}\stackrel{~}{\gamma }^1(_k𝒮)+(1/3)(_k\beta ^k)\stackrel{~}{\gamma }^{ij}\stackrel{~}{\gamma }^1(_j𝒮)`$
$`+(5/6)\beta ^k\stackrel{~}{\gamma }^2\stackrel{~}{\gamma }^{ij}(_k𝒮)(_j𝒮)+(1/2)\beta ^k\stackrel{~}{\gamma }^1(_k\stackrel{~}{\gamma }^{ij})(_j𝒮)+(1/3)\beta ^k\stackrel{~}{\gamma }^1(_j\stackrel{~}{\gamma }^{ji})(_k𝒮).`$
where $`_t^A`$ denotes the part of no replacements, i.e. the terms only use the standard ADM evolution equations in its time derivatives.

From (27)-(31), we understand that all the BSSN evolution equations are adjusted using constraints. This fact will give us the importance of the scaling constraint $`𝒮=0`$ and the tracefree operation $`𝒜=0`$ during the evolution.

As we have pointed out in the case of adjusted ADM systems adjADM ; adjADMsch , certain combinations of adjustments (replacements) in the evolution equations change the eigenvalues of constraint propagation equations drastically. For example, all negative eigenvalues can be negative real by applying Detweiler’s adjustment detweiler  or its simplified version. One common fact we found is that such a case has an adjustment which breaks time reversal parity of the original equation. That is, with a change of time integration direction $`_t_t`$, an adjusted term might become effective if it breaks time reversal symmetry. (This time asymmtric feature was first impremented as a “lambda-system” in BFHR .) Unfortunately, for the case of the BSSN equations, (27)-(31), all the above adjustments keep the time reversal symmetry. So that we can not expect direct decays of constraint violation in the present form. We will give the details on this point later.

## III Constraint propagation analysis in flat space-time

### III.1 Procedures

We start this section overviewing the procedures and our goals. In our series of previous works ronbun2 ; adjADM ; adjADMsch , we have concluded that eigenvalue analysis of the constraint propagation equations are quite useful for explaining or predicting how the constraint violation grows.

Suppose we have a set of dynamical variables $`u^a(x^i,t)`$, and their evolution equations

$$_tu^a=f(u^a,_iu^a,\mathrm{}),$$
(32)
and the (first class) constraints

$$C^\alpha (u^a,_iu^a,\mathrm{})0.$$
(33)
For monitoring the violation of constraints, we propose to investigate the evolution equations of $`C^\alpha `$ (constraint propagation),

$$_tC^\alpha =g(C^\alpha ,_iC^\alpha ,\mathrm{}).$$
(34)
(We do not mean to integrate (34) numerically, but rather to evaluate them analytically in advance.) In order to analyze the contributions of all RHS terms in (34), we propose to reduce (34) in ordinary differential equations by Fourier transformation,

$`_t\widehat{C}^\alpha `$ $`=`$ $`\widehat{g}(\widehat{C}^\alpha )=M^\alpha {}_{\beta }{}^{}\widehat{C}_{}^{\beta },`$ (35)
where $`C(x,t)^\rho =\widehat{C}(k,t)^\rho \mathrm{exp}(ikx)d^3k`$, and then to analyze the set of eigenvalues, say $`\mathrm{\Lambda }_\alpha `$, of the coefficient matrix, $`M^\alpha _\beta `$, in (35). We call $`\mathrm{\Lambda }`$s and $`M^\alpha _\beta `$ the constraint amplification factors (CAFs) of (34) and constraint propagation matrix, respectively. Our guidelines to have ‘better stability’ are that

* If the CAFs have a negative real-part  (the constraints are forced to be diminished), then we see more stable evolution than a system which has a positive real-part.
* If the CAFs have a non-zero imaginary-part  (the constraints are propagating away), then we see more stable evolution than a system which has zero CAFs.

We found heuristically that the system becomes more stable when more $`\mathrm{\Lambda }`$s satisfy the above criteria ronbun1 ; ronbun2 . We note that these guidelines are confirmed numerically for wave propagation in the Maxwell system and in the Ashtekar version of the Einstein system ronbun2 , and also for error propagation in Minkowskii space-time using adjusted ADM systems adjADM . Supporting theorems for above (A) was recently discussed diagCP .

The above features of the constraint propagation, (34), will differ when we modify the original evolution equations. Suppose we add (adjust) the evolution equations using constraints

$$_tu^a=f(u^a,_iu^a,\mathrm{})+F(C^\alpha ,_iC^\alpha ,\mathrm{}),$$
(36)
then (34) will also be modified as

$$_tC^\alpha =g(C^\alpha ,_iC^\alpha ,\mathrm{})+G(C^\alpha ,_iC^\alpha ,\mathrm{}).$$
(37)
Therefore, the problem is how to adjust the evolution equations so that their constraint propagation satisfies the above criteria as much as possible.

### III.2 BSSN constraint propagation equations

Our purpose in this section is to apply the above procedure to the BSSN system. The set of the constraint propagation equations, $`_t(,_i,𝒢^i,𝒜,𝒮)^T`$, turns to be quite long and not elegant (is not a first-order hyperbolic and includes many non-linear terms), and we put them in Appendix. In order to understand the fundamental structure, we hereby show an analysis on the flat space-time background.

For the flat background metric $`g_{\mu \nu }=\eta _{\mu \nu }`$, the first order perturbation equations of (27)-(31) can be written as

$`_t{}_{}{}^{(1)}\varphi `$ $`=`$ $`(1/6){}_{}{}^{(1)}K+(1/6)(\kappa _\varphi 1){}_{}{}^{(1)}𝒜,`$ (38)
$`_t{}_{}{}^{(1)}\stackrel{~}{\gamma }_{ij}^{}`$ $`=`$ $`2{}_{}{}^{(1)}\stackrel{~}{A}_{ij}^{}(2/3)(\kappa _{\stackrel{~}{\gamma }}1)\delta _{ij}{}_{}{}^{(1)}𝒜,`$ (39)
$`_t{}_{}{}^{(1)}K`$ $`=`$ $`(_j_j{}_{}{}^{(1)}\alpha )+(\kappa _{K1}1)_j{}_{}{}^{(1)}𝒢_{}^{j}(\kappa _{K2}1){}_{}{}^{(1)},`$ (40)
$`_t{}_{}{}^{(1)}\stackrel{~}{A}_{ij}^{}`$ $`=`$ $`{}_{}{}^{(1)}(R_{ij}^{BSSN}){}_{}{}^{TF}{}_{}{}^{(1)}(\stackrel{~}{D}_i\stackrel{~}{D}_j\alpha ){}_{}{}^{TF}+(\kappa _{A1}1)\delta _{k(i}(_{j)}{}_{}{}^{(1)}𝒢_{}^{k})(1/3)(\kappa _{A2}1)\delta _{ij}(_k{}_{}{}^{(1)}𝒢_{}^{k}),`$ (41)
$`_t{}_{}{}^{(1)}\stackrel{~}{\mathrm{\Gamma }}_{}^{i}`$ $`=`$ $`(4/3)(_i{}_{}{}^{(1)}K)(2/3)(\kappa _{\stackrel{~}{\mathrm{\Gamma }}1}1)(_i{}_{}{}^{(1)}𝒜)+2(\kappa _{\stackrel{~}{\mathrm{\Gamma }}2}1){}_{}{}^{(1)}_{i}^{},`$ (42)
where we introduced parameters $`\kappa `$s, all $`\kappa =0`$ reproduce no adjustment case from the standard ADM equations, and all $`\kappa =1`$ correspond to the BSSN equations. We express them as

$`\kappa _{adj}`$ $`:=`$ $`(\kappa _\phi ,\kappa _{\stackrel{~}{\gamma }},\kappa _{K1},\kappa _{K2},\kappa _{A1},\kappa _{A2},\kappa _{\stackrel{~}{\mathrm{\Gamma }}1},\kappa _{\stackrel{~}{\mathrm{\Gamma }}2}).`$ (43)
Constraint propagation equations at the first order in the flat space-time, then, become:

$`_t{}_{}{}^{(1)}`$ $`=`$ $`\left(\kappa _{\stackrel{~}{\gamma }}(2/3)\kappa _{\stackrel{~}{\mathrm{\Gamma }}1}(4/3)\kappa _\phi +2\right)_j_j{}_{}{}^{(1)}𝒜+2(\kappa _{\stackrel{~}{\mathrm{\Gamma }}2}1)(_j{}_{}{}^{(1)}_{j}^{}),`$ (44)
$`_t{}_{}{}^{(1)}_{i}^{}`$ $`=`$ $`\left((2/3)\kappa _{K1}+(1/2)\kappa _{A1}(1/3)\kappa _{A2}+(1/2)\right)_i_j{}_{}{}^{(1)}𝒢_{}^{j}`$ (45)
$`+(1/2)\kappa _{A1}_j_j{}_{}{}^{(1)}𝒢_{}^{i}+\left((2/3)\kappa _{K2}(1/2)\right)_i{}_{}{}^{(1)},`$
$`_t{}_{}{}^{(1)}𝒢_{}^{i}`$ $`=`$ $`2\kappa _{\stackrel{~}{\mathrm{\Gamma }}2}{}_{}{}^{(1)}_{i}^{}+((2/3)\kappa _{\stackrel{~}{\mathrm{\Gamma }}1}(1/3)\kappa _{\stackrel{~}{\gamma }})(_i{}_{}{}^{(1)}𝒜),`$ (46)
$`_t{}_{}{}^{(1)}𝒮`$ $`=`$ $`2\kappa _{\stackrel{~}{\gamma }}{}_{}{}^{(1)}𝒜,`$ (47)
$`_t{}_{}{}^{(1)}𝒜`$ $`=`$ $`(\kappa _{A1}\kappa _{A2})(_j{}_{}{}^{(1)}𝒢_{}^{j}).`$ (48)
We will discuss CAFs of (44)-(48).

### III.3 Effect of adjustments

We check CAFs of the BSSN equations in detail. The list of examples is shown also in Table 1. Hereafter we let $`k^2=k_x^2+k_y^2+k_z^2`$ for Fourier wave numbers.

1. The no-adjustment case, $`\kappa _{adj}=`$(all zeros). This is the starting point of the discussion. In this case,

$$CAFs=(0(\times 7),\pm \sqrt{k^2}),$$

i.e., $`(0(\times 7),\pm \text{pure imaginary}\text{(1 pair)})`$. In the standard ADM formulation, which uses $`(\gamma _{ij},K_{ij})`$, CAFs are $`(0,0,\pm \text{Pure Imaginary})`$ adjADM . Therefore if we do not apply adjustments in the BSSN equations the constraint propagation structure is quite similar to that of the standard ADM.
2. For the BSSN equations, $`\kappa _{adj}=`$(all 1s),

$$CAFs=(0(\times 3),\pm \sqrt{k^2}\text{(3 pairs)}),$$

i.e., $`(0(\times 3),\pm \text{Pure Imaginary}\text{(3 pairs)})`$. The number of pure imaginary CAFs is increased over that of No.1, and we conclude this is the advantage of adjustments used in the BSSN equations.
3. No $`𝒮`$-adjustment case. All the numerical experiments so far apply the scaling condition $`𝒮`$ for the conformal factor $`\phi `$. The $`𝒮`$-originated terms appear many places in the BSSN equations (15)-(19), so that we suspect non-zero $`𝒮`$ is a kind of source of the constraint violation. However, since all $`𝒮`$-originated terms do not appear in the flat space-time background analysis, \[no adjusted terms in (38)-(42)\], our analysis is independent of the $`𝒮`$-constraint. (Remark that we do not deny the effect of $`𝒮`$-adjustment in other situation.)
4. No $`𝒜`$-adjustment case. The trace (or traceout) condition for the variables is also considered necessary (e.g. AB2001 ). This can be checked with $`\kappa _{adj}=(\kappa ,\kappa ,1,1,1,1,\kappa ,1)`$, and we get

$$CAFs=(0(\times 3),\pm \sqrt{k^2}\text{(3 pairs)}),$$

independent of $`\kappa `$. Therefore the effect of $`𝒜`$-adjustment is unimportant according to this analysis, i.e. on flat space-time background. (Remark that we do not deny the effect of $`𝒜`$-adjustment in other situation.)
5. No $`𝒢^i`$-adjustment case. The introduction of $`\mathrm{\Gamma }^i`$ is the key in the BSSN system. If we do not apply adjustments by $`𝒢^i`$, ($`\kappa _{adj}=(1,1,0,1,0,0,1,1)`$) then we get

$$CAFs=(0(\times 7),\pm \sqrt{k^2}),$$

which is the same with No.1. That is, adjustments due to $`𝒢^i`$ terms are effective to make a progress from ADM.
6. No $`_i`$-adjustment case. This can be checked with $`\kappa _{adj}=(1,1,1,1,1,1,1,\kappa )`$, and we get

$`CAFs`$ $`=`$ $`(0,\pm \sqrt{\kappa k^2}\text{(2 pairs)},`$
$`\pm \sqrt{k^2(1+4\kappa +|14\kappa |)/6},`$
$`\pm \sqrt{k^2(1+4\kappa |14\kappa |)/6}).`$
If $`\kappa =0`$, then $`(0(\times 7),\pm \sqrt{k^2/3})`$, which is $`(0(\times 7),\pm \text{real value})`$. Interestingly, these real values indicate the existence of the error growing mode together with the decaying mode. Alcubierre et al. potsdam9908  found that the adjustment due to the momentum constraint is crucial for obtaining stability. We think that they picked up this error growing mode. Fortunately at the BSSN limit ($`\kappa =1`$), this error growing mode disappears and turns into a propagation mode.
7. No $``$-adjustment case. The set $`\kappa _{adj}=(1,1,1,\kappa ,1,1,1,1)`$ gives

$$CAFs=(0(\times 3),\pm \sqrt{k^2}\text{(3 pairs)}),$$

independently to $`\kappa `$. Therefore the effect of $``$-adjustment is unimportant according to this analysis, i.e. on flat space-time background. (Remark again that we do not deny the effect of $``$-adjustment in other situation.)

These tests are on the effects of adjustments. We will consider whether much better adjustments are possible in the next section.

We list the above results in Table 1. (Table 1 includes a column of diagonalizability of constraint propagation matrix $`M`$, of which importance was pointed out in diagCP . ) The most characteristic points of the above are No. 5 and No.6 that denote the contribution of the momentum constraint adjustment and the importance of the new variable $`\stackrel{~}{\mathrm{\Gamma }}^i`$. It is quite interesting that the unadjusted BSSN equations (case 1) does not have apparent advantages from the ADM system. As we showed in the case 5 and 6, if we missed a particular adjustment, then the expected stability behavior occationally gets worse than the starting ADM system. Therefore we conclude that the better stability of the BSSN formulation is obtained by their adjustments in the equations, and the combination of the adjustments is in a good balance. That is, a lack of their adjustments might fail to obtain the stability of their system.

## IV Proposals of Improved BSSN systems

In this section, we consider the possibility whether we can obtain a system which has much better properties; whether more pure imaginary CAFs or negative real CAFs.

### IV.1 Heuristic examples

(A) A system which has 8 pure imaginary CAFs:
One direction is to seek a set of equations which make fewer zero CAFs than the standard BSSN case (No.2 in the previous section). Using the same set of adjustments in (38)-(42), CAFs are written in general

$`CAFs`$ $`=`$ $`(0,\pm \sqrt{k^2\kappa _{A1}\kappa _{\stackrel{~}{\mathrm{\Gamma }}2}}\text{(2 pairs)},`$
$`\pm \text{complicated expression},`$
$`\pm \text{complicated expression}).`$
The terms in the first line certainly give four pure imaginary CAFs (two positive and negative real pairs) if $`\kappa _{A1}\kappa _{\stackrel{~}{\mathrm{\Gamma }}2}>0(<0)`$. Keeping this in mind, by choosing $`\kappa _{adj}=(1,1,1,1,1,\kappa ,1,1)`$, we find

$`CAFs`$ $`=`$ $`(0,\pm \sqrt{k^2}\text{(2 pairs)},`$
$`\pm \sqrt{k^2(2+\kappa +|\kappa 4|)/6},`$
$`\pm \sqrt{k^2(2+\kappa |\kappa 4|)/6},).`$
Therefore the adjustment $`\kappa _{adj}=(1,1,1,1,1,4,1,1)`$ gives

$$CAFs=(0,\pm \sqrt{k^2}\text{(4 pairs)}),$$

which is one step advanced from BSSN’s according our guidelines.

We note that such a system can be obtained in many ways, e.g. $`\kappa _{adj}=(0,0,1,0,2,1,0,1/2)`$ also gives four pairs of pure imaginary CAFs.

(B) A system which has negative real CAF:
One criterion to obtain a decaying constraint mode (i.e. an asymptotically constrained system) is to adjust an evolution equation as it breaks time reversal symmetry adjADM ; adjADMsch . For example, we consider an additional adjustment to the BSSN equation as

$$_t\stackrel{~}{\gamma }_{ij}=_t^B\stackrel{~}{\gamma }_{ij}+\kappa _{SD}\alpha \stackrel{~}{\gamma }_{ij},$$
(49)
which is a similar adjustment of the simplified Detweiler-type detweiler  that was discussed in adjADMsch . The first order constraint propagation equations on the flat background space-time become

$`_t{}_{}{}^{(1)}`$ $`=`$ $`_j_j{}_{}{}^{(1)}𝒜(3/2)\kappa _{SD}_j_j{}_{}{}^{(1)},`$
$`_t{}_{}{}^{(1)}_{i}^{}`$ $`=`$ $`(1/6)_i{}_{}{}^{(1)}+(1/2)_j_j{}_{}{}^{(1)}𝒢_{}^{i},`$
$`_t{}_{}{}^{(1)}𝒢_{}^{i}`$ $`=`$ $`_i{}_{}{}^{(1)}𝒜+(1/2)\kappa _{SD}_i{}_{}{}^{(1)}+2{}_{}{}^{(1)}_{i}^{},`$
$`_t{}_{}{}^{(1)}𝒜`$ $`=`$ $`(_j_j{}_{}{}^{(1)}\alpha )^{TF}+({}_{}{}^{(1)}R_{jj}^{BSSN})^{TF},`$
$`_t{}_{}{}^{(1)}𝒮`$ $`=`$ $`2{}_{}{}^{(1)}𝒜+3\kappa _{SD}{}_{}{}^{(1)},`$
where we wrote only additional terms to (44)-(48). The CAFs become

$`CAFs`$ $`=`$ $`(0(\times 2),`$
$`\pm \sqrt{k^2}\text{(3 pairs)},(3/2)k^2\kappa _{SD}),`$
in which the last one becomes negative real if $`\kappa _{SD}<0`$.

(C) Combination of above (A) and (B)
Naturally we next consider both adjustments:

$`_t\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`_t^B\stackrel{~}{\gamma }_{ij}+\kappa _{SD}\alpha \stackrel{~}{\gamma }_{ij}`$ (50)
$`_t\stackrel{~}{A}_{ij}`$ $`=`$ $`_t^B\stackrel{~}{A}_{ij}\kappa _8\alpha e^{4\phi }\stackrel{~}{\gamma }_{ij}_k𝒢^k`$ (51)
where the second one produces the 8 pure imaginary CAFs. The additional terms in the constraint propagation equations (44)-(48) are

$`_t{}_{}{}^{(1)}`$ $`=`$ $`_j_j{}_{}{}^{(1)}𝒜(3/2)\kappa _{SD}_j_j{}_{}{}^{(1)},`$
$`_t{}_{}{}^{(1)}_{i}^{}`$ $`=`$ $`(1/6)_i{}_{}{}^{(1)}+(1/2)_j_j{}_{}{}^{(1)}𝒢_{}^{i}`$
$`\kappa _8_i_k{}_{}{}^{(1)}𝒢_{}^{k},`$
$`_t{}_{}{}^{(1)}𝒢_{}^{i}`$ $`=`$ $`_i{}_{}{}^{(1)}𝒜+(1/2)\kappa _{SD}_i{}_{}{}^{(1)}+2{}_{}{}^{(1)}_{i}^{},`$
$`_t{}_{}{}^{(1)}𝒜`$ $`=`$ $`3\kappa _8_k{}_{}{}^{(1)}𝒢_{}^{k},`$
$`_t{}_{}{}^{(1)}𝒮`$ $`=`$ $`2{}_{}{}^{(1)}𝒜+3\kappa _{SD}{}_{}{}^{(1)}.`$
We then obtain

$`CAFs`$ $`=`$ $`(0,\pm \sqrt{k^2}\text{(3 pairs)},`$
$`(3/4)k^2\kappa _{SD}\pm \sqrt{k^2(\kappa _8+(9/16)k^2\kappa _{SD}^2)})`$
which reproduces case (A) when $`\kappa _{SD}=0,\kappa _8=1`$, and case (B) when $`\kappa _8=0`$. These CAFs can become (0, pure imaginary (3 pairs), complex numbers with a negative real part (1 pair)), with an appropriate combination of $`\kappa _8`$ and $`\kappa _{SD}`$.

### IV.2 Possible adjustments

In order to break time reversal symmetry of the evolution equations adjADM ; adjADMsch ; BFHR , the possible simple adjustments are (1) to add $``$, $`𝒮`$ or $`𝒢^i`$ terms to the equations of $`_t\varphi `$, $`_t\stackrel{~}{\gamma }_{ij}`$, or $`_t\stackrel{~}{\mathrm{\Gamma }}^i`$, or (2) to add $`_i`$ or $`𝒜`$ terms to $`_tK`$ or $`_t\stackrel{~}{A}_{ij}`$. We write them generally, including the above proposal (B), as

$`_t\varphi `$ $`=`$ $`_t^B\varphi +\kappa _\varphi \alpha +\kappa _{\varphi 𝒢}\alpha \stackrel{~}{D}_k𝒢^k,`$ (52)
$`_t\stackrel{~}{\gamma }_{ij}`$ $`=`$ $`_t^B\stackrel{~}{\gamma }_{ij}+\kappa _{SD}\alpha \stackrel{~}{\gamma }_{ij}+\kappa _{\stackrel{~}{\gamma }𝒢1}\alpha \stackrel{~}{\gamma }_{ij}\stackrel{~}{D}_k𝒢^k+\kappa _{\stackrel{~}{\gamma }𝒢2}\alpha \stackrel{~}{\gamma }_{k(i}\stackrel{~}{D}_{j)}𝒢^k+\kappa _{\stackrel{~}{\gamma }𝒮1}\alpha \stackrel{~}{\gamma }_{ij}𝒮+\kappa _{\stackrel{~}{\gamma }𝒮2}\alpha \stackrel{~}{D}_i\stackrel{~}{D}_j𝒮,`$ (53)
$`_tK`$ $`=`$ $`_t^BK+\kappa _K\alpha \stackrel{~}{\gamma }^{jk}(\stackrel{~}{D}_j_k),`$ (54)
$`_t\stackrel{~}{A}_{ij}`$ $`=`$ $`_t^B\stackrel{~}{A}_{ij}+\kappa _{A1}\alpha \stackrel{~}{\gamma }_{ij}(\stackrel{~}{D}^k_k)+\kappa _{A2}\alpha (\stackrel{~}{D}_{(i}_{j)})+\kappa _{A𝒜1}\alpha \stackrel{~}{\gamma }_{ij}𝒜+\kappa _{A𝒜2}\alpha \stackrel{~}{D}_i\stackrel{~}{D}_j𝒜,`$ (55)
$`_t\stackrel{~}{\mathrm{\Gamma }}^i`$ $`=`$ $`_t^B\stackrel{~}{\mathrm{\Gamma }}^i+\kappa _{\stackrel{~}{\mathrm{\Gamma }}}\alpha \stackrel{~}{D}^i+\kappa _{\stackrel{~}{\mathrm{\Gamma }}𝒢1}\alpha 𝒢^i+\kappa _{\stackrel{~}{\mathrm{\Gamma }}𝒢2}\alpha \stackrel{~}{D}^j\stackrel{~}{D}_j𝒢^i+\kappa _{\stackrel{~}{\mathrm{\Gamma }}𝒢3}\alpha \stackrel{~}{D}^i\stackrel{~}{D}_j𝒢^j,`$ (56)
where $`\kappa `$s are possible multipliers (all $`\kappa =0`$ reduce the system the standard BSSN evolution equations).

We show the effects of each terms in Table 2. The CAFs in the table are on the flat space background. We see several terms produce negative real-part in CAFs, which might improve the stability than the previous system. (Table 2 includes again a column of diagonalizability of constraint propagation matrix $`M`$. Diagonalizable ones are expected to reflect the predictions from eigenvalue analysis. That is, the eigenvalue analysis with diagonalizable ones definitely avoid the diverging possibility in constraint propagation when it includes degenerated CAFs. See diagCP ). For the readers convenience, we list up several best candidates here.

(D) A system which has 7 negative CAFs
Simply adding $`\stackrel{~}{D}_{(i}_{j)}`$ term to $`_t\stackrel{~}{A}_{ij}`$ equation, say

$$_t\stackrel{~}{A}_{ij}=_t^{BSSN}\stackrel{~}{A}_{ij}+\kappa _{A2}\alpha (\stackrel{~}{D}_{(i}_{j)})$$
(57)
with $`\kappa _{A2}>0`$, CAFs on the flat background are 7 negative real CAFs.

(E) A system which has 6 negative CAFs
The below two adjustments will make 6 negative real CAFs, while they also produce one positive real CAF (a constraint violating mode). The effectiveness is not clear at this moment, but we think they are worth to be tested in numerical experiments.

(E1) With $`\kappa _{\stackrel{~}{\gamma }𝒢2}<0`$,

$$_t\stackrel{~}{\gamma }_{ij}=_t^{BSSN}\stackrel{~}{\gamma }_{ij}+\kappa _{\stackrel{~}{\gamma }𝒢2}\alpha \stackrel{~}{\gamma }_{k(i}\stackrel{~}{D}_{j)}𝒢^k.$$
(58)
(E2) With $`\kappa _{\stackrel{~}{\mathrm{\Gamma }}𝒢2}<0`$.

$$_t\stackrel{~}{\mathrm{\Gamma }}^i=_t^{BSSN}\stackrel{~}{\mathrm{\Gamma }}^i+\kappa _{\stackrel{~}{\mathrm{\Gamma }}𝒢2}\alpha \stackrel{~}{D}^j\stackrel{~}{D}_j𝒢^i.$$
(59)
## V Concluding remarks

Applying the constraint propagation analysis, we tried to understand why and how the so-called BSSN (Baumgarte-Shapiro-Shibata-Nakamura) re-formulation works better than the standard ADM equations in general relativistic numerical simulations. Our strategy was to evaluate eigenvalues of the constraint propagation equations in their Fourier modes, which method succeeded to explain the stability properties in many other systems in our series of works.

We have studied step-by-step where the replacements in the equations affect and/or newly added constraints work, by checking whether the error of constraints (if it exists) will decay or propagate away. Alcubierre et al potsdam9908  pointed out the importance of the replacement (adjustment) to the evolution equation using the momentum constraint, and our analysis clearly explains why they concluded this is the key. Not only this adjustment, we found, but also other adjustments and other introductions of new constraints also contribute to making the evolution system more stable. We found that if we missed a particular adjustment, then the expected stability behavior occationally gets worse than the ADM system. We further propose other adjustments of the set of equations which may have better features for numerical treatments.

The discussion in this article was only in the flat background space-time, and may not be applicable directly to the general numerical simulations. However, we rather believe that the general fundamental aspects of constraint propagation analysis are already revealed in this article. This is because, for the ADM and its adjusted cases, we found that the better formulations in the flat background are also better in the Schwarzschild space-time, while there are differences on the effective adjusting multipliers or the effective coordinate rangesadjADMsch ; adjADM .

We have not shown any numerical tests here. However, recently, the proposal (B) in §IV was examined numerically using linear wave initial data and confirmed to be effective for controlling the violation of the Hamiltonian constraint with our predicted multiplier signature masarushibata . The systematic numerical comparisons between different formulations are underway mexico , and we expect to have a chance to report them in near future. We are also trying to explain the stability of Laguna-Shoemaker’s implimented BSSN system PabloDeirdre  using the constraint propagation analysis.

There may not be the almighty formulation for any models in numerical relativity, but we believe our guidelines to find a better formulation in a systematic way will contribute a progress of this field. We hope the predictions in this paper will help the community to make further improvements.

## Acknowledgements

HS thanks T. Nakamura and M. Shibata for their comments, and he appreciates the hospitality of the Centre for Gravitational Wave Physics, the Pennsylvania State University, where the part of this work has been done. HS is supported by the special postdoctoral researchers’ program at RIKEN. This work was supported partially by the Grant-in-Aid for Scientific Research Fund of Japan Society of the Promotion of Science, No. 14740179.

## Appendix A Full set of BSSN constraint propagation equations

The constraint propagation equations of the BSSN system can be written as follows.

$`_t=\left((2/3)\alpha K+(2/3)\alpha 𝒜+\beta ^k_k\right)+\left(4e^{4\phi }\alpha (_k\phi )\stackrel{~}{\gamma }^{kj}2e^{4\phi }(_k\alpha )\stackrel{~}{\gamma }^{jk}\right)_j`$ (60)
$`+(2\alpha e^{4\phi }\stackrel{~}{A}^k{}_{j}{}^{}_{k}^{}\alpha e^{4\phi }(_j\stackrel{~}{A}_{kl})\stackrel{~}{\gamma }^{kl}e^{4\phi }(_j\alpha )𝒜e^{4\phi }\beta ^k_k_j`$
$`(1/2)e^{4\phi }\beta ^k\stackrel{~}{\gamma }^1(_j𝒮)_k+(1/6)e^{4\phi }\stackrel{~}{\gamma }^1(_j\beta ^k)(_k𝒮)(2/3)e^{4\phi }(_k\beta ^k)_j)𝒢^j`$
$`+(2\alpha e^{4\phi }\stackrel{~}{\gamma }^1\stackrel{~}{\gamma }^{lk}(_l\phi )𝒜_k+(1/2)\alpha e^{4\phi }\stackrel{~}{\gamma }^1(_l𝒜)\stackrel{~}{\gamma }^{lk}_k+(1/2)e^{4\phi }\stackrel{~}{\gamma }^1(_l\alpha )\stackrel{~}{\gamma }^{lk}𝒜_k`$
$`+(1/2)e^{4\phi }\stackrel{~}{\gamma }^1\beta ^m\stackrel{~}{\gamma }^{lk}_m_l_k(5/4)e^{4\phi }\stackrel{~}{\gamma }^2\beta ^m\stackrel{~}{\gamma }^{lk}(_m𝒮)_l_k+e^{4\phi }\stackrel{~}{\gamma }^1\beta ^m(_m\stackrel{~}{\gamma }^{lk})_l_k`$
$`+(1/2)e^{4\phi }\stackrel{~}{\gamma }^1\beta ^i(_j_i\stackrel{~}{\gamma }^{jk})_k+(3/4)e^{4\phi }\stackrel{~}{\gamma }^3\beta ^i\stackrel{~}{\gamma }^{jk}(_i𝒮)(_j𝒮)_k(3/4)e^{4\phi }\stackrel{~}{\gamma }^2\beta ^i(_i\stackrel{~}{\gamma }^{jk})(_j𝒮)_k`$
$`+(1/3)e^{4\phi }\stackrel{~}{\gamma }^1\stackrel{~}{\gamma }^{pj}(_j\beta ^k)_p_k`$
$`(5/12)e^{4\phi }\stackrel{~}{\gamma }^2\stackrel{~}{\gamma }^{jk}(_k\beta ^i)(_i𝒮)_j+(1/3)e^{4\phi }\stackrel{~}{\gamma }^1(_k\stackrel{~}{\gamma }^{ij})(_j\beta ^k)_i(1/6)e^{4\phi }\stackrel{~}{\gamma }^1\stackrel{~}{\gamma }^{mk}(_k_l\beta ^l)_m)𝒮`$
$`+((4/9)\alpha K𝒜(8/9)\alpha K^2+(4/3)\alpha e^{4\phi }(_i_j\phi )\stackrel{~}{\gamma }^{ij}+(8/3)\alpha e^{4\phi }(_k\phi )(_l\stackrel{~}{\gamma }^{lk})`$
$`+\alpha e^{4\phi }(_j\stackrel{~}{\gamma }^{jk})_k+8\alpha e^{4\phi }\stackrel{~}{\gamma }^{jk}(_j\phi )_k+\alpha e^{4\phi }\stackrel{~}{\gamma }^{jk}_j_k+8e^{4\phi }(_l\alpha )(_k\phi )\stackrel{~}{\gamma }^{lk}+e^{4\phi }(_l\alpha )(_k\stackrel{~}{\gamma }^{lk})`$
$`+2e^{4\phi }(_l\alpha )\stackrel{~}{\gamma }^{lk}_k+e^{4\phi }\stackrel{~}{\gamma }^{lk}(_l_k\alpha ))𝒜,`$
$`_t_i=((1/3)(_i\alpha )+(1/6)_i)+\alpha K_i+(\alpha e^{4\phi }\stackrel{~}{\gamma }^{km}(_k\phi )(_j\stackrel{~}{\gamma }_{mi})(1/2)\alpha e^{4\phi }\stackrel{~}{\mathrm{\Gamma }}_{kl}^m\stackrel{~}{\gamma }^{kl}(_j\stackrel{~}{\gamma }_{mi})`$ (61)
$`+(1/2)\alpha e^{4\phi }\stackrel{~}{\gamma }^{mk}(_k_j\stackrel{~}{\gamma }_{mi})+(1/2)\alpha e^{4\phi }\stackrel{~}{\gamma }^2(_i𝒮)(_j𝒮)(1/4)\alpha e^{4\phi }(_i\stackrel{~}{\gamma }_{kl})(_j\stackrel{~}{\gamma }^{kl})+\alpha e^{4\phi }\stackrel{~}{\gamma }^{km}(_k\phi )\stackrel{~}{\gamma }_{ji}_m`$
$`+\alpha e^{4\phi }(_j\phi )_i(1/2)\alpha e^{4\phi }\stackrel{~}{\mathrm{\Gamma }}_{kl}^m\stackrel{~}{\gamma }^{kl}\stackrel{~}{\gamma }_{ji}_m+\alpha e^{4\phi }\stackrel{~}{\gamma }^{mk}\stackrel{~}{\mathrm{\Gamma }}_{ijk}_m+(1/2)\alpha e^{4\phi }\stackrel{~}{\gamma }^{lk}\stackrel{~}{\gamma }_{ji}_k_l`$
$`+(1/2)e^{4\phi }\stackrel{~}{\gamma }^{mk}(_j\stackrel{~}{\gamma }_{im})(_k\alpha )+(1/2)e^{4\phi }(_j\alpha )_i+(1/2)e^{4\phi }\stackrel{~}{\gamma }^{mk}\stackrel{~}{\gamma }_{ji}(_k\alpha )_m)𝒢^j`$
$`+(\stackrel{~}{A}^k{}_{i}{}^{}(_k\alpha )+(1/9)(_i\alpha )K+(4/9)\alpha (_iK)+(1/9)\alpha K_i\alpha \stackrel{~}{A}^k{}_{i}{}^{}_{k}^{})𝒜,`$
$`_t𝒢^i=2\alpha \stackrel{~}{\gamma }^{ij}_j+((1/2)\beta ^k\stackrel{~}{\gamma }^{il}\stackrel{~}{\gamma }^2(_l𝒮)_k(1/2)\beta ^k\stackrel{~}{\gamma }^{in}(_k\stackrel{~}{\gamma }_{mn})\stackrel{~}{\gamma }^{ml}\stackrel{~}{\gamma }^1_l+(1/2)\beta ^k\stackrel{~}{\gamma }^{il}\stackrel{~}{\gamma }^1_l_k`$ (62)
$`(1/2)(_m\beta ^i)\stackrel{~}{\gamma }^{mk}\stackrel{~}{\gamma }^1_k+(1/3)(_l\beta ^l)\stackrel{~}{\gamma }^{ik}\stackrel{~}{\gamma }^1_k)𝒮+(+4\alpha \stackrel{~}{\gamma }^{ij}(\stackrel{~}{D}_j\phi )\alpha \stackrel{~}{\gamma }^{ij}_j(_k\alpha )\stackrel{~}{\gamma }^{ik})𝒜,`$
$`_t𝒮=+\beta ^k(_k𝒮)2\alpha \stackrel{~}{\gamma }𝒜,`$ (63)
$`_t𝒜=\left(\alpha K+\beta ^k_k\right)𝒜.`$ (64)
The flat background linear order equations, (44)-(48), were obtained from these expression.