# Local condensation of charge-$4e$ superconductivity at a nematic domain wall

Matthias Hecker  School of Physics and Astronomy, University of Minnesota,
Minneapolis 55455 MN, USA  Rafael M. Fernandes  School of Physics and
Astronomy, University of Minnesota, Minneapolis 55455 MN, USA  School of
Physics and Astronomy, University of Minnesota, Minneapolis 55455 MN, USA
School of Physics and Astronomy, University of Minnesota, Minneapolis 55455
MN, USA

(August 31, 2024)

###### Abstract

In the fluctuation regime that precedes the onset of pairing in multi-
component superconductors, such as nematic and chiral superconductors, the
normal state is generally unstable towards the formation of charge-$4e$ order
– an exotic quantum state in which electrons form coherent quartets rather
than Cooper pairs. However, charge-$4e$ order is often suppressed by other
competing composite orders, such as nematics. Importantly, the formation of
nematic domains is unavoidable due to the long-range strains generated,
leading to one-dimensional regions where the competing nematic order is
suppressed. Here, we employ a real-space variational approach to demonstrate
that, in such nematic domain walls, charge-$4e$ order is locally condensed via
a vestigial-order mechanism. We explore the experimental manifestations of
this effect and discuss materials in which it can be potentially observed.

### Introduction.

Shortly after the development of the BCS model for superconductivity, it was
recognized that a gas of bosons could also form a coherent state of pairs of
bosons before the Bose-Einstein condensation of individual particles, provided
that strong enough attractive interactions are present [1, 2]. In condensed
matter systems, pair condensation of bosonic quasiparticles has been studied
in various settings, from biexcitons in semiconductors [3, 4] to two-magnon
bound states in frustrated magnets [5, 6]. A fascinating possibility is the
emergence, in superconducting materials, of a coherent state of pairs of
Cooper pairs, dubbed quartets [7, 8, 9, 10, 11, 12]. Theoretically, a
charge-$4e$ superconducting state is expected to display gapless excitations
[13, 14] and half flux-quantum vortices [10]. Experimentally, charge-$4e$
bound states have been recently invoked to explain puzzling magneto-transport
data in kagome superconductors [15, 16].

In the case of bosonic particles, the state with paired bosons is
thermodynamically stable only when there is more than one bosonic “flavor”
available for condensation (e.g. spin-1 bosons) [17]. This suggests that
“multi-flavor” superconductors are a promising setting to search for
charge-$4e$ superconductivity. Indeed, theoretical proposals for quartet
formation have included spin-3/2 systems [18], spinor condensates [19], multi-
band superconductors [20], pair-density waves [10, 21, 22, 23, 24, 25, 26,
27], and multi-component superconductors [28, 29, 30, 31, 32, 33]. In the
latter case, the superconductor is described by multiple gap functions related
by lattice symmetries,
$\boldsymbol{\Delta}=\left(\Delta_{1},\,\Delta_{2},\cdots\right)$; in group-
theory jargon, $\boldsymbol{\Delta}$ transforms as a multi-dimensional
irreducible representation (IR) of the point group. A broad range of pairing
states belong to this category, including several versions of $p$-wave and
$d$-wave states in tetragonal, hexagonal, and cubic lattices [34, 35]. More
importantly, there is experimental evidence for the realization of multi-
component pairing in various systems of interest, from heavy-fermions [36, 37,
38] to moiré superlattices [39] to doped topological insulators [40, 41, 42].

Charge-$4e$ order emerges in multi-component superconductors via the
condensation of a _complex-valued_ composite order parameter
$\left\langle\boldsymbol{\Delta}\cdot\boldsymbol{\Delta}\right\rangle\neq 0$
(distinct from the _real-valued_ composite
$\left\langle\boldsymbol{\Delta}^{\dagger}\cdot\boldsymbol{\Delta}\right\rangle$),
while the superconducting order parameter itself remains zero,
$\left\langle\boldsymbol{\Delta}\right\rangle=0$, such that the transition
temperature of the composite order, $T_{4e}$, is larger than the
superconducting one, $T_{c}$ [28, 43]. This spontaneous symmetry-breaking,
which is driven by fluctuations and thus not captured by mean-field
approaches, lowers the $U(1)$ gauge symmetry to $Z_{2}$; the latter is further
broken if $\left\langle\boldsymbol{\Delta}\right\rangle$ becomes non-zero. It
is said then that the charge-$4e$ and charge-$2e$ superconducting states are
intertwined, and that the former is a vestigial order of the latter [44, 45].

The main obstacle for the stabilization of charge-$4e$ vestigial order is the
competition with other vestigial phases, most notably nematic and
ferromagnetic. Indeed, besides $U(1)$ symmetry, the ground state of a multi-
component superconductor also breaks either time-reversal (chiral
superconductor) or rotational symmetry (nematic superconductor) [46]. These
additional symmetries can be broken before the onset of superconductivity via
the condensation of _real-valued_ composite order parameters. Large-$N$ and
variational calculations found that the corresponding vestigial nematic and
ferromagnetic orders generally preempt the onset of charge-$4e$ order [47, 48]
– except for the special case of a hexagonal nematic superconductor [28].

[][]

Figure 1: (a) Real-space sketch of our main result, the emergence of
charge-$4e$ order at a nematic domain wall. (b) Phase diagram of the vestigial
orders supported by the nematic superconducting ground state of a tetragonal
system. Here, $u$, $v$ and $w$ are the Landau parameters of Eq. (4). The
different shaded regions indicate which vestigial channel is attractive: none
(green), nematic (blue), and leading nematic with subleading charge-$4e$
(red).

In this work, we investigate whether the local suppression of the leading
nematic or ferromagnetic vestigial order enables the local condensation of
charge-$4e$ order. As a motivation, we first discuss a 2D triplet
superconductor without spin-orbit coupling (SOC), whose multi-component gap
function has a continuous SU(2) symmetry. In this case, the continuous
symmetry globally forbids quasi-long-range order of any composite order
parameter, except for charge-$4e$. We then study the more realistic situation
of discrete multi-component superconductors in tetragonal and hexagonal
systems. While long-range order in the competing vestigial channel is
unavoidable, it is locally suppressed due to the formation of domains.
Focusing on a nematic superconductor on the tetragonal lattice, we employ a
real-space variational approach that treats all composite orders on an equal
footing. We find a wide range of parameters for which charge-$4e$ order is
condensed at the nematic domain walls while the superconducting order
parameter remains uncondensed, as illustrated schematically in Fig. 1(a). We
finish by discussing the experimental manifestations of this local
condensation and candidate materials.

### Global suppression of the competing vestigial phases.

To set the stage, we discuss a special case in which, for symmetry reasons
alone, the only order that can be stabilized is charge-$4e$. Consider a 2D
system in which the electrons experience negligible SOC (e.g. graphene) and,
in addition to spin, have another pseudospin degree of freedom (e.g.
sublattice or valley). Enforcing the pairing state to be momentum-independent
and pseudospin-singlet, the gap function must be spin-triplet and represented
by an order parameter
$\boldsymbol{\Delta}=(\Delta_{1},\,\Delta_{2},\,\Delta_{3})^{T}$ that
transforms as a vector in SU(2) spin space. This is nothing but the
$\boldsymbol{d}$-vector of a triplet superconductor [34, 35], albeit even in
momentum. While there is no indication of valley-singlet spin-triplet
superconductivity in graphene, this type of state has been studied in twisted
moiré systems [49, 50, 51, 52] and Bernal bilayer graphene [32, 53]. The
superconducting action is given by

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{\mathsf{q}}\chi_{\mathsf{q}}^{-1}\left|\boldsymbol{\Delta}_{\mathsf{q}}\right|^{2}+\int_{\mathsf{x}}\Big{(}u\left|\boldsymbol{\Delta}_{\mathsf{x}}\right|^{4}+v\left|\boldsymbol{\Delta}_{\mathsf{x}}\cdot\boldsymbol{\Delta}_{\mathsf{x}}\right|^{2}\Big{)},$
(1)

