# Chaos and Intermittency in the DNLS Equation Describing the Parallel Alfvén Wave Propagation

Authors: Gustavo Krause; Sergio Elaskar; Andrea Costa
Journal: Journal of Astrophysics
 (2014)

---

## Abstract

When the Hall effect is included in the magnetohydrodynamics equations (Hall-MHD model) the wave propagation modes become coupled, but for propagation parallel to the ambient magnetic field the Alfvén mode decouples from the magnetosonic ones, resulting in circularly polarized waves that are described by the derivative nonlinear Schrödinger (DNLS) equation. In this paper, the DNLS equation is numerically solved using spectral methods for the spatial derivatives and a fourth order Runge-Kutta scheme for time integration. Firstly, the nondiffusive DNLS equation is considered to test the validity of the method by verifying the analytical condition of modulational stability. Later, diffusive and excitatory effects are incorporated to compare the numerical results with those obtained by a three-wave truncation model. The results show that different types of attractors can exist depending on the diffusion level: for relatively large damping, there are fixed points for which the truncation model is a good approximation; for low damping, chaotic solutions appear and the three-wave truncation model fails due to the emergence of new nonnegligible modes.


---

## Body

## 1. Introduction

Alfvén waves are one of the most characteristic features of magnetized laboratory and space plasmas. They are driven by different sources, for example, nonuniform plasma parameters, beams of charged particles, and electrostatic and electromagnetic waves [1]. In space plasmas, large MHD amplitude fluctuations with typical proton cyclotron local frequencies were detected on the Earth magnetosphere [2]. Also, nonlinear Alfvén waves were extensively detected in the solar wind [3] and they are believed to be responsible for the turbulent heating of stellar coronas [4]. The comprehension of nonlinear properties of dispersive Alfvén waves is of crucial importance to interpret the abundant amount of low frequency data provided by space plasma observations.On the other hand, the interaction of spatial tethers with the Earth ionosphere and the ambient magnetic field leads to the emission of Alfvén waves forming structures called Alfvén wings. This phenomenon could be applied to produce electric power, generate artificial auroras [5], or to moderate spatial trash [6]. Far from the tethers, a linear analysis could be appropriate [7], but near the conductor intense waves with important nonlinear effects are expected.For the study of the plasma behavior the magnetohydrodynamics (MHD) equations are usually used, but when the frequencies of interest are of the order of the ion-cyclotron frequency or when the characteristic longitudes are comparable with the ion inertial length, these equations have to be extended to include the effect of a finite ion-cyclotron frequency, which is referred to as the “Hall effect”. In the resulting Hall-MHD model, although the magnetosonic and the Alfvén modes can still be identified, they are coupled leading to a dispersive evolution. This allows that, for a homogeneous plasma, nonlinear high-frequency Alfvén waves can propagate parallel to a uniform magnetic field since the nonlinearities are balanced by the dispersive term which is constituted by the Hall effect in the induction equation.Using a two-fluid, quasineutral approximation with electron inertia and current displacement neglected, the MHD-Hall equations normalized by density and magnetic field reference values,B
0 and ρ
0, are [8]
(1)
∂
ρ
∂
t
+
∇
·
(
ρ
u
)
=
0
,
ρ
(
∂
∂
t
+
u
·
∇
)
u
=
-
∇
p
+
V
A
2
(
∇
×
B
)
×
B
,
∂
B
∂
t
=
∇
×
(
u
×
B
)
-
V
A
2
Ω
i
∇
×
[
1
ρ
(
∇
×
B
)
×
B
]
,
where ρ is the plasma density, u the velocity, B the magnetic field, Ω
i the ion-cyclotron frequency, and V
A
=
B
0
/
(
ρ
0
μ
0
)
1
/
2 is the Alfvén speed where μ
0 is the vacuum permeability. The equation of state obeying the polytropic law p
∝
ρ
-
γ, with γ the polytropic coefficient, completes the set of equations.Weak nonlinearities were studied by perturbation theory [8] leading to a Korteweg-de Vries (KdV) equation for the magnetosonic modes and a modified KdV equation for the Alfvén mode when the propagation angle is large enough [9–11]. When propagation parallel (or quasiparallel) to the ambient magnetic field is considered, the MHD waves are degenerated; the Alfvén mode is decoupled and the waves are circularly polarized, being described by the derivative nonlinear Schrödinger (DNLS) equation (e.g., [12–15]) although the case of arbitrary propagation angles is also valid for a high-β plasma [16].The DNLS equation was derived by several authors using different techniques: through the Vlasov equation [12] and using the Hall-MHD equation (e.g., [13–15]). The nondimensional DNLS equation in a frame moving with the Alfvén speed along the magnetic field direction (z) with V
A
≫
c
s (β
=
c
s
/
V
A
≈
0, with c
s the sound speed) results in
(2)
∂
b
∂
t
+
∂
∂
z
(
|
b
|
2
b
)
+
i
∂
2
b
∂
z
2
+
γ
^
b
=
0
,
where the positive sign in the dispersive term indicates that left-hand polarized waves are considered, γ
^ is an appropriate damping/driving linear operator, and the dimensionless variables b, z, and t are defined by
(3)
b
=
B
x
+
i
B
y
2
B
0
,
2
Ω
i
t
⟶
t
,
2
Ω
i
V
A
z
⟶
z
.The advantage of using the DNLS is due to the fact that this equation “belongs” to soliton theory; thus much about its solutions is known although it is nonlinear. On the other hand, soliton theory was used to explain many observations of quasiparallel finite-amplitude MHD waves where the competition between nonlinear steepening and dispersion occurs, such as in the foreshock regions of quasiparallel shocks occurring in connection with the bow shocks of planets or comets [17]. Therefore, the solutions of the DNLS equation can be attributed to specific physical problems [18].Many exact solutions of the DNLS equations are known and stability analysis has been made for some of them [19]. Among these solutions, the circularly polarized wave is of particular interest because the stability of this solution is easily investigated (see Section 3).On the other hand, different analytical approaches were used, which consist of reducing the DNLS equation to a system with low dimensions, either considering stationary waves or truncating the system by including only a finite number of modes. In the first case, the DNLS equation is reduced to a set of three ordinary differential equations where the free variables are the two components of the transverse magnetic field and the phase wave [20]. In this way, a continuous three-dimensional dynamical system is obtained, which would allow retaining the nonlinear evolution of driven conservative and dissipative Alfvén waves that is registered in more complicated high-dimensional models [21]. This approach was extensively used to study the Alfvén intermittent turbulence in space plasmas [22, 23], showing that the onset of Alfvén turbulence can occur via a crisis-induced intermittency.The other analytical approach, consisting in the truncation of the system, was also widely studied. This approximation supposes that the solution is the sum of a finite number of modes, by which a set of ordinary differential equations is obtained, where the order of the system depends on the number of modes used for the truncation. Based on numerical studies where the DNLS was fully integrated, three-wave truncation models with a resonance relation2
k
0
=
k
1
+
k
2 (k
0 is the mother wave) have been proposed. Although the results of the first analysis of this model did not exhibit a chaotic dynamics [21], subsequent works showed that the three-wave truncation model does present chaos through different routes [24–27]. Similar studies were performed in the context of the nonlinear Schrödinger equation (NLS) [28].Numerical studies of the DNLS have been carried out by numerous authors exploring the integrability properties of this equation [21, 27, 29–31]. Most of these works are based on spectral methods assuming periodic boundary conditions, but different integration schemes have been used to solve stability problems.In this paper, a numerical analysis of the driven dissipative DNLS equation for the particular case of three resonant Alfvén waves is presented. The numerical results are compared with those obtained by a three-wave truncation model which was carried out to represent Alfvén wave fronts generated by orbiting conductive tethers interacting with the ambient magnetic field in the ionosphere [5, 26, 27]. This analysis allows establishing the application range of the truncation method by inspecting the power spectrum, in a similar way as the work by Ghosh and Papadopoulos [21]. However, in this case a simpler numerical scheme is used with an aliasing filter that only cancels 1
/
3 of the wavenumber domain, allowing retaining a greater number of effective modes [32]. In addition, new features of the dynamic solutions are obtained, showing that there are more configurations which present bifurcation diagrams exhibiting a Hopf bifurcation that separates the more diffusive region, with stationary solutions, from the less diffusive one, with complex dynamic solutions, where chaotic attractors are developed by different routes (Feigenbaum cascades, intermittency, and crises). Also, new coherence relations are found for the stationary solutions, where the existence of five resonant modes is detected. In this way, the analysis of the present paper extends the results of previous works [21, 27]. On the other hand, another new aspect of this research is the detailed analysis of the stability of the numerical scheme presented here, which not only allows verifying the code, but also permits to study the evolution of the conservative DNLS equation under unstable conditions. These results show that the presence of periodic or chaotic evolutions depends on the parameters of the initial wave.The organization of the paper is as follows. In Section2, a brief description of the numerical method is presented. In Section 3, the numerical method is tested by verifying the analytical conditions of modulational stability in the nondiffusive case. In Section 4, numerical results for the driven diffusive DNLS equation under a three-wave resonant initial condition are shown, comparing the solutions with those of the three-wave truncation model. Finally, in Section 5 the main conclusions of the work are presented.


