text:
# Temperature dependence of the conductivity of the electronic crystal

## Abstract

We study the temperature dependence of the conductivity of the 2D electronic solid. In realistic samples, a domain structure forms in the solid and each domain randomly orients in the absence of the in-plane field. At higher temperature, the electron transport is governed by thermal activation form of $`\sigma _{xx}(T)e^{\mathrm{\Delta }_0/k_BT}`$. The impurities will localize the electron states along the edges of the crystal domains. At sufficient low temperature, another transport mechanism called Mott’s variable range hopping mechanism, similar to that in a disorder insulator takes effect. We show that as the temperature decreases, a crossover from the fixed range hopping of the transport to the variable range hopping of transport in the 2D electron system may be experimentally observed.

It was initially predicted by Wigner that two-dimensional (2D) electrons crystallize into a triangular lattice in the low density limit where the electron-electron interactions dominate over the kinetic energy. In an ideally clean 2D system, the critical $`r_s`$ was presented to be $`37\pm 5`$ from quantum Monte Carlo simulations Tanatar . A strong magnetic field perpendicular to the 2D plane can effectively localize electron wave functions while keeping the kinetic energy controlledYosh . Since this lessens the otherwise severe low-density condition, it is believed that the Wigner crystal (WC) can be stabilized in a sufficiently strong magnetic fieldLau ; Lam ; Kuk ; Santos . Approximate calculationsCM1  have shown that the WC becomes the lowest energy state when the filling factor $`\nu <1/6`$ for $`GaAs/AlGaAs`$ electron system and around $`\nu =1/3`$ for the hole system.

As well-known, while the transport behavior of the WC is characterized by non-linear I-V curves, the temperature dependence of the conductivity of a WC is believed to be normal, namely, it has an ordinary thermal activation formChui ; Li . Moreover, since the impurities pin the electronic crystal, a domain structure forms in realistic samplesSher ; Fil . While the electrons in a domain have an order as they are in the ideal crystal, the orientations of the domains are random. Another role played by the impurities is they localize the electron states along the edges of the crystal domains. The electrons may hop between the edges of the domains. In this case, the transport behavior of the system is similar to that in a disordered insulator. In this work, we will show that the transport of the electrons in the 2D electronic solid obeys a generalized Mott’s variable range hopping theorymott  for a low temperature $`TT_0`$, while it obeys the ordinary thermal activation law for $`TT_0`$. We find that $`T_0`$ is in the experimentally reachable regime if the sample parameters are properly chosen. This implies that one may experimentally observe a different temperature dependence of the conductivity at different temperature regimes. As $`T`$ varies, the conductivity, according to our calculation, may have a crossover from $`\sigma _{xx}(T)e^{A/T^{1/2}}`$ for $`TT_0`$ to a thermal activation form $`\sigma _{xx}(T)e^{\mathrm{\Delta }_0/k_BT}`$.

Pinning of the WC by impurities as a result of breaking of the translational invariance has been widely investigatedKuk ; Millis . In realistic samples a domain structure is formed due to a finite impurity density. The electrons in each domain are ordered as they are in the crystal. ShermanSher  and FilFil  had studied the angular pinning and the domain structure of the electronic crystal mediated by acoustic-phonon in $`IIIV`$ semiconductor. In the absence of the in-plane field, each domain orients randomly, just like the domains in ferromagnets. It can be shown that an in-plane field favors the domains to orient to the same directionyang . Hence the in-plane magnetic field may serve as a tunable means to probe the orientation of the crystal. An ideal electronic crystal is an insulator and the conductivity $`\sigma _{xx}e^{\mathrm{\Delta }_0/k_BT}`$. This thermal activation form of the conductivity implies that the electrons are hopping with a fixed range mechanism. It has been confirmed by experiments with $`\mathrm{\Delta }_0`$ typically of the order 1K jiang . In a realistic domain structure, however, the localized electrons may hop between the edges of the randomly oriented domains. Since the experimentally reachable temperature is as low as 10mK, the variable range hopping mechanism mott  may work in this temperature regime. In the following, we will calculate the electron conductivity according to the different electron hopping mechanisms and determine the characteristic temperature of the crossover region $`T_0`$.

In the usual Anderson localization the envelope of the wave function falls off exponentially as $`\varphi e^{r/\xi }`$, where $`\xi `$ is the localization length, while in the magnetic field the electronic wave function of a perfect system is essentially a Gaussian as $`\varphi e^{r^2/2l_B^2}`$, where $`l_B`$ is the magnetic length. In a slightly disordered system one can think that some of the states will be pinned at certain isolated impurity site. The mixing of these states due to quantum-mechanical tunnelling leads to a simple exponential tail in the wave functionLi . In a strong magnetic field, the electrons condense into a crystal at lower filling factors. The electrons are Coulomb localized. When the temperature is high enough the transport is of the thermal activation form, which implies that the electrons are hopping with a fixed range mechanismChui ; Li . The hopping range is determined by $`R_0=\sqrt{1/\pi n_I}`$, where $`n_I`$ is the impurity density. However, localized states by impurities may exist along the edges of the domains of electronic crystal. When the temperature is sufficiently low such that there is nearly no phonon with energy to assist the electron making the nearest hopping, Mott’s variable hopping mechanismmott  allows the electrons hop a larger distance $`R>R_0`$ to a state which has a smaller energy difference $`\mathrm{\Delta }(R)`$. In turn, the hopping conduction is determined by the typical decay rate of the tails of the wave function. The hopping probability is then given by

$$p\mathrm{exp}[R/\xi \mathrm{\Delta }/k_BT],$$
(1)
where $`R=|\stackrel{}{R}_i\stackrel{}{R}_j|`$ and $`\mathrm{\Delta }`$ is the activation energy.

For non-interacting electrons, Mott hopping with an approximately constant density of states at the Fermi energy gives,

$$\rho _{xx}(T)=\rho _0(T)exp(r_c/\xi )=\rho _0(T)exp(A_0/T)^{1/3},$$
(2)
where $`r_c`$ is a characteristic hopping length, which in this regime is equal to the Mott hopping length. However, this simple treatment does not work for the WC because of the strong Coulomb interaction between electronsPollak . The Coulomb gap depresses the density of states near the Fermi surface efros1 ; efros2 ; YM ; Chan . Efros et al efros  had derived the density of states near the Fermi surface $`N(E)|\mathrm{\Delta }E|=|EE_F|`$. The condition to find one state within a circle of radius $`R`$ is given by

$$\pi R^2N(E_F)\mathrm{\Delta }(R)=1.$$
(3)
Substitute $`N(E_F)`$ by the available states near the Fermi surface, i.e., $`N(E_F)\overline{N}(E_F)=\frac{1}{\mathrm{\Delta }E}_0^{\mathrm{\Delta }E}𝑑EN(E)`$, and note that $`|\mathrm{\Delta }E|1/R`$, we get

$$\mathrm{\Delta }(R)\frac{\mathrm{}v_F}{R},$$
(4)
where $`v_F`$ is the Fermi velocity. Put this $`R`$-dependent energy difference into formula (1) and maximizing $`p`$, one finds the optimal hopping range $`R=\overline{R}`$ and the maximum of the probability are given by

$`\overline{R}^2={\displaystyle \frac{\mathrm{}v_F\xi }{k_BT}},pe^{2\overline{R}/\xi }=e^{A/T^{1/2}},`$ (5)
with $`A=[\frac{4\mathrm{}v_F}{k_B\xi }]^{1/2}`$. The conductivity in the variable range hopping is thenefros1 ; efros2

$$\sigma _{xx}pe^{A/T^{1/2}}.$$
(6)
The characteristic temperature $`T_0`$ above which the fixed range hopping dominates is determined by $`\overline{R}=R_0`$, namely

$`k_BT_0=\pi n_I\mathrm{}v_F\xi =\pi n_I{\displaystyle \frac{\mathrm{}^2\xi }{m_bl_B}}.`$ (7)
In a strong magnetic field the decay length is comparable to the cyclotron radius $`\xi R_c`$Ebert ; Fogler . We find for sample with $`n_I1.0\times 10^8`$cm<sup>-2</sup>, $`T_040`$mK . This temperature is experimentally reachable. We anticipate the different dependence of the conductivity in different temperature regimes can be observed in future experiments.

Now, we briefly discuss the effect of the tilted field. Consider an electron moving on a $`x`$-$`y`$ plane under the influence of a strong magnetic field which is tilted an angle $`\theta `$ to the normal, with $`\stackrel{}{B}=(B\mathrm{tan}\theta ,0,B)`$. The electron is confined in a harmonic potential $`V(z)=\frac{1}{2}m_b\mathrm{\Omega }^2z^2`$ in the $`z`$-direction, where $`m_b`$ is the band mass of the electron and $`\mathrm{\Omega }`$ the characteristic frequency. Such a quantum well has been chosen to deal with many quantum Hall systemsPhillips ; Jungwirth ; yu  to substitute the realistic potential which is either triangular cha  or squaremac . We work in the ”Landau gauge” by choosing the vector potential $`\stackrel{}{A}=\{0,xB_zzB_x,0\}`$. The single particle wave function for the lowest LL are:

$`\varphi _X(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{L_y}}}e^{\frac{iXy}{l_B^2}}\mathrm{\Phi }_0^{\omega _+}((xX)\mathrm{sin}\stackrel{~}{\theta }+z\mathrm{cos}\stackrel{~}{\theta })`$ (8)
$`\times \mathrm{\Phi }_0^\omega _{}((xX)\mathrm{cos}\stackrel{~}{\theta }z\mathrm{sin}\stackrel{~}{\theta }),`$
where $`l_B`$ is the magnetic length and $`X`$ is an integer multiple of $`2\pi l_B^2/L`$. $`\mathrm{\Phi }_0^{\omega _\pm }`$ is the harmonic oscillator wave function in the lowest energy level corresponding to the frequencies $`\omega _\pm `$ and $`\mathrm{tan}\stackrel{~}{\theta }=\frac{\omega _c^2}{\omega _+^2\omega _c^2}\mathrm{tan}\theta `$, with the cyclotron frequency $`\omega _c=eB/m_bc`$. The frequencies $`\omega _\pm `$ are given byyang2

$$\omega _\pm ^2=\frac{1}{2}(\mathrm{\Omega }^2+\frac{\omega _c^2}{\mathrm{cos}^2\theta })\pm \sqrt{\frac{1}{4}(\mathrm{\Omega }^2\frac{\omega _c^2}{\mathrm{cos}^2\theta })^2+\mathrm{\Omega }^2\omega _c^2\mathrm{tan}^2\theta }.$$
(9)
Obviously, the existence of an in-plane field deforms the electron wave function. However, this wave function deformation does not qualitatively change the electron hopping mechanism at a given temperature. The major effect of the tilted field would be on the variation of $`T_0`$. Fig. 1 illustrates the cohesive energy for two typical configurations of the crystal orientation with respective to the in-plane field: the $`[100]`$ or the $`[110]`$ direction parallel to the in-plane field. The energy is always lower for the case of the $`[110]`$ direction parallel to the in-plane field. The applied in-plane field lowers the cohesive energy of the electronic crystal and forces the domains align to the same direction. Thus, the role of the in-plane field is to integrate the domains into larger ones. In this way, the in-plane field causes some of impurities to be irrelevant and therefore reduces the effective impurity density. In determining $`T_0`$ from (7), only the relevant impurities should be counted in. Hence, one can replace $`n_I`$ by an effective impurity density $`n_I(B_{})`$. From eq. (7), we see that $`T_0`$ is sensitive to $`n_I(B_{})`$. Therefore, it is possible to observe the influence of the in-plane magnetic field on the characteristic temperature $`T_0`$ under proper parameters as the tilting angle varies.

In conclusion, the temperature dependence of conductivity is explored for the Wigner crystal in 2-dimensional electrons under a strong magnetic field. We argued that there are domains of electronic crystal in a realistic sample and predicted that the temperature dependence of the transport behavior may be different in different temperature regimes. We found that the conductivity experiences a crossover from the fixed range hopping to the variable range hopping mechanism as the temperature varied. The crossover temperature $`T_0`$ is reachable under present experimental technique. Finally, the possible effect on the crossover region by the in-plane magnetic field is discussed.

This work was supported in part by the NSF of China.

Figure Captions

Figure 1 The cohesive energy of the WC versus the tilting angle for two configurations of the crystal orientation. The real line: $`[110]`$ parallel to the in-plane field (as shown in the inset); The broken line: $`[100]`$ parallel to the in-plane field.


format: text

text:
# 1 Results of DFA2 for the daily maximum temperature records of 18 inner continental stations (∘) and two coastline stations (×): (a) 1 Urumchi, 2 Tomsk, 3 Atbasar, 4 Chita, 5 Olekminsk, 6 Horog, 7 Swerdlowsk, 8 Surgut, 9 Jakutsk, 10 Aleksandrovsk; (b) 1 Huron, 2 Academy, 3 Cheyenne, 4 Gothenburg, 5 Gunnison, 6 Spokane, 7 Winnemucca, 8 Pendleton, 9 Tusc¿on, 10 New York. The estimated minimum distance to the oceans is written in parentheses. The scale of 𝐹(𝑠) is arbitrary. The straight lines in the curves represent the best fits between 𝑠=150 and 2500. For each curve, the variance of the slope is about 0.01. The line at the bottom has a slope of 0.15, corresponding to 𝛼=0.65.

Comment on ”Scaling of Atmosphere and Ocean Temperature Correlations in Observations and Climate Models”

In a recent Letter , Fraedrich and Blender (FB) studied the scaling of atmosphere and ocean temperature. They analyzed the fluctuation functions $`F(s)s^\alpha `$ of monthly temperature records (mostly from grid data) by using the detrended fluctuation analysis (DFA2) and claim that the scaling exponent $`\alpha `$ over the inner continents is equal to $`0.5`$, being characteristic of uncorrelated random sequences. Here we show that this statement is (i) not supported by their own analysis and (ii) disagrees with the analysis of the daily observational data from which the grid monthly data have been derived. We conclude that also for the inner continents, the exponent is between $`0.6`$ and $`0.7`$, similar as for the coastline-stations.

(i) Figure 1a in  shows the representative results of FB for $`F(s)`$ for the inner continental site of Krasnojarsk (observational, grid and model data). Close inspection of the curves shows that none of the observational data approaches the exponent 0.5, but rather yield exponents close to 0.6 or above.

(ii) Figure 1 shows representative results of $`F(s)`$ (DFA2) obtained from daily observational records, for 18 inner continental sites in North America and Asia, where the stations are between 400 and 2600 km away from the ocean. For comparison, we also show the results for two coastline stations. The maximum s-value in each curve is below one quarter of the length of the corresponding record. To facilitate the evaluation of the data, we have divided $`F(s)`$ by $`s^{0.5}`$. A plateau now indicates loss of correlations. One can see clearly that none of the curves approaches a plateau, i. e. an exponent $`\alpha =0.5`$ is never seen. All asymptotic slopes have values above $`0.6`$. The mean value of $`\alpha `$, averaged over all innercontinental stations, is 0.65 $`\pm 0.04`$. There is no remarkable difference between coastline and inner continental stations. Figure 2 shows the dependence of $`\alpha `$ on the shortest distance to the oceans. There is no tendency towards a lower exponent at larger distances.

Finally, we like to comment on the claim that $`\alpha 1`$ for sea surface temperatures. As has been shown in , there is a remarkable crossover at about 1 year in the sea surface temperatures: At small scales, the exponent is significantly larger than 1, $`\alpha 1.3`$, while at large time scales $`\alpha `$ is between $`0.65`$ and $`0.95`$, with the average at $`0.8`$.

This work has been supported by the Deutsche Forschungsgemeinschaft and the Israeli Science Foundation.

A. Bunde1, J. F. Eichner1, S. Havlin2, E. Koscielny-Bunde1,3, H. J. Schellnhuber3,4, and D. Vjushin2
1 Institut für Theoretische Physik III, Universität Giessen, D-35392 Giessen, Germany
2 Minerva Center and Department of Physics, Bar-Ilan University, Ramat-Gan 52900, Israel
3 Potsdam Institute for Climate Impact Research, D-14412 Potsdam, Germany
4 Tyndall Center for Climate Change Research, Norwich, NR4 7JD, England

PACS numbers: 92.10.Fj 89,75,DA 92.70.Gt


text:
# Angular Momentum Transport by Internal Gravity Waves

## 1 Clues to angular momentum transport in low-mass stars

In many locations of the Hertzsprung-Russell diagram, stars exhibit signatures of processes that require challenging modeling beyond the standard stellar theory111By standard we refer to the modeling of non-rotating, non-magnetic stars, in which convection is the only large-scale mixing considered.. In this context, rotation has become a major ingredient of modern models, especially when abundance anomalies have to be accounted for. In order to correctly describe the effects of rotation on stellar structure and evolution as well as on their “byproducts” like surface abundances or chemical yields, special emphasis has to be put on the evolution of the angular momentum distribution; this indeed is the feature that determines the extent and magnitude of rotation-induced mixing in stars.

In recent theoretical developments, the internal rotation law evolves as a result of contraction, expansion, meridional circulation, shear turbulence and mass loss; mixing of chemicals is directly linked to the rotation profile (see e.g Talon 2003). The most sophisticated treatment of these hydrodynamical processes is based mainly on the work by Zahn (1992), Maeder (1995), Talon & Zahn (1997) and Maeder & Zahn (1998). It rests on only one assumption which is that the turbulence sustained by the shear is highly anisotropic. Two free parameters describe the magnitude of the horizontal shear and the erosion of the restoring force due to both the thermal and the mean molecular weight stratifications.

Such a self-consistent treatment has been successfully applied in various parts of the HR diagram. The strength of the theory is that the use of the same free parameters give very satisfactory results for stars over a large range of masses and evolutionary phases. For example, it explains the observed He and N enrichment in main sequence O-type and early B-type stars, in OB supergiants (e.g., Gies & Lambert 1992, Lennon et al. 1991, Herrero et al. 1999) and in A-type supergiants in the SMC (Venn 1999), as well as the B depletion in main sequence B-type stars (e.g., Venn et al. 1996); it also helps reproducing the number ratio of blue to red supergiants in the SMC (e.g., Meylan & Maeder 1982) and the observed WR/O ratios at solar metallicity. We refer to Maeder & Meynet (2000) for a more detailed description of the effects of rotation in massive stars and for relevant references.

