Bistability, multistability and nonreciprocal light propagation in Thue-Morse multilayered structures

The nonlinear properties of quasiperiodic photonic crystals based on the Thue-Morse sequence are investigated. The intrinsic spatial asymmetry of these one-dimensional structures for odd generation numbers results in bistability thresholds which are sensitive to the propagation direction. Along with resonances of perfect transmission, this feature allows to achieve strongly non-reciprocal propagation and to create an all-optical diode. The salient qualitative features of such optical diode action is readily explained through a simple coupled resonator model. The efficiency of a passive scheme, which does not necessitate of an additional short pump signal, is compared to an active scheme, where such a signal is required.


Introduction
The experimental discovery of quasicrystals [1] has given rise to a sudden breakthrough in the area of solid state physics. It had been assumed for a long time that a periodic arrangement of unit cells was a main prerequisite for the specific properties of crystalline matter. The structure of quasicrystals does exhibit long-range order or correlations between distant parts, but at the same time there is no underlying periodicity, in the sense that a shifted copy of the crystal never matches exactly the original one. In fact, quasicrystals represent an intermediate stage between random media and periodic crystals, effectively combining both localization properties as a result of short-range disorder and the presence of band gaps due to long-range correlations [2].
The concept of quasiperiodicity was readily transferred to photonic crystals and proved to be of great value for many practical purposes [3][4][5][6]. Photonic quasicrystals are deterministically generated dielectric structures with a non-periodic modulation of the refractive index. In the one-dimensional (1D) case, they can be formed by stacking together dielectric layers of several different types according to the substitutional sequence under investigation (Cantor, Fibonacci, Rudin-Shapiro, Thue-Morse, etc.) [7]. The Fibonacci sequence is of particular importance, since it leads to the existence of two incommensurable periods in the spatial spectrum of the structure. Such behavior is typical of sequences with a so-called pure point spectrum, which makes the Fibonacci sequence truly quasiperiodic, as a consequence of the appearance of Bragg-like peaks in the spatial spectrum [8]. This property has been demonstrated to be very valuable for nonlinear optical applications, such as third harmonic generation [9]. In fact, the latter is a two-stage process, and for each stage different phase matching conditions should be met. The Fibonacci sequence allows to fulfill both of them simultaneously and on the same crystal [10].
In contrast to the Fibonacci sequence, the Thue-Morse (ThM) sequence possesses a singular continuous spectrum, which is neither continuous nor singular [11][12]. Thus, strictly speaking, the ThM sequence is not quasiperiodic, but rather deterministically aperiodic. To an even different class belongs the Rudin-Shapiro sequence, which shows evidence of a continuous spectrum, analogous to the one exhibited by random sequences [13]. With a slight abuse of terminology, we shall always refer to the above sequences as quasiperiodic, and we shall focus our attention specifically to the ThM sequence, for reasons explained below.
The bandgaps of quasiperiodic multilayered structures, also called 'pseudo band gaps' [14] or 'fractal bandgaps' in the literature [15], often contain resonant states, which can be considered as a manifestation of numerous defects distributed along the structure. Their quality factors are not very large when compared with those of symmetrically placed defects inside conventional Bragg structures, but the mode profiles associated to these resonances are extended in space and, as a consequence, very suitable for the enhancement of nonlinear effects throughout the whole body of the crystal. In order to investigate the specific properties of bistability and multistability in quasicrystals, in the following we shall restrict our attention exclusively to 1D structures. As a representative case, we focus our attention to the ThM aperiodic multilayer, and in making this choice we were guided by the following two reasons. Firstly, the spatial asymmetry of the ThM sequence (of odd generation numbers) interplays with the nonlinearity and makes the transmission sensitive to the propagation direction [16]. Secondly, almost all resonances in ThM quasicrystals are resonances of complete (i.e. 100%) transmission, with some of them being situated well inside the pseudo band gap regions. Taken together, the above two properties provide favorable conditions for the design of an all-optical diode, i.e. a device that shows a strong contrast in the transmission between forward and backward incidence. A proper understanding of bistability phenomena in quasiperiodic crystals is very important for the understanding of the optical diode action described above.
The paper is organized as follows. In section 2, we briefly review the well-known linear properties of ThM quasicrystals and classify their resonances by using the method of trace maps. It is shown that the field profiles at resonance frequencies follow the pattern of Thue-Morse sequence and the classification suggested can be used as a measure of this self-similarity. In section 3, we show that the spatial asymmetry inherent to ThM quasicrystals of odd generation numbers can interplay with Kerrnonlinearity in such a way that transmission become sensitive to the propagation direction. It is very important that the dynamics of this interplay can have a sudden jump due to bistability when the intensity of the incident light changes, and the thresholds for such jumps are different for the incidence from the left and from the right of the multilayer. In section 4, we apply a phenomenological model based on a coupled resonator model to explain this behavior analytically for both bistable and multistable cases. The general formulas describing nonlinear transmission spectra are derived and their relationship to the level of self-similarity in the field profiles is emphasized. In section 5, we propose a design of nonlinear optical diode based on the differences in bistability thresholds between forward and backward propagation. The efficiency of two schemes is compared: passive, when only one pulse is used, and active, when an additional short pump signal is applied to facilitate switching between stable branches of hysteresis. The last section summarizes the results and presents the conclusions. The appendix provides some details about the FDTD method used in the paper for the dynamical simulations of multilayered structures with instantaneous Kerr-nonlinearity. above methods, and the presence of two incommensurable periods in their spatial spectra is implicitly related to the projection of a two-dimensional grid onto a one-dimensional line [18]. However, these two approaches are not equivalent to each other, and the use of substitutional sequences tends to be a more general procedure [19]. In particular, this is applicable to ThM sequences and makes it distinct from strictly quasiperiodic sequences. Sometimes, it is called a deterministic aperiodic sequence, in order to emphasize that it has more disorder than quasiperiodic ones and stands closer to random sequences [20]. On the other hand, the development of resonances when changing generation number in ThM quasicrystals resembles the analogous development in periodic structures to a certain extent and brings these two kinds of crystal closer to each other [21].
Two dielectric layers of different materials are required to compose ThM quasicrystals, which shall be indicated with the letters 'A' and 'B' . They should be arranged in the same way as in literal ThM sequence, which is governed by the following substitutional ('inflation') rules: 'A' 'AB' → , 'B' 'BA' → . Thus, starting from a single layer 0 'A' S = , which is defined to be the ThM quasicrystal of 0th generation or ThM 0 , one obtains 1 'AB' S = , 2 'ABBA' S = , 3 'ABBABAAB' S = and so forth, with each step giving a sequence of generation number increased by one (figure 1a). It can be shown that an additional recurrence relation follows from this definition, which holds for ThM sequences as single blocks: , where n is the generation number. In this notation n S indicates a sequence 'conjugated' to n S , where all letters are interchanged to their opposite as in the rule 'A' 'B' ↔ .
It is interesting to note that there is only a slight difference between the ThM and Bragg (or periodic) sequences as far as inflation rules are concerned. More specifically, Bragg sequences possess the inflation rule 'A' 'AB' → , 'B' 'AB' → . Therefore, the growth rate of these sequences is the same, and the total number of layers is subject to a rapid (exponential) increase as 2 n . The occurrence of layers 'A' and 'B' is also the same for both ThM and Bragg structures for any generation number. Nevertheless, such a small modification of the inflation rules leads to completely distinct transmission and localization properties.
The primary method to compute transmission spectra and field profiles for multilayered structures is the transfer matrix method [22]. If normal incidence is assumed, the electric E and magnetic H fields at the opposite sides of each single layer can be related in the following way: where m d is the thickness of the layer with refractive index m n , and / c Ω = ω is the wavenumber in vacuum. In the representation of equation (1) the imaginary unit is deliberately kept with the magnetic field, in order to show explicitly that transfer matrices of nonabsorbing layered structures can be characterized by real matrix elements. In general, the refractive indices of the materials can depend on frequency, but for the sake of simplicity we will assume that all materials are dispersionless. By multiplying together all transfer matrices of subsequent layers, it is possible to construct the total transfer matrix of the structure and to find the transmission spectrum in terms of its matrix elements Another -more efficient -approach to obtain transfer matrices of ThM quasicrystals is to make use of the following two-step recurrence relations layers 'A ' and 'B' , respectively. The subscripts in curly brackets indicate generation numbers of corresponding ThM quasicrystals. Formulas (5)-(6) take their origin in the recurrence relations of ThM sequence explained above, and stress explicitly the interplay between 'ordinary' n S and 'conjugated' n S counterparts of the sequence. The self-similarity of ThM quasicrystals is naturally embedded in this approach, and this simplifies considerably the computation of transmission spectra for large generation numbers. However, even this modified representation of the transfer matrix method still contains redundant information about the structure, if one is interested only in the analysis of the properties of the resonance frequencies. In this case, it is sufficient to operate with traces of transfer matrices. The trace is defined as 11 22 Tr( ) For ThM quasicrystals of fixed generation numbers, these traces are unchanged under conjugation and for adjacent generation numbers, an additional relation or 'trace map' holds [23]: -5 - The formulas (8)-(9) are valid for 1 n ≥ , and in derivation of (9) the Cayley-Hamilton theorem was applied to prove that any unimodular matrix satisfies the relation where U x is the trace of a unimodular matrix U , and I is the identity matrix. To use the equation for the trace map (7), two initial conditions are required. They can be found directly from the multiplication of transfer matrices having the form given by equation (2): 2 2cos 2 cos 2 ( / / )sin 2 sin 2 where . Several important conclusions can be extracted from the trace map of ThM quasicrystals. The analysis simplifies greatly in absence of absorption, so that the imaginary parts of the refractive indices are negligible. In this case, all traces will be real, and there will be frequencies in the spectrum in correspondence to which 0 n x = . Independent of the value of 1 n x + , this makes 2 2 n x + = , which coincides with the trace of identity matrix. This simply means that frequencies for which 0 n x = correspond to frequencies of perfect transmission. To prove this, an auxiliary matrix should be constructed, the trace of which gives the difference of diagonal elements of the original transfer matrix: 11 22 Tr( ) Therefore, the resonances of complete transmission in ThM quasicrystals can be characterized by identity transfer matrices, and they are preserved from one generation to another. Only additional resonances can appear in new generations, and the overall transmission spectrum of these quasicrystals shows a fractal nature. As generation number increases, the spectrum reveals self-similarity and characteristic trifurcation of resonances [23].
In what follows, we will consider that layers 'A ' and 'B' have the same optical thickness, , and we will introduce the reference wavenumber 0 Ω , defined as 0 / 2 A A n d Ω =π , which corresponds to the quarter wavelength condition. The ratio 0 / Ω Ω will be called 'normalized frequency' and will make the comparison of different structures easier. For example, a normalized frequency 0 / 1 Ω Ω = corresponds to the middle of bandgap region for Bragg crystals, but it is the centre of a wide 'pseudo-pass band' for ThM quasicrystals, a region made of densely located shallow resonances, not always of perfect transmission (figure 2). This band exists only for the specific case of equal optical thicknesses of the layers, and it starts to appear from the second generation, for which formulas (8)-(9) are not valid. The above condition makes them applicable for 0 n ≥ . Moreover, as long as the condition is fulfilled, the transmission spectra of all 1D structures will be translationally invariant along frequency axis,

