We study the properties of nonlinear Bloch waves in a diamond chain waveguide lattice in the presence of a synthetic magnetic flux. In the linear limit, the lattice exhibits a completely flat (wavevector *k*-independent) band structure, resulting in perfect wave localization, known as Aharonov–Bohm caging. We find that in the presence of nonlinearity, the Bloch waves become sensitive to *k*, exhibiting bifurcations and instabilities. Performing numerical beam propagation simulations using the tight-binding model, we show how the instabilities can result in either the spontaneous or controlled formation of localized modes, which are immobile and remain pinned in place due to the synthetic magnetic flux.

## I. INTRODUCTION

Nonlinear waves are general solutions of nonlinear partial differential equations such as the Korteweg–de Vries, cubic nonlinear Schrödinger, sine-Gordon, or Boussinesq equations,^{1} which describe wave processes in diverse fields ranging from nonlinear optics, hydrodynamics, and plasma to biological systems; the theory of general relativity; and the topological theory of knots. Nonlinear wave systems can exhibit singularities, coherent structures, solitary waves, and instabilities. The latter is responsible for the phenomena of wave breaking in hydrodynamics,^{2} filamentation of laser beams in nonlinear media,^{3} wave collapse,^{4} Langmuir waves in plasma,^{5} and ocean turbulence.^{6}

Studies of nonlinear wave dynamics in periodic lattice systems date from the investigations of Fermi, Pasta, Ulam, and Tsingou,^{7,8} who posed a question regarding energy relaxation in a one-dimensional string of nonlinearly interacting particles. More recently, the study of nonlinearities in periodic photonic structures attracts broad interest, offering new routes to manipulate light propagation.^{8–11} Apart from supporting lattice solitons, nonlinearities in lattices can induce coupling between different Bloch modes and bands, giving rise to fundamental phenomena such as modulation instability^{12,13} and spontaneous pattern formation within and/or between different bands.^{14} In short, nonlinear wave dynamics properties can be manipulated by engineering the photonic band structure and vice versa.

Nowadays, a new front for nonlinear wave dynamics in lattices has been opened by recent advances in creating artificial gauge fields for light, which can be used to manipulate photonic band structures to control their dispersion and topological properties,^{15–18} leading to many opportunities for both fundamental studies and device applications, such as the creation of edge states protected against backscattering. The number of studies of nonlinear effects in the presence of artificial gauge fields is growing rapidly,^{19–23} but a general understanding of how nonlinear phenomena differ from those in conventional systems is still far from being complete.

The aim of this study is to shed further light on the effect of artificial gauge fields on nonlinear wave dynamics in photonic lattices, focusing on the diamond chain lattice. Application of synthetic gauge flux to the diamond chain can give rise to a peculiar band structure where all three of its Bloch bands are completely flat (independent of the wavevector *k*), known as an Aharonov–Bohm cage.^{24,25} The fully flat spectrum results in perfect (compact) localization of arbitrary localized excitations of the lattice, which has been observed using fs laser-written waveguide (WG) arrays.^{26–28} Recent theoretical studies indicate that this Aharonov–Bohm caging can persist for localized excitations in the presence of nonlinear interactions.^{29–31}

In this article, we study the effect of nonlinearity on delocalized nonlinear Bloch wave excitations of the Aharonov–Bohm cage lattice, comparing against the properties of the flux-free diamond chain lattice. We find that nonlinearity generically lifts the degeneracy of the flat band Bloch waves, resulting in nonlinear Bloch wave solutions, which are dispersive and have *k*-dependent stability properties. At low intensities, two of the bands are linearly stable for small wavenumbers (in the vicinity of *k* = 0), while at larger wavenumbers, all bands exhibit instabilities, strongest at *k* = *π*/2, which result in the spontaneous formation of localized modes. To complement the linear stability analysis, we carry out numerical simulations of light propagation under the tight-binding approximation. We show how the combination of the synthetic gauge flux and nonlinearity offers the ability to shape the optical field profile using controlled perturbations, enabling the generation of localized nonlinear modes from the nonlinear Bloch waves.

