PHYSICAL REVIEW E 89, 043103 (2014)

Dust-acoustic shocks in strongly coupled dusty plasmas

S. E. Cousens,
1,*
V. V. Yaroshenko,
2
S. Sultana,
3
M. A. Hellberg,
4
F. Verheest,
5,4
and I. Kourakis
1,
1
Centre for Plasma Physics, Department of Physics and Astronomy, Queens University Belfast, BT7 1NN Northern Ireland, United Kingdom
2
Deutsches GeoForschungsZentrum, GFZ, 14473 Potsdam, Germany
3
Department of Physics, Jahangirnagar University, Dhaka-1342, Bangladesh
4
School of Chemistry and Physics, University of KwaZulu-Natal, Durban 4000, South Africa
5
Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281, B-9000 Gent, Belgium
(Received 18 July 2013; published 4 April 2014)
Electrostatic dust-acoustic shock waves are investigated in a viscous, complex plasma consisting of dust
particles, electrons, and ions. The system is modelled using the generalized hydrodynamic equations, with strong
coupling between the dust particles being accounted for by employing the effective electrostatic temperature
approach. Using a reductive perturbation method, it is demonstrated that this model predicts the existence of
weakly nonlinear dust-acoustic shock waves, arising as solutions to Burgerss equation, in which the nonlinear
forces are balanced by dissipative forces, in this case, associated with viscosity. The evolution and stability
of dust-acoustic shocks is investigated via a series of numerical simulations, which conrms our analytical
predictions on the shock characteristics.
DOI: 10.1103/PhysRevE.89.043103 PACS number(s): 52.27.Lw, 52.35.Fp, 52.35.Tc, 52.27.Gr
I. INTRODUCTION
At small amplitude, waves in plasmas may often be
treated as linear phenomena, but as the amplitude of these
waves increases, various processes occur which may not be
described by linear theory. An important nonlinear effect is
that the higher-amplitude parts of a wave travel faster than
the lower-amplitude parts, leading to a steepening of the
wave front. Wave steepening is not the only mechanism to
affect waves in a plasma, and its effects may be mitigated, or
even balanced, by other processes. An example of one such
effect is dispersion, by which various frequency components
in the wave travel at different speeds, causing the wave to
be stretched over time. When the effect of wave steepening
is balanced by dispersion, wave structures known as solitary
waves, or solitons, may be produced. Another process is
energy dissipation or damping of the waves in the medium.
This involves the transfer of energy from the wave, and a
balance between wave steepening and dissipation results in
the formation of shock waves. Various techniques have been
developed to analytically model nonlinear waves, such as the
Sagdeev pseudopotential approach [1] for solitons, and the
reductive perturbation method [2], which may be employed
whenever the amplitude of the nonlinear wave is small.
A very common type of plasma, and one which is of
particular interest to this paper, is a complex or dusty plasma,
which consists of massive dust particles in addition to the
usual electron and ion components [3]. The dust particles
may attain a large net negative charge, since the electrons
typically have a larger thermal speed than the ions. Initially,
the electrons collide with the dust particles more frequently
than the ions. As the dust grains become more negatively
charged, electrons are increasingly repelled, reducing the
electron ux, while at the same time the ion ux increases,
until the equilibrium dust charge is reached when the electron
*
[email protected]

