# Predicting Succinylation Sites in Proteins with Improved Deep Learning
Architecture

Olusola T. Odeyomi and Gergely Zaruba  Department of Electrical Engineering
and Computer Science
Wichita State University
Wichita, Kansas 67260–0083
Email<EMAIL_ADDRESS>
gergely.zaruba @wichita.edu

###### Abstract

Post-translational modifications (PTMs) in proteins occur after the process of
translation. PTMs account for many cellular processes such as deoxyribonucleic
acid (DNA) repair, cell signaling and cell death. One of the recent PTMs is
succinylation. Succinylation modifies lysine residue from $-1$ to $+1$.
Locating succinylation sites using experimental methods, such as mass
spectrometry is very laborious. Hence, computational methods are favored using
machine learning techniques. This paper proposes a deep learning architecture
to predict succinylation sites. The performance of the proposed architecture
is compared to the state-of-the-art deep learning architecture and other
traditional machine learning techniques for succinylation. It is shown from
the performance metrics that the proposed architecture provides a good trade-
off between speed of computation and classification accuracy.

##  I Introduction

Post-translational modification (PTMs) occur in proteins by the addition of a
covalent bond after protein is synthesized from messenger RNA (mRNA) by
ribosomes. PTMs are very important in many cellular processes such as DNA
repair, cell death and cell signaling. Sites in protein that undergo PTMs are
those that can act as a nucleophile in reactions. There are many types of PTMs
such as phosphorylation, glycosylation, succinylation etc. Protein
succinylation is a recently emerged protein post-translational modification on
lysine. Succinylation introduces a large structural moeity (e.g., 100 Da) and
its modifies the net charge of modified lysine residue from -1 to +1. It is
suggested that protein succinylation is important in cellular metabolism and
many other cellular functions. Experimental methods for determining protein
succinylation sites, such as mass spectrometry, are laborious and not cost-
effective. Computational methods are thus favored.

Most existing work that uses computational methods employs machine learning
techniques such as in iSuc-PseAAC [1], iSuc-PseOPT [2], pSuc-Lys [3],
SuccinSite [4], SuccinSite 2.0 [5], GPSuc [6], and PSuccE [7]. The problem
with traditional machine learning methods is the requirement for manual
extraction of features. Also, traditional machine learning methods cannot
handle very large FASTA file reliably. Most deep learning methods in
literature are applied to phosphorylation rather than succinylation, and use
one-hot encoding that does not account for the proper classification of
features. However, there is a recent work that uses a deep learning
architecture for succinylation with state-of-the-art performance. The authors
called this architecture deepSuccinylSite-embedding [8].

Thus, in this paper, the performance of deepSuccinylSite-embedding is compared
with the performance of a proposed deep learning architecture called
DeepSuccinylCNN+LSTM. While deepSuccinylSite uses 2D convolutional network,
the novel architecture uses 1D Convolutional network + LSTM (Long Short-Term
Memory). The motivation for this architecture is comparable performance with
deepSuccinylSite but with faster computational speed. This is very important
when dealing with large FASTA file.

##  II Dataset Preprocessing

The dataset for this computation was obtained from experimentally derived
lysine residue reported in [6] and [7]. Proteins with more than thirty percent
$30\%$ sequence identity were removed using CD-HIT. This left a remainder of
$5009$ succinylation sites, and $53,542$ non-succinylation sites. Out of
these, $4755$ succinylation sites and $50,565$ non-succinylation sites were
used as training set. This meant that $254$ succinylation sites and $2977$
non-succinylation sites were used as independent test set. However, $5$ out of
the $4755$ positive sites were lost because they contained other characters.
The datasets were balanced using undersampling. The final training set had
$4750$ positive and $4750$ negative sites. The final independent test set
contained $254$ positive and $254$ negative sites. The optimal window size was
$33$. If left and right side were less than half the size of the window, then
a pseudo-residue “-” was used to recover all positive sites. Table 1 shows the
size of the positive and negative sites for the training and independent test
sets.

TABLE I: Size of Training and Independent Test Set  Dataset | Positive | Negative
---|---|---
Training | 4750 | 4750
Independent Test | 254 | 254

##  III Input Modeling

The novel deep learning architecture Conv1D+LSTM used keras embedding layer
similar to DeepSuccinylSite. The embedding layer came before the first
convolutional layer. The embedding layer output was a 21-dimensional vector
space. The embedding layer grouped co-occuring items together. Other types of
sophisticated embedding techniques common in natural language processing could
be used, such as Word2Vec or GloVe, but they often come with increased cost of
computation. The training set was divided into $80\%$ training and $20\%$
validation set. Validation was done at every epoch to prevent overfitting.
Checkpointer was used to select the optimal model from the epochs based on
validation accuracy. The protein sequence dataset was in FASTA format. The
input dimension to the convolutional neural network had a dimension of $33$ by
$21$, since the window size was $33$, and the vector size of the output from
the embedding size was $21$. An advantage of using deep learning is the
elimination of handcrafted features, which is often laborious.

##  IV Deep Learning Architecture

The input to the first layer of the deep learning architecture is the
embedding layer. DeepSuccinylSite used a Lambda layer after this embedding
layer, but the Lambda layer was not included in the architecture used in this
paper. Three convolutional layers were used with $64$, $128$ and $256$ number
of neurons respectively. The hidden layers were for hierarchical feature
extraction. Although, DeepSuccinylSite used $2D$ convolutional layers instead
of $1D$ convolutional layers for more feature extraction, it was observed that
it slowed down the speed of the computation in comparison to the improvement
obtained. Hence, $1D$ convolutional layer was used together with LSTM in
DeepSuccinylCNN+LSTM. The kernel size used was $17$, since the PTM site lied
in the $17^{th}$ position. The use of this kernel size removed the need for
paddling. Dropouts of $0.6$ were inserted between the convolutional layers to
avoid overfitting. The activation function was linear rectifier (ReLu) because
it could minimize overfitting and maximize the predictive power of the model.
A single maxpool layer was used followed by the LSTM layer. Three dense layers
with neurons size $256$, $768$ and $1024$ were used. Dropouts of $0.5$, $0.5$
and $0.25$ were used between the dense layers respectively. The optimization
used was Adam optimization, which is an adaptive stochastic gradient
optimization algorithm. Adam optimization inherently possesses the advantages
of adaptive gradient optimization and root-mean-square propagation. Also,
binary cross entropy was used as the loss function because this is a binary
classification problem. The final output layer contained 2 nodes and used
softmax activation function. Figure 1 is a snapshot of the parameters for the
simulation.

Figure 1: A snapshot of the parameters for DeepSuccinylCNN+LSTM architecture

##  V Results

For this simulation, $10-$fold cross validation was used. This means that the
training set was divided into $9$ parts for training and $1$ part for
validation at every epoch. This process is repeated, so that each part is used
as the validation set. The performance metrics for the simulation are Matthew
correlation, specificity and sensitivity. These performance metrics are
defined below.

###  V-A Sensitivity

Sensitivity is a measure of how many positive succinylation sites were
correctly identified. It is given by the (3) below:

$sensitivity=\frac{TP}{TP+FN}$ (1)

where $TP$ refers to the True Positives and $FN$ refers to the False
Negatives.

###  V-B Specificity

Specificity measures the number of negative succinylation sites correctly
identified.

$specificity=\frac{TN}{TN+FP}$ (2)

where $TN$ refers to the True Negatives and $FP$ refers to the False
Positives. The $TP$, $TN$, $FP$ and $FN$ were gotten from the confusion
matrix.

###  V-C Matthew Correlation Coefficient

The Matthew correlation coefficient (MCC) is a measure of the correlation
between the predicted binary class and the observed binary class. It ranges
from $-1$ to $+1$. A coefficient of $-1$ means that there is a total
disagreement between the predicted and the observed binary classifications. A
coefficient of $0$ means that the prediction is random. A coefficient of $1$
means that there is a perfect agreement between the predicted and the observed
binary classifications. The formula for MCC is given as

$MCC=\frac{(TP\times TN)-(FP\times FN)}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$
(3)

###  V-D Receiver Operating Characteristics Curve

The receiver operating characteristic curve (ROC) is a graphical tool to
investigate the discriminatory ability of a binary classifier. It compares
sensitivity to specificity, with the true positive rate on the y-axis, and the
false positive rate on the x-axis. The diagonal line in the plot divides the
ROC space. Points above this diagonal line depicts good classification
results, while points below this line depicts poor classification results.

###  V-E Area Under ROC Curve

The area under the ROC curve (AUC) measures the entire area area under the ROC
curve. It can be defined as the probability that a binary classifier will rank
higher a randomly chosen positive site than a randomly chosen negative site. A
classifier whose predictions are completely wrong has an AUC of $0.0$, while a
classifier whose predictions are completely correct has an AUC of $1.0$. Thus,
the AUC ranges from $0$ to $1$.

The simulation used the same datasets, window size and embedding size to
compare the performance of the state-of-the-art DeepSuccinylSite deep learning
architecture and the proposed deep learning architecture. Both deep learning
architectures were trained on the training set, but their performance were
measured on the independent test set. The performance of both deep learning
architectures is shown in Table II.

TABLE II: DeepSuccinylSite Versus DeepSuccinylCNN+LSTM Prediction Schemes | Sensitivity | Specificity | MCC
---|---|---|---
DeepSuccinylSite | 0.73 | 0.70 | 0.43
DeepSuccinylSiteCNN+LSTM | 0.77 | 0.63 | 0.39

The ROC curves for DeepSuccinylSite and DeepSuccinylCNN+LSTM are both shown in
Figs. 2 and 3 respectively.

Figure 2: ROC curve for DeepSuccinylSite. Figure 3: ROC curve for
DeepSuccinylCNN+LSTM.

It can be observed from Table II that DeepSuccinylSite has a higher MCC and
specificity than DeepSuccinylCNN+LSTM, but a lower sensitivity. However, the
difference is not much. Also, from Figs. 2 and 3, it can also be observed that
DeepSuccinylSite has a higher AUC than DeepSuccinylCNN+LSTM. The main
advantage of DeepSuccinylCNN+LSTM over DeepSuccinylSite is faster computation.
The average time for running an epoch in DeepSuccinylSite was 29 seconds. The
average time for running an epoch in DeepSuccinylLSTM+CNN was 20 seconds.
However, the number of epoch used was $60$ and the batch size is $256$. This
meant that for $60$ epochs, the average running time difference was $540$
seconds, and for $256$ batch size, the average running time difference was
$38.4$ hours, using the same computing resource, i.e., a Intel core i7 HP
computer. Although, this running time difference will be reduced by using a
graphics processing unit (GPU). However, when the size of the datasets
increase, the running time difference will become obvious.