with (bare) superconducting susceptibility $\chi_{\mathsf{q}}^{-1}$ and
$\mathsf{q}=\left(\mathbf{q},\omega_{n}\right)$ denoting momentum and
Matsubara frequency. The interaction part has Landau coefficients $u$ and $v$,
where $\mathsf{x}=\left(\mathbf{x},\tau\right)$ denotes position and imaginary
time. The mean-field phase diagram of this model is well-established [54, 55],
displaying different types of unitary and non-unitary pairing depending on the
sign of $v$.

To describe the vestigial orders, however, it is necessary to go beyond mean-
field [45]; for our purposes, a group-theoretical analysis is sufficient.
There are four symmetry-breaking bilinear combinations of
$\boldsymbol{\Delta}$: two real-valued composites
$\boldsymbol{\Psi}^{(l=1)}=\mathsf{i}\boldsymbol{\Delta}\times\bar{\boldsymbol{\Delta}}$
and
$\Psi_{\mu\nu}^{(l=2)}=\frac{1}{2}\left(\Delta_{\mu}\bar{\Delta}_{\nu}+\Delta_{\nu}\bar{\Delta}_{\mu}\right)-\frac{1}{3}\,\delta_{\mu\nu}\left|\boldsymbol{\Delta}\right|^{2}$
and two complex-valued ones,
$\psi^{(l=0)}=\boldsymbol{\Delta}\cdot\boldsymbol{\Delta}$ and
$\psi_{\mu\nu}^{(l=2)}=\Delta_{\mu}\Delta_{\nu}-\frac{1}{3}\,\delta_{\mu\nu}\left(\boldsymbol{\Delta}\cdot\boldsymbol{\Delta}\right)$,
where $\bar{z}\equiv z^{*}$. The superscript indicates the transformation
properties in SU(2) spin-space, corresponding to a scalar ($l=0$), a vector
($l=1$), or a tensor ($l=2$). Each composite has a clear physical
interpretation: $\boldsymbol{\Psi}^{(l=1)}$ corresponds to time-reversal
symmetry-breaking and, thus, vestigial ferromagnetic order.
$\Psi_{\mu\nu}^{(l=2)}$ is associated with rotational symmetry-breaking in
spin-space, and therefore denotes (spin-)nematic vestigial order [56, 57, 31].
Finally, $\psi^{(l=0)}$ and $\psi_{\mu\nu}^{(l=2)}$ correspond to $s$-wave and
$d$-wave charge-$4e$ vestigial orders, respectively.

Because $\boldsymbol{\Psi}^{(l=1)}$, $\Psi_{\mu\nu}^{(l=2)}$,
$\psi_{\mu\nu}^{(l=2)}$, and
$\boldsymbol{\Delta}\equiv\boldsymbol{\Delta}^{(l=1)}$ itself transform non-
trivially in SU(2) spin-space (i.e. they are at least Heisenberg-type order
parameters), none of them can sustain (quasi-)long-range order at non-zero
temperatures in 2D. On the other hand, because $\psi^{(l=0)}$ is a complex
scalar (i.e. XY-type), it can establish quasi-long-range order through a BKT
transition. Therefore, the only state allowed to develop quasi-long-range
order in this model is the vestigial charge-$4e$ phase. We note that similar
2D and 1D models for triplet superconductors [31, 32, 33] and spinor
condensates [19] have been previously studied and shown to support charge-$4e$
order.

### Local suppression of the competing vestigial phases.

Despite illuminating, the simple model above is not representative of
realistic multi-component superconductors, where either SOC is not negligible
or singlet states are realized. Yet, it highlights an efficient strategy to
realize charge-$4e$ order: suppression of the other, leading, vestigial
phases. While this generally cannot be accomplished globally via Mermin-
Wagner’s theorem, vestigial nematic or ferromagnetic states tend to form
domains to minimize the elastic or magnetic dipolar energies. To explore this
idea, we consider a generic two-component superconducting order parameter
$\boldsymbol{\Delta}=(\Delta_{1},\Delta_{2})$, which could describe
$(p_{x},p_{y})$-wave or $(d_{xz},d_{yz})$-wave states in tetragonal and
hexagonal lattices [34, 35], or $(d_{x^{2}-y^{2}},d_{xy})$-wave pairing in
hexagonal systems and $45^{\circ}$-twisted bilayer tetragonal $d$-wave
superconductors [58]. In contrast to the previous example,
$\boldsymbol{\Delta}$ now transforms as a two-dimensional IR of a discrete
point group. While we will focus on tetragonal ($\mathrm{D_{4h}}$)
superconductors, where $\boldsymbol{\Delta}$ transforms as the IR $E_{g}$ or
$E_{u}$, the conclusions apply to all other cases.

We start by classifying all possible composite order parameters. As shown in
[43], seven different bilinear combinations can be formed. Apart from the
symmetry-preserving bilinear
$\Psi^{A_{1g}}=\boldsymbol{\Delta}^{\dagger}\tau^{0}\boldsymbol{\Delta}$, with
Pauli matrices $\tau^{i}$ acting on the $\boldsymbol{\Delta}$ subspace, there
are three additional real-valued bilinears:

$\displaystyle\Psi^{A_{2g}}$
$\displaystyle=\boldsymbol{\Delta}^{\dagger}\tau^{y}\boldsymbol{\Delta},$
$\displaystyle\Psi^{B_{1g}}$
$\displaystyle=\boldsymbol{\Delta}^{\dagger}\tau^{z}\boldsymbol{\Delta},$
$\displaystyle\Psi^{B_{2g}}$
$\displaystyle=\boldsymbol{\Delta}^{\dagger}\tau^{x}\boldsymbol{\Delta}.$ (2)

Here, $\Psi^{A_{2g}}$ breaks time-reversal-symmetry and causes ferromagnetism,
while $\Psi^{B_{1g}}$ and $\Psi^{B_{2g}}$ break tetragonal symmetry and cause
nematicity. The three complex-valued bilinears are given by

$\displaystyle\psi^{A_{1g}}$
$\displaystyle=\boldsymbol{\Delta}^{T}\tau^{0}\boldsymbol{\Delta},$
$\displaystyle\psi^{B_{1g}}$
$\displaystyle=\boldsymbol{\Delta}^{T}\tau^{z}\boldsymbol{\Delta},$
$\displaystyle\psi^{B_{2g}}$
$\displaystyle=\boldsymbol{\Delta}^{T}\tau^{x}\boldsymbol{\Delta},$ (3)

and describe, respectively, $s$-wave, $d_{x^{2}-y^{2}}$-wave, and
$d_{xy}$-wave charge-$4e$ superconductivity. In our notation, $\Psi^{n}$
($\psi^{n}$) denotes real-valued (complex-valued) bilinears, whereas the
superscript $n$ indicates the IR according to which the composite transforms.
The superconducting action is
$\mathcal{S}=\int_{\mathsf{x}}r_{0}\left|\boldsymbol{\Delta}_{\mathsf{x}}\right|^{2}+\;\mathcal{S}^{\mathrm{grad}}\;+\;\mathcal{S}^{\mathrm{int}}\,,$
where $r_{0}=a_{0}(T-T_{0})$ denotes the bare SC transition ($a_{0},T_{0}>0$)
and $\mathcal{S}^{\mathrm{grad}}$ contains the symmetry-allowed gradient terms
[35]. The interaction part is given by [43]

