Berry phases and the intrinsic thermal Hall effect in high temperature cuprate superconductors

The Bogoliubov quasiparticles move in a practically uniform magnetic field in the vortex state of high temperature cuprate superconductors. Do the quasiparticles experience a Lorentz force when set in motion by an externally applied heat current ${\bf j}_Q$, bending their trajectories and causing the temperature gradient perpendicular to ${\bf j}_Q$ and the applied field ${\bf H}$, or is the thermal Hall effect a consequence of Berry phases as in an intrinsic anomalous Hall effect of a semiconductor/metal with spin-orbit coupling? Here we show that it is the latter, and for the first time, calculate the temperature, ${\bf H}$-field and the $d$-wave pairing gap $\Delta$ dependence of the intrinsic thermal Hall conductivity, $\kappa_{xy}$. We find that the intrinsic contribution to $\kappa_{xy}$ displays a rapid onset with increasing temperature, which compares favourably with existing experiments at high ${\bf H}$-fields on the highest purity samples. This finding may help to settle a much-debated question of the bulk value of the pairing strength in cuprate superconductors in magnetic field.

The low energy excitations of a superconductor, the so called Bogoliubov quasiparticles (qps), are a quantum mechanically coherent superposition of an electron and a hole. As such, they do not carry a fixed electrical charge. In the vortex state of an extreme type-II superconductor, such as high T c cuprates, the externally applied magnetic field H penetrates the sample in the form of flux tubes, whose width, set by the penetration depth λ, is much larger than the intervortex separation set by the magnetic length ℓ H = hc/(eH), for essentially the entire range of H-field applied in a typical experiment. Therefore, the qps move in a practically uniform magnetic (induction) field B. Among the best established physical properties of cuprate superconductors in the vortex state is the √ Hincrease of the low temperature (T ) specific heat when the H-field is applied perpendicular to the CuO 2 planes. This effect was predicted by Volovik 1 based on a semiclassical analysis of the qp motion in the vortex state of d x 2 −y 2 superconductor. Although such analysis successfully captures the overall √ H increase of the qp density of states, a full quantum mechanical solution is necessary in order to investigate non-trivial topological structure of the qp wavefunctions in the vortex state. Here we provide such solution. Using the Berry phases associated with the qp wavefunctions, we calculate, for the first time, the full T , H-field and the d-wave pairing gap ∆ dependence of the intrinsic contribution to the thermal Hall conductivity, thereby shedding light on physics beyond the Volovik effect.
In the thermal Hall effect a small applied energy current j E causes a transverse temperature gradient ∇T , i.e., a non-zero off-diagonal component of the Fourier law j E µ = −κ µν ∇ ν T . This effect has been reported experimentally in (slightly overdoped) YBa 2 Cu 3 O 6.99 samples by Ong and collaborators 2 , and by Zeini et.al. 3 . The measurements extend up to 15Tesla, and down to the lowest temperature of ∼ 10Kelvin. At the largest Hfields, κ xy is positive, decreases with decreasing T . However, linear extrapolation to T → 0 would yield finite negative value of κ xy , which is unphysical because κ xy must vanish in such limit. In fact, the data at the highest fields is consistent with κ xy being very small bellow T onset ∼ 7K, followed by a rapid increase at T onset with an almost linear T -dependence at higher T . New (unpublished) data at lower doping and lower temperature are consistent with such onset.
Initial theoretical arguments of Simon and Lee 4 were based on a continuum model, linearized in the vicinity of the Dirac points. They used Kubo formula to argue for the scaling form κ xy ∼ T 2 F (k B T ℓ H / v) where v = v F ∆/p F , p F and v F being Fermi momentum and velocity respectively and ∆ being the d-wave gap parameter. Based on the same linearized continuum model, Vishwanath 5 argued that the qp spectrum in the vortex state acquires a small gap which grows linearly with H, and that lim T →0 κ xy /T = nπk 2 B /(6 ) with n = 0 or ±2. However, the linearized approximation has been found to be problematic. The qp spectrum is not invariant 6,7 under the large (singular) gauge transformations 8,9 , which were being employed in the calculations 5 of κ xy . Because such large gauge invariance must be present, the validity of the results has been in doubt. A tight-binding lattice regularization was found to remove such problems 6 , and the initial model calculations 10 of κ xy for unrealistically small values of ℓ H , corresponding to over 1700 Tesla, indeed found that for a space inversion symmetric vortex lattice, the qp spectrum is generally gapped, and that lim T →0 κ xy /T = nπk 2 B /(6 ), albeit the values of n ranged from −2 to 12. Later, index theoretic arguments 11 very near half-filling found commensuration based oscillations between n = ±4 and n = 0. Finally, using the linearized model and a number of drastic simplifying assumptions, the weak field κ xy ∼ T √ H was argued be due to thermally excited nodal qps, scattering primarily from impurities, with a small skew component due to vortices 12 .
Here we build on the observation 10 that the intrinsic contribution to κ xy is a direct consequence of the qp wavefunctions acquiring a non-trivial Berry curvature 13 upon adiabatic changes of the vortex crystal momentum, k. Because j E is perpendicular to the thermal Hall gradient, the (non-dissipative) effect described here does not contribute to the entropy production, dS/dt = − d 2 r(∇T /T 2 ) · j E (r). The effect can therefore be addressed without considering irreversible processes. Without making any further simplifying assumptions, apart from treating the qps as non-interacting and the vortices forming a perfect (stationary) Abrikosov lattice, we calculate the intrinsic contribution to κ xy using the tightbinding regularization (1). We find a new scaling collapse of the extensive numerical results which clearly displays the rapid temperature onset at a small fraction of the maximum pairing gap, followed by an almost linear Tdependence.

