# Optical ladder operators in the Glauber-Fock oscillator array

I. Bocanegra<EMAIL_ADDRESS>Instituto Nacional de Astrofísica
Óptica y Electrónica, Calle Luis Enrique Erro No. 1
Santa María Tonantzintla, Puebla, 72840, Mexico  L. Hernández-Sánchez
Instituto Nacional de Astrofísica Óptica y Electrónica, Calle Luis Enrique
Erro No. 1
Santa María Tonantzintla, Puebla, 72840, Mexico  I. Ramos-Prieto  Instituto
Nacional de Astrofísica Óptica y Electrónica, Calle Luis Enrique Erro No. 1
Santa María Tonantzintla, Puebla, 72840, Mexico  F. Soto-Eguibar  Instituto
Nacional de Astrofísica Óptica y Electrónica, Calle Luis Enrique Erro No. 1
Santa María Tonantzintla, Puebla, 72840, Mexico  H. M. Moya-Cessa  Instituto
Nacional de Astrofísica Óptica y Electrónica, Calle Luis Enrique Erro No. 1
Santa María Tonantzintla, Puebla, 72840, Mexico

###### Abstract

In this study, we investigate the stationary states of the Glauber-Fock
oscillator waveguide array. We begin by transforming the associated
Hamiltonian into the form of a quantum harmonic oscillator Hamiltonian,
allowing the implementation of a supersymmetric (SUSY) approach. By
considering the simplest case for the intertwining operator, the optical
ladder operators are straightforwardly constructed and shown to map
eigensolutions into eigensolutions of the corresponding Hamiltonian operator,
in pretty much the same manner as it is done for the quantum harmonic
oscillator case. The ladder of the corresponding (eigen) supermodes is then
easily established.

##  I Introduction

From the point of view of Supersymmetric (SUSY) quantum mechanics [1, 2, 3, 4,
5, 6], if there exists an operator $B$, intertwining the Hamiltonians $H$ and
$\tilde{H}$, this is, satisfying

$BH=\tilde{H}B,$ (1)

then the solutions of the Schrödinger equation associated with $\tilde{H}$ can
be obtained from those corresponding to $H$ and vice versa: the operator $B$
is then called the intertwiner. Actually, the intertwining relation (1) is
often encountered in connection with the so-called Darboux transformation [7,
8, 9, 10, 11, 12] and the factorization method [13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23]. Consequently, supersymmetry, the Darboux transformation and
the factorization method are sometimes treated as equivalent [8] (and are
actually equivalent under certain conditions).

The relation (1) has been exploited to construct (families of) exactly
solvable Hermitian [17, 19, 20] and non-Hermitian Hamiltonians [5, 6, 22, 23],
in the time-independent [2, 8, 11, 21] as well as time-dependent cases [8, 4,
12]. In particular, in the stationary regime, if $B=a$ ($B=a^{\dagger}$), with
$a$ ($a^{\dagger}$) the annihilation (creation) operator of the harmonic
oscillator, and $H=n$, with $n=a^{\dagger}a$ the usual number operator, then
$\tilde{H}=n+1$ ($\tilde{H}=n-1$) [24]. This is, indeed, the basis for the
construction of the solution of the eigenvalue equation of the quantum
harmonic oscillator developed early by Schrödinger [13, 14], and made popular
by Dirac in his book [15].

Besides, due to the formal equivalence between the Schrödinger equation and
the paraxial Helmholtz equation [25] (also between the stationary Schrödinger
equation and the Helmholtz equation), optical waveguides are suitable devices
to observe, study and test quantum phenomena [26]. Therefore, either isolated
waveguides or waveguide arrays (called optical lattices), are susceptible to
be transformed by means of Darboux or SUSY transformations [27, 28, 29, 30,
31, 32, 33, 34, 35], in the Hermitian [27, 28, 29, 30, 32] and non-Hermitian
[31, 34, 35, 22, 23] regimes. Specifically, an optical lattice associated with
a Hamiltonian of the type of the quantum harmonic oscillator can be
“intertwined” with itself, giving rise to the ladder of eigenstates
(supermodes [36, 37]).

In the present work we study the stationary regime of a very particular (semi-
infinite) waveguide array, referred to as the Glauber-Fock (oscillator) array
[38, 39] (also see [40, 41, 36]), and characterized by a non-uniform distance
between adjacent waveguides, as well as a gradient of refractive index
increasing with the waveguide site $k$, $k=0,1,\dots$ (see Fig. 1). The
corresponding Hamiltonian $H$ can be taken to the form of the quantum harmonic
oscillator Hamiltonian, by means of a rather elementary transformation [36].
Then, by considering the simplest case for the intertwiner $B=a$, the
operators mapping eigensolutions into eigensolutions of $H$ are
straightforward to obtain, and the corresponding ladder of supermodes is
easily constructed.

With that in mind, the objective of this paper is twofold. On the one hand,
the optical ladder operators associated with $H$ are constructed in a formal
way, and interpreted as switchers between the stationary (eigen) supermodes of
the Glauber-Fock oscillator array. On the other hand, we set the precedent for
further constructions of (Hermitian and non-Hermitian) waveguide arrays that
can be obtained by considering more general intertwining operators $B$ (in the
stationary as well as the time-dependent regimes).

The general outline is as follows: in section II the system under study is
introduced, as well as the transformation to turn the corresponding
Hamiltonian into the form of the quantum harmonic oscillator Hamiltonian. In
section III, the generic SUSY approach is presented, for an arbitrary
intertwiner operator $B$. In turn, in section IV, the foundations presented in
section III are implemented for the simplest case of the intertwining
operator. This results in the optical ladder operators for the stationary
eigen supermodes of the Glauber-Fock oscillator array. Finally, in section V
the main conclusions are drawn.

Figure 1: In the Glauber-Fock oscillator array, the waveguides are separated
in a non-uniform way, $\Delta x$ is the distance between the first two
waveguides. As the site $k$ of the waveguide increases, the waveguides become
closer and closer, as shown. In addition, as $k$ grows the refractive index of
each waveguide increases as well. Such gradient in refractive index is
schematically shown with darker tones of red as $k\to\infty$.

##  II Glauber-Fock oscillator array

In the Glauber-Fock oscillator array [38, 39], the amplitude of the electric
field propagating in the $k$-th waveguide, $k=0,1\dots$, and here denoted
$c_{k}(z)$, is ruled by

$\operatorname{\textrm{i}}\dot{c}_{k}+\omega
kc_{k}+g\left(\sqrt{k}c_{k-1}+\sqrt{k+1}c_{k+1}\right)=0,$ (2)

