# Nonlinear photonic disclination states

Boquan Ren  Key Laboratory for Physical Electronics and Devices of the
Ministry of Education & Shaanxi Key Lab of Information Photonic Technique,
School of Electronic Science and Engineering, Faculty of Electronic and
Information Engineering, Xi’an Jiaotong University, Xi’an 710049, China
Hongguang Wang  Key Laboratory for Physical Electronics and Devices of the
Ministry of Education & Shaanxi Key Lab of Information Photonic Technique,
School of Electronic Science and Engineering, Faculty of Electronic and
Information Engineering, Xi’an Jiaotong University, Xi’an 710049, China
Yaroslav V. Kartashov  Institute of Spectroscopy, Russian Academy of Sciences,
Troitsk, Moscow, 108840, Russia  Yongdong Li  Key Laboratory for Physical
Electronics and Devices of the Ministry of Education & Shaanxi Key Lab of
Information Photonic Technique, School of Electronic Science and Engineering,
Faculty of Electronic and Information Engineering, Xi’an Jiaotong University,
Xi’an 710049, China  Yiqi Zhang<EMAIL_ADDRESS>Key Laboratory for
Physical Electronics and Devices of the Ministry of Education & Shaanxi Key
Lab of Information Photonic Technique, School of Electronic Science and
Engineering, Faculty of Electronic and Information Engineering, Xi’an Jiaotong
University, Xi’an 710049, China

###### Abstract

Higher-order topological insulators are unusual materials that can support
topologically protected states, whose dimensionality is lower than the
dimensionality of the structure at least by 2. Among the most intriguing
examples of such states are zero-dimensional corner modes existing in two-
dimensional higher-order insulators. In contrast to corner states, recently
discovered disclination states also belong to the class of higher-order
topological states, but are bound to the boundary of the disclination defect
of the higher-order topological insulator and can be predicted using the bulk-
disclination correspondence principle. Here, we present the first example of
the nonlinear photonic disclination state bifurcating from its linear
counterpart in the disclination lattice with a pentagonal or heptagonal core.
We show that nonlinearity allows to tune location of the disclination states
in the bandgap and notably affects their shapes. The structure of the
disclination lattice is crucial for stability of these nonlinear topological
states: for example, disclination states are stable in the heptagonal lattice
and are unstable nearly in the entire gap of the pentagonal lattice. Nonlinear
disclination states reported here are thresholdless and can be excited even at
low powers. Nonlinear zero-energy states coexisting in these structures with
disclination states are also studied. Our results suggest that disclination
lattices can be used in the design of various nonlinear topological functional
devices, while disclination states supported by them may play an important
role in applications, where strong field confinement together with topological
protection are important, such as the design of topological lasers and
enhancement of generation of high harmonics.

##  I Introduction

Topological insulators that behave as conventional insulators in the bulk, but
allow topologically protected currents at their edges are frequently
considered as a new phase of matter [1, 2]. Originating from solid state
physics, the concept of topological insulators nowadays extends to many other
areas, such as photonics [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14], acoustics
[15, 16, 17, 18], physics of atomic systems [19], mechanics [20], and physics
of polariton condensates in microcavities [21, 22, 23], to name just a few.
After experimental realization of the first photonic topological insulator [4]
and first Floquet photonic insulator [5], the interest to rich physics of
photonic topological systems has grown dramatically and nowadays one can talk
about the birth of a new area – topological photonics [24, 25]. Among the
reasons of the blooming of this research direction are not only unusual
physical properties of photonic topological systems, but also richness of
their potential practical applications, such as the development of robust
topological lasers [26, 27, 28, 29, 30, 31, 32], telecommunications and robust
information transfer [33, 34], and design of on-chip topological devices [35,
36].

Formation of topological edge states is usually connected with breakup of
certain symmetries of the Hamiltonian of the system. For example, systems with
broken time-reversal symmetry can support unidirectional edge states, while in
valley-Hall systems topological edge states appear, when inversion symmetry of
the structure is broken [24, 25]. Particular attention is now paid to so-
called higher-order topological insulator (HOTI) phases that usually can be
realized by proper shift of the sites in the unit cell of the structure [37,
38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 52, 53], and whose
distinctive feature is the existence of topologically protected states with
dimensionality at least by 2 lower than the dimensionality of the underlying
structure. For example, two-dimensional (2D) HOTI along with 1D topological
edge states supports also 0D corner modes [54]. By analogy with bulk-boundary
correspondence principle explaining the emergence of the edge states in
conventional topological insulators [25], higher-order bulk-boundary
correspondence was proposed to explain the emergence of the corner states in
HOTIs, which suggests that the corner state is a result of the filling anomaly
associated with the fractional charge [55]. Corner states in HOTIs are
characterized by strong concentration of the light field, beneficial for the
enhancement of nonlinear processes [56] and realization of low-threshold
lasing [57, 58, 59, 60], as well as for construction of high-Q topological
cavities [61, 62].

It was proposed that disclinations in topological crystalline insulators may
lead to the emergence of very unusual type of topological states [63, 64, 65,
66, 67, 68]. The disclination that serves as a crystallographic defect,
disrupts the lattice structure and strongly traps fractional charges,
alongside which higher-order topological states can be predicted. The bulk-
disclination correspondence proposed for such systems illustrates the
importance of fractional charges as a probe of crystalline topology. Very
recently, fractional charges as well as the formation of the higher-order
topological states due to the disclination defects were observed
experimentally in linear systems [69, 70, 71]. These higher-order topological
states are disclination states and they typically localize at the boundary of
the hollow disclination core of the structure. In addition, disclination
defects are useful for vortex state generation [72] and may act as a domain
wall hosting photonic topological edge states [73]. If the disclination defect
does not break the chiral symmetry of the system, midgap zero-energy states
can also form in the system, which are confined to the disclination core [74]
even though there is no fractional charge.

It is well known that one of the advantages of the photonic platform is the
possibility of strong nonlinear response of the system. However, the impact of
nonlinearity on topological systems with disclinations and on topological
disclination states, in particular, have never been explored. At the same
time, the nonlinearity plays a significant role in control and switching of
topological states [75], and brings a variety of intriguing phenomena,
including modulational instability [76] and bistability of the unidirectional
edge states [77], topological transitions [78, 10], formation of topological
solitons [79, 80, 81, 82, 83, 84], and breakdown of topological transport [85,
86]. It is thus important to understand how nonlinear response of the system
can change the properties, localization, and conditions of formation of
topological disclination states.

Figure 1: Disclination lattices and corresponding linear spectra. (a)
Illustration of the principle of construction of the disclination lattices
with pentagonal (b) and heptagonal (c) cores from the honeycomb lattice by
removing or inserting a sector with the Frank angle $\Omega=2\pi/6$,
respectively. Two magnified unit cells displayed above panel (a) indicate
spacing $\ell_{\rm intra}$ and $\ell_{\rm inter}$ between the waveguides when
shift is introduced, while the first Brillouin zone of the original honeycomb
structure with high-symmetric points $(\Gamma,\,\rm K,\,M)$ is shown below
panel (a). The quantity $r=\ell_{\rm intra}/\ell_{\rm inter}$ characterizes
the separation between sites in the unit cell. (d) Spectrum of the
disclination lattice with pentagonal core versus $r$ in which the critical
value $r_{c}=1$ is indicated by the vertical dashed line. Wannier centers (red
dots) are also shown, which are at the center of the unit cell if $r<r_{c}$
and at the edges of the unit cell if $r>r_{c}$. (e) Magnified band structure
with red and green lines indicating disclination and zero-energy states,
respectively. (f, g) Spectrum and its magnified version for the lattice with a
heptagonal core. Figure 2: Examples of the disclination states. (a)
Eigenvalues of the modes of the disclination lattice with a pentagonal core at
$r=1.92$. Yellow, cyan, blue, and red dots correspond to the corner, edge,
zero-energy, and disclination states, respectively. Black dots are bulk
states. (b) Field modulus distributions of exemplary states shown within the
window $-35\leq x,y\leq 35$. Deformed lattice cells are shown with yellow
lines, while lattice sites are indicated by white circles. (c, d) Eigenvalues
of modes and examples of their profiles for the lattice with a heptagonal core
at $r=1.92$. The insets in the bottom panels show phase distributions in the
disclination states (red color represents $\pi$ and blue color $0$) within the
window $-10\leq x,y\leq 10$.

In this work, for the first time to our knowledge, we report on the formation
and stability of the nonlinear disclination states in photonic structures with
different symmetries. Based on the honeycomb lattice, we construct two
different photonic lattices with a pentagonal disclination core and a
heptagonal disclination core by removing or inserting a sector with a Frank
angle of $2\pi/6$, respectively. The chiral symmetry of both these
disclination lattices is broken and fractional charge will appear in the unit
cell around the disclination core if they are in the topological insulator
phase. We present spectra of modes of such disclination lattices and obtain
disclination states in the topological phase. We found that in the presence of
the focusing nonlinearity of the material the nonlinear thresholdless
disclination states can bifurcate from linear ones and that their propagation
constants can shift under the action of nonlinearity into the gap or even
further into allowed band leading to notable modification of the shape of the
state. Linear stability analysis predicts important differences in stability
properties of the nonlinear disclination states in lattices with pentagonal
and heptagonal cores. Importantly, all our results are obtained in the frames
of the continuous model (nonlinear Schrödinger-like equation describing
paraxial light propagation) taking into account all details of shallow
refractive index landscape and couplings between all waveguides forming the
lattice. The nonlinear disclination state opens up many potentials for deep
exploration of the behavior and applications of higher-order topological
states.

Figure 3: Nonlinear disclination and zero-energy state families. (a) Nonlinear
families of topological states in the disclination lattice with a pentagonal
core at $r=1.82$. Red line corresponds to the disclination state, while blue
line corresponds to the zero-energy one. Solid lines correspond to stable
branches, while dashed lines correspond to unstable ones. Gray regions
represent the bulk band (right) and edge band (left), while vertical dashed
line shows band of corner states. (b) Field modulus distributions at different
propagation constant values corresponding to the dots numbered $1\sim 4$ in
panel (a). Panels in the top row show disclination states and those in the
bottom row show zero-energy states. (c, d) Setup as in (a, b), but for the
disclination lattice with a heptagonal core at $r=1.79$.

##  II Disclination lattice, fractional disclination charge and disclination
state

We consider paraxial propagation of a light beam in the material with
imprinted waveguide array and focusing cubic nonlinearity that can be
described by the nonlinear Schrödinger-like equation for the dimensionless
light field amplitude $\psi$:

$i\frac{\partial\psi}{\partial
z}=-\frac{1}{2}\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}\right)\psi-\mathcal{R}(x,y)\psi-|\psi|^{2}\psi.$ (1)

Here, $x,y$ are the transverse coordinates normalized to the characteristic
transverse scale $r_{0}=10\,\mu\rm m$, $z$ is the propagation distance
normalized to the diffraction length $kr_{0}^{2}$, where ${k=2\pi n/\lambda}$
is the wavenumber, ${\lambda=800\,\rm nm}$ is the wavelength, $n=1.45$ is the
refractive index of the material, the function ${\mathcal{R}}(x,y)$ describes
the waveguide array and it can be written as

${\mathcal{R}}(x,y)=p\sum_{m,n}e^{-[(x-x_{m})^{2}+(y-y_{n})^{2}]/\sigma^{2}},$
(2)

where $(x_{m},y_{n})$ are the coordinates of the lattice “sites” (i.e.
waveguides), $\sigma$ is the waveguide width, and $p=k^{2}r_{0}^{2}\delta n/n$
is the depth of the lattice that is proportional to the refractive index
modulation depth $\delta n$. We assume here that the structure that we
consider is fabricated using the fs-laser direct writing technique in fused
silica [87, 88]. The original finite honeycomb lattice that is used for the
construction of disclination structures is displayed in Fig. 1(a). The
horizontal separation between centers of two neighboring honeycomb cells
(indicated by solid hexagons) in this non-deformed structure is $\sqrt{3}a$,
where $a$ determines the size of the unit cell as marked in Fig. 1(a). We
further use the parameters $a=3.2$ (corresponding to $32\,\mu\rm m$),
$\sigma=0.5$ (corresponding to $5\,\mu\rm m$) and ${p=10}$ (${\delta n\sim
1.12\times 10^{-3}}$) representative for fs-laser written waveguide arrays
[89, 90]. In our normalizations dimensionless propagation distance $z=1$
corresponds to $\sim 1.14~{}\textrm{mm}$, while dimensionless intensity
$|\psi|^{2}=1$ corresponds to peak intensity of
$I=n|\psi|^{2}/k^{2}r_{0}^{2}n_{2}\approx 5\times 10^{15}~{}\textrm{W/m}^{2}$.
Similar structures can be also optically induced in photorefractive crystals
[91, 92, 84], where the nonlinearity is saturable. Notice that continuous
nonlinear Schrödinger equation (1) takes into account all details of the
refractive index landscape created in the material and on this reason it
accounts for coupling between all waveguides in the structure as well as for
dynamics of modal fields inside individual waveguides.

