text:
# 𝑑-wave pairing symmetry in cuprate superconductors

## 1 INTRODUCTION

Pairing symmetry in the cuprate superconductors was a subject of intense debate for many years after the discovery of high-temperature superconductivity , primarily because the interpretation of many conventional techniques, such as quasiparticle tunneling, NMR, angle-resolved photoemission, and penetration depth measurements, was model-dependent. Nevertheless, these phase-insensitive techniques have produced a large body of evidence for $`d`$-wave pairing in the cuprates. The recent development of phase-sensitive pairing symmetry test , has yielded compelling evidence for predominantly $`d`$-wave pairing symmetry in a number of optimally doped cuprates . A question naturally arises: How universal is the $`d_{x^2y^2}`$ pairing in cuprate superconductors?

There are numerous theoretical studies suggesting the stability of the $`d`$-wave pair state depends on the details of band structure and the pairing potential . There has also been considerable theoretical studies indicating  that a pure $`d_{x^2y^2}`$ pair state is not stable against the formation of time reversal symmetry breaking states such as $`d_{x^2y^2}+id_{xy}`$ or $`d_{x^2y^2}+is`$, at surfaces, interfaces, near impurities, or below a certain characteristic temperature. On the experimental side, Raman data on Bi2Sr2CaCu2O8+δ and Tl2Ba2CuO8+δ systems as a function of oxygen content indicate that the order parameter has $`d`$-wave symmetry near optimal doping and isotropic pairing in the over-doped regime . There are additional indirect experimental evidences for the fragility of the pure $`d`$-wave state . The following will examine the universality issues based mainly on the results of phase sensitive pairing symmetry experiments. We will conclude with a brief discussion on implications of $`d`$-wave superconductivity in the cuprates.

## 2 PHASE-SENSITIVE PAIRING SYMMETRY TESTS

The sign and magnitude of the pair tunneling critical current Ic across a Josephson junction made with at least one superconductor with unconventional pairing symmetry depends sensitively on the gap function symmetry and relative orientation of the junction electrodes. The ground state of a superconducting loop with an odd number of negative critical currents (termed a “frustrated” or “$`\pi `$” loop) is doubly degenerate and shows spontaneous magnetization of one-half flux quantum ($`\mathrm{\Phi }_0/2=hc/4e=1.035\times 10`$-7$`Gcm^2`$), provided that $`LI_c>>\mathrm{\Phi }_0`$, where $`L`$ is the self-inductance of the ring and $`I_c`$ is the critical current of the weakest junction in the loop . By varying the loop geometry, the presence or absence of the half-integer flux quantum effect can be used for a definitive determination of the pairing symmetry.

The tricrystal pairing symmetry tests use a ring consisting of three crystals with controlled orientation (see Fig. 1) to define the direction of the pair wavefunction.

The magnetic flux threading through the tricrystal ring is measured with a scanning SQUID microscope  for different tricrystal configurations to differentiate between various pairing symmetries. The $`\mathrm{\Phi }_0/2`$ effect is intrinsic in a frustrated geometry, and is observed in ring, disk, and blanket film samples at the tricrystal point.

## 3 UNIVERSALITY of the $`d`$-WAVE PAIR STATE

The establishment of $`d`$-wave pairing in some of the cuprates has imposed well defined constraints on possible models of high-temperature superconductivity, but it does not specify the high-Tc mechanism. To gain further insight into the origin of high-temperature superconductivity, it is important to determine whether the $`d_{x^2y^2}`$ pairing symmetry is universal. One can investigate the universality issues from the following aspects.

### 3.1 Various cuprate systems

Phase-sensitive SQUID interferometry , tricrystal magnetometry , SQUID magnetometry , and single junction interferometry  experiments have demonstrated $`d`$-wave pairing symmetry in YBCO single crystals or $`c`$-axis oriented epitaxial films. In addition, tricrystal magnetometry experiments have demonstrated $`d`$-wave pairing symmetry in various high-Tc hole-doped cuprate systems such as Tl2Ba2CuO6+δ, GdBa2Cu3O7-δ, and Bi2Sr2CaCu2O8+δ. (For a more complete tabulation of the phase-sensitive experiments see Ref. ). More recently, the pairing symmetry of the electron-doped cuprate superconductors Nd1.85Ce0.15CuO4-δ (NCCO) and Pr1.85Ce0.15CuO4-δ (PCCO) has been determined to also be predominantly $`d`$-wave by observing the half-flux quantum effect in $`c`$-axis thick blanket films ($``$1$`\mu `$m thick) epitaxially grown on tricrystal substrates with the configuration of Fig. 1 . Samples with two other geometries, designed to be unfrustrated for a $`d`$-wave superconductor, do not show the half-flux quantum effect. Shown in Fig. 2a is a scanning SQUID microscope image of a NCCO film

deposited on a frustrated tricrystal STO substrate, cooled to 4.2K in nominal zero field. This image was obtained by subtracting an image with an externally applied field of -0.15mG, which stabilized the half-quantum magnetic vortex with the flux pointing down, from an image taken with an applied field of 0.15mG, resulting in a half-quantum magnetic vortex pointing up, and dividing by 2. This results in an image of the half-quantum magnetic vortex with effects from surface roughness, a smoothly varying magnetic background, and the mutual inductance between the SQUID and sample subtracted out. Figure 2b shows fitting of cross-sections through the data of Fig. 2a to expressions for the magnetic signals expected from a Josephson vortex centered at the tricrystal point, with a total flux of $`\mathrm{\Phi }_0/2`$ . The height $`z`$ of the pickup loop was determined by fitting a bulk Abrikosov vortex. The solid lines in Fig. 2b are fits to cross-sectional data for the horizontal and diagonal grain boundaries, assuming a total flux of $`\mathrm{\Phi }_0/2`$ for this vortex, using the Josephson penetration depths $`\lambda _J`$ for the two grain boundaries as the sole fitting parameters. The dashed lines, assuming a flux of $`\mathrm{\Phi }_0`$, are much worse fits. These results indicate that the electron-doped superconductors, just as the hole-doped superconductors, have $`d`$-wave pairing symmetry.

### 3.2 Temperature dependence of $`\mathrm{\Phi }_0/2`$

It is of interest to find out whether $`d`$-wave pairing will change to $`s`$-wave, for example, as a function of temperature. However, the phase-sensitive experiments described above were all performed at liquid helium temperature. With a recently developed variable sample temperature scanning SQUID microscope , the temperature dependence of the half-integer flux quantum effect in tricrystal YBCO was studied. It was found that the total flux of the Josephson vortex at the tricrystal point remains constant (i.e. $`\mathrm{\Phi }_0/2`$) from 0.5K through Tc ($``$90K) . This finding means that the $`d`$-wave pair state dominates at all temperatures below Tc, and that there is little, if any, time-reversal symmetry breaking component in the order parameter, over the entire temperature range. This is consistent with earlier studies, which found less than a few percent of time-reversal symmetry breaking component in phase-sensitive measurements at low temperatures .

### 3.3 The effect of doping

Is there any band structure (doping) effect on pairing symmetry? In the competition between $`s`$-wave and $`d`$-wave channels for high-temperature superconductivity, is there any factor that can stabilize the $`d`$-wave pairing? In the literature, there have been some theoretical studies dealing with these issues: for example, based on a next-near-neighbor pairing interaction, the symmetry of a BCS condensate is predicted to vary as a function of energy band filling and other band parameters . The $`d`$-wave pairing channel is energetically favored in a wide range of doping centered around half-filling. The stability of the $`d`$-wave pair state may be enhanced by the proximity of the Fermi level to the van Hove singularity in the 2D bands of the CuO2 plane . A series of phase-sensitive pairing symmetry tests as a function of doping in a model system such as Tl2Ba2CuO6+δ, HgBa2CuO4+δ, Bi2Sr2CaCu2O8+δ, or Y1-xCaxBa2Cu3O7-δ could test for these effects.