where $\omega\in\mathbb{R}$ is the propagation constant of the first
waveguide, and $g\in\mathbb{R}$ is the coupling between the first to
waveguides; $g$ is proportional to $e^{-\Delta x}$, with $\Delta x$ the
distance between the first two sites (see Fig. 1). Equation (2) can be
alternatively written as

$\operatorname{\textrm{i}}\partial_{z}\ket{\psi(z)}=H\ket{\psi(z)},$ (3)

where

$\ket{\psi(z)}=\sum_{k=0}^{\infty}c_{k}(z)\ket{k},$ (4)

with $c_{j}=\braket{j}{\psi(z)}$ for $j=0,1,\dots$. The set
$\left\\{|k\rangle\right\\}_{k=0,1,\dots}$ represents the Fock basis, which
spans the Hilbert space $\mathcal{H}$. The Hamiltonian in (3) is given by

$H=-\omega n-g(a^{\dagger}+a).$ (5)

The state $\ket{\psi(z)}$ in (4) contains the information of the total
electric field in the array for each $z$. As the Hamiltonian $H$ is time-
independent, the mathematical solution of (3), is simply

$\ket{\psi(z)}=e^{-\operatorname{\textrm{i}}Hz}\ket{\psi(0)}.$ (6)

Here we are particularly interested in the initial condition

$\ket{\psi(0)}=\ket{\psi_{\ell}},\qquad\ell=0,1,\dots,$ (7)

with $\ket{\psi_{\ell}}$ satisfying the eigenvalue equation

$H\ket{\psi_{\ell}}=E_{\ell}\ket{\psi_{\ell}},$ (8)

and such that the solution (6) is stationary, namely,

$\ket{\psi(z)}=e^{-\operatorname{\textrm{i}}Hz}\ket{\psi_{\ell}}=e^{-\operatorname{\textrm{i}}E_{\ell}z}\ket{\psi_{\ell}}.$
(9)

By making the transformations

$\ket{\psi(z)}=D^{\dagger}\left(\frac{g}{\omega}\right)\ket{w(z)},\qquad\ket{\psi_{\ell}}=D^{\dagger}\left(\frac{g}{\omega}\right)\ket{w_{\ell}},$
(10)

where $D(\xi)=\exp(\xi a^{\dagger}-\xi^{*}a)$, $\xi\in\mathbb{C}$, is the
Glauber displacement operator [24], equation (8) and the solution (9) turn,
respectively, into

$\bar{H}\ket{w_{\ell}}=E_{\ell}\ket{w_{\ell}},\qquad\ket{w(z)}=e^{-\operatorname{\textrm{i}}E_{\ell}z}\ket{w_{\ell}},$
(11)

where

$\bar{H}=DHD^{\dagger}=-\omega n+\frac{g^{2}}{\omega},$ (12)

is diagonal. The shortcut notation $D=D\left(\frac{g}{\omega}\right)$,
$D^{\dagger}=D^{\dagger}\left(\frac{g}{\omega}\right)$ is used from now on.
The Hamiltonian (12) has the basic form of the quantum harmonic oscillator
Hamiltonian, therefore it is susceptible to be supersymmetrically transformed.
Next, we set the foundations for a generic SUSY transformation.

##  III Generic SUSY transformation

If one considers an arbitrary operator $B$ intertwining $\bar{H}$ in (12) with
some other Hamiltonian $\tilde{H}$, this is

$B\bar{H}=\tilde{H}B,$ (13)

the equalities in (11) are transformed as

$\tilde{H}\ket{\phi_{\ell}}=E_{\ell}\ket{\phi_{\ell}},\qquad\ket{\phi(z)}=e^{-\operatorname{\textrm{i}}E_{\ell}z}\ket{\phi_{\ell}},$
(14)

where we have defined

$\ket{\phi(z)}=\kappa B\ket{w(z)},\qquad\ket{\phi_{\ell}}=\kappa
B\ket{w_{\ell}},$ (15)

with $\kappa\in\mathbb{C}$ a normalization constant. Therefore, from the first
equality in (12), (13) becomes the intertwining relation between $H$ and
$\tilde{H}$:

$AH=\tilde{H}A,\qquad A=BD,$ (16)

and finally, the solution $\ket{\phi_{\ell}}$ in the second expression of (15)
can be written in terms of the solution $\ket{\psi_{\ell}}$ of the eigenvalue
equation (8), as

$\ket{\phi_{\ell}}=\kappa A\ket{\psi_{\ell}}.$ (17)

In turn, the stationary evolution is given by the second expression in (14).

Therefore, in order to perform the supersymmetric transformation, the
solutions of the eigenvalue equation (8) must be known. The obtention of both
the set of eigenvalues and eigenvectors of (8), $E_{\ell}$ and
$\ket{\psi_{\ell}}$, respectively, are given in the Appendix A.

In what follows, the simplest case for the intertwiner $B$ is considered, the
Hamiltonian $H$ is “intertwined” with itself, giving rise to the
straightforward construction of optical ladder operators connecting solutions
of (8), as well as the corresponding ladder of eigen supermodes.

##  IV Construction of the ladder operators

In this section, the case $B=a$ ($B=a^{\dagger}$) is studied, as we are
interested in the construction of the ladder of (eigen) supermodes associated
with the Hamiltonian (5). Nevertheless, the present discussion is intended to
set the precedent for more general choices of the intertwiner $B$, as
explained in Section III. From (12), it is easy to prove that

$a\bar{H}=(\bar{H}-\omega)a,\qquad
a^{\dagger}\bar{H}=(\bar{H}+\omega)a^{\dagger}.$ (18)

From the point of view of SUSY quantum mechanics [compare (13) with both
expressions in (18)], this means that $a$ ($a^{\dagger}$) intertwines
$\bar{H}$ with itself (see Ref. [16]). From both expressions in (18), and by
using the first equality in (12), we obtain, respectively

$(D^{\dagger}aD)H=(H-\omega)(D^{\dagger}aD)$ (19)

and

$(D^{\dagger}a^{\dagger}D)H=(H+\omega)(D^{\dagger}a^{\dagger}D),$ (20)

where $D^{\dagger}aD=a+\frac{g}{\omega}$ and
$D^{\dagger}a^{\dagger}D=a^{\dagger}+\frac{g}{\omega}$ are the optical ladder
operators connecting the stationary (eigen) supermodes of $H$. By applying the
operator $D^{\dagger}aD$ to the left of equation (8), it is obtained

