text:
# On the derivation of the equations of motion in theories of gravity.

## 1 Introduction

General Relativity (GR) is unique among the class of field theories in the treatment of the equations of motion. The equations of motion of massive particles are completely determined by the field equation.
By comparison, classical electrodynamics postulates, in addition to the Maxwell field equations, the force (Coulomb or Lorentz) as well the form of the inertial term.
The derivation of the equation of motion for massive particles from the GR field equations was obtained in the pioneering work of Einstein, Infeld and Hoffman  of 1939. The particles were treated as singular points of the metric field. The first order approximation of the field equations resulted in the verification of Newton law of attraction between the singular points. In subsequent papers , the higher order approximations to the geodesic line equation (the equations of motion) was calculated. The successive approximations were all derived without additional hypothesis (i.e. the Geodesic Postulate). The method used in the cited papers is known as the EIH-procedure. The derivation is rather formal. A similar approach was adopted by Fock . He, however, studied Einstein field equations with the RHS-terms \- the tensor of energy-momentum of matter. Sternberg  proved the geodesic postulate in a general form for a system that incorporates also Yang-Mills fields. He refers to the establishment of the postulate as to the EIH condition. In his work there is no explicit computation of the trajectories for the $`N`$-body problem. The length and complexity of EIH analysis undoubtedly discouraged further research this area.
The derivability of the equations of motion is attributed to two factors:

* The GR field equations are non-linear (as opposed to the other classical field equations). As a result, the linear combination of 1-point static solutions does not provide a new static solution.
* The system of the GR field equations is subject to additional identities (Bianchi) which reduce the number of independent equations. The conservation laws (which rest upon these identities) play a crucial role in the original EIH-procedure.

Thus two natural questions arise:

* What other non-linear field equations have the property of derivability of the equations of motion?
* Is there a way to embed the moving singularities in a field satisfying the prescribed field equation.

In this paper we deal with the first question and touch upon the second one. First we define a suitable class of the field equations. These are linear in the second order derivatives and quadratic in the first order derivatives both with coefficients that depends only on the field variables. It is also required that the equations are Lorentz invariant. It will be demonstrated in the sequel that the linear part is the agent of inertia and the quadratic part is the agent of the force. Let the field equation have a static spherical symmetric solution (exact or approximate) which is singular at a point. Accordingly, by the Einstein description, this is the field of one particle. The field of a particle moving with constant velocity is taken to be the Lorentz transform of the static field. If it moves on a curved trajectory $`𝝍(t)`$ then the Lorentz transformation with the instantaneous velocity $`\dot{𝝍}`$ is taken. The field of $`N`$ particles is taken to be, approximately, the superposition of the fields of the single particles. When this is inserted to the field equation two singular expressions emerge: one from the linear part and the other one from the quadratic part. By equating the highest order terms of these singularities the strength of the singularity of the solution is reduced. It turns out that this balance equation is Newton’s law of attraction.

## 2 The general form of field equations and of their solutions.

We will require the following for the field equation:

* PDE of the second order,
* Linear in the second order derivatives with coefficients that depend only on the field variables,
* Quadratic in the first order derivatives.

The last condition can be motivated by the following dimensional argument.111It was communicated to us by F. W. Hehl. Let the field variable $`\mathrm{\Phi }`$ will be dimensionless. Thus, the second order derivatives $`\mathrm{\Phi }_{,\mu \nu }`$ have the dimension of $`[l^2]`$. The only polynomial of the first derivatives that has the same dimension is quadratic in the first order derivatives: $`\mathrm{\Phi }_{,\mu }\mathrm{\Phi }_{,\nu }`$. Thus, in order to have a field equation with dimensionless constants we have to consider a PDE, which is linear in the second order derivatives and quadratic in the first order derivatives.
As a consequence we will restrict ourselves to field equations of the form:

$$a\mathrm{\Phi }_{,\mu \nu }b\mathrm{\Phi }_{,\mu }\mathrm{\Phi }_{,\nu }=0,$$
(2.1)
where the coefficients $`a`$ and $`b`$ are the dimensionless functions of the field variable (or constants)

$$a=a(\mathrm{\Phi }),b=b(\mathrm{\Phi }).$$

The covariance condition as well as the existence of an action functional provide certain restrictions on the coefficient $`a(\mathrm{\Phi })`$ and $`b(\mathrm{\Phi })`$. Note that the form (2.1) is still symbolic, since the field $`\mathrm{\Phi }`$ as well as the coefficient functions $`a`$ and $`b`$ can possess the interior or coordinate indices and these indices should be contracted in a proper manner.
Note also that the Einstein equation does have the form (2.1).
Let us turn now to the properties of the possible physical relevant solutions of the field equation (2.1).
We will apply the Einstein model of a particle: it is a smooth solution of the field equation which tends to a constant at infinity and singular at a point. In this model the field pertaining to one particle has only one singularity. In the proper coordinates, the field around the singularity is spherically symmetric and tends to infinity at the singular point. The field pertaining to $`N`$ particles is approximated by the superposition of rigid motions of of the fields pertaining to all the particles. The trajectory of a singularity, which is taken to be the trajectory of the correspondent particle, is affected by the presence of the fields pertaining to the other particles. This means that the field equation should be non-linear (for a linear equation the fields can not affect each other).
Let us also require that a static spherical symmetric solution of (2.1) has a long distance expansion

$$\mathrm{\Phi }=\mathrm{\Phi }_0+\frac{A_1}{r}+\frac{A_2}{r^2}+\mathrm{}+\frac{A_n}{r^n}+\mathrm{},$$
(2.2)
where $`\{A_n,n=1,2,\mathrm{}\}`$ are constants.
It is enough to require the convergence of this expansion for distances greater of some characteristic length of the system. The dimension of $`A_n`$ is $`[l^n]`$. A second order equation usually produces two independent constants of integration. Thus the set of constants can be normalized:

$$A_1=m,A_n=C_nm^n,n=2,3,\mathrm{},$$
(2.3)
where the dimension of the constant $`m`$ is equal to $`[l]`$ and the constants $`C_i`$ are dimensionless.
Einstein derives the equation of motion by considering the integral relations which are implied by the conservation laws. We propose a different procedure to be exhibited in the sequel.

## 3 A Qualitative Description of the Algorithm

In this section we describe a novel algorithm for deriving the equations of motion from the field equation. We consider the general field equation (2.1) without specification of the tensorial nature of the field variable $`\mathrm{\Phi }`$. We do require, however, that (2.1) is Lorentz invariant.
Let the field equation (2.1) have a static spherical symmetric solution $`\mathrm{\Phi }(𝐫𝐫_0)`$ with a singularity located at $`𝐫=𝐫_0`$. Denote the Lorentz transformation based on the velocity $`v`$ by $`L_v`$. Consequently, if $`\mathrm{\Phi }=\mathrm{\Phi }(𝐫𝐫_0)`$ is a time independent solution of (2.1) then $`L_v\mathrm{\Phi }`$ is also a solution of the same equation for an arbitrary Lorentz transformation $`L_v`$. (If $`\mathrm{\Phi }`$ is multi-component like a tensor then $`L_v`$ involves not only a coordinate change but, also, a transformation of components of $`\mathrm{\Phi }`$.) The solution $`L_v\mathrm{\Phi }`$ describes the field of a pointwise singularity moving with a constant velocity $`v`$ on the trajectory $`𝝍=𝐯t`$. Let us try to construct a generalization of Lorentz transformation where the origin moves on a general trajectory $`𝝍=𝝍(t)`$. Denote such a transformation by $`N_\psi `$. The choice

$$N_\psi =L_{\dot{\psi }}$$

is a plausible candidate. $`N_\psi \mathrm{\Phi }`$ is a rigid motion of the field.
Now substitute $`N_\psi \mathrm{\Phi }`$ in (2.1). If $`\dot{\psi }=const`$ then $`N_\psi \mathrm{\Phi }`$ is also a solution of (2.1). If not, then the linear part produces extra terms. These extra terms come from two sources:

* The derivatives of the Lorentz root $`\sqrt{1\dot{\psi }^2}`$. In our consideration the velocity of the particle $`\dot{\psi }`$ as well as its time and spatial derivatives are assumed to be small. It follows that the derivatives of the root are of the form $`𝒪(\ddot{\psi }\dot{\psi })`$. Thus they may be discarded.
* The linear part. Since $`\mathrm{\Phi }`$ is time independent the first derivatives of $`\mathrm{\Phi }`$ are only the spatial ones. Thus, the first derivatives of of $`L_{\dot{\psi }}\mathrm{\Phi }(x)`$ involve spatial derivatives of $`\mathrm{\Phi }`$ multiplied by $`\dot{\psi }`$. The second derivatives of $`L_{\dot{\psi }}\mathrm{\Phi }(x)`$ cancel each other the same way as in a Lorentz transformed solution $`L_v\mathrm{\Phi }(x)`$ with one exception: The second derivative $`\ddot{\psi }`$ multiplied by spatial first derivatives of $`\mathrm{\Phi }`$ do remain. This extra term that come from the linear part is the agent of inertia.
* The quadratic part. It involves only first order derivatives. Consequently the fact that $`\dot{\psi }`$ is variable does not affect it’s form.

Construct now a solution that describes the field of $`N`$ particles, i.e. a field with $`N`$ singular points. It can be approximated by a superposition of 1-singular solutions moving on arbitrary trajectories $`{}_{}{}^{(j)}\psi ={}_{}{}^{(j)}\psi (t)`$

$$\mathrm{\Phi }=\underset{j}{}L_{{}_{}{}^{(j)}\dot{\psi }}{}_{}{}^{(j)}\mathrm{\Phi }(x).$$

Substituting this approximate solution in (2.1) we obtain, in the linear part, only the second derivatives $`{}_{}{}^{(j)}\ddot{\psi }`$ multiplied by the spatial first derivatives of $`\mathrm{\Phi }`$.
Consider the quadratic part. It is composed of the first order derivatives of $`{}_{}{}^{(j)}\mathrm{\Phi }(x)`$ multiplied by the first order derivatives of $`L_{{}_{}{}^{(k)}\mathrm{\Phi }}{}_{}{}^{(k)}\mathrm{\Phi }(x)`$. If $`j=k`$, since $`{}_{}{}^{(k)}\mathrm{\Phi }(x)`$ is a solution of (2.1), these products are cancelled by the linear part operating on $`{}_{}{}^{(k)}\mathrm{\Phi }(x)`$. If $`jk`$ the products will, hopefully, be an approximation to the interaction between the $`j`$-th and the $`k`$-th particles.
Near the $`k`$-th singularity, for the linear part, only the terms coming from $`{}_{}{}^{(k)}\mathrm{\Phi }(x)`$ will be dominant. Likewise, for the quadratic part, only the terms involving the derivatives of $`{}_{}{}^{(k)}\mathrm{\Phi }(x)`$ will be dominant. Equating, near the singularity, of the two terms above (again to the leading order) should, hopefully, result in the Newtonian law of attraction.
Let us summarize the novel algorithm for the derivation of equations of motion.

* Compute a static, spherically symmetric solution of the field equation. It will be singular at the origin. This will be taken to be the field generated by a single particle.
* Move the solution on a trajectory $`𝝍(t)`$ and apply the instantaneous Lorentz transformation based on $`\dot{𝝍}(t)`$.
* Take the field generated by $`n`$ particles to be the superposition of the fields generated by the single particles.
* Compute the leading part of the equation. Hopefully, only terms that involves $`\ddot{𝝍}`$ will be dominant.
* Compute the “force” between the particles by the quadratic part of the equation.
* Equate for each singularity, the highest order terms of the singularities that came from the linear part and the quadratic parts, respectively. This is an equation between the inertial part and the force.

