L F Buchmann
email: buchmann@iesl.forth.gr

G M Nikolopoulos

O Zobay

P Lambropoulos

Passage-time statistics of superradiant light pulses from Bose-Einstein condensates

Keywords: numbers: 03.75.Nt, 67.85.-d, 37.10.Vz, 05.40.-a

We discuss the passage-time statistics of superradiant light pulses generated during the scattering of laser light from an elongated atomic Bose-Einstein condensate. Focusing on the early-stage of the phenomenon, we analyze the corresponding probability distributions and their scaling behaviour with respect to the threshold photon number and the coupling strength. With respect to these parameters, we find quantities which only vary significantly during the transition between the Kapitza Dirac and the Bragg regimes. A possible connection of the present observations to Brownian motion is also discussed.

Introduction

The interaction of optical fields with Bose-Einstein condensates (BECs) has attracted considerable interest over the last years. A rather intriguing phenomenon in this context is the generation of superradiant light pulses accompanied by condensate momentum side-modes, when far off-resonant laser light is scattered from an elongated BEC [1,2,3,4,5,6,7,8]. The process is initially driven by quantum fluctuations in the vacuum electromagnetic field (i.e., spontaneous Rayleigh scattering) [6]. After the first photons have been scattered, bosonic stimulation sets in and additional successive scattering events give rise to a rapid growth of the emitted radiation and of the side-mode populations. Moreover, due to the geometry of the condensate, the gain is largest when the scattered photons leave the condensate along its long axis, in the so called endfire modes. The fact that the process is triggered by quantum fluctuations is expected to be reflected in various macroscopic parameters of the superradiant (SR) light pulses (e.g., peak intensity, shape, etc), which should also fluctuate when the same experiment is repeated.

Throughout this work we focus on the passage-time statistics, which refer to the fluctuations in the time at which the pulse intensity at a given point in space reaches a certain reference intensity. The effect of quantum fluctuations is expected to be rather prominent for short times, where the observed intensity is relatively small. Given that the populations of the endfire modes and the side-modes grow fast in time, the evolution of the system becomes practically deterministic very soon [8] . For our purposes, it therefore suffices to focus on the early stage of the phenomenon, where the BEC is practically undepleted, only first-order atomic side-modes are significantly occupied, and the dynamics of the system are governed by linearized equations of motion. In view of the above arguments, the passage-time statistics at the early stage are expected to represent adequately also the delay-time statistics, which refer to the first maximum of the pulse that appears at later times.

Analogous theoretical as well as experimental studies have been performed in the context of conventional (Dicke) superradiance [9,[START_REF] Benedict | Super-Radiance[END_REF][START_REF] Gross | [END_REF], where excited atoms interacting with the vacuum of the electromagnetic field emit collectively; analytic expressions for the probability distribution of the passage and delay times have been obtained [START_REF] Gross | [END_REF]12]. In the system under consideration, however, condensed atoms may also scatter a photon from the endfire modes back into the laser mode, giving thus rise to recoiling atoms moving against the direction of the applied laser field (backward atomic sidemodes). This essential feature of superradiance from Bose-Einstein condensates is a key difference between this process and conventional superradiance. In this work we investigate how the presence of backward atomic side-modes affects the passage-time statistics, and examine whether the corresponding probability distributions usually discussed in conventional superradiance can also describe accurately the fluctuations in systems involving ultracold atoms. To the best of our knowledge, none of these questions has been addressed either theoretically or experimentally so far . In the following section we present the theoretical framework, while in Sec. 3 we analyze the passage-time statistics and their scaling behaviour. We conclude in Sec. 4.

Theoretical Framework

The system under consideration pertains to a BEC elongated along the z-axis, consisting of N atoms. A linearly polarized laser with frequency ω l = k l c, far detuned from the closest atomic transition by a value of δ, is illuminating the cloud along the x-axis. Due to the coherent nature of the condensate, successive Rayleigh scattering events are strongly correlated and lead to collective superradiant behaviour. As a result of the cigar-shape of the condensate, the gain is largest when the scattered photons leave the condensate along its long axis, traveling in the so called endfire modes with wave vectors k ≈ ±k l e z . A condensate atom can scatter a laser photon into the endfire modes, experiencing a recoil q ≈ (k l e xk). On the other hand, it can also scatter a photon from the endfire modes into the laser mode, in which case its momentum changes by q ≈ (-k l e x + k). These processes lead to the formation of two pairs of atomic side-modes, consisting of counterpropagating atoms with a narrow momentum spread (compared to k l ). Of course, atoms within these side-modes can also scatter photons, thereby acquiring higher momenta, but since we are interested in the early stage of the process, we consider only first order side-modes to be populated.

The Model

Neglecting any coupling between counterpropagating photonic endfire modes, the dynamics in the two halves of the elongated system become independent . We can thus focus on the endfire modes with k ≈ +ke z and consequently on the two sets of atomic side-modes (with central momenta q = ± (k l e x -ke z )) that are coupled to them. Due to the strong confinement along the xand y-axis, we can assume the transverse profiles of the matter field ψ ⊥ (x, y) to be well described by a classical function, independent of the z-coordinate [7]. Assuming a Fresnel number close to unity for the electromagnetic fields, we can apply the same approximation for the transverse part of the radiation field u ⊥ (x, y), effectively reducing the problem to one dimension.

We expand the field operators as

ψ(x, t) = ψ ⊥ (x, y) 1 j=-1 ψj (z, t)e ij(k l x-kz)-iω j t (1a) Ê(+) (x, t) = E 0 2 e y e i(k l x-ω l t) + u ⊥ (x, y) Ê(+) + (z, t)e y , (1b)
where ω ±1 = (k 2 l + k 2 )/2M with M the atomic mass and ω 0 = 0. The matter-wave operator is split up in three parts, describing the two side-modes (j = ±1) and the BEC at rest (j = 0), ψj (z, t) = e iω j t p∈∆ 0

e ipz √ L ĉ-jk+p (t), (2a)
where ∆ 0 is the interval (-k/2, k/2) in k-space, L is the length of the BEC and ĉp annihilates an atom with momentum p. Since the BEC at rest remains practically undepleted it can be treated as a time independent classical function and hence we can set ψ0 (z, t) ≡ ψ 0 (z). Similarly, we expand the endfire mode operator as

Ê(+) + (z, t) = ie i(kz-ωt) p∈∆ 0 ω k+p 2ε 0 1 √ L e ipz e iωt âk+p (t) ≈ ω 2ε 0 e i(kz-ωt) ê+ (z, t), ( 2b
)
where âp is the photon annihilation operator. The frequencies ω k+p are approximated by ω = k/c, the frequency of the scattered photons, and can therefore be taken out of the sum. This approximation is justified by the fact that dominant contributions to the sum come from momenta of order 1/L, which is several orders of magnitude smaller than k.

The one-dimensional field operators satisfy usual boson commutation relations

ψi (z 1 , t), ψ † j (z 2 , t) = δ ij δ ∆ (z 1 -z 2 ), (3a) ê+ (z 1 , t), ê † + (z 2 , t) = δ ∆ (z 1 -z 2 ), (3b) ê+
(z, t 1 ), ê † + (z, t 2 ) = 1 c δ ∆ (t 1 -t 2 ), (3c)
where δ ij denotes the Kronecker delta and δ ∆ (z) is a delta-like distribution with width of order 1/k. Since the fields vary slowly on length scales of the order of 1/k, we will approximate the δ ∆ (z) by actual delta functions. Inserting expansions (2.1) in the Maxwell-Schrödinger equations for the coupled matter-wave and electric fields and applying the slowly-varying-envelope approximation, we find the equations of motion for the operators

∂ ψ † +1 (ξ, τ ) ∂τ = iκê + (ξ, τ )ψ * 0 (ξ), (4a)
∂ ψ -1 (ξ, τ ) ∂τ = -iκê + (ξ, τ )ψ 0 (ξ) -2i ψ -1 (ξ, τ ), ( 4b
) ∂ê + (ξ, τ ) ∂τ + χ ∂ê + (ξ, τ ) ∂ξ = -i κψ 0 (ξ) ψ † +1 (ξ, τ ) +κ ψ -1 (ξ, τ )ψ * 0 (ξ) . (4c)
Here we have defined ψ -1 (ξ, τ ) = e -2iτ ψ-1 (ξ, τ ) and rescaled length and time to dimensionless units

ξ = k l z, τ = 2ω r t, (5)
where ω r = k 2 l /2M. Accordingly, the fields are rescaled as

