# Supplementary Information

## Supplementary Information

###### Contents

  1. I Extended data
    1. I.1 Spectra of all devices
    2. I.2 The in-gap features of the monolayer device
    3. I.3 Temperature dependence of the spectrum
    4. I.4 Parallel Magnetic Field
    5. I.5 Metallic H-TMD Sequence with Additional Materials
  2. II Theoretical Computation of the upper critical field
    1. II.1 Tight-binding model
      1. II.1.1 Without CDW
      2. II.1.2 With the CDW
    2. II.2 Superconductivity – BCS theory
      1. II.2.1 Analytical solution in the case $U_{1}=U_{2}=U_{3}$
      2. II.2.2 Self-consistent determiniation of $T_{c}$
    3. II.3 Superconductivity – Ginzburg-Landau theory
      1. II.3.1 Quadratic terms:
      2. II.3.2 Quartic term:
    4. II.4 Details of computations

###  I Extended data

####  I.1 Spectra of all devices

Figure S1: Tunneling spectra of all devices at base temperature and zero
field.

Fig. S1 shows the base temperature, zero-field tunneling spectra of all TaS2
devices of thicknesses between one to 32 layers. The horizontal axes are all
bias voltage (in units of mV) and the vertical axes are the differential
conductance dI/dV, normalized by the value in the normal regime.

####  I.2 The in-gap features of the monolayer device

The singular measured monolayer devices exhibits clear in-gap features that
disappear at a temperature 0.12$T_{C}\approx$ 0.4 K (see Fig. S2a),
reminiscent of a smaller, second gap, such as the Se-derived gap observed in
NbSe2 Dvir  _et al._ (2018). However a chalcogen-derived second order
parameter cannot be expected in TaS2, since the chalcogen band is not
populated. The identification of the sub-gap feature with a second gap is also
not supported by the fact that up to the highest field measured (8.5T), this
feature closely follows the behavior of the main gap (Fig. S2b), in contrast
to the effect seen in NbSe2 where the two gaps respond qualitatively
differently to magnetic field due to the distinct dimensionality and different
diffusion parameters.

Figure S2: The spectrum of a TaS2 monolayer junction. (a) The spectrum at base
temperature and at T=0.4K, equivalent to $\approx$0.12$T_{C}$, where the inner
feature becomes undetectable. (b) The derivative of the spectrum d2I/dV2 at
base temperature, as a function of in-plane field (only the positive bias is
shown for clarity). The dashed lines indicate constant voltage, highlighting
the main peak and the inner feature.

####  I.3 Temperature dependence of the spectrum

We measure the spectrum of each junction as a function of temperature, ranging
from the cryostat base temperature up to slightly above the temperature at
which the spectrum becomes flat, indicating normal metallic behavior. The
simultaneous effects of the temperature dependence of $\Delta$ and the
smearing of the spectrum due to thermal fluctuations of the order $k_{B}T$
result in filling of the gap and reduction of the height of the quasiparticle
peaks, effectively smearing the spectrum until it becomes metallic. The
spectra of two junctions at various temperatures are shown in Fig. S3. We note
a small mismatch between the continuously measured zero-bias conductance
$G_{0}$ during cooldowns and the values extracted from the full spectra (see
Fig. 2a of the main text). This mismatch likely originates from the small
temperature gradient between the sample and its temperature sensor, expressed
in the data measured during cool-down being shifted slightly to colder
temperatures compared to those extracted from the slower scans in which we
measure a full spectrum at a more stable temperature.

Figure S3: The spectra of (a) a bilayer junction and (b) an 11-layer junction
as a function of temperature, with blue (red) indicating lower (higher)
temperature. The full range of temperatures is indicated in the panels.

####  I.4 Parallel Magnetic Field

Fig. S4 displays the spectra of a bilayer and a 20-layer sample (bulk) under
parallel magnetic field. The effect of Ising SOC is evident in the remarkably
small effect of fields up to 8.5T on the bilayer sample. The bulk sample, on
the other hand, is much more noticeably affected.

Figure S4: Parallel field response of (a) a bilayer and (b) a 20-layer (bulk)
device. Ising SOC is evident in the robustness of the bilayer gap. (c)
Overlaid individual spectra from (a,b), with zero field in blue and high
fields in red. The full range of fields is indicated in the panels.

####  I.5 Metallic H-TMD Sequence with Additional Materials

Fig. S5 displays the $H_{C2}$ vs. $T_{C}$ data presented in Fig. 4 of the main
text along with additional datasets:

  * •

Liquid-ion gated MoS2 Ye  _et al._ (2012), displaying $H_{C2}$ $\propto\Delta$
behavior indicative of dirty-limit superconductivity;

  * •

Liquid-ion gated MoS2 Costanzo  _et al._ (2016), displaying $H_{C2}$
$\propto\Delta^{2}$ behavior but with a different constant of proportionality
than the metallic H-TMD Sequence;

  * •

Data measured on 4Hb-TaS2 Ribak  _et al._ (2020), which is a qualitatively
different superconductor with reports of various exotic phenomena;

  * •

Intercalated TaS2 compounds: TaS2 intercalated with pyridine molecules
(highest $T_{C}$), and misfit compounds of 1H-TaS2 sheets separated by square
layers of PbS (middle) and SnS (lowest).

Figure S5: $H_{C2}$ vs. $T_{C}^{2}$ relation of various TMDs. Notably, the
metallic H-TMD sequence retains $H_{C2}$ $\propto T_{C}^{2}$ with the same
slope, liquid-ion gated MoS2 can display either the dirty-limit result
$H_{C2}$ $\propto T_{C}$ or the clean-limit result in different experiments
(but with a different slope than the metallic ones). Less conventional TaS2
compounds such as intercalated samples or the 4Hb polytype can behave quite
differently. Data adapted from Xi  _et al._ (2016); Yokota  _et al._ (2000);
Tissen  _et al._ (2013); Costanzo  _et al._ (2016); Ye  _et al._ (2012); Ribak
_et al._ (2020).

###  II Theoretical Computation of the upper critical field

####  II.1 Tight-binding model

#####  II.1.1 Without CDW

To model TMD materials in the absence of CDW order, we use the 3-orbital
tight-binding model developed in Ref. Liu  _et al._ (2013). In the basis
${\ket{\psi_{1}}}={\ket{d_{z^{2}}},\ {\ket{\psi_{2}}}=\ket{d_{xy}},\
{\ket{\psi_{3}}}=\ket{d_{x^{2}-y^{2}}}}$, the normal state tight-binding
Hamiltonian is given by

$\displaystyle H_{{\boldsymbol{k}}}=\sigma_{0}\otimes
H_{0{\boldsymbol{k}}}+\frac{\lambda}{2}\sigma_{z}\otimes L_{z},$ (1)

where $L_{z}$ in the spin-orbit-coupling term is the z-component of angular
momentum operator

$\displaystyle L_{z}=\begin{pmatrix}0&0&0\\\ 0&0&-i\\\ 0&i&0\end{pmatrix},$
(2)

and

$\displaystyle H_{0{\boldsymbol{k}}}=\begin{pmatrix}V_{0}&V_{1}&V_{2}\\\
V_{1}^{*}&V_{11}&V_{12}\\\ V_{2}^{*}&V_{12}^{*}&V_{22}\end{pmatrix},$ (3)

where

$\displaystyle V_{0}=$
$\displaystyle\epsilon_{1}+2t_{0}(2\cos\alpha\cos\beta+\cos
2\alpha)+2r_{0}(2\cos 3\alpha\cos\beta+\cos 2\beta)+2u_{0}(2\cos 2\alpha\cos
2\beta+\cos 4\alpha)$ (4) $\displaystyle V_{1}=$
$\displaystyle-2\sqrt{3}t_{2}\sin\alpha\sin\beta+2\left(r_{1}+r_{2}\right)\sin
3\alpha\sin\beta-2\sqrt{3}u_{2}\sin 2\alpha\sin 2\beta$ (5)
$\displaystyle+i\left(2t_{1}\sin\alpha(2\cos\alpha+\cos\beta)+2\left(r_{1}-r_{2}\right)\sin
3\alpha\cos\beta+2u_{1}\sin 2\alpha(2\cos 2\alpha+\cos 2\beta)\right),$ (6)
$\displaystyle V_{2}=$ $\displaystyle 2t_{2}(\cos
2\alpha-\cos\alpha\cos\beta)-\frac{2}{\sqrt{3}}\left(r_{1}+r_{2}\right)(\cos
3\alpha\cos\beta-\cos 2\beta)+2u_{2}(\cos 4\alpha-\cos 2\alpha\cos 2\beta)$
(7)
$\displaystyle+i\left(\sqrt{3}t_{1}\cos\alpha\sin\beta+\frac{2}{\sqrt{3}}\sin\beta\left(r_{1}-r_{2}\right)(\cos
3\alpha+2\cos\beta)+2\sqrt{3}u_{1}\cos 2\alpha\sin 2\beta\right),$
$\displaystyle V_{11}=$
$\displaystyle\epsilon_{2}+\left(t_{11}+3t_{22}\right)\cos\alpha\cos\beta+2t_{11}\cos
2\alpha+4r_{11}\cos 3\alpha\cos\beta+2\left(r_{11}+\sqrt{3}r_{12}\right)\cos
2\beta$ (8) $\displaystyle+\left(u_{11}+3u_{22}\right)\cos 2\alpha\cos
2\beta+2u_{11}\cos 4\alpha,$ $\displaystyle V_{12}=$
$\displaystyle\sqrt{3}\left(t_{22}-t_{11}\right)\sin\alpha\sin\beta+4r_{12}\sin
3\alpha\sin\beta+\sqrt{3}\left(u_{22}-u_{11}\right)\sin 2\alpha\sin 2\beta,$
(9) $\displaystyle+i\left(4t_{12}\sin\alpha(\cos\alpha-\cos\beta)+4u_{12}\sin
2\alpha(\cos 2\alpha-\cos 2\beta)\right),$ (10) $\displaystyle V_{22}=$
$\displaystyle\epsilon_{2}+\left(3t_{11}+t_{22}\right)\cos\alpha\cos\beta+2t_{22}\cos
2\alpha+2r_{11}(2\cos 3\alpha\cos\beta+\cos
2\beta)+\frac{2}{\sqrt{3}}r_{12}(4\cos 3\alpha\cos\beta-\cos 2\beta)$
$\displaystyle+\left(3u_{11}+u_{22}\right)\cos 2\alpha\cos 2\beta+2u_{22}\cos
4\alpha,$ (11)

with $\alpha=k_{x}a/2$ and $\beta=\sqrt{3}k_{y}a/2$, where $a$ is the lattice
constant.

#####  II.1.2 With the CDW

According to Ref. Flicker and van Wezel (2016), in the commensurate 3Q CDW
state the metal atoms’ distortion is

$\displaystyle
u({\boldsymbol{R}}_{a})=\sum_{i=1}^{3}u_{0}\hat{q}_{i}\cos({\boldsymbol{q}}_{i}\cdot{\boldsymbol{R}}_{a}+\phi)\,,$
(12)

where ${\boldsymbol{q}}_{1}=\frac{{\boldsymbol{b}}_{1}}{3}$,
${\boldsymbol{q}}_{2}=\frac{-{\boldsymbol{b}}_{1}+{\boldsymbol{b}}_{2}}{3}$,
${\boldsymbol{q}}_{3}=-\frac{{\boldsymbol{b}}_{2}}{3}$ are three CDW-vectors.
Here ${\boldsymbol{b_{1}}},{\boldsymbol{b_{2}}}$ are primitive vectors in
reciprocal space forming 1BZ of the system in the absence of CDW. In the CDW
phase, the little BZ (lBZ) is formed by vectors ${\boldsymbol{b_{1}}}/3$ and
${\boldsymbol{b_{2}}}/3$.

In principle, one can derive the CDW-term completely from symmetry
considerations. But to avoid mistakes and aid ourselves with something
tangible, we assume that CDW order develops due to electron-phonon interaction
(which is likely trueFlicker and van Wezel (2016)):

$\displaystyle
V_{e-p}=\sum_{a}V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}}-{\boldsymbol{u_{a}}})-\sum_{a}V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})=-\sum_{a}\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot{\boldsymbol{u_{a}}}\,.$ (13)

Consequently, the parameters of the tight-binding model are renoromalized.
This can be accounted for by computing the matrix element for scattering
between ${\boldsymbol{k^{\prime}}}$ and ${\boldsymbol{k}}$ due to the CDW
potential

$\displaystyle-\matrixelement{\psi_{i}({\boldsymbol{k}})}{V_{e-p}}{\psi_{j}({\boldsymbol{k^{\prime}}})}$
$\displaystyle=\sum_{a}\int
d^{d}{\boldsymbol{r}}\psi_{i{\boldsymbol{k}}}({\boldsymbol{r}})^{*}\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot{\boldsymbol{u_{a}}}\psi_{j{\boldsymbol{k^{\prime}}}}({\boldsymbol{r}})$
(14)
$\displaystyle=\sum_{a}\frac{1}{N^{3}}\sum_{{\boldsymbol{R}},{\boldsymbol{R^{\prime}}}}e^{i{\boldsymbol{k^{\prime}}}{\boldsymbol{R^{\prime}}}-i{\boldsymbol{k}}{\boldsymbol{R}}}\int
d^{d}{\boldsymbol{r}}\phi_{i}({\boldsymbol{r}})^{*}\gradient
V_{0}({\boldsymbol{r}}-({\boldsymbol{R_{a}}}-{\boldsymbol{R}}))\cdot{\boldsymbol{u_{a}}}\phi_{j}({\boldsymbol{r}}-({\boldsymbol{R^{\prime}}}-{\boldsymbol{R}}))$
$\displaystyle=\sum_{a}\frac{1}{N^{3}}\sum_{{\boldsymbol{R}},{\boldsymbol{R^{\prime}}}}e^{i{\boldsymbol{k^{\prime}}}{\boldsymbol{R^{\prime}}}-i{\boldsymbol{k}}{\boldsymbol{R}}}\int
d^{d}{\boldsymbol{r}}\phi_{i}({\boldsymbol{r}})^{*}\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot{\boldsymbol{u_{a}}}({\boldsymbol{R_{a}}}+{\boldsymbol{R}})\phi_{j}({\boldsymbol{r}}-({\boldsymbol{R^{\prime}}}-{\boldsymbol{R}})).$