[email protected]; www.kourakis.eu
and ion uxes are equal and the charging process is complete.
An understanding of dusty plasma physics is necessary to
explain various astrophysical phenomena [4] and is important
for many industrial and manufacturing applications [5]. The
introduction of a massive, charged dust component is found to
greatly affect the dynamics of a usual electron-ion plasma. This
includes both a modication to existing linear wave modes, and
the introduction of new ones, such as the dust-acoustic wave
(DAW), in which the pressures of the relatively inertialess
electrons and ions provide the restoring force for the inertial
dust [6,7]. It is interesting that the large size of the dust
particles, as well as their low velocities, allow scientists to
view the propagation of dust-acoustic waves in the laboratory
in real time [8].
Another interesting aspect of dusty plasma physics is
that the dust particles may experience signicant electro-
static forces from their neighbors. The coupling parameter
= Z
2
d
e
2
/4
0
k
B
T
d
n
1/3
d
is the ratio between the Coulomb
potential energy and the thermal energy of the dust and is used
as measure of the relative importance of these electrostatic
forces on the motion of the dust particles. Here Z
d
e, k
B
T
d
, and
n
d
are the charge, kinetic temperature, and number density
of the dust particles, respectively. The high charge number
and low kinetic temperature of dust particles often leads to
coupling parameters that are greater than unity. In such plasmas
the electrostatic forces are of great importance and we class
the dust particles as being strongly coupled [9]. When the
coupling parameter gets sufciently large, crystallization can
occur, and so dusty plasmas provide an excellent opportunity
to study phase transitions, which is another very active topic of
research. In a strongly coupled dusty plasma, viscosity often
plays an important role in the dynamics of the dust particles and
of the properties of the wave modes. Of particular relevance
to this study is that the dissipation of energy through viscous
forces may provide a balance for nonlinear forces that drive
the steepening of the wave front and allowstable dust-acoustic
shocks to form.
Gozadinos et al. [10] presented a model for strongly cou-
pled crystalline dusty plasmas, whereby the effects of strong
1539-3755/2014/89(4)/043103(9) 043103-1 2014 American Physical Society
S. E. COUSENS et al. PHYSICAL REVIEW E 89, 043103 (2014)
coupling are accounted for via an electrostatic pressure. This
model, although originally developed for crystalline plasma
structures, has recently been applied as an approximation
to the equation of state for strongly coupled plasmas near
to the liquid-crystal phase transition. This has included the
study of Bohm sheaths [11], double-layer formation [12,13],
the linear DAW mode [14], and nonlinear solitary wave
structures [15,16]. Yaroshenko et al. [14] found that, by
interpreting the strongly coupled interactions as a consequence
of an electrostatic temperature, the transition to a thermal
mode at high wave numbers for dispersion curves obtained
in previous experimental studies [17] may be accurately
predicted.
In Ref. [18] we theoretically investigated weakly nonlinear
dust-acoustic solitary waves in a complex plasma in which the
dust particles are strongly coupled, with the effects of strong
coupling being accounted for by the electrostatic temperature
approach. To account for local variations in electrostatic
temperature in the vicinity of a wave, we introduced the
concept of electrostatic temperature perturbations. In doing
so, we noted that including the electrostatic temperature
perturbations predicts modications to the sound speed in the
plasma, as well as the amplitude and width of solitary waves.
In this paper, we extend this model to study dust-acoustic
shock waves in a dusty plasma in which the dust particles
are in a strongly coupled, but still uid, state. In particular,
since we use the electrostatic pressure approach of Gozadinos
et al. [10], our results are applicable in the vicinity of the
liquid-crystal phase transition. This state is often observed
in the laboratory, signied by a coexistence of crystalline
structures alongside a viscous uidlike plasma [19]. Adetailed
description of the model, along with the key assumptions
made, is presented in Sec. II. We use the reductive perturbation
method [2], described in Sec. III, to derive both the linear
dispersion relation and Burgerss equation for this system in
Sec. IV. In addition to this, we then present a discussion
on dust-acoustic shocks, as they arise as solutions to the
derived Burgerss equation, for the plasma parameters given
by Rosenberg et al. [20]. In Sec. V, we present a numerical
simulation for dust-acoustic shocks propagating in a plasma
system for which they are not a solution to Burgerss equation.
This numerical study is based on the parameters presented
in Ref. [20] and employs a Runge-Kutta iteration method to
integrate Burgerss equation using a specied shock potential
prole as an initial condition. Finally, in Sec. VI we discuss
the key ndings of this approach.
II. MODEL SYSTEM
In this section, we present a uid model for dust-acoustic
waves in a one-dimensional, three-component plasma consist-
ing of dust particles, electrons, and ions. This will include a
brief discussion on the key assumptions we make in the model
and the equations that result from them.
In any particular dusty plasma system the dust particles
may have a range of sizes, with charge numbers varying both
spatially and temporally. In this paper, we take a simplied
case such that the dust particles each have a constant mass m
d
and a xed charge number Z
d
. The latter means that we do
not include charge exchange processes, which might result in
the dust particles gaining or losing charges, in our analysis.
This commonly used assumption implies that the study is
restricted to plasmas in which the typical (dust-acoustic) time
scale differs signicantly from (and is, in fact, much larger
than) the dust grain charging time. We note, for completeness,
that in some situations (where the grain charging time scale is
comparable with the excited wave time scale) grain charging
processes do have effects on the dynamics of the system [21].
The dust number density, n
d
, is thus conserved, and this is
reected in the continuity equation,
n
d
t
+

x
(n
d
u
d
) = 0, (1)
where u
d
is the dust uid velocity.
We consider three forces on the dust particles. The rst
is the electrostatic force caused by separation of charges in
the plasma. We assume that there are no external electric
elds applied to the system, and so this charge separation
is caused solely by the plasma being perturbed by a passing
wave. We express this force as the gradient of the electrostatic
potential .
The second force we consider is due to the mutual
electrostatic repulsion of similarly charged dust particles. We
model a uid plasma that is close to the liquid-crystal phase
transition, and so to obtain a rst approximation for this force,
we use the equation of state derived by Gozadinos et al. [10]
for dust particles in a crystalline structure. In this model,
the coupling energy per particle, E
0
, is assumed to be of the
isotropic screened Coulomb (Debye-H uckel or Yukawa) form
between the particle and its nearest neighbors, such that
E
0
N
nn
Z
2
d
e
2
4
0
a
exp () , (2)
where N
nn
is the number of nearest neighbors, determined by
the crystalline plasma structure, and = a/
D
is the normal-
ized interparticle distance, dened as the ratio between the
dust interparticle distance a = n
1/3
d
and the Debye screening
length in the plasma,
D
=
_

0
k
B
T
i
T
e
/e
2
(n
i
T
e
+n
e
T
i
). By
assuming that the coupling energy is the dominant source of
free energy, F, in the system such that F N
d
E
0
, where
N
d
is the number of dust particles, Gozadinos et al. [10]
determined the pressure via the thermodynamic relation P

