# Distilling the Knowledge of Romanian BERTs Using Multiple Teachers

###### Abstract

Running large-scale pre-trained language models in computationally constrained
environments remains a challenging problem yet to be addressed, while transfer
learning from these models has become prevalent in Natural Language Processing
tasks. Several solutions, including knowledge distillation, network
quantization, or network pruning have been previously proposed; however, these
approaches focus mostly on the English language, thus widening the gap when
considering low-resource languages. In this work, we introduce three light and
fast versions of distilled BERT models for the Romanian language: Distil-BERT-
base-ro, Distil-RoBERT-base, and DistilMulti-BERT-base-ro. The first two
models resulted from the individual distillation of knowledge from two base
versions of Romanian BERTs available in literature, while the last one was
obtained by distilling their ensemble. To our knowledge, this is the first
attempt to create publicly available Romanian distilled BERT models, which
were thoroughly evaluated on five tasks: part-of-speech tagging, named entity
recognition, sentiment analysis, semantic textual similarity, and dialect
identification. Our experimental results argue that the three distilled models
offer performance comparable to their teachers, while being twice as fast on a
GPU and $\sim$35% smaller. In addition, we further test the similarity between
the predictions of our students versus their teachers by measuring their label
and probability loyalty, together with regression loyalty - a new metric
introduced in this work.

Keywords: Knowledge Distillation, BERT, Romanian language, Multiple Teachers,
Loyalty

languageresourceLanguage Resources

Distilling the Knowledge of Romanian BERTs Using Multiple Teachers

Andrei-Marius Avram1, Darius Catrina*††thanks: *Work done during an interhsip
at the Research Institute for Artificial Intelligence, Romanian Academy.3,
Dumitru-Clementin Cercel2, Mihai Dascalu2,
---
Traian Rebedea2, Vasile Păiș1, Dan Tufiș1
Research Institute for Artificial Intelligence, Romanian Academy1
University Politehnica of Bucharest, Faculty of Automatic Control and
Computers2
Duke University3
{andrei.avram, vasile<EMAIL_ADDRESS>{dumitru.cercel, mihai.dascalu,
<EMAIL_ADDRESS>

Abstract content

## 1\. Introduction

Knowledge transfer from Transformer-based language models [Vaswani et al.,
2017] trained on large amounts of data achieves state-of-the-art results on
most Natural Language Processing (NLP) tasks [Devlin et al., 2019, Liu et al.,
2019, He et al., 2020]. However, the best performing models usually have
billions [Brown et al., 2020] or even trillions [Fedus et al., 2021] of
parameters, making them impractical in certain real-world situations.
Moreover, both training and using these language models usually comes at a
high environmental cost [Strubell et al., 2019].

Several attempts were made to reduce the size of models by distilling their
knowledge [Hinton et al., 2015] accumulated during the pre-training phase
[Sanh et al., 2019], after fine-tuning the model on a specific task [Turc et
al., 2019], or both pre-training and fine-tuning [Jiao et al., 2020]. Other
methods consider shrinking the size of the models by either quantizing their
weights to integer values [Shen et al., 2020], or pruning parts of the neural
network [Brix et al., 2020]. In addition, efficient attention mechanisms were
developed to overcome the quadratic bottleneck in the sequence length of
multi-head attention [Zaheer et al., 2020, Choromanski et al., 2020].

However, the vast majority of these efforts focused on developing English
models, and little attention was paid on increasing the efficiency of pre-
trained models on other languages, with few singular cases of such compressed
models like BERTino [Muffo and Bertino, 2020] for Italian, MBERTA for Arabic
[Alyafeai and Ahmad, 2021], or
GermDistilBERT111https://huggingface.co/distilbert-base-german-cased for
German. As a response to this issue, we focus on Romanian, a language on which
BERT has recently attracted a surge of attention from the local community and
has shown promising results in various areas like dialect identification
[Zaharia et al., 2021, Popa and Ștefănescu, 2020, Zaharia et al., 2020],
document classification [Avram et al., 2021] or satire detection in news
[Rogoz et al., 2021]. Thus, our work introduces three compressed BERT versions
for the Romanian language that were obtained through a distillation process:

  * •

Distil-BERT-base-ro222https://huggingface.co/racai/distilbert-base-romanian-
cased was obtained by distilling the knowledge of BERT-base-ro [Dumitrescu et
al., 2020] using its original training corpus and tokenizer;

  * •

Distil-RoBERT-base333https://huggingface.co/racai/distilbert-base-romanian-
uncased was created from RoBERT-base [Masala et al., 2020] in similar
conditions (i.e., using both original training corpus and tokenizer);

  * •

DistilMulti-BERT-base-ro444https://huggingface.co/racai/distilbert-multi-base-
romanian-cased considered the distillation of the knowledge from an ensemble
consisting of BERT-base-ro and RoBERT-base, while relying on the combined
corpus and coupled with the tokenizer of the former model.

Our three compressed models were further evaluated on five Romanian datasets
and the results showed that they maintained most of the performance of the
original models, while being approximately twice as fast when run on a GPU. In
addition, we also measure the label, probability and regression loyalties
between each of the three distilled models and their teachers, as quantifying
the performance on specific tasks does not show how similar the predictions
between a teacher and a student really are. The models, together with the
distillation and evaluation scripts were open-sourced to improve the
reproducibility of this work555https://github.com/racai-ai/Romanian-
DistilBERT.

The rest of the paper is structured as follows. The next section presents a
series of solutions related to the knowledge distillation of pre-trained
language models. The third section outlines our approach of distilling the
knowledge of Romanian BERTs, whereas the fourth section presents the
evaluation setup and their performance on various Romanian tasks. The fifth
section evaluates the prediction loyalty between each distilled version and
its teacher, while the sixth section evaluates their inference speed. The
final section concludes our work and outlines potential future work.

## 2\. Related Work

Knowledge Distillation [Hinton et al., 2015] is a compression method in which
a smaller framework, the student model, is trained to reduce the loss
$\mathcal{L}_{KD}$ over the soft probabilities predicted by a larger model
(i.e., the teacher):

$\mathcal{L}_{KD}=\sum_{i}t_{i}\cdot log(s_{i})$ (1)

where $t_{i}$ and $s_{i}$ are the probabilities predicted by the teacher and
the student, respectively.

The technique usually uses a temperature parameter $T$ in the softmax function
that controls the smoothness of the distribution given by probabilities
$p_{i}$, known as softmax-temperature:

$p_{i}=\frac{exp(z_{i}/T)}{\sum_{j}exp(z_{j}/T)}$ (2)

where $z_{i}$ is the logit corresponding to the probability at index $i$.

?) introduced the first distilled Transformer-based language model for
English, where the authors created a model that was 1.6x smaller, 2.5x faster,
and retained 97% of the performance of the original BERT [Devlin et al., 2019]
on the General Language Understanding Evaluation (GLUE) benchmark [Wang et
al., 2018]. Soon after, TinyBERT [Jiao et al., 2020] was introduced and it
reduced the size of BERT by 7.5x, improved the inference time by 9.4x, while
maintained 96.8% of the performance of the original model on GLUE; TinyBERT
performed distillation at both the pre-training and fine-tuning stages. Also,
?) introduced MobileBERT, a slightly heavier model compared to TinyBERT - 4.3x
smaller and 5.5x faster than the original BERT, which managed to retain most
of the its knowledge by achieving a GLUE score of 77.7 (0.6 lower than BERT).

One of the first attempts to distill the knowledge of pre-trained language
models for languages other than English was BERTino [Muffo and Bertino, 2020]
for Italian. The distillation corpus was composed of 1.8 billion tokens and
bert-base-italian-xxl-uncased666https://huggingface.co/dbmdz/bert-base-
italian-xxl-uncased was used as the teacher model. The resulting model
retained most of the knowledge from the original model, as its performance was
below with values ranging from 0.29% to 5.15% on the evaluated tasks.

Additionally, a distilled version of Multilingual BERT (mBERT) - LightMBERT -
[Jiao et al., 2021] was developed to reduce the discrepancies between
languages by transferring the cross-lingual capabilities of the original model
into a smaller one. The authors first initialized the student layers with the
bottom layers of the teacher - mBERT. Then, they froze the embedding layer,
and performed distillation only on the other layers. This simple process
allowed them to maintain the cross-lingual generalization capabilities of the
original model.

Nevertheless, the knowledge transferred from a single teacher through
distillation may be limited and even biased, resulting in a student model of
low quality. As such, ?) proposed a knowledge distillation approach from
several teachers in both the pre-training and fine-tuning stages. Their
experimental results showed a significant improvement over single teacher
models like DistilBERT or TinyBERT, and even over other models used as
teachers.

## 3\. Distillation Process

Figure 1 introduces the overall architecture of the distillation process of
all three compressed Romanian BERT variants. The Distil-BERT-base-ro and
Distil-RoBERT-base models were obtained by using BERT-base-ro and RoBERT-base
as teacher, respectively, together with their teacher’s training corpora and
tokenizer. DistilMulti-BERT-base-ro used an ensemble of teachers consisting of
both Romanian BERTs - BERT-base-ro and RoBERT-base, their combined corpora and
the BERT-base-ro tokenizer.

Figure 1: Knowledge distillation process of the three distilled Romanian BERT
versions: Distil-BERT-base-ro (left), Distil-RoBERT-base (middle), and
DistilMulti-BERT-base-ro (right). DistilMulti-BERT-base-ro uses the BERT-base-
ro tokenizer, marked with dotted arrow in the figure.

### 3.1. Corpora

Corpus | Lines | Words | Size
---|---|---|---
RoWiki | 1.5M | 60.5M | 0.4 GB
OPUS | 55.1M | 635M | 3.8 GB
OSCAR | 33.6M | 1725.8M | 11 GB
Total | 90.2M | 2421.3M | 15.2 GB
Table 1: BERT-base-ro training corpus. Corpus | Lines | Words | Size
---|---|---|---
RoWiki | 2M | 50M | 0.3 GB
RoTex | 14M | 240M | 1.5 GB
OSCAR | 87M | 1780M | 10.8 GB
Total | 103M | 2070M | 12.6 GB
Table 2: RoBERT-base variants training corpus.

Tables 1 and 2 summarize the corpora used for training BERT-base-ro and
RoBERT-base, respectively. Both training corpora contain the Romanian
Wikipedia and the Open Super-large Crawled Aggregated coRpus (OSCAR) [Suárez
et al., 2019] datasets, although their final size is not the same due to
preprocessing differences. However, the BERT-base-ro corpus uses the OPUS
corpus [Tiedemann, 2012] while RoBERT-base uses the RoTex custom
collection777https://github.com/aleris/ReadME-RoTex-Corpus-Builder, resulting
in a 2.5 GB difference between the final datasets.

We observed that, in comparison with RoBERT-base training corpus, the BERT-
base-ro training corpus was a bit noisier and contained more artifacts; as
such, we applied several additional rules to ensure a cleaner version by
removing lines in the following conditions: a) lines that contained noise for
diacritics (e.g., ”c?nd” instead of ”când” - eng., ”when”); b) lines with
uncapitalized versions of the most frequent Romanian named entities (e.g.,
”bucurești” instead of ”București” - eng.: ”Bucharest”); and c) lines that
were not detected as written in Romanian by
langdetect888https://github.com/Mimino666/langdetect. Web page parsing
artifacts, such as ”Articolul Anterior” (eng., ”Previous Article”) at the
beginning of the sentence or ”Articolul Următor” (eng., ”Next Article”) at the
end of the sentence, were also removed. In contrast to the original 15.2 GB of
plain text used to train BERT-base-ro, our cleaned version of the corpus
contains 14.6 GB of text and approximately 2.3 billions words (i.e., roughly
5.9% smaller).

The combined corpus used to train DistilMulti-BERT-base-ro was obtained by
merging and deduplicating the two corpora, resulting in a corpus that
contained 25.3 GB of text and 4.1 billion words words.

### 3.2. Teacher Networks

The teacher models used to perform knowledge distillation during pre-training
were the base-cased version of BERT-base-ro [Dumitrescu et al., 2020] for
Distil-BERT-base-ro, RoBERT-base [Masala et al., 2020] for Distil-RoBERT-base,
and an ensemble of the two for DistilMulti-BERT-base-ro. Although, to the best
of our knowledge, two other Romanian BERT models exist, namely RoBERT-small
and RoBERT-large that were introduced in the same work as RoBERT-base [Masala
et al., 2020], we only chose these variants because their parameters
configuration best fit the compression process used to obtain DistilBERT.

### 3.3. Student Networks

All the distilled versions introduced in this work follow the architecture of
DistilBERT, having 6 Transformer encoder layers (i.e., half the layers of
BERT-base-ro), 12 attention heads and a hidden dimension of 768. All follow-up
analyses consider as reference also mBERT [Devlin et al., 2019], the most
accessible multi-lingual BERT-based model that supports Romanian; XLM-RoBERTa
(XLM-R) [Conneau et al., 2020] was not considered due to its size which is one
order of magnitude larger than the other models. Table 3 outlines the size,
the number of parameters, and several layers statistics (i.e., the number of
layers and heads, and the hidden size) for mBERT, RoBERT, BERT-base-ro, and
the distilled models. Both Distil-BERT-base-ro and DistilMulti-BERT-base-ro
have a size of 312 MB and contain 81 millions parameters, reducing the size of
BERT-base-ro by $\sim$35%. Distil-RoBERT-base is slightly smaller, with 72
millions parameters and a size of 282 MB, compressing RoBERT-base by the same
amount999The difference in the number of parameters between Distil-BERT-base-
ro and RoBERT-base is generated by the difference in vocabulary size..