ê+ (ξ, τ ) ≡ 1 √ k l ê+ (z, t) ψj (ξ, τ ) ≡ 1 √ k l ψj (z, t), (6)
and the speed of light becomes χ ≡ ck 2ωr . The effective 1D-coupling is given by κ

= g √ k l L/(2ω r ) with g = |d • e y | 2 E 0 2 δ ω 2 0 L dxdyu ⊥ (x, y)ψ 2 ⊥ (x, y). (7)
For a detailed derivation of (2.1) see [8].

Backwards recoiling atoms are a particular feature of superradiant Rayleigh scattering off condensates. The physical process underlying the backwards modes violates energy conservation by an amount ∆E 4 ω r . Thus, according to Heisenberg uncertainty principle, it can take place only for times shorter than a critical time /∆E which, in our units, is given by [1],

τ c = 0.5. ( 8
)
For such short pulses, and for sufficiently high power one typically observes an X-shaped pattern for the distribution of the atomic side-modes with the initial BEC in the center and the recoiling atoms moving both in and against the direction of the applied laser pulse (Kapitza-Dirac or strong-pulse regime). On the other hand, for weaker pulses with duration longer than τ c , the distribution of the side-modes exhibits a fan pattern, involving mainly forward recoiling atoms (Bragg or weak pulse regime). If we neglect the atomic backwards side-mode altogether, (2.1) become formally equivalent to descriptions of "conventional" superradiance in excited atomic gases [9].

Solutions to the equations of motion

We can use the Laplace transform to find exact solutions to the system (2.1) in terms of the operators evaluated at the boundary of their domain -i.e. at ξ = 0 and τ > 0 and at ξ > 0 and τ = 0. More details on this procedure are given in [8]. In the following we are interested in the radiation field only, whose time evolution is given by

ê+ (ξ, τ ) = τ 0 dτ ê+ (0, τ )F 0,0 (γ ξ,0 , τ -τ ) - iκ χ ξ 0 dξ ψ 0 (ξ ) ψ † +1 (ξ , 0)F 1,0 (γ ξ,ξ , τ) +ψ * 0 (ξ ) ψ -1 (ξ , 0)F 0,1 (γ ξ,ξ , τ) , (9)
where we have neglected retardation effects and introduced

γ ξ,ξ = κ 2 χ [ρ(ξ) -ρ(ξ )], (10)
with

ρ(ξ) = ξ 0 dξ |ψ 0 (ξ )| 2 . The functions F µ,ν (u, v) are defined as F µ,ν (u, v) = L -1 p→v e u/p e -u/(p+2i) p µ (p + 2i) ν , (11)
where L -1 p→v denotes the inverse Laplace transform. Explicit expressions for the functions F µ,ν (u, v) appearing in (9) are given in the appendix. They are combinations of Bessel functions and integrals thereof. For the numerical calculations, we assumed the BEC to consist of N = 10 6 Thomas-Fermi distributed 87 Rb atoms, such that ψ 0 = Θ(z)N 6(Lz -z 2 )/L 3 with L = 130µm. For the incoming laser we chose a rectangular profile and a wavenumber k l = 8.05 × 10 6 m -1 , which results in a dimensionless length of the BEC of Λ ∼ 1000. Coupling strengths are conveniently expressed in terms of the superradiant gain Γ = κ 2 N/χ, whose value separates the two regimes identified by experimental observations of superradiance from condensates [3]. Typically, the weak coupling regime is characterized by g ∼ 10 5 s -1 and Γ ∼ 1, while g ∼ 10 6 s -1 and Γ 1 in the strong coupling regime. In our calculations, we chose Γ = 1 and Γ = 100 for the two regimes.

Averaging over semiclassical trajectories

The presented model is closely related to the mean-field (MF) model presented in [7]. In fact, we arrive at the MF equations with the same approximations by replacing all operators in (2.1) by their expectation value. Due to the generality of the Laplace transform, the same replacement in (9) gives the solutions to the MF model. A detailed comparison between the two models is given in [8]. One of the problems of the MF model is the fact that it cannot describe the spontaneous onset of the superradiance process. In order to start the evolution, one has to "seed" one of the modes with some initial values. In the context of conventional superradiance, Haake et al. introduced the idea of averaging over many semiclassical "trajectories" to obtain quantum results [12]. The MF equations are initially seeded with random variables according to a particular distribution. To obtain a particular quantum expectation value, one has to average over various solutions for the corresponding semiclassical quantity. We have investigated the extension of this idea to superradiant Rayleigh scattering off BECs where, in contrast to conventional superradiance, backwards recoiling atoms are also present.

To see how this works, let us assume we want to calculate the normal-ordered nth order correlation function of the emitted light.

ê † + (ξ, τ ) n [ê + (ξ , τ )] n . (12)
From ( 9) and using the commutation relations as well as the fact that our initial state is the vacuum for all modes, we find the only non-vanishing expectation value involved to be ψ+1 (ξ (1) , 0) . . .

ψ+1 (ξ (n) , 0) ψ † +1 (ξ (n+1) , 0) . . . ψ † +1 (ξ (2n) , 0) , (13)
where the ξ (j) are integrated from 0 to ξ for j = 1, . . . , n and from 0 to ξ for j = n + 1, . . . , 2n. Using the commutation relations, correlation (13) reads

π n j=1 δ(ξ (j) -ξ (n+π(j)) ), (14)
where δ denotes the Dirac delta function and the sum runs over all permutations π of order n. Let us now seed the semiclassical model with [8]

ψ +1 (ξ, 0) = C ξ , (15a)
ψ -1 (ξ, 0) = 0, (15b)
e + (ξ, 0) = 0, ( 15c
)
where C ξ is a random, normally distributed complex variable, with zero mean and variance 1/ √ ∆ξ, with ∆ξ the spatial step of a numerical implementation ‡. The average of the product ψ +1 (ξ (1) , 0) . . . ψ +1 (ξ (n) , 0)ψ * +1 (ξ (n+1) , 0) . . . ψ * +1 (ξ (2n) , 0) over many trajectories (seedings) will effectively converge towards [START_REF] Chhikara | The Inverse Gaussian Distribution[END_REF]. Due to the formal equivalence of the semiclassical solutions to the quantum ones, the average of the product e * + (ξ 1 , τ) n e + (ξ 2 , τ) n will consequently converge towards the quantum expectation value (12). To find correlation functions of other operators, other seedings have to be used, which can be found in an analogous way and are summarized in [8].

As far as densities and populations are concerned, the convergence of the averaging procedure towards the quantum solution is fairly fast. In figure 1 we plot the photon density

n(ξ, τ ) = ê † + (ξ, τ )ê + (ξ, τ ) (16)
averaged over different numbers of trajectories compared to the quantum result.

Typically, one has to average over a couple of thousand trajectories, to obtain well converged, smooth density profiles. Clearly, the averaging procedure is also applicable to correlations of higher order. One can calculate any normal ordered correlation function of the system by averaging over a sufficiently large ensemble of MF trajectories, ‡ In our case 400 spatial gridpoints provided that the operators involved in the correlation function have the same seeding requirements. The convergence, however, becomes rapidly slower with every order added, such that higher order correlations will require a larger number of trajectories. Higher order correlations are related to higher order statistical moments of a particular observable, such as the variance or skew of a particular distribution.

Passage-Time Statistics

We proceed to apply the technique of averaging over semiclassical trajectories to the study of passage-time statistics. First, we will find a simple probability density function which matches the distributions. Following this, we will investigate the behaviour of the distributions with respect to scaling of the two most important parameters of such an experiment and find invariant quantities for both cases.

Distribution of Passage-Times

In view of Sec. 2.3, a single semiclassical trajectory represents a single experiment on superradiance off a given BEC [8,[START_REF] Gross | [END_REF]12]. The process is initially driven by quantum fluctuations; within the MF model this is simulated by seeding ψ +1 with Gaussian white noise. After the first spontaneous scattering events, bosonic stimulation leads to a rapid growth of the emitted radiation and of the side-mode populations. The signature of the quantum fluctuations which started the process, however, is imprinted onto the superradiant light pulses. One way to investigate quantum fluctuations in this context, is by studying the first passage-time statistics, which pertain to the distribution of times at which the intensity of the SR pulse reaches a certain reference intensity. Equivalently, one may look at the times for which the number of photons that have left the condensate volume, given by

