text:
# The Discovery of Time Variability in OH Megamasers

## 1 Introduction

The OH megamaser (OHM) in IRAS 21272+2514 was discovered in August 1999 and reobserved in October 1999 as part of the Arecibo111The Arecibo Observatory is part of the National Astronomy and Ionosphere Center, which is operated by Cornell University under a cooperative agreement with the National Science Foundation. OHM survey (Darling & Giovanelli 2000, 2001, 2002a). Spectra of this strong megamaser spanning 85 days showed a remarkable change in some, but not all, spectral features. This is the first documented detection of variability in an OHM, and it suggests that many OHMs may exhibit variability on time scales of days to months. Observations of variability can segregate maser emission into compact (variable) and extended (quiescent) emission in a similar and complementary manner to very long baseline interferometry (VLBI) observations which show that roughly two-thirds of OHM emission is compact (Lonsdale et al., 1998; Diamond et al., 1999; Pihlström et al., 2001).

IRAS 21272+2514 has a heliocentric optical redshift $`z=0.15102\pm 0.00012`$, while the bulk of the multi-component OH emission is blueshifted with respect to the optical, spanning 849 km s-1 in the rest frame at 10$`\%`$ of peak flux density and centered on $`z_{OH}=0.15021\pm 0.00005`$ (Darling & Giovanelli 2000, 2002b). The integrated OH luminosity is among the highest known at $`L_{OH}=4.57\times 10^3L_{},`$222We assume $`H_{}=75`$ km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0.7. and the far-IR (FIR) luminosity is constrained to be in the range $`L_{FIR}=(4.9`$$`7.8)\times 10^{11}L_{}`$ (the 100 $`\mu `$m flux is undetected by IRAS). The 1.4 GHz flux density is modest, at $`4.4\pm 0.5`$ mJy (Condon et al., 1998), and IRAS 21272+2514 has some IR excess but is broadly consistent with the radio-FIR relationship for star forming galaxies (Darling & Giovanelli 2002a). This object has a Seyfert 2 nucleus (Darling & Giovanelli 2002b).

In this letter, we present follow-up observations of IRAS 21272+2514 designed to confirm the observed variability (§2), we discuss the modulation present in the OH spectral line components (§3), we set limits on the size scales of all observed maser components (§4), we constrain accelerations in the OH spectral lines (§5), and we discuss the implications of variability for understanding OHMs and the environments which produce them (§6).

## 2 Observations and Data Reduction

Observations were performed at Arecibo with the L-narrow and L-wide receivers as indicated in Table 1, nodding on- and off-source for 4 minutes each. We record spectra in 1-sec intervals in two circular polarizations to facilitate radio frequency interference (RFI) excision as described in Darling & Giovanelli (2001). A weights spectrum is produced for each epoch which shows the fraction of RFI-free records used to construct each channel in a spectrum. The weights spectrum spanning all epochs used to construct the average spectrum shown in Figure 1 indicates that while some RFI is present in 1$`\%`$–4$`\%`$ of the records, it is not the source of the observed spectral variation. Final spectra are averages of two circular polarizations which do not significantly differ except in channels with significant RFI.

Observations have been consolidated into five epochs to enhance the spectral signal to noise as indicated in Table 1. Effective modified Julian dates (MJD) are time-weighted averages of observing sessions within each epoch. Although the separation of observations within epochs spans 1–8 days, no variability is observed within epochs. Note that intraepoch variability is poorly constrained due to low signal-to-noise ratio.

Figure 1 shows the residual spectra (each epoch minus the average spectrum) from the five consolidated epochs, the average spectrum over all time, and the weights spectrum of the average. The main emission components are labeled in the average spectrum. The rms noise is 0.17 mJy in the average spectrum and typically 0.35 mJy in the residuals. The residual spectra show significant modulation only in emission line components 1 and 2. Spectra have been scaled to account for broadband changes in sensitivity and calibration as described in §3.

## 3 Spectral Line Modulation

Comparing measurements from different observing epochs can potentially introduce significant spurious time-dependent signals. The observations presented here bracket significant changes in the performance of the telescope, including a gross and fine refiguring of the primary surface, implementation of active control of the pointing with tie-downs anchoring the platform, refined pointing models, and general fine-tuning of the system. Observations were also performed with two different feeds (see Table 1). The solution to the significant changes between observing epochs lies in accurate calibration of the gain for each observing session, in the broadband nature of the changes in telescope performance, and in the complicated, multi-component spectrum of IRAS 21272+2514. Measuring relative changes between spectral features or changes of a given line as a fraction of the whole spectrum can resolve broadband calibration issues. Sources of spurious narrow-band variability (RFI) have been identified by their short timescales and excised.

Measurements of relative variations between spectral features or with respect to the total OH spectrum is complicated by the possibility that all spectral features may vary, and that the variable components make up a significant fraction of the total OH flux. We approach this problem iteratively, first identifying both the variable and the apparently quiescent features, then comparing the fluxes of quiescent features to the total quiescent flux over time to constrain the modulation actually present in each feature. Once stable spectral features have been identified, the variable features are measured with respect to the total quiescent flux to eliminate broadband instrumental variability between epochs.

Figure 2 shows for each spectral component labeled in Figure 1 the integrated flux density, the peak flux density, and the residual change in both scaled by the integrated spectrum minus spectral components 1 and 2. The average peak and integrated flux density of each spectral component is listed in columns (2) and (3) of Table 2. We find that components 1 and 2 show significant modulation of $`10.4\%\pm 1.6\%`$ and $`14.4\%\pm 1.6\%`$, respectively, while the remaining components do not show significant modulation (col.  of Table 2). The source modulation index $`m_s`$ is the rms fluctuation of the integrated line flux density.

A useful diagnostic of the temporal structure of variability is the structure function: $`D(\tau )=[F(t)F(t+\tau )]^2_t`$. The structure functions of lines 3 and 4 are consistent with zero for all 10 time baselines, whereas lines 1 and 2 show large and highly variable values which indicate characteristic timescales less than the shortest sampling interval of 39 days.

## 4 Maser Size Scales

If the observed modulation is intrinsic to the source, then the usual light travel time constraint applies to the amplified continuum source rather than the masing gas if the masing is unsaturated (unless the variability is due to motions in the gas itself). Variation is observed on the shortest time baseline of 39 days, indicating that the variable source is no larger than 0.03 pc or 6.7 $`\mu `$ac. If the masing region is of roughly the same size, it would have a brightness temperature of $`T6\times 10^{13}`$ K. Such a small source would certainly scintillate. Hence, the observed variability is at least partly due to interstellar scintillation (ISS), and we choose the weaker constraint on the angular sizes of the variable masing regions by adopting ISS to interpret the observations.

Walker (1998) provides rough predictions of the scintillation properties of extragalactic sources based on the interstellar medium model of Taylor & Cordes (1993) which we apply to the observed variation in the OH line components of IRAS 21272+2514. The OH line radiation must pass through two scintillating screens, but the dominant influence will be from the electron screen in our own galaxy. The galactic coordinates of 21272+2514, $`(l,b)=(75.7^{},18.4^{})`$, indicate a transition frequency from strong to weak scattering at $`\nu _{}=15\pm 1`$ GHz and a first Fresnel zone at $`\nu _{}`$ of size $`\theta _F=2\mu `$as. At $`\nu =1.45`$ GHz, the scattering should be strong refractive ISS (RISS). For a point source, Walker (1998) predicts a modulation index of $`m_p=27\%`$, a scattering disk 0.34 mas in size, and a scintillation timescale of 10 days. We scale the point source modulation index $`m_p`$ to match the observed modulation $`m_s`$ in each source component to obtain rough estimates of source sizes and scintillation timescales from the scaling relations $`\theta _s/\theta _p=(m_p/m_s)^{6/7}`$ and $`\tau _s/\tau _p=\theta _s/\theta _p`$. Estimates of the angular sizes, RISS timescales, and physical sizes of the OH line components are listed in columns (5)–(7) of Table 2. The two components which show strong modulation have submilliarcsecond angular sizes, are less than 2 pc in size, and are expected to show a characteristic RISS timescale of 17–22 days. This is in good agreement with the observations and our inability to observe changes on short time baselines of 1–8 days. This treatment assumes that all of the emission within a variable spectral component originates from the same compact region. It is also possible that part of the emission is compact and part is extended, so the scintillation predictions are actually upper limits on the source size and the RISS timescales may be as low as 10 days. If the compact variable features are smaller than the scattering disk, then they are expected to show the full point source modulation of 27$`\%`$ and must account for at least 39$`\%`$ and 53$`\%`$ of the emission seen in components 1 and 2, respectively. In other words, no more than roughly half of the emission in components 1 and 2 can be extended. These constraints on the fraction of compact emission present in OHMs are consistent with VLBI observations of nearby OHMs by Diamond et al. (1999) and others who find that roughly two-thirds of the total OH emission is compact.

Components 3–5, which do not show significant modulation, are constrained to be larger than 1 mas or a few parsecs. RISS timescales in these components are constrained by their modulation indices to be larger than 28–51 days. Note that although observations span longer time baselines, the observed limits on modulation provide the stronger constraint — modulation may be present in components 3–5 on timescales of a few months, but the observations lack the sensitivity to detect variations of less than about 4$`\%`$. The fraction of compact emission in components 3, 4, and 5 is less than 15$`\%`$, 16$`\%`$, and 30$`\%`$, respectively.

## 5 Acceleration

Like water masers associated with AGN, one might expect OH megamasers to exhibit accelerations in emission line features. The kinematical structure of major mergers is certainly complicated and dynamic, and it is possible that accelerating masing gas will be observable over the timescales of this study. Water maser features in the best-studied case, NGC 4258, show accelerations of 9.3 km s-1 yr-1 along the line of sight (Herrnstein et al., 1999). Accelerations of this magnitude in OHM line features can be detected by this study if they exist. One might expect that since OHM emission profiles are thought to be an ensemble of many masing regions, then accelerations will be averaged out by the summation. This is especially true if OHMs are produced in highly turbulent or chaotic accumulations of molecular gas. The wide range of line profiles suggests that this may be the case, although there are a number of narrow, physically compact emission lines in the sample. It might be surprising if the lines showing intensity fluctuations do not show significant accelerations.

Overall, frequency/velocity residuals are consistent with the uncertainties in the line peaks of $`\pm 1`$ channel, corresponding to 24 kHz in the observer frame. The mean frequency of the peak of each spectral component is listed in column (8) of Table 2. The associated uncertainties are generally consistent with $`\pm 24/\sqrt{N}`$ kHz where $`N`$ is the number of epochs used in the mean ($`N`$ is listed in parentheses at the end of col.  of Table 2). Acceleration can be measured from the slope of a linear fit to the time series of line peak frequencies and expressed in the rest frame of the source by including a factor $`(1+z)^2`$ to account for time dilation and the narrowing of frequency intervals with redshift. The minimum detectable acceleration observable in an OHM line feature is set by the uncertainty of the best-fit slope obtained from $`N`$ observations with velocity resolution $`\pm 4.4`$ km s-1(observer frame):

$$\sigma _a=\frac{\sigma _v(1+z)}{\sigma _t\sqrt{N1}}=(4.4\text{ km s}\text{-1}\text{ yr}\text{-1})\frac{(1+z)^2}{\sigma _t(\text{yr})\sqrt{N1}}$$
(1)
where $`\sigma _t^2`$ is the variance of the sampling in time and $`\sigma _v`$ is the uncertainty of $`\pm 1`$ channel in identifying the peak of a line. For 21272+2514, with five epochs sampled, $`\sigma _a=2.8`$ km s-1 yr-1. This sets a rough lower limit on the line accelerations which may be measured from the available data.

Table 2 lists the rest frame line accelerations for each component from fits to the time series of line peaks. Uncertainties are obtained from the fits and $`N`$ is listed in parentheses. Except component 2 where the peak could not be identified in one epoch, accelerations are obtained from fits to all data points. Line accelerations are consistent with zero except for line 3 which shows a 3.3 $`\sigma `$ deviation from zero. The acceleration of $`2.6\pm 0.8`$ km s-1 yr-1 for component 3 is very close to the expected lower detection limit of 2.8 km s-1 yr-1. Confirmation of this acceleration requires re-observation, which will provide a longer time baseline. It is worth mention that the peaks of line 3 show a steady decrease in frequency over time with little scatter and no outliers. It is interesting that this line also shows no significant modulation in amplitude.

## 6 Conclusions

