Numerical calculation of the complex berry phase in non-Hermitian systems

We numerically investigate topological phases of periodic lattice systems in tight-binding description under the influence of dissipation. The effects of dissipation are effectively described by $\mathcal{PT}$-symmetric potentials. In this framework we develop a general numerical gauge smoothing procedure to calculate complex Berry phases from the biorthogonal basis of the system's non-Hermitian Hamiltonian. Further, we apply this method to a one-dimensional $\mathcal{PT}$-symmetric lattice system and verify our numerical results by an analytical calculation.


Introduction
Due to their robustness against local defects or disorder topologically protected states as Majorana fermions [1][2][3] are of high value for physical applications such as quantum computation [4].However, no physical system is completely isolated and dissipation can have an important influence on the states [5].Majorana fermions can even be created with the help of dissipative effects [6,7].
Of special importance in this context is the case of balanced gain and loss as described by PT -symmetric complex potentials [8], which has attracted large interest in quantum mechanics [9][10][11][12][13].The stationary Schrödinger equation was solved for scattering solutions [14] and bound states [15], and also questions concerning other symmetries [16,17] as well as the meaning of non-Hermitian Hamiltonians have been discussed [18,19].Their influence on many-particle systems has been studied mainly in the context of Bose-Einstein condensates [20][21][22][23][24][25] but also on lattice systems [26][27][28][29][30][31].In the latter systems it was shown that PT -symmetric complex potentials may eliminate the topologically protected states existing in the same system under complete isolation [26-30, 32, 33].However, in some PT -symmetric potential landscapes they can survive [28][29][30][31][32]34], which has been confirmed in impressive experiments [33,35] In most theoretical studies the topologically protected states have been identified via their property to close the energy gap of the band structure or their localisation at edges or interfaces of the systems [26,[28][29][30]32].The identification and calculation of topological invariants such as the Zak phase [36] known from Hermitian systems leads to new challenges in the case of non-Hermitian operators [27,37,38].This is especially true if the eigenstates are only available numerically.Indeed, in Refs.[27,37,38] all calculations have been done for eigenstates which are analytically available.However, for arbitrary PT -symmetric complex potentials analytical access to the eigenstates is not available and a reliable numerical procedure is required.For the calculation of the invariants an integration of phases over a loop in parameter space is typically necessary.
For example, in the case of the Zak phase, which is applied to one-dimensional systems, this integral runs over the first Brillouin zone.The integrand contains the eigenstates of the Hamiltonian and their first derivatives.In a numerical calculation it is evaluated at discrete points in momentum space and each state possesses an arbitrary global phase spoiling the phase relation.
This is the point our study sets in.In this article we introduce a robust method of calculating the complex Berry phase numerically.It is based on a normalisation of the left-and right-hand eigenvectors with the biorthogonal inner product [39], which reduces to the c product [40] in the case of complex symmetric Hamiltonians.To obtain an unambiguous complex Berry phase we introduce a numerical gauge smoothing procedure.It consists of two parts.First we have to remove the influence of the arbitrary and unconnected global phases of the eigenstates, which is unavoidably attached to them for each point in parameter space.This is achieved by relating the eigenvectors of consecutive steps in parameter space, and then normalising them.With this approach the eigenvectors are not yet single-valued, i.e. the vectors at the starting and end point of the loop possess different phases.These points have to be identified and will be refered as the basepoint of the loop later on.The phase difference between the different left and right states at the basepoint has to be corrected by ensuring the eigenvectors to be identical at the basepoint.
The article is organised as follows.First, we introduce the complex Berry phase in section 2. In section PREPRINT 3 we establish the algorithm of the gauge smoothing procedure for non-Hermitian (and Hermitian) Hamiltonians.An example is presented in section 4, where we apply the previously developed method to a non-Hermitian extension of the Su-Schrieffer-Heeger (SSH) model [41] to calculate its complex Zak phase.

