Asymmetric Wave Propagation Through Saturable Nonlinear Oligomers

In the present paper we consider nonlinear dimers and trimers (more generally, oligomers) embedded within a linear Schr{\"o}dinger lattice where the nonlinear sites are of saturable type. We examine the stationary states of such chains in the form of plane waves, and analytically compute their reflection and transmission coefficients through the nonlinear oligomer, as well as the corresponding rectification factors which clearly illustrate the asymmetry between left and right propagation in such systems. We examine not only the existence but also the dynamical stability of the plane wave states. Lastly, we generalize our numerical considerations to the more physically relevant case of Gaussian initial wavepackets and confirm that the asymmetry in the transmission properties also persists in the case of such wavepackets.


I. INTRODUCTION
In the last two decades, the subject of nonlinear dynamical lattices has gained considerable attraction and interest due to its emergence and relevance to a wide range of diverse applications. These include, among others, arrays of nonlinear-optical waveguides [1], Bose-Einstein condensates (BECs) in periodic potentials [2], micromechanical cantilever arrays [3], Josephson-junction ladders [4], granular crystals of beads interacting through Hertzian contacts [5], layered antiferromagnetic crystals [6], halide-bridged transition metal complexes [7], and dynamical models of the DNA double strand [8].
On the other hand, a specific phenomenon that has been intensely explored in a wide variety of recent studies is that of potentially asymmetric (i.e., non-reciprocal) wave propagation. This has been examined e.g., in asymmetric phonon transmission through a nonlinear interface layer between two very dissimilar crystals [9]. Another example is the proposed thermal diode [10] which, in turn, led to its experimental realization [11]. The optical diode was theoretically suggested in [12] (see also [13] for a setup of unidirectional transmission in photonics crystals) and experimentally achieved in [14]. Such rectification effects have also been proposed in left-handed metamaterials [15], granular crystals [16] and systems with gain-loss bearing so-called PT -symmetry [17,18], among others.
In the present work, we focus on an, arguably simpler, implementation of the diode effect, which is fundamentally due to nonlinearity, as has been presented recently in the work of [19]. There, a linear chain was considered with a pair (or more) of nonlinear sites between the two linear ends of the chain. The nonlinear nature of the dynamics, coupled to a potential asymmetry between the characteristics of the nonlinear sites, was at the heart of the asymmetric propagation observed. The nonlinear sites were modeled as a prototypical system that has arisen in numerous applications, either as a direct model of relevance or as an envelope approximation in the form of the so-called discrete nonlinear Schrödinger (DNLS) equation [20]. While structurally simple, this model incorporates the fundamental characteristics of such lattices, namely diffraction (i.e., a discrete analogue of dispersion) and nonlinearity. The scenario of [19] explores the standard cubic nonlinearity associated with the Kerr effect [1]. However, often in applications, other types of nonlinearities are important as well. For instance, defocusing lithium niobate waveguide arrays exhibit a different type of nonlinearity, namely a saturable, defocusing one due to the photovoltaic effect [21]. In the latter context, dark solitons have been identified not only in regular homogeneous lattices, but also in higher gaps [22], as parts of multi-component soliton complexes (such as dark-bright solitary waves) [23], and even in heterogeneous chains with alternating couplings [24].
In the present work, we combine the experimentally accessible form of the saturable nonlinearity with the asymmetric propagation nonlinear phenomenology of [19]. The model setup is presented in Section II. We consider in Section III exact plane wave solutions by solving the linear parts of the chain and gluing them through the nonlinear saturable "defect" sites. In this way, we identify the asymmetry between left and right transmittivities, due to the nonlinear propagation through the asymmetric dimer, trimer or more generally oligomer (i.e., few site) configuration. We explore the stability of these configurations and generically identify them as unstable. The dynamical evolution of the corresponding instabilities is explored through direct numerical simulations. Finally, more realistic (for experimental purposes) wavepackets of a Gaussian form are also considered in Section IV and the manifestation of the asymmetry in propagation under such initial data is systematically quantified. Section V summarizes our findings and presents our conclusions including some possible directions of future work.

II. THE MODEL
Motivated by the above application of lithium niobate waveguide arrays, we consider a nonlinear Schrödinger type chain with governing equation for γ n ∈ R and φ n ∈ C. Here t plays the role of the (spatial) evolution variable. The saturable nonlinearity term is present in a finite region in the middle of the chain. That is, γ n = 0 only for 1 ≤ n ≤ N . The wave propagates freely i.e., linearly outside of the finite region containing the nonlinearity. The system is Hamiltonian [25] with We will examine the transmission properties of stationary solutions that take the form of plane waves on the linear portions of the lattice. We also explore the linear stability analysis for these solutions and for a number of dynamically unstable scenarios, we evolve the corresponding solutions through direct numerical simulations over the propagation parameter t. Finally, we examine the propagation of more physically applicable localized Gaussian wavepackets and summarize the corresponding transmission properties in connection to the corresponding (potentially observable experimentally) asymmetries.

A. Plane Waves
We seek standing wave solutions by setting φ n = ψ n e −iωt for ω ∈ R. This gives a set of algebraic equations which can be written in the form of a backwards transfer map for ψ n ∈ C independent of t. Following the procedure in [18] and for now assuming k ≥ 0, we begin by assuming solutions in the form of plane waves on the linear portions of the chain. That is, 3 with R 0 , R, T ∈ C representing the incident, reflected and transmitted amplitudes, respectively. The solution (4) solves Eq. (3) for n ∈ {1, . . . , N } only if the wavenumber k satisfies ω = −2 cos(k). Also directly from (4) with n = 0, 1 we have The stationary solution across the whole lattice is then known by the following procedure: given values for γ n for n ∈ {1, . . . , N } we start by specifying values for k and T . Then we compute ψ 0 , ψ 1 , . . . ψ N −1 via Eq. (3) and then R, R 0 by (5). Such a procedure of finding the input as a function of the output is referred to as a "fixed output problem" [26]. Stationary solutions where the amplitude R 0 is incident from the right-hand-side and the wavenumber is taken as −k ≤ 0 can also be formulated in a similar way. In order to avoid swapping the format of (4) (so that R 0 , R would apply on the right and T on left), it is more convenient to leave (4) as-is with positive wave number k and instead flip left-to-right the configuration of the nonlinearities, i.e. {γ 1 , γ 2 , . . . , γ N } → {γ N , γ N −1 , . . . , γ 1 }. In this way the computation of the solution for negative wavenumber is unchanged from the above outline aside from the swap of the order of the γ's. Plots of these plane wave stationary solutions are shown in the next section, where we also address the stability of their t-propagation. In practice we truncate the lattice and refer to its finite length as L.
For a solution that has been determined by the processes described above, we next compute the transmission coefficient τ def.
= |T | 2 /|R 0 | 2 explicitly assuming that T is given. For this purpose it is convenient to write ψ n def.
= T e ikN Ψ n with n = N − l so that l = 0 corresponds to n = N , and incrementing l corresponds to decreasing n. In this notation we have Ψ N = 1 for l = 0 and the value of Ψ n for each subsequent node towards the left is given by rewriting Eq. (3) as See the Appendix where we record a few iterations of Eq. (6). Then by (5) with ψ 0 , ψ 1 computed according to Eq. (6) we have and In the linear case (γ 1 = γ 2 = 0) and in the symmetric case ({γ 1 , γ 2 , . . . , γ N } = {γ N , γ N −1 , . . . , γ 1 } as an ordered set) it is immediately seen that τ is the same for waves incoming from the left or right side. For N = 1 the transmission is always symmetric. We also define here a quantity to measure the asymmetric propagation. We will use the definition of a rectification factor f in the form of where the quantity τ (k, T ) corresponds to transmission of a left-incoming wave with positive wavenumber k ≥ 0 and τ f lip (k, T ) with k ≥ 0 is equivalent to the transmission of a right-incoming wave with negative wavenumber (recall the process described above of keeping k positive while flipping the order of the γ's). This way nonzero values of f in the range [−1, 1] measure the asymmetry of transmission in the system. Symmetry in transmission corresponds to f = 0 and f > 0 corresponds to greater transmission of incident waves originating from the left (transmitted on the right) as compared with incident waves originating from the right (transmitted on the left). Of course, f < 0 corresponds to greater transmission of waves originating from the right. Figure 1 shows plots of the transmission coefficient τ and the rectification factor f as a function of the amplitude T and the wavenumber k of the extended plane wave solutions. We find that whether more is transmitted for waves incoming from the right or left is variable as a function of k and T . Notice that values for γ's are chosen in Figure 1 to be such that γ 1 < γ 2 in the N = 2 case and γ 1 < γ 2 < γ 3 in the N = 3 case. In other words, with increasing γ's from left to right we observe that transmission properties vary with the choice of the parameters T, k. We find that in accordance with our above analysis the N = 1 case is symmetric. Although we do not show an N = 1 analogue  Figure 2 shows the transition in the eigenvalue plots as we move towards brighter regions of the | min(Im(ν)| diagram.
of Figure 1, such plots look similar to Fig. 1 but there is exact symmetry and f = 0 for all T, k. It is interesting to point out here that the rectification factor appears to acquire its largest (absolute) values for k close to π i.e., at the edge of the Brillouin zone. Furthermore, both in the N = 2 and in the N = 3 case, the dependence of f near this value appears to be a non-sign-definite function of T (i.e., different ranges of T values appear to favor propagation in one or the other direction).

B. Stability
In order to analyze spectral stability of stationary states of the form discussed in the previous subsection we write φ n (t) = e −iωt ψ n + ε a n e iνt + b n e −iν * t (10) for a n , b n , ν ∈ C, ε small, and with ψ n being a stationary solution from the previous section. The resultng linear stability equations then read for where G is a sparse matrix with ones on both the super-and sub-diagonals. Given a stationary plane wave solution ψ n and values of γ n which encode the nonlinearity for 1 ≤ n ≤ N , one then calculates the eigenvalues ν in (11). If ν has a negative imaginary part this indicates that the perturbed solution φ n (t) is unstable, as is easily seen by Eq. (10). In practice, one diagonalizes a finite truncation of the matrix in (11), ensuring that the relevant eigenvalues are not affected by the truncation error. In other words, F 1 , F 2 , F 3 , F 4 and G are all L × L matrices and in the matrix Eq. (11) it is now convenient to think of a n and b n as length L column vectors. Furthermore, the Hamiltonian symmetry of the solution ensures that the relevant instability eigenvalues come either in pairs (if ν is imaginary) or in quartets (if ν is genuinely complex).
In Figures 2, 4 we show a plot of min(Im(ν)) as a function of T and γ. We find that an increase in the magnitude of a γ i parameter (with other nearby γ's held fixed) leads to min(Im(ν)) of larger magnitude indicating greater instability. Figures 2, 3, 5, 6 show eigenvector and eigenvalue plots alongside snapshots of φ n (t) to show the behaviour of typical propagation in the t variable of the unstable plane waves. The boundary conditions are calculated according to (4) at t = 0 and evolved by multiplying by e −iωt for t > 0 so as to conform with (10). The unstable plane wave solution, when propagated in the evolution variable, exhibits a few effects: if k > 0 then amplitude leaks over to the right-hand side (to the left if k < 0), and due to the localized instability a peak appears in the center of the lattice. Of course, given the conservation laws of the system, the power n |φ n (t)| 2 and the Hamiltonian H(t) are preserved over t. The figures also show a transition in the eigenvalue plots for unstable solutions. A weak instability (corresponding to dim but nonzero regions of the min(Im(ν)) plots of Figures 2, 4) results in eigenvalue plots in the complex plane where a quartet appears off the real axis; see Figures 2 and 5. As the instability is enhanced for larger values of γ (comparably brighter regions of the min(Im(ν)) plots), the two pairs constituting the quartet merges on the imaginary axis and subsequently split with one pair headed towards zero; see Figures 3,5. For the highest magnitude of instability (brightest regions on the plots of min(Im(ν))) the eigenvalues indicating the instability are in the form of a pair on the imaginary axis; see Figures 3,6. In the examples shown, the instability generically appears to transport power to the right part of the lattice, deforming (decreasing the power of) the corresponding n < 0 portion of the plane wave. On the other hand, critically (per the localized eigenvector of the instability), a localized mode appears to form at the central nonlinear nodes within the domain.
In comparing our results with that of [19], we find that the asymmetry associated with the saturable nonlinearity (presented here) is less pronounced to that of a system with a cubic nonlinearity and a linear potential term (presented in [19]). In the t propagation of extended solutions we can also compare the top plot in Fig. 3 of [19] with our Fig. 7 in which we show plots over space and the propagation parameter. The two systems both experience a concentration of amplitude at the center of the lattice as t moves forward. In the case of the cubic nonlinearity in [19] there are three concentrations of amplitude (two of which are moving). Here we see only the one central concentration of amplitude at the center while sites nearby the center drop amplitude in comparison to the highest points of the initialized state at t = 0. Also, in contrast to [19] where the amplitude concentrations more dramatically rise above the background, here the central concentration of amplitude is more similar to the maximal amplitude of the initialized state.

IV. PROPAGATION OF A GAUSSIAN
Finally, in this section we look at the propagation of a Gaussian wavepacket through the lattice for each of the cases N = 1, 2, 3. While it is less straightforward to prepare the delocalized initial conditions needed for the plane wave solutions of the previous section (whose asymmetric propagation, however, can be analytically quantified), preparing the Gaussian initial data of the present section appears to be considerably more tractable in optical experiments e.g., with lithium niobate waveguide arrays. On the other hand, in this latter setting, we will have to rely on detailed numerical computations of the rectification factor, as this set of initial conditions is less amenable to detailed analytical considerations. The wavepacket considered is given by the equation for starting position index n 0 and width parameter s. We measure the transmission at some value t = t 0 sufficiently large so that the wavepacket has interacted with the nonlinear region and, as a result, some portion of it has been accordingly transmitted through and reflected from the relevant interval. The transmission is then measured by for k > 0 and k < 0, respectively. Then the rectification factor takes the form f = (τ + − τ − )/(τ + + τ − ). The transmission is, of course, equal in both directions (f = 0) in the N = 1 case. In Figure 8 we plot f as a function of γ 1 , γ 2 in the case of N = 2, and as a function of γ 1 , γ 3 with γ 2 fixed in the case of N = 3. In Figure 8 we also show some typical space-t propagation plots for the Gaussian initial data case. Similar to the plane wave solutions case, the rectification factor may be positive or negative as the parameters change. Here the Gaussian in this particular 8 case with k = π/2 has the following property. For parameter values concentrated in the vertical and horizontal bands in the right-side four plots of Fig. 8 more is transmitted if the wave hits the lower γ-value first, in comparison to the wave hitting the larger γ-value first, i.e., encountering the region which is closer to linear is more conducive towards transmission, while encountering the more nonlinear sites at first is more prone to reflection, a feature that seems to be intuitively justified.

V. CONCLUSIONS
In the present work, we considered a lattice setting where embedded in a linear Schrödinger chain was a nonlinear "segment" of the saturable type. Our analytical considerations were focused around plane waves enabling us to analytically compute both the transmittivity and the rectification factor between left-and right-propagating such waves. These features evidence the asymmetric nature of the propagation in a way that is analytically tractable. This asymmetry can be explicitly traced in the nonlinearity of the relevant setup. We also considered the spectral stability of such states in which we observe some effects of low versus high rates of instability. Finally, we considered the asymmetry of propagation of a Gaussian wavepacket. The latter was also clearly evidenced both in the case of a dimer, as well as in that of a trimer, paving the way for the experimental observation of relevant phenomena, such as the enhanced transmission of a wavepacket when encountering a region of increasing, rather than that of a decreasing nonlinear index profile.
There are numerous aspects that may be worthwhile to further explore. In the context of lithium niobate waveguide arrays, it may be relevant to examine settings that involve a genuinely nonlinear lattice but with its central sites bearing a different nonlinearity than the background. This type of "spatial profile" of the nonlinearity coefficient has attracted considerable interest in numerous recent studies as evidenced by the review of [27] and may also be quite experimentally tractable. On the other hand, the vast majority of the present studies on the nonlinearity-induced asymmetry that we are aware of have focused chiefly on one-dimensional configurations. However, it would be both more numerically challenging and also theoretically intriguing to explore scattering of two-dimensional wavepackets from a central (two-dimensional) segment of a lattice which is genuinely nonlinear. Such aspects are currently under investigation and will be reported in future publications.
Finally, we should note that after submission of this manuscript, we were notified of a related work in [28]. While the models analyzed in these works are fairly similar, distinctive features of the present work are (a) that we explored systematically the stability of the extended waves and we studied the outcomes of their dynamical instabilities; and (b) we also considered the effect of the incidence of a Gaussian wave packet. Lastly, (c) although the latter work is restricted to dimers, our considerations here have been provided for general N and examined for N = 1, 2, 3.