$\displaystyle\mathcal{S}^{\mathrm{int}}$
$\displaystyle=\int_{\mathsf{x}}\Big{[}u\,(\Psi_{\mathsf{x}}^{A_{1g}})^{2}+v\,(\Psi_{\mathsf{x}}^{A_{2g}})^{2}+w\,(\Psi_{\mathsf{x}}^{B_{1g}})^{2}\Big{]},$
(4)

and contains three independent interaction parameters $u>0$ and $v,w>-u$. The
mean-field phase diagram in the $\left(\frac{v}{u},\,\frac{w}{u}\right)$
parameter space is well-established, displaying chiral and two types of
nematic superconductivity [35, 10, 46].

The vestigial orders associated with each mean-field ground state were
analyzed in Ref. [43] via a variational approach. The leading vestigial
instability is always that of a real-valued composite (nematic or
ferromagnetic) whereas the vestigial charge-$4e$ orders are always subleading.
While our conclusions hold across the entire phase diagram, hereafter we focus
on the $v>0>w$ region [Fig. 1(b)], where the superconducting ground state is
nematic and the competing vestigial phases are $d_{x^{2}-y^{2}}$ nematic
($\Psi^{B_{1g}}$) and $s$-wave charge-$4e$ ($\psi^{A_{1g}}$). In bulk, the
only vestigial order realized is the nematic one [43]. However, because the
effective interaction in the charge-$4e$ channel is attractive, the system
could gain energy by condensing this mode in regions where nematic order is
suppressed. Due to its Ising-like character, $\Psi^{B_{1g}}$ can sustain long-
range order at non-zero temperatures. However, because of the linear coupling
between $\Psi^{B_{1g}}$ and strain, nematic domains must form to minimize the
elastic energy and accommodate long-range lattice deformations. This opens up
the possibility of $\psi^{A_{1g}}$ condensation at nematic domain walls.

### Charge-4e condensation at the domain wall.

To proceed, we employ a real-space Gaussian variational approach, which treats
all vestigial channels equally, on a 1D grid of length $L$ and $i=1,\dots,N$
sites. The variational ansatz consists of a trial action $\mathcal{S}_{0}$
[59, 60, 61, 43], which in our case is
$\mathcal{S}_{0}=\frac{1}{2}\frac{L}{T}\sum_{i}\hat{\boldsymbol{\Delta}}_{i}^{\dagger}\,G_{i}^{-1}\,\hat{\boldsymbol{\Delta}}_{i}+\mathcal{S}^{\mathrm{grad}}$.
Represented in the Nambu basis
$\hat{\boldsymbol{\Delta}}_{i}=(\boldsymbol{\Delta}_{i},\bar{\boldsymbol{\Delta}}_{i})$,
the local inverse Green’s function

$\displaystyle G_{i}^{-1}$
$\displaystyle=\left(\begin{smallmatrix}R_{i}+\Phi_{i}^{B_{1g}}&\Phi_{i}^{B_{2g}}-\mathsf{i}\Phi_{i}^{A_{2g}}&\phi_{i}^{A_{1g}}+\phi_{i}^{B_{1g}}&\phi_{i}^{B_{2g}}\\\
&R_{i}-\Phi_{i}^{B_{1g}}&\phi_{i}^{B_{2g}}&\phi_{i}^{A_{1g}}-\phi_{i}^{B_{1g}}\\\
&&R_{i}+\Phi_{i}^{B_{1g}}&\Phi_{i}^{B_{2g}}+\mathsf{i}\Phi_{i}^{A_{2g}}\\\
\mathrm{H.c.}&&&R_{i}-\Phi_{i}^{B_{1g}}\end{smallmatrix}\right),$ (5)

contains the real-valued ($\Phi_{i}^{n}$) and complex-valued
$(\phi_{i}^{n},\bar{\phi}_{i}^{n})$ variational parameters, and the mass
renormalization parameter $R_{i}=r_{0}+\Phi_{i}^{A_{1g}}$. Because
$\mathcal{S}_{0}$ is Gaussian, it is straightforward to compute the
variational free energy

$\displaystyle F_{v}$ $\displaystyle=-T\log
Z_{0}+T\langle\mathcal{S}-\mathcal{S}_{0}\rangle_{0},$ (6)

where $Z_{0}=\int
D(\boldsymbol{\Delta},\bar{\boldsymbol{\Delta}})\,e^{-\mathcal{S}_{0}}$ is the
partition function of the trial action. The detailed evaluation of Eq. (6) is
shown in the Supplementary Material (SM). Importantly, the original Landau
coefficients $u$, $v$, and $w$ appear in $F_{v}$ as different combinations in
each symmetry channel, corresponding to effective interactions
$U_{A_{1g}}=3u+v+w$, $U_{A_{2g}}=u+3v-w,$ $U_{B_{1g}}=u-v+3w$,
$U_{B_{2g}}=u-v-w$, $u_{A_{1g}}=u-v+w$, $u_{B_{1g}}=u+v+w$, and
$u_{B_{2g}}=u+v-w$ [43]. An attractive interaction ($U_{n},u_{n}<0$) indicates
a potential instability, signaled by the condensation of the corresponding
variational order parameter ($\Phi_{i}^{n},\,\phi_{i}^{n}$). Importantly, a
vestigial phase only emerges if superconductivity is not immediately triggered
by $\Phi_{i}^{n},\,\phi_{i}^{n}\neq 0$. In the variational approach, this
condition can be verified by confirming that the variational superconducting
susceptibility $\chi$, which is given by a combination of $R_{i}$,
$\Phi_{i}^{n}$, and $\phi_{i}^{n}$ (see SM for details), remains finite.

Figure 2: Variational nematic order parameter ($\Phi^{B_{1g}}$), inverse
superconducting susceptibility ($\chi^{-1}$), and variational charge-$4e$
order parameter ($\phi^{A_{1g}}$) obtained from the numerical minimization of
the variational free energy (6) across two nematic domains (all in units of
the gradient-term stiffness $t_{0}$). Each panel corresponds to a different
temperature, parametrized by $r_{0}\propto T-T_{0}$. The parameters used here
are $w=-0.5u$ and $v=0.1u$ [red cross in Fig. 1(b)]. Here and in Fig. 3, we
set $t_{0}/\sqrt{uT_{0}/L}=20/7$.  Figure 3: Same variational parameters as in
Fig. 2, but evaluated for the parameters $w=-0.05u$ and $v=2.1u$ [purple cross
in Fig. 1(b)], corresponding to a subleading attractive charge-$4e$ channel
($u_{A_{1g}}<0$). Local condensation of $\phi^{A_{1g}}$ is observed at the
nematic domain wall.