## 4 CONCLUDING REMARKS

Phase-sensitive experiments have produced strong evidence for predominantly $`d`$-wave pairing symmetry in hole- and electron-doped cuprate superconductors. These tests also indicate, within the experimental accuracy ($`3\%`$ at low temperatures), the absence of pairing with broken time-reversal symmetry. The universality of $`d`$-wave pairing symmetry in bulk high-temperature superconductors is thus well-established. The possibility of pair states without time-reversal symmetry invariance (e.g. $`d_{x^2y^2}+id_{xy}`$ or $`d_{x^2y^2}+is`$) induced at surfaces, interfaces, or around impurities or defects is being actively investigated (see Ref.  and references therein). Recent observations of a splitting of the zero-bias conductance peak in quasiparticle tunneling  and spontaneous magnetization in $`c`$-axis YBCO thin films  represent possible evidence for such pair states.

Dominant $`d`$-wave pairing symmetry in both the hole- and electron-doped cuprates is expected from single-band models of the CuO2 planes, in which the electron-hole symmetry is an intrinsic property (see e.g. Poole et al. ). From another perspective, the predominance of $`d_{x^2y^2}`$ pairing in cuprates underscores the important role of strong electron correlation in determining the superconducting gap symmetry and other properties. It is experimentally and theoretically well established that a strong on-site Coulomb repulsion is present in all cuprates. Such a strong electron correlation causes the universally-seen Mott transition at half-filling. The strong on-site Coulomb repulsion rules out simple $`s`$-wave pairing.

In future, a systematic study of pairing symmetry as a function of doping, impurities, … is important for a better understanding of high-temperature superconductivity.

The authors wish to thank M. Bhushan, C.C. Chi, A. Gupta, Z.Z. Li, M.B. Ketchen, K.A. Moler, D.M. Newns, H. Raffy, Z.F. Ren, J.Z. Sun, G. Trafas, and J.H. Wang for invaluable assistance in the course of our tricrystal experiments.


format: text

text:
# Renormalization of spectral lineshape and dispersion below 𝑇_𝑐 in Bi2Sr2CaCu2O8+δ

## Abstract

Angle-resolved photoemission (ARPES) data in the superconducting state of Bi2Sr2CaCu2O8+δ show a kink in the dispersion along the zone diagonal, which is related via a Kramers-Krönig analysis to a drop in the low-energy scattering rate. As one moves towards $`(\pi ,0)`$, this kink evolves into a spectral dip. The occurrence of these anomalies in the dispersion and lineshape throughout the zone indicate the presence of a new energy scale in the superconducting state.

The high temperature superconductors exhibit many unusual properties, one of the most striking being the linear temperature dependence of the normal state resistivity. This behavior has been attributed to the presence of a quantum critical point, where the only relevant energy scale is the temperature. However, new energy scales become manifest below $`T_c`$ due to the appearance of the superconducting gap and resulting collective excitations. The effect of these new scales on the ARPES spectral function below $`T_c`$ have been well studied near the $`(\pi ,0)`$ point of the zone . In this Letter we show how these scales manifest themselves in the spectral functions over the entire Brillouin zone.

Remarkably, we find that these effects are manifest even on the zone diagonal where the gap vanishes, with significant changes in both the spectral lineshape and dispersion below $`T_c`$, relative to the normal state (where the nodal points exhibit quantum critical scaling ). Specifically, below $`T_c`$ a kink in the dispersion develops along the diagonal at a finite energy ($``$70 meV). This is accompanied, as required by Kramers-Krönig relations, by a reduction in the linewidth leading to well-defined quasiparticles. As one moves away from the node, the renormalization increases, and the kink in dispersion along the diagonal smoothly evolves into the spectral dip , with the same characteristic energy scale throughout the zone. We suggest that a natural interpretation of all of these spectral renormalizations is in terms of the electron interacting with a collective excitation below $`T_c`$, which is likely that seen directly by neutron scattering .

We begin our analysis by recalling  that, within the impulse approximation, the ARPES intensity for a quasi-two-dimensional system is given by  $`I(𝐤,\omega )=I_0(𝐤)f(\omega )A(𝐤,\omega )`$. Here $`𝐤`$ is the in-plane momentum, $`\omega `$ is the energy of the initial state relative to the chemical potential, $`f`$ is the Fermi function, $`I_0`$ is proportional to the dipole matrix element $`\left|M_{fi}\right|^2`$, and $`A`$ is the one-particle spectral function. Fig. 1 shows data  as a function of $`𝐤`$ and $`\omega `$.

$`A(𝐤,\omega )=(1/\pi )\mathrm{}mG(𝐤,\omega +i0^+)`$ can be written as

$$A(𝐤,\omega )=\frac{1}{\pi }\frac{\left|\mathrm{\Sigma }^{\prime \prime }(𝐤,\omega )\right|}{\left[\omega ϵ_𝐤\mathrm{\Sigma }^{}(𝐤,\omega )\right]^2+\left[\mathrm{\Sigma }^{\prime \prime }(𝐤,\omega )\right]^2}.$$
(1)
where the self-energy $`\mathrm{\Sigma }=\mathrm{\Sigma }^{}+i\mathrm{\Sigma }^{\prime \prime }`$ and $`ϵ_𝐤`$ is the bare dispersion. For $`𝐤`$ near $`k_F`$, and varying normal to the Fermi surface (shown in the inset in Fig. 1), we may write $`ϵ_𝐤v_F^0(kk_F)`$, where both $`k_F(\theta )`$ and the bare Fermi velocity $`v_F^0(\theta )`$ depend in general on the angle $`\theta `$ along the Fermi surface.

In Fig. 2a, we plot the dispersion of the spectral peak above $`T_c`$ obtained from constant $`𝐤`$ scans (energy distribution curves or EDCs), and the peak in momentum obtained from constant $`\omega `$ scans (momentum distribution curves or MDCs)  from data similar to Fig. 1. We find that the EDC and MDC peak dispersions are very different, a consequence of the $`\omega `$ dependence of $`\mathrm{\Sigma }`$. To see this, we note from Eq. (1) that the MDC at fixed $`\omega `$ is a Lorentzian centered at $`k=k_F+\left[\omega \mathrm{\Sigma }^{}(\omega )\right]/v_F^0`$, with a width (HWHM) $`W_M=|\mathrm{\Sigma }^{\prime \prime }(\omega )|/v_F^0`$, provided (i) $`\mathrm{\Sigma }`$ is essentially independent  of $`k`$ normal to the Fermi surface, and (ii) the dipole matrix elements do not vary significantly with $`k`$ over the range of interest. That these two conditions are fulfilled can be seen by the nearly Lorentzian MDC lineshape of the data in Fig. 1b.

On the other hand, in general, the EDC at fixed $`𝐤`$ (Fig. 1c) has a non-Lorentzian lineshape reflecting the non-trivial $`\omega `$-dependence of $`\mathrm{\Sigma }`$, in addition to the Fermi cutoff at low energies. Thus the EDC peak is not given by $`\left[\omega v_F^0(kk_F)\mathrm{\Sigma }^{}(\omega )\right]=0`$ but also involves $`\mathrm{\Sigma }^{\prime \prime }`$, unlike the MDC peak. Further, if the EDC peak is sharp enough, making a Taylor expansion we find that its width (HWHM) is given by $`W_E|\mathrm{\Sigma }^{\prime \prime }(E_k)|/[1\mathrm{\Sigma }^{}/\omega |_{E_k}]`$, where $`E_k`$ is the peak position.