Concentrating on CDW-vector
${\boldsymbol{q}}\in\\{{\boldsymbol{q_{1}}},{\boldsymbol{q_{2}}},{\boldsymbol{q_{3}}}\\}$
such that ${\boldsymbol{k}}={\boldsymbol{k^{\prime}}}+{\boldsymbol{q}}$, and
denoting $\gamma_{a}({\boldsymbol{r}},{\boldsymbol{R}})=\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot{\boldsymbol{u_{a}}}({\boldsymbol{R_{a}}}+{\boldsymbol{R}})$,
$\gamma({\boldsymbol{r}},{\boldsymbol{R}})=\sum_{a}\gamma_{a}({\boldsymbol{r}},{\boldsymbol{R}})$,
we obtain

$\displaystyle-\matrixelement{\psi_{i}({\boldsymbol{k^{\prime}}}+{\boldsymbol{q}})}{V_{e-p}}{\psi_{j}({\boldsymbol{k^{\prime}}})}=\frac{1}{N^{3}}\sum_{{\boldsymbol{R}}}e^{-i{\boldsymbol{q}}{\boldsymbol{R}}}\left(\sum_{{\boldsymbol{R^{\prime}}}}e^{i{\boldsymbol{k^{\prime}}}{\boldsymbol{R^{\prime}}}}\int
d^{d}{\boldsymbol{r}}\phi_{i}({\boldsymbol{r}})^{*}\gamma({\boldsymbol{r}},{\boldsymbol{R}})\phi_{j}({\boldsymbol{r}}-{\boldsymbol{R^{\prime}}}))\right).$
(15)

Finally, denoting
$F_{n}({\boldsymbol{r}})=\sum_{{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{n}}}{\boldsymbol{R^{\prime}}}}\gamma({\boldsymbol{r}},{\boldsymbol{R^{\prime}}})$
and $A^{n}_{ij}({\boldsymbol{R}})=\frac{1}{N^{3}}\int
d^{d}{\boldsymbol{r}}\phi_{i}({\boldsymbol{r}})^{*}F_{n}({\boldsymbol{r}})\phi_{j}({\boldsymbol{r}}-{\boldsymbol{R}}))$,
we obtain

$\displaystyle\Lambda_{{\boldsymbol{q}}_{n}}=\matrixelement{\psi_{i}({\boldsymbol{k}}+{\boldsymbol{q_{i}}})}{V_{e-p}}{\psi_{j}({\boldsymbol{k}})}=-\sum_{{\boldsymbol{R}}}e^{i{\boldsymbol{k}}{\boldsymbol{R}}}A^{n}_{ij}({\boldsymbol{R}}).$
(16)

Now, we will try to establish the relation between
$\Lambda_{{\boldsymbol{q_{i}}}}$, $i=1,2,3$.

First, denoting

$\displaystyle f_{1}^{i}({\boldsymbol{r}})=\sum_{a}\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot\hat{q}_{i}\cos({\boldsymbol{q}}_{i}\cdot{\boldsymbol{R}}_{a}+\phi),$
(17) $\displaystyle f_{2}^{i}({\boldsymbol{r}})=\sum_{a}\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot\hat{q}_{i}\sin({\boldsymbol{q}}_{i}\cdot{\boldsymbol{R}}_{a}+\phi),$

we can write

$\displaystyle
F_{i}({\boldsymbol{r}})=u_{0}\sum_{{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{i}}}{\boldsymbol{R^{\prime}}}}\sum_{j}\left(f_{1}^{j}({\boldsymbol{r}})\cos({\boldsymbol{q_{j}}}\cdot{\boldsymbol{R^{\prime}}})-f_{2}^{j}({\boldsymbol{r}})\sin({\boldsymbol{q_{j}}}\cdot{\boldsymbol{R^{\prime}}})\right).$
(18)

Next, denoting

$\displaystyle
g_{ij}^{1}=\sum_{{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{i}}}{\boldsymbol{R^{\prime}}}}\cos({\boldsymbol{q_{j}}}\cdot{\boldsymbol{R^{\prime}}}),$
(19) $\displaystyle
g_{ij}^{2}=\sum_{{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{i}}}{\boldsymbol{R^{\prime}}}}\sin({\boldsymbol{q_{j}}}\cdot{\boldsymbol{R^{\prime}}}),$

we find, noting that
${\boldsymbol{q_{1}}}+{\boldsymbol{q_{2}}}+{\boldsymbol{q_{3}}}=0$,

$\displaystyle g_{ij}^{1}+ig_{ij}^{2}$
$\displaystyle=\sum_{{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{i}}}{\boldsymbol{R^{\prime}}}}e^{i{\boldsymbol{q_{j}}}\cdot{\boldsymbol{R^{\prime}}}}=N^{3}\delta_{ij},$
(20) $\displaystyle g_{ij}^{1}-ig_{ij}^{2}$
$\displaystyle=\sum_{{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{i}}}{\boldsymbol{R^{\prime}}}}e^{-i{\boldsymbol{q_{j}}}\cdot{\boldsymbol{R^{\prime}}}}=0.$

From which we obtain

$\displaystyle g_{ij}^{1}=\frac{N^{3}}{2}\delta_{ij},\ \
g_{ij}^{2}=-i\frac{N^{3}}{2}\delta_{ij}.$ (21)

Thus, the expression for $F_{i}({\boldsymbol{r}})$ simplifies to

$\displaystyle
F_{i}({\boldsymbol{r}})=u_{0}\frac{N^{3}}{2}\left(f_{1}^{i}+if_{2}^{i}\right)=u_{0}\frac{N^{3}}{2}\sum_{a}\gradient
V_{0}({\boldsymbol{r}}-{\boldsymbol{R_{a}}})\cdot\hat{q}_{i}e^{i{\boldsymbol{q}}_{i}\cdot{\boldsymbol{R}}_{a}+i\phi}.$
(22)

Now, using the symmetry of the system, we can establish the relation between
$A_{ij}^{n}({\boldsymbol{R}})$, $n=1,2,3$. For monolayer, the point group
symmetry is $D_{3h}$, with symmetry elements
$\mathcal{G}\in\\{E,C_{3},\bar{C}_{3},C_{2},C_{2}^{\prime},C_{2}^{\prime\prime},S_{3},\bar{S}_{3},\sigma_{v},\sigma_{v^{\prime}},\sigma_{v^{\prime\prime}}\\}$.

$\displaystyle
A_{ij}^{n}({\boldsymbol{R}})=\matrixelement{\phi_{i}({\boldsymbol{r}})}{F_{n}({\boldsymbol{r}})}{\phi_{j}({\boldsymbol{r}}-{\boldsymbol{R}})}$
(23)

Under transformation operation $\mathcal{G}$:

$\displaystyle
A_{ij}^{n}(\mathcal{G}{\boldsymbol{R}})=\matrixelement{\phi_{i}(\mathcal{G}{\boldsymbol{r}})}{F_{n}(\mathcal{G}{\boldsymbol{r}})}{\phi_{j}(g({\boldsymbol{r}}-{\boldsymbol{R}}))}=\matrixelement{\phi_{k}({\boldsymbol{r}})\Gamma_{\mathcal{G}^{-1},ki}^{r_{i}}}{F_{n}(\mathcal{G}{\boldsymbol{r}})}{\phi_{l}({\boldsymbol{r}}-{\boldsymbol{R}})\Gamma_{\mathcal{G}^{-1},lj}^{r_{j}}},$
(24)

where $r_{i},r_{j}$ stand for the representation of $\phi_{i},\phi_{j}$
functions. The transformation properties of $F_{n}({\boldsymbol{r}})$ are
obtained using Eq. 22 and are specfied in Table 1.

| $C_{3}$ | $\bar{C}3$ | $C_{2}$ | $C_{2}^{\prime}$ | $C_{2}^{\prime\prime}$ | $\sigma_{h}$ | $S_{3}$ | $\bar{S}_{3}$ | $\sigma_{v}$ | $\sigma_{v}^{\prime}$ | $\sigma_{v}^{\prime\prime}$
---|---|---|---|---|---|---|---|---|---|---|---
$F_{1}({\boldsymbol{r}})$ | $F_{3}(C_{3}{\boldsymbol{r}})$ | $F_{2}(\bar{C_{3}}{\boldsymbol{r}})$ | $F_{1}(C_{2}{\boldsymbol{r}})$ | $F_{2}(C_{2}^{\prime}{\boldsymbol{r}})$ | $F_{3}(C_{2}^{\prime\prime}{\boldsymbol{r}})$ | $F_{1}({\boldsymbol{r}})$ | $F_{3}(S_{3}{\boldsymbol{r}})$ | $F_{2}(\bar{S_{3}}{\boldsymbol{r}})$ | $F_{1}(\sigma_{v}{\boldsymbol{r}})$ | $F_{2}(\sigma_{v}^{\prime}{\boldsymbol{r}})$ | $F_{3}(\sigma_{v}^{\prime\prime}{\boldsymbol{r}})$
$F_{2}({\boldsymbol{r}})$ | $F_{1}(C_{3}{\boldsymbol{r}})$ | $F_{3}(\bar{C_{3}}{\boldsymbol{r}})$ | $F_{3}(C_{2}{\boldsymbol{r}})$ | $F_{1}(C_{2}^{\prime}{\boldsymbol{r}})$ | $F_{2}(C_{2}^{\prime\prime}{\boldsymbol{r}})$ | – | – | – | $F_{3}(\sigma_{v}{\boldsymbol{r}})$ | $F_{1}(\sigma_{v}^{\prime}{\boldsymbol{r}})$ | $F_{2}(\sigma_{v}^{\prime\prime}{\boldsymbol{r}})$
$F_{3}({\boldsymbol{r}})$ | $F_{2}(C_{3}{\boldsymbol{r}})$ | $F_{1}(\bar{C_{3}}{\boldsymbol{r}})$ | $F_{2}(C_{2}{\boldsymbol{r}})$ | $F_{3}(C_{2}^{\prime}{\boldsymbol{r}})$ | $F_{1}(C_{2}^{\prime\prime}{\boldsymbol{r}})$ | – | – | – | $F_{2}(\sigma_{v}{\boldsymbol{r}})$ | $F_{3}(\sigma_{v}^{\prime}{\boldsymbol{r}})$ | $F_{1}(\sigma_{v}^{\prime\prime}{\boldsymbol{r}})$
Table 1: Transformations of $F_{i}({\boldsymbol{r}})$, $i=1,2,3$.

Finally, in matrix form, we have

$\displaystyle
A^{n}(g{\boldsymbol{R}})=\Gamma_{g}^{r_{i}}A^{\tilde{n}}(R)\Gamma_{g}^{r_{j}\dagger},$
(25)

where $\tilde{n}$ is found from Table 1. This allows us to symmetry-constrain
the form of the CDW term:

For ${\boldsymbol{R}}=0$, denoting $B^{n}=A^{n}({\boldsymbol{0}})$, $n=1,2,3$,
we then find system of equations

$\displaystyle B^{1}$
$\displaystyle=\Gamma_{\sigma_{\nu}}^{r_{i}}B^{1}\Gamma_{\sigma_{\nu}}^{r_{j}\dagger},\
\ \ \
B^{2}=\Gamma_{\sigma_{\nu}}^{r_{i}}B^{3}\Gamma_{\sigma_{\nu}}^{r_{j}\dagger},$
(26) $\displaystyle B^{1}$
$\displaystyle=\Gamma_{\bar{C_{3}}}^{r_{i}}B^{1}\Gamma_{\bar{C_{3}}}^{r_{j}\dagger},\
\ \ \ B^{2}=\Gamma_{C_{3}}^{r_{i}}B^{3}\Gamma_{C_{3}}^{r_{j}\dagger},$ (27)
$\displaystyle B^{1}$
$\displaystyle=\Gamma_{\sigma_{\nu}^{\prime}}^{r_{i}}B^{2}\Gamma_{\sigma_{\nu}^{\prime}}^{r_{j}\dagger},\
\ \ \
B^{2}=\Gamma_{\sigma_{\nu}^{\prime\prime}}^{r_{i}}B^{2}\Gamma_{\sigma_{\nu}^{\prime\prime}}^{r_{j}\dagger},$
(28) $\displaystyle B^{1}$
$\displaystyle=\Gamma_{C_{3}}^{r_{i}}B^{2}\Gamma_{C_{3}}^{r_{j}\dagger},\ \ \
\
B^{3}=\Gamma_{\sigma_{\nu}^{\prime}}^{r_{i}}B^{3}\Gamma_{\sigma_{\nu}^{\prime}}^{r_{j}\dagger},$
(29) $\displaystyle B^{1}$
$\displaystyle=\Gamma_{\sigma_{\nu}^{\prime\prime}}^{r_{i}}B^{3}\Gamma_{\sigma_{\nu}^{\prime\prime}}^{r_{j}\dagger},$
(30)

Solving this system, we find that

$\displaystyle B^{1}=B^{2}=B^{3}=\begin{pmatrix}\tilde{b}_{1}&0&0\\\
0&\tilde{b}_{2}&0\\\ 0&0&\tilde{b}_{2}\end{pmatrix},$ (32)

where $\tilde{b}_{1}$ and $\tilde{b}_{2}$ are some complex numbers.
Analogously, we can derive similar relations for vectors ${\boldsymbol{R}}$
connecting nearest neighbours, next nearest neigbors, and so on. The free
parameters like $\tilde{b}_{1},\ \tilde{b}_{2}$ in Eq. 32 then could be
determined self-consistently in the BdG-like approach applied to the CDW
Hamiltonian.

