text:
# Freely decaying weak turbulence for sea surface gravity waves

## Abstract

We study numerically the generation of power laws in the framework of weak turbulence theory for surface gravity waves in deep water. Starting from a random wave field, we let the system evolve numerically according to the nonlinear Euler equations for gravity waves in infinitely deep water. In agreement with the theory of Zakharov and Filonenko, we find the formation of a power spectrum characterized by a power law of the form of $`|𝐤|^{2.5}`$.

After the pioneering work by Kolmogorov  on the equilibrium range in the spectrum of an homogeneous and isotropic turbulent flow, there have been a number of studies on cascade processes in many other fields of classical physics such as plasma physics, magnetohydrodynamics and ocean waves. For surface gravity waves the first seminal theoretical work was done by O.M. Phillips in 1958, . Using dimensional arguments, he argued that the frequency spectrum in the inertial range was of the form $`F(\omega )=\alpha g^2\omega ^5`$, where $`\alpha `$ was supposed to be an absolute constant and $`g`$ is gravity. Even though in the introduction of Phillips’ paper it was stated that “a necessary condition for the equilibrium range over a certain part of the spectrum is the appreciable non-linear interactions among these wave-numbers” (from ), his arguments were based on the geometrical features of the free surface elevation. One of his basic assumptions was that the only variable of interest was gravity, while the friction velocity, $`u_{}`$, was not supposed to be involved in the spectral relation, limiting the possibility for a correct dimensional analysis.

Some years later Zakharov and Filonenko  established that in infinite water the direct cascade should produce a power spectrum of the surface elevation of the form $`P(|𝐤|)|𝐤|^{2.5}`$ that corresponds, using the linear dispersion relation in infinite depth, to an $`\omega ^4`$ frequency power spectrum: the result was found as an exact solution of the kinetic wave equation (see ). The theory developed is known as ”weak” or ”wave turbulence” and has many important applications in different fields of physics such as hydrodynamics, plasma physics, nonlinear optics, solid state physics, etc. see . It is called weak turbulence because it deals with resonant interactions among small-amplitude waves. Thus, contrary to fully developed turbulence, it leads to explicit analytical solutions provided some assumptions are made. The first experimental support of the theory for surface gravity waves was made by Toba  who was completely unaware of the paper by Zakharov and Filonenko. He reformulated the Phillips’ equilibrium range law in the following way: $`F(\omega )=\beta gu_{}\omega ^4`$, where $`\beta `$ should now be a universal dimensionless constant. After the work by Toba, successive experimental observation of the $`\omega ^4`$ law have been made by a number of authors, see for example .

Even though there is a consensus on this result, it must be stressed that so far the verification of the theory has never been established from first principles and moreover the mechanisms that lead to the power law $`\omega ^4`$ are not universally recognized: geometrical aspects related to wave breaking, without invoking the nonlinear wave-wave interaction mechanism, are still retained by many oceanographers as fundamental for generating an $`\omega ^4`$ power law. Confirmation of the Zakharov-Filonenko solution to the kinetic equation has being given through numerical simulations of the kinetic wave equation itself  , solving exactly the so called $`S_{nl}`$ term. Nevertheless, it must be underlined that the kinetic equation is derived from the primitive equations of motion under a number of hypotheses (see for example ), therefore it cannot be concluded $`a`$ $`priori`$ that power law solutions of the kinetic equation are also shared by the fully nonlinear wave equations.

One way to verify weak turbulence theory is to perform direct numerical simulations of the primitive equations of motion. The numerical confirmation of the theory for gravity waves propagating on a surface has not been an easy task (for capillary waves see , for one dimensional wave turbulence see ), basically because of the intrinsic difficulties of the computation of the boundary conditions. Different numerical approaches have been used for integrating the fully nonlinear surface gravity waves equations (see  for a review). The numerical methods based on volume formulations show very interesting results, in particular they are capable of modeling in a quite appropriate way wave breaking. Unfortunately they have the disadvantage that they require large computational resources, and therefore are not suitable for long time numerical simulations. For irrotational and inviscid flows boundary formulations are usually preferred: only the surface is discretized reducing the dimension of computation (from three to two). The Higher-Order Spectral Methods (HOS), indeed the method used in our numerical simulations, introduced independently by West et al.  and by Dommermuth et al. , belongs to this second approach (see also the recent work byTanaka ). Very recently three new methods have been proposed as very promising for simulating water waves . Results using these new approaches on turbulent cascades are still to be completed.

In this Letter we establish numerically, using a HOS method, that nonlinear interactions are sufficient for generating power laws in wave spectra; moreover we show that the Zakharov-Filonenko theory is completely consistent with the primitive equations of motion. We consider a system of random waves localized in wave number space and we show how nonlinearities “adjust” the spectrum in agreement with the Zakharov and Filonenko prediction. Numerical work in the case of a forced and dissipated system has been attempted by Willemsen  using what sometimes are called the “Krasitskii equations” (see also ). In order to avoid the effects of external forcing, we considered the case of a freely decaying wave field. If the simulations, as we will see, show the formation of a power law then the conjecture that this power law is caused by geometrical features related to forcing and wave breaking must be excluded, since forcing is absent and wave breaking cannot be taken into account using the numerical method considered. From a physical point of view, the freely decaying case corresponds to the evolution of a swell wave field. A generic wave field is considered at time $`t=0`$ and it is allowed to evolve in a natural way using a high order approximation of the Euler equations. Since numerical computations are limited by the dimension of the grid considered, an artificial dissipation is needed at high wave numbers in order to prevent accumulation of energy and a break down of the numerical code. The fluid is considered inviscid, irrotational and incompressible. Under these conditions the velocity potential $`\varphi (x,y,z,t)`$ satisfies the Laplace’s equation everywhere in the fluid. The boundary conditions are such that the vertical velocity at the bottom is zero and on the free surface the kinematic and dynamic boundary conditions are satisfied for the velocity potential $`\psi (x,y,t)=\varphi (x,y,\eta (x,y,t),t)`$ (we assume that fluid is of infinite depth):

$$\psi _t+g\eta +\frac{1}{2}\left[\psi _x^2+\psi _y^2(\varphi _z|_\eta )^2(1+\eta _x^2+\eta _y^2)\right]=0$$
(1)
$$\eta _t+\psi _x\eta _x+\psi _y\eta _y\varphi _z|_\eta (1+\eta _x^2+\eta _y^2)=0,$$
(2)
The major difficulty for numerical simulations of the system (1)-(2) consists in that we have to compute the derivatives of $`\varphi `$ with respect to $`z`$ on the surface $`\eta `$. This problem can be overcome if we express the velocity potential $`\psi (x,y,t)`$ as a Taylor expansion around $`z=0`$. Inverting asymptotically the expansion one can express $`\varphi _z|_\eta `$ as an expansion of derivatives of $`\psi (x,y,t)`$ that can then be computed using the Fast Fourier Transform, simplifying notably the computation. This is nothing other than a different way for formulating the HOS method. We underline that this is the same approach that has originally been used in  for deriving analytically the equation that is usually known as the “Zakharov equation”. The order of the simulation can be decided a priori and depends on how many terms are retained in the Taylor expansions; in our numerical simulations we considered the expansion necessary to take into account four wave interactions so that we are consistent with the order of the “Zakharov equation”.

A delicate point in our numerical simulations is related to the dissipation of energy at high wave numbers. We remark that this dissipation is completely artificial since we are dealing with a potential flow. Nevertheless we have considered the dissipation phenomenon of the wave field to be similar to the one that takes place in a turbulent flow, i.e. that is mathematically expressed by a Laplacian that operates on the velocity. As is usually done in direct numerical simulations of box turbulent flows, in order to increase the inertial range, we have used a higher order diffusive term. More explicitly on the right hand side of equation (1)-(2), we have added respectively two extra terms: $`\nu (^2)^n\psi `$ and $`\mu (^2)^m\eta `$, where $`\nu `$ and $`\mu `$ represent an artificial viscosity coefficient and $`^2`$ is the horizontal Laplacian. If $`n`$ and $`m`$ are greater than 1 the viscosity is known as “hyperviscosity”.

It has to be noted that, at first sight, one would use a very high power of the Laplacian in order to increase notably the inertial range, unfortunately very high values of $`m`$ and $`n`$ could bring about the “bottleneck effect” , i.e. an accumulation of energy at high wave numbers that could distort the power law expected . In our numerical simulations we used $`\nu =\mu =3\times 10^4`$ and $`n=m=8`$. These values have been selected after some trial and error during the development of the numerical code: because of our limitated number of grid points, smaller values of $`m`$ and $`n`$, would obscure almost completely the inertial range. In our numerical simulations we did not impose any a dissipation at low wave numbers.

In order to prepare the initial wave field it is reasonable to consider a directional spectrum $`S(|𝐤|,\theta )=P(|𝐤|)G(\theta )`$. The directional spreading function $`G(\theta )`$ used here is a cosine-squared function in which only the first lobe (relative to the dominant wave direction) is considered:

$$G(\theta )=\{\begin{array}{cc}\frac{1}{\sigma }cos^2\left(\frac{\pi }{2\sigma }\theta \right)\hfill & if\sigma \theta \sigma \hfill \\ 0\hfill & else\hfill \end{array}$$
(3)
$`\sigma `$ is a parameter that provides a measure of the directional spreading, i.e. as $`\sigma 0`$, the waves become increasingly unidirectional. In our numerical simulations we selected the value of $`\sigma =\pi /2`$. We tried to avoid the complete isotropic case in order to verify if the theory still holds for intermediate values of the spreading. At the same time the selection of a large value of $`\sigma `$ was motivated by the fact that recently it has been found  that, for sufficiently narrow angle of spreading, the Benjamin-Feir instability can be responsible for the formation of freak waves. As a consequence the nonlinear energy transfer could be slightly altered and some corrections to the prediction could be necessary (this very interesting topic is now under investigation and results will be reported in a different paper). $`P(|𝐤|)`$ is chosen to be any localized spectrum. We have performed numerical simulations with a gaussian function or with a “chopped JONSWAP” spectrum (a JONSWAP spectrum with amplitudes equal to zero for frequencies greater than 1.5 times the peak frequency) with random phases. For the case of the gaussian function, wave numbers lower than a selected threshold have been set to zero in order to avoid extremely long large waves. The velocity potential is then computed from the initial wave field using the linear theory. Both gaussian and JONSWAP spectra led to the same results in terms of the turbulent cascade.

Our computation is performed in dimensional units; we have selected the initial spectrum centered at 0.1 Hz, i.e. we are considering 10 seconds waves. The initial steepness computed as $`\epsilon =k_0H_s/2`$ was chosen to be around 0.15 ($`H_s`$ was computed as 4 times the standard deviation of the wave field). The wave field was contained in a square grid (the resolution is $`256\times 256`$) of length $`L=1417.6`$ meters. The time step considered was $`1/50`$ the dominant frequency, i.e. $`\mathrm{\Delta }t=0.2`$ seconds. We have performed our numerical simulations on a 400Mh PC. In Fig. 1 we show the evolution of the wave power spectrum for different real times (t=0, 0.1, 0.5, 1 hours). We see that, as expected, the tail of the spectrum starts to grow. This process seems to be quite fast: as is shown in the figure after a few dominant wave periods some energy is already injected into high wave numbers. The process of adjusting the power law to the “correct” one becomes then very slow, especially for low wave number. This could be due to the frozen turbulent phenomenon , i.e. a condition in which the energy fluxes towards high wave numbers are reduced because of the discretness of the spectrum. Moreover decaying numerical simulations are very time consuming with respect to forced simulations because, as time passes, energy is lost due to viscosity, thus reducing the significant wave height of the wave field and therefore the steepness.

Even though it is not clear from the Log-Log representation in Fig. 1, there is a downshifting of the peak of the spectrum towards lower wave numbers; as a consequences the steepness subsequently decreases over time. The time scale of the nonlinear energy transfer becomes larger and larger. In Fig 2 we show the power spectrum of the surface elevation after 4 hours (the steepnes of the wave field is $`\epsilon 0.07`$). In the same plot we show two power laws $`k^{2.5}`$ and $`k^3`$: the first one seems to better fit the data. In order to be completely sure that the numerical data are in agreement with the prediction of Zakharov and Filonenko, we show in Fig. 3 compensated spectra with different compensation powers: $`z=2.5`$ seems to be the most plausible power law. Thus there seems to be ample evidence from our numerical simulations that the power law is in sufficiently good agreement with the value predicted by the theory.

After the pioneering work by Zakharov and Filonenko the kinetic wave theory has developed further, making available quantitative predictions for other physical observables such as energy fluxes, downshifting of the peak, energy dissipation etc. All these quantities will be examined and results will be reported in future papers. Other questions naturally arise from our results: in HOS simulations the order of the computation depends on how many terms are retained in the Taylor expansion; do higher order terms influence the cascade process? Our computation has been performed in a freely decaying case; could external forcing (especially if anisotropic) influence the power law? And more, what would be the influence of the water depth? These are all questions to be answered in the near future.

Acknowledgements This work was supported by the Office of Naval Research of the United States of America (T. F. Swean, Jr.) and by the U.S. Army Engineer Research and Development Center. Consortium funds and Torino University funds (60 %) are also acknowledged. M. O. was also supported by a Research Contract from the Università di Torino.


format: text

text:
# Cosmic rays above the ankle from Z-bursts

## 1 Introduction

The existence of ultrahigh-energy cosmic rays (UHECR) with energies above the Greisen-Zatsepin-Kuzmin (GZK) cutoff  of about $`5\times 10^{19}`$ eV, presents an outstanding problem . Nucleons and photons with those energies have short attenuation lengths and could only come from distances of 100 Mpc or less , while plausible astrophysical sources for those energetic particles are much farther away.

An elegant and economical solution to this problem, proposed by T. Weiler , consists of the production of the necessary photons and nucleons close to Earth, in the annihilation at the $`Z`$-resonance of ultrahigh-energy neutrinos coming from remote sources, $`\nu _{\mathrm{UHE}}`$, and relic background neutrinos. These events were named “$`Z`$-bursts” by T. Weiler. One of the most appealing features of the “$`Z`$-bursts” scenarios is that the energy scale of $`10^{2021}`$ eV, at which the unexpected events have been detected, is generated naturally given the possible mass range of relic neutrinos. The $`Z`$-resonance occurs when the energy of the incoming $`\nu _{\mathrm{UHE}}`$ is $`E_{\nu _{\mathrm{UHE}}}=E_{Res}`$,

$$E_{Res}=\frac{M_Z^2}{2m_\nu },$$
(1)
where $`m_\nu `$ is the mass of the relic neutrinos. This is the new cutoff of the UHECR energy in these models. It depends inversely on the mass of the relic neutrinos. Since arguments of structure formation in the Universe show $`m_\nu <`$ few eV, then $`E_{\nu _{\mathrm{UHE}}}>10^{21}`$ eV, precisely above the GZK cutoff, as needed.

In this paper we concentrate on a particular “Z-bursts” scenario , in which the relic neutrinos are lighter than 1 eV. These lighter neutrinos, contrary to those in the original scenario, cannot be gravitationally bound, they have a constant density over the whole Universe. In particular, we concentrate on relic neutrinos with mass compatible with Super-Kamiokande results, if neutrino masses are hierarchical (however the results we obtain apply to heavier relic neutrinos, while light enough to not cluster). Super-Kamiokande has provided a strong evidence for the oscillation in atmospheric showers of two neutrino species with masses $`m_1`$, $`m_2`$ and $`\delta m^2`$ = $`m_1^2m_2^2`$ = $`(18)\times 10^3`$ eV, consisting mostly of about equal amounts of $`\nu _\mu `$ and another flavor eigenstate neutrino, $`\nu _\tau `$ or a sterile neutrino . If neutrino masses are hierarchical, as those of the other leptons and quarks, then the heavier of the two oscillating neutrinos, call it $`\nu _{\mathrm{SK}}`$, has a mass $`m_{\mathrm{SK}}=\sqrt{\delta m^2}0.07`$ eV. With $`m_\nu =m_{\mathrm{SK}}`$, the new UHECR cutoff becomes

$$E_{Res}0.6\times 10^{23}eV.$$
(2)
Due to the large multiplicity of the Z-decays, after energy losses in the propagation (as shown in detail here) this value of $`E_{Res}`$ predicts many super GZK UHECR events (many more than with eV relic neutrino masses).

We agree with Farrar and Piran , who have argued that any mechanism accounting for the events beyond the GZK cutoff should also account for the events down to the ankle, including their isotropy and spectral smoothness. We show here that the model we consider can account for the ankle and all the events above it if the position of the ankle is close to that measured by AGASA, $`E_{ankle}=10^{19.0}`$ eV (see , in particular Table V, and references therein). Moreover, a reliable prediction of the model is that most primaries above the ankle should be nucleons up to about $`10^{20.0}`$ eV and photons at higher energies. We also find that nucleons do accumulate at the GZK cutoff energy, which could account for the hint of a slight accumulation seen in the data. Photons become dominant at energies higher than $`10^{20}`$ eV in our model. So these photons are all above the threshold energy (which is about $`5\times 10^{19}`$ eV) to pair produce in the Earth’s magnetic field (which should generate a certain amount of north-south asymmetry in the arrival direction distribution).

Let us return to the issue of the isotropy of the events above the ankle, i.e., the absence of a correlation with the galactic halo. Because the relic neutrinos we assume do not gravitationally cluster, the isotropy of the events reflects the isotropy of the ultrahigh-energy neutrino sources. In particular, relic neutrinos of mass $`m_{\mathrm{SK}}`$ require a large flux of neutrinos with energies of the order of $`10^{23}`$ eV. It is unlikely that active galactic nuclei , neutron stars , or other astrophysical sources could produce such a high energy flux of ultrahigh-energy neutrinos. Topological defects , or superheavy relic particles , could instead easily generate the requisite flux of primary neutrinos (but there are still problems with these sources). For example, with unstable superheavy relic particles, which form part of the cold dark matter, decaying mostly into neutrinos , the directions of UHECR could map the distribution of parent particles (which should coincide with the distribution of cold dark matter) at large redshifts. This is because the initial energy of the $`\nu _{\mathrm{UHE}}`$ produced in the decay needs to be redshifted to the energy of the “$`Z`$-burst” in its way to the Earth. Directional clustering, evident in the existing data , would then reflect the distribution of matter at a certain red shift determined by this process of “cosmological filtering”. Thus, absence of directional correlations with the galactic halo, as well as directional clustering, could be easily accommodated .

Besides, $`\nu _{\mathrm{UHE}}`$ produced by unstable superheavy relic particles would have a spectrum opposite to an astrophysical spectrum, growing rapidly with energy, up to a sharp cutoff at an energy of the order of the parent particle mass. This spectrum has almost no neutrinos at low energy where bounds exist . Most bounds on “$`Z`$-bursts” models (see for example ) assume that the $`\nu _{\mathrm{UHE}}`$ have a typical astrophysical spectrum, decreasing with energy as $`E^\gamma `$, with $`\gamma `$ a number of order one. These bounds do not hold if the $`\nu _{\mathrm{UHE}}`$ spectrum has a very different energy dependence. However, a model for these parent particles is arguably difficult to obtain . Moreover, the EGRET bound on the diffuse low-energy gamma ray flux resulting from the “$`Z`$-bursts” imposes important constraints , which might rule out heavy particles decaying mostly into neutrinos as sources .

In the next section we present our simulations and the resulting spectrum of UHECR. We would like to point out that the main result of this paper, the fact that “$`Z`$-bursts” can account for the ankle and all events above it, does also hold for larger relic neutrino masses, for which the problem of the sources becomes less severe. In fact even if we used $`m_{\mathrm{SK}}`$ here, our considerations apply with trivial changes to other masses for which relic neutrinos are too light to gravitationally cluster. As the relic neutrino mass increases, all the features in the spectrum we find here should move progressively to lower energies.

## 2 Spectrum of UHECR from “$`Z`$-bursts”

We have performed simulations of the photon, nucleon and neutrino fluxes coming from a uniform distribution of “$`Z`$-bursts”, namely $`\nu \overline{\nu }`$ annihilations at the Z pole ($`\nu \overline{\nu }Zp\gamma \mathrm{}`$), with the energy of Eq. 2, corresponding to relic neutrinos of mass $`m_{\mathrm{SK}}`$. The “$`Z`$-bursts” were simulated using PYTHIA 6.125 , and the absorption of photons and nucleons was modeled using energy attenuation lengths provided by Bhattacharjee and Sigl , supplemented by radio-background models by Protheroe and Biermann .

We simulated a uniform distribution of about $`10^7`$ Z particles up to a maximum $`z_{max}=2`$. Even if the shape of the spectrum of nucleons and photons changes somewhat with other Z particle distributions, the main features of the spectrum stated here remain the same.

The decay of the Z bosons through all possible channels was done automatically by the PYTHIA routines , using the default options of this program. For comparison with the other figures we show in Fig. 1 the spectra given by PYTHIA, normalized per Z boson (not including redshifts and energy absorption). The multiplicities that PYTHIA gives per Z-decay are (in each case counting particles and antiparticles) 1.6 nucleons, 17 photons, 15 $`\nu _e`$, 30 $`\nu _\mu `$ and 0.23 $`\nu _\tau `$.

In our simulation, each Z boson generated by PYTHIA was placed on the “event list” of the cascade generator at a randomly selected distance. The cascade of decay products was then boosted. The $`\gamma `$ factor corresponding to an energy $`E_\nu =E_{res}`$ is

$$\gamma =\frac{E_\nu +m_\nu }{M_z}=6.25\times 10^{11}.$$
(3)
The gamma factor of each boost was actually corrected to include the redshift of the decaying Z particle. We then followed the propagation of the nucleons, photons and neutrinos resulting from the boosted Z decays. Neutrinos do not interact in their propagation. Thus, the energy spectra of the three kinds of neutrinos were simply generated by counting the number of particles per energy bin and normalizing this number to the total number of Z particles used.

We included energy absorption for nucleons and photons. Each nucleon or photon was created by PYTHIA in the initial cascade at a fixed position, with fixed energy and direction of motion with respect to the Earth frame of reference. The distance each particle had to travel before reaching Earth was compared with the appropriate attenuation length in space for the particle energy. If the distance was smaller than the attenuation length, the particle was assumed to reach Earth unchanged. In the opposite case, the energy and momentum of the particle were degraded by a factor $`1/e`$ after traveling a distance equal to the attenuation length (and the particle was placed again in the list constituting the cascade at the new position, with the degraded energy and momentum).

This process was continued until all nucleons and photons reached Earth and were counted in the final spectra, or until their energies became too small to be significant, in which case they were simply discarded from the cascade. At this point, the final nucleon and photon spectra were normalized to the total number of Z particles used (the same done with the neutrino spectra). The results are given in Fig. 2, with an approximate fit to the AGASA cosmic ray data .

We used the nucleon energy attenuation length given by Bhattacharjee and Sigl in the Fig. 9 of Ref. , based on results from Ref.  and .

The energy attenuation length of photons is poorly known, due to the uncertainty in the radio background. Using the attenuation length shown in Fig. 11 of Ref. , based on observations of Clark et al. , the resulting photon flux is shown in Figs. 2 and 3 as curve (1). Protheroe and Biermann  produced two models for the radio background which lead to shorter interaction lengths than those based on Clark et al. From the provided interaction lengths we constructed approximate attenuation lengths for the models of Protheroe and Biermann by reducing the attenuation length based on Clark et al. by the difference between the interaction lengths. With the attenuation lengths so constructed we obtained the curves (2) and (3) in Figs. 2 and 3. This is obviously only an approximation, since the mean interaction and energy attenuation lengths do not have exactly the same energy dependence. However, we believe the three curves we obtained give a good representation of the possible range of predicted photon fluxes, in view of the uncertainties related to the energy attenuation of photons in space.

We have arbitrarily used the middle photon flux, curve number (2), when computing the total flux. We added the proton and photon fluxes obtained in this paper to a power law spectrum with slope $`3.23`$ found by AGASA to fit the data below the ankle (from $`10^{17.6}`$ to $`10^{19.0}`$ eV; see Table V of ). Our results depend very little on which of the three photon fluxes we use.

The fit of the AGASA data with our total flux provides the normalization of the photon and nucleon differential fluxes $`F`$, denoted as $`d\varphi /dE`$ in the figures (in this case $`F_{AGASA}`$ = $`10^{14.2}(\mathrm{m}^2\mathrm{sr}\mathrm{s})^1`$ $`F_{PYTHIA}`$ = $`6\times 10^{15}(\mathrm{m}^2\mathrm{sr}\mathrm{s})^1`$ $`F_{PYTHIA}`$), that allows us to determine the (assumed homogeneous and isotropic) flux of ultrahigh-energy neutrinos close to the $`Z`$-resonance energy of Eq. 2 (at some energy between $`E_{Res}/(1+z_{max})=E_{Res}/3`$ and $`E_{Res}`$) to be

$$F_{\nu _{\mathrm{UHE}}}1\times 10^{36}\frac{1}{\mathrm{eV}\mathrm{m}^2\mathrm{sr}\mathrm{s}},$$
(4)
if no lepton asymmetry is assumed in the neutrino background. This flux is shown in Figs. 2 and 3, with the label “$`\nu _{\mathrm{UHE}}`$”. With the level of accuracy of our simulation we can only determine the order of magnitude of this flux. Moreover, for the light neutrinos we are considering, present experimental bounds allow for a lepton asymmetry which could increase the number of relic neutrinos by a factor smaller than 10, reducing the necessary $`F_{\nu _{\mathrm{UHE}}}`$ by the same factor.

We have continued the power law spectrum accounting for the events below the ankle, presumably due to galactic sources, way beyond the ankle, while this contribution may die out at or soon above the ankle. In any event, without adding up the power law, i.e., taking into account only the photon and nucleon fluxes from “$`Z`$-bursts” computed here, we had made a very similar fit , with which we had obtained the same $`F_{\nu _{\mathrm{UHE}}}`$.

As we just mentioned, without making an assumption on the spectrum of $`\nu _{\mathrm{UHE}}`$, the “$`Z`$-bursts” mechanism provides an estimate of the $`\nu _{\mathrm{UHE}}`$ flux only close to the $`Z`$-resonance energy $`E_{Res}`$ (at some energy between $`E_{Res}/(1+z_{max})=E_{Res}/3`$ and $`E_{Res}`$) . It is interesting to point out that, on the basis of that one estimate, the “Goldstone Experiment”(GLUE)  could start testing our model after about 300 hours of observation. This can be seen in Fig. 3, where we show the present limits of this experiment, based on about 30 hours of observation. In Fig. 3 the fluxes of Fig. 2 have been multiplied by $`E^2`$. GLUE searches for lunar radio emissions from interactions of neutrinos above $`10^{19}`$ eV of energy and is expected to accumulate 120 hours of observation later this year.

Figs. 3 and 4 include Akeno data taken from Fig. 23 of Ref. , originally from Ref.  (besides the AGASA data). In Fig. 4 we show our nucleon flux, only our middle photon flux (curve (2) in previous figures), the power law providing a good fit to data below the ankle, and the total flux (the sum of these three components).

We believe that Fig. 4, in which the plotted spectra have been multiplied by $`E^3`$, shows best the change of slope close to the position of the ankle as measured by AGASA, $`E_{ankle}=10^{19.0}`$ eV. Obtaining a value of $`E_{ankle}`$ a factor of 3 smaller, close to $`10^{18.5}`$ eV as measured by Fly’s Eye (see Table V of ) may require to increase the $`z_{max}`$ of the source distribution considerably (which would point to particles as sources), since this would move the lower energy edge of the nucleon flux to lower energies. Alternatively, larger values of the relic neutrino mass may also work to lower $`E_{ankle}`$, since all the features in the spectrum should move to lower energies (even if we used $`m_{\mathrm{SK}}`$ here, our considerations apply with trivial changes to other masses for which relic neutrinos are too light to gravitationally cluster).

In Fig. 4 one can clearly see the enhancement of the predicted flux at the GZK cutoff energy, at about $`5\times 10^{19}`$ eV, due to the accumulation of nucleons, which can also be seen in the AGASA data.

## 3 Conclusions

In this paper we considered a particular “Z-bursts” scenario , in which the relic neutrinos are lighter than 1 eV, and thus do not gravitationally cluster. Using in particular 0.07 eV relic neutrinos, we have shown that “$`Z`$-bursts” may account not only for the super GZK events, but for the “ankle” and all UHECR events above it, including their isotropy and spectral smoothness. In our simulation we found the “ankle” close to $`E_{ankle}=10^{19.0}`$ eV as measured by AGASA. A reliable prediction of the model is that most primaries above the ankle should be nucleons up to about $`10^{20.0}`$ eV and photons at higher energies. Moreover, the nucleons do accumulate at the GZK cutoff energy, which could account for the hint of a slight accumulation seen in the data. The model predicts a new cutoff, which with 0.07 eV relic neutrinos is at $`E_{Res}0.6\times 10^{23}eV`$.

We have not included the effect of extragalactic magnetic fields, thus for the predictions of this paper to be true these fields should be sufficiently small, probably smaller than $`10^9G`$.

Finally let us comment on recent related papers. Photon and nucleon spectra very similar to those presented here are shown in Fig. 2a of Ref.  corresponding to “$`Z`$-bursts” with 0.1 eV relic neutrinos and a different model for the distribution of “$`Z`$-bursts” with redshift up to $`z_{max}=3`$. This model satisfies the EGRET bound on low energy photons (even if with astrophysical sources emitting only UHE neutrinos). This reassures us that the result we present here is robust. The EGRET bound has been computed for various “$`Z`$-bursts” scenarios , including the particular one we concentrated on here , which seems to work well with sources emitting only UHE neutrinos. The most serious problem with these sources may be the electroweak cascading of the UHE neutrinos produced in the decays, as recently claimed in .

Events above the ankle have been previously fitted with “$`Z`$-bursts” products plus an additional hypothetical component of galactic or extragalactic protons in Ref. , varying the relic neutrino mass, in an attempt to provide a determination of this mass using UHECR data. However the normalization and slope of the extra proton flux and the normalization of the nucleon flux from “$`Z`$-bursts” (photons from “$`Z`$-bursts” were neglected) were taken as fitting parameters. This procedure does not make clear if it is actually the “$`Z`$-bursts” which account for the change of slope at the “ankle” and for the events above the “ankle”. Moreover, with this procedure , there is no prediction of the position of the ankle, since this is one of the parameters to be fitted by the sum of the mentioned two fluxes. Here we took instead the measured flux below the ankle, with the measured slope and normalization. We considered it to be of galactic origin, as the correlation with the galactic center of the arrival directions measured by AGASA around $`10^{18}eV`$  seems to indicate. Then we added to this fixed flux the new component due to “$`Z`$-bursts”, which led us to find a prediction for the position of the ankle, and the normalization of the flux of primaries due to “$`Z`$-bursts”.

We believe our flavor of “$`Z`$-bursts” provides a plausible explanation to the puzzle of ultrahigh-energy cosmic rays, not only for the super GZK events, but for the “ankle” and all events above it.

This work was supported in part by the US Department of Energy grant DE-FG03-91ER40662, Task C. G.V. was also supported by an award from Research Corporation. We thank A. Kusenko and S. Nussinov for many valuable discussions, and P. Biermann and E. Roulet for comments and suggestions.


text:
# Voltage-controlled wavelength conversion by terahertz electro-optic modulation in double quantum wells

## Abstract

An undoped double quantum well (DQW) was driven with a terahertz (THz) electric field of frequency $`\omega _{THz}`$ polarized in the growth direction, while simultaneously illuminated with a near-infrared (NIR) laser at frequency $`\omega _{NIR}`$. The intensity of NIR upconverted sidebands $`\omega _{sideband}=\omega _{NIR}+\omega _{THz}`$ was maximized when a dc voltage applied in the growth direction tuned the excitonic states into resonance with both the THz and NIR fields. There was no detectable upconversion far from resonance. The results demonstrate the possibility of using gated DQW devices for all-optical wavelength shifting between optical communication channels separated by up to a few THz.

A basic function in a wavelength division multiplexed (WDM) optical communications network is to switch data between different wavelength channels. All-optical switching in a WDM network requires the ability to shift the frequency of a near-infrared (NIR) carrier of a NIR carrier by several terahertz (THz). This has been accomplished by four-wave mixing of NIR beams in semiconductor optical amplifiers, and three-wave mixing in quasi-phase-matched AlGaAs waveguides.

Generation of THz optical sidebands on a NIR carrier beam has recently been studied in various semiconductor systems. Great flexibility for practical optoelectronic devices that operate at THz frequencies lie with gated double quantum well (DQW) structures, in which two quantum wells (QW) are separated by a thin tunnel barrier. Their intersubband spacings can be engineered to lie in the THz frequency range by adjusting the width of the tunnel barrier. Meanwhile, the optical bandgap can be separately tuned by the width of the individual wells. Large optical nonlinearities can be built into the DQW by using an asymmetric DQW in which carriers feel a non-centrosymmetric confinement potential. Finally, the intersubband spacing can be tuned by applying a DC voltage to the gates.

This paper describes experiments in which a near-infrared (NIR) probe laser beam is modulated at THz frequencies in a gated, asymmetric DQW. The THz field couples to an excitonic intersubband excitation while the NIR field couples to an excitonic interband excitation. Applying a DC voltage to the gates can bring the intersubband transition into resonance with the THz field, and the interband transition into resonance with the NIR field. When these resonance conditions are met the NIR probe is modulated resulting in the emission of optical sidebands which appear at frequencies

$$\omega _{sideband}=\omega _{NIR}+n\omega _{THz}$$
(1)
where $`\omega _{NIR}`$ ($`\omega _{THz}`$) is the frequency of the NIR (THz) beam and $`n=\pm 1,2,\mathrm{}`$.

The sample consisted of an active region containing the DQWs, a distributed Bragg reflector (DBR), and two gates, each consisting of a narrow n-doped quantum-well. The band diagram of the sample is shown in Fig. 1.

The active region consisted 5 periods of DQW, each consisting of a 120 Å GaAs QW and a 100 Å GaAs QW separated by a 25 Å Al0.2Ga0.8As tunnel barrier. Each period is separated by a 300 Å Al0.3Ga0.7As barrier. The active region is thin compared to a NIR wavelength. A thicker active region generates stronger sidebands, but is more difficult to model due to reabsorption, phase-matching, and screening effects.

The dimensions of the DQW were designed so that the two lowest-lying electron subbands were separated by $`10`$ meV ($`2.4`$ THz) at flat-band. This is around the center of the range of high-quality operating frequencies of our THz source, the UCSB Free Electron Laser (FEL). Within the 10 meV energy splitting constraint, the degree of asymmetry was designed to peak the strength of the sideband near flat-band conditions.

The DBR consisted of 15 periods of 689 Å AlAs and 606 Å Al0.3Ga0.7As. It had a low-temperature passband nearly centered on the low-temperature bandgap of the DQW, making it about 95% reflective for the NIR probe beam and sidebands.

The active region was sandwiched between the two gates, each consisting of a Si delta-doped 70 Å QW with carrier density $`1\times 10^{12}`$ cm-2. The gates are separated from the active region by 3000 Å Al0.3Ga0.7As barriers. Since the gate QWs are much narrower than the active DQWs, the gate QWs are transparent to both the THz and NIR beams. A mesa was etched and NiGeAu ohmic contacts annealed to the gate QWs.

The experimental setup is illustrated in Fig. 2. The sample was cooled to 21 K in a closed-cycle He crystat. The THz beam from the FEL is polarized in the growth direction, propagated in the QW plane and focused by a 90o off-axis parabolic mirror (F/1) onto the cleaved edge of the sample.

NIR light from a continuous-wave Ti:sapphire laser is chopped by an acousto-optic modulator into $`25`$ $`\mu `$s pulses which overlap the $`5`$ $`\mu `$s FEL pulses. The vertically-polarized NIR beam propagates normal to the THz beam, and was focused (F/10) at a power density of $``$ 50 W/cm2 to the same small interaction volume in the sample. The reflected beam, sidebands, and photoluminescence (PL) are analyzed by a second polarizer, dispersed by a 0.85 m Raman double-monochromator, and detected by a photomultiplier tube (PMT).

A typical $`n>0`$ sideband spectrum is shown in Fig. 3. The conversion efficiency scale was calibrated by attenuating the NIR laser to match the peak sideband signal, and factoring the loss through the collection and dispersion optics. We present results for the $`n=+1`$ sidebands as a continuous function of gate bias ($`V_g`$) and NIR laser frequency ($`\omega _{NIR}`$) at various THz laser frequencies ($`\omega _{THz}`$).

A sideband voltage scan with $`\omega _{NIR}=1546`$ meV and $`\omega _{THz}=2.5`$ THz is shown in Fig. 4. In this scan, $`V_g`$ is scanned while all other parameters are kept constant so that $`\omega _{detect}=\omega _{NIR}+\omega _{THz}`$. Note that at dc electric fields far from resonance there are no detectable sidebands. This type of scan demonstrates the voltage tunability of the THz-NIR modulation.

The distinct peaks are labelled with excitonic transition assignments derived from a nonlinear susceptibility calculation for THz-NIR mixing due to excitons. The label $`E_\mu HH_\nu X`$ refers to the exciton consisting of an electron in subband $`\mu `$ and a heavy hole in subband $`\nu `$. A double-resonance condition holds when the NIR field is resonant with an exciton and the THz field resonantly couples two excitons.

To understand the full resonant structure of the sideband generation process at each THz frequency, the sideband was measured as a continuous function of both $`\omega _{NIR}`$ and $`V_g`$. By taking a sideband voltage scan at each $`\omega _{NIR}`$ we measured the sideband maps shown in Figs. 5 and 6 at $`\omega _{THz}=`$2.0 THz and 2.5 THz (8.2 meV and 10.4 meV). Again, peaks in the sideband emission occur when the NIR field is resonant with an exciton and the THz field resonantly couples two excitons; they are labelled with excitonic transition assignments.

In summary, we have modulated a NIR laser beam at THz frequencies by driving the excitonic intersubband transition of an asymmetric DQW. A pair of n-doped QW gates which are transparent to both the NIR and THz fields allow the resonant response of the device to be tuned by applying a dc voltage, allowing the device to act as a voltage-controllable wavelength converter. The peak conversion response occurs when the NIR field is resonant with an exciton and the THz field resonantly couples two excitons.

We gratefully thank Christoph Kadow and Art Gossard for their input. This research was funded by NSF-DMR 0070083.


text:
# Demonstration of the Lateral Casimir Force

## Abstract

The lateral Casimir force between a sinusoidally corrugated gold coated plate and large sphere was measured for surface separations between 0.2 $`\mu `$m to 0.3 $`\mu `$m using an atomic force microscope. The measured force shows the required periodicity corresponding to the corrugations. It also exhibits the necessary inverse fourth power distance dependence. The obtained results are shown to be in good agreement with a complete theory taking into account the imperfectness of the boundary metal. This demonstration opens new opportunities for the use of the Casimir effect for lateral translation in microelectromechanical systems.

PACS numbers: 12.20.Fv, 42.50.Lc, 61.16.Ch

The archetypical Casimir force  leads to an attraction between two neutral metal plates placed in vacuum. The force results from the alteration by the metal boundaries of the zero point electromagnetic energy which is present in empty space. Uniquely the Casimir force is independent of the electric charge and other interaction constants. As a result, it is strongly dependent on geometry and topology of the boundary and can be both attractive as well as repulsive (for a review of different aspects of the Casimir effect see monographs ).

The normal Casimir force which leads to an attraction perpendicular to the two surfaces has been demonstrated first between flat plates  and between a flat plate and a lens . Recently extensive experimental research on the normal Casimir force between a large sphere and a flat plate has been performed with increased precision . Also the non-trivial boundary dependence of the Casimir force acting between a large sphere and plate with periodic uniaxial sinusoidal corrugations was demonstrated in Ref. . This has led to extensive theoretical study of the corrections to the Casimir effect due to various factors as finite conductivity of the boundary metal, surface roughness, and nonzero temperature. Also the combined effect of these factors was investigated (new experimental and theoretical developments are presented in review ). Even more importantly, the Casimir effect is finding new applications in fundamental science and engineering. Thus, with the advent of modern unified theories involving compact dimensions, precision measurements of the Casimir force have been used to set limits on the presence of hypothetical forces . With regard to technological applications, both static and dynamic microelectromechanical machines have been created actuated by the normal Casimir force . Also, adhesion and sticking of moving parts in microelectromechanical systems due to the Casimir effect were investigated .

In this letter we report the first demonstration of the lateral Casimir force. It acts between two aligned corrugated surfaces and leads to a mechanical force acting tangential to the surfaces. Similar to the normal Casimir force, the lateral Casimir force also originates from the modifications of electromagnetic zero point oscillations. The possibility of a lateral Casimir force for anisotropic boundaries was investigated theoretically and the harmonic dependence on a correspondent angle was predicted (see  for the plates made of anisotropic material and  for the corrugated plates made of ideal metal). Note that other motional frictional forces between two flat parallel surfaces have been suggested but they would be several orders of magnitude smaller . Here we compare the experimental data with a complete theory taking into account the finite conductivity corrections of the boundary metal. The predicted sinusoidal dependence of the lateral force on the phase shift and inverse fourth power dependence on surface separation is confirmed and good agreement with a complete theory is demonstrated. This demonstration opens new opportunities for the actuation of lateral motion in microelectromechanical systems based entirely on the vacuum effects of quantum electrodynamics.

The regularized zero point energy per unit area for two parallel plates of infinite conductivity a distance $`z`$ apart is given by $`U(z)=\pi ^2\mathrm{}c/(720z^3)`$ . This results in the normal Casimir force per unit area $`F=U/z=\pi ^2\mathrm{}c/(240z^4)`$. As metals have only a finite conductivity, corresponding corrections can be incorporated in terms of the plasma wavelength $`\lambda _p`$ leading to

$$E_{pp}(z)=\frac{\pi ^2\mathrm{}c}{720z^3}\left[1+\underset{n=1}{\overset{4}{}}c_n\left(\frac{\lambda _p}{2\pi z}\right)^n\right],$$
(1)
where the coefficients are $`c_1=4`$, $`c_2=72/5`$,

$`c_3={\displaystyle \frac{320}{7}}\left(1{\displaystyle \frac{\pi ^2}{210}}\right),c_4={\displaystyle \frac{400}{3}}\left(1{\displaystyle \frac{163\pi ^2}{7350}}\right).`$

Eq.(1) is applicable at the separations $`z\lambda _p`$ with an accuracy of about (1–2)%.

Now let us consider plates covered by the longitudinal corrugations described by the functions

$$z_1=A_1\mathrm{sin}(2\pi x/\mathrm{\Lambda }),z_2=z+A_2\mathrm{sin}(2\pi x/\mathrm{\Lambda }+\phi ),$$
(2)
where $`z`$ is the mean separation distance between the two surfaces, $`\mathrm{\Lambda }`$ is the corrugation period, $`A_1`$ and $`A_2`$ are the corrugation amplitudes.

The normal separation distance between two points of the corrugated surfaces can be represented as

$$z_2z_1=z+b\mathrm{cos}(2\pi x/\mathrm{\Lambda }\alpha ),$$
(3)
where

$`b=b(\phi )=\left(A_1^2+A_2^22A_1A_2\mathrm{cos}\phi \right)^{1/2},`$ (4)
$`\mathrm{tan}\alpha =(A_2\mathrm{cos}\phi A_1)/(A_2\mathrm{sin}\phi ).`$ (5)
For corrugated plates with large corrugation period $`(\mathrm{\Lambda }>z)`$ the Casimir energy density can be found by the additive summation of the results obtained for plane plates (see Eq. (1)) using the fact that all the separation distances $`z_2z_1`$ given by Eq. (3) are equally probable

$$E_{pp}^{cor}(z,\phi )=\frac{1}{\mathrm{\Lambda }}_0^\mathrm{\Lambda }E_{pp}(z_2z_1)𝑑x.$$
(6)
Experimentally it is hard to maintain two parallel plates uniformly separated by distances less than a micron. So one of the plates is replaced by a metal sphere of radius $`R`$ where $`Rz`$. Here a sphere of radius $`R=100\mu `$m imprinted by sinusoidal corrugations was used instead of one of the plates (see Fig. 1). As both $`zR`$ and $`\mathrm{\Lambda }R`$, the normal Casimir force can be calculated by use of the Proximity Force Theorem  as $`F^{nor}(z,\phi )=2\pi RE_{pp}^{cor}(z,\phi )`$. The accuracy of this theorem for our parameters is of order 0.2% (note that the recent result  claiming a worse accuracy for the PFT is applicable only to the pure nonrelativistic regime with separations $`z`$ no larger than several nanometers and also small $`R`$). By integrating the normal force with respect to separation, the energy of a corrugated sphere and a plate is obtained. Then, differentiating with respect to the phase shift, one finds the lateral Casimir force

$$F^{lat}(z,\phi )=\frac{4\pi ^2R}{\mathrm{\Lambda }^2}\frac{}{\phi }_z^{\mathrm{}}𝑑y_0^\mathrm{\Lambda }E_{pp}(z_2z_1)𝑑x.$$
(7)
Here all integrals can be calculated explicitly by the use of representation (3) for the separation distance. The result is

$$F^{lat}(z,\phi )=\frac{\pi ^4R\mathrm{}c}{120z^4}\frac{A_1A_2\mathrm{sin}\phi }{\mathrm{\Lambda }(1\beta ^2)^{5/2}}\left[1+\underset{n=1}{\overset{4}{}}c_{n,x}\left(\frac{\lambda _p}{2\pi z}\right)^n\right],$$
(8)
where $`\beta b(\phi )/z`$ and the coefficients are given by

$`c_{1,x}={\displaystyle \frac{4+\beta ^2}{3(1\beta ^2)}}c_1,c_{2,x}={\displaystyle \frac{5(4+3\beta ^2)}{12(1\beta ^2)^2}}c_2,`$ (9)
$`c_{3,x}={\displaystyle \frac{8+12\beta ^2+\beta ^4}{4(1\beta ^2)^3}}c_3,c_{4,x}={\displaystyle \frac{7(8+20\beta ^2+5\beta ^4)}{24(1\beta ^2)^4}}c_4.`$ (10)
There are also corrections due to the nonzero temperature which can be neglected for the separations reported here . It can be shown that the impact of surface roughness on $`F^{lat}`$ is insignificant.

The experiments are performed using a standard Atomic Force Microscope (AFM) at a pressure below 50 mTorr and at room temperature. A schematic diagram of the experiment is shown in Fig. 1. To implement this experiment a diffraction grating with an uniaxial sinusoidal corrugations of period $`\mathrm{\Lambda }=1.2\mu `$m and an amplitude of 90 nm was used as the template. In order to obtain perfect orientation and phase between the two corrugated surfaces a special procedure was developed, where the in situ pressure imprinting of the corrugations from the plate onto the gold coated sphere was done. This imprinting procedure required special adaptation of the cantilever. First a polystyrene sphere was attached to the cantilever tip with conductive silver epoxy. Next a $`<10\mu `$m thick, $`100200\mu `$m wide and 0.5 mm long piece of freshly cleaved mica is attached to the bottom of the sphere with silver epoxy. Then a second polystyrene sphere of 200$`\pm 4\mu `$m diameter was mounted on the tip of mica with the same silver epoxy. This second sphere is used in the measurement. The cantilever (with mica plate and spheres), corrugated plate and a smooth flat plate (polished sapphire) were then coated with about 400 nm of gold in a thermal evaporator. A small section of the corrugated plate is also coated with 100 nm of aluminum (as $`Al`$ is much harder than gold the imprinting of the corrugations on the sphere is done on this surface). The sphere and the plates are mounted as shown in Fig. 1. As the AFM cantilever is only 320 $`\mu `$m long, the addition of the mica plate and the second sphere is necessary to (i) allow access to the central regions of the corrugated plate (ii) reduce the effect of scattered laser light from the edges of the corrugated plate and (iii) increase the sensitivity of the cantilever. Now the sphere is moved over to the region of the corrugated surface coated with $`Al`$ and the corrugations are imprinted on the gold coating of the surface by pressure using the piezos shown. After this imprinting, the sphere is translated over to the gold coated area of the plate. Extreme care to preserve the parallel orientation of the two corrugations is necessary during this translation, as even a misallignment of 0.1 rad can lead to crossovers of the axis of the two corrugations and destruction of any lateral Casimir force. We used optical beams to monitor the orientation of the cantilever holder with an accuracy of 2$`\times 10^3`$ radians during the translation process. The amplitude of the corrugation on the plate and the imprinted one on the sphere were measured using the AFM, to be $`59\pm 7`$nm and $`8\pm 1`$nm respectively. This measurement was made after completion of the experiment.

The experimental arrangement in Fig. 1 is so chosen to increase the sensitivity for detection of lateral forces and suppress the effect of the normal Casimir force on the cantilever. Thus a lateral force tangential to the corrugated sphere surface would result in the usual bending of the cantilever in response to the force. Whereas a force acting normal to the sphere and corrugated plate (from the normal Casimir force) would lead to the torsional deflection (rotation) of the cantilever. The torsional spring constant of this cantilever $`k_{tor}`$ is much greater than the bending spring constant $`k_{ben}`$, making it much more sensitive to detecting the lateral Casimir force, while simultaneously suppressing the effect of the normal Casimir force. The calibration of the cantilever ( $`k_{tor}`$ and $`k_{ben}`$) and the measurement of the residual potentials between the sphere and plate is done by electrostatic means in a manner similar to that in Refs. . Here the sphere is kept grounded and various voltages are applied to the plates. The normal electrostatic force between the corrugated sphere and plate is given by:

$$F^{el}(z,\phi )=\pi R\epsilon _0\frac{(V_1V_0)^2}{z}\frac{1}{\sqrt{1\beta ^2}},$$
(11)
where as before $`z`$ is the total distance between the zero corrugation levels on the corrugated surfaces. $`V_1`$ are the voltages applied on the corrugated plate and $`V_0`$ is the residual potential on the grounded sphere. If $`V_1`$ is applied to the corrugated plate, the electrostatic force on the sphere leads to the torsional rotation of the cantilever. From the electrostatic force at different $`V_1`$ we measure the torsional spring constant $`k_{tor}=0.138\pm 0.005`$N/m and $`V_0=0.135`$V. Similarly if the sphere is brought close to the bottom smooth plate and voltages $`V_1`$ are applied to the bottom plate (here in Eq. (11) $`A_1=A_2=\beta =0`$ due to the smooth surfaces), the electrostatic force leads to the normal bending of the cantilever. From the force at different voltages $`V_1`$ we again measure the normal spring constant $`k_{ben}=0.0052\pm 0.0001`$N/m. Note that $`k_{tor}k_{ben}`$ is required for isolation and detection of the lateral Casimir force. The piezo extension with applied voltage was calibrated by optical interferometry .

Next the lateral Casimir force is measured. The residual potential difference between the corrugated sphere and plate is compensated by application of voltage $`V_0`$ to the corrugated plate. The corrugated plate is mounted on two piezos that allow independent movement of the plate in the vertical and the horizontal directions with help of a $`x`$-piezo and $`z`$-piezo respectively. Movement in the $`x`$ direction with the $`x`$-piezo is necessary to achieve lateral phase shift $`\phi `$ between the corrugations on the sphere and the plate. Independent movement in the $`z`$ direction is necessary for control of the surface separation between the corrugated sphere and plate. Also, the corrugated plate could be mounted with a small but non-zero tilt away from the vertical. This would lead to changes in surface separation during lateral translation of the plate in the $`x`$ direction. To prevent this a small correction voltage is applied to the $`z`$-piezo synchronous with the lateral translation in $`x`$ direction, to keep the surface separation distance constant.

At the start of lateral force measurement, the sphere is brought close to the corrugated plate and the separation distance is kept fixed. To measure the lateral force $`F^{lat}`$ as a function of the phase $`\phi `$, the corrugated plate is moved in the $`x`$-direction in 0.46 nm steps using the $`x`$-piezo and the lateral force is measured at each step. This is repeated 60 times and the average lateral force at each step is recorded. The average lateral force measured is shown as the solid squares in Fig. 2. The sinusoidal oscillations expected from Eq. (8) as a function of the phase difference between the two corrugations are clearly observed. Also the periodicity of the lateral force oscillation is in agreement with corrugation period of the plate. A sine curve fit to the observed data is shown as the solid line and corresponds to an amplitude of 3.2$`\times 10^{13}`$N. From Eq. (8) this corresponds to a separation distance of $`z=221\pm `$2 nm between the two corrugated surfaces. The mean quadratic error of the average lateral force amplitude is $`\sigma _{\overline{A}}=0.15\times 10^{13}`$N. Systematic error is around 5% leading to the absolute error $`\mathrm{\Delta }_A^{(s)}=0.16\times 10^{13}`$N. Using the value of Student coefficient $`t_{0.95,\mathrm{\hspace{0.17em}60}}=2`$ one obtains for the confidence interval $`\mathrm{\Delta }_A=0.46\times 10^{13}`$N with a 95% confidence probability. The resulting precision of the amplitude measurement is around 15%.

Next the separation between the sphere and corrugated plate is increased by 12 nm with the $`z`$-piezo and the measurement is repeated. The average measured amplitude of lateral force of $`2.6\times 10^{13}`$N corresponds from Eq. (8) to $`z=233\pm `$2 nm consistent with the 12 nm increase in the separation distance. Thus the measured lateral force is in agreement with the complete theory taking into account the conductivity corrections. The separation distance is increased in 12 nm steps and the lateral force is measured two more times. The amplitude of the measured forces was found to be consistent with the corresponding separation distances. In Fig. 3 a log-log plot of the amplitude of measured lateral force as a function of the separation distance is shown as solid squares. Here the separation distance of 221 nm determined above is used for the closest point. The remainder of the points are fixed by the 12 nm step increase in the separation distance. A linear fit to the data yields a slope of 4.1$`\pm `$0.2 consistent with the inverse fourth power $`z`$ dependence of the lateral force expected from Eq. (8) (corrections to this dependence are rather small as the maximal value of $`\beta =0.3`$). Note that a lateral force resulting from electrostatic effects would lead to an inverse second power $`z`$ dependence.

In conclusion, we have made the first measurement of the lateral Casimir force. It was measured between a corrugated plate and corrugated sphere using an AFM. Exact allignment of the corrugations on the sphere and plate was achieved by a special imprinting technique. The measured lateral Casimir force is sinusoidal with a periodicity corresponding to the corrugation period. It was shown to have an inverse fourth power distance dependence which is in agreement with the theory taking into account the conductivity corrections. This demonstration of the lateral Casimir force should open new opportunities for the application of the Casimir effect in microelectromechanical systems.

This work is supported by a National Science Foundation Nanoscale Exploratory Research grant and the National Institute for Standards and Technology through a Precision Measurement Grant. G.L.K. and V.M.M. were also supported by CNPq.

Figure 1: Schematic of experimental setup. For clarity, the size of the corrugations have been exaggerated and the additional mica plate and the first sphere attached to the cantilever are not shown. The $`x`$-piezo and $`z`$-piezo are independent.
Figure 2: The average measured lateral Casimir force as a function of the lateral displacement of the corrugated plate is shown as solid squares. The solid line is the best fit sine curve to the data leading to a lateral force amplitude of $`3.2\times 10^{13}`$N.
Figure 3: log-log plot of the measured lateral Casimir force amplitude as a function of the surface separation distance is shown as solid squares. The slope of the straight line fit is 4.1$`\pm `$2.


text:
# Is the mean-field approximation so bad? A simple generalization yelding realistic critical indices for 3D Ising-class systems

## Abstract

Modification of the renormalization-group approach, invoking Stratonovich transformation at each step, is proposed to describe phase transitions in 3D Ising-class systems. The proposed method is closely related to the mean-field approximation. The low-order scheme works well for a wide thermal range, is consistent with a scaling hypothesis and predicts very reasonable values of critical indices.

It is well understood since 1960’s that an essential increase of fluctuations accompanies a second-order phase transition. Exact results for 2D Ising model Ising2D  showed that Landau phenomenology Landau  was too crude to describe the vicinity of the transition point. Ginzburg demonstrated that it is correlations what breaks Landau approach down Ginzburg . He established the parameter which determines the thermal range for the so-called fluctuation region. According to the scaling hypothesis PP , the properties of a system in the fluctuation region are determined by the single quantity $`r_c`$, being the correlation length for fluctuations. This parameter shows a power-law dependence on the external field $`h`$ and on the thermal interval $`t`$ from the transition: $`r_ch^\mu `$, $`r_ct^\nu `$. Values of $`\mu `$ and $`\nu `$ are determined by the universal characteristics, like the dimensionality $`D`$, the number of components of the order parameter, and the type of the interaction in the system (long-range or short-range). Other critical indices can be expressed via $`\mu `$ and $`\nu `$. For instance for the scalar order parameter $`\overline{x}`$ one obtains $`\overline{x}t^\beta ,\beta =\nu D\nu /\mu `$.

The general framework to obtain critical indices is the renormalization-group (RG) approach Wilson ; RevRG ; ERGE . The RG transformation consists of the two stages. First, one integrates out all $`x_k`$ with wave vectors $`k>k_0/s`$, where $`k_0`$ is a cut-off and $`s>1`$ is a parameter. The main assumption of the method is that the potential energy of thus renormalized system can be truncated to certain simplified form, for example to the form of $`\varphi ^4`$ model. The second step is a scale transformation. It is aimed to restore the cut-off for $`k`$ back to $`k_0`$ and preserve the original dispersion law for the long-wavelength modes.

The RG transformation should be applied recursively; its unstable fixed point corresponds to a phase transition. Lyapunov exponents for the vicinity of the fixed point determine the critical indices. The scaling hypothesis requires exactly two of the Lyapunov exponents to be positive, they correspond to the two independent critical indices.

With respect to the value of $`s`$, two modifications of the RG scheme can be distinguished. Wilson’s $`\epsilon `$ expansion Wilson ; RevRG  deals with $`s1`$. The method is formally valid near the upper critical dimension $`D_{up}`$ (for the Ising-class classical systems $`D_{up}=4`$).The ”most divergent“ diagram series can be pointed out in this case. This gives the asymptotic expansion in terms of $`\epsilon D_{up}D`$. In the continuous version of the scheme (so-called exact RG approach) ERGE , $`s`$ is infinitesimally close to unity. Certain decoupling (truncation) should be introduced in this case to solve the integro-differential equation for RG flow. An accuracy is then determined by a form of the decoupling used.

All the approaches mentioned above are very well-developed. They allow to calculate critical indices and universal combinations of amplitudes with a high accuracy. Progress is also achieved in the calculation of non-universal quantities; particularly a crossover between Gaussian (Landau-like) and critical behaviour is analyzed crossover . However RG is applicable only for long-wavelength excitations since they can be described as a fluctuating field. Therefore, in a practical calculation, short-wavelength excitations should be first integrated out by means of some perturbation expansion, and then RG should be applied. The short- and the long-wavelength excitations should be separated ad hoc. Such a two-step procedure is not very elegant. At least from the methodological point of view it is desired to develop a scheme which gives all the information in a uniform manner.

It is interesting to note that calculations are much simpler outside the critical region. In fact, even the mean-field approximation (MFA) is usually sufficient for the qualitative analysis of this region. In MFA, the non-local interaction in the system is replaced by an effect of a self-consistent average field. The MFA predictions for the transition temperature and thermal behaviour of the order parameter in 3D case are quite reasonable MFA . On the other hand MFA gives Landau set of indices and is not very accurate. Note that Landau approach itself has a mean-field nature.

In this paper, we present a simple generalization of MFA, performing well for the 3D Ising-class systems both nearby and far from the transition. Particularly the critical behaviour is reproduced correctly: the values of indices are predicted with a few-percent accuracy. Essentially, the present approach consists in the alternating application of the Stratonovich transformation and the most primitive version of the renormalization-group transformation at a finite value of $`s`$.

We consider a classical scalar field with an anharmonic local potential and a short-range harmonic interaction at temperature $`T`$. The partition function of this system is

$`Z`$ $`=`$ $`{\displaystyle [Dx]e^{V/T}};`$ (1)
$`V`$ $`=`$ $`{\displaystyle U(x_r+\overline{x})𝑑𝐫}+{\displaystyle \underset{k}{}}{\displaystyle \frac{\mathrm{\Omega }_k|x_k|^2}{2}}.`$
Here $`U(x+\overline{x})`$ is an even on-site potential, $`x_k`$ is a Fourier transform of the displacement $`x_r`$ from the average position $`\overline{x}`$; $`\overline{x}0`$ below the transition temperature, and the average of $`x_r`$ equals zero by definition. To make the definitions unambiguous, we put additional conditions $`U(\overline{x})=0;\mathrm{\Omega }_{k=0}=0`$. The wave vector $`k`$ does not exceed the unititary cut-off: $`k^21`$. The case of a short-range interaction corresponds to the asymptotic behaviour

$$\mathrm{\Omega }_k=\omega k^2,k0.$$
(2)
The simplest way to integrate out the modes is to skip $`x_k`$ with $`k>1/s`$. This corresponds to the zero-th approximation in the nonlinearity of $`U`$. The subsequent scaling results in the RG transformation $`ksk,Us^3U,xx/\sqrt{s}`$. The defect of thus applied RG method comes from the the increase of correlations (non-linear part of $`U`$ grows). Therefore, after several iterative transformations, the RG procedure is not valid any more. Although the formal analysis of the fixed point $`U(x)=0`$ can be done, it wrongly gives Landau-like critical indices.

To avoid the problem higher-order approximations and the truncation of the potential are used in known schemes described above (a very clear explanation is given in fermi ). Here we use another procedure: the zero-th approximation is used, but the change of variables is made at each RG step. The change of variables reduces correlations and is aimed to compensate its increase at RG transformation.

At the first stage, we integrate out (omit) all the modes with $`\mathrm{\Omega }_k>\overline{\mathrm{\Omega }}`$, where $`\overline{\mathrm{\Omega }}=\frac{3}{4\pi }\mathrm{\Omega }_kd^3k`$ is an average of $`\mathrm{\Omega }_k`$. For the rest of the modes, the Stratonovich identity is utilized:

$$\mathrm{exp}\left(\frac{(\mathrm{\Omega }_k\overline{\mathrm{\Omega }})|x_k|^2}{T}\right)=\frac{\overline{\mathrm{\Omega }}\mathrm{\Omega }_k}{\pi \overline{\mathrm{\Omega }}^2}\mathrm{exp}\left(\frac{1}{T}\left(2\overline{\mathrm{\Omega }}\mathrm{Re}(x_kf_k^{})+\frac{\overline{\mathrm{\Omega }}\mathrm{\Omega }_k|f_k|^2}{\overline{\mathrm{\Omega }}\mathrm{\Omega }_k}\right)\right)𝑑f_k.$$
(3)
Here the integration over the complex variable $`f_k=\mathrm{}(f_k)+i\mathrm{}(f_k)`$ denotes the integration over the complex plane: $`𝑑f_k_{\mathrm{}}^{\mathrm{}}d\mathrm{}(f_k)_{\mathrm{}}^{\mathrm{}}d\mathrm{}(f_k)`$.

After (3) is substituted into (1), expressions $`_k\mathrm{Re}(x_kf_k^{})=x_rf_r𝑑𝐫`$ and $`_k\overline{\mathrm{\Omega }}|x_k|^2=\overline{\mathrm{\Omega }}x_r^2`$ can be collected in the exponent. Therefore

$$Z[Dx][Df]\mathrm{exp}\left(\frac{1}{T}\left(\left(U(x_r+\overline{x})\mathrm{\Omega }_kf_rx_r+\frac{\overline{\mathrm{\Omega }}x_r^2}{2}\right)𝑑𝐫+\underset{k}{}\frac{\overline{\mathrm{\Omega }}\mathrm{\Omega }_k|f_k|^2}{2(\overline{\mathrm{\Omega }}\mathrm{\Omega }_k)}\right)\right),$$
(4)
where $`f_r`$ is the inverse Fourier transform of $`f_k`$.

Now, we can integrate over $`x_r`$. Let us introduce the function $`F(f)`$ accordingly to the equation:

$$\mathrm{exp}\left(\frac{F(f+\overline{x})F_0}{T}\right)=\mathrm{exp}\left(\frac{1}{T}\left(U(x+\overline{x})+\frac{\overline{\mathrm{\Omega }}(xf)^2}{2}\right)\right)𝑑x,$$
(5)
where $`F_0`$ is a constant, defined from the condition $`F(\overline{x})=0`$.

With thus defined $`F`$, the partition function takes the form

$`Z`$ $``$ $`{\displaystyle [Df]e^{W/T}};`$ (6)
$`W`$ $`=`$ $`{\displaystyle F(f_r+\overline{x})𝑑𝐫}+{\displaystyle \underset{k}{}}{\displaystyle \frac{\overline{\mathrm{\Omega }}\mathrm{\Omega }_k|f_k|^2}{2(\overline{\mathrm{\Omega }}\mathrm{\Omega }_k)}}.`$
Equation (6) is formally very similar to Eq.(1). Function $`F`$ is even; the dispersion at $`k0`$ coincides that of Eq.(1). One can also guess from Eq.(3) that as the average of $`x_r`$ equals zero, the average of $`f_r`$ also vanishes. Indeed, if this average would not be zero, $`f_{k=0}`$ would fluctuate around a macroscopically large average value, and (3) would not be fulfilled.

We argue here that variables $`f_k`$ are much less correlated, than $`x_k`$. Let us consider the parabolic approximation for $`F`$: $`F(\overline{x}+f)\frac{F(\overline{x})}{f}f+\frac{^2F(\overline{x})}{2f^2}f^2`$. Variables $`f_k`$ are uncorrelated in this approximation. Obviously $`<f>`$ vanishes only if $`\frac{F(\overline{x})}{f}=0`$, that gives

$$\frac{\overline{\mathrm{\Omega }}\overline{x}}{T}=\frac{x\mathrm{exp}\left(T^1\left(U(x+\overline{x})+\frac{1}{2}\overline{\mathrm{\Omega }}x^2\right)\right)𝑑x}{\mathrm{exp}\left(T^1\left(U(x+\overline{x})+\frac{1}{2}\overline{\mathrm{\Omega }}x^2\right)\right)𝑑x}.$$
(7)
This is exactly the mean-field equation of state for the system (1). So, even the complete neglect of correlations in $`f_k`$ still gives something reasonable. On the other hand, the parabolic approximation for $`U`$ directly in Eq.(1) gives non-sense: since $`U(x)`$ does not depend on temperature, the system does not show a transition at all in this approximation. Such a comparison convinces us that correlations in $`f_k`$ are less essential, than in $`x_k`$. As it was pointed above we expect that this can compensate an increase of correlations in RG transformation.

Now, we finish the RG transformation, operating with Eq.(6). All $`f_k`$ with $`k>1/s`$ are to be integrated out (omitted in our approximation). After the scale transformation $`ksk,fx/\sqrt{s}`$ is done, the potential energy of the system takes the form

$$U^{(1)}(x_r+\overline{x}^{(1)})𝑑𝐫+\underset{k}{}\frac{\mathrm{\Omega }_k^{(1)}|x_k|^2}{2},$$
(8)
where

$$\mathrm{\Omega }_k^{(1)}=\frac{s^2\overline{\mathrm{\Omega }}\mathrm{\Omega }_k}{\overline{\mathrm{\Omega }}\mathrm{\Omega }_{k/s}},$$
(9)
$$U^{(1)}(x+\overline{x}^{(1)})=s^3F(x/\sqrt{s}+\overline{x}),\overline{x}^{(1)}=\overline{x}\sqrt{s}.$$
(10)
Such a transformation should be applied recursively, like in the standard RG approach. The analysis can be easily performed numerically.

The sequence $`\mathrm{\Omega }_k,\mathrm{\Omega }_k^{(1)},\mathrm{\Omega }_k^{(2)},\mathrm{}`$ is independent of $`U`$ and $`T`$. At every step a new $`\overline{\mathrm{\Omega }}^{(n)}`$ are to be calculated. After several transformations $`\overline{\mathrm{\Omega }}^{(n)}`$ and $`\mathrm{\Omega }_k^{(n)}`$ converge to certain limit. Calculation from formula (9) gives

$$\mathrm{\Omega }_k^{(\mathrm{})}=\frac{\omega \overline{\mathrm{\Omega }}^{(\mathrm{})}(s^21)k^2}{\overline{\mathrm{\Omega }}^{(\mathrm{})}(s^21)\omega k^2};$$
(11)
$`\overline{\mathrm{\Omega }}^{(\mathrm{})}`$ is implicitly defined from the equation $`\overline{\mathrm{\Omega }}^{(\mathrm{})}=\frac{3}{4\pi }\mathrm{\Omega }_k^{(\mathrm{})}d^3k`$.

Now, let us consider the properties of the sequence $`U,U^{(1)},U^{(2)},\mathrm{}`$ at $`\overline{x}=0`$. Transformation $`U^{(n)}U^{(n+1)}`$ has a non-trivial unstable fixed point $`U_f(x)`$. This function is shown in the upper panel of Fig.1. Consider small deviations from this fixed point: $`U^{(n)}(x)=U_f(x)+u_n(x)`$. The linearization of the transformation (10) gives

$$u_{n+1}(y)=u^0+s^3u_n(x)A(x,y)𝑑x;$$
(12)
$$A=\mathrm{exp}\left(\frac{U_f(x)s^3U_f(y)+F+\frac{\overline{\mathrm{\Omega }}}{2}(xy/\sqrt{s})^2}{T}\right).$$

where $`u^0`$ stands to fulfill the condition $`u_{n+1}(0)=0`$. The analysis of this operator shows that among its eigenvalues only two exceed unity. Denote them $`s^{\lambda _e}`$ and $`s^{\lambda _o}`$; they correspond to an even and odd eigenvector, respectively. Thus the theory is consistent with the scaling hypothesis, as it was mentioned in an introductory part.

Exponents $`\lambda _e`$ and $`\lambda _o`$ are the inverse of the critical indices $`\nu `$ and $`\mu `$, respectively RevRG ; fermi . In the ideal case, critical indices should be independent of the particular $`U(x)`$ and $`\mathrm{\Omega }_k`$, as well as of the value of the factor $`s`$. Indeed, after a proper re-scaling of the units, the values of $`T`$ and $`\overline{\mathrm{\Omega }}`$ drop out from Eq.(10). Therefore in our model $`\nu `$ and $`\mu `$ can depend only on $`s`$. Further, it turns out that $`\mu =2/5`$ at any $`s`$, as it should be for the local-potential approximation ERGE . What is about $`\nu `$, it depends on $`s`$ very weakly. The numerical result for $`\nu (s)`$ is shown in the lower panel of Fig.1.

On the other hand it is important to point out that the crude universal properties of the model enter the transformation (10) essentially. The behaviour $`\mathrm{\Omega }_{k0}k^2`$ corresponds to a short-range range interaction. Such behaviour and the dimensionality of the model determine the scale transformation, i.e. powers of $`s`$ in Eq.(10). It is also directly reflected in the formula that the order parameter is a scalar.

For the more detailed analysis it is reasonable to consider an extreme value of $`\nu (s)`$, because in this case small variations of $`s`$ do not affect $`\nu `$. This gives $`s2.5,\nu 0.652`$. The obtained values of indices are collected in Table 1. An agreement with the numerical results for the 3D Ising model is even better than one could expect. Since the scaling hypothesis holds in the model, the values of all other critical indices can be determined and should have a small error bar. As an example, the calculated value of $`\beta `$ is presented in Table 1.

The advantage of the present approach is that it allows to estimate not only critical indices, but also the (non-universal) thermodynamic quantities in a wide thermal range. Here we calculate the dependence of the order parameter $`\overline{x}`$ on temperature for the 3D Ising model. The model is defined as a discrete cubic lattice with $`U(x)=\delta (x1)+\delta (x+1)`$ and the nearest-neighbors interaction $`\mathrm{\Omega }_k=3\mathrm{cos}\pi k_x\mathrm{cos}\pi k_y\mathrm{cos}\pi k_z`$. The initial Brillouin zone is a cube $`1<k_{x,y,z}<1`$; at the first step the value $`\overline{\mathrm{\Omega }}=3`$ is used and all the modes with $`|k|>1/s`$ are integrated. The procedure for further steps is exactly as described above. To simplify the calculation, we put $`\overline{\mathrm{\Omega }}^{(n)}=\overline{\mathrm{\Omega }}^{(\mathrm{})}`$ at these steps; this simplification almost does not change the result obtained. To determine the value of $`\overline{x}`$ we use the condition

$$\frac{}{x}U_n(\overline{x}_n)=0,n\mathrm{}$$
(13)
(compare with the derivation of formula (7)). Numerically, it is enough to calculate $`U_n`$ up to $`n8`$. The result for $`\overline{x}^2(T)`$ is presented in Fig.2. The mean-field and experimental (Monte-Carlo) data are sketched in the same figure for comparison. Far from the transition $`U^{(1)}`$ is almost parabolic, therefore the theory passes into MFA. The critical behaviour occurs at the transition point. Note that Ising model is very non-linear, but the scheme still performs good.

We conclude that the presented low-order scheme describes correctly all the physics of the second-order phase transitions in 3D systems with a short-range interaction and a scalar order parameter. It gives a unified description both of the fluctuation region and of the system far from the transition point. The method does not pretend to be very accurate; its advantage consists in simplicity. In principle, one can improve the accuracy by a higher-order procedure for the integration in RG transformation. An applicability of the method for another universality classes is the subject of further studies.

The work was supported by RFFI foundation (grant 00-02-16253) and by the ”Russian Scientific Schools” program (grant 96-1596476).