text:
# 1 Introduction

## 1 Introduction

The electromagnetic polarizabilities $`\overline{\alpha }`$ and $`\overline{\beta }`$ are fundamental characteristics of the bound system. Their magnitudes depend not only on the quantum numbers of the constituents, but also on the properties of the interaction between these constituents. Therefore, the experimental and theoretical investigation of the electromagnetic polarizabilities are of a great importance. In particular, their prediction and the comparison with experimental data may serve as a sensitive tool for tests of hadron models. Correspondingly, a large number of researchers have been attracted by this fascinating possibility. The electromagnetic polarizabilities can be obtained from the low-energy Compton scattering amplitude. In the lab frame the amplitude of Compton scattering on the compound system of total angular momentum $`S=0,\mathrm{\hspace{0.17em}1}/2`$ up to $`O(\omega ^2)`$ terms reads

$$T=T_{Born}+\overline{\alpha }\omega _1\omega _2\mathit{ϵ}_1\mathit{ϵ}_2^{}+\overline{\beta }(𝐤_1\times \mathit{ϵ}_1)(𝐤_2\times \mathit{ϵ}_2^{}),$$
(1)
where $`\omega _i,𝐤_i,`$ and $`\mathit{ϵ}_i`$ are the energy, momentum, and polarization vector of incoming ($`i=1`$) and outgoing ($`i=2`$) photons ($`\mathrm{}=c=1`$). The contribution $`T_{Born}`$ corresponds to the amplitude of Compton scattering off a point-like particle with spin, mass, charge, and magnetic moment equal to those of the compound system. For spin $`S1`$ the $`O(\omega ^2)`$ part of the Compton scattering amplitude contains additional terms, proportional to quadrupole and higher multipoles of the bound system . In particular, for $`S=1`$ there is a contribution, proportional to the quadrupole moment operator.

The investigation of electromagnetic polarizabilities is interesting not only for systems, bound by the electromagnetic interaction, like atoms, but also for those, bound by strong interaction, such as atomic nuclei  or hadrons . At present there are many different approaches used for the description of the electromagnetic polarizabilities of hadrons: the MIT bag model , the nonrelativistic quark model , the chiral quark model , the chiral soliton model  and the Skyrme model . Here we mentioned only a small part of the publications on these topics (see also review ). Though much effort has been devoted to these calculations, all of them can not be considered as completely satisfactory. In particular, there is a problem in the explanation of the magnitudes of proton and neutron electric polarizabilities within a nonrelativistic quark model. It was derived many years ago  that $`\overline{\alpha }`$ can be represented as a sum

$$\overline{\alpha }=\frac{2}{3}\underset{n0}{}\frac{|n|𝐃|0|^2}{E_nE_0}+\mathrm{\Delta }\alpha =\alpha _{}+\mathrm{\Delta }\alpha ,$$
(2)
where $`𝐃`$ is the internal electric dipole operator, $`|0`$ and $`|n`$ are the ground and excited states in terms of internal coordinates, and $`E_n`$ and $`E_0`$ the corresponding energies. The term $`\mathrm{\Delta }\alpha `$ in $`\overline{\alpha }`$ has a relativistic nature contained and its leading term is equal to

$$\mathrm{\Delta }\alpha =\frac{er_E^2}{3M},$$
(3)
where $`e`$ and $`M`$ are the particle charge and mass, $`r_E`$ is the electric radius defined through the Sachs form factor $`G_E`$111For convenience, our definition of $`r_E^2`$ absorbs a total charge $`e`$ of the system.. The calculation of the quantity $`\alpha _{}`$ in the nonrelativistic quark model without relativistic corrections taken into account leads to the same magnitude of $`\alpha _{}`$ for proton and neutron. Since $`\mathrm{\Delta }\alpha `$ is equal to zero for the neutron but gives a significant contribution to $`\overline{\alpha }`$ for the proton, one has a contradiction between the theoretical prediction of $`\overline{\alpha }`$ for nucleons and their experimental values, since the latter are close to each other. In fact, this approach is not consistent, because there are relativistic corrections to $`\alpha _{}`$ which are of the same order as $`\mathrm{\Delta }\alpha `$.

Starting from second-order perturbation theory one gets the following expression for $`\alpha _{}`$:

$$\alpha _{}=\frac{2}{3}\underset{n0}{}\frac{|n|𝐉|0|^2}{(E_nE_0)^3},$$
(4)
where $`𝐉`$ is the internal electromagnetic current. Using the identity $`𝐉=i[H,𝐃]`$, where $`H`$ is the Hamiltonian, one comes to the form Eq. (2) for $`\alpha _{}`$. The relativistic corrections to $`\alpha _{}`$ in Eq. (4) come from the corrections to wave functions and energies of the ground and excited states, and correction to current $`𝐉`$. Note, that in  (see also ) it was mentioned, that such corrections could exist, but no explicit calculations were made for a realistic system and the importance of these corrections was not realized.

Due to the relation between $`𝐉`$ and $`𝐃`$ it is clear that there also is a relativistic correction to the electric dipole moment operator (see below) which is connected with the appropriate relativistic definition of the center-of-mass coordinate. The neglect of this relativistic correction leads to an incomplete expression for $`\overline{\alpha }`$, and the missing piece which is calculated in the following turns out to be very essential. We expect, that the inclusion of all relativistic corrections allow one to remove the big difference between the predictions of the nonrelativistic quark model for proton and neutron electric polarizabilities due to the difference in $`\mathrm{\Delta }\alpha `$.

The expression for the magnetic polarizability in the nonrelativistic quark model with no exchange and momentum-dependent forces has the form .

$$\overline{\beta }=\frac{2}{3}\underset{n0}{}\frac{|n|𝐌|0|^2}{E_nE_0}\underset{i}{}\frac{e_i^2r_i^2}{6m_i}\frac{𝐃^2}{2M},$$
(5)
where $`𝐌`$ is the internal magnetic dipole operator and the summation in the second term on the r.h.s. is performed over the constituent quarks, $`𝐫_i`$ being the corresponding internal radius vector.

In order to understand the importance of the different relativistic corrections for the polarizabilities, it is useful to consider the example of a system where the relativistic corrections can be obtained analytically. In this paper we calculate the low-energy Compton scattering amplitude for a system of two particles with masses $`m_{1,2}`$ and charges $`e_{1,2}`$, bound by electromagnetic forces. We consider the case $`e_{1,2}^21`$ which provides the validity of the nonrelativistic expansion. We consider in detail the cases of spin $`0`$ and $`1/2`$ of the particles and give the result for general case of arbitrary spins.

## 2 Scattering amplitudes

For the electromagnetic interaction between particles a simple estimate shows that the part $`t`$ of Compton scattering amplitude, proportional to $`\omega ^2`$ has the form:

$$t=\omega ^2a^3\left(c_1+c_2\frac{\epsilon _0}{\mu }+c_3\frac{\epsilon _0^2}{\mu ^2}+\mathrm{}\right),$$
(6)
where $`a=1/(\mu g)`$ is the Bohr radius, $`g=e_1e_2>0`$, $`\mu =m_1m_2/(m_1+m_2)`$ is the reduced mass, $`\epsilon _0=\mu g^2/2`$ is the ground state binding energy in the nonrelativistic approximation, and $`c_i`$ are some quantities, bilinear with respect to $`\mathit{ϵ}_1`$, $`\mathit{ϵ}_2^{}`$ and depending on the ratio of charges and masses. Here, for definiteness, we assume $`\omega =\omega _1`$. The two first terms of this expansion contain the parameter $`g1`$ in the denominators and, therefore, come from the contribution of big distances $`ra`$ (or small momenta $`p1/a=\mu g`$) to the matrix element. These two terms which have no contributions from the Born amplitude are the ones we are going to calculate in this article. Since they are determined by a contribution from big distances (small momenta), it is possible to use the nonrelativistic expansion in the calculations. In fact, the first term is known and is proportional to the electric polarizability $`\alpha _{}`$, calculated in the leading nonrelativistic approximation (see below). Some contributions to the second term are also known, namely, those, containing the magnetic polarizability $`\beta `$ Eq. (5) and the correction $`\mathrm{\Delta }\alpha `$ Eq. (3). These contributions come from the expansion of the photon wave functions over $`\mathrm{𝐤𝐫}\omega /\mu g`$ and from the expansion of the propagator of the system with respect to the photon energy and the center-of-mass kinetic energy of intermediate states. The corresponding results for $`\beta `$ and $`\mathrm{\Delta }\alpha `$ can be obtained using the nonrelativistic Hamiltonian of the system. As was mentioned above, the other source of the contributions to $`c_2`$, which has not been investigated so far, is the relativistic correction to the Hamiltonian of the system and the corresponding corrections to the wave functions, energy levels, and currents. We will obtain the complete result for the second term in the expansion Eq. (6). In the expression for the electromagnetic current we neglect for a while the dependence of the form factors on the momentum transfer. We will take this dependence into account at the consideration of the general case of arbitrary spins.

### 2.1 The system of two spin-$`0`$ particles

Let us consider first the bound state of two spin-$`0`$ particles. In order to calculate the Compton scattering amplitude, it is convenient to put the system into the external electromagnetic field $`𝐀(𝐱,t)`$. In this case the nonrelativistic Hamiltonian has the form

$$\stackrel{~}{H}_{nr}[𝐀]=\frac{𝝅_1^2}{2m_1}+\frac{𝝅_2^2}{2m_2}\frac{g}{|𝐫_1𝐫_2|},$$
(7)
where $`𝝅_i=𝐩_ie_i𝐀(𝐫_i,t)`$. Let us pass to the variables $`𝐫`$ and $`𝐑`$, corresponding to the relative and center-of-mass coordinate:

$$𝐫_1=𝐑+\frac{m_2}{M}𝐫,𝐫_2=𝐑\frac{m_1}{M}𝐫,M=m_1+m_2.$$
(8)
Then, the momenta $`𝐩_i`$ are

$$𝐩_1=\frac{m_1}{M}𝐏+𝐩,𝐩_2=\frac{m_2}{M}𝐏𝐩,$$
(9)
where $`𝐏=i\mathbf{}_R`$ and $`𝐩=i\mathbf{}_r`$. For $`𝐀=0`$ we have

$$\stackrel{~}{H}_{nr}[𝐀=0]=\frac{𝐏^2}{2M}+H_{nr}=\frac{𝐏^2}{2M}+\frac{𝐩^2}{2\mu }\frac{g}{r}.$$
(10)
The first relativistic correction $`H_B`$ (Breit Hamiltonian, see, e.g., ) to Eq. (7) reads

$$\stackrel{~}{H}_B[𝐀]=\frac{1}{8}\left(\frac{\left(𝝅_1^2\right)^2}{m_1^3}+\frac{\left(𝝅_2^2\right)^2}{m_2^3}\right)+\frac{g}{2m_1m_2}\left(\frac{\delta ^{ij}}{r}+\frac{r^ir^j}{r^3}\right)\pi _1^i\pi _2^j.$$
(11)
The first term in Eq. (11) is the correction to the kinetic energy and the second one is the correction due to magnetic quanta exchange, corresponding to the space component of the photon propagator in the Coulomb gauge. If $`𝐀=0`$ then in the center of mass frame where the eigenvalue of the operator $`𝐏`$ is equal to zero we have