N (τ ) = χ τ 0 dτ n(Λ, τ ), (17)
reaches a certain threshold N th . Performing a large number of semiclassical simulations, each one seeded with a new random ψ +1 , and following the evolution of N (τ ) for every trajectory, one can record the corresponding passage-times T . The quantum fluctuations responsible for the start-up of the system are translated into fluctuations of the passage times, whose statistics can provide information about the onset of the phenomenon. Moreover, due to the rapid growth of the mode populations and the deterministic evolution of the system, the passage-time statistics are expected to also reflect the delay-time statistics, which refer to the times at which the first maximum of the SR pulse occurs.

Based on these ideas, Haake et. al. investigated the statistics of passage-(and delay) times of SR light pulses obtained from a conventional superradiant system [12]. A semi-analytic expression was given for the distribution of the passage-times, which was found to be well approximated by the Gumbel probability density function (PDF) [START_REF] Gross | [END_REF]

P G (T ; µ, λ) = λ -1 e -T -µ λ e -exp(-T -µ λ ) . (18)
Their predictions about the standard deviation of the passage-times were found to be in good agreement with experimental data [13].

In contrast to conventional superradiant systems, in superradiance off BECs we also have the presence of backwards atomic side-modes. This is an essential feature of our system and can influence statistical properties of the scattered light, such as the the first-and second-order correlation functions [8]. Here, we will discuss its influence on the passage-time statistics. To this end, we performed a number of numerical simulations using the semiclassical equations of motion. For each trajectory, the equations were seeded in accordance with (2.3), and we followed the time-evolution of the number of outgoing photons N (τ ). As soon as the threshold N th was reached, we recorded the corresponding time τ = T , and started a new simulation with new random seeding. The histograms depicted in figure 2 describe the distribution of the recorded times after a large number of simulations.

In the weak pulse regime, the distributions are narrower (i.e., the variance of the distributions is smaller) than in the strong pulse regime. An intuitive understanding for this broadening can be obtained by recalling that in the strong pulse regime, the growth of the photonic endfire modes is slowed down by the presence of backwards recoiling atoms, as each backscattered atom removes a photon from the endfire modes. Therefore, with respect to the scaled time Γτ , the spread of solutions for N (τ ) corresponding to different seedings, is larger in the strong than in the weak pulse regime. This is illustrated in figure 3, where the evolution of N (τ ) is plotted for two different seedings in the strong-and the weak-pulse regimes.

All the distributions exhibit positive skew, which limits the PDFs appropriate to model the data. Several two-parameter, positively skewed PDFs were fitted to the data according to a least-square fit and assessed by calculating the residual errors as well as comparing mean, variance and skewness of the distributions to the ones of the data. While the Gumbel PDF of conventional superradiance gives a very satisfying fit, it is for all parameter values studied inferior § to a fit with the inverse Gaussian PDF

P IG (T ; µ, λ) = λ 2πT 3 e -λ (T -µ) 2 2µ 2 T , (19)
which, as a further advantage, does not support negative T . If a random variable X is distributed according to (19), the mean and variance of X are given by

X = µ, var(X) = µ 3 /λ. ( 20
)
In figure 2 we show two fits to the numerical data, corresponding to the Gumbel and the inverse Gaussian PDFs. § The assessment of fitted PDFs was performed by calculating the sum of residual errors as well as by checking the agreement of the PDFs mean, variance, skew and kurtosis with the ones of the raw data. Dark(red) histograms show the difference between the two cases, in the two regimes. The curves give least square fits with inverse Gaussian (solid) and Gumbel (dashed) PDFs.

We now turn to the comparison of the passage-time distributions for superradiance off BECs, to the ones obtained for conventional superradiance. As discussed in section 2.1, the latter system can be simulated within our model by neglecting the backwards atomic side-modes in (2.1). The corresponding distributions of the recorded passagetimes are depicted in figures 2 (b) and (d). From the difference of the two histograms, also plotted in figure 2, we see that in the weak pulse regime (Γ = 1), the two distributions are virtually the same [compare figures 2 (a) and (b)]. For weak couplings, the number of backwards recoiled atoms is not big enough to change dramatically the passage time statistics. In the strong pulse regime, however, the distribution for conventional superradiance is narrower than the one for superradiance off BECs, while it is also shifted towards smaller passage times [compare figures 2 (c) and (d)]. Finally, in the weak pulse regime, fits with either of the aforementioned two PDFs agree very well with the distributions of the numerical data [see figure 2 contrary, in the strong-pulse regime, and if backwards recoiled atoms are present, we find noticeable discrepancies between the two fits, particularly for small T . They become more prominent for small thresholds N th .

Scaling with coupling strength

We have seen that our equations of motion can describe conventional superradiance if we set all the couplings to backward atomic side-modes to zero. The resulting reduced set of equations is also expected to be a good approximation in the weak-coupling regime of superradiance off condensates, provided that the time scales of interest are such that the backward atomic side-modes are scarcely populated. Assuming furthermore a homogeneous condensate, analytic solutions for the operators can be found, which moreover are a function of the scaled time Γτ rather than just the time alone [8]. If one multiplies inverse-Gaussian distributed random numbers with a constant factor, the resulting numbers are again inverse-Gaussian distributed, with the mean µ and the shape parameter λ rescaled by the same number. Hence the scaling behaviour is also reflected in the distributions, which can be seen by comparing the distributions of figure 2 (b) and (d), where a scaling factor or 100 is clear. We can therefore extend our previous observation on the scaling behaviour of the passage-time distributions for conventional superradiance, to the context of superradiance off condensates and define a more general "scaled" version of the inverse Gaussian PDF, which describes accurately passage-time distributions for various couplings within the weak-coupling regime,

P WP (τ ; µ, λ) = P IG (τ, µ WP /Γ, λ WP /Γ). (21)
For any N th , the two open parameters µ WP and λ WP are constants characterizing the weak-pulse regime. For N th = 10, their value was estimated to µ WP ≈ 3.2 and λ WP ≈ 16.

In the case of conventional superradiance, it is known that the mean of the passage times scales as a constant divided by the atomic decay rate, while their variance scales as another constant divided by the decay rate squared [START_REF] Gross | [END_REF]. In our model, the quantity corresponding to the decay rate is the superradiant gain Γ, which is directly proportional to the superradiant scattering rate. From ( 20) and ( 21) we therefore find that the inverse Gaussian distribution recovers the same scaling laws as found in [START_REF] Gross | [END_REF].

If we increase the coupling strength and leave the weak-coupling regime, we lose this scaling behaviour, as can be seen by comparing the histograms (a) and (c) in figure 2. Surprisingly, however, we recover (21) deep within the strong pulse regime, albeit with new constant parameters, which for N th = 10 are given by µ SP ≈ 5 and λ SP ≈ 20. To gain further insight into the change of the scaling behaviour of the passage-time distribution with the coupling strength, we calculated numerically the distributions H(ΓT, N th ), of the rescaled passage-times ΓT , for various values of Γ with Γ low ≤ Γ ≤ Γ high , and various thresholds N th . The upper and lower values of Γ were chosen as Γ low = 0.1 and Γ high = 100, so that both regimes are covered. Subsequently, for a given N th , we calculated the normalized overlap max i∈{low,up}

∞ 0 H(ΓT, N th )H(Γ i T, N th )dT 2 ∞ 0 H(ΓT, N th ) 2 dT ∞ 0 H(Γ i T, N th ) 2 dT , (22)
and the corresponding density plot is depicted in figure 4. While the parameters µ WP/SP and λ WP/SP are independent of the coupling, they are nontrivial functions of the photon threshold N th , and their numerical dependence is also depicted in figures 4 (a,c). Quite obviously, the scaling law within the two regimes is very well obeyed by the system. We can therefore conclude that in the early stages of superradiant scattering off condensates, the distribution of passage times in either regime is universally well modelled by the inverse Gaussian PDF and two constant parameters. The transition between the two scaling behaviours begins in the region, where the parameters N th and Γ match up in a way such that T ≈ τ c . The distributions in this transient regime do not differ qualitatively from the ones presented in figure 2, but one cannot define a universal scaling law for them. As can be seen in figure 4 (b), however, this transient regime (dark region) in the parameter space is rather small compared to the ones where the system is unambiguously found in either of the two regimes (white regions). As a final remark, it should be noted that for any threshold N th , the ratio of mean and shape parameter of the resulting distribution is independent of Γ, which by virtue of (20) implies that the ratio of the variance and the squared mean of passage times is independent of Γ.

Scaling with threshold photon number

Having investigated the behaviour of the passage time distributions for varying Γ, it is natural to ask how they behave with a constant coupling but varying the other significant parameter of the system, the photon threshold N th . For any choice of Γ, we find again a simple combination of shape parameter and mean, which remains constant independent of N th . It is the ratio

s = λ µ 2 = T var(T ) , (23)
where the second equality is a simple consequence of (20). Figure 5 (a) depicts a contour plot of the value of s for different couplings and photon thresholds, and reveals the independence of this value of the coupling Γ. Since these results are based on the analysis of a large ensemble of random trajectories, a certain variation is to be expected, and the values of s as a function of Γ are obtained by averaging s over N th . Moreover, The parameter s, averaged over N th , as a function of the coupling Γ (crosses) and the fits (24) in the weak and strong pulse regimes (lines). For every threshold and coupling 10 4 trajectories were evolved.

for each regime, the quantity s has a very simple linear dependence on the coupling Γ, which is illustrated in figure (5) (b). The functions plotted are the fits s WP (Γ) = 1.56Γ (weak pulse regime)

s SP (Γ) = 0.71Γ + 2.86 (strong pulse regime).

Note that the equation for the strong pulse regimes is valid for large values of Γ, such that the offset becomes negligable and s SP is directly proportional to Γ. Hence this relation is consistent with equation 21. Moreover, also in this case the transition between the two fits (24), takes place in a regime where the parameters result in times around critical time τ c = 0.5. Random variables distributed according to the inverse Gaussian PDF and with this very ratio as a constant, also appear in different contexts. Consider the stochastic process

X 0 = 0 (25a) X t = νt + σB t , (25b)
where B t is the standard Brownian motion and ν > 0 is a positive, constant drift. The times at which X t crosses a certain limit α > 0 are distributed according to the inverse Gaussian distribution [START_REF] Chhikara | The Inverse Gaussian Distribution[END_REF] T α = inf{0 < t|X t = α}

∼ P IG T α ; µ = α ν , λ = α 2 σ 2 . ( 26
)
The ratio of shape parameter and squared mean is independent of α, which is analogous to s being independent of the photon threshold N th and we can thus make a connection between Brownian motion and superradiance off BECs. In particular, we can construct a stochastic system analogous to (3.3), but in terms of T and s. Note that the particular choice of corresponding parameters is not unique. From (26) we find that α/ν must correspond to the quantum expectation value T , while

s = ν 2 σ 2 . ( 27
)
We can therefore set the drift ν to unity and find T = α and σ = 1/ √ s. Hence, the times the stochastic system X 0 = 0, (28)

X τ = τ + 1 √ s B τ , (29)
will cross X τ = T will reproduce the distribution of passage-times for a superradiant BEC emitting N th photons. This result by no means implies that superradiance is in fact a physical realization of Brownian motion. In contrast to Brownian motion, superradiance off condensates is a deterministic evolution from random initial conditions. After its initial stage, the system can be described by deterministic equations. This is never the case for Brownian motion.

Conclusion

We have used random seedings to investigate the passage-time statistics of superradiant light pulses emitted from an elongated atomic condensate. In general, the distributions of the passage times are very well approximated by the inverse Gaussian probability density function. The distributions in the Kapitza-Dirac regime, however, appear wider than in the Bragg regime, while they are also shifted towards larger passage times. This is due to the presence of backwards recoiling atoms in the system, which in the Kapitza-Dirac regime lead to a slow growth of end-fire modes. For both regimes the mean of the distributions is inversely proportional to the superradiant gain Γ, whereas their variances scale with Γ -2 . The corresponding proportionality constants are different in the two regimes and for each regime, the ratio of shape parameter λ and mean µ is independent of Γ. For a small region in parameter space, the system is not in either of the two regimes but finds itself in a transient one, where the distribution of the passage times changes smoothly between the two reference (outer most) distributions. This transient region is characterized by timescales around the critical time, at which energy conservation suppresses backwards recoiling atoms. For constant coupling and varying threshold, the ratio λ/µ 2 is constant. The latter behaviour leads to a phenomenological correspondence between superradiance from BECs and Brownian motion with drift. We focused at the very early stage of the phenomenon, but due to the deterministic evolution of the system the observed behaviour can be expected to hold for later times and consequently the distributions are expected to also describe the delay-times of a superradiant pulse from a condensate. Finally we want to stress that the presented results are accessible to experimental verification. Similar experimental studies have been performed for conventional superradiance [START_REF] Gross | [END_REF]13].

