Jean-Marc Luck
email: jean-marc.luck@ipht.fr

An investigation of PT -symmetry breaking in tight-binding chains

We consider non-Hermitian PT -symmetric tight-binding chains where gain/loss optical potentials of equal magnitudes ±iγ are arbitrarily distributed over all sites. The main focus is on the threshold γc beyond which PT -symmetry is broken. This threshold generically falls off as a power of the chain length, whose exponent depends on the configuration of optical potentials, ranging between 1 (for balanced periodic chains) and 2 (for unbalanced periodic chains, where each half of the chain experiences a non-zero mean potential). For random sequences of optical potentials with zero average and finite variance, the threshold is itself a random variable, whose mean value decays with exponent 3/2 and whose fluctuations have a universal distribution. The chains yielding the most robust PT -symmetric phase, i.e., the highest threshold at fixed chain length, are obtained by exact enumeration up to 48 sites. This optimal threshold exhibits an irregular dependence on the chain length, presumably decaying asymptotically with exponent 1, up to logarithmic corrections.

Introduction

Non-Hermitian Hamiltonians involving complex optical potentials have been used in various guises for a long time, in order to model either inelastic scattering or absorption, in analogy with the use of complex refraction indices in optics. Non-Hermitian operators have gained much attention lately, both in classical and quantum physics (see [1] for a recent comprehensive review). Non-Hermitian Hamiltonians generically have complex spectra, with a resulting non-unitary dynamics. It was however emphasized by Bender and collaborators [2,3], generalizing an observation made by Bessis and Zinn-Justin on the Lee-Yang Hamiltonian H = p 2 + ix 3 (see [4]), that a whole class of non-Hermitian Hamiltonians may have real spectra. These are the PT -symmetric Hamiltonians, which are invariant under the combined action of spatial parity (P) and time reversal (T ) (see [5] for a full technical exposition, and [6,7] for more accessible overviews). The generic scenario is the following. Let H be a non-Hermitian PT -symmetric Hamiltonian, obeying [H, PT ] = 0. If H is close to being Hermitian, in the sense that H † -H is small enough, H and PT have a complete set of common eigenvectors, all energy eigenvalues of H are real, and the dynamics is unitary. The system is in a PT -symmetric phase, and PT -symmetry is said to be unbroken. If H † -H gets larger, it may occur that some of the eigenvectors of H are not eigenvectors of PT any more. This is possible since the PT operator is not linear, but anti-linear (see [5] for details). The corresponding energy eigenvalues are complex, so that the norm of a generic state vector blows up exponentially fast in time. PT -symmetry is said to be broken. The number of pairs of complex conjugate energy eigenvalues usually increases with H † -H . The threshold for PT -symmetry breaking corresponds to the very first occurrence of such a pair of complex eigenvalues.

Among the vast variety of PT -symmetric classical and quantum systems that have been investigated so far in more than 3,000 published papers, tight-binding models are definitely among the easiest ones to grasp, at least on the theoretical side. Such models describe e.g. optical arrays consisting of coupled units, some of them having gains and losses. Most studies concern the one-dimensional geometry of tight-binding chains [8,9,10,11,12,13,14,15]. Various cases of finite chains consisting of N = 2M sites, either pristine or disordered, hosting one or several PT -symmetric pairs of non-Hermitian impurities have been considered. Scattering states and transport properties of PT -symmetric chains coupled to an environment through leads have also been considered [16,17,18,19]. Most of these studies have dealt with the situation where all gains and losses have the same magnitude γ, i.e., the optical potentials of all sites with gain (respectively, with loss) read +iγ (respectively, -iγ). The threshold for PT -symmetry breaking is then a single number γ c , whose dependence on various model parameters has been explored in some detail in various circumstances. The presence of a single PT -symmetric pair of sites with gain and loss already yields a rich phenomenology [9,10,11,15]. The situation where each site bears a random gain/loss optical potential ±iγ has also been considered [14]. The dependence of the threshold on the chain length has been found to strongly differ from one situation to another: it reaches a finite limit in the case of a single pair of PT -symmetric impurities at the endpoints of a long pristine chain [9,11], whereas it is exponentially small in M in the insulating regime induced by Anderson localization, i.e., when the on-site energies have strong enough random real parts [8,12,14]. For PT -symmetric gain/loss disorder alone, numerical simulations have revealed that the threshold has a smooth distribution whose mean value falls off slowly as a function of the chain length [14].

The aim of this paper is to put some of the findings recalled just above in a broader perspective, by providing a systematic study of the threshold of PT -symmetric tightbinding chains of length N = 2M where gain/loss optical potentials are distributed over all sites in a periodic, random, or any other fashion. For definiteness, we focus our attention onto a minimal setting, where all gains and losses have the same magnitude γ, whereas the chains are otherwise pristine, in the sense that there are no other sources of inhomogeneity, besides the optical potentials iγε n , with ε n = ±1. This model easily lends itself to an exact enumeration of all possible sequences ε n of reduced optical potentials for a given chain size. In this setting, the threshold for PT -symmetry breaking is a single number γ c , depending only on the arrangement of the reduced optical potentials ε n . Right at this threshold, the spectrum usually exhibits two symmetry-related exceptional points, i.e., two (usually twofold) degenerate eigenvalues at opposite energies (±E c , with E c > 0), or, albeit more rarely, a single exceptional point, i.e., a single (again usually twofold) degenerate eigenvalue at zero energy. More complex patterns of exceptional points may however occur at γ c . We shall meet two examples of such situations in the course of this work. In Section 7 we provide an example with M = 6, where four twofold degeneracies take place simultaneously at non-zero energies given by (7.3). In Appendix B it is shown that for M = 4 with generic optical potentials there are four isolated points in the γ 1 -γ 2 plane for which a fourfold degeneracy occurs at threshold at zero energy.

Throughout this work, the main emphasis is on the scaling behavior of the threshold γ c as a function of the chain length for various classes of configurations of optical potentials. Section 2 contains the theoretical framework of the model, including the transfer-matrix approach and the discrete Riccati formalism. General results are presented in Section 3, including explicit results for all chains up to M = 3 (Section 3.1) and an overview of the behavior of the threshold for various classes of sequences (Section 3.2). More specific results for these classes of sequences are derived in the subsequent sections. This detailed study includes diblock chains and other unbalanced periodic chains (Section 4), alternating chains and other balanced periodic chains (Section 5), random chains (Section 6), and the chains yielding the most robust PT -symmetric phases, i.e., the largest threshold at fixed chain length (Section 7). Section 8 is devoted to a brief discussion. Three appendices contain more technical matters.

Theoretical framework

The model

The focus of this work is on the time-independent PT -symmetric tight-binding equation

ψ n-1 + ψ n+1 + iγε n ψ n = Eψ n (2.1)
on a finite open chain of N = 2M sites, with Dirichlet boundary conditions

ψ 0 = ψ 2M+1 = 0. (2.2)
Each site bears a purely imaginary gain/loss optical potential of magnitude γ, whereas ε n = ±1 is an arbitrary sequence of signs, with ε n = +1 (abbreviated to ε n = +) corresponding to gain, and ε n = -1 (abbreviated to ε n = -) to loss. We always assume that the model is PT -symmetric. This translates to the constraint

ε 2M+1-n = -ε n , (2.3)
imposing that the sequence of signs is of the form

ε 1 , . . . , ε 2M = ε 1 , ε 2 , . . . , ε M , -ε M , . . . , -ε 2 , -ε 1 .
(2.4)

In the following, (ε 1 , . . . , ε 2M ) will be referred to as the sequence defining the chain, and (ε 1 , . . . , ε M ) as the corresponding half-sequence. For a given PT -symmetric sequence, the spectrum of (2.1) consists of 2M energy eigenvalues E a , which can either be real or occur in complex conjugate pairs (see Section 2.3 for more details). This spectrum is invariant under P and T separately, so that one can choose γ ≥ 0 without loss of generality. We are mainly interested in the threshold γ c associated with the occurrence of the first pair of complex conjugate eigenvalues, and especially in the dependence of γ c on the length M and on the nature (periodic, random, etc.) of the half-sequence. It is worth recalling here that only PT -symmetric sequences give rise to a non-zero γ c . This can be shown by means of first-order perturbation theory (see [14], and Appendix A for a detailed proof).

Transfer-matrix approach