This brief study of IRAS 21272+2514 demonstrates the powerful constraints that variability studies can place on the emission regions of OHMs, including the segregation of different spectral components into different physical size scales. Better sampling of intermediate time scales can provide even stronger constraints on the source properties including the physical environments responsible for OHMs. There is mounting evidence that OHMs are composed of both diffuse and compact emission regions with different saturation states (Lonsdale et al., 1998; Diamond et al., 1999). This study obtains constraints on both the size and proportion of compact emission regions: compact features must be smaller than about 2 pc and account for 27$`\%`$–58$`\%`$ of the total OH emission. To account for such high surface brightness emission we suggest that either shocks are involved in the pumping of OH or that these compact regions represent fortuitous geometries which provide enhanced masing path lengths as suggested by Pihlström et al. (2001). Accelerations of the OH lines of IRAS 21272+2514 are generally constrained to be less than 3 km s-1 yr-1, although further observations are required to resolve suggestions of acceleration in one component.

Many thanks indeed to B. Catinella for obtaining the spectra of IRAS 21272+2514 during the MJD 52185 and 52224 epochs and to M. Walker for supplemental scintillation calculations. It is a pleasure to thank the staff of the Arecibo Observatory for excellent support. This work has received support from STScI grant 8373 and NSF grant AST 00-98526.


format: text

text:
# Anomalous scaling of a passive scalar advected by the turbulent velocity field with finite correlation time: Two-loop approximation

## I Introduction

In recent years, considerable progress has been achieved in the understanding of intermittency and anomalous scaling of fluid turbulence. The crucial role in these studies was played by a simple model of a passive scalar quantity advected by a random Gaussian field, white in time and self-similar in space, the so-called Kraichnan’s rapid-change model . There, for the first time the existence of anomalous scaling was established on the basis of a microscopic model  and the corresponding anomalous exponents were calculated within controlled approximations  and a systematic perturbation expansion in a formal small parameter . Detailed review of the recent theoretical research on the passive scalar problem and the bibliography can be found in Ref. .

Within the approach developed in Refs. , nontrivial anomalous exponents are related to “zero modes,” that is, homogeneous solutions of the closed exact differential equations satisfied by the equal-time correlation functions. In this sense, the rapid-change model appears “exactly solvable.”

In a wider context, zero modes can be interpreted as statistical conservation laws of the particle dynamics . The concept of statistical conservation laws appears rather general, being also confirmed in numerical simulations by Refs. , where the passive advection in the two-dimensional Navier–Stokes velocity field  and a shell model of a passive scalar  were studied. This observation is rather intriguing because in those models no closed equations for equal-time quantities can be derived due to the fact that the advecting velocity has a finite correlation time (for a passive field advected by a velocity with given statistics, closed equations can be derived only for different-time correlation functions, and they involve infinite diagrammatic series).

One may thus conclude that breaking the artificial assumption of the time decorrelation of the velocity field is the crucial point .

An important issue related to the effects of the finite correlation time is the universality of the anomalous exponents. It was argued that the exponents may depend on more details of the velocity statistics than only the exponents $`\eta `$ and $`\epsilon `$ . This idea was supported in Refs. , where the case of short but finite correlation time was considered for the special case of a local turnover exponent. In those studies, the anomalous exponents were derived to first order in small correlation time, with Kraichnan’s rapid-change model  or analogous shell model for a scalar field  taken as zeroth order approximations. The exponents obtained appear nonuniversal through the dependence on the correlation time.

In Ref.  and subsequent papers , the field theoretic renormalization group (RG) and operator product expansion (OPE) were applied to the rapid-change model and its descendants. In that approach, anomalous scaling emerges as a consequence of the existence in the model of composite operators with negative scaling dimensions, identified with the anomalous exponents. This allows one to construct a systematic perturbation expansion for the anomalous exponents, analogous to the famous $`\epsilon `$ expansion in the RG theory of critical behavior, and to calculate the exponents to the second  and third  orders. For passively advected vector fields, any calculation of the exponents for higher-order correlations calls for the RG techniques already in the lowest-order approximation .

Besides the calculational efficiency, an important advantage of the RG approach is its relative universality: it is not related to the aforementioned solvability of the rapid-change model and can also be applied to the case of finite correlation time or non-Gaussian advecting field. In Ref.  (see also  for the case of compressible flow) the RG and OPE were applied to the problem of a passive scalar advected by a Gaussian self-similar velocity with finite (and not small) correlation time. The energy spectrum of the velocity in the inertial range has the form $`E(k)k^{12\epsilon }`$, while the correlation time at the wavenumber $`k`$ scales as $`k^{2+\eta }`$. It was shown that, depending on the values of the exponents $`\epsilon `$ and $`\eta `$, the model reveals various types of inertial-range scaling regimes with nontrivial anomalous exponents. For $`\eta >\epsilon `$, they coincide with the exponents of the rapid-change model and depend on the only parameter $`2\epsilon \eta `$, while for $`\epsilon >\eta `$ they coincide with the exponents of the opposite (“quenched” or “frozen”) case and depend only on $`\epsilon `$.

The most interesting case is $`\eta =\epsilon `$, when the exponents can be nonuniversal through the dependence on the correlation time (more precisely, on the ratio $`u`$ of the velocity correlation time and the turnover time of the passive scalar). In the field theoretic language, the nonuniversality of the exponents in this regime is a consequence of the degeneracy of the corresponding fixed point of the RG equations. It agrees with the findings of Refs.  since the borderline $`\eta =\epsilon `$, including the “Kolmogorov” point $`\eta =\epsilon =4/3`$, corresponds to the case of a local turnover exponent. It is also interesting to note that the same relation $`\eta =\epsilon `$ for the boundary between the time-decorrelated and quenched cases is encountered in a model of passive advection by a strongly anisotropic flow, studied in Refs. . It was argued in  that the same boundary will be observed with very general assumptions on the velocity statistics. Although the possibility of the nonuniversality of anomalous exponents for $`\eta =\epsilon `$ was demonstrated by the rigorous RG analysis, the practical calculation by Ref.  has shown that they appear universal (independent of $`u`$) to the first order in $`\epsilon `$ and $`\eta `$: in the one-loop approximation, the anomalous dimensions of the relevant composite operators depend on a combination of the model parameters (couplings) that remains constant along the line of the fixed points. This fact is rather disappointing, because it means that in the one-loop approximation it is impossible to judge how the finite correlation time affects intermittency, in particular, whether the anomalous scaling is enhanced or suppressed in comparison with the rapid-change or quenched limits.

In this paper, we present the anomalous exponents to order $`O(\epsilon ^2)`$ (two-loop approximation) for the most interesting case $`\eta =\epsilon `$, including the exponents of the anisotropic contributions, and study their dependence on $`u`$. \[It is not necessary to separately consider the cases $`\eta >\epsilon `$ ($`\eta <\epsilon `$), because the corresponding exponents are the same as for the rapid-change (quenched) velocity field and can be obtained from the case $`\eta =\epsilon `$ in the limits $`u\mathrm{}`$ ($`u0`$)\].

In Sec. II, we describe our model and its interesting special cases. In Sec. III, we briefly recall the field theoretic formulation, the RG and OPE approach to the model, and the $`O(\epsilon )`$ result for the anomalous exponents . The results of the two-loop calculation are presented and discussed in Sec. IV: in Sec. IV A, we give the anomalous exponents to order $`O(\epsilon ^2)`$ and then discuss them separately for the isotropic (Sec. IV B) and anisotropic (Sec. IV C) contributions.

The main conclusion of the paper can be formulated as follows: the qualitative effect of the finite correlation time on the anomalous scaling depends essentially on the correlation function considered, the value of $`u`$, and the space dimensionality $`d`$. For the low-order structure functions and in low dimensions ($`d=2`$ or 3), the inclusion of finite correlation time enhances the intermittency in comparison with the both limits: the time-decorrelated ($`u=\mathrm{}`$) and time-independent ($`u=0`$) ones. Although the anomalous exponents have a well-defined limits for $`u0`$, they show interesting irregularities in the vicinity of the quenched limit: a rapid falloff when $`u=0`$ increases from zero, with infinite slope for $`d=2`$, with a pronounced minimum for $`u1`$. On the contrary, the behavior in the region of large $`u`$ is smooth, like for the shell model studied in . For higher-order structure functions and large $`d`$, the anomalous scaling is always weaker in comparison with the rapid-change limit and the corresponding (positive) correction is maximal for $`u=0`$ and monotonically decreases to zero as $`u`$ tends to infinity.

## II The model

The advection of a passive scalar field $`\theta (x)\theta (t,𝐱)`$ is described by the stochastic equation

$$_t\theta =\nu _0^2\theta +f,_t_t+v_i_i,$$
(1)
where $`_t/t`$, $`_i/x_i`$, $`\nu _0`$ is the molecular diffusivity coefficient, $`^2`$ is the Laplace operator, $`𝐯(x)\{v_i(x)\}`$ is the divergence-free (owing to the incompressibility) velocity field, and $`ff(x)`$ is an artificial Gaussian random noise with zero mean and correlation function

$$f(x)f(x^{})=C(tt^{},𝐫),𝐫=𝐱𝐱^{}.$$
(2)
The form of the correlator is unessential; it is only important that the function $`C`$ in Eq. (2) decreases rapidly for $`rL`$, where $`L`$ is some integral scale. The noise maintains the steady state of the system and, if $`C`$ depends on the vector $`𝐫`$ and not only its modulus $`r𝐫`$, is a source of large-scale anisotropy. In a more realistic formulation, the noise is replaced by an imposed constant gradient of the scalar field; see e.g. Refs. .

In the real problem, the velocity field satisfies the Navier–Stokes equation. Following Refs. , we assume for $`𝐯(x)`$ in Eq. (1) a Gaussian distribution with zero mean and correlator

$$v_i(x)v_j(x^{})=\frac{d𝐤}{(2\pi )^d}P_{ij}(𝐤)D_v(tt^{},k)\mathrm{exp}[\mathrm{i}𝐤(𝐱𝐱^{})],$$
(3)
where $`P_{ij}(𝐤)\delta _{ij}k_ik_j/k^2`$ is the transverse projector and the function $`D_v(t,k)`$ will be chosen in the form

$$D_v(tt^{},k)=\frac{D_0}{2u_0}\frac{1}{k^{d2+2\epsilon }}\mathrm{exp}[\omega (k)(tt^{})],\omega (k)=u_0\nu _0k^{2\eta }.$$
(4)
Here $`D_0`$ and $`u_0`$ are positive amplitude factors and the positive exponents $`\epsilon `$ and $`\eta `$ play the part of small expansion parameters in the RG theory; see Refs. . It is also convenient to introduce the “coupling constant” $`g_0D_0/\nu _0^2`$ (expansion parameter in the ordinary perturbation theory). The infrared (IR) regularization is provided by the sharp cutoff in all momentum intergals from below at $`k=m`$ with $`m1/L`$.

The model (3), (4) contains two special cases that possess some interest on their own: in the limit $`u_0\mathrm{}`$, $`g_0^{}g_0/u_0^2=\mathrm{const}`$ we arrive at the rapid-change model,

$$D_v(\omega ,k)g_0^{}\nu _0\delta (tt^{})k^{d2\epsilon +\eta },$$
(5)
while the limit $`u_00`$, $`g_0^{\prime \prime }g_0/2u_0=\mathrm{const}`$ corresponds to the case of a quenched (time-independent) velocity field,

$$D_v(\omega ,k)g_0^{\prime \prime }\nu _0^2k^{d+22\epsilon }.$$
(6)
The latter case has a close formal resemblance with the well-known models of random walks in random environment with long-range correlations; see Refs. .

## III Renormalization group and operator product expansion

The RG theory of the model (1)–(4) is presented in Refs.  in detail; below we briefly recall only the necessary information. The stochastic problem (1)–(4) can be cast as a field theory with action functional

$$S(\theta ,\theta ^{},v)=𝐯D_v^1𝐯/2+\theta ^{}D_f\theta ^{}/2+\theta ^{}\left[_t+\nu _0^2\right]\theta ,$$
(7)
where $`\theta ^{}`$ is an auxiliary scalar field and $`D_f`$ and $`D_v`$ are correlators (2) and (3), respectively. In Eq. (7), all the required integrations over $`x=(t,𝐱)`$ and summations over the vector indices are understood.

The model (7) is logarithmic for $`\epsilon =\eta =0`$; the ultraviolet (UV) singularities have the form of poles in various linear combinations of $`\epsilon `$ and $`\eta `$ in the correlation functions. They can be removed by the only counterterm of the form $`\theta ^{}^2\theta `$, which is equivalent to the following multiplicative renormalization of the parameters $`g_0`$, $`u_0`$ and $`\nu _0`$ in the action functional (7):

