format: text

text:
# The Orientation of Jets Relative to Dust Disks in Radio GalaxiesBased partly on observations made with the NASA/ESA Hubble Space Telescope, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

## 1 Introduction

The orientation of a black hole rotation axis relative to its host galaxy rotation axis is an important diagnostic for the structure of the accretion disk and the nuclear region of the galaxy. In the case of spiral galaxies, the simplest picture suggests that they should be aligned, since the gas which fuels the nucleus is expected to originate in the galaxy disk, or if new gas is added to the galaxy, it will end up in the disk. However, Clarke, Kinney & Pringle (1998), Nagar & Wilson (1999) and Kinney et al. (2000, hereafter K00), found that this assumption is wrong for Seyfert galaxies. By comparing the orientation of their radio jets, which are believed to be perpendicular to accretion disks and aligned with the black hole axes, with the major axis of their host galaxies, they showed that the observed distribution can be reproduced by randomly oriented jets.

In the case of radio galaxies the situation is a little different from that of Seyferts. Most of their hosts are ellipticals or S0s and have very little cool gas to feed the nucleus, meaning that the gas is likely to have an external origin, probably from mergers. Birkinshaw & Davies (1985) and Sansom et al. (1987) presented the comparison between the position angle of jets and host galaxy major axes for this class of objects. They found that there is very little or no correlation between the two, a result similar to the one found for Seyfert galaxies. This result is contrasted by the one obtained by Kotanyi & Ekers (1979), who approached the problem using dust lanes to trace the orientation of the gas in these galaxies. They found that radio galaxies have jets aligned nearly perpendicular to these dust lanes, implying a close connection between the gas in these structures and the accretion disk. Möllenhoff, Hummel & Bender (1992) obtained a similar result using ground based data. More recently, van Dokkum & Franx (1995) and de Koff et al. (2000) used HST images to detect dust lanes in radio galaxies. They found a large range of values for the angle between the jets and dust lanes, but suggest that the majority of the jets are still aligned nearly perpendicular to the dust lanes (however, see de Ruiter et al. 2001 for a different result).

Although the previous studies about the orientation of jets in radio galaxies presented interesting results, they faced several problems. For those papers which compared the orientation of jets relative to host galaxies major axes, it would be important to determine if their host galaxies are triaxial or oblate systems. It was shown by Binney (1985) that the intrinsic major axis does not necessarily correspond to the projected one in triaxial galaxies, and consequently the results obtained by Sansom et al. (1987) and Birkinshaw & Davies (1985) could be contaminated to orientation effects. In the case of papers which used the orientation of jets relative to dust lanes, it is assumed that the gas and dust define a preferential plane, thus eliminating the uncertainties of host galaxy projection (Kotanyi & Ekers 1979). However, the inspection of the HST images presented by De Koff et al. (2000) and Capetti et al. (2000) shows that most of the dust lanes in radio galaxies are patchy and not very well organized, suggesting that they may not be dynamically relaxed systems, thus not directly related to the feeding of the nucleus. It is extremely difficult, and even impossible in some of the cases, to determine the orientation and measure a position angle of the dust in most of these galaxies.

In this letter we address the question of the orientation of jets relative to dust structures in radio galaxies using a new approach. Instead of using all radio galaxies with dust lanes, we use only those with well defined dust disks in their nuclei. We also improve on previous analyzes by using the statistical techniques devised by K00 (see below).

## 2 Data and Measurements

In order to avoid the uncertainties faced by previous works, we decided to use only radio galaxies which present a nuclear dust disk. These disks are relaxed structures, settled in the gravitational potential of the galaxy, for which we can determine the position angle of the major axis, as well as the inclination relative to the line of sight, under the assumption that they are circular. Another advantage of using only dust disks is the fact that it allows us to use the same technique developed by K00 for Seyfert galaxies and obtain information about the 3-dimensional distribution of the jets. This technique determines, in a statistical way, which distributions of angles $`\beta `$, the angle between the jet and dust disk rotation axis in 3 dimensions, are the best representations to the data. Most of the previous studies in this field used only the difference between the position angles of dust lanes and jets for their analysis, which are 2-dimensional projections on the plane of the sky and do not give information about the 3 dimensional distribution.

