# Self-organized vortex and antivortex patterns in laser arrays

M. Honari-Latifpour  J. Ding  S. Takei  M.-A. Miri<EMAIL_ADDRESS>Department of Physics, Queens College of the City University of New York,
Queens, New York 11367, USA
Physics Program, The Graduate Center, City University of New York, New York,
New York 10016, USA

###### Abstract

Recently it is shown that dissipatively coupled laser arrays simulate the
classical XY model. We show that phase-locking of laser arrays can give rise
to the spontaneous formation of vortex and antivortex phase patterns that are
analogous to topological defects of the XY model. These patterns are stable
although their formation is less likely in comparison to the ground state
lasing mode. In addition, we show that small ratios of photon to gain lifetime
destabilizes vortex and antivortex phase patterns. These findings are
important for studying topological effects in optics as well as for designing
laser array devices.

††preprint: AIP/123-QED

##  I Introduction

The two-dimensional XY model consists of a lattice of interacting fixed-length
spins that are constrained to rotate in a plane Chaikin, Lubensky, and Witten
(1995). The classical XY model in two spatial dimensions is governed by the
Hamiltonian $\mathcal{H}(\phi_{1},\phi_{2},\cdots,\phi_{n})$:

$\mathcal{H}=\sum_{\langle i,j\rangle}\kappa_{ij}\cos(\phi_{i}-\phi_{j})$ (1)

where, $\phi_{1},\cdots,\phi_{n}$ represent the orientation of the $n$ spins,
and $\kappa_{ij}$ is the interaction for the pair $i$,$j$. This model supports
nontrivial equilibrium spin configurations — a class of topological defects —
known as vortices, which are characterized by the phase of the spins going
through a multiple of $2\pi$ as one traces a loop enclosing the vortex, e.g.,

$\int d\phi=\pm 2\pi$ (2)

for a single vortex/antivortex. In the 2D XY model, vortices have found
applications in various areas of condensed matter physics including superfluid
helium-4 Bishop and Reppy (1980), superconductivity in thin films Beasley,
Mooij, and Orlando (1979); Hebard and Fiory (1980), liquid crystals Pargellis,
Green, and Yurke (1994), and the melting of 2D crystals Halperin and Nelson
(1978); Nelson and Halperin (1979).

Recently, it has been realized that the classical XY model can be optically
simulated with an array of coupled optical oscillators Nixon  _et al._ (2013);
Berloff  _et al._ (2017); Parto  _et al._ (2020); Honari-Latifpour and Miri
(2020). What makes this possible is the random phase of a laser above
oscillation thresholds, which emulates a classical spin confined to a two-
dimensional plane. In addition, dissipative interaction facilitates
synchronization of an array of lasers to a globally-phase-locked state Ding
_et al._ (2019). In this case, one can show that the laser array reaches an
equilibrium phase pattern that locally minimizes a cost function that is
asymptotically equivalent with the classical XY Hamiltonian Honari-Latifpour
and Miri (2020). Accordingly, laser arrays have been utilized for simulating
interesting phenomena related to spin systems such as geometric frustration
Nixon  _et al._ (2013).

In this work, we investigate the formation of vortex and antivortex
singularities as self-organized phase patterns in laser arrays. As shown
schematically in Fig. LABEL:fig1, these patterns are equilibrium states of
laser arrays when reaching a globally-phase-locked state. From a nonlinear
dynamics point of view, these topological defects are fixed-point solutions of
nonlinear dynamical equations governing laser arrays. However, in general,
they can be considered as metastable states with finite basins of attraction
which restricts their formation to proper initial conditions. In addition, we
find that the stability of these states depends critically on the gain level
and lifetime. In the following, after introducing a dynamical model governing
laser arrays, first we numerically investigate the formation of vortices.
Next, we draw a connection between the governing dynamical model and a class
of Ginzburg-Landau systems that are known to support vortex patterns. Finally,
we investigate the stability of the vortex patterns with respect to the
competing time scales of optical cavity and gain decay rates of the lasers.

##  II Model