Model | Layers | Hidden | Heads | Vocab | Size | Params
---|---|---|---|---|---|---
mBERT | 12 | 768 | 12 | 120K | 681 MB | 177M
BERT-base-ro | 12 | 768 | 12 | 50K | 477 MB | 124M
RoBERT-small | 12 | 256 | 8 | 38K | 74 MB | 19M
RoBERT-base | 12 | 768 | 12 | 38K | 441 MB | 114M
RoBERT-large | 24 | 1024 | 16 | 38K | 1.3 GB | 341M
Distil-BERT-base-ro | 6 | 768 | 12 | 50K | 312 MB | 81M
Distil-RoBERT-base | 6 | 768 | 12 | 38K | 282 MB | 72M
DistilMulti-BERT-base-ro | 6 | 768 | 12 | 50K | 312 MB | 81M
Table 3: Model size comparison of mBERT, Romanian BERTs, and our distilled
versions.

Following ?), the total loss function $\mathcal{L}$ used to pre-train Distil-
BERT-base-ro and Distil-RoBERT-base was composed of three parts
$\mathcal{L}_{KD}$, $\mathcal{L}_{MLM}$ and $\mathcal{L}_{COS}$, weighted by
different coefficients that sum to one:

$\mathcal{L}=\lambda_{KD}\mathcal{L}_{KD}+\lambda_{MLM}\mathcal{L}_{MLM}+\lambda_{COS}\mathcal{L}_{COS}$
(3)

where $\mathcal{L}_{KD}$ is the knowledge distillation loss [Hinton et al.,
2015], $\mathcal{L}_{MLM}$ is the masked language modeling (MLM) loss [Devlin
et al., 2019], $\mathcal{L}_{COS}$ is the cosine similarity embedding loss
used to align the hidden states of the teacher and student models [Sanh et
al., 2019], and $\lambda_{KD}$, $\lambda_{MLM}$ and $\lambda_{COS}$ are the
weights of each loss. A higher weight for the knowledge distillation loss
$\lambda_{KD}=0.625$ is considered because the training corpora for the
student network represents a large portion ($\sim$94%) of the pre-training
corpora for the teacher; as such, more of the internal knowledge acquired by
the teachers would be transferred to the students in the process, at the
detriment of not learning much by itself based on the MLM loss.
$\lambda_{MLM}$ and $\lambda_{COS}$ are set to 0.25 and 0.125, respectively.

The knowledge distillation loss $\mathcal{L}_{KD}$ and the cosine similarity
loss $\mathcal{L}_{COS}$ were split for training DistilMulti-BERT-base-ro into
two equally weighted parts corresponding to the losses of each model in the
ensemble :

$\mathcal{L}_{KD}=\frac{\mathcal{L}^{1}_{KD}+\mathcal{L}^{2}_{KD}}{2}$ (4)

$\mathcal{L}_{COS}=\frac{\mathcal{L}^{1}_{COS}+\mathcal{L}^{2}_{COS}}{2}$ (5)

where $\mathcal{L}^{1}_{KD}$, $\mathcal{L}^{1}_{COS}$ are the distillation and
cosine similarity losses corresponding to BERT-base-ro, and
$\mathcal{L}^{2}_{KD}$, $\mathcal{L}^{2}_{COS}$ are the distillation and
cosine similarity losses corresponding to RoBERT-base.

### 3.4. Training Settings

An important practical aspect of knowledge distillation is the initialization
of the student parameters. Following the setup of DistilBERT, we initialize
the students layers with the teacher layers by taking out the first layer of
two101010We initialized DistilMulti-BERT-base-ro with the parameters of BERT-
base-ro. (i.e. first layer, third layer, fifth layer etc.). The distilled
models were trained for 3 epochs with a batch size of 256 and a learning rate
of 5e-4. We also applied a weight decay of 1e-4 and a warm-up for the first 5%
of the entire training process, together with gradient clipping [Pascanu et
al., 2013] for gradient norms surpassing the value of 5. The temperature $T$
from Equation 2 was set to 2. The training process took approximately 30 days
on two GeForce GTX 1080 Ti for each model variant.

Task | Epochs | Batch | Warm | L-Rate
---|---|---|---|---
UPOS | 10 | 16 | 1000 | 1e-4
XPOS | 10 | 16 | 1000 | 4e-5
NER | 15 | 16 | 500 | 5e-5
SAPN | 10 | 16 | 1000 | 3e-5
SAR | 10 | 16 | 1000 | 5e-5
DI | 5 | 8 | 1500 | 5e-5
STS | - | 256 | 0 | 2e-5
Table 4: Hyperparamters used to fine-tune the models on each evaluation task.

## 4\. Task Evaluation

Five Romanian tasks were considered to create a strong evaluation setup for
our distilled models:

  * •

Part-of-Speech (POS) Tagging: Label a sequence of tokens with Universal Part-
of-Speech (UPOS) and eXtended Part-of-Speech (XPOS).

  * •

Named Entity Recognition (NER): Tag a sequence of tokens with
Inside–Outside–Beginning (IOB) labels [Ramshaw and Marcus, 1999].

  * •

Sentiment Analysis (SA): Predict whether a review expresses a positive or
negative sentiment (SAPN), together with its rating (SAR).

  * •

Dialect Identification (DI): Identify the Romanian/Moldavian dialects in news
articles.

  * •

Semantic Textual Similarity (STS): Given a pair of sentences, measure how
semantically similar they are.

The fine-tuning of the models for each task was performed by using the AdamW
optimizer [Loshchilov and Hutter, 2017] and a scheduler that linearly
decreased the learning rate to 0 at the end of the training. The number of
epochs, the batch size, the learning rate, and the warm-up steps used for each
individual task are depicted in Table 4. In line with the recommendations of
the authors of the STS evaluation
scripts111111https://github.com/dumitrescustefan/RO-STS/tree/master/baseline-
models, we did not employ any warm-up steps, a higher batch size of 256 was
considered, and early stopping was used for training instead of a maximum
number of epochs. Each experiment was run 5 times and we report the average
scores in order to mitigate the variation in performance due to the random
initialization of the weights.

### 4.1. Part-of-Speech Tagging

#### Task Description.

The original splits of the Romanian Reference Trees (RRT) corpus [Barbu
Mititelu et al., 2016] from Universal Dependencies (UD) v2.7 were used to
train and evaluate the distilled models on UPOS and XPOS tagging. The
evaluation metric employed to measure the performance on both subtasks was the
macro-averaged F1-score.

#### Methodology.

The output embeddings of the model $E_{i}$ were projected into a tensor of
dimension equal with the number of classes of UPOSes or XPOSes. The
transformation function was a feed-forward layer with weights $W$, bias $b$,
and the $LeakyReLU$ activation function that produced the logits $y_{i}$
corresponding to each embedding (i.e., $y_{i}=LeakyReLU(W^{T}E_{i}+b)$); a
dropout of 0.1 was also applied for regularization. The logits were further
transformed using the $softmax$ function to obtain the output distribution.
The cross-entropy loss between the target labels and the predicted
probabilities was employed as the value to be minimized.

#### Results

Model performance on UPOS and XPOS prediction are outlined in Table 5. The
highest F1-score of 98.07% on UPOS evaluation was achieved by DistilMulti-
BERT-base-ro, outperforming both its teachers and being the second model in
the overall leaderboard behind RoBERT-large. The best XPOS performance was
obtained by Distil-BERT-base-ro with 97.08%, surpassing BERT-base-ro with
96.46%, but falling behind its teacher and the large BERT variant.

An in-depth analysis of Romanian UPOS and XPOS evaluation was performed in
[Păiş et al., 2021], where several basic language processing kits (BLARK) were
tested on RRT. Our distilled models obtained superior results on UPOS and
comparable results with the best models on XPOS.

Model | UPOS | XPOS
---|---|---
mBERT | 97.87 | 96.16
RoBERT-small | 97.43 | 96.05
RoBERT-base | 98.02 | 97.15
RoBERT-large | 98.12 | 97.81
BERT-base-ro | 98.00 | 96.46
Distil-BERT-base-ro | 97.97 | 97.08
Distil-RoBERT-base | 97.12 | 95.79
DistilMulti-BERT-base-ro | 98.07 | 96.83
Table 5: UPOS and XPOS evaluation results on RRT.

### 4.2. Named Entity Recognition

#### Task Description.

NER evaluation was performed on Romanian Named Entity Corpus (RONEC)
[Dumitrescu and Avram, 2020]. Models were evaluated according to ?) and the
macro-averaged F1-scores are reported for the exact matches of the IOB labels.

#### Methodology.

The approach used to fine-tune the models is the same as the one used in POS
tagging, the only difference being the dimension of the output tensor that was
adjusted to match the number of classes found in RONEC.

#### Results

The results on this task are presented in Table 6. The Distil-BERT-base-ro
obtained a strict F1-score of 79.42% which is almost identical to the strict
F1-score obtained by the distilled ensemble - 79.43%, both models
outperforming Distil-RoBERT-base by approximately 0.3% on this metric. Our
compressed models lagged behind the base models by more than 3% on the strict
F1 metric (i.e., RoBERT-large and BERT-base-ro) and by more than 2.5% on the
exact F1-score (i.e., BERT-base-ro), but they managed to achieve a performance
close to the teachers on the other two metrics.

Model | Type | Partial | Strict | Exact
---|---|---|---|---
mBERT | 84.52 | 86.27 | 80.60 | 84.13
RoBERT-small | 83.11 | 84.59 | 78.61 | 82.06
RoBERT-base | 85.92 | 87.21 | 82.05 | 85.14
RoBERT-large | 86.45 | 87.19 | 82.61 | 85.09
BERT-base-ro | 86.21 | 87.84 | 82.54 | 85.88
Distil-BERT-base-ro | 83.83 | 85.73 | 79.42 | 83.35
Distil-RoBERT-base | 85.80 | 87.39 | 79.15 | 83.11
DistilMulti-BERT-base-ro | 85.48 | 87.66 | 79.43 | 83.22
Table 6: NER evaluation results on RONEC.

### 4.3. Sentiment Analysis

#### Task Description

SA was performed on the Large Romanian Sentiment Data Set (LaRoSeDa) [Tache et
al., 2021], a dataset that contains 15,000 reviews written in Romanian, of
which 7,500 are positive and 7,500 negative. We fine-tuned the models to
predict both positive and negative sentiments, as well as the rating of each
review. The models were evaluated by measuring the accuracy and the F1-macro
score for the two prediction subtasks.

#### Methodology

The model was given as input two sentences - the title of the review and the
content of the review - separated by the token [SEP]. We projected the
embedding $C$ of the token [CLS] into a scalar $y$, representing the sentiment
of the review, by using a linear layer with weights $W$ and bias $b$, to which
the $LeakyReLU$ activation function is applied (i.e., $y=\sigma(W^{T}C+b)$).
Then, the $sigmoid$ function is applied to the scalar $y$ to obtain the output
probability.

For the rating prediction, the input was also composed of the title and the
content of the review separated by the [SEP] token. However, the embedding
corresponding to the [CLS] token was projected using a feed-forward neural
network into four dimensions representing the logits of possible ratings.
Also, the $softmax$ function is employed instead of the $sigmoid$ function to
transform the output of the neural network. The models were trained to reduce
the binary cross-entropy for SAPN and the cross-entropy losses for SAR,
respectively.

#### Results

Table 7 depicts the results of the SA task. DistilRo-BERT achieved the best
performance on the positive versus negative sentiment analysis out of all
other distilled models, with 98.20% accuracy and a 98.12% F1-score. For rating
prediction, DistilRo-BERT obtained an accuracy score of 90.14% and a F1-score
of 80.51% surpassing all the other evaluated models. The difference in scores
between the sentiment and the rating prediction for each model might be due to
the noisy labeling of the ratings, as stated by the authors of the dataset
[Tache et al., 2021].

| SAPN | SAR
---|---|---
Model | Acc | F1 | Acc | F1
mBERT | 97.43 | 97.28 | 88.12 | 78.98
RoBERT-small | 97.37 | 97.22 | 87.55 | 77.81
RoBERT-base | 98.30 | 98.20 | 89.69 | 79.40
RoBERT-large | 98.25 | 98.16 | 89.91 | 79.82
BERT-base-ro | 98.07 | 97.94 | 88.45 | 79.61
Distil-BERT-base-ro | 98.20 | 98.12 | 90.14 | 80.51
Distil-RoBERT-base | 98.01 | 97.61 | 88.33 | 79.58
DistilMulti-BERT-base-ro | 98.11 | 97.74 | 89.43 | 79.77
Table 7: SA evaluation results for positive and negative and rating
classification on LaRoSeDa.

### 4.4. Dialect Identification

#### Task Description

The MOldavian and ROmanian Dialectal COrpus (MOROCO) [Butnaru and Ionescu,
2019] is a dataset that contains 33,564 samples of text that were collected
from news, annotated with their dialect and with one of six topics that each
sample belongs to. We evaluated the proposed models in terms of accuracy and
F1-score (macro averaged) on the identification of Romanian versus Moldavian
dialects on MOROCO samples.

#### Methodology

The same approach as in the previously described positive versus negative
sentiment analysis being a binary classification task was considered for this
task. The only slight difference is that the MOROCO samples do not have a
title; as such, no separation of the title and the content was necessary.

#### Results

DI results are outlined in Table 8. As it can be observed, Distil-BERT-base-ro
and DistilMulti-BERT-base-ro achieve better performance than their teachers on
this task (i.e., accuracy/F1-scores of 96.36%/96.31% and 96.26%/96.16%),
outperforming even the large variant of RoBERT.

Model | Acc | F1
---|---|---
mBERT | 95.12 | 95.06
RoBERT-small | 95.48 | 95.41
RoBERT-base | 96.10 | 96.07
RoBERT-large | 96.20 | 96.17
BERT-base-ro | 95.64 | 95.58
Distil-BERT-base-ro | 96.36 | 96.31
Distil-RoBERT-base | 96.13 | 96.11
DistilMulti-BERT-base-ro | 96.26 | 96.18
Table 8: DI evaluation results on MOROCO.