Bistability and multistability in Thue-Morse multilayers
The bistability and multistability properties of Bragg gratings have been studied in depth by a number of authors [24][25][26]. The knowledge of the field profiles is crucial for the proper understanding of bistability and multistability in ThM quasicrystals. In correspondence to the resonance frequencies these field profiles show a self-similar pattern which strongly resembles that of a ThM sequence of smaller generation number (figure 1b). Due to this fact, these resonant states are also called lattice-like states [27]. The level of self-similarity depends on the generation number in correspondence to which a particular resonance first appears. For example, the field profile shown in figure 1b belongs to a ThM quasicrystal of 7th generation. It corresponds to the resonance with normalized frequency 0 / 0.756 Ω Ω = , and, according to the trifurcation diagram shown in figure 2, this resonance is inherited from the 4th generation. Therefore, the field profile mentioned above has the level of similarity equal to 7 4 3 − = , and it is possible to distinguish 3 2 8 = independent blocks. Moreover, these blocks are arranged according to the ThM sequence of 3rd generation. We have found that for any ThM structure with layers of equal optical thickness, the level of self-similarity is maximal for normalized frequencies belonging to the following series: -7 -However, the localization strength is relatively small at these frequencies, and thus the nonlinear response cannot be enhanced significantly, because the electric field profile is almost flat along the structure. For this purpose, resonances located near the edges of pseudo band gaps are more suitable.
In general, the localization strength is inversely proportional to the level of self-similarity, and when the length of the structure increases, the resonant states gradually change from localized to extended ones [28]. The case of Kerr (or cubic) nonlinearity was considered in this work, so that a nonlinear refractive index was additionally taken into account for each type of layer. This leads to an intensity-dependent self-phase modulation, which is able to shift resonance frequencies [25]. The direction of this shift is determined by the sign of the Kerr coefficient. The most interesting cases evidently occur when the sign of the nonlinearity is such that resonances in the spectrum of transmission bend towards the bandgap regions (figures 3-4). In the following, we always choose positive nonlinear coefficients for the two materials, One of the most serious advantages provided by ThM sequence is that for odd generation numbers the corresponding photonic structures are intrinsically asymmetric, and nonlinearity is capable of making transmission sensitive to the propagation direction [16]. This feature is completely absent in nonlinear Bragg structures, where hysteresis curves are the same regardless of the direction of incidence. Although a similar nonreciprocal behavior can be achieved in the framework of linear optics, the latter typically requires making use of magneto-optical media with externally applied static magnetic fields [31] or chiral media such as cholesteric liquid crystals [32]. Moreover, the polarization  state is very important for the correct operation of these linear devices, while the design suggested here is free from such complications.
The level of self-similarity in the field profiles is also related to the type of hysteresis that can be observed near a particular resonance. In presence of several independent localization centers inside the structure, the interplay between them gives rise to multistability. It is related mainly to the fact that ThM structures of higher generation numbers can be decomposed into those of lower ones. Therefore, to single out a particular hysteresis type in its pure form, it is necessary to check that the corresponding resonance does not occur in previous generations.
Two characteristic examples are shown in figures 3-4. The first case (figure 3) corresponds to a resonance with no self-similarity in the field profile and clearly demonstrates bistability of transmission, whereas in the second case (figure 4) two independent localization centers can be distinguished. The hystereses are shown in figure 3b and 4b, the nonlinear transmission spectra are given in figure 3c and 4c. A detailed explanation of this behavior will be given in the next section. Notice that similarly to Bragg structures, the maxima of electric field can be located mostly in the layers of one type, so that in principle it would be sufficient to use only one nonlinear material.