In low-mass stars, the strongest observational constraints come from the so-called Li dip and from helioseismology. The Li dip is a characteristic feature of lithium abundances which is seen both in the field and in open clusters (e.g., Wallerstein et al. 1965, Boesgaard & Tripicco 1986, Balachandran 1995). It corresponds to a narrow region in effective temperature (around $`T_{\mathrm{eff}}`$6650 K) where the surface Li abundances are reduced by up to $``$2.5 dex. In the same region, low-mass stars are beginning to be spun down efficiently early on the main sequence (Boesgaard 1987), most probably via magnetic torquing. In modern rotating stellar models in which the transport of angular momentum and of chemicals by meridional circulation and shear turbulence is self-consistently taken into account, this leads to important lithium (and beryllium) destruction at the correct $`T_{\mathrm{eff}}`$ (Talon & Charbonnel 1998, hereafter TC98; Palacios et al. 2003, hereafter PTCF03). In addition rotation induced mixing inhibits the atomic diffusion in a way that explains the constancy of the CNO abundances within the Li dip (e.g., Varenne & Monier 1998, Takeda et al. 1998). Last but not least, these models do explain the evolution of Li surface abundances in evolved stars that originate from the hot side of the dip. It is worth recalling that all these results for low-mass stars on the blue side of the dip are obtained with the same values for the two above-mentioned parameters as those used in massive stellar models.

For still cooler main sequence stars however, the magnetic torque strengthens as the stellar convective envelope grows. If we assume here that all the momentum transport is assured by the wind-driven meridional circulation as it is on the blue side of the Li-dip, then too much lithium burning is obtained below $`T_{\mathrm{eff}}6550`$ K. Thus the rise of the Li abundance on the red side of the dip can be interpreted as the signature of another mechanism that efficiently transports angular momentum (but not chemicals) in lower mass stars (TC98). As a result, in these objects the magnitude of both meridional circulation and shear turbulence is reduced, as well as the Li depletion due to rotational mixing.

Such a mechanism is also required to shape the Sun’s flat rotation profile (Chaboyer et al. 1995, Matias & Zahn 1998). At the solar age, models relying solely on turbulence and meridional circulation for momentum transport still predict large angular velocity gradients which are not present in the Sun.

Only two mechanisms have been proposed so far to transport angular momentum in addition to the classical hydrodynamical processes in order to enforce the flat rotation profile measured by helioseismology : magnetic fields (Charbonneau & MacGregor 1993, Barnes et al. 1999) and gravity waves (Schatzman 1993, Kumar & Quataert 1997, Zahn, Talon & Matias 1997). As explained in TC98, in order to obtain a consistent picture of rotational mixing in all types of stars, the correct mechanism must become efficient in the center of the Li dip, namely at an effective temperature of $``$6650 K.

While the original propositions for momentum extraction by gravity waves have received much criticism (Gough & McIntyre 1998, Ringot 1998), recent calculations performed by Talon, Kumar & Zahn (2002, hereafter TKZ02) show that, through differential filtering, gravity waves are indeed able to extract momentum from the solar interior. Here we explore the same physics in the stellar mass range around the Li dip in order to verify whether it has the proper mass dependence as described by TC98 (see Fig. 1). Our goal is to present a coherent understanding of rotational mixing in all stars of the main sequence. We will firstly recall the main descriptions that have been made for wave excitation (§ 2) and explain how thermal dissipation of these waves leads to the formation of a filtering shear layer (§ 3). Results for momentum extraction by waves in stars of various masses are then presented (§ 4) and discussed (§ 5).

## 2 Turbulent wave excitation

### 2.1 Two excitation mechanisms

The generation of gravity waves depends on the structure of the stellar convective envelope. On the main sequence, this region is rather large in solar mass stars ($``$30$`\%`$ in radius), while it shrinks as the stellar mass rises, representing only $``$1$`\%`$ in a $`1.45M_{}`$ star lying inside the Li dip. The characteristic convective frequency varies from $`\nu _c`$0.75 $`\mu `$Hz to $`\nu _c`$600 $`\mu `$Hz while the spherical harmonic number $`\mathrm{}`$ corresponding to the largest convective scale varies from $`\mathrm{}_c`$60 to $`\mathrm{}_c`$1000 in the same mass range. These variations are not smooth, rapid changes of convective properties occurring when $`M_{}\mathrm{}>\mathrm{\hspace{0.17em}1.3}M_{}`$. We wish to investigate how the wave spectrum produced by these convection zones varies as a function of stellar mass.

The exact properties of wave spectra remain somewhat uncertain. Excitation can be related to internal stresses that correlate with the mode’s eigenfunction (Goldreich & Kumar 1990, Balmforth 1992, Goldreich et al. 1994). Gravity waves may also be excited by penetration below a convection zone, as observed in laboratory experiments (Townsend 1958) and in numerical simulations (Hurlburt et al. 1986, 1994, Nordlund et al. 1996, Kiraga et al. 2000). In a real star, both sources would contribute to wave generation and they are thus additive.

In this paper, we examine both sources, with emphasis not on their absolute values, but rather on their variation with stellar mass. Our computations are performed for stars between $`1.1`$ and $`1.45M_{}`$ with $`Z=0.02`$ (see Table 1). We use the same code and input physics (eos, opacities, nuclear reactions) as in PTCF03. Wave characteristics are computed on the basis of ZAMS stellar structures. Note that in these objects the depth and structure of the convective envelope does not vary significantly over a main sequence lifetime (at least up to the age of the Hyades).

### 2.2 Excitation by the Reynolds stresses

First, we consider excitation by the Reynolds stresses as described by Goldreich et al. (1994). This description, which has also been used by Kumar, Talon & Zahn (1999, hereafter KTZ99), uses a free parameter that has been calibrated on the solar p-mode spectrum.

The energy flux per unit frequency due to Reynolds stresses $`_{E}^{}{}_{}{}^{R}`$ is then

$`_{E}^{}{}_{}{}^{R}(\mathrm{},\omega )`$ $`=`$ $`{\displaystyle \frac{\omega ^2}{4\pi }}{\displaystyle _{r_c}^R}𝑑r{\displaystyle \frac{\rho ^2}{r^2}}\left[\left({\displaystyle \frac{\xi _r}{r}}\right)^2+\mathrm{}(\mathrm{}+1)\left({\displaystyle \frac{\xi _h}{r}}\right)^2\right]`$ (1)
$`\times \mathrm{exp}\left[h_\omega ^2\mathrm{}(\mathrm{}+1)/2r^2\right]{\displaystyle \frac{v^3L^4}{1+(\omega \tau _L)^{15/2}}},`$
where $`\xi _r`$ and $`[\mathrm{}(\mathrm{}+1)]^{1/2}\xi _h`$ are the radial and horizontal displacement wave-functions which are normalized to unit energy flux just below the convection zone, $`v`$ is the mixing length convective velocity, $`L`$ is the radial size of an energy bearing turbulent eddy, $`\tau _LL/v`$ is the characteristic convective time, and $`h_\omega `$ is the radial size of the largest eddy at $`r`$ with characteristic frequency of $`\omega `$ or greater ($`h_\omega =L\mathrm{min}\{1,(2\omega \tau _L)^{3/2}\}`$). The gravity waves are evanescent in the convection zone, the region where they are excited. Their wave-functions $`\xi _r`$ and $`\xi _h`$ are thus proportional to $`k_r^{1/2}\mathrm{exp}\left(idrk_r\right)`$. The above equation was derived under the assumption that the turbulence spectrum is Kolmogorov. The momentum flux per unit frequency $`_{J}^{}{}_{}{}^{R}`$ is then

$$_{J}^{}{}_{}{}^{R}(m,\mathrm{},\omega )=\frac{m}{\omega }_{E}^{}{}_{}{}^{R}(\mathrm{},\omega ).$$
(2)
The momentum spectra corresponding to the various stellar masses are shown in Fig. 2 where $`_J=4\pi r_{\mathrm{cz}}^2_J`$ is the momentum luminosity; note that some of the lower frequency waves have zero magnitude. This is due to the fact that the corresponding damping is too large to permit them to be waves.

Spectrum characteristics evolve together with the structure of the convection zone which depends on the stellar mass. Low frequency waves disappear when the stellar mass increases. This is related to the fact that the surface convection zone then becomes thinner, leading to a larger thermal diffusivity just below it. This in turn implies stronger damping (see Eq. 3.1) and the disappearance of low frequency waves. Furthermore, as stellar mass rises up to about $`1.35M_{}`$ (i.e., $`T_{\mathrm{eff}}`$6400 K), so does the flux associated with a given frequency. This is related to the rise of the luminosity with mass, and thus, of the energy in convective motions. For even more massive stars, the flux associated with a given frequency remains about constant, but most low frequency waves are damped too rapidly to exist.

Let us recall that these results are obtained for stellar structures on the ZAMS, and that we have found no significant dependence of the wave spectrum with age on the main sequence; it can thus be considered as constant for the duration of the main sequence.

In this description, waves with frequencies up to $`N_c`$ (the Brunt-Väisälä frequency at the base of the convective envelope) can be excited. However, the high-frequency waves will have little impact on momentum evolution; indeed, as frequencies increase, filtering by the shear layer (see §3) becomes less efficient and less differential. Furthermore, damping is so small that it leads to only a small amount of momentum redistribution. Finally, as wave excitation diminishes with frequency, they carry much less momentum overall. These high-frequency waves will form standing waves, the g-modes of helio- and astero-seismology.

### 2.3 Excitation by convective penetration

Greater uncertainties exist on wave generation via convective penetration. The only theoretical estimates that exist have been made by Press (1981) and García López & Spruit (1991), and later used by Zahn et al. (1997). It is the formulation we shall adopt here.

