{'fraction_non_alphanumeric': 0.033947178714201456, 'fraction_numerical': 0.008781702601743234, 'mean_word_length': 4.327107697678478, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A general algorithm for computing bound states in infinite tight-binding systems

We propose a robust and efficient algorithm for computing bound states of infinite tight-binding systems that are made up of a finite scattering region connected to semi-infinite leads. Our method uses wave matching in close analogy to the approaches used to obtain propagating states and scattering matrices. We show that our algorithm is robust in presence of slowly decaying bound states where a diagonalization of a finite system would fail. It also allows to calculate the bound states that can be present in the middle of a continuous spectrum. We apply our technique to quantum billiards and the following topological materials: Majorana states in 1D superconducting nanowires, edge states in the 2D quantum spin Hall phase, and Fermi arcs in 3D Weyl semimetals.


Introduction
Simulating quantum devices that are connected to infinite leads is a commonly occurring problem in quantum nanoelectronics. In particular, finding the propagating states in the leads and coupling them to the device allows to evaluate transport properties such as the conductance of the device [1] . Numerical methods for solving this scattering problem have a long history [2][3][4] . Modern stable algorithms are available in several software packages, e.g. Refs. 5-7. The focus of this work is identifying bound states: states, that are localized inside the central region and decay exponentially in the leads or any other translationally invariant bulk of the material (see Fig. 1). Although these states do not contribute to transport, they frequently are the central object of study. Examples of bound states relevant to current research include the quantum well states in semiconductor heterostructures, surface states in metals, impurity states (phosphorus donors in silicon, nitrogen-vacancy centers in diamond [8] , Shiba states due to magnetic impurities in superconductors), Andreev states [9] in Josephson junctions, and various kinds of protected edge states in topological materials (e.g. Majorana states in superconducting nanowires [10][11][12] , chiral edge states in quantum spin Hall systems or Fermi arcs in Weyl semi-metals). Computing the bound states can also be desirable from a mathematical or technical perspective: the calculations of Feynman diagrams due to electronelectron interactions requires the full basis of the system (see for instance Ref. 13).
Closed-form solutions of the bound state problem exist for several regimes, such as the strong coupling (or short junction) limit of superconducting junctions [14,15] , or the semiclassical approach to impurity bound states [16] .
To solve that problem, one could also use a brute-force method, by cutting out a finite part of the infinite system and diagonalizing the corresponding sub-Hamiltonian. A more evolved numerical method using the boundary element method has been developed in Ref 17, but its drawback is that the scattering region must consist of several homogeneous region and only compute the bound states at negative energies, while our algorithm does not bear these constraints. If the state decays fast enough in the leads and a sufficiently large portion of them has been included, this results in a precise determination of the bound states. This  approach is however not always satisfactory due to its significant computational overhead when the decay length of the bound state diverges. In addition, this brute-force method by itself does not allow to distinguish the bound states from the continuum spectrum when an exact bound state coexists with the continuum spectrum [18] . We present an algorithm for calculating the bound state spectrum directly for the infinite system. Compared to the brute force method, our approach has the following advantages: • It is typically more efficient because it operates with smaller matrices than a truncated finite system.
• It provides an exact asymptotic form of the bound states inside the lead.
• Its performance is not hindered by the presence of slowly decaying modes.
• It allows to reliably distinguish bound states from scattering states including the situation when the scattering states exist at the same energy as a bound state.
The outline of this article is as follows. Sec. 2 introduces the generic model and the formulation of the bound state problem. The algorithm is developed is Sec. 3. The first application, a quantum non-homogeneous discretized billiard, as opposed to continuous ones like in Ref. 19, is presented in Sec. 4 where we study the difference between integrable and chaotic billiards and show that in the former there exist bound states in the continuum (BICs) as in Ref. 20. (Indeed, the ability of our algorithm to isolate bound states from the continuum is one of its strengths.) Sections 5.1, 5.2, and 5.3 continue with further applications: the calculation of edge states for three different kinds of topological phases. These are, respectively, a 1D Majorana bound state in a superconducting wire, a 2D quantum spin Hall phase within the Bernevig-Hughes-Zhang (BHZ) model, and Fermi arcs in a 3D Weyl semimetal.
2 The bound state problem 2.1 General model A typical system of interest is sketched in Fig. 1: a central system of N sr orbitals is connected to one or more leads. The leads themselves are semi-infinite and periodic. They consist of unit cells that contain N t orbitals each and are repeated up to infinity. Without loss of generality, we can assume that there is only one lead: if there are several, they all can be considered as a single effective lead that is made up of disconnected parts. To simplify the notation, we include the first unit cell of the lead in the scattering region. With these conventions, the total Hamiltonian of the system can be written aŝ where the Hamiltonian of the scattering region is a finite but potentially big N sr × N sr matrix H sr , the Hamiltonian H acting within each unit cell of the lead is a N t × N t matrix, while the N t × N t hopping matrix V describes the coupling of neighboring unit cells. The N t × N sr diagonal rectangular matrix P sr is defined as [P sr ] ij = 1 when the site i of the first unit cell of the lead is coupled to the site j of the scattering region and zero otherwise. In the case where the system has several leads, the matrices H and V are block-diagonal with each block corresponding to a single lead. We search for the eigenstatesψ of the matrixĤ tot , i.e. the solutions of For finite-sized problems, this amounts to diagonalizing a finite Hermitian matrix. For infinite systems, however, two cases arise depending on whether the energy E lies within one of the bands of the lead or whether it corresponds to a localized eigenstate. The first case-the scattering problem-has been extensively studied in the literature in various formulations, see Refs. 2,21-23, and can be cast into solving a set of linear equations [5] . Because the scattering problem has a full set of solutions for any E belonging to the continuum spectrum of the lead, the energy E is an input of the calculation. In contrast, in the second case-the bound state problem-the energy is an output of the calculation since localized eigenstates exist only for specific values of E.

Lead modes
In the spirit of the scattering problem let us first analyze the possible modes: states that exist in the lead at a given energy. Since the lead is translationally invariant, any wave function can be decomposed into the eigenstates of the translation operator, where λ = e ik characterizes the momentum in the lead and j labels its unit cells (the index grows with increasing distance from the scattering region). The Schrödinger equation in the lead takes the form or alternatively where we have introduced the Bloch Hamiltonian Introducing ξ ≡ λφ we transform the Eq. (3) the generalized eigenstate problem that can, be solved using standard linear algebra routines (in practice this step requires some care and will be discussed extensively in Ref. 24). Solutions of the Eq. (5) with |λ| = 1 (real values of the momentum k) are propagating while those with |λ| < 1 are evanescent. Finally, modes with |λ| > 1 diverge exponentially with the distance from the scattering region and therefore they play no role in the bound state problem.

Formulation of the bound state problem
Inside the lead an eigenstateψ can be expressed as a superposition of propagating and evanescent modes. We can now formulate the bound state problem as follows: a bound state is an eigenstateψ that has no overlap with the propagating modes, i.e. that is purely evanescent.
To be more specific, we gather all the eigenvalues |λ| < 1 at a given energy E into a diagonal matrix Λ e (E) (of size N e ×N e ) and the corresponding evanescent eigenstates φ into the matrix Φ e (E) (each column of the N t × N e matrix Φ e is a vector φ corresponding to an evanescent state). With these definitions Eq. (3) takes the form Using this notation, an eigenstateψ assumes the following general form in the lead: where the vector q e contains the coefficients of the expansion in terms of the evanescent states. Denoting the subvector ofψ inside the scattering region ψ sr , we arrive at the following equation with unknown (ψ sr , q e ): Here we have reinstated the explicit energy dependence of the mode matrices, and this result is very similar to Eq.D2 of Ref. 25, but the authors did not look for a way to formulate this nonlinear problem in an efficient and practical algorithm. We therefore reduce the bound state problem to finding E, ψ sr , and q e satisfying the Eq. (8).
Eliminating q e from the Eq. (8) yields an equivalent equation where the self-energy Σ is In the remainder, we focus on the formulation Eq. (8) and do not make use of the alternative self-energy formulation. Note that Eq. (10) is only defined when the matrix V is invertible. A better definition of the self-energy is given in the next section.

Numerical algorithm for finding bound states
We now turn to the construction of practical algorithms that use Eq. (8) to calculate the bound states of an infinite system. For completeness and pedagogical purposes, we present three different algorithms of increasing effectiveness. Algorithm III is numerically superior to the other two.

Algorithm I: non-Hermitian formulation
The first algorithm is formulated directly in the non-Hermitian representation of Eq. (8). The first thing to notice is that when the matrix V is not invertible, it admits trivial solutions.
As an example, let us consider One immediately sees that this leads to the last row of Eq. (8) to be made up of zeros, so that the equation admits solutions for any energy E. However, while one can find a non-zero vector q e , it corresponds to a vanishing eigenvalue λ for a vector φ that belongs to the kernel of V (i.e. V φ = 0) so that the actual state given by Eq. (7) vanishes everywhere and is therefore not a true solution.
To remove these spurious solutions, we perform a singular value decomposition of V that we write as We order the matrices D and Λ e such that their vanishing eigenvalues are placed in the lower right part of the matrix. It can be noted that the number of eigenvalues λ = 0 is equal to the dimension of the kernel of V , as seen form Eq. (3). Discarding the zero-valued trailing rows of Eq. (8) and an equal number of columns, that are either made of zeros or correspond to modes associated with λ = 0 contributions, we arrive at whereΦ e ,Λ e ,q e ,D andW † V are the truncated versions of Φ e , Λ e , q e , D and W † V , where the elements that play no role are removed. The matrix in Eq. (11) has the size (N sr + N V ) × (N sr + N e ), where N V is the rank of V . Note that the matrix is square if and only if there are no propagating modes.
The bound state problem is now reduced to finding the values of E for which the left-hand side Q(E) of Eq. (11) is not invertible. However, Q(E) is in general not Hermitian and not even a square matrix. Singular value decomposition provides a solution to this problem: we form the matrix Q † (E)Q(E) and obtain its lowest eigenvalues using a sparse solver (we use the Arpack implementation of the Lanczos algorithm using the shift-invert technique). Since the matrix Q † Q is positive definite its eigenvalues are also positive (see an example in Fig. 2b), we are hence looking for zero singular values usind standard one-dimensional minimization algorithms. They are however less efficient than the root-finding algorithms that we apply in algorithms II and III.
Algorithm I is similar to the algorithm introduced in Refs. 26,27 which maps the bound state problem on the solution of a non-linear eigenvalue problem of a non-Hermitian matrix. There is, however, a notable difference: Refs. 26,27 looks for a vanishing determinant instead of the lowest singular value. This is less numerically efficient since determinant calculations can easily overflow/loose precision for large matrices and they do not exploit the matrix sparsity structure.

Hermitian formulation
To improve on algorithm I, we slightly reformulate Eq. (8) which is in general overcomplete since the matrix on the left-hand side is (N sr + N t ) × (N sr + N e ). By multiplying the second line of Eq. (8) by Λ * e Φ † e we obtain a set of necessary conditions for a bound state to occur: with the explicit dependence on energy E removed for clarity. The matrix on the left-hand side of Eq. (12) now has a square (N sr + N e ) × (N sr + N e ) shape and is moreover Hermitian, a desirable property for numerical purposes. Indeed, as we show in App. B, Eq. (6) leads to and therefore to Hermiticity of Eq. (12). The next step is to get rid of the spurious λ = 0 solutions that can be present in Λ e . We reorder the Λ e matrices so that their vanishing eigenvalues are placed in the lower right part of the matrix. After this reordering the corresponding last rows of Eq. (12) vanish, and we therefore remove them from the matrix. Similarly, by using Eq. (13), we also remove the last columns of the matrix. Introducing truncated quantitiesΦ e ,Λ e , andq e where the zero entries due to the λ = 0 contributions have been disregarded, we finally arrive at the central result of this article, where we show the explicit energy dependence to emphasize the non-linearity of the eigenproblem. Eliminatingq e from the Eq. (14) we also obtain an alternative expression of the self-energy that is defined even when V is not invertible: Combining this expression with Eq. (9) provides an equivalent alternative to the algorithms presented below. We do not pursue this route because there is no a-priori indication of any benefits from using the self-energy formulation, either by solving directly Eq. (9) or using an iterative scheme.

Algorithm II: root finder in the Hermitian formulation
The problem is now set in a Hermitian form suitable for numerical computation. Let us abbreviate Eq. (14) as with and The E-dependence has not been written explicitly in Eq. (17), but the solutionsΛ e (E) and Φ e (E) of Eq. (6) are nontrivial functions of E, which makes the eigenproblem a nonlinear one. We are interested in finding the values of E and the corresponding eigenvectors at which Eq. (16) admits nontrivial solutions. A necessary condition for this is the presence of a zero eigenvalue of H eff (E). For any given E, the eigenvalues of H eff (E) that are close to zero can be computed using a sparse solver (again, we use the Arpack implementation of the Lanczos algorithm using the shift-invert technique). The values of E where any eigenvalue vanishes can be then found using one of the many standard root-finding algorithms for one-dimensional functions.
Once a vanishing eigenvalue has been found, a check is necessary in order to avoid a solution that is not a physical bound state. A trivial example of such a false positive is an eigenstate of H sr at an energy E such that there are no evanescent states. In that case the matrices Φ e and Λ e are empty and Eq. (16) is nothing but the Schrödinger equation for the scattering region alone. Hence, once a candidate solution has been found one checks that V P sr ψ sr − V Φ e q e = 0 (19) to verify that the solution indeed satisfies the original set of equations. Figure 2a shows an example of the eigenvalues of H eff (E) as a function of E for an integrable billiard [the specific Hamiltonian is introduced later in Eq. (23)]. The vanishing eigenvalues correspond to possible bound states solutions, while the arrows indicate those that are true bound states.

Algorithm III: improved convergence thanks to gradient calculation
To improve the convergence of the root finder algorithm we expand H eff (E) up to a linear order in E. Since H eff (E) is a Hermitian matrix, we use first-order perturbation theory to obtain the gradient of its eigenvalues with respect to E, where ε α is an eigenvalue of H eff and ψ eff,α the corresponding eigenvector. Knowing the gradient allows to use root-finding algorithms that converge much faster, such as the Newton-Raphson method. When implementing this algorithm, one needs to evaluate dH eff /dE which amounts to calculating the derivatives dΛ e /dE and dΦ e /dE. These matrices are built from the non-Hermitian generalized eigenproblem (5). To compute these derivatives, we follow Ref. 28 for a general eigenproblem of the form where A and B are the matrices on the left and right-hand side of Eq. (5), x = φ ξ and κ = 1/λ. We assume the eigenvalue κ to be non-degenerate. Since the matrix B does not depend on energy, taking the first derivative of Eq. (21) gives The left-hand side of the above equation can be rewritten as an extended matrix-vector product of the form This system cannot be solved directly since there are 2N t +1 unknowns but only 2N t equations. An additional constraint arises from the normalization of the eigenvector x. We choose to set the value of the largest component of x to unity: writing m for the index of the largest component of x, we have x m ≡ 1 therefore dxm dE = 0 and we can remove the corresponding m-th column from the left-hand side of Eq. (22). We are left with a system which, unless κ corresponds to a degenerate eigenvector, has linearly independent columns [28] and can be solved numerically. Figure 2 shows the eigenvalues of Eq. (14) and the lowest singular values of Eq. (11) for an integrable quantum billiard (to be introduced in the next section) as a function of E. At energies below the band bottom, when there are no propagating modes (left of the vertical dashed line), we observe a one-to-one correspondence between eigenvalues going through zero, zero-valued minima of singular values, as well as eigenstates of a truncated system. Due to the finite size of the truncated system there is a small mismatch between its eigenenergies and the zero crossings/minima. For energies (right of the vertical dashed line), the scattering region is connected to a continuum of states. The transition is marked by a discontinuity in the eigenvalues of Eq. (14). In this regime there exist eigenstates of the truncated system that do not match a vanishing singular value. These states correspond to the continuum and more of them appear as more lead cells are included in the truncated system. The eigenvalues of Eq. (14) have more zeros than the singular values of Eq. (11) because some solutions of Eq. (14) are false positives that can be eliminated using Eq. (19).

Visualization of algorithms I, II and III
The convergence rate of the different approaches is shown in Fig. 3, which confirms that algorithm III is the fastest. Since perfect convergence is achieved after only a few iterations, this algorithm is capable of obtaining the bound states for a cost comparable to that of calculating the propagating modes (scattering matrix).

Application to quantum billiards
We consider a circular (integrable) quantum billiard discretized on a square lattice with nearest neighbor hoppings:    Figures 5a) and 5b).
Here, i, j stands for a summation over nearest neighbours. The summation in the last term is restricted to the scattering region, applying an onsite potential v g inside it, while this potential is equal to 0 in the lead.
We compute the density of states (DOS) of this system using Kwant [5] , and the bound state spectrum using our algorithm. The results are shown in Fig. 4. We observe that there is an energy range in which bound states coexist with the continuum without mixing with it. At a certain energy a second propagating channel opens in the lead and washes out the remaining bound states. Below the bottom of the band, only bound states are present in the system. We emphasize that the bound states inside the continuum (BICs) present in the intermediate energy range 0.03 ≤ E ≤ 0.15 are very difficult to observe with a direct diagonalization of a finite system, since there is no simple way to distinguish eigenstates that originate from the continuum from actual bound states. Two examples of wavefunctions corresponding to particular bound states are shown in Fig. 5.  Fig. 4. State (a) is a bound states that coexists with the continuum. State (b) was selected at an energy below (but very close to) the opening of the first mode. This way, the bound state decays only slowly in the lead and would be costly to compute by diagonalizing a truncated system.
The presence of the BICs is a consequence of the symmetry of the circular billiard and of the corresponding lead: the lowest propagating mode is symmetric with respect to the y-axis while the BICs are antisymmetric, as shown in Fig. 5a), hence the absence of hybridization. On the other hand, the BICs do not appear in the chaotic billiard of Fig. 6 which does not have reflection symmetry.

