{'fraction_non_alphanumeric': 0.05859444149406975, 'fraction_numerical': 0.022747388918392637, 'mean_word_length': 3.692068106312292, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lyapounov exponent and density of states of a one-dimensional non-Hermitian Schroedinger equation

We calculate, using numerical methods, the Lyapounov exponent gamma(E) and the density of states rho(E) at energy E of a one-dimensional non-Hermitian Schroedinger equation with off-diagonal disorder. For the particular case we consider, both gamma(E) and rho(E) depend only on the modulus of E. We find a pronounced maximum of rho(|E|) at energy E=2/sqrt(3), which seems to be linked to the fixed point structure of an associated random map. We show how the density of states rho(E) can be expanded in powers of E. We find rho(|E|) = 1/pi^2 + 4/(3 pi^3) |E|^2 + ... This expansion, which seems to be asymptotic, can be carried out to an arbitrarily high order.


Introduction
It has been realised over the last few years that a number of non-equilibrium problems can be described through the time evolution of a non-Hermitian random Hamiltonian. Some noteworthy applications include the motion of a particle in an imaginary vector potential [1], vortex line pinning in superconductors [2,3], and growth models in population biology [4]. The interest for non-Hermitian random Schrödinger equations increased further [5,6,7,8,9,10,11,12] as it was realised that, at least in some cases, a transition from localised to delocalised states could take place even in one dimension.
The simplest examples considered were one-dimensional tight-binding models for which the wave function satisfies e h ψ n+1 + e −h ψ n−1 + V n ψ n = Eψ n , (1.1) where V n is a random potential at site n and E is the energy. For real h and periodic boundary conditions (n ≡ N + n), it was shown [1,5] that the eigenvalues are located (for a large system size) either on the real axis, or on lines which can be understood from the knowledge of the Lyapounov exponent for h = 0.
Here we consider another simple one-dimensional non-Hermitian Schrödinger equation, first introduced by Feinberg and Zee [3]: e iθn ψ n+1 + e iχn ψ n−1 = Eψ n . (1.2) In this case there is no random potential V n , and the disorder is purely off-diagonal. Both θ n and χ n are uniformly distributed between 0 and 2π, and they are all independent.
A numerical study [3] of the spectrum of this Schrödinger equation suggested that the density of states was roughly uniform in a circle of radius ∼ π/2. We found it challenging to see whether existing tools [13,14,15,16] could be used to calculate this spectrum.
In the present work we mostly study the density of states of Eq. (1.2) by calculating the associated Lyapounov exponent and using the Thouless formula [17] (which remains valid for non-Hermitian one-dimensional models as shown in Section 2).
The calculation of the Lyapounov exponent, and consequently the density of states, can be done either by Monte Carlo methods or by solving numerically an integral equation for the Ricatti variable r n = |ψ n /ψ n−1 |. This is done in Section 3, where we also discuss the possibility of singularities both in the probability distribution of r n and in the density of eigenvalues E of Eq. (1.2).
In Section 4 we develop a method to perturbatively calculate the Lyapounov exponent and the density of states, in powers of the energy E. This approach is an extension of similar calculations done in the past, where the anomalous behaviour of the density of states at the band centre was obtained in the Hermitian case [15,16]. Our expansion indicates that the radius of convergence is E = 0 and that the series is asymptotic.

The Lyapounov exponent and the density of states
An old and powerful way [13,14]  can be recast as where M n is a random 2 × 2 matrix so that for any choice of the energy E and of the wave function ψ 0 and ψ 1 at sites 0 and 1, the wave function ψ n for n ≥ 2 reads If ψ 0 and ψ 1 are chosen arbitrarily (i.e. they are fixed, and not tuned functions of E and of all the θ p and χ p ), one expects that for large N where γ(E) is the largest Lyapounov exponent of the product of random matrices (2.3).
Because the matrices M n are invariant under the transformation M n (E, θ n , χ n ) → M n (Ee iβ , θ n + β, χ n + β) for any β, and because by shifting all the phases θ n and χ n by a constant β one gets an equally likely random sample, it is clear that γ(E) = γ(Ee iβ ).
Thus the Lyapounov exponent depends only on the modulus of E. For similar reasons, the average density ρ(E) of eigenvalues is invariant under the change E → Ee iβ : In fact these two quantities are related by the Thouless formula [17] γ which can be derived as in the Hermitian case by the following argument: Consider a finite system of N sites with Dirichlet boundary conditions (ψ 0 = ψ N +1 = 0). If one fixes ψ 0 = 0 and ψ 1 = 1, the ψ N +1 calculated from the recursion (2.1) is found to be a polynomial of degree N in E, the zeroes of which are the N eigenvalues E α (for Dirichlet boundary conditions): Applying Eq. (2.4) for large N leads to Eq. (2.6).
It turns out that for all E = 0, the Lyapounov exponent γ(E) is strictly positive by the Furstenberg theorem [18,19]. (This will in fact be confirmed by our numerical results, and by the small-E expansion of Section 4.) One can then argue that the density of states is independent of the boundary conditions. Indeed, for general E, Dirichlet (ψ N +1 (E) = ψ 0 (E) = 0) and periodic boundary conditions (ψ N +1 (E) = ψ 1 (E)) are very or, exploiting the rotational symmetry of γ(E) and of ρ(E), that ρ(|E|) = 1 2π For |E| > 2 the density of states vanishes. To see this, consider an eigenfunction ψ n corresponding to an eigenvalue E α . By applying the triangular inequality to Eq. (1.2), one obtains |E α ψ n | ≤ |ψ n−1 | + |ψ n+1 |. Choosing n such that |ψ n | = max (|ψ i |) N i=1 , this implies |E α | ≤ 2. As the density of states ρ(E) vanishes for |E| > 2, by expanding the logarithm in Eq. (2.6) and by using the rotational symmetry (2.5) of ρ(E) one readily finds that (2.10) All our calculations which follow are based on another reformulation of Eq. (1.2).
Introducing the Ricatti variable r n r n = ψ n ψ n−1 (2.11) the problem (2.3) takes on the form of iterating a random map shown in Figure 1 r n+1 = 1 r n + Ee iϕn , (2.12) where ϕ n (which is equal to χ n − θ n + π modulo 2π) is uniformly distributed between 0 and 2π. The Lyapounov exponent is then given by For large n the probability distribution of r n becomes independent of n, and the invariant distribution P (r, E) satisfies (2.14) From the knowledge of this invariant distribution, the calculation of the Lyapounov exponent γ(E) follows from Eq. (2.13), and one has whereas the density of states ρ(E) is given by Eq. (2.9).
All the rest of the paper is devoted to the determination of the invariant distribution P (r, E) solving Eq. (2.14), and to the calculation of the Lyapounov exponent γ(E) as well as the density of states ρ(E).

Numerical study
In this section we numerically determine the invariant distribution P (r, E) solution of Eq. (2.14), the Lyapounov exponent γ(E), and the density of states ρ(E).
According to Eq. (2.9) the computation of the density of states requires both γ(E) and its first two derivatives. Whilst γ(E) is rather easy to determine numerically, a precise estimate of its derivatives with respect to the energy turns out to be far more difficult. We have attacked the problem by two different methods.

The Monte Carlo approach
This first method consists in iterating Eq. (2.12) for a large sample. Typical shapes obtained for the stationary distribution P (r, E) are shown in Figure 2.
First, one notices that depending on the value of E, the support of the invariant measure P (r, E) is either finite of infinite. If one assumes that the support is an interval [r min , r max ], we see from Eq. (2.12) that if r n ∈ [r min , r max ], the requirement that also r n+1 ∈ [r min , r max ] gives conditions that r min and r max should satisfy. After analysing all the possible cases, one finds that the support of the distribution P (r, E) is the whole positive real axis for E < 2, whereas for E > 2 the support is finite, with r min and r max satisfying r min = E − 1/r min and r max = E + 1/r min , that is Second, for large enough E, there are apparent singularities at certain values of r.
These singularities seem to be remarkable points of the maps T + and T − obtained from Eq. (2.12) by choosing ϕ n = 0 or ϕ n = π, i.e.
T + (r) = E + 1 r , For each value r n , it is easy to see that T + (r n ) and T − (r n ) are the two extreme values that r n+1 can take:  iterations.
The fixed point of T + coincides with the most visible singularity of P (r, E) in Figure 2, whereas the iterate T − (r + ) gives the position of another very visible singularity of P (r, E). Other apparent singularities seem to be located at the points T 2 − (r + ) and T + (T − (r + )).
In Appendix A, we show that the singularities in P (r, E) are, if at all present, much weaker than they look on Figure 2. This is confirmed by Figure 3, which is an enlargement of Figure 2.c, where there is a clear rounding of the singularity at r + .
The calculation of the Lyapounov exponent is straightforward by Monte Carlo methods. One just iterates the random recursion (2.12) a large number N of times (typically N = 10 8 ) and one obtains γ(E) by where N ′ is a large enough number of iterations to eliminate the transient effects due to the initial choice of r 0 (here we took N ′ = N/10 which is exceedingly sufficient).
For the density of states one needs Eq. (2.9) to calculate the first two derivatives of To avoid this difficulty we tried to iterate directly the derivatives of r n with respect to E. The recursion (2.12) has the form Now, if r (1) n and r (2) n denote the first and the second derivatives of r n with respect to E, one can obtain (in a complicated form which we do not reproduce here because we wrote it in our programmes through a change of variables) recursion formulae for these derivatives r (1) Iterating these recursions N times yields an estimate for the first two derivatives of γ(E) Whilst the results for γ(E) seem to be alright, the density of states exhibits some irregularities. We tried to understand the origin of these irregularities by changing the length N of the sample, the generator of random numbers, and the way in which Eqs.
sign which more or less compensates the first big contribution. These events (where r n is very small) have a dramatic effect on the accuracy of our results, with the unfortunate consequence that the more we iterate, the more of these events we see, so that the more irregularities are visible.
We unfortunately could not find a satisfactory way of avoiding these difficulties, and so we tried a completely different approach discussed in the next subsection.
In Figure 4.b, we see that the density of states ρ(E) exhibits a non-trivial structure with a pronounced maximum at an energy E ≃ 1.15. Trying to understand the origin of this maximum, we realised that at energy the point T − (r + ) becomes a fixed point of T − . In a more physical language, if one considers an infinite sample for which θ n = 0 and χ n = π for n < 0, θ n = 0 and χ n = 0 for n = 0, (3.10) θ n = π and χ n = 0 for n > 0, one finds that there is a bound state at energy E = 2/ √ 3 corresponding to a wave function ψ n ∝ 3 −|n|/2 in addition to the continuous part of the spectrum at imaginary We believe that the origin of the sharp maximum in ρ(E) at energy E = 2/ √ 3 is due to the existence of these bound states, as in the case of impurity bands for Hermitian problems [19,20]. We did not succeed, however, to determine whether the density of states is analytic at this special energy E = 2/ √ 3, or whether it presents some weak non-analyticity.
Our numerical results for γ(E) in Figure 4.a can be fitted for small E by an expression of the form with k 1 = 0.1595 ± 0.0005 and k 2 = 0.017 ± 0.001. To be more precise, we first find k 1 as the extrapolated interception of γ(E)/E 2 with the ordinate. Then, subtracting the k 1 E 2 term, k 2 can be similarly found from the residual E 4 dependence. Needless to say, this procedure greatly enhances numerical noise at each stage, and so we had to content ourselves finding the first two terms. These are in good agreement with the results k 1 = 1 2π ≃ 0.1592 and k 2 = 1 6π 2 ≃ 0.01689 of the perturbation theory to be developed in Section 4 (see Eq. (4.25)).

Discretisation of the integral equation for P (r, E)
To avoid the problems of the Monte Carlo method, we tried to numerically solve the integral equation (2.14) for P (r, E). To this end the r-variable was discretised by taking 1000 values of r at the points r k = 2k/(2001 − 2k) for 1 ≤ k ≤ 1000. The integral operator (2.14) then becomes a 1000 × 1000 matrix, and the stationary distribution can be obtained by simply iterating a linear system.  We see that both γ(E) and ρ(E) have essentially the same shapes as we have previously encountered in figure 4, using the Monte Carlo method. The irregularities which were visible in Figure 4.b have now disappeared in Figure 6.b. Moreover, the inset in Figure 6.a, showing the asymptotic approach towards the exact result (2.10), indicates that the accuracy of γ(E) has been improved by roughly two orders of magnitude.
In the insets of Figures 4.a and 6.a, the difference γ(E) − log(E) seems to vanish at E ≃ 1.8. We believe that this difference is non-zero up to E = 2 but that it becomes very small (ρ(E) has an essential singularity: it vanishes and all its derivatives vanish at E = 2). The argument that eigenstates exist up to E = 2 can be borrowed from the theory of Lifshitz tails [20,21] in the Hermitian case: in a random sample, one can find arbitrarily large regions where θ n ≃ χ n ≃ 0.

Perturbation theory
In this section we develop a perturbation theory in powers of E to determine the invariant distribution P (r, E) solution of Eq. (2.14). When this is known, we can obtain It is useful to analyse first the small-E expansion of the integral operator which appears in Eq. (2.14). Consider two distributionsQ(r) and Q(r) related bỹ When E < r, this can be rewritten (by integrating over s) as Assuming E ≪ r the integrand can then be expanded as a power series in E, and after performing the integral over ϕ we find that We see that the effect of a non-zero energy is thatQ(r) is not just trivially related to Q 1 r , but also depends on all the derivatives of Q. Quite remarkably, the above expansion can be written in a compact form valid to all orders. Defining the second It is straightforward to see that this agrees with the naïve expansion formula (4.5). The validity of Eq. (4.7) to an arbitrary order in E is established in Appendix B. Note that a simple change of variable r → 1/r transforms Eq. (4.7) intõ E 2n 4 n (n!) 2 L n r 2 Q(r) , (4.8) where the operator L is defined by Coming back to the invariant distribution P (r, E), we see by combining Eqs. (4.7) and (4.8) that it satisfies If we look for a solution P (r, E) which can be expanded as P (r, E) = P 0 (r) + E 2 4 P 1 (r) + E 2 4 2 P 2 (r) + · · · , (4.11) we obtain by equating the two sides of Eq. (4.10), order by order in E, that The term with n = 0 on the right-hand side is just P m+1 (r), same as the left-hand side, so that P m+1 (r) disappears from the equation. Therefore Moving all but the two n = 1 terms, which are the only ones involving P m (r), to the right-hand side we arrive at the equation (4.14) Note that an M operator has been factored out (by using the fact that Lr 2 = r 2 Mr 4 ) from all the terms. Integrating over this, and shifting the summation variables, we finally where A m and B m are a priori arbitrary constants.
When m = 0 we find P 0 (r) = A 0 r log r 1 + r 4 + B 0 r 1 + r 4 . (4.16) By inserting this into Eq. (4.7), with Q(r) andQ(r) replaced by P (r, E), we see that for P 0 (r) to satisfy this equation to zeroth order in E, i.e. P 0 (r) = P 0 (1/r)/r 2 , the constant A 0 should vanish, whereas B 0 could be arbitrary (the problem is linear and B 0 would be fixed by the normalisation of P (r, E)). This leads to P 0 (r) = B 0 r/ (1 + r 4 ). for m > 0, and B 0 is an arbitrary function of E (which can be fixed by requiring that the distribution P (r, E) be normalised).
P 0 = f 1 , In order to compute the contributions to γ(E) and to the normalisation, one needs the following integrals Using (4.18), these are easily evaluated Combining all of this we arrive at the expansion (4.25) In this way we can in principle generate γ(E) to any order. We calculated numerically the coefficients a 2n of the expansion of γ(E) γ(E) = a 2 E 2 + a 4 E 4 + · · · + a 2n E 2n + · · · up to n = 46 in order to estimate the radius of convergence of this series. In Figure 7 we show a log-log plot of |a 2n | −1/(2n) versus 1/n. For n > 10 the data are well fitted by the form |a 2n | −1/(2n) ∝ n −ν with ν ≃ 1.08, shown as a dashed line on the figure.
This result indicates rather convincingly that the large-n limit is zero, and thus that the radius of convergence vanishes. We therefore believe that our expansion is an asymptotic expansion.

Discussion
In this paper we have determined numerically the Lyapounov exponent γ(E) and the density of states ρ(E) of the one-dimensional non-Hermitian Schrödinger equation (1.2) when the phases θ n and χ n are uniformly distributed. We have also developped a method of expanding these quantities in powers of the energy E and obtained γ(E) up to order E 16 . Our procedure can be extended to obtain, in principle, an arbitrarily high order.
The expansion in powers of E is based on the expansion of the integral equation (2.14) satisfied by the invariant distribution P (r, E). Since this expansion is in principle only valid for E ≪ r, there is a possibility that the perturbation series for γ(E) suffers from the fact that when we integrate over r, we use an expression valid for r ≫ E even in the neighborhood of r = 0. One cannot exclude for example that when r and E become comparable (and small), P (r, E) becomes a scaling function of r/E. This could invalidate our expoansion of γ(E). However when we compared our expansion with the results of the simulation we found an excellent agreement for the first two terms and so it is possible that our expansion, a priori valid for r ≫ E, remains valid down to r = 0.
Our numerical results show that the density of states vanishes outside a circle of radius |E| = 2, is non-zero inside this circle |E| < 2, and has a non-trivial shape with a pronounced maximum at an energy E = 2/ One could try to extend our approach to other situations, in particular to cases where the phases θ n and χ n have a non-uniform distribution [3]. In that case, the phase and the modulus of the ratio ψ n+1 /ψ n are in general correlated and one should look for the invariant distribution P (r, E) of the complex variables r n and E. Still, the fact that γ(E = 0) = 0 would remain true, and one could try to expand the invariant distribution in powers of E and use the Thouless formula to calculate ρ(E).
A Nature of the singularities of P (r, E) In this Appendix we discuss the occurence of singularities in the invariant measure P (r, E) at the fixed point r + of the map T + (r) = E + 1 r , when one iterates the random map (2.12) with uniformly distributed ϕ n between 0 and 2π.
We believe that the apparent singularities in P (r, E) seen in Figure 2 are due to events where several successive r n are close to the fixed point r + , and that the nature of the singularity can be understood by analysing only the neighbourhood of r + .
In physical terms, there is a competition between the deterministic part of the map T + (r) which tries to concentrate the distribution on its attractive fixed point r + , and the noise due to ϕ n which tends to broaden the distribution. To analyse the neighbourhood of r + , we consider the linearised problem where a ∈]0, 1[ is a fixed slope, and t n is a random positive number. The variable s n represents the difference r + − r n , when this difference is small. By expanding Eq. (A.1) to second order in ϕ n , it is seen that for Eqs. (A.1) and (A.2) to be equivalent in the neighbourhood of r + , one should choose The essential feature of the distribution of t n is that all the t n are positive and that the distribution Q(t) of t n diverges as Q(t) ∼ t −1/2 as t → 0. If we iterate (A.2) numerically for E = 1.7, that is for a = 0.213851 · · ·, we recover at s = 0 a singularity (see Figure 8.a) which ressembles the one seen on Figure 2.c. Actually, the two distributions are approximately mirror images near their respective fixed points, but this is due to the fact that s n ≃ r + − r n rather than s n ≃ r n − r + .
Trying to analyse the nature of this singularity, we see from Eq. (A.2) that s n can be written as s n = t n−1 − at n−2 + a 2 t n−3 − a 3 t n−4 + a 4 t n−5 − . . . where y n and y ′ n are two positive independent and equally distributed random variables of the form y n = τ n + a 2 τ n−1 + a 4 τ n−2 + a 6 τ n−3 + · · · , y ′ n = τ ′ n + a 2 τ ′ n−1 + a 4 τ ′ n−2 + a 6 τ ′ n−3 + · · · . (A.6) All the τ n and the τ ′ n are independent and distributed according to the same distribution Q(t) as for the t n .
Let us first consider the case where s n is still given by Eq. (A.5) but where the distribution π(y) of the positive variables y is arbitrary. Using Eq. (A.5) one can show that the distribution P (s) of s is given by P (s) = If we assume that π(y) and its first n+1 2 derivatives vanish at the origin, and as y → ∞, it is easy to see by integration by parts that Eqs. (A.9) and (A.10) coincide, so that P (s) has at least n derivatives at 0.
We are now going to argue that since y n is given by Eq. Because Q(t) ∼ t −1/2 for small t, one has for large β For large β, one has e −βy ≃ e −βa 2 y and one finds that to leading order log e −βy ≃ 1 8 log 2 β log a .
(A. 13) We see that e −βy vanishes faster (as log a < 0) than any negative power of β as β → ∞.
As a consequence π(y) and all its derivatives vanish at y = 0.
Therefore we can conclude that all the derivatives of the distribution P (s) are continuous at s = 0. This is confirmed by a magnification of the small-s region, as seen in Figure 8.b, where the rounding of P (s) at s = 0 becomes visible.
Coming back to the stationary distribution P (r, E), we observe the same rounding of the apparent singularity at r + (see Figure 3).

B Justification of Eq. (4.7)
Imagine that s is a random positive variable distributed according to some given distribution Q(s) and let r be given by r = 1 s + Ee iϕ , (B.1)