In bulk, for most of the $\left(\frac{v}{u},\,\frac{w}{u}\right)$ phase
diagram, there are two competing attractive vestigial-order channels,
corresponding to a real-valued composite (nematic/ferromagnetic) and a
complex-valued one ($s$-wave/$d$-wave charge-$4e$). As demonstrated in Ref.
[43], the nematic/ferromagnetic instability is always the leading one in bulk
[Fig. 1(b)]. Our goal is to determine the fate of these phases along a nematic
domain wall. We therefore numerically minimize the free energy (6) for the
$N$-site 1D grid with domain-wall boundary conditions
$\Phi_{1}^{B_{1g}}=-\Phi_{N}^{B_{1g}}=\Phi_{0}^{B_{1g}}$, where
$\Phi_{0}^{B_{1g}}$ is the (self-consistently obtained) bulk value of the
nematic order parameter (see SM for additional details). Consider first the
phase diagram region where the only attractive vestigial-order channel is the
nematic, i.e. $U_{B_{1g}}<0$ but $u_{A_{1g}}>0$ [red cross in Fig. 1(b)]. The
results are shown in Fig. 2. As the control parameter $r_{0}\propto T-T_{0}$
is decreased, a nematic domain emerges below a temperature that coincides with
the bulk nematic critical temperature. As expected, the domain wall becomes
sharper as the temperature is lowered, since the wall width scales as
$t_{0}/|\Phi_{0}^{B_{1g}}|$ where $t_{0}$ is the gradient-term stiffness. At
the domain-wall center, the superconducting susceptibility $\chi$ is
suppressed, consistent with the fact that vestigial order enhances the
superconducting transition. Interestingly, $\chi^{-1}$ has a non-monotonic
spatial dependence, displaying a dip at the domain-wall boundaries, which can
lead to local condensation of superconductivity [yellow region of Fig. 2(c)].
No sign of charge-$4e$ order is observed, consistent with a repulsive
effective interaction in this channel.

Consider now the phase-diagram region where the charge-$4e$ vestigial channel
is attractive, but subleading to the nematic, $U_{B_{1g}}<u_{A_{1g}}<0$
[purple cross in Fig. 1(b)] . As shown in Fig. 3, upon decreasing $r_{0}$,
nematic order emerges first, in two domains. However, in contrast to Fig. 2,
the charge-$4e$ order parameter $\phi^{A_{1g}}$ condenses inside the domain
wall while the superconducting susceptibility $\chi$ remains finite. This is
the main result of the paper. Because the domain wall is one-dimensional, this
should be understood as a local condensation, since phase slips will destroy
long-range order along the wall. Note, even though the $d$-wave charge-$4e$
channel is repulsive in this phase-diagram region, the order parameter
$\phi^{B_{1g}}$ condenses due to the trilinear coupling between
$\Phi^{B_{1g}}$, $\phi^{A_{1g}}$, and $\bar{\phi}^{B_{1g}}$ (see the inset)
[43]. We verified that this behavior is not particular to this set of
parameters, but occurs in any region of the
$\left(\frac{v}{u},\,\frac{w}{u}\right)$ phase diagram where $s$-wave/$d$-wave
vestigial charge-$4e$ order is sufficiently attractive but subleading to
nematic/ferromagnetic vestigial order.

### Discussion.

In this work, we employed a variational approach to demonstrate that
charge-$4e$ order locally condenses at nematic domain walls that emerge in the
normal-state of nematic superconductors upon approaching the superconducting
transition. Before the onset of superconductivity, the system is unstable
towards the formation of vestigial nematic and charge-$4e$ orders arising from
the fluctuation-induced condensation of composite superconducting order
parameters. Because nematic order generally wins over charge-$4e$ order, the
suppression of nematic order along the domain wall enables the system to
further minimize the free energy by condensing charge-$4e$ order. Since an
analogous mechanism also applies for chiral superconductors, our results point
to a wide class of systems – multi-component superconductors – where local
charge-$4e$ order can potentially emerge. Among the materials for which there
is strong experimental evidence for nematic superconductivity, doped
$\mathrm{Bi_{2}Se_{3}}$ [62, 40, 42, 41] and twisted bilayer graphene (TBG)
[39, 63, 64] are the most promising candidates to display this effect. Indeed,
a vestigial nematic phase exists in doped $\mathrm{Bi_{2}Se_{3}}$ [65, 66],
whereas in TBG, normal-state nematic order appears close to the
superconducting dome [39]. There are also several chiral-superconductor
candidates [67], including the widely studied heavy-fermion UPt3 [36, 37, 38].

An important question is how to experimentally detect this effect. Since
charge-$4e$ order emerges at nematic domain walls, local probes such as
scanning tunneling microscopy (STM) and scanning near-field optical microscopy
(SNOM) are ideal. Because the charge-$4e$ state is expected to be gapless [13,
14], its density-of-states (DOS) profile, which can be accessed in a standard
STM measurement, will likely differ from the normal-state DOS only in subtle
ways. On the other hand, Josephson-STM, in which a superconducting tip is used
[68, 69], could provide more direct evidence for quartets. The SNOM technique
has the unique capability of probing the local optical response with a few
nanometers resolution, from which the properties of the local optical
conductivity $\sigma(\omega,\boldsymbol{r})$ can be inferred [70, 71]. Because
the charge-$4e$ state has a non-zero superfluid density [14],
$\mathrm{Im}\,\sigma(\omega,\boldsymbol{r})\sim 1/\omega$ is expected to
emerge at low frequencies near nematic domain walls. While this behavior could
also be due to a charge-$2e$ superconducting filament, only in the charge-$4e$
case this behavior would be accompanied by a gapless DOS, which can be probed
via STM. These results also reveal the tantalizing possibility of using
uniaxial strain to control the charge-$4e$ phase, since beyond a critical
strain value, the sample becomes mono-nematic-domain and local charge-$4e$
order disappears.

###### Acknowledgements.

We thank J. Schmalian, L. Fu, and R. Willa for fruitful discussions, and
particularly Y. Wang for in-depth discussions about the properties of
charge-$4e$ superconductors. This work was supported by the U. S. Department
of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and
Engineering Division, under Award No. DE-SC0020045.

