text:
# Weak Wave Turbulence Scaling Theory for Diffusion and Relative Diffusion in Turbulent Surface Waves

## I Introduction

The study of particles moving on surface waves has shown that the particle motion often is far from being Brownian Ramshankar90 ; Ramshankar91 ; Schroder97a ; Schroder97b ; Hansen97 . Similar results are found in ocean studies Okubo71 . The turbulence observed in surface waves is strongly influenced by the presence of a dispersion relation (for deep-water surface waves) $`\omega _k^2=gk+\sigma k^3/\rho `$ between the wave oscillation frequency and the wave vector $`k`$ dispersion . Here $`g`$ is the gravitational acceleration, $`\sigma `$ is the surface tension coefficient, and $`\rho `$ is the density. Experimentally surface waves can be studied by vertically oscillating a fluid with a free surface. On the surface of the fluid, waves with a clearly discernible wavelength $`\mathrm{\Lambda }`$ are then formed if the amplitude of oscillations exceeds a critical value, the so-called Faraday instability nevolin . The excitation of surface waves in the Faraday system is a parametric effect giving rise to a fundamental wave frequency $`\mathrm{\Omega }`$ equal half the driving frequency. The frequency $`\mathrm{\Omega }`$ and wavelength $`\mathrm{\Lambda }`$ are associated through the dispersion relation. At small frequencies, $`k<k_a=(g\rho /\sigma )^{1/2}`$, the first term in the dispersion relation dominates, and the waves formed are called gravity waves (typical of ocean waves). At large frequencies ($`k>k_a`$) the effect of gravity can be neglected, and the second term dominates. The waves thus formed are called capillary waves. For water, the wave length $`a=2\pi /k_a`$, separating gravity and capillary waves, is about 2 cm.

Several studies have been carried out to measure the diffusivity of particles moving on capillary surface waves generated above the Faraday instability Ramshankar90 ; Ramshankar91 ; Schroder97a ; Schroder97b ; Hansen97 ; Schroder96 . Two particles separated by a distance $`R`$ move more rapidly apart as $`R`$ is increased. The relative diffusivity can be defined as

$$\frac{dR^2}{dt}=2R\delta v,$$
(1)
where $`\delta v`$ is the longitudinal component of the velocity difference. The variance $`\sigma ^2`$ of the relative diffusivity $`dR^2/dt`$ is related to the velocity difference squared,

$`\left({\displaystyle \frac{\sigma }{R}}\right)^2`$ $``$ $`|\delta v(R)|^2.`$ (2)
It has been found experimentally that in capillary waves $`\left(\frac{\sigma }{R}\right)^21a(R/\mathrm{\Lambda })^{1/4}`$, with $`a=0.66`$ Schroder96 .

The purpose of the present investigation is to examine the applicability and limitations of the weak wave turbulence theory in explaining experimental scaling results obtained for the diffusion and relative diffusion of particles moving on turbulent surface waves. Diffusion theories based on random wave fields have likewise been considered RWF , and Herterich and Hasselmann Herterich  particularly considered diffusion by surface gravity waves. For relative diffusion on capillary waves we find the exponent $`1/4`$ as in experiments in the regime $`R>\mathrm{\Lambda }`$. For gravity waves, we find an exponent $`7/61.14`$, a result which is intriguingly close to the value $`1.15`$ found in ocean studies by Okubo Okubo71 . We note that the relative diffusion exponent 4/3 obtained in fully developed turbulence theory previously have been used in the discussion of experimental data Monin75 .

To obtain the diffusivity, the particle displacement $`\mathrm{\Delta }x(t)=x(t_0+t)x(t_0)`$ over a time $`t`$ is measured along an arbitrary axis for many initiation times $`t_0`$. From the resulting distribution the variance $`V(t)=[\mathrm{\Delta }x(t)]^2`$ is found, and the diffusivity $`D`$ can be extracted, $`V(t)=2Dt`$. For Brownian motion the diffusivity $`D`$ is a constant. For deep water capillary waves (we only consider deep-water surface waves) formed at high frequencies, the diffusion has been found experimentally to be anomalous

$$V(t)t^\lambda ,$$
(3)
with $`\lambda 1`$ (super-diffusive) Schroder97a ; Schroder97b ; Hansen97 ; Schroder96 . The exponent is observed to change drastically from a value $`\lambda =1.61.9`$ at length scales below $`\mathrm{\Lambda }`$ to a value $`\lambda =1.01.3`$ at length scales above $`\mathrm{\Lambda }`$. The values of $`\lambda `$ is observed to decrease with increasing drive; at large drive $`\lambda `$ attains the Brownian motion value 1Schroder97a ; Schroder97b . Assuming a connection between spatial and temporal correlations we find for capillary waves $`\lambda =16/9`$ for length scales below $`\mathrm{\Lambda }`$, and $`\lambda =1`$ for length scales above, i.e. values rather close to the experimental values. Ocean studies have been performed using floaters Stommel49 , chemical tracers Okubo71 , and near surface drifters. To what extent these studies relate to surface waves is unclear. Our results do, however, indicate a connection at least for the chemical tracers, where Okubo finds an exponent of 2.34 and from weak wave turbulence theory we find $`8/3`$.

## II Spatial correlations and relative diffusion

Spatial correlations between pairs of particles at $`𝐫`$ and $`𝐫+𝐑`$ are related to the velocity difference squared

$`|\delta v(R)|^2`$ $`=`$ $`|𝐯(𝐫+𝐑)𝐯(𝐫)|^2`$ (4)
$`=`$ $`2[|𝐯(𝐫)|^2𝐯(𝐫+𝐑)𝐯(𝐫)].`$
Following the calculations of Ref. Schroder98 , we have an expression for the space-correlation of the velocity

$`𝒞(R)=𝐯(𝐫)𝐯(𝐫+𝐑)`$ $``$ $`{\displaystyle k\omega _kn_ke^{ikR\mathrm{cos}\theta }𝑑𝐤}`$ (5)
$``$ $`{\displaystyle _0^{\mathrm{}}}k^2\omega _kn_kJ_0(kR)𝑑k,`$
where $`n_k`$ is the isotropic Kolmogorov spectrum Zakharov92 , and $`J_0`$ is the zeroth order Bessel function.

For capillary waves, the dispersion relation is $`\omega _k=(\sigma /\rho )^{1/2}k^\alpha `$ with $`\alpha =3/2`$, and the isotropic Kolmogorov spectrum (three wave interactions) is Zakharov92 ; Zakharov67

$$n_kP^{1/2}\rho ^{3/4}\sigma ^{1/4}k^\beta ,$$
(6)
where $`\beta =17/4`$ and $`P`$ is the energy flux assumed to be constant. From measurements of the spectrum of wave-amplitudes henrya ; henryb  we know that the spectrum levels off below $`K=2\pi /\mathrm{\Lambda }`$. In order to account for this change we introduce a scaling function $`g(k/K)`$ across the drive multiplying $`n_k`$. At $`k<K`$, $`g(k/K)`$ must rise from zero in such a way as to quench the decrease of $`n_k`$. At $`k>K`$, $`g(k/K)`$ can be assumed constant (up to the inverse dissipative scale $`k_d`$ where the theory breaks down). The motion at length scales below the dissipative scale is expected to become ballistic. Also the system size L is a relevant length scale. Obviously the diffusion change character near the system size since particles cannot travel longer distances than L. Also we shall consider the characteristic length scale $`a`$ separating capillary from gravity waves.

The integral for the spatial correlations is given by

$$𝒞(R)k^{2+\alpha \beta }g(k/k_F)J_0(kR)𝑑𝐤$$
(7)
where the specific $`k`$ dependence for $`\omega _k`$ and $`n_k`$ has been inserted. For capillary waves $`2+\alpha \beta =3/4`$. This integral depends explicitly on $`R`$. Consider the value $`K_R2.4/R`$, the lowest $`k`$ value at which $`J_0(kR)=0`$. For $`K_R>K`$, the dominant contribution to the integral comes from the range $`K<k<K_R`$, where the scaling function $`g`$ is constant,

$`𝒞(R)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}k^{2+\alpha \beta }g(k/K)J_0(kR)𝑑k`$ (8)
$``$ $`k^{3+\alpha \beta }|_K^{K_R}R^{3\alpha +\beta }.`$
The result for the velocity difference is