$$\nu _0=\nu Z_\nu ,g_0=g\mu ^{2\epsilon +\eta }Z_g,u_0=u\mu ^\eta Z_u,$$
(8)
where $`g`$, $`u`$, and $`\nu `$ are the renormalized counterparts of the bare parameters, $`\mu `$ is the reference mass in the minimal subtraction (MS) scheme, which we always use in practical calculations, and $`Z_i=Z_i(g,u;d;\epsilon ,\eta )`$ are the renormalization constants satisfying the identities

$$Z_g=Z_\nu ^3,Z_u=Z_\nu ^1.$$
(9)
Fixed points of the corresponding RG equations are found from the requirement that the $`\beta `$ functions

$$\beta _g\stackrel{~}{𝒟}_\mu g=g[2\epsilon \eta +3\gamma _\nu ],\beta _u\stackrel{~}{𝒟}_\mu u=u[\eta +\gamma _\nu ],\gamma _\nu \stackrel{~}{𝒟}_\mu \mathrm{ln}Z_\nu $$
(10)
vanish. Here $`\stackrel{~}{𝒟}_\mu `$ is the operation $`\mu _\mu `$ for fixed bare parameters and the relations between $`\beta `$ functions and the anomalous dimension $`\gamma _\nu `$ result from the definitions and the relation (9).

From Eq. (10) the exact relation $`\beta _g/g3\beta _u/u=2(\eta \epsilon )`$ follows which shows that the $`\beta `$ functions cannot vanish simultaneously for finite values of their arguments, except for the case $`\eta =\epsilon `$. Therefore, to find the fixed points we must set either $`u=\mathrm{}`$ or $`u=0`$ and simultaneously rescale $`g`$ so that the anomalous dimension $`\gamma _\nu `$ remain finite. These two options correspond to the two limits (5) and (6), so that the rapid-change and quenched cases are fixed points of the general model. The analysis shows that the former is IR stable (and thus describes the inertial-range asymptotic behavior) for $`\eta >\epsilon `$ while the latter is stable for $`\eta <\epsilon `$.

Most interesting is the case $`\eta =\epsilon `$, when the $`\beta `$ functions become proportional, and the set $`\beta _g=\beta _u=0`$ reduces to a single equation. As a result, the corresponding fixed point is degenerate: rather than a point, one obtains a line of fixed points in the $`g`$$`u`$ plane. They can be labelled by the value of the parameter $`u`$, which has the meaning of the ratio of the velocity correlation time and the scalar turnover time.

Existence of the IR stable fixed points implies certain scaling properties of various correlation functions at scales larger than the dissipative length $`g_0^{1/3\epsilon }`$. In particular, for the equal-time structure functions

$$𝒮_n(𝐫)=[\theta (t,𝐱)\theta (t,𝐱^{})]^n,𝐫=𝐱𝐱^{}$$
(11)
one obtains

$$𝒮_n(𝐫)=D_0^{n/2}r^{n(1\epsilon /2)}F_n(m𝐫)$$
(12)
(odd structure functions are nontrivial if the correlation function $`vf`$ is nonzero or if a constant gradient of the scalar field is imposed). In the presence of anisotropy the scaling functions $`F_n(m𝐫)`$ can be decomposed in irreducible representations of the $`SO(d)`$ group. In the simplest case of uniaxial anisotropy (which is sufficient to reveal all anomalous exponents) one can write

$$F_n(m𝐫)=P_l(z)F_{nl}(mr),z=(𝐧𝐫)/r,$$
(13)
where $`P_l(z)`$ is the $`l`$th order Gegenbauer polynomial (Legendre polynomial for $`d=3`$) and $`𝐧`$ is a unit vector that determines the distinguished direction.

The leading behavior of the functions $`F_{nl}`$ for $`mr0`$ (inertial range) is found from the corresponding operator product expansion and has the form

$$F_{nl}(mr)^{\mathrm{\Delta }_{nl}},$$
(14)
where the “anomalous exponent” $`\mathrm{\Delta }_{nl}`$ is nothing other than the critical dimension of the irreducible traceless $`l`$th rank tensor composite operator built of $`n`$ fields $`\theta `$ and minimal possible number of derivatives . For $`ln`$ such an operator has the form

$$_{i_1}\theta \mathrm{}_{i_l}\theta (_i\theta _i\theta )^p+\mathrm{},n=l+2p.$$
(15)
Here the dots stand for the appropriate subtractions involving the Kronecker $`\delta `$ symbols, which ensure that the resulting expressions are traceless with respect to any given pair of indices, for example, $`_i\theta _j\theta \delta _{ij}_k\theta _k\theta /d`$. We also note that the numbers $`n`$ and $`l`$ necessarily have the same parity, that is, they can only be simultaneously even or odd.

For the most interesting case of the degenerate fixed point, the dimensions $`\mathrm{\Delta }_{nl}`$ are calculated in the form of series in the only independent exponent $`\epsilon =\eta `$, that is,

$$\mathrm{\Delta }_{nl}=\underset{k=1}{\overset{\mathrm{}}{}}\epsilon ^k\mathrm{\Delta }_{nl}^{(k)}.$$
(16)
In the lowest order one obtains :

$$\mathrm{\Delta }_{nl}^{(1)}=\frac{n(n2)(d1)+\lambda _l(d+1)}{2(d1)(d+2)}$$
(17)
with $`\lambda _ll(d+l2)`$. For $`k2`$, the coefficients $`\mathrm{\Delta }_{nl}^{(k)}`$ depend not only on $`d`$ but also on the parameter $`u`$, the ratio of the velocity correlation time and the scalar turnover time, that labels fixed points in the $`g`$$`u`$ plane (see above).

The reader not interested in the details of practical calculation can skip the end of this section and pass to the result for $`\mathrm{\Delta }_{nl}^{(2)}`$. Calculation of the higher-order coefficients in the $`\epsilon `$ expansions for the rapid-change model is presented in Refs.  in detail. Analogous calculations for the finite correlated case are more difficult in two respects. First, there are more relevant Feynman diagrams in the same order of perturbation theory (for zero correlation time, many diagrams contain closed circuits of retarded propagators and vanish). Second, and more important distinction, is that the diagrams for the finite correlated case involve two different dispersion laws: $`\omega k^2`$ for the scalar and $`\omega k^{2\eta }`$ for the velocity fields. As a result, the calculation, as well as expressions for the renormalization constants, become rather cumbersome already in the lowest (one-loop) approximation; see Refs. .

The latter difficulty can be circumvented as follows. Careful analysis shows that in the MS scheme all the needed anomalous dimensions, $`\gamma _\nu `$ from (10) and $`\gamma _{nl}\stackrel{~}{𝒟}_\mu \mathrm{ln}Z_{nl}`$, in contrast with the respective renormalization constants $`Z_\nu `$ and $`Z_{nl}`$, are independent of the exponents $`\epsilon `$ and $`\eta `$ in the two-loop approximation (for the one-loop approximation this is obvious from the explicit expressions; see ). It is thus sufficient to calculate them for any specific choice of the exponents $`\epsilon `$ and $`\eta `$ that guarantees UV finiteness of the diagrams. The most convenient choice is $`\eta =0`$ and arbitrary $`\epsilon `$: all the diagrams remain finite, the exponents in the aforementioned dispersion laws become identical, and the practical calculations drastically simplify and become feasible.

To avoid possible misunderstandings, it should be emphasized that such an independence is not guaranteed by the renormalizability of the model. The renormalizability in the analytic regularization only guarantees that the renormalization scheme can be chosen such that the correlation functions, along with the coefficients $`\beta `$ and $`\gamma `$ in the RG equations, will be analytic at the origin in the space of two complex variables $`\epsilon `$ and $`\eta `$ . We used another scheme in which the functions $`\beta `$ and $`\gamma `$ are independent of $`\epsilon `$ and $`\eta `$ in the first two orders, which does not exclude nonanalytic dependence on these parameters in higher orders. We expect that in the three-loop approximation nonanalytic constructions like $`(\epsilon +\eta )/(\epsilon +2\eta )`$ will indeed appear in the anomalous dimensions, in particular, due to the necessity to take into account UV finite parts of the two-loop diagrams (with our choice of the sharp IR cutoff in Eq. (4), the one-loop diagrams have no UV finite parts; cf.  for the rapid-change case).

Thus, we conclude that the knowledge of the renormalization constants for $`\eta =0`$ is sufficient to obtain the anomalous dimensions, $`\beta `$ functions, coordinates of the fixed points, and the critical dimensions of composite operators for arbitrary values of $`\eta `$ and $`\epsilon `$, including the most interesting case $`\eta =\epsilon `$, which we always discuss from now on.

## IV Anomalous exponents in the two-loop approximation

### A General expressions

We have performed the complete two-loop calculation of the RG functions (10) and the critical dimensions (16) of the composite operators (15) for arbitrary values of $`n`$, $`l`$, $`d`$, and $`u`$ and obtained the following expression for the second coefficient in expansion (16):

$`\mathrm{\Delta }_{nl}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{(d1)^2(d+2)^2(d+4)}}\{2(d+4)𝒜[n(n2)(d1)+\lambda _l]+(n2)\{6[n(n4)(d1)+3\lambda _l]+`$ (18)
$`+`$ $`9𝒞[n(d+n)(d1)\lambda _l(d+1)]\}\}`$ (19)
with $`\lambda _ll(d+l2)`$. Here and below we denote

$`𝒜`$ $`=`$ $`{\displaystyle \frac{(u11/u)(d+1)}{2(d+2)(1u)}}+{\displaystyle \frac{(d+1)}{2(d+4)(1u)u(1+u)^2}}F_3\left({\displaystyle \frac{1}{(u+1)^2}}\right)+{\displaystyle \frac{2ud(d+2)}{(1d)(1u)}}𝒥(u,d),`$ (21)
$``$ $`=`$ $`{\displaystyle \frac{(d+1)}{3(1u)^2(d+4)}}\left[{\displaystyle \frac{u}{u+1}}F_3\left({\displaystyle \frac{1}{2(u+1)}}\right){\displaystyle \frac{1}{(u+1)^2}}F_3\left({\displaystyle \frac{1}{(u+1)^2}}\right){\displaystyle \frac{u^2}{4}}F_3\left({\displaystyle \frac{1}{4}}\right)\right],`$ (22)
$`𝒞`$ $`=`$ $`{\displaystyle \frac{1}{9(1u)^2}}\{{\displaystyle \frac{3u^2(d1)}{4}}F_2\left({\displaystyle \frac{1}{4}}\right){\displaystyle \frac{u[2d1+u(d2)]}{(u+1)}}F_2\left({\displaystyle \frac{1}{2(u+1)}}\right)+{\displaystyle \frac{[d+1+2u(d2)]}{(u+1)^2}}F_2\left({\displaystyle \frac{1}{(u+1)^2}}\right)`$ (23)
$``$ $`{\displaystyle \frac{u^2(d+1)}{(d+4)}}F_3\left({\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{4u(d+1)}{(u+1)(d+4)}}F_3\left({\displaystyle \frac{1}{2(u+1)}}\right){\displaystyle \frac{4(d+1)}{(u+1)^2(d+4)}}F_3\left({\displaystyle \frac{1}{(u+1)^2}}\right)\}.`$ (24)
We also have denoted $`F_k(x)F(1,1;d/2+k;x)`$ for the hypergeometric series

$`F(a,b;c;z)1+{\displaystyle \frac{ab}{c}}z+{\displaystyle \frac{a(a+1)b(b+1)}{c(c+1)}}{\displaystyle \frac{z^2}{2!}}+\mathrm{}`$

The values of $`F_k`$ entering into (IV A) can be related by the recurrent relation $`(x1)F_2(x)=x(d+2)F_3(x)/(d+4)1`$, but the resulting expressions look more cumbersome and we shall keep both $`F_2`$ and $`F_3`$ in the formulas.

The quantity $`𝒥(u,d)`$ in Eq. (21) can only be expressed in the form of a single convergent integral, suitable for numerical calculation:

$`𝒥(u,d)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(d/2)}{\sqrt{\pi }\mathrm{\Gamma }((d1)/2)}}{\displaystyle _0^1}dz{\displaystyle \frac{(1z^2)^{d/2}}{(u1)^2+4uz^2}}\{z^2(1z^2)\mathrm{ln}\left({\displaystyle \frac{1+u}{2}}\right)z(u1+2z^2)\mathrm{arcsin}z`$ (25)
$``$ $`{\displaystyle \frac{z(1z^2)^{1/2}(1uz^2)}{[2(1+u)z^2]^{1/2}}}\mathrm{arctan}\left[{\displaystyle \frac{z[2(1+u)z^2]^{1/2}}{(1+uz^2)}}\right]\},`$ (26)
where $`\mathrm{\Gamma }(\mathrm{})`$ is the Euler Gamma function.

The quantities (IV A), and hence the dimensions (19), have finite limits for $`u\mathrm{}`$ and $`u0`$. In the first limit, $`\mathrm{\Delta }_{nl}^{(2)}`$ coincides with the known result for the Kraichnan’s rapid-change model (see  for $`l=0`$ and 2 and  for general $`l`$). The O(1/u) correction to the rapid-change limit can be found from the following asymptotic expressions for the coefficients (IV A):

$`𝒜`$ $`=`$ $`{\displaystyle \frac{(d+1)}{2(d+2)}}(1+2/u)+O(1/u^2),`$ (28)
$``$ $`=`$ $`{\displaystyle \frac{(d+1)}{12(d+4)}}F_3\left({\displaystyle \frac{1}{4}}\right)(1+2/u)+O(1/u^2),`$ (29)
$`𝒞`$ $`=`$ $`\left[{\displaystyle \frac{(d1)}{12}}F_2\left({\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{(d+1)}{9(d+4)}}F_3\left({\displaystyle \frac{1}{4}}\right)\right]+{\displaystyle \frac{1}{u}}\left[{\displaystyle \frac{(d1)}{6}}F_2\left({\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{2(d+1)}{9(d+4)}}F_3\left({\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{(d2)}{9}}\right]+O(1/u^2).`$ (30)
The opposite case, $`u=0`$, corresponds to the quenched (time-independent) velocity field. This case was extensively studied in connection with the so-called “random-random walks” (random walks in random environments); see the review paper  and references therein. Our results for the function $`\beta _g`$ from (10) and the corresponding fixed point are in agreement with the two-loop results quoted in Ref.  for the model of random-random walks. To our knowledge, the dimensions of composite operators (15) have not been studied in that context, and below we give the asymptotic expressions for the coefficients (IV A) sufficient to find the dimensions (19) up to order $`O(u)`$ near the quenched limit:

$`𝒜`$ $`=`$ $`{\displaystyle \frac{(d+1)(3d+4)}{2d(d+2)}}+u(d+1)\left\{{\displaystyle \frac{(3d+2)}{2d(d2)}}{\displaystyle \frac{2}{d(d+2)}}+{\displaystyle \frac{\mathrm{ln}2}{(d+4)}}+{\displaystyle \frac{(d+3)}{(d+4)(d+6)}}F_4\left({\displaystyle \frac{1}{2}}\right){\displaystyle \frac{2(d+3)}{(d+4)^2}}\right\},`$ (32)
$``$ $`=`$ $`{\displaystyle \frac{(d+1)}{3(d+2)}}+{\displaystyle \frac{u(d+1)}{3}}\left\{{\displaystyle \frac{1}{(d+4)}}F_3\left({\displaystyle \frac{1}{2}}\right)+{\displaystyle \frac{4}{d(d+2)}}\right\},`$ (33)
$`𝒞`$ $`=`$ $`{\displaystyle \frac{(d+1)(d^2+4)}{9d(d+2)}}+{\displaystyle \frac{u}{9}}\left\{(12d)F_2\left({\displaystyle \frac{1}{2}}\right)+{\displaystyle \frac{4(d+1)}{(d+4)}}F_3\left({\displaystyle \frac{1}{2}}\right){\displaystyle \frac{4(d+1)(d+2)}{d(d2)}}+{\displaystyle \frac{2(d^2+2d+4)}{(d+2)}}\right\},`$ (34)
up to corrections of order $`O(u^2)`$.

It is worth noting that the $`O(u)`$ terms in Eqs. (32) and (34) contain poles in $`(d2)`$ and thus diverge for $`d=2`$. Analysis shows that, for $`d=2`$, the leading correction to the $`u=0`$ result is not analytical in $`u`$ and has the form $`u\mathrm{ln}u`$. Formally, the singularity at $`d=2`$ is explained as follows. Some of the two-loop diagrams contain “energy denominators” of the form $`(𝐤+𝐪)^2+O(u)`$, where $`𝐤`$ and $`𝐪`$ are two independent integration momenta. The numerators contain factors $`(P_{ij}(𝐤)q_iq_j)^2`$ stemming from transverse projectors in the propagators. These factors suppress the singularity at $`𝐤=𝐪`$, occurring in the denominators for $`u=0`$, and ensure the existence of the integrals over $`𝐤`$ and $`𝐪`$, but the “collinear” divergence at $`𝐤=𝐪`$ occurs if the $`O(u)`$ to the denominators is taken into account. Physically, this divergence can be related to a strong resonant interaction between the excitations of the passive scalar field with the opposite momenta $`𝐤=𝐪`$ of equal moduli in two dimensions. We shall see below that this singularity remarkably affects the behavior of the dimensions (19) for the values of $`d`$ much larger than $`d=2`$.

Many studies have been devoted to the analysis of the inertial-range turbulence in the limit $`d\mathrm{}`$ . Our model has no finite “upper critical dimension,” above which anomalous scaling would vanish. Like in the rapid-change case  and, probably, in the Navier-Stokes turbulence , the anomalous scaling disappears at $`d=\mathrm{}`$, but it reveals itself already in the $`O(1/d)`$ approximation. Along with the results  for the scalar rapid-change model, where the $`O(1/d)`$ expression for the anomalous exponents were derived for any $`\epsilon `$, this confirms the importance of the large-$`d`$ expansion for the issue of anomalous scaling in fully developed turbulence.

Straightforward analysis of the expressions (IV A) shows that, for $`d\mathrm{}`$, one has $`𝒜=O(d^0)`$ \[it is important here that $`𝒥(u,d)=O(1/d)`$\], $`=O(d^0)`$ and $`𝒞=O(d)`$, namely,

$$𝒞=\frac{(u+2)(3u+2)}{36(u+1)^2}d+O(d^0).$$
(35)
It then follows that for large $`d`$, the dimension (19) behaves as $`O(1/d^2)`$ and is completely determined by the only contribution with the coefficient $`𝒞`$. This gives

$$\mathrm{\Delta }_{nl}^{(2)}=\frac{(n2)(nl)(u+2)(3u+2)}{4(u+1)^2d^2}+O(1/d^3).$$
(36)
The general expressions (19), (IV A) are rather cumbersome, and in the subsequent sections we shall separately discuss isotropic contribution (even $`n`$, $`l=0`$) and anisotropic ones (general $`n`$, $`l0`$).

### B Isotropic sectors

Expression (19) simplifies for the most important case of the isotropic sector (even $`n`$ and $`l=0`$):

$`\mathrm{\Delta }_{n0}^{(2)}={\displaystyle \frac{n(n2)}{(d1)(d+2)^2(d+4)}}\left\{2(d+4)𝒜+6(n4)+9(d+n)𝒞\right\}.`$ (37)
This gives $`\mathrm{\Delta }_{20}^{(2)}=0`$ in agreement with the exact result $`\mathrm{\Delta }_{20}=0`$ . This means that the second-order structure function is not anomalous. The formal proof is based on certain Schwinger equation, which has the meaning of the energy conservation law; it is almost identical to the analogous proof for the Kraichnan model, given in Ref. .

For the simplest nontrivial case $`n=4`$ one obtains

$`\mathrm{\Delta }_{40}^{(2)}=8(2𝒜+9𝒞)/(d1)(d+2)^2,`$ (38)
that is, the quantity $``$ does not enter into the result. For $`n6`$, all the coefficients (IV A) contribute to the result.

In Fig. 1, we show the behavior of the quantity

$`\zeta _n\left[\mathrm{\Delta }_{n0}^{(2)}\mathrm{\Delta }_{n0}^{(2)}|_{u=\mathrm{}}\right]/n^3`$ (39)
for $`n=4`$, 6, 8 and 20 (from below to above) as a function of $`u`$ for several values of $`d`$. We have subtracted the value of the dimension for the rapid-change case, such that the curves approach zero as $`u\mathrm{}`$, and divided the difference by $`n^3`$, such that the results for different $`n`$’s have the same order of magnitude \[the quantity (19) is a third-order polynomial in $`n`$\]. It is worth noting that, since the leading coefficient (17) is independent of $`u`$, it drops from the difference $`\mathrm{\Delta }_{n0}\mathrm{\Delta }_{n0}|_{u=\mathrm{}}`$ of the exact dimensions, and in the leading order $`O(\epsilon ^2)`$ the latter is proportional to the quantity $`\zeta _n`$ introduced above.

As one can easily see from Fig. 1, the qualitative behavior of $`\zeta _n`$ depends essentially on the values of $`n`$ and $`d`$. For moderate $`n`$ and $`d`$ (e.g., $`n=4`$, 6 and 8 for $`d=2`$ and 3), finite correlation time enhances the intermittency (anomalous dimensions become more negative) in comparison with the both limits: the rapid-change ($`u=\mathrm{}`$) and quenched ($`u=0`$) ones. While the $`O(1/u)`$ correction leads to a smooth decrease of $`\zeta _n`$ for increasing $`1/u`$ from zero (in agreement with the numerical simulation of Ref.  for a shell model near the rapid-change limit), the rapid falloff of $`\zeta _n`$ is observed for increasing $`u`$ from the quenched limit $`u=0`$. As a result of the competition between these two effects, $`\zeta _n`$ is not a monotonous function of $`u`$ and has a pronounced minimum in the interval between 0 and 1.

We recall that the slopes of the functions $`\zeta _n`$ at $`u=0`$ are infinite in two dimensions for all values of $`n`$ due to the presence of poles $`(d2)`$ in the $`O(u)`$ terms in Eqs. (32) and (34); see Sec. IV A. For $`d>2`$, the slopes become finite but they still remain very steep for $`d=3`$ and lead to a rapid falloff of $`\zeta _n`$, as Fig. 1(b) shows. This fact also suggests that the quenched case, in contrast with the rapid-change one, can hardly serve as a good zero-order approximation in studying more realistic models of passive advection by the velocity field with finite correlation time.

If $`n`$ or $`d`$ is large enough, the minimum becomes less pronounced, the behavior of $`\zeta _n`$ becomes more regular (e.g., $`n=20`$ for $`d=2`$) and eventually $`\zeta _n`$ becomes a monotonically decreasing function of $`u`$ ($`n=20`$ for $`d=3`$). For such values of $`n`$ or $`d`$, the function $`\zeta _n`$ approaches zero at $`u\mathrm{}`$ from above. In other words, the $`O(1/u)`$ correction to the rapid-change limit suppresses the intermittency, in contrast with the case of moderate $`n`$.

In the limit of large $`d`$, from Eq. 36 one easily obtains

$$\mathrm{\Delta }_{n0}^{(2)}\mathrm{\Delta }_{n0}^{(2)}|_{u=\mathrm{}}=\frac{n(n2)(2u+1)}{4(u+1)^2d^2}+O(1/d^3)$$
(40)
(note that the $`O(n^3)`$ term in this approximation disappears and $`\mathrm{\Delta }_{n0}^{(2)}`$ becomes only quadratic in $`n`$). One can see that for all values of $`n`$, the difference (40) is positive, decreases monotonically when $`u`$ grows, and approaches zero from above when $`u\mathrm{}`$.

In Ref. , the $`O(1/u)`$ correction to the rapid-change case was derived by the zero-mode techniques in the limit of large $`d`$ and arbitrary (not small) values of $`\epsilon `$, for the case of a local turnover exponent ($`\epsilon =\eta `$). Although the anomalous exponents were shown to be nonuniversal (dependent on $`u`$), our results disagree with Ref.  in two respects. First, due to the universality (independence of $`u`$) of the leading term (17), the ratio (40) is of order $`O(\epsilon ^2)`$ and not $`O(\epsilon )`$. Second, the $`O(1/u)`$ correction in (40) is positive for all $`n`$, while, according to , inclusion of the finite correlation time makes the anomalous exponents more negative in comparison with the rapid-change limit for all $`n`$ and $`\epsilon `$. It is not clear whether this disagreement can be explained by some distinctions between our model (3), (4) and the velocity ensemble employed in . It is possible to show, however, that any modification of the function (4) consistent with the RG analysis performed in Ref.  and section III above, leads to an universal (independent of $`u`$) expression for the leading term in $`\mathrm{\Delta }_{nl}^{(1)}`$, so that the $`O(1/u)`$ correction to the zero-correlated limit remains of order $`O(\epsilon ^2)`$.

