text:
# Nonlinear Stochastic Resonance with subthreshold rectangular pulses

\]http://numerix.us.es

## Abstract

We analyze the phenomenon of nonlinear stochastic resonance (SR) in noisy bistable systems driven by pulsed time periodic forces. The driving force contains, within each period, two pulses of equal constant amplitude and duration but opposite signs. Each pulse starts every half-period and its duration is varied. For subthreshold amplitudes, we study the dependence of the output signal-to-noise ratio (SNR) and the SR gain on the noise strength and the relative duration of the pulses. We find that the SR gains can reach values larger than unity, with maximum values showing a nonmonotonic dependence on the duration of the pulses.

In recent work Casgomprl ; Casgompre , we have carried out a detailed analytical and numerical study of the nonlinear response of a noisy bistable system subject to a subthreshold time periodic driving force. The force considered in those works is such that it remains constant within the duration of each half period while switching its sign every half period. We have focused our attention on the analysis of the dependence of the output signal-to-noise ratio and the corresponding SR gain on the noise strength. The analytical study was based on a two-state approximation amenable to exact treatment Casgomprl . We showed that, for subthreshold input signals of sufficiently long periods, the phenomenon of SR can be accompanied by SR gains larger than unity. This is a genuine characterization of nonlinearity, as SR gains larger than unity are strictly forbidden within a linear response description dykluc95 ; dewbia95 ; casgom03 . The analytical results were corroborated by numerical simulations Casgompre .

In the last few years, Gingl and collaborators ginmak00 ; ginmak01 ; ginmak02  have carried out analog simulations of noisy bistable systems subject to driving forces of the type given in Eq. (1)