### 4.5. Semantic Textual Similarity

#### Task Description

Romanian Semantic Textual Similarity
(RoSTS)121212https://github.com/dumitrescustefan/RO-STS is a dataset that was
obtained by translating the Semantic Textual Similarity (STS)
dataset131313https://ixa2.si.ehu.eus/stswiki/index.php/STSbenchmark and used
to compare the capacity of the employed BERT models in capturing the semantic
similarity between Romanian sentences. The evaluation metrics are the Pearson
and the Spearman coefficients.

#### Methodology

The fine-tuning of the language models was performed by giving as input the
two sentences separated by the [SEP] token, followed by the projection of the
resulted [CLS] token embedding into a scalar $y$ using a linear layer with
weights $W$ and bias $b$, to which the $LeakyRELU$ function is applied (i.e.,
$y=\sigma(W^{T}C+b)$). Then, the $sigmoid$ function is applied to the output
logit. To match with the $y$ interval, the original similarities were
normalized from the $[0,5]$ to $[0,1]$. The models were trained using the mean
squared error loss.

#### Results

Task results are outlined in Table 9. DistilMulti-BERT-base-ro achieved the
highest average Pearson/Spearman scores (80.66%/80.27%) out of all distilled
models, slightly outperforming Distil-BERT-base-ro by approximately 0.1%/0.25%
and Distil-RoBERT-base by approximately 0.85%/0.25%. In comparison with the
teacher models, both Distil-BERT-base-ro and DistilMulti-BERT-base-ro achieved
better results than BERT-base-ro, but lagged behind RoBERT-base by
approximately 0.5% on both metrics and RoBERT-large by approximately 2%/1.5%.

Model | Pears | Spear
---|---|---
mBERT | 76.64 | 76.41
RoBERT-small | 78.24 | 77.84
RoBERT-base | 81.18 | 80.69
RoBERT-large | 82.51 | 81.83
BERT-base-ro | 80.30 | 79.94
Distil-BERT-base-ro | 80.57 | 80.02
Distil-RoBERT-base | 79.80 | 78.82
DistilMulti-BERT-base-ro | 80.66 | 80.27
Table 9: STS evaluation results on RoSTS.

## 5\. Loyalty

?) introduced two new metrics - loyalty and robustness - for measuring the
similarities between a student and a teacher model because just comparing the
evaluation metric for a specific task does not reflect how alike the student
and the teacher model behave. For loyalty, the authors measure similarities
between the labels (L-L) and the probabilities (P-L) predicted by the models
fine-tuned on Multi-Genre Natural Language Inference dataset [Williams et al.,
2018]. However, we choose to evaluate the loyalties of our distilled models on
the test set of LaRoSeDa as there is no natural language inference dataset for
Romanian. In addition, we introduce a new loyalty evaluation metric for
measuring the similarity for regression - regression loyalty (R-L) -, using
the test set of RoSTS as target. Robustness was measured by ?) using the
after-attack accuracy (AA); however, it is not used in this work due to the
lack of high-quality word embeddings that are retrofitted for synonymy141414It
must be noted that a version of synonymy aware word embeddings was introduced
in [Dumitrescu et al., 2018], but they were not open-sourced for public
usage.. The loyalty metrics for MultiDistil-BERT-base-ro were computed using
the average of the metrics computed with each individual teacher.

Model | L-L | P-L | R-L
---|---|---|---
Distil-BERT-base-ro | 87.80 | 74.76 | 94.05
Distil-RoBERT-base | 84.63 | 71.95 | 92.24
DistilMulti-BERT-base-ro | 89.75 | 73.23 | 94.64
Table 10: Label, probability and regression loyalty results between our
distilled models and their teachers.

### 5.1. Label Loyalty

Label Loyalty (L-L) directly measures the similarity between the labels
predicted by two models as the accuracy between the teacher labels acting as
the ground truth versus the student labels:

$LL=Accuracy(label_{t},label_{s})$ (6)

where $label_{t}$ and $label_{s}$ are the labels predicted by the teacher and
the student, respectively.

The highest L-L score was obtained by DistilMulti-BERT-base-ro with 89.75%, as
outlined in Table 10, showing that the distilled ensemble has superior label
alignment with its teachers than each individual distilled model.

### 5.2. Probability Loyalty

Computing the similarities between the output probabilities of a teacher and
student after compression is also important in industrial applications that
focus on confidence and meaningfulness [Guo et al., 2017]. P-L evaluates the
distance between the teacher output probabilities $P_{t}$ and the student
output probabilities $P_{s}$ as:

$PL=1-\sqrt{D_{JS}(P_{t}||P_{s})}$ (7)

where $D_{JS}$ is the Jensen–Shannon divergence, defined using the
Kullback–Leibler divergence over probabilities $X$ $D_{KL}(P||Q)=\sum_{x\in
X}P(x)log(\frac{P(x)}{Q(x)})$ as:

$D_{JS}(P_{t}||P_{s})=\frac{D_{KL}(P_{t}||P_{s})+D_{KL}(P_{s}||P_{t})}{2}$ (8)

The results of P-L evaluation are depicted in the second column of Table 10.
The highest score was obtained by Distil-BERT-base-ro with 74.76%,
outperforming the distilled ensemble model with approximately 1.5%.

### 5.3. Regression Loyalty

The previous two metrics do not assess performance in the case of regression
problems. Thus, this work introduces regression loyalty as the Pearson
Correlation Coefficient [Breese et al., 1998] between the output of the
teacher $pred_{t}$ and the output of the student $pred_{s}$:

$RL=Pearson(pred_{t},pred_{s})$ (9)

The results for R-L are outlined in the last column of Table 10. DistilMulti-
BERT-base-ro achieved the highest regression loyalty with its two teachers (a
94.64% R-L score), followed by Distil-BERT-base-ro (94.05%) and Distil-RoBERT-
base (92.24%).

## 6\. Inference Speed

The inference time of the three distilled versions in comparison with the
other Romanian or multilingual models was evaluated on random sequences with
lengths ranging from 16 tokens to 512 tokens, using both a CPU - Intel
i7-7700K - and a GPU - GeForce GTX 1080 Ti. The models were grouped by their
size into four categories to make the visualization more appealing, namely:
Multilingual (i.e., mBERT), Small (i.e., RoBERT-small), Base (i.e., BERT-base-
ro and RoBERT-base), Large (i.e., RoBERT-large), and Distilled (i.e., Distill-
BERT-base-ro, Distill-RoBERT-base, and DistillMulti-BERT-base-ro). The
inference times are depicted in Figure 2.

The distilled models obtained a significant improvement on the GPU, being
almost twice as fast as the base and small models, and at least three times
faster than the RoBERT-large. It should be noted that the small variant
obtains a similar speed on the GPU for larger sequence lengths due to more
parallelization happening at the level of the input tokens; RoBERT-small
benefits from its reduced size and overcomes its disadvantage in
parallelization by having 12 layers instead of 6. The inference times on the
CPU of the distilled models is also better than those obtained by the base and
large models, but worse than those of RoBERT-small because it has less
parameters.

Figure 2: Inference time on CPU (up) and on GPU (down) of the evaluated
models, grouped by their dimensions: Multilingual (mBERT), Small (RoBERT-
small), Base (BERT-base-ro, RoBERT-base), Large (RoBERT-large) and Distilled
(Distill-BERT-base-ro, Distill-RoBERT-base, DistillMulti-BERT-base-ro).

## 7\. Conclusions

Creating models for low resource languages is an important area of research,
aiming to empower the usage of NLP in various contexts. This work introduces
three distilled models from Romanian BERT models: Distil-BERT-base-ro, Distil-
RoBERT-base, and DistilMulti-BERT-base-ro, that are 35% smaller than their
teacher models (BERT-base-ro and RoBERT-base). Each model was evaluated on
five Romanian NLP tasks: part-of-speech tagging, named entity recognition,
sentiment analysis, dialect identification, and semantic textual similarity.
Our experimental results showed that the distilled models maintain most of the
prediction performance of their original models, while containing only half of
their layers. Moreover, the distilled models even outperformed their teachers
and other evaluated models on several tasks. We further tested the inference
speed of the three distilled models on various sequence lengths, and our
results outlined a reduction of almost 50% on a single GPU in comparison to
their teacher models.

Future work considers the creation of a distilled version using RoBERT-large
as teacher, as well as further reducing the size of the current models by
pruning and/or quantizing their weights. In addition, we intend to fine-tune
the distilled models and annotate the Romanian sub-corpus from the CURLICAT
multilingual corpus covering domains relevant for CEF Digital Service Infra-
structures (DSIs). Besides standard mark-up (lemmatization, POS tagging,
chunking, dependency parsing) annotation will include recognition and labeling
of the occurrences of terms recorded in the Interactive Terminology for Europe
(IATE) 151515https://iate.europa.eu/home, s described in the Curated
Multilingual Language Resources for CEF.AT (CURLICAT)
161616https://curlicat.eu/ project. Finally, the results presented in this
paper are going to be included on the five evaluation tasks in LiRo
[Dumitrescu et al., 2021], a benchmark for Romanian language NLP tasks.

## Acknowledgements

This research was supported by the EC grant INEA/CEF/ICT/A2018/28592472 for
the Action No: 2019-EU-IA-0034 entitled “Curated Multilingual Language
Resources for CEF.AT” (CURLICAT) and by the Romanian Ministry of European
Investments and Projects through the Competitiveness Operational Program (POC)
project ”HOLOTRAIN” (grant no. 29/221, SMIS code: 129077).

## 8\. Bibliographical References