$$\stackrel{~}{H}_B[𝐀=0]|_{𝐏=0}H_B=\frac{1}{8}(\frac{1}{m_1^3}+\frac{1}{m_2^3})\left(𝐩^2\right)^2\frac{g}{2m_1m_2}(\frac{\delta ^{ij}}{r}+\frac{r^ir^j}{r^3})p^ip^j,$$
(12)
The terms, containing the operator $`𝐏`$ in the Hamiltonian, determine the contribution of recoil effect to the Compton scattering amplitude. Within the precision of the present calculations these terms should be taken into account only in the Hamiltonian $`\stackrel{~}{H}_{nr}`$ and can be ommited in $`\stackrel{~}{H}_B`$ (see below). The correction $`\delta \epsilon _0`$ to the ground state energy, related to the Hamiltonian $`H_B`$ reads

$$\delta \epsilon _0=0|H_B|0=g^4\left[\frac{5}{8}\mu ^4\left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)+\frac{\mu ^3}{m_1m_2}\right].$$
(13)
Let us start the calculation of the Compton scattering amplitude with the amplitude $`T_{nr}`$ obtained with the use of the nonrelativistic Hamiltonian Eq. (10). This amplitude can be represented as a sum $`T_{nr}=T_{res}+T_s`$ of resonance and seagull parts. The part $`T_{res}`$ is determined by the second order of perturbation theory with respect to the terms in $`\stackrel{~}{H}_{nr}[𝐀]`$, linear in the vector potential $`𝐀`$. In the lab frame it has the form

$`T_{res}`$ $`=`$ $`\psi _0|\mathrm{exp}[i(𝐤_1𝐤_2)𝐑]\mathit{ϵ}_2^{}[{\displaystyle \frac{e_1}{m_1}}𝐩_1\mathrm{exp}(i𝐤_2𝐫_1)+{\displaystyle \frac{e_2}{m_2}}𝐩_2\mathrm{exp}(i𝐤_2𝐫_2)]\times `$ (14)
$`\times `$ $`[\epsilon _0+\omega _1{\displaystyle \frac{𝐏^2}{2M}}H_{nr}]^1\mathit{ϵ}_1\left[{\displaystyle \frac{e_1}{m_1}}𝐩_1\mathrm{exp}(i𝐤_1𝐫_1)+{\displaystyle \frac{e_2}{m_2}}𝐩_2\mathrm{exp}(i𝐤_1𝐫_2)\right]|\psi _0+`$
$`+`$ $`(\mathit{ϵ}_1\mathit{ϵ}_2^{},\omega _1\omega _2,𝐤_1𝐤_2).`$
Here $`\psi _0(𝐫)=\pi ^{1/2}(\mu g)^{3/2}\mathrm{exp}(\mu gr)`$ is the wave function of the ground state, depending on the relative coordinate $`𝐫`$. The final momentum of the bound system is equal to $`𝐤_1𝐤_2`$. Using the relations Eq. (8) and Eq. (9) and making a simple transformation in order to cancel the exponents containing $`𝐑`$, we obtain

$`T_{res}`$ $`=`$ $`\psi _0|\mathit{ϵ}_2^{}[{\displaystyle \frac{e_1}{m_1}}(𝐩+{\displaystyle \frac{m_1}{M}}𝐤_1)\mathrm{exp}(i{\displaystyle \frac{m_2}{M}}𝐤_2𝐫){\displaystyle \frac{e_2}{m_2}}(𝐩{\displaystyle \frac{m_2}{M}}𝐤_1)\mathrm{exp}(i{\displaystyle \frac{m_1}{M}}𝐤_2𝐫)]\times `$ (15)
$`\times `$ $`G(\omega _1)\mathit{ϵ}_1𝐩\left[{\displaystyle \frac{e_1}{m_1}}\mathrm{exp}(i{\displaystyle \frac{m_2}{M}}𝐤_1𝐫){\displaystyle \frac{e_2}{m_2}}\mathrm{exp}(i{\displaystyle \frac{m_1}{M}}𝐤_1𝐫)\right]|\psi _0+`$
$`+`$ $`(\mathit{ϵ}_1\mathit{ϵ}_2^{},\omega _1\omega _2,𝐤_1𝐤_2).`$
Here $`G(\omega )=[\epsilon _0+\omega \omega ^2/2MH_{nr}]^1`$ is the nonrelativistic propagator of the system in the operator form. The seagull amplitude $`T_s`$ is determined by first order of perturbation theory with respect to the terms in $`\stackrel{~}{H}_{nr}[𝐀]`$ which are quadratic in $`𝐀`$. Similar to Eq. (15), we obtain

$$T_s=\mathit{ϵ}_1\mathit{ϵ}_2^{}\psi _0|\left[\frac{e_1^2}{m_1}\mathrm{exp}(i\frac{m_2}{M}(𝐤_1𝐤_2)𝐫)+\frac{e_2^2}{m_2}\mathrm{exp}(i\frac{m_1}{M}(𝐤_1𝐤_2)𝐫)\right]|\psi _0.$$
(16)
Performing the expansion of Eq. (14) and Eq. (16) with respect to $`𝐤_{1,2}`$ and $`\omega _{1,2}`$ up to quadratic terms and using the relation $`\omega _1\omega _2=(𝐤_1𝐤_2)^2/2M`$ , we obtain

$`T_{nr}`$ $`=`$ $`\mathit{ϵ}_1\mathit{ϵ}_2^{}{\displaystyle \frac{(e_1+e_2)^2}{M}}+\mathit{ϵ}_1\mathit{ϵ}_2^{}\omega ^2\left[{\displaystyle \frac{9}{2\mu g^4}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2+{\displaystyle \frac{e_1+e_2}{Mg^2}}\left({\displaystyle \frac{e_1}{m_1^2}}+{\displaystyle \frac{e_2}{m_2^2}}\right)\right]`$ (17)
$`[\mathit{ϵ}_1\times 𝐤_1][\mathit{ϵ}_2^{}\times 𝐤_2]\left[{\displaystyle \frac{1}{2g^2}}\left({\displaystyle \frac{e_1^2}{m_1^3}}+{\displaystyle \frac{e_2^2}{m_2^3}}\right)+{\displaystyle \frac{3}{2Mg^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\right].`$
There is no need here to distinguish between $`\omega _1`$ and $`\omega _2`$ in the $`O(\omega ^2)`$ term. Therefore, we set $`\omega _1=\omega _2=\omega `$ in Eq. (17). The result Eq. (17) is in agreement with Eq. (2) and Eq. (5), with $`\alpha _{}`$ calculated in the nonrelativistic approximation, since in our model

$`𝐃_{nr}=\mu \left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)𝐫,`$
$`\alpha _{nr}={\displaystyle \frac{2}{3}}\psi _0|𝐃_{nr}G_0𝐃_{nr}|\psi _0={\displaystyle \frac{9}{2\mu g^4}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2,`$
$`\mathrm{\Delta }\alpha ={\displaystyle \frac{e_1+e_2}{3M}}\psi _0|e_1r_1^2+e_2r_2^2|\psi _0={\displaystyle \frac{e_1+e_2}{Mg^2}}\left({\displaystyle \frac{e_1}{m_1^2}}+{\displaystyle \frac{e_2}{m_2^2}}\right),`$ (18)
$`\psi _0|𝐌|\psi _{n0}=0,\psi _0|𝐃^2|\psi _0={\displaystyle \frac{3}{g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2,`$
$`\psi _0|\left[{\displaystyle \frac{e_1^2r_1^2}{m_1}}+{\displaystyle \frac{e_2^2r_2^2}{m_2}}\right]|\psi _0={\displaystyle \frac{3}{g^2}}\left({\displaystyle \frac{e_1^2}{m_1^3}}+{\displaystyle \frac{e_2^2}{m_2^3}}\right).`$
Here $`G_0`$ is the reduced Green function in the operator form:

$$G_0=[\epsilon _0H_{nr}+i0]^1(1|\psi _0\psi _0|).$$
(19)
The details of calculations of different matrix elements, containing the operator $`G_0`$ are presented in Appendix A.

We pass now to the calculation of the relativistic corrections to the electromagnetic polarizabilities, connected with the Breit Hamiltonian $`\stackrel{~}{H}_B`$ Eq. (11). We perform the calculations in the same way as at the derivation of Eq. (14) and Eq. (16), but for the Hamiltonian $`\stackrel{~}{H}=\stackrel{~}{H}_{nr}+\stackrel{~}{H}_B`$. At the calculation of these corrections within our accuracy the terms of the Hamiltonian $`\stackrel{~}{H}_B[𝐀]`$, quadratic in $`𝐀`$ do not contribute to the electromagnetic polarizabilities, i.e. the seagull contribution from $`\stackrel{~}{H}_B[𝐀]`$ is absent. In the corrections to the electromagnetic polarizabilities from the resonance part of the amplitude, we can take the second order of expansion with respect to $`\omega `$ of the operator Green function and put $`𝐤_{1,2}=0`$ elsewhere. This means that within our accuracy we can neglect in $`\stackrel{~}{H}_B[𝐀]`$ the terms containing the total momentum $`𝐏`$ and replace the exponents in the photon wave function by unity. Therefore, the terms in the Hamiltonian $`\stackrel{~}{H}_B[𝐀]`$, linear in $`𝐀`$, can be represented in the form $`𝐀(0)𝐉_B`$, where

$$𝐉_B=\left(\frac{e_1}{m_1^3}\frac{e_2}{m_2^3}\right)\frac{𝐩^2}{2}𝐩\frac{g(e_1e_2)}{2m_1m_2r}\left(𝐩+\frac{𝐫}{r^2}(𝐫𝐩)\right)$$
(20)
is the correction to the operator of the total internal current $`𝐉`$ :

$$𝐉=𝐉_{nr}+𝐉_B=(e_1/m_1e_2/m_2)𝐩+𝐉_B.$$
(21)
Let us now discuss the relativistic correction to the electric dipole moment operator. In the lab frame it is equal to $`𝐃_{tot}=e_1𝐫_1+e_2𝐫_2`$. For the total Hamiltonian of the system the following relation holds

$`i[\stackrel{~}{H}_{nr}[𝐀=0]+\stackrel{~}{H}_B[𝐀=0],𝐃_{tot}]=`$ (22)
$`=e_1{\displaystyle \frac{𝐩_1}{m_1}}(1{\displaystyle \frac{𝐩_1^2}{2m_1^2}}){\displaystyle \frac{ge_1}{2m_1m_2r}}(𝐩_1+{\displaystyle \frac{𝐫}{r^2}}(𝐫𝐩_1))+(12).`$
The r.h.s. of Eq. (22) is nothing but the operator of total current in the lab frame, containing the total momentum $`𝐏`$. This can be verified by differentiating the Hamiltonian $`\stackrel{~}{H}_{nr}[𝐀]+\stackrel{~}{H}_B[𝐀]`$ over $`𝐀`$ at $`𝐀=0`$. Therefore, there are no relativistic corrections to the total electric dipole operator $`𝐃_{tot}`$. The internal electric dipole moment operator is defined as

$$𝐃=𝐃_{tot}(e_1+e_2)𝐑_{cm},$$
(23)
where $`𝐑_{cm}`$ is the center-of-mass vector. This vector is defind in such a way that it satisfies the following relations:

$$[𝐑_{cm},𝐏]=i,i[H_{tot},𝐑_{cm}]=\frac{𝐏}{H_{tot}},$$
(24)
where $`H_{tot}`$ is the total relativistic Hamiltonian of the system, and $`𝐏`$ is the total momentum. Within our accuracy the second relation in Eq. (24) reads:

$$i[\stackrel{~}{H}_{nr}[𝐀=0]+\stackrel{~}{H}_B[𝐀=0],𝐑_{cm}]=\frac{𝐏}{M}\left(1\frac{\stackrel{~}{H}_{nr}[𝐀=0]}{M}\right).$$
(25)
It is known (see, e.g. , ) that there is a relativistic correction to $`𝐑_{cm}`$ in classical electrodynamics. For the case of two particles, interacting due to electromagnetic field, the corresponding operator which satisfies the relations Eq. (24) has the form

