Research ArticlePHYSICS

Evidence of pair-density wave in spin-valley locked systems

See allHide authors and affiliations

Science Advances  29 Mar 2019:
Vol. 5, no. 3, eaat4698
DOI: 10.1126/sciadv.aat4698


Cooper pairs with a finite center-of-mass momentum form a remarkable state in which the superconducting order parameter is modulated periodically in space. Although intense interest in such a “pair-density wave” (PDW) state has emerged due to recent discoveries in high Tc superconductors, there is little theoretical understanding of the mechanism driving this exotic state. The challenge is that many competing states lie close in energy in seemingly simple models, such as the Hubbard model, in the strongly correlated regime. Here, we show that inversion symmetry breaking and the resulting spin-valley locking can promote PDWs over more commonly found spin stripes through frustration against magnetic order. Specifically, we find the first robust evidence for a PDW within density matrix renormalization group simulation of a simple fermionic model. Our results point to a tantalizing possibility in hole-doped group VI transition metal dichalcogenides, with spin-valley locked band structure and moderate correlations.


Recent experimental and theoretical developments have brought a renaissance to the idea of a modulated superconducting state that spontaneously breaks translational symmetry [see (1) and references therein and (29)]. Earlier efforts toward realization of modulated superconductors (10, 11) or toward an interpretation of associated experiments (12) have relied on generating finite-momentum pairing using spin-imbalance under an (effective) magnetic field, in close keeping with the original proposals [Fulde-Ferrell-Larkin-Ovchinnkov (FFLO)] (13, 14). Alternatively, momentum space split, spinless fermions in the context of doped Weyl semimetals (15) and hole-doped transition metal dichalcogenides (TMDs) (16) have been proposed as a platform for modulated superconductors due to pairing within Fermi pockets centered at finite crystal momentum. On the other hand, a modulated paired state proposed for cuprates requires a strong coupling mechanism (2). Such a strong coupling–driven state has been dubbed a pair-density wave (PDW) as a state distinct from FFLO-type superconductors.

The need for a strong coupling mechanism led to a search for the PDW state in numerical simulations. Numerous variational and mean-field studies have shown that PDW-type states are energetically competitive with uniform d-wave superconducting states in generalized t-J models, and it is thought that PDWs may become favorable in the presence of anisotropy (1721). Nevertheless, numerical evidence from the controlled approach of the density matrix renormalization group (DMRG) is lacking within simple fermionic models as the only evidence of a PDW within the DMRG was established in the one-dimensional (1D) Kondo-Heisenberg model (22). One signature difficulty in such a realization is that DMRG calculations on a Hubbard or t-J model on a square lattice with spin-rotation symmetry often find spin and charge stripe ground states instead of the PDW state (23). However, one could hope that frustrating spin order might nudge systems into a PDW state. Here, we turn to a Hubbard model on the frustrated triangular lattice with broken inversion symmetry that captures the hole-doped monolayer group IV TMDs.

Rapidly growing interest in the monolayer group VI TMDs has been fueled by the exotic possibilities driven by spin-orbit coupling (SOC) and lack of centrosymmetry (16, 2429), as well as superconductivity in the n-doped TMDs (2426, 3032). While the symmetry properties of the observed superconducting states remain unknown, the different translationally invariant superconducting channels for the TMDs have been previously classified in mean-field studies (3335). Recently Hsu et al. (16) used a weak-coupling renormalization group (RG) approach to investigate a repulsive interaction-driven pairing mechanism, predicting two topological superconducting instabilities with one of them being a spatially modulated intrapocket state. However, potentially strong correlation effects have largely been neglected despite the fact that the conduction electrons have substantial d-character. In this letter, we use DMRG calculations to study the effects of SOC on superconducting tendencies driven by repulsive interactions.


To capture the spin-valley locked Fermi surfaces that occur in the valence band of group VI TMDs (24, 36, 37) in a one-band model, we consider a nearest-neighbor tight-binding model on a triangular lattice with a staggered, spin-dependent magnetic flux of per plaquette (see Fig. 1A). This spin-dependent flux breaks the C6v symmetry down to C3v while preserving the time reversal symmetry, mimicking the Sz preserving SOC present in the Kane-Mele model and generating two distinct spin-polarized pockets for our Fermi surface. Furthermore, the flux introduces a small amount of anisotropy in the pockets (see Fig. 1B) analogous to that present in real materials such as MoS2 (38). Last, we include an on-site interaction. Hence, our model Hamiltonian isH=ijtij,σciσciσμiσciσciσ+Uinini(1)where tij is the spin-dependent complex nearest-neighbor hopping, μ is the chemical potential, and U is an on-site Hubbard interaction. The dispersion takes the form ofσ(k)=2i[Re(t)cos(δik)+Im(t)sin(δik)σz](2)where δi{x^,12x^+32y^,12x^32y^},σz=± for spin up and down, respectively; we define t such that t = ti+,i;↑, and the lattice spacing has been set to 1.

Fig. 1 Model and Fermi surface.

(A) The spin-dependent staggered flux pattern for one-spin component with ± Φ flux per plaquette. An opposite flux pattern for the other spin component guarantees time-reversal symmetry. The arrows indicate the direction of positive phase hopping. (B) Our Fermi surface with ti+x^,i;=23ei0.3π and μ = 4.6 in the tight-binding model in Eqs. 1 and 2. Here, the spin-valley locked, circular Fermi pockets are evident.

Within this model, we probe the superconducting response to a pair-edge field. We initially consider a moderately attractive interaction with U = −2, where uniform pairing in the A1-irrep is expected. Data for the following statements can be found in sections S1 and S3. First, we note that the phase disorder due to edge field effects quickly disappears upon moving into the bulk (see fig. S1). Inspecting the bond-singlet component of the superconducting order parameter Δijsinglet on the directed nearest-neighbor bonds away from the probe, we find a well-ordered uniform phase structure (fig. S3). Moreover, zooming into the phase structure within a single unit cell, we find the pair field expectation value to be isotropic and definitively s-wave. The uniform and isotropic nature of the order parameter phase distribution is a robust property of our results in the negative U regime, which is insensitive to the profile of the edge fields and occurs for all system sizes studied. The triplet channel behaves analogously, displaying homogeneous f-wave pairing (see fig. S3 in section S3).

Armed with the attractive, U < 0, result that can serve as a reference, we now study the moderately repulsive Hubbard regime with U = +2. Note that estimates for the band structure parameters of MoS2 have an SOC of 0.08 eV, a hopping strength of 0.5 to 1 eV, and an on-site repulsion of the 4d Mo orbitals of 2 to 10 eV (37, 39). Thus, with U/t ≈ 2 for our calculations, we lie at the lower end of these estimates for the effective interaction strength of these materials. Given the difficulty in estimating U/t for actual materials, TMDs represent a realistic material setting for our simulation, at least within the errors of these estimates. Earlier work using two-stage perturbative RG on a similar spin-valley locked model with repulsive U predicted superconductivity in the 2D E representation, where some linear combination of p- and d-wave symmetries occurs because of the lack of inversion symmetry (16). Unexpectedly, our DMRG simulation in this repulsive interaction regime reveals a tendency to break translational symmetry along the length of the cylinder. Specifically, the system forms a modulated paired state, where both the singlet and triplet bond pair order parameters are everywhere real with modulation in their signs (see Fig. 2). From the symmetry perspective, the observed state is analogous to the state proposed by (14). This modulation in the pair amplitude is evident in the plot of Δijsinglet for U = +2 in Fig. 2A, where an anisotropic phase structure within the unit cell is repeated with period 2. A similar unit cell doubling of the phase is seen in the triplet channel (see section S4). We find that this tendency to form a PDW is robust against changes in chemical potential, although the periodicity changes appropriately based on the filling (see section S5). For instance, increasing the chemical potential, μ, from μ = 4.6 to μ = 6.0 enlarges the unit cell by an additional lattice site (see fig. S7). Although there has been much interest in modulated superconducting states, to the best of our knowledge, this is the first report of a strong coupling–driven PDW within DMRG simulations of a simple fermionic model.

Fig. 2 Evidence of PDW oscillations.

(A) Arg (Δijsinglet) for all nearest-neighbors with U = +2 for our 3 by 36 lattice with periodic boundary conditions along the short direction and open boundary conditions along the long direction. For visibility, we truncate the plot so that only the third farthest from the edge field is shown. The line thickness is proportional to the pairing amplitude. (B) We plot the real and imaginary components of Δijsinglet and Δijtriplet for i,j along the middle rung of our lattice to present the phase oscillations.

To further understand this translational symmetry breaking, we examine the oscillations of the superconducting order parameter. Since Δijsinglet is characterized exclusively by π-phase shifts in the bulk, we may project it to real space after a gauge transformation and look at the decay properties of the pair field by plotting it for bonds directed along the middle rung (see Fig. 2B). For attractive interactions, the singlet pairing strength falls off gradually, as expected from the quasi-1D geometry of the system, and exhibits slowly varying oscillations (see fig. S4 in section S3) due to finite size effects induced by the open boundary conditions (40). On the other hand, the pairing profile for the U = +2 simulation shows an additional rapid oscillation centered about zero (see Fig. 2B). These oscillations centered about zero occur for all repulsive Hubbard simulations near U = +2. Thus, although the exact strength of the pairing response and the penetration depth of the edge field appear to have some dependence on the edge field profile and the length of the lattice, the PDW-type behavior reported has been observed for all system sizes and all edge-field types. This plot strongly resembles the plot of the same quantity in the Kondo-Heisenberg model with PDW (22).

