# Superclustering by finding statistically significant separable groups of
optimal gaussian clusters

Oleg I.Berngardt<EMAIL_ADDRESS>
###### Abstract

The paper presents the algorithm for clustering a dataset by grouping the
optimal, from the point of view of the BIC criterion, number of Gaussian
clusters into the optimal, from the point of view of their statistical
separability, superclusters.

The algorithm consists of three stages: representation of the dataset as a
mixture of Gaussian distributions - clusters, which number is determined based
on the minimum of the BIC criterion; using the Mahalanobis distance, to
estimate the distances between the clusters and cluster sizes; combining the
resulting clusters into superclusters using the DBSCAN method by finding its
hyperparameter (maximum distance) providing maximum value of introduced matrix
quality criterion at maximum number of superclusters. The matrix quality
criterion corresponds to the proportion of statistically significant separated
superclusters among all found superclusters.

The algorithm has only one hyperparameter - statistical significance level,
and automatically detects optimal number and shape of superclusters based of
statistical hypothesis testing approach. The algorithm demonstrates a good
results on test datasets in noise and noiseless situations. An essential
advantage of the algorithm is its ability to predict correct supercluster for
new data based on already trained clusterer and perform soft (fuzzy)
clustering. The disadvantages of the algorithm are: its low speed and
stochastic nature of the final clustering. It requires a sufficiently large
dataset for clustering, which is typical for many statistical methods.

###### keywords:

Gaussian mixture , Mahalanobis distance , statistically significant separation
, fuzzy clustering

††journal: Pattern Recognition

[istp]organization=Institute of Solar-Terrestrial physics SB
RAS,addressline=126a, Lermontov Str., city=Irkutsk, postcode=664033,
country=Russia

##  1 Introduction

Data clustering algorithms are usually divided into two big categories:
statistical methods and metric methods. Statistical ones describe the data by
a combination of statistical distribution functions, the basis of metric
methods is the pairwise distance between points. The popular metric method is
Density-based spatial clustering of applications with noise (DBSCAN),
suggested in (Ester et al., 1996), and discussed in (Schubert et al., 2017).
Widely used example of statistical methods is the mixture of normal (gaussian)
distributions (Gaussian Mixture), discussed in (Reynolds, 2015). Often, metric
methods are less resistant to outliers and noise than statistical methods.

Most of the methods have one or more hyperparameters \- for DBSCAN this is the
maximum distance between a pair of objects, for K-means(Arthur and
Vassilvitskii, 2007), for Spectral(Ng et al., 2001) and Gaussian Mixture, this
is the number of classes. One of the main tasks for a researcher is to choose
a clustering method and find the hyperparameters values that provide the best
solution of his clustering problem.

The choice of a clusterer is related to the data model - the clusterer that is
optimal for the researcher should correspond to the data: on the one hand, it
should have the necessary freedom to correctly describe the dataset for any
expert separation of the dataset into a clusters - in terms of the Kleinberg
theorem (Kleinberg, 2002) \- to provide richness. On the other hand it should
correctly label the original unlabeled dataset so that the number and shape of
clusters coincide with what an expert expect from the dataset.

Let us explain the main idea of the algorithm, suggested in the paper. A
dataset consisting of two expert-labeled classes can be separated into an
arbitrary number of clusters, and by many ways.

Let us consider the Rand index (Rand, 1971) at dataset $S$ in classification
(not clusterization) problem:

$RI=\frac{TP+TN}{C_{2}^{N}}$ (1)

where TP (True Positive) is the number of pairs of elements that are in the
same clusters in set $S_{l}$ (dataset $S$, labeled by an expert) and in the
same cluster in set $S_{c}$ (dataset $S$, clustered by an algorithm); TN (True
Negative) is the number of pairs of elements that are in different clusters in
set $S_{l}$ and in different clusters in set $S_{c}$; $C_{2}^{N}$ \- total
number of pairs in dataset $S$.

In this case, pairwise true positive (PWTP) metric can be defined as:

$PWTP=\frac{TP}{C_{2}^{N}}$ (2)

and pairwise true negative (PWTN) metric can be defined as:

$PWTN=\frac{TN}{C_{2}^{N}}$ (3)

Rand Index is sum of PWTP and PWTN. Obviously, the easiest way to achieve the
maximum of PWTN metric, is to put each point into a separate cluster, but PWTP
metric in this variant will become zero. The easiest way to achieve the
maximum of PWTP metric, is to put all the points into the same single cluster,
but PWTN metric in this variant will become zero. In semi-supervised learning
(Lu and Leen, 2004) PWTN and PWTP metrics can be interpreted as quality of
cannot-links and must-links solutions correspondingly.

Let us refer the optimal PWTN algorithm as the algorithm providing maximum
PWTN, and refer the optimal PWTP algorithm as the algorithm providing maximum
PWTP. We suggest that by corresponding grouping clusters produced by optimal
PWTN algorithm we can obtain optimal PWTP algorithm. Adequate classification
can be defined as a clustering that is both optimally PWTP and optimally PWTN
with respect to the teacher’s labeling. It is obvious that this corresponds to
the absolute maximum of RI because in this case the errors of the first and
second kind (FN, FP) are equal to zero. The absolute maximum of both PWTP and
PWTN metrics is reached simultaneously when the number of clusters in the
labeled dataset $S_{l}$ and clustered dataset $S_{c}$ are the same, and their
shapes are also the same. Therefore, the adequate classification problem in
the presence of a labeled dataset $S_{l}$ can be considered as the problem of
achieving RI its absolute maximum value of 1. A good approximation to adequate
classification can be found by supervised learning.

However, in most cases, it seems to be impossible to achieve adequate
clustering by unsupervised learning. We do not know the set $S_{l}$, and
therefore RI will be most likely less than 1 due to the imperfection of our
knowledge about $S_{l}$. Therefore, each clustering model has its own
limitations \- our initial assumptions about the unknown set $S_{l}$: K-Means
searches for centroids, DBSCAN for distance separated sets, Gaussian Mixture
for sets which distributions are close to normal ones.

Approximation by Gaussian distributions is widely used in clustering because
potentially it can provide a high PWTN metric - an infinite sum of Gaussian
distributions can fit (and separate) any arbitrary discrete set of points.

This paper demonstrates an algorithm that uses grouping of gaussian clusters
into superclusters to create a clustering close to adequate one, based on our
own assumptions about the set $S_{l}$, that will be explained later.

##  2 GMSDB algorithm

###  2.1 The idea

We will solve the problem sequentially: at the first stage, we will
approximate our data most accurately by the optimal number of clusters from
the point of view of likelihood (we will construct a clustering with a
potentially high PWTN metric). At the second stage we will group these
clusters into the optimal number of superclusters (we will increase PWTP
metric) from the point of view of the distance between these superclusters.
Optimality in terms of likelihood at the first stage should allow us to best
separate the data clusters from a statistical point of view, and at the second
stage we group too similar (or close) clusters. The problem conceptually is
close to the well-known hierarchical method - agglomerative clustering(Sibson,
1973), but uses only two stages instead of using their sequence.

Important requirements for our algorithm are: robustness to a noise,
interpretability of each of its stages in terms of a statistical approach, and
minimizing the number of hyperparameters.

The first stage can be simply solved when we are not limited by the maximum
number of clusters: we can choose the number of clusters equal to the number
of objects and place each object in its own cluster. If we want to limit the
number of clusters we need to use some kind of clustering as an initial stage.
We choose an approximation by mixture of normal distributions \- Gaussian
Mixture (GM). The advantage of this method is that it is statistical one.
Another advantage - it has a widely used criterion for finding the optimal
number of objects - the Bayesian Information Criterion - BIC, (Schwarz, 1978).
Therefore, we expect that most likely the optimal number of initial clusters
will be less than the number of points (objects) in the dataset. The
disadvantage of this method is its stochastic nature associated with the
iterative Expectation-Maximization (EM) algorithm used for search of GM
parameters, and dependence of EM results on its stochastic initial state
(Reynolds, 2015).

The second stage is close to agglomerative clustering. The main principle of
agglomerative clustering is sequential grouping of the closest clusters. In
agglomerative approach, one should group them until all the clusters are
grouped into a single supercluster, or until necessary number of superclusters
reached, or until optimum of a certain criterion reached.

Therefore two main problems arise: by what principle to combine clusters, and
by what criteria to detect if the clustering becomes the optimal one. There
are many methods for solving the first problem. The simplest method is
sequential one - greedy joining. In this case, we find a pair of closest
clusters, and group them into one (Sibson, 1973). Another common method is
density based distribution (DBSCAN). However, DBSCAN has problems in choosing
value of a hyperparameter - the maximum distance between points in a cluster
($\varepsilon$). So one need to vary (or to use some search method) to find
optimal $\varepsilon$, which slows down the algorithm. We choose DBSCAN method
due to it is easy to implement.

The second problem is the choice of metric \- the distance between clusters.
Since we use GM, one of the native metrics is the Mahalanobis distance
(Mahalanobis, 1936), which scales the space depending on the covariance matrix
of the found Gaussian distribution. The distance between two clusters can also
be calculated using this metric. Our algorithm stops when the distance between
all found superclusters exceeds a given threshold level.

The name of our algorithm (GMSDB) is formed from the three elements it based
on: Gaussian Mixture (GM), statistical approach (S) and DBSCAN (DB). Let us
describe the GMSDB algorithm in details.

###  2.2 Stage 1: Data approximation by the optimal GM

At this stage, we iterate over the number of clusters $N$ until the minimum of
the BIC criterion(Schwarz, 1978) is reached (or until the number of clusters
reaches a high enough value that we have specified). In this stage, we assume
that the set of points $\overrightarrow{x}$ satisfies the distribution with
the probability density $P_{GM}(\overrightarrow{x})$:

$P_{GM}(\overrightarrow{x})=\sum_{n=1}^{N}A_{i}P(\overrightarrow{x}|\theta_{i})$
(4)

$\sum_{n=1}^{N}A_{i}=1$ (5)

Here $P(\overrightarrow{x}|\theta_{i})$ is the normal distribution described
by the set of its unknown parameters $\theta_{i}$ (mean and covariance
matrix), and unknown weights $A_{i}$ of the corresponding distribution in the
sum. N is the number of distributions in the mixture. The optimal
distributions $P_{GM}(\overrightarrow{x})$ and
$P(\overrightarrow{x}|\theta_{i})$, as well as $N=N_{BIC},A_{i},\theta_{i}$,
at which the minimum of the BIC criterion is reached, will be the result of
this stage.