The outline of this article is as follows: Sec. II introduces the tight-binding model for the nonlinear diamond chain in the presence of a synthetic gauge flux, computing its nonlinear Bloch modes and their linear stability. In Sec. III, we validate the linear stability analysis against direct numerical propagation simulations of perturbed nonlinear Bloch waves. Section IV shows how the application of a controlled perturbation to unstable nonlinear Bloch waves can be used to create localized nonlinear modes in a controlled way. Section V concludes this paper.

## II. MODEL AND NONLINEAR BLOCH MODES

We consider light propagation in the magnetic flux-dressed diamond chain with Kerr nonlinearity under the tight-binding (nearest neighbor coupling) approximation, illustrated schematically in Fig. 1(a), described by the discrete nonlinear Schrödinger equation^{29,30}

where *z* is the propagation coordinate (in units of the inter-waveguide coupling strength). We write the wave field as a three component vector with complex amplitudes $\psi n=(an\u2009\u2009bn\u2009\u2009cn)T$, *g* is the Kerr nonlinearity strength, and *ϕ* is the artificial magnetic flux, which can be introduced to *fs* laser-written waveguide arrays by bending the waveguides, introducing auxiliary waveguides, or using higher order waveguide modes carrying orbital angular momentum.^{26–28} The model equations can be derived from the system Hamiltonian *H* = *H*_{L} + *H*_{NL}, where

are the linear and nonlinear energies, respectively, and *N* is the total number of unit cells in the lattice, and we assume periodic boundary conditions. *H* is conserved, along with the total norm $N=\u2211n=1N(|an|2+|bn|2+|cn|2)=1$. We fix $N$ without loss of generality as changes in $N$ are equivalent to rescaling *g*.

In the following, we consider parameter values representative of experiments in fs laser-written waveguide arrays,^{22,26–28,32} where the typical inter-waveguide coupling strength is *κ* = 0.85 cm^{−1} and propagation lengths up to *L* = 10 cm are feasible,^{27} corresponding to a normalized total propagation length of *κL* = 8.5. For waveguides inscribed in fused silica glass, the effective nonlinear coefficient is *γ* = 1.7 cm^{−1} MW^{−1}, and probe beam powers up to *P* = 4 MW are accessible.^{32} Consequently, the normalized nonlinear coefficient *g* = *Pγ*/*κ* can go up to 8. Note that in this study, we are mainly focused on the weakly nonlinear regime *g* < 1, which gives the freedom to potentially increase the coupling strength to observe longer dimensionless propagation lengths.

To obtain the Bloch band structure of Eq. (1), we transform to Fourier space and seek stationary solutions of the form *ψ*_{n}(*z*) = (*a*_{0}, *b*_{0}, *c*_{0})*e*^{i(kn−λz)}, where *k* is the wavevector and *λ* is the propagation constant along *z*. This yields the set of nonlinear algebraic equations

To solve this system, we decompose the site amplitudes into their intensities and phases,

where the *b* site amplitude is fixed by our choice of normalization $N$, and we have set its phase to zero. After straightforward algebraic calculations, we find the expressions for the phases *θ*_{a} = −(*ϕ*/2 + *k*)/2 ± *π*/2, *θ*_{c} = (*ϕ*/2 − *k*)/2 ∓ *π*/2. Substituting these solutions back into Eq. (3), we obtain a set of coupled real polynomial equations that can be solved numerically to obtain the eigenvalues (propagation constants) *λ* = *λ*(*k*, *ϕ*, *g*).

The Bloch wave spectrum in the linear limit *g* = 0 is shown in Figs. 1(c) and 1(d). When *ϕ* ≠ *π*, the spectrum consists of one flat band at *λ* = 0 and two *k*-dependent dispersive bands. Owing to the degeneracy of the flat band, one can construct compact localized eigenmodes as a superposition of its Bloch waves.^{29,30} When *ϕ* = *π*, which corresponds to the Aharonov–Bohm caging, the other two bands also become flat, corresponding to a *k*-independent spectrum *λ* = 0, ±2. In this limit, any excitation of the system can be decomposed in terms of compact localized modes, and therefore, it remains compactly localized for all *z*.