Fourier transforming these oscillations, we find a single dominant wave vector that is approximately twice the Fermi radius of a single pocket, Q = 2kF, suggesting that the finite momentum of our Cooper pairs originates from the interplay between the pockets. This picture is reinforced by probing the effect of increasing the chemical potential (decreasing the pocket radius). Here, a PDW also develops, and we see a shift in wave vector consistent with our change in pocket size (see fig. S7). A schematic of the finite momentum pairing is presented in Fig. 3B. Since our PDW modulations are unidirectional and lie orthogonal to the applied edge field, our results are especially relevant to proximity-induced superconductivity in TMDs.

Fig. 3 Fourier decomposition of the PDW and bond charge order.

(A) Fourier transforms of the PDW and charge bond order. Zero momentum, i.e., constant contributions and decay effects have been removed. (B) Depiction of pairing in momentum space. The regions demarcated by dashed lines are the approximate pairing regions.

Last, we compare oscillations in the singlet pairing strength and in the bond charge density. As the pairing amplitude profile of our U = +2 simulation has net pair amplitude on the whole system due to the edge field and our edge field explicitly breaks gauge symmetry, charge modulation of the same period (Q = 2kF) is anticipated to be driven by the net component and the modulated pairing components (1). From a Ginzburg-Landau perspective, these add terms of the form ρQ(Δ0*Δ+Q+Δ0ΔQ*) and ρQ(Δ+Q+ΔQ*), respectively, to the free energy that account for the development of 2kF charge density wave oscillations as opposed to 4kF oscillations from ρ2QΔQ*Δ+Q, where ρQ and ΔQ here correspond to the density and PDW order parameters with momentum Q. This is clearly seen in the Fourier decomposition of the PDW and the bond charge order presented in Fig. 3A, where both orders are dominated by the same Fourier mode. We remark that while both the attractive and repulsive charge bond densities have oscillations that dip below their mean values, only the superconducting response of the repulsive case has oscillations centered about zero, suggesting that the PDW oscillations are not finite size effects.


In summary, we used the DMRG to study the superconducting tendencies of a repulsive U Hubbard model on a triangular lattice with spin-valley locking. These tendencies were probed by studying the pairing response profile in response to uniform and random pair fields along one edge. Our calculations indicate that the superconducting phase diagram of the model may be more complex than what was revealed by the previous perturbative RG study (16), with translational symmetry breaking superconducting states possibly in competition with a uniform state. The PDW observed breaks translational symmetry with the superconducting order parameter alternating sign. This atypical pairing response may be related to the fact that Ising SOC and the triangular lattice conspire to frustrate any spin order including spin stripe. We fail to find any appreciable spin response to be induced by an Sz-edge field coexistant with our pair field. Here, the moment dies off rapidly away from the edge, leaving no discernible moment and reaching practically zero (5 × 10−5 in units of the applied field) by site 4. Moreover, the fact that the apparent periodicity is not related to any Fermi surface properties suggests a purely strong coupling origin of the observed translational symmetry breaking. It will be interesting to study whether the observed PDW state can be found in a truly 2D setting using a different method such as the density matrix embedding theory, where a long-range order can be observed (41).


The DMRG is a powerful, nonperturbative method for studying strongly interacting systems (22, 4246). It has been used with great success to explore a diverse selection of strongly correlated phenomena highlighted by stripes, spin liquids, and superconductivity (22, 4351). However, since DMRG is quasi-1D in nature, no true long-range order can be seen in the correlations. Thus, to access our system’s superconducting tendencies, we implemented a pair-edge field motivated by the field-pinning approach underlying several earlier studies (22, 43, 52, 53). By biasing the system toward a particular superconducting state and studying the emergent symmetry of the appropriate order parameter in the bulk, one can gauge the model’s propensity for various instabilities.

Even in the absence of the spin-valley locking special to our model, the geometric frustration of the triangular lattice is known to foster exotic phases. While a consensus has emerged that the ground state of the Heisenberg model is a 120° Néel antiferromagnet (54), the Hubbard model has been shown to have tendencies toward spin liquid and chiral d + id superconducting states (55, 56). Note that within the context of the Heisenberg model, frustration can inhibit spin stripes.