$$|\delta v(R)|^21b(KR)^{3\alpha +\beta }1b(KR)^{1/4},$$
(9)
where $`b`$ is a constant of order one note ).

For $`K_R<K`$, the Bessel function oscillations set in below $`K`$ while the integrand is still increasing in size due to the behavior of the scaling function $`g(k/K)`$. Therefore the dominant contribution to the integral comes from a peak centered around $`K`$, it is oscillatory in nature with a vanishing envelope falling off like $`R^{3/2}`$. For large $`R`$ values we essentially have $`|\delta v(R)|^2`$ constant.

For gravity waves the dispersion relation is $`\omega =g^{1/2}k^\alpha `$ with $`\alpha =1/2`$ (four-wave interaction), and we have energy flow towards higher frequencies compared to $`\mathrm{\Omega }`$ giving rise to the isotropic Kolmogorov spectrum Zakharov92 ; Zakharov66

$$n_kP^{1/3}\rho ^{2/3}k^\beta ,$$
(10)
where $`\beta =4`$ and $`P`$ is the energy flux assumed to be constant. Below $`\mathrm{\Omega }`$ a constant wave-number flux $`Q`$ towards lower frequencies yields Zakharov82

$$n_kQ^{1/3}\rho ^{2/3}g^{1/6}k^{\beta +\alpha /3},$$
(11)
with $`\beta \alpha /3=23/6`$.

As above we invoke a scaling function $`g(k/K)`$ describing the change in the behavior of $`n_k`$ as $`k`$ crosses the wavenumber $`K`$ of the drive. According to the form of the spectra we have