There are other traditional machine learning methods for predicting
succinylation sites that utilizes other features apart from the amino-acid
sequence. Such features include the physiochemical properties such as the
Pseudo Amino Acid Composition (PAAC), K-spaced Amino Acid Pairs (AAP) etc.
These machine learning methods are iSuc-PseAAC [1], iSuc-PseOpt [2], PSuc-Lys
[3], SuccinSite [4],SuccinSite2.0 [5] GPSuc[6] and PSuccE [7]. Table III shows
that the performance of these machine learning methods in comparison with
DeepSuccinylLSTM+CNN. The values for sensitivity, specificity and MCC for the
machine learning methods were obtained directly from the published articles.
It can be seen from Table III that DeepSuccinylCNN+LSTM still performs better
than these machine learning methods. Hence, DeepSuccinylCNN+LSTM provides a
good trade-off between speed of computation and performance.

TABLE III: Traditional machine learning Versus DeepSuccinylCNN+LSTM Prediction Schemes | Sensitivity | Specificity | MCC
---|---|---|---
iSuc-PseAAC | 0.12 | 0.89 | 0.01
iSuc-PseOpt | 0.30 | 0.70 | 0.04
Psuc-Lys | 0.22 | 0.83 | 0.04
SuccinSite | 0.37 | 0.88 | 0.20
SuccinSite2.0 | 0.45 | 0.88 | 0.26
GPSuc | 0.50 | 0.88 | 0.30
PSuccE | 0.38 | 0.89 | 0.20
DeepSuccinylSiteCNN+LSTM | 0.77 | 0.63 | 0.39

##  VI Brief Overview of Convolutional Neural Networks

Convolutional neural network (CNN) was inspired by the biological processes in
the visual cortex of a rabbit. Unlike multi-layer perceptron which is fully
connected, the CNN reduces the number of training parameters by sharing
weights. CNN has been successfully for image classification [9]. It has been
successful for automatically extracting features without the need for human
intervention unlike traditional machine learning techniques. However, it still
has the challenge of pose estimation in 3D images.

CNN consists of an input layer, at least one hidden layer, and an output
layer. The input layer interfaces with the object. The hidden layers is
necessary for feature engineering. The hidden layers comprises the
convolutional layer, where a sliding-dot product operation occurs, the ReLu
layer as the activation function, the pooling layer for downsampling and
feature extraction, and the fully connected dense layer that identifies the
object. The output layer does the classification using the softmax function.

##  VII Brief Overview of Long Short Term Memory

Long Short Term Memory (LSTM) is a type of Recurrent Neural Network (RNN) deep
learning architecture useful for processing sequential data such as speech and
text [10]. Also, it is useful for processing time-series data. Unlike CNN
which is a feedforward architecture, LSTM possesses feedback architecture.

LSTM consists of a cell, an input gate, and a forget gate. The cell remembers
values over arbitrary time intervals, while the three gates control the flow
of information in and out of the cell. One advantage of LSTM over traditional
RNN is that it overcomes the challenge of vanishing gradient encountered in
traditional RNNs. The activation function for LSTM is commonly the sigmoid
function.

##  VIII Conclusion and Future Work

Succinylation is a post-translational modification process in proteins that
plays an important role in many cellular functions. Finding the succinylation
sites using experimental methods such as mass spectrometry is very laborious.
Hence, computational methods are favored. Therefore, this paper used an
improved deep learning architecture called DeepSuccinylCNN+LSTM to predict
positive succinylation sites. First, the performance of
DeepSuccinylSiteCNN+LSTM was compared with DeepSuccinylSite, the state-of-the-
art deep learning architecture for succinylation. It was shown that
DeepSuccinylCNN+LSTM performed close to DeepSuccinylSite, but with a better
trade-off in terms of computational speed. Second, the performance of
DeepSuccinylCNN+LSTM was compared with other traditional machine learning
architectures for succinylation in literature. It was shown that
DeepSuccinylCNN+LSTM performed better using some well-known performance
metrics. Hence, it can be concluded that DeepSuccinylCNN+LSTM is better
favored when large datasize is used, judging by its computational speed and
performance accuracy.

This work can be extended to finding better deep learning architecture that
will have the best performance in terms of speed and other performance
metrics. More so, novel embedding technique may be used for better
performance.

## References

  * [1] Y. Xu, Y.-X. Ding, J. Ding, Y.-H. Lei, L.-Y. Wu, and N.-Y. Deng, “isuc-pseaac: predicting lysine succinylation in proteins by incorporating peptide position-specific propensity,” _Scientific reports_ , vol. 5, p. 10184, 2015\.
  * [2] J. Jia, Z. Liu, X. Xiao, B. Liu, and K.-C. Chou, “isuc-pseopt: identifying lysine succinylation sites in proteins by incorporating sequence-coupling effects into pseudo components and optimizing imbalanced training dataset,” _Analytical biochemistry_ , vol. 497, pp. 48–56, 2016.
  * [3] ——, “psuc-lys: predict lysine succinylation sites in proteins with pseaac and ensemble random forest approach,” _Journal of theoretical biology_ , vol. 394, pp. 223–230, 2016.
  * [4] M. M. Hasan, S. Yang, Y. Zhou, and M. N. H. Mollah, “Succinsite: a computational tool for the prediction of protein succinylation sites by exploiting the amino acid patterns and properties,” _Molecular bioSystems_ , vol. 12, no. 3, pp. 786–795, 2016.
  * [5] M. M. Hasan, M. S. Khatun, M. N. H. Mollah, C. Yong, and D. Guo, “A systematic identification of species-specific protein succinylation sites using joint element features information,” _International journal of nanomedicine_ , vol. 12, p. 6303, 2017.
  * [6] M. M. Hasan and H. Kurata, “Gpsuc: Global prediction of generic and species-specific succinylation sites by aggregating multiple sequence features,” _PloS one_ , vol. 13, no. 10, p. e0200283, 2018.
  * [7] Q. Ning, X. Zhao, L. Bao, Z. Ma, and X. Zhao, “Detecting succinylation sites from protein sequences using ensemble support vector machine,” _BMC bioinformatics_ , vol. 19, no. 1, p. 237, 2018.
  * [8] N. Thapa, M. Chaudhari, S. McManus, K. Roy, R. H. Newman, H. Saigo, and D. B. Kc, “Deepsuccinylsite: a deep learning based approach for protein succinylation site prediction,” _BMC bioinformatics_ , vol. 21, pp. 1–10, 2020.
  * [9] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [10] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.


# Fractional Skyrmion molecules in a $\mathbb{C}P^{N-1}$ model

Yutaka Akagi  Department of Physics, Graduate School of Science, The
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  Yuki Amari  BLTP, JINR,
Dubna 141980, Moscow Region, Russia  Department of Mathematical Physics,
Toyama Prefectural University, Kurokawa 5180, Imizu, Toyama, 939-0398, Japan
Sven Bjarke Gudnason  Institute of Contemporary Mathematics, School of
Mathematics and Statistics, Henan University, Kaifeng, Henan 475004, P. R.
China  Muneto Nitta  Department of Physics, and Research and Education Center
for Natural Sciences, Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa
223-8521, Japan  Yakov Shnir  BLTP, JINR, Dubna 141980, Moscow Region, Russia

###### Abstract

We study fractional Skyrmions in a $\mathbb{C}P^{2}$ baby Skyrme model with a
generalization of the easy-plane potential. By numerical methods, we find
stable, metastable, and unstable solutions taking the shapes of molecules.
Various solutions possess discrete symmetries, and the origin of those
symmetries are traced back to congruencies of the fields in homogeneous
coordinates on $\mathbb{C}P^{2}$.

##  I Introduction

Solitons ubiquitously appear in nature from nonlinear media, water, condensed
matter to particle physics, dense stars, and the universe, see e.g. Refs.
Drazin and Johnson (1989); *babelon. Among solitons, topological solitons are
topologically protected stable excitations carrying topological charges, see
e.g. Refs. Manton and Sutcliffe (2004); *Shnir2018. Topological charges are
usually integer numbers since they take values in certain homotopy groups
which are isomorphic to integers, ${\mathbb{Z}}$, or to a discrete group like
${\mathbb{Z}}_{2}$.

However, fractionally quantized topological charges may also exist for a
variety of reasons. An example of a fractional topological charge was first
found in the context of vortices. The so-called Feynman–Onsager’s quantization
condition makes the circulation of superfluids or magnetic flux of
superconductors quantized integers due to the single-valuedness of the wave
function, or more mathematically due to fact that they belong to the first
homotopy group $\pi_{1}$. Nevertheless, fractionally quantized vortices are
ubiquitously allowed in diverse condensed matter systems with multiple
components, such as the A phase of superfluid 3He Salomaa and Volovik (1985);
*salomaaRMP, the uniaxially disordered superfluid 3He Autti  _et al._ (2016);
*Makinen:2018ltj, unconventional superconductors Kee  _et al._ (2000);
*Chung:2007zzc; *Garaud2012a; *Zyuzin2017; *Etter2020; *How2020,
multicomponent superconductors Babaev (2002); *Babaev:2004rm; *Babaev:2007;
*Goryo2007; *Tanaka2007; *Crisan2007; *Guikema2008; *Garaud:2011zk;
*Garaud:2012pn; *Tanaka2017; *Tanaka2018; *Chatterjee:2019jez, spinor Bose-
Einstein condensates (BECs) Leonhardt and Volovik (2000); *Ho:1998zz;
*doi:10.1143/JPSJ.67.1822; *Shinn:2018zde; *Semenoff:2006vv;
*Kobayashi:2008pk; *2010arXiv1001.2072K; *Borgh:2016cco, multicomponent BECs
Son and Stephanov (2002); *Kasamatsu2004; *Eto:2011wp; *Cipriani2013;
*Tylutki2016; *Kasamatsu2016; *MenciaUranga2018; *Eto2018; *Kobayashi:2018ezm;
*Eto:2019uhe; Eto and Nitta (2012); *Eto:2013spa; *Nitta:2013eaa,
${}^{3}P_{2}$ neutron superfluids Masuda and Nitta (2020); *Masaki:2021hmk,
exciton-polariton condensates Rubo (2007); *Lagoudakis974, and nonlinear
optics Pismen (1994); *PhysRevLett.84.634; *Pismen. In high-energy physics,
examples can be found in quantum chromodynamics at high density Balachandran
_et al._ (2006); *Nakano:2007dr; *Eto:2009kg; *Eto:2009bh; *Eto:2009tr;
*Eto:2021nle; *Fujimoto:2020dsa; Eto  _et al._ (2014), two-Higgs doublet
models which are an important and intensely studied part of physics beyond the
Standard Model of elementary particles Dvali and Senjanovic (1993);
*Eto:2018hhg; *Eto:2018tnk, and supersymmetric gauge theories Hanany and Tong
(2003); *Auzzi:2003fs; *Eto:2005yh; *Eto:2006cx; *Tong:2005un; *Eto:2006pg;
*Shifman:2007ce; *Shifman:2009zz; *Eto:2009bg; *Eto:2009bz; *Gudnason:2010yy.