I. MODEL HAMILTONIAN AND DENSITY OF STATES SCALING
We work on a two dimensional square lattice of spacing a -that we set to unity -and the magnetic field perpendicular to it. Our tight-binding Hamiltonian for the d-wave superconductor in the vortex state is This model has been discussed extensively elsewhere 6,11,14 : c r,σ is the electron annihilation operator, the sum over the spin projection σ =↑ or ↓ in the first and the last terms is implicit, and t r,r+δ = −te −iA r,r+δ . The chemical potential and the Zeeman coupling enter via µ ↑(↓) = µ ± h Z . The magnetic flux Φ through an elementary plaquette enters the Peierls factors as A r,r+x = −πyΦ/φ 0 , A r,r+ŷ = πxΦ/φ 0 , where φ 0 = hc/e. The ansatz for the pairing term is ∆ r,r+δ = ∆ δ e iθ(r) exp i 2 r+δ r dl · ∇θ , with the d-wave amplitude ∆x = −∆ŷ = ∆. Vortex positions r j enter through θ(r) which is chosen to be the solution of the (continuum) London's equations ∇ × ∇θ(r) = 2πẑ j δ(r − r j ), ∇ · ∇θ(r) = 0. The closed form solution for θ(r) can be found in the SI. Vortices are arranged within a unit cell as shown in the inset of Fig. 2: each magnetic unit cell, aligned with the tight-binding lattice, is threaded by magnetic flux hc/e and contains a pair of vortices. We study a variety of vortex lattices (VL): when L x = L y the vortices form a square lattice, for L x /L y ≈ √ 3 the lattice is triangular, the intermediate ratio L x /L y ≈ 1.4 yields oblique VL. Performing the singular gauge transformation 6,8,11 turns the hopping and the pairing terms in H periodic with periodicity L x and L y , allowing us to use the Bloch theorem 6,8,11 . As such, we perform the particle-hole transformation by letting c r,↑ , c † r,↓ , where ψ r,σ (k) is periodic in r → r + R and in the sum k x,y ∈ (− π Lx,y , π Lx,y ]. In terms of the new operators (suppressing k on ψ's) where t ↑↑ r,r+δ = t ↓↓ * r,r+δ = t r,r+δ e i 2 θ(r+δ) e − i 2 θ(r) , λ ↑↓ r,r+δ = ∆ r,r+δ e − i 2 θ(r+δ) e − i 2 θ(r) , andμ ↑(↓) = h Z ± µ. Further details, including the treatment of the branch-cuts in e i 2 θ(r) , are provided in SI.
Equation (2) has the form of a Hamiltonian for a fictitious semiconductor/metal with spin-orbit coupling. Although the time-reversal symmetry is explicitly broken in (2), the vector potential does not couple minimally in all the complex hopping parameters 8 , whether they preserve or flip the spin, and as such there is no Lorentz force on ψ's. By direct analogy with the anomalous Hall effect 15 in semiconductors/metals, any intrinsic κ xy is therefore driven by Berry phase mechanisms.
Clearly, for ∆ = 0, the model in Eq.(1) reduces to the well known square lattice Hofstadter "butterfly" The Chern integers ℓ r , determining the electrical conductivity σ (el.) xy = e 2 h ℓ r associated with M filled sub-bands, obey the Diophantine equation 17,18 M = sL 2 H +ℓ r , where s is an integer, and where ℓ r is subject to the restriction |ℓ r | ≤ L 2 H /2. Solving for ℓ r we find that 18 (for even L 2 H ) each of the lowest, and each of the highest, L 2 H /2 − 1 sub-bands carries the Chern number +1. In addition, there is a pair of touching bands just above and just below zero energy, which carries the (large) Chern number 2 − L 2 H , equally divided between the two anomalous sub-bands 18 . Near the band extrema the dispersion is parabolic and the energy gap between the magnetic subbands is There are four gapless Dirac points in the 1 st Brillouin zone, with the velocity anisotropy α D = v F /v ∆ = t/∆, which dominate the low temperature thermodynamics.
Thus, the pairing term in H mixes ∼ 4∆/ ω c = L 2 H /(πα D ) sub-bands near the Fermi level. In the physically relevant situation, each of these magnetic sub-bands mixed by ∆ carries a unit Chern number, thus resembling Landau levels of a continuum problem. Moreover, the number of such occupied sub-bands should be large compared to L 2 H /(πα D ). Therefore, we choose to set µ = 2t; the results presented here are not too sensitive to this choice as long as the pairing term does not mix the anomalous band with the large Chern number, or as long as we are not too close to the band minimum.
We now turn to B = 0 and ∆ = 0. In the semiclassical approximation of Volovik 1 , the spatially varying phase θ(r) induces a finite zero-energy density of states, N (E = 0), which scales as √ H. More generally, the scaling expected 4,8 from the approximate Dirac model is As shown in Fig.2, up to small variations 14 , we indeed recover this scaling in our model when L 2 H /(πα D ) ≫ 1. The integrated den-sity of states per area, per spin, per layer, which is more amenable for comparison when obtained numerically, is seen to follow the scaling

