Roughening of the anharmonic Larkin model
V. H. Purrello J. L. Iguain A. B. Kolton


vpurrello@ifimar-conicet.gob.ar

iguain@mdp.edu.ar
Instituto de Investigaciones Físicas de Mar del Plata (IFIMAR),
Facultad de Ciencias Exactas y Naturales,
Universidad Nacional de Mar del Plata,
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET),
Deán Funes 3350, B7602AYL Mar del Plata, Argentina


koltona@cab.cnea.gov.ar
Centro Atómico Bariloche and Instituto Balseiro,
Comisión Nacional de Energía Atómica (CNEA),
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET),
Universidad Nacional de Cuyo (UNCUYO),
Av. E. Bustillo 9500, R8402AGP San Carlos de Bariloche, Río Negro, Argentina




We study the roughening of d-dimensional directed elastic interfaces subject
to quenched random forces. As in the Larkin model, random forces are considered
constant in the displacement direction and uncorrelated in the perpendicular
direction. The elastic energy density contains an harmonic part, proportional to
(∂_x u)^2, and an anharmonic part, proportional to (∂_x
u)^2n, where u is the displacement field and n>1 an integer. By heuristic
scaling arguments, we obtain the global roughness exponent ζ, the dynamic
exponent z, and the harmonic to anharmonic crossover length scale, for
arbitrary d and n, yielding an upper critical dimension d_c(n)=4n. We find
a precise agreement with numerical calculations in d=1. For the d=1 case we
observe, however, an anomalous “faceted” scaling, with the spectral roughness
exponent ζ_s satisfying ζ_s > ζ > 1 for any finite n>1, hence
invalidating the usual single-exponent scaling for two-point correlation
functions, and the small gradient approximation of the elastic energy density in
the thermodynamic limit. We show that such d=1 case is directly related to a
family of Brownian functionals parameterized by n, ranging from the
random-acceleration model for n=1, to the Lévy arcsine-law problem for n =
∞. Our results may be experimentally relevant for describing the
roughening of nonlinear elastic interfaces in a Matheron-de Marsilly type of
random flow.




Roughening of the anharmonic Larkin model
    A. B. Kolton
    December 30, 2023
=========================================







§ INTRODUCTION


The general study of elastic interfaces in random media is relevant for
understanding generic properties displayed by a variety of experimental systems,
and to successfully classify them into universality
classes <cit.>. Disorder adds pinning, metastability, and
generates macroscopically rough interfaces, while the elasticity tends to
flatten them. When driven by a uniform force F, extended interfaces display a
nontrivial zero-temperature depinning transition from a static to a sliding
regime at a threshold value F_c <cit.>, while for
finite temperatures a collective thermally activated creep motion persists below
the depinning threshold,
F<F_c <cit.>. At zero
force, the interface gets trapped in one of the many available deep metastable
states <cit.>. In all cases, disorder
induces statistically self-affine rough geometries in the putative steady
states.

Universality allows convenient minimalistic models to capture the relevant
disorder-elasticity interplay. A rather minimalistic, yet nontrivial model, is
the one describing the position of a directed elastic interface at a time t as
an univalued scalar displacement field u(x,t), where x ∈ℝ^d,
with d the internal dimension of the manifold (d=0 a particle, d=1 for a
string, d=2 for a sheet, etc.), D=d+1 being the dimension of the space. Such
elastic interfaces do not allow the formation of overhangs nor pinch-off
bubbles. Specifically, we consider an elastic interface, subject to a pinning
force g[u(x,t),x] with a given elastic energy E_el(u). If we assume a
nondriven overdamped relaxational dynamics at zero temperature, then the
equation of motion of the interface results in

    ∂_t u(x,t) = -δ E_el/δ u(x,t) + g[u(x,t),x].

A rather generic disorder is specified by the average over disorder realizations
of the pinning force g(x)=0, and its spatial autocorrelation
function g(x,u)g(x',u')=^2 δ(x-x')Δ_ξ(u-u').
Here, Δ_ξ(u) is some even function of u with ξ denoting its
shortest characteristic length. If we choose Δ_ξ(0)=1, then 
measures the strength of the disorder (see, for instance,
Ref. <cit.>).

The simplest harmonic form E_el = ∫ dx (c_2/2)(∂_x u)^2, with
c_2 ≥ 0 the elastic constant, leads to the celebrated (zero-temperature)
quenched Edwards-Wilkinson (QEW) equation <cit.>:

    ∂_t u(x,t) = c_2 ∂^2_x u + g[u(x,t),x].
 
The main difficulty of Eq. (<ref>) is the nonlinearity and
heterogeneity of the pinning forces g(u,x). This has led Larkin and coworkers
 <cit.> to approach the problem
perturbatively in the disorder such that, at first order, the pinning force can
be simply approximated by an x-dependent but u-independent quenched random
force, g[u(x,t),x] ≈ f(x), with

    f(x)   =    0, f(x)f(x')=^2 δ(x-x').

Replacing g by f in Eq. (<ref>), we obtain the linear
equation

    ∂_t u(x,t) = c_2 ∂^2_x u + f(x),

governing the dynamics of the so-called Larkin model (LM). Since
Eq. (<ref>) is translationally invariant, bulk disorder does not
pin the interface. However, it makes the interface evolve to a rough steady
state.

