Two-dimensional flat-band solitons in superhoneycomb lattices

1 Introduction

Existence of a flat band or of several flat bands in the spectrum of a linear Hamiltonian significantly changes the properties of a system. This is well understood theoretically and verified experimentally in diverse areas of physics ranging from solid-state physics, where flat bands appear in bi-layer graphene, to atomic and optical mono-layered systems, where flat bands arise in spectra of some lattice potentials [1]–[5]. The majority of studies of optical flat-band systems was concentrated on discrete settings, which represent tight-binding limits of the respective continuous models with periodic refractive index landscapes. From the theoretical point of view discrete lattices benefit from the existence of exact flat bands and from the possibility of algorithmic design of structures that feature such bands [6]. From the experimental point of view, discrete models adequately describe light propagation dynamics in arrays of sufficiently deep waveguides [7]–[17] representing also one of the most powerful platforms for the exploration of self-action of light in periodic environment, but produce considerable deviations from observed dynamics in lattices with smooth refractive index landscapes or in arrays with relatively small refractive index contrast. Very recently, photonic structures with multiple flat bands on the basis of microwave resonators were reported [18].

Exactly flat bands in discrete systems have opened opportunities for construction of new types of self-sustained states. Compactons and solitons were shown to exist in both conservative and dissipative discrete arrays [19]. Two-dimensional thresholdless solitons bifurcating from linear compact modes have been found in kagome arrays [20]. Flat-band solitons and nonlinear localized flat-band modes were also reported in sawtooth-like lattices [21], diamond chain lattices [22], [23], Stub lattices [24], [25], Lieb lattices [26]–[29], and octagonal-diamond lattices [30]. A summary of recent theoretical and experimental advances in the area of multidimensional localized structures in optical media, including the formation of different types of lattice solitons in discrete and continuous physical models can be found in Ref. [31].

However, in genuinely continuous systems, whose accurate description requires going beyond tight-binding approximations, exactly flat bands do not exist. Instead, one can realize continuous systems with nearly flat bands, whose width is sufficiently small, but not exactly zero. We emphasize that here the flatness of the band is implied in all directions of the reduced Brillouin zone (BZ). This opens the important and so far unaddressed question about the diversity, existence, and bifurcations of self-sustained nonlinear states in such two-dimensional (2D) continuous waveguiding systems.

Some preliminary conclusions that the behaviour of nonlinear states in such systems can be very unusual and can sharply contrast with behaviour of discrete systems, can be drawn on the basis of recent studies of Wannier solitons in 1D flat-band system describing spin–orbit coupled Bose–Einstein condensates [32] and studies of light propagation in photonic moiré lattices [33], where nearly all higher bands of the optical potential are flat, allowing excitation of 2D solitons practically with zero power threshold [34].

In addition to obvious limitations of the tight-binding (discrete) models that neglect other (non-flat) bands which are necessarily excited in nonlinear systems, there is also a conceptual difference consisting in certain ambiguity of such models. Indeed, deriving a discrete model from a continuous one requires definition of a proper basis. Even restricting the consideration to the most natural choice of a basis of Wannier functions (WFs), the latter are not uniquely defined, what is particularly relevant for two- and three-dimensional systems [35]. The coupling coefficients of the respective one-band discrete approximation are determined by the hopping of the chosen basis functions, and thus are not uniquely defined for a given continuous model, as well. On the other hand, different continuous models may result in the same equation of the tight-binding approximation. This does not allow one to make one-to-one correspondence between the original nonlinear continuous model and its tight-binding approximation, what makes direct study of the continuous models particularly relevant.

In this work we use continuous superhoneycomb lattices that possess a nearly flat band in their linear spectrum, to study the emergence and stability properties of solitons of very different types. Our choice of the superhoneycomb lattice is dictated by its properties which are known due to previous studies [16], [36], [37] allowing us to straightforwardly find the parameter enabling a nearly flat band. By changing the depths of some waveguides it is possible to tune the band structure such that forbidden gap appears either above or below the flat band, as shown in Figure 1 in the following text. This allows to obtain flatband solitons in this continuous system both in the medium with focusing (above the flat band) and defocusing (below the flat band) nonlinearity. In addition the degree of the band flatness that one can achieve for this type of the lattice exceeds band flatness for more traditional Lieb and kagome lattices with the same depth and waveguide spacing. We report families of flat-band soliton and vortex-soliton solutions and show that in continuous lattices, unlike in their discrete counterparts, they are not thresholdless anymore. The thresholds for solitons of very different types, such as fundamental and vortical ones, may be very small and comparable.

Figure 1:

Superhoneycomb lattice and its band structure. (a) Superhoneycomb lattice, where all waveguides have the same depth p = 8.0. The basis vectors of the Bravais lattice v _1,2, the distance d between two sites, and labels of sites within primitive lattice cell are shown. (b) Superhoneycomb lattice with depths of sites A and B p n A = p n B = 8.0 and depths of sites C, D, E p n C = p n D = p n E = 7.8 . We use the notation n _A for positions of A sites, and similar notations n _D , …, n _E for positions of other sites. (c) Superhoneycomb lattice with p n A = p n B = 7.8 and p n C = p n D = p n E = 8.0 . (d–f) Band structures of lattices depicted in panels (a–c), respectively. The flat band present in the spectrum is indicated by red line in (d) and is preserved even when the gap appears in (e, f). Other parameters are d = 2 and w = 0.5.

2 Flat bands in superhoneycomb lattices

The dimensionless envelope Ψ of a paraxial light beam propagating along the z-direction in a 2D optical lattice, described by the function R ( r ) , where r = (x, y), is governed by the nonlinear Schrödinger (NLS) equation

with the linear Hamiltonian

where ∇ = (∂/∂x, ∂/∂y) and σ = ±1 corresponds to the Kerr medium either with a focusing (σ = +1) or defocusing (σ = −1) nonlinearity. In this work we focus on stationary states of the problem described by Equations (1) and (2). Such states have the form Ψ = ψ( r ) exp(ibz), where b is the propagation constant, and the function ψ( r ) describing the transverse profile of the mode, solves the stationary nonlinear problem

While the theory reported below is not restricted to any specific type of optical potential (except the requirement to possess a nearly flat band), in all numerical simulations we employ the superhoneycomb lattice described by

Here n = [n ₁, n ₂] with integers n _1,2 defines the positions of the lattice “sites”: v _n = n ₁ v ₁ + n ₂ v ₂, where v _1,2 are the basis vectors of the Bravais lattice, p _n stands for the depth of each waveguide, and w is the waveguide width. The profile of the lattice for the particular case, when all p _n are equal, p _n = p, is shown in Figure 1(a).

The spectrum of the underlying linear Hamiltonian in Equation (2), is determined by the 2D linear problem

for 2D linear Bloch modes φ _{ν k} ( r ) = e ^{i k ⋅ r} u _{ν k} ( r ), where u _{ν k} ( r ) = u _νk ( r + v _n ), v _n is a lattice vector, ν is the index of the band, and k is the Bloch vector in the first BZ. In Figure 1(d) we show the spectrum of the problem corresponding to the lattice in Figure 1(a) (it was obtained using plane-wave expansion method and details on this method can be find in the Section 7), where a flat band is present, as indicated by the red line. The width of the flat band with index ν given by

is not exactly zero and can be used as a parameter characterizing band flatness. The ideally flat band would have Δ _ν = 0, while for the band shown in Figure 1(d) one has Δ _ν ≈ 2.77 × 10⁻⁵. Notice that Δ _ν increases with decrease of the lattice depth p. We stress that we consider here shallow periodic photonic structures defined by small refractive index modulations. For example, the dimensionless lattice depth of p = k 2 r 0 2 δ n / n ∼ 8 used here corresponds to the refractive index contrast δn ∼ 8 × 10⁻⁴, where k = 2πn/λ is the wavenumber at the wavelength λ = 800 nm, n = 1.45 is the unperturbed refractive index of the material (fused silica, for example), r ₀ = 10 μm is the characteristic transverse scale to which coordinates x, y are normalized. Thus the Equation (1) is ideally suited for accurate description of the paraxial light propagation in such structures [33], [38]–[44].

The existence of gap solitons hinges on the requirement that the flat band is isolated from the rest of the spectrum by a gap, positioned either below or above the band. Examples of such situations are known [36], [37], [45]. To make possible the existence of gap solitons in our case, we notice that each primitive cell of the superhoneycomb lattice has five waveguides, denoted by letters A–D in Figure 1(a), and judiciously detune the depths of two groups of waveguides by making p _n different. Thus, when the waveguides A and B have larger depths than other waveguides, as shown in Figure 1(b), a gap appears above the flat band, while the width of this band practically does not increase and is equal to Δ _ν ≈ 3.54 × 10⁻⁵ [see red line in Figure 1(e)]. Meantime no gap opens below this flat band that touches the next band ν + 1 at the center of the BZ. On the other hand, if the depths of waveguides A and B are reduced in comparison with depths of other three waveguides in the unit cell [Figure 1(c)], then a relatively wide gap opens below the flat band, that in this case keeps its width at Δ _ν ≈ 2.94 × 10⁻⁵, as shown by the red line in Figure 1(f).

3 Two-dimensional Wannier solitons and vortex solitons: numerical study

We are interested in properties of gap solitons, i.e., square integrable solutions, whose propagation constants b belong to a gap adjacent to the flat band. We consider here only static solitons that do not move across the lattice, i.e., their spatial group velocity is zero, while field modulus distribution does not change with propagation distance z. Families of such solutions can be parameterized by the total power P ( b ) = ∫ r ∈ R 2 | Ψ | 2 d r carried by the mode. Let us introduce the notation b _co = max _{k ∈BZ} b _ν ( k ) (b _co = min _{k ∈BZ} b _ν ( k )) in the case of a gap above (below) the flat band corresponding to the case illustrated in Figure 1(d) and (e). In 2D lattices with non-flat bands, 2D solitons carrying a finite total power cannot bifurcate from linear Bloch states in a sense that P(b) does not tend to zero in the limit b → b _co [46] (see also numerical studies in Refs. [47], [48]), and hence, there exists a minimal threshold power P _th which is achieved at a certain b = b _th: P(b) ≥ P _th = P(b _th) (such constraint does not exist in discrete models). The amplitude of a nonlinear solution still vanishes at b → b _co. However, when a solution forms in a gap above or below the flat band, one observes a different picture, and the dependence P(b) acquires qualitatively new properties. We illustrate this first using direct numerical calculation of the families of the flat-band solitons. Such solutions can be found using standard Newton method [49]. To check their stability, we introduce a small-scale perturbation into soliton profiles and model their long-distance propagation according to Equation (1). If the profile of soliton remains unchanged upon propagation, we conclude that it is stable, otherwise it is considered unstable. We also would like to stress that the presence of optical lattice potential suppresses critical collapse that otherwise could occur for solitons in uniform focusing nonlinear medium.

3.1 Focusing case

We start with the focusing case σ = +1 when gap solitons exist above the flat band in the lattice depicted in Figure 1(b). Figure 2(a) shows three different families of solitons with the vertical dashed line representing the location of the flat band. The lowest (green) family corresponds to the fundamental solitons. In a flat-band system such solitons remain surprisingly well localized even very close to the band, i.e., even at relatively small values of detuning δb = b − b _co from the band [see soliton 1 shown in Figure 2(b)]. When δb increases, the localization of soliton moderately increases too [see soliton 2 shown in Figure 2(b)]. Fundamental solitons are stable in the largest part of the gap, but become unstable close to its upper edge. One can use such fundamental solitons for construction of a high-order solitons that can be viewed as coupled fundamental states placed one next to the other. For example, we can use two fundamental solitons ψ ₁ and ψ ₂ to construct an initial guess via ψ ₁(y + δy) − ψ ₂(y − δy) with δy being a shift in y. Using this guess one can obtain dipole soliton via Newton iterations [see soliton 3 in Figure 2(b)]. Dipole solitons also demonstrate stability in the largest part of the gap [see red solid line in Figure 2(a)], becoming unstable only for large enough values of detuning δb [see red dashed line in Figure 2(a)]. We also found a family of vortex solitons shown by the blue line in Figure 2(a). The stability region of vortex solitons, however, is much smaller as compared to that for fundamental or dipole solitons. Representative field modulus and phase distribution in stable vortex soliton 5 taken very close to the flat band is shown in Figure 2(b). The solitons corresponding to dot 4 and dot 6 in Figure 2(b) are more localized in comparison with these corresponding to dot 3 and dot 5.

Figure 2:

The families of solitons and examples of profiles in focusing medium. (a) Flat-band soliton families for focusing nonlinearity. The green, red and blue curves show families of fundamental, dipole, and vortex solitons, respectively. Solid and dashed curves represent stable and unstable solitons. The vertical dashed line at b = 2.4283 indicates the location of the flat band, whose width cannot be discerned on the scale of the figure. The families are terminated at the propagation constant b = 2.429. (b) Field modulus distributions |ψ| in selected flat-band solitons numbered 1–6. Panels in (b) are shown within −15 ≤ x, y ≤ 15 window.

The families of all such solitons share two common features. First, although gap solitons cannot bifurcate from linear Bloch modes (in the above mentioned sense) and exist above a certain power threshold, P _th, this threshold becomes extremely small in flat-band system. We successfully traced the families depicted in Figure 2(a) practically up to the edge of the flat band. Further approaching b _co became technically infeasible, even with very small increments of the propagation constant. Notice that even in such close proximity of the flat band, solitons remain well localized, in particular in comparison with solitons in lattices with dispersive bands that expand dramatically when their propagation constant approaches the edge of the band. Second, in a certain range of propagation constants near the flat band, the power P becomes a linear function of detuning δb, as discussed below (refer to Figure 4).

3.2 Defocusing case

From the physical point of view, a spatial soliton forms due to a delicate balance between diffraction and nonlinearity. Since in the leading order the diffraction is determined by the curvature of the allowed band, it is natural to expect that vanishing curvature of a flat band should result in similarity of properties of bright and dark solitons forming in the gap above and below such a band. This is fully confirmed by Figure 3 that shows soliton families [Figure 3(a)] and shapes [Figure 3(b)] of representative flat-band solitons in defocusing medium (σ = −1), for the lattice depicted in Figure 1(c). Now solitons exist in the region δb < 0, but except for this the obtained P(b) dependencies in Figure 3(a) are remarkably similar to dependencies obtained for focusing case in Figure 2(a) (we use the same color coding for fundamental, dipole, and vortex solitons). Stability properties of soliton families are also similar, i.e., fundamental and dipole solitons in defocusing medium are stable in the largest part of the gap adjacent to the flat band, while vortex solitons are stable only a narrow region near the flat band. Just as in focusing medium analysis of P(b) dependencies near the flat band reveals the presence of very small thresholds (comparable for solitons of different types) also in defocusing medium.

Figure 3:

The families of solitons and examples of profiles in defocusing medium. Color coding of curves is the same as in Figure 2. The vertical dashed line at b = 2.5494 indicates the location of the flat band, whose width cannot be discerned on the scale of the figure. The cut-off propagation constants of the red, green, and blue lines are b = 2.5488, b = 2.549, and b = 2.5486, respectively.

It should be pointed out that while the dependencies P(b) for flat-band solitons are qualitatively similar in focusing and defocusing media (except for the sign of the derivative dP/db for the corresponding families), and the field modulus distributions for solitons shown in Figures 2(b) and 3(b) are also similar, the internal structure of solitons in focusing and defocusing media is different. This is particularly well visible from comparison of solutions 2 and 4 in Figures 2 and 3. One can see that in focusing medium the spots nearest to the central maximum of soliton are out-of-phase with central spot, while in defocusing medium they are in-phase with central spot. Therefore, the phrase structure in tails of flat-band solitons in focusing and defocusing media is clearly different. One can also see such differences from phase and field modulus distributions of vortex solitons in Figures 2(b) and 3(b).

4 Flat-band solitons: understanding the results

4.2 Dependence of power on propagation constant

For the Wannier soliton (16) one obtains that the total power is |a|². This leads to the relation

(17) P ( a , b ) ⋅ χ ( a ) ≈ σ ⋅ δ b .

A peculiarity of this expression is that its right-hand side does not depend on the soliton amplitude: even though both power P and form-factor χ depend on a, their product for the flat-band soliton is determined only by the detuning δb. Furthermore, the approximate expansion in Equation (9) does not imply bifurcation from a linear mode, leaving freedom in the choice of a WF (the phase θ(a, k ) introduced in Equation (6) remains arbitrary) while the exact nonlinear solution for Equation (1) is well defined. These factors not yet accounted in the expansion, strictly speaking, leave open exact applicability conditions for Equation (17). If however one conjecture that a flat-band soliton is approximated by some of WFs well enough, instead of searching for such specific WF, one can use the numerically found nonlinear solution instead of Equation (16), interpreting the result in Equation (17) as the relation

(18) P s ⋅ χ s ≈ σ α ⋅ δ b ,

where P _s and χ _s are the power and the form-factor of the exact soliton, while α is the correction factor of order one.

We have checked the accuracy of formula in Equation (18) for all soliton families depicted in Figures 2 and 3, and obtained a remarkable agreement illustrated in Figure 4. In all cases the accurate linear fit of the families of fundamental, dipole, and vortex solitons was obtained with a universal factor α = 1.11.

Figure 4:

Comparison of flat-band soliton families obtained numerically (see blue curves corresponding to the numerically obtained soliton families shown in Figures 2 and 3) and theoretically (red lines, based on Equation (18)). Panels (a1–a3) correspond, respectively, to families of fundamental, dipole, and vortex solitons in focusing medium. Panels (b1–b3) correspond to fundamental, dipole, and vortex solitons, but in defocusing nonlinear medium.

5 Flat-band multi-soliton solutions

Since flat-band (alias Wannier) solitons are strongly localized and can be stable in considerable part of the gap, it is natural to use them as building blocks for composition of multi-soliton solutions or even of soliton trains. Some solutions of this type are illustrated in Figure 5. Solitons in Figure 5(a), (b) and (e) are obtained in focusing nonlinear medium, while solitons in Figure 5(c), (d) and (f) were found for defocusing nonlinearity. Using the in-phase superposition two fundamental solitons, i.e., ψ ₁(y + δy) + ψ ₂(y − δy) (as opposed to out-of-phase combination yielding dipole solitons presented in Figures 2 and 3) one obtains the 2-soliton (or even) solution depicted in Figure 5(a). We also found various multi-soliton solutions aligned along the x-axis. Thus, in Figure 5(b), we display a 7-soliton solution, which demonstrates that the flat-band solitons exist in abundant forms. Similar higher-order states exist and can be stable in defocusing medium, as shown in Figure 5(c) and (d). While the 7-solitons in Figure 5(b) and (d) are composed of in-phase states, out-of-phase higher-order multi-soliton solutions also exist. To obtain such multi-soliton solutions, we use an initial guess in the Newton method, combining several simple solitons with appropriately selected phase multipliers. An example of such an out-of-phase solution, involving several elements aligned along the y-axis, is shown in Figure 5(e) and (f). The different phase structures of the soliton tails in focusing and defocusing media are also evident in these plots. The solutions presented in Figure 5 can be considered as lattice soliton trains. It is noteworthy that lattice soliton trains have also been reported in non-flat-band systems [52]. However, the important difference from the states presented here is that intentional truncation or partial removal of the train in our case does not lead to significant shape variations of the remaining train upon propagation, whereas it can be quite destructive and lead to the decay of the train in non-flat-band systems. From an experimental point of view, the excitation of multi-soliton states involves creating the proper input, which can be achieved using interferometric techniques and spatial light modulators. Lattices of this type can be easily inscribed using fs-laser writing techniques in fused silica [38], [39], [40], [41], [53], [54], [55].

Figure 5:

Examples of flat-band solitons of other types, including vertical 2-solitons (a, c), in-phase 7-solitons aligned in x (b, d), and out-of-phase 7-solitons aligned in y (e, f). Solutions (a, b) are obtained in focusing nonlinear medium, while solutions (c, d) are obtained in defocusing nonlinear medium. Solutions in (a, c) are shown within −15 ≤ x, y ≤ 15 window, solutions in (b, d) are shown within −30 ≤ x ≤ 30 and −15 ≤ y ≤ 15 window, while solutions in (e, f) are shown within −15 ≤ x ≤ 15 and −44 ≤ y ≤ 44 window.

6 Conclusions

Summarizing, we have described properties of Wannier solitons that can be excited in two-dimensional continuous photonic lattices with a flat band in linear spectrum. Such flat-band solitons exist in both focusing and defocusing media. Their representative feature is that they remain well localized even in close proximity of the flat band. We have found several families of solitons, including fundamental, dipole, multi-soliton, and vortex ones. All these solutions can be stable in focusing or in defocusing medium despite the fact that they form in finite gap adjacent to the flat band. The analytical theory explaining the dependence of power of such solitons on propagation constant is developed for our continuous model. These results illustrate new nonlinear localization scenarios possible in real-world continuous flat-band periodic lattices with the translational symmetry (e.g., the Lieb lattice [7], [8], [14], [56], the kagome lattice [12], [57]–[59], etc.).

7 Methods

The plane-wave expansion method is used to calculate the band structure shown in Figure 1. The primitive vectors of the superhoneycomb lattice are v 1 = [ 2 3 d , 0 ] T and v 2 = [ 3 d , 3 d ] T with the angle between them θ = π/3, as shown in Figure 1(a). The symbol T means the transpose operation. We then apply the following coordinate transformation in the direct and Fourier spaces:

In new coordinate system the primitive lattice vectors become v 1 ′ = [ 2 3 d , 0 ] T and v 2 ′ = [ 0,2 3 d ] T , and they become orthogonal in this frame. In new coordinate system, the Equation (5) can be rewritten as:

where R ′ ( r ′ ) is the superhoneycomb lattice in the frame r ′ after transformation.

We further expand φ and R ′ into the Fourier series:

where c _m,n and ρ _m,n are the Fourier coefficients with m, n are integers, and K m = m π / 3 d . Plugging expansions (25) into Equation (24), after simple algebraic transformations one obtains a series of linear equations:

Rewriting Equation (26) in matrix form and subsequently diagonalizing it, one obtains the eigenvalues b( k ′) (i.e., the band structure) and the corresponding c _m,n that allow to construct the eigenmodes u of the array using the expressions (25).