A numerically exact study of Weyl superconductivity

TL;DR

This paper investigates the effects of interactions on Weyl superconductivity using exact quantum Monte Carlo simulations, revealing unique pairing mechanisms influenced by pseudo-spin and topology.

Contribution

It introduces a numerically exact study of Weyl superconductivity in a modified Hubbard model, highlighting the role of pseudo-spin in pairing and topological charge.

Findings

Pseudo-spin influences pairing mechanisms.

Attractive interactions lead to mixed singlet and triplet pairs.

Provides detailed characterization of exotic pairing in topological systems.

Abstract

We study the interplay of interactions and topology in a pseudo-spin Weyl system, obtained from a minimally modified Hubbard model, using the numerically exact auxiliary-field quantum Monte Carlo method complemented by mean-field theory. We find that the pseudo-spin plays a key role in the pairing mechanism, and its effect is reflected in the structure of the pairing amplitude. An attractive on-site interaction leads to pairing between quasiparticles carrying opposite spin and opposite topological charge, resulting in the formation of real-spin singlet pairs that are a mixture of pseudo-spin singlet and pseudo-spin triplet. Our results provide a detailed characterization of the exotic pairing behavior in this system, and represent an important step towards a more complete understanding of superconductivity in the context of topological band structures, which will help guide searches for topological superconductivity in real materials and ultracold atoms.

Tables

Table 1. Table 1: Comparison of exact diagonalization results with AFQMC results.

$N$		$L_x$		$L_y$		$U$		$v$		$w$		$t_d$		$E_{\text{ED}}$		$E_{\text{QMC}}$
8		4		3		-8.0		0.5		1.3		0.6		-36.64137		-36.6459(80)
8		4		3		-8.0		0.2		1.0		0.3		-34.30740		-34.3099(76)
2		4		2		-12.0		0.2		0.8		0.4		-12.68698		-12.6820(43)
4		4		2		-12.0		0.8		1.0		0.9		-28.69211		-28.6930(70)

Equations

$${\hat{H}}_0$$

$${\hat{H}}_0=\sum _{\mathbf{k}\sigma }{\mathbf{c}}_{\mathbf{k}\sigma }^{†}\mathcal{H}(\mathbf{k}){\mathbf{c}}_{\mathbf{k}\sigma },$$

$$\mathbf{h}(\mathbf{k})=\left(\begin{array}{c}{v}_1+w_1\cos (k_{x}a)\\ w_1\sin (k_{x}a)\\ 0\end{array}\right),$$

$$\left({\chi }_{\mathbf{k}\sigma }^{(-)†},{\chi }_{\mathbf{k}\sigma }^{(+)†}\right)=\left(c_{\mathbf{k}\sigma }^{(\text{A})†},c_{\mathbf{k}\sigma }^{(\text{B})†}\right)\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}-e^{-i{\theta }_{{\mathbf{h}}_{\mathbf{k}}}/2}& & e^{-i{\theta }_{{\mathbf{h}}_{\mathbf{k}}}/2}\\ e^{i{\theta }_{{\mathbf{h}}_{\mathbf{k}}}/2}& & e^{i{\theta }_{{\mathbf{h}}_{\mathbf{k}}}/2}\end{array}\right),$$

$${\hat{H}}_0$$

$${ϵ}_{\mathbf{k}}^{(\pm )}$$

$${\hat{H}}_{\text{I}}=\sum _{i,\alpha }U{n}_{i\uparrow }^{\alpha }{n}_{i\downarrow }^{\alpha },$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{\text{I}},$$

$$|{\Psi }_0⟩\propto \underset{\beta \to \infty }{\lim }{e}^{-\beta \hat{H}}|{\Psi }_T⟩,$$

$${\hat{H}}_{MF}$$

$${\Delta }^{(\alpha \beta )}$$

$${\Delta }^{\text{AA}}={\Delta }^{\text{BB}},{\Delta }^{\text{AB}}={\Delta }^{\text{BA}}.$$

$${\hat{H}}_{MF}$$

$$\Delta$$

$$1$$

$$E_{\mathbf{k}}^{(\pm )}$$

$${\gamma }_{\mathbf{k}\uparrow }^{(\pm )}$$

$${\gamma }_{\mathbf{k}\downarrow }^{(\pm )†}$$

$$u_{\mathbf{k}}^{(\pm )}$$

$$v_{\mathbf{k}}^{(\pm )}$$

$${\mathbf{S}}_{\mathbf{k}}=\frac{1}{2}\sum _{\alpha ,\beta ,\sigma }{\mathbf{\sigma }}_{\alpha ,\beta }{c}_{\mathbf{k}\sigma }^{†(\alpha )}{c}_{\mathbf{k}\sigma }^{(\beta )},$$

$$M_{\mathbf{kq}}^{\mu \nu {\mu }^{\prime }{\nu }^{\prime }}=⟨{\Delta }_{\mathbf{k}}^{†\mu \nu }{\Delta }_{\mathbf{q}}^{{\mu }^{\prime }{\nu }^{\prime }}⟩,$$

$${\Delta }_{\mathbf{k}}^{†\mu \nu }=\frac{1}{\sqrt{2}}\left({\chi }_{\mathbf{k}\uparrow }^{(\mu )†}{\chi }_{-\mathbf{k}\downarrow }^{(\nu )†}-{\chi }_{\mathbf{k}\downarrow }^{(\mu )†}{\chi }_{-\mathbf{k}\uparrow }^{(\nu )†}\right),$$

$${\Psi }_{\mathbf{k}}^p={\psi }_{\mathbf{k}}^{++}+{\psi }_{\mathbf{k}}^{--},$$

