text:
# Remote sensing of bubble clouds in seawater

## 1 Introduction

“The effect of bubbles on the color of the sea may be observed in breaking waves… Where a great many bubbles have been entrained by a breaking wave it is white. But where there are fewer of them it is blue-green or green, brighter than the sea but not as bright as the foamiest parts of the wave. Even after a wave has broken and the water is again quiescent, a pastel green patch often remains, slowly fading into the surrounding sea as the bubbles dissipate. Thus the effect of bubbles on the color of the sea is similar to that of solid particles (Bohren, 1987).” Bubbles within the water and foam on its surface (Bukata, 1995; Frouin et al., 1996; Stramski, 1994) can predominate in determining the radiative transfer properties of the sea surface at higher wind speeds. However, there is a limited knowledge about the radiative transfer properties of bubble clouds, their inherent optical properties (IOP), and their global climatology. Mobley (1994) and Bukata (1995) discuss qualitatively the surface properties of bubble clouds. Frouin et al. (1996) performed spectral reflectance measurements of sea foam at the Scripps Institution of Oceanography Pier. They observed a decrease of the foam reflectance in the near-infrared and proposed that the foam reflectance can not be decoupled from the refelectance by bubbles. Stramski (1994) concentrates on light scattering by submerged bubbles in quiescent seas and shows the scattering coefficient and the backscattering coefficient at 550 nm in comparison with scattering and backscattering coefficients of sea water as estimated from the chlorophyll-based bio-optical models for Case 1 waters. In this exploratory paper, we report on the influence of bubble clouds generated by breaking waves on the remote sensing reflectance and calculate not only the inherent optical, but also apparent optical properties using the radiative transfer model. We show that the optical effects of bubbles on remote sensing of the ocean color are significant. Furthermore, we present a global map of volume fraction of air in water. This map, together with the parameterization of the microphysical properties, shows the significance of bubble clouds on the global albedo of incoming solar energy. By proxy, we show the influence of the bubble clouds on the remote sensing retrieval of organic and inorganic components of the natural waters. It is worth mentioning that the bubble clouds coincide with the upper range of the euphotic zone and will, therefore, contribute to the dynamics of the upper-ocean boundary layer, heat distribution, and sea surface temperature (Thorpe et al., 1992). In fact, our initial motivation for this work was an observation that the asymptotic radiance distribution is established close to the ocean surface in apparent contradiction with theoretical studies (Flatau et al., 1999). Thus, the light field must become diffuse at shallower depths than usually modeled. This leads to search for alternative mechanisms influencing the light distribution. Thus, the importance to light scattering of the bubble clouds goes beyond the remote sensing issues considered in this work.

In the next section, we discuss in more detail the microphysical and morphological properties of bubble clouds, because they have a direct bearing on their optical properties and radiative transfer.

## 2 Physical properties of bubble clouds

### 2.1 Morphology of bubble clouds in natural waters

Individual bubble clouds are generated by breaking waves, persist for several minutes (Thorpe, 1995), and reach to mean depths of about $`4H_s`$, where $`H_s`$ is the significant wave height, but with some clouds extending to about $`6H_s`$. There is evidence that at high wind speeds, separate bubble clouds near the surface coalesce, producing a stratus layer (Thorpe, 1995). Fig 1 is based on a sonograph of Thorpe (1984). The “bubble-stratocumulus” (b-Sc) is often observed by acoustic means (Farmer and Lemon, 1984; Thorpe, 1995). The depth of the b-Sc layer is related to the wind speed and wind variability, but more specifically it is set by larger waves, such us those breaking predominantly in groups (Thorpe, 1995). The “stratus layer” description should not be taken too literally. For sufficiently high winds there will be significant concentrations through out the upper layer, but the variability within this layer can be very high.

### 2.2  Bubble cloud climatology

According to Thorpe et al. (1992) in the absence of precipitation and in wind speeds exceeding about 3 m s<sup>-1</sup>, wave breaking generally provides the dominant source of bubbles. The wind speed $`W_{10}`$ at 10 m above the mean sea surface level is used to parameterize the volume fraction of air in water, $`f=V_{\mathrm{air}}/V_w`$, where $`V_{\mathrm{air}}`$ is the volume of air, and $`V_w`$ is the volume of water. This parameterization is an approximation of a more complex wind-wave relationship.

Figure 2 shows the volume fraction of air in water estimated from the $`W_{10}`$ winds for January of 1992. The volume fraction and the wind speed are assumed to follow a non-linear relationship (Walsh and Mulhearn, 1987)

$$f=f_0W_{10}^{4.4}+f_1$$
(1)
The coefficients $`f_0`$ and $`f_1`$ were calculated at 2m assuming that f= $`10^8`$ for $`W_{10}=6.2\mathrm{m}\mathrm{s}^1`$ and f= $`10^7`$ for $`W_{10}=10.5\mathrm{m}\mathrm{s}^1`$ (Walsh and Mulhearn, 1987). Vagle and Farmer (1992) show that the volume fraction decreases with depth, changing from about $`10^6`$ at 0.3 m to $`10^7`$ at 2.7 m. We base our parameterization on these findings and extrapolate (1) to near-sea-surface depth using exponential fit.

The monthly averages of $`f`$ were obtained by employing 1992 daily surface winds from the NCAR/NCEP (Kalnay et al., 1996) reanalysis project, and averaging the daily volume fractions for each month. Thus, Figure 2 is based on the variability of the wind field on the scale of one day. In the winter of the northern hemisphere, one can observe maxima associated with the midlatitude storm tracks in the Northern Pacific. Cyclogenesis, common in western parts of the oceans during the winter, contributes to mixing and a large bubble cloud volume fraction. This can be observed to the east of the North American continent. The Intertropical Convergence Zone (ITCZ) region, with its associated deep convection, may also be a region of enhanced production of bubbles. The winds of the Southern Ocean have a strong effect on bubble formation during both summer and winter. In subtropical regions, to the west of the continents, the subsidence associated with the descending branch of the Hadley circulation is responsible for the relative minimum in $`f`$. It should be stressed that these results are qualitative and that they can be improved by more detailed breaking wave climatology models (Kraus and Businger, 1994). The data in Fig 2 is indicative of regions where bubble clouds are potentially important in the interpretation of remotely sensed reflectance.

### 2.3 Optical thickness

The size distribution of bubble clouds determines the optical thickness and is, therefore, one of the most critical parameters entering the theory. There are two assumptions which simplify the development here: (a) we consider light scattering in the geometric optics regime for which the ratio proportional to bubble radius to wavelength $`x=2\pi r/\lambda `$ is large and (b) we assume that there is an effective radius $`r_{\mathrm{eff}}`$ which determines the optical properties of the size distribution. Both (a) and (b) are quite probable. The size parameter, $`x=50`$, corresponding to a bubble radius of approximately 5 $`\mu \mathrm{m}`$, is already in the geometric optics regime, and $`r_{\mathrm{eff}}=10\mu \mathrm{m}`$ will satisfy (a). For bubbles with an effective radius of $`r_{\mathrm{eff}}`$, the volume attenuation (equal to scattering for a non-absorbing sphere) can be expressed as

$$b=Q_{\mathrm{sca}}\frac{N}{V_w}s=2\frac{V_{\mathrm{air}}}{V_w}\frac{N}{V_{\mathrm{air}}}s=2f\frac{s}{v}.$$
(2)
Thus

$$b=2\frac{f}{r_{\mathrm{eff}}}$$
(3)
where $`s=\pi r^2`$ is the cross-section of a bubble with radius $`r`$, $`v`$ is the volume of such a bubble. $`N`$ is the number of bubbles in volume $`V_w`$ of water and $`f`$ is the fraction of air in a volume of water. The effective radius is defined as $`r_{\mathrm{eff}}=v/s`$ (Stephens et al., 1990; King et al., 1993; Bricaud and Morel, 1986). Scattering efficiency defines how much of incoming light is being “blocked” by a particle by scattering processes, for large size parameters $`Q_{\mathrm{sca}}`$ tends to 2. This issue is discussed in detail by Bohren and Huffman (1983).

The physical significance of $`f`$ comes from the fact that it is determined by the large scale forcing such as the wind field. Thus, for a given synoptic or climatological setting, the mixing ratio $`f`$ is, to some extent, pre-determined. On the other hand, the effective radius $`r_{\mathrm{eff}}`$ depends on processes of much smaller scale then the large scale. These processes are coagulation, coalescence, coating by organic material, saturation, buoyancy, pressure, etc. Thus, Eq. 3 defines the optical properties of bubble clouds on both the large- and sub-scales. The expression $`b=2f/r_{\mathrm{eff}}`$ holds for polydispersions, and the only difference with a monodispersion is that $`r_{\mathrm{eff}}`$ is defined via distribution averaged $`s`$ and $`v`$.

### 2.4 Effective radius of the size distribution

Numerous observations of bubble size distributions are reported in the literature based on acoustic, photographic, optical, and holographic methods (Wu, 1988a). Akulichev and Bulanov (1987) summarize results from 22 experiments using different techniques. As the origins of bubbles are biological (within the volume and at the bottom), as well as physical (at the rough surface), we may expect large regional and temporal variations of bubble concentration between coastal and open oceanic waters and between plankton bloom or no-bloom conditions (Thorpe et al., 1992).