$H\left[D^{\dagger}aD\ket{\psi_{\ell}}\right]=(E_{\ell}+\omega)\left[D^{\dagger}aD\ket{\psi_{\ell}}\right],$
(21)

where

$(D^{\dagger}aD)\ket{\psi_{\ell}}=e^{\operatorname{\textrm{i}}\omega
z}\sqrt{\ell}\ket{\psi_{\ell-1}}.$ (22)

Similarly, by acting with $D^{\dagger}a^{\dagger}D$ on the left of (8), we
obtain

$H\left[D^{\dagger}a^{\dagger}D\ket{\psi_{\ell}}\right]=(E_{\ell}-\omega)\left[D^{\dagger}a^{\dagger}D\ket{\psi_{\ell}}\right],$
(23)

where

$(D^{\dagger}a^{\dagger}D)\ket{\psi_{\ell}}=e^{-\operatorname{\textrm{i}}\omega
z}\sqrt{\ell+1}\ket{\psi_{\ell+1}}.$ (24)

From expressions (22) and (24) it can be seen that the operators
$D^{\dagger}a^{\dagger}D$ and $D^{\dagger}a^{\dagger}D$ are indeed ladder
operators, switching between eigensolutions of (8), and therefore give the
stationary evolution dictated in (9). The ladder of supermodes can be
straightforwardly constructed. Similar to the quantum harmonic oscillator
case, the $\ell$-th eigenstate (supermode) of the ladder, can be written in
terms of the lowest state $\ket{\psi_{0}}$, as

$\ket{\psi_{\ell}}=\frac{e^{\operatorname{\textrm{i}}\omega\ell
z}}{\sqrt{\ell!}}D^{\dagger}(a^{\dagger})^{\ell}D\ket{\psi_{0}}.$ (25)

Figure 2 shows the effect of the optical annihilation and creation operators
(22) and (24), respectively, for some specific values of the parameters. By
departing from the normalized eigen supermode $\ket{\psi_{4}}$ in Fig. 2
(middle), the (normalized) state $\ket{\psi_{3}}$ in Fig. 2 (up) is obtained
by means of the application of the operator (22). In turn, again departing
form $\ket{\psi_{4}}$, the normalized eigen supermode $\ket{\psi_{5}}$ in Fig.
2 (down) is reached through the operator (24).

Figure 2: Effect of the optical ladder operators (22) and (24) on a given
initial eigen supermode of the Glauber-Fock oscillator array. The chosen
parameters are: $\omega=1.3$, $g=0.29$. By departing from the normalized eigen
supermode $\ket{\psi_{4}}$ (middle), the (also normalized) state
$\ket{\psi_{3}}$ (up) can be obtained by means of the optical annihilation
operator (22), while the normalized state $\ket{\psi_{5}}$ (down) is obtained
through the optical creation operator (24).

In addition, Fig. 3 (up) shows the normalized distributions $c_{k}$
corresponding to the three eigenstates ($\ell=3,4,5$) shown in Fig. 2. The
green, red and blue lines correspond, respectively, to $\ket{\psi_{3}}$ [Fig.
2 (up)], $|\psi_{4}\rangle$ [Fig. 2 (middle)] and $\ket{\psi_{5}}$ [Fig. 2
(down)]. It can be observed that, as $\ell$ increases, the distribution
$c_{k}$ becomes smaller and broader, in agreement with Fig. 2. In turn, Fig. 3
(down) shows $c_{k}$ for a slightly greater coupling $g=0.36$. It is seen
that, as the coupling $g$ grows, the distributions become broader as well.
This can be particularly appreciated in the $\ket{\psi_{5}}$ supermode (blue)
in Fig. 3 (down). In both cases the solid lines correspond to the analytical
results while the the stars correspond to the numerical solutions.

Figure 3: Up: Distributions $c_{k}$ corresponding to the eigen supermodes
$\ket{\psi_{3}}$ (green), $\ket{\psi_{4}}$ (red) and $\ket{\psi_{5}}$ (blue)
in Fig. 2 (up), Fig. 2 (middle) and Fig. 2 (down), respectively. As the reader
can appreciate, as $\ell$ increases ($\ell=3,4,5$), the distribution becomes
smaller and broader. This is in agreement with Fig. 2. Down: Distributions
$c_{k}$ for $g=0.36$. As can be seen, as the coupling $g$ grows, the
distributions become broader as well. This can be particularly observed in the
eigen supermode $\ket{\psi_{5}}$ (blue). The solid lines correspond to the
analytical results while the the stars correspond to the numerical solutions.

##  V Conclusions

The optical ladder operators for the Glauber-Fock oscillator waveguide array
were constructed in the stationary regime. They switch between eigen
supermodes, as expected. In particular, when $g$ is relatively smaller, they
switch between the excited waveguide, as shown in Figure 2. Therefore, the
ladder of eigen supermodes is straightforwardly constructed by departing from
the lowest state, in much the same manner as in quantum harmonic oscillator
case.

In addition, as the ladder operators were constructed by departing from a
supersymmetric (SUSY) approach, a different choice of the intertwiner might
lead to new configurations of waveguide arrays. These might possibly be
associated with non-Hermitian or even with time-dependent Hamiltonians. Work
in these direction is already in progress.

##  Appendix A Solution of the initial eigenvalue equation

By using Eq. (3.1.12) of Ref. [24], from (6) and (12), it is straightforward
to obtain

$\ket{\psi(z)}=D^{\dagger}\left(\frac{g}{\omega}\right)e^{\operatorname{\textrm{i}}z(\omega
n-\frac{g^{2}}{\omega})}D\left(\frac{g}{\omega}\right)\ket{\psi(0)}.$ (26)

By inspection it is seen that, by choosing
$\ket{\psi(0)}=D^{\dagger}\left(\frac{g}{\omega}\right)\ket{\ell}$, with
$\ket{\ell}$ a Fock state satisfying

$n\ket{\ell}=\ell\ket{\ell},$ (27)

the stationary solution $\ket{\psi(z)}=\ket{\psi_{\ell}}$, with

$\ket{\psi_{\ell}}=e^{\operatorname{\textrm{i}}z(\omega
n-\frac{g^{2}}{\omega})}D^{\dagger}\left(\frac{g}{\omega}\right)\ket{\ell},$
(28)

is obtained. In turn, by acting with $H$ on the state (28), the corresponding
spectrum

$E_{\ell}=-\omega\ell+\frac{g^{2}}{\omega},$ (29)