Our sample was selected from a list of radio galaxies with z$`0.1`$ available in NED. A source was considered a radio galaxy if its host morphology is elliptical or S0, with optical absolute magnitude M$`{}_{R}{}^{}21`$ and radio power P$`{}_{1.4GHz}{}^{}10^{22}`$W Hz-1. This list was cross correlated with the HST archive, and the galaxies for which we could find broad band images were inspected for dust disks. We found a total of 20 galaxies with dust disks which can be used for our analysis. Notice that this sample is not complete, like the 3C sample, since it involves galaxies from different surveys (3C, 4C, Bologna, PKS), discovered based on observations at different frequencies, with different sensitivities, but it includes all radio galaxies with known dust disks. If we had used only the galaxies from the 3C survey, the sample would be reduced to too small a number to allow us to perform our analysis. Of the 20 galaxies in the sample, 13 already had measurements of disk major axis position angle and inclination published by other authors. For the remaining galaxies we carried out these measurements by overplotting different ellipses on top of the images, centered at the nucleus, and determined by eye which one gave the best reproduction to the observed dust disk. The uncertainties in these measurements are of the order of 5, similar to the ones involved in the measurements of the position angles of the jets. The smallest scale on which the dust disks can be traced is 0.1′′, the HST resolution, which corresponds to sizes between 20 pc and 150 pc, depending on the distance of the galaxy (we use H$`{}_{0}{}^{}=75`$ km s-1 Mpc-1). The radio measurements from VLA, and in some cases also from VLBI, were obtained from the literature, measured directly on published radio images, and, in some cases, measured on images obtained reducing VLA archival data.

Galaxies like 3C 40 were excluded from the sample, because, even though we could see a dust disk, the galaxy does not present a jet, only the hotspots far away from the nucleus, making it impossible to obtain a robust estimate of the position angle of the extended radio emission. Galaxies like Centaurus A and NGC7626 were not included in our study because their dust disks are warped and we are not able to obtain a reliable value for their inclinations. Similarly, cases like 3C 264, which show a face-on dust disk (Sparks et al. 2000), also had to be eliminated from the sample, since in these cases it is impossible to constrain the position of the jet relative to the disk. The fact that Sparks et al. (2000) find round disks supports our assumption that they are intrinsically round, and will be discussed in detail below.

In Table 1 we present the sample and the measurements. The first column gives the names of the galaxies, followed by their radial velocities in column 2. Columns 3, 4, 5 and 6 present the position angles of the dust disks major axes, the host galaxies major axes, the jets on VLBI scales and on the VLA scales, respectively. Column 7 presents $`\delta `$, the difference between the position angle of the jets measured on VLA images and the dust disks major axes, column 8 gives the inclination ($`i`$) of the dust disks relative to the line of sight and column 9 gives the minimum angle between the jets and the dust disks rotation axes ($`\beta _{min}`$). The morphological types of the host galaxies are given in column 10, and the references from which the measurements were obtained are presented in column 11.

## 3 Results

We first compare the orientation of the dust disks relative to their host galaxies. Figure 1 presents the difference between the position angle of the host galaxies major axes and dust disks major axes for our sample. In most of the cases the dust disk major axis is closely aligned with the host galaxy major axis, with the exception of 3C 402N and 4C–03.43, where they are aligned close to the minor axis. Analyzing these two galaxies in detail we find that 3C 402N is interacting with other galaxies (de Juan et al. 1994), and the isophotes close to the nucleus change position angle by 90 relative to the outer ones, so the dust disk is in fact aligned with the major axis in the nuclear region. In the case of 4C–03.43, we find that the host galaxy is almost face-on which makes it very difficult to accurately measure the position angle of the major axis.

A naive analysis of the result from Figure 1, ignoring projection effects, might suggest that these are triaxial galaxies, since according to Schwarzschild (1979) only orbits around the major and minor axes are stable in these systems. However, it was shown by Binney (1985), Franx, Illingworth & de Zeeuw (1991) and Ryden (1992) that, due to projection effects, the observed major and minor axes, projected on the plane of the sky, do not correspond to the intrinsic ones in triaxial galaxies. In contrast, the distribution seen in Figure 1 is very narrow, with most of the dust disks aligned very close to the host galaxy major axis, similar to the expected distribution for oblate systems. This result indicates that our galaxies either have disks or are oblate axisymmetric systems. In fact, 30% of them are classified as S0 or E/SO (Table 1), indicating that they have disks. The remaining galaxies are classified as ellipticals and we consider them to be oblate systems. Based on this fact we assume that the dust disks are intrinsically circular, and use this assumption to derive their inclinations relative to the line of sight. This agrees with the results from Sparks et al. (2000) who present several radio galaxies with circular dust disk, seen face-on.

The above result is similar to the one obtained by Tran et al. (2001) for normal elliptical galaxies, where they also found that their nuclear dust disks are aligned close to the host galaxies major axes, indicating that they are oblate systems. A complicating factor in the assumption that the dust disks are intrinsically circular have been recently presented by Andersen & Bershady (2002), who showed that disk galaxies are not always exactly round, but have slightly elongated shapes, with an average ellipticity e$`=0.06`$. If we assume that the intrinsic ellipticity of all our dust disks is e$`=0.06`$, we get that in 75% of the cases, those with $`i60^{}`$, the inclinations would be reduced by $`1^{}2^{}`$, which is smaller than the uncertainty in the measurements.