Fig.1 shows examples of the BIC dependence on the number of clusters $N$ and
the resulting splitting of the original dataset into the optimal number
$N_{BIC}$ Gaussian clusters $P(\overrightarrow{x}|\theta_{i})$, as well as
decision regions, defining the boundaries between these clusters. Obviously,
only the first example (’blobs’, shown Fig.1A) is an adequate enough model;
the other cases (Fig.1B-D) are just some clusterings with high enough PWTN.

Figure 1: Clustering points into BIC-optimal mixture of Gaussian
distributions. On the left is a separation by the number of clusters $N_{BIC}$
corresponding to the BIC minimum. The background corresponds to decision
areas, all points of which the algorithm will refer to given cluster. Colors
correspond to different clusters. On the right is the dependence of BIC on the
number of clusters. A) ’medium blobs’; B) ’two horseshoes’; C) ’three nested
rings with noise’; D) ’two nested rings with noise’

###  2.3 Stage 2: Calculation of cluster sizes and intercluster distances

We require the algorithm to be robust to a noise and interpretable from
statistical point of view. So metrics should be carried out not in Euclidean
space, but in probability space. We used a mixture of Gaussian distributions
so the distance between clusters and cluster sizes are calculated using
Mahalanobis distances (Mahalanobis, 1936), native for gaussian statistical
distribution:

$d(\overrightarrow{x},\overrightarrow{y})=\sqrt{\overrightarrow{x}S_{y}^{-1}\overrightarrow{y}^{T}}$
(6)

where $S_{y}^{-1}$ is the matrix inverse to the covariance matrix of the
cluster to which point $\overrightarrow{y}$ belongs.

The distance is widely used in different clustering and classification tasks
(Shiming et al., 2008). Qualitatively, the Mahalanobis distance is the
distance from a point $\overrightarrow{x}$ to an ellipsoid normalized to the
ellipsoid width. Therefore it could be qualitatively interpreted as
t-statistic value(Student, 1908). Thus, the Mahalanobis distance is related to
the likelihood of the hypothesis (t-test) that the point $\overrightarrow{x}$
belongs to a given normal distribution, to which $\overrightarrow{y}$ belongs.
When it reaches a certain threshold value (for example, T=1.96), we can talk
about the rejection of this hypothesis with the corresponding statistical
significance level ($\alpha=0.05$).

Since in the future we will evaluate the quality of clustering from an
approach close to the Silhouette criterion (Rousseeuw, 1987), based on the
concepts of internal (cluster) size and external (intercluster) distance, we
will need to evaluate not only the distance between two clusters, but also
their own sizes. The Mahalanobis distance is usually defined for the case when
both points are from the same distribution. We will use points from different
distributions (clusters). Accordingly, their covariance matrices ($S_{y}$ and
$S_{x}$) in the general case will be different, and the distance calculated
(6) will not be a distance, since it is not symmetrical:

$\sqrt{\overrightarrow{x}S_{y}^{-1}\overrightarrow{y}^{T}}\neq\sqrt{\overrightarrow{y}S_{x}^{-1}\overrightarrow{x}^{T}}$
(7)

In this case, we have 4 statistical distributions of these quasi-metrics:
$p_{ii}(d),p_{jj}(d),p_{ji}(d),p_{ij}(d)$, where the indices i,j correspond to
the points in i-th and j-th clusters, obtained in Stage 1 GM/BIC
clusterization (4). The distributions $p_{ii}(d),p_{jj}(d)$ are the cluster
sizes, and $p_{ij}(d),p_{ji}(d)$ are the intercluster distances taking into
account the j-th and i-th cluster shapes, respectively.

The intercluster distance $R_{ij}$ must be symmetrical ($R_{ij}=R_{ji}$), so
we should form a symmetric functions from $p_{ij}(d),p_{ji}(d)$. As the
average size of the clusters $L_{i},L_{j}$, we choose the median values of the
distance between the points of the corresponding cluster, and as the
intercluster distance - the $\alpha$-th percentile:

$\begin{array}[]{c}L_{i}=median(p_{ii}(x))\\\
R_{ij}=Max(Percentile_{5\%}(p_{ij}),Percentile_{5\%}(p_{ji}))\\\
R_{ii}=0\end{array}$ (8)

The 5th percentile (5%) was chosen for statistical reasons and corresponds to
the lower limit, below which, with a standard statistical confidence level of
95%, the distance between points in the same cluster does not fall down. In
addition, percentiles are more robust to outliers than standard deviations and
means. The maximization between two values is used to symmetrize the distance.
It corresponds to the transition to the coordinate system of the cluster in
which the intercluster distance is higher, and the two clusters can be more
confidently separated. The use of the median as a characteristic cluster size
can also be justified by its resistance to outliers.

Fig.2A-D shows examples of clusters and their corresponding Mahalanobis
distance distributions (6). It can be seen from the figure that the distance
distributions $p_{ij}$ and $p_{ji}$ generally do not coincide, which explains
distance symmetrization (8). The scheme for determining distances is
illustrated in Fig.2E.

Figure 2: A-D) Distributions of points within clusters (different colors
correspond to different clusters) and distributions of Mahalanobis distances
corresponding to them (orange and green correspond to distributions of
intercluster distances $p_{ij}(d),p_{ji}(d)$, red and blue correspond to
cluster size distributions $p_{ii}(d),p_{jj}(d)$). E) An illustration of
determining the cluster size $L_{i}$ and the intercluster distance $R_{ij}$
from the distributions of Mahalanobis distances for two clusters.

###  2.4 Stage 3: Grouping clusters

DBSCAN was used as an algorithm for grouping the clusters into superclusters,
matrix of intercluster distances $R_{ij}$ 8 is used as inter-object distances.
To simplify the enumeration of the iterations over the maximal element
distance in the supercluster $\varepsilon$, the DBSCAN $\varepsilon$
hyperparameter is used as iteration identifier. To do this all unique values
of $\varepsilon$ were ordered in increasing order, and intermediate values
$\overline{\varepsilon}_{i}=(\varepsilon_{i-1}+\varepsilon_{i})/2;\varepsilon_{0}=0$
between the neighbours of this sequence were calculated. The resulting
sequence of $\overline{\varepsilon}_{i}$ is used as iteration indices.
Obviously, at the smallest distance $\overline{\varepsilon}_{1}$ DBSCAN
produces the maximum possible value of clusters $N_{BIC}$ (determined from the
optimality of the BIC criterion at the first stage of the algorithm), and at
the largest $\overline{\varepsilon}_{j}$ DBSCAN produces only one cluster. In
order to speed up algorithm, if iterating over 10 consecutive values of
$\overline{\varepsilon}_{i}$ DBSCAN produces the number of clusters equal to
1, the iteration process over $\overline{\varepsilon}_{i}$ is stopped.

Further calculations of stop criteria require a matrix of distances between
new superclusters $D_{ij}(\overline{\varepsilon})$. The distance
$D_{SC_{i},SC_{j}}(\overline{\varepsilon})$ between two superclusters
$SC_{i},SC_{j}$ is taken to be the minimum distance between clusters
$C_{n},C_{m}$ belonging to these two superclusters :

$D_{SC_{i},SC_{j}}(\overline{\varepsilon})=Min_{m,n:C_{m}\subseteq
SC_{i},C_{n}\subseteq SC_{j}}(R_{mn})$

and as the size $D_{SC_{i},SC_{i}}(\overline{\varepsilon})$ of the
supercluster $SC_{i}$, the minimum size of the cluster belonging to this
supercluster is chosen:

$D_{SC_{i},SC_{i}}(\overline{\varepsilon})=Min_{m:C_{m}\subseteq
SC_{i}}(L_{m})$

In terms of speed, this algorithm is intermediate between a greedy algorithm
(for example, agglomerative) and an exhaustive search. At the same time, at
each iteration of $\overline{\varepsilon}$, it accurately estimates the number
of superclusters produced by this solution. It should be noted that despite
the fact that the DBSCAN method is a metric method, the used distance matrix
$R_{ij}$ (defined in the previous step) is essentially a set of t-statistic
values. Therefore, when performing this stage, we cluster the differences in
the statistical characteristics of clusters, rather than their Cartesian
characteristics.

###  2.5 Stage 4: Stop criteria

An important task is the choice of the optimal superclusters configuration and
their number, which is reduced to the choice of the optimality criterion. It
is obvious that new use of the BIC criterion will not bring any sufficient
results: according to this criterion, we have already chosen the maximum
number of clusters $N_{BIC}$. The other widely used criterion is the
Silhouette criterion. However, it has a problem - it is related with the
concept of subtracting distances (and hence connected with the Euclidean
metric), which makes its statistical interpretation difficult, when distances
are t-statistics. Therefore, we use a more clear criterion from statistical
point of view - the matrix quality criterion (MC).

By analogy Silhouette criterion, to calculate MC the distance matrix
$D_{ij}(\overline{\varepsilon})$ is normalized so its diagonals become 1. This
makes a transformation in which all superclusters are scaled to the same size
$\overline{D}_{ii}(\overline{\varepsilon})=1$:

$\overline{D}_{ij}(\overline{\varepsilon})=D_{ij}(\overline{\varepsilon})/D_{ii}(\overline{\varepsilon})$

Statistically, this can be interpreted as the resulting distances between
superclusters are measured in fractions of their covariances, and can be
interpreted as t-statistics for the superclusters. It is convenient to
interpret them from the point of view of testing hypotheses. So the native
condition for separability of superclusters arises. If for some i-th
supercluster there is exists a j-th supercluster $j\neq i$, the distance
$\overline{D}_{ij}(\overline{\varepsilon})$ to which is below the given
threshold level $\delta D$:

$\overline{D}_{ij}(\overline{\varepsilon}_{i})<\delta D$

the superclusters i-th and j-th coincide, and the number of superclusters can
be reduced by making next iteration $\overline{\varepsilon}_{i+1}$.

Since we initially work with normal distributions, and the distances
$\overline{D}_{ij}(\overline{\varepsilon})$ are interpreted as t-statistics,
we choose threshold value as the value of the t-statistic for the significance
level $\alpha=0.05$ , which is about 1.96. Thus, if two superclusters are
located at a distance $\overline{D}_{ij}(\overline{\varepsilon})$ less than
$\delta D=1.96$, they can be considered similar (intersecting) with a
significance level of 0.05.

Therefore, as a matrix quality criterion (MC), we use the proportion of
superclusters that do not have similar superclusters:

$MC(\overline{\varepsilon})=\frac{|\\{1|Min_{i\neq
j}(\overline{D}_{ij}(\overline{\varepsilon}))>\delta
D\\}|}{|\overline{D}_{ii}(\overline{\varepsilon})|}$ (9)

where $|X|$ is the number of elements.

For each supercluster, we determine whether the supercluster closest to it is
separable from it or not. We divide the resulting number of separable
superclusters by the total number of superclusters.