=
(F/V)
T
d
n
2
d
(E
0
/n
d
)
T
d
. Here V denotes the volume
of the plasma, such that V = N
d
/n
d
. The force exerted on the
dust particles is then the gradient of this electrostatic pressure,
P

/x. Following the work of Yaroshenko et al. [14], we

dene the electrostatic temperature via the equation P

=
n
d
k
B
T

, which yields
k
B
T

=
N
nn
Z
2
d
e
2
12
0
3

n
d
(1 +)exp(). (3)
For the purposes of this paper, we take N
nn
= 12, which is the
value for a face-centered cubic lattice. Since we model a uid
plasma, this value of N
nn
is used only for the sake of analytical
tractability, to proceed with numerical computations, and
should be taken as an approximation for the system.
We note that we have employed the Debye-H uckel (or
Yukawa) model in the derivation of the electrostatic pressure.
This is commonly used in the literature [19,22] and has, for
043103-2
DUST-ACOUSTIC SHOCKS IN STRONGLY COUPLED . . . PHYSICAL REVIEW E 89, 043103 (2014)
instance, proven to be useful in providing results that are
supported by experiments [23]. For the idealized dusty plasma
we describe in this paper, the Debye-H uckel model is therefore
used as an approximation to the interaction potential of the dust
particles.
Third, we introduce a dissipative force, representing mo-
mentumtransfer in the system. In principle, energy dissipation
in our model is not due to collisions with neutrals but rather
is associated with the intrinsic longitudinal viscosity,
l
, of
the dust uid. This transport coefcient (
l
) is connected
to the shear
s
and bulk,
b
viscosities via the relation

l
= (4
s
/3) +
b
. Obtaining these quantities is not a straight-
forward task and, in general, they must be found empirically
by either experimental measurement or molecular dynamic
(MD) simulations [24], in conjunction with, for example, the
Green-Kubo relations [25]. Such simulations are generally
based on the Debye-H uckel model and are applicable in a broad
range of physical systems. Anumber of papers provide tabular
and graphical results obtained via these methods [2630]. For
this study, we note that the bulk viscosity is typically three
orders of magnitude smaller than the shear viscosity [29], and
so we will only consider the shear component in this paper,
such that
l

s
. We obtain this value of
s
from the work of
Sanbonmatsu and Murillo [28], who spatially modulated the
velocity of the particles in a MD simulation and calculated
the shear viscosity from the relaxation time of the modulation
prole for a range of values of and .
The net force on the dust particles is thus a combination
of these three forces, and this is reected in the momentum
equation,
m
d
n
d
_
u
d
t
+u
d
u
d
x
_
= Z
d
en
d

x
+
l

2
u
d
x
2
.
(4)
Finally, Poissons equation gives the electrostatic potential
generated from a given distribution of charges,

x
2
=
e

0
(n
i
n
e
Z
d
n
d
) , (5)
where we have assumed that the ions have a charge number
Z
i
= 1 for simplicity. Since the mass of the electrons and ions
is much less than that of the dust particles, we assume that,
when the dust particles are perturbed due to the passing of the
wave, the electrons and ions instantaneously redistribute them-
selves according to the Maxwell-Boltzmann distribution, such
that n
e
= n
e0
exp (e/k
B
T
e
) and n
i
= n
i0
exp (e/k
B
T
i
).
We note that though we do not include ionization of
neutral particles in our model, the physical reality is that
electrons and ions continuously arrive at the dust particles
and recombine on their surfaces, and so ionization must be
taking place somewhere in the systemto maintain this dynamic
equilibrium. Furthermore, we note that at low gas pressures,
plasma ionization becomes important for ion-acoustic and
dust-acoustic waves as demonstrated by DAngelo [31,32].
This effect is omitted here for simplicity and analytical
tractability.
For notational clarity, we normalize the dynamic variables
appearing in Eqs. (1), (4), and (5). The effective electrostatic
temperature is scaled such that d = T

/T
0
, where T
0
=
Z
2
d
n
d0
T
i
T
e
/(n
i0
T
e
+n
e0
T
i
). The dust density, dust velocity,
and electrostatic potential are then scaled as n = n
d
/n
d0
, u =
u
d
/

k
B
T
0
/m
d
and = /(k
B
T
0
/Z
d
e), respectively. The
spacelike and timelike variables are scaled by the Debye
length,
D0
=
_

0
k
B
T
0
/n
d0
Z
2
d
e
2
, and plasma period
1
pd
=
_

0
m
d
/n
d0
Z
2
d
e
2
, respectively, and the normalized version of
these quantities is denoted with a tilde. The longitudinal
viscosity is normalized such that =
l
/m
d
n
d0

pd

2
D0
. Equa-
tions (1), (4), and (5) thus become
n

t
+

x
(nu) = 0, (6)
n
_
u

t
+u
u
x
_
= n

(nd)
x
+

2
u
x
2
, (7)
and

x
2
= n +

1
exp (
e
)
1
1
exp (
i
) , (8)
with = n
e0
/n
i0
, = T
i
/T
e
,
i
= (1 )/(1 +), and