Figure 2: Equilibrium phase patterns of an array of dissipatively coupled
lasers arranged on a $16\times 16$ square lattice. (a-d) The ground state,
vortex, antivortex and bound vortex-antivortex states for the case of
attractive coupling ($\kappa>0$), associated with a ferromagnetic spin system
[first row], and repulsive coupling ($\kappa<0$), associated with an
antiferromagnetic spin system [second row]. The XY energy level associated
with these phase patterns are shown on top of the panels. Here, $g_{0}=30$,
$\kappa=-1$ for the top row and $\kappa=1$ for the bottom row.

To build a dynamical model governing laser arrays, we first consider an array
of passive and single-mode optical resonators that dissipatively interact Ding
and Miri (2019); Ding  _et al._ (2019). For the sake of simplicity, we assume
all resonators being identical in resonance frequency $\omega_{0}$ and
linewidth $1/\tau_{p}$. Thus, in the framework of the temporal coupled mode
theory Haus (1984), the complex modal amplitude of the electric field in the
$i$th resonator is governed by:

$\dot{a}_{i}(t)=-a_{i}-\gamma_{i}a_{i}-\sum_{j\neq i}\kappa_{ij}a_{j},$ (3)

where, the equations are written within the gauge $a_{i}\rightarrow
a_{i}e^{-i\omega_{0}t}$ and time is normalized to the photon lifetime
$\tau_{p}$. In this relation, $\kappa_{ij}$ represents the rate of dissipative
coupling between the $i$th and $j$th resonators,
$\gamma_{i}=\sum_{j}|\kappa_{ij}|$ is the external loss of the $i$th resonator
due to its coupling with other resonators as demanded by conservation
relations Suh, Wang, and Fan (2004). Here, all coupling coefficients
$\kappa_{ij}$ are normalized to the photon decay rate $1/\tau_{p}$, while the
choice of normalization for the complex amplitudes depends on the gain as
discussed next.

By incorporating a saturable gain mechanism, equations (3) can be modified to
support self-sustained finite stationary solutions of the field amplitudes.
Here, we consider the gain being a dynamical variable as in the so-called
class-B laser model Tredicce  _et al._ (1985). In this model, the gain of a
laser oscillator is driven at a finite pump rate, while it decays linearly for
small field intensities and nonlinearly when the field intensity grows. The
normalized rate equations for the $i$th oscillator can then be written as:

$\displaystyle\dot{a}_{i}(t)$
$\displaystyle=[g_{i}(t)-1]a_{i}-\gamma_{i}a_{i}-\sum_{j\neq
i}\kappa_{ij}a_{j},$ (4a) $\displaystyle\dot{g}_{i}(t)$
$\displaystyle=(\tau_{p}/\tau_{g})[g_{0}-(1+|a_{i}|^{2})g_{i}],$ (4b)

where $g_{i}$ represents the gain of the $i$th oscillator, $g_{0}$ is the pump
parameter, and $1/\tau_{g}$ is the gain decay rate. In these relations both
the field amplitude and gain are dimensionless and the time is normalized to
the photon lifetime $\tau_{p}$. This model has been applied to solid-state
lasers Fabiny  _et al._ (1993). In addition, it can be generalized to model
semiconductor lasers by incorporating the linewidth enhancement factor which
plays an important role in the dynamics Ohtsubo (2012).

Figure 4: The steady state phase pattern of a laser array arranged on a
triangular lattice. Here, $\kappa=-1$ and $g_{0}=30$.

Equations (4) can be greatly simplified when the gain decay rate is much
larger than the photon decay rate, i.e., $1/\tau_{g}\gg 1/\tau_{p}$. In this
case, the gain almost instantaneously follows the dynamics of the field. Thus,
one can adiabatically eliminate the dynamics of the gain, i.e.,
$\dot{g}_{i}(t)\approx 0$, to reach at an instantaneous nonlinear gain
$g_{i}(|a_{i}|)=g_{0}/(1+|a_{i}|^{2})$. In this manner, one reaches at a
reduced model that is suitable for a so-called class-A laser Tredicce  _et
al._ (1985); Mandel and Wolf (1995). Here, we use a polynomial gain term
$g_{i}(a_{i})=g_{0}(1-|a_{i}|^{2})$, which is enough to guarantee the bounded
oscillations of the laser. Thus, we reach at the following reduced model,