In each unit cell of the lattice there are six sites which form a hexagon with
the center in the center of the unit cell. The spacing between sites in the
unit cell can be increased or reduced. This gives rise to the structure with
separation between two neighboring sites in the unit cell denoted as
$\ell_{\rm intra}$ and horizontal separation between two nearest sites from
different cells denoted as $\ell_{\rm inter}$ [93, 94], as illustrated in Fig.
1(a). The case $\ell_{\rm intra}=\ell_{\rm inter}=a$ corresponds to the ideal
non-deformed honeycomb structure. Further, by removing a sector with the Frank
angle $\Omega=2\pi/6$ that is shaded in Fig. 1(a), and gluing the remaining
nodes of the model, i.e. filling the remaining gap by uniformly stretching the
lattice in the circular direction, one obtains a finite disclination lattice
with a pentagonal core [see Fig. 1(b)]. Similarly, one may also insert an
extra sector to the honeycomb lattice and adjust it uniformly in the circular
direction (by compressing the lattice) to obtain a disclination lattice with a
heptagonal core, as shown in Fig. 1(c).

To explore the spectrum of such finite disclination lattices, we neglect the
nonlinear term in Eq. (1) and search for the eigenmodes of the structure in
the form $\psi=u(x,y)e^{ibz}$, with $b$ being the propagation constant
(eigenvalue) and $u(x,y)$ being the modal field. In doing so, we obtain from
Eq. (1) the linear eigenvalue problem

$bu=\frac{1}{2}\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right)u+\mathcal{R}u,$ (3)

that can be solved numerically by using the plane-wave expansion method. In
Fig. 1(d) we present the spectrum of the modes of disclination lattice with
the pentagonal core as a function of parameter $r=\ell_{\rm intra}/\ell_{\rm
inter}$ determining the intra- and inter-cell separation between waveguides.
One finds that topological states appear in the forbidden gap of the spectrum
in the region $r>r_{c}=1$, with $r=r_{c}$ indicated in the figure by the
vertical dashed line. In Fig. 1(e) we plot the magnified part of the spectrum
that is indicated by the dashed rectangle in Fig. 1(d). One can observe the
formation of the disclination states shown by the red lines, and of the “zero-
energy” states (we keep this traditional notion for them that came from tight-
binding models also here despite the fact that their propagation constant can
actually change with $r$ in our continuous system) shown by the green lines.
One can see that disclination states in general are not degenerate, except for
the point $r\approx 1.79$. The spectrum of the disclination lattice with a
heptagonal core is presented in Fig. 1(f), while its magnified version is
shown in Fig. 1(g). The formation of disclination states in this lattice is
obvious too. Notice that localization of disclination states generally
increases with the increase of parameter $r$. The emergence of disclination
states is associated with transition of the lattice into topologically
nontrivial phase when parameter $r$ exceeds $r_{c}$, as described by the
associated topological invariant (see Methods).

Disclination lattices support several types of modes. To illustrate them we
choose a particular value of parameter $r=1.92$ exceeding the value $r_{c}=1$
(i.e. the lattice is in topological insulator phase) and show in Fig. 2(a) for
the lattice with a pentagonal core the eigenvalues of $n=450$ modes with
largest propagation constants. In this figure corner, edge, zero-energy, and
disclination states are highlighted by the colour dots, while bulk states are
indicated by the black dots. Field modulus distributions of some modes are
depicted in Fig. 2(b), where they are superimposed on the lattice profile with
indicated lattice cells (yellow lines) and sites (white circles). One finds
that both corner and zero-energy states form around the outer corners of the
finite lattice, i.e. they belong to the class of 0D states. Notice that in
contrast to corner states, zero-energy ones do not occupy the very corner site
of the lattice. In our sufficiently large lattice there are forty edge states
(cyan dots), ten corner states (yellow dots), and five degenerated zero-energy
states (blue dots) in the lattice with a pentagonal core. Beside these states
that can also be encountered in usual HOTIs [44, 50, 39, 38, 40, 41], there
are also five disclination states associated with the fractional disclination
charge (see Methods), which are labeled by the red dots with numbers $4\sim
8$. We display the first three disclination states (numbered $4\sim 6$) in the
second row of Fig. 2(b). In contrast to all other states, disclination states
localize at the five sites of the disclination core. As mentioned above, five
disclination states are not completely degenerate: only pair of states with
numbers 5 and 6 and pair of states 7 and 8 are degenerate ones. Degeneracy for
such states is lifted because they are localized on compact inner disclination
core, where sites located on the core do feel neighboring sites belonging to
the core and therefore interact. We also depict typical phase distributions
for disclination states in the insets of the second row of Fig. 2(b). The
disclination state numbered 4 with the largest propagation constant includes
five in-phase spots, while phase distributions in all other disclination
states are more complex, with alternating groups of in-phase and out-of-phase
spots.

Figure 4: Stability properties of the disclination and zero-energy states.
Maximal real part of the perturbation growth rate $\beta$ versus propagation
constant $b$ for the disclination state (a, b) and zero-energy state families
(c, d) in the lattice with a pentagonal core (a, c) and a heptagonal core (b,
d). Parameters are the same as in Fig. 3.

Disclination lattice with a heptagonal core also supports corner, edge, zero-
energy, and disclination states, as illustrated in Fig. 2(c) showing
eigenvalues of modes at $r=1.92$. Field modulus distributions in
representative states are illustrated in Fig. 2(d). According to the symmetry
of this disclination lattice, there are fifty-six edge states (cyan dots),
fourteen corner states (yellow dots), seven degenerate zero-energy states
(blue dots), and seven disclination states (red dots). Again, the disclination
states are also not completely degenerate. In this particular case, pair
states numbered 4 and 5, pair states numbered 6 and 7, and pair states
numbered 8 and 9 depicted in Fig. 2(d) are degenerate. Notice that the state
numbered 10 consists of seven in-phase spots. In the following section, we
will investigate nonlinear disclination states bifurcating from the linear
disclination state numbered 4 for the pentagonal case and state numbered 10
for the heptagonal case.

Figure 5: Dynamical excitation of the nonlinear disclination states. (a) The
inverse form-factor $f^{-1}$ of the field distribution versus propagation
distance for input Gaussian beams with different powers injected into single
site of the pentagonal disclination core at $r=2.17$. (b) Field modulus
distributions at different distances. Top and bottom panels correspond to the
black and red curves in (a), respectively. (c, d) Setup as in (a, b), but for
the structure with heptagonal core.

##  III Nonlinear disclination states

We now take into account focusing nonlinearity of the underlying material and
consider nonlinear disclination and zero-energy state families
$\psi=u(x,y)e^{ibz}$ bifurcating from corresponding linear states. Note that
similar phenomena can be expected if defocusing nonlinearity is considered.
Such families can be found using the Newton method and they are displayed in
Fig. 3. Such families are characterized by the dependence of the power
$P=\iint|\psi|^{2}dxdy$ of the nonlinear state on its propagation constant $b$
that is now free parameter defining the family [in Fig. 3 red lines show the
families of the disclination states, while blue lines show families of zero-
energy states]. The case of pentagonal core is illustrated in Fig. 3(a), while
the case of heptagonal core in Fig. 3(c). We found nonlinear states not only
in the forbidden gaps, but also in the bulk band shown as gray region (i.e.
nonlinear can drive propagation constant of nonlinear modes into bulk bands
and continuation of the same family can be found in higher gap, as illustrated
in the figure). One finds that the power $P$ of the nonlinear states increases
with increase of the propagation constant $b$, changing nearly linearly within
forbidden gaps. The shape of the topological disclination or zero-energy state
strongly depends on the location of the propagation constant in the spectrum:
in particular, the state gradually acquires long tails when the propagation
constant penetrates into the band, especially when the slope of $P(b)$
dependence changes and the power starts increasing rapidly. Transformation of
field modulus distributions in the nonlinear disclination (top row) and zero-
energy (bottom row) states is illustrated in Figs. 3(b) and 3(d) for
pentagonal and heptagonal lattices, respectively. We choose four nonlinear
states for each case, and two of them have propagation constants falling into
forbidden gap. The corresponding in-gap states are well-localized. Only the
in-band high-amplitude states corresponding to circles 4 couple with the bulk
modes thereby gradually extending across the entire lattice. Notice that both
disclination and zero-energy nonlinear families cross very narrow regions
associated with the bands of linear edge and corner states [left outermost
narrow gray region and vertical dashed lines in Figs. 3(a) and 3(c)]
practically without coupling with these states, which indicate that they
belong to the class of embedded solitons [95]. In the case of disclination
states this is explained by the fact that they form at the disclination core
located very far from the outer edges or corners of the structure, while in
the case of zero-energy states very weak coupling with corner state can be
explained by different symmetries of these states. Nevertheless, the weak
coupling can result in some cases in the instability of the nonlinear zero-
energy states, as will be shown later.

Among the most important properties of the nonlinear topological states is
their stability. We found that stability properties can be very different in
the lattices with pentagonal and heptagonal cores. We employ linear stability
analysis method to investigate stability of the nonlinear disclination states.
The method is based on substitution of perturbed solutions into governing Eq.
(1), linearization around stationary state and solution of the resulting
linear eigenvalue problem that yields perturbation growth rates $\beta$ (see
Methods). Nonlinear states are stable only if real part of the growth rate
$\beta_{\rm re}=\rm{Re}\\{\beta\\}$ is zero or negative for all possible
perturbations. We calculate perturbation growth rates for all nonlinear
disclination and zero-energy families from Figs. 3(a) and 3(c). Maximal real
part of the perturbation growth rate $\max\\{\beta_{\rm re}\\}$ versus $b$ is
plotted in Figs. 4(a, b) and 4(c, d) for the disclination and the zero-energy
state, respectively. In accordance with this analysis, stable branches were
plotted with solid lines in Figs. 3(a) and 3(c), while unstable ones were
plotted with dashed lines. One finds that the in-gap nonlinear disclination
states in the pentagonal disclination lattice are unstable in the largest part
of the gap, while those in the heptagonal disclination lattice are stable.
Interestingly, when nonlinear disclination states couple with bulk modes in
pentagonal disclination lattice, they become stable. In contrast, in
heptagonal disclination lattice such states destabilize when they enter the
bulk band. As shown in Figs. 4(c) and 4(d) the nonlinear zero-energy states
may be unstable in the vicinity of the edge/corner bands, but they become
completely stable with increase of the propagation constant $b$. When they
penetrate into the bulk band, nonlinear zero-energy states remain stable in
the pentagonal case, but become unstable in the heptagonal case.

Last but not least, we investigate the possibility of excitation of the
nonlinear disclination states by Gaussian input beams. To this end, we launch
a Gaussian beam into one of the sites of the disclination core. By adjusting
its power, we record the corresponding propagation dynamics. To monitor
transformation of the beam during propagation, we calculate the form-factor of
the field distribution at different distances that can be defined as:

$f^{2}=\left.{\iint|\psi|^{4}dxdy}\right/{\left(\iint|\psi|^{2}dxdy\right)^{2}}.$
(4)

The quantity $f^{-1}$ characterizes the width of the wavepacket. We first
investigate the case of the pentagonal disclination lattice. Corresponding
results illustrating excitation dynamics are shown in Figs. 5(a) and 5(b).
When the power of the Gaussian beam is low, $P=0.01$, the beam does not remain
exclusively in the excited site, the power couples between all five
disclination sites resulting in the oscillations of beam width upon
propagation, see black curve in Fig. 5(a) and field modulus distributions at
different distances shown in the first row of Fig. 5(b). In this case the
input Gaussian excites superposition of several disclination states and they
experience subsequent beatings, with excitation remaining always around
disclination core. Increasing the power of the Gaussian beam to $P=0.09$ leads
to strong localization in the only excited site, i.e. excitation of a strongly
localized variant of the nonlinear disclination state, when nonlinearity
suppresses coupling also between sites of the disclination core, see the red
curve in Fig. 5(a) and second row in Fig. 5(b). Similar results can be
obtained if we launch a Gaussian beam into the disclination core of the
heptagonal disclination lattice, as shown in Figs. 5(c) and 5(d). According to
the results of Fig. 5, dynamical excitation of the nonlinear disclination
states can be achieved at very low power levels (for example, nonlinear
disclination state in the heptagonal lattice is excited even if $P=0.0081$),
i.e. they are thresholdless, which is advantageous for potential experimental
implementation.