## References

  * Valatin and Butler [1958] J. Valatin and D. Butler, On the collective properties of a boson system, Nuovo Cimento 10, 37 (1958).
  * Evans and Imry [1969] W. Evans and Y. Imry, On the pairing theory of the bose superfluid, Nuovo Cimento B 63, 155 (1969).
  * Lampert [1958] M. A. Lampert, Mobile and immobile effective-mass-particle complexes in nonmetallic solids, Phys. Rev. Lett. 1, 450 (1958).
  * Moskalenko [1958] S. A. Moskalenko, Opt. Spektrosk 5, 147 (1958).
  * Wang  _et al._ [2017] Z. Wang, A. E. Feiguin, W. Zhu, O. A. Starykh, A. V. Chubukov, and C. D. Batista, Chiral liquid phase of simple quantum magnets, Phys. Rev. B 96, 184409 (2017).
  * Strockoz  _et al._ [2022] J. Strockoz, D. S. Antonenko, D. LaBelle, and J. W. Venderbos, Excitonic instability towards a Potts-nematic quantum paramagnet, arXiv:2211.11739 (2022).
  * Korshunov [1985] S. Korshunov, Two-dimensional superfluid Fermi liquid with $p$-pairing, Zh. Eksp. Teor. Fiz 89, 539 (1985).
  * Volovik [1992] G. E. Volovik,  _Exotic properties of superfluid 3He_, Vol. 1 (World Scientific, 1992).
  * Aligia  _et al._ [2005] A. A. Aligia, A. P. Kampf, and J. Mannhart, Quartet formation at $(100)/(110)$ interfaces of $d$-wave superconductors, Phys. Rev. Lett. 94, 247004 (2005).
  * Berg  _et al._ [2009] E. Berg, E. Fradkin, and S. A. Kivelson, Charge-$4e$ superconductivity from pair-density-wave order in certain high-temperature superconductors, Nature Physics 5, 830 (2009).
  * Moon [2012] E.-G. Moon, Skyrmions with quadratic band touching fermions: A way to achieve charge $4e$ superconductivity, Phys. Rev. B 85, 245123 (2012).
  * Khalaf  _et al._ [2022] E. Khalaf, P. Ledwith, and A. Vishwanath, Symmetry constraints on superconductivity in twisted bilayer graphene: Fractional vortices, $4e$ condensates, or nonunitary pairing, Phys. Rev. B 105, 224508 (2022).
  * Jiang  _et al._ [2017] Y.-F. Jiang, Z.-X. Li, S. A. Kivelson, and H. Yao, Charge-$4e$ superconductors: A Majorana quantum Monte Carlo study, Phys. Rev. B 95, 241103 (2017).
  * Gnezdilov and Wang [2022] N. V. Gnezdilov and Y. Wang, Solvable model for a charge-$4e$ superconductor, Phys. Rev. B 106, 094508 (2022).
  * Ge  _et al._ [2022] J. Ge, P. Wang, Y. Xing, Q. Yin, H. Lei, Z. Wang, and J. Wang, Discovery of charge-$4e$ and charge-$6e$ superconductivity in kagome superconductor $CsV_{3}Sb_{5}$, arXiv:2201.10352 (2022).
  * Varma and Wang [2023] C. M. Varma and Z. Wang, Extended superconducting fluctuation region and 6e and 4e flux-quantization in a kagome compound with a normal state of 3q-order, arXiv:2307.00448 (2023).
  * Noziéres and Saint James [1982] P. Noziéres and D. Saint James, Particle vs. pair condensation in attractive Bose liquids, Journal de Physique 43, 1133 (1982).
  * Wu [2005] C. Wu, Competing orders in one-dimensional spin-$3/2$ fermionic systems, Phys. Rev. Lett. 95, 266404 (2005).
  * Mukerjee  _et al._ [2006] S. Mukerjee, C. Xu, and J. E. Moore, Topological defects and the superfluid transition of the $s=1$ spinor condensate in two dimensions, Phys. Rev. Lett. 97, 120406 (2006).
  * Herland  _et al._ [2010] E. V. Herland, E. Babaev, and A. Sudbø, Phase transitions in a three dimensional $U(1)\times U(1)$ lattice London superconductor: Metallic superfluid and charge-$4e$ superconducting states, Phys. Rev. B 82, 134511 (2010).
  * Radzihovsky and Vishwanath [2009] L. Radzihovsky and A. Vishwanath, Quantum liquid crystals in an imbalanced Fermi gas: Fluctuations and fractional vortices in Larkin-Ovchinnikov states, Phys. Rev. Lett. 103, 010404 (2009).
  * Radzihovsky [2011] L. Radzihovsky, Fluctuations and phase transitions in Larkin-Ovchinnikov liquid-crystal states of a population-imbalanced resonant Fermi gas, Phys. Rev. A 84, 023611 (2011).
  * Agterberg  _et al._ [2011] D. F. Agterberg, M. Geracie, and H. Tsunetsugu, Conventional and charge-six superfluids from melting hexagonal Fulde-Ferrell-Larkin-Ovchinnikov phases in two dimensions, Phys. Rev. B 84, 014513 (2011).
  * Agterberg  _et al._ [2020] D. F. Agterberg, J. C. S. Davis, S. D. Edkins, E. Fradkin, D. J. V. Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, The physics of pair-density waves: Cuprate superconductors and beyond, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * Shaffer  _et al._ [2023] D. Shaffer, F. J. Burnell, and R. M. Fernandes, Weak-coupling theory of pair density wave instabilities in transition metal dichalcogenides, Phys. Rev. B 107, 224516 (2023).
  * Yu [2023] Y. Yu, Nondegenerate surface pair density wave in the kagome superconductor ${\mathrm{csv}}_{3}{\mathrm{sb}}_{5}$: Application to vestigial orders, Phys. Rev. B 108, 054517 (2023).
  * Wu and Wang [2023] Y.-M. Wu and Y. Wang, $d$-wave charge-$4e$ superconductivity from fluctuating pair density waves, arXiv:2303.17631 (2023).
  * Fernandes and Fu [2021] R. M. Fernandes and L. Fu, Charge-$4e$ superconductivity from multicomponent nematic pairing: Application to twisted bilayer graphene, Phys. Rev. Lett. 127, 047001 (2021).
  * Jian  _et al._ [2021] S.-K. Jian, Y. Huang, and H. Yao, Charge-$4e$ superconductivity from nematic superconductors in two and three dimensions., Phys. Rev. Lett. 127, 227001 (2021).
  * Zeng  _et al._ [2021] M. Zeng, L.-H. Hu, H.-Y. Hu, Y.-Z. You, and C. Wu, Phase-fluctuation induced time-reversal symmetry breaking normal state, arXiv:2102.06158 (2021).
  * Rampp  _et al._ [2022] M. A. Rampp, E. J. König, and J. Schmalian, Topologically enabled superconductivity, Phys. Rev. Lett. 129, 077001 (2022).
  * Curtis  _et al._ [2023] J. B. Curtis, N. R. Poniatowski, Y. Xie, A. Yacoby, E. Demler, and P. Narang, Stabilizing fluctuating spin-triplet superconductivity in graphene via induced spin-orbit coupling, Phys. Rev. Lett. 130, 196001 (2023).
  * Poduval and Scheurer [2023] P. P. Poduval and M. S. Scheurer, Vestigial singlet pairing in a fluctuating magnetic triplet superconductor: Applications to graphene moiré systems, arXiv:2301.01344 (2023).
  * Annett [1990] J. F. Annett, Symmetry of the order parameter for high-temperature superconductivity, Advances in Physics 39, 83 (1990).
  * Sigrist and Ueda [1991] M. Sigrist and K. Ueda, Phenomenological theory of unconventional superconductivity, Rev. Mod. Phys. 63, 239 (1991).
  * Sauls [1994] J. Sauls, The order parameter for the superconducting phases of $UPt_{3}$, Advances in Physics 43, 113 (1994).
  * Schemm  _et al._ [2014] E. Schemm, W. Gannon, C. Wishne, W. Halperin, and A. Kapitulnik, Observation of broken time-reversal symmetry in the heavy-fermion superconductor $UPt_{3}$, Science 345, 190 (2014).
  * Avers  _et al._ [2020] K. E. Avers, W. J. Gannon, S. J. Kuhn, W. P. Halperin, J. Sauls, L. DeBeer-Schmitt, C. Dewhurst, J. Gavilano, G. Nagy, U. Gasser, _et al._ , Broken time-reversal symmetry in the topological superconductor $UPt_{3}$, Nature Physics 16, 531 (2020).
  * Cao  _et al._ [2021] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, Science 372, 264 (2021).
  * Matano  _et al._ [2016] K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-q. Zheng, Spin-rotation symmetry breaking in the superconducting state of $Cu_{x}Bi_{2}Se_{3}$, Nature Physics 12, 852 (2016).
  * Pan  _et al._ [2016] Y. Pan, A. Nikitin, G. Araizi, Y. Huang, Y. Matsushita, T. Naka, and A. De Visser, Rotational symmetry breaking in the topological superconductor $Sr_{x}Bi_{2}Se_{3}$ probed by upper-critical field experiments, Scientific reports 6, 28632 (2016).
  * Asaba  _et al._ [2017] T. Asaba, B. J. Lawson, C. Tinsman, L. Chen, P. Corbae, G. Li, Y. Qiu, Y. S. Hor, L. Fu, and L. Li, Rotational symmetry breaking in a trigonal superconductor $Nb$-doped $Bi_{2}Se_{3}$, Phys. Rev. X 7, 011009 (2017).
  * Hecker  _et al._ [2023] M. Hecker, R. Willa, J. Schmalian, and R. M. Fernandes, Cascade of vestigial orders in two-component superconductors: Nematic, ferromagnetic, $s$-wave charge-$4e$, and $d$-wave charge-$4e$ states, Phys. Rev. B 107, 224503 (2023).
  * Fradkin  _et al._ [2015] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Colloquium: Theory of intertwined orders in high temperature superconductors, Rev. Mod. Phys. 87, 457 (2015).
  * Fernandes  _et al._ [2019] R. M. Fernandes, P. P. Orth, and J. Schmalian, Intertwined vestigial order in quantum materials: Nematicity and beyond, Annual Review of Condensed Matter Physics 10, 133 (2019).
  * Gali and Fernandes [2022] V. Gali and R. M. Fernandes, Role of electromagnetic gauge-field fluctuations in the selection between chiral and nematic superconductivity, Phys. Rev. B 106, 094509 (2022).
  * Hecker and Schmalian [2018] M. Hecker and J. Schmalian, Vestigial nematic order and superconductivity in the doped topological insulator $Cu_{x}Bi_{2}Se_{3}$, npj Quantum Materials 3, 1 (2018).
  * Willa [2020] R. Willa, Symmetry-mixed bound-state order, Phys. Rev. B 102, 180503 (2020).
  * Venderbos and Fernandes [2018] J. W. F. Venderbos and R. M. Fernandes, Correlations and electronic order in a two-orbital honeycomb lattice model for twisted bilayer graphene, Phys. Rev. B 98, 245103 (2018).
  * Scheurer and Samajdar [2020] M. S. Scheurer and R. Samajdar, Pairing in graphene-based moiré superlattices, Phys. Rev. Research 2, 033062 (2020).
  * Wang  _et al._ [2021] Y. Wang, J. Kang, and R. M. Fernandes, Topological and nematic superconductivity mediated by ferro-su(4) fluctuations in twisted bilayer graphene, Phys. Rev. B 103, 024506 (2021).
  * Lake  _et al._ [2022] E. Lake, A. S. Patri, and T. Senthil, Pairing symmetry of twisted bilayer graphene: A phenomenological synthesis, Phys. Rev. B 106, 104506 (2022).
  * Dong  _et al._ [2023] Z. Dong, A. V. Chubukov, and L. Levitov, Transformer spin-triplet superconductivity at the onset of isospin order in bilayer graphene, Phys. Rev. B 107, 174512 (2023).
  * Volovik and Gor’kov [1985] G. Volovik and L. Gor’kov, Superconducting classes in heavy-fermion systems, Zh. Eksp. Teor. Fiz 88, 1412 (1985).
  * Blagoeva  _et al._ [1990] E. J. Blagoeva, G. Busiello, L. De Cesare, Y. T. Millev, I. Rabuffo, and D. I. Uzunov, Fluctuation-induced first-order transitions in unconventional superconductors, Phys. Rev. B 42, 6124 (1990).
  * Blume and Hsieh [1969] M. Blume and Y. Hsieh, Biquadratic exchange and quadrupolar ordering, J. Appl. Phys. 40, 1249 (1969).
  * Andreev and Grishchuk [1984] A. Andreev and I. Grishchuk, Spin nematics, Sov. Phys. JETP 60, 267 (1984).
  * Haenel  _et al._ [2022] R. Haenel, T. Tummuru, and M. Franz, Incoherent tunneling and topological superconductivity in twisted cuprate bilayers, Phys. Rev. B 106, 104505 (2022).
  * Moshe and Zinn-Justin [2003] M. Moshe and J. Zinn-Justin, Quantum field theory in the large N limit: a review, Physics Reports 385, 69 (2003).
  * Fischer and Berg [2016] M. H. Fischer and E. Berg, Fluctuation and strain effects in a chiral $p$-wave superconductor, Phys. Rev. B 93, 054501 (2016).
  * Nie  _et al._ [2017] L. Nie, A. V. Maharaj, E. Fradkin, and S. A. Kivelson, Vestigial nematicity from spin and/or charge order in the cuprates, Phys. Rev. B 96, 085142 (2017).
  * Fu [2014] L. Fu, Odd-parity topological superconductor with nematic order: Application to $Cu_{x}Bi_{2}Se_{3}$, Phys. Rev. B 90, 100509 (2014).
  * Kozii  _et al._ [2019] V. Kozii, H. Isobe, J. W. F. Venderbos, and L. Fu, Nematic superconductivity stabilized by density wave fluctuations: Possible application to twisted bilayer graphene, Phys. Rev. B 99, 144507 (2019).
  * Chichinadze  _et al._ [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Nematic superconductivity in twisted bilayer graphene, Phys. Rev. B 101, 224513 (2020).
  * Sun  _et al._ [2019] Y. Sun, S. Kittaka, T. Sakakibara, K. Machida, J. Wang, J. Wen, X. Xing, Z. Shi, and T. Tamegai, Quasiparticle evidence for the nematic state above ${T}_{\mathrm{c}}$ in $Sr_{x}Bi_{2}Se_{3}$, Phys. Rev. Lett. 123, 027002 (2019).
  * Cho  _et al._ [2020] C.-w. Cho, J. Shen, J. Lyu, O. Atanov, Q. Chen, S. H. Lee, Y. San Hor, D. J. Gawryluk, E. Pomjakushina, M. Bartkowiak, M. Hecker, J. Schmalian, and R. Lortz, $Z_{3}$-vestigial nematic order due to superconducting fluctuations in the doped topological insulator $Nb_{x}Bi_{2}Se_{3}$ and $Cu_{x}Bi_{2}Se_{3}$, Nature communications 11, 1 (2020).
  * Ghosh  _et al._ [2020] S. K. Ghosh, M. Smidman, T. Shang, J. F. Annett, A. D. Hillier, J. Quintanilla, and H. Yuan, Recent progress on superconductors with time-reversal symmetry breaking, Journal of Physics: Condensed Matter 33, 033001 (2020).
  * Naaman  _et al._ [2001] O. Naaman, W. Teizer, and R. C. Dynes, Fluctuation dominated josephson tunneling with a scanning tunneling microscope, Phys. Rev. Lett. 87, 097004 (2001).
  * Bastiaans  _et al._ [2019] K. M. Bastiaans, D. Cho, D. Chatzopoulos, M. Leeuwenhoek, C. Koks, and M. P. Allan, Imaging doubled shot noise in a josephson scanning tunneling microscope, Phys. Rev. B 100, 104506 (2019).
  * Yang  _et al._ [2013] H. U. Yang, E. Hebestreit, E. E. Josberger, and M. B. Raschke, A cryogenic scattering-type scanning near-field optical microscope, Review of Scientific Instruments 84, 023701 (2013).
  * McLeod  _et al._ [2017] A. S. McLeod, E. van Heumen, J. G. Ramirez, S. Wang, T. Saerbeck, S. Guenon, M. Goldflam, L. Anderegg, P. Kelly, A. Mueller, M. K. Liu, I. K. Schuller, and D. N. Basov, Nanotextured phase coexistence in the correlated insulator $V_{2}O_{3}$, Nature Physics 13, 80 (2017).

Supplementary Material: Local condensation of charge-$4e$ superconductivity at
a nematic domain wall Matthias Hecker Rafael M. Fernandes August 31, 2024

Supplementary Material: Local condensation of charge-$4e$ superconductivity at
a nematic domain wall

##  SI Derivation of the variational free energy

Here, we derive the variational free energy associated with the ansatz
$\mathcal{S}_{0}$ of the main text, i.e. we evaluate equation (3) of the main
text. For convenience, we first repeat the setup of the problem and the
notations introduced in the main text. Expressing the two-component
superconducting order parameter $\boldsymbol{\Delta}=(\Delta_{1},\Delta_{2})$
in terms of the four-component Nambu basis
$\hat{\boldsymbol{\Delta}}=\left(\boldsymbol{\Delta},\bar{\boldsymbol{\Delta}}\right)$,
we can rewrite the real-valued and complex-valued bilinear combinations as

$\displaystyle\Psi^{n}$
$\displaystyle=\hat{\boldsymbol{\Delta}}^{\dagger}M^{n}\hat{\boldsymbol{\Delta}},$
$\displaystyle\psi^{n}$
$\displaystyle=\hat{\boldsymbol{\Delta}}^{\dagger}m^{n}\hat{\boldsymbol{\Delta}}.$
(S1)

Here, we defined the matrices

$\displaystyle M^{A_{1g}}$ $\displaystyle=\tau^{0}\sigma^{0}/2,$
$\displaystyle m^{A_{1g}}$ $\displaystyle=\tau^{0}\sigma^{-},$ $\displaystyle
M^{A_{2g}}$ $\displaystyle=\tau^{y}\sigma^{z}/2,$ $\displaystyle M^{B_{1g}}$
$\displaystyle=\tau^{z}\sigma^{0}/2,$ $\displaystyle m^{B_{1g}}$
$\displaystyle=\tau^{z}\sigma^{-},$ $\displaystyle M^{B_{2g}}$
$\displaystyle=\tau^{x}\sigma^{0}/2,$ $\displaystyle m^{B_{2g}}$
$\displaystyle=\tau^{x}\sigma^{-},$ (S2)

with $\sigma^{\pm}=(\sigma^{x}\pm\mathsf{i}\sigma^{y})/2$ and $\tau^{i}$,
$\sigma^{i}$ acting respectively on the internal superconducting subspace and
on the Nambu space. For our specific setting of a one-dimensional grid of
length $L$ with lattice sites labeled $i,j=1,\dots,N$, the superconducting
action becomes

$\displaystyle\mathcal{S}$
$\displaystyle=\frac{L}{T}\sum_{i,j}\boldsymbol{\Delta}_{i}^{\dagger}\left[r_{0}\delta_{ij}+\frac{1}{2}f_{ij}^{0}\right]\tau^{0}\boldsymbol{\Delta}_{j}\;+\;\mathcal{S}^{\mathrm{int}}\,,$

where we introduced the gradient term

$\displaystyle\mathcal{S}^{\mathrm{grad}}$
$\displaystyle=\frac{L}{2T}\sum_{i,j}\boldsymbol{\Delta}_{i}^{\dagger}\tau^{0}f_{ij}^{0}\boldsymbol{\Delta}_{j},$
(S3)

described in terms of the hopping function
$f_{ij}^{0}=\frac{t_{0}}{2}\left(2\delta_{ij}-\delta_{i,j+1}-\delta_{i,j-1}\right)$
and the stiffness parameter $t_{0}>0$. Recall that $r_{0}=a_{0}(T-T_{0})$
denotes the bare superconducting transition temperature with $a_{0},T_{0}>0$.
The interaction part is given by [43] (see also Eq. (4) in main text)

$\displaystyle\mathcal{S}^{\mathrm{int}}$
$\displaystyle=\frac{L}{T}\sum_{i}\Big{[}u\,(\Psi_{i}^{A_{1g}})^{2}+v\,(\Psi_{i}^{A_{2g}})^{2}+w\,(\Psi_{i}^{B_{1g}})^{2}\Big{]},$
(S4)

where the Landau parameters satisfy the conditions $u>0$ and $v,w>-u$, in
order for the action to be bounded.

As discussed in the main text, within the Gaussian variational approach, we
choose a trial action

$\displaystyle\mathcal{S}_{0}$
$\displaystyle=\frac{1}{2}\frac{L}{T}\sum_{i,j}\hat{\boldsymbol{\Delta}}_{i}^{\dagger}\,\mathcal{G}_{i,j}^{-1}\,\hat{\boldsymbol{\Delta}}_{j},$
(S5)

that is characterized by the inverse Green’s function

$\displaystyle\mathcal{G}_{ij}^{-1}$
$\displaystyle=G_{i}^{-1}\,\delta_{ij}+f_{ij}^{0}M^{A_{1g}},$ (S6)
$\displaystyle G_{i}^{-1}$
$\displaystyle=2R_{i}M^{A_{1g}}+2\sum_{n\in\mathbb{G}_{\mathbb{R}}}\Phi_{i}^{n}M^{n}+\sum_{n\in\mathbb{G}_{\mathbb{C}}}\left(\bar{\phi}_{i}^{n}m^{n}+\mathrm{h.c.}\right),$
(S7)

which contains all variational parameters. Recall that we use $\Phi_{i}^{n}$
for real-valued variational composite order parameters and
$(\phi_{i}^{n},\bar{\phi}_{i}^{n})$ for the complex-valued ones, where $n$
denotes the irreducible representation (IR) according to which the composite
transforms. Moreover, we also define the mass renormalization parameter
$R_{i}=r_{0}+\Phi_{i}^{A_{1g}}$. For convenience of notation, we introduce the
IR sets $\mathbb{G}_{\mathbb{R}}=\\{A_{2g},B_{1g},B_{2g}\\}$ and
$\mathbb{G}_{\mathbb{C}}=\\{A_{1g},B_{1g},B_{2g}\\}$ in Eq. (S7). Note that
the local inverse Green’s function (S7) is identical to Eq. (5) in the main
text.

Since the trial action (S5) is Gaussian, it is straightforward to evaluate the
variational free energy [Eq. (6) in the main text]

$\displaystyle F_{v}$ $\displaystyle=-T\log
Z_{0}+T\langle\mathcal{S}-\mathcal{S}_{0}\rangle_{0},$ (S8)

with $Z_{0}=\int
D(\boldsymbol{\Delta},\bar{\boldsymbol{\Delta}})\,e^{-\mathcal{S}_{0}}$ being
the partition function of the trial action, see for example Ref. [43] for
technical details. In the real-space representation (S5), it is convenient to
promote the $4$-component field $\hat{\boldsymbol{\Delta}}_{i}$ to a
$(4N)$-component field via
$\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\boldsymbol{\Delta}}}}=\sum_{i}\hat{P}_{i}\hat{\boldsymbol{\Delta}}_{i}$,
where the projector $\hat{P}_{i}$ is a $(4N)\times 4$ dimensional matrix whose
elements are either $0$ or $1$, and
$\hat{P}_{i}^{T}\hat{P}_{j}=\delta_{ij}\mathbbm{1}_{4}$. The resulting
variational free energy $F_{v}$, up to an unimportant constant, is given by