The transfer-matrix approach is a very useful tool in the theory of one-dimensional disordered systems [START_REF] Crisanti | Products of Random Matrices in Statistical Physics Springer Series in Solid-State Sciences[END_REF][START_REF] Luck | Systèmes désordonnés unidimensionnels[END_REF][START_REF] Pendry | [END_REF]23,[START_REF] Comtet | Impurity models and products of random matrices Stochastic Processes and Random Matrices[END_REF]. In the present situation, the tight-binding equation (2.1) can be recast as

ψ n+1 ψ n = T n ψ n ψ n-1 , ψ n-1 ψ n = T n ψ n ψ n+1 , (2.5)
where the transfer matrix T n associated with site n reads

T n = E -iγε n -1 1 0 . (2.6)
Iterating the recursion (2.5) along the half-sequence and along its PT -symmetric partner, using the Dirichlet boundary conditions (2.2) and the PT -symmetry constraint (2.3), we respectively obtain ‡

ψ M+1 ψ M = U M ψ 1 0 , ψ M ψ M+1 = U ⋆ M ψ 2M 0 , (2.7)
with

U M = T M . . . T 1 = a M b M c M d M , U ⋆ M = T M+1 . . . T 2M = a ⋆ M b ⋆ M c ⋆ M d ⋆ M . (2.8)
We have therefore

ψ M = c M ψ 1 = a ⋆ M ψ 2M , ψ M+1 = a M ψ 1 = c ⋆ M ψ 2M .
(2.9)

The quantization condition which determines the spectrum of (2.1) therefore reads

a M a ⋆ M = c M c ⋆ M , i.e., |a M | = |c M |.
(2.10)

We have c n = a n-1 , whereas the entries a n obey the linear recursion

a n = (E -iγε n )a n-1 -a n-2 , (2.11)
with a 0 = 1 and a -1 = 0. ‡ Throughout this paper a star denotes complex conjugation.

Discrete Riccati formalism

Discrete Riccati variables are another classic tool of the theory of one-dimensional disordered systems [START_REF] Crisanti | Products of Random Matrices in Statistical Physics Springer Series in Solid-State Sciences[END_REF][START_REF] Luck | Systèmes désordonnés unidimensionnels[END_REF][START_REF] Pendry | [END_REF]23,[START_REF] Comtet | Impurity models and products of random matrices Stochastic Processes and Random Matrices[END_REF]. In the present situation, they read

R n = a n a n-1 .
(2.12)

These ratios obey the non-linear recursion

R n = E -iγε n - 1 R n-1 , (2.13)
known as a discrete Riccati equation, with initial value R 0 = ∞. In terms of these variables, the quantization condition (2.10) reads

R M R ⋆ M = 1, i.e., |R M | = 1. (2.14)
The above condition yields after reduction a polynomial equation of the form P (E, γ) = 0, where the characteristic polynomial P has degree 2M in E and in γ, has real coefficients, and is even in each of the variables E and γ (see examples in Section 3.1). These properties have several consequences.

For a generic value of γ, the spectrum of (2.1) consists of 2M distinct energy eigenvalues, which are the zeros of P . If E is an eigenvalue, so are E ⋆ , -E and -E ⋆ . These four numbers degenerate into two either if E is real (and so E ⋆ = E) or if E is imaginary (and so E ⋆ = -E). They degenerate into one if E = 0.

At the threshold γ c , there are usually two (usually twofold) degenerate eigenvalues at opposite real energies (±E c , with E c > 0), or, albeit more rarely, a single (usually twofold) degenerate eigenvalue at zero energy. As announced in the Introduction, more complex degeneracy patterns may take place at γ c .

The recursion (2.13) implies

|R n | + 1 |R n-1 | ≥ E 2 + γ 2 (2.15)
whenever E is real. If the r.h.s. of the above inequality is larger than 2, it can be shown by recursion that |R n | > 1 for all n, so that the quantization condition (2.14) cannot be satisfied. In other words, all real eigenvalues of (2.1) obey

|E| ≤ 4 -γ 2 , (2.16)
independently of the length M . The real spectrum is therefore reduced with respect to that of a free tight-binding particle (|E| ≤ 2). It shrinks down to zero as γ → 2. In particular, all sequences have γ c < 2.

General results

Shortest chains

As already said, the spectrum of (2.1) is invariant under P and T separately. As a consequence, besides choosing γ ≥ 0, one can fix ε 1 = +, and so ε 2M = -, without loss of generality. There are therefore 2 M-1 essentially different PT -symmetric sign sequences on a chain of length N = 2M , each of them being encoded in a half-sequence of length M starting with +. We begin our investigation with a detailed study of the shortest chains. There are 1, 2 and 4 different half-sequences of lengths M = 1, 2 and 3, i.e., chain lengths N = 2M = 2, 4 and 6. These sequences are successively investigated below. Figure 1 shows their real eigenvalues in the E-γ plane for γ ≥ 0. Blue symbols show the thresholds γ c and the corresponding degenerate eigenvalues. These plots exhibit a variety of topologically different patterns, hinting at the complexity of the problem.

-

2 -1 0 1 2 E 0 0.5 1 1.5 2 γ + + -- (2)
-2 -1 0 1 2 E 0 0.5 1 1.5 2 γ + -+ - (3)
Figure 1. Real eigenvalues of the seven shortest chains in the E-γ plane (see text). The full PT -symmetric sequence is given in each case. All plots are at the same scale in order to facilitate the comparison between them. Blue symbols: thresholds γc and corresponding degenerate eigenvalues.

(1) The half-sequence (+) is the simplest of all. It yields the full PT -symmetric sequence (+-). Its characteristic polynomial is

P = E 2 + γ 2 -1. (3.1)
The energy eigenvalues E = ± 1 -γ 2 are real for γ < 1 and imaginary for γ > 1.

At the threshold

γ c = 1 (3.2)
there is a twofold degenerate eigenvalue at zero energy.

(2) The half-sequence (++) has characteristic polynomial

P = E 4 + (2γ 2 -3)E 2 + γ 4 + γ 2 + 1. (3.3)
At the threshold

γ c = √ 5 4 ≈ 0.559 016 (3.4)
there are degenerate energy eigenvalues at ±E c , with

E c = √ 19 4 ≈ 1.089 724. ( 3.5)
(3) The half-sequence (+-) has characteristic polynomial

P = E 4 + (2γ 2 -3)E 2 + γ 4 -3γ 2 + 1. (3.6)
At the threshold

γ c = √ 5 -1 2 ≈ 0.618 033 (3.7)
there is a twofold degenerate eigenvalue at zero energy.

(4) The half-sequence (+ + +) has characteristic polynomial

P = E 6 + (3γ 2 -5)E 4 + (3γ 4 -2γ 2 + 6)E 2 + γ 6 + 3γ 4 + 2γ 2 -1.(3.8)
At the threshold γ c ≈ 0.322 142, (3.9)

obeying 2048γ 6 c + 560γ 4 c + 392γ 2 c -49 = 0, (3.10)
there are degenerate energy eigenvalues at ±E c , with

E c ≈ 1.506 895, (3.11)
obeying

2048E 6 c -7216E 4 c + 6632E 2 c -1831 = 0. ( 3.12)
(5) The half-sequence (+ + -) has characteristic polynomial

P = E 6 + (3γ 2 -5)E 4 + (3γ 4 -6γ 2 + 6)E 2 + γ 6 -γ 4 -2γ 2 -1.(3.13)
At the threshold

γ c = 1 2 (3.14)
there are degenerate energy eigenvalues at ±E c , with

E c = √ 3 2 ≈ 0.866 025. (3.15) (6)
The half-sequence (+ -+) has characteristic polynomial

P = E 6 + (3γ 2 -5)E 4 + (3γ 4 -10γ 2 + 6)E 2 + γ 6 -5γ 4 + 6γ 2 -1.(3.16)
At the threshold

γ c = 2 sin π 14 ≈ 0.445 041 (3.17)
there is a twofold degenerate eigenvalue at zero energy. (7) The half-sequence (+ --) has characteristic polynomial

P = E 6 + (3γ 2 -5)E 4 + (3γ 4 -6γ 2 + 6)E 2 + γ 6 -γ 4 + 2γ 2 -1.(3.18)
At the threshold

γ c ≈ 0.534 036, (3.19)
obeying