In terms of $p_{value}$ (due to $\overline{D}_{ij}$ is interpreted as
t-statistic), the MC criterion will be:

$MC(\overline{\varepsilon})=\frac{|\\{1|Max_{i\neq
j}(p_{value,ij}(\overline{\varepsilon}))<\alpha\\}|}{|p_{value,ii}(\overline{\varepsilon})|}$
(10)

and has clear statistical interpretation: maximal value 1 of this criterion
demonstrates that every supercluster differs from others with significance
level at least $\alpha$.

Obviously, $MC(\overline{\varepsilon})$ lies within $[0..1]$, and reaches its
maximum value of 1 when we find the clustering with statistically significant
separable superclusters with significance level $\alpha$.

Fig.3A-B shows dependence of MC criteria on iteration number and number of
superclusters. Fig.3C-J shows the shapes of superclusters tested during
training. The Fig.3A-B demonstrates the need to stop at the
$\overline{\varepsilon}$ stage, at which $MC(\overline{\varepsilon})$ becomes
1 : after this the number of superclusters only decreases, and the value of
the $MC(\overline{\varepsilon})$ criterion does not change. From the example
shown in the figure, it can be seen that for the first time MC reaches the
value of 1 at the 39th iteration, which leads to the splitting of the dataset
into 4 superclusters we expected \- three nested rings and a noise. Further
iterations do not improve the clustering, and leave MC unchanged.

Figure 3: Dependence of the MC (A), number of superclusters (B) and
superclustering during training (C-J) as function of iteration number using
’three nested rings with noise’ dataset as an example.

##  3 Testing

At first we test the algorithm in noise-free situations. The algorithm was
tested on artificial datasets, the code of which is given in (Berngardt,
2023). Some of these datasets are identical to those used for sklearn library
(Pedregosa et al., 2011) clustering tests (https://scikit-
learn.org/stable/modules/clustering.html). All the datasets were generated by
various transformations from a normal and uniform distributions, so datasets
have varying cluster shapes and point densities within clusters.

Fig.4 shows the results of such testing. The ’three grains’ (Fig.4A), large
(Fig.4B) and small (Fig.4C) ’blobs’, which are useful for testing clustering
with statistical models and models based on centroids and distances between
them; ’three horseshoes’ (Fig.4D) and ’two horseshoes’ (Fig.4E), which are
useful for testing clustering with distance-based models, and ’three nested
rings’ (Fig.4F) and ’two nested rings’ (Fig.4G).

The results of the GMSDB algorithm and the Agglomerative algorithm (Sibson,
1973) was compared in the following way. The maximum possible number of
clusters is 50, which exceeds $N_{BIC}$ for any dataset. For agglomerative
clustering, the distance is Euclidean, and the distance between clusters is
calculated by their nearest elements.

The figure shows that in the case of well-separated ’grains’ (Fig.4A) and
’small blobs’ (Fig.4C) the algorithms work very similarly and the results are
the same as expected \- 3 and 4 main classes correspondingly. It can be seen
that the agglomerative algorithm leaves artifacts on the cluster boundaries.
The GMSDB algorithm produces no artifacts, and accurately guesses the number
of superclusters. A similar situation occurs in the case of ’three horseshoes’
and ’two horseshoes’ (Fig.4D,E) - agglomerative clustering guesses the main
body of the cluster quite well, but the edges with clustering artifacts become
even larger. The GMSDB algorithm continues to correctly guess the number and
shape of superclusters. ’Big blobs’ (Fig.4B), where the distance between
clusters is of the order of their size, is problematic for both algorithms,
but GMSDB guesses a little closer than the agglomerative one. The biggest
difference is in ’three nested rings’ and ’two nested rings’ (Fig.4F,G) - the
agglomerative one cannot adequately separate them - as the radius increases
(and the density of points decreases), the ring is divided into an increasing
number of clusters. Unlike the agglomerative algorithm, GMSDB adequately
clusterizes these datasets.

Figure 4: Comparison of clustering of different identical datasets by GMSDB
and agglomerative algorithms (minimal interobject distance as intercluster
distance, metric - euclidean). Noise-free situation. The maximum allowed
number of clusters is 50. From left to right - the original dataset,
clustering by the GMSDB algorithm and the found optimal number of
superclusters, decision regions of the trained GMSDB algorithm, clustering by
the agglomerative algorithm. Bottom left corner in second column images - the
number of superclusters detected by GMSDB.

Examples of testing at noisy datasets are shown in fig.5. Here the situation
is much different. The impact of noise on ’middle blobs’ clustering (Fig.5A)
makes classification more difficult with the agglomerative method, but still
allows clustering with the GMSDB algorithm. ’Two horseshoes’ (Fig.5B), ’two
nested rings’ (Fig.5D) and ’three nested rings’ (Fig.5C) clustered by
agglomerative algorithm in the presence of noise cluster more adequately than
without noise.

The last test is ’two snakes’ (Fig.5E) \- the result of the t-SNE transform
(van der Maaten and Hinton, 2008) of the ’two noisy horseshoes’ dataset
(Fig.5B). Two snakes - a dataset of two clusters and noise, which has a
complex shape and an inhomogeneous noise-like structure inside and outside
clusters. The distances between two clusters in different regions are
different and comparable with the characteristic sizes of the clusters
themselves. This dataset cannot be adequately clustered by both algorithms
(agglomerative and GMSDB), but it can be easily seen that the GMSDB algorithm
finds fewer superclusters (ten) than the agglomerative algorithm and detects
two largest superclusters corresponding to ’snakes’ sufficiently well.

Figure 5: Comparison of clustering of different identical datasets by GMSDB
and agglomerative algorithms (minimal interobject distance as intercluster
distance, metric - euclidean). Noisy situation. The maximum allowed number of
clusters is 50. From left to right - the original dataset, superclustering by
the GMSDB algorithm and the found optimal number of superclusters, decision
regions of the trained GMSDB algorithm, clustering by the agglomerative
algorithm. Top right corner in second column images - the number of
superclusters detected by GMSDB.

##  4 Discussion

According to Kleinberg’s impossibility theorem (Kleinberg, 2002), there is no
clustering algorithm that is simultaneously scale-invariant, consistent, and
rich. Obviously, the proposed GMSDB algorithm is unrich \- it does not
enumerate all clustering/superclustering variants. Let us evaluate how it
differs from the rich one in order to understand what possible clusterings it
is limited to. Obviously, if the number of Gaussian distributions is equal to
the number of points in the dataset, then the algorithm ’find the parameters
of all Gaussian distributions, and then combine them to obtain a clustering
labeled by expert’ is almost rich - for discrete points (we usually work with
such datasets), it is enough to choose Gaussian distributions much narrower
than the minimum distance between points. Compared to the full algorithm in
the GMSDB algorithm we limit the number of normal distributions and their
shape: their shape is determined from the EM algorithm by the maximum
likelihood condition, and the number is determined from the minimum of the BIC
criterion. The possible combinations of clusters (superclusters) in the GMSDB
algorithm are also limited - there should be the maximum number of
superclusters separable by the Mahalanobis distance with a significance level
$\alpha$. This limits the set of possible clusterings that the algorithm is
able to determine adequately. In practice, it is not able to correctly
separate non-connected or significantly overlapping clusters, as well as
clusters in which there are not enough elements for the EM algorithm to work
or which are ’not likehood enough’, and could be clustered by different way
with more likehood in terms of Gaussian Mixture.

The high quality of the algorithm is demonstrated in noiseless (Fig.4) and
noisy (Fig.5) situations compared to the agglomerative algorithm.

The advantages of the GMSDB method include also its statistical formulation,
which makes it possible to estimate the probabilities of points belonging to a
particular supercluster, and use it to classify new data both in the hard
clustering mode and in the soft (fuzzy) clustering mode. Indeed, it is known
that a mixture of Gaussians is capable of producing soft clustering - to
estimate for each point the probability of its belonging
$P_{C_{j}}(\overrightarrow{x})$ to one or another Gaussian cluster $C_{j}$.
Each supercluster $SC_{i}$ is a union of the original Gaussian clusters and
the sets of clusters in different superclusters do not intersect. So the
probabilities of any new point belonging to supercluster
$P_{SC_{i}}(\overrightarrow{x})$ can be calculated from the probabilities of
belonging this point to the clusters $C_{j}\subseteq SC_{i}$ that make up this
supercluster:

$P_{SC_{i}}(\overrightarrow{x})=\sum_{j:C_{j}\subseteq
SC_{i}}P_{C_{j}}(\overrightarrow{x})$ (11)

Obviously, the resulting expression (11) is a probability, since it satisfies
Kolmogoroff’s definition of probability (Kolmogoroff, 1933): it can be shown
that if $P_{C_{j}}(\overrightarrow{x})$ is a probability, then for any
$\overrightarrow{x}$:

$\sum_{i}P_{SC_{i}}(\overrightarrow{x})=1$

$P_{SC_{i}\cup SC_{j},i\neq
j}(\overrightarrow{x})=P_{SC_{i}}(\overrightarrow{x})+P_{SC_{j}}(\overrightarrow{x})$

Table 1 shows the time of clustering datasets by the GMSDB and Agglomerative
algorithms (in seconds, by an identical computer), as well as the Rand Index
(RI) determined from the initial labeling of the test datasets during their
creation. In noisy cases, noise was marked as a separate class.

From Table 1 it can be seen that the GMSDB algorithm is most often more
adequate than the agglomerative algorithm, especially on noise-free or very
complex datasets (for example, ’two snakes’), and also more adequate separates
spherical datasets (’medium blobs’) in presence of noise. Higher RI values
indicate higher adequacy.

The main disadvantage of the algorithm is its speed, shown in Table 1. It can
be seen from the table that GMSDB is two orders of magnitude slower than the
agglomerative algorithm, so its use is recommended when one takes into account
its low speed, in the case when the adequacy of clustering is very important,
and there is much enough data for its operation. The algorithm is implemented
in Python, so using fast search algorithms or other programming languages for
these tasks will improve its performance. Its current variant (Berngardt,
2023) has implementations for the speed optimization by a kind of dichotomy
search for stage 1 and Monte-Carlo search for stage 3, significantly improving
clustering speed.

