# Péclet-number dependence of optimal mixing strategies identified using
multiscale norms

Conor Heffernan  Department of Applied Mathematics and Theoretical Physics,
Wilberforce Road, Cambridge, CB3 0WA, Cambridgeshire, United Kingdom  Colm-
cille P. Caulfield  BP Institute, Bullard Laboratories, Madingley Road,
Cambridge, CB3 0EZ, Cambridgeshire, United Kingdom

###### Abstract

The optimization of the mixing of a passive scalar at finite Péclet number
$Pe=Uh/\kappa$ (where $U,h$ are characteristic velocity and length scales and
$\kappa$ is the scalar diffusivity) is relevant to many significant flow
challenges across science and engineering. While much work has focused on
identifying flow structures conducive to mixing for flows with various values
of $Pe$, there has been relatively little attention paid to how the underlying
structure of initial scalar distribution affects the mixing achieved. In this
study we focus on two problems of interest investigating this issue. Our
methods employ a nonlinear direct-adjoint looping (DAL) method to compute
fluid velocity fields which optimize a multiscale norm (representing the
‘mixedness’ of our scalar) at a finite target time. First, we investigate how
the structure of optimal initial velocity perturbations and the subsequent
mixing changes between initially rectilinear ‘stripes’ of scalar and disc-like
‘drops’. We find that the ensuing stirring of the initial velocity
perturbations varies considerably depending on the geometry of the initial
scalar distribution. Secondly, we examine the case of lattices of multiple
initial ‘drops’ of scalar and investigate how the structure of optimal
perturbations varies with appropriately scaled Péclet number defined in terms
of the drop scale rather than the domain scale. We find that the
characteristic structure of the optimal initial velocity perturbation we
observe for a single drop is upheld as the number of drops and $Pe$ increase.
However, the characteristic vortex structure  and associated mixing exhibits
some nonlocal variability, suggesting that rescaling to a local $Pe$ will not
capture all the significant flow dynamics.

###### keywords:

keyword one , keyword two

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

††journal: Nonlinear Physics D

Direct-adjoint-looping method used to compute fluid flows to optimize mixing.

Questions of interest relate to the structure of the underlying passive
scalar.

Qualitative differences between a rectilinear and disc geometry are found.

The disc geometry problem is scaled up and we find the base structure is
preserved.

##  1 Introduction

A robust understanding of fluid mixing is of huge importance in the study of
various physical phenomena. Irreversible scalar mixing is characterized by a
relatively large scale ‘stirring’ combined with smaller scale diffusion to
homogenize an initially inhomogeneous distribution of a scalar field [1]. One
of the many difficulties of studying irreversible mixing lies in the fact that
there is no unified theory that rigorously defines it. For example, an
interesting question from a mathematical viewpoint is how exactly does one
define the ‘mixedness’ of a substance? Intuitively, the variance of a mean
zero passive scalar $\theta$ on a domain $\Omega$, defined as

$\displaystyle||\theta(\cdot,t)||^{2}_{L^{2}}$
$\displaystyle=\frac{1}{\mu(\Omega)}\int_{\Omega}|\theta(x,t)|^{2}\,dx,$ (1)

where $\mu$ denotes Lebesgue measure (area or volume) is a natural measure of
the mixedness of a scalar field. In the particularly simple situation when the
domain is a 2D torus $\Omega=\mathbb{T}^{2}$, it can be shown that Equation 1
reduces to

$\displaystyle||\theta(\cdot,t)||^{2}_{L^{2}}$
$\displaystyle=\sum_{\mathbf{k}}|\hat{\theta}_{\mathbf{k}}|^{2},$ (2)

where $\hat{\theta}_{\mathbf{k}}$ are the Fourier coefficients of $\theta$ and
the sum is taken over all wavenumbers $\mathbf{k}$.

While this measure is intuitive, there are non-trivial issues which need to be
addressed if the scalar field satisfies an advection-diffusion equation with
appropriate boundary conditions and a divergence-free velocity field
${\mathbf{u}}$ defined with periodic boundary conditions with spatial period
$L$ defined at points $(x,y)\in\Omega$

$\displaystyle\frac{\partial\theta}{\partial t}+\mathbf{u}\cdot\nabla\theta$
$\displaystyle=\frac{1}{Pe}\nabla^{2}\theta,$ (3)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,$ (4)
$\displaystyle\theta(x,y,t)$ $\displaystyle=\theta(x+L,y,t),$ (5)
$\displaystyle\theta(x,y,t)$ $\displaystyle=\theta(x,y+L,t).$ (6)

First in the limit when the diffusivity goes to zero, or equivalently when
$Pe\rightarrow\infty$, the variance is actually a conserved quantity. Even
when $Pe$ is finite yet large, the decay of the variance is generically
expected to be relatively slow, which introduces computational challenges to
flow optimization calculations [2]. A multiscale measure has been developed to
overcome this problem, based around the use of Sobolev norms of negative index
$-s$ [3]:

$\displaystyle||\theta(\cdot,t)||^{2}_{H^{-s}}$
$\displaystyle=\sum_{\mathbf{k}\neq\mathbf{0}}|\mathbf{k}|^{-2s}|\hat{\theta}_{\mathbf{k}}|^{2}.$
(7)

This measure is similar to Equation 2 but weights wavenumbers differently
(depending on $s$), with smaller wavenumbers (corresponding to larger scale
structures) leading to larger values of this measure. This measure is commonly
referred to as a ‘mix-norm’ [4], as if a scalar field is ‘mixed’ (in the
intuitive sense of having no large scale coherent structure) then the mix-norm
is expected to be small. A mix-norm is also advantageous as it is consistent
with the mathematically rigorous ergodic mixing and any choice of $s>0$ is
consistent with the theory [5].

Mixing quantification is a necessary component for mixing optimization subject
to a finite initial energy constraint. One successful method in particular has
been to use a cost functional-based approach and maximize the energy growth at
a target time $T$ [6, 7]. To solve this type of optimization problem a
‘direct-adjoint-looping’ (DAL) method has been developed [8, 9, 10]. The
algorithm computes initial flows $\mathbf{u}$ which optimize a given cost
functional at a prescribed target time $T$. It is advantageous in that it
enforces the fully nonlinear Navier-Stokes equations as a constraint of the
problem. A benchmark of its success has been established by producing optimal
perturbations which effectively homogenize a scalar field undergoing mixing.
Interestingly, optimal perturbations which minimize the mix-norm (with given
index) at relatively short target times can be excellent proxies for
perturbations which minimize variance at longer target times (and hence ensure
thorough mixing). Furthermore, for a given initial energy, such mix-norm-
minimizing perturbations are much more efficient at mixing than perturbations
which are chosen to maximize (perturbation) energy growth [2, 11, 12]. More
recently, the DAL method has also been used to investigate how variations in
the mix-norm index $s$ affects the ensuing mixing dynamics in a toroidal
geometry of the associated ‘optimal’ perturbations [13, 14].

The connections (and contrasts) between variance and mix-norm based strategies
has been investigated for flows with passive and dynamic (i.e. the flow is
density-stratified in a gravitational field) scalars [2, 11] as well as at
high and low values of finite Péclet number in three-dimensional flows [12].
All these studies have demonstrated that the mix-norm is indeed a useful proxy
for computing optimal mixing. A common theme of the previously considered
mixing optimization problems has been to consider relatively simple
rectilinear structures for the initial passive (or dynamic) scalar, typically
with one or two interfaces between areas of high and low concentration, with
$\theta\approx 1$ and $\theta\approx-1$ respectively such that the spatial
mean is zero. Such considerations were natural first steps in demonstrating
that the DAL method can identify optimal perturbations, and also that the use
of the mix-norm is appropriate and computationally efficient.

Now that we have confidence in the DAL method, in this paper we build on such
previous studies [2, 11, 12, 13, 14] to investigate mixing optimization for
initial scalar distributions with more detailed structure. We wish to gain
insight into the particular fluid-dynamical processes which lead to mixing for
different structures in the initial scalar distribution. In particular, we
attempt to answer two questions. First, we ask the question:

Q1: what differences may be found in mixing results for rectilinear ‘stripes’
of scalar as compared to disc-like ‘drops’?.

To answer this question, we perform a comparative study between an initial
scalar distribution of a rectilinear stripe used previously [13, 14] and an
initial scalar distribution of one disc-like drop. We find that there is a
qualitative difference between flows associated with these initial structures,
with the stripe-like structure undergoing more mixing compared with the drop-
like structure for the same initial perturbation energy. For both structures,
the optimal initial perturbation takes the form of multiple vortices that
localize along the interfaces (i.e. the locations where $\theta$ vary rapidly)
and hence enhance mixing. These vortices have a characteristic scale
significantly smaller than both the stripe and the (single) drop.

Motivated by the generic observation of vortices localized at interfaces, the
second question we ask is:

Q2: at what scale, if any, do these vortices cease to align along drop
interfaces for smaller and smaller drops and if they do cease to appear what
structures replace them?

In other words, if the initial structure is a lattice of sufficiently many and
hence sufficiently small drops, does the initial optimizing velocity
perturbation (in the form of interfacially-localized vortices) change from
what has been found in [2, 11, 12, 13, 14] and if so, what sort of mixing
dynamics and initial optimal perturbation(s) replace it? The distribution and
scales present in the mixing will be characterized by the number of drops in a
(square) lattice which will be given as $4^{n-1}$ for $n\in\\{1,2,...\\}$.
Indeed, an associated question is:

Q2a: Is it possible to identify an appropriate self-similar rescaling of the
flow dynamics, in particular in terms of a local Péclet number in terms of the
drop scale rather than the domain scale?

As we are always considering the flow on a two-dimensional torus, the square
lattices effectively allow us to consider whether non-local mixing processes
associated for example with broken symmetries in the mixing from drop to drop
can actually lead to enhanced mixing. To address these questions, the rest of
the paper is organised as follows. In section 2, we briefly review the DAL
method, and describe the various initial scalar distributions which we
consider. We then present and interpret our results in section 3, and finally
draw our conclusions in section 4.

##  2 Methods

There are four key control parameters for the mixing optimization problems
which we consider. Two are directly associated with the optimization problem:
the particular choice of the mix-norm index $s$ used in the objective
functional and the target time $T$ of the optimization problem over which the
mix-norm is minimized. The other two are associated with the properties of the
flow: the Péclet number $Pe$ and the Reynolds number $Re$, defined as

