Weyl-Majorana solenoid

A Weyl semimetal wire with an axial magnetization has metallic surface states (Fermi arcs) winding along its perimeter, connecting bulk Weyl cones of opposite topological charge (Berry curvature). We investigate what happens to this"Weyl solenoid"if the wire is covered with a superconductor, by determining the dispersion relation of the surface modes propagating along the wire. Coupling to the superconductor breaks up the Fermi arc into a pair of Majorana modes, separated by an energy gap. Upon variation of the coupling strength along the wire there is a gap inversion that traps the Majorana fermions.


I. INTRODUCTION
A three-dimensional Weyl semimetal has topological features that are lacking in its two-dimensional counterpart, graphene [1][2][3]. One striking feature is the appearance of surface states, in Fermi arcs connecting Weyl cones of opposite topological charge (Chern number or Berry curvature) [4]. Unlike the surface states of a topological insulator, which are the only source of metallic conduction, the Fermi arcs at the surface compete with the Weyl cones in the bulk when it comes to transport properties. Quantum oscillations in the magnetoresistance are one example of an effect where the Fermi arcs play a prominent role [5,6], the chiral magnetic effect without Landau levels is another example [7].
An interesting way to differentiate surface from bulk is to bring the Weyl semimetal into contact with a superconductor. While the Weyl cones in the bulk remain largely unaffected, the surface states acquire the mixed electron-hole character of a charge-neutral Bogoliubov quasiparticle -a Majorana fermion [8][9][10][11][12][13]. Here we investigate this proximity effect in the nanowire geometry of Fig. 1, in which an axial magnetization causes the surface modes to spiral along the wire, essentially forming a solenoid on the nanoscale [7]. We study the dispersion relation of the Majorana modes and identify a mechanism to trap the quasiparticles at a specified location along the wire.
In the next section we identify the pair of Z 2 quantum numbers ν, κ that label the four surface modes in a given orbital subband. The electron-hole index ν is generic for any surface state where electrons and holes are coupled by Andreev reflection [14][15][16]. The connectivity index κ is specific for the Fermi arcs, it distinguishes whether the surface state reconnects in the bulk with the Weyl cone at positive or negative energy. In Sec. III we construct the 4 × 4 matrix Hamiltonian in the ν, κ basis, constrained by particle-hole symmetry, as an effective low-energy description of the two-dimensional surface modes.
We then proceed in Sec. IV with a numerical calculation of the three-dimensional band structure of a mi- croscopic model Hamiltonian. The unexpected feature revealed by this simulation is a gap inversion, visible in the band structure as a level crossing between two surface modes with the same connectivity index. The gap inversion can be controlled by variation of the tunnel coupling between the semimetal and the superconductor. At the domain wall where the gap changes sign, a chargeneutral quasiparticle is trapped -as we demonstrate numerically and explain within the context of the effective surface Hamiltonian in Sec. V. In Sec. VI we study the same gap inversion analytically, via a mode-matching calculation. In the concluding Sec. VII we comment on the relation of the gap inversion to the flow of Berry curvature in the Brillouin zone.  Fig. 1b, calculated from the tight-binding model described in the text [17]. In panel a) there is only the Weyl semimetal, in panel b) the superconducting contacts have been added. Inversion symmetry has not been broken, so the spectrum has ±pz symmetry, in addition to the particlehole symmetry E(pz) = −E(−pz). In the slab geometry the transverse wave vector ky is a good quantum number, and to make the figure less crowded only subbands at a single value of ky are shown. (The Fermi arcs in panel a) are approximately at ±v φ sin ky.) The superconductor breaks up the two Dirac fermion surface modes in panel a) into four Majorana fermion modes in panel b), labeled by a pair of indices κ, ν = ±1. The Majorana modes are nearly charge-neutral, as indicated by the color scale (with electron charge +e).