In spite of its simplicity, the LM has many interesting properties and is a
fundamental model in the theory of disordered elastic systems. The linear
Eq. (<ref>) is known to correctly describe the behavior of the
QEW interfaces, which evolve according to the nonlinear
Eq. (<ref>), at scales smaller than the so-called Larkin length
L_c ∼ (c ξ/)^2/4-d, for d < 4.  Above L_c the LM
solution crossovers to the so-called random-manifold regime. Beyond the LM
description, further crossovers at larger scales are still possible, depending
on the properties of Δ_ξ(u), and also on the drive
F <cit.>. In any
case, the LM is relevant to estimate the fundamental physical units for
expressing the large-scale solution of the full QEW problem. Quite remarkably,
L_c can be related to the macroscopic critical depinning threshold F_c ∼/L_c^d/2, the elementary pinning energy barrier U_c = F_c L_c^dξ, and the finite size crossover to
metastability <cit.>. Above
L_c, in the random-manifold regime, pinning and elastic energies scale in the
same way as ∼ U_c (L/L_c)^θ, with θ a characteristic
exponent, and the interface global width as ξ(L/L_c)^ζ, with ζ a
roughness exponent. It is worth also noting that the LM is  relevant per se if
the pinning forces have a large correlation length in the direction of interface
displacements, such that L_c is large compared to system size. This can be
achieved if ξ is very large, the system is elastically stiff or the disorder
very weak. Such situations can arise experimentally, displaying a finite-size
pinning crossover <cit.>.

Interestingly, the physics of the LM is not only relevant for describing small
length scales. Pinning forces such as f(x) can be also generated via
coarse-graining at large length scales L ≫ L_c and dominate over other
forces in driven disordered elastic systems. This happens if the pinning force
correlator Δ_ξ(u) has a periodicity <cit.>. One simple
example is a QEW interface, Eq. (<ref>), in a disordered medium
with periodic boundary conditions in the direction of displacements, and an
additional driving force F. Another less trivial example, is a one-dimensional
elastic chain of particles (or interfaces) with average separation a_0, driven
in a completely uncorrelated random potential <cit.>. Such
models are often used to study friction <cit.>, but the
physics is also relevant for charge density waves
depinning <cit.>. In these cases, the large-scale
roughness is described by the solution of Eq. (<ref>).

The LM being a fundamental model in the theory of disordered elastic systems, it
is worth analyzing how their properties change under the influence of additional
physically motivated terms in Eq. (<ref>). In this paper, we
focus in the influence of anharmonic corrections to the elasticity. To motivate
the introduction of such corrections, we note that the LM, being linear, can be
easily solved for a particular realization of f(x) in any dimension d, and
averages over realizations are straightforward. Much universal information can
be obtained by studying the critical nonstationary relaxational dynamics, from a
flat initial condition at the origin of displacements, i.e., u(x,t=0)=0. Since
this process is dominated by a single dynamical growing length-scale l(t), for
long enough times before global equilibration.

We define the structure factor of the manifold as S(q,t)≡|û(q,t)|^2, where  is the Fourier transform of u(x,t). Using that
û(q,t=0)=0 we get,

    S(q,t) ∝^2/c_2^2 q^4(1-e^-c_2q^2 t)^2.

Then, the structure factor satisfies the general scaling S(q,t)∼
q^-(d+2ζ)G(q t^1/z), with G(y)=y^2(ζ-ζ_s) for y ≫ 1,
and G(y)=y^d+2 ζ for y ≪ 1  <cit.>. In our case,
it is easy to check that the so-called global roughness exponent is
ζ=(4-d)/2, and coincides with the spectral roughness exponent
ζ_s, while the so-called dynamical exponent, related to the growing
length-scale l(t) ∼ t^1/z, is z=2. In the steady state, roughly
reached at times t such that l(t) ∼ L, S(q,t→∞)∼
q^-(d+2ζ). For an interface of size L, this gives the global squared
width with respect to the center of mass position W^2 ≡[u(x,t)-v_cmt]^2=∫ dq S(q,t) ∼ L^2ζ, with
v_cm = L^-d∫ dx f(x) ∼ L^-d/2 the finite-size residual
center of mass velocity. The interface then becomes macroscopically self-affine,
with exponent ζ. That is, the rescaling x' → b x and u' →
b^ζ u leads to a statistically equivalent interface.

Interestingly, in the d=1 LM, we have ζ_s=ζ>1. This situation, which
has been called the super-rough case in Ref. <cit.>,
has a physical peculiarity. It implies that the harmonic elastic approximation
in Eq. (<ref>) to the local elastic couplings must break down in
the thermodynamic limit, since W/L ∼ L^ζ-1→∞ as L→∞. Local elongations are thus not  bounded in the thermodynamic limit. A
similar situation occurs for the roughness exponent at the depinning threshold
for the driven QEW model, where ζ≈
1.25 <cit.>. To remedy this situation, in
Ref. <cit.> ad hoc anharmonic corrections to the elastic energy
were introduced, such that

    E_el(u) = ∫ dx [c_2/2 (∂_x u)^2 +
            c_2/2 (∂_x u)^2],

with =2,3,4..., and constants c_2>0, c_2>0. It is worth
noting that this kind of correction is Hamiltonian, convex, and being a
correction to the elasticity only, translational invariant. In particular, note
that the presence of the anharmonic correction breaks the tilt-symmetry of the
full QEW equation, since both g(u,x) and the harmonic elasticity, in contrast
with the (∂_x u)^2 term for n>1, are statistically invariant
by the transformation u → u-s x_α, with α denoting any of the
d internal directions and s the parameter measuring the tilt deformation.