Table 1: The execution time for clustering identical datasets by the GMSDB and Agglomerative algorithms (in seconds), as well as the Rand Index (RI) determined from the initial labeling of test datasets during their creation (in noisy cases, noise is marked as a separate class). Bold indicates the best RI result for given dataset. | GMSDB |  Agglomerative algorithm
---|---|---
Dataset | Stage 1 | Stage 2 |  Stages 3,4 | Total | RI | time | RI
Noise-free datasets (Fig.4)
Grains | 12.16 | 14.03 |  0.002 | 26.20 | 1.000 | 0.44 | 0.993
Big blobs | 12.74 | 10.57 |  0.008 | 23.32 | 0.759 | 0.46 | 0.523
Small blobs | 14.64 | 10.01 |  0.002 | 24.65 | 1.000 | 0.455 | 0.998
Three horseshoes | 33.69 | 34.32 |  0.34 | 68.34 | 1.000 | 0.72 | 0.964
Two horseshoes | 24.02 | 12.85 |  0.10 | 36.98 | 1.000 | 0.36 | 0.935
Three nested rings | 3.24 | 2.17 |  0.15 | 5.56 | 1.000 | 0.05 | 0.881
Two nested rings | 2.76 | 0.91 |  0.02 | 3.69 | 1.000 | 0.025 | 0.797
Noisy datasets (Fig.5)
Medium blobs | 14.39 | 14.14 |  0.01 | 28.54 | 0.986 | 0.53 | 0.527
Two horseshoes | 28.98 | 13.74 |  0.11 | 42.83 | 0.995 | 0.38 | 0.985
Three nested rings | 6.50 | 2.52 |  0.12 | 9.14 | 0.953 | 0.04 | 0.952
Two nested rings | 3.79 | 1.06 |  0.02 | 4.87 | 0.882 | 0.03 | 0.911
Two snakes | 21.75 | 14.67 |  1.17 | 37.59 | 0.743 | 0.40 | 0.630

The next drawback is the presence of a control hyperparameter - the
significance level $\alpha$, but to obtain the results described above, we
used its standard statistical value 0.05.

The last disadvantage of the method is the stochasticity of the results \- the
optimal separation of the same dataset can change from run to run. This is due
to the algorithm is based on the EM implementation of GM, which is based on an
iterative algorithm and depends on the initial conditions.

An interesting feature of the method is the consequence of its stochasticity
\- in some cases, with some probability, it is able to separate even slightly
intersecting clusters. Two examples and their probabilities are shown in
Fig.6. This effect is possibly related to the features of the Mahalanobis
distance and their dependence not only on the distances between the clusters,
but also on the orientation of their main axes. As it can be seen in Fig.2A-D,
sometimes it allows to separate clusters with different orientations of the
main axes in the area of their intersection .

Figure 6: The effect of stochasticity to clustering in the case of
intersecting clusters and the probability of occurrence of these cases (left).
Top left corner in second column images - the number of superclusters detected
by GMSDB.

Fig.7 shows the examples of clustering real images from the (USC-SIPI, )
dataset: Madril,Pirate and Peppers (Fig.7A). For testing, we reduced all the
images to 256x256 size, converted them to grayscale color, and create a
dataset from three-dimensional vectors corresponding to two coordinates of a
pixel and its intensity. The resulting dataset was clustered by the
Agglomerative algorithm (Fig.7B) and the GMSDB algorithm with different
significance levels - 0.05, 0.10, 0.15 (Fig.7C-E). It can be seen from the
figure that the agglomerative algorithm is poorly adapted to distinguish
clusters in such images, and GMSDB in some cases is able to detect well-
interpreted clusters in an image. From the figure it can be seen that GMSDB
with the standard significance level of 0.05 is poorly suited for clustering
such images, while GSMDB with significance levels of 0.10 and 0.15 quite
successfully perform some kind of clustering. This may be due to smoothness of
real images, leading to the difficulty of separating the superclusters
confidently. Increasing the level of significance makes it possible to
consider more similar superclusters as different.

Figure 7: Examples of using GMSDB for superclustering images
(Madril,Pirate,Peppers), converted to greyscale 256x256. A)original image;
B)Agglomerative clustering results; C-E) GMSDB superclustering with different
significance levels $\alpha$: 0.05, 0.10, and 0.15 correspondingly.

##  5 Conclusion

The paper presents a statistical algorithm for GMSDB adequate clustering a
dataset by an optimal grouping of inadequate gaussian models.

The algorithm performs two sequential operations - approximates the set of
points under study by the optimal superposition of normal distributions, each
of them, taken separately, is generally not an adequate cluster in the data.
To create a model of adequate clusters, the algorithm groups the resulting
Gaussian distributions into the maximum number of separable superclusters
(with given statistical significance level). Thus, the algorithm clusters the
dataset into statistically different superclusters, each of which is described
by a superposition of normal distributions.

The algorithm consists of three main stages: approximation of the dataset by a
mixture of Gaussian distributions, the number of which is determined based on
the minimum of the BIC criterion; using the Mahalanobis distance, estimates
the statistical distances between clusters and cluster sizes; combines the
resulting clusters into superclusters using the DBSCAN method by sequential
iterations over its hyperparameter (maximum distance $\varepsilon$) from
smaller to larger values; stops the iterations over the hyperparameter when
the matrix quality criterion MC introduced by us reaches its maximum value of
1. The matrix quality criterion calculates the proportion of statistically
significant separable superclusters among all found superclusters.

Although the algorithm is computationally very slow, it shows good results on
several test datasets, in noise and noiseless situations, as well as for image
clustering. An essential advantage of the method is the ability to predict a
cluster for new data based on an already trained clusterer and the possibility
of soft (fuzzy) clustering. The Python source code of the clusterer is
available in (Berngardt, 2023).

## Acknowledgements

The work has been done under financial support of the Ministry of Science and
Higher Education of the Russian Federation (Subsidy No.075-GZC3569278).

## References

  * Arthur and Vassilvitskii (2007) Arthur, D., Vassilvitskii, S., 2007\.  K-Means++: The Advantages of Careful Seeding, in: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, Society for Industrial and Applied Mathematics, USA. p. 1027–1035.  doi:10.5555/1283383.1283494.
  * Berngardt (2023) Berngardt, O., 2023.  berng/gmsdb: v1.0.  URL: https://zenodo.org/record/8320328, doi:10.5281/ZENODO.8320328.
  * Ester et al. (1996) Ester, M., Kriegel, H.P., Sander, J., Xu, X., 1996.  A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise, in: Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, AAAI Press. p. 226–231.  doi:10.5555/3001460.3001507.
  * Kleinberg (2002) Kleinberg, J., 2002.  An Impossibility Theorem for Clustering, in: Proceedings of the 15th International Conference on Neural Information Processing Systems, MIT Press, Cambridge, MA, USA. pp. 463–470.
  * Kolmogoroff (1933) Kolmogoroff, A., 1933.  Grundbegriffe der Wahrscheinlichkeitsrechnung.  Springer Berlin Heidelberg.  doi:10.1007/978-3-642-49888-6.
  * Lu and Leen (2004) Lu, Z., Leen, T., 2004.  Semi-supervised Learning with Penalized Probabilistic Clustering, in: Saul, L., Weiss, Y., Bottou, L. (Eds.), Advances in Neural Information Processing Systems, MIT Press. p. n/a.
  * van der Maaten and Hinton (2008) van der Maaten, L., Hinton, G., 2008\.  Visualizing Data using t-SNE.  Journal of Machine Learning Research 9, 2579–2605.
  * Mahalanobis (1936) Mahalanobis, P., 1936.  On the Generalised Distance in Statistics.  Proceedings of the National Institute of Sciences of India 2, 49–55.
  * Ng et al. (2001) Ng, A., Jordan, M., Weiss, Y., 2001.  On Spectral Clustering: Analysis and an algorithm, in: Dietterich, T., Becker, S., Ghahramani, Z. (Eds.), Advances in Neural Information Processing Systems, MIT Press. p. n/a.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011\.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825–2830.
  * Rand (1971) Rand, W.M., 1971.  Objective Criteria for the Evaluation of Clustering Methods.  Journal of the American Statistical Association 66, 846–850.  doi:10.1080/01621459.1971.10482356.
  * Reynolds (2015) Reynolds, D., 2015.  Gaussian Mixture Models. Springer US, Boston, MA. chapter n/a.  pp. 827–832.  doi:10.1007/978-1-4899-7488-4\\_196.
  * Rousseeuw (1987) Rousseeuw, P.J., 1987.  Silhouettes: A graphical aid to the interpretation and validation of cluster analysis.  Journal of Computational and Applied Mathematics 20, 53–65.  doi:10.1016/0377-0427(87)90125-7.
  * Schubert et al. (2017) Schubert, E., Sander, J., Ester, M., Kriegel, H.P., Xu, X., 2017.  DBSCAN Revisited, Revisited: Why and How You Should (Still) Use DBSCAN.  ACM Trans. Database Syst. 42.  doi:10.1145/3068335.
  * Schwarz (1978) Schwarz, G., 1978.  Estimating the Dimension of a Model.  The Annals of Statistics 6, 461 – 464.  doi:10.1214/aos/1176344136.
  * Shiming et al. (2008) Shiming, X., Feiping, N., Changshui, Z., 2008.  Learning a Mahalanobis distance metric for data clustering and classification.  Pattern Recognition 41, 3600–3612.  doi:10.1016/j.patcog.2008.05.018.
  * Sibson (1973) Sibson, R., 1973.  SLINK: An optimally efficient algorithm for the single-link cluster method.  The Computer Journal 16, 30–34.  doi:10.1093/comjnl/16.1.30.
  * Student (1908) Student, 1908.  The probable error of a mean.  Biometrika 6, 1\.  doi:10.2307/2331554.
  * (19) USC-SIPI, .  The USC-SIPI Image Database.  https://sipi.usc.edu/database/database.php.


# Modulated dust-acoustic wave packets in an opposite polarity dusty plasma
system

S. Jahana, N. A. Chowdhury, A. Mannan, and A. A. Mamun  Department of Physics,
Jahangirnagar University, Savar, Dhaka-1342, Bangladesh.
∗Email<EMAIL_ADDRESS>
###### Abstract

The nonlinear propagation of the dust-acoustic bright and dark envelope
solitons in an opposite polarity dusty plasma system (composed of non-
extensive $q$-distributed electrons, isothermal ions, and positively as well
as negatively charged warm dust) have been theoretically investigated. The
reductive perturbation method (which is valid for a small, but finite
amplitude limit) is employed to derive the nonlinear Schrödinger equation. Two
types of modes (namely, fast and slow dust-acoustic (DA) modes) have been
observed. The conditions for the modulational instability (MI) and its growth
rate in the unstable regime of the DA waves are significantly modified by the
effects of non-extensive electrons, dust mass, temperatures of different
plasma species, etc. The implications of the obtained results from our current
investigation in space (viz. Jupiters magnetosphere, upper mesosphere, and
comets tails) and laboratory (e.g. direct current and radio-frequency
discharges, plasma processing reactors, fusion plasma devices, and solid-fuel
combustion products, etc.) dusty plasmas are briefly discussed.