Figure 3 presents a comparison of bubble spectra under breaking waves and quiescent sea. Recently published observations, using laser holography near the ocean surface, have shown that the densities of 10 to 15$`\mu `$m radius bubbles can be as high as $`10^6`$ (per cubic meter per micron radius increment) within 3 m of the surface of quiescent seas (O’Hern et al., 1988). These results are plotted as solid squares connected with a vertical solid line. The majority of bubbles injected into the surface layers of natural waters is unstable, either dissolving due to enhanced surface tension and hydrostatic pressures or rising to the air-water interface where the bubbles break (Johnson and Wangersky, 1987). However, bubbles with long residence times, i.e. stable microbubbles, have been observed. For example Medwin (1977) observed nearly $`2.5\times 10^6`$of bubbles per cubic meter in the radius range $`18355`$ $`\mu `$m for small wind speeds. One of the stabilization mechanisms (Mulhearn, 1981; Johnson and Wangersky, 1987) assumes that the surfactant material is a natural degradation product of chlorophyll, present in almost all photosynthesizing algae. Isao et al. (1990) have observed very large populations of neutrally buoyant particles with radii between $`0.11\mu `$m. Johnson and Wangersky (1987) and Thorpe et al. (1992) proposed another stabilization mechanism based on monolayers of adsorbed particles. Numerical modeling (Thorpe et al., 1992) can be used to study the effects of water temperature, dissolved gas saturation levels, and particulate concentrations on the size distribution of subsurface bubbles. The results of such numerical models provide additional evidence for the existence of a small size bubble fraction which is not adequately measured by acoustic or photographic techniques. The dashed line on Fig. 3 presents the mean concentration in the model steady-state for the water temperature of 0C. It can be seen that maximum concentration is around $`15\mu `$m and it is 2.5 orders of magnitude larger than for 100$`\mu `$m bubbles. Other results presented (Figure 3) are those of Johnson and Cooke (1979) observations at 4m in wind speed 11-13$`\mathrm{m}\mathrm{s}^1`$ (open squares). In situ acoustic measurements of microbubbles at sea by Medwin (1977) are plotted as solid triangles. The solid triangles joined by a solid line are concentrations at 4m depth and 3.3m$`s^1`$ wind speed. These spectra were obtained on August 7, 1975 in Monterey Bay. The solid triangles are for midafternoon, February 10-16, 1965 at Mission Bay, San Diego, 3m below the surface, and in 1.7-2.8$`\mathrm{m}\mathrm{s}^1`$ winds. The open triangles are from Baldy (1988) and include data at 30cm depth with wind and swell and at 25cm with wind only. They are based on extensive laboratory experiments. The solid hexagons are data from Medwin and Breitz (1989) acoustic measurements obtained in the open sea at 25cm depth under the water surface during 12$`\mathrm{m}\mathrm{s}^1`$ winds under spilling breakers. The literature reviewed here and encapsulated in Figure 3 is a mixture of descriptions of the effect of active wave breaking, and of stabilised microbubbles observed largely in coastal situations. Currently, it is not clear how to parameterise the stabilized bubbles.

From results such as those presented in Fig. 3 in the case of “transient,” open ocean bubbles, it can be estimated that the size distribution follows the power law dependence $`n(r)r^a`$ and $`a4`$ (Walsh and Mulhearn, 1987; Wu, 1988a). Even though small microbubbles may not contribute to the total mass, they may be important for the light scattering. Therefore, it is of interest to estimate the contribution of small bubbles to the optical thickness. Assuming that microbubbles are spherical ($`V_{\mathrm{air}}=4/3\pi r^3N`$) we can show that the optical thickness ($`\tau =bh`$) of a layer with geometrical thickness $`h`$ is

$$bhV_{\mathrm{air}}^{2/3}N^{1/3}$$
(4)
Both $`h`$ and $`V_{\mathrm{air}}`$ are assumed to be fixed. It can be seen from Eq. 4 that the contribution to optical thickness by very small particles will be the same as that by very large particles if their concentration varies as $`N^3`$ (or steeper). Indeed is the case (Walsh and Mulhearn, 1987; Wu, 1988a).

What remains to be defined is the effective radius $`r_{\mathrm{eff}}`$. On the basis of measurements, estimates of small particle fraction, existence of background 1$`\mu `$m microbubbles, modeling predictions, and steep slope of microbubble size distributions, we decided to use $`r_{\mathrm{eff}}=10\mu \mathrm{m}`$ as a typical “radiative response radius.” This choice does not exclude existence of larger or smaller particles. The real value may be between $`1\mu `$m for stabilized particles and $`50\mu `$m for the open ocean and will depend on many environmental factors such as storm passage, wind speed, swell, wind variability, phytoplankton concentration, water temperature, gas saturation, and other properties. The 10-15 fold increase in the size of effective radius or similar decrease in the air volume fraction will reduce the importance of air bubbles to a very small effect. It may be instructive to calculate the scattering coefficient $`b`$ for typical size distribution of bubbles in water. In such case we have

$$b=Q_{\mathrm{sca}}\pi r^2\frac{dN(r)}{V_w}$$
(5)
or $`b=2fs/v`$ where $`s=\pi r^2𝑑N(r)`$ and $`v=4/3\pi r^3𝑑N(r)`$, and $`dN(r)`$ is the number of bubbles between $`r`$ and $`r+dr`$ in a volume of water $`V_w`$. For typical size distribution of bubbles in water $`dN(r)/dr1/r^4`$ we have

$$r_{\mathrm{eff}}=\frac{4}{3}\mathrm{ln}(r_1/r_0)/\left(\frac{1}{r_0}\frac{1}{r_1}\right)$$
(6)
Consider $`r_1=150`$ and $`r_0=10`$ micrometers. This gives $`r_{\mathrm{eff}}3.6r_0`$ which shows that the choice of small bubble cut-off is important for the bubbles’ optical properties. However, the choice of this cutoff is non-trivial because the spectrum of small bubbles is not understood well at present.

We close this section with some general comment about the effective radius. It is not a directly measurable quantity and, in essence, it defines how dispersed the given amount of mass is. Scattering of incoming solar radiation is sensitive to total projected surface rather than to total mass. For this reason the effective radius is commonly used in radiation calculations. However, it should be stressed that the effective radius is a semi-inherent optical property because it carries information not only about the size itself, but also about the orientation of particles, their morphology, coating, size distribution, or departure from a spherical shape. In addition, estimates of effective radius, as used in satellite remote sensing, often contain bias due to unrealistic assumptions about other optical properties such as optical thickness, leakage of photons due to horizontal transfer, wavelengths, or technique employed in retrieval. In that sense the effective radius is also used (or abused) as a semi-apparent optical property.

## 3 Results

### 3.1 Numerical model

The numerical radiative transfer model used in this study is a slightly modified version of the Hydrolight 3.0 code (Mobley, 1994; Mobley et al., 1994). In brief, this model computes from first principles the radiance distribution within, and leaving, any plane-parallel water body. Input to the model consists of the absorbing and scattering properties of the water body, the nature of the wind-blown sea surface and of the bottom of the water column, and the sun and sky radiance incident on the sea surface.

Pure sea water absorption and scattering coefficients are determined from the data of Pope and Fry (1997). 35 model wavebands were specified to cover the 400-700 nm region with a typical resolution of 10nm. The water column was specified as infinitely deep. Up to 62 depth layers, extending to 50 meters, were specified with a resolution of 0.5 m close to the surface. A clear sky was assumed, but the diffuse sky radiance was included. Three or four component systems were considered, consisting of pure water, particulates with or without bubble clouds, and dissolved organic matter. The spectral absorption of dissolved organic matter was defined as

$$a(\lambda )=a(\lambda _0)\mathrm{exp}[0.014(\lambda \lambda _0)]$$
(7)
where $`a(\lambda _0)=0.1`$ $`\mathrm{m}^1`$, $`\lambda _0=440`$ nm (Bricaud et al., 1981). The phase function of phytoplankton was defined as an average of Petzold’s clear ocean, coastal ocean, and turbid harbor cases (Mobley, 1994; Tyler, 1977). The bubble cloud phase function was calculated for the real relative refractive index $`m=0.75`$ as an average for the size parameter range between $`x=100300`$ with a resolution $`\mathrm{\Delta }x=1`$ using Wiscombe (Fiedler-Ferrari et al., 1991) code. Figure 4 compares Petzold and bubble phase functions.

The main reason for the size distribution average was to remove transient spikes. This normalized phase function is scaled by $`b`$ calculated from Eq. 3. The volume scattering and absorption coefficients for particulates were determined from the Gordon-Morel model (Mobley, 1994; Gordon and Morel, 1983; Morel, 1988)

$$b_p=\frac{550}{\lambda }0.3C^{0.62}$$
(8)
$$a^{\mathrm{Case1}}(\lambda )=[a_w(\lambda )+0.06a_c^{}(\lambda )C^{0.65}][1+0.2\mathrm{exp}(0.014(\lambda 440))]$$
(9)
Here $`a_w(\lambda )`$ is the absorption coefficient of pure water, $`a_c^{}`$ is the chlorophyll-specific absorption coefficient, and $`C`$ is the chlorophyll concentration in $`\mathrm{mg}\mathrm{m}^3`$ (Mobley, 1994). The chlorophyll concentration was set as constant (well-mixed) with depth and equal to 0.8$`\mathrm{mg}\mathrm{m}^3`$ or 0.08$`\mathrm{mg}\mathrm{m}^3`$. We used Case 1 water parameterization but assumed that not all dissolved organic material is correlated with the chlorophyll concentration. In an apparent contradiction, the wind speed which defined surface reflectance and transmittance functions due to the wind-blown water surface was set to 0. The reason for this was to estimate the effect on remote sensing properties of the sub-surface bubble clouds. However, we investigated the sensitivity of the reflectance to change in the wind speed between 0 and 10 $`\mathrm{ms}^1`$, and the effect was small compared to the influence of bubbles. The azimuth direction was divided into 24 equally spaced sectors, the zenith-nadir range was divided into 20 equally spaced sectors. The profile of the bubble cloud volume fraction was determined by the expression