For large , Eq. (<ref>) is equivalent to
impose a hard-constraint to local elongations, an usual modeling of directed
polymers <cit.>. The proposed nonquadratic term succeeds to
save the elastic approximation in the thermodynamic limit at the depinning
transition of the QEW model, with a new (“physical”) roughness exponent ζ
= ζ_s ≈ 0.63 for all >1 <cit.>. Quite
remarkably, this value is in perfect agreement with the roughness at the
depinning transition of the paradigmatic quenched Kardar-Parisi-Zhang (QKPZ)
model <cit.>, which is inherently a nonequilibrium effective equation
that cannot be derived directly from a Hamiltonian or free energy. The
anharmonic model hence allows us to study a model with a Hamiltonian, and a
well-defined equilibrium state at F=0, which nevertheless spontaneously
generates the ubiquitous Kardar-Parisi-Zhang (KPZ) <cit.> term, when it
is driven by a force F. Thus, one may ask whether the breakdown of the small
gradient approximation in the elastic energy of the d=1 Larkin model
[Eq. (<ref>)] is, similarly, protected by introducing a nonlinear
elasticity of the form proposed in Eq. (<ref>), and how
ζ and ζ_s would change upon its addition. Solving the resulting
“anharmonic Larkin model” would also allow us to find a possibly modified
Larkin length L_c and related quantities, which are fundamental to estimate
both the critical force and the crossover length to the random manifold regime,
for the anharmonic depinning model defined by Eqs. (<ref>)
and (<ref>). The depinning transition of such model was
analyzed in Ref. <cit.>.

Motivated by the above phenomenology, in this work, we study the Larkin model
with anharmonic elasticity, for a general dimension d and ≥ 2. We
obtain the global roughness exponent as a function of  and d, and
also describe the crossover from harmonic to anharmonic regimes of roughness,
when two terms, one with =1 and another with >1 coexist. For
the special d=1 case, where the small elastic deformation approximation is
compromised, we show that, unlike the depinning model, the anharmonic correction
is never able to reduce ζ below unity, even in the large
 almost hard-constraint limit. Moreover, we show that for all
>1 the d=1 interface displays anomalous scaling properties, with
ζ_s > ζ≥ 1, the so-called “faceted regime” in
Ref. <cit.>. Interestingly, we show also that this
anomalous case is closely connected to an -parameterized family of
Brownian functionals which interpolate between the random-acceleration process
for =1, to the arcsine law Lévy stochastic process, for
=∞. For d>1, however, we find ζ≤ 1.

This article is organized as follows: in Sec. <ref> we describe the
anharmonic Larkin model, the observables of interest, and the methods. In
Sec. <ref> we derive, via heuristic arguments, the global
roughness and dynamical exponents, and also the harmonic-anharmonic crossover
length. In Sec. <ref>, we numerically test these predictions for
d=1, solving both the relaxational dynamics of an interface as a function of
time, and the statics, for two different boundary conditions. In
Sec. <ref>, we discuss the relation for the anomalous d=1 case
with a family of Brownian functionals. Finally, in Sec. <ref> we
present our conclusions along with new open questions, and we suggest some
possible applications for our results.



§ MODEL AND OBSERVABLES


We consider the anharmonic Larkin model using Eq. (<ref>) in the
Larkin approximation g(u,x) → f(x), and Eq. (<ref>)
for the elastic energy. The resulting equation of motion reads

    ∂_t u(x,t) =
    c_2 ∂^2_x u + c_2∂_x 
    [(∂_x u)^2-1] + f(x),

where >1. We will call Eq. (<ref>) the
“anharmonic Larkin model” (ALM). Fig. <ref> schematically
represents the model for the particular d=2 case.



Let us consider an interface of linear dimension L. The time-dependent global
width for the relaxation process W(L,t) is defined by the mean-squared
fluctuations around the center of mass,

    W^2(L,t) = 1/L^d∫ dx  (u-v_cmt)^2,

where  [from
Eq. (<ref>)] is the finite-size disorder-dependent but
constant center of mass velocity, that may be observed for periodic or free
boundary conditions. It is also useful to compute the time-dependent structure
factor S(q,t,L), as the interface displacement spatial power spectrum,

    S(q,t,L) = |û(q,t)|^2,

with û(q,t) the Fourier transform of u(x,t) as defined just before
Eq. (<ref>).

Solving Eq. (<ref>) from a flat initial condition
 allows us to extract several universal exponents
characterizing the generically expected critical relaxational dynamics. Such
dynamics are expected to be controlled by a single growing correlation length
scale,

    l(t) ∼ t^1/z,

with z defining the dynamical exponent. Such scaling is expected to be valid
at times larger than the microscopic time scale, but smaller than the time at
which the system becomes completely correlated, i.e., when l(t) ∼ L.
However, from W^2(L,t), and S(q,t,L) we can define in principle three
different roughness exponents (z, ζ, and ζ_s) characterizing the
random geometry below l(t) <cit.>.

The global roughness exponent ζ is defined from

    W(t,L) ∼ t^ζ/zW̃[L/l(t)],

with

    W̃(x) ∼
        x^ζ    if  x ≪ 1 
    const    if  x ≫ 1,

or more directly from the stationary interface width

    W(t→∞,L) ∼ L^ζ.


Both the global ζ and the spectral ζ_s roughness exponents can be
obtained from the expected scaling

    S(q,t) = q^-(2ζ+1) S̃(qt^1/z),

where S̃(x) is the scaling function and has the general form

    S̃(x) ∼
        x^2(ζ-ζ_s),      if  x ≫ 1 
    
        x^2ζ+1,      if  x ≪ 1.

The stationary limit is reached when the correlation length l(t) reaches L,
at times of order L^z. Thus,

    S(q) ∼ q^-(2ζ_s+1) L^2(ζ-ζ_s).


