Displacement-Field-Driven Transition between Superconductivity and Valley Ferromagnetism in Transition Metal Dichalcogenides

TL;DR

This paper presents a theoretical mechanism for a displacement-field-controlled transition between superconductivity and valley ferromagnetism in 2D hexagonal systems, explaining recent experimental observations in twisted bilayer WSe2.

Contribution

It introduces a simple VHS-only model capturing the transition driven by displacement field D, predicting a switch from chiral superconductivity to valley ferromagnetism.

Findings

Weak D favors chiral d/p-wave superconductivity

Strong D induces valley ferromagnetic phase

Transition is controlled by inter-VHS interactions

Abstract

Recent experiments have observed transitions between superconductivity and correlated magnetism in twisted bilayer WSe$_2$ near van-Hove fillings, driven by the displacement field $D$. Motivated by the experiment, we theoretically propose a general mechanism for a $D$-controlled transition between superconductivity and ferromagnetism in two-dimensional (2D) spin-orbit-coupled hexagonal systems, where van Hove singularities (VHS) lie on the Fermi level. We show that such a transition can be naturally captured by a simple VHS-only model without Fermi surface details, where the inter-VHS interactions that govern the Fermi surface instabilities is controlled by $D$ through the band projection of screened Coulomb interaction. By treating this simple model with renormalization group technique beyond mean-field level, we find that a chiral $d/p$-wave superconductivity naturally dominates under a weak displacement field $D<D_c$. At a stronger displacement field $D>D_c$, a \textit{valley ferromagnetic phase} (vFM) takes over, which is spatially non-uniform due to valley-modulated magnetization. Finally, we discuss generic conditions for the predicted superconductivity-to-ferromagnetism transition to take place in the rich family of few-layer hexagonal van der Waals material systems. Taking twisted bilayer WSe$_2$ as a case study, we discuss experimental detections that can falsify our prediction.

Equations

$$h_{\uparrow }(\text{r})=\left(\begin{array}{cc}-\frac{{\hslash }^2(\text{k}-{\mathbf{\kappa }}_t{)}^2}{2m^{\ast }}+{\Delta }_t(\text{r})& {\Delta }_{\text{T}}(\text{r})\\ {\Delta }_{\text{T}}^{†}(\text{r})& -\frac{{\hslash }^2(\text{k}-{\mathbf{\kappa }}_b{)}^2}{2m^{\ast }}+{\Delta }_b(\text{r})\end{array}\right),$$

$$H_{\text{int}}=\frac{1}{2A}\sum _{l{l}^{\prime }=t,b}\sum _{\text{k}}{V}_{l{l}^{\prime }}(k){\rho }_l(\text{k}){\rho }_{l^{\prime }}(-\text{k}),$$

$${\rho }_{l,\tau }(\text{q})=\sum _{s=\uparrow ,\downarrow }\sum _{\text{p}\in \text{mBZ}}{\Lambda }_{sl\tau }(\text{p},\text{p}+\text{q})c_{s\tau }^{†}(\text{p})c_{s\tau }(\text{p}+\text{q}),$$

$$\displaystyle H_{\text{patch}}=\sum_{n=1}^{3}\Big{[}\tilde{g}_{2}c^{\dagger}_{n}c_{n}c^{\dagger}_{\bar{n}}c_{\bar{n}}+\sum_{m\neq n}\Big{(}\tilde{g}_{3}c^{\dagger}_{m}c_{n}c^{\dagger}_{\bar{m}}c_{\bar{n}}$$

$$\displaystyle\quad+\frac{1}{2}\tilde{g}_{6}(c^{\dagger}_{n}c_{n}c^{\dagger}_{{m}}c_{{m}}+c^{\dagger}_{\bar{n}}c_{\bar{n}}c^{\dagger}_{\bar{m}}c_{\bar{m}})+\tilde{g}_{6}^{\prime}c^{\dagger}_{n}c_{n}c^{\dagger}_{\bar{m}}c_{\bar{m}}\Big{)}\Big{]},\;\;$$