e
=
i
. This is a modied version of the well-known
generalized hydrodynamic model [3335], with the inclusion
of the dust electrostatic pressure term accounting for strong-
coupling effects.
Since we are interested in weakly nonlinear waves and there
are no external electric elds applied to the system, we assume
that the normalized potential is small, such that 1. As a
consequence of this, we may expand the exponential functions
appearing in Eq. (8) such that exp(j) = 1 +j +(j)
2
/2 +
O(
3
). For the purpose of deriving Burgerss equation for
this system, we only require up to the
2
term. With that
approximation, Poissons equation simplies to

x
2
= (n 1) +c
1
+c
2

2
+O(
3
), (9)
where
c
1
= 1, c
2
=
(1 )(1
2
)
2(1 +)
2
. (10)
III. REDUCTIVE PERTURBATION METHOD
Since we are dealing with weak nonlinearities, we may lin-
earize the equations using the reductive perturbation method,
whereby we expand the dynamical quantities in terms of a
small parameter . We then expand the dust number density n,
the dust uid velocity u, and the electrostatic potential such
that
n = 1 +n
1
+
2
n
2
+ , (11)
u = u
1
+
2
u
2
+ , (12)
and
=
1
+
2

2
+ . (13)
As mentioned in Sec. I, T

is also a dynamically varying

quantity, depending on both and n
d
. For notational clarity,
we expand the normalized electrostatic temperature, d =
T

/T
0
separately, such that
d = d
0
+d
1
+
2
d
2
+ , (14)
043103-3
S. E. COUSENS et al. PHYSICAL REVIEW E 89, 043103 (2014)
with the perturbations d
1
and d
2
being
d
1
= d
11
n
1
+d
12

1
(15)
and
d
2
= d
21
n
2
+d
22

2
+d
23
n
2
1
+d
24
n
1

1
+d
25

2
1
. (16)
The expressions for d
ij
are presented in the Appendix, and a
derivation of these quantities is given in Ref. [18].
For the sake of rigor, it may be noted that, since the
viscosity coefcient is a function of both and (see,
e.g., Ref. [24]), its value may in principle be perturbed by
the passing shock. Therefore, one might consider
0
+

1
+
2

2
+ . However, in deriving Burgerss equation
in Sec. IV below, we only require orders up to
3
for the
viscosity term in the momentum equation. Given our choice
of stretched coordinates, only the equilibrium term of the
viscosity contributes to terms of this order, so all higher-order
perturbations of are not required for the present study, hence,
from this point forward we take
0
.
IV. RESULTS
A. Linear
To obtain the linear dispersion relation for this system,
we expand the state variables in Eqs. (6), (7), and (9) in
a power series of , as described above, with terms
2
or
higher neglected. By then assuming oscillatory solutions to the
perturbed quantities, such that /

t i and / x i

k,
we obtain the linear dispersion relation

2
+i

k
2
=

k
2
1 +

k
2
+

k
2
, (17)
in which we have
= 1 d
12
, = d
0
+d
11
. (18)
It may be deduced from their form in the Appendix that the
electrostatic temperature coefcients d
11
and d
12
are positive
quantities for dusty plasmas for which both and are
less than unity. For the parameters considered in this paper,
these two conditions are satised. From this observation, we
see that in this case 1 and 0. For weakly coupled
dusty plasmas (for which
0
1 and 1) the electrostatic
temperature terms tend to zero, in which case = 1 and = 0.
For the purposes of the following discussion, we separate
the dispersion relation into its real and imaginary components
by rst setting =
R
+i
I
. Thus, by Eq. (17), we ob-
tain
2
R

2
I

I

k =

k
2
/(1 +

k
2
) +

k
2
and
R
(2
I
+

k
2
) = 0 as the real and imaginary components, respectively.
For nonzero real frequency, the imaginary component reduces
to
I
=

k
2
/2. This reects the energy dissipation associ-
ated with the viscosity.
We rst consider the nondissipative case ( = 0), which
leads to the dispersion relation
2
=

k
2
/(1 +

k
2
) +

k
2

2
0
. From this, it is found immediately that for both weak
and strong-coupling models an acoustic mode is obtained for
small

k, with the real frequency being proportional to

k, such
that
R
( +)
1/2

k. While in the strongly coupled case

R
remains proportional to

k for large

k, with
R

1/2

k, in the
weakly coupled model the frequency saturates at =
pd
, the
dust plasma frequency for the weakly coupled mode.
The introduction of nite dissipation leads to
2
R
=
2
0
(

k)

k
4
/4, and so for

k 1 this is seen to have only a small effect
on the dispersion, which is still described by an acoustic curve
of slightly lower value. However, an interesting feature of the
damped case ( > 0) is that it follows directly from Eq. (16)
that for

k greater than some critical value

k
cr
, the real frequency

R
= 0. This is given by

k
cr
=
_
(
2
4) +
_
(
2
+4)
2
+16
2
)
2
2
_
1/2
(19)
and applies to both the weakly and strongly coupled cases.
Further, for both cases, the damping (|
I
|) increases as

k
2
(viz.
I
=

k
2
/2) for

k <

k
cr
. Physically, this corresponds
to a nonpropagating damped mode (further algebraic details
omitted here). Note that the limit of the above expression for
vanishing is innity, thus the condition for the nonpropagat-
ing mode to occur (

k >

k
crit
) is never satised in this case, as
expected.
As