$\displaystyle F_{v}$
$\displaystyle=\frac{T}{2}\log\det\Big{(}\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}^{-1}\Big{)}+T\\!\sum_{i}\\!\Big{\\{}2\big{[}r_{0}-R_{i}\\!+\tilde{U}_{A_{1g}}G_{ii}^{A_{1g}}\big{]}G_{ii}^{A_{1g}}$
$\displaystyle-2\sum_{n\in\mathbb{G}_{\mathbb{R}}}\left(\Phi_{i}^{n}\\!-\tilde{U}_{n}G_{ii}^{n}\right)G_{ii}^{n}-\sum_{n\in\mathbb{G}_{\mathbb{C}}}\big{[}\left(\phi_{i}^{n}-\tilde{u}_{n}g_{ii}^{n}\right)\bar{g}_{ii}^{n}+\mathrm{c.c.}\big{]}\\!\Big{\\}}.$
(S9)

Here, $\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}$ is
the inverse of
$\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}^{-1}=\sum_{i,j}\hat{P}_{i}\mathcal{G}_{ij}^{-1}\hat{P}_{j}^{T}$
and $G_{ij}^{n}$ is given by the decomposition of
$\mathcal{G}_{ij}=\hat{P}_{i}^{T}\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}\hat{P}_{j}$
onto the symmetry channels according to