$\dot{a}_{i}(t)=[g_{0}(1-|a_{i}|^{2})-1]a_{i}-\gamma_{i}a_{i}-\sum_{j\neq
i}\kappa_{ij}a_{j},$ (5)

which is accurate as long as the photon decay rate is smaller than the decay
rates of the atomic degrees of freedom that give result to the gain. In this
work, first we focus on simulating equations (5). Next, we return to equations
(4) and discuss the effect of the gain dynamics.

The dynamical model of Eq. (5) admits a Lyapunov function $F$ such that
$\dot{a}_{i}=-\partial{F}/\partial{a_{i}^{*}}$, where Honari-Latifpour and
Miri (2020)

$F=\sum_{i}[-(g_{0}-1-\gamma_{i})|a_{i}|^{2}+\frac{g_{0}}{2}|a_{i}|^{4}]+\frac{1}{2}\sum_{i,j}a_{i}^{*}\kappa_{ij}a_{j}.$
(6)

The fixed point solutions of the dynamical system of Eq. (5) are local minima
of this Lyapunov function. In a recent work we showed that the diagonal term
in this cost function behaves like a penalty term that tends to force all
oscillators to a constant amplitude in the large gain limit, $g_{0}\gg 1$
Honari-Latifpour and Miri (2020). Thus, considering the stationary state
solution of the oscillators as $a_{i}=\sqrt{I_{i}}e^{i\phi_{i}}$ by enforcing
the condition of $I_{i}=I_{0}$, the cost function of Eq. (6) reduces to the XY
Hamiltonian of Eq. (1).

##  III Formation of Topological Defects

In the following, we investigate self-organization of topological defects by
numerically simulating the dynamical model of Eqs. (5). We first consider an
array of lasers arranged on a square lattice with uniform nearest-neighbor
coupling of strength $\kappa$. Here, the gain is assumed to be large
($g_{0}\gg 1$), so that steady-state amplitudes become nearly uniform and the
phase pattern obeys the XY Hamiltonian. Figure 2 depicts the equilibrium phase
patterns obtained by simulating a square lattice arrangement of $16\times 16$
oscillators for both scenarios of attractive ($\kappa<0$) and repulsive
($\kappa>0$) coupling, which are respectively associated with ferromagnetic
and antiferromagnetic cases. In both cases, different stable patterns are
observed, including the ground state, isolated vortex and antivortex states,
and paired vortex-antivortex states. The XY energy levels associated with
these equilibrium phase patterns are listed in Fig. 2, which shows higher
energy levels for the topological defects. The difference between the XY
energy of the vortex (Fig. 2(b)) and the ground state (Fig. 2(a)) is
comparable with the approximate formula $\Delta E=\pi\kappa\ln{L}$, where $L$
is the lattice length.

The stability of these fixed point solutions is directly evaluated through the
Jacobian matrix of the dynamical system of Eqs. (5):

$\mathrm{J}=g_{0}\begin{pmatrix}\mathrm{diag}(1-2\bar{\mathbf{a}}\odot\bar{\mathbf{a}}^{*})-\mathrm{I}-\mathrm{Q}&-\mathrm{diag}(\bar{\mathbf{a}}\odot\bar{\mathbf{a}})\\\
-\mathrm{diag}(\bar{\mathbf{a}}^{*}\odot\bar{\mathbf{a}}^{*})&\mathrm{diag}(1-2\bar{\mathbf{a}}\odot\bar{\mathbf{a}}^{*})-\mathrm{I}-\mathrm{Q}\end{pmatrix}.$
(7)

In this relation, $\bar{\mathbf{a}}=(\bar{a}_{1},\cdots,\bar{a}_{n})^{t}$ is
the stationary state, $\mathrm{I}$ represents the identity matrix,
$\mathrm{Q}$ is the coupling matrix, where $q_{ij}=\kappa_{ij}$ and
$q_{ii}=\sum_{j}|\kappa_{ij}|$, $\odot$ represents the entry-wise product, and
$\mathrm{diag}$ creates a diagonal matrix of a given vector. The Jacobian
matrix turns out to be a negative semi-definite matrix for all cases shown in
Fig. 2.

