## Abstract

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 *T*_{c} 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.

## INTRODUCTION

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 (*2*–*9*)]. 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 (*17*–*21*). 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*, *24*–*29*), as well as superconductivity in the n-doped TMDs (*24*–*26*, *30*–*32*). 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 (*33*–*35*). 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.

## RESULTS

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 C_{6v} symmetry down to C_{3v} while preserving the time reversal symmetry, mimicking the *S*_{z} 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 MoS_{2} (*38*). Last, we include an on-site interaction. Hence, our model Hamiltonian is*t*_{ij,σ} 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*t* such that *t = t*_{i+x̂,i;↑}, and the lattice spacing has been set to 1.

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 *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 MoS_{2} 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 *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.

To further understand this translational symmetry breaking, we examine the oscillations of the superconducting order parameter. Since *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* = 2*k*_{F}, 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.

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* = 2*k*_{F}) 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 *k*_{F} charge density wave oscillations as opposed to 4*k*_{F} oscillations from _{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.

## DISCUSSION

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 *S*_{z}-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*).

## MATERIALS AND METHODS

The DMRG is a powerful, nonperturbative method for studying strongly interacting systems (*22*, *42*–*46*). It has been used with great success to explore a diverse selection of strongly correlated phenomena highlighted by stripes, spin liquids, and superconductivity (*22*, *43*–*51*). 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 *S*_{z} preserving SOC allows for the mixing of *S*_{z} = 0 singlet and triplet states, i.e., the bond pair order parameter

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_{ij} chosen to transform under the A1-irrep or be random for *i* ≠ *j* 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 C_{3v} 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.

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 *S*_{z} = 0 quantum number. *S*_{z} and fermion parity, N_{2}, 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 *S*_{z} = N_{2} = 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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/3/eaat4698/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).