$$f(z)=f_0(f_1/f_0)^p$$
(10)
where $`f_0=f(z_0)=10^6`$, $`f_1=f(z_1)=10^7`$, $`z_0=0`$ m, $`z_1=8`$ m, and $`p=(zz_0)/(z_1z_0)`$. Our choice of the speed of attenuation is perhaps unreasonably gradual, except for the case of stable microbubbles. On the other hand this choice is of secondary importance for the remote sensing properties of bubble clouds which are dominated by the surface volume fraction. In-water asymptotic light field is discussed by Flatau et al. (1999). The scattering was conservative. The volume attenuation was determined from the asymptotic expression (3). The probability distribution function for scattering was discretized and was set as constant throughout the depth range. We calculated the asymmetry parameter to be $`g0.85`$. This parameter defines the probability of photon scattering towards the forward or backward hemisphere. It is equal to 1 if all photons are scattered forward and to -1 if all photons are scattered backward.

### 3.2 Remote sensing reflectance

The reference model runs were performed with the 3 component system of phytoplankton, pure water, and DOM with a constant chlorophyll profile. Two cases were computed (but only one is presented), representative of clear coastal ($`C=0.8`$) and oceanic ($`C=0.08`$) water. The remote sensed reflectance is defined as

$$R_{\mathrm{rs}}(\lambda )=\frac{L_w(\lambda )}{E_d(\lambda )}$$
(11)
where $`E_d`$ is the downwelling irradiance onto the sea surface, and $`L_w(\lambda )`$ is the upwelling water-leaving radiance. The remote-sensing radiance is a measure (Mobley, 1994) of how much of the downwelling light that is incident onto the water surface is returned into the zenith direction. The remote sensing reflectances are plotted in Fig. 5. The physical importance of the remote-sensing reflectance is evident from asymptotic theories which relate $`R_{\mathrm{rs}}`$ to inherent optical properties (Zaneveld, 1995)

$$R_{\mathrm{rs}}\frac{\beta (\pi \theta _m)}{a}$$
(12)
where $`\beta `$ is the phase function, $`\theta _m`$ is related to the sun zenith angle, and $`a`$ is the absorption coefficient. Thus, $`R_{\mathrm{rs}}`$ is approximately proportional to the probability of back- or side-scattering, and inversly proportional to the absorption of water column. Figure 5 shows the remote-sensing reflectance for the 3- and 4-component system with and without ocean bubbles but with constant pigment amount. The total single scattering albedo is strongly influenced by scattering from air bubbles. This leads to enhanced reflectance at all wavelengths. The results for $`C=0.08`$ (not presented) show even larger sensitivity. In Fig 5 the gray rectangles indicate bands (wavelenghts) which are used by the current ocean color satellite instrument (SeaWiFS). It is of interest to comment on the remote sensing of pigments and bubble cloud retrievals. Consider an algorithm based on the ratio of remote-sensing reflectance and define the ratio of remote-sensing reflectances without (Chl) and with (Chl+b) microbubbles as

$$\mathrm{ratio}(\lambda )=R_{\mathrm{rs}}^{\mathrm{Chl}}/R_{\mathrm{rs}}^{\mathrm{Chl}+\mathrm{b}}.$$
(13)
Figure 6 shows $`\mathrm{ratio}(\lambda )`$. Performance of the pigment algorithms based on the ratio of reflectances will depend on $`\mathrm{ratio}(\lambda _1)/\mathrm{ratio}(\lambda _2)`$. Submerged microbubble clouds seem to be wavelength-selective and even the ratio algorithms may require slight systematic correction. Given the increased sensitivity of the current generation of ocean color instruments, the absolute value of the radiances at the top of the atmosphere can be used for pigment retrievals.

## 4 Summary

Our calculations indicate that the optical effects of submerged microbubbles on the remote sensing reflectance of the ocean color are significant. These results are of importance for the retrievals of pigments from the ocean color measurements and for studies of the energetics of the ocean mixed layer. We provide information on how to reduce the systematic error due to microbubbles in pigment retrieval schemes via the $`\mathrm{ratio}(\lambda )`$. We also derive apparent optical property of bubbles - remote sensing reflectance - for the whole solar spectrum. This AOP is directly observable by the satellites and remote sensors. We expect that these and similar AOPs will have to be invoked in the case of hyperspectral retrievals for Case 2 waters where the signals from minerals, bubbles, chlorophyll, and dissolved organic material (CDOM) are not well correlated. New algorithms for current satellite instruments such as MODIS and SeaWIFS should employ this information.

We also present a global map of the volume fraction of air in water derived from daily wind speed data. We expect that such a map can be improved by knowing the day-to-day variability of the wind-wave relationship and better estimates of the volume fraction.

The paper is exploratory. Therefore, it is perhaps worth playing advocatus diaboli and speculate why the bubble clouds may not be important, at least in current satellite ocean color retrieval practice. Here are some reasons: (1) The high wind and clouds are correlated. This masks (bias) the effect of bubbles as observed from the satellites; (2) Both whitecaps (Wu, 1988b; Gordon and Wang, 1994; Frouin et al., 1996) and bubble clouds are correlated via their dependence on wind speed. Therefore, our results, as well as the hypothesis of Frouin et al. (1996), indicate that reflectance of foam has to be considered together with the reflectance due to bubble clouds. On the other hand, there are cases in which strong wind is not correlated with clouds. For example, the cross equatorial flow during the summer monsoon in the southern Indian Ocean is strong but the ITCZ position is in the northern hemisphere.

It is interesting to note that stabilized, coated microbubbles are hypothesized to be correlated to phytoplankton and CDOM concentrations; we need parameterization of this process. The optical properties of the first several meters below the surface are difficult to measure and are often removed from data due to experimental problems such as ship shadow or wave activity. This is the region where more detailed studies are needed.

## Acknowledgements

P. J. Flatau was supported in part by the Office of Naval Research Young Investigator Program and NASA SIMBIOS program. M. Flatau acknowledges NOAA/UCAR Global Climate Change Fellowship and J. R. V. Zaneveld acknowledges support of the Environmental Optics program of the Office of Naval Research and the Biogeochemistry program of NASA. C. D. Mobley acknowledges support of the Environmental Optics program of the Office of Naval Research, which also supported in part the development of the Hydrolight model.

## References


format: text

text:
# 1 Introduction

## 1 Introduction

In the past 3 years the Rossi X-ray Timing Explorer has discovered kilohertz quasi-periodic oscillations (kHz QPOs) in the persistent flux of 19 low-mass X-ray binaries (LMXBs). The power density spectra of most of these sources show twin kHz peaks. Initial results seemed to indicate that the frequency separation between the twin kHz peaks, $`\mathrm{\Delta }\nu =\nu _2\nu _1`$, remained constant even as the peaks gradually moved up and down in frequency (typically over a range of several hundred Hz) as a function of time. But thanks to precise measurements of the frequencies of the kHz QPOs, we know now that in several sources $`\mathrm{\Delta }\nu `$ is not constant, but it decreases as $`\nu _1`$ and $`\nu _2`$ increase.

In this chapter I describe a new technique that we have been using in the past few years to get precise measurements of the frequency separation of the kHz QPOs in some LMXBs. My plan is to show how this technique (that we call “shift-and-add”) works, and to present some of the results we obtained using it. It is not my purpose, however, to discuss here the details of the kHz QPO phenomenon, as this subject is extensively covered elsewhere in this book (see the chapter by van der Klis).

## 2 The Average Power Spectrum

The standard Fourier techniques in X-ray timing are described in detail in van der Klis (1989a); let us just recall here the definition of the power spectrum. A time series $`x(t)`$ is first divided into $`N`$ time intervals of length $`\tau `$, and for each interval the Fourier transform is calculated as:

$$X_n(\nu )=_{t_0+n\tau }^{t_0+(n+1)\tau }x(t)e^{2\pi i\nu t}𝑑t,n=0,1,2,\mathrm{},N,$$
(1)
where $`i=\sqrt{1}`$, and $`t_0`$ is some initial time. The power spectrum of $`x(t)`$ is defined as the average of the square of the Fourier transforms, $`P(\nu )=<X_n^2(\nu )>`$. The average is computed to improve the signal-to-noise ($`S/N`$) ratio of the power spectrum (see van der Klis 1989a for details). This definition assumes that the time series $`x(t)`$ is stationary, i.e. that all time segment of length $`\tau `$ have the same statistical properties.

As I already mentioned, one property of the kHz QPOs is that their frequencies vary with time (therefore the underlying time series is non-stationary). This is shown in Figure 1, for the source 4U 1728–34. On the left panel I show a dynamical power spectrum, where the x-axis represents time (the origin of time in this case is 1997 October 1 at 06:09 UTC), the y-axis is the Fourier frequency, and the different gray levels indicate the power in each frequency bin. The dark feature that crosses the Figure almost horizontally is one of the two QPOs detected in this source. The right panel shows power spectra ($`P(\nu )`$ vs. Fourier frequency) at some selected times: The same QPO is seen moving between $`700`$ and 900 Hz.

Obviously, the average power spectrum of the whole observation will show several peaks, at the frequencies at which the QPO stays for a while during the observation. This is shown in Figure 2 left panel, where the average power spectrum has been calculated as described above. However, as the QPO is strong enough to allow us to measure it in short time segments, we can think of “shifting” the frequency axis of each individual power spectrum before we calculate the average so that the new power spectrum is $`\stackrel{~}{P}(\nu )=<X_n^2(\nu \nu _{0,n})>`$. Here $`\nu _{0,n}`$ is the frequency shift that we need to apply to the power spectrum of each segment if we want to align the strong QPO at the same frequency before we calculate the average. Notice that this is equivalent to multiplying the function $`x(t)`$ by $`e^{2\pi i\nu _{0,n}t}`$ before calculating the integral in eq.(1).

The right panel in Figure 2 shows the result of such procedure (which, of course, also implies measuring $`\nu _{0,n}`$ for each segment; certainly this is extra work, but it comes with some rewards, as I will explain in §3). If it were only for the strong QPO in the average power spectrum (this QPO is stronger in the new power spectrum than in the original one) we would have not gained much by this procedure. Although there is now only one sharp peak, this shifted-and-averaged power spectrum provides no new information about the strong QPO. However, there is a second (weaker, but yet significant) QPO in the shifted power spectrum, $`350`$ Hz above the strong QPO. Although there seems to be evidence for such a weak QPO in the original power spectrum (Fig. 2, left panel), there it was not significant enough.

