Hofstadter butterflies in nonlinear Harper lattices, and their optical realizations

The ubiquitous Hofstadter butterfly describes a variety of systems characterized by incommensurable periodicities, ranging from Bloch electrons in magnetic fields and the quantum Hall effect to cold atoms in optical lattices and more. Here, we introduce nonlinearity into the underlying (Harper) model and study the nonlinear spectra and the corresponding extended eigenmodes of nonlinear quasiperiodic systems. We show that the spectra of the nonlinear eigenmodes form deformed versions of the Hofstadter butterfly and demonstrate that the modes can be classified into two families: nonlinear modes that are a ‘continuation’ of the linear modes of the system and new nonlinear modes that have no counterparts in the linear spectrum. Finally, we propose an optical realization of the linear and nonlinear Harper models in transversely modulated waveguide arrays, where these Hofstadter butterflies can be observed. This work is relevant to a variety of other branches of physics beyond optics, such as disorder-induced localization in ultracold bosonic gases, localization transition processes in disordered lattices, and more.

2 the energy carried by the wavepacket spread among the linear modes of the system [4,5,37]? What are the nonlinear modes of the system-both localized and extended, and what are their stability properties [3,6,7,38]. These questions become especially intriguing when the system has two (or more) incommensurable periodicities, a scenario appearing in many areas of physics. Perhaps the best known case of a linear wave system with two incommensurable periodicities is the one described by the Harper model, whose set of eigenvalues comprise the Hofstadter butterfly. The ubiquitous Hofstadter butterfly is known to describe a variety of incommensurate systems ranging from Bloch electrons in uniform magnetic fields [8]- [13] and cold atoms in optical lattices [14], to the Quantum Hall effect [15]- [17], and microwaves [18] and acoustic wave systems [19]. The purpose of this paper is to address the behavior of these generic systems under the action of nonlinearity.
The answers for many of the questions posed above are intimately connected to the underlying linear spectrum of the system. For example, the spreading of a localized input beam can be ballistic (if the system is linear and perfectly periodic [20]), diffusive (if the system contains disorder [21]), or suppressed altogether (i.e. localized [22]- [24]), depending on the linear spectrum of the structure. Another example is the existence of nonlinear modes. A localized nonlinear mode, a soliton of propagation constant β, can exist in a system only when the system has an exponentially decaying (and usually, unbounded) linear mode that has the same propagation constant β (eigenvalue) as the soliton. As such, the propagation constant of a soliton typically resides in the gaps of the linear spectrum (either in the semi-infinite gap beyond the first band, or in a gap between two bands). Specifically, closely examining the extended nonlinear modes of a periodic system reveals that the nonlinear spectrum (the collection of the eigenvalues of all the possible nonlinear extended modes) may consist of two parts: modes that are 'continuations' of the linear modes of the system and new modes that have no counterparts in the linear spectrum. For example, the linear spectrum of a continuous linear periodic system consists of bands separated by gaps. However, the spectrum of the corresponding nonlinear problem may involve, in addition to deformed bands, also a swallow-tail structure [25]- [27].
In the context of one-dimensional lattices, an interesting case arises when the potential is described by two incommensurable periodicities. In the linear domain and under the tightbinding approximation, this fundamental problem is described by the Harper equation [8]: ψ n+1 + 2 cos (2π αn + θ ) ψ n + ψ n−1 = Eψ n . (1) As first recognized by Azbel [9], this same equation can be used to describe the spectra of electrons in a two-dimensional square periodic potential under the action of a uniform perpendicular magnetic field. In this context, ψ n is the amplitude of the local mode at lattice site n, α is a dimensionless parameter-the magnetic flux through a unit cell in units of magnetic flux quanta, θ represents a transverse wave number and the eigenvalue E is the energy. The set of eigenvalues of this equation (for all values of the parameters α and θ ) was subsequently shown by Douglas Hofstadter [10] to have a fractal structure that bears his name-the Hofstadter butterfly. It was later proved rigorously that the spectrum is a set of zero Lebesgue measures [11] for (almost all) irrational α. Evidently, adding a nonlinear term to the Harper equation could change the physics of this system dramatically. So far, this problem has remained largely unexplored: only the dynamics of localized wavepackets in a nonlinear version of this equation and in a related quasiperiodic polaron model have been considered [28]- [30]. Contradictory conclusions of these papers were recently settled in [2], where the anomalous diffusion of the Harper model at criticality was shown to be destructed due to weak nonlinearity. Naturally, the nonlinear Harper model raises several basic questions. How does the nonlinearity affect the linear spectrum of nonlinear eigenmodes-the Hofstadter butterfly? Is the butterfly still fractal, under nonlinear conditions, or not? How does the butterfly vary as a function of nonlinearity, and how would the sign of the nonlinearity (self-focusing/attractive versus self-defocusing/repulsive) affect the butterfly? Finally, in what physical settings can these phenomena be observed?
In this paper, we show that the sets of eigenvalues of the nonlinear Harper model give rise to deformed versions of the Hofstadter butterfly, for either sign of the nonlinearity (self-focusing and self-defocusing). We find that the extended modes of the nonlinear Harper equation are classified into two families: modes that are 'continuations' of the linear modes of the system and new modes that have no origin in the linear spectrum. Finally, we propose to observe the linear and nonlinear Hofstadter butterfly spectra in transversely modulated optical waveguide arrays. Such waveguide lattices can provide a realistic platform upon which one could study for the first time not only the static properties but also the dynamic behavior of nonlinear Harper lattices. We emphasize that this work can shed light on other areas of physics dealing with wave behavior in nonlinear random lattices. These include, for example, localization induced by disorder in ultracold bosonic gases experiencing the action of purely random potentials and bichromatic optical lattices [31], and the onset of localization dynamics in Aubry-André systems [32,39].
Consider the nonlinear eigenvalue equation associated with the nonlinear Harper model [2,30]: In the context of optics, equation (2) describes the propagation of an optical monochromatic field in a Kerr nonlinear array of coupled waveguides, provided that the lattice satisfies Harper's rule (involving two incommensurable periodicities). In order to study the spectrum of the extended stationary states of this equation, we compute the extended nonlinear eigenmodes and eigenvalues for a large set of rational and irrational values of α, for many values of θ and for different maximal intensities: I max = max(|ψ n | 2 ). Clearly, since the equation is nonlinear, its eigenvalues and eigenmodes depend on the maximal intensity. The computation is carried out using two methods. For a rational α= p/q, where p and q are relatively prime and taking the boundary condition to be ψ n+q = exp(ikq)ψ n , equation (2) reduces to a nonlinear q × q eigenvalue matrix equation, which we solve using an iterative self-consistency method. We also employ a second method, appropriate for both rational and irrational α's-we guess values for E, ψ 0 and ψ 1 and then use equation (2) to compute the series ψ 2 ,ψ 3 , . . . , ψ N up to a large N . If the series does not diverge, then a nonlinear mode and its eigenvalue E are obtained. The identified nonlinear modes can be classified into two categories: (i) modes resulting as a 'continuation' of linear modes of the system and hence admit similar symmetries related to the periodicity of the underlying lattice, and (ii) modes that are unique to the nonlinear system and do not admit any periodicity. While the linear modes are quasiperiodic with several peaks in the Fourier spectrum (see figure 1(a) showing a linear mode with E = −2.902 11 at the upper edge of the lowest band for α = 1/5, and its spectrum in figure 1(b)), the nonlinear modes show much more diverse behavior. For example, nonlinear modes that happen to be 'extensions' of linear Harper modes possess similar quasiperiodic properties (figures 1(c) and (d)). On the other hand, modes with a much broader spectrum that is still centered around a few peaks (figures 1(e) and (f)) and chaotic modes whose spectrum is almost uniform (figures 1(g) and (h)) are also possible. This type of chaotic modes is unique to the nonlinear Harper system, as it appears that it cannot be obtained as a natural 'continuation' from a linear mode. To check the stability of the extended modes, we have used two methods: analytic linearization and the numerical beam propagation method where the modes are propagated with added noise. All the modes that we have checked were found to be stable. The eigenvalue of a nonlinear mode depends on the intensity of the mode; hence the spectrum of the nonlinear Harper model depends on the maximal intensity. Plotting the eigenvalues of the nonlinear system as a function of α, one gets the nonlinear counterparts of the Hofstadter butterfly. Figure 2 depicts the linear Hofstadter butterfly (a) and the nonlinear Hofstadter butterflies for three different maximal intensities (b)-(d), for a self-focusing nonlinearity. As the intensity is increased, the diagram goes up and deforms, losing its reflection symmetry about the E = 0 line while keeping its reflection symmetry about the α = 1/2 line.
For a self-defocusing nonlinearity, the spectrum is exactly the negative of the spectrum of the focusing nonlinear Harper model. That is, if ψ n is a nonlinear eigenmode with eigenvalue E of the focusing nonlinear Harper model for a specific α and a specific maximal intensity, then ϕ n = (−1) n ψ n is an eigenfunction of the defocusing nonlinear Harper model with the eigenvalue -E, the same α, the same maximal intensity, and with the transverse wavenumber changing as θ → θ + π . This does not imply that the nonlinear dynamics is the same. For example, a stable eigenmode in a focusing nonlinear system may become unstable in a defocusing system.
Next, we propose an incommensurate optical structure where both the linear and nonlinear Hofstadter butterflies can be observed. In this system, the linear Hofstadter fractal is formed from the propagation constants of the corresponding extended linear modes, whereas the propagation constants of its nonlinear extended modes produce the nonlinear Hofstadter butterfly. The proposed optical realization is a nonlinear optical waveguide array, where the widths of the individual channels and their separations are transversely modulated. By varying 6 these parameters, one can ensure that light propagation in this structure is described in the tight-binding approximation by the nonlinear Harper model. This structure could be fabricated in Fe-doped LiNbO 3 using Ti in-diffusion (similar to the structure used in [33] to demonstrate discrete modulation instability), where the nonlinearity originates from the large photovoltaic-photorefractive effect [34]. The propagation constants of the eigenmodes of such a structure can be measured using the prism method demonstrated in [35].
In order to choose the proper widths and separations for the waveguide channels, we first need an appropriate model for the modulated waveguide array. We then develop its tight-binding analogue using coupled mode theory. Starting with the normalized non-dimensional nonlinear Schrodinger equation: where A is the slowly varying amplitude of the electric component of the optical field, the potential n(x) is the refractive index profile of the array, and the nonlinear term represents the optical Kerr effect. We then describe the field A in terms of the eigenmodes of the individual waveguides, where ψ m is the normalized ground state of the mth waveguide, with eigenvalue β n , and g n ≡ ψ 4 n dx is a site-dependent effective nonlinear coefficient. By inserting this expression into equation (3) and by considering only nearest-neighbor coupling effects, we obtain the following equation: i ∂ E n ∂z + (β n + C n )E n + K n+1,n g n g n+1 E n+1 + K n−1,n g n g n−1 E n−1 + |E n | 2 E n = 0.
Here C n ≡ f n ψ 2 n dx is a modification of the propagation constant due to the presence of other waveguide channels, K n±1,n ≡ f n±1 ψ n±1 ψ n dx is the coupling between waveguide site n and its nearest neighbors n ± 1, and f m (x) = n(x) − n m (x) is the index profile of the waveguide array when n m (x)-the profile of the mth waveguide-is removed. Equation (2) contains an on-site potential term β n + C n , and for each waveguide there are two different coupling coefficients to its nearest neighbors. Comparing this equation, which models our structure, to the nonlinear Harper system, one can see that the nonlinear Harper equation is in fact a special case of our equation. If we judiciously choose the waveguides' parameters so that (β n + C n ) = 2 cos(2π αn), and K n±1,n √ g n /g n±1 = 1 for every n, then the propagation of light in our structure would be governed by the nonlinear Harper equation. We find that, for all practical purposes, this is possible, since we can find parameters for which the coupling values K n±1,n √ g n /g n±1 differ only slightly from unity, with the difference being equal to α sin(2παn). In our optimization of the structure's parameters, we restrict ourselves to varying the width and separation of each waveguide without changing the refractive index. Changing the refractive indices of the waveguides would give more degrees of freedom, so equation (4) would correspond more accurately to the nonlinear Harper model. Since light propagation in such a structure is governed by the nonlinear Harper model, measuring the linear propagation constants of the linear (or low-intensity) modes in this structure will directly provide an optical realization of the Hofstadter butterfly. On the other hand, measuring the nonlinear (or highintensity) propagation constants will lead to the observation of nonlinear Hofstadter butterflies.