$${\Psi }_{\mathbf{k}}^p={\psi }_{\mathbf{k}}^{\text{AA}}+{\psi }_{\mathbf{k}}^{{\text{AB}}_t}+{\psi }_{\mathbf{k}}^{{\text{AB}}_s}+{\psi }_{\mathbf{k}}^{\text{BB}},$$

$${\Delta }_{\mathbf{k}}^{†\text{AA}}=\frac{1}{\sqrt{2}}\left(c_{\mathbf{k}\uparrow }^{(†\text{A})}{c}_{-\mathbf{k}\downarrow }^{(†\text{A})}-c_{\mathbf{k}\downarrow }^{†(\text{A})}{c}_{-\mathbf{k}\uparrow }^{†(\text{A})}\right);$$

$${\Delta }_{\mathbf{k}}^{†\text{BB}}=\frac{1}{\sqrt{2}}\left(c_{\mathbf{k}\uparrow }^{†(\text{B})}{c}_{-\mathbf{k}\downarrow }^{†(\text{B})}-c_{\mathbf{k}\downarrow }^{†(\text{B})}{c}_{-\mathbf{k}\uparrow }^{†(\text{B})}\right);$$

$${\Delta }_{\mathbf{k}}^{†{\text{AB}}_t}=\frac{1}{2}\{\left(c_{\mathbf{k}\uparrow }^{†(\text{A})}{c}_{-\mathbf{k}\downarrow }^{†(\text{B})}-c_{\mathbf{k}\downarrow }^{†(\text{A})}{c}_{-\mathbf{k}\uparrow }^{†(\text{B})}\right)$$

$$+\left(c_{\mathbf{k}\uparrow }^{†(\text{B})}{c}_{-\mathbf{k}\downarrow }^{†(\text{A})}-c_{\mathbf{k}\downarrow }^{†(\text{B})}{c}_{-\mathbf{k}\uparrow }^{†(\text{A})}\right)\};$$

$${\Delta }_{\mathbf{k}}^{†{\text{AB}}_s}=\frac{1}{2}\{\left(c_{\mathbf{k}\uparrow }^{†(\text{A})}{c}_{-\mathbf{k}\downarrow }^{†(\text{B})}-c_{\mathbf{k}\downarrow }^{†(\text{A})}{c}_{-\mathbf{k}\uparrow }^{†(\text{B})}\right)$$

$$-\left(c_{\mathbf{k}\uparrow }^{†(\text{B})}{c}_{-\mathbf{k}\downarrow }^{†(\text{A})}-c_{\mathbf{k}\downarrow }^{†(\text{B})}{c}_{-\mathbf{k}\uparrow }^{†(\text{A})}\right)\}.$$

$${\mathbf{h}}_{\mathrm{linear}}(\mathbf{k})=\left(\begin{array}{c}{a}_1{k}_y\\ a_2{k}_x\\ 0\end{array}\right),$$

$${\mathbf{A}}_{\pm }(\mathbf{k})=\frac{a_1{a}_2}{2a_1^2{k}_y^2+2a_2^2{k}_x^2}\left[\begin{array}{c}{k}_y\\ -k_x\\ 0\end{array}\right].$$

$${\mathbf{h}}_{\mathrm{quadratic}}(\mathbf{k})=\left(\begin{array}{c}-b_1{k}_x^2+b_2{k}_y^2\\ 2b_1{k}_x\\ 0\end{array}\right),$$

$$e^{-\beta \hat{H}}={\left(e^{-\delta \tau \hat{H}}\right)}^m,$$

Full text

I Introduction

Pairing and superconductivity have been focal points of condensed matter physics for several decades, whereas the comparatively modern discovery of topological materials has generated intense recent interest in the role of topology in condensed matter systems. The microscopic, many-body origins of superconductivity are the subject of an expansive body of theoretical and computational work, while exhaustive schemes have been developed to classify topological systems based on symmetries of the Hamiltonian and band structure. There is a fundamental conceptual division between these efforts; pairing is a many-body behavior that only emerges in the presence of electron interactions, while topological behavior is typically characterized by the non-interacting one-body picture. A successful unification of these concepts is a crucial step towards a more complete treatment of many intriguing problems, including topological superconductivity. Given the enticing potential applications of these ideas in quantum computing and information, a detailed, quantitative description is a compelling open challenge.

Since the inception of the Hubbard model, quantum lattice models have served as a testbed for theories of pairing and superconductivity. One prominent example being the repulsive Hubbard model, which has been at the center of a persistent and ongoing effort to understand the origin and mechanisms of superconductivity in the cuprates. These models are broadly relevant across condensed matter, nuclear, and atomic physics, and most recently in the context of ultracold atoms, which offer the possibility of high-accuracy and finely tunable experimental realizations of a variety of lattice models Cheuk et al. (2016a); Miyake et al. (2013).

While pairing has been a central theme of condensed matter research for over sixty years, the last decade has seen the emergence of new classes of materials with topological character, in particular, the class of Weyl semimetals, that have revolutionized condensed matter physics. These materials are especially interesting because they become superconducting at sufficiently low temperature Guguchia et al. (2017); Rhodes et al. (2017, 2019). The presence of superconductivity in these Weyl systems prompts an intriguing question: to what extent is the superconductivity due to Weyl fermions that seem to be present near the Fermi energy?