To understand why we detect this QPO in the shifted power spectrum, while in the standard power spectrum it was not significant enough, we need to recall that the $`S/N`$-ratio of the QPO is proportional to $`(T/W)^{1/2}`$, where $`T=N\times \tau `$ is the total length of the observation, and $`W`$ is the width of the QPO (see van der Klis 1989b). If the frequency separation between the 2 QPOs is roughly constant111Although it is not really constant, variations in $`\mathrm{\Delta }\nu `$ are smaller than changes in the centroid frequencies of the QPOs. (a close inspection to Fig. 1 left panel shows that there might be such weak QPO in the dynamic power spectrum, “following” the strong QPO as this one moves in frequency), the shift applied to the power spectra of the individual segments to align the strong QPO also aligns the weak QPO. This makes this new QPO narrower, and therefore more significant, in $`\stackrel{~}{P}(\nu )`$ than in $`P(\nu )`$.

This is how we first discovered the second QPO during 1996 outburst of 4U 1608–52 (Méndez et al. 1998a), which we afterwards detected again (this time without requiring the sensitivity-enhancing technique described above) during another outburst in 1998 (Méndez et al. 1998b).

Not surprisingly, this same technique can be used to obtain much better measurements of the frequency separation between the kHz QPOs. This is because the error in the centroid frequency of a QPO is proportional to the width of the QPO, and inversely proportional to its S/N-ratio (see Downs & Reichley 1983 for a similar situation when measuring the arrival time of individual pulses in radio pulsars). As described above, the weak QPO is usually narrower in the shifted power spectrum, $`\stackrel{~}{P}(\nu )`$, than in the standard power spectrum, $`P(\nu )`$, and therefore its frequency can be measured more accurately in $`\stackrel{~}{P}(\nu )`$ than in $`P(\nu )`$.

One way of measuring $`\mathrm{\Delta }\nu `$ using the shifted power spectrum is the following: Divide the data into $`N`$ segments of length $`\tau `$, and for each segment calculate the power spectrum. Measure the frequency of the strong QPO, $`\nu _1`$, in each of these $`N`$ power spectra, and group the data in $`M`$ sets such that $`\nu _1`$ is more or less constant within each set. This ensures that $`\mathrm{\Delta }\nu `$ is more or less constant within each set222Notice that there is a trade-off here: the choice of small frequency intervals for $`\nu _1`$ ensures small changes in $`\mathrm{\Delta }\nu `$, preventing artificial broadening of the weak QPO when the power spectra are shifted; but it also reduces the amount of data that are averaged in each set, and therefore decreases the $`S/N`$-ratio of the weak QPO at $`\nu _2`$. (early results in Sco X–1 showed that $`\mathrm{\Delta }\nu `$ decreased monotonically as $`\nu _1,\nu _2`$ increased; van der Klis et al. 1997). For each set, align the power spectra of the individual time segments using $`\nu _1`$ as a reference, and combine them to produce an average power spectrum. The result is $`M`$ (aligned) power spectra for which both kHz QPOs are as narrow as possible. For each of these $`M`$ power spectra measure $`\nu _1`$ and $`\nu _2`$, calculate $`\mathrm{\Delta }\nu =\nu _2\nu _1`$, and plot $`\mathrm{\Delta }\nu `$ vs. the average value of $`\nu _1`$ in each set.

Figure 3 shows the results of applying such procedure to the kHz QPOs in 4U 1608–52, 4U 1728–34, and Sco X–1. The implications of these results (specially in the case of 4U 1728–34) to the existing models for the kHz QPOs are presented elsewhere in this book (see the chapter by van der Klis), and will not be discussed here.

## 3 Other Results

As I mentioned in §2, there is a reward for the extra work of measuring $`\nu _1`$ in the individual time segments: We can combine these measurements with, for instance, the source count rate or X-ray colors in those same segments. One example of this is shown in Figure 4. The left panel shows a plot of $`\nu _1`$ vs. the $`216`$ keV count rate for 4U 1608-52. Each point represents 128 s of data. For the same data used in the plot on the left panel, the right panel shows the relation between $`\nu _1`$ and an X-ray color (see the caption of this figure for the definition of the color).

Interesting in this Figure is the complex dependence of $`\nu _1`$ upon X-ray count rate, which is usually assumed to be a good measure of mass accretion rate $`\dot{M}`$, and how this complexity is reduced to a single track in the frequency vs. hard color diagram. As before, the interested reader can find an extensive discussion of this subject elsewhere in this book. I just want to say here that these are “side” results of the technique described in §2.

## 4 Conclusion

I presented a new technique aimed at increasing the sensitivity to weak moving QPOs in the power density spectra of LMXBs. I showed some examples of the results obtained by applying this technique to data obtained with the Rossi X-ray Timing Explorer. In particular, this technique provides more precise measurements of the frequencies of the kHz QPOs, which can be used to set constraints on the models so far proposed to explain this new phenomenon.

###### Acknowledgements.

This work was supported by the Netherlands Research School for Astronomy (NOVA), the Netherlands Organization for Scientific Research (NWO) under contract number 614-51-002, the Leids Kerkhoven-Bosscha Fonds (LKBF), and the NWO Spinoza grant 08-0 to E.P.J. van den Heuvel. MM is a fellow of the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.


text:
# Self-organization in systems of self-propelled particles

Self-organization and pattern formation in systems of self-propelled entities are ubiquitous in nature. Examples can be found in a variety of fields and include animal aggregation , traffic patterns  and cell migration . Recently, the problem of flocking, in which a large number of moving particles (e.g. fish or birds) remain coherent over long times and distances, has attracted considerable attention. A simulation of a simple numerical model by Vicsek et al.  in which each particle has a constant identical speed and a direction of motion that is determined by the average direction of its neighbors, revealed that an ordered phase exists, even in the present of noise and disorder. A subsequent analysis of a continuum model by Toner and Tu  investigated this ordered phase further and derived conditions for its stability.

These models have in common that the flocks have infinite extent and, in simulations, fill the entire computational box. In reality, however, flocks have a finite size, with its density dropping sharply at the edge of the flock . In this Letter we present a discrete model consisting of self-propelled interacting particles that obey simple rules. We show that self-organization in our model leads to coherent localized states in one dimension (1D) and in two dimensions (2D) that are stable in the presence of noise and disorder. Furthermore, we present a continuum version of our discrete model which is obtained by coarse-grain averaging the discrete equations. The continuum flock solutions in 1D agree very well with the discrete solutions and are characterized by having a finite extent with densities that drop off sharply at the edges. In 2D, we focus on a vortex state in which the particles are rotating around a common center and show that the discrete model and the continuum model agree well.

Our discrete particle model consists of $`N`$ particles with mass $`m_i`$, position $`\stackrel{}{x}_i`$ and velocity $`\stackrel{}{v}_i`$. Each particle experiences a self-propelling force $`\stackrel{}{f_i}`$ with fixed magnitude $`\alpha `$. To prevent the particles from reaching large speeds, a friction force with coefficient $`\beta `$ is introduced. In addition, each particle is subject to an attractive force which is characterized by an interaction range $`l_a`$. This force is responsible for the aggregation of the particles and corresponds in animal aggregates to an awareness of the position of surrounding animals. To prevent a collapse of the aggregate, a shorter-range repulsive force is introduced with interaction range $`l_r`$. Thus, the governing equations for each particle is

$`m_i_t\stackrel{}{v}_i=\alpha \widehat{f_i}\beta \stackrel{}{v}_i\stackrel{}{}U`$ (1)
$`_t\stackrel{}{x}_i=\stackrel{}{v}_i`$ (2)
We have checked that our qualitative results are independent of the explicit form of the interaction potential and we have chosen here an exponentially decaying interaction:

$`U={\displaystyle \underset{ji}{}}C_a\mathrm{exp}(|\stackrel{}{x}_i\stackrel{}{x}_j|/l_a)`$ (3)
$`{\displaystyle \underset{ji}{}}C_r\mathrm{exp}(|\stackrel{}{x}_i\stackrel{}{x}_j|/l_r)`$ (4)
where $`C_a`$, $`C_r`$ determine the strength of the attractive and repulsive force respectively. The direction of the self-propelling force can be chosen along the instantaneous velocity vector or, similar to the numerical model of Ref. , can be determined by aligning it with the average velocity direction of the neighboring particles:

$`\stackrel{}{f}_i=\stackrel{}{v}_i`$ $`\mathrm{without}\mathrm{averaging}`$ (5)
$`\stackrel{}{f}_i={\displaystyle \underset{ji}{}}\stackrel{}{v}_j\mathrm{exp}(|\stackrel{}{x}_i\stackrel{}{x}_j|/l_c)`$ $`\mathrm{with}\mathrm{averaging}`$ (6)
where $`l_c`$ is a correlation length.

Let us now present our simulation results of the discrete model. The model was integrated by solving Eqs. 1,2 using a simple Euler integration routine with timestep $`\mathrm{\Delta }t=0.2`$. The simplest coherent localized solution in our model is a 1D flock in which all particles move with constant velocity $`v=\alpha /\beta `$. An example of this solution is presented in Fig. 1 where we have plotted, as solid circles, the density defined as $`\rho _i=2/(x(i+1)x(i1))`$ as a function of the position of the particle. The density can be seen to drop abruptly to zero at the edge of the flock. We have checked that this solution is stable in the presence of moderate amounts of noise (added to Eq. 2) and of disorder in the parameters.

Further simulations revealed that increasing the number of particles does not change the shape of the density function and that the total size of the flock reaches a constant value. This is illustrated in Fig. 2 where the dashed line represents the size of the flock as a function of the number of particles. This obviously unrealistic behavior of the model is due to the soft-core repulsive force which allows the particles to be very close. Our model can easily be extended to incorporate hard-core repulsive forces. In fact, a flock with a size that scales approximately linearly with the total number of particles can be obtained by introducing a hard-core repulsive force in Eq. 1 . The specific form of the hard-core potential is not important and we have chosen

