T Gu

D S Dean

Universal time-dependent dispersion properties for diffusion in a one-dimensional critically tilted potential

We consider the time-dependent dispersion properties of overdamped tracer particles diffusing in a onedimensional periodic potential under the influence of an additional constant tilting force F . The system is studied in the region where the force is close to the critical value F c at which the barriers separating neighboring potential wells disappear. We show that, when F crosses the critical value, the shape of the mean-square displacement (MSD) curves is strongly modified. We identify a diffusive regime at intermediate-time scales with an effective diffusion coefficient which is much larger than the late-time diffusion coefficient for F > F c , whereas for F < F c the late-time and intermediate-time diffusive regimes are indistinguishable. Explicit asymptotic regimes for the MSD curves are identified at all time scales.

I. INTRODUCTION

In a variety of physical systems, the motion of tracer particles can be described by Fokker-Planck equations or their associated stochastic differential or Langevin equations [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF][START_REF]Stochastic Differential Equations[END_REF][START_REF] Gardiner | Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences[END_REF]. In such systems, the motion results from the combined action of deterministic, the so called drift, and stochastic forces whose amplitudes are given by a diffusivity or local diffusion constant. In a periodic heterogeneous medium in which the local transport coefficients are constant in time but spatially periodic, the late-time large scale dynamics can be characterized by a mean velocity and an effective diffusion tensor which characterize the mean drift and the spatial extension of a cloud of initially close tracer particles. These effective transport coefficients can be widely different from typical microscopic ones and are important in a number of phenomena, such as mixing, pollutant spreading, or chemical reactions [START_REF] Le Borgne | [END_REF][5][6][7][8].

The late-time effective transport coefficients have theoretically been characterized in a number of systems. For example, they can be calculated in the case of tracer particles diffusing in incompressible hydrodynamic flows [9][10][11][12][13][14] or in porous media [15][16][17][18][19][20] using homogenization theory. In another context, dispersion properties were derived using methods of statistical physics for particles diffusing in periodic potentials with uniform molecular diffusivity [21][22][23][START_REF] Zwanzig | Proc. Natl. Acad. Sci[END_REF][START_REF] Gennes | [END_REF][26]. The problem of dispersion in one-dimensional (1D) systems also was investigated at length [27][28][29][30][31][32][33][34][35][36][37][38][39]. A remarkable prediction in the case of diffusion in a one-dimensional tilted potential is a huge increase in the effective diffusivity when the tilting force approaches a critical value [27], a phenomenon observed in various experiments [40][41][42][43] which was recently used to estimate the energy barrier opposing the steps of a rotary molecular motor [44]. More recently, analytical results demonstrated how an external force can influence dispersion in periodic systems in higher dimensions [45].

These results, however, characterize only the late-time effective diffusivity, whereas dispersion can also be characterized by the time-dependent mean-square displacement (MSD) which is routinely measured for example in single particle tracking experiments. The general theory developed in Ref. [45] shows that, starting from steady state initial conditions, the average drift is independent of the time whereas the MSD actually evolves in time, eventually attaining the late-time diffusive limit. The approach to the diffusive limit has been calculated in equilibrium systems [46,47]. Apart from approximate forms of the distribution of particles [48,49] in tilted potentials, little is known about the temporal evolution of dispersion coefficients in general nonequilibrium periodic media.

Recently, a very general formalism was proposed to calculate the MSD in a wide class of nonequilibrium systems [45,50]. This formalism examines the full temporal behavior of dispersion, and here we use it to analytically calculate the time-dependent MSD of overdamped particles diffusing in a 1D tilted periodic potential. We derive exact and explicit asymptotic expressions for the dispersion at various time scales in the case where the external force is close to its critical value (at which the barriers between successive potential wells disappear). In earlier studies [27,29], it has been shown that, when the tilting force is close to its critical value, the tracer particles spend most of their time in a very narrow window of positions and thus the effective diffusivity takes a universal form depending only on the properties of the potential near these positions. Here we will see that this property carries over to the time-dependent MSD, which admits simple and universal forms at various time scales.

The outline of the paper is as follows. The model is briefly introduced in Sec. II. In Sec. III, we derive a formula for the MSD in terms of first-passage time (FPT) densities. This formula is analyzed asymptotically in Sec. IV where we show that the shape of the MSD curves show a remarkable change when F crosses its critical value. Our predictions are in excellent agreement with results obtained by simulating the stochastic trajectories.

II. MODEL AND QUANTITIES OF INTEREST

We consider the motion of an overdamped tracer particle of position X(t) at time t in a 1D space, moving in a periodic potential V (x) (of period L with x as the spatial coordinate) and subject to an additional external tilting force F at finite temperature T [see Fig. 1(a)]. The overdamped dynamics 1 FIG. 1. (a) Sketch of the problem investigated in this paper: We calculate the dispersion properties of particles diffusing in periodic tilted potentials. When the noise is weak and the force is close to its critical value, the tracer particles spend most of their time in narrow regions, termed the slow regions, of typical size l * much smaller than period L. (b) Value of the stationary PDF P s (x) inside the slow region obtained by evaluating Eqs. (19) and (20) for ε = -5 (F below the critical force F c ), ε = 0 (F = F c ), and ε = 5 (F > F c ). The inset: same curves on a semilogarithmic scale.

satisfies the force balance equation, which in Langevin form reads

