text:
# The skewed multifractal random walk with applications to option smiles

## 1 Introduction

Volatility clustering in financial markets is a well known phenomenon. More surprising is the fact that volatility correlations are found to be long-ranged, and cannot be characterized by a single correlation time. Recent empirical works  suggest that the volatility correlation function of various assets (in particular stocks) actually decays as a power-law of the time lag, with a rather small exponent. Furthermore, the volatility fluctuations are found to be close to (but not exactly) log-normal . These observations have lead several authors  to propose an interesting class of multifractal models, where the log-volatility is a Gaussian random variable with a correlation function that decays in time as a logarithm. An explicit continuous time construction of such a process was proposed and studied in , and was showed to exhibit strict multifractal properties, in the sense that the even moments of the log-price difference scale with a non trivial power of the time lag (more precise statements will be given below). This means in particular that the kurtosis of the process decreases only very slowly with the time lag, in contrast with most simple models of stochastic volatility, where the kurtosis drops exponentially with time. This model is therefore of interest for option pricing, because it is consistent with smiles that flatten only very slowly with time . The ‘multifractal random walk’ (mrw) constructed in  is however explicitly symmetric, in the sense that all the odd moments of the log-price difference are strictly zero. The same is true of the ‘cascade’ construction of Mandelbrot, Fisher and Calvet  (see also  for an interesting causal extension of this model). On the other hand, an important stylized fact of the time series of stocks and stock indices is the so called ‘leverage effect’: past price returns and future volatilities are negatively correlated . This effect was documented quantitatively in . This effect is particularly strong for stock indices, and induces a significant (negative) skewness in the distribution of price returns. The aim of this paper is to study a generalization of the mrw which accounts for the leverage effect and the corresponding skewness. A particularly interesting class of models preserves the multiscaling property of the mrw and extends it to odd moments. We discuss several financial applications of these models, in particular to option pricing. In the presence of non zero skewness, the volatility smile itself becomes skewed. Approximate theories, based on a cumulant expansion, are also discussed in this context.

## 2 The Multifractal Random Walk (mrw)

Let $`(X_t,t0)`$ be a stochastic process, with stationary increments. We note by $`\delta _{\mathrm{}}X_t`$ the increments $`X_{t+\mathrm{}}X_t`$ and by $`M(q,\mathrm{})`$ the $`qth`$ moment $`E(|\delta _{\mathrm{}}X_t|^q)`$. $`(X_t,t0)`$ is said to be a fractal process if

$$M(q,\mathrm{})=C_q\mathrm{}^{\zeta _q}.$$
(1)
When the function $`\zeta _q`$ is linear in $`q`$, one speaks of monofractal process. This is the case of the brownian motion for which $`\zeta _q=q/2`$, and more generally of self-similar processes , for which the following equality in law holds

$$\delta _b\mathrm{}X_t=b^H\delta _{\mathrm{}}X_t,$$

where $`b`$ is an arbitrary scaling factor. When the function $`\zeta _q`$ is non linear in $`q`$, one speaks of multifractal process. In this case, as emphasized in , Eq. (1) can in fact only hold in a certain ‘scaling regime’, $`\mathrm{}T`$. The construction of  is based on the following discretized process:

$$X_{\mathrm{\Delta }t}(t)=\underset{k=1}{\overset{\frac{t}{\mathrm{\Delta }t}}{}}\delta X_{\mathrm{\Delta }t}[k];\delta X_{\mathrm{\Delta }t}[k]=ϵ_{\mathrm{\Delta }t}[k]e^{\omega _{\mathrm{\Delta }t}[k]}.$$

For financial applications, $`X_{\mathrm{\Delta }t}(t)`$ can be thought of as the logarithm of the price at time $`t`$. The quantities $`ϵ_{\mathrm{\Delta }t}[k]`$ and $`\omega _{\mathrm{\Delta }t}[k]`$ are independent Gaussian variables with the following covariance structure:

$`E(ϵ_{\mathrm{\Delta }t}[j]ϵ_{\mathrm{\Delta }t}[k])`$ $`=`$ $`\sigma ^2\mathrm{\Delta }t\delta _{j,k}`$
$`E(\omega _{\mathrm{\Delta }t}[j]\omega _{\mathrm{\Delta }t}[k])`$ $`=`$ $`\lambda ^2\mathrm{ln}\rho _{\mathrm{\Delta }t}[|jk|]`$
where

$$\rho _{\mathrm{\Delta }t}[\mathrm{}]=\{\begin{array}{cc}\frac{T}{(|\mathrm{}|+1)\mathrm{\Delta }t}\hfill & \text{when}|\mathrm{}|\frac{T}{\mathrm{\Delta }t}1\hfill \\ 1\hfill & \text{otherwise}\hfill \end{array}$$

$`T`$ is the integral time beyond which the multifractal scaling ceases to hold. $`ϵ_{\mathrm{\Delta }t}[k]`$ is of mean $`0`$ and for the variance of the process to converge when $`\mathrm{\Delta }t0`$, one has to choose:

$$E(\omega _{\mathrm{\Delta }t}[k])=Var(\omega _{\mathrm{\Delta }t}[k])=\lambda ^2\mathrm{ln}\frac{T}{\mathrm{\Delta }t}.$$

Under these conditions, it was shown in  that the even moments of the process are, in the limit $`\mathrm{\Delta }t0`$, given by:

$$M(2p,\mathrm{})=K_{2p}\left(\frac{\mathrm{}}{T}\right)^{p2p(p1)\lambda ^2}.$$

with

$$K_{2p}=(2p1)!!(\sigma ^2T)^p_0^1𝑑u_1\mathrm{}_0^1𝑑u_p\underset{i<j}{}|u_iu_j|^{4\lambda ^2}.$$

The explicit computation of the above integral can be found in , that shows that the moments are finite only for $`q<q^{}=2+1/2\lambda ^2`$, beyond which these moments are infinite for all $`\mathrm{}`$. \>From the definition of $`\zeta _q`$, we obtain:

$$\zeta _q=\frac{q}{2}(1(q2)\lambda ^2),q=2p.$$

One can also compute the correlation of the local volatility, given by $`\sigma [k]\sigma e^{\omega _{\mathrm{\Delta }t}[k]}`$. It is easily shown that this correlation decays as a power-law of the time lag :

$$E(\sigma [k]^q\sigma [k+\mathrm{}]^q)E(\sigma [k]^q)^2\mathrm{}^{q^2\lambda ^2}.$$

Empirically, $`\lambda ^2`$ is found to be rather small, in the range $`0.020.1`$ .

## 3 A Skewed Multifractal Random Walk (smrw)

### 3.1 Definition

We now generalize the construction of  in order to account for the leverage effect, where the volatility is correlated with past price returns. We first consider the following discretized model (we omit the $`\mathrm{\Delta }t`$ in $`ϵ`$ and $`\omega `$ for simplicity):

$$\delta X_{\mathrm{\Delta }t}[k]=ϵ[k]e^{\stackrel{~}{\omega }[k]}\stackrel{~}{\omega }[k]\omega [k]\underset{i<k}{}K(i,k)ϵ[i],$$

where $`K(i,k)`$ is a certain kernel describing how the sign of the return at time $`i`$ affects the (log)-volatility at a later time $`k`$. We think of $`K(i,k)`$ as being positive, and therefore the minus sign accounts for the sign of the leverage effect. In order to preserve the multiscaling properties of the process, we choose the kernel $`K(i,j)`$ to decay as a power-law:

$$K(i,j)=\frac{K_0}{(ji)^\alpha \mathrm{\Delta }t^\beta }(j>i).$$

We now compute exactly the first moments of this process and then give an approximate formula (to first order in $`\mathrm{\Delta }t`$) for the higher moments.

### 3.2 The second moment

Let $`t=n\mathrm{\Delta }t`$, with $`1nN=\frac{T}{\mathrm{\Delta }t}`$. We find, using the notation $`\mathrm{}E_{ϵ,\omega }(\mathrm{})`$:

$`X_t^2`$ $`=`$ $`({\displaystyle \underset{i=0}{\overset{n1}{}}}\delta X_{\mathrm{\Delta }t}[i])^2={\displaystyle \underset{i=0}{\overset{n1}{}}}ϵ[i]^2e^{2\stackrel{~}{\omega }[i]}`$
$`=`$ $`\sigma ^2\mathrm{\Delta }t{\displaystyle \underset{i=0}{\overset{n1}{}}}e^{2\omega [i]}e^{2_{k<i}\frac{K_0}{(ik)^\alpha \mathrm{\Delta }t^\beta }ϵ[k]}`$
$`=`$ $`\sigma ^2\mathrm{\Delta }t{\displaystyle \underset{i=0}{\overset{n1}{}}}e^{2K_0^2\sigma ^2\mathrm{\Delta }t^{12\beta }_{k<i}\frac{1}{(ik)^{2\alpha }}}`$
We now make the following assumptions:

$`\alpha >{\displaystyle \frac{1}{2}}`$
$`\sigma ^2K_0^2\mathrm{\Delta }t^{12\beta }1`$ (2)
which, in particular, requires that $`\beta \frac{1}{2}`$ when $`\mathrm{\Delta }t0`$. The conditions (2) ensure that the sum $`_{k<i}\frac{1}{(ik)^{2\alpha }}=_{\mathrm{}=1}^{\mathrm{}}\frac{1}{\mathrm{}^{2\alpha }}`$ converges and that we can replace, to lowest order in $`K_0`$, the exponential terms in the expression for $`X_t^2`$, by $`1`$. Under these assumptions, we simply have:

$$X_t^2=\sigma ^2t,$$

as for the simple Brownian motion, or for the symmetric mrw considered in .

### 3.3 The third moment

The third moment, related to the skewness, explicitly reads:

$`X_t^3`$ $`=`$ $`{\displaystyle \underset{0i,j,k<n}{}}ϵ[i]ϵ[j]ϵ[k]e^{\stackrel{~}{\omega }[i]}e^{\stackrel{~}{\omega }[j]}e^{\stackrel{~}{\omega }[k]}`$
$`=`$ $`3{\displaystyle \underset{0i<j<n}{}}ϵ[i]ϵ[j]^2e^{\stackrel{~}{\omega }[i]+2\stackrel{~}{\omega }[j]}`$
$`=`$ $`3\sigma ^2\mathrm{\Delta }t{\displaystyle \underset{0i<j<n}{}}e^{\omega [i]+2\omega [j]}ϵ[i]e^{2K(i,j)ϵ[i]}`$
$`e^{_{k<i}\left(K(k,i)+2K(k,j)\right)ϵ[k]}e^{2_{i<k<j}K(k,j)ϵ[k]}`$
Within the conditions (2), the exponential terms (i.e. $`e^{K(k,j)ϵ_k}`$) can again be set to 1 to first order, provided the sums converge, which is the case whenever $`\alpha >1/2`$. Using the equality $`ϵe^{\lambda ϵ}=\lambda \sigma ^2e^{\frac{\lambda ^2\sigma ^2}{2}}`$, we find

$`X_t^3`$ $``$ $`6(\sigma ^2\mathrm{\Delta }t)^2{\displaystyle \underset{0i<j<n}{}}K(i,j)\rho _0^{\frac{\lambda ^2}{2}}\rho _{i,j}^{2\lambda ^2}`$ (3)
$`=`$ $`6K_0\sigma ^4T^{\frac{3\lambda ^2}{2}}\mathrm{\Delta }t^{\alpha \beta +\frac{\lambda ^2}{2}}{\displaystyle \underset{0i<j<n}{}}{\displaystyle \frac{\mathrm{\Delta }t}{((ji)\mathrm{\Delta }t)^\alpha }}{\displaystyle \frac{\mathrm{\Delta }t}{((ji+1)\mathrm{\Delta }t)^{2\lambda ^2}}}`$
$`\stackrel{\mathrm{\Delta }t0}{}`$ $`6K_0\sigma ^4T^{\frac{3\lambda ^2}{2}}\mathrm{\Delta }t^{\alpha \beta +\frac{\lambda ^2}{2}}{\displaystyle \underset{0<u<v<t}{}}{\displaystyle \frac{dudv}{|uv|^{\alpha +2\lambda ^2}}}`$
One easily proves that, under the condition $`\nu _2=\alpha +2\lambda ^2<1`$,

$$\underset{0<u<v<t}{}\frac{dudv}{|uv|^{\alpha +2\lambda ^2}}=\frac{t^{2\nu _2}}{(2\nu _2)(1\nu _2)},$$

so that our final result is:

$$X_t^3\frac{6K_0\sigma ^4T^{\frac{3\lambda ^2}{2}}}{(2\nu _2)(1\nu _2)}\mathrm{\Delta }t^\mu t^{2\nu _2},$$
(4)
where $`\mu =\alpha \beta +\frac{\lambda ^2}{2}`$. Therefore, the third moment grows as a power of the time lag $`t`$, with an exponent $`\zeta _3=2\nu _2`$. Due to the conditions (2), $`\mu `$ is always positive and therefore $`X_t^3`$ vanishes in the limit $`\mathrm{\Delta }t0`$. In other words, the skewness of the continuous-limit process disappears. In practice, however, the elementary time scale $`\mathrm{\Delta }t`$, beyond which price increments are uncorrelated, is equal to several seconds, even in futures markets, and to several minutes in stocks markets. \>From a theoretical perspective, however, it would be interesting to construct a multifractal process which remains skewed in the continuous time limit. We will come back to this point in section 4.

### 3.4 Higher order moments

As $`q`$ grows, the algebra of the explicit computation of the $`qth`$ moment rapidly becomes messy, and the exact result cannot be computed. However, as we show in appendix, one can still compute the dominating term in the limit $`\mathrm{\Delta }t0`$.

#### 3.4.1 Odd moments

We want to evaluate the generic moment:

$$X_t^{2p+1}=\underset{0i_1,\mathrm{},i_{2p+1}<n}{}ϵ[i_1]\mathrm{}ϵ[i_{2p+1}]e^{\stackrel{~}{\omega }[i_1]}\mathrm{}e^{\stackrel{~}{\omega }[i_{2p+1}]}.$$
(5)
The sum over $`2p+1`$ variables contains many different terms but, as we justify in the Appendix, the dominating one is:

$$𝒟_{2p+1}=\underset{0j_1,\mathrm{},j_{p+1}<n}{}ϵ[j_1]ϵ[j_2]^2\mathrm{}ϵ[j_{p+1}]^2e^{\stackrel{~}{\omega }[j_1]+2\stackrel{~}{\omega }[j_2]+\mathrm{}+2\stackrel{~}{\omega }[j_{p+1}]}$$
(6)
where all the variables are paired two by two but one. For small $`K_0`$, we finally find:

$$X_t^{2p+1}_{2p+1}K_0\mathrm{\Delta }t^\mu t^{p+12p^2\lambda ^2\alpha },$$
(7)
where $`_{2p+1}`$ is a positive numerical prefactor, and $`\mu `$ was defined above. Therefore, we find that all odd moments scale as a power of the time lag $`t`$, with $`\zeta _q=p+12p^2\lambda ^2\alpha `$ ($`q=2p+1`$). In this sense, the model considered is multifractal for odd moments as well.

#### 3.4.2 Even moments

We now want to evaluate the generic even moment

$$X_t^{2p}=\underset{0i_1,\mathrm{},i_{2p}<n}{}ϵ[i_1]\mathrm{}ϵ[i_{2p}]e^{\stackrel{~}{\omega }[i_1]+\mathrm{}+\stackrel{~}{\omega }[i_{2p}]}$$

Using similar arguments, detailed in the appendix, we can show that the dominating term is simply:

$$𝒟_{2p}=\underset{0j_1,\mathrm{},j_p<n}{}ϵ[j_1]^2\mathrm{}ϵ[j_p]^2e^{2\stackrel{~}{\omega }[j_1]+\mathrm{}+2\stackrel{~}{\omega }[j_p]}$$

and that we retrieve, for $`K_0`$ small, the result of :

$$X_t^{2p}_{2p}t^{p2p(p1)\lambda ^2}.$$

Hence, we find that a power-law leverage correlation function, introduced as we have proposed, does not affect even moments (to lowest order in $`K_0`$), but give to the odd moments a multifractal behaviour.

## 4 An extension to slowly decaying kernels

As we have noticed above, a condition for the convergence of sums appearing in the exponentials is $`\alpha >1/2`$. A way to relax this is to add a large time cut-off to the power-law kernel $`K(i,j)`$. More precisely, we set:

$$K(i,j)=\frac{K_0e^{\mathrm{\Gamma }(ji)}}{(ji)^\alpha \mathrm{\Delta }t^\beta }(j>i),$$

where $`\mathrm{\Gamma }\mathrm{\Delta }t/T1`$. We now consider the case condition $`\alpha <1/2`$, but impose:

$$\mathrm{\Gamma }^{2\alpha 1}K_0^2\mathrm{\Delta }t^{12\beta }1.$$
(8)
This ensures that the sum in the exponential terms is convergent and remains small. Following the same route as above, the third moment reads

$$X_t^3\stackrel{\mathrm{\Delta }t0}{}6K_0\sigma ^4T^{\frac{3\lambda ^2}{2}}\mathrm{\Delta }t^\mu \underset{0<u<v<t}{}\frac{e^{\frac{|uv|}{T}}dudv}{|uv|^{\alpha +2\lambda ^2}}.$$

Note that the double integral over $`u,v`$ converges provided $`\alpha +2\lambda ^2<1`$. Under the condition $`tT`$, the exponential term can be dropped and we recover the same result as above, Eq. (4), with $`\mu =\alpha \beta +\frac{\lambda ^2}{2}`$. Unfortunately, (8) with $`\mathrm{\Gamma }\mathrm{\Delta }t/T`$ impose that $`\alpha \beta 0`$ in the limit $`\mathrm{\Delta }t0`$. Therefore, the exponent $`\mu `$ is strictly positive and the skewness again disappears in the continuous time limit. On the other hand, the possibility of choosing $`\alpha <1/2`$ leads to the interesting possibility of a skewness that grows with time. Indeed, the third moment (given by Eq. (4)), divided by the third power of the volatility, scales as $`t^{1/2\nu _2}`$, with $`\nu _2=\alpha +2\lambda ^2`$. Small values of $`\alpha `$ therefore leads to a growing skewness (at least for $`tT`$), whereas $`\alpha >1/2`$ necessarily leads to a skewness that decays with time.

## 5 Numerical results

### 5.1 Methodology

#### 5.1.1 Simulation

In order to generate correlated Gaussian random variables (the $`\omega `$’s), we follow a quite standard procedure explained in , based on the utilization of the Fast Fourier Transform .
The major difficulty is to efficiently compute the skewness convolution $`_{i<k}K(i,k)ϵ[i]`$. We first generate $`2n`$ independent Gaussian variables $`\eta [i]`$, and periodize our series as:

$$ϵ[i]=\{\begin{array}{cc}\eta [i]\hfill & \text{if }0i<k\hfill \\ \eta [i+2n]\hfill & \text{if }k2ni<0\hfill \end{array}$$

and we define

$$\overline{K}(i,k)=\{\begin{array}{cc}K(i,k)\hfill & \text{if }0<kin\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$

We can now evaluate the skewness convolution by $`_{i=k2n}^{k1}\overline{K}(i,k)ϵ[i]`$, again using Fast Fourier Transforms. The use of $`n`$ ‘pseudo’ past variables and the cut-off imposed to $`\overline{K}(i,k)`$ are essential to avoid spurious correlations between past volatilities and future returns. The systematic use of the Fast Fourier Transform makes the whole procedure quite efficient and allows us to simulate long series of our process (in the applications we used series of $`2^{13}`$ variables but longer series can be obtained).

#### 5.1.2 Evaluation of the moments

We have evaluated the moments by averaging over a large number of realizations of the process. This is not such an easy task, because the quantities we deal with can be small and the convergence can be very slow due to the existence of long-range correlations in the variables $`\{\omega [i]\}`$, and the slow power-law decay of the kernel $`K(i,j)`$. In order to eliminate spurious skewness effects and to be more accurate, we have systematically averaged over the processes generated by the variables $`\{ϵ[i],\omega [i]\}`$ and by the variables $`\{ϵ[i],\omega [i]\}`$. In the absence of asymmetry (i.e. when $`K(i,j)0`$), the second process is exactly the “mirror image” of the first one and therefore all odd moments are strictly zero. Conversely in the presence of asymmetry, the second process is no longer the “mirror image” of the first one (see Fig 1), and one finds non zero odd moments. By proceeding this way, we are sure that any non zero result can indeed be attributed to the leverage kernel.

### 5.2 Some results

The results we present are based on 10 000 independent realizations of the process. A total of $`2^{13}`$ variables were used for each series. The parameters we used are $`T=4000\mathrm{\Delta }t`$, $`\lambda =0.175`$, corresponding to $`\lambda ^20.03`$. We chose two different values of the kernel exponent: $`\alpha =0.3`$ with $`\beta =\alpha +\lambda ^2/2`$ and $`\alpha =0.6`$ with $`\beta =0.45`$. The strength of the leverage effect is taken to be $`K_0=0.1`$, small enough to trust our approximate formulas. We have evaluated the different moments for times ranging from $`t=10\mathrm{\Delta }t`$ to $`t=1300\mathrm{\Delta }t`$, i.e. in a region where $`\mathrm{\Delta }ttT`$. We have plotted the logarithm of the different moments with as a function of the logarithm of the time lag. We expect to obtain a straight line with a slope given by $`\zeta _q`$.

We indeed observe a nice power-law scaling for the second and fourth moment of the process (Fig. 2), with a perfect adequation between the theoretical prediction and the simulations. The agreement is again very good for the third moment (Fig. 3), in particular for $`\alpha =0.3`$. There is a systematic difference for $`\alpha =0.6`$, that tends to disappear for large time lags. This is due to the fact that corrections to our asymptotic formulas tend to be stronger when $`\alpha `$ increases. Indeed, in the above calculations (3), we have set:

$$\underset{0i<j<n}{}\frac{\mathrm{\Delta }t}{((ji)\mathrm{\Delta }t)^\alpha }\frac{\mathrm{\Delta }t}{((ji+1)\mathrm{\Delta }t)^{2\lambda ^2}}\underset{0<u<v<t}{}\frac{dudv}{|uv|^{\alpha +2\lambda ^2}}.$$

Although this is correct in the limit $`\mathrm{\Delta }t/t0`$, there are corrections of the order of $`(\mathrm{\Delta }t/t)^{1\nu _2}`$, which persist when $`\nu _21`$. Let us note that this finite time correction is very small for the even moments, since in this case $`\nu _2`$ is replaced by $`4\lambda ^2`$, which is very small (see Fig. 4).

## 6 Comparison with empirical data

The comparison between the predictions of the (symmetric) mrw and real financial data was discussed in details in . In particular, the way to calibrate the parameters $`\lambda ^2`$ and $`T`$ is to study the correlation function of the log-volatility, that should behave as the log of the time lag up to time $`T`$, with a slope which is precisely $`\lambda ^2`$. While $`\lambda ^2`$ is quite well determined (its typical value is 0.03), the ‘integral time’ $`T`$ is much less precisely calibrated (it is found to be $`2`$ years). The elementary time interval $`\mathrm{\Delta }t`$ is given by the chosen discretization of the time series, provided it is larger than the correlation time of the returns (otherwise the assumption that $`ϵ[i]ϵ[j]=\delta _{i,j}`$ is not justified). For example, for daily data, $`\mathrm{\Delta }t=1`$ day.

As far as skewness is concerned, the following ‘leverage’ correlation function was proposed and studied in :

$$(i,j)=\frac{\delta X[i]\delta X[j]^2}{\left(\delta X[k]^2\right)^2}i<j.$$

Using the above definition we find, in the limit $`K_0^2\sigma ^2\mathrm{\Delta }t1`$, the following relation:

$$(i,j)2\left(\frac{T}{\mathrm{\Delta }t}\right)^{\frac{3\lambda ^2}{2}}\frac{K(i,j)}{|ij|^{2\lambda ^2}},$$
(9)
that allows in principle to calibrate $`K(i,j)`$ using empirical data, once $`\lambda ^2`$ and $`T`$ are fixed. (Note that only the combination $`K_0T^{\frac{3\lambda ^2}{2}}`$ appears in the above formulae, so that the uncertainty on $`T`$ does not transpire on the skewness). The leverage correlation function is however found to be, in the case of stock indices, close to a pure exponential in $`|ij|`$, with a decay time of ten days or so , in contrast with Eq. (9). For individual stocks, however, the decay time is much longer. The introduction of well defined time scale would ruin the scaling properties of the mrw process. We choose to preserve the scaling properties of the process for three reasons: (a) a strict multifractal skewed process is an interesting process in its own right, that may have applications outside finance (for example in turbulence ), (b) the scaling properties of the smrw might enable one to find exact analytic formulas for, e.g. probability distributions, and (c) the slow decay of the leverage function allows one to obtain long time persistent skews, and is useful in the context of option pricing (see section 8 below). In other words, even if the historical time series do not exhibit a long ranged leverage correlation function, the implied distributions, consistent with option smiles, might do so.

## 7 Volatility asymmetry

In , an interesting time asymmetry in financial series was detected using a wavelet analysis (see also  for related work). According to these authors, “the volatility at large scales influences causally (in the future) the volatility at shorter scales”, whereas the converse effect is much weaker. This finding is actually deeply related to the leverage effect  and our model, with its extended time correlations between returns and volatility, should be able to reproduce it, as we show now. One can quantify the effect reported in  by computing a correlation between volatilities at different scales $`n`$ and $`m`$, defined as:

$`C(n,m)`$ $`=`$ $`\left({\displaystyle \underset{i=0}{\overset{n1}{}}}\delta X[i]\right)^2\left({\displaystyle \underset{j=n}{\overset{n+m}{}}}\delta X[j]\right)^2`$
$`=`$ $`{\displaystyle \underset{i=0}{\overset{n1}{}}}{\displaystyle \underset{j=n}{\overset{n+m}{}}}\delta X[i]^2\delta X[j]^2+2{\displaystyle \underset{0i_1<i_2<n}{}}{\displaystyle \underset{j=n}{\overset{n+m}{}}}\delta X[i_1]\delta X[i_2]\delta X[j]^2.`$
Although this quantity is not exactly the one considered in , one expects that it behaves very similarly. The first term $`C_1`$ is easy to compute and is equal, in the limit $`\mathrm{\Delta }t0`$, to:

$`C_1`$ $`=`$ $`\sigma ^4T^{4\lambda ^2}{\displaystyle _{u=0}^\tau }{\displaystyle _{v=\tau }^{\tau +\tau ^{}}}{\displaystyle \frac{dudv}{|uv|^{4\lambda ^2}}}`$
$`=`$ $`{\displaystyle \frac{\sigma ^4T^{4\lambda ^2}}{(14\lambda ^2)(24\lambda ^2)}}[(\tau +\tau ^{})^{24\lambda ^2}\tau ^{24\lambda ^2}\tau _{}^{}{}_{}{}^{24\lambda ^2}],`$
where $`\tau =n\mathrm{\Delta }t`$ and $`\tau ^{}=m\mathrm{\Delta }t`$. This result is symmetric in $`\tau `$ and $`\tau ^{}`$, and therefore cannot explain the asymmetry. The second term, on the other hand, is equal to:

$`C_2`$ $`=`$ $`2\sigma ^6K_0^2T^{4\lambda ^2}\mathrm{\Delta }t^{\lambda ^2+2\alpha 2\beta }{\displaystyle _{0u_1<u_2<\tau v<\tau +\tau ^{}}}{\displaystyle \frac{2}{(vu_2)^\alpha }}`$
$`\left[{\displaystyle \frac{1}{(u_2u_1)^\alpha }}+{\displaystyle \frac{2}{(vu_1)^\alpha }}\right]{\displaystyle \frac{du_1du_2dv}{(u_2u_1)^{\lambda ^2}(vu_1)^{2\lambda ^2}(vu_2)^{2\lambda ^2}}}.`$
We now change variables and set:

$`u_1`$ $`=`$ $`x_1\tau `$
$`u_2`$ $`=`$ $`x_2\tau `$
$`v`$ $`=`$ $`\tau +y\tau ^{},`$
and assume first that $`\tau ^{}=ϵ\tau `$ with $`ϵ1`$, i.e. we compare future short scale volatilities with past large scale volatilities. One finds:

$$C_2K_0^2ϵ\tau ^{32\alpha 5\lambda ^2}.$$

Now, the opposite case where $`\tau =ϵ\tau ^{}`$ with $`ϵ1`$, i.e. where we compare future large scale volatilities with past small scale volatilities leads to:

$$C_2K_0^2ϵ^{2\alpha \lambda ^2}\tau ^{32\alpha 5\lambda ^2}.$$

The ratio of the second indicator to the first is therefore equal to $`ϵ^{1\alpha \lambda ^2}`$, which goes to zero provided that $`\alpha +\lambda ^2<1`$. In this case, the asymmetry indeed has the sign found in , that is, the correlation between past volatility at short scale and future volatility at large scale is indeed weaker than the correlation between past volatility at large scale and future volatility at short scale.

## 8 An application to option pricing: skewed volatility smile

It is a well known fact that Black-Scholes hypotheses are not satisfied in practice. In particular, fat tail effects and volatility clustering are absent in the Black-Scholes world. This has important consequences on option pricing (particularly for exotic options) and hedging. One particularly clear symptom is the so called volatility smile: implied volatility from option prices is not constant but varies across both strike and maturity, defining a volatility surface. Although this effect is widely known and studied by both academics and practitioners, a detailed understanding of the implied volatility surface is still one of the major challenge of modern finance.
Several authors  have given simple approximate formulas for option prices (or for the resulting volatility smile) when the asset is modelled by an arbitrary stochastic process, which capture the basic mechanisms that lead to non trivial smiles. The idea is to consider that the distribution of the price increments is weakly non Gaussian. Using a truncated cumulant expansion , one finds corrections to the constant volatility case in terms of the skewness and kurtosis of the underlying process (see also  for equivalent results but in a different language). For example, when the log-price difference between times $`0`$ and $`t`$ has a variance $`\sigma _t`$, skewness $`\kappa _{3,t}`$ and kurtosis $`\kappa _{4,t}`$, the volatility smile $`v_t`$ is approximately given by <sup>1</sup><sup>1</sup>1The formula given in  is in fact slightly different in that it neglects terms in $`\sigma _t`$ compared to the moneyness $`d_t`$, which we have kept below.:

$`v_t`$ $`=`$ $`\sigma _0\left[1+{\displaystyle \frac{\kappa _{3,t}}{3!}}(2\sigma _td_t){\displaystyle \frac{\kappa _{4,t}}{4!}}(1d_t^2+3d_t\sigma _t3\sigma _t^2)\right]`$ (11)
$`d_t`$ $`=`$ $`{\displaystyle \frac{\mathrm{log}(S_0/K)+rt+\sigma _t^2/2}{\sigma _t}},`$
where $`S_0`$ is the current price of the underlying, $`\sigma _0`$ is (true) volatility, $`K`$ the strike and $`r`$ the interest rate. In the limit where $`\sigma _t,rt1`$ and $`S_0KS_0`$, the above formula boils down to the formula obtained in  in the context of an additive model:

$`v_t`$ $`=`$ $`\sigma _0\left[1{\displaystyle \frac{\kappa _{3,t}}{3!}}d_t{\displaystyle \frac{\kappa _{4,t}}{4!}}(1d_t^2)\right]`$ (12)
$`d_t`$ $`=`$ $`{\displaystyle \frac{S_0K}{S_0\sigma _t}},`$
Since our process displays both anomalous skewness and kurtosis, we expect to observe this smile effect. The above formula is interesting as it gives an explicit expression of the option smiles for the smrw in terms of its known skewness and kurtosis. It is therefore important to compare the result (11) with an exact smile obtained by Monte-Carlo simulations (Fig.5), in order to test the accuracy of the above formula. We model the returns of the underlying by a geometric smrw:

$$\frac{\delta S_t}{S_t}=\mu \delta t+\sigma _0\delta X_t$$
(13)
where $`X_t`$ is a smrw. We then evaluate the price of an (european) option as an unconditional average of the final payoff $`(S_TK)_+=\mathrm{max}(S_TK,0)`$ over many realizations of our process with $`\mu =r`$. This assumes that

* the option is hedged using the Black-Scholes $`\mathrm{\Delta }`$ , which is a good approximation for the minimum variance hedging strategies. In this case, the ‘risk neutral’ process is indeed given by Eq. (13). Other hedging schemes are in fact possible, and would lead to small corrections to the option price, but this is beyond the aim of the present analysis, which is to test Eq. (11) in the framework within which it is established.
* one does not attempt to measure the instantaneous volatility and to condition the paths to start from this measured value. In other words, we study the unconditional, time independent volatility surface. Correspondingly, both the price and the hedge are independent of the current level of volatility.

In Fig.5 we have drawn the volatility smiles given by the Monte-Carlo simulation and the theoretical smile given by formula (11), with the theoretical skewness $`\kappa _{3,t}^{the}`$ and kurtosis $`\kappa _{4,t}^{the}`$ obtained from our above analysis:

$`\kappa _{3,t}^{the}`$ $`=`$ $`{\displaystyle \frac{6K_0\sigma T^{\frac{3\lambda ^2}{2}}t^{\frac{1}{2}\nu _2}}{(2\nu _2)(1\nu _2)}}\mathrm{\Delta }t^\mu `$ (14)
$`\kappa _{4,t}^{the}`$ $`=`$ $`{\displaystyle \frac{6T^{4\lambda ^2}t^{4\lambda ^2}}{(24\lambda ^2)(14\lambda ^2)}}3.`$ (15)
As already mentioned in section 4, the skewness grows with time when $`\alpha <\frac{1}{2}`$.<sup>2</sup><sup>2</sup>2At least in the regime $`tT`$. For $`tT`$, the multifractal features disappear and both the skewness and kurtosis vanish, as the process converges to a Gaussian. This can be of great interest for financial applications since it is often observed that the skewness of the option smile can persist for large maturities. We clearly see in (Fig.5) that the non Gaussian nature of our process induces an asymmetric smile and, visually, the agreement between simulations and theory is satisfying for maturities $`0.5,\mathrm{\hspace{0.17em}0.75},\mathrm{\hspace{0.17em}1}`$ and strikes between $`80`$ and $`120`$. The case of the shortest maturity is less convincing, as expected from a cumulant expansion, that is in principle only valid for long enough maturities. In order to be more precise we then have performed a parametric fit of these curves using Eq. (11) to obtain implied values of the skewness ($`\kappa _3^{imp}`$) and the kurtosis ($`\kappa _4^{imp}`$). We have chosen to fit the curves over the whole interval ($`[75:125]`$) of strikes. It is not obvious that this is the best choice because the implied volatility for deep in-the-money options (small strike) is known to be very sensitive to errors (because of a small Vega). We have also computed the empirical skewness and kurtosis from the simulations. We report these results in Tab.1 and Tab.2.

The agreement between the theoretical and implied values is only fair, even if the order of magnitude is correct. The implied kurtosis is smaller than the theoretical one for small maturities with a better agreement for large maturities, whereas the implied skewness tends to depart from its theoretical value for large maturities. Note that the skewness increases with maturity, at variance with most theoretical models in which it decreases as $`1/\sqrt{t}`$. Note also that for a fixed range of strikes, the range of relative moneyness decreases with maturity and the skewness becomes more difficult to ascertain precisely. In order to track the origin of the observed discrepancy<sup>3</sup><sup>3</sup>3Note that the agreement between empirical and theoretical volatility is much better (see Fig. 5)., we have directly tested the cumulant expansion on the cumulative distribution of our process. We show in Fig. 6 the difference between the numerical cumulative distribution of the smrw and the Gaussian distribution for two different time lags, $`t=0.25`$ and $`t=1`$, for $`\alpha =0.6`$.

We see that the cumulant expansion truncated at second order is only qualitatively correct in our case, since we observe systematic deviations between the empirical values and the predictions. This is particularly true when one uses the exact theoretical cumulants, and suggests that higher order cumulants cannot be neglected. This is in fact expected, since higher order cumulants, for example even ones, only decay with time as $`\kappa _{2p,t}t^{2p(p1)\lambda ^2}`$ for the smrw, instead of $`t^{p/2}`$ for iid returns. For example, $`\kappa _{6,t}/\kappa _{4,t}t^{8\lambda ^2}t^{0.25}`$ in our case. (Actually, cumulants even diverge beyond a certain order $`p^{}=1+1/4\lambda ^2`$, see ). Of course, the fit of the cumulative distribution with the implied values extracted from the option smile is better, since these values try to correct the inadequacy of the truncated expansion. The conclusion of this study is that a simple cumulant expansion for the smrw only gives a qualitative description of option smiles. In order to get quantitative results, one should use Monte-Carlo paths. However, the scaling properties of the model might allow one to find exact analytic solutions in some cases. Finally, we want to show many typical smile shapes can be obtained with our process by tuning its degree of asymmetry (measured by the parameter $`K_0`$). We have shown in (Fig.7) three smiles obtained for different values of $`K_0`$: $`0.1,0.5,1`$. In this example a more pronounced skewness can be obtained. The term structure of the skewness and of the curvature of the smile can be chosen by playing with the parameters $`\alpha `$ and $`\lambda ^2`$.

## 9 Conclusion and prospects

In this paper, we have generalized the construction of the multifractal random walk (mrw) of  to take into account the asymmetric character of the financial returns. Our first motivation was to introduce in this class of models the observed correlation between past returns and future volatility, in such a way that the scale invariance properties of the mrw are preserved. This correlation leads to a power-law skewness of the returns, and is in fact deeply connected with the so-called “harch” effect, or “causal cascade” observed in . Such a process can be very useful in finance since it captures most of the stylized facts of the returns: non Gaussian shape at short scales, volatility clustering, long range dependence, and, with the present work, tunable skewness. Possible applications in finance are numerous: pricing of options, risk management, etc. As illustrated above, the model considered here leads to very versatile shapes of volatility surfaces with anomalous term structure. In particular, one can account for smiles with skewness that is constant or even increases with maturity. This would not be possible without a slowly decaying return-volatility (leverage) correlation. From a more academic point of view, we have not been able to construct a skewed multifractal process which remains skewed in the continuous time limit. Whether this is a fundamental limitation, or that alternative constructions are possible, is an open question. It might be possible to marry the retarded volatility model of  with the mrw and obtain a skewed process in the continuous time limit. Another promising path would be to study the precise connections between the present model and the recent “causal cascade” construction of Calvet and Fisher, that extends the work presented in  and might allow to introduce the skew in a satisfactory way. We hope to answer these questions in future work .
Acknowledgments: We thank E. Bacry and J-F. Muzy for discussions at an early stage of this work, and especially J-F. Muzy for having pointed out to us the link with the “causal cascade” effect. We also thank L. Calvet for comments about this work, A. Matacz and M. Potters for many discussions on the leverage effect and on volatility smiles.

## Technical appendix

The computation of the $`2p+1^{th}`$ moment of $`X_t`$ involves terms of the form

$$\underset{0i_1,\mathrm{},i_q<n}{}ϵ[i_1]^{\gamma _1}\mathrm{}ϵ[i_q]^{\gamma _q}e^{\gamma _1\stackrel{~}{\omega }[i_1]+\mathrm{}+\gamma _q\stackrel{~}{\omega }[i_q]},$$

with

$`\gamma _i`$ $``$ $`𝒩^{}`$
$`{\displaystyle \underset{i}{}}\gamma _i`$ $`=`$ $`2p+1`$
Writing more explicitly $`\stackrel{~}{\omega }[i]`$, we can expand these terms as:

$$\underset{0i_1,\mathrm{},i_q<n}{}ϵ[i_1]^{\gamma _1}e^{𝒦_1ϵ[i_1]}\mathrm{}ϵ[i_q]^{\gamma _q}e^{𝒦_qϵ[i_q]}e^{\gamma _1\omega [i_1]+\mathrm{}+\gamma _q\omega [i_q]}\underset{\mathrm{}<i<n,ii_1,\mathrm{},i_q}{}e^{𝒦_iϵ[i]}.$$

The exponential terms $`ϵ[i]^{\gamma _i}e^{𝒦_𝒾ϵ[i]}`$ follow from the definition of the $`\stackrel{~}{\omega }_i`$ and $`𝒦_i`$ are complicated prefactors. The above expression looks quite involved but as we search for the leading term, we can use (2) in order to simplify it. In particular, to the first order in $`K_0`$, $`e^{𝒦_iϵ_i}1`$. We also have

$`ϵ[i]^{\gamma _i}e^{𝒦_𝒾ϵ[i]}`$ $``$ $`ϵ[i]^{\gamma _i}(1+𝒦_𝒾ϵ[i]+\mathrm{}`$
$``$ $`\{\begin{array}{cc}ϵ[i]^{\gamma _i}\hfill & \text{if }\gamma _i\text{ is even}\hfill \\ 𝒦_𝒾ϵ[i]^{1+\gamma _i}\hfill & \text{if }\gamma _i\text{ is odd}\hfill \end{array}`$
Due to (2), we deduce that the lowest order term in $`K_0`$ is found when a maximum number of $`\gamma _i`$ are even, that is all $`\gamma _i`$ except one (because we look at the $`2p+1^{th}`$ moment) are even. Finally, using basic properties of the Gaussian vector,{$`\omega _i`$}, we find that:

$$e^{\gamma _1\omega [i_1]+\mathrm{}+\gamma _q\omega [i_q]}=\underset{1k<\mathrm{}q}{}\left(\frac{T}{(1+|i_ki_{\mathrm{}}|)\mathrm{\Delta }t}\right)^{\gamma _k\gamma _{\mathrm{}}\lambda ^2}\left(\frac{T}{\mathrm{\Delta }t}\right)^{\lambda ^2(\frac{1}{2}{\scriptscriptstyle \gamma _i^2}{\scriptscriptstyle \gamma _i})}.$$

We want to know which combination of the $`\{\gamma _i\}`$ leads to the dominating term. One way to do this is to study the scaling with $`\mathrm{\Delta }t`$ of the term corresponding to a given combination of $`\{\gamma _i\}`$. Combining the previous computations, we find that the exponent of $`\mathrm{\Delta }t`$ is:

$$\mathrm{\Delta }t^{1+(\alpha \beta )+\lambda ^2+p(1+2\lambda ^2)q\frac{1}{2}\lambda ^2_{i=1}^q\gamma _i^2}.$$

In the limit $`\mathrm{\Delta }t0`$, only the smallest exponent will contribute. The only parameters we can adjust are $`q`$ and $`_{i=1}^q\gamma _i^2`$, with the constraint $`_{i=1}^q\gamma _i=2p+1`$. The two extrema are:

* $`q=2,\gamma _{i_1}=2p,\gamma _{i_2}=1`$ which corresponds to an exponent:

$$\mu _2=(\alpha \beta )+\frac{\lambda ^2}{2}+(p1)(12p\lambda ^2).$$
* $`q=p+1,\gamma _{i_1}=\mathrm{}=\gamma _{i_p}=2,\gamma _{i_{p+1}}=1`$ which corresponds to an exponent:

$$\mu _{p+1}=(\alpha \beta )+\frac{\lambda ^2}{2}.$$

The difference $`\mu _2\mu _{p+1}=(p1)(12p\lambda ^2)`$ is positive as soon as $`p<p^{}=\frac{1}{2\lambda ^2}`$. Therefore, we can conclude that the behaviour of the first $`2p^{}+1`$ moments is dominated by only one term, corresponding to pairing all points except one. A similar discussion can be made for even moments as well. The above argument is clearer in the mrw model , where the calculation is easier and can be made without any further assumptions. In this model, the odd moments are zero and the even moments have a scaling in $`\mathrm{\Delta }t`$ given by the exponent $`p(1+2\lambda ^2)q\frac{\lambda ^2}{2}_{i=1}^q\gamma _i^2`$. We immediately see that if $`q=p\text{ and }\gamma _i=2`$, this exponent is 0, which is the condition for the existence of a non trivial continuous limit .


format: text

text:
# A simple model for dynamical melting of moving vortex lattices interacting with periodic pinning

## I introduction

The study of dynamical phases of lattice structures moving in a periodic medium has received a great deal of attention lately due to its relevance to physical systems such as adsorbed atomic layers in boundary lubrication , vortices in Josephson Junction Arrays (JJA) and in films with artificial defect lattices.

Numerical studies of driven vortices in clean superconducting films interacting with a periodic lattice of columnar defects (CD) and in JJA  find dynamical melting of a moving vortex-lattice (VL) into a moving vortex liquid. In these systems the moving vortices order in a VL at low temperatures when the driving force magnitude is sufficiently large. Melting of the VL is found to take place at a temperature that, for a given vortex density, decreases with the driving force magnitude, approaching the equilibrium VL melting temperature when this magnitude approaches infinity. These studies also show that the vortices center of mass (CM) velocity is in general not parallel to the driving force. This leads to anisotropic V-I curves, with voltages both longitudinal and transverse to the direction of the applied current.

In this paper we show that dynamical melting lines and anisotropic V-I curves similar to those found in the above mentioned numerical studies can be obtained using a simple model for the vortex dynamics.

In driven vortex systems interacting with random or periodic pinning, moving vortex lattices are expected to occur at low temperatures and in the limit of very large CM velocities, because the pinning potential is averaged in the direction of motion, as first pointed out by Schmid and Hauger . For periodic pinning, Ref. shows that in this limit the moving vortex lattice can be commensurate or incommensurate with the periodic pinning potential, depending on the direction of motion and on the vortex density. Commensurate lattices occur only for motion along directions of high symmetry of the pinning potential, such as , ,  and \[-1,1\] for square symmetry, if the vortex density is not too high . For motion along other directions, the moving vortex lattices are incommensurate. Recent numerical simulations find that both moving commensurate and incommensurate lattices may undergo dynamical melting into a moving liquid. However, dynamical melting of commensurate lattices is much more complicated than that of the incommensurate ones. The reason is that the moving commensurate lattices are pinned in channels running along the direction of motion (transverse pinning), and melting into a liquid occurs simultaneously with depinning from the channels. In this paper we consider only dynamical melting of incommensurate lattices.

Dynamical melting of vortex lattices is still poorly understood. It is unclear from the simulation results reported so far if the change from lattice to liquid spatial symmetries is a true phase transition, with accompanying singularities in correlation functions, or merely a crossover. In this paper we do not offer any insight into this fundamental question. Our aim is to estimate the dynamic melting temperature based on a well know method for the equilibrium melting temperature, namely Lindemann’s criterion. This method was proposed in a classical paper by Houghton, Pelcovitz and Sudbø. They estimate the equilibrium VL melting temperature for clean bulk superconductors by applying Lindemann’s criterion to the mean-square vortex-line displacement caused by thermal fluctuations. A simple model to estimate the same temperature was proposed later by Frey, Nelson and Fisher . It considers a single vortex line trapped in a harmonic potential representing its interaction with the other vortex lines in the lattice, assumed straight and fixed at their equilibrium lattice positions. We refer to this model here as the equilibrium cage model. The equilibrium melting temperature is estimated by applying Lindemann’s criterion to the mean-square thermal fluctuations of the caged vortex line. The result agrees with that of Ref., if the harmonic potential is properly chosen. This simple idea can also be applied to two-dimensional VL melting in thin films. In this case the equilibrium cage model consists of a point vortex tied to a spring, representing its interaction with the other point vortices, assumed fixed at their lattice positions. The melting temperature is obtained by applying Lindemann’s criterion to the caged vortex mean square displacement due to thermal fluctuations . The resulting melting temperature agrees with the one obtained from dislocation unbinding theory, if the spring constant is properly chosen. Our model, which we call the dynamical cage model, generalizes the two-dimensional equilibrium one for a moving vortex lattice. The dynamical melting temperature is estimated by applying Lindemann’s criterion to the mean-square vortex displacements caused by thermal fluctuations and by motion in the periodic pinning potential. Application of Lindemann’s criterion to estimate the dynamical melting line of vortices in bulk superconductors interacting with random pinning was first proposed by Koshelev and Vinokur.

We start from Langevin’s equations describing two-dimensional vortices at temperature $`T`$, interacting between themselves and with a periodic pinning potential, and driven by an external force. We consider vortex motion at large CM velocities and make the following assumptions: i) The moving vortices order in a lattice, incommensurate with the periodic pinning potential. ii) The vortex lattice flows elastically. This means that there are no dislocations and that vortex-vortex interactions can be approximated by harmonic forces. We then propose a mean-field like treatment of the resulting equations of motion. This approximation reduces the problem to that of a single vortex tied to a spring, both moving with the CM velocity, and interacting with the periodic pinning potential. We find that the dynamical cage model is governed by two equations of motion. One for the vortex displacement in the frame moving with the vortex lattice CM velocity (CM frame), and another relating the driving force and the CM velocity. The vortex displacement equation is found to be nonlinear, unlike that for the equilibrium cage model, due to the vortex-pinning interaction. These equations are solved by a perturbation method, similar to that introduced by Schmid and Hauger  for weak random pinning forces. It consists in expanding the time-dependent periodic pinning force (in the CM frame) in powers of the vortex displacement. We find that the perturbation expansion is valid at large CM velocities for both weak and strong pinning. From these solutions we obtain expressions for the vortex mean square displacement and for the relationship between the CM velocity and the driving force. The dynamical melting temperature is obtained by applying Lindemann’s criterion to the caged vortex mean-square displacement. These results are applied to a square pinning array, and compared to numerical simulation ones.

This paper is organized as follows. In Sec. II we derive the dynamical cage model equations of motion and their perturbation theory solution. In Sec. III. we apply these results to strong pinning by a square array of pinning centers. The conclusions of the paper are presented in Sec. IV. Details of the perturbation calculation are given in the Appendix.

## II Dynamical Cage Model

We consider $`N_v`$ two-dimensional vortices at temperature $`T`$ interacting between themselves and with a periodic potential produced by an array of pinning centers, and driven by a force $`𝐟_d`$. The equations of motion are

$$\eta \frac{d𝐫_j}{dt}=𝐟_d+𝐅_j^{\text{v-v}}+𝐅^{\text{v-p}}(𝐫_j)+𝚪_j,$$
(1)
where $`j=1,\mathrm{},N_v`$, $`\eta `$ is the friction coefficient,

$$𝐅_j^{\text{v-v}}=\underset{ij=1}{\overset{N_v}{}}_jU^{\text{v-v}}(𝐫_j𝐫_i),$$
(2)
is the force of interaction with other vortices, $`U^{\text{v-v}}(𝐫)`$ being the vortex-vortex interaction potential in two dimensions,

$$𝐅^{\text{v-p}}(𝐫_j)=\underset{𝐐}{}(i𝐐)U_𝐐^{\text{v-p}}e^{i𝐐𝐫_j},$$
(3)
is the force of interaction with the periodic pinning array, $`𝐐`$ denotes the pinning array reciprocal lattice vectors, and $`U_𝐐^{\text{v-p}}`$ is the Fourier transform of the vortex-single pin interaction potential, and $`𝚪_j`$ is the random force, satisfying

$$\mathrm{\Gamma }_{j\alpha }(t)\mathrm{\Gamma }_{l\beta }(t^{})=2k_BT\eta \delta _{j,l}\delta _{\alpha ,\beta }\delta (tt^{}),$$
(4)
where $`\mathrm{}`$ denotes average over the random force distribution.

We assume that the vortices center of mass (CM) moves with constant velocity $`𝐯`$, defined, as usual, by

$$𝐯=\frac{1}{N_v}\underset{j=1}{\overset{N_v}{}}\frac{d𝐫_j(t)}{dt}.$$
(5)
Using Eq. (1), we find that

$$\eta 𝐯=𝐟_d+\frac{1}{N_v}\underset{j=1}{\overset{N_v}{}}𝐅^{\text{v-p}}(𝐫_j(t)).$$
(6)
To derive Eq. (6) we use the fact that for large $`N_v`$, the random force term is negligible, since $`_{i=1}^{N_v}𝚪_j\sqrt{N_v}`$. According to Eq. (6), in order that $`𝐯`$ is constant $`𝐟_d`$ must depend on time as well as on the random force. To obtain physical results we average Eq. (6) over time and over the random force distribution, that is

$$\eta 𝐯=𝐅_d+\frac{1}{N_v}\underset{i=1}{\overset{N_v}{}}\frac{1}{\tau }_0^\tau 𝐅^{\text{v-p}}\left(𝐫_j(t)\right),$$
(7)
where $`𝐅_d`$, a constant vector, denotes the average of $`𝐟_d`$ over time and over the random force distribution and $`\tau `$ is a time large compared with the characteristic times in Eq. (1). We interpret $`𝐅_d`$ as the force due to the applied current, so that Eqs. (7) give the V-I curves.

Now we consider a vortex-lattice flowing elastically and transform to the CM frame. We write the vortex positions as

$$𝐫_j(t)=𝐑_j+𝐯t+𝐮_j(t),$$
(8)
where $`𝐑_j`$ and $`𝐮_j(t)`$ ($`j=1,\mathrm{},N_v`$) are, respectively, the vortex equilibrium positions and displacements from equilibrium in the CM frame. In elastic flow the displacements $`𝐮_j(t)`$ are small and the vortex-vortex interactions can be approximated by harmonic forces. Our dynamical cage model considers a particular vortex, and substitute its elastic interaction with the other vortices in the lattice by a spring of constant $`\kappa `$. The equation of motion for this vortex in the CM frame is found to be

$`\eta {\displaystyle \frac{d𝐮_l(t)}{dt}}`$ $`=`$ $`\kappa 𝐮_l(t)+𝐅^{\text{v-p}}(𝐑_l+𝐮_l(t)+𝐯t)`$ (10)
$`+𝐟_d\eta 𝐯+𝚪_l.`$
The relationship between the spring constant and the vortex-vortex interaction potential will be presented shortly. Our dynamical cage model approach is similar to the mean-field theory introduced by Fisher for sliding charge-density waves interacting with random pinning.

Our objective here is to use Eq. (10), together with Eq. (7), to calculate the vortex mean-square displacement required for Lindemann’s criterion, and to obtain the $`𝐯`$ vs. $`𝐅_d`$ relationship (V-I curves). First we note that Eq. (10) depends on the particular vortex $`l`$, through the $`𝐅^{\text{v-p}}`$\- term. Since $`𝐅^{\text{v-p}}`$ is periodic in the pinning-array lattice, its argument in Eq. (10) can be reduced, at any instant $`t`$, to a position within the pinning-array primitive unit cell. As mentioned in Sec. I, we consider only moving incommensurate vortex lattices. In this case the reduced vortex positions are uniformly distributed over the unit cell area. Accordingly, we define the vortex mean square displacement as

$$u^2=\frac{1}{N_v}\underset{j=1}{\overset{N_v}{}}\frac{1}{\tau }_0^\tau 𝐮_j(t)^2.$$
(11)
According to the above considerations, our dynamical cage model is described by the set of equations

$`\eta {\displaystyle \frac{d𝐮_l(t)}{dt}}`$ $`=`$ $`\kappa 𝐮_l(t)+𝐅^{\text{v-p}}(𝐑_l+𝐮_l(t)+𝐯t)`$ (13)
$`{\displaystyle \frac{1}{N_v}}{\displaystyle \underset{j=1}{\overset{N_v}{}}}𝐅^{\text{v-p}}(𝐑_j+𝐮_j(t)+𝐯t)+𝚪_l,`$
and

$$\eta 𝐯=𝐅_d+\frac{1}{N_v}\underset{j=1}{\overset{N_v}{}}\frac{1}{\tau }_0^\tau 𝐅^{\text{v-p}}(𝐑_j+𝐮_𝐯t).$$
(14)
To obtain $`𝐮_l(t)`$ we solve Eq. (13) by the perturbation theory method introduced in Ref., which consists in expanding $`𝐅^{\text{v-p}}`$ in powers of $`𝐮_i(t)`$ and solving the resulting equations by iteration. We keep only terms to first order. The relationship between $`𝐯`$ and $`𝐅_d`$ are obtained from Eq. (14) by a similar expansion. The details of these calculation are given in the Appendix.

For the mean square displacement we obtain to leading order

$`u^2`$ $`=`$ $`{\displaystyle \frac{2k_BT}{\kappa }}+{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{Q^2|U_𝐐|^2}{\eta ^2(𝐐𝐯)^2+\kappa ^2}}`$ (16)
$`+{\displaystyle \frac{2k_BT}{\kappa }}{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{Q^4|U_𝐐|^2}{\eta ^2(𝐐𝐯)^2+4\kappa ^2}}.`$
To the same order of approximation the $`𝐯`$ vs. $`𝐅_d`$ relationship is

$`𝐯`$ $`=`$ $`{\displaystyle \frac{𝐅_d}{\eta }}{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{𝐐Q^2|U_𝐐|^2}{\eta ^2(𝐐𝐯)^2+\kappa ^2}}𝐐𝐯`$ (18)
$`{\displaystyle \frac{k_BT}{\kappa }}{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{𝐐Q^4|U_𝐐|^2}{\eta ^2(𝐐𝐯)^2+4\kappa ^2}}𝐐𝐯.`$
According to Lindemann’s criterion the VL melts when $`u^2=c_L^2a_v^2`$, where $`a_v`$ is the VL lattice parameter. Thus, our dynamical cage model only makes sense as long as $`u^2<c_L^2a_v^2`$. This requires that in Eq. (16) both $`k_BT/\kappa `$, and the velocity dependent terms are sufficiently small compared to $`a_v^2`$. The latter are so for large enough $`v`$. To state more precisely the conditions under which the perturbation expansion in Eq. (16) is valid, we restrict our considerations to periodic pinning forces that can be represented by Eq. (3) with only a few non-zero $`U_𝐐`$. Typically those for $`𝐐`$ corresponding to nearest and next-nearest neighbor reciprocal lattice points. In these cases $`Q2\pi /a_p`$, $`a_p`$ being the pinning array lattice parameter. We also assume that the number of vortices is comparable to the number of pins, so that $`a_va_p`$. Under these assumptions, the condition for the smallness of both velocity dependent terms in Eq. (16) depends on the relative strengths of the pinning and elastic forces. For weak pinning, defined as $`|U_𝐐|\kappa a_p^2/2\pi `$, these terms are always small, independent of $`𝐯`$, so that the effects of vortex

motion on $`u^2`$ are negligible. For strong pinning, that is $`|U_𝐐|\kappa a_p^2/2\pi `$, the velocity dependent terms in Eq. (16) are small for all $`𝐐`$, if $`\eta v\mathrm{cos}\vartheta _𝐐\kappa a_p/2\pi `$, where $`\vartheta _𝐐`$ is the angle between the particular $`𝐐`$ and $`𝐯`$. Similar considerations apply to the perturbation expansion for $`𝐯`$, Eq. (18). For strong pinning, the perturbation expansion breaks down if, for some $`𝐐`$ entering the sums in Eq. (16), $`𝐐𝐯=0`$ ($`\mathrm{cos}\vartheta _𝐐=0`$). In this case the velocity dependent terms in Eq. (16) are not small. This breakdown can also be seen from Eq. (10). If $`𝐐𝐯=0`$ the pining force $`𝐅^{\text{v-p}}(𝐑_l+𝐮_l(t)+𝐯t)`$ has a static component, periodic in the direction perpendicular to $`𝐯`$, which gives rise to a static vortex displacement comparable to $`a_v`$ for strong pinning. Such a displacement indicates that our assumption that the moving VL is incommensurate with the pinning array lattice breaks down for motion along the directions for which $`𝐐𝐯=0`$. There are only a few directions of motion for which $`𝐐𝐯=0`$ since, as mentioned above, only a few $`𝐐`$ enter the sums in Eq. (16). We note that the terms in sum in Eq. (18) with $`𝐐𝐯=0`$ do not contribute to $`𝐯`$, so that there is no breakdown in the perturbation expansion for this quantity.

In the limit of very large $`v`$ ($`v\mathrm{}`$) our dynamical cage model reduces to the static one, since vortex motion averages out the pinning potential (we assume $`𝐐𝐯0`$). In this limit $`u^2`$ reduces to the equilibrium mean square displacement of a vortex tied to a spring, $`u^2=\frac{2k_BT}{\kappa }`$. Application of Lindemann’s criterion gives for the melting temperature $`k_BT_m=\kappa a_v^2c_L^2/2`$. To compare this estimate of $`T_m`$ with the dislocation unbinding theory one, $`\kappa `$ is identified with $`\kappa =(^2U^{vv}(r)/r^2)_{r=a_v}`$, where $`U^{vv}(r)=(\varphi _0^2/8\pi \mathrm{\Lambda })\mathrm{ln}(r/\xi )`$ is the 2D vortex-vortex interaction potential ($`\varphi _0`$ is the flux quantum and $`\mathrm{\Lambda }`$ the film effective penetration depth ). The result is $`\kappa =\varphi _0^2/8\pi \mathrm{\Lambda }a_v^2`$ and $`k_BT_m=c_L^2\varphi _0^2/16\pi \mathrm{\Lambda }`$. This melting temperature agrees with the dislocation unbinding theory one  if $`c_L0.12`$. From here on we adopt for $`\kappa `$ and $`c_L`$ the above mentioned values.

Application of Lindemann’s criterion to Eq. (16) gives the velocity dependent melting temperature $`T_{dm}(𝐯)`$ as

$`k_BT_{dm}(𝐯)={\displaystyle \frac{\frac{\kappa a_v^2c_L^2}{2}k_BT_{pp}(𝐯)}{1+_𝐐\frac{Q^4|U_𝐐|^2}{\eta ^2(𝐐𝐯)^2+4\kappa ^2}}},`$ (19)
where

$$k_BT_{pp}(𝐯)\frac{\kappa }{2}\underset{𝐐}{}\frac{Q^2|U_𝐐|^2}{\eta ^2(𝐐𝐯)^2+\kappa ^2}$$
(20)
is identified with an effective temperature, resulting from vortex vibrations due to the periodic pinning force, that adds to the thermodynamic one in the expression for $`u^2`$. The temperature $`T_{pp}`$ is similar to the ’shaking temperature’ introduced by Koshelev and Vinokur  for random pinning. It arises from a $`T`$-independent contribution to $`u^2`$, the second term in the right hand side of Eq. (16). There is also a $`T`$-dependent contribution to $`u^2`$ in Eq. (16) that cannot be identified with an effective temperature. However, this term leads to a small correction in $`T_{dm}(𝐯)`$, since the denominator in Eq. (19) is close to one. Thus, the dynamical melting condition is essentially that $`k_B(T_{dm}(𝐯)+T_{pp}(𝐯))=\kappa a_v^2c_L^2/2`$. The $`𝐯`$ dependence of $`T_{dm}(𝐯)`$, Eq. (19), is such that $`T_{dm}(𝐯)`$ decreases with decreasing $`v`$.

The V-I curves follow from the relationship between $`𝐯`$ and $`𝐅_d`$, Eq. (18). In the $`v\mathrm{}`$ limit $`𝐯=𝐅_d/\eta `$. For finite $`v`$, vortex motion is no longer along the drive direction, leading to anisotropic V-I curves. The $`𝐯`$ vs. $`𝐅_d`$ relationship also allows us to obtain the melting temperature as a function of the driving force, which is more usual.

In the next section we study in detail the dynamical cage model predictions by applying them to a typical periodic pinning potential and carrying out numerical calculations.

## III  square pinning-lattice

We consider a square pinning array with lattice parameter $`a_p`$, and assume that the vortex-pinning interaction force is given by Eq. (3) with the $`U_𝐐^{\text{v-p}}`$ chosen as

$`U_𝐐^{\text{v-p}}=\{\begin{array}{cc}U_1\hfill & \text{for }𝐐=\pm \frac{2\pi }{a_p}\widehat{x}\text{}\pm \frac{2\pi }{a_p}\widehat{y}\text{;}\hfill \\ U_2\hfill & \text{for }𝐐=\pm \frac{2\pi }{a_p}(\widehat{x}\pm \widehat{y})\text{;}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}`$ (24)
The $`x`$ and $`y`$ axis are along the pinning lattice $`[1,0]`$, $`[0,1]`$ directions, respectively, as shown in Fig. 1(a). We choose $`U_2=U_1/2`$. This gives the pinning potential shown in Fig. 1(b). We set $`U_1/\kappa a_p^2=1/2\pi =0.16`$. As discussed in Sec. II, this value corresponds to strong pinning.

We choose to present the numerical results for the expressions derived in Sec. II as functions of the driving force $`𝐅_d`$, instead of $`𝐯`$. To do so we invert Eq. (18) by an iterative method. We start at a high value of the driving force magnitude, $`F_d=50\kappa a_p`$, and a given orientation of $`𝐅_d`$ with respect to the $`x`$-axis, $`\alpha `$ (Fig. 1(a)), and approximate the solution by the result for infinite drive $`𝐯=𝐅_d/\eta `$. Then $`F_d`$ is decreased by small steps, keeping $`\alpha `$ fixed. At each step, $`𝐯`$ in the right hand side of Eq. (18) is replaced by the $`𝐯`$ obtained in the previous step. From the $`𝐯`$ vs. $`𝐅_d`$ relationship thus obtained, the $`u`$ and $`T_{dm}`$ vs. $`F_d`$ curves follow using Eqs. (16) and (19). The results are shown in Figs. 2-5. In the figures shown here the following units are used. i) Force: $`F_{max}=8\pi U_1/a_p=`$ magnitude of the maximum force exerted by the square pinning potential defined in Eqs. (24). ii) Velocities: $`F_{max}/\eta `$. iii) Temperatures: $`T_m=`$ VL equilibrium melting temperature (see Sec. II).

In Fig. 2 we show $`𝐯`$ and $`u`$ as functions of $`F_d`$ for several $`\alpha `$ at $`T=0`$. The horizontal lines in the bottom panel of Figs. 2(a) and 2(b) indicate the $`u=c_La_v`$ lines for a vortex density corresponding to two vortices per pin ($`B=2B_\varphi `$), with $`c_L=0.12`$, at which the VL melts, according to Lindemann’s criterion. The region of validity of our dynamical cage model is below this line. The $`v`$ vs. $`F_d`$ curves in the top panels of Figs. 2(a) and 2(b) show that the CM velocity is not parallel to the driving force, which leads to anisotropic V-I curves. This anisotropy is more clearly seen by plotting the direction of motion, $`\theta `$, as a function of the driving force magnitude, as shown in Fig. 3, for fixed directions of drive, $`\alpha `$. These plots show that, depending on $`\alpha `$, the moving VL is attracted towards the  or  directions. This behavior is similar to that seen in recent numerical simulations . Above the $`c_L=0.12`$ lines in the bottom panel of Figs. 2(a) and 2(b) we find a discontinuous jump in $`u`$, accompanied by the trapping of moving caged vortex along  for $`0<\alpha <32.5^o`$ \[Fig. 2(b)\], or along  for $`32.5^o<\alpha <45^o`$ \[Fig. 2(a)\]. This trapping is similar to transverse pinning observed in numerical simulations and JJA experiments. This is interesting because it indicates that, even though our dynamical cage model is not strictly valid in the region where trapping occurs, it contains physical ingredients capable of describing the phenomenon. We also note that the jump in $`u`$ occurs because for vortex motion along  and  there are terms in the sums in Eq. (16) for which $`𝐐𝐯=0`$. According to the discussion in Sec. II, the perturbation expansion for $`u`$ breaks down and a large contribution to $`u`$ results from the terms in the sum in Eq. (16) with $`𝐐𝐯=0`$.

The finite-$`T`$ behavior of the results is shown in Fig. 4. We find that the $`𝐯`$ vs. $`F_d`$ curves change little with $`T`$ as illustrated in Fig. 4(a). We also find that in $`u^2`$ the $`T`$-dependent term in Eq. (16) gives a small contribution. This is so because the $`T`$-dependence in $`u`$ and $`𝐯`$ results from the third terms in the right hand side of Eqs. (16) and (18), which are, essentially, the product of two quantities which are small in the region of validity of our dynamical cage model.

The dynamical melting lines $`T_{dm}`$ vs. $`F_d`$ are shown in the top panel of Fig. 5 for several $`\alpha `$. In the bottom panel, the dynamical phase diagram $`F_d`$ vs $`\alpha `$ for $`T=0`$ is shown. The $`F_d`$ vs. $`\alpha `$ dynamical melting line in the bottom panel corresponds to the intersections of the $`T_{dm}`$ vs. $`F_d`$ curves shown in Fig. 5 top with the $`T=0`$ line. We find that for other temperatures $`T<T_m`$ the $`F_d`$ vs $`\alpha `$ phase diagrams are similar.

## IV conclusion

In conclusion then, we introduce a simple dynamical model for vortex lattices flowing elastically in the presence of periodic pinning. The model is solved by a perturbation method that allows the calculation of the dynamical melting curves, using Lindemann’s criterion,

and of the V-I curves. This approach is valid for periodic pinning potentials in general, and predicts that the moving vortex-lattice melting temperature varies significantly with the CM velocity only for strong pinning. The perturbation solution is applied to a square pinning lattice. The results show qualitative agreement with numerical simulation ones. Namely: i) the predicted dynamical melting temperature approaches the thermal equilibrium one in the limit of very large center of mass velocities, and decreases with the velocity. The model predicts that this decrease is, essentially, as $`v^2`$. ii) The center of mass velocity is not parallel to the driving force in general, leading to anisotropic V-I curves. The vortex lattice is attracted towards the  or  direction, depending on the direction of drive. Transverse pinning, with the vortex lattice motion restricted to the  or  directions, is found outside the region where the model is strictly valid.

For other periodic pinning potentials, and for strong pinning, we expect similar dynamical melting lines and anisotropic V-I curves.

It is significant that the model predicts anisotropic V-I curves so closely resembling the numerical simulation ones. These curves follow from the relationship between the CM velocity and the driving force, Eq. (13). The only approximation used in our solution of this equation is the mean-field treatment of Langevin’s equations for the whole vortex lattice. Our results for the V-I curves show that this approximation is justifiable.

The model can be also applied along the lines described here to study the effects of weak random potentials, either from material defects or from imperfections in the pinning array, on the melting lines and V-I curves. It can also be generalized to study transverse pinning more rigorously. Work along these lines is under way and will be reported elsewhere.

###### Acknowledgements.

Research supported in part by CNPq, FAPERJ, FUJB and FACEPE. We thank J. A. Aguiar and L. R. E. Cabral for stimulating discussions and useful suggestions.

## A

To obtain the perturbation expansion of Eq. (13) we write it, using Eq. (3), as

$`\eta {\displaystyle \frac{d𝐮_l(t)}{dt}}`$ $`=`$ $`\kappa 𝐮_l(t)+{\displaystyle \underset{𝐐}{}}𝐅_𝐐^{\text{v-p}}e^{i\omega _𝐐t}\times `$ (A2)
$`{\displaystyle \frac{1}{N_v}}{\displaystyle \underset{j=1}{\overset{N_v}{}}}(N_v\delta _{j,l}1)e^{i𝐐[𝐑_j+𝐮_j(t)]}+𝚪_l,`$
where $`𝐅_𝐐^{\text{v-p}}i𝐐U_𝐐^{\text{v-p}}`$ and $`\omega _𝐐=𝐐𝐯`$. Expanding $`e^{i𝐐𝐮_j}`$ to first order, and carrying out the Fourier transformation in time, assuming periodic boundary conditions in the time interval $`(0,\tau )`$, justifiable for large $`\tau `$ ($`\tau a_p/v`$), we find

$`𝐮_l(\omega )`$ $`=`$ $`𝐮_l^T(\omega )+{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{𝐅_𝐐^{\text{v-p}}}{i\eta \omega +\kappa }}{\displaystyle \frac{1}{N_v}}{\displaystyle \underset{j=1}{\overset{N_v}{}}}(N_v\delta _{j,l}1)\times `$ (A4)
$`e^{i𝐐𝐑_j}\left[\delta _{\omega ,\omega _𝐐}+i𝐐𝐮_j(\omega +\omega _𝐐)\right],`$
where $`𝐮_l^T(\omega )`$ is the thermal displacement

$$𝐮_l^T(\omega )=\frac{\mathrm{\Gamma }_l(\omega )}{i\eta \omega +\kappa },$$
(A5)
$`\mathrm{\Gamma }_l(\omega )`$ is the Fourier transform of the random force and $`\delta _K`$ denotes the Kronecker delta. The first order solution for $`𝐮_l(\omega )`$ is obtained by neglecting the linear term in Eq. (A4). We find

$$𝐮_l^{(1)}(\omega )=𝐮_l^T(\omega )+\underset{𝐐}{}\frac{𝐅_𝐐^{\text{v-p}}e^{i𝐐𝐑_l}}{i\eta \omega _𝐐+\kappa }\delta _{\omega ,\omega _𝐐}.$$
(A6)
To derive Eq. (A6) we use $`_{j=1}^{N_v}e^{i𝐐𝐑_j}=\delta _{𝐐,0}`$, since the vortex and pinning lattices are incommensurate.

The second order solution is obtained by substituting $`𝐮_j`$ in the linear term in the left hand side of Eq. (A4) by the first order solution. The result is

$`𝐮_l^{(2)}(\omega )`$ $`=`$ $`𝐮_l^{(1)}(\omega )+{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{𝐅_𝐐^{\text{v-p}}}{i\eta \omega +\kappa }}\times `$ (A9)
$`{\displaystyle \frac{1}{N_v}}{\displaystyle \underset{j=1}{\overset{N_v}{}}}(N_v\delta _{j,l}1)e^{i𝐐𝐑_j}i𝐐𝐮_j^T(\omega +\omega _𝐐)+`$
$`{\displaystyle \underset{𝐐,𝐐^{}}{}}{\displaystyle \frac{𝐅_𝐐^{\text{v-p}}i𝐐𝐅_𝐐^{}^{\text{v-p}}[e^{i𝐊𝐑_l}\delta _{𝐐,𝐐^{}}]}{(i\eta \omega +\kappa )(i\eta \omega _𝐐^{}+\kappa )}}\delta _{\omega ,\omega _𝐊},`$
where $`\omega _𝐐^{}=𝐐^{}𝐯`$, $`\omega _𝐊=𝐊𝐯`$, and $`𝐊=𝐐+𝐐^{}`$. With this solution, one can calculate the mean square displacement for a particular vortex $`l`$, $`u_l^2=\frac{1}{\tau }_0^\tau |𝐮_l(t)|^2`$, which up to the second order in the pinning potential is expressed by:

$`u_l^2`$ $`=`$ $`{\displaystyle \frac{2k_BT}{\kappa }}+{\displaystyle \underset{𝐐,𝐐^{}}{}}{\displaystyle \frac{𝐅_𝐐^{\text{v-p}}𝐅_𝐐^{}^{\text{v-p}}e^{i𝐊𝐑_l}}{\eta ^2\omega _𝐐^2+\kappa ^2}}\delta _{\omega _𝐐,\omega _𝐐^{}}+`$ (A13)
$`{\displaystyle \frac{2k_BT}{\kappa }}{\displaystyle \frac{1}{N_v}}{\displaystyle \underset{𝐐}{}}{\displaystyle \frac{Q^2𝐅_𝐐^{\text{v-p}}𝐅_𝐐^{}^{\text{v-p}}}{\eta ^2\omega _𝐐^2+4\kappa ^2}}\times `$
$`\left[(N_v2)e^{i𝐊𝐑_l}+\delta _{𝐐,𝐐^{}}\right].`$
Taking the average over the $`N_v`$ vortices and assuming $`N_v1`$, Eq. (16) is derived.

A similar procedure is used to calculate the driving force dependence of $`𝐯`$. Expanding Eq. (7) in the form

$`\eta 𝐯`$ $`=`$ $`𝐅_d+{\displaystyle \frac{1}{N_v}}{\displaystyle \underset{j=1}{\overset{N_v}{}}}{\displaystyle \frac{1}{\tau }}\times `$ (A15)
$`{\displaystyle _0^\tau }{\displaystyle \underset{𝐐}{}}𝐅_𝐐^{\text{v-p}}e^{i𝐐𝐑_j}e^{i\omega _𝐐t}\left[1+𝐐𝐮_j(t)\right],`$
and using the second order solution for $`u_l(t)`$ we obtain Eq. (18).