k is increased through

k
cr
, two zero-frequency, strongly
damped modes appear. Both the phase and group velocities
vanish for such zero-frequency modes, and no energy is trans-
mitted. They do not appear to have any physical signicance.
In particular, as they occur for short wavelengths, they have no
effect on the Burgerss discussion that follows, as that assumes
small

k. We note that the zero-frequency modes arise from the
imaginary termin the dispersion relation and thus also occur in
the case of other dissipation mechanisms such as dust-neutral
collisions.
This behavior of the wave modes is illustrated in Fig. 1,
in which we present the real and imaginary parts of the
frequency as functions of wave number, both with and without
the dissipative term. Figure 1(a) shows the weakly coupled
case, with = 1 and = 0, while in Fig. 1(b) we have
included the effects of strong coupling. The parameters used
are taken fromRosenberg et al. [20] and yield values = 0.79
and = 1.2. Note that, contrary to Ref. [20], we do not
assume that the dust particles are characterized by high kinetic
temperatures of T
d
1080 eV but use a standard assumption
that the dust is in thermal equilibrium with ions, T
d
T
i
.
The value of the viscosity was obtained by rst calculating
the equilibrium coupling parameter, , and the normalized
interparticle distance,
0
, for this system. These values were
then compared to Figs. 1 and 2 in Ref. [28] to obtain a value
for the normalized shear viscosity. Multiplying this value by
4/3, and assuming that the bulk viscosity is negligible, we then
obtain the normalized longitudinal viscosity, .
The nondissipative case is represented by the dashed blue
curves in the two gures, which conrm the discussion of
the undamped acoustic wave outlined above for the two
models under consideration. The continuous red curves show
the dispersion in the presence of viscosity. One notes that
the dissipation has a minimal effect on the long-wavelength
dispersion of the wave in both gures. We also see that both the
strongly and weakly coupled models have dispersion curves
of similar shape, with an acoustic mode for small

k, and then
a attening of the curve, followed by a reduction in real
frequency as

k is increased, until the zero-frequency mode
is found for all

k >

k
cr
. However, it is seen that the critical
wave number

k
cr
is signicantly larger for the strongly coupled
043103-4
DUST-ACOUSTIC SHOCKS IN STRONGLY COUPLED . . . PHYSICAL REVIEW E 89, 043103 (2014)
FIG. 1. (Color online) Real (positive, upper curves) and imagi-
nary (negative, lower curves) parts of the linear dispersion relation
for dust-acoustic waves in both (a) weakly and (b) strongly coupled
dusty plasmas. Here we have n
i0
= 5 10
14
m
3
, n
e0
= 3.3
10
14
m
3
, n
d0
= 1 10
11
m
3
, Z
d
= 1700, k
B
T
i
= 0.05 eV, and
k
B
T
e
= 2.5 eV.
model than for the weakly coupled case in Fig. 1(a), that is, the
nite frequency mode can be found over a wider wavelength
range, although it is increasingly damped as

k is increased
(dotted red curves).
B. Nonlinear shocks
To obtain Burgerss equation, we use the moving and
stretched reference frame such that

= ( x v

t ), =
2

t , (20)
where v is the linear phase velocity of the DAW or the sound
speed of the plasma. Then, to lowest order, the expansion of
Eqs. (6)(9) leads to the relations
v =
_
1 +d
0
+d
11
d
12
, (21)
n
1
=
1
, u
1
= v
1
. (22)
To next order, one may obtain Burgerss equation for this
system such that

1

+

A
1

=

C

2
, (23)
in which we have

A =
(1 +2v
2
+2c
2
+2 )
2v
,

C =

2
, (24)
where is dened in Eq. (18) and is dened as
= d
11
d
12
+d
23
d
24
+d
25
. (25)
Equation (23) admits shock solutions of the form

1
(

, ) =
m
_
1 tanh
_

__
, (26)
where =



U ,
m
=

U/

Ais the amplitude of the shock,

= 2

C/

U is the normalized width, and

U is the normalized
velocity of the nonlinear wave in the moving reference frame.
Given the form of Eq. (26), it is clear that both the
amplitudes and widths of the shocks are dependent on

U, the
speed of the shock relative to the sound speed. In particular,
in any specied plasma system, faster shocks will have larger
amplitudes and narrower widths than slower shocks. Also, it
is interesting to note that the quantity
m

is independent of
the speed of the shock in any particular system.
In a laboratory situation, one is likely to measure both the
sound speed and the shock speed directly. However, from a
theoretical perspective, an underestimation of the sound speed
in a dusty plasma may occur if one does not account for the
effects of strong coupling in the model. For example, for the
normalization we have chosen, the sound speed in a weakly
coupled dusty plasma system, that is, without including the
electrostatic temperature terms, is scaled to unity. Using the
typical plasma parameters displayed in Table I in Sec. V
(System A), the sound speed in our strongly coupled model is
predicted to be v 1.41. This then has to be taken into account
when choosing the value of

U for the shocks. Furthermore,
since the lower bound for the speed of shock waves in the
laboratory frame is v, this model predicts that shocks do not
occur with a laboratory speed V such that 1 V v, which
are predicted in a weakly coupled plasma.
In addition to