##  I Lead Paragraph

The nonlinear propagation of the dust-acoustic bright and dark envelope
solitons in an opposite polarity dusty plasma system (containing isothermal
ions, non-extensive $q$-distributed electrons, and positively as well as
negatively charged warm dust) has been studied to identify some new features
(e.g. conditions for the modulational instabiliuty of the dust-acoustic waves
(DAWs), growth rate, formation of the bright and dark envelope solitons, etc.)
of the DAWs. These new features should be useful for understanding the
nonlinear electrostatic disturbances in space and laboratory dusty plasmas.

##  II Introduction

Now-a-days, a great deal of interest has been devoted to the new and
fascinating field of dusty plasma due to their existence in space Shukla2002 ;
Verheest1996a ; Goertz1986 ; Whipple1985  and massive applications in
laboratory plasmas Mamun1996a ; Rao1990 ; Angelo1990 . The dusty plasma is
generally considered to be an ensemble of dust particles, free electrons, and
ions, which has been found widely in planetary rings, asteroid zones Rao1990 ,
cometary tails, magnetosphere as well as the lower part of the earth’s
atmosphere Goertz1986 ; Whipple1985 . However, the co-existence of opposite
polarity dusty in plasmas, introduces a new dusty plasma model called
“opposite polarity dusty plasma” (OPDP) Mamun2016b , whose main species are
positively and negatively charged warm massive dust Mamun2016b . The exclusive
property of this OPDP, which makes it completely unique from other plasmas
(viz. electron ion and electron-positron plasmas) is that the ratio of the
size of positively charged dust to that of negatively charged dust can be
smaller Chow1993  or larger Mendis1994a  or equal to unity Mendis1994a . Dust
species can be charged either positively or negatively, because of different
charging processes, like photo-ionization, secondary electron emission,
thermionic emission, and collection of plasma particles (electrons and ions),
etc. The existence of positively charged dust has been observed in different
regions of space (e.g. Jupiter’s magnetosphere Horányi1993a , cometary tails
Mendis1994a ; Mendis1991b ; Horányi1996b ; Chow1993 , upper mesosphere
Havnes1996 , etc). There are three main processes by which dust grains become
positively charged, namely, Fortov1998 ; Mamun2008c  (a) secondary emission of
electrons from the surface of the dust grains; (b) thermionic emission induced
by the radiative heating; and (c) photoemission in the presence of a flux of
ultraviolet photons. There is also direct evidence for the co-existence of
both positively and negatively charged dust particles in different regions of
space plasmas (viz. cometary tails Horányi1996b ; Mendis1991b , upper
mesosphere Havnes1996 , Jupiter’s magnetosphere Horányi1996b ; Horányi1993a ;
Mamun2008c ; Chowdhury2017a , etc.) and laboratory devices (e.g. direct
current and radio-frequency discharges, plasma processing reactors, fusion
plasma devices, solid-fuel combustion products, etc Shukla2002 ; Akhter2012 ).
In space and astrophysical sectors, if the plasma particles move very fast
compared to their thermal velocities Renyi1955  then the Maxwellian
distribution is no longer valid. For that reason, Tsallis proposed the non-
extensive statistics Tsallis1988 ; Driouch2017 , which is the generalisation
of Boltzmann-Gibbs-Shannon entropy. The importance of Tsallis statistics lies
with describing the system of long range interactions, such as, plasmas and
dusty plasmas Driouch2017 .

The researchers are focused on wave dynamics, specifically, dust-acoustic (DA)
waves (DAWs), dust-acoustic rogue waves (DARWs), and dust ion-acoustic waves
(DIAWs) in understanding electrostatic density perturbations and potential
structures (viz. soliton, shock, vortices, and rogue waves Chowdhury2017a ;
Dubinov2009 ; Gill2010 ; Tasnim2013 ). The DAWs Rao1990  have low phase
velocity (where the inertia is provided by the dust mass and restoring force
is provided by the thermal pressure of electrons and ions Chowdhury2017a ) in
comparison with the electron and ion thermal velocities. The research on
modulational instability (MI) of DAWs in nonlinear and dispersive mediums have
been increasing significantly due to their existence in astrophysics, space
physics Chowdhury2018b  as well as in application in many laboratory
situations Misra2006 . A great number of researchers have used the nonlinear
Schrödinger (NLS) equation, which governs the dynamics of the DAWs and the
formation of the envelope solitons Misra2006 . Sayed and Mamun Sayed2007
studied the finite solitary potential structures that exist in OPDP. El-
Taibany El-Taibany2013  examined the DAWs in inhomogeneous four-component
dusty plasmas with opposite charge polarity dust grains and observed that only
compressive soliton is created corresponding to fast DAWs velocities. To the
best knowledge of the authors, no attempt has been made on MI, and
corresponding dark and bright envelope solitons associated with the DAWs
containing non-extensive $q$-distributed electrons, isothermal ions, and
positively and negatively charged warm massive dust. Therefore, in our
article, we will derive the NLS equation by employing the reductive
perturbation method and examine the conditions for the MI of the DAWs (in
which inertia is provided by the dust mass and restoring force is provided
from the thermal pressure of non-extensive $q$-distributed electrons and
isothermal ions).

The manuscript is organized as the following fashion: The governing equations
of our considered plasma model are stated in Sec. III. The NLS equation is
derived in Sec. IV. The stability of DAWs and envelope solitons are examined
in Sec. V and VI, respectively. Finally, the summary of our discussion is
provided in Sec. VII.

##  III Model Equations

In this paper, we consider a collisionless, fully ionized, unmagnetized four
component dusty plasma system formed of non-extensive q-distributed electrons
(charge $-e$, mass $m_{e}$), isothermal ions (charge $+e$; mass $m_{i}$) and
inertial warm negatively charged dust particles (charge $q_{1}=-z_{1}e$, mass
$m_{1}$) as well as positively charged warm dust particles (charge
$q_{2}=+z_{2}e$; mass $m_{2}$); where $z_{1}$ ($z_{2}$) is the charge state of
the negatively (positively) charged warm dust particles. The negatively and
positively charged warm dust particles can be displayed by continuity and
momentum equations, respectively, as given below:

$\displaystyle\frac{\partial n_{1}}{\partial t}+\frac{\partial}{\partial
x}(n_{1}u_{1})=0,$ (1) $\displaystyle\frac{\partial u_{1}}{\partial
t}+u_{1}\frac{\partial u_{1}}{\partial
x}=\frac{z_{1}e}{m_{1}}\frac{\partial\varphi}{\partial
x}-\frac{1}{m_{1}n_{1}}\frac{\partial p_{1}}{\partial x},$ (2)
$\displaystyle\frac{\partial n_{2}}{\partial t}+\frac{\partial}{\partial
x}(n_{2}u_{2})=0,$ (3) $\displaystyle\frac{\partial u_{2}}{\partial
t}+u_{2}\frac{\partial u_{2}}{\partial
x}=-\frac{z_{2}e}{m_{2}}\frac{\partial\varphi}{\partial
x}-\frac{1}{m_{2}n_{2}}\frac{\partial p_{2}}{\partial x},\ $ (4)

where $n_{1}(n_{2})$ is the number densities of the negatively (positively)
charged warm dust particles; $t(x)$ is the time (space) variable;
$u_{1}(u_{2})$ is the fluid speed of the negatively (positively) charged warm
dust species; $e$ is the magnitude of the charge of the electron; $\varphi$ is
the electrostatic wave potential; and $p_{1}$ ($p_{2}$) is the adiabatic
pressure of the negatively (positively) charged warm dust particles. The
system is closed through Poisson’s equation, stated as

$\displaystyle\frac{\partial^{2}\varphi}{\partial x^{2}}=4\pi
e(n_{e}-n_{i}+z_{1}n_{1}-z_{2}n_{2}),$ (5)

where $n_{i}$ and $n_{e}$ are, respectively, the ion and electron number
densities. The quasi-neutrality condition at equilibrium can be shown as

$\displaystyle n_{i0}+z_{2}n_{20}=n_{e0}+z_{1}n_{10},$ (6)

where $n_{i0}$, $n_{20}$, $n_{e0}$, and $n_{10}$ are the equilibrium number
densities of the isothermal ions, positively charged warm dust particles, non-
extensive $q$-distributed electrons, and negatively charged warm dust
particles, respectively. Now, in terms of normalized variables, namely,
$N_{1}=n_{1}/n_{10}$, $N_{2}=n_{2}/n_{20}$, $U_{1}=u_{1}/C_{d1}$ (with
$C_{d1}$ being the sound speed of the negatively charged warm dust particles),
$U_{2}=u_{2}/C_{d1}$; $\phi=e\varphi/T_{i}$ (with $T_{i}$ being the
temperature of the isothermal ion); $T=t\omega_{pd1}$ (with $\omega_{pd1}$
being the plasma frequency of the negatively charged warm dust particles),
$X=x/\lambda_{Dd1}$ (with $\lambda_{Dd1}$ being the Debye length of the
negatively charged warm dust species); $C_{d1}=(z_{1}T_{i}/m_{1})^{1/2}$,
$\omega_{pd1}=(4\pi e^{2}z_{1}^{2}n_{10}/m_{1})^{1/2}$;
$\lambda_{Dd1}=(T_{i}/4\pi e^{2}z_{1}n_{10})^{1/2}$;
$p_{1}=p_{10}(n_{1}/n_{10})^{\gamma}$ [with $p_{10}$ being the equilibrium
adiabatic pressure of the negatively charged warm dust particles and
$\gamma=(N+2)/N$, where $N$ is the degree of freedom, for one-dimensional
case, $N=1$ so that $\gamma=3$], $p_{10}=n_{10}T_{1}$ (with $T_{1}$ being the
temperature of the negatively charged warm dust particles);
$p_{2}=p_{20}(n_{2}/n_{20})^{3}$ (with $p_{20}$ being the equilibrium
adiabatic pressure of the positively charged warm dust particles), and
$p_{20}=n_{20}T_{2}$ (with $T_{2}$ being the temperature of the positively
charged warm dust particles). After normalization, Eqs. (1)$-$(5) becomes

$\displaystyle\frac{\partial N_{1}}{\partial T}+\frac{\partial}{\partial
X}(N_{1}U_{1})=0,$ (7) $\displaystyle\frac{\partial U_{1}}{\partial
T}+U_{1}\frac{\partial U_{1}}{\partial X}+3\sigma_{1}N_{1}\frac{\partial
N_{1}}{\partial X}=\frac{\partial\phi}{\partial X},$ (8)
$\displaystyle\frac{\partial N_{2}}{\partial T}+\frac{\partial}{\partial
X}(N_{2}U_{2})=0,$ (9) $\displaystyle\frac{\partial U_{2}}{\partial
T}+U_{2}\frac{\partial U_{2}}{\partial X}+3\sigma_{2}N_{2}\frac{\partial
N_{2}}{\partial X}=-\alpha\frac{\partial\phi}{\partial X},$ (10)
$\displaystyle\frac{\partial^{2}\phi}{\partial
X^{2}}=(\mu_{i}+\beta-1)N_{e}-\mu_{i}N_{i}+N_{1}-\beta N_{2},$ (11)