ζ ∂ t X(t) = -V (X(t)) + F + ξ (t), ( 1
)
where ζ is the frictional drag coefficient and the thermal fluctuating forces ξ (t) have zero mean white noise Gaussian statistics characterized by the correlation function ξ (t)ξ (t ) = 2k B T ζ δ(tt ). We denote D 0 = k B T /ζ as the local molecular diffusivity of the tracer particle, and we define the drift field,

u(x) = ζ -1 [-∂ x V (x) + F ]. (2)
Equivalently, the process can be described by the Fokker-Planck equation,

∂ t P = -∂ x [u(x)P ] + D 0 ∂ 2 x P , (3)
where P (x,t) is the probability density function of particles at positions x at time t.

In this paper, we aim to calculate the time-dependent dispersion quantified by the MSD function ψ(t) defined as

ψ(t) = [X(t) -X(0)] 2 -X(t) -X(0) 2 , ( 4
)
where • • • denotes ensemble averaging over realizations of the white noise. Note that we assume here that the system has reached a steady state at time t = 0 (in the sense that the probability distribution over a unit cell of one period is the stationary one).

III. GENERAL EXPRESSION OF THE MEAN-SQUARE DISPLACEMENT

Our starting point is the following Kubo formula, derived in Refs. [45,50]:

ψ(s) = 2D 0 s 2 - 2 s 2 L 0 dx L 0 dx 0 u(x)u * (x 0 ) × P (x,s|x 0 ) - P s (x) s P s (x 0 ). (5)
Here, ψ(s) is the Laplace transform of the MSD, i.e., ψ(s) = ∞ 0 ψ(t)e -st dt, whereas P (x,t|x 0 ) is the propagator of the process mod L (with periodic boundary conditions), i.e., the probability density function for the tracer particle at x (mod L) at t given an initial position x 0 (also mod L), and P is its temporal Laplace transform. Moreover, P s (x) = lim t→∞ P (x,t|x 0 ) is the probability density function of the position (mod L) in the steady state. Note that P s is not the equilibrium-Boltzmann distribution which is only applicable for finite systems with reflecting or confining rather than periodic boundary conditions. Here it describes a nonequilibrium steady state and is characterized by a nonzero flux J s , which has both a convective and a diffusive (Fickian) component,

J s = u(x)P s (x) -D 0 ∂ x P s (x). (6)
Note that J s must be constant in a 1D problem. Finally, in Eq. ( 5), u * (x) represents the drift of the time-reversed stochastic process [50],

u * (x) = u(x) -2J s /P s (x). (7)
The first step of the present analysis consists of expressing the MSD in terms of FPT densities rather than propagators; this will greatly simplify the asymptotic analysis in the next sections. Consider f (x,t|x 0 ) as the probability density of reaching position x (mod L) for the first time at t, starting from the initial position x 0 . The propagators and the FPT densities can be linked by the following, well known, renewal equation [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF]:

P (x,t|x 0 ) = t 0 dt f (x,t |x 0 )P (x,t -t |x). (8)
Physically, the above equation states that, if a particle reaches x at t, it also means that x was reached for the first time at some earlier instant t and that the particle subsequently reached x again in a time tt . The Laplace transform of Eq. ( 8) is as follows:

P (x,s|x 0 ) = f (x,s|x 0 ) P (x,s|x). ( 9
)
We also consider the FPT density averaged over initial conditions,

f s (x,t) ≡ L 0 dx 0 f (x,t|x 0 )P s (x 0 ), (10)
and we remark that averaging (9) over the stationary distribution for x 0 gives

P s (x) s = fs (x,s) P (x,s|x). ( 11
)
Using the above expressions, Eq. ( 5) for the MSD becomes

ψ(s) = 2D 0 s 2 - 2 s 3 L 0 dx L 0 dx 0 u(x)u * (x 0 )P s (x 0 )P s (x) × f (x,s|x 0 ) fs (x,s) -1 . ( 12
)
This expression for the MSD is well adapted for the asymptotic analysis of dispersion presented in the next section.

IV. DISPERSION IN A CRITICAL TILTED POTENTIAL AT DIFFERENT TIME SCALES

A. Regions of fast and slow motion

We now focus on the case where the external force is very close to the critical tilt force F c = max[V (x)] [Fig. 1(a)]. When F F c , the late-time effective diffusivity varies as D 1/3 0 [27,29] and can therefore be much larger than the molecular diffusivity D 0 when the latter is small. Our goal here is to predict the approach to this diffusive limit, which is universal in the sense that it does not depend on the detailed shape of the potential. Following the notations of Refs. [27,29], we introduce the parameter μ defined as

μ = -V (0)/6, ( 13
)
where we chose the origin of coordinates such that V (x) is maximum at x = 0. Note that we assume here that the potential admits a third derivative; more general cases could be treated with our approach but are not considered here for simplicity. We also introduce a parameter ε to measure the distance to the critical force,

ε = F -F c μ 1/3 (k B T ) 2/3 . ( 14
)
The limit considered here is that of weak noise k B T → 0 while keeping the parameter ε constant. In this limit, since D 0 ∝ k B T , the convective terms dominate over the diffusive terms in Eq. ( 3) almost everywhere except in a narrow region of characteristic size l * located around x = 0. In this region, the convective flux is J c = uP s ζ -1 μx 2 P s whereas the diffusive flux is on the order of J d = D 0 ∇P s D 0 P s /l * . The fluxes J c and J d are on the same order of magnitude when x ∼ l * if we set the characteristic size l * to

l * = (k B T /μ) 1/3 . ( 15
)
The characteristic time in the inner region is

τ * = (l * ) 2 D 0 = ζ μ 2/3 (k B T ) 1/3 . ( 16
)
The time τ * diverges in the small temperature limit, and we therefore refer to the region around x = 0 as a slow region as opposed to the outer region called the fast region.