Next, we consider the properties of the weakly nonlinear Bloch modes in Figs. 1(e) and 1(f). In the flux-free case, *ϕ* = 0, the nonlinear bands keep the form of their linear counterparts, being only shifted upward in the presence of weak *g* > 0 with respect to the linear bands. Thus, the middle band remains flat and touches the dispersive bands at the BZ boundaries. On the other hand, when *ϕ* = *π*, the presence of nonlinearity immediately makes all the Bloch bands dispersive, with the central band more strongly affected than the outer bands. Thus, nonlinearity lifts the degeneracy of the flat bands of the Aharonov–Bohm cage.

Figure 2 illustrates the modal properties as a function of the nonlinearity strength *g*. At small *g*, the deviation of the magnitude of the corresponding nonlinear eigenvalues from the linear ones is proportional to the nonlinearity strength. At stronger *g*, we observe the formation of new mode branches due to nonlinear bifurcations. New branches bifurcate from the upper band for strong *g* > 1 in the flux-free lattice and from all bands in the *π* flux lattice. The upper band bifurcation occurs at *g* = 4 at the neighborhood of *k* = 0 and moves toward smaller values of *g* as *k* increases. In general, bifurcations occur for smaller values of *g* in the *π* flux lattice, compared to the flux-free lattice.

We investigate the mode stability by applying the standard linear stability analysis,^{30} which characterizes the evolution of small perturbations with wavevector *p* to the nonlinear Bloch waves. The imaginary part of the perturbation mode eigenvalues *λ*_{p} determines the growth rate of the corresponding perturbation; any positive Im(*λ*_{p}) implies that the nonlinear Bloch wave is unstable. Further details are provided in the Appendix.

The growth rate of the most unstable perturbation mode is plotted in Fig. 2. The lower band in the flux-free lattice is stable, while the other bands are unstable, with instability growth rates insensitive to *k* in the region of weak nonlinearity. On the other hand, in the *π* flux lattice, the mode stability becomes strongly *k*-dependent. The two outer bands are stable in the vicinity of *k* = 0, while the middle band is unstable throughout the Brillouin zone. The bifurcations of the outer bands give rise to new stable modes. They also appear as loops around |*k*| = *π* in Fig. 1(e).

In summary, we have found that the nonlinear Bloch waves in the *π* flux lattice are sensitive to both *k* and *g*, exhibiting *k*-dependent linear stability properties and nonlinear bifurcations. Thus, nonlinearity lifts the degeneracy of the Bloch wave spectrum in the Aharonov–Bohm cage lattice.

## III. DYNAMICS OF NONLINEAR BLOCH WAVES

To validate the stability properties discussed in Sec. II, we now consider the propagation dynamics of the nonlinear Bloch waves. We solve the discrete nonlinear Schrödinger equation (1) numerically. We use the fourth order Runge–Kutta method, with perturbed nonlinear Bloch waves of the form Ψ_{n} + *δψ*_{n} as the initial state at *z* = 0. Here, *δψ*_{n} = *ξ*_{n} exp(*i θ*_{n}), with *ξ*_{n} and *θ*_{n} being random amplitudes and phases, respectively.

To characterize the development of instabilities, we use the inverse participation ratio, ⟨*P*⟩, averaged with respect to different realizations of the initial random perturbation field and integrated over *z*, which is a standard measure for characterizing the localization of a wave field,

where *ℓ* is the propagation length and ⟨..⟩ denotes the ensemble averaging, and |*ψ*_{n}|^{2} = |*a*_{n}|^{2} + |*b*_{n}|^{2} + |*c*_{n}|^{2}. The minimum value of ⟨*P*⟩ (≈1) is attained for beams with a spatially uniform intensity. Over the course of the instability development, the fragmentation of the nonlinear Bloch wave will result in an increase in *P*. Thus, ⟨*P*⟩ provides an independent estimate of the strength of the instabilities.