This description is based on the matching of wave pressure fluctuations with those of turbulent convection. Furthermore, it takes into account the combination of turbulent eddies of a given size to produce larger fluctuations. The range of horizontal scales available is thus

$$0<\mathrm{}<\mathrm{}_u\mathrm{with}\mathrm{}_u=\mathrm{}_c\left(\frac{\omega }{\omega _c}\right)^{3/2}$$
(3)
where $`\mathrm{}_c`$ is the spherical harmonic number associated with the largest convective scale and $`\omega _c`$ is the corresponding convective frequency. Turbulence is assumed to follow a Kolmogorov spectrum, all frequencies with $`\omega \omega _c`$ thus being available. The associated energy flux $`_{E}^{}{}_{}{}^{P}`$ is

$$_{E}^{}{}_{}{}^{P}(\mathrm{},\omega )=\frac{\rho v_c^3}{2}\frac{\omega _c^3}{\omega ^2}\frac{\mathrm{}}{\mathrm{}_c}\frac{1}{N}\left(1\frac{\omega ^2}{N^2}\right)^2$$
(4)
and the corresponding angular momentum flux is still given by Eq. (2).

### 2.4 A word of caution

Let us stress that the fluxes calculated here are somewhat uncertain, especially in the second case. Indeed, it is well known that the “structure” of convection, which contains e.g. plumes that travel down across the whole convection zone, is not well reproduced by the mixing length model even though the convective velocities are more or less correct (e.g. Hurlburt et al. 1986). However, their differential properties should have the proper dependence, and this will be our focus in the rest of this paper.

## 3 Formation of a shear layer

### 3.1 Differential filtering of the waves by a shear layer

As described by several authors (Gough & McIntyre 1998, Ringot 1998, KTZ99, Kim & MacGregor 2001), the most obvious feature of momentum transport by gravity waves is the formation of a double peaked shear layer just below the convection zone (see Fig. 3).

Gough & McIntyre (1998) and Ringot (1998) argued that such a layer would prevent waves from propagating beyond. Indeed, local radiative and viscous damping is greatly dependent on the local frequency and increases when that frequency diminishes (this property leads to the formation of the double peaked shear layer). The local momentum luminosity integrated over the whole spectrum writes

$$_J(r)=\underset{\sigma ,\mathrm{},m}{}_{J}^{}{}_{\mathrm{},m}{}^{}\left(r_{\mathrm{zc}}\right)\mathrm{exp}\left[\tau (r,\sigma ,\mathrm{})\right].$$
(5)
The local amplitude $`\mathrm{exp}\left[\tau (r,\sigma ,\mathrm{})\right]`$ depends on the integrated damping due to thermal diffusion $`K_T`$ and (turbulent) viscosity $`\nu _t`$

$`\tau (r,\sigma ,\mathrm{})=`$
$`[\mathrm{}(\mathrm{}+1)]^{\frac{3}{2}}{\displaystyle _r^{r_c}}\left(K_T+\nu _t\right){\displaystyle \frac{NN_T^2}{\sigma ^4}}\left({\displaystyle \frac{N^2}{N^2\sigma ^2}}\right)^{\frac{1}{2}}{\displaystyle \frac{\mathrm{d}r}{r^3}}`$
where $`N^2=N_T^2+N_\mu ^2`$ is the Brunt-Väisälä frequency, $`N_T^2`$ is its thermal part and $`N_\mu ^2`$ is due to the mean molecular weight stratification (Zahn et al. 1997).

The prograde peak222The prograde peak corresponds to the layer that rotates more rapidly than the convection zone, while the retrograde peak designates the layer that rotates more slowly. thus filters prograde waves, while the retrograde peak filters retrograde waves. If the peaks were infinite in height, not a single wave could travel through it. This is however not the case, as the magnitude of the shear layer is self regulated by shear turbulence333Baroclinic instabilities could also set in, and would have the same effect..

As described by TKZ02, if no initial differential rotation is present, the average wave momentum luminosity that traverses the shear layer is null, and there is no net effect on the interior. However, if differential rotation is initially present, as in the case of low mass stars that are braked by a magnetic torque, the average magnitude of the two peaks is not equal; in the case of interest (i.e. with the surface convection zone rotating more slowly than the interior), the prograde peak is always larger, leading to differential filtering that favors the penetration of retrograde waves. These waves are then damped in the whole radiative region that they spin down by depositing negative momentum.

### 3.2 Dynamics of the shear layer

To study the effect of gravity waves in low mass stars, one must thus be able not only to evaluate the global momentum luminosity in waves, but also the dynamics of the shear layer and the resulting differential filtering. This is addressed here by performing numerical simulations of angular momentum evolution. We study a few cycles and evaluate the resulting net momentum luminosity for various masses and various differential rotations.

Angular momentum evolves under the action of gravity waves and shear turbulence according to

$$\rho \frac{\mathrm{d}}{\mathrm{d}t}\left[r^2\mathrm{\Omega }\right]=\frac{1}{r^2}\frac{}{r}\left[\rho \nu _tr^4\frac{\mathrm{\Omega }}{r}\right]\frac{3}{8\pi }\frac{1}{r^2}\frac{}{r}_J(r)$$
(7)
where $`\rho `$ is the density, $`\mathrm{\Omega }`$ the angular velocity and $`\nu _t`$ the turbulent viscosity. The local momentum luminosity $`_J(r)`$ is integrated over the whole wave spectrum. Equation (7) is solved throughout the radiation zone. The upper boundary condition expresses conservation of momentum of the star as a whole, with the convection zone rotating as a solid body. On their way to the core, waves are reflected when their local frequency $`\sigma `$ equals the Brunt-Väisälä frequency

$$\sigma (r,m)\omega m\left[\mathrm{\Omega }(r)\mathrm{\Omega }_c\right]=N.$$
(8)
Waves then deposit more momentum as they travel back to the surface convection zone. The high-frequency waves of lowest degree considered in this study can be reflected up to 5000 times before they are completely damped; these contribute little to momentum redistribution444As the frequency rises, damping diminishes and standing waves may form..

Figure 3 shows the oscillation cycle of the shear layer for our $`1.2M_{}`$ star. Each model (i.e., each stellar mass) has a different characteristic time-scale for this oscillation which depends mainly on the net wave flux in the high degree, low frequency wave. It is thus shortest for masses around $`1.21.3M_{}`$.555Let us note that the very different time-scale found in TKZ02 for the Sun is related not only to the smaller flux but also to the reduced number of frequencies and degrees used in order to be able to perform the solar calculation over the required $`10^7`$ time-steps.

## 4 Momentum extraction

### 4.1 Fluxes

In all our stellar models, we follow the detailed time evolution of the shear layer according to Eq. (7) over several cycles, for two values of the initial differential rotation (a very small differential rotation of $`\delta \mathrm{\Omega }=0.00005\mu \mathrm{Hz}`$ over $`0.05R_{}`$, and a larger differential rotation $`\delta \mathrm{\Omega }=0.05\mu \mathrm{Hz}`$ over $`0.05R_{}`$). We calculate the net luminosity below the filtering shear layer, at $`r=R_{\mathrm{cz}}0.03R_{}`$. Resulting instantaneous and mean luminosities are detailed for the $`1.2M_{}`$ model in Fig. 4. The change of sign of the instantaneous luminosity illustrates the varying asymmetry of the double peak shear layer. In the case of very small initial differential rotation, the curve is symmetric and the mean luminosity is always close to zero. On the other hand, if the initial differential rotation is large, retrograde waves are favored, as illustrated by the negative mean luminosity. It is this net mean luminosity below the shear layer that gives rise to a net momentum extraction in the interior.

In order to obtain differential filtering, it is important to include several different modes;

* High-$`\mathrm{}`$ modes, that damp close to the convection zone (see Eq. 3.1);
* Low-$`\mathrm{}`$ modes, that penetrate deeply.

Furthermore, it is the lowest frequency waves (which undergo the largest amount of differential damping) that give rise to momentum extraction. While wave frequencies range from $`0`$ to $`N`$, for each considered mass, there is a maximum frequency for waves that contribute to differential filtering. That maximum frequency has been determined for each individual mass (see Table 1). In the most massive models considered here, which correspond to stars inside the dip, this maximum frequency is much higher than in lower mass stars. Indeed, as can be seen from the wave spectrum (Fig. 2), in those more massive models, high-$`\mathrm{}`$ modes of low frequency experience too much damping to travel beyond a single vertical wavelength in the vertical. One must thus include higher wave frequencies in order to capture high-$`\mathrm{}`$ modes that produce the oscillating double peaked shear layer.

The García López & Spruit model for wave excitation by convective penetration leads to the generation of larger frequency waves, leading to a limited amount of extraction. This is especially true for the more massive stars presented in this study and is thus less relevant to our discussion here666Let us remind that both generating sources are actually additive.. For that reason, that model has been ignored in detailed calculations.

The time evolution of the resulting filtering in all stellar masses considered is shown in Fig. 5. The abscissa represents the time iteration, each mass having its own characteristic time-scale for oscillation. The results are thus shown in terms of this characteristic time, and not in a real physical time sequence. The net luminosity slightly increases with the stellar mass up to $`1.2M_{}`$, and then dramatically decreases as one moves to stellar masses inside the Li dip. For the largest mass considered here ($`1.44`$ and $`1.45M_{}`$), the process of differential filtering does not lead to momentum extraction anymore (the net luminosity is actually positive although very small for these stars on the blue side of the dip777This reverse effect is probably linked to the fact that thermal damping dominates just below the convection zone. Then, the low-$`\mathrm{}`$, high frequency waves that are included and that travel back and forth over 100 times, are damped mostly there. However, retrograde waves that have a larger frequency in the interior are less damped there than the prograde waves, leading to an effect opposite to what is observed in other stars.).