7
Before closing, we would like to discuss an important feature of the nonlinear Hofstadter butterfly: its fractal dimension. Since all fractal systems are characterized by their non-integer dimension, it is compelling to ask how does the fractal dimension of the nonlinear butterfly vary when the nonlinearity is increased. This question, however, turns out to be difficult to answer conclusively. Basically, the linear Hofstadter butterfly is multi-fractal, that is, it has fractal dimensions that depend on the scale. Likewise, the nonlinear butterfly is also multi-fractal. In the nonlinear regime, these scale-dependent fractal dimensions vary in a non-uniform fashion as the nonlinearity is increased. Our figure 2 highlights this feature nicely: as the maximum intensity (nonlinearity) is increased, the spectrum becomes denser in some places and less dense in others. It is therefore difficult to characterize the fractal dimension of the nonlinear butterfly, as the nonlinearity is increased, because the fractals vary in a non-uniform fashion. Undoubtedly, this is a very intriguing issue for future research. In the same vein, it would be very interesting to study transport in weakly nonlinear systems characterized by the nonlinear Hofstadter model. More specifically, recent experiments have shown that a weak attractive nonlinearity enhances localization in periodic systems containing disorder [23]. Similar experiments with quasicrystals containing disorder reveal that disorder can enhance transport in linear quasi-crystals, while nonlinearity may also play an important role-especially when the evolution is short range [36]. Clearly, the results presented here open up many new intriguing questions.
In summary, we have demonstrated that the nonlinear Harper model supports stable extended nonlinear modes with varying structural properties. In all cases, the resulting energy diagrams lead to deformed versions of the Hofstadter butterfly. We have proposed an optical realization based on a modulated nonlinear waveguide array that can be used to experimentally study the properties of linear and nonlinear Harper models. This model offers a wide scope for theoretical and experimental research as to its nonlinear modes (solitons), their dynamics and much more. Especially intriguing is the question of how do the fractal dimensions of the (multi-fractal) nonlinear butterfly vary as the nonlinearity is increased. We envision a variety of nonlinear phenomena where the dynamics in a Hofstadter-type waveguide array (made of two incommensurable periodicities) would be fundamentally different from that in a periodic array, such as irreversible (possibly chaotic) propagation dynamics, unique shock waves and much more.