## 2. The Method

The DNLS equation, as type-KdV equations, admits that periodic boundary conditions are used when the initial conditions are periodic [33]. Therefore, boundary conditions of the form
(4)
b
(
z
,
t
)
=
b
(
z
+
L
,
t
)
are appropriate to implement numerical simulations [19], where L is the domain size.Spectral methods based on trigonometric functions automatically and individually satisfy the periodic boundary conditions. The discrete Fourier expansion allows computing the Fourier coefficients of a functionb
(
z
) defined in the domain (
-
L
/
2
,
L
/
2
] and known at points
(5)
z
j
=
-
L
2
+
L
N
j
,
j
=
1
,
…
,
N
,
where N is the number of grid points. The discrete Fourier transform and its inverse are
(6)
ℱ
[
b
j
]
=
b
¯
k
=
1
N
∑
j
=
1
N
‍
b
j
exp
⁡
(
-
i
2
π
n
L
z
j
)
,
ℱ
-
1
[
b
¯
k
]
=
b
j
=
∑
n
=
-
(
N
/
2
)
+
1
N
/
2
‍
b
¯
k
exp
⁡
(
-
i
2
π
n
L
z
j
)
,
where the notations b
j
=
b
(
z
j
) and b
¯
k
=
b
¯
(
k
) are used. The wavenumbers k are given by
(7)
k
=
2
π
n
L
,
n
=
-
N
2
+
1
,
…
,
N
2
.Equations (6) define an exact transform pair between the N grid values b
j and the N discrete Fourier coefficients b
¯
k. The transformation can be performed from one representation to the other without loss of information [34]. Note that from (6) the derivatives at the grid points are given by
(8)
d
n
b
j
d
z
n
=
ℱ
-
1
[
(
i
k
)
n
b
¯
k
]
,
and then we can write the DNLS equation as
(9)
∂
b
∂
t
=
ℱ
-
1
[
-
i
k
ℱ
[
|
b
|
2
b
]
+
i
k
2
ℱ
[
b
]
]
-
γ
^
b
;
thus, the time derivative can be computed for a defined damping/driving operator γ
^. In this way, the time integration can be explicitly performed using a fourth order Runge-Kutta scheme.It should be noted that theN truncation of the series introduces the “aliasing” error, which is associated with the difference between the discrete and the continuous Fourier expansion. This error can cause numerical instabilities in the time integration of nonlinear equations [35]. The usual strategy employed to avoid these effects is to remove the “aliased” modes by applying an “all-or-nothing” filter where the Fourier coefficients external to a certain portion of the central domain are canceled. In this paper, contrary to previous works where half of the modes are preserved [21, 27], the 2
/
3 of the central domain is used, allowing to extend the amount of effective modes used in the simulations [32]. The subsequent results show that this is a good implementation.


### 2.1. Modeling of the Diffusive Term