Given the finite attraction basins of the phase patterns shown in Fig. 2,
formation of these patterns depends on the initial conditions. In a realistic
scenario, one can assume that the initial conditions are uniformly sampled
from the high-dimensional phase space of the laser array. Thus, as local
minima of the Lyapunov function, topological defects are expected to have
lower chances of formation. To explore this aspect, we simulated the dynamical
system for a large ensemble of random initial conditions. In these
simulations, the initial phases were drawn from a uniform random distribution
within the range $[0,2\pi]$, while the initial amplitudes were taken to be
small values. The resulting equilibrium energy distribution is shown in Fig.
LABEL:fig3 for both cases of ferromagnetic and antiferromagnetic systems. As
this figure clearly indicates, topological defects form at much lower
probabilities compared to the ground state.

The vortex and antivortex phase patterns can also form in the triangular
lattice geometry as shown in the simulations of Fig. 4. In this lattice
geometry, the antiferromagnetic case is more complex due to the geometric
frustration. The XY energy distributions associated with the equilibrium phase
patterns of the ferromagnetic and antiferromagnetic cases are depicted in Fig.
LABEL:fig5. This figure indicates the complexity of the antiferromagnetic
system due to its large number of higher energy states.

It is worth noting that the transient dynamics of the laser array reveals
unstable topological defects that disappear by reaching the lattice boundaries
or by collapsing of vortex-antivortex pairs. An exemplary dynamics of a laser
network is shown in Fig. LABEL:fig6. To visualize the transient dynamics, the
phase pattern in depicted at intermediate time scales. It is worth noting that
similar behavior was reported in a reduced model based on Kuramoto phase
oscillators, which showed rapid decay of a large number of transient
topological defects Mahler  _et al._ (2019).

##  IV Discrete Ginzburg Landau Equation

In order to provide insight to the formation vortices in laser arrays, we draw
a connection between the lattice model of Eq. (5) with its continuum
counterpart, which turns out to be the well-known Ginzburg-Landau Equation
(GLE). To understand this analogy, we focus our attention to the case of a 2D
square lattice arrangement of lasers with uniform nearest neighbor coupling
$\kappa$, and consider the ferromagnetic case ($\kappa<0$). By considering a
pair of integer indices $(m,n)$ for describing the horizontal and vertical
coordinates of a square lattice, one can rewrite Eqs. (5) as:

$\dot{a}_{m,n}(t)=[g_{0}(1-|a_{m,n}|^{2})-1]a_{m,n}-\kappa\mathrm{L}a_{m,n}$
(8)

where, $\mathrm{L}$ is the discrete Laplacian operator on a square lattice
graph that acts as
$\mathrm{L}a_{m,n}=a_{m-1,n}+a_{m+1,n}+a_{m,n-1}+a_{m,n+1}-4a_{m,n}$. This
relation is clearly in a finite difference form. To construct the continuum
counterpart of this relation, we use $(m,n)\rightarrow(x,y)$,
$a_{m,n}(t)\rightarrow\psi(x,y,t)$, and $\mathrm{L}\rightarrow\nabla^{2}$,
which results in the Ginzburg-Landau Equation (GLE):

$\dot{\psi}(x,y,t)=[g_{0}(1-|\psi|^{2})-1]\psi-\kappa\nabla^{2}\psi.$ (9)

Given that all coefficients are real-valued, despite the fact that $\psi$ is
complex, this equation is often referred to as the real Ginzburg-Landau
equation Aranson and Kramer (2002).

The dynamical equation (9) can be written in terms of a variation of a
functional $\mathcal{F}[\psi]$, as
$\dot{\psi}=-\delta\mathcal{F}/\delta\psi^{*}$. The functional $\mathcal{F}$
is found to be:

$\mathcal{F}=\int{dxdy\left[-(g_{0}-1)|\psi|^{2}+\frac{g_{0}}{2}|\psi|^{4}+|\nabla\psi|^{2}\right]},$
(10)