In Fig. 3, we compare the maximal perturbation growth rate obtained from the linear stability analysis with the time-averaged participation number for the three bands. In general, the dynamical calculations confirm the results of the linear stability analysis. The eigenvalues and instability growth rates are symmetric about *k* = *π*/2. We observe stable propagation of nonlinear waves from the outer bands of the *π*-flux diamond chain lattice at (and in the neighborhood of) the Brillouin zone center *k* = 0 and boundary *k* = *π*. Moreover, in the presence of weak nonlinearity, the small perturbation growth rates result in the persistence of the nonlinear Bloch waves over the short propagation distances accessible in experiments, characterized by ⟨*P*⟩ ≈ 1. As *k* → *π*/2, the perturbation growth rate increases for all three bands, and ⟨*P*⟩ attains its maximal value. Figures 4(a), 4(c), and 4(e) quantify this *k*-dependence of ⟨*P*⟩ for different nonlinearity strengths.

In Figs. 4(b), 4(d), and 4(f), the dynamics of the weakly unstable nonlinear modes for the three different bands in the *π* flux lattice with *g* = 0.3 are illustrated. We observe that the final stage of the dynamics exhibits self-trapping, with the formation of localized intensity peaks, with random positions dictated by the initial perturbation. The localization is reflected in the integral participation ratio profile [Figs. 4(a), 4(c), and 4(e)] averaged over different initial random perturbation field realizations. Interestingly, the intensity maxima are not independent of *z* but exhibit persistent periodic oscillations. This is in contrast to regular (non-flat band) nonlinear lattices, which typically exhibit stationary solitons under the modulational instability.

Upon closer inspection, we find that the localized breathing structures have a transverse (*n*) periodicity of four lattice cells, with *z*-periodicity proportional to the energy difference between the flat bands. Thus, the peculiarity of the nonlinear Bloch wave propagation in the *π* flux weakly nonlinear diamond chain lattice is its ability to spontaneously form localized periodic breathers. The *b* “hub” sites play a pivotal role in the formation of these breathers due to their ability to block wave spreading even in the presence of local nonlinearities.^{29,30} Thus, after the development of the instability, the resulting localized peaks remain confined (caged) between two *a* or *c* sites. In the flux-free lattice *ϕ* = 0, we do not observe the formation of periodic breathers.

## IV. CONTROLLED GENERATION OF LOCALIZED NONLINEAR MODES

The results of Sec. III indicate the ability of the Aharonov–Bohm cage to trap the random localized breathers generated by the instability. Here, we perform a numerical experiment demonstrating the possibility to seed the modulational instability in order to generate breathers in a controlled way by perturbing the initial nonlinear Bloch wave. This provides a mechanism to create a structured light beam from a uniform intensity Bloch wave using nonlinearity.

We inject the nonlinear wave corresponding to the mode from the lowest band in the lattice with *N* = 40, *g* = 0.3, which is characterized by the amplitude $\Psi n=(\u22121/2,1/2,0)$. As in Sec. III, we include a small random perturbation to account for imperfections of the input laser beam, which are inevitably present in experiment. In addition, we kick the input nonlinear Bloch wave by a weak spatially localized wavepacket with *k* = 0 [*ξ* in Fig. 1(b)], *δξ*_{n} = (*ξ*_{a,n}, *ξ*_{b,n}, *ξ*_{c,n})exp(*i η*), where *ξ*_{j,n} are randomly distributed amplitudes and *η* is a fixed phase. Such a kick can be generated by injecting the perturbing beam at normal incidence (*k* = 0) into the lattice after passing it through an amplitude mask. The spatial region in which the kick is applied is indicated by arrows in Figs. 5(a) and 5(b).

Without the initial kick, an irregular pattern of light spots was formed during the propagation, as shown in Fig. 4(d). By turning on the kick, an intensity modulation with a period of four lattice cells is provided to seed the instability, which grows rapidly and generates a regular array of intensity peaks, while the intensity in the un-kicked region remains much weaker, as illustrated in Fig. 5(a). We repeat the simulations for the lattice in the absence of the synthetic magnetic flux in Fig. 5(b). Here, instead of the formation of a periodic wave train, the envelope spreads beyond the initially perturbed region, and small irregular oscillations persist across the entire lattice. This indicates the decisive role of the synthetic magnetic flux and Aharonov–Bohm caging in driving the structured light formation.