Figure 2 presents the distribution of $`\delta `$’s, obtained using the jet directions measured on VLA scales. We prefer to use the measurements based on VLA images because they are available for all galaxies in the sample. Although the VLBI measurements would trace the jets in regions much closer to the nucleus, they are available only for 8 galaxies. The position angles obtained from VLBI differ by no more than 7 from the VLA results, and making use of the VLBI measurements does not significantly affect our conclusions. The inspection of this Figure shows that there is a large range of $`\delta `$’s in our sample from 19 to 89, suggesting that there is not a good alignment of the jets perpendicular to the dust disks. As explained above, the $`\delta `$ distribution is a 2-dimensional projection on the plane of the sky, and we are actually interested on the orientation of jets in 3-dimensions, distribution of angles $`\beta `$, which is obtained using the techniques developed by K00.

We give here a brief description of the geometry of the problem, and a more detailed description can be seen in K00. We assume the nucleus to be centered at the origin of the reference frame, with the dust disk lying in the XY-plane, the major axis along the X-axis, the rotation axis along the Z-axis, and the dust disk inclined relative to the line of sight by an angle $`i`$. The dust disk major axis and the jet are projected onto the sky plane, and the difference between their position angles is $`\delta `$. The angle between the jet and the galaxy rotation axis is $`\beta `$. It was shown by Clarke et al. (1998) that for a given galaxy with a measured pair of angles $`i`$ and $`\delta `$, the jet is determined to lie on two great circles, which are mirror reflections of each other. These great circles are drawn on a sphere centered at the origin of the system, the nucleus. A particular property of the geometry of this problem is the fact that all points in a great circle lie at angles $`\beta \beta _{min}`$, meaning that for a given pair ($`i`$,$`\delta `$) the jet will never be at an angle closer than $`\beta _{min}=\mathrm{sin}^1(\mathrm{sin}i\mathrm{cos}\delta )`$ relative to the disk rotation axis, but can always be as large as 90, along the equator. In order words, for a dust disk of inclination $`i`$, there is a large range of $`\beta `$ angles which can reproduce an observed value of $`\delta `$.

The values of $`\beta _{min}`$ for our sample are presented in Table 1 and the cumulative distribution of these angles is presented in Figure 3. This distribution varies smoothly and shows a large range of values. In the case of some radio galaxies, like NGC 315, the jet could be very well aligned with the dust disk rotation axis, with $`\beta _{min}=1^{}`$. However, in most of the cases the jet has to be misaligned by at least 20, with 30% of the sample having $`\beta _{min}>30^{}`$. The largest value of $`\beta _{min}`$ corresponds to 53 for UGC 12064. The uncertainties in the values of $`\beta _{min}`$, calculated using propagation of errors in $`i`$ and $`\delta `$, is of the order of 4$`{}_{}{}^{}5^{}`$.

We would like to be able to determine the value of $`\beta `$ for each one of the galaxies. Unfortunately, with the exception of a small number of cases (Pringle et al. 1999), we cannot determine this value directly, but rather a range of acceptable values. The distribution of $`\beta _{min}`$’s already indicates that the jets do not lie preferentially perpendicular to the dust disks. In order to improve on this result we used the method described by K00 to test whether a homogeneous distribution of jets over polar caps of different sizes, in the range 0$`{}_{}{}^{}\beta \beta _2`$, can represent the observed $`\delta `$ distribution. This technique assumes a homogeneous distribution of $`\beta `$’s and uses, for each one of the galaxies, the information about $`i`$ and $`\delta `$. This distribution is integrated between 0$`\delta \delta _k`$, where $`\delta _k`$ is the observed value of $`\delta `$ for the galaxy, to obtain the centile point $`c_k`$. The final distribution of $`c_k`$’s is then compared with a homogeneous distribution in the interval  using the Kolmogorov-Smirnov test (see K00 for a detailed explanation of the technique).

We tested which distributions of jets can represent the observed data, by varying $`\beta _2`$ between 55 and 90. The lower limit (55) was chosen based on the fact that the largest value of $`\beta _{min}`$ is $`53^{}`$, so any distribution of jets with $`\beta _2<53^{}`$ would not include all the galaxies available in the sample. We find that polar caps with $`\beta _2`$ between $`55^{}`$ and $`77^{}`$, which correspond to regions between 40% and 78% of the parameter space, respectively, are good representations to the data. The KS test shows that the hypothesis that the observed data is drawn from a homogeneous distribution of jets over a polar cap with $`\beta _2=55^{}`$ is consistent at the 94% level, while a polar cap with $`\beta _2=77^{}`$ is consistent at the 5% level. For polar caps with $`\beta _2>77^{}`$ the KS test gives increasingly smaller probabilities of the models agreeing with the observations, reaching a value of 0.9% for jets distributed over the entire sphere ($`0^{}\beta 90^{}`$). Loosely speaking, our finding is that, if the jets are uniformly distributed over angles between 0 and $`\beta _2`$ relative to the disk axis, then at the 2$`\sigma `$ level $`\beta _2`$ lies between 55 and 77.