## References

  * Alyafeai and Ahmad, 2021 Alyafeai, Z. and Ahmad, I.  (2021).  Arabic compact language modelling for resource limited devices.  In Proceedings of the Sixth Arabic Natural Language Processing Workshop, pages 53–59.
  * Avram et al., 2021 Avram, A.-M., Păiș, V., and Tufis, D. I.  (2021).  Pyeurovoc: A tool for multilingual legal document classification with eurovoc descriptors.  In Proceedings of the International Conference on Recent Advances in Natural Language Processing (RANLP 2021), pages 92–101.
  * Barbu Mititelu et al., 2016 Barbu Mititelu, V., Ion, R., Simionescu, R., Irimia, E., and Perez, C.-A.  (2016).  The romanian treebank annotated according to universal dependencies.  In Proceedings of the tenth international conference on natural language processing (hrtal2016).
  * Breese et al., 1998 Breese, J. S., Heckerman, D., and Kadie, C.  (1998).  Empirical analysis of predictive algorithms for collaborative filtering.  In Proceedings of the Fourteenth conference on Uncertainty in artificial intelligence, pages 43–52.
  * Brix et al., 2020 Brix, C., Bahar, P., and Ney, H.  (2020).  Successfully applying the stabilized lottery ticket hypothesis to the transformer architecture.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 3909–3915.
  * Brown et al., 2020 Brown, T. B., Mann, B., Ryder, N., Subbiah, M., Kaplan, J., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., et al.  (2020).  Language models are few-shot learners.  arXiv preprint arXiv:2005.14165.
  * Butnaru and Ionescu, 2019 Butnaru, A. and Ionescu, R. T.  (2019).  Moroco: The moldavian and romanian dialectal corpus.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 688–698.
  * Choromanski et al., 2020 Choromanski, K., Likhosherstov, V., Dohan, D., Song, X., Gane, A., Sarlos, T., Hawkins, P., Davis, J., Mohiuddin, A., Kaiser, L., et al.  (2020).  Rethinking attention with performers.  arXiv preprint arXiv:2009.14794.
  * Conneau et al., 2020 Conneau, A., Khandelwal, K., Goyal, N., Chaudhary, V., Wenzek, G., Guzmán, F., Grave, É., Ott, M., Zettlemoyer, L., and Stoyanov, V.  (2020).  Unsupervised cross-lingual representation learning at scale.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 8440–8451.
  * Devlin et al., 2019 Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  (2019).  Bert: Pre-training of deep bidirectional transformers for language understanding.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186.
  * Dumitrescu and Avram, 2020 Dumitrescu, S. D. and Avram, A.-M.  (2020).  Introducing ronec-the romanian named entity corpus.  In Proceedings of The 12th Language Resources and Evaluation Conference, pages 4436–4443.
  * Dumitrescu et al., 2018 Dumitrescu, S. D., Avram, A. M., Morogan, L., and Toma, S.-A.  (2018).  Rowordnet–a python api for the romanian wordnet.  In 2018 10th International Conference on Electronics, Computers and Artificial Intelligence (ECAI), pages 1–6. IEEE.
  * Dumitrescu et al., 2020 Dumitrescu, S., Avram, A.-M., and Pyysalo, S.  (2020).  The birth of romanian bert.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: Findings, pages 4324–4328.
  * Dumitrescu et al., 2021 Dumitrescu, S. D., Rebeja, P., Lorincz, B., Gaman, M., Avram, A., Ilie, M., Pruteanu, A., Stan, A., Rosia, L., Iacobescu, C., Morogan, L., Dima, G., Marchidan, G., Rebedea, T., Chitez, M., Yogatama, D., Ruder, S., Ionescu, R. T., Pascanu, R., and Patraucean, V.  (2021).  Liro: Benchmark and leaderboard for romanian language tasks.  In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 1).
  * Fedus et al., 2021 Fedus, W., Zoph, B., and Shazeer, N.  (2021).  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity.  arXiv preprint arXiv:2101.03961.
  * Guo et al., 2017 Guo, C., Pleiss, G., Sun, Y., and Weinberger, K. Q.  (2017).  On calibration of modern neural networks.  In International Conference on Machine Learning, pages 1321–1330. PMLR.
  * He et al., 2020 He, P., Liu, X., Gao, J., and Chen, W.  (2020).  Deberta: Decoding-enhanced bert with disentangled attention.  arXiv preprint arXiv:2006.03654.
  * Hinton et al., 2015 Hinton, G., Vinyals, O., and Dean, J.  (2015).  Distilling the knowledge in a neural network.  In NIPS Deep Learning and Representation Learning Workshop.
  * Jiao et al., 2020 Jiao, X., Yin, Y., Shang, L., Jiang, X., Chen, X., Li, L., Wang, F., and Liu, Q.  (2020).  Tinybert: Distilling bert for natural language understanding.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: Findings, pages 4163–4174.
  * Jiao et al., 2021 Jiao, X., Yin, Y., Shang, L., Jiang, X., Chen, X., Li, L., Wang, F., and Liu, Q.  (2021).  Lightmbert: A simple yet effective method for multilingual bert distillation.  arXiv preprint arXiv:2103.06418.
  * Jin et al., 2020 Jin, D., Jin, Z., Zhou, J. T., and Szolovits, P.  (2020).  Is bert really robust? a strong baseline for natural language attack on text classification and entailment.  In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 8018–8025.
  * Liu et al., 2019 Liu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D., Levy, O., Lewis, M., Zettlemoyer, L., and Stoyanov, V.  (2019).  Roberta: A robustly optimized bert pretraining approach.  arXiv preprint arXiv:1907.11692.
  * Loshchilov and Hutter, 2017 Loshchilov, I. and Hutter, F.  (2017).  Decoupled weight decay regularization.  arXiv preprint arXiv:1711.05101.
  * Masala et al., 2020 Masala, M., Ruseti, S., and Dascalu, M.  (2020).  Robert–a romanian bert model.  In Proceedings of the 28th International Conference on Computational Linguistics, pages 6626–6637.
  * Muffo and Bertino, 2020 Muffo, M. and Bertino, E.  (2020).  Bertino: An italian distilbert model.  In CLiC-it.
  * Pascanu et al., 2013 Pascanu, R., Mikolov, T., and Bengio, Y.  (2013).  On the difficulty of training recurrent neural networks.  In International conference on machine learning, pages 1310–1318. PMLR.
  * Popa and Ștefănescu, 2020 Popa, C. and Ștefănescu, V.  (2020).  Applying multilingual and monolingual transformer-based models for dialect identification.  In Proceedings of the 7th Workshop on NLP for Similar Languages, Varieties and Dialects, pages 193–201.
  * Păiş et al., 2021 Păiş, V., Ion, R., Avram, A.-M., Mitrofan, M., and Tufiș, D.  (2021).  In-depth evaluation of Romanian natural language processing pipelines.  Romanian Journal of Information Science and Technology (ROMJIST), 24(4):384–401.
  * Ramshaw and Marcus, 1999 Ramshaw, L. A. and Marcus, M. P.  (1999).  Text chunking using transformation-based learning.  In Natural language processing using very large corpora, pages 157–176. Springer.
  * Rogoz et al., 2021 Rogoz, A.-C., Gaman, M., and Ionescu, R. T.  (2021).  Saroco: Detecting satire in a novel romanian corpus of news articles.  arXiv preprint arXiv:2105.06456.
  * Sanh et al., 2019 Sanh, V., Debut, L., Chaumond, J., and Wolf, T.  (2019).  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.  arXiv preprint arXiv:1910.01108.
  * Segura-Bedmar et al., 2013 Segura-Bedmar, I., Martínez, P., and Herrero-Zazo, M.  (2013).  Semeval-2013 task 9: Extraction of drug-drug interactions from biomedical texts (ddiextraction 2013).  In Second Joint Conference on Lexical and Computational Semantics (*SEM), Volume 2: Proceedings of the Seventh International Workshop on Semantic Evaluation (SemEval 2013), pages 341–350.
  * Shen et al., 2020 Shen, S., Dong, Z., Ye, J., Ma, L., Yao, Z., Gholami, A., Mahoney, M. W., and Keutzer, K.  (2020).  Q-bert: Hessian based ultra low precision quantization of bert.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 8815–8821.
  * Strubell et al., 2019 Strubell, E., Ganesh, A., and McCallum, A.  (2019).  Energy and policy considerations for deep learning in nlp.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 3645–3650.
  * Suárez et al., 2019 Suárez, P. J. O., Sagot, B., and Romary, L.  (2019).  Asynchronous pipeline for processing huge corpora on medium to low resource infrastructures.  In 7th Workshop on the Challenges in the Management of Large Corpora (CMLC-7). Leibniz-Institut für Deutsche Sprache.
  * Sun et al., 2020 Sun, Z., Yu, H., Song, X., Liu, R., Yang, Y., and Zhou, D.  (2020).  Mobilebert: a compact task-agnostic bert for resource-limited devices.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 2158–2170.
  * Tache et al., 2021 Tache, A., Mihaela, G., and Ionescu, R. T.  (2021).  Clustering word embeddings with self-organizing maps. application on laroseda–a large romanian sentiment data set.  In Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume, pages 949–956.
  * Tiedemann, 2012 Tiedemann, J.  (2012).  Parallel data, tools and interfaces in opus.  In Eighth International Conference on Language Resources and Evaluation (LREC’12), pages 2214–2218.
  * Turc et al., 2019 Turc, I., Chang, M.-W., Lee, K., and Toutanova, K.  (2019).  Well-read students learn better: On the importance of pre-training compact models.  arXiv preprint arXiv:1908.08962.
  * Vaswani et al., 2017 Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  (2017).  Attention is all you need.  In Advances in neural information processing systems, pages 5998–6008.
  * Wang et al., 2018 Wang, A., Singh, A., Michael, J., Hill, F., Levy, O., and Bowman, S.  (2018).  Glue: A multi-task benchmark and analysis platform for natural language understanding.  In Proceedings of the 2018 EMNLP Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP, pages 353–355.
  * Williams et al., 2018 Williams, A., Nangia, N., and Bowman, S.  (2018).  A broad-coverage challenge corpus for sentence understanding through inference.  In Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers), pages 1112–1122.
  * Wu et al., 2021 Wu, C., Wu, F., and Huang, Y.  (2021).  One teacher is enough? pre-trained language model distillation from multiple teachers.  arXiv preprint arXiv:2106.01023.
  * Xu et al., 2021 Xu, C., Zhou, W., Ge, T., Xu, K., McAuley, J., and Wei, F.  (2021).  Beyond preserved accuracy: Evaluating loyalty and robustness of bert compression.  arXiv preprint arXiv:2109.03228.
  * Zaharia et al., 2020 Zaharia, G.-E., Avram, A.-M., Cercel, D.-C., and Rebedea, T.  (2020).  Exploring the power of romanian bert for dialect identification.  In Proceedings of the 7th Workshop on NLP for Similar Languages, Varieties and Dialects, pages 232–241.
  * Zaharia et al., 2021 Zaharia, G.-E., Avram, A.-M., Cercel, D.-C., and Rebedea, T.  (2021).  Dialect identification through adversarial learning and knowledge distillation on romanian bert.  In Proceedings of the Eighth Workshop on NLP for Similar Languages, Varieties and Dialects, pages 113–119.
  * Zaheer et al., 2020 Zaheer, M., Guruganesh, G., Dubey, K. A., Ainslie, J., Alberti, C., Ontanon, S., Pham, P., Ravula, A., Wang, Q., Yang, L., et al.  (2020).  Big bird: Transformers for longer sequences.  In NeurIPS.


# Connected and disconnected stable regions of solitons of nonlinear
Schrödinger equation with $\mathcal{PT}$-symmetric potential

Niladri Ghosh  Email<EMAIL_ADDRESS>Department of Mathematics,
University of Kalyani, Kalyani - 741235, WB, India.  Amiya Das  Email:
<EMAIL_ADDRESS>Department of Mathematics, University of Kalyani, Kalyani
- 741235, WB, India.  Debraj Nath  Corresponding author. Email:
<EMAIL_ADDRESS>Department of Mathematics, Vivekananda College, Kolkata -
700063, WB, India.

###### Abstract

We have considered cubic nonlinear Schrödinger equation along with
supersymmetric $\mathcal{PT}$ like potential and obtained exact stationary
solutions in terms of bright and brigh-dark interacting solitons. The
$\mathcal{PT}$ broken and $\mathcal{PT}$ unbroken regions are demonstrated
also depicted. Connected and disconnected stable regions of bright and dark
solitons are examined incorporating linear stability analysis validated by
direct numerical simulations. Moreover, the strength of stability has been
illustrated through excitations of bright and dark solitons.

##  I Introduction

A fundamental axiom of quantum mechanics is that every physical observable is
associated with a real spectrum and thus operators representing such physical
observables must be Hermitian. In quantum mechanics, in order for the energy
levels to be real and the theory to be probability conserving (or unitary
evolution), it is usually assumed that the Hamiltonian (Schrödinger) operator
be Hermitian. In the case of the Hamiltonian operator, this requirement gives
rise to real energy levels and the theory guaranteing conservation of
probability. In recent years, in a series of papers by Bender and coworkers
[1], a considerable attention has been shown in a weaker version of the
Hermiticity axiom in which many non-Hermitian Hamiltonians may also lead to
entirely real spectra provided they possess something known as $\mathcal{PT}$
(parity-time) symmetry. The fundamental idea of $\mathcal{PT}$ symmetric
quantum mechanics is to replace the concept of Hermitian Hamiltonian by weaker
condition that it possess space-time reflection symmetry
($\mathcal{PT}$-symmetry). The linear space-reflection operator $\mathcal{P}$,
responsible for spatial reflection, is defined through the operations
$p\to-p\,,x\to-x$, while the anti-linear time-reversal operator $\mathcal{T}$
leads to $p\to-p\,,x\to x$ and to complex conjugation
$\textrm{i}\to-\textrm{i}$. Dynamical systems are known to be $\mathcal{PT}$
symmetric if they remain invariant under the combined parity $\mathcal{P}$ and
time-reversal $\mathcal{T}$ transformation. A necessary condition for a
Hamiltonian to be $\mathcal{PT}$ symmetric is that the potential function
$V_{PT}=V(x)+\textrm{i}W(x)$ should satisfy the condition
$V_{PT}(x)=V^{*}_{PT}(-x)$, with $*$ denoting the complex conjugation. A
$\mathcal{PT}$ symmetric system exhibits entirely real spectra as long as the
robustness of the imaginary component of the complex potential is less than a
certain threshold value. The corresponding system undergoes a drastic phase
transition also known as spontaneous $\mathcal{PT}$ symmetry breaking [2] on
exceeding the threshold value. The $\mathcal{PT}$ symmetric wave propagation
in optical structures can be replaced by the complex refractive index and
balanced gain-loss profile which was judiciously predicted through the
experimental execution in optical metamaterials, synthetic photonic lattices,
microring lasers and optically induced atomic lattices etc [3]. Though the
theory was initially originated in quantum mechanics, but later on the concept
of $\mathcal{PT}$ symmetry has also found applications in different fields,
e.g., optics [4, 5, 6, 7, 8], electronics [9], Bose-Einstein condensation [10,
11, 12], metamaterials [13, 14, 15, 16, 17], etc. Furthermore $\mathcal{PT}$
symmetry has been realized experimentally [18, 19, 20, 21].

In the last decades, the theory of existence of nonlinear localized modes in
$\mathcal{PT}$ symmetric potential and their linear stability has been studied
very promptly [22, 23]. Several interesting potentials for example Scarf-II
potential [24], harmonic potential [25], Rosen-Morse potential [26], Gaussian
potential [27, 28, 29, 30, 31, 32], sextic anharmonic double-well potential
[33], time-dependent harmonic-Gaussian potential [34] etc are considered for
study.

The nonlinear Schrödinger equation (NLSE) is one of the pioneer mathematical
model which arises in diversified physical systems. This equation is a key
model describing wave processes in plasma physics, Bose-Einstein condensates
(BEC), gravitational models for quantum mechanics, wave propagation in
biological and geological systems [35] and nonlinear optics [36]. In some of
these fields and many others, the NLS equation appears as an asymptotic limit
for a slowly varying dispersive wave envelope propagating in a nonlinear
medium. It also plays an important role in describing full spatiotemporal
optical solitons or light bullets in the theory of nonlinear optics [36].
During the past few years, exact solutions of NLSE with $\mathcal{PT}$
symmetric potentials and their stability have been studied by many authors
since they are useful in different contexts [26, 37, 38, 39, 40, 41, 42, 43,
44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56]. In Ref.[54], the exact
stationary solutions of derivative nonlinear Schrödinger equation have been
obtained in presence of a $\mathcal{PT}$ symmetric potential as a sum of
super-Gaussian and parabolic potentials. The theory of supersymmetry and
deformed supersymmetry are used in [55, 56] to retrieve exact analytical
localized solutions in context of a number of complex $\mathcal{PT}$ symmetric
potentials and supersymmetric potentials along with power law nonlinearity.

As we know that, a pair of potentials $V_{\pm}(x)$ are said to be
supersymmetric, if there exists a function $U(x)$ such that,
$V_{\pm}(x)=U^{2}(x)\pm U^{\prime}(x)$ [57]. In this paper, we have considered
$(a,v_{0},w_{0},U(x))$, where $U(x)=U_{R}(x)+\textrm{i}U_{I}$ such that,

$\begin{array}[]{l}V=v_{0}\,Re(U^{2}-U^{\prime})=v_{0}\left(U_{R}^{2}-U_{I}^{2}-U_{R}^{\prime}\right),\\\
W=w_{0}\,Im\left(U^{2}-U^{\prime}\right)=w_{0}\,\left(2U_{R}U_{I}-U_{I}^{\prime}\right),\end{array}$
(1)

