text:
# Re 𝐴₀ and Re 𝐴₂ from Quenched Lattice QCD

## 1 INTRODUCTION

At energies below the electroweak scale the weak interactions are described by local four-fermi operators multiplied by effective coupling constants, the Wilson coefficients. The formal framework to achieve this is the operator product expansion (OPE) which allows one to separate the calculation of a physical amplitude into two distinct parts: the short distance (perturbative) calculation of the Wilson coefficients and the long distance (generally non-perturbative) calculation of the hadronic matrix elements of the operators $`Q_i`$. We calculate on the lattice $`K\pi `$ and $`K0`$. This allows us to calculate the low energy constants in chiral perturbation  which, after incorporating the non-perturbative renormalization factors are then translated into $`K\pi \pi `$ matrix elements. The CP-PACS collaboration has also presented a very similar calculation at this meeting .

## 2 DETAILS OF THE SIMULATIONS

We have used the Wilson gauge action, quenched, at $`\beta =6.0`$ on a $`16^3\times 32`$ lattice which corresponds to an inverse lattice spacing $`a^1=1.922\mathrm{GeV}`$. The domain wall fermion height $`M_5=1.8`$ and fifth dimension $`L_s=16`$ give a residual symmetry breaking $`m_{\mathrm{res}}=0.001244\mathrm{MeV}`$ ; 400 configurations separated by 10000 heat-bath sweeps were used in this analysis. $`K\pi `$ matrix elements were calculated in the $`SU(3)_{\mathrm{flavor}}`$ limit for 5 light quark masses $`m_f=0.01,0.02,0.03,0.04,0.05`$. Since the $`K0`$ matrix elements vanish in the $`SU(3)`$ limit these matrix elements were calculated with non-degenerate quark propagators for 10 mass combinations subject to the constraint $`m_sm_d=0.01,0.02,0.03,0.04`$. We have also calculated the so called eye diagrams with an active charm quark for $`m_f=0.1,0.2,0.3,0.4`$ ( the physical charm quark is around 0.5). However, the analysis for charm-in is still in progress; in this presentation we concentrate on the case with 3-active flavors wherein charm is integrated out assuming it is very heavy. The calculation took about 4 months on 800 Gflops (peak). Quark propagators were calculated using the conjugate gradient method with a stopping residual of $`10^8`$ with periodic and anti-periodic boundary conditions which amounts to doubling the lattice size in time direction. The two wall source propagators at $`t_K=5`$ and $`t_\pi =27`$ were fixed to Coulomb gauge. For eye diagrams we employed random wall sources spread over time slices $`t=1417`$ with 2 hits per configuration. Dividing the three-point correlation functions by the wall-wall pseudoscalar-pseudoscalar correlation function yields the desired matrix elements up to a factor of $`2m_\pi `$ which is determined from a covariant fit to the wall-point two-point function in the range $`t=1220`$ for each mass.

## 3 CALCULATION OF LOW ENERGY CONSTANTS

Since our results unambiguously show that Re$`A_0`$ and Re$`A_2`$ come essentially from the current-current operators (recall these have the largest Wilson coefficients) we will concentrate on these operators from now on. Quenched chiral perturbation theory predicts

$$m_\pi ^2=a_\pi (m_f+m_{\mathrm{res}})\left[1\delta \mathrm{ln}\left(\frac{a_\pi (m_f+m_{\mathrm{res}})}{\mathrm{\Lambda }_{Q\chi pt}^2}\right)\right]$$

We find a quenched chiral logarithm coefficient $`\delta =0.05(2)`$ which has a negligible contribution in our matrix element calculation. Unlike the quenched chiral logarithms, the conventional logarithms coming from quenched chiral perturbation theory induce large corrections to the $`\mathrm{\Delta }I=3/2`$ $`K\pi `$ matrix element as can be seen in Figure 1. We fit these amplitudes to

$$b_1^{(27,1)}m_M^2[1(\delta +\frac{6m_M^2}{(4\pi f)^2})\mathrm{ln}(\frac{m_M^2}{\mathrm{\Lambda }^2})]+b_2^{(27,1)}m_M^4$$

where $`m_M^2=a_\pi (m_f+m_{res})`$, $`f=137\mathrm{MeV}`$, $`\mathrm{\Lambda }=1\mathrm{GeV}`$. The conventional chiral logarithm $`m_M^2\mathrm{ln}(m_M^2/\mathrm{\Lambda }^2)`$ is almost linear over the mass range we have used so the fitting routine cannot distinguish this term from the linear term if we leave the coefficient of the logarithm as a free parameter. Since the large coefficient -6 of the logarithm makes the contribution of this term comparable to the contribution of the linear term omitting this term would change $`b_1^{(27,1)}`$ by almost a factor of two. The quenched chiral log contribution is very small.

$`\mathrm{\Delta }I=1/2`$ $`K\pi `$ matrix elements mix with $`\overline{s}d`$ with a power divergent coefficient $`(m_s+m_d)/a^2`$. We define a subtracted matrix element $`\pi ^+|Q_i^{(8,1)}|K^+_{sub}`$ by

$$\pi ^+|Q_i^{(8,1)}|K^++\eta _i(m_s+m_d)\pi ^+|\overline{s}d|K^+$$

where $`\eta _i`$ is obtained from a linear fit to $`\frac{0|Q_i|K}{0|\overline{s}\gamma ^5d|K}`$. For an explanation of this subtraction of the power divergence we refer the reader to . The quenched chiral perturbation theory corrections to $`\pi ^+|Q_i^{1/2}|K^+_{sub}`$ are not known. As can be seen in Figure 2 our data is consistent with a linear fit $`c_{0,i}+c_{1,i}m_f`$ with the slope $`c_{1,i}`$ determining the low energy constants $`\alpha _i^{1/2}`$ and the intercept $`c_{0,i}`$ arising from residual chiral symmetry breaking.

## 4  Re $`A_0`$ AND Re $`A_2`$

We use chiral perturbation theory to compute the lattice $`K\pi \pi `$ matrix elements. Using non-perturbative Z factors we obtain the continuum matrix elements which are then multiplied by Wilson coefficients to yield the physical amplitudes. We present an extrapolation to the Kaon mass scale to lowest order in chiral perturbation theory and a second extrapolation which includes one loop logarithmic effects. We multiply the pseudoscalar masses by $`\xi `$ so that for $`\xi 0`$ the chiral perturbation theory extrapolation is increasingly accurate but we need the extrapolation at $`\xi =1`$, the physical point. In figure 3 we present Re$`A_0`$ and Re$`A_2`$ as a function of the parameter $`\xi `$. The chiral logarithm correction for Re$`A_0`$ is large (about $`40\%`$). In addition one expects a large correction (not included here) coming from the tree level $`O(p^4)`$ terms necessary to cancel the dependence on the chiral perturbation theory scale $`\mathrm{\Lambda }_{\chi PT}`$.

If the Z factors and the Wilson coefficients were calculated to all orders in perturbation theory the physical amplitudes that we calculate would not depend on the scale $`\mu `$ where the transition between the lattice and the continuum operators is made. To a good approximation this is what we find, even though at $`\mu =1.51\mathrm{GeV}`$ one expects non-perturbative effects in the Z factors and at $`\mu =3.02\mathrm{GeV}`$ the discretization errors may be large. In Figure 4 we present the ratio Re$`A_0`$/Re$`A_2`$, the so called $`\mathrm{\Delta }I=1/2`$ rule which shows a large enhancement in the $`I=0`$ channel in accord with experiment (note, the chiral logarithm corrections largely cancel in the ratio so a large enhancement is seen for both extrapolation choices). The residual scale dependence in the physical amplitudes is slight (see Figure 4).

## 5 CONCLUSIONS

In conclusion, Re$`A_0`$, Re$`A_2`$ and especially the ratio Re$`A_0`$/Re$`A_2`$ were found reasonably close to the experimental values. We see this as an important success of the lattice method. However there were a number of major approximations in our calculation, the hardest to quantify is the use of quenched QCD. Also the chiral logarithms in quenched $`K\pi `$, $`\mathrm{\Delta }I=1/2`$ are not known and we have included only the logarithmic portion of the next-to-leading-order, 1-loop corrections in $`K\pi \pi `$ extrapolations.


format: text

text:
# Fluctuating-friction molecular motors

## 1 Introduction

Driven by the modern tools of molecular biology and by the opening perspectives of nanotechnology, there is currently a strong interest into understanding all the mechanisms by which systems at a molecular scale can efficiently convert chemical or light energy into mechanical energy. Such systems are often called “molecular motors”, and are at the root of biological processes such as muscle contraction, cell motility, and several intracellular transport processes . Contrary to ordinary engines, these molecular motors are conceived to work at a single temperature. Moreover they essentially exploit Brownian motion, converting its random behaviour into an ordered directional motion by means of some physical mechanism related with energy dissipation. An example of such a mechanism is the so-called “Brownian ratchet effect” . Artificial Brownian ratchets at microscopic scales have been proposed for a variety of technological applications, and in some cases they have been experimentally demonstrated . Very recently, the first steps toward engineering truly molecular motors have been taken .