Here we introduce a simple lattice fermion model, which is a straightforward extension of the two-dimensional (2D) Hubbard model, with Weyl quasipartlcles near the Fermi surface that interact via on-site attraction. Other models of Weyl systems have been studied either without interactions Yang and Zhang (2016) or at the mean-field level Meng and Balents (2012); Bednik et al. (2015), but these mean-field approaches are approximative and their results can be unreliable. Exact diagonalization, while free of approximations, is limited to very small sized systems, and is therefore unsuitable for the detection of long range correlations and pairing. Quantum Monte Carlo methods, based on the truly underlying degrees of freedom, the electrons, have demonstrated unique capability in the treatment of strongly-correlated many-body systems, including systems with fermionic pairing Shi et al. (2015a, 2016); Rosenberg et al. (2017, 2019); Vitali et al. (2016). However, many lattice models with topological character suffer from the fermion sign problem, so very few unbiased many-body studies of strongly-correlated topological systems have been done Otsuka et al. (2018); Li et al. (2017); Hung et al. (2014); Hohenadler et al. (2011), and a quantitative description of pairing in these systems remains an important goal.

With a minimal modification to the well-known attractive Hubbard model, we obtain a Hamiltonian whose one-body term describes a Weyl system, even in the absence of spin-orbit coupling. The one-body part is a set of stacked 1D SSH chains Su et al. (1979) coupled via an inter-chain hopping term to form a 2D lattice. The addition of interactions, which render the model analytically intractable and exceptionally computationally challenging, leads to the emergence of several exotic behaviors including pairing between Weyl quasiparticles. We find that with suitable hopping parameters the system supports a phase with pairs composed of quasiparticles from the same pseudo-helicity branch, but of opposite topological charge (chirality). We study the model using the cutting-edge auxiliary-field quantum Monte Carlo (AFQMC) method, complemented by mean-field theory calculations, for a variety of system sizes to obtain good control over the finite-size effects. The AFQMC and mean-field approaches are in reasonable agreement. The Hamiltonian preserves time reversal symmetry, which guarantees that our AFQMC calculations are free of the sign problem Wu and Zhang (2005), meaning that these results are numerically exact and offer a uniquely accurate, detailed quantitative description of pairing behavior in a strongly-correlated topological system. In addition to demonstrating the capability of the AFQMC technique to treat interacting topological systems, this work can serve as a guide to ongoing experimental efforts in topological superconducting materials.

The quantitative picture we develop yields several observations about the nature of pairing in this topological model. We find that the attractive interaction leads to spin singlet pairing of Weyl quasiparticles of the same helicity but opposite topological charge into a BCS-like superconductivity, which emerges at a small but finite interaction strength. We also examine in detail the pairing mechanism and its effects on the spatial variation of the pairing amplitude. While the constituents of each pair have non-trivial topological character, the attractive on-site interaction leads to a topological-charge-neutral form of superconductivity. We also discuss pairing mechanisms capable of producing topological superconductivity in this Weyl system Li and Haldane (2018), which is a direction of future research.

The remainder of the paper is organized as follows, In section II we introduce our lattice model. In section III we describe our approach, which combines AFQMC and mean-field theory calculations. We present our results in section IV, and conclude with a discussion of experimental implications, routes towards topological superconductivity, and an outlook in section V.

II The Model

Our model has the following lattice Hamiltonian,

[TABLE]

where the operator $c^{(\textmd{A})\dagger}_{{n,m}}$ creates an electron on the A site of the unit cell at position $\mathbf{R}=na\,\hat{x}+mb\,\hat{y}$. Each unit cell is composed of an A site and a B site. The intra- and inter-unit-cell hopping strengths in the $\hat{x}$-direction are controlled by the parameters $v$ and $w$, respectively, and the diagonal hopping strength is given by $t_{d}$.

In momentum space the Hamiltonian can be written,

[TABLE]

with the vector $\mathbf{c}^{\dagger}_{\mathbf{k}\sigma}=(c^{(\textmd{A})\dagger}_{\mathbf{k}\sigma},c^{(\textmd{B})\dagger}_{\mathbf{k}\sigma})$, and the matrix $\mathcal{H}(\mathbf{k})=-\mathbf{h}(\mathbf{k})\cdot\boldsymbol{\sigma}$. This has the form of a Weyl Hamiltonian, where $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector of the Pauli matrices in the sublattice (or pseudo-spin) basis, and

[TABLE]

with $v_{1}=v+2t_{d}\cos(k_{y}b)$, and $w_{1}=w+2t_{d}\cos(k_{y}b)$. We plot the lattice geometry and band structure in Fig. 1. For $v+w<4t_{d}$ there are two Weyl nodes at $\mathbf{k}=(0,\pm k_{N})$, with $k_{N}=\cos^{-1}(-(v+w)/4t_{d})$. These nodes are protected by the combined inversion (A $\rightarrow$ B, $\mathbf{k}\rightarrow-\mathbf{k}$) and time-reversal symmetry of the lattice Hamiltonian. For a more detailed characterization of the topological features of the model see Appendix A.

The connection between the pseudo-spin basis and the diagonal pseudo-helicity basis is defined by the unitary transformation,

[TABLE]

which introduces a new set of creation operators, $\chi^{(\pm)\dagger}_{\mathbf{k}\sigma}$, with the angle, $\theta_{\mathbf{h}_{\mathbf{k}}}\equiv\tan^{-1}(h^{y}_{\mathbf{k}}/h^{x}_{\mathbf{k}})$. In the pseudo-helicity basis the tight-binding Hamiltonian takes the form:

[TABLE]

This simple lattice model with a topological band structure provides an ideal setting in which to study the interplay of topology and superconductivity at the many-body level. With this motivation in mind, we proceed by introducing interactions into the system in the form of an on-site attractive Hubbard term,

[TABLE]

where $i$ labels the unit cell, $\alpha=\textrm{A},\textrm{B}$, and the interaction parameter $U<0$. Our full Hamiltonian is now,

[TABLE]

with $\hat{H}_{0}$ defined by Eq. (2), and $\hat{H}_{I}$ defined by Eq. (7).

