Non-Abelian Majorana fermions in topological $s$-wave Fermi superfluids

By solving the Bogoliubov--De Gennes equations analytically, we derive the fermionic zero-modes satisfying the Majorana property that exist in vortices of a two-dimensional $s$-wave Fermi superfluid with spin-orbit coupling and Zeeman spin-splitting. The Majorana zero-mode becomes normalisable and exponentially localised to the vicinity of the vortex core when the superfluid is topologically non-trivial. We calculate the energy splitting due to Majorana hybridisation and identify that the $s$-wave Majorana vortices obey non-Abelian statistics.

Introduction -In two dimensions, the quasi-particle excitations of topologically ordered systems are generally non-Abelian anyons [1]. Possible realization of non-Abelian statistics has been studied in connection with the ν = 5/2 fractional quantum Hall (FQH) state [2] and the vortex state of chiral p x + ip y superconductors and superfluids [3,4]. Majorana fermions in the cores of superfluid vortex excitations have been actively considered in the context of p-wave pairing [5,6], but can also occur in an s-wave superconductor coupled by the proximity effect to a topological insulator [7]. Majorana vortices in the topological s-wave superconductor are non-Abelian, and in the same topological class [1] with the Moore-Read Pfaffian FQH state [8,9], p x + ip y superconductors [10,11], and the gapped non-Abelian spin liquid phase of the Kitaev model [12]. A key benefit of non-Abelian exchange statistics is that the Majorana vortices can be used to realize braiding operations. Noncommutative braiding, which amounts to exchanging two Majorana vortices adiabatically [13], is a key ingredient needed for performing quantum logic operations on a fault-tolerant topological quantum computer [14].
In ultra-cold atom experiments [15][16][17], however, the densities and temperatures are low, and the scattering between fermions takes place typically in the s-wave channel [18]. Experimentally, the ultra-cold atomic swave Fermi gas is a highly flexible quantum many-body system whose interactions, spin balance, and trapping geometries can be tuned nearly arbitrarily [19,20]. In the mean-field picture, it is possible to achieve band inversion leading to a topological phase [1,21] by combining the effects of two-dimensional spin-orbit (SO) coupling [22,23], and spin imbalance [24,25]. Recently, tunable 2D SO coupling was demonstrated experimentally [26]. These advances open up a promising perspective for creating a highly-controlled atomic topological superfluid [27,28] in the laboratory in the near future.
Despite the promising recent experimental progress, Majorana zero-modes in vortices of an s-wave Fermi superfluid with SO coupling and a Zeeman field remain relatively poorly understood. Here, we derive analytically the Majorana vortex zero-mode, and use * lauri.toikka@gmail.com it to calculate the energy splitting due to inter-vortex tunnelling. The tunnelling generally lifts the zeromode degeneracy [29], and has been studied for the Moore-Read state [30], Kitaev's honeycomb model [31], and p-wave superconductors [32][33][34]. Hybridisation of Majorana fermions in dense vortex lattices gives rise to a band structure [35][36][37]. Of particular interest for experimentally controlled quantum simulation are topologically non-trivial bands of Majoranas [38], and the possibility of flattening them [39].
Superfluidity in spin-balanced Fermi gases results from the formation of Cooper pairs, bound states of two fermions around the Fermi surface with opposite momenta +k and −k. For spin-imbalanced pairing [41], theoretical predictions such as FFLO [42,43] and Sarma phases [44] together with deformed Fermi surface super-fluidity [45] have been presented. In the FFLO phase the Cooper pairs have a non-zero center-of-mass momentum which gives rise to a periodic spatial modulation of the order parameter and also possibly of density. While spin-imbalance can thus influence the pairing and subsequently the dynamics of topological defects by frustrating the Cooper pairing, in experiments the excess nonsuperfluid particles are typically spatially separated from the completely paired BCS superfluid [25,46,47].
In typical current experimental conditions in ultra-cold atoms, the pairs form in the s-wave state with zero angular momentum. The Pauli principle then necessitates that the spin state be a singlet. Since the Cooper pairs have no angular momentum and no net spin, the relative pair wavefunction can be fully characterised by a single complex amplitude ∆(r, t).
Ignoring now all the other quantum correlations apart from the pair correlations captured by the order parameter, the mean-field Hamiltonian describing quasi-particle excitations in the superfluid Fermi gas without spin-rotational invariance reads is the Bogoliubov-De Gennes (BdG) matrix representation. In what follows, we are not interested in explicit time dependence and therefore drop the labels.
In the mean-field picture, we replace the many-body problem (1) with the effective single-particle Hamiltonian (2) parametrised by the pair potential (order parameter) ∆(r). The pair potential together with the chemical potential is determined by optimising the effective Hamiltonian (2) such that it minimises the total free energy. The result for the pair potential is and for the atomic densities where f ν = 1/ e Eν /(kBT ) + 1 is the Fermi distribution for occupations and T is the temperature. We set T = 0. We have expressed the order parameter in terms of the amplitudes u and v, defined by the Bogoliubov-Valatin transformationsψ σ (r) = ν u ν,σ (r)γ ν + v * ν,σ (r)γ † ν , γ † ν = dr σ u ν,σ (r)ψ † σ (r) + v ν,σ (r)ψ σ (r) , which diagonalise the effective single-particle mean-field Hamiltonian.
Symmetries of the Hamiltonian -To obtain the static properties of Majorana fermions at the vortex cores in a topological s-wave Fermi superfluid, we need to solve the fundamental eigenvalue equation for the zero-mode ν = 0 with E 0 = 0. When considering equations of the form (5), it is instructive to first understand the symmetries of the underlying Hamiltonian. The BdG Hamiltonian H(r, t) possesses a particle-hole symmetry (PHS) defined by {Ξ s−wave , H(r, t)} = 0 with Ξ s−wave = e iθ τ x ⊗ 1 σ K, where K is the complex conjugation operator (in momentum space, K additionally flips the sign of the momentum), τ y (σ z ) is the Pauli matrix in particle-hole (spin) space, and θ ∈ R. As a result of the particle-hole symmetry, if Additionally, if and only if ∆(r, t) ∈ R and h = 0, the BdG Hamiltonian possesses a time-reversal symmetry defined by [T, H(r, t)] = 0, where T = e iθ 1 τ ⊗ σ y K (T 2 = −1, for half-integer spin), and θ is arbitrary. We have defined 1 τ to be the identity matrix in the particle-hole space, while σ y is the Pauli matrix in spin space.
A straightforward calculation shows that nondegenerate zero-modes (that is, the modes H(r)h 0 = E 0 h 0 with E 0 = 0), must always satisfy the important symmetry which is the so-called Majorana property making the zero-mode special. The Majorana property ensures that the quasi-particle operatorŝ satisfyγ † 0 =γ 0 ; they represent Majorana fermions. Let us temporarily assume that time-reversal symmetry is respected, i.e. ∆(r) ∈ R and h = 0. If T is a solution of Eq. (5) with energy E ν , then T h ν is a solution with the same energy E ν . This is the pair of the zero-energy Majorana mode. With time-reversal symmetry, the Majorana probability densities are identically overlapping in space: We need to break the time-reversal symmetry to separate them spatially, for example, by having a vortex in the order parameter.
A broken time-reversal symmetry T , a broken spinrotational symmetry (due to spin-orbit coupling), and an unbroken PHS Ξ s−wave result in the BdG Hamiltonian H(r) belonging to the symmetry class D in the Altland-Zirnbauer classification [48]. According to Ref. [49], point defects (such as the vortices in the order parameter considered here) with a Hamiltonian that belongs to the class D are associated with a Z 2 -valued topological invariant. This means that in the topologically non-trivial regime the number of exact Majorana zero-modes is given by W mod 2, where W is the sum over all vortex winding numbers, with the rest hybridising in pairs forming a band structure. Generally, when multiple vortices are brought close together, the Majorana zero-modes hybridize into a band structure, leading to complex fermion states at positive and negative energy [35][36][37]. We calculate the energy splitting in Eq. (18). For periodic vortex lattices, such as the square and the triangular lattice, the resulting low-energy theory is typically gapped and topologically non-trivial [36].
Eigenvalue equation for the zero-mode -We now include a vortex in the order parameter, where r and ϕ are polar coordinates centred on the vortex. Energy considerations require that ∆(0) = 0 at the vortex core. The integer denotes the vorticity, and ∆(r) is a real function of r that vanishes at r = 0. In what follows, we will solve Eq. (5) analytically for the zero-mode h 0 that exists in the core of the vortex (8).
For convenience, let us first apply the unitary transformation U s , which transforms the basis as shuffles the zeros in the rows and columns of the Hamiltonian H(r) [Eq.
(2)] to give where D σ = diag(K σσ , −K σσ ) and Setting D σ = 0 decouples the spins, and the Hamiltonian reduces to a 2 × 2 BdG matrix corresponding to the Fu-Kane model at neutrality, whose Majorana modes have been considered [50]. This is the low-energy regime near µ σ = 0, which is the topological transition point in typical p-wave systems. The Majorana zero-modes have also been considered for the Fu-Kane model at µ = 0 [33].
However, here we consider the full 4 × 4 structure of Eq. (5):K When looking for zero-modes, the system (11) can be simplified by using the symmetry property (6), that is, we seek simultaneous eigenstates of the PHS. There are two possibilities for θ, namely θ = θ 0 , θ 0 + π, which both satisfy the Majorana relation. This means that with the final solution we must allow for both possibilities. In both cases, we obtain a coupled 2×2 system, which must be solved for u 0,↑ and u 0,↓ . In polar coordinates, the spin-orbit coupling terms Observing the azimuthal symmetry, let us try a separable ansatz with angular momentum eigenstates, with the assumption that ∆ has no r-dependence. Physically, this approximation treats the vortices as points with only a phase profile. Substitution into the system (12) shows that we can eliminate the angular dependence by taking m The only requirement here is that ∈ Z. Explicitly, the coupled system (12) then reads Below, we solve analytically the coupled system (14) for F ↑ (r) and F ↓ (r). Series solution for F σ (r) -The vortex core (r = 0) is an irregular singular point of Eq. (14) where numerical methods are inherently unstable. An obvious solution is to start the numerical integration at a point r 0 > 0 thus avoiding the singularity, but in this case we need to know accurately the initial condition at the rather arbitrary point r 0 . It is therefore highly desirable to have analytic solutions, at least around the vortex core. In what follows we solve the system (14) using the Wasow method [51], a generalisation of the Fröbenius series method for coupled systems of differential equations with irregular singular points.
The details for solving the system (14) are shown in Appendix A. For definitess, we set = +1. The result is a series solution for F σ (r) around the origin that can be evaluated analytically to arbitrary order, and reads While F ↑ (0) can be finite, the boundary condition F ↓ (0) = 0 is forced to keep F ↓ finite at the origin, which results from the series solution for F ↓ starting with the power r −1 whereas that for F ↑ starts with r 0 . That one spin component is zero and the other finite at the vortex core has been observed in all numerical studies of the Majorana zero-mode in s-wave Fermi gases [52,53]. In principle the system (14) needs four initial conditions, but requiring that F ↓ be finite at the origin consumes two of them leaving only the two initial conditions F ↑ (0) and F ↓ (0) in the solution (15).
Existence of the Majorana zero-mode -The existence of the analytical series solutions (15) alone does not guarantee the existence of a topologically protected Majorana zero-mode. The condition for the physical existence of the zero-mode is that it remain normalisable as r → ∞.
In the limit r → ∞, the kinetic energy and terms in 1/r can be dropped in system (14) to obtain Physically, we are interested in only the solutions for system (16) that are normalisable. For example, if λ > 0 and ∆ < 0, the only normalisable asymptotic solution is where C is a constant. Even then, the solution (17) is normalisable only when h > μ 2 + |∆| 2 . The value at equality corresponds to the critical Zeeman field marking the phase transition to the topological regime with a uniform order parameter, ∆(r) = ∆ [54], which here is associated with the boundedness of the Majorana zero-mode at infinity. However, it is not straightforward to find initial conditions F ↑ (0) and F ↓ (0) such that the zero-mode has the desired large-r asymptotics.
In principle,μ and ∆ must be obtained selfconsistently in conjunction with solving the BdG eigenvalue equation. Considering a uniform vortex-free system, we can solve Eqs.  [54]. From Eq. (17), therefore, we know that with these parameters there is an exponentially decaying asymptotic solution, and we have performed a direct search for F ↑ (0) and F ↓ (0) that match with this asymptotic behaviour.
We retain terms upto O(r 10 ), and use the analytical solution (15) to evaluate accurate initial conditions for a numerical integration of the system (14). The truncated analytical series solution agrees with the numerical integration for r 0.5 (Fig. 1). The analytic series expansion can be developed to arbitrary precision. As expected, the numerical solution becomes unstable as the vortex core is approached, but more importantly it can be used to study the boundedness of the mode as r → ∞. We find a bounded zero-mode in vortices with unit positive circulation with the specific parameter values above for 3.488885 < F ↑ (0) < 3.488886 and F ↓ (0) = −12.518. The full solution including azimuthal dependence is then given by Eq. (13).
Majorana energy splitting -We now use the Majorana zero-mode to calculate the energy splitting in an s-wave Fermi gas that results from hybridisation between the Majorana modes of two vortices located at points r 1 and r 2 . Then, ∆(r) ≈ ∆ i (r) near vortex i (i = 1, 2), where ∆ i (r) is the single-vortex pair function (8) for vortex i. Introduction of a single-vortex background phase θ i , ∆ i (r) = ∆e i( iϕi+θi) with ϕ i the azimuthal angle measured with respect to vortex i, amounts to the Majorana mode changing by u 0,σ (r − r i ) → e i θ i 2 u 0,σ (r − r i ), leaving invariant the system (12).
In the tight-binding approximation, the Hamiltonian H ∆ describing hopping between the two vortices, with a straightforward generalisation to lattices of Majorana vortices, is given by H ∆ = i t 12γ0,1γ0,2 , wherê γ 0,i is the Majorana mode (7) at vortex i, and t 12 = g (2) 0 |H (Us) |g is the overlap integral that gives the energy splitting. Here g the Majorana zero-mode centered at vortex i.
, by definition of the zero-mode at vortex 1 we obtain Compared with p-wave pairing [55], under Majorana ex-change t 12 is now symmetric with respect to θ 1 − θ 2 , but still anti-symmetric due to the spin degree of freedom. It was pointed out by Fujimoto [1] that the existence of the Majorana zero-mode with SO coupling does not automatically guarantee their non-Abelian statistics. Majorana interchange 1 ↔ 2 gives t 12 = −t 21 , and the U(1) gauge transformation θ i of the Majorana fermion has the important property that when θ i changes from 0 → 2π, the Majorana changes sign,γ 0,i → −γ 0,i . Braiding of vortices i and j changes the superfluid phase at one vortex by 2π amounting toγ 0,i →γ 0,j ,γ 0,j → −γ 0,i , and it was shown by Ivanov [10] that this property together with quantisation of i gives rise to non-Abelian statistics.
Conclusions -We have derived analytically the Majorana zero-mode in a topological Fermi superfluid with s-wave pairing, two-dimensional spin-orbit coupling, and a Zeeman field. We find an exponentially localised Majorana zero-mode at the vortex core only in the topologically non-trivial regime of the superfluid. We find that in the s-wave Fermi superfluid the Majorana fermions obey non-Abelian exchange statistics, and the energy splitting due to Majorana hybridisation is determined by both spin sectors. Knowing the Majorana zero-mode analytically paves the way for studies of quantum many-body correlations and simulation of topological quantum matter in lattices of Majorana vortices in a new experimentally well-controlled setup.