Overall, for the commensurate CDW, the Hamiltonian is

$\displaystyle\hat{H}=\sum_{\boldsymbol{k}}c_{\boldsymbol{k}}^{\dagger}H_{\boldsymbol{k}}c_{\boldsymbol{k}}+\sum_{n,{\boldsymbol{k}}}\left(c_{{\boldsymbol{k}}+{\boldsymbol{q_{n}}}}^{\dagger}\Lambda^{n}_{\boldsymbol{k}}c_{\boldsymbol{k}}+\mathrm{cc}.\right),$
(33)

where
$c_{\boldsymbol{k}}=\left(c_{{\boldsymbol{k}}z^{2}\uparrow},c_{{\boldsymbol{k}}xy\uparrow},c_{{\boldsymbol{k}}x^{2}-y^{2}\uparrow},c_{{\boldsymbol{k}}z^{2}\downarrow},c_{{\boldsymbol{k}}xy\downarrow},c_{{\boldsymbol{k}}x^{2}-y^{2}\downarrow}\right)$,
and $c^{\dagger}_{{\boldsymbol{k}}\alpha s}$ creates an electron with momentum
${\boldsymbol{k}}$ in orbital $\alpha$ and spin $s$. In the space of vectors
$\\{c_{\boldsymbol{k}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}},c_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern 1.5mu}}}}\\}$, the CDW
Hamiltonian takes form

$\displaystyle
H_{\mathrm{CDW}}({\boldsymbol{k}})=\begin{pmatrix}H_{\boldsymbol{k}}&\Lambda^{1}_{\boldsymbol{k}}&\Lambda^{2}_{\boldsymbol{k}}&\Lambda^{3}_{\boldsymbol{k}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}&0&0\\\
\Lambda^{1\dagger}_{\boldsymbol{k}}&H_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&0&0&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern 1.5mu}}}}\\\
\Lambda^{2\dagger}_{\boldsymbol{k}}&0&H_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&0&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern 1.5mu}}}}\\\
\Lambda^{3\dagger}_{\boldsymbol{k}}&0&0&H_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern 1.5mu}}}}\\\
\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&H_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}&0&0&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}\\\
\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&0&H_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}&0&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}\\\
\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{1}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{2}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{3}}}}&0&0&H_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}\\\
0&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}&H_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern 1.5mu}}}}&0\\\
0&\Lambda^{3}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{1}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{2}_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern
1.5mu}}}}&\Lambda^{1\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{12}}}}&\Lambda^{3\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{13}}}}&\Lambda^{2\dagger}_{{\boldsymbol{k}}+{\boldsymbol{q_{23}}}}&0&H_{{\boldsymbol{k}}+{\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern 1.5mu}}}}\end{pmatrix},$ (34)

where ${\boldsymbol{k}}$ is restricted to lBZ,
${\boldsymbol{k}}\in\mathrm{lBZ}$, and we used notation
${\boldsymbol{q_{12}}}={\boldsymbol{q_{1}}}+{\boldsymbol{q_{2}}},\
{\boldsymbol{q_{13}}}={\boldsymbol{q_{1}}}+{\boldsymbol{q_{3}}},\
{\boldsymbol{q_{23}}}={\boldsymbol{q_{2}}}+{\boldsymbol{q_{3}}},\
{\boldsymbol{q_{\mkern 1.5mu\overline{\mkern-1.5mu12\mkern-1.5mu}\mkern
1.5mu}}}={\boldsymbol{q_{1}}}-{\boldsymbol{q_{2}}},\ {\boldsymbol{q_{\mkern
1.5mu\overline{\mkern-1.5mu13\mkern-1.5mu}\mkern
1.5mu}}}={\boldsymbol{q_{1}}}-{\boldsymbol{q_{3}}}$.

The values of the CDW gaps $\Lambda^{i}_{\boldsymbol{k}}$ can be determined
self-consistently. However, expecting that results of interest will not depend
on the fine details of the CDW-term (an argument supporting this will be
present in the next section), we use the functional form determined in Flicker
and van Wezel (2016) for NbSe2 in the band basis and adapt it to our approach.
We set

$\displaystyle\Lambda^{1}_{\boldsymbol{k}}$
$\displaystyle=\Lambda^{2}_{\boldsymbol{k}}=\Lambda^{3}_{\boldsymbol{k}}=b\begin{pmatrix}\Lambda({\boldsymbol{k}})&0&0\\\
0&\Lambda({\boldsymbol{k}})&0\\\ 0&0&\Lambda({\boldsymbol{k}})\end{pmatrix},$
(35)

where

$\displaystyle\Lambda({\boldsymbol{k}})$ $\displaystyle=c_{1}(2\cos(\alpha
k_{x})\cos(\beta k_{y})+\cos(\alpha k_{x}))+c_{2}(2\cos(3\alpha
k_{x})\cos(\beta k_{y})+\cos(2\beta k_{y}))+c_{3}(2\cos(2\alpha
k_{x})\cos(2\beta k_{y})+\cos(4\alpha k_{x}))$ (36)
$\displaystyle+c_{4}(\cos(5\alpha k_{x})\cos(\beta k_{y})+\cos(4\alpha
k_{x})\cos(2\beta k_{y})+\cos(\alpha k_{x})\cos(3\beta
k_{y}))+c_{5}(2\cos(3\alpha k_{x})\cos(3\beta k_{y})+\cos(6\alpha
k_{x}))+c_{0},$ $\displaystyle\alpha$ $\displaystyle=\frac{a}{2},\
\beta=\frac{\sqrt{3}a}{2}\,.$ (37)

The parameter $b$ is considered as a phenomenological parameter hereafter. It
is used to control the overall strength of the CDW-term. Following Flicker and
van Wezel (2016), where the constants $c_{i}$ are determined self-
consistently, we determine $c_{i}$ from the condition that
$\Lambda({\boldsymbol{k}})$ is not zero only at one high-symmetry point,
namely at $K$ we require $\Lambda({\boldsymbol{k}})=1.25$. We find

$\displaystyle c_{0}$ $\displaystyle=0.135;\ c_{1}=-0.114;$ (38)
$\displaystyle c_{2}$ $\displaystyle=0.116;\ c_{3}=-0.060;$ $\displaystyle
c_{4}$ $\displaystyle=-0.095;\ c_{5}=0.141;$

We will use the same form of the CDW term for TaS2 as well. We plot
$\Lambda({\boldsymbol{k}})$ function with superimposed Fermi surfaces for
NbSe2 and TaS2 in Fig. S6.

(a)

(b)

Figure S6: Density plot of $\Lambda({\boldsymbol{k}})$ in 1BZ with
superimposed FS for NbSe2 (a) and TaS2 (b).

####  II.2 Superconductivity – BCS theory

Now we add interaction that leads to $s$-wave superconductivity. Using group
theory arguments as outlined in Ref. Möckli and Khodas (2018), one can
conclude that allowed $s$-wave correlations (in real space) are

$\displaystyle\expectationvalue{\Psi_{1,i}}=\expectationvalue{c_{iz^{2},\uparrow}c_{iz^{2},\downarrow}},\
\mathrm{and}\
\expectationvalue{\Psi_{2,i}}=\expectationvalue{c_{ixy,\uparrow}c_{ixy,\downarrow}}+\expectationvalue{c_{ix^{2}-y^{2},\uparrow}c_{ix^{2}-y^{2},\downarrow}},$
(39)

where $i$ is the index indicating the location in the lattice. We then can
write the pairing contact interaction in the form

$\displaystyle\hat{H}_{\mathrm{I}}=-\sum_{i}\left[U_{1}\Psi_{1,i}^{\dagger}\Psi_{1,i}+U_{2}\Psi_{2,i}^{\dagger}\Psi_{2,i}+U_{3}\left(\Psi_{1,i}^{\dagger}\Psi_{2,i}+\Psi_{2,i}^{\dagger}\Psi_{1,i}\right)\right]$
(40)

Correspondingly, in the momentum space, restricting interaction to most
relevant channel (i.e., the one corresponding to zero momentum of the center
of mass of a pair), we write down effective interaction in the form

$\displaystyle\hat{H}_{\mathrm{I}}=-\sum_{k,k^{\prime}}\left[U_{1}\Psi_{1,{\boldsymbol{k}}}^{\dagger}\Psi_{1,{\boldsymbol{k^{\prime}}}}+U_{2}\Psi_{2,{\boldsymbol{k}}}^{\dagger}\Psi_{2,{\boldsymbol{k^{\prime}}}}+U_{3}\left(\Psi_{1,{\boldsymbol{k}}}^{\dagger}\Psi_{2,{\boldsymbol{k^{\prime}}}}+\Psi_{2,{\boldsymbol{k}}}^{\dagger}\Psi_{1,{\boldsymbol{k^{\prime}}}}\right)\right],$
(41)

where
$\Psi_{1,{\boldsymbol{k}}}=c_{z^{2}{\boldsymbol{k}},\uparrow}c_{z^{2}-{\boldsymbol{k}},\downarrow}$
and
$\Psi_{2,{\boldsymbol{k}}}=c_{xy{\boldsymbol{k}},\uparrow}c_{xy-{\boldsymbol{k}},\downarrow}+c_{x^{2}-y^{2}{\boldsymbol{k}},\uparrow}c_{x^{2}-y^{2}-{\boldsymbol{k}},\downarrow}$.

Denoting
$\Delta_{i}=U_{i}\expectationvalue{c_{i{\boldsymbol{p}},\uparrow}c_{i-{\boldsymbol{p}},\downarrow}}$,
where $i=1,2,3$ correspond to $z^{2},xy,x^{2}-y^{2}$, respectively, in the
mean-field treatment the interaction takes form

$\displaystyle\hat{H}_{\mathrm{I,MF}}=\sum_{\boldsymbol{k}}\left(\Delta_{1}+\frac{U_{3}}{U_{2}}(\Delta_{2}+\Delta_{3})\right)c_{1{\boldsymbol{k}}\uparrow}^{\dagger}c_{1-{\boldsymbol{k}}\downarrow}^{\dagger}+\left(\Delta_{2}+\Delta_{3}+\frac{U_{3}}{U_{1}}\Delta_{1}\right)\left(c_{2{\boldsymbol{k}}\uparrow}^{\dagger}c_{2-{\boldsymbol{k}}\downarrow}^{\dagger}+c_{3{\boldsymbol{k}}\uparrow}^{\dagger}c_{3-{\boldsymbol{k}}\downarrow}^{\dagger}\right).$
(42)

To proceed, we use the BdG formalism

$\displaystyle
H_{\mathrm{BdG}}({\boldsymbol{k}})=\frac{1}{2}\begin{pmatrix}H({\boldsymbol{k}})&\hat{\Delta}\\\
\hat{\Delta}^{\dagger}&-H(-{\boldsymbol{k}})^{T}\end{pmatrix},\ \
\mathrm{and}\ \
\hat{H}=\sum_{{\boldsymbol{k}}}\Psi_{\boldsymbol{k}}^{\dagger}H_{BdG}({\boldsymbol{k}})\Psi_{\boldsymbol{k}},$
(43)

where $\hat{\Delta}$ is the pairing matrix and $\Psi_{{\boldsymbol{k}}}$ is
the fermion destruction operator in the Nambu space. In the case of no CDW,

$\displaystyle\Psi_{\boldsymbol{k}}$
$\displaystyle=\left(c_{1{\boldsymbol{k}}\uparrow},c_{2{\boldsymbol{k}}\uparrow},c_{3{\boldsymbol{k}}\uparrow},c_{1{\boldsymbol{k}}\downarrow},c_{2{\boldsymbol{k}}\downarrow},c_{3{\boldsymbol{k}}\downarrow},c_{1-{\boldsymbol{k}}\uparrow}^{\dagger},c_{2-{\boldsymbol{k}}\uparrow}^{\dagger},c_{3-{\boldsymbol{k}}\uparrow}^{\dagger},c_{1-{\boldsymbol{k}}\downarrow}^{\dagger},c_{2-{\boldsymbol{k}}\downarrow}^{\dagger},c_{3-{\boldsymbol{k}}\downarrow}^{\dagger}\right)^{T},$
(44) $\displaystyle\hat{\Delta}$ $\displaystyle=i\sigma_{y}\otimes\Delta_{o},\
\ \mathrm{where}\ \
\Delta_{o}=\begin{pmatrix}\Delta_{1}+\frac{U_{3}}{U_{2}}(\Delta_{2}+\Delta_{3})&0&0\\\
0&\Delta_{2}+\Delta_{3}+\frac{U_{3}}{U_{1}}\Delta_{1}&0\\\
0&0&\Delta_{2}+\Delta_{3}+\frac{U_{3}}{U_{1}}\Delta_{1}\end{pmatrix},$

where $\sigma_{y}$ is the Pauli matrix.

When the CDW is present, the Nambu space is formed by

$\displaystyle\Psi_{\boldsymbol{k}}=\left(c_{{\boldsymbol{k_{1}}}\uparrow},...,c_{{\boldsymbol{k_{9}}}\uparrow},c_{{\boldsymbol{k_{1}}}\downarrow},...,c_{{\boldsymbol{k_{9}}}\downarrow},c_{-{\boldsymbol{k_{1}}}\uparrow}^{\dagger
T},...,c_{-{\boldsymbol{k_{9}}}\uparrow}^{\dagger
T},c_{-{\boldsymbol{k_{1}}}\downarrow}^{\dagger
T},...,c_{-{\boldsymbol{k_{9}}}\downarrow}^{\dagger T}\right)^{T},$ (45)