In (2), γ
^ is used to represent an appropriate damping/driving linear operator, which is assumed to have the characteristic diffusive form
(10)
γ
^
=
-
η
∂
2
∂
z
2
,
where η is the damping coefficient. Taking into account this definition and (8), it is obtained for the generic sample point z
j that
(11)
γ
^
j
=
-
η
k
j
2
.This operator is equivalent to consider a resistive damping model, where the dissipation is proportional tok
j
2 [24]:
(12)
|
γ
^
j
|
≈
V
A
2
k
j
2
2
Ω
e
Ω
i
τ
c
=
η
k
j
2
,
where Ω
e is the electron-cyclotron frequency and τ
c is the characteristic Braginskii collision time.To represent the excitation of the mother wave (identified withk
0), the Dirac function δ is used:
(13)
γ
^
j
=
γ
g
δ
(
k
j
-
k
0
)
-
η
k
j
2
,
where γ
g gives the excitation level of the wave. For k
j other than k
0, the Dirac function is zero and (13) reduces to (10).With the latest definitions, the real and the imaginary parts of (9) are
(14)
Re
[
∂
b
j
∂
t
]
=
ℱ
-
1
[
k
j
Im
⁡
[
ℱ
[
|
b
j
|
2
b
j
]
]
-
k
j
2
Im
⁡
[
b
¯
k
j
]
]
+
ℱ
-
1
[
(
γ
g
δ
(
k
j
-
k
0
)
-
η
k
j
2
)
Re
[
b
¯
k
j
]
]
,
Im
⁡
[
∂
b
j
∂
t
]
=
ℱ
-
1
[
-
k
j
Re
[
ℱ
[
|
b
j
|
2
b
j
]
]
+
k
j
2
Re
[
b
¯
k
j
]
]
+
ℱ
-
1
[
(
γ
g
δ
(
k
j
-
k
0
)
-
η
k
j
2
)
Im
⁡
[
b
¯
k
j
]
]
,
and the time derivative of the DNLS equation is finally obtained.


## 2.1. Modeling of the Diffusive Term

In (2), γ
^ is used to represent an appropriate damping/driving linear operator, which is assumed to have the characteristic diffusive form
(10)
γ
^
=
-
η
∂
2
∂
z
2
,
where η is the damping coefficient. Taking into account this definition and (8), it is obtained for the generic sample point z
j that
(11)
γ
^
j
=
-
η
k
j
2
.This operator is equivalent to consider a resistive damping model, where the dissipation is proportional tok
j
2 [24]:
(12)
|
γ
^
j
|
≈
V
A
2
k
j
2
2
Ω
e
Ω
i
τ
c
=
η
k
j
2
,
where Ω
e is the electron-cyclotron frequency and τ
c is the characteristic Braginskii collision time.To represent the excitation of the mother wave (identified withk
0), the Dirac function δ is used:
(13)
γ
^
j
=
γ
g
δ
(
k
j
-
k
0
)
-
η
k
j
2
,
where γ
g gives the excitation level of the wave. For k
j other than k
0, the Dirac function is zero and (13) reduces to (10).With the latest definitions, the real and the imaginary parts of (9) are
(14)
Re
[
∂
b
j
∂
t
]
=
ℱ
-
1
[
k
j
Im
⁡
[
ℱ
[
|
b
j
|
2
b
j
]
]
-
k
j
2
Im
⁡
[
b
¯
k
j
]
]
+
ℱ
-
1
[
(
γ
g
δ
(
k
j
-
k
0
)
-
η
k
j
2
)
Re
[
b
¯
k
j
]
]
,
Im
⁡
[
∂
b
j
∂
t
]
=
ℱ
-
1
[
-
k
j
Re
[
ℱ
[
|
b
j
|
2
b
j
]
]
+
k
j
2
Re
[
b
¯
k
j
]
]
+
ℱ
-
1
[
(
γ
g
δ
(
k
j
-
k
0
)
-
η
k
j
2
)
Im
⁡
[
b
¯
k
j
]
]
,
and the time derivative of the DNLS equation is finally obtained.


## 3. Verification of the Numerical Method