We see that it is much simpler to interpret the MDC peak positions, and thus focus on the change in the MDC dispersion going from the normal (N) to the superconducting (SC) state shown in Fig. 2b. The striking feature of Fig. 2b is the development of a kink in the dispersion below $`T_c`$. At fixed $`\omega `$ let the dispersion change from $`k_N`$ to $`k_{SC}`$. Using $`v_F^0(k_Nk_{SC})=\mathrm{\Sigma }_{SC}^{}(\omega )\mathrm{\Sigma }_N^{}(\omega )`$, we directly obtain the change in real part of $`\mathrm{\Sigma }`$ plotted in Fig. 2c. The Kramers-Krönig transformation of $`\mathrm{\Sigma }_{SC}^{}\mathrm{\Sigma }_N^{}`$ then yields $`\mathrm{\Sigma }_N^{\prime \prime }\mathrm{\Sigma }_{SC}^{\prime \prime }`$, plotted in Fig. 2d, which shows that $`|\mathrm{\Sigma }_{SC}^{\prime \prime }|`$ is smaller than $`|\mathrm{\Sigma }_N^{\prime \prime }|`$ at low energies.

We compare these results in Fig. 3a with the $`W_M=|\mathrm{\Sigma }^{\prime \prime }|/v_F^0`$ estimated directly from the MDC Lorentzian linewidths. The normal state curve was obtained from a linear fit to the corresponding MDC width data points in Fig. 3a, and then the data from Fig. 2d was added to it to generate the low temperature curve. We are thus able to make a quantitative connection between the appearance of a kink in the (MDC) dispersion below $`T_c`$ and a drop in the low energy scattering rate in the SC state relative to the normal state, which leads to the appearance of quasiparticles below $`T_c`$ . We emphasize that we have estimated these $`T`$-dependent changes in the complex self-energy without making fits to the EDC lineshape, thus avoiding the problem of modeling the $`\omega `$ dependence of $`\mathrm{\Sigma }`$ and the extrinsic background.

In Fig. 3b, we plot the EDC width obtained as explained in  from Fig. 3d. As an interesting exercise, we present in Fig. 3c the ratio of this EDC width to the MDC width of Fig. 3a (dotted lines), and compare it to the renormalized MDC velocities, $`1/vdk/d\omega `$, obtained directly by numerical differentiation of Fig. 2b (solid lines). We note that only for a sufficiently narrow EDC lineshape is the ratio $`W_E/W_Mv_F^0/[1\mathrm{\Sigma }^{}/\omega ]=v_F`$. Interestingly, only in the SC state below the kink energy do these two quantities agree, which implies that only in this case does one have a Fermi liquid.

Similar kinks in the dispersion have been seen by ARPES in normal metals due to the electron-phonon interaction. Phonons cannot be the cause here, since our kink disappears above $`T_c`$. Rather, our effect is suggestive of coupling to an electronic collective excitation which only appears below $`T_c`$.

We now study how the lineshape and dispersion evolve as we move along the Fermi surface. Away from the node a quantitative analysis (like the one above) becomes more complicated  and will be presented in a later publication. Here, we will simply present the data. In Fig. 4, we plot raw (2D) data as obtained from our detector for a series of cuts parallel to the $`MY`$ direction (normal state in left panels, superconducting state in middle panels). We start from the bottom row that corresponds to a cut close to the node and reveals the same kink described above. As we move towards $`(\pi ,0)`$, the dispersion kink (middle panels) becomes more pronounced and at around $`k_x`$=0.55 develops into a break separating the faster dispersing high energy part of the spectrum from the slower dispersing low energy part. This break leads to the appearance of two features in the EDCs, shown in the right panels of Fig. 4. Further towards $`(\pi ,0)`$, the low energy feature, the quasiparticle peak, becomes almost dispersionless. At the $`(\pi ,0)`$ point, this break effect becomes the most pronounced, giving rise to the well known peak/dip/hump  in the EDC. We note that there is a continuous evolution in the zone from kink to break, and these features all occur at exactly the same energy.

The above evolution is suggestive of the self-energy becoming stronger as the $`(\pi ,0)`$ point is approached. This can be quantified from the observed change in the dispersion. In Fig. 3(e) we plot the ratio of the EDC dispersion slope above and below the kink energy at various points along the Fermi surface obtained from middle panels of Fig. 4. Near the node, this ratio is around 2, but becomes large near the $`(\pi ,0)`$ point because of the nearly dispersionless quasiparticle peak.

The lineshape also indicates that the self-energy is larger near $`(\pi ,0)`$, as is evident in Fig. 5. Along the diagonal, there is a gentle reduction in $`\mathrm{\Sigma }^{\prime \prime }`$ at low energies, as shown in Fig. 3 (a) and (b), with an onset at the dispersion kink energy scale. In contrast, near the $`(\pi ,0)`$ point there must be a very rapid change in $`\mathrm{\Sigma }^{\prime \prime }`$ in order to produce a spectral dip, as quantified in Refs. . Despite these differences, it is important to note that these changes take place throughout the zone at the same characteristic energy scale (vertical line in Fig. 5).

As discussed in Ref.  the near-$`(\pi ,0)`$ ARPES spectra can be naturally explained in terms of the interaction of the electron with a collective mode of electronic origin which only exists below $`T_c`$. It was further speculated that this mode was the neutron resonance , an interpretation which received further support from Ref.  where the doping dependence of ARPES spectra were examined. Here we have shown that dispersion and lineshape anomalies have a continuous evolution throughout the zone and are characterized by a single energy scale. This leads us to suggest that the same electron-mode interaction determines the superconducting lineshape and dispersion at all points in the zone, including the nodal direction . In essence, there is a suppression of the low energy scattering rate below the finite energy of the mode. Of course, since the neutron mode is characterized by a $`(\pi ,\pi )`$ wavevector, one would expect its effect on the lineshape to be much stronger at points in the zone which are spanned by $`(\pi ,\pi )`$ , as observed here.

In summary, we have shown by a simple, self-consistent analysis based on general properties of the spectral function and self-energy, that Bi2Sr2CaCu2O8+δ shows a dispersion renormalization along the zone diagonal which is directly related to a drop in the low energy scattering rate below $`T_c`$. The anomalies in the dispersion and lineshape evolve smoothly as one moves from the zone diagonal to the zone corner, but always show the same characteristic energy scale. We suggest that this suppression of the scattering rate below $`T_c`$ at all points in the Brillouin zone is due to the presence of a gap and a finite energy collective mode, which we identify with the magnetic resonance observed by neutron scattering.

MR would like to thank P.D. Johnson for discussions. This work was supported by the NSF DMR 9974401, the U.S. DOE, Basic Energy Sciences, under contract W-31-109-ENG-38, the CREST of JST, and the Ministry of Education, Science, and Culture of Japan. The Synchrotron Radiation Center is supported by NSF DMR 9212658. JM is supported by the Swiss National Science Foundation, and MR in part by the Indian DST through the Swarnajayanti scheme.

Note added: After completion of this work, we became aware of related work by Bogdanov et al., cond-mat/0004349.


text:
# Toda equations for Hurwitz numbers

## 1 Introduction

### 1.1