Nonlinear sigma (NL$\sigma$) models admit lumps in $d=2+1$ or sigma model
instantons in $d=2+0$ (or $d=1+1$) dimensions, characterized by the second
homotopy group $\pi_{2}$ Polyakov and Belavin (1975). The gauged
$\mathbb{C}P^{1}$ model (or O(3) model) admits fractionally quantized lumps
Schroers (1995); *Schroers:1996zy; *Baptista:2004rk; *Nitta:2011um; a single
lump is decomposed into two fractional lumps which are simultaneously gauged
vortices with vortex winding numbers that sum to zero. Furthermore, the gauged
$\mathbb{C}P^{N-1}$ model, with the maximal torus action $U(1)^{N-1}$ gauged,
admits $1/N$ quantized fractional lumps Schroers (1995); *Schroers:1996zy;
*Baptista:2004rk; *Nitta:2011um, and a single lump is decomposed into $N$ of
$1/N$ fractional lumps which are simultaneously $U(1)^{N-1}$ gauged vortices
with vortex winding numbers that sum to zero. However, no explicit numerical
solutions have been constructed for $N>2$ yet and although their fractionality
is understood, the interactions between their constituents could be
complicated. Similarly, in the baby Skyrme model, which is the
$\mathbb{C}P^{1}$ model with an easy-plane potential and the baby Skyrme term,
a single baby Skyrmion splits into two fractional Skyrmions with global vortex
winding numbers that sum to zero Jaykka and Speight (2010);
*Kobayashi:2013aja; *Kobayashi:2013wra. In this case, these fractional
Skyrmions form a molecule of a certain size determined by the competition of a
long-range attraction and a short-range repulsion.111 The $\mathbb{C}P^{1}$
model can be formulated as an Abelian-Higgs model with two complex scalars, in
which lumps are replaced by semilocal vortices Vachaspati and Achucarro
(1991). An easy-plane potential, in this case, splits the semilocal vortex
into a molecule of two fractional vortices with global vortex charges that sum
to zero. This provides a stabilization mechanism of the semilocal vortex which
is unstable in the type-II region of parameter space Eto  _et al._ (2016).
Such fractional Skyrmions also appear in chiral magnets described by the
$\mathbb{C}P^{1}$ model with the Dzyaloshinskii-Moriya term and an easy-plane
potential Lin  _et al._ (2015), and are sometimes called merons. A
generalization to Skyrmions in $d=3+1$ characterized by the third homotopy
group, $\pi_{3}$, is also known; in this case, a single Skyrmion is decomposed
into two fractional Skyrmions with global monopole charges that sum to zero,
forming a molecule Gudnason and Nitta (2015).

The examples described above possess a fractionalization mechanism that is due
to the introduction of a suitable potential term. There are also other
mechanisms of fractionalization of topological charges. One mechanism is to
introduce a twisted boundary condition along a compactified direction. A
single $SU(N)$ Yang-Mills instanton on ${\mathbb{R}}^{3}\times S^{1}$ with a
twisted boundary condition along $S^{1}$ is decomposed into $N$ fractional
instantons with induced monopole charges that sum to zero Unsal (2009), while
a single $\mathbb{C}P^{N-1}$ NL$\sigma$ model instanton (Skyrmion instanton in
the $SU(N)$ principal chiral model) on ${\mathbb{R}}^{1}\times S^{1}$
(${\mathbb{R}}^{2}\times S^{1}$) with a twisted boundary condition along
$S^{1}$ is decomposed into $N$ fractional instantons with induced domain wall
(global vortex) charges that sum to zero Eto  _et al._ (2005a); *Eto:2006mz;
*Dunne:2012ae; *Dunne:2012zk; *Misumi:2014jua; *Misumi:2014bsa (Nitta (2015)).
Another mechanism is to introduce a suitable modification of the kinetic term.
The so-called modified XY model, consisting of a half($1/N$)-periodic nearest-
neighbor interaction term in addition to the usual nearest-neighbor
interaction term, admits a molecule of two ($N$) half($1/N$)-quantized
vortices Carpenter and Chalker (1989); *Kobayashi:2019sus; *Kobayashi:2019npl.

In this paper, we explore molecules of fractional Skyrmions in a
$(2+1)$-dimensional $\mathbb{C}P^{N-1}$ baby Skyrme model. We find that
fractional Skyrmions with $1/N$ topological charges appear as the minimum
constituents and they constitute molecules. In particular, when there are two
or more molecules, they form very nontrivial bound states as real molecules of
atoms. The Lagrangian density that we consider is given by

$\mathcal{L}=\frac{M^{2}}{2}{\partial}_{\mu}n^{a}{\partial}^{\mu}n^{a}-\frac{\kappa^{2}}{4}(f_{abc}n^{a}{\partial}_{\mu}n^{b}{\partial}_{\nu}n^{c})^{2}-V,$
(1)

where $M$ and $\kappa$ are real coupling constants, $V$ is a potential term,
$f_{abc}$ are the structure constants of $SU(N)$ given by
$f_{abc}=-\frac{i}{4}\mathrm{Tr}\left(\lambda_{a}\left[\lambda_{b},\lambda_{c}\right]\right)$
with $SU(N)$ generators $\lambda_{a}$ subject to the normalization
$\mathrm{Tr}(\lambda_{a}\lambda_{b})=2\delta_{ab}$. Here, the fields $n^{a}$
($a=1,2,...,N^{2}-1$) satisfy the following $N^{2}$ constraints

$n^{a}n^{a}=\frac{2(N-1)}{N},\qquad n^{a}=\frac{N}{2(N-2)}d_{abc}n^{b}n^{c},$
(2)

where
$d_{abc}=\frac{1}{4}\mathrm{Tr}\left(\lambda_{a}\\{\lambda_{b},\lambda_{c}\\}\right)$,
which imply the third-order Casimir identity,
$d_{abc}n^{a}n^{b}n^{c}=4(N-1)(N-2)/N^{2}$. The first term in Eq. (1) is the
Lagrangian density of the $\mathbb{C}P^{N-1}$ NL$\sigma$ model, also known as
the ($\mathbb{C}P^{N-1}$) Dirichlet term.

The second term of Eq. (1) is a Skyrme term introduced as a stabilizer against
shrinkage of soliton configurations, and is given by the square of the
topological current222Different types of Skyrme terms have been studied in
Refs. Ferreira and Klimas (2010); *Ferreira:2011ja; *Amari:2015sva;
*Amari:2016ynl. These Skyrme terms have a rich mathematical structure but are
somewhat intricate. Therefore, for simplicity, we will concentrate on the
Skyrme term of Eq. (1)., 333$\mathbb{C}P^{2}$ soliton configurations can also
be stabilized by the Dzyaloshinskii-Moriya type interaction Akagi  _et al._
(2021). . This term can be viewed as the $\mathbb{C}P^{N-1}$ counterpart of
the Skyrme term in the $SU(N)/U(1)^{N-1}$ Skyrme-Faddeev model Faddeev and
Niemi (1999); *Amari:2018gbq, and also a relativistic, $SU(3)$ generalization
of the so-called chiral-chiral interaction Grytsiuk  _et al._ (2020) in the
continuum limit. This term has also been used to stabilize a higher-
dimensional soliton solution in Ref. Radu  _et al._ (2013).

As for potential terms for $(2+1)$-dimensional $\mathbb{C}P^{N-1}$ models444In
$3+1$-dimensional spacetime, the so-called V-shape potential has also been
used to construct compact $Q$-solitons Klimas and Livramento (2017);
*Klimas:2018ywv; *Sawado:2020ncc; *Sawado:2021rsc., a generalization of the
easy-axis anisotropic potential, which appears for example in the multi-
layered Josephson junction Nitta (2012); *Fujimori:2016tmw, has been studied
in Refs. Abraham (1992); *Abraham:1992vb; *Gauntlett:2000ib; *Isozumi:2004jc;
*Isozumi:2004va; *Eto:2005cp; *Eto:2020cys. Here, we consider instead a
generalization of the easy-plane anisotropy term Jaykka and Speight (2010),
also known as the XY-type of potential, of the form

$V=\sum_{l=1}^{N-1}\mu^{2}\left(n^{l(l+2)}\right)^{2},$ (3)

where $\mu^{2}$ is a positive constant. This potential is a potential for the
field components corresponding to the maximal Abelian torus $U(1)^{N-1}\subset
SU(N)$. When $N=2$, Eq. (3) is just the easy-plane anisotropy term, which
generates molecules of half Skyrmions Jaykka and Speight (2010);
*Kobayashi:2013aja; *Kobayashi:2013wra. For $N=2,3$, this potential
(accompanying the Dirichlet term) appears in a system of BECs with
pseudospin-$(N-1)/2$, or in a mixture of $N$ species of bosons Graß  _et al._
(2014). We expect that for all $N>1$, the potential emerges from such a
pseudospin system. This potential is thus a simple and motivated potential
giving rise to fractional Skyrmion molecules, but modifications of such type
of potential may also possibly give rise to fractional Skyrmions.

In this paper, we numerically construct molecules of fractional Skyrmions with
$1\leq Q\leq 4$ in the $N=3$ case, where each constituent possesses a
topological charge $1/3$. For $Q=1$, we obtain a stable solution that is a
bound state of three fractional Skyrmions forming a triangle, as well as an
unstable solution that is a chain of fractional Skyrmions. For $Q=2$, we find
a stable solution whose energy density is $\mathbb{Z}_{3}$ axially symmetric,
as well as two unstable solutions composed of fractional Skyrmions that take
the shape of a hexagon and a chain, respectively. Moreover, we obtain a
metastable state which can be regarded either as a bound state of two $Q=1$
stable solutions or as a parallelogram of two $Q=1$ chains. For $Q=3$, we find
a stable solution of diamond shape, three metastable states, and two unstable
solutions. In the $Q=4$ case, we show only a stable solution and a metastable
solution, whose energy is quite close to that of the stable solution. In this
large topological charge sector, we expect that there exist many metastable
and unstable solutions, and we will not attempt to find all of them here.

This paper is organized as follows. In the next section, we introduce the
mathematical structure of the model in detail. In Sec. III, we show numerical
solutions of molecules of fractional Skyrmions and discuss their properties.
Finally, Sec. IV concludes the paper with a short discussion.

##  II The model

In this paper, we investigate static solutions of the $\mathbb{C}P^{2}$ baby
Skyrme type model (1) with the potential (3). It follows that the energy of
the model is given by $E=E_{2}+E_{4}+E_{0}$ with

$\displaystyle E_{2}=\frac{M^{2}}{2}\int
d^{2}x{\partial}_{i}n^{a}{\partial}_{i}n^{a},$ (4) $\displaystyle
E_{4}=\frac{\kappa^{2}}{4}\int
d^{2}x(f_{abc}n^{a}{\partial}_{i}n^{b}{\partial}_{j}n^{c})^{2},$ (5)
$\displaystyle E_{0}=\mu^{2}\int
d^{2}x\sum_{l=1}^{N-1}\left(n^{l(l+2)}\right)^{2},$ (6)

where the subscript $i$ of $E_{i}$ ($i=2,4,0$) indicates the number of
derivatives that the functional contains. Since the total energy includes both
terms with higher and with lower order of derivatives than the spatial
dimensions (i.e. the Skyrme term and the potential, respectively), we can
clearly evade Derrick’s no-go theorem, and Derrick’s scaling argument tells us
that the solutions should satisfy the scaling relation, $E_{4}/E_{0}=1$.