Figure 1 .

 1 Figure 1. (Colour online) Example of the photon density within the BEC volume obtained through averaging over randomly seeded MF solutions. The quantum solution is shown in solid black, while averages over 20 and 2000 trajectories are shown in red dashed and blue dotted lines. Parameters: Γ = 1, and τ = 2.

Figure 2 .

 2 Figure 2. (Colour online) Passage-time statistics of light pulses in superradiance from BECs. Histograms (a) and (c) show the distribution of the times at which the number of photons that have left the BEC reaches the threshold N th = 10, for (a) Γ = 1, and (c) Γ = 100. The corresponding histograms with neglected backwards recoiling atoms (equivalent to conventional superradiance), are plotted with inverted ordinate (b,d).Dark(red) histograms show the difference between the two cases, in the two regimes. The curves give least square fits with inverse Gaussian (solid) and Gumbel (dashed) PDFs.

Figure 3 .

 3 Figure 3. (Colour online) Evolution of the number of photons in the SR light pulse for two different seeds, as a function of scaled time in the strong and weak coupling regime. The vertical lines indicate the wider spread due to the slower growth in the strong coupling regime.

Figure 4 .

 4 Figure 4. (Colour online) (a) Shape parameters of the inverse gaussian PDFs as a function of photon threshold N th . (b) Density plot of the overlap (22), for various couplings Γ and thresholds N th . White corresponds to a value of unity. Note the logarithmic scale for Γ. The contour lines indicate the points for which the expectation time for T is equal to 0.5,1 and 2 respectively. Each histogram contains 10 4 trajectories. (c) Means µ of the scaled passage times as a function of photon threshold N th .

Figure 5 .

 5 Figure 5. (Colour online) (a) Contour plot of the value of s as defined in (23). (b)The parameter s, averaged over N th , as a function of the coupling Γ (crosses) and the fits (24) in the weak and strong pulse regimes (lines). For every threshold and coupling 10 4 trajectories were evolved.

Appendix A. Inverse Laplace Transforms

Using elementary techniques for Laplace inversion, one finds

with J i and I i the ith Bessel function of the first kind and the ith modified Bessel function respectively [8].


O M Braun

M Peyrard
email: michel.peyrard@ens-lyon.fr

Dependence of kinetic friction on velocity: Master equation approach

Keywords: numbers: 81.40.Pq, 46.55.+d, 61.72.Hh

We investigate the velocity dependence of kinetic friction with a model which makes minimal assumptions on the actual mechanism of friction so that it can be applied at many scales provided the system involves multi-contact friction. Using a recently developed master equation approach we investigate the influence of two concurrent processes. First, at a nonzero temperature thermal fluctuations allow an activated breaking of contacts which are still below the threshold. As a result, the friction force monotonically increases with velocity. Second, the aging of contacts leads to a decrease of the friction force with velocity. Aging effects include two aspects: the delay in contact formation and aging of a contact itself, i.e., the change of its characteristics with the duration of stationary contact. All these processes are considered simultaneously with the master equation approach, giving a complete dependence of the kinetic friction force on the driving velocity and system temperature, provided the interface parameters are known.

I. INTRODUCTION