Figure 6: Wannier center (red dots) and spectral charge per unit cell in the
disclination lattice with a heptagonal core. (a) Topological insulator phase
with $r>r_{c}$. (b) Normal insulator phase with $r<r_{c}$. The unit cell with
same color has same spectral charge.

##  IV Conclusion

Summarizing, we have investigated nonlinear disclination states in
disclination lattices with a pentagonal/heptagonal cores that are constructed
by removing/inserting a sector with the Frank angle $2\pi/6$ from/into the
honeycomb lattice. The disclination lattice enters into topological insulator
phase upon variation of the separation among lattice sites in each unit cell,
this is accompanied by the shift of the Wannier centers to the boundaries of
the unit cell. In this phase, there appears fractional charge in the unit cell
around the disclination core due to filling anomaly, which predicts the
existence of disclination states according to the bulk-disclination
correspondence. We obtained rich nonlinear disclination state families
bifurcating from such linear topological states. We have confirmed that such
states are practically thresholdless in contrast to nonlinear self-sustained
states in periodic photonic lattices that always require certain minimal power
for their formation. The interesting property of the system considered here is
that disclination states coexist in it with zero-energy, corner and edge
states. Despite their coexistence, the bifurcation of nonlinear disclination
states into the gap occurs in such a way that they do not show any noticeable
coupling with other states (for example, with zero-energy or corner states)
present in the topological gap, even though their propagation constants may
overlap, as one can see in Fig. 3. Moreover, while the phenomenon of
bifurcation of localized gap states from defects is known even in trivial
periodic structures[96], we would like to stress here that in our case, in
contrast to trivial inclusion of local defect, the deformation of the lattice
after removal or addition of the sector is global and affects the entire
structure, in a sense that the structure becomes aperiodic, but still it can
possess the gap in the spectrum, where topological disclination states may
form. We have also found that stability of disclination states crucially
depends on the geometry of the lattice: nonlinear disclination states in the
lattice with pentagonal core are unstable in the largest part of the gap,
while nonlinear disclination states in the lattice with heptagonal core are
instead stable in the entire gap. Our results not only introduce the first
example of the nonlinear topological disclination states in photonic system,
but they also open a pathway to enhancement of light-matter interactions.
Additionally, we believe that the work will inspire investigations of
disclination states in non-Hermitian systems, realization of topological
lasing in such states, and design of compact on-chip photonic devices.

##  V Methods

###  V.1 Topological invariant

The topological invariant that can be associated with the band gap of the
disclination lattice [66, 67, 69, 70, 68] is defined by

$\chi=(\chi_{\rm M},\chi_{\rm K}),$ (5)

where the high symmetry indicators are $\chi_{\rm M}=\\#{\rm
M}_{1}^{(2)}-\\#\Gamma_{1}^{(2)}$ and $\chi_{\rm K}=\\#{\rm
K}_{1}^{(3)}-\\#\Gamma_{1}^{(3)}$. Here $\\#\Pi_{q}^{(n)}$ is the number of
bands below the bandgap at a high-symmetry point $\Pi=\Gamma,\,\rm M,\,K$ [see
the inset in Fig. 1(a)] with the eigenvalue of the $C_{n}$ rotation matrix
$e^{i2\pi(q-1)/n}~{}(q=1,\cdots,n)$. It can be shown for the disclination
lattices with pentagonal and heptagonal cores that $\chi=(2,\,0)$ for
$r>r_{c}$ and $\chi=(0,\,0)$ for $r<r_{c}$, that explains the appearance of
topologically nontrivial states. In addition, the bulk-disclination
correspondence establishes the link between the fractional disclination charge
and the localized states bound to the disclination core [70, 68]. The
fractional disclination charge $\mathcal{Q}$ is determined by the topological
index of the gap [66, 67, 69, 70, 68] and can be written as

$\mathcal{Q}=\frac{\Omega}{2\pi}\left(\frac{3}{2}\chi_{\rm M}-\chi_{\rm
K}\right)\,{\rm modulo}\,1,$ (6)

which means that the fractional charge is $\mathcal{Q}=1/2$ for $r>r_{c}$ and
$\mathcal{Q}=0$ for $r<r_{c}$. Here $\Omega=2\pi/6$ is the Frank angle. The
fractional disclination charge can be also obtained by counting the number of
the Wannier centers. For example, as shown in the insets in spectra presented
in Figs. 1(d) and 1(f), the Wannier centers are located at the edges of the
unit cell if $r>r_{c}$ and at its center if $r<r_{c}$. In Fig. 6, we show the
Wannier center and the spectral charge in each unit cell in the disclination
lattice with a heptagonal core. The unit cell around the disclination core has
five bulk Wannier centers if $r>r_{c}$, which give a $5/2$ spectral charge per
unit cell, as shown in Fig. 6(a). While if $r<r_{c}$, the spectral charge per
unit cell around the disclination core is $3$ [see Fig. 6(b)].

###  V.2 Linear stability analysis method

To perform linear stability analysis of the obtained nonlinear disclination
states, we write perturbed solution in the form

$\psi=\left[u(x,y)+v(x,y)e^{\beta z}+w^{*}(x,y)e^{\beta^{*}z}\right]e^{ibz},$
(7)

where $v,w\ll 1$ are small perturbations, and $\beta$ is the perturbation
growth rate that can be complex. Inserting Eq. (7) into Eq. (1) and
linearizing it around stationary solution $u$, we arrive at the linear
eigenvalue problem:

$\begin{split}i\beta v=\,&-\frac{1}{2}\left(\frac{\partial}{\partial
x^{2}}+\frac{\partial}{\partial
y^{2}}\right)v-({\mathcal{R}}-b)v-2|u|^{2}v-|u|^{2}w,\\\ i\beta
w=\,&+\frac{1}{2}\left(\frac{\partial}{\partial
x^{2}}+\frac{\partial}{\partial
y^{2}}\right)w+({\mathcal{R}}-b)w+2|u|^{2}w+|u|^{2}v.\\\ \end{split}$ (8)

Solving the problem (8) using standard eigenvalue solver, we obtain the
dependence of the perturbation growth rate $\beta$ (and associated
perturbation profiles $v,w$) on the propagation constant $b$ for a given
family of the nonlinear states. If $\beta_{\rm{re}}\leq 0$ for all possible
perturbations, corresponding nonlinear state $u$ is linearly stable and in the
presence of small perturbations it will exhibit only small amplitude
oscillations upon evolution, otherwise, if at least one perturbation mode has
${\beta_{\rm{re}}>0}$, the nonlinear state is unstable and will decay in the
course of propagation.

Acknowledgements
This work was supported by the National Natural Science Foundation of China
(12074308, U1537210), the Russian Science Foundation (21-12-00096), and the
Fundamental Research Funds for the Central Universities (xzy022022058).

Data availability
The data that supports the results within this paper are available from the
corresponding authors upon reasonable request.

Code availability
The codes used to obtain the results described in this paper are available
from the corresponding authors upon reasonable request.

Competing interests The authors declare no competing interests.