U, the amplitude of a shock wave is dependent
on the

A coefcient in Burgerss equation. In this paper, we
take a dynamical form of the electrostatic temperature, and
by using this we obtain the value of

A as dened in Eq. (24).
This is the same coefcient that was found in Ref. [18], in
which solitons were studied. If we instead had neglected the
perturbations and assumed T

= T
0
, we would obtain

A =
3 +2c
2
+2d
0
2

1 +d
0
, (27)
15 10 5 0 5 10 15
0.2
0.15
0.1
0.05
0
/
D0

0
T
*
T
*
(n
d
,)
T
*
T
*0
T
*
=0
FIG. 2. (Color online) Shock solution to Burgerss equation for
the three different models. Here we have n
i0
= 5 10
14
m
3
, n
e0
=
3.3 10
14
m
3
, n
d0
= 1 10
11
m
3
, Z
d
= 1700, k
B
T
i
= 0.05 eV,
k
B
T
e
= 2.5 eV, 0.23, and

U = 0.1.
043103-5
S. E. COUSENS et al. PHYSICAL REVIEW E 89, 043103 (2014)
15 10 5 0 5 10 15
0.12
0.1
0.08
0.06
0.04
0.02
0
/
D0

0
= 0.1
0.23
= 0.5
FIG. 3. (Color online) Shock solution to Burgerss equation for
different dissipation coefcients. Here we have n
i0
= 5 10
14
m
3
,
n
e0
= 3.3 10
14
m
3
, n
d0
= 1 10
11
m
3
, Z
d
= 1700, k
B
T
i
=
0.05 eV, k
B
T
e
= 2.5 eV, and

U = 0.1.
and if we neglect the electrostatic temperature terms com-
pletely, we obtain

A =
3 +2c
2
2
. (28)
For the typical dusty plasma parameters we have used, there is
a signicant difference in the amplitude of the shocks between
the three models, as shown in Fig. 2.
The normalized widths of the shocks are only dependent
on the normalized viscosity coefcient , such that higher
viscocities lead to wider shocks, which is consistent across all
three models. Shock solutions for typical plasma parameters
and various values for are shown in Fig. 3. This suggests
that measuring the thickness of shocks in a plasma may be a
possible method of determining the viscosity of the plasma.
V. NUMERICAL INVESTIGATION
In this section, we present a simple numerical investigation
of the evolution of dust-acoustic shocks propagating in plasma
TABLE I. Parameters for the two regions in this numerical study.
Parameter Region A Region B
Plasma parameters
Z
d
1700 1400
m
d
(kg) 1 10
15
1 10
15
n
d0
(m
3
) 1 10
11
1.2 10
11
n
i0
(m
3
) 5 10
14
2.625 10
14
n
e0
(m
3
) 3.3 10
14
9.45 10
13
k
B
T
i
(eV) 0.05 0.05
k
B
T
e
(eV) 2.5 2.5
Viscosity calculations
390 280

0
2.92 1.98
0.23 0.31
Burgerss equation coefcients
A 1.82 1.51
C 0.11 0.16
Initial conditions

U 0.1 0.1

m
0.055 0.066
2.27 3.15
20 10 0 10 20 30 40
0.12
0.1
0.08
0.06
0.04
0.02
0
/
D0

0
= 0
= 50
=100
FIG. 4. (Color online) Shock originating in Region A, propagat-
ing in Region B.
conditions that are different to those in which they formed.
Physically, these changing conditions may occur as the shock
propagates along the discharge chamber, for example.
Again, we use the parameters presented by Rosenberg
et al. [20], but we also have an additional region of higher
dust density. For notational convenience, we denote the set of
parameters fromRosenberg et al. [20] Region Aand the region
of greater dust density Region B. We then use a Runge-Kutta
method to numerically integrate Burgerss equation for Region
B using the shock solution of Region A as an initial condition
and vice versa. The parameters for these two regions, as well
as the initial conditions used in the simulations, are presented
in Table I. The normalized longitudinal viscosity, , was
determined for each system using the method described in
Sec. IVA. The various scaling parameters,
D0
and
pd
, for
example, differ for each region since the plasma parameters
FIG. 5. (Color online) Contour plots of a shock which is a
solution to Burgerss equation for Region A. Plot (a) shows that
the shock is stable in Region A, whereas plot (b) shows that if this
shock is propagating in Region B its width increases with time.
043103-6
DUST-ACOUSTIC SHOCKS IN STRONGLY COUPLED . . . PHYSICAL REVIEW E 89, 043103 (2014)
20 10 0 10 20 30 40
0.3
0.2
0.1
0
/
D0

0
= 0
= 50
= 100
FIG. 6. (Color online) Shock originating in Region B, propagat-
ing in Region A.
have been modied. To directly compare the shocks arising
fromthese two sets of parameters, it is more convenient to rst
convert them to the same normalizations. For the simulations
in this section, we have converted the relevent parameters
for Region B into the scaling of Region A, such that the
shock parameters for Region B become
m
0.13 and
4.35, while its Burgerss coefcients become

A 1.00
and