II. INTRINSIC CONTRIBUTION TO κxy
Having established the connection with known results, we now move to the main focus of this paper. It has long been known that for B = 0 the thermal transport coefficients cannot be calculated using the Kubo formula alone 10,19-22 . Rather, κ xy is given by 10,21 where f (ξ) = 1/ e ξ/kB T + 1 , and  The double sum over m and n is to be performed subject to the stated constrain. It is well known, thatσ xy (ξ) = m Ω m (ξ), where Ω m (ξ) is the integral over the Berry curvature of the m th magnetic sub-band over the parts of the magnetic Brillouin zone which are below the energy ξ. If the sub-band is fully occupied then the integral over the Berry curvature is the first Chern number, i.e., an integer 10 .
Our main result is shown in Fig. 2. The ratio κ xy (T, ∆, H)/κ xy (T, 0, H), in the regime where κ xy (T, 0, H) is linear in T and L 2 H /(πα D ) ≫ 1, is seen to collapse quite well onto a single curve for each value of the Dirac cone anisotropy α D , and approach 1 with increasing T /∆, independent of the vortex lattice geometry. The rapid onset of κ xy with increasing temperature can be understood from the ξ-dependence of σ xy (ξ) and the formula in Eq.(3). The quantityσ xy (ξ), shown in Fig. 3, has a simple physical interpretation: it is proportional to the quasiparticle contribution to the T = 0 spin Hall conductivity 10 if all the qp subbands below the energy ξ are occupied. In addition, for ∆ = 0 it can be related to the usual electrical Hall conductivity σ (el.) . We see that for ∆ = 0, the quantity ξ can be interpreted as a fictitious Zeeman energy splitting. Therefore, the loss of the total Chern number in the minority band is largely compensated by the gain in the majority band, and for ∆ = 0,σ xy (ξ) is essentially independent of ξ. For a range of ξ's near zero, its value is approximately 2ℓ r , i.e., the solution of the mentioned Diophantine equation determining the normal state electrical conductivity at µ (see 1/α D =0 curve in Fig. 3 for small ξ). This is expected, since an analogous formula to Eq.(3) holds in the normal state as well 21 , in which case it relates the electrical Hall conductivity to the thermal Hall conductivity. For ∆ = 0, theσ xy (ξ) still reaches the value ∼ 2ℓ r , i.e., of order L 2 H , when ξ ≫ ∆. That is because for ξ large compared to ∆, the contribution to the Chern numbers comes predominantly from the qp bands which are weakly affected by the pairing term. On the other hand, when ξ ≪ ∆,σ xy (ξ) becomes small, non-universal and its value oscillates around zero. This trend is displayed in Fig. 3, where the ξ dependence ofσ xy (ξ) is shown for various values of ∆. Such precipitous drop near ξ = 0 is a consequence of the qp states near zero energy being an almost equal superposition of the electron and a hole, ridding the sub-bands of the Berry curvature. It is fully consistent with the previous result 11 where the value ofσ xy (0), small compared with L 2 H , was obtained near µ = 0.
Convoluting such ξ-dependence ofσ xy (ξ) with the thermal factor in Eq.(3) at low T , will result in a vanishingly small κ xy . As the temperature increases, the thermal function broadens, and the rapid increase ofσ xy (ξ) with ξ is mirrored by the rapid increase of κ xy with T . In Fig. 3, we restrict the k B T ≪ t guaranteeing that in the normal state (∆ = 0), κ xy is linear in T . The magnetic field dependence follows readily from the above discussion as well. Since for ξ ≫ ∆,σ xy (ξ) ∼ L 2 H ∼ 1/H, we find that the H-field dependence of κ xy (∆ = 0) is entirely determined by the H-dependence of κ xy (∆ = 0) ∼ 1/H. This is shown in Fig. 3. We find the dependence of κ xy on the Zeeman coupling, h Z , to be barely observable when h Z ≪ ∆. (see Fig. 3 in the SI). This is due to the fact that, with the Zeeman term, we convolute the average ofσ xy (ξ ± h Z ) instead ofσ xy (ξ). For any h Z ≪ ∆, the main effect of the averaging is to smear the fluctuations inσ xy (ξ), while the overall shape of the curve remains unchanged. Interestingly, the analogy with the anomalous Hall effect in semiconductors suggests that, in the presence of quenched disorder either in the form of vortex lattice imperfections or impurity scattering, there may be regimes where skew scattering or side jump play a role.
The importance and the novelty of our results stem from the role that the overall d-wave pairing strength, ∆, plays in determining the temperature scale at which κ xy rapidly increases from a negligibly small value at low T . It may therefore provide a means of measuring ∆ in magnetic field via a bulk transport measurement.

