text:
# Orbital-resolved Soft X-Ray Spectroscopy in NaV2O5

## I Introduction

Orbital degrees of freedom play a crucial role in many technologically interesting materials, such as superconductors, magnetic materials, and other strongly correlated compounds. For instance, the spin-orbit coupling effects determine whether the magnetic crystalline anisotropy in magnetic thin films favors magnetic moments normal or parallel to the surface of the sample. Such important effects motivate many investigations that probe orbitals in real materials. Soft x-ray spectroscopy (SXS)  has a unique capability to reveal spin and orbital information specific to one element. By employing x-ray magnetic circular dichroism, one can measure spin and orbital moments; using x-ray magnetic linear dichroism, one can image antiferromagnetic domains at sub-micrometer scales. Despite these tremendous achievements, imaging atomic orbitals has not often been done. The first attempt at the C $`K`$ edge of graphite was done by Rosenberg et al. A more elaborate study was performed a decade later with an intense synchrotron source. On surfaces, the first investigation was made by Nilsson et al. They used $`s`$-polarized soft x-rays to selectively excite $`\pi `$ orbitals of CO adsorbed on Ni and used $`p`$-polarized light at off-normal angles to excite both $`\sigma `$ and $`\pi `$ orbitals, from which the orbital information could be extracted. However, obtaining such valuable orbital information in bulk materials is still a challenge for SXS.

In this work, we consider $`\mathrm{NaV}_2\mathrm{O}_5`$, a fascinating transition-metal compound that features three inequivalent oxygen sites: O1, O2 and O3 \[cf Fig. 1(a)\]. At 34 K, this material undergoes a spin-Peierls-like phase transition, which has attracted great interest. $`\mathrm{NaV}_2\mathrm{O}_5`$ is orthorhombic with symmetry of $`P_{mmn}`$ and lattice constants $`a`$=11.316 $`\mathrm{\AA }`$, $`b`$=3.611 $`\mathrm{\AA }`$ and $`c`$=4.797 $`\mathrm{\AA }`$. It is a ladder structure compound along the $`b`$ axis with V-O1-V bonds as its rungs extending along the $`a`$ axis. The neighboring ladders along the $`b`$ axis are connected by shifting a half-unit cell along the $`b`$ axis. These ladders form quasi-two dimensional layers between which the sodium atoms are intercalated. One O1, three oxygen O2 and one apical O3 atoms form a pyramid inside which one vanadium atom is situated. These three different oxygen atoms present a challenge to site-specific SXS, because one desires to obtain information, not only about different orbitals of atoms of different elements, but also about different orbitals on the same site.

Specifically, we report orbital-resolved x-ray absorption measurements, using the total electron yield method (TEY), on $`\mathrm{NaV}_2\mathrm{O}_5`$ at the O $`K`$ (1s) and V $`L_3`$ (2p3/2) edges. By changing the angle $`\theta `$ between the electric field-polarization direction and certain crystal axes of $`\mathrm{NaV}_2\mathrm{O}_5`$, we can determine characteristic orbital energies, not only for orbitals on inequivalent O sites, but also for different orbitals on the same O site. With the aid of theoretical calculations, we are able to pinpoint the origin of the angular dependence of the O $`K`$ edge absorption spectrum, distinguish different contributions thereto, and accurately determine the crucial energy ordering of the different orbitals from three inequivalent oxygen atoms. Combining different scans along different axes, we are able to project out spatial orientations of oxygen orbitals in real space. Thus, our study highlights an important capability of SXS that has not been explored in great detail.

Our experiment reveals a surprising difference between the vanadium and oxygen edges. At the V $`L_3`$ edge, the spectra show a very small dependence on $`\theta `$, but at the O $`K`$ edge there is a large dependence. Our calculations are scalar-relativistic, so that we model the V $`L_{2,3}`$ manifold without spin-orbit effects (including $`p_{1/2}p_{3/2}`$ splitting) taken into account. However, this description should still be able to address many aspects of a polarization dependence, which depends on states near the Fermi level. The theory shows that, subject to the dipole selection rule $`\mathrm{\Delta }l=\pm 1`$ and $`\mathrm{\Delta }m=0,\pm 1`$, transitions from $`s`$ orbitals only have one initial state, i.e., $`m=0`$, while transitions from $`p`$ orbitals have three initial states, i.e., $`m=0,\pm 1`$. Polarization dependences for absorption by electrons in each of the initial states compensate each other and ultimately lead to a very small polarization dependence. This explains the observed trend for the V $`L_3`$ (2p3/2) edge.

The remainder of this article is arranged as follows. A brief review of our experimental measurements is described in Section II. In Section III, we present our theoretical scheme. The results and discussions are shown in Sec. IV. Finally, we conclude our paper in Section V. In a companion paper, we analyze emission spectra and experimentally determine the occupied $`p`$ density of states from each oxygen site. This paper is devoted entirely to the analysis of absorption spectra.

## II Experimental measurements

Single crystals of $`\mathrm{NaV}_2\mathrm{O}_5`$ were prepared at the Oak Ridge National Laboratory. The crystals grow as rectangular platelets with the c axis normal to the plate surface, the a axis along the short side of the rectangle and b along the long side, as shown in Fig. 1(b). The samples were cleaved along the c axis. The experimental measurements were done at Beamline 8.0 of the Advanced Light Source at the Ernest Orlando Lawrence Berkeley National Laboratory. The geometry of the measurements is indicated in Fig. 1(b). Horizontally polarized x-rays from an undulator and monochromator were incident on the sample. In the TEY detection mode, the x-ray absorption spectrum is measured by collecting the electrons that escape from the sample. Two samples were used in the experiment, and mounted at right angles to each other. Rotation of the sample holders scanned the polarization vector of the incident x-rays in the a-c plane for one sample and in b-c plane for the second sample.

$`\mathrm{NaV}_2\mathrm{O}_5`$ is an ideal system, in which to reveal a polarization dependence, because its structure is highly anisotropic. In Fig. 2, we show the normalized absorption spectra of a polarization scan from along the b axis to along the c axis, specifically from $`\theta =7.5^{}`$ to $`\theta =75^{}`$ with a step of 7.5. For each angle, the normalization is done with respect to the total area between 514 eV and 536 eV because this energy region covers both the V-L and O-K edges.

We note that the TEY method has a finite sampling depth of about 100 $`\mathrm{\AA }`$. This suggests that our electron yield arises from a local spot of the sample. To illuminate the sample in the same place during measurement and to ensure data consistency, we adjusted the sample position after each change of its orientation angle. We found that data obtained using this technique are reproducible within 5 %.

We first consider the vanadium $`L_3`$ edge around 518 eV \[Fig. 2(a)\]. When we scan the polarization from along the b axis to along the c axis, the spectra change very little. We observed a similar behavior when we scanned from along the a axis to along the c axis (not shown). However, for the oxygen $`K`$ edge, we see a totally different picture \[Fig. 2 (b)\]. Increasing the angle $`\theta `$ from 7.5 to 75 strongly changes the spectra. Three main peaks A, B and C, located at 530.6 eV, 532.1 eV and 532.8 eV, respectively, can be clearly identified. Peak A’s intensity decreases sharply and eventually disappears as the polarization moves to along the c axis, which is normal to the sample surface. In the meantime, intensities of peaks B and C increase strongly and become the two main peaks on the spectra. To have a clear view of the intensity change, in Fig. 3 we plot the intensities versus angle $`\theta `$ for those three peaks. Their variations can be nicely fitted to a $`\mathrm{cos}2\theta `$ dependence, which is the simplest possible variation (see the solid lines in Fig. 3).

## III Theoretical framework

Our theoretical calculation to simulate the experiment is based on pseudopotential, local-density-approximation (LDA) calculations. Electron states are found by solving the self-consistent Kohn-Sham (KS) equations,

$$\widehat{H}\mathrm{\Psi }_{n𝐤}(𝐫)=\{\frac{\mathrm{}^2}{2m}^2+\widehat{V}_{\mathrm{eff}}\}\mathrm{\Psi }_{n𝐤}(𝐫)=E_{n𝐤}\mathrm{\Psi }_{n𝐤}(𝐫)$$
(1)
Here $`\mathrm{\Psi }_{n𝐤}`$ is the Bloch state for an electron in band $`n`$ at crystal wave-vector $`𝐤`$. The effective potential is given by

$$\widehat{V}_{\mathrm{eff}}=V_{\mathrm{ion}}+V_\mathrm{H}(𝐫)+V_{\mathrm{xc}}(𝐫).$$
(2)
The three terms account for the pseudopotential interaction of electrons with ion cores $`V_{\mathrm{ion}}`$, the Hartree potential $`V_\mathrm{H}(𝐫)`$ that is obtained from the electron charge density, given by

$$n(𝐫)=\underset{n𝐤}{\overset{\mathrm{occ}^{}\mathrm{d}.}{}}|\mathrm{\Psi }_{n𝐤}(𝐫)|^2$$
(3)
and exchange and correlation effects $`V_{\mathrm{xc}}`$. We use the Ceperley-Alder form for the exchange term $`V_{\mathrm{xc}}=\delta E_{\mathrm{xc}}[n]/\delta n(𝐫)`$. The ionic potential $`V_{\mathrm{ion}}`$ is a norm-conserving pseudopotential computed within the Hartree-Fock approximation, and with core-valence correlation effects included using the core-polarization-potential method. As basis sets, we use plane waves or optimized basis functions.

To compute an x-ray absorption spectrum, we use the expression,

$$S(\omega )2\pi \underset{f}{}|f|𝐩𝐀|i|^2\delta (\omega +E_iE_f)$$
(4)
where $`𝐩𝐀`$ is the electron-photon interaction approximated by a dipole operator; $`|i`$ is the initial core state with energy $`E_i`$, $`|f`$ is the final state with energy $`E_f`$, and $`\omega `$ is the incident photon energy. In our present analysis, which mainly deals with determining the symmetries of orbitals, we neglect electron-core hole interactions in the final state.

## IV Results and discussions

### A Oxygen $`K`$ (1s) edge absorption spectrum

In the experiment, we always rotated our polarization direction from along the b axis to along the c axis in the way already described. Theoretically, we did calculations in the same way. In Fig. 4(a), we show theoretical x-ray absorption spectra at the O $`K`$ edge, where three peaks, A, B and C, are also clearly observed with different intensities as a function of polarization. The spectra were calculated with $`12\times 12\times 12`$ k-point meshes. The results are found well-converged. The angular dependence of these peaks is in good agreement with the experimental observations: Peak A decreases as the polarization moves away from along the b axis, whereas both peaks B and C increase. The change also follows a cosine function of $`\theta `$, consistent with the experimental observation. Note that the experimental measurements cannot directly determine which oxygen orbitals contribute to these main peaks, whereas calculations should provide important insight into this.

Figures 4(b-d) show the respective contributions from three inequivalent oxygen sites. To simplify our notation, we let Om(Pn) denote the $`n`$th orbital of oxygen Om that participates the absorption process. That is, the local, partial density of states for such an orbital can play a significant role in absorption as a function of excitation energy, and this role can be strongest at some “characteristic” energy. We first consider peaks B and C. By comparing Fig. 4(a) with Figs. 4(b), 4(c) and 4(d), we infer that peaks B and C in Fig. 4(a) and Fig. 2(b) originate from O1(P3), O2(P4), O2(P5), and O3(P5) orbitals, whose characteristic orbital energies (expressed relative to the O 1s level) are 533.0 eV, 533.1 eV, 533.7 eV and 533.1 eV, respectively.