Let us discuss the (continuous) symmetries of the model. For this purpose, it
is convenient to parametrize the field as

$n^{a}=-\sqrt{\frac{N-1}{2N}}\mathrm{Tr}\left(\lambda_{a}U\lambda_{N^{2}-1}U^{\dagger}\right),$
(7)

where $U\in U(N)$. When the potential is absent, i.e. $\mu^{2}=0$, the energy
is invariant under the global $U(N)$ transformation $U\to gU$ with $g\in
U(N)$, which is equivalent to
$n^{a}\to(n^{b}/2)\mathrm{Tr}(\lambda_{a}g\lambda_{b}g^{\dagger})$. This
symmetry is broken by the potential to $U(1)^{N}$, which is the diagonal
subgroup of $U(N)$. In addition, the entire energy functional possesses the
hidden local symmetry given by $U\to Uh$ with $h\in SU(N-1)\times U(2)$
commuting with $\lambda_{N^{2}-1}$, which keeps $n^{a}$ intact.

We are in particular interested in finite energy static solutions. In order
for the static energy to be finite, the fields $n^{a}$ must decay to constant
values at spatial infinity (especially for $a=3,8,...,N^{2}-1$, $n^{a}$ must
vanish at spatial infinity). This condition effectively compactifies physical
space, $\mathbb{R}^{2}$, to the 2-sphere, i.e.
$\mathbb{R}^{2}\cup\\{\infty\\}\simeq S^{2}$, and therefore a static
configuration with finite energy corresponds to a map from $S^{2}$ to
$\mathbb{C}P^{N-1}$. Such configurations are characterized by the topological
charge

$Q=\frac{1}{8\pi}\int{\mathrm{d}}^{2}x~{}\varepsilon^{jk}f_{abc}n^{a}{\partial}_{j}n^{b}{\partial}_{k}n^{c},$
(8)

which is an element of the second homotopy group
$\pi_{2}(\mathbb{C}P^{N-1})=\mathbb{Z}\ni Q$. In this paper, we study
fractional Skyrmions, but they always appear as molecules possessing a total
topological charge being an integer.

The fields $n^{a}$ enable us to interpret the physical meaning of the
components in the Lagrangian or energy functional. However, they do not appear
to be the best option for constructing solutions because of their complicated
nonlinear constraints given in Eq. (2). Thus, we introduce the homogeneous
coordinates of $\mathbb{C}P^{N-1}$, $\\{Z^{\alpha},\bar{Z}^{\alpha}\\}$
($\alpha=1,2,...,N$) which satisfy the constraint
$\sum_{\alpha=1}^{N}|Z^{\alpha}|^{2}=1$. In terms of the homogeneous
coordinates, $n^{a}$ is given by

$n^{a}=Z^{\dagger}\lambda_{a}Z,$ (9)

with $Z=(Z^{1},Z^{2},...,Z^{N})^{\rm T}$, where T stands for the transpose
operation. It follows that $Z$ is an $N$-component complex unit vector, i.e.,
$Z^{\dagger}Z=1$. One can show that the above parametrization satisfies all
the constraints of Eq. (2) by using the identity
$(\lambda_{a})_{\alpha\beta}(\lambda_{a})_{\alpha^{\prime}\beta^{\prime}}=2(\delta_{\alpha\beta^{\prime}}\delta_{\alpha^{\prime}\beta}-\delta_{\alpha\beta}\delta_{\alpha^{\prime}\beta^{\prime}}/N)$.
This parametrization is equivalent to Eq. (7) with
$U=(Y_{1},Y_{2},...,Y_{N-1},Z)$, where $Y_{\alpha}$ are vectors forming a
complete basis in $\mathbb{C}^{N}$ along with $Z$. In terms of $Z$, the energy
and the topological charge are given by

$\displaystyle E=$ $\displaystyle\int
d^{2}x\left[2M^{2}(D_{i}Z)^{\dagger}D_{i}Z\right.$
$\displaystyle\qquad\qquad+\kappa^{2}\\{(D_{i}Z)^{\dagger}D_{j}Z-(D_{j}Z)^{\dagger}D_{i}Z\\}^{2}$
(10)
$\displaystyle\left.\qquad\qquad\qquad\quad+\frac{2\mu^{2}}{N}\sum_{\alpha<\beta}\left(|Z^{\alpha}|^{2}-|Z^{\beta}|^{2}\right)^{2}\right],$
$\displaystyle Q=$ $\displaystyle-\frac{i}{2\pi}\int
d^{2}x~{}\varepsilon^{jk}(D_{j}Z)^{\dagger}D_{k}Z,$ (11)

where $D_{i}Z={\partial}_{i}Z-(Z^{\dagger}{\partial}_{i}Z)Z$. In the next
section, we construct molecules of fractional Skyrmions by optimizing the
energy subject to the constraint $Z^{\dagger}Z=1$.

Figure 1: Toric diagram of $\mathbb{C}P^{2}$. Table 1: The numerical values of
the energy and topological charge of stable ($\circledcirc$), metastable
($\bigcirc$), and unstable ($\bigtriangleup$) solutions.

$\begin{array}[]{ccccccccc}\hline\cr Q&\text{label
(shape)}&\text{stability}&E&E_{2}&E_{4}&E_{0}&E_{4}/E_{0}&Q_{\rm Num}\\\
\hline\cr\hline\cr
1&\text{triangle}&\circledcirc&83.53&19.58&32.07&31.88&1.006&1.000\\\
&\text{chain}&\bigtriangleup&84.47&21.02&31.86&31.58&1.001&1.000\\\ \hline\cr
2&\text{curly triangle}&\circledcirc&163.11&36.09&63.59&63.43&1.002&2.000\\\
&\text{hexagon}&\bigtriangleup&163.14&36.56&63.35&63.24&1.002&2.000\\\
&\text{parallelogram}&\bigcirc&164.28&37.01&63.80&63.47&1.005&2.000\\\
&\text{chain}&\bigtriangleup&167.28&40.80&63.76&62.72&1.017&2.000\\\ \hline\cr
3&\text{3A (diamond)}&\circledcirc&243.20&51.71&95.80&95.69&1.001&3.000\\\
&\text{3B (circle)}&\bigcirc&243.59&54.21&94.73&94.65&1.000&3.000\\\
&\text{3C}&\bigcirc&243.74&52.19&95.86&95.69&1.002&3.000\\\
&\text{3D}&\bigcirc&245.05&54.02&95.68&95.35&1.003&3.000\\\
&\text{3E}&\bigtriangleup&246.41&56.55&95.03&94.83&1.002&3.000\\\ &\text{3F
(chain)}&\bigtriangleup&248.91&59.47&95.26&94.18&1.011&2.999\\\ \hline\cr
4&\text{4A}&\circledcirc&323.27&67.61&127.92&127.75&1.001&3.999\\\
&\text{4B}&\bigcirc&323.27&67.42&128.01&127.84&1.001&3.999\\\
\hline\cr\end{array}$

We introduce a tool to understand the geometrical aspects of the solutions. In
the idea of toric geometry, the manifold $\mathbb{C}P^{N-1}$ is represented by
an $(N-1)$-simplex over which an $(N-1)$-torus $T^{N-1}$ is fibered Leung and
Vafa (1998); Eto  _et al._ (2014). Each $(N-2)$-simplex composing the
$(N-1)$-simplex describes $\mathbb{C}P^{N-2}$, and therefore one of the
$T^{1}$ fibers shrink on those. At the vertices of the $(N-1)$-simplex, the
whole $T^{N-1}$ fiber shrinks. Fig. 1 is an example of the representation,
known as the toric diagram, for $\mathbb{C}P^{2}$. On the inside of the
triangle, there is a $T^{2}$ fiber corresponding to the $U(1)^{3}$ action on
the phase of the homogeneous coordinates $Z^{\alpha}\to
e^{i\theta_{\alpha}}Z^{\alpha}$ modulo the trivial $U(1)$ action, i.e.,
$\theta_{1}=\theta_{2}=\theta_{3}$. On the edges, one of the
$|Z^{\alpha}|^{2}$ vanishes, and therefore a $T^{1}$ fiber shrinks. The
fractional Skyrmions we study are interpreted as global vortices whose cores
locate where some corresponding fiber shrinks. In addition, on the toric
diagram for $\mathbb{C}P^{2}$, any straight line from a vertex to the opposite
edge represents a $\mathbb{C}P^{1}$ submanifold. If a mapping of the value
$(n^{3},n^{8})$ of a configuration onto the toric diagram describes a straight
line, the configuration is a $\mathbb{C}P^{1}$ embedding in $\mathbb{C}P^{2}$.
On the other hand, when the mapping is two-dimensional, the configuration is a
full map to $\mathbb{C}P^{2}$.

##  III Numerical solutions

We numerically construct molecules of fractional Skyrmions with $1\leq Q\leq
4$. We restrict ourselves to the $N=3$ case but conjecture that the potential
is sufficient for the model to possess fractional Skyrmion molecules for all
$N\geq 2$. We recognize $Z$ as the fundamental field and perform optimization
of the energy under the constraint $Z^{\dagger}Z=1$. Here, the optimization is
done by numerically solving the equations of motion

$\begin{split}&\frac{{\partial}{n^{s}}}{{\partial}\mathrm{Re}Z^{\alpha}}\frac{\delta
E}{\delta n^{s}}-\omega\mathrm{Re}Z^{\alpha}=0,\\\
&\frac{{\partial}{n^{s}}}{{\partial}\mathrm{Im}Z^{\alpha}}\frac{\delta
E}{\delta n^{s}}-\omega\mathrm{Im}Z^{\alpha}=0,\end{split}$ (12)

where $\omega$ is a Lagrange multiplier. The functional derivative of $E$ with
respect to $n^{s}$, $\delta E/\delta n^{s}$, is written as