In a recent paper , R. Pandharipande showed that the conjectural Toda equation for the Gromov-Witten potential $`F_^1`$ of the Riemann sphere $`^1`$ implies a Toda equation for the generating function of the Hurwitz numbers $`H_{g,d}`$. Here $`H_{g,d}`$ is the number of degree $`d`$ connected coverings of $`^1`$ by a curve of genus $`g`$ which have simple ramifications over some $`2g+2d2`$ fixed points of $`^1`$. This Toda equation results in a simple recurrence relation for the numbers $`H_{g,d}`$.

The purpose of this paper is to prove that a more general generating function, where arbitrary ramifications are allowed over two given points of $`^1`$ (for example, over $`0,\mathrm{}^1`$) satisfies an entire hierarchy of PDE’s, namely the Toda lattice hierarchy of Ueno and Takasaki . This proves and generalizes the result of Pandharipande.

We will call these more general numbers of ramified coverings the double Hurwitz numbers.

### 1.2

Our approach is quite straightforward: we remark that the generating function for double Hurwitz numbers is almost by definition a certain matrix element in the infinite wedge space. It is a well known result of the Kyoto school that such matrix elements are $`\tau `$-functions of integrable hierarchies.

These Toda equations should probably be not too surprising from the physical point of view. Physicists consider the enumeration of coverings of $`^1`$ as a sort of string theory on $`^1`$, see for example , and relate it to a number of other models where the occurrence of integrable hierarchies is common.

### 1.3

It is known, see  and also , that ordinary Hurwitz numbers, that is, the numbers of coverings where the ramification over $`\mathrm{}^1`$ is arbitrary and all other ramifications are simple, can be expressed as certain Hodge integrals over moduli spaces of curves. It is an interesting problem to transfer the Toda equation to the Hodge integrals side. In particular, it would be nice to reproduce the KdV equations or Virasoro constraints, see for example , for just the Gromov-Witten potential of a point, which corresponds to Hodge integrals with no $`\lambda `$-classes.

It is also interesting and challenging to find an expression for the double Hurwitz numbers in terms of the GW theory of $`^1`$. The conjectural Toda equation for the GW potential of $`^1`$, see , may be an indication of such a connection.

### 1.4

The Toda lattice equations is an infinite number of PDE’s satisfied by the generating function for coverings, each resulting is some recurrence relations for the double Hurwitz numbers. These relations uniquely determine the double Hurwitz numbers.

In particular, the lowest equation of the Toda hierarchy, see Sections 2.8 and 2.9, describes what will happen if we add an extra sheet to our covering leaving the nontrivial part of the monodromy the same.

Many recurrent relations for Hurwitz numbers are known in the literature, see for example  and references therein, but as R. Pandharipande writes in , the form of the Toda recurrence relations is “certainly the simplest recursion for the Hurwitz numbers known to the author” . It is absolutely beyond the scope of the present note to give any sort of survey on the vast literature on Hurwitz numbers and their applications.

### 1.5

I want to thank Rahul Pandharipande for helpful correspondence.

## 2 Coverings and Toda equations

### 2.1

Let $`S(d)`$ denote the symmetric group and fix some conjugacy classes

$$C_1,\mathrm{},C_sS(d).$$

Let $`\mathrm{Cov}_d(C_1,\mathrm{},C_s)`$ denote the weighted number of $`d`$-fold coverings of $`^1`$ ramified over $`s`$ fixed points of $`^1`$ with monodromies in the conjugacy classes $`C_1,\mathrm{},C_s`$. The weight of any covering in $`\mathrm{Cov}_d`$ is the reciprocal of the order of its group of automorphisms. This automorphisms group is the same as the centralizer of all monodromies in $`S(d)`$. Note that $`\mathrm{Cov}_d`$ counts all coverings, that is, including the disconnected ones.

For this weighted number there is the following formula, see for example , which goes back to Burnside

$$\mathrm{Cov}_d(C_1,\mathrm{},C_s)=\underset{|\lambda |=d}{}\left(\frac{dim\lambda }{d!}\right)^2\underset{i=1}{\overset{s}{}}f_{C_i}(\lambda ),$$
(1)
where the sum is over all partitions $`\lambda `$ of $`d`$, $`dim\lambda `$ is the dimension of the corresponding representation of the symmetric group $`S(d)`$, and

$$f_{C_i}(\lambda )=|C_i|\frac{\chi ^\lambda (C_i)}{dim\lambda }.$$
(2)
Here $`|C_i|`$ is the cardinality of the conjugacy class $`C_i`$ and $`\chi ^\lambda (C_i)`$ is the value of the irreducible character $`\chi ^\lambda `$ on any element of $`C_i`$. The expression (2) is in fact a polynomial in $`\lambda `$ in a certain natural sense, see for example .

### 2.2

Recall that the Schur functions $`s_\lambda `$ are certain distinguished symmetric functions associated to any partition $`\lambda `$. In terms of the power-sum symmetric functions $`p_1,p_2,\mathrm{}`$ they have the following expression, see ,

$$s_\lambda (p_1,p_2,\mathrm{})=\frac{1}{d!}\underset{|\mu |=d}{}\chi ^\lambda (C_\mu )|C_\mu |p_\mu ,$$

where the summation is over all partitions $`\mu `$ of $`d=|\lambda |`$, $`C_\mu `$ is the conjugacy class corresponding to the partition $`\mu `$, and $`p_\mu =p_{\mu _i}`$. We will treat the $`p_k`$’s as independent variables.

Now let $`P=(p_1,p_2,\mathrm{})`$ and $`P^{}=(p_1^{},p_2^{},\mathrm{})`$ be two sets of variables. We have the following equality:

$$\begin{array}{c}\underset{d=0}{\overset{\mathrm{}}{}}q^d\underset{|\mu |=|\nu |=d}{}p_\mu p_\nu ^{}\mathrm{Cov}_d(C_\mu ,C_\nu ,C_1,\mathrm{},C_s)=\hfill \\ \hfill \underset{\lambda }{}q^{|\lambda |}s_\lambda (P)s_\lambda (P^{})\underset{i=1}{\overset{s}{}}f_{C_i}(\lambda ),\end{array}$$
(3)
where $`q`$ is a new variable which keeps track of the degree $`d`$ of the covering. The variables $`q`$ is redundant since for every term its degree $`d`$ in $`q`$ is equal to both its degree in $`P`$ and in $`P^{}`$, where

$$\mathrm{deg}p_k=\mathrm{deg}p_k^{}=k,$$

but it is convenient to keep it.

### 2.3

Let $`C_{(2)}`$ be the conjugacy class of a transposition. Introduce the following generating function

$$\tau (P,P^{},\beta ,q)=\underset{d,b,\mu ,\nu }{}q^d\beta ^bp_\mu p_\nu ^{}\mathrm{Cov}_d(C_\mu ,C_\nu ,\underset{b\text{ times}}{\underset{}{C_{(2)},\mathrm{},C_{(2)}}})/b!.$$

This is the generating function for degree $`d`$ possibly disconnected coverings of $`^1`$ such that the ramification over $`0,\mathrm{}^1`$ can be fixed arbitrarily and the other $`b`$ ramifications are simple. It is a general property of the exponential generating functions, see e.g. Chapter 3 in , that

$$\mathrm{log}\tau (P,P^{},\beta ,q)=H(P,P^{},\beta ,q),$$

where $`H(P,P^{},\beta ,q)`$ is the generating function for connected coverings

$$H(P,P^{},\beta ,q)=\underset{d,b,\mu ,\nu }{}q^d\beta ^bp_\mu p_\nu ^{}\mathrm{Hur}_{d,b}(\mu ,\nu )/b!,$$