II. CONNECTIVITY INDEX OF SURFACE FERMI ARCS
The geometry under consideration is shown in Fig.  1. A Weyl semimetal wire oriented along the z-axis is covered by a superconductor. We include a thin insulating layer between the superconductor and the Weyl semimetal, forming a tunnel barrier. A magnetization in the z-direction breaks time-reversal symmetry and separates the Weyl cones along the p z momentum direction in the Brillouin zone. (Induced superconductivity in the presence of time-reversal symmetry, with minimally four Weyl points, has a different phenomenology [13].) The surface states connecting the Weyl cones are chiral, circulating with velocity v φ in a direction set by the magnetization. If inversion symmetry is broken the surface states also spiral with velocity v z along the wire [7].
As shown in Fig. 2, resulting from a model calculation described in Sec. IV, at the interface with a superconductor the surface spectrum is drastically modified. We seek an effective Hamiltonian that describes this proximity effect on the Fermi arcs.
The first question we have to address is which pairs of states are coupled by the superconducting pair potential ∆. In the bulk spectrum the answer is well known [8,12]: Superconductivity couples electrons in a Weyl cone of positive Berry curvature to holes in a Weyl cone of negative Berry curvature, and vice versa. To decide this question for the surface states, we assign to each Fermi arc a "connectivity index" κ = ±1, depending on whether it reconnects in the bulk with the Weyl cone at positive or negative energy. Inspection of Fig. 2 shows that ∆ predominantly couples Fermi arcs with same κ, pushing them apart, without removing the crossing between states of opposite κ.
More explicitly, in a slab geometry we can identify κ = sign k y and in a cylindrical wire geometry we would have κ = sign p φ . The coupling of states with different κ is then forbidden by (translational or rotational) symmetry. More generally, in the absence of any symmetry, the sign of κ = ±1 says whether the Fermi arc connects with the Weyl cone at ±E, and thus identifies which pairs of Fermi arcs are predominantly coupled by ∆.

III. EFFECTIVE SURFACE HAMILTONIAN
The superconducting proximity effect is governed by the Bogoliubov-De Gennes (BdG) Hamiltonian, describing the coupling of electrons and holes by the pair potential. In the numerical simulations we will work with the BdG Hamiltonian in a 3D microscopic model. For analytical insight we aim for an effective 2D description involving only surface modes.
Each orbital subband n is associated with four Majorana modes, labeled by a pair of Z 2 indices κ, ν. (See Fig.  2.) The connectivity index κ = ± identifies the connectivity of the surface mode (with the Weyl cone at positive or negative energy), the electron-hole index ν = ± identifies the pair of Majorana fermions that form a Dirac fermion. The corresponding BdG Hamiltonian H n is a 4×4 matrix with p z -dependent elements. In what follows we omit the subband index n for ease of notation.
The fundamental symmetry of the BdG Hamiltonian is particle-hole symmetry, with Pauli matrices κ α and ν α acting, respectively on the connectivity and electron-hole degree of freedom (α = 1, 2, 3 → x, y, z and α = 0 for the unit matrix). The operation of particle-hole conjugation squares to +1, which places the system in symmetry class D [19] -this is the appropriate symmetry class in the absence of timereversal and spin-rotation symmetry.
If we neglect the mixing by disorder of surface states with opposite connectivity index κ = ±, the 4 × 4 matrix H decouples into two blocks H ± related by particle-hole symmetry, The 2 × 2 matrices H ± can be decomposed into Pauli matrices, with real p z -dependent coefficients D α . Diagonalization of the Hamiltonian (3.2) gives the dispersion relation E κ,ν (p z ) of the four Majorana modes in the n-th subband, is an even function of p z while each of the functions D 1 , D 2 , D 3 has a definite parity (even or odd).