$`U_{hc}={\displaystyle \underset{ji}{}}C_{hc}(|x_ix_j|l_{hc})^5|x_ix_j|l_{hc}`$ (8)
$`U_{hc}=0|x_ix_j|>l_{hc}`$ (9)
In Fig. 2 we show, as a solid line, the size of the flock in the presence of this additional repulsive force as a function of $`N`$ while in the inset we show the corresponding density of the flock for different $`N`$. Naturally, the force has only an effect when the interparticle separations are smaller than $`l_{hc}`$. Hence, for small $`N`$ the flock solution is unaffected by the additional force. As $`N`$ is increased and the interparticle spacing becomes smaller than $`l_{hc}`$, the particles in the center of the flock are pushed apart. For large $`N`$, the flock reaches a constant density in its center and its size scales linearly with $`N`$.

Let us now turn to 2D where we have obtained several flocking states. One, not shown here, is the equivalent of our 1D flock: all particles are moving in the same direction with $`|\stackrel{}{v_i}|=\alpha /\beta `$. The particles arrange themselves on a disk and this aggregate is stable under small amounts of noise and disorder. The solution which we will focus on here consists of a vortex state where the particles rotate around a common center and which is common in fish schools , bacterial colonies  and amoeba aggregates . This solution has been observed previously in models of self-propelling particles but only in the presence of a confining boundary  or when a rotational chemotaxis term is invoked .

In our model, the vortex solution can be obtained from a wide variety of initial conditions including one in which all particles are randomly placed on a disk with speed $`\alpha /\beta `$ and random initial velocities. A typical evolution of the particles starting from this initial condition is shown in Fig. 3. This figure was obtained in the absence of the velocity averaging term and illustrates that in this case some particles move clockwise while the others rotate counter clockwise. When the velocity averaging term is included the final vortex state consists of all particles turning the same way. An example of this case is shown in the inset of Fig. 4. In both cases, the speed of each particle was found to be sharply peaked around $`|\stackrel{}{v}|=\alpha /\beta `$

The average size of the vortex remains constant in time as shown in Fig. 4 where we have plotted the average density (obtained by averaging over $`10^6`$ iterations) of the vortex structure. Fig. 4 displays several remarkable features. First, there is a well-defined core which remains void even for extended simulation runs. Using different parameter values however, one can also produce a vortex without a core. Second, as in the one dimensional flock, the density does not decay smoothly to zero at the edges. Instead, it increases at both the inner and outer edge of the aggregate and then drops abruptly. Qualitatively similar vortex solutions were found when an additional hard-core repulsion like the one discussed above is added.

As in the case of traffic models , it is useful to develop a continuum version of our model. To this end, we simply coarse-grain average the equations which results in, after dividing by a common factor of $`\rho `$,

$`_t\stackrel{}{v}+(\stackrel{}{v}\stackrel{}{})\stackrel{}{v}=\alpha \widehat{f}\beta \stackrel{}{v}+\stackrel{}{G}`$ (10)
together with the conservation equation for the density $`\rho `$

$`_t\rho +\stackrel{}{}(\stackrel{}{v}\rho )=0.`$ (11)
The interaction force is given by

$`\stackrel{}{G}(\stackrel{}{x})={\displaystyle \rho (\stackrel{}{x}^{})\stackrel{}{}U(\stackrel{}{x},\stackrel{}{x}^{})𝑑\stackrel{}{x}^{}}`$ (12)
and the self-propulsive force direction is given by either

$`\stackrel{}{f}(\stackrel{}{x})={\displaystyle \rho (\stackrel{}{x}^{})\stackrel{}{v}(\stackrel{}{x}^{})exp(|\stackrel{}{x}\stackrel{}{x}^{}|/l_d)𝑑\stackrel{}{x}^{}}`$ (13)
in the velocity averaging case or simply by $`\widehat{f}=\stackrel{}{v}/|\stackrel{}{v}|`$ otherwise.

A comparison between the discrete model and the continuum model can be carried out for the solutions presented here. For a 1D flock, we have simply $`\widehat{f_i}=1`$, and the solution in the continuum model is given by $`G=0`$ or

$$\rho (x^{})U(x,x^{})𝑑x^{}=D$$
(14)
where $`D`$ is a constant determining the total number of particles. Since the sought after solution has a finite extent with a discontinuity at the edge where the density drops to zero we discretize the integral using $`M`$ points and discretization step $`\mathrm{\Delta }x`$. The last point corresponds to the edge of the flock. The resulting linear set of $`M`$ equations for $`\rho `$ is easily solved using standard linear algebra packages. and $`\mathrm{\Delta }x`$ was varied until the slope at the center of the flock vanished. The result, with $`D`$ chosen such that $`\rho (x)𝑑x=N`$, is shown in Fig. 1 as a solid line. The density profile in the continuum model is discontinuous at the edge and agrees well with the profile obtained in the discrete model. Note that since the equations are linear in $`\rho `$ it is not surprising to find that the density profile of the discrete model did not change as the number of particles is increased. Clearly,the simple coarse-grained averaging procedure is not adequate for the hard-core potential case, where higher order terms in the density are important.

The continuum model can also be used to find the vortex solution. To this end, we use the fact that all particles undergo approximately a circular motion with constant speed ($`\alpha /\beta `$). Thus, a continuum vortex solution can be found by requiring that the force $`\stackrel{}{G}`$ is centripetal:

$$_0^{2\pi }𝑑\varphi _0^{\mathrm{}}𝑑r\rho (r)U(r,\varphi )=D(\alpha /\beta )^2\mathrm{ln}(r)$$
(15)
After performing the integration over $`\varphi `$, the remaining integral was discretized as in the 1D case. The resulting matrix was solved for $`\rho (r)`$ and used in a Newton solver that searched for the size of the hole and the overall size of the vortex (i.e. discretization $`\mathrm{\Delta }x`$) with a condition for smoothness of the solution at both discontuous edges. This condition simply consisted of requiring that the first and last point can be obtained by linear interpolation from its two neighboring points. We have checked that the solutions we obtained are converged by increasing the number of discretization points from 80 to 1480. In Fig. 4 we compare the discrete solution to the one obtained by Eq. 15 where the integration constant $`D`$ was varied until $`\rho (\stackrel{}{x})𝑑\stackrel{}{x}=N`$. Again, the continuum profile is discontinuous at the edges and the agreement between the continuum profile and the discrete profile is very good.

The continuum equation can be used to explore the (large) parameter space more efficiently. An example of such an exploration is shown in Fig. 5 where we plot two different solutions found by our Newton solver for the same model parameters and total number of particles ($`N`$=200) but with different integration constant $`D`$. Preliminary simulations of the discrete model show that the solution with the larger core is unstable and a formal stability analysis of the continuum solution will be carried out in the future.

In this Letter we have presented a model for localized aggregates and flocks. Our model contains very simple rules and can be straightforwardly extended to incorporate additional and different types of interactions. For example, the forces that maintain bacterial aggregates are believed to be short-range adhesion forces together with a short-range hard-core repulsion. The investigation of these types of interactions will be the subject of future work. Finally, it would be interesting to compare our results to animal flocks. Unfortunately, such a comparison is currently not possible since not enough quantitative data is available.

We acknowledge useful conversations with E. Ben-Jacob and W.F. Loomis. The work of HL and WJR was supported in part by NSF DBI-95-12809. IC acknowledges support from The Colton Scholarships and a Israeli-US Binational Science Foundation BSF grant.


text:
# Constraining New Physics with Vertex Corrections

## Abstract

We discuss how new physics not encompassed within the $`STU`$ oblique correction framework can be constrained from precision electroweak measurements via vertex corrections to $`Z`$–pole observables.

## 1. Limitations of the Oblique Correction Analysis

In constructing models beyond the Standard Model (SM) such as technicolor or supersymmetric theories, it is important to test their viability by looking for possible conflicts with currently available experimental data. A popular procedure has been the $`STU`$ oblique correction analysis of Ref. . However, the method can only be applied to theories which satisfy three conditions, namely (1) the electroweak gauge group is the standard $`SU(2)_L\times U(1)_Y`$, (2) vertex and box corrections from new physics are negligible, and (3) the scale of new physics is large compared to the electroweak scale. Furthermore, the oblique correction analysis requires (4) the model to be completely specified since any particle with electroweak gauge quantum numbers will show up in a vacuum polarization loop. Many theories of interest do not fall into this category. For instance, in the topcolor assisted technicolor model of Ref. , (1) the electroweak gauge group is $`SU(2)_L\times U(1)\times U(1)`$ resulting in an extra $`Z^{}`$ boson, (2) there are potentially large corrections to the $`Zf\overline{f}`$ vertex for third generation fermions coming from $`Z^{}`$ and coloron exchange, (3) the scale of new physics may be as low as a TeV, and (4) the technisector of the theory remains unspecified.

## 2. Vertex Corrections at the $`Z`$–pole

The solution to our predicament lies in the observation that a large number of electroweak observables have been measured accurately at LEP and SLD on the $`Z`$–pole, and that the majority of them are just ratios of coupling constants of the $`Z`$ to quarks and leptons. This means that they are only sensitive to oblique corrections through the effective value of $`\mathrm{sin}^2\theta _W`$ at that particular energy scale. Therefore, one can use just one of the observables, $`A_{\mathrm{LR}}`$ say, to fix the value of $`\mathrm{sin}^2\theta _W`$, and use it to predict all the rest for the SM. Any deviations of the data from those predictions must be due to vertex corrections from new physics.

This method of comparing SM predictions to $`Z`$–pole electroweak observables solves all of the problems mentioned above since (1) mixing of the $`Z`$ with extra gauge bosons can be treated as a vertex correction, (2) vertex corrections need not be small while box corrections are still naturally suppressed on the $`Z`$–pole, (3) the new physics scale need not be large since only observables at one scale are used in the analysis, and (4) the theory need not be specified completely since the method is blind to oblique corrections.