Let us briefly derive the stationary probability density and steady state flux; although these quantities are known in the literature, we show below that they play an important role in the dispersion properties. Since each particle spends most of its time in one of the slow regions, the steady state stationary density P s (x) of positions (mod L) is localized in these regions where it satisfies

-ζ -1 ∂ x [(3μx 2 + F -F c )P s ] + D 0 ∂ 2 x P s = 0, (17)
where we have used a Taylor expansion of the potential V (x) at next to leading order around x = 0. We introduce the dimensionless position y = x/ l * , probability density Ps (y) = P s l * , and flux Js = τ * J s . Using these notations, integrating Eq. ( 17) and comparing with Eq. ( 6) leads to

-[(3y 2 + ε) Ps ] + ∂ y Ps = -Js . ( 18
)
Assuming that Ps (y) vanishes for y → ∞ (that is, for x l * ), the solution of this equation is

Ps (y) = Js ∞ y du e -u 3 -εu+y 3 +εy . ( 19
)
The normalization imposes that 1

= L/2 -L/2 dx P s (x)
∞ -∞ dy Ps (y) and thus leads to the identification of the dimensionless current,

J -1 s = ∞ -∞ dy ∞ 0 du e -(u+y) 3 -εu+y 3 . ( 20
)
The rescaled stationary PDF in the slow region is represented in Fig. 1(b) where one observes the transition from narrow distributions shifted in the region x < 0 when ε < 0 (that is, for forces F < F c ) to broader and more centered distributions when F > F c . The dimensionless flux is represented in Fig. 2(a), it almost vanishes for F < F c , and then significantly increases when F > F c . In the outer region |x| l * , the convective term takes over the diffusive term in Eq. ( 6), and we simply obtain

P s (x) J s /u(x), (21)
physically this means that, in the regions where drift dominates, the probability of presence is inversely proportional to the speed. We note that the expressions in the outer and the inner regions [Eqs. (19) and ( 21)] can be matched in the region l * |x| L where both approximations give P s (x) J s ζ /(3μx 2 + F -F c ).

B. MSD at intermediate-time scales: ballistic and diffusive regimes

We now investigate the properties of the MSD ψ(t) at different time scales. First, at very short times, t → 0 (or s → ∞), we have ψ(t) = 2D 0 t where the effective diffusivity is exactly equal to the molecular diffusivity. We do not discuss this regime any further and consider now longer times with the additional condition that t τ * (or equivalently sτ * 1). At this time scale, one can consider that the events of crossing the slow region or escaping from it take an infinite time. In this regime, when x,x 0 are in the fast region with 0 < x 0 < x < L, we can approximate the FPT to x from x 0 by the duration of the deterministic trajectory that links x 0 and x, f (x,t|x 0 ) δ t -

x x 0 dy u(y) . ( 22
)
In the opposite direction where x 0 > x, reaching x from x 0 requires crossing the slow region, which takes a time infinite compared to the time scale considered here. Hence f (x,t|x 0 ) is simply approximated by 0 in this case. Taking all this into account, the FPT density in Laplace space reads f (x,s|x 0 ) exp -s

x x 0 dy u(y) , if 0 < x 0 < x < L, 0, if 0 < x < x 0 < L. ( 23
)
We now calculate fs (s,x), the average FPT density to x, with stationary initial conditions, fs (s,x)

x x - dx 0 P s (x 0 ) f (s,x|x 0 ). ( 24
)
Here the upper integration limit has been set to x because if x 0 > x the FPT from x 0 to x is infinite at this time scale. In turn, the lower integration limit has been set to x -, which is a position such that l * x - L (its precise value will not change the result, see below). The reason for this is that, if x 0 ∼ l * is in the slow region, its escape time is ∼ τ * and therefore almost no trajectory can bring it to x at the considered time scale. For x 0 l * , we can approximate P s (x 0 ) J s /u(x 0 ) [see Eq. ( 21)] and therefore, fs (s,x)

x x - dx 0 J s u(x 0 ) exp -s x x 0 dy u(y) . ( 25
)
The integration over x 0 is now straightforward, fs (s,x) J s s , (26) and the result does not depend on the lower bound of the integral x -as soon as x - L. Furthermore, when x,x 0 are in the fast regions, the approximation u(x)P s (x) = J s = -u * (x 0 )P s (x 0 ) holds [cf. Eq. ( 21)]. Hence, inserting (23) and (26) into Eq. ( 12) and keeping only the dominant contribution to the MSD (coming from the x,x 0 in the fast region), the expression of ψ near criticality is considerably simplified,

ψ(s) 2J 2 s s 3 L 0 dx x 0 dx 0 s exp -s x x 0 dy u(y) J s -1 . ( 27
)
Here we investigate the regime t τ * or equivalently sτ *

1. Since J s ∼ 1/τ * , we obtain s/J s 1 and

ψ(s) 2J s s 2 L 0 dx x 0 dx 0 exp -s x x 0 dy u(y) . ( 28
)
The above equation describes the MSD at time scales smaller than τ * . Let us consider its asymptotics. First, for s → ∞, the terms x x 0 essentially contribute in the integral, and therefore

ψ(s) 2J s s 2 L 0 dx x 0 dx 0 exp -s(x -x 0 ) u(x) . ( 29
)
Performing the integral, we obtain

ψ(s) 2J s s 3 L 0 dx u(x) = 2J s F c L ζ s 3 . ( 30
)
If we invert the Laplace transform, we obtain