In order to test the numerical scheme, the nondiffusive DNLS equation (γ
^
=
0) is analyzed to verify the accomplishment of the analytical conditions of modulational stability, which can be established for an initial magnetic field condition
(15)
b
(
z
,
0
)
=
b
0
=
A
0
exp
⁡
(
i
k
0
z
)
,
where the initial wavenumber is
(16)
k
0
=
2
π
n
0
L
,
n
0 being an integer, L the domain size, and A
0 the initial amplitude.The DNLS modulational instability analysis is performed using the modulation of the constant amplitude wave train expressed by(17)
b
(
z
,
t
)
=
A
0
exp
⁡
[
i
(
k
0
z
-
ω
t
)
]
[
1
+
ɛ
(
z
,
t
)
]
,
where b and ɛ are complex periodic functions in the range [
-
L
/
2
,
L
/
2
], |
ɛ
|
2
≪
1 and ω
=
A
0
2
k
0
-
k
0
2. According to the positive sign in the dispersive term of (2), only the left-hand polarized waves are considered (k
0
>
0). On the other hand, the right-hand polarized waves (k
0
<
0) are unconditionally stable for parallel modulation [30]. These conditions are due to the left-hand polarized waves leading to the resonance of the ion cyclotron mode at higher values of k, while the right-hand polarized branch connects to the whistler mode without resonance [18].Expandingɛ
(
z
,
t
) in a Fourier base gives
(18)
ε
(
z
,
t
)
=
∑
j
=
-
∞
∞
‍
ε
^
j
(
t
)
exp
⁡
(
i
ℓ
j
z
)
,
ε
^
j
(
t
)
=
ε
^
j
(
0
)
exp
⁡
(
i
λ
j
t
)
,
where ℓ
j
=
2
π
j
/
L is the modulational wavenumber. Replacing this expression in the DNLS equation reads [36]
(19)
λ
j
=
i
2
(
A
0
2
-
k
0
)
ℓ
j
±
|
ℓ
j
|
(
2
k
0
-
A
0
2
)
A
0
2
-
ℓ
j
2
.From (19), the instability conditions for the left-hand polarized waves are obtained:
(20)
2
k
0
<
A
0
2
⟶
Marginally
stable
,
2
k
0
>
A
0
2
⟶
Unstable
.To numerically verify this condition the parameterE
k
0 is defined, which represents the ratio between the energy carried by the initial wave and the total energy of the system in the Fourier space:
(21)
E
k
0
=
|
b
¯
k
0
|
∑
j
=
1
N
|
b
¯
k
j
|
,
where b
¯
k
0 is the discrete Fourier transform of the initial wave, that is, the amplitude of wave k
0, N is the total number of modes used in the simulation, and k
j indicates the wavenumber. E
k
0
=
1 implies that the energy is concentrated in the initial wave. When the instability is triggered, the energy initially concentrated in k
0 is transferred to other modes and E
k
0
<
1. Figure 1 shows, for different initial wavenumbers, the triggering time of the instability as a function of the initial amplitude A
0.Instability time as a function of the initial amplitudeA
0 for different initial wavenumbers k
0
=
2
π
n
0
/
L. Dashed lines indicate the stability limits for each configuration, (20). The stable region is bounded by the dashed line and A
0
→
∞ (for clarity, the figure is divided into even and odd n
0 cases).
(a)
(b)Figure1 indicates that the instability condition A
0
<
2
k
0 is always verified. For larger amplitudes, the stability was confirmed for initial amplitude values twice the limit amplitude and a simulation time of 50000 units. The time step used in the simulations is τ
=
1
0
-
3, with a grid of N
=
256 points. These values were established searching for the minimum computational time for which conditions (20) and the stability of the right-hand polarized waves were verified.This analysis not only allows verifying the stability of the numerical scheme but also permits to study the evolution of the unstable solutions, which is not possible by the theoretical analysis. Figure2 shows the parameter E
k
0 as a function of time for different initial wavenumbers k
0 and amplitudes A
0. In addition, the energy distribution is also shown in order to know how the initial energy (curve |
b
¯
k
|
0) is transferred when E
k
0 reaches its minimum value (curve |
b
¯
k
|
min
⁡).Left panels in (a) and (b):E
k
0 evolution as a function of time for different initial wave amplitudes A
0. Right panels in (a) and (b): energy distribution among the different modes for the initial condition |
b
-
k
|
0 and the corresponding to the minimum E
k
0 value (|
b
-
k
|
min
⁡).
(a)
n
0
=
5
(b)
n
0
=
20It can be seen that the evolution ofE
k
0 is similar in all cases independent of the initial wavenumber k
0. This evolution can be divided in two different behaviors depending on the initial amplitude A
0. For small initial amplitudes, quasiperiodic solutions are registered in which the initial wave recovers its energy and E
k
0 has a cyclic evolution between a maximum (E
k
0
≈
1) and a minimum value. This minimum value is associated with the amplitude A
0, being smaller with shorter periods for larger amplitudes. Concerning the energy transfer, it is observed that it occurs mainly among few modes, the initial wave (mother wave k
0) and the two adjacent modes (daughter waves) with wavenumbers k
1 and k
2 that satisfy the resonance relation 2
k
0
=
k
1
+
k
2.For large amplitudes, the regular solution is lost and a spatio-temporal chaotic evolution is produced, which is independent of the initial wavenumber. Unlike the previous case, the initial wave does not recover its initial energy and a random distribution among the different modes occurs.The transition between the quasiperiodic and the chaotic behavior occurs forA
0 values where the results are a combination of both solutions. In the transition range, the initial wave partially recovers the initial energy with a cyclic performance as in the first case, but a descendant behavior of the maximum energy occurs until the irregular trend is established.The initial amplitude values associated with the different behaviors depend on the initial wavenumber. From the figures, it is noted that higher initial modes have periodic behavior ranges smaller with bounded transition intervals. On the contrary, low initial modes have periodic behavior at higher amplitudes and the transition is smoother.


## 4. Results of the DNLS Equation with Diffusive and Excitatory Effects

In this section the results for the driven diffusive DNLS equation are presented, where a three-wave initial condition under a resonant relation is used. The wavek
0 (mother wave) is linearly unstable and the others are damped. Such simplified system was extensively studied since it can be useful in making a preliminary analysis of Alfvén waves generated by tethers [24, 26, 27]. On the other hand, the system driven by a one unstable mode in the Fourier spectrum was also widely used as a physical analogy of a warm-low-density field-aligned ion beam creating growth along a finite bandwidth of modes through the electromagnetic ion-beam cyclotron instability [21].The three-wave initial condition is(22)
b
(
z
,
0
)
=
∑
j
=
0
2
‍
a
j
exp
⁡
(
i
k
j
z
)
,
with wavenumbers k
1
=
(
1
-
Δ
)
k
0 and k
2
=
(
1
+
Δ
)
k
0 that satisfy the resonance relation 2
k
0
=
k
1
+
k
2 and the periodic boundary conditions, (4).The numerical results are compared with a three-wave truncation model under resonant interaction. The truncation model is obtained looking for solutions of the DNLS in the form(23)
b
(
z
,
t
)
=
∑
j
=
0
2
‍
a
j
exp
⁡
[
i
(
k
j
z
-
ω
j
t
+
ϕ
j
)
]
,
involving three modes satisfying the resonance relation and the dispersion relation ω
j
=
∓
k
j
2. Substituting (23) in the DNLS and only considering terms involving k
0, k
1, and k
2, the truncation equations can be derived to obtain the amplitudes a
0, a
1, and a
2 and the relative phase θ
=
ν
t
+
ϕ
1
+
ϕ
2
-
2
ϕ
0 with ν
=
2
ω
0
-
ω
1
-
ω
2
=
±
2
Δ
2, a frequency mismatch. The upper (lower) sign corresponds to the left-hand (right-hand) polarized waves [26, 27].The comparison between both methods is done considering the amplitudes of the resonant wavesa
0, a
1, and a
2 and the energy in the whole domain
(24)
E
m
=
∫
-
L
/
2
L
/
2
‍
|
b
|
2
d
z
.
In the case of the truncation model, where the only non-zero amplitudes are a
0, a
1, and a
2, the last equation results in:
(25)
E
m
=
∫
-
L
/
2
L
/
2
‍
|
b
|
2
d
z
=
(
a
0
2
+
a
1
2
+
a
2
2
)
L
.The analysis is carried out considering attractors with a mother wave with wavenumber given byn
0
=
9 and an excitation γ
^
0
=
0.02. The attractors are numerically obtained by the simulation for certain damping coefficient values depending on the wavenumbers of the daughter waves. It should be noted that in this analysis only long term dynamic solutions are considered; therefore, transient solutions as those due to chaotic saddles are not presented. However, it is mentioned that transient chaos appears in almost all the studied configurations, either in the case of stationary solutions or before the convergence to limit cycles. A very detailed study of this type of nonattracting chaotic sets can be found at [22, 23, 28].To clarify the analysis, the results are separated into stationary solutions and dynamic solutions as follows.