Non-separable resonances
In order to provide a qualitative analytical understanding for the nonlinear behavior of ThM quasicrystals, we apply the coupled resonator model, which has been recognized as a useful tool for describing the dynamical properties of single-mode cavities in 2D photonic crystal waveguides [33][34][35]. Being a phenomenological theory based on general physical concepts like conservation of energy and time-reversal symmetry [36], it is not restricted to the above 2D systems and with little modifications we adapt it to 1D structures. Taking into account the intrinsic spatial asymmetry of resonances, the main equations of this method can be written in the following form In essence, they represent an extension of the scattering matrix method to the time domain. The first equation (17) describes the temporal evolution of a resonant mode in the vicinity of a resonant frequency 0 ω . The amplitude of this mode A is normalized in such a way that 2 | | A gives the energy inside the cavity. Energy is carried to the cavity from the two ports or channels, which are denoted as U and V , located on the left and the right sides, respectively (figure 5a). The amplitudes of the incident signals are denoted as u + and v + . They are normalized in such a way that 2 | | u + and 2 | | v + give the power flow. The interaction between these signals and the cavity is described by the two coupling coefficients u κ and v κ , respectively. Since the cavity not only accumulates the incident energy, but may also partly absorb it or return it back to the channels, there is a damping constant γ , which characterizes the efficiency of these processes. It can be related to the quality factor of the resonator 0 /(2 ) Q ω γ = .
As to the influence of nonlinearity, it results in the shift of the resonant frequency. In equation (17), we introduce a characteristic energy 0 P that allows to take this shift into account [37], which can be either positive or negative depending on the common sign of the Kerr coefficients in the multilayered structure.
The second equation (18) describes the output from the cavity. The right hand side of this equation consists of two parts. The first part characterises the background reflection from the cavity as if the amplitude of the mode inside it were zero. In fact, it has exactly the form of a scattering matrix. The background transmission coefficient t is equal for both directions even for nonsymmetric structures, while the reflection coefficients in general have different phases and are denoted as u r and v r . As to the second part of equation (18), it is proportional to the amplitude of the mode, and thus provides the corrections caused by the interaction of incident signals with the cavity. Assuming that the system operates at a fixed frequency, it is possible to derive from equations (17), (18) an explicit expression for the nonlinear transmission spectrum along the forward (backward) propagation direction: ) is the output intensity, 0 ( )/ = − δ ω ω γ is the normalized detuning from the resonance frequency, and The typical form of nonlinear transmission spectrum following from formula (19).
In the case of single Dirac-delta-like nonlinear layer in the system, formula (19) coincides with a previously investigated exact analytical solution [38]. However, taking into account a continuously distributed nonlinearity gives rise to a new feature: if the structure is spatially asymmetric, the characteristic intensity can depend on the direction of incidence, as it is evident from equation (20).
Since the formula (19) is equivalent to a third-order polynomial in the output intensity , out u v I , by definition it fails to describe resonances with multistability, but it gives a good approximation for resonances with a strictly bistable response (compare for instance figure 5b with figure 3b-c).

Separable resonances
In the more complex case, when the field profiles consist of several independent parts, one simply needs to increase the number of equations in the coupled resonator model. There should be exactly one set of equations (17)-(18) for each separable resonance. For example, the case shown in figure 4 can be described by the following system of equations The system (21)-(24) can be readily solved in the frequency domain. We found the following formula, which describes the nonlinear transmission spectra with multistable behaviour: Expression (25) contains two parameters, a ξ and b ξ , which can be calculated given the output signal v − (or u − ) by using the following ladder: The typical nonlinear transmission spectra calculated by using formula (25) are shown in figure 6bc. The main parameter which determines the shape of nonlinear transmission spectra is the phase mismatch a b + ϕ ϕ . It describes the efficiency of interaction between separable resonances and is responsible for the splitting of initially identical resonant frequencies even in linear case. The only difference between figures 6b and 6c is in the phase mismatch, and the two well-separated resonances shown in the first case become undistinguishable for small intensities in the second case (see gray lines). As the intensity of incident signal increases, this degeneracy disappears giving rise to multistability (see red and blue lines for the forward and backward incidence, respectively). Qualitatively, the nonlinear transmission spectra shown in figure 6c corresponds to that in figure 4c.
It is worth noting that if 0 , the formula (25) reduces to which means that the two resonances do not interact with each other and the resulting nonlinear transmission spectra reveals only a bistable behaviour.