Almost three centuries ago Charles Coulomb (1736-1806) discovered that kinetic friction does not depend on the sliding velocity [START_REF] Dowson | History of tribology[END_REF]. Later, more careful experiments showed that this law is only approximately valid [START_REF] Persson | Sliding friction: Physical principles and applications[END_REF][START_REF] Robbins | Handbook of modern tribology[END_REF][START_REF] Robbins | Jamming and rheology: Constrained dynamics on microscopic and macroscopic scales[END_REF][START_REF] Persson | [END_REF][6][7]. Friction does depend on the sliding velocity, but this dependence is far from universal: some measurements find an increase when velocity increases, while others find a decay [8][9][10] or even a more complex non-monotonous behavior [START_REF] Persson | Sliding friction: Physical principles and applications[END_REF]. A logarithmic dependence, often quoted, has been found for two extreme scales, friction at the tip of an AFM (see for instance [8,9,[11][12][13]) or at the scale of a fault in the earth crust [10], but it is often only approximate and observed in a fairly narrow velocity range. Therefore understanding the velocity dependence of kinetic friction is still an open problem, and what makes it difficult is that several phenomena contribute, the thermal depining of contacts, their aging, and the delay in contact formation.

In this study we investigate the velocity dependence of friction with a model that includes these three contributions and makes minimal assumptions on the actual mechanism of friction so that it can be applied at many scales provided the system involves multi-contact friction. Our aim is to elucidate the respective role of these three contributions to the velocity dependence of friction and to provide analytical treatments in some limits, or simple numerical approaches that allow the investigation of a velocity range that may span many orders of magni-tude.

At the most fundamental level multi-contact friction can be described as resulting a succession of breaking and formation of local contacts which possess a distribution of breaking thresholds. This viewpoint was first applied to describe earthquakes [14,15] and then adopted to friction by Persson [16].

We recently developed a master equation (ME) approach to describe the breaking and attachment events [17,18]. It splits the analysis in two independent parts: (i) the calculation of the friction force, given by the master equation provided the statistical properties of the contacts are known, and (ii) the study of the properties of the contact themselves, which is system dependent. This method is very general and allows us to calculate the velocity dependence of friction, which results from the interplay of two concurrent processes. First, at a nonzero temperature thermal fluctuations allow an activated breaking of contacts which are still below their mechanical breaking threshold. This phenomenon leads to a monotonic increase of the friction force F with the velocity v. Second, the aging of contacts [19][20][21] leads to a decrease of the friction force with velocity. It includes two processes: the delay in contact formation, i.e., time lag between contact breaking and re-making [18,[20][21][22][23], and the aging of a contact itself, i.e., the change of its characteristics with the time of stationary contact. To incorporate the latter effect, the master equation must be completed by an equation for the evolution of static thresholds.

In earlier studies [17,18] we considered thermal and aging effects separately, to set up the method. However, to relate the results to experiments, both contributions must be taken into account simultaneously. This is the aim of the present paper, which is organized as follows. Section II is a brief review of the mas-ter equation approach. Section III discusses temperature effects. Whereas our earlier work [18] focussed on time-dependent phenomena to analyze stick-slip, here we concentrate on the steady-state case (constant velocity). This allows us to proceed further and derive explicit expressions for the influence of temperature. Then Sec. IV introduces the second effect, the aging of the contacts. It first summarizes the method introduced earlier and its main results [18], which only considered the T = 0 case, and then studies the combined influence of aging and temperature fluctuations. Section V adds the influence of the delay in contact formation after breaking, to get the full picture, allowing us to compute the velocity dependence of friction. Section VI discusses all those results in the context of experimental data. The difficulty to apply the theory to actual experiments is to properly assess the values of the parameters that enter in the theoretical expressions, and not simply try to fit experimental curves, which would not be very significant owing to the number of parameters which are involved. Therefore Sec. VI focusses on this assessment. Finally, Sec. VII concludes the paper with a discussion of perspectives for its further development.

II. MASTER EQUATION

The earthquake (EQ) model is the most generic model for friction due to multiple contacts at an interface. The sliding interface is treated as a set of N c "contacts" which deform elastically with the average rigidity k. The contacts represent, for example, asperities for the interface of rough surfaces [24], or patches of lubricant or its domains ("solid islands" [25]) in the case of lubricated friction. The ith contact connects the slider and the bottom substrate through a spring of elastic constant k i . When the slider moves, the position of the contact point changes, and the contact's spring elongates or shortens, so that the slider experiences the force -F = f i from the interface, where f i = k i x i and x i (t) is the spring length. The contacts are coupled frictionally to the slider. Namely, as long as the force |f i | is below a certain threshold f si (corresponding to the onset of plastic flow of the entangled asperity, or to local shear-induced melting of the boundary lubrication layer), this contact i moves together with the slider. When the force exceeds the threshold, the contact breaks, and then re-attaches again in the unstressed state after some delay time τ . Thus with every contact we may associate the threshold value f si , which takes random values from a distribution P c (f ) having a mean value f s . The spring constants are related to the threshold forces by the relationship k i = k (f si /f s ) 1/2 , because the value of the static threshold is proportional to the area A i of the given contact, while the transverse rigidity k i is proportional to contact's size, k i ∝ √ A i . When a contact is formed again (re-attached to the slider), new values for its parameters have to be assigned.

Rather than studying the evolution of the EQ model by numerical simulation it is possible to describe it analytically [17,18]. Let P c (x) be the normalized probability distribution of values of the thresholds x si at which contacts break; it is coupled with the distribution of threshold forces by the relationship P c (x) dx = P c (f ) df . To describe the evolution of the model, we introduce the distribution Q(x; X) of the stretchings x i when the bottom of the solid block is at a position X. Let us consider a small displacement ∆X > 0 of the bottom of the sliding block.

It induces a variation of the stretching x i of the contacts which has the same value ∆X for all contacts (here we neglect the elastic deformation of the block). The displacement X leads to three kinds of changes in the distribution Q(x; X): first, there is a shift due to the global increase of the stretching of the asperities; second, some contacts break because their stretching exceeds the maximum value that they can withstand; and third, those broken contacts form again, at a lower stretching, after a slip at the scale of the asperities, which locally reduces the tension within the corresponding asperities. These three contributions can be written as a master equation for Q(x; X):

∂ ∂x + ∂ ∂X + P (x) Q(x; X) = R(x) Γ(X) , (1)
where P (x) ∆X describes the fraction of contacts that break when the slider position changes from X to X + ∆X. At zero temperature P (x) is coupled with the threshold distribution P c (x) by the relationship [17,18]

P (x) = P c (x)/J c (x) , J c (x) = ∞ x dξP c (ξ) . (2)
The function Γ(X) in Eq. ( 1) describes the contacts that form again after breaking,

Γ(X) = ∞ -∞ dξ P (ξ) Q(ξ; X) (3)
(the delay time is neglected at this stage), and R(x) is the (normalized) distribution of stretchings for newborn contacts. Then, the friction force is given by

F (X) = N c k ∞ -∞ dx x Q(x; X) . (4)
The evolution of the system in the quasi-static limit where inertia effects are neglected shows that, in the long term, the initial distribution approaches a stationary distribution Q s (x) and the total force F becomes independent of X. This statement is valid for any distribution P c (x) except for the singular case of P c (x) = δ(x -x s ).

In the present work we concentrate on the steady state (smooth sliding). In what follows we use R(x) = δ(x) for simplicity. The steady-state solution of Eq. (1) is

Q(x) = Θ(x)E P (x)/C[P ] , (5)
where Θ(x) is the Heaviside step function (Θ(x) = 1 for x ≥ 0 and 0 otherwise),

E P (x) = exp[-U (x)], U (x) =
x 0 dξ P (ξ), and

C[P ] = ∞ 0 dx E P (x). Note also that, in the steady state, Γ = 1/C[P ] , (6) because
∞ 0 dξ P (ξ)E P (ξ) = ∞ 0 dU e -U = 1.
The distribution P c (f ) can be estimated for the contact of rough surfaces [18,24] as well as for the contact of polycrystal substrates [18,26]: its general shape may be approximated by the function

P c (f ) ∝ f n exp(-f /f * ) , (7)
where n ≥ 0 depends on the nature of the interface.

Then, the distribution P c (x) can be related to the distribution P c (f s ) of the static friction force thresholds of the contacts. If a given contact has an area A, then it is characterized by the static friction threshold f s ∝ A and the (shear) elastic constant k ∝ √ A (assuming that the linear size of the contact and its height are of the same order of magnitude, see [16] and Appendix A in Ref. [18]). The displacement threshold for the given contact is

x s = f s /k, so that f s ∝ x 2 s , or df s /dx s ∝ x s . Then, using P c (x s ) dx s = P c (f s ) df s , we obtain P c (x s ) ∝ x s P c [f s (x s )], or P c (x) ∝ x 1+2n exp(-x 2 /x 2 * ) , (8)
where x * may be estimated from experiments as N c kx * ≈ F s . In the SFA/B (surface force apparatus/balance) experiments, where the sliding surfaces are made of mica, the interface may be atomically flat over a macroscopic area. But even in this case the lubricant film cannot be ideally homogeneous throughout the whole contact area -it should be split into domains, e.g., with different orientation, because this will lower the system free energy due to the increase of entropy. Domains of different orientations have different values for the thresholds f si , i.e., they play the same role as asperities in the contact of rough surfaces.

For the normalized distribution of static thresholds given by Eq. ( 8) with n = 1,

P c (x) = (2/x * ) u 3 e -u 2 , where u ≡ x/x * , (9)
we can express the steady-state solution of the master equation analytically. In this case

J c (x) = (1 + u 2 ) e -u 2 , (10)
so that at zero temperature we have

P (x) = (2/x * ) u 3 /(1 + u 2 ) , (11)
U (x) = u 2 -ln(1 + u 2 ) , (12)
E P (x) = J c (x) = (1 + u 2 ) e -u 2 , (13)
C[P ] = x * /C 0 , where C 0 = 4 3 √ π ≈ 0.752, (14)
Q(x) = (C 0 /x * )(1 + u 2 ) e -u 2 , u ≥ 0 (15)
and the kinetic friction is

f k ≡ F k / (N c k) = f k0 ≡ C 0 x * . (16)
The ME formalism described above can be extended to take into account various generalizations of the EQ model, such as temperature effects and contact aging, which are examined in the following sections.

III. NONZERO TEMPERATURE

Temperature effects enter in the ME formalism through their effect on the fraction of contacts that break per unit displacement of the sliding block, P (x), because thermal fluctuations allow an activated breaking of any contact which is still below the threshold [16,18,[20][21][22][23]. For a sliding at velocity v so that X = vt, the thermally activated jumps can be incorporated in the master equation, if we use, instead of the zero-temperature breaking fraction density P (x), an expression P T (x) defined by (see [18])

P T (x) = P (x) + H(x) , (17)
where the temperature contribution is given by

H(x) = ω v e kx 2 /2kB T ∞ x dξ P c (ξ) e -kξ 2 /2kB T (18)
for "soft" contacts or by

H(x) = ω v ∞ x dξ P c (ξ) 1 - x ξ 1 2 e -kξ 2 (1-x ξ ) 3 2 /2kB T (19
) in the case of "stiff" contacts which have a deep pinning potential so that their breaking only occurs with a significant probability when their stretching is close to the threshold. Here ω is the attempt frequency of contact breaking, ω ∼ 10 10 s -1 according to Refs. [16,21].

For concreteness, in what follows we assume that the contacts are soft, Eq. ( 18), and we select n = 1 in Eq. ( 8), so that P c (x) is given by Eq. ( 9).

At a nonzero temperature the total rate of contact breaking, Eq. ( 17), is equal to P T (x) = P (x)+(ω/v) h(x), where for the soft contacts

h(x) = 1 + (1 + b) u 2 (1 + b) 2 e -u 2 (20) with b(T ) = kx 2 * 2k B T . ( 21
)
The condition b = 1 defines a crossover temperature

k B T * = 1 2 kx 2 * . (22)
Then, a straightforward integration gives the function

U T (x) = x 0 dξ P T (ξ) = U (x) + ∆U (x), where ∆U (x) = S 0 (v, T ) erf(u) -S 1 (T ) u e -u 2 (23)
with

S 0 (v, T ) = ωx * C 0 v (1 + b/3) (1 + b) 2 (24)
and

S 1 (T ) = C 0 2 (1 + b) (1 + b/3) . ( 25
)
The coefficient S 1 (T ) weakly changes with temperature from S 1 (0

) = 2/ √ π ≈ 1.128 to S 1 (∞) = C 0 /2 ≈ 0.
376. On the other hand, the coefficient S 0 (v, T ) determines whether the effect of temperature is essential or not. The temperature-induced breaking of contacts is essential at low driving velocities only, when S 0 (v, T ) ≫ 1. Thus, the equation S 0 (v * , T ) = 1 defines the crossover velocity:

v * (T ) = ωx * C 0 [1 + b(T )/3] [1 + b(T )] 2 . ( 26
)
We see that v * monotonically increases with temperature as v * (T ) ≈ 0.443 ωx * T /T * at T ≪ T * and approaches the maximal value v * ≈ 1.33 ωx * at T ≫ T * .

Then, E P T (x) = e -UT (x) = (1 + u 2 ) e -u 2 e -∆U(x) , and we can find the kinetic friction force:

f k (v, T ) = ∞ 0 dx xE P T (x) ∞ 0 dx E P T (x) . (27)
At a low driving velocity, v ≪ v * , we may put ∆U (x) ≈ S 2 u, where

S 2 (v, T ) = ωx * v (1 + b) 2 , (28)
and Eq. ( 27) leads to

f k ≈ x * /S 2 = (v/ω)(1 + b) 2 . ( 29
)
A linear dependence of the kinetic friction on the driving velocity at low velocities corresponds to the creep motion due to temperature activated breaking of contacts and was predicted in several earlier studies [16,[21][22][23], although our approach allowed us to derive it rigorously.

The dependence (29) may be interpreted as an effective "viscosity" of the confined interface: At a high velocity, v ≫ v * , when e -∆U(x) ≈ 1 -

η = f k v = 1 ω 1 + kx 2 * 2k B T 2 . ( 30
)
S 2 u + 1 2 (S 2 u) 2 , we obtain C ≈ 3 √ π/4 -S 2 + (5 √ π/16)S 2 2 so that f k ≈ f k0 1 -C 1 S 2 + C 2 S 2 2 , (31)
where 31) agrees qualitatively with that found by Persson [16] in the case of b ≫ 1.

C 1 = 5 √ π/8-C 0 ≈ 0.356 and C 2 = 16/9π-1/3 ≈ 0.233. Equation (
Approximate expressions [START_REF] Klein | [END_REF]31) together with the numerical integration of Eq. ( 27) are presented in Fig. 1. Also we showed a logarithmic fitting which operates in a narrow interval of velocities near the crossover velocity only. Persson [16] showed that the logarithmic dependence may be obtained analytically, only if the P c (x) distribution has a sharp cutoff at some x = x s as, e.g., in simplified versions of the EQ model with

P c (x) = δ(x -x s ).
Although we cannot obtain analytical results for the stiff contacts, we calculated the f k (v) dependences numerically (see Fig. 2), which shows that the effect remains qualitatively the same.

IV. AGING OF CONTACTS

The aging of contacts was considered in our work [18] where, however, we ignored the temperature-induced breaking of contacts discussed above in Sec. III. When aging is taken into account the master equation for Q(x, X) must be completed by an equation for the evolution of P c (x), which in turns affects P (x). Let the newborn contacts be characterized by a distribution P ci (x), while at t → ∞, due to aging the distribution P c (x) approaches a final distribution P cf (x). If we assume that the evolution of P c (x) corresponds to a stochastic process, then it should be described by a Smoluchowsky equation

∂P c ∂t = D LP c , where L ≡ ∂ ∂x B(x) + ∂ ∂x , (32)
in which the "diffusion" parameter D describes the rate of aging, B(x) = d Ū (x)/dx, and the "potential" Ū (x) determines the final distribution, P cf (x) ∝ exp -Ū (x) , so that we can write

B(x) = -[dP cf (x)/dx] /P cf (x) . (33)
However, because the contacts continuously break and form again when the substrate moves, this introduces two extra contributions in the equation determining ∂P c /∂X in addition to the pure aging effect described by Eq. (32): a term P (x; X)Q(x; X) takes into account the contacts that break, while their reappearance with the threshold distribution P ci (x) gives rise to the second extra term in the equation. Thus, the evolution of P c is described by the equation

∂P c (x; X) ∂X -D v LP c (x; X) + P (x; X) Q(x; X) = P ci (x) Γ(X) , (34)
where D v = D/v, and v = dX(t)/dt is the driving velocity. Finally, we come to the set of equations (1-3, 34).

For the steady-state regime, Eq. (34) reduces to

D v C[P ] LP c (x) = P (x) E P (x) -P ci (x) , (35)
where we used Eqs. ( 5) and (6). Taking also into account the identity P (x) E P (x) = P c (x) [18], we finally come to the equation

D v C[P ] LP c (x) = P c (x) -P ci (x) . (36)
It was shown [18] that the kinetic friction monotonically decreases with the driving velocity as F k (v) - (black solid), 10 -2 (magenta dash-dotted), 10 -3 (blue shortdotted), and D = 0 (red dotted curve). k = 1, ω = 1, kBT = 1; the initial and final Pc(x) distributions are given by Eq. ( 9) with x * i = 0.1 and x * f = 1 correspondingly. Dashed curves shows the dependences at T = 0.

F k (0) ∝ -v/D in the low-velocity limit and F k (v) -F k (∞) ∝ D/v in the high-velocity case. One may expect that at low velocities this decreasing will compensate the friction increasing due to temperature induced jumps. The problem, however, is more involved. When the temperature effects are incorporated, Eq. (36) for the function P c (x) in the steady state takes the form

D v C[P T ] LP c (x) = P c (x) -P ci (x) . ( 37
)
Numerical solutions of Eq. (37) are presented in Figs. 3 and4a: the initial increase of the kinetic friction F with the driving velocity v due to the temperature activated breaking of contacts is followed by the decrease of F due to contacts aging. Figure 4b shows also the dependence of the effective "viscosity" η = f k /v on the driving velocity v. It is constant at low velocity and then decreases; the latter may be approximately fitted by a power law η(v) ∝ v -α with the exponent α changing from 1.5 to 1 as D decreases.

Using the definition (32) of the operator L and Eq. (33) for the function B(x), the l.h.s. of Eq. (37) may be rewritten as

D v C LP c (x) = D v C d dx P cf (x) d dx P c (x) P cf (x) , (38)
while the r.h.s. of Eq. (37) may be presented as

P c (x) -P ci (x) = - d dx [J c (x) -J ci ] , (39)
where J ci (x) = we can find the first integral of Eq. ( 37):

D v C[P T ] P cf (x) d dx P c (x) P cf (x) = J ci (x) -J c (x) . (40)
Integration of Eq. ( 40) leads to an integral equation for the function P c (x):

P c (x) = P cf (x) 1 + v DC[P T ] x 0 J ci (ξ) -J c (ξ) P cf (ξ) dξ .
(41) Substituting J c (x) ≈ J cf (x) = ∞ x dξ P cf (ξ) into the r.h.s. of Eq. ( 41), one may analytically find the lowvelocity behavior of the kinetic friction, for example, the decrease of f k with v for T = 0. At a nonzero temperature, however, aging does not affect the low-velocity behavior [START_REF] Klein | [END_REF] and only reduces the interval of velocities where Eq. ( 29) is valid, as demonstrated in Figs. 3 and4. Indeed, at v → 0 and T > 0 the main contribution to P T (x) comes from the function H(x) ∝ ω/v in Eq. ( 17), which only weakly depends on P c (x).

The limit D v = D/v → 0 may be studied with the help of Eq. ( 37) by substituting P c (x) ≈ P ci (x) into its left-hand side. For the function (9), this approach gives

P c (x) -P ci (x) ≈ - 16D v C[P ] x * i x 2 if u i e -u 2 i 1 - 1 2 u 2 i ( 42
)
and

J c (x) -J ci (x) ≈ 4D v C[P ] x 2 if u 4 i e -u 2 i , ( 43
)
where u i = x/x * i and (black dash-dotted) and 10 -3 (red dotted curve) as compared with approximate expressions (thin solid curves). k = 1, ω = 1, kBT = 1; the initial and final Pc(x) distributions are given by Eq. ( 9) with x * i = 0.1 and x * f = 1 correspondingly.

1 x 2 if = 1 x 2 * i - 1 x 2 * f . ( 44
Then, taking the corresponding integrals, we obtain to first order in v -1

C ≈ x * i C -1 0 -S 2 + C 12 S 3 (45)
and

f k ≈ x 2 * i C 1 -(C 0 + C 1 )S 2 + C 12 S 3 , (46)
where

S 3 (v) = 4D v C[P ]/x 2 if , (47)
C 12 ≈ 1.324 and C 12 ≈ 1.789 are numerical constants. A comparison of the exact and approximate expressions is shown in Fig. 5.

V. DELAY IN CONTACT FORMATION

Finally, let us take into account the delay in contact formation following the work of Schallamach [22]. Let τ be the delay time, N be the total number of contacts, N c be the number of coupled (pinned) contacts, and N f = N -N c be the number of detached (sliding) contacts. The fraction of contacts that detach per unit displacement of the sliding block is Γ(v, T ) = dx P (x) Q(x), i.e., when the slider shifts by ∆X, the number of detached contacts changes by N c Γ∆X, so that N f = Γvτ N c . Using N c + N f = N , we obtain N c = N/(1 + Γvτ ) and N f = N Γvτ /(1 + Γvτ ). If we define x = 1/Γ and v d = x/τ , we can write and kBT = 1 (other parameters as in Fig. 3). Inset shows the same in log-log scale.

N c = N 1 + v/v d and N f = N v/v d 1 + v/v d . (48)
The coupled contacts produce the force f k defined above by the steady-state solution of the master equation. The combined dependence which incorporates temperature effects, aging and delay in contact formation, is shown in Fig. 6 for different values of the parameter v d .

However, above we assumed that the sliding contacts experience zero friction, while these contacts may experience a viscous friction force f l = η l v, where η l corresponds to the (bulk) viscosity of the liquid lubricant. In this case the kinetic friction should be additionally multiplied by a factor

β(v) = 1 + v 2 /v η v d , where v η = f k /η l (v η ≫ v d ).
Such a correction may be expected at huge velocities only, e.g., for v ∼ 1 m/s. In this case the function F k (v), after decreasing, reaches a minimum at a velocity v 0 ≈ (v η v d ) 1/2 , and then increases according to a law f k (v) ∝ η l v. Note that the viscous friction which comes from the excitation of phonons in the substrates, as shown in MD simulation [7], may also depend on the velocity, e.g., as η l ∝ v 4 .

VI. MAKING THE LINK WITH EXPERIMENTS

For a real system, the results presented in the previous sections allow the calculation of the kinetic friction force F k (v, T ) provided the parameters of the model are known. In this section section we examine how they can be evaluated from experiments.

The contact parameters k and ω may be estimated with the help of elastic theory [START_REF] Landau | Theory of Elasticity[END_REF]. Let us assume that a contact has a cylinder shape of height h (the thickness of the interface) and radius r c , so that it is characterized by the section A i = πr 2 c , the (geometrical) inertial momentum I = πr 4 c /4, a mass density ρ and a Young modulus E. If the cylinder foot is fixed and a force ∆f is applied to its top, the latter will be shifted on the distance ∆x = ∆f h 3 /3EI (the problem of bending pivot, see Sec. 20, example 3 in Ref. [START_REF] Landau | Theory of Elasticity[END_REF]). Thus, the effective elastic constant of the contact is k = ∆f /∆x = 3EI/h 3 . The minimal frequency of bending vibration of the pivot with one fixed end and one free end, is given by ω ≈ (3.52/h 2 )(EI/ρA i ) 1/2 (see Sec. 25, example 6 in Ref. [START_REF] Landau | Theory of Elasticity[END_REF]).

Next, let a be the average distance between the contacts, so that the total area of the interface is A = N a 2 , and introduce the dimensionless parameter γ c = r c /a (γ c < 0.5). The threshold distance x * may be estimated as follows. At the beginning, when all contacts are in the unstressed state, the maximal force the slider may sustain is equal to F * ≈ N kx * (this force corresponds to the first large stick spike in the F (t) dependence at the beginning of stick-slip motion at low driving). Thus, we obtain that kx * ≈ a 2 σ * , where σ * = F * /A is the maximal shear stress.

Let us consider a contact of two rough surfaces and assume that a = h = r c . Then we obtain

ω ≈ 1.76 E/ρ r c (49)
for the attempt frequency,

k = (3π/4)Er c (50)
for the contact elasticity, and

x * = r 2 c σ * /k (51)
for the threshold distance. For steel substrates we may take ρ = 10 4 kg/m 3 the mass density, E = 2 × 10 11 N/m 2 for the Young modulus, and σ c = 10 9 N/m 2 for the plasticity threshold. Assuming that σ * = σ c and r c ≈ 1 µm, we find that ω ≈ 7.9 × 10 9 s -1 , k ≈ 4.7 × 10 5 N/m, x * ≈ 2.1 × 10 -9 m, b ≈ 2.7 × 10 8 for room temperature (i.e., b ≫ 1), so that the crossover velocity is quite low, v * ≈ 0.03 µm/s. Now let us consider a lubricated system, e.g., the one with a few OMCTS layers as studied by Klein [29] and Bureau [27], and assume that the lubricant consists of solidified islands which melt under stress as proposed by Persson [25]. In this case, instead of using the Young modulus, let us assume that x * = r c ; this allows us to find the parameter EI = ah 3 σ * /3γ c . Then, the elastic constant is k = aσ * /γ c , the attempt frequency is

ω ≈ 1.15 σ * /ahργ 3 c , (52)
the parameter b is given by

b ≈ γ c a 3 σ * /2k B T , (53)
and in the case of b ≫ 1 the crossover velocity is For a four-layer OMCTS film [START_REF] Klein | [END_REF] one may take ρ = 956 kg/m 3 , h ≈ 3.5 × 10 -9 m, F * ≈ 2 × 10 -5 N and A ≈ 10 -10 m 2 so that σ * ≈ 2×10 5 Pa. Assuming γ c = 0.5 and a ≈ 1 µm, we obtain for room temperature, k B T = 4 × 10 -21 J, that ω ≈ 8 × 10 8 s -1 and b ≈ 1.25 × 10 7 , i.e. this system is in the low-temperature limit too, although the crossover velocity is much higher than for rough surfaces, v * ≈ 16 µm/s.

v * ≈ γ c ωa/3C 0 b ≈ k B T / a 5 hρσ * γ 3 c . (54)
Moreover, we may calculate the dependence f k (v) for different thicknesses of the lubricant film. If the film consists of n l layers, then the film thickness is h = n l d, where d ≈ 8.75 Å is the diameter of the OMCTS molecule. Let us assume that the maximal shear stress exponentially decreases with the number of layers according to the results of MD simulation [7], σ * = σ 0 e -βn l , where β ∼ 1 is a numerical constant. Taking σ 0 = 4 × 10 6 N/m 2 and β = 1.5, we obtain the dependences of the shear stress σ = F k /A on the shear rate γ = v/h shown in Fig. 7, which may be compared with the experimental dependences (Fig. 2a) of Bureau [27].

Note that our approach may overestimate the value of the crossover velocity v * . First, the crossover will occur earlier if the delay and/or aging effects play a significant role. Besides, at low temperatures the stiff contacts lead to higher "viscosity" and lower values of v * than the soft contacts considered above (see Fig. 2). Second, we completely ignored the elastic interaction between the contacts. If the latter would be incorporated, a breaking of one contact may stimulate neighboring contacts to break as well, i.e., the value of the parameter a should describe such a cooperative "contact" size which may be much larger than those of individual ones.

Giving a quantitative evaluation of the influence of aging on the velocity dependence of the friction coefficient is harder than for the temperature dependence due to insufficient experimental data. Aging appears to cause a decrease of friction as velocity increases, and thus, when such a behavior is observed experimentally [9,13] it can be considered as a strong indication of the presence of aging. Our analysis indicates that the combined effect of temperature and aging leads to a maximum in the friction coefficient versus velocity. Therefore, when aging is manifested by a decreasing friction versus velocity, extending the experiments to lower velocities and temperatures might detect the maximum and thus provide some quantitative data to evaluate the aging parameters.

Although the aim of our work was to find the dependence of the kinetic friction on the driving velocity, our approach allows us to find the dependence on temperature as well. However, the behavior of a real tribological system is more involved, because all parameters may depend on temperature T in a general case. For example, the delay time τ may exponentially depend on T if the formation of a new contact is an activated process [22]; the same may be true for the aging rate D. In this case one may obtain a nonmonotonic temperature dependence of friction with, e.g., a peak at cryogenic temperatures [21].

VII. CONCLUSION

In this study we determined the dependence of the kinetic friction force in the smooth sliding regime on the driving velocity. In a general case, the friction linearly increases with the velocity (this creep motion may be interpreted as an effective "viscosity" of the confined film), passes through a maximum and then decreases due to delay/aging effects. The decay may be followed by a new growth in friction in the case of liquid lubricant. Estimation showed that for the contact of rough surfaces, the initial growth of friction should occur at quite low velocities, v ≪ 0.1 µm/s, so that for typical velocities the friction is independent on velocity in agreement with the Coulomb law. However, for the case of lubricated friction with a thin lubricant film which solidifies due to compression, the f k (v) dependence is essential, and the linear dependence may stay valid up to velocities v ∼ 10 ÷ 10 3 µm/s. At higher velocities the growth saturates and the f k (v) dependence may be fitted by a logarithmic law. The latter velocity interval is narrow if the distribution of static thresholds is wide; the logarithmic law may be found analytically for a wide interval of velocities when the thresholds are approximately identical, i.e., for the singular distribution P c (x) = δ(x -x s ).

We emphasize that our approach is only valid for a system with many contacts, for example, N > 20 at least [21]. When the contact is due to a single atom as it may occur in the AFM/FFM devices, the friction can be accurately described by the Prandtl-Tomlinson model and should follow the logarithmic f k (v) dependence, f k (v) ∝ (ln v/v 0 ) 2/3 [30]. But if the AFM/FFM tip is not too sharp so that the contact is due to more than one atom, the logarithmic dependence is only ap-proximate and, moreover, for some systems the friction may decrease with the velocity which has to be attributed to the aging/delay effects [8,9,31].

In this work we had in mind that contacts correspond to real asperities in the case of the contact of rough surfaces or to "solid islands" for the lubricated interface. However, the ME approach also operates when the contact is due to long molecules which are attached by their ends to both substrates. Such a system was first studied by Schallamach [22] and then further investigated by Filippov et al. [20], Srinivasan and Walcott [23], and Barel et al. [21]. Note that when all molecules are identical, they are characterized by the same static threshold, i.e., this system is close to the singular one, where the logarithmic f k (v) dependence has to have a wide interval of operation.

Finally, let us discuss restrictions of our approach. First of all, we assumed the somehow idealized case of wearless friction; wearing may mask the predicted dependences. Besides, the interface is heated during sliding; this effect is hard to describe analytically as well as to control experimentally. Then, we did not estimated the delay/aging parameters; moreover, these parameters, e.g., the delay time τ , may depend on the driving velocity v. Besides, we assumed the simplest mechanism of aging described by the Smoluchowsky equation, while the real situation may be more involved, e.g., it may correspond to the Lifshitz-Slözov mechanism [18]. Also, we assumed that the reformed contacts appear in the unstressed state, R(x) = δ(x) in Eq. ( 1), which may not be the case in real systems.

The most important issue, however, is the incorporation of the elastic interaction between the contacts as well as elastic deformation of substrates at sliding. This point certainly deserves a detailed investigation and is the topic of our future work.

FIG. 1 :

 1 FIG. 1: (Color online): The friction force f k as a function of the driving velocity. Dash-dotted blue line shows the low-v approximation (29), dashed red line shows the high-v approximation (31), and dotted magenta line shows a logarithmic fitting. k = 1, ω = 1, x * = 1, kBT = 1.

FIG. 2 :

 2 FIG. 2: (Color online): The friction force f k as a function of the driving velocity v for soft (dotted) and stiff (solid curves) contacts at low temperature kBT = 0.01 (blue) and high temperature kBT = 1 (black). k = 1, ω = 1, x * = 1.

FIG. 3 :

 3 FIG. 3: (Color online):The kinetic friction force f k as a function of the driving velocity v for different values of the aging rate: D = ∞ (red dotted), 1 (blue short-dashed), 10 -1 (black solid), 10 -2 (magenta dash-dotted), 10 -3 (blue shortdotted), and D = 0 (red dotted curve). k = 1, ω = 1, kBT = 1; the initial and final Pc(x) distributions are given by Eq. (9) with x * i = 0.1 and x * f = 1 correspondingly. Dashed curves shows the dependences at T = 0.

∞xFIG. 4 :

 4 FIG. 4: (Color online): (a) The same as in Fig. 3 in loglog scale. (b) The dependences of the effective "viscosity" η = f k /v on the velocity v (dashed lines show power-law fits).

)FIG. 5 :

 5 FIG. 5: (Color online): The kinetic friction force f k as a function of the driving velocity v for different values of the aging rate D = 1 (magenta dashed), 10 -1 (blue solid), 10 -2(black dash-dotted) and 10 -3 (red dotted curve) as compared with approximate expressions (thin solid curves). k = 1, ω = 1, kBT = 1; the initial and final Pc(x) distributions are given by Eq. (9) with x * i = 0.1 and x * f = 1 correspondingly.

FIG. 6 :

 6 FIG. 6: (Color online): The kinetic friction force f k as a function of the driving velocity v for different values of the delay time: v d = ∞ (black solid), 10 (red dot-dashed), 1 (magenta dotted) and 0.1 (blue dashed curve) for D = 0.1 and kBT = 1 (other parameters as in Fig.3). Inset shows the same in log-log scale.

1 )FIG. 7 :

 17 FIG. 7: The shear stress σ as a function of the shear rate γ for different values of the OMCTS film thickness from n l = 2 to 6 monolayers. Inset: the crossover velocity v * as a function of the number of layers.

Acknowledgments

This work was supported by CNRS-Ukraine PICS grant No. 5421.