IV. NUMERICAL SIMULATION OF A MICROSCOPIC MODEL
We now turn to a microscopic model of a Weyl semimetal in contact with a superconductor, which we solve numerically. The Weyl semimetal has BdG Hamiltonian with chemical potential µ and charge operator The Pauli matrices σ α and τ α refer to spin and orbital degrees of freedom, respectively, while ν α acts on the electron-hole index. The momentum k varies over the Brillouin zone |k α | < π of a simple cubic lattice (lattice constant a 0 ≡ 1). This is a model of a layered material in the Bi 2 Se 3 family [20], with weak coupling t z < t in the z-direction, perpendicular to the layers in the x-y plane. The particle-hole symmetry relation is The magnetization term ∝ β breaks time-reversal symmetry, H W (k) = σ y H * W (−k)σ y . Inversion symmetry, The Weyl semimetal is in contact with a spin-singlet s-wave superconductor, with Hamiltonian There are different chemical potentials in the Weyl semimetal, µ = µ W , and in the superconductor, µ = µ S . At the NS interface we include an electrostatic potential barrier of width d barrier , raising µ to a value µ B ≡ U barrier . The resulting spatial profile µ(x) is shown in Fig. 3. We consider the two geometries shown in Fig. 1, a wire geometry and a computationally more efficient slab geometry [21]. In each case there is translational invariance along the z-direction. In the slab geometry there is in addition translational invariance in the y-direction, so the modes are labeled by a continuous quantum number k y [22].
The dispersion relation in the slab geometry is shown in Fig. 2. The mode crossings at nonzero p z appear because modes with different connectivity index κ are uncoupled in the absence of disorder. In Fig. 4 we show a different type of crossing, near p z = 0 between modes with the same κ, induced by variation of the tunnel barrier height. This crossing appears generically when we vary interface parameters, in Fig. 5 we show that it persists at nonzero chemical potential µ = µ W in the Weyl semimetal [23]. Inversion symmetry breaking by a nonzero λ moves the crossing point away from p z = 0, but does not destroy it. The wire geometry gives similar results, see Fig. 6.
To model this effect in the framework of the surface Hamiltonian (3.3), we take a momentum-independent complex off-diagonal potential D 1 − iD 2 ≡ ∆ with amplitude ∆ 0 = c(U barrier − U c ) that crosses zero at some critical barrier height U c . Inversion symmetry imposes a definite parity on the real diagonal potential D 3 ≡ µ(p z ), such that even a small admixture of an odd-parity component enforces µ(0) = 0 when λ = 0. If we take µ(p z ) = c λ + c p z the dispersion relation (3.4) in the pair of modes with κ = +1 has the form (4.5) The dashed curves in Fig. 4 are fits to this functional form, with λ = 0 and a quartic D 0 (p z ). The qualitative behavior agrees reasonably well.

V. QUASIPARTICLE TRAPPING BY GAP INVERSION
The gap inversion of Fig. 4 can be used to trap a quasiparticle by varying the tunnel barrier height U barrier (z) (by means of a variation in the thickness of the insulating layer), from a value above the critical strength U c to a value below U c . A demonstration of this effect in the slab geometry is shown in Fig. 7, where we  [17], except for the tunnel barrier height U barrier , which is varied to tune through the gap inversion. The dashed curves are fits [18] to the dispersion (4.5) from the effective surface Hamiltonian.
In terms of the surface Hamiltonian, the quasiparticle trapping is described by the Schrödinger equation (5.1) We take a real ∆(z) = c(U barrier (z) − U c ) and, respectively, an even and odd p z -dependence of D 0 and µ = c p z -consistent with inversion symmetry. If we neglect quadratic terms in D 0 we have a matrix differential equation of first order, Let ∆(z)/c vary from a positive value for z < 0 and z > L to a negative value in the interval 0 < z < L. For sufficiently large L we can consider the domain wall at z = 0 separately from the one at z = L. At energy E = ±D 0 (0) there is a bound state at z = 0 with wave function This should be a decaying function of |z|, so ψ ± (0) = (1, ±i) is an eigenstate of ν y with eigenvalue ±1. Fig. 7 shows that the bound state is a charge-neutral quasiparticle. There is one state at energy +D 0 (0) and a second state at −D 0 (0), but because the BdG equation doubles the spectrum only a single Majorana fermion is trapped at z = 0. A second Majorana fermion is trapped at z = L. All of this is for a single orbital mode n. We have found numerically that the critical barrier height U c is weakly n-dependent, so a domain wall traps one Majorana fermion per orbital subband.  Fig. 2b [17], but we took a nonzero µW = 0.05 t0 (notice the displacement of electron and hole bands in the bulk Weyl cones) in order to demonstrate that the level crossing does not require a vanishing chemical potential. The level crossing also persists if inversion symmetry is broken by a nonzero λ = 0.05 t0, but the crossing point is displaced away from pz = 0 (compare black and red curves in panel b, at pz = 0 and pz = −6 · 10 −4 /a0).