where $\sigma_{1}=T_{1}/z_{1}T_{i}$, $\sigma_{2}=m_{1}T_{2}/z_{1}m_{2}T_{i}$,
$\alpha=m_{1}z_{2}/m_{2}z_{1}$, $\beta=z_{2}n_{20}/z_{1}n_{10}$, and
$\mu_{i}=n_{i0}/z_{1}n_{10}$. The number densities of the non-extensive
q-distributed Chowdhury2017a  electron can be given by the following
normalized equation

$\displaystyle N_{e}=[1+(q-1)\delta\phi]^{\frac{q+1}{2(q-1)}},$ (12)

where $\delta=T_{i}/T_{e}$ and $q$ is the non-extensive parameter describing
the degree of non-extensivity, i.e., $q=1$ indicates the Maxwellian
distribution, whereas $q<1$ refers to the super-extensivity, and the opposite
condition $q>1$ corresponds to the sub-extensivity Chowdhury2017c . The number
densities of the iso-thermally distributed Chowdhury2017a  ion can be
represented as

$\displaystyle N_{i}=\mbox{exp}(-\phi).$ (13)

Now, substituting Eqs. (12) and (13) into Eq. (11), and extending up to third
order in $\phi$, we can write

$\displaystyle\frac{\partial^{2}\phi}{\partial X^{2}}+\beta
N_{2}-N_{1}=\beta-1+\gamma_{1}\phi+\gamma_{2}\phi^{2}+\gamma_{3}\phi^{3}+..,$
(14)

where

$\displaystyle\gamma_{1}=\frac{(\beta+\mu_{i}-1)(q+1)\delta}{2}+\mu_{i},$
$\displaystyle\gamma_{2}=\frac{(\beta+\mu_{i}-1)(q+1)(3-q)\delta^{2}}{8}-\frac{\mu_{i}}{2},$
$\displaystyle\gamma_{3}=\frac{(\beta+\mu_{i}-1)(q+1)(q-3)(3q-5)\delta^{3}}{48}+\frac{\mu_{i}}{6}.\
$

The left hand side of Eq. (14), is the contribution of electron and ion
species.

##  IV Derivation of the NLS equation

We will use the reductive perturbation method to derive the NLS equation to
know the modulation of the DAWs. Let us consider, the stretched co-ordinate
Chowdhury2017a  as

$\displaystyle\xi=\epsilon(X-V_{g}T),$ (15) $\displaystyle\tau=\epsilon^{2}T,$
(16)

where $V_{g}$ is the envelope group velocity and $\epsilon$ ($0<\epsilon<1$)
is a small (real) parameter. The dependent variables Chowdhury2017a  can be
written as

$\displaystyle
N_{1}=1+\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l=-\infty}^{\infty}N_{1l}^{(m)}(\xi,\tau)~{}\mbox{exp}(il\Upsilon),$
(17) $\displaystyle
U_{1}=\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l=-\infty}^{\infty}U_{1l}^{(m)}(\xi,\tau)~{}\mbox{exp}(il\Upsilon),$
(18) $\displaystyle
N_{2}=1+\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l=-\infty}^{\infty}N_{2l}^{(m)}(\xi,\tau)~{}\mbox{exp}(il\Upsilon),$
(19) $\displaystyle
U_{2}=\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l=-\infty}^{\infty}U_{2l}^{(m)}(\xi,\tau)~{}\mbox{exp}(il\Upsilon),$
(20)
$\displaystyle\phi=\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l=-\infty}^{\infty}\phi_{l}^{(m)}(\xi,\tau)~{}\mbox{exp}(il\Upsilon),$
(21)

where $\Upsilon=kX-\omega T$ and $k$($\omega$) is the carrier wave number
(frequency). The derivative operators in the above equations are considered as
follows:

$\displaystyle\frac{\partial}{\partial T}\rightarrow\frac{\partial}{\partial
T}-\epsilon
V_{g}\frac{\partial}{\partial\xi}+\epsilon^{2}\frac{\partial}{\partial\tau},$
(22) $\displaystyle\frac{\partial}{\partial
X}\rightarrow\frac{\partial}{\partial
X}+\epsilon\frac{\partial}{\partial\xi}.$ (23)

Now, by substituting Eqs. (15)$-$(23) into Eqs. (7)$-$(10), and (14), and
collecting power term of $\epsilon$, the first order approximation ($m=1$)
with the first harmonic ($l=1$) provides the following relation

$\displaystyle ikU_{11}^{(1)}-i\omega N_{11}^{(1)}=0,$ (24) $\displaystyle
ik\lambda N_{11}^{(1)}-i\omega U_{11}^{(1)}-ik\phi_{1}^{(1)}=0,$ (25)
$\displaystyle ikU_{21}^{(2)}-i\omega N_{21}^{(2)}=0,$ (26) $\displaystyle
ik\theta N_{21}^{(1)}+ik\alpha\phi_{1}^{(1)}-i\omega U_{21}^{(1)}=0,$ (27)
$\displaystyle\beta
N_{21}^{(1)}-N_{11}^{(1)}-k^{2}\phi_{1}^{(1)}-\gamma_{1}\phi_{1}^{(1)}=0,$
(28)

where $\lambda=3\sigma_{1}$ and $\theta=3\sigma_{2}$. Now, these equations
detracted to the following pattern

$\displaystyle N_{11}^{(1)}=\frac{k^{2}}{S}\phi_{1}^{(1)},$ (29)
$\displaystyle U_{11}^{(1)}=\frac{k\omega}{S}\phi_{1}^{(1)},$ (30)
$\displaystyle N_{21}^{(2)}=\frac{k^{2}\alpha}{A}\phi_{1}^{(1)},$ (31)
$\displaystyle U_{21}^{(1)}=\frac{k\omega\alpha}{A}\phi_{1}^{(1)},$ (32)

where $A=\omega^{2}-\theta k^{2}$ and $S=\lambda k^{2}-\omega^{2}$. Therefore,
the dispersion relation for the DAWs can be written as

$\displaystyle\omega^{2}=\frac{k^{2}G\pm k^{2}\sqrt{G^{2}-4HM}}{2H},$ (33)

where $G=(\theta k^{2}+\lambda
k^{2}+\theta\gamma_{1}+\lambda\gamma_{1}+\alpha\beta+1),H=(k^{2}+\gamma_{1})$,
and $M=(\theta\lambda
k^{2}+\theta\gamma_{1}\lambda+\theta+\alpha\beta\lambda)$. The condition
$G^{2}>4HM$ must be satisfied in order to obtain real and positive values of
$\omega$. Normally, two types of DA modes exist, namely, fast ($\omega_{f}$)
and slow ($\omega_{s}$) DA modes according to the positive and negative sign
of the Eq. (33). Now, we have studied the dispersion properties by depicting
$\omega$ with $k$ in Figs. 2 and 2, which clearly indicates that (a) the value
of $\omega_{f}$ increases exponentially with the increasing values of $z_{2}$
for fixed value of $z_{1}$, $n_{20}$, and $n_{10}$ although, a saturation
region begins after a certain value of $k$ (as is clearly shown in Fig. 2);
(b) on the other side, the magnitude of $\omega_{s}$ linearly increases with
the increase of $z_{2}$ for the fixed value of $z_{1}$, $n_{20}$, and $n_{10}$
(as is clearly shown in Fig. 2). It is important to mention that in fast DA
mode, both positive and negative warm dust species oscillate in phase with
electrons and ions, whereas, in slow DA mode, only one of the inertial massive
dust components oscillate in phase with electrons and ions, but the other
species are in anti-phase with them Chowdhury2017a . Next, with the help of
second-order ($m=2$ with $l=1$) equations, we obtain the expression of $V_{g}$
like that

$\displaystyle V_{g}=\frac{F_{1}}{2\omega k(A^{2}+\alpha\beta S^{2})},$ (34)
$\displaystyle F_{1}=\lambda A^{2}k^{2}+\alpha\beta\theta
k^{2}S^{2}+\omega^{2}A^{2}+\alpha\beta\omega^{2}S^{2}$
$\displaystyle-2A^{2}S^{2}-SA^{2}+\alpha\beta AS^{2}.\ $

Now, the second-harmonic mode of the carrier comes from nonlinear self
interaction caused by the components ($l=2$) for the second order ($m=2$)
compressed equations in the following form

$\displaystyle N_{12}^{(2)}=C_{1}|\phi_{1}^{(1)}|^{2},$ (35) $\displaystyle
U_{12}^{(2)}=C_{2}|\phi_{1}^{(1)}|^{2},$ (36) $\displaystyle
N_{22}^{(2)}=C_{3}|\phi_{1}^{(1)}|^{2},$ (37) $\displaystyle
U_{22}^{(2)}=C_{4}|\phi_{1}^{(1)}|^{2},$ (38)
$\displaystyle\phi_{2}^{(2)}=C_{5}|\phi_{1}^{(1)}|^{2},$ (39)

where

$\displaystyle C_{1}=\frac{2C_{5}k^{2}S^{2}-(3\omega^{2}k^{4}+\lambda
k^{6})}{2S^{3}},$ $\displaystyle C_{2}=\frac{\omega C_{1}S^{2}-\omega
k^{4}}{kS^{2}},$ $\displaystyle
C_{3}=\frac{3\alpha^{2}\omega^{2}k^{4}+\theta\alpha^{2}k^{6}+2\alpha
C_{5}A^{2}k^{2}}{2A^{3}},$ $\displaystyle C_{4}=\frac{\omega
C_{3}A^{2}-\omega\alpha^{2}k^{4}}{kA^{2}},$ $\displaystyle
C_{5}=\frac{F_{2}}{F_{3}},$ $\displaystyle F_{2}=3\omega^{2}A^{3}k^{4}+\lambda
A^{3}k^{6}-2\gamma_{2}A^{3}S^{3}+3\beta\alpha^{2}\omega^{2}S^{3}k^{4}$
$\displaystyle\hskip 22.76228pt+\beta\theta\alpha^{2}S^{3}k^{6},$
$\displaystyle
F_{3}=2k^{2}S^{2}A^{3}+8k^{2}A^{3}S^{3}+2\gamma_{1}A^{3}S^{3}-2\alpha\beta
A^{2}k^{2}S^{3}.$