III Methods

III.1 Auxiliary-field quantum Monte Carlo

Strongly correlated many-body systems are a well-known theoretical and computational challenge. One method that has demonstrated considerable accuracy in the treatment of these types of systems is auxiliary-field quantum Monte Carlo (AFQMC) Sugiyama and Koonin (1986); Zhang (2013); Assaad (2002). The method has been widely applied to both quantum chemistry Motta and Zhang (2018); Landinez Borda et al. (2019) and model Hamiltonians Qin et al. (2016); Vitali et al. (2019, 2016); Shi et al. (2015a, 2016); Sun et al. (2013). In sign-problem-free cases, such as the attractive spin-unpolarized Hubbard model, AFQMC calculations are numerically exact. The technique can be used to calculate an array of ground state properties and provide a reliable quantitative many-body description of the pairing, spin, and charge behaviors of strongly correlated systems.

The AFQMC algorithm is built on the idea of imaginary-time projection, by which the many-body ground-state, $|\Psi_{0}\rangle$, of a Hamiltonian, $\hat{H}$, can be computed via application of a many-body projection operator to an initial trial wave function, $\left|\Psi_{T}\right>$,

[TABLE]

where $\beta$ is imaginary time, and we require $\left<\Psi_{T}\vphantom{\Psi_{0}}|\Psi_{0}\vphantom{\Psi_{T}}\right>\neq 0$ (the trial wave function must have non-vanishing overlap with the many-body ground-state). In the AFQMC framework, the many-body projection operator is decomposed into a set of one-body operators coupled to auxiliary fields. The projection process in Eq. (9) is then recast as a path integral in auxiliary-field space that can be evaluated using Monte Carlo techniques. See Appendix B for details on the AFQMC method. Our calculations treat periodic lattices with over 700 sites and 650 electrons to provide a systematic and high-accuracy characterization of the ground-state properties of the system.

III.2 Bogoliubov-Valatin-de Gennes-BCS theory

To complement the AFQMC approach outlined above, we have performed a set of calculations within the BCS mean-field theory (MFT) framework. The many-body Hamiltonian in Eq. (8) can be written in quadratic form using an appropriate mean-field decomposition,

[TABLE]

with,

[TABLE]

where $\alpha,\beta$ take the values A and B, $\mu$ is the chemical potential, $\hat{N}$ counts the total number of electrons, and $N_{c}$ is the number of unit cells. The Hamiltonian preserves inversion symmetry (${\bf k}\to-{\bf k}$ and A $\to$ B), which, taking into account the summation over $\mathbf{k}$ in Eq. (11), imposes the following identities:

[TABLE]

In the basis of the helicity eigenstates:

[TABLE]

with $\epsilon^{(\pm)}_{\bf k}$ defined in Eq. (6) and,

[TABLE]

The identities in Eq. (12) guarantee that pairing only occurs between quasiparticles from the same pseudo-helicity branch. Therefore, in this basis, the mean-field Hamiltonian is block-diagonal, and the gap equation takes the form of the standard BCS gap equation:

[TABLE]

with the Bogoliubov-Valatin quasiparticles given by the well-known transformation:

[TABLE]

where,

[TABLE]

IV Results

IV.1 Mean-field theory calibration

As outlined in the preceding section, our approach combines cutting-edge many-body AFQMC calculations with mean-field theory. Our model is free of the fermion sign problem, which means that the AFQMC results are numerically exact and can be used to calibrate the mean-field theory on finite-size systems. We begin this section with a comparison of AFQMC and MFT calculations of the pairing amplitude (defined in Sec. IV.3) across several parameter sets. This comparison provides a qualitative validation of the mean-field description, which is an important complement to the AFQMC analysis, and also offers an estimate of any finite-size effects in the AFQMC calculations. We find that the mean-field treatment correctly captures the prominent features of the AFQMC results, and that the finite-size effects are generally small. We illustrate below the quality of agreement between the MFT and AFQMC results with several representative examples; additional examples are given in Appendix D.

In Fig. 2 we present a comparison of results for the pairing amplitude in the pseudo-helicity basis from AFQMC and MFT. There is good qualitative agreement between the calculations, which show that pairing occurs primarily in the vicinity of the non-interacting Fermi surface, and that there is relatively little occupation of momentum states above the non-interacting Fermi surface. The MFT result shows a smoother momentum distribution, with occupation spreading above the non-interacting Fermi surface towards the node at $\mathbf{k}=(0,\pm k_{N})$, and a correspondingly smooth pairing amplitude, but the essential features of the AFQMC result are evident in the MFT calculation. As we illustrate in Fig. 3 this qualitative agreement holds for increased interaction strength and different values of the hopping asymmetry, $|v-w|$. Having established a reasonable level of agreement between the AFQMC and MFT results, we can extend our discussion to systems directly in the thermodynamic limit with a certain degree of confidence. Because of their size, these systems are computationally inaccessible to AFQMC, but they are manageable with MFT.

IV.2 Pseudo-spin distribution and bond-order

Several features that emerge from the lattice Hamiltonian inform our investigation of the pairing behavior in this system. We first consider the pseudo-spin distribution, presented in the left panel of Fig. 5. We plot the pseudo-spin direction, $(\langle S^{x}_{\mathbf{k}}\rangle,\langle S^{y}_{\mathbf{k}}\rangle)$, calculated from the many-body expectation values of the operators,

[TABLE]

