text:
# The surface carbon and nitrogen abundances in models of ultra metal-poor stars

## 1 Introduction

Spectroscopic observations of ultra metal-poor (UMP) stars have disclosed that carbon and nitrogen are significantly overabundant with respect to Fe, compared to typical halo stars which show \[C/Fe\]$``$\[N/Fe\]$``$0.0. According to Beers (1999) and Rossi et al. (1999) about 10% of the stars with \[Fe/H\]$``$$``$2.5 display \[C/Fe\]$``$1.0; this fraction rises to about 25% when \[Fe/H\]$``$$``$3.0. Also N appears to be overabundant by the same amount, while $`\alpha `$-elements show a ratio with respect to Fe that is typical of halo stars, i.e. \[$`\alpha `$/Fe\]$``$0.4. As discussed, e.g., by Norris et al. (1997a), it is difficult to explain these abundance ratios in terms of binary star evolution or ejecta from very massive Population III supernovae, which in particular produce few nitrogen but vast amounts of oxygen (Heger & Woosley, 2002), so that the origin of the surface abundances of UMP stars is still shrouded in mystery.

Concerning the nature of the UMPs, different ideas may be followed, the first one being that they are true Pop III stars — i.e. of initially zero metallicity — with the observed heavy elements (Fe, etc.) resulting from pollution of just the envelope by other stars such as supernovae. In two recent papers (Weiss et al. 2000, Paper I, and Schlattl et al. 2001, Paper II) we have thus discussed in detail the evolution of low-mass zero-metal stars, including the effect of atomic diffusion, mass loss and surface metal pollution. One important result was that the small entropy barrier between the H- and He-rich regions allows the He-flash driven convective zone to penetrate the overlying H-rich layers at the tip of the red giant branch. The consequent inward migration of protons into high-temperature regions leads to an H-shell flash, which further increases the extension of the central convective region. At the late stages of this phase the convective envelope deepens and merges with the inner convective zone. The initially metal-free surface is therefore enriched by a large amount of matter processed in He- and H-burning reactions. In particular, very high C and N abundances follow. We denote this process of surface C and N enrichment as HElium Flash induced Mixing (HEFM).

This scenario — already suggested by Fujimoto et al. (2000) — could potentially explain the anomalous abundance pattern at the surface of UMP stars, since via this mechanism high abundances of C and N are produced in the post He-flash phases. In Paper II we have shown that an initially metal-free star of 0.82$`M_{}`$, polluted by 0.0003$`M_{}`$ of $`Z`$=0.02 material, undergoes HEFM and reproduces in its post He-flash phase approximately luminosity and surface gravity of two of the best studied UMP stars, namely CS 22892-052 and CS 22957-027. The stellar mass was chosen such that the age on the RGB of 13.7 Gyr is compatible with the current estimates of the age of the universe, and with the age of CS 22892-052 estimated by Cowan et al. (1997) from nuclear chronology (15.6$`\pm `$4.6 Gyr). The amount and composition of the polluting material was adjusted such that the observed surface \[Fe/H\] abundance (\[Fe/H\]$``$$``$3) was obtained and, as the result of the HEFM, \[C/Fe\] and \[N/Fe\] ratios considerably larger than zero resulted. However, the values of \[C/Fe\] and \[N/Fe\] predicted by the models are about 2 orders of magnitude higher than the observed ones; i.e., we predict \[C/Fe\]$``$\[N/Fe\]$``$4, while observations of these stars yield \[C/Fe\]$``$\[N/Fe\]$``$1–2 (see Table 1). We note, however, that higher N overabundances are observed, e.g., in CS 22949-037 (Norris et al., 2001).

An alternative explanation for the C and N enhancement in UMPs is provided by HEFM during the thermally pulsating asymptotic giant-branch phase (TP-AGB) which occurs in metal-free stars with $`M`$$``$1$`M_{}`$ (Cassisi et al., 1996; Fujimoto et al., 2000). However, age, surface gravity, and effective temperature of these AGB stars are more difficult to reconcile with those of the observed objects, which are therefore unlikely to be TP-AGB stars themselves. The idea in this case is that the observed low-mass UMP objects form binary systems with more massive AGB stars, which have transferred part of their C- and N-enriched envelope to their companion (Fujimoto et al., 2000). A potential problem of this hypothesis is that HEFM during the TP-AGB of intermediate-mass stars may enhance the surface oxygen content similar to carbon. This feature can be found in models of Siess et al. (2002), which however do not agree in this respect with computations of Chieffi et al. (2001). An oxygen enhancement of \[O/Fe\]$``$1–2 might not necessarily be in contradiction to observations, as only very weak upper limits on the \[O/Fe\] ratio exist (Sneden et al., 1996; Norris et al., 1997a). But the main problem of this scenario is that not all UMPs are in binary systems, and the orbital periods of those who are, are inconsistent with an AGB mass transfer paradigm (Preston & Sneden, 2001). Moreover, the origin of Fe in these stars is as uncertain as in our favoured scenario. It is presently not clear whether intermediate-mass stars with \[Fe/H\]$``$$``$3 undergo HEFM during the TP-AGB or rather evolve like ordinary Pop II objects. In the latter case the surface C, N, O, and Fe abundances of UMP stars would then be the result of mass transfer from their companion and the accretion of further material, e.g., a nearby supernova.

In summary, HEFM during the major He flash is presently the most promising astrophysical explanation for the high \[C,N/Fe\] ratios observed in UMP stars despite their failure of reproducing the absolute value of the C and N enrichment. Therefore we are studying to what extend HEFM occurs also in low-mass stars with initial $`Z`$$`>`$0, and whether those objects are showing the observed \[C,N/Fe\] ratios of CS 22892-052 and CS 22957-027. This extension of our model assumption is warranted by the results about star formation in primordial environments, which agree on the fact that the formation of (very) massive objects is strongly favoured in the absence of metals (Abel et al., 2002; Bromm et al., 2002). If this indeed is the case, the UMPs would constitute the extreme end of Pop II, and had formed out of material already bearing the signature of individual supernova events (Shigeyama & Tsujimoto, 1998).

The final alternative for the origin of the carbon-rich UMPs is that they were already born with the abundances they now show. This scenario would receive support if single dwarfs of similar composition should definitely be identified. While one shifts the problem of identifying the source of the peculiar abundances to other stars, it is worthwhile, nevertheless, to investigate the evolution of such C-rich stars. For example, they are expected to reach standard RGB-tip luminosities and low gravities, which could be compared with those of observed objects. Stars of extremely low metallicity ($`Z=0`$ in particular), have a definitely shorter RGB, in contrast.

In addition to the three plausible models for the origin of the C-rich UMPs, we also investigate in this paper a physical problem associated with the HEFM. Due to the large uncertainties still affecting the treatment of convection in stellar models, we will assess the sensitivity of the HEFM process to different assumptions about the extension and mixing efficiency of convective regions in the stellar interiors.

The paper is structured as follows. Section 2 summarizes briefly the model computational techniques, input physics, and evolution of low-mass stars undergoing HEFM. Section 3 discusses the occurrence of HEFM and resulting surface abundances for different initial metallicities and various parameterizations of the stellar convective regions; Sect. 4 summarizes the results.

## 2 Modelling and evolution of stars undergoing HEFM during the He flash

All calculations (as the ones in Paper II) have been performed using the Garching stellar evolution code (Weiss & Schlattl 2000). The numerical details are presented in Paper II and will not be repeated here. The important feature of our code is the ability to follow in detail the evolution through the He flash considering simultaneously mixing and burning processes. Mixing is described by means of a time-dependent algorithm that treats convection as a fast diffusive process, where the diffusion constant is proportional to the convective velocity obtained from a convection theory. We have discussed in Paper II that the convective velocities predicted by both mixing-length and Canuto & Mazzitelli’s (1991) theory lead to very similar results. In general, arbitrary changes by up to 2 orders of magnitude in the mixing speed do not affect appreciably the HEFM process.

Atomic diffusion is treated according to Thoul et al. (1994), opacities are from Iglesias & Rogers (1996) and Alexander & Ferguson (1994), neutrino energy losses from Munakata et al. (1985), and the equation of state (EOS) is a Saha EOS plus a simplified description of the degenerate electron gas in the core regions following Kippenhahn & Weigert (1990). Our reaction network follows the abundances of H, $`{}_{}{}^{3}\mathrm{He}_{}^{}`$, $`{}_{}{}^{4}\mathrm{He}_{}^{}`$, $`{}_{}{}^{12}\mathrm{C}_{}^{}`$, $`{}_{}{}^{13}\mathrm{C}_{}^{}`$, $`{}_{}{}^{14}\mathrm{N}_{}^{}`$, $`{}_{}{}^{15}\mathrm{N}_{}^{}`$, $`{}_{}{}^{16}\mathrm{O}_{}^{}`$, $`{}_{}{}^{17}\mathrm{O}_{}^{}`$, without assuming a priori equilibrium compositions for any of these chemical elements. The reaction rates for the $`pp`$\- and CNO-burning have been taken from Adelberger et al. (1998), while the $`{}_{}{}^{12}\mathrm{C}_{}^{}(\alpha ,\gamma ){}_{}{}^{16}\mathrm{O}_{}^{}`$-rate agrees with Caughlan et al.’s (1985) value. The latter one is, in the relevant temperature range ($`0.1<T/(10^9\mathrm{K})<0.2`$), about a factor of 2 higher than the rates of Caughlan & Fowler (1988) and about 30–50% higher than the most recent value of Kunz et al. (2002).

Figure 1, as a reference, displays the evolutionary track of the model with 0.82$`M_{}`$, $`Z`$=0, and surface chemical pollution mentioned in the introduction (similar to model B1 of Paper II), which experiences the HEFM.

The 0.0003$`M_{}`$ of $`Z`$=0.02 polluting material are instantaneously accreted before the star reaches the zero-age main sequence; this approach maximizes the effect of accretion and diffusion effects. After this episode the star evolves along the main sequence producing energy through the $`pp`$-chain. Diffusion is never able to bring the surface metals deep enough to contribute to the burning.

At the sub-giant branch the core increases its temperature and density, and gradually 3-$`\alpha `$ reactions are setting in. In stars more massive than $`0.9M_{}`$ the central H content is not fully exhausted when sufficient C is produced ($`X({}_{}{}^{12}\mathrm{C}_{}^{})10^{10}`$) to gain considerable nuclear energy from the CNO-cycle. As a consequence a thermal runaway through the CNO-cycle ensues, causing a characteristic blue loop in the H-R diagram. The consequent expansion of the inner regions reduces the efficiency of the 3-$`\alpha `$ process and CNO-cycle, terminating the runaway which lasts only about $`10^7`$ years. In less massive models, as the ones considered in this work, the central H is already exhausted before sufficient C can be produced, and thus no blue-loop occurs on the sub-giant branch.