is reached.

## Acknowledgments

I. Bocanegra acknowledges CONAHCyT (México) for financial support through the
postdoctoral fellowship 711878 and projects A1-S-24569 and CF19-304307. He is
also grateful to IPN (México) for supplementary economical support through the
project SIP20232237. L. Hernández Sánchez also thanks the Instituto Nacional
de Astrofísica, Óptica y Electrónica (INAOE) and the Consejo Nacional de
Humanidades, Ciencias y Tecnologías (CONAHCyT) for the PhD scholarship awarded
(No. CVU: 736710).

## References

  * [1] Witten E. “Dynamical breaking of supersymmetry”, Nucl. Phys. B 188, 513–554 (1981).
  * [2] Andrianov A. A., Borisov N. V., Ioffe M. V. and Eides M. I., “Supersymmetric mechanics: a new look at the equivalence of quantum systems”, Theor. Math. Phys. 61, 965–972 (1984).
  * [3] Cooper V., Khare A. and Sukhatme U., “Supersymmetry and quantum mechanics”, Phys. Rep. 251 267–385 (1995).
  * [4] Bagrov V. G. and Samsonov B., “Supersymmetry of a nonstationary Schrödinger equation”, Phys. Lett. A 210, 60–64 (1996).
  * [5] Zelaya K., Cruz y Cruz S. and Rosas-Ortiz O., “On the construction of non-Hermitian Hamiltonians with all-real spectra through supersymmetric algorithms”, on P. Kielanowsky A. Odzijewicz E. Previato Geometric methods in Physics XXXVIII, Trends in Mathematics, Springer WS AG2020, 283–292 (2020).
  * [6] Cruz y Cruz S., Romero-Osnaya A. and Rosas-Ortiz O., “Balanced gain-and-loss optical waveguides. Exact solutions for guided modes in Susy-QM”, Symmetry 13, 1583 (2021).
  * [7] Darboux G. “Sur une proposition relative aux équations linéaires”, C. R. Acad. Sci. Paris 94 1456–9 (1882).
  * [8] Bagrov V. G. and Samsonov B. F., “Darboux transformation, factorization and supersymmetry in one-dimensional quantum mechanics”, Theor. Math. Phys. 104, 1051–1060 (1995).
  * [9] Bagrov V. G., Samsonov B. F. and Shekoyan L. A., “Darboux transformation for the nonsteady Schrödinger equation”, Russian Physics Journal 38, 706–712 (1995).
  * [10] Bagrov V. G. and Samsonov B. F., “Darboux transformation of the Schrödinger equation”, Phys. Part. Nucl. 28 374–397 (1997).
  * [11] Fernández D. J., Nieto L. M. and Rosas-Ortiz O., “Distorted Heisenberg algebra and coherent states for isospectral oscillator Hamiltonians”, J. Phys. A.: Math. Gen. 28 2693–2708 (1995)
  * [12] Cruz y Cruz S., Razo R. Rosas-Ortiz O. and Zelaya K., “Coherent states for exactly solvable time-dependent oscillators generated by Darboux transformations”, Phys. Scr. 95 044009-1–26 (2020).
  * [13] Schrödinger E., “A method of determining quantum-mechanical eingenvalues and eigenfunctions”, Proc. R. Ir. Acad. A 46, 9–16 (1940).
  * [14] Schrödinger E., “Further studies on solving eigenvalue problems by factorization”, Proc. R. Ir. Acad. A 46, 183–206 (1940).
  * [15] Dirac P. A. M., “The principles of quantum mechanics”, 2nd Ed., Oxford: Clarendon (1935).
  * [16] Infeld L. and Hull T. E., “The factorization method”, Rev. Mod. Phys. 23, 21–68 (1951).
  * [17] Mielnik B., “Factorization method and new potentials with the oscillator spectrum”, J. Math. Phys. 25, 3387 (1984).
  * [18] Mielnik B. and Rosas-Ortiz O., “Factorization: little or great algorithm”, J. Phys. A.: Math. Gen. 37 10007–10035 (2004).
  * [19] Kuru S. and Negro J., “Factorization of one-dimensional classical systems”, Ann. Phys. 323 413–431 (2008).
  * [20] Cruz y Cruz S., Kuru S. and Negro J., “Classical motion and coherent states for Pöschl-teller potentials”, Phys. Lett. A: 372 1391–1405 (2008).
  * [21] Cruz y Cruz S. and Santiago-Cruz C., “Position dependent mass Scarf Hamiltonians generated via the Riccati equation”, Math. Meth. Appl. Sci. 42, 4909–4924 (2019).
  * [22] Bocanegra I. and Cruz y Cruz S., “Classes of balanced gain-and-loss waveguides as non-Hermitian potential hierarchies”, Symmetry 14, 432 (2022).
  * [23] Bocanegra I., “New families of complex hyperbolic-secant refractive-index profiles through the factorization method”, J. Phys.: Conf. Ser. 2448, 012015 (2023).
  * [24] Louisell W. H., “Quantum statistical properties of radiation”, Wiley: New York, (1990).
  * [25] Marte M. A. M. and Stenholm S., “Paraxial light and atom optics: The optical Schrödinger equation and beyond”, Phys. Rev. A 56, 2940–2953 (1997).
  * [26] Longhi S., “Quantum-optical analogies using photonic structures”
  * [27] Miri M. A., Heinrich M. El-Ganainy R. and Christodoulides D. N., “Supersymmetric optical structures”, Phys. Rev.Lett. 110, 233902 (2013).
  * [28] Miri M. A., Heinrich M. and Christodoulides D. N., “SUSY-inspired one-dimensional transformation optics”, Optica 1, 89–95 (2014).
  * [29] Heinrich M., Miri M. A., Stützer S., El-Ganainy R., Nolte S., Szameit A. and Christodoulides D. N., “Supersymmetric mode converters”, Nat. Commun. 5, 1–7 (2014).
  * [30] Heinrich M., Miri M. A., Stützer S., Nolte S., Christodoulides D. N. and Szameit A., “Observation of supersymmetric scattering in photonic lattices”, Opt. Lett. 39, 6130–6133 (2014).
  * [31] Midya B., “Supersymmetry-generated one-way-invisible PT-symmetric optical crystals”, Phys. Rev. A 89, 032116 (2014).
  * [32] Zúñiga-Segundo A., Rodríguez-Lara B. M., Fernández D. J. and Moya-Cessa H. M., “Jacobi photonic lattices and their SUSY partners”, Opt. Exp. 22, 987–994 (2014).
  * [33] El-Ganainy R., Ge L., Khajavikhan M. and Christodoulides D. N., “Supersymmetric laser arrays”, Phys. Rev. A 92 033818 (2015).
  * [34] Contreras-Astorga A. and Jakubský V., “Photonic systems with two-dimensional landscapes of complex refractive index via time-dependent supersymmetry”, Phys. Rev. A 99, 053812 (2019).
  * [35] Queraltó G., Kremer M., Maczewsky L. J., Heinrich M., Mompart J., Ahufinger V. and Szameit A., “Topological state engineering via supersymmetric transformations”, Communication Physics 3, 49 (2020).
  * [36] Bocanegra I. and Moya-Cessa H. M., “Non-Hermitian transport in Glauber-Fock optical lattices”, https://doi.org/10.48550/arXiv.2307.06425 (2023).
  * [37] Bocanegra I. and Moya-Cessa H. M., “Non-Hermitian propagation in equally-spaced waveguide arrays”, https://doi.org/10.48550/arXiv.2307.06952 (2023).
  * [38] Perez-Leija A., Keil R., Szameit A., Abouraddy A. F., Moya-Cessa H. M. and Christodoulides D. N., “Tailoring the correlation and anticorrelation behavior of path-entangled photons in Glauber-Fock oscillator lattices”, Phys. Rev A 85, 013848 (2012).
  * [39] Keil R., Perez-Leija A., Aleahmad P., Moya-Cessa H., Nolte S., Christodoulides D. N. and Szameit A., “Observation of Bloch-like revivals in semi-infinite Glauber-Fock photonic lattices”, Opt. Lett. 37 (18), 3801–3 (2012).
  * [40] Perez-Leija A., Moya-Cessa H. M., Szameit A. and Christodoulides D. N., “Glauber-Fock photonic lattices”, Opt. Lett. 35 (14), 2409–2411 (2010).
  * [41] Keil R., Perez-Leija A., Dreisow F., Heinrich M., Moya-Cessa H., Nolte S., Christodoulides D. N. and Szameit A., “Classical analogue of displaced Fock states and quantum correlations in Glauber-Fock photonic lattices”, Phys. Rev. Lett. 107, 103601-1–5 (2011).