Nonlinear optical diode action in Thue-Morse multilayers
The difference in bistability thresholds can be small, but it creates favorable conditions for unidirectional propagation, so that ThM quasicrystals can be used as all-optical diodes [16]. Similar to electronic circuits, these devices are indispensable, when there is a need to suppress the flow of light in one direction or to avoid problems caused by unwanted reflections. Various types of optical diodes have been proposed and realized, which can work both in linear [31][32] and nonlinear regime [39][40][41][42]. The primary figure of merit, which determines the efficiency of this device, is the contrast ratio f b / C T T = between transmission along the forward ( f T ) and the backward ( b T ) directions. C can be very large in ThM structures due to the fact that this device can operate on two different hysteresis branches depending on propagation direction (figure 3b).
It is possible to use the intensities of both up and down transitions to achieve strongly nonreciprocal transmission. In the former case, the scheme works in passive mode, but the maximum value of transmission is limited [43]. In the latter case, the transmission can be almost perfect, but the scheme requires an additional short pump signal in order to switch to the higher stable branch of hysteresis [44].
The FDTD method was applied to demonstrate the switching dynamics in time domain (figure 7). To maintain the second order accuracy of the Yee-scheme in case of multilayered structures, nonuniform spatial mesh was used with electric field nodes aligned to the boundaries between layers (see Appendix). Since the structure can accumulate and release energy, the sum of reflected (given by red dotted line in figure 7) and transmitted intensities (blue dashed line) should not necessarily give the incident one (gray solid line). This constraint becomes valid only after all the transient dynamics is over. The input intensity was set first to the value, which is sufficiently large for the up switching in the forward direction, but at the same time is smaller than corresponding value for the backward direction. The maximum transmission obtained was f 65% T = with the contrast ratio 7.2 C = , which is quite close to the theoretical limit of max 9 C = following from (31) in the assumption of passive scheme [43]. In the second stage, the input intensity was set just above the down switching threshold for the forward direction, which ensures transmission of almost f 100% T = in this direction, while in the opposite direction transmission will be strongly suppressed due to the presence of the pseudo band gap. The contrast ratio achieved was 20 C = , but it is not limited in principle. The major drawback of this scheme is that it is necessary to force the system to switch to the upper branch of hysteresis whether by applying an auxiliary pump signal, or by temporarily increasing the input intensity.

Conclusions and future work
The nonlinear properties of quasicrystals based on ThM sequence were investigated. It was shown that the interplay between the intrinsic spatial asymmetry of these structures for odd generation numbers and Kerr nonlinearity makes the switching thresholds induced by bistability and multistability sensitive to the propagation direction. The role of self-similarity was emphasized to explain the shape of hysteresis curves observed near a particular resonance, and conditions necessary to achieve highly non-reciprocal propagation were formulated. Coupled resonator model was used as a phenomenological model. FDTD simulations have been used to confirm the results obtained by the nonlinear transfer matrix method, and to show how ThM multilayers can exhibit an effective optical diode action, both in passive and active operation. To separate incident and reflected waves from the structure, the so-called Total-Field Scattered-Field technique is used on the left (and of the right) of the structure Smooth envelops for switching on and off monochromatic signals are defined by splines of 6th order, which ensure the continuity of fields up to the second derivative [47] ( / ) for 0 , ( ) 1 for , where T is the total time when the signal is nonzero, τ is the duration of the switching, and