where $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector of the Pauli matrices in pseudo-spin space. On the right side of Fig. 5 we show the momentum space “effective field” $\bf{h}({\bf k})$ given by Eq. (3). Notice that the pseudo-spin is parallel to the effective field, even in the presence of interactions. This implies that the Weyl quasiparticles in the system experience this field, and the pairing mechanism should be consistent with this effect; in general, different forms of superconducting order can compete and the ground state represents the most energetically favorable pairing state in the presence of this pseudo-spin distribution. We elaborate on this connection between the pairing mechanism and the pseudo-spin distribution later in this section.

We next study the effect of the asymmetric hopping along the $\hat{x}$-direction, which results in an oscillation of the bond density, defined by the operator $\rho^{\mu}_{\ell}=c^{\dagger}_{\mathbf{r}_{\ell}}c_{\mathbf{r}_{\ell}+\hat{\mu}}+\textrm{h.c.}$, that measures the density on the $\ell$-th bond along the $\mu$-direction, where the first site of the bond is located at position $\mathbf{r}_{\ell}$. We compute the bond-density correlation function, $\langle\rho^{\mu}_{\ell}\rho^{\nu}_{m}\rangle$, for several values of hopping asymmetry and interaction strength. A typical example is presented in Fig. 4. We observe a clear oscillation along the $\hat{x}$-direction that is proportional to the hopping asymmetry (see upper panel in Fig. 4(b)). A similar oscillation, with smaller amplitude, is evident along the diagonal directions (lower panel in Fig. 4(b)). The magnitude of this order is relatively insensitive to the strength of the interaction, which suggests that the behavior is a primarily one-body effect, though as we illustrate later, the hopping asymmetry responsible for this bond-density wave is also intricately related to the pairing behavior in real and momentum space.

IV.3 Pairing from a topological band structure

We now present a detailed picture of the pairing behavior. Our description focuses on the pairing amplitude, which is the eigenstate corresponding to the leading eigenvalue of the two-body density matrix Yang (1962). We define the elements of this matrix in the pseudo-helicity basis as:

[TABLE]

where,

[TABLE]

and $\mu,\nu=\pm$. We note that with an attractive on-site interaction, in the absence of spin-orbit coupling, there is no spin-triplet pairing. Additionally, pairing between pseudo-helicity quasiparticles from different branches ($\mu\neq\nu$ above) is identically zero, so the pairing matrix can be written in $2N_{c}\times 2N_{c}$ form. In this basis the pairing amplitude is,

[TABLE]

where $\psi^{+{}+}_{\mathbf{k}}$ ($\psi^{-{}-}_{\mathbf{k}}$) corresponds to pairs formed from quasiparticles in the upper (lower) pseudo-helicity branch.

In the pseudo-spin basis the pairing amplitude can be written,

[TABLE]

which corresponds to the following set of pair creation operators,

[TABLE]

As in the pseudo-helicity basis, all of these operators create spin singlets. The first three create pseudo-spin triplet pairs, and the fourth creates a pseudo-spin singlet.

In the remainder of this section, we analyze in quantitative detail the effects of both interaction and hopping asymmetry on pairing in this system. All the results we present are for densities slightly below half-filling, $n\approx 0.94$, and interaction strengths consistent with the energy scales in relevant materials Rhodes et al. (2017). The Hamiltonian has particle-hole symmetry, which establishes a direct mapping between electron-doped and hole-doped systems, and implies that the results are symmetric about half-filling.

Figure 6 summarizes the behavior of the pairing amplitude in the pseudo-helicity basis versus interaction strength and hopping asymmetry (see Appendix E for other parameter values). This quantity is directly related to the gap function, which can be measured via STM or Bogoliubov quasiparticle interference Allan et al. (2013); Sprau et al. (2017); Wang and Lee (2003). We see significant amplitude for pairing between quasiparticles from the lower helicity band, mostly confined to the vicinity of the non-interacting Fermi surface. Some pairing is evident away from the Fermi surface, which is a reflection of the underlying band structure. The regions where the pairing amplitude has significant magnitude away from the Fermi surface correspond to states near the edges of the Brillouin zone that surround the deep minimum of the lower helicity band. These states are close in energy to the Fermi energy (see Fig. 1) and therefore participate in pairing. Interestingly, the AFQMC calculations reveal a small amplitude for pairing in the upper helicity band, with a peak near the nodal point (Fig. 6(b)). This type of pairing only appears in the mean-field picture for interaction strengths $|U|\gtrsim 2.0$. As the hopping asymmetry decreases, the curvature of the elliptical Fermi surface decreases, which brings the two sides of the Fermi surface closer together. The effect of this change to the geometry of the Fermi surface on the pairing behavior is evident in the magnitude of the pairing amplitude along the lower Fermi surface (Fig. 6(c)). As the hopping asymmetry decreases, there is increased pairing inside the Fermi surface, close to the nodal point, and the pairing is of equal magnitude along either side of the Fermi surface. The pairing amplitude shows similar quantitative changes with increasing interaction strength, which leads to larger values of the pairing amplitude along the Fermi surface. The mean-field description in the thermodynamic limit is consistent with the behavior seen in the AFQMC calculations on finite-size systems (see Fig. 7). Here we find that the pairing amplitude is essentially constant along the non-interacting Fermi surface, and that the pairing amplitude away from the Fermi surface is closely connected to the underlying band structure. As the hopping asymmetry decreases, the helicity bands become flatter and the pairing amplitude away from the Fermi surface increases while the pairing amplitude along the Fermi surface decreases.