$\displaystyle\frac{\delta E}{\delta n^{s}}\equiv\frac{{\partial}{\cal
E}}{{\partial}n^{s}}-{\partial}_{k}\frac{{\partial}{\cal
E}}{{\partial}({\partial}_{k}n^{s})}$
$\displaystyle=-M^{2}{\partial}_{i}^{2}n^{s}+\frac{\kappa^{2}}{2}\left(2{\partial}_{i}F_{ij}R^{s}_{j}+3F_{ij}{\partial}_{i}R^{s}_{j}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\mu^{2}(\delta_{s,3}n^{3}+\delta_{s,8}n^{8}),$
(13)

where ${\cal E}$ denotes the energy density and

$F_{ij}=f_{abc}n^{a}{\partial}_{i}n^{b}{\partial}_{j}n^{c},\quad
R^{a}_{j}=f_{abc}n^{b}{\partial}_{j}n^{c}.$ (14)

Replacing spatial derivatives by their 4th-order finite difference
approximations, we explore the solutions of Eq. (12) using a nonlinear
conjugate gradients method, with the augmented Lagrangian method.

To prepare initial configurations, we write

$Z=\frac{\left(1,u_{1},u_{2}\right)^{T}}{\sqrt{1+|u_{1}|^{2}+|u_{2}|^{2}}},$
(15)

where $u_{1}$ and $u_{2}$ are complex scalar fields, or inhomogeneous
coordinates of $\mathbb{C}P^{2}$. Finiteness of the energy requires that both
$n^{3}$ and $n^{8}$ vanish at spatial infinity, which corresponds to
$|u_{1}|^{2}$ and $|u_{2}|^{2}$ both approaching unity. Guided by this
condition, we employ initial configurations of a single lump configuration
with topological charge $Q$ of the form

$u_{1}=1+\frac{r_{1}}{(x+iy)^{Q}},\qquad u_{2}=1-\frac{r_{2}}{(x+iy)^{Q}},$
(16)

where $r_{1}$ and $r_{2}$ are nonzero real constants. In order to check the
stability of solutions obtained using the initial configuration (16) and to
construct metastable solutions that cannot be obtained when we use Eq. (16),
we also conduct simulations employing initial configurations with $Q$
separated integer Skyrmion states

$\begin{split}&u_{1}=1+\sum_{k=1}^{Q}\frac{r^{(1)}_{k}e^{i\theta^{(1)}_{k}}}{x-a_{k}+i(y-b_{k})},\\\
&u_{2}=1+\sum_{k=1}^{Q}\frac{r^{(2)}_{k}e^{i\theta^{(2)}_{k}}}{x-a_{k}+i(y-b_{k})},\end{split}$
(17)

where $r^{(j)}_{k}$ are positive constants, and $\theta^{(j)}_{k}$, $a_{k}$
and $b_{k}$ are real constants. Since the initial configurations (16) and (17)
are both holomorphic functions, they are BPS solutions of the
$\mathbb{C}P^{2}$ NL$\sigma$ model, i.e. in the $\kappa=\mu=0$ case.

Figure 2: The energy density, the norm of the homogeneous coordinates, and the
toric diagram of the $Q=1$ solutions with the coupling constants
$(M^{2},\kappa^{2},\mu^{2})=(1,16,10)$. The dots in the toric diagram
represent the values $(n^{3},n^{8})$ on each lattice point of the numerical
solution, and the color of the background shows the value of the potential.
Figure 3: The energy density, the norm of the homogeneous coordinates, and the
toric diagram of the $Q=2$ solutions. For more details, see the caption of
Fig. 2.

Figure 4: The energy density (left) and toric diagram (right) of the solutions
with $Q=3$ and $Q=4$.

We performed numerical simulations with $201^{2}$ lattice sites for $Q=1,2$
and $301^{2}$ lattice sites for $Q=3,4$, fixing the lattice spacing to
$\Delta=0.2$ and the coupling constants as
$(M^{2},\kappa^{2},\mu^{2})=(1,16,10)$. We show the energy and topological
charge of the solutions we obtain in Table 1. As shown in the table, the
numerical values of the topological charge of the solutions, $Q_{\rm Num}$,
are close to their desired integer values, and the Derrick scaling relation,
i.e. $E_{4}/E_{0}=1$, is satisfied very well too. We show the energy
distribution, the norm of the homogeneous coordinates, and a scatter plot of
$(n^{3},n^{8})$ on the toric diagram of the $Q=1$ solutions in Fig. 2. The
solutions with $Q=2$ are similarly given in Fig. 3, and for the $Q=3$ and
$Q=4$ solutions, we show only the energy densities and the toric diagrams in
Fig. 4, due to a large number of solutions. The topological charge
distribution is similar to that of the energy, so we do not present it to
avoid redundancy. We will adopt the following nomenclature for describing the
stability of the solutions, based on quantum mechanical language: stable
solutions are the global minima of the energy functional in solution space;
metastable are local minima and unstable solutions have at least one run-away
direction where the solution classically can fall into a lower-energy state.
In the following, we discuss the properties of the solutions with each
topological charge in order.

For $Q=1$, we find molecules of three fractional Skyrmions of the same size,
forming an equilateral triangle, as the stable state. In addition, unlike the
$\mathbb{C}P^{1}$ case Kobayashi and Nitta (2013) there exists an unstable
configuration, a chain of fractional Skyrmions. For both solutions, the dots
in the toric diagram in Fig. 3 scatter like a trefoil, and hence the solutions
are genuine $\mathbb{C}P^{2}$ configurations. The symmetry of the trefoils
comes from the congruency of the norm of the homogeneous coordinates: the
trefoil for the triangle solution has $\mathbb{Z}_{3}$ symmetry because all
$|Z^{\alpha}|^{2}$ are congruent to each other; on the other hand, that of the
chain is not $\mathbb{Z}_{3}$ but $\mathbb{Z}_{2}$ symmetric, because
$|Z^{1}|^{2}$ is not congruent to both $|Z^{2}|^{2}$ and $|Z^{3}|^{2}$. As
expected from the discussion of the toric diagram, one can observe in Fig. 2
that the position of Skyrmions, i.e., the energy peaks coincide with the zeros
of the homogeneous coordinates, in the same way as the positions of fractional
vortices in multicomponent BECs correspond to zeros of the condensate wave
functions Kobayashi and Nitta (2014); Eto and Nitta (2012, 2013); Nitta  _et
al._ (2013).

For $Q=2$, we obtain several molecules of six fractional Skyrmions as
(meta)stable or unstable configurations. Among them, the lowest energy
configuration is the curly triangle. The second-lowest is the hexagon, which
finally decays into the curly triangle configuration after a large number of
iterations and so is unstable. The energy density of the hexagon has literally
$\mathbb{Z}_{6}$ axial symmetry, which stems from the mirror symmetry of
$|Z^{\alpha}|^{2}$ and the congruence relation on the norm of the homogeneous
coordinates. The norm of the homogeneous coordinates of the curly triangle
type are also congruent with each other, and the congruency provides the
$\mathbb{Z}_{3}$ symmetry of its energy density. In addition to these two
configurations, we have a metastable solution that is a parallelogram state
and an unstable solution, i.e. a chain of Skyrmions. Their energy density is
$\mathbb{Z}_{2}$ axially symmetric. One can regard the parallelogram state as
a bound state of two $Q=1$ triangles or as a parallelogram of two $Q=1$
chains. On the other hand, the $Q=2$ chain is a linearly connected state of
two $Q=1$ chains. We can observe from Fig. 3 that two Skyrmions, corresponding
to zeros in the same component $|Z^{\alpha}|^{2}$, repel each other, and those
in different components attract.

In the $Q=3$ case, the stable solution is a configuration of diamond shape.
Like the $Q=2$ case, the most symmetric configuration, a molecule of nine
fractional Skyrmions sitting on a circle, has slightly higher energy than the
lowest energy state, although this state is not unstable but only metastable.
Moreover, we find two further metastable and two unstable states.

For $Q=4$, the situation becomes a bit different from the $Q=3$ case; a
molecule of fractional Skyrmions sitting on a circle does not appear as a
metastable configuration; and there appear two states, 4A (a parallelogram)
and 4B (a trapezoid), whose total energies are almost equal to one another.

As a generalization of the results shown above, it is natural to believe that
in general, the potential will allow for molecules of $3Q$ fractional
Skyrmions. In addition, we conjecture that the potential (3) equips the model
with solutions which are molecules of $NQ$ factional Skyrmions for $N>3$.

##  IV Summary and discussion

In this paper, we have studied fractional Skyrmions in a $\mathbb{C}P^{2}$
baby Skyrme model made of the standard kinetic term, the Skyrme term
generalized to $\mathbb{C}P^{2}$ and a generalization of the easy-plane
potential. The formulation of the model in terms of homogeneous coordinates is
the most convenient for numerical calculations, whereas a corresponding
inhomogeneous coordinate representation proved useful for generating initial
conditions. The solutions we have found are all molecules of $NQ$ (with $N=3$)
fractional Skyrmions with each constituent carrying $1/N$ of the topological
charge $Q$. Furthermore, all the solutions are genuine $\mathbb{C}P^{2}$
solutions, which we have shown by generating scatter plots on the toric
diagram of $\mathbb{C}P^{2}$ of the numerical solutions. For $Q=1$, the stable
solution is a triangular composition of fractional Skyrmions, and the unstable
solution is a chain that eventually will bend over and collapse into the
triangle. For $Q=2$, the stable solution has $\mathbb{Z}_{3}$ symmetry, and we
have denoted it as the curly triangle solution. A parallelogram turns out to
be a metastable state, whereas a hexagon and a chain made of two $Q=1$ chains
are both unstable states. We have also found stable, metastable, and unstable
solutions for $Q=3,4$.

While we have studied configurations of only lower $Q$ in this paper, those of
higher $Q$ remain as a future problem. In particular, revealing the structure
at sufficiently large $Q$, is the most important problem. In the case of the
$\mathbb{C}P^{1}$ model, the ground state is a square lattice for sufficiently
large $Q$ Jaykka and Speight (2010); *Kobayashi:2013aja; *Kobayashi:2013wra.
This situation also occurs in chiral magnets with an easy-plane potential Lin
_et al._ (2015). We thus expect that in our case of the $\mathbb{C}P^{2}$
model, the ground state, for sufficiently large $Q$, is a triangular lattice.

A further natural continuation of this study would be to see if the $1/N$
fractional Skyrmions persist for the $\mathbb{C}P^{N-1}$ model with our
generalization of the easy-plane potential, as we have conjectured. In
particular, it is interesting to determine whether the ground states of $Q=1$
configurations are square, pentagon, hexagon, or more generally polygons for
$N=4,5,6$ and so on. For sufficiently large $Q$, what kind of ground state is
realized is a very interesting question to explore. Since a plane can be
filled by tiling a square, triangular, or hexagonal lattice, this may
correspond to $N=2(4),3$, or $6$, respectively. It is, however, impossible for
general $N$ to fill a plane periodically. It is also interesting to study
fractional Skyrmions in the flag manifold NL$\sigma$ model Affleck  _et al._
(2021); *Ueda_2016; *Amari:2017qnb with a Skyrme term and an easy-plane type
potential. It would also be interesting to study generalizations of the
potential considered in this paper and classify which types of potentials give
rise to fractional solitons. We leave these issues for future studies.

The fractional Skyrmions studied here may very well exist also in systems
directly applicable to condensed matter physics. The most imminent and
promising situation is where the model is modified to use a generalization of
the Dzyaloshinskii-Moriya term, instead of the Skyrme term. Such a model is
expected to be realizable in a mixture of $N$ species of bosons with
artificial gauge potentials Graß  _et al._ (2014).

Another direction is to investigate whether this model admits stable Hopfions
in $d=3+1$. Since Hopfions are topologically unstable in $\mathbb{C}P^{N-1}$
models due to $\pi_{3}(\mathbb{C}P^{N-1})=0$ except for $N=2$, this is a
highly nontrivial dynamical question. If they are stable in a certain
parameter region, they would be fractional Hopfion molecules, which generalize
the fractional Hopfions in the usual Faddeev-Skyrme model with a symmetry-
breaking potential Samoilenka and Shnir (2017), possibly of the form of closed
fractional lump-strings linking each other.

### Acknowledgments

Y. Akagi thanks Hosho Katsura for a helpful discussion. The work of Y. Akagi
is supported by JSPS KAKENHI Grant No. JP20K14411 and JSPS Grant-in-Aid for
Scientific Research on Innovative Areas “Quantum Liquid Crystals” (KAKENHI
Grant No. JP20H05154). S. B. G. thanks the Outstanding Talent Program of Henan
University for partial support. The work of S. B. G. is supported by the
National Natural Science Foundation of China (Grants No. 11675223 and No.
12071111). The work of M.N. is supported in part by JSPS Grant-in-Aid for
Scientific Research (KAKENHI Grant No. 18H01217). Ya. S. gratefully
acknowledges support by the Ministry of Education of Russian Federation,
project FEWF- 2020-0003. The computations in this paper were run on the
”GOVORUN” cluster supported by the LIT, JINR.

## References

  * Drazin and Johnson (1989) P. G. Drazin and R. S. Johnson,  _Solitons: an introduction_ (Cambridge university press, 1989).
  * Babelon  _et al._ (2003) O. Babelon, D. Bernard, and M. Talon,  _Introduction to Classical Integrable Systems_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2003).
  * Manton and Sutcliffe (2004) N. S. Manton and P. M. Sutcliffe,  _Topological solitons._ , Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2004).
  * Shnir (2018) Y. M. Shnir,  _Topological and Non-Topological Solitons in Scalar Field Theories_ , Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2018).
  * Salomaa and Volovik (1985) M. M. Salomaa and G. E. Volovik, Phys. Rev. Lett. 55, 1184 (1985).
  * Salomaa and Volovik (1987) M. M. Salomaa and G. E. Volovik, Rev. Mod. Phys. 59, 533 (1987).
  * Autti  _et al._ (2016) S. Autti, V. V. Dmitriev, J. T. Mäkinen, A. A. Soldatov, G. E. Volovik, A. N. Yudin, V. V. Zavjalov, and V. B. Eltsov, Phys. Rev. Lett. 117, 255301 (2016), arXiv:1508.02197 .
  * Mäkinen  _et al._ (2019) J. T. Mäkinen, V. V. Dmitriev, J. Nissinen, J. Rysti, G. E. Volovik, A. N. Yudin, K. Zhang, and V. B. Eltsov, Nature Commun. 10, 237 (2019), arXiv:1807.04328 .
  * Kee  _et al._ (2000) H.-Y. Kee, Y. B. Kim, and K. Maki, Phys. Rev. B 62, R9275 (2000).
  * Chung  _et al._ (2007) S. B. Chung, H. Bluhm, and E.-A. Kim, Phys. Rev. Lett. 99, 197002 (2007), arXiv:0705.2660 .
  * Garaud and Babaev (2012) J. Garaud and E. Babaev, Phys. Rev. B 86, 060514 (2012).
  * Zyuzin  _et al._ (2017) A. A. Zyuzin, J. Garaud, and E. Babaev, Phys. Rev. Lett. 119, 167001 (2017), arXiv:1705.01718 .
  * Etter  _et al._ (2020) S. B. Etter, W. Huang, and M. Sigrist, New J. Phys. 22, 093038 (2020), arXiv:2006.00911 .
  * How and Yip (2020) P. T. How and S.-K. Yip, Phys. Rev. Res. 2, 043192 (2020), arXiv:2005.03366 .
  * Babaev (2002) E. Babaev, Phys. Rev. Lett. 89, 067001 (2002), cond-mat/0111192 .
  * Babaev  _et al._ (2004) E. Babaev, A. Sudbo, and N. W. Ashcroft, Nature 431, 666 (2004), cond-mat/0410408 .
  * Babaev and Ashcroft (2007) E. Babaev and N. W. Ashcroft, Nature Phys. 3, 530 (2007), arXiv:0706.2411 .
  * Goryo  _et al._ (2007) J. Goryo, S. Soma, and H. Matsukawa, Europhys. Lett. 80, 17002 (2007), cond-mat/0608015 .
  * Tanaka  _et al._ (2007) Y. Tanaka, A. Crisan, D. D. Shivagan, A. Iyo, K. Tokiwa, and T. Watanabe, Jpn. J. Appl. Phys. 46, 134 (2007).
  * Crisan  _et al._ (2007) A. Crisan, Y. Tanaka, D. D. Shivagan, A. Iyo, L. Cosereanu, K. Tokiwa, and T. Watanabe, Jpn. J. Appl. Phys. 46, L451 (2007).
  * Guikema  _et al._ (2008) J. W. Guikema, H. Bluhm, D. A. Bonn, R. Liang, W. N. Hardy, and K. A. Moler, Phys. Rev. B 77, 104515 (2008), arXiv:0802.2129 .
  * Garaud  _et al._ (2011) J. Garaud, J. Carlstrom, and E. Babaev, Phys. Rev. Lett. 107, 197001 (2011), arXiv:1107.0995 .
  * Garaud  _et al._ (2013) J. Garaud, J. Carlström, E. Babaev, and M. Speight, Phys. Rev. B 87, 014507 (2013), arXiv:1211.4342 .
  * Tanaka  _et al._ (2017) Y. Tanaka, H. Yamamori, T. Yanagisawa, T. Nishio, and S. Arisawa, Phys. C Supercond. its Appl. 538, 12 (2017).
  * Tanaka  _et al._ (2018) Y. Tanaka, H. Yamamori, T. Yanagisawa, T. Nishio, and S. Arisawa, Phys. C Supercond. its Appl. 548, 44 (2018).
  * Chatterjee  _et al._ (2020) C. Chatterjee, S. B. Gudnason, and M. Nitta, JHEP 04, 109 (2020), arXiv:1912.02685 .
  * Leonhardt and Volovik (2000) U. Leonhardt and G. E. Volovik, Pisma Zh. Eksp. Teor. Fiz. 72, 66 (2000), cond-mat/0003428 .
  * Ho (1998) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998), cond-mat/9803231 .
  * Ohmi and Machida (1998) T. Ohmi and K. Machida, J. Phys. Soc. Jpn. 67, 1822 (1998), cond-mat/9803160 .
  * Shinn and Fischer (2018) S.-H. Shinn and U. R. Fischer, Phys. Rev. A 98, 053602 (2018), arXiv:1805.05623 .
  * Semenoff and Zhou (2007) G. W. Semenoff and F. Zhou, Phys. Rev. Lett. 98, 100401 (2007), cond-mat/0610162 .
  * Kobayashi  _et al._ (2009) M. Kobayashi, Y. Kawaguchi, M. Nitta, and M. Ueda, Phys. Rev. Lett. 103, 115301 (2009), arXiv:0810.5441 .
  * Kawaguchi and Ueda (2012) Y. Kawaguchi and M. Ueda, Phys. Rep. 520, 253 (2012), arXiv:1001.2072 .
  * Borgh and Ruostekoski (2016) M. O. Borgh and J. Ruostekoski, Phys. Rev. Lett. 117, 275302 (2016), [Erratum: Phys.Rev.Lett. 118, 129901 (2017)], arXiv:1611.09735 .
  * Son and Stephanov (2002) D. T. Son and M. A. Stephanov, Phys. Rev. A 65, 063621 (2002), cond-mat/0103451 .
  * Kasamatsu  _et al._ (2004) K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. Lett. 93, 250406 (2004), cond-mat/0406150 .
  * Eto  _et al._ (2011) M. Eto, K. Kasamatsu, M. Nitta, H. Takeuchi, and M. Tsubota, Phys. Rev. A 83, 063603 (2011), arXiv:1103.6144 .
  * Cipriani and Nitta (2013) M. Cipriani and M. Nitta, Phys. Rev. Lett. 111, 170401 (2013), arXiv:1303.2592 .
  * Tylutki  _et al._ (2016) M. Tylutki, L. P. Pitaevskii, A. Recati, and S. Stringari, Phys. Rev. A 93, 043623 (2016), arXiv:1601.03695 .
  * Kasamatsu  _et al._ (2016) K. Kasamatsu, M. Eto, and M. Nitta, Phys. Rev. A 93, 013615 (2016), arXiv:1510.00139 .
  * Mencia Uranga and Lamacraft (2018) B. Mencia Uranga and A. Lamacraft, Phys. Rev. A 97, 043609 (2018), arXiv:1801.09912 .
  * Eto and Nitta (2018) M. Eto and M. Nitta, Phys. Rev. A 97, 023613 (2018), arXiv:1702.04892 .
  * Kobayashi  _et al._ (2019) M. Kobayashi, M. Eto, and M. Nitta, Phys. Rev. Lett. 123, 075303 (2019), arXiv:1802.08763 .
  * Eto  _et al._ (2020a) M. Eto, K. Ikeno, and M. Nitta, Phys. Rev. Res. 2, 033373 (2020a), arXiv:1912.09014 .
  * Eto and Nitta (2012) M. Eto and M. Nitta, Phys. Rev. A 85, 053645 (2012), arXiv:1201.0343 .
  * Eto and Nitta (2013) M. Eto and M. Nitta, Europhys. Lett. 103, 60006 (2013), arXiv:1303.6048 .
  * Nitta  _et al._ (2013) M. Nitta, M. Eto, and M. Cipriani, J. Low Temp. Phys. 175, 177 (2013), arXiv:1307.4312 .
  * Masuda and Nitta (2020) K. Masuda and M. Nitta, PTEP 2020, 013D01 (2020), arXiv:1602.07050 .
  * Masaki  _et al._ (2021) Y. Masaki, T. Mizushima, and M. Nitta, (2021), arXiv:2107.02448 .
  * Rubo (2007) Y. G. Rubo, Phys. Rev. Lett. 99, 106401 (2007), arXiv:0704.2721 .
  * Lagoudakis  _et al._ (2009) K. G. Lagoudakis, T. Ostatnický, A. V. Kavokin, Y. G. Rubo, R. André, and B. Deveaud-Plédran, Science 326, 974 (2009), arXiv:0911.2555 .
  * Pismen (1994) L. M. Pismen, Phys. Rev. Lett. 72, 2557 (1994).
  * Aranson and Pismen (2000) I. S. Aranson and L. M. Pismen, Phys. Rev. Lett. 84, 634 (2000), cond-mat/9909108 .
  * Pismen (1999) L. Pismen,  _Vortices in Nonlinear Fields: From Liquid Crystals to Superfluids, from Non-Equilibrium Patterns to Cosmic Strings_ , International Series of Monographs on Physics (Clarendon Press, 1999).
  * Balachandran  _et al._ (2006) A. P. Balachandran, S. Digal, and T. Matsuura, Phys. Rev. D 73, 074009 (2006), hep-ph/0509276 .
  * Nakano  _et al._ (2008) E. Nakano, M. Nitta, and T. Matsuura, Phys. Rev. D 78, 045002 (2008), arXiv:0708.4096 .
  * Eto and Nitta (2009) M. Eto and M. Nitta, Phys. Rev. D 80, 125007 (2009), arXiv:0907.1278 .
  * Eto  _et al._ (2009a) M. Eto, E. Nakano, and M. Nitta, Phys. Rev. D 80, 125011 (2009a), arXiv:0908.4470 .
  * Eto  _et al._ (2010) M. Eto, M. Nitta, and N. Yamamoto, Phys. Rev. Lett. 104, 161601 (2010), arXiv:0912.1352 .
  * Eto and Nitta (2021) M. Eto and M. Nitta, (2021), arXiv:2103.13011 .
  * Fujimoto and Nitta (2021) Y. Fujimoto and M. Nitta, Phys. Rev. D 103, 054002 (2021), arXiv:2011.09947 .
  * Eto  _et al._ (2014) M. Eto, Y. Hirono, M. Nitta, and S. Yasui, PTEP 2014, 012D01 (2014), arXiv:1308.1535 .
  * Dvali and Senjanovic (1993) G. R. Dvali and G. Senjanovic, Phys. Rev. Lett. 71, 2376 (1993), hep-ph/9305278 .
  * Eto  _et al._ (2018a) M. Eto, M. Kurachi, and M. Nitta, Phys. Lett. B 785, 447 (2018a), arXiv:1803.04662 .
  * Eto  _et al._ (2018b) M. Eto, M. Kurachi, and M. Nitta, JHEP 08, 195 (2018b), arXiv:1805.07015 .
  * Hanany and Tong (2003) A. Hanany and D. Tong, JHEP 07, 037 (2003), hep-th/0306150 .
  * Auzzi  _et al._ (2003) R. Auzzi, S. Bolognesi, J. Evslin, K. Konishi, and A. Yung, Nucl. Phys. B 673, 187 (2003), hep-th/0307287 .
  * Eto  _et al._ (2006a) M. Eto, Y. Isozumi, M. Nitta, K. Ohashi, and N. Sakai, Phys. Rev. Lett. 96, 161601 (2006a), hep-th/0511088 .
  * Eto  _et al._ (2006b) M. Eto, K. Konishi, G. Marmorini, M. Nitta, K. Ohashi, W. Vinci, and N. Yokoi, Phys. Rev. D 74, 065021 (2006b), hep-th/0607070 .
  * Tong (2005) D. Tong, in  _Theoretical Advanced Study Institute in Elementary Particle Physics: Many Dimensions of String Theory_ (2005) hep-th/0509216 .
  * Eto  _et al._ (2006c) M. Eto, Y. Isozumi, M. Nitta, K. Ohashi, and N. Sakai, J. Phys. A 39, R315 (2006c), hep-th/0602170 .
  * Shifman and Yung (2007) M. Shifman and A. Yung, Rev. Mod. Phys. 79, 1139 (2007), hep-th/0703267 .
  * Shifman and Yung (2009) M. Shifman and A. Yung,  _Supersymmetric solitons_ , Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2009).
  * Eto  _et al._ (2009b) M. Eto, T. Fujimori, S. B. Gudnason, K. Konishi, T. Nagashima, M. Nitta, K. Ohashi, and W. Vinci, JHEP 06, 004 (2009b), arXiv:0903.4471 .
  * Eto  _et al._ (2009c) M. Eto, T. Fujimori, S. B. Gudnason, K. Konishi, T. Nagashima, M. Nitta, K. Ohashi, and W. Vinci, Phys. Rev. D 80, 045018 (2009c), arXiv:0905.3540 .
  * Gudnason (2010) S. B. Gudnason, Nucl. Phys. B 840, 160 (2010), arXiv:1005.0557 .
  * Polyakov and Belavin (1975) A. M. Polyakov and A. A. Belavin, JETP Lett. 22, 245 (1975).
  * Schroers (1995) B. J. Schroers, Phys. Lett. B 356, 291 (1995), hep-th/9506004 .
  * Schroers (1996) B. J. Schroers, Nucl. Phys. B 475, 440 (1996), hep-th/9603101 .
  * Baptista (2006) J. M. Baptista, Commun. Math. Phys. 261, 161 (2006), math/0411517 .
  * Nitta and Vinci (2012) M. Nitta and W. Vinci, J. Phys. A 45, 175401 (2012), arXiv:1108.5742 .
  * Jaykka and Speight (2010) J. Jaykka and M. Speight, Phys. Rev. D 82, 125030 (2010), arXiv:1010.2217 .
  * Kobayashi and Nitta (2013) M. Kobayashi and M. Nitta, Phys. Rev. D 87, 125013 (2013), arXiv:1307.0242 .
  * Kobayashi and Nitta (2014) M. Kobayashi and M. Nitta, J. Low Temp. Phys. 175, 208 (2014), arXiv:1307.1345 .
  * Vachaspati and Achucarro (1991) T. Vachaspati and A. Achucarro, Phys. Rev. D 44, 3067 (1991).
  * Eto  _et al._ (2016) M. Eto, M. Nitta, and K. Sakurai, JHEP 10, 048 (2016), arXiv:1608.03516 .
  * Lin  _et al._ (2015) S.-Z. Lin, A. Saxena, and C. D. Batista, Phys. Rev. B 91, 224407 (2015), arXiv:1406.1422 .
  * Gudnason and Nitta (2015) S. B. Gudnason and M. Nitta, Phys. Rev. D 91, 085040 (2015), arXiv:1502.06596 .
  * Unsal (2009) M. Unsal, Phys. Rev. D 80, 065001 (2009), arXiv:0709.3269 .
  * Eto  _et al._ (2005a) M. Eto, Y. Isozumi, M. Nitta, K. Ohashi, and N. Sakai, Phys. Rev. D 72, 025011 (2005a), hep-th/0412048 .
  * Eto  _et al._ (2006d) M. Eto, T. Fujimori, Y. Isozumi, M. Nitta, K. Ohashi, K. Ohta, and N. Sakai, Phys. Rev. D 73, 085008 (2006d), hep-th/0601181 .
  * Dunne and Unsal (2012) G. V. Dunne and M. Unsal, JHEP 11, 170 (2012), arXiv:1210.2423 .
  * Dunne and Ünsal (2013) G. V. Dunne and M. Ünsal, Phys. Rev. D 87, 025015 (2013), arXiv:1210.3646 .
  * Misumi  _et al._ (2014) T. Misumi, M. Nitta, and N. Sakai, JHEP 06, 164 (2014), arXiv:1404.7225 .
  * Misumi  _et al._ (2015) T. Misumi, M. Nitta, and N. Sakai, PTEP 2015, 033B02 (2015), arXiv:1409.3444 .
  * Nitta (2015) M. Nitta, JHEP 08, 063 (2015), arXiv:1503.06336 .
  * Carpenter and Chalker (1989) D. B. Carpenter and J. T. Chalker, Journal of Physics: Condensed Matter 1, 4907 (1989).
  * Kobayashi  _et al._ (2020) M. Kobayashi, G. Fejős, C. Chatterjee, and M. Nitta, Phys. Rev. Res. 2, 013081 (2020), arXiv:1908.11087 .
  * Kobayashi and Nitta (2020) M. Kobayashi and M. Nitta, Phys. Rev. D 101, 085003 (2020), arXiv:1912.09456 .
  * Ferreira and Klimas (2010) L. A. Ferreira and P. Klimas, JHEP 10, 008 (2010), arXiv:1007.1667 .
  * Ferreira  _et al._ (2011) L. A. Ferreira, P. Klimas, and W. J. Zakrzewski, JHEP 12, 098 (2011), arXiv:1111.2338 .
  * Amari  _et al._ (2015) Y. Amari, P. Klimas, N. Sawado, and Y. Tamaki, Phys. Rev. D 92, 045007 (2015), arXiv:1504.02848 .
  * Amari  _et al._ (2016) Y. Amari, P. Klimas, and N. Sawado, Phys. Rev. D 94, 025032 (2016), arXiv:1604.06125 .
  * Akagi  _et al._ (2021) Y. Akagi, Y. Amari, N. Sawado, and Y. Shnir, Phys. Rev. D 103, 065008 (2021), arXiv:2101.10566 .
  * Faddeev and Niemi (1999) L. D. Faddeev and A. J. Niemi, Phys. Lett. B 449, 214 (1999), hep-th/9812090 .
  * Amari and Sawado (2018a) Y. Amari and N. Sawado, Phys. Lett. B 784, 294 (2018a), arXiv:1805.10008 .
  * Grytsiuk  _et al._ (2020) S. Grytsiuk, J.-P. Hanke, M. Hoffmann, J. Bouaziz, O. Gomonay, G. Bihlmayer, S. Lounis, Y. Mokrousov, and S. Blügel, Nat. Commun. 11, 1 (2020), arXiv:1904.02369 .
  * Radu  _et al._ (2013) E. Radu, D. H. Tchrakian, and Y. Yang, Nucl. Phys. B 875, 388 (2013), arXiv:1305.4784 .
  * Klimas and Livramento (2017) P. Klimas and L. R. Livramento, Phys. Rev. D 96, 016001 (2017), arXiv:1704.01132 .
  * Klimas  _et al._ (2019) P. Klimas, N. Sawado, and S. Yanai, Phys. Rev. D 99, 045015 (2019), arXiv:1812.08363 .
  * Sawado and Yanai (2020) N. Sawado and S. Yanai, Phys. Rev. D 102, 045007 (2020), arXiv:2006.03424 .
  * Sawado and Yanai (2021) N. Sawado and S. Yanai, Phys. Rev. D 103, 125018 (2021), arXiv:2103.05877 .
  * Nitta (2012) M. Nitta, Phys. Rev. D 86, 125004 (2012), arXiv:1207.6958 .
  * Fujimori  _et al._ (2016) T. Fujimori, H. Iida, and M. Nitta, Phys. Rev. B 94, 104504 (2016), arXiv:1604.08103 .
  * Abraham (1992) E. Abraham, Phys. Lett. B 278, 291 (1992).
  * Abraham and Townsend (1992) E. R. C. Abraham and P. K. Townsend, Phys. Lett. B 291, 85 (1992).
  * Gauntlett  _et al._ (2001) J. P. Gauntlett, D. Tong, and P. K. Townsend, Phys. Rev. D 64, 025010 (2001), hep-th/0012178 .
  * Isozumi  _et al._ (2004a) Y. Isozumi, M. Nitta, K. Ohashi, and N. Sakai, Phys. Rev. Lett. 93, 161601 (2004a), hep-th/0404198 .
  * Isozumi  _et al._ (2004b) Y. Isozumi, M. Nitta, K. Ohashi, and N. Sakai, Phys. Rev. D 70, 125014 (2004b), hep-th/0405194 .
  * Eto  _et al._ (2005b) M. Eto, Y. Isozumi, M. Nitta, K. Ohashi, and N. Sakai, Phys. Rev. D 72, 085004 (2005b), hep-th/0506135 .
  * Eto  _et al._ (2020b) M. Eto, M. Kawaguchi, M. Nitta, and R. Sasaki, Phys. Rev. D 101, 105020 (2020b), arXiv:2003.13520 .
  * Graß  _et al._ (2014) T. Graß, R. W. Chhajlany, C. A. Muschik, and M. Lewenstein, Phys. Rev. B 90, 195127 (2014), arXiv:1408.0769 .
  * Leung and Vafa (1998) N. C. Leung and C. Vafa, Adv. Theor. Math. Phys. 2, 91 (1998), hep-th/9711013 .
  * Affleck  _et al._ (2021) I. Affleck, D. Bykov, and K. Wamer, (2021), arXiv:2101.11638 .
  * Ueda  _et al._ (2016) H. T. Ueda, Y. Akagi, and N. Shannon, Phys. Rev. A 93, 021606(R) (2016), arXiv:1511.06515 .
  * Amari and Sawado (2018b) Y. Amari and N. Sawado, Phys. Rev. D 97, 065012 (2018b), arXiv:1711.00933 .
  * Samoilenka and Shnir (2017) A. Samoilenka and Y. Shnir, JHEP 09, 029 (2017), arXiv:1707.06608 .