C 0.27.
First, we consider the propagation of a dust-acoustic shock
that is initially a solution to the Burgerss equation of Region A,
propagating in Region B. The increase in the number density
of dust particles in the plasma for Region B has the effect of
reducing
0
due to both a decrease in the interparticle spacing
and an increase in the screening length since both the electron
and ion densities are reduced. For the calculated value of for
each region, this reduction of
0
results in a signicant increase
in the viscosity of the plasma and thus also an increase in the
width of the shock solutions.
Figure 4 shows the result of this numerical experiment,
showing an increase in the width of the shock over time. The
speed of the shock is also reduced, which may be seen by a
comparison of Figs. 5(a) and 5(b).
Conversely, a shock that is a solution to Region B
propagating in Region A is found to experience a reduction in
its width over time, as illustrated in Fig. 6. Furthermore, the
speed of the shock increases when it propagates in Region A.
This is shown in Figs. 7(a) and 7(b).
VI. CONCLUSION
To summarize, we have theoretically investigated dust-
acoustic shocks in a dusty plasma using the phenomenological
approach of Gozadinos et al. [10], which was originally
developed for crystalline plasmas but has since been applied
successfully in the study of strongly coupled plasmas near
the liquid-crystal phase transition. For analytical tractability,
we have assumed weak-amplitude excitations and proceeded
by adopting a perturbative technique. Although technically
limited to small amplitudes, the qualitative results obtained in
our investigation can be extended to larger amplitudes (e.g., as
initial conditions for numerical investigations) and may thus
provide insight on the salient features of shock structures in
strongly coupled dusty plasmas. Based on the linear dispersion
relation and the weakly nonlinear Burgerss equation, which
includes dissipative effects associated with viscosity, for the
FIG. 7. (Color online) Contour plots of a shock which is a
solution to Burgerss equation for Region B. Plot (a) shows that
the shock is stable in Region B, whereas plot (b) shows that if this
shock is propagating in Region A its width decreases with time.
typical parameter values that we have used, the results from
this model suggest the following, in comparison to the weakly
coupled model:
Linear dust-acoustic waves:
(i) Strong coupling mitigates the effects of dissipation in
the dusty plasma such that larger wave numbers and higher-
frequency waves may exist.
(ii) In the long-wavelength limit, strong coupling leads to
an increase of the sound speed in the plasma.
Weakly nonlinear dust-acoustic shocks:
(i) The increase in the sound speed due to strong coupling
implies a larger minimum velocity for observable dust-
acoustic shocks.
(ii) For a specied shock speed in the laboratory frame,
the sound speed excess,

U, is reduced due to the increase in
predicted sound speed.
(iii) For a specied sound speed excess,

U, strong coupling
reduces the amplitude of the predicted shocks due to an
increase in the magnitude of the nonlinearity coefcient,

A.
(iv) In this model, the width of the shocks is not altered by
strong-coupling effects. The width of the shocks is, however,
directly proportional to the viscosity coefcient, . This
suggests that measurements of the shock thickness, together
with the simulations, may be employed to obtain more accurate
values of the viscosity of a dusty plasma.
Note added in proof. Recently, we were alerted to a
related article, namely, M. Shahmansouri and M. Rezaei [36].
Although the article is based on a similar model to ours, we
note that the similarity is limited to the algebraic structures.
Neither the way that article proceeds, nor the actual research
outcomes, overlap our work here. In addition, we note that they
have apparently used an alternative denition of the viscosity
043103-7
S. E. COUSENS et al. PHYSICAL REVIEW E 89, 043103 (2014)
coefcient , which has different dimensions to those of the
viscosity used here and elsewhere, e.g., Refs. [15,3335].
ACKNOWLEDGMENTS
M.A.H. thanks the National Research Foundation of South
Africa for partial support under Grant No. UID 85408, and
acknowledges that opinions, ndings, and conclusions or
recommendations expressed in this publication are those of
the authors, and that the NRF accepts no liability whatsoever
in this regard.
APPENDIX: ELECTROSTATIC TEMPERATURE
COEFFICIENTS
We take the expansion of the normalized electrostatic
temperature up to second order, such that
d = d
0
+d
1
+
2
d
2
, (A1)
where the equilibriumelectrostatic temperature, in normalized
units, is given as
d
0
= T
0
/T
0
, (A2)
with
k
B
T
0
=
N
nn
Z
2
d
e
2
12
0
3

n
d0
(1 +
0
) exp(
0
), (A3)
where
0
= 1/
3

n
d0

D0
.
The perturbations d
1
and d
2
are
d
1
= d
11
n
1
+d
12

1
(A4)
and
d
2
= d
21
n
2
+d
22

2
+d
23
n
2
1
+d
24
n
1

1
+d
25

2
1
, (A5)
d
11
= d
21
=
d
0
3
1 +
0
+
2
0
1 +
0
, (A6)
d
12
= d
22
= d
0
c
2

2
0
1 +
0
, (A7)
d
23
=
d
0
18

3
0
3
2
0
2
0
2
1 +
0
, (A8)
d
24
=
d
0
3
c
2

2
0
(
0
1)
(1 +
0
)
, (A9)
d
25
=
d
0
2
_
3c
3
c
2
2