$$F(t)=\{\begin{array}{cc}\hfill A:& 0t<t_c\hfill \\ \hfill A:& \frac{T}{2}t<\frac{T}{2}+t_c\hfill \\ \hfill 0:& \text{otherwise}\hfill \end{array}$$
(1)
It is convenient to introduce the parameter $`r=2t_c/T`$, measuring the fraction of a period during which this driving force has a nonvanishing value (the parameter $`r`$ in the present paper corresponds exactly to what Gingl and coworkers term “duty cyle”). In ginmak00 ; ginmak01 ; ginmak02 , the SNR and the SR gain for subthreshold amplitude input signals with $`r0.3`$ were studied. These authors find SR gains larger than unity, and, also, that increasing the $`r`$ value lowers the SNR gain. They rationalize their observations by noting that the input SNR increases as $`r`$ increases, while the output SNR is less sensitive to the value of $`r`$. The case studied by us in Casgomprl ; Casgompre  corresponds to the largest possible value of the parameter $`r`$, namely, $`r=1`$. It seems therefore of interest to extend our analysis to input signals with $`r<1`$ in order to compare with the predictions of Gingl et al.

Let us consider a system characterized by a single degree of freedom, $`x`$, whose dynamics (in dimensionless units) is governed by the Langevin equation

$$\dot{x}(t)=U^{}[x(t),t]+\xi (t),$$
(2)
where $`\xi (t)`$ is a Gaussian white noise of zero mean with $`\xi (t)\xi (s)=2D\delta (ts)`$, and $`U^{}(x,t)`$ represents the force stemming from the time-dependent, archetype bistable potential

$$U(x,t)=\frac{x^4}{4}\frac{x^2}{2}F(t)x,$$
(3)
with $`F(t)`$ given by Eq. (1). The one-time correlation function is defined as

$$C(\tau )=\frac{1}{T}_0^Tdtx(t+\tau )x(t)_{\mathrm{}}.$$
(4)
It can be written exactly as the sum of two contributions: a coherent part, $`C_{coh}(\tau )`$, which is periodic in $`\tau `$ with period $`T`$, and an incoherent part, $`C_{incoh}(\tau )`$, which decays to $`0`$ for large values of $`\tau `$ and reflects the correlation of the output fluctuations about its average (the noisy part of the output). The coherent part $`C_{coh}(\tau )`$ is given by RMP ; PRA91

$$C_{coh}(\tau )=\frac{1}{T}_0^Tdtx(t+\tau )_{\mathrm{}}x(t)_{\mathrm{}},$$
(5)
and $`C_{incoh}(\tau )`$ is obtained from the difference of Eq. (4) and Eq. (5). In the expressions above, the susbcript indicates that the averages are to be evaluated in the limit $`t\mathrm{}`$.

According to McNamara and Wiesenfeld McNWie89 , the output SNR, $`R_{out}`$, is defined in terms of the Fourier transform of the coherent and incoherent parts of $`C(\tau )`$ as

$$R_{out}=\frac{lim_{ϵ0+}_{\mathrm{\Omega }ϵ}^{\mathrm{\Omega }+ϵ}d\omega \stackrel{~}{C}(\omega )}{\stackrel{~}{C}_{incoh}(\mathrm{\Omega })},$$
(6)
where $`\stackrel{~}{H}(\omega )`$ denotes the Fourier cosine transform of $`H(\tau )`$, i.e., $`\stackrel{~}{H}(\omega )=2/\pi _0^{\mathrm{}}d\tau H(\tau )\mathrm{cos}(\omega \tau )`$. Note that this definition of the output SNR differs by a factor $`2`$, stemming from the same contribution at $`\omega =\mathrm{\Omega }`$, from the definitions used in earlier works RMP ; PRA91 . The periodicity of the coherent part gives rise to delta peaks in the spectrum. Thus, the only contribution to the numerator in Eq. (6) stems from the coherent part of the correlation function. The output SNR can then be expressed as

$$R_{out}=\frac{Q_u}{Q_l},$$
(7)
where

$$Q_u=\frac{2}{T}_0^Td\tau C_{coh}(\tau )\mathrm{cos}(\mathrm{\Omega }\tau ),$$
(8)
and

$$Q_l=\frac{2}{\pi }_0^{\mathrm{}}d\tau C_{incoh}(\tau )\mathrm{cos}(\mathrm{\Omega }\tau ).$$
(9)
A nonmonotonic behavior of the SNR with the strength of the noise is a signature of the phenomenon of SR.

The signal-to-noise ratio for an input signal $`F(t)+\xi (t)`$ is given by

$$R_{inp}=\frac{\frac{1}{2}(f_1^2+g_1^2)}{\frac{2}{\pi }D}.$$
(10)
where

$$f_1=\frac{2A}{\pi }\mathrm{sin}\mathrm{\Omega }t_c$$
(11)
and

$$g_1=\frac{2A}{\pi }\left(1\mathrm{cos}\mathrm{\Omega }t_c\right)$$
(12)
It is then clear that for fixed values of $`A`$ and $`D`$, $`R_{inp}`$ increases with $`r`$ as pointed out in Ref. ginmak02 . Also, for given $`A`$ and $`r`$, $`R_{inp}`$ decreases as $`D`$ increases.

The SR gain, $`G`$, is defined as the ratio of the SNR of the output over that of the input; namely,

$$G=\frac{R_{out}}{R_{inp}}.$$
(13)
SR gain values larger than 1 have been obtained in driven nondynamical systems ginmak00 , in stochastic resonators with static nonlinearities driven by square pulses chapeau , or in noisy bistable systems driven by superthreshold input sinusoidal frequencies saga . The existence of SR gains with values larger than 1 indicates a truly nonlinear SR.

Although the two-state approximation introduced in Ref. Casgomprl  can, in principle, be extended to analyze systems driven by input signals with $`r<1`$, the analytical expressions obtained are too cumbersome to be of practical value. Thus, in the present work, we rely on the numerical treatment of the Langevin equation, Eq. (2), following the procedure detailed in Ref. casgom03 .

In Fig. 1, we depict the behavior of $`R_{out}`$ with the noise strength $`D`$, for input signals of the type given in Eq. (1), with subthreshold amplitude $`A=0.35`$, fundamental frequency $`\mathrm{\Omega }=0.0024`$ and $`r=0.1,0.4,0.7,0.95,0.98,1`$. For all values of $`r`$, the characteristic nonmonotonic behavior of the SNR with $`D`$ is obtained. As $`r`$ increases, the maximum value of $`R_{out}`$ increases. Namely, the longer the potential remains asymmetric during each half-cycle, the larger the maximum height in the SNR is. Therefore, $`R_{out}`$ is quite sensitive to the duration of the pulses within each half-cycle.

In Fig. 2 the behavior of $`G`$ with the noise strength $`D`$ is depicted for the same values of $`r`$ as in the previous figure. For all values of $`r`$ considered, there exists a range of noise values such that $`G`$ is larger than unity. The peak value of the SR gain has a nonmonotonic behavior with $`r`$. At the lowest value of $`r`$ considered, the SR gain has a rather large peak value. Then, as $`r`$ is increased, the peak of the SR gain decreases, in agreement with the observations in ginmak01 ; ginmak02 . As the duration of the pulses gets larger so that $`r`$ is closer to $`1`$, the tendency of the SR gain maximum reverses and a considerable increase in the maximum is observed.

The just mentioned features can be rationalized by noting that the SR gain depends on $`R_{out}`$ and on $`R_{in}`$. As indicated above, for fixed values of the noise strength, $`D`$, and amplitude, $`A`$, the input $`R_{in}`$ always increases with $`r`$. Also, for given values of $`A`$ and $`r`$, $`R_{in}`$ decreases monotonically with $`D`$. On the other hand, the results depicted in Fig. 1 indicate that there are two main effects on the location of the maximum of $`R_{out}`$ as $`r`$ increases. First, as noted before, the maximum height increases as $`r`$ increases. Second, the maxima appear at increasingly large values of $`D`$ as the duration of the pulses increases. This second effect manifests itself clearly for pulses of sufficiently long duration, namely for $`r`$ values larger than $`r0.9`$, while it is almost unnoticeable for smaller bvalues of $`r`$. For low values of $`r`$ (let us say $`r=0.1`$), even though the peak of $`R_{out}`$ is the smallest one appearing in Fig. 1, the corresponding value of $`R_{in}`$ is so small (due to the smallness of $`r`$) that the SR gain reaches the large values depicted in Fig. 2. As $`r`$ increases, the height of the $`R_{out}`$ maximum also increases, but appearing at an approximately constant value of the noise strength. Thus, the increase of $`R_{in}`$ with $`r`$ counterbalances the increase of $`R_{out}`$ in such a way that the SR gain decreases. Finally, for long duration pulses, the shift to the right of the $`R_{out}`$ maximum and the large increase in its height are the cause of the increase in the maximum gain observed in Fig. 2.

In conclusion, we observe that the behavior of the SR gain for pulses of relatively short duration is basically a consequence of the behavior of $`R_{in}`$, in agreement with the observation of Gingl et al. ginmak00 ; ginmak01 ; ginmak02 . On the other hand, for $`r`$ close to 1, the behavior of the SR gain is dominated by the output SNR.

###### Acknowledgements.

We acknowledge the support of the Dirección General de Enseñanza Superior of Spain (BFM2002-03822) and the Junta de Andalucía.


format: text

text:
# Strong Gravitational Lensing and Dark Energy Complementarity

## I Introduction

Dark energy poses a fundamental challenge to our understanding of the universe. The acceleration of the cosmic expansion discovered through the Type Ia supernovae distance-redshift relation can be interpreted in terms of a new component of the energy density possessing a substantially negative pressure. Further observations from the cosmic microwave background radiation indicate the spatial geometry of the universe is flat; in combination with the supernova measurements this implies that the unknown “dark energy” comprises roughly 70% of the total density, in concordance with large scale structure data indicating that matter contributes approximately 30% of the critical density.

Such a weight of dark energy leads to its dominance of the expansion, causing acceleration, restricting the growth of large scale structure, and holding the key to the fate of the universe. But apart from its rough magnitude, its nature is almost unknown – whether it arises from the physics of the high energy vacuum, a scalar field, extra dimensional or “beyond Einstein” gravitational effects, etc. One way to gain clues to the underlying physical mechanism is to characterize the behavior of dark energy in terms of its equation of state ratio as a function of redshift (essentially a time parameter), $`w(z)`$. This is conventionally interpreted in terms of the ratio of the dark energy pressure to energy density, but can also be used as an effective parameterization to treat generalizations of the cosmology framework of the Friedmann equations of expansion linjen  or simply in terms of the expansion rate and acceleration itself lingrav .

To obtain the clearest focus on the class of physics responsible for the accelerating universe, we seek the tightest constraints on the equation of state. This must come not merely from high statistical precision of a probe, but from robust control of systematic uncertainties. The greatest accuracy and confidence in the measurements will come from independent crosschecks and complementarity between methods of probing the cosmology. Many studies have considered such complementarity between probes (e.g. fhlt ; linap ; seo ; linbo ; wl ; cluster ), with a promising future for next generation surveys carrying out such measurements. However, all the methods of supernova distances, cosmic microwave background power spectrum, weak gravitational lensing, cluster counts, and baryon oscillations possess a similar fundamental dependence on the equation of state through the Hubble parameter, or expansion rate.

This paper investigates whether a truly complementary probe exists that has nearly orthogonal dependence to the previous ones in the plane of the equation of state parameters of value today, $`w_0`$, vs. time variation, $`w_a`$. Such a probe, if practical, would offer a valuable contribution to uncovering fundamental physics and deserve further consideration among next generation experiments. Section II considers the characteristics such a method would possess and identifies two promising candidates related to strong gravitational lensing. In §III we analyze the sensitivity of these probes to the cosmological parameters and the constraints and complementarity they offer. In the conclusion we summarize the prospects for strong lensing as a cosmological probe and discuss some issues regarding surveys and systematic uncertainties.

## II Complementarity in the Equation of State Plane

Astronomical observations involving distances and volumes all follow from the metric in a simple, kinematic way weinberg . For convenience, the dependence can be written in terms of the conformal distance interval

$$d\eta =dt/a(t)=dz/H(z),$$
(1)
where $`a(t)=(1+z)^1`$ is the scale factor as a function of proper time, equivalently parametrized in terms of redshift $`z`$, and $`H=\dot{a}/a`$ is the Hubble parameter. So in observing some source at redshift $`z`$ the cosmological information is carried by the quantity $`H`$.

Varieties of distance observations – angular diameter, luminosity, proper motion – all contain the same information. One exception is the parallax distance (cf. fpoc , §3.2) which involves not only $`H`$ but the spatial curvature $`k`$ as well; however this is not practical for cosmology and CMB measurements strongly indicate a flat universe. Another subtlety involves exotic models that break the thermodynamic, or reciprocity, relation between the various distances thermo ; tolman ; bassett  by violation of Liouville’s theorem, but this arises from photon properties and not cosmology.

Volume elements, and hence numbers of sources, are built up out of distances and so similarly involve the Hubble parameter. From the Friedmann equation the relation between the expansion rate and the dark energy equation of state is

$$H^2/H_0^2=\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_{DE}e^{3{\scriptscriptstyle d\mathrm{ln}(1+z)[1+w(z)]}},$$
(2)
where $`H_0`$ is the present value of the Hubble parameter (the Hubble constant), and for a flat universe $`\mathrm{\Omega }_{DE}=1\mathrm{\Omega }_m`$. As discussed in linjen  this can be written more generally as defining an effective equation of state

$$w(z)=1+\frac{1}{3}\frac{d\mathrm{ln}\delta H^2/H_0^2}{d\mathrm{ln}(1+z)},$$
(3)
where $`\delta H^2`$ encodes our ignorance of the right hand side of eq. 2 after the first, matter density term.

Making $`w(z)`$ more positive (holding $`H_0`$ and $`\mathrm{\Omega }_m`$ fixed) increases the expansion rate in the past, and so decreases the acceleration (since the rate today is fixed). Distances will be less, sources appear brighter, etc. If we try to characterize the nature of the dark energy in the simplest, nontrivial way, through the value of its equation of state today and its increasing or decreasing positivity into the past, we see that cosmological measurements, depending only on $`H`$, don’t really care in a gross sense how $`w(z)`$ was more positive – either increasing the present value $`w_0`$ or increasing the time variation $`w_adw/d\mathrm{ln}(1+z)`$ does the trick. Thus, for all such probes the quantities $`w_0`$, $`w_a`$ will be anticorrelated; an increase in one can be (at least partially) offset by a decrease in the other. If we plot constraints from astronomical data in the $`w_0w_a`$ plane (marginalizing over any other parameters), the probability contours will show a degeneracy direction tilted counterclockwise from the vertical. This is well known and illustrated for a number of different probes in hutcoo ; note also that a “bare” $`H`$ of course has the same dependence linbo .

However, this has the fundamental implication that complementarity between methods can only be partial in this physically key parameter plane; different data sets cannot be orthogonal. This restricts the ability to constrain the equation of state and understand the nature of dark energy. What would be valuable is construction of a cosmological probe whose sensitivity lies clockwise of vertical, into the unexplored half-plane of positive correlation between $`w_0`$ and $`w_a`$. Note that this argument does not depend on the specifics of how $`w_0`$, $`w_a`$ are defined, but for concreteness we adopt the parametrization

$$w(z)=w_0+w_a(1a)=w_0+w_az/(1+z),$$
(4)
successful in fitting smoothly varying dark energy models linprl ; lin10217 . A characteristic time variation of the equation of state is given by $`w^{}dw/d\mathrm{ln}(1+z)|_{z=1}=w_a/2`$.

Figure 1 illustrates this anticorrelation for the distance-redshift probe (dashed ellipse). A hypothetical probe that involves positive correlation between the equation of state parameters (dotted ellipse) can give nearly orthogonal constraints and so the joint probability contour (solid ellipse) between the two types of probe tightly characterizes the equation of state. Furthermore, since reduced data quality, e.g. from systematic uncertainties, often lengthen the confidence contour along the degeneracy direction (long dashed ellipse), an orthogonal probe immunizes against systematic errors. However, degeneracies with other cosmological parameters interfere with both these desiderata, as discussed in §III.

The growth of large scale structure, or equivalently the evolution of gravitational potentials, and hence, say, the number of galaxy clusters of a certain mass, might seem to offer a different dependence than distances. But the growth equation also has $`H`$ as its predominant ingredient and indeed in hutcoo  we see that the growth factor still possesses the anticorrelation. An exception might occur for models where the dark energy (or alternative gravity theory) is inhomogeneous itself and can act as a source to the matter density perturbations, but this is not expected to occur for scalar field dark energy on subhorizon scales caldwell ; dave .

If simple measures of distances or products of distances all give anticorrelation, yet they involve differing levels of sensitivity. This provides the hope that opposing these quantities through ratios might break the pattern of degeneracy. Indeed this was found in linap  on consideration of the cosmic shear test111A note on the name: this probe involves the observed shearing of a sphere due to the properties of spacetime itself. Despite a perfectly isotropic space, a sphere will appear sheared because the spacetime is not isotropic – directions along and perpendicular to the null geodesics differ. This shear is due to cosmic properties, unlike shears in gravitational lensing due to anisotropies not in spacetime but in space, though the latter is unfortunately sometimes called cosmic shear. originated by alpac . Here the key quantity was $`H(z)d_A(z)H(z)𝑑z/H(z)`$, where $`d_A`$ is the angular diameter distance. For a certain redshift range, $`z1.32.3`$ (for $`\mathrm{\Omega }_m=0.3`$), the competition between the two ingredients causes a positive correlation between $`w_0`$ and $`w_a`$. A critical pitfall however was that systematic uncertainties could blow up the constraint contour in a preferred direction such that the resulting contour actually exhibited an anticorrelation (compare Figs. 5 and 8 of linap ). Another indication of positive correlation was found between discrete values of $`w(z)`$ in neighboring redshift bins by songknox  for ratios of distances entering weak gravitational lensing calculations.

Here we investigate distance ratios that appear in strong gravitational lensing measurements. Specifically, we consider $`r_{ls}/r_s`$, related to the angular separation $`\mathrm{\Delta }\theta `$ between multiple images of a source, $`r_{ls}/(r_lr_s)`$, related to the time delay $`\mathrm{\Delta }t`$ between the images, and $`r_lr_{ls}/r_s`$, connected with the length scale of caustic properties, or the cross section, of the lensing. Here $`r_l`$ is the comoving distance to the lens, $`r_s`$ to the source, and $`r_{ls}`$ between the source and lens; in a flat universe $`r_{ls}=r_sr_l`$.

Figure 2 illustrates the sensitivity of each of these ratios to the cosmological parameters as a function of redshift. We have taken the redshift to represent the lens redshift $`z_l`$ and for convenience fixed the source redshift to be $`z_s=2z_l`$ (since that gives roughly the greatest probability of effective lensing). Generally the greatest sensitivity is to the magnitude of the matter (or dark energy) density, then to the present value of the dark energy equation of state ratio, and least to the time variation of the equation of state.

Of most interest to our investigation are the crossings from negative to positive sensitivity, i.e. at some redshift the distance ratio switches from diminishing to growing as we change the value of a parameter $`\theta =\{\mathrm{\Omega }_m,w_0,w_a\}`$ (increasing it, say). Since these crossings occur at different redshifts if the parameter is $`w_0`$ than if it is $`w_a`$, then the correlation between these quantities can shift from negative to positive. We see that positive correlation occurs for $`r_{ls}/r_s`$ for $`z_l0.91.7`$ and for $`r_{ls}/(r_lr_s)`$ for $`z_l00.6`$; zero crossing in equation of state variables and hence positive correlation does not occur for $`r_lr_{ls}/r_s`$. (The exact values will depend on the fiducial cosmology, here taken to be $`\mathrm{\Omega }_m=0.3`$, $`w_0=1`$, $`w_a=0`$.)

For the ratio $`r_{ls}/r_s`$ we expect rapid evolution in the $`w_0w_a`$ degeneracy direction as the contour will rotate from horizontal at $`z=0.9`$ (no sensitivity to $`w_0`$, so lying parallel to the $`w_0`$-axis) to vertical at $`z=1.7`$ (no sensitivity to $`w_a`$). Fig. 3 shows this behavior in a “flower” plot, with idealized precision to better illustrate the degeneracy direction. A strong lensing survey measuring image separations from lenses within this redshift range should possess good internal complementarity and the ability to constrain the dark energy equation of state. Furthermore, in this range strong lensing will have near orthogonality with supernova distance and weak gravitational lensing surveys and so be a valuable complementary probe.

The ratio $`r_{ls}/(r_lr_s)`$ appears somewhat less promising. Its sensitivity to $`w_a`$ is much weaker in the key redshift range and the evolution of the correlation between $`w_0`$ and $`w_a`$ much less dramatic. In particular the contour never tilts far from the vertical. This means it is not as complementary either internally or with other cosmological probes. But it does possess three interesting characteristics: its overall sensitivity to the equation of state is reasonably large at $`z>1`$, the positive correlation region is at low redshifts where observations are easier (though there is less volume, hence fewer lens systems), and it has an odd null sensitivity with respect to $`\mathrm{\Omega }_m`$ at $`z=0.35`$. This last is where its sensitivity to $`w_a`$ is maximal in the positive correlation region and eases degeneracy with the matter density and hence reduces the need for a tight external prior on it. So time delay measurements should not be wholly neglected as a possibly useful probe.

## III Constraints on Cosmological Parameters

So far we have only concentrated on the degeneracy direction. What about the actual magnitude of constraints capable of being imposed on the cosmological parameters? If we consider one parameter $`\theta `$ at a time, fixing the others, we can obtain a lower limit on the parameter estimation uncertainty. This is given by

$$\delta \theta \left|\frac{[\mathrm{Distance}\mathrm{Ratio}]}{\theta }\right|^1\delta [\mathrm{Distance}\mathrm{Ratio}].$$
(5)
In Fig. 2 we took a fractional measurement error of 1%. Thus, for example, the sensitivity for $`r_{ls}/r_s`$ with respect to the parameter $`w_0`$ at $`z=0.3`$ is -6.5 so the best possible estimate of $`w_0`$ from such a single, 1% measurement is 1/6.5=0.15. That is, the unmarginalized uncertainty is $`\sigma (w_0)=0.15`$. This will be degraded by degeneracies with other parameters, or by worse precision in the measurement, and improved by additional measurements, subject to some systematics floor.

In fact, degeneracies with additional parameters wash out the simple orthogonality of Fig. 1. The parameter phase space is actually 3-dimensional (from $`\mathrm{\Omega }_m`$) or higher and the joint probabilities do not merely trace the intersection of the contours in the equation of state plane. This weakening of apparent complementarity also means that the magnitude of the sensitivity with redshift matters as well as the degree of correlation, and this can alter the favored redshift range. All these effects need to be taken into account.

Here we seek to establish what improvements the complementarity of strong lensing data with other probes such as distance-redshift or weak lensing data can realistically provide. Next generation precision in all three of these come in one package: the proposed Supernova/Acceleration Probe (SNAP: snap ). With some $`10^7`$ galaxies observed in its deep survey (down to AB mag $`R=30`$) and $`10^8`$ in its wide survey (to AB mag $`R=27.5`$), the canonical estimate of one strong lens system per thousand galaxies provides an anticipated wealth of data. Other valuable data sets will come from LSST lsst  and LOFAR lofar .

### III.1 Image Separation

If we consider an optimistic strong lensing measurement precision of 1% in a redshift bin of 0.1, whether through statistical means, e.g. 100 lens systems of 10% precision, or as a systematic floor, we find that the parameter uncertainties are large if we do not a priori fix some of the variables. Even employing the orthogonal combination of measurements of the strong lensing separation variable $`S=r_{ls}/r_s`$ at $`z=0.9`$ and $`z=1.7`$ (cf. Fig. 3), the estimations are $`\sigma (w_0)=0.15`$, $`\sigma (w_a)=0.53`$ even with $`\mathrm{\Omega }_m`$ fixed, and 0.28, 0.93 with a prior on $`\mathrm{\Omega }_m`$ of 0.03. The lack of sensitivity of the strong lensing separation to the dark energy equation of state and, as mentioned above, the degeneracy with the other parameter, $`\mathrm{\Omega }_m`$, prevent as well its complementarity with other probes from being as effective as might have been hoped from §II.

Figure 4 shows contours in the $`w_0w_a`$ plane, now marginalized over $`\mathrm{\Omega }_m`$, for various combinations of probes. This considers 1% measurements at $`z=0.9`$, 1.3, 1.7. The improvements in parameter estimation remain modest: about 15% in $`w_0`$ and $`w_a`$ when added to SNAP supernova distance data, 19% and 3% respectively when added to the supernova plus weak lensing data set, 25% and 13% when added to supernova plus Planck CMB data, and 17% and 6% when added to all three data sets. If strong lensing $`S`$ measurements were available from $`z=0.11.7`$, say (recall this is the lens redshift, with the source assumed to lie at twice this redshift), the improvement to the triple data set increases to 30% and 15%. If the measurement floor for the strong lensing was instead at 2%, then the constraints in conjunction with the other three probes only gain due to strong lensing by 6% and 2% for $`z=0.9`$, 1.3, 1.7 or 17% and 7% for $`z=0.11.7`$.

### III.2 Time Delay

The other strong lensing quantity we consider is $`T=r_{ls}/(r_lr_s)`$, related to the time delay between images. This introduces an extra parameter – the Hubble constant $`H_0`$, which we marginalize over. Again, measurements of strong lensing by itself provide only poor estimates of the cosmological parameters, even for combinations of redshifts where the degeneracy directions are complementary. The insensitivity to $`\mathrm{\Omega }_m`$ at low redshift noted in §II does not help since the overall weak dependence and degeneracy of the equation of state variables prevents precision constraints.

In complementarity with supernova distances, strong lensing $`T`$ measurements of 1% precision at $`z=0.10.6`$ provide improvement by 21% and 9% in $`w_0`$ and $`w_a`$. Adding them to all three other data sets helps by 13% and 9%, or 22% and 15% if the strong lensing measurements extend from $`z=0.11.7`$. Weakening the precision to 2% adjusts this last case to 7% and 5%. Note that no external prior on $`H_0`$ is required; it is determined from the data to better than 1%.

### III.3 Image Separation and Time Delay

Simultaneous addition of both separation and time delay information allows further improvements. Relative to the supernova only data, these are 36% and 41% for 1% strong lensing data at $`z=0.1,0.3,\mathrm{}1.7`$; added to all three other probes they yield tighter constraints by 31% and 17%. With only 2% precision the improvements in the latter case are 15% and 7%. We show the parameter constraints in Fig. 5 for the three probes without strong lensing, with strong lensing for the fiducial redshift bins, for a reduced strong lensing redshift range, and for the weaker 2% strong lensing precision.

The highest redshift bins are not that powerful but it is important to include the $`z=0.1`$ bin to constrain $`H_0`$, which is strongly degenerate with $`w_0`$ and $`w_a`$. Even with $`H_0`$ known to better than 1% from the data, the degeneracy still plays an important role since the strong lensing time delay probe is so much more sensitive to $`H_0`$ than to the equation of state.

As usual, taking the fiducial cosmology to be that of a cosmological constant tends to underestimate both the sensitivity and complementarity of the probes linbo . For comparison Fig. 5 includes contours for a supergravity inspired dark energy model with time varying equation of state: S(UGRA) with $`w_0=0.82`$, $`w_a=0.58`$. Here strong lensing offers more dramatic improvements, by 49% in estimating $`w_0`$ and 54% in $`w_a`$ for 1% precision and 29% and 22% for 2% precision. Even at 2% precision the final uncertainty from the four probes in complementarity achieves an impressive constraint on the time variation of $`\sigma (w^{})=0.05`$.

## IV Conclusion

While strong lensing does not completely fulfill the promise of orthogonal constraints on the dark energy equation of state that would provide vastly improved parameter estimates, it does offer some complementarity, especially for models with time varying equation of state. Furthermore, the data resources required will be innate within the deep and wide optical and near infrared surveys of SNAP; abundant strong lensing data should also come from LOFAR at radio wavelengths and both strong and weak lensing from LSST in the optical. We find that strong lensing adds reasonable value to dark energy constraints when the measurements reach the 1% precision level. This is in accord with other calculations (e.g. yamamoto ; lewisibata ). Precision arises from statistics subject to a systematics floor.

The challenge in the next generation will not be the search for sufficient strong lensing data, it will be achieving systematics limits sufficiently low to allow the statistical wealth and the value of strong lensing complementarity to come into play. Currently 1% sounds rather optimistic for a bound to systematic uncertainties. The main contribution is likely to come from our ignorance of the lensing mass model. For example, within a singular isothermal sphere model 1% measurement of the distance ratio $`S`$ requires knowledge of the velocity dispersion to 0.5% dhk . Cross correlation of different sources with the same lens to remove systematics, à la bernjain  for weak lensing, would not work, both because of the rarity of strong lensing and because sources at different distances and positions would probe different parts of the lensing mass distribution.

One might hope to constrain the lensing mass distribution by making use of the simultaneous weak lensing information that SNAP or LSST would provide, but this will be of little use according to dalal . Basically the weak lensing measurements of shears average over a larger range of scales (or multipoles) than the convergence contribution from the mass to the strong lensing. It is like trying to locate a pin with thick gloves. Nongaussian effects also need to be incorporated in a more rigorous treatment of strong lensing. Relaxation of the simplification that sources lie at redshifts $`z_s=2z_l`$ is unlikely to change significantly the results shown here.

The advantage of image separations and time delays is that they are well defined and used as markers of the cosmic geometry. One could also use the statistics of lensing systems, e.g. the number of lenses out to some redshift (e.g. sarbu ; oguri ; kuhlen ), but then one must deal with astrophysical scatter and observational bias and incompleteness. Furthermore, lensing counts depend on mass density growth and volume factors that take us further from our original goal of finding orthogonal probes in the equation of state plane.

Since systematics are the limiting factor, we can attempt to seek out special lensing systems that ameliorate these uncertainties, such as strongly lensed calibrated candles (e.g. Type Ia supernovae holz ; ogurik ), time evolving, differentially amplified sources (e.g. Type II supernovae goobar ; wagoner ), or relatively clean Einstein cross or ring images. With the wealth of future data such subsampling may be practical.

While interrelations between cosmological parameters entering into distances and the expansion rate prevent clean orthogonality between cosmological probes, the various methods retain sufficient complementarity to both crosscheck each other and further tighten constraints on the nature of dark energy. If the systematics challenge can be met, in special systems at least, strong lensing merits further consideration in joining the cosmological toolbox of next generation probes.

## Acknowledgments

I thank Dragan Huterer for helpful comments. This work has been supported in part by the Director, Office of Science, Department of Energy under grant DE-AC03-76SF00098.


text:
# Three-dimensional phase-field simulations of directional solidification

## I Introduction

Microstructural pattern formation during solidification has been intensely investigated for decades Kurz , and yet progress remains to be made on many fascinating questions. A central one is pattern selection: For a given material under given processing conditions, what is the microstructural pattern that will emerge, and how is it selected out of all other possible structures? This is an important issue both from an applied and a fundamental viewpoint. For the physicist, solidification microstructures are a classic example of pattern formation out of equilibrium, and one would like to have a fundamental principle that is capable to predict the final patterns, as the minimization of an appropriate thermodynamic potential does for equilibrium situations. However, as of yet, no such general selection principle is known for out of equilibrium situations. While empirical rules for solidification microstructures, such as the maximum growth speed or the minimum undercooling criteria, are often very useful in practice, it is often unclear by what detailed mechanisms the final pattern is attained dynamically. Therefore, it is important to understand pattern selection more thoroughly starting from the equations for the evolution of the solidification front.

Recently, modern methods of statistical physics have allowed theorists to make considerable progress along this line. In particular, the phase-field method has emerged as a powerful tool to simulate microstructure evolution. Its main advantage over more traditional methods for front motion is that it avoids an explicit tracking of the interfaces with the help of one or several phase fields, continuous auxiliary fields which distinguish between thermodynamic phases. This makes it possible to simulate the full dynamics of complex morphologies both in two and three dimensions. While up to the mid-1990s, this method had remained a qualitative tool, recent progress has made it possible to treat quantitatively some occurrences of the classic free-boundary problem of solidification Karma98 ; Karma01 ; Folch03  and to match the results with theory and experiments Karma00 .

Here, we present three-dimensional phase-field simulations of the directional solidification of dilute and eutectic binary alloys. Both cases are characterized by the existence of a continuous family of periodic steady-state solutions with different spatial periodicity (spacing): cells for dilute, lamellae or rods for eutectic alloys. It has been demonstrated both in thin-sample experiments Akamatsu98 ; Ginibre97  and two-dimensional simulations Kopcz96 ; Karma96 ; Dresden  that no sharp pattern selection occurs: the final spacing depends on the growth history, and a range of stable spacings is observed. Outside of this range of spacings, steady-state solutions generally exist, but they are unstable, which explains why they cannot be dynamically selected. Therefore, the limits of stability determine the observable spacings. The different instabilities and their onsets have been determined numerically for periodic two-dimensional states Kopcz96 ; Karma96 ; we extend this work to three dimensions, where new instabilities can arise.

## II Dilute Binary Alloy

Thin-sample directional solidification offers a unique opportunity to study the influence of crystalline anisotropy on the interface dynamics experimentally. The anisotropy in the sample plane is determined by the orientation of the three-dimensional crystal with respect to that plane. If a crystalline (100) axis is normal to the sample plane, the in-plane anisotropy is unchanged with respect to the three-dimensional situation; in contrast, if the crystal is oriented with a (111) axis normal to that plane, it is close to zero. It was observed that shallow cells are stable only for anisotropic crystals; in the absence of anisotropy, cell splitting and elimination events occur continuously, and the front as a whole never reaches a steady state Akamatsu98 . This is in agreement with a numerical study in two dimensions by the boundary integral technique Kopcz96 , which indicated that the range of stable cellular states strongly depends on anisotropy. In addition, this study identified the relevant instability modes that limit the stable range: elimination of every other cell on the short-spacing side, and a period-doubling oscillatory mode on the large-spacing side.

We have extended this investigation to three dimensions using a phase-field model that has equal diffusivities in both phases, and a constant concentration jump (partition coefficient $`k=1`$). While this is not very realistic for alloys, the qualitative behavior is similar for a one-sided model Karma01 , and the implementation PlappCells  is particularly simple. Simulations that started from a flat steady-state interface with a small random perturbation are shown in Fig. 1. Clearly, the crystalline anisotropy has a dramatic effect: whereas the interface spontaneously organizes into an array of hexagonal cells for sufficiently strong anisotropy, the evolution remains unsteady up to the end of the run for only slightly lower anisotropy.

To determine the stability boundaries for regular arrays, we generated states with well-defined spacings by placing six cells in simulation boxes of appropriate size and aspect ratio, and changing control parameters until an instability occurred. It turns out that the spatial structure of the relevant instability modes is determined by the symmetry of the underlying steady-state pattern. A hexagonal lattice can be decomposed into three equivalent sublattices of $`\sqrt{3}`$ times larger spacing. Instabilities occur when the symmetry between the three sublattices is broken. For small spacings, one out of the three sublattices is eliminated. For large spacings, the sublattices start to oscillate, with a phase difference of about $`2\pi /3`$ between each other. In both cases, once the first cell eliminations or splittings occur, the system follows an unsteady evolution as in Fig. 1b. The stability boundaries are plotted in Fig. 2. The stability limit of the oscillatory mode shifts to larger spacings with increasing anisotropy. This explains the results of the larger simulations in Fig. 1: Only when a sufficiently large range of stable spacings is available, the system is able to “find” a steady state within a reasonable simulation time.

## III Eutectic Alloy

Coupled eutectic growth leads to lamellar or fibrous microstructures. Detailed studies for lamellae in two dimensions by thin-sample experiments Ginibre97  and by dynamic boundary-integral simulations Karma96  have shown that the range of stable spacings is limited on the short-spacing side by a long-wavelength lamella elimination instability, and on the large-spacing side by oscillatory instabilities of once or twice the initial spacing. In contrast, very little is known about pattern stability in three dimensions, since no reliable numerical method was available that could handle complex composite microstructures.

The phase-field method would seem ideally suited to overcome this difficulty. However, its use has been hampered until recently by the wide range of length scales involved, from an interface which is diffuse on a nanometric scale up to a diffusion layer of several millimeters for slow solidification. The problem is that the variation of the phase fields through the diffuse interfaces has to be resolved by the discretization, which would result in huge simulation boxes. Fortunately, quantitative phase-field models offer the possibility to use an interface thickness $`W`$ much larger than that of real solidification fronts, by ensuring that results are independent of $`W`$. The only limitation is then that $`W`$ has to remain about an order of magnitude smaller than the radius of curvature of the interface Karma98 ; Karma01 .

We have recently developed a quantitative model for two-phase solidification with vanishing solute diffusion in the two solid phases Folch03 . As a test, we have performed simulations in two dimensions at the eutectic composition of a model alloy that has a completely symmetric phase diagram (generalization to other compositions and arbitrary phase diagrams is straightforward Folch03 ), and compared the results to a boundary integral calculation with the code of Ref. Karma96 . As can be seen in Fig. 3, our model produces precise results that are independent of the interface thickness for $`\lambda /W64`$, where $`\lambda `$ is the lamellar spacing, whereas those of a generic model make appreciable errors and depend on $`W`$.

In Fig. 4, we show a three-dimensional simulation of an unstable lamella with reflecting boundary conditions on all sides. The lateral dimensions are $`56\times 168`$ grid points, and the diffusion field along the growth direction is calculated using a multi–grid method with coarser and coarser grids away from the interface. Undulations develop with time, and finally the lamella breaks up, leaving behind a maze of the two phases (not shown). This evolution can be further analyzed to characterize the instability. First, the positions of the triple lines are extracted. Then, we calculate the position of the mid-line and the local width of the lamella. These two variables correspond to two different instability modes: the former to an undulation of the whole lamella, the latter to a pinching mechanism. After a Fourier transform, it can be checked that modes of different wave number $`k`$ along the lamella grow exponentially with a $`k`$-dependent growth rate $`\omega (k)`$. The extracted stability spectrum is plotted in Fig. 5. The “undulating” instability grows faster, and indeed it is this mode that shows up in Fig. 4. It should be noted that this three-dimensional instability occurs for parameters where the two-dimensional system is stable both in simulations and experiments. Therefore, it can already be concluded from this preliminary simulation that the range of stable spacings is considerably smaller in three dimensions than in two.

## IV Conclusion

Phase-field modeling offers a unique opportunity to investigate pattern stability and pattern selection in three dimensions. Here, we have focused on results for generic phase diagrams, but the method can be carried over to the study of particular materials; simulations for the alloy of Refs. Akamatsu98 ; Ginibre97  are currently in progress. This will enable us to carry out critical comparisons of simulations, experiments, and theory, which constitutes a particularly promising way to elucidate many long-standing questions concerning microstructural pattern formation during solidification.

## V Acknowledgments

We thank Alain Karma for many useful discussions, and for the boundary integral code. This work was supported by the Centre National d’Etudes Spatiales, France.


text:
# Large scale numerical simulations of “ultrametric” long-range depinning

## I Introduction

The depinning of an elastic interface in a random environment gives a common theoretical framework to describe physical phenomena as various as the advance of a magnetic wall, the propagation of a fracture front or the wetting of a disordered surface (see e.g. Ref. Leschhorn et al. (1997); Kardar (1998); Fisher (1998) for a recent review). The richness of the physics encountered in these different situations results from the same feature. The disorder of the environment which tends to anchor the interface competes with an elastic-like term which tends to keep the interface smooth. The tuning of an external driving force allows to go through a critical transition. Below threshold the interface can only advance over a finite distance before stopping in a blocked conformation. Above threshold the interface can move freely and acquire a finite velocity. At threshold the system is characterized by a set of universal critical exponents.

In the case of over-damped dynamics the motion can be described by the following stochastic equation:

$$\mu _th(x,t)=F_{ext}(t)+F_{el}(x,t)+\gamma [x,h(x,t)]$$
(1)
where $`h(x,t)`$ denotes the position of the front, $`F_{ext}(t)`$ the external driving force, $`F_{el}(x,t)`$ the elastic term due to the distortion of the front and $`\gamma [x,h(x,t)]`$ the frozen disorder. Depending on the physical phenomenon considered the elastic term can take different forms. In the case of magnetic walls or of wetting in a Hele-Shaw cell, the interactions to be taken into account are short ranged and at first order the elastic force can be estimated by a simple Laplacian term: $`F_{el}(x,t)=^2h(x,t)`$. In the case of the advance or receding of a triple line on a disordered substrate Joanny and de Gennes (1984) or of the propagation of fracture front long range interactions have to be considered. In the latter case elastic interactions are mediated via the bulk along the whole front. A first order perturbation analysis for the front roughness givesGao and Rice (1989):

$$F_{el}(x,t)=\frac{1}{\pi }𝑑x^{}\frac{h(x,t)h(x^{},t)}{|xx^{}|^2}$$
(2)
This long range elastic string model has been widely studied over the last ten years. In particular the interface is shown to exhibit a self-affine roughness: the width $`w`$ of the interface scales with the system size $`L`$ as $`wL^\zeta `$. The different numerical works Dong et al. (1993); Måløy and Schmittbuhl (1997); Tanguy et al. (1998); Rosso and Krauth (2002) performed give estimations of the roughness exponent $`\zeta `$ spreading in the interval $`[0.340.40]`$. The latter exponent and more generally the universality class of the model depend strongly on the long range nature of the kernel in Eq. (2). For instance, if $`F_{el}^2h`$, then $`\zeta =1.25`$. The most recent results obtained by Rosso and Krauth Rosso and Krauth (2002) give a value $`\zeta 0.39`$ significantly larger than the theoretical prediction $`\zeta =1/3`$ obtained by one-loop calculations of renormalization group techniqueNarayan and Fisher (1993); Ertaş and Kardar (1994) and equally smaller than the recent two-loop estimationChauve et al. (2001) $`\zeta 0.47`$. Note however that these values are not consistent with the results $`\zeta 0.50.6`$ obtained experimentally for interfacial fractureMåløy and Schmittbuhl (1997) or wettingPrevost et al. (1999); Moulinet et al. (2002). Therefore it is of uttermost importance to have an accurate determination of these critical exponents, and thus to be able to study large system sizes.

In the following we present an ultrametric version of an extremal model of depinning (see Ref. Tanguy et al. (1998) for a detailed study of the original model with Euclidean metric). The complexity of the elementary step of calculation is shown to scale with the system size $`L`$ as $`\mathrm{log}_2L`$ instead of $`L`$ in the original model. This allows us to perform simulations on systems of size $`L=2^{20}10^6`$, which corresponds to a gain of two to three orders of magnitude compared with other published works based on the Euclidean metrics. The universality class of the model is shown to be unchanged. Beyond the numerical acceleration, this ultrametric model of depinning may thus also serve as a starting point for a real space renormalization analysis.

The paper is organized as follows: In a first part we recall the definition of the original model, we then define the ultrametric version and present the main features of the new model. In the second part we give results of simulations performed on large systems and focus on the critical properties of this ultrametric model of depinning. We conclude that it lies in the same universality class as its euclidean version. The details of the numerical implementation of the algorithm are finally given in an appendix.

## II Extremal driving of a depinning front

Various numerical techniques can be used to study the depinning phenomenon in the vicinity of the critical threshold. Early works used a direct integration of the equation above via EulerDong et al. (1993) or Runge-Kutta Zapperi et al. (2000) schemes. Recently Rosso and Krauth Rosso and Krauth (2002) developed an iterative algorithm to determine the blocked conformations corresponding to a given constant forcing. They used periodic boundary conditions both along the front and in the direction of propagation and the critical threshold is reached when the last blocked conformation starts moving.

In the case of strong pinning the advance of the front proceeds by successive local instabilities. This avalanche behavior is characteristic of the motion of a depinning front. This property can be exploited to develop an efficient algorithm describing the motion of the front close to the critical threshold. Instead of driving the front at a constant external force, it consists of tuning the latter at the exact value such that one and only one site can depin at a time. So doing, the sequence of depinning events is preserved. The implementation of this extremal dynamics which has been used since 1992 in various interface growth models Zaitsev (1992); Sneppen (1992); Bak and Sneppen (1993); Leschhorn and Tang (1994); Roux and Hansen (1994); Paczuski et al. (1995) is straightforward. At each iteration step one needs to identify the weakest site, to advance it up to the next trap and to update the long range elastic forces due to the change of front conformation. The latter operation scales with the size of the front. The great advantage of this method is that the system remain constantly at the edge of the critical behavior and it is not necessary to tune the external driving force (see e.g. Ref. Dickman et al. (2000) for a discussion on the use of extremal dynamics to reach the critical state and more generally on the link between Self Organized Criticality and classical critical transitions).

Based on the above described extremal modeling, we now detail the way to implement the model numerically. The front is discretized along a regular horizontal grid; $`i[0,L1]`$ and $`h_i`$ are the coordinates along the front and in the direction of propagation respectively. Traps of random depths $`\gamma _i`$ are randomly distributed along the direction of propagation. The distortion of the front induces elastic forces $`f_i^{el}`$ via a Green function $`G_{ij}=G(r_{ij})`$ where $`r_{ij}`$ is the euclidean distance separating two sites $`i`$ and $`j`$ along the front. In the case of a fracture front, a discretized version of the elastic redistribution function is such that

$$G_{ij}_{ij}|ij|^2,G_{ii}=\underset{ij}{}G_{ij}.$$
(3)
Let us consider a given conformation of the front. For each site $`(i,h_i)`$ located in a trap of depth $`\gamma _i`$, we can define a local depinning threshold $`s_i=\gamma _if_i^{el}`$: this site depins as soon as the external driving force $`F`$ overcomes the threshold $`F>s_i`$. The depinning threshold $`s(t)`$ of the front conformation obtained at iteration $`t`$ thus corresponds to the minimal external force to be exerted so that at least one site of the front can depin: $`s(t)=\mathrm{min}_is_i`$. Finally the critical threshold $`s^{}`$ above which the front can freely propagate is $`s^{}=\mathrm{max}_ts(t)`$. The basic rule of the extremal driving consists simply at each iteration step $`t`$ of simply tuning the external force at exactly the value of the depinning threshold of the current front conformation: $`F(t)=s(t)`$. Once identified, the extremal site is advanced up to the next trap, the elastic forces are updated to take into account this local displacement and the new value $`s(t+1)`$ of the front depinning threshold is evaluated. We summarize below the elementary steps of the algorithm used to run the model and we estimate their complexity respectively to the size $`L`$ of the system.

| A | initialization $`h_i0,f_i^{el}0,\gamma _irnd`$ | $`[L]`$ |
| --- | --- | --- |
|  | $`i[0L],L=2^n`$ |  |
| B | identification of the extremal site $`i_0`$ | $`[L]`$ |
|  | such that $`\gamma _{i_0}F_{i_0}=\mathrm{min}_i(\gamma _if_i^{el})`$ |  |
| C | advance of the extremal site | $`[1]`$ |
|  | $`\delta h_{i_0}rnd,h_{i_0}h_{i_0}+\delta h_{i_0}`$ |  |
| D | update of the trap depth $`\gamma _{i_0}rnd`$ | $`[1]`$ |
| E | update of the elastic forces $`f_i^{el}f_i^{el}+G_{ii_0}\delta h_{i_0}`$ | $`[L]`$ |
|  | where $`G_{ij}|ij|^2`$ |  |
| F | back to step B |  |

where $`rnd`$ stands for a random number and $`ab`$ for the assignment of value $`b`$ to variable $`a`$. Except the first initialization step, the two limiting steps are the identification of the extremal site and the update of elastic forces along the front which both scale linearly with the system size $`L`$. This sequence of elementary steps is then iterated $`T`$ times to obtain statistical averages of the quantities of interest.

## III Ultrametric depinning

We now turn to the presentation of the ultrametric model. The basic rules of the extremal model remain identical but the redistribution of the elastic forces. Instead of using the natural euclidean distance along the front we use an ultrametric distance. The structure of the algorithm stays roughly similar to the previous one but steps B and E are shown to be characterized by a complexity in $`\mathrm{log}_2L`$ instead of $`L`$.

The most natural structure to be used in the context of a model with ultrametric distance is a dyadic tree. Let us first build such a tree whose final leaves are the $`L`$ sites of the depinning front. As illustrated on Fig. 1, the simplest ultrametric distance between two sites $`i`$ and $`j`$ is the number of branches that composes the shortest path on the tree between the two sites. This is exactly twice the depth $`m`$ of the nearest common ancestor of these two sites. Different choices of distance are possible based on the tree structure. In the following we use a definition that preserves the scaling of the Euclidean distance:

$$d_u(i,j)=2^m1$$
(4)
An important characteristic of this ultrametric distance is its degeneracy. Namely there is one point at distance $`d=2^11=1`$, two points at distance $`d=2^21=3`$ and $`2^{p1}`$ points at distance $`d=2^p1`$. For a set of $`L=2^n`$ points one thus counts only $`n=\mathrm{log}_2L`$ different values for the distance between two points of the set. There lies the main advantage of the choice of an ultrametric distance from the computational complexity point of view. The expression of the elastic Green function then derives directly from the definition of the new distance:

$`G_{ij}=G\left[d_u(i,j)\right]{\displaystyle \frac{1}{d_u(i,j)^2}}`$ (5)
$`G_{ii}={\displaystyle \underset{m=1}{\overset{\mathrm{log}_2L}{}}}2^{m1}G\left[2^m1\right]`$
Using this definition, we can easily accelerate the update of elastic forces (step E). As illustrated on Fig. 2 instead of updating $`L`$ sites, we can update only $`\mathrm{log}_2L`$ subtrees corresponding to sites located at the same ultrametric distance of the extremal site. A similar gain can be obtained on step B for the determination of the extremal site. Technical details regarding the numerical implementation of the algorithm are developed in the appendix. The basic steps of the algorithm then scale with $`\mathrm{log}_2L`$ instead of $`L`$. The price to pay for this numerical advantage is the loss of the translational invariance. We show below that it does not affect the universality class of the model.

The simulations have been performed on systems of sizes up to $`L=2^{20}`$. The numerical runs were performed over a large duration $`T`$ respective to the natural correlation time of the system $`\tau _LL^z`$. In the case of the largest system $`L=2^{20}`$ we used $`T=2.510^{10}100\tau _L`$.

## IV Characterization of the critical state: numerical results

The propagation of depinning fronts at the critical threshold exhibits a rich phenomenology. The front presents a self-affine roughness: its width $`w(\mathrm{\Delta }x)`$ measured over a distance $`\mathrm{\Delta }x`$ scales as $`w(\mathrm{\Delta }x)\mathrm{\Delta }x^\zeta `$, where $`\zeta `$ is the roughness exponent.

The dynamics of depinning is characterized by an avalanche behavior. In the framework of an extremal dynamics this can be described by the distribution of the distances $`r`$ between two successive depinning sites: $`P(r)r^a`$. Generalizing this distribution for sites corresponding to depinning events separated by a given time lag allows to obtain in addition the dynamic exponent $`z`$ which characterizes the spreading of the avalanches. Namely the lateral extension $`\xi `$ of an avalanche of duration $`\mathrm{\Delta }t`$ scales as $`\xi \mathrm{\Delta }t^{1/z}`$.

Another quantity of interest is the external force $`s`$ needed to depin a given conformation of the front. It can be shown Skoe et al. (2002); Vandembroucq and Roux (2003) that the distribution $`𝒬(s)`$ of these front depinning forces exhibits a singular behavior close to the critical threshold $`s^{}`$: $`𝒬(s)(s^{}s)^\mu `$.

In the following we present simulations of ultrametric depinning performed on large systems (up to $`L=2^{20}`$). We recover all critical features described above with exponents numerically indistinguishable from their counterparts in the Euclidean version of the model.

### IV.1 Self-affine roughness

Various statistical roughness estimators can be used to characterize the self-affine properties of a rough front. Consider for example the standard deviation $`\sigma (\mathrm{\Delta }x)`$ of the height differences between points separated by a distance $`\mathrm{\Delta }x`$. A self-affine front obeys $`\sigma (\mathrm{\Delta }x)\mathrm{\Delta }x^\zeta `$. Similarly the width $`w(L)`$ of a front of length $`L`$ (i.e. the standard deviation of the height distribution along the front) scales as $`w(L)L^\zeta `$. Fourier or wavelet transforms are also of standard use.

In the context of this study we can also design a “wavelet like” roughness estimator which exploits the natural tree structure associated to the ultrametric distance. At the level $`\mathrm{}=n=\mathrm{log}_2L`$ we define $`\omega ^2(n)`$ as the variance of the height difference between nearest neighbors:

$$\omega ^2(n)=\left(h_ih_j\right)^2_{d_u(i,j)=2^n1}$$
(6)
At the upper level the height of a node is thus simply defined as the arithmetic average of its two ancestors: $`h^{(\mathrm{}1)}(i)=[h^{(\mathrm{})}(2i)h^{(\mathrm{})}(2i+1)]/2`$ and the corresponding variance $`w(m)`$ is computed. This sequence is iterated up to the root of the tree. At each level $`m`$ corresponds an ultrametric distance $`m=n\mathrm{}+1`$ and we have

$$\omega ^2(m)2^{2\zeta (m1)},\mathrm{or}\omega (m)d_u(m)^\zeta $$
(7)
We present now numerical results obtained for these various roughness estimators. The simulations have been performed on systems of sizes up to $`L=2^{20}`$. In Fig. 3 and 4 we show the scaling behaviors obtained for the wavelet roughness estimator $`\omega (d_u)`$ and the width $`w(L)`$ of the interface. We obtain perfect power law behaviors over six decades and we only show here in logarithmic scale the residuals after normalization by a power law of exponent $`\zeta `$. This procedure is a very sensitive way of detecting deviations from a power law. Note that previous published works deal with $`\mathrm{log}_2L10`$.

We present simulations performed with two kinds of distributions for the trap depths, uniform and Gaussian (respectively denoted by the subscripts $`a`$ and $`b`$ in the figures).

An estimation of the roughness exponent can be extracted from each individual set of data. Note that the fluctuations of this estimate due to the choice of the nature of disorder or the roughness estimator are larger than the deviation form the power law behavior itself. All results are presented here with the central value $`\zeta =0.391`$. This value slightly underestimates the results obtained from the width estimator and slightly overestimates the results obtained from the wavelet estimator. A conservative estimate of the roughness exponent thus appears to be:

$$\zeta =0.391\pm 0.005$$
(8)
### IV.2 Avalanche behavior – Dynamic exponent

The avalanche behavior is characteristic of the dynamics of the front close to threshold forcing. Although an extremal driving does not allow to recover the real dynamics of the front, avalanches associated to a given level of the driving force $`F`$ can be reconstructed from the history of the extremal force signal $`s(t)`$. An avalanche thus consists of a continuous series of depinning events such that $`s(t)<F`$. Instead of reconstructing these avalanches it is classical in the framework of extremal models Furuberg et al. (1998) to work directly from the extremal force signal. Let us consider a time lag $`\mathrm{\Delta }t`$ (a number of iterations), we introduce the distance $`\mathrm{\Delta }x`$ along the front between the sites depinning at $`t`$ and $`t+\mathrm{\Delta }t`$ respectively. It appears Tanguy et al. (1998) that the distributions of these distances $`\mathrm{\Delta }x`$ at fixed time lad $`\mathrm{\Delta }t`$ can be rescaled on a universal form:

$$P(\mathrm{\Delta }x;\mathrm{\Delta }t)=\frac{1}{\mathrm{\Delta }x^a}\psi \left(\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t^{1/z}}\right)$$
(9)
where $`\psi (u)u^a`$ for $`u1`$ and $`\psi (u)\mathrm{cste}`$ for $`u1`$. The exponent $`a`$ is well approximated by the exponent of elastic kernel of Eq. (2) $`a2`$ and the dynamic exponent $`z`$ can be related to the roughness exponent: if a sequence of $`\mathrm{\Delta }t`$ depinning events spreads over a distance $`\mathrm{\Delta }x`$ along the front, the knowledge of the roughness of the front over $`\mathrm{\Delta }x`$ leads to $`\mathrm{\Delta }t\mathrm{\Delta }x\mathrm{\Delta }x^\zeta `$ thus $`z=1+\zeta `$Tanguy et al. (1998); Krishnamurthy et al. (2000).