$${\tilde{g}}_2={\tilde{g}}_6^{\prime }=V_{ll}(0)-{\Lambda }_1(D)(V_{ll}(0)-V_{l\bar{l}}(0)),$$

$${\tilde{g}}_3=V_{ll}(Q){\Lambda }_2(D)+V_{l\bar{l}}(Q){\Lambda }_3(D),$$

$${\tilde{g}}_6={\tilde{g}}_2-{\tilde{g}}_3,$$

$${\Lambda }_1(D)=|{\Lambda }_{\uparrow ,t}({\text{P}}_1,{\text{P}}_1){\Lambda }_{\uparrow ,b}({\text{P}}_1,{\text{P}}_1)|,$$

$${\Lambda }_2(D)=|{\Lambda }_{\uparrow ,t}({\text{P}}_1,{\text{P}}_2){|}^2+|{\Lambda }_{\uparrow ,b}({\text{P}}_1,{\text{P}}_2){|}^2,$$

$${\Lambda }_3(D)=2|{\Lambda }_{\uparrow ,t}({\text{P}}_1,{\text{P}}_2){\Lambda }_{\uparrow ,b}({\text{P}}_1,{\text{P}}_2)|.$$

$$\frac{d{g}_i}{dy}=\sum _{lm}{d}_{lm}{g}_l{g}_m,i,l,m=2,3,6,6^{\prime }.$$

$${\beta }_{d/p\text{-SC}}=G_2-G_3,$$

$${\beta }_{\text{vFM}}=-d(y)(G_2-2G_6+2G_6^{\prime }),$$

$${\int }_{{\text{p}}^{\prime }\in \text{FS}}\frac{d{\hat{\text{p}}}^{\prime }}{|{\text{v}}_{{\text{p}}^{\prime }}|}\Gamma (\text{p},{\text{p}}^{\prime })\Delta ({\text{p}}^{\prime })=\lambda \Delta (\text{p}),$$

$$\displaystyle\delta H_{\text{vFM}}=\Delta_{\text{vFM}}\Big{(}\sum_{\textbf{p}\in K}c_{\uparrow}^{\dagger}(\textbf{p})c_{\uparrow}(\textbf{p})-\sum_{\textbf{p}\in K^{\prime}}c_{\downarrow}^{\dagger}(\textbf{p})c_{\downarrow}(\textbf{p})\Big{)}.\;\;$$

Full text

Introduction— Two-dimensional superconductivity (SC) and its transitions into other symmetry-broken phases have recently been observed in a variety of van der Waals materials, tunable by external fields [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], carrier filling [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22], and other experimental parameters [23, 24]. Besides graphene-based systems, robust SC has recently been demonstrated in moiré transition metal dichalcogenides (TMDs) with strong spin-valley locking [25]. In particular, SC was seen in twisted bilayer WSe2 (tWSe2) at different twist angles, where transitions were observed from superconducting to other correlated phases upon increasing an external displacement field [26, 27]. The natures of the superconducting state [28, 29, 30, 31, 32, 33, 34, 35, 36, 37] and the nearby correlated phases [38, 39, 40] have been extensively studied, where the latter were widely believed to exhibit magnetism. For instance, recent magnetic dichroism measurement has found a ferromagnetic phase close to the Van Hove (VH) filling in the $\theta=3.2^{\circ}$ system [41], along with theories suggesting close competition among ferro- and antiferromagnetic phases [42]. On the other hand, antiferromagnetic orders were proposed by mean-field [43] and functional renormalization group (RG) [44] theories for the $\theta=5^{\circ}$ system, while definitive experimental evidence for the magnetic order remains to be reported. The rapid experimental development and the rich variety of few-layer TMD systems thus call for generic mechanisms for SC-to-magnetism transitions that do not rely on system details.

