Transverse quasi-modes in periodic potentials

We discuss how to find the quasi-modes of the Schrödinger equation when the potential is periodic in time. Our method confirms that the profile of the time-independent (continuous) component of the wavefunction obeys an effective Schrödinger equation where the potential is given by the Kapitza term plus the temporal mean of the original potential, a result originally found by Rahav et al (2003 Phys. Rev. Lett. 91 110404). We then find closed-form expression for the higher order corrections to the quasi-modes, showing how the generic quasi-mode undergoes periodic temporal oscillations and a non-flat phase profile. Validity of our theoretical results is verified against full numerical simulations of the Schrödinger equation. Our findings can be applied both to quantum mechanics and light propagation in the paraxial regime.


Introduction
The interaction of waves with periodic potential is a basic problem in physics. For example, the Rabi model, i.e. the problem of a two-state quantum system interacting with a periodic perturbation [1], is used to describe several fundamental processes, such as the spin dynamics in magnetic fields (used for example in the nuclear magnetic resonance) [2,3], electromagnetic interaction with matter [4], stabilization of repulsive Bose-Einstein condensate [5], directional coupler in integrated optics [6], Josephson junction [7], and many others, all of them featuring the Rabi oscillations [8]. In the last years, the problem of waves in periodic Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. potential beyond the perturbative regime is of strong interest, including optical traps for cold atoms [9,10] and higher harmonic generation in nonlinear optics [11]. For example, the possibility to reduce a full 3D nonlinear Schrödinger equation with inhomogeneous coefficients to a simplified 1D version has been demonstrated [12].
An exact formal solution for the Schrödinger equation in time periodic potential has been first found by Shirley in 1965 [13]. By using the Floquet-Bloch theorem, Shirley showed that the quasi-modes of the Schrödinger equation are eigenfunctions of a time-independent infinite-dimensional matrix, the latter showing non-diagonal terms due to the coupling between different energies coming from the potential modulation. Existence of quasi-modes, i.e. periodic functions replicating themselves at each period, has been proved in infinite Hilbert space some years later [14]. Several approximated approaches, mainly including perturbative series expansions [15], have been exploited in literature. Between the others, we cite the Magnus expansion [16], gauge transformation [17,18] and multiple scale expansion [19].
A counter-intuitive phenomenon related with periodic potential is the confinement of particles in periodic potential encompassing a vanishing temporal average [20]. This effect, originally discovered by Kapitza while studying an inverted forced classical pendulum [21,22], originates from the influence of fast varying scales on slow scales [18]. A time-independent potential proportional to the square of the local (in time and space) force arises [20,23,24], providing a way to confine particles (Paul trap [25]) and waves [10], or even realizing scatteringless potential using non-Hermitian potentials [26].
In this article we discuss a new gauge transformation to describe wave propagation in the presence of a periodic potential, with particular attention to the Kapitza-based confinement [18,20,23,24]. With respect to the approach presented in [18,27], here we employ a different gauge transformation, chosen so that the solutions satisfy certain physical conditions corresponding to a breathing solution in time. First, we show that, in the presence of a nonvanishing average potential, the overall time-independent component of the effective potential is given by the sum of the average plus the effective potential. Then we provide a form for the Kapitza term, account for the influence of any temporal scale on the continuous component. We also provide a set of ordinary differential equation capable to describe every harmonic component of the quasi-mode, including both the amplitude and the phase terms. Finally, we confirmed our theoretical results by using full numerical simulations of the time-dependent Schrödinger equation.

The mathematical model
Let us consider a Schrödinger equation with a time varying potential V(x, t) [13,18] Let us also suppose that V(x, t) is periodic in time with a period T In terms of Fourier series, we can write In other words, the Hamiltonian operator Ĥ depends explicitly on time t, Ĥ =p 2 2m +V(x, t). The goal is to investigate the existence of confined (along x) quasi-modes of (1), defined as time-periodic wave functions with no appreciable mean spreading while evolving along t over temporal intervals of several periods.