On the Required Secure Key Rate for Quantum-Secured Optical Channels

Farzam Toudeh-Fallah1 , Robert Keys1 , Dave Atkinson1

1 Quantum Communication and Photonic Systems, Ciena Corporation

###### Abstract

The current maturity of the quantum-secured optical data channels based on the
Quantum Key Distribution (QKD) technology is at the deployment level in metro
environments, while R&D efforts are also being conducted towards long-distance
deployments. A great deal of research has been conducted on the achievable
Secure Key Rate (SKR) for quantum channels. However, one of the major
questions for network operators is the required SKR for the deployment of
quantum-secured channels. This article addresses this question by defining the
required SKR for quantum-secured optical channels and provides guidelines
towards optimizing this parameter.

##  1 Introduction

The current key distribution methods utilize algorithms to encrypt the secret
symmetric keys to be distributed amongst the communicating parties. These keys
are then used to encrypt and decrypt the messages exchanged between the
parties on the communication channels. However, these algorithm-based
approaches will be vulnerable to attacks by quantum computers [1, 2]. In
contrast, Quantum Key Distribution (QKD) relies on the principles of quantum
mechanics to distribute quantum-secured secret keys, backed by theoretical
proof of unconditional security[3, 4, 5]. Consequently, quantum-secured data
channels based on the QKD technology have attracted a global interest.
Instituting such channels over long-distances faces technological challenges
and is the subject of world-wide research efforts. However, the maturity level
of the short-distance quantum-secured channels is at a deployment level in
metro environments[6].

