text:
# Classical Limit of Demagnetization in a Field Gradient

## 1 Introduction

The problem of a system possessing magnetization in the presence of a magnetic field having a spatial gradient to its strength was considered classically by Carr and Purcell , Torrey , and Abragam . The rate of precession of the transverse component of a particle’s magnetic moment depends on where it is located in the field, but thermal motion and collisions smear this out. This was taken into account in  by carrying out a random walk, and in  and  by adding a diffusion term to the equation for the time rate of change of the magnetization. It leads to a reduction of its transverse component which is exponential in the cube of the time. Stejskal and Tanner  took the same approach but generalised the problem by permitting the gradient to be time dependent. In the approximate solution obtained there  the reduction of the magnetization has a more complicated time dependence. A classical treatment of NMR experiments involving spatial variation of the magnetization, including the effect of RF pulses, is given in Sodickson and Cory . We shall examine the diminution of transverse magnetization as a quantum mechanics problem to determine under what conditions the classical approach described above is justified.

An isolated spin 1/2 particle is described by a wave function the spin part of which can be written as a linear combination of two basis states, which can be taken to be eigenstates of the $`z`$-component of spin with eigenvalues $`\pm 1/2`$. Any given values of the relative magnitudes and phases of the coefficients of those two states represents a particle whose spin is pointing in some definite direction in 3-dimensional space. In the presence of a constant magnetic field the spin rotates coherently about the field.

The theoretical description of the process of interest is more complex and proceeds in two steps, precisely following the order in which the experimental conditions are created. Focussing on a single particle, it is located in a heat bath, which means there is a classical probability distribution for its momentum . In addition, the method of preparation (partially polarized) creates a probability distribution for its spin direction . At this first stage already the momentum and spin variables for the particle must be described by a density matrix . In the presence of a constant magnetic field the probability distribution of the particle’s spin rotates coherently about the field direction, just as for the case of an isolated particle in a pure state. At this stage the central interaction of the particle with the other particles, which are in the same heat bath, has no effect on the spin.

In the second step a gradient to the strength of the field is turned on. The momentum distribution of the particles becomes distorted away from thermal equilibrium because a magnetic moment in the direction of the field experiences a force in the direction of the gradient. Now the 2-body interaction does have an effect. The expectation value of the transverse component of the particle’s spin is determined by the off-diagonal elements of the density matrix in spin space, and we shall see that the rate at which those elements decrease with time is affected by the interplay between the two-body interaction and the field gradient.

The effect of (anti)symmetrization on transport properties was studied by Lhuillier and Laloë , and included initially in the present work. Those papers were not concerned with the effects of an applied magnetic field; and there appears to be a technical difference between their treatment of the effect of 2-body scattering on the 1-body density matrix and that given in Snider , which we follow. This is discussed in the text.

In the Theory section we formulate the problem in terms of the 1-body Wigner distribution functions, and the coupled integro-partial differential equations that they satisfy are obtained as well as the expressions for the expectation value of the particle’s spin in terms of those functions. Some important parameters are identified, including the momentum associated with the magnetic field gradient. Provided the gradient is sufficiently small and a number of other inequalities are satisfied, the equations are solved in Section 3 by treating the gradient of the magnetic field as a perturbation, to second order in its strength. The successive terms of the solution are gathered together in Section 4, where the result is compared with the classical solution. The approximations encountered along the way are summarized in Section 5.

A word about terminology. The bulk of the paper is concerned with the quantum mechanical calculation of the expectation value of the spin of a particle. The classical calculation with which it is compared speaks about magnetization, which is the magnetic moment per unit volume. The two quantities are simply related: multiplying the expectation value of the spin by the magnetic moment of a particle and by the number of particles per unit volume gives the magnetization since all the particles are identical.

## 2 Theory

The Hamiltonian for the system is taken to be a sum of 1-body terms, $`H^{(1)}`$, plus a sum of 2-body potentials $`V`$. $`H^{(1)}`$ contains the kinetic energy plus the interaction energy of the particle’s magnetic moment with the applied magnetic field, $`B_z`$, taken to be in the $`z`$-direction. For spin-1/2 particles this becomes

$$H^{(1)}=\frac{p^2}{2M}\frac{\gamma }{2}\mathrm{}\sigma _zB_z$$
(1)
where $`M`$ is the mass of the particles, $`\gamma `$ is the gyromagnetic ratio and $`\sigma _z`$ is the Pauli spin matrix. With the gradient of the magnetic field strength in the direction of the unit vector $`𝐮`$, $`B_z`$ can be written

$$B_z=B_0+𝐮𝐱G(t).$$
(2)
$`B_0`$ is the constant part of the field and $`G(t)`$ represents the magnitude of the gradient of the field.

We follow the development in Snider , but with the inclusion of (anti)symmetrization. There the density matrix for an N particle system is reduced by taking successive traces down to the 1-body density matrix, $`\rho `$. It satisfies the equation

$$\frac{\rho (1)}{t}=\frac{i}{\mathrm{}}\left([H^{(1)}(1),\rho (1)]+\underset{(2)}{Tr}[V(1,2),\rho ^{(2)}(1,2)]\right)$$
(3)
where the arguments of the operators refer to particles number 1 and 2 and $`\rho ^{(2)}`$ is the 2-body density matrix. (A superscript (1) is omitted on the 1-body $`\rho `$.) It is indicated in the second term within the parenthesis in (3) that the trace is to be taken over the variables of particle 2. The hierarchy of equations is cut off by neglecting the terms in the equation for $`\rho ^{(2)}/t`$ which involve $`\rho ^{(3)}`$, arguing that in a dilute gas when particles 1 and 2 are interacting the interaction with a third particle can be neglected . With $`n`$ being the number of particles per unit volume and $`a`$ the range of the two-body potential, this requires that $`ad=n^{1/3}`$. Normalising to the conditions found in the experiment reported by Schmidt et al.  involving a partially polarized gas of 3He atoms, which was performed at room temperature and a pressure of 7 atmospheres,

$$\frac{a}{d}=0.13\left(\frac{P}{7atm}\frac{293^{}}{T}\right)^{1/3}\left(\frac{a}{2.410^{10}m}\right)$$
(4)
Although the ratio $`a/d`$ may be on the borderline of being sufficiently small for that particular experiment , we shall neglect the possibility of triple collisions, and will continue to take $`a=2.410^{10}`$ m as an estimate of the range of the He-He potential when examining the size of other relevant parameters.

We assume that prior to a collision between particles 1 and 2 they are uncorrelated, and hence $`\rho ^{(2)}`$ is just the product of the two 1-body density matrices; this corresponds to the Boltzmann property for the classical distribution functions. Using scattering theory the effect of a collision can be obtained as

$$\rho ^{(2)}(1,2)=\mathrm{\Omega }^{(+)}\rho (1)\rho (2)\mathrm{\Omega }^{(+)}$$
(5)
where $`\mathrm{\Omega }^{(+)}`$ is the wave operator . When this expression for $`\rho ^{(2)}`$ is inserted into (3) it becomes an equation for the 1-body density matrix $`\rho `$ alone. We shall seek a solution of this equation to second order in the strength of the gradient $`G(t)`$. In Lhuillier and Laloë  instead of the wave matrix $`\mathrm{\Omega }^{(+)}`$ in Eq.(5) the $`S`$-matrix is used. This relates the density matrix in the distant past, before the collision, to its value in the distant future after the collision is completed. But we are interested in the time development of $`\rho `$ during a collision  and choose the wave matrix, therefore. Our numerical results would be different had we chosen the $`S`$-matrix.

In the experiment reported by Schmidt et al. , a gas of 3He atoms is partially polarized in the $`z`$-direction and fairly uniform spatially. At $`t=0`$ it is tipped through an angle $`\theta `$ thereby producing a component of spin in the plane transverse to the direction of the magnetic field as well as a component in the $`z`$-direction. It is the diminution of the transverse component with time that was studied, and which we wish to evaluate. The most general form of the density matrix for spin-1/2 particles can be written

$`\rho `$ $`=`$ $`\rho _II+𝝆𝝈`$ (6)
$`=`$ $`\rho _II+\rho _z\sigma _z+{\displaystyle \frac{1}{2}}(\rho _+\sigma _{}+\rho _{}\sigma _+)`$
where for any operator, $`A_\pm =A_x\pm iA_y`$, and the $`\sigma _i`$ are the Pauli spin matrices. The four operators $`\rho _I`$, $`\rho _z`$, $`\rho _+`$, and $`\rho _{}=\rho _+^{}`$ operate only in ordinary (position or momentum) space.

It is convenient to work with the Wigner distribution function  generalized to a 2x2 matrix to include the spin degree of freedom. This generalization was also made in reference  . It is defined by

$`f_{mm^{}}(𝐱,𝐩,t)`$ $`=`$ $`{\displaystyle \frac{1}{(\mathrm{}\pi )^3}}{\displaystyle \frac{d^3q}{(2\pi \mathrm{})^3}<𝐩𝐪,m|\rho |𝐩+𝐪,m^{}>\mathrm{exp}(2i𝐱𝐪/\mathrm{})}`$ (7)
$`=`$ $`{\displaystyle \frac{1}{(\mathrm{}\pi )^3}}{\displaystyle d^3y<𝐱𝐲,m|\rho |𝐱+𝐲,m^{}>\mathrm{exp}(2i𝐩𝐲/\mathrm{})}`$ (8)
where the indices $`m`$ and $`m^{}`$ take the values $`\pm 1/2`$. It follows from the fact that $`\rho `$ is hermitean that the 2x2 matrix $`f`$ is also. We also need the inverse relation

$$<𝐩,m|\rho |𝐩^{},m^{}>=(2\pi \mathrm{})^3d^3xf_{mm^{}}(𝐱,(𝐩+𝐩^{})/2,t)\mathrm{exp}[i(𝐩^{}𝐩)𝐱/\mathrm{}].$$
(9)
Taking the partial derivative of Eq.(7) or (8) with respect to time and using Eq.(3) gives

$$\frac{f_{mm^{}}}{t}=S_{mm^{}}^K+S_{mm^{}}^\gamma +S_{mm^{}}^V$$
(10)
where $`S^K`$ arises from the kinetic energy term in $`H^{(1)}`$, $`S^\gamma `$ from the magnetic energy term in $`H^{(1)}`$, and $`S^V`$ from the 2-body potential energy term. $`S^K`$ is the well known flow term

$$S_{mm^{}}^K=\frac{1}{M}𝐩\mathbf{}f_{mm^{}}.$$
(11)
To obtain $`S^\gamma `$ requires evaluation of the commutator $`[\sigma _zB_z,\rho ]`$. Making use of Eqs. (2) and (6) yields

$`S_{1/2}^\gamma `$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}\mathrm{}G(t)𝐮\mathbf{}_pf_{1/2}`$ (12)
$`S_{1/2}^\gamma `$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}\mathrm{}G(t)𝐮\mathbf{}_pf_{1/2}`$ (13)
$`S_+^\gamma `$ $`=`$ $`i\gamma B_zf_+`$ (14)
where we have introduced abbreviated notations for all the $`S`$ terms as follows: $`S_{1/2}S_{1/\mathrm{2\; 1}/2}`$ ; $`S_{1/2}S_{1/21/2}`$ ; and $`S_+S_{1/\mathrm{2\; 1}/2}`$; and the same definitions apply to the components of the Wigner distribution functions $`f`$. We shall be particularly interested in the function $`f_+`$, which involves the off-diagonal elements of $`\rho `$ in the 2-dimensional spin space. In conjunction with Eq.(10) for the time rate of change of the distribution functions, the interpretation of Eqs. (12) and (13) is as follows. A magnetic moment having a component in the direction of (or opposite to) a magnetic field whose strength has a spatial gradient will experience a force in the direction of (or opposite to) that gradient. Eq.(14) says that a magnetic moment transverse to the direction of the magnetic field undergoes Larmor precession about the field.

With the density matrix normalized to Tr($`\rho )=N`$, the total number of particles, it is straightforward to obtain the following results for $`n`$, the number of particles per unit volume, and $`<𝝈>`$, the expectation value of the particle’s spin at position $`𝐱`$ and time t

$$n=d^3p[f_{1/2}(𝐱,𝐩)+f_{1/2}(𝐱,𝐩)]$$
(15)
and

$`<\sigma _z>`$ $`=`$ $`{\displaystyle \frac{d^3p[f_{1/2}(𝐱,𝐩)f_{1/2}(𝐱,𝐩)]}{d^3p[f_{1/2}(𝐱,𝐩)+f_{1/2}(𝐱,𝐩)]}}`$ (16)
$`<\sigma _+>`$ $`=`$ $`{\displaystyle \frac{2d^3pf_+(𝐱,𝐩,t)}{d^3p[f_{1/2}(𝐱,𝐩)+f_{1/2}(𝐱,𝐩)]}}.`$ (17)
Once these quantities have been obtained the net magnetic moment per unit volume, $`𝐌`$, can be obtained simply by multiplying $`<𝝈>`$ by the magnetic moment of a single particle, $`\gamma \mathrm{}/2`$, and the number of particles per unit volume, $`n`$.

We now turn to $`S^V`$, the collision term in Eq.(10). Introducing the spin and momentum variables for particle 2, over which the trace indicated in Eq.(3) will be taken, and the antisymmetrization operator ($`I+ϵP_{12}`$) where $`P_{12}`$ permutes all the variables of particles 1 and 2 and $`ϵ=1`$ for fermions, $`S^V`$ can be written

$`S_{mm^{}}^V(𝐱,𝐩)={\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{1}{(\mathrm{}\pi )^3}}{\displaystyle \underset{m_2}{}}{\displaystyle \frac{d^3q}{(2\pi \mathrm{})^3}\frac{d^3p_2}{(2\pi \mathrm{})^3}}`$ (18)
$`<𝐩𝐪,m;𝐩_\mathrm{𝟐},m_2|(I+ϵP_{12})[V,\rho ^{(2)}(1,2)]|𝐩+𝐪,m^{};𝐩_2,m_2>`$
$`\mathrm{exp}(2i𝐱𝐪/\mathrm{}).`$
It is sufficient to antisymmetrize just one of the states because $`V`$ and $`\rho ^{(2)}`$ are symmetric under the interchange of the two particles. In the Appendix we will develop the collision term making two approximations. The first is to assume that the 2-particle potential energy $`V`$ consists purely of a central interaction. For example, at a separation of 2.4 Å the moment-moment interaction energy of two 3He atoms is $`10^8`$ mev, whereas the central interaction is $`1`$ mev. The second approximation amounts to assuming that the gradient of the magnetic field strength is sufficiently small; this comes about as follows.

Examination of Eqs.(10) through (14), which represent the contributions to the time derivative of the respective distribution functions, shows that only $`S_+^\gamma `$ contains dependence on position via the gradient term in the magnetic field, $`B_z`$. Neglecting the collision terms $`S^V`$ for the moment, would lead to the following conclusions. With the particle density $`n`$ and the expectation value of the spin initially independent of position, $`n`$ and $`<\sigma _z>`$ remain so since $`f_{\pm 1/2}`$ are independent of position. $`n`$ is, of course, also independent of time, and $`<\sigma _z>`$ as well since $`\sigma _z`$ commutes with the Hamiltonian. On the other hand $`f_+`$ would contain a factor $`\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))]`$ as the only position dependence in the problem, where

$$F(t)=_0^t𝑑t^{}G(t^{}).$$
(19)
We show below that all these conclusions remain valid with collisions and write, therefore,

$$f_+(𝐱,𝐩,t)=g_+(𝐩)\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))].$$
(20)
The product $`\gamma F(t)`$ is the wave number associated with the magnetic field gradient ; in that reference it is designated as $`k`$.

The consequences for the matrix elements of $`\rho `$ can be seen as follows. From Eq.(9) those matrix elements with $`m=m^{}`$, since they arise from $`f_{\pm 1/2}`$ which do not depend on position, are local in momentum space,

$$<𝐩,m|\rho |𝐩^{},m>=(2\pi \mathrm{})^3f_m(𝐩)(2\pi \mathrm{})^3\delta ^{(3)}(𝐩^{}𝐩).$$
(21)
On the other hand, for the off-diagonal matrix elements we have

$`<𝐩,1/2|\rho |𝐩^{},1/2>=(2\pi \mathrm{})^3g_+((𝐩+𝐩^{})/2)\mathrm{exp}(i\gamma B_0t)`$ (22)
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩^{}𝐩\gamma \mathrm{}𝐮F(t)).`$
Position dependence of a distribution function is sometimes neglected in the collision term on the assumption that the quantity that is varying (the strength of the magnetic field in the present case) changes little over the distance of a mean free path . But the $`\delta `$-function in Eq.(22) contains another relevant quantity, $`\gamma \mathrm{}𝐮F(t)`$, that has nothing to do with the mean free path. It is a momentum that arises from the oppositely directed forces felt by particles with their moments in the direction of the field, or opposite to it, due to the field gradient. After inserting two complete sets of 2-body spin and momentum states into Eq.(18) matrix elements of the type $`<𝐩,1/2|\rho |𝐩^{},1/2>`$ appear linearly in the collision term $`S_+^V`$ and bilinearly (or not at all) in $`S_{\pm 1/2}^V`$. When the $`\delta `$-function from Eq.(22) is combined with the two momentum conserving $`\delta `$-functions from the $`t`$\- and $`\mathrm{\Omega }^{(+)}`$-matrix elements \[see (A-1) and (A-2)\] it has the following effects. In Eq.(18) for $`S_{\pm 1/2}^V`$ the integration variable $`𝐪`$ is required to be zero; but in $`S_+^V`$ it becomes $`\gamma \mathrm{}F(t)𝐮/2`$. Consequently $`S_+^V`$ acquires a factor $`\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))]`$, and as a result every term in Eq.(10) for $`f_+/t`$ has that factor. This confirms the statements above that $`f_{\pm 1/2}`$ are independent of position and $`f_+`$ has the structure shown in Eq.(20). No assumption about the magnitude of $`F(t)`$ has been made at this point.

If the magnitude of $`\gamma \mathrm{}F(t)`$ were not sufficiently small compared to the important momenta in the problem there would be significant quantum mechanical consequences for the collision terms of the displaced argument of the $`\delta `$-function in Eq.(22). The arguments of $`g_+`$ and the arguments of the $`t`$-matrix elements would be shifted by amounts of order $`\gamma \mathrm{}F(t)𝐮`$. We shall delay discussion of the magnitude of the effect on the collision terms due to the altered $`t`$-matrix elements, but note here that it would even lead to the presence of non-energy-conserving matrix elements.

The effect of the variation in $`g_+`$ can be seen directly. Although the distribution functions are of the Fermi-Dirac type, we shall be interested in conditions of negligible degeneracy, $`n\overline{\lambda }^31`$, where $`\overline{\lambda }`$ is the mean wavelength of the particles. In that case we can use the Maxwell-Boltzmann distribution, and its momentum derivative is of order $`(p/\overline{p}^2)`$ times the function itself, where $`\overline{p}`$ is the mean momentum at temperature $`T`$. The fractional shift in the magnitude of $`g_+`$, therefore, is of order $`\gamma \mathrm{}F(t)/\overline{p}`$. We shall estimate it by normalising to the experimental conditions found in reference . The gyromagnetic ratio for 3He is 204 MHz/Tesla, and the maximum value of $`F(t)`$ in that experiment, designated $`F_M`$, is 3.5 $`10^7`$ Tesla-sec/cm. Estimating $`\overline{p}`$ by the mean value of the momentum at temperature $`T`$, we obtain

$$\frac{\gamma \mathrm{}F(t)}{\overline{p}}9.510^8\left(\frac{\gamma }{\gamma _{{}_{}{}^{3}\mathrm{He}}}\right)\left(\frac{F(t)}{F_M}\right)\left(\frac{M_{{}_{}{}^{3}\mathrm{He}}}{M}\frac{293^{}}{T}\right)^{1/2}.$$
(23)
For any other set of experimental parameters it is only necessary to insert the appropriate values of the gyromagnetic ratio $`\gamma `$, the temperature $`T`$, the molecular mass $`M`$, and the integrated strength of the gradient of the magnetic field, $`F(t)`$. Note that a larger strength of the field gradient or a smaller temperature would increase the ratio in Eq.(23), but the variation with $`T`$ is slow.

For the remainder of the paper we shall assume that the shift in the arguments of $`g_+`$ and the $`t`$-matrix elements can be neglected, and this leads to the expressions for the collision terms presented in Eqs.(A-17) and (A-18). Eqs.(10) then become the following set of coupled nonlinear integro- partial differential equations.

$$\frac{f_{1/2}(𝐩)}{t}=\frac{\gamma }{2}\mathrm{}G(t)𝐮\mathbf{}_pf_{1/2}(𝐩)+S_{1/2}^V(𝐩)$$
(24)
$$\frac{f_{1/2}(𝐩)}{t}=\frac{\gamma }{2}\mathrm{}G(t)𝐮\mathbf{}_pf_{1/2}(𝐩)+S_{1/2}^V(𝐩)$$
(25)
and in the equation for $`f_+/t`$ the derivative of the exponential with respect to time cancels the term $`S_+^\gamma `$ leaving

$$\frac{g_+(𝐩)}{t}=i\frac{\gamma }{M}𝐮𝐩F(t)g_+(𝐩)+\overline{S}_+^V(𝐩).$$
(26)
## 3 Perturbation Expansion

Our approach for solving these equations, consistent with the approximation described above, is to treat the strength of the field gradient, $`G(t)`$, as a perturbation. Since it is multiplied by the gyromagnetic ratio $`\gamma `$ we shall expand the distribution functions up to second order in $`\gamma `$,

$$f_{\pm 1/2}=f_{\pm 1/2}^{(0)}+\gamma f_{\pm 1/2}^{(1)}+\gamma ^2f_{\pm 1/2}^{(2)}+O(\gamma ^3)$$
(27)
and

$$g_+=g_+^{(0)}+\gamma g_+^{(1)}+\gamma ^2g_+^{(2)}+O(\gamma ^3).$$
(28)
This procedure converts the collision terms $`S^V`$ from being bilinear in the functions $`f`$ and $`g`$ to being linear in the successive approximations $`f^{(i)}`$ and $`g^{(i)}`$.

### Zeroth Order

The terms of order $`\gamma ^0`$ in Eqs.(24) - (26) represent the initial condition in which the probability of finding the particle with a given momentum is given by a thermal distribution, and the expectation value of its spin is specified; both quantities are initially independent of position. Those distribution functions can be written

$$f_{\pm 1/2}^{(0)}(𝐩)=n_{\pm 1/2}B^{(M)}(p)$$
(29)
and

$$g_+^{(0)}(𝐩)=n_+B^{(M)}(p),$$
(30)
where

$$B^{(M)}(p)=(\beta ^{(M)}/\pi )^{3/2}\mathrm{exp}(\beta ^{(M)}p^2)$$
(31)
is the the Maxwell-Boltzmann distribution for a particle of mass $`M`$, and $`\beta ^{(M)}=(1/2MkT)`$. As expected from conservation of energy , when these zeroth order distribution functions are inserted in the collision integrals they can be shown to vanish, confirming that Eqs.(24) - (26) are satisfied to this order. The unitarity relation gets used here and in the subsequent development \[see Appendix\]. The significance of the normalization constants can be seen from Eqs.(15) to (17), since they determine the initial values of the number density, $`n`$, and the expectation value of the spin, $`<𝝈>`$: $`n=n_{1/2}+n_{1/2},<\sigma _z>=(n_{1/2}n_{1/2})/n`$, and $`<\sigma _+>=n_+/n`$.

### First Order

At the next order in $`\gamma `$ the distribution functions are distorted in the direction $`𝐮`$ of the gradient, so we write

$$f_{\pm 1/2}^{(1)}(𝐩)=n_{\pm 1/2}B^{(M)}(p)𝐮𝐩h_{\pm 1/2}^{(1)}$$
(32)
and

$$g_+^{(1)}(𝐩)=n_+B^{(M)}(p)𝐮𝐩h_+^{(1)}.$$
(33)
In principle the functions $`h^{(1)}`$ could depend on the magnitude of the momentum, and a systematic expansion in a complete set of Sonine polynomials, for example , could be undertaken. We shall take the $`h^{(1)}`$ functions to be independent of $`p`$, thereby limiting that expansion to just the lowest term. Note that these first order terms in the distribution functions contribute nothing to the number density and the expectation value of the spin, Eqs. (15) - (17), because their integrals over all momenta vanish. As expected if $`<\sigma _z>0`$ the particles acquire a non-zero expectation value of momentum, and it can be shown to be proportional to the product of $`(n_{1/2}n_{1/2})/n`$ and $`\gamma \mathrm{}𝐮F(t)`$. The latter factor is the same one we required to be small compared to thermal momenta in the discussion above.

Now pick out of Eqs. (24) - (26) all the terms of first order in $`\gamma `$. Since the functions $`h^{(1)}`$ depend only on the time, these become coupled ordinary first order linear differential equations

$$n_{1/2}B^{(M)}(p)𝐮𝐩\frac{dh_{1/2}^{(1)}}{dt}=n_{1/2}\mathrm{}\beta ^{(M)}B^{(M)}(p)𝐮𝐩G(t)+S_{1/2}^{V(1)}$$
(34)
$$n_{1/2}B^{(M)}(p)𝐮𝐩\frac{dh_{1/2}^{(1)}}{dt}=n_{1/2}\mathrm{}\beta ^{(M)}B^{(M)}(p)𝐮𝐩G(t)+S_{1/2}^{V(1)}$$
(35)
and

$$n_+B^{(M)}(p)𝐮𝐩\frac{dh_+^{(1)}}{dt}=i\frac{1}{M}n_+B^{(M)}(p)𝐮𝐩F(t)+\overline{S}_+^{V(1)}.$$
(36)
The expressions for the first order collision terms $`S^{V(1)}`$ are given in the Appendix, where the above equations are developed into Eqs.(A-24) to (A-26).

There are three quantities in those equations, $`I_U,I_I,`$ and $`I_\pi `$ that depend upon the details of the 2-body potential energy; see Eq.(A-23). $`I_U`$, which is the average over a Maxwell-Boltzmann distribution of the cube of the momentum multiplied by the transport cross section, enters into the direct terms. $`I_I`$ and $`I_\pi `$, which involve the $`t`$-matrix elements themselves and not just their absolute square, appear only in the exchange terms, i.e., those proportional to $`ϵ`$, which arise from (anti)symmetrization of the wave functions. To study the relative importance of the latter it is necessary to specify the 2-body potential energy function for the problem of interest and calculate the $`t`$-matrix that it produces. For He the interatomic potential becomes strongly repulsive at distances less than 2.4 Å , so we shall use a hard sphere potential as a guide; and there are two limiting cases for which analytic approximations to the scattering amplitude can be obtained: $`ka/\mathrm{}1`$ and $`ka/\mathrm{}1`$ where k is the relative momentum of the colliding particles and $`a`$ is the radius of the hard sphere potential. Numerical calculations of the quantum mechanical scattering amplitude in the short wavelength limit  are in qualitative and even semi-quantitative agreement with the classical scattering of a scalar wave on a rigid sphere, as presented in Morse and Feshbach , and also by the scattering of an electromagnetic wave by a perfectly conducting sphere, given in Jackson . That $`t`$-matrix takes the approximate form

$$t(k,\theta )\frac{2\pi \mathrm{}^2}{\mu }\frac{a}{2}\left(i\frac{ka}{\mathrm{}}\frac{(1+\mathrm{cos}\theta )J_1(\frac{ka}{\mathrm{}}\mathrm{sin}\theta )}{\frac{ka}{\mathrm{}}\mathrm{sin}\theta }+\mathrm{exp}[2i\frac{ka}{\mathrm{}}\mathrm{sin}\frac{\theta }{2}]\right)(ka/\mathrm{}1)$$
(37)
where the first term arises from scattering by the shadow side and the second term from the illuminated side; $`\mu `$ is the reduced mass. The pure imaginary shadow scattering amplitude dominates in the small angle region $`\theta \stackrel{<}{}\mathrm{}/ka`$, and oscillates rapidly with decreasing amplitude for larger angles.

Normalizing once again to the conditions in reference  gives

$$\frac{\overline{k}a}{\mathrm{}}=12.5\left(\frac{M}{M_{{}_{}{}^{3}\mathrm{He}}}\frac{T}{293^{}}\right)^{1/2}\frac{a}{2.4\mathrm{\AA }}$$
(38)
where we have chosen 2.4 Å as the hard sphere radius for a gas of $`{}_{}{}^{3}\mathrm{He}`$ atoms. Although it would require a very low temperature to get into the long wavelength limit, even a moderately low temperature, $`10^{}`$, would introduce significant corrections to the $`t`$-matrix as well as the first order collision terms, $`S^{V(1)}`$. In particular, contributions to $`S^{V(1)}`$ arising from antisymmetrization of the wave functions would not be negligible and this would produce important modifications to the remainder of the paper. Keeping that in mind we now estimate the magnitudes of $`I_U,I_I`$ and $`I_\pi `$ in Eq.(A-23) using the short wavelength expression (37). It should be noted that this limit, $`\overline{\lambda }/a1`$, together with the low density assumption, $`na^31`$, imply $`n\overline{\lambda }^31`$, which is the negligible degeneracy assumption made earlier. Using Eq.(37) we can also verify the assumption made earlier that the shift in the momentum (or momentum transfer) argument of the $`t`$-matrix element by the amount $`\gamma \mathrm{}F(t)`$ has a negligible effect for the conditions of reference .

Although the shadow scattering and the scattering from the illuminated side each contribute $`\pi a^2`$ to the total cross section, the factor ($`1\mathrm{cos}\theta `$) in Eq.(A-22) for the transport cross section eliminates the shadow contribution, hence $`\sigma _U(k)=\pi a^2`$; and from Eq.(A-23) $`I_U=\pi \overline{k^3}a^2`$, where $`\overline{k^3}`$ is the average of $`k^3`$ over a Maxwell-Boltzmann distribution for a particle having mass $`\mu `$. It takes a considerable amount of algebra to show that when the short wavelength expression Eq.(37) for the scattering amplitude is inserted into Eqs.(A-13), (A-22) and (A-23), $`I_I`$ is smaller than $`I_U`$ by $`O[(\overline{k}a/\mathrm{})^2]`$. $`I_\pi `$ is even smaller, being $`O(\mathrm{exp}[(\overline{k}a/\mathrm{})^2])`$. Therefore all the terms proportional to $`ϵ`$ in Eqs.(A-24)-(A-26) can be neglected.

In the short wavelength limit Eqs.(A-24)-(A-26) for the time derivative of the first order functions $`h^{(1)}`$ become

$$n_{1/2}\frac{dh_{1/2}^{(1)}}{dt}=n_{1/2}\mathrm{}\beta ^{(M)}G(t)+\frac{4\beta ^{(M)}}{3M}n_{1/2}n_{1/2}[h_{1/2}^{(1)}h_{1/2}^{(1)}]I_U$$
(39)
$$n_{1/2}\frac{dh_{1/2}^{(1)}}{dt}=n_{1/2}\mathrm{}\beta ^{(M)}G(t)\frac{4\beta ^{(M)}}{3M}n_{1/2}n_{1/2}[h_{1/2}^{(1)}h_{1/2}^{(1)}]I_U$$
(40)
and

$$\frac{dh_+^{(1)}}{dt}=i\frac{1}{M}F(t)+\frac{4\beta ^{(M)}}{3M}[n_{1/2}h_{1/2}^{(1)}+n_{1/2}h_{1/2}^{(1)}(n_{1/2}+n_{1/2})h_+^{(1)}]I_U.$$
(41)
Summing the first two equations gives

$$n_{1/2}h_{1/2}^{(1)}+n_{1/2}h_{1/2}^{(1)}=(n_{1/2}n_{1/2})\mathrm{}\beta ^{(M)}F(t)$$
(42)
where $`F(t)`$ is the time integral of the field gradient, Eq.(19). There is no need to solve for $`h_{1/2}^{(1)}`$ and $`h_{1/2}^{(1)}`$ separately because only the combination shown in Eq.(42) is needed in the equation for $`dh_+^{(1)}/dt`$; and furthermore the terms of first order in $`\gamma `$ do not contribute to the magnetization.

The solutions for the real and imaginary parts of $`h_+^{(1)}`$ are

$$\mathrm{Re}h_+^{(1)}=\frac{4\beta ^{(M)}}{3M}(n_{1/2}n_{1/2})\mathrm{}\beta ^{(M)}I_U\mathrm{exp}(\alpha t)_0^t𝑑t^{}F(t^{})\mathrm{exp}(\alpha t^{})$$
(43)
and

$$\mathrm{Im}h_+^{(1)}=\frac{1}{M}\mathrm{exp}(\alpha t)_0^t𝑑t^{}F(t^{})\mathrm{exp}(\alpha t^{})$$
(44)
where

$$\alpha =\frac{4\beta ^{(M)}}{3M}nI_U.$$
(45)
$`n=n_{1/2}+n_{1/2}`$ is the total number of particles per unit volume.

$`\alpha `$ can also be expressed in terms of the diffusion constant, $`D`$, for a related problem. Suppose there is no magnetic field and no transverse component of the magnetization, but there is a spatial gradient to the $`z`$-component of magnetization. This reduces to a standard 2-component problem with equal masses (spin up and spin down); and if the diffusion constant is calculated with the same approximation as above, namely keeping only the leading term in an expansion in Sonine polynomials, then

$$\alpha =\frac{kT}{MD}.$$
(46)
The qualitative significance of $`\alpha ^1`$ can be seen from Eq.(45) or its expression in terms of the diffusion constant: it is the mean free time between collisions. Normalising to the conditions in reference

$$\alpha =3.510^{10}\mathrm{sec}^1\left(\frac{P}{7\mathrm{a}\mathrm{t}\mathrm{m}}\right)\left(\frac{M_{{}_{}{}^{3}\mathrm{He}}}{M}\frac{293^{}}{T}\right)^{1/2}\left(\frac{a}{2.410^{10}\mathrm{m}}\right)^2.$$
(47)
Provided $`G(t)`$, and therefore $`F(t)`$, vary on a time scale much longer than $`\alpha ^1`$, the integrals in Eqs.(43) and (44) are well approximated for $`t>>\alpha ^1`$ by

$$\mathrm{exp}(\alpha t)_0^t𝑑t^{}F(t^{})\mathrm{exp}(\alpha t^{})\frac{F(t)}{\alpha }.$$
(48)
This leads to

$$\mathrm{Re}h_+^{(1)}=\frac{n_{1/2}n_{1/2}}{n}\mathrm{}\beta ^{(M)}F(t)$$
(49)
and

$$\mathrm{Im}h_+^{(1)}=\frac{1}{M\alpha }F(t)$$
(50)
For conditions similar to those of reference  $`\mathrm{Re}h_+^{(1)}`$ is much smaller than $`\mathrm{Im}h_+^{(1)}`$, their ratio being, at most, that of the mean wavelength to the mean free path; this will be small if the previous assumptions that $`\overline{\lambda }a`$ and $`ad`$ are both satisfied, and we shall neglect Re$`h_+^{(1)}`$. It vanishes exactly if there is no $`z`$-component to the expectation value of the spin.

### Second Order

There is a restriction on the form of a distortion of the distribution functions $`f_{\pm 1/2}`$ that is of second order in $`\gamma `$ (and therefore also of second order in the direction vector $`𝐮)`$, which arises from the requirement that it not alter the values of the particle number density, $`n`$, and the expectation value of the $`z`$-component of spin, $`<\sigma _z>`$, both of which are fixed by the terms of zero order. Since these quantities are obtained from Eqs.(15) and (16), the integral over all momenta of the distortions must vanish. The only allowed form is

$$f_{\pm 1/2}^{(2)}(𝐩)=n_{\pm 1/2}B^{(M)}(p)[𝐮𝐮2\beta ^{(M)}(𝐮𝐩)^2]h_{\pm 1/2}^{(2)}.$$
(51)
Indeed, such a form arises naturally from a second application of the operator $`𝐮_p`$ from Eq.(24) on $`f_{\pm 1/2}^{(0)}`$. For $`g_+^{(2)}`$ we write

$$g_+^{(2)}=n_+B^{(M)}(p)(𝐮𝐩)^2h_+^{(2)}(t),$$
(52)
which arises similarly from a second application of $`𝐮𝐩`$ from Eq.(26) on $`g_+^{(0)}`$ .

Now pick out of Eqs.(24)-(26) for the time derivatives of the distribution functions the terms of second order in $`\gamma `$. When the equations for $`f_{\pm 1/2}^{(2)}/t`$ are integrated over all $`𝐩`$ they become 0=0. The equation for the time derivative of $`g_+^{(2)}`$ becomes

$$n_+B^{(M)}(p)(𝐮𝐩)^2\frac{dh_+^{(2)}}{dt}=i\frac{1}{M}n_+B^{(M)}(p)(𝐮𝐩)^2F(t)h_+^{(1)}+\overline{S}_+^V(𝐩)$$
(53)
Integrating this equation over all $`𝐩`$ and changing the variables of integration in the collision term to ($`𝐏,𝐤`$) that term can be shown to contribute nothing; hence

$$h_+^{(2)}(t)=i\frac{1}{M}_0^t𝑑tF(t)h_+^{(1)}(t).$$
(54)
Inserting the expression for Im$`h_+^{(1)}`$ from Eq.(50) gives

$$h_+^{(2)}(t)=\frac{1}{M^2\alpha }_0^t𝑑tF^2(t).$$
(55)
## 4 Time Dependence of $`<𝝈>`$

We now gather together all the contributions to $`<𝝈>`$, the expectation value of the particle’s spin, up to second order in the strength of the gradient of the magnetic field, by inserting into equations (16) and (17) the results from (20), (27) -(31) , (52) and (55).

$$<\sigma _z>=<\sigma _z>(t=0)$$
(56)
and

$`<\sigma _+>=<\sigma _+>(t=0)\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))]`$ (57)
$`\left[1\gamma ^2D{\displaystyle _0^t}𝑑tF^2(t)+O(\gamma ^4)\right],`$
where we have made use of the expression for the diffusion constant $`D`$ from Eq.(46). Eq.(57) represents the rate of decrease of the expectation value of the transverse component of the particle’s spin to second order in the strength of the gradient of the magnetic field.

How does this solution for $`<𝝈>`$ compare with the procedure followed in references - mentioned at the beginning of the paper? This consists of adding a classical diffusion term to the equation for the time rate of change of the magnetization, i.e.,

$$\frac{𝐌}{t}=\gamma (𝐌\times 𝐁)+D^2𝐌.$$
(58)
Recall that $`𝐌`$ is given by $`<𝝈>`$ multiplied by the magnetic moment of a particle, $`\gamma \mathrm{}/2`$, and the number of particles per unit volume. With $`𝐁`$ given in Eq.(2) and the initial value of $`𝐌`$ independent of postion, the exact solution of Eq.(58) is that $`M_z`$ remains at its initial value; and $`M_+`$ is given by

$$M_+(𝐱,t)=M_+(t=0)\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))]\mathrm{exp}\left[\gamma ^2D_0^tF^2(t)𝑑t\right].$$
(59)
The final exponential factor is found in references  and . If that exponential is expanded for small argument, to second order in $`\gamma `$ it gives the corresponding term from Eq.(57).

It was stated in the Introduction that the rate of decrease of the transverse component of magnetization is affected by both thermal motion and the interaction of a particle with its environment. By turning off the 2-body potential in Eq.(26), thereby setting the collision terms $`S^V`$ to zero, the rate due just to the random velocity distribution can be found, and the exact result is

$`M_+(𝐱,t)=M_+(t=0)\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))]\mathrm{exp}\left[\gamma ^2{\displaystyle \frac{kT}{M}}\left({\displaystyle _0^t}F(t)𝑑t\right)^2\right]`$ (60)
(no collisions)
Using expression (46) for the diffusion constant shows that the ratio of the arguments of the exponentials in Eqs.(60) and (59) is qualitatively given by $`\alpha T`$ where $`\alpha ^1`$ is the mean free time between collisions and $`T`$ is the duration of the gradient. For the conditions of the experiment in reference  $`\alpha T10^7`$. It is not surprising that the transverse component of magnetization decreases more rapidly without collisions because the particles sample a much larger range of magnetic field strengths in a given time than they would if experiencing collisions, and consequently their precession rates cover a larger range.

For the conditions of the experiment in reference  the quantity $`\gamma ^2_0^TF^2(t)𝑑t`$ has the value 0.032 sec/mm2, where $`T`$ is the time at which the gradient is turned off. Using the experimental values  of the ratio $`|M_+(T)/M_+(0)|`$ of the magnitudes of the transverse magnetization (at various points along the length of the cylinder) after the gradient is turned off to its initial value, and comparing with Eq.(59), the value of the 3He self diffusion coefficient is found to be $`D`$=15.9 mm2/sec . The discrepancy with the values shown in Fig.5a of reference  is due to a numerical error in that reference . Although the statistical error in the value of $`D`$ is small, only 0.4 mm2/sec, there is a much larger systematic error of almost 25 percent because the cell was prepared initially at much higher than room temperature .

We also evaluated $`D`$ from its relation to $`\alpha `$, Eq.(46), using the estimate given in (47) based on a hard sphere potential with radius $`a`$=2.4 Å, and find $`D`$= 23.1 mm2/sec. This falls between the results obtained in previous experiments  and , after correcting for pressure and isotopic mass dependence, and that reported in .

## 5 Summary

The system under consideration consists of spin-1/2 particles in a magnetic field with a time-varying gradient to its strength, and with its initial magnetization, taken to be spatially uniform, having components both transverse to and along the field direction. We now summarize the line of reasoning that led to the expression (57) for the time development of the expectation value of the transverse component of spin, $`<\sigma _+>`$, starting with Eq.(3) for the time derivative of $`\rho `$, the 1-body density matrix. It is the off-diagonal elements of $`\rho `$ in the 2-dimensional spin space, or equivalently the Wigner distribution function $`f_+`$, that determine $`<\sigma _+>`$.

A number of parameters play a role in this development and certain inequalities among them were assumed, in succession. As each of these relations arose we compared the relevant ratios with their values in the experiment reported by Schmidt et al. , expressed in a way that allows quick comparison with any other set of experimental parameters. For a low density system of spin 1/2 particles, in the equation for the 2-body density matrix we neglected 3-body collisions, and this requires that the mean spacing between particles, $`d`$, be large compare to the range of the force, $`a`$. We also assumed that prior to a collision the particles are uncorrelated, so that the 2-body density matrix is the product of the respective 1-body density matrices. Quantum scattering theory is then used to find the effect of a collision.

After expressing the 1-body density matrix elements in terms of Wigner distribution functions, $`f`$ , which form a 2x2 matrix in spin space, a set of coupled nonlinear integro- partial differential equations (10) for the $`f`$’s is obtained. On the right hand side are three contributions to their time rate of change: $`S^K`$ arising from the kinetic energy; $`S^\gamma `$ from the magnetic energy; and $`S^V`$ from the 2-body potential energy terms in the Hamiltonian. The first two terms are simple to write down and are given in Eqs.(11) -(14); the expressions for the collision terms $`S^V`$ are obtained in the Appendix, Eqs.(A-17) and (A-18). Two approximations were made to obtain those collision terms. The first is to assume that the 2-body potential is purely central, and this is satisfied extremely well for 3He atoms. The second assumes that the gradient of the magnetic field strength is sufficiently small; this arises from the spatial dependence of the phase factor that governs the Larmor precession, Eq.(20), where $`F(t)`$ is the time integral of field gradient. We expected and verified that this is the only spatial dependence of the distribution functions. The product $`\gamma F(t)`$ plays the role of a wave number, and leads to the off diagonal element (in spin space) of the density matrix being non-local in momentum space by the amount $`\gamma \mathrm{}F(t)`$; see Eq.(22). This nonlocality shifts the arguments of $`g_+`$ and the $`t`$-matrix elements in the collision terms by amounts of the same order, and even leads to non-energy-conserving matrix elements. If these shifts were sizable they would represent important quantum mechanical effects. The magnitude of the effect of the shift on $`g_+`$ depends on the ratio of $`\gamma \mathrm{}F(t)`$ to the mean thermal momenta, and in Eq.(23) this is shown to be very small for the conditions of the experiment in reference . Using a hard sphere potential as a guide, this same small ratio makes the effect of the shift in the arguments of the $`t`$-matrix elements negligible as well.

To solve the resulting equations (24) - (26) we expanded the distribution functions in a perturbation series to second order in the strength of the field gradient, $`G(t)`$, using the gyromagnetic ratio $`\gamma `$ as the expansion coefficient since it multiplies $`G`$; this makes the equations linear in each order. The only distortions away from Maxwell-Boltzmann distributions that were used are of the form $`𝐮𝐩`$ in first order and $`𝐮𝐮`$ and $`(𝐮𝐩)^2`$ in second order, where $`𝐮`$ is the direction of the gradient. Omitting additional dependence on $`p^2`$ is equivalent to keeping only the leading term in an expansion in Sonine polynomials .

To first order in $`\gamma `$ the equations (A-24) - (A-26) were written still containing the effect of antisymmetization of the wave functions on the collision terms; these are the terms containing $`ϵ`$ and they involve, through the quantities $`I_I`$ and $`I_\pi `$ defined in (A-23), not just the absolute square of the $`t`$-matrix elements but also their real and imaginary parts. For a hard sphere potential in the short wavelength limit, however, the exchange terms are smaller than the direct terms by $`O(\overline{\lambda }/a)^2`$, where $`a`$ is the hard sphere radius, and we expect this to be true for a more general potential with $`a`$ the range of the force. But note from Eq.(38) that for $`{}_{}{}^{3}\mathrm{He}`$ at room temperature the ratio $`a/\overline{\lambda }`$ is only 12.5 with $`a`$ = 2.4 Å. At a significantly lower temperature corrections due to the exchange terms would become important; they were neglected in the remainder of the paper. The direct terms, on the other hand, manifest themselves through the single quantity $`I_U`$, which is the mean of the product of the cube of the momentum with the transport cross section $`\sigma _U`$, averaged over a Maxwell-Boltzmann distribution for a particle having the reduced mass.

Provided the field gradient $`G(t)`$ varies on a time scale much longer than the mean free time between collisions, $`\alpha ^1`$, the solution for the first order distortion of the transverse Wigner distribution function is given in Eqs.(49) and (50). $`\alpha `$ is simply related to the self diffusion constant $`D`$ for 3He, which can be more readily realized in a different experiment with no field gradient but rather a spatial gradient to the $`z`$-component of magnetization; this relation is given in Eq.(46). The solution for the second order contribution to the transverse distribution function follows readily from the first order term, and is given in Eq.(55).

If all the restrictions mentioned above are satisfied, then the expression for the rate of decrease of the transverse component of spin, to second order in the strength of the gradient, is given in Eq.(57). Comparison with the classical treatment described in the Introduction, the exact solution for which is given in Eq.(59), shows agreement to second order in the strength of the gradient. It was also shown that the rate would be much greater if there were no collisions, because the particles would sample a much larger range of magnetic field strengths in a given time and consequently their precession rates would also.

It is not possible to make a good comparison between the magnitude of the 3He self diffusion coefficient estimated from a hard sphere potential using Eqs.(46) and (47), and the experimental result in  because the temperature varied over a large range . The theoretical estimate of $`D`$= 23.1 mm2/sec falls between the values obtained in that experiment and that obtained in previous experiments .

## 6 Acknowledgments

I am indebted to David Schmidt for introducing me to this subject and for many informative discussions, and for valuable suggestions regarding the manuscript. I also want to thank Emil Motolla and Wojciech Zurek for helpful conversations.

## Appendix A Appendix

Here we work out the expressions for the collision terms $`S^V`$. After substituting Eq.(5) for $`\rho ^{(2)}(1,2)`$ in Eq.(18) the elastic scattering $`t`$-matrix is introduced via the relations  $`\mathrm{\Omega }^{(+)}=(I+G^{(+)}t)`$ and $`V\mathrm{\Omega }^{(+)}=t`$, where $`G^{(+)}`$ is the Green’s function. Since the potential conserves the total momentum and the $`t`$-matrix elements depend only on the relative momenta of the two particles, they can be written

$`<𝐩_A,m_A;𝐩_B,m_B|t|𝐩_{}^{}{}_{A}{}^{},m_A^{};𝐩_{}^{}{}_{B}{}^{},m_B^{}>=\delta _{m_Am_A^{}}\delta _{m_Bm_B^{}}`$ (A-1)
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩_A+𝐩_B(𝐩_{}^{}{}_{A}{}^{}+𝐩_{}^{}{}_{B}{}^{}))t(𝐤_{AB},𝐤_{}^{}{}_{AB}{}^{})`$
where $`𝐤_{AB}=(𝐩_A𝐩_B)/2`$. The corresponding equation for the matrix element of $`\mathrm{\Omega }^{(+)}`$ is

$`<𝐩_A,m_A;𝐩_B,m_B|\mathrm{\Omega }^{(+)}|𝐩_{}^{}{}_{A}{}^{},m_A^{};𝐩_{}^{}{}_{B}{}^{},m_B^{}>=\delta _{m_Am_A^{}}\delta _{m_Bm_B^{}}`$ (A-2)
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩_A+𝐩_B(𝐩_{}^{}{}_{A}{}^{}+𝐩_B^{}))\mathrm{\Omega }^{(+)}(𝐤_{AB},𝐤_{AB}^{})`$
where

$$\mathrm{\Omega }^{(+)}(𝐤_{AB},𝐤_{AB}^{})=[(2\pi \mathrm{})^3\delta ^{(3)}(𝐩_A𝐩_A^{})+t(𝐤_{AB},𝐤_{AB}^{})/(E_{AB}^{}+i\tau E_{AB})]$$
(A-3)
with $`E_{AB}`$ the sum of the 1-body energies of particles A and B. (It is understood that the limit $`\tau 0`$ is to be taken.) Since the magnetic energy term in $`H^{(1)}`$ is the same in the primed and unprimed states only the kinetic energy contributes to the denominator in Eq.(A-3). Hence $`E_{AB}`$ can be taken to be $`k_{AB}^2/2\mu `$ where $`\mu `$, the reduced mass, is $`M/2`$.

As described in the text we are assuming that the magnitude of $`\gamma \mathrm{}F(t)`$ is sufficiently small that the shift in the argument of the $`\delta `$-function in Eq.(22) can be neglected. Hence we shall write

$$<𝐩,m|\rho |𝐩^{},m^{}>=(2\pi \mathrm{})^3\delta ^{(3)}(𝐩𝐩^{})\rho _{mm^{}}(𝐩)$$
(A-4)
where

$$\rho _{mm}(𝐩)=(2\pi \mathrm{})^3f_m(𝐩)$$
(A-5)
and

$$\rho _{1/\mathrm{2\; 1}/2}(𝐩)=(2\pi \mathrm{})^3g_+(𝐩)$$
(A-6)
When two sets of intermediate states are introduced into Eq.(18) two of the integrations can be carried out immediately by making use of Eqs.(A-1), (A-2), and (A-4), giving

$`S_{mm^{}}^V(𝐱,𝐩)={\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \frac{1}{(2\pi \mathrm{})^3}}{\displaystyle \underset{m_2}{}}{\displaystyle \frac{d^3p_2}{(2\pi \mathrm{})^3}\frac{d^3p^{}}{(2\pi \mathrm{})^3}\frac{d^3p_2^{}}{(2\pi \mathrm{})^3}}`$ (A-7)
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩+𝐩_\mathrm{𝟐}(𝐩^{}+𝐩_\mathrm{𝟐}^{}))`$
$`\{[t(𝐤,𝐤^{})\mathrm{\Omega }^{(+)}(𝐤,𝐤^{})\mathrm{\Omega }^{(+)}(𝐤,𝐤^{})t^{}(𝐤,𝐤^{})]\rho _{mm^{}}(𝐩^{})\rho _{m_2m_2}(𝐩_2^{})`$
$`+ϵ[t(𝐤,𝐤^{})\mathrm{\Omega }^{(+)}(𝐤,𝐤^{})\mathrm{\Omega }^{(+)}(𝐤,𝐤^{})t^{}(𝐤,𝐤^{})]\rho _{m_2m^{}}(𝐩^{})\rho _{mm_2}(𝐩_2^{})\}`$
After substituting for $`\mathrm{\Omega }^{(+)}`$ from Eq.(A-3), for the terms linear in the $`t`$-matrix elements the integrations over the primed variables can be carried out, and for the terms bilinear in $`t`$ the energy denominators can be combined, giving

$`S_{mm^{}}^V(𝐱,𝐩)={\displaystyle \frac{i}{\mathrm{}}}\{{\displaystyle \frac{1}{(2\pi \mathrm{})^3}}{\displaystyle \underset{m_2}{}}{\displaystyle }{\displaystyle \frac{d^3p_2}{(2\pi \mathrm{})^3}}([t(𝐤,𝐤)t^{}(𝐤,𝐤)]\rho _{mm^{}}(𝐩)\rho _{m_2m_2}(𝐩_\mathrm{𝟐})`$ (A-8)
$`+ϵ[t(𝐤,𝐤)\rho _{m_2m^{}}(𝐩)\rho _{mm_2}(𝐩_\mathrm{𝟐})t^{}(𝐤,𝐤)\rho _{m_2m^{}}(𝐩_\mathrm{𝟐})\rho _{mm_2}(𝐩)])`$
$`+{\displaystyle \frac{1}{(2\pi \mathrm{})^3}}{\displaystyle \underset{m_2}{}}{\displaystyle \frac{d^3p_2}{(2\pi \mathrm{})^3}\frac{d^3p^{}}{(2\pi \mathrm{})^3}\frac{d^3p_2^{}}{(2\pi \mathrm{})^3}}`$
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩+𝐩_\mathrm{𝟐}(𝐩^{}+𝐩_\mathrm{𝟐}^{}))2\pi i\delta (E^{}E)`$
$`[|t(𝐤,𝐤^{})|^2\rho _{mm^{}}(𝐩^{})\rho _{m_2m_2}(𝐩_\mathrm{𝟐}^{})+ϵt(𝐤,𝐤^{})t^{}(𝐤,𝐤^{})\rho _{m_2m^{}}(𝐩^{})\rho _{mm_2}(𝐩_\mathrm{𝟐}^{})]\}`$
This establishes the fact that only energy conserving matrix elements contribute to the collision terms, so $`t(𝐤,𝐤^{})`$ can be written as $`t(k,\theta )`$.

Note that any integral of the form

$`I(𝐩,𝐩_2)={\displaystyle \frac{d^3p^{}}{(2\pi \mathrm{})^3}\frac{d^3p_2^{}}{(2\pi \mathrm{})^3}}`$ (A-9)
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩+𝐩_\mathrm{𝟐}(𝐩^{}+𝐩_\mathrm{𝟐}^{}))2\pi \delta (E^{}E)F(𝐩,𝐩_\mathrm{𝟐},𝐩^{},𝐩_\mathrm{𝟐}^{})`$
can be rewritten by changing the variables of integration to ($`𝐏^{},𝐤^{}`$), where $`𝐏^{}=𝐩^{}+𝐩_\mathrm{𝟐}^{}`$ and $`𝐤^{}=(𝐩^{}𝐩_\mathrm{𝟐}^{})/\mathrm{𝟐}`$. After carrying out the integrations over $`𝐏^{}`$ and the magnitude of $`𝐤^{}`$, using $`E^{}E=(k^2k^2)/2\mu `$, $`I`$ becomes

$$I(𝐩,𝐩_2)=2\pi \frac{\mu kd\mathrm{\Omega }^{}}{(2\pi \mathrm{})^3}F(𝐩,𝐩_2,\mathrm{\Omega }^{})$$
(A-10)
where $`\mathrm{\Omega }^{}`$ stands for two angles specifying the direction of $`𝐤^{}`$. The combination $`\mu kd\mathrm{\Omega }^{}/(2\pi \mathrm{})^3`$ is the density of final states.

Making use of the relation between the differential scattering cross section and the $`t`$-matrix elements gives

$$v\frac{d\sigma _U}{d\mathrm{\Omega }^{}}=\frac{2\pi }{\mathrm{}}|t(𝐤,𝐤^{})|^2\frac{\mu k}{(2\pi \mathrm{})^3}$$
(A-11)
$$v\frac{d\sigma _A}{d\mathrm{\Omega }^{}}=\frac{2\pi }{\mathrm{}}\left|[t(𝐤,𝐤^{})+ϵt(𝐤,𝐤^{})]/\sqrt{2}\right|^2\frac{\mu k}{(2\pi \mathrm{})^3}$$
(A-12)
and

$$v\frac{d\sigma _I}{d\mathrm{\Omega }^{}}=\frac{2\pi }{\mathrm{}}t(𝐤,𝐤^{})t^{}(𝐤,𝐤^{})\frac{\mu k}{(2\pi \mathrm{})^3}$$
(A-13)
where $`\sigma _U`$ and $`\sigma _A`$ are, respectively, the unsymmetrized and antisymmetrized cross sections. $`\sigma _I`$ is not actually a cross section, just a convenient notation for the interference term; $`v=k/\mu `$ is the relative velocity of the colliding particles.

Further progress is made by making use of rotational invariance of the $`t`$-matrix, and unitarity of the $`S`$-matrix, which takes the form

$`2i\mathrm{Im}t(𝐤,𝐤)=i{\displaystyle \frac{d^3p^{}}{(2\pi \mathrm{})^3}\frac{d^3p_2^{}}{(2\pi \mathrm{})^3}|t(𝐤,𝐤^{})|^2}`$ (A-14)
$`(2\pi \mathrm{})^3\delta ^{(3)}(𝐩+𝐩_\mathrm{𝟐}(𝐩^{}+𝐩_\mathrm{𝟐}^{}))2\pi \delta (E^{}E),`$
and from the expressions above this becomes

$$2i\mathrm{Im}t(𝐤,𝐤)=i\mathrm{}𝑑\mathrm{\Omega }^{}v\frac{d\sigma _U}{d\mathrm{\Omega }^{}}.$$
(A-15)
The corresponding equation for $`t(𝐤,𝐤)`$ is

$$2i\mathrm{Im}t(𝐤,𝐤)=i\mathrm{}𝑑\mathrm{\Omega }^{}v\frac{d\sigma _I}{d\mathrm{\Omega }^{}}.$$
(A-16)
In terms of these quantities the separate cases of $`S^V`$ from Eq.(A-8) can be obtained by carrying out the summation over $`m_2`$ and rearranging some terms.

$`S_{1/2}^V(𝐩)={\displaystyle }d^3p_2d\mathrm{\Omega }^{}\{[f_{1/2}(𝐩^{})f_{1/2}(𝐩_2^{})f_{1/2}(𝐩)f_{1/2}(𝐩_2)]v{\displaystyle \frac{d\sigma _U}{d\mathrm{\Omega }^{}}}`$ (A-17)
$`+[f_{1/2}(𝐩^{})f_{1/2}(𝐩_2^{})f_{1/2}(𝐩)f_{1/2}(𝐩_2)]v{\displaystyle \frac{d\sigma _A}{d\mathrm{\Omega }^{}}}+ϵg_+(𝐩^{})g_+^{}(𝐩_2^{})v{\displaystyle \frac{d\sigma _I}{d\mathrm{\Omega }^{}}}\}`$
$`ϵ{\displaystyle \frac{i}{\mathrm{}}}{\displaystyle }d^3p_2[t(𝐤,𝐤)g_+(𝐩)g_+^{}(𝐩_\mathrm{𝟐})c.c.]`$
The corresponding equation for $`S_{1/2}^V`$ is obtained from (A-17) by making the interchanges $`f_{1/2}f_{1/2}`$ and $`g_+g_+^{}`$. Finally, the expression for $`\overline{S}_+^V`$ becomes

$`\overline{S}_+^V(𝐩)=\mathrm{exp}[i\gamma (B_0t+𝐮𝐱F(t))]S_+^V(𝐱,𝐩)={\displaystyle d^3p_2𝑑\mathrm{\Omega }^{}}`$ (A-18)
$`\{(g_+(𝐩^{})[f_{1/2}(𝐩_\mathrm{𝟐}^{})+f_{1/2}(𝐩_\mathrm{𝟐}^{})]g_+(𝐩)[f_{1/2}(𝐩_\mathrm{𝟐})+f_{1/2}(𝐩_\mathrm{𝟐})])v{\displaystyle \frac{d\sigma _U}{d\mathrm{\Omega }^{}}}`$
$`+ϵ[g_+(𝐩^{})f_{1/2}(𝐩_\mathrm{𝟐}^{})+f_{1/2}(𝐩^{})g_+(𝐩_\mathrm{𝟐}^{})]v{\displaystyle \frac{d\sigma _I}{d\mathrm{\Omega }^{}}}\}`$
$`ϵ{\displaystyle \frac{i}{\mathrm{}}}{\displaystyle }d^3p_2(t(𝐤,𝐤)[g_+(𝐩)f_{1/2}(𝐩_2)+f_{1/2}(𝐩)g_+(𝐩_2)]`$
$`t^{}(𝐤,𝐤)[g_+(𝐩_2)f_{1/2}(𝐩)+f_{1/2}(𝐩_2)g_+(𝐩)])`$
We now extract the collision terms of first order in $`\gamma `$, to be used in Eqs.(34) - (36), by inserting Eqs.(32) and (33) in the above. It is helpful to express each individual momentum variable in terms of the total and relative momenta $`𝐏,𝐤`$ and $`𝐤^{}`$. When this is done the terms in $`S_{\pm 1/2}^{V(1)}`$ involving $`d\sigma _A/d\mathrm{\Omega }^{}`$ vanish by conservation of momentum; and making use of unitarity, Eq.(A-16), the terms involving the real part of $`h_+^{(1)}`$ cancel, leaving

$`S_{1/2}^{V(1)}(𝐩)=𝐮{\displaystyle }d^3p_2d\mathrm{\Omega }^{}B^{(M)}(p)B^{(M)}(p_2)\{n_{1/2}n_{1/2}(𝐩^{}𝐩)(h_{1/2}^{(1)}h_{1/2}^{(1)})v{\displaystyle \frac{d\sigma _U}{d\mathrm{\Omega }^{}}}`$ (A-19)
$`+ϵ|n_+|^22i𝐤^{}\mathrm{Im}h_+^{(1)}v{\displaystyle \frac{d\sigma _I}{d\mathrm{\Omega }^{}}}\}`$
$`+ϵ|n_+|^2{\displaystyle \frac{4}{\mathrm{}}}\mathrm{Im}h_+^{(1)}𝐮{\displaystyle d^3p_2B^{(M)}(p)B^{(M)}(p_2)𝐤\mathrm{Re}t(𝐤,𝐤)}.`$
Furthermore, $`S_{1/2}^{V(1)}(𝐩)`$ = $`S_{1/2}^{V(1)}(𝐩)`$. Following similar reasoning leads to

$`\overline{S}_+^{V(1)}(𝐩)=𝐮n_+{\displaystyle d^3p_2𝑑\mathrm{\Omega }^{}B^{(M)}(p)B^{(M)}(p_2)}`$ (A-20)
$`\{(𝐩^{}𝐩)[h_+^{(1)}(n_{1/2}+n_{1/2})(n_{1/2}h_{1/2}^{(1)}+n_{1/2}h_{1/2}^{(1)})]v{\displaystyle \frac{d\sigma _U}{d\mathrm{\Omega }^{}}}`$
$`+ϵ𝐤^{}[(n_{1/2}h_{1/2}^{(1)}n_{1/2}h_{1/2}^{(1)})(n_{1/2}n_{1/2})h_+^{(1)}]v{\displaystyle \frac{d\sigma _I}{d\mathrm{\Omega }^{}}}\}`$
$`iϵ{\displaystyle \frac{2}{\mathrm{}}}n_+[(n_{1/2}h_{1/2}^{(1)}n_{1/2}h_{1/2}^{(1)})(n_{1/2}n_{1/2})h_+^{(1)}]`$
$`𝐮{\displaystyle d^3p_2B^{(M)}(p)B^{(M)}(p_2)𝐤\mathrm{Re}t(𝐤,𝐤)}`$
When these expressions for the first order collision terms are inserted into Eqs.(34) - (36) we see that every term is of the form $`𝐮`$ some vector, where $`𝐮`$ is an arbitrary direction in space; therefore those equations must be valid as vector equations with the $`𝐮`$ omitted. A necessary consequence can be obtained by forming $`𝐩`$ those equations and integrating over all $`𝐩`$. At the same time we make the following observations about the collision terms. Change the variables of integration from $`(𝐩,𝐩_\mathrm{𝟐})`$ to $`(𝐏,𝐤)`$, where $`𝐏=𝐩+𝐩_\mathrm{𝟐}`$ and $`𝐤=(𝐩𝐩_\mathrm{𝟐})/2`$, and make use of

$$B^{(M)}(p)B^{(M)}(p_2)=B^{(2M)}(P)B^{(\mu )}(k)$$
(A-21)
which follows from conservation of energy. For each fixed value of $`𝐤`$ choose its direction as the polar axis for the integration over the angles of $`𝐤^{}`$, so that $`d\mathrm{\Omega }^{}`$ becomes $`d\mathrm{\Omega }=2\pi \mathrm{sin}\theta d\theta `$ where $`\theta `$ is the angle between $`𝐤`$ and $`𝐤^{}`$. We also define the following quantities.

$`\sigma _U(k)`$ $`=`$ $`{\displaystyle _0^\pi }𝑑\mathrm{\Omega }(1\mathrm{cos}\theta ){\displaystyle \frac{d\sigma _U}{d\mathrm{\Omega }}}`$
$`\sigma _I(k)`$ $`=`$ $`{\displaystyle _0^\pi }𝑑\mathrm{\Omega }\mathrm{cos}\theta {\displaystyle \frac{d\sigma _I}{d\mathrm{\Omega }}}`$ (A-22)
and

$`I_U`$ $`=`$ $`{\displaystyle d^3kk^3B^{(\mu )}(k)\sigma _U(k)}=<k^3\sigma _U(k)>^{(\mu )}`$
$`I_I`$ $`=`$ $`{\displaystyle d^3kk^3B^{(\mu )}(k)\sigma _I(k)}=<k^3\sigma _I(k)>^{(\mu )}`$
$`I_\pi `$ $`=`$ $`{\displaystyle d^3kk^2B^{(\mu )}(k)\mathrm{Re}t(k,\theta =\pi )}`$ (A-23)
From the definition of $`d\sigma _I/d\mathrm{\Omega }`$, Eq.(A-13), it is seen that $`\sigma _I(k)`$ and $`I_I`$ are pure imaginary, whereas $`I_U`$ and $`I_\pi `$ are real. $`\sigma _U(k)`$ is called the transport cross section  and $`I_U`$ is the mean value of $`k^3\sigma _U(k)`$ averaged over a Maxwell-Boltzmann distribution for a particle of mass $`\mu `$.

Eqs.(34) - (36) now become

$`n_{1/2}{\displaystyle \frac{dh_{1/2}^{(1)}}{dt}}=n_{1/2}\mathrm{}\beta ^{(M)}G(t)`$ (A-24)
$`+{\displaystyle \frac{4\beta ^{(M)}}{3M}}\left\{n_{1/2}n_{1/2}[h_{1/2}^{(1)}h_{1/2}^{(1)}]I_U2ϵ|n_+|^2\mathrm{Im}h_+^{(1)}\mathrm{Im}I_I\right\}`$
$`+ϵ{\displaystyle \frac{8\beta ^{(M)}}{3\mathrm{}}}|n_+|^2\mathrm{Im}h_+^{(1)}I_\pi ,`$
$`n_{1/2}{\displaystyle \frac{dh_{1/2}^{(1)}}{dt}}=n_{1/2}\mathrm{}\beta ^{(M)}G(t)`$ (A-25)
$`{\displaystyle \frac{4\beta ^{(M)}}{3M}}\left\{n_{1/2}n_{1/2}[h_{1/2}^{(1)}h_{1/2}^{(1)}]I_U2ϵ|n_+|^2\mathrm{Im}h_+^{(1)}\mathrm{Im}I_I\right\}`$
$`ϵ{\displaystyle \frac{8\beta ^{(M)}}{3\mathrm{}}}|n_+|^2\mathrm{Im}h_+^{(1)}I_\pi `$
and

$`{\displaystyle \frac{dh_+^{(1)}}{dt}}=i{\displaystyle \frac{1}{M}}F(t)+{\displaystyle \frac{4\beta ^{(M)}}{3M}}\{[n_{1/2}h_{1/2}^{(1)}+n_{1/2}h_{1/2}^{(1)}(n_{1/2}+n_{1/2})h_+^{(1)}]I_U`$ (A-26)
$`+ϵ[(n_{1/2}h_{1/2}^{(1)}n_{1/2}h_{1/2}^{(1)})(n_{1/2}n_{1/2})h_+^{(1)}]I_I\}`$
$`iϵ{\displaystyle \frac{4\beta ^{(M)}}{3\mathrm{}}}[(n_{1/2}h_{1/2}^{(1)}n_{1/2}h_{1/2}^{(1)})(n_{1/2}n_{1/2})h_+^{(1)}]I_\pi .`$
These are four coupled ordinary linear differential equations for the four unknown functions $`h_{\pm 1/2}^{(1)}`$ and the real and imaginary parts of $`h_+^{(1)}`$.


format: text

text:
# Gravitational lensing in the strong field limit

## I Introduction

Gravitational lensing is one of the first applications of General Relativity ever studied Einstein . Firstly it was recognized in the light deflection by sun, secondly in lensing of quasars by foreground galaxies, then in the formation of giant arcs in galaxy clusters and finally in galactic microlensing. Now it is an ordinary phenomenon in the panorama of astronomical observations (see SEF  for a complete treatment and references therein).

The full theory of gravitational lensing has been developed following the scheme of the weak field approximation and, in this formulation, it has been successfully employed to explain all the physical observations.

In the last years, however, the scientific community is starting to look at this phenomenon from the opposite point of view, opening a strong field perspective. Viergutz Vie  made a semi-analytical investigation about geodesics in Kerr geometry; in Ref. Bar  the appearance of a black hole in front of a uniform background was studied; Falcke, Melia and Agol FMA  considered the emission of the accretion flow as source. Virbhadra & Ellis VirEll  showed that a source behind a Schwarzschild black hole would produce one set of infinite relativistic images on each side of the black hole. These images are produced when a light ray with small impact parameter winds one or several times around the black hole before emerging. Later on, by an alternative formulation of the problem, Frittelli, Kling & Newman FKN  attained an exact lens equation, giving integral expressions for its solutions, and compared their results to those by Virbhadra & Ellis. The same problem has been investigated by Bozza et al. in Ref. BCIS , where a strong field limit was first defined in Schwarzschild black hole lensing and used to find the position and the characteristics of all the relativistic images analytically. Eiroa, Romero and Torres ERT  applied the same technique to a Reissner - Nordstrom black hole. Recently, in another work VirEll2 , Virbhadra & Ellis distinguished the main features of gravitational lensing by normal black holes and by naked singularities, analyzing the Janis, Newman, Winicour metric. They remarked the importance of these studies in providing a test for the cosmic censorship hypothesis.

The reason for such an interest in gravitational lensing in strong fields is that by the properties of the relativistic images it may be possible to investigate the regions immediately outside of the event horizon. High resolution imaging of black holes by VLBI VLBI  could be able to detect the relativistic images and retrieve information about strong fields stored within these new observables. Moreover, since alternative theories of gravitation must agree with GR in the weak field limit, in order to show deviations from GR it is necessary to probe strong fields in some way. Indeed, deviation of light rays in strong fields is one of the most promising grounds where a theory of gravitation can be tested in its full form.

Of course, the study of null geodesics in strong fields is not easy and up to now it has always been carried out using numerical techniques. An analytical treatment would enlighten the dependence of the observables on the parameters of the system, allow easy checks about the detectability of the images and open the way to comparisons between the results in different metrics. In Ref. BCIS , a new way to expand the deflection angle in the Schwarzschild metric was suggested. The deflection angle near its divergence was approximated by its leading order and its first regular term and then plugged into the lens equation. In this way, very simple and reliable analytical formulae were derived for the relativistic images and their main features.

So, as the weak field limit takes the first order deviation from Minkowski, the strong field limit starts from complete capture of the photon and takes the leading order in the divergence of the deflection angle.

The strong field limit of Ref. BCIS  was only developed in Schwarzschild spacetime. In this paper, we provide a general method to extend the strong field limit to a generic static spherically symmetric spacetime. Our method is universal and can be applied to any spacetime in any theory of gravitation, provided that photons satisfy the standard geodesics equation. The parameters of the strong field limit expansion are directly connected with the observables, providing an effective tool to discriminate among different metrics. In Sect. 2, we state the problem and carry out the strong field limit of the deflection angle. In Sect. 3, we apply the method to some simple metrics: Schwarzschild, Reissner-Nordstrom and Janis-Newman-Winicour black hole, discussing their differences with reference to the gravitational lensing phenomenology. In Sect. 4, we establish a connection between the strong field limit coefficients and the relativistic images, analyzing the case of the black hole at the center of our galaxy as a concrete example where our results can be tested. Finally, Sect. 5 contains the summary.

## II Strong field expansion of the deflection angle

A generic spherically symmetric spacetime has line-element Wei

$$ds^2=A(x)dt^2B(x)dx^2C(x)\left(d\theta ^2+\mathrm{sin}^2\theta \varphi ^2\right).$$
(1)
where

$$\begin{array}{c}A(x)^x\mathrm{}1\frac{2M}{x}\\ \\ B(x)^x\mathrm{}1+\frac{2M}{x}\\ \\ C(x)^x\mathrm{}x^2,\end{array}$$
(2)
in order to correctly match the weak gravitational field far from the lensing object.

We require that the equation

$$\frac{C^{}(x)}{C(x)}=\frac{A^{}(x)}{A(x)}.$$
(3)
admits at least one positive solution. We shall call the largest root of Eq. (3) the radius of the photon sphere $`x_m`$ (for an alternative definition of the photon sphere, see CVE ). $`A`$, $`B`$, $`C`$, $`A^{}`$ and $`C^{}`$ must be positive for $`x>x_m`$

For metrics expressed in standard coordinates ($`C(x)=x^2`$) a sufficient condition for the existence of $`x_m`$ is the presence of a static limit (a radius $`x_s`$ such that $`A(x_s)=0`$). Our strong field expansion takes the photon sphere as the starting point. In our study, therefore, we shall not consider naked singularities without a photon sphere. For a numerical study of this situation, see Ref. VirEll2 .

A photon incoming from infinity with some impact parameter $`u`$, will be deviated while approaching the black hole. It will reach a minimum distance $`x_0`$ and then emerge in another direction. Of course, the approach phase is identical to the departure phase, with the time reversed.

By the conservation of the angular momentum, the closest approach distance is related to the impact parameter $`u`$ by

$$u=\sqrt{\frac{C_0}{A_0}}.$$
(4)
where the subscript $`0`$ indicates that the function is evaluated at $`x_0`$.

From the geodesics equation it is easy to extract the quantity

$$\frac{d\varphi }{dx}=\frac{\sqrt{B}}{\sqrt{C}\sqrt{\frac{C}{C_0}\frac{A_0}{A}1}}$$
(5)
which gives the angular shift of the photon as a function of the distance from the center (see Wei  for the complete derivation).

The deflection angle can then be calculated as a function of the closest approach:

$`\alpha (x_0)=I(x_0)\pi `$ (6)
$`I(x_0)={\displaystyle \underset{x_0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{2\sqrt{B}dx}{\sqrt{C}\sqrt{\frac{C}{C_0}\frac{A_0}{A}1}}}.`$ (7)
It is easy to check that for a vanishing gravitational field ($`A=B=1`$, $`C=x^2`$) $`\alpha (x_0)`$ identically vanishes. In the weak field limit, the integrand is expanded to the first order in the gravitational potential. This limit is no longer valid when the closest approach distance significantly differs from the impact parameter (which, by Eq. (4), means that $`A(x_0)`$ significantly differs from 1 and/or $`C(x_0)`$ from $`x_0^2`$, i.e. the photon passes in a strong gravitational field).

When we decrease the impact parameter (and consequently $`x_0`$), the deflection angle increases. At some point, the deflection angle will exceed $`2\pi `$, resulting in a complete loop around the black hole. Decreasing $`u`$ further, the photon will wind several times before emerging. Finally, for $`x_0=x_m`$ (see Sect. II.1), corresponding to an impact parameter $`u=u_m`$, the deflection angle diverges and the photon is captured (see Fig. 1).

In this section, we will show that this divergence is logarithmic for all spherically symmetric metrics. Our purpose is to get an analytical expansion of the deflection angle close to the divergence in the form

$$\alpha (x_0)=a\mathrm{log}\left(\frac{x_0}{x_m}1\right)+b+O\left(x_0x_m\right),$$
(8)
where all the coefficients depend on the metric functions evaluated at $`x_m`$.

As the angular separation of the image from the lens is $`\theta =\frac{u}{D_{OL}}`$, where $`D_{OL}`$ is the distance between the lens and the observer, we need to express the deflection angle in terms of this variable. So, we shall finally transform Eq.(8) to

$$\alpha (\theta )=\overline{a}\mathrm{log}\left(\frac{\theta D_{OL}}{u_m}1\right)+\overline{b}+O\left(uu_m\right)$$
(9)
which we define as the strong field limit of the deflection angle. The rest of this section is devoted to the calculation of the two coefficients $`a`$ and $`b`$ in a first step and then of $`\overline{a}`$ and $`\overline{b}`$.

### II.1 Divergent term of the deflection angle

We define two new variables

$`y=A(x)`$ (10)
$`z={\displaystyle \frac{yy_0}{1y_0}}`$ (11)
where $`y_0=A_0`$. The integral (7) in the deflection angle becomes

$`I(x_0)={\displaystyle \underset{0}{\overset{1}{}}}R(z,x_0)f(z,x_0)𝑑z`$ (12)
$`R(z,x_0)={\displaystyle \frac{2\sqrt{By}}{CA^{}}}\left(1y_0\right)\sqrt{C_0}`$ (13)
$`f(z,x_0)={\displaystyle \frac{1}{\sqrt{y_0\left[\left(1y_0\right)z+y_0\right]\frac{C_0}{C}}}}`$ (14)
where all functions without the subscript $`0`$ are evaluated at $`x=A^1\left[\left(1y_0\right)z+y_0\right]`$.

The function $`R(z,x_0)`$ is regular for all values of $`z`$ and $`x_0`$, while $`f(z,x_0)`$ diverges for $`z0`$. To find out the order of divergence of the integrand, we expand the argument of the square root in $`f(z,x_0)`$ to the second order in $`z`$:

$`f`$ $`(z,x_0)f_0(z,x_0)={\displaystyle \frac{1}{\sqrt{\alpha z+\beta z^2}}}`$ (15)
$`\alpha `$ $`={\displaystyle \frac{1y_0}{C_0A_0^{}}}\left(C_0^{}y_0C_0A_0^{}\right)`$ (17)
$`\beta `$ $`={\displaystyle \frac{\left(1y_0\right)^2}{2C_0^2A_{0}^{}{}_{}{}^{3}}}[2C_0C_0^{}A_{0}^{}{}_{}{}^{2}+`$
$`+(C_0C_0^{\prime \prime }2C_{0}^{}{}_{}{}^{2})y_0A_0^{}C_0C_0^{}y_0A_0^{\prime \prime }].`$
When $`\alpha `$ is non zero, the leading order of the divergence in $`f_0`$ is $`z^{1/2}`$, which can be integrated to give a finite result. When $`\alpha `$ vanishes, the divergence is $`z^1`$ which makes the integral diverge. Examining the form of $`\alpha `$, we see that it vanishes at $`x_0=x_m`$, with $`x_m`$ defined by Eq. (3). Each photon having $`x_0<x_m`$ is captured by the central object and cannot emerge back.

To solve the integral (12), we split it into two pieces

$$I(x_0)=I_D(x_0)+I_R(x_0),$$
(18)
where

$$I_D(x_0)=\underset{0}{\overset{1}{}}R(0,x_m)f_0(z,x_0)𝑑z$$
(19)
contains the divergence and

$$I_R(x_0)=\underset{0}{\overset{1}{}}g(z,x_0)𝑑z$$
(20)
$$g(z,x_0)=R(z,x_0)f(z,x_0)R(0,x_m)f_0(z,x_0).$$
(21)
is the original integral with the divergence subtracted. We shall solve both integrals separately and then sum up their results to rebuild the deflection angle. Here we deal with $`I_D`$ and its divergence, while in the next subsection we shall verify that $`I_R`$ is indeed regular.

The integral $`I_D(x_0)`$ can be solved exactly, giving

$$I_D(x_0)=R(0,x_m)\frac{2}{\sqrt{\beta }}\mathrm{log}\frac{\sqrt{\beta }+\sqrt{\alpha +\beta }}{\sqrt{\alpha }}.$$
(22)
Since we are interested in the terms up to $`O(x_0x_m)`$, we expand $`\alpha `$ as

$$\alpha =\frac{2\beta _mA_m^{}}{1y_m}\left(x_0x_m\right)+O\left(x_0x_m\right)^2,$$
(23)
where

$$\beta _m=\beta |_{x_0=x_m}=\frac{C_m\left(1y_m\right)^2\left(C_m^{\prime \prime }y_mC_mA^{\prime \prime }(x_m)\right)}{2y_m^2C_{m}^{}{}_{}{}^{2}}$$
(24)
and substitute in $`I_D(x_0)`$. Rearranging all terms, we find

$$I_D(x_0)=a\mathrm{log}\left(\frac{x_0}{x_m}1\right)+b_D+O(x_0x_m),$$
(25)
$`a={\displaystyle \frac{R(0,x_m)}{\sqrt{\beta _m}}}`$ (26)
$`b_D={\displaystyle \frac{R(0,x_m)}{\sqrt{\beta _m}}}\mathrm{log}{\displaystyle \frac{2(1y_m)}{A_m^{}x_m}}.`$ (27)
Eq. (25) yields the leading order in the divergence of the deflection angle, which is logarithmic, as anticipated before. The coefficient $`a`$ of Eq. (8) is then given by Eq. (26).

### II.2 Regular term of the deflection angle

In order to find the correct coefficient $`b`$ in Eq. (8), we have to add to the term $`b_D`$ coming from Eq. (27), an analogous term coming from the regular part of the original integral, defined by Eq. (20).

We can expand $`I_R(x_0)`$ in powers of $`(x_0x_m)`$

$$I_R(x_0)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{n!}(x_0x_m)^n\underset{0}{\overset{1}{}}\frac{^ng}{x_0^n}|_{x_0=x_m}dz$$
(28)
and evaluate the single coefficients.

If we had not subtracted the singular part from $`R(z,x_0)f(z,x_0)`$, we would have an infinite coefficient for $`n=0`$, while all other coefficients would be finite. However, the function $`g(z,x_0)`$ is regular in $`z=0,x_0=x_m`$ as can be easily checked by a power expansion, recalling that $`\alpha _m=0`$.

Since we are interested to terms up to $`O(x_0x_m)`$, we will just retain the $`n=0`$ term

$$I_R(x_0)=\underset{0}{\overset{1}{}}g(z,x_m)𝑑z+O(x_0x_m)$$
(29)
and then

$$b_R=I_R(x_m)$$
(30)
is the term we need to add to $`b_D`$ in order to get the regular coefficient. Recalling also the term $`\pi `$ in the deflection angle, we have