Several interesting features of the pairing behavior emerge in the pseudo-spin basis. We explore the connection between the pseudo-spin degree of freedom and pair formation in Fig. 8. We see from the pseudo-spin distribution that there is a $2\pi$ vortex-like rotation of the pseudo-spin around the upper node at $\mathbf{k}=(0,k_{N})$ and a corresponding rotation with opposite vorticity around the lower node at $\mathbf{k}=(0,-k_{N})$. This behavior is related to the pairing mechanism, sketched by the circles and connected arrows, which illustrate the formation of zero-momentum pairs along the inner and outer edges of the Fermi surface. Each pair is a real-spin singlet, as well as a mixture of pseudo-spin singlet and pseudo-spin triplet. Pairs formed by electrons in states on the outer edges of the Fermi surface have a large net pseudo-spin triplet component because the electrons in these states have nearly parallel pseudo-spin, and consequently this component of the pairing amplitude has a large magnitude. The same effect diminishes the pairing amplitude in the pseudo-spin singlet sector along the outer edges of the Fermi surface. Along the inner edges, the net pseudo-spin triplet content of the pair is smaller and goes to zero before changing orientation, which leads to a node in the triplet sector of the pairing amplitude. There is also pairing in the pseudo-spin singlet sector along the inner edges of the Fermi surface. This sector of the pairing amplitude exhibits a node at $k_{x}=0$, corresponding to the momenta at which the pair becomes a net pseudo-spin triplet.

We present a typical example of the pairing amplitude in the pseudo-spin basis in Fig. 9. The left panel plots the AA-sector of the pairing amplitude. This component has significant amplitude along the Fermi surface, and near the edges of the BZ. As in the pseudo-helicity basis, for finite-size systems, the pairing amplitude is larger along the outer edges of the Fermi surface relative to the inner edges. In the middle panel we show the pseudo-spin triplet and singlet sectors, which illustrate the effects of the pseudo-spin degree of freedom in the formation of pairs, as described above. While $\psi^{\textmd{AB}_{t}}_{\mathbf{k}}$ has a large amplitude along the outer edges of the Fermi surface, $\psi^{\textmd{AB}_{s}}_{\mathbf{k}}$ is zero, whereas along the inner edges, both components are non-zero, with $\psi^{\textmd{AB}_{s}}_{\mathbf{k}}$ having a node at $k_{x}=0$. The right panel, which shows each sector of the pairing amplitude along the path in momentum space drawn in Fig. 2, provides another illustration of the nodal structure of the different components, in particular, the node in $\psi^{\textmd{AB}_{t}}_{\mathbf{k}}$ between the inner and outer edges of the Fermi surface.

This exotic pairing behavior also emerges from MFT calculations in the thermodynamic limit, an example of which we present in Fig. 10. We see a clear display of the nodal structures of both the pseudo-spin triplet and singlet sectors of the pairing amplitude. The triplet component shows a sign change along the inner edges of the Fermi surface (left panel of Fig. 10), and the singlet component shows a nodal line, dividing the BZ along the $k_{y}$-axis (middle panel). These nodal structures appear to be largely insensitive to the magnitude of the hopping asymmetry (provided the asymmetry is non-zero).

We find similarly rich pairing properties in real-space. Figure 11 presents the real-space structure of the pseudo-helicity pairing amplitude, given by the Fourier transform of $\psi^{-{}-}_{\mathbf{k}}$. There is a prominent peak at the origin, reflecting the large on-site component of the pair. Away from the origin there are a set of smaller peaks separated by nodes along the $x$ and $y$ axes. We observe that the pairing amplitude becomes more spatially extended as the hopping asymmetry decreases. This behavior is evident in the upper panel of Fig. 11(b) and the inset, which shows $\psi^{-{}-}_{\mathbf{r}}$ along the $x$-axis at decreasing values of $|v-w|$. The central peak is reduced for smaller values of the hopping asymmetry, while the additional peaks along the axis grow. The lower panel illustrates that the central peak of the pairing amplitude increases with interaction strength, but otherwise the qualitative behavior does not change over the range of interaction strengths we have considered. These behaviors are also evident in the MFT description in the thermodynamic limit (Fig. 12).

V Discussion and Conclusion

The discovery of superconducting Weyl materials has ushered in a new era in condensed matter physics. These materials have a number of exciting potential applications in quantum computing and quantum information that have motivated an intense effort to understand their properties. An accurate quantitative description of these systems at the many-body level remains an essential goal. Such a description presents a unique challenge, demanding a unified treatment of topology and strong correlations. In this work we provide a detailed, high-accuracy characterization of the rich pairing behaviors of a strongly-correlated topological system. We find that both the spin and pseudo-spin (sublattice) degrees of freedom play an important role, leading to a pairing amplitude with multiple components that have different spatial behaviors. These behaviors, as well as the pairing mechanism, are connected to the pseudo-spin distribution and bond-density order. We observe that in the pseudo-helicity basis pairing occurs only between quasiparticles from the same pseudo-helicity band carrying opposite topological charge, which results in topological-charge-neutral pairs.

This topologically neutral form of superconductivity is a consequence of the on-site attractive interaction, which favors the formation of spin-singlet pairs with zero net topological charge. The ongoing search for Majorana fermions has indicated that higher order pairing symmetries, for instance $p$-wave pairing, may be an essential ingredient in topological superconductors. Simple modifications to this model, such as the addition of spin-orbit coupling, or non-local interactions, can induce pairing mechanisms with these different symmetries. As we have demonstrated, the AFQMC method is well-suited to treat these types of topological systems, which makes understanding the origins and mechanisms of topological superconductivity at the many-body level a promising direction for future research.