In this work, we propose a mechanism for a displacement-field $D$-driven transition between SC and ferromagnetism in two-dimensional spin-orbit-coupled hexagonal systems along the Van Hove filling curve $n_{V}(D)$ (see Fig. 1a). In general, such a system can exhibit three VH singularities per valley on the Fermi level with appropriate gating (see Fig. 1b,c), where these six VH momenta evolve on the Fermi surface (FS) with $D$. For instance, in $\theta=5^{\circ}$ tWSe2, the theoretically determined VH filling $n_{V}(D)$ was found to trace through the observed SC phase and the nearby magnetic phase, which agrees well with the Hall resistivity measurement [27]. Due to the diverging density of states at the VH momenta, a SC-to-magnetism transition in such spin-orbit-coupled hexagonal systems, if exists, should be driven by the evolution of these VHS under $D$, regardless of the Fermi surface details.

Motivated by this hypothesis, within a simple model of six-VHS, we examine how the symmetry-allowed inter-VHS interactions receive a displacement field $D$-dependence from screened Coulomb interactions and renormalize in the long-wavelength limit. We show that among all possible instabilities, the renormalized inter-VHS interactions drive a mixed $p$- and $d$-wave chiral SC at $D<D_{c}$ and spatially non-uniform valley ferromagnetism at $D>D_{c}$, where the displacement field $D$ triggers a Stoner-like transition. We expect that such a transition can generically occur in the family of few-layer TMD systems and other two-dimensional (2D) hexagonal spin-orbit-coupled systems near the VH filling, in the weak to intermediate coupling regime.

Microscopic model— We start from a well-known continuum model for general twisted homobilayer TMDs, where the topmost valence bands in the Moire Brillouin zones (mBZ) at valley $K$ and $K$’ are spin-up and -down, respectively, due to the Ising spin-orbit coupling (see Fig. 1b). The non-interacting model reads $H_{0}=\sum_{\sigma=\uparrow,\downarrow}\int d^{2}\textbf{r}\psi_{\sigma}^{\dagger}(\textbf{r})h_{\sigma}(\textbf{r})\psi_{\sigma}(\textbf{r})$ [46], where $\psi_{\sigma}^{\dagger}(\textbf{r})$ creates an electron with spin $\sigma$ at spatial position r, and the spin-up Hamiltonian at valley $K$ is given by

[TABLE]

where $h_{\downarrow}(\textbf{r})$ at valley $K^{\prime}$ is related to Eq. 1 by time-reversal symmetry. In Eq. 1, k is the in-plane crystalline momentum measured from $\Gamma$-point in Fig. 1b, $l=t,b$ denotes the top and bottom layer, while $\Delta_{t,b}(\textbf{r})=\pm\frac{eV_{z}}{2}+2v\sum_{j=1,3,5}\cos(\textbf{g}_{j}\cdot\textbf{r}\pm\psi)$ and $\Delta_{\text{T}}(\textbf{r})=w(1+e^{-i\textbf{g}_{2}\cdot\textbf{r}}+e^{-i\textbf{g}_{3}\cdot\textbf{r}})$ describe the layer-dependent potentials and inter-layer tunneling, respectively. The potential difference, $V_{z}$ is directly proportional to the displacement field $D$ through $V_{z}=D\times d/\epsilon$ where $d$ is the layer separation and $\epsilon$ is the dielectric constant. Here, $\textbf{g}_{j=1,2,3}$ are the reciprocal lattice vectors of mBZ that can be obtained by rotating $\textbf{g}_{1}=\frac{4\pi\theta}{\sqrt{3}a_{0}}\hat{x}$ counterclockwise by an angle of $(j-1)\pi/3$, $j=2,3$. The resulting electronic dispersion generally exhibits three VHS at momenta $\textbf{P}_{1,2,3}$ between the mBZ corners $\bm{\kappa}_{t,b}$ in valley $K$, and another three at $\textbf{P}_{\bar{1},\bar{2},\bar{3}}=-\textbf{P}_{1,2,3}$ in valley $K^{\prime}$ [28] (see Fig. 1c). When the electronic density $n$ is tuned to the Van Hove filling $n_{V}(D)$, at which the VH points sit on the FS, the density of state (DOS) becomes logarithmically divergent so that the FS becomes susceptible to a wide variety of instabilities under even weak-to-intermediate interactions.