Topological materials
We now present applications to topological materials where our algorithm is used to illustrate the bulk-boundary correspondence. In these three examples, we consider infinite systems of different dimensionality that are cut at x = 0 to form, respectively, a semi-infinite wire (with 0D Majorana), a half plane (with 1D edge state), and half space (with a 2D Fermi arc surface state). In all these examples we consider a lattice termination boundary with H sr = H, and P sr = 1.

One-dimensional superconducting nanowire, Majorana bound states
The starting point for the topological superconducting nanowire is the Bloch Hamiltonian from Ref. 11: where the Pauli matrices (τ 0 , τ x , τ y , τ z ) and (σ 0 , σ x , σ y , σ z ) act in the Nambu (electron-hole) space and spin space, respectively. The parameter ∆ is the superconducting gap, µ the chemical potential, α so the strength of the spin-orbit coupling and B the Zeeman magnetic field. Equation (24) corresponds to a tight-binding Hamiltonian with N t = 4, In Fig. 7 we compare diagonalization of a finite system with the output of our algorithm. The upper panel shows the spectrum of a superconducting wire of finite length as a function of the chemical potential. The region with a dense set of energies corresponds to the continuum while the Majorana bound states are isolated. The lower panel displays the bound states as calculated with our algorithm and the DOS obtained with Kwant [5] . Both panels match perfectly, but the bound state algorithm has the advantage of working directly in the thermodynamic limit which enables it to distinguish bound states from the continuum. In Fig. 8 we compute the wave function ψ † sr ψ sr at the boundary of the topological region. Close to the quantum phase transition between the topological and non-topological phases the brute force diagonalization exhibits finite size corrections.
Another application of the bound state algorithm to the same Hamiltonian, illustrated in Fig. 9, is the calculation of the Andreev bound states of a Josephson junction. The hopping term between the left superconductor and the normal metal is described by Eqs. (25) and (26) but acquires a phase difference ϕ such that The normal section in the center consists of one site where the gap ∆ = 0 and an onsite potential V N τ z is added. The Fig. 9 shows the energy of the Andreev bound states as a function of the phase difference ϕ in the topological (a) and trivial (b) case. As expected, one recovers the 4π and 2π periodicity of the respective spectrum.