which, can be considered as a counterpart of the Lyapunov function of Eq. (6).
In an amplitude and phase representation
$\psi(r,\phi,t)=R(r,\phi,t)\exp[i\Phi(\phi(r,\phi,t))]$, Eq. (9) can be
rewritten as:

$\displaystyle\dot{R}(r,\phi,t)$
$\displaystyle=[g_{0}(1-R^{2})-1]R-\kappa(\nabla^{2}-|\nabla\Phi|^{2})R,$
(11a) $\displaystyle R\dot{\Phi}(r,\phi,t)$ $\displaystyle=-\kappa(2\nabla
R\cdot\nabla\Phi+R\nabla^{2}\Phi).$ (11b)

These equations are a special class of a reaction-diffusion system, named
$\lambda-\omega$ systems Aranson and Kramer (2002). It is shown that these
equations support stable single-arm spiral wave solutions of the form
$R(r,\phi,t)=R(r)$ and $\Phi(r,\phi,t)=\phi$, while multi-arm spiral waves are
unstable Hagan (1982). Driven from this intuition, one would expect discrete
counterparts of the spiral patterns in lattices of coupled lasers. An obvious
deviation of the lattice model from the continuum case is the finite amplitude
of the field at the center of the vortex. In the continuum scenario, the
vanishing amplitude at the phase singularity guarantees a well-defined field.
In the case of the lattice model, however, there is no such requirement. This
is clearly irrelevant in the case of the lattice model, although the
amplitudes still tend to be lower at the lattice sites that are located near
the vortex center. Nevertheless, the amplitudes become uniform in the large
gain limit ($g_{0}\rightarrow\infty$), where the laser system approaches the
XY model.

##  V Vortices in the Dynamic Gain Regime

The results presented in Figs. 2-LABEL:fig6 are based on the complex amplitude
model with an instantaneous nonlinear saturable gain described by Eqs. (5).
Here, we consider the more realistic scenario that involves gain as a
dynamical variable as described through Eqs. (4). The gain lifetime
$\tau_{g}$, which in reality is dictated by the laser gain material, is found
to play a critical role in the dynamics and in formation of stationary
topological defects. As mentioned earlier, when the ratio of the photon
lifetime over the gain lifetime is large ($\tau_{p}/\tau_{g}\gg 1$), the gain
dynamics can be effectively disregarded. Thus, in this regime, stable vortex
and antivortex patterns are expected. It remains to investigate the stability
of these topological defects when the photon to gain lifetime ratio becomes
small ($\tau_{p}/\tau_{g}\sim 1$ or $\tau_{p}/\tau_{g}\ll 1$).

To explore this aspect, we first simulated the dynamics of a small array of
lasers for different values of $\tau_{p}/\tau_{g}$ while using identical
initial conditions. The results are depicted in Figs. LABEL:fig7. As these
exemplary simulations indicate, by decreasing the ratio of
$\tau_{p}/\tau_{g}$, the transient dynamics becomes more complicated,
involving effects like self-pulsations that are known processes in class-B
lasers. On the other hand, these figures show that the vortex state does not
form for smaller ratios of $\tau_{p}/\tau_{g}$ (Figs. LABEL:fig7(b,c)).

This can be explained through the dynamics of the gain. According to Figs.
LABEL:fig7(b,c), by increasing $\tau_{g}$, the gain increases slowly and
during this process, the laser array finds the time to escape from patterns
associated with higher XY energy levels and to stabilize into the ground state
phase pattern. This is in agreement with our recent study based on the reduced
model, while the gain was adiabatically increased to avoid trapping into local
minima Honari-Latifpour and Miri (2020). Here, the dynamics allows the gain to
automatically increase over a time scale governed by $\tau_{g}$. To
systematically explore this aspect, the laser array was simulated with an
ensemble of random initial conditions and for three different time scale
ratios. In addition, both cases of rectangular and triangular lattices were
considered. The XY energy distributions of the associated steady-state phase
patterns are shown in Fig. 8. As expected, for $\tau_{p}/\tau_{g}\gg 1$, the
equilibrium energy distribution is similar to the case of instantaneous gain
(Figs. LABEL:fig3(b) and LABEL:fig5(b)). As this ratio increases, however, the
occurrence of higher energy patterns drops significantly.