We consider intra- and inter-layer screened Coulomb interactions projected onto the topmost band

[TABLE]

where $A$ is the system area. In Eq. 2, $V_{ll}(k)$ and $V_{l\bar{l}}(k)$ are the intra- and inter-layer Coulomb potentials at momentum $\textbf{k}=\textbf{q}+\textbf{G}$, respectively, where q lies in the first mBZ, and G is the moire reciprocal lattice vector 111By solving the Poisson equation with double-gate boundary conditions, we obtain the $k$-dependence of $V_{ll^{\prime}}(k)$ in Fig. 1d (see Supplementary Material (SM) Section I [60]). Since $V_{ll^{\prime}}(k)$ decays fast with $k$, we will neglect the $G\neq 0$ components for simplicity.. Next, the density operator for layer $l$ can be expressed as $\rho_{l}=\sum_{\tau}\rho_{l,\tau}$ in terms of the band $\tau$-resolved density

[TABLE]

where we will project $\rho_{l,\tau}$ onto the topmost band and omit the band label $\tau$ in the following. Here, $c_{s}(\textbf{p})$ annihilates spin-$s$ electrons at momentum p, and $\Lambda_{sl}(\textbf{p},\textbf{p}+\textbf{q})=\langle u_{sl}(\textbf{p})|u_{sl}(\textbf{p}+\textbf{q})\rangle$ is given by the overlap between Bloch wavefunctions, which satisfies the normalization condition $\sum_{l=t,b}\Lambda_{sl}(\textbf{p},\textbf{p})=1$ for all p and $s$. Note that although the Coulomb potentials $V_{ll^{\prime}}$ does not depend on the displacement field $D$, the interaction $H_{\text{int}}$ in Eq. 2 still receives an explicit $D$-dependence through the wavefunction overlap $\Lambda_{sl}(D)$, leading to a $D$-driven phase transition.

Patch model — To investigate the FS instabilities at VH-fillings $n_{V}(D)$ at different displacement field strengths $D$, we first simplify the model $H=H_{0}+H_{\text{int}}$ by restricting the momentum k to only the six patches centered at VH points with sizes $k_{\Lambda}\ll G_{0}$ (see Fig. 1c), rather than the full mBZs. This patch approximation is justified by the diverging DOS at the VH points. For the low-energy electrons lying within the patches, the screened Coulomb interactions $H_{\text{int}}$ in Eq. 2 now simplifies to only four inequivalent inter-patch interactions due to the spin-valley locking [28]

[TABLE]

where $c^{\dagger}_{n}$ creates spin-up electrons in patch $n=1,2,3$ at valley $K$, and $c^{\dagger}_{\bar{n}}$ creates spin-down electrons in patch $\bar{n}$ at valley $K^{\prime}$. The strengths of the intra-valley (inter-valley) density-density interactions $\tilde{g}_{6}$ ($\tilde{g}_{2}$ and $\tilde{g}_{6}^{\prime}$) as well as the inter-valley scattering $\tilde{g}_{3}$ in Eq. 4 are related to the Coulomb potentials in Eq. 2 as

[TABLE]

where $Q\equiv Q_{nm}=|\textbf{P}_{n}-\textbf{P}_{m}|$ for all patches $n$ and $m$, and

[TABLE]

Using Eqs. 5 and 6, we numerically calculate the $D$-dependence of the bare interactions $\tilde{g}_{i}$ in Fig. 2a. In particular, since $\Lambda_{1(2,3)}$ are decreasing (increasing) functions of $D$ and ${V}_{ll}(0)-{V}_{l\bar{l}}(0)>0$, we find that $\tilde{g}_{6}^{\prime}$ and $\tilde{g}_{2}$ decrease with $D$ while $\tilde{g}_{3}$ and $\tilde{g}_{6}$ increase with $D$ (see Fig. 2a).