Regardless of the occurrence of this CNO flash, all low-mass models settle on the RGB and evolve towards higher luminosities. While on the sub-giant branch their main energy source remains the $`pp`$-chain, the contribution of the CNO-cycle is increasing gradually with luminosity as more carbon is created at the inner tail of the H-burning shell by 3-$`\alpha `$ reactions. In a 0.8$`M_{}`$ star about 50% of the H-burning energy star is produced by the CNO-cycle at the tip of the RGB.

At this stage the mass difference between the location of the maximum energy generation in the core and the He-core boundary is $``$0.29$`M_{}`$. Soon after the onset of the He flash a convective shell develops, which reaches the H-rich matter in the envelope about 1 month after the flash (left arrow in Fig. 2); as soon as H is carried into the interior hotter convective region, it starts to burn at a very high rate. Due to this extra energy input the upper boundary of the convective shell immediately moves closer to the surface (thus ingesting even more protons), while the He-burning rate is significantly reduced. The single inner convective region splits into two zones, one for the He- and one for the H-burning. The rapidly weakening He-burning shell can no longer support the underlying layers against contraction, and the released gravothermal energy forms a further convective region which disappears again after a few 1,000 years555This convective region appeared also in the calculations presented in Paper II, but was not shown in the corresponding figure due do insufficient plot data..

About 500 yr after the H ingestion the convective envelope deepens and merges with the convective region above the H-burning region (horizontal arrow in Fig. 2). A huge amount of matter processed through H and He burning is brought to the surface enriched in He, C (produced during the He burning) and N (produced by the CNO cycle at the expenses of C). The resulting surface \[C/Fe\] and \[N/Fe\] ratios are, respectively, 4.1 and 4.3, while the surface oxygen abundance is practically unaltered with respect to the initial value (Table 2).

As a result of this dredge up of heavy elements (which happens on timescales of a few weeks) the envelope opacity is increased significantly; this produces an abrupt discontinuity in the effective temperature visible in Fig. 1. At this stage, with a practically ceased He burning, the star behaves as a newborn RGB star, climbing its own RGB. The H shell produces the energy needed to support the star. Another He flash ensues at the end of this second RGB phase, this time with a smaller and less degenerate core. The flash is weaker and no HEFM happens. The subsequent evolution does not show any peculiar feature.

As a general rule, all the changes in the physical inputs and/or initial conditions that can contribute to an increase of the electron degeneracy and thus cause the location of the He-flash ignition to move closer to the border of the He core favour the occurrence of HEFM. Changes in the relative location as small as $``$0.02$`M_{}`$ can make a difference between models undergoing HEFM and models which do not experience it. In Paper II we have shown how an increase of the initial He abundance, or an increase of the initial mass, or the inclusion of heavy element pollution at the stellar surface are disfavouring the onset of HEFM. On the other hand, the inclusion of element diffusion favours this process, while the inclusion of mass loss from the stellar surface does not affect at all the onset of the HEFM.

## 3 The influence of initial metallicity and convection

In Paper II we have investigated the dependence of the HEFM process on several input-physics parameters and various assumptions adopted in computing stellar models; however, we have considered only stars with initial zero metallicity, plus some eventual metal pollution at the surface. It is worthwhile to investigate the dependence of the HEFM on the value of the initial stellar metallicity, given the unknown composition of these stars at the time of formation, as outlined in the introduction.

Because of the uncertainty in the treatment of convection in stellar structures, it is also important to verify how much the HEFM process is sensitive to different assumptions about the extension and mixing efficiency of convective regions in the stellar interior. In particular, we considered various amounts of overshooting from the canonical formal convective boundaries fixed by the Schwarzschild criterion.

In this section we will discuss these effects in theoretical models representative of CS 22892-052 and CS 22957-027.

### 3.1 Models with $`Z`$$`>`$0

Models with initial $`Y`$=0.23 and increasing $`Z`$ have been computed in order to derive the value of the maximum initial metallicity which produces HEFM in a star with 0.82$`M_{}`$. The heavy element mixture has been considered to be $`\alpha `$-enhanced, with \[$`\alpha `$/Fe\]=0.4; carbon and nitrogen were scaled solar. We obtain HEFM only for $`Z<10^7`$, corresponding to \[Fe/H\]$`<`$$``$5.6. Since all observed \[Fe/H\]-abundances in UMPs so far are higher than $`4`$, this already excludes the possibility that the C/N-anomalies were produced in stars with initial Fe-abundances as observed today. Due to the mixing of H-rich material into the H-depleted interior during the HEFM, the surface H abundance drops by about a factor of two, and the upper limit for the surface Fe abundance after this phase is \[Fe/H\]$`<`$$``$5.3, about 2 orders of magnitude smaller than the observed metallicities of CS 22892-052 and CS 22957-027. Whenever HEFM happens, the \[C/Fe\] and \[N/Fe\] ratios in the post He-flash phase are always of the order of 4, still about 2 orders of magnitude higher than the observed values.

This upper value for the metallicity is slightly increased if we take into account the revised plasma-neutrino rates by Haft et al. (1994). In this case the limiting metallicity is $`Z=10^6`$, corresponding to \[Fe/H\]=$``$4.3 after the HEFM, i.e., 1 dex higher than in case of the old neutrino losses. The reason is that the degeneracy of the core increases and the He flash starts more off-centre. We also wish to notice here that, in spite of their different input physics and cruder treatment of the coupling between convection and burning, Fujimoto et al. (2000) found a similar upper limit for the metallicity of a 0.80$`M_{}`$ star undergoing HEFM, and similar values of the surface \[C,N/Fe\] ratios.

We have also tested the effect of implementing a new EOS by Irwin (2000), which covers the entire stellar structure in all relevant evolutionary phases, and closely reproduces the OPAL EOS (Rogers et al., 1996) in the common validity range (see, e.g., Schlattl, 2002). The inclusion of this new EOS in the models does not modify the value of the largest critical heavy elements abundance for which the HEFM occurs. However, with this equation of state the stellar ages are slightly reduced (by about 0.5 Gyr) with respect to our reference EOS, and therefore one has to choose a slightly lower mass to reproduce the observations, which might favour the occurrence of HEFM.

Therefore, we have computed a, what we call, ‘best model’ with an initial \[Fe/H\] equal to the observed one, and a mass which yields an age at the RGB tip well matching the age estimated for CS 22892-052 and CS 22957-027; we considered $`M`$=0.78$`M_{}`$, $`Z`$=$`2.0\times 10^5`$ (which corresponds, including the $`\alpha `$ enhancement by 0.4 dex, to \[Fe/H\]=$``$3.3), neutrino emission following Haft et al. (1994), Irwin’s (2000) EOS, plus atomic diffusion of H, He, and metals (which further reduces the ages by about 0.5 Gyr). The corresponding evolutionary track is shown in Fig. 1, too (solid line). In spite of the reduced stellar mass the effect of the initial metallicity dominates, and *no* HEFM occurs. The location of the maximum energy generation at the flash is 0.30$`M_{}`$ away from the boundary of the He core (solid line in Fig. 3), a distance about 0.01$`M_{}`$ larger than in the cases when HEFM occurs. We computed also the evolution of stars of 0.75$`M_{}`$ and 0.82$`M_{}`$ with the same initial composition as the 0.78$`M_{}`$ star to span the age range determined for CS 22892-052 from nuclear chronology, not obtaining HEFM in any case.

### 3.2 Models including overshooting

Models including overshooting have been computed in order to check its influence on both the onset of HEFM and the amount of C and N dredged up to the surface. Since overshooting increases the extension of convective regions, the occurrence of HEFM is in principle favoured. Our overshooting description follows Blöcker et al. (1998) and it is modelled as an exponential diffusive process. The diffusion constant of the overshoot region decays exponentially outside the Schwarzschild convective boundary, starting from the value assumed by the convective diffusion coefficient at the convective boundary; the decay length is $`F`$ times the pressure scale height at the convective boundary. Hydrodynamical simulations of shallow convective envelopes in A stars by Freytag et al. (1996) predict $`F`$=0.25$`\pm `$0.05, while for the overshooting from the convective envelopes of AGB stars Herwig et al. (1997) took $`F=0.02`$.

We have computed a series of models (adopting the same input physics as in Paper II) with different values for $`F`$, considering $`M`$=0.82$`M_{}`$, $`Y`$=0.23 and $`Z`$=$`2.0\times 10^5`$ which corresponds to \[Fe/H\]=$``$3.3, approximately the \[Fe/H\] value determined for the two UMP stars under scrutiny. In these models the simpler EOS (Saha-EOS completed with degenerate electron gas in the deep interior) and the neutrino losses of Munakata et al. (1985) have been used.

We have first considered overshooting from all the convective boundaries, starting from the main-sequence phase. As a general result, we found that the inclusion of overshooting during the main sequence increases the lifetime of the small convective core which develops at the beginning of the main sequence, when the $`{}_{}{}^{3}\mathrm{He}_{}^{}`$ abundance gradually reaches nuclear equilibrium in the centre; this, in turn, causes an increase of the main-sequence lifetime and a decreased electron degeneracy in the He core along the RGB, which disfavours the occurrence of HEFM. Moreover, the increased lifetime would demand a larger stellar mass in order to obtain not too high stellar ages, and that would further reduce the probability for HEFM.

Models with $`M`$=0.82$`M_{}`$ and F$`=`$0.03 do not show HEFM, and are continuing their evolution through the horizontal branch to the TP-AGB (dash-dotted line in Fig. 1). During this phase the inter-shell convective zone is able to penetrate the H-rich envelope, similar to the HEFM process (Cassisi et al., 1996; Fujimoto et al., 2000). The surface C and N abundances rise again by about 4 dex, while O is increased by about 3 dex. Hence, an O enrichment about one order of magnitude smaller than for C is expected in these stars, which implies that \[O/Fe\] should be at most 1 for the stars in Table 1. Since no oxygen abundance could be determined in these stars, but only an upper limit of \[O/Fe\]$``$0.6 (Sneden et al., 1996, for CS 22982-052), this scenario could provide a possible explanation for C-rich UMPs, too. However, surface gravity and temperature of these AGB stars are too low to be in agreement with the stars under scrutiny (Fig. 1).

Higher values of $`F`$ do not lead to HEFM. In fact, for $`F=0.1`$ there is no He flash, but quiescent He-burning ignition and of course no HEFM.

In case overshooting is included only from the RGB phase onwards, the stellar lifetime is not affected, and one needs at least $`F=0.17`$ to get HEFM. This value is about 30% smaller as the one obtained by Freytag et al. (1996) for A-star envelopes, but appears to be quite high compared to that needed in AGB stars (Herwig et al., 1997); moreover, switching on overshooting only during the RGB phase appears to be entirely ad hoc; in addition, the resulting \[C/Fe\] and \[N/Fe\] ratios after the HEFM are again of the order of 4, i.e., too high with respect to the observations.

### 3.3 Models with reduced mixing efficiency