To further characterize the wave localization following the development of the instability, we consider the spatial coherence of the final state characterized by the spatial correlation function

where ⟨..⟩ denotes averaging over different realizations of the initial uniform random perturbations and |*ψ*_{n}⟩ is a three component vector (Sec. II). |*C*_{n,m}| = 1 indicates perfect coherence between the two cells, with relative phase independent of the initial random perturbation, while |*C*_{n,m}| ≈ 0 indicates a random field profile sensitive to the initial perturbation.

Figure 5(c) shows the spatial correlations between the middle unit cell (*n* = 20) and the other unit cells. In the case of the *π* flux diamond chain without the initial kick, there is low coherence between different unit cells, whose amplitudes are sensitive to the local initial perturbations [blue line in Fig. 5(c)]. On the other hand, when the kick is applied, strong spatial correlations persist within the kicked region (red line), indicating the caging of the generated localized structures. In the flux-free case, the spatial correlations remain high and do not decay with the separation between the sites (orange line), indicating the preservation of long-range order.

To confirm this interpretation, we compare the averages and standard deviation (SD) of the final intensity obtained for an ensemble of 1000 random initial perturbations for the flux and flux-free cases in Fig. 5(e). In general, the flux-free case *ϕ* = 0 is characterized by a lower SD than the *ϕ* = *π* flux case due to its slower development of the modulation instability for the selected set of parameters. However, in this case, the SD level is the same for both kicked and non-kicked areas, while in the *ϕ* = *π* case, the values of SD are significantly reduced in the kicked region as compared to the non-kicked one.

The peculiarity of the nonlinear train dynamics can be summarized as follows: The light field is modulated by slight kick/tilt at *z* = 0, which drives the restoration of the cage-shaped periodicity. Thus, the hump states *b* preserve their pivotal role in limiting the spreading of energy in the Aharonov–Bohm cage lattice, resulting in a characteristic four-cell spatial periodicity, shown in Fig. 5(d). This four-cell periodicity originates from the periodicity of the initial field with wavenumber *k* = *π*. By changing the phase of the kicking beam, the positions of the cage centers can be shifted from one *b* site to the neighboring sites, preserving the overall spatial periodicity. The center of each hump is therefore in each fourth *b* site in the train, and the localized train profile can be interpreted as the result of the interplay between the linear light trapping (Aharonov–Bohm caging) and the nonlinear localization (self-trapping). This opens interesting opportunities for the experimental design of localized wave trains in waveguide arrays.

One way to generate the requisite initial field profiles in experiments is using a spatial light modulator (SLM),^{33,34} which enables manipulation of both the amplitude and phase profile of the probe beam. Alternatively, to achieve the required intensity modulation, one could shift the initial position (*z* = 0) of the c-sublattice waveguides slightly deeper into the glass sample such that the input probe beam only couples strongly into the a and b sublattices. Then, the relative phase between the *a* and *b* sublattices can be controlled via the angle of the input beam or by adjusting the optical path length of the *a* and *b* sublattices at the input facets.

The perturbation can be created either from the noise inevitably present in experiments or by simultaneous injection of a second perturbation field *δψ*_{n}, as illustrated schematically in Fig. 1(b). After propagation through the sample, the beam intensity profile may be measured either in the near or far-field.

## V. CONCLUSIONS

Managing the light propagation in periodic media by nonlinearity is one of the main issues in photonics. We have studied the dynamics of nonlinear waves in the diamond chain waveguide lattice in the presence of a synthetic magnetic flux and nonlinearity, which exhibits perfect wave localization (Aharonov–Bohm caging) in the linear limit. We considered the persistence of the lattice’s Bloch waves in the presence of nonlinearity. Applying the linear stability analysis and direct numerical simulations of the nonlinear Bloch wave propagation, we have shown how the dispersionless (*k*-independent propagation constant) Bloch waves in the linear limit become strongly *k*-dependent in the presence of nonlinearity; some ranges of wavenumbers exhibit stable propagation within certain ranges of the nonlinearity strength, while others are always unstable. In the weakly nonlinear regime, we have demonstrated the generation of transient nonlinear wave trains as a result of the modulational instability and the way they can be created in a controlled way by suitably perturbing the nonlinear Bloch waves.