Renormalization group analysis — The inter-patch interactions in Eq. 4 can lead to a wide variety of FS instabilities even in the weak-to-intermediate-coupling regime since the DOS diverges at VH filling. To identify the dominant FS instability at different displacement field strengths $D$, we perform a renormalization group (RG) analysis, which unbiasly considers the full set of instabilities and identifies the one that is driven by the most relevant interaction in the low-energy limit [48, 49, 50, 51, 52]. This approach was widely applied to various few-layer Van der Waals systems in the weak-to-intermediate coupling regime [53, 54, 55, 56, 57, 58, 54, 59]. In particular, the framework for general 2D hexagonal systems without spin degeneracy was proposed in Ref. [28], which consists of three steps. We will briefly review this established method in the following, where details are included in the SM section II for completeness [60].

In Step 1, we examine the evolution of the inter-patch coupling constants $g_{i}(y)$ as a function of an inverse energy scale $y\equiv\nu_{0}\text{ln}^{2}(\Lambda/E)/2$ 222The prefactor $\nu_{0}$ depends on the specific dispersion obtained from Eq. 1. For the RG calculation, we take the phenomenological value $\nu_{0}=\frac{1}{2\pi W}$ for the bandwidth $W$., which flows from a higher energy scale $\Lambda$ at $y=0$, corresponding to the patch size $k_{\Lambda}$, towards the low-energy limit at $y\rightarrow\infty$. Such evolution is governed by one-loop RG equations with a general form of

[TABLE]

Here, $g_{i}\equiv\nu_{0}\tilde{g}_{i}$, and the factor $d_{lm}=\frac{1}{\nu_{0}}\frac{d}{dy}\Pi(\textbf{q})$ captures the evolution of the non-interacting susceptibilities $\Pi(\textbf{q})$ at momenta $\textbf{q}=0,\textbf{Q}_{nm}$, which physically correspond to the density of states or the particle-hole and particle-particle nesting degrees between patches $n$ and $m$, respectively. This set of RG equations is capable of capturing weak-to-intermediate coupling physics in systems with no spin-degeneracy and no perfectly nested FS by including the subleading ln-divergent terms on top of the leading ln2 terms [28].

In Step 2, we identify the interaction $\Gamma_{j}$ that supports each FS instability $j$. We consider all symmetry-allowed candidate order parameters $\Delta_{j}$, including zero-momentum superconductivity $\Delta_{j}(\textbf{p})=\langle c^{\dagger}_{\uparrow}(\textbf{p})c^{\dagger}_{\downarrow}(-\textbf{p})\rangle$, pair-density waves $\langle c^{\dagger}_{\uparrow}(\textbf{p})c^{\dagger}_{\uparrow}(-\textbf{p}+\textbf{Q}_{nm})\rangle$, $\langle c^{\dagger}_{\uparrow}(\textbf{p})c^{\dagger}_{\downarrow}(-\textbf{p}+\textbf{Q}_{n\bar{m}})\rangle$, density waves $\langle c^{\dagger}_{s}(\textbf{p})\sigma^{a}_{ss^{\prime}}c_{s^{\prime}}(\textbf{p}+\textbf{q})\rangle$, as well as uniform spin and charge orders $\langle c^{\dagger}_{s}(\textbf{p})\sigma^{a}_{ss^{\prime}}c_{s^{\prime}}(\textbf{p})\rangle$, where $\sigma^{a}$ with $a=0,1,2,3$ are the spin Pauli matrices. The corresponding interaction $\Gamma_{j}(y)=\sum_{ji}\eta_{ji}g_{i}(y)$ that drives each instability $j$ is given by a certain linear combination of the inter-patch interactions $g_{i}$, which can be systematically identified diagrammatically (see SM Section II [60]).

In Step 3, we identify the most relevant FS instability $j$ by identifying the strongest tendency $\beta_{j}\equiv d_{j}(y)\Gamma_{j}(y)$ as the inverse energy $y$ flows towards a critical low-energy limit $y_{c}$, quantified by the product of the driving interaction $\Gamma_{j}$ and the factor $d_{j}(y)=\frac{1}{\nu_{0}}\frac{d}{dy}\Pi(\textbf{q})$ that captures the ratio between the non-interacting susceptibilities of instability $j$ and superconductivity. This tendency determines the susceptibility $\chi_{j}\sim(y_{c}-y)^{2\beta_{j}+1}$ of instability $j$, where the dominant instability has the most negative exponent $2\beta_{j}+1<0$.