At zero temperature, the nonstationary scaling behavior from
Eqs. (<ref>), (<ref>), (<ref>), and
(<ref>), and also the stationary scaling described by
Eq. (<ref>), are verified by the QEW [Eq. (<ref>)]
at or above the depinning threshold F ≥ F_c <cit.>,
and in particular by the harmonic Larkin model [Eq. (<ref>)] for
any force, as shown in Eq. (<ref>). In all these cases, a single
growing correlation length controls the relaxation toward a unique self-affine
stationary state, without memory of the initial condition. Different roughness
exponents are obtained in each case. For instance, the d=1 QEW equation at
F_c has ζ =  ζ_s ≈ 1.25 and z ≈
1.433 <cit.>, while at F>F_c, it crossovers to the
Edwards-Wilkinson exponents ζ = ζ_s = 1/2 and z=2. The harmonic
Larkin model, however, as discussed in Sec. <ref>, has
ζ=ζ_s=(4-d)/2 and z=2. In the following sections, we will show that
the ALM also displays the same scaling forms, and we will obtain their exponents
ζ, ζ_s, and z for =1,2,...,∞.



§ SCALING ARGUMENTS


The dynamical and roughness exponents of the EW, the LM and other linear models
can be successfully obtained by simple scaling arguments (see
Ref. <cit.> for many more relevant examples). Although such approach
may fail in general (for instance, KPZ, QKPZ, or QEW equations require
renormalization group calculations instead), we will employ the same naive
methodology for the ALM nonlinear dynamics in
Eq. (<ref>). For simplicity, at first we consider only
the nonlinear elastic term. Putting c_2=0, Eq. (<ref>)
reduces to

    ∂_t u(x,t) = c_2∂_x
      [(∂_x u)^2-1] + f(x).


Hence, to obtain the scaling exponents, we follow the standard naive procedure,
which consists in rescaling  space and time as

    x    → x' ≡ b x,
    
      u    → u' ≡ b^ζ u,
    
      t    → t' ≡ b^z t.

Then, making the strong assumption that the elasticity and pinning parameters
are not changed by rescaling, after this transformation,
Eq. (<ref>) results in

    b^ζ-z∂_t u =
    
      b^(ζ-1)(2-1)-1 c_2∂_x
        [(∂_x u)^2-1] + b^-d/2 f,

where we have used that the random uncorrelated pinning forces scale as a
d-dimensional random walk.

We finally require the resulting equation must be invariant under the
transformation, which leads to the Flory-like exponents,

    ζ(d,)    = 4-d/4-2, 
    
      z(d,)    = ζ(d,) + d/2,

for =1,2,...,∞, and general dimension d. 

As expected, for =1, we recover the LM exponents, ζ=(4-d)/2 and
z=2, and the known upper critical dimension d_c(=1)=4. From the
prediction in Eqs. (<ref>) and (<ref>), it is
worth noting the following:


  (1) The upper critical dimension increases with , as
d_c()=4.

  (2) In the hard-constraint case,  corresponding to , we get, for a fixed dimension d, that ζ→ 1, and z →
1+d/2.

  (3) The anharmonic elastic energy density scales as
L^2[ζ(d,)-1], thus for >1 it is thermodynamically
bounded only for d>2, since ζ(d > 2,) < 1. The d=2 case is
always marginal, ζ(d=2,) = 1, while the d<2 case diverges, since
ζ(d < 2,) > 1.

  (4) The spectral roughness exponent ζ_s does not come out from the
arguments.


If the predicted exponents are valid, then we arrive to the striking conclusion
that the unbounded local displacements predicted for the d=1 LM remain
unbounded, in spite of the anharmonic elasticity, for any value of .
Therefore, the anharmonic correction to the elasticity can not save the small
local deformation assumption behind the gradient expansion of the elastic energy
density, even in the n →∞ hard-constraint limit, in the large-size
limit. This is in sharp contrast to what happens for the d=1 QEW equation at
depinning, where the elastic approximation is saved by the very same correction,
making the depinning roughness exponent to change from ζ_dep≈ 1.25 to the n-independent (n>1) value ζ_dep≈
0.63 <cit.> [Note that if the ALM is considered as a
short-scale approximation of the QEW with anharmonic
corrections <cit.>, it applies only below the corresponding finite
Larkin length, calculated in the Appendix. In this case, the elastic
approximation is not necessarily compromised in the large scale limit in spite
of the super-rough (anharmonic) Larkin regime.].

Let us now consider both the usual harmonic elasticity term and the nonlinear
correction together, i.e., c_2 > 0 and c_2 > 0, for a given
>1, in Eq. (<ref>). In this case, the harmonic
term should dominate the behavior for small distortions, so there might be an
anharmonic crossover length L_anh between two regimes of
roughness, from the harmonic LM (c_2>0, c_2=0) to the previously
analyzed purely anharmonic ALM (c_2=0, c_2>0) universality classes.
Heuristically, we can propose that, for a length-scale l>L_anh, the
purely anharmonic ALM displacement behaves as u ∼ u_anh
(l/L_anh)^ζ(d,); while for l<L_anh, we have the
harmonic result u ∼ (/c_2) l^ζ(d,=1) = (/c_2)
l^(4-d)/2. Sharply matching these two regimes at L_anh implies

    u_anh=/c_2 L_anh^(4-d)/2.

We propose that the crossover occurs when the two elastic
energy terms are equally important, c_2
(u_anh/L_anh)^2=(c_2/)
(u_anh/L_anh)^2, thus,

    u_anh =
        ( c_2/c_2)^1/2(-1) 
      L_anh
                 = /c_2 L_anh^(4-d)/2,

which leads to the crossover length,

    L_anh = (c_2/)^2/2-d ( c_2/c_2)^1/(2-d)(-1).