When considering peak A in Figs. 4(a) and 2(b), we notice that all three inequivalent atoms contribute but are favored in different ranges of energy. Figure 4(b) illustrates the main contributions from O1. It is important to note that only O1(P1) appears at 530.2 eV, which helps us to resolve different orbital contributions (see next paragraph). The peak intensity decreases as the polarization direction deviates from the b axis, in a similar way as peak A in Fig. 2(b) or Fig. 4(a). O1(P1) is immediately followed by O3(P1) at 530.6 eV with an energy difference between them of about 0.4 eV, \[cf. Fig. 4 (d)\]. The O3(P1) intensity changes in a similar manner as does that for O1(P1). After the O3(P1) orbital, the O2(P3) orbital gradually appears around 531 eV to 532 eV \[cf. Fig. 4(c)\]. Because of their unique positions, O2 orbitals are distorted, and not along one particular axis. Thus, the angular dependence of their contribution to absorption is weaker than for orbitals on O1 and O3.

These results show that the measured peak A has contributions from the O1(P1), O3(P1), and part of O2(P3) orbitals. This could prove crucial for our tentative explanation of emission spectra that is given subsequently in Ref. 16. In particular, because our theory shows that these three orbitals have different energies, it is possible to excite different orbitals at different energies and to resolve them experimentally. This is exactly what we did. Because the energy increases in the order of O1, O3 and O2, we can selectively excite these atoms. When we excite at the O1(P1) orbital edge, we only see one peak in the emission spectrum, because O2 and O3 atoms are not excited, though in the normal fluorescence spectrum, we expect to see two main peaks. Therefore, with theoretical guidance, we demonstrate that we are able to resolve these different orbitals associated not only with different sites such as O1(P1), O3(P1), O2(P3) but also with the same site such as O1(P1) and O1(P3).

By scanning polarization from along the a axis to along the b axis and from along the a axis to along the c axis, we can project out orbital orientations in real space. Putting together all the results, we can show some of the orbitals in the real space in Fig. 5. One sees that there are two kinds of O3 orbitals: one is along the a axis, the other is along the c axis, and one O1 orbital along the b axis, two O2 orbitals along the b and a axes. Note that, unlike charge-density plots for orbitals, the orbitals presented here are experimentally accessible. Thus, we demonstrate that the soft x-ray spectroscopy has the capability of resolving orbitals from the same site.

### B Vanadium $`L_3`$ (2p3/2) edge

The x-ray absorption spectrum at the V $`L_3`$ edge is rather unique and shows a very weak angular dependence. At room temperature, there is only one inequivalent vanadium atom. In our scalar-relativistic calculations, excitation of V 2p electrons can involve three electron initial states: $`m=0,\pm 1`$. Our absorption spectrum intrinsically involves a sum over excitations with the indicated initial states, and the experimentally observed spectrum is an analogous sum at the V $`L_3`$ edge. We first show the results of the first channel in Fig. 6(a). We note that with the incident polarization turning away from the b axis, the whole spectrum clearly changes. If this were the dominant excitation channel, then one would see a strong angular dependence in experiment. However, our experiment showed otherwise, which implies that we should take into account the other channels. In Fig. 6(b), we plot the change because of the $`l=1`$, $`m=\pm 1`$ channel. Without spin-orbit coupling, the spectra for both $`l=1,m=+1`$ and $`l=1,m=1`$ channels are same. Two prominent features are: (1) the mean peak strongly overlaps the spectra of the $`l=1`$, $`m=0`$ channel; (2) the portion showing strong angular dependence is on the right-hand side, which is opposite to that in Fig. 6(a). Because of this compensating effect, the total absorption spectrum does not show a strong angular dependence as one can see in Fig. 6(c). This explains our experimental observations.

Now comparing this with the oxygen $`K`$-edge spectra where only one excitation channel, $`l=0`$ and $`m=0`$, participates in the process, we clearly see that appropriate summation over the multiple excitation channels (i.e., corresponding to different electron initial states) obscures the polarization dependence, which ultimately leads to the big difference between the role of polarization at the V $`L_3`$ edge vs the O $`K`$ edge. This finding demonstrates the importance of understanding selection rules in x-ray spectroscopy, and it might help one optimize the polarization direction to reveal the maximal information from x-ray spectroscopy measurements in the future. In particular, it would be interesting to ask whether the present polarization effects for the various components of the V $`L_3`$ spectrum could be resolved by circular dichroism. However, the difficulty is that even with the 2p core spin-orbit coupling, in the absence of spin-polarization of the V-d bands, $`\mathrm{\Delta }m=\pm 1`$ transitions should also be the same.

## V Conclusions

In conclusion, we measured angular dependence of x-ray absorption spectra at the vanadium $`L_3`$ edge and oxygen $`K`$ edge. Rotating the polarization from along the b axis to along the c axis, the spectra at these two edges behave quite differently. There is very little change at the V $`L_3`$ edge, while a dramatic change is observed at the O $`K`$ edge. Such a large change enables us to selectively excite different orbitals belonging to the same site. Thus, we are able to resolve orbitals for the same element using soft x-ray spectroscopy. Our theoretical calculation provides the detailed information to understand the experimental spectra. It clarifies how the predominant first peak at the oxygen $`K`$ edge is from three different orbitals of three different oxygen sites (O1, O2 and O3) and determines the crucial energy ordering of the orbitals on these sites. This also helps explain the emission spectra. Our calculations also show how a compensating effect can arise in the polarization dependence of various excitation channels for the V $`L_3`$ edge, so that the total absorption spectrum at the edges exhibits a small polarization dependence. After we submitted this paper, we noticed three recent papers  that intensively debated whether the anomalous scattering factors in the x-ray diffraction technique could be used to probe the charge ordering in $`\mathrm{NaV}_2\mathrm{O}_5`$. Our study shows that if excited at the V $`L_3`$ edge, the orbital or structural effects are small. Thus, one might explain the anomalous scattering as resulting from the charge ordering, provided that the scattering intensity of the charge ordering is larger than that of the orbital ordering.

Acknowledgment

This research is supported by NSF DMR-9801804 and NSF DMR-0072998. Oak Ridge National Laboratory is managed by UT-Battelle, LLC, for the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. Measurements were carried out at the Advanced Light Source at the Ernest Orlando Lawrence Berkeley National Laboratory, supported U.S. Department of Energy Contract DE-A003-76SF00098. We would like to acknowledge Christen Halloy and the staff at the Joint Institute of Computational Sciences (JICS) at the University of Tennessee at Knoxville, where part of calculation has been done.

Mailing address

#Present address

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6


format: text

text:
# Magnetoconvection and dynamo coefficients:

## 1 Introduction

Magnetic pumping refers to transport of magnetic fields in convective layers which does not result from bulk motion. The existence of magnetic pumping has been known for several decades now. Drobyshevski & Yuferev (drobyshevski74 (1974)) demonstrated in a numerical investigation that magnetic flux in a convective layer is expelled downward if the downflows form connected lanes that separate isolated upflowing regions, an effect known since as topological pumping (cf. also Moffatt moffatt78 (1978), p. 70; Petrovay petrovay91 (1991)). A different but related effect is turbulent pumping. For the two-dimensional case, Drobyshevski (drobyshevski77 (1977)) showed that a poloidal magnetic field is transported against the gradient of the density. Analytical investigations based on the first order smoothing approximation (FOSA) suggest that turbulent pumping also occurs in the presence of a gradient in the intensity of the turbulence, due to the diamagnetic effect (Rädler raedler68 (1968), Roberts & Soward roberts75 (1975), Krause & Rädler krause80 (1980), Kichatinov & Rüdiger kichatinov92 (1992)), or as a result of rotationally induced anisotropies (Kichatinov kichatinov91 (1991)), or through magnetic buoyancy of small-scale fields (Kichatinov & Pipin kichatinov93 (1993)). Numerical simulations of stratified convection showed that dynamo-generated magnetic fields are pumped downward (Nordlund et al. nordlund92 (1992), Brandenburg et al. brandenburg96 (1996)). Tobias et al. (tobias98 (1998), tobias01 (2001)) and Dorch & Nordlund (dorch01 (2001)) performed detailed numerical experiments of magnetic pumping in box simulations of magnetoconvection. They did not, however, measure the pumping velocity which plays a role in mean-field electrodynamics.

Although in most investigations the possible significance of magnetic pumping for stellar dynamo action is stressed, it was in fact rarely incorporated in models. This is true partly because early mean-field solar dynamo models did not appear to ‘need’ any pumping effect since they were able to reproduce the large-scale solar magnetic field rather well using the simplest assumption that there are just alpha effect and differential rotation (for the alpha effect see, e.g., Parker parker55 (1955), or Krause & Rädler krause80 (1980)). Hence the pumping effect, along with many other dynamo coefficients, was thought to be an unnecessary complication. First indications for the usefulness of turbulent pumping came from an $`\alpha \mathrm{\Omega }`$-dynamo simulation of Brandenburg et al. (brandenburg92 (1992)), where downward pumping of the mean magnetic field was found to favour an equatorward migration of magnetic field in the presence of mainly latitudinal differential rotation.

A major problem in solar dynamo theory was the storage of magnetic flux for sufficiently long times (see, e.g., Moreno-Insertis et al. morenoinsertis92 (1992)). It was believed that magnetic buoyancy would remove magnetic field quickly from the convection zone. However, numerical simulations of dynamo action in turbulent stratified convection indicated quite clearly that the effects of magnetic buoyancy can be offset by the downward pumping (Brandenburg & Tuominen brandenburg91 (1991), Nordlund et al. nordlund92 (1992), Brandenburg et al. brandenburg96 (1996), Sect. 4.4). Thus the vertical pumping effect may be important because sufficient amplification of the toroidal magnetic fields generated by the solar dynamo can be achieved only in or near the stably stratified overshoot layer below the convection zone.

Apart from the vertical direction, pumping effects can also exist in other directions. So far though, non-radial pumping effects have not received much attention, although they have interesting consequences for stellar dynamos, as is discussed further below. Their existence has been predicted on the basis of symmetry considerations within the framework of mean-field electrodynamics (Krause & Rädler krause80 (1980)), and some analytical results based on FOSA are presented in Kichatinov (kichatinov91 (1991)) and Rüdiger & Kichatinov (ruediger93 (1993)). No numerical studies of non-radial pumping effects have previously been carried out, as far as we know. The present investigation provides the first numerical evidence for their existence, as well as quantitative results for the corresponding transport coefficients.

The study of pumping effects is part of an effort to calculate dynamo coefficients through numerical simulations. This work is motivated by the idea that progress in understanding stellar dynamo action can be achieved only through a combination of approaches which include numerical simulation as well as mean-field modelling. On the one hand, fully global MHD calculations of stellar dynamos have not been successful so far because computer capacities are insufficient to capture the parameter regime of stellar convection zones. Local MHD simulations of small subvolumes offer the best possibility to approach solar parameters. On the other hand, observations of the large-scale solar magnetic field (solar cycle, butterfly diagram, Hale’s polarity laws) suggest that some form of mean field description is possible. Local numerical investigations of aspects of stellar dynamos may provide useful input information for the construction of mean-field models. Also, the study of dynamo coefficients can be a powerful diagnostic tool to isolate from numerical simulations the various effects that play a role in stellar dynamos and to quantify them. In order to understand the role of dynamo coefficients, consider the averaged MHD induction equation:

$$\frac{\overline{B}}{t}=\mathrm{curl}\left(\overline{U}\times \overline{B}+\eta \mathrm{curl}\overline{B}\right),$$
(1)
Here $`\overline{B}`$ is the mean magnetic field, $`\overline{U}`$ is the mean flow, $`\eta `$ is the magnetic diffusivity, and $`=\overline{u\times b}`$ is the electromotive force (EMF) which depends on $`u=U\overline{U}`$, the fluctuating component of the flow, and $`b=B\overline{B}`$, the fluctuating component of the magnetic field. Mean quantities are either ensemble averages or azimuthal averages. In the present case of a box with periodic boundary conditions in the horizontal directions (see below), horizontal averages can be used instead, and in practice, there is additional temporal averaging. The dynamo coefficients appear in an expansion of the electromotive force in terms of spatial derivatives of the mean magnetic field. In Cartesian coordinates, and assuming the summation convention, this may be written as

$$_i=\alpha _{ij}\overline{B}_j+\beta _{ijk}NewA_j\overline{B}_k+\mathrm{}.$$
(2)
The coefficients in the expansion, of which only the first two are shown, are pseudo tensors of increasing rank. In what follows, we consider only the first term, i.e. we assume that the mean magnetic field in the simulation domain is approximately independent of spatial coordinates. Commonly one separates the $`\alpha `$-tensor into symmetric and antisymmetric parts, i.e.

$$\alpha _{ij}^\mathrm{S}=(\alpha _{ij}+\alpha _{ji})/2;\alpha _{ij}^\mathrm{A}=(\alpha _{ij}\alpha _{ji})/2.$$
(3)
The latter gives rise to an advection term, $`\alpha ^\mathrm{A}\overline{B}=\gamma \times \overline{B}`$, where

$$\gamma =\left(\begin{array}{c}\alpha _{yz}^\mathrm{A}\\ \alpha _{zx}^\mathrm{A}\\ \alpha _{xy}^\mathrm{A}\end{array}\right).$$
(4)
This term describes the pumping of all components of the magnetic field with equal speed, and we shall refer to it as the ‘general pumping effect’. Note that the pumping effect as it is used here refers to the zeroth-order terms in the expansion of the EMF only. Turbulent diffusion ($`\beta `$) and higher order contributions can also give rise to transport phenomena of the mean magnetic field, but they cannot be expressed by advection terms and are to be treated separately. One may also ask whether there are circumstances under which different components of the magnetic field have different advection speeds and directions. A necessary condition that must be fulfilled when addressing this question is that the magnetic field can be decomposed in divergence-free components. Only then can one write down separate equations for the magnetic flux of each field component through a surface moving with a velocity $`w`$, say, thereby allowing one to distinguish advection from other processes. This holds, for example, for a decomposition using toroidal and poloidal components (Kichatinov kichatinov91 (1991)). But unlike Alfvén’s theorem in ideal MHD (where $`w=U`$ for all divergence-free components), it is generally impossible in mean-field electrodynamics to define a moving surface in which flux is conserved. That is due to the presence of the alpha effect and turbulent diffusion. Essentially, what can be achieved is to define a velocity that optimizes flux conservation for a given divergence-free component. For example, in the case of $`\overline{B}_{\mathrm{tor}}`$, one sets $`w_{\mathrm{tor}}=\overline{U}+\gamma ^{(\mathrm{tor})}`$, where $`\gamma ^{(\mathrm{tor})}`$ denotes the pumping velocity of the toroidal field (even though it is itself poloidal).

As a result of the averaging over the horizontal coordinates, our mean magnetic field is of the form $`\overline{B}=[\overline{B}_x(z),\overline{B}_y(z),\overline{B}_z]`$, so that all three Cartesian components are divergence free. If one interprets this field as the mean field of an axisymmetric spherical dynamo, then $`\overline{B}_xe_x+\overline{B}_ze_z`$ corresponds to the poloidal component and $`\overline{B}_ye_y`$ to the toroidal. The meaning of the tensor components $`\alpha _{ij}`$ may be grasped most easily by writing

$$=\alpha ^\mathrm{D}\overline{B}+\gamma ^{(x)}\times \overline{B}_x+\gamma ^{(y)}\times \overline{B}_y+\gamma ^{(z)}\times \overline{B}_z,$$
(5)
where $`\overline{B}_x=\overline{B}_xe_x`$, etc. Here, $`\alpha ^\mathrm{D}`$ denotes the diagonal part of $`\alpha `$, which is responsible for maintaining a large-scale magnetic field in solar-type stars. The off-diagonal terms of the $`\alpha `$ tensor represent various contributions to the pumping effect. As shown by Eq. (5), they are equivalent to an effective advection of the mean magnetic field, where

$`\gamma ^{(x)}`$ $`=`$ $`\gamma +\left(\begin{array}{c}0\\ \alpha _{xz}^\mathrm{S}\\ \alpha _{xy}^\mathrm{S}\end{array}\right),\gamma ^{(y)}=\gamma +\left(\begin{array}{c}\alpha _{yz}^\mathrm{S}\\ 0\\ \alpha _{xy}^\mathrm{S}\end{array}\right),`$
$`\gamma ^{(z)}`$ $`=`$ $`\gamma +\left(\begin{array}{c}\alpha _{yz}^\mathrm{S}\\ \alpha _{xz}^\mathrm{S}\\ 0\end{array}\right)`$ (6)
are the pumping velocities for the respective Cartesian components of the mean field. This representation shows that the pumping velocity can, in principle, be different for each component, as was pointed out by Kichatinov (kichatinov91 (1991)), and that the differences are equal to the off-diagonal, symmetric components of the $`\alpha `$-tensor.

In general, the nine components $`\alpha _{ij}`$ are independent. At the poles however, the rotation axis coincides with the direction of gravity, resulting in isotropy in the horizontal plane. Since the only remaining preferential direction is the vertical one, it follows from elementary symmetry considerations that the pumping velocity cannot have horizontal components, i.e. $`\alpha _{xz}=\alpha _{zx}=\alpha _{yz}=\alpha _{zy}=0`$, and that the pumping effect acts in identical ways on $`\overline{B}_x`$ and $`\overline{B}_y`$, i.e. $`\alpha _{xy}=\alpha _{yx}`$. By the same token $`\alpha _{xx}=\alpha _{yy}`$. Hence at the poles the EMF is given by $`=\alpha _{xx}\overline{B}_\mathrm{H}+\alpha _{zz}\overline{B}_z+\gamma \times \overline{B}`$, where $`\overline{B}_\mathrm{H}=[\overline{B}_x,\overline{B}_y,0]`$ is the horizontal mean field and $`\gamma =\alpha _{xy}e_z`$ is the pumping velocity. At the equator, similar arguments lead to the conclusion that the pumping velocity in the $`x`$-direction vanishes, i.e. $`\alpha _{yz}=\alpha _{zy}=0`$, and that $`\alpha ^\mathrm{D}=0`$. These same components are antisymmetric with respect to the equatorial plane, while the other ones are symmetric.

The remainder of the paper is structured as follows. After introducing the model (Sect. 2), we explain the method by which the dynamo coefficients are calculated (Sect. 3). In Sect. 4 we present results on the pumping effects, their dependence on various parameters, and we make a qualitative comparison with analytical results, where they are available. Briefly, some results on the diagonal alpha coefficients will be mentioned. Finally, the results and their possible implications for stellar dynamos are discussed (Sect. 5).

## 2 Model

For a detailed description of the model and for the full equations and boundary conditions, the reader is referred to Ossendrijver et al. (ossendrijver01 (2001); hereafter cited as Paper I). The simulation domain consists of a rectangular box which is defined on a Cartesian grid, with $`x`$ and $`y`$ denoting the two horizontal coordinates (corresponding to latitude and longitude in spherical coordinates), and $`z`$ denoting depth, respectively (Fig. 1). Gravity is directed in the positive $`z`$-direction. The box rotates about an axis in the $`xz`$-plane, $`\mathrm{\Omega }=\mathrm{\Omega }(e_x\mathrm{sin}\theta +e_z\mathrm{cos}\theta )`$, where $`\theta `$ is the angle between the rotation vector and the $`z`$-axis, or the angular distance to the south pole.

From top to bottom, the simulation domain consists of a thin cooling layer, a convectively unstable layer, and a stably stratified layer with overshooting convection. The consecutive boundaries of the layers are located at depths $`z_1`$, $`z_2`$, $`z_3`$ and $`z_4`$, where $`z_2=0`$ and $`z_3=d`$ in all cases. The horizontal extent of the box, $`L\times L`$, is chosen such that there are always at least a few convective cells within the box. The governing equations are those describing magnetic induction, mass continuity, and the balance of momentum and energy. Solutions are governed by the following independent dimensionless parameters. The Prandtl number and the magnetic Prandtl number are defined as

$$\text{Pr}=\nu /\chi _0,\text{Pm}=\nu /\eta ,$$
(7)
where $`\chi _0=\kappa _2/\gamma _{\mathrm{ad}}\rho _0`$ denotes a reference value of the radiative diffusivity for region 2. In order to avoid possible confusion with the pumping effect we write for the ratio of specific heats $`\gamma _{\mathrm{ad}}=5/3`$. The parameter $`\xi _0`$ determines the internal energy at the top, $`e_0`$, and hence the pressure scale height there, i.e.

$$\xi _0=(\gamma _{\mathrm{ad}}1)e_0/gd.$$
(8)
The initial thermal structure of the three regions is characterized by polytropic indices $`m_i=(1NewA_i)/NewA_i`$, where $`NewA_i=(d\mathrm{ln}T/d\mathrm{ln}p)_i`$ denotes the radiative temperature gradient ($`i=1,2,3`$). The adiabatic temperature gradient is $`NewA_{\mathrm{ad}}=(\gamma _{\mathrm{ad}}1)/\gamma _{\mathrm{ad}}`$. A measure for instability is provided by the (initial) superadiabaticity, $`\delta _i=NewA_iNewA_{\mathrm{ad}}`$, which is positive in an unstably stratified medium. The reference Rayleigh number is defined as

$$\text{Ra}=\frac{d^4g\delta _2}{\nu \chi _0H_{ph}},$$
(9)
where $`H_{ph}=\xi _0d+0.5d/(m_2+1)`$ is the pressure scale height in the middle of the unstable layer, and the parameters $`\delta _2`$ and $`H_{ph}`$ refer to the initial stratification. The Taylor number,

$$\text{Ta}=(2\mathrm{\Omega }d^2/\nu )^2,$$
(10)
measures the importance of rotation relative to viscous dissipation. Finally, $`\sigma _0`$ represents the rate at which internal energy adjusts to the specified reference value $`e_0`$. The remaining parameters are secondary. The Coriolis number, or inverse Rossby number, measures the importance of the Coriolis force:

$$\text{Co}=2\mathrm{\Omega }\tau .$$
(11)
Here $`\tau =d/u_{\mathrm{rms}}`$ is an estimate of the correlation time, and $`u_{\mathrm{rms}}`$ is the rms velocity defined in a suitable way (e.g., by averaging over time and over the bulk of the convection zone in the box). In the initial state, the radiative energy flux, $`F_{\mathrm{rad},z}=\kappa de/dz`$, is assumed to be constant throughout the domain. This determines the radiative conductivities in the three regions according to $`\kappa _i/\kappa _2=(m_i+1)/(m_2+1)`$. In order to obtain a continuous function, $`\kappa `$ is smoothed across thin intermediate layers between the three regions. The initial stratification, obtained from a simplified mixing length model, is calculated with the constraint $`\rho (z_3)=\rho _0`$. Kinetic and magnetic Reynolds numbers are defined as

$$\text{Re}=u_{\mathrm{rms}}d/\nu ,\text{Rm}=u_{\mathrm{rms}}d/\eta .$$
(12)
Dimensionless units are introduced by setting $`d=\rho _0=g=\mu _0=1`$, where $`\mu _0`$ is the vacuum permeability.

## 3 Method

Initially, a simulation is performed without a magnetic field until convection has attained a statistically stationary state. Then at most three numerical experiments are done with imposed weak homogeneous magnetic fields in orthogonal directions. For each run, the EMF and the mean magnetic field are calculated, and by combining the results, the nine coefficients $`\alpha _{ij}`$ are determined through the relation $`=\alpha \overline{B}`$. A necessary condition for the results to be valid is that the statistical properties of the flow are independent of the imposed magnetic field, i.e. that no additional anisotropies are introduced. This is ensured because during the simulation, the field strength remains small compared to the equipartition field strength, so that the Lorentz force is weak.

As was pointed out in Paper I, the determination of dynamo coefficients with sufficient accuracy requires both spatial averaging and temporal averaging over at least several turnover times. Compared to Paper I, we employ an improved method that gives more accurate results than can be obtained by performing an average over the horizontal coordinates plus a simple time average. The method aims to minimize contributions to the EMF (2) due to first and higher order derivatives of the mean magnetic field. Such contributions can arise because with the exception of $`\overline{B}_z`$, mean quantities, defined as horizontal averages, can have a dependence on $`z`$. Thus, a few turnover times after the introduction of a homogeneous magnetic field, the mean field is characterized by rather steep gradients at several locations. The only exception is $`\overline{B}_z`$, which by virtue of the periodic boundary conditions in the horizontal directions must remain unchanged. In the ideal situation where the $`z`$-derivatives of the mean field vanish, the only contributing coefficients are $`\alpha _{ij}`$, and they can be calculated exactly by performing a series of (at most) three runs using imposed magnetic fields with orthogonal orientations. Note that for $`\theta =0\mathrm{°}`$ or $`180\mathrm{°}`$ a pumping effect exists only in the $`z`$-direction; it can be measured from a single run with an imposed field in the $`x`$ or $`y`$ directions. In the presence of significant nonvanishing gradients on the other hand, one could account for the higher-order contributions by doing more experiments using imposed magnetic fields that have gradients and calculating the corresponding dynamo coefficients. Even from one experiment one can calculate the additional diffusive contributions if the mean field is varying in space and time. This allowed Brandenburg & Sokoloff (brandenburg02 (2002)) to estimate all components of the $`\alpha `$ and $`\beta `$ tensors (Eq. 2) from local simulations of accretion discs which feature migratory dynamo cycles.

Instead of following this very cumbersome recipe, we have proceeded by periodically resetting the magnetic field in the entire box to the same initial homogeneous field at times $`t_0+iT`$, where $`i=0,1,\mathrm{},N`$, thereby allowing the field to evolve for only a limited duration $`T`$. In this manner, a sequence of $`N`$ time slices is obtained which, to good approximation, are independent realisations. By averaging a quantity, say $`A`$, over the horizontal coordinates and over the time slices according to

$$\overline{A}(z,t)=\frac{1}{NL^2}\underset{i=0}{\overset{N1}{}}𝑑x𝑑yA(x,y,z,t_0+t+iT),$$
(13)
where $`0t<T`$, one obtains something equivalent to an ensemble average of a time slice. The ‘error’ in $`\overline{A}(z,t)`$ is estimated from the standard deviation, $`\sigma `$, obtained from all the time slices according to $`\sigma _{\mathrm{mean}}=\sigma /\sqrt{N}`$. In order to obtain significant results $`N`$ must exceed a minimal value which depends on several factors such as the magnitude of the coefficient that one tries to measure, the amount of spatial averaging, and the level of turbulence; for the present runs, it was found that $`N`$ should be at least $`5`$.

The most delicate part of the method consists of identifying an optimal subinterval within the mean time slice, say $`t_1tt_2`$, that provides the most reliable estimate of the dynamo coefficients. The best estimate will be defined as the time average over this subinterval. Two conditions must ideally be fulfilled. First, the EMF must have attained the equilibrium value given by Eq. (2). This condition is not trivially fulfilled, because $`=\overline{u\times b}=0`$ at times $`t_0+iT`$, in apparent contradiction with Eq. (2). Obviously the paradox arises because Eq. (2) is not valid instantaneously from the moment when a homogeneous magnetic field is introduced, since it takes a while for the magnetic fluctuations to develop. Thus $`t_1`$ must be chosen sufficiently late. Secondly, $`t_2`$ must be chosen such that the gradients of the mean field are still numerically small compared to the imposed field itself divided by, say, the thickness of the unstable layer. If this condition is satisfied, then $`\alpha _{ij}`$ is the dominant contribution to the EMF, provided that $`\alpha _{ij}`$ does not vanish for reasons of symmetry or other. Both conditions can be met if the EMF evolves sufficiently more rapidly than does the mean magnetic field. In practice, the time scales are rather similar in the present simulations, and both conditions are at best marginally fulfilled.

The method is illustrated in Figs. 2 and 3 for run A45X (see Table 1 for the parameters of this run). Fig. 2 shows the $`z`$-component of the EMF at a fixed depth in the box, defined as an average over $`x`$ and $`y`$, as a function of time. Although the temporal evolution varies strongly between the different initialisations, the additional averaging over the subinterval $`[t_1,t_2]`$ and over the subsequent realisations leads to a well-converged result. This may be compared with Fig. 2 of Cattaneo & Hughes (cattaneo96 (1996); henceforth referred to as CH96), which suggests that, for a weak initital magnetic field, convergence is achieved only after more than about $`50`$ turnover times. However, there are important differences compared with their calculation. First, the dynamics in CH96 are very different, because the flow results from external forcing rather than moderately turbulent convection, and because the Lorentz force dominates, whereas it is negligible in our case. We have smaller fluctuations in the horizontal averages, partly because the length scale of convection is always at least a few times smaller than the box size, whereas they are similar in CH96. Also, Re is smaller in comparison with CH96, and therefore the level of fluctuations. Note that in order to compare Re in Table 1 with the Reynolds number, $`U/\nu k_\mathrm{f}=100`$, of CH96, where $`k_\mathrm{f}`$ is the wave number of the forcing, Re should be multiplied with a factor of about $`1/2\pi `$. Secondly, the averaging in CH96 starts from an initial state where there is no alpha effect, which exists only after $`30`$ turnover times. This introduces a transient in the cumulative time average, which is absent in our case. As a result of these differences, we obtain convergence already if $`N6`$, even though this corresponds to a total averaging time, $`N(t_2t_1)`$, of a few turnover times only. Fig. 3, top, shows mean time slices for the $`z`$-component of the EMF at different depths in the box. After an initial phase of approximately linear growth which, depending on the depth, lasts about $`48`$ units, the EMF levels off to the equilibrium value. After that, it continues to evolve more slowly as it adapts to the changing mean magnetic field in accordance with Eq. (2). The lower figure shows $`\overline{B}_x`$ as a function of depth at consecutive moments. Starting from the initial condition, the mean magnetic field evolves while developing increasing gradients, especially near the upper and lower boundaries of the convective layer. In the time interval $`[5,10]`$ the magnitude of the gradients in the bulk of the unstable layer is $`d/dz\mathrm{ln}\overline{B}_x0.5`$, and in the interval $`[20,25]`$ it is about $`1`$. It follows that the best available time interval for calculating dynamo coefficients from run A45X is roughly between $`5`$ and $`10`$.

As soon as the EMF levels off to its equilibrium value (2), the mean magnetic field already begins to develop considerable gradients at some depths. Nevertheless, the results are less contaminated with higher-order contributions than those obtained by performing a simple time average accross a long simulation in which the magnetic field was initialized only once. An indication of the error in the coefficients is provided by the variation of the coefficients within the optimal time interval as shown by the shading in Figs. 48. In order to assess the possible difference between the results of the improved method and those of the less accurate method of Paper I, we have in a few cases also calculated $`\alpha _{ij}`$ using the time interval $`[t_1,t_2]=[20,25]`$. Although this choice of $`[t_1,t_2]`$ yields the correct sign and order of magnitude, the higher-order contributions can amount to $`40\%`$ in some cases. All results presented below were therefore derived using the earliest possible interval, so that the influence of the higher-order contributions is minimized. As is shown in Fig. 3, bottom, the gradients of the mean field can be particularly strong near the top and the bottom of the unstable layer, so that at these locations the results should be treated with some caution. To avoid misunderstandings, we emphasize that the accuracy discussed here is distinct from the error introduced if the time sequence is too short to average out the fluctuations.

## 4 Results

The numerical simulations are grouped together in several series named A–E, each of which may contain several experiments performed at different latitudes and having different orientations of the imposed field. Each run is referred to by a code such as A45X, which identifies it as belonging to series A, having $`\theta =45\mathrm{°}`$ and an imposed magnetic field in the $`x`$-direction. The parameters of each series are summarized in Table 1. They may be characterized as follows: series A is moderately turbulent and has a rotation rate appropriate for the base of the solar convection zone. Series B has the same parameters as A except for rotation, which is slower. Series C has the highest Rayleigh number and the most turbulent convection. Series D has a Rayleigh number just above marginal stability. Accordingly, the flow is laminar, and the Reynolds number is small. It has the slowest rotation of all. Series E is identical to D, except for the rotation rate. In series A, $`\theta `$ was varied between $`0\mathrm{°}`$ and $`90\mathrm{°}`$ with a stepsize of $`15\mathrm{°}`$; all other simulations were done for $`\theta =45\mathrm{°}`$. The optimal time interval, $`[t_1,t_2]`$, was determined to be $`[5,10]`$ for series A and B, $`[4,7]`$ for C, and $`[15,20]`$ for D and E. The value of $`N`$ was in the range between $`6`$ and $`21`$. Although the total amount of time used in the averaging, $`N(t_2t_1)`$, corresponds to several turnover times only, the combined spatial and temporal averaging as defined by Eq. (13) was in most cases sufficient to yield significant results. This is illustrated by the convergence of the cumulative average shown in Fig. 2, by the well-defined shapes of the curves of $`_z`$ and $`\overline{B}_x`$ in Fig. 3, and by the narrowness of most of the shaded areas in Figs. 48.

In the following subsections we present the main results of the simulations. As will become apparent, many of the results are new in the sense that apart from the pumping effect in the vertical direction, none of the effects have been numerically investigated before as far as we know. Also, this is the first numerical investigation in which pumping effects are quantified in terms of gamma coefficients. For that reason, and in order to compare the various pumping effects, it is worthwile to consider also the vertical pumping effect, although it has already received much recent attention in the literature. Results for the coefficients $`\alpha _{ij}`$ of all series are shown in Figs. 48 in the form of the three Cartesian components of $`\gamma `$ and the independent elements of $`\alpha _{ij}^\mathrm{S}`$. From these, all components of $`\alpha _{ij}`$ can, if necessary, be derived using Eqs (3) and (6). The black curves in Figs. 48 are averages of a given coefficient over the time interval $`[t_1,t_2]`$. The shaded areas indicate the range of values attained at individual times within the interval $`[t_1,t_2]`$, and are meant to give an error indication. Some confirmation that the error estimates are consistent, provided that they are used conservatively, can be seen in series A90, where $`\gamma _x,\alpha _{yz}^\mathrm{S},\alpha _{xx},\alpha _{yy}`$ and $`\alpha _{zz}`$ are known to vanish on symmetry grounds (Figs. 4, 6 and 7). For these coefficients the shaded area is not always consistent with zero. This departure gives an additional estimate of the expected error. We reiterate that results for $`\alpha _{ij}`$ near the top and the bottom of the unstable layer should be treated with caution because of the enhanced possibility of contamination with higher-order terms.

### 4.1 Latitudinal pumping effect

In the latitudinal direction, the simulations yield a weak net general advection, measured by $`\gamma _x`$, which is equatorward ($`\gamma _x>0`$) in the bulk of the unstable layer and can be poleward near the upper surface and in the overshoot layer (Figs. 45). At the pole and at the equator the effect vanishes; it is strongest near $`\theta =60\mathrm{°}`$. The symmetric contribution $`\alpha _{yz}^\mathrm{S}`$ is about as weak as $`\gamma _x`$, but shows little systematic behaviour as a function of $`\theta `$ (Fig. 6). Both increase with increasing rotation rate, since they are rotationally induced (compare D45 and E45). As can be inferred from the relations $`\gamma _x^{(y)}=\gamma _x+\alpha _{yz}^\mathrm{S}`$ and $`\gamma _x^{(z)}=\gamma _x\alpha _{yz}^\mathrm{S}`$ (Eq. 6), the $`y`$ and $`z`$ components of the mean field have different latitudinal velocities (Fig. 9). In all simulations, $`\overline{B}_y`$ tends to be pumped towards the equator ($`\gamma _x^{(y)}>0`$) in the bulk of the unstable layer. The situation with regard to the component $`\overline{B}_z`$ is less clear. It is generally pumped towards the pole ($`\gamma _x^{(z)}<0`$) in the lower part of the unstable layer, but the behaviour in the bulk of the unstable layer appears to vary with latitude. The average over the unstable layer turns out to be mostly positive (equatorward) for both components (Fig. 9). The lowest-order contribution to latitudinal pumping in the expansion of the EMF presented by, e.g., Krause & Rädler (krause80 (1980)), p. 193, can be written as $`\mathrm{\Delta }\alpha _{ij}=A\mathrm{\Omega }_iϵ_{jkl}g_k\mathrm{\Omega }_l+B\mathrm{\Omega }_jϵ_{ikl}g_k\mathrm{\Omega }_l`$, where $`A`$ and $`B`$ are true scalars, i.e. even functions of $`\mathrm{cos}\theta `$. Thus, the latitudinal pumping is at least of second order in $`\mathrm{\Omega }`$. One mechanism that may explain the effect is provided by the rotational anisotropies studied by Kichatinov (kichatinov91 (1991)). Indeed, the general tendency of his result is that the toroidal magnetic field, which for an axisymmetric mean field corresponds to our $`\overline{B}_y`$, is pumped towards the equator, while the poloidal field (i.e. $`\overline{B}_x`$ and $`\overline{B}_z`$) is pumped towards the poles (cf. his Eq. 3.26).

### 4.2 Longitudinal pumping effect

Surprisingly, the longitudinal pumping effect had so far not been identified in simulations of magnetoconvection, although it is often the strongest pumping effect except if the box is situated near the pole (Figs. 45, column 2). The predominantly negative sign of $`\gamma _y`$ signifies that the mean field is advected in the retrograde direction within the bulk of the unstable layer and the overshoot layer. In most cases, there is a thin layer near the top of the box where the field is pumped in the prograde direction. Unlike the vertical pumping effect (see below), here the symmetric part of $`\alpha `$-tensor dominates. From the relations $`\gamma _y^{(x)}=\gamma _y\alpha _{xz}^\mathrm{S}`$ and $`\gamma _y^{(z)}=\gamma _y+\alpha _{xz}^\mathrm{S}`$ (Eq. 6) it follows that the longitudinal pumping velocities of $`\overline{B}_x`$ and $`\overline{B}_z`$ can be opposite. Since $`\alpha _{xz}^\mathrm{S}`$ is positive in the bulk of the unstable layer (Fig. 6, column 3), one finds that $`\overline{B}_x`$ is pumped mainly in the retrograde direction, but $`\overline{B}_z`$ mainly progradely near the top of the unstable layer, while retrogradely near the bottom. This behaviour is reflected in the volume averages (Fig. 9). The longitudinal pumping effect is strongly dependent on latitude; it vanishes at the pole and peaks at the equator. The lowest-order term which contributes to longitudinal pumping (Krause & Rädler krause80 (1980), p. 193) is $`\mathrm{\Delta }\alpha _{ij}=Ag_i\mathrm{\Omega }_j+Bg_j\mathrm{\Omega }_i`$, where $`A`$ and $`B`$ are true scalars. Clearly, the effect requires both gravity and rotation. The antisymmetric part, corresponding to general longitudinal pumping, has been computed by Rüdiger & Kichatinov (ruediger93 (1993)) in the form of $`\alpha _3(\mathrm{\Omega }_iG_j\mathrm{\Omega }_jG_i)`$, where the coefficient $`\alpha _3`$ proved to be positive-definite. Here $`G`$ is the gradient vector of the turbulence intensity, i.e. $`G=\mathrm{ln}u_{\mathrm{rms}}`$. One obtains

$$\gamma _y\mathrm{\Omega }\mathrm{sin}\theta \frac{d\overline{u^2}}{dz},$$
(14)
while for the vertical pumping the simple relation

$$\gamma _z\frac{d\overline{u^2}}{dz}$$
(15)
results, with only a very weak dependence on the Coriolis number. Thus $`\gamma _y`$ and $`\gamma _z`$ should always satisfy the relation $`\gamma _y\gamma _z<0`$, which is indeed the case in all the presented simulations. According to Eqs (14) and (15), the nodes which exist in all profiles of $`\gamma _y`$ and $`\gamma _z`$ should be very close together, which is really true in Figs. 45. For completeness we also note that longitudinal pumping cannot be attributed to the type of rotational anisotropies studied by Kichatinov (kichatinov91 (1991)), since the resulting pumping effects exist only in the radial and latitudinal directions.

### 4.3 Vertical pumping effect

The main features of the vertical pumping effect are as follows. The direction of vertical pumping is downward ($`\gamma _z>0`$) in the bulk of the unstable layer and in the overshoot layer below it (Figs. 45). Near the top of the box there is a thin layer where the pumping is directed upwards ($`\gamma _z<0`$). The strongest effect is observed in the laminar cases (D45 and E45), where upward and downward speeds can amount to about $`15\%`$ of the rms convective velocity (Table 1). In the more turbulent runs, the maximum speeds are smaller. Since the top layer where the pumping is upward is less pronounced in the turbulent runs, the decrease is less apparent in the volume-averaged coefficient (Fig. 9) because there is less cancellation. There is little dependence on rotation (compare A45 with B45 and D45 with E45) or latitude (Figs. 4 and 9). Also the rms convective velocity has a depth dependence that is similar in all runs, namely a plateau at the top of the domain and below that a monotonic decrease with depth. This in combination with the robustness of the vertical pumping, in particular its virtual independence of rotation, suggests an explanation in terms of ’turbulent diamagnetism’, i.e. by transport down the gradient of the rms convective velocity (Eq. 15). Analytical expressions for the diamagnetic pumping effect derived under the FOSA assumption were given by Kichatinov & Rüdiger (kichatinov92 (1992)).

If the Rayleigh number is increased from $`10^6`$ (series B) to $`410^6`$ (series C), so that the convection becomes more strongly turbulent, $`|\gamma _z|`$ decreases. This would appear to be at odds with the result of Tobias et al. (tobias01 (2001)), who found more efficient pumping for increased Rayleigh numbers. But they measure the pumping effect through the ratio of magnetic fluxes in the unstable layer and the overshoot layer. This ratio is determined both by transport properties and storage properties of the magnetic field, i.e. by the balance between pumping, turbulent diffusion in the convection zone, and molecular diffusion in the overshoot layer ($`\eta `$). If $`\eta `$ is decreased, flux concentrations decay more slowly, so that the storage capacity is increased. Since the increase of the Rayleigh number in Tobias et al. (tobias01 (2001)) is accompanied by a decrease of $`\eta `$, the flux ratio increases, even if the pumping is somewhat less efficient.

In series A it can be seen that the symmetric contribution $`\alpha _{xy}^\mathrm{S}`$ is positive and that it increases if $`\theta `$ increases from $`15\mathrm{°}`$ to $`75\mathrm{°}`$ (Fig. 6). Using the relations $`\gamma _z^{(x)}=\gamma _z+\alpha _{xy}^\mathrm{S}`$ and $`\gamma _z^{(y)}=\gamma _z\alpha _{xy}^\mathrm{S}`$ (Eq. 6) one can infer from this that the pumping speeds for $`\overline{B}_x`$ and $`\overline{B}_y`$ deviate increasingly from one another: the downward pumping of $`\overline{B}_x`$ becomes faster, and the top layer where $`\overline{B}_y`$ is pumped upward becomes deeper; this is reflected in Fig. 9. A comparison of runs B45 with A45 and D45 with E45 shows that $`\alpha _{xy}^\mathrm{S}`$ increases with increasing rotation rate, while $`\gamma _z`$ is hardly affected (Figs. 46). In run D45, which has the weakest rotation of all, the value of $`\alpha _{xy}^\mathrm{S}`$ is not above the noise level, and the same is probably the case for run B45. From the symmetry considerations of Krause & Rädler (krause80 (1980)), p. 193, it is inferred that the lowest-order contribution to $`\gamma _z`$ stems from the term $`\mathrm{\Delta }\alpha _{ij}=Aϵ_{ikj}g_k`$, while that to $`\alpha _{xy}^\mathrm{S}`$ stems from $`\mathrm{\Delta }\alpha _{ij}=A\mathrm{\Omega }_iϵ_{jkl}g_k\mathrm{\Omega }_l+B\mathrm{\Omega }_jϵ_{ikl}g_k\mathrm{\Omega }_l`$, where $`A`$ and $`B`$ are true scalars. The latter term also gives rise to latitudinal pumping, as mentioned above. Unlike $`\gamma _z`$, the coefficient $`\alpha _{xy}^\mathrm{S}`$ is rotationally induced, and this is confirmed by the simulations.

### 4.4 Alpha effect

As a byproduct of the determination of the pumping effects, we obtain results for the diagonal part of the $`\alpha `$-tensor. These results are an extension to what was reported in Paper I, where all simulations were done for $`\text{Ra}=\mathrm{50\hspace{0.17em}000}`$ and $`\theta =0\mathrm{°}`$. The sign and magnitude of $`\alpha _{yy}`$ are of special interest, since this component is responsible for the generation of the large-scale magnetic field in solar-type stars. The depth dependence of $`\alpha _{xx}`$ and $`\alpha _{yy}`$ has a common shape in all series, provided there is sufficient rotation, namely a negative sign in the bulk of the convection zone, and a positive sign in the overshooting layer (Figs. 78). These features, which are expected for the southern hemisphere, are explained, e.g., in Paper I. Series A (Fig. 7) provides detailed results on the latitude dependence. The amplitudes of $`\alpha _{xx}`$ and $`\alpha _{yy}`$ increase with angular distance from the equator up to a point close to the south pole, more or less consistent with the commonly assumed $`\mathrm{cos}\theta `$-function (Fig. 9). At the pole, there appears to be an unexpected dip. The amplitudes of $`\alpha _{xx}`$ and $`\alpha _{yy}`$ are approximately the same near the pole – for $`\theta =0\mathrm{°}`$ they are equal on symmetry grounds – but closer to the equator $`\alpha _{yy}`$ is somewhat larger (especially E45; Fig. 8). Without rotation, the diagonal $`\alpha `$-terms vanish (the results in Fig. 8 for D45, which has the weakest rotation, are probably not significantly different from zero), and the increase of $`\alpha _{xx}`$ and $`\alpha _{yy}`$ with rotation (compare B45 with A45 and D45 with E45 in Fig. 8) confirms the results of Paper I. The same can be said for $`\alpha _{zz}`$, which for $`\text{Co}2`$ has the opposite sign of $`\alpha _{xx}`$ (tentatively visible for B45 in Fig. 8 but clear in Paper I). If rotation increases beyond a certain point, $`\alpha _{zz}`$ as a function of depth developes multiple sign changes (cf. series A), and its amplitude falls behind that of $`\alpha _{xx}`$, even more so if one also averages over depth (Fig. 9). This is the rotational quenching of the vertical alpha effect reported in Paper I.

## 5 Discussion

We have performed numerical simulations of magnetoconvection in a Cartesian box with a rotating frame of reference which is meant to represent a section from the lower part of a stellar convection zone. Through an inversion of the relation between the mean magnetic field and the electromotive force, we have measured pumping effects in all three directions as well as the symmetric elements of the alpha tensor. For the first time, the existence of non-radial pumping effects is hereby established numerically, and quantified in terms of gamma coefficients. Previously, their existence was known only on the basis of symmetry considerations in the framework of mean-field electrodynamics (Krause & Rädler krause80 (1980), p. 193), and analytical results based on the first-order smoothing approximation (Kichatinov kichatinov91 (1991); Rüdiger & Kichatinov ruediger93 (1993)). Such analytical results are available only for some of the coefficients that appear in the expansion of the EMF given by Krause & Rädler (krause80 (1980)). Moreover, the validity of results based on FOSA is questionable for real and simulated stellar convection zones alike, and they can be readily evaluated only for simple turbulence models but not for actual convective flows. For these reasons, it is difficult at present to explain all of the phenomena that are found in the simulations, although we did observe that several qualitative features of the expressions derived using FOSA are confirmed by the simulations, as was the case in Paper I.

The Rayleigh number was varied between $`10^4`$, corresponding to laminar convection, and $`410^6`$, corresponding to mildly turbulent convection. Within this parameter range, we found remarkably consistent behaviour of the dynamo coefficients. Clearly, the regime of stellar convection, which is characterized by Rayleigh numbers many orders of magnitude larger, remains far out of reach because of the limited spatial resolution. The main resulting difference between the simulations and stellar convection is as follows. In the lower part of stellar convection zones, the superadiabaticity is very small (typically $`\delta 10^6`$) and, according to the standard mixing-length expressions (e.g., Stix stix89 (1989), p. 198), this results in velocities of the order of $`50`$ m s<sup>-1</sup>, or $`\text{Ma}10^4`$. In the simulations, convection is less turbulent. After the onset of convection, the superadiabaticity in the unstable layer decreases from its initial value, $`\delta _2=0.1`$, to a value of at least about $`10^2`$ (series C). The convective velocities are therefore much larger than in the Sun, if expressed in m s<sup>-1</sup> or as a Mach number. We may illustrate this for series A, where the unstable layer covers $`3.7`$ pressure scale heights (Table 1). Counting upwards from the bottom of the solar convection zone, this corresponds to about $`\mathrm{90\hspace{0.17em}000}`$ km. Setting $`g540`$ m s<sup>-2</sup>, one would obtain $`\sqrt{dg}210^5`$ m s<sup>-1</sup> for the unit of velocity. Thus the rms convective velocity of series A (Table 1) would correspond to about $`210^4`$ m s<sup>-1</sup>, i.e. approximately $`400`$ times larger than the solar value. About the same ratio holds for the Mach numbers. Only simulations with a much higher resolution and much smaller entropy gradients can improve this ratio. However, the rotation rate was chosen such that the Coriolis number, i.e. the relative importance of the Coriolis force in the equation of motion, is within the appropriate range. Also, it is to be expected that beyond some Rayleigh number, a further decrease of the molecular viscosity or the magnetic diffusivity affects only the smallest scales while having little or no effect on the large scale. From the consistency of the results one may therefore suspect that the current simulations already provide some clues about the mean-field transport properties of stellar convection, if properly scaled.

In order to assess the astrophysical relevance of the results, the pumping speeds should therefore be compared with the rms convective velocity for the unstable layer (Table 1) rather than be directly expressed in m s<sup>-1</sup>. The latitudinal pumping effect for the toroidal mean field ($`\overline{B}_y`$) was found to be mainly equatorward with speeds up to about $`0.1u_{\mathrm{rms}}`$ in the turbulent runs for $`\theta =45\mathrm{°}`$, while the volume averages ($`\gamma _x^{(y)}=\gamma _x+\alpha _{yz}^\mathrm{S}`$ in Fig. 9) are less than about $`0.03u_{\mathrm{rms}}`$. The equatorward motion in the solar butterfly diagram corresponds to a speed of about $`1`$ m s<sup>-1</sup>, so that one may conclude that the latitudinal pumping effect of the toroidal mean field has the same order of magnitude. We therefore point out the possibility that the equatorward motion of solar magnetic activity is not the result of a dynamo wave or meridional bulk motion, but of latitudinal pumping of the toroidal mean magnetic field. Results for the latitudinal pumping of the radial mean field do not show consistent behaviour; the poleward motion at the top of the box found in some runs amounts to at most $`0.1u_{\mathrm{rms}}`$, but this should be treated with caution because of contamination with higher-order terms due to the strong gradients in the mean field near the top of the box. Since the simulations are not meant to represent the photosphere, the relevance of this result for the observed poleward motion of poloidal fields at the solar surface is in any case hard to establish. We merely point out the possibility that the observed motion is the result of pumping of the deep-seated global poloidal field rather than of meridional bulk surface motion. The longitudinal pumping attains maximal speeds of about $`0.2u_{\mathrm{rms}}`$ in the retrograde direction. Again, the mainly prograde pumping of the poloidal mean field near the top of the box should be treated with caution because of the possibility of strong higher-order contributions. Perhaps, if carried over to the surface of the Sun, it could explain the observed motion of the active longitudes in the solar wind relative to solar surface rotation (Ruzmaikin et al. ruzmaikin01 (2001)). We also note that preliminary calculations of a mean-field model indicate that the inclusion of longitudinal pumping in an $`\alpha ^2`$-type dynamo can lead to oscillating solutions, which are rare for such models. This can be relevant for explaining magnetic cycles in stars which do not have differential rotation (Rüdiger et al. ruediger02a (2002)). The downward pumping in the turbulent runs attains speeds of up to about $`0.1u_{\mathrm{rms}}`$, and the volume average about half of that. In the lower part of the solar convection zone this translates to a few m s<sup>-1</sup>. Such speeds may be sufficient to overcome the magnetic buoyancy of flux generated in the convection zone and transport the flux downward into the overshoot layer. This layer coincides more or less with the tachocline, a region where the rotation rate changes rather suddenly with depth. Through stretching of field lines, the shearing motions can generate from poloidal field components a strong field of predominantly toroidal orientation, as is required in order to account for various features of the large-scale solar magnetic field.

The origin of the poloidal magnetic field and its locus of generation still remain an open question. According to one likely scenario, the interface dynamo (Parker parker93 (1993), Ossendrijver & Hoyng ossendrijver97 (1997), Charbonneau & MacGregor charbonneau97 (1997)), the poloidal field is generated by an alpha effect in the lower half of the convection zone. Our results for the coefficient $`\alpha _{yy}`$, which corresponds to $`\alpha _{\varphi \varphi }`$ in spherical coordinates, are of particular interest for interface-type stellar dynamo models. We should stress however that it cannot be decided from the present simulations whether the large-scale magnetic field of solar-type stars is generated by the convective alpha effect investigated here or by a different alpha effect such as that resulting from the buoyancy instability of magnetic flux tubes in the overshoot layer. We confirm that the convective alpha effect is strongest near the poles. For the Sun, this should, in combination with the differential rotation near the poles, as inferred from helioseismic inversions, result in strong dynamo action near the poles (Rüdiger & Brandenburg ruediger95 (1995)), which is not observed – sunspots are concentrated in belts that extend to a latitude of at most $`40\mathrm{°}`$ on either side of the equator. This well-known problem may perhaps be resolved by taking into account meridional circulation, and by noting that the buoyancy instability of toroidal magnetic flux tubes anchored in the overshoot layer is reduced or suppressed at high latitudes (Ferriz-Mas & Schüssler ferrizmas95 (1995)). The sign of $`\alpha _{yy}`$ confirms what is generally found in MHD simulations (e.g. Paper I), namely a negative sign in the bulk of the convection zone on the southern hemisphere, i.e. a positive sign on the northern hemisphere, and a sign change at the bottom of the convection zone. This distribution of signs puts a tough constraint on theories of the solar dynamo. Generally, mean-field dynamo models that incorporate solar differential rotation, in particular a positive $`d\mathrm{\Omega }/dr`$ near the equator near the base of the convection zone, produce a butterfly diagram with the required equatorward migration of the magnetic activity belts only in combination with a negative $`\alpha _{\varphi \varphi }`$ on the northern hemisphere. There are various ways in which this dilemma could be resolved; perhaps the relevant dynamo effect is located at the very bottom of the convection zone where the sign has changed. If $`\alpha _{\varphi \varphi }`$ is positive on the northern hemisphere, equatorward migration can be obtained if a suitable meridional circulation is included (Choudhuri et al. choudhuri95 (1995); Durney durney96 (1996); Küker et al. kueker01 (2001); Rüdiger et al. ruediger01 (2001); Rüdiger & Arlt ruediger02 (2002)). As mentioned above, the present results suggest the alternative possibility of latitudinal pumping. This suggestion needs however to be verified in mean-field models that take into account all relevant dynamo coefficients as well as solar differential rotation.