*Displacement-field-driven phase transition — * We now perform the RG analysis on the patch model $H_{\text{patch}}$ to show how a displacement field $D$ can drive a superconductivity-magnetism transition in generic 2D hexagonal systems at Van Hove filling. Within the framework of patch RG analysis, the dependence of $D$ explicitly enters the initial conditions $\tilde{g}_{i}\equiv g_{i}(y=0)/\nu_{0}$ of the RG differential equations in Step 1. As shown in Fig. 2a and b, the initial conditions $\tilde{g}_{2,3,6,6^{\prime}}$ exhibit different field-dependence through the Bloch-function overlaps at $y=0$ (see Eqs. 5, 6) so that the signs and relevance of $g_{i}(y)$ in the low-energy limit $y\rightarrow\infty$, and thus the instability tendency $\beta_{j}$, can qualitatively alter as $D$ increases. A change in the dominant instability $j$ can thus occur when $D$ reaches some critical value $D_{c}$, hence a displacement-field driven phase transition.

By numerically calculating the tendencies $\beta_{j}$ for all instabilities $j$ considered in Step 2, we find that the two leading instabilities are the $d/p$-wave superconductivity ($d/p$-SC) and the valley ferromagnetism (vFM) when the parameters for tWSe2 are used. Specifically, their corresponding tendencies $\beta_{j}$ are given by [28]

[TABLE]

where $d(y)=\frac{d}{dy}\Pi_{\text{ph}}(0)$. In Fig. 2b, we show the competition between these two leading instabilities quantified by $2\beta_{j}+1$, where we find that the $d/p$-SC phase dominates at a weaker field $D<D_{c}$, whereas the vFM phase is favored at a stronger field $D>D_{c}$ through a Stoner-like transition.

We now discuss how the $D$-dependence in the initial conditions at a higher-energy $y=0$ in Fig. 2a drives the transition from $d/p$-SC to vFM in the low-energy limit in Fig. 2b. Under a weak field $D<D_{c}$ where SC dominates, the sign and relevance of the inter-patch interactions $g_{i}$ are shown in the RG flows in Fig. 2c. The density-density interaction $g_{2}$ between opposite patches $\textbf{P}_{n}$ and $\textbf{P}_{\bar{n}}$ is the key interaction that drives the Kohn-Luttinger-like SC. Although starting out as a repulsion at the higher energy scale $y=0$, $g_{2}$ becomes attractive due to the particle-hole fluctuations from other interactions in the low-energy limit. Moreover, we find that the unconventional $d/p$-wave pairing dominates over the conventional pairing driven by $G_{2}+2G_{3}$ (See SM Section II [60]) due to the repulsively relevant inter-patch scattering $g_{3}$, which favors momentum-dependent order parameters.

As the displacement field grows stronger to $D>D_{c}$, the typical RG flows of $g_{i}(y)$ qualitatively change from Fig. 2c to Fig. 2d due to the seemingly mild changes in the initial condition, which now favor vFM over superconductivity. There are two key changes in the RG flows that lead to the phase transition: (1) An inter-valley density density interaction $g_{2}$ flips from a relevant attraction that drives SC into a repulsion that drives FM instead. Second, the intra-valley density-density interaction $g_{6}$ flips from a relevant repulsion that hurts vFM into an attraction that supports vFM (see Eq. 8). These two changes collectively drive a Stoner-like transition when the displacement field $D$ becomes stronger.

Topological superconductivity at $D<D_{c}$ — Our RG results at a weaker field finds that the superconductivity favors two degenerate nodal $d/p$-wave pairing gaps (see Fig. 3a), which form a two-dimensional irreducible representation (irrep) in the relevant point group $C_{3v}$. Due to energetic considerations [36, 52], these two degenerate gaps are known to prefer a chiral linear combination so that the FS can be fully gapped (See SM Section II [60]). We therefore expect a chiral $d/p$-wave SC that is topological with Chern number 2.