Secure Key Rate (SKR) (bits per second) is a measure of the key generation
rate by QKD systems. Although several studies have been conducted on the
achievable SKR for QKD systems [7, 8, 9, 10], a major question raised by the
network operators is the required SKR and how it relates to their specific
operational environments. This article addresses this question with the intent
of fostering a pragmatic approach towards defining this requirement for such
deployments.

##  2 Quantum-Secured Optical Channels

Figure 1 depicts a generic model for establishing quantum-secured optical data
channels.

Fig. 1: A generic model representing quantum-secured optical channels

As shown in this figure, QKD technology utilizes two channels (quantum and
classical) to generate symmetric secret keys in both communication sites.
These keys are then utilized to encrypt the optical data channels provisioned
between the two sites using a symmetric encryption methodology such as
AES-256, hence establishing quantum-secured optical channels. On some QKD
systems, the quantum channel is established on the C band, in which case
mostly a dedicated fiber will be required to prevent any interference, such as
Raman Scattering, from the optical channels. However, some QKD systems
establish the quantum channel on the O band, which allows the quantum channel
to co-exist on the same fiber as the optical data channels. This study is
independent of how quantum channel is established, its wavelength, encoding,
protocol, or any other attributes.

##  3 Defining the Required Secure Key Rate

In the process of deploying quantum-secured optical channels, network
operators must first establish the Required Secure Key Rate (RSKR). In high-
capacity operational links, the entire C band, or in some cases, the entire
C+L bands are filled with DWDM optical data channels and a deployment model
must be able to quantum-secure the entire link. Therefore, the question arises
on the Required SKR to satisfy this model. This requirement is dictated by the
contributing factors on the optical data channels, as described in this
section.