Correspondence and requests for materials should be addressed to Y.Z.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, “Colloquium: Topological insulators,” Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, 1057 (2011).
  * Haldane and Raghu [2008] F. D. M. Haldane and S. Raghu, “Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry,” Phys. Rev. Lett. 100, 013904 (2008).
  * Wang  _et al._ [2009] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, “Observation of unidirectional backscattering-immune topological electromagnetic states,” Nature 461, 772 (2009).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, “Photonic Floquet topological insulators,” Nature 496, 196 (2013).
  * Lindner  _et al._ [2011] N. H. Lindner, G. Refael, and V. Galitski, “Floquet topological insulator in semiconductor quantum wells,” Nat. Phys. 7, 490 (2011).
  * Hafezi  _et al._ [2011] M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, “Robust optical delay lines with topological protection,” Nat. Phys. 7, 907 (2011).
  * Stützer  _et al._ [2018] S. Stützer, Y. Plotnik, Y. Lumer, P. Titum, N. H. Lindner, M. Segev, M. C. Rechtsman, and A. Szameit, “Photonic topological Anderson insulators,” Nature 560, 461 (2018).
  * Yang  _et al._ [2019] Y. Yang, Z. Gao, H. Xue, L. Zhang, M. He, Z. Yang, R. Singh, Y. Chong, B. Zhang, and H. Chen, “Realization of a three-dimensional photonic topological insulator,” Nature 565, 622 (2019).
  * Maczewsky  _et al._ [2020] L. J. Maczewsky, M. Heinrich, M. Kremer, S. K. Ivanov, M. Ehrhardt, F. Martinez, Y. V. Kartashov, V. V. Konotop, L. Torner, D. Bauer, and A. Szameit, “Nonlinearity-induced photonic topological insulator,” Science 370, 701 (2020).
  * Biesenthal  _et al._ [2022] T. Biesenthal, L. J. Maczewsky, Z. Yang, M. Kremer, M. Segev, A. Szameit, and M. Heinrich, “Fractal photonic topological insulators,” Science 376, eabm2842 (2022).
  * Pyrialakos  _et al._ [2022] G. G. Pyrialakos, J. Beck, M. Heinrich, L. J. Maczewsky, N. V. Kantartzis, M. Khajavikhan, A. Szameit, and D. N. Christodoulides, “Bimorphic Floquet topological insulators,” Nat. Mater. 21, 634 (2022).
  * Barik  _et al._ [2018] S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi, and E. Waks, “A topological quantum optics interface,” Science 359, 666 (2018).
  * Mittal  _et al._ [2018] S. Mittal, E. A. Goldschmidt, and M. Hafezi, “A topological source of quantum light,” Nature 561, 502 (2018).
  * Yang  _et al._ [2015] Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, and B. Zhang, “Topological acoustics,” Phys. Rev. Lett. 114, 114301 (2015).
  * He  _et al._ [2016] C. He, X. Ni, H. Ge, X.-C. Sun, Y.-B. Chen, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, “Acoustic topological insulator and robust one-way sound transport,” Nat. Phys. 12, 1124 (2016).
  * Lu  _et al._ [2017] J. Lu, C. Qiu, L. Ye, X. Fan, M. Ke, F. Zhang, and Z. Liu, “Observation of topological valley transport of sound in sonic crystals,” Nat. Phys. 13, 369 (2017).
  * Ma  _et al._ [2019] G. Ma, M. Xiao, and C. T. Chan, “Topological phases in acoustic and mechanical systems,” Nat. Rev. Phys. 1, 281 (2019).
  * Jotzu  _et al._ [2014] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, “Experimental realisation of the topological Haldane model,” Nature 515, 237 (2014).
  * Huber [2016] S. D. Huber, “Topological mechanics,” Nat. Phys. 12, 621 (2016).
  * Nalitov  _et al._ [2015] A. V. Nalitov, D. D. Solnyshkov, and G. Malpuech, “Polariton $\mathbb{Z}$ topological insulator,” Phys. Rev. Lett. 114, 116401 (2015).
  * Klembt  _et al._ [2018] S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, R. Ge, M. A. Bandres, M. Emmerling, L. Worschech, T. C. H. Liew, M. Segev, C. Schneider, and S. Höfling, “Exciton-polariton topological insulator,” Nature 562, 552 (2018).
  * Liu  _et al._ [2020] W. Liu, Z. Ji, Y. Wang, G. Modi, M. Hwang, B. Zheng, V. J. Sorger, A. Pan, and R. Agarwal, “Generation of helical topological exciton-polaritons,” Science 370, 600 (2020).
  * Lu  _et al._ [2014] L. Lu, J. D. Joannopoulos, and M. Soljačić, “Topological photonics,” Nat. Photon. 8, 821 (2014).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, “Topological photonics,” Rev. Mod. Phys. 91, 015006 (2019).
  * St-Jean  _et al._ [2017] P. St-Jean, V. Goblot, E. Galopin, A. Lemaître, T. Ozawa, L. Le Gratiet, I. Sagnes, J. Bloch, and A. Amo, “Lasing in topological edge states of a one-dimensional lattice,” Nat. Photon. 11, 651 (2017).
  * Bandres  _et al._ [2018] M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, “Topological insulator laser: Experiments,” Science 359, eaar4005 (2018).
  * Bahari  _et al._ [2017] B. Bahari, A. Ndao, F. Vallini, A. El Amili, Y. Fainman, and B. Kanté, “Nonreciprocal lasing in topological cavities of arbitrary geometries,” Science 358, 636 (2017).
  * Parto  _et al._ [2018] M. Parto, S. Wittek, H. Hodaei, G. Harari, M. A. Bandres, J. Ren, M. C. Rechtsman, M. Segev, D. N. Christodoulides, and M. Khajavikhan, “Edge-mode lasing in 1D topological active arrays,” Phys. Rev. Lett. 120, 113901 (2018).
  * Zhao  _et al._ [2018] H. Zhao, P. Miao, M. H. Teimourpour, S. Malzard, R. El-Ganainy, H. Schomerus, and L. Feng, “Topological hybrid silicon microlasers,” Nat. Commun. 9, 981 (2018).
  * Kartashov and Skryabin [2019] Y. V. Kartashov and D. V. Skryabin, “Two-dimensional topological polariton laser,” Phys. Rev. Lett. 122, 083902 (2019).
  * Zeng  _et al._ [2020] Y. Zeng, U. Chattopadhyay, B. Zhu, B. Qiang, J. Li, Y. Jin, L. Li, A. G. Davies, E. H. Linfield, B. Zhang, Y. Chong, and Q. J. Wang, “Electrically pumped topological laser with valley edge modes,” Nature 578, 246 (2020).
  * Shalaev  _et al._ [2019] M. I. Shalaev, W. Walasik, A. Tsukernik, Y. Xu, and N. M. Litchinitser, “Robust topologically protected transport in photonic crystals at telecommunication wavelengths,” Nat. Nanotech. 14, 31 (2019).
  * Yang  _et al._ [2020] Y. Yang, Y. Yamagami, X. Yu, P. Pitchappa, J. Webber, B. Zhang, M. Fujita, T. Nagatsuma, and R. Singh, “Terahertz topological photonics for on-chip communication,” Nat. Photon. 14, 446 (2020).
  * Chen  _et al._ [2021] Y. Chen, X.-T. He, Y.-J. Cheng, H.-Y. Qiu, L.-T. Feng, M. Zhang, D.-X. Dai, G.-C. Guo, J.-W. Dong, and X.-F. Ren, “Topologically protected valley-dependent quantum photonic circuits,” Phys. Rev. Lett. 126, 230503 (2021).
  * Yang  _et al._ [2022] Q. Yang, D. Wang, S. Kruk, M. Liu, I. Kravchenko, J. Han, Y. Kivshar, and I. Shadrivov, “Topology-empowered membrane devices for terahertz photonics,” Adv. Photon. 4, 046002 (2022).
  * Peterson  _et al._ [2018] C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, “A quantized microwave quadrupole insulator with topologically protected corner states,” Nature 555, 346 (2018).
  * Xue  _et al._ [2019] H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, “Acoustic higher-order topological insulator on a kagome lattice,” Nat. Mater. 18, 108 (2019).
  * Ni  _et al._ [2019] X. Ni, M. Weiner, A. Alù, and A. B. Khanikaev, “Observation of higher-order topological acoustic states protected by generalized chiral symmetry,” Nat. Mater. 18, 113 (2019).
  * Chen  _et al._ [2019] X.-D. Chen, W.-M. Deng, F.-L. Shi, F.-L. Zhao, M. Chen, and J.-W. Dong, “Direct observation of corner states in second-order topological photonic crystal slabs,” Phys. Rev. Lett. 122, 233902 (2019).
  * Xie  _et al._ [2019] B.-Y. Xie, G.-X. Su, H.-F. Wang, H. Su, X.-P. Shen, P. Zhan, M.-H. Lu, Z.-L. Wang, and Y.-F. Chen, “Visualization of higher-order topological insulating phases in two-dimensional dielectric photonic crystals,” Phys. Rev. Lett. 122, 233903 (2019).
  * Mittal  _et al._ [2019] S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, “Photonic quadrupole topological phases,” Nat. Photon. 13, 692 (2019).
  * Serra-Garcia  _et al._ [2018] M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, “Observation of a phononic quadrupole topological insulator,” Nature 555, 342 (2018).
  * Noh  _et al._ [2018] J. Noh, W. A. Benalcazar, S. Huang, M. J. Collins, K. P. Chen, T. L. Hughes, and M. C. Rechtsman, “Topological protection of photonic mid-gap defect modes,” Nat. Photon. 12, 408 (2018).
  * Wang  _et al._ [2021a] H.-X. Wang, L. Liang, B. Jiang, J. Hu, X. Lu, and J.-H. Jiang, “Higher-order topological phases in tunable ${C_{3}}$ symmetric photonic crystals,” Photon. Res. 9, 1854 (2021a).
  * Imhof  _et al._ [2018] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, “Topolectrical-circuit realization of topological corner modes,” Nat. Phys. 14, 925 (2018).
  * Zhang  _et al._ [2019] X. Zhang, H.-X. Wang, Z.-K. Lin, Y. Tian, B. Xie, M.-H. Lu, Y.-F. Chen, and J.-H. Jiang, “Second-order topology and multidimensional topological transitions in sonic crystals,” Nat. Phys. 15, 582 (2019).
  * Zangeneh-Nejad and Fleury [2019] F. Zangeneh-Nejad and R. Fleury, “Nonlinear second-order topological insulators,” Phys. Rev. Lett. 123, 053902 (2019).
  * Qi  _et al._ [2020] Y. Qi, C. Qiu, M. Xiao, H. He, M. Ke, and Z. Liu, “Acoustic realization of quadrupole topological insulators,” Phys. Rev. Lett. 124, 206601 (2020).
  * El Hassan  _et al._ [2019] A. El Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, and M. Bourennane, “Corner states of light in photonic waveguides,” Nat. Photon. 13, 697 (2019).
  * Xu  _et al._ [2021] C. Xu, Z.-G. Chen, G. Zhang, G. Ma, and Y. Wu, “Multi-dimensional wave steering with higher-order topological phononic crystal,” Sci. Bull. 66, 1740 (2021).
  * Li  _et al._ [2022a] J. Li, Q. Mo, J.-H. Jiang, and Z. Yang, “Higher-order topological phase in an acoustic fractal lattice,” Sci. Bull. 67, 2040 (2022a).
  * Zheng  _et al._ [2022] S. Zheng, X. Man, Z.-L. Kong, Z.-K. Lin, G. Duan, N. Chen, D. Yu, J.-H. Jiang, and B. Xia, “Observation of fractal higher-order topological states in acoustic metamaterials,” Sci. Bull. 67, 2069 (2022).
  * Xie  _et al._ [2021] B. Xie, H.-X. Wang, X. Zhang, P. Zhan, J.-H. Jiang, M. Lu, and Y. Chen, “Higher-order band topology,” Nat. Rev. Phys. 3, 520 (2021).
  * Peterson  _et al._ [2020] C. W. Peterson, T. Li, W. A. Benalcazar, T. L. Hughes, and G. Bahl, “A fractional corner anomaly reveals higher-order topology,” Science 368, 1114 (2020).
  * Kruk  _et al._ [2021] S. S. Kruk, W. Gao, D.-Y. Choi, T. Zentgraf, S. Zhang, and Y. Kivshar, “Nonlinear imaging of nanoscale topological corner states,” Nano Lett. 21, 4592 (2021).
  * Zhang  _et al._ [2020] W. Zhang, X. Xie, H. Hao, J. Dang, S. Xiao, S. Shi, H. Ni, Z. Niu, C. Wang, K. Jin, X. Zhang, and X. Xu, “Low-threshold topological nanolasers based on the second-order corner state,” Light Sci. Appl. 9, 109 (2020).
  * Kim  _et al._ [2020] H.-R. Kim, M.-S. Hwang, D. Smirnova, K.-Y. Jeong, Y. Kivshar, and H.-G. Park, “Multipolar lasing modes from topological corner states,” Nat. Commun. 11, 5758 (2020).
  * Han  _et al._ [2020] C. Han, M. Kang, and H. Jeon, “Lasing at multidimensional topological states in a two-dimensional photonic crystal structure,” ACS Photon. 7, 2027 (2020).
  * Zhong  _et al._ [2021a] H. Zhong, Y. V. Kartashov, A. Szameit, Y. Li, C. Liu, and Y. Zhang, “Theory of topological corner state laser in Kagome waveguide arrays,” APL Photon. 6, 040802 (2021a).
  * Ota  _et al._ [2019] Y. Ota, F. Liu, R. Katsumi, K. Watanabe, K. Wakabayashi, Y. Arakawa, and S. Iwamoto, “Photonic crystal nanocavity based on a topological corner state,” Optica 6, 786 (2019).
  * Xie  _et al._ [2020] X. Xie, W. Zhang, X. He, S. Wu, J. Dang, K. Peng, F. Song, L. Yang, H. Ni, Z. Niu, C. Wang, K. Jin, X. Zhang, and X. Xu, “Cavity quantum electrodynamics with second-order topological corner state,” Laser Photon. Rev. 14, 1900425 (2020).
  * Rüegg and Lin [2013] A. Rüegg and C. Lin, “Bound states of conical singularities in graphene-based topological insulators,” Phys. Rev. Lett. 110, 046401 (2013).
  * Teo and Hughes [2013] J. C. Y. Teo and T. L. Hughes, “Existence of Majorana-Fermion bound states on disclinations and the classification of topological crystalline superconductors in two dimensions,” Phys. Rev. Lett. 111, 047006 (2013).
  * Benalcazar  _et al._ [2014] W. A. Benalcazar, J. C. Y. Teo, and T. L. Hughes, “Classification of two-dimensional topological crystalline superconductors and Majorana bound states at disclinations,” Phys. Rev. B 89, 224503 (2014).
  * Benalcazar  _et al._ [2019] W. A. Benalcazar, T. Li, and T. L. Hughes, “Quantization of fractional corner charge in ${C}_{n}$-symmetric higher-order topological crystalline insulators,” Phys. Rev. B 99, 245151 (2019).
  * Li  _et al._ [2020] T. Li, P. Zhu, W. A. Benalcazar, and T. L. Hughes, “Fractional disclination charge in two-dimensional ${C}_{n}$-symmetric topological crystalline insulators,” Phys. Rev. B 101, 115115 (2020).
  * Wu  _et al._ [2021] S. Wu, B. Jiang, Y. Liu, and J.-H. Jiang, “All-dielectric photonic crystal with unconventional higher-order topology,” Photon. Res. 9, 668 (2021).
  * Peterson  _et al._ [2021] C. W. Peterson, T. Li, W. Jiang, T. L. Hughes, and G. Bahl, “Trapped fractional charges at bulk defects in topological insulators,” Nature 589, 376 (2021).
  * Liu  _et al._ [2021] Y. Liu, S. Leung, F.-F. Li, Z.-K. Lin, X. Tao, Y. Poo, and J.-H. Jiang, “Bulk-disclination correspondence in topological crystalline insulators,” Nature 589, 381 (2021).
  * Chen  _et al._ [2022] Y. Chen, Y. Yin, Z.-K. Lin, Z.-H. Zheng, Y. Liu, J. Li, J.-H. Jiang, and H. Chen, “Observation of topological $p$-orbital disclination states in non-Euclidean acoustic metamaterials,” Phys. Rev. Lett. 129, 154301 (2022).
  * Wang  _et al._ [2021b] Q. Wang, Y. Ge, H.-x. Sun, H. Xue, D. Jia, Y.-j. Guan, S.-q. Yuan, B. Zhang, and Y. D. Chong, “Vortex states in an acoustic Weyl crystal with a topological lattice defect,” Nat. Commun. 12, 3654 (2021b).
  * Wang  _et al._ [2020] Q. Wang, H. Xue, B. Zhang, and Y. D. Chong, “Observation of protected photonic edge states induced by real-space topological lattice defects,” Phys. Rev. Lett. 124, 243602 (2020).
  * Deng  _et al._ [2022] Y. Deng, W. A. Benalcazar, Z.-G. Chen, M. Oudich, G. Ma, and Y. Jing, “Observation of degenerate zero-energy topological states at disclinations in an acoustic lattice,” Phys. Rev. Lett. 128, 174301 (2022).
  * Smirnova  _et al._ [2020] D. Smirnova, D. Leykam, Y. Chong, and Y. Kivshar, “Nonlinear topological photonics,” Appl. Phys. Rev. 7, 021306 (2020).
  * Kartashov and Skryabin [2016] Y. V. Kartashov and D. V. Skryabin, “Modulational instability and solitary waves in polariton topological insulators,” Optica 3, 1228 (2016).
  * Kartashov and Skryabin [2017] Y. V. Kartashov and D. V. Skryabin, “Bistable topological insulator with exciton-polaritons,” Phys. Rev. Lett. 119, 253904 (2017).
  * Hadad  _et al._ [2018] Y. Hadad, J. C. Soric, A. B. Khanikaev, and A. Alù, “Self-induced topological protection in nonlinear circuit arrays,” Nat. Electron. 1, 178 (2018).
  * Lumer  _et al._ [2013] Y. Lumer, Y. Plotnik, M. C. Rechtsman, and M. Segev, “Self-localized states in photonic topological insulators,” Phys. Rev. Lett. 111, 243905 (2013).
  * Mukherjee and Rechtsman [2020] S. Mukherjee and M. C. Rechtsman, “Observation of Floquet solitons in a topological bandgap,” Science 368, 856 (2020).
  * Leykam and Chong [2016] D. Leykam and Y. D. Chong, “Edge solitons in nonlinear-photonic topological insulators,” Phys. Rev. Lett. 117, 143901 (2016).
  * Ablowitz and Cole [2017] M. J. Ablowitz and J. T. Cole, “Tight-binding methods for general longitudinally driven photonic lattices: Edge states and solitons,” Phys. Rev. A 96, 043868 (2017).
  * Ivanov  _et al._ [2020] S. K. Ivanov, Y. V. Kartashov, A. Szameit, L. Torner, and V. V. Konotop, “Vector topological edge solitons in Floquet insulators,” ACS Photon. 7, 735 (2020).
  * Zhong  _et al._ [2021b] H. Zhong, S. Xia, Y. Zhang, Y. Li, D. Song, C. Liu, and Z. Chen, “Nonlinear topological valley Hall edge states arising from type-II Dirac cones,” Adv. Photon. 3, 056001 (2021b).
  * Jürgensen  _et al._ [2021] M. Jürgensen, S. Mukherjee, and M. C. Rechtsman, “Quantized nonlinear Thouless pumping,” Nature 596, 63 (2021).
  * Fu  _et al._ [2022] Q. Fu, P. Wang, Y. V. Kartashov, V. V. Konotop, and F. Ye, “Nonlinear Thouless pumping: Solitons and transport breakdown,” Phys. Rev. Lett. 128, 154101 (2022).
  * Kirsch  _et al._ [2021] M. S. Kirsch, Y. Zhang, M. Kremer, L. J. Maczewsky, S. K. Ivanov, Y. V. Kartashov, L. Torner, D. Bauer, A. Szameit, and M. Heinrich, “Nonlinear second-order photonic topological insulators,” Nat. Phys. 17, 995 (2021).
  * Kartashov  _et al._ [2022] Y. V. Kartashov, A. A. Arkhipova, S. A. Zhuravitskii, N. N. Skryabin, I. V. Dyakonov, A. A. Kalinkin, S. P. Kulik, V. O. Kompanets, S. V. Chekalin, L. Torner, and V. N. Zadkov, “Observation of edge solitons in topological trimer arrays,” Phys. Rev. Lett. 128, 093901 (2022).
  * Tan  _et al._ [2021] D. Tan, Z. Wang, B. Xu, and J. Qiu, “Photonic circuits written by femtosecond laser in glass: improved fabrication and recent progress in photonic devices,” Adv. Photon. 3, 024002 (2021).
  * Li  _et al._ [2022b] L. Li, W. Kong, and F. Chen, “Femtosecond laser-inscribed optical waveguides in dielectric crystals: a concise review and recent advances,” Adv. Photon. 4, 024002 (2022b).
  * Fu  _et al._ [2020] Q. Fu, P. Wang, C. Huang, Y. V. Kartashov, L. Torner, V. V. Konotop, and F. Ye, “Optical soliton formation controlled by angle twisting in photonic moiré lattices,” Nat. Photon. 14, 663 (2020).
  * Hu  _et al._ [2021] Z. Hu, D. Bongiovanni, D. Jukić, E. Jajtić, S. Xia, D. Song, J. Xu, R. Morandotti, H. Buljan, and Z. Chen, “Nonlinear control of photonic higher-order topological bound states in the continuum,” Light Sci. Appl. 10, 164 (2021).
  * Wu and Hu [2015] L.-H. Wu and X. Hu, “Scheme for achieving a topological photonic crystal by using dielectric material,” Phys. Rev. Lett. 114, 223901 (2015).
  * Yves  _et al._ [2017] S. Yves, R. Fleury, T. Berthelot, M. Fink, F. Lemoult, and G. Lerosey, “Crystalline metamaterials for topological properties at subwavelength scales,” Nat. Commun. 8, 16023 (2017).
  * Champneys  _et al._ [2001] A. Champneys, B. Malomed, J. Yang, and D. Kaup, “Embedded solitons: solitary waves in resonance with the linear spectrum,” Phys. D 152–153, 340 (2001).
  * Sukhorukov  _et al._ [2001] A. A. Sukhorukov, Y. S. Kivshar, O. Bang, J. J. Rasmussen, and P. L. Christiansen, “Nonlinearity and disorder: Classification and stability of nonlinear impurity modes,” Phys. Rev. E 63, 036601 (2001).