This means that, for  the combined elasticities, we would have, in the
stationary limit of a large system L ≫ L_anh, a crossover behavior
as a function of the length-scale in the correlation functions. For instance,
the mean-square width [Eq. (<ref>)] is expected to satisfy

    W^2(L,t→∞) = L^2ζ(d,n=1) w(L/L_anh),

with w(x)=x^2[ζ(d,n>1)-ζ(d,n=1)], for x ≫ 1; and , for x ≪ 1.

In the following two sections, we numerically check  all these heuristic scaling
predictions, for the special d=1 case.



§ NUMERICAL RESULTS FOR D=1


To test the validity and robustness of the scaling predictions, we have found
convenient to analyze separately the static solution, using appropriate boundary
conditions. With these numerical methods, we are able to access all the
exponents ζ, ζ_s, and z, and the crossover length
L_anh.



 §.§ Static solution


We study the static geometric properties of the ALM,  by solving
Eq. (<ref>) in the stationary limit, where the elastic
and
random forces exactly balance, i.e.,

    c_2 ∂^2_x u + c_2∂_x [(∂_x u)^2-1] =
      -f(x).




  §.§.§ Global roughness exponent

Since the anharmonic term dominates the geometry at large length scales, we
first focus on the purely anharmonic ALM by fixing c_2=0 in
Eq. (<ref>). This is then straightforward to solve:

    u(x) =  -∫_0^x dx'
            [∫_0^x'dx”f(x”)/c_2-(∂_x
                  u(0))^2-1]^1/2-1,

where we have used the particular boundary condition u(0,t)=0. As elastic
interactions are short ranged, this particular boundary condition does not alter
the large-scale scaling properties we are interested in.

To numerically evaluate Eq. (<ref>), we simply discretize the
d=1 interface in segments of size δ x, and perform the sums

    u_i = -∑_j=0^i-1[∑_k=0^j f_k/c_2]^1/2-1,
       for  i>0,

with u_0=0 fixed, and f_0 ≡ -∑_i=1^L f_i; taking care of the
statics of the whole polymer of size L in the discrete description, and
accounting for the term  in the continuum
Eq. (<ref>) [Note that the constraint force f_0,
needed to fix the boundary conditions in the first displacement u_0 of the
polymer, is not described by the bulk random forces, with the statistical
properties in Eq. (<ref>).]. Note that, without loss of
generality, we have taken δ x=1 exploiting that the pinning force is
completely uncorrelated in the internal coordinate. We have drawn f_k from a
uniform distribution with zero mean and variance ^2 = 1/12.

Equation (<ref>) can  then be solved very efficiently for large
system sizes, by using parallel random number generators and parallel prefix-sum
algorithms implemented in massively parallel coprocessors, such as graphic
cards. As the same scaling analysis leading to the prediction in
Eq. (<ref>) applies to Eq. (<ref>) or
Eq. (<ref>), regardless of the  fixed-end condition u_0=0, we
expect to get ζ(d=1,) from the numerical evaluation of
Eq. (<ref>). Since the system is not translational invariant along
x, it is convenient to directly measure the last-point scaling, for which we
expect u_L^2∼ L^2ζ(d=1,), for large enough L.

In Fig. <ref> we show there is an excellent agreement between
predicted behavior and the obtained results for various values of , for
system sizes larger than around 10^2, and up to 10^8 string elements.





  §.§.§ Harmonic to anharmonic crossover

To study the crossover between the LM and the ALM, we now consider the static
condition Eq. (<ref>), with c_2 ≠ 0 and c_2≠
0, i.e., by including both the harmonic and the anharmonic terms. Integrating
it, we get

    c_2 ∂_x u + c_2 (∂_x u)^2-1   =   
      c_2 ∂_x u(0) + c_2 [∂_x u(0)]^2-1
       -   ∫_0^x  dx” f(x”) ,

which is an odd (2 - 1)-degree polynomial in ∂_x u. Since
c_2>0 and c_2>0, it has only one real root for each x. By
integrating this
x-dependent root we obtain u(x). We perform this task numerically, for the
discretized
interface u_i, as before; using that c_2 ∂_x u(0) + c_2[∂_x u(0)]^2-1∼ -f_0 in this mixed case. Finally, to test
the prediction for L_anh [Eq. (<ref>)], we use the
expected crossover scaling for the polymer end-point,

    u_L^2∼ L_anh^4-d h(L/L_anh),

with h(x)∼ x^4-d, for x ≪ 1; and h(x) ∼ x^2ζ(d=1,),
for x ≫ 1.

We have tested the prediction for L_anh from
Eq. (<ref>), for the particular d=1 and =2 case,
where Eq. (<ref>) becomes an analytically solvable cubic
polynomial. The prediction of Eq. (<ref>) for this case is

    L_anh = (c_2/)^2 2 c_2/c_4
      = 2 c_2^2 ^-2 c_4^-1.

In Fig. <ref>, we show the results for the rescaled families
of curves, u_L^2/L_anh^3 versus L/L_anh,
combining c_2 = 1  with various values for c_4, from c_4 = 1 to c_4=
0.0001. The master curve strongly supports the validity of
Eq. (<ref>).





 §.§ Dynamic Solution


To test dynamical scaling (i.e., involving the time variable), we have performed
numerical simulations of Eq. (<ref>) in d=1 with
periodic boundary conditions, and we have averaged the results over many
disorder realizations. We have implemented this by using a spatial finite
difference scheme, and we have solved the resulting system of equations
following the standard Euler method. If the discretization is δ x, such
that u_i ≡ u(x=iδ x), for i=1,2,…,L, then
Eq. (<ref>) can be approximated by

    ∂_t u_i    =
      c_2/(δ x)^3(u_i+1+u_i-1-2 u_i) 
       + c_2/(δ x)^2+1[
      (u_i+1-u_i)^2-1 -
      (u_i-u_i-1)^2-1] 
       + f_i,