Within this field, a specific line of investigation is aimed at developing physical models that capture the essential features of these molecular systems but that are simple enough to be easily understood and studied . Many models are focused on a very simple ideal system: a point-like particle undergoing overdamped Brownian motion under the effect of molecular friction and of a periodic energy potential $`U(x)`$, where $`x`$ is the particle coordinate. At a given constant absolute temperature $`T`$, friction and Brownian diffusion are characterized by a single quantity that can be taken to be the friction coefficient $`\eta `$ or equivalently the diffusion constant $`D=kT/\eta `$, where $`k`$ is Boltzmann constant. In this case, it is well known that even if $`U(x)`$ is asymmetric for $`xx`$, as for example in the case of a saw-tooth potential, no directional long-range motion can be induced unless the system is driven out of thermodynamic equilibrium by some mechanism .

One of the simplest proposed mechanisms to induce directional motion is a random switching of the system between two internal states characterized by different potentials, say $`U_1(x)`$ and $`U_2(x)`$  (see also  for a detailed two-state model of kinesin). In other words, the particle is subject to a fluctuating potential, where the fluctuations are described as sudden switches from $`U_1(x)`$ to $`U_2(x)`$ and backwards. These switches are not completely random, but obey to stochastic laws in which detailed balance and hence thermodynamic equilibrium is broken. For example, one must assign the probability per unit time $`I_{ij}(x)`$ of having a nonthermal forced transition from state $`i`$ to state $`j`$, associated with some input of free energy. Moreover, there is the probability per unit time of having spontaneous transitions associated with thermal equilibrium and therefore obeying to detailed balance. One may prove for this system that the correlated energy potential fluctuations associated with the random transitions may result into a nonzero average force $`F`$ acting on the particle and therefore into its drift at a constant average velocity $`v=F/\eta `$ . This long-range motion is the mechanical output of the motor, whereas the (chemical or optical) free energy input arises from the nonthermal state transitions. The direction of this motion is dictated by the asymmetry of one or both the potentials $`U_i(x)`$. However, a directional motion can be induced also by asymmetric transition rates .

The possibility of a difference in the friction coefficient (or diffusion constant) experienced by the particle in the two internal states, i.e., $`\eta _1\eta _2`$ (or $`D_1D_2`$), was considered in previous works but only with a passive role. In particular, if taken alone, this difference cannot give rise to a nonzero average force and velocity of the motor. In other words a fluctuation of friction in the presence of a constant force potential cannot power the motor.

In this paper we show that a stochastic fluctuation of friction induced by asymmetric transition rates actually can give rise to a nonzero average force $`F`$. Unexpectedly, however, this steady average force does not generate any long-range motion of the particle, i.e., $`v=0`$. More precisely, it can be shown that the average friction force experienced by the particle, owing to the fluctuations in the friction coefficient, does not vanish even for $`v=0`$. Therefore, a balance between the average nonzero force arising from the (constant) potential $`U(x)`$ and the average friction force is established at $`v=0`$. At first sight it seems therefore that there is no possibility at all of obtaining motion and mechanical work out of Brownian motion with fluctuating friction coefficient. We show in the following that this is not completely true.

Let us consider explicitly the two “external” bodies that interact with the particle and provide respectively the potential forces and the friction forces acting on the particle. In the following, we call these two bodies the “sources” of the potential and friction forces. Usually, in discussing molecular motor models, these bodies are considered as large systems that do not move. However, this picture is not always realistic. For example, in biological actomyosin motors the actin (“thin”) filaments (usually considered as the source of potential forces in the framework of these models) are displaced relative to myosin (“thick”) filaments by the action of many myosin motors. It is this displacement the actual output of the motor. Therefore, we are led naturally to consider also the motion of these “large” bodies, besides the particle itself, in analyzing the mechanical output of these systems. The whole motor system must actually be seen then as a three-body dynamical system.

Now, the bodies generating the opposite friction and potential forces acting on the particle are also, by reaction, experiencing opposite nonvanishing average forces. Therefore these two large bodies can be set in motion relative to each other, unless they are blocked by some constraint or they actually coincide. In other words, the motor particle undergoing Brownian motion with a fluctuating friction coefficient can induce a long-range relative displacement of the two bodies with which it is interacting. In the process, the particle itself may also undergo a directional long-range displacement.

This paper is organized as follows. In the next Section we describe a simple model of a three body-system that converts chemical (or optical) free energy into mechanical work by exploiting the Brownian motion of a particle that experiences stochastic sudden fluctuations of its friction coefficient between two values. In section III we solve analytically the equations of this system for a specific choice of the potential landscape and transition probabilities. We also performed some numerical studies to analyze more general cases. In section IV, in order to show how the idea of a fluctuating-friction molecular motor can be useful in the context of biological systems, we apply our model to the actomyosin system. Some concluding remarks are given in section V.

## 2 The model

Consider three bodies denoted as A, B, and C, corresponding to the motor active unit (the “particle”), the potential source, and the friction source, respectively (see figure 1). We start working in a frame of reference $``$ (we may consider it “fixed”, for the sake of simplicity) in which C is motionless while A and B move, their motion being described by coordinates $`x_A(t)`$ and $`x_B(t)`$. The motion of B is supposed to be deterministic, i.e., its Brownian fluctuations are neglected, for example because it is a comparatively larger body. For the time being, we assume that B moves at a constant velocity $`V_B`$, i.e., $`x_B=V_Bt`$. The motion of the small particle A is instead stochastic, due to Brownian fluctuations. Moreover, we assume that it is overdamped, i.e., inertia of A can be neglected. Therefore, to each internal state $`i=1,2`$ of the motor we may associate an evolving probability density $`f_i(x_A,t)`$. The interaction between A and B is described by a potential $`U_i(x_A,x_B)=U_i(x_A,V_Bt)`$ that is time-dependent, because B is moving. It is therefore convenient to switch to another inertial frame of reference $`_B`$ that is comoving with B, thus introducing a new relative coordinate $`x=x_Ax_B`$ in terms of which the potential $`U_i(x)`$ is stationary.

The stochastic dynamics of body A can then be described by the following pair of coupled Smoluchowski equations:

$$\frac{f_i}{t}+\frac{J_i}{x}=W_i(x),$$
(1)
where $`J_i(x)`$ is the probability current in the frame $`_B`$ and $`W_i(x)`$ is the net rate of transitions to state $`i`$ at position $`x`$. $`J_i`$ is given by

$$J_i=D_i\left(\frac{f_i}{x}+\frac{f_i}{kT}\frac{U_i}{x}\right)f_iV_B,$$
(2)
where $`D_i=kT/\eta _i`$ is the diffusion constant in state $`i`$. The first two terms in equation (2) are the standard diffusion and drift currents, respectively. The last term appears because we are working in a frame of reference, $`_B`$, that is moving with respect to the source of friction (normally, overdamped Brownian motion is studied in the privileged reference frame of the fluid providing the friction). Formally, it can be derived starting from equations (1) written in the fixed frame $``$, by applying the Galilei coordinate transformation $`x=x_AV_Bt`$ to the $`f_i/t`$ term. Physically, it describes the drift current induced by the drag force generated by the friction source body C that – in the frame $`_B`$ – is moving at constant speed $`V_B`$. This term would arise, for example, when describing the Brownian motion of a particle in a viscous fluid that is flowing at constant speed $`V_B`$.

For comparing the contributions of fluctuating friction and fluctuating potential, we allow here for both of them to depend on the internal state. For the transition rates we take

$$W_2(x)=W_1(x)=I(x)f_1(x)\frac{f_2(x)}{\tau },$$
(3)
with $`\tau `$ being the lifetime of state 2 (state 1 is assumed to be stable, while $`1/\tau `$ is the rate of spontaneous transitions $`21`$). The latter follows from assuming that state transitions are local, i.e., do not involve a variation of $`x`$, and thermal-induced transitions from level 1 to level 2 are negligible (because the level energy difference is much larger than $`kT`$). Note moreover that we are assuming $`I(x)`$ to depend only on the relative coordinate $`x=x_Ax_B`$.

Thus far we have made no hypothesis on the periodicity of $`U(x)`$ and $`I(x)`$. Actually, this is not a strictly necessary condition in our model, as we will discuss later. However, we assume in the following that the functions $`U_i(x)`$ and $`I(x)`$ are periodic, owing to the extended periodic structure of body B. Denoting then by $`L`$ the half-period in the variable $`x`$, we may limit our solution to the interval $`x[L,L]`$ by imposing the periodic boundary conditions $`f_i(L)=f_i(L)`$ and $`J_i(L)=J_i(L)`$. The distribution functions $`f_i`$ are normalized by the condition $`_i_L^Lf_i𝑑x=1`$.

In the stationary regime, when $`f_i/t=0`$, the condition $`W_1=W_2`$ combined with equations (1) implies that the total current $`J_t=J_1(x)+J_2(x)`$ is independent of $`x`$. This first-integral can be exploited to reduce equations (1) and (2) to a set of three first-order differential equations in the unknown functions $`f_1(x),f_2(x)`$, and $`J_2(x)`$, with $`J_t`$ playing the role of an eigenvalue. Then, the three remaining periodic boundary conditions together with the normalization condition uniquely determine the unknown functions and the eigenvalue $`J_t`$. In this way, the problem can be easily solved numerically for any functional shape of $`U_i(x)`$ and $`I(x)`$ and for a given velocity $`V_B`$. For the case of piecewise linear functions shown in figure 2, the problem can also be solved analytically, as discussed in section III.

Once solved the stationary random-walk problem, we are interested in computing the average force $`F_{AB}`$ that B exerts on A, given by

$$F_{AB}=\underset{i}{}_L^L\frac{dU_i}{dx}f_i𝑑x,$$
(4)
and the average friction force $`F_{AC}`$ arising in the interaction with C, given by