## 4 Discussion

The above results clearly show that jets are not closely aligned perpendicular to dust disks. They can be represented by homogeneous distributions over large polar caps ($`55^{}\beta _277^{}`$), encompassing 40% to 78% of the parameter space. The smaller polar caps agree in part with the suggestion by Sparks et al. (2000) that jets are aligned at 30$`{}_{}{}^{}40^{}`$ from the rotation axis, but a larger parameter space cannot be ruled out by the current data. The models show a probability between 5% and 0.9% that the jets uniformly populate a region extending more than $`77^{}`$ from the pole. This implies a deficit of jets that are closely aligned along dust disks, as one might expect since they would either destroy the dust disk, or not be able to propagate through it.

The reason why other papers found a better alignment between dust structures (lanes and disks) and radio jets is not very clear. Some of the previous studies were based on small samples (less than 10 galaxies), where selection effects and small number statistics may have influenced the results. Something else which might have influenced the results presented by other papers is the fact that most of them based their conclusions only on the difference between the position angles of the jets and dust structures. As mentioned above, these measurements are only the projection of the intrinsic 3 dimensional distribution on the plane of the sky. They can have a much tighter distribution than the intrinsic one, creating the impression of a better alignment perpendicular to the dust lanes.

As is also the case in Seyfert galaxies, the result we obtained contradicts the simplest expectations that the black hole should have its angular momentum vector aligned to that of the host galaxy. In their study, K00 presented a series of mechanisms to explain how the observed misalignment could happen, some of which can be used to explain radio galaxies. Misaligned inflow of gas towards the nucleus of the galaxy could, in principle, explain the results. Since we are dealing with early type galaxies, the origin of gas to feed the nucleus is probably external from the galaxy, from mergers with gas rich galaxies. Under this scenario, the observed jet would be the result of a previous merger, while the observed dust disk is the result of a more recent one, which had a different orientation, and will produce jets oriented in a different direction. There are several evidence which point against this scenario. First, we showed that the dust disks are closely aligned with the host galaxy major and minor axes, meaning that they had enough time to settle into the gravitational potential of the galaxy and are not the result of recent accretion. Second, the comparison between the orientation of jets on small (VLBI) and large (VLA) scales also indicates that their directions are not regulated by mergers. For 8 of the galaxies we found high resolution VLBI continuum data, which show an alignment better than 7 with the VLA scale jets, showing their long stability over time, from 1 pc to 10 kpc scales.

Since the feeding of the nucleus by a misaligned inflow of gas does not seem to explain the observations, the misalignment has to occur between the inner region of the dust disk, which corresponds to dimensions between 20 pc and 150 pc for the galaxies studied here, and the point where the jet is produced (a few Schwarzshild radii). This can be achieved by the warping of the accretion disk by self-irradiation instability (Pringle 1996, 1997; Maloney, Begelman & Pringle 1996), by the Bardeen-Petterson effect in the case the galaxy has a rapidly rotating black hole, or by a misaligned gravitational potential in the nuclear region of the galaxy. However, some of these mechanisms may present some problems in radio galaxies.

Warping of the accretion disk by self-irradiation instability, possibly enhanced by mhd wind warping, could explain the observed misalignment. In this case the black hole should not be rotating rapidly, as suggested by Blandford (1990) to be the case in radio galaxies. According to Pringle (1997), a rapidly rotating black hole stabilizes the inner region of the accretion disk and can inhibit the self-irradiation instability.

Warping of the accretion disk by the Bardeen-Petterson effect can explain misaligned dust disk and jet axes if the black hole is rotating rapidly, with a rotation axis different from that of the galaxy and the dust disk. This can result from the merger between two black holes (Wilson & Colbert 1995), originated from a galaxy merger, but the time scale should be long enough to allow the gas and dust to settle in the galaxy potential. It is also possible that the black hole with misaligned axis relative to the host galaxy was created by a previous merger, while the gas and dust we see are the result of a more recent interaction. Here the rapidly rotating black hole forces the precession and alignment of the inner regions of the accretion disk with its spin axis. One problem with this scenario, pointed out by Rees (1978) and more recently by Natarajan & Pringle (1998), is that the black hole exerts a torque on the accretion disk, which will also exert a torque on the black hole, eventually forcing the alignment between the black hole and the source of gas feeding it, in this case the dust disk. The timescale for this alignment is somewhat uncertain, but if it is relatively fast it might contradict the fact that we observe the jets to be aligned from the VLBI to the VLA scales.