$v_{0}$ and $w_{0}$ are amplitudes of real and imaginary parts of the
potential, with $w_{0}\neq 0$ and $a$ is a real constant, which will be
conneted to solutions of a NLSE. Then the functions $U_{R}$ and $U_{I}$
satisfy $\mathcal{PT}$-symmetric conditions

$U_{R}(x)=-U_{R}(-x),~{}~{}~{}~{}U_{I}(x)=U_{I}(-x),$ (2)

and a cubic nonlinear Schrödinger equation has exact solutions for deformed
supersymmetric potential $V(x)+iW(x)$. The optical solitons are classified
into two forms based on the sign of group velocity dispersion. Solitons in
negative dispersion regime are known as bright solitons (BS) and solitons in
positive dispersion regime are called dark solitons (DS). Dark soliton can be
stressed as a localized pulse appearing as an acute dip on the background of
continuous wave. On the other hand, bright soliton appears as a crest above
the continuous wave. Recently, in a series of paper [58], the existence of
dark solitons are verified experimentally in the area of nonlinear fiber
optics, plasma, BEC and waveguide arrays. In most of the cases, cubic NLSE has
either a bright or dark soliton. In this paper, we will find two solutions of
a cubic NLSE and they will be the ground and first excited states of linear
Schrödinger equiation. In this paper, we will investigate stability analysis
of BS and DS of a cubic nonlinear Schrödinger equation with deformed
supersymmetric potential. Then we will find stable regions of BS and DS and
regions will be verified by different numerical methods. Finally, regions will
be checked by excitations of BS and DS of a cubic time dependent NLSE.

The rest of the paper is structured as follows. In Sec. II, we will obtain
exact stationary solutions in terms of BS and bright-dark soliton interaction
of a cubic nonlinear Schrödinger equation along with complex $\mathcal{PT}$
symmetric potential. In Sec. III, linear stability analysis is illustrated
incorporating direct numerical simulations, which assures about stable
behavior of BS and DS depending on potential parameters. Also, excitations of
nonlinear BS and DS of a time dependent NLSE and some numerical results have
been reported in Sec. IV. Finally, conclusions have been drawn in Sec. V.

##  II Methodology

Let us consider a cubic NLSE

$i\Psi_{t}=-\Psi_{xx}+\left(V(x)+\textrm{i}W(x)\right)\Psi+g|\Psi|^{2}\Psi$
(3)

with a complex $\mathcal{PT}$-symmetric potential $V(x)+\textrm{i}W(x)$, where
$g$ is the nonlinearity parameter. To solve Eq.(3) a transformation

$\Psi(x,t)=e^{-\textrm{i}\mu t}\psi(x)$ (4)

is considered, where $\mu$ is a real propagation constant. Then from Eq.(3)
one can obtain

$-\frac{d^{2}\psi}{dx^{2}}+(V+\textrm{i}W)\psi+g|\psi|^{2}\psi=\mu\psi$ (5)

and it has a solution of the form
$C_{0}\,e^{-\int(U_{R}+\textrm{i}\,a\,U_{I})dx}$, $a\in\mathbb{R}$. In
particular if, $(a,v_{0},w_{0},U)=(1,1,1,U)$ and $g=0$, then the NLSE reduces
to a LSE with exact supersymmetric potentials with zero energy normalized
solution $e^{-\int(U_{R}+\textrm{i}\,U_{I})dx}$. In this paper we will find
solutions of NLSE (5) as well as (3), where $(a,v_{0},w_{0})\neq(1,1,1)$.

###  II.1 Solution 1

A solution of Eq.(5) is taken to be of the form

$\psi_{0}=C_{0}\,e^{-\int(U_{R}+\textrm{i}\,a\,U_{I})dx},~{}~{}\mu=\mu_{0}$
(6)

where $C_{0}$ is the amplitude and
$\lim\limits_{|x|\rightarrow\infty}|\psi_{0}(x)|=0$. Now substituting Eq.(6)
into Eq.(5) and equating real and imaginary parts, we obtain [55]

$\begin{array}[]{l}\left(v_{0}-1\right)\left(U_{R}^{2}-U_{R}^{{}^{\prime}}\right)+\left(a^{2}-v_{0}\right)U_{I}^{2}\\\
+gC_{0}^{2}\,e^{-2\int U_{R}(x)\,dx}=\mu_{0},\\\
\left(2U_{R}U_{I}-U_{I}^{{}^{\prime}}\right)\left(w_{0}-a\right)=0.\end{array}$
(7)

If $2U_{R}U_{I}-U^{\prime}_{I}=0$, then $\mathcal{PT}$-symmetric potential
becomes real potential. In this paper we have considered $w_{0}=a$ and
therefore $(a,v_{0},w_{0},U)$ exists which satisfy relations of (7) and
solution (6) esists.

###  II.2 Solution 2

If possible let, there be another solution of the form

$\begin{array}[]{l}\psi_{1}(x)=C_{1}\,\displaystyle
e^{-\displaystyle\int\left(U_{R}+\textrm{i}\,a\,U_{I}\right)\,dx}F(x),\mu=\mu_{1},\end{array}$
(8)

where $\lim\limits_{|x|\rightarrow\infty}|\psi_{1}(x)|\neq 0$. Using Eq.(8)
into Eq.(5), we obtain

$\begin{array}[]{ll}-F^{\prime\prime}+2\left(U_{R}+\textrm{i}\,a\,U_{I}\right)F^{\prime}\\\
+g\left\\{C_{1}^{2}|F|^{2}-C_{0}^{2}\right\\}\,e^{-2\displaystyle\int
U_{R}\,dx}F=\left(\mu_{1}-\mu_{0}\right)F.\end{array}$ (9)

Eq.(9) is the additional condition for existence of solution of Eq.(8). In
particular, $g=0$, then conditions are

$\begin{array}[]{l}\left(v_{0}-1\right)\left(U_{R}^{2}-U_{R}^{{}^{\prime}}\right)+\left(a^{2}-v_{0}\right)U_{I}^{2}=\mu_{0},\\\
\left(2U_{R}U_{I}-U_{I}^{{}^{\prime}}\right)\left(w_{0}-a\right)=0,\\\
-F^{\prime\prime}+2\left(U_{R}+i\,a\,U_{I}\right)F^{\prime}=\left(\mu_{1}-\mu_{0}\right)F.\end{array}$
(10)

Eq.(10) has the solution $(a,v_{0},w_{0})=(1,1,1)$ and $F$ satisfies the
relation

$-F^{\prime\prime}(x)+2U(x)F^{\prime}(x)=\mu_{1}\,F(x),$ (11)

which has infinitely many solutions depends on $U(x)$. In this paper, we will
find solution of (9) by an example. Then the NLSE (3) has two solutions of the
forms (6) and (8). The solution (13) is a BS and the solution (14) may be BS
or DS.

###  II.3 Example

Let us choose $U_{R}(x)=\tanh x,~{}~{}U_{I}(x)=\beta\rm{sech}\,x$, where
$\beta$ is a real constant. Then $\mathcal{PT}$-symmetric potential
$V(x)+iW(x)$ is defined by

$\begin{array}[]{ll}V(x)=v_{0}-v_{0}(2+\beta^{2})\,sech^{2}x,\\\
W(x)=3w_{0}\beta\,sech\,x\,\tanh x.\end{array}$ (12)

This potential is known as $\mathcal{PT}$-symmetric Scarf potential [57]. In
Fig. 1, we have plotted real and imaginary parts of the potential, for
different values of parameters and the parameters are taken from
$\mathcal{PT}$ broken and unbroken region of the linear Schrödinger equation
with this complex potential. In Sec. II.4, the $\mathcal{PT}$ broken and
unbroken symmetry are discussed. It is to be noted that, this potential was
considered within the context of cubic NLSE [4, 37, 38, 40, 26, 46, 50, 52,
59, 55, 56].

Figure 1:  Plot of real and imaginary parts of $\mathcal{PT}$-symmetric
potential (12). Parameters are (A) black lines
$v_{0}=0.133,\beta=1,w_{0}=-0.218$ and blue lines
$v_{0}=0.063,\beta=1,w_{0}=0.148$; (B) black lines
$v_{0}=0.2,\beta=1,w_{0}=0.11$ and blue lines $v_{0}=0.4,\beta=1,w_{0}=-0.08$;
(C) black lines $v_{0}=0.5,\beta=0.4,w_{0}=-4$ and blue lines
$v_{0}=0.5,\beta=1.3,w_{0}=2$; (D) black lines
$v_{0}=0.5,\beta=-0.5,w_{0}=0.2$ and blue lines
$v_{0}=0.5,\beta=-1.4,w_{0}=-0.016$.

Then solution 1 is given by

$\Psi_{0}(x,t)=C_{0}\,\rm{sech}\,x\,e^{-\textrm{i}w_{0}\beta\tan^{-1}\left(\sinh
x\right)-\textrm{i}\mu_{0}t},$ (13)

and solution 2 is given by

$\Psi_{1}(x,t)=C_{1}\left(\gamma\rm{sech}\,x+\textrm{i}\delta\tanh
x\right)e^{-\textrm{i}\,w_{0}\beta\tan^{-1}\left(\sinh
x\right)-\textrm{i}\mu_{1}t},$ (14)

where

$~{}~{}\mu_{0}=v_{0}-1,~{}~{}C_{0}=\left\\{\frac{v_{0}(2+\beta^{2})-w_{0}^{2}\beta^{2}-2}{g}\right\\}^{\frac{1}{2}},$
(15)

and

$\begin{array}[]{ll}\mu_{1}=v_{0}+\displaystyle\frac{v_{0}(2+\beta^{2})-w_{0}^{2}\beta^{2}-2}{4\beta^{2}w_{0}^{2}-1},\\\
C_{1}=\displaystyle\left\\{\frac{v_{0}(2+\beta^{2})-w_{0}^{2}\beta^{2}-2}{g\delta^{2}(4\beta^{2}w_{0}^{2}-1)}\right\\}^{\frac{1}{2}},\\\
\gamma=-2\beta\delta w_{0}.\end{array}$ (16)

We have plotted real and imaginary parts of solution (13) in Fig. 2 and of
solution (14) in Fig. 3 for different set of parameters. It is to be noted
that, for Fig. 2 and Fig. 3, all parameters are taken from $\mathcal{PT}$
broken stable, $\mathcal{PT}$ broken unstable and $\mathcal{PT}$ unbroken
stable, $\mathcal{PT}$ unbroken unstable regions of Fig. 6, 7 and 8. We
observe that, solution (13) is a BS and solution (14) is a BS, if
$w_{0}^{2}>\frac{1}{4\beta^{2}}$ and a DS, if
$w_{0}^{2}<\frac{1}{4\beta^{2}}$. Therefore, the phase transition between BS
and DS is occured for the solution (14) at the points on the curve
$w_{0}^{2}=\frac{1}{4\beta^{2}}$. Now, we will focus on LSE with complex
potential (12).

Figure 2:  Plot of real and imaginary parts of (13). The parameters are taken
from (A) stable $\mathcal{PT}$ broken region: for black lines
$v_{0}=0.133,w_{0}=-0.218$ and for blue lines $v_{0}=0.063,w_{0}=0.148$; (B)
stable $\mathcal{PT}$ unbroken region: for black lines $v_{0}=0.2,w_{0}=0.11$
and for blue lines $v_{0}=0.4,w_{0}=-0.08$; (C) unstable $\mathcal{PT}$ broken
region: for black lines $\beta=0.4,w_{0}=-4.0$ and for blue lines
$\beta=1.3,w_{0}=2.0$, (D) stable $\mathcal{PT}$ unbroken region: for black
lines $\beta=-0.5,w_{0}=0.2$ and for blue lines $\beta=-1.4,w_{0}=-0.016$. All
regions are shown in Fig. 6. Figure 3:  Plot of real and imaginary parts of
(14). The parameters are taken from (A) unstable $\mathcal{PT}$ broken region:
(black lines) $v_{0}=-0.2,w_{0}=0.40$ and (blue lines)
$v_{0}=0.4,w_{0}=-0.45$; (B) stable $\mathcal{PT}$ unbroken region: (black
lines) $v_{0}=0.20,w_{0}=-0.16$ and (blue lines) $v_{0}=0.65,w_{0}=0.27$; (C)
unstable $\mathcal{PT}$ broken region: (black lines) $\beta=0.23,w_{0}=-2.00$
and (blue lines) $\beta=0.29,w_{0}=1.70$, (D) stable $\mathcal{PT}$ unbroken
region: (black lines) $\beta=-0.10,w_{0}=1.40$ and (blue lines)
$\beta=-0.85,w_{0}=-0.13$ . All the regions are shown in Fig. 8.

###  II.4 Spectrum of linear Schrödinger equation

We first consider the linear spectrum of $L$ with $\mathcal{PT}$-symmetric
potential (12) as

$\begin{array}[]{l}L\psi(x)=\mu\psi,\\\
L\equiv-\frac{d^{2}}{dx^{2}}+V(x)+\textrm{i}W(x),\end{array}$ (17)

where $\mu$ is the eigenvalue and $\psi(x)$ is the corresponding
eigenfunction. For the bound states
$\lim\limits_{x\rightarrow\infty}|\psi(x)|\rightarrow 0$. Now spectrum of the
linear operator $L$ with complex $\mathcal{PT}$-symmetric potential (12) may
be real or complex conjugates. Then, one can find $v_{0},w_{0},\beta$, for
which linear operator has real spectrum and it is possible, if [60]

$3|w_{0}\beta|\leq\frac{1}{4}+v_{0}(2+\beta^{2}).$ (18)