The changeover from the behavior typical to low spatial dimensions to the behavior described by Eq. (40) also produces interesting patterns, as illustrated by Fig. 1 for $`d=5`$ and 10.

### C Anisotropic sectors

Let us turn to the analysis of anisotropic contributions in the structure functions (12), (13), described by the dimensions (16) with $`l0`$. We recall that such contributions appear in the inertial-range expression (13) if the forcing (2) is chosen to be anisotropic, or a constant gradient of the scalar field is imposed.

An important property of the first-order result (17) is that for any fixed $`n`$, the quantity $`\mathrm{\Delta }_{nl}^{(1)}`$ increases monotonically with $`l`$ . One can say that the exponents, associated with tensor composite operators (15), exhibit a kind of hierarchy related to the degree of anisotropy: the less is the rank $`l`$, the less is the dimension and, consequently, the more important is the corresponding contribution to the inertial-range expression (13). The leading terms in the even structure functions (11) are given by the scalar operators (15) with $`l=0`$, that is, they are the same as in the model with isotropic forcing (we recall that $`n`$ and $`l`$ should be simultaneously even or odd).

This behavior is in agreement with the existing phenomenological ideas, according to which the anisotropy introduced at large scales by the forcing (boundary conditions, geometry of an obstacle etc.) dies out when the energy is transferred down to smaller scales owing to the cascade mechanism . The hierarchy of anisotropic contributions appears rather universal, being also observed for a vector (magnetic) field, advected by the Kraichnan velocity ensemble , for the scalar advected by the two-dimensional Navier–Stokes velocity field  and for the turbulent velocity field itself .

Nevertheless, the anisotropy survives in the inertial range and reveals itself in dimensionless ratios involving odd structure functions,

$$_k𝒮_{2k+1}/𝒮_2^{k+1/2}.$$
(41)
For a number of models it was shown that the skewness factor $`_1`$ decreases down the scales but slower than predicted by phenomenological theories , while the higher-order odd ratios (hyperskewness $`_2`$ etc.) increase, thus signalling of persistent small-scale anisotropy . Due to the aforementioned hierarchy of the dimensions (17), the leading terms in the odd structure functions (11) in our model are determined by the vector operators (15) with $`l=1`$, and it is easy to check the above statements from the explicit expression (17).

Of course, the $`O(\epsilon ^2)`$ contribution (19) cannot change all the above properties of the dimensions $`\mathrm{\Delta }_{nl}`$, determined by their leading $`O(\epsilon )`$ contribution (17), as far as small values of $`\epsilon `$ are concerned. However, since the dependence on $`u`$ occurs only in the $`O(\epsilon ^2)`$ contribution, it should be taken into account if one wishes to discuss how finite correlation time affects the hierarchy of the dimensions or the behavior of the dimensionless ratios.

Consider the effects of the finite correlation time on the hierarchy of the anisotropic contributions. To this aim, consider the difference of the coefficients (19) for a fixed value of $`n`$ and two neighboring subsequent values of $`l`$:

$`\mathrm{\Delta }_{n,l+2}^{(2)}\mathrm{\Delta }_{nl}^{(2)}={\displaystyle \frac{2(2l+d)𝒦_n(d,u)}{(d1)^2(d+2)^2(d+4)}},𝒦_n(d,u)\left\{2(d+4)𝒜+9(n2)\left[2(d+1)𝒞\right]\right\}.`$ (42)
(we recall that for a fixed $`n`$, all possible values of $`l`$ are either even or odd, so that the subsequent values of $`l`$ differ by 2). It is clear from Eq. (42) that the sign and the dependence on $`u`$ of the whole expression is determined by the behavior of the function $`𝒦_n(d,u)`$.

In Fig. 2, we plot the quantity $`𝒦_n(d,u)`$ as a function of $`u`$ for $`n=2`$, 4, 6 and 20 (from above to below) for $`d=2`$ (a) and $`d=3`$ (b). The function is always negative for all the cases studied and increases monotonically with $`u`$. This behavior persists in the limit of large $`d`$, as follows from the asymptotical expression (36).

We thus conclude that the $`O(\epsilon ^2)`$ contribution in the exact dimension (16) “tries to cope” with the hierarchy, set by the $`O(\epsilon )`$ term, for all values of $`n`$, $`l`$, $`d`$ and $`u`$; this effect is at its strongest for $`u=0`$ and weakens monotonically as $`u`$ increases from 0 to $`\mathrm{}`$.

Now let us turn to the dimensionless ratios $`_k`$ in (41). From the discussion below Eq. (41) and asymptotic representation (12) one obtains the power-like inertial-range asymptotic expression $`_k(mr)^{\mathrm{\Delta }_{2k+1,1}}`$ with $`\mathrm{\Delta }_{2k+1,1}`$ from Eq. (16) \[we recall that in our model $`\mathrm{\Delta }_{2,0}=0`$; see Sec. IV B\]. Due to the $`u`$-independence of the first-order answer (17), the $`O(\epsilon )`$ contribution in the exponent $`\mathrm{\Delta }_{2k+1,1}=\epsilon (d+24k^2)/2(d+2)+O(\epsilon ^2)`$ coincides with its analog for the Kraichnan model; see Refs.  for $`k=1`$ and  for general $`k`$. It completely determines the qualitative behavior of the quantities $`_k`$: for $`k=1`$ one has $`\mathrm{\Delta }_{3,1}>0`$ and the skewness factor $`_1`$ decreases with $`mr`$, while for $`k1`$ one has $`\mathrm{\Delta }_{2k+1,1}<0`$ and the higher-order ratios $`_k`$ increase for $`mr0`$.

In Fig. 3, we show the behavior of the second-order correction $`\mathrm{\Delta }_{2k+1,1}^{(2)}`$, obtained from the general formula Eq. (19) and divided by $`(2k+1)^3`$ like for the even dimensions, as a function of $`u`$ for $`k=1`$, 2, 3 and 4 (from below to above) for $`d=2`$ (a) and $`d=3`$ (b).

One can see that the effect of the $`O(\epsilon ^2)`$ correction on the inertial-range behavior of the ratios $`_k`$ is different for different $`d`$, $`k`$ and $`u`$. In two dimensions, the corrections are negative for all $`u`$ and moderate $`k`$: the decay of the skewness factor $`_1`$ for $`mr0`$ appears even slower than indicated by the $`O(\epsilon )`$ approximation, while the growth of the ratios with $`k2`$ becomes faster.

In three dimensions, the correction is negative for $`k=1`$ so that the decay $`_1`$ for $`mr0`$ is also slower than in the $`O(\epsilon )`$ approximation. For $`k=2`$, the correction is negative for small $`u`$ (so that the growth of the hyperskewness factor $`_2`$ for $`mr0`$ is faster than in the first-order approximation), but it changes its sign for some finite value of $`u`$ and the growth of $`_2`$ slows down. For $`k3`$, the corrections are negative for all $`u`$ and the growth of the corresponding higher-order ratios $`_k`$ appears slower than predicted by the $`O(\epsilon )`$ expression. One thus may conclude that for $`d=3`$, with the exception of the $`k=2`$ case, the effect of the second-order term is opposite to the tendency set by the first-order approximation.

For $`k`$ large enough and any $`d`$, the behavior of the quantities $`\mathrm{\Delta }_{2k+1,1}^{(2)}`$ becomes similar to that of the even dimensions $`\mathrm{\Delta }_{2k,0}^{(2)}`$ discussed in Sec. IV B: they decrease monotonically as $`u`$ increases, comparatively fast for small $`u`$ (due to the singularity in the slope for $`d=2`$; see Sec. IV A) and rather slow when $`u`$ becomes large enough. This follows from the fact that the $`l`$-independent contribution in the general expression $`\mathrm{\Delta }_{nl}^{(2)}`$ behaves as $`O(n^3)`$ for $`n\mathrm{}`$, while its $`l`$-dependent contribution behaves only as $`O(n)`$; see Eq. (19).

We also note that for moderate $`k`$, the quantities $`\mathrm{\Delta }_{2k+1,1}^{(2)}`$ show a nonmonotonous dependence on $`u`$ in the region of small $`u`$ and in this respect they also resemble the even dimensions; see Fig. 1 and the discussion in Sec. IV B.

## V Conclusion

We have applied the RG and OPE methods to a simple model of a passive scalar quantity advected by the synthetic Gaussian velocity field with a given self-similar covariance with finite correlation time. The structure functions of the scalar field exhibit inertial-range anomalous scaling behavior, as a consequence of the existence in the model of composite operators with negative scaling dimensions, identified with anomalous exponents.

For the special case of a local turnover exponent, the anomalous exponents are nonuniversal through the dependence on a dimensionless parameter $`u`$ that has the meaning of the ratio of the velocity correlation time and the scalar turnover time. The universality reveals itself only in the second order of the $`\epsilon `$ expansion, and we have derived the exponents to order $`O(\epsilon ^2)`$, including anisotropic contributions.

It is shown that, for isotropic contributions, the qualitative effect of finite correlation time depends essentially on $`n`$, the order of the structure function, and the space dimensionality $`d`$. For moderate $`n`$ and $`d`$, finite correlation time enhances the intermittency in comparison with the both limits: the rapid-change ($`u=\mathrm{}`$) and quenched ($`u=0`$) ones. The $`O(\epsilon ^2)`$ term shows a highly nontrivial behavior in the vicinity of the quenched limit: a rapid falloff when $`u=0`$ increases from zero, with infinite derivative at $`u=0`$ for $`d=2`$, with a pronounced minimum for $`u1`$. This irregularity shows that the time-independent advecting field can hardly be a reasonable approximation in studying more realistic models of passive advection by the velocity field with finite correlation time. The behavior near the opposite limit, $`u=\mathrm{}`$, is smooth in agreement with the existing simulation for a shell model .

The behavior changes remarkably when $`n`$ and/or $`d`$ become large enough: the correction to the limit $`u=\mathrm{}`$ due to finite correlation time is positive for all $`u`$ (that is, the anomalous scaling is suppressed in comparison with the rapid-change case), it is maximal for $`u=0`$ and monotonically decreases to zero as $`u`$ tends to infinity.

In the anisotropic sectors, the $`O(\epsilon ^2)`$ terms diminish the hierarchy revealed by the first-order terms for all values of the parameters $`n`$, $`l`$ and $`d`$; this effect is maximal at $`u=0`$ and decreases monotonically with $`1/u`$.

The effect of the $`O(\epsilon ^2)`$ corrections on the inertial-range behavior of the dimensionless ratios involving odd-order structure functions depends on $`d`$. For $`d=2`$ and moderate $`k`$ these corrections are negative; the decay of the skewness factor $`_1`$ for $`mr0`$ is slower while the growth of the higher-order ratios $`_k`$ with $`k2`$ is faster than indicated by the $`O(\epsilon )`$ approximation by Refs. . For $`d=3`$, the effect is, for most cases, opposite to the tendency set by the first-order approximation: the decay of the skewness factor slows down as well as the growth of the higher-order ratios.

Our analysis has been confined with the region of small $`\epsilon `$, where the results obtained within the $`\epsilon `$ expansion are internally consistent and undoubtedly reliable \[we recall again that, although the leading terms of the anomalous exponents are of order $`O(\epsilon )`$, the leading terms in which the effects of finite correlation time occur are of order $`O(\epsilon ^2)`$\]. We do not discuss here the serious issue of validity of the $`\epsilon `$ expansions for finite $`\epsilon =O(1)`$. One can think that, in our model, the natural region of validity of the $`\epsilon `$ expansion is restricted by the value $`\epsilon =1/2`$, where the velocity field acquires negative critical dimension (along with all its powers) and new IR singularities, related to the well-known sweeping effects, occur in the diagrams; see the discussion in Ref. . (It should be noted, however, that such singularities do not necessarily lead to a changeover in the inertial-range behavior, as shown in  for the special case of the structure function $`𝒮_2`$ for $`u=0`$). On the other hand, $`\epsilon =1/2`$ can be regarded as the upper bound of the range of validity of the model itself: the lack of Galilean covariance becomes a serious drawback of the synthetic Gaussian velocity ensemble when the sweeping effects become important. The next important step should be the analytical derivation of anomalous exponents of a passive scalar advected by the Galilean covariant velocity field; this work is now in progress.

###### Acknowledgements.