## 4 Einsteinian gravity

Let us apply the procedure described above to the Einsteinian gravity. Our first step is to construct an approximated dynamical $`N`$-particle solution from the static 1-particle solution. The Einstein field equation in vacuo is

$$R_{ik}=0.$$
(4.1)
Consider the Schwarzschild solution in isotropic coordinates

$$ds^2=\left(\frac{1m/4r}{1+m/4r}\right)^2dt^2(1+m/4r)^4(dx^2+dy^2+dz^2).$$
(4.2)
We will use the general diagonal metric of the form

$$ds^2=e^{2f}dt^2e^{2g}(dx^2+dy^2+dz^2).$$
(4.3)
The static solution (4.2) has the long-distance expansion

$`f`$ $`=`$ $`{\displaystyle \frac{m}{2r}}+𝒪\left({\displaystyle \frac{m}{r}}\right)^3`$ (4.4)
$`g`$ $`=`$ $`{\displaystyle \frac{m}{2r}}{\displaystyle \frac{1}{16}}\left({\displaystyle \frac{m}{r}}\right)^2+𝒪\left({\displaystyle \frac{m}{r}}\right)^3.`$ (4.5)
Hence, up to $`𝒪\left(\frac{m}{r}\right)^3`$, the following relation holds

$$g=f\frac{1}{4}f^2.$$
(4.6)
We will seek the time depended solution of the field equation (4.1) in the diagonal form (4.3) with the functions $`f=f(x^i,t)`$ and $`g=g(x^i,t)`$.
Let us assume that the static relation (4.6) holds for the time dependent metric as well.
Thus, in order to determine the metric up to $`𝒪(\frac{m}{r})^3`$, it is enough to find the function $`f`$. This way, by simplifying the assumptions above, the solution of (4.1) is, approximately, dependent on one function $`f`$. If this is needed we can look at the scalar covariant equation

$$R=0.$$
(4.7)
Let us first insert the diagonal metric (4.3) in (4.7). The calculations result in

$$R=6e^{2f}(\ddot{g}\dot{f}\dot{g}+2\dot{g}^2)+2e^{2g}\left(\mathrm{}f+2\mathrm{}g+(f)^2+fg+(g)^2\right)$$
(4.8)
($`\dot{f}`$ is used for for the time derivative of $`f`$).
Insert the relation (4.6) to this equation. The result is

$$R=3e^{2f}\left((2f)\ddot{f}(7+5f+f^2)\dot{f}^2\right)+e^{2f+1/2f^2}\left(2(1+f)\mathrm{}f+(f+\frac{1}{2}f^2)(f)^2\right).$$
(4.9)
Estimate the terms of this equation. Note, first, that the the term containing $`\mathrm{}f`$ vanishes for the static as well as the time dependent solution obtained from it by Lorentz type transformations. As for the gradient term, it is estimated by

$$f(f)^2\frac{1}{r^2}𝒪(\frac{m}{r})^3.$$

(Note that the equation does not include the term $`\frac{1}{r^2}𝒪(\frac{m}{r})^2`$ which may come from the expression $`(f)^2`$.) Thus, up to $`\frac{1}{r^2}𝒪(\frac{m}{r})^3`$, the curvature scalar is

$$R=6\ddot{f}2\mathrm{}f.$$
(4.10)
Be rescaling the coordinates we obtain the Lorentz invariant equation

$$R=2\text{ }\text{ }\text{ }\text{ }f=0.$$
(4.11)
Recall that the truly dynamical variable are the metric components determined by the function $`\phi =e^{2f}`$. The equation for this function is

$$\text{ }\text{ }\text{ }\text{ }\phi =\frac{1}{\phi }(\phi )^2.$$
(4.12)
Note that the function $`\phi =1+\frac{m}{R}𝒪(1)`$. Hence up to $`\frac{1}{r^2}𝒪(\frac{m}{r})^3`$ the equation (4.12) takes the form

$$\text{ }\text{ }\text{ }\text{ }\phi =(\phi )^2.$$
(4.13)
This scalar equation is Lorentz invariant and non-linear. In the next section it will be shown that the algorithm is applicable for this equation.

## 5 A Non-linear scalar model

Let a flat Minkowskian $`4D`$-space with a metric $`\eta _{ab}=\eta ^{ab}=diag(1,1,1,1)`$ and a scalar field $`f=f(x)`$ be given. Consider, first, the wave equation - a linear second order equation for the field $`f`$

$$\text{ }\text{ }\text{ }\text{ }f=ddf=0.$$
(5.1)
It has a unique 1-singular at the origin spherical symmetric and asymptotically zero solution

$$f=\frac{m}{R}.$$
(5.2)
It is natural to interpret this solution as the field of a pointwise body located at the origin. The linearity of the field equation (5.1) yields the existence of $`N`$-singular asymptotically zero solution

$$f=\underset{i=1}{\overset{n}{}}\frac{m_i}{|𝑹_i|},𝑹_i=𝐫𝐫_i.$$
(5.3)
The solution (5.3) have to be interpreted as a field produced by a static configuration of bodies of masses $`m_i`$ located at the points $`𝐫_i`$. This interpretation is, of course, not physical because it completely ignores the interaction between the particles. One may attempt to apply now an instantaneous Lorentz transformation (with a time dependent velocity) of the static solution in order to produce a solution with interactions. For that, let the singular point move on an arbitrary trajectory $`𝝍(t)`$. Consider the function

$$f=\underset{i}{}\frac{m_i}{|𝑹_i|},$$
(5.4)
where

$$𝑹_i=(𝐫𝐫_i)\alpha _i\dot{𝝍}_i(\dot{𝝍}_i,(𝐫𝐫_i))\beta _i𝝍_i,$$
(5.5)
and the Lorentz parameters $`\alpha _i`$ and $`\beta _i`$ are functions only of $`|\dot{𝝍}_i|^2`$

$$\beta _i=\frac{1}{\sqrt{1|\dot{𝝍}_i|^2}},\alpha _i=\frac{1}{|\dot{𝝍}_i|^2}\left(1\frac{1}{\sqrt{1|\dot{𝝍}_i|^2}}\right)$$
(5.6)
This is the form of Lorentz transformation for an arbitrary velocity (Cf. Appendix A).
For $`𝝍_i=\dot{𝝍}_it`$, at the point $`𝑹_i=0`$

$$𝐫𝐫_i=\dot{𝝍}_i\left(\alpha _i(\dot{𝝍}_i,(𝐫𝐫_i))+\beta _it\right).$$
(5.7)
So

$$(\dot{𝝍}_i,(𝐫𝐫_i))=|\dot{𝝍}_i|^2\left(\alpha _i(\dot{𝝍}_i,(𝐫𝐫_i))+\beta _it\right).$$
(5.8)
Or

$$(\dot{𝝍}_i,(𝐫𝐫_i))=\frac{|\dot{𝝍}_i|^2\beta _it}{1\alpha _i|\dot{𝝍}_i|^2}.$$
(5.9)
Substituting this relation in (5.7) we obtain the equation motion of the singularity as

$$𝐫𝐫_i=\frac{\beta _i}{1\alpha _i|\dot{𝝍}_i|^2}\dot{𝝍}_it=\dot{𝝍}_it,$$
(5.10)
which is the motion of a free particle.
Calculate approximately the d’Alembertian of the function (5.4), in particular, omit the time derivatives of $`\alpha _i`$ and $`\beta _i`$ to obtain (Cf. Appendix B)

$$\text{ }\text{ }\text{ }\text{ }f=\underset{i=1}{\overset{n}{}}\frac{m_i\beta _i}{R_i^3}(\ddot{𝝍}_i,𝑹_i)+𝒪(\ddot{𝝍}\dot{𝝍}).$$
(5.11)
Thus, for the field equation (5.1) $`f`$ is a solution if and only if $`\ddot{𝝍}_i=0`$.
Consequently, the linear equation describes a free inertial motion of an arbitrary system of singularities.
Consider now a nonlinear covariant field equation

$$ddf=kdfdf,$$
(5.12)
or, in coordinate form

$$\text{ }\text{ }\text{ }\text{ }f=k\eta ^{ab}f_{,a}f_b.$$
(5.13)
It is easy to see that the the field equation (5.13) can be transformed to the field equation (5.1) by the redefinition of the scalar field

$$\phi =e^{kf}.$$
(5.14)
If the field $`f`$ satisfies the equation (5.13) the new field $`\phi `$ satisfies (5.1). The transformation (5.14) can be used in order to derive a close form solution to the nonlinear equation (5.13).

$$f=\frac{1}{k}\mathrm{ln}\left(1+k\frac{m}{|𝑹|}\right).$$
(5.15)
In order to have a singularity only at the origin we have to require the constants $`k`$ and $`m`$ to be positive.
The solution (5.15) is singular at the origin and vanishes at the infinity. It’s long-distance expansion begins with $`m/r`$. In the limit case $`k0`$ the solution (5.15) approaches the Newtonian potential $`f\frac{m}{r}`$.
Because of the transformation (5.14) the equation (5.13) has also a static $`N`$-singular points solution

$$f=\frac{1}{k}\mathrm{ln}\left(1+k\underset{i=1}{\overset{N}{}}\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}\right),$$
(5.16)
where $`{}_{}{}^{(i)}𝑹=𝐫𝐫_i`$ and $`𝐫_i`$ are the radius vectors of the $`N`$ fixed points.
The solution (5.16) is static. Let us look for an approximate dynamic one. In order to construct an approximate solution take the superposition of the fields

$$f=\frac{1}{k}\underset{i=1}{\overset{N}{}}\mathrm{ln}\left(1+k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}\right).$$
(5.17)
Let us apply the Lorentz transformation with variable velocities. Now the vector $`𝑹`$ c.f. (5.5) is time dependent. Calculate the linear part of the equation (Cf. Appendix C) to obtain

$$\text{ }\text{ }\text{ }\text{ }f=\underset{i=1}{\overset{n}{}}\left(\frac{{}_{}{}^{(i)}m{}_{}{}^{(i)}\beta }{{}_{}{}^{(i)}R_{}^{3}}\frac{({}_{}{}^{(i)}𝑹,{}_{}{}^{(i)}\ddot{𝝍})}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}\frac{{}_{}{}^{(i)}m_{}^{2}k}{{}_{}{}^{(i)}R_{}^{4}}\frac{1}{(1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R})^2}\right).$$
(5.18)
By Appendix D the nonlinear part is:

$`k\eta ^{ab}f_{,a}f_{,b}`$ $`=`$ $`k{\displaystyle \underset{i}{}}{\displaystyle \frac{{}_{}{}^{(i)}m_{}^{2}k}{{}_{}{}^{(i)}R_{}^{4}}}{\displaystyle \frac{1}{(1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R})^2}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}}{\displaystyle \frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}}[({}_{}{}^{(i)}𝑹,{}_{}{}^{(j)}𝑹)+`$
$`({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}𝑹)({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\beta {}_{}{}^{(j)}\beta +{}_{}{}^{(i)}\alpha +{}_{}{}^{(j)}\alpha {}_{}{}^{(i)}\alpha {}_{}{}^{(j)}\alpha ({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}\dot{𝝍}))].`$
The field equation takes the form