This scaling is recovered in the framework of ultrametric depinning, where we measured $`P(d_u,\mathrm{\Delta }t)`$. We check on Fig. 5 that after rescaling all distributions collapse on a unique master curve. The rescaling was obtained with the value $`z=1.39`$, for large arguments the behavior of $`\psi `$ is well approximated by a power law of exponent $`a=2`$.

## V Scaling functions

In addition to critical exponents, the pinned state is also characterized by scaling functionsRacz (2002). In the specific case of interface dynamics, this means that the critical properties of the front are described by a universal exponent and a universal function describing the fluctuations of the (rescaled) width of the interface. This property has been evidenced for various growth models (Edwards-Wilkinson, Kardar-Parisi-Zhang…)Racz (2002) and has been recently applied to the case of depinning interfaces Rosso et al. (2003). Note that in the latter case, beyond the interface width, the technique can be used to characterize other fluctuating quantities. In particular, the distribution of the depinning threshold of a finite elastic line under extremal driving can be shown to be universalSkoe et al. (2002); Vandembroucq et al. (2003). In the present study we show that the choice of an ultrametric distance does not affect these universal distributions. More precisely the fluctuations of interface width and depinning threshold are shown to be described by universal functions and these functions appear to be very close to or identical to their counterparts obtained in the framework of a depinning model with Euclidean distance.