## 3. Applications

This method has been used successfully to constrain the size of possible vertex corrections to the $`Zb\overline{b}`$ vertex , and to constrain the topcolor assisted technicolor model mentioned above . We have also applied the technique to constrain: an R–parity violating extension to the MSSM , Higgs masses in a general two Higgs doublet model , and the fundamental Planck scale of a model with large extra dimensions . Details are provided in the respective references.

## Acknowledgements

This work was supported in part by the U.S. Department of Energy, grant DE–FG05–92–ER40709, Task A.


text:
# Phase Transition and Acoustic Localization in Arrays of Air-Bubbles in Water

## Figure Captions

Averaged transmission as a function of $`ka`$.

Left column: The phase diagram for the two dimensional phase vectors defined in the text for one random distribution of bubbles. Right column: The acoustic energy distribution or transmission, averaged for all directions, as a function of distance away from the source; the distance is scaled by the sample size $`R`$. (Original in Color)

The acoustic band structure for a simple cubic array of air-bubbles. The complete band gap lies between the two dashed horizontal lines. The plot is rendered in terms of the non-dimensional $`ka`$ as a function of Bloch vectors.

The phase diagram for the two dimensional phase vectors and the energy distribution, averaged for all directions, for the cubic array of the bubbles which has the same bubble volume fraction as the random case. The lattice constant $`d=1/(\beta )^{1/3}a`$. (Original in Color)


text:
# Comment about the Letter entitled “Scalar fields as dark matter in spiral galaxies”

## Acknowledgments

This work was in part supported by DGAPA-UNAM grant No IN121298 and by CONACyT grant 32272-E.


text:
# Local properties of extended self-similarity in 3D turbulence

## Abstract

Using a generalization of extended self-similarity we have studied local scaling properties of incompressible homogeneous isotropic 3D turbulence in a direct numerical simulation. We have found that these properties are consistent with lognormal-like behavior of the velocity increments with moderate amplitudes for space scales $`r`$ beginning from Kolmogorov length $`\eta `$ up to the largest scales, and in the whole range of the Reynolds numbers: $`50R_\lambda 459`$. The locally determined intermittency exponent $`\mu (r)`$ varies with $`r`$; it has a maximum at scale $`r=14\eta `$, independent of $`R_\lambda `$.

So-called extended self-similarity (ESS) in incompressible turbulent flows is intensively studied in recent years (see, for instance,  and references therein). The ESS implies scaling relation between moments of different order. For example, for absolute value of longitudinal velocity increments over separation $`r`$ in the inertial range the ESS means

$$u_r^qu_r^p^{\rho (q)/\rho (p)},$$
(1)
where the scaling exponent $`\rho (q)`$ is some function of $`q`$. It is shown in numerous experiments and numerical simulations that the range of applicability of the ESS is substantially larger than that for ordinary self-similarity, and the ESS can exist even for situations where the ordinary self-similarity cannot be observed at all.

The aim of the present Letter is twofold. The first is to study local properties of ESS defined as

$$u_r^pu_r^3^{\zeta _p(r)},$$
(2)
where $`\zeta _p(r)`$ are depending on $`r`$. The local approach allows to develop some old ideas. At past a representation $`\zeta _p=p/3\mu p(p3)/18`$ derived from a lognormal model, was shown to hold for the values averaged over an inertial interval of certain extension. Here we will show that the representation holds locally (i.e. with $`\zeta _p(r)`$ and $`\mu (r)`$ depending on $`r`$) in the above defined sense, ranging from the dissipative Kolmogorov length $`\eta `$ to integral scale for any observed Reynolds number in direct numerical simulations up to $`R_\lambda =459`$.

The second aim is to examine the implication of the variation of the local intermittency exponent $`\mu (r)`$. It is often mentioned that there is no characteristic length in turbulence, so that the structure functions obey a power law in the inertial region and the associated scaling exponents are independent of scale. However, a recent study of turbulence reveals the existence of the structures . A question, then, naturally arises as to whether there may be any characteristic length ascribed to the structures . The peculiar variation of $`\mu (r)`$ with respect to $`r`$, observed in the Letter, indicates that there is a certain length affecting the ESS.

The She-Leveque model  is also very popular in the last years. For this model the local exponent is expressed in a general form:

$$\zeta _p(r)=\frac{p}{3}(1\gamma )+\frac{\gamma }{1\beta }\left[1\beta ^{p/3}\right],$$
(3)
where $`\beta `$ and $`\gamma `$ can be also supposed to depend on $`r`$; $`\beta (r)`$ and $`\gamma (r)`$ could be in principle evaluated by plotting the local exponent $`\zeta _p(r)`$ against $`p`$ with $`r`$ fixed. Since here we focus the lower order structure functions which are computed more reliably than the higher order ones, the comparison of the data with the formula is done using a simpler lognormal model with one parameter. If one could compute the higher order structure functions, the comparison should be made with the She-Leveque model.

Let us begin with the derivation of useful formulae on the basis of the lognormal distribution of the dissipation rate $`\epsilon _r`$ averaged over spheres of radius $`r`$:

$$P(\epsilon _r)=\frac{\epsilon _r^1}{\sqrt{2\pi \sigma ^2}}\mathrm{exp}\left(\frac{(\mathrm{ln}\epsilon _ra)^2}{2\sigma ^2}\right),$$
(4)
from which we obtain

$$\frac{\epsilon _r^q}{\epsilon _r^q}=e^{\sigma ^2q(q1)/2}$$
(5)
that results in a parameter-independent type of ESS  of turbulent energy dissipation

$$\frac{\epsilon _r^q}{\epsilon _r^q}=\left(\frac{\epsilon _r^p}{\epsilon _r^p}\right)^{\frac{q(q1)}{p(p1)}}.$$
(6)
What is the equivalent relation for the velocity increment? According to the refined similarity method  $`u_r(r\epsilon _r)^{1/3}`$ in the inertial region, while $`u_rr\epsilon _r^{1/2}`$ in the dissipative region. If we take regions of scale $`r`$ between the inertial and dissipative regions, some are described by the former relation, while the others are by the latter one, so that $`u_r`$ is related to $`\epsilon _r`$ in a probabilistic way. For the sake of simplicity we assume a general relation

$$u_r=c(r)\epsilon _r^{1/\alpha (r)}$$
(7)
in a mean sense. Here $`\alpha (r)`$ is a function of $`r`$ and coefficient $`c(r)`$ does not necessarily scale with $`r`$. Substituting (7) into (6) yields

$$\frac{u_r^q}{u_r^\alpha ^{q/\alpha }}=\left(\frac{u_r^p}{u_r^\alpha ^{p/\alpha }}\right)^{\frac{q(q\alpha )}{p(p\alpha )}}.$$
(8)
Note that this expresses how the structure functions of different order are related to each other with value of $`r`$ fixed. This relation can be considered as a functional equation and a solution to this equation is

$$u_r^qu_r^{\alpha (r)}^{q/\alpha (r)+b(r)q(q\alpha (r))},$$
(9)
where $`b(r)`$ is an arbitrary function of $`r`$. Relation (9) is a generalization of ordinary ESS . The difference between ordinary ESS and relation (9) is that parameters used in (9) can depend on $`r`$. Therefore we will call this type of ESS as extended local self-similarity (ELSS).

To be consistent with the present data processing we express the $`q`$th order structure function in terms of 3rd order structure function as (2). Making use of (9), we are led to

$$\zeta _q(r)=\frac{q}{3}\frac{\mu (r)}{18}q(q3),$$
(10)
where

$$\mu (r)=\frac{6b(r)\alpha (r)}{1+\alpha (r)(3\alpha (r))b(r)}.$$
(11)
In the frame of ELSS the exponent $`\zeta _q`$ depends also on $`r`$, and below we compare (2) and (10) with data of DNS for different values of $`r`$. It should be noted that the ELSS expression (10) holds for any value of $`\alpha (r)`$, so that the expression can be compared with the data for any scale separation without paying attention to which region is being considered.

We have performed a series of direct numerical simulations (DNS’s) of incompressible homogeneous isotropic turbulence using a resolution up to $`1024^3`$. Reynolds numbers range from 50 to 459 . The random force is statistically homogeneous, isotropic and Gaussian white, and applied to the band $`1\stackrel{<}{}k\stackrel{<}{}3`$ in which the forcing spectrum is constant. The code uses the pseudo spectral method and the 4th order Runge-Kutta-Gill one. Initial conditions are Gaussian random velocity fields with the energy spectrum $`E(k)k^4\mathrm{exp}(2(k/k_0)^2)`$, and the resolution is $`N=256^3`$ for $`R_\lambda =69`$, $`N=512^3`$ for $`R_\lambda =125,176,259`$, and $`N=1024^3`$ for $`R_\lambda =374,459`$. After about two eddy turnover times all the turbulent fields attained statistically steady states, which were confirmed by observing the time evolution of the total energy and enstrophy, and the skewness of the longitudinal velocity derivative. For $`R_\lambda =459`$ run, the Reynolds number was gradually increased through two steady states. The condition $`k_{\mathrm{max}}\eta >1`$ for the resolution of DNS is satisfied for most runs, but that of $`R_\lambda =459`$ is slightly less than unity. The statistical averages were taken as the time average over tens of turnover times for lower Reynolds numbers and over a few turnover times for the higher Reynolds numbers; over 45 samples during 2.9 eddy turnover times for $`R_\lambda =374`$ and over 31 samples during 1.4 eddy turnover times $`R_\lambda =459`$. Computations with $`R_\lambda 259`$ have been done using a vector parallel machine with 16 processors, Fujitsu VPP700E at RIKEN, and those for higher $`R_\lambda `$, using Fujitsu VPP5000/56 with 32 processors at Nagoya University Computation Center.

Now turn to the data analysis. Figure 1 is a plot of $`u_r^2/(\overline{\epsilon }\eta )^{2/3}`$ against $`r/\eta `$ for various values of Reynolds number, where $`\overline{\epsilon }`$ is the average dissipation rate, and $`\eta =(\nu ^3/\overline{\epsilon })^{1/4}`$ with molecular viscosity $`\nu `$. Here a straight solid line proportional to $`r^{2/3}`$ is inserted. It is remarkable that all data points collapse on a single line in the dissipative region, which indicates that all simulations are carried out with the good resolution at small scales. Although the slope of $`u_r^2`$ could be estimated for large Reynolds numbers as seen from Fig.1, the scaling exponents of higher order structure functions as well as low order ones for small Reynolds numbers can be reliably evaluated only on the basis of the ESS method, i.e. by plotting $`u_r^p`$ against $`u_r^3`$ .

In order to know the $`r`$-dependence of $`\zeta _p(r)`$ for various Reynolds numbers, we prepare Fig.2, in which $`\zeta _p(r)`$ with $`p=4,6,8`$ are depicted for $`R_\lambda =69,125,259,374`$ and $`459`$. (The eighth order structure function is confirmed to converge statistically.) Note that the data for scales larger than integral scales are not shown, because a universal property of turbulence is not expected in those data. It is remarkable that there is a dip at about $`r/\eta 10`$, and that it grows in depth with Reynolds number. The exception is the case $`R_\lambda =69`$, where a dip does not appear. As the scale increases beyond the dip, $`\zeta _4(r)`$ and $`\zeta _6(r)`$ tend to approach constant values, although the corresponding data for $`R_\lambda =259`$ behave in a slightly different way from other cases. For $`p=8`$ the situation is the same as for $`p=6`$, but the variation is larger. It is of interest to notice that for the largest Reynolds number $`459`$ the flat region is observed in the interval $`100\stackrel{<}{}r/\eta \stackrel{<}{}300`$, which may be identified with the inertial region. For smaller Reynolds numbers it is a little difficult to find the flat region. For $`R_\lambda =69`$ we see the flat region in the interval $`10\stackrel{<}{}r/\eta \stackrel{<}{}30`$, and the corresponding slopes for $`p4`$ are larger than those for $`R_\lambda =459`$. However, the flat region at smaller scales for $`R_\lambda =69`$ is different from one at larger scales for $`R_\lambda =459`$.

Before we compare the formula (10) with the data, the validity of the lognormal distribution of the velocity increment should be ensured. The pdf of the intermediate amplitudes is certainly lognormal. In order to determine the range of the lognormality, we calculated a peak position $`u_r^{}(p)`$ of $`u_r^pP(u_r)`$ at $`r/\eta =24`$ for $`R_\lambda =121`$. The pdf is satisfactorily fitted by a lognormal curve for $`u_r^{}(p=1.5)u_ru_r^{}(p=6)`$. On the other hand, the pdf of $`\epsilon _r`$ is found to be lognormal in much wider interval. If we employ the same notation as above, the lognormality of $`\epsilon _r`$ holds at least in the interval of $`\epsilon _r^{}(p=4)\epsilon _r\epsilon _r^{}(p=6)`$ for the same Reynolds number. If the refined similarity hypothesis $`u_r^3r\epsilon _r`$ is assumed to hold for any amplitude, the corresponding pdf of $`u_r`$ should be lognormal for $`u_r^{}(p=12)u_ru_r^{}(p=18)`$, which is much wider than the observed lognormal interval. The reason for the discrepancy is that the refined similarity holds only for the intermediate amplitudes of $`\epsilon _r`$ in agreement with the observation. Hence, the use of the lognormal expression for the exponent (10) is completely justified for intermediate values of $`p`$.

In order to analyze a nature of the $`r`$-dependence of the local slope $`\zeta _p(r)`$, we rewrite (10) in the following form:

$$\frac{\zeta _p}{p}=\left(\frac{1}{3}+\frac{\mu (r)}{6}\right)\frac{\mu (r)}{18}p.$$
(12)
What is the range of $`p`$? To decide the range we calculated $`\zeta _p(r)`$ for various values of $`p`$ at units of 0.1 at $`r/\eta =24`$ for $`R_\lambda =121`$, and plotted $`\zeta _p/p`$ against $`p`$. Although such a plot is not given here, $`\zeta _p/p`$ is fitted by a straight line for $`1p6`$. The curve is deviated upward from the straight fitting line for $`p7`$ and downward for $`0<p<1`$. Hence the comparison will be made in the range $`1p6`$, but the intermittency coefficient $`\mu (r)`$ estimated there plays a significant role to represent even the whole intermittency effect.

Figure 3a shows a curve $`\zeta _p(r)/p`$ vs. $`p`$ for several values of $`r/\eta `$ with $`R_\lambda =69`$. Straight lines in this figure are the best fit lines for $`1p6`$. The data points for $`p=7`$, which are not included for the comparison, are slightly deviated from the lognormal lines as mentioned above. Intermittency index, $`\mu (r)`$, can be calculated from this figure using the slope of the fitting straight lines or the intersection point of the fitting straight lines with vertical axis (cf.(12)). The calculated values of $`\mu (r)`$ are shown in the inset of the Fig.3a ($`\mu _1`$ corresponds to calculations using the ’slope’ method, whereas $`\mu _2`$ are given by the ’intersection’ method.) In the plot of $`\mu (r)`$ Taylor microscale $`\lambda `$ and an integral scale $`L`$ are marked for convenience. (We have confirmed the well-known prediction  $`\lambda /\eta =15^{1/4}R_\lambda ^{1/2}`$ and $`L/\eta R_\lambda ^{3/2}`$.) For Reynolds numbers $`259`$ and $`459`$ we have obtained similar pictures ( Figs.3b and 3c). Other Reynolds numbers give the same result.

As seen from Figs.3a to 3c, $`\mu (r)`$ substantially depends on $`r`$ with a typical “two-maxima” shape. Even for $`R_\lambda =69`$ in Fig.3a there is a small peak around $`r/\eta 10`$. Local maximum of the $`\mu (r)`$ at smaller scales exhibits a few interesting properties. (I) For $`R_\lambda =50`$ the maximum is not observed, and the first appearance of the maximum occurs at $`R_\lambda `$ between 50 and 69. (II) A position of this maximum (normalized by $`\eta `$) is independent of Reynolds number, and takes value $`\mathrm{}_114\eta `$, which is actually the beginning of the dissipative region. (III) Value of $`\mu (\mathrm{}_1)`$ scales with $`R_\lambda `$ as $`\mu (\mathrm{}_1)R_\lambda ^{0.26}`$. (IV) The scale $`\lambda `$ is located on the right descending hill.

The flat region where $`\mu (r)`$ is constant appears in between $`r/\eta 100`$ and $`r/\eta 300`$ for the highest Reynolds number $`R_\lambda =459`$ (see Fig.3c). In such an interval $`\mu (r)0.25`$. This value is consistent with those known in literature for observations corresponding to very large Reynolds number .

It should be noted that strong tube-like vortices are believed to be sparsely distributed in space in fully developed turbulence . Those vortices are frequently assumed as Burgers’ vortex with mean radius $`10\eta `$ . The energy dissipation takes place strongly around those vortices. Therefore extreme intermittency of the energy dissipation at this scale is consistent with the observation of the maximum of $`\mu (r)`$ at scale $`r=14\eta `$. On the other hand, the usual scaling region, where $`\zeta _p(r)`$ as well as $`\mu (r)`$ are independent of scale in a certain interval of scale, can be seen only for largest Reynolds number, $`R_\lambda =459`$ in the present work. Such a scaling region is located for $`r>\lambda `$; in our simulation the scaling region starts from approximately $`2.5\lambda `$.

The work of T.G. presented here was supported by a Grant-in Aid for Scientific Research (C-2) 12640118 from the Japan Society for the Promotion of Science. We are very grateful to RIKEN Computer Center and Nagoya University Computer Center for their support. T. Ochiai at NIT is also acknowledged for his assistance in numerical computation involved in this work. One of the authors A.B. is grateful to C.H. Gibson and K.R. Sreenivasan for numerous discussions on the problem.


text:
# On Batchelor passive advection by a finite-time correlated random velocity field

## I Introduction

Magnetic field passively advected by an incompressible velocity field obeys the induction equation

$`_tB^i+v^k_kB^iB^k_kv^i=\eta \mathrm{\Delta }B^i,`$ (1)
where $`\eta `$ is the diffusivity. We will assume that the diffusivity is negligibly small (large Prandtl number), and the magnetic field is therefore frozen into the fluid. The incompressibility of the velocity field is not relevant for further calculations; for the compressible case, Eq. (1) holds for $`B^i/\rho `$, where $`\rho `$ is the density of the fluid. The random velocity field is assumed to be Gaussian in the Eulerian frame:

$`v^i(𝐱,t)=0,`$ (2)
$`v^i(𝐱,t)v^j(𝐱^{},t^{})=\kappa ^{ij}(𝐱𝐱^{},tt^{}),`$ (3)
where $`\kappa ^{ij}(𝐱,t)`$ is an arbitrary function regular at $`𝐱=0`$, this is the so-called Batchelor regime . The solution of the Eq. (1) is

$`B^i(x,t)=B_0^i+{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \frac{\stackrel{~}{v}^i(y,\tau )}{y^k}}\text{d}\tau B_0^k=B_0^i+{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \frac{\stackrel{~}{v}^i\mathbf{(}y(x,t),\tau \mathbf{)}}{x^k}}\text{d}\tau B^k\mathbf{(}x(y,t),t\mathbf{)},`$ (4)
where $`y^i=x^i(y,0)`$ and $`B_0^i=B^i(y,0)`$ are initial values of the fluid-particle coordinate and of the magnetic field. The object $`\stackrel{~}{v}^i(y,\tau )`$ in the integrand is the Lagrangian velocity field. The Eulerian velocity has the form $`v^i(x,\tau )=\stackrel{~}{v}^i\mathbf{(}y(x,\tau ),\tau \mathbf{)}`$. To make use of the formula (2), we need to “stop” the moving point $`x^i=x^i(y,t)`$ in the function $`\stackrel{~}{v}^i\mathbf{(}y(x,\tau ),\tau \mathbf{)}`$. This is done differently in the Eulerian and Lagrangian frames (see Fig.1).