A. Hamiltonian with spatially dependent coefficients
To analytically substantiate our numerical findings we have performed a mode-matching calculation in the slab geometry of Fig. 1b, matching electron and hole modes in the normal (N) region 0 < x < W to Bogoliubov quasiparticles in the superconducting (S) regions x < 0, x > W . This procedure can be greatly simplified if we choose a single BdG Hamiltonian H with x-dependent coefficients, rather than the different H W and H S of Sec. IV -the former choice is a less realistic model of an SNS junction than the latter, but as we will see the results are essentially equivalent.
Our starting point is therefore the Hamiltonian We will compare our analytical mode-matching calculation to a numerical solution of the discretized Hamiltonian (6.1). For this analytics, but not for the numerics, we make one further simplification, which is to linearize the Hamiltonian in the transverse momentum component k x , so that the mode-matching calculation requires the solution of a set of first order differential equation in x. We thus replace sin k x → k x and replace the mass term (6.2) bỹ m(k y , k z ) = m 0 + t(1 − cos k y ) + t z (1 − cos k z ). (6.3)

B. First-order decoupling of the mode-matching equations
The Schrödinger equation Hψ = Eψ produces 8 coupled differential equations, and an attempt at direct solution produces unwieldy results. Our approach is to partially decouple these by suitable unitary transformations of H. We take the inversion symmetry breaking strength λ and chemical potential µ as small parameters and seek a decoupling up to corrections of first or second order in λ, µ.
For a first-order decoupling we rotate the ν x and τ x spinors by the unitaries The rotation angles θ, φ are x and k z -dependent, cos θ = −(t z /∆ eff ) sin k z , sin θ = ∆/∆ eff , (6.5a) cos φ = ∆ eff /M, sin φ =m/M, (6.5b) ∆ eff (x) = ∆ 2 (x) + t 2 z sin 2 k z , (6.5c) Notice that cos θ → −sign k z for ∆ → 0. We can avoid this discontinuity at k z = 0 by keeping a small nonzero ∆ in the normal region. The transformed Hamiltonian, − µ cos θ ν z τ 0 σ 0 − µ sin θ cos φ ν x τ z σ z − µ sin θ sin φ ν x τ x σ 0 + λ sin θ ν x τ 0 σ z + λ cos θ cos φ ν z τ z σ 0 + λ cos θ sin φ ν z τ x σ z , (6.6) is diagonal in the ν and τ degrees of freedom up to corrections of first order in λ, µ, and up to a boundary potential V b (x) resulting from the commutator of k x = −i∂/∂x and the x-dependent superconducting gap ∆(x) at the NS interface. In this section we discard the boundary potential, to simplify the calculations -we will fully include it in the Appendix. The term ∝ µν x τ x σ 0 in the Hamiltonian (6.6) can be made diagonal in ν and τ with the unitary transformation The four blocks in the shift matrix P 3 [with (P 3 ) 3 = 1] refer to the ν degree of freedom. The transformed Hamiltonian is The symbol δ keeps track of the order in λ, µ of the diagonal ("diag") and off-diagonal ("offdiag") blocks.

C. Second-order decoupling via Schrieffer-Wolff transformation
The Schrieffer-Wolff transformation with Hermitian off-diagonal matrix δS given by [δS, H diag ] = iδH offdiag , (6.10) removes the off-diagonal blocks up to corrections of second order in δ: The solution of Eq. (6.10) is [25] The Schrieffer-Wolff matrix δS contributes terms of order δ 2 to the energy spectrum, which is given by the eigenvalues of H diag + δH diag + δH SW with

D. Dispersion relation of the surface modes
The mode-matching calculation at energy E with the Hamiltonian H diag + δH diag (not yet including the Schrieffer-Wolff correction) now involves four uncoupled differential equations, labeled by ν, τ ∈ {−1, +1}, for a two-component spinor ψ(x): We solve this for piecewise constant coefficients. For the normal (N) region at 0 < x < W we choose ∆ = ∆ N , µ = µ N , (6.15a) and for the superconducting (S) region at x < 0 and x > W we choose ∆ = ∆ S , µ = µ S , (6.15b) demanding continuity of ψ(x) at x = 0, W . We keep a finite pair potential ∆ N in the normal region to avoid the discontinuity at p z = 0 noted in Sec. VI B.
To obtain the dispersion relation at a single NS interface we may take W → ∞ and match decaying wave functions at both sides of the interface at x = 0. Such a bound surface state is possible if M ν − β has the opposite sign in N and S, which requires ν = +1 (since β and M are both positive). We denote M ≡ M N in N and M ≡ M S in S, and similarly denote (6.16) The sign ± accounts for the quantum number τ in Eq. (6.14).
For a surface state we need M N −β < −|U ± N |, M S −β > |U ± S | in some interval of E, k y , k z around zero. Solution of Eq. (6.14) gives the wave function profile , for x < 0, (6.18) with inverse decay lengths on the normal and superconducting sides of the NS interface. The amplitudes C N and C S are to be adjusted so that ψ(x) is continuous at x = 0. By requiring that the matrix of coefficients of the mode-matching equations has vanishing determinant, we arrive at the dispersion relation of the surface modes, (6.20) discarding terms of second order in µ, λ. The level crossing at k z = 0, for a given k y , happens for m 0 = t(cos k y −1). The corresponding charge expectation value Q = −e∂E/∂µ is one order in µ, λ less accurate than the energy. In Fig. 8 we compare the numerical diagonalization of the Hamiltonian (6.1) with the analytical mode matching calculation. Unlike the comparison in Fig. 4, here there is not a single fit parameter. The agreement is excellent for the energy, somewhat less for the average charge.