A promising experimental platform for observing the nonlinear Bloch wave dynamics studied here is offered by fs laser-written waveguide arrays.^{26–28} The schemes for creating photonic Aharonov–Bohm cages based on the dynamic modulation^{22} or auxiliary waveguides^{27} can be used in the nonlinear regime as long as the nonlinear energy shifts remain small compared to the modulation frequency or auxiliary waveguide detuning, which are typically much stronger than the effective inter-site coupling strength. On the other hand, in the case of the orbital angular momentum mode-based method of Ref. 28, our tight-binding model is not likely to be applicable due to the occurrence of two degenerate copies of the AB cage Hamiltonian, which will become coupled due to the Kerr nonlinearity.

For our study, we considered a specific choice of the synthetic gauge potential. While the linear and nonlinear eigenvalue spectra are independent of the gauge choice (provided the gauge-invariant flux remains the same), for a fixed initial perturbation, the details of the modulation instability dynamics can be sensitive to the choice of gauge, as has recently been theoretically considered in Ref. 35. Thus, depending on the specific gauge realized in experiments, one may need to adjust the form of the controlled perturbation in order to create the nonlinear wave trains.

## AUTHORS’ CONTRIBUTIONS

N.C. contributed equally to this work.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## ACKNOWLEDGMENTS

This research was supported by the Institute for Basic Science in Korea (Grant Nos. IBS-R024-D1 and IBS-R024-Y1). N.C. would like to thank the China Scholarship Council (CSC) for financial support. A.M. acknowledges support from the Ministry of Education, Science and Technological Development of the Republic of Serbia, Grant No. 451-03-68/2020-14/200017. D.L. and D.G.A. acknowledge support from the National Research Foundation, Prime Minister’s Office, Singapore; the Ministry of Education, Singapore, under the Research Centres of Excellence programme; and the Polisimulator project co-financed by Greece and the EU Regional Development Fund. S.K. acknowledges support from the NSFC, Grant Nos. 11674026 and 11974053.

### APPENDIX: LINEAR STABILITY ANALYSIS

The initial development of the modulation instability can be studied using the linear stability analysis,^{30} which can be used to calculate the growth rate of small perturbations to a nonlinear mode. It starts by considering the evolution of a small perturbation added to a stationary state,

where $\Psi n=(an\u2009\u2009bn\u2009\u2009cn)T$ is the stationary mode profile with propagation constant *λ* and $\delta \psi n(z)=[\delta an\u2009\u2009\delta bn\u2009\u2009\delta cn]T$ is a weak perturbation. By substituting Eq. (A1) into the model equations [Eq. (1)] and linearizing the nonlinear terms for *g*|*δψ*_{n}|^{2} ≪ 1, a set of linear ordinary differential equations for the perturbations can be obtained,

where *H*_{L} was given by Eq. (2). The solution to this equation can be decomposed into perturbation modes of the form

where *ν*_{n}, *η*_{n} are three component vectors, *p* is the wavevector of the small perturbations, and *λ*_{p} are the complex eigenvalues indicating the stability properties of small perturbations. Collecting terms with the same time dependence and taking the complex conjugate of one set of terms, we obtain the linear stability eigenvalue problem $\lambda p\delta \psi =M\u0303\delta \psi $, where

while the eigenvalue matrix $M\u0303$ can be explicitly written as

and the other parameters were defined in the main text in Sec. II. The imaginary part of *λ*_{p} is proportional to the instability growth rate of the perturbation. The pure imaginary *λ*_{p} denotes the exponential instability, while the complex *λ*_{p} with nonzero imaginary part denotes the oscillatory instability.