# Modeling Turbulent Flows with LSTM Neural Network

Hugo D. Pasinato<EMAIL_ADDRESS>Nicolás F. Moguilner Rhe

In this study, we explore the application of an artificial recurrent neural
network (RNN) called Long Short-Term Memory (LSTM) as an alternative to a
turbulent Reynolds-Averaged Navier-Stokes (RANS) model. The LSTM models are
utilized to predict the shear Reynolds stress in developed and developing
turbulent channel flows. We conduct comparative analyses, comparing the LSTM
results propagated through computational fluid dynamics (CFD) simulations with
the outcomes from the $\kappa-\epsilon$ model and data acquired from direct
numerical simulation (DNS). These analyses demonstrate a good performance of
the LSTM approach.

##  1 Introduction

Fluid flow predictions are of utmost importance in various fields of
engineering, including mechanical, chemical, and aeronautical, to achieve
diverse objectives such as design optimization and equipment enhancement. Many
of these engineering applications involve turbulent flows, which require the
use of Reynolds-averaged Navier-Stokes (RANS) modeling. RANS is a commonly
employed technique for high-velocity scenarios with high Reynolds numbers. It
involves utilizing the Reynolds averaging method to derive the time-averaged
general equations of fluid motion, known as the RANS equations (Pope (2000)).
These equations, expressed in dimensionless form using tensor notation, are
the conservation of mass

$\frac{\partial U_{i}}{\partial x_{i}}=0,$ (1)

and the momentum equations

$\frac{\partial U_{i}}{\partial t}+\frac{\partial(U_{j}U_{i})}{\partial
x_{j}}\;=\;-\;\frac{\partial P}{\partial
x_{i}}\;+\;\;\frac{1}{Re_{\tau}}\frac{\partial^{2}U_{i}}{\partial
x_{j}\partial x_{j}}-\frac{\partial\langle u_{j}u_{i}\rangle}{\partial
x_{j}},$ (2)

where $U_{i}$ (equivalent to U, V, and W), is the mean velocity and $P$ the
mean pressure, $x$ and $t$ are the space and time coordinates, respectively,
and $R_{\tau}$ is the friction Reynolds number $u_{\tau}\delta/\nu$. All
variables are non-dimensionalized using the friction velocity, $u_{\tau}$ and
half the distance between walls, $\delta$, which represent the characteristic
velocity and length scales. Furthermore, when the original Navier-Stokes
equations are time-averaged, a tensor arises consisting of new unknowns, known
as the Reynolds stress tensor, $\langle u_{j}u_{i}\rangle$. Given that this
tensor is symmetric, for solving these equations there are 6 new unknowns, in
addition to $U_{i}$ and P.

Previous task involves incorporating a procedure or method for computing these
6 new unknowns (the velocities U, V, and W and pressure P are found using the
conservation equations of momentum and mass). These techniques are called RANS
modeling for turbulence, as commented above, with the Kappa-Epsilon model
being an example of them. In general, these models employ a pseudo viscosity
that is determined based on the influence of turbulence. This approach draws
an analogy between the behavior of a Newtonian fluid and that of a turbulent
flow. By incorporating this pseudo viscosity, the models aim to capture the
characteristics and dynamics of turbulent flows in a manner similar to how a
Newtonian fluid behaves.

RANS models depend on empirically-determined constants that are calibrated
using experimental turbulent flows or more accurate techniques such as Large
Eddy Simulation (LES) or Direct Numerical Simulation (DNS). However, when
these models are utilized beyond the range of their calibrated constants,
their predictions may lack accuracy. This represents the primary challenge
associated with RANS techniques.

Following a prolonged period of stagnation, RANS modeling is now entering a
new phase with the integration of Artificial Intelligence (AI) techniques. AI
involves a diverse range of data-driven algorithms, ranging from well-known
methods like linear regression to more advanced concepts like neural networks.
With the availability of high-fidelity data, AI allows for the capturing of
potentially complex relationships between turbulence mean-flow characteristics
(features) and modeling terms (predictions)

In the last years AI has been widely used in RANS models, employing a variety
of methods for addressing a wide range of problems (Milano and Koumoutsakos
(2002); Lecun et al. (2015); Ling and Templeton (2015); Ling et al. (2016a)).
One such method that has gained attention for its flexibility and precision is
deep learning (DL), which involves transforming input features through
multiple layers of nonlinear interaction (Ling et al. (2016b)).

Turbulence, however, is a highly complex phenomenon that often proves
challenging to accurately model using conventional AI techniques like deep
learning (DL). One of the primary reasons for this complexity lies in the non-
local effects of turbulence. The Reynolds stresses or turbulence stresses not
only depend on the flow characteristics at the specific location where they
are calculated but also on the flow characteristics in the surrounding
regions. This interdependency poses a significant bottleneck in turbulence
modeling. To address these non-linear problems more effectively, recurrent
neural networks (RNN) are better suited due to their ability to capture
spatial or temporal dependencies.

In particular, Long Short-Term Memory (LSTM) recurrent neural network is
currently regarded as one of the most interesting types of neural networks,
with potential for effectively capturing the non-locality of turbulence. LSTMs
have found widespread applications in various fields, including language
modeling, sentiment analysis, machine translation, speech recognition, and
time series forecasting. In engineering, LSTMs are particularly useful for
predicting trends and patterns in time series data, such as stock prices,
weather patterns, electricity demand, anomaly detection in sensor data,
industrial machinery health monitoring, and equipment failure prediction,
among others (Hochriter and Schmidhuber (1997); Sutskever et al. (2014)).
LSTMs excel in tasks where data has a sequential nature and long-term
dependencies need to be captured, making them a promising approach for
turbulence modeling.

Figure 1: Mean flow and Reynolds stress dimensionless longitudinal profiles
from DNS data for a channel flow with $Re_{\tau}=300$, perturbed with an
adverse pressure gradient at a narrow region of wide $W^{+}=220$ at the buffer
region. Vertical solid lines are the perturbation region limits, and the 3
vertical broken-lines, show middle, $W^{+}$ and $5W^{+}$ distance from
perturbation region. Solid line, $\partial U^{+}/\partial y^{+}$;
$-.-.-$$\partial P^{+}/\partial y^{+}$; $---$, $20\times\partial
U^{+}/\partial x^{+}$; $\bullet\cdot\bullet\cdot\bullet$, $20\times\partial
V^{+}/\partial x^{+}$; $\Box\cdot\Box\cdot\Box$ , $50\times\langle
uv\rangle^{+}$; $\diamond\cdot\diamond\cdot\diamond$$10\times\langle
uu\rangle^{+}$ .

In this study, an application of LSTM is performed to predict Reynolds
stresses from data generated by Direct Numerical Simulation (DNS). The
following section discusses aspects related to the Reynolds stress prediction
using RNN; section 3 presents and discusses some of the LSTM models tested in
the study; section 4 gives some numerical details of the RANS simulations;
section 5 discusses and shows the results of the neural network fitting and
propagation of predictions through RANS simulations, and finally section 6
lists some conclusions.