Figure 8: The XY energy distribution associated with the equilibrium phase
pattern of laser arrays governed by the dynamical equations (4) for different
photon to gain lifetime ratios. Top: The square lattice for (a)
$\tau_{p}/\tau_{g}=10$, (b) $1.0$, and (c) $0.01$. Bottom: The triangular
lattice for (a) $\tau_{p}/\tau_{g}=10$, (b) $1.0$, and (c) $0.01$. Here,
$\kappa=-1$ and $g_{0}=30$, while the lattice size is the same as in Figs. 2
and 4 for the top and bottom rows, respectively.

These results clearly suggest that smaller ratios of $\tau_{p}/\tau_{g}$
destabilize the vortex and antivortex patterns. However, it is quite
interesting to note that these patterns pass the linear stability test in the
both regimes of small and large time scale ratios. To show this, first we
consider the fixed point solutions of Eqs. (4), i.e.,
$(a_{i}(t),g_{i}(t))=(\bar{a}_{i},\bar{g}_{i})$. The fixed points are governed
by a set of nonlinear algebraic equations:

$\displaystyle[\bar{g}_{i}-1]\bar{a}_{i}-\gamma_{i}\bar{a}_{i}-\sum_{j\neq
i}\kappa_{ij}\bar{a}_{j}=0,$ (12a)
$\displaystyle[g_{0}-(1+|\bar{a}_{i}|^{2})\bar{g}_{i}]=0,$ (12b)

which, are independent from the time scale ratio $\tau_{p}/\tau_{g}$. Thus, as
an example, the vortex state found in Fig. LABEL:fig7(a) for
$\tau_{p}/\tau_{g}=10$ is also a valid solution of the system explored in Fig.
LABEL:fig7(c) for $\tau_{p}/\tau_{g}=0.1$. However, while numerical
simulations indicate that the vortex state is stable in the former case, it
remains to investigate its stability in the latter scenario. In this case
again, it is straightforward to perform a linear stability analysis through
the Jacobian matrix:

$\mathrm{J}=\begin{pmatrix}\mathrm{diag}(\bar{\mathbf{g}})-\mathrm{I}-\mathrm{Q}&\mathrm{diag}(\bar{\mathbf{a}})&\mathbf{0}\\\
-\tau\mathrm{diag}(\bar{\mathbf{a}}\odot\bar{\mathbf{g}})^{*}&-\tau(\mathrm{I}+\mathrm{diag}(\bar{\mathbf{a}}\odot\bar{\mathbf{a}}^{*}))&-\tau\mathrm{diag}(\bar{\mathbf{a}}\odot\bar{\mathbf{g}})\\\
\mathbf{0}&\mathrm{diag}(\bar{\mathbf{a}})^{*}&\mathrm{diag}(\bar{\mathbf{g}})-\mathrm{I}-\mathrm{Q}\end{pmatrix}$
(13)

where, $\bar{\mathbf{a}}=(\bar{a}_{1},\cdots,\bar{a}_{n})^{t}$,
$\bar{\mathbf{g}}=(\bar{g}_{1},\cdots,\bar{g}_{n})^{t}$, and
$\tau=\tau_{p}/\tau_{g}$. By considering the vortex pattern of Fig.
LABEL:fig7(a) as the fixed point solution, the eigenvalues of the Jacobian
matrix were numerically calculated for a range of values of
$\tau_{p}/\tau_{g}$. Figure 9 depicts the real part of all eigenvalues versus
$\tau_{p}/\tau_{g}$, while the lead eigenvalue, i.e., the one with the largest
real part, is highlighted in red. Clearly, all eigenvalues exhibit negative
real parts in the entire range of the time scale ratio, which shows linear
stability of the vortex state in all regimes. However, according to this
figure, the eigenvalues undergo a square-root transition and increase toward
zero for small values of $\tau_{p}/\tau_{g}$. This clearly suggests the
evolution of the vortex from a stable to a metastable state that can disappear
with a small perturbation, which is in agreement with the simulations.