The authors are thankful M. Hnatich, A. Kupiainen, P. Muratore Ginanneschi, M. Yu. Nalimov, A. N. Vasil’ev, and A. Vulpiani for discussions. The work was supported by the Nordic Grant for Network Cooperation with the Baltic Countries and Northwest Russia No. FIN-18/2001. N.V.A. and L.Ts.A. were also supported by the program “Universities of Russia” and the GRACENAS Grant No. E00-3-24. N.V.A. and L.Ts.A. acknowledge the Department of Physical Sciences of the University of Helsinki for kind hospitality.


text:
# Heat current in parametric quantum pump

## Acknowledgments

We gratefully acknowledge support by a RGC grant from the SAR Government of Hong Kong under grant number HKU 7091/01P.


text:
# A study of correlation between the oscillation amplitude and stellar parameters for 𝛿 Scuti stars in open clusters

## 1 Introduction

The $`\delta `$ Scuti stars are variables located in the lower part of the Cepheid instability strip with spectral types from A2 to F0. These pulsators, when belonging to the Main Sequence, seem particularly suitable for determining the extent of the convective core and internal rotation rate, and thereby for probing poorly understood hydrodynamical processes occurring in deep stellar interiors. Within the last decade great efforts have been made in developing the seismology of $`\delta `$ Scuti stars (Breger, 2000; Handler, 2000). However, several aspects of the pulsating behavior of these stars (e.g. the number and the determination of the excited modes) within the instability strip are not completely understood.

Several statistical works investigated potential correlations between the oscillation amplitude and different stellar parameters. Particularly, a correlation between the amplitude, period and absolute magnitude of low amplitude $`\delta `$ Scuti stars has been established (Antonello et al., 1981). In a preliminary work for $`\delta `$ Scuti stars in the Praesepe cluster, Suárez et al. (2001) have suggested the importance of taking into account the effect of fast rotation as proposed by Michel et al. (1999) in such studies.

Here we extend this work to 5 open clusters, in order to estimate the impact of the mentioned technique in the studies of correlation. In addition, these tools allow us to obtain a rough estimation of several stellar parameters like the angle of inclination of the star or the rotation rate, which are not usually accessible. We present a study of single correlation between the parameters mentioned herebefore and the oscillation amplitude. For values indicating a reasonable degree of correlation, we also search for a possible functional dependence.

This paper is organized as follows: in Section 2 we present the observational material and the main characteristics of the selected open clusters. Details on the computation of models are given in Section 3. In Section 4, we introduce the stellar parameters corrected for the effect of fast rotation. Section 5 is dedicated to describe the technique of correlation and the analysis of the results. Finally, we bring up our conclusions in section 6.

## 2 Observational data

The oscillation amplitude (hereafter $`A_{\mathrm{osc}}`$), as classically, is taken to be the highest peak-to-peak measurement in the observed light curve. Values for $`A_{\mathrm{osc}}`$ used here are taken from Li & Michel (1999) and references therein. The interpretation of this parameter is not straightforward; it can be considered as a maximum constructive interference of the observed modes. Thus, in such an interpretation, and for a given angle of inclination, it can be taken as an indicator of the energy of these modes.With such an estimation, it is evident that short time series might induce an underestimation of the amplitude, especially in presence of beating phenomena due to close frequencies. In order to take this effect into account, we associate with each amplitude estimation a weight corresponding to the length of the time series. In Table 5 are listed such values where number 1 corresponds to effective time series shorter than 4 hours; 2, between 4 and 20 hours; 3, between 20 and 200 hours, and finally 4 corresponds to time series longer than 200 hours.

The five young open clusters considered are: Praesepe, Pleiades, Hyades, Coma Ber and $`\alpha `$-Persei, with ages between $`90`$ and $`700`$ Myr. Detailed information concerning metallicity, age, distance modulus and reddening are listed in Table 1. Apparent magnitudes and colour indexes are taken from the Rufener’s cataloge (Rufener, 1999).

For $`\alpha `$-Persei, the estimation of the global metallicity $`[M/H]`$ by spectroscopic measurements is taken from Pinsonneault et al. (1998). For the other four clusters, we used metallicities derived by M. Grenon from Geneva photometry of single stars with spectral type in the range F4-K3 (Grenon, 2000; Robichon et al., 1999). To correct the visual absorption due to interstellar matter, the relation

$`A_\mathrm{v}=3.6E(B_2V_1)`$

was used. The intrinsic colours in the Geneva photometric system were obtained from Nicolet (1981).

The distance moduli for Coma ($`4.70\pm 0.04`$ mag) and Praesepe ($`6.28\pm 0.12`$ mag) are derived from mean cluster parallaxes computed using *Hipparcos*  intermediate data (Robichon et al., 1999). In the case of Pleiades, the estimation of the distance modulus is still a subject of discussion. The *Hipparcos* values ($`5.36\pm 0.06`$) and the one obtained by Pinsonneault et al. (1998), $`5.60\pm 0.06`$ mag illustrate the range of values found in the literature (Stello & Nissen, 2001). For the Hyades, we have also access from literature, to the individual parallaxes of our target $`\delta `$ Scuti stars (de Bruijne et al., 2001) which are listed in Table 2. For the other stars we use an estimation of $`3.33\pm 0.01`$ mag by Perryman et al. (1998) from *Hipparcos* observations. Finally for $`\alpha `$-Persei, we used a distance modulus of $`6.23\pm 0.06`$ from Pinsonneault et al. (1998).

The projected velocity ($`V\mathrm{sin}i`$) values listed in Table 5 are taken from Uesugi & Fukuda (1982).

## 3 Stellar models and isochrones

The stellar evolutive models have been computed using the evolution code CESAM (Morel, 1997). We consider input physics appropriate to the mass range covered by $`\delta `$ Scuti stars . Metallicity values ($`Z`$) used to compute models, are deduced from $`[M/H]`$ assuming $`Y_{\mathrm{pr}}=0.235`$ and $`Z_{\mathrm{pr}}=0`$ as helium and heavy elements primordial concentrations, and an enrichment ratio of $`\mathrm{\Delta }Y/\mathrm{\Delta }Z=2`$. We work with OPAL radiative opacity tables (Iglesias & Rogers, 1996). For the atmosphere reconstruction, the Eddington’s $`T(\tau )`$ law is considered. Convection is described according to the classical mixing-length theory, with a solar calibrated mixing length parameter $`l_m=1.61\mathrm{H}p`$, where $`\mathrm{H}p`$ is the local pressure scale-height. For all the models, we take into account an overshooting of the mixed convective core over a distance $`\alpha _{\mathrm{ov}}=0.2\mathrm{H}p`$ following the prescription of Schaller et al. (1992).

Sets of evolutionary sequences were computed for masses between $`1M_{\mathrm{}}`$ and $`5M_{\mathrm{}}`$ from the zero-age Main Sequence to the subgiant branch. With representative sequences for each cluster, we obtain the corresponding isochrones using the Geneva isochrone program. These isochrones were then transformed into a $`M_\mathrm{v}(B_2V_1)`$ diagram using the calibration of Kunzli et al. (1997) for $`(B_2V_1)`$ and the calibration of Schmidt-Kaler (1982) for $`M_\mathrm{v}`$.

In Fig. 1 we present colour-magnitude (hereafter CM) diagrams with the observed data for each cluster. Continuous lines correspond to their associated isochrones, given the distances and ages listed in Table 1. Potential binarity and the effect of fast rotation discussed in Section 4 are expected to induce systematic shifts toward higher luminosity and lower effective temperature, compared to single non-rotating stars. Thus, in order to compensate both effects, the fit of isochrones has been made adjusting them to the bottom envelope of the cluster in the CM diagrams.

In the case of Pleiades, we use our best fit corresponding to a distance modulus of $`5.5`$ mag, which is approximately half way between the main sequence fitting result and *Hipparcos* distance (see Section 2). The corresponding age found is $`130`$ Myr.

## 4 Stellar parameters corrected for the effect of fast rotation

We will not enter here in the details of the method, already described in Pérez Hernández et al. (1999). Following this paper, we apply the method to selected $`\delta `$ Scuti stars for each of our five clusters. In Fig. 2 we illustrate such corrections for the five clusters. For a given $`V\mathrm{sin}i`$, segments represent the potential position of the star, in absence of rotation, on the CM diagram. The position of this non-rotating copartners varies with the angle of inclination, $`i`$, and the rotation rate, $`\omega `$. This last dimensionless quantity is defined by

$$\omega \frac{\mathrm{\Omega }}{\stackrel{~}{\mathrm{\Omega }}_c}\text{with}\stackrel{~}{\mathrm{\Omega }}_\mathrm{c}^2\frac{8GM}{27R_{\mathrm{p}}^{}{}_{}{}^{3}}$$

where $`\mathrm{\Omega }`$ is the angular rotational velocity of the star; $`\stackrel{~}{\mathrm{\Omega }}_\mathrm{c}`$ is the angular rotational velocity that a star with the same mass ($`M`$) and polar radius ($`R_\mathrm{p}`$) would have if the centrifugal force balanced the gravitational attraction at the equator. The non-rotating counterparts are thus, computed as a function of the angle of inclination, from $`i=90^{}`$ ($`\omega _{\mathrm{min}}`$) to $`i_{\mathrm{min}}`$ obtained for $`\omega _{\mathrm{max}}=0.95`$.

For each selected $`\delta `$ Scuti star we compute two associated segments corresponding to corrections assuming an estimation of the error in $`V\mathrm{sin}i`$ of $`\pm 10\%`$. Finally, models at the intersection of segments with the isochrone will be our non-rotating counterparts, giving a set of corrected photometric parameters, like absolute magnitudes ($`M_\mathrm{v}`$) and ($`B_2V_1`$) colour indices. In addition to this, the method gives us the possibility of obtaining an estimation of other stellar parameters like the angle of inclination of the star, the radial velocity, and the rotation rate, given in Table 6.

## 5 Searching for correlations with the oscillation amplitude

### 5.1 The errors

We search for correlations between, on one side, the $`A_{\mathrm{osc}}`$ parameter, and on the other side, stellar parameters essentially obtained after correcting for the effect of fast rotation. Here we discuss the error on these parameters.

For $`A_{\mathrm{osc}}`$, we estimate $`0.005`$ and $`0.01`$ as minimum and maximum typical errors and we calculate error bars assuming a linear relation with the weights listed in Table 5.

Considering the stellar parameters, we have to deal with different error sources: input photometric indexes, $`V\mathrm{sin}i`$ values, distances, models and calibrations (isochrones), and of course, the method to correct parameters for the effect of rotation.

The photometric input parameters, basically the $`m_V`$ and $`B_2V_1`$, have typical observational errors around $`10^3`$ magnitudes. This is between 10 to 100 times smaller than the rotation effect found for these parameters (Pérez Hernández et al., 1999; Suárez et al., 2001). The error in the estimation of $`V\mathrm{sin}i`$ (estimated to 10$`\%`$) is probably a dominant factor here. It can induce errors up to 0.15 and 0.02 mag in $`M_\mathrm{v}`$ and $`B_2V_1`$ corrected values respectively. The error on the fit of the isochrones (including those on the distance determination and on the calibration between $`B_2V_1`$ and $`T_{\mathrm{eff}}`$) can be of the same order than the one associated to $`V\mathrm{sin}i`$. However they are expected to be systematic, at least for objects within the same cluster.

The errors associated more specifically to the models and the correction for fast rotation are difficult to estimate, because we lack elements to compare with. Although for $`M_\mathrm{v}`$ and $`B_2V_1`$, they are expected to be systematic to a good extent, for the rest of estimated stellar parameter we cannot guarantee it. In particular for $`i`$, $`V`$ and $`\omega `$, we can only pretend to obtain rough estimations, which will be indicated with a ().

We thus, consider that the error on $`V\mathrm{sin}i`$ is representative of the non-systematic errors here, being the most important source of errors when searching for correlations. In addition, other observational aspects like unknown binarity, differences with the metallicity of the cluster and possible not detected differential extinction inside the cluster, can have an influence.

For corrected parameters in Table 6, the error bars correspond to values obtained using $`\pm 10\%`$ in the observed $`V\mathrm{sin}i`$. In Fig. 3 we present the different diagrams with the parameters corrected for rotation versus $`A_{\mathrm{osc}}`$, with the error bars for both.

### 5.2 The correlation technique

We computed the Spearman rank-order correlation coefficient, $`r_\mathrm{s}`$. This coefficient is preferred to the classic Pearson one (also computed as a complement) as it is more robust to outliers and does not presuppose a linear relation. For the error made in the correlation coefficient, we use the following expression from standard statistics,