From the computations performed until now, it appears that one of the main problem to the HEFM scenario is the too high C and N enrichment of the surface. One solution to obtain lower surface \[C,N/Fe\] ratios would be to have a smaller overlap between the convective region developing at the flash, and the overlying H-rich region. Smaller overlap implies less protons and therefore less processed C and N. This smaller overlap could be achieved by strongly reducing the diffusion coefficient associated to the convective mixing. Using standard mixing-length velocities, the crossing time of the convective shell is of the order of a few hours. A reduction of the “convective” diffusion coefficient is performed by simply reducing the multiplicative factor that links the coefficient itself to the convective velocity.

The upper boundary of the He-flash driven convective shell is at first order dictated by the energy of the flash and the location of the ignition point. However, after some protons are engulfed the increased energy production (due to the additional H burning) pushes the upper boundary further up, thus increasing the number of proton ingested. If the mixing efficiency is reduced, protons are mixed less deep and less energy is produced, with a smaller increase of the convective shell extension. This would cause a smaller final number of protons ingested and a lower N production.

In order to explore this scenario in more detail, we recomputed the polluted $`Z`$=0 model of Sect. 2, which experiences HEFM, decreasing the “convective” diffusion coefficient by a factor of $`2\times 10^4`$. The extension of the convective regions are displayed in Fig. 4. The upper boundary of the H-flash driven convective shell extends less into H-rich layers, and the location of the H-burning region is less deep; the thickness of the H convective shell before the dredge up is also thinner, implying a reduction of the amount of C dredged up with respect to the case with standard convective mixing. The final surface abundance ratios \[C/Fe\] and \[N/Fe\] are in this case 2.3 and 3.2 (Table 2), respectively, closer to the observed values.

The surface $`{}_{}{}^{12}\mathrm{C}_{}^{}`$/$`{}_{}{}^{13}\mathrm{C}_{}^{}`$ ratio is also affected, being now increased to 5.6 compared to 4.5 which we obtain by using the standard mixing efficiency. This increased value is in slightly better agreement with the observed values of about 10. In spite of the much reduced mixing efficiency the dredge up still happens fast, on timescales of the order of $`10^3`$ years.

## 4 Discussion

The results shown in the previous section highlight the fact that it is very difficult to reproduce the surface abundances of CS 22892-052 and CS 22957-027 from single-star evolution, at least in the case the observed surface abundances pattern is not primordial. To investigate this alternative, we have computed additionally the evolution of an 0.82$`M_{}`$ star, with the same physics of our ‘best model’ and $`Y`$=0.23, but with $`Z`$=$`6\times 10^4`$ and \[C/Fe\]=\[N/Fe\]=2 (implying \[Fe/H\]=$``$3.3), in agreement with the observed abundances. The evolutionary track is shown in Fig. 1, while the He-core boundary before and at He ignition is plotted in Fig. 3. As expected, this model does not experience HEFM (the He flash starts very deep in the core), so that the surface \[C/Fe\] and \[N/Fe\] ratios are not modified with respect to the initial values. Note that the effective temperature along the RGB is higher than for the case with HEFM, because of the much lower C-abundance. The location of the track fits reasonably well the observed gravities and effective temperatures, and thus, this model would be able to reproduce the observed properties of CS 22957-027. However, as discussed in Sect. 1, if this were a realistic model for the observed object, the source for the high carbon and nitrogen abundances still remains unknown, as it is not clear how to pollute the interstellar medium with matter showing these anomalous abundance ratios.

Thus, assuming that the observed surface composition is not primordial, the HEFM is presently the only way to get a surface C and N enhancement without invoking some kind of mass transfer from binary companions. From our tests in Sect. 2 there appear to be two main problems with the HEFM scenario:

Firstly, too much of C and N is transported to the surface after the HEFM, and secondly, no HEFM appears to occur for an initial metallicity that can match the observed \[Fe/H\] of the stars under scrutiny, unless surface heavy element pollution has been efficient, as discussed in Paper II.

Once again, we want to point out that, in spite of their different input physics and cruder treatment of the coupling between convection and burning, Fujimoto et al. (2000) found similar upper limits for the metallicity of stars undergoing HEFM, and similar values of the \[C,N/Fe\] ratio for stellar mass and metallicity values similar to our ‘best model’.

Concerning the appearance of HEFM in stars with \[Fe/H\] in agreement with observations, one possible solution would be, as discussed in Sect.3.3, a completely ad hoc choice of overshooting efficiency. Another possibility is related to the H depletion during the dredge up of the HEFM products. During this phase the convective envelope reaches deeper regions involved in the HEFM, which are devoid of H. According to our models the surface \[Fe/H\] increases by $``$0.3 dex because of the dilution of H. If the drop in hydrogen abundance is about five times higher, a surface \[Fe/H\] of about $``$3.3 could be obtained with an initial \[Fe/H\]=$``$4.3 corresponding to the upper limit of getting HEFM.

In order to achieve a larger H depletion, a smaller envelope thickness (in mass) is needed, i.e., a smaller H reservoir. We estimate that right before the flash ignition an envelope of only about 0.03$`M_{}`$ is needed in order to reproduce the observed post-HEFM \[Fe/H\] surface values with our ’best model’; with this envelope mass the star would still be able to ignite He close to its RGB location (see, e.g., Castellani & Castellani 1993), i.e., the effective temperature and gravity after the HEFM would still be in accordance with CS 22892-052 and CS 22957-027. Taking into account the age of the star of about 14 Gyr, which fixes its initial mass to be about 0.8$`M_{}`$, it would imply that about 0.27$`M_{}`$ are lost during the evolution, a value which is not much higher than the average mass loss experienced by globular cluster stars (about 0.20$`M_{}`$). The lack of metals in the outer layers of UMPs, however, should cause a strong decrease of mass-loss efficiency in comparison with more metal-rich stars, if radiative driven winds are the main source of mass loss. Nevertheless, due to our poor knowledge about mass-loss processes in RGB stars, this possibility cannot be completely ruled out.

###### Acknowledgements.

We warmly thank A. Irwin for allowing us to use his new equation of state as well as for help in managing his routines. We also thank him for discussions and suggestions about the use of his EOS in the evolutionary code. We are grateful to T. Beers for his continuing interest and support of this work. One of us (S.C.) acknowledges the financial support by MURST/Cofin2000 under the project ”Stellar Observables of Cosmological Relevance”. H.S. has been supported by a Marie Curie Fellowship of the European Community programme “Human Potential” under contract number HPMF-CT-2000-00951. A.W. thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work.


format: text

text:
# The magnetic response of disordered metallic rings: the large contribution of the far levels

## Abstract

We calculate the orbital linear magnetic response of disordered metallic rings to an Aharonov-Bohm flux using the BCS model for attractive electron-electron interaction. The contribution of all levels including those up to a high energy cutoff results in a much larger value than previously obtained using the local interaction model. The possible relevance of our results to the resolution of the discrepancy between the experimental and theoretical values for the ensemble-averaged persistent currents in these systems is discussed.

One of the remarkable phenomena of mesoscopic physics is the existence of equilibrium persistent currents in small normal metal rings, in the regime where the elastic mean free path $`l`$ is much smaller than the ring’s circumference $`L`$. This was both predicted theoretically  and observed experimentally . While the phenomenon and the periodicity in the flux are well understood by now, the magnitude of the ensemble-averaged persistent current found experimentally is much larger than that obtained theoretically using the model of noninteracting electrons. An attempt to account for this discrepancy by the inclusion of electron-electron (“e-e”) interactions indeed increased the theoretical value, but still came short by a factor of about $`5`$.

Along with the large value of the ensemble-averaged persistent current observed, its sign poses significant questions as well. While in the first experiment the sign of the ensemble-averaged persistent current was tentatively identified as diamagnetic, it was not determined definitively (see Ref. ). Later experiments reported a predominantly diamagnetic ensemble-averaged persistent current. Theoretically, the paramagnetic/diamagnetic sign follows for repulsive/attractive effective e-e interactions (the interplay between the renormalized repulsive e-e interaction and the phonon mediated attractive interaction determines the sign of the effective interaction). Thus, the diamagnetic sign of the persistent currents observed in Gold, Copper, and Silver suggests that these materials are very weak superconductors.

The model used to treat both signs of the interaction was that of local e-e interaction. The discrepancy between the theoretical values and the experimental results exists for repulsive, as well as attractive interactions . For the repulsive case, the higher order renormalization which reduces the first order result is responsible for this discrepancy. For the attractive interaction case, the higher order renormalization does not reduce, and even somewhat increases the first order result. However, the restriction on the bare interaction constant, due to the experimental fact that Gold, Copper and Silver are not superconducting in observed temperatures, limits the value of the theoretical result in this case. For the magnetic response (the derivative of the persistent current at zero flux), there exist a further logarithm in the numerator, irrespective of the above renormalization of the effective interaction. Here we show that the use of the BCS Hamiltonian to model attractive e-e interaction gives rise to a much larger cutoff of this logarithm, and hence to a significantly larger value for the magnetic response to flux. The enhancement of the magnetic response is a result of the contribution of high energy levels (to be called “far levels”) originating from pairing correlations which persist up to the energy of the Debye frequency $`\omega _\mathrm{D}`$. Thus, correlations on an energy scale much larger than the Thouless energy affect the magnetic response.

The contribution of the far levels was investigated in connection to superconductivity in small grains . There, the correlations of levels much further than the superconducting gap from the Fermi energy $`E_\mathrm{F}`$, up to $`\omega _\mathrm{D}`$ prove to be significant. Using an exact solution of the reduced BCS Hamiltonian , it was shown that this contribution results in a much larger condensation energy than that given by the BCS value, in a wide parameter regime in which superconducting correlations are well developed. The far levels also affect single particle properties of a superconducting grain, such as the Matveev-Larkin parameter.

The difference between the local interaction model and the BCS Hamiltonian lies in the $`q`$ dependence of the interaction in each, where $`q`$ is the sum of the incoming (outgoing) momenta of the pair scattering. While the local interaction model assumes that the bare interaction does not depend on any momenta, and specifically on $`q`$, the BCS Hamiltonian assumes a sharp cutoff at small total momentum, and considers the interaction to be a delta function at $`q=0`$. While for physical properties related to the superconducting phase of the system, near and below $`T_\mathrm{c}`$ (such as $`T_\mathrm{c}`$ itself) the two models give similar results, we show that for the magnetic response in the perturbative regime (equivalent to $`TT_\mathrm{c}`$) the results in the two models differ substantially. We discuss the relation between the two models, and the effect of relaxing the BCS assumption by taking the interaction to have finite width in $`q`$.

Recently it was suggested that the largeness of the observed current is due to AC noise  or interactions of the electrons with two level systems , and a relation to dephasing was suggested (one should also keep in mind that the latter interactions may lead to an additional attractive interaction). We, however, consider $`T=0`$, and do not include dephasing (we consider the dephasing length $`L_\varphi `$ to be larger than all relevant lengths). A finite $`L_\varphi `$ will result in a suppression of the persistent current, exponential in $`L/L_\varphi `$.