Complex Berry phase
Topological phases of closed one-dimensional periodic lattice systems are characterised by the Zak phase [36], which is the Berry phase [42] picked up by the eigenstate when it is transported once along the Brillouin zone.In the presence of an antiunitary symmetry these phases are quantised [43] and can be related to the winding number of a vector n(k) determining the Bloch Hamiltonian where σ denotes the vector of Pauli matrices and k is the wave number parametrising the Brillouin zone, which acts as parameter space.In this case the Zak phase characterises the system's topological phase.This concept can be generalised to dissipative systems effectively described by a PT -symmetric non-Hermitian Hamiltonian H(α).The complex Berry phase γ n of a biorthogonal pair of eigenvectors χ n | and |φ n of H(α), follows from the lowest order of the adiabatic approximation of the time evolution of a state in parameter space [44].Here C is a loop in parameter space and α = (α 1 , ..., α i , ...) are its coordinates.We consider PT -symmetric non-Hermitian Hamiltonians of the form where H h denotes its Hermitian and H nh represents its PT -symmetric non-Hermitian part.The non-Hermitian part is a complex potential modelling the gain and loss of particles.The complex Berry phase γ n arising from the periodic modulation of states in the parameter space of a PT -symmetric one-dimensional system cannot be related to a real winding number calculated from eq. ( 1).Hence, the calculation of the complex Berry phase requires the determination of gauge-smoothed eigenvector pairs along the loop C in parameter space allowing for the evaluation of eq. ( 2).
Here, it is important to note that the argumentation of Hatsugai [43] for the quantisation of Berry phases of Hermitian Hamiltonians can be extended on complex Berry phases of non-Hermitian PT -symmetric Hamiltonians in the case of unbroken PT symmetry.One finds the real part of the complex Berry phase to take values 0 or π modulo 2π.Thus, a strict quantisation is still present and the PT symmetry protects the topological phases occurring in such systems.

Numerical gauge smoothing
In this section we present a numerical procedure to determine the left and right eigenvectors χ n | and |φ n in an appropriate smoothed gauge to compute complex Berry phases on a discretised loop C = (α 1 , ..., α j , ..., α M = α 1 ) in parameter space.This is necessary for the evaluation of integrals of the form in eq. ( 2).Typically the left and right eigenvectors have to be calculated independently, and each of them has an arbitrary global phase.The biorthogonal normalisation condition [39], chooses one arbitrary global phase for each α j .This is sufficient if only products or matrix elements of eigenstates belonging to the same point in parameter space are required.However, for numerical derivatives used in eq. ( 2) the remaining global phases of successive steps in α j along the loop C can spoil the complex Zak phase.A fixation of the phase between consecutive steps that does not distort the desired result is required.
Starting point for the gauge smoothing procedure are the left and right handed versions of the timeindependent Schrödinger equation, defining a set of natural left and right basis states χ n (α)| and |φ n (α) .These equations are solved for every point α j of the discretised loop C in parameter space providing the eigenvalues E n (α j ) and the unnormalised states of a biorthogonal basis { χ n (α j )|, |φ n (α j ) } of the Hamiltonian H(α j ) at each point α j .Here the basis states are determined up to the aforementioned arbitrary phases.
To smooth the gauge within the loop in parameter space with basepoint α 1 one chooses an arbitrary global phase.It is most convenient to do this for the basepoint.The corresponding eigenstates are normalised according to the conditions (4a) and (4b).The following two-stage procedure transfers the choice of the global phase at α 1 onto the other basis states along the loop C.
First one modifies the phases of the states χ n (α j )| and |φ n (α j ) iteratively by followed by a normalisation of the states according to eqs.(4a) and (4b).Eqs.(6a) and (6b) relate the vectors of step j to those of step j − 1 by ensuring no. / Numerical Calculation of the Complex Berry Phase which is a valid condition in the continuous limit.The normalisation conditions (4a) and (4b) ensure that the basis states now fulfil χ m (α j )|φ n (α j ) = δ mn (8) for j ∈ {1, ..., M } and for all n and m.
As a result of the first step the arbitrary global phases have been removed.Only one arbitrary phase is left, which has no influence, since it is identical for all right eigenvectors and its complex conjugate for all left eigenvectors.However, the biorthogonal basis following from the procedure so far is not singlevalued in the parameter space.In particular, the vectors at the starting and end point of the loop are not identical.For the calculation of a Berry phase a continuous single-valued phase function is essential [42], and thus has to be established.
To this end in the second step one adjusts the basis states such that they are the same at the starting and the end point of the loop.This can be achieved by compensating the phase difference between the states at the basepoint, χ n (α 1 )| and χ n (α M )|, respectively, as well as |φ n (α 1 ) and |φ n (α M ) .This remains true for single vector components.Therefore we calculate the phase difference of the first non-vanishing component p of the left basis states χ n (α 1 )| and where ϕ n,j = arg χ n (α j )| p is the argument of component p of the left eigenvector at the point α j in parameter space and X denotes the sum of directed crossings of the phase ϕ n,j over the borders of the standard interval [−π, π).Starting with X = 0 we increase X by one for every jump of ϕ n,j from −π to π and subtract 1 for the opposite direction.The states of the biorthogonal basis can then finally be gauge-smoothed by multiplying the states at α j by a phase factor according to for j ∈ {1, ..., M }, where f ∆ϕn (x) is any "smooth" real valued continuous function with z ∈ Z.Its explicit form is not critical since it only has to correct the total phase change over the whole range of the loop.However, a linear progression of the phase correction from step to step turns out to be a good choice.
It should be mentioned that in case of a degeneracy of the eigenvalue at α j the solution of eqs.(5a) and (5b) yields an arbitrary linear combination of eigenvectors of the degenerate eigenspace.To find the correct eigenvectors χ n (α j )| and |φ n (α j ) one can apply a biorthogonal Gram-Schmidt algorithm [45].If α j−1 is a point neighbouring the degeneracy one tries to find a linear combination of the vectors of the left degenerate eigenspace fulfilling χ m (α j )|φ n (α j−1 ) ≈ δ mn (12) and then chooses the right eigenvectors such that Alternatively one can treat the real and imaginary parts of the degenerate eigenvector components as "smooth" functions.Then the eigenvector components at degeneracy points can be predicted by fitting a spline to the vector components at neighbouring points α l .An approximation to the correct eigenvectors of the degenerate eigenspace can be determined by a linear combination of the obtained vectors of the degenerate eigenspace such that they fit best to the prediction.Hermitian Hamiltonians can be treated as a special case, in which the left eigenvector fulfils