###### Acknowledgements.

M.O. acknowledges financial support by the Deutsche Forschungsgemeinschaft. The calculations were carried out largely on the Origin 2000 of the Rechenzentrum of the University of Freiburg i. Br..


text:
# 1 Introduction

## 1 Introduction

Separation of variables for the Hamilton–Jacobi equation is a classical topic that is still very much investigated and has connections with several important research fields, such as stationary reductions of soliton equations (see  and references therein), algebraic completely integrable systems , Riemannian geometry , Bäcklund transformations and Baxter’s $`Q`$-operator , and string theory .

The aim of this paper is to present, in the simple example of the Neumann system, the main ideas of a new approach  to separation of variables. This point of view is based on the geometry of bi-Hamiltonian manifolds, and has been successfully applied to the stationary reductions of the KdV hierarchy  and to Toda lattices .

The Neumann system is a well-known and very much studied mechanical system (see, e.g., ), given by a (mass 1) particle moving on the (unit) sphere under the influence of a quadratic potential $`V(x,y,z)=\frac{1}{2}(a_1x^2+a_2y^2+a_3z^3)`$, where $`a_1<a_2<a_3`$ and $`(x,y,z)`$ are Cartesian coordinates whose origin coincides with the center of the sphere. In 1859, Carl Neumann showed that the Hamilton–Jacobi equation of this system is additively separable in the so-called spheroconical (or elliptical spherical) coordinates $`(\lambda _1,\lambda _2)`$, given by

$$\frac{x^2}{\lambda a_1}+\frac{y^2}{\lambda a_2}+\frac{z^2}{\lambda a_3}=\frac{(\lambda \lambda _1)(\lambda \lambda _2)}{(\lambda a_1)(\lambda a_2)(\lambda a_3)},$$

with $`a_1<\lambda _1<a_2<\lambda _2<a_3`$. Once they have been introduced, one can check that they are separation coordinates, but the problem we want to address is to find out these coordinates in a systematic way. It is known that the separation variables of the Neumann system can be supplied by a Lax representation, but then the problem is to find such a representation for a given system, which is in general a quite difficult task. We will show that a careful study of the Neumann system will allow us to deduce in a natural way the spheroconical coordinates as separation variables.

The plan of the paper is the following. In Section 2 we will introduce the idea (due to A. Nijenhuis) of the geometrization of a coordinate system by means of a suitable tensor field $`L`$. This allows one to find a kind of compatibility condition between $`L`$ and a given Hamiltonian $`H`$, ensuring that the coordinates induced by $`L`$ separate the Hamilton–Jacobi equation associated with $`H`$. This condition is an intrinsic (i.e., coordinate-free) form of the classical Levi-Civita separability conditions. Before applying this results to the Neumann system in Section 4, we will comment on their meaning in the theory of bi-Hamiltonian manifolds in Section 3. Section 5 is devoted to two additional results on the Neumann system, that can be easily obtained from our standpoint (using some facts highlighted in ): A simple rule for the construction of the integrals of motion and the extension of the bi-Hamiltonian structure from the cotangent bundle $`T^{}S^2`$ of the sphere to a 5-dimensional manifold, in order to give the Neumann system a bi-Hamiltonian formulation, which is missing in $`T^{}S^2`$. The final section contains some concluding remarks.

Acknowledgments. The results presented in this paper have been obtained in collaboration with Gregorio Falqui and Franco Magri, to whom I am very grateful. I wish to thank also Claudio Bartocci, Sergio Benenti, and Giovanni Rastelli for useful discussions, and the organizers of NEEDS 2001 for the opportunity to present these results there and for the nice atmosphere at the conference. This work has been partially supported by the Italian M.I.U.R. under the research project Geometry of Integrable Systems.

## 2 The search for separation coordinates

Let $`𝒬`$ be an $`n`$-dimensional manifold, that can be thought of as the configuration space of a mechanical system (the sphere, for the Neumann system), and let $`H`$ be a function on $`T^{}𝒬`$ (the Hamiltonian). In this section we will present a method to look for the variables in which the Hamilton–Jacobi equation for $`H`$ separates. This strategy will prove to be quite efficient in the case of the Neumann system.

The first important point, due to Nijenhuis , is the idea of the geometrization of a coordinate system. If $`\{q^i\}`$ are local coordinates on an open subset $`U𝒬`$, then the tensor field $`L`$ of type $`(1,1)`$ defined by

$$L\frac{}{q^i}=q^i\frac{}{q^i}$$

has vanishing Nijenhuis torsion, i.e.,

$$[LX,LY]L[LX,Y]L[X,LY]+L^2[X,Y]=0$$

for every pair $`(X,Y)`$ of vector fields on $`U`$. Viceversa, a tensor field $`L`$ of type $`(1,1)`$ whose torsion is zero induces local coordinates in a neighborhood of any point where the eigenvalues $`(\lambda _1,\mathrm{},\lambda _n)`$ of $`L`$ are distinct. In many cases these eigenvalues are functionally independent, so that

$$L\frac{}{\lambda _i}=\lambda _i\frac{}{\lambda _i},$$

i.e., they can be chosen as the coordinates associated with $`L`$. In the following we will suppose to be in this situation, that is, we will look for a tensor field $`L`$ with vanishing torsion and functionally independent eigenvalues.

Once we have replaced the coordinate system with the geometric object $`L`$, it is quite natural to look for a “compatibility condition” between $`L`$ and a function $`HC^{\mathrm{}}(T^{}𝒬)`$ entailing that the coordinates given by $`L`$ are separation variables for $`H`$. This condition can be obtained by lifting $`L`$ from $`𝒬`$ to $`T^{}𝒬`$ by means of a procedure called “complete lifting” . It gives rise to a torsionless tensor field $`N`$ of type $`(1,1)`$ on $`T^{}𝒬`$, described in fibered coordinates $`(q^i,p_i)`$ as

$$N\frac{}{q^i}=L_i^j\frac{}{q^j}+\left(\frac{L_j^k}{q^i}\frac{L_i^k}{q^j}\right)p_k\frac{}{p_j},N\frac{}{p_i}=L_j^i\frac{}{p_j}.$$

Now that we have $`H`$ and $`N`$ on the same manifold $`T^{}𝒬`$, we can obtain the separability condition we are looking for as follows. First, we use $`H`$ and $`N`$ to construct the 2-form $`\omega _H=d\left(N^{}dH\right)`$, where $`N^{}`$ is the adjoint of $`N`$. Then we consider the Hamiltonian vector field $`X_H`$ and the vector fields $`NX_H`$, $`N^2X_H`$, …obtained by iteration, and we suppose that the distribution $`D_H`$ spanned by $`X_H`$, $`NX_H`$,…, $`N^{n1}X_H`$ is $`n`$-dimensional (i.e., that these vector fields are linearly independent at any point).

###### Theorem 1

In the above-mentioned hypotheses, the coordinates associated with $`L`$ are separation variables for $`H`$ if and only if the 2-form $`\omega _H`$ annihilates the distribution $`D_H`$:

$$\omega _H|_{D_H}=0.$$
(1)
The proof consists in writing the conditions

$$\omega _H(N^kX_H,N^lX_H)=0,\text{for }k,l=0,\mathrm{},n1\text{,}$$
(2)
in the canonical coordinates $`(\lambda _1,\mathrm{},\lambda _n,\mu _1,\mathrm{},\mu _n)`$, where the $`\lambda _i`$ are the eigenvalues of $`L`$ and the $`\mu _i`$ are their conjugate momenta. Then, taking into account that $`N`$ is diagonal in these coordinates, it is not difficult to realize that the equations (2) are equivalent to the classical Levi-Civita separability conditions (, p. 208).

Therefore, we have found the separability condition between $`L`$ and $`H`$. It is simply the vanishing of a suitable 2-form on the distribution generated from $`X_H`$ by means of the complete lift $`N`$ of $`L`$. We stress that this condition is a concise and, above all, intrinsic form of the Levi-Civita equations. This means that one can check the separability of a given Hamiltonian $`H`$ in the coordinates associated with a tensor field $`L`$ before computing these coordinates. Moreover, one can impose condition (1) on $`L`$ to search for separation coordinates for $`H`$. In Section 4 we will show how to exploit this fact in order to systematically deduce the separability of the Neumann system.

## 3 The bi-Hamiltonian meaning of the separability condition

Before applying the results of the previous section to the Neumann system, let us make some comments on their geometrical meaning.

1. Let $`𝒬`$, $`L`$, and $`N`$ be as in Section 2. Then the cotangent bundle of $`𝒬`$ is a bi-Hamiltonian manifold . The first Poisson bracket is the canonical one,

$$\{F,G\}=\omega (X_F,X_G),$$

while the second one is given by

$$\{F,G\}^{}=\omega (NX_F,X_G).$$

If $`P`$ (resp. $`P^{}`$) is the Poisson tensor of $`\{,\}`$ (resp. $`\{,\}^{}`$), this means that $`P^{}=NP`$.

2. Assume that $`H`$ fulfills the separability conditions (1). Then the distribution $`D_H`$ is integrable, so that there (locally) exist $`n`$ independent functions $`H_1,\mathrm{},H_n`$ which are constant on the leaves of $`D_H`$. They are (local) first integrals for $`X_H`$. Notice however that the Hamiltonian vector field associated with the separable Hamiltonian $`H`$ is not (in general, and in the particular example of the Neumann system) bi-Hamiltonian.

3. The distribution $`D_H`$ is bi-Lagrangian, that is,

$$\{H_i,H_j\}=0,\{H_i,H_j\}^{}=0\text{for all }i,j\text{.}$$

As explained in , this condition characterizes the (integrable) Hamiltonian systems that are separable in the coordinates induced by $`L`$. We remark in passing that bi-Lagrangian foliations play an important role in the study of special Kähler manifolds .

4. As we have already said in the previous section, the separation coordinates $`(\lambda _i,\mu _i)`$ are canonical for $`\omega `$ (this is obvious) and for $`N`$,

$$N\frac{}{\lambda _i}=\lambda _i\frac{}{\lambda _i},N\frac{}{\mu _i}=\lambda _i\frac{}{\mu _i}.$$

They are often called Darboux-Nijenhuis coordinates (since $`N`$ is sometimes called the Nijenhuis tensor) and have been used as separation coordinates in, e.g., .

## 4 The case of the Neumann system

In this section we will exploit the separability condition (1) to find the variables of separation for the Neumann system. To this aim, we will look for a torsionless tensor field $`L`$ of type $`(1,1)`$ on $`S^2`$ such that its complete lift $`N`$ on $`T^{}S^2`$ satisfies

$$\omega _H(X_H,NX_H)=0,$$
(3)
where $`\omega _H=d(N^{}dH)`$ and $`H`$ is the Neumann Hamiltonian. Due to the form of the constraint and of the potential, it is quite natural to perform the computations in the local coordinates $`X=x^2`$, $`Y=y^2`$, parametrizing every connected component of the subset obtained by removing from $`S^2`$ its intersections with the coordinate planes. In these coodinates the Hamiltonian of the Neumann system has the form:

$$\begin{array}{cc}\hfill H=& \mathrm{\hspace{0.17em}2}\left[X(1X)p_X^22XYp_Xp_Y+Y(1Y)p_Y^2\right]\hfill \\ & +\frac{1}{2}(a_1a_3)X+\frac{1}{2}(a_2a_3)Y,\hfill \end{array}$$
(4)
where $`(p_X,p_Y)`$ are the conjugate momenta of $`(X,Y)`$. The unknown tensor field $`L`$ can be written as

$$\begin{array}{cc}& L^{}(dX)=AdX+BdY\hfill \\ & L^{}(dY)=CdX+DdY,\hfill \end{array}$$
(5)
where $`(A,B,C,D)`$ are functions of the coordinates $`(X,Y)`$ that must satisfy two additional conditions. The first one is that the torsion of $`L`$ has to vanish, that is,

$$\begin{array}{cc}& B\left(C_YD_X\right)BA_XDA_Y+AA_Y+CB_Y=0\hfill \\ & C\left(B_XA_Y\right)AD_XCD_Y+BC_X+DD_X=0.\hfill \end{array}$$

The second one is the independency of the eigenvalues of $`L`$. Our aim is to find $`(A,B,C,D)`$ in such a way to verify also the separability conditions (3). Thus we need to compute the complete lift $`N`$ of $`L`$, which turns out to be given by

$$\begin{array}{cc}\hfill N^{}dX& =AdX+BdY\hfill \\ \hfill N^{}dY& =CdX+DdY\hfill \\ \hfill N^{}dp_X& =\left[p_X\left(B_XA_Y\right)+p_Y\left(D_XC_Y\right)\right]dY+Adp_X+Cdp_Y\hfill \\ \hfill N^{}dp_Y& =\left[p_X\left(B_XA_Y\right)+p_Y\left(D_XC_Y\right)\right]dX+Bdp_X+Ddp_Y.\hfill \end{array}$$
(6)
Then (3) becomes a differential equation where the coordinates $`(X,Y)`$, the momenta $`(p_X,p_Y)`$, and the unknown functions $`(A,B,C,D)`$ and their derivatives appear polynomially. This suggests to seek for a solution which also depends polynomially on the coordinates. Trying the simplest solution, one finds

$$\begin{array}{cc}\hfill A& =(a_3a_1)X+a_1\hfill \\ \hfill B& =(a_3a_2)X\hfill \\ \hfill C& =(a_3a_1)Y\hfill \\ \hfill D& =(a_3a_2)Y+a_2,\hfill \end{array}$$
(7)
which leads to the tensor field $`L`$ defined by

$$\begin{array}{cc}\hfill L^{}(dX)& =a_1dX+Xd[(a_3a_1)X+(a_3a_2)Y]\hfill \\ \hfill L^{}(dY)& =a_2dY+Yd[(a_3a_1)X+(a_3a_2)Y].\hfill \end{array}$$
(8)
It can be checked that $`L`$ can be extended to the whole sphere. Indeed, it is the conformal Killing tensor associated with the spheroconical coordinates . They are the eigenvalues $`(\lambda _1,\lambda _2)`$ of $`L`$, since

$`det(\lambda IL)`$ $`=\lambda ^2\lambda [(a_3a_1)X+(a_3a_2)Y+(a_1+a_2)]`$
$`+a_2(a_3a_1)X+a_1(a_3a_2)Y+a_1a_2`$
$`=\lambda ^2(x^2+y^2+z^2)`$
$`\lambda [(a_3a_1)x^2+(a_3a_2)y^2+(a_1+a_2)(x^2+y^2+z^2)]`$
$`+a_2(a_3a_1)x^2+a_1(a_3a_2)y^2+a_1a_2(x^2+y^2+z^2)`$
$`=\lambda ^2(x^2+y^2+z^2)\lambda [(a_2+a_3)x^2+(a_1+a_3)y^2+(a_1+a_2)z^2]`$
$`+[a_2a_3x^2+a_1a_3y^2+a_1a_2z^2]`$
$`=(\lambda a_2)(\lambda a_3)x^2+(\lambda a_1)(\lambda a_3)y^2+(\lambda a_1)(\lambda a_2)z^2`$
$`=\left[{\displaystyle \frac{x^2}{\lambda a_1}}+{\displaystyle \frac{y^2}{\lambda a_2}}+{\displaystyle \frac{z^2}{\lambda a_3}}\right](\lambda a_1)(\lambda a_2)(\lambda a_3).`$
Hence,

$$\frac{(\lambda \lambda _1)(\lambda \lambda _2)}{(\lambda a_1)(\lambda a_2)(\lambda a_3)}=\frac{x^2}{\lambda a_1}+\frac{y^2}{\lambda a_2}+\frac{z^2}{\lambda a_3},$$

meaning that the eigenvalues of $`L`$ are the spheroconical cooordinates on $`S^2`$. We have thus deduced the usual separation coordinates of the Neumann system simply using the separability condition (1) discussed in Section 2.

## 5 Integrals of motion for the Neumann system

The aim of this section is to show that in our setting one can easily describe other interesting features of the Neumann system, such as a simple construction of the integrals of motion and the existence of a bi-Hamiltonian formulation in a suitable extension of the phase space. The crucial point is that the Hamiltonian $`H`$ of the Neumann system and the complete lift $`N`$ of the tensor field (8) satisfy the stronger condition

$$d(N^{}dH)=dp_1dH,\text{where }p_1=\text{tr}L=\lambda _1+\lambda _2\text{,}$$
(9)
implying the separability condition (1). This fact has important consequences, that we are going to show in the general setting considered in Sections 2 and 3.

Let $`L`$ be a torsionless tensor field of type $`(1,1)`$, with functionally independent eigenvalues, on an $`n`$-dimensional manifold $`𝒬`$, so that we can endow the symplectic manifold $`(T^{}𝒬,\{,\})`$ with the additional Poisson bracket

$$\{F,G\}^{}=\omega (NX_F,X_G),$$

using the complete lift $`N`$ of $`L`$. Let

$$det(\lambda IL)=\lambda ^n\left(p_1\lambda ^{n1}+p_2\lambda ^{n2}+\mathrm{}+p_n\right)$$

be the characteristic polynomial of $`L`$, and suppose that $`HC^{\mathrm{}}(T^{}𝒬)`$ satisfies

$$d(N^{}dH)=dp_1dH.$$
(10)
Then one can easily show that $`\omega _H=d(N^{}dH)`$ vanishes on the distribution $`D_H`$, so that the separability condition (1) is fulfilled and $`H`$ is separable in the coordinates associated with $`L`$. We know from Section 3 that there are (local) integrals of motion for $`X_H`$. If (10) holds, they can be easily found by writing it in the form $`d(N^{}dHp_1dH)=0`$, and choosing $`H_2`$ such that $`dH_2=N^{}dHp_1dH`$. Then (as shown in  in the case where $`H`$ is quadratic in the momenta) also the 1-form $`N^{}dH_2p_2dH`$ is closed, and the process can be iterated, so that the integrals of motion for $`X_H`$ can be determined from the recursion relations

$$\begin{array}{cc}& dH_2=N^{}dHp_1dH\hfill \\ & dH_3=N^{}dH_2p_2dH\hfill \\ & \mathrm{}\hfill \\ & dH_n=N^{}dH_{n1}p_{n1}dH.\hfill \end{array}$$
(11)
Moreover, we have that

$$0=N^{}dH_np_ndH.$$
(12)
For the Neumann system the global existence of the function $`H_2`$ is garanteed by the fact that $`T^{}S^2`$ is simply connected. One finds that, in local coordinates,

$`H_2`$ $`=2(a_1p_{Y}^{}{}_{}{}^{2}Y+a_2p_{X}^{}{}_{}{}^{2}X)(X+Y1)2a_3XY(p_Yp_X)^2`$
$`{\displaystyle \frac{1}{2}}a_2(a_1a_3)X{\displaystyle \frac{1}{2}}a_1(a_2a_3)Y,`$
which coincides with the constant of motion provided by the Lax matrix. Notice that the Hamiltonian vector field $`X_H`$ is not bi-Hamiltonian on $`T^{}𝒬`$. Nevertheless, from (12) we can conclude that

$$X_H=\{,H\}=\frac{1}{p_n}\{,H_n\}^{}.$$

In the terminology of , this means that $`X_H`$ is a quasi-bi-Hamiltonian vector field (see also ). Next we will show that we can obtain a bi-Hamiltonian representation of $`X_H`$ on the extended phase space $`:=T^{}𝒬\times `$, following . The first step is the extension of the bi-Hamiltonian structure from $`T^{}𝒬`$ to $``$. If $`c`$ is a coordinate in the “additional dimension” and $`FC^{\mathrm{}}(T^{}𝒬)`$, then the new Poisson brackets are defined as

$`\{F,c\}_{}`$ $`=`$ $`0`$
$`\{F,c\}_{}^{}`$ $`=`$ $`\{F,H\}c\{F,p_1\}.`$
They endow $``$ with a bi-Hamiltonian structure. The second step is to notice that the recursion relations (11) on the $`H_i`$ become the usual Lenard relations on the functions $`\widehat{H}_i`$ defined on $``$ as

$$\widehat{H}_0=c,\widehat{H}_i=Hcp_i\text{for }i=1,\mathrm{},n\text{.}$$

Indeed, we have

$$\begin{array}{ccc}\hfill \{\widehat{H}_0,\}_{}& =& 0\hfill \\ \hfill \{\widehat{H}_1,\}_{}& =& \{\widehat{H}_0,\}_{}^{}\hfill \\ & \mathrm{}& \\ \hfill \{\widehat{H}_n,\}_{}& =& \{\widehat{H}_{n1},\}_{}^{}\hfill \\ \hfill 0& =& \{\widehat{H}_n,\}_{}^{}.\hfill \end{array}$$

In the Neumann case, the restriction to $`c=0`$ of the bi-Hamiltonian vector field $`\{\widehat{H}_1,\}_{}=\{\widehat{H}_0,\}_{}^{}`$ is the Neumann vector field. Thus we have shown that the Neumann system admits a bi-Hamiltonian formulation in an extended phase space, recovering in this way a result of .

## 6 Concluding remarks

1. The extension process of a chain of the form (11) into a Lenard chain is the opposite of the reduction technique presented in , where a suitable quotienting produces a chain like (11) from a Lenard chain.

2. All the results presented here can be extended to the $`n`$-dimensional Neumann system. It would be interesting to compare our approach with the one based on the Lax representation.

3. The condition (9) appears also in , in the case where $`(𝒬,g)`$ is a Riemannian manifold and $`H=\frac{1}{2}g^{ij}(q)p_ip_j+V(q)`$. The authors show that (9) implies that $`L`$ is a conformal Killing tensor of $`g`$ with vanishing torsion. Notice however that in our approach to separability the Riemannian structure of $`𝒬`$ plays no role and the important objects are bi-Lagrangian foliations on bi-Hamiltonian manifolds.

4. The idea of extending the phase space in order to obtain a bi-Hamiltonian formulation of a given (Hamiltonian) vector field goes back, to the best of our knowledge, to , where it has been used also in the context of separation of variables.