We consider a quasi one dimensional disordered ring penetrated by a constant Aharonov-Bohm (AB) flux $`\mathrm{\Phi }`$ in its middle. We calculate the magnetic response, denoted by $`F_2`$, to first order in the e-e interaction, for attractive interaction using the reduced BCS Hamiltonian.

The result we obtain with logarithmic accuracy is

$$F_2\frac{dI}{d\mathrm{\Phi }}|_{\mathrm{\Phi }=0}=\frac{8\pi \lambda E_{\mathrm{Th}}}{\mathrm{\Phi }_0^2}\mathrm{ln}\frac{E_{\mathrm{co}}}{d},$$
(1)
where $`\lambda `$ is the dimensionless pairing parameter \[see Eq. (7)\], $`d`$ is the level spacing, $`E_{\mathrm{Th}}=l^2/(s\tau L^2)`$ is the Thouless energy, which is the inverse of the time to circulate the ring ($`s=1,2,3`$ is the effective dimension of the ring for diffusive motion), $`\mathrm{\Phi }_0=hc/(2e)`$ and $`\mathrm{}`$ denotes ensemble averaging. The upper logarithmic cutoff $`E_{\mathrm{co}}\mathrm{min}[\omega _\mathrm{D},1/\tau ]`$ represents the fact that energies (measured from the $`E_\mathrm{F}`$) further than $`E_{\mathrm{Th}}`$, up to $`E_{\mathrm{co}}`$ (far levels) contribute to the magnetic response ($`\tau `$ is the elastic mean free time). Their contribution enhances $`F_2`$ by about an order of magnitude due to the much larger logarithmic cutoff as compared to $`E_{\mathrm{Th}}`$ of the known result  \[for $`1/(\tau E_{\mathrm{Th}})10^4`$ the ratio of the logarithms is roughly $`10`$\].

In order to obtain Eq. (1) we consider the general Hamiltonian

$$H=\underset{\alpha }{}𝑑r\psi _\alpha ^{}(r)\left[\frac{1}{2m}(\stackrel{}{P}\frac{e}{c}\stackrel{}{A})^2+U(r)\right]\psi _\alpha (r)+H_{\mathrm{ee}},$$
(2)
where $`U(r)`$ is the external potential which includes the disorder, and $`H_{\mathrm{ee}}`$ represents the e-e interaction. In the London gauge $`\stackrel{}{A}=\mathrm{\Phi }/(2\pi \rho )\widehat{\varphi }`$ where $`\rho `$ is the distance from the origin and $`\widehat{\varphi }`$ is in the clockwise direction of the ring. The free energy of the system and the persistent current are flux dependent, and related by $`I(\mathrm{\Phi })=dF/d\mathrm{\Phi }`$. By time reversal symmetry $`I(0)=0`$, and for small flux

$$F(\mathrm{\Phi })=F_0\frac{1}{2}F_2\mathrm{\Phi }^2+\mathrm{};\frac{dI}{d\mathrm{\Phi }}|_{\mathrm{\Phi }=0}=F_2.$$
(3)
We now turn to the calculation of $`F_2`$ to first order in the interaction. Perturbative analysis of the reduced BCS attractive interaction \[Eq. (7) below\] is valid, at $`T=0`$, for $`\lambda <1/\mathrm{ln}(\omega _\mathrm{D}/d)`$ and at finite temperature ($`T>d`$) for $`TT_\mathrm{c}`$. We take as the unperturbed Hamiltonian

$$H_0=\underset{\alpha }{}𝑑r\psi _\alpha ^{}(r)\left[\frac{P^2}{2m}+U(r)\right]\psi _\alpha (r)$$
(4)
and the magnetic field together with the e-e interaction as perturbation

$`H_I=H_{\mathrm{ee}}+`$ (5)
$`{\displaystyle \underset{\alpha }{}}{\displaystyle 𝑑r\psi _\alpha ^{}(r)}`$ $`\left[{\displaystyle \frac{e}{2mc}}(\stackrel{}{P}\stackrel{}{A}+\stackrel{}{A}\stackrel{}{P})+{\displaystyle \frac{e^2}{2mc^2}}A^2\right]\psi _\alpha (r)`$ $`.`$ (6)
We calculate perturbatively the energy to third order in $`H_\mathrm{I}`$, and consider only the terms which are second order in the flux and first order in the interaction.

We denote by $`|i`$ the eigenstates of the noninteracting electrons in the disordered ring without magnetic field. In this basis, within the reduced BCS model, we obtain

$$H_\mathrm{I}^{\mathrm{BCS}}=\underset{ij\alpha }{}\frac{e\mathrm{\Phi }}{mcL}P_{ij}c_{i\alpha }^{}c_{j\alpha }+\lambda d\underset{ij}{\overset{}{}}c_i^{}c_i^{}c_jc_j.$$
(7)
Here $`c_i`$ destroys an electron in the state $`|i`$ with wavefunction $`\chi _i(r)`$ and $`P_{ij}=i|P_{}|j`$ is the matrix element of the momentum parallel to the ring’s direction. The second sum is restricted to levels within $`\omega _\mathrm{D}`$ of $`E_\mathrm{F}`$. Note that contrary to the usual convention, $`\lambda `$ is negative for attractive interaction. We assume that the width of the ring is much smaller than its radius . The $`A^2`$ term is then interaction independent. We choose the $`\chi _i`$’s to be real, and then $`P_{ij}`$ is pure imaginary and $`P_{ii}=0`$.

Using third order perturbation theory  we find that for $`T=0`$ to first order in $`\lambda `$

$$F_2=8\left(\frac{e}{mcL}\right)^2\lambda d\underset{\begin{array}{c}\\ 0<i<\omega _\mathrm{D}\\ \omega _\mathrm{D}<j<0\end{array}}{}\frac{|P_{ij}|^2}{\omega _{ij}^2},$$
(8)
where $`0<i<\omega _\mathrm{D}`$ represents states whose energy is between $`E_\mathrm{F}`$ and $`E_\mathrm{F}+\omega _\mathrm{D}`$, and $`\omega _{ij}=ϵ_iϵ_j`$. Performing disorder averaging over $`|P_{ij}|^2`$ we obtain

$$|P_{ij}|^2=\frac{p_\mathrm{F}^2d\tau }{\pi (1+\omega _{ij}^2\tau ^2)s},$$
(9)
which is roughly constant for $`\omega _{ij}<1/\tau `$ and zero for $`\omega _{ij}>1/\tau `$. In Eq. (8), for energies smaller than $`\omega _\mathrm{D}`$ and $`1/\tau `$ there is a double sum and a second power of energy in the denominator. The result is therefore a logarithm with an upper cutoff which is given by $`\mathrm{min}[\omega _\mathrm{D},1/\tau ]`$. We evaluate the sum in Eq. (8) taking for simplicity the noninteracting spectrum to be equally spaced, and obtain Eq. (1).

This equation differs from the result of Ambegaokar and Eckern (AE) by the upper cutoff of the logarithm, being $`\mathrm{min}[\omega _\mathrm{D},1/\tau ]`$ in comparison to $`E_{\mathrm{Th}}`$ in Ref. . We now explain the origin of this difference in some detail. AE consider the same Hamiltonian as in our Eq. (2), only with local e-e interaction $`V(r_1r_2)=\stackrel{~}{\lambda }N(0)^1\delta (r_1r_2)`$. In terms of the noninteracting eigenstates $`H_{\mathrm{ee}}`$ is then given by

$$H_{\mathrm{ee}}^{\mathrm{local}}=\frac{1}{2}\underset{ijkl,\sigma }{}V_{ijkl}c_{i\sigma }^{}c_{j\overline{\sigma }}^{}c_{k\overline{\sigma }}c_{l\sigma },$$
(10)
where $`V_{ijkl}=\stackrel{~}{\lambda }N(0)^1\chi _i(r)\chi _j(r)\chi _k(r)\chi _l(r)𝑑r`$ and $`\sigma \overline{\sigma }`$ due to the interaction being local. The dependence of the energy on flux, to first order in the interaction, can be written as

$$\mathrm{\Delta }F=\underset{i,j<0}{}V_{ijji}(\mathrm{\Phi })=\underset{ij<0}{}V_{ijji}(\mathrm{\Phi })+\underset{i<0}{}V_{iiii}(\mathrm{\Phi }).$$
(11)
We now show that the second, diagonal part contribution to the magnetic response corresponds to the BCS result as given in Eq. (1). $`\mathrm{\Delta }F`$ can be expanded to second order in the flux to obtain $`F_2`$ for the local interaction model by following the same procedure done above for the reduced BCS Hamiltonian, with $`H_{\mathrm{ee}}`$ given by Eq. (10). The result is

$$F_2=4\left(\frac{e}{mcL}\right)^2\underset{j<0,i>0}{}\frac{|P_{ij}|^2}{\omega _{ij}^2}\underset{l<0}{}(V_{jllj}V_{illi}).$$
(12)
The diagonal contribution to this result is given by the $`l=i`$ terms, and since to leading order $`V_{iiii}V_{ijji}=2\stackrel{~}{\lambda }d`$ we obtain for the diagonal contribution in the local interaction model the same expression as Eq. (8), without the explicit cutoff at $`\omega _\mathrm{D}`$. Due to Eq. (9), the diagonal contribution gives Eq. (1) with $`E_{\mathrm{co}}=1/\tau `$. Therefore, the BCS approximation is equivalent to assuming that the flux dependence of the offdiagonal matrix elements is small, and can be neglected. In the case of the local interaction model this is not the case, and the contribution of the offdiagonal elements to the magnetic response is significant, and opposite in sign to that of the diagonal element. This results in a partial cancellation, and therefore a reduction of the high logarithmic cutoff in Eq. (1).

Another way to understand the relation between the two models is to consider the $`q`$ dependence of the interaction in both (where $`q`$ is the sum of the incoming momenta). AE derive their result as a sum over all $`q`$’s (see Eqs. \[(12)-(14)\] of Ref. ), and the result within the BCS Hamiltonian \[Eq. (1)\] corresponds to their $`q=0`$ term. For the local interaction model, the contribution to the magnetic response of the high energies, above $`E_{\mathrm{Th}}`$ exactly cancels between equal magnitude and opposite signs of the $`q=0`$ term and the sum of all $`q0`$ terms. *This makes the $`q`$ independent assumption for the bare interaction crucial*. The existence of excess interaction at small total momentum $`q`$ would thus significantly affect the result. The BCS interaction assumes just that, the existence of excess interaction at $`q=0`$.