Application to a one-dimensional lattice system
In this section we apply the gauge smoothing procedure developed in section 3 to a PT -symmetric onedimensional lattice system to calculate the complex Zak phase where the parameter space is given by the discretised Brillouin zone BZ and k is the wave number.
As an example we consider the SSH model [41] with N lattice sites, tunnelling amplitudes t + and t − , and creation (annihilation) operators of spinless fermions We apply an alternating non-Hermitian PT -symmetric potential of the form where Γ denotes the parameter of gain and loss.The PT -symmetric Hamiltonian describing this model (sketched in Figure 1) is given by To evaluate eq. ( 14) we need to represent this Hamiltonian in the reciprocal space, where the Brillouin zone M. Wagner et al. acts as parameter space.This is done by rewriting the Hamiltonian with creation and annihilation operators in the reciprocal space in the limit N → ∞,

PREPRINT
where the sum runs over discrete values of k in steps of k = 2π/N and the annihilation operator of an electron with wave number k is given by with r n = an and the lattice spacing a.The matrix occurring in eq. ( 18) is the Bloch Hamiltonian H(k) of the system, which can be decomposed into the Pauli matrices, with a coefficient vector n and the Pauli vector σ.
From this form the energy eigenvalues can be obtained explicitly, In the limit Γ → 0 the Hamiltonian from eq. ( 18) reproduces the Hermitian SSH model, which possesses time-reversal, reflection, particle-hole, and a chiral symmetry (represented by σ 3 ).The introduction of a PT -symmetric non-Hermitian on-site potential Γ breaks these symmetries.The non-Hermitian Bloch Hamiltonian is invariant under the combined action of the parity and the time inversion operator.Further the particle-hole symmetry is broken because the sources (sinks) of an electron correspond to sinks (sources) of holes.The non-Hermitian system is therefore symmetric under the action of the combination of the parity and the charge conjugation operator.Further it has no chiral symmetry Λ = a 0 σ 0 + a • σ because a chiral symmetry would fulfil with a coefficient vector a which is independent of the value of k , the 2 × 2 identity matrix σ 0 , and the anti-commutation relations {σ i , σ j } = 2δ ij σ 0 of the Pauli matrices.Therefore one finds a 0 = 0, and thus which cannot be satisfied for a constant vector a because the vector n(k) rotates on a cylindrical surface as k runs through the Brillouin zone.Hence the non-Hermitian Hamiltonian in eq. ( 18) does not possess a chiral symmetry.However, this does not mean there is no quantised real part of the Zak phase since its quantisation is ensured by the argument of Hatsugai [43] in the PT -unbroken parameter regime as mentioned above.At the critical point Γ = 0 the system reproduces the Hermitian SSH model, which possesses the previously mentioned symmetries.For Γ < 0 the particle sinks and sources are interchanged leading to a spatially reflected system with the same general properties as the system with Γ > 0.
From the Bloch Hamiltonian (cf.eq. ( 18)) the complex Zak phase can be calculated following the steps explained in section 3. We choose (c.f.eq.(11a)) which is the most simple function fulfilling the conditions (11b). Figure 2 illustrates the gauge smoothing process using the first non-vanishing component p = 1 of the left basis states as an example.The component χ 1 (α j )| 1 of the unnormalised left eigenvector is shown in Figure 2 (a).One can see two different phase branches as a result of the numerical diagonalisation and a constant modulus.After the gauge smoothing and normalisation according to the first step described in section 3 as shown in Figure 2 (b) the modulus varies with the wave number k and there is only one phase branch left, but the basis is not yet single-valued as there is still a phase difference at the boundaries of the Brillouin zone.In this example the factor X = −1 has to be used as ϕ 1,j jumps from −π to π.After the second step the component χ 1 (α j )| 1 of the final left eigenvector is the same at the Brillouin zone boundaries in Figure 2 (c 4a) and (4b) (here ∆ φ1 = ϕM,1 − ϕ1,1 and X = −1 c.f. eq. ( 9)) the phase is smooth within the Brillouin zone but discontinuous at its boundaries, and the modulus varies with k.(c) After the gauge smoothing process the phase difference ∆ φ1 is compensated and the phase is continuous and smooth in the whole Brillouin zone.
is used to describe the difference between the two tunnelling amplitudes t + and t − , with the mean value of the tunnelling amplitude t and the dimerization strength ∆.
To verify our results we compare them with the analytical ones derived in [37], where q = t + /t − is the ratio of the tunnelling amplitudes, η = Γ/(2t − ) and ν = 4q/((q + 1) 2 − η 2 ) and µ = 4q/(q + 1) 2 .K(ν) and Π(µ, ν) are elliptic integrals of first and third kind, The numerical calculations perfectly reproduce the analytical results.The grey shaded areas in Figure 3 17) in dependence of the control parameter θ with t = 1, ∆ = 1/2 and Γ = 1 (all values in a.u.) compared to the analytical result (real and imaginary part are represented by a solid black line coinciding with the numerical results) following from eq. ( 26).The grey shaded area marks the PT -broken phase of the Bloch Hamiltonian.
mark the values of θ for which the system is in a PTbroken phase, for all other values of θ the system is in a PT -unbroken phase.In the PT -symmetric regime the real part of the complex Zak phase is either 0 or π modulo 2π and can be used to characterise the topological phase of the system.

Conclusion
We developed a numerical method to determine a gauge-smoothed biorthogonal basis of eigenstates of a PT -symmetric non-Hermitian Hamiltonian required for complex Zak phases.It is also applicable to Hermitian systems.In the course of this we removed the arbitrary and unconnected global phases of the biorthogonal eigenstates of the PT -symmetric Hamiltonian at each point in parameter space and made the basis single-valued.This allows for the calculation of the complex Berry phase by explicitly evaluating eq.
(2) even in complicated lattice systems for which no analytical access to the eigenstates is approachable.We demonstrated the action of the gauge smoothing method by applying the developed algorithm to a PTsymmetric extension of the SSH model.An excellent agreement of the numerical and analytical results was found.In future work this provides the basis for the identification of the Zak phase as topological invariant in many-body systems with arbitrary PT -symmetric complex potentials.

Figure 1 .
Figure 1.(Colour online) Sketch of the SSH model with N lattice sites subject to the alternating imaginary potential U .The minus (plus) sign marks a negative (positive) imaginary potential corresponding to particle sinks (sources).
).The two complex Zak phases γ 1 and γ 2 following from the eigenvector pairs of H(k) are shown in Figures 3 (a) and (b), where the control parameter θ

Figure 2 .
Figure 2. (Colour online) First component of the left handed eigenvector χ1(αj)|1 in dependence of the wave number k with t = 1, ∆ = 1/2, Γ = 1 and θ ≈ 0.3π : (a) Before the steps described in eqs.(6a) and (6b) one can identify two different phase branches (blue line) and a constant modulus (filled red dots).(b) After the steps characterised by eqs.(4a) and (4b) (here ∆ φ1 = ϕM,1 − ϕ1,1 and X = −1 c.f. eq.(9)) the phase is smooth within the Brillouin zone but discontinuous at its boundaries, and the modulus varies with k.(c) After the gauge smoothing process the phase difference ∆ φ1 is compensated and the phase is continuous and smooth in the whole Brillouin zone.

Figure 3 .
Figure 3. (Colour online) Numerical results of the real part (filled red dots) and the imaginary part (open blue circles) of the complex Zak phases γ1 (a) and γ2 (b) following from the Hamiltonian given in eq.(17) in dependence of the control parameter θ with t = 1, ∆ = 1/2 and Γ = 1 (all values in a.u.) compared to the analytical result (real and imaginary part are represented by a solid black line coinciding with the numerical results) following from eq. (26).The grey shaded area marks the PT -broken phase of the Bloch Hamiltonian.