$\displaystyle\mathcal{G}_{ij}$
$\displaystyle=2\\!\\!\\!\sum_{n\in\\{A_{1g},\mathbb{G}_{\mathbb{R}}\\}}\\!\\!G_{ij}^{n}M^{n}+\\!\\!\sum_{n\in\mathbb{G}_{\mathbb{C}}}\\!\\!\Big{[}g_{ij}^{n}(m^{n})^{\dagger}\\!+\bar{g}_{ij}^{n}m^{n}\Big{]}.$
(S10)

Finally, in Eq. (S9), the effective interaction parameters
$\left\\{\tilde{U}_{n},\tilde{u}_{n}\right\\}=\frac{T}{L}\left\\{U_{n},u_{n}\right\\}$
in each symmetry channel are given by

$\displaystyle U_{A_{1g}}$ $\displaystyle=3u+v+w,$ $\displaystyle u_{A_{1g}}$
$\displaystyle=u-v+w,$ $\displaystyle U_{A_{2g}}$ $\displaystyle=u+3v-w,$
$\displaystyle U_{B_{1g}}$ $\displaystyle=u-v+3w,$ $\displaystyle u_{B_{1g}}$
$\displaystyle=u+v+w,$ $\displaystyle U_{B_{2g}}$ $\displaystyle=u-v-w,$
$\displaystyle u_{B_{2g}}$ $\displaystyle=u+v-w,$ (S11)