We emphasize that due to the lack of inversion symmetry, even and odd pairing components can coexist (57). Thus the Sz preserving SOC allows for the mixing of Sz = 0 singlet and triplet states, i.e., the bond pair order parameter Δij=cicj should be an admixture of singlet and triplet componentsΔijsinglet=cicjcicjΔijtriplet=cicj+cicj(3)

Note that since our system lacks translational symmetry due to open boundary conditions, the pairing symmetry is not constrained to transform under a single irreducible representation. Nevertheless, the real-space structure of these bond-centered order parameters provides insight into the nature of the dominant pairing state.

We carried out our DMRG simulations on a cylinder with 3 unit cells in the periodic direction and 12, 18, 24, and 36 unit cells in the nonperiodic direction. The width is sufficiently large to sample both types of pockets in the Fermi surface but not so large as to make the DMRG prohibitively expensive for our available computational resources. We kept the band structure fixed and explored the effects of varying U. We investigated the superconducting susceptibility by applying a pair field along one edge, as illustrated above in Fig. 4. To reveal any inherent preferences for a particular superconducting channel, we considered two different phase structures for the edge field: a uniform field described by the A1-irrep and a random field, that isΔijedge=Vcicjeiϕij+h.c.(4)with the phase ϕij chosen to transform under the A1-irrep or be random for ij and 0 for i = j. We remark that all results presented in this paper have been shown to be independent of the phase structure of the edge field applied (see section S1). The strength of the pair field was fixed to be V = 0.1, about an order of magnitude less than the hopping amplitude, which is consistent with that used in previous studies (22, 52, 53). While the cylindricity of our geometry and the addition of a pair field break the C3v symmetry of the lattice and translational invariance, if there is a well-defined structure to order parameter in the bulk, then we expect to gain information about the inherent superconducting tendencies through the real-space structure of the order parameters in Eq. 3.

Fig. 4 Lattice and edge field.

A depiction of our lattice. It is periodic in the short direction with three unit cells and has open boundaries in the long direction. The ellipses on the right signify that multiple lengths are studied: L = 12, 18, 24, 36. The edge field, shown as red lines, is a pair field of the form given in Eq. 4. The nearest-neighbor hopping structure for spin up is also shown with the spin down hopping structure being the complex conjugate of that shown above.

Our DMRG calculation used the iTensor library developed by Stoudenmire and White (58). We performed up to 14 sweeps with a final bond dimension of M = 2500. This is sufficient to obtain energy convergence to O(10−7) for our repulsive calculations and O(10−10) for attractive calculations (see section S2). We focus exclusively on the interpocket instabilities of our model, where we may exploit the conservation of the Sz = 0 quantum number. Sz and fermion parity, N2, are conserved quantum numbers, but the U(1) particle number symmetry is broken by the pair field. As our starting point, we constructed an matrix product state that randomly samples the Sz = N2 = 0 sector of our Hilbert space picking 10 to 30 states from each even particle number sector with the number of states depending on the system size. In addition, we used exact diagonalization (ED) to ensure correctness in the DMRG simulations with the energetics from the two methods agreeing to within machine precision for small 3 by 3 systems where ED is computationally tractable.


Supplementary material for this article is available at

Section S1. Phase structure near the pair field

Section S2. DMRG convergence

Section S3. Phase structure for attractive interactions

Section S4. Triplet phase structure for repulsive interactions

Section S5. Effects of chemical potential on the PDW

Fig. S1. Singlet phase structure near the edge field.

Fig. S2. DMRG convergence.

Fig. S3. Phase structure for attractive interactions.

Fig. S4. Amplitude for attractive interactions.

Fig. S5. Triplet phase structure for repulsive interactions.

Fig. S6. Effect of chemical potential on the PDW phase structure.

Fig. S7. Effect of chemical potential on the PDW dominant Fourier mode.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank L. Balents, E. Berg, G. Chan, Y.-T. Hsu, S. Kivelson, K. T. Law, P. Lee, S. White, and H. Yao for the useful discussions. Funding: E.-A.K. acknowledges Simons Fellow in Theoretical Physics Award #392182 and DOE support under Award DE-SC0010313. E.-A.K. is grateful to the hospitality of the Kavali Institute of Theoretical Physics supported by NSF under grant no. NSF PHY-1125915, where this work was completed. J.V. acknowledges support by the NSF [Platform for the Accelerated Realization, Analysis, and Discovery of Interface Materials (PARADIM)] under Cooperative Agreement No. DMR-1539918 and in part by NSF DMR-1308089. Author contributions: J.V. performed the DMRG simulations. E.-A.K. was responsible for the initial conception of the problem. Both J.V. and E.-A.K. contributed to the interpretation of results and the writing of the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article