$$b=\pi +b_D+b_R.$$
(31)
The coefficient $`b_R`$ can be easily evaluated numerically for all metrics, since the integrand has no divergence. However, in many cases it is also possible to build a completely analytical formula for $`b_R`$ as well. In fact, in Schwarzschild metric, the integral (29) is solved exactly (see Sect. III.1). Then, in most metrics, we can expand Eq. (29) in powers of their parameters, starting from the Schwarzschild limit, and evaluate each coefficient separately. This is what we shall do for Reissner-Nordstrom and Janis-Newman-Winicour black holes (see Sect. III).

### II.3 From $`\alpha (x_0)`$ to $`\alpha (\theta )`$

From Eq. (4), we see that the minimum impact parameter is

$$u_m=\sqrt{\frac{C_m}{y_m}}.$$
(32)
Expanding Eq. (4), we find

$`uu_m=c\left(x_0x_m\right)^2`$ (33)
$`c={\displaystyle \frac{C_m^{\prime \prime }y_mC_mA_m^{\prime \prime }}{4\sqrt{y_m^3C_m}}}=\beta _m\sqrt{{\displaystyle \frac{y_m}{C_m^3}}}{\displaystyle \frac{C_{m}^{}{}_{}{}^{2}}{2\left(1y_m\right)^2}}.`$ (34)
Using this relation, we can write the deflection angle as a function of $`\theta `$:

$`\alpha (\theta )=\overline{a}\mathrm{log}\left({\displaystyle \frac{\theta D_{OL}}{u_m}}1\right)+\overline{b}`$ (35)
$`\overline{a}={\displaystyle \frac{a}{2}}={\displaystyle \frac{R(0,x_m)}{2\sqrt{\beta _m}}}`$ (36)
$`\overline{b}=b+{\displaystyle \frac{a}{2}}\mathrm{log}{\displaystyle \frac{cx_m^2}{u_m}}=\pi +b_R+\overline{a}\mathrm{log}{\displaystyle \frac{2\beta _m}{y_m}}.`$ (37)
This concludes our general discussion of the form of the deflection angle in the strong field limit. Even if the proof is somewhat tricky, the application to concrete cases is very straightforward, as we shall see in Sect. III. In fact, once we write a metric, it is sufficient to:

* solve Eq. (3) to find $`x_m`$;
* write $`\beta _m`$ from Eq. (24) and $`R(0,x_m)`$ from Eq. (13);
* compute $`b_R`$ from Eq. (30) numerically or by a proper expansion in the parameters of the metric;
* compute the coefficients $`u_m`$, $`\overline{a}`$ and $`\overline{b}`$ from Eqs. (32), (36) and (37), respectively.

The crucial step is the calculation of $`b_R`$, since it is the only integral involved in the whole procedure.

## III Applications

In this section, we apply the general formulation of the strong field limit to three simple examples. First, we shall revisit Schwarzschild spacetime, whose strong field limit gravitational lensing has already been studied in Ref. BCIS . Then we will apply our method to Reissner-Nordstrom spacetime, which was explored numerically in Ref. ERT . Finally, as an example of extended theory of gravitation, we shall consider a Janis, Newman, Winicour black hole, whose strong field limit expansion has not been investigated so far. In all three cases, we shall derive analytical formulae for the strong field limit coefficients, in order to analyze the functional dependences of the deflection angle on the parameters of the metrics.

### III.1 Schwarzschild lensing

This is the simplest spherically symmetric metric describing the outer solution for a black hole. It only depends on the mass of the central object (by Birkhoff’s theorem). It is convenient to define the Schwarzschild radius $`x_s=2M`$ as the measure of distances; then, in standard coordinates, the functions in the metric become

$`A(x)=1{\displaystyle \frac{1}{x}}`$ (38)
$`B(x)=\left(1{\displaystyle \frac{1}{x}}\right)^1`$ (39)
$`C(x)=x^2,`$ (40)
which obviously satisfy all hypotheses required in Sect. II, with static limit $`x_s=1`$.

The two functions $`R(z,x_0)`$ and $`f(z,x_0)`$ read

$`R(z,x_0)=2`$ (41)
$`f(z,x_0)={\displaystyle \frac{1}{\sqrt{\left(2\frac{3}{x_0}\right)z+\left(\frac{3}{x_0}1\right)z^2\frac{z^3}{x_0}}}}.`$ (42)
From Eqs. (17), (17), or directly from the expansion of the denominator of $`f`$, we read off the coefficients $`\alpha `$ and $`\beta `$

$`\alpha =2{\displaystyle \frac{3}{x_0}}`$ (43)
$`\beta ={\displaystyle \frac{3}{x_0}}1.`$ (44)
The equation $`\alpha =0`$ defines the radius of the photon sphere

$$x_m=\frac{3}{2}.$$
(45)
Consequently, $`\beta _m=1`$.

In this simple case, it is possible to solve the integral in Eq. (30) exactly and write the regular term in the deflection angle

$$b_R=2\mathrm{log}\left[6\left(2\sqrt{3}\right)\right]=0.9496.$$
(46)
From Eqs. (36) and (37), we derive the coefficients $`\overline{a}`$, $`\overline{b}`$ and $`u_m`$ of the deflection angle

$`\overline{a}=1`$ (47)
$`\overline{b}=\pi +b_R+\mathrm{log}6=0.4002`$ (48)
$`u_m={\displaystyle \frac{3\sqrt{3}}{2}}.`$ (49)
Then the Schwarzschild deflection angle, in the strong field limit, is

$$\alpha (\theta )=\mathrm{log}\left(\frac{2\theta D_{OL}}{3\sqrt{3}}1\right)+\mathrm{log}\left[216\left(74\sqrt{3}\right)\right]\pi ,$$
(50)
as obtained in Ref. BCIS .

In this simple case, we can compare the exact deflection angle $`\alpha _{ex}(\theta )`$, calculated numerically and the strong field limit $`\alpha _{app}(\theta )`$. The most external image appears where $`\alpha (\theta )`$ falls below $`2\pi `$. This happens at $`uu_m=0.003264`$. Here the discrepancy between $`\alpha _{ex}`$ and $`\alpha _{app}`$ is about $`0.06\%`$, which corresponds to an error in the position of the outer image of the order $`0.4\%`$. With such a high accuracy, we are encouraged to take the Schwarzschild strong field limit as the starting point for successive series expansions to evaluate $`b_R`$ in more advanced metrics.

### III.2 Reissner-Nordstrom lensing

Reissner-Nordstrom metric describes the gravitational field of a spherically symmetric massive object endowed with an electric charge $`q`$. The metric functions in standard coordinates are

$`A(x)=1{\displaystyle \frac{1}{x}}+{\displaystyle \frac{q^2}{x^2}}`$ (51)
$`B(x)=\left(1{\displaystyle \frac{1}{x}}+{\displaystyle \frac{q^2}{x^2}}\right)^1`$ (52)
$`C(x)=x^2.`$ (53)
They satisfy the hypotheses required in Sect. II, only when $`q\frac{3}{4\sqrt{2}}`$. However, beyond the critical value $`q=0.5`$, there is no event horizon and causality violations appear Carter ; HawEll . We shall restrict to $`q<0.5`$.

The coefficients $`\alpha `$ and $`\beta `$ are

$`\alpha =\left(2{\displaystyle \frac{3}{x_0}}+{\displaystyle \frac{4q^2}{x_0^2}}\right){\displaystyle \frac{x_0q^2}{x_02q^2}}`$ (54)
$`\beta =\left({\displaystyle \frac{3}{x_0}}1{\displaystyle \frac{9q^2}{x_0^2}}+{\displaystyle \frac{8q^4}{x_0^3}}\right){\displaystyle \frac{x_0\left(x_0q^2\right)^2}{\left(x_0^32q^2\right)^3}}`$ (55)
which reduce to the Schwarzschild coefficients when $`q0`$. From the equation $`\alpha =0`$, we derive the radius of the photon sphere

$$x_m=\frac{3}{4}\left(1+\sqrt{1\frac{32q^2}{9}}\right),$$
(56)
which yields

$`\beta _m`$ $`=[9+32q^2144q^4+512q^6+\sqrt{932q^2}`$ (57)
$`(3+16q^280q^4)\left]\right[8(q4q^3)]^2.`$
The regular term $`b_R`$ cannot be calculated analytically, but we can expand the integrand in (29) in powers of $`q`$ and evaluate the single coefficients. We get

$$b_R=b_{R,0}+b_{R,2}q^2+O(q^4),$$
(58)
where $`b_{R,0}`$ is the value of the coefficient for an uncharged black hole, calculated in the previous subsection and given by Eq. (46); the correction is quadratic in the charge $`q`$ of the black hole. Its coefficient is

$$b_{R,2}=\frac{8}{9}\left\{\sqrt{3}4+\mathrm{log}\left[6\left(2\sqrt{3}\right)\right]\right\}=1.5939.$$
(59)
It is very easy to calculate further terms in the expansion of $`b_R`$, deriving analytical formulae which prove to be very accurate even for large values of $`q`$.

Following Sect. II, we calculate the coefficients for the formula of the deflection angle

$`\overline{a}`$ $`={\displaystyle \frac{x_m\sqrt{x_m2q^2}}{\sqrt{\left(3x_m\right)x_m^29q^2x_m+8q^4}}}`$ (61)
$`\overline{b}`$ $`=\pi +b_R+\overline{a}\mathrm{log}2(x_mq^2)^2`$
$`{\displaystyle \frac{\left[\left(3x_m\right)x_m^29q^2x_m+8q^4\right]}{\left(x_m2q^2\right)^3\left(x_m^2x_m+q^2\right)}}`$
$`u_m`$ $`={\displaystyle \frac{\left(3+\sqrt{932q^2}\right)^2}{4\sqrt{2}\sqrt{38q^2+\sqrt{932q^2}}}}.`$ (62)
In ERT , the coefficients $`a`$ and $`b`$ were calculated numerically. Here, by our general method, we have been able to derive the coefficient $`\overline{a}`$ for the logarithmic divergence exactly and find a formula for $`\overline{b}`$ which is valid up to second order in $`q`$, indicating the way to extend it to an arbitrary order.

We notice that the radius of the photon sphere $`x_m`$ decreases as the charge increases, but becomes imaginary only for $`q=\frac{3}{4\sqrt{2}}>\frac{1}{2}`$; i.e., even when there is no horizon, a photon can be captured by the gravitational field of a hypothetical object with a charge greater than the critical value.

In Fig. 2, we evaluate the deflection angle at $`u=u_m+0.003`$ as a function of $`q`$. The plot shows that, just using the first correction to $`b_R`$, we get an excellent approximation to the deflection angle. In fact, we see that up to $`q=0.3`$, the error in the position of the outer image, calculated using the strong field limit, stays below $`4\%`$.

Finally, in Fig. 3, we plot the coefficients of the strong field limit as functions of $`q`$, calculating the full $`\overline{b}`$ numerically. We see that $`\overline{a}`$ and $`\overline{b}`$ deviate from the corresponding Schwarzschild values as the charge increases. As we shall see in Sect. IV, the strong field limit coefficients are directly connected to the observables. It is then possible, in principle, to distinguish a Reissner-Nordstrom black hole from a Schwarzschild black hole, using strong field gravitational lensing.

### III.3 Janis-Newman-Winicour lensing

The spherically symmetric solution to the Einstein massless scalar equations ($`R_{\mu \nu }=\mathrm{\Phi }_{,\mu }\mathrm{\Phi }_{,\nu }`$, $`\mathrm{\Phi }_{;\mu }^{;\mu }`$) can be written in Janis-Newman-Winicour (JNW) coordinates Brans

$`A(x)=\left(1{\displaystyle \frac{1}{x}}\right)^\gamma `$ (63)
$`B(x)=\left(1{\displaystyle \frac{1}{x}}\right)^\gamma `$ (64)
$`C(x)=\left(1{\displaystyle \frac{1}{x}}\right)^{1\gamma }x^2`$ (65)
$`\mathrm{\Phi }(x)={\displaystyle \frac{q}{2\sqrt{M^2+q^2}}}\mathrm{log}\left(1{\displaystyle \frac{1}{x}}\right)`$ (66)
$`\gamma ={\displaystyle \frac{M}{\sqrt{M^2+q^2}}}.`$ (67)
where all distances are measured in terms of $`x_s=2\sqrt{M^2+q^2}`$ and $`q`$ is the scalar charge of the central object. This metric admits a photon sphere external to the static limit when $`\gamma >\frac{1}{2}`$, i.e. when $`q<M`$. We shall thus restrict our investigations to objects with scalar charge lower than their mass. In Ref. VirEll2 , the gravitational lensing of this object was investigated numerically even when $`q>M`$. In this situation, it was shown that a drastically different and interesting phenomenology shows up.

As in the previous cases, we compute all the coefficients, taking into account that our metric is not written in standard coordinates. The coefficients $`\alpha `$ and $`\beta `$ are

$`\alpha `$ $`=\left(2{\displaystyle \frac{2\gamma +1}{x_0}}\right){\displaystyle \frac{1}{\gamma x_0^{\gamma 1}}}\left[x_0^\gamma \left(x_01\right)^\gamma \right]`$ (69)
$`\beta `$ $`={\displaystyle \frac{\left(2\gamma +1\right)\left(\gamma +1\right)2x_0\left(3\gamma +1\right)+2x_0^2}{\left(x_01\right)^\gamma }}`$
$`{\displaystyle \frac{\left[x_0^\gamma \left(x_01\right)^\gamma \right]^2}{2\gamma ^2x_0^\gamma }},`$
which easily reduce to Eqs. (43) and (44) when $`\gamma =1`$. From the Eq. $`\alpha =0`$, we derive the radius of the photon sphere

$$x_m=\frac{2\gamma +1}{2}$$
(70)
and then