Figure 6 represents the average luminosity after 500 time steps as a function of effective temperature. On the same figure are shown the limits established by TC98 in terms of momentum extraction required in order to consistently understand the Li dip in terms of rotational mixing.

### 4.2 Coriolis force

In order to fully describe wave transport in rotating stars, we must also take into account the influence of the Coriolis force on waves. Its first order effect is to change the horizontal structure of the modes, confining them closer to the equator. The maximum co-latitude of propagation $`\theta _{\mathrm{crit}}`$ is given by the condition

$$\sigma ^2=4\mathrm{\Omega }^2\mathrm{cos}^2\theta _{\mathrm{crit}}$$
(9)
(KTZ99). For a given rotating velocity, this condition implies that only a fraction

$$𝒮=_{\theta _{\mathrm{crit}}}^{\pi /2}\mathrm{sin}\theta \mathrm{d}\theta $$
(10)
of the surface will support waves and the corresponding luminosity should be diminished accordingly. Table 2 gives the critical angle for propagation and the corresponding efficiency for a typical frequency of $`1\mu `$Hz. The rotation velocities used here correspond, for each stellar mass, to the mean velocity observed at the age of the Hyades. Let us note that the fraction of the surface $`𝒮`$ varies linearly with $`\nu `$, and that the whole surface will support waves when $`\nu 2\nu _{\mathrm{rot}}`$.

Exact calculations can be made by reducing the wave luminosity corresponding to each frequency by the proper coefficient. This has been done for 3 models ($`1.2`$, $`1.3`$ and $`1.4M_{}`$). The amplitude of low frequency waves (that dominate the formation of the shear layer) being reduced, including the Coriolis force leads to a somewhat smaller double peaked shear layer, that oscillates on a slightly longer period. Furthermore, differential filtering corresponding to momentum extraction is also dominated by rather low frequency waves. This leads to a smaller efficiency of momentum extraction by waves (see Fig. 6). However, the overall mass dependence is not affected.

## 5 Discussion

Several abundance anomalies are best explained in terms of large-scale mixing. As a natural physical parameter (beyond mass, chemical composition and age), rotation is an important factor in understanding various features which are not predicted from classical stellar models (see Maeder & Meynet 2000 for details). Furthermore, if one is to built a consistent model of rotational mixing, it has to apply to all stars.

Recent improvements of the description of the Eddington-Sweet meridional circulation by Zahn (1992) and Maeder & Zahn (1998) have been applied to several kinds of stars. These studies show that, with the same free parameters, it is possible to explain abundance anomalies in B (Talon et al. 1997) and O stars (Maeder & Meynet 2000) as well as the blue side of the Li dip (TC98, PTCF03) and Li anomalies in sub-giants (Charbonnel & Talon 1999).

However, in order to understand the helioseismic rotation profile, some extra mechanism for momentum transport is required in low mass stars (Matias & Zahn 1998). TC98 argued that, for rotational mixing to be responsible for the Li dip, this mechanism has to become efficient at an effective temperature around $`65006600`$K, i.e., in a region where the stellar convective envelope becomes substantial. If the same description is to be applied to all stars, there must be some physical reason for the lack of efficiency of this mechanism in hotter stars. The goal of this paper was to establish whether gravity waves do have the proper effective temperature dependence.

The calculations we performed here show that the built up of a double peaked shear layer allowing differential filtering of the waves is possible in low-mass main sequence stars up to an effective temperature of $`6700`$K. The net momentum luminosity below the shear layer increases with stellar mass up to $`1.2M_{}`$ (i.e., $`T_{\mathrm{eff}}59005950`$K); it then drops in hotter stars, even becoming slightly positive for stars more massive than $`1.44M_{}`$ (i.e., $`T_{\mathrm{eff}}>6700`$K). While the quantitative results shown here are based on the excitation model by Goldreich et al. (1994), the García López & Spruit (1991) model would lead to an even stronger decrease in efficiency. Indeed, as can be seen in Fig. 2, in the more massive models, there are simply no low frequency waves, implying an even more abrupt decrease in the efficiency of differential filtering.

Let us discuss a final issue. As can be seen from Fig. 1, the dispersion in lithium abundances is somewhat larger on the cold side of the dip ($`T_{\mathrm{eff}}60006500`$K) than in cooler stars. A very straightworward explanation of this arises within our framework. As we have shown indeed, for a given stellar mass a higher rotation velocity implies a stronger reduction of the net momentum luminosity in waves due to the Coriolis force. As a consequence, the efficiency of angular momentum extraction by the waves is somewhat lower within faster rotators, leading to a larger efficiency of the hydrodynamical processes, and thus to stronger Li destruction. In this context, the larger lithium dispersion measured in stars laying immediatly on the right of the dip can be interpreted as a consequence of the larger dispersion in rotational velocities that these stars do exhibit in comparison with the less massive ones (Gaigé 1993). Detailed calculations would have to be performed to check whether this effect could be large enough to explain the observed dispersion.

All our results clearly indicate that momentum transport by gravity waves has the proper mass dependence to be the required process in low-mass main sequence stars on the cold side of the Li dip.

Note that a similar mass dependence is not expected if momentum transport is dominated by magnetic fields. There is another test which could permit to differentiate the magnitude of transport by magnetic field compared to transport by waves which applies to the Sun. Certain early helioseismic inversions (Elsworth et al. 1995, Corbard et al. 1997) suggested that the central region of the Sun (below $`0.4R_{}`$) could rotate somewhat slower than the rest of the radiative region. These results seem to be confirmed by recent SOHO data (Couvidat et al. 2003). TKZ02 explained how such a feature could be produced by gravity waves. Both arguments favor momentum transport by gravity waves with respect to magnetic field.

We are now in a position to present a coherent picture of rotational mixing in main sequence stars of all masses:
On the hot side of the Li dip and in more massive stars, the transport of angular momentum and of chemicals by meridional circulation and shear instabilities do explain the Li as well as the He and CNO patterns. In lower mass stars, gravity waves dominate the transport of angular momentum, thereby reducing the magnitude of meridional circulation and shears and shaping the Li pattern on the cold side of the dip. Within this framework, we thus predict that Pop I main sequence stars with initial masses lower than $`1.4M_{}`$ must be quasi-solid body rotators, as the Sun is.

###### Acknowledgements.

We acknowledge financial support from the French Programme National de Physique Stellaire (PNPS). We thank the Réseau québécois de calcul de haute performance (RQCHP) and the Centre informatique national de l’enseignement supérieur (CINES) for useful computational resources.


text:
# Modern nuclear force predictions for the neutron-deuteron scattering lengths

## I Introduction

It has been observed a long time ago that the nd scattering length for total three-nucleon spin S=1/2 ($`{}_{}{}^{2}a_{nd}^{}`$) is correlated to the $`{}_{}{}^{3}H`$ binding energy ($`E_{{}_{}{}^{3}H}`$). This correlation is known as the Phillips line . Indeed, calculations years later based on simplistic or more realistic NN model forces (see ) yielded quite a few results for the $`{}_{}{}^{3}H`$ binding energy and the $`{}_{}{}^{2}a_{nd}^{}`$ scattering length, which lie on or rather close to a line in the two-dimensional plane spanned by $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$. Also 3N forces of the $`2\pi `$-exchange type have been added. In  it was found that this line passes well through the experimental point.