##  2 Reynolds stress prediction with LSTM RNN

Name | $R_{\tau}$ | Perturbation | Parameter | Purpose
---|---|---|---|---
Developed | 300 |  |  | train. test.
BRE3006 | 300 | blowing | $v^{+}=0.6$ | train. test.
SRE3006 | 300 | suction | $v^{+}=0.6$ | train. test.
APGRE30025 | 300 | adv. p. grad. step | $\partial P^{+}/\partial x^{+}=+0.25$ | train. test.
FPGRE30025 | 300 | fav. p. grad. step | $\partial P^{+}/\partial x^{+}=-0.25$ | train. test.
BRE3004 | 300 | blowing | $v^{+}=0.4$ | val.
SRE3004 | 300 | suction | $v^{+}=0.4$ | val.
APG100045 | 600 | adv. p. grad. step | $\partial P^{+}/\partial x^{+}=+0.40$ | val.
FPG100045 | 600 | fab. p. grad. step | $\partial P^{+}/\partial x^{+}=-0.40$ | val.
Developed | 1000 |  |  | val.
Table 1: Data set for training, testing or validation. $v^{+}$ is the
dimensionless wall-normal velocity in the slot at the wall; $\partial
P^{+}/\partial x^{+}$ is the pressure gradient step in the perturbation
region.

One of the most challenging aspects of turbulence is its non-local nature.
Let’s consider a scenario where we have a statistically stationary turbulent
flow (the mean flow characteristics remain constant with time). When this flow
is disturbed by an immersed body, for example, the primary effect is on the
pressure field, mainly upstream and around the body. These pressure field
changes subsequently influence the mean velocity field around the object,
resulting in an increase in turbulence primarily downstream of the body. In
other words, the modifications in mean pressure and velocity, and turbulence
fields, occur in different regions of the flow, emphasizing the non-local
characteristics of turbulent flows.

Figure 1 shows the longitudinal profiles of mean flow and Reynolds stress in a
channel flow perturbed by an adverse pressure gradient in a narrow region near
the wall. This figure reveals that the maximum values of the shear Reynolds
stress $\langle uv\rangle$ occurs downstream of the perturbation region, which
is clearly downstream where the maximum changes in the mean pressure and
velocity are observed.

In the context of RANS turbulence models, such as the $\kappa-\epsilon$ model
(Pope (2000)), the inclusion of non-local effects is achieved through the
convection and diffusion of kinetic energy ($\kappa$) and turbulence
dissipation ($\epsilon$). However, incorporating these non-local effects
becomes more challenging when using a regression model in deep learning (DL).
In this case, it is more appropriate to utilize a recurrent neural network
such as LSTM.

Another aspect of an NN model for Reynolds stress is its universality. The
question arises regarding the level of universality and complexity the model
should possess, or conversely, how specialized and simple it can be.
Developing a universal model with an NN would necessitate a large number of
adjustable parameters to capture all possible flow characteristics.
Furthermore, fine-tuning an NN with a high parameter count is exceedingly
challenging, and integrating a neural network with a large number of
parameters into a computational fluid dynamics (CFD) software presents a
significant obstacle due to the substantial computational demands it entails.

Consequently, this study aims to utilize an RNN with a reduced parameter count
for rapid prediction in a RANS simulation. To achieve this, the LSTM neural
network was trained, tested, and validated, using data generated from Direct
Numerical Simulations (DNS) of fully developed or perturbed (developing)
turbulent flows. All the flows considered in this study are statistically 2D
channel flows, comprising the following scenarios: a) developed flows; b)
flows perturbed with wall injection; c) flows perturbed with wall suction; d)
flows perturbed with an adverse pressure gradient in a localized region; and
e) flows perturbed with a favorable pressure gradient in a localized region.
Each perturbed flow initially represents a developed flow, with a small region
of width W+=220 being perturbed. To generate boundary conditions for a
perturbed flow, two parallel DNS simulations were conducted, with the first
simulation employing periodic boundary conditions. For more detailed
information, please refer to the bibliography for additional insights
(Pasinato (2012)). Table 1 provides a comprehensive list of the flows utilized
in the study, accompanied by pertinent details such as perturbation parameters
and their respective purposes.

Given that the turbulent flows studied are statistically stationary, the
neural network is trained using space-sequences instead of time-sequences.
Each space-sequence consists of 64 consecutive values along the longitudinal
direction, with an approximate dimensionless distance of $X^{+}\simeq 41$
between each value. Consequently, the complete sequence covers a distance of
approximately $X^{+}\simeq 2620$, representing the entire physical domain. In
other words, if for a DNS with $R_{\tau}=300$ a numerical grid with $N=254$
nods is used in the longitudinal direction, for each spatial sequence N/4
nodes of the DNS solution are employed at a fixed $y^{+}$ distance from the
wall.

To ensure a diverse representation during training, the full sequences are
randomly shuffled. Additionally, the values are normalized between 0 and 1,
using the global maximum and minimum values of the data set. This
normalization process enables consistent scaling across different variables.
Finally, the entire data set is exclusively divided into training and testing
data, with the testing data comprising $20\%$ of the complete data set.

##  3 LSTM models

Name | Parameters | Learning rate | Features | Predictions
---|---|---|---|---
L110DL | 531 | 0.025-0.001 |  $S_{11}^{+}$; $S_{12}^{+}$ | $\langle uv\rangle^{+}$
L15L25DL | 386 | 0.025-0.001 |  $S_{11}^{+}$; $S_{12}^{+}$ | $\langle uv\rangle^{+}$
L15DL-Y | 186 | 0.025-0.001 |  $S_{11}^{+}$; $S_{12}^{+};Y$ | $\langle uv\rangle^{+}$
L110DL-Y | 571 | 0.025-0.001 |  $S_{11}^{+}$; $S_{12}^{+};Y$ | $\langle uv\rangle^{+}$
Table 2: LSTM tested models.

Figure 2: Prediction of $\langle uv\rangle^{+}$ with the LSTM model L110DL-Y.
Solid line, $S_{12}^{+}$; $---$, $S_{11}^{+}$; $-.-.-$, $y/\delta$;
$\Box\cdot\Box\cdot\Box$ , Prediction; $\circ\cdot\circ\cdot\circ\;$, DNS
data. The data set was normalized in the 0-1 range with the global maximum and
minimum.

Figure 3: Mean velocity for a developed channel flow with $Re_{\tau}=300$ with
$\langle uv\rangle^{+}$ predicted from the LSTM model L110DL-Y. Solid line,
$U^{+}$; $+\cdot+\cdot+$, Logarithmic profile $1/0.41\;ln(y^{+})+5.5$.

Figure 4: Shear stresses for a developed channel flow with $Re_{\tau}=300$
with $\langle uv\rangle^{+}$ predicted with the LSTM model L110DL-Y. Filled
diamonds, total stresses; $\Box\cdot\Box\cdot\Box$ , $\langle uv\rangle^{+}$;
$\circ\cdot\circ\cdot\circ\;$, molecular stress $dU^{+}/dy^{+}$.

All LSTM models were implemented in two different numerical codes: one
developed in-house using modern Fortran, and the other implemented in Python
using the Keras-TensorFlow libraries (Abadi et al. (2015)). The porpouse of
the Fortran code was to gain insight into the forward and backward algorithms
of LSTM and compare them with the Keras version ( Chen (2018)).

The LSTM, like other deep learning models, has three hyper-parameters: a) the
number of stacked layers (or stacked LSTM cells), b) the number of memory
units in each LSTM layer, and c) the learning rate. In this study, these
hyper-parameters were selected through trial and error.

It is worth mentioning that the TensorFlow libraries offer a wide range of
optimization methods, including dynamic learning-rate techniques that speed up
convergence. However, both the Fortran and Keras codes faced challenges when
optimizing models with many stacked layers (5 or more layers) and numerous
memory units. To simplify the models and reduce the number of parameters, a
maximum of 3 stacked LSTM layers with a maximum of 10 memory units were used
for all tested models. Finally the best predicted solution were from models
with 1 LSTM cell and 10 memory units for the data set used.

Table 2 presents some of the LSTM models tested in this study. The name
$L15L25DL$ indicates an LSTM with 2 stacked layers, each having 5 memory
units, and a dense output layer. The components $S_{11}^{+}$ and $S_{12}^{+}$
represent the dimensionless mean rate-of-strain-tensor $S_{ij}$, which are
calculated as $S_{ij}^{+}=1/2(\partial U^{+}_{i}/\partial x^{+}_{j}+\partial
U^{+}_{j}/\partial x^{+}_{i})$. Additionally, $Y$ denotes the dimensionless
wall distance ($y/\delta$). It is important to note that the axes in $S{ij}$
are non-dimensionalized differently from the distance to the wall $Y$.

Although the primary objective of this study was to predict the shear Reynolds
stress $\langle uv\rangle$, some of the models were developed to predict the
normal Reynolds stress $\langle uu\rangle$, as an initial step, followed by an
LSTM model to predict $\langle uv\rangle$ based on the predicted $\langle
uu\rangle$. In all of these cases, the optimization of the LSTM models for
$\langle uu\rangle$ yielded superior and faster results compared to the LSTM
models for $\langle uv\rangle$.

The previous architecture of one RNN to predict $\langle uu\rangle$ and a
second RNN to predict $\langle uv\rangle$ is based on the non-linearity of
turbulence. In the context of channel turbulent flows, it is well-established
that turbulence acquires energy from the mean flow first through $\langle
uu\rangle$. Subsequently, the instantaneous pressure field transfers a portion
of this energy to the $\langle vv\rangle$ and $\langle ww\rangle$ components
(Tennekes and Lumley (1972)). As consequence as depicted in Fig. 1, the
modification of $\langle uu\rangle$ occurs promptly after the mean velocity
and pressure field are perturbed, while the alteration of $\langle uv\rangle$
experiences a certain delay. No results are presented here of the LSTM model
of $\langle uv\rangle$ based on a previous prediction of $\langle uu\rangle$.

##  4 Numerical details of the RANS simulations

The RANS simulations with the propagation of the LSTM predictions and the
$\kappa-\epsilon$ model were computed with a physical domain of
$3\pi\delta\times 2\delta\times 4/3\pi\delta$, and a numerical grid of
$64\times 64\times 64$ ($\Delta x^{+}=44$, $\Delta z^{+}=19.6$, $\Delta
y^{+}_{max}=26.4$, and $\Delta y^{+}_{min}=1$.), and the time step
$0.15\nu/u_{\tau}^{2}$. In the wall-normal direction a non-uniform mesh
distributed with the hyperbolic tangent function is used and the expansion
ratio is adjusted to ensure that the $y^{+}$ of the first cell center is equal
to 1. Furthermore, a van Driest function near the wall is used with the
$\kappa-\epsilon$ model. Considering the low Reynolds number flows employed in
this study, it is believed that the wall modeling approach adopted is
adequate.

The RANS simulation of perturbed flow follows a similar approach to DNS.
Initially, a fully developed flow is considered, and a small region of width
W+=220 is perturbed. To generate boundary conditions for the perturbed flow, a
second simulation is conducted in parallel with periodic boundary conditions,
representing the fully developed flow. In other words, at every time-step, the
input boundary condition for the perturbed flow is extracted from the middle
of the computational domain of the fully developed flow.

##  5 Results and Discussion

Figure 5: Comparison of half of the skin-friction for a channel flow with
$Re_{\tau}=300$ perturbed with blowing with an injection dimensionless
velocity $0.40$ from a narrow slot at the wall, with $\langle uv\rangle^{+}$
predicted with the LSTM L110DL-Y model. Solid line, LSTM; $---$,
$\kappa-\epsilon$ model; $\Box\cdot\Box\cdot\Box$ , DNS.

In this section, we present the results of the LSTM predictions and their
integration with a CFD code for propagation. The term ”propagation” refers to
the utilization of the LSTM RNN as a substitute for a RANS turbulence model.
At each time-step, sequences are formed and the LSTM model is employed to
calculate the Reynolds stress $\langle uv\rangle$ at every node of the grid.
While we provide some comparisons between the LSTM prediction and results
obtained from the $\kappa-\epsilon$ model (which computes all shear and normal
Reynolds stresses), the RANS simulations with LSTM predictions only consider
the shear Reynolds stress $\langle uv\rangle$. The remaining Reynolds stresses
were taken equal to zero.