These statistical distributions are however sensitive to the boundary conditions (periodic boundary conditions vs isolated system)Antal et al. (2002). In the following we use periodic boundary conditions. In the ultrametric case the re-summation over all replicas induces an additional mean field contribution $`1/2L^2`$ equally shared by all $`G_{ij},ij`$.

### V.1 Universal width fluctuations

Following Ref. Rosso et al. (2003) we study the distribution of the rescaled width $`w^2/w^2`$ where $`w^2`$ is the temporal average of the width of a depinning front of finite extent $`L`$. In Fig. 6 we present results obtained for various system sizes in both cases of Euclidean and ultrametric distance. We observe that all rescaled distributions collapse onto a master curve whatever the size of the front and the Euclidean or hierarchical metric.

### V.2 Universal depinning force fluctuations

As developed above, the extremal dynamics gives a direct access to the fluctuations of the driving force needed to depin the front site by site. Most of these fluctuations simply correspond to the depinning of near neighbors (which receive the largest contributions of the elastic redistribution) and are highly sensitive to the details of the pinning force disorder. The other part of these fluctuations corresponds to depinning events taking place at a larger distance $`d`$ from the previous depinning site. In other terms, the front can be regarded as pinned over the distance $`d`$ between the two successive extremal sites. Conditioning the depinning force distribution to this distance $`d`$ between successive extremal sites thus allows to define distributions of size dependent effective thresholds. (Note that in the context of fracture, it can be seen as a distribution of effective toughnessCharles et al. (2003)). An effective threshold over a distance $`d`$ can also be seen as the force needed to entirely depin a front of size $`d`$. It can actually be shown Skoe et al. (2002); Vandembroucq et al. (2003) that as $`d`$ increases, the distributions $`P(s;d)`$ tend to peak and approach the critical threshold $`s^{}`$. This behavior is recovered in our ultrametric variant as illustrated on Fig. 7 where we plotted the global distribution $`P(s)`$ and the contributions $`P(s;d_u)`$ corresponding to the different ultrametric distances.