where
$c_{\boldsymbol{k_{1}s}}=\left(c_{\boldsymbol{1k_{1}s}},c_{\boldsymbol{2k_{1}s}},c_{\boldsymbol{3k_{1}s}}\right)$
with $s$ being a spin index, ${\boldsymbol{k_{1}}}\in\mathrm{lBZ}$ and
${\boldsymbol{k_{2}}},...,{\boldsymbol{k_{9}}}$ vectors belong to domains
$\mathrm{lBZ}_{2},...,\mathrm{lBZ}_{9}$, respectively, in
${\boldsymbol{k}}$-space formed by shifting $\mathrm{lBZ}$ by linear
combinations of ${\boldsymbol{q_{1}}},{\boldsymbol{q_{2}}},\ \mathrm{and}\
{\boldsymbol{q_{3}}}$, such that the union of these domains with
$\mathrm{lBZ}$ covers the domain equivalent to $\mathrm{1BZ}$. We note that in
this case $H(-{\boldsymbol{k}})$ in Eq. 34 is not simply formed by taking
${\boldsymbol{k}}\rightarrow-{\boldsymbol{k}}$ in $H({\boldsymbol{k}})$ since,
for example,
$-{\boldsymbol{k}}+{\boldsymbol{q_{1}}}\neq-({\boldsymbol{k}}+{\boldsymbol{q_{1}}})$:
proper rearrangement of blocks is necessary.

In principle, one can work with such defined Hamiltonian numerically. But the
size of the matrices is quite large, and practice showed that straightforward
computations with this Hamiltonian defined in $k$-space takes too much time.
However, the special case $U_{1}=U_{2}=U_{3}$ allows significant analytical
advancement.

#####  II.2.1 Analytical solution in the case $U_{1}=U_{2}=U_{3}$

Let us first rewrite Eq. 43 more explicitly expanding in spin space:

$\displaystyle
H_{\mathrm{BdG}}({\boldsymbol{k}})=\frac{1}{2}\begin{pmatrix}H_{\uparrow}({\boldsymbol{k}})&0&0&\hat{\Delta}_{c}\\\
0&H_{\downarrow}({\boldsymbol{k}})&-\hat{\Delta}_{c}&0\\\
0&-\hat{\Delta}_{c}^{\dagger}&-H_{\downarrow}({\boldsymbol{k}})&0\\\
\hat{\Delta}_{c}^{\dagger}&0&0&-H_{\uparrow}({\boldsymbol{k}})\end{pmatrix},$
(46)

where $\Delta_{c}=I_{9}\otimes\Delta_{0}$. Here we also used the time-reversal
symmetry relations to write down the bottom right block. In the special case
$U_{1}=U_{2}=U_{3}$, $\Delta_{c}$ becomes proportional to the identity matrix,
$\Delta_{c}=\Delta\times I_{27}$, where
$\Delta=\Delta_{1}+\Delta_{2}+\Delta_{3}$. In this case, in the band basis
there is no inter-band pairings. Indeed, let

$\displaystyle
U_{\boldsymbol{k}}=\begin{pmatrix}U_{{\boldsymbol{k}}\uparrow}&0\\\
0&U_{{\boldsymbol{k}}\downarrow}\end{pmatrix}$ (47)

be the unitary operator which rotates $H({\boldsymbol{k}})$ into the band
basis ${\boldsymbol{a}}_{\boldsymbol{k}}=U{\boldsymbol{c}}_{\boldsymbol{k}}$.
Then

$\displaystyle\begin{pmatrix}0&\hat{\Delta}_{c}\\\
-\hat{\Delta}_{c}&0\end{pmatrix}\rightarrow\begin{pmatrix}U_{{\boldsymbol{k}}\uparrow}&0\\\
0&U_{{\boldsymbol{k}}\downarrow}\end{pmatrix}\begin{pmatrix}0&\hat{\Delta}_{c}\\\
-\hat{\Delta}_{c}&0\end{pmatrix}\begin{pmatrix}U_{-{\boldsymbol{k}}\uparrow}^{T}&0\\\
0&U_{-{\boldsymbol{k}}\downarrow}^{T}\end{pmatrix}=\begin{pmatrix}0&U_{{\boldsymbol{k}}\uparrow}\hat{\Delta}_{c}U_{{\boldsymbol{k}}\uparrow}^{\dagger}\\\
-U_{{\boldsymbol{k}}\downarrow}\hat{\Delta}_{c}U_{{\boldsymbol{k}}\downarrow}^{\dagger}&0\end{pmatrix}=\begin{pmatrix}0&\hat{\Delta}_{c}\\\
-\hat{\Delta}_{c}&0\end{pmatrix}.$ (48)

In the band basis, the BdG Hamiltonian takes form

$\displaystyle
H_{\mathrm{BdG}}({\boldsymbol{k}})=\frac{1}{2}\begin{pmatrix}\hat{\epsilon}_{\uparrow}({\boldsymbol{k}})&0&0&\hat{\Delta}_{c}\\\
0&\hat{\epsilon}_{\downarrow}({\boldsymbol{k}})&-\hat{\Delta}_{c}&0\\\
0&-\hat{\Delta}_{c}^{\dagger}&-\hat{\epsilon}_{\downarrow}({\boldsymbol{k}})&0\\\
\hat{\Delta}_{c}^{\dagger}&0&0&-\hat{\epsilon}_{\uparrow}({\boldsymbol{k}})\end{pmatrix},$
(49)

where $\hat{\epsilon}_{\uparrow/\downarrow}({\boldsymbol{k}})$ is the diagonal
matrix of eigenvalues $\epsilon_{i\uparrow/\downarrow}({\boldsymbol{k}})$
corresponding to spin up/down. Permuting the basis, we can rewrite the BdG
Hamiltonian in a more convenient way

$\displaystyle
H_{\mathrm{BdG}}({\boldsymbol{k}})=\frac{1}{2}\begin{pmatrix}\hat{\epsilon}_{\uparrow}({\boldsymbol{k}})&\hat{\Delta}_{c}&0&0\\\
\hat{\Delta}_{c}^{\dagger}&-\hat{\epsilon}_{\uparrow}({\boldsymbol{k}})&0&0\\\
0&0&\hat{\epsilon}_{\downarrow}({\boldsymbol{k}})&-\hat{\Delta}_{c}\\\
0&0&-\hat{\Delta}_{c}^{\dagger}&-\hat{\epsilon}_{\downarrow}({\boldsymbol{k}})\end{pmatrix}$
(50)

Squaring the Hamiltonian in Eq. 50, we can find the eigen-energies

$\displaystyle
E_{ns}({\boldsymbol{k}})^{2}=\epsilon_{ns}({\boldsymbol{k}})^{2}+\Delta^{2},\
\ (n=1,...,27).$ (51)

To find the eigenvectors, we have to solve

$\displaystyle\begin{pmatrix}\hat{\epsilon}_{s}({\boldsymbol{k}})&s\hat{\Delta}_{c}\\\
s\hat{\Delta}_{c}&-\hat{\epsilon}_{s}({\boldsymbol{k}})\end{pmatrix}\begin{pmatrix}{\boldsymbol{u}}\\\
{\boldsymbol{v}}\end{pmatrix}=E_{n}\begin{pmatrix}{\boldsymbol{u}}\\\
{\boldsymbol{v}}\end{pmatrix},$ (52)

where $s=\pm$ correspond to spin $\uparrow,\downarrow$. The solution is

$\displaystyle u_{i}$
$\displaystyle=\frac{\epsilon_{ns}({\boldsymbol{k}})+E_{n}}{\left(\left(\epsilon_{ns}({\boldsymbol{k}})+E_{n}\right)^{2}+\Delta^{2}\right)^{\frac{1}{2}}}\delta_{in},$
(53) $\displaystyle v_{i}$
$\displaystyle=\frac{\Delta}{\left(\left(\epsilon_{ns}({\boldsymbol{k}})+E_{n}\right)^{2}+\Delta^{2}\right)^{\frac{1}{2}}}\delta_{in},\
\ (i=1,...,27).$

We can write the transformation to the basis in which the BdG Hamiltonian Eq.
50 is diagonal in the form

$\displaystyle a_{{\boldsymbol{k}}is}$
$\displaystyle=u_{{\boldsymbol{k}}is}^{j\tau}\gamma_{{\boldsymbol{k}}j\tau}+v_{{\boldsymbol{k}}is}^{j\tau}\gamma_{-{\boldsymbol{k}}j\tau}^{\dagger},$
(54) $\displaystyle a_{-{\boldsymbol{k}}is}^{\dagger}$
$\displaystyle=v_{-{\boldsymbol{k}}is}^{j\tau*}\gamma_{{\boldsymbol{k}}j\tau}+u_{-{\boldsymbol{k}}is}^{j\tau*}\gamma_{-{\boldsymbol{k}}j\tau}^{\dagger},$

with

$\displaystyle u_{{\boldsymbol{k}}is}^{j\tau}$
$\displaystyle=\frac{\epsilon_{is}({\boldsymbol{k}})+E_{is}}{\left(\left(\epsilon_{is}({\boldsymbol{k}})+E_{is}\right)^{2}+\Delta^{2}\right)^{\frac{1}{2}}}\delta_{s\tau}\delta_{ij},$
(55) $\displaystyle v_{{\boldsymbol{k}}is}^{j\tau}$
$\displaystyle=\frac{\Delta}{\left(\left(\epsilon_{i\tau}(-{\boldsymbol{k}})+E_{i\tau}\right)^{2}+\Delta^{2}\right)^{\frac{1}{2}}}(1-2\delta_{\downarrow,\tau})\delta_{s\bar{\tau}}\delta_{ij},$

where $i,j=1,...,27$ are band indices, and $s,\tau$ are spin indices. In the
next subsection we proceed calculating $T_{c}$ in preparation for the final
goal – the ratio of ${H_{c2}}/{\Delta^{2}}$.

#####  II.2.2 Self-consistent determiniation of $T_{c}$

The self-consistent equations from which $T_{c}$ is determined are

$\displaystyle\Delta_{\alpha}=-U_{\alpha}\sum_{{\boldsymbol{k}}\in\mathrm{BZ}}\expectationvalue{c_{\alpha-{\boldsymbol{k}}\downarrow}c_{\alpha{\boldsymbol{k}}\uparrow}}=-U_{\alpha}\sum_{i}\sum_{{\boldsymbol{k_{i}}}\in\mathrm{lBZ_{i}}}\expectationvalue{c_{\alpha-{\boldsymbol{k_{i}}}\downarrow}c_{\alpha{\boldsymbol{k_{i}}}\uparrow}},\
\ \alpha=1,2,3.$ (56)

One can show that

$\displaystyle\Delta=\sum_{\alpha}\Delta_{\alpha}=-U_{1}\sum_{{\boldsymbol{k}}\in\mathrm{lBZ}}\sum_{i=1}^{27}\expectationvalue{a_{-{\boldsymbol{k}}i\downarrow}a_{{\boldsymbol{k}}i\uparrow}}.$
(57)

Using Eq. 54, we find

$\displaystyle-
U_{1}\sum_{{\boldsymbol{k}}\in\mathrm{lBZ}}\expectationvalue{a_{-{\boldsymbol{k}}i\downarrow}a_{{\boldsymbol{k}}i\uparrow}}=U_{1}\sum_{{\boldsymbol{k}}}\expectationvalue{\sum_{j\tau}v_{{\boldsymbol{k}}i\uparrow}^{j\tau}u_{-{\boldsymbol{k}}i\downarrow}^{j\tau}\gamma_{-{\boldsymbol{k}}j\tau}^{\dagger}\gamma_{-{\boldsymbol{k}}j\tau}+\sum_{j\tau}u_{{\boldsymbol{k}}i\uparrow}^{j\tau}v_{-{\boldsymbol{k}}i\downarrow}^{j\tau}\gamma_{{\boldsymbol{k}}j\tau}\gamma_{{\boldsymbol{k}}j\tau}^{\dagger}}.$
(58)

Using
$\expectationvalue{\gamma_{{\boldsymbol{k}}j\tau}^{\dagger}\gamma_{{\boldsymbol{k}}j\tau}}=n_{F}(E_{j\tau}({\boldsymbol{k}}))$
and
$\\{a_{{\boldsymbol{k}}is},a_{{\boldsymbol{k^{\prime}}}j\tau}^{\dagger}\\}=\delta_{ij}\delta_{s\tau}\delta_{{\boldsymbol{k}}{\boldsymbol{k^{\prime}}}}$,
we find

$\displaystyle\Delta=U_{1}\sum_{{\boldsymbol{k}},ij\tau}v_{{\boldsymbol{k}}i\uparrow}^{j\tau}u_{-{\boldsymbol{k}}i\downarrow}^{j\tau}\left(n_{F}(E_{j\bar{\tau}}({\boldsymbol{k}}))+n_{F}(E_{j\tau}({\boldsymbol{k}}))-1\right).$
(59)

And eventually