Figure 9: The real part of the Jacobian matrix eigenvalues versus the photon
to gain lifetime ratio $\tau_{p}/\tau_{g}$. These eigenvalues represent the
linear stability test for a vortex pattern (shown as inset) in a $6\times 6$
rectangular lattice array. The largest eigenvalue is highlighted in red.

##  VI Conclusion

In summary, we investigated the formation of vortex and antivortex phase
patterns in laser arrays. We showed that for large gains these patterns exist
in direct analogy with topological defects of the XY model. However, their
stability depends critically on the ratio of photon to gain lifetime
$\tau_{p}/\tau_{g}$. An important finding was that large gain lifetimes
($\tau_{g}$) destabilize topological defects and force the system into the
ground state.

The presence of self-organized vortex and antivortex phase patterns in laser
arrays could have practical applications given that these patterns can emit
optical vortex beams. It is worth noting that both solid-state and
semiconductor lasers are considered class-B lasers given that in these systems
the fluorescence lifetime is by several orders of magnitudes larger than the
photon lifetime ($\tau_{p}/\tau_{g}\ll 1$). A potential route for enforcing
such laser arrays to create vortex phase patterns is to fix the phase of the
boundary elements which could be done by seeding through a master laser. In
this manner, vortex patterns can be enforced by fixing the topological charge
of the array through boundary elements. However, it remains to investigate
this aspect given that pinning the topological charge in the nonlinear
dynamical systems discussed here can result in wandering of a vortex and
potentially other instabilities.

Finally, it is worth mentioning that the formation of bound vortex-antivortex
states suggests an analogy with the Berezinskii–Kosterlitz–Thouless (BKT)
phase transition. In two spatial dimensions (2D), the theorem by Mermin and
Wagner precludes conventional long-range order at any finite temperature in
systems with continuous symmetry and short-range interactions Mermin and
Wagner (1966). Nevertheless, certain 2D systems can exhibit signs of a BKT
transition from a high-temperature disordered state, with short-range
correlations, to a quasi-ordered state, with algebraic correlations, below a
critical temperatureKosterlitz and Thouless (1973); Kosterlitz (1974). When
viewed in terms of the vortices, the BKT transition is a vortex-unbinding
transition. These topological objects, while tightly bound in pairs in the
low-temperature state, unbind and freely proliferate above the critical
temperature, causing the system to melt its quasi-order. While the present
work suggests such a transition, a rigorous investigation of this aspect
requires a non-equilibrium thermodynamic treatment of the laser array which
could be the subject of future studies.

###### Acknowledgements.

We thank Dr. Morrel H. Cohen for bringing the subject of the BKT phase
transition to our attention.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## references