Preliminary considerations leading to our gauge transformation
Solutions of (1) take different forms according to the properties of the commutator between Hamiltonians considered at different times [16]. If Hamiltonians at different times commute, the solution is given by [22] ψ( When the adiabatic theorem holds valid (i.e. slow variations in Ĥ ), eigenvectors of equation (4) are eigenstates of the energy, but with an eigenvalue E depending on time t. Conversely, when the Hamiltonians at different time do not commute, wavefunction evolution is governed by a Dyson series [14], implying in general a much more complex coupling between different instantaneous energy eigenstates at different time. The commutator Ĉ for equation (1) Let us first suppose the simplest case when temporal and spatial terms in the Hamiltonian are separated, i.e. V(x, t) = U(x) + W(t). From equation (5) it is Ĉ = 0: the Hamiltonian at different time commute, and the energy eigenstates are thus being ψ n t the energy eigenstate in absence of temporal modulation, i.e. when W(t) = 0 at every instant t. The energy eigenvalue is then E (D) When Ĉ does not vanish, energy eigenstates in the form (6) are no more solution of the original Schrödinger equation: evolution is now governed by a cumbersome Dyson series [14,16]. Problems arise from the fact that spatial shifts (i.e. motion along x) and temporal shifts (i.e. motion in time) do not commute anymore. As in the Magnus approach [28], the idea is to conserve an exponential-like solution by adding a dependence on x on the exponent in equation (6) [17]. Moreover, due to the temporal periodicity of V(x, t), we guess that such exponent, in some complicated way related to the wave energy, is periodic, in agreement with the Bloch-Floquet theorem. Thus, we look for modes in the following form where t 0 is an arbitrary initial time. Each function E n (x,t) represent the nth harmonic. In general, the spatial (x) and the temporal (t) dependence of the functions E n cannot be factorized, for example when the temporal spectrum of the external potential varies in space. Equation (7) corresponds to a gauge transformation. A gauge transformation is also used in [18,27], where the transformation is constructed so that the Hamiltonian is time-independent in the new framework. Our gauge transformation expressed by the ansatz (7) has been devised following a different approach: we chose it so that the solutions directly satisfy some physical properties we expect from the quasi-modes. Our position thus allows to find simpler forms for the oscillating part of the quasi-mode, which is the main aim of the current paper. The basic properties we expect from a quasi-mode are the following: • Solution (7) features an energy varying in time, i.e. i ∂ t ψ/ψ = n=0 E n (x, t).
• When (E n ) = 0 for ∀n, |ψ| 2 = |g(x)| 2 does not vary with the propagation coordinate t, the exponential term affecting only the evolution of the phase distribution. • For (E n ) = 0, the intensity profile |ψ| 2 changes with time. If the functions E n are also time periodic, ψ represents the most generic type of quasi-mode [13], i.e. breathing solutions. Noteworthy, such breathing behavior is inherently different from breathing due to the interference between non-degenerate energy eigenmodes. • When (E n ) = 0 and E n depends on x, the transverse profile of the quasi-mode ψ changes during propagation. • When E n is periodic in time the ansatz (7) is compatible with the Bloch-Floquet theorem.
Let us briefly discuss the orthogonality between two different quasi-modes ψ. Defining the scalar product c ij between ψ i and ψ j reads Even if g i and g j are orthogonal to each other in a given instant, the scalar product c ij will be non-vanishing at later instant t owing to the dependence on x of the exponential.

The master equations
Substituting the ansatz (7) into equation (1), we find the modified eigenvalue problem Equations (9) form an infinite discrete set of equations, corresponding to each harmonic E n of the wave pseudo-energy. The third term on the rhs of equation (9) is nonlinear with respect to E n , thus introducing the coupling between different harmonics E n . An important role is played by the initial time t 0 [29]: to fully address its influence on the wave propagation, we define the functions thus yielding n t t0

Continuous component
Let us start by the CW component, i.e. equation (9) for n = 0. We find Passing from equation (9) to (12), we split the dependence of the definite temporal integrals from the time t and the initial time t 0 by using equations (10) and (11).
In the case of localized quasi-modes, all the secular terms in equation (9) directly proportional to t should vanish. This means that Hence, in agreement with the assumption of waves with average value invariant in propagation, the CW component E 0 does not depend on x. Furthermore, the explicit presence of the imaginary unit in the eigenvalue problem (12) implies, in general, a non-flat phase profile for the mode amplitude at the lowest order g(x). This inhomogeneous phase can be eliminated making the gauge transformation Summarizing, equation (12) can be rewritten as Hence, the problem of finding the profile ϕ(x) is equivalent to the computation of the energy eigenstates for a quantum particle subject to a potential V eff independent from time [18], whose expression is The second term on the rhs of equation (16) is due to the Kapitza effect [21]: all the harmonics E n affects the modal profile ϕ(x), i.e. the continuous component. We also note that, if is a purely real function [20]. According to equation (16) and differently from previous works on the topic [18,20,23,27], in our gauge the time-independent effective potential does not depend only on the external potential V(x, t), but also on how the quasi-mode oscillates in time. We have thus found that a strongly oscillating quasi-mode changes its own effective potential. From a physical point of view, this effective nonlinearity arises from the coupling induced by V(x, t) between different temporal frequencies. Similar effective nonlinearities appear when the Madelung transformation is applied to the Schrödinger equation [30]. When such variations are strong enough to substantially deform the shape of V eff (x), confinement can be lost. In other words, our gauge provides an intuitive and simple interpretation for the absence of localized quasimodes for certain parameters of the external potential V(x, t) [23,25].
If period T is small enough so that terms containing temporal integral in equation (9) can be neglected, we obtain for n = 0 Substitution into equation (16) provides which is the Kapitza potential summed to the CW component of the potential. In fact, assuming with f n independent on x (differently from the general case pictured by V n (x)) and where is assumed to be a real function. Equation (18) then turns into in agreement with [20,23]. Via equations (7) and (14) the quasi-modes are where is the dynamic (time-dependent) part of the potential. This result is a generalization of [17], where the case of a sinusoidal temporal modulations is considered. In this limit the time-modulation of the potential is too fast with respect to the diffraction length of the wave: intensity profile remains unvaried along evolution. The dynamic part of the potential has a two-fold role: it affects the mode profile g(x) via the Kapitza effect and modulates, both in space and time, the phase distribution of the field.

Higher-order harmonics
In the previous section we found the equations obeyed by the continuous component of the quasi-mode. In this section we will instead focus on the behavior of E n (x, t) for |n| > 0, the latter being nonetheless necessary to find the effective potential (18) acting on the CW component. Inspired by equation (17), we will suppose that the pseudo-energies E n can be factorized, i.e. we assume that E n (x, t) = e n (x) exp (−2πint/T). Noteworthy, a quasi-mode breaking the periodic symmetry imposed by the external potential V(x, t) would imply a spontaneous symmetry breaking with respect to the temporal variable, similarly to the recently introduced time crystals [31], but actually accessible solely in the nonlinear regime. Equation (9) then provides It is convenient to recast equation (21) as a function of ϕ(x). Temporal integral in equations (10) and (11) can now be easily computed [18] Using (14) and the definition of I 1 (22), equation (21) finally yields

Pseudo-energy functions
To proceed further, we now assume that V(x, t) = f (t)U(x) for the sake of simplicity. Our approach is to write e n (x) as a power series of the modulation period T [18]. Mathematically we set Substitution of ansatz (25) into (21) provides, up to the order T 2 Accordingly, terms depending on the Planck constant h start to contribute to the effective potential equation (18) at the order O(T 3 ), i.e. quantum contributions to the Kapitza potential becomes important for long enough modulation periods. We would like to stress that, with the gauge transformation used in [27], quantum corrections to the classical Hamiltonian are proportional to T 4 , thus remarking the differences between the two methods. When the temporal average of fn n e −2πin t 0 T = f (t)dt t0 . Thus, given that the temporal integral of f (t) is a periodic function oscillating around zero, it is always possible to choose an initial time t 0 so that the last term in equation (28) vanishes. At the following order in T it is not more possible to set I 1 (x, t 0 ) identically zero across all the axis x. From a physical point of view, the phase surface of each harmonics will be flat for different values of the time t.
Re-arranging we get

First-order pseudo-energies
Let us start by discussing the lowest-order terms e n (x). According to our expansion in powers of the period T, the functions e (1) n (x) share the same profile along x, regardless of the index n. This means that we can set where the spatial profile e (1) (x) is given by and the weight of the nth harmonic is Equation (29), together with equation (32), shows clearly that the functions E n (x,t) share the same transverse profile up to the order T (conversely, the shapes versus x of e (2) n (x) depend on the parameter n, see equation (42) below). According to equation (15), ϕ(x) is a complex function, thus equation (30) states that the first-order correction to the energy is also complex, in general. Function e (1) (x) can be made real using the arbitrariness in the choice of the initial time t 0 . In fact, the phase of g(x) is flat across x if we set Condition (35) cannot be satisfied for every x in the general case. But if we limit to the component e (0) , equation (35) reads that is, the same condition required to cancel out the last term in equation (28).

Second-order pseudo-energies
We now discuss the behavior of the components of e n (x) proportional to the square of the modulation period. To simplify the notation we define Last row of equation (31) now reads If we introduce the x-independent coefficients Equation (39) can be recast as where we defined s n = l =0,n f l f n−l l(n−l) . The terms depending on c (2) n couples different harmonics, resulting in a wave oscillating at different frequencies with respect to the one associated with the external field V(x, t). This effect can be associated with the pseudo-nonlinearity stemming from the gauge transformation, as discussed after equation (16) for V eff . Figure 1 summarizes the behavior of the introduced function versus the transverse coordinate x. For simplicity, in figure 1 the first order for V eff provided by equation (19) is plotted, i.e. the effective nonlinearity is neglected. To fix the ideas, we considered U(x) to be Gaussian-shaped, U(x) = U 0 exp −(x 2 /w 2 0 ) and sinusoidal modulation in time, i.e. f (t) = cos (2πt/T). We also normalized equation (1) by introducing the normalized transverse coordinate x = x/w 0 . Then, by dividing both the sides of equation (1) by , we define the normalized dispersion coefficient D N = / 2mw 2 0 and the normalized potential U N = U/ . Accordingly, the adimensional normalized effective potential is V eff / ( D N ). To describe the spatial behavior of e (1) (x) and e (2) n (x), we introduce the normalized function

Pseudo-energies for a Gaussian-shaped potential
The potential V eff is M-shaped, presenting a potential well around x = 0. Nonetheless, given that the value of the effective potential is the same in x = 0 and for |x| → ∞, only leaky modes (i.e. finite lifetime states in quantum mechanics) are supported [32]. The fundamental mode (i.e. bell-shaped and presenting the lowest attenuation coefficient) can be found by considering only the central lobe of the potential [23], i.e. considering for |x| > x max a constant potential of value equal to the maximum of V eff , where x max is the coordinate corresponding to the maximum of V eff . Solutions for ϕ plotted in figure 1 are computed numerically after discretization the Schrödinger equation and then solving the eigenvalue problem of the found finite matrix. The attenuation coefficient of the two modes plotted in figure 1 is very small, given that the spatial overlap between the mode and the lateral lobes is negligible [32]. Due to a higher potential barrier, the quasi-mode gets narrower in the transverse direction as the temporal period T increases. At the same time, the amplitude of both the functions M(x ) and F N (x ) grows up, in turn corresponding to a greater breathing for the quasi-mode in each single period. The nature of this time-dependent modulation on the wavefunction profile will be discussed in detail in the next section.

Wavefunction
According to equation (7), the wavefunction of the quasi-mode is The potential minimum U 0 is fixed and equal to −1. The quasi-mode is computed by considering only the central lobe of V eff , see [23].
The component D (x, t)dt corresponds to the firstorder solution given by equation (20). Here, we are interested in the form assumed by as the period T changes. Due to the definitions (44) and (45), the complete wavefunction is ψ = ψ (0) ψ (1) ψ (2) .

First-order correction.
First point to address is the behavior of the imaginary part in the exponential terms in equations (44) and (45): in fact, a real exponent means a change in the intensity profile of the wave. Regarding equation (44), from equation (34) it is straightforward , thus providing where we used the property ∠a n = ∠f n + π stemming directly from equation (34). Equation (46) states that, if the imaginary part of e (1) (x) vanishes, only a periodic variation on the wave amplitude takes place, without affecting the phase distribution. Computing the integral we obtain Summarizing, if we develop E n (x,t) up to the linear terms in T, quasi-modes are given by where By definition P(t) is a periodic function featuring a null temporal average: thus, in this limit P(t) and f (t) feature the same discrete spectrum, but with the nth harmonic being modulated by a factor n −2 , plus an overall temporal shift of π. Stated otherwise, the modulus of the field follows adiabatically the external potential V(x, t) due to the smallness of the dispersion length with respect to the modulation period. Accordingly, the function ϕ(x) represent the averaged value of the intensity profile, whereas R(x, t) accounts for periodic fluctuations in the intensity distribution due to the time-variable potential. The sign of the function M(x) changes across the spatial transverse dimension x (see figure 1), thus resulting in R(x, t) larger or smaller than one across the wave cross-section for a given time. Thus, equation (48) correctly predicts a transverse motion of energy from the edges (center) to the center (edges) during the temporal intervals when the instantaneous potential V(x, t) is attracting (repelling). From equation (48) we can evince the phase delay between the oscillations of the phase term e − i V (0) D (x,t)dt = e iφ (0) (x,t) and of the periodic amplitude modulation given by R(x, t). For the sake of simplicity, we consider a sinusoidal modulation in time by setting f (t) = cos 2πt T . Integrating we find that Indeed, φ (0) represents the curved wavefront oscillating in time and providing the energy flux described above: the amplitude R(x, t) oscillates in phase with the external potential V(x, t), whereas the phase front curvature given by φ (0) is in advance of π/2. The wave evolution in one single period for a bell-shaped U(x) is depicted in figure 2.
The wave is the narrowest in t = 0, corresponding to the minimum in time of the potential (i.e. maximum attractive force). In t = T/4 and t = 3T/4 the instantaneous potential is vanishing and R(x, t) = 1, thus the wave profile corresponds to the CW component computed from equation (15), i.e. ψ = ϕ. The wavefront in these points reaches its maximum warping in absolute value, that is, the field variations in time are the fastest. At t = T/2, thus in correspondence to the maximum values in time of the potential V(x, t), the quasi-modes reach the maximum width. The curvature of the wavefront is vanishing in t = 0 and t = T/2, where the absolute value of the potential achieves its maximum. Figure 2 also shows the typical evolution of the periodic wavefunction given by equation (48) versus the modulation period T. For short periods, variations in the intensity profile in propagation are negligible. For longer periods T, the quasi-mode becomes narrower, increasing the tendency of the wave to diffract. Owing also to the shorter modulation period, the wave profile now undergoes a strong breathing in propagation. Thus, the width oscillation of the mode increases as the modulation period becomes longer. In the plotted example, for T = 60 µs light is bell-shaped everywhere. For longer periods, more complex profiles, encompassing a dip in the center, appear when the field achieves the maximum transverse size. To understand if this is an artifact coming from our approximation, middle row of figure 2 plots the power (or particle density for matter waves) P-defined as P = |ψ| 2 dx-versus the propagation time t. The power changes in propagation, and the amplitude of the oscillation increases with the modulation period T. In fact, equations (48)-(50) hold valid only in the limit of small T. Accordingly, the power trend shows a periodicity of T/2, implying the presence of higher-order terms in the expansion of e n (x) given by equation (25). Eventually, our approximation ceases to be valid for T = 70 µs: this result will be later confirmed using full numerical simulations of equation (1).

Second-order correction.
Let us now turn discuss the wave component ψ (2) . We start by addressing the symmetry properties of e (2) n (x) with respect to the transformation n → −n. First term on the rhs of equation (42) obeys the relationships b ; the same symmetry is possessed by the second term c (2) n . Calling ψ (2) bn and ψ (2) cn the terms depending on b (2) n and c (2) n respectively, equation (45) provides Equations (51) and (52) contribute both to the overall quasi-mode as a pure phase factor, thus not affecting the amplitude of the field.

Higher-order corrections on the effective potential
We now substitute higher order terms of the pseudo-energies E n (x,t) into the Kapitza-like potential V eff given by equation (16). Halting to the terms proportional to T 4 we find (53) In equation (53) the terms proportional to T 3 and T 4 are of quantum origin due to their dependence on , differently from the term ∝ T 2 , indeed corresponding to the Kapitza potential in the classical case [21]. Additionally, the terms starting from T 3 contain the CW profile ϕ itself, i.e. in our model the effective potential becomes nonlinear. The effective nonlinearity depends on log g 2 , thus on transverse width of the profile (it is log ϕ 2 ≈ −2x 2 /w 2 being w the width of ϕ, the former expression being exact if ϕ(x) is a fundamental Gaussian mode): the wave propagation does not depend on the amplitude of the wave itself, in accordance with the absence of nonlinearity in equation (1).

Comparison with numerical simulation in the case of sinusoidal potential
The simplest case of periodic potential is a sinusoidal modulation, f (t) = cos(2πt/T) and f n = 0.5 [δ n,1 + δ n, −1 ] in the real and in the transformed domain, respectively. We note that s n does not vanish only for |n| = 2, thus the last term in equation (53) is null. The effective potential (53) now reads Behavior of the effective potential V eff versus x for three values of the period T is shown in figure 3. First, we notice that the higher-order terms feature a transverse shape more complex (i.e. more lobes) than the dominant term. For short periods (T = 10 µs), the contributions proportional to T 3 and T 4 are negligible. For intermediate values (T = 40 µs), the higher-order contributions are about 5% of the dominant term. Despite that, their effect consists solely in a slight widening of the central lobe. For larger periods (T = 70 µs) a bigger distortion on the central lobe appear, inducing an appreciable change in the shape of the potential around |x/w 0 | ≈ 1.

Numerical simulations
We now check the accuracy of our theoretical results by direct comparison with numerical simulations of equation (1) in the case of a sinusoidal modulation in time by setting f (t) = cos 2πt T + α 0 . The quasi-mode ϕ calculated from equation (15) corresponds to the intensity profile of the wave in 2πt T + α 0 = π 2 . Thus, in the simulations we fixed α 0 = π/2 and, in agreement with equation (48), we used ϕ times e −iφ (0) as the initial condition in t = 0. Figure 4 shows the behavior of the wave over several periods, for different values of the modulation period and a Gaussian shaped potential featuring U 0 = −1. The quasi-mode computed from equation (17) propagates almost unvaried for short modulation periods. As period increases, the quasi-mode becomes narrower, and at the same time the breathing oscillations increase in amplitude. Comparison with the theoretical predictions (middle and bottom rows) confirm a good agreement between theory and numerical simulations for not too long periods T: for T = 70 µs the diffraction becomes too large, and the quasi-mode rapidly decays in propagation due to the large losses. This is in good quantitative agreement with figure 3, where the overall effective potential is appreciably distorted when T = 70 µs. Figure 5 compares the solutions on longer intervals. We first compute the maximum for each instant t for both the approaches. We then plot the behavior of the maxima versus the time t for four different periods T. The matching is almost perfect up to T = 30 µs. Appreciable discrepancies (about 10% with respect to the amplitude of the field; the error in the oscillating component varies between 30% and 100%) appear when T = 50 µs. For longer periods, the distance between  the two curves rapidly increases due to the O(T 4 ) terms neglected in the theory. For T = 60 µs the field distribution remains periodic, despite an appreciable amount of losses at the interface due to the relevance of the higher order terms.

Conclusions
In this paper we discussed the propagation of waves obeying the Schrödinger equation in a periodic potential. We started by making a special ansatz for the wave shape driven by physical considerations about the property of the system. Using this ansatz, we first demonstrated a general formula for the Kapitza potential, including the case when the temporal average of the potential is non-vanishing. We then analytically computed the higher-order harmonics of the field as a power series of the modulation period T, including both amplitude and phase terms. These higher orders allow also the computation of the Kaptiza potential as a power expansion of the modulation period [18], thus addressing the validity range of the model. For bell-shaped potential wells, our results predict a quasi-invariant profile for the fundamental quasi-mode when the modulation period is small. For longer periods, the Kapitza potential becomes stronger, but at the same the temporal oscillations of the wavefunction grow up. When the period reaches a threshold value, quasi-confinement is no more observed. Our results have been confirmed by direct numerical simulations, the latter confirming also the dependence from the modulation period of the accuracy of the theoretical model. These findings can find applications in several fields, ranging from light confinement of particles (both in the classical [33] and quantum regime [34]), to periodic gratings and segmented waveguides in integrated optics [35,36], Floquet periodic structures [37,38], nonlinear periodic systems [39] and highharmonic generation in nonlinear optics [11].

Acknowledgments
The Author thanks DFG for funding (GRK 2101) and Dr C P Jisha for the careful reading of the manuscript.