$\displaystyle\Delta=U_{1}\sum_{{\boldsymbol{k}},i}\frac{\Delta}{2E_{i\downarrow}({\boldsymbol{k}})}\left(1-n_{F}(E_{i\uparrow}({\boldsymbol{k}}))-n_{F}(E_{i\downarrow}({\boldsymbol{k}})\right).$
(60)

At $T_{c}$, this equation reduces to

$\displaystyle\frac{1}{U_{1}}=\sum_{{\boldsymbol{k}},i}\frac{1}{2\absolutevalue{\epsilon_{i\downarrow}({\boldsymbol{k}})}}\left(1-n_{F}(\absolutevalue{\epsilon_{i\uparrow}({\boldsymbol{k}})})-n_{F}(\absolutevalue{\epsilon_{i\downarrow}({\boldsymbol{k}})})\right).$
(61)

Note that energies $\epsilon_{is}$ are measured with respect to the chemical
potential.

We compute the dependence of $T_{c}$ on the CDW strength $b_{1}$ numerically.
These dependencies are plotted in Fig. S7.

(a)

(b)

Figure S7: $T_{c}$ as a function of CDW strength $b_{1}$ for (a) NbSe2 and (b)
TaS2.

Mainly, for a fixed value of $U_{1}$ the value of $T_{c}$ is controlled by the
density of states $\mathcal{N}(\epsilon_{F})$. The increase of $T_{c}$ for
small values of $b_{1}$ happens due to the chemical potential, $\mu$, falls
onto the van Hove singularity developing due to CDW order. The dependence of
$\mu$ on the CDW strength $b_{1}$ is plotted in Fig. S8, and the density of
states for some values of $b_{1}$ are plotted in Figs. S9 and S10 for NbSe2
and TaS2, respectively. In Fig. S11, we plot the smeared DOS for NbSe2 and
TaS2.

(a)

(b)

Figure S8: $\mu$ as a function of CDW strength $b_{1}$ for (a) NbSe2 and (b)
TaS2.

(a)

(b)

Figure S9: DOS for NbSe2 for CDW strength (a) $b_{1}=0$ eV and (b)
$b_{1}=0.25$ eV.

(a)

(b)

Figure S10: DOS for TaS2 for CDW strength (a) $b_{1}=0$ eV and (b) $b_{1}=0.3$
eV.

(a)

(b)

Figure S11: Smeared DOS for (a) NbSe2 at $b_{1}=0$ eV (blue) and $b_{1}=0.25$
eV (yellow) and for (b) TaS2 at $b_{1}=0$ eV (yellow) and $b_{1}=0.3$ eV
(blue).

####  II.3 Superconductivity – Ginzburg-Landau theory

Our goal is to compute the ratio of $H_{c2}/\Delta_{0}^{2}$ at the temperature
at which the experimental measurements were conducted (which is essentially
$T=0$). However, we will do this only approximately: Our line of attack is to
compute $H_{c2}/\Delta^{2}$ at $T_{c}$ from the Ginzburg-Landau theory

$H_{c2}(T)\sim{\phi_{0}\over
2\pi\xi_{GL}^{2}}(1-T/T_{c})\,\,,\,\Delta(T)\sim\Delta_{0}\sqrt{1-T/T_{c}}\;,\;T\to
T_{c}$

and extrapolate to zero temperature. Namely, estimate the above equation at
$T=0$, $H_{c2}(0)$. In actuality, this is an upper bound on this ratio, as it
curves and saturates to a smaller value.

The form of the GL theory, which is consistent with the point group of the
system, is given by

$\displaystyle S=\int
d^{d}xd\tau\left(\frac{1}{2\eta}\absolutevalue{\left(-i\gradient-2eA\right)\Delta}^{2}+r(T)\absolutevalue{\Delta}^{2}+\frac{u(T)}{2}\absolutevalue{\Delta}^{4}+...\right)$
(62)

Having $H_{c2}={\phi_{0}}/{2\pi\xi^{2}}$, $\Delta^{2}=-{r(T)}/{u(T)}$,
$\xi^{2}=-{1}/{2\eta r(T)}=\xi_{GL}^{2}(1-T/T_{c})$, we find that the ratio of
the upper-critical field to $\Delta^{2}$, in terms of the GL parameters, is
given by

$\displaystyle\frac{H_{c2}}{\Delta^{2}}=\frac{\phi_{0}}{\pi}\eta u.$ (63)

Thus, our next (and final) goal is to compute these two coefficients from the
microscopic model.

#####  II.3.1 Quadratic terms:

The coefficient of the quadratic terms in Eq. 62 are obtained from

$\displaystyle\Tr\left(G_{0}\hat{\Delta}\right)^{2}=-2\Tr\left(G_{e\uparrow{\boldsymbol{p}}}\hat{\Delta}_{c}G_{h\downarrow{\boldsymbol{p-k}}}\hat{\Delta}_{c}+G_{e\downarrow{\boldsymbol{p}}}\hat{\Delta}_{c}G_{h\uparrow{\boldsymbol{p-k}}}\hat{\Delta}_{c}\right),$
(64)

where

$\displaystyle
G_{es{\boldsymbol{p}}}=\left(-i\omega+H_{\mathrm{CDW},s}({\boldsymbol{p}})\right)^{-1},\
G_{hs{\boldsymbol{p}}}=\left(-i\omega-H^{T}_{\mathrm{CDW},s}(-{\boldsymbol{p}})\right)^{-1},$
(65)

where $s=\uparrow,\downarrow$, and $G_{0}$ is the normal part of the Green’s
function of the BDG Hamiltonian

$\displaystyle G_{0}=2G_{e}\oplus G_{h},$ (66)

where

$\displaystyle G_{e}=G_{e\uparrow}\oplus G_{e\downarrow},\
G_{h}=G_{h\uparrow}\oplus G_{h\downarrow}.$ (67)

Recall that for $U_{1}=U_{2}=U_{3}$, $\hat{\Delta}_{c}=I_{27}\otimes\Delta$.
Then Eq. 64 simplifies to

$\displaystyle\Tr\left(G_{0}\hat{\Delta}\right)^{2}=-2\Tr\left(G_{e\uparrow{\boldsymbol{p}}}G_{h\downarrow{\boldsymbol{p-k}}}+G_{e\downarrow{\boldsymbol{p}}}G_{h\uparrow{\boldsymbol{p-k}}}\right)\Delta^{2}.$
(68)

Denoting
$\Pi_{s}({\boldsymbol{k}})=\Tr\left(G_{es{\boldsymbol{p}}}G_{h\bar{s}{\boldsymbol{p-k}}}\right)$
and switching to the band-basis (recall that $\hat{\Delta}_{c}$ is invariant
under this transformation in the particular case we consider), we can write

$\displaystyle\Pi_{s}({\boldsymbol{k}})=\frac{T}{L^{3}}\sum_{\omega_{n}}\sum_{{\boldsymbol{p}}}\sum_{i}\frac{1}{-i\omega_{n}+\epsilon_{i{\boldsymbol{p}}s}}\frac{1}{-i\omega_{n}-\epsilon_{i{\boldsymbol{p}}-{\boldsymbol{k}}s}}\approx\frac{T}{L^{3}}\sum_{\omega_{n}}\sum_{{\boldsymbol{p}}}\sum_{i}\frac{1}{-i\omega_{n}+\epsilon_{i{\boldsymbol{p}}s}}\frac{1}{-i\omega_{n}-\epsilon_{i{\boldsymbol{p}}s}+{\boldsymbol{v}}_{i}({\boldsymbol{p}})\cdot{\boldsymbol{k}}},$
(69)

and

$\displaystyle\frac{1}{\eta}=-\left(\frac{\partial^{2}\Pi_{\uparrow}({\boldsymbol{0}})}{\partial
k_{x}^{2}}+\frac{\partial^{2}\Pi_{\downarrow}({\boldsymbol{0}})}{\partial
k_{x}^{2}}\right).$ (70)

Expanding in small ${\boldsymbol{k}}$ and switching from the sum to the
integration over the momentum

$\displaystyle\frac{1}{2}\frac{\partial^{2}\Pi_{s}({\boldsymbol{0}})}{\partial
k_{x}^{2}}\approx\frac{1}{2}T\sum_{i}\sum_{\omega_{n}}\int
pdpd\theta\frac{1}{-i\omega_{n}+\epsilon_{i{\boldsymbol{p}}s}}\frac{v_{i}(p,\theta)^{2}\cos^{2}\alpha_{i}(\theta)}{\left(-i\omega_{n}-\epsilon_{i{\boldsymbol{p}}s}\right)^{3}},$
(71)

where $\alpha_{i}(\theta)$ is the angle between ${\boldsymbol{v}}_{i}$ and
${\boldsymbol{k}}$.

One can show that
$\frac{1}{2}\frac{\partial^{2}\Pi_{s}({\boldsymbol{k}})}{\partial
k_{x}^{2}}=\frac{1}{2}\frac{\partial^{2}\Pi_{s}({\boldsymbol{k}})}{\partial
k_{y}^{2}}$, and then we can write

$\displaystyle\frac{\partial^{2}\Pi_{s}({\boldsymbol{0}})}{\partial
k_{x}^{2}}=\frac{1}{2}\frac{\partial^{2}\Pi_{s}({\boldsymbol{0}})}{\partial
k_{x}^{2}}+\frac{1}{2}\frac{\partial^{2}\Pi_{s}({\boldsymbol{0}})}{\partial
k_{y}^{2}}=\frac{1}{2}T\sum_{i}\sum_{\omega_{n}}\int
pdpd\theta\frac{1}{-i\omega_{n}+\epsilon_{i{\boldsymbol{p}}s}}\frac{v_{i}(p,\theta)^{2}}{\left(-i\omega_{n}-\epsilon_{i{\boldsymbol{p}}s}\right)^{3}}.$
(72)

Assuming $v_{i}(p,\theta)$ does not depend much on the absolute value of the
momentum ${\boldsymbol{p}}$ such that $v_{i}(p,\theta)\approx v_{i}(\theta)$
is a good approximation, we can rewrite the last expression using
$\theta$-dependent density of states per angle introduced as follows

$\displaystyle\int d^{2}{\boldsymbol{k}}f(\epsilon_{\boldsymbol{k}},\theta)$
$\displaystyle=\int kdk\int d\theta f(\epsilon_{\boldsymbol{k}},\theta)=\int
kdk\int d\theta\int d\epsilon\delta(\epsilon-\epsilon_{k})f(\epsilon,\theta)$
(73) $\displaystyle=\int d\epsilon\int d\theta\int kdk\
\delta(\epsilon-\epsilon_{k})f(\epsilon,\theta)=\int d\epsilon\int
d\theta\mathcal{N}_{\theta}(\epsilon)f(\epsilon,\theta),$

where $f(\epsilon_{\boldsymbol{k}},\theta)$ is an arbitrary function. The
approximation $v_{i}(p,\theta)\approx v_{i}(\theta)$ is not valid in the
vicinity of a Van Hove singularity.

Thus

$\displaystyle\frac{\partial^{2}\Pi_{s}({\boldsymbol{0}})}{\partial
k_{x}^{2}}$
$\displaystyle=\frac{1}{2}T\sum_{i}\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon_{i}\int
d\theta\frac{\mathcal{N}_{is\theta}(\epsilon_{i})}{-i\omega_{n}+\epsilon_{i}}\frac{v_{i}(\theta)^{2}}{\left(-i\omega_{n}-\epsilon_{i}\right)^{3}}\approx\frac{1}{2}T\sum_{i}\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon_{i}\int
d\theta\frac{\mathcal{N}_{is\theta}(0)}{-i\omega_{n}+\epsilon_{i}}\frac{v_{i}(\theta)^{2}}{\left(-i\omega_{n}-\epsilon_{i}\right)^{3}}$
(74)
$\displaystyle=\frac{1}{2}T\sum_{i}\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon_{i}\frac{\mathcal{N}_{is}(0)}{-i\omega_{n}+\epsilon_{i}}\frac{\mkern
1.5mu\overline{\mkern-1.5muv_{i}^{2}\mkern-1.5mu}\mkern
1.5mu}{\left(-i\omega_{n}-\epsilon_{i}\right)^{3}}=\frac{1}{2}T\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon\frac{\mathcal{N}_{s}(0)}{-i\omega_{n}+\epsilon}\frac{\mkern
1.5mu\overline{\mkern-1.5muv^{2}\mkern-1.5mu}\mkern
1.5mu}{\left(-i\omega_{n}-\epsilon\right)^{3}},$

where $\mathcal{N}_{is\theta}$ is the DOS per angle associated with the band
$i$ and spin $s$, and we introduced average velocities associated with band
$i$ via $\int
d\theta\mathcal{N}_{is\theta}(0)v_{i}(\theta)^{2}=\mathcal{N}_{is}(0)\mkern
1.5mu\overline{\mkern-1.5muv_{i}^{2}\mkern-1.5mu}\mkern 1.5mu$ and overall
average velocity $\sum_{i}\mathcal{N}_{is}(0)\mkern
1.5mu\overline{\mkern-1.5muv_{i}^{2}\mkern-1.5mu}\mkern
1.5mu=\mathcal{N}_{s}(0)\mkern
1.5mu\overline{\mkern-1.5muv^{2}\mkern-1.5mu}\mkern 1.5mu$. We can proceed
calculating the integral and eventually obtaining the Gorkov’s result, but
$\mkern 1.5mu\overline{\mkern-1.5muv^{2}\mkern-1.5mu}\mkern 1.5mu$ has to be
computed numerically.

Eq. 74 can be used for gaining some insights into the effect of CDW if Van
Hove singularities are sufficiently far away from the chemical potential
level. In any case, we compute $\eta$ via numerically differentiating Eq. 69.

#####  II.3.2 Quartic term:

The quartic term is obtained from the static uniform part, $S$, of

$\displaystyle\Tr\left(G_{0}\hat{\Delta}\right)^{4}$
$\displaystyle=\Tr\\{G_{e{\boldsymbol{p}}}\hat{\Delta}({\boldsymbol{k_{1}}})G_{h{\boldsymbol{p}}-{\boldsymbol{k_{1}}}}\hat{\Delta}^{\dagger}(k_{2}))G_{e{\boldsymbol{p}}-{\boldsymbol{k_{1}}}-{\boldsymbol{k_{2}}}}\hat{\Delta}({\boldsymbol{k_{3}}})G_{h{\boldsymbol{p}}-{\boldsymbol{k_{1}}}-{\boldsymbol{k_{2}}}-{\boldsymbol{k_{3}}}}\hat{\Delta}^{\dagger}(k_{4}))$
(75)
$\displaystyle+G_{h{\boldsymbol{p}}}\hat{\Delta}({\boldsymbol{k_{1}}})G_{e{\boldsymbol{p}}-{\boldsymbol{k_{1}}}}\hat{\Delta}^{\dagger}(k_{2}))G_{h{\boldsymbol{p}}-{\boldsymbol{k_{1}}}-{\boldsymbol{k_{2}}}}\hat{\Delta}({\boldsymbol{k_{3}}})G_{e{\boldsymbol{p}}-{\boldsymbol{k_{1}}}-{\boldsymbol{k_{2}}}-{\boldsymbol{k_{3}}}}\hat{\Delta}^{\dagger}(k_{4}))\bigg{\\}}.$

Using band basis, we find