After that, we consider the expression for ($m=3,l=0$) and ($m=2,l=0$), which
leads the zeroth harmonic modes. Finally, we get

$\displaystyle N_{10}^{(2)}=C_{6}|\phi_{1}^{(1)}|^{2},$ (40) $\displaystyle
U_{10}^{(2)}=C_{7}|\phi_{1}^{(1)}|^{2},$ (41) $\displaystyle
N_{20}^{(2)}=C_{8}|\phi_{1}^{(1)}|^{2},$ (42) $\displaystyle
U_{20}^{(2)}=C_{9}|\phi_{1}^{(1)}|^{2},$ (43)
$\displaystyle\phi_{0}^{(2)}=C_{10}|\phi_{1}^{(1)}|^{2},$ (44)

where

$\displaystyle C_{6}=\frac{2\omega V_{g}k^{3}+\lambda
k^{4}+\omega^{2}k^{2}-C_{10}S^{2}}{S^{2}(V_{g}^{2}-\lambda)},$ $\displaystyle
C_{7}=\frac{C_{6}V_{g}S^{2}-2\omega k^{3}}{S^{2}},$ $\displaystyle
C_{8}=\frac{2\omega
V_{g}\alpha^{2}k^{3}+\theta\alpha^{2}k^{4}+\alpha^{2}\omega^{2}k^{2}+\alpha
C_{10}A^{2}}{A^{2}(V_{g}^{2}-\theta)},$ $\displaystyle
C_{9}=\frac{C_{8}V_{g}A^{2}-2\omega\alpha^{2}k^{3}}{A^{2}},$ $\displaystyle
C_{10}=\frac{F_{4}}{F_{5}},$ $\displaystyle F_{4}=(A^{2}V_{g}^{2}-\theta
A^{2})(2\omega V_{g}k^{3}+\lambda k^{4}+\omega^{2}k^{2})$ $\displaystyle\hskip
5.69046pt+(V_{g}^{2}-\lambda)(2\gamma_{2}A^{2}S^{2}V_{g}^{2}-2\gamma_{2}\theta
A^{2}S^{2})$ $\displaystyle\hskip 5.69046pt-(2\omega
V_{g}\alpha^{2}k^{3}+\theta\alpha^{2}k^{4}+\alpha^{2}\omega^{2}k^{2})$
$\displaystyle\hskip 5.69046pt(\beta S^{2}V_{g}^{2}-\beta\lambda S^{2}),$
$\displaystyle F_{5}=\alpha\beta A^{2}S^{2}V_{g}^{2}-\alpha\beta\lambda
A^{2}S^{2}+A^{2}S^{2}V_{g}^{2}-\theta A^{2}S^{2}$ $\displaystyle\hskip
5.69046pt-(V_{g}^{2}-\lambda)(\gamma_{1}A^{2}S^{2}V_{g}^{2}-\gamma_{1}\theta
A^{2}S^{2}).$

Later, with the help of Eqs. (29)$-$(44), the third harmonic ($m=3$) with
($l=1$) modes can be accomplished. Thus, we can write the NLS equation as;

$\displaystyle
i\frac{\partial\Phi}{\partial\tau}+P\frac{\partial^{2}\Phi}{\partial\xi^{2}}+Q|\Phi|^{2}\Phi=0,$
(45)

where $\Phi=\phi_{1}^{(1)}$ for simplicity and $P$ ($Q$) is the dispersion
(nonlinear) coefficient, and is written by

$\displaystyle P=\frac{F_{6}}{2\omega ASk^{2}(A^{2}+\alpha\beta S^{2})},\ $

and

$\displaystyle Q=\frac{F_{7}}{2\omega k^{2}(A^{2}+\alpha\beta S^{2})},\ $

where

$\displaystyle F_{6}=(\omega V_{g}A^{3}-\lambda kA^{3})(\lambda k^{3}-2\omega
V_{g}k^{2}+k\omega^{2}-kS)$ $\displaystyle\hskip
11.38092pt+(kV_{g}A^{3}-\omega A^{3})(\lambda\omega
k^{2}-2kV_{g}\omega^{2}+\omega^{3}-kSV_{g})$ $\displaystyle\hskip
11.38092pt-(\theta k^{3}-2\omega V_{g}k^{2}+k\omega^{2}+Ak)(\alpha\beta\omega
V_{g}S^{3}$ $\displaystyle\hskip 11.38092pt-\alpha\beta\theta
kS^{3})-A^{3}S^{3}-(\alpha\beta kV_{g}S^{3}-\alpha\beta\omega S^{3})$
$\displaystyle\hskip 11.38092pt(\theta\omega
k^{2}-2kV_{g}\omega^{2}+\omega^{3}+AkV_{g}),\ $ $\displaystyle
F_{7}=2\gamma_{2}C_{5}A^{2}S^{2}+2\gamma_{2}C_{10}A^{2}S^{2}+3\gamma_{3}A^{2}S^{2}$
$\displaystyle-2\omega C_{2}A^{2}k^{3}-2\omega
C_{7}A^{2}k^{3}-2\alpha\beta\omega C_{4}S^{2}k^{3}$
$\displaystyle-2\alpha\beta\omega
C_{9}S^{2}k^{3}-C_{1}\omega^{2}A^{2}k^{2}-C_{6}\omega^{2}A^{2}k^{2}$
$\displaystyle-\lambda C_{1}A^{2}k^{4}-\lambda C_{6}A^{2}k^{4}-\alpha\beta
C_{3}\omega^{2}k^{2}S^{2}$ $\displaystyle-\alpha\beta
C_{8}\omega^{2}k^{2}S^{2}-\alpha\beta\theta C_{3}S^{2}k^{4}-\alpha\beta\theta
C_{8}S^{2}k^{4}.$

Figure 1: The variation of $\omega_{f}$ with $k$ for different values of
$\beta$; along with $\alpha=2.0$, $\delta=0.3$, $\mu_{i}=0.4$,
$\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, and $q=1.5$.

Figure 2: The variation of $\omega_{s}$ with $k$ for different values of
$\beta$; along with $\alpha=2.0$, $\delta=0.3$, $\mu_{i}=0.4$,
$\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, and $q=1.5$.

Figure 3: The variation of $P/Q$ with $k$ for different values of $\mu_{i}$;
along with $\alpha=2.0$, $\beta=0.07$, $\delta=0.3$, $\sigma_{1}=0.0001$,
$\sigma_{2}=0.001$, $\omega_{f}$, and $q=1.5$.

Figure 4: The variation of $P/Q$ with $k$ for different values of $\alpha$;
along with $\beta=0.07$, $\delta=0.3$, $\mu_{i}=0.4$, $\sigma_{1}=0.0001$,
$\sigma_{2}=0.001$, $\omega_{s}$, and $q=1.5$.

Figure 5: The variation of MI growth rate $(\Gamma)$ with $\tilde{k}$ for
different values of $\alpha$; along with $\beta=0.07$, $\delta=0.3$,
$\mu_{i}=0.4$, $\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\omega_{f}$, $q=1.0$,
$\Phi_{0}=0.5$, and $k=0.3$.

Figure 6: The variation of MI growth rate $(\Gamma)$ with $\tilde{k}$ for
different values of $\alpha$; along with $\beta=0.07$, $\delta=0.3$,
$\mu_{i}=0.4$, $\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\omega_{f}$, $q=1.5$,
$\Phi_{0}=0.5$, and $k=0.3$.

Figure 7: The variation of MI growth rate $(\Gamma)$ with $\tilde{k}$ for
different values of $\alpha$; along with $\beta=0.07$, $\delta=0.3$,
$\mu_{i}=0.4$, $\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\omega_{f}$,
$q=-0.6$, $\Phi_{0}=0.5$, and $k=0.3$.

Figure 8: The variation of MI growth rate $(\Gamma)$ with $\tilde{k}$ for $q$
(within sub-extensive range); along with $\beta=0.07$, $\delta=0.3$,
$\mu_{i}=0.4$, $\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\omega_{f}$,
$\Phi_{0}=0.5$, and $k=0.3$.

Figure 9: The variation of MI growth rate $(\Gamma)$ with $\tilde{k}$ for $q$
(within super-extensive range); along with $\beta=0.07$, $\delta=0.3$,
$\mu_{i}=0.4$, $\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\omega_{f}$,
$\Phi_{0}=0.5$, and $k=0.3$.

Figure 10: The variation of $Re(\Phi)$ with $\xi$ for bright envelope
solitons; along with $\alpha=2.0$, $\beta=0.07$, $\delta=0.3$, $\mu_{i}=0.4$,
$\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\tau=0$, $\psi_{0}=0.008$,
$\omega_{f}$, $q=1.5$, $k=0.3$, $\Omega_{0}=0.4$, and $U=0.4$.

Figure 11: The variation of $Re(\Phi)$ with $\xi$ for dark envelope solitons;
along with $\alpha=2.0$, $\beta=0.07$, $\delta=0.3$, $\mu_{i}=0.4$,
$\sigma_{1}=0.0001$, $\sigma_{2}=0.001$, $\tau=0$, $\psi_{0}=0.008$,
$\omega_{f}$, $q=1.5$, $k=0.1$, $\Omega_{0}=0.4$, and $U=0.4$.

##  V Stability of DAWs

To study the MI of DAWs, we consider the linear solution of the NLS equation
(45) in the form $\Phi=\hat{\Phi}e^{iQ|\hat{\Phi}|^{2}\tau}+c.c$ (c.c denotes
the complex conjugate), where
$\hat{\Phi}=\hat{\Phi}_{0}+\epsilon\hat{\Phi}_{1}$ and
$\hat{\Phi}_{1}=\hat{\Phi}_{1,0}e^{i(\tilde{k}\xi-\tilde{\omega}\tau)}$+c.c
(the perturbed wave number $\tilde{k}$ and the frequency $\tilde{\omega}$ are
different from $k$ and $\omega$). Now, by substituting these values in Eq.
(45) the following nonlinear dispersion relation can be obtained as
Sultana2011 ; Kourakis2005 ; Schamel2002 ; Fedele2002

$\displaystyle\tilde{\omega}^{2}=P^{2}\tilde{k}^{2}\left(\tilde{k}^{2}-\frac{2|\hat{\Phi}_{0}|^{2}}{P/Q}\right).$
(46)