Acknowledgments
We wish to thank Prof. Ong for sending us their unpublished data. We also acknowledge helpful conversations with Prof. Kun Yang, Dr. Ashot Melikyan, and support from the NSF CAREER award (

I. THE HAMILTONIAN AND ITS SYMMETRIES
We concentrate on the Hamiltonian as presented in Eq.
(2) in the main text which we restate here: Here the Franz-Tesanovic singular gauge transformation has been performed 1-3 , so the the hopping amplitudes are given as The branch cut variable z (2) r,r+δ is equal to +1 on all links except those crossed by the branch cut (green wavy line in the upper left panel of Fig. 1) where it is equal to −1. The pairing strength ∆x = +∆, ∆ŷ = −∆ has the d x 2 −y 2 -wave symmetry. Per construction, z (2) r,r+δ is periodic under the magnetic unit cell translation. The phase θ(r) is a solution to the London equation, ∇ × ∇θ(r) = 2πẑ j δ(r − r j ). With the periodic arrangement of vortices, as shown in Fig. 1, the solution is Here z = x + iy, and the summation over j accounts only for the two vortex positions z j = x j + iy j . σ(z; ω, ω ′ ) is the Weierstrass σ function with periods ω = L x and ω ′ = iL y . Constants γ and v 0 are determined by the boundary conditions. In order to ensure that the superfluid velocity v ∼ 1 2 ∇θ − e c A satisfies periodic boundary conditions, where ζ(. . . ; ω, ω ′ ) is the Weierstrass ζ function. For square lattices γ = 0. Setting v 0 = π/L 2 H , makes the overall current in the system vanish. The former boundary condition ensures that the hopping amplitudes t in Eq. (1) are periodic.
The Eq. (1) can be rewritten as i L Let us relabel ψ r,σ (k) = χ L−r,σ (−k) in Eq. (1) so that we can write it as As implied by the notation, we introduced k ′ = −k and r ′ = L − r − δ. The term in the first line of Eq. (9) comes from the "h.c." term of the first line in Eq. (1), and the similar is true for the third line. From Eq. (8) we find that the couplings in Eq. (9) obey The Hamiltonian Eq. (9) is therefore describing the same system as Eq. (1) apart from the branch cut which has been moved by the inversion operation. This Hamiltonian therefore corresponds to the situation in the upper right panel on Fig. 1. The branch cut may be restored to its original position by a gauge transformation, γ , which changes sign across the border of the rectangular area delineated by the original and the new branch cut and shown in gray in Fig.  1, r ∈ means that site r is inside the rectangular area. After the gauge transformation Eq. (9) becomes This procedure therefore defines a unitary transformation corresponding to {i|L} followed by gauge transformation γ , whose matrix elements are The Kronecker delta δ r,L−r ′ is defined modulo magnetic unit cell translations. Eq. (13) is then written as Comparing the second line with Eq. (7) we prove thatĤ(−k) = U † γ {i|L} Ĥ (k)U γ {i|L} which was our aim. Another operation preserving the vortex arrangement is vertical mirror perpendicular to the y-axis and passing through the coordinate center, {σ y |0}. When applied to Eq. (5) this operation implies The negative sign in the exponent on the right hand side reflects the fact that a vertical mirror reverses the sign of the magnetic field, therefore we should follow {σ y |0} by the time reversal (TR) operation Θ which acts as a conjugation. We then relabel ψ r,σ (k) = χ σ y r,σ (−σ y k) * = χ σ y r,σ (σ x k) * , and repeat the procedure done for the inversion operation. The result is analogous to Eq. (9): the Hamiltonian for χ(k) fields corresponds to the system with branch cuts repositioned by {σ y |0} operation, as shown in the lower left panel of Fig. 1. A gauge transformation γ x intended to restore the branch cut to its original position must change sign on the vertical line defined by the old and new portions of the branch cut. For example, it could multiply χ r,σ (k) by −1 on the sites left of that line and by +1 on the sites right of the line. However, such a gauge transformation is not periodic unless followed by a phase prefactor: The symmetries of a triangular VL combined with optional gauge transformations (γ) and/or TR (Θ) which relate the single particle Bloch Hamiltonians at different k points. The matrix elements of the corresponding (anti-)unitary transformations are listed in the next column, K representing the complex conjugation. The next column gives k ′ related to k through the operation, and the last column the color corresponding to the area of the magnetic BZ in Fig. 2.
Here, x = r ·x and 0 ≤ x < L x . The anti-unitary transformation following from this set of operation is therefore defined by its matrix elements and it establishes thatĤ Following the same procedure, we show that vertical mirror {σ x |0} followed by a TR and e i πŷ Lx ·r maps the original problem to the one with the branch cuts shown in the lower right panel of Fig. 1 In the case of a VL which is not a square lattice, all other symmetries, eight in total, are generated from the three we have presented. These are enumerated in Table I. When we introduce the methods for the calculation of the spin Hall conductivity,σ xy (ξ), in the next section, we use these symmetries to prove that the contributions toσ xy (ξ) from two different symmetry related k's are equal. This ultimately speeds up the calculation by reducing the integration domain from the full magnetic BZ to only a small representative fraction thereof.
The special case of the square lattice, where L x = L y , possesses additional symmetries. A vertical-diagonal mirror which contains A and B vortices, followed by a TR, Θ{σ Y |0}, preserves the VL arrangement. Since such a transformation changes the position of the branch cut, it must be followed by γ in order to restore the branch cuts. This is not sufficient to recover the original Hamiltonian, due to the nature of the d x 2 −y 2 -wave gap. This operation changes the sign of the off-diagonal hopping terms, λ's, and must be combined with a charge gauge transformation which restores the sign on λ's, γ λ = e iπσ 3 = iσ 3 . Only then we can define an anti-unitary operation, such that The combination of this operation with the previously defined eight operations yields another eight operation and the total of sixteen point group operations for the case of a square VL.
In the absence of the Zeeman term, h Z = 0, another symmetry, which is independent of the VL arrangement, relates the single particle Hamiltonian at k and −k: As a result, for each eigenstate |m; k with energy E m (k), there is another state |m ′ ; −k = σ 2 |m; k with energy The magnetic Brillouin zone (BZ). The single particle Hamiltonians at k's of different colors are related by symmetries tabulated in Table I. For example, k in white and k ′ = −k in gray area have their respective single particle Hamiltonians related by γ {i|L} symmetry.

II. SPIN-HALL CONDUCTIVITY
Thermal Hall conductivity at finite temperature is obtained by convolution of spin Hall conductivity,σ s µν (ξ), given by Eq. (3) in the main text which we repeat here: Here, f (ξ) = 1/(e ξ/(kB T ) + 1). The key ingredient in our work is therefore the spin Hall conductivity and here we present methods to calculate it and discuss their efficiency.
In the main text we have the following formula for the transverse spin Hall conductivity (also Ref. 3), The velocity operator isV and V mn µ (k) = m; k|V µ (k)|n; k . Notice that the transverse spin Hall conductivity is a topological property of the system:σ s xy (ξ = 0) is proportional to the sum of first Chern number over the occupied bands 3 . The Eq. (26) can be understood asσ xy (ξ) = 1 being the Berry curvature summed over all states that are below ξ at k. We will present several formulas for calculating F k (ξ) and prove their equivalence. In the next subsection we will explain the best strategies to calculateσ xy (ξ) based on these formulae. The double summation over all occupied/unoccupied states in Eq. (28) requires iterating over both m and n which, when naively implemented, is inefficient and does not have a clear physical connection to the topological properties of individual bands. This can be addressed by transforming this equation into a sum over individual states with eigenvalues below ξ, where, for brevity, we write ∂m k ∂kµ = ∂ ∂kµ |m; k . When defining the derivatives of the eigenstates in the previous equation, we have to be mindful about the freedom to choose the overall phase of each individual eigenstate. Transformation |m; k → |m; k e iαm(k) , where α m (k) is an arbitrary, but smooth, function for each band m, leaves Eq. (29) invariant: terms proportional to ∂ x α m (k)∂ y α m (k) arise, but mutually cancel. In practice, however, the appearance of such terms leads to numerical instability of Eq. (29). In order to evaluate the wave-function derivatives, one needs to diagonalize the Hamiltonian at k as well as in some vicinity of that point. This leads to large terms appearing in this expression. While these are supposed to mutually cancel, they introduce numerical errors comparable to the value of the integrand.
One way to circumvent the problem is to use projectors in place of the wave functions 4 . Projectors, P (l) (k) = |l; k l; k|, do not carry the ambiguous phase factor, which is why the previously discussed problem does not occur during the differentiation. Here we derive expressions equivalent to Eq. (29). The first step is to evaluate V mn µ . Differentiating single particle Hamiltonian we obtain From here it follows that assuming that we are interested in terms with m = n where the first term from Eq. (30) vanishes. We use the following shorthand notation, ∂P (l) /∂k µ = P (l) ,µ . We also drop k in the notation from now on. A useful relation at this step, and later in the derivation, follows from differentiating P (l) P (l ′ ) = δ ll ′ P (l) identity: Setting l = l ′ , we write Eq. (31) as eliminating the summation over l. Differentiating the identity resolution, we obtain another useful identity l P (l) ,µ = 0, which allows us to write the second term in Eq. (33) as and The Berry curvature now becomes The summation over states with E n > ξ can be substituted by a summation over states below ξ using the identity resolution, This yields Let us first show that the second term is zero. Notice that this term can be cast as Terms in the sum with m = m ′ vanish due to the cyclical property of the trace. We are left only with terms where m = m ′ . For those, we use the properties of projector operators and make the following transformations ,y Notice that the first and the last line are equal up to a minus sign, therefore this expression must be equal to zero. We are left with the first term in Eq. (39), that is, Each term in the sum is the Berry curvature of m'th band at point k in the BZ. An equivalent formula can be found in Ref. 4. While Eq. (42) does not suffer from the numerical instability the same way Eq. (29), we found that it may become unstable at, and in the vicinity of, k's where two bands are degenerate. When this is the case there is an SU(2) freedom in choosing the respective projectors for the two degenerate bands at k. Just like before, this freedom results in arbitrarily large terms that are supposed to cancel, but instead introduce an error of the order of F k (ξ). We found an alternative, but related, numerically more stable expression to remedy this problem. Instead of calculating Berry curvature of individual bands, here we define projector onto all states below energy ξ, Utilizing the identities for projector operators show above, it can be shown that Eq. (42) is equivalent to This expression does not suffer from any of the previously discussed numerical instabilities.

A. Strategies for computing F k (ξ)
Calculating the spin Hall conductivity, translates into doing a Riemann sum over the Brillouin zone (BZ) with F k (ξ) as the integrand. The method we developed for the Riemann summation is presented in the next section. Here we concentrate on the methods for calculating F k (ξ) and analyze their computational complexity. The calculation of F k (ξ) requires computationally costly operations, diagonalization of large matrices (dimension 2L 2 H ), construction of projectors, summation over a large number of states, etc. Here we outline our algorithm for performing these operations in the most efficient way, using either Eq. (26), Eq. (42), or Eq. (44) The convolution in Eq. (25) requires the knowledge ofσ xy (ξ) at arbitrary values of ξ; therefore it is desirable that we are able to return the value of F k (ξ) quickly for any ξ. An important observation is that F k (ξ) is a constant as a function of ξ unless ξ coincides with one of the Hamiltonian eigenvalues E l ; there F k (ξ) exhibits a jump. This is true regardless of which definition of F k (ξ) we use. We rely on this property to create an ordered lookup table where we store pairs of values (E l , F k (E l + 0 + )), i.e., the list of eigenvalues and F k at ξ's just above each eigenvalue. Once the lookup table has been created, finding F k (ξ) requires a search for the highest E l < ξ in the lookup table and reading off the corresponding F (when ξ is below the bottom of the band F k (ξ) = 0). Given that the table is ordered, the complexity of such an operation is O(log L H ).
What about the computational complexity of the lookup table generation? The diagonalization of the Hamiltonian at k is the first step regardless the method of finding F k (ξ). Its computational cost is O(L 6 H ) as we need the entire spectrum as well as all the eigenstates. In the case when we are supposed to find the derivatives of the projectors, we need to perform several diagonalizations in the vicinity of the k point and find the projectors there to. This adds an overall constant factor to the computational cost, but does not change its asymptotic behaviour.
For each eigenvalue E l , building the corresponding projector P (E l + 0 + ) requires computation of each element of the projector matrix, which is a dense matrix, thus O(L 4 H ) operations. While it may seem that for E l 's far from both the bottom and the top of the band one needs to repeat this l ∼ O(L 2 H ) times, once for each band below E l , this can be avoided by deriving P (E l + 0 + ) from already constructed P (E l−1 + 0 + ), Calculating the derivatives, as mentioned before, requires this operation to be repeated several times in the vicinity of k. While adding a constant prefactor, these operations do not change the computational complexity of this step asymptotically. The most computationally expensive step comes from evaluating the product of matrices in Eq. (44). Naive matrix product operation has a computational complexity proportional to the cube of the matrix size, in our case O(L 6 H ). Having to repeat this for each entry in the lookup table, we conclude that the construction of the lookup table for F k (ξ) has the computational complexity which scales as O(L 8 H ). Methods for matrix multiplication that scale with lower power of the matrix size 5,6 are known. Overlooking the fact that these algorithms offer only a negligible ( 15%) speedup over the standard library methods 7 for typical matrix sizes we use, we note that even the most sophisticated matrix-product algorithm presently known 8 cannot reduce the asymptotic computational complexity of the lookup table creating below O(L 6.745 H ). The theoretical limit on the matrix multiplication complexity 9 sets the lower bound for the lookup table creation with projectors to O(L 6 H log L H ). We conclude that this is the lower bound for the computational complexity of the lookup tables for F k (ξ) based on either Eq. (42) and Eq. (44). In practice, however, we found that the execution time of the lookup table computation scaled as O(L 8−ǫ H ) with ǫ < 1. Here we show that we can do much better using the original, double sum, formula, Eq. (28). Naively, the integrand F k (ξ) in that equation contains summation over two indices, m and n, and both sums have O(L 2 H ) terms unless ξ is close to the bottom or the top of the band. Further, finding each V mn µ for a general V µ matrix would require summing O(L 4 H ) terms. In total, these estimates imply that the creation of a lookup table has O(L 10 H ) computational complexity. We reduce the complexity to O(L 6 H ) using the following observations. First notice that the difference between two adjacent F 's in the lookup table is given by The most of the terms in the difference cancel and we are left with only 2L 2 H − 1 terms to sum over. Therefore, if we build the lookup table by adding Eq. (46) to F k (E l−1 + 0 + ) to obtain F k (E l + 0 + ), for each lookup Next, we observe that V µ are sparse matrices, each having only O(L 2 H ) non-zero elements. This follows from the fact that they derive from the single particle HamiltonianĤ(k) which is itself a sparse matrix. The computation of V mn µ therefore requires a summation of only O(L 2 H ) terms, thereby improving the computational time by another O(L 2 H ) factor. The computation of V mn µ 's can further be optimized by noticing that, once l is fixed, we can find V µ |l , store those two vectors, and use them to compute V ml µ and V ln µ faster. Overall, we find that the computational complexity of the lookup table creation is O(L 6 H ) if Eq. (26) is used in conjunction with these improvements. A further advantage of this method over those using the projectors is that, for each F k , the diagonalization of the Hamiltonian has to be performed only once, at k point.

B. Symmetries and F k (ξ)
Now we show that, if two Bloch Hamiltonians,Ĥ(k) andĤ(k ′ ) are related by an unitary or anti-unitary transformation, then F k (ξ) = F k ′ (ξ) for all ξ. This property allows us to calculate the contributions toσ xy (ξ) coming only from a fraction of the magnetic BZ which contains symmetry inequivalent k points only; the contributions from the remainder of the integration domain is then proportional to what we have calculated. In Fig. 2 such an area is represented by a white triangle for the case of a non-square VL. Since this area is one eight of the original magnetic BZ, the use of symmetries speeds up our calculation by a factor of 8. When the vortices form a square lattice, the integration domain is further reduced by a factor of two and the speed of the calculation is correspondingly improved.
Let us consider, for illustration, γ {i|L} operation in Table I for which k ′ = −k. We have |m; −k = U † γ {i|L} |m; k = U † |m; k with U listed in the table and its argument suppressed. We also have,Ĥ(−k) = U †Ĥ (k)U . Since U is k-independent we havê Here we used the fact that the derivative w.r.t. q of an even function in q is an odd function. Combining these identities we find that Eq. (28), evaluated at −k, becomes We used E l (−k) = E l (k) and dropped the k-argument accordingly. The proof is similar for other unitary symmetry operations.
To illustrate how to equate F k (ξ) for k and k ′ related by an anti-unitary symmetry, let us consider Θ{σ x |0} in Table  I. For this operation k ′ = πŷ Ly + σ y k and we use k ′ for brevity. Since this operation is anti-unitary, |m; k ′ = U * |m; k * . As before,Ĥ(k ′ ) = U †Ĥ (k) * U which impliesV x (k ′ ) = −U †V x (k) * U andV y (k ′ ) = U †V y (k) * U . Combining these we find We have just shown that F πŷ Ly +σ y k (ξ) = F k (ξ). The proof is analogous for all other anti-unitary operations. From here it follows that F k ′ (ξ) = F k (ξ) for eight k ′ corresponding to each operation in Table I. In Fig. 2 we colour different sections of the magnetic BZ according to these symmetries. It is sufficient to calculate the Berry curvature in one triangle, say the white one, to be able to find the Berry curvature in the entire magnetic BZ by applying these symmetries. Hence, in order to find the spin Hall conductivity we only need to concentrate on the white triangle and multiply the result we obtain by 8.
In the case of the square VL, the diagonal mirrors and the symmetries generated through these fold the area of the integration further by a factor of two, since, in this case we can prove that F (kx,ky) (ξ) = F (ky,kx) (ξ). When this is the case, the integration domain is one half of the white triangle, and this domain is defined by corners at (0, 0), (π/2ℓ x , π/2ℓ y ), and (π/2ℓ x , −π/2ℓ y ).
The symmetry of the Hamiltonian expressed in Eq. (24), since unitary, implies that where E l 's are the eigenvalues at k and, thanks to this symmetry, −E l 's are the eigenvalues at −k. Given that F k (ξ) is a constant unless ξ coincides with an eigenvalue at k, the right hand side of Eq. (50), we can rewrite that as The expression on the left hand side contributes toσ xy (ξ = E l ), whereas the one on the right hand side contributes toσ xy (ξ = −E l ) up to the minus sign. Therefore, we conclude thatσ xy (ξ) is an even function of ξ, When the VL has the inversion symmetry, as it is the case for us, this symmetry also implies that F k (ξ) is an even function of ξ. We therefore calculate F k (ξ) only for negative ξ's and from there infer the values for positive ξ, thereby reducing the computation time by another factor of two.

C. Strategies for computingσxy(ξ)
Recalling the equation for the spin Hall conductivity from the previous section,σ xy (ξ) = 1 2π d 2 k F k (ξ), we see thatσ xy (ξ) can be approximated and calculated by expressing it as a Riemann sum, where the integration domain (in our case a fraction of the magnetic BZ defined previously) is partitioned into s areas, each with area A i . For each area, the Berry curvature is calculated at some point k i within that area. The partitioning of the integration domain is performed using adaptive mesh integration method. The details of the method will be presented elsewhere, it suffices to know that this method starts with a coarse mesh of k-points and then refines the mesh, i.e., increases the density of k-points in the areas where the absolute errors are the largest. In this manner the overall error is minimized using the fewest number, s, of k-points where F k (ξ) is calculated. The calculation of F k (ξ) at s points requires O(sL 6 H ) operations using the method presented previously and this is the most time consuming part of our calculation. From here, we could store the lookup tables for F at each k i and invoke Eq. (53) whenever we needσ xy (ξ). Each such call would require O(s log L H ) time. For us 10 3 s 3 ·10 3 , and this method might slow down the calculation of the spin-Hall conductivity significantly. Therefore, we constructed another lookup table holding values of ξ and the correspondingσ xy (ξ).
We rely on the fact that all the terms in Eq. (53) are constant as a function of ξ, except for ξ'w which coincide with any of the eigenvalues E m (k i ); there the sum exhibits a discontinuous jump originating only from F ki (ξ). The lookup table therefore is an ordered list of pairs (ξ,σ xy (ξ + 0 + )). For each eigenvalue E m (k i ) there is a pair with ξ = E m (k i ) andσ xy (ξ + 0 + ) is the value of the spin-Hall conductivity at ξ infinitesimally above ξ. Such a lookup table contains 2sL 2 H . Since it is ordered, findingσ xy (ξ) for any given ξ takes O(log s) + 2O(log L H ) time which is much faster than the naive implementation mentioned before.
An naive way to build the lookup table is to get all the eigenvalues E m (k i ), and findσ xy (E m (k i ) + 0 + ) according to Eq. (53). Since there are 2sL 2 H distinct eigenvalues, such an approach has computational complexity O(s 2 L 2 H log L H ). While this calculation step would not be the most time consuming, we nevertheless found a much more efficient algorithm to generate the lookup table forσ xy (ξ). We first notice that, for any two consecutive eigenvalues, E m (k j ) ≤ E m ′ (k j ′ ) (here k j and k j ′ may or may not be the same), the difference of the spin-Hall conductivities corresponding to these is given bỹ Since there are no other eigenvalues between E m (k j ) and E m ′ (k j ′ ), all F xy 's except F k j ′ (ξ) are constant and therefore cancel. The change between two consecutiveσ xy (ξ)'s in the lookup table is therefore calculated by updating the contribution from k j ′ only. This requires only two operations. We further optimize the processing of all the eigenvalues in the correct, sorted, order by utilizing priority queue, a data structure in which both inserting (pushing) new elements and returning the smallest one (popping) take time proportional to the logarithm of the number of elements in the structure. We do that by initially building the queue by pushing the smallest eigenvalue at each k i . In each iteration the smallest eigenvalue, E m (k i ), in the queue is popped and replaced by the next eigenvalue at the same k i (until the eigenvalues are exhausted). Since the eigenvalues at each k i are sorted, we are guaranteed to pop all the eigenvalues in the ascending order. Whenever an eigenvalue is popped, we find the spin-Hall conductivityσ xy (E m (k i ) + 0 + ) according to Eq. (54) and add a new entry to the lookup table. Since each iteration has computational complexity O(log s), the total computational complexity for building the entire lookup table forσ xy (ξ) using this algorithm is O(sL 2 H log s).

III. κxy IN THE PRESENCE OF ZEEMAN EFFECT
In Fig. 2 in the main text we present the thermal Hall conductivity with no Zeeman term present. The justification for this lies in the fact that when the result is plotted with the Zeeman term h Z < ∆ we found no observable change in κ xy . We illustrate this in Fig. 3 where we plot the particular curve of Fig. 2, the one onto which other data collapse (α D = 7.1), with the Zeeman term included. The Zeeman term strength is set to h Z /∆ = 0, 0.071, 0.143, and 0.357.