$\displaystyle
S=2\Tr\left\\{G_{e\uparrow{\boldsymbol{p}}}G_{h\downarrow{\boldsymbol{p}}}G_{e\uparrow{\boldsymbol{p}}}G_{h\downarrow{\boldsymbol{p}}}+G_{e\downarrow{\boldsymbol{p}}}G_{h\uparrow{\boldsymbol{p}}}G_{e\downarrow{\boldsymbol{p}}}G_{h\uparrow{\boldsymbol{p}}}\right\\}\absolutevalue{\Delta}^{4}.$
(76)

Denoting
$\Sigma_{s}=\Tr\left\\{G_{es{\boldsymbol{p}}}G_{h\bar{s}{\boldsymbol{p}}}G_{es{\boldsymbol{p}}}G_{h\bar{s}{\boldsymbol{p}}}\right\\}$,
we have

$\displaystyle u=4(\Sigma_{\uparrow}+\Sigma_{\downarrow})=8\Sigma_{\uparrow},$
(77)

where the last equality holds due to TRS.

In the band basis

$\displaystyle\Sigma_{s}$
$\displaystyle=\frac{T}{L^{3}}\sum_{i}\sum_{\omega_{n}}\sum_{{\boldsymbol{p}}}\frac{1}{\left(-i\omega_{n}+\epsilon_{i{\boldsymbol{p}}s}\right)^{2}}\frac{1}{\left(-i\omega_{n}-\epsilon_{i{\boldsymbol{p}}s}\right)^{2}}$
(78)
$\displaystyle=\frac{T}{L^{3}}\sum_{i}\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon_{is}\frac{\mathcal{N}_{is}(\epsilon_{is})}{\left(-i\omega_{n}+\epsilon_{is}\right)^{2}}\frac{1}{\left(-i\omega_{n}-\epsilon_{is}\right)^{2}}=\frac{T}{L^{3}}\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon\frac{\mathcal{N}_{s}(\epsilon)}{\left(-i\omega_{n}+\epsilon\right)^{2}}\frac{1}{\left(-i\omega_{n}-\epsilon\right)^{2}}$
$\displaystyle\approx\frac{T}{L^{3}}\sum_{\omega_{n}}\int_{-\omega_{D}}^{\omega_{D}}d\epsilon\frac{\mathcal{N}_{s}(0)}{\left(-i\omega_{n}+\epsilon\right)^{2}}\frac{1}{\left(-i\omega_{n}-\epsilon\right)^{2}}.$

Completing this integral and obtaining the Gor’kov’s result, we compute $u$
numerically from the first line of equation Eq. 78. Finally, we plot the
computed $H_{c2}/\Delta^{2}$ for NbSe2 and TaS2 in Fig. S12.

Figure S12: $H_{c2}/\Delta^{2}$ vs $T_{c}$ for TaS2 (blue) and NbSe2 (orange).

####  II.4 Details of computations

In numerical computations, we utilize $C_{3}$ symmetry of the problem and
generate Monkhorst-Pack grid Monkhorst and Pack (1976) for the one third of
the lBZ with $L(L+1)$ points. For computing of all values, we used $L=600$.

The DOS were calculated by counting number of states $N_{E}$ in the intervals
of energy $(E,E+\Delta E)$ with $\Delta E=0.001$ eV from $E_{0}=-1$ eV to
$E_{f}=1$ eV. The smeared DOS were calculated from the computed ones by
smearing the delta functions in
$N_{E}=N_{E}\int\delta\left(\epsilon-(E+\frac{\Delta E}{2})\right)d\epsilon$
into normal distribution functions with dispersion $\sigma=0.01$ eV.

The shift of the chemical potential for non-zero CDW strength was calculated
by preserving number of occupied states on the same grid with $L=600$.

In computing critical temperatures, we set $T_{c}=9.8\ (3.4)$ K for NbSe2
(TaS2) for zero CDW strength, $b_{1}=0$, and compute the corresponding $U_{1}$
from Eq. 61. For a non-zero CDW strength, $b_{1}\neq 0$, we then use the
computed $U_{1}$ to find $T_{c}$ from Eq. 61.

In computing $\eta$ from Eq. 70, we use the central difference formula for the
second derivative with $\Delta k_{x}=10^{-9}\frac{3\sqrt{3}}{4}v_{1}$, where
$v_{1}$ is the length of the primitive vector in the reciprocal lattice that
defines lBZ.

### References

  * Dvir  _et al._ (2018) Tom Dvir, Freek Massee, Lotan Attias, Maxim Khodas, Marco Aprili, Charis HL Quay, and Hadar Steinberg, “Spectroscopy of bulk and few-layer superconducting nbse2 with van der waals tunnel junctions,” Nature communications 9, 598 (2018).
  * Ye  _et al._ (2012) JT Ye, Yj J Zhang, R Akashi, Ms S Bahramy, R Arita, and Y Iwasa, “Superconducting dome in a gate-tuned band insulator,” Science 338, 1193–1196 (2012).
  * Costanzo  _et al._ (2016) Davide Costanzo, Sanghyun Jo, Helmuth Berger, and Alberto F. Morpurgo, “Gate-induced superconductivity in atomically thin MoS2 crystals,” Nature Nanotechnology 11, 339–344 (2016).
  * Ribak  _et al._ (2020) A Ribak, R Majlin Skiff, M Mograbi, PK Rout, MH Fischer, J Ruhman, K Chashka, Y Dagan, and A Kanigel, “Chiral superconductivity in the alternate stacking compound 4hb-tas2,” Science advances 6, eaax9480 (2020).
  * Xi  _et al._ (2016) Xiaoxiang Xi, Zefang Wang, Weiwei Zhao, Ju-Hyun Park, Kam Tuen Law, Helmuth Berger, László Forró, Jie Shan, and Kin Fai Mak, “Ising pairing in superconducting nbse2 atomic layers,” Nature Physics 12, 139–143 (2016).
  * Yokota  _et al._ (2000) Ken-ichi Yokota, Gouo Kurata, Tomohiro Matsui, and Hiroshi Fukuyama, “Superconductivity in the quasi-two-dimensional conductor 2h-tase2,” Physica B: condensed matter 284, 551–552 (2000).
  * Tissen  _et al._ (2013) VG Tissen, MR Osorio, Jean-Pascal Brison, NM Nemes, Mar García-Hernández, Laurent Cario, Pierre Rodiere, Sebastian Vieira, and Hermann Suderow, “Pressure dependence of superconducting critical temperature and upper critical field of 2h-nbs2,” Physical Review B 87, 134502 (2013).
  * Liu  _et al._ (2013) Gui-Bin Liu, Wen-Yu Shan, Yugui Yao, Wang Yao, and Di Xiao, “Three-band tight-binding model for monolayers of group-vib transition metal dichalcogenides,” Phys. Rev. B 88, 085433 (2013).
  * Flicker and van Wezel (2016) Felix Flicker and Jasper van Wezel, “Charge order in nbse2,” Physical Review B 94, 235135 (2016).
  * Möckli and Khodas (2018) David Möckli and Maxim Khodas, “Robust parity-mixed superconductivity in disordered monolayer transition metal dichalcogenides,” Physical Review B 98, 144518 (2018).
  * Monkhorst and Pack (1976) Hendrik J Monkhorst and James D Pack, “Special points for brillouin-zone integrations,” Physical review B 13, 5188 (1976).


# Exploring Top-Quark Signatures of Heavy Flavor-Violating Scalars at the LHC
with Parametrized Neural Networks

Alexandre Alves<EMAIL_ADDRESS>Departamento de Física, Universidade
Federal de São Paulo, UNIFESP, 09972-270, Diadema-SP, Brazil  Eduardo da Silva
Almeida<EMAIL_ADDRESS>Departamento de Física do Estado Sólido,
Universidade Federal da Bahia,
UFBA, 40170-115, Salvador-BA, Brazil  Alex G. Dias<EMAIL_ADDRESS>Centro de Ciências Naturais e Humanas, Universidade Federal do ABC,
UFABC, 09210-580, Santo André-SP, Brazil  Diego S. V. Gonçalves
<EMAIL_ADDRESS>Centro de Ciências Naturais e Humanas, Universidade
Federal do ABC,
UFABC, 09210-580, Santo André-SP, Brazil

###### Abstract

In this work, we study flavor-violating scalars (flavons) in a range of large
masses that have not been explored previously. We model the interactions with
an effective field theory formulation where the flavon is heavier than the top
quark. In addition, we assume that the flavon only couples to fermions of the
Standard Model in a flavor-changing way. As the flavon couples strongly to top
quarks, same-sign and opposite-sign top quark pair signals can be explored in
the search for those particles. Using parametrized neural networks, we show
that it is possible to probe flavons with masses in the 200–1600 GeV range
through their interactions with a top quark plus up and charm quarks for
effective couplings of order $10^{-2}$ TeV-1 at the 14 TeV High-Luminosity
LHC.

##  I Introduction

Hypothetical scalar particles with flavor-violating interactions are part of
the particle spectra of several models proposed to tackle problems that do not
find a solution within the Standard Model (SM). Among the models predicting
such particles are the ones based on the Froggatt-Nielsen mechanism Froggatt
and Nielsen (1979), where the hierarchical flavor structure of the SM fermions
originates from suppressed effective Yukawa couplings given by powers of the
ratio of the condensation scale of a scalar field breaking the global flavor
$U(1)_{F}$ symmetry – one that acts differently in each fermion multiplet of
the gauge symmetry – and a high-energy scale. The particle excitation of the
pseudo-Nanbu-Goldstone boson field of this symmetry breaking is the
$\it{flavon}$ and has flavor-violating interactions. Axions and axion-like
particles Di Luzio  _et al._ (2020); Jaeckel and Ringwald (2010) can also have
flavor-violating interactions when their associated symmetry acts similarly to
a flavor symmetry Ema  _et al._ (2017); Björkeroth  _et al._ (2018); Bauer
_et al._ (2022); Jho  _et al._ (2022); Dias  _et al._ (2021). As is typical of
these particles, their flavor-violating interactions with matter are
proportional to the mass of the fermions involved but suppressed by the new
physics energy scale. Especially for axion-like particles and flavons111In the
typical, axion models the suppression scale is too high Di Luzio  _et al._
(2017). However, it is possible to constrain axion flavor interactions in
meson decays Björkeroth  _et al._ (2018)., the suppression is such that it
would still allow for signals of their production at the LHC Bauer  _et al._
(2017).

In this work, we investigate signals of heavy flavons at the LHC 14 TeV under
the hypothesis that their main interactions violate flavor (the maximally
violating flavons) with perturbative couplings proportional to the quark
masses. We explore the effective top quark-flavon couplings which are not
still severely constrained, different from the flavor-violating couplings
involving the down-type quarks. We assume that flavor-conserving couplings of
the flavon to the fermions, as well as to photons, gluons, and electroweak
gauge boson pairs, are relatively small compared to the considered flavor-
violating ones.

Our analysis focuses on flavon masses higher than 200 GeV, which decay
predominantly into a top quark plus an up/charm quark. We identify a peculiar
phenomenology of these new particles where final states contain both same-sign
and opposite-sign top quark pairs, impacting the signal/background separation.
Same-charge tops signals from flavor-violating scalars in two-Higgs doublets
models have been searched by ATLAS Aad  _et al._ (2023) and CMS Hayrapetyan
_et al._ (2024).

Machine learning (ML) techniques Mehta  _et al._ (2019); Carleo  _et al._
(2019); Bourilkov (2020); Radovic  _et al._ (2018) are powerful tools in new
physics searches where signal-to-background ratios are expected to be too
small to produce statistically significant signals. We chose to employ a
parametrized neural network (pNN) Baldi  _et al._ (2016); Anzalone  _et al._
(2022) for the signal-background classification task to boost the discovery
and exclusion limits of our effective field theory (EFT). The pNN option is
adequate if we want the neural networks to learn to separate signals from
backgrounds as a continuous function of the flavon mass parameter.

Our main result is that flavons in the considered mass interval can be
discovered, or excluded, in the 14 TeV High-Luminosity LHC (HL-LHC) for
effective couplings up to order $10^{-2}$ TeV-1, extending considerably the
search prospects for these particles. The findings of the present analysis
represent a complementary study of previous works Alves  _et al._ (2024);
Cornella  _et al._ (2020) for flavons lighter than the top quark.

In section II, we present the effective field theory formulation of the
interactions; in section III, the signal and background simulations are
discussed; the analysis strategy and the neural network training are discussed
in section IV; our results are presented in section V, while the conclusions
can be found in section VI.

##  II Effective Lagrangian

We assume the flavon field originates from the spontaneous symmetry breaking
of a global $U(1)_{F}$ flavor symmetry as in the Froggatt-Nielsen scheme
Froggatt and Nielsen (1979). The main interactions of the flavon with fermions
arise from effective operators like
$y_{ij}(\frac{\phi}{M})^{n_{f}}\bar{F}_{iL}\Phi f^{\prime}_{jR}$, with
$y_{ij}$ couplings that in principle could be of order one; $M$ a high-energy
scale;
$\phi(x)=(\rho(x)+\sqrt{2}\langle\phi\rangle)e^{ia(x)/\Lambda}/\sqrt{2}$ a
scalar singlet field of the SM, hosting the flavon field $a$ and having vacuum
expectation value $\langle\phi\rangle$ breaking the $U(1)_{F}$ symmetry;
$F_{iL}$ the fermion doublet of left-handed fields and $f^{\prime}_{jR}$ the
right-handed singlet fermion fields; $\Phi=H$ (or $i\sigma_{2}H^{*}$ depending
on the fermion $f^{\prime}_{jR}$) the Higgs doublet, and $n_{f}$ integer
numbers fixed by the flavor symmetry. The scale $\Lambda$ may depend on the
vacuum expectation values of other scalar fields charged under
$U(1)_{F}$222This is, in fact, the flavon decay constant and can depend on the
$U(1)_{F}$ charge $c_{i}$ of the scalar fields $\varphi_{i}$ getting vacuum
expectation value as
$\Lambda=\sqrt{2(\sum_{i}c_{i}^{2}\langle\varphi_{i}\rangle^{2})}$.. With the
vacuum expectation values of $H$ and $\phi$, we have that the flavon
interactions with SM fermions are described by the effective Lagrangian (e.g.
Ema  _et al._ (2017))