Although our RG result was obtained at VH filling, the finding of chiral $d/p$-wave pairing gap is robust even away from the VH filling in spin-orbit coupled hexagonal systems, assuming that the screened Coulomb interaction is stronger than electron-phonon coupling. This is because such systems typically has a point group of $C_{3v}$, where the inversion symmetry is broken by the spin-orbit coupling so that the even- and odd-parity gaps can mix. There are thus only two possible fully gapped superconducting states described by the irreps of $C_{3v}$: the fully gapped $s/f$-wave gap in the trivial irrep A, and our case of chiral $d/p$-wave gap in the two-dimensional irrep E. Due to the $s$-wave component in the former possibility, screened Coulomb interaction generally favors the chiral $d/p$-wave gap.

To further demonstrate the stability of the chiral $d/p$-wave SC away from the Van Hove filling $n_{V}(D)$, in Fig. 3b,c we calculate the gap function $\Delta(\textbf{p})$ on the full FS and the critical temperature $T_{c}(n)/T_{c}(n_{V})$ along the black dotted line in Fig. 1a by solving the linearized gap equation self-consistently

[TABLE]

where $\textbf{v}_{\textbf{p}}$ is the Fermi velocity, the interaction vertex $\Gamma(\textbf{p},\textbf{p}^{\prime})$ is set to be the dressed Hubbard interaction up to one-loop correction [62], and the effective pairing interaction strength $\lambda$ is obtained by solving for the eigenvalues of $\Gamma(\textbf{p},\textbf{p}^{\prime})$, which determines the critical temperature $T_{c}$. In Fig. 3b, we find that the energetically favored gap functions remain in the $E$ irrep when moving away from the VH filling, where we show the two degenerate $d/p$-wave gap functions on the FS at the representative point labeled in Fig. 1a. Moreover, we show that this $d/p$-wave pairing gap remains stable away from the VH filling, as signaled by the smallest effective interaction $\lambda_{\text{min}}<0$ being negative (see Fig. 3c). This effective pairing interaction determines the critical temperature as $T_{c}\propto e^{-1/|\lambda_{\text{min}}|}$, which peaks at the VH filling $n_{V}$ as expected from the diverging DOS.

Valley ferromagnetism at $D>D_{c}$— The vFM phase we find at a stronger field $D>D_{c}$ is characterized by the order parameter (see Fig. 3d)

[TABLE]

Due to the strong Ising SOC, this ferromagnetic order parameter pushes the topmost bands at valley $K$ and $K^{\prime}$ up and down, respectively (see Fig. 3e), resulting in a vFM state. Intuitively, this order parameter acts as a time-reversal-breaking pseudo-magnetic field that enhances (suppresses) the Ising SOC strength at valley $K$ ($K^{\prime}$), leading to a finite magnetization that is spatially non-uniform (see Fig. 3f).

We now comment on the difference and competition between the vFM state we find and antiferromagnetic (AFM) states that were often proposed in prior works for twisted bilayer TMDs [43, 44]. Although both vFM and AFM order parameters are spatially non-uniform, vFM has a finite magnetization while AFM has zero magnetization. In the vFM state, the magnetization originates from the imbalance between spin densities between the spin-up density at valley $K$ and the spin-down density at valley $K^{\prime}$, and thus spatially modulates at a wavelength $\propto\pm 1/(2K)$. This vFM instability can naturally occur when a spin ferromagnetic order develops in an Ising spin-orbit-coupled metal with multiple valleys. The competition between the vFM and AFM instabilities can be quantified by the tendency $\beta_{j}\sim\Pi_{j}\Gamma_{j}$, where $\Pi_{j}$ and $\Gamma_{j}$ are the corresponding non-interacting susceptibility and the driving interaction of instability $j$ at the critical inverse energy scale $y_{c}$, respectively [63]. As shown in Fig. 4, although the non-interacting susceptibility of the AFM state $\Pi_{j}=\Pi_{\text{ph}}(\textbf{Q}_{n\bar{m}})$ is slightly larger than that of the vFM state $\Pi_{\text{ph}}(0)$ at a larger displacement field $D>D_{c}$, the driving interaction of the vFM state $\Gamma_{\text{vFM}}$ is much larger than the AFM interaction $\Gamma_{\text{AFM}}$ so that the system tilts towards the vFM state in this close competition.

