# Secondary frequency control stabilizing the voltage

Eder Batista Tchawou Tchuisseu  Institute of Thermomechanics, Academy of
Science of the Czech Republic, 18200 Prague 8, Czech Republic
<EMAIL_ADDRESS>Eric-Donald Dongmo  Department of Mechanical Engineering,
College of Technology, University of Buea, Po. Box 63, Buea, Cameroon
Laboratory of Modeling and Simulation in Bio-Engineering and Prototypes,
University of Yaounde 1  Pavel Procházka  Institute of Thermomechanics,
Academy of Science of the Czech Republic, 18200 Prague 8, Czech Republic
<EMAIL_ADDRESS>Paul Woafo  Laboratory of Modeling and Simulation in Bio-
Engineering and Prototypes, University of Yaounde 1  Pere Colet  Instituto de
Física Interdisciplinar y Sistemas Complejos, IFISC (CSIC-UIB), Campus
Universitat Illes Balears, E-07122 Palma de Mallorca, Spain  Benjamin Schäfer
<EMAIL_ADDRESS>School of Mathematical Sciences, Queen Mary
University of London, London E1 4NS, United Kingdom  Faculty of Science and
Technology, Norwegian University of Life Sciences, 1432 Ås, Norway

# Secondary frequency control stabilizing voltage dynamics

Eder Batista Tchawou Tchuisseu  Institute of Thermomechanics, Academy of
Science of the Czech Republic, 18200 Prague 8, Czech Republic
<EMAIL_ADDRESS>Eric-Donald Dongmo  Department of Mechanical Engineering,
College of Technology, University of Buea, Po. Box 63, Buea, Cameroon
Laboratory of Modeling and Simulation in Bio-Engineering and Prototypes,
University of Yaounde 1  Pavel Procházka  Institute of Thermomechanics,
Academy of Science of the Czech Republic, 18200 Prague 8, Czech Republic
<EMAIL_ADDRESS>Paul Woafo  Laboratory of Modeling and Simulation in Bio-
Engineering and Prototypes, University of Yaounde 1  Pere Colet  Instituto de
Física Interdisciplinar y Sistemas Complejos, IFISC (CSIC-UIB), Campus
Universitat Illes Balears, E-07122 Palma de Mallorca, Spain  Benjamin Schäfer
<EMAIL_ADDRESS>School of Mathematical Sciences, Queen Mary
University of London, London E1 4NS, United Kingdom  Faculty of Science and
Technology, Norwegian University of Life Sciences, 1432 Ås, Norway

###### Abstract

The ongoing energy transition challenges the stability of the electrical power
system. Stable operation of the electrical power grid requires both the
voltage (amplitude) and the frequency to stay within operational bounds. While
much research has focused on frequency dynamics and stability, the voltage
dynamics has been neglected. Here, we study frequency and voltage stability in
the case of simple networks via linear stability and bulk analysis. In
particular, our linear stability analysis of the network shows that the
frequency secondary control guarantees the stability of a particular electric
network. Even more interesting, while we only consider secondary frequency
control, we observe a stabilizing effect on the voltage dynamics, especially
in our numerical bulk analysis.

##  I Introduction

The need for good and stable electricity is a current and urgent quest in our
society [1, 2]. Electricity is generated by the conversion of a primary source
of energy such are mechanical, chemical, nuclear or thermal to electrical
energy. To power small devices, energy harvesting systems are commonly used
[3, 4], while for powering cities or countries synchronous generators and
renewable energy sources are needed [5, 6, 7]. The infrastructure connecting
such generators and consumers of electricity is called the electric power
grid. The traditional power system stability remains an important task to be
achieved by electric utilities in order to ensure a good electric quality and
supply security for the consumers [5, 6, 8, 9, 10, 11, 12]. Power imbalances
are one of the principal causes of grid instabilities, which can lead the
network to blackouts. In fact, any power imbalance induces a variation of the
frequency and the voltage (amplitude) of the electric grid.

In the literature, much work has been devoted to propose controllers that
stabilize the grid when facing any power imbalance. Many such controllers
mainly focused on either controlling the frequency [5, 11, 7, 13] of the grid
or the voltage [14] of the power grid. For the frequency control, the governor
of the power plant is often used through the load frequency control and the
automatic generation control, mostly known as the primary, secondary and
tertiary frequency controls [7, 15, 16, 17]. The voltage on the other hand is
controlled through the automatic voltage regulator, which ensures that the
voltage is kept within an admissible range. Therefore, many studies only
focused on either the control of the frequency or the control of the voltage
and rarely both.

There are three leading models [18] to mathematically describe the power grid,
which are: the effective network; the structure-preserving model and the
synchronous motor model. In this paper, we are using the synchronous motor
model, where each node in the network ("the motor") can be considered as an
aggregate of generators or consumers (e.g. a small region or large city). The
voltage sources of these machines are usually considered constant, such that
the dynamic of the power system is reduced to its frequency, hence phase
dynamics [18, 5, 19, 20].

Within this article, we investigate the stability of the high-voltage
transmission system while including the voltage dynamics. This model has been
studied in the electric network research community e.g. by Katrin
Schmietendorf et al. in [21], Sabine Auer et al. [22] and Florian Döfler et
al. [23], all stressing the need to include voltage dynamics even on
transmission level. However, these works do not consider any type of control
applied to the variables of the network (phase, frequencies and voltages). But
this control is essential for the global understanding of the power system
dynamics and stability. For example, it has been shown that applying secondary
frequency control to the well-known second order electric network model, can
deeply modify the dynamics of a network [1, 15, 5]. Thus, we aim to couple the
frequency secondary control to the electric network considering the voltage
dynamics. The results obtained including the voltage dynamics will be compared
with the ones obtained with the classical uncontrolled model.

The rest of this paper is organized as follows: Section 2 provides a
mathematical model of the power grid considered in this work as networks of
synchronous machines controlled each by frequency controllers. Based on linear
stability analysis and bulk dynamics, we also quantify the stability of the
power grid. Numerical analyses are presented in Section 3 for small (N=2) and
larger (N>2) networks. Finally, the paper is concluded in Section 4.

##  II Mathematical model and stability analysis

###  II.1 Mathematical model

The electric network is modeled as coupled synchronous machines combined with
their respective frequency controllers (secondary) and by considering the
voltage dynamics, as follows:

$\left\\{\begin{array}[]{ll}\dot{\mathbf{\theta}}_{i}=\omega_{i}\\\
\dot{\omega}_{i}=-\alpha_{i}\dot{\mathbf{\theta}}_{i}+P_{i}^{*}-\sum\limits_{j=1}^{N}E_{q,i}^{{}^{\prime}}B_{ij}E^{{}^{\prime}}_{q,j}\sin(\theta_{i}-\theta_{j})+u_{i}\\\
T^{{}^{\prime}}_{{d,i}_{0}}\dot{\overline{E}}_{q,i}^{{}^{\prime}}=E_{f,i}-E_{q,i}^{{}^{\prime}}+(X_{d,i}-X_{d,i}^{{}^{\prime}})\sum\limits_{j=1}^{N}B_{ij}E^{{}^{\prime}}_{q,j}\cos(\theta_{i}-\theta_{j})\\\
\tau_{g_{i}}\dot{u}_{i}=-u_{i}-\gamma_{i}\theta_{i}-\beta_{i}\omega_{i},\end{array}\right.$
(1)

where $i\in{1,...,N}$ denotes the node index in the network, $\theta$ the
voltage phase angle, $\omega$ the frequency, $E$ the voltage amplitude and $u$
the control. $\alpha$ is a damping constant, $P^{*}$ the power
consumed/generated at a node, $B$ the susceptibility matrix, while $T$ and
$\tau_{g}$ are time constants, $E_{f}$ is the rotor’s field voltage, and
$X_{d}$ and $X_{d}^{{}^{\prime}}$ are voltage dynamic parameters [21]. The
nodes (islanded power grids or synchronous machines), which compose the
network, are assumed to be all to all connected. For simplicity we assume that
the frequency controller acts instantaneously ($\tau_{g}$=0), such that it is
described as the proportional derivative control given by:

$u_{i}=-\gamma_{i}\theta_{i}-\beta_{i}\omega_{i}$. The derivative term
$\beta_{i}\dot{\mathbf{\theta}}_{i}$ can be absorbed into the damping term
$\alpha_{i}\dot{\mathbf{\theta}}_{i}$ of the swing equation. One can then also
rewrite the set of equation 1 by Eq. 2 substituting
$(X_{d,i}-X_{d,i}^{{}^{\prime}})$ by $X_{i}$, $T^{{}^{\prime}}_{{d,i}_{0}}$
and $E_{q,i}^{{}^{\prime}}$ by respectively $T_{{d,i}}$ and $E_{i}$ without
loss of generality

$\left\\{\begin{array}[]{ll}\dot{\mathbf{\theta}}_{i}=\omega_{i}\\\
\dot{\omega}_{i}=-\alpha_{i}\dot{\mathbf{\theta}}_{i}-\gamma_{i}\theta_{i}+P_{i}^{*}-\sum\limits_{j=1}^{N}E_{i}B_{ij}E_{j}\sin(\theta_{i}-\theta_{j})\\\
T_{{d,i}}\dot{\overline{E}}_{i}=E_{f,i}-E_{i}+X_{i}\sum\limits_{j=1}^{N}B_{ij}E_{j}\cos(\theta_{i}-\theta_{j}).\\\
\end{array}\right.$ (2)

###  II.2 Linear stability analysis

The stable operation of the electric network requires that the frequency,
voltage and difference of phases between connected nodes is constant, meaning
that the network is operating in a synchronous regime. Such system is said to
be linearly stable if, subjected to a small perturbation, it regains its
stable operation, in the case of a power grid its synchronous state.
Konstantin Sharafutdinov et al. have extensively studied the stability of such
electric network model considering the voltage dynamics as described in Eq.
(2), but without secondary control. Thus, based on their work and mainly on
the necessary and sufficient conditions of the uncontrolled system to be
linearly stable, the effects of the secondary control on the stability of the
network will be investigated.

Thus, to analyze the stability of the system with respect to small
perturbation, we linearize Eq. 2 around a steady state
$(\theta_{i}^{*},\omega_{i}^{*},E_{i}^{*})$. We denote small perturbation
around the steady state as $\theta_{i}=\theta_{i}^{*}+\delta\theta_{i}$,
$\omega_{i}=\omega_{i}^{*}+\delta\omega_{i}$, $E_{i}=E_{i}^{*}+\delta E_{i}$.
The linearization of the equation 2 leads to the Eq. 3, where
$\mathbf{X_{1}}$, $\mathbf{X_{2}}$ and $\mathbf{X_{3}}$ are n-dimensional
vectors of $\delta\theta_{i}$, $\delta\omega_{i}$ and $\delta E_{i}$.

$\left\\{\begin{array}[]{ll}\dot{\mathbf{X}}_{\mathbf{{}_{1}}}=\mathbf{X_{2}},\\\
\dot{\mathbf{X}}_{\mathbf{{}_{2}}}=-(\mathbf{P}+\mathbf{\Gamma})\mathbf{X_{1}}-\mathbf{A}\mathbf{X_{2}}-\mathbf{\Lambda}\mathbf{X_{3}}\\\
\dot{\mathbf{X}}_{\mathbf{{}_{3}}}=\mathbf{T}^{-1}\mathbf{\chi}\mathbf{\Lambda}\mathbf{X_{1}}+\mathbf{T}^{-1}(\mathbf{\chi}\mathbf{C}-\mathds{1})\mathbf{X_{3}},\end{array}\right.$
(3)

where $\mathbf{T}^{-1}$, $\mathbf{\Gamma}$, $\mathbf{A}$ and $\mathbf{\chi}$
are diagonal matrices with elements $T_{ii}$= $\frac{1}{T_{i}}$,
$\mathbf{\Gamma_{ii}}$= $\gamma_{i}$,$A_{ii}$= $\alpha_{i}$ and
$\chi_{ii}$=$X_{i}$ respectively, representing the relaxation time of the
transient voltage dynamics matrix,the control, the damping and the transient
reactance matrix of the synchronous machine. Matrices $\mathbf{P}$,
$\mathbf{\Lambda}$, $\mathbf{C}$, $\mathbf{\Gamma}$ $\in$ $\mathds{R}^{n\times
n}$ whose the elements are respectively defined as follows:

$P_{ij}=\left\\{\begin{array}[]{ll}-{E_{i}}^{*}B_{ij}{E_{j}}^{*}\cos({\theta_{i}}^{*}-{\theta_{j}}^{*}),&i\neq
j,\\\ \sum\limits_{l\neq
i}^{N}{E_{i}}^{*}B_{il}{E_{l}}^{*}\cos({\theta_{l}}^{*}-{\theta_{i}}^{*}),&i=j,\end{array}\right.$
(4)
${\Lambda}_{ij}=\left\\{\begin{array}[]{ll}{E_{j}}^{*}B_{ij}\sin({\theta_{i}}^{*}-{\theta_{j}}^{*}),&i\neq
j,\\\
-\sum\limits_{l=i}^{N}{E_{l}}^{*}B_{il}\sin({\theta_{l}}^{*}-{\theta_{i}}^{*}),&i=j,\end{array}\right.$
(5) ${C}_{ij}=B_{ij}\cos({\theta_{i}}^{*}-{\theta_{j}}^{*}).$ (6)

This later set of equations can be rewritten into the following compact form

$\frac{\text{d}}{\text{d}t}\left({\begin{array}[]{ccc}\mathbf{X_{1}}\\\
\mathbf{X_{2}}\\\
\mathbf{X_{3}}\end{array}}\right)=\underbrace{\left[{\begin{array}[]{ccc}\mathbf{0}&\mathds{1}&\mathbf{0}\\\
-(\mathbf{P}+\mathbf{\Gamma})&-\mathbf{A}&-\mathbf{\Lambda}\\\
\mathbf{T}^{-1}\mathbf{\chi}\mathbf{\Lambda}&\mathbf{0}&\mathbf{T}^{-1}(\mathbf{\chi}\mathbf{C}-\mathds{1})\end{array}}\right]}_{:=\mathbf{J}}\left[{\begin{array}[]{cc}\mathbf{X_{1}}\\\
\mathbf{X_{2}}\\\ \mathbf{X_{3}}\end{array}}\right].$ (7)

The evaluation of the Jacobian matrix $\mathbf{J}$ around an existing fixed
point $(\theta_{i}^{*},\omega_{i}^{*},E_{i}^{*})$ determines the stability of
the system around this fixed point. That is done by computing the eigenvalues
$\mu_{i}$ associated to this fixed point. The system is then linearly stable
if the maximal real part of the eigenvalues is negative. In this case,
perturbed, the system will regain its fixed point. But, if at least one
eigenvalue has a positive real part, then, the system is said to be linearly
unstable. This previous analysis is possible once we can determine the steady
state of the system and therefore compute the eigenvalues.
In most of the cases, the fixed points can be difficult to determine, hence
complicating the determination of eigenvalues. In such cases, the
determination of condition of existence of the fixed point is necessary and
the stability conditions of the system can be derived through some
mathematical formulations.
Sharafutdinov et al in [24] have provided for the uncontrolled system two
sufficient and necessary stability conditions grouped into two propositions
(proposition I and II in [24]), which both have to hold to have a linearly
stable fixed point. In fact, it is worth noting that the Jacobian matrix
$\mathbf{J}$ has one eigenvector $(1,0,0)$ with eigenvalue $\mu_{1}=0$, which
corresponds to a global phase shift of the synchronous machine [25, 24]. This
particular case, which does not have any physical implications for the
stability, is excluded from the stability analysis. Thereby we reduce the
space of possible states to $S_{\bot}$, which is the space perpendicular to
the solution $\theta_{i}^{*}$+c(1,0,0)/c $\in\mathds{R}$ and defined by

$S_{\bot}=\left\\{(\mathbf{X_{1}},\mathbf{X_{2}},\mathbf{X_{3}})\in\mathds{R}^{3n}/\mathbf{1}\mathbf{X_{1}}=0\\\
\right\\}.$ (8)

The Proposition I in [24] states that a given steady state
$(\theta_{i}^{*},\omega_{i}^{*},E_{i}^{*})$ is linearly stable if and only if

  * •

given the space
$S_{\bot}^{1}=\left\\{(\mathbf{X_{1}},\mathbf{X_{2}})\in\mathds{R}^{2N}/\mathbf{1}\mathbf{X_{1}}=0\\\
\right\\}$, the matrix $\mathbf{P+\Gamma}$ is positive definite on
$S_{\bot}^{1}$, furthermore,

  * •

the matrix
$\mathbf{C}-\mathbf{\chi}^{-1}+\mathbf{\Lambda}\mathbf{P}^{+}\mathbf{\Lambda}^{T}$
is negative definite.

Where $T$ denotes the transpose of a matrix, which is identical to its inverse
if it is an orthogonal matrix. In addition, the matrix $\mathbf{P}^{+}$,
represents the Moore-Penrose pseudo-inverse of $\mathbf{P}$. From this
proposition I, the only difference appearing in comparison to the study made
in Sharafutdinov et al in [24] is the matrix $\mathbf{\Gamma}$ in the sum
$\mathbf{P+\Gamma}$. Thus, from the proposition I, it is clear that the
secondary control represented by the matrix $\mathbf{\Gamma}$ can improve the
stability of the network with the necessary condition to be positive definite.
In our case, where the nodes are assumed to have the same system parameters
such are $\alpha_{i}=\alpha$, $\gamma_{i}=\gamma$, $T_{d,i}=T_{d}$ and
$E_{f,i}=E_{f}$. The matrix $\mathbf{\Gamma}$ is a diagonal matrix, which
elements $\gamma$ are always positive. Thus, all the eigenvalues of the
symmetric matrix $\mathbf{\Gamma}$ are positive, hence, the matrix
$\mathbf{\Gamma}$ is positive definite (See Appendix A).

###  II.3 Bulk dynamics

In order to analyze the stability of the network, we simplify the analysis by
focusing on the ensemble dynamics, i.e. we consider the _bulk_ or average
dynamics of the network. As in the previous stability analysis, we wish to
evaluate the impact of the secondary control through its parameter $\gamma$ on
the dynamics of the system. Thus, we first consider a system with a constant
voltage, which is a well known and studied case. Next, we include the voltage
dynamics and analyze the effects of the secondary control on the voltage.

####  II.3.1 Constant voltage

Let’ s consider Eq. (2) and assume the voltage to be constant. The resulting
equation is the simple Kuramoto model with secondary control, discussed for
example in [5]:

$\left\\{\begin{array}[]{ll}\dot{\mathbf{\theta}}_{i}=\omega_{i},\\\
\dot{\omega}_{i}=-\alpha_{i}\dot{\mathbf{\theta}}_{i}-\gamma_{i}\theta_{i}+P_{i}^{*}-\sum\limits_{j=1}^{N}E_{i}B_{ij}E_{j}\sin(\theta_{i}-\theta_{j}).\\\
\end{array}\right.$ (9)

Now, we take the average of the variables of Eq. (9) over the number of nodes
$n$ that form the network. We further assume that the node parameters
$\alpha_{i}=\alpha$, $\gamma_{i}=\gamma$ are identical for all the nodes.
Thus, we obtain the following differential equations:

$\left\\{\begin{array}[]{ll}\dot{\overline{\theta}}=\overline{\omega},\\\
\dot{\overline{\omega}}=-\alpha\dot{\overline{\theta}}-\gamma\overline{\theta}+\frac{1}{N}\sum_{i=1}^{N}{P_{i}^{*}}.\\\
\end{array}\right.$ (10)

Eq. (10) can easily be solved and the average value of the principal variable
are obtained as follows:

$\overline{\theta}(t)=C_{1}\exp{r_{1}t}+C_{2}\exp{r_{2}t}+\frac{\sum_{i=1}^{N}{P_{i}^{*}}}{N\gamma},$

where $C_{1}$ and $C_{2}$ are constants which are determined by the initial
conditions and $r_{1}$ and $r_{2}$ are expressed as follows:

$\left\\{\begin{array}[]{ll}r_{1}=-\frac{\alpha}{2}+\frac{\sqrt{\alpha^{2}-4\gamma}}{2},\\\
r_{2}=-\frac{\alpha}{2}-\frac{\sqrt{\alpha^{2}-4\gamma}}{2}.\\\
\end{array}\right.$ (11)

Thus, it is clear that for time tending to infinity and for any value of the
secondary control parameter $\gamma$ ($\gamma>0$), the average value of the
angle of rotation of each node in the network is constant and the average
value of the frequency is zero:

$\overline{\theta}(t)=\frac{\sum_{i=1}^{N}{P_{i}^{*}}}{N\gamma},\quad\bar{\omega}=0.$

As we can see, the average of the angle $\overline{\theta}$ for the controlled
system depends on the control parameter $\gamma$, the balance of the network
as well as the size of the network. For a balanced network, i.e. with
$\sum_{i=1}^{N}{P_{i}^{*}}=0$, the mean value of the angles is zero. For an
imbalanced and fixed network, the mean angle is inversely proportional to the
secondary control parameter, thus for a large control, the average of the
angles tends to zero.
For an uncontrolled electric network ($\gamma=0$), the average of the angles
is time varying and given by the following expression:

$\overline{\theta}(t)=D_{1}+D_{2}\exp{(-\alpha
t)}+\frac{\sum_{i=1}^{N}{P_{i}^{*}}}{N\alpha}t,$

where $D_{1}=-D_{2}=-{\sum_{i=1}^{N}{{P_{i}^{*}}}/{N\alpha^{2}}}$ are
constants determined by the initial conditions. Thus, the mean value of the
frequency of the network as a function of time $t$ is given by the Eq. (12):

$\overline{\omega}(t)=-\alpha{D_{2}}\exp{(-\alpha
t)}+\frac{\sum_{i=1}^{N}{P_{i}^{*}}}{N\alpha},$ (12)

this basically corresponds to the frequency deviation for an uncontrolled and
imbalanced network.

This imbalance can be easily absorbed in large networks (large $N$) and those
with large primary control (given by the $\alpha$ parameter).

We have so far shown that, when the voltage is constant, any imbalance in the
network induces the deviation of the mean frequency, which is constant and
non-zero for a network without control. This implies that the average value of
the angle is linearly increasing with the time. For the same network, but in
presence of secondary control in all the nodes, one notices that the average
frequency in the network is always zero after a long period of time and the
average angular in the network is a constant value proportional to the
disturbance and inversely proportional to the size of the network and the
secondary control parameter. We have therefore seen the effects of the
secondary control in such an electric network in which the voltage dynamics is
not considered. But how does the voltage dynamics change these results?

####  II.3.2 Dynamical voltage

In this part, we investigate the effects of the secondary control on the
dynamics of the voltage in an electric network. Thus, we consider the full
equation (2), summing up all three equations and dividing them by the total
number of nodes. Again, we assume that the parameters of the nodes are
constant and identical for all the nodes, obtaining the following equations:

$\left\\{\begin{array}[]{lll}\dot{\overline{\theta}}=\overline{\omega}\\\
\dot{\overline{\omega}}=-\alpha\dot{\overline{\theta}}-\gamma\overline{\theta}+\frac{1}{N}\sum_{i=1}^{N}{P_{i}^{*}}\\\
T_{d}\dot{\overline{E}}=E_{f}-\overline{E}+\frac{X}{N}\sum_{i=1}^{N}\sum\limits_{j=1}^{N}B_{ij}E_{j}\cos(\theta_{i}-\theta_{j})\end{array}\right.$
(13)

As shown in Eq. (13), the equations describing the dynamics of the mean
voltage (the third equation) contain the difference of phases of connected
nodes as argument of a cosine function. The presence of this term makes the
determination of the mean voltage fixed point very difficult. Nevertheless,
one can find the range of variation of $\overline{E}$ considering the
following relations:

$\begin{cases}B_{ij}=B_{ji}=\left\\{\begin{array}[]{ll}B_{0}&\mbox{if }i=j\\\
B_{1}&\mbox{else}\end{array}\right.\\\ -1\leq\cos{(\theta_{i}-\theta_{j})}\leq
1.\par\end{cases}$

Thus, the term with cosine in Eq.(13) is approximated by the following
relation:

$\begin{cases}-X(B_{0}+(N-1)B_{1})\overline{E}\leq\frac{X}{N}\sum_{i=1}^{N}\sum\limits_{j=1}^{N}B_{ij}E_{j}\cos(\theta_{i}-\theta_{j}),\\\
\frac{X}{N}\sum_{i=1}^{N}\sum\limits_{j=1}^{N}B_{ij}E_{j}\cos(\theta_{i}-\theta_{j})\leq
X(B_{0}+(n-1)B_{1})\overline{E}.\end{cases}$ (14)

The range of variation of the mean value can then be determined, solving the
following inequalities:

$\begin{cases}T_{d}\dot{\overline{E}}+\overline{E}-E_{f}\leq X[\
(B_{0}+(N-1)B_{1})]\ \overline{E},\\\
T_{d}\dot{\overline{E}}+\overline{E}-E_{f}\geq X[\ (B_{0}-(N-1)B_{1})]\
\overline{E}.\end{cases}$ (15)

Thus, the average of the voltage is bounded as follows:

$\begin{cases}\overline{E(t)}\leq
E_{0}\exp{(\frac{-(1-\sigma_{1})t}{T_{d}})}+\frac{E_{f}}{T_{d}(1-\sigma_{1})}(\
1-\exp{(\frac{-(1-\sigma_{1})t}{T_{d}})})\,\\\ \overline{E(t)}\geq
E_{0}\exp{(\frac{-(1+\sigma_{2})t}{T_{d}})}+\frac{E_{f}}{T_{d}(1+\sigma_{2})}(\
1-\exp{(\frac{-(1+\sigma_{2})t}{T_{d}})})\ ,\end{cases}$ (16)

where $E_{0}$ is the initial value of the mean voltage:
$\overline{E}(t=0)=E_{0}$, and $\sigma_{2}$ and $\sigma_{1}$ are given by:

$\begin{cases}\sigma_{1}=X[\ (B_{0}+(N-1)B_{1})],\ \\\ \sigma_{2}=-X[\
(B_{0}+(N-1)B_{1})].\ \end{cases}$

Evaluating all previous inequalities, we observe that the dynamics of the mean
voltage is strictly related to the size of the network. Thus, according to Eq.
(16), the mean voltage will be bounded between real values if and only if:

$\begin{cases}1-\sigma_{1}=1-X[\ (B_{0}+(N-1)B_{1})]\geq 0,\\\
1-\sigma_{2}=1+X[\ (B_{0}-(N-1)B_{1})]\geq 0.\\\ \end{cases}$ (17)

Thus the stability of the electric voltage is achieved if the network size and
the network coupling are chosen in the range given by the following
inequality:

$1-X\frac{1-B_{0}}{B_{1}}\leq N\leq 1+X\frac{1-B_{0}}{B_{1}}$ (18)

In the present study, the values of the system’s parameters are set as: $X=1$,
$B_{0}=-0.8$ and $B_{1}=1$, following literature values [21, 22]. Thus, the
mean voltage of the network of size $N$ will be bounded if and only if:

${\color[rgb]{0,0,0}N=2}$ (19)

Indeed, let’s consider an electric network with the same parameters as the
studied case and composed of $n$ nodes, where each node represents a small
region or a city. The result above states that an extension of the current
network composed of $N=2$ nodes by connecting it with additional nodes leads
to the instability of the new extended network in terms of voltage.

Thus, from this bulk analysis, it appears that the secondary frequency control
has explicitly no effects on the stabilization of the mean voltage, but does
have an effect when it comes to stabilizing the mean frequency of the network.
This implies therefore the necessity of having another form of control for the
mean voltage in larger electric networks.

##  III Numerical analysis

In the previous section, we analysed the linear stability as well as the bulk
dynamics analytically. In this part, we present a numerical analysis for N=2
and $N\geq 2$ nodes with the objective to evaluate the effects of the
frequency secondary control ($\gamma$) on the network, thereby complementing
the analytical results.

###  III.1 2-node system

The two node system considered here consists of a generator node of power
$P_{1}>0$ and a consumer node of power $P_{2}<0$. The two nodes represent
simply two interconnected power grids (or synchronous machines) with equal
control parameters $\gamma_{1}=\gamma_{2}=\gamma$ and are described by the Eq.
(2) and parameters, such as $T_{d}$, inspired by literature values [21, 22].
Initially, without disturbance, the two systems evolve till reaching a stable
final state, in which the variables of each node tend to a steady state as
shown in Fig.1 without control (first column in black) and with control
(second column in red).

Figure 1: An unperturbed system converges to a steady state with and without
control. [Description of figure] parameters of simulation are the following:
$\alpha_{1}=\alpha_{2}=0.2$, $T_{d,1}=T_{d,2}=2$,$E_{f,1}=E_{f,2}=1.$ ,
$E_{1}(0)=E_{2}(0)=1.14$, $\omega_{1}(0)=\omega_{2}(0)=50$ ,
$\theta_{1}(0)=\theta_{2}(0)=0$, $P_{dist}$=0, $\gamma_{1}=\gamma_{2}=0$
(first column), $\gamma_{1}=\gamma_{2}=1$ (second column)

One observes that in both cases the network is synchronized, meaning that the
frequency at each node is equal to $\omega_{syn}=50$ $Hz$ and the difference
of phases between the connected node tends to a constant value. The final
states reached by the voltages in the controlled case are sightly greater than
the the ones of the uncontrolled case. This slight increase results from the
effects of the secondary control parameter at each node, which affects the
dynamics of the voltage through the angles. In fact, at a constant voltage,
the secondary control is basically the governor which activates the online or
offline substations, thereby increasing or reducing the power generated in
order to balance the power in the system ($\omega=0$) [7, 5, 26]. Thus, once
the frequency is brought back to its nominal value, the phases in each node
tend to constant values, which basically implies a constant difference of
phases between connected nodes, hence constant voltage. Without control, the
two nodes of the system evolve until sharing a constant amount of power and
reach a power balanced state.

Figure 2: Secondary control stabilizes a 2 node system. [Description of
figure] parameters of simulation are the following:
$\alpha_{1}=\alpha_{2}=0.2$, $T_{d,1}=T_{d,2}=2$,$E_{f,1}=E_{f,2}=1.$,
$E_{1}(0)=E_{2}(0)=1.14$, $\omega_{1}(0)=\omega_{2}(0)=50$,
$\theta_{1}(0)=\theta_{2}(0)=0$, $P_{dist}$=1, $\gamma_{1}=\gamma_{2}=0$
(first column), $\gamma_{1}=\gamma_{2}=1$ (second column)

Next, we consider a perturbation to the power system: We assume that from time
$t=40$ $s$ to $t=42$ $s$ , the power $P_{1}$ at node 1 experiences a gradual
increase of power from 0 to $P_{dist}=1$, resulting in a sudden increase of
its corresponding frequency, as shown in Fig. 2 (dashed lines) in the
controlled (red) and uncontrolled (black) case. This perturbation leads to an
instability of the uncontrolled case, visible by the emergence of periodic
oscillation in the frequency and the voltages around states far away from
their original states,while the corresponding phases are rising in opposite
directions. On the other hand, the variables in the controlled system (plot in
red in Fig.2) converge to new steady states for the voltages and phases, while
the corresponding frequencies return to their synchronous states as before the
perturbation. Thus, the secondary control not only stabilized the frequencies
but also stabilized the voltage, thereby avoiding the system to drop into an
unstable regime. Thus, the secondary frequency control acts as a damping for
the voltages. This is a very similar behavior as the primary frequency control
has on the frequency through the damping coefficient $\alpha$ on the frequency
of the uncontrolled ($\gamma=0$) network. Hence, we can consider secondary
frequency control acts somehow as a primary control for the voltage.

Figure 3: Uncontrolled networks display frequency instability and large
voltage oscillations. Illustration of average of the phase (a), frequency
deviation (b) and the voltage (c) as a function of the time for different size
of the network, with the corresponding zoom of the voltage dynamics given at
the bottom of the figure. Parameter values: $\gamma=0$, $\alpha=0.2$,
$Pdist=1$, $X=1$, $E_{f}=1$, $E_{i}(t=0)=1.14$, $T_{d}=1$, $B_{0}=-0.8$ and
$B_{1}=1$.

Figure 4: Secondary control stabilizes the frequency and voltage following
initial oscillations. Illustration of average of the phase (a), frequency
deviation (b) and the voltage (c) as a function of the time for different size
of the network, with the corresponding zoom of the voltage dynamics given at
the bottom of the figure. Parameter values: $\gamma=1$, $\alpha=0.2$,
$Pdist=1$

###  III.2 Bulk dynamics

Moving towards the $N>2$, case, we now present the numerical results
complementing the analytical bulk analysis. Thus, the question to be answered
is the following: how does the secondary control affect the mean voltage
dynamics in networks? Thus, we consider the electric network, whose dynamics
is governed by the Eq. (2) and we analyse the dynamics of the mean voltage and
frequency for different network sizes. We consider for these purposes that we
have an interconnected network constituted of $n$ nodes, which can represent
$n$ interconnected isolated power grids or synchronous machines. So far, we
assumed that the nodes are all-to-all-coupled, each node having equal power
$P_{i}$ in absolute value, and the network without perturbation is power
balanced ($\sum_{i}{P_{i}=0}$). Such all-to-all-coupled networks naturally
emerge after Kron reduction of any network topology [27]. In all the studied
scenarios, the system is perturbed from the time $t=40$ $s$ to $t=42$ $s$ by
gradually increasing the power at a single node from $0$ to $P_{dist}$.

First, we consider the uncontrolled ($\gamma=0$) case and plot the deviation
of the mean phase, frequency and voltage as function of time for different
network size, and without secondary control in Fig. 3. We observe a clear
agreement with the analytical predictions given by the expression of the
frequency deviation in Eq. (12), meaning that the mean frequency for a long
time tends to a constant. Since we assumed that the primary control
parameter/the inertia of each power plant is constant, we observe that the
mean of the frequency deviation decreases when the size of the network
increases. The corresponding curve of the average voltage is plotted at the
third column of the first row. Complementing the analytical results, we
observe large transient voltage dynamics, which increase in amplitude with
network size. To better observe the dynamics of the mean value of the voltage
towards the end of our simulation window, we provide zooms in the lower row of
Fig. 3 for each size $n$ of the network. For the two nodes system, the average
frequency tends to a constant value. The mean voltage conversely is
oscillating around a stable value $\overline{E}\approx 0.55$ $Volt$. For
larger values of $n$ on the other hand, the mean voltage is fluctuating and
displays large periodic peaks. These regular peaks persist throughout the
running time of our simulation and are observed for $N$ = $10$ and $N$ = $20$.
For $N$ = $50$, after the transient, the mean voltage mostly fluctuates around
a constant value and only displays one peak, likely pointing to a longer
periodicity of these peaks due to the larger network. The small fluctuations
observed are probably due to the difference of phases existing between the
connected nodes. In this uncontrolled studied case, the system is clearly not
a synchronous stable electric network.

Now, we include secondary control ($\gamma>0$) and repeat the same simulations
as before. Fig. 4 (top) shows the time evolution of the mean frequency and the
mean voltage of the considered network in presence of the secondary control.
First, we notice that the mean frequency deviation tend to zero for every
network size, as we predicted in our linear stability analysis. Hence, we
focus on the mean voltage dynamics, for which we could not derive an equality
from the Eq.(13) but only constrained its range of fluctuations. Thus, only
these numerical simulation can tell us how the voltage evolves with time.
Similar to the frequency, after an initial transient phase, the mean voltage
tends to a constant value for all the considered sizes of the network.

Again, to highlight the voltage dynamics towards the end of our simulation
window, we provide a zoom of the voltages in the second row of figure 4. For
the two node system, we observe that the mean voltage tends to steady state,
and it is not oscillating as in the corresponding uncontrolled case. Thus, the
secondary control clearly stabilizes the voltage. In fact, without secondary
control, the phase of each node is increasing continuously with the time (see
Fig. 2), leading to the oscillation of the voltage. The secondary frequency
control on other hand stabilizes the phases, hence the cosine function in the
expression of each voltage becomes a constant, leading to constant voltage at
each node. For larger networks, it may happen that due to the power flow
between connected nodes, the phases slightly vary, leading at some times to
the fluctuations of the voltage. This justifies some rare peaks observed in
the mean voltage. Nevertheless, the magnitude of the voltage and its peaks in
the controlled network mostly remains lower than the one in the uncontrolled
network. In addition, we observe that the steady state reached in the end
decreases when the network size increases.

###  III.3 Relaxation time

We have shown in the previous analysis how the dynamics of the phase, the
frequency, the voltage as well as their corresponding mean values evolve with
the secondary parameter $\gamma$, once the network is perturbed. In fact, once
perturbed from its stable state, the variables of the system vary until they
reach new steady states, often different from their original steady states
without perturbation. Now, we investigate how long it takes the system to
relax to its (new) steady state, by computing the return time or relaxation
time. According to the bulk dynamics of the network, the return time of the
frequency and phase average depends strongly on the damping of the system, and
slightly on the secondary control parameter, which impacts mostly the
oscillatory regime of these variables. If the return time of the mean
frequency and phase can be obtained explicitly, this cannot so easily be
achieved for the mean value of the voltage. Hence, we again use numerical
computations to obtain the return time for different control parameters.
In particular, we consider the previously described 2 nodes system, including
voltage dynamics. For this system, we aim to determine the evolution of the
return time as a function of the secondary control parameter $\gamma$. The
network is perturbed as previously, by gradually decreasing the power at the
node $1$ from $P_{1}$ to $P_{1}+P_{dist}$ ($P_{dist}=-1$) from time t=$40s$ to
t=$42s$ and all other parameters also remain unchanged. The relaxation time or
return time is numerically computed by recording the time taken by the voltage
at a node to regain a steady state after perturbation. Thus, we define by
$\overline{E}(t)$ the mean voltage at the time $t$ and $\overline{E}(t-T)$ the
mean at the time $t=t-T$, where $T$ is a characteristic time. We define also
by $\xi$ the numerical tolerance. The mean voltage is considered stable after
perturbation of the network if and only if
$|\overline{E}(t)-\overline{E}(t-T)|\leq$ $\xi$. Note that we focus here on
the relaxation and therefore of the voltage, not the frequency, as we are
interested in quantifying the impact of frequency control on voltage
stability. The so computed return times decrease with increasing secondary
control parameter $\gamma$, see Fig. 5. First, we note that for any non-zero
value of the secondary control parameter $\gamma$, the system after
perturbation always regains a steady state, i.e. the secondary frequency
control guarantees stability and return. Secondly, we observe that the return
time is reduced by increasing the control amplitude, approaching zero for
sufficiently large control $\gamma$. This means that the secondary frequency
control is strong enough to immediately compensate the introduced disturbance
and the system never leaves its original fixed point.

Figure 5: Secondary control stabilizes the system by reducing the return time.
We plot the return time of the mean voltage when the control parameter
increases in the two node system. The parameters of simulation are the
following: $\alpha_{1}=\alpha_{2}=0.2$, $T_{d,1}=T_{d,2}=2$,
$E_{f,1}=E_{f,2}=1$, $P_{dist}=1$ $N=2$ nodes.  Figure 6: Exploring
heterogeneous parameters and non-all-to-all-coupling. Illustration of the
topology (a), the average of the frequency deviation (b) and the voltage (c)
as a function of the time for the controlled (black) and controlled (red)
network. Parameter values:
$P=[-0.4,0.53,-0.51,0.56,0.52,0.48,-0.55,-0.45,0.491,0.509,-0.482,-0.518,-0.46,-0.64,0.42,0.58,-0.5,0.5,0.35,-0.45]$,
primary control $\alpha_{i}=0.2|P_{i}|$, secondary control
$\gamma_{i}=|P_{i}|$, while $P_{dist}$, $X$ and $B$ are chosen as in the other
simulations.

###  III.4 Case study: Heterogeneous parameters and network influence

So far we have assumed networks that are all-to-all-coupled and have
homogeneous parameters, i.e. identical loads, generation and control values.
We used this simplification to derive analytical results. To demonstrate that
our results are in principle also applicable to more general systems we
consider one non-all-to-all-coupled system with heterogeneous parameters, see
Fig. 6a. In particular, we simulate the dynamical behavior of 20 nodes
connected to a common bus with power values $P$ randomly drawn between -0.7
and 0.7. We assume nodes with large absolute value of $P$ to provide more
control power and hence set the control values proportional to the absolute
power value: $\alpha_{i}=0.2|P_{i}|$, secondary control $\gamma_{i}=|P_{i}|$,
where $|...|$ denotes the absolute value. Analogue to the earlier analysis, we
observe that the average frequency can only be restored once secondary control
is used, see Fig. 6b. Meanwhile, the voltage dynamics is not controlled and
oscillates (Fig. 6c) since we have more than two nodes, consistent with our
analytical results. How exactly non-homogeneous parameters and network
topology affect the dynamics and controlability of both frequency and voltage
is beyond the scope of this study.

##  IV Conclusion

This article has investigated the effects of the secondary frequency control
on the voltage dynamics of electric networks of different sizes. We considered
an simple networks, both all-to-all-coupled as well as a bus topology. Our
analytical linear stability analysis of the network has shown that the
secondary control can guarantee the stability of the network. In addition,
considering the network as a simplified bulk, we have demonstrated that the
stability of the mean phase and frequency are independent of the mean voltage
of the network. On the other hand, the mean voltage does depend on the nodes’
phases. The different numerical simulations computed for the perturbed network
and in presence of the secondary control have shown that, the secondary
control actually plays the role of a primary control for the voltage. The
frequency secondary control after perturbation stops the variation of the
voltage and stabilizes it to a new steady value different to the original one.

Our results showcase how voltage stability and secondary frequency control
should be considered in other power system stability analyses: Including only
primary frequency control might stabilize the frequency but does not guarantee
any voltage stability. Hence, the voltage should then explicitly be considered
when assessing stability in power systems. Meanwhile, if secondary frequency
control is applied, the voltage stability is no longer an immediate concern
and might be neglected, especially if only the short-term stability is of
interest.

In the future, it would be interesting to complement the secondary frequency
control, which acts as an effective "primary voltage control" by a "secondary
voltage control" to bring the voltage back within its operational boundaries.
Furthermore, how precisely network topology and heterogeneous parameters
affect the voltage and frequency stability and return times still remains a
mostly open question for now. Similarly, the choice of specific parameters,
such as $T_{d}$ and $X_{i}$ should be investigated further. Finally, our
linear stability and return time analysis could be supplemented by a detailed
analysis of the basin of attraction [28].

### Acknowledgments

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 840825.
E.B.T.T. and P.P. acknowledge funding from the European Regional Development
Fund under Grant No. CZ.02.1.02/0.0/0.0/15003/0000493, "CENDYNMAT - Centre of
Excellence for Nonlinear Dynamics Behaviour of Advanced Materials in
Engineering" and the Academy of Sciences CR under Grant Strategy AV21, VP03:
"Efficient energy conversion and storage, Vibrodiagnostics of rotating blades
of rotary machines in power engineering.
E.B.T.T., D.G., and P.C. acknowledge funding from the Ministerio de Ciencia e
Innovación (Spain), the Agencia Estatal de Investigación (AEI, Spain), and the
Fondo Europeo de DesarrolloRegional (FEDER, EU) under Grant No. PACSS
(RTI2018-093732-B-C22) and the Maria de Maeztu program for Units of Excellence
in R&D (No. MDM-2017-0711). E.B.T.T. also acknowledges the fellowship from the
AEI and MINEICO, Spain under the FPI program(No. FIS2015-63628-CZ-Z-R).

### Competing interests

The authors declare no competing interests.

## References

  * Machowski et al. [2011] J. Machowski, J. Bialek, and J. Bumby, _Power System Dynamics: Stability and Control_ (John Wiley & Sons, 2011).
  * Mancarella et al. [2021] P. Mancarella, J. Moriarty, A. Philpott, A. Veraart, S. Zachary, and B. Zwart, _Introduction: the mathematics of energy systems_ (2021).
  * Tchawou and Woafo [2014] E. Tchawou and P. Woafo, Nonlinear Engineering 3, 89 (2014).
  * Tékam et al. [2014] G. O. Tékam, E. T. Tchuisseu, C. K. Kwuimy, and P. Woafo, Nonlinear Dynamics 76, 1561 (2014).
  * Tchuisseu et al. [2018] E. B. T. Tchuisseu, D. Gomila, P. Colet, D. Witthaut, M. Timme, and B. Schäfer, New Journal of Physics 20, 083005 (2018).
  * Witthaut and Timme [2012] D. Witthaut and M. Timme, New Journal of Physics 14, 083036 (2012).
  * Tchuisseu et al. [2017] E. T. Tchuisseu, D. Gomila, D. Brunner, and P. Colet, Physical Review E 96, 022302 (2017).
  * Witthaut and Timme [2013] D. Witthaut and M. Timme, The European Physical Journal B 86, 1 (2013).
  * Schäfer et al. [2015] B. Schäfer, M. Matthiae, M. Timme, and D. Witthaut, New Journal of Physics 17, 015002 (2015).
  * Dongmo and Woafo [2015] E. D. Dongmo and P. Woafo, The European Physical Journal B 88, 1 (2015).
  * Dongmo et al. [2017] E. D. Dongmo, P. Colet, and P. Woafo, The European Physical Journal B 90, 6 (2017).
  * Filatrella et al. [2008] G. Filatrella, A. H. Nielsen, and N. F. Pedersen, The European Physical Journal B 61, 485 (2008).
  * Gorjão et al. [2020] L. R. Gorjão, M. Anvari, H. Kantz, C. Beck, D. Witthaut, M. Timme, and B. Schäfer, IEEE Access 8, 43082 (2020).
  * Sun et al. [2019] H. Sun, Q. Guo, J. Qi, V. Ajjarapu, R. Bravo, J. Chow, Z. Li, R. Moghe, E. Nasr-Azadani, U. Tamrakar, et al., IEEE Transactions on Power Systems 34, 2790 (2019).
  * Weitenberg et al. [2017] E. Weitenberg, Y. Jiang, C. Zhao, E. Mallada, C. De Persis, and F. Dörfler, arXiv preprint arXiv:1711.07332 (2017).
  * Tyloo and Jacquod [2020] M. Tyloo and P. Jacquod, IEEE Control Systems Letters 5, 929 (2020).
  * Böttcher et al. [2020] P. C. Böttcher, A. Otto, S. Kettemann, and C. Agert, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 013122 (2020).
  * Nishikawa and Motter [2015] T. Nishikawa and A. E. Motter, New Journal of Physics 17, 015012 (2015).
  * Schäfer et al. [2016] B. Schäfer, C. Grabow, S. Auer, J. Kurths, D. Witthaut, and M. Timme, The European Physical Journal Special Topics 225, 569 (2016).
  * Schäfer et al. [2017] B. Schäfer, M. Matthiae, X. Zhang, M. Rohden, M. Timme, and D. Witthaut, Physical Review E 95, 060203 (2017).
  * Schmietendorf et al. [2014] K. Schmietendorf, J. Peinke, R. Friedrich, and O. Kamps, Euroean Physical Journal Special Topics 223, 2577 (2014).
  * Auer et al. [2016] S. Auer, K. Kleis, P. Schultz, J. Kurths, and F. Hellmann, Euroean Physical Journal Special Topics 225, 609 (2016).
  * Dörfler and Bullo [2014] F. Dörfler and F. Bullo, Automatica 50, 1539 (2014).
  * Sharafutdinov et al. [2018] K. Sharafutdinov, L. Rydin Gorjão, M. Matthiae, T. Faulwasser, and D. Witthaut, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 033117 (2018).
  * Manik et al. [2014] D. Manik, D. Witthaut, B. Schäfer, M. Matthiae, A. Sorge, M. Rohden, E. Katifori, and M. Timme, The European Physical Journal Special Topics 223, 2527 (2014).
  * Tchuisseu et al. [2019] E. T. Tchuisseu, D. Gomila, and P. Colet, International Journal of Electrical Power & Energy Systems 108, 145 (2019).
  * Dorfler and Bullo [2013] F. Dorfler and F. Bullo, Circuits and Systems I: Regular Papers, IEEE Transactions on 60, 150 (2013).
  * Hellmann et al. [2016] F. Hellmann, P. Schultz, C. Grabow, J. Heitzig, and J. Kurths, Scientific Reports 6 (2016).
  * Horn [1962] A. Horn, Pacific Journal of Mathematics 12, 225 (1962).

##  Appendix A Eigenvalues of sum of matrices

To compute the eigenvalues of the sum of two matrices, as used in this paper,
we follow this simple line of reasoning: The eigenvalue equation of matrix $A$
are obtained by:

$\det{\left(A-\lambda I\right)}=0,$ (20)

with identity matrix $I$ and where $\lambda_{1},\lambda_{2},...,\lambda_{N}$
are the eigenvalues of $A$. Now suppose we add a diagonal matrix to $A$, i.e.
$B=A+\gamma I$, then the eigenvalue equation reads

$\displaystyle\det{\left(B-\tilde{\lambda}I\right)}$ $\displaystyle=$
$\displaystyle 0,$ (21) $\displaystyle\det{\left(A+\gamma
I-\tilde{\lambda}I\right)}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\det{\left(A-\left(\tilde{\lambda}-\gamma\right)I\right)}$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\det{\left(A-\lambda
I\right)}$ $\displaystyle=$ $\displaystyle 0,$

i.e. the eigenvalues $\tilde{\lambda}$ of $B=A+\gamma I$ are given as
$\tilde{\lambda}_{k}=\lambda_{k}+\gamma$. So if the eigenvalues are ordered as
$\lambda_{1}\leq\lambda_{2}\leq...$ and $\lambda_{1}<0$ is the smallest
eigenvalue, we can choose any $\gamma>\lambda_{1}$ so that $B=A+\gamma I$ has
only positive eigenvalues, i.e. it the matrix $B$ is positive definite.

The interested reader might consult more general results on the eigenvalues of
two Hermitian matrices [29].


# Stable patterns in the Lugiato-Lefever equation with a confined vortex pump

Shatrughna Kumar1 Wesley B. Cardoso2 Boris A. Malomed1,3 1Department of
Physical Electronics, School of Electrical Engineering, Faculty of
Engineering, and Center for Light-Matter Interaction, Tel Aviv University, Tel
Aviv 69978, Israel  2Instituto de Física, Universidade Federal de Goiás,
74.690-900, Goiânia, Goiás, Brazil  3Instituto de Alta Investigación,
Universidad de Tarapacá, Casilla 7D, Arica, Chile

###### Abstract

We introduce a model of a passive optical cavity based on a novel variety of
the two-dimensional Lugiato-Lefever equation, with a localized pump carrying
intrinsic vorticity $S$, and the cubic or cubic-quintic nonlinearity. Up to
$S=5$, stable confined vortex-ring states (vortex pixels) are produced by
means of a variational approximation and in a numerical form. Surprisingly,
vast stability areas of the vortex states are found, for both the self-
focusing and defocusing signs of the nonlinearity, in the plane of the pump
and loss parameters. When the vortex-rings are unstable, they are destroyed by
azimuthal perturbations which break the axial symmetry. The results suggest
new possibilities for mode manipulations in passive nonlinear photonic media
by means of appropriately designed pump beams.

##  I Introduction and the model

Optical solitons are a broad class of self-trapped states maintained by the
interplay of nonlinearity and dispersion or diffraction in diverse photonic
media [1, 2]. In addition to that, dissipative optical solitons are supported
by the equilibrium of loss and gain or pump, concomitant to the nonlinearity-
dispersion/diffraction balance [3, 4]. Dissipative solitons have been studied
in detail, theoretically and experimentally, in active setups, with the loss
compensated by local gain (essentially, provided by lasing), being modelled by
one- and two-dimensional (1D and 2D) equations of the complex Ginzburg-Landau
(CGL) type [5, 6].

In passive nonlinear optical cavities, the losses are balanced by the pump
field supplied by external laser beams, with the appropriate models provided
by the Lugiato-Lefever (LL) equations [7]. This setting was also studied in
the 1D and 2D forms [8, 9, 10, 11]. Widely applied in nonlinear optics,
equations of the LL type play a crucial role in understanding fundamental
phenomena such as the modulation instability (MI) and pattern formation in
dissipative environments [8]-[24]. The relevance of these models extends to
the exploration of complex dynamics of various nonlinear photonic modes,
tremendously important applications being the generation of Kerr solitons and
frequency combs in passive cavities [12]-[24], as well as the generation of
terahertz radiation [25]. In addition to rectilinear cavity resonators,
circular ones can be used too [26]. In many cases, they operate in the
whispering-gallery regime [27]-[32].

In most cases, solutions of the one- and two-dimensional LL equations are
looked for under the action of the spatially uniform pump, which approximately
corresponds to the usual experimental setup. However, the use of localized
(focused) pump beams is possible too, which makes it relevant to consider LL
equations with the respective shape of the pump terms. In fact, truly
localized optical modes in the cavities can be created only in this case,
otherwise the uniform pump supports nonzero background of the optical field.
In particular, exact analytical solutions of the LL equations with the 1D pump
represented by the delta-function, and approximate solutions maintained by the
2D pump in the form of a Gaussian were reported in Ref. [33]. In Ref. [26],
the LL equation for the ring resonator with localized pump and loss terms
produced nonlinear resonances leading to multistability of nonlinear modes and
coexisting solitons, that are associated with spectrally distinct frequency
combs.

Further, solutions for fully localized robust pixels with zero background were
produced by the 2D LL equation, incorporating the spatially uniform pump,
self-focusing or defocusing cubic nonlinearity, and a tight confining
harmonic-oscillator potential [34]. Additionally, this model with a vorticity-
carrying pump gives rise to stable vortex pixels. In particular, in the case
of the self-defocusing sign of the nonlinear term, the pixels with zero
vorticity and ones with vorticity $S=1$ were predicted analytically by means
of the Thomas-Fermi (TF) approximation.

In this work, we introduce the 2D LL equation for complex amplitude field
$u\left(x,y,t\right)$ of the light field, with the cubic or cubic-quintic
nonlinearity:

$\frac{\partial u}{\partial t}=-\alpha
u+\frac{i}{2}\nabla^{2}u+i\sigma\left(\left|u\right|^{2}-\eta^{2}\right)u-ig\left(\left|u\right|^{4}-\eta^{4}\right)u+f(r)e^{iS\theta},$
(1)

and a confined pump beam represented by factor $f(r)$. Here $i$ is the
imaginary unit, $\alpha>0$ is the loss parameter, $\sigma=+1$ and $-1$
corresponds, respectively, to the self-focusing and defocusing Kerr (cubic)
nonlinearity, $g>0$ or $g<0$ represents the self-defocusing or focusing
quintic nonlinearity (which often occurs in optical media [35, 36, 37], in
addition to the cubic term), parameter $\eta$ defines the cavity mismatch,
which is $\sigma\eta^{2}-g\eta^{4}$ (the coefficient multiplying the linear
term $\sim iu$), in terms of the linearized LL equation, and

$f(r)=f_{0}r^{S}\exp\left(-r^{2}/W^{2}\right)$ (2)

written in terms of polar coordinates $\left(r,\theta\right)$, corresponds to
the confined pump beam with real amplitude $f_{0}$, radial width $W$, and
integer vorticity $S\geq 1$. Vortex beams, shaped by the passage of the usual
laser beam through an appropriate phase mask, are available in the experiment
[37].

Equation (1) is written in the scaled form. All figures are plotted below in
the same notation. In physical units, $r=1$ and $t=1$ normally correspond to
$\char 126\relax 50$ $\mathrm{\mu}$m and $\sim 50$ ps, respectively. Then, the
typical width $W=2$, considered below, corresponds to the pump beam with
diameter $\sim 100$ $\mathrm{\mu}$m, which an experimentally relevant value.
Accordingly, the characteristic evolution time in simulations presented below,
$\grave{t}\sim 100$, corresponds to times $\char 126\relax 5$ ns.

Stationary solutions of Eq. (1) are characterized by values of the total power
(alias norm),

$P=\int_{-\infty}^{+\infty}dx\int_{-\infty}^{+\infty}dy~{}\left|u\left(x,y\right)\right|^{2}\equiv
2\pi\int_{0}^{\infty}\left|u(r,t)\right|^{2}rdr,$ (3)

and angular momentum,

$M=i\int_{-\infty}^{+\infty}dx\int_{-\infty}^{+\infty}dyu^{\ast}\left(y\partial_{x}u-x\partial_{y}u\right)dxdy$
(4)

(with $\ast$ standing for the complex conjugate), even if the power and
angular momentum are not dynamical invariants of the dissipative equation (1).
In the case of the axisymmetric solutions with vorticity $S$, i.e.,
$u\left(x,y\right)=u(r)e^{iS\theta}$ [38, 39, 40], the expressions for the
power and angular momentum are simplified:

$P=2\pi\int_{0}^{\infty}\left|u(r)\right|^{2}rdr,~{}M=SP.$ (5)

Our objective is to construct _stable_ ring-shaped vortex solitons
(representing vortex pixels, in terms of plausible applications), as localized
solutions of Eq. (1) with the same $S$ as in the pump term (2). The stability
is a challenging problem, as it is well known from the work with models based
on the nonlinear Schrödinger and CGL equations that (in the absence of a tight
confining potential) vortex-ring solitons are normally vulnerable to the
splitting instability. In the case of a narrow ring shape, the splitting
instability may be considered as quasi-one-dimensional MI of the ring against
azimuthal perturbations which break its axial symmetry [40, 2]. The azimuthal
MI is driven by the self-focusing nonlinearity, and inhibited by the self-
defocusing.

To produce stationary solutions for the vortex solitons in an approximate
analytical form (parallel to the numerical solution), we employ a variational
approximation (VA). Our results identify regions of the existence and
stability of the vortex solitons with $1\leq S\leq 4$ in the space of
parameters of Eqs. (1) and (2) (in particular, in the plane of
$\left(f_{0},\alpha\right)$) for both signs of the cubic nonlinearity,
$\sigma=\pm 1$, while the mismatch parameter is fixed to be $\eta=1$ by dint
of scaling. The stability areas are vast, provided that the loss coefficient
$\alpha$ is, roughly speaking, not too small. A majority of the results are
produced for the pure cubic model, with $g=0$, but the effect of the quintic
term, with $g\neq 0$, is considered too. Quite surprisingly, a stability area
for the vortices with $S\leq 3$ is found even in the case of $\sigma=+1,g<0$,
when both the cubic and quintic terms are self-focusing, which usually implies
strong propensity to the azimuthal instability of the vortex rings [40].

The rest of the paper is structured as follows. The analytical approach, based
on the appropriate VA, is presented in Sec. II. An asymptotic expression for
the tail of the vortex solitons, decaying at $r\rightarrow\infty$, is found
too in that section. Systematically produced numerical results for the shape
and stability of the vortex solitons are collected (and compared to the VA
predictions) in Sec. III. The paper is concluded by Sec. IV.

##  II Analytical considerations

This section summarizes the analytical part of the work and results produced
by this part. Two directions of the analytical considerations for the present
model are possible: the investigation of the decaying “tails" of the localized
stationary states, in the framework of the linearized model, and detailed
development of VA for the full model, including the nonlinear terms.

###  II.1 Asymptotic forms of the vortex solitons

Direct consideration of the linearized version of Eqs. (1) and (2) readily
produces an explicit result for the soliton’s tail decaying at
$r\rightarrow\infty$:

$u(r,\theta)\approx(i/2)W^{4}r^{S-2}\exp\left(-r^{2}/W^{2}+iS\theta\right),$
(6)

with the power of the pre-Gaussian factor, $r^{S-2}$, which is lower than that
in the pump term, $r^{S}$. Due to this feature, the asymptotic expression (6)
formally predicts a maximum of local power $\left|u(r)\right|^{2}$ at $S>2$,
at

$r^{2}=r_{\max}^{2}\equiv\sqrt{(S/2-1)}W.$ (7)

A local maximum is indeed observed in numerically found radial profiles of all
vortex solitons (see Figs. 2, 4, 7(b), 9(b), 10, and 11 below, for $S=1$, $4$,
$2$, $2$, $1$, and $3$, respectively). In fact, for these cases Eq. (7)
predicts values of $r_{\max}$ which are smaller by a factor $\simeq 0.6$ than
the actually observed positions of the maxima. The discrepancy is explained by
the fact that the asymptotic expression (6) is valid at values of $r$ which
are essentially larger than $r_{\max}$.

In a looser form, one can try to construct an asymptotic approximation for the
solution at moderately large $r$, by adopting the ansatz which follows the
functional form of the pump term (2), viz.,

$u(r,\theta)\approx\mathrm{u}(r)r^{S}\exp\left(-r^{2}/W^{2}+iS\theta\right),$
(8)

where $\mathrm{u}(r)$ is a complex slowly varying function, in comparison with
those which are explicitly present in ansatz (8). Substituting the ansatz in
Eq. (1) and omitting derivatives of the slowly varying function, one can
develop an approach which is akin to the TF approximation applied to the model
with the tight trapping potential in Ref. [33]. The result for the linearized
version of Eq. (1), which implies a small amplitude of the mode pinned to the
pump beam, is

$\mathrm{u}(r)=f_{0}\left[\alpha-i\left(\frac{2r^{2}}{W^{4}}-\frac{2(S+1)}{W^{2}}-\sigma+g\right)\right]^{-1},$
(9)

where, as said above, $\eta=1$ is substituted. In the limit of
$r\rightarrow\infty$, Eqs. (9) and (8) carry over into the asymptotically
rigorous expression (6). On the other hand, Eq. (9) predicts a maximum of the
local power at

$r^{2}=\left(r_{\max}^{2}\right)_{\mathrm{TF}}=\left(S+1\right)W^{2}+\frac{\sigma-g}{2}W^{4},$
(10)

cf. Eq. (7). One may expect that the prediction of the local maximum of the
vortex soliton at point (10) is valid when it yields values of
$\left(r_{\max}^{2}\right)_{\mathrm{TF}}$ which are large enough, i.e., if $S$
and $W$ are relatively large. Indeed, the comparison with the numerically
produced profiles of the vortex solitons, displayed below in Figs. 6 and 10,
demonstrates that Eq. (10) predicts, relatively accurately,
$\left(r_{\max}\right)_{\mathrm{TF}}=4$ for $S=5$, $W=2$, $\sigma=-1$, and
$g=0$. However, the prediction given by Eq. (10) is not accurate for $S=1$ and
$2$.

Finally, in the limit of $r\rightarrow 0$, the asymptotic form of the solution
is simple, $u(r,\theta)\approx u_{0}r^{S}$, but constant $u_{0}$ cannot be
found explicitly, as it depends on the global structure of the vortex-soliton
solution. In particular, the crude TF approximation given by Eq. (9) yields
$u_{0}=f_{0}\left[\alpha+i\left(\frac{2(S+1)}{W^{2}}+\sigma-g\right)\right]^{-1}$.

###  II.2 The variational approximation (VA)

A consistent global analytical fit for the vortex solitons may be provided by
VA, based on the Lagrangian of the underlying equation [2]. While Eq. (1),
which includes the linear dissipative term, does not have a Lagrangian
structure, it can be converted into an appropriate form by the substitution
suggested by Ref. [41], which absorbs the dissipative term:

$u\left(r,\theta,t\right)=U(r,t)e^{iS\theta-\alpha t},$ (11)

producing the following time-dependent equation for complex function $U(r,t)$,
where, as said above, we set $\eta=1$ by means of scaling:

$\frac{\partial U}{\partial t}=\frac{i}{2}\left(\frac{\partial^{2}}{\partial
r^{2}}+\frac{1}{r}\frac{\partial}{\partial
r}-\frac{S^{2}}{r^{2}}\right)U+i\sigma\left(\left|U\right|^{2}e^{-2\alpha
t}-1\right)U-ig\left(\left|U\right|^{4}e^{-4\alpha t}-1\right)U+f(r)e^{\alpha
t}.$ (12)

The real Lagrangian which precisely produces the time-dependent equation (12)
is

$\displaystyle L=\int_{0}^{\infty}\left[\frac{i}{2}\left(\frac{\partial
U^{\ast}}{\partial t}U-U^{\ast}\frac{\partial U}{\partial
t}\right)+\frac{1}{2}\left|\frac{\partial U}{\partial
r}\right|^{2}+\left(\frac{S^{2}}{2r^{2}}+\sigma-g\right)|U|^{2}\right.$
$\displaystyle-\left.\frac{\sigma}{2}e^{-2\alpha
t}|U|^{4}+\frac{g}{3}e^{-4\alpha t}|U|^{6}+if(r)e^{\alpha
t}\left(U^{\ast}-U\right)\right]rdr.$ (13)

The simplest ansatz which may be used as the basis for VA follows the form of
the pump term (2):

$U(r,t)=U_{0}e^{i\phi}r^{S}\exp\left(-\frac{r^{2}}{W^{2}}+\alpha t\right),$
(14)

where variational parameters $U_{0}$ and $\phi$ are the real amplitude and
phase shift of the solution with respect to the pump. Power (3) for this
ansatz is

$P_{S}=\pi\Gamma(S+1)\left(\frac{W^{2}}{2}\right)^{S+1}U_{0}^{2},$ (15)

where $\Gamma(S+1)\equiv S!$ is the Gamma-function, and the time-dependent
factors, $\exp\left(\pm\alpha t\right)$, mutually cancel when relations (11)
and (14) are substituted in expression (3). Note that the local power
$|U(r)|^{2}$, corresponding to ansatz (14), attains it maximum at
$r^{2}=SW^{2}/2$.

The substitution of ansatz (14) in Lagrangian (13) and straightforward
integration yields the respective VA Lagrangian,

$\displaystyle L_{\mathrm{VA}}$ $\displaystyle=$ $\displaystyle
U_{0}e^{2\alpha
t}\left\\{\left[6^{-(3S+2)}\Gamma\left(3S+1\right)gW^{2(3S+1)}\right]U_{0}^{5}-\left[2^{-4(S+1)}\Gamma\left(2S+1\right)\sigma
W^{2(2S+1)}\right]U_{0}^{3}\right.$ (16) $\displaystyle+$
$\displaystyle\left.2^{-(S+2)}\Gamma\left(S+1\right)W^{2S}\left[\left(\frac{d\phi}{dt}+(\sigma-g)\right)W^{2}+\left(S+1\right)\right]U_{0}\right.$
$\displaystyle+$
$\displaystyle\left.2^{-(S+1)}\Gamma\left(S+1\right)W^{2(S+1)}f_{0}\sin\phi\right\\}.$

Then, the Euler-Lagrange equations for $U_{0}$ and $\phi$ are obtained as

$\partial L_{\mathrm{VA}}/\partial U_{0}=\delta L_{\mathrm{VA}}/\delta\phi=0,$
(17)

where $\delta/\delta\phi$ stands for the variational derivative. Taking into
regard that Lagrangian (16) must be substituted in the respective action,
$\int L_{\mathrm{VA}}dt$, and then the action must be actually subjected to
the variation, one should apply the time differentiation to factor $e^{2\alpha
t}$ in Lagrangian (16), while deducing the appropriate form of $\delta
L_{\mathrm{VA}}/\delta\phi$ in Eq. (17). Once the Euler-Lagrange equations
(17) have been derived , we consider their stationary (fixed-point) solutions
by setting $dU_{0}/dt=d\phi/dt=0$, which yields

$U_{0}\alpha-f_{0}\cos(\phi)=0,$ (18) $\displaystyle
6^{-(3S+1)}\Gamma\left(3S+1\right)gW^{2(2S+1)}U_{0}^{5}-2^{-2(2S+1)}\Gamma\left(2S+1\right)\sigma
W^{2S+2}U_{0}^{3}$ $\displaystyle+$
$\displaystyle+2^{-(S+1)}\Gamma\left(S+1\right)\left\\{f_{0}W^{2}\sin\phi+\left[(\sigma-g)W^{2}+\left(S+1\right)\right]U_{0}\right\\}$
$\displaystyle=$ $\displaystyle 0.$ (19)

It is relevant to mention that the evolution equation for the power (3), which
follows from Eq. (12), is

$\frac{dP}{dt}=-2\alpha
P+4\pi\int_{0}^{\infty}f(r)\text{Re}\left\\{U(r,t)\right\\}e^{-\alpha t}rdr.$
(20)

The stationary states must satisfy the balance condition, $dP/dt=0$. Then, the
substitution of _ansatz_ (14) and expression (2) for the pump in this
condition yields a simple relation,

$\cos\phi=\frac{\alpha U_{0}}{f_{0}},$ (21)

which is identical to Eq. (18). In particular, Eq. (21) implies that, for the
fixed pump’s amplitude $f_{0}$, the amplitude of the established localized
pattern cannot exceed the maximum value, which corresponds to $\phi=0$ in Eq.
(21):

$U_{0}\leq\left(U_{0}\right)_{\max}=f_{0}/\alpha.$ (22)

###  II.3 VA for the cubic ($g=0$) and quintic ($g\rightarrow\infty$) models

First, we aim to predict stationary states, as solutions of Eqs. (18) and
(19), for the model with cubic-only nonlinearity, i.e., $g=0$, under the
assumption that the loss and pump terms in Eq. (1) may be considered as small
perturbations. In the lowest approximation, i.e., dropping the small term
$\sim f_{0}$ in Eq. (19), one obtains a relatively simple expression which
predicts the squared amplitude of the vortex soliton,

$\left(U_{0}^{2}\right)_{\mathrm{VA}}^{(g=0)}=\frac{2^{2S+1}}{(2S-1)!!}W^{-2S}\left[1+\sigma\left(S+1\right)W^{-2}\right],$
(23)

the respective power (15) of the underlying ansatz (14) being

$P_{\mathrm{VA}}^{(g=0)}=\frac{2^{S}S!\pi}{(2S-1)!!}\left[W^{2}+\sigma\left(S+1\right)\right].$
(24)

Note that expressions (23) and (24) are always meaningful for the self-
focusing sign of the cubic nonlinearity, $\sigma=+1$, while in the case of
defocusing, $\sigma=-1$, the expressions are meaningful if they are positive,
which imposes a restriction on the width of the Gaussian pump: it must be
broad enough, viz.,

$W^{2}>S+1.$ (25)

The dependence of the power given by Eq. (24) on the pump’s squared width
$W^{2}$ for three different values of the vorticity, $S=1,2,3$, is plotted in
Figs. 1(a) and (b) for the self-focusing and defocusing signs of the cubic
term, i.e., $\sigma=+1$ and $-1$, respectively. Note that in Fig. 1(b) for
$\sigma=-1$, there is no solution in the region in which condition (25) does
not hold, hence the VA solution does not exist.

Figure 1: The power of the VA solution, in the case of $g=0$ (no quintic
nonlinearity), vs. the pump’s squared width, as given by Eq. (24), which
neglects weak effects of the pump and loss (small $f_{0}$ and $\alpha$), for
the self-focusing ($\sigma=+1$) in (a) and defocusing ($\sigma=-1$) in (b)
signs of the cubic term. The black curve with circles, the red one with
squares, and the blue one with triangles pertain to vorticities $S=1$, $2$,
and $3$, respectively.

In the limit of the dominant quintic nonlinearity, i.e.,
$g\rightarrow\pm\infty$, opposite to the pure cubic model considered above,
the asymptotic solution of Eq. (19) is

$\left(U_{0}^{2}\right)_{\mathrm{VA}}^{(g\rightarrow\pm\infty)}\approx
2^{S}\sqrt{3^{3S+1}\frac{S!}{(3S)!}}W^{-2S}$ (26)

(in which the large coefficient $g$ cancels out), the respective expression
for power (15) being

$P_{\mathrm{VA}}^{(g\rightarrow\pm\infty)}=\frac{3^{(3S+1)/2}\pi\left(S!\right)^{3/2}}{2\sqrt{(3S)!}}W^{2},$
(27)

cf. Eqs. (23) and (24).

In the following section, we report results of numerical solution of Eq. (1),
comparing them to solutions of the full system of the VA equations (18) and
(19), which include effects of the pump and loss terms.

##  III Numerical results

Simulations of Eq. (1) were conducted by means of the split-step pseudo-
spectral algorithm. The solution procedure started from the zero input, and
was running until convergence to an apparently stable stationary profile (if
this outcome of the evolution was possible). This profile was then compared to
its VA counterpart, produced by a numerical solution of Eqs. (18) and (19)
with the same values of parameters $\alpha$, $\sigma=\pm 1$, $g$, and $f_{0}$,
$W$, $S$ (see Eq. (2)). The results are presented below by varying, severally,
loss $\alpha$, vorticity $S$, the pump’s width $W$ and strength $f_{0}$, and,
eventually, the quintic coefficient $g$. The findings are eventually
summarized in the form of stability charts plotted in Fig. 12.

###  III.1 Variation of the loss parameter $\alpha$

In Fig. 2(a) we display the cross-section (drawn through $y=0$) of the
variational and numerical solutions for the stable vortex solitons obtained
with $\alpha=0.5$, $1.0$, and $2.0$, while the other parameter are fixed as
$\sigma=+1$ (the self-focusing cubic nonlinearity), $g=0$, $f_{0}=1$, $W=2$,
and $S=1$. The accuracy of the VA-predicted solutions presented in Fig. 2(a)
is characterized by the relative power difference from their numerically found
counterparts, which is $5.1\%$, $3.1\%$, and $0.5\%$ for

$\alpha=0.5,1.0,\mathrm{~{}}2.0,$ (28)

respectively. Thus, the VA accuracy improves with the increase of $\alpha$.

Similar results for the self-defocusing nonlinearity, $\sigma=-1$, are
presented in Fig. 2(b), which shows an essentially larger discrepancy between
the VA and numerical solutions, viz., $18.9\%$, $17.4\%$, and $3.4\%$ for the
same set (28) of values of the loss parameter, other coefficients being the
same as in Fig. 2(a). The larger discrepancy is explained by the fact that
localized (bright-soliton) modes are not naturally maintained by the self-
defocusing, hence the ansatz (14), which is natural for the self-trapped
solitons in the case of the self-focusing, is not accurate enough for
$\sigma=-1$. In the same vein, it is natural that, in the latter case, the
discrepancy is more salient for stronger nonlinearity, i.e., smaller $\alpha$,
which makes the respective amplitude higher.

Figure 2: The comparison between cross sections (drawn through $y=0$) of the
VA solutions and their numerically found counterparts (dashed black and solid
blue lines, respectively) for different values of the loss parameter $\alpha$
in Eq. (1), taken from set (28). Panels (a) and (b) pertain to the self-
focusing ($\sigma=+1$) and defocusing ($\sigma=-1$) signs of the cubic
nonlinearity, respectively. The other parameters in Eqs. (1) and ( 2) are
fixed as $g=0$, $\eta=1$, $f_{0}=1$, $W=2$, and $S=1$.

There is a critical value of $\alpha$ below which the vortex solitons are
unstable. As an example, Fig. 3 shows the VA-predicted and numerically
produced solutions for $\alpha=0.2$ and $\sigma=+1$ (the self-focusing
nonlinearity). The observed picture may be understood as a result of the
above-mentioned azimuthal MI which breaks the axial symmetry of the vortex
soliton. More examples of the instability of this type are displayed below.
For the values of other parameters fixed as in Fig. 3, the instability
boundary is $\alpha_{\mathrm{crit}}\approx 0.35$. The stabilizing effect of
the loss at $\alpha>\alpha_{\mathrm{crit}}$ is a natural feature. On the other
hand, the increase of $\alpha$ leads to decrease of the soliton’s amplitude,
as seen in Fig. 2.

In the case of the self-defocusing ($\sigma=-1$), all the numerically found
vortex modes are stable, at least, at $\alpha\geq 0.1$, although the
discrepancy in the values of the power between these solutions and their VA
counterparts is very large at small $\alpha$, exceeding $75\%$ at
$\alpha=0.1$. As mentioned above, the growing discrepancy is explained by the
increase of the soliton’s amplitude with the decrease of $\alpha$. At still
smaller values of $\alpha$, the relaxation of the evolving numerical solution
toward the stationary state is very slow, which makes it difficult to identify
the stability.

Figure 3: Profiles of $|u|^{2}$ produced by VA (a) and numerical solution at
$t=100$ (b) in the case of the self-focusing ($\sigma=+1$) for $\alpha=0.2$.
Other parameter are the same as in Fig. 2(a).

###  III.2 Variation of the pump’s vorticity $S$

To analyze effects of the winding number (vorticity) $S$, we here fix $g=0$
(the pure cubic nonlinearity) and set $\alpha=1$, $f_{0}=1$, $W=2$ in Eqs. (1)
and (2). In the self-focusing case ($\sigma=+1$), the numerically produced
solutions are stable for $S=1$ and $2$, and unstable for $S\geq 3$. In the
former case, the power difference between the VA and numerical solutions is
$3.1\%$ and $4.9\%$ for $S=1$ and $2$, respectively, i.e., the VA remains a
relatively accurate approximation in this case.

In the self-defocusing case ($\sigma=-1$), considering the same values of the
other parameters as used above, the numerical solution produces stable vortex
solitons at least until $S=5$. For the same reason as mentioned above, the
accuracy of VA is much lower for $\sigma=-1$ than for the self-focusing case
($\sigma=+1$), with the respective discrepancies in the power values being
$17.4\%$, $6.7\%$, $24.1\%$, $29.0\%$, and $29.6\%$ for $S=1$, $2$, $3$, $4$,
and $5$, respectively.

For the cogent verification of the stability of the localized vortices in the
case of the self-defocusing, we have also checked it for smallest value of the
loss parameter considered in this work, viz., $\alpha=0.1$, again for $S=1$,
$2$, $3$, $4$, and $5$, and the above-mentioned values of the other
coefficients, i.e., $f_{0}=1$ and $W=2$. Naturally, the discrepancy between
the VA and numerical findings is still higher in this case, being $75\%$,
$72.5\%$, $65.7\%$, $55.2\%$, and $41.1\%$ for $S=1$, $2$, $3$, $4$, and $5$,
respectively. The result is illustrated by Fig. 4 for a relatively large
vorticity, $S=4$. In particular, the pattern of $|u(x,y)|^{2}$ and the
corresponding cross section, displayed in Figs. 4(b) and (c), respectively,
exhibit the established vortex structure and background “garbage" produced by
the evolution.

Figure 4: (a) The VA-predicted pattern and (b) the corresponding result of the
direct simulation of Eq. (1) with $\sigma=-1$ and $g=0$ (cubic self-
defocusing) at $t=100$, initiated by the zero input at $t=0$, for
$\alpha=0.1$, $f_{0}=1$, $W=2$, and vorticity $S=4$ in the pump term (2). (c)
The respective cross sections drawn through $y=0$. (d) The evolution of the
total power $P$ (see Eq. (3) of the numerical solution in the course of the
simulation.

###  III.3 Variation of the pump’s width $W$

To address the effects of the variation of parameter $W$ in Eq. (2), we here
fix $g=0$, $\alpha=1$, $f_{0}=1$, and $S=1$. In Fig. 5(a), the power of the
VA-predicted and numerically found stable vortex-soliton solutions is plotted
as a function of $W$ for both the self-focusing and defocusing cases, i.e.,
$\sigma=+1$ and $\sigma=-1$, respectively. In the former case, the azimuthal
MI sets in at $W\geq 2.7$, see an example in Fig. 5(b) for $W=2.75$.

Figure 5: (a) The power of the VA-predicted vortex-soliton solutions (solid
blue and dashed black lines, pertaining to the self-focusing, $\sigma=+1$, and
self-defocusing, $\sigma=-1$, cubic nonlinearitry, respectively) and their
numerically found counterparts (red circles and green squares pertaining to
the self-focusing and self-defocusing nonlinearitry, respectively) vs. the
pump’s width $W$. Other parameters are $g=0$, $\alpha=1$, $f_{0}=1$, and
$S=1$. (b) The profile produced, at $t=100$, by the numerically generated
unstable solution in the case of the self-focusing, $\sigma=+1$, with
$W=2.75$.

In the self-focusing case, $\sigma=+1$, the azimuthal MI for the solitons with
higher vorticities, $S=2$, $3$, $4$, or $5$, sets in at $W\geq 2.1$, $1.7$,
$1.6$, and $1.4$, respectively. In the self-defocusing case, no
existence/stability boundary was found for the vortex modes with $S=1$, $S=2$,
and $S=3$ (at least, up to $W=5$). At higher values of the vorticity, the
localized vortices do not exist, in the defocusing case, at $W>3.5$ and
$W>2.5$, for $S=4$, and $S=5$, respectively.

###  III.4 Variation of the pump’s strength $f_{0}$

Effects of the variation of $f_{0}$ are reported here, fixing other parameters
as $g=0$, $\alpha=1$, and $W=2$. In the case of the cubic self-focusing,
$\sigma=+1$, the vortex soliton with $S=1$ are subject to MI at $f_{0}\geq
1.6$. As a typical example, in Fig. 6 we display the VA-predicted solution
alongside the result of the numerical simulations for $f_{0}=1.7$. For higher
vorticities, $S=2$, $3$, $4$, and $5$, the azimuthal instability sets in at
$f_{0}\geq 1.1$, $0.6$, $0.3$, and $0.08$, respectively. Naturally, the narrow
vortex rings with large values of $S$ are much more vulnerable to the quasi-
one-dimensional azimuthal MI.

The power of the vortex solitons with $S=1$ and $2$, as produced by the VA and
numerical solution, is plotted vs. the pump amplitude $f_{0}$ in Fig. 7(a). As
an example, Fig. 7(b) showcases an example of the cross-section profile of the
vortex soliton with $S=2$, demonstrating the reliability of the VA prediction.
In the range of $f_{0}\leq 1$, the highest relative difference in the power
between the numerical and variational solutions cases is $5.5\%$ and $7.5\%$
for $S=1$ and $S=2$, respectively.

Figure 6: (a) The VA-predicted profile of $|u|^{2}$ in the self-focusing case,
with parameters $g=0$, $\sigma=+1$, $\alpha=1$, $f_{0}=1.7$, $W=2$, and $S=1$.
(b) The unstable solution, produced, at $t=100$, by the simulations of Eq. (1)
for the same parameters.

Figure 7: (a) The power versus $f_{0}$ for the confined vortex modes in the
self-focusing case ($\sigma=+1$). The VA solutions for $S=1$ and $2$ are shown
by solid blue and dashed black lines, respectively. The corresponding
numerical solutions are represented by red circles and green squares,
respectively. Recall that the numerical solutions are stable, in this case, at
$f_{0}<1.6$ and $f_{0}<1.1$, for $S=1$ and $2$, respectively. (b) The VA-
predicted and numerically obtained (the dashed red and solid blue lines,
respectively) profiles of the stable solution with $S=2$ and $f_{0}=0.4$,
drawn as cross sections through $y=0$ . The other parameters are $g=0$,
$\sigma=+1$, $\alpha=1$, $W=2$.

It is relevant to mention that the “traditional” azimuthal instability of
vortex-ring solitons with winding number $S$ demonstrates fission of the
original axially symmetric shape into a set consisting of a large number
$N\geq S$ of symmetrically placed localized fragments [40, 2], while the above
examples, displayed in Figs. 3(b), 5(b), and 6(b), demonstrate the appearance
of a single bright fragment and a “garbage cloud” distributed along the
original ring. At larger values of $f_{0}$, our simulations produce examples
of the “clean” fragmentation, viz., with $N=4$ produced by the unstable vortex
rings with $S=1$ in Fig. 8(a), $N=5$ produced by $S=2$ and $3$ in (b) and (d),
$N=7$ by $S=2$ and $3$ in (c) and (e), and $N=8$ by $S=3$ in (f). These
outcomes of the instability development are observed at the same evolution
time $t=100$ as in Figs. 3(b), 5(b), and 6(b). The gradual increase of the
number of the fragments on $S$ is explained by the dependence of the azimuthal
index of the fastest growing eigenmode of the breaking instability on the
underlying winding number $S$, which is a generic property of vortex solitons
[40, 2].

Figure 8: Examples of the fission of unstable vortex-ring solitons produced by
simulations of the LL equation (1) with $\sigma=+1$ (cubic self-focusing),
$g=0$ (no quintic nonlinearity), $\alpha=2$, $W=2$, and $\eta=1$. Each plot
displays the result of the numerical simulations at time $t=100$. Values of
the initial vorticity and pump’s strength are indicated in panels.

In self-defocusing case, $\sigma=-1$, a summary of the results produced by the
VA and numerical solution for the stable vortex solitons with $S=1$ and $2$,
in the form of the dependence of their power on $f_{0}$, is produced in Fig.
9(a) (cf. Fig. 7(a) for $\sigma=+1$). Naturally, the VA-numerical discrepancy
increases with the growth of the pump’s strength, $f_{0}$, see an example in
Fig. 9(b). Unlike the case of $\sigma=+1$, in the case of the self-defocusing
the vortex modes with $S\leq 5$ remain stable, at least, up to $f_{0}=5$
(here, we do not consider the case of $S>5$).

Figure 9: (a) The power versus $f_{0}$ for the vortex modes in the self-
defocusing case ($\sigma=-1$). The VA solutions for $S=1$ and $2$ are shown by
solid blue and dashed black lines, respectively. The corresponding numerical
solutions are presented by red circles and green squares, respectively. (b)
The VA-predicted and numerically obtained (the dashed red and solid blue
lines, respectively) profiles of the solution with $S=2$ and $f_{0}=2$, drawn
as the cross sections through $y=0$ . The other parameters are $g=0$,
$\sigma=-1$, $\alpha=1$, $W=2$.

###  III.5 Influence of the quintic coefficient $g$

In the above analysis, the quintic term was dropped in the LL equation (1),
setting $g=0$. To examine the impact of this term, we first address the case
shown above in Fig. 3, which demonstrated that the vortex soliton with $S=1$,
as a solution to Eqs. (1) and (2) with $g=0$, $\sigma=+1$ and $f_{0}=1$,
$W=2$, was unstable if the loss parameter fell below the critical value,
$\alpha=0.35$. We have found that, adding to Eq. (1) the quintic term with
either $g=-1$ or $g=+1$ (the self-focusing or defocusing quintic nonlinearity,
respectively) leads to the _stabilization_ of the vortex mode displayed in
Fig. 3, which was unstable in the absence of the quintic term. The
stabilization of the soliton by the quintic self-defocusing is a natural fact.
More surprising is the possibility to provide the stabilization by the self-
focusing quintic nonlinearity because, in most cases, the inclusion of such a
term gives rise to the supercritical collapse in 2D, making all solitons
strongly unstable [40, 2]. However, it is concluded from the stability charts
displayed below in Fig. 12 that the stabilizing effect of the quintic self-
focusing occurs only at moderately small powers, for which the the quintic
term is not a clearly dominant one. In the general case, solitons stability
regions naturally shrink in Figs. 12 under the action of the quintic self-
focusing quintic term, with $g<0$.

In the presence of the quintic term, the comparison of the numerically found
stabilized vortex-soliton profiles with their VA counterparts, whose
parameters are produced by a numerical solution of Eqs. (18) and (19), is
presented in Fig. 10. Similar to the results for the LL equation with the
cubic-only nonlinearity, the VA is essentially more accurate in the case of
the self-focusing sign of the quintic term ($g<0$) than in the opposite case,
$g>0$. In particular, in the case shown in Fig. 10, the power-measured
discrepancy for $g=-1$ and $+1$ is, respectively, $4\%$ and $64\%$. An
explanation for this observation is provided by the fact that the soliton’s
amplitude is much higher in the latter case.

Figure 10: Cross-section profiles (drawn through $y=0$) for stable vortex
solitns with $S=1$, produced by Eq. (1) with the quintic self-focusing, $g=-1$
(a) or defocusing, $g=+1$ (b) term. In both cases, the cubic self-focusing
cubic term, with $\sigma=+1$, is present. The numerically found profiles and
their VA-produced counterparts are displayed, respectively, by the solid blue
and dashed lines. Other parameter are $\alpha=0.2$, $\eta=1$, and $f_{0}=1$,
$W=2$.

Another noteworthy finding is the stabilization of higher-vorticity solitons
by the quintic term. For instance, it was shown above that, for parameters
$\alpha=1$, $\eta=1$, $f_{0}=1$, $W=2$, and $\sigma=+1$ (the cubic self-
focusing) all vortex solitons with $S\geq 3$, produced by Eq. (1), were
unstable. Now, we demonstrate that the soliton with $S=3$ is stabilized by
adding the self-defocusing quintic term with a small coefficient, just
$g=0.1$, see Fig. 11(b). As a counter-intuitive effect, the stabilization of
the same soliton by the self-focusing quintic term is possible too, but the
necessary coefficient is large, $g=-6$, see Fig. 11(a) (recall that $g=-1$ is
sufficient for the stabilization of the vortex soliton with $S=1$ and
$\alpha=0.2$ in Fig. 10(a)). Nevertheless, similar to what is said above, the
set of Figs. 12(c,f,i) demonstrates the natural shrinkage of the stability
area under the action of the quintic self-focusing. For the stabilized vortex
modes shown in Fig. 11(a), in the case when both the cubic and quintic terms
are self-focusing, the relative power-measured discrepancy between the
numerical and VA-predicted solutions is very small, $\approx 0.7\%$, while in
the presence of the weak quintic self-defocusing in Fig. 11(b) the discrepancy
is $6.6\%$.

Figure 11: Cross-section profiles (drawn through $y=0$) of vortex solitons
with $S=3$, stabilized by the self-focusing (a) or defocusing (b) quintic
tterm in Eq. (1), with the respective coefficient $g=-6$ or $g=0.1$, other
parameters being $\sigma=-1$ (cubic self-focusing), $\alpha=1$, $\eta=1$, and
$f_{0}=1$, $W=2$ in Eq.( 2). The numerically found solutions and their VA-
predicted counterparts are plotted by the solid blue and dashed red lines,
respectively.

###  III.6 Stability charts in the parameter space

The numerical results produced in this work are summarized in the form of
stability areas plotted in Fig. 12 in the parameter plane
$\left(f_{0},\alpha\right)$, for the vortex-soliton families with winding
numbers $S=1$, $2$, $3$, and three values of the quintic coefficient, $g=-1$,
$0$, $+1$, while the cubic term is self-focusing, $\sigma=+1$, and the width
of the pump beam is fixed, $W=2$. In addition to that, stability charts
corresponding to the combination of the cubic self-defocusing ($\sigma=-1$)
and quintic focusing ($g=-1$), also for $S=1,2,3$, are plotted in Fig. 13.

The choice of the parameter plane $\left(f_{0},\alpha\right)$ in the stability
diagrams displayed in Fig. 12 is relevant, as the strength of the pump beam,
$f_{0}$, and loss coefficient, $\alpha$,are amenable to accurate adjustment in
the experiment (in particular, $\alpha$ may be tuned by partially compensating
the background loss of the optical cavity by a spatially uniform pump, taken
separately from the confined pump beam). As seen in all panels of Figs. 12 and
13, the increase of $\alpha$ naturally provides effective stabilization of the
vortex modes, while none of them may be stable at $\alpha=0$, in agreement
with the known properties of vortex-soliton solutions of the 2D nonlinear
Schrödinger equation with the cubic and/or cubic-quintic nonlinearity [40, 2].
The apparent destabilization of the vortices with the increase of the pump’s
amplitude $f_{0}$ is explained by the ensuing enhancement of the destabilizing
nonlinearity. Other natural features exhibited by Figs. 12 are the general
stabilizing/destabilizing effect of the quintic self-defocusing/focusing (as
discussed above), and expansion of the splitting-instability area with the
increase of $S$ (the latter feature is also exhibited by Fig. 13). The latter
finding is natural too, as larger $S$ makes the ring-shaped mode closer to the
quasi-1D shape (see, in particular, Fig. 8), which facilitates the onset of
the above-mentioned azimuthal MI (modulational instability).

Figure 12: Stability areas for families of the vortex solitons with winding
numbers $S=1$, $2,$ and $3$, in the plain of the loss coefficient ($\alpha$)
and amplitude of the pump beam ($f_{0}$), for three different values of the
quintic coefficient, $g=-1,0,+1$ (recall that $g<0$ and $g>0$ correspond,
respectively, to the self-focusing and defocusing). Other parameters of Eqs.
(1) and (2) are $\sigma=+1$ (the self-focusing cubic term), $\eta=1$, and
$W=2$. Figure 13: The same as in Figs. 12(a-c), but for $\sigma=-1$ and $g=-1$
(the cubic self-defocusing and quintic focusing terms).

Thus, the inference is that the instability mode which determines the boundary
of the stability areas in Figs. 12 and 13 is the breaking of the axial
symmetry of the vortex rings by azimuthal perturbations, as shown above, in
particular, in Figs. 3(b), 5(b), and 6(b). The destabilization through the
spontaneous splitting of the rings into symmetric sets of fragments (see Fig.
8) occurs deeply inside the instability area, i.e., at larger values of
$f_{0}$.

##  IV Conclusion

We have introduced the two-dimensional LL (Lugiato-Lefever) equation including
the self-focusing or defocusing cubic or cubic-quintic nonlinearity and the
confined pump with embedded vorticity (winding number), $S\leq 5$. Stable
states in the form of vortex solitons (rings) for these values of $S$ are
obtained, in parallel, in the semi-analytical form by means of the VA
(variational approximation) and numerically, by means of systematic
simulations of the LL equation starting from the zero input. The VA provides
much more accurate results in the case of the self-focusing nonlinearity than
for the defocusing system. Stability areas of the vortex solitons with
$S=1,2,3$ are identified in the plane of experimentally relevant parameters,
viz., the pump amplitude and loss coefficient, for the self-focusing and
defocusing signs of the cubic and quintic terms. Stability boundaries for the
vortex rings are determined by the onset of the azimuthal instability which
breaks their axial symmetry. These findings suggest new possibilities for the
design of tightly confined robust optical modes, such as vortex pixels.

As an extension of this work, it may be interesting to construct solutions
pinned to a symmetric pair of pump beams with or without intrinsic vorticity.
In this context, it is possible to consider the beam pair with identical or
opposite vorticities. In the case of the self-focusing sign of the
nonlinearity, one may expect onset of spontaneous breaking of the symmetry in
the dual-pump configuration. Results for this setup will be reported
elsewhere.

## Acknowledgments

We thank Prof. Branko Dragovich for invitation to submit the paper to the
Special Issue of Symmetry on the topic of “Selected Papers on Nonlinear
Dynamics".

The work of S.K. and B.A.M. is supported, in part, by the Israel Science
Foundation through grant No. 1695/22. W.B.C. acknowledges the financial
support of the Brazilian agency CNPq (grant #306105/2022-5). This work was
also performed as a part of program #465469/2014-0 of the Brazilian National
Institute of Science and Technology (INCT) for Quantum Information.

## References

  * [1] Y. S. Kivshar and G. Agrawal, Optical Solitons: From Fibers to Photonic Crystals (Elsevier Science, 2003).
  * [2] B. A. Malomed, Multidimensional Solitons (American Institute of Physics Publishing, Melville, NY, 2022).
  * [3] N. N. Rosanov, Spatial Hysteresis and Optical Patterns (Springer: Berlin, 2002).
  * [4] Dissipative Optical Solitons, ed. by M. F. S. Ferreira (Springer Nature Switzerland AG, Cham, Switzerland, 2022).
  * [5] P. Grelu and N. Akhmediev, Dissipative solitons for mode-locked lasers, Nature Phot. 6, 84-92 (2012).
  * [6] M. X. Jiang, X. N. Wang, Q. Ouyang, and H. Zhang, Spatiotemporal chaos control with a target wave in the complex Ginzburg-Landau equation system, Phys. Rev. E 69, 56202 (2004).
  * [7] L. A. Lugiato and R. Lefever, Spatial Dissipative Structures in Passive Optical Systems, Phys. Rev. Lett. 58, 2209 (1987).
  * [8] M. Tlidi and K. Panajotov, Two-dimensional dissipative rogue waves due to time-delayed feedback in cavity nonlinear optics, Chaos 27, 013119 (2017).
  * [9] K. Panajotov, M. G. Clerc, and M. Tlidi, Spatiotemporal chaos and two-dimensional dissipative rogue waves in Lugiato-Lefever model, Eur. Phys. J. B 71, 176 (2017).
  * [10] M. Tlidi and M. Taki, Rogue waves in nonlinear optics, Adv. Opt. Photon. 14, 87-147 (2022).
  * [11] Y. Sun, P. Parra-Rivas, F. Mangini, and S. Wabnitz, Multidimensional localized states in externally driven Kerr cavities with a parabolic spatiotemporal potential: a dimensional connection, arXiv:2401.15689.
  * [12] G. J. de Valcárcel and K. Staliunas, Phase-bistable Kerr cavity solitons and patterns, Phys. Rev. A 87, 043802 (2013).
  * [13] S. Coen, H. G. Randle, T. Sylvestre, and M. Erkintalo, Modeling of octave-spanning Kerr frequency combs using a generalized mean-field Lugiato-Lefever model, Opt. Lett. 38, 37-39 (2013).
  * [14] M. R. E. Lamont, Y. Okawachi, and A. L. Gaeta, Route to stabilized ultrabroadband microresonator-based frequency combs, Opt. Lett. 38, 3478-3481 (2013).
  * [15] C. Godey, I. V. Balakireva, A. Coillet, and Y. K. Chembo, Stability analysis of the spatiotemporal Lugiato-Lefever model for Kerr optical frequency combs in the anomalous and normal dispersion regimes, Phys. Rev. A 89, 063814 (2014).
  * [16] V. E. Lobanov, G. Lihachev, T. J. Kippenberg, and M. L. Gorodetsky, Frequency combs and platicons in optical microresonators with normal GVD, Opt. Exp. 23, 7713-7721 (2015).
  * [17] M. Karpov, H. Guo, A. Kordts, V. Brasch, M. H. P. Pfeiffer, M. Zervas, M. Geiselmann, and T. J. Kippenberg, Raman Self-Frequency Shift of Dissipative Kerr Solitons in an Optical Microresonator, Phys. Rev. Lett. 116, 103902 (2016).
  * [18] F. Copie, M. Conforti, A. Kudlinski, A. Mussot, and S. Trillo, Competing Turing and Faraday Instabilities in longitudinally modulated passive resonators, Phys. Rev. Lett. 116, 143901 (2016).
  * [19] P. Parra-Rivas, D. Gomila, P. Colet, and L. Gelens, Interaction of solitons and the formation of bound states in the generalized Lugiato-Lefever equation, Eur. Phys. J. D 71, 198 (2017).
  * [20] M. G. Clerc, M. A. Ferré, S. Coulibaly, R. G. Rojas, and M. Tlidi, Chimera-like states in an array of coupled-waveguide resonators, Opt. Lett. 42, 2906-2909 (2017).
  * [21] B. Garbin, Y. Wang, S. G. Murdoch, G.-L. Oppo, S. Coen, and M. Erkintalo, Experimental and numerical investigations of switching wave dynamics in a normally dispersive fibre ring resonator, Eur. Phys. J. D 71, 240 (2017).
  * [22] Q. Li, T. C. Briles, D. A. Westly, T. E. Drake, J. R. Stone, B. R. Ilic, S. A. Diddams, S. B. Papp, and K. Srinivasan, Stably accessing octave-spanning microresonator frequency combs in the soliton regime, Optica 4, 193-203 (2017).
  * [23] L. A. Lugiato, F. Prati, M. L. Gorodetsky, and T. J. Kippenberg, From the Lugiato–Lefever equation to microresonator-based soliton Kerr frequency combs, Phil. Trans. R. Soc. A. 376, 20180113 (2018).
  * [24] X. Dong, C. Spiess, V. G. Bucklew, and W. H. Renninger, Chirped-pulsed Kerr solitons in the Lugiato-Lefever equation with spectral filtering, Phys. Rev. Research 3, 033252 (2021).
  * [25] S.-W. Huang, J. Yang, S.-H. Yang, M. Yu, D.-L. Kwong, T. Zelevinsky, M. Jarrahi, and C. W. Wong, Globally stable microresonator Turing pattern formation for coherent high-power THz radiation on-chip, Phys. Rev. X 7, 041002 (2017).
  * [26] Y. V. Kartashov, O. Alexander, and D. V. Skryabin, Multistability and coexisting soliton combs in ring resonators: the Lugiato-Lefever approach, Opt. Express 25, 11550-11555 (2017).
  * [27] A. Coillet, I. Balakireva, R. Henriet, K. Saleh, L. Larger, J. M. Dudley, C. R. Menyuk, Y. K. Chembo, Azimuthal Turing patterns, bright and dark cavity solitons in Kerr combs generated with whispering-gallery-mode resonators, IEEE Photon. J. 5, 6100409-6100409 (2013).
  * [28] Y. K. Chembo and C. R. Menyuk, Spatiotemporal Lugiato-Lefever formalism for Kerr-comb generation in whispering-gallery-mode resonators, Phys. Rev. A 87, 053852 (2013).
  * [29] H. Taheri, A. A. Eftekhar, K. Wiesenfeld, and A. Adibi, Soliton Formation in Whispering-Gallery-Mode Resonators via Input Phase Modulation, IEEE Photon. J. 7, 2200309 (2015).
  * [30] Y.-Y. Wang, M.-M. Li, G.-Q. Zhou, Y. Fan, and X.-J. Lai, Rotating vortex-like soliton in a whispering gallery mode microresonator, Eur. Phys. J. Plus 134, 161 (2019).
  * [31] T. Daugey, C. Billet, J. Dudley, J.-M. Merolla, and Y. K. Chembo, Kerr optical frequency comb generation using whispering-gallery-mode resonators in the pulsed-pump regime, Phys. Rev. A 103, 023521 (2021).
  * [32] Q.-H. Cao, K.-L. Geng, B-W. Zhu, Y.-Y. Wang, and C.-Q. Dai, Scalar vortex solitons and vector dipole solitons in whispering gallery mode optical microresonators, Chaos, Sol. & Fract. 166, 112895 (2023).
  * [33] W. B. Cardoso, L. Salasnich, and B. A. Malomed, Localized solutions of Lugiato-Lefever equations with focused pump, Sci. Rep. 7, 16876 (2017).
  * [34] W. B. Cardoso, L. Salasnich, and B. A. Malomed, Zero-dimensional limit of the two-dimensional Lugiato-Lefever equation. Eur. Phys. J. D 71, 112 (2017).
  * [35] M. Quiroga-Teixeiro and H. Michinel, Stable azimuthal stationary state in quintic nonlinear optical media, J. Opt. Soc. Am. B 14, 2004-2009 (1997).
  * [36] G. Boudebs, S. Cherukulappurath, H. Leblond, J. Troles, F. Smektala, and F. Sanchez, Experimental and theoretical study of higher-order nonlinearities in chalcogenide glasses, Opt. Commun. 219, 427-432 (2003).
  * [37] A. S. Reyna and C. B. de Araújo, High-order optical nonlinearities in plasmonic nanocomposites – a review, Adv. Opt. Phot. 9, 720-774 (2017).
  * [38] D. L. Andrews, symmetry and quantum features in optical vortices, Symmetry 13, 1368 (2021).
  * [39] A. Ramaniuk, N. V. Hung, M. Giersig, K. Kempa, V. V. Konotop, and M. Trippenbach, Vortex creation without stirring in coupled ring resonators with gain and loss, Symmetry 10, 195 (2018).
  * [40] B. A. Malomed, (INVITED) Vortex solitons: Old results and new perspectives, Physica D 399, 108-137 (2019).
  * [41] R. K. Bullough, A. P. Fordy, and S. V. Manakov, Adiabatic invariants theory of near-integrable systems with damping, Phys. Lett. A 91, 98-100 (1982).