As mentioned earlier, the purpose of using the LSTM RNN in this study is not
to capture the universal behavior of turbulence, but rather to predict
statistical stationary 2D turbulent flows specifically for channel flow. The
LSTM model serves as a proof of concept for predicting Reynolds stresses.
Consequently, the turbulent flows utilized for validating the LSTM model share
certain similarities with the flows employed for model adjustment and testing.

For instance, all the turbulent flows used for model adjustment, testing, and
validation pertain to channel flows. These flows encompass developed or
perturbed developing flows, characterized by low or medium friction Reynolds
numbers $R_{\tau}$ (300-1000). It is worth noting that the mean velocity and
normal Reynolds stresses, along with other statistical measures, exhibit low
Reynolds number effects in these flows. The developing flows represent
turbulent flows that have been perturbed from the wall or in small regions in
the buffer region.

Figure 6: Comparison of half of the skin-friction for a channel flow with
$Re_{\tau}=300$ perturbed with suction from a narrow slot at the wall, with
$\langle uv\rangle^{+}$ predicted with the LSTM L110DL-Y model. Solid line,
LSTM; $---$, $\kappa-\epsilon$ model; $\Box\cdot\Box\cdot\Box$ , DNS.

Figure 7: Comparison of the shear Reynolds stress for a channel flow with
$Re_{\tau}=300$ perturbed with suction from a narrow slot at the wall, with
$\langle uv\rangle^{+}$ predicted from the L110DL-Y model. $y^{+}=5$, solid
line, LSTM; $\circ\cdot\circ\cdot\circ\;$DNS; $y^{+}=30$, $---$, LSTM;
$\diamond\cdot\diamond\cdot\diamond$, DNS.

Prediction of the shear Reynolds stress without propagation: Figure 2
illustrates the LSTM-based prediction of the shear Reynolds stress $\langle
uv\rangle^{+}$ utilizing $S_{11}^{+}$, $S_{12}^{+}$, and Y($=y/\delta$) as
input parameters, with L110DL-Y model. The LSTM model employed for this
prediction follows a sequence-by-sequence (seqxseq) architecture (Sutskever et
al. (2014)). To ensure consistent scaling, all data points have been
normalized within the range of 0-1, utilizing the global minimum and maximum
values across the entire dataset. These normalized values are subsequently
transformed back to their original, real-world scale using the same minimum
and maximum values.

Propagation of predicted shear Reynolds stress in developed flow: In fully
developed turbulent flows, the momentum equation (2) in the longitudinal
direction can be expressed in dimensionless form using the characteristic
velocity scale $u_{\tau}$ and boundary layer thickness $\delta$, as mentioned
earlier. It takes the following form:

$\frac{dU^{+}}{dy^{+}}\;=\;\langle uv\rangle^{+}+(1-\frac{y^{+}}{R_{\tau}})$
(3)

where the pressure gradient has been substituted with the wall shear stress
(Tennekes and Lumley (1972)). Equation (3) illustrates the equilibrium between
the Reynolds shear stress and the wall-normal gradient of the longitudinal
mean velocity for a fully developed turbulent flow. This equilibrium is
depicted in the plots of $U^{+}$, $\langle uv\rangle^{+}$, and $dU^{+}/dy^{+}$
as functions of the dimensionless distance to the wall, $y^{+}$, in Figs. 3
and 4. Specifically, the agreement of the $U^{+}$ profile with the logarithmic
law of the wall in the logarithmic sub-region, shown in Fig. 3, indicates that
the predicted value of $\langle uv\rangle^{+}$ is accurate.

Thus, as observed in these figures, the mean velocity exhibits a logarithmic
sub-region with appropriate values for the given Reynolds number, indicating
the characteristic velocity profile of a developed turbulent flow.
Furthermore, the shear stresses, including both molecular and turbulent
components, demonstrate a state of perfect equilibrium under developed flow
conditions. This implies that the LSTM predictions successfully capture the
essential features of developed turbulent flows.

In order to evaluate the stability of the results obtained from the LSTM
model, a series of tests were conducted using flows with initial conditions
significantly distant from a fully developed state. The objective of these
tests was to assess the capability of the LSTM model to converge towards a
stable solution.

Propagation of predicted shear Reynolds stress in developing flows: In this
section the results of the LSTM model L110DL-Y propagated in developing
turbulent channel flow are presented. The first case is a blowing case in a
channel with a dimensionless injection velocity of $0.40$ (the set of data for
training included an injection case with a dimensionless injection velocity of
$0.60$). The results of the propagated LSTM predictions are shown in
comparison with $\kappa-\epsilon$ model and data from DNS, Figs. 5 to 7.

Figures 5 and 6 clearly demonstrate the LSTM model’s superior performance in
computing the skin friction when compared to the traditional $\kappa-\epsilon$
model. Additionally, Figure 7 illustrates that the propagated shear Reynolds
stress follows the general trend observed in the DNS data.

A RANS prediction depends on various factors, such as the numerical grid,
software precision, and, of course, the RANS model itself. In this study, all
these aspects were reasonably implemented. Based on the presented comparisons
between the LSTM RNN and the $\kappa-\epsilon$ model predictions, it can be
concluded that LSTM-RANS exhibits superior performance. However, it is
essential to note that the LSTM-RANS simulations with the L110DL-Y model used
here required approximately 50% more time for computation than
$\kappa-\epsilon$ simulations.

Regarding the time-processing aspect of LSTM-RANS, it is worth mentioning that
a more comprehensive study should be conducted. This broader investigation
should encompass the use of different RANS models and LSTM networks with
varying numbers of parameters to draw more definitive conclusions.

In summary, the study demonstrates that LSTM-RANS outperforms the traditional
$\kappa-\epsilon$ model, but further research is needed to explore its
efficiency and effectiveness fully.

##  6 Conclusions

In this study, we introduce the application of an artificial Recurrent Neural
Network (RNN) known as LSTM as a promising alternative to a turbulent RANS
model. The primary objective was to utilize the LSTM model for predicting the
shear Reynolds stress in both developed and developing turbulent channel
flows, and subsequently propagate this prediction.

To assess its performance, we conduct a comparative analysis, where we compare
the LSTM results, propagated through CFD simulations, with the outcomes from
the traditional $\kappa-\epsilon$ model and data obtained from DNS.
Remarkably, the LSTM approach demonstrates a strong performance in these
analyses, showing its potential as an effective technique for modeling the
shear Reynolds stress in turbulent flows.

It is essential to note that this study serves as a proof of concept, and the
validation process utilizes 2D statistical turbulent flows as the training and
test data sets. Despite this limitation, the results presented here indicate
the LSTM RNN’s efficacy in modeling the shear Reynolds stresses accurately in
turbulent flows.

Overall, the findings from this study open up exciting possibilities for
further exploring the capabilities of LSTM models as alternatives to
conventional RANS models in the context of turbulence prediction and modeling.

## References

  * Pope [2000] S. Pope.  _Turbulent Flow_.  University Press, Cambridge, Cambridge, UK, 1st. edition, 2000.
  * Milano and Koumoutsakos [2002] M. Milano and P. Koumoutsakos.  Neural network modeling for near wall turbulent flow.  _J. Comp. Physics_ , pages 1–12, 2002.  doi: 10.1006/jcph.2002.7146.
  * Lecun et al. [2015] Y. Lecun, Y. Bengio, and G. Hinton.  Deep learning.  _Nature_ , 521, 2015.  doi: 10.1038/nature14539.
  * Ling and Templeton [2015] J. Ling and J. Templeton.  Evaluation of machine learning algorithms for prediction of regions of high Reynolds-averaged Navier-Stokes uncertainty.  _Phys. Fluids_ , 27:99, 2015.
  * Ling et al. [2016a] J. Ling, Jones R., and J. Templeton.  Machine learning strategies for systems with invariance properties.  _J. Comput. Phys_ , pages 22–35, 2016a.
  * Ling et al. [2016b] J. Ling, Kurzawski A., and Templeton J. D.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _J. Fluid Mechanics_ , 807:155–166, 2016b.
  * Hochriter and Schmidhuber [1997] S. Hochriter and J. Schmidhuber.  Long short term memory.  _Neural Computation_ , 9:1–32, 1997.
  * Sutskever et al. [2014] I. Sutskever, Vinyals O., and Quoc V. Le.  Sequence to sequence learning with neural networks.  _arXiv_ , pages 1–9, 2014.  doi: org/10.48550/arXiv.1409.3215.
  * Pasinato [2012] H. D. Pasinato.  Dissimilarity of turbulent fluxes of momentum and heat in perturbed turbulent flows.  _J. Heat Transfer-ASME_ , 135:1–12, 2012.
  * Abadi et al. [2015] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Chen [2018] G. Chen.  A gentle tutorial of recurrent neural network with error backpropagation.  pages 1–10, 2018.
  * Tennekes and Lumley [1972] H. Tennekes and J. L. Lumley.  _A First Course in Turbulence_.  MIT Press, Cambridge USA, 1972.


# Numerical solution to the time-independent Gross-Pitaevskii equation

Tsogbayar Tsednee  Banzragch Tsednee  Tsookhuu Khinayat  Institute of Physics
and Technology, Mongolian Academy of Sciences, Peace Ave 54B, 13330
Ulaanbaatar, Mongolia

###### Abstract

We solve the time-independent Gross-Pitaevskii equation modeling the Bose-
Einstein condensate trapped in an anistropic harmonic potential using a
pseudospectral method. Numerically obtained values for an energy and a
chemical potential for the condensate with positive and negative scattering
length have been compared with those from the literature. The results show
that they are in good agreement when an atomic interaction is not too strong.

Bose-Einstein condensate, nonlinear Schrödinger equation, chemical potential,
rubidium, lithium,

###### pacs:

Valid PACS appear here

##  I Introduction

When the thermal de Broglie wavelength exceeds the mean spacing between
identical boson-particles, bosons are stimulated by presence of other bosons
in the lowest energy state to occupy that state as well, resulting in
macroscopic occupation of a single quantum state Bose24 ; Einstein25 . This
phenomenon is named the Bose-Einstein condensation and the condensate that
forms constitutes a macroscopic quantum mechanical object. This theoretical
prediction had been confirmed experimentally 70 years later, particularly for
${}^{87}\mbox{Rb}$ Anderson95 , ${}^{7}\mbox{Li}$ Bradley95  and
${}^{23}\mbox{Na}$ Davis95 . The vapors of alkali atoms employed in the
experiments are very dilute, so one can expect that the two-body collision
accounting for by the knowledge of the $s$-wave scattering length might be
dominate. This also implies that the Gross-Pitaevskii theory Pitaevskii61 ;
Gross61  for weakly interacting bosons can be suitable for the system which
can be simulated to be confined in an isotropic Edwards95 ; Ruprecht95  and an
anisotropic Dalfovo96 ; Baym96  traps.

In this work we solve the time-independent Gross-Pitaevskii equation (GPE) for
$N$ alkali atoms in an anisotropic trap. We compute the condensate wave
function at $T=0$ for bosons interacting through positive and negative
scattering lengths and obtain the chemical potential and energy as a function
of $N$. Numerical method we choose to solve the GPE is a pseudospectral
method, which we had applied successfully in the past Tsogbayar13 .

The paper is organized as follows. In Section 2, we show the formalism of the
Gross-Pitaevskii theory for the anisotropic trap. In Section 3 we give a brief
discussion of a pseudospectral approach for the 3D problem. In Section 4 we
present the numerical results for the two cases of positive ${}^{87}\mbox{Rb}$
and negative ${}^{7}\mbox{Li}$ scattering lengths. Then a conclusion follows.

##  II Gross-Pitaevskii theory for trapped bosons

The mean field theory for a dilute assembly of bosons at $T=0$ results in an
effective nonlinear Schrödinger equation for the condensate’s wave function.
This equation, the Gross-Pitaevskii, nonlinear Schrödinger equation for
condensed bosons has a form:

$\displaystyle i\hbar\frac{\partial\Psi(\mathbf{r},t)}{\partial
t}=\Big{(}-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{ext}+\frac{4\pi\hbar^{2}aN}{m}|\Psi(\mathbf{r},t)|^{2}\Big{)}\Psi(\mathbf{r},t).$
(1)