The warping of the accretion disk could also be due to a misaligned gravitational potential around the nucleus. Since black holes seem to be surrounded by a nuclear star cluster (Gebhardt et al. 1996; Faber et al. 1997; van der Marel 1999), if the star cluster is axisymmetric and has its symmetry axis misaligned relative to that of the galaxy, it would make the disk precess and force its alignment with the symmetry axis of the cluster. This effect may not be very important in radio galaxies. Considering that their black holes have masses of 108 to 109 M, their radius of influence is of the order of 100 pc. A globular cluster would need to have a similar radius and be axisymmetric, in order to induce the precession of the accretion disk.

We would like to thank J. Ulvestad and R. Antonucci for helpful comments. JEP is grateful for continued support from the STScI visitor program. This work was partially supported by the NASA grant AR-8375.01-97A. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jep Propulsion Laboratory, Caltech, under contract with NASA.


text:
# Enhancement of the electron electric dipole moment in gadolinium 3+

## I introduction

A violation of the combined symmetry of charge conjugation (C) and parity (P) has been discovered in the decay of the $`K^0`$ meson about 40 years ago . The exact origin of this symmetry violation remains an enigma, although the standard model of electroweak interactions can describe these processes phenomenologically. It has also been proposed by Sakharov  that the matter-antimatter asymmetry observed in our universe could have arisen from a CP-violating interaction active at an early stage of the big bang. The CP-violation implies a time-reversal (T) asymmetry and hence violation of the combined TP-symmetry, because there are strong reasons to believe that the combined CPT-symmetry should not be violated . An electric dipole moment (EDM) of a system in a stationary quantum state indicates a violation of the TP-symmetry; this is why searches for EDMs of elementary particles, atoms and molecules are a very important approach to the studies of violations of fundamental symmetries . In the present work we concentrate on the EDM of the electron $`d_e`$.

At present the best limitation on $`d_e`$ comes from the Berkeley experiment with an atomic thallium beam , $`d_e<1.610^{27}e\mathrm{cm}`$. There are some ideas for improving the sensitivity. One way of improvement is an experiment with metastable levels of PbO molecules . A breakthrough could be achieved in solid state experiments. This idea was suggested by Shapiro in 1968 . The application of strong electric fields to electrons bound within a solid would align the electric dipole moments of these electrons. This should lead to a simultaneous alignment of the electron spins; the magnetic field arising from this alignment could be detected experimentally. An experiment of this kind has been performed with nickel-zinc ferrite . However, due to experimental limitations, the result was not very impressive. Interest to the approach has been renewed recently due to suggestions of Lamoreaux  and Hunter  to perform similar experiments with gadolinium gallium garnet, Gd3Ga5O12, and gadolinium iron garnet, Gd3Fe5O12, employing new experimental techniques. First estimates of sensitivity promise to improve the current upper limit on the electron EDM by at least three orders of magnitude, depending on the experimental setup. A thorough analysis of the problem requires the calculation of the enhancement coefficient for the electron EDM in the Gd3+ ion. We perform this calculation in the present work.

## II Single particle contribution

The Gd3+ ion has a nucleus with charge $`Z=64`$ and 61 electrons. The 7 electrons in the outer shell occupy the $`4f`$ orbitals. So the shell is half filled, and hence the total orbital angular momentum is zero, $`L=0`$, and total spin $`S=7/2`$. Let us consider the state with maximum z-projection of the spin. In the representation of second quantization the ground state wave function is of the form

$$|gs=f_3^{}f_2^{}f_1^{}f_0^{}f_1^{}f_2^{}f_3^{}|0,$$
(1)
where $`f_{m\sigma }^{}`$ is the creation operator for a $`4f`$ electron with spin $`\sigma `$ and z-projection of orbital angular momentum $`l_z=m`$. The TP-odd interaction of the electron EDM with the electric field $`𝑬`$ is of the form, see e.g. Ref.

$$V_d=d_e\gamma _0𝚺\mathbf{}𝑬,$$
(2)
where $`\gamma _0`$ and $`𝚺=\gamma _0\gamma _5𝜸`$ are Dirac $`\gamma `$-matrices. Because of Schiff’s theorem  it is crucially important to account for very complex many-body screening effects, when working with the Hamiltonian (2); technically this means that the many-body perturbation theory practically is not convergent. The standard way  to avoid this complication is to split the Hamiltonian into two terms: $`V_d=d_e\gamma _0𝚺\mathbf{}𝑬=d_e𝚺\mathbf{}𝑬d_e(\gamma _01)𝚺\mathbf{}𝑬`$. Then due to Schiff’s theorem  the contribution of the first term to the EDM of the ion is identically zero, so one can reduce the interaction

$$V_dV_d^r=d_e(\gamma _01)𝚺\mathbf{}𝑬.$$
(3)
Perturbation theory with this operator is reasonably convergent.

In leading order of single particle perturbation theory the EDM of the ion $`D_{\mathrm{sp}}`$ is given by diagrams shown in Fig.1.

The corresponding formula reads

$$D_{\mathrm{sp}}=2\underset{mn}{}\frac{4f_m|er\mathrm{cos}\theta |nn|V_d^r|4f_m}{E_{4f}E_n}=K_{\mathrm{sp}}d_e,$$
(4)
where the EDM enhancement coefficient is given by