Here $`\mathrm{Hur}_{d,b}(\mu ,\nu )`$ is the weighted number of *connected* degree $`d`$ coverings of $`^1`$ with monodromy around $`0,\mathrm{}^1`$ being $`\mu `$ and $`\nu `$, respectively, and $`b`$ additional simple ramifications. The genus of such a covering is

$$g=(b+2\mathrm{}(\mu )\mathrm{}(\nu ))/2,$$

where $`\mathrm{}(\mu )`$ is the number of parts of $`\mu `$. The weight of each covering in $`\mathrm{Hur}_{d,b}(\mu ,\nu )`$ is the reciprocal of the order of its automorphism group. We call the numbers $`\mathrm{Hur}_{d,b}(\mu ,\nu )`$ the *double Hurwitz numbers*.

We also point out that if when one of ramifications $`\mu `$ and $`\nu `$ is trivial and the covering is connected of degree $`d>2`$ then the automorphism group is trivial and so $`H=\mathrm{log}\tau `$ counts any such covering with weight $`1`$.

### 2.4

Denote by $`f_2`$ the polynomial (2) corresponding to the class $`C_{(2)}`$ of a transposition. From (3) we have

$$\tau (P,P^{},\beta ,q)=\underset{\lambda }{}q^{|\lambda |}e^{\beta f_2(\lambda )}s_\lambda (P)s_\lambda (P^{}).$$
(4)
Explicitly, the polynomial $`f_2`$ is given by

$$f_2=\frac{1}{2}\underset{i}{}\left[\left(\lambda _ii+\frac{1}{2}\right)^2\left(i+\frac{1}{2}\right)^2\right].$$
(5)
Remark that the upper limit of summation can be taken to be $`+\mathrm{}`$ in the above sum since $`\lambda _i=0`$ for all sufficiently large $`i`$. Using some standard combinatorics, see for example Section 5 of , one can rewrite the formula (5) as follows

$$f_2=\underset{k𝔖(\lambda )_+}{}\frac{k^2}{2}\underset{k𝔖(\lambda )_{}}{}\frac{k^2}{2}$$
(6)
where

$$𝔖(\lambda )=\{\lambda _ii+\frac{1}{2}\}+\frac{1}{2}$$

and, by definition, for any $`S+\frac{1}{2}`$ we set

$$S_+=S\left(_0\frac{1}{2}\right),S_{}=\left(_0\frac{1}{2}\right)S.$$

### 2.5

The sum (4) admits a straightforward interpretation as a certain matrix element in the infinite wedge space. It will be convenient for us to use the notation of , see especially the Appendix to  for a summary of the infinite wedge basics. A good general reference on the infinite wedge space is Chapter 14 of .

Note that the power-sum variables $`P`$ and $`P^{}`$ are related to the variables $`t`$ and $`t^{}`$ used in  by

$$t_k=\frac{p_k}{k},t_k^{}=\frac{p_k^{}}{k},k=1,2,\mathrm{},$$

We will use the following abbreviations

$$\mathrm{\Gamma }_+=\mathrm{\Gamma }_+(p_1,\frac{p_2}{2},\frac{p_3}{3},\mathrm{}),\mathrm{\Gamma }_{}=\mathrm{\Gamma }_{}(p_1^{},\frac{p_2^{}}{2},\frac{p_3^{}}{3},\mathrm{}).$$

The formula (4) becomes

$$\tau (P,P^{},\beta ,q)=(\mathrm{\Gamma }_+q^He^{\beta F_2}\mathrm{\Gamma }_{}v_{\mathrm{}},v_{\mathrm{}}),$$
(7)
where $`v_{\mathrm{}}`$ is the vacuum vector of the infinite wedge space, $`H`$ is the energy operator, and $`F_2`$ is the following operator

$$F_2=\underset{k+\frac{1}{2}}{}\frac{k^2}{2}:\psi _k\psi _k^{}:,$$

where the colons denote the normal ordering. It is clear from (6) that the operator $`F_2`$ acts as follows

$$F_2v_\lambda =f_2(\lambda )v_\lambda ,$$

where

$$v_\lambda =\underset{¯}{s_1}\underset{¯}{s_2}\underset{¯}{s_3}\mathrm{},\{s_1,s_2,\mathrm{}\}=𝔖(\lambda ),$$

is the vector corresponding to the partition $`\lambda `$.

### 2.6

Any operator of the form $`\xi =\psi _i\psi _j^{}`$ satisfies the relation

$$[\xi 1+1\xi ,\mathrm{\Omega }]=0,\mathrm{\Omega }=\psi _k\psi _k^{},$$

Since the operators $`q^H`$ and $`e^{\beta F_2}`$ are exponentials of sums of such operators $`\xi `$, we have

$$[q^He^{\beta F_2}q^He^{\beta F_2},\mathrm{\Omega }]=0.$$

It follows that the sequence

$$\tau _n=(\mathrm{\Gamma }_+q^He^{\beta F_2}\mathrm{\Gamma }_{}v_n,v_n),n,$$
(8)
where

$$v_n=\underset{¯}{n\frac{1}{2}}\underset{¯}{n\frac{3}{2}}\underset{¯}{n\frac{5}{2}}\mathrm{},$$

is the vacuum vector in the charge $`n`$ subspace, is a sequence of $`\tau `$-functions for the Toda lattice hierarchy. Our old $`\tau (P,P^{},\beta ,q)`$ is the $`\tau _0`$ term of this sequence.

### 2.7

We will now show that all terms of the sequence (8) can be expressed in terms of $`\tau (P,P^{},\beta ,q)`$. We have $`v_n=R^nv_{\mathrm{}}`$, where $`R`$ is the translation operator, and also $`[R,\mathrm{\Gamma }_\pm ]=0`$. We compute

$`R^nF_2R^n`$ $`={\displaystyle \underset{k>0}{}}{\displaystyle \frac{k^2}{2}}\psi _{kn}\psi _{kn}^{}{\displaystyle \underset{k<0}{}}{\displaystyle \frac{k^2}{2}}\psi _{kn}^{}\psi _{kn}`$
$`={\displaystyle \underset{k+\frac{1}{2}}{}}{\displaystyle \frac{(k+n)^2}{2}}:\psi _k\psi _k^{}:+{\displaystyle \underset{k=1/2}{\overset{n1/2}{}}}{\displaystyle \frac{k^2}{2}}`$
$`=F_2+nH+{\displaystyle \frac{n^2}{2}}C+{\displaystyle \frac{n(4n^21)}{24}},`$
where $`H`$ and $`C`$ are the energy and charge operators. Similarly,

$$R^nHR^n=H+nC+\frac{n^2}{2}.$$

Since the charge operator $`C`$ commutes with $`\mathrm{\Gamma }_\pm `$ and $`Cv_{\mathrm{}}=0`$, we have

$`\tau _n`$ $`=(\mathrm{\Gamma }_+q^He^{\beta F_2}\mathrm{\Gamma }_{}v_n,v_n)`$
$`=(\mathrm{\Gamma }_+R^nq^He^{\beta F_2}R^n\mathrm{\Gamma }_{}v_{\mathrm{}},v_{\mathrm{}})`$
$`=q^{n^2/2}e^{n(4n^21)\beta /24}(\mathrm{\Gamma }_+(e^{n\beta }q)^He^{\beta F_2}\mathrm{\Gamma }_{}v_{\mathrm{}},v_{\mathrm{}})`$
$`=q^{n^2/2}e^{n(4n^21)\beta /24}\tau (P,P^{},\beta ,e^{n\beta }q).`$
Substituting this expression into the Hirota bilinear equations for the Toda lattice hierarchy, we obtain the following