A physical justification for taking a $`q`$ dependent interaction can be obtained, for the attractive phonon mediated interaction, from the usual restriction that all (incoming and outgoing) scattered states are within $`\omega _\mathrm{D}`$ of $`E_\mathrm{F}`$ (see e.g. Ref. ). This restriction implies a significant $`q`$ dependence on the scale of $`q_c=\omega _\mathrm{D}/v_\mathrm{F}`$. One can therefore take a $`q`$ dependent attractive interaction in the form of a step function with width $`q_c`$. For $`q_c<2\pi /L`$ one finds that only the $`q=0`$ term survives, and the BCS result for the magnetic response is recovered.

For metals such as Gold, Copper and Silver, this condition requires the circumference of the ring to be smaller than $`0.20.4\mu m`$, which is an order of magnitude smaller than the relevant experimental lengths. The effects of relaxing this condition, i.e. having $`q_c2\pi /L`$, as well as the effect of the $`q`$ dependence on the persistent current itself will be considered elsewhere.

The central result of this Letter is the large interaction correction to the derivative of the persistent current at zero flux, within the BCS model, as is given in Eq. (1). Physically, this is due to the large contribution of the far levels, up to $`\mathrm{min}[\omega _\mathrm{D},1/\tau ]`$ from the Fermi level.

Our calculations were done for $`T=0`$. At finite temperature $`T<\mathrm{min}[\omega _\mathrm{D},1/\tau ]`$ (but neglecting dephasing, assuming $`L_\varphi L`$) the magnetic response would be given by Eq. (1) with $`d`$ replaced by $`T`$. Therefore, on top of the large magnetic response, we predict a weak, logarithmic, temperature dependence of the magnetic response up to $`T\mathrm{min}[\omega _\mathrm{D},1/\tau ]`$. The weak temperature dependence is due to the addition of the many small contributions, all having the same sign, of the levels up to the large energy cutoff (see also Refs.).

The local interaction model and the reduced BCS Hamiltonian are two models used to describe the effective e-e interaction resulting from the retarded electron-phonon interaction. Though different, these two models give similar results for many properties of the superconductor near and below $`T_\mathrm{c}`$. The robustness of the physics near and below the transition temperature makes the differences between the two models irrelevant. However, the magnetic response in the perturbative regime is a more subtle property, that distinguishes between the two models.

Our result for the reduced BCS Hamiltonian also follows by expanding Richardson’s exact solution. This has the analytical merit of being a first order expansion of an exact solution, as well as an additional viewpoint on the contribution of the far levels, and will be given elsewhere.

We would like to thank O. Agam and I. Aleiner for discussing with us the special role of the $`q=0`$ term in the calculation of the persistent current within the local interaction model. We benefited from valuable discussions with N. Argaman, V. Kravtsov, D. Orgad, A. Punnoose, Y. Aharonov, B. Altshuler, R. Berkovits, H. Bouchiat, P. Brouwer, J. von Delft, A. Finkel’stein, Y. Gefen, D. E. Khmelnitskii, M. Khodas, F. von Oppen, Z. Ovadyahu, R. Smith, U. Smilansky, R. A. Webb, and H. A. Weidenmueller. M.S. is thankful for the support by the Lady Davis fund. This work was supported by a Center of Excellence of the Israel Science Foundation (ISF), Jerusalem and by the German Federal Ministry of Education and Research (BMBF) within the Framework of the German-Israeli Project Cooperation (DIP), by the German-Israeli Foundation(GIF), and by ISF grant number 160-01/1.


text:
# 𝑑+𝑖𝑠 vs 𝑑+𝑖𝑑' Time Reversal Symmetry Breaking States in Finite Size Systems

## Abstract

We report self-consistent quasiclassical calculations of spontaneous currents and magnetic moments in finite size unconventional superconducting systems, namely: (i) in isolated $`d`$-wave superconductor islands where, in addition to the dominant order parameter (with a $`d_{x^2y^2}`$ symmetry), a subdominant order parameter of $`s`$ or $`d_{xy}`$ symmetry is added; (ii) in grain boundary junctions between two arbitrarily oriented $`d`$-wave superconductors, and between a $`d`$-wave and an $`s`$-wave superconductor. We show that the profile of the spontaneous current density and the magnetic field distribution depend on the time-reversal symmetry breaking properties of the system. For the $`d_{x^2y^2}+id_{xy}`$ state, vortices appear near the edges of the finite size systems. We associate these vortices with the chiral nature of the mixed order parameter. The method developed here is quite general, and can be used for predicting properties of any finite size superconducting system.

The possibility of the existence of mixed order parameter symmetries in unconventional superconductors has been suggested and recently investigated both theoretically  and experimentally . A significant feature of the mixed symmetry states is that they may produce spontaneous currents and magnetic moments which can be measured using appropriate experimental techniques. High $`T_c`$ cuprates present an especially interesting class of materials for studying mixed symmetry states. In these systems, recently developed technology  provides an opportunity to fabricate different structures with controllable characteristics, such as $`\pi `$-junctions  (junctions with equilibrium phase difference equal to $`\pi `$), submicron size $`\varphi _0`$-junctions  (with an equilibrium phase difference $`\varphi _0`$ which is different from 0 or $`\pi `$ ), $`\pi `$-SQUIDs  and $`\pi /2`$-SQUIDs , and superconducting qubit prototypes . Therefore, a quantitative prediction of properties of such finite size unconventional superconducting systems becomes extremely important for both experiments and technology.

The mixed symmetry states can be realized near a surface of a $`d`$-wave superconductor or at a grain boundary junction between two different superconductors. In the grain boundary junction case, the symmetry of the mixed state is dictated by the symmetries of the order parameters on both sides of the junction (proximity effect). On the other hand, at a surface, a subdominant order parameter may appear due to an attractive interaction potential in the corresponding channel. This may involve a second order phase transition below the superconducting transition temperature . The mixed symmetry states that have been suggested for $`d`$-wave superconductors are $`d_{x^2y^2}+id_{xy}`$ and $`d_{x^2y^2}+is`$  (we denote them $`d+id^{}`$ and $`d+is`$, respectively). Which of these states is realized near the surface of a $`d`$-wave superconductor is still an open question. An identification of this state would provide useful information about microscopic interactions in the superconductor. The $`d+id^{}`$ and $`d+is`$ states are inherently different in the sense that the Cooper pairs in the former case have intrinsic magnetic moment, while in the latter, they do not . Therefore, it is natural to expect differences in some measurable quantities, such as spontaneous current and magnetic field.

In this article, we perform self-consistent calculations of spontaneous current densities and magnetic field distributions for finite size systems with different mixed symmetries. We show that these experimentally observable characteristics are very dependent on the time-reversal symmetry properties of the system. For the $`d+id^{}`$ state, we find that vortices appear, and we connect their appearance with the chiral nature of the order parameter in this state.

Most general approaches to calculation of currents in superconducting systems are based on Gorkov equations for Green’s functions of the superconductor. It is widely accepted that this “mean field” approach is valid on a phenomenological level , independent of further developments in the first principles’ theory of high $`T_c`$ superconductivity. In a quasiclassical limit, these equations (which are also called Eilenberger equations, see Ref. ) can be solved by integrating along the classical trajectories of the quasiparticles (with boundary conditions at infinity). In finite size systems, where the trajectories undergo multiple reflections from surfaces and interfaces, the problem of definition of the boundary conditions becomes complicated. We demonstrate how a modification of the Eilenberger formalism based on the Schopohl-Maki transformation  can be used for stable numerical calculations in finite size two-dimensional (2D) systems.

The Eilenberger equations for quasiclassical Green’s functions can be written as

$$𝐯_F\widehat{g}+[\omega \widehat{\tau }_3+\widehat{\mathrm{\Delta }},\widehat{g}]=\widehat{0},\widehat{g}^{2}=\widehat{1},$$
(1)
where $`\omega `$ is the Matsubara frequency and

$`\widehat{\tau }_3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\widehat{g}=\left(\begin{array}{cc}g& f\\ f^{}& g\end{array}\right),\widehat{\mathrm{\Delta }}=\left(\begin{array}{cc}0& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& 0\end{array}\right).`$

The matrix Green’s function $`\widehat{g}`$ and the superconducting order parameter $`\mathrm{\Delta }`$ are both functions of the Fermi velocity $`𝐯_F`$ and the position $`𝐫`$. $`\mathrm{\Delta }`$ is determined by the self-consistency equation, which in the two-dimensional (2D) case can be written as

$`\mathrm{\Delta }(\theta )=2\pi N(0)T{\displaystyle \underset{\omega >0}{}}V_{\theta \theta ^{}}f(\theta ^{})_\theta ^{}`$ (2)
where $`\theta `$ is the angle between $`𝐯_F`$ and the $`x`$-axis, $`V_{\theta \theta ^{}}`$ is the interaction potential, $`N(0)`$ is the density of states at the Fermi surface, and $`\mathrm{}_\theta `$ represents averaging over $`\theta `$. We also have $`\mathrm{\Delta }^{}=\mathrm{\Delta }^{}`$, which holds for any singlet order parameter symmetry. Generally, it is possible to obtain an order parameter which is a mixture of several terms with different symmetries, e.g., $`\mathrm{\Delta }=\mathrm{\Delta }_{x^2y^2}+\mathrm{\Delta }_{xy}+\mathrm{\Delta }_s`$, where $`\mathrm{\Delta }_{x^2y^2}=\mathrm{\Delta }_1\mathrm{cos}2\theta ,`$ $`\mathrm{\Delta }_{xy}=\mathrm{\Delta }_2\mathrm{sin}2\theta `$, and $`\mathrm{\Delta }_s`$ are the dominant $`d_{x^2y^2}`$ component, the subdominant $`d_{xy}`$, and the $`s`$ components of the order parameter, respectively. The corresponding interaction potential,

$$V_{\theta \theta ^{}}=V_{d1}\mathrm{cos}2\theta \mathrm{cos}2\theta ^{}+V_{d2}\mathrm{sin}2\theta \mathrm{sin}2\theta ^{}+V_s,$$
(3)
should be substituted in the self-consistency equation (2) in order to assure a self-consistency of the solution . The current density $`𝐣(𝐫)`$ is given by

$`𝐣=4\pi ieN(0)T{\displaystyle \underset{\omega >0}{}}𝐯_Fg_\theta .`$ (4)
For numerical calculations, it is convenient to parameterize the quasiclassical Green’s functions via

$$g=\frac{1ab}{1+ab},f=\frac{2a}{1+ab},f^{}=\frac{2b}{1+ab}.$$
(5)
Functions $`a`$ and $`b`$ satisfy two independent, but nonlinear, (Riccati) equations:

$`𝐯_Fa`$ $`=`$ $`\mathrm{\Delta }\mathrm{\Delta }^{}a^22\omega a,`$ (6)
$`𝐯_Fb`$ $`=`$ $`\mathrm{\Delta }^{}\mathrm{\Delta }b^22\omega b.`$ (7)
From these equations it follows that $`a(𝐯_F)=b^{}(𝐯_F)`$ and $`b(𝐯_F)=a^{}(𝐯_F)`$. The solutions of these equations can be interpreted as trajectories of quasiparticles. One should integrate these equations along all possible trajectories and perform the summation over the trajectories to calculate the current. To find $`a`$ and $`b`$ along the trajectories, one needs to use boundary conditions at the ends of the trajectories. In infinite systems, one usually assumes that all trajectories go deep into the bulk of the superconductor, i.e., one can use the bulk solutions