E. Effective surface Hamiltonian
In Sec. III we constructed an effective surface Hamiltonian by relying only on particle-hole symmetry. As an alternative route, we present here a derivation starting from the model Hamiltonian (6.8).
The motion perpendicular to the NS interface at x = 0 is governed by the reduced Hamiltonian with neglect of the terms ∝ µ, λ as well as the k y and k z -dependent terms for motion parallel to the interface. The wave function profile ψ(x) at E = 0, decays for x → −∞ (inside the superconducting region) because of the term ∝ M (−∞) > β and for x → +∞ (inside the Weyl semimetal region) because of the term ∝ β > M (∞). This two-sided decay is ensured if ψ(0) is an eigenstate with eigenvalue +1 of both ν 0 τ 0 σ y and ν z τ 0 σ y . The resulting eigenspace has rank two. The 2×2 effective surface Hamiltonian H eff for motion parallel to the surface is obtained by projecting H onto this two-dimensional eigenspace, resulting in The corresponding charge operator is momentum dependent, In this effective surface description the energy scales ∆ and µ should be regarded as weighted averages of the x-dependent parameters from Eq. (6.15).
The two surface modes have opposite charge Q ± = ±e (1 − ∆ 2 /M 2 ) 1/2 and dispersion relation E ± (k z ) = t sin k y − (∆ 2 +m 2 (k y , k z ) + t 2 z sin 2 k z ) −1/2 × λt z sin k z ± µ m 2 (k y , k z ) + t 2 z sin 2 k z , (6.26) representing the spiraling surface Fermi arc illustrated in Fig. 1. The ± index corresponds to the ν index of Sec. III, the κ index is taken care of by the sign of sin k y . The We interpret m eff as the effective coupling strength of the surface state to the superconductor, and as the parameter that in the microscopic model of Sec. IV is varied by varying U barrier . The level crossing then happens when m eff = 0. At the level crossing the excitations are charge neutral.
We may include the Schrieffer-Wolff correction, by projecting δH SW from Eq. (6.13) onto the surface eigenspace. The result is a correction of order δ 2 to the effective surface Hamiltonian, The dominant effect of this correction is to shift the level crossing away from k z = 0 to k z = −(λ/β)(t/t z ) sin k y .

VII. CONCLUSION
In summary, we have investigated the superconducting proximity effect on the dispersion relation of surface modes in a Weyl-Majorana solenoid -a Weyl semimetal nanowire with an axial magnetization covered by a superconductor. The surface Fermi arc connecting bulk Weyl cones is broken up into nearly charge-neutral Majorana modes. We have identified a "connectivity index" that determines between which pair of modes a gap is opened by the superconductor.
We have discovered that the sign of the induced gap can be inverted by variation of the tunnel coupling strength between the semimetal and the superconductor. A domain wall separating segments of the nanowire with opposite sign of the gap traps a charge-neutral quasiparticle. This bound Majorana fermion is not at zero energy, so it should not be confused with the Majorana zero-modes in semiconductor nanowires [27][28][29]. The gap inversion is studied for a 3D model Hamiltonian, both numerically in a tight-binding formulation, and analytically via mode matching at the normal-superconductor interface. Further insight is obtained by an effective 2D surface Hamiltonian.
In closing we remark on a global aspect of the gap inversion in terms of the flow of Berry curvature (topological charge) in the Brillouin zone [30]. The minimal number of two Weyl cones in a Weyl semimetal with broken time-reversal symmetry is doubled if we include the electron-hole degree of freedom. The sign of the Berry curvature at a given point in the Brillouin zone is not changed by the doubling [8], so the Fermi arc connecting Weyl cones of opposite Berry curvature must still run across the Brillouin zone -but now it has a choice: it may connect cones of the same or opposite electrical charge. If we inspect Fig. 4 we see that the Fermi arcs always connect Weyl cones of the same electrical charge (coded blue or red), except at the gap inversion point. At the critical tunnel barrier height U barrier = U c the Majorana surface modes connect bulk states of opposite electrical charge (from blue to red).
In Fig. 4 the anomalous connection by Fermi arcs of Weyl cones of opposite electrical charge and opposite topological charge happens only at an isolated point in parameter space, because the superconductivity is induced only at the surface of the Weyl semimetal. By inducing superconductivity throughout the bulk (for example, using the heterostructure approach of Ref. 8) one should be able to stabilize the anomalous connection in an entire region of parameter space. We expect an anomalous Josephson effect to develop in the Weyl-Majorana solenoid as a result of this topologically nontrivial connection.
where we abbreviated For simplicity we omitted V b (x) from the mode-matching calculations and the derivation of the effective surface Hamiltonian in Sec. VI. In the following we include it in the calculation, resulting in an improved agreement of the analytics with the numerics but without simple closed-form expressions as Eqs. (6.20) and (6.21). The step-function variation of the pair potential ∆(x) at the NS interfaces x = 0, W produces a delta-function boundary potential. Let us focus on the interface at x = 0, with ∆ = ∆ N for x > 0 and ∆ = ∆ S for x < 0. Because of the boundary potential, the wave function does not vary continuously across the NS interface. Instead, the wave functions at the two sides of the interface x = 0 are related by the transfer matrix, ψ(0 + ) = e iMNS ψ(0 − ), where the angle α is given by the integral Note that at the level crossing point we have m z = 0 hence α = 0, so the level crossing itself is not affected by the boundary potential. As explained in Sec. VI E, to obtain the effective surface Hamiltonian we impose a two-sided decay of the wave function, by demanding that ψ is an eigenstate with eigenvalue +1 of ν 0 τ 0 σ y in S and of ν z τ 0 σ y in N. The former condition can be rewritten as a boundary condition in N, Note that U b and ν z τ 0 σ y commute, so they can be diagonalized simultaneously. The rank two eigenspace of eigenvalue +1 is spanned by the vectors v 1 = (0, 0, sin α, i sin α, 1 − cos α, −i + i cos α, 0, 0) , v 2 = (sin α, i sin α, 0, 0, 0, 0, 1 − cos α, −i + i cos α) .
The Hamiltonian projected onto this eigenspace is H eff = τ 0 t sin k y − (γ/M )(λτ 0 t z sin k z − µτ z m z ), γ = cos α + (∆/m z ) sin α, FIG. 9: Colored data points: Energy spectrum (color scale as in Fig. 2) and average charge obtained from a numerical diagonalization of the discretized Hamiltonian (6.1). The parameters are the same as in Fig. 8. The black dashed curves result from the mode-matching calculations including the boundary potential and the full Hamiltonian (with the off-diagonal terms).
where the x-dependent gap ∆(x) in the full Hamiltonian has been replaced by a spatial average∆, andM = (m 2 z + ∆ 2 ) 1/2 .
Comparison with Eq. (6.24) shows that the effect of the boundary potential is to renormalize the parameters λ and µ by a factor γ. For ∆ S m z we have The full mode-matching calculation of Sec. VI D is also modified by the new boundary condition. Since Eq. (A3) mixes the ν and τ indices, we can no longer use the blockdiagonalization of the Hamiltonian to simplify the mode matching, and we could not find a closed-form solution analogous to Eqs. (6.20) and (6.21). Including both the diagonal and off-diagonal terms in the Hamiltonian (6.8) we find the energy and charge expectation value shown in Fig. 9 (dashed curves). The solid curves are the numerical solution of the tight-binding model. Comparison with Fig. 8, where we did not include the bound-ary potential and discarded off-diagonal ν, τ terms in the Hamiltonian, shows little difference in the energy but an improved agreement in the charge.