###### Theorem.

The function $`\tau (P,P^{},\beta ,q)`$ satisfies the following Hirota-type bilinear equations:

$$q^{m+1}e^{m(m+1)\beta /2}[z^{1m}]e^{2{\scriptscriptstyle ns_n/z^n}}\mathrm{}\mathrm{}=[z^{m+1}]e^{2{\scriptscriptstyle ns_n^{}/z^n}}\mathrm{}\mathrm{},$$
(9)
where $`m`$, $`z`$ is a formal variable, $`[z^{m+1}]`$ denotes the coefficient of $`z^{m+1}`$, $`S=(s_1,s_2,\mathrm{})`$ and $`S^{}=(s_1^{},\mathrm{})`$ are two arbitrary sequences,

$`\mathrm{}`$ $`=\tau (P+S,P^{}+S^{}+\stackrel{}{z},\beta ,e^{(m+1)\beta }q),`$ $`\mathrm{}`$ $`=\tau (P+S\stackrel{}{z},P^{}+S^{},\beta ,e^{m\beta }q),`$
$`\mathrm{}`$ $`=\tau (PS,P^{}S^{}\stackrel{}{z},\beta ,e^\beta q),`$ $`\mathrm{}`$ $`=\tau (PS+\stackrel{}{z},P^{}S^{},\beta ,q),`$
and $`\stackrel{}{z}=(z,z^2,z^3,\mathrm{})`$ .

Remark that, as pointed out earlier,

$`\tau (P,P^{},\beta ,q)`$ $`=\tau (qp_1,q^2p_2,q^3p_3,\mathrm{},P^{},\beta ,1)`$
$`=\tau (P,qp_1^{},q^2p_2^{},q^3p_3^{},\mathrm{},\beta ,1).`$
### 2.8

In particular, setting $`m=0`$ and taking the coefficient of $`s_1`$ gives

$$\tau \frac{^2\tau }{p_1p_1^{}}\frac{\tau }{p_1}\frac{\tau }{p_1^{}}=q\tau (\mathrm{},e^\beta q)\tau (\mathrm{},e^\beta q).$$

This is equivalent to

$$\frac{^2}{p_1p_1^{}}\mathrm{log}\tau =q\frac{\tau (\mathrm{},e^\beta q)\tau (\mathrm{},e^\beta q)}{\tau ^2}.$$
(10)
Recall that

$$\mathrm{log}\tau =H$$

is the generating function for the connected coverings.

The combinatorial meaning of the operator $`{\displaystyle \frac{^2}{p_1p_1^{}}}`$ is that it peels off a sheet of the covering which is unramified over the special points $`0,\mathrm{}^1`$. The equation (10) gives thus a recursive procedure for the computation of the double Hurwitz numbers with a fixed nontrivial part of the monodromy around $`0,\mathrm{}`$.

### 2.9

In the particular case when there is no ramification over $`0,\mathrm{}`$, that is, when,

$$p_2=p_3=\mathrm{}=p_2^{}=p_3^{}=\mathrm{}=0$$

the function $`\tau `$ depends on the variables $`p_1`$, $`p_1^{}`$, and $`q`$ only via their product $`qp_1p_1^{}`$. If we set

$$u=\mathrm{log}qp_1p_1^{}$$

then the equation (10) becomes

$$e^u\frac{^2}{u^2}\mathrm{log}\tau =\frac{\tau (u+\beta ,\beta )\tau (u\beta ,\beta )}{\tau ^2},$$

which is equivalent to the Toda equation derived in .


text:
# Temporal variations of the rotation rate in the solar interior

## 1 Introduction

Helioseismic data allow us to probe the rotation rate in the solar interior as a function of radius and latitude (Thompson et al. 1996; Schou et al. 1998). With the accumulation of Global Oscillations Network Group (GONG) and Michelson Doppler Imager (MDI) data over the last five years, it has also become possible to study the temporal variation in the rotation rate and other properties of the solar interior.

It is well known that frequencies of solar oscillations vary with time (Elsworth et al. 1990; Libbrecht & Woodard 1990; Dziembowski et al. 1998; Bhatnagar, Jain & Tripathy 1999; Howe, Komm & Hill 1999). These variations are also known to be correlated with solar activity. However, these frequency changes appear to result from variations in solar structure close to the surface; there is no clear evidence for structural changes in deep interior (Basu & Antia 2000a). Inversions for the rotation rate in solar interior on the other hand, have confirmed that the temporal variations in the rotation rate penetrate to somewhat deeper layers (Howe et al. 2000a; Toomre et al. 2000). Furthermore, the changing pattern of solar rotation rate in interior also agrees with the torsional oscillations observed at the solar surface (Howard & LaBonte 1980). While most of these studies are restricted to the outer part of the convection zone, the seat of the solar dynamo is generally believed to be near the base of the convection zone and rotation is believed to play an important role in the operation of the solar dynamo. This is also the region where the tachocline is located (Kosovichev 1996; Basu 1997). Thus one needs to look for possible changes in rotation rate in this region. Recently, Howe et al. (2000b) have reported a 1.3 year periodicity in variation of equatorial rotation rate at a radial distance of $`0.72R_{}`$. It is not clear if the period of 1.3 years is associated with solar cycle variations. Considering the fact that we have only 4 years of data, this period is also uncomfortably close to the orbital period of the Earth and one may expect systematic errors on this period being introduced in the data. If this periodicity is confirmed, it may have implications for the theories of solar dynamo. Hence, this needs to be checked by an independent analysis.

In this work we use oscillation frequencies obtained from GONG observations covering a period from May 1995 to June 1999, to study possible temporal variations in the rotation rate as a function of radius and latitude. We use both 1.5d and 2d inversion using the Regularized Least Squares (RLS) technique (Antia, Basu & Chitre 1998) to determine the rotation rate from the observed splitting coefficients.

## 2 The data and technique

We use data for GONG months 1–44 to determine the rotation rate in the solar interior. Each of these data sets covers a period of 108 days and is centered on dates 36 days apart. Thus there is significant overlap between neighboring data sets. Each GONG month covers a period of 36 days with month 1, starting on May 7, 1995 and month 44 ending on September 6, 1999. We have used 42 data sets centered on GONG months 2–43 to study the temporal variations in the rotation rate. The GONG project determines the frequencies of each mode $`\nu _{n\mathrm{}m}`$ (Hill et al. 1996) where $`n,\mathrm{},m`$ are respectively, the radial order, degree and azimuthal order of the mode. The frequencies of all observed modes in an $`n,\mathrm{}`$ multiplet are fitted to an expansion in terms of orthogonal polynomials (Ritzwoller & Lavely 1991) in $`m`$,

$$\nu _{n\mathrm{}m}=\nu _n\mathrm{}+\underset{j=1}{\overset{j_{\mathrm{max}}}{}}c_j(n,\mathrm{})𝒫_j^{(\mathrm{})}(m),$$
(1)
where $`c_j`$ are the splitting coefficients which are determined by least squares fit to the frequencies $`\nu _{n\mathrm{}m}`$, and the number of terms $`j_{\mathrm{max}}`$ is set such that $`j_{\mathrm{max}}2\mathrm{}`$. For this work, we have fitted splitting coefficients up to $`c_{16}`$; beyond this there is very little signal in the fitted coefficients, and because of the orthogonality of these polynomials, inclusion of higher coefficients does not make any difference to the fitted values of the lower order coefficients. The odd-order coefficients $`c_1,c_3,\mathrm{}`$ are determined by the rotation rate in solar interior and hence these can be used to determine rotation rate as a function of latitude and radius. The splitting coefficients are sensitive only to the north-south symmetric component of rotation rate and hence that is the only component that can be determined through inversions. We therefore restrict ourselves to this symmetric component. The anti-symmetric component in the outer layers can be studied using other techniques such as ring diagram analysis (e.g., Basu, Antia & Tripathy 1999; Basu & Antia 2000b). We use splitting coefficients for modes with $`1.5<\nu <3.5`$ mHz and $`\mathrm{}150`$ only to perform the inversions of the different GONG data sets.