$$a_\pm =\frac{\mathrm{\Delta }}{\omega \pm \mathrm{\Omega }},b_\pm =\frac{\mathrm{\Delta }^{}}{\omega \pm \mathrm{\Omega }},$$
(8)
with $`\mathrm{\Omega }=\sqrt{\omega ^2+|\mathrm{\Delta }|^2}`$, as the boundary conditions at infinity. In the case of finite size systems where multiple reflections from surfaces and interfaces are possible, stability of the numerical procedure for calculating the current is not obvious because of complications with choosing a proper boundary condition for a given trajectory. Nevertheless, we show that a stable numerical procedure can still be developed.

Numerical integration for $`a`$ ($`b`$) in Eq. (7) should be taken in the direction of $`𝐯_F`$ (-$`𝐯_F`$) to ensure stability. When $`\mathrm{\Delta }`$ is a constant, the solution for $`a`$ can be written analytically,

$`a_f`$ $`=`$ $`a_++{\displaystyle \frac{a_ia_+}{1+\frac{\mathrm{\Delta }^{}}{\mathrm{\Omega }}(a_ia_{})e^{\mathrm{\Omega }\tau }\mathrm{sinh}\mathrm{\Omega }\tau }}`$ (9)
$``$ $`a_++{\displaystyle \frac{\mathrm{\Omega }}{\mathrm{\Delta }^{}}}\left({\displaystyle \frac{a_ia_+}{a_ia_{}}}\right)e^{2\mathrm{\Omega }\tau }(\mathrm{if}\mathrm{\Omega }\tau 1),`$ (10)
where $`a_i`$ and $`a_f`$ are the values of $`a`$ at the initial ($`𝐫_i`$) and final ($`𝐫_f`$) points of the trajectory, and $`\tau =|𝐫_f𝐫_i|/v_F`$ is the migration time between the initial and final points. It is clear that the solution for $`a`$ relaxes to the bulk value $`a_+`$ at the distance $`L=𝐯_F/2\mathrm{\Omega }`$ which is of the order of the coherence length $`\xi _0`$. In other words, when the quasiparticle moves away from the initial point at a distance of a few $`\xi _0`$’s, any information about the initial point $`a_i`$ is lost. The same argument is valid for the function $`b`$.

Let us now consider a system with a restricted geometry. After integrating over a few $`\xi _0`$’s (considering the reflections), $`a_f`$ becomes almost independent of $`a_i`$. This solution corresponds to a simple exponential relaxation of the functions $`a`$ and $`b`$ to their local “steady-state” values (which can be different from the bulk values) defined by the local values of the order parameter. Therefore, to find $`\mathrm{\Delta }`$ at a given point, one does not need the values of $`\mathrm{\Delta }`$ at distances larger than several $`\xi _0`$ along the trajectory. Such a “relaxational” property of the $`a`$ and $`b`$ functions significantly simplifies the numerical solution of the self-consistent 2D problem. All the “multiple-reflection history” of a trajectory becomes a moot point, and for practical calculations, this trajectory can be cut at a distance of a few $`\xi _0`$’s from a given point. We used this observation in our self-consistent calculations. For the “cutting” distance, we choose $`10\xi _0`$$`20\xi _0`$. We set the bulk values of $`a`$ and $`b`$ ($`a_+`$ and $`b_+`$) calculated at the initial point of the “truncated” trajectory as the boundary conditions (such a choice does not affect the final results because the system has no memory) and integrate along the trajectory until we get to the point where we calculate the current. We found that the results are very stable and do not depend on the value of the “cutting” distance.

To compare the spontaneous currents and magnetic fields in the $`d+is`$ and $`d+id^{}`$ states, we performed self-consistent calculations of the order parameter in small size triangular regions of clean $`d`$-wave superconductor with different subdominant order parameters. We assume specular boundaries with the length of the longer edge of the triangles being $`30\xi _0`$. Figs. 1 and 2 show the spontaneous current and magnetic field distributions in systems with the two different symmetries. We imposed the desired subdominant symmetry by introducing an additional attractive interaction potential in the corresponding channel \[$`V_{d2}`$ or $`V_s`$ in Eq. (3)\]. The subdominant critical temperatures are taken to be $`0.5T_c`$ in both cases. The orientation of the dominant order parameter makes a $`45^{}`$ angle with respect to the longer (right) edge of the triangle (see the figures). As a result, only this boundary of the triangle is pair breaking; the quasiparticles face a different sign of the dominant order parameter after reflection from the surface . The suppression of the main order parameter at the surface allows the subdominant order parameter to appear near the surface.

The spontaneous current and magnetic field distributions are evidently different in Figs. 1 and 2. In the first case, the current makes a counterclockwise loop throughout the region, while in the second case, the direction of rotation of the current changes throughout the triangle. Especially, appearance of the two vortices in Fig. 2 reflects the chiral nature of the $`d+id^{}`$ symmetry of the order parameter . There is no phase winding around these vortices and therefore no flux quantization is expected. In fact, the flux trapped in these vortices is much smaller than a flux quantum. The magnitude of the magnetic field at the maximum positions is of the order of $`10^4`$$`10^3G`$ in both cases. This value agrees with the magnitude of the magnetic field reported in Ref. . In the $`d+id^{}`$ case, the magnetic field is strongly peaked at the vortices. The magnetic field in other spatial points is noticeable only in the vicinity of the pair-breaking edge (within a few coherence lengths) and is almost one order of magnitude weaker than the field at the vortex peaks. This is a result of the superscreening effect which happens when the spontaneous current is due to the intrinsic angular moment of the order parameter . The total flux generated by this magnetic field is very small and difficult to measure. This is in agreement with the magnitude of the flux, calculated using a different technique, in Ref. . In the $`d+is`$ case, on the other hand, the peak is spread over a wider region. In fact, in larger systems, the size of this region could be of the order of the penetration depth. This is because in the $`d+is`$ case, the superscreening effect is absent and the only mechanism to return the current is the Meissner effect, which happens over the length scale of a penetration depth. The flux generated by this magnetic field can be large and, therefore, measurable. This suggests that the flux measured in Ref.  can be generated by a $`d+is`$ symmetry breaking state and not by a $`d+id^{}`$ one.

The spontaneous magnetic moments of mesoscopic islands with $`d+is`$ or $`d+id^{}`$ symmetry can be in principle used as qubits. Since the “up” and “down” states are degenerate eigenstates of the effective Hamiltonian, tunneling between them becomes only possible in the presence of external magnetic field, field gradient, etc., depending on the multipole structure of the spontaneous current distribution. In order to determine the tunneling amplitude, a three-dimensional picture should be considered, which is beyond the scope of this publication.

We also have calculated the spontaneous current and magnetic field distributions in grain boundary junctions. The geometry of these systems consists of a square of size $`=30\xi _0\times 30\xi _0`$, divided into two equal parts separated by a grain boundary junction. The right half of the square is a clean $`d`$-wave superconductor with a $`45^{}`$ orientation of the order parameter with respect to the boundaries. The left half, on the other hand, can be either a clean $`s`$-wave superconductor ($`s`$-$`d`$ junction, Fig. 3), or a clean $`d`$-wave superconductor with a $`0^{}`$-orientation ($`d`$-$`d`$ junction, Fig. 4). The grain boundary is taken to be perfectly transparent with no roughness or faceting . Although this choice of grain boundary does not exactly correspond to the reality, study of it may provide useful information for real systems. To be able to compare the two systems, we assumed that the $`s`$-wave and $`d`$-wave superconductors have the same transition temperatures. We did not introduce any subdominant order parameter here. However, we introduced an additional phase difference of $`\mathrm{\Delta }\varphi =\pi /2`$ between the two sides. Such a choice corresponds to the equilibrium phase difference of the junction at which the total current passing through the junction is zero . Calculations are done at $`T=0.1T_c`$.

The results of spontaneous current density and magnetic field distributions are shown in Fig. 3 for the $`s`$-$`d`$ junction, and Fig. 4 for the $`d`$-$`d`$ junction. Note that the current distribution is not symmetric with respect to the grain boundary. For the $`s`$-$`d`$ junction of Fig. 3, the current has a maximum at the grain boundary and returns through the bulk of the superconductors on both sides. If the system is large, the size of the region where the spontaneous current is non-zero, should be of order of the penetration depth. For the $`d`$-$`d`$ junction (Fig. 4), on the other hand, the current changes the direction just within a few coherence length from the boundary, again due to the superscreening effect . Notice also that near the edges of the system, on the left side of the grain boundary junction ($`0^{}`$-orientation of the order parameter), the current returns along the diagonal, whereas on the right side ($`45^{}`$-orientation) it forms two small vortices and antivortices. These vortices, which are absent in the $`s`$-$`d`$ case, again reflect the chiral nature of the $`d+id^{}`$ symmetry; although the subdominant order parameter is absent in this case, the correlation functions convey this symmetry near the grain boundary junction due to the proximity effect . The magnetic field distributions, displayed in Figs. 3 and 4, is of the order of $`10^4G`$ at the positions of the maxima, in both $`s`$-$`d`$ and $`d`$-$`d`$ junctions. In the case of $`d`$-$`d`$ junction (Fig. 4), magnetic field is peaked at the location of vortices. Away from the vortices, the magnetic field is localized near the grain boundary junction and almost one order of magnitude smaller. The flux generated by such a magnetic field is very small. Thus, observation of large flux as in Ref.  can not be associated with such an effect. In real systems, due to the faceting effect , some $`\pi `$-loops may exist along the grain boundary junction. They can produce large fluxes, of the same order as observed in Ref. .

It is important to emphasize here that the existence of the spontaneous current at the grain boundary junctions does not depend on the presence of a subdominant order parameter (unlike in the previous case, unless if we consider other effects ). Thus, experiments similar to that of Ref.  (that only probes the mixing of the symmetry of the order parameter) can not exclude the possibility of the existence of spontaneous current. Indeed, addition of a subdominant order parameter will suppress the spontaneous current at the boundary (see Refs. ).

One should note that in the presence of magnetic field, the procedure described here is not valid, because a path dependent phase will be accumulated to $`a`$ and $`b`$ functions, and the relaxation mechanism along the trajectory may no longer hold. This does not mean that we cannot use our approach to calculate currents and fields in realistic systems without calculating the magnetic field due to spontaneous currents self-consistently. Indeed, it is not difficult to show that the corresponding Doppler shift is rather small. We saw that for systems with size $`L30\xi _05\times 10^4`$Å, spontaneous magnetic field is always smaller than $`10^3`$G. Therefore, the induced superfluid momentum is at most $`p_s(e/c)HL10^{26}\mathrm{g}\mathrm{cm}/\mathrm{s}`$. If we take $`\mathrm{}v_F1`$ eV$``$Å for the Fermi velocity in high-Tc cuprates, we obtain $`p_sv_F`$ 0.1 $`\mu `$eV $``$ 1 mK, which is negligible compared to other energy scales in the problem. However, these restrictions may be important if the system is placed in a strong external magnetic field.