Quantum spin Hall effect
We continue with the two-dimensional BHZ model for the quantum spin Hall effect [29] : where h(k) = (k) + d(k) · σ, x + k 2 y )).
The vector σ = (σ x , σ y , σ z ) contains the Pauli matrices. The constants A, B, C, D and M are material parameters. The tight-binding model is two-dimensional with the onsite Hamiltonian and hopping matrices Aσ y along the two directions. We apply the Bloch theorem in the x-direction, and compute the bound state spectrum and the DOS as a function of k x . We calculate the bound state of the effective quasi-onedimensional system (parametrized by k x ) given by for which our algorithm can be applied directly. The results are shown in Fig. 10 where the DOS in the topological insulating phase is shown together with the positions of the bound states. As expected, the helical edge states appear in the gap.

Fermi arcs in Weyl semimetals
We conclude with a last example that uses the same procedure for a three-dimensional Weyl semimetal whose Bloch Hamiltonian is given by [30] where µ(k) = µ 0 + t(2 − cos(k x ) − cos(k y )) + t z (1 − cos(k z ))). This Hamiltonian corresponds to the 3D tight-binding model Once again, applying the Bloch theorem in the x-direction and using k y and k z as parameters, so that which can be used directly in our algorithm. The results are shown in Fig. 11 and 12. As expected, we obtain the Fermi arcs that couple the two Weyl points. Again, our algorithm allows to study a single surface of the Weyl semimetal, in contrast with the finite size calculations (shown in the lower panel of Fig. 11) where states belonging to the two surfaces hybridize.   Figure 12: Same as the upper panel of Fig. 11 with an alternative view: the energy E = 0.3 is fixed and the DOS is plotted as a function of (k y , k z ).