$$F_{AC}=\underset{i}{}\frac{kT}{D_i}_L^L\left(J_i+f_iV_B\right)𝑑x$$
(5)
(the latter is easily understood by considering that $`J_i+f_iV_B`$ is the current in the fixed frame $``$, i.e., the current relative to body C). By exploiting equations (2), it is straightforward to prove that $`F_{AB}+F_{AC}=0`$, reflecting the balance of all forces acting on the element A of the motor. Due to this balance, the particle A does not acquire any momentum (that it could not accumulate, since its inertia is negligible). On the other hand, by interacting with B and C, A mediates a continuous transfer of momentum from the friction source body C to the potential source body B. This transfer is expressed by the reaction force $`F_{BA}=F_{AB}=F_{AC}`$ acting on B. In the following, for brevity, we denote the latter just by $`F=F_{BA}`$. A convenient expression of $`F`$, obtained from equation (5), is the following:

$$F(V_B)=\frac{2kTL}{D_1}\left[J_t+\left(\frac{D_1}{D_2}1\right)\overline{J}_2\right]\overline{\eta }V_B,$$
(6)
where

$$\overline{J}_2=\frac{1}{2L}_L^LJ_2(x)𝑑x$$
(7)
is the average current in state 2 and

$$\overline{\eta }=kT_L^L\left(\frac{f_1(x)}{D_1}+\frac{f_2(x)}{D_2}\right)𝑑x$$
(8)
is the average friction coefficient.

The average relative velocity $`v`$ of the active element A in the frame $`_B`$ can be evaluated as $`v=2LJ_t`$. If we set $`V_B=0`$ and there are no potential fluctuations, i.e., $`U_1(x)=U_2(x)=U(x)`$, then $`v=J_t=0`$, for any functional form of $`U(x)`$. To prove it, let us introduce the function

$$g(x)=\left[D_1f_1(x)+D_2f_2(x)\right]e^{U(x)/kT}.$$
(9)
The periodic boundary conditions and the periodicity of $`U(x)`$ imply that $`g(L)=g(L)`$. However, from equations (2) we find the relation

$$\frac{dg}{dx}=J_te^{U(x)/kT},$$
(10)
which after an integration gives

$$J_t=\frac{g(L)g(L)}{_L^Le^{U(x)/kT}𝑑x}=0.$$
(11)
This means that there can be no long range displacement of the particle A with respect to the potential source for $`V_B=0`$, as anticipated in the Introduction. From equation (6) one can see, however, that the vanishing of $`J_t`$ does not imply the vanishing of the force $`F`$ if $`D_1D_2`$, as the average current $`\overline{J}_2`$ does not necessarily vanishes for $`V_B=0`$. In fact, as illustrated graphically in figure 3, the system may achieve a stationary nonequilibrium state in which there is a continuous cycling between state 1, in which the particle in the average moves one way, and state 2, in which in the average the particle moves the opposite way by the same amount. The net particle displacement vanishes but, owing to the different friction coefficients in the two states, the average friction force does not. The actual achievement of this stationary state is proved in section III for the particular choice of $`U(x)`$ and $`I(x)`$ shown in figure 2. In general, we found numerically that $`F`$ does not vanish if the nonthermal transition probability $`I(x)`$ is described by a periodic function that is shifted with respect to the potential $`U(x)`$, so that the resulting overall transition rates $`W_i(x)`$ are asymmetric. Note that if one sets $`V_B0`$, then $`J_t0`$, i.e., particle A will acquire a directional motion relative to the other two bodies.

Let us now come back to body B. Thus far we have simply assumed that B is moving at constant speed $`V_B`$, and we have thus determined the average force $`F=F(V_B)`$ that the motor develops on B itself. Now we can use this information to write a self-consistent equation for the (deterministic) motion of B at steady-state,

$$m_B\frac{dV_B}{dt}=F(V_B)\eta _pV_BF_{\mathrm{ext}}=0,$$
(12)
where $`\eta _p`$ is an additional friction coefficient associated with the motion of body B, and $`F_{\mathrm{ext}}`$ is the external load. In particular, setting $`V_B=0`$ in equation (12) we obtain the “stalling” force $`F_{\mathrm{ext}}=F(0)`$. The maximum velocity $`V_{\mathrm{max}}`$ is instead obtained by solving the force-balance equation with $`F_{\mathrm{ext}}=0`$.

Note that we are assuming here that body B responds only to the average force $`F(V_B)`$, and not to the instantaneous fluctuating force exchanged with A. As we said, this assumption is justified for example if B is much larger and therefore slower than A, so that it responds only to its time-averaged motion and ignores its fast Brownian fluctuations. Another case that justifies this approach occurs when B interacts simultaneously with many replicas of A and therefore it responds only to their total force (proportional to the average force of a single motor). However, it should be noted also that this simplifying assumption of our model is not essential to our main conclusions and that qualitatively similar results are to be expected when Brownian fluctuations of B are taken into account.

One final comment is in order. At first sight, one could be led to believe that by assuming from the start a nonzero velocity for body B we have artificially introduced a bias in the random-walk of A, and that it is only this bias that finally leads to the directional motion. That this is not the case is proved by the fact that the internal force $`F(V_B)`$ does not vanish for $`V_B=0`$. It is this force that defines the long-range motion direction. And it is this internal force that, if not counteracted by a stalling external load, eventually sets body B in motion.

## 3 Solutions

Let us now describe some specific solutions of our model. It is convenient to use dimensionless quantities, obtained from the corresponding dimensional ones by using $`L,L^2/D_2,kT`$ as units of length, time, and energy, respectively. We denote a dimensionless quantity by a tilde sign placed above the corresponding dimensional symbol. In particular $`\stackrel{~}{I}(x)=I(x)L^2/D_2`$, $`\stackrel{~}{F}=FL/kT`$, $`\stackrel{~}{D}_1=D_1/D_2`$, and $`\stackrel{~}{V}_B=V_BL/D_2`$. For definiteness, let us specialize to the case of triangular-wave potentials and square-wave transition probability shown in figure 2, i.e., $`\stackrel{~}{U}_i(\stackrel{~}{x})=[2\theta (\stackrel{~}{x})1]\stackrel{~}{u}_i\stackrel{~}{x}`$ and $`\stackrel{~}{I}(\stackrel{~}{x})=\stackrel{~}{I}_0\theta (\stackrel{~}{x})`$, where $`\theta (\stackrel{~}{x})`$ is the unitary step function. In this case all results are analytical. The general expression of $`\stackrel{~}{F}`$ is omitted here for brevity, but we give its linear limit obtained for small excitation probability $`\stackrel{~}{I}_0`$ and velocity $`\stackrel{~}{V}_B`$. The stalling force is

$$\stackrel{~}{F}(0)=\frac{2}{\stackrel{~}{D}_1}\stackrel{~}{J}_t+2\left(\frac{1}{\stackrel{~}{D}_1}1\right)\stackrel{~}{\overline{J}}_2,$$
(13)
with

$`\stackrel{~}{J}_t`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{I}_0\stackrel{~}{u}_d\stackrel{~}{u}_1^2e^{\stackrel{~}{u}_1}}{4\mathrm{\Delta }^2\left(1e^{\stackrel{~}{u}_1}\right)^2}}\left[\mathrm{\Delta }+{\displaystyle \frac{\mathrm{cosh}(\stackrel{~}{u}_1\stackrel{~}{u}_2/2)\mathrm{cosh}\delta }{\mathrm{sinh}\delta /(2\delta )}}\right]`$
$`\stackrel{~}{\overline{J}}_2`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{I}_0\stackrel{~}{u}_1\stackrel{~}{u}_2\stackrel{~}{\tau }}{4\mathrm{\Delta }}}\left[{\displaystyle \frac{1}{\stackrel{~}{u}_1\stackrel{~}{\tau }}}+{\displaystyle \frac{\stackrel{~}{u}_2}{2}}{\displaystyle \frac{\mathrm{cosh}\delta \mathrm{cosh}\frac{\stackrel{~}{u}_1}{2}\mathrm{cosh}\frac{\stackrel{~}{u}_d}{2}}{\mathrm{sinh}\frac{\stackrel{~}{u}_1}{2}\mathrm{sinh}\delta /\delta }}\right]`$
where $`\stackrel{~}{u}_d=\stackrel{~}{u}_2\stackrel{~}{u}_1`$, $`\delta =\sqrt{\stackrel{~}{u}_2^2/4+1/\stackrel{~}{\tau }}`$, and $`\mathrm{\Delta }=\stackrel{~}{u}_d\stackrel{~}{u}_1+1/\stackrel{~}{\tau }`$.

For $`U_1=U_2`$ and $`D_1D_2`$ (friction fluctuations only, no potential fluctuations), $`\stackrel{~}{J}_t=0`$ and the stalling force is proportional to $`\stackrel{~}{D}_11=(D_1D_2)/D_2`$. Therefore, everything else being fixed, the direction and strength of the force and of the ensuing motion is determined by the variation of diffusion constant (or friction coefficient) induced by the state fluctuations.

The nonlinear behaviour of the stalling force $`\stackrel{~}{F}(0)`$ versus pumping intensity $`\stackrel{~}{I}_0`$ and lifetime $`\stackrel{~}{\tau }`$ is illustrated in figure 4. The graph shows that there is an optimal $`\stackrel{~}{\tau }`$, i.e., an optimal ratio between the excited state lifetime and diffusion time. This optimal value decreases for increasing $`\stackrel{~}{I}_0`$. This behaviour is similar to that exhibited by the fluctuating-potential models .