## References

  * Chaikin, Lubensky, and Witten (1995) P. M. Chaikin, T. C. Lubensky, and T. A. Witten,  _Principles of condensed matter physics_ , Vol. 10 (Cambridge university press Cambridge, 1995).
  * Bishop and Reppy (1980) D. Bishop and J. Reppy, “Study of the superfluid transition in two-dimensional he 4 films,” Physical Review B 22, 5171 (1980).
  * Beasley, Mooij, and Orlando (1979) M. Beasley, J. Mooij, and T. Orlando, “Possibility of vortex-antivortex pair dissociation in two-dimensional superconductors,” Physical Review Letters 42, 1165 (1979).
  * Hebard and Fiory (1980) A. Hebard and A. Fiory, “Evidence for the kosterlitz-thouless transition in thin superconducting aluminum films,” Physical Review Letters 44, 291 (1980).
  * Pargellis, Green, and Yurke (1994) A. N. Pargellis, S. Green, and B. Yurke, “Planar xy-model dynamics in a nematic liquid crystal system,” Physical Review E 49, 4250 (1994).
  * Halperin and Nelson (1978) B. Halperin and D. R. Nelson, “Theory of two-dimensional melting,” Physical Review Letters 41, 121 (1978).
  * Nelson and Halperin (1979) D. R. Nelson and B. Halperin, “Dislocation-mediated melting in two dimensions,” Physical Review B 19, 2457 (1979).
  * Nixon  _et al._ (2013) M. Nixon, E. Ronen, A. A. Friesem, and N. Davidson, “Observing geometric frustration with thousands of coupled lasers,” Phys. Rev. Lett. 110, 184102 (2013).
  * Berloff  _et al._ (2017) N. G. Berloff, M. Silva, K. Kalinin, A. Askitopoulos, J. D. Töpfer, P. Cilibrizzi, W. Langbein, and P. G. Lagoudakis, “Realizing the classical XY Hamiltonian in polariton simulators,” Nature Materials 16, 1120–1126 (2017).
  * Parto  _et al._ (2020) M. Parto, W. Hayenga, A. Marandi, D. N. Christodoulides, and M. Khajavikhan, “Realizing spin Hamiltonians in nanoscale active photonic lattices,” Nature Materials (2020), 10.1038/s41563-020-0635-6.
  * Honari-Latifpour and Miri (2020) M. Honari-Latifpour and M.-A. Miri, “Mapping the $xy$ hamiltonian onto a network of coupled lasers,” Phys. Rev. Research 2, 043335 (2020).
  * Ding  _et al._ (2019) J. Ding, I. Belykh, A. Marandi, and M.-A. Miri, “Dispersive versus dissipative coupling for frequency synchronization in lasers,” Phys. Rev. Applied 12, 054039 (2019).
  * Ding and Miri (2019) J. Ding and M.-A. Miri, “Mode discrimination in dissipatively coupled laser arrays,” Optics letters 44, 5021–5024 (2019).
  * Haus (1984) H. Haus, “Waves and fields in optoelectronics.” PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, NJ 07632, USA, 1984, 402 (1984).
  * Suh, Wang, and Fan (2004) W. Suh, Z. Wang, and S. Fan, “Temporal coupled-mode theory and the presence of non-orthogonal modes in lossless multimode cavities,” IEEE Journal of Quantum Electronics 40, 1511–1518 (2004).
  * Tredicce  _et al._ (1985) J. R. Tredicce, F. T. Arecchi, G. L. Lippi, and G. P. Puccioni, “Instabilities in lasers with an injected signal,” JOSA B 2, 173–183 (1985).
  * Fabiny  _et al._ (1993) L. Fabiny, P. Colet, R. Roy, and D. Lenstra, “Coherence and phase dynamics of spatially coupled solid-state lasers,” Physical Review A 47, 4287 (1993).
  * Ohtsubo (2012) J. Ohtsubo,  _Semiconductor lasers: stability, instability and chaos_ , Vol. 111 (Springer, 2012).
  * Mandel and Wolf (1995) L. Mandel and E. Wolf,  _Optical coherence and quantum optics_ (Cambridge University Press, 1995).
  * Mahler  _et al._ (2019) S. Mahler, V. Pal, C. Tradonsky, R. Chriki, A. A. Friesem, and N. Davidson, “Dynamics of dissipative topological defects in coupled phase oscillators,” Journal of Physics B: Atomic, Molecular and Optical Physics 52, 205401 (2019).
  * Aranson and Kramer (2002) I. S. Aranson and L. Kramer, “The world of the complex ginzburg-landau equation,” Reviews of modern physics 74, 99 (2002).
  * Hagan (1982) P. S. Hagan, “Spiral waves in reaction-diffusion equations,” SIAM journal on applied mathematics 42, 762–786 (1982).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one-or two-dimensional isotropic heisenberg models,” Physical Review Letters 17, 1133 (1966).
  * Kosterlitz and Thouless (1973) J. Kosterlitz and D. Thouless, “J. phys. c: Solid state phys.” (1973).
  * Kosterlitz (1974) J. Kosterlitz, “The critical properties of the two-dimensional xy model,” Journal of Physics C: Solid State Physics 7, 1046 (1974).
  * Lee  _et al._ (2010) T. E. Lee, H. Tam, G. Refael, J. L. Rogers, and M. Cross, “Vortices and the entrainment transition in the two-dimensional kuramoto model,” Physical Review E 82, 036202 (2010).