Perturbation Theory for PT-Symmetric Sinusoidal Optical Lattices at the Symmetry-Breaking Threshold

The $PT$ symmetric potential $V_0[\cos(2\pi x/a)+i\lambda\sin(2\pi x/a)]$ has a completely real spectrum for $\lambda\le 1$, and begins to develop complex eigenvalues for $\lambda>1$. At the symmetry-breaking threshold $\lambda=1$ some of the eigenvectors become degenerate, giving rise to a Jordan-block structure for each degenerate eigenvector. In general this is expected to give rise to a secular growth in the amplitude of the wave. However, it has been shown in a recent paper by Longhi, by numerical simulation and by the use of perturbation theory, that for an initial wave packet this growth is suppressed, giving instead a constant maximum amplitude. We revisit this problem by developing the perturbation theory further. We verify that the results found by Longhi persist to second order, and with different input wave packets we are able to see the seeds in perturbation theory of the phenomenon of birefringence first discovered by Makris et al.


I. INTRODUCTION
The study of quantum mechanical Hamiltonians that are P T -symmetric but not Hermitian [1]- [6] has recently found an unexpected application in classical optics [7]- [15], due to the fact that in the paraxial approximation the equation of propagation of an electromagnetic wave in a medium is formally identical to the Schrödinger equation, but with different interpretations for the symbols appearing therein. It turns out that propagation through such a medium exhibits many new and interesting properties, such as power oscillations and birefringence. The equation of propagation takes the form where ψ(x, z) represents the envelope function of the amplitude of the electric field, z is a scaled propagation distance, and V (x) is the optical potential, proportional to the variation in the refractive index of the material through which the wave is passing. A complex V corresponds to a complex refractive index, whose imaginary part represents either loss or gain. In principle the loss and gain regions can be carefully configured so that V is P T symmetric, that is V * (x) = V (−x). There is also a non-linear version of this equation, arising from sufficiently intense beams, where there is an additional term proportional to |ψ| 2 ψ. However, for the purposes of this paper we shall limit ourselves to the linear case.
A model system exemplifying some of the novel features of beam propagation in P T -symmetric optical lattices uses the sinusoidal potential This model has been studied numerically and theoretically in Refs. [9,12,13]. The propagation in z of the amplitude ψ(x, z) is governed by the analogue Schrödinger equation1, which for an eigenstate of H, with eigenvalue β and z-dependence ψ ∝ e −iβz reduces to the eigenvalue equation It turns out that these eigenvalues are real for λ ≤ 1, which corresponds to unbroken P T symmetry, where the eigenfunctions respect the (anti-linear) symmetry of the Hamiltonian. Above λ = 1 complex eigenvalues begin to appear, and indeed above λ ≈ 1.77687 all the eigenvalues are complex [16]. Clearly one would expect oscillatory behaviour of the amplitude below the threshold at λ = 1 and exponential behaviour above the threshold, but the precise form of the evolution at λ = 1 is less obvious. At first sight one would expect linear growth because of the appearance of Jordan blocks associated with the degenerate eigenvalues that merge at that value of λ, but, as Longhi [12] has emphasized, this behaviour can be significantly modified depending on the nature of the initial wave packet.
In a previous paper [17] we approached this problem by explicitly constructing the Bloch wave-functions and the associated Jordan functions corresponding to the degenerate eigenvalues and then using the method of stationary states to construct the z-dependence. We found that the explicit linear dependence arising from the Jordan associated functions is indeed cancelled by the combined contributions from the non-degenerate wave-functions and were able to understand how this cancellation came about. In the present paper we approach the problem from a different point of view by revisiting the complementary perturbative calculation of Longhi [12]. In Section 2 we briefly recapitulate how the spectrum and eigenfunctions are calculated. Then in Section 3, which forms the main body of the paper, we give an explicit expression for the first-order contribution and carry out the second-order calculation in detail. This enables us to investigate the saturation phenomenon for a variety of different inputs. Finally, in Section 4 we give a brief discussion of our results.

II. BAND STRUCTURE AT THRESHOLD
At the threshold λ = 1, the potential V in Eq. (2) becomes the complex exponential This is a form of the Bessel equation, as is seen by the substitution y = y 0 exp(iπx/a), where where q 2 = β(a/π) 2 . Thus the spectrum is that of a free massive particle, shown in the reduced zone scheme in Fig. 1, and for q ≡ ka/π not an integer the solutions ψ k (x) = I q (y) and ψ −k (x) = I −q (y) are linearly independent, and have exactly the correct periodicity, ψ k (x + a) = e ika ψ k (x), to be the Bloch wave-functions. It is important to note, however, that because the original potential is P T -symmetric rather than Hermitian, these functions are not orthogonal in the usual sense, but rather with respect to the P T inner product, namely However, for q = n, a non-zero integer, I n (y) and I −n (y) are no longer independent. In that case the Bloch eigenfunctions do not form a complete set, and we must search for other functions, still with the same periodicity, to supplement them. These are the Jordan associated functions, which we denote by ϕ k (x) ≡ χ n (y). They may be defined as derivatives of the eigenfunctions with respect to β, and satisfy the generalized eigenvalue equation The crucial feature of the Jordan functions is that because of this latter equation they naturally give rise to linear growth in z, provided that they are excited: However, as was found numerically in Ref. [12], and explored further in Ref. [17], this natural linear growth may become saturated due to the contributions of neighbouring Bloch functions, which are closely correlated with those of the Jordan functions.