Conclusion
We have derived a robust and efficient algorithm to compute the energies and the wavefunctions of the bound states of infinite quasi one-dimensional systems described by general tightbinding Hamiltonians. We have applied our approach to a variety of systems and shown that it can efficiently calculate the bound states, including the situations where other approaches would fail or become computationally prohibitive.
This algorithm can be used on its own to obtain e.g. the Josephson relation of superconductingnormal-superconducting junctions, the properties of quantum wells or characterize topological materials and their interfaces. It may also be used to obtain the properties of perfect surfaces or interfaces (in arbitrary topological, metallic or superconducting materials) that can then be used as the starting point for further calculations.
Funding information Financial support from US Office of Naval Research (ONR) is gratefully acknowledged. AA and MW were supported by the Netherlands Organisation for Scientific Research (NWO/OCW). AA was also supported by the ERC Starting Grant 638760.

A Normalization of the bound state
Bound states should be correctly normalized: However, our algorithm does not ensure this normalization automatically. For a given set (ψ sr , q e,α ) we find that We recognize a geometric series and arrive at with the matrix N defined as Eq. (37) can now be used to normalize the bound state wave functions to unity.

B Proof of Eq. (13)
To prove Eq. (13), we begin by multiplying Eq. (6) by Φ † e : The complex conjugate of the above equation reads Now, multiplying Eq. (39) by Λ * e on the left and Eq. (40) by Λ e on the right, we arrive after substracting one equation from the other at Since we are dealing with evanescent states, we can simplify by 1 − λ * α λ β = 0 and arrive at which is essentially Eq. (13).