16γ 3 c + 16γ 2 c -7 = 0, (3.20)
there are degenerate energy eigenvalues at ±E c , with

E c ≈ 1.418 436, (3.21) obeying 256E 6 c -768E 4 c + 576E 2 c -135 = 0. (3.22)

Panoramic overview

We now turn to an overview of the behavior of the threshold γ c for PT -symmetry breaking for various kinds of configurations (periodic, random, etc.) of the optical potentials. We have run a computer routine which enumerates all the 2 M-1 halfsequences of length M starting with ε 1 = +, up to a maximal length M = 24, i.e., N = 48, and determines for each sequence accurate numerical values of the threshold and of the corresponding degenerate energy eigenvalues. The most notable emerging features are illustrated in Figure 2, and will be studied in detail in the forthcoming sections. In the regime of large chains, thresholds fall off to zero, generically as a power law of the form

γ c ∼ 1 M β . (3.23)
Different classes of sequences yield different decay exponents β. The smallest of all γ c , decaying with exponent β = 2, is reached for the diblock chains having a constant half-sequence (see black dataset and Section 4.1). The same scaling holds for all unbalanced periodic sequences with a non-zero mean value (see Section 4.2). Alternating half-sequences, as well as all balanced periodic ones with zero mean value, have thresholds falling off with exponent β = 1 (see green dataset and Section 5). The mean value of γ c over all half-sequences of given length scales with exponent β = 3/2. The same scaling holds for typical random sequences, i.e., for almost all half-sequences (see red dataset and Section 6). Finally, the largest threshold at fixed M has an irregular dependence on M and presumably also scales asymptotically as 1/M , up to logarithmic corrections (see blue dataset and Section 7).

Diblock and other unbalanced periodic chains

Diblock chains

This section is devoted to diblock chains, which are the most unbalanced of all PTsymmetric chains. On a diblock chain of of length N = 2M , the M sites with gain are to the left, and the M sites with loss are to the right. In other terms, the halfsequence is constant, i.e., ε 1 = . . . = ε M = +, and so ε M+1 = . . . = ε 2M = -. It is expected that the PT -symmetric phase is the most fragile on these diblock chains. This is corroborated by the observation made in Section 3.2 that these chains yield the smallest threshold at fixed chain length (see black dataset in Figure 2). The latter observation will now, in turn, be made quantitative by analytical means. Instead of using a straightforward approach borrowed from elementary Quantum Mechanics, i.e., looking for eigenfunctions of (2.1) as linear combinations of plane waves with complex wavenumbers in each block, we prefer to have recourse to the Riccati formalism introduced in Section 2.3, since the latter approach will prove more powerful in more complex situations.

The recursion (2.13) obeyed by the Riccati variables does not depend explicitly on n and reads

R n = 2 cosh µ - 1 R n-1 (n = 1, . . . , M ), (4.1)
where

µ = η + iq (4.2) is a complex variable such that 2 cosh µ = E -iγ, (4.3)
i.e., E = 2 cos q cosh η, γ = -2 sin q sinh η.

(4.4)

The fixed points of (4.1) are R = e ±µ . Setting

Y n = R n -e -µ R n -e µ , (4.5)
the recursion (4.1) boils down to the multiplication by a constant, i.e.,

Y n = e 2µ Y n-1 . (4.6)
The initial value R 0 = ∞ yields Y 0 = 1, hence Y n = e 2nµ , and finally

R n = sinh(n + 1)µ sinh nµ . (4.7)
The quantization condition (2.14) therefore reads

sinh(M + 1)µ sinh(M + 1)µ ⋆ = sinh M µ sinh M µ ⋆ , (4.8)
i.e., sin q sin(2M + 1)q + sinh η sinh(2M + 1)η = 0. (4.9)

The above equations are exact for any finite M . We are mainly interested in the scaling behavior of γ c for long chains. It turns out that the most unstable energy eigenvalues are close to the band edges (E → ±2). Let us focus our attention onto the upper band edge (E → 2, i.e., q → 0) and introduce the scaling variables

x = (2M + 1)q, g = (2M + 1) 2 4 γ. (4.10)
We have then (2M + 1)η ≈ (-2g/x), and so the quantization condition (4.9) becomes

sin x ≈ - 2g x 2 sinh 2g x , (4.11)
up to corrections of relative order 1/(2M + 1) 2 .

In the absence of optical potentials (g = 0), the r.h.s. of (4.11) vanishes. We thus obtain

x ≈ aπ (a = 1, 2, . . .), (4.12)
in perfect agreement with (A.2). As the rescaled magnitude g of gains and losses increases, the rescaled spectrum is more and more distorted, especially so for the very first values of a, as shown in Figure 3. The blue line (g = 0) intersects the sine function at integer values of x/π (see (4.12)). The green curve (g = 2) yields a distorted albeit still entirely real spectrum. The red one, corresponding to

g c ≈ 4.475 308, ( 4.13)
is tangent to the sine function at

x c ≈ 5.331 598 (4.14) (blue symbol), where the lowest pair of eigenvalues (a = 1 and 2) merges before becoming complex. Other mergings between a = 3 and 4, a = 5 and 6, and so on, take place at higher well-defined critical values of g.

0 1 2 3 4 5
x/π In brief, for long diblock chains of length N = 2M , the threshold for PTsymmetry breaking scales as

γ c ≈ 4g c (2M + 1) 2 ∼ g c M 2 . (4.15)
This provides a proof of the scaling law announced in Section 3.2, with exponent β = 2, together with a prediction for the amplitude g c . The first two pairs of degenerate eigenvalues sit very near the band edges, at ±E c , with

2 -E c ≈ q 2 -η 2 ≈ 4σ c (2M + 1) 2 ∼ σ c M 2 , (4.16)
with

σ c = x 2 c 4 - g 2 c
x 2 c ≈ 6.401 903. (4.17)

Both estimates (4.15) and (4.16) hold up to corrections of relative order 1/(2M + 1) 2 .

Figure 4 provides a quantitative check of the prediction (4.15), including the order of magnitude of the corrections. The combination (2M + 1) 2 γ c /4 is plotted against 1/(2M + 1) 2 , up to M = 50. The blue line has the theoretical intercept g c given in (4.13).

Unbalanced periodic chains

We now turn to unbalanced periodic chains. A PT -symmetric chain is said to be periodic if the half-sequence ε 1 , . . . , ε M is periodic, i.e., ε n = ε n+L for some period L ≥ 1 and all n = 1, . . . , M -L. The motif (ε 1 , . . . ε L ) is called the unit cell of the periodic sequence. In this section we make the assumption that the periodic sequence is unbalanced, i.e., its mean value ε is non-zero. § In the present case, we have where L + and L -are the numbers of plus and minus signs (i.e., of sites with gain and loss) in the unit cell, with