# Localized states with non-trivial symmetries: localized labyrinthine
patterns

M.G. Clerc  Departamento de Física and Millennium Institute for Research in
Optics, FCFM, Universidad de Chile, Casilla 487-3, Santiago, Chile  S.
Echeverría-Alar  Departamento de Física and Millennium Institute for Research
in Optics, FCFM, Universidad de Chile, Casilla 487-3, Santiago, Chile  M.
Tlidi  Faculté des Sciences, Université libre de Bruxelles (U.L.B), CP. 231,
1050 Brussels, Belgium

###### Abstract

The formation of self-organized patterns and localized states are ubiquitous
in Nature. Localized states containing trivial symmetries such as stripes,
hexagons, or squares have been profusely studied. Disordered patterns with
non-trivial symmetries such as labyrinthine patterns are observed in different
physical contexts. Here we report stable localized disordered patterns in
spatially extended dissipative systems. These 2D and 3D localized structures
consist of an isolated labyrinth embedded in a homogeneous steady state. Their
partial bifurcation diagram allows us to explain this phenomenon as a
manifestation of a pinning-depinning transition. We illustrate our findings on
Swift-Hohenberg type of equations and other well-established models for plant
ecology, nonlinear optics, and reaction-diffusion systems.

Spatiotemporal patterning resulting from a symmetry-breaking instability is a
central issue in almost all driven far from equilibrium systems [1, 2, 3].
Localized structures, dissipative solitons, and localized patterns belong to
this field of research. They consist of one or more regions in one state
surrounded by a region in a qualitatively different state [5, 6, 4, 7, 8].
Spatial localization appears not only in nonlinear systems but can occur in
linear ones such as Anderson localization that arises in inhomogeneous systems
[9]. Localized states appear in other classes of experimentally relevant
systems such as nonlinear optics and photonics. Spatial localized patterns
possess potential applications to all-optical control of light, optical
storage, and information processing [6, 8].