ψ(t) = J s F c L ζ t 2 , ( 31
)
and thus the regime identified here is a ballistic regime with a MSD ψ(t) αt 2 . The nontrivial coefficient α identified in Eq. ( 31) is the product of two velocities; the first velocity is the average velocity J s L, whereras the second velocity is the characteristic velocity in the fast region F c /ζ . Note the nontrivial temperature dependence [as (k B T ) 1/3 ] of the coefficient of the MSD in this ballistic regime. Now, let us take the small s limit of Eq. ( 28), for which

ψ(s) = 2J s s 2 L 0 dx x L 2 J s s 2 . ( 32
)
Inverting the Laplace transform leads to ψ(t) = 2D int t with D int as an effective diffusion coefficient at intermediate times,

D int = L 2 J s /2. ( 33
)
Thus, the analysis reveals the existence of a diffusive regime at intermediate-time scales. The effective diffusivity D int is proportional to the flux of particles and therefore increases significantly when the force becomes larger than the critical force [that is, for increasing ε, see Fig. 2(a)]. When t τ * , the motion becomes effectively diffusive with an effective diffusion coefficient D ∞ = (L 2 /τ * )G(ε), where G is the same dimensionless function identified in Refs. [27,29]. It is represented in Fig. 2(b) and shows a maximum for ε 0, leading to the giant enhancement of diffusivity at the critical force. For completeness we provide in Appendix A a derivation of D ∞ within our formalism, which is an alternative to the approach of Refs. [27,29].

The ratio of the late-time over the intermediate-time diffusivity D ∞ /D int is shown in Fig. 2(c). It is almost equal to 1 for forces below the critical force (ε < 0) but then vanishes for larger values of ε. This means that, when the critical force is reached, the shape of the MSD curves changes drastically. For F < F c , one observes a direct transition between a ballistic and the long-time ballistic regime. When F > F c , one observes the intermediate diffusive regime with a diffusivity larger than the effective diffusivity, which translates by an overshoot of the MSD and an apparent regime of subdiffusion.

If we summarize all the results, we obtain

ψ(t) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 2D 0 t, if t t 2 , J s vLt 2 if t 2 t L/v, 2D int t, if L/v t τ * , 2D ∞ t, if τ * t, (34)
with v = F c /ζ as the velocity in the fast region and t 2 = k B T /(F c LJ s ). In order to check these predictions, we performed stochastic simulations of the Langevin equation (1) in the case of a sine potential (see Appendix B for details on the simulation algorithm). The results presented in Fig. 3 confirm the validity of the asymptotic regimes described by Eq. ( 34)

FIG. 3. MSD curves for the diffusion in a sine potential V (x) = V 0 sin(2πx/L) for (a) a tilting force below the critical force (ε = -2), (b) at the critical force (ε = 0), and (c) above the critical force (ε = 2). Different values of the noise amplitude are indicated in the legend. Symbols: results of the simulations of the Langevin equation (1) (see Appendix B for details). Lines: results of the numerical evaluation of the Kubo formula (5). Bold lines: asymptotic results of Eq. ( 34) (with no free parameters). Note that the regime ψ = 2D 0 t is not represented but is present at short times. Units of length and times are L and ζ L 2 /V 0 , respectively.

for tilting forces that are below, above, or equal to the critical force.

V. CONCLUSION

In this paper we have studied the time-dependent dispersion properties of particles diffusing in a near critically tilted one-dimensional periodic potential. We have derived explicit asymptotic expressions for the MSD at different time scales [Eq. (34)]. The approach to the late-time diffusive limit depends only on a small number of parameters which characterize the potential in particular regions where the dynamics is slow. The approach to the diffusive limit is therefore universal in the sense that it does not depend on the details of the potential shape.

It has been proposed that the giant increase in the late-time diffusivity that occurs at the critical force [27] can be used to estimate barriers of potential energy [40,44]. Here we have shown that the time-dependent dispersion properties are strongly modified when the force becomes larger than F c , in particular, we have found the presence of a second diffusive regime at intermediate-time scales. The strong difference between these two diffusion coefficients could be used as another signature of the effect of the crossing of the critical force.

In this paper, we have quantified the transition between the short-time regime of molecular diffusion and the late-time effective diffusion. The mean-square displacement between these regimes can be regarded as "anomalous" in the sense that it is a nonlinear function of time. Anomalous diffusion can have a variety of origins in different physical systems and may occur in the late-time regime for example for particles undergoing jumps whose sizes or durations follow large distributions [51] when the distribution of energy barriers leads to diverging mean occupation times in the local energy minima or when the convective velocity field has long range correlations [52]. Anomalous diffusion also arises in fractal media [START_REF] Ben-Avraham | Diffusion and Reactions in Fractals and Disordered Systems[END_REF] when considering time-dependent microscopic diffusivities [START_REF] Jeon | [END_REF] or when the tracer trajectory results from a collective dynamics, such as in polymer systems [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Panja | [END_REF] or complex fluids [57,58]. This paper where we study out-of-equilibrium tracer particles in periodic media is an example where one can entirely characterize the anomalous time-dependent dispersion properties at all intermediate-time scales between the molecular diffusion regime and the final late-time regime, which is however one of normal diffusion. mean and variance 1, as the time step, and x i as the position at time t = i . Each trajectory was simulated during a time t m , the position of the particle at the end of a trajectory being used as the initial position for the next trajectory. The MSD was then computed with Eq. ( 4) by averaging over distinct trajectories. For each parameter, was chosen small enough so that we could observe the regime ψ(t) 2D 0 t for small times (not shown in Fig. 3) and we carefully controlled that MSD curves, obtained with different 's, overlapped. In order to ensure that different runs are independent, the time t m was chosen large enough to be located in the late-time diffusive regime. The number of runs was always larger than 8000 for each parameter set. In the late-time regime, where X(t) -X(0) is Gaussian distributed, the 95% confidence intervals can be evaluated to (1 ± √ 8/N run )ψ so that the precision on ψ is approximately ±3%.

FIG. 2 .

 2 FIG. 2. (a) Value of the dimensionless flux Js = J s τ * for near critical forces, which is also proportional to the effective diffusion coefficient in the intermediate-time regime [see Eq. (33)]. (b) Rescaled late-time effective diffusivity D ∞ , deduced from Eq. (A9). (c) Ratio of the late-time effective diffusivity over the intermediatetime effective diffusivity.

APPENDIX A: LATE-TIME EFFECTIVE DIFFUSIVITY

In this appendix we briefly derive an expression for the late-time effective diffusivity D ∞ . If we expand the temporal Laplace transform of the FPT densities for s → 0 we obtain f (x,s|x 0 ) 1τ

where τ (x|x 0 ) = ∞ 0 dt f (x,t|x 0 )t is the mean first-passage time (MFPT) to reach the position x mod L starting from x 0 , whereas τ s (x) is the MFPT to x with x 0 averaged over stationary initial conditions. Using Eq. ( 12), we then find that ψ 2D ∞ /s 2 where the effective diffusivity D ∞ reads

For forces close to the critical force, this expression can be simplified considerably. The dominant contribution to this integral comes from the values of x,x 0 that are outside the slow region in which the diffusive component of the flux is negligible, and thus u(x)P s (x) = J s = -u * (x 0 )P s (x 0 ) [cf.

Eq. ( 21)]. Furthermore, taking x and x 0 in [0,L], we realize that τ (x|x 0 ) is negligible when x 0 < x (because the convection brings the particle to x almost immediately), whereas for x < x 0 this time is equal to the time to cross the slow region, which is exactly the inverse of the flux J s .

We define τ e as the time to escape the slow region starting from stationary initial conditions, and we note that τ s (x) τ e is independent of x. Following these considerations, and neglecting the term D 0 in Eq. (A3), we obtain

Now let us consider τ ∞ (x 0 ), the mean time to escape the slow region, starting from x 0 from which τ e = τ ∞ (x 0 ) can be deduced. We pose y 0 = x 0 /l * and τ ∞ (x 0 ) = τ * τ∞ . Then, the dimensionless MFPT satisfies the backward equation [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF][START_REF] Gardiner | Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences[END_REF],

Noting that ∂ y 0 [e y 3 0 +εy 0 ∂ y 0 τ∞ ] = -e y 3 0 +εy 0 , (A6) the differential equation (A5) can be integrated twice, leading to

where we took into account the condition that τ∞ (y 0 ) vanishes at infinity. Now, the average escape time from the slow region can be calculated by averaging over P s given in Eq. ( 19), leading to

Taking y 0 → -∞ in (A7) we can check that the average time to cross the slow region is also 1/J s . Therefore, we obtain at the end,

where the dimensionless function G(ε) is deduced from Eqs. ( 20) and (A8), is represented in Fig. 2, and is in excellent agreement with the results of Refs. [27,29].

APPENDIX B: DETAILS ON SIMULATIONS

We performed numerical simulations of the Langevin equation (1) by using the algorithm x i+1 = x i + u(x i ) + √ 2 D 0 w i with w i as a Gaussian random variable of zero


Amalya Minasyan

Clément Trovato

Jérôme Degert

Eric Freysz

Etienne Brasselet
email: etienne.brasselet@ubordeaux.fr

Emmanuel Abraham
email: emmanuel.abraham@ubordeaux.fr

Geometric phase shaping of terahertz vortex beams

Geometric phase shaping of terahertz vortex beams

We propose a topological beam-shaping strategy of terahertz (THz) beams using geometric phase elements made of space-variant birefringent slabs. Quasi-monochromatic THz vortex beams are produced and characterized both in amplitude and phase from the reconstructed realtime two-dimensional imaging of the electric field. Nonseparable superpositions of such vortex beams are also obtained and characterized by two-dimensional polarimetric analysis. These results emphasize the versatility of the spin-orbit electromagnetic toolbox to prepare on-demand structured light endowed with polarization-controlled orbital angular momentum content in the THz domain, which should find many uses in future THz technologies. The unique properties of terahertz (THz) radiation, such as good penetration through various nonconductive materials, low scattering, free space propagation, low photon energy, and broad spectral bandwidth, encouraged the development of THz photonics [1]. More precisely, the initial developments of THz optics mainly consisted of engineer sources and detec tors, which subsequently implied the need for THz optical el ements such as lenses, splitters, polarizers, and retarders [2]. Nowadays, there is still plenty of room for technological im provements, for instance, regarding THz beam shaping. A par ticular situation corresponds to topological beam shaping that refers to tailor made structured fields endowed with phase and/ or polarization singularities. Such beams are now common in the visible domain and have already found a huge number of applications (e.g., in sensing, microscopy and astronomical im aging, trapping and manipulating of matter, information and communication technologies), whereas fewer works have been reported so far in the THz domain.

A basic ingredient of structured light corresponds to phase singularities (or "vortices") that are location in space, where the phase is undefined. Such singularities are locally described in a plane by a pure phase factor of the form expilφ, with l as an integer and φ as the usual azimuthal angle. Many techniques have been implemented to engineer beams carrying on axis phase singularities, usually called vortex beams, in the visible domain and, quite naturally, a straightforward strategy con sisted of performing mapping to the THz domain. One can mention the use of refractive three dimensional spiral phase plates, with an early demonstration in the sub THz domain 20 years ago [3], and more recent realizations covering the THz domain, for instance, at 0.3 THz [4], 0.62 THz [5], and in the range from 2 to 4 THz [6]. Flat THz vortex beam generators have also been realized, for instance, at 0.75 THz using an array of wavelength size V shaped antennas [7], at 0.31 THz [8], in the range from 1 to 1.6 THz [9], and at 2.1 THz [10] using a computer generated holograms.

Vector beams, which refer to light beams with space variant polarization states in the transverse plane, represent another kind of structured beam that has attracted interest in THz op tics. So far, the most studied cases are radially and azimuthally polarized beams. Radially polarized beams have been produced by using different kinds of radial polarizers, such as a tapered waveguide at 0.1 THz [11] and azimuthally patterned sub wavelength dichroic gratings in the range from 0.3 to 0.9 THz [12]. Concentric and radial antennas have been also proposed to generate radial and azimuthal polarization states, respectively, in the range from 0.5 to 2 THz [13]. On the other hand, various nonlinear processes have also been explored to generate THz vector beams [14 16]. Recently, the broadband generation of radially polarized THz beams has been shown in the range from 0.75 to 2 THz [17] by optical rectification of radially polarized near infrared pulses. Moreover, broadband THz vortex beams of topological charge l 1 have been extracted by selective polarization filtering of a generated poly chromatic THz vector beam [17].

Although various aspects of singular THz optics have been already explored, there is not yet a well established topological shaping approach enabling simple yet efficient generation of THz beams with arbitrary orbital angular momentum states. Here, we propose to exploit the geometric phase (see Ref. [18] for an overview in the context of spin orbit interactions of light) in order to achieve polarization controlled THz vortex beams and superposition them from "regular" THz beams free from any singularities. This is done by using a single element, namely, a space variant birefringent slab endowed with an azi muthally varying in plane optical axis orientation. Indeed, this allows versatile interplay between the spin (polarization degree of freedom) and orbital (spatial degree of freedom) electromag netic angular momentum. By doing so, this work extends to the THz domain's previous use of geometric phase optical elements [19] in the visible domain [20] and should find many uses in future THz technologies.

The experimental setup is presented in Fig. 1. An amplified Ti:Sapphire femtosecond laser source centered at an 800 nm wavelength with a 1 kHz repetition rate, 1.5 mJ pulse energy, and 50 fs pulse duration is split into pump and probe beams for THz field generation and detection. Regular broadband linearly polarized THz pulses over the typical range from 0.1 to 2.5 THz are prepared by optical rectification of amplified laser pulses in a (110) cut ZnTe nonlinear crystal [21], and a quasi monochromatic beam is obtained by placing a 1 THz bandpass filter (200 GHz full width at half maximum bandwidth) after the ZnTe crystal. Then, topological shaping of THz pulses ba sically relies on the use of a half wave quartz birefringent retarder whose optical axis orientation angle ψ in the transverse plane has the form ψ qφ (q half integer), as introduced in the framework of optical vortex generation in the visible do main [22]. The effect of such an optical element (often called a q plate) on an incident field can be easily grasped from the Jones calculation by applying the φ dependent Jones matrix J to an incident, fully polarized monochromatic plane wave de scribed by the Jones vector E in . Indeed, neglecting diffraction, this gives an output field E out J • E in , which is expressed on the x; y basis as

E out cos 2ψ sin 2ψ sin 2ψ cos 2ψ • a b ; (1)
where the complex number a; b defines the incident polarization state. In this work, we consider two situations: (i) incident circular polarization, namely a; b 1∕ 2 p 1; iσ, where σ 1 refers to left/right polarization handedness and (ii) incident linear polarization, namely, a; b cos θ; sin θ where θ refers to polarization state azi muth. Without lack of generality, we restrict our study to the case q 1∕2, with the subsequent generalization to arbitrary values of q being straightforward. In addition, we practically use a 8 segmented half wave retarder for a 1 THz frequency made of quartz crystal (from Tydex; see Fig. 1), instead of a smoothly varying azimuthal profile for the angle ψ, which does not affect the results. The case (i) is carried out by preparing an incident circularly polarized THz beam with helicity σ, which is done by inserting a 1 THz quarter wave retarder after the bandpass filter [see in set (a) of Fig. 1]. One gets from Eq. ( 1)

E out 1 2 p 1 iσ e iσφ ; (2)
namely, a circularly polarized vortex beam with a helicity con trolled topological charge l σ. The real time generated THz vortex electric field is fully reconstructed in two dimensions using a pump probe detection scheme. This is made using a high speed 1 kHz CMOS camera coupled to an electro optical detection owing to a second (110) cut ZnTe crystal, where the THz beam and delayed probe laser pulse interact [23,24] via a second order nonlinear Pockels effect (see Fig. 1).

The resulting image of the modulated probe light intensity is then read out on the camera, from which the time dependent THz electric field amplitude is deduced. Indeed, through the electro optic effect, the probe pulse experiences a polarization modulation that is proportional to the strength of the THz electric field. This polarization modulation is measured by sending the probe beam through an analyzer; then, an objective lens images the crystal plane onto the CMOS camera [25].

Note that the signal is obtained by setting the analyzed THz field polarization state to linear, which is made by placing another quarter wave retarder after the geometric phase element, which does not affect the sought after spiraling phase spatial distribution of Eq. ( 2). Images are acquired at a 500 Hz frame rate via dynamic subtraction of the background [26].

The signal to noise ratio is eventually improved by averaging 500 successive frames, which corresponds to a 1 s integration time. The temporal evolution of the THz electric field pulse is then reconstructed by controlling the time delay between the THz and probe pulses via an optical delay line (see Fig. 1). The results are presented in Fig. 2 at 1 THz. Indeed, we can directly use the broadband THz pulses (0.1 2.5 THz), since the recorded data contain both the amplitude and phase information at all frequencies. This implies that transforming the temporal information into spec tral information can be numerically done by a fast Fourier transform. By doing so, the amplitude and phase spatial distri bution data are obtained for each frequency. A demonstration is illustrated in Fig. 3, where the intensity and phase transverse patterns at 1 THz extracted from the polychromatic data are shown in panels (a) and (b), respectively, for σ 1. Again, one thus unambiguously concludes that the generation of a spin controlled unit charge THz vortex, whose linear azimuthal dependence of the phase along the circle line in Fig. 3(a) is quantitatively shown in Fig. 3(c), is in good agree ment with the theory, Φ σφ. In addition, we note that a broadband THz vortex could be easily obtained by mere use of achromatic quarter wave retarders instead of monochromatic ones; however, this would be at the expense of the overall vortex generation efficiency. Indeed, the frequency dependent Jones matrix Jν is expressed as

J cos Δ 2 i sin Δ 2 cos 2ψ i sin 2ψ sin Δ 2 i sin 2ψ sin Δ 2 cos Δ 2 i sin Δ 2 cos 2ψ ; (3)
where Δν πν∕ν 0 is the frequency dependent birefringent phase retardation, ν 0 being the frequency associated with half wave retardance, here with ν 0 1 THz, which assumes frequency independent birefringence, which is a fair approxi mation. One thus gets

E out 1 2 p cos Δ 2 1 iσ i sin Δ 2 1 iσ e iσφ ; (4)
which shows that the σ polarized output component carries a phase singularity with topological charge l σ whatever the frequency ν, up to a frequency dependent intensity modulation factor sin 2 πν∕2ν 0 . Then, any quantity of interest when dealing with vortex generation from polychromatic fields can be readily derived from the latter expression.

Then we explore the case (ii), which is done by preparing an incident linearly polarized THz beam with an azimuth angle θ. One gets from Eq. ( 1)

E out cosφ θ sinφ θ ; (5)
namely, a vector beam characterized by an inhomogeneous lin ear polarization state having an azimuth angle Ψ φ θ. In particular, a radially polarized beam is expected when θ 0, whereas an azimuthally polarized beam is expected when θ π∕2. In practice, THz vector beam generation is demon strated from the spatially resolved azimuth angle Ψ of the electric field. From the definition Ψ arctanE y ∕E x , it is therefore enough to place a linear polarizer after the segmented half wave retarder, as shown in inset (b) of We note that such a direct determination of the electric field orientation is an option that one cannot afford in the visible region of the near infrared domain, where the electric field vec tor can hardly be measured. In both cases, we observe a total increase of 2π of the polarization azimuth angle ψ. Moreover, this angle is equal to 0 if the polar angle φ is 0 or π∕2, respec tively, for radial and azimuthal beams. Finally, the angular dependence of polarization azimuths for radial and azimuthal beams along the circle lines in Figs. Interestingly, the approach by Imai et al. that consists of extracting a vortex beam from a vector beam by projection onto a circular polarization state [17] can be grasped by rewriting Eq. ( 5) as

E out 1 2 1 i e -iφ 1 i e iφ ; ( 6
)
where the ± sign refers to θ 0; π∕2. Indeed, this corre sponds to a nonseparable spin orbit state (i.e., neither the polarization nor the spatial degree of freedom can be factorized) that consists of the superposition of contra circularly polarized vortex fields with opposite unit topological charges. circular polarization filtering of this electric field will thus give a vortex beam with a topological charge l σ. In our case, where the birefringent phase retardation is Δ π, we have access either to spin controlled THz vortex beams or azimuth controlled THz vector beams without need for post polarization filtering.

The broadband features merely depend on the use of achro matic retarders, which is not a serious technological issue. More over, the present demonstration in the particular case where q 1∕2 can be straightforwardly generalized to any half integer value of the parameter q, which would lead to the generation of THz vortex beams with topological charge l 2σq in the case (i) and THz vector beams of order l 2q in the case (ii), similar to what happens in the well documented visible domain.

In conclusion, THz vortex beams and vector beams were generated using geometric phase elements. A high speed CMOS camera combined with EO detection was used to char acterize the beams by employing a simple and direct measure ment of the THz electric field. The possibility of generating broadband vortex beams with adapted optical elements and their characterization by means of a fast Fourier transform was also shown. Finally, we note recent theoretical developments toward alternative nonlinear strategies to generate THz twisted beams [27,28], calling for experimental demonstrations. This is likely to complement the present linear spin orbit THz toolbox to prepare on demand structured light endowed with polariza tion controlled orbital angular momentum content.

Funding. Agence Nationale de la Recherche (ANR) (ANR 10 IDEX 03 02).

  OCIS codes: (110.6795) Terahertz imaging; (080.4865) Optical vortices; (260.5430) Polarization; (060.5060) Phase modulation.

for σ 1 ,

 1 where THz intensity transverse patterns are obtained by averaging the THz pulse intensity for different time delays [Figs. 2(a) and 2(b)]. Instantaneous amplitude transverse patterns are shown in Figs. 2(c) and 2(d) over half a THz period duration (i.e., T ≃ 1 ps, that is ≃300 μm free space delay). As expected, the doughnut shaped intensity patterns shown in Figs.2(a) and 2(b) are reminiscent of the on axis phase singularity that is associated with a value of zero of the THz intensity. On the other hand, the σ dependent spatiotemporal spiraling structure of the electric field shown in Figs.2(c) and 2(d) corresponds to the presence of a phase singularity with a topological charge l σ, since a σ dependent 2π rotation of the pattern per wavelength is observed. This can also be observed in the media file Visualization 1, in which the left part corresponds to σ 1 [i.e., Fig.2(d)] and the right part to σ 1 [i.e., Fig.2(c)].It is noteworthy that the above spin controlled THz vortex generation does not require us to prepare incident THz pulses

Fig. 1 .

 1 Fig. 1. Experimental setup for the generation and detection of THz vortex beams [case (i), see inset (a)] and THz vector beams [case (ii), see inset (b)]. See text for details.

Fig. 1 ,

 1 and to image the vertical and horizontal output beam electric field compo nents, E y and E x . For the individual measurements of E x and E y , the polarization of the laser probe beam was kept horizon tal, and, for a given position of the THz polarizer, the ZnTe detection crystal was rotated in order to maximize the electro optical signal. The corresponding data are shown in Figs. 4(a) and 4(b) for θ 0 and Figs. 5(a) and 5(b) for θ π∕2, from which the Ψ map is evaluated, as shown in Figs. 4(c) and 5(c).

  4(c) and 5(c) are plotted in Figs. 4(d) and 5(d). They are in good agreement with theoreti cal linear fits. Distributions of polarization vectors are recon structed from the polarization azimuth angle values for radial and azimuthal beams, as shown in the insets of Figs. 4(d) and 5(d), respectively.

Fig. 2 .

 2 Fig. 2. (a) and (b) Intensity patterns of the 1 THz vortex beam ob tained in case (i) for σ 1. (c) and (d) Temporal evolution of the THz electric field transverse distribution revealing phase singularity of topological charge l σ. Scale bar: 10 mm (see Visualization 1).

Fig. 3 .

 3 Fig. 3. Intensity (a) and phase (b) distributions of the generated vortex beam at 1 THz extracted from the polychromatic field proc essed by the geometric phase vortex generator where σ 1; see text for details. (c) Azimuthal angle dependence of the phase along the circle line shown in panel (a). Markers: experimental data. Solid line: theory.

Fig. 5 .

 5 Fig. 5. Analysis of azimuthally polarized THz beam (θ π∕2). Panel descriptions are similar to those of Fig. 4.

Fig. 4 .

 4 Fig. 4. Analysis of radially polarized THz beam (θ 0). Transverse distributions of the vertical (E y ) and horizontal (E x ) electric field com ponents are shown in panels (a) and (b), respectively. (c) Spatial dis tribution of the azimuth angle Ψ arctanE y ∕E x . (d) Dependence of the polarization azimuth Ψ of the radial vector beam: circles, experimental data along the circle line shown in panel (c). Markers: experimental data; Solid line: theory.

Acknowledgment. This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of "the Invest ments for the future" Programme IdEx Bordeaux LAPHIA.




  12 normal French speakers (6 women & 6 men)  4 consonants with voicing & manner contrast  7 phonations types : spontaneous (baseline), low or high, soft or loud & breathy or pressed  Control of intra-speaker difference to baseline for pitch (F0), intensity (SLP) & quality (H1-H2)  24 trains of 8 syllables, like "pafavabapafavaba"  For each, selection of the 4 central consonants  Measure of the relative fundamental frequency on the 10 cycles of the vowel offset & onset Abstract #2575950 Camille F Robieux & Christine Meunier Aix-Marseille Univ, CNRS, LPL, Aix-en-Provence, France (camille.robieux@lpl-aix.fr) Relative fundamental frequency of voiced and voiceless stops and fricatives for different phonation types in normal French speakers This study was conducted thanks to the financial support of the region Provence-Alpes-Côte d'Azur and the compagny SQ-Lab Low-soft-breathy phonation Δ F0 (ST) Δ SPL (dB) Δ H1-H2 (dB)



  Measure of 10 periods (time between 2 pulses) at vowel offset & onset  Calcul of the instantaneous fundamental frequencies (F0) = 1/period  Conversion of all frequencies from Hertz (Hz) to Semitones (ST) Calcul of the relative fundamental frequency (RFF) to the farthest point from the consonant (cycle 1 in vowel offset & cycle 10 in vowel onset)  RFF varies with vocal effort: higher RFF/lower effort in LSB phonation (+ voiceless consonants) & lower RFF/higher effort in HLP phonation (+ voiced consonants)  Possible to decrease but not to increase the vocal effort for voiced consonants & possible to increase but not to decrease the vocal effort for voiceless consonants  In spontaneous phonation: production of voiced consonants at their maximum vocal effort but production of voiceless consonants at their minimum vocal RFF lower for voiced than for voiceless stops & fricatives at vowel onset  RFF for voiced consonants: (except /b/ at onset: LSB > Baseline > HLP) Low-soft-breathy phonation > Baseline = High-loud-pressed phonation  RFF for voiceless consonant: (except /p/ at offset: LSB > Baseline = HLP) Low-soft-breathy phonation = Baseline > High-loud-pressed phonation LSB phonation Baseline HLP phonation