0
_
2
0
1 +
0
, (A10)
where c
2
is dened by Eq. (10), and
c
3
=
(1 )
2
(1 +
3
)
6(1 +)
3
. (A11)
[1] R. Z. Sagdeev, Reviews of Plasma Physics, Vol. 4, edited by
M. A. Leontovich (Consultants Bureau, New York, 1966),
p. 23.
[2] H. Washimi and T. Taniuti, Phys. Rev. Lett. 17, 996 (1966).
[3] P. K. Shukla and A. A. Mamun, Introduction to Dusty Plasma
Physics (Institute of Physics, Bristol, 2002).
[4] F. Verheest, Waves in Dusty Space Plasmas (Kluwer Academic,
Dordrecht, 2000).
[5] R. L. Merlino, Plasma Physics Applied, edited by C. Grabbe
(Transworld Research, Kerala, 2006), p. 74.
[6] N. N. Rao, P. K. Shukla, and M. Y. Yu, Planet. Space Sci. 38,
543 (1990).
[7] P. K. Shukla and V. P. Silin, Phys. Scripta 45, 508
(1992).
[8] A. Barkan, R. L. Merlino, and N. DAngelo, Phys. Plasmas 2,
3563 (1995).
[9] H. Ikezi, Phys. Fluids 29, 1764 (1986).
[10] G. Gozadinos, A. V. Ivlev, and J. P. Boeuf, New J. Phys. 5, 32
(2003).
[11] V. V. Yaroshenko, F. Verheest, H. M. Thomas, and G. E. Morll,
New J. Phys. 11, 073013 (2009).
[12] V. V. Yaroshenko, V. Nosenko, M. A. Hellberg, F. Verheest,
H. M. Thomas, and G. E. Morll, New J. Phys. 12, 073038
(2010).
[13] V. V. Yaroshenko, M. H. Thoma, H. M. Thomas, and G. E.
Morll, IEEE Trans. Plasma Sci. 38, 869 (2010).
[14] V. V. Yaroshenko, V. Nosenko, and G. E. Morll, Phys. Plasmas
17, 103709 (2010).
[15] A. A. Mamun and P. K. Shukla, Europhys. Lett. 87, 55001
(2009).
[16] A. A. Mamun, K. S. Ashra, and P. K. Shukla, Phys. Rev. E 82,
026405 (2010).
[17] V. Nosenko, S. K. Zhdanov, S.-H. Kim, J. Heinrich, R. L.
Merlino, and G. E. Morll, Europhys. Lett. 88, 65001 (2009).
[18] S. E. Cousens, S. Sultana, I. Kourakis, V. V. Yaroshenko,
F. Verheest, and M. A. Hellberg, Phys. Rev. E 86, 066404
(2012).
[19] S. Khrapak, D. Samsonov, G. Morll, H. Thomas,
V. Yaroshenko, H. Rothermel, T. Hagl, V. Fortov, A. Nefedov,
V. Molotkov, O. Petrov, A. Lipaev, A. Ivanov, and Y. Baturin,
Phys. Plasmas 10, 1 (2003).
[20] M. Rosenberg, E. Thomas, Jr., and R. L. Merlino, Phys. Plasmas
15, 073701 (2008).
[21] P. K. Shukla, Phys. Plasmas 8, 1791 (2001).
[22] V. E. Fortov, A. V. Ivlev, S. A. Khrapak, A. G. Khrapak, and
G. E. Morll, Phys. Rep. 421, 1 (2005).
[23] U. Konopka, G. E. Morll, and L. Ratke, Phys. Rev. Lett. 84,
891 (2000).
[24] A. Budea, A. Derzsi, P. Hartmann, and Z. Donko, Contrib.
Plasma Phys. 52, 194 (2012).
[25] J. P. Hansen and I. McDonald, Theory of Simple Liquids
(Academic Press, New York, 1986).
043103-8
DUST-ACOUSTIC SHOCKS IN STRONGLY COUPLED . . . PHYSICAL REVIEW E 89, 043103 (2014)
[26] S. Ichimaru and S. Tanaka, Phys. Rev. Lett. 56, 2815
(1986).
[27] M. S. Murillo, Phys. Rev. E 62, 4115 (2000).
[28] K. Y. Sanbonmatsu and M. S. Murillo, Phys. Rev. Lett. 86, 1215
(2001).
[29] G. Salin and J.-M. Caillol, Phys. Rev. Lett. 88, 065002 (2002).
[30] T. Saigo and S. Hamaguchi, Phys. Plasmas 9, 1210 (2002).
[31] N. DAngelo, Phys. Plasmas 4, 3422 (1997).
[32] N. DAngelo, Phys. Plasmas 5, 3155 (1998).
[33] P. K. Kaw and A. Sen, Phys. Plasmas 5, 3552 (1998).
[34] A. A. Mamun, P. K. Shukla, and T. Farid, Phys. Plasmas 7, 2329
(2000).
[35] P. K. Shukla and A. A. Mamun, IEEE Trans. Plasma Sci. 29,
221 (2001).
[36] M. Shahmansouri and M. Rezaei, Astrophys. Space Sci. 351,
197 (2014).
043103-9