$`\sigma (r_\mathrm{s})^2={\displaystyle \frac{(1r_\mathrm{s})^2}{n1}}\left(1+{\displaystyle \frac{11r_{\mathrm{s}}^{}{}_{}{}^{2}}{2n}}\right)+𝒪(n^3)`$

where n is the dimension of the sample. In addition to this, we computed, for the Spearman’s coefficient, the parameter p which gives the probability of deviation from the null-hypothesis for this coefficient.

On the other hand, for the $`r_\mathrm{s}`$ values indicating a reasonable degree of correlation, we also search for a possible functional dependence. In such cases, we have used an orthogonal distance regression routine (Press et al., 1989, see) which combines data sets with errors in both variables.

### 5.3 Results and discussion

We have applied the statistical methods described hereabove to all parameters corrected for rotation computed in Section 4. In Table 3 are listed the coefficients obtained for correlations between these parameters and the $`A_{\mathrm{osc}}`$.

We now analyze the results obtained for the absolute magnitude and the colour index. These two photometric parameters have been generally used in the literature for statistical works on $`\delta `$ Scuti stars . Particularly, in a multi-correlate analysis made by Antonello et al. (1981), a significant individual correlation between $`A_{\mathrm{osc}}`$ and $`M_\mathrm{v}`$ is found. No strong correlation is found here neither for $`M_\mathrm{v}`$ nor for $`B_2V_1`$, however we notice that the correlation coefficient is approximately two times better for corrected values of $`M_\mathrm{v}`$ than for non-corrected ones ($`r_\mathrm{s}=0.26`$). They are of the same order in the case of $`B_2V_1`$ ($`r_\mathrm{s}0.25`$). We thus conclude, following (Suárez et al., 2001), that it is worthwhile to take into account such a correction in this kind of study. We attribute the low value obtained for the coefficient in the case of the magnitude to the fact that our sample is still too limited.

A less classical parameter for such studies is the projected rotational velocity, $`V\mathrm{sin}i`$, for which an unexpected significant correlation is found ($`r_\mathrm{s}=0.76`$). This parameter has been considered before in Li & Michel (1999), however no correlation was calculated. This striking result has two possible explanations: a physical effect (linked to $`V^{}`$ or maybe rather $`\omega ^{}`$) and/or a geometric effect (associated with $`i^{}`$). Making use of the different stellar parameters we have access to, we tried to go further to understand the origin of this correlation. We can notice that there is no significant correlation neither for the $`V^{}`$ parameter nor for the $`\omega ^{}`$ parameter.

On the other hand, for $`i^{}`$, another significant correlation is found ($`r_\mathrm{s}=0.75`$) such that, the higher is the value of $`A_{\mathrm{osc}}`$ the larger is $`i^{}`$. These results suggest that the correlation observed for $`V\mathrm{sin}i`$ has a pure geometric origin. At this point, we tried to figure out if a possible bias could be induced by the method used to correct parameters for rotation. We can only figure out a possible influence of our limited sample of stars. In addition, the correlation found between $`V\mathrm{sin}i`$ and $`A_{\mathrm{osc}}`$ does not involve the correction method.

In the framework of the classical description of the oscillations, which attributes to each mode a given spherical harmonic, the preceding result would first confirm that the variation of flux observed in $`\delta `$ Scuti light curves are not dominated by the radial modes, which visibility coefficients are independent of $`i^{}`$. This is in agreement with the results obtained by Hernández et al. (1998) who showed that they could not find more than two radial modes per star, in oscillation spectra, for a group of selected stars in the Praesepe cluster.

The $`A_{\mathrm{osc}}`$ value being considered as the maximum constructive interference of the observed oscillation modes, the present results would suggest that, at least, some modes with visibility coefficients increasing with $`i^{}`$are favored in the $`\delta `$ Scuti stars oscillations compared with other modes.

As can be seen in Pesnell (1985), for instance, modes corresponding to this description are the *sectoral* modes ($`m=\pm \mathrm{}`$), in opposition to the *zonal* modes ($`(\mathrm{}=1,m=0)`$,$`(\mathrm{}=2,m=0)`$) which visibility coefficient is low for $`i0`$, and $`(\mathrm{}=2,m=\pm 1)`$ which visibility coefficient is low for $`i0`$ and for $`i90`$. In this context, our results would suggest that the $`\delta `$ Scuti oscillations are dominated, at least to some extent, by modes with $`m=\pm \mathrm{}`$. If confirmed, this interpretation would constitute an important new element in the problem of mode identification of $`\delta `$ Scuti stars .

However, we want to stress that the classical spherical harmonics system might not be the best context to interpret this result. It is well known now, that $`\delta `$ Scuti stars are fast rotators (this is confirmed here where, for our sample of $`\delta `$ Scuti stars , large $`V`$ values are found) and that the effect of rotation affects their structure to a point where the description in terms of spherical harmonics might become questionable. An alternative representation is being developed by Lignieres (Lignières, 2001) by means of a non-perturbative numerical method. Their results obtained for homogeneous ellipsoids show modes with geometric horizontal features very different from classical spherical harmonics. These results are preliminary and have to be applied to models closer to stellar internal structure, but they might be indicative of relevant geometric effects.

For projected velocity and $`i^{}`$, we have also searched for a functional dependence. In Table 4 their linear fit parameters are listed. We have to keep in mind that the coefficients of these linear fits might be sensitive to eventual systematic errors, which were neglected when searching for correlations (see Section 5). While good confidence values ($`q`$) are obtained, we can appreciate from values of $`\chi ^2`$, that scatter is important as well as errors in fit coefficients. This situation can be improved by increasing the number of target stars.

Finally, no correlation is found for the mass.

## 6 Conclusions

For a sample of 17 $`\delta `$ Scuti stars , we report here the results of a search for correlations between the oscillation amplitude and a set of estimated fundamental parameters.

In this process we confirm the interest of applying correction for fast rotation, as suggested by (Suárez et al., 2001) for this kind of study.

We found a significant correlation coefficient for the projected velocity ($`V\mathrm{sin}i`$) and show that it is due to the strong correlation between the oscillation amplitude and an estimation of the angle of inclination of the star, $`i^{}`$. We observe that the oscillation amplitudes increase with $`i^{}`$. This striking result is discussed in terms of the visibility of modes. In the context of spherical harmonics, this suggests that the oscillations of $`\delta `$ Scuti stars would follow a selection rule in favor of *sectoral* modes (i.e. $`m=\pm \mathrm{}`$ modes).

An alternative interpretation might turn up with future developments of the work initiated by Lignieres (Lignières, 2001).

A future increase of our sample of $`\delta `$ Scuti stars in clusters would help to refine the present results.


text:
# Mesoscopic effects in superconductor-ferromagnet-superconductor junctions

## Abstract

We show that at zero temperature the supercurrent through the superconductor - ferromagnetic metal - superconductor junctions does not decay exponentially with the thickness $`L`$ of the junction. At large $`L`$ it has a random sample-specific sign which can change with a change in temperature. In the case of mesoscopic junctions the phase of the order parameter in the ground state is a random sample-specific quantity. In the case of junctions of large area the ground state phase difference is $`\pm \pi /2`$.

This work has been motivated by recent experiments  on the superconductor-metallic ferromagnet-superconductor junctions where the ground state of the system with the superconducting phase difference equal to $`\pi `$ has been observed. It is well known that the sign of the critical supercurrent of pure SFS junctions oscillates with the width $`L`$ of the ferromagnetic region . This is due to the difference in Fermi momentum of electrons with different spins at the Fermi energy causing the superconducting wave function in the ferromagnetic region to oscillate with the characteristic distance $`\mathrm{}/(k_{}k_{})`$. Here $`k_{},k_{}`$ are Fermi momenta of up and down spins, which are different because of the finite exchange spin-splitting energy $`I`$ in ferromagnets.

It is important to note that in pure junctions the characteristic distance of oscillations is inversely proportional to $`I`$ and that at zero temperature the modulus of the critical current does not decay exponentially with $`L`$.

In the opposite limit of disordered ferromagnets $`LL_I=\sqrt{\frac{D}{I}}l`$ the average critical current $`I_c`$ decays exponentially with $`L`$ and with temperature $`T`$

$$I_cRe\mathrm{exp}(\frac{L}{\xi _F});\xi _F=\sqrt{\frac{\mathrm{}D}{2(\pi T+iI)}}$$
(1)
(where brackets $``$ stand for averaging over the random realizations of the scattering potential and $`D`$ is the the electron diffusion coefficient in the ferromagnet. It is interesting however, that in addition to the exponential decay, the average critical current oscillates as a function of $`L`$. It also changes sign as a function of $`T`$, provided $`D/L^2\stackrel{}{<}I`$.

The value of $`I_c`$ can be negative which means that in the ground state of the junction the phase difference $`\varphi `$ of the superconducting order parameter is $`\pi `$ rather than zero. The experimentally measured critical current in this case is the absolute value $`|I_c|`$.

The oscillations of the average critical current have been observed experimentally . We would like to note, however, that according to Eq.1 , these oscillations can be observed only in the case when the exchange spin-splitting energy in the ferromagnet is relatively small so that the characteristic distance of oscillations is large. This limits significantly the choice of ferromagnetic metals which can be used in the junctions to observe this effect.

In this paper we show that the exponential decay of the supercurrent with $`L`$ in Eq.1 originates from the averaging procedure. Before averaging over the impurity configurations at $`LL_I`$ the supercurrent $`I_s(\varphi ,L)`$ has a random sample-specific sign while its modulus does not decay exponentially. The practical consequence of this conclusion is that the Josephson effect survives in the case of ferromagnets with large $`I`$ when $`LL_I`$. We would like to mention that this feature is a particular case of a general statement that the Friedel oscillations in disordered metals do not decay exponentially at zero temperature .

Below we discuss the limit of thick SFS junctions $`LL_I`$ when their superconducting properties are determined by mesoscopic effects. We show that in this case the critical current of the junctions does not depend on the spin splitting $`I`$. The ground state of such a junction is, generally speaking, doubly degenerate with the phase difference having a random sample-specific modulus distributed in the interval ($`0,\pi `$). We also show that the critical current undergoes random oscillations as a function of $`T`$. Thus the junctions with $`LL_I`$ should exhibit the same sequence of effects as the junctions with $`L<L_I`$.

The energy of the Josephson junction $`E(\varphi )`$ is an even and periodic function of the phase of the order parameter. It can be represented in the form

$$E(\varphi )=\underset{n=1}{\overset{\mathrm{}}{}}E_n\mathrm{cos}n\varphi $$
(2)
while the current through the contact is determined by the relation $`J=2e\frac{dE}{d\varphi }`$. The coefficients $`E_n`$ are random sample-specific quantities. For $`LL_I`$ all average coefficients $`E_n\mathrm{exp}(L/L_I)`$ are exponentially small. In this case the typical values of the coefficients $`E_n`$ can be estimated from their variances. To simplify the analysis we consider the case when the superconductors and the ferromagnet are separated by tunneling barriers of small transparency. We will show below that at $`D/L^2T`$ and $`L\xi `$ we have

$$(E_1)^2=S\xi ^2(\frac{g}{8\pi \nu _0D})^4(\frac{D}{2\pi ^2L^2})^2$$
(3)
$$(E_2)^2=S\xi ^2(\frac{g}{8\pi \nu _0D})^8\frac{1}{(4\pi )^3}\frac{\xi ^2D^2}{L^2}$$
(4)
(where $`S`$ is the area of the junction, $`\nu _0`$ is the density of states per spin in the ferromagnet, $`\xi `$ is the zero-temperature coherence length in the superconductor and $`g`$ is the conductance per unit area of the surface between the ferromagnet and the superconductor. Eqs. 3,4 correspond to the diagrams shown in Figures 1.b),c). In the case when the tunneling transmission coefficient of the insulator between the superconductor and the ferromagnet is small, the typical amplitude $`((E_2)^2)^{1/2}`$ of the second harmonic contains an additional power of $`g`$ compared to the amplitude $`((E_1)^2)^{1/2}`$ of the first harmonic and therefore for typical samples of low $`g`$, $`E(\varphi )`$ is well approximated by the first harmonic. Since $`((E_1)^2)^{1/2}`$ does not decay exponentially with $`L`$, it is much larger than $`E_1`$. This means that $`E_1`$ has a random sign. In the cases when $`E_1<0`$ the ground state of the junction corresponds to $`\varphi _{GS}=\pi `$.