For nonzero velocity $`\stackrel{~}{V}_B`$, the (linearized) force is

$$\stackrel{~}{F}(\stackrel{~}{V}_B)=\stackrel{~}{F}(0)\stackrel{~}{\eta }_{\mathrm{eff}}\stackrel{~}{V}_B,$$
(14)
where

$$\stackrel{~}{\eta }_{\mathrm{eff}}=\frac{1}{\stackrel{~}{D}_1}\left[1\frac{\stackrel{~}{u}_1^2e^{\stackrel{~}{u}_1}}{(1e^{\stackrel{~}{u}_1})^2}\right].$$
(15)
Therefore the maximum velocity of the motor is

$$\stackrel{~}{V}_{\mathrm{max}}=\frac{\stackrel{~}{F}(0)}{(\stackrel{~}{\eta }_{\mathrm{eff}}+\stackrel{~}{\eta }_p)}.$$
(16)
For $`U_1U_2`$ and $`D_1=D_2`$ (potential fluctuations only, no friction fluctuations), a nonzero current $`\stackrel{~}{J}_t`$ is established, i.e., a long-range displacement of body A with respect to body B occurs. If body B is also fixed (i.e., it is linked to body C), this case reduces to the “standard” one discussed for example in reference . The average speed of body A will then be $`\stackrel{~}{v}=2\stackrel{~}{J}_t=\stackrel{~}{D}_1\stackrel{~}{F}(0)`$.

An interesting difference between the two cases of fluctuating potential and fluctuating friction is given by the behaviour of $`\stackrel{~}{F}(0)`$ for large potential energies, i.e., $`\stackrel{~}{u}_i=u_i/kT1`$. Indeed, when $`D_1=D_2`$ the force shows a thermal-activated behaviour, with $`\stackrel{~}{F}e^{\stackrel{~}{u}_m}`$, where $`\stackrel{~}{u}_m`$ is the minimum of $`\stackrel{~}{u}_1`$ and $`\stackrel{~}{u}_2`$. On the contrary, when $`D_1D_2`$ and $`u_1=u_2`$ a term survives that decays only as $`\stackrel{~}{F}1/\stackrel{~}{u}_1`$. This behaviour is evident in figure 5. It can be explained by considering that for $`D_1=D_2`$ the force is proportional to the total average current $`\stackrel{~}{J}_t`$. To contribute to $`\stackrel{~}{J}_t`$ the motor moving element must overcome the potential barrier by acquiring sufficient thermal energy and the probability for this to happen is proportional to $`e^{\stackrel{~}{u}_i}`$. In the case of $`D_1D_2`$, the force acquires another term proportional to the average current $`\stackrel{~}{\overline{J}}_2`$ that need not vanish for $`\stackrel{~}{J}_t=0`$. It corresponds to a circulation of the particle that moves preferentially forward when in state 1 and backward when in state 2. However, the motor does not have to reach the potential maximum and therefore to be thermally activated. This cycling between the two states and the resulting currents are schematically described in figure 3.

## 4 Biological applications: actomyosin motor

We did not develop this model having in mind a specific biological application. Nonetheless, it is possible that basic features of some biological molecular mechanism be captured by our model. The basic working mechanism of important biological molecular motors is still under debate , so that there is room for proposing new models and offering different perspectives of this problem. Therefore we show in the following how our model can be applied to describing a biological molecular motor. However, we stress that we are not attempting here to build a “realistic” model, our purpose being only to illustrate the basic idea.

To be specific, let us consider the actin-myosin II system, that is the molecular motor providing the chemomechanical coupling in muscle contraction. When a muscle is contracted, thin filaments of actin and thick filaments of myosin are displaced relative to each other by the concerted action of many myosin “heads”. A myosin head is a molecular unit linked to the myosin filament through a flexible molecular “arm” (or “neck”) and interacting via weak intermolecular interactions with the actin filament . A schematic picture of the actomyosin system is in figure 6.

In reference  the actomyosin system is modeled by associating the myosin head with the particle (body A), the actin filament with the source of potential (body B), and the surrounding fluid or the actin filament itself with the friction source (body C). The myosin filament is treated as an additional passive “load” dragged on by the moving myosin head. The potential $`U(x)`$ describing the adsorption interaction of the head with the actin filament is periodic and asymmetric reflecting the underlying periodic and asymmetric structure of the filament. Attaching-detaching of the head to/from the actin filament are described as sudden changes in the potential landscape.

Here we propose a different set of correspondence rules in order to apply our model to the actomyosin system, as shown in figure 6. We still identify the Brownian particle (body A) with the myosin head. However, we model here the weak interaction between A and the actin filament as molecular friction, thus identifying the actin filament with body C. In a limit case, the low friction state allows almost free sliding of A on C, while the high friction state impedes any sliding. The myosin filament plays instead the role of body B. The potential $`U(x)`$ is used here to model the deformation energy of the molecular arm connecting the myosin filament with the head. The function $`U(x)`$ therefore needs not be periodic in this case, as the limited range of arm deformation constrains the particle A into a single potential well. The periodic potential formalism can be still used, however, by setting infinitely high energy barriers at $`\pm L`$.

Adenosinetriphosphate (ATP) binding on the myosin head and its subsequent hydrolysis triggers the variation of molecular friction coefficient of the myosin head (from high to low friction). By means of some suitable mechanism (specific intermolecular interactions), only when the conformation of the myosin arm is bent in a given direction, for example negative $`x`$ (figure 6b), the molecule enzymatic activity is large and ATP binding and hydrolysis may occur. This corresponds roughly to having $`I(x)>0`$ for $`x<0`$ (figure 6b) and $`I(x)=0`$ for $`x>0`$ (figure 6a), where the rate $`I`$ is proportional to ATP concentration. After an average time $`\tau `$, the hydrolysis reaction has gone through its full cycle and the motor switches back to a high friction state.

Given these correspondence rules, our model predicts a continuous directional motion of the myosin motor (head plus filament, i.e., bodies A and B) with respect to the actin filament (body C) (of course this is a relative motion: actually none of the three bodies can be considered fixed). The general behaviour of the motor as a function of external load and ATP concentration is in qualitative agreement with experimental data .

We emphasize that no direct causal link between the ATP hydrolysis and conformational transformations is introduced in our model, i.e., the conformational changes, described in our model by the continuous variable $`x`$, are not chemically-driven but are governed solely by Brownian motion. The asymmetric transition rates allow the elastic energy stored during the Brownian motion to be converted into directed long-range motion. In particular the main “force-generating step” occurs when the system happens to switch back to the high friction state for $`x>0`$ (figure 6a), relaxing then its deformation energy $`U`$ by moving back to $`x0`$. This step may be considered a sort of “Brownian power-stroke”. In contrast, our model assumes no “chemically-driven power-stroke”. In this sense, it differs from many other models of biological motors. A power-stroke is often modeled either as a sudden change of the elastic constant characterizing the potential , occurring as a consequence of ATP-induced state transitions, or as free sliding on an asymmetric potential as in Brownian ratchet models . On the other hand, the force-generation mechanism of our model in the limit of very high friction in state 2 and very low friction in state 1 becomes similar to that of the old model by A. F. Huxley . A similar sort of “Brownian power-stroke” is assumed also in references , in connection with the modeling of kinesin.

## 5 Conclusions

The model presented in this paper was actually developed for a specific application, namely the recently discovered strong enhancement of the light-induced molecular reorientation taking place in dye-doped nematic liquid crystals . A more detailed treatment of our modeling of this phenomenon has been already published elsewhere , and here we limit ourselves to a brief discussion. The relationship between this optical effect and molecular motors was first pointed out by Palffy-Muhoray and Weinan, without introducing however the possibility of an active role of fluctuating friction . In our model, dye molecules play the role of motor particles (body A) and the nematic molecular director plays the role of potential source (body B). The coordinate $`x`$ is here the angle between the dye molecule orientation and the molecular director, so that the motor is a rotary one. The role of the friction source C is played by the translational (centre-of-mass) degrees of freedom of liquid crystal molecules themselves. Light powers the motor by continuously promoting internal electronic transitions in the dye molecules, and the final output is a torque (corresponding to the force $`F`$) acting on the molecular director and eventually causing its reorientation. There is substantial experimental evidence  that the rotational friction coefficients of dye molecules in the two electronic states involved in the transitions are quite different, so that the fluctuating-friction mechanism is indeed at work. This proves that this hypothetical mechanism is actually realistic, at least in one specific example.

In this paper, our main goal was to illustrate in a general abstract form, independent of any specific application, the idea that a stochastic fluctuation of a kinetic coefficient, such as friction, mobility, or diffusion constant, can be an effective mechanism for converting chemical (or light) energy into directional motion and mechanical work at the molecular scale, at which Brownian fluctuations dominate. The possibility of extracting work from Brownian motion by means of a suitable modulation of kinetic properties, as opposed to equilibrium potential forces, appears to be new, at least within the field of molecular motors.

We believe that this concept can be profitable in several other fields, besides the nonlinear optics one. In particular, as discussed in section IV, many biological motors are still waiting for a detailed modeling, and any new idea can be very useful in the quest for complete understanding of these complex systems. Similarly, in the field of nanotechnology, the possibility of driving transport via a suitable modulation of the microscopic kinetic coefficients can be an interesting design concept.