which agree with the expressions found in the bulk case [43].

##  SII details of the free energy minimization

In this section, we present a few technical details related to the
minimization of the free energy (S9). As emphasized in the main text, we model
the domain wall through the boundary conditions
$\Phi_{1}^{B_{1g}}=-\Phi_{N}^{B_{1g}}=\Phi_{0}^{B_{1g}}$, $R_{1}=R_{N}=R_{0}$,
and all other fields $\Phi_{1,N}^{n}=\phi_{1,N}^{n}=0$, where
$\Phi_{0}^{B_{1g}}$ and $R_{0}$ are the corresponding bulk values.

The free energy (S9) effectively depends only on the parameters $T$,
$\hat{t}_{0}\equiv t_{0}/\sqrt{uT/L}$, $\\{v,w\\}/u$, and
$\\{r_{0},R_{i},\Phi_{i}^{n},\phi_{i}^{n}\\}/t_{0}$. In the spirit of the
Ginzburg-Landau approach, we assume the most important temperature dependence
to be in $r_{0}=a_{0}(T-T_{0})$, and we set $T=T_{0}$ elsewhere, with $T_{0}$
denoting the bare SC transition temperature. To facilitate the minimization of
the free energy (S9), we supply the minimizer with the gradient expressions
given by

$\displaystyle\frac{\partial F_{v}}{\partial X_{i}}$
$\displaystyle=T\\!\sum_{j}\\!\Bigg{\\{}V_{j}^{A_{1g}}\frac{\partial
G_{jj}^{A_{1g}}}{\partial
X_{i}}+\sum_{n\in\mathbb{G}_{\mathbb{R}}}V_{j}^{n}\frac{\partial
G_{jj}^{n}}{\partial X_{i}}$
$\displaystyle\quad+\sum_{n\in\mathbb{G}_{\mathbb{C}}}\left(v_{j}^{n}\frac{\partial
g_{jj}^{n}}{\partial X_{i}}+\mathrm{c.c.}\right)\Bigg{\\}},$ (S12)

where $X_{i}\in\\{R_{i},\Phi_{i}^{n},\phi_{i}^{n},\bar{\phi}_{i}^{n}\\}$ can
be any of the variational parameters. In writing Eq. (S12), we exploited the
fact that the partial derivatives $\frac{\partial F_{v}}{\partial
X_{i}}\big{|}_{G^{n},g^{n}}=0$ vanish [43]. Moreover, we defined

$\displaystyle V_{j}^{A_{1g}}$
$\displaystyle=2\big{[}r_{0}-R_{j}\\!+2\tilde{U}_{A_{1g}}G_{jj}^{A_{1g}}\big{]},$
$\displaystyle v_{j}^{n}$
$\displaystyle=-\left(\bar{\phi}_{j}^{n}-2\tilde{u}_{n}\bar{g}_{jj}^{n}\right),$
$\displaystyle V_{j}^{n}$
$\displaystyle=-2\left(\Phi_{j}^{n}\\!-2\tilde{U}_{n}G_{jj}^{n}\right).$ (S13)

To determine the remaining derivatives in Eq. (S12), we use the Green’s
function (S10) to identify

$\displaystyle G_{jj}^{n}$
$\displaystyle=\frac{1}{2}\mathrm{tr}\left[\mathcal{G}_{jj}M^{n}\right],$
$\displaystyle g_{jj}^{n}$
$\displaystyle=\frac{1}{2}\mathrm{tr}\left[\mathcal{G}_{jj}m^{n}\right].$
(S14)

Then, using the introduced projector matrix $\hat{P}_{i}$, we obtain the
relationship

$\displaystyle\frac{\partial\mathcal{G}_{jj}}{\partial X_{i}}$
$\displaystyle=-\hat{P}_{j}^{T}\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}\frac{\partial\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}^{-1}}{\partial
X_{i}}\hat{\vphantom{\rule{1.0pt}{6.57643pt}}\smash{\hat{\mathcal{G}}}}\hat{P}_{j}=-\sum_{i_{1}i_{2}}\mathcal{G}_{ji_{1}}\frac{\partial\mathcal{G}_{i_{1}i_{2}}^{-1}}{\partial
X_{i}}\mathcal{G}_{i_{2}j},$ (S15)

which leads to the expressions

$\displaystyle\frac{\partial\mathcal{G}_{jj}}{\partial R_{i}}$
$\displaystyle=-2\mathcal{G}_{ji}M^{A_{1g}}\mathcal{G}_{ij},$
$\displaystyle\frac{\partial\mathcal{G}_{jj}}{\partial\Phi_{i}^{n}}$
$\displaystyle=-2\mathcal{G}_{ji}M^{n}\mathcal{G}_{ij},$
$\displaystyle\frac{\partial\mathcal{G}_{jj}}{\partial\bar{\phi}_{i}^{n}}$
$\displaystyle=-\mathcal{G}_{ji}m^{n}\mathcal{G}_{ij}.$ (S16)

It is now straightforward to compute the derivatives in Eq. (S14) by using Eq.
(S16); we find, for instance

$\displaystyle\frac{\partial G_{jj}^{n}}{\partial R_{i}}$
$\displaystyle=-\mathrm{tr}\left[\mathcal{G}_{ji}M^{A_{1g}}\mathcal{G}_{ij}M^{n}\right],$
$\displaystyle\frac{\partial g_{jj}^{n}}{\partial R_{i}}$
$\displaystyle=-\mathrm{tr}\left[\mathcal{G}_{ji}M^{A_{1g}}\mathcal{G}_{ij}m^{n}\right],$
(S17)

and similar expressions for the other variational parameters.