text:
# Contents

## 1 Introduction

The orbital path of a test particle freely falling in the gravitational field of a central body of mass $`M`$ and angular momentum $`J`$ is affected by various post-Newtonian general relativistic effects of order $`𝒪(c^2)`$ like the Einstein gravitoelectric secular precession of the argument of pericentre111For an explanation of the Keplerian orbital elements of the orbit of a test particle see, e.g.,  and Figure 1. $`\omega `$  and the Lense-Thirring gravitomagnetic secular precessions of the longitude of the ascending node $`\mathrm{\Omega }`$ and the argument of the pericentre $`\omega `$ .

The Einstein precession is now known at a $`10^4`$ level of relative accuracy, e.g., from Solar System planetary motions by analyzing more than 280,000 position observations (1913-2003) of different types including radiometric observations of planets and spacecrafts, CCD astrometric observations of the outer planets and their satellites, meridian transits and photographic observations . On the contrary, to date there are not yet clear and undisputable direct tests of the various effects induced by the gravitomagnetic field. In April 2004 the Gravity Probe B (GP-B) mission  has been launched. Its goal is to measure, among other things, the gravitomagnetic precessions of the spins of four superconducting gyroscopes  carried onboard with a claimed accuracy of 1$`\%`$ or better. The Lense-Thirring effect on the orbits of the LAGEOS and LAGEOS II Earth’s artificial satellites has been experimentally checked for the first time some years ago with a claimed accuracy of the order of 20$`\%`$ , but, at present, there is not yet a full consensus on the real accuracy obtained in that test .

Recently, it has been shown that also the orbital period of a test particle is affected by the post-Newtonian gravitoelectromagnetic forces. The gravitomagnetic correction to the period of the (equatorial) right ascension222It is nothing but the azimuth angle $`\varphi `$ of a spherical coordinate system in a frame whose origin is in the center of mass of Earth, the $`\{x,y\}`$ plane coincides with the Earth equatorial plane and the $`x`$ axis points toward the Vernal Equinox $``$ (see Figure 1). In satellite dynamics it is one of the direct observable quantities .  $`\alpha `$, for $`e=0`$, $`i=0^{}`$, where $`e`$ and $`i`$ are the eccentricity and the inclination, respectively, of the test particle’s orbit, has been worked out in . The case of more general orbits has been treated in  ($`e=0`$, $`i0^{}`$) and  ($`e0`$, $`i0^{}`$). The gravitoelectric correction has been treated in . Let us consider a pair of counter-orbiting satellites, conventionally denoted as $`(+)`$ and $`()`$, following identical circular and equatorial orbits along opposite directions around a central spinning body of mass $`M`$ and proper angular momentum $`J`$; while the gravitoelectric correction to the Keplerian period has the same sign for both the directions of motion, similar to the classical perturbing corrections of gravitational origin (oblateness of the central mass, tides, N-body interactions), the gravitomagnetic one changes its sign if the motion is reversed. This allows, at least in principle, to single out the latter one by measuring the difference of the periods of the counter-orbiting satellites. Some preliminary error analyses can be found in .

As we will show here, it turns out that also the mean longitudes333$``$ is the mean anomaly of the test particle’s orbit. $`l`$, along with the Keplerian orbital elements, comes out from the machinery of the data reduction process performed by the orbit determination softwares like GEODYN II and UTOPIA. $`l=+\mathrm{\Omega }+\omega `$ are affected, among other perturbations, by the post-Newtonian general relativistic gravitoelectromagnetic forces. Also in this case the gravitomagnetic correction is sensitive to the direction of motion, contrary to the gravitoelectric one.

In this paper we will investigate the possibility of singling out the gravitomagnetic effect on the mean longitudes by means of a suitable space-based mission in the Earth space environment. We will mainly deal with a number of competing classical effects which could represent a source of systematic errors and bias. No measurements error budget is carried out.

## 2 The Gravitomagnetic Effect on the Mean Longitude

From the Gauss perturbative equations  it turns out that the rate equation for $`l`$, for small but finite values of $`e`$ and $`i`$, is

$$\frac{dl}{dt}n\frac{2}{na}A_R\left(\frac{r}{a}\right)+\left(\frac{i^2}{2}\right)\frac{d\mathrm{\Omega }}{dt}+\left(\frac{e^2}{2}\right)\frac{d\omega }{dt},$$
(1)
where $`a`$ is the satellite’s semimajor axis and $`n=(GMa^3)^{\frac{1}{2}}`$ is the Keplerian mean motion; $`A_R`$ is the radial component of the perturbing acceleration $`𝐚_{\mathrm{pert}}`$. The currently available technologies allows to insert Earth artificial satellites in orbits with444This value holds for initial orbital injection. Later, various perturbations, mainly due to the terrestrial gravitational field, affect $`e`$ which can reach values up to $`5\times 10^3`$.  $`e10^3`$ and $`i0.01^{}`$; then, the third and fourth terms in eq.(1), which are proportional to $`i`$ and $`e`$ because $`d\mathrm{\Omega }/dt1/\mathrm{sin}i`$ and $`d\omega /dt1/e`$, can be neglected and the equation

$$\frac{dl}{dt}n\frac{2}{na}A_R$$
(2)
can be used instead of eq.(1) to a good level of approximation.

According to , the radial components of the post-Newtonian general relativistic gravitoelectromagnetic accelerations, for generic orbits around a central spinning body, are

$`A_R^{(\mathrm{GE})}`$ $`=`$ $`{\displaystyle \frac{(GM)^2(1+e\mathrm{cos}f)^2}{c^2a^3(1e^2)^3}}(3+2e\mathrm{cos}fe^2+4e^2\mathrm{sin}^2f),`$ (3)
$`A_R^{(\mathrm{GM})}`$ $`=`$ $`{\displaystyle \frac{2nGJ\mathrm{cos}i}{c^2a^2(1e^2)^{\frac{7}{2}}}}(1+e\mathrm{cos}f)^4,`$ (4)
where $`f`$ is the true anomaly. They are induced by the post-Newtonian general relativistic gravitoelectromagnetic fields555In the weak-field and slow-motion approximation of the General Theory of Relativity the equations of motion of a test particle are
$$\frac{d^2𝐫}{dt^2}=𝐄_\mathrm{g}𝐄^{(\mathrm{GE})}2\frac{𝐯}{c}\times 𝐁^{(\mathrm{GM})},$$
(5) where $`𝐄_\mathrm{g}=GM𝐫/r^3`$ is the usual Newtonian monopole term.

$`𝐄^{(\mathrm{GE})}`$ $`=`$ $`{\displaystyle \frac{GM}{c^2r^3}}\left({\displaystyle \frac{4GM}{r}}𝐯^2\right)𝕣{\displaystyle \frac{4GM}{c^2r^2}}(𝕣𝐯)𝐯,`$ (6)
$`𝐁^{(\mathrm{GM})}`$ $`=`$ $`{\displaystyle \frac{GJ}{cr^3}}[\widehat{𝐉}3(\widehat{𝐉}\widehat{𝐫})\widehat{𝐫}].`$ (7)
From eq.(3)-eq.(4) it can be noted that, while the gravitoelectric acceleration is insensitive to the sense of motion of the test particle along its orbit, it is not so for the gravitomagnetic acceleration due to its dependence on $`\mathrm{cos}i`$. Indeed, according to

$`\mathrm{sin}i`$ $`=`$ $`{\displaystyle \frac{(h_x^2+h_y^2)^{\frac{1}{2}}}{\left|𝐡\right|}},`$ (8)
$`\mathrm{cos}i`$ $`=`$ $`{\displaystyle \frac{h_z}{\left|𝐡\right|}},`$ (9)
where $`𝐡`$ is the orbital angular momentum per unit mass whose components change sign when $`𝐯𝐯`$. By reversing the sign of the velocity vector one obtains ($`+`$ and $``$ denote the pro- and retrograde orbits, respectively)

$`\mathrm{sin}i^{()}`$ $`=`$ $`\mathrm{sin}i^{(+)},`$ (10)
$`\mathrm{cos}i^{()}`$ $`=`$ $`\mathrm{cos}i^{(+)},`$ (11)
from which it follows

$$i^{()}=180^{}i^{(+)}.$$
(12)
By inserting eq.(3)-eq.(4) in eq.(2), it can be obtained for the time $`P_l`$ required to $`l`$ for passing from 0 to 2$`\pi `$

$`P_l^{(\pm )}:=P_l^{(0)}+P_l^{(\mathrm{GE})}+P_l^{(\mathrm{GM})}={\displaystyle \frac{2\pi }{n}}+{\displaystyle \frac{12\pi (GMr_0)^{\frac{1}{2}}}{c^2}}\pm {\displaystyle \frac{8\pi J}{c^2M}}.`$ (13)
Eq.(13) yields for identical orbits followed in opposite directions,

$$\mathrm{\Delta }P_l:=P_l^{(+)}P_l^{()}=\frac{16\pi J}{c^2M},$$
(14)
which amounts to $`5\times 10^7`$ s for Earth and is four times larger than the corresponding effect for the right ascension.

## 3 Sources of Systematic Errors

Eq.(14) would be valid only if the Keplerian periods (and the various classical and post-Newtonian gravitoelectric perturbative corrections to them) of the two satellites were exactly equal; this condition, however, cannot be achieved due to the unavoidable orbital injection errors. Then, eq.(14) has also to account for the difference induced in the Keplerian periods and the various classical and post-Newtonian gravitoelectric corrections to them, e.g., by the difference $`d`$ in the semimajor axes of the two satellites. The differences in $`i`$, which do not affect the Keplerian periods and the post-Newtonian gravitoelectric correction, would have, instead, an impact on the classical perturbative terms which cannot be neglected.

In general, the difference between the mean longitude periods of the two counter-rotating satellites can be written as

$$\mathrm{\Delta }P_l=\mathrm{\Delta }P^{(0)}+\mathrm{\Delta }P_l^{\mathrm{class}}+\mathrm{\Delta }P_l^{(\mathrm{GE})}+\mathrm{\Delta }P_l^{(\mathrm{GM})}.$$
(15)
Over many orbital revolutions, say $`N`$, the accuracy in determining the gravitomagnetic time shift becomes

$$\delta [\mathrm{\Delta }P^{(\mathrm{GM})}]\frac{\delta (\mathrm{\Delta }P)^{\mathrm{exp}}}{N}+\delta [\mathrm{\Delta }P^{(0)}]+\delta [\mathrm{\Delta }P_l^{\mathrm{class}}]+\delta [\mathrm{\Delta }P_l^{\mathrm{GE}}],$$
(16)
where $`\delta (\mathrm{\Delta }P_l)^{\mathrm{exp}}`$ is the experimental error in the difference of the obtained $`P_l^{(\pm )}`$ over $`N`$ revolutions. After a sufficiently high number of orbital revolutions it should be possible to make the term $`\delta (\mathrm{\Delta }P_l)^{\mathrm{exp}}/N`$ smaller than $`\mathrm{\Delta }P_l^{(\mathrm{GM})}`$. In order to get an estimate, let us calculate the angular shift corresponding to $`\mathrm{\Delta }P_l^{(\mathrm{GM})}`$ over an orbital revolution by using

$$\frac{\mathrm{\Delta }l}{\overline{l}}\frac{\mathrm{\Delta }P_l^{(\mathrm{GM})}}{P_l^{(0)}}.$$
(17)
For666In order just to fix the ideas, we shall consider the semimajor axis of the existing SLR ETALON satellites. We will see that high altitudes contribute to reduce the impact of many aliasing perturbations.  $`r_0=25498`$ km and $`P_l^{(0)}=4.05200895378\times 10^4`$ s, $`\overline{l}=2\pi `$ the gravitomagnetic shift amounts to $`\mathrm{\Delta }l=1\times 10^2`$ milliarcseconds (mas in the following; 1 mas = 4.8$`\times 10^9`$ rad), over an angular span. Since the accuracy with which the orientation of the axes of the International Terrestrial Reference777Its realization, the International Terrestrial Reference Frame (ITRF), is obtained by estimates of the coordinates and velocities of a set of observing stations on the Earth . System (ITRS) is of the order888See on the Internet http://hpiers.obspm.fr/iers/bul/bulb/explanatory.html of 3 mas, after 300 orbital revolutions or 140 days, the gravitomagnetic time shift would reach this sensitivity cutoff.

### 3.1 The Impact of the Imperfect Cancellation of the Keplerian Periods

As we will see, the major limiting factor in measuring the gravitomagnetic time shift of interest is the difference of the Keplerian orbital periods

$$\mathrm{\Delta }P^{(0)}=\frac{2\pi [r_0^{\frac{3}{2}}(r_0+d)^{\frac{3}{2}}]}{(GM)^{\frac{1}{2}}}\frac{3\pi dr_0^{\frac{1}{2}}}{(GM)^{\frac{1}{2}}},$$
(18)
where $`r_0`$ represents the nominal value of the semimajor axis of the two satellites. The shift of eq.(18) cannot be made smaller than eq.(14), in practice, by choosing suitably the orbital geometry of the satellites. Indeed, from eq.(18) follows that the relation $`dr_{0}^{}{}_{}{}^{\frac{1}{2}}1.059\times 10^3`$ cm$`^{\frac{3}{2}}`$ should be fulfilled; for $`r_010^9`$ cm it would imply $`d3\times 10^2`$ cm. Then, in order to be able to measure the relativistic effect of interest, which accumulates during the orbital revolution, the difference should be subtracted from the data provided that its error $`\delta [\mathrm{\Delta }P^{(0)}]`$, which is present at every orbital revolution and is due to the uncertainties in the Earth $`GM`$, $`r_0`$ and $`d`$, is smaller than the gravitomagnetic time shift. This error is given by

$$\delta [\mathrm{\Delta }P^{(0)}]\left|\frac{3\pi dr_0^{\frac{1}{2}}}{2(GM)^{\frac{3}{2}}}\right|\times \delta (GM)+\left|\frac{3\pi d}{2(GMr_0)^{\frac{1}{2}}}\right|\times (\delta r_0)+\left|\frac{3\pi r_0^{\frac{1}{2}}}{(GM)^{\frac{1}{2}}}\right|\times (\delta d)$$
(19)
and by assuming $`r_0=25498`$ km, $`(\delta r_0)^{\mathrm{exp}}=1`$ cm and $`\delta (GM)=8\times 10^{11}`$ cm3 s-2 , eq.(19) yields

$$\delta [\mathrm{\Delta }P^{(0)}](2.8\times 10^{14}\mathrm{cm}^1\mathrm{s})\times d+(2.3\times 10^5\mathrm{cm}^1\mathrm{s})\times (\delta d).$$
(20)
Eq.(20) tells us that the error due to the uncertainty in999Note that the uncertainty in the Earth’s $`GM`$ is a not a totally negligible limiting factor. Indeed, for an intersatellite separation of a few kilometers the induced bias is of the order of $`10^8`$ s. If only one satellite was considered the systematic error in the Keplerian period due to $`GM`$ $`\pi (a/GM)^{\frac{3}{2}}\delta (GM)`$ would amount to $`4\times 10^5`$ s, as pointed out in the preliminary evaluations of .  $`GM`$ and, especially, $`r_0`$ is less than the gravitomagnetic effect, while $`\delta d`$ should be at the level of $`2\times 10^2`$ cm. It seems to be impossible to meet this very stringent requirement with the current SLR technology due to many measurement errors (station errors, random errors in precession, nutation and Earth rotation, observation errors). However, an intersatellite tracking approach could yield better results. A level of $`(\delta d)^{\mathrm{exp}}`$ of the order of $`10^3`$ cm is currently available with the K-band Ranging (KBR) intersatellite tracking technology used for the GRACE mission . On the other hand, it must be noted that it would not be easy to use this approach for satellites following opposite directions. The GRACE twins revolve in the same sense.

In regard to $`(\delta d)^{\mathrm{systematic}}`$, thinking about a pair of completely passive, spherical LAGEOS–like satellites, it should be pointed out that for orbits with $`ei0`$, for which eq.(14) holds, many non–gravitational perturbations101010On the semimajor axis there are no long–term gravitational perturbations. affecting the semimajor axis vanish ; the remaining ones could be strongly constrained by constructing the two satellites very carefully with regard to their geometrical and physical properties. Some figures will be helpful. We will adopt the physical properties of the existing LAGEOS satellites. The atmospheric drag induces a decrease in the semimajor axis of $`\mathrm{\Delta }a=a^2\varrho C_DS/m`$ over one orbital revolution, where $`C_D`$ is the satellite drag coefficient, $`S/m`$ is the satellite area-to-mass ratio and $`\varrho `$ is the atmospheric density at the satellite altitude. For $`C_D4.9`$, $`S/m=7\times 10^3`$ cm2 g-1 and $`\varrho =8.4\times 10^{21}`$ g cm-3 (estimate of the atmospheric neutral density at LAGEOS altitude; for $`r_0=25498`$ km it should be smaller, of course) the decrease in the semimajor axis amounts to 4$`\times 10^4`$ cm over one orbital revolution. The impact of direct solar radiation pressure on the semimajor axis of a circular orbit vanishes, also when eclipses effects are accounted for. The nominal Poynting-Robertson effect on the semimajor axis of a circular, equatorial orbit ($`r_0=25498`$ km) is of the order of $`10^3`$ cm over one orbital revolution. However, it is linearly proportional to the satellite reflectivity coefficient $`C_R`$ for which a 0.5–1$`\%`$ mismodelling can be assumed. The semimajor axis of a circular orbit is neither affected neither the albedo over one orbital revolution. The same holds also for the IR Earth radiation pressure and the solar Yarkovsky-Schach effect (by neglecting the eclipses effects). The Earth Yarkovsky-Rubincam effect would affect the semimajor axis of an orbit with $`i0^{}`$ by means of a nominal secular trend less than $`1\times 10^1`$ cm ($`r_0=25498`$ km) over one orbital revolution111111It would vanish if the satellite spin axis was aligned with the $`z`$ axis of an Earth-centered inertial frame.. However, a 20-25$`\%`$ mismodelling can be assumed on it. The direct effect of the terrestrial magnetic field on the semimajor axis of a circular orbit vanishes over one orbital revolution.

Another possibility could be the use of more complex (and more expensive) spacecraft endowed with a drag–free apparatus; it would be helpful in suitably constraining $`(\delta r_0)^{\mathrm{systematic}}`$ and, especially, $`(\delta d)^{\mathrm{systematic}}`$. Indeed, for orbits with semimajor axis of the order of 109 cm the unperturbed Keplerian period is of the order of $`P^{(0)}4\times 10^4`$ s. If the maximum admissible error in knowing the position of the satellites is of the order of $`2\times 10^2`$ cm then the maximum disturbing acceleration that would not mask the gravitomagnetic clock effect, over one orbital revolution, is $`\delta d/[P^{(0)}]^21\times 10^{11}`$ cm s-2. But the drag–free technologies currently under development for LISA  and OPTIS  missions should allow to cancel out the accelerations of non–gravitational origin down to $`3\times 10^{13}`$ cm s-2 Hz$`^{\frac{1}{2}}`$ and 10-12 cm s-2 Hz$`^{\frac{1}{2}}`$, respectively. In our case, for $`r_0=25498`$ km and an orbital frequency of $`2.4\times 10^5`$ Hz, a $`1\times 10^9`$ cm s-2 Hz$`^{\frac{1}{2}}`$ drag-free level would be sufficient. More precisely, from the Gauss perturbative rate equation of the semimajor axis it turns out that, for circular orbits, the effects induced on $`r_0`$ by periodically time-varying non–gravitational along-track accelerations with frequencies multiple of the orbital frequency vanish over one orbital revolution. For constant along-track accelerations $`A_T`$ the situation is different: over one orbital revolution (corresponding to a frequency of the order of 10-5 Hz) their effect on $`r_0`$ would amount to

$$\delta r_0^{(\mathrm{NG})}=4\pi \left(\frac{r_0^3}{GM}\right)A_T.$$
(21)
In order to make $`\delta r_0^{(\mathrm{NG})}=2\times 10^2`$ cm the maximum value for $`A_T`$ would be $`3.8\times 10^{11}`$ cm s-2 in a frequency range with $`10^5`$ Hz as upper limit. It is a difficult but not impossible limit to be obtained with the drag-free technologies which are currently under development.

### 3.2 The Impact of the Imperfect Cancellation of the Post-Newtonian Gravitoelectric Periods

The perturbative correction to $`P_l`$ induced by the post-Newtonian gravitoelectric acceleration, given in eq.(14), amounts to

$$P_l^{(\mathrm{GE})}=\frac{12\pi (GMr_0)^{\frac{1}{2}}}{c^2}.$$
(22)
For an Earth high-orbit satellite it is of the order of 10-5 s. The orbital injection errors in $`r_0`$ would yield

$$\mathrm{\Delta }P_l^{(\mathrm{GE})}\frac{6\pi }{c^2}\left(\frac{GM}{r_0}\right)^{\frac{1}{2}}d,$$
(23)
which amounts to 10-9 s for $`r_0=25498`$ km and $`d=5`$ km. This result justifies, a posteriori, our choice of neglecting the small corrections to eq.(22) induced by the small, but finite, eccentricity of the orbit.

### 3.3 The Impact of the Classical Gravitational Perturbations

In this section the perturbations of gravitational origin on the mean longitude $`l`$ and their impact on the measurement of $`\mathrm{\Delta }P_l^{(\mathrm{GM})}`$ are investigated.

#### 3.3.1 The Geopotential

Since we are interested in effects which are averaged over many orbital revolutions, only the zonal harmonics of geopotential, which induces secular perturbations on $`\mathrm{\Omega },\omega `$ and $``$, will be considered. In the following the Gaussian approach together with eq.(2) will be used.

Let us write the zonal part of geopotential as

$$\mathrm{\Phi }^{\mathrm{zon}}:=\mathrm{\Phi }^{(0)}+\mathrm{\Delta }\mathrm{\Phi }^{\mathrm{zon}}=\frac{GM}{r}GM\underset{\mathrm{}=2}{\overset{\mathrm{}}{}}\frac{R_{}^{\mathrm{}}}{r^{\mathrm{}+1}}J_{\mathrm{}}P_{\mathrm{}}(\mathrm{cos}\theta ).$$
(24)
According to (6.98b) of , the radial component of the perturbing acceleration, in spherical coordinates, for an even zonal harmonic of degree $`\mathrm{}`$, is

$$A_R^{(J_{\mathrm{}})}=\frac{\mathrm{\Delta }\mathrm{\Phi }^{\mathrm{zon}}}{r}=(\mathrm{}+1)\frac{GMR_{}^{\mathrm{}}}{r^{\mathrm{}+2}}J_{\mathrm{}}P_{\mathrm{}}(\mathrm{cos}\theta ).$$
(25)
It is easy to see that the odd zonal harmonics of geopotential do not affect the mean longitude of a satellite in an equatorial orbit. Indeed, it is well known that the odd-degree Legendre polynomials vanish for $`\theta =90^{}`$ . The perturbing radial accelerations due to the even zonal harmonics are

$`A_R^{(J_2)}`$ $`=`$ $`{\displaystyle \frac{3GMR_{}^2J_2}{2r^4}},`$ (26)
$`A_R^{(J_4)}`$ $`=`$ $`{\displaystyle \frac{15GMR_{}^4J_4}{8r^6}},`$ (27)
$`A_R^{(J_6)}`$ $`=`$ $`{\displaystyle \frac{35GMR_{}^6J_6}{16r^8}},`$ (28)
$`A_R^{(J_8)}`$ $`=`$ $`{\displaystyle \frac{315GMR_{}^8J_8}{128r^{10}}}.`$ (29)
Inserting eq.(26)-eq.(29) in eq.(2) yields

$$P_l^{\mathrm{obl}}=\frac{6\pi R_{}^2J_2}{(GMr_0)^{\frac{1}{2}}}+\frac{15\pi R_{}^4J_4}{(GMr_0^5)^{\frac{1}{2}}}\frac{35\pi R_{}^6J_6}{4(GMr_0^9)^{\frac{1}{2}}}+\frac{315\pi R_{}^8J_8}{32(GMr_0^{13})^{\frac{1}{2}}}\mathrm{}=$$

$$=8.2341063\mathrm{s}1.9266\times 10^3\mathrm{s}2.34\times 10^5\mathrm{s}+5\times 10^7\mathrm{s}+\mathrm{}$$
(30)
for an orbit radius $`r_0=25498`$ km. This shows that, for such a radius, just the first four even zonal harmonics are of the same order of magnitude of $`P_l^{(\mathrm{GM})}10^7`$ s.

Let us now evaluate the impact of the even zonal harmonics the geopotential on the measurement of $`\mathrm{\Delta }P_l^{(\mathrm{GM})}`$ for a pair of counter-orbiting satellites whose orbits differ by a small amount $`dr_0`$ in radius. With

$$\frac{1}{(r_0)^{\frac{n}{2}}}\frac{1}{(r_0+d)^{\frac{n}{2}}}\frac{nd}{2(r_0)^{\frac{n+2}{2}}},$$
(31)
the difference in the perturbing terms due to the even zonal harmonics are, for $`d=5`$ km and $`r_0=25498`$ km

$`\mathrm{\Delta }P_l^{(J_2)}`$ $`=`$ $`{\displaystyle \frac{3\pi R_{}^2J_2d}{(GMr_0^3)^{\frac{1}{2}}}}=8.068\times 10^4\mathrm{s},`$ (32)
$`\mathrm{\Delta }P_l^{(J_4)}`$ $`=`$ $`{\displaystyle \frac{75\pi R_{}^4J_4d}{2(GMr_0^7)^{\frac{1}{2}}}}=9\times 10^7\mathrm{s},`$ (33)
$`\mathrm{\Delta }P_l^{(J_6)}`$ $`=`$ $`{\displaystyle \frac{315\pi R_{}^6J_6d}{8(GMr_0^{11})^{\frac{1}{2}}}}=2\times 10^8\mathrm{s},`$ (34)
$`\mathrm{\Delta }P_l^{(J_8)}`$ $`=`$ $`{\displaystyle \frac{4095\pi R_{}^8J_8d}{64(GMr_0^{15})^{\frac{1}{2}}}}=6\times 10^{10}\mathrm{s}.`$ (35)
This means that only $`\mathrm{\Delta }P_l^{(J_2)}`$ and $`\mathrm{\Delta }P_l^{(J_4)}`$ are relevant. The uncertainty in the knowledge of the Earth $`GM`$ and of the even zonal harmonics yield

$`\delta [\mathrm{\Delta }P_l^{(J_2)}]_{GM}`$ $`=`$ $`{\displaystyle \frac{3\pi R_{}^2J_2d}{2(GMr_0)^{\frac{3}{2}}}}\times \delta (GM)8\times 10^{13}\mathrm{s},`$ (36)
$`\delta [\mathrm{\Delta }P_l^{(J_2)}]_{J_2}`$ $`=`$ $`{\displaystyle \frac{3\pi R_{}^2d}{(GMr_0^3)^{\frac{1}{2}}}}\times (\delta J_2)3\times 10^{11}\mathrm{s},`$ (37)
$`\delta [\mathrm{\Delta }P_l^{(J_4)}]_{GM}`$ $`=`$ $`{\displaystyle \frac{75\pi R_{}^4J_4d}{4(GM)^{\frac{3}{2}}(r_0)^{\frac{7}{2}}}}\times \delta (GM)9\times 10^{16}\mathrm{s},`$ (38)
$`\delta [\mathrm{\Delta }P_l^{(J_4)}]_{J_4}`$ $`=`$ $`{\displaystyle \frac{75\pi R_{}^4d}{2(GMr_0^7)^{\frac{1}{2}}}}\times (\delta J_4)3\times 10^{11}\mathrm{s},`$ (39)
For $`\delta (GM)`$ and $`\delta J_{\mathrm{}}`$ the values of the IERS convention  and of the EIGEN-2 Earth gravity model , respectively, have been used.

Another source of error is represented by the uncertainty in the knowledge of the separation $`d`$ between the two orbits and of their radius $`r_0`$. Their impact is, by assuming $`\delta d\delta r_01`$ cm

$`\delta [\mathrm{\Delta }P_l^{(J_2)}]_{r_0}`$ $`=`$ $`{\displaystyle \frac{9\pi R_{}^2J_2d}{2(GM)^{\frac{1}{2}}(r_0)^{\frac{5}{2}}}}\times (\delta r_0)4\times 10^{13}\mathrm{s},`$ (40)
$`\delta [\mathrm{\Delta }P_l^{(J_2)}]_d`$ $`=`$ $`{\displaystyle \frac{3\pi R_{}^2J_2}{(GMr_0^3)^{\frac{1}{2}}}}\times (\delta d)1\times 10^9\mathrm{s},`$ (41)
$`\delta [\mathrm{\Delta }P_l^{(J_4)}]_{r_0}`$ $`=`$ $`{\displaystyle \frac{525\pi R_{}^4J_4d}{4(GM)^{\frac{1}{2}}(r_0)^{\frac{9}{2}}}}\times (\delta r_0)1\times 10^{15}\mathrm{s},`$ (42)
$`\delta [\mathrm{\Delta }P_l^{(J_4)}]_d`$ $`=`$ $`{\displaystyle \frac{75\pi R_{}^4J_4}{2(GMr_0^7)^{\frac{1}{2}}}}\times (\delta d)2\times 10^{12}\mathrm{s},`$ (43)
These results clearly show that, for the same semimajor axis of, say, the ETALON SLR satellites and for not too stringent requirements on the orbital injection errors for the separation between the two orbits, the impact of the geopotential on the measurement of $`\mathrm{\Delta }P_l^{(\mathrm{GM})}`$ is negligible.

Finally, let us consider what could be the impact of the secular variations of the even zonal harmonics $`\dot{J}_{\mathrm{}}`$. For the first even zonal harmonic the effective $`\dot{J}_2^{\mathrm{eff}}\dot{J}_2+0.371\dot{J}_4+0.079\dot{J}_6+0.006\dot{J}_80.003\dot{J}_{10}\mathrm{}`$, whose magnitude is of the order of $`(2.6\pm 0.3)\times 10^{11}`$ yr-1, can be considered. The insertion of such value in eq.(32) yields $`\mathrm{\Delta }P_l^{(\dot{J}_2^{\mathrm{eff}})}=1\times 10^{11}`$ s over one year, so that it can be concluded that its effect can be safely neglected.

#### 3.3.2 The Tides

Another source of potential systematic error is represented by the solid Earth and ocean tides .

For a constituent of given frequency $`\mathrm{f}`$ of the solid Earth tidal spectrum, which is the more effective in perturbing the Earth satellites’orbits, the perturbation of degree $`\mathrm{}`$ and order $`m`$ induced on $`l`$ is

$$\frac{dl}{dt}=n+A_\mathrm{}m(GM)^{\frac{1}{2}}R_{}^\mathrm{}1r_0^{(\frac{2\mathrm{}+3}{2})}2(\mathrm{}+1)\underset{p=0}{\overset{\mathrm{}}{}}\underset{q=\mathrm{}}{\overset{+\mathrm{}}{}}F_{\mathrm{}mp}G_{\mathrm{}pq}k_\mathrm{}m^{(0)}H_{\mathrm{}}^m\mathrm{cos}\gamma _{\mathrm{f}\mathrm{}mpq},$$
(44)
where $`A_\mathrm{}m=[(2\mathrm{}+1)(\mathrm{}m)!/(4\pi )(\mathrm{}+m)!]^{\frac{1}{2}}`$, $`F_{\mathrm{}mp}(i)`$ and $`G_{\mathrm{}pq}(e)`$ are the inclination and eccentricity functions, respectively , $`k_\mathrm{}m^{(0)}`$ are the Love numbers, $`H_{\mathrm{}}^m`$ are the tidal heights and $`\gamma _{\mathrm{f}\mathrm{}mpq}`$ is built up with the satellite’s orbital elements $`\mathrm{\Omega }`$, $`\omega `$ and $``$, the lunisolar variables and the lag angle $`\delta _{\mathrm{f}\mathrm{}m}`$ . Note that the dependence on $`F_{\mathrm{}mp}`$ and $`G_{\mathrm{}pq}`$ is the same also for the ocean tidal perturbations.

From an inspection of the explicit expressions of the inclination and eccentricity functions  and from the condition $`\mathrm{}2p+q=0`$, which must be fulfilled for the perturbations averaged over an orbital revolution, it turns out that the $`\mathrm{}=3`$ part of the entire tidal spectrum does not affect the mean longitude of a satellite with $`i=e=0`$. With regard to the $`\mathrm{}=2`$ constituents, only the long-period zonal ($`m=0`$) tides induce non vanishing perturbations on $`P_l`$ for circular and equatorial orbits. Among them, the most powerful is by far the 18.6-year lunar tide. Its effect on $`P_l`$ is given by

$$P_l^{(18.6\mathrm{yr})}=\frac{3(5\pi )^{\frac{1}{2}}R_{}k_{20}^{(0)}(\overline{\mathrm{f}})H_2^0(\overline{\mathrm{f}})}{(GMr_0)^{\frac{1}{2}}}\mathrm{cos}\gamma _{\overline{\mathrm{f}}2010}.$$
(45)
For $`r_0=25498`$ km it amounts to 3$`\times 10^6`$ s. The difference in the orbits of the two satellites would yield

$$\mathrm{\Delta }P_l^{(18.6\mathrm{yr})}=\frac{3(5\pi )^{\frac{1}{2}}R_{}k_{20}^{(0)}(\overline{\mathrm{f}})H_2^0(\overline{\mathrm{f}})d}{2(GMr_0^3)^{\frac{1}{2}}}\mathrm{cos}\gamma _{\overline{\mathrm{f}}2010}6\times 10^{10}\mathrm{s}$$
(46)
for $`d=5`$ km and $`r_0=25498`$ km. As can be seen, $`\mathrm{\Delta }P_l^{(\mathrm{GM})}\mathrm{\Delta }P_l^{(18.6\mathrm{yr})}`$, so that the impact of the direct tidal perturbations on the measurement of the gravitomagnetic clock effect on $`l`$ can be neglected.

The cross-coupling between the static Earth oblateness and the Earth tides  would induce the following perturbations on the mean longitude

$$\frac{dl}{dt}=n\frac{21J_2GMR_{}^{\mathrm{}+1}A_\mathrm{}m_{p=0}^{\mathrm{}}_{q=\mathrm{}}^+\mathrm{}F_{\mathrm{}mp}G_{\mathrm{}pq}[(\mathrm{}2p)m]k_\mathrm{}m^{(0)}H_{\mathrm{}}^m}{2r_0^{\mathrm{}+5}\mathrm{f}}\mathrm{cos}\gamma _{\mathrm{f}\mathrm{}mpq}.$$
(47)
From an inspection of the explicit expressions of the inclination functions (see Appendix A) and from the condition $`\mathrm{}2p=0`$, which holds for even $`\mathrm{}`$, it turns out that the $`\mathrm{}=2`$ part of the indirect tidal spectrum, which is the most powerful, has no effect on $`l`$ for an equatorial orbit. More precisely, $`F_{201}`$ does not vanish for $`i=0^{}`$, but in this case, $`\mathrm{}2p=m=0`$. The other $`\mathrm{}=2`$ inclination functions vanish for equatorial orbits. The same holds also for the $`\mathrm{}=3`$ constituents for which $`\mathrm{}2p0`$. It turns out that $`F_{311}0`$ for equatorial orbits; in this case $`\mathrm{}2p=m=1`$. The other $`\mathrm{}=3`$ inclination functions vanish for $`i=0^{}`$.

#### 3.3.3 The Impact of the Errors in the Inclinations

Until now we have propagated the uncertainties in the orbit radius by assuming for the inclination and the eccentricity their nominal values $`i=e=0`$. Now we wish to fix $`r_0`$ and propagate the errors in the inclination . As pointed out before, the uncertainties in $`i`$ affect the classical perturbative corrections to the Keplerian periods. The effects of the Earth oblateness are the most prominent; thus, we will investigate them in order to establish upper bounds in the admissible separation $`I`$ between the two orbital planes.

The radial component of the perturbing acceleration due to $`J_2`$ for a generic orbit reads

$$A_R^{(J_2)}=\frac{3GMR_{}^2J_2}{4r^4}\left[3\mathrm{cos}2(\omega +f)+3\mathrm{cos}^2i3\mathrm{cos}^2i\mathrm{cos}2(\omega +f)1\right].$$
(48)
Neglecting all terms of order $`𝒪(e^2)`$, eq.(48) in eq.(2) yields

$$P_l^{(J_2)}=\frac{3\pi R_{}^2J_2}{(GMr_0)^{\frac{1}{2}}}(13\mathrm{cos}^2i).$$
(49)
Then, for the counter-orbiting satellites along identical (almost) circular orbits with inclinations $`i^{(+)}i`$ and $`i^{()}=(180^{}i)+I`$, with $`I/i1`$, from $`\mathrm{cos}^2i1+i^2`$,

$$\mathrm{\Delta }P_l^{(J_2)}\frac{18\pi R_{}^2J_2}{(GMr_0)^{\frac{1}{2}}}(iI).$$
(50)
For $`r_0=25498`$ km and $`i=0.01^{}`$, $`\mathrm{\Delta }P_l^{(J_2)}\mathrm{\Delta }P_l^{(\mathrm{GM})}`$ if $`I0.006^{}=1\times 10^4`$ rad. It is worth noting that the current technology does allow to obtain equatorial orbits tilted by $`0.01^{}`$ to the equator, e.g., for many geostationary satellites. Many orbital perturbations affect the inclination with shifts $`\mathrm{\Delta }i`$ which must be kept smaller than $`I`$. The static part of geopotential does not induce long-term perturbations on the inclinations. The tesseral and sectorial solid Earth and ocean tides induces long-period perturbations on $`i`$ of the order of a few tens of mas. The general relativistic gravitoelectric de Sitter, or geodetic, precession  induces a secular variation of the inclination of 84 mas yr-1. The inclination is sensitive to the non-gravitational perturbations, but the drag-free apparatus could resolve this problem. Suffices it to say that the uncancelled non-conservative perturbations on $`i`$ for the LAGEOS satellites amount to a few tens of mas yr-1 . It is also important to note that departures of $`i`$ from their nominal values of the order of $`I0.0001^{}`$ are feasible with the current technologies for, e.g., the planned GP-B mission  and the current GRACE mission121212 http://www.csr.utexas.edu/grace/newsletter/2002/august2002.html. It turns out that the requirements on the semimajor axis and the eccentricity are far less demanding than those on the inclination .

#### 3.3.4 The N-Body Gravitational Perturbations

Another source of perturbations on an Earth satellite’s mean longitude is represented by the gravitational effect induced by the major bodies of the Solar System (Sun, Moon, Jupiter, other planets, the asteroids). Let us calculate the effects induced by some of them.

The perturbative effect of the planet of mass $`m^{^{}}`$ on the satellite of mass $`m`$ is given by131313Here we adopt the Lagrangian approach; $``$ is the disturbing function and the (approximate) equation for the rate of the mean longitude is
$$\frac{dl}{dt}=n\frac{2}{na}\frac{}{a}.$$
(51)

$$_{\mathrm{planets}}=Gm^{^{}}\left(\frac{1}{|𝐫𝐫^{^{}}|}\frac{𝐫𝐫^{^{}}}{r^{}_{}{}^{}3}\right).$$
(52)
It turns out that the second term in eq.(52) does not induces secular perturbations.

After expressing the first term of eq.(52) in terms of the orbital elements of the satellite and the planet and averaging it over one period of the mean longitudes $`l`$ and $`l^{^{}}`$ one finds for the largest contribution which does not contain the terms of second order in the eccentricities and the inclinations

$$P_l^{(\mathrm{planets})}=\frac{4\pi Gm^{^{}}}{n^3a^{}_{}{}^{}3}.$$
(53)
It is important to notice that eq.(53) does not depend on the sense of motion along the orbits, so that it cancels, in principle, in $`\mathrm{\Delta }P_l`$. Instead, the difference $`d`$ in the semimajor axes of the two satellites would induce an aliasing effect

$$\mathrm{\Delta }P_l^{(\mathrm{planets})}\frac{18\pi m^{^{}}da^{\frac{7}{2}}}{G^{\frac{1}{2}}a^{}_{}{}^{}3M_{}^{\frac{3}{2}}}.$$
(54)
The nominal values of eq.(54) for the Sun and the Moon amount to 1.178$`\times 10^4`$ s and $`2.565\times 10^4`$ s, respectively. However, the uncertainties in $`Gm^{^{}}`$ yield

$`\delta [\mathrm{\Delta }P_l^{(\mathrm{Sun})}]_{GM_{}}`$ $`=`$ $`7\times 10^{15}\mathrm{s},`$ (55)
$`\delta [\mathrm{\Delta }P_l^{(\mathrm{Moon})}]_{Gm_{\mathrm{Moon}}}`$ $`=`$ $`6\times 10^{11}\mathrm{s},`$ (56)
where we have used $`\delta (GM)_{}=8\times 10^9`$ m3 s-2  and $`\delta (Gm)_{\mathrm{Moon}}=1.2\times 10^6`$ m3 s-2 . The errors induced by the uncertainties in $`d`$, $`a`$ and $`GM_{}`$ turn out to be negligible. The corrections $`\mathrm{\Delta }P_l^{(\mathrm{planets})}`$ for Venus, Mars, Jupiter and Ceres, whose mass amounts to approximately one third of the total mass of asteroids $`m^{\mathrm{asteroids}}=2.3\times 10^{24}`$ g, are negligible because their nominal values are $`10^9`$ s. This is particularly important for Jupiter, whose contribution would amount to $`\mathrm{\Delta }P_l^{(\mathrm{Jup})}=1\times 10^9`$ s, because its sidereal orbital period, with respect to the Sun, amount to almost 11 years.

It can be shown that the indirect effect induced on $`P^{(0)}`$ by the perturbations on $`a`$ can be neglected. Indeed, there are no N-body secular perturbations on $`a`$ and, most importantly141414After all, short-periodic terms in the planetary longitudes would result in semi-secular signatures on the time-scales of an Earth satellite., their amplitudes are proportional to factors like $`e^{|h_1|}e^{{}_{}{}^{}|h_2|}\mathrm{sin}^{|h_3|}i\mathrm{sin}^{|h_4|}i^{^{}}`$, where $`h_1,h_2,h_3`$ and $`h_4`$ are integer numbers  constrained by the conditions $`|h_1|+|h_2|+|h_3|+|h_4|2`$ and, for the short–period perturbations, $`h_1+h_2+h_3+h_40`$.

### 3.4 The Impact of the Non-Gravitational Perturbations

From eq.(2) it turns out that, in general, periodically time-varying radial accelerations with frequencies multiple of the orbital frequency do not affect the mean longitude over an orbital revolution, at least at order $`𝒪(e)`$.

The situation is different for radial accelerations which can be considered constant in time over an orbital revolution. Let us consider an acceleration of non–gravitational origin whose radial component $`A_R`$ is constant over the timescale of one orbital revolution. From eq.(2) it turns out that

$$P_l^{(\mathrm{NG})}=4\pi \left(\frac{r_0^7}{GM^3}\right)^{\frac{1}{2}}A_R,$$
(57)
from which it follows

$$\mathrm{\Delta }P_l^{(\mathrm{NG})}=14\pi d\left(\frac{r_0^5}{GM^3}\right)^{\frac{1}{2}}A_R.$$
(58)
The maximum value of $`A_R`$ which makes $`\mathrm{\Delta }P_l^{(\mathrm{NG})}\mathrm{\Delta }P_l^{(\mathrm{GM})}`$ is $`A_R^{\mathrm{max}}=6\times 10^7`$ cm s-2 for $`d=5`$ km and $`r_0=25498`$ km. It is not a too stringent constraint. For LAGEOS, which is completely passive, the largest non–gravitational acceleration, i.e. the direct solar radiation pressure, amounts to almost $`4\times 10^7`$ cm s-2; a drag-free apparatus with not too stringent performances–well far from the $`10^{13}`$ cm s-2 Hz$`^{\frac{1}{2}}`$ level of LISA–could meet in a relatively easy way such requirement. For $`\delta d\delta r_01`$ cm, $`\delta (GM)=8\times 10^{11}`$ cm3 s-2 and $`A_R=6\times 10^7`$ cm s-2 it turns out

$`\delta [\mathrm{\Delta }P_l^{(\mathrm{NG})}]_d`$ $`=`$ $`14\pi \left({\displaystyle \frac{r_0^5}{GM^3}}\right)^{\frac{1}{2}}A_R\times (\delta d)=1\times 10^{12}\mathrm{s},`$ (59)
$`\delta [\mathrm{\Delta }P_l^{(\mathrm{NG})}]_{r_0}`$ $`=`$ $`35\pi d\left({\displaystyle \frac{r_0}{GM}}\right)^{\frac{3}{2}}A_R\times (\delta r_0)=5\times 10^{16}\mathrm{s},`$ (60)
$`\delta [\mathrm{\Delta }P_l^{(\mathrm{NG})}]_{GM}`$ $`=`$ $`21\pi d\left({\displaystyle \frac{r_0}{GM}}\right)^{\frac{5}{2}}A_R\times \delta (GM)=1\times 10^{15}\mathrm{s}.`$ (61)
These results show that the direct impact of the non-gravitational accelerations on the measurement of the gravitomagnetic time shift on $`l`$ could be made negligible with a drag-free apparatus of relatively modest performances, contrary to the indirect effects on the difference of the Keplerian periods which, as seen before, would require a much more effective drag–free cancellation in the frequency range $`[0,10^5\mathrm{Hz}]`$.

## 4 Conclusions

In this paper we have examined the possibility of measuring the gravitomagnetic clock effect on the mean longitudes of a pair of counter-rotating satellites following almost identical circular equatorial orbits in the gravitational field of Earth. While the gravitomagnetic signature depends only on the Earth parameters, the aliasing classical effects depend also on the orbital geometry of the satellites, so that it is possible to choose it suitably in order to reduce their impact on the measurement of the post-Newtonian effect. In this respect, our choice of a nominal value of 25498 km for the semimajor axis yields good results.

From the point of view of the observational accuracy, one source of error would be the uncertainty with which the axes of the ITRS are known. However, it turns out that the level of the gravitomagnetic effect can be reached over a sufficiently high number of orbital revolutions.

In regard to the systematic errors, the major limiting factors are

* The unavoidable imperfect cancellation of the Keplerian periods, which yields a 10-2 cm constraint in knowing the difference $`d`$ between the semimajor axes of the satellites
* The required accuracy in the knowledge of the inclinations $`i`$ of the satellites in presence of not exactly equatorial orbits. For $`i=0.01^{}`$ the difference between the inclinations of the two orbital planes $`Ii^{(+)}i^{()}`$ should be less than 0.006.

A pair of drag-free ($`10^9`$ cm s-2 Hz$`^{\frac{1}{2}}`$) spacecraft endowed with an intersatellite tracking apparatus might allow to meet the stringent requirements of such a mission whose realization seems to be very difficult although, perhaps, not completely impossible with the present-day or forthcoming space technologies.

## Acknowledgments

L.I. is grateful to L. Guerriero for his support while at Bari. L.I. and H. L. gratefully thanks the 2nd referee for his/her useful and important efforts for improving the paper.


format: text

text:
# Cluster Transformation Coefficients for Structure and Dynamics Calculations in n-Particle Systems: Atoms, Nuclei, and Quarks

## 1 Introduction

In the Heisenberg’s picture the time evolution of a system of particles is described by a commutator equation which involves the n-body Hamiltonian and the creation operators of the n-body ground- and excited-modes. However, the excitations of the reference vacuum, resulting from the scattering of particles from the vacuum to higher states (the particle-hole and particle-particle excitations), are completely neglected in this formulation, namely, the time evolution of the n-body-modes is described by a linearized Equation of Motion (EOM)  which involves only valence particles. Recently the original Heisenberg formulation of the n-body-dynamics has been generalized within the framework of the Dynamic-Correlation Model (DCM)  for fermions, the Boson Dynamic-Correlation Model (BDCM)  for bosons, and the Superconductive Dynamic- Correlation Model (SDCM) which describes superconductive- and polarization-effects . In these models, the inclusion of the structure and the dynamics of odd/even quasi-particles into the calculations of the excitations of the model-vacuum led to the modification of the formal strcture of the original Heisenberg’s commutator equation. The new dynamics system is characterized by a set of coupled commutator equations which involve simultaneously the excitations of the valence quasi-particles and the excitations of the Intrinsic-Supercunductive-Vacuum States (ISVSs) i.e.: mixed-mode states formed by coupling the valence quasi-particles to the excitations of the vacuum (particle-particle and particle-hole). As in Ref. , the resulting mixed-mode states are classified in terms of the following mixed-mode wave functions: (a) valence particles coupled to ISVSs formed by include particle-hole, particle-particle and hole-hole vacuum-excitation modes that have the same parity of the valence particles; (b) valence particles and ISVSs formed by exciting particle-hole, particle-particle and hole-hole having a parity different from that of the valence particles. The superconductive vacuum states (b) are not considered by perturbative theories, although they may be associated to the creation of virtual particles  giving important contribution to dynamic theories. The superconductive vacuum excitation modes are important mainly at high densities , although the overlap of the particle-particle excitation modes may influence, due to the strong Pauli blocking effects in the ISVSs, the particle-hole excitation mechanism also at medium energies. In this paper we disregard the superconductive effects, which will be the subject of a paper in preparation, and discuss only the effect of the particle-hole excitation mechanism on ground state properties of atomic, nuclear, and quark systems. In principle, the system of commutator equations may be solved perturbatively by means of substitution method which consists in inserting the higher order commutators into the previous calculated commutators. However, such perturbative method is usually not convergent for systems of strong interacting particles and, therefore, will not be discussed here. We will study nonperturbative solutions based on the following: (a) a linearization ansatz motivated by the consideration that in the low energy domain only few particle excitations are energetically representative while the higher excitation modes (with admixture probability less than 0.1 %) may be omitted from the model space; (b) cluster transformation coefficients which provide an exact solution for the complex n-body matrix elements that are input of the commutator equations. The wave functions, solutions of the calculated collective eigenvalue equation are classified as in Refs.  in terms of Configuration Mixing Wave Functions (CMWFs). In Section II, we calculate the commutator chain for one particle-, two particle-, and three particle- creation operators. The generalization of the model to hole- and (particle-hole)-valence operators is not given explicitly because of easy extrapolation from the present particle formulation. The commutator chain is then linearized and converted into an eigenvalue problem which is solved calculating the n-body matrix elements within the cluster transformation coefficients. In Section III, solutions of the nonlinear eigenvalue equation for one valence particle/hole are obtained and the HFS constants of hydrogen- and lithium-like heavy ions are calculated. In Section IV, the commutator chain is solved for two and three valence particles and the resulting CMWFs are used to calculate the matter distribution of halo nuclei. In Section V, an application to the three dressed quark systems is discussed in order to provide new theoretical data for the polarizability of the proton. For a measurement of the proton polarization, novel ultra/intense lasers may be useful, like the PHELIX-petawatt laser presently under construction at GSI, have to be used. A direct way for a precision measurement seems to be in reach by the latest improvements of energetic multi-MeV photon sources using laser backscattering at electron storage rings.

## 2 The nonlinear commutator chain

The nonlinear commutator chain is an extension of the Heisenberg’s equation. It allows us to address the situation where valence clusters and core clusters are almost energetically degenerated and may, therefore, coexist. We introduce valence systems formed of either neutrons or protons. and define the valence states as

$`|\mathrm{\Gamma }_J(\{n\}particles)={\displaystyle \underset{\alpha _n}{}}X_{\alpha _nJ_n}^{(n)}N_{\alpha _nJ_n}^{(n)}A_n^{}(\alpha _nJ_n;J)|0,`$ (1)
where $`n=1,2,3`$ and $`J`$ is the total spin. The $`N_{\alpha _nJ_n}^{(n)}`$ are the normalization constants, $`X_{\alpha _nJ_n}^{(n)}`$ the mode amplitudes, and $`\alpha _1`$ denotes the other quantum numbers of the valence particles. The creation operator is defined by

$`A_1^{}(\alpha _1j_1;J)a_{j_1}^{}`$ (2)
for one valence particle (fermion);

$`A_2^{}(\alpha _2J;J)(a_{j_1}^{}a_{j_2}^{})^J`$ (3)
for one valence pair (boson);

$`A_3^{}(\alpha _3J_1;J)((a_{j_1}^{}a_{j_2}^{})^{J_1}a_{j_3}^{})^J`$ (4)
for three valence particles (fermion). Hence, in Eq. (1), $`J_n=j_1`$ for one particle, $`J_n=J`$ for two particles, and $`J_n=J_1`$ for three particles.
In the literature, various approaches, as reported in Ref. , have been proposed for including core excitations. In Shell-Model calculations the residual interaction between the valence particles excites the valence pairs to higher single particle states, leaving the vacuum in the $`J=0`$ ground state. In Ref.  the core excitations have been included in the Shell-Model calculations through introducing the coupling between the valence particles and the collective $`J=2^+`$ core-excited states. In the DCM and BDCM the core excitation is included through coupling the valence fermionic/bosonic states, Eqs. (23) to intrinsic bosonic states corresponding to particle-hole excitations of the nuclear core. In this paper, we consider the following mixed-mode fermionic/bosonic states: a) valence fermionic/bosonic states coupled to the dynamic particle-hole states of normal parity; b) valence fermionic/bosonic states coupled to the dynamic particle-hole states of non-normal parity. The particle-hole coupling is implemented through the two-body force $`V`$, such that

$`H={\displaystyle \underset{\alpha }{}}ϵ_\alpha a_\alpha ^{}a_\alpha +{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha \beta \gamma \delta }{}}v_{\alpha \beta \gamma \delta }a_\alpha ^{}a_\beta ^{}a_\delta a_\gamma =H_0+V,`$ (5)
where $`v_{\alpha \beta \gamma \delta }`$ are the matrix elements of the realistic two-body potential which has two parts: a central part and a tensor part. The tensor component of the realistic two-body potential shapes the many-body Hamiltonian and in particular its long tail acting between the valence fermions/bosons and causing simultaneously the excitation of the particles to high shell-model states and the deformation of the nuclear core.

The mixed-mode $`\{123\}particles\{1^{}2^{}\}bosons`$ states, denoted $`\mathrm{\Gamma }_J^{},`$ can therefore be expanded as follows:

$`|\mathrm{\Gamma }_J^{}(\{n=123particles\}\{1^{}2^{}\}bosons)=[{\displaystyle \underset{\alpha _nJ_n}{}}X_{\alpha _nJ_n}^{(n)}N_{\alpha _n;J_n}^{(n)}A_n^{}(\alpha _nJ_n;J)`$
$`+{\displaystyle \underset{\alpha _{n+1^{}}J_nJ_{n+1^{}}}{}}X_{\alpha _{n+1^{}}J_nJ_{n+1^{}}}^{(n+1^{})}N_{\alpha _{n+1^{}}J_nJ_{n+1^{}}}^{(n+1^{})}A_{n+1^{}}^{}(\alpha _{n+1^{}}(J_nJ_{n+1^{}};J))`$
$`+{\displaystyle \underset{\alpha _{n+2^{}}J_nJ_{n+1^{}}J_{n+2^{}}}{}}X_{\alpha _{n+2^{}}J_nJ_{n+1^{}}J_{n+2^{}}}^{(n+2^{})}N_{\alpha _{n+2^{}}J_nJ_{n+1^{}}J_{n+2^{}}}^{(n+2^{})}A_{n+2^{}}^{}(\alpha _{n+2^{}}(J_nJ_{n+1^{}}J_{n+2^{}};J))]|0,`$ (6)
where $`1,2,3`$ labels the particles and $`1^{},2^{}`$ label the bosons (particle-hole pairs). The operators

$`A_{n+1^{}}^{}(\alpha _{n+1^{}}(J_1J_{n+1^{}});J)([A_n^{}(\alpha _n(J_n;J)(a_{j_{n+1}}^{}b_{j_{n+2}}^{})^{J_{n+1^{}}}]^{J_{(n,n+1^{})}})^J,`$ (7)
$`A_{n+2^{}}^{}(\alpha _{n+2^{}}(J_nJ_{n+1^{}}J_{n+2^{}});J)([A_n^{}(\alpha _n(J_n;J))`$
$`((a_{j_{2n+1}}^{}b_{j_{2n+2}}^{})^{J_{n+1^{}}}(a_{j_{2n+3}}^{}b_{j_{2n+4}}^{})^{J_{n+2^{}}})^{J_{(n+1^{},n+2^{})}}]^{J_n,J_{(n+1^{},n+2^{})}})^J`$ (8)
contain consequently the hole creators $`b_j^{}`$. To conclude, in the DCM and BDCM one starts with the valence fermions/bosonic states of Eqs. (24) and constructs subsequently mixed-mode nuclear states by including components having additional bosons formed by the particle-hole pairs of the core excitations. The resulting nuclear states are then classified in terms of (CMWFs) of increasing degrees of complexity (number of particle-hole pairs) Ref. . The basic dynamic equations of the models are the commutator equations that involve the nuclear many-body Hamiltonian $`H`$ and the operators $`A_n^{}(n=1,2,3)`$. After a lengthy but elementary algebra, we obtain the following results:

(a) Commutator equation for $`\{n=1,2,3\}`$ states:

$`[H,A_n^{}(\alpha _nJ_n;J)]={\displaystyle \underset{\beta _nJ_n^{}}{}}A_n(\alpha _nJ_n;J)HA_n^{}(\beta _nJ_n^{};J)A_n^{}(\beta _nJ_n^{};J)`$
$`+{\displaystyle \underset{\beta _{n+1^{}}J_n^{}J_{n+1^{}}^{}}{}}A_n(\alpha _nJ_n;J)HA_{n+1^{}}^{}(\beta _{n+1^{}}J_n^{}J_{n+1^{}}^{};J)`$
$`\times A_{n+1^{}}^{}(\beta _{n+1}J_n^{}J_{n+1^{}}^{};J).`$ (9)
(b) and Commutator equation for $`\{n=1,2,3;n^{}=1^{}\}`$ states

$`[H,A_{n+1^{}}^{}(\alpha _{n+1^{}}J_nJ_{n+1^{}};J)]`$
$`={\displaystyle \underset{\beta _{n+1^{}}J_n^{}J_{n+1^{}}^{}}{}}A_{n+1^{}}(\alpha _{n+1^{}}J_n^{}J_{n+1^{}};J)HA_{n+1^{}}^{}(\beta _{n+1^{}}J_n^{}J_{n+1^{}}^{};J)`$
$`\times A_{n+1^{}}^{}(\beta _{n+1^{}}J_n^{}J_{n+1^{}}^{};J)`$
$`+{\displaystyle \underset{\beta _{n+2^{}}J_n^{}J_{n+1^{}}^{}J_{n+2^{}}^{})}{}}A_{n+1^{}}(\alpha _{n+1^{}}J_nJ_{n+1^{}};J)HA_{n+2^{}}^{}(\beta _{n+2^{}}J_n^{}J_{n+1^{}}^{}J_{n+2^{}}^{});J)`$
$`\times A_{n+2^{}}^{}(\beta _{n+2^{}}J_n^{}J_{n+1^{}}^{}J_{n+2^{}}^{};J).`$ (10)
The commutator equations are used to obtain the solutions of DCM with use of the Equations of Motion (EOM) method . In this latter method, one looks for an operator $`C_m^{}`$ such that $`C_m^{}|0=|m`$ with $`H|m=E_m|m`$ and $`H|0=E_0|0`$ $`(E_00)`$. Here $`|m`$ denotes the excited state and $`|0`$ the correlated ground state. One then has the operator identity $`[H,C_m^{}]=E_mC_m^{}`$. Upon replacing $`C^{}`$ by $`A^{}`$ and using the results in Eqs. (910) for the l.h.s. of the above operator identity, one obtains a set of equations which, after linearization, can be transformed into a system of eigenvalue equations. As one can see, the linearization consists in approximating the higher-order diagrams by an effective term. The solutions of the linearized commutator equations, can, therefore, be regarded as eigenvalues of a model fermionic/bosonic Hamiltonian.

The valence particles become, therefore, the dressed solutions of the following collective nonlinear Hamiltonian :

$`H_{coll}=(T_{coll}+V_{coll})={\displaystyle \underset{i}{\overset{N}{}}}E_{coll}^i𝒜_i^{}(\alpha _iJ_1\mathrm{}J_i;J)𝒜_i(\alpha _iJ_1\mathrm{}J_i;J)`$ (11)
with

$`𝒜_i^{}(\alpha _i(J_1\mathrm{}J_i);J)=[X_{\alpha _iJ_i\mathrm{}J_i}^iA_i^{}(\alpha _iJ_i\mathrm{}J_i;J)`$
$`+X_{\alpha _{i+1^{}}J_i\mathrm{}J_{i+1^{}}}A_{i+1^{}}^{}(\alpha _{i+1^{}}J_1\mathrm{}J_{i+1^{}};J)+X_{\alpha _{i+2^{}}J_i\mathrm{}J_{i+2^{}}}A_{i+2^{}}^{}(\alpha _{i+2^{}}J_i\mathrm{}J_{i+2^{}};J)]`$ (12)
where the $`X_\alpha `$ are the calculated mixed-mode amplitudes.

The input to solving the commutator equations are the matrix elements of the $`nparticle`$ configuration mixing wave functions (CMWFs). The latter can be easily calculated by use of the cluster-factorization method of Ref. . The method is guided by the observation that the EOM connect the $`n`$particle states to the $`(n1)`$-particle states. In the DCM, the $`nn1`$ reduction is achieved by factorizing out 1 boson, which can be $`\{2p\}`$, $`\{1p1h\}`$, or $`\{2h\}`$. In the following, we exemplify the method by considering the $`\{4p2h\}`$ $`(n=3)`$ parent configuration. For the sake of simpler notation, we will not write the detailed coupling of the creation operators but write the wave functions $`\mathrm{\Psi }`$. There are three active pairs associated with the $`\mathrm{\Phi }_{JM}^3(\alpha _3J_1J_2J_3)`$ configuration. To go to the $`(n=2)`$ (two active pairs) configuration, we can factorize out either $`\{1p1h\}`$, $`\{2p\}`$, or $`\{2h\}`$ to arrive at $`\{3p1h\}`$, $`\{2p2h\}`$, or $`\{4p\}`$, respectively. This leads to the following expansion in terms of $`\{3p1h\}`$-CMWFs with coordinates $`\{\alpha _2\}`$, $`\{2p2h\}`$-CMWFs with coordinates $`\{\lambda _2\}`$, and $`\{4p\}`$-CMWFs with coordinates $`\{ϵ_2\}`$:

$`|\mathrm{\Phi }_J^{4p2h}(\alpha _3J_1J_2J_3)=`$
$`={\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \underset{\alpha _2\overline{\alpha }_2J_{k_1}J_{k_2}J_rJ_s}{\overset{3}{}}}T_J(\alpha _3J_1J_2J_3|\}\alpha _2J_r\overline{\alpha }_2J_s)[|\mathrm{\Psi }_{J_r}^{3p1h}(\alpha _2J_{k_1}J_{k_2})|\mathrm{\Psi }_{J_s}^{1p1h}(\overline{\alpha }_2J_s)]^J`$
$`+{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \underset{\lambda _2\overline{\lambda }_2J_{k_1}J_{k_2}J_rJ_s}{\overset{3}{}}}V_J(\alpha _3J_1J_2J_3|\}\lambda _2J_r\overline{\lambda }_2J_s)[|\mathrm{\Psi }_{J_r}^{2p2h}(\lambda _2J_{k_1}J_{k_2})|\mathrm{\Psi }_{J_s}^{2p}(\overline{\lambda }_2J_s)]^J`$
$`+{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \underset{ϵ_2\overline{ϵ}_2J_{k_1}J_{k_2}J_rJ_s}{\overset{3}{}}}Z_J(\alpha _3J_1J_2J_3|\}ϵ_2J_r\overline{ϵ}_2J_s)[|\mathrm{\Psi }_{J_r}^{4p}(ϵ_2J_{k_1}J_{k_2})|\mathrm{\Psi }_{J_s}^{2h}(\overline{ϵ}_2J_s)]^J.`$ (13)
The cluster transformation coefficients in Eq. (13), i.e: the $`{}_{}{}^{3}T_{J}^{}(\alpha _3J_1J_2J_3|\}\alpha _2J_r\overline{\alpha }_2J_s)`$, the $`{}_{}{}^{3}V_{J}^{}(\alpha _3J_1J_2J_3|\}\lambda _2J_r\overline{\lambda }_2J_s)`$, and the $`{}_{}{}^{3}Z_{J}^{}(\alpha _3J_1J_2J_3|\}ϵ_2J_r\overline{ϵ}_1J_s)`$ are then calculated, as in Ref. , by reducing the $`SU_{2J+1}(3)`$ representations carried by the wave functions on the right-hand side of Eq. (7), i.e., by diagonalizing the Casimir operators  of the $`SU_{2J+1}(3)`$ groups in the basis states of Eq. (13). Hence, the coefficients $`{}_{}{}^{3}T_{J}^{}(\alpha _3J_1J_2J_3|\}\alpha _2J_r\overline{\alpha }_2J_s)`$ are eigenvalues of the following matrix:

$`({}_{}{}^{3}T_{J}^{}(\alpha _3J_1J_2J_3|\}\alpha _2J_r\overline{\alpha }_2J_s))^{}{}_{}{}^{3}T_{J}^{}(\alpha _3J_1J_2J_3|\}\alpha _2J_r\overline{\alpha }_2J_s)=`$
$`={\displaystyle \underset{kJ_iJ_i^{}J_r^2J_r^3}{}}(1)^{J_i+J_i^{}+J_r^2+J_r^3+J_s+J_s^{}}(\widehat{k})^{1/2}\widehat{J}_r\widehat{J}_r^{}`$
$`\left\{\begin{array}{ccc}J_r& J& J_s\\ J& J_r& k\end{array}\right\}\left\{\begin{array}{ccc}J_r^{}& J& J_s^{}\\ J& J_r^{}& k\end{array}\right\}\left\{\begin{array}{ccc}J_i& J_r& J_r^2\\ J_r& J_i& k\end{array}\right\}\left\{\begin{array}{ccc}J_i^{}& J_r^{}& J_r^3\\ J_r^{}& J_i^{}& k\end{array}\right\}`$ (22)
$`({}_{}{}^{2}T_{J_r}^{k}(\alpha _2J_{k_1}J_{k_2}|\}\alpha _1J_i\overline{\alpha }_1J_r^2))^{}{}_{}{}^{2}T_{J_r}^{k}(\alpha _2J_{k_1}J_{k_2}|\}\alpha _1J_i\overline{\alpha }_1J_r^2)`$
$`({}_{}{}^{2}T_{J_r^{}}^{k}(\beta _2J_{k_3}J_{k_4}|\}\beta _1J_i^{}\overline{\beta }_1J_r^3))^{}{}_{}{}^{2}T_{J_r^{}}^{k}(\beta _2J_{k_3}J_{k_4}|\}\beta _1J_i^{}\overline{\beta }_1J_r^3).`$ (23)
In Eq. (23), the sum over $`J_{k_1}`$, $`J_{k_2}`$, $`J_{k_3}`$, $`J_{k_4}`$, $`\alpha _1`$, $`\overline{\alpha _1}`$, $`\beta _1`$, and $`\overline{\beta _1}`$ has not been written explicitly but is understood. The 6$`J`$ coefficients are defined according to Ref. . The same procedure has been applied to calculate the $`{}_{}{}^{3}V`$ and $`{}_{}{}^{3}Z`$ cluster-transformation coefficients. The cluster-transformation coefficients of the $`{}_{}{}^{4}T`$ are then expanded in terms of $`{}_{}{}^{3}T`$ using Eq. (23) in the limit $`{}_{}{}^{4}T^3T`$ and $`{}_{}{}^{3}T^2T`$. In tables 1 to 5 the transformation coefficients are exemplified. The above recursive procedure is equally valid for cluster transformation coefficients involving fermionic CMWFs. Within the cluster transformation coefficients introduced for the fermionic/bosonic CMWFs, the matrix elements needed for the EOM can be readily obtained:

The matrix elements for two particles are

$`(a)\mathrm{\Phi }_J^{2p}|V|\mathrm{\Phi }_J^{2p}=\mathrm{\Phi }_J^{2p}|G(E)|\mathrm{\Phi }_J^{2p}`$
$`(b)\mathrm{\Psi }_J^{1p1h}|V|\mathrm{\Psi }_J^{1p1h}=\mathrm{\Psi }_J^{1p1h}|V_{phen}|\mathrm{\Psi }_J^{1p1h}`$ (24)
In Eq. (24), G(E) is the realistic G-matrix, and $`V_{phen}`$ is a phenomenological potential. The matrix elements for:

$`(c)\{3p\},\{2p1h\},\{2h1h\},\{3p1h\},\{4p\},\{4p1h\},\{4p2h\}`$ (25)
are given symbolically in Fig. 1.

## 3 HFS in hydrogen- and lithium-like ions

With the reduction of the CMWFs and the factorization of the n-body matrix elements at our disposal we have performed calculations for the hyperfine splitting (HFS) of heavy ions . In table 6 the DCM results are compared to ground-state HFS splittings calculated for a point nucleus and with results which take the finite spatial distribution of the nuclear charge into account (Breit-Schawlow correction). Additionally, QED radiative corrections  are listed and combined with the DCM results. The pure DCM splittings agree remarkably well with the experimental values, which are given in the last row of the table, while the wavelengths obtained after adding the QED contributions show a systematic shift to larger wavelengths. To clarify this still open point and in order to obtain further informations about the hyperfine interaction of the nucleus with the electronic cloud of high-$`Z`$ ions , we have also performed calculations for the HFS of the lithium-like ions. The result of the calculation  of the HFS for 209Bi80 is given and compared with other theories in table 7. In table 7 the experimental result of Ref.  is also given. A new experiment will be performed at GSI (2003) in order to localize the resonance within an even smaller error. The boiling of the QED vacuum terms given in table 7, i.e. the electron-positrons contributions up to now not considered by other theories, may of course also generate a double counting effect with the QED contributions in lithium-like ions . In the DCM the mixed-mode (b) states generate contributions to the HFS that may be responsible for a double counting effect with the QED perturbative calculation. The terms (b) are not considered by any perturbative calculation performed for the magnetic structure of nuclei.

## 4 Elastic proton scattering on exotic nuclei

The transformation coefficients find also their applications in the calculation of the matter distribution of halo nuclei . Here the calculation are performed for 6Li and 11Li considering a large dimensional space and treating consistently the interaction of the valence particles with the core excitation. For 11Li, in order to reproduce the three-body force we have performed calculations with three dressed neutrons , which are solutions of the symbolic eigenvalue equation given in Fig. 2. The results of the calculations are given in Fig. 3 and Fig. 4. In Fig. 3, we plot: (left) the experimental matter distribution calculated in Ref.  with the BDCM and the experimental matter distribution; and (right) the p-6Li experimental scattering cross sections  compared with the results obtained in Ref.  for the theoretical cross section. In Fig. 4, we plot: (left) the calculated  and the experimental matter distributions; (right) the calculated p-11Li scattering cross sections  compared with the experimental one . For 6Li the calculated and the experimental matter radii given in Fig. 3 are in good agreement with the experimental matter radius given in the figure. The calculated matter radius of 11Li is 3.64 fm  also in good agreement with the experimental value of 3.62(19) fm . The BDCM and the DCM results obtained for the scattering cross sections are in good agreement with the data of Ref. . For 11Li the agreement is better for momentum transfer smaller than .04 (GeV/c)2, as for momentum transfer between 0.04-0.05 (GeV/c)2. This is probably due to high-lying core states, which are not included in the present DCM calculations, but are of increasing importance at higher momentum transfer. The consideration of these core states will increase drastically the dimension of the CMWFs; in the calculations presented here already about two hundred components were used. However, such calculations presently under investigation become easily feasiable with the use of cluster transformation coefficients.

## 5 Proton polarization

The scattering cross section of low-energy photons by a structureless charged system is given by the Thomson formula which derives the photon scattering cross section only in terms of linear contributions which are proportional to the frequency of the incoming photon. A second class of terms which are linear in the frequency of the incident photon and take care of the magnetic moment and spin of the target have also been included as corrective factor to the Thomson scattering. Additionally to these linear effects, however, the cross section get contributions from terms proportional to the square of the photon frequency. These terms, are the polarization terms of Ref.  and modify the cross section accordingly to the square bracket of the following equation:

$`{\displaystyle \frac{\delta \sigma }{\delta \mathrm{\Omega }}}=({\displaystyle \frac{\delta \sigma }{\delta \mathrm{\Omega }}})_{point}+\left[r\left({\displaystyle \frac{E_{sc}}{E_i}}\right)^2{\displaystyle \frac{E_{sc}E_i}{(\mathrm{}c)^2}}\left({\displaystyle \frac{\alpha +\beta }{2}}(1+\mathrm{cos}\theta )^2+{\displaystyle \frac{\alpha \beta }{2}}(1\mathrm{cos}\theta )^2\right)\right]`$ (26)
In Eq. (26), the first term $`(\frac{\delta \sigma }{\delta \mathrm{\Omega }})_{point}`$ give the scattering of photons from a point-like hadron corrected with the size-spin contributions, r is the classical radius of the proton, $`E_i`$ is the incidental photon energy, $`E_{sc}`$ the photon scattering energy, $`\theta `$ the scattering angle, $`\alpha `$ the electric- and $`\beta `$ the magnetic-polarizability ı.e.: terms quadratic in the photon frequency. In Fig. 5, we plot the polarization cross section using for $`\alpha `$ and $`\beta `$ the values calculated with the model independent dispersion sum rule in Ref. :

$`(\alpha +\beta )=(14.2\pm 0.03)\times 10^{49}m^3`$ (27)
These values are around those measured for the proton electric- and magnetic-polarization in Ref. . The plotted polarization terms are ca. twelve order of magnitude smaller then the Compton term. Under this consideration with the present parameters of the PHELIX-laser in mind it is difficult to perform measurements to obtain the proton polarizability within a better resolution. However, additionally to the polarization effect produced by the energy of the incoming photon field, an intrinsic polarization effect  may contribute to the scattering cross section. This would result from effects analogous to the boiling of the QED terms. The three quarks are exciting from the vacuum in a nonlinear mechanism the $`q\overline{q}`$ condensate. The three dressed quarks are also eigenvalues of the equation symbolically represented in Fig. 2. Calculations are in progress using for the quark the harmonic oscillator base proposed in Ref. . The QCD boiling terms should then give rise to a radial distribution with an halo shape as calculated in Fig. 4 for 11Li. De facto the intrinsic polarization effect increases the polarizability of the proton and may produce a sizable effect for the present parameters of PHELIX-laser. Therefore, considering this new source of polarization, experiments can reveal unexpected results.

## 6 Conclusion

The linearization approximation and the cluster transformation coefficients are important tools in solving nonlinear systems such as those characterizing n dressed particles: the linearization ansatz defines a collective hamiltonian that is selfconsistently solvable for the model eigenvalues, while the transformation coefficients provide easy computation for the complex n-body matrix elements which are the input to the collective eigenvalue equation. Within the linearization ansatz and the transformation coefficients, new theoretical insight have been obtained for the physics of interacting particles which coexist with the excitations of the model vacuum.

## References


text:
# Doing numerical cosmology with the Cactus code

## 1 Introduction

The core of modern cosmology is the theory of General Relativity, (-), a 4-dimensional theory involving one dimension of time and three of space, having as field equations the Einstein equations :

$$R_{ij}\frac{1}{2}g_{ij}R+\lambda g_{ij}=\frac{8\pi G}{c^4}T_{ij}$$
(1)
where $`\lambda `$ is the cosmological constant, $`R_{ij}`$ the Ricci tensor, $`R`$ the Ricci scalar, $`g_{ij}`$ the space-time metric, $`T_{ij}`$ the stress-energy tensor, G the gravitational constant, $`c`$ the speed of light and $`i,j=0,1,2,3`$. Numerical Relativity (NR) is concerned with the study of numerical solutions of the Einstein’s equations for the gravitational field (), and when we apply these equations to cosmological models we have numerical cosmology as a branch of NR. Einstein equations are an extremely complicated system of coupled, non-linear, partial differential equations and solving them numerically makes enormous demands on the processing power and memory of a computer. The Cactus code () has been mainly designed as a computational toolkit (freely available for the scientific community) for simulating different systems of partial differential equations, as Einstein equations are The code can be used to simulate fully 3-dimensional systems with strong gravitational fields: collapsing gravitational waves, colliding black-holes, neutron stars, and other violent astrophysical processes generating gravitational waves.

As a branch of NR, numerical cosmology can use the Cactus code for his purposes. Although Einstein field equations reduce to a small number of equations in most cosmological models (-), trying to simulate more complex models (as models with cosmological constant or having matter fields coupled with gravity) it becomes difficult to compose special numerical codes for every specific case. Thus our purpose was to develop a piece of code to be used inside a general numerical machinery for solving Einstein equations, as the Cactus code is. This article is a report on how we extended the Cactus code for numerical cosmology. In this purpose we developed a new thorn of the Cactus code (“thorns” are called the computer code packages embedded in Cactus code for some specific application). The starting point of this new thorn, called “Cosmo” was an existent thorn of the Cactus code (the “Exact” thorn) previously reported in () for testing the Cactus code on exact solutions of the Einstein equations.

We investigated some cosmological models being exact solutions of the Einstein equations in our previous article (). Here the scale factor of the universe ($`R(t)`$) and thus the components of the metric tensor $`g_{ij}`$ are known analytic function of time, given in the code at the initial time as initial data. But in our new Cosmo thorn, the main point is that we have to use the general Robertson-Walker metric, having the scale factor of the universe ($`R(t)`$) as a unknown function to be evolved during the time evolution of the Einstein equations through the Cactus code. Thus the only job to be done by the Cosmo thorn is to provide good initial data (that means at the initial time $`t_0`$ considered the actual time of the universe, for example) and then let the Cactus code to solve the Einstein equations through the time. In addition it is necessary to prescribe convenient boundary values for the functions and variables involved. This is necessary because the Cactus code has a limited boundary methods implemented, not appropriate for our numerical cosmology purposes.

The article is organized as follows : the next section 2, presents some general facts about the way we built the Cosmo thorn, how we solved the problems exposed above and the necessary theoretical cosmology background. Section 3 is dedicated to a matter dominated universe model based on Robertson-Walker metric (), with pressure-less matter - “dust” and without cosmological constant. Models with cosmological constant are presented in the next section 4 where we concentrated again only on models having matter as a pressure-less dust. A wide larger class of possible initial data is possible including the so called “coasting cosmologies”, models not having as a starting point a Big Bang. In all cases studied and reported in this article we investigated the convergence and the stability of the Cactus code, in the same way we presented in detail in .

As a major conclusion we have to point out the good second order convergent and stable behavior of the code, even in the cases we studied both forward and backward time evolutions. Future perspective are open for developing this thorn for inflationary models. It is in our view to study those models having one or more scalar field coupled with gravity to control the behavior of the model and to simulate an accelerated late time evolution, actually called “cosmic acceleration” being well proved experimentally in the last years. Thus together with other similar thorns being in preparation at AEI, we shall have a new class of applications for the Cactus code, a major instrument for new studies in numerical cosmology and in cosmology, in general.

Being included in the list of thorns of the Cactus code, actually the files for our thorn are available by request, directly at the author’s e-mail address and in the near future it will be included in the CVS repository of the Cactus code (see ) for free download.

The simulations were performed mainly on single processor machines, using both the AEI computer network (SGI or Dec machines with UNIX operating system) and a Pentium III machine with a UNIX FreeBSD operating system at the West University of Timişoara. Some of the simulations, with similar results, were also done on the Origin 2000 supercomputer at the AEI.

Through this article and in the Cactus code we shall use geometrical units with $`G=c=1`$.

## 2 The Cactus code, the Robertson-Walker metric and the Cosmo Thorn

As we mentioned earlier, in modern cosmology we are using the Robertson-Walker metric (RW), namely (-)

$$ds^2=c^2dt^2+R(t)^2\left[\frac{dr^2}{1kr^2}+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)\right]$$
(2)
or in isotropic coordinates :

$$ds^2=c^2dt^2+\frac{R(t)^2}{(1+\frac{k}{4}r^2)^2}\left[dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)\right]$$
(3)
where

$`r^{}={\displaystyle \frac{2r}{1+\sqrt{1kr^2}}}`$
as a generic metric for describing the dynamics of the universe through the Einstein equations (1). Here $`k`$ is a constant with arbitrary value, positive (for closed universes), negative (for open universe) and zero for flat universes. Usually, this constant is taken $`1`$, $`1`$ or $`0`$ respectively. When we consider all the matter in the universe as a perfect fluid, the stress-energy tensor may be given as

$$T_{ij}=(ϵ+p)u_iu_j+pg_{ij}$$
(4)
and is necessary to prescribe a state equation, i.e. a relation between the pressure $`p`$ and the density of the universe $`\rho `$. For matter dominated universes usually the pressure is $`p=0`$ and this will be the case we concentrated in our simulations.

As our main purpose was to build a new thorn for the Cactus code, generically named from now one “The Cosmo thorn” for evolving the RW metric above in all three cases for $`k=1,0`$ and $`+1`$. The main problem was to to ”convince” Cactus to evolve, without having the time behavior of $`R(t)`$ function (actually, $`R(t)`$ is a function of the thorn). Although we used the same strategy as in the above mentioned Exact thorn (see ) here, not having an analytic function of time for $`R(t)`$ the solution was to provide Cactus code all initial data he needs, namely $`R(t)`$ and $`\dot{R}(t)=dR(t)/dt`$ at the initial time $`t_0`$ we choose ! As a result, new routines were added for calculating the extrinsic curvature  components at the initial time, namely :

$$K_{\alpha \beta }=2R(t)^2\frac{\dot{R}(t)}{R(t)}\delta _{\alpha \beta }=2R(t)^2H(t)\delta _{\alpha \beta }$$
(5)
where $`H(t)`$ is the Hubble constant function (!). Note that we have here as lapse function, $`N=1`$ and a shift vector as $`N_\alpha =0`$, $`\alpha ,\beta =1,2,3`$ as always in RW cosmology.

Another item we concentrated on was the boundary problem : - Cactus code and the thorns which solve numerically the Einstein equations (namely ADM\_BSSN and ADM - see ) have not implemented proper boundary conditions for the RW metric. In our previous article  on Exact thorn we used “flat” boundary condition but this is not appropiate here ! To solve this problem we used a simple hint : $`R(t)`$ and $`H(t)`$ are functions of time only, thus they are constant on the numerical grid at one time, having the same value inside the grid and on the boundary. As a result, we added a new boundary method through a new value of the parameter bound of the ADM\_BSSN thorn (see  and ), called external which activates some new routines we composed in this purpose, having as a generic name recover. In these routines the spatial metric components $`g_{\alpha \beta }`$ and the extrinsic curvature components $`K_{\alpha \beta }`$ on the boundaries are “recovered” from the values of the $`R(t)`$ and $`H(t)`$ on the interior points by calculating them using equations (2) and (5). Also these new boundary conditions are necessary during the intermediate steps of integration used by ADM\_BSSN thorn in the so-called “IterativeCN.F” routine (where a three steps Cranck-Nicholson integration method is applied). For this purpose we slightly modified the “IterativeCN.F” routine for “injecting” our specific “recover” boundary values. Actually this was possible through a facility of the Cactus code which makes possible to call different routines from different thorns using special “include” commands (see ).

## 3 Friedman-Robertson-Walker (FRW) cosmologies

By “FRW cosmologies” we mean those cosmological models based on the above RW metric and being solutions of the Einstein equations (1) without cosmological constant ($`\lambda =0`$) - see for example () and (). Restricting ourselves only to matter dominated universes (i.e. $`p(t)=0`$) as we mentioned earlier, and solving the conservation law of the perfect fluid which emulates the matter, we have for the stress-energy tensor components

$$T_{\alpha \beta }=0\text{ and }T_{00}=\frac{c^2C}{R(t)^3}$$
(6)
where $`C`$ is a constant depending on the the mass density of the universe at the initial time $`t_0`$. Solving the Einstein equations for this constant and at the initial time $`t_0`$ where $`R(t_0)=R_0`$ and $`H(t_0)=H_0`$, we have

$$C=\mathrm{\Omega }_0\rho _cR_0^3=2q_0\rho _cR_0^3=q_0c^2\frac{3H_0^2R_0^3}{4\pi GR(t)^3}$$
(7)
where $`q_0`$ is the deceleration parameter at the initial time, i.e.

$$q(t)=\frac{\ddot{R}(t)}{H(t)^2R(t)}$$
(8)
and $`\mathrm{\Omega }_0`$ is the so-called density factor, defined as

$$\mathrm{\Omega }_0=\frac{\rho _0}{\rho _c}$$
(9)
$`\rho _0`$ being the initial density of the universe and $`\rho _c`$ the critical density (i.e. the density of a flat three-dimensional universe where $`k=0`$). It is proven () that we have :

$$\mathrm{\Omega }_0=2q_0\text{ and }\rho _c=\frac{3H_0^2}{8\pi G}$$
(10)
In order to provide initial data for the Cactus code we need to solve the Einstein equations at the initial time for the initial Hubble constant $`H_0`$ (to calculate the values of the extrinsic curvature through equation (5)) and for the initial scale factor $`R_0`$. It is known that we have here three solutions, depending on the geometric factor $`k=0,1,1`$. Thus, for $`k=0`$ (flat three-dimensional spacetime) where

$`q_0={\displaystyle \frac{1}{2}}\text{ and }\mathrm{\Omega }_0=1`$
we have (see LABEL:8 eq. 4.37)

$$H_0=\frac{2}{3t_0}$$
(11)
and we shall take $`R_0=1`$, as a value for calculating, for now one the relative scale factor, i.e. $`R(t)/R_0`$ which will be our numerical output. For the case of closed universe, $`k=1`$ we have

$`q_0>{\displaystyle \frac{1}{2}}\text{ and }\mathrm{\Omega }_0>1`$
In this case (as in the next one for $`k=1`$) we do not have an analytic solution for the scale factor, but for the initial time $`t_0`$ we can easily obtain that (LABEL:8 eq. 4.51) :

$$H_0=\frac{q_0}{(2q_01)^{3/2}}\left[cos^1\left(\frac{1q_0}{q_0}\right)\frac{\sqrt{2q_01}}{q_0}\right]\frac{1}{t_0}$$
(12)
and the initial scale factor is :

$$R_0=\frac{c}{H_0\sqrt{2q_01}}$$
(13)
Finally, for the open universes ($`k=1`$) where we have

$`q_0<{\displaystyle \frac{1}{2}}\text{ and }\mathrm{\Omega }_0<1`$
As it is shown ( eq. 4.64), here we have :

$$H_0=\frac{q_0}{(12q_0)^{3/2}}\left[\frac{\sqrt{12q_0}}{q_0}ln\left(\frac{1q_0+\sqrt{12q_0}}{q_0}\right)\right]\frac{1}{t_0}$$
(14)
and the initial scale factor is :

$$R_0=\frac{c}{H_0\sqrt{12q_0}}$$
(15)
With these results reminded, we can now proceed to explain how we produced our simulations. Taking as initial data the initial time ($`t_0`$) in billion years, being an internal parameter of the Cactus code, namely CCTK\_intitial\_time, we then have all the variables (in our geometrical units with $`c=G=1`$, of course) rescaled in billion years. We also have the initial value of the deceleration parameter, $`q_0`$ denoted in the code as robson\_q. As output we have the scale factor, more precisely the relative scale factor $`R(t)/R_0`$ (denoted in the code as raza ) and the Hubble parameter function H(t) (denoted in the code as hubble. For evolving these two functions we produced two internal routines of the thorn, which mainly calculates their values from the general output values provided by the Cactus code and the ADM\_BSSN thorn (as the extrinsic curvature and the three-dimensional metric components).

Here we shall present the results we obtained on evolving a matter dominated universe (i.e pressure-less matter, $`p=0`$) based on the RW metric. Some explanations are necessary in order to clarity our figures ( see also ()). Here and in the next figures, “normalized” means the L2 norm of the respective function calculated using all the values of the function on the computational grid at one time and the respective output Cactus files are denoted with \_nm2 extension.

Cactus code is evolving the Einstein equations using the 3+1 decomposition of space-time (both in ADM or BSSN evolution methods -see () and ()). Thus the Einstein equations can be split in two groups : the dynamic equations (for the time derivatives of the three-dimensional metric and extrinsic curvature) and the constraint equations (the Hamiltonian constraint and the Momentum constraint) - see  and . The constraint equations are satisfied during all time evolution of the system. Thus one of the main tests on the convergence in the Cactus code is given by the time behavior of the Hamiltonian constraint (an output of the Cactus code through a thorn called ADMConstraints, namely ham). In some of the next figures (here and in the next sections) we show the convergence of the L2 norm of the Hamiltonian constraint for different number of iterations, using two different resolutions on the grid, one with $`14^3`$ and the second with $`28^3`$ points (both grids cover the same region of spacetime, so the grid with more points has a smaller value of the finite difference interval $`\mathrm{\Delta }x`$). Notice that the Hamiltonian constraint for a true solution of the Einstein equations should be equal to zero. Finite differencing errors imply that the numerical solution will have a non-vanishing value of the Hamiltonian constraint. For a consistent finite difference approximation of the Einstein equations, we should expect the Hamiltonian constraint to approach zero as the resolution is increased. For a second order approximation, the value of the Hamiltonian constraint should go down by a factor of four when the resolution is doubled. This was the method to test the convergence of the code in all examples presented in this article. We have obtained good second order convergence in all examples as is shown in the next figures.

First of our figures (1 and 2) are pointing out the convergence of the code for different examples of FRW cosmologies. The time is scaled in billion years (as for all figures in this section) starting with an initial time of 20 billion years. As it is obvious from these figures a good second order convergence is visible, even for long time evolutions. For the case of the closed universe ($`k=1`$) we had to restrict our simulations only till around 9000 iterations (for a grid with $`14^3`$ points) because in this case the universe is evolving to a Big Crunch (opposite to the Big Bang) as it is pointed out in the right panel of the next figure (3). Here the scale factor is evolving to a maximum having the value

$`R_{max}={\displaystyle \frac{2q_0}{(2q_01)^{3/2}}}{\displaystyle \frac{c}{H_0}}`$
twice of the actual value $`R_0`$ (we used here $`q_0=1`$). For spatially flat universes we used, of course, $`q_0=0.5`$ and for open universes we had $`q_0=0.25`$.

Figure (3) contains also in his left panel the evolution of the Hubble parameter function in time. We plotted here, for convenience the modulus of this function being calculated through the L2 norm as is output by the code.

The same comment is valuable for the left panel of the figure (4) where we can point out also the coincidence of the minimum zero value of the Hubble parameter with the maximum of the scale factor from the previous figure (3). The right panel of (4) is dedicated to the simultaneous plot of the three cases, pointing out the differences between the geometries of the three models.

Last of our figures from this section, (5) contains the plots of the relative scale factor time evolution for different closed matter dominated FRW universes, having different deceleration parameter (i.e. different density parameters). It is visible here that as the deceleration parameter (density parameter) is increasing so is decreasing the time-life of the closed universe between a Big-Bang and a Big-Crunch.

## 4 FRW cosmologies with cosmological constant

Here we shall deal with cosmologies based on the above RW metric but having a non-vanishing cosmological constant. In spite of the failure of the the initial tentative of Einstein to use the cosmological constant as a trick to fit general relativity with, at that time presumed, static universe, the cosmological constant is again used in modern cosmology in order to establish more realistic models of the universe, as the inflationary models and, recently, the so called cosmic acceleration. Early attempts (see ()) on using numerical routines for solving the Einstein equations for these models were reported. Because our goal is to implement in the Cactus code a thorn for doing similar simulations, but on a larger scale and fully three-dimensional we shall refer to these previous results only for comparing our results. Meanwhile we shall use here the same notations and scaling of the parameters as in () for prescribing the initial data sets for the Cactus code.

This time I will proceed in a slightly different way than in the previous section, scaling all our variables and parameters in terms of the initial value of the Hubble parameter function $`H_0`$, which is not completely free in our universe, having been measured to within a factor of 2. Thus for the cosmological constant $`\lambda `$ for instance, we use $`\lambda /H_0^2`$ as a basic parameter, denoted in the code with robson\_l. Of course, for dimensional scaling reasons we prescribed at the initial time $`H_0=1`$. Thus, for example, at the output we have the relative value of the Hubble function $`H(t)/H_0`$. Meanwhile we need a parameter to describe the density of the universe, and this time we used $`\mathrm{\Omega }_0`$ instead of the deceleration parameter (remark that here we have no more the relation $`\mathrm{\Omega }_0=2q_0`$). We have also the Gaussian curvature $`K(t)`$ defined as

$`K(t)={\displaystyle \frac{k}{R(t)^2}}`$
which will be, of course identical to the geometric factor $`k`$ at the initial time when we take the initial scale factor $`R_0=1`$ as we done in our simulations. Following the line indicated in () we can define the density factor as

$$\mathrm{\Omega }(t)=\frac{8\pi G}{3}\frac{\rho (t)}{H(t)^2}\text{ and at the initial time }\mathrm{\Omega }_0=\frac{8\pi G}{3}\frac{\rho _0}{H_0^2}$$
(16)
where $`\rho _0`$ is the initial density of the matter in the universe. Now solving the Einstein equations at the initial time and after some algebraic manipulations we obtain (for $`c=G=1`$, $`R_0=1`$) :

$$K_0=k=H_0^2\left(\mathrm{\Omega }1+\frac{1}{3}\frac{\lambda }{H_0^2}\right)$$
(17)
and for the stress-energy tensor (where we included also the term containing the cosmological factor, as we did in ()) :

$$T_{ij}=T_{ij}^0\frac{c^4}{8\pi G}\lambda g_{ij}\text{ where }T_{00}^0=\frac{3}{8\pi G}\mathrm{\Omega }_0H_0^2\frac{R_0^3}{R(t)^3}\text{ ; }T_{\alpha \beta }^0=0\alpha ,\beta =1,2,3$$
(18)
Now our numerical results are pointed out in the next figures. First we show, in figure (6) the time evolution of the scale factor for some cosmological models, having the density factor $`\mathrm{\Omega }`$ being $`0.1`$ (left panel) and $`1.0`$ (right panel). The models are differentiated by their respective cosmological constant (more precisely the factor $`\lambda /H_0^2`$). For the time scale we used here Hubble time units, i.e. the time is expressed as $`t/t_H`$ where $`t_H`$ is the Hubble time : $`t_H=H_0^1`$. In these simulations, as in the next ones we started from the actual time (taken here $`t_0=0`$) and we investigated the time behavior of the model in both directions of time, forward, as well as backward in time.

Next figure (7) shows, in his left panel some special cases of FRW models with cosmological constant, where the scale factor shows off a different evolution : namely it starts somewhere back in time with a certain finite value, it decreases in time, then it has a period (longer or shorter) of constant value, evolving finally as an expanding model (of De-Sitter type, for example). These models, having a period of constant value of the scale factor, are sometimes called “coasting cosmologies” and are investigated as an alternative (realistic or not) to the standard Big-Bang model, not being generated from a Big-Bang, as can be seen from our simulations. The right panel of the figure (7) contains the time evolution of the Hubble constant function (the L2 norm, as it is output by the Cactus code, otherwise the values after the null minimum value being negative ones) for a model with $`\mathrm{\Omega }=0.1`$ and $`\lambda /H_0^2=5`$.

We must point out that we performed here also several convergence tests, as in the previous section, obtaining again good second order convergence. As an example we shall present in the last figure (8) the convergence test for $`14^3`$ and $`28^3`$ points on the grid for a model having again $`\mathrm{\Omega }=0.1`$ and $`\lambda /H_0^2=5`$ for an evolution forward in time (left panel) and backward in time (right panel).

## Acknowledgments

The author is deeply indebted to Prof. E. Seidel for patience, continuous support and encouragement, and to prof. Schutz for his understanding and help. The friendly hospitality and environment the author experienced at AEI during his visits is kindly acknowledged.


text:
# Generation of Large Photon-Number Cat States using Linear Optics and Quantum MemoryTo appear in the Proceedings of the Sixth International Conference on Quantum Communication, Measurement and Computing (QCMC’02), Boston, July 22-26, 2002.

## Abstract

A recursive method for producing path-entangled states of light is presented. These states may find applications in quantum lithography and high-precision interferometric measurements. The required resources are single-photon sources, linear optics components, and photodetectors. Adding a quantum memory greatly enhances the yield in comparison with the previously known schemes.

Entangled states play a fundamental role in quantum information processing as they are basic ingredients in various tasks such as quantum teleportation, quantum key distribution, or quantum computing . Entanglement has also been shown to be a useful resource for high-precision frequency measurements or quantum lithography . In this paper, we will focus on the entangled states needed in the latter applications, that is, photon-number path-entangled states or *cat states*. These states are of the form

$$|N_+_{a,b}=\frac{1}{\sqrt{2N!}}(a^N+b^N)|0,$$
(1)
where $`a^{}`$ and $`b^{}`$ denote the bosonic creation operators for two modes of the electromagnetic field, and $`|0`$ denotes the vacuum state. Deterministic schemes to produce cat states using $`\chi ^{(3)}`$ nonlinearity have been proposed . Unfortunately, the presently available $`\chi ^{(3)}`$ couplings are about $`10^{16}`$ times weaker than the required values, which makes these schemes out of the scope of current technology. Recently, however, alternative schemes that only require photon-number state sources, linear optics components (beam splitters and phase shifters) and photodetectors have been proposed for $`N=4`$ first, and then for an arbitrary photon number $`N`$. These schemes are unfortunately *probabilistic*, and the probability $`p(N_+)`$ of a successful state preparation scales exponentially poorly with $`N`$ for all of these schemes. Typically, if $`N`$ is even, $`p(N_+)`$ scales as $`c^N`$ with the constant $`c=\sqrt{2}e`$ for the scheme  and $`c=e`$ in . This is so because these schemes work by adding (or subtracting) the photons one by one (or two by two), each step having a non-unity probability of success. Such a multiplicative process (the yield is a product of success probabilities) naturally scales exponentially.

In this paper, we will demonstrate that using a *tree-like* structure for the linear optics circuit makes it possible to significantly enhance the yield. We will describe a recursive scheme that generates a state of the form of Eq. (1) when $`N`$ is a power of $`2`$ with an approximate yield of

$$(4\sqrt{\pi })^{\mathrm{log}_2N}N^{(1\mathrm{log}_2N)/4}.$$
(2)
Thus, in contrast with the previous schemes, the scaling is *sub-exponential* in $`N`$, which exploits the fact that only $`\mathrm{log}_2N`$ recurrence levels are needed in the tree. As we will see, the price to pay is that a quantum memory (or fast optical switches) is required in addition to the other optical components.

The basic operation in our method is the transformation $`T_N`$, which probabilistically transforms two cat states of $`N`$ photons into one cat state containing $`2N`$ photons, that is, $`T_N:|N_+_{a,b}|N_+_{c,d}|2N_+_{a,b}`$ with probability $`p(T_N)`$. The transformation $`T_N`$ is realized as follows. First, change the phase of mode $`d`$ as $`de^{i\pi /N}d`$. Next, combine the modes $`a`$ and $`c`$ at a $`50`$:$`50`$ beam splitter and modes $`b`$ and $`d`$ at another $`50`$:$`50`$ beam splitter. Defining the action of a $`50`$:$`50`$ beam splitter on two impinging modes $`a`$ and $`b`$ as $`a(aib)/\sqrt{2}`$ and $`b(ia+b)/\sqrt{2}`$, we obtain the 4-mode state

$$\frac{1}{2^{N+1}N!}\left((a^{}+ic^{})^N+(b^{}+id^{})^N\right)\left((ia^{}+c^{})^N(ib^{}+d^{})^N\right)|0.$$
(3)
Finally, measure the number of photons contained in modes $`c`$ and $`d`$. The transformation succeeds if *no* photon is detected in these two modes, in which case we get the (unnormalized) state

$$\frac{1}{2^{N+1}N!}\left((a^{})^{2N}(b^{})^{2N}\right)|0.$$
(4)
Transforming the phase of mode $`b`$ as $`be^{i\pi /2N}b`$ then yields the desired state $`|2N_+_{a,b}`$ with a success probability $`p(T_N)=\frac{(2N)!}{2^{2N+1}(N!)^2}`$. Upon using the Stirling formula, we get

$$p(T_N)\frac{1}{\sqrt{4\pi N}}\left(1+O(1/N)\right).$$
(5)
Remarkably, $`p(T_N)`$ does not decrease exponentially with $`N`$ as if the photons had been added one by one, but only as $`1/\sqrt{N}`$, so the successive steps of our iterative process keep a reasonably good efficiency as $`N`$ increases.

Consider first a naive scheme that does not use a quantum memory. Starting from $`N`$ one-photon states $`|1_+`$, which are easily obtainable with a single-photon source and a beam splitter, one probabilistically produces $`N/2`$ states $`|2_+`$, then $`N/4`$ states $`|4_+`$, and so on until one gets the desired state $`|N_+`$. What is the overall probability $`p(N_+)`$ to produce $`|N_+`$ with this procedure? Solving the recurrence equation $`p(2N_+)=p(N_+)^2p(T_N)`$, we get

$$p(N_+)=\frac{2N!}{(2N)^N}\sqrt{8\pi N}(2e)^N\left(1+O(1/N)\right).$$
(6)
The reason why this success probability decreases exponentially with $`N`$ is that we require all the $`(N1)`$ individual transformations $`T`$ at all recurrence levels to succeed simultaneously. Since each transformation $`T`$ has a non-unity probability of success, the probability that they all succeed is exponentially small.

Let us now show that we can strongly improve on this by using a quantum memory. We start with $`|1_+^{M_1}`$ where the number $`M_1`$ of states we need to prepare initially will be determined below. Suppose that, after a few steps, we have obtained the state $`|n_+^{M_n}`$, assuming $`M_n`$ is an even number. We carry out the transformation $`T_n`$ on each pair of states, which produces on average $`M_{2n}`$ states $`|2n_+`$ where $`M_{2n}=\frac{M_n}{2}P(T_n).`$ We then discard the states for which the transformation failed and only keep those for which it succeeded (there are, on average, $`M_{2n}`$ such states). We use these states to construct the state $`|4n_+`$, and so on. If we want to produce $`M_N`$ cat states $`|N_+`$ (where $`N`$ is a power of 2) with a probability of order one, then, solving the above recurrence using Eq. (5), we find that the average number $`M_1`$ of states $`|1_+`$ that are needed initially is

$$M_1(4\sqrt{\pi })^{\mathrm{log}_2N}N^{(\mathrm{log}_2N1)/4}\left(M_N+𝒪(\sqrt{M_N})\right),$$
(7)
which immediately implies Eq. (2). This shows that the resources (i.e., the number of single photons, of operations $`T_n`$, of beam splitters, of photodetections, etc.) needed for generating a cat state increase only sub-exponentially with $`N`$. This recursive scheme thus requires much less operations than the naive scheme, but the dynamics is now conditional on whether the operations $`T_n`$ succeed or not. Such a conditional dynamics requires fast optical switches in order to store the photons in a loop, or, alternatively, a quantum memory in which the photon states can be stored. These extra technological resources make our scheme much more demanding than those based only on linear optics, single photon sources, and photodetectors only , but $`\chi ^{(3)}`$ couplings are still not required. Note that we may start our recursion from two-photon cat states $`|2_+`$ (instead of $`|1_+`$) since they can be deterministically generated from two one-photon states using a beam splitter, resulting in an extra factor of 4 in the yield Eq. (2).

Our scheme is sensitive to various imperfections such as the losses in the beam splitters, imperfect mode-matching, etc. Let us focus on one of them, namely the non-unity detector efficiency $`\eta <1`$. The probability that the detector does not click when it has absorbed $`j`$ photons is given by $`(1\eta )^j`$. At first sight, it seems that the main source of errors will be due to the events where no photon is detected in modes $`c`$ and $`d`$ while a photon has actually been absorbed in either one of the detectors. By expanding Eq. (3) in a series in $`c^{}`$ and $`d^{}`$, we find that the probability that there is a single photon in mode $`c`$ or $`d`$ is $`p(c^{}c+d^{}d=1)=\frac{2N}{2^{2N}}`$. This probability is thus exponentially smaller than $`p(T_N)`$, so this event can in fact be neglected in good approximation. Next, if we calculate the probability that two photons are in mode $`c`$ or $`d`$, $`p(c^{}c+d^{}d=2)=\frac{(2N2)!}{2^{2N+1}(N1)!^2}\frac{P(T_N)}{4}`$, we find that, for large $`N`$, the probability that one decides that $`T_N`$ succeeded while it actually failed is approximately given by $`(1\eta )^2/4`$. So we conclude that our scheme is robust with respect to detector inefficiencies up to the second order in $`(1\eta )`$.

We thank Stefano Pironio and JiangFeng Du for contributing to stimulating discussions and the European Science Foundation for its financial support. S.I. acknowledges support from the FRIA foundation (Belgium). S.M. is Research Associate of the FNRS (Belgium). This work was partly funded by the European Union under the project EQUIP (IST-FET programme).