$$𝐑_{cm}=𝐑+\frac{1}{2M}\left(\{𝐫_1,\frac{𝐩_1^2}{2m_1}\frac{g}{2r}\}+\{𝐫_2,\frac{𝐩_2^2}{2m_2}\frac{g}{2r}\}\{𝐑,\frac{𝐩_1^2}{2m_1}+\frac{𝐩_2^2}{2m_2}\frac{g}{r}\}\right).$$
(26)
Here we took into account the first relativistic correction and use the notation $`\{a,b\}=ab+ba`$. In terms of the variables $`𝐫`$ and $`𝐩`$ (see Eq. (8) and Eq. (9)) we obtain :

$$𝐑_{cm}=𝐑+\frac{(m_2m_1)}{2M^2}\left(\{𝐫,H_{nr}\}+g\frac{𝐫}{r}\right),$$
(27)
where the term proportional to the total momentum $`𝐏`$ is omitted. Substituting this expression into Eq. (23) we obtain the relativistic correction to the internal electric dipole moment:

$$𝐃_B=\frac{(e_1+e_2)(m_1m_2)}{2M^2}\left(\{𝐫,H_{nr}\}+g\frac{𝐫}{r}\right).$$
(28)
Note that, as should be the case, the operator of total internal current $`𝐉=𝐉_{nr}+𝐉_B`$ satisfies within our accuracy the relation

$$𝐉=i[H_{nr}+H_B,𝐃_{nr}+𝐃_B],$$
(29)
where $`𝐃_{nr}`$ is defined in Eq. (2.1), and $`H_{nr}+H_B`$ is the internal part of Hamiltonian (see Eq. (10) and Eq. (12)).

Let us write down now the corrections to the $`O(\omega ^2)`$ term of Compton scattering amplitude, related to the Breit Hamiltonian and the corresponding current. The correction due to $`𝐉_B`$ reads

$$t_c=\omega ^2\psi _0|[\mathit{ϵ}_2^{}𝐉_BG_0^3\mathit{ϵ}_1𝐉_{nr}+\mathit{ϵ}_2^{}𝐉_{nr}G_0^3\mathit{ϵ}_1𝐉_B]|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$
(30)
The $`O(\omega ^2)`$ correction to the amplitude, connected with the expansion of the propagator with respect to $`H_B`$, has the form

$$t_p=\omega ^2\psi _0|\mathit{ϵ}_2^{}𝐉_{nr}[G_0^2H_BG_0^2+G_0^3H_BG_0+G_0H_BG_0^3]\mathit{ϵ}_1𝐉_{nr}|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$

The contribution due to the correction to wave function is

$$t_w=\omega ^2\psi _0|[\mathit{ϵ}_2^{}𝐉_{nr}G_0^3\mathit{ϵ}_1𝐉_{nr}G_0H_B+H_BG_0\mathit{ϵ}_2^{}𝐉_{nr}G_0^3\mathit{ϵ}_1𝐉_{nr}]|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$

At last, the contribution corresponding to the correction to the ground state energy reads:

$$t_e=3\omega ^2\delta \epsilon _0\psi _0|\mathit{ϵ}_2^{}𝐉_{nr}G_0^4\mathit{ϵ}_1𝐉_{nr}|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$

In order to calculate the matrix elements in Eq. (30)-Eq. (2.1) it is convenient to use the following relations (see Appendix A):

$`G_0𝐫\psi _0={\displaystyle \frac{𝐫}{g}}\left({\displaystyle \frac{r}{2}}+a\right)\psi _0,G_0r𝐫\psi _0={\displaystyle \frac{𝐫}{g}}\left({\displaystyle \frac{r^2}{3}}+{\displaystyle \frac{5ar}{6}}+{\displaystyle \frac{5a^2}{3}}\right)\psi _0,`$
$`G_0r^2\psi _0={\displaystyle \frac{1}{g}}\left({\displaystyle \frac{r^3}{3}}+ar^3{\displaystyle \frac{11a^3}{2}}\right)\psi _0,`$ (31)
$`G_0r^3\psi _0={\displaystyle \frac{1}{g}}\left({\displaystyle \frac{r^4}{4}}+{\displaystyle \frac{5ar^3}{6}}+{\displaystyle \frac{5a^2r^2}{2}}{\displaystyle \frac{155a^4}{8}}\right)\psi _0,`$
where $`a=1/\mu g`$. Using (2.1) and also the relation $`𝐩=i\mu [H_{nr},𝐫]`$ , we get the following results for the contributions Eq. (30)-Eq. (2.1)

$`t_c`$ $`=`$ $`{\displaystyle \frac{2\mu \omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}}{9g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)\left[20\left({\displaystyle \frac{e_1}{m_1^3}}{\displaystyle \frac{e_2}{m_2^3}}\right)+{\displaystyle \frac{89}{4}}{\displaystyle \frac{(e_1e_2)}{\mu m_1m_2}}\right],`$ (32)
$`t_p`$ $`=`$ $`{\displaystyle \frac{\mu ^2\omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}}{g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[{\displaystyle \frac{1061}{288}}\left({\displaystyle \frac{1}{m_1^3}}+{\displaystyle \frac{1}{m_2^3}}\right)+{\displaystyle \frac{25}{3\mu m_1m_2}}\right],`$
$`t_w`$ $`=`$ $`{\displaystyle \frac{\mu ^2\omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}}{g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{1}{m_1^3}}+{\displaystyle \frac{1}{m_2^3}}\right)+{\displaystyle \frac{14}{9\mu m_1m_2}}\right],`$
$`t_e`$ $`=`$ $`{\displaystyle \frac{\mu ^2\omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}}{g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2{\displaystyle \frac{129}{4}}\left[{\displaystyle \frac{5}{8}}\left({\displaystyle \frac{1}{m_1^3}}+{\displaystyle \frac{1}{m_2^3}}\right)+{\displaystyle \frac{1}{\mu m_1m_2}}\right].`$
Representing the sum of all contributions in Eq. (32) as $`\omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}\alpha _B`$ we obtain the following result for $`\alpha _B`$ :

$$\alpha _B=\frac{1}{g^2}\left(\frac{e_1}{m_1}\frac{e_2}{m_2}\right)^2\left(\frac{121}{6\mu }\frac{113}{4M}\right)\frac{(e_1+e_2)}{Mg^2}\left(\frac{e_1}{m_1}\frac{e_2}{m_2}\right)\frac{(m_1m_2)}{2m_1m_2}.$$
(33)
If one starts the calculation of $`\alpha _{}`$ from Eq. (2), then it is easy to check that the first term in Eq. (33) corresponds to the sum of corrections due to modification of wave function, propagator, and ground state energy. The second term in Eq. (33) corresponds to the contribution due to the relativistic correction to the electric dipole moment Eq. (28). Therefore, this correction appears due to the correct description of the center-of-mass motion. It is seen that the first term in Eq. (33) has the same dependence on charges as $`\alpha _{nr}`$ while the second term is similar to $`\mathrm{\Delta }\alpha `$ , Eq. (2.1). Taking a sum of $`\alpha _B`$ , $`\alpha _{nr}`$, and $`\mathrm{\Delta }\alpha `$, we come to the following result for $`\overline{\alpha }`$ for the system under consideration:

$`\overline{\alpha }`$ $`=`$ $`{\displaystyle \frac{1}{\mu g^4}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[{\displaystyle \frac{9}{2}}g^2\left({\displaystyle \frac{121}{6}}{\displaystyle \frac{113\mu }{4M}}\right)\right]+`$ (34)
$`+{\displaystyle \frac{(e_1+e_2)}{Mg^2}}\left[{\displaystyle \frac{3}{2}}\left({\displaystyle \frac{e_1}{m_1^2}}+{\displaystyle \frac{e_2}{m_2^2}}\right){\displaystyle \frac{(e_1+e_2)}{2m_1m_2}}\right].`$
Thus, the relativistic corrections to $`\alpha _{}`$ has reduced to a renormalization of $`\alpha _{nr}`$ and essentially to a modification of $`\mathrm{\Delta }\alpha `$ , Eq. (3). One can expect that the last statement is valid not only for the system under consideration. Indeed, due to the definition of $`𝐃`$, the correction to $`\overline{\alpha }`$ related to the modification of $`𝐑_{cm}`$ is proportional to $`Q/M`$ , where Q is the total charge of the system, and, therefore, has the same structure as $`\mathrm{\Delta }\alpha `$.

As a nontrivial test of our method of calculation we checked the fulfilment of the low-energy theorem for the Compton scattering amplitude. At $`\omega =0`$ this amplitude should have the form

$$T(\omega =0)=\mathit{ϵ}_1\mathit{ϵ}_2^{}\frac{(e_1+e_2)^2}{E_0}\mathit{ϵ}_1\mathit{ϵ}_2^{}\frac{(e_1+e_2)^2}{M+\epsilon _0}\mathit{ϵ}_1\mathit{ϵ}_2^{}\frac{(e_1+e_2)^2}{M}\left(1\frac{\epsilon _0}{M}\right),$$
(35)
where $`E_0=M+\epsilon _0`$ is the mass of the system. It is interesting, that the term in r.h.s. of Eq. (35), proportional to the nonrelativistic energy $`\epsilon _0`$, appears as a contribution of terms from the Breit Hamiltonian $`\stackrel{~}{H}_B[𝐀]`$, which we checked by explicit calculations (see Appendix B).

### 2.2 The system of a spin-$`0`$ particle and a spin-$`1/2`$ particle

Let the first particle have the spin $`1/2`$ and the second particle have the spin $`0`$. Then we should add the term

$$\delta \stackrel{~}{H}_{nr}[𝐀]=\frac{e_1(1+\kappa _1)}{m_1}𝐬_1𝐇$$
(36)
to the nonrelativistic Hamiltonian $`\stackrel{~}{H}_{nr}[𝐀]`$, Eq. (7). Here $`𝐇`$ is the external magnetic field, $`𝐬_1=𝝈_1/2`$ is the spin operator of the first particle, and $`\kappa _1`$ is its anomalous magnetic moment in units $`e_1/2m_1`$. There is also some additional contribution $`\delta \stackrel{~}{H}_B[𝐀]`$ to $`\stackrel{~}{H}_B[𝐀]`$, Eq. (12) (see, e.g. ). The terms of $`\delta \stackrel{~}{H}_B[𝐀]`$ linear in $`𝐬_1`$ as well as Eq. (36) determine the $`O(\omega )`$ terms of the Compton amplitude. These terms are well-known and follow, together with the $`\omega `$-independent term, from the low-energy theorem . As it was explained in the previous subsection, it is sufficient within our accuracy to account for the Breit Hamiltonian only in the long-wave limit, i.e. at $`\omega _{1,2}=0`$ in order to obtain the $`O(\omega ^2)`$ terms of Compton amplitude. In this limit the Hamiltonian $`\delta \stackrel{~}{H}_B[𝐀]`$ reads

$$\delta \stackrel{~}{H}_B[𝐀]=\frac{e_1e_2(1+2\kappa _1)}{2m_1^2}\left(\pi \delta (𝐫)+\frac{1}{r^3}𝐬_1(𝐫\times 𝝅_1)\right)+\frac{e_1e_2(1+\kappa _1)}{m_1m_2r^3}𝐬_1(𝐫\times 𝝅_2).$$
(37)
The explicit calculation shows that the contribution of $`\delta \stackrel{~}{H}_{nr}[𝐀]`$ given by Eq. (36) and the terms in Eq. (37) linear in $`\sigma _1`$ do not lead to any contributions to $`c_{1,2}`$ in Eq. (6), i.e. they can be neglected in the calculation of polarizabilities within our accuracy. In particular, there are no terms linear in the spin in the quantities $`c_{1,2}`$, which is in agreement with the general conclusion on the absence of terms $`O(\omega ^2)`$ linear in spin in the non-Born part of the Compton amplitude . The only term which should be taken into account in addition to those considered in the previous subsection, is the spin-independent term in Eq. (37) (Darwin term):

$$\delta _DH_B=\frac{\pi g(1+2\kappa _1)}{2m_1^2}\delta (𝐫).$$
(38)
It follows from Eq. (37) that there is no correction to the current associated with the Hamiltonian $`\delta _DH_B`$. Using the expressions (2.1-2.1) with the replacement $`H_B\delta _DH_B`$ and the relations Eq. (2.1) we obtain

$`\delta _Dt_c=0,\delta _Dt_w={\displaystyle \frac{5\mu \omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}}{8m_1^2g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2(1+2\kappa _1),`$ (39)
$`\delta _Dt_p=0,\delta _Dt_e={\displaystyle \frac{129\mu \omega ^2\mathit{ϵ}_1\mathit{ϵ}_2^{}}{8m_1^2g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2(1+2\kappa _1).`$
As a result, the correction to the electric polarizability associated with the Breit Hamiltonian in the system of spin-$`0`$ and spin-$`1/2`$ will be the sum of $`\alpha _B`$ Eq. (33) and

$$\delta \alpha _B=\frac{31\mu }{2m_1^2g^2}\left(\frac{e_1}{m_1}\frac{e_2}{m_2}\right)^2(1+2\kappa _1).$$
(40)
### 2.3 The system of two spin-$`1/2`$ particles

In the case of two spin-$`1/2`$ particles it is necessary to account for two Darvin terms in addition to the Breit Hamiltonian Eq. (12), corresponding to both particles

$$\delta _DH_B=\frac{\pi g(1+2\kappa _1)}{2m_1^2}\delta (𝐫)+\frac{\pi g(1+2\kappa _2)}{2m_2^2}\delta (𝐫)$$
(41)
and the Hamiltonian, corresponding to spin-spin interaction :

$$\delta _sH_B=\frac{g(1+\kappa _1)(1+\kappa _2)}{m_1m_2}\left[\frac{3(𝐧𝐬_1)(𝐧𝐬_2)𝐬_1𝐬_2}{r^3}+\frac{8\pi }{3}\delta (𝐫)𝐬_1𝐬_2\right],$$
(42)
with $`𝐧=𝐫/r`$. It is more convenient to rewrite $`\delta _sH_B`$ in terms of the total spin operator $`𝐒=𝐬_1+𝐬_2`$:

$$\delta _sH_B=\frac{g(1+\kappa _1)(1+\kappa _2)}{2m_1m_2}\left[\frac{3n_in_jQ_{ij}}{2r^3}+4\pi \delta (𝐫)\left(\frac{2}{3}𝐒^21\right)\right],$$
(43)
where the operator $`Q_{ij}`$, quadratic in $`𝐒`$, is equal to

$$Q_{ij}=S_iS_j+S_jS_i\frac{2}{3}𝐒^2\delta _{ij}.$$
(44)
Note that in such a system as positronium it is necessary to add the contribution of the annihilation diagram, which results in the replacement $`(2𝐒^2/31)(7𝐒^2/61)`$ in the coefficient of the $`\delta `$-function in Eq. (43) (of course, in this case $`m_1=m_2`$, $`e_1=e_2`$). As in the previous subsection, the terms proportional to the $`\delta `$-function in Eq. (41) and Eq. (43) give the contribution $`\delta \alpha _B`$, which should be added to $`\alpha __B`$, Eq. (33):

$$\delta \alpha _B=\frac{31\mu }{2g^2}\left(\frac{e_1}{m_1}\frac{e_2}{m_2}\right)^2\left[\frac{1+2\kappa _1}{m_1^2}+\frac{1+2\kappa _2}{m_2^2}+\frac{4(1+\kappa _1)(1+\kappa _2)}{m_1m_2}\left(\frac{2}{3}S(S+1)1\right)\right].$$
(45)
Here we replaced $`𝐒^2`$ by its eigenvalue $`S(S+1)`$, where $`S=0,1`$ is the total spin of the system. The term in Eq. (45) containing the tensor operator $`3(𝐧𝐒)^2𝐒^2`$ determines the contribution to the $`O(\omega ^2)`$ part of the Compton amplitude, which has the form

$$t_{(tensor)}=\omega ^2\alpha _Tϵ_1^iϵ_2^jQ_{ij},$$
(46)
where $`\mathrm{}`$ denotes the averaging over the spin part of the wave function. Of course, $`t_{(tensor)}`$ vanishes if $`S=0`$. Since there is no correction to the current or to the energy of the ground state due to the tensor part of $`\delta _sH_B`$, the contributions to $`\alpha _T`$ come only from the corrections to the propagator and to the wave function. Using Eq. (2.1), Eq. (2.1), and the relations (see Appendix A)

$`G_0(3r_ir_jr^2\delta _{ij})|\psi _0`$ $`=`$ $`{\displaystyle \frac{3r_ir_jr^2\delta _{ij}}{g}}\left({\displaystyle \frac{r}{3}}+{\displaystyle \frac{a}{2}}\right)|\psi _0`$ (47)
$`G_0(3r_ir_jr^2\delta _{ij})r|\psi _0`$ $`=`$ $`{\displaystyle \frac{3r_ir_jr^2\delta _{ij}}{g}}\left({\displaystyle \frac{7a^2}{8}}+{\displaystyle \frac{7ar}{12}}+{\displaystyle \frac{r^2}{4}}\right)|\psi _0,`$
we obtain

$$\alpha _T=\frac{47(1+\kappa _1)(1+\kappa _2)}{40Mg^2}\left(\frac{e_1}{m_1}\frac{e_2}{m_2}\right)^2.$$
(48)
In the system of two spin-1/2 particles there is a big paramagnetic contribution to the magnetic polarizability from the first term in Eq. (5). The main contribution corresponds to the transition from the ground state with the total spin $`S=0`$ to the state with $`S=1`$, with both states having the same angular momenta, $`l=0`$, and radial quantum nambers $`n_r=0`$ (hyperfine splitting). Representing the spin part of the magnetic moment operator in the form

$$𝐌_s=f_1𝐬_1+f_2𝐬_2,f_i=\frac{e_i(1+\kappa _i)}{m_i},$$

and using Eq. (42), we obtain:

$$\beta _1=\frac{3(f_1f_2)^2a^3}{16f_1f_2},a=\frac{1}{\mu g}.$$
(49)
As was pointed out in the previous section, for positronium it is necessary to change the coefficient of $`\delta `$-function in Eq. (43). As a result, the contribution of the first term in Eq. (5) to the magnetic polarizability of positronium is:

$$\beta _1=\pm \frac{3}{7}a^3,$$
(50)
where upper sign corresponds to parapositronium ($`S=0`$), and lower sign to orthopositronium ($`S=1`$).

### 2.4 The system of two particles with arbitrary spins

Let the particles have the spins $`s_{1,2}`$ and magnetic moments $`\mu _{1,2}`$ which we represent in the form

$$\mu _a=\frac{e_as_a}{m_a}(1+\kappa _a),a=1,2.$$
(51)
The electromagnetic current for each particle has the form (see, e.g., )

$$j_\mu =\overline{\psi }(p^{})\left[F_e\frac{p_\nu +p_\nu ^{}}{2m}+\frac{G_m}{2m}\mathrm{\Sigma }_{\mu \nu }q^\nu \right]\psi (p),$$
(52)
where $`q=p^{}p`$. The operator $`\mathrm{\Sigma }_{\mu \nu }`$ is a generalisation of the corresponding matrix for spin $`1/2`$. The indices numerating the particles have been omitted. The quantities $`F_e`$ and $`G_m`$ depend on $`q^2`$ and $`(s_\mu q^\mu )^2`$, where $`s_\mu `$ is the 4-vector of the spin operator. This quantities are normalized as follows:

$$F_e(q=0)=1,G_m(q=0)=1+\kappa .$$
(53)
If we neglect the $`q`$-dependence of the form factors, then, in addition to the Breit Hamiltonian for two spin-$`0`$ particles, Eq. (11), it is necessary to take into account the Hamiltonian Eq. (42) (with the corresponding spin operators) and two other contributions . Namely, the Darwin Hamiltonian

$$\delta _DH_B=\underset{a=1,2}{}\frac{2\pi g}{3m_a^2}(1+2\kappa _a)(s_a+\zeta _a)\delta (𝐫),$$
(54)
$`\zeta =0`$ for integer spin and $`\zeta =1/4`$ otherwise, and the term containing the quadrupole moments of the particles:

$$\delta _QH_B=\underset{a=1,2}{}\frac{g(1+2\kappa _a)\xi _a}{2m_a^2r^3}(3(𝐧𝐬_a)^2𝐬_a^2),$$
(55)
$`\xi =1/(2s1)`$ for integer spin and $`\xi =1/(2s)`$ otherwise. It is clear that all matrix elements can be calculated in the same way as in the previous subsection. The averaging over the spin variables can be done using the following relations

$`SS_z|\left[s_{1i}s_{2j}+s_{2i}s_{1j}{\displaystyle \frac{2}{3}}\delta _{ij}𝐬_1𝐬_2\right]|S,S_z^{}`$ $`=`$ $`A(S,s_1,s_2)SS_z|Q_{ij}|S,S_z^{},`$ (56)
$`SS_z|\left[s_{1i}s_{1j}+s_{1i}s_{1j}{\displaystyle \frac{2}{3}}\delta _{ij}𝐬_1^2\right]|S,S_z^{}`$ $`=`$ $`B(S,s_1,s_2)SS_z|Q_{ij}|S,S_z^{}`$
where $`𝐒=𝐬_1+𝐬_2`$ is the total spin operator, $`Q_{ij}`$ is defined in Eq. (44), and for $`S1`$

$`A(S,s_1,s_2)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }^2+2\mathrm{\Lambda }(\lambda _1+\lambda _2)3(\lambda _1\lambda _2)^2}{2\mathrm{\Lambda }(4\mathrm{\Lambda }3)}}`$ (57)
$`B(S,s_1,s_2)`$ $`=`$ $`{\displaystyle \frac{3\mathrm{\Lambda }^2+\mathrm{\Lambda }(2\lambda _16\lambda _23)+3(\lambda _1\lambda _2)(\lambda _1\lambda _21)}{2\mathrm{\Lambda }(4\mathrm{\Lambda }3)}}.`$
Here $`\mathrm{\Lambda }=S(S+1)`$, $`\lambda _{1,2}=s_{1,2}(s_{1,2}+1)`$ are the eigenvalues of the operators $`𝐒^2`$ and $`𝐬_{1,2}^2`$, respectively. For $`S=0,\mathrm{\hspace{0.17em}1}/2`$ we put $`A=B=0`$. As a result, we obtain the following generalization of Eq. (45) to the case of arbitrary spins:

$`\delta \alpha _B`$ $`=`$ $`{\displaystyle \frac{62\mu }{3g^2}}({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}})^2[{\displaystyle \frac{1+2\kappa _1}{m_1^2}}(s_1+\zeta _1)+{\displaystyle \frac{1+2\kappa _2}{m_2^2}}(s_2+\zeta _2)+`$
$`+{\displaystyle \frac{2(1+\kappa _1)(1+\kappa _2)}{m_1m_2}}(\mathrm{\Lambda }\lambda _1\lambda _2)].`$
The generalization of Eq. (48) is

$`\alpha _T`$ $`=`$ $`{\displaystyle \frac{47\mu }{40g^2}}({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}})^2[{\displaystyle \frac{2(1+\kappa _1)(1+\kappa _2)}{m_1m_2}}A(S,s_1,s_2)+`$
$`+{\displaystyle \frac{1+2\kappa _1}{m_1^2}}\xi _1B(S,s_1,s_2)+{\displaystyle \frac{1+2\kappa _2}{m_2^2}}\xi _2B(S,s_2,s_1)].`$
Let us now take into account the $`q`$-dependence of the electromagnetic form factors of the constituents defined in Eq. (52). We assume, that the scale of variation of these form factors are much larger than the typical momentum transfer $`\mu g`$. In other words, the characteristic size of each constituent is much smaller than the size of the whole system $`a=1/\mu g`$. In this case it is sufficient, within our accuracy, to take $`G_m=1+\kappa `$ and to expand the form factor $`F_e`$ up to quadratic in $`q`$ terms:

$$F_e(q^2,(s^\mu q_\mu )^2)1\frac{r_e^2𝐪^2}{6}+\frac{r_s^2(𝐬𝐪)^2}{2},$$
(60)
where $`r_{e,s}^2`$ are some constants. Multiplying the $`O(𝐪^2)`$ terms in this expression by $`4\pi g/𝐪^2`$ and performing the Fourier transform, we obtain the additional terms in the Hamiltonian

$$\delta _fH_B=\underset{a=1,2}{}\left[\frac{2\pi g}{3}\left(r_{ea}^2r_{sa}^2𝐬_a^2\right)\delta (𝐫)+gr_{sa}^2\frac{3(𝐬_a𝐧)^2𝐬_a^2}{2r^3}\right].$$
(61)
Since the terms in this Hamiltonian have the same structure as above, it is easy to write down the result for the corresponding corrections to polarizabilities:

$`\delta _f\alpha _B`$ $`=`$ $`{\displaystyle \frac{62\mu }{3g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2(r_{e1}^2r_{s1}^2\lambda _1+r_{e2}^2r_{s2}^2\lambda _2),`$ (62)
$`\delta _f\alpha _T`$ $`=`$ $`{\displaystyle \frac{47\mu }{40g^2}}\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[r_{s1}^2B(S,s_1,s_2)+r_{s2}^2B(S,s_2,s_1)\right].`$
If the parameters of the form factors $`r_{e,s}^21/m^2a^2`$, then the contributions Eq. (62) to the polarizabilities are of the same order as $`\alpha _B`$. The first relativistic correction to the Compton scattering amplitude at $`\omega =0`$, Eq. (35), is proportional to $`\epsilon _0=\mu g^2/2`$ and is independent of the spins of the constituents. Then, the correction to the amplitude at $`\omega =0`$ connected with spin-dependent terms in Breit Hamiltonian as well as the Darwin terms (also having the spin origin) should vanish. This statement was checked explicitly (see Appendix B).

Let us consider now the paramagnetic contribution to the magnetic polarizability from the first term in Eq. (5). Let $`s_1s_2`$. Then, the total spin of the ground state is $`S=s_1s_2`$, and the main contribution corresponds to the transition from the ground state to the state with $`S=s_1s_2+1`$, with both states having the same angular momentum, $`l=0`$, and radial quantum namber $`n_r=0`$ (hyperfine splitting). A simple explicit calculation leads to:

$$\beta _1=\frac{(f_1f_2)^2s_2(s_1+1)a^3}{4f_1f_2(s_1s_2+1)^2}.$$
(63)
This term should be added to the diamagnetic contribution $`\beta _{dia}`$ (see Eq. (17)) :

$$\beta _{dia}=\frac{1}{2g^2}\left(\frac{e_1^2}{m_1^3}+\frac{e_2^2}{m_2^3}\right)\frac{3}{2Mg^2}\left(\frac{e_1}{m_1}\frac{e_2}{m_2}\right)^2.$$
(64)
## 3 Conclusion

We have obtained the complete result for the first relativistic corrections to the electromagnetic polarizabilities, including the tensor part which exists for the total spin $`S1`$. We demonstrated that, within our accuracy, this tensor part contains the quadrupole moment of the system and no any higher multipoles. For the system of two spinless particles it is easy to check that the total relativistic correction Eqs. (2.1) and (33) is negative at arbitrary masses and charges. In the general case of non-zero spins and arbitrary anomalous magnetic moments the relativistic correction $`\mathrm{\Delta }\alpha +\alpha _B+\delta \alpha _B`$, where $`\delta \alpha _B`$ is given by Eq. (2.4), can be positive. It is interesting to consider some special cases. The first of them is a hydrogen-like ion. In this case $`e_1=e`$, $`e_2=Ze`$, and $`m_2m_1`$. In the limit $`m_2\mathrm{}`$ the result for electric polarizabilities is independent of the spin and magnetic moment of the nucleus. Neglecting also the anomalous magnetic moment $`\kappa _1`$ of the electron, we obtain from Eq. (17), Eq. (2.1), Eq. (33), and Eq. (45)

$$\overline{\alpha }=\frac{9}{2m^3\alpha _{em}^3Z^4}\frac{14}{3m^3\alpha _{em}Z^2},$$
(65)
where $`\alpha _{em}=e^2=1/137`$ is the fine-structure constant. Note that in this limit the correction $`\mathrm{\Delta }\alpha `$ , Eq. (3), vanishes. The result Eq. (65) is in agreement with that obtained in  with the use of the reduced Green function of the Dirac equation for an electron in a Coulomb field. For the magnetic polarizability at $`s_2=0`$ we have

$$\overline{\beta }=\beta _{dia}=\frac{1}{2m^3\alpha _{em}Z^2}.$$
(66)
For $`s_2=1/2`$ in the limit $`m_2m_1`$ there is a very big contribution from the paramagnetic part of the magnetic polarizability, Eq. (49). In the Compton scattering amplitude this contribution should be taken into account only for photon energies $`\omega `$ much smaller than the energy $`E_{hf}\alpha _{em}^4m_1^2/m_2`$ of the hyperfine splitting. For $`\mu g^2\omega E_{hf}`$ the paramagnetic contribution should be omitted.

Another interesting example is positronium. As we mentioned above, in this case it is necessary to replace $`(2𝐒^2/31)(7𝐒^2/61)`$ in the coefficient of the $`\delta `$-function in Eq. (43) due to the contribution of the annihilation diagram. Putting $`m_1=m_2=m`$, $`e_1=e_2=e`$, and $`\kappa _1=\kappa _2=0`$, we obtain $`\mathrm{\Delta }\alpha =0`$ and the complete result for the polarizabilities

$`\overline{\alpha }`$ $`=`$ $`{\displaystyle \frac{36}{(m\alpha _{em})^3}}+{\displaystyle \frac{1}{6m^3\alpha _{em}}}\times \{{\displaystyle \genfrac{}{}{0pt}{}{1001\text{for }S=0}{735\text{for }S=1}}`$
$`\overline{\beta }`$ $`=`$ $`(1)^S{\displaystyle \frac{3}{56m^3\alpha _{em}^3}}{\displaystyle \frac{4}{m^3\alpha _{em}}}.`$ (67)
As in the previous case, for photon energy $`\omega \alpha _{em}^4m`$ the paramagnetic contribution should be omitted in the Compton scattering amplitude.

For $`S=1`$ (orthopositronium) we also have the tensor polarizability

$$\alpha _T=\frac{47}{20m^3\alpha _{em}}.$$
(68)
Thus, we have shown that the complete set of the first relativistic corrections differs essentially from the commonly used term $`\mathrm{\Delta }\alpha `$. We suppose that for the electromagnetic polarizabilities of hadrons investigated within the constituent quark model an analogous situation may be found.

## Acknowledgments

A.I.M. thanks the members of the Zweites Physikalisches Institut at the University of Göttingen for the warm hospitality during his stay, when a part of this work has been done. One of the authors (M.S.) thanks Deutsche Forschungsgemeinschaft for the suport of this work thought he grants Schu222 and 436RUS113/510.

## Appendix A

In this Appendix we derive the formulas for the result of the action of the operator $`G_0`$, Eq. (19), on the wave function $`|\psi _0`$, multiplied by some polynomial of $`𝐫`$. More precisely, we obtain the expression for $`G_0Y_{lm}(𝐫/r)r^n|\psi _0`$ in the form of the product $`Y_{lm}(𝐫/r)P(r)|\psi _0`$,where $`P(r)`$ is some polynomial. Since the Hamiltonian $`H_{nr}`$ commutes with the operator of angular momentum $`𝐥=𝐫\times 𝐩`$, we can make the following transformation:

$$G_0Y_{lm}(𝐫/r)r^n|\psi _0=Y_{lm}(𝐫/r)G_0^{(l)}(r^n\delta _{l\mathrm{\hspace{0.17em}0}}r^n)|\psi _0,$$
(69)
where $`G_0^{(l)}=[\epsilon _0H_{nr}^{(l)}+i0]^1`$, $`H_{nr}^{(l)}=(2\mu r)^1_r^2r+l(l+1)/(2\mu r^2)g/r`$ is the radial Hamiltonian with the angular momentum $`l`$, and

$$r^n=\psi _0|r^n|\psi _0=\frac{(n+2)!}{2^{n+1}}a^n,$$
(70)
where $`a=1/\mu g`$. In the derivation of Eq. (69) we used the identity

$$(1|\psi _0\psi _0|)Y_{lm}(𝐫/r)r^n|\psi _0=Y_{lm}(𝐫/r)(r^n\delta _{l\mathrm{\hspace{0.17em}0}}r^n)|\psi _0.$$

For our purposes it is sufficient to consider the cases $`nl1`$ for $`l0`$ and $`n1`$ for $`l=0`$. It is easy to check that in these cases one can represent the result of action of $`G_0^{(l)}`$ in r.h.s. of Eq. (69) in the form

$$G_0^{(l)}(r^n\delta _{l\mathrm{\hspace{0.17em}0}}r^n)|\psi _0=\underset{k=0}{\overset{\mathrm{}}{}}C_kr^k|\psi _0,$$
(71)
where $`C_k`$ are some constants to be found. Acting on both sides of this equation with the operator $`\epsilon _0H_{nr}^{(l)}`$ and collecting the coefficients with different powers of $`r`$, we obtain

$`r^n\delta _{l\mathrm{\hspace{0.17em}0}}r^n`$ $`=`$ $`{\displaystyle \frac{l(l+1)}{2\mu }}C_0r^2{\displaystyle \frac{(l1)(l+2)}{2\mu }}C_1r^1+`$
$`+{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{(kl+2)(k+l+3)}{2\mu }}C_{k+2}g(k+1)C_{k+1}\right)r^k.`$
From this relation we can find the coefficients $`C_i`$. For the case $`nl1`$, $`l0`$, we finally obtain

$$G_0Y_{lm}(𝐫/r)r^n|\psi _0=Y_{lm}(𝐫/r)\frac{\left(nl+1\right)!\left(n+l+2\right)!}{g(2/a)^{n+1}\left(n+1\right)!}\underset{k=l}{\overset{n+1}{}}\frac{\left(k1\right)!(2r/a)^k}{\left(kl\right)!\left(k+l+1\right)!}|\psi _0.$$
(73)
For the case $`n1`$, $`l=0`$ we have

$$G_0r^n|\psi _0=\frac{(n+2)!}{g(2/a)^{n+1}}\underset{k=2}{\overset{n+1}{}}\frac{1}{k}\left(\frac{(2r/a)^k}{(k+1)!}\frac{k+2}{2}\right)|\psi _0.$$
(74)
Using the formulas Eq. (Appendix A) and Eq. (73) one can easily calculate all matrix elements needed.

## Appendix B

In this Appendix we check the fulfilment of the low energy theorem. Namely, we reproduce the two first terms of the expansion with respect to $`\epsilon _0/M`$ of the Compton scattering amplitude at $`\omega =0`$:

$$T(\omega =0)\mathit{ϵ}_1\mathit{ϵ}_2^{}\frac{(e_1+e_2)^2}{M}\left(1\frac{\epsilon _0}{M}\right).$$
(75)
In fact, the first term is contained in Eq. (17). In order to obtain the second term, we have to take into account the corrections to the current, seagull, wave function, propagator, and energy due to the Breit Hamiltonian $`\stackrel{~}{H}_B[𝐀]`$. The contribution to the amplitude at $`\omega =0`$ due to $`𝐉_B`$ reads

$$T_c=\psi _0|[\mathit{ϵ}_2^{}𝐉_BG_0\mathit{ϵ}_1𝐉_{nr}+\mathit{ϵ}_2^{}𝐉_{nr}G_0\mathit{ϵ}_1𝐉_B]|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$
(76)
The contribution due to the correction to seagull (the terms in $`\stackrel{~}{H}_B[𝐀]`$ being quadratic in $`𝐀`$) reads

$`T_s`$ $`=`$ $`\psi _0|\{({\displaystyle \frac{e_1^2}{m_1^3}}+{\displaystyle \frac{e_2^2}{m_2^3}})[(\mathit{ϵ}_1𝐩)(\mathit{ϵ}_2^{}𝐩)+(\mathit{ϵ}_1\mathit{ϵ}_2^{}){\displaystyle \frac{𝐩^2}{2}}]+`$
$`+{\displaystyle \frac{g^2}{m_1m_2}}[{\displaystyle \frac{\mathit{ϵ}_1\mathit{ϵ}_2^{}}{r}}+{\displaystyle \frac{(\mathit{ϵ}_1𝐫)(\mathit{ϵ}_2^{}𝐫)}{r^3}}]\}|\psi _0.`$
The contribution connected with the expansion of propagator with respect to $`H_B`$ has the form