$$g(k/k_F)\{\begin{array}{cc}1\hfill & \text{if }K<k<k_a,k_d\hfill \\ k^{\alpha /3}\hfill & \text{if }k_L<k<K.\hfill \end{array}$$
(12)
Other corrections are relevant at length scales outside the above regime. The motion at sufficiently large length scales ($`k<k_L`$) may be expected to become Brownian, as also found in studies of motion of drifters near the surface of the oceans Okubo71 .

For the velocity difference squared we have

$`|\delta v(R)|^2`$ $``$ $`{\displaystyle _0^{\mathrm{}}}k^{2+\alpha \beta }g(k/K)[1J_0(kR)]𝑑k`$ (13)
$``$ $`{\displaystyle _{K_R}^{\mathrm{}}}k^{3/2}g(k/K)𝑑k.`$
For $`K_R>K`$, the scaling function is constant in the integration region above $`K_R`$, and we have

$$|\delta v(R)|^2K_R^{1/2}(KR)^{1/2}.$$
(14)
For $`K_R<K`$ the part of the integration, where $`g(k/K)k^{1/6}`$, gives the $`R`$ dependent contribution, and

$$|\delta v(R)|^2+constK_R^{1/3}(KR)^{1/3}.$$
(15)
The growth of $`|\delta v(R)|^2`$ stops at length scale $`L`$. We note that the exponent 1/3 for the velocity difference squared implies an exponent 7/6 for the relative diffusivity measure $`\sigma `$. It is interesting that Okubo Okubo71  finds the relative diffusion exponent to be 1.15 for oceanic diffusion, although we should emphasize that the relative diffusion exponent 4/3 obtained in fully developed turbulence theory has generally been applied to explain experimental data Monin75 .

## III Temporal correlations and capillary waves

Next we consider single particle diffusion in capillary waves. For arbitrary times, the mean-square displacement $`r^2=|r(t_0+t)r(t_0)|^2`$ is expressed in terms of the Lagrangian velocity correlation function $`𝒞(\tau )=𝐯(𝐫(t))𝐯(𝐫(t+\tau ))`$ as hansen

$$r^2=2t_0^t(1\tau /t)𝒞(\tau )𝑑\tau .$$
(16)
In the Taylor limit of very small $`t`$, $`𝒞(\tau )v_L^2`$, where $`v_L^2`$ is the mean-square Lagrangian velocity, to be determined by sampling along particle orbits. In this case

$$r^2v_L^2t^2.$$
(17)
For incompressible flows, the mean-square velocity obtained by Eulerian sampling, $`v_E^2`$, is identical to the mean-square velocity, $`v_L^2`$, obtained from Lagrangian measurements Tennekes72 . The subscripts on $`v^2`$ will be omitted from here on. For large times, $`t\tau _L`$, where $`\tau _L=\frac{1}{v^2}_0^{\mathrm{}}𝒞(\tau )𝑑\tau `$ is the Lagrangian integral time scale, Eq. (16) gives a diffusion-like dispersion

$$r^22v^2\tau _Lt.$$
(18)
First we consider capillary waves. If we invoke a scaling form $`r^2t^\lambda `$ at $`t<T`$, then Eq. (16) gives

$$𝒞(\tau )\tau ^{\lambda 2}.$$
(19)
On the other hand, invoking the same scaling $`r^2t^\lambda `$ and using Eq. (8) yields

$$𝒞(\tau )(\tau ^{\frac{\lambda }{2}})^{\frac{1}{4}}.$$
(20)
The result is a self-consistency condition for the exponent $`\lambda `$. We find

$$\lambda =\frac{4}{5+\alpha \beta }=\frac{16}{9}.$$
(21)
This value obtained for $`\lambda `$ is in good agreement with experimental observations which vary $`\pm `$10% from this value Schroder97a ; Schroder97b . From this the exponent for the time correlation function becomes $`2/9`$.

Above it has tacitly been assumed that the time scale $`t`$ is larger than the one corresponding to the dissipative frequency $`\omega _d`$ corresponding to the inverse dissipative scale $`k_d`$. Experimentally, the sampled diffusion distances at small times may approach the dissipative scale at low values of the drive. Below this scale $`\lambda `$ has the ballistic value 2, and the observed slope (1.8–1.9) slightly larger than the theoretical value 16/9 may be an effect hereof.

We have also assumed that the particles can follow the fluid velocity field, i.e. if $`𝐯`$ is the velocity field of the particles and $`𝐮`$ is that of the fluid, our assumption is that $`𝐯=𝐮`$, or at least that the weak ansatz $`𝒞(\tau )=𝐮(t+\tau )𝐮(t)`$ is fulfilled. However, at large driving amplitudes the particles may not follow the fluid velocity. As an example, assume that corrections for the velocity correlation takes the following form

$$𝐯(t)𝐯(t+\tau )\tau ^\delta 𝐮(t)𝐮(t+\tau )$$
(22)
where $`𝐯`$ and $`𝐮`$ are the velocity fields of the particles and of the flow respectively. Following the calculation above, we then find

$$\lambda =\frac{42\delta }{5+\alpha \beta }=\frac{168\delta }{9},$$
(23)
which may explain why a slightly smaller exponent than 16/9 is observed for higher driving amplitudes (e.g. $`\delta =0.2`$ gives an exponent $`\lambda =1.6`$).

Equivalently, the self-consistent condition for $`\lambda `$ (Eq. (21)), can be obtained in terms of the energy spectrum

$$E_kk^\nu ,$$
(24)
where the energy per area is

$$E=E_𝐤𝑑𝐤=E_k𝑑k=k\omega _kn_k𝑑k,$$
(25)
we have the relationship $`\nu =(1+\alpha \beta )`$, and

$$\lambda =\frac{4}{4\nu }.$$
(26)
The capillary wave exponent is $`\nu =7/4`$.

For $`t>T`$, the time correlation $`𝒞(\tau )`$ is oscillatory and therefore negligible for $`t>T`$. In this case, Eq. (18) applies

$$r^2t,$$
(27)
and $`\lambda =1`$. The Brownian motion value obtained for $`\lambda `$ at times $`\tau `$ above $`T`$ is in agreement with experimental observations for large drives. The slightly larger values of $`\lambda `$ at lower drives may arise from corrections to scaling. Experimentally, at sufficiently large times (above those considered in Schroder97a ; Schroder97b ), the sampled diffusion distances are limited by the systems size $`L`$ and eventually one must have $`\lambda =0`$.

## IV Temporal correlations and gravity waves

For gravity waves, we can follow a similar assumption as the one proposed for capillary waves. However, in contrast to the situation for capillary waves, the scaling is given for $`|𝐯(t)𝐯(t+\tau )|^2`$. Thus, for $`\tau <T`$, we have from Eq. (14)

$$|𝐯(t)𝐯(t+\tau )|^2\tau ^{\lambda /4},$$
(28)
and for $`\tau >T`$, we get from Eq. (15)

$$|𝐯(t)𝐯(t+\tau )|^2\tau ^{\lambda /6}.$$
(29)
Assuming

$$\sigma \frac{R^2}{t},$$
(30)
the exponents can now be calculated using Eq. (2) from the self-consistent condition

$$t^{\lambda 1}\{\begin{array}{cc}t^{5\lambda /8}\hfill & \text{small }t\hfill \\ t^{7\lambda /12}\hfill & \text{large }t\hfill \end{array}$$
(31)
The values obtained for $`\lambda `$ are $`8/3`$ for small $`t`$ and $`12/5`$ for large $`t`$. Interestingly, a similar situation is observed in ocean studies where $`\lambda =2.34`$ is found in several dye studies collected and analyzed by Okubo Okubo71 .

## V Conclusions

Good agreement is found between theoretical and experimental scaling results for the diffusion and relative diffusion of particles on weakly turbulent surface waves. In the case of relative diffusion, for capillary waves the exponent -1/4 for the space correlation function is recovered. For gravity waves, a scaling behavior is derived for the velocity difference squared with the exponent changing from 1/2 at small distances to 1/3 at large distances. This corresponds respectively to scaling exponents 5/4 and 7/6 for the relative diffusivity. Intriguingly, Okubo finds a relative diffusion exponent of 1.15 for oceanic diffusion Okubo71 . This differs from the relative diffusion exponent 4/3 obtained in fully developed turbulence theory, previously used in the discussion of experimental data Monin75 .

For capillary waves the theory yields a crossover in the diffusivity at the wavelength $`\mathrm{\Lambda }`$ with a change in diffusion exponent $`\lambda `$ from $`\lambda =16/91.78`$ to $`\lambda =1`$, in good agreement with experimental findings. Particles suspended on the fluid surface may not follow the fluid far above the Faraday instability, and we suggest that this may explain the small discrepancy between weak turbulence theory and experimental results in this regime. For gravity waves the weak turbulence theory gives a diffusion exponent of $`12/5=2.4`$ for long times which compares favoritely with the exponent $`2.34`$ found by Okubo for dye diffusion. To what extent the weak wave turbulence theory really explains the oceanic diffusivity is unclear but it is fascinating to speculate that wave motion taking place mostly on a short scale could be responsible for the long time diffusion in oceans. At very long time scales, i.e. above the Lagrangian time scale, ordinary Brownian motion is expected.

## VI Acknowledgments

We acknowledge financial support from the Danish Natural Sciences Research Council and CORE.


format: text

text:
# Scaling in a Nonconservative Earthquake Model of Self-Organised Criticality

## I Introduction

Many dynamical phenomena in nature are intermittent. This “bursty” dynamics may be related to an underlying complex state, often characterised by long range correlations in space and time. For example, the crust of the earth alternates long periods of relative quiescence with burst of activity (earthquakes), which have a wide range of possible sizes. The behaviour of earthquakes is described by the empirical Gutenberg-Richter (GR) law , where the distribution of energy dissipated in earthquake events is a power law over many orders of magnitude in energy. The GR scaling extends from the smallest measurable earthquakes, which are equivalent to a truck passing by, to the most disastrous that have been recorded. Similar scale free behaviour of bursts is observed in vastly different kinds of physical systems such as flux motion through disordered type II superconductors placed in a magnetic field , or in granular piles, under some conditions , etc.

Self-organised criticality (SOC)  has been proposed as a general dynamical principal behind the observed complex behaviour of many natural phenomena. It refers to the fundamental property of slowly driven, extended systems to organise, over a sufficiently long transient period, into a dynamical critical state which lacks any characteristic time or length scale. The amplitude of the response of the system to an external perturbation follows a power law distribution. A number of simple models have been developed to test the applicability of SOC to a variety of complex interacting dynamical systems, such as sand piles and earthquakes (for a review see e.g. ref. ).

One of the basic theoretical problems is to identify robust, and thus physically relevant mechanisms for SOC to emerge, and to define a space of parameters and dynamical processes where SOC is a stable feature. Although it has been proposed that the presence of conservation laws (e.g. sand grains being transported in a sand pile) or special symmetries was necessary for SOC , many known examples of physical phenomena and some models have been found where no apparent conservation law or special symmetry exists. These include, besides earthquakes, biological evolution, forest fires, epidemics, (possibly) solar flares, (possibly) reconnection events in the magnetotail, etc. . In contrast to conservative systems, the mechanisms for SOC in non-conservative systems are not very well established.

A nonconservative SOC model that in recent years has attracted much attention is the so called OFC model . The OFC model is a simplified lattice representation of a spring-block model for earthquake dynamics which was originally introduced by Burridge and Knopoff . The Burridge-Knopoff model can be schematized as a two dimensional network of blocks interconnected by springs. All blocks are subject to an external driving force, which pulls them, and to a static friction, which opposes their motion. In the OFC model, each site of a lattice is associated with a continuous variable, which represent the force acting on a block. A slow driving is applied to the system by increasing uniformly and simultaneously the forces of all the elements. When the force on a site exceeds some threshold value (the maximal static friction), the site relaxes and distributes part of its force to nearest neighbour sites. Each such discharge event is accompanied by a local loss in accumulated force from the system. This conceptually simple and seemingly numerically tractable model reproduces some of the qualitative phenomenology of the statistics of earthquake events such as power law behaviour over a range of sizes, and intermittency or clustering of large events . Extensions of the model have been recently developed which reproduce to some extent Omori’s law and other temporal patterns associated with earthquakes .

In the context of nonconservative models, the OFC model is of particular interest as it is possible to directly control the level of conservation of the dynamics through a parameter $`\alpha `$. Early analysis on relatively small systems indicated that the OFC model exhibited SOC, in the sense that earthquakes in the steady state obeyed finite size scaling (FSS) when the system size was varied . However, the critical coefficients obtained using the FSS ansatz were found to be nonuniversal. In particular the exponents characterising the power law distributions appeared to vary with both the dissipation parameter, $`\alpha `$, and the form of the boundary conditions. This would have been in sharp contrast to the usual fixed point picture of critical phenomena where most microscopic details are irrelevant. Moreover some apparent critical exponents obtained using FSS violated physical bounds , putting some doubt on the existence of criticality in the model. Recently it was shown using a multiscaling analysis of large-scale simulations that, actually, the avalanche size distribution has a universal power law behaviour, independent of the dissipation parameter and for different boundary conditions, but that the cutoff in the power-law distribution does not behave according to FSS . In larger systems, proportionally more of the force can be dissipated in the largest events that occur, and the cutoff function becomes sharper and sharper as the system size increases.

Departures from standard FSS have been reported for various SOC models as, for example, the sandpile model , the Drossel-Schwabl (DS) forest fire model  and the Zhang model . In this paper we address the question of the origin of the breaking of FSS and its relation to the mechanism responsible for SOC in the OFC model. In particular we test the implicit assumption behind the FSS hypothesis that a finite systems behaves as a subsystem of a larger system. The paper is organised as follows. In the section II we describe in some detail the model. In section III we present the results of our numerical study relative to two different type of probability distributions for earthquake sizes. The first distribution concerns earthquakes which are localised within a given subsystem. We show that this subset of earthquakes exhibits ordinary FSS as long as the linear extent, $`\lambda `$, of the subsystem is sufficiently small compared to the linear extent, $`L`$, of the entire system. The second distribution groups earthquakes according to the position of their starting site relative to the boundaries of the system. From this investigation, we deduce that FSS is violated due to large events initiated in a region near the boundary. Finally, in section IV we discuss our results and draw some conclusions.

## II The model

We consider a two-dimensional square lattice of $`L\times L`$ sites. To each site $`i`$ of the lattice we associate a continuous variable $`F_i`$, which initially take some random values between zero and a threshold value $`F_{th}`$. The dynamics proceeds then indefinitely. In the limit of infinite time scale separation between the slow driving and the (almost) instantaneous earthquake process, the dynamics is:

1. Uniform drive: all forces $`F_i`$ are increased at the same rate, until one of them reaches the value $`F_{th}`$.
2. Earthquake: when a site becomes unstable (i.e $`F_iF_{th}`$), the uniform driving is stopped and the system evolves according to the following local relaxation rule

$$F_iF_{th}\{\begin{array}{c}F_i0\hfill \\ F_{nn}F_{nn}+\alpha F_i\hfill \end{array}$$
(1)
until there are no more unstable sites. In eq. (1), the subscript “nn” stands for the four nearest neighbours to site $`i`$.

Since only a fraction, $`4\alpha `$, of the force is redistributed in each relaxation event (toppling), the model is nonconservative for $`\alpha <1/4`$. In the following we concentrate on this case, i.e. $`0<\alpha <1/4`$.

To completely define the model we need to specify the boundary conditions. Boundaries are believed to play a crucial role for the observation of critical behaviour in the OFC model. It has been suggested that they act as inhomogeneities which frustrate the natural tendency of the model to order into a periodic state . Indeed, for sufficiently small values of the conservation parameter $`\alpha `$ ($`\alpha <\alpha _c0.18`$), a system with periodic boundary conditions quickly reaches an exactly periodic state with only earthquakes of size one. For larger values of $`\alpha `$ the situation is slightly more complicated with multiple topplings involved in a single avalanche, but the avalanches are still localised and criticality is not observed . A system with open boundaries is prevented from reaching a periodic state because boundary sites have fewer neighbours and therefore cycle at a different frequency from bulk sites. Middleton and Tang suggested that the inhomogeneity created by the boundaries propagates into the bulk of the system, developing, in this way, long range spatial correlations. They named this mechanism “marginal synchronisation” or phase locking . In accordance with previous studies, therefore, we consider open boundary conditions. If a boundary (corner) site topples, an extra amount $`\alpha F_i`$ ($`2\alpha F_i`$) is simply lost by the system.

## III Results: probability distributions

After a sufficiently long transient time, the system reaches a stationary state. Several statistical properties can be used to characterise this state. Most previous studies of the OFC model have focused on the behaviour of the probability distribution of earthquake sizes, $`P_L(s)`$, where $`L`$ is the size of the system and $`s`$ is the total number of topplings events during an earthquake . We choose instead to analyse the behaviour of different distributions for avalanches sizes, which distinguish between earthquakes according to the region of the lattice involved (e.g. bulk or boundary) and the coordinates of the triggering site (see figure 1). This investigation is particularly relevant for the OFC model in view of the strong inhomogeneity in the spatial distribution of avalanches . According to ref. , for example, large avalanches are localised near the boundary (at least for $`\alpha <\alpha _c`$). As a minor technical remark, we point out that we exclude from our data avalanches which involve only one site ($`s=1`$) as they appear to behave according to their own statistics . As we are mainly interested in asymptotically large earthquakes, this does not alter our conclusions.

We consider a subsystem of linear extent $`\lambda `$ centred in a system of size $`L`$. The first distribution we introduce, $`P_{conf}(\lambda ,L,s)`$, is the normalised distribution of earthquake sizes restricted to earthquakes which are confined entirely within the subsystem (e.g. avalanche (a) in fig. 1). The model is driven according to its usual dynamics but only those particular earthquakes are counted. According to our definition, the case $`\lambda =L`$ corresponds to the distribution of avalanches which do not reach the boundary of the system. As shown in fig. 2, the distribution $`P_{conf}(\lambda ,L,s)`$ becomes independent of $`L`$, if $`L`$ is considerably larger then $`\lambda `$ (approximately $`L2\lambda `$). When $`L`$ approaches $`\lambda `$, this is no longer the case and the cutoff in the distribution is pushed to larger sizes. Although we have shown in Fig. 2 only the distributions for $`\alpha =0.18`$ and $`\lambda =32,64`$, analogous considerations apply to different values of $`\alpha `$ and for different sizes, $`\lambda `$. Since for a generic $`L`$, $`P_{conf}(\lambda ,L,s)P_{conf}(\lambda ,\lambda ,s)`$, a small portion of a large system is substantially different from a finite system of the same size, contrary to what happens in equilibrium critical phenomena. A similar observation was made in ref.  for the DS forest fire model. In the following, we denote with $`P_{conf}(\lambda ,s)`$ the distribution $`P_{conf}(\lambda ,L,s)`$ in the limit where the distribution does not appear to depend on $`L`$. In order to determine numerically these distributions, for each value of $`\lambda `$ we have simulated (for at least $`210^9`$ earthquakes) a system of size $`L=2\lambda `$. The dependence on $`L`$ of $`P_{conf}`$ can in this case be safely neglected. With this choice, the accuracy of the measures and the range of scales investigated are optimised, within our computational limits.

In Fig. 3 we report a FSS collapse of $`P_{conf}(\lambda ,s)`$ for different values of $`\alpha `$ . Contrary to the entire distribution of earthquake sizes, $`P_L(s)`$, we observe that $`P_{conf}(\lambda ,s)`$ satisfies the FSS hypothesis, i.e. $`P_{conf}(\lambda ,s)\lambda ^\beta f(s/\lambda ^D)`$, with universal critical coefficients. The curve corresponding to $`\alpha =0.15`$ and $`\lambda =256`$ shows some noisy behaviour, due to the difficulties in collecting good statistics in this case. Indeed by decreasing $`\alpha `$, the relative fraction of earthquakes in the bulk of the system (with size $`s>1`$) diminishes. Nonetheless, there is no evident sign that FSS is violated in this case. The critical exponents used in the fit of Fig. 3 are $`\beta =3.6`$ and $`D=2`$, independent of the dissipation parameter $`\alpha `$. The value of the histogram exponent $`\tau =\beta /D1.8`$ we obtain is the same as that found for $`P_L(s)`$ . In addition, the value of $`D`$ we find corresponds to the largest possible value for the entire distribution ($`D_{max}`$ in ref. ), as it can be shown that non-conservation requires $`D2`$ .

The scaling behaviour of $`P_{conf}`$ appear to be reasonably robust with respect to translating the subsystem within the entire system; in fig. 4 we report a FSS plot for the subsystem placed on a boundary and on a corner of the system for $`\alpha =0.18`$. While the FSS collapse for the subsystem placed on the boundary is rather good, some deviations from FSS are observed in the cut-off region for the case of the subsystem placed in the corner. We believe this behaviour can be ascribed to the enhanced boundary effects in the latter case (two side of the subsystem are boundary sides instead of only one) and would disappear if larger (sub)systems could be studied. This picture is confirmed by choosing different $`\alpha `$ values: for the subsystem in the corner, deviations from FSS are more pronounced for $`\alpha =0.21`$ and are absent for $`\alpha =0.15`$. For the subsystem on the boundary, instead, the quality of FSS collapse is rather convincing in all cases.

We introduce next the distributions $`P_<(\lambda ,L,s)`$ and $`P_>(\lambda ,L,s)`$. These are the normalised distribution of earthquakes which start respectively within ($`P_<`$) and outside ($`P_>`$) the subsystem of size $`\lambda `$, irrespective of whether they stay in or go out of the subsystem (see fig. 1). The only difference between these two distributions, $`P_<`$ and $`P_>`$, is the location of the site that triggers the avalanche. We observe numerically that the distributions $`P_<(\lambda ,L,s)`$ and $`P_>(\lambda ,L,s)`$ become independent of $`\lambda `$ respectively in the limit $`\lambda L`$ and $`\lambda L`$. As an example we report in Fig. 5 the behaviour of $`P_<(\lambda ,L,s)`$ and $`P_>(\lambda ,L,s)`$ for $`\alpha =0.18`$, $`L=256`$ and for various $`\lambda `$. For simplicity, in the following we denote with $`P_<(L,s)`$ and $`P_>(L,s)`$ the distributions in the limit where they do not depend on $`\lambda `$.

We consider therefore two centred subsystems of linear extent $`\lambda _2>\lambda _1`$, such that the above conditions are satisfied. More specifically, we choose $`\lambda _1=\frac{3}{16}L`$ and $`\lambda _2=\frac{7}{8}L`$. In this case $`P_>`$ corresponds to the subset of earthquakes which are initiated in some “boundary” region and $`P_<`$ corresponds to the subset of earthquakes which are initiated within some “bulk” region. In fig. 6 we report a FSS scaling plot both for $`P_>`$ and $`P_<`$. In this figure we choose $`D=2`$ as the maximum allowed value. It is clear that the boundary distribution $`P_>`$ cannot be collapsed according to the FSS ansatz. In fact it develops a sharper and sharper cutoff which changes shape and which has an excess of large events (the cutoff moves towards right for increasing $`L`$). The bulk distribution $`P_<`$ instead does not develop a noticeably sharper cutoff and does not appear to change its shape. It may possibly be collapsed according to the FSS ansatz. Consistent with the results for $`P_{conf}`$ and for $`P_L`$, the power law exponent for $`P_<`$ is $`\tau =\beta _1/D1.8`$.

The above numerical analysis indicates that the large events which violate FSS are triggered by sites in a boundary region. Indeed the behaviour of the cutoff for the collapsed probability distributions $`P_>`$ and $`P_L`$ is very similar (see Fig. 1 in ref. ). Although large earthquakes are focussed mainly toward the boundary, as suggested in ref. , they occur also in the bulk of the system, even for low values of $`\alpha `$, as can be deduced from Fig. 3 and Fig. 6. Moreover, we do not observe any significant qualitative change in the behaviour of the system around $`\alpha =\alpha _c0.18`$, as claimed in ref. .

## IV Discussion and conclusions

Similarly to other SOC systems , the nonconservative OFC model shows relevant deviations from simple FSS . In this paper we have investigated the origin of this phenomenon, finding that FSS in the OFC model is violated because of large, system wide earthquakes. In fact, we have found that earthquakes localised within a given subsystem do obey ordinary FSS, with universal critical exponents, independently of whether the subsystem is placed at the centre or on a boundary of the system. The value of the power law exponent, $`\tau 1.8`$, for the “confined” distribution agrees with the one for the entire distribution. We have shown, moreover, that the probability distribution for earthquakes initiated in a boundary region do not obey FSS, because of an “excess” of large events. This would result in an apparent exponent $`D>2`$ which is not allowed in the nonconservative case. On the other hand, the probability distribution for earthquakes starting in the bulk of the system is compatible with a FSS hypothesis. In particular, the critical exponent is $`D=2`$ in this case, indicating that large earthquakes responsible for breaking finite size scaling are initiated predominantly near the boundary of the system.

Self-organised criticality in the OFC model has been ascribed to a mechanism of “marginal synchronisation” . A system with open boundaries becomes almost synchronised by an invasion process where spatial correlations develop from the boundaries. It was suggested that sites close to the boundaries start to organise themselves first, building up long range correlations. The critical region grows with time, until, in the stationary state, it invades the whole lattice . Our findings on the large events occurring at the boundary seem to indicate that the effect of synchronisation is stronger for boundary sites than for bulk sites. This view is supported also by the “on screen” observation that large earthquakes tend to be triggered repetitively by the same sites over a long time scale (a result which seems to be confirmed also by the study in ref. ).

S.L. acknowledges financial support from EPSRC (UK), Grant No. GR/M10823/01.


text:
# Addendum to “COULOMB INTERACTION DOES NOT SPREAD INSTANTANEOUSLY”

## Abstract

We describe the results of experimental research by configurations changed with respect to the original experimental configuration . In this way it is demonstrated that the basic part of the signal registered is due to the Coulomb electric field and does not originate from the transversal electromagnetic wave.

Received 2001

PACS: 03.50.-z, 03.50.De

The basic doubt regarding the reliability of our results is due to the fact that it is not known which part of the signal obtained by us is due to the Coulomb interaction, and which part is due to the presence of a regular transversal electromagnetic wave. The theoretical consideration is difficult because of the complexity of the experimental construction.

For this reason, it is more appropriate to determine this by means of experiment.

In our Coulomb electric field generator (Fig. 1) there are three elements which emit transversal electromagnetic waves:

FIG. 1. Regular experimental configuration.

> 1. the currents which flow along the surface of the spheres of the Van de Graaf generators
>
> 2. the spark, which jumps between the electrodes of the discharger
>
> 3. the currents which flow along the discharge cable

During our experiments, the discharge cable was electrically and magnetically shielded. For its construction, the discharge cable was heavy coaxial cable. The discharge current flowed along the central conductor. As an electric shield, its shield with a thickness of approximately 1 mm was used, which is much greater than the depth of the skin layer (0.03 mm) for the resonant frequency of 14 Mhz. For the magnetic shield, a 5 mm thick iron sectioned tube was used. The coaxial cable was placed inside the iron tube without electric connection between them. The discharger was not shielded. The shielding of the spheres of the Van de Graaf generators is not recommended as the signal due to the Coulomb electric field would be lost.

Consequently, it is necessary to evaluate the contribution of the transversal electromagnetic waves emitted both by the discharger and by the Van de Graaf generator’s spheres. In practice, it is only necessary to consider the influence of the active Van de Graaf generator as the passive generator is considerably further away from the antenna 2. For this reason, the signal originating from the passive Van de Graaf generator reaches antenna 2 (the measuring antenna) with a delay greater than 60 ns. That is, the said signal reaches antenna 2 after measuring and does not influence the experimental result.

The general idea of the proposed experiments is simple. If the initial direction of the currents that flow through the discharger and along the sphere’s surface changes, then the sign of the first front of the electromagnetic wave emitted by them must change. Consequently, if the predominant part of the signal is due to the electromagnetic wave, the sign of the first flank of the signal obtained by the antennas will change. The amplitude will change simultaneously.

With this purpose, an orifice will be perforated in the head of the active Van de Graaf generator. A thick ceramic wall tube, which does not allow electric rupture up to a voltage of 150 kV (Fig.2), is placed in the orifice.

FIG. 2. Modified head of the Van de Graaf active generator.

A conductor, which is connected electrically to the other pole of the sphere, is placed along the axis of the tube. As we will demonstrate further on, a construction of this type allows for a reversal of the currents along the surface of the sphere. It is important to point out that the electromagnetic wave emitted from the current, which passes through the central conductor, is isolated from the metallic walls of the sphere.

Experiments have been carried out with four configurations. In fig.1 we show a regular configuration; in figures 3, 4 and 5, the configurations are modified. The direction of the initial movement of the electrons is indicated by arrows in all of the figures.

FIG. 3. Inverse direction configuration of the electric current which passes along the surface of the head of the Van de Graaf active generator.

FIG. 4. Inverse direction configuration of the electric current which passes through the discharger.

FIG. 5. Inverse direction configuration of the electric current which passe through the discharger and along the surface of the head of the Van de Graaf active generator.

Each configuration has been investigated in three regimes:

> a. the discharge cable is not shielded.
>
> b. the discharge cable is a coaxial cable. The discharge current passes through the central thread. The cable shield is not connected to earth.
>
> c. as in b., but the shield is connected to earth.

In this way we can evaluate the influence of the discharge cable screen on the processes. As we can see from the figures, the general flow of the electrons is always in the same direction - from the passive Van de Graaf generator towards the active Van de Graaf generator. Because of this, the signal due to the Coulomb electric field in the case of the four configurations does not undergo significant changes. Consequently, if the basic part of the signal is due to the Coulomb interaction, we should observe, in the case of the four configurations, an impulse with negative first flank (Fig.6).

FIG. 6. Real form of the signal obtained by all of the configurations.

At the same time, in two configurations the direction of the current in the discharger changes in the opposite way (Fig.4, Fig.5). In an analogous manner, the direction of the currents which flow along the surface of the sphere changes in two configurations (Fig.3, Fig.5). Consequently, if our signal is due to the electromagnetic waves emitted, or to the discharger, or the sphere, or both together, at least in one configuration we should observe a change in the sign of the first flank on the oscilloscope screen (Fig. 7).

FIG. 7. Supposed form of the signal generated from the transversal electromagnetic wave by changed configuration.

As the experiments are carried out we can see that ALWAYS, regardless of the configuration and measurement regime, an impulse as pointed out in Fig.6 is obtained.

The proportion between the signal amplitudes, obtained in the configurations presented in figures 3, 4 and 5, taking the signal amplitude of the regular configuration as reference point (Fig. 1), has been investigated. The result is that the relative deviation of the amplitudes does not exceed 5%.

Consequently, we can arrive at the definite conclusion that the dominant part of our signal is due to the Coulomb interaction. Consequently, if the Coulomb interaction is propagated instantaneously, we should record this fact. However, in reality we only record the propagation of a wave. In this manner, we can consider that it has been demonstrated experimentally that the Coulomb electric field is propagated like a Coulomb wave (i.e. not instantaneously).


text:
# Momentum transferred to a trapped Bose-Einstein condensate by stimulated light scattering

## I Introduction

Recent experiments based on two-photon Bragg spectroscopy  have given first valuable information on the density-density response function of a trapped Bose-Einstein condensate, thereby providing an important opportunity to understand the role of quantum correlations and to check theoretical predictions on the dynamics of nonuniform condensates at both low and high excitation energy . In these experiments two laser beams impinge upon a condensate, whose atoms can undergo stimulated light scattering events by absorbing a photon from one of the beams and emitting into the other. The difference in the wave vectors of the beams defines the momentum transfer in a single scattering event, $`\mathrm{}𝐪`$ , while the frequency difference defines the corresponding energy transfer, $`\mathrm{}\omega `$. Both the values of $`𝐪`$ and $`\omega `$ can be tuned by changing the angle between the two beams and varying their frequency difference. For small values of $`𝐪`$ the system is excited in the phonon regime and the response is usefully detected by measuring the net momentum imparted to the gas. The latter is measured by observing the center of mass motion of the sample after releasing the trap.

Since the momentum acquired by the condensate is the physical quantity measured in these experiments, it is important to understand what are the underlying mechanisms characterizing its behavior and what can be learnt from a careful study of its time dependence. In this work we investigate this problem by using the Gross-Pitaevskii (GP) theory for the order parameter of the condensate, with special emphasis on the phonon (low $`q`$) regime. First we present results of simulations based on the numerical integration of the time dependent GP equation. Then we analyze the same results from the viewpoint of linear response theory, testing the accuracy of the local density approximation for the dynamic structure factor. Finally we compare our results with the experimental data obtained from two-photon Bragg spectroscopy .

## II Time dependent Gross-Pitaevskii equation

The interaction of the atoms with the laser field is described by the Hamiltonian

$$H_{\mathrm{Bragg}}=\frac{V}{2}\left(\delta \rho _𝐪^{}e^{i\omega t}+\delta \rho _𝐪e^{i\omega t}\right)$$
(1)
where $`\delta \rho _𝐪=\rho _𝐪\rho _𝐪`$ is the fluctuation of the density operator

$$\rho _𝐪=\underset{j=1}{\overset{N}{}}e^{i𝐪r_j},$$
(2)
while the strength of the perturbation, $`V`$, is proportional to the intensity of the laser beam and can be related to the two-photon Rabi frequency . The validity of expression (1) is not limited to the regime of small perturbations, but can also be used to explore non linear effects, including, for example, the coherent splitting of the condensate .

We simulate the time evolution of the system by solving the Gross-Pitaevskii (GP) equation

$$i\mathrm{}\frac{}{t}\mathrm{\Phi }(𝐫,t)=\left[\frac{\mathrm{}^2}{2m}^2+V_{\mathrm{ho}}(𝐫)+g|\mathrm{\Phi }|^2+\vartheta (t)V\mathrm{cos}(qz\omega t)\right]\mathrm{\Phi }(𝐫,t)$$
(3)
for the order parameter of the gas, $`\mathrm{\Phi }(𝐫,t)`$. The harmonic confinement is represented by the potential $`V_{\mathrm{ho}}(𝐫)=(m/2)[\omega _{}^2r_{}^2+\omega _z^2z^2]`$, with $`r_{}^2=x^2+y^2`$, while the coupling constant in the mean field term, $`g=4\pi \mathrm{}^2a/m`$, is fixed by the $`s`$-wave scattering length $`a`$. The last term in the r.h.s. comes from Eq. (1), while $`\vartheta (t)`$ is the Heaviside step function. The GP equation (3) is expected to provide an accurate description of the system at sufficiently low temperature . We numerically solve it for a system of sodium atoms confined by a cigar-shaped potential, with $`\omega _{}=8\omega _z`$. For the Thomas-Fermi parameter $`15Na/a_{}`$ we choose the value $`10^4`$. Here $`a_{}=\sqrt{\mathrm{}/m\omega _{}}`$ is the oscillator length in the transverse direction. Actually, by a proper rescaling of the spatial coordinates (in units of $`a_{}`$) and energies (in units of $`\mathrm{}\omega _{}`$), the solutions of the time-dependent GP equation (3) depend only on the Thomas-Fermi parameter, the ratio $`\omega _z/\omega _{}`$ and, of course, the parameters characterizing the Bragg potential. Finally, we take $`𝐪`$ along $`z`$ to preserve the axial symmetry of the problem as in the MIT experiment . In order to select the phonon regime, one has to choose a wave vector $`q`$ such that $`q\xi <1`$, where $`\xi `$ is the healing length. The latter is defined as $`\xi =[8\pi na]^{1/2}`$, where $`n`$ is the density evaluated in the center of the trap. In our calculation, a typical value is $`q=3a_{}^1`$, with $`q\xi =0.42`$.

The numerical integration of the time dependent GP equation (3) is performed by using a method previously developed in , which is suitable for axially symmetric condensates. The equation is mapped on a two-dimensional grid of points and solved by using a combination of Fast Fourier and Crank-Nicholson integration algorithms. From the solution of Eq. (3), one can easily calculate net momentum imparted to the condensate by direct integration of the current density associated with the order parameter

$$P_z(t)=\frac{\mathrm{}}{2i}𝑑𝐫\mathrm{\Phi }^{}(𝐫,t)_z\mathrm{\Phi }(𝐫,t)+\mathrm{c}.\mathrm{c}..$$
(4)
In Fig. 1 we show the results obtained for $`P_z(t)`$, having fixed the frequency difference of the Bragg pulses to be $`\omega =13\omega _{}`$. We use different values of the strength $`V`$ in order to check whether we are in the appropriate limit of linear response to the external perturbation. In the linear regime the momentum $`P_z`$ of the condensate should depend quadratically on the intensity $`V`$ of the Bragg pulse. We have verified that, by taking $`V`$ smaller than $`0.2\mathrm{}\omega _{}`$, the linear regime is ensured for times shorter than $`0.2(2\pi /\omega _{})`$. We have solved the time dependent GP equation also in the case of a strong Bragg pulse ($`V>10\mathrm{}\omega _{}`$). For such a pulse the change in the momentum distribution of the condensate is huge  and our simulation shows that, in this case, the momentum transfer stays systematically below the value predicted by linear theory.

In Fig. 2 the net momentum transfer is shown as a function of the frequency $`\omega `$ for two different choices of the time duration $`\tau _B`$ of the Bragg pulse. The figure shows that the shape of the signal exhibits a significant time dependence, becoming narrower by increasing $`\tau _B`$.

## III Linear response theory

In order to better understand the behavior of the momentum transfer as a function of $`t`$ and $`\omega `$, we use the formalism of linear response function, valid for weak intensities of the Bragg pulse.

First, we write the Heisenberg equation of motion for the momentum $`P_z=_{j=1}^Np_j^z`$:

$$\frac{dP_z(t)}{dt}=\frac{1}{i\mathrm{}}[\underset{j=1}{\overset{N}{}}p_j^z,H_{\mathrm{tot}}]$$
(5)
where the Hamiltonian of the system is

$$H_{\mathrm{tot}}=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2m}+\underset{i=1}{\overset{N}{}}V_{\mathrm{ext}}(𝐫_i)+g\underset{i<j}{}\delta (𝐫_i𝐫_j)+\vartheta (t)H_{\mathrm{Bragg}}.$$
(6)
The commutator can be explicitly evaluated, obtaining the exact equation

$$\frac{dP_z(t)}{dt}=m\omega _z^2Z\frac{iqV}{2}\left(\delta \rho _𝐪^{}e^{i\omega t}\delta \rho _𝐪e^{i\omega t}\right).$$
(7)
The first term on the r.h.s originates from the confining oscillator potential, while the second one from the Bragg perturbation. The quantity $`Z=_{j=1}^Nz_j`$ is the expectation value of the $`z`$-th component of the center of mass coordinate. Notice that the commutator in Eq. (5) is not affected by the kinetic energy term, nor by two-body interactions which are translational invariant quantities.

Now, one can use the linear response theory for the density fluctuations. A convenient way consists in writing the time dependent perturbation as

$$\vartheta (t)H_{\mathrm{Bragg}}=\frac{iV}{4\pi }\delta \rho _𝐪^{}𝑑\omega ^{}\frac{e^{i\omega ^{}t}}{\omega ^{}\omega +i\eta }+H.c.,$$
(8)
which follows from the Fourier representation of the step function. This perturbation induces density fluctuations in the form

$$\delta \rho _𝐪=\frac{Ve^{i\omega t}}{2\mathrm{}}𝑑\omega ^{}\left[S(𝐪,\omega ^{})S(𝐪,\omega ^{})\right]\frac{e^{i(\omega \omega ^{})t}1}{\omega \omega ^{}},$$
(9)
where

$$S(𝐪,\omega )=\frac{1}{𝒵}\underset{m,n}{}e^{\beta E_m}m\delta \rho _𝐪n^2\delta (\omega \omega _{mn})$$
(10)
is the dynamic structure factor of the system, with $`\beta =1/kT`$, $`\mathrm{}\omega _{mn}=E_mE_n`$ and $`𝒵`$ is the usual partition function. In deriving this result, one assumes that the fluctuation $`\delta \rho _𝐪`$ induced by the $`e^{i\omega ^{}t}`$ component of the perturbation (8) oscillates like $`e^{i\omega ^{}t}`$, so neglecting the term oscillating like $`e^{i\omega ^{}t}`$. The latter exactly vanishes in a uniform body and, in general, is exponentially small if the value of $`q`$ is much larger than the inverse of the size of the system. In the following we will always neglect this contribution.

Inserting the density fluctuation (9) into Eq. (7), one gets the following equation for the rate of momentum transfer at $`t>0`$:

$$\frac{dP_z(t)}{dt}=m\omega _z^2Z+\frac{2q}{\mathrm{}}\left(\frac{V}{2}\right)^2𝑑\omega ^{}\left[S(𝐪,\omega ^{})S(𝐪,\omega ^{})\right]\frac{\mathrm{sin}\left[(\omega \omega ^{})t\right]}{\omega \omega ^{}}.$$
(11)
This equation explicitly contains the effect of the confining potential through the first term in the r.h.s.; this term is small for times short compared to the oscillator period $`2\pi /\omega _z`$. We have verified that for the times considered in our simulation the effect of the external potential in the equation for the momentum rate can be safely ignored.

For very short times the expansion of (11) yields the model independent result

$$\frac{dP_z(t)}{dt}=\frac{\omega V^2q^3t^3}{3m}$$
(12)
where we have made use of both the $`f`$-sum rule $`𝑑\omega \omega S(𝐪,\omega )=N\mathrm{}q^2/2m`$ and the relationships

$$𝑑\omega \left[S(𝐪,\omega )S(𝐪,\omega )\right]=𝑑\omega \omega ^2\left[S(𝐪,\omega )S(𝐪,\omega )\right]=0.$$
(13)
A closed equation for $`P_z(t)`$ can be obtained by taking the time derivative of (11) and using the exact equation $`dZ(t)/dt=P_z(t)/m`$. If the duration of the pulse is short compared to the oscillator time, but at the same time large compared to the inverse of the frequency of the applied field, then the equation for the momentum rate approaches the golden rule result:

$$\frac{dP_z(t)}{dt}=q\left(\frac{V}{2}\right)^2\frac{2\pi }{\mathrm{}}\left[S(𝐪,\omega )S(𝐪,\omega )\right].$$
(14)
Although restrictive, the conditions $`\omega t1`$ and $`\omega _zt1`$, needed to derive result (14), are compatible. In principle, even in the phonon regime where the excitation energy $`\mathrm{}\omega `$ should be smaller than the chemical potential $`\mu `$, the two conditions can be simultaneously satisfied if one works in the Thomas-Fermi regime where $`\mu \mathrm{}\omega _{}`$. In practice, we have found difficulties in exploring the large time domain because of the occurrence of significant nonlinear effects.

Analogous results, valid in the linear regime, can be obtained for the energy rate. In this case one finds the standard result of perturbation theory,

$$\frac{dE(t)}{dt}=\frac{V^2}{2\mathrm{}}𝑑\omega ^{}\omega ^{}\left[S(𝐪,\omega ^{})S(𝐪,\omega ^{})\right]\frac{\mathrm{sin}\left[(\omega \omega ^{})t\right]}{\omega \omega ^{}},$$
(15)
which yields the golden rule result

$$\frac{dE(t)}{dt}=\omega \left(\frac{V}{2}\right)^2\frac{2\pi }{\mathrm{}}\left[S(𝐪,\omega )S(𝐪,\omega )\right]$$
(16)
in the large $`t`$ limit. From the experimental viewpoint it is very difficult to extract the net energy transfer since it would require a very high precision in the determination of the release energy. Furthermore, the release energy does not coincide with the total energy of the sample which includes also the contribution of the confining potential.

It is worth noticing that result (11) for the momentum rate, like Eq. (15) for the energy rate, is sensitive to the difference $`[S(𝐪,\omega )S(𝐪,\omega )]`$ rather than to the dynamic structure factor itself. This difference characterizes the imaginary part of the density-density response function

$$\mathrm{Im}\left(\chi _𝐪(\omega )\right)=(\pi /\mathrm{})\left[S(𝐪,\omega )S(𝐪,\omega )\right]$$
(17)
and follows from the fact that atoms can scatter by absorbing a photon from either of the laser beams. This represents an important difference with respect to other scattering experiments (like, for example, neutron scattering from helium) where, by detecting the scattered probe, one instead measures directly the dynamic structure factor. The dynamic structure factor and the imaginary part of $`\chi _𝐪(\omega )`$ have the same behavior at $`T=0`$ since in this case $`S(𝐪,\omega )`$ is zero for negative $`\omega `$. They instead differ at finite temperatures if $`k_BT`$ is of the order or higher than the excitation energy $`\mathrm{}\omega `$. Actually the difference (17) significantly suppresses the thermal effects exhibited by the dynamic structure factor so that, by measuring $`\chi _𝐪(\omega )`$, one has an easy access to the zero temperature value of $`S(𝐪,\omega )`$. For this reason, even if experiments are carried out at temperatures which do not satisfy the condition $`kT\mathrm{}\omega `$, one can safely restrict the theoretical analysis to the simpler $`T=0`$ case, provided the temperature is sufficiently low to ignore the effects due to the thermal depletion of the condensate.

## IV Local Density Approximation

If the value of the wave vector $`q`$ is larger than the inverse of the size of the condensate the response of the system can be obtained, in first approximation, using a local density approximation (LDA). In fact, in this case the effects of discretization in the excitation spectrum can be safely ignored and the system can be treated as a locally uniform medium. In our simulation we have $`qR_{}=21`$ and $`qR_z=168`$ and consequently this approximation should be well satisfied. In the local density approximation the dynamic structure factor is given by the analytic expression

$$S_{\mathrm{LDA}}(𝐪,\omega )=\frac{15\mathrm{}}{8}\frac{(\mathrm{}^2\omega ^2E_r^2)}{E_r\mu ^2}\left[1\frac{(\mathrm{}^2\omega ^2E_r^2)}{2E_r\mu }\right]^{1/2}$$
(18)
where

$$E_r=\frac{\mathrm{}^2q^2}{2m}$$
(19)
is the free recoil energy and $`\mu =\mathrm{}^2/(2m\xi ^2)`$ is the Thomas-Fermi value of the chemical potential. In deriving result (18) one has to average the Bogoliubov expression $`S(𝐪,\omega )=[\mathrm{}^2q^2/(2mϵ(q))]\delta (\omega ϵ(q)/\mathrm{})`$ over the Thomas-Fermi density profile $`n(𝐫)=(1/g)[\mu V_{\mathrm{ext}}(𝐫)]`$, by evaluating the Bogoliubov excitation spectrum

$$ϵ(q)=\left[\frac{\mathrm{}^2q^2}{2m}\left(\frac{\mathrm{}^2q^2}{2m}+gn(𝐫)\right)\right]^{1/2}$$
(20)
at the corresponding density. In order to apply the LDA, the momentum transfer $`\mathrm{}q`$ should not be however too large, because this approximation ignores the Doppler effect associated with the spreading of the momentum distribution of the condensate which is expected to become the leading effect at very large values of $`q`$ . This happens for values of $`q`$ much larger than the inverse of the healing length, a situation that is not considered in the present work.

Equation (18) shows that, differently from the case of a uniform gas, the dynamic structure factor of a trapped condensate is no longer a delta function, its value being different from zero in the interval $`E_r<E<E_r[1+2\mu /E_r]^{1/2}`$. The value $`E=E_r`$ corresponds to the excitation energy in the region near the border where the gas is extremely dilute and hence noninteracting. The value $`E=E_r[1+2\mu /E_r]^{1/2}`$ is instead the excitation energy of a Bogoliubov gas evaluated at the central density. Notice that the LDA expression (18) for $`S(𝐪,\omega )`$ does not depend on the direction of the vector $`𝐪`$ even in the presence of a deformed trap.

Results (11) and (18) allows one to evaluate the time dependence of the momentum imparted to the condensate in the linear regime. In Figs. 1 and 2 we compare the results of this local density approximation with the ones of the numerical integration of the GP equation. The agreement is excellent and proves that the LDA is quite adequate to describe the response of the system in the conditions considered in our analysis. We have also carried out the same comparison for a smaller value of momentum transfer ($`q=1a_{}^1`$) finding a similar good agreement.

The successful comparison between the predictions of the LDA and the numerical simulation suggests that equations (18) and (11) are indeed useful for the analysis of the experimental results of Ref  (see Fig. 3). The values of the parameters used in the experiment ($`q=4.3a_{}^1`$, $`\omega _z/\omega _{}=0.12`$, $`q\xi =0.39`$) do not differ significantly from the ones used in our simulation, except for the value of the chemical potential which is significantly larger ($`\mu =60\mathrm{}\omega _{}`$). Because of this, the numerical simulation of the MIT experiment would be a highly time consuming calculation, beyond the purpose of this work. On the other hand, in the above experimental conditions the applicability of the LDA should be even better, because of the increase of the value of $`qR_{}`$ and of $`qR_z`$.

In Fig. 3 we compare the LDA prediction (solid line) with the available experimental data (points). The latter are given in arbitrary units, since the value of $`V`$ used in  is not available, and hence the comparison is limited to the shape and position of the peak. The LDA curve is obtained by integrating Eq. (11) between $`t=0`$ and $`t=\tau _B`$ and using expression (18) for the dynamic structure factor. The agreement is reasonably good. It is also interesting to compare the LDA prediction with the golden rule (14) integrated over the same time interval (dashed line). The LDA curve is broadened and quenched as a result of the finite duration of the Bragg pulse, which enters Eq. (11) through the sinusoidal factor in the integrand. In the golden rule (14) the same factor is replaced with a delta function, but the figure shows that the effect of this replacement is significant and can not be neglected.

## V Conclusions

In conclusion we have shown that the net momentum imparted to a trapped condensate by light scattering can be calculated by solving the time dependent Gross-Pitaevskii equation in a regime which is significant for current experiments. In the linear response limit, this quantity is directly connected with the dynamic structure factor for which the local density approximation (LDA) turns out to be a reliable approximation in the explored range of $`q`$’s. The connection between the dynamic structure factor and the quantities measured in stimulated light scattering experiments in the diverse accessible regimes deserves further investigations (see also Ref. ). The present work shows that, within appropriate limits, a quantitative analysis of experimental data is indeed possible. More accurate and systematic measurements, including the explicit dependence of the momentum imparted to the condensate on the duration of the pulse, would allow for a better test of the theoretical predictions for the dynamic structure factor of trapped Bose-Einstein condensates.

###### Acknowledgements.

This work is supported by the Ministero dell’Università e della Ricerca Scientifica. F.D. would like to thank the Dipartimento di Fisica di Trento for the hospitality.


text:
# Valley Splitting in Si-Inversion Layers at Low Magnetic Fields

## I Appendix: Energy spectrum in Si-inversion layer at low magnetic fields

The results on the valley splitting presented above shed a light on the energy spectrum of the two-valley electron system in (100)Si-inversion layers in low fields. We plotted in Fig. 5 the energy for three lowest Landau levels, $`N=0,1,2`$, according to Eqs. (1) and (2) (from the main section). In high fields, the cyclotron splitting is the largest and the sequence of energy levels corresponds to the inset to Fig. 1. Counting from the lowest Fermi energy, the 1st gap in the energy spectrum (filling factor $`\nu =1`$) corresponds to valley splitting; $`\nu =2`$ is for the reduced Zeeman splitting, $`\nu =3`$ is again for the valley splitting and $`\nu =4`$ is the reduced cyclotron gap. This picture was verified in numerous experiments .

Due to the nonzero valley splitting at $`B=0`$, as magnetic fields decreases, the energy levels start crossing each other at $`B1`$ T . In the region of crossing, this single-electron picture fails and the repulsion between levels should be taken into account. There is almost no experimental data available on the details of the energy spectrum in Si in low fields. It is well known only  that in low fields/low density regime quantum oscillations are missing for $`\nu =3,4,5`$ and remain pronounced for $`\nu =1,2`$ and for $`\nu =6`$ and 10; this result was recently confirmed in ref. .

In fields $`B<0.7`$ T, the sequence of energy levels changes as shown in the inset to Fig. 5. For example, at $`B=0.6`$ T, the 1st splitting, $`\nu =1`$, corresponds to the Zeeman gap, $`\nu =2`$ is the cyclotron gap (reduced by the spin flip energy), $`\nu =3`$ is a combination of the valley and cyclotron splitting etc. In this simplified picture, at some field values ($`B=0.45`$ T and 0.25 T in the inset to Fig. 5) a number of energy levels coincide. We don’t think though that the simplified single-electron picture of Fig. 5 may predict exact location of the level crossing regions.

The identification of each energy splitting in low field/low density regime is not transparent and requires more thorough theoretical calculations of the spectrum and much more detailed experiments. As an example, on the $`nB`$ plane in Fig. 6 we reproduce from Ref.  a Landau fan diagram of oscillations in the low density limit; similar data were reported recently in Ref. . In Fig. 6, the $`\rho _{xx}`$-minima at $`\nu =6`$ can be traced without interruption down to the field $`B=0.55T`$, through the ‘metal-insulator’ boundary, however the origin of this splitting may change from ‘Zeeman gap’ at high fields to any other one, e.g. ‘valley gap’ or ‘cyclotron gap’ at low fields. By now, there is no firm physical background to identify $`\nu =6`$ oscillation in low fields with spin splitting. For example, one can see from the inset to Fig. 5, that in the single-electron approximation the $`\nu =6`$ gap at $`B=0.6`$ T corresponds to the transition with changing valley (but is reduced by the spin flip energy) and $`\nu =4`$ gap is the combined gap with changing both, the cyclotron number and the spin projection.

In Ref. , from the estimated ratio of the two gaps, at $`\nu =6`$ (which authors of Ref.  assumed to be the spin-gap) and at $`\nu =4`$ (was assumed to be the cyclotron gap), a conclusion was drawn on the unexpectedly strong enhancement of the $`g^{}`$-factor at low electron densities. Taking into account vanishing $`\mathrm{\Delta }_c`$ in the crossing region and ill-defined splittings at low field, we find this argument unjustified. Recent direct measurements of the Zeeman energy in low fields vs carrier density  show that the Zeeman energy increases smoothly as density decreases down to $`n_c`$, in agreement with the anticipated Fermi liquid renormalization and with recent results by Okamoto et al. , without any unexpected deviations or divergence in the range of densities $`n=(190)\times 10^{11}`$cm<sup>-2</sup>.

It is noteworthy that for the weak field interference pattern which we analyzed above, the interpretation of each oscillation is not important and does not influence on the conclusion about the valley origin of the interference and on the results of our analysis.