The typical elastic force fluctuations over a distance $`d`$ can moreover be estimated: $`\delta f^{el}d^{(1\zeta )}`$. Both the width $`\sigma _s(d)`$ and the gap to the critical threshold $`\delta _s(d)=s^{}s(d)`$ of the center of these distribution actually follow this scaling. As developed in Ref. Skoe et al. (2002); Vandembroucq et al. (2003) this property allows for a precise extrapolation of the critical threshold by extrapolation of the linear relationship between $`s(d)`$ and $`\sigma _s(d)`$ up to the force value canceling $`\sigma _s`$. Moreover this rescaling results in a collapse of all distribution over a single master curve:

$$P(s;d)=d^{1\zeta }\chi \left[d^{1\zeta }(s^{}s)\right].$$
(10)
We see on Fig. 8 that the choice of an ultrametric distance slightly changes the shape of the distribution obtained after rescaling. This may mean that these conditional distributions are more sensitive to the boundary conditions and the loss of translation invariance induced by the ultrametric model than a macroscopic quantity such as the width of the interface.

## VI Conclusion

We developed a depinning model based on the use of an ultrametric distance. This choice allows to reduce the complexity of an elementary step of computation to $`\mathrm{log}_2L`$ with respect to the system size $`L`$. We performed numerical simulations on systems of size $`10^6`$ more than two orders of magnitude larger than in previously published results. We propose an estimate of the roughness exponent $`\zeta =0.390\pm 0.01`$ consistent with the values measured for the Euclidean counterpart of the model. Moreover we obtained scaling functions either identical or very close in both cases. Therefore the choice of an ultrametric distance appears not to affect the universality class of the depinning transition.