$$\beta _m=\frac{\left[\left(2\gamma +1\right)^\gamma \left(2\gamma 1\right)^\gamma \right]^2}{4\gamma ^2\left(4\gamma ^21\right)^{2\gamma 1}}$$
(71)
In the same way as in Reissner-Nordstrom metric, increasing the charge of the central object, the photons are allowed to get closer to the centre. When $`\gamma =0.5`$, $`x_m=x_s`$ and we should change our coordinates frame to perform any significant study.

Here as well we cannot solve the integral (29) exactly. However, we can expand the integrand in powers of $`(\gamma 1)`$ to get

$$b_R=b_R^00.1199(\gamma 1)+O(\gamma 1)^2$$
(72)
with $`b_R^0`$ given by Eq. (46).

Finally we compute the coefficients of the strong field limit

$`\overline{a}`$ $`=1`$ (74)
$`\overline{b}`$ $`=\pi +b_R+`$
$`+2\mathrm{log}{\displaystyle \frac{\left[\left(2\gamma +1\right)^\gamma \left(2\gamma 1\right)^\gamma \right]^2\left(2\gamma +1\right)}{2\gamma ^2\left(2\gamma 1\right)^{2\gamma 1}}}`$
$`u_m`$ $`={\displaystyle \frac{\left(2\gamma +1\right)^{\gamma +\frac{1}{2}}}{2\left(2\gamma 1\right)^{\gamma \frac{1}{2}}}}`$ (75)
Surprisingly, the leading order coefficient $`\overline{a}`$ is the same as in the Schwarzschild case, in spite of the fact that the radius of the photon sphere has changed.

Once we fix $`uu_m=0.003`$, we see in Fig. 4 that the deflection angle decreases as we increase the charge (decrease $`\gamma `$) until $`\gamma `$ reaches the value $`0.6`$, then the deflection angle increases again. The strong field limit, with $`b_R`$ truncated to first order in $`\gamma `$, at $`\gamma =0.7`$ is precise up to $`4\%`$ in the determination of the outer image.

In Fig. 5, we plot the coefficients of the strong field limit. $`\overline{a}`$ is constant, $`\overline{b}`$ increases as $`\gamma `$ decreases, but $`u_m`$ decreases enough to make the deflection angle decrease at constant $`uu_m`$ until $`\gamma `$ reaches the value 0.6. Comparing with Fig. 3, we can observe that a JNW charge has completely different effects on the strong field limit coefficients than an electric charge and can be identified without confusion.

## IV Obervables in the strong field limit

In Sect. II we have proved that the strong field limit approximation can be used to obtain a simple and reliable formula for the deflection angle, which contains a logarithmic and a constant term. Now we plug the formula (9) into the lens equation and establish direct relations between the position and the magnification of the relativistic images and the deflection angle, calculated according to the strong field limit.

The lens equation in the strong field limit was derived in Ref. BCIS . It reads

$$\beta =\theta \frac{D_{LS}}{D_{OS}}\mathrm{\Delta }\alpha _n.$$
(76)
where $`D_{LS}`$ is the distance between the lens and the source, $`D_{OS}=D_{OL}+D_{LS}`$, $`\beta `$ is the angular separation between the source and the lens, $`\theta `$ is the angular separation between the lens and the image, $`\mathrm{\Delta }\alpha _n=\alpha (\theta )2n\pi `$ is the offset of the deflection angle, once we subtract all the loops done by the photon.

To pass from the deflection angle $`\alpha (\theta )`$ to the offset $`\mathrm{\Delta }\alpha _n`$, we need to find the values $`\theta _n^0`$ such that $`\alpha (\theta _n^0)=2n\pi `$. Solving Eq. (35) with $`\alpha (\theta )=2n\pi `$, we find

$`\theta _n^0={\displaystyle \frac{u_m}{D_{OL}}}\left(1+e_n\right)`$ (77)
$`e_n=e^{\frac{\overline{b}2n\pi }{\overline{a}}}.`$ (78)
The offset $`\mathrm{\Delta }\alpha _n`$ is then found expanding $`\alpha (\theta )`$ around $`\theta =\theta _n^0`$. Letting $`\mathrm{\Delta }\theta _n=\theta \theta _n^0`$, we find

$$\mathrm{\Delta }\alpha _n=\frac{\overline{a}D_{OL}}{u_me_n}\mathrm{\Delta }\theta _n.$$
(79)
The lens equation becomes

$$\beta =\left(\theta _n^0+\mathrm{\Delta }\theta _n\right)+\left(\frac{\overline{a}D_{OL}}{u_me_n}\frac{D_{LS}}{D_{OS}}\right)\mathrm{\Delta }\theta _n.$$
(80)
Now we derive the position of all relativistic images, their magnification and the critical curves of the lens.

The second term in the r.h.s. of Eq. (80) is negligible when compared to the last one (since $`u_mD_{OL}`$). Immediately, we find

$$\theta _n=\theta _n^0+\frac{u_me_n\left(\beta \theta _n^0\right)D_{OS}}{\overline{a}D_{LS}D_{OL}}.$$
(81)
where the correction to $`\theta _n^0`$ is much smaller than $`\theta _n^0`$.

This formula is valid both for the images on the same side of the source and the images on the opposite side. In fact, to find the latter, it is sufficient to take a negative $`\beta `$ in Eq. (81).

Finally, we have expressed the position of the relativistic images in terms of the coefficients $`\overline{a}`$, $`\overline{b}`$ and $`u_m`$. If we manage to determine these coefficients from the observation of the relativistic images, we are rewarded with information about the parameters of the black hole stored in them.

The critical curves are just Einstein rings corresponding to a source perfectly aligned with the lens. Their radius is obtained putting $`\beta =0`$ in (81).

Another important source of information is the magnification of the images which is the inverse of the Jacobian evaluated at the position of the image. For simplicity, we approximate the position of the images by $`\theta _n^0`$, recalling that the correction provided by Eq. (81) is negligible.

$$\mu _n=\frac{1}{\frac{\beta }{\theta }\frac{\beta }{\theta }}|_{\theta _n^0}.$$
(82)
We have

$$\frac{\beta }{\theta }|_{\theta _n^0}=1+\frac{\overline{a}D_{OL}}{u_me_n}\frac{D_{LS}}{D_{OS}}.$$
(83)
where the first term is small compared to the second and can be neglected.

Finally, the magnification is

$$\mu _n=\frac{1}{\left|detJ|_{\theta _n^0}\right|}=\frac{\theta _n^0}{\beta \frac{\beta }{\theta }|_{\theta _n^0}}=e_n\frac{u_m^2\left(1+e_n\right)D_{OS}}{\overline{a}\beta D_{OL}^2D_{LS}},$$
(84)
which decreases very quickly in $`n`$.

Formulae (77) and (84) relate the position and the magnification to the strong field limit coefficients. The successive step is to solve the inversion problem, i.e. we have to find out the most effective way to go back from measured positions and fluxes to the strong field limit coefficients, which carry the information about the nature of the black hole.

The minimum impact parameter can be simply obtained as

$$u_m=D_{OL}\theta _{\mathrm{}}$$
(85)
where $`\theta _{\mathrm{}}`$ represents the asymptotic position approached by a set of images, obtained by Eq. (77) in the limit $`n\mathrm{}`$.

To obtain the coefficients $`\overline{a}`$ and $`\overline{b}`$, we need to separate at least the outermost image from all the others. We shall thus consider the simplest situation where only the outermost image $`\theta _1`$ is resolved as a single image, while all the remaining ones are packed together at $`\theta _{\mathrm{}}`$.

Our observables will thus be

$`s=\theta _1\theta _{\mathrm{}}`$ (86)
$`r={\displaystyle \frac{\mu _1}{\underset{n=2}{\overset{\mathrm{}}{}}\mu _n}},`$ (87)
which respectively represent the separation between the first image and the others, and the ratio between the flux of the first image and the flux coming from all the other images.

The sum of the fluxes of all the set of relativistic images except the first is

$$\underset{n=2}{\overset{\mathrm{}}{}}\mu _n=\frac{u_m^2D_{OS}e^{\frac{\overline{b}}{\overline{a}}}}{\overline{a}\beta D_{OL}^2D_{LS}}\frac{e^{\frac{4\pi }{\overline{a}}}+e^{\frac{2\pi }{\overline{a}}}+e^{\frac{\overline{b}}{\overline{a}}}}{e^{\frac{4\pi }{\overline{a}}}1}.$$
(88)
We can notice that $`e^{\frac{2\pi }{\overline{a}}}1`$ and $`e^{\frac{\overline{b}}{\overline{a}}}`$ is of order one, since $`\overline{a}`$ and $`\overline{b}`$ are of order one too. Using these simple observations, we can simplify our formulae to have

$`s=\theta _{\mathrm{}}e^{\frac{\overline{b}}{\overline{a}}\frac{2\pi }{\overline{a}}}`$ (89)
$`r=e^{\frac{2\pi }{\overline{a}}}.`$ (90)
These two formulae can be easily inverted to give

$`\overline{a}={\displaystyle \frac{2\pi }{\mathrm{log}r}}`$ (91)
$`\overline{b}=\overline{a}\mathrm{log}\left({\displaystyle \frac{rs}{\theta _{\mathrm{}}}}\right).`$ (92)
Finally, just measuring an angular separation and a flux ratio, we are able to reconstruct the full strong field limit expansion of the deflection angle for the observed gravitational lens.

The coefficients of the strong field limit are constrained by the characteristics of the metric to be well precise real numbers. Thus, if strong field gravitational lensing will be detected, comparing the experimental coefficients with the theoretical expectations, calculated according to different models, we are able to identify the nature of the lensing black hole unambiguously.

### IV.1 A numerical example: lensing by the galactic supermassive black hole

It is significant to consider a realistic case where we can discuss the instrumental sensitivity required to detect relativistic images and possibly distinguish between different black holes through the reconstruction of the strong field limit coefficients.

The centre of our galaxy is believed to host a black hole with mass $`M=2.8\times 10^6M_{}`$ Richstone . The lensing of a background source by this supermassive black hole was discussed in detail by Virbhadra & Ellis VirEll . Taking $`D_{OL}=8.5`$kpc as the distance between the sun and the center of Galaxy, they found that the separation between each set of relativistic images with respect to the central lens would be $`\theta _{\mathrm{}}17`$ microarcsecs. In principle, such a resolution is reachable by actual VLBI projects, but we must be aware that the disturbances intrinsic in such observations (mainly due to extinction and emission by accreting matter), would make the identification of the relativistic images very difficult, as already pointed out in Ref. VirEll .

In Tab. 1, following this line, we estimate the quantities we need for a complete strong field limit reconstruction in different situations, starting from a simple Schwarzschild black hole and then going to black holes with different values of the electric charge and JNW charge.

Looking at Tab. 1, we can make some considerations of different order. Indeed, the easiest parameter to evaluate is the minimum impact parameter $`u_m`$, since a microarcsec resolution is reachable in the next years. This information alone, can already distinguish between a Schwarzschild or other types of geometry. In fact, since the total mass and the distance of the black hole are known to a reasonable accuracy (and possibly will be even better fixed in the next years), an $`u_m`$ smaller than predicted would signal that the structure of spacetime close to the central black hole is not described by Schwarzschild solution. On the other hand, if $`u_m`$ is compatible within experimental uncertainties with Schwarzschild case, we would set an upper limit for the parameters describing other black holes, such as an electric or a scalar charge.

In a second extent, to fit all the strong field limit coefficients into any black hole model, we need to separate at least the outermost relativistic image from the others. We see that this can be done only increasing optical resolution at least by two orders of magnitude with respect to actual observational projects. Therefore, with these numbers, it seems that we are forced to wait for further technological developments. However, given the evolution rate of astronomical facilities in the last twenty years, it is not unthinkable that these two orders of magnitude will be reached within a not so far future.

Black holes are also present in the bulge of other galaxies. As far as we have investigated, in the best cases, the instrumental resolutions needed for strong field gravitational lensing are about the same as for the Milky Way black hole. So, the determination of $`u_m`$ by next future observations would possibly confirm or disprove Schwarzschild geometry for several extra-galactic black holes too. The complete determination of the strong field limit coefficients instead remains a long term project, unless more favourable astrophysical situations emerge.

## V Summary

Gravitational lensing is undoubtedly a potential powerful tool for the investigation of strong fields. By general arguments we have shown that the deflection angle diverges logarithmically as we approach the photon sphere. We have drawn a general method to compute the coefficient of the leading order divergent term and the first regular term. When the latter cannot be calculated analytically, we have seen that it can be well approximated by a simple series expansion starting from Schwarzschild spacetime.

We have applied our method to Schwarzschild, Reissner-Nordstrom and Janis-Newman-Winicour black hole, explicitly calculating and plotting the strong field limit coefficients.

Of course, it is possible to apply the strong field limit, in the form given in this paper, to any spherically symmetric metric representing a black hole. In this way, it is possible to compare the gravitational lensing behaviour of these objects in different theories of gravitation. In principle, the extension to non-spherically symmetric and rotating black holes is possible. However, the dependence of the deflection angle on more than one variable can put severe obstacles in the way of analytic solutions of the problem. Nevertheless, this is indeed another important point which needs to be investigated to complete the picture of black hole lensing.

Differences in the deflection angle are immediately reflected on the relativistic images. If the mass and the distance of the lens is known, then the detection of any set of relativistic images would immediately check the Schwarzschild geometry. VLBI should be able to provide an observational answer, if the relativistic images are not hidden behind environmental noise.

Furthermore, if the outermost image is resolved from the others, it is then possible to fully reconstruct the strong field limit coefficients and select a precise black hole model. Our present observational facilities are not so far from the required resolutions, which, for the galactic black hole, are of the order of 0.01 microarcsecs. As a long term project the detection of the outermost image stands as a very interesting, non-trivial challenge for future technology.

Strong field limit represents an important step in the construction of a robust theoretical scheme connecting the gravitational lensing with the strong field properties. By a simple and reliable expansion, it clarifies the whole phenomenology and the differences between various models that should be expected in the appearance of relativistic images. If these so elusive features will be detected, we will finally have a way to effectively discriminate between alternative theories of gravitation and grow in our knowledge of spacetime.

###### Acknowledgements.

I would like to thank Gaetano Lambiase and Salvatore Capozziello for helpful comments and discussions.