Here $\Psi(\mathbf{r},t)$ is the Bose-Einstein condensate (BEC) wave function,
(also called the order parameter), $m$ is the mass of boson, $V_{ext}$ is an
external confining potential (trap), $a$ is the $s$-wave scattering length and
$N$ is the number of bosons in the condensate.

A stationary solution $\Psi(\mathbf{r},t)=e^{i\mu t/\hbar}\psi(\mathbf{r})$
obeys

$\displaystyle\Big{[}-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{ext}(\mathbf{r})+V_{mf}(\mathbf{r})\Big{]}\psi(\mathbf{r})=\mu\psi(\mathbf{r}),$
(2)

where the mean-field (mf) potential is
$V_{mf}=\frac{4\pi\hbar^{2}aN}{m}|\psi(\mathbf{r})|^{2}$. Once this equation
is solved the chemical potential $\mu$ is known and the free energy can be
calculated using

$\displaystyle E=\mu-\frac{1}{2}\langle
V_{mf}\rangle=\int\psi(\mathbf{r})^{\ast}\Big{(}-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{ext}(\mathbf{r})+\langle
V_{mf}\rangle\Big{)}\psi(\mathbf{r})d\mathbf{r}.$ (3)

Since $\psi$ and $V_{mf}$ in equation (2) depend on each other, the GP
equation must be solved self-consistently. One first uses an initial guess for
the wave function $\psi$ to calculate $V_{mf}$ using equation (4). This value
is then employed in equation (2) to obtain a new $\psi$, which is then used to
calculate $V_{mf}$ again. This process is repeated until self-consistency is
reached.

In our calculation we use a following anisotropic harmonic oscillator
potential:

$\displaystyle
V_{ext}(x,y,z)=\frac{m}{2}\omega^{2}_{x}x^{2}+\frac{m}{2}\omega^{2}_{y}y^{2}+\frac{m}{2}\omega^{2}_{z}z^{2}.$
(4)

By introducing the standard lengths $a_{\perp}=(\hbar/m\omega_{\perp})^{1/2}$
and $a_{z}=(\hbar/m\omega_{z})^{1/2}$, we can rescale the spatial coordinate,
the energy, and the wave function as $\mathbf{r}=a_{\perp}\mathbf{r}_{1}$,
$E=\hbar\omega_{\perp}E_{1}$ and
$\psi(\mathbf{r})=\sqrt{N/a^{3}_{\perp}}\psi_{1}(\mathbf{r}_{1})$, with
$\omega_{x}=\omega_{y}=\omega_{\perp}$. Here the wave function $\psi_{1}$ is
normalized to $1$. With help of the introduced asymmetry parameter
$\lambda=\omega_{z}/\omega_{\perp}$ and the quantity $g=4\pi aN/a_{\perp}$,
the time-independent GP equation (2) can be written as:

$\displaystyle\Big{[}-\frac{1}{2}\nabla^{2}_{1}+\frac{(x^{2}_{1}+y^{2}_{1})}{2}+\frac{\lambda^{2}z^{2}_{1}}{2}+g|\psi(\mathbf{r}_{1})|^{2}\Big{]}\psi(\mathbf{r}_{1})=\mu\psi(\mathbf{r}_{1}),$
(5)

##  III Numerical prodecure

In our calculation we use the Legendre-pseudospectral method Tsogbayar13 . In
terms of this approach, function $\psi(r)$ can be expressed as:

$\displaystyle\psi(r)\approx\psi_{N_{r}}(r)=\sum^{N_{r}}_{i=0}\psi(r_{i})g_{i}[x(r)],$
(6)

where $g_{i}[x(r)]$ is a cardinal function given with

$\displaystyle
g_{i}[x(r)]=-\frac{1}{N_{r}(N_{r}+1)P_{N_{r}}(x_{i})}\frac{(1-x^{2})P^{{}^{\prime}}_{N_{r}}(x)}{x-x_{i}}$
(7)

and $g_{i}(x_{j})=\delta_{ij}$. Here $N_{r}$ is a number of grid point along
$r(x)=a+(1+x)(b-a)/2,\,x\in[-1,1]$ and $r\in[a,b]$ with a length parameters
$(a,b)$. Here the Legendre-Gauss-Lobatto grid pints $x_{i}$ are determined as
the roots of the first derivative of the Legendre polynomial $P_{N}(x)$ with
respect to $x$, $P^{\prime}_{N_{r}}(x_{i})=0,\,\,i=0,\ldots,N_{r}$. In the
approach, the Laplace operator $\nabla^{2}$ can be approximated with a
differentiation matrix $d_{ij}$ Tsogbayar13 . So, for the 3D calculation, we
can approximate $\nabla^{2}_{x}+\nabla^{2}_{y}+\nabla^{2}_{z}\approx
I_{zz}\otimes(d^{2}_{xx}\otimes I_{yy})+I_{zz}\otimes(I_{xx}\otimes
d^{2}_{yy})+(I_{xx}\otimes I_{yy})\otimes d^{2}_{zz}$. Here $I$ is the unit
matrix, and $\otimes$ expresses the Kroneckor (tensor) product Tsogbayar13 .
In our numerical calculation we use $a_{x}=a_{y}=a_{z}=a=-5$,
$b_{x}=b_{y}=b_{z}=b=5$ in units of $a_{\perp}$, and $N_{x}=N_{y}=N_{z}=32$.

##  IV Results and discussion

As an example of atoms with repulsive interaction, we choose
${}^{87}\mbox{Rb}$, as in the experiment of Ref. Anderson95 . In our
calculation, all values of the physical parameters are taken from Ref.
Dalfovo96 : the $s$-wave triplet-spin scattering length, $a_{s}=100a_{0}$
where $a_{0}$ is the Bohr radius; the asymmetry parameter of the experimental
trap is $\lambda=\omega_{z}/\omega_{\perp}=\sqrt{8}$; the axial frequency
$\omega_{z}/2\pi=220\,Hz$; the corresponding characteristic length is
$a_{\perp}=1.222\times 10^{-12}\,cm$ and the ratio between the scattering and
the oscillator lengths is $a/a_{\perp}=4.33\times 10^{-3}$. In our calculation
number of grid points is $N_{x}=N_{y}=N_{z}=24$, and results are independent
on this number. Table 1 shows the excess chemical potential and energy per
particle for three values of $N=100,1000$ and $1000$, and our calculated
values are close to those in Ref. Dalfovo96 , which had been obtained with a
direct minimization approach combined with an imaginary time technique. Both
quantities are expressed in units of $\hbar\omega_{\perp}$.

Table 1:  Results for the ground state of ${}^{87}\mbox{Rb}$ atoms in a trap with $\lambda=\sqrt{8}$. Chemical potential and energy are in unit of $\hbar\omega_{\perp}$. A number of grid point is $N_{x}=N_{y}=N_{z}=24$.  | $N=100$ | $N=1000$ | $N=5000$
---|---|---|---
| $\mu$ | $E$ | $\mu$ | $E$ | $\mu$ | $E$
This work | 2.88 | 2.67 | 4.77 | 3.84 | 8.15 | 6.13
Dalfovo96  | 2.88 | 2.66 | 4.77 | 3.84 | 8.14 | 6.12

In Figure 1 we show plots of the wave function along the $x$ (panel a) and the
$z$ axis (panel b) for four values of $N$. When $N$ increases the repulsion
among the atoms tends to lower the central density, and expands the cloud of
the atoms towards region where the trapping potential is higher. This results
in increase of an energy per particle.

Figure 1:  The ground-state wave function for ${}^{87}\mbox{Rb}$ along the $x$
axis (a) and along the $z$-axis (b). Distances are in units of $a_{\perp}$.
Blue, green, red and black lines corresponds to $N=100,500,2000$ and $5000$,
in descending order of central density. Cross are taken from Ref. Dalfovo96 .

As an example of atoms with an attractive interaction we choose
${}^{7}\mbox{Li}$, as in the experiment of Ref. Bradley95 . In the
calculation, we use $a_{s}=-27\,a_{0}$; $\omega_{z}=2\pi\times 117$ Hz;
$a_{\perp}=2.972\times 10^{-4}\,$cm; $|a|/a_{\perp}=0.48\times 10^{-3}$;
$\omega_{\perp}=2\pi\times 163\,$Hz, and
$\lambda=\omega_{z}/\omega_{\perp}=0.72$. Table 2 shows numerical values of
chemical potential and energy for the ${}^{7}\mbox{Li}$ ground state in unit
of $\hbar\omega_{\perp}$ for three values of $N$.

Table 2: Same results as shown in Table 1, but for ${}^{7}\mbox{Li}$ atoms in a trap with $\lambda=0.72$. | $N=100$ | $N=500$ | $N=1000$
---|---|---|---
| $\mu$ | $E$ | $\mu$ | $E$ | $\mu$ | $E$
This work | $1.33$ | $1.34$ | $1.17$ | $1.27$ | $0.86$ | $1.16$

Figure 2a presents the ground state wave function for the ${}^{7}\mbox{Li}$
atom for three values of $N$. In this plot, the central density of cloud
increases rapidly with $N$ since more attractive potential energy is added.
Fig. 2b shows two-dimensional decsription of ground state wave function
$\psi(x,y,0)$ for $N=1000$.

Figure 2:  Same plots as shown in Figure 1, but for the ${}^{7}\mbox{Li}$
ground state along the $x$ axis (a), and blue, red and black lines correspond
to $N=100,500$ and $1000$, in ascending order of central density (a). Panel
(b) shows an interpolated two-dimensional plot for $N=1000$. Crosses are taken
from Ref. Dalfovo96 .

##  V Conclusion

In this paper we have solved the time-independent Gross-Pitaevskii equation,
non-linear eigenvalue problem, for a dilute gas of alkali atoms in an
anisotropic traps using the pseudospectral method. The ground state wave
function for the condensate with repulsive $(^{87}\mbox{Rb})$ and attractive
(${}^{7}\mbox{Li}$) behaviors at $T=0$ has been obtained and natures of these
wave function depending on number of particle have been discussed. Chemical
potential and energy per particle for the $(^{87}\mbox{Rb})$ and
(${}^{7}\mbox{Li}$) condensates for different values of $N$ have been
presented as well.

## References

## References

  * (1) S.N. Bose, Planck’s Law and Light Quantum Hypothesis, Z. Phys., 26 (2004), 178-181.
  * (2) A. Einstein, Quantentheorie des einatomigen idealen Gases, Sitzber. Kgl. Preuss. Akad. Wiss. 1925 (1925), 3-14.
  * (3) M.H. Anderson, J. R. Ensher, M.R. Matthews, C.E. Wieman, and E.A. Cornell, Observation of Bose-Einstein Condensation in a Dilute Atomic Vapor, Science 269, (1995), 198-201.
  * (4) C.C. Bradley, C.A. Sackett, J.J. Tollett, and R.G. Hulet, Evidence of Bose-Einstein Condensation in an Atomic Gas with Attractive Interactions, Phys. Rev. Lett. 75, (1995), 1687-1690.
  * (5) K.B. Davis, M.-O. Mewes, M.R. Andrews, N.J. van Druten, D.S. Durfee, D.M. Kurn, and W. Ketterle, Bose-Einstein Condensation in a Gas of Sodium Atoms, Phys. Rev. Lett. 75, (1995), 3969-3973.
  * (6) L.P. Pitaevskii, Vortex lines in the imperfect Bose-gas, Zh. Eksp. Teor. Fiz. 40, (1961), 646-649.
  * (7) E.P. Gross, Structure of a quantized vortex in boson systems, Nuovo Cimento 20, (1961), 454-477.
  * (8) M. Edwards and K. Burnett, Numerical solution of the nonlinear Schrödinger equation for small samples of trapped neutral atoms, Phys. Rev. A 51, (1995), 1382-1386.
  * (9) P.A. Ruprecht, M.J. Holland, K. Burnett, and M. Edwards, Time-dependent solution of the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms, Phys. Rev. A 51, (1995), 4704-4711.
  * (10) F. Dalfovo and S. Stringari, Bosons in anisotropic traps: Ground state and vortices, Phys. Rev. A 53, (1996), 2477-2485.
  * (11) G. Baym and C. Pethick, Ground-State Properties of Magnetically Trapped Bose-Condensed Rubidium Gas, Phys. Rev. Lett. 76, (1996), 6-9.
  * (12) Ts. Tsogbayar, M Horbatsch, Calculation of Stark resonance parameters for the hydrogen molecular ion in a static electric field, J. Phys. B 46, (2013), 085004 (8pp).