$$T_p=\psi _0|(\mathit{ϵ}_2^{}𝐉_{nr}G_0H_BG_0\mathit{ϵ}_1𝐉_{nr}|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$

The contribution due to the correction to wave function is

$$T_w=\psi _0|[\mathit{ϵ}_2^{}𝐉_{nr}G_0\mathit{ϵ}_1𝐉_{nr}G_0H_B+H_BG_0\mathit{ϵ}_2^{}𝐉_{nr}G_0\mathit{ϵ}_1𝐉_{nr}]|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$

At last, the contribution corresponding to the correction to the ground state energy reads:

$$T_e=\delta \epsilon _0\psi _0|\mathit{ϵ}_2^{}𝐉_{nr}G_0^2\mathit{ϵ}_1𝐉_{nr}|\psi _0+(\mathit{ϵ}_1\mathit{ϵ}_2^{}).$$

Using the results of Appendix A, we obtain the following expressions for the corrections:

$`T_c`$ $`=`$ $`{\displaystyle \frac{\mathit{ϵ}_1\mathit{ϵ}_2^{}}{3}}g^2\mu ^2\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)\left[5\mu \left({\displaystyle \frac{e_1}{m_{1}^{}{}_{}{}^{3}}}{\displaystyle \frac{e_2}{m_{2}^{}{}_{}{}^{3}}}\right)+4{\displaystyle \frac{e_1e_2}{m_1m_2}}\right],`$ (78)
$`T_s`$ $`=`$ $`{\displaystyle \frac{\mathit{ϵ}_1\mathit{ϵ}_2^{}}{6}}g^2\mu \left[5\mu \left({\displaystyle \frac{e_{1}^{}{}_{}{}^{2}}{m_{1}^{}{}_{}{}^{3}}}+{\displaystyle \frac{e_{2}^{}{}_{}{}^{2}}{m_{2}^{}{}_{}{}^{3}}}\right)+{\displaystyle \frac{8g}{m_1m_2}}\right],`$
$`T_p`$ $`=`$ $`{\displaystyle \frac{\mathit{ϵ}_1\mathit{ϵ}_2^{}}{12}}g^2\mu ^3\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[7\mu \left({\displaystyle \frac{1}{m_{1}^{}{}_{}{}^{3}}}+{\displaystyle \frac{1}{m_{2}^{}{}_{}{}^{3}}}\right)+{\displaystyle \frac{12}{m_1m_2}}\right],`$
$`T_w`$ $`=`$ $`{\displaystyle \frac{\mathit{ϵ}_1\mathit{ϵ}_2^{}}{6}}g^2\mu ^3\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[9\mu \left({\displaystyle \frac{1}{m_{1}^{}{}_{}{}^{3}}}+{\displaystyle \frac{1}{m_{2}^{}{}_{}{}^{3}}}\right)+{\displaystyle \frac{14}{m_1m_2}}\right],`$
$`T_e`$ $`=`$ $`{\displaystyle \frac{\mathit{ϵ}_1\mathit{ϵ}_2^{}}{4}}g^2\mu ^3\left({\displaystyle \frac{e_1}{m_1}}{\displaystyle \frac{e_2}{m_2}}\right)^2\left[5\mu \left({\displaystyle \frac{1}{m_{1}^{}{}_{}{}^{3}}}+{\displaystyle \frac{1}{m_{2}^{}{}_{}{}^{3}}}\right)+{\displaystyle \frac{8}{m_1m_2}}\right].`$
Summing up these contributions, we get

$$T_B=\mathit{ϵ}_1\mathit{ϵ}_2^{}\frac{(e_1+e_2)^2}{2M^2}\mu g^2,$$
(79)
which is the second term in Eq. (75). Let us consider now the contribution to the Compton amplitude at $`\omega =0`$, connected with the spin-dependent terms and the Darwin terms in Breit Hamiltonian. Note that all these terms are proportional to either $`\delta (𝐫)`$ or to the operator $`(3n_in_j\delta _{ij})/r^3`$. The terms $`\delta (𝐫)`$ give the contributions to $`T_w`$ and $`T_e`$. Using the results of Appendix A, it is easy to show that the sum of these two contributions is zero. The terms $`(3n_in_j\delta _{ij})/r^3`$ give the contributions to $`T_w`$ and $`T_p`$. Again, direct calculations show that they also cancel each other. Therefore, we proved, that the first relativistic correction to the Compton amplitude at $`\omega =0`$ is spin-independent, which is in agreement with the low-energy theorem.


format: text

text:
# Adaptive evolution on neutral networks

## 1 Introduction

Sudden bursts of adaptive activity which punctuate long periods of stagnancy seem to be a common observation in evolving systems. Such *epochal evolution* , has been found in the fossil record , in evolving bacteria , the evolution of tRNA structures  or artificial systems such as digital organisms  and evolutionary optimization . The most complete theoretical analysis of epochal evolution has probably been presented in a series of papers by van Nimwegen and coworkers . The general picture is as follows. A population can easily climb to the nearest local optimum, but escape from there only with difficulty. Once a local optimum has been reached, the population is trapped and experiences a metastable equilibrium. With a relatively low probability, the population can discover a portal genotype leading to the next local optimum, i.e., a genotype with a higher fitness than what is currently present in the population. Once a portal is discovered, the population moves quickly away from its current peak, and towards the new peak. There, it settles down again in equilibrium, until the next portal is discovered.

The above description focuses on the dynamics between local optima, but not on the dynamics at one local optimum. However, in the presence of neutrality, i.e., when a number of genotypes share the same identical replication rate, the dynamics around one local optimum can be quite intriguing. A population does not drift over a set of neutral genotypes, like a set of random walkers, but has a quite different, even to some extent deterministic dynamics. In a completely flat fitness landscape, for example, a population does not assume a Gaussian distribution with a variance that increases over time, as would be expected from a simple diffusion process. Rather, the population stays clustered together and the cluster moves about as a whole . If neutral and deleterious genotypes are mixed, i.e., if each genotype has as direct neighbors in genotype space both neutral and deleterious genotypes, then the population moves to the regions of the most connected neutral genotypes, as long as the mutation rate is small .

In this paper, we study the dynamics of a finite but large population of asexually replicating genetic sequences on a fitness landscape that contains both neutral and strongly deleterious genotypes. We show that the neutral genotypes naturally decompose into disjunct sets, and that an evolving population can be trapped within such sets. Moreover, different sets yield different reproductive success for the populations residing on them. When a population discovers a set of genotypes with higher reproductive success, the population moves over to that set. For the average fitness of the population, such transitions are reflected in a stepwise increase, exactly as it is observed in standard scenarios of epochal evolution. However, the increase observed here is not due to the discovery of faster replicating genotypes, but solely to the discovery of genotypes with increased robustness against mutations. Our results are valid for arbitrary mutation rates, and they generalize the previous findings of van Nimwegen et al. .

## 2 Theory

We assume that there exist two classes of genotypes, those with a relatively high replication rate $`\sigma `$, and the ones with a much lower replication rate. For reasons of simplicity, we assume the latter to have replication rate 0. This assumption is quite common in the literature . Our analysis is similar to the standard treatment of the quasispecies model (see for example Schuster and Swetina , Eigen et al. , Wilke et al. ) and to the work of van Nimwegen et al. . Let the vector $`𝒂`$ define the set of neutral genotypes, i.e., the genotypes with replication rate $`\sigma `$:

$$a_i=\{\begin{array}{cc}1& \text{if }i\text{ has replication rate }\sigma \text{,}\hfill \\ 0& \text{else,}\hfill \end{array}$$
(1)
where $`i`$ runs over all possible genotypes. We assume a discrete time model, and write the average population fitness as $`f`$. In equilibrium, we have

$$𝒙=\frac{\sigma }{f}\text{Q}\text{A}𝒙,$$
(2)
where $`𝒙`$ is the vector of concentrations, the diagonal matrix $`\text{A}=\mathrm{diag}(a_0,a_1,\mathrm{})`$ contains the set of neutral genotypes, and the matrix Q defines the mutation probabilities between different genotypes, i.e., genotype $`j`$ mutates into genotype $`i`$ with probability $`Q_{ij}`$.

In the following, we assume that the genotypes can be represented as sequences of length $`l`$ over an alphabet of $`A`$ different symbols. Moreover, we assume a uniform copy fidelity $`q`$ per symbol. That means, the $`l`$ symbols in a sequence mutate independently from each other, and the substitution probability in one generation is $`1q`$ for each symbol. With this assumption, the mutation matrix can be written as

$$Q_{ij}=q^l\left(\frac{1q}{q(A1)}\right)^{d(i,j)},$$
(3)
where $`d(i,j)`$ is the Hamming distance between sequences $`i`$ and $`j`$. It is useful to introduce the reduced mutation rate $`\stackrel{~}{\mu }`$,

$$\stackrel{~}{\mu }=\frac{1q}{q(A1)},$$
(4)
which allows us to write Q as a sum of matrices,

$$\text{Q}=q^l\underset{k=0}{\overset{l}{}}\stackrel{~}{\mu }^k\text{D}^{(k)}.$$
(5)
The matrices $`\text{D}^{(k)}`$ define the connection graphs at Hamming distance $`k`$ in sequence space, i.e.,

$$D_{ij}^{(k)}=\{\begin{array}{cc}1& \text{if }d(i,j)=k\text{,}\hfill \\ 0& \text{else.}\hfill \end{array}$$
(6)
We insert (5) into (2), and obtain

$$𝒙=\frac{\sigma q^l}{f}\underset{k=0}{\overset{l}{}}\stackrel{~}{\mu }^k\text{D}^{(k)}\text{A}𝒙.$$
(7)
It is useful to introduce the matrices

$$\text{G}^{(k)}=\text{A}\text{D}^{(k)}\text{A}.$$
(8)
These matrices define the connection graphs at Hamming distance $`k`$ for the neutral genotypes. We now disregard all non-neutral sequences, and introduce the concentration vector $`𝒑`$, which holds the concentrations of all neutral sequences. The total number of neutral sequences is then $`P=_ip_i`$. Moreover, we assume all columns and rows corresponding to non-neutral sequences to be deleted from the matrices $`\text{G}^{(k)}`$. From equation (7), we obtain the eigenvalue equation

$$\left(\frac{f}{\sigma q^l}1\right)𝒑=\left(\underset{k=1}{\overset{l}{}}\stackrel{~}{\mu }^k\text{G}^{(k)}\right)𝒑.$$
(9)
Consequently, the equilibrium state of the population is fully determined by the matrix

G $`={\displaystyle \underset{k=1}{\overset{l}{}}}\stackrel{~}{\mu }^{k1}\text{G}^{(k)}`$
$`=\text{G}^{(1)}+\stackrel{~}{\mu }\text{G}^{(2)}+\stackrel{~}{\mu }^2\text{G}^{(3)}\mathrm{}`$ (10)
In the following, we will call matrices such as G *generalized connection matrices*. The difference to the normal connection matrices $`\text{G}^{(k)}`$ as defined in equation (8) is that generalized connection matrices may contain powers of $`\stackrel{~}{\mu }`$, i.e., they define a connection graph with weighted edges, whereas for the graphs defined by $`\text{G}^{(k)}`$, all edges have the same weight.

We can obtain further insight from relating the average population fitness $`f`$ to the average fraction of neutral offspring, $`\nu `$. Under the assumption that the non-neutral sequences have a vanishing replication rate (this assumption is equivalent to neglecting back mutations), we have

$$\nu =\frac{f}{\sigma }.$$
(11)
The probability $`\nu _i`$ for a single sequence $`i`$ to have a neutral genotype as offspring is given by

$`\nu _i`$ $`={\displaystyle \underset{j}{}}a_jQ_{ji}`$
$`=q^l{\displaystyle \underset{k=0}{\overset{l}{}}}\stackrel{~}{\mu }^k{\displaystyle \underset{j}{}}a_jD_{ji}^{(k)}`$
$`=q^l\left(1+{\displaystyle \underset{k=1}{\overset{l}{}}}\stackrel{~}{\mu }^kd_i^{(k)}\right),`$ (12)
where $`d_i^{(k)}`$ gives the number of neutral neighbors at Hamming distance $`k`$ of sequence $`i`$. When we take the average over all viable sequences in the population on both sides of equation (2), we arrive at

$$\nu =q^l\left(1+\underset{k=1}{\overset{l}{}}\stackrel{~}{\mu }^kd^{(k)}\right).$$
(13)
The quantities $`d^{(k)}`$ give the average number of neutral neighbors at Hamming distance $`k`$ in the population. They can be expressed as

$$d^{(k)}=\frac{1}{P}\underset{i}{}p_i\underset{j}{}G_{ij}^{(k)}.$$
(14)
With equation (11), we can rewrite equation (13) as

$$\underset{k=1}{\overset{l}{}}\stackrel{~}{\mu }^kd^{(k)}=\frac{f}{\sigma q^l}1.$$
(15)
We notice that the right-hand-side of (15) is identical to the factor in front of $`𝒑`$ on the left-hand-side of (9). Therefore, we arrive at

$$\left(\underset{k=1}{\overset{l}{}}\stackrel{~}{\mu }^{k1}d^{(k)}\right)𝒑=\left(\underset{k=1}{\overset{l}{}}\stackrel{~}{\mu }^{k1}\text{G}^{(k)}\right)𝒑.$$
(16)
In the limit of a small mutation rate, (16) becomes

$$d^{(1)}𝒑=\text{G}^{(1)}𝒑.$$
(17)
In that case, the equilibrium distribution $`𝒑`$ depends solely on the one-mutant connections of the neutral sequences, and the population neutrality $`d^{(1)}`$ is given by the spectral radius of the first-order connection matrix. Equation (17) corresponds to the result of van Nimwegen et al. . Thus, we find that equation (16) is the generalization of that result to arbitrary mutation rates. For larger mutation rates, the equilibrium distribution is influenced by the higher order connection matrices, and the population neutrality may deviate from the spectral radius of $`\text{G}^{(1)}`$.

### 2.1 Separate neutral networks

Following the conventional nomenclature in the literature , we call a set of neutral sequences that can be transversed by one-point mutations a *neutral network*. The full set of neutral sequences in genotype space will in general decompose into several disjunct such neutral networks. Assume there are $`n`$ disjunct neutral networks. By reordering sequences, we can arrange G in block-matrix form

$$\text{G}=\left(\begin{array}{ccccc}\text{G}_1& & & & \\ & \text{G}_2& & 𝒪(\stackrel{~}{\mu })& \\ & & \text{G}_3& & \\ & 𝒪(\stackrel{~}{\mu })& & \mathrm{}& \\ & & & & \text{G}_n\end{array}\right),$$
(18)
where the matrices $`\text{G}_i`$ are the generalized connection matrices for the different neutral networks, and all off-diagonal terms are at least of the order of $`\stackrel{~}{\mu }`$. Let us further assume that the matrices $`\text{G}_i`$ are ordered with descending spectral radius, i.e., the spectral radius of $`\text{G}_i`$ is larger than the one of $`\text{G}_{i+1}`$ for all $`i`$. For a *finite* population initially residing on network $`n`$, we can then expect the following dynamics. If the mutation rate is small, the population will equilibrate within network $`n`$. The discovery of a sequence that is part of a neutral network $`i<n`$ is very unlikely, due to the small off-diagonal terms. However, eventually such a sequence will be discovered. Since a progeny of that sequence has a smaller probability to fall off its neutral network ($`\text{G}_i`$ has a larger spectral radius than $`\text{G}_n`$), the sequences on network $`i`$ will have a higher reproductive success than the sequences on network $`n`$, and the population will move over to the newly discovered network. There, the population will equilibrate, until the next higher-connected network is discovered. If the mutation rate is large, on the other hand, the off-diagonal terms cannot be considered small anymore. In that case, the discovery of higher connected regions is much more likely, and a population will move straight away into the more densely connected regions of the genotype space. In Sec. 4, we present examples for both of these behaviors.

The above considerations show that the decomposition of the full set of neutral sequences into neutral networks according to the definition given at the beginning of this subsection is somewhat arbitrary. Depending on the mutation rate, it might be justified, for example, to disregard contributions of order $`𝒪(\stackrel{~}{\mu }^2)`$ to the matrix G, but not the ones of order $`𝒪(\stackrel{~}{\mu })`$. In that case, it would be more natural to group the sequences into sets that can be transversed by a combination of one- or two-point mutations. In general, the neutral sequences should therefore be subdivided into sets such that two arbitrary sequences of two disjunct sets are at least a Hamming distance $`k`$ apart, where $`k`$ represents the smallest number of simultaneous mutations that can be considered rare at the respective mutation rate. For the remainder of this paper, we will understand the term *neutral network* in this more general sense.

## 3 A simple exactly solvable landscape

In this section, we study a simple example landscape for which the matrix G can be diagonalized exactly, to all orders in $`\stackrel{~}{\mu }`$. We consider binary sequences of length $`l=2n`$, and break them down into $`n`$ pairs of bits. For each pair, we assume that there are three states (00, 01, 10) which are neutral, and one state (11) which is lethal. Therefore, a sequence which contains at least one pair for which both bits are set to 1 has a fitness of 0, and all other sequences have a fitness of $`\sigma `$. For simplicity, we set $`\sigma =1`$. We will refer to this landscape as the *Neutral Bitpairs* landscape.

For a single pair, the matrix $`\text{G}_1`$ reads (in this section, the subscript $`i`$ in $`\text{G}_i`$ indicates the number of pairs we are considering)

$$\text{G}_1=\left(\begin{array}{ccc}0& 1& 1\\ 1& 0& \stackrel{~}{\mu }\\ 1& \stackrel{~}{\mu }& 0\end{array}\right),$$
(19)
and its largest eigenvalue is $`\lambda _1=(\stackrel{~}{\mu }+\sqrt{8+\stackrel{~}{\mu }^2})/2`$. For a sequence with $`2`$ pairs, the corresponding matrix $`\text{G}_2`$ can be written as a tensor product (see Rumschitzki , Dress and Rumschitzki , Wilke ),

$$\text{G}_2=\stackrel{~}{\mu }^1\left[(\mathrm{𝟏}+\stackrel{~}{\mu }\text{G}_1)(\mathrm{𝟏}+\stackrel{~}{\mu }\text{G}_1)\mathrm{𝟏}\right].$$
(20)
The symbol $`\mathrm{𝟏}`$ stands for the identity matrix in the appropriate matrix space. In general, for a sequence consisting of $`n`$ pairs, we can define the matrix $`\text{G}_n`$ recursively,

$$\text{G}_n=\stackrel{~}{\mu }^1\left[(\mathrm{𝟏}+\stackrel{~}{\mu }\text{G}_1)(\mathrm{𝟏}+\stackrel{~}{\mu }\text{G}_{n1})\mathrm{𝟏}\right].$$
(21)
As a consequence, the largest eigenvalue of $`\text{G}_n`$ reads

$`\lambda _n`$ $`=\stackrel{~}{\mu }^1[(1+\stackrel{~}{\mu }\lambda _1)^n1],`$
$`=\stackrel{~}{\mu }^1[(1+{\displaystyle \frac{\stackrel{~}{\mu }^2}{2}}+{\displaystyle \frac{\stackrel{~}{\mu }}{2}}\sqrt{8+\stackrel{~}{\mu }^2})^n1].`$ (22)
The average population fitness in this landscape follows from equation (9). We obtain

$$f=q^{2n}(1+\frac{\stackrel{~}{\mu }^2}{2}+\frac{\stackrel{~}{\mu }}{2}\sqrt{8+\stackrel{~}{\mu }^2})^n.$$
(23)
For small mutation rates, we can write the average fitness as

$$f=q^{2n}(1+\sqrt{2}n\stackrel{~}{\mu }),$$
(24)
and the average population neutrality becomes

$$d^{(1)}=\sqrt{2}n.$$
(25)
Let us now compare the full solution for the average population fitness to the approximation given by van Nimwegen et al. ,

$$f=1\mu \left(1\frac{d^{(1)}}{l(A1)}\right),$$
(26)
where $`\mu `$ is the genomic mutation rate,

$$\mu =l(1q).$$
(27)
Equation (26) follows from equation (24) if we disregard all terms of the order $`𝒪\left((1q)^2\right)`$ or higher. Figure 1a shows the exact solution for the average fitness $`f`$ equation (23) and the approximation equation (26) as a function of the genomic mutation rate $`\mu `$. For comparison, we have also displayed results from numerical simulations with a genetic algorithm. As was to be expected, the approximation works well for $`\mu 1`$, but breaks down for higher mutation rates. For $`\mu 1`$, the approximation significantly underestimates the average population fitness. Interestingly, this goes along with a decrease in the population neutrality (Fig. 1b). What happens is the following. In general, a population moves to that region of the genotype space where the probability of neutral offspring is maximized for the given mutation rate. Clearly, when the mutation rate is low, only the immediate neighbors influence that probability. If, however, the mutation rate is high, such that offspring with two or even more mutations become common, the immediate neighbors lose their importance. In the extreme case of a per symbol copy fidelity of $`q=0.5`$, the probability of giving birth to a viable sequence becomes identical for all sequences. Hence, in that extreme we can expect the population neutrality to coincide with the network neutrality $`\nu `$ (the network neutrality is the average number of neutral neighbors of a viable sequence), which amounts to $`\nu =\frac{2}{3}l`$ in our case. At the same time, most of the average population fitness does not stem from the one-mutants anymore, but from mutants further away, which explains why the approximation equation (26) underestimates the true average fitness.

## 4 The Neutral Staircase landscape

The fitness landscape that we have studied in the previous section contains only a single large neutral network, and hence epochal dynamics as predicted in Sec. 2.1 cannot be observed in that landscape. However, with a small modification, we can create a landscape that possesses the required properties. We subdivide a binary sequence into $`b`$ blocks of length $`2n`$, with a set of $`k`$ bits in between each block (Fig. 2). The total length of the sequence is thus $`l=2bn+k(b1)`$. The blocks can be *active* or *inactive*. An active block has properties similar to the sequences of the previous section. The block is subdivided into $`n`$ pairs. If any of those $`n`$ pairs are in the configuration 11, the fitness of the whole sequence is zero. Otherwise, the fitness is not affected by that block. The rightmost block is always active. A block further to the left becomes active if the block immediately to its right is active, and all $`k`$ bits between the two blocks are set to 1. Finally, any bit to the left of the leftmost active block that is set to 1 results in fitness zero for the whole sequence. We call this landscape the *Neutral Staircase* landscape, in analogy to the *Royal Staircase* introduced by van Nimwegen and Crutchfield . Note that the Neutral Staircase differs from the Royal Staircase in an important aspect: all sequences have either fitness 0 or fitness 1. No higher fitness genotypes can be discovered, and the population’s dynamics is determined only by the topology of the neutral sequences.

The analytical treatment of the Neutral Staircase landscape is straightforward for $`k>1`$. The neutral sequences decompose into $`b`$ neutral networks, one for each possible number of blocks that can be active. Within one neutral network, the average population fitness is readily available from (23). Namely, if $`i`$ counts the number of active blocks, we have

$$f_i=q^l(1+\frac{\stackrel{~}{\mu }^2}{2}+\frac{\stackrel{~}{\mu }}{2}\sqrt{8+\stackrel{~}{\mu }^2})^{ni}.$$
(28)
Likewise, it is simple to derive the probability with which a new block is discovered. First, we note that the probability with which an offspring sequence remains on the neutral network $`i`$ is given by $`\nu _i=f_i`$, according to (11). Hence, a single offspring sequence ends up on the next neutral network with probability $`f_i\stackrel{~}{\mu }^kq^{2n}(1q^2\stackrel{~}{\mu }^2)^n`$ (we need $`k`$ extra mutations to set the $`k`$ bits left of the leftmost active block to one, but we are allowed some miscopies in the newly activated block). In a finite population of $`N`$ sequences, there are on average $`Nf_i`$ sequences that can give rise to offspring. A sequence that belongs to the network $`i+1`$ is therefore created in one generation with probability

$$P_{\mathrm{crea},i}=1[1f_i\stackrel{~}{\mu }^kq^{2n}(1q^2\stackrel{~}{\mu }^2)^n]^{Nf_i}.$$
(29)
Next, we are interested in the probability of fixation of that newly discovered sequence, $`\pi _i`$. For large populations, $`\pi _i`$ can be approximated by $`2s`$ , where $`s`$ is the selective advantage of the newly discovered network. In our case, we obtain thus

$$\pi _i=2f_{i+1}/f_i1.$$
(30)
The probability of a transition from network $`i`$ to network $`i+1`$, $`P(ii+1)`$, is given by $`P(ii+1)=P_{\mathrm{crea},i}\pi _i`$, and hence the average epoch length $`\tau _i`$ follows as

$$\tau _i=\frac{1}{P(ii+1)}=\frac{1}{P_{\mathrm{crea},i}\pi _i}.$$
(31)
The time to convergence $`\tau _{\mathrm{conv}}`$, i.e., the time until the optimum sequence distribution has been found, is given by

$$\tau _{\mathrm{conv}}=\underset{k=1}{\overset{b1}{}}\tau _k.$$
(32)
Equations (28), (31), and (32) capture most of the immediately observable quantities in an evolving population in the Neutral Staircase landscape. We have tested their applicability with numerical simulations of a genetic algorithm, and have found good agreement. In Fig. 3, we compare the average fitness in the different epochs, $`f_i`$, to the simulation results. Figure 3a shows that the population fitness fluctuates around the predicted value during the metastable equilibrium, and that transitions between two equilibria happen fast. In Fig. 3b, we have made the entropic barrier between two networks smaller (we have reduced $`k`$). Then, all neutral sequences actually form only a single large neutral network, and the population transitions immediately to the highest possible level of the average fitness.

In Fig. 4, we compare the time to convergence $`\tau _{\mathrm{conv}}`$ with simulation results. We find that for $`\mu 3`$, our analytical prediction agrees well with the simulations, whereas for larger $`\mu `$, the prediction fails to capture the observed behavior. The origin of this failure is the following. We have derived equation (32) under the assumption that it is hard to discover the next higher network. Now, when $`\mu `$ becomes of the order of $`k`$ (we have $`k=3`$ in Fig. 4), mutations that flip all $`k`$ bits at once to activate the next block become likely, and hence this assumption fails. As in the case of Fig. 3b, the neutral sequences then form a single gigantic network, and an expression derived under the assumption of disjunct sub-networks must break down.

## 5 Discussion

The analysis presented in the preceding sections has shown that the topology of the neutral sequences in genotype space has an important influence on the dynamics of a population. It is therefore not justified to regard the evolution of a population of neutral genotypes as a simple diffusion process, unless either the population size or the mutation rate are very small (when the product of the population size, $`N`$ and the genomic mutation rate, $`\mu `$, is much smaller than one, the population as a whole performs essentially a random walk on the neutral network .) We have demonstrated that evolution on neutral networks can lead to the same kind of epochal dynamics that previously was thought to be caused solely by the discovery of genotypes of higher fitness. Of course, our results do not imply that sudden transitions in a population are never caused by such discoveries. Normally, a transition that we observe in a population evolving in an unknown fitness landscape will be due to the discovery of a faster replicating genotype. However, sometimes we may observe a transition to a higher average fitness without an increase (or even with a decrease, see below) in the fitness of the dominant or fastest replicating genotype. Such a transition is then due to the effects described in this paper.

In most cases, selection acts first and foremost on replication rates. However, if all viable genotypes are identical in terms of their replication rate, as is the case in the fitness landscapes we have studied in this paper, the next important quantity selection acts upon is the probability with which genotypes have viable offspring. This probability is more a property of a set of genotypes than of a single genotype, because a genotype that itself has high reproductive success, but produces mainly offspring with poor reproductive success, will ultimately have only a small number of progeny. Therefore, the selective pressure we have described here acts solely on clouds of mutants, in an extreme form of quasispecies-like selection .

Throughout this paper, we have assumed that all sequences have the same replication rate $`\sigma `$. It is natural to ask what happens when sequences with different replication rates are present. The simplest such situation occurs when all sequences within a single neutral network $`i`$ have the same replication rate $`\sigma _i`$. In that case, we can take the analysis given in Sec. 2.1 one step further. The equilibrium fitness $`f_i`$ within a network $`i`$ is then given by

$$f_i=\sigma _iq^l(1+\stackrel{~}{\mu }\rho _i),$$
(33)
where $`\rho _i`$ is the spectral radius of the connection matrix of the network, $`\text{G}_i`$. A population will of course try to move to the particular network that gives the maximum equilibrium fitness. For a small mutation rate $`\stackrel{~}{\mu }`$, it is clear that the network with the largest $`\sigma _i`$ yields the maximum equilibrium fitness. If, however, the mutation rate is large, such that $`\stackrel{~}{\mu }\rho _i`$ becomes of the order of unity or even exceeds that value, then depending on the distribution of the $`\rho _i`$’s and the $`\sigma _i`$’s over the different networks, a population may actually transition from a network $`j`$ with a higher replication rate $`\sigma _j`$ to a network $`i`$ with a smaller replication rate $`\sigma _i`$, while at the same time increasing its average fitness. A similar effect was already predicted by , although they investigated peaks with different support from slightly deleterious mutations, rather than neutral networks with different connection densities. In recent work, Ofria et al.  have observed this selection for neutrality empirically in digital organisms. In a large number of evolution experiments with different mutation rates, they observed that the number of neutral neighbors of the dominant genotype (a crude but practical measure for the spectral radius of the connection matrix) would increase at higher mutation rates, at the expense of the digital organisms’ replication rate. Similarly, Wilke et al.  showed that digital organisms with a vastly inferior replication rate could outcompete seemingly superior digital organisms at high mutation rates if the slower replicating organisms had a higher robustness against mutations.

It is interesting to relate our results for the average population fitness to the mutational load $`L`$, which is defined as

$$L=1\frac{f}{\sigma }.$$
(34)
From equation (9), we obtain

$$L=1q^l(1+\stackrel{~}{\mu }\rho ),$$
(35)
where $`\rho `$ is the spectral radius of the connection matrix G. As was already noted by van Nimwegen et al. , the load can deviate significantly from Haldane’s result $`L=\mu `$  if neutrality is present \[note that equation (35) becomes identical to Haldane’s result in the absence of neutrality ($`\rho =0`$) and in the limit of a small mutation rate\]. More importantly, as a generalization of Haldane’s result, it is often cited that the genetic load in an asexual population is independent of the fitness landscape, and therefore also of epistasis, and that it is equal to

$$L=1e^\mu .$$
(36)
(see e.g. Kondrashov , Charlesworth , derivation given by Kimura and Maruyama , Crow ). This result, however, holds only in the absence of neutral mutants. If neutral mutants are present, then the topology of the neutral genotypes in the genotype space is coupled to the mutational load, by virtue of $`\rho `$ in equation (35). Since the topology is also coupled to the type of epistasis that we observe , epistasis must have an influence on the mutational load, even in asexual populations, as long as neutral mutations can occur. However, equation (36) is the basis of the deterministic mutational hypothesis of the evolution of sex : if the mutational load is independent of epistasis for asexual populations, but depends strongly on the sign of epistasis for sexual populations, then for certain types of epistasis the mutational load of a sexual population may be much smaller than that of an asexual population. Because of the differences between equations (35) and (36), the deterministic mutational hypothesis could break down if neutrality were taken into account. As a consequence, this hypothesis should be reconsidered for the case of landscapes with neutrality.

## 6 Conclusions

We have shown that adaptive evolution can take place in the complete absence of what is ordinarily understood as advantageous genotypes. Even if the fitnesses of all viable genotypes are completely identical does selection favor particular regions of genotype space over others. What makes the difference is the density of neutral sequences. In a region of genotype space with a higher density of neutral sequences, chances are higher that a mutated offspring is neutral rather than deleterious. Therefore, neutral sequences in such a region have a higher robustness against mutation, and hence a higher reproductive success. This gives them sufficient selective advantage to outcompete sequences from a less densely connected region. The transitions between different such regions will often occur in sudden jumps, followed by relatively long periods of stasis. Evolution on neutral networks alone can thus lead to epochal dynamics observed in so many natural and artificial evolving systems.

## Acknowledgements

We thank Chris Adami for carefully reading this manuscript. This work was supported by the NSF under contract DEB-9981397.