with u_0 ≡ u_L and u_L+1≡ u_1. The d-dimensional
generalization is straightforward. We draw the random forces f_i from a
uniform distribution with zero mean and variance ^2 = 1/12, and we
solve Eq. (<ref>) starting from a flat
configuration
u_i(t=0) = 0.



  §.§.§ Time dependent structure factor

The structure factor S(q,t,L) allows us to obtain ζ, and in particular,
the dynamical exponent z and the spectral roughness exponent ζ_s, from
the expected scaling of Eq. (<ref>).

In Fig. <ref> we show the corresponding S(q,t) for c_2
= 1, c_2=0, L=262144, and different times t. The main panel
shows the raw data (displayed in the inset) rescaled using
Eq. (<ref>). The result is perfectly consistent with the scaling
and confirms the exponents ζ=3/2 and z=2 of the analytical solution of
Eqs. (<ref>) and (<ref>) for d=1 and
=1.



To contrast with the pure harmonic case considered before, we now consider the
purely anharmonic case with c_2=0 and c_2=1 using =2. In
Fig. <ref>, we show the corresponding S(q,t) for
different times. Interestingly, we can observe that large q behavior does not
follow a master curve, but there is a downward drift by increasing times, which
differs from the harmonic result shown in the inset of
Fig. <ref>. Moreover, fitting the large q power-law
behavior for a long time t, yields a different power-law exponent than fitting
the power-law envelope for all the times. In Fig. <ref> we
show that, nevertheless, the scaling of Eq. (<ref>) works
perfectly
by using the global roughness exponent ζ(d=1,=2) = 7/6 and the
dynamical exponent z(d=1,=2) =5/3, both predicted by dimensional
analysis. However, in spite of the perfect collapse using the predictions of the
scaling arguments of Eqs. (<ref>) and
 (<ref>)
for d=1 and =2, the master curve displays a clearly larger and well
defined roughness exponent. From Eqs. (<ref>) and
 (<ref>), we can identify this new exponent with the
spectral exponent ζ_s(d=1,=2)≈ 1.39.







  §.§.§ Critical exponents as a function of 

Repeating the above procedure, we have  fitted z, ζ, and ζ_s for
different values of , from =2 to =100. In
Fig. <ref> we show the corresponding results, where it can be
appreciated that numerical simulations perfectly agree, within the error
bars [We use different criteria to estimate the error bars for each
exponent. In the case of ζ, we estimate the uncertainty fitting the data
of Fig. <ref>. The error bar corresponding to ζ_s is obtained
by fitting the slope in the rescaled structure factors, shown in
Fig. <ref>. For z we make a conservative estimation of
its uncertainty by determining the range where a good collapse can be seen.],
with the predicted exponents ζ(d=1,) and z(d=1,) (shown by
solid lines), from Eqs. (<ref>) and (<ref>).
In all cases, the fitted spectral exponent is found to be ζ_s(d=1,)
> ζ(d=1,). We also observe that all exponents tend to converge to a
well defined finite value in the large  limit.  The results also show a
roughly constant difference, ζ_s(d=1,)-ζ(d=1,) ≈
0.23 for >1, i.e., for all the anharmonic cases. This is very
interesting, as we do not have any analytical prediction for
ζ_s(d=1,).





§ DISCUSSIONS


We find an excellent agreement of the general scaling predictions of
Sec. <ref> with static and dynamical numerical simulations
of a d=1 interface as a function of the anharmonicity parameter .
These predictions include the global roughness exponent ζ, the dynamical
exponent z, and the harmonic to anharmonic crossover length
L_anh. We also obtain numerically the spectral roughness exponent
ζ_s, for which no analytical prediction is available, and we find that it
is a different exponent, at least for d=1 and n>1. In the following sections
we discuss possible connections of the ALM with other problems, and then
we discuss the anomalous scaling.



 §.§ Family of Brownian functionals


It is interesting to note that the static solution u_L of
Eq. (<ref>) is directly related to a family of Brownian
functionals <cit.> parametrized by ,

    u_L = -∫_0^L dx sign[B(x)] |B(x)|^1/2-1,

with the Brownian path B(x) given by

    B(x)=∫_0^L dx  f(x).

From the statistical properties assumed for the elementary pinning forces
[Eq. (<ref>)], B(x) is clearly a Wiener process.

Let us now analyze the family of Brownian functionals. The =1 case,

    u_L =∫_0^L dx  B(x),

simply corresponds to the random-acceleration process (see, for instance,
Refs. <cit.>). This well-known
process has a Gaussian end-point distribution with a particular super-extensive
variance,

    P(u_L,=1)    =   exp[-u_L^2/2σ_L^2]/√(2πσ_L^2),
    
    σ_L^2    ∼    L^3.

Hence, u_L ∼ L^3/2 is in perfect agreement with the scaling result u_L
∼ L^ζ, with the global roughness exponent
 predicted in
Eq. (<ref>).

On the other extreme of the ALM family, when →∞ we get

    u_L =∫ dx sign[B(x)],

which corresponds to the so-called arcsine-law Lévy process <cit.>.
We have also an exact result for this case,

    P(u_L, =∞)=1/π1/√((L-u_L)(u_L+L)).

Its variance grows as ∼ L^2 implying that , which is
also in perfect agreement with our predicted
 in Eq. (<ref>).

For intermediate cases, 1<<∞, P(u_L,) is, to the best our
knowledge, not known analytically. In Fig. <ref> we show
P(u_L,) obtained numerically from Eq. (<ref>), for several
values of . These end-point distributions  characterize the family of
Brownian functionals describing the d=1 ALM. We observe that, in general, all
>1 distributions are symmetric but non-Gaussian, with twin peaks far
from the origin that grow with increasing , and diverge when →∞, as the Lévy arcsine process is approached.



Different methods, such as the Feynman-Kac formula (see, for instance,
Ref. <cit.>), could be used in principle to get
P(u_L;) analytically for general , and to compare with our
numerical results of Fig. <ref>.



 §.§ Relation to variants of the KPZ equation

Following the connection made in Ref. <cit.>, between the QKPZ
equation and the QEW equation with anharmonic elasticity at depinning, one may
ask whether the ALM roughness properties can be related to other nonlinear
roughening equations. At this respect, it is worth noting that the KPZ equation
with temporally correlated noise,

    ∂_t u = c_2 ∂_x^2 u + λ (∂_x u)^2 + η(x,t),

with the noise field such that ⟨η(x,t) ⟩=0, and ⟨η(x,t)η(x',t') ⟩ = |x-x'|^2ψ-d |t-t'|^2ϕ-1, was found to
satisfy the scaling relation <cit.>,

    z(1+2ϕ)-2ζ-d+2ψ=0.

Interestingly, for quenched spatially uncorrelated
noise <cit.> (ϕ=1/2 and ψ=0), we get
2z-2ζ-d=0, which is identical to Eq. (<ref>). This can be
related to the fact that Eq. (<ref>) is derived from
the nonrenormalization of the temporal component of the noise correlator, while
Eq. (<ref>) was derived assuming, analogously, the no-scaling of
the random force amplitude. Since the global roughness and dynamic exponents of
the ALM are parametrized by n, however, one can ask whether a value of n
exists such that the ALM and the quenched noise KPZ share the same exponents
ζ, ζ_s, and z.

At this respect, numerical simulations of temporally correlated KPZ for d=1
and λ>0 report, for the  ϕ=1/2 quenched-noise limit, ζ≈ 1.06, z ≈ 1.54 <cit.> and  ζ≈ 1.07, z
≈ 1.15 <cit.>. Exponents are thus roughly consistent with
the n ≈ 4–5 ALM, where ζ≈ 1.06 ± 0.01 and z ≈
1.56. Interestingly, in Ref. <cit.> it is shown that configurations
display facets, with an spectral exponent ζ_s ≈ 1.5, which is
slightly larger than the ζ_s ≈ 1.29 ± 0.01 we find for n ≈
4–5 ALM. More recent results for the KPZ with temporally correlated noise
approaching the limit ϕ=1/2, report ζ_s ≈ 1.3, showing a closer
agreement with the n ≈ 4–5 universality class of the
ALM <cit.>. The quenched noise KPZ with λ<0, on the other 
hand, has been associated (though for more general random forces) to faceted
interfaces <cit.>, similar to the ones we show for d=1 in
Fig. <ref>. The possible connection of the ALM for general n with
n-dependent extensions of Eq. (<ref>) remains an interesting
open question.

Finally, it is also worth pointing out that the d=1 connection of the ALM with
Brownian functionals of the previous section also suggests, from a different
perspective, a connection with other KPZ variants. The Feynman-Kac formula is
indeed connected to a family of quantum Hamiltonians describing a single
particle in a family of potentials pU(X) with p a parameter (see, for
instance, Ref. <cit.>). The quantum propagator satisfying each
imaginary-time Schrödinger-like equation can then be mapped, by a suitable
“Cole-Hopf” transformation, to a particular forced KPZ
equation <cit.>. From Eq. (<ref>) we have
U(X)=sign[X] |X|^1/2-1 so the potential is
parametrized by both p and . It would be interesting to test these
connections further.





 §.§ Two different roughness exponents

The ALM displays interesting anomalous scaling properties. In particular, in the
steady state, the spectral roughness exponent does not coincide with the global
roughness exponent (ζ_s>ζ for >1) in d=1. This is in sharp
contrast to what happens for the LM and many other linear surface growth models,
such as the EW equation, where interfaces are characterized by a single
roughness exponent, and the usual Family-Vicsek scaling applies <cit.>.
It is worth noting, however, that also nonlinear models such as the KPZ equation
with additive noise <cit.>, or the zero temperature QEW
equation at the depinning threshold <cit.>, seem to display a
normal Family-Vicsek scaling of the two point correlations functions, such that
ζ≈ζ_s, with ζ_s ≈ 1/2 in the former and ζ_s
≈ 1.25 in the latter, in d=1. Therefore, nonlinearity, either in the
elasticity or in the disorder, is not sufficient to produce the anomaly, even
combining it with super-rough behavior ζ_s>1, as it occurs in the QEW at
depinning. We also note that breaking the tilt symmetry (as for instance in the
KPZ equation), or breaking the statistical tilt symmetry (as for instance in the
d=1 QKPZ equation at depinning where ζ_s ≈ζ≈
0.63 <cit.>), does not seem to be sufficient for
producing the anomaly either.

However, other models are known to display the ζ≠ζ_s anomalous
behavior. The anomalous case with ζ_s<1, so called “intrinsic,” has
been, for instance, observed <cit.>. The case
ζ_s>1 has been observed, however, in the Sneppen model A of self organized
criticality <cit.>. Interestingly, in Ref.
<cit.>, the ζ≠ζ_s >1 behavior has been
associated to generic “faceted” stationary interfaces, and a simple model of
random ± 1 slopes used to illustrate the association. In
Fig. <ref>, we show that this is indeed qualitatively the case for
the ALM family parameterized by  in d=1. As far as we can tell, the
generic  origin of the ζ≠ζ_s behavior, either in the ζ_s<1
(“intrinsic”) or ζ_s>1 (“faceted”) remains open, however.

We speculate that the observed anomalous roughness property ζ≠ζ_s
of the d=1 ALM should be related to the shape of P(u_L,) for
>1. Indeed, only for the LM (=1) the distribution is Gaussian
and ζ=ζ_s. For >1 it is non-Gaussian and presents two twin
peaks and, as expected, the larger and sharper they are, the more faceted the
individual configurations look like. Since the non-Gaussian bimodal shape of
P(u_L,) presents  scale invariance, the associated faceted structure
should be also scale invariant and thus visible in the structure factor shape.
This makes us conjecture that ζ_s(d=1,) must be somehow coded in
P(u_L,). It would be useful to test this hypothesis for d=1 by
analytical calculations using, for instance, Brownian functionals techniques.

For d>1, however, it would be useful to go beyond heuristic methods and use
the Gaussian variational method (see, for instance, Ref. <cit.>),
or renormalization group techniques, to calculate both ζ and particularly
ζ_s, which we were unable to estimate by simple scaling arguments. One
important question is when ζ_s ≠ζ, and to know if ζ_s is
universal. In other words, the question is whether ζ_s is needed to make a
finer classification into (anomalous) universality classes.



§ CONCLUSIONS


We have studied the roughness and dynamical properties of the anharmonic Larkin
model. By Flory-like scaling arguments, we have obtained the global ζ
roughness exponent, the dynamic exponent z, and the harmonic to anharmonic
crossover length scale for arbitrary both d and n. An excellent agreement is
obtained by comparing with numerical calculations in d=1, and we have showed
that, in this case, the model directly relates to a family of Brownian
functionals parameterized by n; ranging from the random-acceleration model
(n=1) to the Lévy arcsine-law problem (n = ∞).  It would be
interesting to exploit this connection, not only to find analytical methods of
solution, but to link or map the d=1 ALM model to other problems.

From the analytical and numerical calculations, we have found an intriguing
anomalous scaling for the d=1 case. The spectral roughness exponent ζ_s
is different from the global ζ, and ζ_s > ζ >1, for n>1. On
the one hand, this implies that the small gradient expansion for the elastic
energy is always compromised in the thermodynamic limit, even in the  n ≫ 1
hard-constraint limit. On the other hand, that we need two roughness exponents,
already at the level of the two point correlation functions, to describe the
self-affine geometry of the interface. It would be interesting to see if the
“faceted” scaling anomaly ζ_s ≠ζ persists at larger dimensions,
and in particular, to find an analytical estimate of ζ_s as a function of
dimension d and the anharmonicity parameter n using, for instance,
variational or renormalization group approaches. This would allows us to test
the universality of ζ_s, and to point out the mechanism originating the
anomaly in general cases.

Finally, although the Larkin model is usually a local approximation for more
realistic models of disordered elastic systems with many metastable states, it
would be nevertheless interesting to check some of our results experimentally.
This is in principle  possible in systems with anisotropically correlated random
forces. Indeed, interesting anomalous scaling was found for interfaces with
“columnar noise,” both theoretically <cit.> and experimentally
<cit.>. Finite systems with either stiff elastic couplings or very
weak disorder can also comply with the Larkin approximation for disorder at the
relevant scales. Larkin random forces could be also spontaneously generated, by
coarse-graining, well beyond the Larkin length in large driven periodic systems
such as elastic chains <cit.>. In all these cases, the
additional anharmonicity needed for realizing the ALM described here may arise
from some nonlinear local interaction breaking the tilt symmetry of the clean
(i.e. nondisordered) system. Directed polymers or membranes in a layered
Matteron-de Marsily scalar flow field  <cit.> may be a realization
of the ALM schematically represented in Fig. <ref>.


We thank S. Bustingorry, A. Rosso,  C. Texier, and T. Giamarchi for useful and
motivating discussions. This work was partly supported by Grants No.
PICT2016-0069/FONCyT and No. UNCuyo C017, from Argentina. We have used Mendieta
Cluster from CCAD-UNC and GPGPU-CAB Cluster from CAB, which are part of
SNCAD-MinCyT, Argentina.




§ LARKIN LENGTH WITH PURE NON-HARMONIC ELASTICITY


If we see Eq. (<ref>) as a short-scale approximation of
the full model of Eq. (<ref>) with
Eq. (<ref>), then it is useful to get the corresponding
(“anharmonic”) Larkin length.

By balancing the elastic force with the pinning force of a piece of linear size
l, we get c_2 u^2-1 l^-2(2-1)-1≈
l^-d/2. We can define L_c(d,) as the length l corresponding to
u=ξ,

    ξ = (/c_2)^1/2-1
          L_c^ζ(d,),

so we obtain

    L_c(d,) = [ξ^2-1(c_2/)]^2/(4-d).

Thus, for =1, we recover the well-known harmonic result
L_c(d,1)=(ξ c_2/)^2/(4-d). For =∞,
however, L_c →ξ. The length L_c(d,n) marks the crossover to the
random-manifold regime of the purely anharmonic interface. While the Larkin
regime is described by the exponents ζ(d,n)
[Eq. (<ref>)] and z(d,n) [Eq. (<ref>)],
regardless of the presence of the driving force F, the random-manifold
exponents do change with F. The equilibrium random manifold exponents (present
at equilibrium or in the creep regime at intermediate
scales <cit.>) are expected to be the same as those for
the harmonic elasticity or QEW equation, while the depinning random manifold
exponents (present at depinning or in the creep regime at large scales for
vanishing velocities  <cit.>) coincide with the ones of
the QKPZ equation <cit.>.