Our model offers a simple, fundamental description of the pairing properties that emerge in the recently discovered class of superconducting Weyl materials. Modern cold atom experiments also offer an ideal platform to simulate and study lattice models based on these and other materials. There have already been experimental realizations of the Hubbard model Jördens et al. (2008); Brown et al. (2017); Cheuk et al. (2016b); Parsons et al. (2016), as well as more exotic variations, including Weyl systems Atala et al. (2013); Miyake et al. (2013); Aidelsburger et al. (2013); Dubček et al. (2015). These experiments offer a clean, highly tunable setting in which to explore the intersection of topology and interaction Cooper et al. (2019). Our results provide important guidance to this next generation of experiments exploring the intersection of topological band structures and interaction, and in turn these experiments can serve as a testbed for numerical many-body approaches probing interaction effects in strongly-correlated topological systems.

In summary, we have presented an illustrative and quantitative description of the pairing properties of an interacting topological system. Our model contains a pair of Weyl nodes that generate a dipole-like pseudo-spin distribution in momentum space. This distribution exhibits a strong effect on the nature of the pairing amplitude, which represents the most energetically favorable paired state in the presence of the distribution. As a consequence, we observe pairing between quasiparticles that carry opposite topological charge, resulting in a topological-charge-neutral pairing amplitude composed of real-spin singlets that are a mixture of pseudo-spin singlet and triplet. We provide a thorough characterization of these exotic pairing behaviors using a combination of numerically exact AFQMC calculations and mean-field theory in the thermodynamic limit. Our results demonstrate the power of the AFQMC method to treat strongly-interacting topologically non-trivial systems at the many-body level. This is an important step towards a more complete understanding of pairing and strong correlations in the context of topological band structures, a subject that has captivated the condensed matter community given the potentially impactful applications of topological superconductivity across quantum computing and quantum information.

VI Acknowledgments

This work was supported in part by the U.S. National High Magnetic Field Laboratory, which is funded by NSF DMR-1157490 and the State of Florida. Computations were performed in part using the facilities at the College of William and Mary which were provided by contributions from the National Science Foundation, the Commonwealth of Virginia Equipment Trust Fund and the Office of Naval Research.

Appendix A Topological characterization of the non-interacting model: Berry phase and Fermi arcs

In the following section we provide a detailed description of the topological features of our model. As alluded to in the main text, the lattice Hamiltonian given in Eq. (3) takes the form of a Weyl Hamiltonian with a linear dispersion in the vicinity of the band touching points. These nodal points occur at $\mathbf{k}=(0,\pm k_{N})$, with $k_{N}=\cos^{-1}(-(v+w)/4t_{d})$, provided $v+w<4t_{d}$. This Weyl form emerges from a Taylor expansion of $\mathbf{h}(\mathbf{k})$ around the nodal points, which yields,

[TABLE]

with $a_{1}=\pm 4t_{d}\mathrm{sin}(k_{N})$, and $a_{2}=(w-v)/2$. In this case, to leading order in $k_{x}$ and $k_{y}$, the dispersion around each node is linear. In the $k_{x}$-$k_{y}$ plane, the dispersion is elliptical, with the degree of anisotropy determined by the ratio of $a_{1}$ to $a_{2}$.

The Berry potential provides another measure of the topological character of this band structure. In terms of the pseudo-helicity eigenstates, $|\pm\rangle$, the Berry potential vector $\mathbf{A}_{\pm}(\mathbf{k})=\langle\pm|\boldsymbol{\nabla}_{\mathbf{k}}|\pm\rangle$ in the linear regime is given by:

[TABLE]

As noted in the main text, the Hamiltonian preserves both time-reversal and inversion symmetry, which leads to a vanishing Berry curvature at all points in the BZ, with the exception of the nodal points, where the Berry curvature diverges. Despite this property of the Berry curvature, the Berry potential and the Berry phase around either node are non-zero. Figure 13 plots the vector field $\mathbf{A}_{-}(\mathbf{k})$ around a single Weyl point for two different values of anisotropy in $a_{1}$ and $a_{2}$. The vector field around the other Weyl point has equal magnitude but opposite vorticity. The top panel of Fig. 13 corresponds to a small anisotropy, while the bottom panel corresponds to a large anisotropy. Note that when $a_{1}=a_{2}$, the vortex pattern of the Berry potential is isotropic. An isotropic Berry potential and its corresponding dispersion with Weyl nodes, can be engineered by an appropriate choice of hopping parameters, for example, setting $a_{1}=a_{2}=1$ gives, $w=2+v$ and $t_{d}=\sqrt{1+4(1+v)^{2}}/4$. The integral of the Berry potential, $\mathbf{A}_{\pm}(\mathbf{k})$, along any path that encloses a single node gives the Berry phase, which, for the two Weyl nodes in this model is equal to $\pm\pi$.

For systems with $v+w\geq 4t_{d}$, the two nodes merge at $\mathbf{k}=(0,\pm\pi)$ and annihilate, which opens a gap in the spectrum without breaking any symmetries of the Hamiltonian. When the equality above holds (i.e. $v+w=4t_{d}$) there is a single nodal point, however the dispersion around this point is quadratic. In this situation the low energy Hamiltonian is given by:

[TABLE]

where, $b_{1}=w-2t_{d}$, and $b_{2}=2t_{d}$. To leading order in $k_{x}$ the dispersion in this case is linear, however, unlike Eq. (27), which is linear in both $k_{x}$ and $k_{y}$, this dispersion is quadratic to leading order in $k_{y}$.

Finally, we consider the surface spectrum and Fermi arc states. In Fig. 14 we present the surface spectrum, calculated for a system with open boundary conditions along either the $\hat{x}$- or $\hat{y}$-direction. Figure 14 plots the dispersion along the $k_{y}$-direction for a ribbon of finite length in the $\hat{x}$-direction. The two Weyl nodes at $k_{y}=\pm k_{N}$ are connected by dispersionless one dimensional modes, the Fermi arc states, which are localized on the two edges of the ribbon. Such one dimensional edge modes are inherited from the zero dimensional edge states of the SSH chain, which is the building block of our lattice. These edge states are absent for a ribbon terminated along the $\hat{y}$-direction, as seen in Fig. 14. However, unlike the SSH chain which hosts edge states only when the lattice terminates on the bond with stronger hopping Asbóth et al. (2016), in our model, edge states are present irrespective of the hopping strength along the terminating bond. The strength of the terminating bond does however determine the orientation of the edge states, as illustrated in Fig. 14.