ε = 1 L L n=1 ε n = L + -L - L , ( 4
L + + L -= L.
For the time being, let us follow a heuristic line of thought. If ε is positive, on a spatial scale much larger than the period L, the left half of the chain has a roughly homogeneous gain, whereas the right half has a roughly homogeneous loss, and vice versa if ε is negative. Any unbalanced periodic sequence therefore amounts to an effective diblock sequence, with γ being replaced by the effective value γ eff = γε. This approach predicts that the threshold of a long unbalanced periodic chain is asymptotically given by the product of the expression (4.15) for diblock chains and of an amplification factor

f = 1 |ε| = L |L + -L -| , (4.19)
depending only on the contents of the unit cell of the sequence. This reads

γ c ≈ 4f g c (2M + 1) 2 ∼ f g c M 2 . (4.20)
Consistently, diblock chains themselves can be viewed as unbalanced periodic ones with L = L + = 1, and so f = 1. Allowing larger and larger periods, the amplification factor f may become either arbitrarily close to unity or arbitrarily large: if L ≥ 3 and

L ± = 1, we have f = L/(L -2); if L is odd and L ± = (L -1)/2, we have f = L.
The result (4.20) turns out to be asymptotically correct, in spite of its heuristic derivation. This will be shown in a broader setting in Section 6.1, by means of a continuum Riccati approach.

The non-trivial unbalanced chain with shortest period has L = 3 and unit cell (+ + -), hence f = 3. Figure 5 provides a quantitative check of the prediction (4.20) on this example. Here again, the combination (2M + 1) 2 γ c /4 is plotted against 1/(2M +1) 2 , up to M = 50. Data for finite chains converge to the predicted limit value 3g c ≈ 13.425 (horizontal blue line), oscillating with the period L = 3 of the underlying sequence of gains and losses. Oscillations are damped linearly on the chosen scale. , where γc is the threshold of unbalanced chains with unit cell (+ + -). Horizontal blue line: predicted limit value 3gc ≈ 13.425.

Alternating and other balanced periodic chains

Alternating chains

In this section we consider alternating chains, where sites with gains and losses alternate along the chain, so that ε n = (-1) n . If the chain length N = 2M is even, this configuration of optical potentials obeys PT -symmetry. Alternating chains are periodic with period L = 2 and balanced, in the sense that ε = 0, so that the analysis of Section 4.2 does not apply.

Our goal is to derive the exact threshold γ c for any finite chain length N = 2M (see (5.11)). The analysis again relies on the Riccati formalism. The recursion (2.13) takes two different forms according to the parity of n:

R 2p+1 = E + iγ - 1 R 2p , R 2p+2 = E -iγ - 1 R 2p+1 . (5.1)
These two equations can be combined, yielding

R 2p+2 = (E 2 + γ 2 -1)R 2p -E + iγ (E + iγ)R 2p -1 . (5.2)
Introducing the parametrization

E + iγ = 2 cos q e iθ , (5.3)
i.e., E = 2 cos q cos θ, γ = 2 cos q sin θ,

and so

E 2 + γ 2 = 4 cos 2 q, (5.5)
the fixed points of (5.2) are R = e -i(θ±q) . Setting

Y n = R n -e -i(θ+q) R n -e -i(θ-q) , (5.6)
the recursion (5.2) boils down to

Y 2p+2 = e 4iq Y 2p .
(5.7)

The initial value R 0 = ∞ yields Y 0 = 1, hence Y 2p = e 4ipq , and finally R 2p = e -iθ sin(2p + 1)q sin 2pq , R 2p+1 = e iθ sin(2p + 2)q sin(2p + 1)q .

(5.8)

The quantization condition (2.14) reads sin(2M + 1)q = 0, (5.9)

irrespective of θ and of the parity of M . This equation is considerably simpler than its counterpart (4.9) for diblock chains. It yields the explicit result

q = kπ 2M + 1 (k = 1, . . . , M ).
(5.10)

For each quantized value of q, labelled by the integer k, a pair of eigenvalues merges at zero energy as θ → π/2, whereas γ reaches the critical value γ k = 2 cos q k . The pattern of real eigenvalues in the E-γ plane therefore has the form of a rainbow of M nested arches (see panels (1), ( 3) and ( 6) of Figure 1). The first merging event corresponds to the smallest critical value, i.e., k = M . We thus obtain the exact expression

γ c = 2 sin π 2(2M + 1) (5.11)
for the threshold of alternating PT -symmetric chains of any finite length. We recover known values of γ c for the first values of M , namely

γ c = 2 sin π/6 = 1 for M = 1 (see (3.2)), γ c = 2 sin π/10 = ( √ 5 -1)/2 for M = 2 (see (3. 7
)), and γ c = 2 sin π/14 for M = 3 (see (3.17)). For long alternating chains, the threshold falls off as

γ c ≈ π 2M .
(5.12)

This law of decay with exponent β = 1 was announced in Section 3.2.

One example with period four

We now turn to an example of balanced periodic chains where pairs of sites with gains and losses form an alternating pattern. These chains have period L = 4 and unit cell (+ + --). The analysis again relies on the Riccati formalism. The recursion (2.13) reads

R 4p+1 = E -iγ - 1 R 4p , R 4p+2 = E -iγ - 1 R 4p+1 , R 4p+3 = E + iγ - 1 R 4p+2 , R 4p+4 = E + iγ - 1 R 4p+3 . (5.13)
These four equations can be combined into a recursion expressing R 4p+4 as a function of R 4p , whose expression is more intricate than (5.2). Anticipating that the threshold corresponds to a merging of eigenvalues around E c = √ 2 (see the discussion below (5.34)), we set E = √ 2 + a, and simplify the recursion in the scaling regime where both a and γ are small. We thus obtain

R 4p+4 = (1 -2 √ 2 a)R 4p + 2iγ + 4a 1 + 2 √ 2 a + (2iγ -4a)R 4p . (5.14)
Introducing parameters w, φ and θ such that

a = 1 2 w cos φ, γ = w sin φ, cos θ = cos φ √ 2 , (5.15)
with 0 < φ < π and π/4 < θ < 3π/4, the fixed points of (5.14) are R = e i(φ±θ) . Setting

Y 4p = R 4p -e i(φ-θ) R 4p -e i(φ+θ) , (5.16)
the recursion (5.14) boils down to

Y 4p+4 = e -2iλ Y 4p , (5.17)
with tan λ = 2w sin θ.

(5.18)

The initial value R 0 = ∞ yields Y 0 = 1, hence Y 4p = e -2ipλ , and finally

R 4p = e iφ sin(pλ -θ) sin pλ , R 4p+1 = √ 2 -e -iφ sin pλ sin(pλ -θ) . (5.19)
For M = 4p (and also for M = 4p + 2), the quantization condition (2.14) reads

sin(pλ -θ) = ± sin pλ, (5.20)
yielding

λ = θ + kπ 2p , (5.21)
with k integer. The threshold is expected to correspond to the smallest value of λ, i.e., k = 0. It can be shown that the rescaled combinations

G = M γ, x = M a = M (E - √ 2), (5.22)
asymptotically only depend on θ, according to

G = θ √ -cos 2θ sin θ , x = θ cot θ √ 2 .
(5.23)

The threshold γ c corresponds to the point in the x-G plane where the curve given parametrically by (5.23) has a horizontal tangent, as in Figure 1, i.e., to the maximum of the function G(θ). As θ is varied between π/4 and 3π/4, G reaches its maximum G 1 for θ 1 such that θ 1 = cos 2θ 1 tan θ 1 , i.e., θ 1 ≈ 1.937 183, where we have

G = (θ + α) √ -cos 2θ sin θ , x = (θ + α) cot θ √ 2 .
(5.28)

The threshold is obtained by the maximization procedure applied to (5.23) To sum up, for large chains with even M , the threshold scales as

γ c ≈ G 1 M , (5.30)
and corresponds to two pairs of degenerate eigenvalues at ±E c , with

E c ≈ √ 2 + x 1 M , (5.31)
where G 1 and x 1 are given in (5.24). For large chains with odd M , the threshold scales as

γ c ≈ G 2 M , (5.32)
and corresponds to two pairs of degenerate eigenvalues at ±E c , with

E c ≈ √ 2 + x 2 M , (5.33)
where G 2 and x 2 are given in (5.29).

Figure 6 provides a quantitative check of the above predictions. The product M γ c is plotted against 1/M , up to M = 100. The upper (resp. lower) dataset corresponds to odd (resp. even) values of M . The blue lines have the theoretical intercepts G 1 and G 2 (see (5.24)

Balanced periodic chains

Let us now turn to balanced periodic chains with even period L and ε = 0 (see (4.18)).

Only a few general features can be stated in the present situation. This is in strong contrast with unbalanced chains, investigated in Section 4.2, where the quantitative prediction (4.20) has been derived in full generality.

The Fourier transform of a periodic chain with period L is supported by all integer multiples of the reciprocal period 2π/L. This singles out L -1 resonant energies inside the band of the free tight-binding chain, i.e., E a = 2 cos aπ L (a = 1, . . . , L -1).

(5.34)

For long balanced periodic chains (M ≫ 1), it is expected that the energy E c at which the first merging of eigenvalues takes place approaches one of those resonant energies.

The convention E c ≥ 0 selects a = 1, . . . , L/2, the last case yielding E c = 0. For the alternating chain (see Section 5.1), L = 2 and a = 1 correctly predict E c = 0. For the chain studied in Section 5.2, L = 4 and a = 1 correctly predict E c = √ 2, since a = 2 is ruled out by symmetry. For chains with larger periods L, there is no systematic way of predicting the resonant energy E c . There are actually examples where E c keeps varying periodically with M mod L.

The threshold for PT -symmetry breaking is expected to decay as

γ c ≈ G c M (5.35)
for all balanced periodic chains. For the alternating chain (see Section 5.1), we have G c = π/2 (see (5.12)). For the chain with period four studied in Section 5.2, G c alternates between the two numbers G 1 and G 2 , according to the parity of M (see (5.24), (5.29)). For chains with larger periods, the amplitude G c is generically expected to vary periodically with M mod L. Many questions remain open, including whether the value π/2 obtained for alternating chains is a lower bound for G c and whether G c may become arbitrarily large.

Random chains

This section is devoted to random PT -symmetric chains, where each symbol of the half-sequence (ε 1 , . . . , ε M ) is taken to be ε n = ± with equal probabilities. Equivalently, the 2 M different half-sequences are considered as equally probable. The threshold γ c for PT -symmetry breaking is now a random variable, as it depends on the chain under scrutiny. This situation has been studied by numerical means in [14].

Continuum Riccati formalism

The forthcoming analysis again relies on the Riccati formalism introduced in Section 2.3. We are mainly interested in the behavior of the threshold γ c for long random chains. In this regime, γ c is typically small, whereas the first mergings concern pairs of eigenvalues close to the band edges (E → ±2).

Let us focus for definiteness our attention onto the upper band edge (E → 2). In the vicinity of the latter, setting E = 2 -q 2 , in the presence of small optical potentials iγ n = iγε n , it is advantageous to consider the small variable

Z n = R n -1, (6.1)
in terms of which the recursion (2.13) reads

Z n = -q 2 -iγ n + Z n-1 1 + Z n-1 = -q 2 -iγ n + Z n-1 -Z 2 n-1 + • • • (6.2)
This expansion suggests to use a continuum formalism, considering Z(n) as a continuous function of the real variable n, obeying the Riccati differential equation dZ dn = -q 2 -iγ(n) -Z 2 , (6.3) with initial condition

Z(0) = ∞. (6.4)
In spite of this, the function Z(n) is small all over the relevant range. In particular, if the optical potential is constant, setting µ 2 = -q 2 -iγ, the solution to (6.3) reads

Z = µ cotanh nµ. (6.5)
To leading order as µ → 0, for fixed values of the product nµ, this is in agreement with (4.7), which can be recast as

R n = cosh µ + sinh µ cotanh nµ = 1 + µ cotanh nµ + 1 2 µ 2 + • • • (6.
6) Before we pursue with random chains, let us come back for a while to unbalanced periodic chains, investigated in Section 4.2. If the sequence ε n is periodic with period L and has a non-zero mean value ε, the continuous function γ(n) has the same period and mean value γ = γε. In the situation of interest, namely sequences of length M ≫ L, the function Z(n) varies typically on a scale M . It is therefore slowly varying at the scale of one period, so that it is legitimate to replace γ(n) by its mean value γ. This puts the heuristic line of thought of Section 4.2 on a firmer basis.

Predictions for random chains

For long random PT -symmetric chains, the function γ(n) becomes in the continuum limit a Gaussian white noise with covariance

γ(m)γ(n) = γ 2 δ(m -n).
(6.7)

The continuum Riccati equation ( 6.3) for the function Z(n) has to be integrated from n = 0, with initial condition (6.4), to n = M , the length of the half-sequence, where the quantization condition (2.14) translates to

Re Z(M ) = 0. (6.8)

The Riccati differential equation driven by a real Gaussian white noise, i.e., dZ dx = -q 2 -V (x) -Z 2 , (6.9) with V (x)V (y) = w 2 δ(x -y), (6.10) has a long history, dating back to pioneering studies in the 1960s of Anderson localization in a one-dimensional white-noise potential [START_REF] Frisch | [END_REF]26]. The Lyapunov exponent Γ = Z(∞), (6.11) representing the inverse localization length of the model on an infinite line, obeys a scaling law of the form Γ ≈ 1 2 w 2 1/3 F (x), x = 1 2 w 2 -2/3 q 2 , (6.12)

in the low-energy scaling regime where q and w are simultaneously small. The scaling function F (x) is known in terms of Airy functions:

F (x) = Ai(x)Ai ′ (x) + Bi(x)Bi ′ (x) Ai(x) 2 + Bi(x) 2 . ( 6.13)
The result (6.12) is commonly referred to as anomalous band-edge scaling. It has been rederived several times, both in continuum and in lattice models [27,28,29]. Several works have been devoted to the fluctuations of the process Z(x) and of closely related ones around their stationary mean value [30,31,32]. The mathematics involved there is quite intricate. Fortunately enough, predictions of interest concerning the statistics of the threshold γ c of long random chains can be derived from the sole consideration of scale invariance. The continuous Riccati equations (6.3), (6.9) driven by noises obeying (6.7), (6.10) are invariant under the scaling transformation

(n, x) → 1 λ (n, x), (Z, q) → λ(Z, q), (γ, w) → λ 3/2 (γ, w). (6.14)
For the problem of Anderson localization, the scale invariance of the Riccati equation (6.9) dictates the scaling form (6.12) of the Lyapunov exponent. Of course, the determination of the exact form (6.13) of the scaling function F (x) requires a more advanced analysis.

In the present context, the full scale invariance of the Riccati equation ( 6.3), including its boundary conditions (6.4) at n = 0 and (6.8) at n = M , implies that the mean value of the threshold of long random chains scales as

γ c ≈ A M 3/2 , (6.15)
with a decay exponent β = 3/2, as announced in Section 3.2, and that the reduced random variable

ξ = γ c γ c (6.16)
has a well-defined limiting distribution with density f (ξ), such that ξ = 1. It should be clear from the above reasoning that the latter distribution is universal, in the sense that it would hold for any symmetric, i.e., even distribution of gain/loss optical potentials γ n with finite variance. The determination of the amplitude A and of the distribution f (ξ) by analytical means remains a difficult open problem. We therefore have to rely on numerical analysis. Figure 7 shows the product M 3/2 γ c , plotted against 1/M . Red circles show the outcome of a systematic enumeration of all half-chains, available up to M = 24. These data were already presented in Section 3.2 (see red dataset in Figure 2). Blue squares show the outcome of numerical simulations with 10 6 samples for M = 30, 40 and 50. The black curve is a quadratic fit to all data points with intercept 5.660, confirming thus the scaling law (6.15) to a high degree of accuracy, and yielding the estimate A ≈ 5.7.

(6.17 Figure 8 shows the reduced combination of moments

Q = γ 2 c (γ c ) 2 , (6.18)
plotted against 1/M , with the same conventions as in Figure 7. The black curve is a quadratic fit to all data points with intercept 1.351. The observed convergence to a well-defined limit provides a first quantitative piece of information on the limiting distribution f (ξ), namely an estimate of its second moment,

ξ 2 ≈ 1.
Q Figure 8
. Plot of the quantity Q defined in (6.18) against 1/M . Same conventions as in Figure 7.

Figure 9 shows a histogram plot of the distribution of the reduced random variable ξ defined in (6.16). Data are based on numerical simulations with 10 6 samples for M = 40 and 50 (see legend). Every second point of each dataset is plotted alternatively. A good collapse is observed, implying that finite-size corrections are small for the chosen chain lengths, and so that the plotted data give an accurate representation of the limiting distribution f (ξ). The black curve shows a very accurate four-parameter fit of the form

f (ξ) ∼ exp - a ξ 2 -b -cξ -dξ 2 . (6.20)
This functional form is meant to capture approximately the asymptotic behavior of the exact distribution at small and large ξ. For a fixed length M , the smallest threshold is that of diblock sequences (see Section 4.1). We have γ c ∼ 1/M 2 , and so ξ ∼ 1/ √ M , with probability 1/2 M , suggesting f (ξ) ∼ exp(-a/ξ 2 ). The largest threshold, to be investigated in Section 7, scales as γ c ∼ (ln M )/M , and so ξ ∼ √ M ln M , again with probability 1/2 M , suggesting f (ξ) ∼ exp(-dξ 2 ), up to logarithmic corrections. It should however be clear that the exact distribution f (ξ) is expected to be infinitely more complex than (6.20).

Most robust PT -symmetric phase

In this section we consider the chains yielding the most robust PT -symmetric phase, i.e., the highest threshold γ c at fixed length M . We describe in some detail the outcome of the systematic investigation described in Section 3.2, yielding exact data up to M = 24, i.e., N = 48. The value of the threshold, as well as many other characteristics of these most robust chains, exhibit an irregular dependence on the length M of the half-sequence. Let us begin with the value of the highest threshold itself. It has been shown in Section 5 that the thresholds of balanced periodic sequences fall off as γ c ≈ G c /M (see (5.35)), where the amplitude G c has a complicated dependence on the length and unit cell of the sequence. Roughly speaking, the highest threshold γ c at fixed length M must be larger than all those thresholds for all periods L ≤ M . As it turns out, γ c is larger than 1/M by an irregularly varying but altogether very slowly growing factor. This is demonstrated in Figure 10, showing a plot of M times the highest threshold γ c against ln M . These data were already presented in Section 3.2 (see blue dataset in Figure 2). The red line with slope 1.82 shows an acceptable least-square fit to the data, suggesting the behavior Table 1 provides a list of all half-sequences (ε 1 , . . . , ε M ) yielding the highest threshold shown in Figure 10, for all lengths up to M = 24, with the convention ε 1 = + (see Section 3.1). These sequences can be viewed as consisting of alternating domains (+ -+ -. . .) separated by domain walls, in the words of magnetism, denoted by vertical bars. The number of domain walls and their positions optimize the value of the threshold. We have not succeeded in identifying any pattern in those optimal positions.

γ c ≈ C ln M M , C ≈ 1.8. ( 7
In order to illustrate the complexity of the dependence of γ c on the positions of the domain walls, we consider in Appendix C alternating chains having a single pair of PT -related domain walls between sites K and K + 1 and between sites 2M -K and 2M -K + 1. Already in this situation, the threshold is found to have an intricate asymptotic dependence on the ratio α = K/M , involving four different regimes. This threshold is always enhanced by the presence of the domain walls, with respect to that of pristine alternating chains (see (5.12)). The enhancement factor takes its maximal value 3 for α = 1/3. Some degeneracies can be seen in Table 1. For M = 6, three different halfsequences yield the maximal γ c = 1/2, with mergings occurring at different critical energies (see below). For M = 16 and M = 19, two different half-sequences yield the same maximal threshold and mergings at the same critical energies. The threshold γ c generically corresponds to the appearance of two twofold degenerate eigenvalues at opposite energies ±E c , where E c > 0 has itself an irregular dependence on M (data

M ε 1 , . . . , ε M 1 + 2 +- 3 + -|- 4 +| + -+ 5 + -| -+- 6    + -+| + -|- + -| -+ -+ +| + -+ -+ 7 + -| -+ -+- 8 + -| -+ -+ -+ 9 + -| -+ -+ -+|+ 10 + -+| + -+ -+ -+ 11 + -+| + -+ -+ -+ - 12 +| + -+ -| -+ -+ -+- 13 + -+ -+ -+| + -+ -| -+ 14 + -+ -+ -| -+ -+| + -| -+ 15 + -+| + -| -+ -+ -+| + -+ - 16 +| + -+ -+ -+ -| -| -+| + -+ - +| + -+ -+ -+ -| -+| + | + -+ - 17 + -+ -+ -+| + -| -+| + -+ -| -+ 18 + -+| + -| -+ -+| + -+ -+ -| -+- 19 +| + -+ -+ -+ -+ -| -| -+| + -+ -+ +| + -+ -+ -+ -+ -| -+| + | + -+ -+ 20 +| + -+ -+ -| -+| + -+ -| -+| + -+ | + - 21 + -| -+| + -| -+| + -+ -| -+ -+| + -+ -+ 22 + -+| + -+ -+ -+ -+ -| -+| + -+ -| -+|+ 23 + -+ -+ -+| + -| -+| + -| -+ -+| + -+ -| -+ 24 +| + -+ -+ -| -+ -+ -+ -+ -+| + -| -+ -+|+ Table 1.
Half-sequences yielding the highest thresholds γc at fixed length M , up to M = 24. Vertical bars: domain walls between alternating domains. not shown). For M = 2, 7, 8, 12 and 21, there is a single twofold degenerate eigenvalue at zero energy.

The case M = 6 is special in several regards. Three distinct half-sequences, listed in Table 1, yield the highest threshold γ c = 1/2. This simple value allows for an analytical study of the corresponding energy spectra. For the first half-sequence (+ -+ + --), the first mergings occur at ±E c , with

E c = √ 7 2 ≈ 1.322 875. (7.2)
For the second half-sequence (+ --+ -+), four mergings occur simultaneously at ±E c1 and ±E c2 , with

E c1 = 5 -2 √ 5 2 ≈ 0.363 271, E c2 = 5 + 2 √ 5 2 ≈ 1.538 841. (7.3)
For the third half-sequence (+ + -+ -+), which is the only unbalanced one, the first mergings occur at ±E c , with

E c = √ 3 2 ≈ 0.866 025. (7.4)

Discussion

This paper has been devoted to non-Hermitian PT -symmetric tight-binding chains of N = 2M sites, where gain/loss optical potentials of equal magnitudes, of the form iγε n with ε n = ±1, are distributed over all sites in a periodic, random, or any other fashion.

We have provided a systematic study of the threshold γ c for PT -symmetry breaking, with the main emphasis being on the dependence of γ c on the chain length for various classes of configurations of optical potentials. Our main findings are summarized in the panoramic overview of the asymptotic behavior of γ c for various classes of chains given in Section 3.2 (see Figure 2). The threshold generically falls off as a power of the chain length, i.e.,

γ c ∼ 1 M β . (8.1)
The decay exponent β depends on the configuration of optical potentials. The largest exponent β = 2 is reached for diblock chains, and more generally for all unbalanced periodic chains, i.e., periodic configurations of optical potentials where each half of the chain has on average either gain or loss. The smallest exponent β = 1 is attained for alternating chains, and more generally for all balanced periodic chains, i.e., periodic configurations of optical potentials whose spatial average over one period vanishes. For random sequences of optical potentials with zero average and finite variance, the threshold γ c is itself a random variable, whose mean value decays with exponent β = 3/2, whereas the fluctuations around this mean value obey a universal distribution f (ξ).

These results have been corroborated and made quantitative by detailed analytical studies performed in Sections 4 to 6. This analysis relies on concepts and techniques from the physics of one-dimensional disordered systems, including the transfer-matrix approach and the discrete Riccati formalism.

The above results can be put in a broader context as follows. Consider the sum

S = ε 1 + • • • + ε M (8.2)
of all reduced optical potentials over the left half of the chain, and assume that this sum grows as

S ∼ M ω , (8.3)
with a wandering exponent in the range 0 ≤ ω ≤ 1. The three classes of configurations of optical potentials recalled above respectively yield β = 2, 1 and 3/2, whereas they clearly obey (8.3) with ω = 1, 0 and 1/2. This observation suggests the hyperscaling relation

β = ω + 1. (8.4)
This relation provides a quantitative form of the rule of thumb according to which the more the sites with gain and loss are well mixed along the chain, the more robust the corresponding PT -symmetric phase. A large class of deterministic aperiodic sequences have non-trivial wandering exponents ω, besides the three classical values recalled above [33,34,[START_REF] Queffélec | Substitution Dynamical Systems[END_REF]36]. The relation (8.4) certainly holds in those cases as well, thus predicting a whole zoo of non-trivial asymptotic decay exponents β.

The chains yielding the most robust PT -symmetric phase, i.e., the highest threshold γ c at fixed chain length, have been investigated by means of a systematic enumeration up to M = 24, i.e., N = 48 (see Section 7). Both the optimal γ c itself and the sequences responsible for it exhibit an irregular dependence on M . The highest threshold presumably falls off asymptotically as (ln M )/M . These logarithmic corrections somehow represent an exception to the above rule of thumb. The chains yielding the optimal thresholds indeed contain several domain walls, i.e., pairs of consecutive sites with gain or loss. They are therefore slightly less well mixed than pristine alternating ones, whereas the corresponding thresholds are larger.

To close, it would certainly be interesting to extend at least some aspects of the present study to more complex situations, including e.g. a periodic, quasiperiodic or random modulation of on-site energies and/or hopping rates, or quasi-one-dimensional geometries such as ladders, exhibiting either internal band gaps or flat bands.

The forthcoming analysis, including notations, closely follows Section 5.1, which was devoted to pristine alternating chains. For definiteness we restrict ourselves to the case where M and K are even. The expressions (5.8) of the Riccati variables still hold for n = 1, . . . , K. We have therefore in particular R K = e -iθ sin(K + 1)q sin Kq .

(C.2)

For n = K + 1, . . . , M , the recursion obeyed by the Riccati variables R n is obtained by changing γ and θ into their opposites in (5.1) and (5.2). In particular, the variables

Y n = R n -e i(θ-q)
R n -e i(θ+q) (C.3) still obey the recursion (5.7). We have therefore

Y M = e 2iq(M-K) Y K . (C.4)
Using the three above equations, we obtain

R M = e iθ c M -c M+2 + i tan θ (c M + c M+2 -2c M-2K ) c M-1 -c M+1 + i tan θ (c M-1 + c M+1 -2c M-2K-1 )
, (C.5)

with the shorthand notation c n = cos nq. Finally, the quantization condition (2.14) simplifies to 4 sin 2 θ sin Kq sin(K + 1)q sin(2M -2K)q + sin 2 q sin(2M + 1)q = 0. (C.6)

The above equation is exact for any even integers M and K.

For long chains, where both M and K are large, the threshold for PT -symmetry breaking scales as

γ c ≈ G c (α) M , (C.7)
where the numerator has a rich dependence on the ratio

α = K M , (C.8)
as announced in Section 7. The amplitude G c (α) is always larger than that of pristine alternating chains (see (5.12)), which is recovered for α → 0 and α → 1. The effect of the domain walls indeed disappears in both limits, where they either attain the endpoints of the chain or annihilate each other in the middle. The threshold amplification is maximal for α → (1/3) -, where we have G c ((1/3) -)/G c (0) = 3.

Figure C1 shows a plot of G c (α) against α. Numerical results for the thresholds for M = 60 and all possible positions of the domain walls are compared to the exact expression of G c (α), involving four different regimes, to be successively investigated below.

• Regime (1) holds for 0 < α < 1/3. The amplitude increases from G c (0) = π/2 ≈ 1.570 796 to G c ((1/3) -) = 3π/2 ≈ 4.712 388.

• Regime (2) holds for 1/3 < α < α 2 . At α = 1/3 the amplitude drops discontinuously to G c ((1/3) + ) = 3π √ 3/4 ≈ 4.081 048. We have α 2 ≈ 0.583 417 and G c (α 2 ) ≈ 2.697 312.

• Regime (3) holds for α 2 < α < α 3 . The amplitude is continuous at α 2 and at α 3 , with α 3 ≈ 0.776 365 and G c (α 3 ) ≈ 2.841 880. • Regime (4) holds for α 3 < α < 1. The amplitude decreases from G c (α 3 ) to G c (1) = π/2 ≈ 1.570 796.

Analytical expressions of the amplitude G c (α) in each regime can be derived by zooming onto the exact quantization condition (C.6) in appropriate ranges.

In Regimes (1) and ( 4), the first merging takes place at zero energy. This corresponds to θ = π/2, whereas q = π/2 -η, and so γ c ≈ 2η, with η ∼ 1/M . The quantization condition (C.6) asymptotically becomes cos(1 -2α)G c = 0.

(C.9)

The amplitude G c (α) is the smallest positive solution to the above equation. We thus obtain the explicit expressions

G c (α) = π 2(1 -2α)
in Regime (1), (C.10)

G c (α) = π 2(2α -1)
in Regime (4). (C.11)

In Regime (2), the first mergings take place near zero energy. This corresponds to q = π/2 -η with η ∼ 1/M and θ variable, and so E c ≈ 2η cos θ and γ c ≈ 2η sin θ. Introducing the product x = 2M η, the quantization condition (C.6) asymptotically becomes 2 sin 2 θ sin αx sin(1 -α)x + cos x = 0.

(C.12)

The threshold is obtained by the maximization procedure applied to (5.23), yielding G c (α) = max In Regime (3), the first mergings take place near the band edges. This corresponds to both q and θ scaling as 1/M . Introducing the product x = M q, we obtain the asymptotic prediction G c (α) = max where the maximum is taken over x in the range π/2 < x < π at fixed α. The amplitude takes its minimal value G c (2/3) ≈ 2.493 125 for α = 2/3.

Figure 2 .

 2 Figure2. Log-log plot of various thresholds γc as a function of the length M of the half-sequence, ranging from 2 to 24. From bottom to top: Black: smallest γc at fixed M , reached for diblock chains. Red: mean value of γc over all half-sequences of length M . Green: γc for alternating halfsequences. Blue: largest γc at fixed M . Dashed lines have the theoretical slopes β = 2, 3/2 and 1, to be derived in the forthcoming sections.

Figure 3 .

 3 Figure 3. Graphical illustration of the threshold mechanism for the rescaled quantization condition (4.11) near the upper edge of the spectrum of long diblock chains (see text).

4 Figure 5 .

 45 Figure 5. Red symbols: combination (2M + 1) 2 γc/4 plotted against 1/(2M + 1) 2 , where γc is the threshold of unbalanced chains with unit cell (+ + -). Horizontal blue line: predicted limit value 3gc ≈ 13.425.

G 1 ≈ 1 For M = 4p+ 1 (

 111 and also for M = 4p+ 3), the quantization condition (2.14) reads sin 2 (pλ -θ) + sin 2 pλ -4 cos θ sin pλ sin(pλ -θ) α < π. The quantities introduced in (5.22) now read

Figure 6 .

 6 Figure 6. Red symbols: product M γc for chains with unit cell (+ + --), plotted against 1/M up to M = 100. The upper (resp. lower) dataset corresponds to odd (resp. even) values of M . The blue lines have the theoretical intercepts G1 and G2 (see (5.24), (5.29)).

Figure 7 .

 7 Figure 7. Plot of the product M 3/2 γc against 1/M . Red circles: systematic enumeration of all half-chains up to M = 24. Blue squares: numerical simulations with 10 6 samples. Black curve: quadratic fit to all data points.

Figure 9 .

 9 Figure 9. Histogram plot of limiting distribution f (ξ) of reduced random variable ξ defined in (6.16), based on numerical simulations with 10 6 samples for M = 40 and 50 (see legend). Black curve: fit of the form (6.20).

Figure 10 .

 10 Figure 10. Plot of the product M γc against ln M , where γc is the highest threshold at fixed length M , up to M = 24. Red line with slope 1.82: least-square fit.

Figure C1 .

 C1 Figure C1. Plot of the amplitude Gc(α) against the reduced position α of the left domain wall. Full curves: analytical asymptotic predictions (C.10), (C.13), (C.17), (C.11) in the four regimes. Symbols: data for M = 60 and all possible positions of the domain walls.

  is taken over x in the range π/2 < x < 3π/2 at fixed α. The amplitude takes its minimal value G c (1/2) ≈ 2.529 882 for α = 1/2, i.e., when the left domain wall is in the middle of the half-sequence. It has a square-root endpoint singularity at α = 1/3, below which Regime (1) takes over. Setting

Acknowledgments

It is a pleasure to thank Philippe Di Francesco for fruitful exchanges.

Appendix A. First-order perturbation theory

In this appendix we derive within the framework of first-order perturbation theory the spectrum of the tight-binding equation

with Dirichlet boundary conditions, for an arbitrary sequence of purely imaginary optical potentials γ 1 , . . . , γ 2M . Our goal is to show that all energy levels remain real within this framework if and only if the sequence γ n is PT -symmetric.

In the absence of optical potentials, the non-degenerate energy eigenvalues and normalized eigenstates read

To first order in the optical potentials, the energy eigenvalues read E a + iΓ a , with

all the sums in the parentheses in the last line of (A.4) vanish, and therefore all the firstorder imaginary parts Γ a vanish, suggesting the existence of a non-zero threshold γ c for PT -symmetry breaking.

The proof of the reciprocal proceeds as follows. The first-order imaginary parts (A.4) obey the mirror symmetry Γ 2M+1-a = Γ a , so that it is sufficient to consider a = 1, . . . , M . Requesting Γ a = 0 for a = 1, . . . , M yields a system of M linear equations for the M unknown quantities (γ n + γ 2M+1-n ) for n = 1, . . . , M . The corresponding determinant,

with

is non-zero. It can indeed be evaluated by using

). This is useful since the entries of the matrix A 2 have simple expressions [37,38]. Expanding the sines into complex exponentials and performing the sums, we indeed obtain after some algebra

The determinant of a matrix of size M ×M with entries λ+µδ ab reads µ M-1 (M λ+µ).

We thus obtain

The sign of D M is not predicted by this approach. It depends periodically on M mod 4.

To sum up, Γ a = 0 for all a implies γ n +γ 2M+1-n = 0 for all n, i.e., PT -symmetry. This completes the proof.

Appendix B. Chain of length four with arbitrary gains and losses

In this appendix we consider the tight-binding equation

on a PT -symmetric chain of length N = 2M = 4, with arbitrary imaginary optical potentials, obeying

The first two optical potentials γ 1 and γ 2 are kept as free parameters. The quantization condition (2.14) is expressed by the characteristic polynomial B1.

1). (B.3) This polynomial has a rich bifurcation diagram shown in Figure

There is a twofold degenerate eigenvalue at zero energy on the red curves, i.e., for

There are two twofold degenerate eigenvalues at ±E c = 0 (either real or imaginary) on the green curves, i.e., for

The red and green curves have four contact points (blue symbols), at and their opposites, where the characteristic polynomial has a fourfold degenerate eigenvalue at zero energy. These four points are cocyclic, obeying γ 2 1 + γ 2 2 = 3. The above curves define various regions in the γ 1 -γ 2 plane. The structure of the energy spectrum is given in each region, with R, I and C denoting real, imaginary and generic complex eigenvalues. For instance, 2R+2I means two opposite non-zero real eigenvalues and two opposite non-zero purely imaginary ones.

PT -symmetry is unbroken whenever all eigenvalues are real, i.e., in the central region marked 4R. This quadrangular region is bounded by two red arcs and two green ones meeting at the four blue contact points. The predictions (3.4) and (3.7) for the thresholds of binary chains with half-sequences (++) and (+-) can be respectively recovered by intersecting one of the green arcs with the bisector line of equation γ 2 = γ 1 and of the red ones with the bisector line of equation γ 2 = -γ 1 .

Appendix C. Alternating chains with a single pair of domain walls

This appendix is devoted to almost pristine alternating chains of length N = 2M having a single pair of domain walls at PT -related positions between sites K and K + 1 and between sites 2M -K and 2M -K + 1. The corresponding half-sequence reads


Olivier Moranne
email: olivier.moranne@chu-nimes.fr

Jean-Paul Cristol

With Single Plasma Exchange, a Better Understanding of the Potential Clinical Effects of Albumin Replacement Is Required

Keywords: Albumin, Single plasma exchange, Physiology

Therapeutic plasmapheresis (TP) can be either nonselective (single plasma exchange [SPE]), semi-selective (double filtration plasmapheresis) or selective (immunoadsorption) [START_REF] Sanchez | The selective therapeutic apheresis procedures: the Selective Therapeutic Apheresis Procedures[END_REF]. SPE requires replacement fluid infusion with human serum albumin (HSA) to prevent hypotension but could be suggested for treatment infusion, just as fresh frozen plasma is used in thrombotic thrombocytopenic purpura [START_REF] Mcleod | Plasma and plasma derivatives in therapeutic plasmapheresis: plasma for therapeutic apheresis[END_REF]. Indeed, each SPE session requires infusion of 60 g to 120 g of albumin according to the exchanged plasma volume and replacement fluid protocol. Moreover, the exchanged plasma volume depends on weight, hematocrit, and clinical indication. Then HSA volume and concentration could be isotonic and isovolemic or depending to baseline patient's blood albumin concentration and volemia. As HSA has several functions, we need to consider the potential benefits and adverse effects of albumin exchange when comparing different TP procedures. It is therefore crucial to understand the metabolism of albumin, its binding capacities, antioxidant power, and posttranslational modifications.

HSA, the most abundant protein in plasma, is produced by hepatocytes at a rate of 10-15 g per day and about 5% per hour is exchanged between the vascular and extracellular space. Thus, the vascular albumin pool is estimated at 118 g (40%) and the extracellular pool at 242g. Albumin has one of the longest protein half-lives (13-18 days). The tertiary structure is characterized by three homologous domains (DI, DII, DIII) connected via long flexible loops. Two subdomains, known as Sudlow I (DIIA) and 2 (DIIIA), serve as binding sites for lowmolecular-weight ligands [START_REF] Bern | The role of albumin receptors in regulation of albumin homeostasis: implications for drug delivery[END_REF]. Moreover, albumin is the major pool of thiol in plasma with 17 disulfide bridges and one-free cysteine in position 34.

Albumin is the major determinant in plasma's oncotic pressure, mainly due to its molecular mass but also its negative charge which attracts sodium [START_REF] Fanali | Human serum albumin: from bench to bedside[END_REF]. Albumin is also a primary transporter of several physiologic or pathologic substances (fatty acids, hormones, metals [iron, copper], heme, bilirubin, hippuric acid and protein-bound uremic toxins like indoles and phenols, amyloid beta protein, thiols) but also medication (antibiotics, antiepileptics, anticoagulants) [START_REF] Fanali | Human serum albumin: from bench to bedside[END_REF][START_REF] Ghuman | Structural basis of the drug-binding specificity of human serum albumin[END_REF]. Albumin is not simply a sponge: it may also an antioxidant activity [START_REF] Halliwell | Albumin: an important extracellular antioxidant?[END_REF][START_REF] Quinlan | Albumin influences total plasma antioxidant capacity favorably in patients with acute lung injury[END_REF][START_REF] Bourdon | Differential effects of cysteine and methionine residues in the antioxidant activity of human serum albumin[END_REF], anti-infectious effects [START_REF] Schneider | A pilot study on continuous infusion of 4% albumin in critically ill patients: impact on nosocomial infection via a reduction mechanism for oxidized substrates[END_REF] and interactions with the endothelium affecting vascular integrity and coagulation [START_REF] Van De Wouw | Albumin is an interface between blood plasma and cell membrane, and not just a sponge[END_REF]. We know the potential clinical implications of HAS infusion with SPE in several conditions like thyrotoxicosis [START_REF] Tańska | Therapeutic plasma exchange with albumin as a valuable method of preparing thyrotoxic patients for a life-saving thyroidectomy[END_REF], liver failure [START_REF] Kapoor | Therapeutic plasma exchange in liver failure[END_REF], Alzheimer [START_REF] Boada | Neuropsychological, neuropsychiatric, and quality-of-life assessments in Alzheimer's disease patients treated with plasma exchange with albumin replacement from the randomized AMBAR study[END_REF], medication intoxication [START_REF] Shi | MARS: optimistic therapy method in fulminant hepatic failure secondary to cytotoxic mushroom poisoning-a case report[END_REF], and another ongoing study [15]. Furthermore, posttranslational albumin modifications, which could be detected using liquid chromatography coupled to tandem mass spectrometry, are observed in diabetes, end-stage renal disease, or age-related disease. These modifications impair HSA physiological properties and may reduce HSA ability to bind ligands, decrease drug bioavailability, change colloid oncotic pressure, enzymatic properties, or impair antioxidant capacity mainly due to cysteine and thiol oxidation [16]. Moreover, modified albumin gains new properties including neutrophil stimulation and enhances reactive oxygen production leading to inflammation, endothelial dysfunction and, ultimately, cardiovascular disease [16]. Furthermore, HSA infusion may cause adverse events with the risk of contamination by infectious agents and secondary deficiency of substances fixed to albumin like medication, vitamins, hormones, and metals [17].

Thus, HSA exchange may be interesting but there is little information available on the potential modifications induced by therapeutic albumin preparations including the redox state of different pharmaceutical products [18] and substances fixed on HSA after ethanol treatment [19] except for binding metals which are controlled by national instance. Therefore, when SPE is associated with HSA infusion, we must consider the potential effects of albumin exchange and gain a better understanding by evaluating its clinical effects, including adverse effects. HSA is currently a limited, expensive resource which, in future, could be replaced by synthetic/recombinant albumin [START_REF] Fanali | Human serum albumin: from bench to bedside[END_REF]20].

Acknowledgments

We thank Teresa Sawyers, Medical Writer at the BESPIM, Nîmes University Hospital, France, for her expertise in reviewing and editing this paper.

Funding Sources

The authors wish to thank Nîmes University Hospital, France, for its structural, human, and financial support through the award obtained by our team during the internal call for tenders "Thématiques émergentes" (emerging themes).

Conflict of Interest Statement

The authors have no conflicts of interest to declare.

Author Contributions

O.M. developed the idea and O.M. and J.-P.C. wrote the manuscript.