Fig. 4 exhibits $\mathcal{PT}$ symmetry broken and unbroken phase transition
of the complex potential (12). The region on and above the surface of Fig. 4
is $\mathcal{PT}$ unbroken and the region below the surface is $\mathcal{PT}$
broken.

Figure 4:  Plot of $\mathcal{PT}$-symmetric phase diagram of $v_{0}$ of the
linear operator (17) with respect to $w_{0}$ and $\beta$.

The imaginary parts of some eigen values of the linear operator $L$ (17) are
shown in Fig. 5, in $(w_{0},v_{0})$ plane for $\beta=1$ and in $(w_{0},\beta)$
plane for $v_{0}=0.5$. From Fig. 4 one can find a suitable set of parametrs
for which the potential is $\mathcal{PT}$ unbroken and the spectrum of
operator $L$ 17 are real. First of all under this set of parameters, we will
investigate the linear stability analysis of solutions (13) and (14) and we
denote this point by $P_{1}$ and then we will invesigate further for another
set of parameters, which is denoted by $P_{2}$. Next, we will investigate
their stability analysis for collection of points joing $P_{1}$ and $P_{2}$ by
a curve in $(w_{0},v_{0})$ plane for fixed $\beta$ and in $(w_{0},\beta)$
plane for fixed $v_{0}$ in the next Sec. III.

Figure 5:  Plot of imaginary parts of some eigen values of (17), for (A)-(D)
$v_{0}=0.5$ and for (E)-(H) $\beta=1$.

##  III Linear stability analysis

One of the most important property of solutions of NLSE is their stability.
Here we will examine linear stability of solutions of a NLSE. To this end,
solution in Eq.(4) is given small perturbations $v(x),w(x)$ and is taken to be
of the form [61, 62, 63, 64, 65]

$\begin{array}[]{l}\Psi(x,t)=\\{\psi(x)+\left[\widetilde{v}(x)+\widetilde{w}(x)\right]e^{\lambda
t}\\\
+\displaystyle\left[\widetilde{v}^{*}(x)-\widetilde{w}^{*}(x)\right]e^{\lambda^{*}t}\\}e^{\textrm{i}\mu
t},\end{array}$ (19)

where the superscript $*$ represents the complex conjugation, $|v|,|w|\ll 1$
are infinitesimal normal mode perturbation eigenfunctions, which may grow
during propagation with the perturbation growth rate $\lambda$. Now
substituting the above expression in Eq.(3) and linearizing, one gets the
following coupled set of eigenvalue problem:

$\textrm{i}\left(\begin{array}[]{cc}h_{0}&\nabla^{2}+h_{1}\\\
\nabla^{2}+h_{2}&h_{3}\end{array}\right)\left(\begin{array}[]{c}\widetilde{v}\\\
\widetilde{w}\end{array}\right)=\lambda\left(\begin{array}[]{c}\widetilde{v}\\\
\widetilde{w}\end{array}\right),$ (20)

where

$\begin{array}[]{l}h_{0}=\displaystyle\frac{g}{2}\left(\psi^{2}-\psi^{*2}\right)+i~{}W(x),\\\
h_{1}=-\mu+V(x)+g|\psi|^{2}+\displaystyle
g\left\\{|\psi|^{2}-\frac{1}{2}\left(\psi^{2}+\psi^{*2}\right)\right\\},\\\
h_{2}=-\mu+V(x)+g|\psi|^{2}+\displaystyle
g\left\\{|\psi|^{2}+\frac{1}{2}\left(\psi^{2}+\psi^{*2}\right)\right\\},\\\
h_{3}=-\displaystyle\frac{g}{2}\left(\psi^{2}-\psi^{*2}\right)+i~{}W(x).\end{array}$
(21)

The linear stability is determined by the nature of the eigenvalue problem
(20). If there exists any $\lambda$ with a positive real part, then the
perturbed solution will grow exponentially with $t$ resulting in unstable
mode. On the other hand, mode is completely stable only when, real parts of
$\lambda$ are not positive. The BS and DS are satisfying boundary conditions

$\begin{array}[]{l}|\Psi_{0}(x,t)|\rightarrow
0,|x|\rightarrow\infty,~{}t>0,\\\ \frac{\partial}{\partial
x}\Psi_{0}(x,t)\rightarrow 0,|x|\rightarrow\infty,~{}t>0,\end{array}$ (22)
$\begin{array}[]{l}|\Psi_{1}(x,t)|\rightarrow
C_{1}|\delta|,|x|\rightarrow\infty,~{}t>0,w_{0}^{2}<\frac{1}{4\beta^{2}},\\\
\frac{\partial}{\partial x}\Psi_{1}(x,t)\rightarrow
0,|x|\rightarrow\infty,~{}t>0.\end{array}$ (23)

In this paper, we will apply the Fourier collocation method (FCM) and finite
difference method (FDM) for finding eigenvalues of Eq.(20) for BS and DS
respectively. The simulation of these two solutions are obtained by
pseudospectral method and finite difference method. Using the definition of
power and the VK statbility condition [66], one can say that, BS defined in
(13) is stable, if $g>0$ and unstable, if $g<0$, for the real potential i.e.,
when $w_{0}=0$. The stable region of BS and DS can not be defined for
arbitrary values of $v_{0},w_{0}$ and $\beta$. Therefore, we will find stable
regions in $(w_{0},v_{0})$ plane for fixed $\beta$ and in $(w_{0},\beta)$
plane for fixed $v_{0}$ by applying some numerical techniques.

###  III.1 Numerical method

All the special grid points are defined by $x_{j}=-L+\delta xj$,
$j=1,2,3,...,N+1$ ($L$ being the half-width), where $\delta x=2L/N$ is taken
to be the lattice spacing (resolution). The left and right boundary points are
denoted by $j=1$ and $j=N+1$, respectively. The boundary conditions for BS and
DS at the end points are given by

$\begin{array}[]{ll}\displaystyle|\Psi_{0}(\pm
L,t)|=0,~{}\displaystyle\left[\partial_{x}\Psi_{0}(x,t)\right]_{x=\pm
L}=0,\end{array}$ (24)

and

$\begin{array}[]{l}\displaystyle|\Psi_{1}(\pm
L,t)|=C_{1}|\delta|,\displaystyle\left[\partial_{x}\Psi_{1}(x,t)\right]_{x=\pm
L}=0,w_{0}^{2}<\frac{1}{4\beta^{2}}.\end{array}$ (25)

In this paper, we will consider two finite difference methods for discretizing
Eqs. (3), (24) and (25). For BS (13), we will apply (i) speudospectral method
for second-order spatial derivatives and 4th order Runge-Kutta (RK4) for the
temporal derivative and (ii) second-order central-difference formula for
second-order spatial derivatives and Crank-Nicholson finite difference method
for the temporal derivative. Similarly, for DS we will apply (iii) second-
order central-difference formula for second-order spatial derivatives and 4th
order Runge-Kutta (RK4) for the temporal derivative and (iv) second-order
central-difference formula for second-order spatial derivatives and Crank-
Nicholson finite difference method for the temporal derivative [67]. Then, a
solution will be considered a stable mode, if none of eigen frequencies
$\lambda=\lambda_{R}+i\lambda_{I}$ has a positive real part $\lambda_{R}$.
Finally, the stable regions will be obtained from the eigenvalue problem (20)
and they will be checked by direct dynamical evolution of the NLSE (3) forward
with respect to $t$. To validate linear stability, we will add random noise
perturbation for stable mode add eigen vector
$(\widetilde{v},\widetilde{w})^{T}$ of (20), which corresponds to largest
positive real part $\lambda_{R}$ for unstable mode to initial solution. The
random numbers are generated from $(0,1)$.

###  III.2 Stability of BS

We now consider two cases of $g$. In the first one, $g=-1$, the existence of
the solution (13) of the NLSE (5) in $(w_{0},v_{0})$ and $(w_{0},\beta)$
planes are shown in Fig. 6 (A) -(B). The space $(w_{0},v_{0})$ is divided into
seven regions $I_{i},~{}i=1-5,~{}J_{1}$ and $J_{2}$. In region $I_{1}$,
solution (13) does not exist. In addition, $\mathcal{PT}$ symmetric phase
transition is shown in Fig. 6 (A). In green region $(I_{5}\cup J_{1}\cup
J_{2})$ linear operator $L$ (17) is $\mathcal{PT}$ broken and in white region
$(\cup_{i=1}^{4}I_{i})$ $L$ is $\mathcal{PT}$ unbroken. The stable region for
solution (13) is shown in Fig. 6 (A) and it is bounded by blue and red curves,
which is equal to $I_{3}\cup J_{1}\cup J_{2}$. From this figure it is clear
that $I_{3}$ is $\mathcal{PT}$ unbroken stable, whereas $J_{1}\cup J_{2}$ is
$\mathcal{PT}$ broken stable, which is very small with respect to $I_{3}$ and
total stable region $I_{3}\cup J_{1}\cup J_{2}$ is connected. The
$(w_{0},\beta)$ space is divided into nine regions $(I_{i},~{}i=1-9)$ and they
are shown in Fig. 6 (B). From this figure one can see that, solution (13) does
not exists in $I_{3}\cup I_{9}$. The linear operator (17) is $\mathcal{PT}$
broken in green region $(I_{1}\cup I_{5})$, whereas it is unbroken in white
region $\left((\cup_{i=2}^{4}I_{i})\cup(\cup_{j=6}^{9}I_{j})\right)$. It is
interesting to observe that solution (13) is stable in a connected region
$I_{7}$, which entirely lies in $\mathcal{PT}$ unbroken and there is no stable
$\mathcal{PT}$ broken mode. [64, 65].

Figure 6: (A) $\mathcal{PT}$-symmetric phase transitions of the linear
operator (17) with potential (12) and stable-unstable regions of the solution
(13) in $(w_{0},v_{0})$ plane, for $g=-1,\beta=1.0$. (B) The phase transitions
for the same linear operator in $(w_{0},\beta)$ plane and stable-unstable mode
of the solution (13) in $(w_{0},\beta)$ plane, for $g=-1,v_{0}=0.5$.

In the second case $g=1$, the existence of BS in $(w_{0},v_{0})$ and
$(w_{0},\beta)$ planes are shown in Fig. 7 (A) and (C). From Fig. 7 (A) we see
that, BS does not exist below the red curve and the operator $L$ of (17) is
broken in green region, which is entirely below the red curve. BS is stable in
two disjoint regions $J_{1}$ and $J_{2}\subset I_{1}$ and they are shown in
Fig. 7 (B). From Fig. 7 (C) we see that, BS exists in $I_{2}$ which is the
upward interior of the red parabola and in $I_{5}$ which is the downward
interior of red parabola. The stable regions exist in four small parts such as
$J_{3}$, $J_{4}\subset I_{2}$ and $J_{5}$, $J_{6}\subset I_{5}$ and they are
shown in Fig. 7 (D) and (E). All parts are connected but
$\displaystyle\cup_{i=3}^{6}J_{i}$ is disconneted. For $g=1$, stable region in
both planes are very small, which agree with [46].

Figure 7: (A) Phase transition of the linear operator (17) with potential (12)
and stable-unstable mode of the solution (13) in $(w_{0},v_{0})$ plane for
$g=1,\beta=1$. (C) The phase transitions for the same linear operator and
stable-unstable mode of the solution (13) in $(w_{0},\beta)$ plane, for
$g=1,v_{0}=0.5$.

###  III.3 Stability of DS

For $g=-1$ the DS is stable in small region in $(w_{0},v_{0})$ and
$(w_{0},\beta)$ planes satisfying ${w_{0}}^{2}<\frac{1}{4\beta^{2}}$ and some
additional conditions. In particular, $\beta=1,\delta=1$, we have found a
stable region in $(w_{0},v_{0})$ plane, such that $\frac{2+w_{0}^{2}}{3}\leq
v_{0}\leq\frac{2+w_{0}^{2}}{3}+\epsilon$, where $0\leq\epsilon<0.1$. We now
consider the second case $g=1$, DS exists in a domain which lies between two
vertical lines $w_{0}^{2}<\frac{1}{4\beta^{2}}$ and a paraboic curve
$v_{0}<\frac{w_{0}^{2}\beta^{2}+2}{\beta^{2}+2}$ in $(w_{0},v_{0})$ plane. In
particular, for $\beta=1$, $\delta=1$, the region is shown in Fig. 8 (A) and
the linear operator is broken in green region $\cup_{i=5}^{7}I_{i}$. Stable
modes of DS are marked by _blue dots, $I_{3}$_ and _magenta dots, $I_{4}$_
represent unstable modes. It is to be found that, stable modes are exist
within $\mathcal{PT}$ unbroken region. We did not find stable mode in broken
region after several trials. Similarly, existence of (14) and $\mathcal{PT}$
broken-unbroken phase transition shown in Fig. 8 (B) in $(w_{0},\beta)$ plane,
for $v_{0}=0.5,\delta=1$. The space $(w_{0},\beta)$ is divided into eleven
regions $I_{i},~{}i=1-11$. The linear operator $L$ is broken in $I_{5}\cup
I_{6}\cup I_{9}\cup I_{10}$ and unbroken in
$\displaystyle\cup_{i=1}^{4}I_{i}\cup I_{7}\cup I_{8}\cup I_{11}$ and DS
exists in a region bounded by red curves. Stable modes of DS are marked by
_blue dots, $I_{8}$_ and unstable in $\displaystyle\cup_{i=1}^{7}I_{i}$. The
DS has stable modes only in $\mathcal{PT}$ unbroken region. To confirm all
these stable modes of BS and DS, we will investigate their excitations in the
next section IV.