In all cases, the relevance of this idea is further enhanced if one takes into account the high structural sensitivity of kinetic coefficients occurring commonly in activated systems, where small changes in the activation energy can lead to huge variations of the kinetic coefficients.

We close by adding that the fluctuating-friction motor idea is not limited to stochastic fluctuations. Just as in the case of fluctuating-potential motors, a periodic modulation in time of the friction coefficient is also expected to induce similar phenomena. The modulation period will then have to be close to the typical diffusion time of the particle, in order to be effective, leading to a new kind of stochastic resonance phenomenon.

This work was realized within the framework of the European Brite-Euram Network LC-PHOTONET.

## References


text:
# Fractal analysis of Sampled Profiles: Systematic Study

## I Introduction

The characterization of interfaces and of the mechanisms underlying their formation and evolution is a subject of paramount importance for a broad variety of phenomena such as crystal growth, rock fracture, biological growth, vapor deposition, surface erosion by ion sputtering, cluster assembling, etc … (Barabasi and Stanley (1995); Mandelbrot (1982); Family and Landau (1984); Meakin (1998); Pietronero and Tosatti (1985) and references therein). Since the pioneering work of B.B. Mandelbrot, fractal geometry has been widely used as a model to describe these physical systems that are too disordered to be studied with other mathematical tools but that still hold a sort of “order” in a scale-invariance sense Barabasi and Stanley (1995); Mandelbrot (1982); Vicsek (1992). In particular, the growth of interfaces resulting from the irreversible addition of subunits from outside (vapor deposition of thin films, low energy cluster beam deposition, etc…) shows a typical asymmetric scale invariance, because of the existence of a privileged direction (e.g. the direction of growth) Kardar et al. (1986); Witten and Sander (1983); Family and Vicsek (1985); Messier and Yehoda (1985); Meakin et al. (1986); Meakin (1987, 1998); Lichter and Chen (1986); Family (1990); Krim and Palasantzas (1995); Csahok and Vicsek (1992); Bales et al. (1990); Muller-Pfeiffer et al. (1990); Sit et al. (1999); Yoona et al. (1999); Milani et al. (2001a, b); Buzio et al. (2000). These interfaces belong to the class of self-affine fractals and they can be described either by the fractal dimension $`D`$ or by the well-known Hurst exponent $`H`$ Peitgen,H.-O. and Saupe (1989); Meakin (1989); Voss (1989); Schmittbuhl et al. (1995a); Mandelbrot (1985, 1986). If these systems are the result of a temporally evolving process, they usually show also a time scale-invariance described by the exponent $`\beta `$ Vicsek (1992); Barabasi and Stanley (1995). Because of the close relationship between the scaling exponent(s) and the fundamental mechanisms leading to scale invariance, universality classes can be defined Vicsek (1992); Barabasi and Stanley (1995). An accurate knowledge of $`H`$ (and $`\beta `$) is required to identify the universality class of the system and to give a deep insight on the underlying formation processes.

The possibility of characterizing the topography of an interface in a dimension range from the nanometer up to several tens of microns, in a relatively simple and quick way by Atomic Force Microscopy (AFM) and Scanning Tunneling Microscopy (STM) Binnig et al. (1982, 1986) has stimulated an upsurge of experimental report claiming for self-affine structures (see Refs. Malcai et al. (1997); Avnir et al. (1998) and references therein). The abundance of experimental characterization of different systems and the limited sampling capability of the scanning probe microscopies (SPM) prompted at the attention of many authors the need of an accurate methodological approach to the determination of the exponent $`H`$ and of its error Deng et al. (1999a); Tate (1998), realistically considering the consequences of the finite sampling inherent to SPM. Typical sampling with an AFM or a STM is 256 or 512 points per line, for a maximum of 512 lines. Most of the results published in the late eighties and early nineties were based upon 256x256-point data-sheets, or even smaller ones (see list of references in Ref. Malcai et al. (1997)). Commercially available SPMs offer today a maximum of 512x512-point resolution, and homemade instruments hardly go beyond this value.

Many authors have questioned the reliability of the measurement of the Hurst exponent from a poorly sampled profile Dubuc et al. (1989); Schmittbuhl et al. (1995b); Mehrabi et al. (1997); Schmittbuhl et al. (1996). In order to quantify the influence of the sampling on the determination of $`H`$, a numerical analysis can be performed on artificial self-affine profiles, generated with a specific algorithm, with a fixed number of points $`L`$ and known Hurst exponent $`H_{in}`$. The “true” exponents ($`H_{in}`$) are then compared with the ones measured directly from the generated profiles ($`H_{out}`$). Usually a sensible discrepancy between the measured $`H_{out}`$ and the expected $`H_{in}`$ is found Dubuc et al. (1989); Mehrabi et al. (1997); Schmittbuhl et al. (1996). The discrepancy is not uniform but depends on the value of $`H_{in}`$. As one would expect, the discrepancy is globally dependent on the number $`L`$ and it approaches zero for large values of $`L`$. In particular, for $`L<1000`$ the sampling effect is of great importance since the discrepancy can be of the order of the exponent itself (100% relative error) Schmittbuhl et al. (1995b). Dubuc et al. have reported that even for values of $`L`$ as high as 16384, the discrepancy is still significant Dubuc et al. (1989).

Although the problem of sampling has been clearly addressed and discussed, quite surprisingly a systematic analysis of the problem, considering different generation algorithms, is still lacking. The dependence of the sampling effect on $`L`$ has been investigated Dubuc et al. (1989); Schmittbuhl et al. (1995b) and also many different methods for the measurement of $`H_{out}`$ have been considered for different values of $`H_{in}`$ in the range $`[0.11]`$ Dubuc et al. (1989); Schmittbuhl et al. (1995b); Mehrabi et al. (1997); Schmittbuhl et al. (1996). However, either only one single generation algorithm has been used Schmittbuhl et al. (1995b); Schmittbuhl et al. (1996), or the results from different generation algorithms have not been compared Mehrabi et al. (1997). We believe that this comparison is of fundamental importance.

Indeed profiles from different generation algorithms can be considered as different self-affine objects sampled in $`L`$ points. For a fixed value of $`H_{in}`$, these objects would all have the same fractal dimension if they were sampled with an infinite number of points. The fundamental question at this point is whether the discrepancy of $`H_{out}`$ from $`H_{in}`$, for a finite value of $`L`$, is the same for every self-affine object (i.e. for every generation algorithm). Only an analysis that considers different self-affine objects has a statistical validity and allows a reliable interpretation of the results. Up to now the results obtained in literature from a single generation algorithm did non allow a discussion of the nature of the aforementioned discrepancy, which has been interpreted as an uncontrollable error affecting the analysis of sampled profiles. The main conclusion drawn by these authors is the non-reliability of results obtained from profiles with less than 1024 sampling points Schmittbuhl et al. (1995b).

Our aim is to achieve a deeper understanding of the effects of sampling in order to answer the question whether the measurement of the Hurst exponent with a poor number of sampling points is reliable or not. This point is crucial both for future analysis of self-affine profiles and for a correct interpretation of the results already present in literature.

From a more general point of view, fractality is characterized by the repetition of somehow similar structures at all length scales and can be described in its major properties by a single number: the fractal dimension $`D`$ Falconer (1990); Mandelbrot (1982). Any finite sampling of a fractal object poses both an upper and a lower cut-off to this scale invariance. It has been shown that these cut-offs introduce a deviation in $`D`$ and the sampled object has a dimension different from the one of the underlying continuous object Dubuc et al. (1989); Schmittbuhl et al. (1996); Mehrabi et al. (1997). However, it is still unknown whether the sampling influences in a different way different objects characterized by the same ideal dimension, thus breaking the sort of universality that makes a fractal be identified by its dimension only.

In this paper we present a systematic analysis considering together all the generation algorithms found in literature. The aim of our analysis is to understand whether the discrepancy of the measured $`H_{out}`$ for a fixed $`L`$ and for every generation algorithm is completely random or has a universal dependence on $`H_{in}`$. The latter observation can be interpreted as a reminiscence of the fact that a fractal object is completely characterized by its dimension 111As pointed out by Voss et al. in Refs. Voss (1986a, b) this is not completely true, because fractal objects have also properties that do not depend on the fractal dimension only, such as lacunarity. However, the main statistical properties of a fractal object are characterized by its dimension Mandelbrot (1982); Falconer (1990). The distinction is of crucial importance because in the case of universal dependence of $`H_{out}`$ on $`H_{in}`$ , one can empirically correct the discrepancy of the measured exponents from the “true” ones. Some authors independently suggested to use directly the $`H_{out}`$ vs. $`H_{in}`$ curves as correction, but they considered only one generation algorithm without discussing the universal character that these curves must have in order to be utilized for any self-affine object Deng et al. (1999a).

Conversely, on the basis of our analysis, we will interpret the discrepancy in terms of two distinct contributions: a universal deviation and a random dispersion. We will propose a powerful method to correct the universal deviation and we will discuss the nature of the dispersion, which is due to both statistical fluctuations and an intrinsic sampling effect. The latter turns out to be a sort of systematic error that cannot be corrected unless one knows the generation algorithm that produced the self-affine object. In the case of generic self-affine profiles which have not been generated by a specific algorithm, such as experimental profiles, the above arguments no longer hold. A new procedure to quantify the intrinsic error in the measurement of the Hurst exponent of generic self-affine profiles is thus needed.

On these basis, we will discuss the effect of sampling on the reliability of the fractal analysis of poorly sampled self-affine profiles, focusing on both the deviation and the dispersion of the measured exponents from the ideal ones, showing that the conclusions drawn by Schmittbuhl et al. that “… a system size less than 1024 can hardly be studied seriously, unless one has some independent way of assessing the self-affine character of the profiles and very large statistical sampling” were too restrictive Schmittbuhl et al. (1995b). Moreover, we will point out that the estimate of the intrinsic error is essential for a correct classification of a process in terms of universality classes. In fact, in order to distinguish exponents belonging to different classes, it is necessary to quantify the error on the measurement. Up to now, the statistical error or the error of the fit have been used to quantify the error on the measurement of $`H`$ Krim et al. (1993); Deng et al. (1999b); Iwasaki and Yoshinobu (1993). Both the statistical error and the error of the linear fit can be made very small, if a large number of profiles are averaged. However, if the measurement is likely to be affected by more subtle intrinsic errors, such as the aforementioned dispersion due to the sampling, considering only the statistical error may lead to serious misleading. The intrinsic error in many cases may indeed be much larger than the statistical one.

In the following sections we will present a systematic analysis of synthetic self-affine profiles with the aim of both achieving a deep understanding of the effects of sampling and providing the experimentalists of a reliable tool for the fractal analysis of surfaces and interfaces. To this purpose we have developed a new automated fitting protocol in order to avoid any arbitrariness in the measurement. With this new methodology we will study the effects of sampling, enlightening the main characteristics of the deviation and the dispersion of the measured exponents. We will present a new powerful method to correct the deviation of $`H_{out}`$ and to estimate the error of the measurement. Finally, we will apply our empirical correction procedure to 512-point profiles created with the directed percolation (DP) algorithm Buldyrev et al. (1992). This system provides a simple benchmark to test our protocol and allows noticing the opportunity of the correction.

## II The Automated Fitting Protocol

Self-affine systems occurring in nature are usually profiles or surfaces. In order to measure their Hurst exponents the 2+1 dimensional case of surfaces is usually reduced to 1+1 dimensions, considering the intersection of the surface with a normal plane. The particular case of in-plane anisotropy results in a dependence of $`H`$ on the orientation of the plane with respect to the surface Dubuc et al. (1989); Schmittbuhl et al. (1995b); Falconer (1990); Barabasi and Stanley (1995).

Once we have scaled down the analysis to 1+1 dimensions, the following general properties characterize a self-affine profile. If $`h(x)`$ is the height of the profile in the position $`x`$, the orthogonal anisotropy can be expressed by the scaling relationship:

$$h(\lambda x)=\lambda ^Hh(x)$$
(1)
where $`H(0,1)`$ is the Hurst exponent, $`\lambda `$ is a positive scaling factor and the equation holds in a statistical sense Barabasi and Stanley (1995); Stanley (1971). The fractal dimension $`D`$ of the profile is related to the Hurst exponent by the equation $`D=2H`$ while the dimension of the surface is $`D=3H`$ Mandelbrot (1986); Moreira et al. (1994). The lower is $`H`$, the more space invasive is the surface. In most of the physical self-affine surfaces, the scale invariance does not extend to all length scales but there is an upper cut-off above which the surface is no longer correlated. The length at which this cut-off appears is defined as the correlation length $`\xi `$ Barabasi and Stanley (1995); Malcai et al. (1997). In the present analysis, we consider only profiles whose correlation length (expressed in number of points) is equal to their length $`L`$. To this purpose we have carefully studied each generation algorithm in order to grant the condition $`\xi =L`$. For this reason we were often forced to generate very long profiles and to consider only their central portion Makse et al. (1996); Mehrabi et al. (1997); Simonsen and Hansen (1999). The usual procedure to measure the Hurst exponent of a self-affine profile $`h(x)`$ is to calculate appropriate statistical functions from the whole profile. These functions of analysis (AFs) show a typical power law behavior on self-affine profiles:

$$AF[h(),k]=ck^{f(H)}$$
(2)
where $`c`$ is a constant, $`k`$ is a variable indicating the resolution at which the profile $`h`$ is analyzed (typically a frequency or a spatial/temporal separation), and $`f(H)`$ is a simple function of the Hurst exponent $`H`$ Mehrabi et al. (1997); Barabasi and Stanley (1995); Yang et al. (1997); Moreira et al. (1994); Peng et al. (1994); Press et al. (1986); Simonsen et al. (1998). The power law behavior of the AF is then fitted in a log-log plot in order to calculate the exponent $`H`$. In the analysis of statistical self-affine profiles there are random fluctuations super-imposed to this power law behavior. The signal-to-noise ratio of these fluctuations is scale-dependent, the AFs being calculated as averages of statistical quantities at different length scales Barabasi and Stanley (1995). To reduce this noise, the average of the AFs obtained from $`N`$ independent profiles is usually taken before the execution of the linear fit. However, while small-scale fluctuations are easily smoothed, larger scale fluctuations converge very slowly.

The identification of the linear region in the analysis of the AFs is a puzzling point. Windowing saturation is present at length scales comparable with the profile length depending on the nature of the profiles Yang et al. (1997). This results in a departure from the power law behavior to a constant value. Moreover, the degradation of the fractality due to the sampling causes a diversion of the AFs from their ideal power law behavior. This produces both a discrepancy of the measured Hurst exponent from the ideal value (a change of the slope in the log-log plot) and a shortening of the linear region as shown in Fig. 1.

Here, the presence of curved regions is clearly visible. It can be seen that this anomalous behavior is not localized at length scales close to the length of the profile, but involves also the shortest length scales especially for values of $`H`$ close to zero. It is important to notice that this effect is not due to experimental conditions, such as the finite size of the SPM scanning probe. Thus it is necessary, in particular for small values of $`H`$, to chose a linear region instead of fitting the whole function. The methods proposed in the literature to identify the linear region (e.g. the consecutive slopes method Bourbonnais et al. (1991); Barabasi and Stanley (1995), correlation index method Barlow (1989), the coefficient of determination method Hamburger et al. (1996) and the “fractal measure” method Li et al. (1999)) are usually based on an arbitrary (human) choice. This is particularly delicate since the curvature in the AFs can be so small, if compared to the statistical noise, that it is hard to distinguish the correct linear region. Because of this reason, we think that the proposed methods suffer of a high degree of arbitrariness. Moreover, all these methods make no distinction between a straight line with statistical noise and a slightly curved line.

Due to the previous arguments and since no universally accepted fitting procedure is available in literature, we were prompted to develop an automated fitting protocol (AFP) with two purposes: to reduce as much as possible the effects of the curved regions on the measured exponent, and to define a standard algorithm for the choice of the linear region, eliminating, as much as possible, any arbitrariness. This is very important for the reliability of the results, in particular for the comparison of different generation algorithms. Moreover, the automation of the fitting procedure is essential to perform a systematic analysis. In fact, in order to have good statistics, a large number of AFs must be calculated and fitted.

In our procedure, that is an implementation of the consecutive slopes algorithm Barabasi and Stanley (1995), the curve to be fitted is divided in many portions of the same length $`\mathrm{}`$ (in number of points) and each of them is considered separately. A linear and a cubic fit are performed on each portion. Comparing the mean distance of the linear fit from the portion to the mean distance of the cubic from the linear fit, we evaluate whether the portion is almost linear with uncorrelated noise or it presents a definite curvature. Obviously, the distinction is not immediate and we have to set a threshold to separate the two cases through a parameter in the fitting procedure. The use of a parameter is common to other methods (see for example the coefficient of determination method used in Ref. Hamburger et al. (1996)). Once the fitting parameter is set, our procedure is able to decide automatically whether the portion is “curved” or “linear”. Only the “linear” portions are then considered. They undergo a straight-line-fit analysis through which the slopes and their errors are determined. A distribution of the slopes weighted with the values of the errors is then built (see Fig. 2a) and its main peak position and width are measured. We do not consider here the presence of more than one linear region with different slopes. Thus, there is a well-defined main peak in the distribution. We have extended our procedure also to the case of more than one linear region, but this extension is out of the scopes of this article.

The procedure described above is repeated varying the length $`\mathrm{}`$ of the portions from a minimum value $`\mathrm{}_{min}`$ up to the length of the curve. The results are then shown in a plot of the peak position (i.e. a slope value) versus the length of the portion, with the peak widths as error bars (see Fig. 2b).

If the analyzed curve presents a linear region, this plot shows a plateau for $`\mathrm{}`$ ranging from $`\mathrm{}_{min}`$ to the length of the whole linear region. This plateau is usually very easy to be identified because of the distinction between linear and curved portions. In fact, portions of length larger than the length of the whole linear region are considered curved portions and discarded. Thus, the plot usually drops to zero at the end of the plateau. Eventually, through an average and a standard deviation, we obtain the final slope value and its fitting error, while the length of the plateau gives the length of the linear region. In conclusion, our AFP is able to identify not only the slope of the linear region but also its length. We have tested our AFP before its application to the systematic analysis and we have found that the measured Hurst exponent is widely independent of the fitting parameter 222As Hamburger et al. did in Ref. Hamburger et al. (1996), we have not chosen a value of the parameter, being this choice not relevant to the purposes of the article. Nevertheless, we have tested a wide range of reasonable values of the parameter in order to be sure of the repeatability of the results obtained.. Conversely, the length of the linear region strongly depends upon the value of the parameter and must be considered only an internal parameter of the analysis and not a direct measurement of the scale invariance range.

## III Numerical Analysis