## II Lagrangian frame

In the Lagrangian frame, $`𝐲`$ is fixed and we need to express $`\stackrel{~}{v}^i(y,\tau )`$ in terms of $`v^i(y,\tau )`$. The general expression is $`\stackrel{~}{v}^i(y,\tau )=v^i(y^k+\underset{0}{\overset{\tau }{}}\stackrel{~}{v}^k(y,\tau ^{}),\tau )`$. For small $`t`$, we expand this expression up to the second order in $`\underset{0}{\overset{\tau }{}}\stackrel{~}{v}^i(y,\tau ^{})`$, and iterate once with respect to $`\stackrel{~}{v}`$. We obtain

$`\stackrel{~}{v}^i(\tau )=v^i(\tau )+v_l^i(\tau ){\displaystyle \underset{0}{\overset{\tau }{}}}v^l(\tau ^{})+v_l^i(\tau ){\displaystyle \underset{0}{\overset{\tau }{}}}v_j^l(\tau ^{}){\displaystyle \underset{0}{\overset{\tau ^{}}{}}}v^j(\tau ^{\prime \prime })+{\displaystyle \frac{1}{2}}v_{lj}^i(\tau ){\displaystyle \underset{0}{\overset{\tau }{}}}v^l(\tau ^{}){\displaystyle \underset{0}{\overset{\tau }{}}}v^j(\tau ^{\prime \prime })+\mathrm{},`$ (5)
where all the velocities are taken at the point $`𝐲`$, and we use the short-hand notation $`v_l^i(\tau )v^i(y,\tau )/y^l`$ and $`v_{lj}^i(\tau )^2v^i(y,\tau )/y^ly^j`$.

To obtain the equation for the second moment of $`B^i`$, we substitute the expansion (5) into the formula (4), raise the latter to the second power and average using formula (2). We assume that $`t`$ is much smaller than the inverse growth rate of the second moment, and that the correlation time $`\tau _c`$ is much smaller than $`t`$. Assuming isotropic and spatially homogeneous initial distribution, $`B_0^i=0`$, $`B_0^iB_0^j=\frac{1}{d}\delta ^{ij}H_2(0)`$, we find the equation for the second-order moment $`H_2(t)=|𝐁|^2(t)`$ in the Lagrangian frame. For the incompressible velocity field, it takes the form

$`{\displaystyle \frac{\text{d}}{\text{d}t}}H_2={\displaystyle \frac{1}{d}}\kappa _{jj}^{ii}H_2{\displaystyle \frac{1}{24d}}\tau _c(\kappa _{jjlm}^{ii}\kappa ^{lm}8\kappa _{jl}^{im}\kappa _{jm}^{il}+10\kappa _{lm}^{ii}\kappa _{jj}^{lm})H_2.`$ (6)
To derive this equation we assumed that $`\kappa (y,t)\kappa (y)T(t)`$, where $`T(t)`$ is a $`\delta `$ function smeared over the correlation time $`\tau _c`$. In particular, we have chosen the “box” regularization: $`T(t)=1/\tau _c`$ for $`t[\tau _c/2,\tau _c/2]`$, and $`T(t)0`$ otherwise. The coefficients in the correction terms in the formula (6) are given by the following integrals:

$`I_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \underset{0}{\overset{t}{}}}T(\tau _1\tau _2){\displaystyle \underset{\tau _2}{\overset{\tau _1}{}}}{\displaystyle \underset{\tau _1}{\overset{\tau _2}{}}}T(\tau ^{}\tau ^{\prime \prime })={\displaystyle \frac{1}{24}}t\tau _c+\mathrm{},`$ (7)
$`I_2`$ $`=`$ $`2{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \underset{0}{\overset{\tau _1}{}}}{\displaystyle \underset{0}{\overset{\tau _2}{}}}T(\tau _2\tau ^{})T(\tau _1\tau ^{\prime \prime })={\displaystyle \frac{1}{3}}t\tau _c+\mathrm{},`$ (8)
$`I_3`$ $`=`$ $`{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \underset{0}{\overset{t}{}}}T(\tau _1\tau _2){\displaystyle \underset{0}{\overset{\tau _1}{}}}{\displaystyle \underset{0}{\overset{\tau _2}{}}}T(\tau ^{}\tau ^{\prime \prime })={\displaystyle \frac{t^2}{2}}{\displaystyle \frac{5}{12}}t\tau _c+\mathrm{}.`$ (9)
All the derivatives of $`\kappa (y)`$ in formula (6) are taken at $`y=0`$; the subscripts denote derivatives with respect to the corresponding components of $`𝐲`$, and $`d`$ is the space dimension.

One can easily check that the first-order $`\tau _c`$ correction to the growth rate is negative. It is important to note that this correction is not universal since the integrals (7-9) depend on the chosen regularization $`T(t)`$ of the $`\delta `$ function.

## III Eulerian frame

In the Eulerian frame, the point $`𝐱`$ is fixed in the formula (4), and we need to express $`\stackrel{~}{v}^i\mathbf{(}y(x,t),\tau \mathbf{)}`$ in terms of $`v^i(x,\tau )`$. By analogy with the expression (5), we write $`\stackrel{~}{v}^i(y,\tau )=v^i(x^k\underset{\tau }{\overset{t}{}}\stackrel{~}{v}^k(y,\tau ^{}),\tau )`$. The expansion now takes the form

$`\stackrel{~}{v}^i(\tau )=v^i(\tau )v_l^i(\tau ){\displaystyle \underset{\tau }{\overset{t}{}}}v^l(\tau ^{})+v_l^i(\tau ){\displaystyle \underset{\tau }{\overset{t}{}}}v_j^l(\tau ^{}){\displaystyle \underset{\tau ^{}}{\overset{t}{}}}v^j(\tau ^{\prime \prime })+{\displaystyle \frac{1}{2}}v_{lj}^i(\tau ){\displaystyle \underset{\tau }{\overset{t}{}}}v^l(\tau ^{}){\displaystyle \underset{\tau }{\overset{t}{}}}v^j(\tau ^{\prime \prime })+\mathrm{},`$ (10)
where we use the short-hand notation $`v_l^i(\tau )v^i(x,\tau )/x^l`$. To calculate the second moment of $`B^i`$, we have to use the expression (4):

$`B^i(x,t)=B_0^k\left(\delta _k^i{\displaystyle \underset{0}{\overset{t}{}}}{\displaystyle \frac{\stackrel{~}{v}^i(y,\tau )}{x^k}}\text{d}\tau \right)^1,`$ (11)
that should be expanded up to the fourth order in $`\underset{0}{\overset{t}{}}\frac{\stackrel{~}{v}^i(y,\tau )}{x^k}\text{d}\tau `$. Then, we substitute the expansion (10) for $`\stackrel{~}{v}`$. Finally, assuming isotropic and spatially homogeneous initial distribution, $`B_0^i=0`$, $`B_0^iB_0^j=\frac{1}{d}\delta ^{ij}H_2(0)`$, and making use of (2), we can find the equation for $`H_2(t)`$. We do not present the detailed calculation here since for the incompressible velocity field, the answer coincides with the Lagrangian case (6), as it should.

The appearance of the term $`\kappa _{jjlm}^{ii}\kappa ^{lm}`$ in Eq. (6) shows that for the finite correlation time $`\tau _c`$, the velocity field $`v^i(x,t)`$ cannot be treated as linear in the calculation of the growth rate. This is obvious since a Lagrangian particle is swept by the finite distance $`v^i\tau _c`$ before the velocity field becomes decorrelated. The appearance of $`\kappa (0)`$ (or, equivalently, $`u_{rms}`$) is a manifestation of the absence of the Galilean invariance. For the velocity field with finite correlation time, such invariance is broken as it can be seen from Eq. (2).

## IV Discussion

The Batchelor advection, i.e. the advection by a smooth velocity field is a good model for many physical problems. It is usually a valid approximation when one is interested in the correlators of the advected fields on the scales $`l_{adv}`$ that are much smaller than the viscous scale of the velocity field $`l_v`$.

An additional simplification arises when the correlation time of the velocity field $`\tau _c`$ is much smaller than the inverse growth rates of the advected fields. Mathematically, such a limit is described by the $`\delta `$-correlated in time velocity field. In such a case, the growth rates depend only on $`\kappa ^{\prime \prime }(0)`$, as a consequence of the scale invariance. In this case one can substitute the linear velocity field $`v^i(x,t)=\sigma _k^i(t)x^k`$, $`\sigma _l^i(t)\sigma _m^j(0)=\kappa _{lm}^{ij}(0)\delta (t)`$, for the real field, and thus reduce the field-theoretical problem to the quantum mechanics of the fluctuating matrix $`\sigma _k^i`$(t); for details and references, see .

We have demonstrated that the possibility of such a reduction is an artifact of the $`\delta `$-correlated velocity field, and have presented a simple method for calculating the finite-time corrections to the growth rates. The criterion of the applicability of the quantum mechanical reduction is not only $`l_{adv}/l_v1`$, but also $`\tau _cv/x1`$. For $`\tau _c0`$, the universal growth rates do not exist, they are determined by the statistics of the velocity field on the integral scale and by the form of the short-time regularization. This fact is sometimes overlooked in the literature; in the present note we have tried to clarify the question. The general formalism allowing to systematically find the finite-time corrections to the dynamo growth rates, has been developed in  by a different method. A simple discussion of the stochastic calculus applied to the dynamo problem can be found in .

I am grateful to Alexander Schekochihin for many valuable discussions and comments on both the substance and the style of the paper, and to Michael Chertkov for important conversations.