To summarize, we described self-consistent calculations of some equilibrium properties in finite size superconducting systems. We calculated distributions of spontaneous current and magnetic field in different small samples in the presence of mixed order parameter symmetries. The nature of the mixed symmetry state in all cases affects the shape of the current and magnetic field distributions. In particular, the chiral nature of the $`d+id^{}`$ states exhibits itself through the appearance of vortices close to the edges of the system. The vortices are absent in the $`d+is`$ cases. The method described here is quite general and can be applied to any 2D geometry with proper boundary conditions as long as external magnetic field is not present. The self-generated magnetic field by the spontaneous currents, however, is usually very small, i.e., it can be neglected in most of the practically important cases.

Many discussions with A. Maassen van den Brink, G. Rose, and A. Yu. Smirnov are greatly acknowledged.


text:
# 1 Introduction

## 1 Introduction

In regard to General Relativity (GR), the experience shows that designing and effectively performing experiments aimed to test it is very difficult and often expensive. This is particularly true for its linearized weak–field and slow–motion approximation which is adequate in the Solar System. The relativistic effects quite often fall below or lie just at the edge of the experimental sensitivity of various techniques which could be used. Moreover, there are lots of biases and aliasing effects induced by a host of classical phenomena which pose severe constraints to the obtainable accuracy in space-based experiments. In many cases the systematic errors are larger than the standard statistical errors resulting from the least-squares procedures with which the relativistic effects are commonly extracted from the data.

So, it is of the utmost importance

* To design experiments capable of effectively testing some post-Newtonian effects not yet directly checked at relatively low expense of time and money
* To assess as more reliably and clearly as possible the error budget so to single out the various aliasing competing forces and to evaluate the impact of their systematic errors
* To investigate more thoroughly and reanalyze the role played by the already existing or proposed experiments in order to correctly evaluate their importance

### 1.1 The importance of the present Lense-Thirring LAGEOS experiment

Nowadays, a first attempt to measure the Lense-Thirring drag of the orbit of a test body \[Ciufolini and Wheeler, 1995\] in the gravitational field of the Earth, by using the following linear combination of orbital residuals of the rates of the node $`\delta \dot{\mathrm{\Omega }}`$ of LAGEOS and LAGEOS II and the perigee $`\delta \dot{\omega }`$ of LAGEOS II \[Ciufolini, 1996\]

$$\delta \dot{\mathrm{\Omega }}^\mathrm{I}+c_1\delta \dot{\mathrm{\Omega }}^{\mathrm{II}}+c_2\delta \dot{\omega }^{\mathrm{II}}\mu _{\mathrm{LT}}60.2,$$
($`1`$)
has been reported in \[Ciufolini et al., 1998; Ciufolini, 2002\]. In eq.($`1`$) $`c_10.295`$, $`c_20.35`$ are coefficients suitably designed in order to cancel out every mismodelled contributions of the first two even ($`l=2,4`$) zonal ($`m=0`$) harmonics $`J_2`$ and $`J_4`$ of the multipolar expansion of the non-spherical Earth’s gravitational field and are built up with the orbital parameters of LAGEOS and LAGEOS II (see section 3 for further details). Finally, $`\mu _{\mathrm{LT}}`$ is the parameter in terms of which the Lense-Thirring drag is expressed: it is 1 in GR and 0 in Newtonian mechanics. GR predicts for eq.($`1`$) a linear trend with a slope of 60.2 milliarcseconds per year (mas/y). Eq.($`1`$) would allow to obtain a claimed total experimental error of the order of $`20\%`$ \[Ciufolini et al., 1998; Ciufolini, 2002\]. The important GP–B mission \[Everitt et al., 2001\], which is aimed to the detection of a complementary gravitomagnetic effect on freely falling gyroscopes and which is scheduled to fly at the beginning of 2003, would reach a claimed accuracy level of 1$`\%`$. It is worth noting that the strategy of the combined residuals of LAGEOS and LAGEOS II could allow to constraining strongly the hypothesis of a Yukawa–like fifth force \[Iorio, 2002a\] and to measure in the field of the Earth the gravitoelectric relativistic perigee advance \[Iorio et al., 2002a\].

What is the real importance of the observable of eq.($`1`$)? Is it the only possibility we have at present in order to measure the Lense-Thirring drag of the orbit of a test body in the terrestrial space environment \[Casotto et al., 1990\]? Would it be possible to adopt alternative combinations involving orbital elements of other SLR satellites? We will try to answer these questions in the following.

In Table 1 we quote the orbital parameters of the existing SLR geodetic satellites and of the proposed LARES. In it $`a,e,i`$ and $`n`$ are the semimajor axis, the eccentricity, the inclination and the Keplerian mean motion, respectively.

The paper is organized as follows. In section 2 we analyze the systematic errors affecting such measurement. In section 3 we look for alternative combinations including the orbital residuals of other existing SLR satellites. Section 4 is devoted to the conclusions.

## 2 The systematic errors

The systematic errors affecting such measurement are of two types

* Long-periodic aliasing harmonics of gravitational (solid Earth and ocean tides \[Iorio, 2001; Iorio and Pavlis, 2001; Pavlis and Iorio, 2002\]) and non-gravitational (direct solar radiation pressure, Earth albedo, direct Earth IR radiation pressure, thermal effects \[Lucchesi, 2001; 2002\]) origin. If their periods $`P`$ are shorter than the observational time span $`T_{\mathrm{obs}}`$ they can be viewed as empirically fitted quantities and can be removed from the signal. If $`P=kT_{\mathrm{obs}},k=1,2,\mathrm{}`$ they average out. If $`P>T_{\mathrm{obs}}`$, as is the case for the insidious 18.6-year lunar tide acting on the node and the perigee, the $`K_1,l=3,m=1,p=1,q=1`$ ocean tide, whose period amounts to 5.07 years, acting on the perigee of LAGEOS II, and a harmonic constituent due to the direct solar radiation pressure, whose period amounts to 11.6 years, acting on the perigee of LAGEOS II as well111Notice that the period of the 18.6–year tide depends only on the luni–solar variables, while the periods of the $`K_1`$ tide and of the direct solar radiation pressure harmonic depend on the orbital geometry of the satellite., they may alias the relativistic trend acting as superimposed linear signals. However, the 18.6-year tide is an even zonal perturbation which can be canceled out by designing suitably the observable, as done with eq.($`1`$). The impact of the other two mismodelled harmonics is reduced by the small value of the coefficient $`c_2`$ with which the perigee of LAGEOS II enters the combination of eq.($`1`$)
* Secular aliasing trends of gravitational (static even zonal harmonics of the geopotential \[Kaula, 1966\]) and non-gravitational (some tiny thermal effects \[Lucchesi, 2001b\]) origin. The impact of the mismodelled even zonal harmonics is far the most relevant in the error budget. Indeed, the systematic error $`\delta \mu _{\mathrm{LT}}`$ induced by them cannot be suitably handled or removed from the signal. One can only assess as more reliably as possible its impact on the measurement of $`\mu _{\mathrm{LT}}`$.

Eq.($`1`$) cancels out the first two even zonal harmonics: the effect of the remaining higher degree mismodelled harmonics $`\delta J_6,\delta J_8,\mathrm{}`$ amounts to 12.9$`\%`$ according to the covariance matrix of EGM96 Earth gravity model \[Lemoine et al., 1998\] up to degree $`l=20`$222In regard to this point, it should be mentioned that, at present, there is no full consensus in the scientific community on the reliability of such estimate which might be rather optimistic \[Ries et al., 1998\]. Indeed, it would not be entirely correct to automatically extend the validity of the covariance matrix of EGM96, which is based on a multi–year average that spans the 1970, 1980 and early 1990 decades, to any particular time span like that, e.g., of the LAGEOS–LAGEOS II analysis which extends from the middle to the end of the 1990 decade. Indeed, there would not be assurance that the errors in the even zonal harmonics of the geopotential during the time of the LAGEOS–LAGEOS II experiment remained correlated exactly as in the EGM96 covariance matrix, in view of the various secular, seasonal and stochastic variations that we know occur in the terrestrial gravitational field and that have been neglected in the EGM96 solution..

## 3 The search for alternative combinations

A possible strategy for improving the accuracy of the present–day Lense–Thirring LAGEOS–LAGEOS II experiment consists of suitable combinations of the orbital residuals $`\delta \dot{\mathrm{\Omega }}`$ and $`\delta \dot{\omega }`$ of the rates of the nodes and the perigees of different SLR satellites. Such combinations can be written in the form

$$\underset{i=1}{\overset{N}{}}c_if_i=X_{\mathrm{GR}}\mu _{\mathrm{GR}},$$
($`2`$)
in which the coefficients $`c_i`$ are, in general, suitably built up with the orbital parameters of the satellites entering the combinations, the $`f_i`$ are the residuals $`\delta \dot{\mathrm{\Omega }},\delta \dot{\omega }`$ of the rates of the nodes and the perigees of the satellites entering the combination, $`X_{\mathrm{GR}}`$ is the slope, in mas/y, of the general relativistic trend of interest and $`\mu _{\mathrm{GR}}`$ is the solve–for parameter, to be determined by means of usual least–square procedures, which accounts for the general relativistic effect. For example, in the case of the Lense–Thirring LAGEOS–LAGEOS II experiment \[Ciufolini, 1996\] $`X_{\mathrm{LT}}=60.2`$ mas/y. More precisely, the combinations of eq.($`2`$) are obtained in the following way. The equations for the residuals of the rates of the $`N`$ chosen orbital elements are written down, so to obtain a non homogeneous algebraic linear system of $`N`$ equations in $`N`$ unknowns. They are $`\mu _{\mathrm{GR}}`$ and the first $`N1`$ mismodelled spherical harmonics coefficients $`\delta J_l`$ in terms of which the residual rates are expressed \[Iorio, 2002c\]. The coefficients $`c_i`$ and, consequently, $`X_{\mathrm{GR}}`$ are obtained by solving for $`\mu _{\mathrm{GR}}`$ the system of equations. So, the coefficients $`c_i`$ are calculated in order to cancel out the contributions of the first $`N1`$ even zonal mismodelled harmonics which represent the major source of uncertainty in the Lense–Thirring precessions. The coefficients $`c_i`$ can be either constant 333In general, the coefficient of the first orbital element entering a given combination is equal to 1, as for the combinations in \[Ciufolini, 1996; Iorio, 2002b; Iorio et al., 2002a; 2002b\].  or depend on the orbital elements of the satellites entering the combinations.