Localized patterns involving trivial symmetries such as stripes, hexagons, or
squares have been abundantly discussed and is by now fairly well understood,
including their respective snaking bifurcation diagrams [10, 7]. Indeed, these
localized patterns involve few Fourier modes. However, localized patterns with
non-trivial symmetries have neither experimentally observed nor documented,
nor theoretically predicted. An example of this type of patterning phenomenon
is referred to as localized labyrinthine patterns (LLP). They are observed in
population biology, such as in vegetation populations, on the skin of animals
or human bodies [cf. Fig. 1]. All these examples show an area, which is not
necessarily circular, containing complex spatial structures, a labyrinth, and
surrounded by a uniform state. In the vegetation populations, this intrigued
phenomenon seems to be stationary [see the Supplementary Material [20].
Extended labyrinthine patterns refer to 2D or more dimensional dissipative
structures characterized by a circular or spherical powderlike spectrum
globally [12], and they exhibit a short-range order with a single Fourier mode
[13]. A power spectrum with a powdered ring (sphere) structure is the main
characteristic of patterns with non-trivial symmetries.

Figure 1:  Snapshots of localized labyrinthine patterns in natural systems.
(a) Irregular distribution of vegetation embedded in a uniform vegetated cover
observed in central Cameroon using Google Earth software (with ground
coordinate 3∘58’22.70” N 12∘19’05.84” E). (b) Mamilaria Elongata Cristata with
a contorted tissue (courtesy of David Stang). It is a localized structure in
the bare soil background. (c) Smudgeds composed of stripes and spots in the
skin dorsum of Dendropsophus ozzyi species [11] (reproduced with permission
from the copyright holder). (d) Skin lessions of Tinea imbricata disease
(courtesy of Michael Marks).

In this Letter, we account for the formation of localized patterns with non-
trivial symmetries in well established models from ecology, optics, and
reaction-diffusion systems. We illustrate and investigate this phenomenon
using a Swift-Hohenberg equation (SHE) [14], which constitutes a well-known
paradigm in the study of spatial periodic or localized patterns in spatially
extended systems [2]. We show that this model supports static and stable LLP.
Considering adequate initial conditions, LLP are generated in the coexistence
region between extended labyrinth and homogenous state. We draw the partial
bifurcation diagram showing the stability domain of LLP and their pinning-
depinning transitions, where the localized labyrinth exists as a robust
solution. Free energy allows us to study the relative stability analysis. We
show numerical evidence of stable three-dimensional localized labyrinthine
patterns. Further within the pinning range of parameters, three LLP with
different sizes are generated for a fixed value of the system parameters.

The SHE reads ([14])

$\partial_{t}u=\epsilon u-u^{3}-\nu\nabla^{2}u-\nabla^{4}u,$ (1)

where the real order parameter $u=u(x,y,z,t)$ is an excess scalar field
variable measuring the deviation from criticality, $\epsilon$ is the control
parameter, and $\nu$ the (anti) diffusion coefficient for (positive) negative
value. The cubic term accounts for the nonlinear response of the system under
study. The Laplace operator
$\nabla^{2}=\partial_{xx}+\partial_{yy}+\partial_{zz}$ acts in the
$(x,y,z)$-Euclidean space, and $t$ is time. The last term on the right-hand
side, the bi-Laplacian, stands for hyper-diffusion. Note that Eq. (1) can also
be used to describe two-dimensional systems, where the Laplacian, bi-Laplacian
and the order parameter $u$ are defined in the $(x,y)$-Euclidean space. The
model equation (1) can be rewritten in a variational form as
$\partial_{t}u=-\delta{\cal{F}}/\delta u$, where ${\cal{F}}$ is a Lyapunov
functional or a free energy

${\cal{F}}=\int\frac{dxdydz}{2}\left(-\epsilon
u^{2}+\frac{u^{4}}{2}-\nu(\nabla u)^{2}+(\nabla^{2}u)^{2}\right).$ (2)

The variational structure of the SHE (1) indicates that only stationary
solutions such as uniform states, spatially periodic, or localized patterns
are possible. The SHE (1) is a well-known paradigm for the study of periodic
and localized patterns, and has been first derived in hydrodynamics [14], and
later in other fields of natural science, such as chemistry [15], and
nonlinear optics [16]. In the last decade, it has been established that Eq.
(1) had already been constructed by Alan Turing, but was unpublished [17].

Other real SHE has been derived for out of equilibrium systems [18, 19]

$\partial_{t}u=\eta-\epsilon u-u^{3}-(\nu-bu)\nabla^{2}u-\nabla^{4}u-c(\nabla
u)^{2},$ (3)

where $\epsilon$ and $\eta$ are control parameters, $\nu$, $b$ are diffusion
parameters, and $c$ is the nonlinear advection strength. The presence of
nonlinear diffusion and nonlinear advection terms, $u\nabla^{2}u$ and $(\nabla
u)^{2}$, render Eq. (3) non-variational. In general, this equation does not
admit a Lyapunov functional.

Figure 2: Stationary localized labyrinthine patterns obtained in different
pattern forming models: (a) SHE (1) ($\epsilon=1.17$, $\nu=1$). (b)
generalized SHE (3) ($\epsilon=0.2$, $\nu=1$, $\eta=-0.06$, $b=0.1$, $c=0.1$).
(c) Non-local vegetation. (d) Passive diffractive resonator. (e) reaction-
diffusion. See the Supplementary Material for more details of models and
parameters used in (c), (d), and (e) [20]. The right-upper inset in (a) shows
the powderlike ring specturm of the LLP in the SHE model (1). All the
localized structures fulfill the definition of labyrinthine patterns (see
Supplementary Material [20] for details).

Whether a SHE model is variational or not, numerical simulations of both
models, Eqs. (1,3), with periodic boundary conditions show evidence of stable
stationary localized labyrinthine patterns [see Figs. 2(a) and 2(b)].

To obtain localized labyrinthine patterns, the initial conditions consist of a
circular area of a stable labyrinthine pattern in the center of the simulation
box, embedded in a uniform background. The evolution towards equilibrium
starts with a quick adjustment of the interface mediated by the curvature of
the stripe patterns; then, there is an accommodation of the stripe patterns in
the bulk. To the end, some retraction of stripes in the interface takes place
(cf. Video 1 and the stabilization of LLP in the Supplementary Material [20]).
The localized region is not perfectly circular, containing finite segments of
deformed stripes separated by spots of the same width. These finite-size
stripes can be interconnected or not. They support all stripe orientations
along the motionless interface separating the labyrinth to the homogeneous
steady state, as shown in Fig. 2. Besides, the formation of LLP in the above
scalar model equations in the form of SHE, additional models are also
considered that are experimentally relevant: (i) a generic interaction
redistribution model describing vegetation pattern formation which is an
integrodifferential model equation. This simple modeling approach based on the
interplay between short-range and long-range interactions governing plant
communities captures localized labyrinthine pattern as shown in Fig. 2(c).
(ii) broad area photonics devices such as nonlinear resonators subjected to a
coherent injected beam [see Fig. 2(d)]. In this case, the resulting equation
is a complex Ginzburg-Landau type equation. Finally, a reaction-diffusion
model for chemical dynamics, also supports LLP as shown in Fig. 2(e). The
description of these models and the values of the parameters are provided in
the Supplementary Materials [20]. Similar solutions when using Dirichlet and
Neumann boundary conditions are observed (see the Supplementary Material [20]
for details).

Figure 3: Three dimensional localized labyrinthine pattern solution of Eq.
(3). (a) Colormap of the full simulation, (b) and (c) are colormap slices of
the localized labyrinthine pattern. (d) Isosurface of the localized
labyrinthine pattern with $u=1.3$. Parameters are $\epsilon=4.2$, $\nu=5$,
$\eta=-6.8$, $b=0$, and $c=0$. The mesh integration is $40\times 40\times 40$.

The 2D LLP are robust structures in 2D systems in the various natural system,
as shown in Fig. 1 and in the Supplementary Material [20]. It has been shown
that the Swift-Hohenberg equation supports 3D extended patterns with trivial
symmetries such as lamellae, body-centered cubic crystals, hexagonally packet
cylinders [21, 22, 23, 24], and localized patterns [23, 24, 25]. Recently,
clusters of three-dimensional bullets forming a localized crystal with trivial
symmetry have been reported [26]. We extend this analysis to 3D non-trivial
symmetry patterns and we show the existence of stable 3D localized
labyrinthine patterns. They consist of finite size curved and connected tubes
embedded in a homogeneous background. The width of the tubes is half of the
critical wavelength at the symmetry-breaking instability. They are obtained by
numerical simulations of the generalized SHE Eq. (3) with Neumann boundary
conditions along $x$, $y$, and $z$ directions. Figure 3 shows a typical 3D
localized labyrinthine pattern.

The homogeneous steady states $u_{s}=0$ and $u_{s\pm}=\pm\epsilon^{1/2}$
solutions of Eq. (1) undergo symmetry-breaking instabilities at
$\epsilon_{c1}=-\nu^{2}/4$ and $\epsilon_{c2}=\nu^{2}/8$. At both critical
bifurcations points the critical wavelength is
$\lambda_{c}=2\pi/k_{c}=2\sqrt{2}\pi/\sqrt{\nu}$. Indeed, when the linear
coefficient of the Laplacian is negative; $\nu>0$, the spontaneous pattern
formation process becomes possible thanks to the appearance of a finite band
of linearly unstable Fourier modes that triggers the appearance of spatially
periodic patterns. The upper cutoff is affected by the bi-Laplacian term,
which is always stabilizing for short distances, since dispersion is an
efficient mixing mechanism. Numerical simulations of the bidimensional Eq. (1)
in the neighborhood of the critical point $\epsilon=\epsilon_{c1}$ indicate
the emergence of extended patterns, as shown in Fig. 4(a). When increasing the
control parameter, the sequence of symmetry-breaking transitions fingerprint,
glassy, and scurfy labyrinthine patterns are observed [13].

Figure 4:  Bifurcation diagram of homogenous solutions (a) and relative
stability analysis of LLP (b) in SHE (1), with $\nu=1$. The uniform $u(x,y)=0$
state suffers a Turing instability at $\epsilon_{c1}=-0.25$. The uppermost
curve (blue) shows the maximum $u(x,y)$ of the different equilibrium patterns
[stripe (I), fingerprint-type labyrinth (II), glassy labyrinth (III), and
scurfy labyrinth (IV)]. At $\epsilon_{c2}=0.125$, the non-zero homogenous
states (HSS-green curves) are unstable to in-homogenous perturbations. In the
narrow region limited by $\epsilon_{p}^{-}=1.16$ and $\epsilon_{p}^{+}=1.19$,
where the extended labyrinthine pattern and the uniform solutions coexist, the
existence of localized labyrinthine patterns (V) is possible. Free energy
${\cal{F}}$ given by Eq. (2) is computed for the extended labyrinths (EL), the
stripe pattern (SP), the homogenous states (HSS), and the localized
labyrinthine patterns (LLP) near the pinning region.

Spatial confinement leading to the formation of localized patterns with
nontrivial symmetry occurs in parameter space ($\epsilon>\epsilon_{c2}$),
where extended labyrinthine patterns coexist with a homogeneous steady state.
Within this hysteresis loop, there generally exists a so-called pinning range
of parameters [27], delimited by $\epsilon^{\pm}_{p}$, in which robust LLP can
be observed [cf. Fig. 4(a)]. Taking advantage of the variational structure of
Eq. (1), we address the problem of the relative stability analysis. We
evaluate numerically $\mathcal{F}$, associated with uniform states, extended
labyrinth, and LLP. Figure 4(b) summarizes the results. These equilibria
correspond to a local or global minimum of Lyapunov functional $\mathcal{F}$
given by Eq. (2). From Fig. 4(b), we see that the LLP are more robust than the
extended labyrinth but less stable than the homogeneous steady state. The
localized labyrinth solutions are stable only in the pinning range of
parameters $\epsilon^{-}_{p}<\epsilon<\epsilon^{+}_{p}$. The localized
labyrinths are stationary and robust patterns since their localized area never
expands despite diffusion and never shrinks despite nonlinearity and
dissipation.

Figure 5:  Three stable branches of LLP and pinning-depinning transitions in
the SHE model (1). (a) Plot of $||u||^{2}$ for three initial conditions (i,
ii, and iii) with different sizes (see Fig. 3 in the Supplementary Material
[20] for details). The upward (downward) triangles account for the increasing
(decreasing) of $\epsilon$, starting from $\epsilon_{ic}=1.17$. (b) Temporal
sequence of the pinning-depinning transition when crossing
$\epsilon_{p}^{-}$=1.16, and (c) $\epsilon_{p}^{+}=1.19$.

Localized patterns with trivial symmetry (stripes and hexagons) are organized
into a complex diagram following a homoclinic snaking bifurcation [7]. This
type of diagram is obtained by a continuation method. The symmetries of the
localized structures are relevant for the convergence of this algorithm.
However, in the case of LLP there is a lack of continuation algorithms to
characterize the full bifurcation diagram. To figure out the existence region
of stable LLP, we performed direct numerical simulations of Eq. (1). Figure
5(a) summarizes the results, where we plot

$||u||^{2}=\dfrac{1}{L_{x}L_{y}}\int_{0}^{L_{x}}\int_{0}^{L_{y}}(u(x,y)-u_{s+})^{2}dxdy,$
(4)

as a function of the bifurcation parameter. The full bifurcation diagram can
be complex, so we display only three branches of LLP obtained with different
initial conditions shown in the insets (i), (ii), and (iii) of Fig. 5(a). The
maximum amplitude of the three LLP is the same, but they have different size.
Varying $\epsilon$ from these initial conditions, we obtain the three branches
shown in Fig. 5(a). Whatever the initial condition, when increasing the
bifurcation parameter the LLP decrease in size, mediated by shrinking of
fingers and accommodation of defects [see the Supplementary Material [20] for
details]. All LLP disappear close to $\epsilon>\epsilon_{p}^{+}$ and the
system exhibit a transition towards a mixture of circular localized peaks and
dips. Figure 5(b) illustrates this transition, during which we observe the
contraction of fingers, which transform to circular peaks or dips. This
process correspond to the inverse of the invagination of localized structures
[28]. Starting from the initial conditions shown in the insets (i), (ii), and
(iii) of Fig. 5(a) and decreasing the bifurcation parameter, we observe an
increase in the size of the LLP. Further decreasing
$\epsilon<\epsilon_{p}^{-}$, we observe a transition to an extended
fingerprint-like labyrinthine pattern. This depinning transition mediated by
front propagation has the tendency to reduce the number of circular spots and
dips, and enhance the invagination process as illustrated in Fig. 5(c).

By varying the control parameter within the pinning region delimited by
$\epsilon_{p}^{-}$ and $\epsilon_{p}^{+}$, we see that for a fixed $\epsilon$,
and near $\epsilon=\epsilon_{p}^{-}$ the sizes of the coexisting LLP are
different. However, close to $\epsilon=\epsilon_{p}^{+}$ the system reaches
more or less the same size. We stress that the position of LLP and their size
depends on the initial conditions, and the maximum of the coexisting LLP is
essentially constant for fixed values of the system parameters. The number of
coexisting LLP with different sizes can be much larger than the three branches
shown in the bifurcation diagram displayed in Fig. 5(a).

The localized patterns with trivial symmetries have a well established
bifurcation diagram based on continuation methods. However, when dealing with
localized patterns with non trivial symmetries, there are no available
algorithms for the continuation to handle this problem. Whether the localized
labyrinthine patterns present a homoclinic snaking bifurcation diagram or not
remains an open question. The plausibility of spatial varying parameters can
be responsible for complex localized patterns. However, our result opens a
novel possibility of localized patterns with non-trivial symmetries even in
homogenous and isotropic systems.

###### Acknowledgements.

MGC thanks for the financial support of ANID–Millennium Science Initiative
Program–ICN17_012 (MIRO) and FONDECYT projects 1210353. S.E.-A. acknowledges
the financial support of ANID by Beca Doctorado Nacional 2020-21201376. A part
of this work was supported by the ”Laboratoire Associé International”
University of Lille - ULB on ”Self-organisation of light and extreme events”
(LAI-ALLURE). M.T acknowledges financial support from the Fonds de la
Recherche Scientifique FNRS under Grant CDR no. 35333527 ”Semiconductor
optical comb generator”.

## References

  * [1] P. Glansdorff and I. Prigogine, Thermodynamic theory of structure, stability and fluctuations (Wiley-Interscience, 1971).
  * [2] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * [3] J. D. Murray, Mathematical biology: I. An introduction, (Springer Science & Business Media, 2007).
  * [4] O. Descalzi, M. Clerc, S. Residori, and G. Assanto, Localized States in Physics: Solitons and Patterns (Springer, New York, 2011).
  * [5] H. G. Purwins, H. U. Bödeker, and Sh. Amiranashvili, Dissipative solitons, Adv. Phys. 59, 485 (2010).
  * [6] M. Tlidi, K. Staliunas, K. Panajotov, A.G. Vladimirov, and M.G. Clerc, Localized structures in dissipative media: from optics to plant ecology, Phil. Trans. R. Soc. A 372, 20140101 (2014).
  * [7] E. Knobloch, Spatial Localization in Dissipative Systems, Annu. Rev. Condens. Matter Phys. 6, 325 (2015).
  * [8] L. Lugiato, F. Prati, and M. Brambilla, Nonlinear Optical Systems (Cambridge University Press, 2015).
  * [9] P. W. Anderson, Absence of Diffusion in Certain Random Lattices, Phys. Rev. 109, 492 (1958).
  * [10] P. Coullet, Localized patterns and fronts in nonequilibrium systems, Int. J. Bifurcation Chaos Appl. Sci. Eng. 12, 2445 (2002).
  * [11] V. G. D. Orrico, A new ”Bat-Voiced” species Dendrophosus Fitzinger, 1843 (Anura, Hyliade) from the Amazon Basin, Brazil, Zootaxa 3881, 4 (2014).
  * [12] M. Le Berre, E. Ressayre, A. Tallet, Y. Pomeau, and L. Di Menza, Example of a chaotic crystal: The labyrinth, Phys. Rev. E 66, 026203 (2002).
  * [13] S. Echeverria-Alar and M. G. Clerc, Labyrinthine patterns transitions, Phys. Rev. Research 2, 042036(R) (2020).
  * [14] J. Swift and P.C. Hohenberg, Hydrodynamic fluctuations at the convective instability, Phys. Rev. A 15, 319 (1977).
  * [15] M’F Hilali, G. Dewel, and P. Borckmans, Subharmonic and strong resonances through coupling with a zero mode, Phys. Lett. A 217, 263 (1996).
  * [16] M. Tlidi, P. Mandel, and R. Lefever, Localized structures and localized patterns in optical bistability, Phys. Rev. Lett. 73, 640 (1994).
  * [17] J. H. Dawes, After 1952: The later development of Alan Turing’s ideas on the mathematics of pattern formation, Hist. Math. 43, 49 (2016).
  * [18] M. G. Clerc, A. Petrossian, and S. Residori, Bouncing localized structures in a liquid-crystal light-valve experiment, Phys. Rev. E 71, 015205(R) (2005).
  * [19] G. Kozyreff and M. Tlidi, Nonvariational real Swift-Hohenberg equation for biological, chemical, and optical systems, Chaos 17, 037103 (2007).
  * [20] See Supplementary Material for details on the complexity of localized labyrinthine patterns in the models described here (Eqs. 1 and 3) and in [20] (vegetation, nonlinear optics, reaction-diffusion). Video 1 shows the stabilization process of a localized labyrinthine pattern in the SHE model 1. Details on numerical simulations of LLP in the SHE model 1 with Neumann and Dirichlet boundary conditions. The growth mechanism along one stable branch of LLP is described.
  * [21] U. Thiele, A.J. Archer, M.J. Robbins, H. Gomez, and E. Knobloch, Localized states in the conserved Swift-Hohenberg equation with cubic nonlinearity, Phys. Rev. E, 87, 042915 (2013).
  * [22] H. G. Lee, A semi-analytical Fourier spectral method for the Swift–Hohenberg equation. Computers & Mathematics with Applications 74, 1885 (2017).
  * [23] K. Staliunas, Three-dimensional Turing structures and spatial solitons in optical parametric oscillators, Phys. Rev. Lett. 81, 81 (1998).
  * [24] M. Tlidi and P. Mandel, Three-dimensional optical crystals and localized structures in cavity second harmonic generation, Phys. Rev. Lett. 83, 4995 (1999).
  * [25] I. Bordeu and M. G. Clerc, Rodlike localized structure in isotropic pattern-forming systems, Phys. Rev. E 92, 042915 (2015).
  * [26] S. S. Gopalakrishnan, K. Panajotov, M. Taki, and M. Tlidi, Dissipative Light Bullets in Kerr Cavities: Multistability, Clustering, and Rogue Waves, Phys. Rev. Lett. 126, 153902 (2021)
  * [27] Y. Pomeau, Front motion, metastability and subcritical bifurcations in hydrodynamics, Physica (Amsterdam) 23D, 3 (1986).
  * [28] I. Bordeu, M. G. Clerc, R. Lefever, and M. Tlidi. From localized spots to the formation of invaginated labyrinthine structures in a Swift–Hohenberg model, Comm. in Nonlinear Sci. and Num. Sim. 29, 482 (2015).