Traditionally (e.g., Snodgrass 1984), the latitudinal dependence of the solar surface rotation rate is expressed as a second degree polynomial in $`\mathrm{cos}^2\theta `$, $`\theta `$ being the colatitude. This is the “smooth” part of the rotation rate. The torsional oscillations at the solar surface have been studied using the residuals left after subtracting the mean rotation rate obtained from this second degree polynomial fit (Howard & LaBonte 1980). Since the signal in these residuals is rather weak other prescriptions for calculating the “mean” rotation rate have also been tried to study the sensitivity of resulting pattern to these changes (LaBonte & Howard 1982; Snodgrass 1992). The polynomial approximation to the rotation rate is essentially equivalent to using splitting coefficients $`c_1,c_3`$, and $`c_5`$ only for inversions. Thus we can determine the smooth part of rotation rate by using only these three coefficients and take a mean over the period for which data are available. Since helioseismic data are available only for a small fraction of the magnetic cycle, the average will not represent a true mean of rotation rate over long time intervals, but that is the best we can do with the available data sets. This mean can be subtracted from the rotation rate at any given epoch to get the time varying component of the rotation rate. Thus we can write

$$\delta \mathrm{\Omega }(r,\theta ,t)=\mathrm{\Omega }(r,\theta ,t)\mathrm{\Omega }_s(r,\theta ,t),$$
(2)
where $`\mathrm{\Omega }`$ is the rotation rate at any given epoch and $`\mathrm{\Omega }_s`$ is the smooth part of rotation rate as a function of latitude. Here the angular brackets denote average over the time duration for which data are available. This time varying component $`\delta \mathrm{\Omega }`$ is generally called the zonal flow. We will refer to this time varying component variously as the residual in rotation rate or the zonal flow. It may be noted that this residual cannot be directly compared with the observations at solar surface as different mean rotation would have been subtracted in the two cases.

We use both 1.5d and 2d RLS techniques as described by Antia et al. (1998) to invert the splitting coefficients to determine the rotation rate as a function of radius and latitude for each data set. The advantage of the 1.5d inversion is that each coefficient is treated separately and hence we can combine the results from required coefficients as per the need for smooth part and full rotation rate. In the 2d inversion technique one has to do separate inversions to calculate these parts independently. As a result, for 2d inversion it is more convenient to take the temporal mean of the entire rotation rate and subtract it from the rotation rate at each epoch to get the residual velocity. This prescription has also been used in literature (e.g., LaBonte & Howard 1982; Howe et al. 2000a) to isolate temporal variation in rotation rate. Unless otherwise stated, all results in this paper have been obtained by subtracting the mean of only the smooth component of rotation rate. The use of two independent inversion techniques enable us to test the results for possible systematic errors in inversion.

Observations from the Earth or the first Lagrange point, where the SOHO satellite is located, will only give the synodic rotation rate, which needs to be converted to sidereal value. This conversion may introduce some error, which will have periodic variations on the time-scale of one year. There may be other systematic errors in inverted rotation rate, but most of these may be independent of time and may get cancelled when the mean rotation rate is subtracted.

## 3 Results

The zonal flows obtained by using the 1.5d inversion technique are shown in Fig. 1 as a function of latitude for a few data sets and depths. These results were obtained by subtracting the average of the smooth rotation rate (as obtained from coefficients $`c_1`$, $`c_3`$ and $`c_5`$) from the the full rotation rate (eq. 2). It is clear that the residuals are significant and the time-variations in the layers immediately below the surface are well above the noise level. It also appears that the magnitude of $`\delta \mathrm{\Omega }`$ increases with solar activity and the peaks in the latitudinal profile tend to shift towards the equator with time.

The zonal flow pattern is more clearly seen in Fig. 2, which shows the contours of constant $`\delta \mathrm{\Omega }`$ at a few selected depths as a function of latitude and time. At any given time the flow pattern has alternating bands in latitude where the rotation rate is either faster or slower than the average. These bands appear to move towards the equator, and at some time after the solar activity minimum the equatorial region changes from faster to slower rotation. This pattern is similar to the well known torsional oscillations observed at the solar surface (Howard & LaBonte 1980). A similar pattern has been seen in other global helioseismic analyses too (Schou 1999; Howe et al. 2000a; Toomre et al. 2000), and also in the near surface layers studied using local-helioseismology methods (Basu & Antia 2000b). There is a good agreement between different independent results obtained in the outer layers of the Sun.

Unfortunately, the prescription for isolating the time-dependent component from the full rotation rate is not unique and in many cases the mean of the entire rotation rate (not just the smooth part) is subtracted to obtain the residual. A problem with this prescription is that the mean depends on the time duration over which the averaging is done, unless the data covers a very long period. On the other hand, the smooth component of rotation rate shows much less temporal variation and hence the mean of the smooth part is relatively insensitive to duration of data. The variation of the smooth component of the flow with time can be seen in Fig. 3. This figure shows the residuals obtained by subtracting the time-averaged smooth component from the smooth component at each epoch. The time variation in this case is much smaller than the variation in the full rotation rate. For example, at low latitudes the full rotation rate varies by about $`\pm 2`$ nHz, while the smooth component varies by only $`\pm 0.5`$ nHz. This clearly shows that most of the temporal variation is in the higher-order splitting coefficients. However, the temporal variation seen in Fig. 3 is still larger than the error estimates and hence is statistically significant. Of course, the resulting pattern is different from that seen in Fig. 2 with full rotation profile as smooth component by definition cannot yield more than two zeros in each hemisphere. There will also be significant distortion due to the fact that the mean over a short period of 4 years is subtracted to get the residuals. The variation seen in Fig. 3 is consistent with that seen in smooth component of surface rotation rate (LaBonte & Howard 1982) in the latitude range shown in Fig. 3.

In order to test the sensitivity of results to varying prescriptions for isolating the time-dependent component, we also compute the residuals where the mean of the full rotation rate is subtracted and the results at $`r=0.98R_{}`$ are shown in Fig. 4. These results can be compared with the upper left hand panel in Fig. 2. It can be seen that although the general pattern is similar in both cases, there are significant differences in the details. In Fig. 2, which shows the results obtained by subtracting only the mean of the smooth component of rotation, the amplitude of residual appears to increase with activity, whereas this increase is not seen in the pattern when the mean of the full rotation rate is subtracted. In order to understand the difference between the two prescriptions for calculating zonal flow component of rotation rate, we show in Fig. 5 the difference in temporal averages of the full rotation rate and that of the smooth component at a few selected depths. This difference will have to be added/subtracted to convert from one definition to the other. It is clear that since the difference itself shows latitudinal variations with amplitude comparable to the zonal flow pattern, there will be a significant difference between the residuals when different averages are subtracted. Fig. 4 also compares the results at $`r=0.98R_{}`$ obtained using different inversion techniques. We find that the zonal flow pattern obtained using 2d RLS inversions agree reasonably well with those obtained with 1.5d RLS, though there are differences at high latitudes. Figs. 8,9 also compare the 1.5d and 2d RLS results and once again it is clear that there is reasonable agreement between the two sets of results. Apart from this we have also varied the regularization used in each inversion technique to check the sensitivity of results to choice of smoothing. It is found that results are not particularly sensitive to the choice of smoothing used in the inversions. Similarly, the results are also insensitive to the choice of modes (frequency range) used for the inversions.