The correlation function $`E_1E_2`$ contains an additional power of $`g`$ with respect to $`(E_2)^2`$ and can be neglected which means that $`E_1`$ and $`E_2`$ are random uncorrelated quantities. In rare samples where the amplitude of the first harmonic is small $`(E_1E_2)`$ one has to take into account the second harmonic in Eq.2. In this case the ground state is doubly degenerate and corresponds to the phase difference

$$\varphi _{GS}\pm \mathrm{arccos}\frac{E_1}{4E_2}$$
(5)
Therefore the sample-specific random absolute value of the ground state phase difference is in the interval $`0<|\varphi _{GS}|<\pi /2`$, provided $`|E_1|<4|E_2|`$. The probability of such an event is of order

$$\sqrt{\frac{(E_2)^2}{(E_1)^2}}\xi L(\frac{g}{8\pi \nu _0D})^2$$
(6)
We would like to mention that in the case of a transparent superconductor-ferromagnet boundary $`\varphi _{GS}`$ is a random sample-specific quantity of order one.

Let us now discuss the temperature dependence of the critical current $`I_c(T)`$. To do so we calculate the correlation function

$$E_1(T_1)E_1(T_2)\mathrm{exp}(2L\sqrt{\pi (T_1+T_2)/D})$$
(7)
It follows from Eq.7 that the quantities $`(E_1(T))^2\mathrm{exp}(2^{3/2}L\sqrt{T/D})`$ and $`(E_1(T)E_1(0)\mathrm{exp}(2L\sqrt{T/D})`$ decay with different rates as $`T`$ increases. This indicates that in addition to an exponential decay, the quantity $`E_1(T)`$ exhibits random sample-specific oscillations with a period of order $`T^{}D/L^2`$.

The results presented above were obtained in the approximation when the variations of the phase of the order parameter along the superconductor-ferromagnet surface are neglected. This is a good approximation for the samples of small area. Below we show that in the samples of large area the possibility of spatial fluctuations of the order parameter phase along the superconductor- ferromagnet surface leads to an average critical current which is proportional to the area of the junction and does not decay exponentially even in the limit when $`LL_I`$. In this case the ground state of the junction is doubly degenerate and

$$\varphi _{GS}=\pm \pi /2$$
(8)
This can be shown using an expression for the effective Josephson energy

$$E=E_c(𝝆)\mathrm{cos}\varphi (𝝆)𝑑𝝆+\frac{N_s}{2m}𝑑𝐫(\varphi (𝐫)^2)$$
(9)
where $`𝝆`$ is the coordinate along the surface between the superconductor and the ferromagnet, the $`𝐫`$-integration is performed in the bulk of the superconductors, $`m`$ is the electron mass, and $`N_s`$ is the superfluid density in superconductors. Eq.9 is valid on a scale larger than $`L`$ along the surface. The random function $`E_c(𝐲)`$ is characterized by its average $`E_c(𝝆)=0`$ and a quickly decaying (at $`|𝝆𝝆^{\mathbf{}}|>L`$) correlation function $`E_c(𝝆)E_c(𝝆^{\mathbf{}})`$ $`E_c(𝝆\mathbf{}𝝆^{\mathbf{}})E_c(0)(E_1(S=L^2))^2/L^4`$. The phase difference $`\varphi (𝝆)=\varphi +\delta \varphi (𝝆)`$ is a random function of $`𝝆`$. At small $`E_1`$ $`\delta \varphi (𝝆)1`$. Minimizing Eq.9 with respect to $`\delta \varphi (𝝆)`$ we get an expression for the effective energy per unit area

$$E_{eff}=E_1^{eff}\mathrm{sin}^2\varphi $$
(10)
(where $`E_1^{eff}\frac{m}{N_s\xi L^2}(E_1(S=L^2))^2`$), determining the average critical current density as $`I_c=2eE_1^{eff}`$ and the average phase difference in the ground state as given by Eq.8.

The results in Eqs.3,4 were calculated microscopically describing the SFS junction by a Hamiltonian of the form

$$\widehat{H}=\widehat{H}_{BCS}+\widehat{H}_T+\widehat{H}_F$$
(11)
where $`\widehat{H}_{BCS}`$ is the Hamiltonian of superconducting leads, the Hamiltonian

$$\widehat{H}_T=t\underset{i=1,2;\alpha }{}_{S_i}d𝝆_i(\mathrm{\Psi }_i^+(𝝆_i,z_i;\alpha )\mathrm{\Psi }_F(𝝆_i,z_i;\alpha )+h.c)$$
(12)
describes tunneling between superconductors labeled by $`i=1,2`$ and the ferromagnetic metal ($`z_1=0,z_2=L`$), labeled by index $`F`$, and $`\alpha `$ is the spin index. The integration is taken over the surfaces between superconductors and the ferromagnetic metal. The last term in Eq.11 corresponds to the disordered ferromagnetic metal unperturbed by the presence of superconductors, where spin up and down bands are split by the exchange field $`I`$

$$\widehat{H}_F=\widehat{H}_0+I𝑑𝐫\mathrm{\Psi }_F^+(𝐫;\alpha )\sigma _{\alpha \beta }^z\mathrm{\Psi }_F(𝐫;\beta ).$$
(13)
Here $`\widehat{H}_0`$ is the Hamiltonian of noninteracting electrons which contains the operators of the kinetic energy and a random field $`U\left(𝐫\right)`$ . We assume that the random potential is white-noise correlated so $`U\left(𝐫\right)=0`$ and $`U\left(𝐫\right)U\left(𝐫^{}\right)=\left(2\pi \nu _0\tau \right)^1\delta \left(𝐫𝐫^{}\right)`$ (where $`\tau `$ is the mean free time).

In the lowest order in tunneling through the superconductor-ferromagnet boundary we get

$$E_1=\frac{T}{2}\underset{ϵ_k;\alpha }{}t^4_{S_1;S_2}d𝝆_1d𝝆_2d𝝆_1^{}d𝝆_2^{}\{F^+(ϵ_k;𝝆_1,0;𝝆_1^{},0)G_{\alpha ;\alpha }(ϵ_k;𝝆_1^{},0;𝝆_2^{},L)F(ϵ_k;𝝆_2^{},L;𝝆_2,L)G_{\alpha ;\alpha }(ϵ_k;𝝆_2,L;𝝆_1,0)+h.c\}$$
(14)
where $`ϵ_k=\pi (2k+1)T`$ is the Matzubara frequency and $`k=0,1,\mathrm{}`$ is an integer. We use the usual diagrammatic technique for averaging the products of electron Green’s functions . Diagrams for correlation functions $`\left(E_1\right)^2`$ and $`\left(E_2\right)^2`$ are shown in Fig.1.b),c). There are two important blocks in these diagrams. The first one corresponds to diffuson and cooperon ladder made of single particle Green’s functions in the ferromagnet. In the case of a large $`I`$ only diffusons and cooperons for parallel spins survive. They are equal to each other and the same as in the absence of the external magnetic field, satisfying the equation

$$\left(D^2+\left|ϵ_kϵ_k^{}\right|\right)C(\left|ϵ_kϵ_k^{}\right|;𝐫,𝐫^{})=\theta \left(ϵ_kϵ_k^{}\right)\delta \left(𝐫𝐫^{}\right)$$
(15)
We would also like to mention that for $`I\tau 1`$ we can neglect contributions of diagrams shown in Fig.1.d).

The second block is the ladder made of the anomalous Green’s functions $`F(ϵ_k;𝐫,𝐫^{})`$ in the superconductor. This average can be expressed in terms of the averaged product of the advanced and retarded Green’s function in the normal metal. Since

$$F(ϵ_k;𝐫,𝐫^{})=i\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{dE}{2\pi }\frac{\mathrm{\Delta }}{\mathrm{\Omega }^2+E^2}\left(G^R(E,𝐫,𝐫^{})G^A(E,𝐫,𝐫^{})\right)$$
(16)
(where $`\mathrm{\Omega }^2\mathrm{\Delta }^2+ϵ_k^2`$), the averaging $`F(ϵ_k;𝐫,𝐫^{})F(ϵ_k^{};𝐫,𝐫^{})`$ is equivalent to the averaging of Green’s functions in the normal metal. The diffusion propagator $`F(ϵ_k;𝐫,𝐫^{})F(ϵ_k^{};𝐫,𝐫^{})`$ $`P(\mathrm{\Omega },\mathrm{\Omega }^{};𝐫,𝐫^{})`$ obeys the diffusion equation

$$\left(D^2+\mathrm{\Omega }+\mathrm{\Omega }^{}\right)P(\mathrm{\Omega },\mathrm{\Omega }^{};𝐫,𝐫^{})=\frac{1}{2}\frac{\mathrm{\Delta }}{\mathrm{\Omega }}\frac{\mathrm{\Delta }}{\mathrm{\Omega }^{}}\delta \left(𝐫𝐫^{}\right)$$
(17)
The result of integration of four electron Green’s functions over the surface in diagrams shown in Fig.1.b),c) is estimated as $`\tau ^2g`$ , where $`g`$ is the dimensional tunneling conductance per unit area of the boundary.

The solution of Eq.15 satisfying the boundary conditions $`{\displaystyle \frac{d}{dz}}C(\left|ϵ_kϵ_k^{}\right|;𝐫,𝐫^{})=0`$ at the superconductor - ferromagnet metal surfaces is given by ( for $`z>z^{}`$ )

$$C(\left|ϵ_kϵ_k^{}\right|;𝐫,𝐫^{})=\theta \left(ϵ_kϵ_k^{}\right)\frac{dq_xdq_y}{(2\pi )^2}\frac{\mathrm{cosh}(Q(zL/2))\mathrm{cosh}(Q(z^{}+L/2))}{DQ\mathrm{sinh}QL}e^{iq_x(xx^{})+iq_y(yy^{})}$$
(18)
where $`Q=\sqrt{q^2+{\displaystyle \frac{\left|ϵ_kϵ_k^{}\right|}{D}}}`$ . Substituting this expression into diagrams shown in Fig.1.b)c) we get

$$\left(E_1\right)^2S\xi ^2(\frac{g}{8\pi \nu _0D})^4T^2\underset{ϵ_k;ϵ_k^{}}{}\theta \left(ϵ_kϵ_k^{}\right)\frac{d^2𝐪}{(2\pi )^2}\frac{1}{\left(Q\mathrm{sinh}QL\right)^2}$$
(19)
$$\left(E_2\right)^2S\xi ^4D^2(\frac{g}{8\pi \nu _0D})^8T^2\underset{ϵ_k;ϵ_k^{}}{}\theta \left(ϵ_kϵ_k^{}\right)\frac{d^2𝐪}{\left(2\pi \right)^2}\frac{1}{\left(Q\mathrm{sinh}QL\right)^4}$$
(20)
Calculating integrals in Eqs.19,20 we arrive at Eqs.3-4.

To calculate the temperature dependence $`E_1\left(T_1\right)E_1\left(T_2\right)`$ we should substitute $`ϵ_k=\pi \left(2k+1\right)T_1`$ and $`ϵ_k^{}=\pi \left(2k+1\right)T_2`$ into Eq.19 . Then for $`L_{T_1}=(D/T_1)^{1/2},L_{T_2}=(D/T_1)^{1/2}<L`$ we get $`\left(\mathrm{sinh}QL\right)^2=4\mathrm{exp}[2L(q^2+(\pi T_1+\pi T_2)/D)^{1/2}]`$ and finally we arrive to Eq.7.

In conclusion we have shown that the critical current of a mesoscopic SFS junction at small temperatures does not decay exponentially with the ferromagnet thickness. It has a random sign, which changes with temperature. The ground state phase difference of the junction is a random quantity $`0<\varphi _{GS}<\pi /2`$. In the case of junctions of large area the phase difference is $`\varphi _{GS}=\pi /2`$. Let us estimate the typical value of the critical current using Eq.3. Taking, for example, iron as a ferromagnet with $`L\xi =3\times 10^6`$ cm, the area of the surface $`S10^7`$ cm2, the elastic mean free path $`l10^6`$ cm and assuming that the transmission coefficient through the superconductor-ferromagnet boundary is of order one we get $`I_c10^6`$ A. The estimate based on Eq.3 scales with the junction area as $`\sqrt{S}`$ and is valid at small $`S`$. For junctions of large area, Eq.10 implies that the critical current is proportional to $`S`$.

We express our thanks to M. Gershenzon, A.D.Kent and Z. Radović for valuable discussions.

This work was supported by Division of Material Sciences, U.S.National Science Foundation under Contract No. DMR-9205144 and (ZA) by Russian Fund for Fundamental Research 01-02-17794.