In recent years chiral perturbation theory and effective theories have been applied to nuclear physics. In the pionless formulation , adequate for extreme low energy phenomena, it has been shown that $`{}_{}{}^{3}H`$ can be energetically stabilized only if a 3N contact force is introduced (see however ). In the two lowest orders of that framework there is just one parameter connected to that 3N force. Thus both quantities, $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$, depend on that one parameter and are therefore correlated though the line does not hit the experimental point. In higher orders additional parameters show up and the correlation is broken, which makes the two quantities independent. The same observation was made in an approach based on chiral perturbation theory  which includes explicitely the pion degrees of freedom. In the next-to-next-to-leading-order (NNLO), 3N forces occur the first time and they depend on two parameters. This makes $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$ independent and the Phillips line correlation is broken. In fact in that framework the two experimental values are used to fix those two parameters of the 3N force . Thus we find it interesting to ask, whether the conventional and high precision NN forces AV18 , CD Bonn 2000 , Nijm I, Nijm II, and Nijm 93  alone or in combination with the two most popular 3N force models , Urbana IX  and TM99  (an updated Tucson-Melbourne $`2\pi `$-exchange 3NF  modified in view of chiral symmetry) lead to a strict correlation between $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$ or whether that Phillips line correlation is also broken. Further we ask whether the NN and 3NF combinations adjusted to $`E_{{}_{}{}^{3}H}`$ (or may be only one of them) also describe $`{}_{}{}^{2}a_{nd}^{}`$. One more reason to confront $`{}_{}{}^{2}a_{nd}^{}`$ to state-of-the-art calculations is the recent appearance of a precision neutron interferometric measurements of the nd coherent scattering length ($`b_{nd}`$.

The coherent scattering length $`b_{nd}`$ depends in addition to $`{}_{}{}^{2}a_{nd}^{}`$ also on the second s-wave scattering length for the state of total 3N spin S=3/2 , $`{}_{}{}^{4}a_{nd}^{}`$ . Because of the Pauli principle this quantity is supposed not to be sensitive to short range details of nuclear forces. We want to investigate also that quantity in the light of modern nuclear forces.

Then we would like to add two more investigations. Charge-symmetry breaking (CSB) in the strong NN forces is mostly pronounced in the states $`{}_{}{}^{1}S_{0}^{}`$, where the scattering length for the neutron-neutron (nn), $`a_{nn}`$, and proton-proton (pp), $`a_{pp}`$, systems are different. However the value for $`a_{nn}`$ is still under debate . Therefore we would like to present results, where the nn forces are replaced by the (strong) pp forces. This will provide some insight into the magnitudes of the shifts in $`{}_{}{}^{2}a_{nd}^{}`$ caused by small changes in $`a_{nn}`$.

The other investigation is due to effects on the scattering lengths and $`E_{3H}`$ caused by electromagnetic interactions, mostly due to magnetic moment interactions (MMI). MMI is a relativistic effect and including only that specific force is of course inconsistent, since other relativistic effects are not taken into account (see for instance  and references therein). But it is interesting to see its separate effect on $`{}_{}{}^{2}a_{nd}^{}`$ and $`{}_{}{}^{4}a_{nd}^{}`$ (its effect on the binding energy of $`{}_{}{}^{3}H`$ and $`{}_{}{}^{3}H`$e is known and older results have been reconfirmed recently ). Here we hit some ”defects” in current NN force models. The NN potentials CD Bonn 2000, Nijm I, II, and Nijm 93 are fitted directly to the NN data without taking electromagnetic interactions (EMI) into account (of course the point Coulomb force in case of the pp system has been included). Therefore those strong forces include the effects of the MMI’s (and further electromagnetic corrections). If we want to see then the separate effects of the EMI’s we have to subtract from those model NN forces the EMI’s and compare to results without that subtraction. In case of AV18 the strong force plus EMI’s have been fitted to the data. Thus the force free of EMI’s is just the strong AV18 force alone. In this manner by comparing results with and without electromagnetic interactions one can see their effects.

The paper is organised as follows. In Section II the theoretical formulation is briefly outlined. The results are given in Section III, and we end with a summary and an outlook in Section IV. More technical details are deferred to an appendix.

## II Formulation

We use the Faddeev scheme. Including a three-nucleon force a convenient basic formulation for part of the $`ndn+n+p`$ break-up amplitude $`T`$  is the integral equation

$`T`$ $`=`$ $`tP\varphi +(1+tG_0)V_4^{(1)}(1+P)\varphi +tPG_0T+(1+tG_0)V_4^{(1)}(1+P)G_0T.`$ (1)
The driving term contains the NN $`t`$-operator $`t`$, permutation operators $`P`$, the free 3N propagator $`G_0`$ and part of the 3N force, $`V_4^{(1)}`$. Any 3N force can be split into three pieces, where for instance the first piece is symmetrical under exchange of particles 2 and 3 , the second under 3-1 exchange etc. Thus the quantity $`V_4^{(1)}`$ is the part symmetrical under 2-3 exchange like the $`t`$-operator $`t`$, which is supposed to act on the pair 23. Finally, $`\varphi `$ is the initial channel state composed of the deuteron state and a momentum eigenstate of the projectile neutron.That integral equation can be solved precisely in a partial wave decomposition in momentum space. For details see .

The operator $`U`$ for elastic scattering is given in terms of the amplitude $`T`$ by quadrature as follows

$`U`$ $`=`$ $`PG_0^1\varphi +PT+V_4^{(1)}(1+P)\varphi +V_4^{(1)}(1+P)G_0T.`$ (2)
We want to solve the integral equation (1) directly at the threshold of nd scattering. This is for zero initial relative momentum $`\stackrel{}{q}_0`$ of the projectile and will then lead directly to the scattering length. For the convenience of the reader we briefly sketch the necessary steps . Our partial wave momentum space basis is denoted by $`|pq\alpha >`$, where p and q are the magnitudes of standard Jacobi momenta and $`\alpha `$ a string of angular momentum and isospin quantum numbers (see ). Then for $`\stackrel{}{q}_0`$ in z- direction we define the auxiliary amplitude

$`U_{\alpha ,\lambda I}(p,q)`$ $`=`$ $`{\displaystyle \underset{m,m_d}{}}{\displaystyle \frac{\sqrt{4\pi \widehat{\lambda }}}{\widehat{J}}}(\lambda 0{\displaystyle \frac{1}{2}}m|Im)(j_dm_dIm|Jm_d+m)<pq\alpha |U|\varphi >`$ (3)
for the projectile nucleon with orbital angular momentum $`\lambda `$ ($`\widehat{\lambda }2\lambda +1`$) and total angular momentum $`I`$ combined with the deuteron total angular momentum $`j_d=1`$ to total 3N angular momentum $`J`$. Out of that amplitude one obtains the partial wave projected nd elastic scattering amplitude

$`U_{\lambda ^{}I^{},\lambda I}^J`$ $`=`$ $`{\displaystyle \underset{l^{}}{}}{\displaystyle p_{}^{}{}_{}{}^{2}𝑑p^{}\varphi _l^{}(p^{})U_{\alpha _{}^{}{}_{d}{}^{},\lambda I}(p^{},q_0)}.`$ (4)
where $`\varphi _l(p)`$ are the s- and d-wave components of the deuteron and $`\alpha _d`$ contains the deuteron quantum numbers.

Finally, the projectile spin and the deuteron spin can be combined to the total spin $`\mathrm{\Sigma }`$ and one obtains

$`U_{\lambda ^{}\mathrm{\Sigma }^{},\lambda \mathrm{\Sigma }}^J`$ $`=`$ $`{\displaystyle \underset{I,I^{}}{}}\sqrt{\widehat{\mathrm{\Sigma }}^{}\widehat{I}^{}}()^{JI^{}}\left\{{\displaystyle \genfrac{}{}{0pt}{}{\lambda ^{}}{j_d}}{\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}}{J}}{\displaystyle \genfrac{}{}{0pt}{}{I^{}}{\mathrm{\Sigma }^{}}}\right\}\sqrt{\widehat{\mathrm{\Sigma }}\widehat{I}}()^{JI}\left\{{\displaystyle \genfrac{}{}{0pt}{}{\lambda }{j_d}}{\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}}{J}}{\displaystyle \genfrac{}{}{0pt}{}{I}{\mathrm{\Sigma }}}\right\}U_{\lambda ^{}I^{},\lambda I}^J.`$ (5)
The $`S`$\- matrix element is then given as

$`S_{\lambda ^{}\mathrm{\Sigma }^{},\lambda \mathrm{\Sigma }}^J`$ $`=`$ $`\delta _{\lambda ^{}\lambda }\delta _{\mathrm{\Sigma }^{}\mathrm{\Sigma }}i{\displaystyle \frac{4\pi }{3}}mq_0(i)^{\lambda ^{}\lambda }U_{\lambda ^{}\mathrm{\Sigma }^{},\lambda \mathrm{\Sigma }}^J`$ (6)
leading to the doublet and quartet scattering lengths for $`q_0=0`$

$`{}_{}{}^{2}a_{nd}^{}`$ $`=`$ $`{\displaystyle \frac{2\pi }{3}}mU_{0\frac{1}{2},0\frac{1}{2}}^{\frac{1}{2}}`$ (7)
$`{}_{}{}^{4}a_{nd}^{}`$ $`=`$ $`{\displaystyle \frac{2\pi }{3}}mU_{0\frac{3}{2},0\frac{3}{2}}^{\frac{3}{2}}.`$ (8)
One also defines a coherent scattering length $`b_{nd}`$ as

$`b_{nd}`$ $`=`$ $`{\displaystyle \frac{m_n+m_d}{m_d}}[({\displaystyle \frac{1}{3}}){}_{}{}^{2}a_{nd}^{}+({\displaystyle \frac{2}{3}}){}_{}{}^{4}a_{nd}^{}]`$ (9)
We defer the special form of the Faddeev integral equation (1) at $`q_0=0`$ to the appendix. It is free of singularities and therefore as easily solved as a bound state problem. Also the explicit form of the elastic amplitude for $`q_0=0`$ is given there.

## III Results

We used the NN forces CD Bonn 2000, AV18, Nijm I, II and Nijm 93 alone or various combinations with the three-nucleon forces Urbana IX, the older Tucson-Melbourne force and the modified one TM99. When we combine the Urbana IX 3NF with CD Bonn 2000 the strength of the repulsive part of this 3NF has been reduced by multiplying it with the factor 0.812 in order to get the proper $`E_{{}_{}{}^{3}H}`$.

Due to their nonnegligible influence on the nd scattering lengths, we took special care of electromagnetic interactions. In the case of the AV18 potential it is clear how to separate the strong AV18 force from the electromagnetic parts because both are well defined and added together in fitting the total force to the NN data. In case of the np system the EMI’s are given in Eqs. (11), (12), and (15) of Ref.  and in Eq. (16) for the nn system (for the np system we did not include the very small class IV charge asymmetric force $`\stackrel{}{L}1/2(\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j)`$. Also we neglected the energy dependence of the $`\alpha ^{}`$). This is different for the CD Bonn 2000 and Nijmegen potentials, which were fitted directly to the NN data without adding to them electromagnetic interactions, with the exception of the point Coulomb force in case of the pp system. Therefore to define the strong forces in the particular NN system one needs to subtract the corresponding EMI, which we assumed to be given as in ref. . To be precise for the np system we subtract from the np CD Bonn and the Nijmegen forces the np EMI’s as defined above. Similarly for the nn system we subtract from CD Bonn the MMI as defined above. Since we also want to see the effect of replacing the strong nn force by the strong pp force we have also to define the strong pp CD Bonn and Nijmegen forces. In this case we subtract from those potentials the pp EMI’s as given in Eqs. (3)-(8) of Ref.  without the leading 1 in $`F_c(r)`$ from Eq. (10) of Ref. , since the point Coulomb force has been taken into account for those potentials fitting to the pp data.

Before we report on our results we give some comments on our numerical accuracy. As usual the partial wave decomposition is truncated at a certain total two-body angular momentum $`j_{max}`$. Fig. 1 documents the convergence of $`{}_{}{}^{2}a_{nd}^{}`$ as function of $`j_{max}`$. This shows that we reached an accuracy of about three digits. This is for CD Bonn  alone and valid also for the other NN forces alone. Adding a three-nucleon force we were limited to $`j_{max}=5`$ due to computer resources. Nevertheless, as Fig. 2 documents, the convergence reached for $`{}_{}{}^{2}a_{nd}^{}`$ is two digits. In case of $`{}_{}{}^{4}a_{nd}^{}`$ with NN forces alone we reach 4 digits convergence and including a 3N force an accuracy close to that. This is documented in Fig. 3. The other numerical ingredients (discretization in the momenta) are safely under control. In all calculations we took into account charge dependence of the NN forces using a simple “$`\frac{2}{3}t_{pp(nn)}+\frac{1}{3}t_{np}`$” rule to generate t-matrices in isospin $`t=1`$ 2N states . The total isospin $`T=3/2`$ 3N states have been neglected . We checked that their inclusion does not change $`{}_{}{}^{4}a_{nd}^{}`$ up to the fifth digit and the change of $`{}_{}{}^{2}a_{nd}^{}`$ is of the order of $`0.1\%`$. The triton binding energies have been obtained using $`j_{max}=6`$. They are accurate to $`2`$ keV.

As an overview we show all our results for $`{}_{}{}^{2}a_{nd}^{}`$ and $`E_{{}_{}{}^{3}H}`$ in Fig. 4. We see a group of results in the right half of the figure based on NN forces alone and another group close to the experimental area including 3N forces. We performed several investigations. First we take CD Bonn 2000 as it is (fitted to the NN data) and use the np-nn force combinations appropriate to the nd system. In this way the EMI’s in the np and nn systems are effectively included inside the strong forces. In case of AV18 we keep all electromagnetic corrections as in  except the energy dependence of $`\alpha ^{}`$ (MMI’s for the nn system and the MMI’s together with the one photon Coulomb term $`V_{C_1}(np)`$ for the np system). The corresponding two predictions are shown as stars in Fig. 4.

Secondly we want to see the effect of replacing the strong nn forces by the strong pp forces. The difference between nn and pp strong forces is mostly located in the different scattering lengths $`a_{nn}`$ and $`a_{pp}`$ (strong) and will therefore give some information how changes in $`a_{nn}`$ will show up in changes of $`{}_{}{}^{2}a_{nd}^{}`$. Since thereby we do not want to change the EMI’s we keep in case of AV18 the nn MMI. For the CD Bonn and the Nijmegen potentials the strong pp potentials are defined as above and the nn MMI (as for AV18) is added. The results are shown as the 5 open circles in Fig 4. As seen in case of CD Bonn 2000 and AV18 the 3H binding energy is decreased and $`a_{nd}`$ increased. For the Nijmegen potentials such a comparison is not done, since no nn forces have been introduced. The effects on the 3N binding energies are known. These two first investigations provide theoretical predictions for the nd scattering lengths and triton binding energy including all electromagnetic interactions similarly as is the case for the measured values.

In the next two investigations we address the effects of the electromagnetic interactions themselves by switching them off while generating theoretical predictions. For the AV18 potential we take just the np-nn and np-pp strong force combinations alone, while in the cases of CD Bonn 2000, Nijm I, II, and 93 we only use the corresponding strong forces stripped off from the EMI’s, as described above. The resulting theoretical predictions are shown as pluses and squares in Fig. 4 for the np-nn and np-pp combinations, respectively. Again the binding energy is decreased and $`a_{nd}`$ increased.

The individual results of these four investigations are summarized in Fig. 4 also as dashed (np-nn with EMI’s), dotted (np-pp with EMI’s), solid (np-nn), and dashed-dotted (np-pp) straight lines fitted in a $`\chi ^2`$ sense. These lines include also the corresponding results including 3NF’s (see below). We see a small shift of the lines under exchanges of nn versus pp forces, but a more significant shift if the electromagnetic forces are switched off. Though the two curves (dashed and dotted) for the cases when the electromagnetic forces are added come close to the experimental range spanned by the uncertainty in $`{}_{}{}^{2}a_{nd}^{}`$, they miss it clearly. When the electromagnetic forces are switched off the np-nn (solid) and np-pp (dashed-dotted) lines go through the experimental point well inside the $`{}_{}{}^{2}a_{nd}^{}`$ error bar.

Now we want to regard our results in more detail as displayed in Table I and in the inset of Fig. 4. The theory has to be finally compared to the experimental values, which are $`{}_{}{}^{2}a_{nd}^{}=(0.65\pm 0.04)`$ fm , $`{}_{}{}^{4}a_{nd}^{}=(6.35\pm 0.02)`$ fm , and $`b_{nd}=(6.669\pm 0.003)`$ fm .

The results in Table I are grouped into NN force predictions only and selected combinations with the 3N forces TM, TM99 and Urbana IX. For each potential or potential combination we show the results for the various scattering lengths and the 3H binding energies. This is given for the np-pp NN forces, with (without) EMI’s in the first (second) row. For AV18 and CD Bonn 2000 we also show the results for np-nn forces with (third row) and without (fourth row) electromagnetic interactions. Note that in case of the np-nn forces including EMI’s (as described above) the combinations with TM99 and Urbana IX are fitted well to the experimental value -8.48 MeV of the 3H binding energy. This is also the case for Nijm I and II, which, however, refers to the np-pp forces. For the older TM 3N force we did not perform a precise (re)fit and the results are only included in view of investigating, whether a straight line correlation between $`{}_{}{}^{2}a_{nd}^{}`$ and $`E_{{}_{}{}^{3}H}`$ exists. A glance at Fig. 4 tells that the individual results scatter around the four straight lines. Thus obviously no straight line correlation exist (this has been known before, though for some older calculations the numerical accuracy might be not sufficient to give a reliable judgement).

Let us now concentrate on the group of results with 3N forces. These are displayed in the inset of Fig. 4. We see four results (stars) for the np-nn forces including TM99 or Urbana IX, where the binding energy has been exactly fitted but where the $`{}_{}{}^{2}a_{nd}^{}`$ value is too small. These are the results achieved under the supposedly most realistic assumptions in this paper. If one switches off the electromagnetic interaction (pluses) the binding energy increases and interestingly $`{}_{}{}^{2}a_{nd}^{}`$ moves to larger values. Regarding all results, the inclusion of the electromagnetic force in our studies shows that they cause shifts of up to about 40 keV less binding energy and of up to about 0.04 fm decrease in $`{}_{}{}^{2}a_{nd}^{}`$. In no case studied the experimental value of $`{}_{}{}^{2}a_{nd}^{}`$ is reproduced for np-nn and np-pp strong forces combined with different 3NF’s with exception of np-pp AV18 combined with TM 3NF, for which the theoretical prediction lies at the lower limit of the error bar.

As one learned from the approach in chiral perturbation theory , where two parameters are needed to fix the short range 3N forces at NNLO and consequently two 3N observables to adjust them, one could foresee that the straight lines in Fig. 4 could only by accident pass through the experimental region. For the conventional forces used in this paper, one can think of additional 3N force diagrams (the most obvious one the $`\pi \rho `$ exchange) where a sufficient number of parameters would be available to fit both, $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$.

Going back to Table I we see that $`{}_{}{}^{4}a_{nd}^{}`$ sticks always close to the value 6.34 for the np-pp and np-nn NN force choices, without or with 3N forces and with EMI’s included. This is well within the experimental $`{}_{}{}^{4}a_{nd}^{}`$ error bar. Interestingly, the electromagnetic interactions increase $`{}_{}{}^{4}a_{nd}^{}`$ in nearly all cases by about 0.02 and the pure strong force predictions lie always outside the experimental error bar.

Finally, one can confront theory to the very precisely known experimental value of the coherent scattering length $`b_{nd}`$ . Clearly the supposedly most realistic dynamics (nn-np NN forces plus TM99 or Urbana IX 3NF’s) misses that value. As can be seen from Table I, when electromagnetic interactions are included the np-pp force combination reaches the experimental value in case of the AV18 and CD Bonn 2000 potentials combined with Urbana IX and AV18 with TM99. However, this agreement is accidental and caused by the corresponding decrease in $`{}_{}{}^{3}H`$ binding.

## IV Summary and Outlook

A recently performed precise neutron interferometric measurement of the nd coherent neutron scattering length  and a planned precision measurement of the doublet nd scattering length  stimulated us to investigate the theoretical predictions of that quantity for the high precision NN forces CD Bonn 2000, AV18, Nijm I, II and 93 in combination with currently popular 3N force models. These are the modified $`2\pi `$-exchange Tucson-Melbourne (TM99) and the Urbana IX 3N forces. We have chosen several NN and 3NF combinations, which are separately adjusted to the $`{}_{}{}^{3}H`$ binding energy. For NN forces alone with and without EMI’s we recovered the approximate correlation between $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$, but the scatter around a thought straight line (Phillips line) inside the band spanned by the 4 lines in Fig. 4 is quite strong. Adding 3N forces shifts the values into the neighbourhood of the experimental range of $`{}_{}{}^{2}a_{nd}^{}`$, but misses the experimental value including its error bar in all cases, when electromagnetic forces are included. The inset of Fig. 4 clearly shows that for equal or nearly equal 3H binding energies $`{}_{}{}^{2}a_{nd}^{}`$ can vary significantly and vice versa.

Thus one has to conclude that $`{}_{}{}^{2}a_{nd}^{}`$ has to be considered as a low energy observable, which is independent from the $`{}_{}{}^{3}H`$ binding energy. This observation has been found before in approaches based on pure effective field theory (pionless formulation) and on chiral perturbation theory (including pion degrees of freedom). Thus in future investigations, adjusting both observables, $`E_{{}_{}{}^{3}H}`$ and $`{}_{}{}^{2}a_{nd}^{}`$, for conventional nuclear forces will require more flexibility in the choice of 3N forces. Adding more mechanisms (on top of the $`2\pi `$-exchange) for 3N forces should be no obstacle. This is a step already performed in the effective theory approaches .

We also investigated the effects on $`{}_{}{}^{2}a_{nd}^{}`$ resulting from electromagnetic interactions given in . The effects on $`{}_{}{}^{2}a_{nd}^{}`$ and even $`{}_{}{}^{4}a_{nd}^{}`$ are noticeable. For $`{}_{}{}^{2}a_{nd}^{}`$ including the electromagnetic interactions reduces its value by up to 0.04 fm. It is interesting to note that $`{}_{}{}^{4}a_{nd}^{}`$ is perfectly stable under all exchanges of nuclear forces studied in this paper but the electromagnetic interactions affect its value, though only in the 3’rd digit. However, only when EMI’s are included the experimental value is reproduced.

The effects of adding the electromagnetic interactions on the $`{}_{}{}^{3}H`$ binding energy are well known and can reach shifts of up to 40 keV less binding energy.

Due to the failure to describe $`{}_{}{}^{2}a_{nd}^{}`$ also the recently newly measured coherent scattering length $`b_{nd}`$ cannot be reproduced theoretically. The good reproduction of $`{}_{}{}^{4}a_{nd}^{}`$ by all interactions and the small error bar of the coherent scattering length suggests that the value of the doublet nd scattering length might be somewhat smaller than the presently one, namely around 0.63 fm, different from the present experimental value of $`{}_{}{}^{2}a_{nd}^{}`$. This strongly calls for a new, more precise measurement.

Since the scattering lengths are (extreme) low energy observables, it appears that the mentioned effective theory approaches are the most adequate ones. Because one works there below a certain momentum cut off, which is smaller than the nucleon mass, they allow also to incorporate relativistic effects in a well defined and converged manner. Also 3N forces appear in those approaches in a well organised manner, according to a certain power counting scheme, and are consistent to the NN forces. In other words, one can take into account all these subtle effects, relativity, 3N forces, isospin breaking, in a well controlled and systematic manner. In conventional approaches on the other hand, which include a lot of phenomenological parametrisations and where no momentum cut-off is used, a reliable treatment of relativistic effects poses still a problem and the choices of 3N force mechanisms are also quite open. Therefore in conventional approaches physically reliable predictions to $`{}_{}{}^{2}a_{nd}^{}`$ will very likely remain a challenge for quite some time.

## V Appendix

This appendix summarizes various expressions exactly at the nd threshold $`q_0=0`$. The first part of the driving term in Eq.(1) turns out to be

$`<pq\alpha |tP|\varphi >`$ $`=`$ $`\delta _{\lambda _0,0}{\displaystyle \underset{l_\alpha ^{}l_0I_0}{}}<pl_\alpha |t^\alpha ({\displaystyle \frac{3}{4m}}q^2)|{\displaystyle \frac{q}{2}}l_\alpha ^{}>\phi _{\alpha _0}(q)`$ (10)
$`(1m_dI_0m_n|JM)(\lambda _00{\displaystyle \frac{1}{2}}m_n|I_0m_n)\sqrt{{\displaystyle \frac{\widehat{\lambda }_0}{4\pi }}}.`$ (11)
The quantities with index 0 refer to the initial state.

The kernel applied on $`T`$ is given as

$`<pq\alpha |tG_0PT`$ $`=`$ $`{\displaystyle q_{}^{\prime \prime }{}_{}{}^{2}𝑑q^{\prime \prime }}`$ (12)
$`{\displaystyle _1^{+1}}𝑑x{\displaystyle \frac{\underset{l_\alpha ^{}}{}\frac{m}{qq^{\prime \prime }}\frac{t^{l_\alpha ,l_\alpha ^{}}(p,\pi _1;\frac{3}{4m}q^2)}{\pi _1^{l_\alpha ^{}}}\underset{\alpha ^{\prime \prime }}{}\frac{G_{\overline{\alpha }\alpha ^{\prime \prime }}(q,q^{\prime \prime },x)}{\pi _2^{l_{\alpha ^{\prime \prime }}}}<\pi _2q^{\prime \prime }\alpha ^{\prime \prime }|T}{x_0x}}`$ (13)
with

$`x_0`$ $``$ $`{\displaystyle \frac{k_d^2q^{\prime \prime 2}q^2}{qq^{\prime \prime }}}`$ (14)
and $`\overline{\alpha }`$ contains the same quantum numbers as $`\alpha `$ with the exception of $`l_\alpha `$ replaced by $`l_\alpha ^{}`$.

For our notation see . The deuteron binding energy is written as $`(k_d^2/m)`$. The remaining parts related to $`V_4^{(1)}`$ can be worked out correspondingly and can be found in . Evaluating the elastic scattering amplitude one needs it at $`q=q_0`$ (see Eq.(4)). Therefore the point $`q=q_0=0`$ was included. Then Eq.(13) simplifies to

$`<pq=0\alpha |tG_0PT`$ $`=`$ $`2m\delta _{\lambda _\alpha ,0}{\displaystyle q_{}^{\prime \prime }{}_{}{}^{2}𝑑q^{\prime \prime }}`$ (15)
$`{\displaystyle \underset{l_\alpha ^{}}{}}{\displaystyle \frac{t^{l_\alpha ,l_\alpha ^{}}(p,q^{\prime \prime };0)\underset{\alpha ^{\prime \prime }}{}2^{l_{\alpha ^{\prime \prime }}}g_{\overline{\alpha }\alpha ^{\prime \prime }}^{0l_\alpha ^{}0l_{\alpha ^{\prime \prime }}0}<\frac{1}{2}q^{\prime \prime }q^{\prime \prime }\alpha ^{\prime \prime }|T}{k_d^2q^{\prime \prime 2}}}.`$ (16)
One ends up with the elastic scattering amplitude at threshold:

$`U_{\lambda ^{}I^{},\lambda I}^{J^\pi }`$ $`=`$ $`{\displaystyle \frac{2k_d^2}{m}}\delta _{\lambda ,0}\delta _{\lambda ^{},0}g_{\alpha _d^{}\alpha _d}^{00000}{\displaystyle \underset{l,l^{}}{}}{\displaystyle \frac{\phi _l^{}(p)}{p^l^{}}}|_{p=0}{\displaystyle \frac{\phi _l(p)}{p^l}}|_{p=0}`$ (17)
$`+`$ $`\delta _{\lambda ^{},0}{\displaystyle \underset{l^{},\alpha ^{\prime \prime }}{}}2^{l^{\prime \prime }+1}g_{\alpha _d^{}\alpha ^{\prime \prime }}^{0l^{}0l^{\prime \prime }0}{\displaystyle q_{}^{\prime \prime }{}_{}{}^{2}𝑑q^{\prime \prime }\phi _l^{}(q^{\prime \prime })<\frac{1}{2}q^{\prime \prime }q^{\prime \prime }\alpha ^{\prime \prime }|T}`$ (18)
$`+`$ $`{\displaystyle \underset{l^{}}{}}{\displaystyle p_{}^{}{}_{}{}^{2}𝑑p^{}\phi _l^{}(p^{})\{V_4^{(1)}(1+P)\varphi +V_4^{(1)}(1+P)G_0T\}_{\alpha _d^{},\lambda I}(p^{},0)}.`$ (19)
The geometrical coefficients $`g_{\alpha \alpha ^{}}^{kl_1l_2l_{}^{}{}_{1}{}^{}l_{}^{}{}_{2}{}^{}}`$ arise from the permutation operator P and are given by Eq.(A19) in .

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft (J.G.,R.S.), by the NFS grant No. PHY0070858, and by the NATO grant PST.CLG.978943. R.S. acknowledges financial support of the Foundation for Polish Science. The authors would like to thank Dr. E. Epelbaum for very constructive discussions. The numerical calculations have been performed on the Cray T90, SV1 and T3E of the NIC in Jülich, Germany.