As mentioned earlier, from Figs. 1 and 2, we find that the amplitude of the zonal flow increases with solar activity. In fact, the flow pattern is almost absent near the solar activity minimum. In order to study this variation, in Fig. 6 we show the average magnitude of this component over the latitude range of 0–60 as a function of time. This figure also shows the mean radio flux at 10.7 cm during the time interval covered by each data set as obtained from Solar Geophysical data web page (www.ngdc.noaa.gov/stp/stp.html) of the US National Geophysical Data Center. The 10.7 cm flux from the Sun is known to track solar activity. It is clear that there is a reasonable correlation between the mean magnitude of the zonal flow in layers immediately below the solar surface and the 10.7 cm flux. A similar correlation is found if we use the maximum magnitude of zonal flow velocity instead of the mean magnitude, or when other solar activity indices like the daily sunspot number are used instead of the 10.7 cm flux. This correlation is not particularly sensitive to the range of latitude used in averaging. Further, the magnitude of mean zonal flow velocity decreases with depth. Similar correlation between the mean magnitude of zonal flow velocity and solar activity has been seen in results obtained using ring diagram technique (Basu & Antia 2000b).

At depths exceeding $`0.1R_{}`$ the correlation between the zonal flow and the 10.7 cm flux becomes markedly weaker and there appears to be no correlation at a depth of $`0.2R_{}`$. To show this result quantitatively, in Fig. 7 we plot the correlation coefficient between the average 10.7 cm flux over the period observations were made and the latitudinally averaged magnitude of $`\delta \mathrm{\Omega }`$ as a function of radius. This figure shows that the correlation between the zonal flow and solar activity does not extend to the bottom of the convection zone. It probably also implies that the flow pattern itself does not penetrate through the convection zone. Similar conclusions can be drawn from Fig. 2, where the flow pattern is seen to change significantly below a radius of $`0.9R_{}`$. The correlation coefficient shows a sudden drop around the radius of $`0.89R_{}`$ and this could be considered to be the depth to which flow pattern penetrates. There is some difficulty in estimating the penetration-depth of the zonal flows to high precision because the errors in inversion increase with depth and it becomes increasingly difficult to detect any pattern by just looking at the results at different depths. But the sudden drop in the correlation coefficient around $`0.89R_{}`$ is difficult to explain as being caused by an increase in the errors. Neither the errors, nor the resolution drops suddenly below $`0.89R_{}`$ and hence we would expect a slower drop due to these factors. We can thus conclude that the zonal flow pattern penetrates to a depth of about $`0.1R_{}`$.

As mentioned earlier, if instead of subtracting the mean rotation rate obtained from only the coefficients $`c_1,c_3,c_5`$, we subtract the mean of full rotation rate, the amplitude of residuals doesn’t show any increase with activity. This can also be seen in Fig. 7 which also shows the correlation coefficient for this case. This is due to the fact that the higher order coefficients $`c_7,c_9,\mathrm{}`$ show an increase in amplitude with activity.

In order to study the radial behavior of the zonal flow, we have plotted the flow in Fig. 8 as a function of radius and time at a few selected latitudes. It can be seen that the pattern is quite different below about $`r=0.9R_{}`$. This is just another indication that the flow pattern persists only in the outer layers of the Sun. At all latitudes there appears to be a tendency of the bands of faster or slower rotation to move up in radius with time. Though the significance of this upward movement is not clear as it is not seen at all times. Around the tachocline and below that there appear to be some pronounced oscillations with time at all latitudes. It is not clear if these are real or arise due to the regularization used in the inversion.

The oscillations in the region around the tachocline have also been studied by Howe et al. (2000b) who claim to see a 1.3 year periodicity in the equatorial rotation rate at a depth of $`0.72R_{}`$. In order to verify this periodicity, in Figs. 9 and 10 we show the residuals in the rotation rate at a few selected depths and latitudes as a function of time. For clarity the error-bars are shown only for the 1.5d RLS inversions. The estimated errors in 2d inversions are comparable to those in 1.5d RLS at $`r=0.95R_{}`$ and $`0.90R_{}`$, while in deeper layers (Fig. 10) the error-bars on the results of the 2d inversions are about a factor of 2 lower than those of the 1.5d inversions. In order to make a proper comparison with Howe et al. (2000b) these residuals have been obtained by subtracting the time average of the full rotation rate from the rotation rate obtained from each data set. It is quite clear that there is a systematic variation with time in the rotation rate of the outer layers (which is clearly seen in Fig. 2 also though that was obtained by subtracting just the smooth part of the mean rotation rate). In deeper layers there are some variations which do not seem to be systematic. In order to test whether there are any periodic variations, we take a discrete Fourier transform of the data for the equator at $`r=0.72R_{}`$. The resulting amplitude is shown in Fig. 11. There is no prominent peak in this spectrum at a period of 1.3 yr, or at any other period. Similar spectrum is obtained for other latitudes and depths in this region. Some of the apparent variations seen in the region around the tachocline can arise from the fact that in this region the rotation rate changes very rapidly with radius. Inversions tend to yield some average of the rotation rate in the neighborhood, and furthermore, since each data set has slightly different error estimates and mode set, the averaging would be different and one can expect to get some spurious variations.

## 4 Conclusions

Inversions to determine the solar rotation rate with data from GONG months 1–42, which mostly cover the rising phase of solar cycle 23, show a clear change in the rotation rate with solar activity. The zonal flow, which we have defined as residual in the rotation rate obtained by subtracting a time-average of a three-term fit to the rotation rate, shows bands of faster and slower rotation. These bands appear to move towards the equator with time. The magnitude of this time varying component is about 2 nHz. This pattern is similar to that found by Howe et al. (2000c) using helioseismic data from GONG and MDI projects. The amplitude of the flow pattern in our results is also similar to those in Howe et al. (2000c) as can be seen by comparing their Fig. 3 with our Fig. 9.

The mean magnitude of the pattern in the outer parts of the solar convection zone appears to increase with solar activity and is reasonably well correlated with the mean 10.7 cm flux which tracks solar activity. The correlation is seen up to a depth of about $`0.1R_{}`$ (70 Mm). This appears to be depth to which the zonal flow pattern penetrates. Below this depth the temporal variations in the rotation rate do not appear to be significant. This inferred depth of zonal flow pattern is consistent with estimate of at least 60 Mm obtained by Howe et al. (2000c) using independent data sets and inversion techniques. Clearly, there is good agreement between different inversion results in the outer convection zone. There appears to be some hint that the flow pattern moves upwards with time, however, the statistical significance of this is as yet unclear.

We do not find any evidence of periodic changes in the solar rotation rate near the base of the solar convection zone. We find that the temporal variations at this depth are random in nature and probably represent statistical fluctuations. There is no clear periodic signal in the residual rotation rate with period less than 2 years at any depth. The present data length of 4 years does not allow us to study periodicities with longer time period. Some of this fluctuations may also arise from systematic errors in data.

This work utilizes data obtained by the Global Oscillation Network Group (GONG) project, managed by the National Solar Observatory, a Division of the National Optical Astronomy Observatories, which is operated by AURA, Inc. under a cooperative agreement with the National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofisico de Canarias, and Cerro Tololo Interamerican Observatory.