Figure 8: (A) Phase transition of the linear operator $L$ (17) with potential
(12) and stable-unstable mode of the solution (14) in $(w_{0},v_{0})$ plane,
$g=1,\beta=1,\delta=1$. (B) Phase transition of the same linear operator and
stable-unstable mode of the solution (14) in $(w_{0},\beta)$ plane, for
$g=1,v_{0}=0.5,\delta=1$.

##  IV Excitation

Now, we will focus on nature of solutions of time dependent NLSE with complex
$\mathcal{PT}$-symmetric potential

$i\,\Psi_{t}=-\Psi_{xx}+\left(V(x,t)+i\,W(x,t)\right)\Psi+g|\Psi|^{2}\Psi,$
(26)

where $V(x,t),W(x,t)$ are given by Eq.(12) with $v_{0}\rightarrow v_{0}(t)$
and $w_{0}\rightarrow w_{0}(t)$ satisfy [68]

$\left\\{v_{0},w_{0},\beta\right\\}(t)=\left\\{\begin{array}[]{ll}\left(\left\\{v_{02},w_{02},\beta_{02}\right\\}-\left\\{v_{01},w_{01},\beta_{01}\right\\}\right)f(t)\\\
+\left\\{v_{01},w_{01},\beta_{01}\right\\},0\leq t\leq 5000\\\
\left\\{v_{02},w_{02},\beta_{02}\right\\},t>5000\end{array}\right.,$ (27)

where $f$ is a real valued function of $t$. It is easy to verify that,
solutions $(\ref{solpsi0})and(\ref{solpsi1})$ do not satisfy Eq.(26) but they
satisfy for initial time $t=0$ and $t\geq 5000$. In this paper we have
considered $f(t)=\sin\frac{\pi t}{10000}$, $0\leq t\leq 5000$.

Figure 9:  (A) Shows time dependent potential amplitudes $w_{0},v_{0}$ for
different curves. (B) Shows that stable and unstable regions of BS and four
curves $C_{i},~{}i=1-4$. Excitations of stable nonlinear localized state of
(26) with intial condition (13), for $g=-1$, $\beta=1$, (C)
$v_{01}=0.062,v_{02}=0.1335,w_{01}=0.15,w_{02}=-0.22$
(D)$v_{01}=0.1335,v_{02}=0.1335,w_{01}=-0.02,w_{02}=-0.22$ (E)
$v_{01}=0.08,v_{02}=0.53,w_{01}=-0.08,w_{02}=0.022$ (F)
$v_{01}=0.08,v_{02}=0.08,w_{01}=-0.08,w_{02}=0.022$.

###  IV.1 Excitation of BS

We have displayed excitations of BS of (26), with initial condition (13) in
Fig. 9 (C) -(F), for $g=-1$ along four curves $C_{1}$ \- $C_{4}$. Each curve
has initial point $P_{1}(w_{01},v_{01})$ and final point
$P_{2}(w_{02},v_{02})$. The time dependent potential amplitudes
$\min\\{v_{01},v_{02}\\}\leq v_{0}(t)\leq\max\\{v_{01},v_{02}\\}$ and
$\min\\{w_{01},w_{02}\\}\leq w_{0}(t)\leq\max\\{w_{01},w_{02}\\}$ satisfy the
Eq.(27) and they are plotted in Fig. 9 (A) with respect to $t$ and the
corresponding locus of $(w_{0}(t),v_{0}(t))$ for $0\leq t\leq 5000$, represent
curves $C_{j}~{},j=1,2,3,4$ and they are shown in Fig. 9 (B). The curve
$C_{1}:P_{2}(-0.22,0.1335)\nwarrow P_{1}(0.15,0.062)$ is started from
$P_{1}(\in J_{2},~{}\mbox{broken~{}region})$ to $P_{2}(\in
J_{1},~{}\mbox{broken~{}region})$ via unbroken region $(I_{3})$, whereas
$C_{2}:P_{2}(-0.22,0.1335)\longleftarrow P_{1}(-0.02,0.1335)$ is started from
$P_{1}(\in I_{3},~{}\mbox{unbroken~{}region})$ to $P_{2}(\in
J_{1},~{}\mbox{broken~{}region})$. Therefore, excitation Fig. 9 (C) started
from a broken point passes through unbroken region and then eneded to a broken
point, whereas, excitation Fig. 9 (D) started from an unbroken point and then
going to broken region. On the otherhand Fig. 9 (E) and (F) are excitations
along $C_{3}:P_{1}(-0.08,0.08)\nearrow P_{2}(0.022,0.53)$ and
$C_{4}:P_{1}(-0.08,0.08)\longrightarrow P_{2}(0.022,0.08)$ and they lie within
unbroken region $I_{3}$. One can see that, curves $C_{1}$, $C_{3}$ are look
like oblique whereas, $C_{2}$, $C_{4}$ are like horizontal but their
directions are different and all excitations in $(w_{0},v_{0})$ plane are
stable. We have checked that, BS is stable for all $(w_{0},v_{0})\in I_{3}\cup
J_{1}\cup J_{2}$ and its excitation is stable along any curve in any direction
which lies in $I_{3}\cup J_{1}\cup J_{2}$.

Figure 10:  Excitations of stable nonlinear localized state of (26) with
initial condition (13), for $g=-1$, $v_{0}=0.5$, (A)
$\beta_{01}=-1.30,\beta_{02}=1.30,w_{01}=-0.017,w_{02}=0.017$ (B)
$\beta_{01}=-1.30,\beta_{02}=1.30,w_{01}=-0.017,w_{02}=-0.017$ (C)
$\beta_{01}=-0.02,\beta_{02}=0.02,w_{01}=-4.7,w_{02}=4.7$ (D)
$\beta_{01}=-0.02,\beta_{02}=-0.02,w_{01}=-4.7,w_{02}=4.7$.

Similarly, in $(w_{0},\beta)$ plane for $g=-1$, $v_{0}=0.5$ excitations of
(26) with initial solution (13) are shown in Figs. 10 (A) -(D) along four
curves $(C_{5}:P_{1}(-0.017,-1.3)\nearrow
P_{2}(0.017,1.3);~{}C_{6}:P_{1}(-0.017,-1.3)\uparrow
P_{2}(-0.017,1.3);~{}C_{7}:P_{1}(-4.7,-0.02)\nearrow
P_{2}(4.7,0.02);~{}C_{8}:P_{1}(-4.7,-0.02)\rightarrow P_{2}(4.7,-0.02))$. From
Figs. 10 (A) -(D), one can see that BS and its excitations are stable. We
observe that, excitation of BS of (26) with intial state (13) is stable along
any curve in any direction which lies in $\mathcal{PT}$ unbroken region
$I_{7}$ of Fig. 6 (B) in $(w_{0},\beta_{0})$ plane.

Figure 11:  Excitations of nonlinear localized state of (26) with initial
condition (13) in $(w_{0},v_{0})$ plane for $g=1$, $\beta=1$, (A)
$v_{01}=0.731,v_{02}=0.731,w_{01}=0.024,w_{02}=-0.024$; (B)
$v_{01}=0.731,v_{02}=0.731,w_{01}=0.024,w_{02}=0.076$.

Figure 12:  Excitations of stable nonlinear localized state of (26) with
initila condition (13) in $(w_{0},\beta)$ plane for $g=1$, $v_{0}=0.5$, (A)
$\beta_{01}=-1.550,\beta_{02}=-1.550,w_{01}=-0.016,w_{02}=0.016$; (B)
$w_{01}=-0.015,\beta_{01}=1.43,w_{02}=0.017,\beta_{02}=1.44$

Similarly, in $(w_{0},v_{0})$ for $g=1$, excitations of BS with same initial
condition (13) along $C_{9}:P_{2}(-0.024,0.731)\leftarrow
P_{1}(0.024,0.731)~{}(\in J_{1})$ and $C_{10}:P_{1}(0.024,0.731)~{}(\in
J_{1})\rightarrow P_{2}(0.076,0.731)$ are shown in Fig. 11 (A) and (B)
respectively. The curve $C_{9}$ lies in the stable region $J_{1}$ of Fig. 7
(B) and the curve $C_{10}$ is started from a stable point of $J_{1}$ and to an
unstable point of $I_{1}$ of Fig. 7 (B). For $C_{9}$ BS is stable and for
$C_{10}$ it is unstable and we observe that, BS is stable for all
$(w_{0},v_{0})\in J_{1}\cup J_{2}$, where $J_{1}$ and $J_{2}$ are connected
regions but $J_{1}\cup J_{2}$ is disconnected.

Similarly, in $(w_{0},\beta)$ plane for $g=1,$ $(w_{0},\beta)$ plane
excitations of BS of (26) along $C_{11}:P_{1}(-0.016,-1.55)~{}(\in
J_{4})\rightarrow P_{2}(0.016,-1.55)~{}(\in J_{4})$ and
$C_{12}:P_{1}(-0.015,1.43)\in J_{6}\nearrow P_{2}(0.017,1.44)\in I_{2}$ curves
are shown in Fig. 12 (A) and (B) respectively. Here, $C_{11}$ lies in a stable
region $J_{4}$ of Fig. 7 (D) and for $C_{12}$, some part lies in a stable
region $J_{6}$ and the remaining part lies in unstable region $I_{2}$ of Fig.
7 (E). Once can sees that, BS is stable for $C_{11}$ and it is unstable for
$C_{12}$. Similarly, BS is stable in $\cup_{i=3}^{6}J_{i}$.

###  IV.2 Excitation of DS

In this case, time dependent potential amplitudes $w_{0}(t)$ and $v_{0}(t)$
are shown in Fig. 13 (A), for four curves
$C_{13}:P_{1}(-0.2,0.25)\longrightarrow
P_{2}(0.2,0.25);~{}C_{14}:P_{1}(-0.2,0.65)\longrightarrow
P_{2}(0.2,0.65);~{}C_{15}:P_{1}(-0.2,0.25)\uparrow
P_{2}(-0.2,0.32);~{}C_{16}:P_{1}(-0.2,0.25)\nearrow P_{2}(0.2,0.65)$, and
curves are shown in Fig. 13 (B) and excitations of DS of (26) with intial
condition (14) along these curves are shown in Fig. 13 (C) - (F), for
$g=\beta=\delta=1$ in $(w_{0},v_{0})$ plane. One can see that, curves
$C_{13}$, $C_{14}$ and $C_{15}$ are passed through stable modes (_blue dots_)
but $C_{16}$ passes through some stable (_blue dots_) and unstable modes
(_magenta dots_) and DS is stable along $C_{13}$, $C_{14}$, $C_{15}$ and
unstable along $C_{16}$. Therefore, DS is stable in a disconneted
$\mathcal{PT}$ unbroken region in $(w_{0},v_{0})$ plane.

Similarly, excitations of DS of (26) in $(w_{0},\beta)$ plane are evaluated
along $(C_{17}:P_{1}(-0.1,-0.85)\longrightarrow
P_{2}(0.1,-0.85);~{}C_{18}:P_{1}(-0.21,-1.4)\longrightarrow
P_{2}(0.21,-1.4);~{}C_{19}:P_{1}(0.1,0.85)\nearrow
P_{2}(0.2,1.15);~{}C_{20}:P_{1}(-0.1,-0.85)\searrow P_{2}(0.21,-1.4))$ are
shown in Fig. 14 (C) - (F), for $g=1,v_{0}=0.5,\delta=1$ and curves are shown
in Fig. 14 (B) and the corresponding $w_{0}(t),\beta(t)$ are shown in Fig. 14
(A). One can see that DS is stable along $C_{17}$, $C_{18}$, $C_{19}$ and
unstable for $C_{20}$, where the curve $C_{20}$ passes through stable and
unstable modes of DS. Theofore, DS is stable in a disconnected $\mathcal{PT}$
unbroken region in $(w_{0},\beta)$ plane.

Figure 13:  (A) Shows time dependent potential amplitudes $w_{0},v_{0}$ for
different curves. (B) Shows that stable (_blue dots_) and unstable (_magenta
dots_) modes of DS and four curves $C_{i},~{}i=13-16$. Excitations of stable-
unstable nonlinear delocalized state of (26) with initial condition (14), for
$g=1$, $\beta=1,\delta=1$, (C)
$v_{01}=0.25,v_{02}=0.25,w_{01}=-0.20,w_{02}=0.20$; (D)
$v_{01}=0.65,v_{02}=0.65,w_{01}=-0.20,w_{02}=0.20$; (E)
$v_{01}=0.25,v_{02}=0.32,w_{01}=-0.20,w_{02}=-0.20$; (F)
$v_{01}=0.25,v_{02}=0.65,w_{01}=-0.20,w_{02}=0.20$.

Figure 14:  (A) Shows time dependent potential amplitudes $w_{0},v_{0}$ for
different curves. (B) Shows that stable (_blue dots_) and unstable (_magenta
dots_) modes of DS and four curves $C_{i},~{}i=17-20$. Excitations of stable-
unstable nonlinear delocalized state of (26) with initial condition (14), for
$g=1$, $v_{0}=0.5,\delta=1$, (C)
$\beta_{01}=-0.85,\beta_{02}=-0.85,w_{01}=-0.10,w_{02}=0.10$; (D)
$\beta_{01}=-1.40,\beta_{02}=-1.40,w_{01}=-0.21,w_{02}=0.21$; (E)
$\beta_{01}=0.85,\beta_{02}=1.15,w_{01}=0.10,w_{02}=0.20$; (F)
$\beta_{01}=-0.85,\beta_{02}=-1.40,w_{01}=-0.10,w_{02}=0.21$.

##  V Conclusion

In this paper, a cubic NLSE with $\mathcal{PT}$-symmetric potential has been
considered. A bright soliton and a bright-dark soliton interaction are
generated in terms of potential amplitudes $v_{0}$ ,$w_{0}$ and $\beta$.
Depending on these parameters, bright-dark soliton interaction reduces to a BS
if $w_{0}^{2}>\frac{1}{4\beta^{2}}$ and a DS, if
$w_{0}^{2}<\frac{1}{4\beta^{2}}$. A phase transition occurs between BS and DS.
In addition, $\mathcal{PT}$ broken-unbroken phase transition is shown with
respect to $w_{0},v_{0},\beta$ and also it is shown in $(w_{0},v_{0})$ and
$(w_{0},\beta)$ planes. The linear stability analysis of BS and DS are
investigated and obtained stable regions of BS in $\mathcal{PT}$ broken
unbroken region and of DS in $\mathcal{PT}$ unbroken region. We have found
that, (i) BS is stable in two connected regions $I_{3}\cup J_{1}\cup J_{2}$
(see Fig. 6 (A)) in $(w_{0},v_{0})$ plane and $I_{7}$ (see Fig. 6 (B)) in
$(w_{0},\beta)$ plane for $g=-1$. (ii) BS is stable in two disconnected
regions $J_{1}\cup J_{2}$ (see Fig. 7 (B)) in $(w_{0},v_{0})$ plane and
$\cup_{i=3}^{6}J_{i}$ (see Fig. 7 (D), (E)) in $(w_{0},\beta)$ plane for
$g=1$. (iii) DS is stable in a small connected region in $(w_{0},v_{0})$ plane
such that $\frac{2+w_{0}^{2}}{3}\leq v_{0}\leq\frac{2+w_{0}^{2}}{3}+\epsilon$,
where $0\leq\epsilon<0.1$, ${w_{0}}^{2}<\frac{1}{4}$, for
$g=-1,\beta=\delta=1$. (iv) DS is stable in two disconnected regions $I_{3}$
(see Fig. 8 (A)) in $(w_{0},v_{0})$ plane and in $I_{8}$ (see Fig. 8 (B))
$(w_{0},v_{0})$ plane for $g=1$. (v) Solution (14) is always unstable for
$w_{0}^{2}>\frac{1}{4\beta^{2}}$. All these regions are verified by eigenvalue
problem (20), then verified by direct numerical simulation, after that they
are verified by exacitaions of BS and DS of time dependent nonlinear
Schrödinger equation with time dependent potential amplitutes. The curves are
taken with different length in different directions, which are defined in
their respective domain. For examples only twenty excitaions are shown along
twenty curves $C_{i},~{}i=1-20$ and just twelve curves $C_{i},~{}i=1-4,13-20$
are shown among them.

## Acknowledgement

Debraj Nath dedicates this article to the memory of his kind brother, late Raj
Kumar Nath. Amiya Das gratefully acknowledges financial support from SERB-DST,
Govt. of India (EEQ/2017/000150) and DST PURSE-II University of Kalyani.

## References

  * [1] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80 (1998) 5243. C. M. Bender, D. C. Brody, H. F. Jones, Phys. Rev. Lett. 89 (2002) 270401. C. M. Bender, S. Boettcher, P. N. Meisinger, J. Math. Phys. 40 (1999) 2201. C. M. Bender, Rep. Prog. Phys. 70 (2007) 947.
  * [2] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, Z. H. Musslimani, Phys. Rev. A 81 (2010) 063807.
  * [3] C. E. Ruter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, D. Kip, Nat. Phys. 6 (2010) 192. A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, D. N. Christodoulides, Phys. Rev. Lett. 103 (2009) 093902. Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, D. N. Christodoulides, Phys. Rev. Lett. 106 (2011) 213901.
  * [4] Z. H. Musslimani, K. G. Makris, R. El-Ganainy and D. N. Christodoulides, Phys. Rev. Lett. 100 (2008) 030402.
  * [5] M. V. Berry, J. Phys. A 41 (2008) 244007.
  * [6] S. Longhi, Phys. Rev. Lett. 103 (2009) 123601.
  * [7] O. Bendix, R. Fleischmann, T. Kottos and B. Shapiro, J. Phys. A. 43 (2010) 265305.
  * [8] K. G. Makris, Z. H. Musslimani, D. N. Christodoulides and S. Rotter, Nat. Commun. 6 (2014) 7257.
  * [9] J. Schindler, Z. Lin, J. M. Lee, H. Ramezani, F. M. Ellis and T. Kotto, H. Ramezani, F. M. Ellis and T. Kotto, J. Phys. A. 45 (2012) 444029.
  * [10] H. Cartarius and G. Wunner, Phys. Rev. A 86 (2012) 013612.
  * [11] E.-M. Graefe, J. Phys. A 45 (2012) 444015.
  * [12] D. Dast, D. Haag, H. Cartarius, J. Main and G. Wunner, J. Phys. A. 46 (2013) 375301.
  * [13] G. Castaldi, S. Savoia, V. Galdi, A. Al‘u and N. Engheta, Phys. Rev. Lett. 110 (2013) 173901.
  * [14] N. Lazarides and G. P. Tsironis, Phys. Rev. Lett. 110 (2013) 053901.
  * [15] M. Kang, F. Liu and J. Li, Phys. Rev. A 87 (2013) 053824.
  * [16] H. Alaeian and J. A. Dionne, Phys. Rev. A 89 (2014) 033829.
  * [17] G. P. Tsironis and N. Lazarides, App. Phys. 115 (2014) 449.
  * [18] T. Kottos, Nat. Phys. 6 (2010) 166.
  * [19] C. E. Rüter, K. G. Markis, R. El-Ganainy, D. N. Christodoulides, M. Segev and D. Kip, Nat. Phys. 6 (2010) 192.
  * [20] L. Feng et al., Nature Mater. 12 (2013) 108.
  * [21] J. Schindler, A. Li, M. C. Zheng, F. M. Ellis and T. Kottos, Phys. Rev. A 84 (2011) 040101(R).
  * [22] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, Phys. Rev. Lett. 100 (2008) 103904 . A. Ruschhaupt, F. Delgado, J. G. Muga, J. Phys. A: Math. Gen. 38 (2005) L171. R. El-Ganainy, K. G. Makris, D. N. Christodoulides, Z. H. Musslimani, Opt. Lett. 32 (2007) 2632. M. V. Berry, J. Phys. A: Math. Theor. 41 (2008) 244007. S. Klaiman, U. Gunther, N. Moiseyev, Phys. Rev. Lett. 101 (2008) 080402. S. Longhi, ibid. 103 (2009) 123601. S. Longhi, Phys. Rev. A 82 (2010) 031801(R). H. Ramezani, D. N. Christodoulides, V. Kovanis, I. Vitebskiy, T. Kottos, Phys. Rev. Lett. 109 (2012) 033902. M. A. Miri, A. Regensburger, U. Peschel, D. N. Christodoulides, Phys. Rev. A 86 (2012) 023807.
  * [23] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, D. N. Christodoulides, Phys. Rev. Lett. 103 (2009) 093902. C. E. Reuter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, D. Kip, Nat. Phys. 6 (2010) 192. A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, U. Peschel, Nature 488 (2012) 167.
  * [24] A. Khare, S. M. Al-Marzoug, H. Bahlouli, Phys. Lett. A 376 (2012) 2880.
  * [25] D. A. Zezyulin, V. V. Konotop, Phys. Rev. A 85 (2012) 043840.
  * [26] B. Midya, R. Roychoudhury, Phys. Rev. A 87 (2013) 045803.
  * [27] S. Hu, X. Ma, D. Lu, Z. Yang, Y. Zheng, W. Hu, Phys. Rev. A 84 (2011) 043818.
  * [28] S. Hu, G. Liang, S. Cai, D. Lu, Q. Guo, W. Hu, arXiv:1203.1862v1 (2012)
  * [29] C. P. Jisha, L. Devassy, A. Alberucci, V. C. Kuriakose, Phys. Rev. A 90 (2014) 043855.
  * [30] B. Midya, Nonlinear Dyn. 79 (2015) 409.
  * [31] P. Li, B. Liu, L. Li, D. Mihalache, Rom. J. Phys. 61 (2016) 577.
  * [32] P. Li, L. Li, D. Mihalache, arXiv:1804.02629v1 [physics.optics] (2018).
  * [33] Z. C. Wen, Z. Yan, Phys. Lett. A 379 (2015) 2025.
  * [34] Y. Chen, Z. Yan, Commun. Nonlin. Sc. Numer. Simulat. 57 (2018) 34.
  * [35] N. N. Akhmediev and A. Ankiewicz, Solitons: Nonlinear pulses and beams (Chapman and Hall, London, 1997)
  * [36] Y. S. Kivshar and G. P. Agrawal, Optical Solitons: From Fibers to Photonic Crystals (Academic, San Diego, 2003).
  * [37] Z. H. Musslimani, K. G. Makris, R. El-Ganainy and D. N. Christodoulides, J. Phys. A 41 (2008) 244019.
  * [38] K. G. Makris, R. El-Ganainy, D. N. Christodoulides and Z. H. Musslimani, Int. J. Theo. Phys. 50 (2011) 1019.
  * [39] F. Kh. Abdullaev, V. V. Konotop, M. Salerno and A. V. Yulin, Phys. Rev. E 82 (2010) 056606.
  * [40] A. Khare, S. M. Al-Marzoug and H. Bahlouli, Phys. Lett. A 376 (2012) 2880.
  * [41] V. Achilleos, P. G. Kevrekidis, D. J. Frantzeskakis and R. Carrectero-González, Phys. Rev. A 86 (2012) 013808.
  * [42] A. Khare and A. Saxena, Phys. Lett. A 377 (2013) 2761.
  * [43] M. Salerno, J. Geom. Symm. Phys. 32 (2013) 25.
  * [44] T. Mayteevarunyoo, B. A. Malomed and A. Reoksabutr, Phys. Rev. E 88 (2013) 022919.
  * [45] H. Xu, P. G. Kevrekidis, Q. Zhou, D. J. Frantzeskakis, V. Achilleos and R. Carretero-González, Rom. J. Phys. 59 (2014) 185.
  * [46] B. Midya and R. Roychoudhury, Ann. Phys. 341 (2014) 12.
  * [47] H. Chen, S. Hu and L. Qi, Opt. Commun. 331 (2014) 139.
  * [48] V. V. Konotop and D. A. Zezyulin, Opt. Express 39 (2014) 553.
  * [49] A. K. Sarma, M.-A. Miri, Z. H. Musslimani and D. N. Christodoulides, Phys. Rev. A 89 (2014) 052918.
  * [50] C.-Q. Dai and Y. -Y.Wang, Opt. Commun. 315 (2014) 303.
  * [51] I. Göskel, N. Antar and I. Bakirtas, Opt. Commun. 354 (2015) 277.
  * [52] Z. Yan, Z.Wen and C. Hang, Phys. Rev. E 92 (2015) 022913.
  * [53] D. Nath, B. Roy, R. Roychoudhury, Chaos, Solitons & Fractals 81 (2015) 91.
  * [54] A. Das, N. Ghosh and D. Nath, Phys. Lett. A 384 (2020) 126681.
  * [55] D. Nath and P. Roy, J. Nonlinear Optic. Phys. Mat. 25 (2016) 1650036.
  * [56] D. Nath and P. Roy, Phys. of Part. and Nuclei Lett. 14 (2017) 347.
  * [57] F. Cooper, A. Khare and U. Sukhatme, Supersymmetry in Quantum Mechanics, World Scientific (2001).
  * [58] Y.S. Kivshar, B. Luther-Davies, Phys. Rep. 298 (1998) 81197. R. Heidemann, S. Zhdanov, R. Stterlin, H.M. Thomas, G.E. Morfill, Phys. Rev. Lett. 102 (2009) 135002. T. Xu, B. Tian, L.L. Li, X. Lv, C. Zhang, Phys. Plasmas 15 (2008) 102307. T. Tsurumi, M. Wadati, J. Phys. Soc. Jpn. 97 (1998) 2294. M. Wadati, N. Tsuchida, J. Phys. Soc. Jpn. 75 (2006) 014301. D. Mandelik, R. Morandotti, J.S. Aitchison, Y. Silberberg, Phys. Rev. Lett. 92 (2004) 093904.
  * [59] Z. Shi, X. Jiang, X. Zhu and H. Li, Phys. Rev. A 84 (2011) 053855.
  * [60] Z. Ahmed, Phys. Lett. A 282 (2001) 343.
  * [61] J. Yang, _Nonlinear Waves in Integrable and Nonintegrable Systems_ (SIAM, Philadelphia, 2010).
  * [62] S. Nixon, L. Ge and J. Yang, Phys. Rev. A. 85 (2012) 023822.
  * [63] D. A. Zezyulin, Y. V. Kartashov and V. V. Konotop, Euro. Phys. Lett. 96 (2011) 64003.
  * [64] D. Nath, B. Roy, R. Roychoudhury, Opt. Commun. 393 (2017) 224.
  * [65] D. Nath, N. Saha, B. Roy, Eur. Phys. J. Plus 133 (2018) 504.
  * [66] M. Vakhitov and A. Kolokolov, Radiophys. Quantum Electron. 16 (1973) 783.
  * [67] W. Bao, Q. Tang, and Z. Xu, J. Comput. Phys. 235 (2013) 423. D. Nath, Y. Gao, R.B. Mareeswaran, T. Kanna, and B. Roy, Chaos 27 (2017) 123102.
  * [68] Z. Yan, Y. Chen, Chaos 27 (2017) 073114. Z. Yan, Z. Wen, V. V. Konotop, Phys. Rev. A 92 (2015) 023821. Y. Chen and Zhenya Yan, Phys. Rev. E 95 (2017) 012205. Y. Chen and Z. Yan, Scientific Reports 6:23478 DOI: 10.1038/srep23478, (2016); Y. Chen, Z. Yan, D. Mihalache and B. A. Malomed, Scientific Reports 7: 1257 DOI:10.1038/s41598-017-01401-3, (2017).