Beyond the numerical efficiency, this model may thus serve as a starting point for a real space renormalization analysis. Indeed a coarse grained picture of the model preserves exactly the same structure with “dressed” thresholds and front advances, and thus the solution of a simple $`L=2`$ front model may open the way to an analytic determination of critical exponents and universal scaling functions.

## Appendix A Numerical implementation of the model

In this appendix we describe the algorithm used to run this ultrametric model of depinning. As above mentioned the structure of the algorithm is roughly similar to the original one but steps B and E are characterized by a complexity in $`\mathrm{log}_2L`$ instead of $`L`$.

Let us first consider the redistribution of elastic forces that take place after a depinning event. The increment of force at a given site $`j`$ only depends on its distance $`d_u(i_0,j)=2^m1`$ to the extremal site $`i_0`$. Actually $`2^{m1}`$ different sites are at this same distance from the extremal site $`i_0`$ and form a subtree (see Fig. 2). Instead of updating sequentially the elastic force increment on every sites it is possible to update the elastic force acting on the whole subtree. To this aim we define the force $`\phi _{i,\mathrm{}}^{el}`$ as the force transported by the $`i^{th}`$ branch of the $`\mathrm{}^{th}`$ level of the tree. The first level ($`\mathrm{}=1`$) here corresponds to the two branches at the root and the level ($`\mathrm{}=\mathrm{log}_2L`$) to the branches pointing toward the $`L`$ sites of the front (the leaves of the tree). Using this definition the elastic force acting on site $`i`$ is nothing but the sum of the forces transported by its ancestor branches:

$$f_i^{el}=\underset{\mathrm{}=1}{\overset{n=\mathrm{log}_2L}{}}\phi _{\frac{i}{2^n\mathrm{}},\mathrm{}}^{el}$$
(11)
In the hierarchical algorithm the step E thus consists of updating the branches of the force tree affected by the redistribution of elastic forces : $`\phi _{j_0,\mathrm{}}^{el}=\phi _{j_0,\mathrm{}}^{el}+G_{i_0,m}\delta h_{i_0},j_0=\{i_0/2^n\mathrm{}\}`$ where $`\{i\}`$ is the operation of exchanging the parity bit of $`i`$.

In the same spirit we can build a tree to determine the value and location of the extremal site. The determination of the extremal site requires a priori $`L1`$ operations of elementary comparisons to find $`s^{}=\mathrm{min}_i(\gamma _if_i^{el})`$. After each depinning event, the elastic forces are updated on every sites so that the same $`L1`$ operations have to be performed at each iteration. In the ultrametric model the situation is slightly different because the elastic force on a site can be written as a sum of the forces transmitted by the branches of the tree (see Eq. 11). As a consequence two sites $`i`$ and $`j`$ share all force components acting on branches above their nearest common ancestor.

$`f_i^{el}`$ $`=`$ $`\mathrm{\Phi }_m+{\displaystyle \underset{\mathrm{}=m+1}{\overset{n=\mathrm{log}_2L}{}}}\phi _{\frac{i}{2^n\mathrm{}},\mathrm{}}^{el}`$ (12)
$`f_j^{el}`$ $`=`$ $`\mathrm{\Phi }_m+{\displaystyle \underset{\mathrm{}=m+1}{\overset{n=\mathrm{log}_2L}{}}}\phi _{\frac{j}{2^n\mathrm{}},\mathrm{}}^{el}`$ (13)
$`\mathrm{\Phi }_m`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=1}{\overset{m}{}}}\phi _{\frac{i}{2^n\mathrm{}},\mathrm{}}^{el}={\displaystyle \underset{\mathrm{}=1}{\overset{m}{}}}\phi _{\frac{j}{2^n\mathrm{}},\mathrm{}}^{el}`$ (14)
We can then perform the comparison operations level by level going up the tree of elastic forces. To be more specific let us define the following hierarchical structure. At level $`n=\mathrm{log}_2L`$ we define

$$\sigma _{i,n}=\gamma _i,\alpha _{i,n}=i$$
(15)
then we proceed iteratively up to the root of the tree:

$`\sigma _{i,p1}`$ $`=`$ $`\mathrm{min}(\sigma _{2i,p}+\phi _{2i,p}^{el},\sigma _{2i+1,p}+\phi _{2i+1,p}^{el})`$ (16)
$$\begin{array}{ccccc}\alpha _{i,p1}\hfill & =\hfill & \alpha _{2i,p}\hfill & \mathrm{if}\hfill & \sigma _{2i,p}+\phi _{2i,p}^{el}<\sigma _{2i+1,p}+\phi _{2i+1,p}^{el}\hfill \\ & =\hfill & \alpha _{2i+1,p}\hfill & \mathrm{if}\hfill & \sigma _{2i,p}+\phi _{2i,p}^{el}>\sigma _{2i+1,p}+\phi _{2i+1,p}^{el}\hfill \end{array}$$
(17)
The location and the value of the extremal site are thus given by

$$s^{}=\sigma _{1,0},i_0=\alpha _{1,0}$$
(18)
Without prior knowledge, the computation for the determination of the extremal site of a given conformation is simply the sum of all comparisons at each level of the tree: $`_1^n2^{p1}=2^n1=L1`$. This is exactly the same result as in the standard case. Consider now the situation after the depinning event: the trap depth $`t_{i0}`$ is updated at the extremal site and only $`n=\mathrm{log}_2L`$ branches of the tree are altered by the elastic force distribution. A sequence of $`n`$ comparisons thus allows to find the new extremal site: we start at the former extremal site $`i_0`$,

$$p=n:j_0=i_0,\sigma _{j_0}=t_{j_0},\alpha _{j_0}=i_0$$
(19)
then for $`p=n`$ to $`p=1`$ we proceed iteratively to update the tree and determine the value $`\sigma _{1,0}`$ and location $`\alpha _{1,0}`$ of the new extremal site.

$$j_p=j_0/2^{np},k_p=\{j_p\}$$
(20)
$$\begin{array}{ccc}\sigma _{j_{p1},p1}\hfill & =\hfill & \mathrm{min}[\sigma _{j_p,p}+\phi _{j_p,p}^{el},\sigma k_p,p+\phi _{k_p,p}^{el}]\hfill \end{array}$$
(21)
$$\begin{array}{ccccc}\alpha _{j_{p1},p1}\hfill & =\hfill & \alpha _{j_p,p}\hfill & \mathrm{if}\hfill & \sigma _{j_{p1},p1}=\sigma _{j_p,p}\hfill \\ & =\hfill & \alpha _{j_p,p}\hfill & \mathrm{if}\hfill & \sigma _{j_{p1},p1}=\sigma _{k_p,p}\hfill \end{array}$$
(22)


text:
# Langevin equation with scale-dependent noise

(March 21, 2003)

## Abstract

A new wavelet based technique for the perturbative solution of the Langevin equation is proposed. It is shown that for the random force acting in a limited band of scales the proposed method directly leads to a finite result with no renormalization required. The one-loop contribution to the Kardar-Parisi-Zhang equation Green function for the interface growth is calculated as an example.

The Langevin equation is one of the most general approximations for the evolution of a dynamical system affected by fluctuating environment. It arises in the description of magnetic at the presence of magnetic field fluctuations, in the description of hydrodynamic turbulence, in stochastic quantization problems, in the description of interface growth and in a large variety of other problems .

In the most general form the Langevin equation can be casted as

$$\frac{\varphi (𝒙,t)}{t}=U[\varphi (𝒙,t)]+\eta (𝒙,t),\eta (x)\eta (x^{})=D(x,x^{}),$$
(1)
where $`U[\varphi ]`$ is the nonlinear interaction potential, $`\eta (t,𝒙)`$ is the Gaussian random noise, which standss for the fluctuations of the environment.The Minkovski-like $`(d+1)`$-dimensional notation $`x(𝒙,t),k(𝒌,\omega )`$ is used hereafter.

The standard way to solve the Langevin equation (1) is to introduce the small parameter $`\lambda `$ in the interaction potential $`U`$, and then solve the system iteratively in each order of the perturbative expansion. The averaging over the Gaussian random force $`\eta `$ consists in evaluation of the pair correlators $`\eta \eta `$. The procedure is simplified by the assumption of the Gaussian statistics of the random noise which allows one to take into account only even order correlators of the random noise: all terms containing the odd number of $`\eta `$ are equal to zero. The diagram technique for the iterative solution of the Langevin equation is often called the Wyld diagram technique . Similarly to that in quantum field theory, it requires elimination of loop divergences by renormalization group methods .

The structure of divergences arising in the perturbative solution of the Langevin equation depends on a particular type of the random force correlation function. Most approximations use $`\delta `$-correlated random force of the form

$$\eta (k_1)\eta (k_2)=(2\pi )^{d+1}\delta (k_1+k_2)D(𝒌_2).$$
(2)
The spatial part of the correlation function $`D(𝒌)`$ is assumed either to be a constant or to have a power behavior ($`D(𝒌)|𝒌|^\beta `$).

From the physical standpoint, a random force acting in a limited range of scales is often desirable. For instance, a limited band forcing is the case for the stirred hydrodynamic turbulence and stochastic interface growth. In this paper we propose a novel method to describe the limited band stochastic forcing. It is shown, that for a narrow band forcing, an appropriate choice of the random force correlation functions yields a theory which is free of loop divergences. The proposed method preserves the whole structure of the perturbation expansion, with only the space of functions being changed. In the limiting case of the scale independent forcing (2) all common results are preserved.

To study the dynamical system separately at each scale, following , instead of the usual space of the random functions $`f(x,)(\mathrm{\Omega },𝒜,P)`$, where $`f(x)\mathrm{L}^2(\mathrm{R}^\mathrm{d})`$ for each given realization of the random process, we turn to the multi-scale representation provided by the continuous wavelet transform

$$W_\psi (a,𝒃,)=|a|^{\frac{d}{2}}\overline{\psi \left(\frac{𝒙𝒃}{a}\right)}f(𝒙,)d^dx.$$
(3)
The wavelet transform is performed here only in the spatial argument, but not in its temporal argument of $`x`$. This is because the structure of divergences and other important properties of the considered processes are determined by their spatial localization.

The existence an uniqueness of the inverse wavelet transform

$$f(𝒙,)=C_\psi ^1|a|^{\frac{d}{2}}\psi \left(\frac{𝒙𝒃}{a}\right)W_\psi (a,𝒃,)\frac{dad𝒃}{a^{d+1}}$$
(4)
is provided by the admissibility condition imposed on the basic wavelet

$$C_\psi =S_d^1\frac{|\widehat{\psi }(𝒌)|^2}{|k|^d}d^dk=\frac{|\widehat{\psi }(a𝒌)|^2}{a}𝑑a<\mathrm{}.$$
(5)
Factor $`S_d`$, the area of a unit sphere in $`d`$ dimensions, stands for the rotationally symmetric wavelet $`\psi (𝒙)=\psi (|𝒙|)`$.

The use of the wavelet images instead of the original stochastic processes provides an extra analytical opportunities. What is important for us, is the possibility of construction of more than one set of random functions $`W(a,𝒃,)`$, the images of which have coinciding correlation functions in the space of $`f(x,)`$.

It is easy to check, that the image of the correlation function of the processes

$$\widehat{W}(a_1,k_1)\widehat{W}(a_2,k_2)=C_\psi (2\pi )^d\delta ^d(k_1+k_2)a_1^{d+1}\delta (a_1a_2)D_0,$$

coincides with that of the white noise

$$\begin{array}{ccc}\widehat{f}(k_1)\widehat{f}(k_2)\hfill & =& (2\pi )^dD_0\delta ^d(k_1+k_2)\hfill \\ \widehat{W}(a_1,k_1)\widehat{W}(a_2,k_2)\hfill & =& (2\pi )^dD_0\delta ^d(k_1+k_2)(a_1a_2)^{d/2}\overline{\widehat{\psi }(a_1k_1)\widehat{\psi }(a_2k_2)}.\hfill \end{array}$$

Therefore, modeling the random force in the space of wavelet coefficients, we can provide a narrow band forcing keeping intact all required properties in the common space.

As an example, let us consider the well known Kardar-Parisi-Zhang model of interface growth :

$$\dot{Z}\nu \mathrm{\Delta }Z=\frac{\lambda }{2}(Z)^2+\eta .$$
(6)
Substituting the wavelet transform in spatial argument

$$Z(x)=C_\psi ^1\mathrm{exp}(ı(𝒌𝒙k_0t))a^{\frac{d}{2}}\widehat{\psi }(a𝒌)\widehat{Z}(a,k)\frac{d^{d+1}k}{(2\pi )^{d+1}}\frac{da}{a^{d+1}}$$
(7)
into the equation (6), with the random force of the form

$$\begin{array}{ccc}\widehat{\eta }(a_1,k_1)\widehat{\eta }(a_2,k_2)\hfill & =& C_\psi (2\pi )^{d+1}\delta ^{d+1}(k_1+k_2)a_1^{d+1}\delta (a_1a_2)D(a_2,𝒌_2),\hfill \\ \widehat{\eta }(a,k)\hfill & =& 0.\hfill \end{array}$$
(8)
after straightforward calculations, we get an integral equation

$$\begin{array}{ccc}(ı\omega +\nu 𝒌^2)\widehat{Z}(a,k)\hfill & =& \widehat{\eta }(a,k)\frac{\lambda }{2}a^{\frac{d}{2}}\overline{\widehat{\psi }(a𝒌)}C_\psi ^2(a_1a_2)^{\frac{d}{2}}\widehat{\psi }(a_1𝒌_1)\widehat{\psi }(a_2(𝒌𝒌_1))\hfill \\ & & 𝒌_1(𝒌𝒌_1)\widehat{Z}(a_1,k_1)\widehat{Z}(a_2,kk_1)\frac{d^{d+1}k_1}{(2\pi )^{d+1}}\frac{da_1}{a_{1}^{}{}_{}{}^{d+1}}\frac{da_2}{a_{2}^{}{}_{}{}^{d+1}},\hfill \end{array}$$

from where, the one loop approximation for the Green function is yielded

$$G(k)=G_0(k)\lambda ^2G_0^2(k)\frac{d^{d+1}k_1}{(2\pi )^{d+1}}\mathrm{\Delta }(k_1)𝒌_1(𝒌𝒌_1)|G_0(k_1)|^2𝒌𝒌_1G_0(kk_1)+O(\lambda ^4)$$
(9)
where $`G_0^1(k)=ı\omega +\nu 𝒌^2`$ is the zeroth order approximation for the Green function, and

$$\mathrm{\Delta }(k)C_\psi ^1\frac{da}{a}|\widehat{\psi }(a𝒌)|^2D(a,𝒌)$$
(10)
is the scale averaged effective force correlation function. The Green function obtained with the random force (8), does not depend on scale $`\widehat{Z}(a,k)=G(k)\widehat{\eta }(a,k).`$

Similarly, for the pair correlation function $`ZZ`$, we get

$`C(a_i,a_f,k)`$ $`=`$ $`{\displaystyle \frac{\lambda ^2}{2}}|G_0(k)|^2\overline{\widehat{\psi }(a_i𝒌)\widehat{\psi }(a_f𝒌)}`$
$`{\displaystyle \frac{d^{d+1}k_1}{(2\pi )^{d+1}}|G_0(k_1)|^2|G_0(kk_1)|^2[𝒌_1(𝒌𝒌_1)]^2\mathrm{\Delta }(k_1)\mathrm{\Delta }(kk_1)}.`$
For the random stirring which does not depend on scale, the integration in equation (10), after substitution $`k_1=k_1^{}+\frac{k}{2}`$ in both equations (9)and (Langevin equation with scale-dependent noise), leads to the known result .

Let us consider a single band forcing

$$D(a,𝒌)=\delta (aa_0)D(𝒌)$$
(12)
and the “Mexican hat” taken as the basic wavelet

$$\widehat{\psi }(𝒌)=(2\pi )^{d/2}(ı𝒌)^2\mathrm{exp}(𝒌^2/2),C_\psi =(2\pi )^d.$$
(13)
Substituting (12) and (13) into the equation (9), after integration over the frequency, in the leading order in the small parameter $`x=|𝒌|/|𝒌_1|`$, we get the contribution to the Green function ($`d>2`$):

$$G(k)=G_0(k)+\lambda ^2G_0^2(k)\frac{S_d}{(2\pi )^d}\frac{a_0^3k^2}{\nu ^2}\frac{d2}{8d}_0^{\mathrm{}}D(𝒒)e^{(a_0𝒒)^2}q^{d+1}𝑑q+O(\lambda ^4).$$
(14)
For constant $`D(𝒒)=D_0`$, the obtained contribution to the Green function is finite and does not require any further renormalization. In the limit $`\omega ,𝒌0`$, the one loop contribution to the surface tension $`\nu `$, which stems from (14), is equal to

$$\nu _{eff}=\nu \left[1\frac{\lambda ^2}{\nu ^3k^2}\frac{d2}{16d}\frac{S_d}{(2\pi )^d}a_0^{1d}D_0\mathrm{\Gamma }\left(1+\frac{d}{2}\right)+O(\lambda ^4)\right].$$
(15)
We have to note that the equations (14,15) were derived from one loop integral (9) after substitution $`k_1=q+k/2`$, followed by perturbation expansion in small parameter $`x=k/q`$, . For concrete correlation functions the direct numeric evaluation of the integral (9) may be more adequate. Similar calculations can be done for other interaction potentials. For instance, for the quadratic interaction $`\frac{\lambda ^2}{2}Z^2`$, the corresponding equations for (9) and (Langevin equation with scale-dependent noise) differ from that obtained above only by the absence of the scalar products of the wave vectors in each vertex.

Concerning the contribution of the higher orders of the perturbation expansion, we should say that for the basic wavelets $`\widehat{\psi }(k)`$ localized in $`k`$-space and limited band noise $`D(a,k)`$, the effective coupling constant $`\overline{\lambda }`$, which is the actual parameter of the perturbation expansion , can be made small by decreasing the noise amplitude. For instance, for the basic wavelets from the family $`\widehat{\psi }(k)=(2\pi )^{d/2}(ık)^n\mathrm{exp}(k^2/2)`$ and the noise correlator $`D(a,k)=D_0\delta (aa_0)`$, the effective coupling constant is

$$\overline{\lambda }^2=\frac{\lambda ^2}{\nu ^3}\frac{D_0}{a_0}.$$

The author is grateful to N.Antonov and V.Priezzhev for useful discussions.