C Degenerate eigenvalues
In some cases, the solutions λ of Eq. (3) can be degenerate, as in the quantum spin Hall model of Sec. (5.2) where the degeneracies arise because of the two species of spin. A set of degenerate eigenvectors is not uniquely defined, as any linear combination is still a valid solution of the eigenproblem (3). This means that the matrix Φ e that was introduced in Sec. 2.3 can be replaced by Φ e ≡ Φ e T , where T is an invertible matrix equal to unity except in the block corresponding to degenerate eigenvalues. This modification results in an uncertainty of the l.h.s. of Eq. (14) since does not have the same eigenvalues as H eff , unless T is unitary. (We use the fact that Λ e commutes with T .) This leads to a problem during the root-finding phase of the algorithms of Sec. 3.1 and 3.3: in a naive implementation H eff is evaluated multiple times for different E, each time with an effectively random T . The resulting fluctuations, shown in Fig. 13(a), are incompatible with efficient root-finder routines. The algorithm presented in 3.4 is not considered here as the derivatives computed in that section are not well defined for degenerate eigenvalues.
It is important to note that the transformation from H eff (E) to H eff (E) does not affect the zero-valued eigenvalues -the only ones with a physical meaning. Furthermore, replacing Φ by ΦT also changes q e by T −1 q e , so that the final wavefunction Ψ(j) in Eq. (7) remains unchanged. This is why the algorithms of Sec. 3 remain correct, but require a fluctuationtolerant root-finder when implemented naively. In practice, it is preferable to fix a unique spectrum such that an efficient standard root-finder can be used. The matrices Φ e cannot be fixed in a unique way, but one can choose the column vectors such that the degenerate ones are orthogonal. To this end, one performs the QR decomposition Φ e = QR, where Q is an orthogonal matrix and R an upper triangular matrix, and substitutes in Eq. (17): The resulting smooth eigenvalues are shown in Fig. 13(b).
Orthogonalizing Φ e forces the matrix T to be unitary, as we can understand using the following geometrical argument. Any superposition on two vectors e 1 and e 2 of the (x, y) plane is still a valid basis as long as they are not collinear. If one forces the vectors to be perpendicular to each other, then the only transformation left is a rotation, in other words a unitary transformation.