ACKNOWLEDGMENTS
I would like to thank Andreas Läuchli for fruitful discussions. This work was supported by the Austrian Academy of Sciences (P7050-029-011).
Introducing the first derivatives x = F ↑ and y = F ↓ as auxiliary variables, we can trivially regroup the system (14) for x and y in terms of x, y, F ↑ , F ↓ . This gives the equivalent matrix equation where where but this approach suffers from the irregular singular point at the origin as well.
Near r = 0, we can always write where M ν are constant 4 × 4 matrices holomorphic at r = 0 for which the series converges component-wise in a neighbourhood of the origin. Here g = 2, and M ν = 0 for only ν = 0, 1, 2.

Solution for the eigenvalue equation (14)
Generalising the vector u into the matrix X, and transforming x = 1/r changes Eq. (A6) into where q = g − 2 and with the only non-zero matrices being The matrix M(x) is holomorphic at x = ∞ meaning that there exists a convergent expansion of the form (A8) for sufficiently large x 0 such that |x| > x 0 . If the integer q + 1 > 0, the singular point is irregular, and if q + 1 = 0, the singular point is regular. For us q = 0. Our goal is to reduce Eq. (A7) (where q = 0) through formal transformations into a system with q = −1, which corresponds to a regular singular point at the origin x = ∞, and can be solved in terms of more standard Fröbenius series ansatz methods.
To form our starting point, we transform M 0 into the Jordan canonical form (JCF) (we define the JCF as having the 1's on the superdiagonal). The matrix M 0 is brought to the JCF by the non-singular constant simi-larity transformation P (0) , where for definiteness we have fixed = +1. A 0 is now in the canonical Jordan block diagonal form where H j (j = 1, 2, 3) are shifting matrices where H 1 is of dimension two and H 2 , H 3 are of dimension 1. The same transformation is applied to all the matrices, defining the new starting point where q = 0, and a. Prepare for a shearing transformation that simplifies the problem The transformation Y (x) = K(x)Z(x), where the matrix K(x) is holomorphic and has a non-vanishing determinant at x = ∞ changes Eq. (A12) into with q ≥ 0 and We seek series solutions Using Eq. (A16), substitution of the series ansätze into Eq. (A15) and comparison of like coefficients results in the recursion relation where the last term in Eq. (A17b) is absent for ν −q−1 < 0, that is, for ν < 1. The recursion relation is of the form If all the eigenvalues of A 0 are distinct, then B(x) will be diagonal and the problem is uncoupled and easily solved. If at least two eigenvalues are distinct, we can take all B ν (ν > 0) zero or block-diagonal. However, this is not possible here because A 0 has only one distinct eigenvalue, zero.
Instead, we partition each Eq. (A18) into blocks of the same order as the Jordan blocks H j for A 0 , and call these blocks K jk ν with j, k = 1, 2, . . . , s. For us s = 3. Then each relation (A18) corresponds to s 2 = 9 relations where H and K are shifting matrices of orders h and k respectively, and M is a h×k matrix whose first h−1 rows are given constant vectors while the entries α 1 , α 2 , . . . , α k of the last row are variables. It can be proven that the numbers α 1 , α 2 , . . . , α k can be determined uniquely in such a way that Eq. (A20) is solvable for the h×k matrix X. We now apply this result to Eq. (A19). Without loss of generality we take all rows but the last of B jk ν to be zero. Then, the series is determined by solving the recursion relations (A19) successively for ν = 1, 2, . . .. While it will in general be divergent, it can be proven that it is the asymptotic expansion of some holomorphic matrix function K(x) for sufficiently large x 0 such that |x| > x 0 . The point of the transformation K is to obtain a differential equation (A14) such that the matrices B, by construction, satisfy the following three properties: (i) H s (the H j are shifting matrices); and, most importantly, (iii) that the only non-zero entries in B ν for ν > 0 occur in the rows corresponding to the last rows of the Jordan blocks H k (k = 1, 2, . . . , s) in the representation of A 0 . The transformation K induces zeros into the matrices preparing them for a shearing transformation. The result is an asymptotic series solution valid at x → ∞ such that We work until ν = ν max . The final series solutions will then be exact upto ν = ν max − 1 for F ↑ (r) and ν = ν max − 2 for F ↓ (r). For spin-↑ the derivative (component 2) will agree with a direct derivative of component 1 upto ν = ν max − 2 and for spin-↓ the derivative (component 4) will agree with a direct derivative of component 3 upto ν = ν max − 3.

b. Simplify the problem by a shearing transformation
Having obtained the matrix B, Eq. (A14) is ready for the shearing transformation with a temporarily unknown positive parameter ξ, which takes Eq. (A14) into where the matrix C(x) The matrix elements c jk (j, k = 1, . . . , 4, n = 4) of C read Here δ jk is the Kronecker delta, and b jk are the matrix elements of B. The parameter ξ must be chosen appropriately to induce, where possible, non-zero elements below the main diagonal into the leading order matrix C 0 . Above the main diagonal it is equal to B 0 . The parameter ξ must be chosen judiciously as follows.
where b jkν = 0 and each positive integer α jk ≥ 1 except for the special elements with a 1 from a shifting matrix in the Jordan decomposition for B 0 for which α jk = 0. Since we could not diagonalise B 0 and instead obtained a Jordan matrix, at least one such special element is present. Before the shearing transformation (ξ = 0) the special elements have the lowest α jk , and the purpose of the transformation is to add non-zero elements on or below the diagonal to B 0 by a suitable choice of ξ. After the shearing transformation the expansions of non-zero offdiagonal elements c jk (δ jk = 0) begin with the power −α jk + ξ(j − k). The special elements on the superdiagonal begin with the power x −ξ . There exists a smallest rational ξ 0 = q (1) /p (1) > 0 with q (1) , p (1) coprime, for which the special elements have the same leading power as an element below the main diagonal, ξ 0 = α jk − ξ 0 (j − k) for some j, k < j. We take ξ = ξ 0 in the shearing transformation; the result of the above process is ξ = q (1) /p (1) with q (1) = 1, p (1) = 2 coprime. Fractional powers thus unavoidably appear in the shearing transformation (A23).
In descending powers of x −1/p (1) , the matrix C will begin with the power x −q (1) /p (1) . The matrix x ξ C as x → ∞ has at least one non-zero entry on or below the main diagonal and above the main diagonal it is equal to B 0 : We remove the fractional powers by multiplying both sides by x ξ and introducing the new independent variable x 1 such that to obtain from Eq. (A24) the differential equation where The branch of the multi-valued function x p (1) can be chosen freely.
While new elements have appeared on and below the main diagonal in D 0 , the upper triangular part above the main diagonal, in particular the superdiagonal, of D 0 is equal to that of B 0 by construction. This completes the purpose of the shearing transformation.
If h < 0, the problem would be solved because then either the singular point would be regular (h = −1), or there would not be any singular point. If h ≥ 0 and D 0 has at least two distinct eigenvalues, then the problem reduces to a set of similar problems of lower order. However, the eigenvalues of D 0 are all equal to zero, and D 0 is nilpotent. It can be proven [51] that if the process outlined above is repeatedly carried out, eventually we obtain a nilpotent D 0 that is itself a shifting matrix i.e. s = 1.

c. Obtain an equation with a regular singular point
The entire process must be reapplied, starting from the JCF for D 0 . Another shearing transformation with the parameter 1/2 is required. Repeating the process once more, the third shearing transformation comes with an integer-valued shearing parameter of 1, which means that via such a chain of transformations we have reduced the coupled system (A12), to the system and x 1 = x 2 2 /4. We have reached our goal that we set below Eq. (A7): Eq. (A32) has only a regular singular point at infinity (or the vortex core in terms of r). This equation can be solved with more standard methods.
Through a sequence of standard double transformations that raise the eigenvalues of a Jordan block of A (3) 0 by an integer followed by a non-singular constant similarity transformation of A (3) 0 to JCF, we obtain such that no eigenvalues of the leading matrix A Let us now transform back to radial coordinates with r 2 = 1/x 2 (here x 1 = x 2 2 /4) to obtain from Eq. (A34) the equation where A (7) (r 2 ) = −A (6) (1/r 2 ) and A (7) (0) is holomorphic.

A
0 K 0 − K 0 A 0 K (7) ν − K (7) ν A ν−s − A Given that the matrix A 0 in the convergent expansion A (7) = ∞ ν=0 A (7) ν r ν 2 has no eigenvalues that differ from each other by positive integers, there exists a formal convergent series K (7) (r 2 ) with K 0 = I such that the formal transformation Y (7) = K (7) (r 2 )Υ reduces the differential equation (A35) to the form of Eq. (A37) such that B (7) ν = 0 for all ν > 0. Therefore, the recursion relations (A38) can be solved for K (7) (r 2 ) such that B (7) ν = 0 for ν > 0. The solution for Y (7) (r 2 ) is then given by Eq. (A36). Here a set of four linearly independent vectors forms the fundamental matrix solution Υ(r 2 ), and Υ(r 0 ) is the fundamental matrix solution at r 2 = r 0 . We can set the auxiliary parameter r 0 = 1 and choose Υ(r 0 ) in such a way as to choose the boundary conditions for F σ (r) at r = 0. The general (vector) solution is then given by linear combinations of the columns. The logarithms of the fundamental solution matrix end up in columns 2 and 4 making them either divergent or zero at the origin r = 0. The columns 1 and 3, in contrast, can be chosen to remain well-behaved and finite at the vortex core r = 0. We do not include the columns 2 and 4 in what follows, that is, we take a linear combination only of columns 1 and 3.
Performing all the transformations carried out in an inverse order, we can compute the matrix solution X for the original problem (14) from knowing Υ(r 2 ). A judicious choice for Υ(1) (with r 0 = 1) then gives the explicit solution shown in Eq. (15) in the main text.