Experimental detections and material systems — The proposed vFM state can be experimentally distinguished from an AFM state by detecting its finite magnetization, e.g. using magnetic dichroism, as was done in $2^{\circ}$ to $3^{\circ}$ tWSe2 [41]. Furthermore, the valley-dependent band structures in Fig. 3e could be measured by angle-resolved photoemission spectroscopy (ARPES), whereas the non-uniform spin density in Fig. 3f could be detected by spin-dependent local probes, such as spin-dependent scanning tunneling microscope or scanning superconducting quantum interference device (SQUID). Finally, with transport measurements, ferromagnetism and AFM can be distinguished from their Curie–Weiss behaviors in temperature dependence [41, 9]. The resistivity of the vFM state depends on the filling due to the relation between the reconstructed low-energy band structure and chemical potential, which could reach an insulating state when the vFM gap is comparable to the bandwidth. As shown in Fig. 3e, the vFM state becomes strictly insulating at filling $n=1\sim n_{V}(D=40\text{meV})$, and the longitudinal resistivity may slowly decrease as moving away from $n=1$.

The nature of the superconducting state can be experimentally examined by several ways. First of all, the hypothesis that the superconductivity is driven by multiple VHS on the FS can be validated by an optical absorption measurement in the superconducting state, signaled by a pronounced peak at frequency $\omega\sim 2\Delta$ exceeding the disorder background [64]. The anisotropy in the current-induced optical conductivity Re$\sigma_{ii}(\omega)$ between $i=x$ and $y$ can further determine whether these VHS are conventional and higher-order, with logarithmic and power-law diverging DOS, respectively [64]. For the pairing symmetry, the conventional $s/f$-wave and unconventional $d/p$-wave order parameters can be distinguished with phase-sensitive measurements [65].

Finally, we comment on the material platforms that could realize the proposed phase transition as well as the valley ferromagnetism and chiral $d/p$-wave superconductivity on the two sides of the transition. One plausible system that exhibits a similar displacement-field-driven superconducting to magnetic phase transition is the $5^{\circ}$-tWSe2 [27], from which we adapt the material parameters for our numerical calculations. We propose that the observed SC to magnetism transition in Ref. [27] can be a chiral $d/p$-wave SC-to-vFM transition mediated by inter-VH interactions, if the FS away from VH points does not exhibit perfect nesting, as we find using the continuum model $H_{0}$ [46]. While previous theoretical works have proposed magnetic states with large momentum transfer, such as 120-degree antiferromagnetism [43] and inter-valley coherent states [44], definitive experimental evidence for the magnetic order parameter remains elusive. We thus urge near-future detections to examine our proposal.

Besides the displacement field driven transitions, we expect that a similar Stoner-like transition into the vFM phase could also be driven by other experimental knobs that can effectively tune the interaction strengths among VH points near the VH filling, such as the twist angle. Recent experimental studies in different smaller-angle tWSe2 systems have reported ferromagnetism [41, 66], offering a promising platform to test our proposed vFM state via the suggested measurements. Finally, our results highlight that such a Stoner-like transition into a vFM state and a nearby chiral $d/p$-wave SC are generic weak-to-intermediate coupling phenomenon that can occur in other twisted or untwisted spin-orbit-coupled few-layer systems near the VH filling away from perfect nesting.

Acknowledgment— The authors are grateful for useful discussions with Fengcheng Wu, Ming Xie, and Daniel Kaplan. H.-J. Y. and Y.-T. H. acknowledge the support from National Science Foundation Grant No. DMR-2238748. H.-J.Y. also acknowledges the support from the Society of Science Fellows Postdoctoral Program in the College of Science. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. This research was also supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP).