1 Introduction↩︎

When a quantum particle in a periodic potential is subjected to a weak linear force, Bloch oscillations [1] can be observed. This fundamental phenomenon can be interpreted as the adiabatic motion of the particle in the reciprocal space, during which it follows the periodic dispersion relation representative for the band-gap spectrum of periodic potentials. Being ubiquitous across different areas of physics, Bloch oscillations have been experimentally observed in condensed matter [2]–[5], atomic [6]–[13], Dahan1996?, and diverse optical [14]–[19] systems.

In recent years, there has been a growing interest in the investigation of various dynamical phenomena in structured aperiodic potentials, conducted using setups similar to those employed for periodic media. For the one-dimensional (1D) quasi-periodic potentials, prominent examples include the evolution of atoms in optical lattices [20], [21], nontrivial beam dynamics in arrays of optical waveguides [22] and induced lattices in photosensitive materials [23]. A natural question that arises is: what would be the dynamics of a wavepacket in such quasi-periodic systems under the influence of a weak gradient? When a quasi-periodic potential is formed by two sublattices of substantially different depths, such that one of the sublattices can be treated as a perturbation, the tight-binding approximation can be applied to the deeper lattice, giving rise to discrete models where the incommensurate lattice acts as a modulation. This approach leads to the well-known Aubry-André model [24]. The study of Bloch oscillations induced by an additional weak external field in such 1D aperiodic systems has been reported in Refs. [25], [26]. In the fully continuous 1D models with incommensurate lattices, damping of Bloch oscillations was explored in [27], where the aperiodic lattice was also approximated as a deep periodic potential with a weak incommensurate perturbation. Experimental studies of atomic Bose-Einstein condensates (BECs) in a bichromatic lattice with incommensurate periods were performed in [20], where a reduction of damping due to interatomic interactions was observed. The theory of Bloch–Landau–Zener (BLZ) oscillations of localized modes in 1D quasiperiodic lattices with comparable potential depths was developed in [28].

If a quasi-periodic potential is generated by two (or more) incommensurate lattices of comparable amplitudes, neither the tight-binding approach nor the perturbative treatment of one of the constituent sublattices is applicable. In this case, the Hamiltonian spectrum no longer exhibits a band-gap structure, but instead acquires fractal characteristics [29], thereby losing one of the key features of periodic systems typically employed to interpret Bloch oscillations. Furthermore, when the potential is not shallow, two main types of eigenstates emerge: spatially localized and extended ones, separated by a mobility edge (ME) [30]–[35] (even more general situations where several MEs, as well as critical states can exist, will not be addressed here). Nevertheless, even under these conditions, a weak linear gradient can still induce oscillatory dynamics, although the nature of these oscillations differs from those in periodic systems. For the 1D localized states, it was shown in Ref. [28] that the oscillatory dynamics arises from tunneling among (quasi-)resonant eigenstates, and in this sense it closely resembles the phenomenon of Landau-Zener tunneling [36], LandauTun1932?, which has been experimentally studied in systems of cold atoms loaded in periodic lattices [37]–[41]. Consequently, the resulting dynamics have been termed Bloch-Landau-Zener (BLZ) oscillations, the term that will be used below.

The physics underlying Bloch oscillations in 1D periodic structures and BLZ dynamics in 1D quasi-periodic systems differ qualitatively, resulting in markedly distinct dynamical behaviors [28]. For instance, the amplitude and frequency of BLZ oscillations depend on both the initial conditions and the applied gradient. Unlike standard Bloch oscillations, BLZ behavior can involve more than two eigenstates and may feature multiple characteristic frequencies of oscillations. Increasing linear gradient may even suppress BLZ oscillations — a phenomenon that never occurs in periodic systems. Notably, BLZ oscillations also exhibit a remarkable robustness against weak nonlinearities of both types, attractive and repulsive, in contrast to conventional Bloch oscillations.

Bloch oscillations in 2D periodic lattices (continuous, discrete, and semi-discrete) have also been considered in earlier theoretical studies [42]–[53], Khomeriki2010?, and observed experimentally in optically induced lattices [54] and in phononic crystals [55]. The dynamics of such oscillations were touched in [56], where a sophisticated aperiodic potential was modeled using tight-binding approximation, as well as in the bilayer discrete moiré lattice [57], where localized modes were not considered. Generally, BLZ dynamics in continuous 2D aperiodic systems remain largely unexplored, the physical mechanism has not been described, and the modes participating in these processes have not been identified.

The goal of this paper is to extend the theory of BLZ oscillations to 2D aperiodic continuous systems, focusing on the propagation of light beams through photonic moiré arrays. Although the analysis is carried out in the context of photonic systems, the underlying model is broadly applicable, e.g., to cold atoms and BECs in moiré potentials. There are several features of the 2D aperiodic system with gradient, making it qualitatively different from 1D aperiodic systems. These include eventual level crossing of localized states, a much larger variety of the eigenstates of the governing Hamiltonian (say, including states localized along only one spatial direction), coupling scenarios, as well as the presence of additional control parameters like the twist angle of a moiré lattice and the orientation of the applied weak gradient. Nonlinearity in 2D systems can lead to a wider range of instabilities, including collapse or enhanced diffraction.

BLZ oscillations described in our work are readily observable experimentally, as photonic moiré lattices can be created using well-established methods. These include photorefractive crystals [58]–[60], femtosecond-laser written waveguide arrays [61], and photonic crystals with moiré patterns produced via holographic fabrication techniques [62]. The holographic method also enables the creation of a gradient in photonic moiré lattices through the addition of an auxiliary sublattice [63]. A refractive index gradient can also be induced in waveguide arrays via thermo-optic response by applying a thermal gradient to a sample [64], [65].

The paper is organized as follows. Section 2 introduces the problem and describes the system under consideration. In Section 3, we derive the selection rule governing BLZ oscillations, which involves conditions of spatial proximity and quasi-resonance between the participating states. This section also presents examples of periodic energy transfer between such states. Section 4 examines the impact of focusing and defocusing nonlinearities on BLZ dynamics. Finally, in Section 5, we discuss several aspects related to the experimental feasibility of observing 2D BLZ oscillations.

2 Statement of the problem↩︎

We consider the propagation of a paraxial optical beam along the \(z\) direction in the material with shallow transverse modulation of the refractive index, governed by the dimensionless Schrödinger equation for the light field amplitude \(\Psi\) (for normalizations see [59]): \[\begin{align} \label{main} i\frac{\partial \Psi}{\partial z}=H_{{\boldsymbol{\alpha}}}\Psi , \end{align}\tag{1}\] where the Hamiltonian is given by \[\begin{align} \label{Hamilt} H_{{\boldsymbol{\alpha}}}= -\frac{1}{2}\nabla^2-U({\boldsymbol{r}})-{\boldsymbol{\alpha}}\cdot{\boldsymbol{r}} , \end{align}\tag{2}\] where \({\boldsymbol{r}}=(x,y)\) and \(\nabla=(\partial_x,\partial_y)\). The optical potential \(U({\boldsymbol{r}})\) represents an incommensurate (i.e., aperiodic) moiré lattice [58]–[60] induced by two sublattices. These sublattices are mutually rotated by a non-Pythagorean angle \(\theta=\pi/6\), implemented via the 2D rotation operator \(R(\theta)\), and shifted relative each other by the displacement vector \(\boldsymbol{s}\): \[\begin{align} \label{pot95total} U({\boldsymbol{r}})=p_1V(R(\theta){\boldsymbol{r}}+\boldsymbol{s})+p_2V({\boldsymbol{r}}). \end{align}\tag{3}\] Here \(p_{1,2}\) are the depths of the sublattices (it is assumed that \(|V({\boldsymbol{r}})|\leq 1\)). The linear gradient is characterized by the vector \({\boldsymbol{\alpha}}=(\alpha_x,\alpha_y)\), which is considered small enough, i.e., \(\alpha=|{\boldsymbol{\alpha}}|\ll p_{1,2}\). The shift \(\boldsymbol{s}\) is introduced to break the discrete rotational symmetry of the lattice. Neither the choice of a non-Pythagorean angle, nor the shift \(\boldsymbol{s}\) are essential for the consideration, provided the respective lattice is incommensurate phase (say, \(\boldsymbol{s}\) can acquire zero value).

The respective stationary eigenvalue problem, obtained by using the ansatz \(\Psi=e^{ibz}\varphi({\boldsymbol{r}})\), where \(b\) is the propagation constant, has the form \[\begin{align} \label{Halpha} H_{{\boldsymbol{\alpha}}}\varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}=-b_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})} \varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}. \end{align}\tag{4}\] We adopt a vector notation for labeling the eigenstates \(\boldsymbol{n}=(n_1,n_2)\) to account for possible degeneracies of the eigenvalues \(b_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}\). Such degeneracies may originate from discrete rotational symmetry (if \(\boldsymbol{s}=(0,0)\)) or occur accidentally [see, e.g., Fig. 2 below]. The first index \(n_1\) numbers the eigenvalues in the descending order of their absolute values, while the second index \(n_2\) numbers eigenfunctions within the respective degenerate subspace. We also explicitly indicate the dependence of the eigenfunctions on the gradient \(\boldsymbol{\alpha}\). In the absence of the gradient, i.e., for \({\boldsymbol{\alpha}}=(0,0)\), one obtains the eigenvalue problem \[\begin{align} H_0\varphi_{\boldsymbol{n}}^{(0)}=-b_{\boldsymbol{n}}^{(0)} \varphi_{\boldsymbol{n}}^{(0)}, \quad H_0=-\frac{1}{2}\nabla^2-U({\boldsymbol{r}}) , \end{align}\] which describes modes of a conventional aperiodic moiré lattice.

Strictly speaking, the spectrum of the Hamiltonian \(H_{{\boldsymbol{\alpha}}}\) with \(\alpha>0\), considered in the entire space \(\mathbb{R}^2\), is expected to be continuous, implying the absence of square-integrable eigenstates (see, e.g., [66]–[68] for rigorous results in 1D). Nevertheless, it is relevant from both experimental and theoretical perspectives to consider a moiré lattice with a large, but finite area \(S\), with zero boundary conditions for eigenmodes. In such a statement, one deals with a discrete spectrum of \(H_{{\boldsymbol{\alpha}}}\), which justifies the use of a discrete two-component index \(\boldsymbol{n}\). On the other hand, as we will see below from a thorough numerical analysis, the moiré lattice with nonzero \({\boldsymbol{\alpha}}\) supports spatially localized states, for which the choice of the specific boundary conditions is not relevant for observed dynamics. The respective eigenmodes can be normalized: \[\begin{align} \int_{S}|\varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}|^2d{\boldsymbol{r}}=1, \end{align}\] where the integration is carried out over the whole lattice area \(S\).

As it was shown in previous numerical and experimental studies of moiré lattices (see e.g. [58]–[60], [69]), when the lattice is sufficiently deep, the eigenmodes of \(H_0\) with largest eigenvalues \(b_{\boldsymbol{n}}^{(0)}\) (the lower states in quantum-mechanical applications) become spatially localized. The indicator of localization of such modes is their form-factor \([\chi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}]^{1/2}\) where \[\begin{align} \chi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}=\langle[\varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}]^2, [\varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}]^2\rangle \end{align}\] (it is also known as the inverse participation ratio), and \(\langle f,g\rangle :=\int_Sf^*gd{\boldsymbol{r}}\). The larger the value of \(\chi\), the stronger the localization. In particular, modes with the largest eigenvalues tend to localize within a small region of the total domain, such that \[\begin{align} \chi_{\boldsymbol{n}}^{-1} \ll S. \end{align}\]

The parameters of the lattice, at which localized modes emerge for a given gradient, correspond to the localization-delocalization transition (LDT). Increase of the lattice depth above this transition leads (in a general situation) to an increase in the number of localized modes. Typically, when the parameters are above the LDT, one can identify a mobility edge, \(b_{\rm ME}^{(\boldsymbol{\alpha})}\), which can be defined as the largest propagation constant of an extended state, such that all modes with \(b_{\boldsymbol{n}}^{(\boldsymbol{\alpha})}>b_{\rm ME}^{(\boldsymbol{\alpha})}\) in the area \(S\) are localized. Although this is not a strict definition, the latter requires information about the modes defined in the entire space \(\mathbb{R}^2\), it is practically useful for characterizing BLZ oscillations, the focus of this work.

The lattice under consideration and the properties of its highest eigenmodes are illustrated in Fig. 1. Fig. 1 (a) shows the lattice profile within a large but finite domain \(S\), which is used in all simulations presented below. We consider an experimentally relevant configuration in which each sublattice \(V\) is composed of Gaussian waveguides arranged in a square array: \[\begin{align} \label{Gaussain} V({\boldsymbol{r}})=-\sum_{{\boldsymbol{m}}\in\mathbb{Q}^2} e^{-({\boldsymbol{r}}-\boldsymbol{m}d)^2/w^2}, \end{align}\tag{5}\] where \(\mathbb{Q}^2\) is the space of 2D vectors with integer components defining the Bravais lattice with period \(d\) and \(w\) is the waveguide radius. Two sublattices are rotated by the angle \(\theta=\pi/6\) corresponding to the incommensurate (aperiodic) moiré structure.

Figure 1 (b) shows positions of localized modes emerging in the lattice without tilt, \(\boldsymbol{\alpha}=0\), characterized by “centers of mass” coordinates defined by \[\begin{align} \label{COM95mode} {\boldsymbol{r}}_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}= \langle \varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}, {\boldsymbol{r}}\varphi_{\boldsymbol{n}}^{({\boldsymbol{\alpha}})}\rangle. \end{align}\tag{6}\] The degree of localization (form-factor) of different modes in this figure is indicated by the color of the corresponding dot. We observe a nearly homogeneous, although irregular, distribution of the localized modes over the whole area of the moiré lattice (this is similar to the behavior of modes in 1D quasi-periodic arrays [28], [70]). Considering this fact, in panels (a) and (b) of Fig. 1 we outline a central area delimited by the circumference of the radius \(\rho=30\), shown as a white solid line in (a) and a black dotted line in (b), inside which the modes used for dynamical simulations are located (this constraint ensures that the boundaries do not affect the dynamics, but it does not have any physical significance, as this area can be further expanded by increasing the integration window). Below, we refer to this domain as the effective area, denoting it as \(S_{\rm eff}\). It is worth emphasizing that its area is \(S= 4900\), while the smallest form-factor shown in Fig. 1 (b) corresponds to \(1/\chi\approx 11.11\), i.e., these modes are well localized within this area (specific examples are shown in the figures below).

Close inspection of the localization properties of the modes reveals that in incommensurate moiré lattices, in contrast to 1D quasi-periodic lattices, the ME is not sharp. Instead, delocalization of modes occurs relatively smoothly as their propagation constant varies. This is illustrated in Fig. 1 (c), where black dots represent the propagation constants \(\beta\) of the modes confined within the effective area \(S_{\rm eff}\), and red squares show corresponding form-factors \(\chi^{1/2}\) for each mode (the notation \(\beta\) for the propagation constant will be explained below). There are no large gaps in the propagation constant spectrum, to which the ME could belong, which also contrasts with the case of 1D quasi-periodic potentials. Roughly speaking, the 200 modes with the largest propagation constants within the effective area can be considered as well localized. The consideration below will be restricted, namely, to these modes. It is important to stress that we consider here only the eigenmodes of \(H_0\) that remain localized in the lattice with gradient described by \(H_{\boldsymbol{\alpha}}\) (recall that \(\alpha\ll 1\)).

Thus, the main goal of this paper is to address the following fundamental question: How does an eigenmode of the Hamiltonian \(H_0\) without a gradient (\(\alpha=0\)), considered as an initial condition for Eq. (1 ), evolve under influence of a weak gradient (\(0<\alpha\ll 1\))? For a linear system, this directly describes the evolution of an initial input beam under the influence of the gradient term \({\boldsymbol{r}}\cdot{\boldsymbol{\alpha}}\).

3 Multi-mode dynamics↩︎

3.1 Localized bases of eigenfunctions↩︎

In the 2D case, the Hamiltonians \(H_0\) and \(H_{\boldsymbol{\alpha}}\) may exhibit degeneracies in their spectra. Moreover, even when there are no degeneracies in the spectrum of \(H_0\) and all modes can be ordered by the magnitude of their eigenvalues at \(\boldsymbol{\alpha}=0\), there is no guarantee that this ordering will be preserved at nonzero \(\boldsymbol{\alpha}\) [see the examples in Fig. 2 (a) and (b)]. To avoid this ambiguity, we introduce a more convenient labeling of localized modes via a single scalar index \(i\). Specifically, we define \(\beta_i^{(0)}=b_{\boldsymbol{n}}^{(0)}\) and set \(i_1 < i_2\) if \(b_{\boldsymbol{n}_1}^{(0)}>b_{\boldsymbol{n}_2}^{(0)}\) (respectively \(\beta_{i_1}^{(0)}>\beta_{i_2}^{(0)}\)). Note that according to the physical meaning of the vector index \(\boldsymbol{n}\) described above, if no degeneracies occur, then one has a simple relation \(i=n_1\). Within degenerate subspaces (i.e., for modes with identical eigenvalues), the order of indices can be chosen arbitrarily. The resulting set of propagation constants \(\beta_i^{(0)}\) is shown in Fig. 1 (c) for \(\alpha=0\) as black dots. The corresponding eigenfunctions and their centers of mass are denoted as \(\varphi_{\boldsymbol{n}}^{(0)}=\phi_{i}^{(0)}\) and \({\boldsymbol{r}}_{\boldsymbol{n}}^{(0)}={\boldsymbol{r}}_{i}^{(0)}\), respectively. The index \(i\) defined in this way is used in Figs. 1 (b) and (c) to label several representative modes employed below to illustrate different dynamical regimes. Importantly, we restrict this labeling to the modes located within the effective area delimited by the circumference, thereby identifying only the modes that are relevant (employed) in subsequent dynamics. This restriction, however, does not imply any theoretical limitation or ambiguity. Indeed, if the integration domain \(S\) is enlarged to \(S'\) or even goes to infinity, many additional modes with eigenvalues above the ME will appear. These new modes may alter the global indexing, but they are typically localized outside the original domain \(S\), i.e., in \(S' \setminus S\). As a result, they do not affect the dynamics of the modes confined well within \(S\). Even potential deformations of modes near the boundary \(\partial S\) have negligible influence on the interior states. This behavior is consistent with previous observations in 2D systems (e.g., [58]) and is discussed in greater detail for the 1D case in [28], [70], [71]. Thus, increasing the integration domain (assuming it is initially large enough) only affects the absolute numbering of the modes, without altering their physical properties, spatial positions, or ordering of eigenvalues within the region of interest.

We emphasize that the ordering of \(\beta_i^{({\boldsymbol{\alpha}})}\) at a given \({\boldsymbol{\alpha}}\) can differ from the one introduced for \(\beta_i^{(0)}\) due to possible level crossing, i.e., it may happen that \(\beta_{{i}_1}^{({\boldsymbol{\alpha}})}<\beta_{{i}_2}^{({\boldsymbol{\alpha}})}\) (or \(b_{\boldsymbol{n}_1}^{({\boldsymbol{\alpha}})}<b_{\boldsymbol{n}_2}^{({\boldsymbol{\alpha}})}\)) for \(i_1<i_2\). An example is provided by the yellow and red branches in Fig. 2 (a), which display the dependence of the propagation constants on \(\alpha\). These two branches undergo two level crossings within the shown interval of \(\alpha\). The yellow branch originates from the state with index 127, while the red branch originates from mode 129. Accordingly, in our labeling scheme, we associate all propagation constants and modes along the yellow branch with index 127, and those along the red branch with index 129, for all values of \(\alpha\). Numerous additional level crossings can also be observed in the gray branches.

Let \(N\) be the number of the localized modes of \(H_0\) (already taking the degeneracy into account) which remain localized in the presence of a nonzero gradient. Considering increasing gradient \((0,0)\to{\boldsymbol{\alpha}}\) as a deformation, one can trace the result of this deformation upon which localized mode \(\phi_{i}^{(0)}\), where \(i=1,...,N\), is transformed into the eigenmode \(\phi_{i}^{({\boldsymbol{\alpha}})}\) of \(H_{{\boldsymbol{\alpha}}}\). Figure 2 (a) shows the dependence of \(b_i^{({\boldsymbol{\alpha}})}\) on \(\alpha\) for several modes with propagation constants in the range \(1.25\)–\(1.45\). Further, we assume that the gradient in polar coordinates is characterized by the angle \(\gamma\) and strength \(\alpha\), i.e., \({\boldsymbol{\alpha}}= (\alpha \cos\gamma, \alpha \sin\gamma)\), and Fig. 2 (a) corresponds to the case of \(\gamma = \pi/3\). Several representative curves in this figure, highlighted in color, correspond to the modes for which the dynamics will be analyzed below. Notably, the modes with indices \(124\) (blue branch) and \(129\) (red branch) exhibit an avoided crossing near \(\alpha \approx 0.003\). In this region, resonant interaction is expected if either of these modes \(\phi_{124}^{(0)}\) or \(\phi_{129}^{(0)}\) is used as the initial condition for propagation in the system with \(\alpha \approx 0.003\). As we will show below, the resonant interaction of two (or more) modes requires yet another condition, namely, a non-negligible coupling of these modes. For the selected modes, this condition is generally not met even at their level crossing with other modes within the effective area. Figure 2 (e) illustrates how the spatial profiles of modes \(124\) and \(129\) change across the avoided crossing. Initially, at \(\alpha = 0\), the blue and red modes have distinct profiles and correspond to points 1 and 4, respectively. Near the avoided crossing (points 2 and 5), both states undergo significant deformation and their amplitude distributions \(|\phi_{124}^{({\boldsymbol{\alpha}})}|\) and \(|\phi_{129}^{({\boldsymbol{\alpha}})}|\) become similar. After the avoided crossing (points 3 and 6), the modes from blue and red branches resemble the red and blue modes at \(\alpha = 0\), respectively, demonstrating an “exchange” of spatial structure. This behavior is further confirmed in Fig. 2 (c), which plots the center-of-mass coordinates \({\boldsymbol{r}}_{i}^{(\alpha)}=\left( x_c,y_c \right)\) of the two modes as functions of \(\alpha\), showing a smooth exchange of position consistent with the avoided crossing scenario.

A more intricate transformation of the states is observed for branches \(127\) (yellow), \(135\) (magenta), and \(137\) (green). As shown in Fig. 2 (a), the magenta branch undergoes two avoided crossings: the first near \(\alpha \approx 0.001\) with the green branch, and the second near \(\alpha \approx 0.005\) with the yellow one. This suggests that during deformation in the parameter space, resonant interactions with different modes may occur at different values of \(\alpha\). The corresponding trajectories of the centers of mass are presented in Fig. 2 (d), also revealing a more complex behavior: mode \(135\) first exchanges its position with mode \(137\), and later with mode \(127\) for larger \(\alpha\). This is consistent with the behavior of the respective propagation constants.

One of the key factors in the 2D case is the direction of the gradient determined by angle \(\gamma\), as it can dramatically alter propagation dynamics. This sensitivity of the modes to the gradient direction is demonstrated in Fig. 2 (b), showing dependence of eigenvalues \(\beta_i^{({\boldsymbol{\alpha}})}\) on \(\alpha\) for angle \(\gamma = -\pi/5\). In comparison with the previous case, the change in gradient direction leads to a complete restructuring of the branches in the spectrum for \(\alpha > 0\). For instance, the blue branch, which previously had a positive slope, now exhibits a negative slope and undergoes two avoided crossings with branches originating from states \(116\) and \(129\). These interactions are absent for the previous angle \(\gamma = \pi/3\). As we will show below, this change in spectral structure has direct consequences for the dynamics: varying the gradient direction can completely change the set of resonantly interacting modes, even when the initial condition remains the same.

If a given localized mode, obtained for \(\alpha=0\), is used at the input in the lattice with a nonzero, but small gradient \(\alpha\), it may, in practice, excite only a few localized modes of the lattice. Thus, to address the central question of this work, one has to determine the modes of \(H_{\boldsymbol{\alpha}}\) that can be efficiently excited by such input. Clearly, due to localization, any newly excited mode should have its center of mass sufficiently close to that of the input mode (see more details below). This observation represents an additional justification for restricting the numerical analysis to the effective area, which in Fig. 1 and Fig. 2 (e) is delimited by the circumference of the radius \(\rho=30\). All modes with appreciable overlap that can participate in the dynamics are confined to the effective area. It is worth noting that the centers of mass of all extended modes, which are located near the origin, also lie within this region. However, such modes are not excited in a system with a weak gradient due to a relatively large difference between their propagation constants and propagation constants of the selected localized states and, therefore, will not be considered further.

Now we introduce two bases \({\boldsymbol{\phi}}^{(0)}\) and \({\boldsymbol{\phi}}^{({\boldsymbol{\alpha}})}\): \[\begin{align} {\boldsymbol{\phi}}^{(0)}=(\phi_1^{(0)}, ...,\phi_N^{(0)})^T \,\, and\,\, {\boldsymbol{\phi}}^{({\boldsymbol{\alpha}})}=(\phi_1^{({\boldsymbol{\alpha}})}, ...,\phi_N^{({\boldsymbol{\alpha}})})^T \end{align}\] (\(T\) stands for transpose) which are considered complete in the subspace of localized modes in the effective area. These bases are connected by the unitary transformation \[\begin{align} \label{unitary} {\boldsymbol{\phi}}^{({\boldsymbol{\alpha}})}=\boldsymbol{S}_{{\boldsymbol{\alpha}}} {\boldsymbol{\phi}}^{(0)}, \qquad S_{ij}^{({\boldsymbol{\alpha}})} =\langle\phi_{j}^{(0)},\phi_i^{({\boldsymbol{\alpha}})}\rangle. \end{align}\tag{7}\] Respectively, one can expand \[\begin{align} \label{expansion} \Psi=\sum_{i=1}^{N}c_i^{(0)}(z)\phi_{i}^{(0)}({\boldsymbol{r}}) =\sum_{i=1}^{N}c_i^{({\boldsymbol{\alpha}})}(z)\phi_{i}^{({\boldsymbol{\alpha}})}({\boldsymbol{r}}) \end{align}\tag{8}\] and introduce the center of mass of the field: \[\begin{align} \label{bR} {\boldsymbol{R}}=\langle \Psi,{\boldsymbol{r}}\Psi\rangle ={\boldsymbol{c}}_0^\dagger{\boldsymbol{R}}^{(0)}{\boldsymbol{c}}_0={\boldsymbol{c}}_{{\boldsymbol{\alpha}}}^\dagger{\boldsymbol{R}}^{({\boldsymbol{\alpha}})}{\boldsymbol{c}}_{{\boldsymbol{\alpha}}}, \end{align}\tag{9}\] where \({\boldsymbol{R}}^{({\boldsymbol{\alpha}})}\) is the matrix with (vector) elements \[\begin{align} \label{r95elem} {\boldsymbol{r}}_{ji}^{({\boldsymbol{\alpha}})}=\langle\phi_{j}^{({\boldsymbol{\alpha}})},{\boldsymbol{r}} \phi_i^{({\boldsymbol{\alpha}})}\rangle, \end{align}\tag{10}\] and \({\boldsymbol{c}}_{{\boldsymbol{\alpha}}}=(c_1^{({\boldsymbol{\alpha}})}, ...,c_N^{({\boldsymbol{\alpha}})})^T\). It follows from (7 ) that \[\begin{align} \label{cSc} {\boldsymbol{c}}^{(0)}= \boldsymbol{S}_{{\boldsymbol{\alpha}}} {\boldsymbol{c}}^{({\boldsymbol{\alpha}})}. \end{align}\tag{11}\]

The evolution of the coefficients is determined by the equations \[\begin{align} \label{dc950} i\frac{d{\boldsymbol{c}}_{0}}{dz}=&-({\boldsymbol{B}}_{0}+{\boldsymbol{\alpha}}\cdot{\boldsymbol{R}}^{(0)}){\boldsymbol{c}}_{0} \end{align}\tag{12}\] in the original basis and \[\begin{align} \label{dc95al} i\frac{d{\boldsymbol{c}}_{{\boldsymbol{\alpha}}}}{dz}=&-{\boldsymbol{B}}_{{\boldsymbol{\alpha}}}{\boldsymbol{c}}_{{\boldsymbol{\alpha}}} \end{align}\tag{13}\] in the basis \(\phi_i^{({\boldsymbol{\alpha}})}\), where \({\boldsymbol{B}}_{{\boldsymbol{\alpha}}}=diag(\beta_1^{({\boldsymbol{\alpha}})},..., \beta_N^{({\boldsymbol{\alpha}})})\). Equation (13 ) is readily solved. Finally, combining (9 ) and (11 ) with the solution of (13 ), one obtains the formula for the evolution of the center of mass: \[\begin{align} \label{COM} {\boldsymbol{R}}(z)={\boldsymbol{c}}_{0}^\dagger(0) \boldsymbol{S}_{{\boldsymbol{\alpha}}}^\dagger {\boldsymbol{\Lambda}}_{{\boldsymbol{\alpha}}}^\dagger(z)\boldsymbol{S}_{{\boldsymbol{\alpha}}}{\boldsymbol{R}}^{(0)} \boldsymbol{S}_{{\boldsymbol{\alpha}}}^\dagger{\boldsymbol{\Lambda}}_{{\boldsymbol{\alpha}}}(z)\boldsymbol{S}_{{\boldsymbol{\alpha}}}{\boldsymbol{c}}_{0}(0), \end{align}\tag{14}\] where \({\boldsymbol{\Lambda}}_{{\boldsymbol{\alpha}}}(z)=diag(e^{i\beta_1^{({\boldsymbol{\alpha}})}z},...,e^{i\beta_N^{({\boldsymbol{\alpha}})}z}).\) We emphasize that the only approximation made so far was the absence of excited modes below the ME and beyond the effective area.

3.2 Selection rule↩︎

Generally, a given input beam can excite multiple localized modes. Thus, to pursue the main objective outlined above, we should first address the following auxiliary task: given an initial state that is an eigenmode of the Hamiltonian \(H_0\) without a gradient (\(\alpha=0\)), determine which modes will be significantly excited for \(0<\alpha\ll 1\) at \(z>0\). This task is addressed below in this subsection.

The initial condition that corresponds to excitation of \(j\)’s mode of \(H_0\) is \(c_i^{(0)}(0)=\delta_{ij}\). Respectively, the amplitude of the state \(\phi_{i}^{({\boldsymbol{\alpha}})}(0)\) at \(z=0\) is determined by \(c_i^{({\boldsymbol{\alpha}})}(0)=S_{ij}^{({\boldsymbol{\alpha}})}\). The modes for which \(c_i^{({\boldsymbol{\alpha}})}(0)\) is negligible are regarded as unexcited. To quantify this, we introduce a parameter \(0<\Delta\ll 1\) and consider only the modes \(j\) with non-negligible projections \[\begin{align} \label{position} |S_{ij}^{({\boldsymbol{\alpha}})}| >\Delta. \end{align}\tag{15}\] Following [28], we call this requirement as spatial proximity conditions. In numerical simulations shown below for identification of the interacting modes, we set \(\Delta= 0.1\). Only the modes satisfying (15 ) can effectively participate in the energy exchange with the input mode \(\phi_j^{(0)}\), while all other modes, even being excited, have negligible weights and do not affect the dynamics. In Fig. 3 (a) and (b), we illustrate the dependence of elements \(S_{ij}^{({\boldsymbol{\alpha}})}\) for all modes \(\phi_i^{({\boldsymbol{\alpha}})}\) from the effective area on gradient \(\alpha\) for the two input modes \(j=124\) (panel (a)) and \(j=137\) (panel (b)). In panel (a), one observes that \(S_{124,124}^{({\boldsymbol{\alpha}})}\) (blue line) remains close to 1 for small \(\alpha\), then decreases monotonically as \(\alpha\) increases. In contrast, \(S_{129,124}^{({\boldsymbol{\alpha}})}\) (red line) is negligible at \(\boldsymbol{\alpha} = 0\), but rises significantly, reaching a value close to \(1\) at \(\alpha = 0.01\). The two curves intersect and take on similar values around the avoided crossing in the spectrum at approximately \(\alpha \approx 0.003\), indicating strong modal overlap. This behavior reflects how the evolving modes for different \(\alpha\) overlap with the initial mode \(\phi_{124}^{(0)}\). At \(\alpha = 0\), the centers of mass of modes \(\phi_{124}^{(0)}\) and \(\phi_{129}^{(0)}\) are spatially separated, as shown in Fig. 2 (c) and (e) at points 1 and 4. As \(\alpha\) increases, the mode profiles deform and their centers of mass approach, as seen at points 2 and 5, leading to comparable values of \(S_{i,124}^{({\boldsymbol{\alpha}})}\) for \(i=124\) and \(i=129\). With further increase in \(\alpha\), the centers of mass separate again (points 3 and 6), and only \(\phi_{129}^{({\boldsymbol{\alpha}})}\) retains significant overlap with the initial mode. For all other localized modes (shown in gray), we obtain \(|S_{{ i,124}}^{({\boldsymbol{\alpha}})}|<0.1\) for all \(\alpha<0.01\) in Fig. 3 (a) and \(|S_{{ i,137}}^{({\boldsymbol{\alpha}})}|<0.1\) for almost all \(\alpha<0.01\) (except for extremely short intervals, which are omitted from our consideration due to their narrowness) in Fig. 3 (b). We also note that in Fig. 3 (a) at \(\alpha\gtrsim 0.005\) the mode 129 has a larger overlap integral with input mode 124, than mode 124 with the same index.

Besides the elements of the overlap matrix \(\boldsymbol{S}_{{\boldsymbol{\alpha}}}\), the evolution of the center of mass \({\boldsymbol{R}}(z)\) of the wavepacket along \(z\) axis, is determined by the exponential factors \(\exp[i(\beta_{i}^{({\boldsymbol{\alpha}})}-\beta_{j}^{({\boldsymbol{\alpha}})})z]\) stemming from the \({\boldsymbol{\Lambda}}_{{\boldsymbol{\alpha}}}(z)\) matrices in expression (14 ). Thus, in addition to (15 ), effectively interacting modes must obey sufficiently close propagation constants, i.e., for \(i\neq j\), one requires \[\begin{align} \label{reson1} |\beta_{i}^{({\boldsymbol{\alpha}})}-\beta_{j}^{({\boldsymbol{\alpha}})}|\ll |\beta_{i}^{({\boldsymbol{\alpha}})}|, |\beta_{j}^{({\boldsymbol{\alpha}})}|. \end{align}\tag{16}\] To express this condition through the characteristics of the modes of the system without gradient (i.e., at \(\alpha=0\), as it is required by the formulated setting of the problem) we recall that the bases are related by the unitary transformation (7 ) and hence, the \(\beta_{i}^{({\boldsymbol{\alpha}})}\) are the eigenvalues also of the matrix \(-({\boldsymbol{B}}^{(0)}+{\boldsymbol{\alpha}}\cdot{\boldsymbol{R}}^{(0)})\) whose diagonal elements are denoted by \[\begin{align} \tilde{\beta}_i^{({\boldsymbol{\alpha}})}=\beta_i^{(0)}+{\boldsymbol{\alpha}}\cdot{\boldsymbol{r}}_{ii}^{(0)}. \end{align}\] For several modes (explored below) the dependencies \(\tilde{\beta}_i^{({\boldsymbol{\alpha}})}\) versus \(\alpha\) are illustrated by the solid lines in Figs. 3 (c) and (d).

Importantly, \(\tilde{\beta}_i^{({\boldsymbol{\alpha}})}\) are determined by the eigenvalues of the untilted lattice and by the gradient. Now we must relate them to the propagation constants \(\beta_{j}^{({\boldsymbol{\alpha}})}\) determining (quasi-)resonances (16 ). To this end, we recall that according to the Gershgorin circle theorem Gershgorin? for each \(\beta_{i}^{({\boldsymbol{\alpha}})}\) one can find \(\tilde{\beta}_{j}^{({\boldsymbol{\alpha}})}\) such that \[\begin{align} \label{energy0} |\beta_{i}^{({\boldsymbol{\alpha}})}-\tilde{\beta}_{j}^{({\boldsymbol{\alpha}})}|\leq \rho_j^{({\boldsymbol{\alpha}})} , \end{align}\tag{17}\] where \[\begin{align} \label{GR} \rho_j^{({\boldsymbol{\alpha}})}=\sum_{j'=1\atop j'\neq j}^{N}|{\boldsymbol{\alpha}}\cdot {\boldsymbol{r}}_{jj'}^{(0)}| \end{align}\tag{18}\] are the Gershgorin radii. Since we consider real spectra, as a matter of fact, \(\rho_j^{({\boldsymbol{\alpha}})}\) are half-lengths of the intervals. In Figs. 3 (c) and (d), they are illustrated by shaded regions around the respective \(\tilde{\beta}_{j}^{({\boldsymbol{\alpha}})}\) lines. Note that generally speaking, \(i\) can be different from \(j\). Thus, if we require that the respective Gershgorin intervals overlap for a given gradient \({\boldsymbol{\alpha}}\) and, i.e., if \[\begin{align} \label{energy} |\tilde{\beta}_{j}^{({\boldsymbol{\alpha}})}-\tilde{\beta}_{j'}^{({\boldsymbol{\alpha}})}|\leq \rho_j^{({\boldsymbol{\alpha}})}+\rho_{j'}^{({\boldsymbol{\alpha}})} \end{align}\tag{19}\] for the modes \(i\) and \(i'\) of \(H_{{\boldsymbol{\alpha}}}\), then one can find the respective modes \(j\) and \(j'\) of \(H_0\) such that \[\begin{align} \label{resonance2} |\beta_{i}^{({\boldsymbol{\alpha}})}-\beta_{i'}^{({\boldsymbol{\alpha}})}|\leq 2 (\rho_{j}^{({\boldsymbol{\alpha}})} +\rho_{j'}^{({\boldsymbol{\alpha}})}) . \end{align}\tag{20}\] Thus, for the gradients small enough, one can always find resonant modes, satisfying (16 ).

Strictly speaking, the derived conditions do not solve the formulated task yet, because having considered modes of \(H_0\) with indices \(j\) and \(j'\), we have not determined yet the modes \(\phi_i^{(0)}\) and \(\phi_{i'}^{(0)}\), which will interact resonantly: formally they can be different from the modes \(\phi_j^{(0)}\) and \(\phi_{j'}^{(0)}\), which give origin to Gershgorin intervals. To complete this last step, we recall that in the limit \(\alpha\to 0\) the Gershgorin intervals collapse to points \(\beta_i^{(0)}\) and remain non-overlapped for sufficiently small \(\alpha\) [see examples in Figs. 3 (c) and (d)]. Due to continuous dependence of \(\beta_i^{({\boldsymbol{\alpha}})}\) on \({\boldsymbol{\alpha}}\), it follows that \(\beta_i^{({\boldsymbol{\alpha}})}\) belongs to the Gershgorin interval originated at \(\beta_i^{(0)}\) even when the intervals start to overlap. Thus, for sufficiently small \(\alpha\), one has \(i=j\) and \(i=j'\). In Fig. 3 (c) and (d), we illustrate quasi-resonances and Gershgorin intervals for several modes.

The condition (20 ) invoking the input modes of the Hamiltonian \(H_0\), can be viewed as a quasi-resonance in the “energy” space, representing the proximity of the propagation constants. Conditions (15 ) and (20 ) verified simultaneously constitute the selection rule determining the modes of \(H_0\) in the coordinate and energy spaces, that may interact resonantly with a given input mode. In other words, these conditions determine modes that can be involved in the energy exchange upon propagation. Each of those conditions alone is not sufficient for efficient coupling between modes.

Finally, we note that although in the 2D case the exact resonance at the level crossing \(\beta_{i}^{({\boldsymbol{\alpha}})}=\beta_{i'}^{({\boldsymbol{\alpha}})}\) is possible [examples can be found in Figs. 1 (a) and (b)], in a generic situation it will not result in resonant interactions of the modes, because of negligible hopping integrals \(\langle\phi_j^{(0)},\phi_{j'}^{(\alpha)}\rangle\). In a general situation, the quasi-resonance occurs at avoided crossings [like the ones shown in Fig. 1 (a) and (b)].

3.3 Examples of BLZ oscillations↩︎

3.3.1 Two-mode dynamics↩︎

Now we turn to the examples of different propagation regimes and BLZ oscillations. For the illustration purposes we choose the input mode \(\phi_{124}^{(0)}\) (the mode labeled as point 1 in Fig. 2) for the lattice with \(\gamma=\pi/3\). At sufficiently small \(\alpha\lesssim 0.001\) one cannot encounter other modes satisfying both the proximity condition (15 ) [Fig. 3 (a)] and quasi-resonance condition (20 ) [Fig. 3 (c)]. For this region of parameters, in spite of the change of the propagation constant (blue line in Fig. 2) due to the gradient, the mode exhibits no significant dynamics beyond weak amplitude fluctuations (not shown here). With further increase of gradient \(\alpha\) both spatial and quasi-resonant conditions for efficient arise between modes 124 and 129, as evidenced by the avoided level crossing, see modes corresponding to points 2 and 5 in Fig. 2 (a) indicated by the vertical dashed line in Fig. 3 (a) and (c). Now, large amplitude (in the sense of amount of the transferred power) BLZ oscillations involving these two modes occur. The results of the direct numerical simulations of these oscillations are shown in Fig. 4 (a). The shown oscillations persist over long propagation distances. This extreme robustness of the BLZ oscillations arises from suppressed diffraction in the incommensurate moiré lattices (observed earlier in different experimental settings [58], [59]). The period of the oscillations in theory is determined by the formula \(Z_{ij}=2\pi/|\beta_i^{({\boldsymbol{\alpha}})}-\beta_j^{({\boldsymbol{\alpha}})}|\). Estimating this quantity for \(\beta_{124}^{({\boldsymbol{\alpha}})}\) and \(\beta_{129}^{({\boldsymbol{\alpha}})}\) at \(\alpha=0.003\) (see Fig. 2 (a)) one obtains \(Z\approx 761.8\) what matches well the period observed in direct simulations [Fig. 4 (a)]. In the dynamical figures, the prediction of the analytical theory is shown as dashed lines and closely matches numerical results.

A remarkable feature of BLZ oscillations in moiré lattices is that increasing the gradient does not necessarily intensify power exchange between modes; on the contrary, it can have the opposite effect, weakening the coupling and thereby damping the beam’s oscillations. Indeed, our numerical simulations (not shown here) confirm that a further increase of the gradient results in much weaker power exchange between modes 124 and 129. Note that at \(\alpha=0.01\) the Gershgorin radii, although still overlapping [Fig. 3 (c)], are relatively large such that the quasi-resonant condition (16 ) is no longer satisfied. In this regime, the only noticeable impact of increasing the gradient is a higher frequency of the low‑amplitude oscillations of the mode intensity, instead of the enhancement of power transfer.

As it was mentioned above, the direction of the gradient plays a crucial role. Fig. 4 (b) shows the dynamics for the same initial state \(\phi_{124}^{(0)}\) in a system with the same gradient \(\alpha = 0.003\), but with angle \(\gamma = -\pi/5\). As evident from the figure, the energy is now transferred to a completely different mode, specifically, to the one labeled 116, demonstrating the strong dependence of the dynamics on the gradient direction. Note that while all modes within the effective area were accounted, in the upper panels of Fig. 4 (a) and Fig. 4 (b) one can distinguish contributions from only three modes (with two modes dominating the dynamics in each case): the input mode 124 and resonantly interacting mode 129 at \(\gamma=\pi/3\) [panel (a)] and 109 at \(\gamma=-\pi/5\) [panel (b)]. This reveals the effect of the direction of the gradient on the BLZ oscillations, which is shown by the arrows in the leftmost snapshots of the dynamics. Remarkably, one observes that the positions of the modes in resonance do not have any visual connection with the direction of the gradient, which is in full agreement with the theoretical predictions since the selection rules do not depend on the angle \(\gamma\) in a direct way. It is important to note that changing the orientation of the gradient \(\gamma\) can not only alter the resonant state involved in power exchange, but in certain directions, it can completely suppress this exchange, causing the initial state to propagate with minimal oscillations. For example, at a fixed \(\alpha=0.003\), the gradient angle can be tuned such that only weak amplitude oscillations occur, without any significant mode coupling.

3.3.2 Oscillations involving several modes↩︎

An even more intricate dependence of the dynamics on the gradient \(\alpha\) emerges when different modes satisfy the selection rule at different values of \(\alpha\). To illustrate this, consider the input mode \(\phi_{137}^{(0)}\) propagating in the lattice with angle of gradient \(\gamma=\pi/3\). The change of its propagation constant upon change of the gradient is shown by the green curve in Fig. 2 (a), while the corresponding coupling coefficients \(S_{i,137}^{({\boldsymbol{\alpha}})}\) and Gershgorin intervals are presented in Fig. 3 (b) and (d), respectively. As these plots indicate, the selection rule is fulfilled for mode 135 around \(\alpha \approx 0.001\). With increasing \(\alpha\), this resonance condition is no longer met for mode 135 but becomes valid for mode 127 near \(\alpha \approx 0.005\). Beyond this, up to \(\alpha = 0.01\), no modes satisfy the selection criteria. This behavior is confirmed by direct propagation simulations of mode \(\phi_{137}^{(0)}\) under different gradients. At \(\alpha = 0.001\), energy transfer between modes 137 and 135 is observed [Fig. 5 (a)]. As the gradient increases to \(\alpha = 0.005\), the interaction with mode 135 is suppressed, and coupling with mode 127 appears [Fig. 5 (b)]. Once again, increasing the gradient does not necessarily enhance mode coupling; in fact, for \(\alpha = 0.01\), energy exchange is almost completely suppressed. This is evident from the snapshots shown in Fig. 5 (c), where the amplitude profile remains nearly unchanged during the propagation.

3.3.3 Multimode oscillations for a different sublattice orientation↩︎

The described physics of BLZ oscillations is confirmed by numerical simulations performed for other sublattice orientations. Figure 6 illustrates the evolution of the mode weights \(|c_n|^2\) of all linear modes at \(|{\boldsymbol{\alpha}}| = 0\) for the system with \(\gamma = \pi/3\) and \(|{\boldsymbol{\alpha}}| = 0.003\), obtained for the angle between sublattices \(\theta = \pi/5\). The initial state is the mode \(n = 101\). For this mode, the spatial proximity condition \(|S_{i,101}^{({\boldsymbol{\alpha}})}| > \Delta\) is satisfied only for \(i = 102\), \(105\), and \(106\), while overlap with other modes is negligibly small. For these three states, the spectral proximity condition is fulfilled as well. Consistently, Fig. 6 demonstrates pronounced power exchange among modes \(101\), \(102\), and \(105\), with weaker coupling to mode \(106\), whereas interactions with all other modes remain suppressed. The lower panels of Fig. 6 display the field modulus distributions at different propagation distances \(z\), marked by the vertical dashed lines in the top panel. These snapshots indicate that the power transfer occurs predominantly along a single spatial line, which is noticeably misaligned with the gradient direction \(\gamma\). Finally, the analytical predictions for weights \(|c_n|^2\) dynamics (dashed curves) are in excellent agreement with the numerical simulations.

3.3.4 Oscillatory dynamics in a perturbed system↩︎

To examine the structural robustness of the BLZ dynamics we simulated propagation of a beam in a moiré lattice where the waveguide depths of the sublattices were randomly distributed within the interval \([p_j(1-\delta_p), p_j(1+\delta_p)]\), \(j=1,2\), and the waveguide positions in both the \(x\) and \(y\) directions were shifted by random values in the range \([-\delta_r, \delta_r]\). We used parameters from Fig. 4 (a) (i.e., \(\gamma = \pi/3\), \(|{\boldsymbol{\alpha}}| = 0.003\), and \(\theta = \pi/3\)) simulated now for \(\delta_p = 0.02\) and \(\delta_r = 0.01\). We used the same input mode, \(n = 124\), of the unperturbed lattice without gradient, \(|{\boldsymbol{\alpha}}|=0\) as in Fig. 4 (a). In Fig. 7 we show the dynamics of the beam projections \(|c_i|^2\) on the eigenmodes of the unperturbed lattice, with different line styles representing distinct disorder realizations. Although the oscillation period and amplitude are moderately modified compared to the unperturbed case, the mode composition remains unchanged, and oscillations persist even when the excitation is not an eigenmode of the perturbed system. We emphasize that for each realization of disorder the eigenmodes of the perturbed system can be determined, by analogy with the unperturbed moiré lattice, and the spectrum remains discrete. Thus, our theoretical framework can then be used to predict the corresponding mode composition and fully describe the BLZ dynamics. The selection rule still determines the set of excited modes, however, a thorough analytical description requires a full numerical analysis of eigenmodes of the randomized lattice. It should be mentioned that for \(|{\boldsymbol{\alpha}}| = 0.003\), the dominant modes exhibit dynamics comparable to those of the unperturbed system, whereas for larger values of \(|{\boldsymbol{\alpha}}|\), the effect of disorder becomes more pronounced, occasionally inducing power transfer to additional modes not involved in the unperturbed dynamics.

4 Nonlinear BLZ oscillations↩︎

Now we turn to the effect of weak nonlinearity of the BLZ oscillations governed by the model \[\begin{align} \label{main-nonlin} i\frac{\partial \Psi}{\partial z}=H_{{\boldsymbol{\alpha}}}\Psi+g|\Psi|^2\Psi , \end{align}\tag{21}\] where \(g=1\) and \(g=-1\) describe defocusing and focusing nonlinearities, respectively, while \(H_{{\boldsymbol{\alpha}}}\) is the Hamiltonian introduced in (2 ). Since our consideration is limited to the localized modes, neglecting excitation of extended states, in the nonlinear case, one can still search for a solution of (21 ) in the form of the expansion over the basis \(\phi_i^{(0)}\) [see (8 )]. For sufficiently small input intensities, where only localized modes are excited at \(z=0\), this condition can be expressed as \(U=\sum_{j=1}^N|c_j|^2\ll 1\). In the leading order, substituting the ansatz into Eq. (21 ) yields the following system for the modal amplitudes: \[\begin{align} \label{dc95095nonlin} i\frac{d{c}_j}{dz}=-\beta_j^{(0)}c_j-\sum_{k=1}^N({\boldsymbol{\alpha}}\cdot{\boldsymbol{R}}^{(0)})_{jk}{c}_k+\chi_j|{c}_j|^2{c}_j , \end{align}\tag{22}\] where \(\chi_j=g\langle[\phi_j^{(0)}]^2,[\phi_j^{(0)}]^2\rangle\) is the effective nonlinearity (to shorten the notations we do not use the upper index \(0\) in the coefficients \(c_j(z)\), since only the basis of \(H_0\) is considered below). In the system (22 ) we neglected all the terms proportional to \(\langle \phi_{i_1}^{(0)}\phi_{i_2}^{(0)},\phi_{i_3}^{(0)},\phi_{i_4}^{(0)}\rangle\), except the terms where all coefficients are equal \(i_1=\cdots=i_4\), i.e., \(\chi_j\). This is justified by the fact that such overlap integrals are significant only for modes that are spatially close to each other (see, e.g., the 1D numerical examples in [71]). Meanwhile, due to the selection rule described above, the modes resonantly interacting at nonzero gradient are typically located relatively far apart at \(\alpha=0\) [see Fig. 1 (b), Fig. 2 (d), and Fig. 3 (a),(b)]. We can thus conclude that the nonlinear term in Eq. (21 ) is expressed predominantly through the self‑action of the modes. The model describes the dynamics well for small input powers, specified above (not shown here).

As an example of light propagation in a nonlinear medium with a refractive index gradient at larger nonlinearities, we consider the initial mode \(\phi_{137}^{(0)}\) in a system with \(\alpha = 0.005\) and \(\gamma = \pi/3\). The same configuration was used for the linear dynamics in Fig. 5 (b). Nonlinear evolution for different input powers \(U=|c_j|^2\) is shown in Fig. 8, with the left column [Figs. 8 (a)–(c)] corresponding to a focusing medium and the right column [Figs. 8 (d)–(f)] to a defocusing one. In the linear case, power transfer occurred between modes 137 and 127. However, under weak nonlinearity (e.g., \(U=0.15\) for focusing, \(U=0.20\) for defocusing), this transfer is suppressed, as evidenced in Figs. 8 (a) and (d): only minor variations in modal weights occur, and the field profiles remain nearly unchanged. Notably, model (22 ) predicts suppression of oscillations for these modes 137 and 127 as power increases; however, at this relatively high power level, it no longer accurately captures the amplitude and period of these low-amplitude dynamics. For lower initial powers \(U\) (not shown here), the predictions of this model remain reasonably accurate. As the input power increases, nonlinear effects can bring new modes into resonance. For instance, at \(U=0.60\) in the focusing regime [Fig. 8 (b)], coupling emerges primarily with mode 114, and to a lesser extent with other modes, resulting in irregular dynamics. In the defocusing case at the same power [Fig. 8 (e)], modes 152 and 176 are excited instead, both having lower propagation constants than the initial mode. At even higher input powers [Figs. 8 (c) and (f)], the system exhibits increasingly complex dynamics with simultaneous coupling to multiple modes, resulting in broad energy redistribution across the spectrum.

Thus, nonlinearities can not only suppress BLZ oscillations, but also trigger new mode couplings. Moreover, BLZ oscillations exhibit remarkable robustness against weak nonlinear effects: the few-mode character of the dynamics is preserved, although the period and amplitude of oscillations can be noticeably altered. Notably, despite the broader range of instabilities typical in two-dimensional systems, no such instabilities are observed in our case.

5 Discussion↩︎

We now address several issues related to the experimental feasibility of observing two-dimensional BLZ oscillations. As mentioned in the Introduction, the system described by Eqs. (1 )-(3 ) can be realized in photorefractive crystals using method of optical induction [54] that also allows creation of moiré potentials at will [58]. In such systems, transverse refractive index gradients can be realized by illumination of the sample with broad beam with linearly varying intensity. Similar configurations can also be created using highly developed technology of fs-laser inscription of waveguide arrays [72] that was already successfully applied to creation of moire arrays [61]. In these systems, the transverse refractive index gradient can be emulated by bending of waveguides along the parabolic trajectory with distance \(z\) [18], [73].

In the last case, the dimensionless parameters used in the model are as follows. The transverse coordinates are defined by \(\boldsymbol{r}=\tilde{\boldsymbol{r}}/r_0\), where \(r_0=10\) \(\mu\text{m}\) (hereafter tilde indicates the respective variables in the physical units) is the characteristic transverse scale and the dimensionless propagation distance \(z=\tilde{z}/L_d\) is normalized to the diffraction length \(L_d=kr_0^2\approx1.14\) \(\text{mm}\), where \(k=2\pi n/\lambda\) is the wavenumber at typical wavelength \(\lambda=800\) \(\text{nm}\), \(n\approx1.45\) is the unperturbed refractive index of transparent dielectric material, where moiré lattice is inscribed (for example, fused silica). The dimensionless field amplitude is related to the real field amplitude \(\mathcal{E}\) via \(\Psi=\left(k^2r_0^2n_2/n\right)^{1/2}\mathcal{E}\), where the nonlinear refractive index (typical for fused silica) \(n_2\approx2.7\times10^{-20}\) \(\text{m}^2/\text{W}\). The dimensionless depth of sublattices forming moire lattice is given by \(p_{1,2}=k^2r_0^2\delta n_{1,2}/n\), where \(\delta n_{1,2}\) is the refractive index modulation depth in each sublattice. The representative values \(p_{1,2}=4\) used here correspond to a refractive index modulation depth of \(\delta n_{1,2}\approx4.5\times10^{-4}\). The dimensionless width of individual waveguides \(w=0.5\) corresponds to \(5\) \(\mu\textrm{m}\), while the period \(d=2.5\) of both sublattices corresponds \(25\) \(\mu\textrm{m}\). For experiments with low-power beams the samples with lengths up to \(20-30\) cm can be fabricated, that correspond to dimensionless propagation distances of \(z\approx180-270\).

It is relevant to notice that modern experimental techniques in optical waveguiding systems allow registration of not only output intensity distributions, but also of the intensity distributions at intermediate cross-sections inside the sample or even the entire light propagation dynamics inside the sample. Moreover, developed holographic techniques and spatial light modulators (SLMs) allow to create practically any desirable field distribution at the input face of the sample that can be focused into selected waveguides/location of the structure to excite a predetermined set of modes and to observe theoretically predicted dynamics.

Regarding the stability of BLZ oscillations, which is crucial for their experimental observation, deviations of the initial beam from the exact eigenmodes do not lead to instabilities, even in the weakly nonlinear regime, as demonstrated by the simulations above (obviously, no instabilities can arise in the linear case when modes do not interact with each other). Concerning structural stability, i.e., deviations of the potential from the exact form used in the numerical simulations (illustrated in Sec. 3.3.4), no instability is observed due to the inherent nature of quasi-periodic systems and their discrete spectrum (for a finite specimen). Indeed, any small perturbation of \(V(\boldsymbol{r})\) can only shift the discrete spectrum slightly and cause weak changes in the position and shape of the mode. However, such perturbations do not affect the selection rule itself. Thus, one may observe quantitative variations in the period and amplitude of the oscillations, but not their destruction.

Another important consideration is the preparation of the input signal, which depends on the experimental configuration. In the case of optical beams, the primary focus of this study, devices such as spatial light modulators enable the generation of nearly arbitrary input profiles. The programmability of a phase mask allows precise control over the spacing between excitation spots, their size, individual intensities, and relative phases. Although various methods exist to create specific input conditions, one of the simplest approaches is a single-site excitation. To assess the robustness and experimental feasibility of BLZ oscillations, we consider the parameters \(\gamma = \pi/3\), \(|\alpha| = 0.005\), and \(\theta = \pi/3\), with a single-site excitation at the waveguide corresponding to the maximum intensity of the mode \(n = 137\). This mode served as the initial condition for the dynamics depicted in Fig. 5 (b). The results of the single-site excitation are presented in Fig. 9. Although the oscillation amplitudes of the mode weights \(|c_n|^2\) are slightly reduced due to the radiation, both the oscillation period and the modal composition remain almost unchanged relative to the case of mode-driven excitation. These results demonstrate that BLZ oscillations can be reliably excited in an optical system using a single-beam input.

6 Conclusion↩︎

To conclude, we have developed a theory of two-dimensional BLZ oscillations for localized modes in aperiodic lattices, illustrating the results with a light beam propagating in the incommensurate moiré lattice with a transverse refractive index gradient. Since modes with propagation constants (energies) above (below) the mobility edge are localized in the coordinate space but possess different propagation constants (energies), the linear dynamics is governed by the simultaneous spatial tunneling and Landau-Zener tunneling in the energy space. This dynamics is highly sensitive to the initial conditions and may involve several resonantly coupled modes determined by the selection rule. In the 2D setting, in contrast to their 1D counterpart, the selection rule depends on the angle of the gradient, which can result in resonantly coupled modes being misaligned with the direction of the linear gradient. The dependence on the gradient arises from the spatial localization of the modes, which further explains why quasi-resonances at avoided crossings are more relevant for the dynamics than the level crossings that are typical in 2D settings. Furthermore, BLZ oscillations may not appear at small gradients or become suppressed with the increase of the gradient, if no modes satisfying the selection rule are found. Weak nonlinearity, whether attractive or repulsive, preserves the few-mode character of the dynamics, although it can significantly modify the period and amplitude of the oscillations.

The reported theory allows direct generalization, including different types of aperiodic potentials, like quasicrystals, or described by two-dimensional almost-periodic functions, all of which have been successfully realized in laboratory settings.

Although our analysis focuses on optical systems, the underlying theory is directly applicable to other experimentally available settings such as BECs [74], where optical lattices in general are highly tunable [75], [76], or acoustic waves  [77] in systems governed by moiré, or more generally, quasi-periodic potentials.

Acknowledgements↩︎

SKI has received funding from the European Union through the Program Fondo Social Europeo Plus 2021-2027(FSE+) of the Valencian Community (Generalitat Valenciana CIAPOS/2023/329). YVK was supported by the Russian Science Foundation (grant 24-12-00167) and partially by the research project FFUU-2024-0003 of the Institute of Spectroscopy of the Russian Academy of Sciences. VVK was supported by the Fundação para a Ciência e Tecnologia under the projects 2023.13176.PEX (DOI
https://doi.org/10.54499/2023.13176.PEX) and by national funds, under the Unit CFTC - Centro de Física Teórica e Computacional, reference UID/00618/2023, financing period 2025-2029.

Data availability statement↩︎

All data that support the findings of this study are included within the article or available upon reasonable request from the authors.

Conflict of Interest↩︎

The authors declare no conflict of interest

References↩︎