![FIG. 4. (Color online) (a) The ratio between the diagonal and off diagonal elements of the Hessian |H0|2|H1| as a function of r . (b)–(e) Coefficients Al,B,Cl,Dl defined in Eq. (11) as a function of r . The solid blue lines depict the coefficient values for l = 1, whereas the black doted lines for l = 2 and the red dashed ones for l = 3. The cyan (shaded) area marks the region of applicability of the RWA as derived from (a). The vertical lines correspond to the values of r we study with the DNLS model. Note that three of them namely r2,r3 and r4 are the same as those studied in the previous section [Fig. 1(b)] whereas r45 lies between r4 and r5. The insets provide a zoom into the region of interest.](/figures/fig-4-color-online-a-the-ratio-between-the-diagonal-and-off-3l26h38g.webp)
PHYSICAL REVIEW E 92, 042905 (2015)
Dynamics of nonlinear excitations of helically confined charges
A. V. Zampetaki,
1
J. Stockhofe,
1
and P. Schmelcher
1,2
1
Zentrum f
¨
ur Optische Quantentechnologien, Universit
¨
at Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany
2
The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg, Germany
(Received 22 July 2015; published 2 October 2015)
We explore the long-time dynamics of a system of identical charged particles trapped on a closed helix. This
system has recently been found to exhibit an unconventional deformation of the linear spectrum when tuning
the helix radius. Here we show that the same geometrical parameter can affect significantly also the dynamical
behavior of an initially broad excitation for long times. In particular, for small values of the radius, the excitation
disperses into the whole crystal whereas within a specific narrow regime of larger radii the excitation self-focuses,
assuming finally a localized form. Beyond this regime, the excitation defocuses and the dispersion gradually
increases again. We analyze this geometrically controlled nonlinear behavior using an effective discrete nonlinear
Schr
¨
odinger model, which allows us among others to identify a number of breatherlike excitations.
DOI: 10.1103/PhysRevE.92.042905 PACS number(s): 05.45.−a, 37.10.Ty, 37.90.+j, 45.90.+t
I. INTRODUCTION
Whereas the harmonic approximation of interactions pro-
vides valuable information about the stability and the prop-
agation of small amplitude excitations in crystals formed
by interacting particles, their real-time dynamics as well as
their thermal and transport properties are typically subject to
some degree of nonlinearity [1]. Among the most prominent
manifestations of such a nonlinearity are the self-focusing
or self-trapping [2–4] of initial wave packet excitations and
the existence of nonspreading excitations such as breathers
and kinks [5–7]. For discrete systems, a prototype equation
incorporating these features is the s o-called discrete nonlinear
Schr
¨
odinger (DNLS) equation consisting of a linear (disper-
sive) coupling and a cubic nonlinear term [8], used to model
plenty of systems ranging from coupled optical waveguides
[9–11] and Bose-Einstein condensates [12–14] to transport in
DNA molecules [15–17]. The standard spatial arrangement of
sites in most of such one-dimensional (1D) studies is that of a
straight equidistant chain in which the coupling (hopping) is
restricted to nearest neighbors.
Nontrivial lattice geometries for 1D discrete nonlinear
systems have also been studied and have been found to
lead to intriguing new phenomena owing to the interplay
between geometry and nonlinearity. In particular, in curved
1D lattices embedded in a 2D space the bending can act
as a trap of excitations and induce a symmetry breaking of
nontopological solitons [18]. For a 3D space, the helicoidal
lattice structure, such as that of DNA molecules, is found
to enhance the existence and stability of discrete breathers
[19,20]. Furthermore, a curved geometry has been proven
to induce nonlinearity in systems where the underlying
interactions are harmonic [21,22].
In the present work we examine the interplay between
nonlinearity and geometry in a system of identical charged
particles, confined on a curved 1D manifold embedded in
the 3D space, namely a closed (toroidal) helix. In a previous
work [23], we have shown that in such a system, a tuning
of the geometry controlled by the helix radius, leads to an
unconventional deformation of the phononic band structure
including a regime of strong degeneracy. As a consequence,
the propagation of small amplitude localized excitations is
affected significantly and a specific geometry exists at which
the excitations remain localized up to long times. A natural
question therefore arises as to what would be the long time
dynamics of a general excitation and whether there is some
geometrically controllable degree of nonlinearity inherent in
the system, which can alter the propagation characteristics.
We provide an answer to this question by studying the
time propagation of an initially broad excitation on the crystal
of charges. We find that for values of the helix radius far
from the degeneracy regime, the excitation initially spreads
with multiple subsequent revivals due to the closed shape
of the crystal. Within the degeneracy regime, however, the
initial excitation focuses in the course of propagation reaching
finally a rather localized state, serving as a hallmark of the
existing nonlinearity. In order to quantify this nonlinearity we
construct an effective DNLS model with additional nonlocal
nonlinear terms [24,25]. Such a model is found to capture
qualitatively well the localization and dispersion features of the
original dynamics, providing a deeper insight into the observed
effect. Even more, it gives us the opportunity to identify some
discrete breatherlike excitations at the degenerate geometry,
thus adding to the dynamical picture.
The structure of this work is as follows. In Sec. II we de-
scribe our system, commenting also on its linearized behavior.
In Sec. III we present our results for the time evolution of an
excitation in the crystal f or different geometries. In Sec. IV
we construct an effective DNLS model for our system at the
geometries of interest and in Sec. V we use it to identify
some breatherlike excitations. Finally, Sec. VI contains our
conclusions.
II. SETUP AND LINEARIZATION
We consider a system of N identical charges of mass m
0
,
which interact via repulsive Coulomb interactions and are
confined to move on a 1D toroidal helix, parametrized as
r(u) =
⎛
⎝
(
R + r cos(u)
)
cos(au)
(
R + r cos(u)
)
sin(au)
r sin(u)
⎞
⎠
,u∈ [0,2Mπ]. (1)
In Eq. (1) R denotes the major radius of the torus [Fig. 1(a)],
h is the helix pitch, and r the radius of the helix (minor
1539-3755/2015/92(4)/042905(10) 042905-1 ©2015 American Physical Society
A. V. ZAMPETAKI, J. STOCKHOFE, AND P. SCHMELCHER PHYSICAL REVIEW E 92, 042905 (2015)
R
h
r
FIG. 1. (Color online) (a) Equidistant configuration of ions con-
fined on the toroidal helix for ν =
1
2
and N = 8. The yellow arrows
indicate the initial velocities of the particles. (b) Highest (solid blue
line) and lowest (solid red line) frequencies of the linearization
spectrum around the equilibrium equidistant configuration as a
function of r for N = 60 particles. The black dots and the empty
diamonds refer to the frequencies corresponding to the center of
mass and the out of phase mode respectively. The vertical lines mark
the radii of the helix we use in our calculations. The small insets
depict the form of the respective vibrational band structures ω(k)at
the corresponding values of r. (c) Initial local energy E
n
profile as a
function of the particle index n.
radius of the torus), whereas a =
1
M
stands for the inverse
number of windings M =
2πR
h
. The total effective interaction
potential, which results from the constrained motion of the
charges on the helical manifold, reads V (u
1
,u
2
,...u
N
) =
1
2
N
i,j=1,i=j
λ
|r(u
i
)−r(u
j
)|
, where u
j
denotes the coordinate of
particle j and λ is the coupling constant characterizing the
standard Coulomb interactions.
Then the Lagrangian of the system in terms of the u
i
coordinates is given by
L({u
i
,
˙
u
i
}) =
1
2
m
0
2
i=1
|∂
u
i
r(u
i
)|
2
˙
u
2
i
−
1
2
N
i,j=1,i=j
λ
|r(u
i
) − r(u
j
)|
, (2)
where r(u) refers to the parametrization of the toroidal
helix given in Eq. (1). Note that the geometry of the
constraint manifold enters the Lagrangian of Eq. (2) in both
the interaction and the kinetic energy, due to the position-
dependent factor |∂
u
i
r(u
i
)|
2
. If desired, the latter factor can
be removed by transforming to arc-length parametrization
s(u) =
|∂
u
r(u)|du, resulting in the familiar second time
derivative terms in the Euler-Lagrange equations of motion
for the s
i
(t) = s[u
i
(t)] [26], at the cost, however, of losing the
explicit analytical form for the interaction energy.
We choose dimensionless units by scaling all our physical
quantities (e.g., position x, time t, and energy E) with λ, m
0
,
and 2h/π as follows:
˜
x =
xπ
2h
,
˜
t = t
λπ
3
8m
0
h
3
,
˜
E =
2Eh
λπ
,
˜
m
0
= 1,
˜
λ = 1,
omitting in the following the tilde for simplicity.
At commensurate fillings, i.e., M = nN, n = 1,2,...with
the filling factor being ν = 1/n 1, it is found that for
values of the helix radius r up to a critical point r
c
the
ground-state configuration of such a system is the equidistant
polygonic configuration u
(0)
j
= 2(j −1)πn [Fig. 1(a)]. This
configuration loses its stability at r
c
undergoing a zigzag
bifurcation [23].
We focus in this work on the dynamical behavior of charged
particles confined on the toroidal helix in the region r<r
c
,
where the ground state is still the polygonic one. We have
shown in Ref. [23] that in such a region the linear spectrum of
the system changes dramatically with tuning the radius of the
helix r, a fact that crucially affects the propagation of small
amplitude localized excitations. Specifically, it was found that
the width of the linear spectrum decreases as one approaches a
point r
d
of strong degeneracy from below and increases again
beyond that point, while interchanging the character between
the eigenmodes corresponding to the highest and the lowest
frequencies [Fig. 1(b)]. In fact, since the degeneracy is not
complete, it is better to refer to a degeneracy regime within
which the inversion of the spectrum is gradually achieved
while its width remains small [Fig. 1(b) (inset)]. We consider
in this work six different geometries, each corresponding to a
different value of r, covering all the regions with a qualitatively
different linear spectrum [Fig. 1(b)] from the ring limit (r
0
= 0)
to the degeneracy (r
3
,r
4
) and the inversion (r
5
)regime.We
focus on the case of half filling ν =
1
2
for N = 60 particles.
III. TIME PROPAGATION OF A GAUSSIAN EXCITATION
In this section we present and discuss the dynamical
response of our system to an initial excitation. Although the
physical results are in principle independent of the exact
character of this excitation and the means used for its
quantification, the determination of both is essential for the
illustration and the theoretical description of our findings.
Dealing with classical systems and seeking an excitation
measure whose total amount is conserved in time, the natural
choice is a (to be defined) energy distribution associated
with each particle, referred to hereafter as local energy E
n
.
Whereas the kinetic energy K consists of parts allocated
to each individual particle, the potential energy cannot be
uniquely partitioned, yielding different definitions of local
energy [27–29]. Aiming for them to be strictly positive for
all possible excitations (a considerably nontrivial requirement
for systems with Coulomb interactions), we define our local
energies E
n
in a rather unconventional way, focusing on a
042905-2
DYNAMICS OF NONLINEAR EXCITATIONS OF . . . PHYSICAL REVIEW E 92, 042905 (2015)
FIG. 2. (Color online) (a)–(f) Time evolution of the initial Gaussian excitation presented in Fig. 1(c) for N = 60, ν = 1/2 for increasing r
corresponding to the points (a) r
0
,(b)r
1
,(c)r
2
,(d)r
3
,(e)r
4
, and (f) r
5
marked in Fig. 1(b). Colors encode the values of local energy E
n
for each
particle n and time t.Forthesamevaluesofr panels (g)–(l) depict the time evolution of the normalized participation ratio P of the excitation.
The solid blue lines are the results from our numerical simulations corresponding to (a)–(f), whereas the dashed red lines correspond to the
results for a harmonic approximation of the potential.
positive decomposition of the harmonic interaction term [29].
Our complete definition and a more detailed discussion of local
energies are provided in the Appendix.
We start with an initially broad excitation of a Gaussian
profile in terms of local energies [Fig. 1(c)]. Since the
local energies E
n
depend trivially on the particles’ velocities
(contrary to what is the case for the particles’ positions), the
most straightforward way to obtain such a Gaussian local
energy profile is by exciting the particles with a suitable
velocity distribution. Of course, for a given local energy profile
the magnitude of such a velocity distribution can be uniquely
determined, but there is a freedom in the direction of the
velocity for each particle. We choose here all the velocities
to point in the positive direction [Fig. 1(a)].
The ensuing dynamics is shown in Figs. 2(a)–2(f).Obvi-
ously, the time evolution of the Gaussian excitation possesses
a drastic dependence on the geometry, controlled by the
helix radius r. In particular, for r<r
d
[Figs. 2(a), 2(b)]the
excitation spreads into the whole crystal and refocuses almost
periodically at the time instants when the left and the right
propagating parts of the excitation meet and superimpose at the
diametrically opposite point of the closed helix. As discussed
in Ref. [23] the spreading velocity decreases as r is increased,
following the width of the linear spectrum [Fig. 1(b)]. As the
width becomes smaller the features of the time evolution alter
significantly. Already at the point r
2
the excitation does not
spread any more into the crystal, but it alternately focuses
and again defocuses to its initial shape [Fig. 2(c)]. Even more
surprisingly, within the degeneracy region [Figs. 2(d), 2(e)]
the initial excitation undergoes a focusing after some time
scale t
F
[t
F
≈ 4000 for Fig. 2(d), t
F
≈ 6000 for Fig. 2(e)].
Subsequently, the wave packet loses its smooth envelope
and fragments into a number of highly localized excitations.
Depending on r, the routes towards such a localized state
can be different, with the wave packet evolving initially one
central peak [Fig. 2(d)] or two side peaks [Fig. 2(e)]. Another
interesting feature within the degeneracy regime is that the
reflection symmetry of the initial excitation profile can break
042905-3
A. V. ZAMPETAKI, J. STOCKHOFE, AND P. SCHMELCHER PHYSICAL REVIEW E 92, 042905 (2015)
in the course of propagation, attaining after some time a
significantly asymmetric form [Fig. 2(d)]. The direction of the
asymmetry depends on the direction of the initial velocities of
the particles. Beyond the degeneracy [Fig. 2(f)] the spreading
of the excitation into the crystal reappears with a periodic
refocusing but the propagation pattern is much different owing
to the inverted form of the vibrational band structure.
In order to quantify the degree of focusing or localization
of the excitation, we examine the time dependence of the
normalized participation ratio
P =
1
N
N
n=1
E
n
2
N
n=1
E
2
n
. (3)
Evidently, this quantity can take values between 1/N and
1, with P = 1 signifying the case of a completely extended
excitation where the energy is equipartitioned between all the
particles and P = 1/N marking the opposite case of a fully
localized excitation in a single particle. Note that the local
energies E
n
in (3) should be non-negative for the definition to
make sense.
Our results, presented in Figs. 2(g)–2(l), support our
discussion above. Especially the focusing of the excitation
after t
F
is evident in Figs. 2(j)–2(k). However, the subsequent
drop in P is much stronger in Fig. 2(j) than in Fig. 2(k),in
line with the observation that at r = r
3
the final localized state
consists of less excited particles [Fig. 2(d)] than at r = r
4
[Fig. 2(e)].
In Ref. [23] it was demonstrated that a small localized
initial excitation does not spread significantly for short times
at the degeneracy point, in contrast to the behavior for other
geometries. This fact was understood solely by an inspection
of the linearization spectrum. Here, however, the situation is
different. Not only t he complete absence of spreading, but
especially the existence of self-focusing calls for an account
of the underlying nonlinearity. This is further emphasized and
supported by Figs. 2(g)–2(l) where the results of the prop-
agation within the harmonic approximation of the potential
are also displayed. As long as the total amplitude of the initial
excitation is small enough, the harmonic approximation works
well. As the amplitude is increased this approximation will
start to fail, and nonlinear effects are expected to show up.
The results of the present work suggest that except for the
amplitude, also the geometry allows to control the importance
of the nonlinearity. Specifically, for the given amplitude and for
geometries far from the degeneracy regime [Figs. 2(g), 2(h),
2(l)] the harmonic approximation qualitatively reproduces the
exact time evolution of the participation ratio, although, as
should be expected, there are quantitative deviations. In con-
trast, close to and within the degeneracy regime the harmonic
approximation fails completely, predicting a spreading and an
extended form of the excitation, instead of localization. This
makes it clear that regarding the focusing, we indeed encounter
a nonlinear phenomenon.
Complementary information about the spreading or
localization of the Gaussian excitation can be obtained from
its time evolution in the reciprocal space, i.e.. in terms of
the wave numbers k ∈ [−
π
s
,
π
s
], where s denotes the arc
length interparticle distance of the equidistant ground-state
configuration. To this extent, we examine how the discrete
FIG. 3. (Color online) Absolute value of the discrete Fourier
transform of the local energy excitation profile as a function of time
and wave number k. (a)–(f) correspond to increasing r, ranging from
r
0
to r
5
as indicated in Fig. 1(b).
Fourier transform of the local energy profile evolves in the
reciprocal space [Figs. 3(a)–3(f)]. The initial excitation, being
overall extended in the coordinate space [Fig. 1(c)], appears
to be rather localized around the wave vector k = 0inthek
space. It remains l ocalized as expected for r<r
d
[Figs. 3(a),
3(b)], but when approaching r
d
more wave vectors in the
vicinity of k = 0 become excited [Fig. 3(c)]. At the degenerate
geometries [Figs. 3(d), 3(e) ] the excitation expands rapidly
in the reciprocal space, populating after the characteristic
time t
F
almost all the wave numbers k. The space localized
solution consists therefore of most of the k modes, with the
population of the initial k = 0 mode being dominant. At
r>r
d
[Fig. 3(f)] the excitation remains, as for small r,inthe
narrow vicinity of the k = 0 mode.
Before proceeding with our study of the nonlinear behavior,
let us note that for the results presented here we have used a
rather small initial excitation with a total energy of the order of
1% of the ground-state energy per particle E
GS
/N. For larger
amplitudes the self-focusing can occur also for smaller radii,
i.e., at a greater distance from the degeneracy region.
IV. EFFECTIVE NONLINEAR MODEL
The dynamics analyzed in the previous section is charac-
terized by a self-focusing process of excitations for the case
042905-4
DYNAMICS OF NONLINEAR EXCITATIONS OF . . . PHYSICAL REVIEW E 92, 042905 (2015)
of degenerate geometries and therefore suggests a prominent
role of the nonlinearity. We aim in this section to identify
and quantify the leading nonlinear terms as well as to derive
a DNLS effective model in the region of degeneracy. As a
result we will, among others, gain insight into the excitation
amplitudes and the time scale t
F
for l ocalization.
A. Dominant nonlinear terms
Since the initial excitation is small enough we attempt to
identify the dominant anharmonic terms by expanding the
potential around the equilibrium configuration {u
(0)
} up to
fourth order. It is advantageous to do so in the arc length
parametrization s so that the final Euler-Lagrange equations
of motion and particularly the kinetic terms would assume the
standard form [ 26]. To this purpose we calculate the matrices
H
ij
=
∂
2
V
∂s
i
∂s
j
{s
(0)
}
,G
ij
=
∂
3
V
∂s
2
i
∂s
j
{s
(0)
}
(4)
corresponding to the Hessian and the matrix of the third
derivatives respectively. Terms involving more than two
different positions are equal to zero, since the total potential
V is a sum of exclusively two-body potential terms.
For the fourth derivative terms we define two further
matrices as
M
ij
=
∂
4
V
∂s
2
i
∂s
2
j
{s
(0)
}
,Q
ij
=
∂
4
V
∂s
3
i
∂s
j
{s
(0)
}
(5)
for i = j and
M
ii
= Q
ii
=
1
2
∂
4
V
∂s
4
i
{s
(0)
}
,
separating the derivatives of the same order in s
i
, s
j
from those
with a different order and splitting the diagonal derivative
terms (involving differentiation with respect to a single
position) in two. The calculation of the above derivatives in
the arc length parametrization can be carried out by using the
respective derivatives in the u coordinate space as well as the
known relations for derivatives of inverse functions. Under
these considerations and denoting s
j
− s
(0)
j
= x
j
the potential
reads
V ≈E
GS
+
1
2
N
i,j=1
x
i
H
ij
x
j
+
1
6
N
i,j=1
x
2
i
G
ij
x
j
+
1
24
N
i,j=1
x
2
i
M
ij
x
2
j
+
1
24
N
i,j=1
x
3
i
Q
ij
x
j
(6)
leading to the equations of motion
¨
x
n
=−
j=n
H
nj
x
j
− H
nn
x
n
−
1
3
j=n
x
n
G
nj
x
j
−
1
6
j=n
G
jn
x
2
j
−
1
6
j=n
x
n
M
nj
x
2
j
−
1
3
M
nn
x
3
n
−
1
8
j=n
x
2
n
Q
nj
x
j
−
1
24
j=n
Q
jn
x
3
j
, (7)
which consist of harmonic, quadratic, and third-order nonlin-
ear terms. All the matrices appearing in these equations are
symmetric except for the matrix G
ij
, relating to the quadratic
force terms, which is fully antisymmetric. This can affect
significantly the symmetry of the expected solutions. If the
quadratic nonlinear force terms are ignored then the equations
of motion (7) possess the symmetry x
n
→−x
n
which also
results in symmetric excitations keeping their symmetry in the
course of propagation. The quadratic force terms, however,
break the reflection s ymmetry and allow for an asymmetric
evolution of initially symmetric excitations as the one observed
in Fig. 2(d) .
Let us mention at this point that in the ring limit r = r
0
= 0,
where a separation of the center of mass holds, the matrices
involved in Eq. (7) are not independent but obey the relations
H
nn
=−
j = n
H
nj
,M
nn
=−
1
8
j = n
Q
nj
,M
nj
=−
3
4
Q
nj
yielding
¨
x
n
=−
j=n
H
nj
(x
j
− x
n
) −
1
6
j=n
G
jn
(x
j
− x
n
)
2
−
1
24
j=n
Q
nj
(x
j
− x
n
)
3
, (8)
which for only nearest-neighbor (NN) couplings leads to a
Fermi-Pasta-Ulam kind of equations of motion [30], with both
quadratic and cubic nonlinear interactions.
Before proceeding, we note that the matrix elements
H
ij
,G
ij
,M
ij
,Q
ij
depend, due to the symmetry of the ground-
state configuration, only on the index difference m = i − j.
Therefore, when referring to these elements in the following
we will use the notation H
m
,G
m
,M
m
,Q
m
.
B. DNLS model
The effect of localization of t he initial wave packet
occurs, as we have observed in Sec. III, in the regime close
to degeneracy. There, the diagonal terms of the Hessian
provide the dominant contribution to the linear spectrum as
the off-diagonal ones are very small. Since |H
0
|2|H
1
|
[Fig. 4(a)], one can use in that regime the so-called rotating
wave approximation (RWA) [31], assuming that the position
coordinate can be described as
x
n
(t) =
n
(t)e
−iω
0
t
+
∗
n
(t)e
iω
0
t
, (9)
with
n
(t) a slowly varying amplitude and ω
2
0
= H
0
yielding
a fast oscillating phase e
±iω
0
t
. Apart from the requirement
of a weak dispersion, a condition for a sufficiently weak
nonlinearity has also to be satisfied [31], namely |H
0
|
1
3
|M
0
|max[x
i
(0)]
2
with max[x
i
(0)] the maximum initial dis-
placement (or respectively momentum) of a single particle.
This criterion is satisfied as well in our case, since the initial
conditions we have used lead to |H
0
| of about 100 times larger
than
1
3
|M
0
|max[x
i
(0)]
2
.
Using the ansatz (9) and making the assumption |
d
n
dt
|
ω
0
|
n
|, as well as neglecting the rapidly oscillating terms with
042905-5