$\mathcal{L}\supset-i\frac{a}{\Lambda}\sum_{f,i,j}\overline{f_{i}}\left[(m_{j}-m_{i})v_{ij}+(m_{j}+m_{i})a_{ij}\gamma_{5}\right]f_{j}\,.$
(1)

where $m_{i}$ is the mass of the mass eigenstate fermion field $f_{i}$;
$a_{ij}$ and $v_{ij}$ the axial and vectorial coefficients, respectively.

Depending on the coefficients $a_{ij}$ and $v_{ij}$, given the mass
proportionality in Eq. (1), the most intense interactions of the flavon would
involve the charged fermions of the third generation. Since $a$ is a
pseudoscalar field, the flavor-violating (FV) couplings $v_{ij}$ imply parity
symmetry violation. However, the CP symmetry is preserved if $v_{ij}$ and
$a_{ij}$ are real, which is the case we will study.

It is also supposed that the flavon has a mass that could be generated by some
explicitly breaking mechanism of the $U(1)_{F}$ symmetry through an anomalous
feature, as in the case of the axion, or directly by an potential term like
$V=({\lambda}\phi^{n}/{M^{n-4}}+h.c.)$. For our purposes the origin of this
mechanism is not relevant so that we simply assume a mass $m_{a}$ for the
flavon.

Our analysis will be based on the following hypotheses concerning the
couplings: the flavon diagonal couplings with fermions, vector bosons, as well
as the non-diagonal couplings with the down-type quarks are all negligible;
the dominant couplings with quarks are the non-diagonal ones involving the top
quark plus charm or up quark; the couplings with leptons play a subdominant
role as the flavon couples to the fermion mass. This way, the set of
parameters we consider is the following:

$\left\\{\begin{array}[]{cc}m_{a},\Lambda&\hbox{flavon mass and new physics
scale}\\\ v_{tc},a_{tc},v_{tu},a_{tu}&\hbox{non-diagonal top-quark
couplings}\end{array}\right.$ (2)

These couplings avoid dangerous flavor-changing neutral currents at tree level
and strong experimental constraints on flavor-diagonal interactions.

To restrict the parameters space further, we will also assume universal flavon
couplings to quarks, that is, $v_{tc}=v_{tu}\equiv v_{tq}$, and
$a_{tc}=a_{tu}\equiv a_{tq}$. Moreover, the cross sections will depend only on
the ratio between $c_{tq}\equiv\sqrt{v_{tq}^{2}+a_{tq}^{2}}$ and $\Lambda$,
thus our searches can be performed on the $c_{tq}/\Lambda$ versus the flavon
mass, $m_{a}$.

Our analysis is restricted to flavons with masses above 200 GeV, which prefer
to decay into top plus up or charm quarks. Even though the effective flavon
couplings to down-type quarks are assumed to be negligibly small, flavons can
decay to $a\rightarrow bq$ with $q=\bar{d},\bar{s}$ at 1-loop. However, this
decay is suppressed compared to the top-quark decay channel since it arises
from next-to-next leading order corrections Alves  _et al._ (2024).

##  III SIGNAL AND BACKGROUND SIMULATIONS

We simulated signals and backgrounds in leading order at the 14 TeV LHC with
MadGraph5 Alwall  _et al._ (2014); Frederix  _et al._ (2018) using the UFO
Degrande  _et al._ (2012) files generated by FeynRules Alloul  _et al._
(2014), where our EFT was implemented. The parton shower and detector
simulations were performed with default settings of PYTHIA 8.3 Bierlich  _et
al._ (2022) and Delphes3 de Favereau  _et al._ (2014), respectively.

$a$$q$$t$$g$$t$$a$$q$$a$$\bar{q}^{\prime}$$t$$t$$q$$t$$\bar{q}^{\prime}$$a$

Figure 1: Feynman diagrams for top pair signals. The left and middle diagrams
depict the contributions from single and double production of flavons,
respectively, while the rightmost diagram shows the virtual flavon
contribution in the $t$-channel. Here, $q=u,c$.

The signal we consider involves the production of 2 top quarks with the same
or opposite charges and up to 2 non-$b$ jets after the decay of flavons into a
top quark plus a $c$ or $u$ quark. We assume that flavons with masses higher
than $m_{t}+m_{c}$ decay exclusively and equally into the dominating flavor-
violating channels, $a\to t\bar{u},t\bar{c},\bar{t}u,\bar{t}c$. The signal
receives the following contributions:

  * (a)

flavon in the $t$-channel, as shown at the right panel of Fig. 1. In this
case, two top quarks of any combination of charges are produced, $tt$,
$\bar{t}\bar{t}$, $t\bar{t}$. After decay, the top quarks produce two $b$ jets
and same-sign (SS) and opposite-sign (OS) leptons ($\ell=e,\mu$);

  * (b)

flavon single production, as shown in the left panel of Fig. 1. Again, we have
OS and SS top quark pairs but with a non-$b$ jet from the flavon decay;

  * (c)

double flavon production as shown in the middle panel of Fig. 1. Now, we have
OS and SS top quark pairs along with two non-$b$ jets from the flavon decay.

Summarizing, we propose to look for heavy flavons in the channels

$\displaystyle pp\to tt,\;\bar{t}\bar{t}+\hbox{jets}\to
bb(\bar{b}\bar{b})+\ell^{\pm}\ell^{\pm}+\not\\!\\!E_{T}+\hbox{jets}$ (3)
$\displaystyle pp\to t\bar{t}+\hbox{jets}\to
b\bar{b}+\ell^{\pm}\ell^{\mp}+\not\\!\\!E_{T}+\hbox{jets}\;.$ (4)

The production cross sections, in pb, for all three channels are shown in Fig.
2. For the range of masses we are interested in, the single production
dominates, followed by the $t$-channel contribution and the double flavon
production, which is suppressed by the mass of the flavon.

The main background sources for the OS leptons signals are $t\bar{t}$ as the
dominating one, next-to-next-to-leading order Drell-Yan plus two $b$ jets
($\ell\ell bb$), and single top $tW$ as subdominant. Same-sign lepton signals
have small reducible backgrounds from misidentifying leptons and jets from the
three background sources described. The only irreducible background is the
vector boson fusion (VBF) channel $qq^{\prime}\to W^{\pm}W^{\pm}+bb$, which is
CKM suppressed. Moreover, in this case, the bottom jets are produced in the
forward regions of the detector, where $b$-tagging is less efficient, while
the signals tend to be produced in the central regions of the detector. For
these reasons, we neglect VBF events but collect all events with two SS
leptons from the background sources described previously.

Figure 2: Production cross sections of all the three production modes
considered in this work: pair production, $pp\to aa$; single production,
$pp\to at$; $t$-channel contribution, $pp\to tt$. In this last case, all
combinations of top quark charges were included.

##  IV ANALYSIS Strategy AND PARAMETRIZED NEURAL NETWORKs

###  IV.1 Basic Requirements

Before training the neural networks, we impose kinematic cuts on both SS and
OS signals and backgrounds. We exclude particles that do not meet the
following selection

$\displaystyle\text{at least 1 $b$-jet}\quad\text{and}\quad\text{exactly 2
leptons},$ $\displaystyle
p_{T_{\ell}}>10\;\mbox{GeV},\;p_{T_{j,b}}>20\;\mbox{GeV},\;\left|\eta_{\ell}\right|<2.5,\;\left|\eta_{j,b}\right|<2.5,$
$\displaystyle
M_{\ell\ell}>100\;\mbox{GeV},\;\not{E}_{T}>20\;$\mathrm{G}\mathrm{e}\mathrm{V}$\;.$
(5)

Requiring at least one $b$-jet efficiently raises the overall $b$-tag to $\sim
70$%. On the other hand, requiring two leptons drops the efficiency of
signals, especially for heavy masses, as it gets more difficult to find
isolated leptons. The cut on the invariant mass of two leptons,
$M_{\ell\ell}$, helps to tame the large Drell-Yan background.

In Table 1, we display the cross sections after imposing the cuts of Eq.
(IV.1). The background rates are still much larger than the signals. We could
have selected events with same-sign leptons to clean up almost all those
backgrounds at the cost of half of the signal events, but because it is easy
to separate SS and OS events, we can use the OS events if we are able to
classify signals versus backgrounds in this case.

$t\bar{t}$ (pb) |  $tW$ (pb) |  $\ell\ell bb$ (pb) |  $200$ GeV (fb) |  $1000$ GeV (fb) |  $1400$ GeV (fb)
---|---|---|---|---|---
7.3 | 0.76 | 0.07 | 45.1 | 1.1 | 0.2
Table 1: The cross sections after cuts of Eq. (IV.1) for backgrounds, in pb,
and signals corresponding to flavon masses of 200, 1000, and 1400 GeV, in fb.
The signal cross sections were calculated with $\Lambda=1$ TeV.

As we will show, the neural networks can efficiently distinguish between OS
signals and backgrounds, while SS signals have very little contamination. Of
course, the neural networks can also easily distinguish between OS and SS
events. However, if we mix them together and frame the problem as OS+SS
signals against backgrounds, we lose statistical significance as the number of
SS signal events is much smaller than $t\bar{t}$ events even after pNN
classification. Instead, we can calculate the statistical significance of the
SS and OS signals separately after the pNN classification, combine them in
quadrature, and profit from both channels.

###  IV.2 Training and Validation of the parametrized Neural Networks

Our goal is to find the regions of the two-dimensional parameter space $m_{a}$
versus $c_{tq}/\Lambda$ where the HL-LHC can exclude, at the 95% CL, or
discover the flavons. We trained neural networks to separate signals and
backgrounds in order to maximize the signal significance for a given mass
hypothesis. The coupling just rescales the remaining number of events after
the neural net separation, but the mass impacts the distribution of the
kinematic variables of the signals. In principle, a neural net trained for a
given flavon mass might perform poorly in identifying signals of other masses,
what would force us to train a large number of algorithms to account for the
change in the flavon mass parameter. Instead, it was demonstrated in Baldi
_et al._ (2016); Anzalone  _et al._ (2022) that a neural network parametrized
by a small fixed grid of mass parameters can learn how to interpolate between
the masses of the grid used to train the algorithm. In our case, we trained
our pNN with signal events of the mass grid: 200, 500, 800, 1100, and 1400
GeV. We then used the model to classify events of other eight different
masses. We checked that the pNN learns to distinguish the intermediate masses
correctly by training dedicated neural nets for a couple of intermediate
masses. The conditioning of the model to the flavon mass is done by adding a
column to the dataset for the mass of the signal event. In the case of
backgrounds, we randomly choose a mass value from a uniform distribution in
the mass range of the grid.

Figure 3: Examples of kinematic distributions used as high-level variables in
the pNN training. For smaller $m_{a}$ all distributions are similar to the
$t\bar{t}$ and $tW$ backgrounds. They become more distinct as $m_{a}$ grows,
which partly explains why the pNN identifies signals with higher $m_{a}$ with
greater accuracy.

The variables used to represent the events might facilitate the performance of
the neural networks to separate the signal and background classes. Many high-
level features can be constructed with the four-momenta of the leptons,
bottom, and light-quark jets. We tested many configurations and found a very
good representation of our events with the kinematic variables listed below:

  * •

the transverse momentum of the two leading $b$-jets, the two leptons, and the
two leading non $b$-jets:
$p_{T_{\ell_{1,2}}},p_{T_{b_{1,2}}},p_{T_{j_{1,2}}}$,

  * •

the missing transverse momentum, $\not\\!\\!p_{T}$,

  * •

the energy, at the laboratory center-of-mass frame, of the two leading
$b$-jets, and the two leptons: $E_{\ell_{1,2}}$, $E_{b_{1,2}}$,

  * •

the pseudorapidity of the two leading $b$-jets, the two leptons, and the
leading non $b$-jet: $\eta_{\ell_{1,2}}$, $\eta_{b_{1,2}}$, $\eta_{j_{1}}$,

  * •

the invariant masses: $M_{\ell_{1}\ell_{2}}$, $M_{b_{1}b_{2}}$,
$M_{\ell_{1}b_{1}}$, $M_{\ell_{1}b_{2}}$, $M_{\ell_{2}b_{1}}$,
$M_{\ell_{2}b_{2}}$, $M_{\ell_{1}j_{1}}$, and
${\not\\!\\!M}_{\ell_{1}b_{1}j_{1}}$, the mass of the system
$p_{\ell_{1}}+p_{b_{1}}+p_{j_{1}}+\not\\!\\!p$, where
$\not\\!\\!p=(\sqrt{\not\\!\\!p_{x}^{2}+\not\\!\\!p_{y}^{2}},\not\\!\\!p_{x},\not\\!\\!p_{y},0)$
in the lab frame,

  * •

the difference in the azimuthal angle of the vectors
$p_{\ell_{1}}+p_{\ell_{2}}$ and $p_{b_{1}}+p_{b_{2}}$:
$\Delta\phi_{\ell_{1}\ell_{2}b_{1}b_{2}}=\Delta\phi_{\ell_{1}+\ell_{2},b_{1}+b_{2}}$,

  * •

the invariant angular distances, $\Delta
R_{i,j}=\sqrt{(\eta_{i}-\eta_{j})^{2}+(\phi_{i}-\phi_{j})^{2}}$: $\Delta
R_{\ell_{1}\ell_{2}}$, $\Delta R_{\ell_{1}b_{1}}$, $\Delta R_{\ell_{1}b_{2}}$,
$\Delta R_{\ell_{1}j_{1}}$, $\Delta R_{\ell_{2}b_{1}}$, $\Delta
R_{\ell_{2}b_{2}}$, $\Delta R_{b_{1}b_{2}}$,

  * •

the number of identified $\mu^{+}$, $\mu^{-}$, $e^{+}$, $e^{-}$, $b$-jets and
non-$b$-jets in the final state.

Some of these variable distributions are shown in Fig. 3.

Concerning the algorithm training, we separated the data into three parts:
training, validation, and test sets. After every epoch – the passing of the
whole training data through the algorithm for updating its weights – we
checked the discerning power of the neural network by examining the training
and validation metrics for overfitting.

A Bayesian-inspired optimization of neural network hyperparameters was done
with the Hyperopt package Bergstra  _et al._ (2013). The hyperparameters
considered were the initial learning rate, batch size, number of hidden
layers, activation functions, dropout rate, and the number of neurons at each
hidden layer. The best set we found is shown in Tab. 2. The minimization of
the binary cross-entropy function was carried out with Adam Kingma and Ba
(2017). To avoid overfitting, an early stopping criterion of 25 epochs without
improvement in the loss function on the validation set was adopted with a
total number of epochs set to 200.

hidden layers | neurons | activation | drop-out | learning rate | batch size
---|---|---|---|---|---
3 | 80 | ELU | 0.265 | 0.00156 | 80
Table 2: Hyperparameters and neural net architecture used to train the pNN.
The hyperparameters were optimized with Hyperopt Bergstra  _et al._ (2013).

To guarantee a reliable classification, a balance in the signal-to-background
ratio in the samples of Monte Carlo events must be kept when training the pNN,
according to Ref. Anzalone  _et al._ (2022). The training was performed with
270000 events evenly distributed among signals from the grid of masses and the
three background sources.

The performance of the pNN model was tested using a 10-fold cross-validation,
splitting the training set containing backgrounds and the signals events
corresponding to the grid masses into two datasets: one to train the pNN and
another set held out to compute the model’s accuracy. The mean accuracy across
the 10 validation sets was 89%.

Figure 4: Output scores of signal(red) and background(blue) on the pNN
classification task. The higher flavon masses are more easily distinguished
once their kinematic distributions get more distinct from the backgrounds. We
also show the performance of the algorithm for each background separately. We
also show the corresponding accuracy (acc), defined as the fraction of correct
identifications in the sample.

After training and cross-validation, we used the best cross-validation model
in the test data, which contains signal events with masses from 200 to 1400
GeV in steps of 100 GeV, 9000 events each, and 9000 background events in
total. The test set contains signal events with masses absent in the training
set. Nonetheless, the performance of the pNN on these events is almost as good
as on the grid masses used for training.

The neural network predicts an output score corresponding to the chance of an
event being tagged as a signal. Fig. 4 shows the distribution of the signal
(red) and background (blue) events as functions of the classifier output
score. Lighter flavons are more difficult to discriminate from the background,
as could have been anticipated from the distributions in Fig. 3. In Fig. 5, we
depict the algorithm performance for each test mass in terms of the signal
efficiency, $\varepsilon_{S}$, and background rejection, $1/\varepsilon_{B}$,
for events that survive a given threshold score cut. For a fixed signal
efficiency, the background rejection increases from 200 GeV flavons to 1000
GeV, reaching a maximum and saturating, showing no improvement for larger
masses. The neural network presents an exquisite discerning power, keeping
only $\sim(10^{-3})10^{-4}$ of background events for a signal efficiency of
50% in the (worst)best mass cases. The accuracy (acc) also increases from
light to heavy flavon masses as can be seen from Fig. (4).

Figure 5: The background rejection, $1/\varepsilon_{B}$, against the signal
efficiency, $\varepsilon_{S}$, after cutting on the pNN output scores.

##  V RESULTS

Because the background contamination in SS samples is expected to be rather
low, we first performed a simple cut analysis requiring two same-sign leptons
in the event alongside the cuts of Eq. (IV.1). In this case, we avoid the
large backgrounds to opposite-sign leptons at the cost of halving the number
of signals. The strategy pays off from the perspective of a cut-based
analysis, and the results for 95% CL and $5\sigma$ discovery regions in the
coupling versus mass plane are shown in Fig. 6. The statistical significances
were calculated using the Asimov formula, assuming systematic uncertainties in
the background rates Cowan  _et al._ (2011) assuming an integrated luminosity
of 3 ab-1. Notice the stability of the SS cut-based results against systematic
uncertainties due to the large $S/B$ ratio after cuts.

The neural network allows us to efficiently distinguish between signals and
background irrespective of their lepton charges. The charges of leptons can be
used later to separate the events into OS and SS classes to calculate their
statistical significances separately and then combine them in quadrature,
raising the overall significance compared to the SS significance alone.

After training the pNN, as described in the previous section, we determined
the number of surviving signals and background events by placing cuts on the
test set scores. The cut on the pNN score was scanned to achieve the highest
SS and OS statistical significances, provided they satisfy the minimum number
of events compatible (2 events for 95% CL exclusion and 25 events for
5$\sigma$ discovery). The final significance of the signal is the quadrature
of both SS and OS significances, as discussed.

Figure 6: The red(purple) lines represent the 95%($5\sigma$) CL limits on the
coupling and flavon mass parameters for a luminosity of 3 ab-1. The
dashed(solid) lines depict limits based on a cut-and-count on SS signals (pNN
on SS+OS signals) analysis. The shaded upper right region delimits the region
where EFT is valid. At the left(right) panel, the systematic uncertainty in
the background normalization is 5%(50%).

We plot the resulting 95% CL exclusion and $5\sigma$ discovery regions for pNN
classification in Fig. 6 also assuming an integrated luminosity of 3 ab-1. The
pNN approach probes higher cutoff scales than the SS cut-based approach for
each mass utilizing both SS and OS events. Roughly speaking, pNNs can probe
couplings 50% smaller than the cut-based SS strategy. The advantage faints as
the flavon mass increases, but it is still better up to 1.6 TeV mass, where
the EFT validity is expected to be broken. The region $m_{a}>\Lambda$,
highlighted in magenta in Fig. 6, refers to the validity limit of the EFT.
Effective couplings of order $10^{-2}$ TeV-1 can be probed for flavon masses
up to $\sim 600$ GeV.

Another important feature of our results is the reduced sensitivity to
systematic uncertainties on the background rates. As we see in Fig. 6, the
limits barely change from the 5% to 50% systematics scenario due to the
exquisite job of the neural network in separating signals from the background
with high signal-to-background ratios. These results extend those for
$m_{a}\lesssim 200$ GeV where the flavon prefers to decay into tau leptons
Alves  _et al._ (2024). In combination, they cover the heavy flavon mass range
from 10 to 1600 GeV, complementing the results for lighter flavon masses,
$m_{a}\lesssim 10$ GeV Cornella  _et al._ (2020).

##  VI CONCLUSIONS

Flavor-violating phenomena at colliders are an interesting area for new
physics searches. Contrary to light quarks and leptons, flavor-violating
interactions involving top quarks and up or charm quarks are much less
constrained by searches for flavor-changing neutral currents. Models like the
Froggatt-Nielsen mechanism, which aims to explain the Yukawa hierarchy,
predict the existence of an axion-like particle that can mediate such flavor-
violating interactions, which we call flavons. We assume these scalars
predominantly change the flavor of the fermions they interact with, avoiding
constraints from flavor-conserving searches.

This work explored a distinctive phenomenology of flavons at the 14 TeV LHC.
Pair and single flavon production, plus virtual $t$-channel contributions,
give rise to pairs of opposite and same-charge top quarks. The decay of same-
charge tops to same-sign leptons leads to almost background-free signals at
the expense of halving the number of signal events. A simple cut-and-count
strategy is able to probe flavon couplings of the order of $10^{-2}$ TeV-1.
Instead of wasting all the opposite-charge top quarks, we trained a neural
network to clean up the signals from $t\bar{t}$ and other backgrounds. We
chose to work with parameterized neural networks to get a smooth classifier as
a function of the flavon mass. After the classification, we separated the
samples according to the lepton pairs charges and combined the SS and OS
signal significances in quadrature, thus profiting from all flavons produced.
The combined analyses extended the limits by 50% compared to the cut analysis
based on SS tops only. Moreover, the combined analysis turned out to be very
insensitive to systematic uncertainties in background rates.

With 3 ab-1, we estimate that the 14 TeV LHC can exclude 200(1600) GeV flavons
for couplings as weak as $c_{tq}/\Lambda\sim 3(50)\times 10^{-2}$ TeV-1. The
discovery prospects are also good, 200(1600) GeV flavons can be discovered for
couplings $c_{tq}/\Lambda\gtrsim 6(70)\times 10^{-2}$ TeV-1. Our findings
extend previous searches for light flavons up to 10 GeV Cornella  _et al._
(2020), and heavy flavons with masses up to $\sim 200$ GeV where they might
prefer to decay into tau leptons Alves  _et al._ (2024).

Acknowledgments: This study was financed in part by Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq), via the Grants No.
307317/2021-8 (A. A.) and No. 308280/2023-7 (A. G. D.), and in part by the
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) –
Finance Code 001 (D. S. V. G.). A. A. also acknowledges support from the
FAPESP (No. 2021/01089-1) Grant. E. d. S. A. partially thanks FAPESP for its
support in the initial steps of this work (Grant No. 2022/07185-5).