$Pe=\frac{Uh}{\kappa},\ Re=\frac{Uh}{\nu},$ (8)

where $U,h$ are characteristic velocity and length scales and $\nu,\kappa$ are
the kinematic viscosity and scalar diffusivity respectively. As we discuss
further below, when we consider initial lattices of disc-like drops of scalar,
care needs to be taken in the choice of the appropriate length scale $h$,
given in terms of the individual drop rather than the domain. This also has
implications for the appropriate definition of the target time $T$, as $Pe$
may also be interpreted as the ratio of the characteristic diffusive time
scale $h^{2}/\kappa$ to the advective time scale $h/U$. We return to this key
point below.

In [2, 11, 12, 14], the Schmidt number $Sc=\nu/\kappa$ was kept fixed at
unity, and so $Pe=Re$. To assess the different behaviour of flows with either
initial stripes or drops, and thus to answer the first question Q1 posed in
the Introduction, we will allow $Pe$ to differ from $Re$. For the second
question Q2, where we are interested in whether the mixing is inherently local
or exhibits meaningful differences between single drops and lattices of drops
(when appropriate rescaling of the key parameters is done), for simplicity we
fix $Re=Pe$. We also fix the initial perturbation energy density, non-
dimensionalized using the characteristic velocity scale $U$. Solutions of the
optimization algorithm corresponding to these parameters will be denoted by
$OA(s,T,Re,Pe)$, with subscripts to distinguish between the striped (SG) and
disc (DG) initial scalar distributions. For comparison between the striped and
disc geometry, we consider four different choices of $Re$ and $Pe$: $Re=100$
with either $Pe=50$ or $Pe=500$; and $Pe=100$ with either $Re=50$ or $Re=500$.
For either the single stripe or the single disc-like drop, we also vary the
index of the mix-norm $s=0.5,1,2$. For the various values of mix-norm index,
we choose the target time $T=s$. We make this choice due to our previous
observation in [14] that increasing either $T$ or $s$ leads to a qualitatively
similar change in both the structure of the optimal initial perturbation and
the subsequent time evolution of the vorticity distribution, in that
increasing either $T$ or $s$ favours larger characteristic vortices in the
initial optimal perturbation. Therefore, keeping the ratio between these two
control parameters leads to some simplification of the problem while still
retaining the key dynamics of interest. When considering the drop lattice
problem, we restrict attention to $s=1,T=0.1,0.2,0.4,0.8$ and
$Pe=400,800,1600,3200$.

We use the nonlinear DAL method [10] to compute the initial velocity field
which will minimize the value of the mix-norm at the prescribed target time
$T$. The flow takes place in a 2D torus of length $2\pi$ with $x$ and $y$
denoting the horizontal and vertical directions respectively. The velocity
field $\mathbf{u}=(u,v)$ and pressure $p$ are governed by the incompressible
Navier-Stokes equations and the passive scalar field $\theta$ is governed by a
conventional advection-diffusion equation. Therefore, the non-dimensionalized
equations governing the evolution of these variables are:

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}$ $\displaystyle=-\nabla
p+Re^{-1}\nabla^{2}\mathbf{u},$ (9) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=0,$ (10) $\displaystyle\frac{\partial\theta}{\partial
t}+\mathbf{u}\cdot{\nabla}\theta$ $\displaystyle=Pe^{-1}\nabla^{2}\theta,$
(11)

where $Pe$ and $Re$ denote the Péclet number and Reynolds number respectively
and are defined by Equation (8).

As the scalar field is passive, we may solve Equations 9 and 10 and Equation
11 separately. In the case of the striped geometry, we choose

$\displaystyle\theta_{SG}(\mathbf{x},0):=\theta_{0}(\mathbf{x})$
$\displaystyle=\tanh\left(6\left(x-\frac{\pi}{2}\right)\right)-\tanh\left(6\left(x-\frac{3\pi}{2}\right)\right)-1.$
(12)

This corresponds to a smooth zero-mean scalar distribution with a vertical
stripe, centred at $x=\pi$ of width $\pi$ of positive $\theta\simeq 1$,
bordered by stripes of negative $\theta\simeq-1$. For the disc ‘drop’
geometry, we choose

$\displaystyle\theta_{DG}^{n}(\mathbf{x},0):=\begin{cases}1,&\quad\text{ if
}\mathbf{x}\in C_{n}\\\ -1,&\quad\text{ if }\mathbf{x}\not\in
C_{n}\end{cases}$ (13)

where $C_{n}$ denotes the set of $4^{n-1}$ discs with centres evenly
distributed in the vertical and horizontal directions on a lattice and radius
given by

$\displaystyle r_{n}$ $\displaystyle=\sqrt{\frac{2\pi}{4^{n-1}}}.$ (14)

In our study, we consider $n=1,2,3,4.$ We also require an initial condition
for the velocity field $\mathbf{u_{0}}=\mathbf{u}(\mathbf{x},0)$ in order to
solve the system. For the very first loop of the DAL method, we set
$\mathbf{u_{0}}$ to be random noise. After each iteration of the loop, we
update $\mathbf{u_{0}}$ to give us our initial condition to evolve the system
of Equations 9, 10 and 11.

In order to identify the initial perturbation which optimizes the mixing of
the initially inhomogeneous fluid, we seek to minimize an objective functional
subject to a constraint on the kinetic energy of the initial perturbation:

$\displaystyle||\mathbf{u_{0}}||^{2}_{L^{2}(\Omega)}$
$\displaystyle=2e_{0}\mu(\Omega),$ (15)

where $e_{0}=0.03$ is the non-dimensional perturbation energy density and
$\mu$ denotes the ‘volume’ (i.e. the area) of the flow domain $\Omega$. As has
been previously demonstrated in [13, 14], this perturbation energy density is
chosen as an appropriate intermediate value. It is both sufficiently large to
allow for the identification of non-trivial initial flow structures, and yet
sufficiently small so that mixing is not too rapid and so there is still the
possibility to distinguish the efficacy for mixing of different initial
perturbation structures.

We define the objective functional as

$\displaystyle\mathcal{J}(\theta(T))=\frac{1}{2}||\theta(\mathbf{x},T)||_{H^{-s}(\Omega)}^{2},$
(16)

i.e. (half) the value of the Sobolev norm of (negative) index $-s$ (which we
refer to as the mix-norm of index $s$) at the target time $T$. We may then
define the constrained optimization problem of interest as

$\displaystyle\text{argmin }\mathcal{J}(\theta(T))\text{ subject to
}||\mathbf{u_{0}}||^{2}_{L^{2}(\Omega)}=2e_{0}\mu(\Omega),$ (17)

where $\\{\mathbf{u},\theta\\}$ solve the system of Equations 9, 10 and 11.
The initial condition $\mathbf{u}_{0}$ does not appear explicitly in the
objective functional, but nevertheless it affects $\mathcal{J}$ through the
evolution of the flow variables which are constrained by the system of
Equations 9, 10 and 11. These constraints are imposed of course by the use of
appropriate Lagrange multipliers, the spatially and temporally evolving so-
called adjoint variables denoted by
$\\{\mathbf{u}^{\dagger},p^{\dagger},\theta^{\dagger}\\}=\\{\mathbf{v},q,\eta\\}$.
This is explained in detail in (for example) [13] and we follow their approach
here. We may define a Lagrangian as

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{J}(\theta(T))-\Sigma_{I\in\\{NS,AD,C,IC\\}}\quad\mathcal{J}_{I},$
(18)

where

$\displaystyle\mathcal{J_{NS}}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\mathbf{v}\cdot\left(\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}+\nabla
p-Re^{-1}\nabla^{2}\mathbf{u}\right),$ (19) $\displaystyle\mathcal{J_{AD}}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\eta\left(\frac{\partial\theta}{\partial
t}+\mathbf{u}\cdot{\nabla}\theta-Pe^{-1}\nabla^{2}\theta\right),$ (20)
$\displaystyle\mathcal{J_{C}}$
$\displaystyle=\int_{0}^{T}\int_{\Omega}q\nabla\cdot\mathbf{u},$ (21)
$\displaystyle\mathcal{J_{IC}}$
$\displaystyle=\int_{\Omega}\mathbf{v_{0}}\cdot(\mathbf{u}(\mathbf{x},0)-\mathbf{u_{0}}).$
(22)

Variation with respect to adjoint variables yields Equations 9, 10 and 11.
Similarly, variation with respect to the direct variables
$\\{\mathbf{u},p,\theta\\}$ results in the so-called ‘adjoint Navier-Stokes’
equations

$\displaystyle\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{v}$ $\displaystyle=-\nabla
q-Re^{-1}\nabla^{2}{\mathbf{v}}+\eta\nabla\theta,$ (23)
$\displaystyle\nabla\cdot\mathbf{v}$ $\displaystyle=0,$ (24)
$\displaystyle\frac{\partial\eta}{\partial t}+\mathbf{u}\cdot\nabla\eta$
$\displaystyle=-Pe^{-1}\nabla^{2}\eta.$ (25)

At $t=T,0$ we also produce the following terminal and initial conditions

$\displaystyle\mathbf{v}(\mathbf{x},T)$ $\displaystyle=0,$ (26)
$\displaystyle\eta(\mathbf{x},T)$
$\displaystyle=\sum_{\mathbf{k}\neq\mathbf{0}}|\mathbf{k}|^{-2s}Re\\{\hat{\theta}_{\mathbf{k}}(T)\exp{(i\mathbf{k}\cdot\mathbf{x})}\\},$
(27) $\displaystyle\mathbf{v_{0}}$ $\displaystyle=\mathbf{v}(\mathbf{x},0),$
(28) $\displaystyle\nabla_{\mathbf{u_{0}}}\mathcal{L}$
$\displaystyle=\mathbf{v_{0}},$ (29)

where $\hat{\theta}_{\mathbf{k}}$ are the Fourier coefficients of $\theta$ and
$Re\\{\cdot\\}$ denotes taking the real part. Due to the negative diffusion
terms $-Re^{-1}\nabla^{2}\mathbf{v}$ and $-Pe^{-1}\nabla^{2}\eta$, Equations
23, 24 and 25 must be integrated backwards in time to avoid numerical
instability. These equations are then integrated backwards from $t=T$ to
$t=0$, thus forming a ‘direct-adjoint loop’. Using a numerical technique from
[15] and with $\mathbf{u_{0}}^{n}$ and $\mathbf{v_{0}}^{n}$ denoting the
direct and adjoint velocities at $t=0$ after $n$ loops of this direct-adjoint-
looping (DAL) method, the updated guess $\mathbf{u_{n+1}}$ can be calculated
by

$\mathbf{u_{0}}^{n+1}=\cos(\phi)\mathbf{u_{0}}^{n}+\sin(\phi)\mathbf{w_{0}}^{n},$

where $\mathbf{w_{0}}^{n}$ denotes the scaled (by the energy constraint)
adjoint velocity $\mathbf{v_{0}}^{n}$ projected onto the hypersurface
tangential to the energy hypersphere at $\mathbf{u_{0}}^{n}$, as described in
detail in [16]. The angle of rotation $\phi$ is calculated by using a
backtracking line search [17]. This looping procedure is repeated until
convergence has been reached as measured by the normalized residual $r$,
defined by

$r=\frac{||\nabla_{\mathbf{u_{0}}}\mathcal{L}^{\perp}||_{L^{2}(\Omega)}^{2}}{||\nabla_{\mathbf{u_{0}}}\mathcal{L}||_{L^{2}(\Omega)}^{2}},$

where the symbol $\perp$ denotes projection onto the hyperplane tangential to
the energy hypersurface. Since the energy is fixed, a small residual ($r\sim
O(10^{-3})$) implies the gradient can only change by varying its magnitude
which is not permissible due to the (explicitly imposed) energy constraint.

The direct and adjoint equations are solved with a $4^{th}$-order mixed Crank-
Nicholson Runge-Kutta scheme with incompressibility enforced through a
fractional step method [2, 18]. Simulations were performed using a
discretization of $N=256$ points in both the $x$ and $y$ directions (we use
$N=512$ for the $Pe=3200$ calculation in Section 3.2). For presenting the
results graphically, the mix-norm and variance are scaled by the evolution of
the purely diffusive passive scalar defined as

$\displaystyle M_{s}(t)$
$\displaystyle=\frac{||\theta(\mathbf{x},t)||^{2}_{H^{-s}(\Omega)}}{||\theta_{d}(\mathbf{x},t)||^{2}_{H^{-s}(\Omega)}},$
(30) $\displaystyle V(t)$
$\displaystyle=\frac{||\theta(\mathbf{x},t)||^{2}_{L^{2}(\Omega)}}{||\theta_{d}(\mathbf{x},t)||^{2}_{L^{2}(\Omega)}},$
(31)

where $\theta_{d}$ is the solution of Equation 11 but with the advective term
dropped, and the dependence of the (scaled) mix-norm on the index $s$ is
labelled by the subscript. Scaled in this way, it is possible to identify the
extent to which fluid motions, and hence non-zero advection, affect the
reduction in mix-norm and variance over time.

##  3 Results

###  3.1 Striped Geometry (SG) versus Disc Geometry (DG) Mixing

In previous studies of nonlinear mix-norm optimization the Schmidt number is
often kept constant and so $Re=Pe$. In [14], we suggested varying both Péclet
number and Reynolds number individually rather than keeping them equal, which
we do here. Initially, we wish to address Q1, and so we wish to investigate
how the optimal perturbations change in structure and mixing performance
depending on the geometry of the passive scalar’s initial distribution.

Figure 1: Time evolution of scaled mix-norm $M_{s}(t)$ and scaled variance
$V_{s}(t)$, as defined in Equations 30 and 31, for flows initially seeded with
perturbations that minimize mix-norms (for a variety of values of index $s$)
with $s=T=0.5$ shown in red, $s=T=1$ shown in blue, and $s=T=2$ shown in
black. The solid lines denote the ‘SG’ flows which are optimized to mix the
striped geometry and the dashed lines denote the ‘DG’ flows optimized to mix
the disc geometry. Figures in the left-hand column correspond to flows with
varying $Pe$ and fixed $Re=100$ while figures in the right-hand column
correspond to flows with varying $Re$ and fixed $Pe=100$.

Figure 2: Time evolution of passive scalar and vorticity for DG flows for
$s=T=2$ and $Re=100,Pe=50$ (left columns) and $Re=500,Pe=100$ (right columns).
Snapshots are taken at $t=0,5,10,20$.

We consider two such geometries here. The first is the striped geometry (SG)
which consists of three layers, one on top of the other. The central stripe
indicates the region where $\theta=1$ and the two side stripes indicate
$\theta=-1$. The second disc geometry (DG) we consider is a disc-like ‘drop’
with radius $r=\sqrt{2\pi}$ (to ensure zero mean). The interior of the disc
corresponds to $\theta=1$ and the surrounding region corresponds to
$\theta=-1$.

Figure 1 show the evolution of scaled mix-norm and variance respectively with
solid lines denoting the striped geometry (SG) flows and dashed lines denoting
the disc (or drop) geometry (DG) flows.

Figure 3: Initial vorticity distribution of the optimal perturbations
identified for mix-norm minimization for flows with $Pe=50$ and $Re\equiv
100$. We keep $s=T$ and show $s=0.5$ (left column), $s=1$ (centre column) and
$s=2$ (right column). Top row denotes DG and bottom row denotes SG.

For flows with varying $Pe$ and fixed $Re$ (left-hand column in Figure 1,
there is a qualitative difference between SG and DG flows. Similarly to the
perturbations in [14], the scaled mix-norms decay rapidly and then approach
pure diffusion (and so the curves become close to horizontal). The degree of
decay depends on $s$ and $T$ and also the initial scalar distribution. SG
flows (solid lines) typically exhibit a larger initial decay yet a smaller
continued decay compared to the DG flows (dashed lines). This behaviour is
consistent with the case of varying $Re$ and fixed $Pe$ shown in the right-
hand column. In both cases of $Re$, the more rapid initial decay means that SG
flows typically mix more thoroughly than DG flows at both target time and
beyond.

Figure 4: Initial vorticity distribution of the optimal perturbations
identified for mix-norm minimization for flows with $Pe=500$ and $Re\equiv
100$. We keep $s=T$ and show $s=0.5$ (left column), $s=1$ (centre column) and
$s=2$ (right column). Top row denotes DG and bottom row denotes SG.

Figure 5: Initial vorticity distribution of the optimal perturbations
identified for mix-norm minimization for flows with $Re=50$ and $Pe\equiv
100$. We keep $s=T$ and show $s=0.5$ (left column), $s=1$ (centre column) and
$s=2$ (right column). Top row denotes DG and bottom row denotes SG.

The qualitative difference between high and low Péclet number cases,
conjectured in [14], appears actually to be more dependent on variations in
the Reynolds number for SG flows (solid lines in Figure 1), as it is clear
that the difference between flows with $Pe=50$ and $Pe=500$ for $Re=100$ (left
column) is not as marked as the difference between flows with $Re=50$ and
$Re=500$ for $Pe=100$ . However, there does appear to be a more significant
qualitative difference with $Pe$ for the DG flows. The scaled variance for the
DG flows (dashed lines) decays markedly more slowly than for the SG flows when
$Pe=500$, whereas for the lower $Pe$ case both flows undergo a similar
evolution.

Figure 6: Initial vorticity distribution of the optimal perturbations
identified for mix-norm minimization for flows with $Re=500$ and $Pe\equiv
100$. We keep $s=T$ and show $s=0.5$ (left column), $s=1$ (centre column) and
$s=2$ (right column). Top row denotes DG and bottom row denotes SG.

We now attempt to understand the differences in mixing behaviours of SG flows
and DG flows, for the various control parameters we have considered. We first
examine the differences between a low $Re$ and a high $Re$ DG flow. In Figure
2, we plot the time evolution of the vorticity and scalar fields for
$OA_{DG}(2,2,100,500)$ and $OA_{DG}(2,2,50,100)$. The initial structures
consist of elongated vortices along the interface. The difference between
these structures is that at higher $Re$ the vortices have perturbed corners
similar to those observed in [14] at the same index and Reynolds number. We
see from their evolution in Figure 2 that the lower $Re$ field dissipates its
energy more quickly; the vortices turn on their axis but retain a more rigid
structure. However, in the case of the higher Reynolds number the vortices
move freely throughout the domain with slower energy dissipation. Higher $Re$
flows more effectively mix the passive scalar as they are successful in fully
tearing apart the disc whereas lower $Re$ flows do not do this as successfully
and rely more on diffusion for full homogenization. As can be seen from the
top rows of Figures 3, 4, 5 and 6, increasing target time $T$ increases the
characteristic length scales of these vortices, and the perturbed corners are
more strongly associated with higher $Re$ than with higher $Pe$. As identified
in [14], these perturbed vortices seem to lead to better mixing which also
appears to hold here. This implies that at an intermediate value of mix-norm
index $s\sim 1-2$ with high $Re$ this initial structure brings about a better
mixture.

We now analyze the perturbations associated with the SG flows. Analogously to
above, in Figure 7 we plot the time evolution of the vorticity and scalar
fields for $OA_{SG}(2,2,100,50)$ and $OA_{SG}(2,2,500,100)$. Again, the
initial perturbation structure consists of vortices aligned along the
interfaces of the three layers. Similarly to the DG flow, we observe that at
low $Re$ the initial vortices do not undergo much change in their position as
the flow evolves in time. The scale set by the mix-norm manages to break the
central stripe effectively which leads to a good homogenization of the passive
scalar by $t=20$. At higher $Re$ the perturbations manage to break the central
stripe in the same time scale. However, the vortices are not fixed in position
and cascade into different scales while moving around the domain. The energy
dissipation in this case occurs over a longer timescale but both generate good
mixing. Furthermore, as is apparent from the bottom rows of Figures 3, 4, 5
and 6, increasing target time $T$ again increases the characteristic length
scales of these vortices. Also, increasing $Re$ once again appears to lead to
more changes in the fine structure of these vortices than increasing $Pe$
does.

Fundamentally, increasing $Re$ in both DG and SG flows causes a distortion in
the otherwise smooth shape of vortices. This distortion then manifests itself
by a cascade and free movement throughout the domain. A similar comparison can
be made at low $Re$ between the geometries as the vortices stay fixed and
dissipate energy more quickly. There is, however a distinct difference between
the two classes of flow. As shown in Figure 1, it is clear that the stirring
action in SG flows is superior to that which occurs in DG flows. A simple
explanation may be framed in terms of the characteristic length scale
($d_{c}$) between regions where initially $\theta\simeq-1$ either side of the
‘central’ region where $\theta\simeq 1$. For SG flows, $d_{c}\simeq\pi$, the
width of the central ‘stripe’ where $\theta=1$, while for (single drop) DG
flows this length is nontrivially larger, as $d_{c}\simeq
2r_{1}=\sqrt{8\pi}\approx 5.0133$, the initial drop diameter.

Figure 7: Time evolution of passive scalar (left) and vorticity (right) for SG
for $s=T=2$ and $Re=100,Pe=50$ (left two columns) and $Re=500,Pe=100$ (right
two columns). Snapshots are taken at $t=0,5,10,20$.

As was observed in [14], different parameter choices lead to different
symmetry/asymmetry structure and this does appear to have a role in how the
central region (with $\theta\simeq 1$) is distorted. However, the difference
in how symmetry is characterized is not the same for SG flows and DG flows.
Symmetry in SG flows is characterized by vortices on the rectilinear
(vertical) interfaces between the different stripes, whereas for DG flows
there is rotational symmetry around the edge of the central drop.

###  3.2 The Drop Lattice Problem

We now turn our attention solely to the DG flow, and analyze the scenario
where there is not one ‘drop’ but many, arranged in a square lattice. The size
and position of these drops will be given by Equation 13 (with mean zero). The
problem of interest is to investigate the mixing dynamics as $n$ increases and
investigate if increasing the number of drops leads to an appropriately self-
similar scaled version of the previously considered DG flow with a single
drop.

Periodic boundary conditions on one drop (as considered in Section 3.1) imply
an infinite lattice of drops next to one another with imposed symmetry. The
question of interest is to investigate the behaviour if we increase the number
of drops and decrease their size within the same $\Omega=[0,2\pi]^{2}$ domain.
We expect that as we scale up the number of drops we observe the same
structure around the edge of the disc for each of the smaller counterparts
(perhaps rotated around the interface). The central question is whether the
optimization algorithm converges onto other solutions which diverge from the
$n=1$ case, allowing for non-local interactions associated with broken
symmetries between neighbouring drops and how the behaviour of these solutions
compare qualitatively with one another. It is of interest to study how the
potential introduction of finer scales set by the Péclet number can change the
dynamics of the drop homogenization by comparing behaviour at different times,
when the time is scaled appropriately, as discussed below.

Figure 8: Time evolution of scaled mix-norm $M_{s}(t)$ and variance
$V_{s}(t)$, as defined in Equation 30 and Equation 31 respectively, for flows
initially seeded with perturbations that minimize mix-norms corresponding to
$n=1,2,3,4$ for the drop lattice problem.

Figure 9: Initial vorticity distribution of the optimal perturbations
identified for mix-norm minimization for flows with $4^{n-1}$ initial drops,
$Pe=50\times 2^{n-1}$ and target times $T_{n}=0.8\times 2^{-(n-1)}$ for
$n=1,2,3,4$.

For this problem we fix $s=1$ and we choose $T$ to scale appropriately (since
we are changing scales). For the flow with a single drop, the characteristic
length scale $h$ used in the definitions of $Pe$ and $Re$ in Equation 8 is
appropriately related to the domain scale. However, as more drops are added,
the appropriate length scale for a lattice of $4^{n-1}$ drops is
$h_{n}=h/2^{n-1}$. Assuming that the velocity scale is kept constant, (and
that the physical properties of the fluid $\nu$ and $\kappa$ are also kept
constant) we then expect the local Péclet number and Reynolds number for an
individual disc also to be smaller by the factor
$h_{n}/h=2^{-(n-1)},n=1,2,3,4$.

We set for a single drop $Pe=400$ and $T=0.8$ and then we increase and
decrease respectively by a factor of $2$ as we increase from $n=1$ to $n=4$
(corresponding to 1 drop and 64 drops respectively). In problems of this type
and those of [2, 11, 12, 13, 14, 16], the primary focus of study is the
stirring action due to the optimal mixing strategies. Therefore, the decision
to scale in this manner is to investigate dynamics on the local advective
timescale $t_{n}=h_{n}/U$ rather than the diffusive timescale
$t_{n}=h_{n}^{2}/\kappa$.

Underlying this rescaling is the assumption that the dynamics are entirely
local, and that the key dynamics are dominated by advection. Therefore, there
are two ways in which there could be different behaviour even for the scaled
cases. First, there could be nonlocal effects due to broken symmetry for
example which could modify the mixing dynamics. Secondly, particularly at
later times, when the flow is dominated by diffusive processes, the diffusive
time scale could become dominant. However, in this study, we choose to focus
solely on the former.

Following this scaling approach, we choose $T=T_{4}=0.1$ for the $n=4$ case
(with 64 drops) and increase by a factor of 2 as we go down through $n=3,2,1$
so that for the base case (with one drop) we have $T_{1}=0.8$. In the left
column of Figure 8 we plot the scaled mix-norm and variance for the evolution
of the passive scalar with $n=1,2,3,4$ (shown in Figure 9) under the action of
its corresponding optimal flow field. We also plot in the right column of
Figure 8 shifted mix-norm and variance with the time scale ‘stretched’. The
behaviour we observe at the $Pe=400$ solution happens over a longer timescale
than in the case of $Pe=3200$ and so we need to rescale in order to compare
dynamics. As can be seen from Figure 8, there is a slight discrepancy in that
the $n=4$ case undergoes a much quicker decay. An explanation of this lies in
the examination of the vorticity of the solutions shown in Figure 10. We
observe that as we increase the number of drops in our computational domain we
do not get a repeat of the $n=1$ structure. This change becomes clear passing
through $n=2,3$ and is fully realized when $n=4$. As can be seen from Figure
10, the $n=1$ case is asymmetric and disordered which results in a slower rate
of decay. The $n=2$ case retains the corners observed for $n=1$ but this
structure begins to break down into a more rigid structure as observed for
$n=3,4$. After this transitionary phase, we oberve a symmetric solution for
the $n=4$ case which leads to the most rapid decay. This would appear to
suggest that the algorithm ‘catches’ some smaller scale structures which cause
a slow drift away from the $n=1$ case to something qualitatively different at
$n=4$. Indeed, the structure appears to diverge away from what our base case
suggests.

The emergence of these finer scales leads to a more rapid homogenization
(particularly at $Pe=3200$) and thus we observe that (after a rescale of time)
we do not see the exact same mix-norm and variance decay between solutions. It
appears that as we increase $Pe$, the optimization algorithm exploits the
range of smaller scales now available as set by $Pe$ and uses them to perturb
the base case in order to lead to a better mixture. Interestingly, these
smaller scale structures are dominated by advective processes, as at later
times the (rescaled) $n=4$ case is observed actually to homogenize more
slowly, particularly relative to the purely diffusive behaviour (i.e. the
denominator in Equation 30). However, care needs to be taken not to over-
interpret the late time behaviour of flows significantly after the imposed
target time ($T=0.8$ on the rescaled plots), as it is only after that time
that the $n=4$ case starts (apparently) to under-perform relative to the other
cases.

We observe a very similar behaviour between solutions as Péclet number is
increased. While there are some finer scales present as we increase $Pe$, they
do appear to follow the same qualitative behaviour as the base case. We can
see that as time progresses, the vorticity field deforms the interface and
homogenizes the disc region into a thin stripe-like shape. This implies that
the mixing dynamics behave in the same way as we scale up, even with the
introduction of finer scales for higher $Pe$. Therefore, while there is some
slight divergence away from the base case for lower $Pe$, the solution
essentially behaves the same way as can be seen in the evolution of solutions
in Figure 10. In fact, for higher $Pe$ (and higher $Re$ since Schmidt number
is fixed at $Sc=1$) we see that the algorithm scales up the base behaviour
even as we approach a highly disordered regime. This suggests that as one
exits an essentially laminar regime and finer scales are introduced, the
underlying mixing dynamics still hold.

To conclude this section, we comment on the emergence of a symmetric structure
in the drop lattice at $Pe=3200$. In Figure 11 we show the evolution of the
passive scalar and vorticity. The evolution shows that the structure is
perfectly symmetric, with each local drop deforming in the same way. This is
in stark contrast to the lower $Pe$ cases which stretches and folds the drops
asymmetrically but maintains essentially the same qualitative homogenization
as shown in Figure 10. It appears that as we scale up, the finer scales tend
to impose a rigid structure on the drop lattice. We conjecture that as $Pe$
increases, the optimal mixing field favours a rigid symmetric structure.
However, if we compare this with the $Pe=800$ case, for example, we observe
that the lower number of drops present leads to an asymmetric lattice of
vortices. This leads to disordered mixing dynamics but as we can see from
Figure 8, this asymmetry both at unscaled and rescaled time decays more
rapidly than in the $Pe=3200$ solution which suggests that the symmetry may
not necessarily be preferable for optimal mixing.

Figure 10: Time evolution of passive scalar (left) and vorticity (right) for
optimal perturbations that minimize the mix-norm with $s=1$ for drop lattice
problems with $Pe=200,n=2,T=0.4$ (left columns) and $Pe=800,n=3,T=0.2$ (right
columns). Snapshots are taken at $t=0,5,10,20$.

Figure 11: Time evolution of passive scalar (top row) and vorticity (bottom
row) for the drop lattice problem with $Pe=3200,n=4,s=1,T=0.1$. Snapshots are
taken at $t=0,5,10,20$.

##  4 Conclusion

In this paper, we have studied two problems relating to the Péclet number
dependence of optimal mixing strategies. These problems were studied using the
nonlinear direct-adjoint-looping method [9, 10] to compute optimal flows which
minimize the mix-norm at a prescribed target time. The central theme of both
of the problems investigated concerned how optimal mixing relates to the
structure of the underlying passive scalar.

In the first problem, we compared optimal mixing on the geometry of a stripe,
as in [13, 14], with that of a single drop. We found similarities between both
regimes, such as the asymmetric alignment of vortices along the interfaces and
the scale of the vortices being set by the mix-norm index. However, we also
found discrepancies between the two geometries. The stripe geometry underwent
much more stirring than the disc geometry, though both homogenized reasonably
well at similar times. We also varied the Schmidt number and so consequently,
as suggested by [14], we had $Pe\neq Re$ in this problem. It would appear from
this study that the Reynolds number controls the scale whereas the Péclet
number controls the shape and symmetry of the vortices.

In the second problem (with $Pe=Re$), we examined the effect on mixing
dynamics of a lattice of many drops distributed throughout the computational
domain. The motivation for this problem was to determine whether the vortex
alignment structures observed so far in [2, 11, 12, 14] scale as $Pe$ is
increased accordingly with the number of drops present in the initial scalar.
Indeed, this turns out to be the case as one approaches the turbulent regime
but with the slight difference of smaller scales being introduced to perturb
the $n=1$ case but ultimately leading to the same types of dynamics. While the
structure of the vortices does indeed change as $n$ is increased, we did not
find a critical scale at which the alignment structure breaks down. Future
studies could look at increasing drops and $Pe$ until such a scale is found.

Looking to future directions, the stripe and disc geometries appear to provide
bounds for the stirring action and it may be of interest to investigate other
structures which may or may not sharpen these bounds. Also, as we found
different behaviour for the drop lattices, particularly for the $n=4$ case, it
is clear that careful consideration needs to be paid to ensure that small-
scale processes are accessible, particularly at higher Péclet number.

Contributions: All authors contributed equally to the present article.

Interests: We declare we have no competing interests.

Funding: This work was done as part of a PhD project made possible by the
National University of Ireland’s Travelling Studentship in the Sciences and
their support is gratefully acknowledged by the authors.

Acknowledgements: The authors wish to thank the two anonymous reviewers whose
comments have greatly improved the presentation and quality of this
manuscript.

Data Accessibility: Accompanying animations for the plots may be found here:
https://data.mendeley.com/datasets/n6ctnkmndx/draft?a=06be2dfb-662f-4e9b-96fb-a2f081a5ca8e

Source code is based on a modified version of the solver found in [14].

## References

  * [1] C. Eckhart, An analysis of the stirring and mixing processes in incompressible fluids., J. Mar. Res. 7 (1948) 265–275.
  * [2] D. P. G. Foures, C. P. Caulfield, P. J. Schmid, Optimal mixing in two-dimensional plane poiseuille flow at finite péclet number, J. Fluid Mech. 748 (2014) 241–277.
  * [3] G. Mathew, I. Mezić, L. Petzold, A multiscale measure for mixing, Physica D: Nonlinear Phenomena 211 (2005) 23–46.
  * [4] J.-L. Thiffeault, Using multiscale norms to quantify mixing and transport, Nonlinearity 25 (2012) R1.
  * [5] Z. Lin, J.-L. Thiffeault, C. Doering, Optimal stirring strategies for passive scalar mixing, J. Fluid Mech. 675 (2011) 465–476.
  * [6] O. Aamo, M. Krstic, Flow Control by Feedback: Stabilization and Mixing, Communications and Control Engineering, Springer London, 2013.
  * [7] O. M. Aamo, M. Krstić, T. R. Bewley, Control of mixing by boundary feedback in 2d channel flow, Autom. 39 (2003) 1597–1606.
  * [8] S. Cherubini, P. De Palma, J.-C. Robinet, A. Bottaro, Rapid path to transition via nonlinear localized optimal perturbations in a boundary-layer flow, Phys. Rev. E 82 (2010) 066–302.
  * [9] C. C. T. Pringle, R. R. Kerswell, Using nonlinear transient growth to construct the minimal seed for shear flow turbulence, Phys. Rev. Lett. 105 (2010) 154–502.
  * [10] R. Kerswell, Nonlinear nonmodal stability theory, Annual Review of Fluid Mechanics 50 (1) (2018) 319–345.
  * [11] F. Marcotte, C. P. Caulfield, Optimal mixing in two-dimensional stratified plane poiseuille flow at finite péclet and richardson numbers, J. Fluid Mech. 853 (2018) 359–385.
  * [12] L. Vermach, C. P. Caulfield, Optimal mixing in three-dimensional plane poiseuille flow at high péclet number, J. Fluid Mech. 850 (2018) 875–923.
  * [13] L. Vermach, Using multiscale norms for nonlinear mixing optimisation, PhD Thesis (2016).
  * [14] C. Heffernan, C. Caulfield, Robust and efficient identification of optimal mixing perturbations using proxy multiscale measures, Phil. Trans. A to appear (2022).
  * [15] S. Douglas, S. Amari, S. Kung, Gradient adaptation under unit-norm constraints, 1998, pp. 144–147, proceedings of the 1998 9th IEEE SP Workshop on Statistical Signal and Array Processing ; Conference date: 14-09-1998 Through 16-09-1998.
  * [16] D. P. G. Foures, C. P. Caulfield, P. J. Schmid, Localization of flow structures using $\infty$ -norm optimization, J. Fluid Mech. 729 (2013) 672–701.
  * [17] J. E. Dennis, Jr., R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Vol. 16 of Classics in Applied Mathematics, SIAM, Philadelphia, PA, USA, 1996.
  * [18] I. Danaila, P. Joly, S. Kaber, M. Postel, An Introduction to Scientific Computing: Twelve Computational Projects Solved with MATLAB, Springer New York, 2008.


# Dust-ion-acoustic shock waves in magnetized plasma having super-thermal
electrons

T. Yeashna1,∗, R.K. Shikha1,∗∗, N.A. Chowdhury2,∗∗∗, A. Mannan1,†, S.
Sultana1,‡, A.A. Mamun1,§ 1 Department of Physics, Jahangirnagar University,
Savar, Dhaka-1342, Bangladesh
2 Plasma Physics Division, Atomic Energy Centre, Dhaka-1000, Bangladesh
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The propagation of dust-ion-acoustic shock waves (DIASHWs) in a three-
component magnetized plasma having inertialess super-thermal electrons,
inertial warm positive ions and negative dust grains has been investigated. A
Burgers’ equation is derived by employing the reductive perturbation method.
Under consideration of inertial warm positive ions and negative dust grains,
both positive and negative shock structures are numerically observed in the
presence of super-thermal electrons. The effects of oblique angle ($\delta$),
spectral index ($\kappa$), kinematic viscosity ($\eta$), number density and
charge state of the plasma species on the formation of the DIASHWs are
examined. It is found that the positive and negative shock wave potentials
increase with the oblique angle. It is also observed that the magnitude of the
amplitude of positive and negative shock waves is not affected by the
variation of the kinematic viscosity of plasma species but the steepness of
the positive and negative shock waves decreases with kinematic viscosity of
plasma species. The implications of our findings in space and laboratory
plasmas are briefly discussed.

###### keywords:

Dusty plasma , $\kappa$-distribution , Shock waves , Burgers’ equation.

††journal: “Plasma Physics Reports”

##  1 Introduction

Dust grains are omnipresent components in most of the astrophysical
environments, viz., interstellar clouds [1, 2, 3, 4], circumstellar clouds [5,
6, 7], interplanetary space [8, 9, 10], Earth’s magnetosphere [11], Saturn’s B
ring [12] and magnetosphere [13, 14], Jupiter’s magnetosphere [14, 15],
cometary tails [16, 17], and also in laboratory plasmas [18, 19]. In a
realistic dusty plasma medium (DPM), dust grains can collide frequently with
lighter electrons than the ions, and due to these frequent collisions, dust
grains are generally negatively charged [20, 21] massive object by collecting
electrons from the surrounding environments [22]. It was observed that dust
grains are responsible to change the dynamics of the DPM [23, 24, 25, 26, 27].

It is now well established that the Earth’s magnetosphere plasma sheet [28,
29, 30, 31], the solar wind [32, 33], Jupiter’s magnetosphere [34], and
Saturn’s magnetosphere [35] plasma contain highly energetic particles which
can exhibit high energy tails. Maxwellian distribution is only appropriate for
describing iso-thermal particles but it fails to describe the dynamics of
these highly energetic particles. The super-thermal kappa ($\kappa$)
distribution was first introduced by Vasyliunas [28] for congruent description
of these energetic particles. The parameter $\kappa$ in the $\kappa$
distribution is called the spectral index which describes the deviation of the
plasma particles from the thermally equilibrium state. The
$\kappa$-distribution reduces to the ordinary Maxwellian distribution for
large values of $\kappa$ (i.e., $\kappa\rightarrow\infty$). Shah et al. [36]
investigated the ion-acoustic (IA) shock waves (IASHWs) in the presence of
super-thermal electrons, and observed that the height of the positive shock
profile increases with an increase in the value of spectral index of
electrons. Adnan et al. [37] studied small amplitude IA solitons in a
magnetized plasma having $\kappa$-distributed electrons, and found that the
amplitude of the positive solitons increases with the electrons super-
thermality. Haider et al. [38] analysed dust-ion-acoustic (DIA) solitary waves
(DIASWs) and DIA shock waves (DIASHWs) in the presence of super-thermal
electrons, and demonstrated that positive potential associated with DIASWs
increases with the increase of super-thermality of the electrons.

Shock wave mostly arises in nature due to balance between nonlinear and
dissipative forces [39]. The collision between charged and neutral particles,
kinematic viscosity of the medium components, and Landau damping are to be
responsible for dissipation [40]. The presence of kinematic viscosity plays a
major role in generating nonlinear waves [41]. Bansal et al. [42] analysed
IASHWs in DPM having super-thermal electrons, and observed that the steepness
of shock structures decreases with the increase of kinematic viscosity of
plasma species. Borah et al. [43] studied DIASWs and DIASHWs in a three-
component DPM, and reported that the steepness of the shock structures reduces
with the increase of the kinematic viscosity of the plasma component.

A number of authors considered an external magnetic field to investigate
nonlinear electrostatic shock [41, 44] and solitary [45] waves in different
plasma medium. Hossen et al. [41] studied the electrostatic shock structures
in magnetized DPM, and found that the magnitude of the positive and negative
shock profiles increases with the oblique angle ($\delta$) which arises due to
the external magnetic field. El-Monier et al. [45] investigated nonlinear IA
solitary structures in a three-component magnetized plasma and highlighted
that the positive potential increases with the increase of oblique angle.
Kamran et al. [44] analysed DIASHWs in a three-component plasma composed of
$\kappa$-distributed electrons, mobile ions, and stationary dust grains, and
observed that the magnitude of the amplitude of negative shocks increases with
the increase of electron number density. To the best knowledge of the authors,
no attempt has been made to study the DIASHWs in a three-component magnetized
DPM by considering kinematic viscosities of both inertial warm positive ions
and negative dust grains, and inertialess super-thermal electrons. The aim of
the present investigation is, therefore, to derive Burgers’ equation and
investigate DIASHWs in a three-component magnetized DPM, and to observe the
effects of various plasma parameters on the configuration of DIASHWs.

This manuscript is organized as follows: The governing equations are described
in section 2. The Burgers’ equation is derived in section 3. Numerical
analysis is reported in section 4. A brief conclusion of our present work is
provided in section 5.

##  2 Governing Equations

We present a simplified fluid model for DIASHWs in a magnetized, three-
component DPM consisting of inertial negatively charged dust particles,
positively charged warm ions, and inertialess $\kappa$-distributed electrons.
An external magnetic field $\mathbf{B}_{0}$ has been considered in the system
directed along the $z$-axis signifying $\mathbf{B}_{0}=B_{0}\hat{z}$ where
$B_{0}$ and $\hat{z}$ represent the strength of the external magnetic field
and unit vector directed along the $z$-axis, respectively. The wave
propagation vector is assumed to produce an oblique angle ($\delta$) with the
external magnetic field. The dynamics of the magnetized DPM is governed by the
following equations:

$\displaystyle\frac{\partial\tilde{n}_{d}}{\partial\tilde{t}}+\acute{\nabla}\cdot(\tilde{n}_{d}\tilde{u}_{d})=0,$
(1)
$\displaystyle\frac{\partial\tilde{u}_{d}}{\partial\tilde{t}}+(\tilde{u}_{d}\cdot\acute{\nabla})\tilde{u}_{d}=\frac{Z_{d}e}{m_{d}}\acute{\nabla}\tilde{\psi}-\frac{Z_{d}eB_{0}}{m_{d}}(\tilde{u}_{d}\times\hat{z})+\tilde{\eta}_{d}\acute{\nabla}^{2}\tilde{u}_{d},$
(2)
$\displaystyle\frac{\partial\tilde{n}_{i}}{\partial\tilde{t}}+\acute{\nabla}\cdot(\tilde{n}_{i}\tilde{u}_{i})=0,$
(3)
$\displaystyle\frac{\partial\tilde{u}_{i}}{\partial\tilde{t}}+(\tilde{u}_{i}\cdot\acute{\nabla})\tilde{u}_{i}=-\frac{Z_{i}e}{m_{i}}\acute{\nabla}\tilde{\psi}+\frac{Z_{i}eB_{0}}{m_{i}}(\tilde{u}_{i}\times\hat{z})$
$\displaystyle\hskip
85.35826pt-\frac{1}{m_{i}\tilde{n}_{i}}\acute{\nabla}P_{i}+\tilde{\eta}_{i}\acute{\nabla}^{2}\tilde{u}_{i},$
(4) $\displaystyle\acute{\nabla}^{2}\tilde{\psi}=4\pi
e[\tilde{n}_{e}-Z_{i}\tilde{n}_{i}+Z_{d}\tilde{n}_{d}],\ $ (5)

where $\tilde{n}_{d}$ ($\tilde{n}_{i}$) is the dust (ion) number density,
$m_{d}$ ($m_{i}$) is the dust (ion) mass, $Z_{d}$ ($Z_{i}$) is the charge
state of the dust (ion), $e$ is the magnitude of electron charge,
$\tilde{u}_{d}$ ($\tilde{u}_{i}$) is the dust (ion) fluid velocity,
$\tilde{\eta}_{d}=\mu_{d}/m_{d}\tilde{n}_{d}$
($\tilde{\eta}_{i}=\mu_{i}/m_{i}\tilde{n}_{i}$) is the kinematic viscosity of
the dust (ion), $P_{i}$ is the pressure of positive ion, and $\tilde{\psi}$
represents the electrostatic wave potential. Now, we are introducing
normalized parameters, namely, $n_{d}\rightarrow\tilde{n}_{d}/n_{d0}$,
$n_{i}\rightarrow\tilde{n}_{i}/n_{i0}$, and
$n_{e}\rightarrow\tilde{n}_{e}/n_{e0}$, where $n_{d0}$, $n_{i0}$, and $n_{e0}$
are the equilibrium number densities of the negative dust grains, positive
ions, and electrons, respectively; $u_{d}\rightarrow\tilde{u}_{d}/C_{i}$,
$u_{i}\rightarrow\tilde{u}_{i}/C_{i}$ [where
$C_{i}=(Z_{i}k_{B}T_{e}/m_{i})^{1/2}$, $k_{B}$ being the Boltzmann constant,
and $T_{e}$ being temperature of the electron];
$\psi\rightarrow\tilde{\psi}e/k_{B}T_{e}$; $t=\tilde{t}/\omega_{p}^{-1}$
[where $\omega_{p}^{-1}=(m_{i}/4\pi e^{2}Z_{i}^{2}n_{i0})^{1/2}$];
$\nabla=\acute{\nabla}/\lambda_{D}$ [where $\lambda_{D}=(k_{B}T_{e}/4\pi
e^{2}Z_{i}n_{i0})^{1/2}$]. The pressure term of the positive ions can be
recognized as $P_{i}=P_{i0}(\tilde{n}_{i}/n_{i0})^{\gamma}$ with
$P_{i0}=n_{i0}k_{B}T_{i}$ being the equilibrium pressure of the positive ions,
and $T_{i}$ is the temperature of warm positive ions, and $\gamma=(N+2)/N$
(where $N$ is the degree of freedom and for three-dimensional case $N=3$, then
$\gamma=5/3$). For simplicity, we have considered
($\tilde{\eta}_{d}\approx\tilde{\eta}_{i}=\eta$), and $\eta$ is normalized by
$\omega_{p}\lambda_{D}^{2}$. The quasi-neutrality condition at equilibrium for
our plasma model can be written as $n_{e0}+Z_{d}n_{d0}\approx Z_{i}n_{i0}$.
Equations (1)$-$(5) can be expressed in the normalized form as:

$\displaystyle\frac{\partial n_{d}}{\partial t}+\nabla\cdot(n_{d}u_{d})=0,$
(6) $\displaystyle\frac{\partial u_{d}}{\partial
t}+(u_{d}\cdot\nabla)u_{d}=\alpha_{1}\nabla\psi-\alpha_{1}\Omega_{c}(u_{d}\times\hat{z})+\eta\nabla^{2}u_{d},$
(7) $\displaystyle\frac{\partial n_{i}}{\partial
t}+\nabla\cdot(n_{i}u_{i})=0,$ (8) $\displaystyle\frac{\partial
u_{i}}{\partial
t}+(u_{i}\cdot\nabla)u_{i}=-\nabla\psi+\Omega_{c}(u_{i}\times\hat{z})-\alpha_{2}\nabla
n_{i}^{\gamma-1}+\eta\nabla^{2}u_{i},$ (9)
$\displaystyle\nabla^{2}\psi=n_{d}(1-\mu_{e})-n_{i}+n_{e}\mu_{e},\ $ (10)

other plasma parameters are defined as $\alpha_{1}=Z_{d}m_{i}/Z_{i}m_{d}$,
$\alpha_{2}=\gamma T_{i}/(\gamma-1)Z_{i}T_{e}$, $\mu_{e}=n_{e0}/Z_{i}n_{i0}$,
and $\Omega_{c}=\omega_{ci}/\omega_{P}$ [where
$\omega_{ci}=Z_{i}eB_{0}/m_{i}$]. Now, the expression for the number density
of electrons following the $\kappa$-distribution can be written as [38, 37]

$\displaystyle
n_{e}=\left[1-\frac{\psi}{\kappa-3/2}\right]^{-\frac{\kappa+1}{2}},\ $ (11)

where the parameter $\kappa$ represents the non-thermal properties of the
electrons. Now, by expanding Eq. (11) up to third order in $\psi$, and
substituting in Eq. (10), we can write

$\displaystyle\nabla^{2}\psi=\mu_{e}+n_{d}(1-\mu_{e})-n_{i}+\sigma_{1}\psi$
$\displaystyle+\sigma_{2}\psi^{2}+\sigma_{3}\psi^{3}+\cdot\cdot\cdot,\ $ (12)

where

$\displaystyle\sigma_{1}=\mu_{e}[(\kappa+1)/2(\kappa-3/2)],$
$\displaystyle\sigma_{2}=\mu_{e}[(\kappa+1)(\kappa+3)/8(\kappa-3/2)^{2}],$
$\displaystyle\sigma_{3}=\mu_{e}[(\kappa+1)(\kappa+3)(\kappa+5)/48(\kappa-3/2)^{3}].\
$

We note that the terms containing $\sigma_{1}$, $\sigma_{2}$, and $\sigma_{3}$
are the contribution of $\kappa$-distributed electrons.

##  3 Derivation of the Burgers’ equation

To derive the Burgers’ equation for the DIASHWs propagating in a magnetized
plasma, we are going to employ reductive perturbation method [46]. First we
introduce the stretched co-ordinates [41, 47]

$\displaystyle\xi=\epsilon(l_{x}x+l_{y}y+l_{z}z-v_{p}t),$ (13)
$\displaystyle\tau={\epsilon}^{2}t,\ $ (14)

where $v_{p}$ express the phase speed and $\epsilon$ represent a small
parameter ($0<\epsilon<1$). The $l_{x}$, $l_{y}$, and $l_{z}$ (i.e.,
$l_{x}^{2}+l_{y}^{2}+l_{z}^{2}=1$) are the directional cosines of the wave
vector $k$ along $x$, $y$, and $z$-axes, respectively. Then, the dependent
variables can be expressed in power series of $\epsilon$ as

$\displaystyle n_{d}=1+\epsilon
n_{d}^{(1)}+\epsilon^{2}n_{d}^{(2)}+\epsilon^{3}n_{d}^{(3)}+\cdot\cdot\cdot,$
(15) $\displaystyle n_{i}=1+\epsilon
n_{i}^{(1)}+\epsilon^{2}n_{i}^{(2)}+\epsilon^{3}n_{i}^{(3)}+\cdot\cdot\cdot,$
(16) $\displaystyle
u_{dx,y}=\epsilon^{2}u_{dx,y}^{(1)}+\epsilon^{3}u_{dx,y}^{(2)}+\cdot\cdot\cdot,$
(17) $\displaystyle
u_{ix,y}=\epsilon^{2}u_{ix,y}^{(1)}+\epsilon^{3}u_{ix,y}^{(2)}+\cdot\cdot\cdot,$
(18) $\displaystyle u_{dz}=\epsilon
u_{dz}^{(1)}+\epsilon^{2}u_{dz}^{(2)}+\cdot\cdot\cdot,$ (19) $\displaystyle
u_{iz}=\epsilon u_{iz}^{(1)}+\epsilon^{2}u_{iz}^{(2)}+\cdot\cdot\cdot,$ (20)
$\displaystyle\psi=\epsilon\psi^{(1)}+\epsilon^{2}\psi^{(2)}+\cdot\cdot\cdot.\
$ (21)

Now, by substituting Eqs. (13)$-$(21) in Eqs. (6)$-$(9) and (12), we obtain a
set of first-order equations in the following form

$\displaystyle n_{d}^{(1)}=-\frac{\alpha_{1}l_{z}^{2}}{v_{p}^{2}}\psi^{(1)},$
(22) $\displaystyle u_{dz}^{(1)}=-\frac{\alpha_{1}l_{z}}{v_{p}}\psi^{(1)},$
(23) $\displaystyle
n_{i}^{(1)}=\frac{3l_{z}^{2}}{\big{(}3v_{p}^{2}-2\alpha_{2}l_{z}^{2}\big{)}}\psi^{(1)},$
(24) $\displaystyle
u_{iz}^{(1)}=\frac{3v_{p}l_{z}}{\big{(}3v_{p}^{2}-2\alpha_{2}l_{z}^{2}\big{)}}\psi^{(1)}.\
$ (25)

Now, the phase speed of DIASHWs can be written as

$\displaystyle v_{p}\equiv
v_{p+}=l_{z}\sqrt{{\frac{-a_{1}+\sqrt{a_{1}^{2}-12\sigma_{1}a_{2}}}{6\sigma_{1}}}},$
(26) $\displaystyle v_{p}\equiv
v_{p-}=l_{z}\sqrt{{\frac{-a_{1}-\sqrt{a_{1}^{2}-12\sigma_{1}a_{2}}}{6\sigma_{1}}}},\
$ (27)

where $a_{1}=3\mu_{e}\alpha_{1}-2\alpha_{2}\sigma_{1}-3-3\alpha_{1}$ and
$a_{2}=2\alpha_{1}\alpha_{2}+2\mu_{e}\alpha_{1}\alpha_{2}$. The $x$ and
$y$-components of the first-order momentum equations can be manifested as

$\displaystyle
u_{dx}^{(1)}=-\frac{l_{y}}{\Omega_{c}}\frac{\partial\psi^{(1)}}{\partial\xi},$
(28) $\displaystyle
u_{dy}^{(1)}=\frac{l_{x}}{\Omega_{c}}\frac{\partial\psi^{(1)}}{\partial\xi},$
(29) $\displaystyle
u_{ix}^{(1)}=-\frac{3l_{y}v_{p}^{2}}{\Omega_{c}(3v_{p}^{2}-2\alpha_{2}l_{z}^{2})}\frac{\partial\psi^{(1)}}{\partial\xi},$
(30) $\displaystyle
u_{iy}^{(1)}=\frac{3l_{x}v_{p}^{2}}{\Omega_{c}(3v_{p}^{2}-2\alpha_{2}l_{z}^{2})}\frac{\partial\psi^{(1)}}{\partial\xi}.\
$ (31)

Figure 1: The variation of nonlinear coefficient $A$ with $\mu_{e}$. Other
plasma parameters are $\alpha_{1}=10^{-3}$, $\alpha_{2}=0.03$,
$\delta=10^{\circ}$, $\kappa=2$, and $v_{p}\equiv v_{p+}$.

Now, by taking the next higher-order terms, the equation of continuity,
momentum equation, and Poisson’s equation can be written as

$\displaystyle\frac{\partial n_{d}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
n_{d}^{(2)}}{\partial\xi}+l_{x}\frac{\partial
u_{dx}^{(1)}}{\partial\xi}+l_{y}\frac{\partial
u_{dy}^{(1)}}{\partial\xi}+l_{z}\frac{\partial u_{dz}^{(2)}}{\partial\xi}$
$\displaystyle\hskip
51.21504pt+l_{z}\frac{\partial}{\partial\xi}\Big{(}n_{d}^{(1)}u_{dz}^{(1)}\Big{)}=0,$
(32) $\displaystyle\frac{\partial
u_{dz}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
u_{dz}^{(2)}}{\partial\xi}+l_{z}u_{dz}^{(1)}\frac{\partial
u_{dz}^{(1)}}{\partial\xi}-\alpha_{1}l_{z}\frac{\partial\psi^{(2)}}{\partial\xi}-\eta\frac{\partial^{2}u_{hz}^{(1)}}{\partial\xi^{2}}=0,$
(33) $\displaystyle\frac{\partial
n_{i}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
n_{i}^{(2)}}{\partial\xi}+l_{x}\frac{\partial
u_{ix}^{(1)}}{\partial\xi}+l_{y}\frac{\partial
u_{iy}^{(1)}}{\partial\xi}+l_{z}\frac{\partial u_{iz}^{(2)}}{\partial\xi}$
$\displaystyle\hskip
51.21504pt+l_{z}\frac{\partial}{\partial\xi}\Big{(}n_{i}^{(1)}u_{iz}^{(1)}\Big{)}=0,$
(34) $\displaystyle\frac{\partial
u_{iz}^{(1)}}{\partial\tau}-v_{p}\frac{\partial
u_{iz}^{(2)}}{\partial\xi}+l_{z}u_{iz}^{(1)}\frac{\partial
u_{iz}^{(1)}}{\partial\xi}+l_{z}\frac{\partial\psi^{(2)}}{\partial\xi}$
$\displaystyle\hskip
14.22636pt+\alpha_{2}l_{z}\Biggr{[}\frac{2}{3}\frac{\partial
n_{i}^{(2)}}{\partial\xi}-\frac{1}{9}\frac{\partial
n_{i}^{(1)^{2}}}{\partial\xi}\Biggr{]}-\eta\frac{\partial^{2}u_{iz}^{(1)}}{\partial\xi^{2}}=0,$
(35) $\displaystyle
n_{i}^{(2)}=(1-\mu_{e})n_{d}^{(2)}+\sigma_{1}\psi^{(2)}+\sigma_{2}\psi^{(1)^{2}}.\
$ (36)

Finally, the next higher-order terms of Eqs. (6)-(9) and (12), with the help
of Eqs. (22)-(36), can provide the Burgers’ equation

$\displaystyle\frac{\partial\Psi}{\partial\tau}+A\Psi\frac{\partial\Psi}{\partial\xi}=C\frac{\partial^{2}\Psi}{\partial\xi^{2}},\
$ (37)

where $\Psi=\psi^{(1)}$ is used for simplicity. The nonlinear coefficient
($A$) and dissipative coefficient ($C$) are represented, respectively, as

$\displaystyle
A=\frac{2\sigma_{2}s_{1}^{3}v_{p}^{4}+3\alpha_{1}^{2}s_{1}^{3}l_{z}^{4}(1-\mu_{e})-81v_{p}^{6}l_{z}^{4}-F_{1}}{2\alpha_{2}v_{p}l_{z}^{2}s_{1}^{3}(\mu_{e}-1)+F_{2}},$
(38) $\displaystyle C=\frac{\eta}{2},\ $ (39)

where $F_{1}=6\alpha_{2}s_{1}v_{p}^{4}l_{z}^{6}$,
$F_{2}=-18s_{1}l_{z}^{2}v_{p}^{5}$, and
$s_{1}=3v_{p}^{2}-2\alpha_{2}l_{z}^{2}$. To obtain stationary shock wave
solution of this Burgers’ equation, we consider a frame of reference that
advances with shock speed $U_{0}$. The space ($\zeta$) and time ($\tau$)
coordinates in such frame are expressed as $\zeta=\xi-U_{0}\tau^{\prime}$ and
$\tau=\tau^{\prime}$. These allow us to write the stationary shock wave
solution as [48, 49]

$\displaystyle\Psi=\Psi_{m}\Big{[}1-\tanh\bigg{(}\frac{\zeta}{\Delta}\bigg{)}\Big{]},\
$ (40)

where the amplitude $\Psi_{m}$ and the width $\Delta$ is given by

$\displaystyle\Psi_{m}=\frac{U_{0}}{A},~{}~{}\text{and}~{}~{}\Delta=\frac{2C}{U_{0}}.\
$ (41)

It can be seen from Eq. (33) that the amplitude of the shock wave becomes
infinity corresponding to the value of $A=0$ (due to $C>0$ and $U_{0}>0$), and
this refers to that our theory (specially, the RPM) is only valid for the
small amplitude waves.

Figure 2: The variation of $\Psi$ with $\zeta$ for different values of
$\delta$ under consideration of $\mu_{e}>\mu_{ec}$. Other plasma parameters
are $\alpha_{1}=10^{-3}$, $\alpha_{2}=0.03$, $\eta=0.3$, $\kappa=2$,
$\mu_{e}=0.8$, $U_{0}=0.2$, and $v_{p}\equiv v_{p+}$. Figure 3: The variation
of $\Psi$ with $\zeta$ for different values of $\delta$ under consideration of
$\mu_{e}<\mu_{ec}$. Other plasma parameters are $\alpha_{1}=10^{-3}$,
$\alpha_{2}=0.03$, $\eta=0.3$, $\kappa=2$, $\mu_{e}=0.4$, $U_{0}=0.2$, and
$v_{p}\equiv v_{p+}$. Figure 4: The variation of $\Psi$ with $\zeta$ for
different values of $\kappa$ under consideration of $\mu_{e}>\mu_{ec}$. Other
plasma parameters are $\alpha_{1}=10^{-3}$, $\alpha_{2}=0.03$,
$\delta=10^{\circ}$, $\eta=0.3$, $\mu_{e}=0.8$, $U_{0}=0.2$, and $v_{p}\equiv
v_{p+}$.

##  4 Numerical analysis

The balance between nonlinearity and dissipation leads to generate DIASHWs in
a three-component magnetized DPM. We have numerically analyzed the variation
of $A$ with $\mu_{e}$ in Fig. 1, and it is obvious from this figure that (a)
$A$ can be negative, zero, and positive depending on the values of $\mu_{e}$;
(b) the value of $\mu_{e}$ for which $A$ becomes zero is known as critical
value of $\mu_{e}$ (i.e., $\mu_{ec}$), and the $\mu_{ec}$ for our present
analysis is almost $0.55$; and (c) the parametric regimes for the formation of
positive (i.e., $\psi>0$) and negative (i.e., $\psi<0$) potential shock
structures can be found corresponding to $A>0$ and $A<0$.

It is clear from Figs. 2 and 3 that (a) with the increase of the oblique angle
($\delta$), the magnitude of the amplitude of the positive and negative shock
structures increases, and this result agrees with the result of Hossen et al.
[41]; (b) the magnitude of the negative potential is always greater than the
positive potential for same plasma parameters. So, the oblique angle enhances
the height of the potential shock structures.

Figure 5: The variation of $\Psi$ with $\zeta$ for different values of $\eta$
under consideration of $\mu_{e}>\mu_{ec}$. Other plasma parameters are
$\alpha_{1}=10^{-3}$, $\alpha_{2}=0.03$, $\delta=10^{\circ}$, $\kappa=2$,
$\mu_{e}=0.8$, $U_{0}=0.2$, and $v_{p}\equiv v_{p+}$. Figure 6: The variation
of $\Psi$ with $\zeta$ for different values of $\mu_{e}$ under consideration
of $\mu_{e}<\mu_{ec}$. Other plasma parameters are $\alpha_{1}=10^{-3}$,
$\alpha_{2}=0.03$, $\delta=10^{\circ}$, $\eta=0.3$, $\kappa=2$, $U_{0}=0.2$,
and $v_{p}\equiv v_{p+}$.

Figure 4 displays the variation of positive shock profile with spectral index
$\kappa$, and it is evident from the figure that (a) shock height is amplified
for higher values of $\kappa$, and this result agrees with the previous work
of Ref. [36]; (b) positive shock height changes abruptly with the increase of
the value of $\kappa$.

It is evident from Fig. 5 that there are some specific correlations between
the dust-ion kinematic viscosity on the positive (under the consideration
$\mu_{e}>\mu_{ec}$) shock profiles. It is really interesting that the
magnitude of the amplitude of positive shock profiles is not affected by the
variation of the dust-ion kinematic viscosity but the steepness of the shock
profile decreases with the increase of dust-ion kinematic viscosity, and this
result is in good agreement with the result of Bansal et al. [42] and Borah et
al. [43].

Figure 6 illustrates the effects of the equilibrium density ratio of electrons
to ions (via $\mu_{e}$) on plasma shock structures. The numerical analysis
exhibits amplification in the magnitude of the amplitude of negative shock
profile for higher values of $\mu_{e}$ which fully agrees with the result of
Kamran et al. [44]. However, the variation of the value of $\mu_{e}$
drastically changes the height of the shock profile.

##  5 Conclusion

We have studied DIASHWs in a three-component magnetized DPM by considering
kinematic viscosities of both negative dust and positive ion species, and
inertialess super-thermal electrons. The reductive perturbation method is used
to derive the Burgers’ equation. The results that have been found from our
investigation can be summarized as follows:

  * 1.

The parametric regimes for the formation of positive (i.e., $\psi>0$) and
negative (i.e., $\psi<0$) potential shock structures can be found
corresponding to $A>0$ and $A<0$.

  * 2.

The magnitude of the amplitude of positive and negative shock structures
increases with the oblique angle ($\delta$).

  * 3.

The magnitude of the amplitude of positive and negative shock profiles is not
affected by the variation of the dust-ion kinematic viscosity ($\eta$) but the
steepness of the shock profile decreases with dust-ion kinematic viscosity
($\eta$).

It should be noted here that the gravitational effect is of great importance
for DPM but it is beyond the scope of our present work. In future for better
understanding, someone can investigate the nonlinear propagation in a three-
component DPM by considering the gravitational effect. The results of our
present investigation will be useful in understanding the nonlinear phenomena
both in astrophysical environments such as interstellar clouds [1, 2, 3],
circumstellar clouds [5, 6, 7], interplanetary space [8, 9, 10], Earth’s
magnetosphere[11], Saturn’s B ring [12] and magnetosphere [13, 14], Jupiter’s
magnetosphere [14, 15], cometary tails [16, 17] and in laboratory plasmas [18,
19].

## References

  * [1] L. Spitzer, Physical Processes in the Interstellar Medium, (Wiley, New York, 1978).
  * [2] B.D. Savage and J.S. Mathis, Ann. Rev. Astron. Astrophys. 17, 73 (1979).
  * [3] E. Herbst, Annu. Rev. Phys. Chem. 46, 27 (1995).
  * [4] M.H. Rahman,et al., Phys. Plasmas 25, 102118 (2018); R.K. Shikha, et al., Eur. Phys. J. D 73, 177 (2019); S.K. Paul, et al., Pramana J. Phys. 94, 58 (2020); S. Jahan, et al., Commun. Theor. Phys. 71, 327 (2019); N.A. Chowdhury, et al., Phys. plasmas 24, 113701 (2017); M.H. Rahman, et al., Chin. J. Phys. 56, 2061 (2018); D.M.S. Zaman, et al., High Temp. 58, 789 (2020); S. Jahan, et al., Plasma Phys. Rep. 46, 90 (2020).
  * [5] V. Ossenkopf, et al., Astron. Astrophys. 261, 567 (1992).
  * [6] S.M. Andrews and J.P. Williams, Astrophys. J. 631, 1134 (2005).
  * [7] C. Dijkstra, et al., Astrophy. J. 633, L133 (2005).
  * [8] I.D.R. Mackinnon and F.J.M. Rietmeijer, Rev. Geophys. 25, 7 (1987).
  * [9] S.A. Sandford and J.P. Bradley, Icarus 82, 146 (1989).
  * [10] J.C. Liou and H.A. Zook, Icarus 124, 429 (1996).
  * [11] M. Horanyi, et al., Astrophys. Space Sci. 144, 215 (1988).
  * [12] T.G. Northrop and J.R. Hill, J. Geophys. Res. 88, A8 (1983).
  * [13] D.H. Humes, J. Geophys. Res. 85, 5841 (1980).
  * [14] T.G. Northrop and J.R. Hill, J. Geophys. Res. 88, A1 (1983).
  * [15] D.H. Humes, et al., J. Geophys. Res. 79, 25 (1974).
  * [16] E.P. Mazets, et al., Astron. Astrophys. 187, 699 (1987).
  * [17] D.A. Gurnett, et al., Geophys. Res. Lett. 13, 3 (1986).
  * [18] G.S. Selwyn, et al., J. Vac. Sci. Technol. A, 7, 2758 (1989).
  * [19] L. Boufendi and A. Bouchoule, Plas. Sour. Sci. Technol. 11, A211 (2002).
  * [20] T.G. Northrop and J.R. Hill, J. Geophys. Res. 87, A8 (1982).
  * [21] R.L. Merlino, et al., Phys. Plasmas 5, 1607 (1998).
  * [22] P.K. Shukla, Phys. Script. 45, 504 (1992).
  * [23] S. Ghosh, et al., Phys. Lett. A 274, 162 (2000).
  * [24] Y. Nakamura, Phys. Plas. 9, 440 (2002).
  * [25] A.A. Mamun and K.S. Ashrafi, Phys. Rev. E 82, 026405 (2010).
  * [26] D. Samsonov, et al., Phys. Rev. Lett. 92, 25 (2004).
  * [27] M. Bacha, et al., Phys. Rev. E 85, 056413 (2012).
  * [28] V.M. Vasyliunas, J. Geophys. Res. 73, 9 (1968).
  * [29] A.T.Y. Lui and S.M. Krimigis, Geophys. Res. Lett. 10, 13 (1983).
  * [30] D.J. Williams, et al., Geophys. Res. Lett. 15, 4 (1988).
  * [31] S.P. Christon, et al., J. Geophys. Res. 93, 2562 (1988).
  * [32] B.A. Shrauner and W. C. Feldman, J. Plasma Phys. 17, 123 (1977).
  * [33] J.T. Gosling, et al., J. Geophys. Res. 86, 547 (1981).
  * [34] N. Divine and H.B. Garrett, J. Geophys. Res. 88, 6889 (1983).
  * [35] T.P. Armstrong, et al., J. Geophys. Res. 88, 8893 (1983).
  * [36] A. Shah, et al., Phys. Plasmas 19, 032302 (2012).
  * [37] M. Adnan, et al., Advan. Spa. Res. 53, 845 (2014).
  * [38] M.M. Haider and A. Nahar, Z. Naturforsch. 72, 627 (2017).
  * [39] Y.B. Zel’dovich and Y.P. Raizer, Physics of Shock Waves and High- Temperature Hydrodynamic Phenomena (Academic, New York, 1967).
  * [40] B. Sahu, et al., Phys. Plasmas. 21, 103701 (2014).
  * [41] M.M. Hossen, et al., High Energy Density Phys. 24, 9 (2017).
  * [42] S. Bansal, et al., Eur. Phys. J. D 74, 236 (2020).
  * [43] P. Borah, et al., Phys. Plasmas. 23, 103706 (2016).
  * [44] M. Kamran, et al., Resul. Phys. 21, 103808 (2021) .
  * [45] S.Y. El-Monier and A. Atteya, AIP Adv. 9, 045306 (2019).
  * [46] M. Hassan, et al., Commun. Theor. Phys. 71, 1017 (2019); N.A. Chowdhury, et al., Plasma Phys. Rep. 45, 459 (2019); N.A. Chowdhury, et al., Contrib. Plasma Phys. 58, 870 (2018); S. Banik, et al., Eur. Phys. J. D 75, 43 (2021); N.A. Chowdhury, et al., Chaos 27, 093105 (2017); N.A. Chowdhury, et al., Vacuum 147, 31 (2018); N. Ahmed, et al., Chaos 28, 123107 (2018).
  * [47] H. Washimi and T. Taniuti, Phys. Rev. Lett. 17, 996 (1966).
  * [48] V.I. Karpman, Nonlinear Waves in Dispersive Media, (Pergamon Press, Oxford, 1975).
  * [49] A. Hasegawa, Plasma Instabilities and Nonlinear Effects, (Springer-Verlag, Berlin, 1975).