III. PERTURBATION THEORY
The analysis of Ref. [17] approached the problem from one point of view, in which the interplay between the contributions of the Bloch eigenfunctions and the Jordan associated functions was made explicit. A complementary way of looking at things, which does not separate these two contributions, is to use the perturbative expansion, which instead emphasizes the contributions of the free propagation and the corrections brought about by the potential. The general framework for an expansion of ψ(x, z) in powers of V 0 , namely ψ(x, z) = ∞ r=0 V r 0 ψ r (x, z), has been given in Ref. [12], along with an approximate form of the first-order term ψ 1 (x, z) for the case q 0 = −1 and w large. In this section we generalize this calculation by obtaining analytic expressions for both the first-and second-order terms for general q 0 and w. Of course, this can only be used as a guide because there is no guarantee that the expansion converges, nor that the large-z behaviour of the complete amplitude can be extracted from the behaviour of the truncated series. We will take as our input a Gaussian profile of the form with offset k 0 and width w. The zeroth-order term, ψ 0 (x, z), is just the freely-propagating wave-packet while the first-order term, ψ 1 (x, z), is given by wheref is the Fourier transform of f (x) of Eq. (8) and k B = 2π/a is the width of the first Brillouin zone. We can reverse the order of integration in the expression for ψ 1 , performing the (Gaussian) k integration first, to obtain The y integration is then also a Gaussian integration over a finite range, giving the result where For the purposes of considering large w, it is convenient to rewrite these in the form The case k 0 = − 1 2 k B , i.e. q 0 = −1, is clearly very special, since in this case the second terms in the contributions to η 1 and η 2 vanish, so that we get the simple expressions In that case, as long as w 2 ≫ 4z the arguments may be treated as effectively real, and each error function behaves like a sign function of its argument (see Fig. 2(a)), so that the sum of the two behaves like the step function θ(k B z−|x|). This is the function Φ(x/(k B z)) of Ref. [12]. In this case the qualitative features of the perturbative calculation are in complete agreement with the spreading of the wave-function in Fig. 3 of that paper, and the saturation 1 of ψ max . However, in this treatment there is of course no mention of whether or not any Jordan functions are excited. In (a) we plot erf(4 + y/w) + erf(4 − y/w) with w = 80. In (b) we plot w e −w 2 |erf(4 + y/w − iw) + erf(4 − y/w + iw)| The same expressions in Eqs. (13) and (15) can also be used for the cases q 0 = 0 and q 0 = 1 in the limit of large w. In each case the arguments η 1 and η 2 now have a large imaginary part. In that situation the modulus of the erf has a narrow peak where the real part vanishes (see Fig. 2(b)). Thus the result consists of two narrow rays, which are centered on x = 2k 0 z and x = 2(k B + k 0 )z. Here we have the seeds of the birefringence first observed in Ref. [9]. For the case k 0 = q 0 = 0, the two rays are centered on x = 0 and x = 2k B z, while for the case q 0 = 1, or k 0 = 1 2 k B , the two rays are centered on x = k B z and x = 3k B z. We now go on to second-order perturbation theory to investigate the behaviour of ψ 2 (x, z), which is given by [12] Again the k integration is a Gaussian, which leaves finite-range Gaussian integrations over ξ and η. It is convenient to change the integration variable ξ to y ≡ 2η − ξ. The integration over y then yields the expression where we have written k 0 = −( 1 2 k B + ∆), and a and b are given by .
The final η integrations are then of the form Thus in principle ψ 2 is expressible in terms of eight error functions. However, it turns out in practice that only six are involved.
In the case k 0 = −k B /2, or q 0 = −1, when ∆ = 0, the arguments of the error functions are such as to give a plateau in |ψ 2 | between x = −3k B z and x = −k B z, i.e. a widening of the beam to the left of ψ 0 , and a much smaller peak, centered around x = 3k B z, that is, a second weak beam to the right. More importantly, the second-order contribution again shows no sign of the linear growth 2 in z naïvely expected from the excitation of Jordan associated functions.
In the other case, q 0 = 0, corresponding to ∆ = −k B /2, there are three peaks, centered around x = 0, x = −2k B z and x = 4k B z, representing a further splitting of the initial beam. Both cases are illustrated in Figure 3.

IV. DISCUSSION
Of course one needs to treat the results of perturbation theory with caution. In general terms we have no proof that the perturbation series converges, and in particular the asymptotic behaviour in z of a few terms of the series does not necessarily give the correct asymptotic behaviour of the entire sum. Nonetheless this series does appear to give reliable results. For the parameters used by Longhi in Ref. [12], with V 0 = 0.2, first-order perturbation theory already reproduces the numerical results very well, and the second-order term gives only an extremely small correction.
In all of our calculations we have not allowed z to become too large, restricting it by the condition z ≪ w 2 /4, in which case saturation is an inbuilt feature of perturbation theory. In fact it was shown in Ref. [17] using the method of stationary states that if one goes to much larger values of z the amplitude so calculated begins to grow again but this ultimate resumption of linear growth is not physical, because it corresponds to the situation where the beam has widened beyond lateral limits of the optical lattice.
It is an elegant feature of the perturbative expansion that the different types of possible behaviour of the beam -spreading or splitting into two or more beams -arise from very simple properties of the error function depending crucially on the offset k 0 . In the first case the arguments are essentially real, and the error functions behave like sign functions, while in the second case there is a large imaginary part, and the moduli of the error functions behave instead like narrow peaks.