Appendix B The auxiliary-field quantum Monte Carlo Method

In this section we provide a concise overview of the AFQMC method intended to highlight the essential features of the technique, and refer interested readers to several pedagogical presentations of the formalism for additional details Zhang (2013); Assaad (2002); Sugiyama and Koonin (1986).

Building on our discussion in the main text, the AFQMC algorithm relies on the projection process defined by Eq. (9). This long imaginary-time projection interval is then divided into $m=\beta/\delta\tau$ time slices,

[TABLE]

which establishes an iterative procedure to obtain the limit in Eq. (9),

[TABLE]

with $\left|\Psi^{(0)}\right>=\left|\Psi_{T}\right>$.

Writing our Hamiltonian in the general form, $\hat{H}=\hat{K}+\hat{V}$, where $\hat{K}$ refers to the one-body terms and $\hat{V}$ refers to the two-body terms, we proceed by applying the Trotter-Suzuki decomposition Trotter (1959); Suzuki (1976), followed by a Hubbard-Stratonovich transformation Hubbard (1959):

[TABLE]

[TABLE]

with $\alpha=$ A, B and $\gamma$ is defined by $\cosh(\gamma)=\exp(-\delta\tau\,U/2)$. Here we have chosen to decouple in the charge channel Hirsch (1983), though decompositions in other channels, such as spin or pairing, exist Rosenberg et al. (2019), and the choice of decomposition can affect the efficiency of the simulation Shi and Zhang (2013). This procedure yields the following form for the propagator,

[TABLE]

where $\mathbf{x}=\{x_{1},x_{2},\ldots,x_{N_{S}}\}$ is a set of auxiliary fields at a given time slice, $p(\mathbf{x})$ is a normalized probability density function, and $\hat{B}(\mathbf{x})$ is a one-body operator. The set of auxiliary fields, $\mathbf{x}$, has dimension $N_{S}$ equal to the size of the single-particle basis, which in this case is double the number of unit cells. Note that we have written the propagator in a more general form, as an integral over continuous auxiliary field variables; for the discrete charge decomposition in Eq. (33) this integral is replaced by a summation over discrete auxiliary field variables. With this choice of Hubbard-Stratonovich transformation $p(\mathbf{x})$ is uniform, and,

[TABLE]

with $\hat{b}_{i}(x^{\alpha}_{i})\equiv\exp\left[{(\gamma x^{\alpha}_{i}-\delta\tau U/2)(n^{\alpha}_{i\uparrow}+n^{\alpha}_{i\downarrow}-1)}\right]$. The many-body propagator is now composed of one-body operators with the fermions in external auxiliary fields. The integration over auxiliary field configurations recovers the two-body interactions.

The many-body, ground-state expectation value of an observable $\hat{O}$ is calculated according to,

[TABLE]

Using the definition of the propagator in Eq. (34), the denominator in Eq. (36) can be written,

[TABLE]

where,

[TABLE]

and we have introduced the notation,

[TABLE]

In the above, $\mathbf{x}^{(\ell)}$ represents an auxiliary field configuration at time slice $\ell$, and the collection of auxiliary fields $\mathbf{X}=\{\mathbf{x}^{(1)},\mathbf{x}^{(2)},\ldots,\mathbf{x}^{(m)}\}$ defines a path in auxiliary field space.

The expectation value in Eq. (36) can now be recast as a path integral in auxiliary-field space,

[TABLE]

with,

[TABLE]

This integral can be evaluated using standard Monte Carlo techniques, such as the Metropolis algorithm, which samples auxiliary-fields from $\mathcal{W}(\mathbf{X})$ to obtain a Monte Carlo estimate of the expectation value in Eq. (39). To accelerate the sampling procedure we employ a dynamic force bias Shi et al. (2015b); Zhang (2013), which improves the acceptance ratio and consequently the efficiency of the algorithm. In addition, we remove the infinite variance problem using the bridge link method Shi and Zhang (2016).

Appendix C Calibration of AFQMC results against exact diagonalization

In table 1 we present a comparison of exact diagonalization and AFQMC calculations of the total energy.

In Fig. 15 we present a comparison of exact diagonalization and AFQMC calculations of different correlation functions.

Appendix D Additional comparison of AFQMC and mean-field results

In this section we provide results to supplement the comparisons between AFQMC and MFT calculations shown in the main text. Figure 16 shows a similar level of qualitative agreement, especially in the limit of small hopping asymmetry. Generally, the AFQMC result has larger magnitude near the non-interacting Fermi surface, while the MFT result is larger away from the Fermi surface. We see similar qualitative agreement in Fig. 17, where both methods show a difference in the magnitude of the pairing amplitude along the outer edges of the Fermi surface relative to the inner edges at large values of hopping asymmetry. This difference decreases with decreasing hopping asymmetry.

Appendix E Additional AFQMC results for pairing amplitude

We provide below a survey of results from AFQMC for the pairing amplitude versus interaction strength and hopping asymmetry. In real-space we observe that the central peak of the pairing amplitude grows with interaction strength, but otherwise the qualitative behavior is relatively insensitive to interaction strength over the range of interaction strengths we consider. For smaller values of hopping asymmetry, the pairing amplitude is more spatially extended and isotropic, and still relatively insensitive to the interaction strength.