It is clear from Eq. (46) that the DAWs will be modulationally stable
(unstable) in the range of values of $\tilde{k}$ in which $P/Q$ is negative
(positive). When $P/Q\rightarrow\pm\infty$, the corresponding value of $k$
($=k_{c}$) is called the critical or threshold wave number ($k_{c}$) for the
onset of MI. The variation of $P/Q$ with $k$ for $\mu_{i}$ and $\alpha$ are
shown in Figs. 4 and 4, respectively, which clearly indicate that (a) the
value of $k_{c}$ increases with the increase of $n_{i0}$ for fixed value of
$z_{1}$ and $n_{10}$; (b) on the other hand, $k_{c}$ value decreases with the
increase of $m_{2}$ for fixed value of $m_{1}$, $z_{2}$, and $z_{1}$. It also
can be written from Eq. (46) that the growth rate ($\Gamma$) of the
modulationally unstable region for the DAWs [when $P/Q>0$ and
$\tilde{k}<\tilde{k}_{c}=(2Q|\hat{\Phi}_{0}|^{2}/P)^{1/2}$] as

$\displaystyle\Gamma=|P|\tilde{k}^{2}\left(\frac{\tilde{k}^{2}_{c}}{\tilde{k}^{2}}-1\right)^{1/2}.$
(47)

Moreover, we have graphically shown how the $\Gamma$ varies with $\tilde{k}$
of different values of $\alpha$ and $q$ in Figs. 7$-$9.

It is obvious from Figs. 7$-$7 that (a) the maximum value of $\Gamma$
increases with the increase of the values of $z_{2}$ for fixed value of $q=1$
(Maxwellian limit), $z_{1}$, $m_{1}$, and $m_{2}$ (as is clearly shown in Fig.
7); (b) the magnitude of $\Gamma$ also increases with the increase of $z_{2}$
for fixed value of $q=1.5$ (subextensive range), $z_{1}$, $m_{1}$, and $m_{2}$
(as is clearly shown in Fig. 7); (c) Moreover, $\Gamma$ increases with the
increase of $z_{2}$ for fixed value of $q=-0.6$ (super-extensive range),
$z_{1}$, $m_{1}$, and $m_{2}$ (as clearly shown in Fig. 7). The physics of
this result is that, since the nonlinearity increases with the increase of the
values of $\alpha$, the growth rate of DAWs increase.

It can be observed from Figs. 9 and 9 that the maximum value of $\Gamma$
increases (decreases) with the decrease (increase) of the values of $q>1$
($q<1$) (as is clearly shown in Fig. 9 and 9). It can be concluded here that
the variation of $\Gamma$ with respect to $\tilde{k}$ is independent on the
sign of the $q$.

##  VI Envelope solitons

There are two types of envelope solitonic solutions exist, namely, bright and
dark envelope solitons, depending on the sign of the coefficients $P$ and $Q$.

###  VI.1 Bright envelope solitons

When $P/Q>0$, the expression of the bright envelope solitonic solution of Eq.
(45) can be written as Sultana2011 ; Kourakis2005 ; Schamel2002 ; Fedele2002

$\displaystyle\Phi(\xi,\tau)=\left[\psi_{0}~{}\mbox{sech}^{2}\left(\frac{\xi-U\tau}{W}\right)\right]^{1/2}$
$\displaystyle\times\exp\left[\frac{i}{2P}\left\\{U\xi+\left(\Omega_{0}-\frac{U^{2}}{2}\right)\tau\right\\}\right],$
(48)

where $\psi_{0}$ indicates the envelope amplitude, $U$ is the travelling speed
of the localized pulse, $W$ is the pulse width, which can be written as
$W=(2P\psi_{0}/Q)^{1/2}$, and $\Omega_{0}$ is the oscillating frequency for
$U=0$. The bright envelope soliton is depicted in Fig. 11. The amplitude of
the bright envelope solitons remains constant but the width of the bright
envelope solitons increases as we increase the value of the $q$-distributed
electron temperature $T_{e}$, for fixed value of ion temperature $T_{i}$ (via
$\delta$).

###  VI.2 Dark envelope solitons

As we know before that the condition for dark envelope soliton is $P/Q<0$. So,
the dark envelope soliton solution of Eq. (45) can be written as Sultana2011 ;
Kourakis2005 ; Schamel2002 ; Fedele2002

$\displaystyle\Phi(\xi,\tau)=\left[\psi_{0}~{}\mbox{tanh}^{2}\left(\frac{\xi-U\tau}{W}\right)\right]^{1/2}$
$\displaystyle\times\exp\left[\frac{i}{2P}\left\\{U\xi-\left(\frac{U^{2}}{2}-2PQ\psi_{0}\right)\tau\right\\}\right].$
(49)

The dark envelope soliton, obtained from Eq. (49) is shown in Fig. 11.

##  VII Discussion

We have studied an unmagnetized realistic space dusty plasma system consists
of non-extensive $q$-distributed electrons, isothermal ions, positively
charged warm dust particles as well as negatively charged warm dust particles.
The reductive perturbation method is used to derive the NLS equation. The
results that have been found from our investigation can be summarized as
follows

  1. 1.

The fast DA mode increases exponentially with the increasing values of $z_{2}$
for fixed value of $z_{1}$, $n_{20}$ (via $\beta$), and $n_{10}$. On the other
hand, the slow DA mode linearly increases with the increase of $z_{2}$ for the
fixed value of $z_{1}$, $n_{20}$, and $n_{10}$ (via $\beta$).

  2. 2.

The DAWs is modulationally stable (unstable) in the range of values of $k$ in
which the ratio $P/Q$ is $P/Q<0$ ($P/Q>0$).

  3. 3.

The value of $k_{c}$ increases with the increase of $n_{i0}$ for fixed value
of $z_{1}$ and $n_{10}$ (via $\mu_{i}$). On the other hand, $k_{c}$ value
decreases with the increase of $m_{2}$ for fixed value of $m_{1}$, $z_{2}$,
and $z_{1}$ (via $\alpha$).

  4. 4.

The value of $\Gamma$ increases with the decrease of $z_{2}$ for fixed value
of $z_{1}$, $m_{1}$, $m_{2}$, and $q$ (where the values of $q$ lies in the
range of $q>1$, $q=1$, and $q<1$, respectively).

  5. 5.

The maximum value of $\Gamma$ increases (decreases) with the decrease
(increase) of the values of $q>1$ ($q<1$). So, the growth rate is independent
on the sign of the $q$.

  6. 6.

The amplitude of the bright envelope solitons remains constant but the width
of the bright envelope solitons increases as we increase the value of the
$q$-distributed electron temperature $T_{e}$, for fixed value of ion
temperature $T_{i}$ (via $\delta$).

The results of our present investigation will be useful in understanding the
nonlinear phenomena both in space (e.g. Jupiters magnetosphere Horányi1996b ;
Horányi1993a ; Mamun2008c ; Chowdhury2017a , upper mesosphere Havnes1996 ,
comets tails Horányi1996b ; Mendis1991b , etc.) and laboratory (viz. direct
current and radio-frequency discharges, plasma processing reactors, fusion
plasma devices Shukla2002 , solid-fuel combustion products Shukla2002 , etc.)
plasma system containing nonextensive $q$-distributed electrons, isothermal
ions, negatively and positively charged warm dust.

## References

  * (1) P. K. Shukla and A. A. Mamun, Introdustion to Dusty Plasma Physics (Institute of Physics, Bristol, 2002).
  * (2) F. Verheest, Space Sci. Rev. 77, 267 (1996).
  * (3) C. K. Goertz, Rev. Geophys. 27, 271 (1989).
  * (4) E. C. Whipple, T. G. Northrop, and D. A. Mendis, J. Geophys. Res. 90, 7405 (1985).
  * (5) A. Mamun, R. A. Cairns, and P. K. Shukla, Phys. Plasmas 3, 702 (1996).
  * (6) N. N. Rao, P. K. Shukla, and M. Y. Yu, Planet. Space. Sci. 38, 543 (1990).
  * (7) N. D’ Angelo, Planet. Space. Sci. 38, 1143 (1990).
  * (8) A. A. Mamun and R. Schlickeiser, Phys. Plasmas 23, 034502 (2016).
  * (9) V. W. Chow, D. A. Mendis, and M. Rosenberg, J. Geophys. Res. 98, 19065 (1993).
  * (10) D. A. Mendis and M. Rosenberg, Annu. Rev. Astron. Astrophys. 32, 419 (1994).
  * (11) M. Horányi, G. E. Morfill, and E. Grün, Nature (London) 363, 144 (1993).
  * (12) D. A. Mendis and M. Horányi, in Cometary Plasma Processes, Monograph 61 (American Geophysical Union, Washington, DC, 1991).
  * (13) M. Horányi, Annu. Rev. Astron. Astrophys. 34, 383 1996.
  * (14) O. Havnes, J. Trøim, T. Blix et al., J. Geophys. Res. 101, 10839 (1996).
  * (15) V. E. Fortov et al., J. Exp. Theor. Phys. 87, 1087 (1998)
  * (16) A. A. Mamun, Phys. Rev. E 77, 026406 (2008).
  * (17) N. A. Chowdhury, A. Mannan, and A. A. Mamun, Phys. Plasmas 24, 113701 (2017).
  * (18) T. Akhter, M. M. Hossain, and A. A. Mamun, Phys. Plasmas 19, 093707 (2012).
  * (19) A. Renyi, Acta Math. Acad. Sci. Hung. 6, 285 (1955).
  * (20) C. Tsallis, J. Stat. Phys. 52, 479 (1988).
  * (21) I. Driouch and H. Chatei, Eur. Phys. J. D 9, 71 (2017).
  * (22) A. E. Dubinov, Plasma Phys. Rep. 35, 991 (2009).
  * (23) T. S. Gill, A. S. Bains, and C. Bedi, Phys. Plasmas 17, 013701 (2010).
  * (24) I. Tasnim, M. M. Masud, M. Asaduzzaman, and A. A. Mamun, Chaos 23, 013147 (2013).
  * (25) N. A. Chowdhury, M. M. Hasan, A. Mannan, and A. A. Mamun, Vacuum 147, 31-37 (2018).
  * (26) A.P. Misra and A. R. Chowdhury, Eur. Phys. J. D 39, 49–57 (2006).
  * (27) F. Sayed and A. A. Mamun, Phys. Plasmas 14, 014501 (2007).
  * (28) W. F. El-Taibany, Phys. Plasmas 20, 093701 (2013).
  * (29) N. A. Chowdhury, A. Mannan, M. M. Hasan, and A. A. Mamun, Chaos 27, 093105 (2017).
  * (30) S. Sultana and I. Kourakis, Plasma Phys. Control. Fusion 53, 045003 (2011).
  * (31) I. Kourakis and P.K. Sukla, Nonlinear Proc. Geophys. 12, 407 (2005).
  * (32) R. Fedele and H. Schamel, Eur. Phys. J. B 27, 313 (2002).
  * (33) R. Fedele, Phys. Scr. 65, 502 (2002).