$$K_{\mathrm{sp}}=2\underset{mn}{}\frac{4f_m|(r/a_B)\mathrm{cos}\theta |nn|(ea_B/d_e)V_d^r|4f_m}{E_{4f}E_n}.$$
(5)
Here $`e`$ is the electron charge and $`a_B`$ is the Bohr radius. We first consider $`5d`$ and $`5g`$ intermediate states, as the contribution of other states is much less. The E1 matrix elements are of the form

$`4f_m|(r/a_B)\mathrm{cos}\theta |5d_m`$ $`=`$ $`\sqrt{{\displaystyle \frac{9m^2}{35}}}r_{4f,5d},`$ (6)
$`4f_m|(r/a_B)\mathrm{cos}\theta |5g_m`$ $`=`$ $`\sqrt{{\displaystyle \frac{16m^2}{63}}}r_{4f,5g},`$ (7)
where $`r_{ik}`$ is the E1-transition radial integral expressed in atomic units (Bohr radius),

$$r_{ik}=_0^{\mathrm{}}R_i(r)R_k(r)r^3dr.$$
(8)
$`R_i(r)`$ is the radial wave function of the corresponding orbital, $`_0^{\mathrm{}}R_i^2(r)r^2𝑑r=1`$.

The matrix element of $`V_d^r`$ has been calculated earlier, see e.g. Ref.. In jj coupling scheme, and using a semiclassical approximation, it reads

$$n^{}jl\pm 1|(ea_B/d_e)V_d^r|njl=\frac{4Z^3\alpha ^2}{\gamma (4\gamma ^21)}\frac{Z_{\mathrm{eff}}^2}{(\nu ^{}\nu )^{3/2}}E_0,$$
(9)
where $`\alpha =1/137.036`$ is the fine structure constant, $`\gamma =\sqrt{(j+1/2)^2(Z\alpha )^2}`$, $`E_0=2\mathrm{Ry}=27.2\mathrm{eV}=2\times 109,737\mathrm{cm}^1`$ is the atomic energy unit, $`Z_{\mathrm{eff}}=4`$ is the effective charge for electron motion at large distances, and finally the effective principal quantum numbers $`\nu `$ and $`\nu ^{}`$ are defined by the electron energy with respect to the ionization limit, $`E=Z_{\mathrm{eff}}^2E_0/(2\nu ^2)`$. Our numerical calculations in Hartree-Fock approximation agree with the semiclassical result (9) within ten per cent. The ground state wave function (1) is given in LS coupling scheme, therefore we have to rewrite the operator $`V_d^r`$ in this representation. A simple calculation gives

$`5d_m|(ea_B/d_e)V_d^r|4f_m`$ $`=`$ $`\sqrt{{\displaystyle \frac{9m^2}{35}}}5d_{5/2}|(ea_B/d_e)V_d^r|4f_{5/2},`$ (10)
$`5g_m|(ea_B/d_e)V_d^r|4f_m`$ $`=`$ $`\sqrt{{\displaystyle \frac{16m^2}{35}}}5g_{7/2}|(ea_B/d_e)V_d^r|4f_{7/2}.`$ (11)
Combining equations (5),(6),(9), and (10) and performing the summation over $`m`$, we find the following expression for the single particle contribution to the enhancement coefficient

$$K_{\mathrm{sp}}=8Z^3\alpha ^2Z_{\mathrm{eff}}^2\left(\frac{r_{4f,5d}}{\gamma _{5/2}(4\gamma _{5/2}^21)(\nu _{4f}\nu _{5d})^{3/2}}\frac{1}{(E_{4f}E_{5d})}+\frac{4}{3}\frac{r_{4f,5g}}{\gamma _{7/2}(4\gamma _{7/2}^21)(\nu _{4f}\nu _{5g})^{3/2}}\frac{1}{(E_{4f}E_{5g})}\right).$$
(12)
Note that from here on we skip the atomic energy unit $`E_0`$ in all equations, assuming that all energies are expressed in units of $`E_0`$. To find values of the parameters that appear in this equation, we have performed a Hartree-Fock calculation for the Gd3+ ion. It is known that a Hartree-Fock calculation for an open shell is not a uniquely defined procedure. In our calculation we used the following averaging: we assumed that the shell is fully occupied, but the occupation number of each single particle $`f`$-orbital is 1/2. This is a crude approximation, it gives reasonable values of energy levels and radial integrals, but one cannot rely on this calculation as far as energy splittings are concerned. As a result of the calculation we obtained the following values of energies and radial integrals:

$`E_{4f}=1.65E_{5d}=1.20E_{5g}=0.32`$ (13)
$`r_{4f,5d}=0.63r_{4f,5g}=0.088`$ (14)
Where available, experimental data should be used, but unfortunately the experimental data on this particular ion is scarce. Only the value of the $`4f`$ energy level (ionization limit) is known, see Ref. ,