Now we expose how to calculate the systematic error due to the mismodelled even zonal harmonics of the geopotential for the combinations involving the residuals of the nodes and the perigees of various satellites.

In general, if we have an observable $`q`$ which is a function $`q=q(x_j)`$, $`j=1,2\mathrm{}M`$ of $`M`$ correlated parameters $`x_j`$ the error in it is given by

$$\delta q=\left[\underset{j=1}{\overset{M}{}}\left(\frac{q}{x_j}\right)^2\sigma _j^2+2\underset{hk=1}{\overset{M}{}}\left(\frac{q}{x_h}\right)\left(\frac{q}{x_k}\right)\sigma _{hk}^2\right]^{\frac{1}{2}}$$
($`3`$)
in which $`\sigma _j^2C_{jj}`$ and $`\sigma _{hk}^2C_{hk}`$ where $`\{C_{hk}\}`$ is the square matrix of covariance of the parameters $`x_j`$.

In our case the observable $`q`$ is any residuals’ combination

$$q=\underset{i=1}{\overset{N}{}}c_if_i(x_j),j=1,2\mathrm{}10,$$
($`4`$)
where $`x_j,j=1,2\mathrm{}10`$ are the even zonal geopotential’s coefficients $`J_2,J_4\mathrm{}J_{20}`$. Since

$$\frac{q}{x_j}=\underset{i=1}{\overset{N}{}}c_i\frac{f_i}{x_j},j=1,2\mathrm{}10,$$
($`5`$)
by putting eq.($`5`$) in eq.($`3`$) one obtains, in mas/y

$$\delta q=\left[\underset{j=1}{\overset{10}{}}\left(\underset{i=1}{\overset{N}{}}c_i\frac{f_i}{x_j}\right)^2\sigma _j^2+2\underset{hk=1}{\overset{10}{}}\left(\underset{i=1}{\overset{N}{}}c_i\frac{f_i}{x_h}\right)\left(\underset{i=1}{\overset{N}{}}c_i\frac{f_i}{x_k}\right)\sigma _{hk}^2\right]^{\frac{1}{2}}.$$
($`6`$)
The percent error, for a given general relativistic trend and for a given combination, is obtained by taking the ratio of eq.($`6`$) to the slope in mas yr-1 of the general relativistic trend for the residual combination considered.

The validity of eq.($`6`$) has been checked by calculating with it and the covariance matrix of the EGM96 gravity model up to degree $`l=20`$ the systematic error due to the even zonal harmonics of the geopotential of the Lense–Thirring LAGEOS–LAGEOS II experiment; indeed the result

$$\delta \mu _{\mathrm{LT}}=13\%\mu _{\mathrm{LT}}$$
($`7`$)
claimed in \[Ciufolini et al., 1998\] has been obtained again444It may be interesting to notice that if, with a more conservative approach, we consider only the diagonal part of the EGM96 covariance matrix, up to degree $`l=20`$, the systematic error due to the mismodelled even zonal harmonics of the geopotential for the Lense–Thirring LAGEOS–LAGEOS II experiment amounts to 46.5$`\%`$ \[Iorio, 2002c\].. For the systematic errors due to the even zonal harmonics of the geopotential of alternative proposed gravitomagnetic and gravitoelectric experiments, according to EGM96, see \[Iorio, 2002b; Iorio et al., 2002a; 2002b\]. It is worth noticing that, since the orbits of the LAGEOS satellites are almost insensitive to the geopotential’s terms of degree higher than $`l=20`$, the estimates based on the covariance matrix of the EGM96 gravity model up to degree $`l=20`$ should be considered rather reliable although the higher degree terms of EGM96 might be determined with a low accuracy.

It should also be pointed out that the evaluations of the systematic error due to geopotential based on the approach of the combined orbital residuals should be free from some uncertainties due to possible secular, seasonal and stochastic effects. Indeed, even putting aside the fact that most of the secular and seasonal variations of the geopotential are concentrated just in the first even zonal harmonics, if we cancel out as many even zonal harmonics as possible, the uncertainties in the evaluation of the error based on the remaining correlated even zonal harmonics of higher degree should be greatly reduced irrespectively of the chosen time span.

A very important point to stress is that the forthcoming new data on the Earth’s gravitational field by CHAMP, which has been launched in July 2000, and GRACE, which has been launched in March 2002, will have a great impact on the reduction of the systematic error due to the mismodelled part of geopotential and on the reliability of the estimates based on it.

### 3.1 The combinations without the perigee of LAGEOS II

In regard to alternative combinations, we first explore the possibility of reducing $`\delta \mu _{\mathrm{LT}}`$ by adopting different residual combinations obtained by substituting the perigee of LAGEOS II555The orbits of the other SLR satellites are far less elliptical than that of LAGEOS II (except for Starlette), so that their perigees are not good observables. It turns out that the inclusion of the perigee of Starlette would not yield any improvements. with one or more nodes of Ajisai, Starlette, Stella, WESTSPAC -1, ETALON -1 and ETALON -2. The advantages of using only the nodes would be

* No odd zonal geopotential harmonics
* No odd degree ocean tidal perturbations
* Non -gravitational perturbations of smaller magnitude and fairly well modelled
* Atmospheric drag negligible not only for the LAGEOS and the ETALON satellites, but also for Starlette and Ajisai orbiting at lower altitudes
* Smaller rms on the experimental residuals

The exclusion of the perigee of LAGEOS II leads to two classes of residual combinations.

* Combinations without the ETALON satellites. They should be discarded because the error due to the uncancelled higher degree even zonal harmonics of the geopotential, according to EGM96, is itself far larger than 20$`\%`$, which is the estimated total error in the LAGEOS–LAGEOS II Lense–Thirring measurement. It is so because of the fact that the satellites orbiting at lower altitudes than LAGEOS and LAGEOS II are more sensitive to the higher degree even zonal harmonics of the geopotential.
* Combinations with the ETALON satellites. Due to their higher altitude with respect to the LAGEOS satellites $`\delta \mu _{\mathrm{LT}}`$ falls, in some cases, below 1$`\%`$. However, they should not be considered because the coefficients of the nodes of the ETALON satellites are, in these cases, of the order of $`10^3`$ so that they greatly amplify every perturbations acting on such elements. This is particularly true for the solid Earth tesseral $`K_1`$ tide. Indeed, it turns out that its $`l=2,m=1,p=1,q=0`$ constituent induces on ETALON–1 a perturbation with a nominal amplitude of 1,216.57 mas and period of 10,880.8 days while for ETALON- 2 we have a nominal amplitude of 1,269.35 mas and a period of 11,130.1 days. Even if the solid Earth tides are known at a 0.5$`\%`$ level \[Iorio, 2001\], the effect of the too large weighing coefficients and the very long periods of such tidal perturbations would induce insidious mismodelled linear trends over reasonable observational time spans of few years. However, it must pointed out that a practical useful use of the SLR data of the ETALON satellites would be problematic mainly because the Russian satellites are, at present, very poorly tracked. Moreover, they wobble and their center of mass is not well defined. Last but not least, the accuracy with which they were constructed and the pre launch ground controls are not comparable to those of LAGEOS satellites.

In conclusion, even if very appealing, the possibility of using residual combinations involving only the nodes of the existing geodetic SLR satellites should be rejected.

### 3.2 The combinations with the perigee of LAGEOS II

It turns out that it is possible to obtain admissible results only by involving the perigee of LAGEOS II.

* Combinations with ETALON -1. We obtained seventeen combinations which, in principle, could be considered because $`\delta \mu _{\mathrm{LT}}`$ $`<`$ 20$`\%`$ for them. However, it turns out that the combinations including the perigee of LAGEOS II and the node of ETALON–1 cannot represent a genuine improvement with respect to the present day LAGEOS–LAGEOS II Lense–Thirring experiment because $`\delta \mu _{\mathrm{LT}}13\%`$ for them. The systematic errors due to the time–dependent non–gravitational perturbations would make the total error budget unfavorable with respect to the current experiment.
* Combinations without ETALON -1. We obtained eight combinations involving only the LAGEOS satellites and the lower altitude satellites. It turns out that the combination

$$\delta \dot{\mathrm{\Omega }}^\mathrm{I}+c_1\delta \dot{\mathrm{\Omega }}^{\mathrm{II}}+c_2\delta \dot{\mathrm{\Omega }}^{\mathrm{Aj}}+c_3\delta \dot{\omega }^{\mathrm{II}}\mu _{\mathrm{LT}}61.2,$$
($`8`$)
with $`c_1=0.443,c_2=0.0275,c_3=0.341`$ is the only one that could lead to a real, slight improvement of the error budget666It may be interesting to notice that if, with a more conservative approach, we consider the diagonal part only of the EGM96 covariance matrix, up to degree $`l=20`$, the systematic error due to the mismodelled even zonal harmonics of the geopotential for such alternative combination amounts to 64.2$`\%`$ \[Iorio, 2002c\]. \[Iorio, 2002b\]. Indeed

+ The systematic error due to the even zonal harmonics amount to $`\delta \mu _{\mathrm{LT}}10.7\%`$ since eq.($`8`$) cancels out the effect of $`\delta J_6`$ as well
    + The estimates relative to the periodic time–varying perturbations for eq.($`1`$) would be valid also in this case because the coefficients of the node and the perigee of LAGEOS II are almost the same
    + The additional perturbations induced by Ajisai are damped since the coefficient of its node amounts only to -0.0275
    + It turns out that the tidal perturbations on the node of Ajisai have periods not longer than few years and the non -gravitational perturbations are fairly well studied \[Sengoku et al., 1995; 1997\]
    + Ajisai, contrary to the ETALON satellites, is a well known and accurately tracked geodetic satellite

## 4 Conclusions

The combination of eq.($`1`$) plays a really important role in the context of the efforts aimed to the detection of the general relativistic Lense-Thirring drag in the gravitational field of the Earth. Indeed, at present, it represents a relatively fast, cheap, reliable and almost unique way for measuring this effect in the terrestrial space environment by means of the analysis of the orbital data of the existing laser–tracked satellites. This should also justify and encourage the efforts aimed to improve the error budget and the knowledge of the aliasing classical forces affecting such measurement. However, according to the covariance matrix of the EGM96 model up to degree $`l=20`$, the alternative combination of eq.($`8`$) might enforce and slightly improve the experimental accuracy of the measurement.

At the present–day level of knowledge of the perturbing forces of the terrestrial space environment, only the not yet approved LAGEOS–LARES mission, whose originally proposed configuration has been recently modified \[Iorio et al., 2002b\], could yield a notable, genuine improvement in measuring the Lense–Thirring effect with SLR.

Moreover, within few years the new data on the geopotential from the CHAMP and GRACE missions should be available. They should improve our knowledge, among other things, of the spherical harmonics coefficients of the geopotential, and, consequently, the systematic gravitational errors in all the proposed Lense–Thirring measurements should reduce as well.