$`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{{}_{}{}^{(i)}m{}_{}{}^{(i)}\beta }{{}_{}{}^{(i)}R_{}^{3}}}{\displaystyle \frac{({}_{}{}^{(i)}𝑹,{}_{}{}^{(i)}\ddot{𝝍})}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}}=k{\displaystyle \underset{ij}{}}{\displaystyle \frac{\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}}{\displaystyle \frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}}[({}_{}{}^{(i)}𝑹,{}_{}{}^{(j)}𝑹)+`$
$`({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}𝑹)({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\beta {}_{}{}^{(j)}\beta +{}_{}{}^{(i)}\alpha +{}_{}{}^{(j)}\alpha {}_{}{}^{(i)}\alpha {}_{}{}^{(j)}\alpha ({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}\dot{𝝍}))]+𝒪(\ddot{\psi }\dot{\psi }).`$ (5.19)
For the approximation of slow motions it is

$$\underset{i=1}{\overset{n}{}}\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}\frac{({}_{}{}^{(i)}𝑹,{}_{}{}^{(i)}\ddot{𝝍})}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}=k\underset{ij}{}\frac{\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}\frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}({}_{}{}^{(i)}𝑹,{}_{}{}^{(j)}𝑹)+𝒪(\ddot{\psi }\dot{\psi })+\frac{1}{R^2}𝒪(\frac{m}{R})^2.$$
(5.20)
The two sides of this equation are functions of the variable $`x`$. Consider the $`p`$-th singularity. Take an arbitrary point $`x`$ close to this singularity. It follows that

$${}_{}{}^{(p)}𝑹0\text{ and }{}_{}{}^{(i)}𝑹\overline{R}_{ip}\text{ for }ip,$$
(5.21)
where $`\overline{R}_{ip}`$ is a vector from the point $`i`$ to the point $`p`$.
In the LHS of the equation (5.20) there is one singular term

$$\frac{{}_{}{}^{(p)}m}{{}_{}{}^{(p)}R_{}^{3}}\frac{({}_{}{}^{(p)}𝑹,{}_{}{}^{(p)}\ddot{𝝍})}{1k\frac{{}_{}{}^{(p)}m}{{}_{}{}^{(p)}R}}.$$
(5.22)
The singular term in the RHS of (5.20) is

$$k\frac{\frac{{}_{}{}^{(p)}m}{{}_{}{}^{(p)}R_{}^{3}}}{1k\frac{{}_{}{}^{(p)}m}{{}_{}{}^{(p)}R}}\underset{jp}{}\frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}({}_{}{}^{(p)}𝑹,{}_{}{}^{(j)}𝑹).$$
(5.23)
The terms (5.22) and (5.23) are $`𝒪(R^2)`$ near the singularity. When these are inserted to the RHS and LHS of (5.20), respectively, the remainder will be $`𝒪(R^1)`$ if only if:

$$({}_{}{}^{(p)}𝑹,{}_{}{}^{(p)}\ddot{𝝍})=k\underset{jp}{}\frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}({}_{}{}^{(p)}𝑹,{}_{}{}^{(j)}𝑹)$$
(5.24)
This is the only way to diminish the strength of the singularity at (5.20).
Take into account that the point $`x`$ is still arbitrary. Hence (5.24) is valid only if

$${}_{}{}^{(p)}\ddot{𝝍}=k\underset{jp}{}\frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}{}_{}{}^{(j)}𝑹$$
(5.25)
For the limiting values in (5.21)

$${}_{}{}^{(p)}\ddot{𝝍}=k\underset{jp}{}\frac{{}_{}{}^{(j)}m𝑹_{jp}}{R_{jp}^3}\frac{1}{1k\frac{{}_{}{}^{(j)}m}{R_{jp}}}$$
(5.26)
The second fraction differs from 1, significantly, only for small distances comparable to the Schwarzschild radius $`r=km`$. Thus, it can be neglected.
It follows that

$${}_{}{}^{(p)}\ddot{𝝍}=k\underset{jp}{}\frac{{}_{}{}^{(j)}m𝑹_{jp}}{R_{jp}^3}$$
(5.27)
For a system of two singular points

$${}_{}{}^{(1)}\ddot{𝝍}=k\frac{{}_{}{}^{(2)}m𝑹_{21}}{R_{21}^3}$$
(5.28)
For $`k<0`$ (5.27) and (5.28) result in attraction between the particles. The absolute value of $`k`$ is unimportant, since it amounts to the rescaling of the mass.
This way Newton’s law is obtained.

## 6 The Embedding problem

The trajectories obtained are approximate. At this point, two avenues are open. The first one, which is adopted by EIH is to get higher order approximations to the trajectories. This procedure is also used in the PPN approach. By these methods, the successive approximations become highly singular near the particle trajectories.
The second avenue is to embed the singularities in a field satisfying the field equations. For that purpose, the successive approximations should add regular terms (and, possible, low order singular terms) near the trajectories. In this article we just touch upon the method to be developed. Consider a first order solution i.e. the trajectories are taken to be fixed and the calculations is done only for the first order. For the scalar model the desired field will be a solution $`\phi `$ of (5.13),

$$\phi =e^{f+g},$$
(6.1)
where $`f`$ is defined by (5.17) and $`g`$ is regular at the vicinity of all the singularities (“multi Green function”). Let us present a formal argument that shows that the construction is possible. Take, approximately $`f=𝒪(\frac{m}{r})`$It follows from the computations in the previous section that

$$(\text{ }\text{ }\text{ }\text{ }f\eta ^{ab}f_{,a}f_{,b})=mF$$
(6.2)
(5.13) is satisfied when

$$\text{ }\text{ }\text{ }\text{ }g2\eta ^{ab}f_{,a}g_{,b}+\eta ^{ab}g_{,a}g_{,b}=mF.$$
(6.3)
Newton’s law holds, thus $`F=𝒪(\frac{1}{r})`$. Consider for $`F`$ two examples:

$$F=\frac{m}{r}+\text{regular terms}$$

or

$$F=\frac{m<r,a>^2}{r^3}+\text{regular terms}$$

In the first case, take $`g^{(1)}=\frac{1}{2}cr^2`$, in the second case ( a is a constant vector) take $`g^{(1)}=\frac{1}{2}c<r,a>r`$. In both cases the remainder in the left hand side of (5.13) will be regular.

## Appendix A Lorentz transformations

The Lorentz transformations are usually written in a very special case when the axes of two reference systems are parallel and one of the systems moves relative to the other with a velocity parallel to the $`x`$-axis 222Recall that we use the system of units with $`c=1`$

$$x=\alpha (\stackrel{~}{x}+v\stackrel{~}{t}),y=\stackrel{~}{y},z=\stackrel{~}{z},t=\alpha (\stackrel{~}{t}+v\stackrel{~}{x}),$$
(A.1)
with the Lorentz parameter

$$\alpha =(1v^2)^{\frac{1}{2}}.$$
(A.2)
It is well known that the Lorentz transformations are non-commutative. Consequently, a general transformation with an arbitrary directed vector of velocity can not be generated as a successive application of three orthogonal transformations (relative to the axes).
In the sequel, a formula for Lorentz transformation for general velocity vector is exhibited. It is hard to believe that such formula does not exists in the literature. The authors, however, could not find a reference for it.
Consider a reference system moving with an arbitrary directed velocity $`𝐯`$ with the axes parallel to the corresponding axes of a rest reference system. Consider a vector $`𝐫`$ to an arbitrary point in space. The projection of the vector $`𝐫`$ on the direction $`𝐯`$ will be

$$P_𝐯𝐫=\frac{𝐯(𝐯,𝐫)}{v^2}.$$
(A.3)
Expand the vector $`𝐫`$ as

$$𝐫=P_𝐯𝐫+N_𝐯𝐫.$$
(A.4)
For a Lorentz transformation

$`P_𝐯\stackrel{~}{𝐫}`$ $`=`$ $`\alpha (P_𝐯𝐫+𝐯t),`$ (A.5)
$`N_𝐯\stackrel{~}{𝐫}`$ $`=`$ $`N_𝐯𝐫.`$ (A.6)
Consequently, the transform of the spatial coordinates is

$$\stackrel{~}{𝐫}=\alpha (P_𝐯𝐫+𝐯t)+N_𝐯$$
(A.7)
or, explicitly

$$\stackrel{~}{𝐫}=𝐫+𝐯\left(\alpha t\frac{(1\alpha )}{v^2}(𝐫,𝐯)\right).$$
(A.8)
The change of the time coordinate is also governed only by the tangential part of the vector $`𝐫`$

$$\stackrel{~}{t}=\alpha (t+|P_𝐯𝐫|v)=\alpha \left(t+(P_𝐯𝐫,𝐯)\right)$$
(A.9)
or explicitly

$$\stackrel{~}{t}=\alpha \left(t+(𝐫,𝐯)\right).$$
(A.10)
In the special case of a motion parallel to the axis $`x`$ the relations (A.8), (A.10) reduce to the ordinary form of Lorentz transformation (A.1).
Therefore an arbitrary Lorentz transformation takes the form

$$\{\begin{array}{cc}& \stackrel{~}{t}=\alpha \left(t+(𝐫,𝐯)\right)\hfill \\ & \stackrel{~}{𝐫}=𝐫+𝐯\left(\alpha t\frac{(1\alpha )}{v^2}(𝐫,𝐯)\right)\hfill \end{array}$$
(A.11)
Taking the derivative of (A.11) the law of transformations of the differentials is obtained

$$\{\begin{array}{cc}& d\stackrel{~}{t}=\alpha \left(dtv_jdx^j\right)\hfill \\ & d\stackrel{~}{x}^i=dx^i+v^i\left(\alpha dt+\frac{(1\alpha )}{v^2}v_jdx^j\right)\hfill \end{array}$$
(A.12)
This can be represented in a matrix form by

$$d\stackrel{~}{x}=Adx<===>d\stackrel{~}{x}^\mu =A_{}^{\mu }{}_{\nu }{}^{}dx^\nu $$
(A.13)
with

$$A=\left(\begin{array}{cccc}\alpha & \alpha v^1& \alpha v^2& \alpha v^3\\ \alpha v^1& 1+(1\alpha )\frac{v_1^2}{v^2}& (1\alpha )\frac{v_1v_2}{v^2}& (1\alpha )\frac{v_1v_3}{v^2}\\ \alpha v^2& (1\alpha )\frac{v_1v_2}{v^2}& 1+(1\alpha )\frac{v_2^2}{v^2}& (1\alpha )\frac{v_1v_3}{v^2}\\ \alpha v^3& (1\alpha )\frac{v_1v_3}{v^2}& (1\alpha )\frac{v_2v_3}{v^2}& 1+(1\alpha )\frac{v_3^2}{v^2}\end{array}\right)$$
(A.14)
For our purposes it is enough to have the second order $`(\frac{1}{c^2})`$ approximation of the transforms (A.8), (A.10). Use

$$\alpha =\frac{1}{\sqrt{1v^2}}=1+\frac{1}{2}v^2+O(v^4)$$
(A.15)
to obtain

$$\{\begin{array}{cc}& \stackrel{~}{𝐫}=𝐫+𝐯t+𝐯(𝐫,𝐯)\hfill \\ & \stackrel{~}{t}=(1\frac{1}{2}v^2)t+(𝐫,𝐯).\hfill \end{array}$$
(A.16)
Consequently the approximation of the transformation law for the differentials is

$$\{\begin{array}{cc}& d\stackrel{~}{t}=(1+\frac{1}{2}v^2)dtv_jdx^j\hfill \\ & d\stackrel{~}{x}^i=dx^i+v^i\left(dt\frac{1}{2}v_jdx^j\right).\hfill \end{array}$$
(A.17)
Thus the matrix of transformation is

$$A=\left(\begin{array}{cccc}1+\frac{1}{2}v^2& v^1& v^2& v^3\\ v^1& 1\frac{1}{2}v_1^2& \frac{1}{2}v_1v_2& \frac{1}{2}v_1v_3\\ v^2& \frac{1}{2}v_1v_2& 1\frac{1}{2}v_2^2& \frac{1}{2}v_1v_3\\ v^3& \frac{1}{2}v_1v_3& \frac{1}{2}v_2v_3& 1\frac{1}{2}v_3^2.\end{array}\right).$$
(A.18)
## Appendix B 1-point ansatz

Calculate the d’Alembertian of the field $`f`$ of a moving particle $`f`$:

$$f=\frac{m}{R},$$
(B.1)
where

$$𝑹=(𝐫𝐫_0)\alpha \dot{\overline{\psi }}(\dot{𝝍},(𝐫𝐫_0))\beta \overline{\psi },$$
(B.2)
with $`\overline{\psi }=\overline{\psi }(t)`$ and $`\alpha `$ and $`\beta `$ are functions of $`|\dot{𝝍}|^2`$:

$$\beta =\frac{1}{\sqrt{1|\dot{𝝍}|^2}},\alpha =\frac{1}{|\dot{𝝍}|^2}\left(1\frac{1}{\sqrt{1|\dot{𝝍}|^2}}\right)$$
(B.3)
By straightforward calculation

$`𝑹_t`$ $`=`$ $`2\alpha ^{}\left((\dot{𝝍},\ddot{𝝍})\dot{𝝍}(\dot{𝝍},𝐫𝐫_0)\right)\alpha \left(\ddot{𝝍}(\dot{𝝍},𝐫𝐫_0)\right)`$ (B.4)
$`\alpha \left(\dot{𝝍}(\ddot{𝝍},𝐫𝐫_0)\right)2\beta ^{}\left((\dot{𝝍},\ddot{𝝍})\overline{\psi }\right)\beta \dot{𝝍}`$
and

$`𝑹_{tt}`$ $`=`$ $`4\alpha ^{\prime \prime }\left((\dot{𝝍},\ddot{𝝍})^2(\dot{𝝍},(𝐫𝐫_0))\right)2\alpha ^{}\left(|\ddot{𝝍}|^2(\dot{𝝍},(𝐫𝐫_0))\right)`$ (B.5)
$`2\alpha ^{}\left((\dot{𝝍},\dot{\ddot{𝝍}})(\dot{𝝍},(𝐫𝐫_0))\right)2\alpha ^{}\left((\dot{𝝍},\ddot{𝝍})(\ddot{𝝍},(𝐫𝐫_0))\right)`$
$`4\beta ^{\prime \prime }\left((\dot{𝝍},\ddot{𝝍})^2𝝍\right)2\beta ^{}\left(|\ddot{𝝍}|^2𝝍\right)2\beta ^{}\left((\dot{𝝍},\dot{\ddot{𝝍}})𝝍\right)2\beta ^{}\left((\dot{𝝍},\ddot{𝝍})\dot{𝝍}\right)`$
$`2\beta ^{}\left((\dot{𝝍},\ddot{𝝍})\dot{𝝍}\right)\beta \ddot{𝝍}.`$
Thus

$`f_t`$ $`=`$ $`{\displaystyle \frac{m}{R^3}}(𝑹_t,𝑹)={\displaystyle \frac{m}{R^3}}(2\alpha ^{}\left((\dot{𝝍},\ddot{𝝍})(\dot{𝝍},𝑹)(\dot{𝝍},𝐫𝐫_0)\right)`$ (B.6)
$`\alpha \left((\ddot{𝝍},𝑹)(\dot{𝝍},𝐫𝐫_0)\right)\alpha \left((\dot{𝝍},R)(\ddot{𝝍},𝐫𝐫_0)\right)`$
$`2\beta ^{}\left((\dot{𝝍},\ddot{𝝍})(\overline{\psi },𝑹)\right)\beta \left((\dot{𝝍},R)\right))`$
In the first approximation (for velocities that are small with respect to the speed of the light) we can take

$$𝑹_t=\beta \dot{𝝍},𝑹_{tt}=\beta \ddot{𝝍}.$$
(B.7)
Consequently,

$$f_t=\frac{m}{R^3}\beta (\dot{𝝍},R)$$
(B.8)
In general, the derivatives of $`\alpha `$ and $`\beta `$ contribute terms that are quadratic in $`\dot{𝝍}`$ and its derivatives. Hence, they can be neglected. Thus, the second time derivative, to the same accuracy, is

$`f_{tt}`$ $`=`$ $`{\displaystyle \frac{m\beta }{R^3}}\left((\ddot{𝝍},𝑹)+(\dot{𝝍},𝑹_t)\right)3{\displaystyle \frac{m\beta }{R^5}}(\dot{𝝍},𝑹)(𝑹_t,𝑹)`$ (B.9)
Substitute

$$𝑹_t=\beta \dot{𝝍}$$
(B.10)
to get

$`f_{tt}`$ $`=`$ $`{\displaystyle \frac{m\beta }{R^3}}\left((\ddot{𝝍},𝑹)\beta (\dot{𝝍},\dot{𝝍})\right)+3{\displaystyle \frac{m\beta ^2}{R^5}}(\dot{𝝍},𝑹)^2`$ (B.11)
$`=`$ $`{\displaystyle \frac{m\beta }{R^3}}(\ddot{𝝍},𝑹)+{\displaystyle \frac{m\beta ^2}{R^5}}\left(3(\dot{𝝍},𝑹)^2R^2|\dot{𝝍}|^2\right).`$
As for the spatial derivatives

$$𝑹_x=𝒆_1\alpha \dot{\psi }_1\dot{𝝍},𝑹_{xx}=0$$
(B.12)
$$f_x=\frac{m}{R^3}(𝑹_x,𝑹)$$
(B.13)
$$f_{xx}=\frac{m}{R^3}(𝑹_x,𝑹_x)+3\frac{m}{R^5}(𝑹_x,𝑹)^2$$
(B.14)
$$\mathrm{}f=3\frac{m}{R^5}\left((𝑹_x,𝑹)^2+(𝑹_y,𝑹)^2+(𝑹_z,𝑹)^2\right)\frac{m}{R^3}\left(R_x^2+R_y^2+R_z^2\right)$$
(B.15)
Since

$$(𝑹_x,𝑹)=R_1\alpha \dot{\psi }_1(\dot{𝝍},R)$$
(B.16)
we get

$$(𝑹_x,𝑹)^2+(𝑹_y,𝑹)^2+(𝑹_z,𝑹)^2=R^22\alpha (\dot{𝝍},R)^2+\alpha ^2|\dot{𝝍}|^2(\dot{𝝍},R)^2$$
(B.17)
and

$$R_x^2+R_y^2+R_z^2=32\alpha |\dot{𝝍}|^2+\alpha ^2|\dot{𝝍}|^4$$
(B.18)
Thus the Laplacian is

$$\mathrm{}f=\frac{m}{R^5}(\alpha ^2|\dot{𝝍}|^22\alpha )\left(3(\dot{𝝍},R)^2|\dot{𝝍}|^2R^2\right)$$
(B.19)
And, for the d’Alembertian

$$\text{ }\text{ }\text{ }\text{ }f=\frac{m\beta }{R^3}(\ddot{𝝍},𝑹)+\frac{m}{R^5}(\beta ^2+2\alpha \alpha ^2|\dot{𝝍}|^2)\left(3(\dot{𝝍},𝑹)^2R^2|\dot{𝝍}|^2\right)$$
(B.20)
Using the expressions for the functions $`\alpha ,\beta `$ we obtain, approximately,

$$\text{ }\text{ }\text{ }\text{ }f=\frac{m\beta }{R^3}(\ddot{𝝍},𝑹)$$
(B.21)
## Appendix C The leading part

From the calculations above

$$𝑹_t=\beta \dot{𝝍},𝑹_{tt}=\beta \ddot{𝝍}.$$
(C.1)
Thus

$$f_t=\frac{m}{R^3}\frac{(𝑹,𝑹_t)}{1k\frac{m}{R}}=\frac{m\beta }{R^3}\frac{(𝑹,\dot{𝝍})}{1k\frac{m}{R}}.$$
(C.2)
And

$`f_{tt}={\displaystyle \frac{m\beta }{R^5}}{\displaystyle \frac{3\beta (𝑹,\dot{𝝍})^2+R^2(𝑹,\ddot{𝝍})\beta R^2|\dot{𝝍}|^2}{1k\frac{m}{R}}}+{\displaystyle \frac{m^2\beta ^2k}{R^6}}{\displaystyle \frac{(𝑹,\dot{𝝍})^2}{(1k\frac{m}{R})^2}}.`$ (C.3)
As for the spatial derivatives we have

$$𝑹_x=𝒆_1\alpha \dot{𝝍}(\dot{𝝍},𝒆_1),𝑹_{xx}=0,$$
(C.4)
where $`e_1`$ is a unit vector on the $`x`$ axis.

$$f_x=\frac{m}{R^3}\frac{(𝑹,𝑹_x)}{1k\frac{m}{R}}$$
(C.5)
$$f_{xx}=3\frac{m}{R^5}\frac{(𝑹,𝑹_x)^2}{1k\frac{m}{R}}\frac{m}{R^3}\frac{(𝑹_x,𝑹_x)}{1k\frac{m}{R}}+\frac{km^2}{R^6}\frac{(𝑹,𝑹_x)^2}{(1k\frac{m}{R})^2}.$$
(C.6)
Thus

$`\mathrm{}f`$ $`=`$ $`3{\displaystyle \frac{m}{R^5}}{\displaystyle \frac{(𝑹,𝑹_x)^2+(𝑹,𝑹_y)^2+(𝑹,𝑹_z)^2}{1k\frac{m}{R}}}`$ (C.7)
$`{\displaystyle \frac{m}{R^3}}{\displaystyle \frac{(𝑹_x,𝑹_x)+(𝑹_y,𝑹_y)+(𝑹_z,𝑹_z)}{1k\frac{m}{R}}}+`$
$`{\displaystyle \frac{km^2}{R^6}}{\displaystyle \frac{(𝑹,𝑹_x)^2+(𝑹,𝑹_y)^2+(𝑹,𝑹_z)^2}{(1k\frac{m}{R})^2}}`$
Substitute the value of $`𝑹_x`$ to get

$`\mathrm{}f`$ $`=`$ $`3{\displaystyle \frac{m}{R^5}}{\displaystyle \frac{R^22\alpha (\dot{𝝍},𝑹)^2+\alpha ^2|\dot{𝝍}|^2(\dot{𝝍},𝑹)^2}{1k\frac{m}{R}}}{\displaystyle \frac{m}{R^5}}{\displaystyle \frac{32\alpha |\dot{𝝍}|^2+\alpha ^2|\dot{𝝍}|^4}{1k\frac{m}{R}}}+`$ (C.8)
$`{\displaystyle \frac{km^2}{R^6}}{\displaystyle \frac{R^22\alpha (\dot{𝝍},𝑹)^2+\alpha ^2|\dot{𝝍}|^2(\dot{𝝍},𝑹)^2}{(1k\frac{m}{R})^2}}`$
Thus the second order l.h.s. of the equation is

$`\text{ }\text{ }\text{ }f`$ $`=`$ $`{\displaystyle \frac{m\beta }{R^5}}{\displaystyle \frac{3\beta (𝑹,\dot{𝝍})^2+R^2(𝑹,\ddot{𝝍})\beta R^2|\dot{𝝍}|^2}{1k\frac{m}{R}}}+{\displaystyle \frac{m^2\beta ^2k}{R^6}}{\displaystyle \frac{(𝑹,\dot{𝝍})^2}{(1k\frac{m}{R})^2}}`$ (C.9)
$`3{\displaystyle \frac{m}{R^5}}{\displaystyle \frac{R^22\alpha (\dot{𝝍},𝑹)^2+\alpha ^2|\dot{𝝍}|^2(\dot{𝝍},𝑹)^2}{1k\frac{m}{R}}}+{\displaystyle \frac{m}{R^3}}{\displaystyle \frac{32\alpha |\dot{𝝍}|^2+\alpha ^2|\dot{𝝍}|^4}{1k\frac{m}{R}}}`$
$`{\displaystyle \frac{km^2}{R^6}}{\displaystyle \frac{R^22\alpha (\dot{𝝍},𝑹)^2+\alpha ^2|\dot{𝝍}|^2(\dot{𝝍},𝑹)^2}{(1k\frac{m}{R})^2}}`$
$`=`$ $`{\displaystyle \frac{m\beta }{R^3}}{\displaystyle \frac{(𝑹,\ddot{𝝍})}{1k\frac{m}{R}}}+{\displaystyle \frac{m}{R^5}}(\beta ^2+2\alpha \alpha ^2|\dot{𝝍}|^2){\displaystyle \frac{(3(𝑹,\dot{𝝍})^2R^2|\dot{𝝍}|^2)}{1k\frac{m}{R}}}+`$
$`{\displaystyle \frac{m^2k}{R^6}}{\displaystyle \frac{R^2+(\beta ^2+2\alpha \alpha ^2|\dot{𝝍}|^2)(𝑹,\dot{𝝍})^2}{(1k\frac{m}{R})^2}}`$
Using the relation $`\beta ^2=\alpha ^2|\dot{𝝍}|^22\alpha `$ we obtain

$$\text{ }\text{ }\text{ }\text{ }f=\frac{m\beta }{R^3}\frac{(𝑹,\ddot{𝝍})}{1k\frac{m}{R}}\frac{m^2k}{R^4}\frac{1}{(1k\frac{m}{R})^2}$$
(C.10)
As for the quadratic r.h.s.

$`\eta ^{ab}f_{,a}f_{,b}`$ $`=`$ $`{\displaystyle \frac{m^2\beta ^2}{R^6}}{\displaystyle \frac{(𝑹,\dot{𝝍})^2}{(1k\frac{m}{R})^2}}{\displaystyle \frac{m^2}{R^6}}{\displaystyle \frac{(𝑹,𝑹_x)^2+(𝑹,𝑹_y)^2+(𝑹,𝑹_z)^2}{(1k\frac{m}{R})^2}}`$ (C.11)
$`=`$ $`{\displaystyle \frac{m^2\beta ^2}{R^6}}{\displaystyle \frac{(𝑹,\dot{𝝍})^2}{(1k\frac{m}{R})^2}}{\displaystyle \frac{m^2}{R^6}}{\displaystyle \frac{R^22\alpha (\dot{𝝍},𝑹)^2+\alpha ^2|\dot{𝝍}|^2(\dot{𝝍},𝑹)^2}{(1k\frac{m}{R})^2}}`$
$`=`$ $`{\displaystyle \frac{m^2}{R^4}}{\displaystyle \frac{1}{(1k\frac{m}{R})^2}}+{\displaystyle \frac{m^2}{R^6}}(\beta ^2+2\alpha \alpha ^2|\dot{𝝍}|^2){\displaystyle \frac{(𝑹,\dot{𝝍})^2}{(1k\frac{m}{R})^2}}`$
Thus, to the first order,

$$\eta ^{ab}f_{,a}f_{,b}=\frac{m^2}{R^4}\frac{1}{(1k\frac{m}{R})^2}$$
(C.12)
## Appendix D The quadratic part

Calculate the nonlinear part

$$f_t=\underset{i=1}{\overset{n}{}}\left(\frac{({}_{}{}^{(i)}𝑹_{t}^{},{}_{}{}^{(i)}𝑹)}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}\right)=\underset{i=1}{\overset{n}{}}\left(\frac{({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}\frac{{}_{}{}^{(i)}m{}_{}{}^{(i)}\beta }{{}_{}{}^{(i)}R_{}^{3}}\right).$$
(D.1)
Thus

$$(f_t)^2=\underset{i,j=1}{\overset{n}{}}\frac{\frac{{}_{}{}^{(i)}m{}_{}{}^{(i)}\beta }{{}_{}{}^{(i)}R_{}^{3}}}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}\frac{\frac{{}_{}{}^{(j)}m{}_{}{}^{(j)}\beta }{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}\dot{𝝍},{}_{}{}^{(j)}𝑹).$$
(D.2)
As for the spatial derivatives

$$f_x=\underset{i=1}{\overset{n}{}}\frac{({}_{}{}^{(i)}𝑹_{x}^{},{}_{}{}^{(i)}𝑹)}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}$$
(D.3)
$`(f,f)`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{n}{}}}{\displaystyle \frac{\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}}{\displaystyle \frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}}(({}_{}{}^{(i)}𝑹_{x}^{},{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}𝑹_{x}^{},{}_{}{}^{(j)}𝑹)+`$ (D.4)
$`({}_{}{}^{(i)}𝑹_{y}^{},{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}𝑹_{y}^{},{}_{}{}^{(j)}𝑹)+({}_{}{}^{(i)}𝑹_{z}^{},{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}𝑹_{z}^{},{}_{}{}^{(j)}𝑹))`$
Using the relation

$${}_{}{}^{(i)}𝑹_{x}^{}=𝒆_1{}_{}{}^{(i)}\alpha {}_{}{}^{(i)}\dot{𝝍}({}_{}{}^{(i)}\dot{𝝍},𝒆_1)$$
(D.5)
and writing

$$({}_{}{}^{(i)}𝑹_{x}^{},{}_{}{}^{(i)}𝑹)=(𝒆_1,{}_{}{}^{(i)}𝑹){}_{}{}^{(i)}\alpha ({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\dot{𝝍},𝒆_1)$$
(D.6)
we obtain

$`\left(({}_{}{}^{(i)}𝑹_{x}^{},{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}𝑹_{x}^{},{}_{}{}^{(j)}𝑹)\right)=\left((𝒆_1,{}_{}{}^{(i)}𝑹)(𝒆_1,{}_{}{}^{(j)}𝑹)\right)`$
$`{}_{}{}^{(j)}\alpha \left((𝒆_1,{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}\dot{𝝍},{}_{}{}^{(j)}𝑹)({}_{}{}^{(j)}\dot{𝝍},𝒆_1)\right)`$
$`{}_{}{}^{(i)}\alpha \left((𝒆_1,{}_{}{}^{(j)}𝑹)({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\dot{𝝍},𝒆_1)\right)+`$
$`{}_{}{}^{(i)}\alpha {}_{}{}^{(j)}\alpha ({}_{}{}^{(i)}(\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\dot{𝝍},𝒆_1)({}_{}{}^{(j)}\dot{𝝍},{}_{}{}^{(j)}𝑹)({}_{}{}^{(j)}\dot{𝝍},𝒆_1)).`$ (D.7)
Thus the brackets in (D.4) are

$`\left(\right)`$ $`=`$ $`({}_{}{}^{(i)}𝑹,{}_{}{}^{(j)}𝑹){}_{}{}^{(j)}\alpha \left(({}_{}{}^{(j)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(j)}\dot{𝝍},{}_{}{}^{(j)}𝑹)\right)`$ (D.8)
$`{}_{}{}^{(i)}\alpha \left(({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}𝑹)({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)\right)+`$
$`{}_{}{}^{(i)}\alpha {}_{}{}^{(j)}\alpha \left(({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}\dot{𝝍})({}_{}{}^{(j)}\dot{𝝍},{}_{}{}^{(j)}𝑹)\right).`$
Consequently, the r.h.s. of the field equation is

$`k\eta ^{ab}f_{,a}f_{,b}`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{n}{}}}{\displaystyle \frac{\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R_{}^{3}}}{1k\frac{{}_{}{}^{(i)}m}{{}_{}{}^{(i)}R}}}{\displaystyle \frac{\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R_{}^{3}}}{1k\frac{{}_{}{}^{(j)}m}{{}_{}{}^{(j)}R}}}[({}_{}{}^{(i)}𝑹,{}_{}{}^{(j)}𝑹)+`$
$`({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}𝑹)({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(i)}𝑹)({}_{}{}^{(i)}\beta {}_{}{}^{(j)}\beta +{}_{}{}^{(i)}\alpha +{}_{}{}^{(j)}\alpha {}_{}{}^{(i)}\alpha {}_{}{}^{(j)}\alpha ({}_{}{}^{(i)}\dot{𝝍},{}_{}{}^{(j)}\dot{𝝍}))]`$


format: text

text:
# 1 Introduction

## 1 Introduction

Observations of double-mode Cepheids provide many information about evolution and structure of massive stars. Parameters of these stars that can be measured with high precision, like periods, period ratios and Fourier coefficients of light curve decomposition, can give us important information not only on the basic physical parameters of these stars but also on chemical composition or distance to them.

Only 14 double-mode Cepheids are known in the Galaxy (Pardo and Poretti 1997). In this sample only one object, CO Aur, pulsates simultaneously in the first and the second overtone (FO/SO) modes. The remaining Galactic double-mode Cepheids have the fundamental mode and the first overtone excited (FU/FO). Cepheids pulsating in two radial modes are difficult to discover, because of the large number of datapoints required to detect and characterize them. Gravitational microlensing surveys provide unique observational material for search for double mode Cepheids. Regular, long term photometric observations of high stellar density regions of the sky, are ideal for examining variable stars, including double-mode Cepheids.

Forty five double-mode Cepheids were identified in the Large Magellanic Clouds by the MACHO microlensing team (Alcock et al. 1995). This sample, increased later to 75 objects (Alcock et al. 1999), allowed detailed study of properties of these objects. Based on the Fourier decomposition of FO/SO double-mode Cepheids light curves, Alcock et al. (1999) characterized properties of the second overtone pulsations, what led to detection of potential candidates for the single-mode second overtone Cepheids (Alcock et al. 1999, Udalski et al. 1999b).

The second phase of the Optical Gravitational Lensing Experiment (OGLE-II) microlensing survey started in January 1997. Since then the Large and Small Magellanic Clouds have been observed regularly, practically on every clear night. Observations were made through the BVI filters, very closely reproducing the standard BVI system. More than three years of monitoring allowed to detect and characterize thousands of Cepheids in both Magellanic Clouds including large sample of double-mode Cepheids.

In the previous papers of this series we presented analysis of 93 double-mode Cepheids detected by the OGLE team in the SMC (Udalski et al. 1999a), discovery of 13 candidates for Cepheids pulsating solely in the second overtone (Udalski et al. 1999b), analysis of the Period-Luminosity and Period-Luminosity-Color relations of Cepheids from the LMC and SMC (Udalski et al. 1999c) and Catalogs of Cepheids in the LMC (Udalski et al. 1999d) and SMC (Udalski et al. 1999e).

This paper, the final of this series, completes the sample of Cepheids detected in Magellanic Clouds during the second phase of the OGLE project. We present here 76 double-mode Cepheids detected in the LMC. Nineteen of these objects are stars pulsating in the fundamental mode and the first overtone, 57 are the first and the second overtone pulsators. Forty two double-mode Cepheids from our sample were previously detected by the MACHO team (Alcock et al. 1995, 1999), the remaining stars are detected for the first time. We provide all basic observational parameters of detected stars. Their photometry is available from the OGLE Internet archive.

## 2 Observations

All observations presented in this paper were carried out during the second phase of the OGLE experiment with the 1.3-m Warsaw telescope at the Las Campanas Observatory, Chile, which is operated by the Carnegie Institution of Washington. The telescope was equipped with the ”first generation” camera with a SITe $`2048\times 2048`$ CCD detector working in the drift-scan mode. The pixel size was 24 $`\mu `$m giving the 0.417 arcsec/pixel scale. Observations of the LMC were performed in the ”slow” reading mode of the CCD detector with the gain 3.8 e-/ADU and readout noise about 5.4 e-. Details of the instrumentation setup can be found in Udalski, Kubiak and Szymański (1997).

Observations of the LMC started on January 6, 1997. 11 driftscan fields, each covering $`14.2\times 57`$ arcmins on the sky, were observed during the first months of 1997. Additional 10 fields were added in October 1997 increasing the total observed area of the LMC to about 4.5 square degree. In this paper we present data collected up to the end of May 2000. Observations were obtained in the standard BVI-bands with the majority of measurements made in the I-band. The reader is referred to Udalski et al. (2000) for more details about methods of data reduction, tests on quality of photometric data, astrometry, location of observed fields etc.

## 3 Selection of Double-Mode Cepheids

The search for Cepheids in the 21 OGLE LMC fields is described in detail in Udalski et al. (1999d). The variable stars classified as Cepheids (about 1500 objects in total) were then searched for double-mode objects.

The selection of double-mode Cepheids was performed using identical two stage algorithm as applied for the SMC double-mode Cepheid search (Udalski et al. 1999a). In the first stage the observations of each Cepheid candidate were folded with the detected period, the light curve was fitted by high order polynomial and subtracted. The residuals were then searched for additional periodic signal and, if detected, such a candidate was selected for further analysis. Then a histogram of the ratio of the shorter to the longer period of double-mode Cepheid candidates was constructed. It exhibited two clear sharp peaks corresponding to the ratio of the first overtone to the fundamental period, $`0.72`$, and the second to the first overtone period, $`0.805`$, in good agreement with Alcock et al. (1995, 1999). The list of the double-mode Cepheid candidates from this search included stars having the period ratio within $`\pm 0.02`$ from these values.

The second, final search for double-mode Cepheids was performed using the Clean algorithm of period determination (Roberts, Lehár and Dreher 1987). All 1500 objects from the LMC Cepheid candidate list were analyzed with the Clean period search algorithm. Having well established limits for the period ratio of double-mode Cepheids from the preliminary analysis, only those objects which exhibited suitable period ratio ($`\pm 0.015`$) between the highest peak in the power spectrum and one of the next four strongest peaks were further analyzed. The final list of the double-mode Cepheid candidates presented in this paper was obtained after careful visual inspection of the Clean power spectra of each object.

## 4 Double-Mode Cepheids in the LMC

Tables 1 and 2 list all double-mode Cepheids detected in the central area of the LMC. They contain 81 entries but only 76 objects. Five stars are located in the overlapping regions between fields and they were discovered independently in each field. Fundamental and first overtone mode pulsators (FU/FO) are listed in Table 1 while Table 2 includes objects pulsating in the first and second overtones (FO/SO). Basic parameters of each star: right ascension and declination (J2000), the intensity-mean I-band magnitude, $`(BV)`$ and $`(VI)`$ colors, both periods and their ratio are provided. Accuracy of periods is about $`410^5P`$. Finding charts for all objects are presented in Appendix A. The size of the I-band subframes is $`60\times 60`$ arcsec; North is up and East to the left.

Appendices B and C show the light curves of FU/FO and FO/SO pulsators, respectively. The first and second columns in each Appendix contain original photometric data folded with the shorter and longer periods while the remaining columns show variability attributed to each mode after subtraction of the other period variability approximated by Fourier series of fifth order. For objects revealing also periodicity equal to the sum and/or difference of both mode frequencies and having an amplitude larger than twice the formal error – such terms were also subtracted from the original data. BVI photometry of all objects is available from the OGLE Internet archive (see Section 6).

Completeness of the sample is determined by completeness of the variable star search in the OGLE databases and efficiency of double-mode Cepheid detection algorithm. OGLE Cepheid catalog of the LMC was estimated to be complete in more than 96% (Udalski et al. 1999d). Completeness of the detection algorithm can be assessed by comparison of results obtained in the preliminary and final (Clean) searches. More than 90% objects in both lists are common suggesting good completeness of the search.

As a test of completeness we cross-identified double mode Cepheids reported by MACHO (Alcock et al. 1995, 1999). 39 out of 42 objects which are located in the OGLE fields were detected during our search. Two of the remaining objects have been misclassified, and have not entered to the list of Cepheids. One star – LMC\_SC4 168269 – has a faint, close companion, and the noise in the power spectrum was greater than the peak of its second pulsation mode. All double-mode Cepheids found by the MACHO team are marked by the letter ’M’ in the last column of Tables 1 and 2.

## 5 Discussion

Seventy six double-mode Cepheids were identified during the presented search in the 4.5 square degree area in the central part of the LMC. Nineteen objects pulsate simultaneously in the fundamental mode and first overtone while 57 objects in the first and second overtones. In the following Subsections we present the basic observational properties of detected sample.

### 5.1 Period Ratio in Double-Mode Cepheids

The ratio of the two periods of double-mode Cepheids is expected to be dependent on metallicity. This dependence is clearly seen in Fig. 1, which presents the period ratio of the FU/FO pulsators plotted as a function of the fundamental mode period. Beside of the LMC double-mode Cepheids we present there the period ratios of FU/FO Cepheids from the SMC (Udalski et al. 1999a) and the Galactic FU/FO Cepheids (Pardo and Poretti 1997). The metal contents of Cepheids in the Galaxy, LMC and SMC are approximately $`Z=0.02`$, 0.008 and 0.004, respectively. It can be seen, that the lower metallicity – the higher the FU/FO period ratio of double-mode Cepheids. The best linear fit for the LMC FU/FO pulsators is given by the equation:

$`P_{\mathrm{FO}}/P_{\mathrm{FU}}=0.734`$ $``$ $`0.035\times \mathrm{log}P_{\mathrm{FU}},`$ (1)
$`0.002`$ $`0.004`$
Period ratios for FO/SO Cepheids in LMC are almost identical with the ratios of the SMC objects. Fig. 2 presents these values in similar fashion as Fig. 1. The only one known FO/SO Galactic Cepheid (CO Aur) is also marked on the diagram. Our sample of the LMC double-mode Cepheids is numerous enough to note that the period ratio-period relation for FO/SO Cepheids is non-linear. The best square fit is as follows:

$$P_{\mathrm{SO}}/P_{\mathrm{FO}}=0.8040.020\times \mathrm{log}P_{\mathrm{FO}}0.067\times \mathrm{log}^2P_{\mathrm{FO}}.$$
$`(2)`$
### 5.2 Fourier Decomposition of Light Curves of Double-Mode Cepheids

Fourier decomposition of light curves of pulsating stars has been widely used for analyzing their properties (Simon and Lee 1981). In the case of Cepheids the ratio of amplitudes of the first harmonic and the fundamental period, $`R_{21}=A_2/A_1`$, and phase difference, $`\varphi _{21}=\varphi _22\varphi _1`$ are particularly useful. Both allow to distinguish between fundamental mode and first overtone pulsators.

Fig. 3 presents the $`R_{21}`$ vs. $`\mathrm{log}P`$ and $`\varphi _{21}`$ vs. $`\mathrm{log}P`$ diagrams constructed for more than 1300 single-mode Cepheids (small dots) taken from the Catalog of Cepheids from the LMC (Udalski et al. 1999d). The $`R_{21}`$ vs. $`\mathrm{log}P`$ diagram shows the characteristic, well separated ”V-shape” sequences for Cepheids pulsating in the fundamental mode and the first overtone. In the similar diagram $`\varphi _{21}`$ vs. $`\mathrm{log}P`$ the sequences for both modes of pulsation are also well defined but the separation is smaller and in some ranges of periods they overlap.

We decomposed the light curves of double-mode Cepheids to the sum of two Fourier series of fifth order corresponding to both periodicities including the terms with periodicities equal to the sum and difference of mode frequencies when their amplitudes were larger than twice the formal errors. Then we calculated $`R_{21}`$ and $`\varphi _{21}`$ for both pulsating modes. They are listed in Tables 1 and 2.

Fig. 3 shows positions of the FU/FO Cepheids on the $`R_{21}`$ vs. $`\mathrm{log}P`$ and $`\varphi _{21}`$ vs. $`\mathrm{log}P`$ diagrams. The values for the fundamental mode pulsation are plotted with large filled dots while for the first overtone mode with open circles. Objects with non-significant first harmonic amplitude, $`A_2`$, (i.e., with almost sinusoidal light curve) have $`R_{21}=0`$ and their $`\varphi _{21}`$ is not defined. Similar diagrams for FO/SO double-mode Cepheids are shown in Fig. 4. Values of $`R_{21}`$ and $`\varphi _{21}`$ for the first overtone pulsation are plotted with large filled dots and with open circles for the second overtone pulsations.

Figs. 3 and 4 are similar to analogous diagrams plotted for the SMC double-mode Cepheids (Udalski et al. 1999a). In both types of double-mode Cepheids the values of $`R_{21}`$ and $`\varphi _{21}`$ of first overtone pulsations fall on the sequences of the single-mode first overtone Cepheids.

Fundamental mode pulsations in the FU/FO Cepheids have much smaller values of the $`R_{21}`$ parameter compared to the single-mode fundamental mode pulsators. They fall practically on the sequence of the first overtone pulsators. This means that the light curves of the fundamental mode pulsations in double-mode Cepheids are more sinusoidal than in the single-mode Cepheids of that type. However, it is clearly seen that $`\varphi _{21}`$ of the fundamental mode pulsations in double-mode Cepheids fall in most cases on the single-mode fundamental mode Cepheid sequence.

Second overtone pulsations in the FO/SO Cepheids are generally small amplitude, almost sinusoidal variations (small $`R_{21}`$). Their $`\varphi _{21}`$ values are usually larger than the first overtone values for similar periodicities (Fig. 4). Identical conclusions on the behavior of second overtone pulsations in double-mode Ceheids were presented by Alcock et al. (1999).

### 5.3 Color-Magnitude Diagram and Colors of Double-Mode Cepheids

Fig. 5 presents the color-magnitude diagram of subfield 2 of the LMC\_SC3 field corrected for the mean $`E(BV)=0.120`$ reddening in this direction (Udalski et al. 1999d). Tiny dots correspond to the field stars. Larger dots indicate positions of the single-mode Cepheids: fundamental mode – darker dots, first overtone – lighter dots. Positions of FU/FO double-mode Cepheids are indicated by large filled circles while FO/SO objects by star symbols.

Fig. 6 shows distribution of color indices $`(VI)_0`$ of single-mode FU and FO Cepheids and double-mode FU/FO and FO/SO pulsators in the LMC. $`E(BV)`$ reddening values from Table 2 of Udalski et al. (1999d) were used to deredden the photometry presented in Tables 1 and 2. The width of the bin is 0.03 mag. Thick solid and dotted lines correspond to the single-mode fundamental and first overtone Cepheids while thin solid and dotted lines to the FU/FO and FO/SO double-mode pulsators. All histograms were fitted with Gaussians which fit well the observed color distributions. In the case of single-mode fundamental mode Cepheids there is a small excess of red objects.

The mean $`(VI)_0`$ color and the standard deviation of its distribution are: (0.425, 0.07), (0.496, 0.08), (0.513, 0.03), (0.594, 0.08) respectively for FO/SO, FO, FU/FO and FU Cepheids. Color indices (temperature) distribution depends on the type of pulsations. The single-mode first overtone Cepheids are on average by about 0.1 mag bluer than the fundamental mode pulsators. As one could expect the FU/FO double-mode Cepheids have $`VI`$ color distribution in between the first and fundamental mode distributions of single-mode stars. The color distribution of FO/SO double-mode Cepheids resembles that of the single-mode first overtone stars but it is shifted bluewards.

## 6 Data Availability

The BVI photometry of the LMC double-mode Cepheids is available to the astronomical community in the electronic form from the OGLE archive:

http://www.astrouw.edu.pl/~ogle

ftp://sirius.astrouw.edu.pl/ogle/ogle2/var\_stars/lmc/cep/dmcep/

or its US mirror

http://bulge.princeton.edu/~ogle
ftp://bulge.princeton.edu/ogle/ogle2//var\_stars/lmc/cep/dmcep/

Acknowledgements. The paper was partly supported by the Polish KBN grant 2P03D00814 to A. Udalski and 2P03D00916 to M. Szymański. Partial support for the OGLE project was provided with the NSF grant AST-9820314 to B. Paczyński.

## REFERENCES

* Alcock, C. et al.  1995, Astron. J., 109, 1652.
* Alcock, C. et al.  1999, Astrophys. J., 511, 185.
* Pardo, I., and Poretti, E. 1997, Astron. Astrophys., 324, 121.
* Roberts, D.H, Lehár, J., and Dreher, J.W. 1987, Astron. J., 93, 968.
* Simon, N.R., and Lee, A.S. 1981, Astrophys. J., 248, 291.
* Udalski, A., Kubiak, M., and Szymański, M. 1997, Acta Astron., 47, 319.
* Udalski, A., Soszyński, I., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1999a, Acta Astron., 49, 1.
* Udalski, A., Soszyński, I., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1999b, Acta Astron., 49, 45.
* Udalski, A., Szymański, M., Kubiak, M., Pietrzyński, G., Soszyński, I., Woźniak, P., and Żebruń, K. 1999c, Acta Astron., 49, 201.
* Udalski, A., Soszyński, I., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1999d, Acta Astron., 49, 223.
* Udalski, A., Soszyński, I., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1999e, Acta Astron., 49, 437.
* Udalski, A., Szymański, M., Kubiak, M., Pietrzyński, G., Soszyński, I., Woźniak, P., and Żebruń, K. 2000, Acta Astron., 50, 307.


text:
# Entangled state based on nonorthogonal state

## I INTRODUCTION

Entanglement and its information theoretic aspects have been studied by many authors . For a pure entangled state of a bipartite system $`|\rho _{AB}`$, a measure of entanglement is defined as

$$E(|\rho _{AB})=\mathrm{Tr}_A\rho _A\mathrm{log}\rho _A,\rho _A=\mathrm{Tr}_B|\rho _{AB}\rho _{AB}|,$$
(1)
which is called as ”entropy of entanglement”. This quantity enjoys two kinds of information theoretic interpretations. One of them is entanglement of formation which means the asymptotic number $`k`$ of standard singlet required to locally prepare faithfully $`n`$ identical copies of a system in bipartite state $`|\rho _{AB}`$ for very large $`k`$ and $`n`$. Other is distillable entanglement which means the asymptotic number of singlets $`k`$ that can be distilled from $`n`$ identical copies of $`|\rho _{AB}`$. In particular, it satisfies

$$\underset{n,k\mathrm{}}{lim}\frac{k}{n}=E(|\rho _{AB}).$$
(2)
Explicit expressions for $`E(|\rho _{AB})`$ is only known in the case of two qubit systems . In fact, it is given as

$$E(|\rho _{AB})=H[\frac{1}{2}(1+\sqrt{1C(\rho _{AB})^2})]$$
(3)
where $`H[x]`$ is the entropy function and $`C(\rho _{AB})`$ is the ”concurrence” defined by $`C(\rho _{AB})=|\rho _{AB}|\stackrel{~}{\rho }_{AB}|`$ with $`|\stackrel{~}{\rho }_{AB}=\sigma |\rho _{AB}^{}`$. The similar analytic formulas for mixed states of qubits is also obtained for the properly defined entanglement of formation . In this paper, we study properties of entangled states based on nonorthogonal states such as coherent states. An implementation scheme for manipulating such states are also discussed.

## II QUASI BELL STATE

### A General definition

Let us define the entangled state based on nonorthogonal states $`|\psi _1`$ and $`|\psi _2`$ such as $`\psi _1|\psi _2=\kappa `$ and $`\psi _2|\psi _1=\kappa ^{}`$. They can be described by

$`\{\begin{array}{ccc}|\mathrm{\Psi }_1\hfill & =& h_1(|\psi _1_A|\psi _2_B+|\psi _2_A|\psi _1_B)\hfill \\ |\mathrm{\Psi }_2\hfill & =& h_2(|\psi _1_A|\psi _2_B|\psi _2_A|\psi _1_B)\hfill \\ |\mathrm{\Psi }_3\hfill & =& h_3(|\psi _1_A|\psi _1_B+|\psi _2_A|\psi _2_B)\hfill \\ |\mathrm{\Psi }_4\hfill & =& h_4(|\psi _1_A|\psi _1_B|\psi _2_A|\psi _2_B)\hfill \end{array}`$ (8)
where $`\{h_i\}`$ are normalized constant:$`h_1=h_3=1/\sqrt{2(1+\kappa ^2)}`$, $`h_2=h_4=1/\sqrt{2(1\kappa ^2)}`$. They are not orthogonal each other. Here, if $`\kappa =\kappa ^{}`$, then the Gram matrix of them becomes very simple as follows:

$$G=\left(\begin{array}{cccc}1& 0& D& 0\\ 0& 1& 0& 0\\ D& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$$
(9)
where $`D=\frac{2\kappa }{1+\kappa ^2}`$. If the basic states are orthogonal, then they are Bell states. Let us discuss the entropy of entanglement for the above states. We, first, calculate the reduced density operators of the quasi Bell state. They are $`\rho _{A}^{}{}_{}{}^{(1)}=\rho _{A}^{}{}_{}{}^{(3)}`$ and $`\rho _{A}^{}{}_{}{}^{(2)}=\rho _{A}^{}{}_{}{}^{(4)}`$. Their concrete forms are

$$\rho _{A}^{}{}_{}{}^{(1)}=\frac{1}{2(1+\kappa ^2)}\{|\psi _1_A\psi _1|+\kappa |\psi _1_A\psi _2|+\kappa |\psi _2_A\psi _1|+|\psi _2_A\psi _2\}$$
(10)
$$\rho _{A}^{}{}_{}{}^{(2)}=\frac{1}{2(1\kappa ^2)}\{|\psi _1_A\psi _1|\kappa |\psi _1_A\psi _2|\kappa |\psi _2_A\psi _1|+|\psi _2_A\psi _2|\}$$
(11)
The eigenvalues of the above density operators $`\rho _{A}^{}{}_{}{}^{(1)}`$(or $`\rho _{A}^{}{}_{}{}^{(3)}`$) are given as follows by using the Gram matrix elements $`C_{ij}=|\mathrm{\Psi }_i|\mathrm{\Psi }_j|`$ of Eq. (9):

$$\lambda _{1/1}=\frac{(1+\kappa )^2}{2(1+\kappa ^2)}=\frac{1+C_{13}}{2},\lambda _{2/1}=\frac{(1\kappa )^2}{2(1+\kappa ^2)}=\frac{1C_{13}}{2}$$
(12)
and for $`\rho _{A}^{}{}_{}{}^{(2)}`$(or $`\rho _{A}^{}{}_{}{}^{(4)}`$), we have

$$\lambda _{1/2}=\frac{1}{2}=\frac{1+C_{24}}{2},\lambda _{2/2}=\frac{1}{2}=\frac{1C_{24}}{2}.$$
(13)
The entropy of entanglement is then

$$E(|\mathrm{\Psi }_1)=E(|\mathrm{\Psi }_3)=\frac{1+C_{13}}{2}\mathrm{log}\frac{1+C_{13}}{2}\frac{1C_{13}}{2}\mathrm{log}\frac{1C_{13}}{2}$$
(14)
and $`E(|\mathrm{\Psi }_2)=E(|\mathrm{\Psi }_4)=1`$. Thus $`|\mathrm{\Psi }_2`$ and $`|\mathrm{\Psi }_4`$ have perfect entanglement, even though the enatangled states consist of nonorthogonal state in each subsystem. These results are true for arbitrary nonorthogonal states with $`\psi _1|\psi _2`$ $`=`$ $`\psi _2|\psi _1`$ $`=`$ $`\kappa `$ and do not depend on the physical dimension of the systems.

### B Generation of quasi Bell states

It is well known that an entangled state can be generated by Walsh-Hadamard gate and CN gate. That is, when the input state for control bit is a superposition state generated by Walsh-Hadamard gate, the output state of the CN gate is an entangled state. The W-H gate(Walsh-Hadamard transformation) is described by

$$U_{WH}=\mathrm{exp}\{\theta (|01||10|)\}$$
(15)
and the unitary operator for the control NOT is

$$U_{CN}=|0_C0|I_T+|1_C1|(|0_T1|+|1_T0|)$$
(16)
The function of the control NOT(CN gate) is as follows:

$`\{\begin{array}{ccc}|0_C|0_T\hfill & & U_{CN}|0_C|0_T=|0_C|0_T\hfill \\ |0_C|1_T\hfill & & U_{CN}|0_C|1_T=|0_C|1_T\hfill \\ |1_C|0_T\hfill & & U_{CN}|1_C|0_T=|1_C|1_T\hfill \\ |1_C|1_T\hfill & & U_{CN}|1_C|1_T=|1_C|0_T\hfill \end{array}`$ (21)
where $`C`$ and $`T`$ mean control mode, and target mode, respectively.

On the two state space spanned by nonorthogonal states: $`|\psi _1`$ and $`|\psi _2`$, we can consider general scheme to manipulate the quasi Bell states. Let us define the orthonormal basis

$`|\psi _e`$ $`=`$ $`(|\psi _1+|\psi _2)/\sqrt{2(1+\kappa )},`$ (22)
$`|\psi _o`$ $`=`$ $`(|\psi _1|\psi _2)/\sqrt{2(1\kappa )}.`$ (23)
They play a role of qubit basis $`\{|0,|1\}`$. We can then go along with quantum logic operations on qubit systems. In terms of the basis defined by Eq. (22), and Eq. (23), the required gates are

$$U_{WH}=\mathrm{exp}\{\theta (|\psi _e\psi _o||\psi _o\psi _e|)\}$$
(24)
$`U_{CN}`$ $`=`$ $`|\psi _e_C\psi _e|I_T`$ (25)
$`+`$ $`|\psi _o_C\psi _o|\left(|\psi _e_T\psi _o|+|\psi _o_T\psi _e|\right)`$ (26)
The W-H gate acts on the input superposition state as follows:

$`|\psi _e_C`$ $``$ $`U_{WH}|\psi _e_C`$ (27)
$`=`$ $`|\psi _e_C+|\psi _o_C`$ (28)
Thus we have

$`K(|\psi _e_C`$ $`+`$ $`|\psi _o_C)|\psi _e_T`$ (29)
$``$ $`U_{CN}K\left(|\psi _e_C+|\psi _o_C\right)|\psi _e_T`$ (30)
$`=`$ $`K\left(|\psi _e_C|\psi _e_T+|\psi _o_C|\psi _o_T\right)`$ (31)
$`=`$ $`h_3\left(|\psi _1_C|\psi _1_T+|\psi _2_C|\psi _2_T\right)`$ (32)
where $`K`$ is the normalized constant. The final state is one of quasi Bell state. Thus we have quasi Bell states based on such operations. If we use the coherent states as the basic states, then the above gates correspond to bosonic gates whose realization is discussed in the later section.

### C General case

Here let us consider the general pure entangled state of nonorthogonal state. They can be described as follows:

$`\{\begin{array}{ccc}|\mathrm{\Psi }_1\hfill & =& g_1(\beta |\psi _1_A|\psi _2_B+\sqrt{1\beta ^2}|\psi _2_A|\psi _1_B)\hfill \\ |\mathrm{\Psi }_2\hfill & =& g_2(\beta |\psi _1_A|\psi _2_B\sqrt{1\beta ^2}|\psi _2_A|\psi _1_B)\hfill \\ |\mathrm{\Psi }_3\hfill & =& g_3(\beta |\psi _1_A|\psi _1_B+\sqrt{1\beta ^2}|\psi _2_A|\psi _2_B)\hfill \\ |\mathrm{\Psi }_4\hfill & =& g_4(\beta |\psi _1_A|\psi _1_B\sqrt{1\beta ^2}|\psi _2_A|\psi _2_B)\hfill \end{array}`$ (37)
where $`g_i`$ is normalized constant, and $`\beta `$ is real number. Since all the elements of the Gram matrix is not zero, they are not orthogonal states. The reduced density operators in this case become $`\rho _{A}^{}{}_{}{}^{(1)}=\rho _{A}^{}{}_{}{}^{(3)}`$ and $`\rho _{A}^{}{}_{}{}^{(2)}=\rho _{A}^{}{}_{}{}^{(4)}`$, and

$`\rho _A^{(1)}`$ $`=`$ $`k_1\{\beta ^2|\psi _1_A\psi _1|+\kappa \beta \sqrt{1\beta ^2}|\psi _1_A\psi _2`$ (38)
$`+`$ $`\kappa \beta \sqrt{1\beta ^2}|\psi _2_A\psi _1|+(1\beta ^2)|\psi _2_A\psi _2|\}`$ (39)
$`\rho _A^{(2)}`$ $`=`$ $`k_2\{\beta ^2|\psi _1_A\psi _1|\kappa \beta \sqrt{1\beta ^2}|\psi _1_A\psi _2|`$ (40)
$``$ $`\kappa \beta \sqrt{1\beta ^2}|\psi _2_A\psi _1|+(1\beta ^2)|\psi _2_A\psi _2|\}`$ (41)
where $`k_i=\frac{1}{1\pm 2\kappa ^2\beta \sqrt{1\beta ^2}}`$ is normalized constant. The analysis in this case is also easy. Here we again assumed that $`\psi _1|\psi _2`$ $`=`$ $`\psi _2|\psi _1`$ $`=`$ $`\kappa `$.

## III QUASI BELL STATES OF COHERENT STATES

Let us consider the binary coherent states of a bosonic mode $`\{|\alpha ,|\alpha \}`$, where $`(\kappa =\alpha |\alpha =\mathrm{e}^{2|\alpha |^2})`$. Then the quasi Bell states are

$`\{\begin{array}{ccc}|\mathrm{\Psi }_1\hfill & =& h_1(|\alpha _A|\alpha _B+|\alpha _A|\alpha _B)\hfill \\ |\mathrm{\Psi }_2\hfill & =& h_2(|\alpha _A|\alpha _B|\alpha _A|\alpha _B)\hfill \\ |\mathrm{\Psi }_3\hfill & =& h_3(|\alpha _A|\alpha _B+|\alpha _A|\alpha _B)\hfill \\ |\mathrm{\Psi }_4\hfill & =& h_4(|\alpha _A|\alpha _B|\alpha _A|\alpha _B)\hfill \end{array}`$ (46)
where $`\alpha `$ is coherent amplitude of light field. The average photon numbers of the reduced states are

$$n_A^{(1)}=\frac{(1\kappa ^2)}{(1+\kappa ^2)}|\alpha |^2,n_A^{(2)}=\frac{(1+\kappa ^2)}{(1\kappa ^2)}|\alpha |^2$$
(47)
Thus the quasi Bell states can have arbitrary photon, and approach to the Bell states as $`|\alpha |\mathrm{}`$. We mention the characteristic function of quasi Bell states defined as

$`C(\xi ,\eta )`$ $`=`$ $`\mathrm{Tr}\left[|\mathrm{\Psi }\mathrm{\Psi }|\mathrm{exp}(\xi \mathrm{a}_A^{})\mathrm{exp}(\xi ^{}\mathrm{a}_A)\mathrm{exp}(\eta \mathrm{a}_B^{})\mathrm{exp}(\eta ^{}\mathrm{a}_B)\right]`$ (48)
$`\times `$ $`\mathrm{exp}\{(|\xi |^2+|\eta |^2)/2\}`$ (49)
where $`a`$ and $`a^{}`$ are the annihilation and creation operators, respectively. They are actually

$`C(\xi ,\eta |i=1,2)`$ $`=`$ $`h_{i}^{}{}_{}{}^{2}\mathrm{exp}\{(|\xi |^2+|\eta |^2)/2\}\{\mathrm{exp}(A_1B_1)\alpha `$ (50)
$`+`$ $`\mathrm{exp}(A_1+B_1)\alpha \pm \mathrm{exp}(A_2B_2)\alpha `$ (51)
$`\pm `$ $`\mathrm{exp}(A_2+B_2)\alpha \}`$ (52)
$`C(\xi ,\eta |i=3,4)`$ $`=`$ $`h_{i}^{}{}_{}{}^{2}\mathrm{exp}\{(|\xi |^2+|\eta |^2)/2\}\{\mathrm{exp}(A_1+B_1)\alpha `$ (53)
$`+`$ $`\mathrm{exp}(A_1B_1)\alpha \pm \mathrm{exp}(A_2+B_2)\alpha `$ (54)
$`\pm `$ $`\mathrm{exp}(A_2B_2)\alpha \}`$ (55)
where $`A_1=(\xi \xi ^{}),A_2=(\xi +\xi ^{}),B_1=(\eta \eta ^{}),B_2=(\eta +\eta ^{})`$. It is worthy to mention that the quasi Bell states do not belong the Gaussian state in contrast to that two mode squeezed state does so.

## IV PHYSICAL REALIZATION

In order to manipulate the quasi Bell states of bosonic coherent states, one needs quantum gates acting on a state space spanned by the relevant coherent states. A convenient basis is the even and odd coherent states. Let us denote them as $`\{|e,|o\}`$ hereafter. It would be much difficult to realize quantum gates for these macroscopic qubits. Cochrane, Milburn, and Munro proposed a physical model for such gates . In their model, a CN gate is made by applying an $`H`$ gate to the mode $`a`$, the target bit, then coupling the target and the control (the mode $`b`$), and finally applying another $`H`$ gate to the target again. This CN gate operation is actually valid in a certain limited case of the coherent state amplitude $`\alpha `$ and the classical field amplitude $`\gamma `$. However, their model is indeed indicative. If the $`H`$ gate is capable of generating the Schrödinger cat state of coherent state, that is, includes an appropriate nonlinear Hamiltonian instead of the linear interaction of Eq. , then their scheme provides the universal bosonic gates. (This was also recognized by Tatsuta et al. ) Concerning the two bit interaction, the cross Kerr effect always suffices. Thus the problem is to synthesize the appropriate nonlinear Hamiltonian.

One bit gate operations for macroscopic bosonic qubits can be represented by rotations on the two-state space. For example, the Hadamard gate is represented by

$$\widehat{U}_H=i\mathrm{exp}[i\frac{\pi }{2}\widehat{Q}]\mathrm{exp}[\frac{\pi }{4}\widehat{P}],$$
(56)
$$\widehat{P}=|eo||oe|,\widehat{Q}=|ee||oo|,$$
(57)
The physical process corresponding to $`\widehat{U}_H`$ includes essentially multiphoton nonlinear precess. The corresponding Hamiltonians were studied by Sasaki and Hirota .

The case of our interest is when the amplitude $`|\alpha |`$ is small, in which the assumption used in . Then the Hamiltonian $`\widehat{P}`$ and $`i\widehat{Q}`$ can be effected by the nonlinear Hamiltonian including finite number of nonlinearity. For mathematical convenience, we consider the Hamiltonian for $`\stackrel{~}{U}_H=\widehat{D}(\alpha )\widehat{U}_H\widehat{D}(\alpha )`$. First we introduce a cut off photon number $`M`$ for a weak coherent state such that its photon number distribution in $`n>M`$ becomes negligibly small. Second define

$$\widehat{P}_M=2\sqrt{1c_0^2}(\underset{l=0}{\overset{M}{}}\frac{(\widehat{a}^{})^l\widehat{a}^l}{l!}\underset{n=1}{\overset{M}{}}d_n\frac{\widehat{a}^n}{\sqrt{n!}}\mathrm{h}.\mathrm{c}.),$$
(58)
$$\widehat{Q}_M=4c_0\underset{l=0}{\overset{M}{}}\frac{(\widehat{a}^{})^l\widehat{a}^l}{l!}+2\sqrt{1c_0^2}(\underset{l=0}{\overset{M}{}}\frac{(\widehat{a}^{})^l\widehat{a}^l}{l!}\underset{n=1}{\overset{M}{}}d_n\frac{\widehat{a}^n}{\sqrt{n!}}+\mathrm{h}.\mathrm{c}.),$$
(59)
where

$$c_n=\mathrm{e}^{2\alpha ^2}\frac{(2\alpha )^2}{\sqrt{n!}},d_n=\frac{c_n}{\sqrt{_{n=1}^M}c_n^2}.$$
(60)
Then $`\stackrel{~}{U}_H`$ can be represented

$$\stackrel{~}{U}_H=i\mathrm{exp}[i\frac{\pi }{2}\widehat{Q}_M]\mathrm{exp}[\frac{\pi }{4}\widehat{P}_M]+O(\delta _M),$$
(61)
where

$$\delta _M=1\frac{\underset{n=1}{\overset{M}{}}c_n^2}{1c_0^2}.$$
(62)
The Hamiltonian of $`\widehat{P}_M`$ and $`\widehat{Q}_M`$ still seem to be unrealistic. One possible way to make realistic is to decompose them into a cascade process of lower order nonlinear processes. In fact, as suggested by Harel and Akulin , and Lloyd and Braunstein, it is possible in principle to synthesize the required unitary dynamics by lower order nonlinear Hamiltonians. In particular, it can be shown that the nonlinearity up to third order and the cross Kerr nonlinearity suffice to implement the universal gates for macroscopic bosonic qubits.

## V CONCLUSION

We investigated properties of entangled states based on nonorthogonal state such as coherent states. Implementation of quantum gates for such macroscopic qubits was suggested. We would like to find more simple generation method of such a quasi Bell states.

###### Acknowledgements.

We are grateful to M.Ban, S.Barnett, C.Bennett, S. J. van Enk, C.Fuchs, A.Holevo, K.Kato, R.Jozsa, M.Osaki, and P.Shor for helpful discussions.