$$E_{4f}=355000\mathrm{cm}^1=1.62E_0.$$
(15)
It agrees well with (13). From (15) and (13) one finds $`\nu _{4f}=2.22`$, $`\nu _{5d}=2.58`$ and $`\nu _{5g}=5.0`$.

The most important is the $`E_{4f}E_{5d}`$ energy splitting. The accuracy of the present Hartree-Fock calculation is not sufficient to determine this splitting. There is reliable experimental data for the energy levels of Gd2+, Gd1+, and Eu2+, see Ref.. Naive extrapolation of the splitting from these ions gives

$$A:E_{5d}E_{4f}40,000\mathrm{cm}^10.18E_0.$$
(16)
On the other hand there is experimental data for Gd3+  that indicates

$$B:E_{5d}E_{4f}100,000\mathrm{cm}^10.45E_0.$$
(17)
Unfortunately Ref. does not contain identifications of all the possible levels, therefore we cannot quite rely on the data. For this reason we will present two estimates of $`K`$: one for the case (A), see (16), and another for the case (B), see (17). Experimental and/or theoretical determination of the $`E_{5d}E_{4f}`$ splitting would be the most important to improve the accuracy of the calculation of the enhancement coefficient $`K`$.

Substituting values of the parameters listed above into eq. (12) we find the contribution of the $`5d`$ intermediate state to the EDM enhancement coefficient:

$`A:K_{5\mathrm{d}}=4.5,`$ (18)
$`B:K_{5\mathrm{d}}=1.8`$ (19)
The contribution of the $`5g`$ intermediate state is very small. We also estimate the contribution of the higher $`d`$-levels (mainly the continuous spectrum) as

$$A,B:K_{\mathrm{nd},\mathrm{n}>5}1.$$
(20)
Altogether this gives the following value of the single particle contribution to the EDM enhancement coefficient

$`A:K_{\mathrm{sp}}=5.5,`$ (21)
$`B:K_{\mathrm{sp}}=2.8`$ (22)
## III Many-body corrections

In the situation with cesium or any other atom with valent $`s`$\- or $`p`$-electrons  the single particle estimate is satisfactory. However here we have $`f`$-electrons that have a very small wave function in the vicinity of the nucleus. As a result the single particle contribution is strongly suppressed. Technically this suppression is reflected in eq. (9); the matrix element is proportional to $`1/\gamma ^3`$, and for $`f_{5/2}`$-electrons $`\gamma 3`$, while for $`s_{1/2}`$\- or $`p_{1/2}`$-electrons $`\gamma 0.88`$. Therefore it is very important to estimate the many-body corrections. The leading many-body corrections to the EDM enhancement coefficient $`K_{\mathrm{mb}}`$ are shown in diagrams Fig.2, where the wavy line denotes residual Coulomb interaction

$$V_C=\frac{1}{|𝒓_𝒊𝒓_𝒋|}=\underset{k=0}{\overset{\mathrm{}}{}}\underset{q=k}{\overset{k}{}}\frac{4\pi }{2k+1}\frac{r_<^k}{r_>^{k+1}}Y_{kq}^{}(\mathrm{\Omega }_i)Y_{kq}(\mathrm{\Omega }_j).$$
(23)
We only account for diagrams with momentum of the Coulomb quantum $`k`$ not higher than 2. The contribution of each diagram from Fig.2 must be doubled, because an opposite order of operators is also possible.

The formula for the many-body contribution to the enhancement coefficient due to the diagrams in Fig.2 reads

$`K_{\mathrm{mb}}=2{\displaystyle \underset{m}{}}`$ $`({\displaystyle \frac{4f_m|(r/a_B)\mathrm{cos}\theta |5d_m5p_{1/2}|5d_m|V_C|6s|4f_m6s|(ea_B/d_e)V_d^r|5p_{1/2}}{\left(E_{4f}E_{5d}\right)\left(E_{5p}E_{6s}\right)}}`$ (26)
$`+{\displaystyle \frac{5s|(ea_B/d_e)V_d^r|6p_{1/2}4f_m|(r/a_B)\mathrm{cos}\theta |5d_m6p_{1/2}|5d_m|V_C|5s|4f_m}{\left(E_{5s}E_{6p}\right)\left(E_{4f}+E_{5s}E_{5d}E_{6p}\right)}}`$
$`+{\displaystyle \frac{4f_m|(r/a_B)\mathrm{cos}\theta |5d_m5s|(ea_B/d_e)V_d^r|6p_{1/2}6p_{1/2}|5d_m|V_C|5s|4f_m}{\left(E_{4f}E_{5d}\right)\left(E_{4f}+E_{5s}E_{5d}E_{6p}\right)}}).`$
We have only included $`sp_{1/2}`$ matrix elements of $`V_d^r`$ and only the intermediate states involving $`5d`$ electrons which give the main contribution. All the diagrams in Fig.2 are exchange ones, this is why the sign in eq.(26) is negative.