First, we define the Required Secure Key Rate Per Data Channel, which is the
minimum key generation rate required for quantum-securing a data channel,
designated as $RSKR_{DC}$ (bits per second). This requirement depends on the
secure key length (bits) and how fast the key is refreshed (per second) on
each channel, as defined in equation 1.

$RSKR_{DC}=L_{k}R_{k}$ (1)

where:
$RSKR_{DC}$: Required Secure Key Rate (RSKR) per data channel (bits/sec)
$L_{k}$: Required key length associated with data channel k (bits)
$R_{k}$: Required key refresh rate associated with data channel k (per second)

A link usually consists of several multiplexed data channels (such as DWDM
fiber links) and in order to calculate the Required Secure Key Rate for the
entire link, the total number of data channels multiplexed on the link has to
be taken into account as follows:

$RSKR_{L}=\sum_{k=1}^{N_{DC}}L_{k}R_{k}$ (2)

where:
$RSKR_{L}$: Required Secure Key Rate (RSKR) per link (bits/sec)
$N_{DC}$: Total number of channels per link

Equation 2 is a general form applied to any link format with different key
length and key refresh rate per channel, fix-grid or flex-grid, or any other
format and provides the minimum SKR required for quantum-securing the entire
operational band on a link. In the case of all channels on a link utilizing
the same key length and key refresh rate, equation 2 turns into the following
form:

$RSKR_{L}=L_{k}R_{k}N_{DC}$ (3)

Equation 3 clearly shows the dependency of the RSKR per link not only on the
key length and the key refresh rate, but also on the number of data channels
per link. Therefore, by reducing the number of data channels on a link, one
can reduce the RSKR per link accordingly. Equation 3 could be applied to both
fix-grid and flex-grid formatted links. However, to study an important
consequence of this equation, let us consider the case of fix-grid DWDM-
formatted links. In this case, equation 3 will result in:

$RSKR_{L}=L_{K}R_{K}\biggl{\lfloor}\frac{\delta_{B}}{\delta_{C}}\biggr{\rfloor}$
(4)

where:
$\delta_{B}$: Total spectrum of the operational band per link (Hz)
$\delta_{C}$: Spectral occupancy per data channel (Hz)

Note that in equation 4 $\lfloor\hskip 2.84544pt\rfloor$ denotes the floor
function.

##  4 Impact of Optical Data Channel Configuration on the Link RSKR

According to equation 4, a shift to higher spectral occupancy per data channel
$(\delta_{C})$ would result into a drop in the RSKR per link. An interesting
case study is the increase in spectral occupancy towards achieving higher
channel capacity. The relation between the channel capacity and spectral
occupancy is a function of several factors, such as signaling rate and
encoding mechanism. Based on equation 4, Figure 2 shows the dependency between
the RSKR per link and the channel capacity for C band-only and C+L band links.

Fig. 2: Required SKR Per Link vs Channel Capacity (at minimum spectral
occupancy)

The channel capacity in Figure 2 is based on the minimum spectral occupancy
and represents the experimental data on some of the current deployed optical
data channels up to 800 Gbps and is extrapolated to higher channel capacities
as specified in Table 1.

Channel Capacity (Gbps) |  Minimum Spectral Occupancy (GHz)
---|---
100 |  50
400 |  75
800 |  112.5
1600 |  225
Table 1: Channel Capacity and the Corresponding Minimum Spectral Occupancy

In addition, the following spectrum ranges were considered in Figure 2:

C band: 191.35 - 196.10 THz

C+L band: 184.49 - 196.10 THz

It should be mentioned that Figure 2 is based on: 1) utilizing AES-256 as the
encryption method for the optical data channels (at present the most widely
used method in operational networks) and 2) a key refresh rate of one per
second.

The analysis presented in Figure 2 reveals an important observation regarding
the Required SKR per link. In both cases (C-band and C+L band), a decline in
the Required SKR per link is observed as a function of increase in the optical
data channel capacity (as it increases the spectral occupancy per channel).
Therefore, according to this analysis, data aggregation to produce channels
with higher throughput reduces the Require SKR per link and could be utilized
as an optimization factor towards deployment of quantum-secured optical data
channels. This correlation is more pronounced in the C+L band scenario.

It should be noted that the above-mentioned analysis is independent of the
quantum channel used to generate the QKD keys. For instance, in addition to
terrestrial quantum channels (established on any wavelength), this analysis is
also applicable to the case where optical data channels are deployed on
terrestrial fiber links, while the quantum channel is satellite-based (free-
space in general).

Note that in addition to channel aggregation to higher capacity data channels,
any other requirement on the optical data channels that causes an increase in
the channel spectral occupancy will result in a decline in the Required SKR
per link. Increasing the signaling rate towards achieving longer distance for
a fixed channel capacity is another example.

##  5 Conclusion

In this article the Required Secure Key Rate (RKSR) was defined at both data
channel and link levels for the deployment of quantum-secured optical
channels. This is one of the fundamental requirements to be considered for
such deployments. Upon defining the RSKR, the optical data channel spectral
occupancy was identified as one of the most important contributing factors in
determining this requirement. Accordingly, an RSKR optimization methodology
based on data channel aggregation was discussed.

##  6 Acknowledgements

The authors would like to thank Charles Laperle and David Fasken for their
review of this article and their comments.

## References

  * [1] M. Mosca and M. Piani, ”Quantum Threat Timeline Report 2022”, Global Risk Institute, December 2022
  * [2] C. Gidney and M. Ekera, ”How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits”, Arxiv:1905.09749v3, 2021
  * [3] P. Shor and J. Preskill, ”Simple Proof of Security of the BB84 Quantum Key Distribution Protocol”, Phys Rev Lett 85, 441, 2000
  * [4] H.-K. Lo and H. Chau, ”Unconditional Security of Quantum Key Distribution Over Arbitrarily Long Distances”, Science, vol. 283, no. 5410, 1999
  * [5] D. Mayers, ”Unconditional Security in Quantum Cryptography”, Journal of the ACM, vol. 48, no. 3, 2001
  * [6] P. Sandle, ”BT and Toshiba trial first commercial quantum-secured network”, Reuters, 26 April 2022. [Online]. Available: https://www.reuters.com/business/bt-toshiba-trial-first-commercial-quantum-secured-network-2022-04-26/
  * [7] A. Shields, ”Performance Limits for Quantum Key Distribution Networks”, ITU-T Workshop, Shanghai, 2019
  * [8] M. Takeoka, S. Guha and M. Wilde, ”Fundamental rate-loss tradeoff for optical quantum key distribution”, Nature Communications 5, 5235, 2014.
  * [9] S. Pirandola, R. Laurenza, C. Ottaviani and L. Banchi, ”Fundamental limits of repeaterless quantum communications”, Nat Commun 8, 15043, 2017.
  * [10] X. Wang, S. Guo, P. Wang, W. Liu and Y. Li, ”Realistic rate–distance limit of continuous-variable quantum key distribution”, Optics Express, vol. 27, no. 9, 2019.