### 4.1. Stationary Solutions

When intermediate and strong diffusion levels (η
≳
0.05) are considered, a set of attractors is obtained which consist of stationary traveling waves with most of the energy concentrated in the three initial modes, such that its amplitudes and the energy can be considered as stable fixed points. Table 1 shows the damping range where each attractor is reached.Table 1
Range of existence of fixed points as a function of the damping coefficient withγ
^
0
=
0.02.
Attractor
Damping range
n
1
n
0
n
2
A
8
η
≥
0.04
1
9
17
A
7
η
≥
0.08
2
9
16
A
6
η
≥
0.05
3
9
15
A
5
η
≥
0.04
4
9
14
A
4
0.05
≤
η
≤
0.15
5
9
13
A
3
0.03
≤
η
≤
0.13
6
9
12To obtain the attractors, The attractors are obtained considering the convergence of the energyE
m and the amplitudes of the resonant modes, a
0, a
1, and a
2. For coefficients η less than the indicated values, in some cases the solutions diverge (E
m
→
∞) and the attractor disappears. In other cases, the lower limit is a supercritical Hopf bifurcation that gives rise to periodic solutions which lead to complex attractor, as explained in the next section.The attractors of Table1 are classified as stable fixed points; however, this feature depends on the integration time that is considered. For a maximum time less than 5000 units, the attractors can be still assumed as stable fixed points, but for longer times it is found that, independent of the daughter wavenumbers n
1 and n
2, all configurations converge to the attractor A
8 via a jump. This is, the energy of the daughter waves is transferred to the modes given by n
1
=
1 and n
2
=
17, and the system evolves towards the mentioned attractor. Figure 3 shows the energy evolution for different initial configurations. Note that the convergence to the attractor A
8 occurs from a pseudo-stable position. Similar behaviors can be observed for the evolution the of mother wave amplitude a
0, while for the daughter waves the amplitudes a
1 and a
2 rapidly converge to zero as the energy jumps to the waves of A
8. The convergence of the solutions to this attractor is due to the strong dissipation, for which the system is obligated to put the energy in the unstable mode and in the mode n
1
=
1 that has the lower diffusion, according to the resistive damping model. Consequently, the solution can be affected by the boundary conditions [27].Energy evolutionE
m for different configurations of resonant initial waves with n
0
=
9, η
=
0.12, and γ
^
0
=
0.02 (in configurations with n
1
>
4 and n
2
<
14, the convergence to the attractor A
8 occurs with small transitory stability intervals and thus they are not considered in the analysis).
(a)
(b)
(c)
(d)In Figures4 and 5, the results for attractors A
8 and A
5 are shown. The results of the others attractors are very similar; therefore, the next comments are valid for all attractors presented in Table 1. In the figures, the energy E
m and the amplitudes a
0, a
1, and a
2 are shown as a function of the damping coefficient η, where circles indicate the numerical simulations and the full lines the results of the truncation model. In addition, it also shows the Fourier energy spectrum of the solution, that is, the amplitude of each mode as a function of the wavenumber k, highlighting with dashed lines the resonant waves. In this way, it can be known how the energy is distributed among the different modes when the attractor is formed. It should be noted that, for attractors other than A
8, the comparison between the numerical and the truncation results is only valid at times before the jump to A
8, since the truncation model only considers the three initial waves.AttractorA
8 (n
0
=
9, n
1
=
1, and n
2
=
17 with γ
^
0
=
0.02).
(a)
(b)AttractorA
5 (n
0
=
9, n
1
=
4, and n
2
=
14 with γ
^
0
=
0.02).
(a)
(b)Observing Figures4 and 5, it can be seen that both methods are in good accordance with most of the damping domain, excepting the region η
≲
0.15, where the results are different for each method. The energy spectrum shows that, in addition to the discrepancy in the evaluation of the amplitude of the daughter waves, the differences between the results of E
m are also due to the emergence of two new modes with nonnegligible energy values. These new modes identified as k
3 and k
4, being k
3
<
k
1
<
k
0
<
k
2
<
k
4, satisfy the resonance relation 2
k
0
=
k
1
+
k
2
=
k
3
+
k
4 and an additional relation k
4
-
k
3
=
2
(
k
2
-
k
1
). This behavior is registered in all the attractors of Table 1 when intermediate values of the damping coefficient are used. On the other hand, when larger values of η are considered, there are no new modes with relevant energy levels and the results are similar for both methods.


### 4.2. Dynamic Solutions

In this section, attractors corresponding to small damping coefficients are presented. In these cases a series of complex solutions is obtained, which are formed by limit cycles, chaotic attractors, and processes of intermittency and crises.Table2 shows the damping range for existence of dynamic solutions for the indicated initial wave configurations. These solutions are only captured by the numerical simulations, while the truncation model mainly converges to spurious fixed points. Wave configurations of attractors A
8, A
7, A
5, and A
3 not included in Table 2 have no dynamic solutions, but they diverge for small damping values (E
m
→
∞).Table 2
Range of existence of dynamic solutions as a function of the damping coefficient withγ
^
0
=
0.02.
Attractor
Damping range
n
1
n
0
n
2
A
6
0.0256
≤
η
≤
0.05
3
9
15
A
4
0.0162
≤
η
≤
0.05
5
9
13
A
2
0.0163
≤
η
≤
0.05
7
9
11In Figure6, bifurcation diagrams of the maximums of the energy and the amplitudes are plotted versus the damping coefficient η for attractor A
2. In this way, a one-period limit cycle for a given damping coefficient is identified as a single point in the bifurcation diagram, a two-period limit cycle as two points, and so on, while for a chaotic attractor the solution is broad.Bifurcation diagram of attractorA
2 (n
0
=
9, n
1
=
7, and n
2
=
11 with γ
^
0
=
0.02). Points indicate the numerical simulations. The truncation solutions are indicated by the full line (branch of fixed points) and the dashed line (maximum of the periodic solution).
(a)
(b)
(c)
(d)According to the numerical results, the attractorA
2 starts at η
c
≈
0.05 with a supercritical Hopf bifurcation where the fixed points solutions that rapidly converge to attractor A
8 become stable periodic solutions modulated in time for η
≲
η
c. In η
≈
0.0467, there is a period doubling bifurcation followed by a Feigenbaum sequence ending in a chaotic attractor. The periodic solution emerges again at η
≈
0.0325 through what appears to be an intermittent process. Then a new Feigenbaum sequence continues, which starts at η
≈
0.029 ending in a new chaotic attractor. By means of another similar sequence, the periodic solution returns to finally culminate in a chaotic attractor that ends the periodic solutions at η
≈
0.016. Finally, the chaotic attractor is suddenly destroyed at η
≈
0.0163 where a boundary crisis appears to take place.With respect to the solutions of the truncation model, it can be seen that the method does not capture the chaotic attractors, but it predicts the existence of a branch of stable fixed points that ends in a Hopf bifurcation atη
≈
0.023, where periodic solutions take place.The discrepancies between both methods are due to the energy transfer, since the energy initially contained in the three resonant waves is transferred to other modes which are not considered in the three-wave truncation model. This feature confirms the results obtained in the analysis of stationary solutions, where it was shown that for small damping new additional modes emerge; although in this case the energy spectrum is broader than only five modes.Figures7 and 8 show the bifurcation diagrams for attractors A
4 and A
6, respectively. For A
4 similar results to attractor A
2 can be observed, where the same Feigenbaum sequences and chaotic attractors at similar damping ranges are registered; although in this case the results of the truncation model are rather worse than those obtained for A
2. On the other hand, in the bifurcation diagram of the daughter waves can be observed what appears to be an attractor merging crisis. With respect to attractor A
6, the above comments are valid for damping coefficients η
≳
0.0256, for which the attractor behaves similarly to attractor A
4. But in this case the boundary crisis take place for η
≈
0.0256 when the attractor is suddenly destroyed.Bifurcation diagram of attractorA
4 (n
0
=
9, n
1
=
5, and n
2
=
13 with γ
^
0
=
0.02).
(a)
(b)
(c)
(d)Bifurcation diagram of attractorA
6 (n
0
=
9, n
1
=
3, and n
2
=
15 with γ
^
0
=
0.02).
(a)
(b)
(c)
(d)To complete the analysis, Figure9 is presented in order to visualize the intermittent process by which the periodic solution is restored at η
c
≈
0.0325. In the figure the time evolution of the mother wave is shown, where it can be seen that for a minimum increment of the damping coefficient the periodic solution is replaced by pseudoperiodic evolutions that are interrupted by chaotic explosions, responding to the characteristic shape of the intermittent phenomenon. According to the bifurcation diagrams, a tangent bifurcation appears to take place at η
c, which indicates the presence of type-I Pomeau-Manneville intermittency [37]. In such a phenomenon, the emergence of chaotic solutions is related to the coalescence of two fixed points in a tangent bifurcation that disappear when the control parameter is changed, generating a narrow channel through which the orbits slowly evolve forming the characteristic laminar phases of the intermittency phenomenon. To visualize this situation, Figure 10 shows a one-dimensional Poincaré map built with the values of the peaks of a
0 for η
=
0.032507. This figure shows the presence of a narrow channel between the map and the bisectrix, evidencing the process of type-I intermittency. A detailed analysis of this phenomenon was recently carried out, where a new methodology to obtain the statistical properties of type-I intermittency with discontinuous reinjection probability density is presented for the solutions of the DNLS equation [38].Time evolution of the mother wave for attractorA
2. (a) Periodic solution. (b) Chaos.
(a)
n
=
0.032506
(b)
n
=
0.032507Figure 10
Poincaré map of the peaks ofa
0 for η
=
0.032507 in the region of the tangent point a
0
≈
0.242 where the phenomenon of type-I intermittency is observed.


## 4.1. Stationary Solutions

When intermediate and strong diffusion levels (η
≳
0.05) are considered, a set of attractors is obtained which consist of stationary traveling waves with most of the energy concentrated in the three initial modes, such that its amplitudes and the energy can be considered as stable fixed points. Table 1 shows the damping range where each attractor is reached.Table 1
Range of existence of fixed points as a function of the damping coefficient withγ
^
0
=
0.02.
Attractor
Damping range
n
1
n
0
n
2
A
8
η
≥
0.04
1
9
17
A
7
η
≥
0.08
2
9
16
A
6
η
≥
0.05
3
9
15
A
5
η
≥
0.04
4
9
14
A
4
0.05
≤
η
≤
0.15
5
9
13
A
3
0.03
≤
η
≤
0.13
6
9
12To obtain the attractors, The attractors are obtained considering the convergence of the energyE
m and the amplitudes of the resonant modes, a
0, a
1, and a
2. For coefficients η less than the indicated values, in some cases the solutions diverge (E
m
→
∞) and the attractor disappears. In other cases, the lower limit is a supercritical Hopf bifurcation that gives rise to periodic solutions which lead to complex attractor, as explained in the next section.The attractors of Table1 are classified as stable fixed points; however, this feature depends on the integration time that is considered. For a maximum time less than 5000 units, the attractors can be still assumed as stable fixed points, but for longer times it is found that, independent of the daughter wavenumbers n
1 and n
2, all configurations converge to the attractor A
8 via a jump. This is, the energy of the daughter waves is transferred to the modes given by n
1
=
1 and n
2
=
17, and the system evolves towards the mentioned attractor. Figure 3 shows the energy evolution for different initial configurations. Note that the convergence to the attractor A
8 occurs from a pseudo-stable position. Similar behaviors can be observed for the evolution the of mother wave amplitude a
0, while for the daughter waves the amplitudes a
1 and a
2 rapidly converge to zero as the energy jumps to the waves of A
8. The convergence of the solutions to this attractor is due to the strong dissipation, for which the system is obligated to put the energy in the unstable mode and in the mode n
1
=
1 that has the lower diffusion, according to the resistive damping model. Consequently, the solution can be affected by the boundary conditions [27].Energy evolutionE
m for different configurations of resonant initial waves with n
0
=
9, η
=
0.12, and γ
^
0
=
0.02 (in configurations with n
1
>
4 and n
2
<
14, the convergence to the attractor A
8 occurs with small transitory stability intervals and thus they are not considered in the analysis).
(a)
(b)
(c)
(d)In Figures4 and 5, the results for attractors A
8 and A
5 are shown. The results of the others attractors are very similar; therefore, the next comments are valid for all attractors presented in Table 1. In the figures, the energy E
m and the amplitudes a
0, a
1, and a
2 are shown as a function of the damping coefficient η, where circles indicate the numerical simulations and the full lines the results of the truncation model. In addition, it also shows the Fourier energy spectrum of the solution, that is, the amplitude of each mode as a function of the wavenumber k, highlighting with dashed lines the resonant waves. In this way, it can be known how the energy is distributed among the different modes when the attractor is formed. It should be noted that, for attractors other than A
8, the comparison between the numerical and the truncation results is only valid at times before the jump to A
8, since the truncation model only considers the three initial waves.AttractorA
8 (n
0
=
9, n
1
=
1, and n
2
=
17 with γ
^
0
=
0.02).
(a)
(b)AttractorA
5 (n
0
=
9, n
1
=
4, and n
2
=
14 with γ
^
0
=
0.02).
(a)
(b)Observing Figures4 and 5, it can be seen that both methods are in good accordance with most of the damping domain, excepting the region η
≲
0.15, where the results are different for each method. The energy spectrum shows that, in addition to the discrepancy in the evaluation of the amplitude of the daughter waves, the differences between the results of E
m are also due to the emergence of two new modes with nonnegligible energy values. These new modes identified as k
3 and k
4, being k
3
<
k
1
<
k
0
<
k
2
<
k
4, satisfy the resonance relation 2
k
0
=
k
1
+
k
2
=
k
3
+
k
4 and an additional relation k
4
-
k
3
=
2
(
k
2
-
k
1
). This behavior is registered in all the attractors of Table 1 when intermediate values of the damping coefficient are used. On the other hand, when larger values of η are considered, there are no new modes with relevant energy levels and the results are similar for both methods.


## 4.2. Dynamic Solutions

In this section, attractors corresponding to small damping coefficients are presented. In these cases a series of complex solutions is obtained, which are formed by limit cycles, chaotic attractors, and processes of intermittency and crises.Table2 shows the damping range for existence of dynamic solutions for the indicated initial wave configurations. These solutions are only captured by the numerical simulations, while the truncation model mainly converges to spurious fixed points. Wave configurations of attractors A
8, A
7, A
5, and A
3 not included in Table 2 have no dynamic solutions, but they diverge for small damping values (E
m
→
∞).Table 2
Range of existence of dynamic solutions as a function of the damping coefficient withγ
^
0
=
0.02.
Attractor
Damping range
n
1
n
0
n
2
A
6
0.0256
≤
η
≤
0.05
3
9
15
A
4
0.0162
≤
η
≤
0.05
5
9
13
A
2
0.0163
≤
η
≤
0.05
7
9
11In Figure6, bifurcation diagrams of the maximums of the energy and the amplitudes are plotted versus the damping coefficient η for attractor A
2. In this way, a one-period limit cycle for a given damping coefficient is identified as a single point in the bifurcation diagram, a two-period limit cycle as two points, and so on, while for a chaotic attractor the solution is broad.Bifurcation diagram of attractorA
2 (n
0
=
9, n
1
=
7, and n
2
=
11 with γ
^
0
=
0.02). Points indicate the numerical simulations. The truncation solutions are indicated by the full line (branch of fixed points) and the dashed line (maximum of the periodic solution).
(a)
(b)
(c)
(d)According to the numerical results, the attractorA
2 starts at η
c
≈
0.05 with a supercritical Hopf bifurcation where the fixed points solutions that rapidly converge to attractor A
8 become stable periodic solutions modulated in time for η
≲
η
c. In η
≈
0.0467, there is a period doubling bifurcation followed by a Feigenbaum sequence ending in a chaotic attractor. The periodic solution emerges again at η
≈
0.0325 through what appears to be an intermittent process. Then a new Feigenbaum sequence continues, which starts at η
≈
0.029 ending in a new chaotic attractor. By means of another similar sequence, the periodic solution returns to finally culminate in a chaotic attractor that ends the periodic solutions at η
≈
0.016. Finally, the chaotic attractor is suddenly destroyed at η
≈
0.0163 where a boundary crisis appears to take place.With respect to the solutions of the truncation model, it can be seen that the method does not capture the chaotic attractors, but it predicts the existence of a branch of stable fixed points that ends in a Hopf bifurcation atη
≈
0.023, where periodic solutions take place.The discrepancies between both methods are due to the energy transfer, since the energy initially contained in the three resonant waves is transferred to other modes which are not considered in the three-wave truncation model. This feature confirms the results obtained in the analysis of stationary solutions, where it was shown that for small damping new additional modes emerge; although in this case the energy spectrum is broader than only five modes.Figures7 and 8 show the bifurcation diagrams for attractors A
4 and A
6, respectively. For A
4 similar results to attractor A
2 can be observed, where the same Feigenbaum sequences and chaotic attractors at similar damping ranges are registered; although in this case the results of the truncation model are rather worse than those obtained for A
2. On the other hand, in the bifurcation diagram of the daughter waves can be observed what appears to be an attractor merging crisis. With respect to attractor A
6, the above comments are valid for damping coefficients η
≳
0.0256, for which the attractor behaves similarly to attractor A
4. But in this case the boundary crisis take place for η
≈
0.0256 when the attractor is suddenly destroyed.Bifurcation diagram of attractorA
4 (n
0
=
9, n
1
=
5, and n
2
=
13 with γ
^
0
=
0.02).
(a)
(b)
(c)
(d)Bifurcation diagram of attractorA
6 (n
0
=
9, n
1
=
3, and n
2
=
15 with γ
^
0
=
0.02).
(a)
(b)
(c)
(d)To complete the analysis, Figure9 is presented in order to visualize the intermittent process by which the periodic solution is restored at η
c
≈
0.0325. In the figure the time evolution of the mother wave is shown, where it can be seen that for a minimum increment of the damping coefficient the periodic solution is replaced by pseudoperiodic evolutions that are interrupted by chaotic explosions, responding to the characteristic shape of the intermittent phenomenon. According to the bifurcation diagrams, a tangent bifurcation appears to take place at η
c, which indicates the presence of type-I Pomeau-Manneville intermittency [37]. In such a phenomenon, the emergence of chaotic solutions is related to the coalescence of two fixed points in a tangent bifurcation that disappear when the control parameter is changed, generating a narrow channel through which the orbits slowly evolve forming the characteristic laminar phases of the intermittency phenomenon. To visualize this situation, Figure 10 shows a one-dimensional Poincaré map built with the values of the peaks of a
0 for η
=
0.032507. This figure shows the presence of a narrow channel between the map and the bisectrix, evidencing the process of type-I intermittency. A detailed analysis of this phenomenon was recently carried out, where a new methodology to obtain the statistical properties of type-I intermittency with discontinuous reinjection probability density is presented for the solutions of the DNLS equation [38].Time evolution of the mother wave for attractorA
2. (a) Periodic solution. (b) Chaos.
(a)
n
=
0.032506
(b)
n
=
0.032507Figure 10
Poincaré map of the peaks ofa
0 for η
=
0.032507 in the region of the tangent point a
0
≈
0.242 where the phenomenon of type-I intermittency is observed.


## 5. Conclusions

In this work, the derivative nonlinear Schrödinger (DNLS) equation for periodic boundary conditions has been numerically solved using spectral methods based on Fourier expansions for the spatial derivatives and a fourth order Runge-Kutta scheme for time integration. The left-hand polarized Alfvén waves in a homogeneous uniform plasma withβ
≈
0 were considered in the analysis.The DNLS equation is extensively used to describe Alfvén waves propagating parallel (or quasiparallel) to the magnetic field direction. Alfvén waves are ubiquitous in space and play an important role in the solar wind, interstellar medium, turbulent heating in magnetosphere, and among others [39, 40]. In these astrophysical processes the Hall effect becomes relevant; thus the magnetohydrodynamics (MHD) equations have to be modified, resulting in a computationally demanding problem. However, for the description of the nonlinear dynamics of large amplitude Alfvén waves traveling along to the background magnetic field, the solution of the full Hall-MHD equations can be replaced by the analysis of the DNLS equation, which represents a large reduction in the computational cost when numerical methods are used. The validity of this idealized model will always depend on the validity of the assumptions carried out during its derivation, in addition to the appropriateness of the use of a fluid-like damping model and a one-wave driver term [8, 18, 20]. Having into account these circumstances, the DNLS equation is widely used to explain many phenomena of astrophysics, such as the Alfvén turbulence [22, 23], the solar wind [41], and the wave-wave interactions [25].The first analysis of this work was carried out to validate the numerical scheme by testing the analytical conditions of modulational stability in the nondiffusive case with a one-wave initial condition. The stability conditions were adequately reproduced by the code. In addition, the instability time could be calculated and a study of the evolution of the unstable configurations was carried out, features not predicted by previous theoretical analysis. In this way the numerical results showed that the wavenumber of the initial wave defines the instability time of the system (higher modes destabilize faster), while the amplitude governs the shape of the evolution: small amplitudes evolve as periodic solutions where the initial wave recovers its energy in each cycle, while greater amplitudes produce spatio-temporal chaos where the energy is randomly distributed among all the modes with an irregular evolution.In the study of the driven dissipative DNLS equation a three-wave initial condition was considered, with the mother wave (k
0) linearly unstable and the remaining waves damped, with the daughter waves (k
1 and k
2) satisfying the resonant condition (2
k
0
=
k
1
+
k
2). This analysis showed that the diffusion level of the system defines the kind of solutions. For relatively large damping a set of attractors consisting in pseudo-stable fixed points is found, which converge to the less diffusive wave configuration via a jump for sufficiently long times. Before the jump, each attractor behaves as a conventional fixed point. In these cases, for larger damping levels, the energy remains contained in the three initial modes, while for smaller diffusion it is distributed to two new modes k
3 and k
4 (being k
3
<
k
1
<
k
0
<
k
2
<
k
4) that always satisfy the resonant condition 2
k
0
=
k
1
+
k
2
=
k
3
+
k
4 and are related to the daughter waves by k
4
-
k
3
=
2
(
k
2
-
k
1
). The emergence of these new modes produces the errors in the results of the three-wave truncation model, which can be considered a good approximation for large damping but fails when the new modes become relevant.When very small damping coefficients are used, the numerical solutions drastically change resulting in a complex dynamics. The bifurcation diagrams of these solutions exhibit limit cycles, chaotic attractors, and processes of intermittency and crisis with the change of the damping coefficient. On the other hand, the three-wave truncation model cannot capture the dynamic solution and predicts a spurious branch of fixed points instead. This is because the power spectrum forη
→
0 is broad and consequently the truncation model is unable to represent the dynamics of the system. As a result, this simplified model cannot be used for small diffusion levels.These results allow concluding that the numerical scheme proposed in this paper is a good option to solve the driven dissipative DNLS equation, while the three-wave truncation model is a good approach when relatively large damping levels are considered, but for small diffusion the energy transfer does not allow the use of the simplified model.Finally, according to the results, it is concluded that the numerical method presented in this work is a good option to solve the driven dissipative DNLS equation, which has been properly validated and by which new results were obtained.



---

*Source: 812052-2014-04-14.xml*