## References

  * Froggatt and Nielsen (1979) C. D. Froggatt and H. B. Nielsen, Nucl. Phys. B 147, 277 (1979).
  * Di Luzio  _et al._ (2020) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Jaeckel and Ringwald (2010) J. Jaeckel and A. Ringwald, Ann. Rev. Nucl. Part. Sci. 60, 405 (2010), arXiv:1002.0329 [hep-ph] .
  * Ema  _et al._ (2017) Y. Ema, K. Hamaguchi, T. Moroi, and K. Nakayama, JHEP 01, 096 (2017), arXiv:1612.05492 [hep-ph] .
  * Björkeroth  _et al._ (2018) F. Björkeroth, E. J. Chun, and S. F. King, JHEP 08, 117 (2018), arXiv:1806.00660 [hep-ph] .
  * Bauer  _et al._ (2022) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 09, 056 (2022), arXiv:2110.10698 [hep-ph] .
  * Jho  _et al._ (2022) Y. Jho, S. Knapen, and D. Redigolo, JHEP 10, 029 (2022), arXiv:2203.11222 [hep-ph] .
  * Dias  _et al._ (2021) A. G. Dias, J. Leite, and D. S. V. Gonçalves, Phys. Rev. D 104, 075014 (2021), arXiv:2106.07518 [hep-ph] .
  * Di Luzio  _et al._ (2017) L. Di Luzio, F. Mescia, and E. Nardi, Phys.Rev. D96, 075003 (2017), arXiv:1705.05370 [hep-ph] .
  * Bauer  _et al._ (2017) M. Bauer, M. Neubert, and A. Thamm, JHEP 12, 044 (2017), arXiv:1708.00443 [hep-ph] .
  * Aad  _et al._ (2023) G. Aad _et al._ (ATLAS), JHEP 12, 081 (2023), arXiv:2307.14759 [hep-ex] .
  * Hayrapetyan  _et al._ (2024) A. Hayrapetyan _et al._ (CMS), Phys. Lett. B 850, 138478 (2024), arXiv:2311.03261 [hep-ex] .
  * Mehta  _et al._ (2019) P. Mehta, M. Bukov, C.-H. Wang, A. G. R. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, Phys. Rept. 810, 1 (2019), arXiv:1803.08823 [physics.comp-ph] .
  * Carleo  _et al._ (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019), arXiv:1903.10563 [physics.comp-ph] .
  * Bourilkov (2020) D. Bourilkov, Int. J. Mod. Phys. A 34, 1930019 (2020), arXiv:1912.08245 [physics.data-an] .
  * Radovic  _et al._ (2018) A. Radovic, M. Williams, D. Rousseau, M. Kagan, D. Bonacorsi, A. Himmel, A. Aurisano, K. Terao, and T. Wongjirad, Nature 560, 41 (2018).
  * Baldi  _et al._ (2016) P. Baldi, K. Cranmer, T. Faucett, P. Sadowski, and D. Whiteson, The European Physical Journal C 76 (2016), 10.1140/epjc/s10052-016-4099-4.
  * Anzalone  _et al._ (2022) L. Anzalone, T. Diotalevi, and D. Bonacorsi, Machine Learning: Science and Technology 3, 035017 (2022).
  * Alves  _et al._ (2024) A. Alves, A. G. Dias, E. da Silva Almeida, and D. S. V. Gonçalves, Phys. Rev. D 109, 055013 (2024), arXiv:2306.14957 [hep-ph] .
  * Cornella  _et al._ (2020) C. Cornella, P. Paradisi, and O. Sumensari, JHEP 01, 158 (2020), arXiv:1911.06279 [hep-ph] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Frederix  _et al._ (2018) R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H. S. Shao, and M. Zaro, JHEP 07, 185 (2018), [Erratum: JHEP 11, 085 (2021)], arXiv:1804.10017 [hep-ph] .
  * Degrande  _et al._ (2012) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, and T. Reiter, Comput. Phys. Commun. 183, 1201 (2012), arXiv:1108.2040 [hep-ph] .
  * Alloul  _et al._ (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), arXiv:1310.1921 [hep-ph] .
  * Bierlich  _et al._ (2022) C. Bierlich _et al._ ,  (2022), 10.21468/SciPostPhysCodeb.8, arXiv:2203.11601 [hep-ph] .
  * de Favereau  _et al._ (2014) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi (DELPHES 3), JHEP 02, 057 (2014), arXiv:1307.6346 [hep-ex] .
  * Bergstra  _et al._ (2013) J. Bergstra, D. Yamins, and D. D. Cox,  28 (2013).
  * Kingma and Ba (2017) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” (2017), arXiv:1412.6980 [cs.LG] .
  * Cowan  _et al._ (2011) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, The European Physical Journal C 71 (2011), 10.1140/epjc/s10052-011-1554-0.