With all the generation algorithms published in literature we have created sampled self-affine profiles with known fractal dimension $`D=2H`$. We have varied the exponent $`H`$ between 0.1 and 1 and we have focused on the value $`L=512`$ sampling points (the best sampling obtainable with most of the SPMs). We discuss also different values of $`L`$ up to 16384. Because there exists only a few algorithms that generate exactly self-affine profiles, we have used algorithms that generate statistically self-affine profiles, which are more difficult to handle but closer to reproduce natural physical systems. The algorithms we have used are known in literature as: the random midpoint displacement Schmittbuhl et al. (1995b); Voss (1986a), the random addition algorithm Feder (1988); Peitgen,H.-O. and Saupe (1989), the fractional Brownian motion Feder (1988), the Weierstrass-Mandelbrot function Lopez et al. (1994); Berry and Lewis (1980), the inverse Fourier transform method Voss (1986a) and a variation of the independent cut method Falconer (1990). For the measurement of the Hurst exponent of self-affine profiles we have used the height-height correlation function $`C_2`$ Yang et al. (1997) and the root mean square variable bandwidth with fit subtraction method Moreira et al. (1994); Peng et al. (1994). The value of $`H_{out}`$ has been calculated from the slope in the log-log plot of the average over $`N`$ statistically independent AFs, measured with our AFP.

The results are expressed in terms of $`H_{out}`$ vs. $`H_{in}`$ plots. Each plot is characteristic of a single AF and generation algorithm and it represents the relationship between the measured Hurst exponent $`H_{out}`$, calculated from the average of $`N`$ AFs, and the nominal exponent $`H_{in}`$ of the profile. Grouping the $`H_{out}`$ vs. $`H_{in}`$ plots obtained using the same AF for all the generation algorithms, the dispersion of the $`H_{out}`$ values comes to evidence.

In Fig. 3 we show the $`H_{out}`$ vs. $`H_{in}`$ graphs obtained from $`N=500`$, $`L=512`$ profiles, as explained in the previous section. We show separately in Figs. 3(a) and 3(b) the different AFs used. Since the profiles are statistically self-affine, the measured $`H_{out}`$ are subject to a statistical error that is inversely related to $`N`$ Deng et al. (1999b). In order to characterize the dependence of this statistical error on the number $`N`$ of averaged AFs, we let $`N`$ vary from 1 to 50 using the same profiles considered in Fig. 3. With these values of $`N`$ we have repeated the numerical analysis (i.e. calculation of the AFs, averaging and application of the AFP) and we have extracted a standard deviation $`\sigma _N`$ of the measured exponents.

In Fig. 4 we show the $`H_{out}`$ vs. $`H_{in}`$ graphs, analogous to those in Fig. 3, with the calculated error bars (twice the standard deviation $`\sigma _N`$), for a few values of $`N`$. We present the results for a single AF (the root mean square variable bandwidth with fit subtraction), the results for the other AFs being similar.

In Fig. 5 we show three $`H_{out}`$ vs. $`H_{in}`$ graphs obtained respectively with $`N=500`$, $`L=512`$ profiles, $`N=50`$, $`L=4096`$ profiles and $`N=15`$, $`L=16384`$ profiles. Again, we present only one AF (the height-height correlation function $`C_2`$).

## IV Results and Discussion: Deviation and Dispersion from the Ideal Behavior

Ideal continuous fractal profiles are statistically characterized by their fractal dimension (universality) and their $`H_{out}`$ vs. $`H_{in}`$ graphs are straight lines Barabasi and Stanley (1995); Feder (1988); Falconer (1990).

In Fig. 3 a deviation from the ideal behavior is observed for both the AFs. It turns out that the sampling of a profile affects in a different way different methods of analysis. The deviation from the ideal behavior has been already observed in literature (for example, see Ref. Schmittbuhl et al. (1995b)) and our results are in good agreement with the previous ones.

Moreover, within the same method of analysis we observe that the different generation algorithms give significantly different $`H_{out}`$ vs. $`H_{in}`$ plots. This dispersion is pointed out here for the first time because different generation algorithms are considered together. The significance of the dispersion can be inferred from the characterization of the statistical error of the measured exponent discussed hereafter.

In Fig. 4 we show that for $`N>25`$ and $`H_{out}<0.3`$ the error bars of $`H_{out}`$ for different generation algorithms hardly overlap. This fact suggests that the statistical error is not the only reason of the differences between the $`H_{out}`$ vs. $`H_{in}`$ plots shown in Fig. 3.

In Fig. 6 we plot the statistical error $`\sigma _N`$ times the square root of $`N`$ vs. $`N`$. For $`N10`$ the curves approach a constant value according to the relationship between the standard deviation of independent, normally distributed measurements and the standard deviation of the mean upon $`N`$ measurements:

$$\sigma _N=\frac{\sigma }{\sqrt{N}}$$
(3)
This result shows that the AFP and the averaging of the AFs do commute. The assessment of this property is non-trivial due to the complexity of the AFP. Thus, we extrapolate the statistical error of the measured exponents in Fig. 3 ($`N=500`$) using Eq. (3) where $`\sigma `$ is extracted from the plateau in Fig. 6. Overestimating $`\sigma `$ with the value 0.16 we obtain $`\sigma _{500}=0.007`$. This value produces an error bar in Fig. 3 as small as the symbol used to mark the data. A direct calculation of $`\sigma _{500}`$, obtained averaging AFs calculated on groups of $`N=500`$ profiles for every $`H[0.1,1]`$ and for every generation algorithm, fitting and extracting a mean value and a standard deviation of $`H`$, would have required a huge and time consuming calculation.

These results suggest that the observed dispersion between the $`H_{out}`$ vs. $`H_{in}`$ curves for different generation algorithms is an intrinsic effect of the sampling, depending only on the number of sampling points $`L`$. This fact has an important consequence on a fractal analysis of experimental surfaces. While looking at a real sample, we do not know what kind of “algorithm” has generated the surface. This introduces an uncertainty on its real fractal dimension independent of the statistical error. Thus, there is an intrinsic upper limit to the precision of the measurement of the exponent. It is useless to strengthen the statistics once the number of acquired profiles makes the statistical error smaller than the intrinsic dispersion.

In Fig. 5 we see that as $`L`$ increases both the deviation and the dispersion decrease in agreement with their expected vanishing in the limit of $`L`$ going to infinity Schmittbuhl et al. (1995b). This is also an a posteriori proof of the correctness of both the generation algorithms and the methods of analysis.

Our interpretation of these effects is that the sampling of a self-affine profile lessens its fractality in such a way that it is no longer characterized universally by its fractal dimension (or Hurst exponent). While for a continuous self-affine profile the relationship $`H_{out}=H_{in}`$ holds, for sampled profiles we can see that different AFs produce different $`H_{out}`$ vs. $`H_{in}`$ plots from the same sampled fractal profile. Considering instead a single AF, our results show that sampled fractal profiles generated with different generation algorithms but with the same ideal dimension give different measured Hurst exponents.

However, Fig. 5 clearly shows that the lessening of fractality of a profile is rather a continuous process than a sharp transition: the poorer is the sampling, the worse are the deviation and the dispersion. In Figs. 5 and 3 we observe that the lessening of fractality acts in a similar way on profiles generated with different algorithms. The common trend of the $`H_{out}`$ vs. $`H_{in}`$ curves obtained from different generation algorithms is interpreted as a consequence of the universality of fractal objects.

It is then reasonable to assume the existence for every AF of a universal region in the $`H_{out}`$-$`H_{in}`$ plane containing all the $`H_{out}`$ vs. $`H_{in}`$ plots obtained with every possible generation algorithm. This region, approximately identifiable with the envelope of the $`H_{out}`$ vs. $`H_{in}`$ plots, has a width that depends on the number of sampling points and approaches the 1-dimensional $`H_{out}=H_{in}`$ ideal curve for very large values of $`L`$. We expect that, given any continuous self-affine profile with a Hurst exponent $`H_{in}`$ and given the exponent $`H_{out}`$ measured from an $`L`$-point sampling of the continuous profile, the pair ($`H_{in}`$,$`H_{out}`$) belongs to the universal region of the corresponding graph (specific for every AF and number of sampling points $`L`$). Provided a good characterization of the aforementioned regions (i.e. using as many generation algorithms as possible), we can use them to generate calibration graphs for every $`L`$ and AF describing the relationship between the measured $`H_{out}`$ and the true value $`H_{in}`$.

To produce the calibration graphs we proceed as follows. First of all, we make two general assumptions in order to take quantitatively into account the problem of measuring the Hurst exponent of a sampled profile. We assume that the $`H_{out}`$ values corresponding to the same $`H_{in}`$ are normally distributed around a mean $`H_{out}`$, and we assume also that the values obtained with the available generation algorithms are a random sampling of the gaussian distribution. We then measure the average and the standard deviation of the dispersed $`H_{out}`$ values corresponding to each $`H_{in}`$ separately. Thus we obtain a sampling of the functions describing the dependence of $`H_{out}`$ and $`\sigma _{H_{out}}`$ from $`H_{in}`$. With an interpolation algorithm using smooth functions, we derive the curve representing the relationship between $`H_{out}`$ and $`H_{in}`$. We also derive the pair of curves corresponding to $`H_{out}+n\sigma _{H_{out}}`$ and $`H_{out}n\sigma _{H_{out}}`$ vs. $`H_{out}`$ which define the $`n`$-th confidence level. For every value of $`H_{out}`$ it is possible to find the confidence interval of $`H_{in}`$ for any given confidence level. The resulting graphs for $`L=512`$ are shown in Fig. 7.

These calibration graphs allow to take into account the deviation and the dispersion due to the sampling. A similar method has been independently proposed in Ref. Deng et al. (1999a) even though the analysis was limited to a single generation algorithm and the discussion on the reliability of the calibration regions together with the intrinsic dispersion were completely neglected.

Using the calibration graphs it is possible to measure the Hurst exponent of poorly sampled profiles correcting for the first time the deviation due to the sampling and providing a reasonable estimate of the error on a confidence level basis. The quantification of the error is of paramount importance, as pointed out in the introduction, since many authors estimated the error from the precision of the linear fit Krim et al. (1993); Iwasaki and Yoshinobu (1993) or from the standard deviation of the measured exponents Deng et al. (1999b). Our results show that they usually underestimated the true error.

## V Application of the calibration graphs to the study of directed percolation numerical profiles

We have applied our procedure to the 1+1 dimensional directed percolation (DP) model, described by S.V. Buldyrev et al. Buldyrev et al. (1992). This model mimics the paper wetting process by a fluid. The resulting pinned interface is self-affine with exponent $`H0.63`$.

We have analyzed $`N=30`$, $`L=16384`$ DP profiles with the height-height correlation function (h.-h. corr) and the variable bandwidth with fit subtraction (vbw), using the automated fitting protocol to measure the Hurst exponents. The results are shown in the second column of Tab. 1. We have not calculated the statistical error (see Section IV) because it would have been excessively time consuming. Thus, the error shown is simply the error of the fit calculated with the AFP.

The values of the measured exponents $`H_{out}^{16384}`$ are significantly lower than the ones predicted by the DP model, suggesting that a correction is needed even in the case of profiles of $`L=16384`$ points, which are widely considered as continuous.

We have then analyzed $`N=1000`$, $`L=512`$ profiles extracted from the $`L=16384`$ profiles. We have applied the correction procedure based on the calibration graphs shown in Fig. 7 to the exponents measured with the AFP. In the third column of Tab. 1, the uncorrected measured exponents ($`H_{out}^{512}`$) are shown. The error is calculated as the root mean square (rms) value of the statistical error $`\sigma _{1000}`$ (evaluated as explained in Section IV) and the error of the fit calculated with the AFP. In the fourth column, the confidence intervals corresponding to the $`68\%`$ probability for the “true” exponents are shown ($`H_{in}^{512}(68\%)`$).

The results summarized in Tab. 1 allow to notice the effectiveness of the calibration graphs in the analysis of self-affine profiles when the effects of sampling are non negligible. In the example reported here, the poor sampling causes a discrepancy of about $`4\%`$ between the measured exponents and the theoretical one for DP profiles. After the correction with the calibration graphs, the expected value $`H_{in}0.63`$ is consistent with the confidence intervals of the three AFs. Moreover, the intrinsic error due to the dispersion (about half the width of the confidence interval) turns out to be usually one order of magnitude larger than the aforementioned rms error.

In conclusion, our calibration graphs have allowed to correct the deviation and to quantify the intrinsic error of the Hurst exponent of poorly sampled ($`L=512`$) DP profiles.

## VI Conclusions

We have carried out a systematic analysis in order to achieve a deeper understanding of the effects of sampling on the measurement of the Hurst exponent of self-affine profiles. This is a crucial point for the assessment of the reliability of fractal analysis of experimental profiles, such as topographic profiles of growing thin films and interfaces acquired with a Scanning Probe Microscope. We have pointed out that some of the steps leading to the measurement of the Hurst exponent have been only superficially discussed, although worth of deeper attention. We have focused on the quantification of the effects of sampling and possibly on their correction, allowing a more reliable identification of the universality class of growth.

In order to perform such a quantitative analysis we have developed a new automated fitting protocol that allows to remove the ambiguity in the choice of the region for the linear fit of the analysis functions. This point is usually underestimated in the published experimental literature, and appears to be a significant source of error in the whole analysis. Moreover, an automated protocol sensibly reduces the time required for the fitting of a large number of noisy curves, allowing a higher statistics. With our automated fitting protocol we have systematically investigated synthetic self-affine profiles generated with all the generation algorithms found in literature using different method of analysis.

The systematic analysis presented in this paper has been carried out on 1+1 dimensional profiles and we have not considered 2-dimensional methods of analysis (e.g. see Deng et al. (1999a); Krim et al. (1993)). However, it is reasonable to suppose that even in this case the effects of sampling cannot be neglected, and the conclusions drawn in Ref. Deng et al. (1999a) are probably incorrect. The similarity between Fig. 1 in Ref. Deng et al. (1999a) and the analogous results presented in this paper (see the variable bandwidth analysis of profiles generated with the random midpoint displacement shown in Fig. 3c) suggests that conclusions very close to those presented here can be drawn also in the 2-dimensional case.

Studying the discrepancy between the measured Hurst exponent $`H_{out}`$ and the “true” one ($`H_{in}`$) for synthetic self-affine profiles with $`L=512`$ points, we have shown that the main effects of sampling are a deviation of the $`H_{out}`$ vs. $`H_{in}`$ plots from the ideal behavior and a dispersion of the exponents calculated from different generation algorithms. Both these effects smoothly reduce with increasing values of $`L`$. The deviation turns out to be universal in the sense that the trend of the $`H_{out}`$ vs. $`H_{in}`$ curves is common to all of the generation algorithms, depending only on the number of sampling points and on the function used in the analysis. We propose that this behavior is reminiscent of the fact that a fractal object is completely characterized by its dimension and therefore the deviation can be at least empirically corrected. The dispersion instead has to be considered as an intrinsic error due to the sampling, but for the very special case of profiles whose generation algorithm allows to build their specific $`H_{out}`$ vs. $`H_{in}`$ plot. This dispersion error must be quantitatively taken into account since it cannot be reduced with an increase in the statistics but only with an increase in the number of sampling points.

The existence of an intrinsic dispersion error in the measurement of the Hurst exponent that depends only on the number of sampling points is very important. In fact, this intrinsic error easily overwhelms the statistical error for poorly sampled profiles. It is definitely clear that a reliable result cannot be based on the consideration of the statistical error only. Moreover, the dispersion poses an upper limit to the precision in the measurement of the Hurst exponent of sampled profiles. It becomes useless to increase the statistics once the statistical error has been made reasonably smaller than the intrinsic one. This is particularly important in an experimental analysis because it usually reduces significantly the number of profiles that have to be acquired, making the analysis much less time consuming.

Thanks to our systematic analysis, we have built, for each method of analysis, a calibration graph representing the region of the $`H_{out}`$-$`H_{in}`$ plane where the true exponents fall within a given confidence level. We have originally proposed to use these graphs as a reliable empirical method to correct the measured value of the Hurst exponent of a poorly sampled profile and to estimate its intrinsic sampling error. The reliability of the calibration graphs is based on two assumptions:

1. The measured exponents for all the possible self-affine profiles, with the same “true” exponent $`H_{in}`$ and with the same number of sampling points, are normally distributed;
2. The numerical generation algorithms known in literature provide a statistically reliable sample of all the possible self-affine profiles.

Even though we have found just six generation algorithms in literature, we believe that they still allow to obtain reasonable results or at least the only ones obtainable to date. These results represent a step forward to a reliable fractal analysis of both numerical and experimental profiles and to the individuation of the universality classes in the study of the evolution of many different systems.

In conclusion, we have demonstrated that a reliable measurement of the Hurst exponent of poorly sampled self-affine profiles is possible, provided that the measured $`H_{out}`$ is corrected of its deviation and that the sampling error is quantitatively taken into account. We have thus given strength to experimental analyses, since the numerical results reported in literature to date led to the conclusion that the analysis of self-affine profiles sampled with less than 1000 points is not reliable Schmittbuhl et al. (1995b). Even with the great improvement introduced by the use of the calibration graphs in the analysis of self-affine profiles, we definitely agree with Schmittbuhl et al. in pointing out that the comparison of the results obtained with different method of analysis is of fundamental importance Schmittbuhl et al. (1995b). Furthermore, we shortly comment on the common experimental procedure of connecting AFs calculated from profiles acquired with different scan sizes Krim et al. (1993); Iwasaki and Yoshinobu (1993); Fang et al. (1997). This connection allows investigating a wider range of length scales with a limited number of sampling points and makes the measurement more reliable. However, the deviation and dispersion are not influenced by this procedure, since they depend only on the number of sampling points of the profiles on which the AFs are calculated.

The AFP and the calibration graphs have been tested on numerically generated 1+1 dimensional directed percolation (DP) profiles, which have provided a benchmark to check our protocol. We have shown that for $`L=512`$ profiles a correction is needed and the calibration graphs allow to recover the theoretical value of $`H`$ predicted by the DP model. We have also shown that a correction is needed even for the $`L=16384`$ profiles, which are widely considered as continuous.

Our results provide a powerful tool for the accurate extraction of the Hurst exponent from poorly sampled profiles, and for the quantification of the error in the measurement. This is of paramount importance for experimentalists who study the scale invariance of surfaces and interfaces by Scanning Probe Microscopy or other techniques, with the aim of identifying the underlying universality classes. The huge amount of experimental results published in the past two decades about the fractality of many interfaces can be now analyzed under a new light.

###### Acknowledgements.

We thank E.H. Roman and G. Benedek for discussions. Financial support from MURST under the project COFIN99 is acknowledged.