Matrix elements of the Coulomb interactions are of the form

$`5p_{1/2}|5d_m|V_C|6s|4f_m`$ $`=`$ $`{\displaystyle \frac{1}{5}}\sqrt{{\displaystyle \frac{9m^2}{35}}}F^{(2)}(5p,4f;5d,6s)`$ (27)
$`6p_{1/2}|5d_m|V_C|5s|4f_m`$ $`=`$ $`{\displaystyle \frac{1}{5}}\sqrt{{\displaystyle \frac{9m^2}{35}}}F^{(2)}(6p,4f;5d,5s)`$ (28)
where $`F^{(k)}(i,j;h,l)`$ is the usual Coulomb radial integral,

$$F^{(k)}(i,j;h,l)=_0^{\mathrm{}}_0^{\mathrm{}}\frac{r_<^k}{r_>^{k+1}}R_i(r_1)R_j(r_1)R_h(r_2)R_l(r_2)r_1^2r_2^2dr_1dr_2,$$
(29)
expressed in atomic units.

Using matrix elements (6), (9) and (27), and performing the summation over $`m`$ in (26), we obtain the following expression for the leading many-body correction to the enhancement coefficient

$`K_{\mathrm{mb}}`$ $`=`$ $`{\displaystyle \frac{8Z^3\alpha ^2Z_{\mathrm{eff}}^2}{5\gamma _{1/2}(4\gamma _{1/2}^21)}}r_{4f,5d}\{{\displaystyle \frac{F^{(2)}(5p,4f;5d,6s)}{(\nu _{6s}\nu _{5p})^{3/2}}}{\displaystyle \frac{1}{(E_{4f}E_{5d})}}{\displaystyle \frac{1}{(E_{5p}E_{6s})}}`$ (31)
$`+{\displaystyle \frac{F^{(2)}(6p,4f;5d,5s)}{(\nu _{5s}\nu _{6p})^{3/2}}}{\displaystyle \frac{1}{(E_{4f}E_{5d}+E_{5s}E_{6p})}}({\displaystyle \frac{1}{(E_{4f}E_{5d})}}+{\displaystyle \frac{1}{(E_{5s}E_{6p})}})\}`$
Like the parameters presented in (13) the energy levels and radial integrals in this formula have been calculated using a Hartree-Fock method

$`E_{5s}=3.15E_{5p}=2.26E_{6s}=1.03E_{6p}=0.84`$ (32)
$`F^{(2)}(5p,4f;5d,6s)=0.028F^{(2)}(6p,4f;5d,5s)=0.023.`$ (33)
The corresponding effective principal quantum numbers are $`\nu _{5s}=1.59`$, $`\nu _{5p}=1.88`$, $`\nu _{6s}=2.78`$, $`\nu _{6p}=3.09`$. Using the $`E_{5d}E_{4f}`$ energy splitting given in (16) and (17) and substituting all the parameters into eq. (31), we find the many-body correction due to the $`6s`$-, $`6p`$-, and $`5d`$\- intermediate states, see Fig.2

$`A:K_{\mathrm{mb1}}0.7,`$ (34)
$`B:K_{\mathrm{mb1}}0.3.`$ (35)
One should also perform the summation over higher $`s`$-, $`p`$-, and $`d`$\- intermediate states. We estimate this contribution as

$$A,B:K_{\mathrm{mb2}}0.2$$
(36)
Combining (34) and (36) one finds the leading many-body correction to the enhancement coefficient

$`A:K_{\mathrm{mb}}0.9,`$ (37)
$`B:K_{\mathrm{mb}}0.5.`$ (38)
It is substantially smaller than the single particle contribution (21), so the many-body perturbation theory is convergent. The final result for the electron EDM enhancement coefficient, $`K=K_{\mathrm{sp}}+K_{\mathrm{mb}}`$, reads

$`A:K6.4,`$ (39)
$`B:K3.3.`$ (40)
We recall that the case A corresponds to the energy splitting (16), and the case B corresponds to (17).

## IV conclusion

We have calculated the electron EDM enhancement coefficient in the Gd3+ ion. The single particle contribution as well as the leading many-body corrections have been taken into account. The result is $`K4.9\pm 1.6`$. The main reason for such a large uncertainty lies in the unknown energy splitting $`E_{4f}E_{5d}`$. Experimental and/or theoretical determination of the splitting in Gd3+ would be the most important step for improving the accuracy of the enhancement coefficient calculation.

###### Acknowledgements.

We are grateful to S. K. Lamoreaux, L. R. Hunter, and D. Budker, who attracted our attention to the problem. We would also like to thank them, as well as W. R. Johnson, for stimulating discussions and interest in the work.