# Designing Phase Sensitive Probes of Monopole Superconducting Order

Grayson R. Frazier  Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, Maryland 21218, USA  Junjia Zhang  Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, Maryland 21218,
USA  Junyi Zhang  Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, Maryland 21218, USA  Xinyu Sun  Department of Physics
and Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA  Yi Li
Department of Physics and Astronomy, Johns Hopkins University, Baltimore,
Maryland 21218, USA

(July 26, 2024)

###### Abstract

Distinct from familiar $s$-, $p$-, or $d$-wave pairings, the monopole
superconducting order represents a novel class of pairing order arising from
nontrivial monopole charge of the Cooper pair. In the weak-coupling regime,
this order can emerge when pairing occurs between Fermi surfaces with
different Chern numbers in, for example, doped Weyl semimetal systems.
However, the phase of monopole pairing order is not well-defined over an
entire Fermi surface, making it challenging to design experiments sensitive to
both its symmetry and topology. To address this, we propose a scheme based on
symmetry and topological principles to identify this elusive pairing order
through a set of phase-sensitive Josephson experiments. By examining the
discrepancy between global and local angular momentum of the pairing order, we
can unveil the monopole charge of the pairing order. We demonstrate the
proposed probe of monopole pairing order through analytic and numerical
studies of Josephson coupling in models of monopole superconductor junctions.
This work opens a promising avenue to uncover the unique topological
properties of monopole pairing orders and to distinguish them from known
pairing orders based on spherical harmonic symmetry.

##  I Introduction

Monopole superconducting order [1] is an exotic family of three-dimensional
topological superconducting orders. This pairing order arises from the
nontrivial geometric phase of Cooper pairs [2] formed between two Fermi
surfaces with different Chern numbers, which can occur in weak-coupling
superconducting materials supporting Weyl Fermi surfaces [3, 4, 5, 6]. For
example, when superconductivity is developed in a three-dimensional inversion
symmetric magnetic Weyl system, intrinsically or by proximity, monopole
harmonic pairing exists as the zero center-of-mass momentum pairing between
two inversion-related Fermi surfaces enclosing Weyl points of opposite
chiralities [1, 7, 8]. Unlike many other known topological superconductors,
where only Bogoliubov-de Gennes (BdG) quasi-particle states cannot be
adiabatically connected to a Bose-Einstein condensate [9, 10], the pairing
order of a monopole superconductor also fails to be globally well-defined over
the entire Fermi surface due to topological obstruction in the $\mathrm{U}(1)$
pairing phase. The topological obstruction in the pairing order also
fundamentally changes the symmetry representation: monopole pairing order is
no longer describable by familiar symmetry in terms of $s$-, $p$-, and
$d$-waves based on spherical harmonics but is rather characterized by the so-
called monopole harmonic functions [11, 12, 1, 7]. Moreover, the existence of
pairing nodes in monopole superconducting order is independent of specific
pairing mechanisms and is instead dictated by nonzero monopole charge of
Cooper pairs, which results in unconventional pairing with nontrivial winding
of pairing phase over an entire Fermi surface [2, 1].

It has been known that, in the presence of a magnetic monopole in real space,
the orbital angular momentum of a charged particle is no longer conserved.
However, the total angular momentum, which includes the angular momentum of
the electromagnetic field, remains conserved in both classical and quantum
scatterings of charged particles [13, 14]. In a monopole superconductor,
Cooper pairs scatter in the presence of a momentum space “magnetic” monopole.
The total angular momentum of the pairing order, which we refer to as the
“global angular momentum,” includes the contribution from the pair Berry flux
and is a conserved physical observable that characterizes the global symmetry
of the pairing order. However, even though the global angular momentum is
conserved, the pairing order is not well-defined globally, i.e. for any choice
of gauge, there exists at least a singular point at the Fermi surface at which
the pairing phase cannot be well-defined. Only away from the singularity is
the pairing order locally well-defined and able to be described in terms of
spherical harmonic functions, which determine the “local angular momentum” of
the pairing order. The shift, or say, the discrepancy, between the conserved
global angular momentum and local angular momentum is a unique signature of
the monopole charge of the pairing order.

To detect superconducting pairing symmetries, phase-sensitive experiments have
long played an essential role [15, 16, 17, 18, 19, 20, 21]. For example, the
$d_{x^{2}-y^{2}}$ pairing symmetry of a high-temperature cuprate
superconductor has been revealed by a corner Josephson junction. The relative
sign difference of the pairing order along the $a$ and $b$ axes of a cuprate
superconductor leads to a destructive Fraunhofer interference pattern at zero
magnetic flux [20]. For $p$-wave topological superconductors, Majorana modes
localized at boundaries have been proposed to contribute to single particle
tunneling and the $4\pi$-periodic Josephson current phase relation [22, 23,
24, 25, 26, 27, 28, 29]. For monopole pairing order, the pairing phase is not
well-defined over an entire Fermi surface, making it challenging to design
experiments sensitive to both its symmetry and topology.

We propose a phase-sensitive approach to address the outstanding question of
experimentally identifying the monopole superconducting order. Guided by
symmetry and topological principles, we develop the design of a set of
Josephson junctions that can extract the shift between the global and local
angular momentum of the monopole pairing order, thereby determining the pair
monopole charge. Specifically, we investigate two classes of Josephson
junctions and study the Josephson currents in these junctions by combining
linear response and numerical tight-binding calculations. The first class
consists of junctions aligned along the common high-symmetry rotational axis
between a monopole superconductor and a superconductor with known spherical
harmonic pairing symmetry. This junction probes the global angular momentum of
the monopole pairing order, as only when the angular momentum components at
two sides of the junction are identical does the system exhibit nonzero first-
order Josephson coupling. The second class consists of junctions between
identical monopole superconductors, probing the local pairing symmetry at the
momenta conserved at the junction interface. Our investigation provides
guiding principles for designing phase-sensitive experiments to probe the
exotic monopole pairing order and distinguish it from all known
superconducting orders with spherical harmonic symmetry.

This paper is organized as follows. In Sec. II, we discuss the first-order
Josephson coupling between two uniform superconductors based on linear
response theory, initially considering spin-orbit interactions at the
interface and then extending to spin-orbit coupling in the bulk, emphasizing
the role of the symmetry of superconducting orders in terms of a scattering
form factor in this microscopic formulation. In Sec. III, we focus on
Josephson junctions involving monopole superconductors, first deriving the
form factor discussed in Sec. II based on linear response theory, and then
analyzing the symmetry and topological principles which are independent of
microscopic tunneling processes. We propose designing Josephson junctions
aligned along or perpendicular to the high-symmetry rotational axis of the
monopole superconductor to reveal the shift between the global and local
angular momentum of the pairing order and extract the pair monopole charge. In
Sec. IV, we design Josephson junctions to extract the global angular momentum
of the monopole superconducting order, presenting results from both linear
response theory and numerical results of a tight-binding model. Sec. V
explores the design of Josephson junctions to probe the local angular momentum
component by considering junctions between identical monopole superconductors.
Finally, Sec. VI addresses the presence of spin-orbit interactions at the
junction interface, proposing additional junction designs to distinguish
monopole superconducting order from chiral spherical harmonic pairing orders.

##  II First order Josephson coupling in the presence of spin-orbit coupling

In this section, we employ linear response theory to derive the first order
Josephson coupling [30, 31] in a superconductor-insulator-superconductor (SIS)
junction in the presence of spin-orbit coupling in the bulk of the
superconductors. Initially, we examine the scenario where only spin-orbit
interactions occur at the junction interface, assuming the absence of spin-
orbit coupling in the bulk. We emphasize the role of a form factor, which
determines the symmetry selection rules for the first order Josephson current.
Then, we additionally consider spin-orbit coupling in the bulk of
superconductors. In the weak-coupling regime, by projecting to states near
Fermi surfaces, we demonstrate that the bulk spin-orbit coupling can be
treated as an effective interface spin-orbit interaction. Using this
formalism, we can obtain the first order Josephson current in the presence of
spin-orbit interactions in the bulk superconductors or at the interface.

We restrict our consideration to the Josephson effects between two uniform
superconductors separated by a thin insulating barrier, as described by the
following Hamiltonian,

$H_{\mathrm{JJ}}=H_{\mathrm{BdG},L}+H_{\mathrm{BdG},R}+H_{\mathrm{link}}.$ (1)

Here, we use subscripts $L$ and $R$ to denote the two sides of a junction,
irrespective of its spatial orientation. The Bogoliubov-de-Gennes (BdG)
Hamiltonians $H_{\mathrm{BdG},\alpha=L,R}$ of the superconductors on opposite
sides of the junction generally take the form

$\displaystyle H_{\mathrm{BdG},\alpha}=$ (2)
$\displaystyle{\sum_{\mathbf{k}}}^{\prime}\Psi^{\dagger}_{\alpha}(\mathbf{k})\left(\begin{array}[]{cc}\mathcal{H}_{\mathrm{kin},\alpha}(\mathbf{k})&\Delta_{\alpha}(\mathbf{k})e^{i\phi_{\alpha}}\\\
\Delta^{\dagger}_{\alpha}(\mathbf{k})e^{-i\phi_{\alpha}}&-\mathcal{H}^{T}_{\mathrm{kin},\alpha}(-\mathbf{k})\end{array}\right)\Psi_{\alpha}(\mathbf{k}),$

where ${\sum_{\mathbf{k}}}^{\prime}$ denotes the sum over half the Brillouin
zone and
$\Psi_{\alpha}(\mathbf{k})=(c_{\alpha,\mathbf{k},\uparrow},c_{\alpha,\mathbf{k},\downarrow},c_{\alpha,-\mathbf{k},\uparrow}^{\dagger},c_{\alpha,-\mathbf{k},\downarrow}^{\dagger})^{\mathrm{T}}$
denotes the annihilation operator of Nambu spinor. Here,
$c_{\alpha=L(R),\mathbf{k},\sigma}$ annihilates a single-particle state with
momentum $\mathbf{k}$ and spin $\sigma=\uparrow,\downarrow$ at the left
(right) side of the junction. $\mathcal{H}_{\mathrm{kin},\alpha}$ is the
kinetic Hamiltonian kernel defined in the particle-hole channel, and
$\Delta_{\alpha}e^{i\phi_{\alpha}}$ is the superconducting pairing matrix
defined in the particle-particle channel with the overall $\mathrm{U}(1)$
pairing phase $e^{i\phi_{\alpha}}$ separated out here. $H_{\mathrm{link}}$
describes the single particle tunnelings across the insulating junction link
and takes the form of the Bardeen-Josephson Hamiltonian [32, 33, 15] as
follows,

$H_{\mathrm{link}}=\sum_{\mathbf{k},\mathbf{k}^{\prime}}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}T_{\sigma,\sigma^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})c^{\dagger}_{R,\mathbf{k},\sigma}c_{L,\mathbf{k}^{\prime},\sigma^{\prime}}+\mathrm{h.c.}$
(3)

The tunneling amplitudes
$T_{\sigma,\sigma^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})$ form a $2\times
2$ matrix $T(\mathbf{k},\mathbf{k}^{\prime})$ which can be decomposed into
spin-independent and spin-dependent tunneling amplitudes as
$T(\mathbf{k},\mathbf{k}^{\prime})=T_{0}(\mathbf{k},\mathbf{k}^{\prime})\sigma_{0}+{\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})\cdot\bm{\sigma}}$,
where $\sigma_{0}$ and $\sigma_{i=x,y,z}$ are the identity and Pauli matrices
in spin space. $\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})$, the spin-
dependent tunneling amplitude, results from the spin-orbit interactions at the
junction interface. The results are also generalizable to systems with
additional internal degrees of freedom, like valley or orbital.

We now derive the first order Josephson current in this SIS junction at zero
bias, where the dissipationless DC Josephson current running in thermal
equilibrium is the only contributor to the tunneling current through the
junction. Using linear response theory, to first order in the perturbation
$H_{\mathrm{link}}$ in the expansion of the thermal weight $e^{-\beta
H_{\mathrm{link}}}$ with $\beta=1/(k_{B}T)$, the DC Josephson current is given
by

$\displaystyle I_{J}(\phi)=\mathrm{lim}_{\delta\rightarrow 0^{+}}$
$\displaystyle\frac{2e}{\hbar}\mathrm{Im}\Big{\\{}\frac{1}{\beta}\sum_{ip_{n}}\sum_{\mathbf{k},\mathbf{k}^{\prime}}\mathrm{Tr}[\mathcal{F}_{L}(\mathbf{k},ip_{n}-i\delta)$
$\displaystyle
T^{\mathrm{T}}(-\mathbf{k},-\mathbf{k}^{\prime})\mathcal{F}^{*}_{R}(\mathbf{k},ip_{n})T(\mathbf{k},\mathbf{k}^{\prime})]\Big{\\}},$
(4)

where $\phi=\phi_{R}-\phi_{L}$ is the overall U($1$) phase difference of
superconducting pairing across the junction, the trace is taken over spin
space, and $p_{n}=(2n+1)\pi/\beta$ with $n\in\mathbb{Z}$ are fermion Matsubara
frequencies. $\mathcal{F}_{\alpha}(\mathbf{k},ip_{n})$ is the anomalous
Green’s function of the unperturbed BdG Hamiltonian $H_{\mathrm{BdG},\alpha}$
in Eq. (2). Generally, the total first order Josephson current can be
expressed as

$I_{J}(\phi)=I_{c}\sin(\phi),$ (5)

in which $I_{c}$ is the Josephson critical current at zero bias.

When the band Hamiltonian $\mathcal{H}_{\mathrm{kin},\alpha}$ has no spin
dependence and when the pairing is unitary,
$\Delta^{\dagger}\Delta\propto\sigma_{0}$, the anomalous Green’s function
$\mathcal{F}_{\alpha}(\mathbf{k},ip_{n})$ can be simplified to
$-e^{i\phi_{\alpha}}\Delta_{\alpha}(\mathbf{k})/(p_{n}^{2}+E^{2}_{\alpha,\mathbf{k}})$.
Here,
$E_{\alpha,\mathbf{k}}=\sqrt{\xi^{2}_{\alpha,\mathbf{k}}+\mathrm{Tr}(\Delta_{\alpha}^{\dagger}(\mathbf{k})\Delta_{\alpha}(\mathbf{k}))/2}$
is the dispersion of the unperturbed BdG quasiparticle on the $\alpha$ side of
the junction with $\xi_{\alpha,\mathbf{k}}$ being the spin-degenerate band
dispersion. The form of $\mathcal{F}(\mathbf{k},ip_{n})$ for more general
pairings, including non-unitary spin triplet pairing, is shown in Appendix
A.1.

The trace part of Eq. (4) can be regarded as the following product of three
distinct factors,

$\displaystyle\frac{1}{\beta}\sum_{ip_{n}}\mathrm{Tr}[\cdots]=e^{-i\phi}w\left(E_{L,\mathbf{k}^{\prime}},E_{R,\mathbf{k}};\beta\right)\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime}).$
(6)

Here, the first term $e^{-i\phi}$ gives rise to $2\pi$-periodicity with
respect to the pairing phase difference $\phi$ in the Josephson current
$I_{J}(\phi)$. This reflects the first order process corresponding to a single
Cooper pair tunneling through the junction barrier. The second term
$w(E_{L;\mathbf{k}^{\prime}},E_{R,\mathbf{k}};\beta)\equiv\frac{1}{\beta}\sum_{ip_{n}}(p_{n}^{2}+E^{2}_{L,\mathbf{k}^{\prime}})^{-1}(p_{n}^{2}+E^{2}_{R,\mathbf{k}})^{-1}$
includes the temperature dependence and contribution from the unperturbed BdG
quasiparticle energy dispersions, $E_{\alpha,\mathbf{k}}$. After frequency
summation, it becomes

$\displaystyle w(E_{L;\mathbf{k}^{\prime}},E_{R,\mathbf{k}};\beta)$ (7)
$\displaystyle=$
$\displaystyle\frac{1}{4E_{L,\mathbf{k}^{\prime}}E_{R,\mathbf{k}}}\mathrm{lim}_{\omega\rightarrow
0^{+}}[(n_{F}(E_{L,\mathbf{k}^{\prime}})-n_{F}(E_{R,\mathbf{k}}))\times$
$\displaystyle(\frac{1}{E_{L,\mathbf{k}^{\prime}}-E_{R,\mathbf{k}}-i\omega}+\frac{1}{E_{L,\mathbf{k}^{\prime}}-E_{R,\mathbf{k}}+i\omega})$
$\displaystyle+(1-n_{F}(E_{L,\mathbf{k}^{\prime}})-n_{F}(E_{R,\mathbf{k}}))\times$
$\displaystyle(\frac{1}{E_{L,\mathbf{k}^{\prime}}+E_{R,\mathbf{k}}+i\omega}+\frac{1}{E_{L,\mathbf{k}^{\prime}}+E_{R,\mathbf{k}}-i\omega})],$

where $n_{F}$ is the Fermi-Dirac distribution.

The last and the most essential term is the form factor
$\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})$, defined as

$\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})\equiv\mathrm{Tr}\Big{[}{\Delta_{L}(\mathbf{k}^{\prime})}{T^{\mathrm{T}}(-\mathbf{k},-\mathbf{k}^{\prime})}{\Delta_{R}^{*}(\mathbf{k})}{T(\mathbf{k},\mathbf{k}^{\prime})}\Big{]}.$
(8)

Though the magnitude of the form factor depends on the microscopic details of
tunneling at the junction interface via $T(\mathbf{k},\mathbf{k}^{\prime})$,
whether it vanishes or not is sensitive to the superconducting orders at the
two sides of the junction, $\Delta_{L}$ and $\Delta_{R}$, as well as the
symmetry of the junction geometry. Therefore, the form factor
$\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})$ gives rise to symmetry
selection rules which determine whether the first order Josephson current
vanishes.

To demonstrate the symmetry selection rule provided by the form factor,
consider a Josephson junction between a spin singlet and a spin triplet
superconductor without spin-orbit coupling in the bulk superconductors but
with spin-dependent tunneling at the junction interface. The pairing matrices
of superconductors at two sides of the junction generally take the form
$\Delta_{L}(\mathbf{k})=d_{L,0}(\mathbf{k})i\sigma_{y}$ and
$\Delta_{R}(\mathbf{k})=(\mathbf{d}_{R}(\mathbf{k})\cdot\bm{\sigma})i\sigma_{y}$,
with $d_{L,0}(\mathbf{k})$ the spatial part of the singlet pairing order and
$\mathbf{d}_{R}(\mathbf{k})$ the $d$-vector of the triplet pairing order [34,
35]. Employing trace identities of Pauli matrices, we simplify the form factor
for this singlet-triplet junction, denoted as
$\mathfrak{F}_{\mathrm{sing}-\mathrm{trip}}$, to the following form,

$\displaystyle\mathfrak{F}_{\mathrm{sing}-\mathrm{trip}}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle=$ $\displaystyle
2d_{L,0}(\mathbf{k}^{\prime})\\{T_{0}(-\mathbf{k},-\mathbf{k}^{\prime})\mathbf{d}^{*}_{R}(\mathbf{k})\cdot\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})$
(9)
$\displaystyle-i\mathbf{d}^{*}_{R}(\mathbf{k})\cdot[\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})\times\mathbf{T}(-\mathbf{k},-\mathbf{k}^{\prime})]$
$\displaystyle-
T_{0}(\mathbf{k},\mathbf{k}^{\prime})\mathbf{d}^{*}_{R}(\mathbf{k})\cdot\mathbf{T}(-\mathbf{k},-\mathbf{k}^{\prime})\\}.$

Consequently, in the absence of spin-orbit interaction at junction interface,
$\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=0$, the form factor vanishes and
forbids first order Josephson coupling. This reflects the symmetry principle
that when pairing orders are in different total spin channels, the first order
Josephson coupling vanishes in the system where spin and orbital angular
momentum are separately conserved.

Further consider Rashba-type spin-orbit interaction at the junction interface,
as will be used later in Sec. VI, which has the form
$\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=T_{\mathrm{SO}}(\hat{\mathbf{n}}\times\hat{\mathbf{k}})\delta_{\mathbf{k},\mathbf{k}^{\prime}}$,
with $\hat{\mathbf{n}}$ being the junction orientation. The form factor
reduces to

$\displaystyle\mathfrak{F}_{\mathrm{sing}-\mathrm{trip}}(\mathbf{k},\mathbf{k}^{\prime})=4T_{0}T_{\mathrm{SO}}d_{L,0}(\mathbf{k}^{\prime})[\hat{\mathbf{n}}\cdot(\hat{\mathbf{k}}\times\mathbf{d}_{R}^{*}(\mathbf{k}))]\delta_{\mathbf{k},\mathbf{k}^{\prime}},$
(10)

where we take $T_{0}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}$, a constant, for
simplicity. Together with Eq. (4), the first order Josephson current takes the
form,
$I_{J}=(8e/\hbar)T_{0}T_{\mathrm{SO}}\mathrm{Im}[e^{-i\phi}\sum_{\mathbf{k}}w(E_{L,\mathbf{k}},E_{R,\mathbf{k}};\beta)d_{L,0}(\mathbf{k})\hat{\mathbf{n}}\cdot(\hat{\mathbf{k}}\times\mathbf{d}_{R}^{*}(\mathbf{k}))]$,
which agrees with the result obtained by Geshkenbein and Larkin [36]. Then,
when $d_{L,0}(\mathbf{k}^{\prime})=\Delta_{L,0}$ for an $s$-wave
superconductor,
$\mathbf{d}_{R}(\mathbf{k})=\Delta_{R,0}(0,0,(k_{x}+ik_{y})/k_{F})$ for a
chiral $p$-wave triplet superconductor with total spin $s_{z}=0$, and the
junction is oriented along $\hat{\mathbf{n}}=\hat{z}\parallel\mathbf{d}_{R}$,
the form factor in Eq. (10) vanishes, forbidding first order Josephson
coupling even though there is nonvanishing Rashba spin-orbit interaction. This
result can be understood from the following spin symmetry selection rules. In
the spin-independent tunneling channel, as the singlet and triplet
superconductors have different total spins, their first order Josephson
coupling is zero. In the spin-dependent tunneling channel, the junction
geometry determines the form of Rashba spin-orbit interaction as consisting of
only in-plane spins which change $s_{z}$, the total spin-$z$ component of the
Cooper pair, by $1$. However, $s_{z}$ of the pairing orders at two sides of
the junction are identical, $s_{z}=0$. Hence, the first order Josephson
coupling in this example vanishes.

In this work, to take into account spin-orbit couplings in bulk topological
superconductors, we introduce a convenient formalism to treat the spin-orbit
coupling in terms of effective spin-dependent tunneling at the junction
barrier. We consider the tunneling of band eigenstates through the junction,
where $H_{\mathrm{link}}$ in Eq. (3) now takes the form

$H_{\mathrm{link}}=\sum_{\mathbf{k},s;\mathbf{k}^{\prime},s^{\prime}}T_{s,s^{\prime}}^{(b)}(\mathbf{k},\mathbf{k}^{\prime})\psi^{\dagger}_{R,\mathbf{k},s}\psi_{L,\mathbf{k}^{\prime},s^{\prime}}+\mathrm{h.c.}$
(11)

Here, $\psi_{\alpha,\mathbf{k},s}$ with $\alpha=L(R)$ annihilates a band
eigenstate with band index $s=\pm$ and momentum $\mathbf{k}$ at the left
(right) of the junction. For a generic two-band Hamiltonian,
$\psi_{\alpha,\mathbf{k},s}$ is related to the annihilation operators in spin
basis, $c_{\alpha,\mathbf{k},\sigma}$, by a unitary transformation,
$\psi_{\alpha,\mathbf{k},s}=\sum_{\sigma=\uparrow,\downarrow}[U^{\dagger}_{\alpha}(\mathbf{k})]_{s,\sigma}c_{\alpha,\mathbf{k},\sigma}$,
where $U_{\alpha}(\mathbf{k})$ encodes the spin-orbit coupling in band
Hamiltonian $\mathcal{H}_{\mathrm{kin},\alpha}(\mathbf{k})$. It is convenient
to define the complex four-vector
$(u_{\alpha,0}(\mathbf{k}),\mathbf{u}_{\alpha}(\mathbf{k}))$, such that
$U_{\alpha}(\mathbf{k})=u_{\alpha,0}(\mathbf{k})\sigma_{0}+\mathbf{u}_{\alpha}(\mathbf{k})\cdot\bm{\sigma}$,
with $u_{\alpha,i=0,\cdots,3}(\mathbf{k})$ being complex functions.
Correspondingly, the tunneling matrix defined in the band representation,
denoted by the superscript “$b$”, is given by
$T^{(b)}(\mathbf{k},\mathbf{k}^{\prime})=U^{\dagger}_{R}(\mathbf{k})T(\mathbf{k},\mathbf{k}^{\prime})U_{L}(\mathbf{k}^{\prime})$
with $T(\mathbf{k},\mathbf{k}^{\prime})$ defined below Eq. (3) in the spin
basis. Therefore, even in the absence of bare spin-orbit interaction at the
junction interface, $\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=0$, the
effective tunneling in the band eigenbasis takes the following form,

$T^{(b)}(\mathbf{k},\mathbf{k}^{\prime})=T^{(b)}_{0}(\mathbf{k},\mathbf{k}^{\prime})+\mathbf{T}^{(b)}(\mathbf{k},\mathbf{k}^{\prime})\cdot\bm{\sigma},$
(12)

where
$T^{(b)}_{0}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}(\mathbf{k},\mathbf{k}^{\prime})(u^{*}_{R,0}(\mathbf{k})u_{L,0}(\mathbf{k}^{\prime})+\mathbf{u}^{*}_{R}(\mathbf{k})\cdot\mathbf{u}_{L}(\mathbf{k}^{\prime}))$
and
$\mathbf{T}^{(b)}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}(\mathbf{k},\mathbf{k}^{\prime})(u^{*}_{R,0}(\mathbf{k})\mathbf{u}_{L}(\mathbf{k}^{\prime})+\mathbf{u}^{*}_{R}(\mathbf{k}){u}_{L,0}(\mathbf{k}^{\prime})+i\mathbf{u}^{*}_{R}(\mathbf{k})\times\mathbf{u}_{L}(\mathbf{k}^{\prime}))$.
For most parts of the manuscript, we consider the case of zero bare interface
spin-orbit interaction for simplicity. We discuss the results when the bare
interface spin-orbit coupling $\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})\neq
0$ in Sec. VI.

Furthermore, in a weak-coupling superconductor, we can project its pairing
order to low energy states near the Fermi surfaces participating in Cooper
pairing. For pairing order already expressed in the band representation
$\Delta_{\alpha}^{(b)}(\mathbf{k})=U_{\alpha}^{\dagger}(\mathbf{k})\Delta_{\alpha}(\mathbf{k})U^{*}_{\alpha}(-\mathbf{k})$,
we capture the key ingredients of the low energy pairing by focusing on the
Fermi surface projected pairing order in this representation, defining
$\Delta^{(bp)}_{\alpha}(\mathbf{k})\equiv
P_{+}^{(b)}\Delta_{\alpha}^{(b)}(\mathbf{k})P_{+}^{(b)}$, where
$P_{+}^{(b)}=\mathrm{diag}(1,0)$ is the projection operator in band
representation. Here, without loss of generality, we have ordered the band
basis so that the states at the Fermi surface of interest correspond to the
first eigenvalues in the particle and hole parts (see Appendix B). The
superscript “$bp$” denotes the band diagonal basis with bands projected to the
Fermi surfaces of interest. The anomalous Green’s function of the
corresponding BdG Hamiltonian in this band projected basis takes a simple
form,
$\mathcal{F}^{(bp)}_{\alpha}(\mathbf{k},ip_{n})\approx-e^{i\phi_{\alpha}}\Delta_{\alpha}^{(bp)}(\mathbf{k})/(p_{n}^{2}+E^{(bp)2}_{\alpha,\mathbf{k}})$.
This form is analogous to the simplified form of
$\mathcal{F}_{\alpha}(\mathbf{k},ip_{n})$ in the case of spin degenerate band
and unitary pairing. However, here, the energy of the BdG quasiparticle is
given by
$E^{(bp)}_{\alpha,\mathbf{k}}=\sqrt{\xi_{\alpha,\mathbf{k}}^{(bp)2}+\mathrm{Tr}(\Delta^{(bp)\dagger}_{\alpha}(\mathbf{k})\Delta^{(bp)}_{\alpha}(\mathbf{k}))}$,
where $\xi^{(bp)}_{\alpha,\mathbf{k}}$ is the dispersion of the band that
forms the Fermi surface. The expression of the form factor in Eq. (8) now
reduces to
$\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})=\mathrm{Tr}[\Delta_{L}^{(bp)}(\mathbf{k}^{\prime})T^{(b)\mathrm{T}}(-\mathbf{k},-\mathbf{k}^{\prime})\Delta_{R}^{(bp)*}(\mathbf{k})T^{(b)}(\mathbf{k},\mathbf{k}^{\prime})]$,
where the pairing and tunneling matrices are in the projected band basis. We
use this band projected form factor to analyze Josephson junctions of monopole
superconductors, as demonstrated in the following section.

##  III Josephson coupling in junctions of monopole superconductors

Having discussed the form of the first order Josephson current in a general
SIS junction, our focus turns to junctions that include monopole
superconductors. This discussion unfolds in two different approaches. We first
derive the form factor that determines the Josephson current using the
microscopic derivation based on linear response theory, as discussed in Sec.
II. We work in the weak-coupling regime, projecting the pairing order to low
energy states that comprise the Fermi surface. Secondly, we supplement the
first approach with a discussion based on the underlying symmetry and
topological principles, which are essential in designing a set of Josephson
junctions to probe the monopole superconducting order.

###  III.1 Microscopic derivation of Josephson coupling for monopole
superconductors

We first extend the results from linear response theory in Sec. II to
determine the first order Josephson current for junctions that include
monopole superconductors. Monopole pairing order is an exotic class of
topologically obstructed pairing order that arises from the nontrivial Berry
phase of Cooper pairs. This, for example, can arise when pairing occurs
between two Fermi surfaces $\mathrm{FS}_{1}$ and $\mathrm{FS}_{2}$ with
different Chern numbers, $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ respectively,
in doped time-reversal symmetry broken Weyl semimetals in three dimensions
[1]. In the weak-coupling regime, we describe this system using mean-field BdG
Hamiltonian in Eq. (2), with its Nambu spinor modified to
$\Psi(\mathbf{k})=(c_{1;\mathbf{k},\uparrow},c_{1;\mathbf{k},\downarrow},c_{2;-\mathbf{k},\uparrow}^{\dagger},c_{2;-\mathbf{k},\downarrow}^{\dagger})^{\mathrm{T}}$
to describe the inter-Fermi surface pairing
$\Delta_{\mathrm{inter}}(\mathbf{k})$ between parity-related $\mathrm{FS}_{1}$
and $\mathrm{FS}_{2}$, where $c_{i;\mathbf{k},\sigma}$ denotes the
annihilation operator of an electron on $\mathrm{FS}_{i}$ ($i=1,2$) with
momentum $\mathbf{k}$ and spin $\sigma=\uparrow,\downarrow$. Correspondingly,
the BdG Hamiltonian kernel takes the following form,

$\mathcal{H}_{\mathrm{BdG}}(\mathbf{k})=\left(\begin{array}[]{cc}\mathcal{H}_{\mathrm{kin},1}(\mathbf{k})&\Delta_{\mathrm{inter}}(\mathbf{k})\\\
\Delta_{\mathrm{inter}}^{\dagger}(\mathbf{k})&-\mathcal{H}_{\mathrm{kin},2}^{\mathrm{T}}(-\mathbf{k})\end{array}\right).$
(13)

Here, we suppress the subscript $\alpha=L,R$ when only discussing the bulk
superconductor residing at one side of the junction in this and later
sections.

We use the following two-band models to describe the two topological Fermi
surfaces with nontrivial Chern numbers $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$:

$\displaystyle\mathcal{H}_{\mathrm{kin},1}(\mathbf{k}=\mathbf{K}_{0}+\tilde{\mathbf{k}})$
$\displaystyle=\mathbf{h}_{1}(\tilde{\mathbf{k}})\cdot\bm{\sigma}-\mu_{1},$
(14a)
$\displaystyle\mathcal{H}_{\mathrm{kin},2}(\mathbf{k}=-\mathbf{K}_{0}+\tilde{\mathbf{k}})$
$\displaystyle=\mathbf{h}_{2}(\tilde{\mathbf{k}})\cdot\bm{\sigma}-\mu_{2},$
(14b)

where $\tilde{\mathbf{k}}$ is the momentum relative to the corresponding Weyl
nodes at $\pm\mathbf{K}_{0}$ and $\mathbf{h}_{i}(\tilde{\mathbf{k}})$ denotes
the pseudo “magnetic” field in momentum space which determines the chirality
of Weyl node at $\mathbf{K}_{0}$ ($-\mathbf{K}_{0}$) for $i=1$ ($i=2$). There
are helical eigenstates $\chi_{i,s}(\tilde{\mathbf{k}})$ of
$\mathcal{H}_{\mathrm{kin},i}$ in which $i=1,2$ refers to the valley index and
$s=\pm$ denotes the band index. Here, we use $s=+$ to denote the band at the
Fermi level. The Fermi surfaces have Chern numbers
$\mathcal{C}_{i}={(1/2\pi)\oint_{\mathrm{FS}_{i}}\bm{\Omega}_{i}(\tilde{\mathbf{k}})\cdot
d\mathbf{S}_{\tilde{\mathbf{k}}}}$, where
$\bm{\Omega}_{i}(\tilde{\mathbf{k}})=\bm{\nabla}_{\tilde{\mathbf{k}}}\times\langle\chi_{i,+}(\tilde{\mathbf{k}})|i\bm{\nabla}_{\tilde{\mathbf{k}}}|\chi_{i,+}(\tilde{\mathbf{k}})\rangle$
is the corresponding single particle Berry curvature.

The inter-Fermi surface pairing,
$\Delta_{\mathrm{inter};\sigma,\sigma^{\prime}}(\mathbf{k})\sim\langle
c_{2;-\mathbf{k},\sigma^{\prime}}c_{1;\mathbf{k},\sigma}\rangle$ in the
spin-$\uparrow,\downarrow$ representation, describes the zero center-of-mass
momentum pairing, since FS1 and FS2 are related by parity. In the weak-
coupling regime, superconducting pairing occurs at low energy close to the
Fermi surfaces; hence, we express the pairing in the helical band
representation as
$\Delta_{\mathrm{inter};s,s^{\prime}}^{(b)}(\tilde{\mathbf{k}})\sim\langle\psi_{2;-\mathbf{\tilde{k}},s^{\prime}}\psi_{1;\mathbf{\tilde{k}},s}\rangle$,
where
$\psi_{1;\mathbf{\tilde{k}},s}=\sum_{\sigma=\uparrow,\downarrow}\chi_{1;s,\sigma}^{*}(\tilde{\mathbf{k}})c_{1;\tilde{\mathbf{k}}+\mathbf{K}_{0},\sigma}$
and
$\psi_{2;\mathbf{\tilde{k}},s}=\sum_{\sigma=\uparrow,\downarrow}\chi_{2;s,\sigma}^{*}(\tilde{\mathbf{k}})c_{2;\tilde{\mathbf{k}}-\mathbf{K}_{0},\sigma}$
are annihilation operators of helical band eigenstates of
$\mathcal{H}_{\mathrm{kin},1}$ and $\mathcal{H}_{\mathrm{kin},2}$
respectively. When the Chern numbers of two Fermi surfaces are different, the
inter-Fermi surface pairing inherits the band topology of single-particle
states at $\mathbf{k}$ and $-\mathbf{k}$ in a nontrivial way, and the Cooper
pair scatters in the presence of a magnetic monopole in momentum space given
by the pair monopole charge, $q_{p}=(\mathcal{C}_{1}-\mathcal{C}_{2})/2$.
Consequently, this leads to a topologically obstructed superconducting pairing
order that cannot be well defined over an entire Fermi surface. Instead, the
pairing order is described by monopole harmonics [11, 37] that all belong to
the topological sector characterized by the pair monopole charge $q_{p}$ [1].

Upon projecting to the helical Fermi surfaces in the weak-coupling limit, the
effective pairing in the helical band-diagonal Nambu basis
$(\psi_{1;\tilde{\mathbf{k}},+},\psi^{\dagger}_{2;-\tilde{\mathbf{k}},+})^{T}$
is

$\Delta^{(bp)}_{\mathrm{inter}}(\tilde{\mathbf{k}})=\left(\begin{array}[]{cc}\Delta_{\mathrm{MSC}}^{(q_{p})}(\tilde{\mathbf{k}})&0\\\
0&0\end{array}\right),$ (15)

in which

${\Delta}_{\mathrm{MSC}}^{(q_{p})}(\tilde{\mathbf{k}})=\left\langle\chi_{1,+}(\tilde{\mathbf{k}})\middle|\Delta_{\mathrm{inter}}(\mathbf{k})\middle|\chi_{2,+}^{*}(-\mathbf{\tilde{k}})\right\rangle.$
(16)

In the pseudospin space formed by the helical band eigenbasis, we define an
effective $d$-vector for the monopole pairing order as

$\mathbf{d}^{(bp)}_{\mathrm{inter}}(\tilde{\mathbf{k}})=-\frac{1}{2}\Delta_{\mathrm{MSC}}^{(q_{p})}(\tilde{\mathbf{k}})(1,i,0)^{\mathrm{T}},$
(17)

with
$\Delta_{\mathrm{inter}}^{(bp)}(\tilde{\mathbf{k}})\equiv\mathbf{d}^{(bp)}_{\mathrm{inter}}(\tilde{\mathbf{k}})\cdot\bm{\sigma}(i\sigma_{y})$.
In the following sections, we will use the effective pairing order and the
symmetry principles encoded in the form factor,
$\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})$, to extract the symmetry and
topology of the monopole superconducting order.

To describe a monopole superconductor with pair monopole charge $q_{p}=-1$, we
consider the pairing between Fermi surfaces of Chern numbers
$\mathcal{C}_{1}=-1$ and $\mathcal{C}_{2}=1$. The Hamiltonian in Eq. (13)
takes a specific form below, with
$\mathbf{h}^{(\nu=1)}_{1}(\tilde{\mathbf{k}})=\hbar
v_{F}(\tilde{k}_{x},\tilde{k}_{y},\tilde{k}_{z})$ and
$\mathbf{h}^{(\nu=-1)}_{2}(\tilde{\mathbf{k}})=\hbar
v_{F}(\tilde{k}_{x},\tilde{k}_{y},-\tilde{k}_{z})$ defined near the Weyl nodes
at $\pm\mathbf{K}_{0}=(0,0,\pm K_{0,z})$ respectively,

$\displaystyle\mathcal{H}_{\mathrm{BdG}}^{(q_{p}=-1)}(\mathbf{k})=\left(\begin{array}[]{cc}\mathcal{H}_{\mathrm{kin}}^{(\nu=1)}({\mathbf{k}})&\Delta_{\mathrm{inter}}(\mathbf{k})\\\
\Delta_{\mathrm{inter}}^{\dagger}(\mathbf{k})&-[\mathcal{H}_{\mathrm{kin}}^{(\nu=-1)}(-\mathbf{k})]^{\mathrm{T}}\end{array}\right).$
(18)

Here, the band Hamiltonians are given by

$\displaystyle\mathcal{H}_{\mathrm{kin}}^{(\nu=1)}(\mathbf{k})$
$\displaystyle=\mathbf{h}^{(\nu=1)}_{1}(\tilde{\mathbf{k}})\cdot\bm{\sigma}-\mu,$
(19a) $\displaystyle\mathcal{H}_{\mathrm{kin}}^{(\nu=-1)}(\mathbf{k})$
$\displaystyle=\mathbf{h}^{(\nu=-1)}_{2}(\tilde{\mathbf{k}})\cdot\bm{\sigma}-\mu,$
(19b)

with the superscript $\nu$ denoting the chiralities of the Weyl nodes. $v_{F}$
is the Fermi velocity, and $\mu>0$ is chosen so that we have disjoint Fermi
surfaces which are related by parity via
$\mathcal{H}^{(\nu=+1)}_{1}(\mathbf{k})=\sigma_{z}\mathcal{H}^{(\nu=-1)}_{2}(-\mathbf{k})\sigma_{z}$.
The helical eigenstate of $\mathcal{H}^{{(\nu=1)}}_{\mathrm{kin}}$ in Eq. (19)
describing $\mathrm{FS}_{1}$ takes the form

$\displaystyle\chi_{1,+}(\tilde{\mathbf{k}})$
$\displaystyle=\left(\cos\frac{\theta_{\tilde{k}}}{2}e^{i(\lambda-1)\varphi_{\tilde{k}}/2},\sin\frac{\theta_{\tilde{k}}}{2}e^{i(\lambda+1)\varphi_{\tilde{k}}/2}\right)$
(20)
$\displaystyle=\sqrt{2\pi}\left(\mathcal{Y}_{\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}}),-\mathcal{Y}_{\frac{1}{2};\frac{1}{2},+\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\right)^{T},$

and is related to states on $\mathrm{FS}_{2}$ by
$\chi_{2,+}(\tilde{\mathbf{k}})=\sigma_{z}\chi_{1,+}(-\tilde{\mathbf{k}})$.
Here, $\bm{\Omega}_{\tilde{{k}}}=(\theta_{\tilde{k}},\varphi_{\tilde{k}})$,
with $\theta_{\tilde{k}}$ and $\varphi_{\tilde{k}}$ being the polar and
azimuthal angles of the spherical Fermi surfaces.
$\mathcal{Y}_{q;l,l_{z}}(\bm{\Omega}_{\tilde{k}})$ is the monopole harmonic
function [11, 37] (reviewed in Appendix C), now defined in momentum space. The
monopole harmonic function is labelled by the monopole charge $q$, the
eigenvalue of angular momentum $l$, and the eigenvalue of $z$-component
angular momentum $l_{z}$. Above, $\lambda=\pm 1$ corresponds to two different
choices of gauge. The upper (lower) sign denotes the choice in which the state
$\chi_{1,+}(\tilde{\mathbf{k}})$ is well-defined locally near the north
(south) pole of $\mathrm{FS}_{1}$, while $\chi_{2,+}(\tilde{\mathbf{k}})$ is
well-defined near the south (north) pole of $\mathrm{FS}_{2}$, with the Dirac
strings located at the respective antipodal points of these two Fermi
surfaces. The Chern numbers of $\chi_{1,+}$ on $\mathrm{FS}_{1}$ and
$\chi_{2,+}$ on $\mathrm{FS}_{2}$ are $\mathcal{C}_{1}=-1$ and
$\mathcal{C}_{2}=+1$ respectively. The unitary transformations that
diagonalize the Weyl band Hamiltonians $\mathcal{H}_{\mathrm{kin}}^{(\nu=\pm
1)}$ in Eq. (19) take the form

$U^{(\pm\mathbf{K}_{0})}(\tilde{\mathbf{k}})=\sqrt{2\pi}\left(\begin{array}[]{cc}\mathcal{Y}_{\pm{\frac{1}{2}};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})&\mathcal{Y}_{{\mp\frac{1}{2}};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\\\
\mp\mathcal{Y}_{\pm{\frac{1}{2}};\frac{1}{2},+\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})&\mp\mathcal{Y}_{{\mp\frac{1}{2}};\frac{1}{2},+\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\end{array}\right).$
(21)

As introduced below Eq. (11), the unitary transformations can be written
conveniently in terms of complex four-vectors, where here, the components are
given by

$\displaystyle u_{0}^{(\pm\mathbf{K}_{0})}(\pm\tilde{\mathbf{k}})$
$\displaystyle=\sqrt{\frac{\pi}{2}}\left[\mathcal{Y}_{\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\mp\mathcal{Y}^{*}_{\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\right],$
$\displaystyle u_{1}^{(\pm\mathbf{K}_{0})}(\pm\tilde{\mathbf{k}})$
$\displaystyle=\sqrt{\frac{\pi}{2}}\left[\mp\mathcal{Y}_{\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})-\mathcal{Y}^{*}_{\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\right],$
$\displaystyle u_{2}^{(\pm\mathbf{K}_{0})}(\pm\tilde{\mathbf{k}})$
$\displaystyle=i\sqrt{\frac{\pi}{2}}\left[\pm\mathcal{Y}_{\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})-\mathcal{Y}^{*}_{\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\right],$
$\displaystyle u_{3}^{(\pm\mathbf{K}_{0})}(\pm\tilde{\mathbf{k}})$
$\displaystyle=\sqrt{\frac{\pi}{2}}\left[\mathcal{Y}_{\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\pm\mathcal{Y}^{*}_{\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\right],$
(22)

in which
$\mathcal{Y}^{*}_{q;l,l_{z}}=(-1)^{q+l_{z}}\mathcal{Y}_{-q;l,-l_{z}}$.

For a general form of inter-Fermi-surface pairing
$[\Delta_{\mathrm{inter}}(\mathbf{k})]_{\sigma,\sigma^{\prime}}$, the
effective pairing order takes the form

$\displaystyle\Delta_{\mathrm{MSC}}^{(q_{p}=-1)}=\sqrt{\frac{4\pi}{3}}\mathcal{Y}_{-1;1,1}(\bm{\Omega}_{\tilde{\mathbf{k}}})[\Delta_{\mathrm{inter}}(\mathbf{k})]_{\uparrow,\uparrow}$
(23)
$\displaystyle-\sqrt{\frac{2\pi}{3}}\mathcal{Y}_{-1;1,0}(\bm{\Omega}_{\tilde{\mathbf{k}}})\Big{(}[\Delta_{\mathrm{inter}}(\mathbf{k})]_{\uparrow,\downarrow}-[\Delta_{\mathrm{inter}}(\mathbf{k})]_{\downarrow,\uparrow}\Big{)}.$
$\displaystyle-\sqrt{\frac{4\pi}{3}}\mathcal{Y}_{-1;1,+1}(\bm{\Omega}_{\tilde{\mathbf{k}}})[\Delta_{\mathrm{inter}}(\mathbf{k})]_{\downarrow,\downarrow}.$

Here, the spin pairing channels in addition to the momentum-dependence of
$\Delta_{\mathrm{inter}}(\mathbf{k})$ determine the rotational irreps in the
decomposition of $\Delta_{\mathrm{MSC}}^{(q_{p}=-1)}$ but do not influence the
underlying topology. When the system has rotational symmetry, $l_{z}$ is a
conserved quantity and further constrains the pairing channels. Globally, the
monopole superconducting order transforms under rotation $R_{z}$ according to
its conserved angular momentum $l_{z}=l_{z,\mathrm{glob}}$. However, the
pairing order is a singular representation of $l_{z}$ and thus cannot be well-
defined globally over the Fermi surface under a single gauge. In other words,
due to $\mathrm{U}(1)$ obstruction of the pairing phase, the monopole pairing
order can no longer be described globally by conventional spherical harmonics.
Though, depending on the choice of gauge, $\lambda=\pm 1$, the $q_{p}=-1$
monopole pairing can be described locally with spherical harmonic symmetry,
with $l_{z,\mathrm{loc}}=l_{z,\mathrm{glob}}+q_{p}\lambda$ denoting the local
angular momentum. Here, $l_{z,\mathrm{loc}}$ plays an essential role in the
underlying topology of the system in momentum space. Namely, the shift in
angular momentum reflects the nontrivial pair monopole charge
$|q_{p}|=|l_{z,\mathrm{loc}}-l_{z,\mathrm{glob}}|,$ which, in addition to the
discrepancy between $l_{z,\mathrm{loc}}$ near the gap nodes, is unique to
monopole harmonic symmetry. This shift between global and local angular
momentum is a general feature of monopole superconducting order, which we
shall later exploit in order to distinguish it from other pairing orders with
spherical harmonic symmetry.

As an example, when there is inter-Fermi surface $s$-wave pairing,
$\Delta_{\mathrm{inter}}=\Delta_{0}i\sigma_{y}$, the effective low-energy
pairing in Eq. (23) resides in the lowest possible angular momentum $l=1$
channel as follows,

$\displaystyle{\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}}(\tilde{\mathbf{k}})$
$\displaystyle=-\Delta_{0}\sqrt{\frac{8\pi}{3}}\mathcal{Y}_{-1;1,0}(\theta_{\tilde{k}},\varphi_{\tilde{k}})$
(24)
$\displaystyle=-4\pi\Delta_{0}\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},+\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})$
$\displaystyle=-\Delta_{0}\sin\theta_{\tilde{k}}e^{-i\lambda\varphi_{\tilde{k}}}.$

Here, $\mathcal{H}_{\mathrm{BdG}}^{(q_{p}=-1)}(\mathbf{k})$ is invariant under
rotation $R_{z}$, and the low-energy pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1)}$, if it does not break the symmetry
spontaneously, is also invariant under $R_{z}$ and characterized by good
quantum number $l_{z}=l_{z,\mathrm{glob}}=0$ that labels the trivial irrep.
However, for given choice of gauge $\lambda=\pm 1$, the
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ pairing order locally has well-
defined $p_{x}-\lambda ip_{y}$ symmetry with angular momentum
$l_{z,\mathrm{loc}}=-\lambda$ near the gap node at the north ($\lambda=+1$) or
south ($\lambda=-1$) pole of $\mathrm{FS}_{1}$. As $\tilde{k}_{z}$ is varied
between the gap nodes, the local symmetry changes from $p_{x}-ip_{y}$ to
$p_{x}+ip_{y}$, as evident in the group velocity of the surface modes (see
Appendix E). The shift between the global and local angular momenta is
indicative of the pair monopole charge, $q_{p}=-1$. Moreover, this shift is
generally true for any form of $\Delta_{\mathrm{inter}}$ and therefore for a
general $\Delta_{\mathrm{MSC}}^{(q_{p}=-1)}$ in Eq. (23).

We now derive the form factor in the first order Josephson coupling between a
spin singlet superconductor with a spin-independent band Hamiltonian at the
left side and the $q_{p}=-1$ monopole superconductor at the right side, in
which there is only bare spin-independent tunneling at the barrier. In the
band-projected basis, the monopole pairing order is described by an effective
$d$-vector shown in Eq. (17). The form factor hence takes the form

$\displaystyle\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}(\mathbf{k},\mathbf{k}^{\prime})=$
(25) $\displaystyle
2{d_{L,0}(\mathbf{k}^{\prime})}{T^{(b)}_{0}(-\mathbf{k},-\mathbf{k}^{\prime})}[{\mathbf{d}^{(bp)*}_{R}(\mathbf{k})}\cdot{\mathbf{T}^{(b)}(\mathbf{k},\mathbf{k}^{\prime})}]$
$\displaystyle-2{d_{L,0}(\mathbf{k}^{\prime})}{T_{0}^{(b)}(\mathbf{k},\mathbf{k}^{\prime})}[{\mathbf{d}^{(bp)*}_{R}(\mathbf{k})}\cdot{\mathbf{T}^{(b)}(-\mathbf{k},-\mathbf{k}^{\prime})}]$
$\displaystyle-2i{d_{L,0}(\mathbf{k}^{\prime})}{\mathbf{d}^{(bp)*}_{R}(\mathbf{k})}\cdot[{\mathbf{T}^{(b)}(\mathbf{k},\mathbf{k}^{\prime})}\times{\mathbf{T}^{(b)}(-\mathbf{k},-\mathbf{k}^{\prime})}],$

which is analogous to that of Eq. (9), only now the effective $d$-vector and
tunneling matrices are evaluated in the projected helical band basis. The
tunneling amplitudes that describe the tunneling between the states at the
left, in the spin-$\uparrow$, $\downarrow$ representation, and the states at
the right, in the helical band representation, are given by

$\displaystyle T^{(b)}_{0}(\pm\mathbf{k},\pm\mathbf{k}^{\prime})$
$\displaystyle=T_{0}(\pm\mathbf{k},\pm\mathbf{k}^{\prime})[u_{R,0}^{(\pm\mathbf{K}_{0})}(\pm\tilde{\mathbf{k}})]^{*}$
(26a) $\displaystyle\mathbf{T}^{(b)}(\pm\mathbf{k},\pm\mathbf{k}^{\prime})$
$\displaystyle=T_{0}(\pm\mathbf{k},\pm\mathbf{k}^{\prime})[\mathbf{u}_{R}^{(\pm\mathbf{K}_{0})}(\pm\tilde{\mathbf{k}})]^{*}.$
(26b)

Above,
$(u_{R,0}^{(\pm\mathbf{K}_{0})}(\tilde{\mathbf{k}}),\mathbf{u}_{R}^{(\pm\mathbf{K}_{0})}(\tilde{\mathbf{k}}))$
are the complex four-vectors defined in Eq. (22).

Consequently, the form factor describing the first order Josephson coupling
between a $q_{p}=-1$ monopole superconductor and a singlet superconductor
simplifies to

$\displaystyle\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}^{(q_{p}=-1)}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-4\pi{\Delta_{L,0}(\mathbf{k}^{\prime})}[{\Delta_{R,\mathrm{MSC}}^{{(q_{p}=-1)}}(\tilde{\mathbf{k}})}]^{*}$
(27)
$\displaystyle\times{T_{0}(-\mathbf{k},-\mathbf{k}^{\prime})}{T_{0}(\mathbf{k},\mathbf{k}^{\prime})}\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}}),$

in which $\Delta_{L,0}(\mathbf{k}^{\prime})=d_{L,0}(\mathbf{k}^{\prime})$,
$\Delta_{R,\mathrm{MSC}}^{{(q_{p}=-1)}}(\tilde{\mathbf{k}})$ is the band-
projected effective monopole harmonic pairing order in Eq. (23) characterized
by the pair monopole charge $q_{p}=-1$, and
$\tilde{\mathbf{k}}=\mathbf{k}-\mathbf{K}_{0}$ is the wave vector with respect
to the Weyl node. Despite the form factor expressed in Eq. (27) including
multiple monopole harmonic functions,
$\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}$ is always non-singular, as it is
related to the Josephson current, a physical observable. By employing the
addition theorem of monopole harmonic functions [37] (see also Appendix C),
which states that the product of two monopole harmonic functions with charges
$q_{1}$ and $q_{2}$ can be expressed as a sum of monopole harmonic functions
with monopole charge $q_{3}=q_{1}+q_{2}$, we find that the form factor in Eq.
(27) has zero total monopole charge. In other words, it can be represented
using standard spherical harmonic functions without singularities over the
Fermi surface.

Similarly, we derive the form factor for a junction between the same
$q_{p}=-1$ monopole superconductor and a standard spin triplet superconductor
as follows, which is non-singular for the same reason as that of
$\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}^{(q_{p}=-1)}$,

$\displaystyle\mathfrak{F}_{\mathrm{trip}-\mathrm{MSC}}^{(q_{p}=-1)}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-2\pi{T_{0}(-\mathbf{k},-\mathbf{k}^{\prime})}{T_{0}(\mathbf{k},\mathbf{k}^{\prime})}\times$
$\displaystyle{\Delta_{R,\mathrm{MSC}}^{{(q_{p}=-1)}^{*}}(\tilde{\mathbf{k}})}[(d_{L,x}(\mathbf{k}^{\prime})-id_{L,y}(\mathbf{k}^{\prime}))\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}^{2}(\bm{\Omega}_{\tilde{k}})$
$\displaystyle+(d_{L,x}(\mathbf{k}^{\prime})+id_{L,y}(\mathbf{k}^{\prime}))\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},\frac{1}{2}}^{2}(\bm{\Omega}_{\tilde{k}})].$
(28)

It depends on, $\mathbf{d}_{L}(\mathbf{k}^{\prime})$, the $d$-vector which
encodes the spin structure of the triplet pairing order at the left. In
contrast, the form factor
$\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}^{(q_{p}=-1)}$ in Eq. (27) takes a
simpler form, as the symmetry of the singlet pairing order at the left is
primarily determined by its symmetry in momentum space. For the microscopic
models studied later in Sec. IV and VI, to identify the momentum space
topology of the monopole pairing order, we will focus on junctions with a spin
singlet superconductor for simplicity.

###  III.2 Symmetry and topological principles to probe monopole
superconducting order

Symmetry principles have long been crucial in the design of Josephson
experiments to reveal the pairing symmetry of unconventional superconductors
[36, 16, 20, 21]. We now propose symmetry and topological principles to guide
the design of a set of Josephson junctions which together can identify
monopole superconducting order. This approach complements the discussion in
Sec. III.1, which was mainly based on microscopic derivations of Josephson
current.

Figure 1:  Josephson junction designs to probe monopole superconducting order,
shown schematically in a hybrid picture of momentum and real space. The
junction barrier (gray, striped) is shown in real space, separating two
superconductors at the left (L) and right (R). Band Fermi surfaces are shown
in momentum space at either side of the junction. (a) Junction design to probe
the global angular momentum, consisting of a junction oriented along the
rotational axis ($z$ axis) between a monopole superconductor (cyan) and other
known superconductor (orange). (b) Junction design to probe the local angular
momentum of the pairing order in momentum space. The junction is between two
identical monopole superconductors and is oriented perpendicular to the
rotational $z$-axis along which the gap nodes lie. For conserved momentum
$k_{z}$, it can be treated as an effective 2D junction (yellow).

Monopole pairing order with pair monopole charge $q_{p}$ generally takes the
form of a superposition of monopole harmonic functions
$\mathcal{Y}_{q_{p},l,l_{\hat{\mathbf{n}},\mathrm{glob}}}$ for different
integer values of $l\geq|q_{p}|$ and $|l_{\hat{\mathbf{n}}}|\leq l$. For
$q_{p}\neq 0$, the monopole harmonic functions exhibit nodes pinned to the
rotation quantization axis, $\hat{\mathbf{n}}$, which we often take as
$\hat{z}$ without loss of generality. Furthermore, as discussed in Sec. III.1,
local angular momenta differ from global angular momentum by the pair monopole
charge,

$l_{\hat{\mathbf{n}},\mathrm{loc}}=l_{\hat{\mathbf{n}},\mathrm{glob}}+\lambda
q_{p}.$ (29)

Here, $\lambda=\pm 1$ corresponds to different choices of gauge (see Appendix
C). To identify this exotic pairing order, we design a set of two Josephson
junctions which respectively probe the global and local angular momentum of
the pairing order. The nonzero shift between $l_{z,\mathrm{glob}}$ and
$l_{z,\mathrm{loc}}$ indicates the pair monopole charge.

(I) Globally, the monopole superconducting order remains a representation of
the rotation group but resides in a topological sector characterized by its
pair monopole charge $q_{p}$. For two superconductors that transform
identically under a given symmetry operation, which is also a symmetry of the
junction, the Josephson energy phase relationship is invariant under this
symmetry transformation. Hence, if rotation $R_{z}$ with respect to the $z$
axis is a symmetry of the junction, then the first order Josephson coupling is
allowed only if the superconducting orders at both sides of the junction
transform in the same way under $R_{z}$, i.e. they both form an irrep of
$R_{z}$ labeled by the same $l_{z}$. Monopole superconducting order given by,
for example, $\mathcal{Y}_{q_{p};l,l_{z,\mathrm{glob}}}$ transforms under
rotation $R_{z}$ according to its globally conserved angular momentum,
$l_{z,\mathrm{glob}}$. As such, to extract $l_{z,\mathrm{glob}}$, we consider
a junction which is oriented along a rotational $z$ axis, as shown in Fig. 1
(a). By placing a monopole superconductor in junction with another
superconductor with known pairing symmetry, one can extract the global angular
momentum $l_{z,\mathrm{glob}}$ based on the existence of first order
$2\pi$-periodic Josephson current.

(II) Locally in momentum space, the symmetry of monopole superconductor
pairing order can be described by regular spherical harmonic symmetry. In the
local representation, the pairing phase possesses nontrivial winding (see
Appendix C), corresponding to chiral pairing order characterized by angular
momentum $l_{z,\mathrm{loc}}$. Between the gap nodes pinned to the rotational
axis, the local chirality changes by twice the pair monopole charge, as shown
in the model in Sec. III.1. To probe $l_{z,\mathrm{loc}}$, we consider the
design shown in Fig. 1 (b), in which two identical monopole superconductors
form a Josephson junction oriented normal to the $z$ axis. For given conserved
momentum $k_{z}$ along the rotational axis, the junction can be regarded
effectively as a junction between two identical two-dimensional
superconductors whose pairing orders are characterized by local angular
momentum $l_{z,\mathrm{loc}}$. Consequently, the total Josephson current of
the proposed junction is a superposition of contributions from these planar
junctions, which can give signatures, for example the periodicity and phase
shift of the current phase relation, that can indicate the local pairing
angular momentum $l_{z,\mathrm{loc}}$.

In the following two sections, we analyze a set of Josephson junctions guided
by the above design principles to probe the pairing symmetry of a prototypical
monopole superconductor arising in a magnetic doped Weyl semimetal with inter-
Fermi surface pairing in Eq. (18). In Sec. IV and Sec. V, we propose junctions
to probe the global angular momentum $l_{z,\mathrm{glob}}$ and local angular
momentum $l_{z,\mathrm{loc}}$ respectively. Together, these junctions can be
used to extract the monopole charge,
$|q_{p}|=|l_{z,\mathrm{glob}}-l_{z,\mathrm{loc}}|$, and provide a foundation
to guide efforts to experimentally uncover the monopole superconducting order.

##  IV Josephson junctions to probe $l_{z,\mathrm{glob}}$ of the monopole
superconducting order

Following the symmetry principles presented in Sec. III.2, we design Josephson
junctions to extract the global angular momentum of monopole superconducting
order, $l_{z,\mathrm{glob}}$. We first examine a Josephson junction oriented
along the $z$ direction between a monopole superconductor with pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ and a conventional $s$-wave
superconductor, demonstrating the nonvanishing first order Josephson current.
After, in Sec. IV.2, we study the Josephson coupling for the monopole pairing
order $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}$ in the same topological
sector but with different global angular momentum, and in Sec. IV.3, monopole
pairing order $\Delta_{\mathrm{MSC}}^{(q_{p}=3,l_{z}=0)}$ in a different
topological sector of $q_{p}=3$.

###  IV.1 First order Josephson current between $s$-wave superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor

####  IV.1.1 Method 1: Linear Response Theory

We study Josephson effects between an $s$-wave superconductor and monopole
superconductor with pairing order $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$
in a junction along the $z$ direction, shown schematically in Fig. 2 (a).
Following the symmetry principles, as discussed in Sec. III.2, the first order
Josephson coupling in this system is allowed due to the pairing order at two
sides of the junction being characterized by the same global angular momentum
$l_{z}=0$. To show this, we consider the continuum model of the monopole
superconductor introduced in Eq. (18), for which the effective pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ in Eq. (24) is invariant under
rotation $R_{z}$, having global angular momentum $l_{z}=0$. Following the
derivation based on linear response theory in Sec. III.1, we demonstrate
nonzero first order Josephson coupling and the resulting critical current.

For a junction between an $s$-wave superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor, the form
factor in Eq. (27) can be simplified to

$\displaystyle\mathfrak{F}_{s-\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=\frac{8\pi}{3}\Delta_{L,0}\Delta_{R,0}|Y_{1,1}(\bm{\Omega}_{\tilde{k}})|^{2}T_{0}^{2}\delta_{\mathbf{k}_{\parallel},\mathbf{k}^{\prime}_{\parallel}},$
(30)

where $\Delta_{L,0}$ is the pairing amplitude of the $s$-wave superconductor,
and $\Delta_{R,0}$ is the amplitude of the inter-Fermi surface pairing in the
monopole superconductor. We model the tunneling amplitude through the junction
barrier as a spin-independent delta function that preserves in-plane momenta
with
$T_{0}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}\delta_{\mathbf{k}_{\parallel},\mathbf{k}^{\prime}_{\parallel}}$
and $\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=0$. Here, $T_{0}$ is a
constant, and
$\mathbf{k}_{\parallel}=\mathbf{k}_{\parallel}^{\prime}=\tilde{\mathbf{k}}_{\parallel}$
is the wave vector in the $k_{x}k_{y}$ plane, which is conserved at the
barrier. Already, the nonvanishing first order Josephson coupling can be seen
from the form factor $\mathfrak{F}_{s-\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$,
which encodes the symmetry of the junction as well as the pairing orders and
itself is invariant under rotation $R_{z}$ and parity.

With the form factor and tunneling matrix above, the Josephson current in Eq.
(4) evaluates to

$\displaystyle
I_{J}(\phi)=\frac{4e}{h}AT_{0}^{2}\Delta_{L,0}\Delta_{R,0}\sin\phi\int\mathrm{d}{k}_{\parallel}\
{k}_{\parallel}J({k}_{\parallel}),$ (31)

in which

$\displaystyle
J({k}_{\parallel})=k_{\parallel}^{2}\frac{\ell_{L}\ell_{R}}{(2\pi)^{2}}\int$
$\displaystyle\frac{\mathrm{d}\tilde{k}_{z}\mathrm{d}k^{\prime}_{z}}{{k}_{\parallel}^{2}+\tilde{k}_{z}^{2}}w(E_{L,\mathbf{k}^{\prime}},E_{R,\tilde{\mathbf{k}}}^{(bp)};\beta)$
(32)

describes the contribution to the Josephson current for conserved transverse
momentum, ${k}_{\parallel}$. Above, $A$ is the surface area of the interface,
and $\ell_{L}$ and $\ell_{R}$ are the system lengths along $z$ direction at
the left and right sides respectively. We consider the Josephson current at
zero temperature, for which
$w(E_{L,\mathbf{k}^{\prime}},E_{R,\tilde{\mathbf{k}}}^{(bp)};\beta)=[2E_{L,\mathbf{k}^{\prime}}E_{R,\tilde{\mathbf{k}}}^{(bp)}(E_{L,\mathbf{k}^{\prime}}^{2}-E_{R,\tilde{\mathbf{k}}}^{(bp)2})]^{-1}$
and is positive-definite. The dispersions of the BdG quasiparticles for the
$s$-wave superconductor and monopole superconductor are given by
$E_{L,\mathbf{k}^{\prime}}=\sqrt{\xi_{L,\mathbf{k}^{\prime}}^{2}+|\Delta_{L,0}|^{2}}$
and
$E_{R,\tilde{\mathbf{k}}}^{(bp)}=\sqrt{\xi_{R,\tilde{\mathbf{k}}}^{(bp)2}+\left|\Delta_{R,0}\right|^{2}{\tilde{k}^{2}_{\parallel}}/{\tilde{k}^{2}}}$
respectively. Here,
$\xi_{L,\mathbf{k}^{\prime}}=\hbar^{2}({k^{\prime}_{z}}^{2}+{k^{\prime}}_{\parallel}^{2})/2m-\mu_{L}$
and $\xi^{(bp)}_{R,\tilde{\mathbf{k}}}=\hbar
v_{F}\sqrt{\tilde{k}_{z}^{2}+\tilde{k}_{\parallel}^{2}}-\mu_{R}$ are the
respective band dispersions. The Fermi surfaces of
$\xi_{L,\mathbf{k}^{\prime}}$ and $\xi_{R,\tilde{\mathbf{k}}}^{(bp)}$ are
shown projected to the $k_{x}k_{y}$ plane in Fig. 2 (b).

Figure 2:  (a) Schematic of the Josephson junction oriented along $z$
direction between an $s$-wave superconductor and a
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor. (b) Bulk
Fermi surface projections on the $k_{x}k_{y}$ plane for $s$-wave
superconductor (hatched) and monopole superconductor (solid) at two sides of
the junction in (a). (c) Relative magnitude of the contributions to the first
order Josephson current phase relation for different transverse momentum,
$k_{\parallel}$, corresponding to Eq. (32). (d) First order Josephson current
between the monopole superconductor and $s$-wave superconductor at zero
temperature and zero bias. Here, $v_{F}=0.5$, $\mu_{R}=0.1$,
$\tilde{k}_{F}=0.2$, and $\mu_{L}=0.8$.
$\Delta_{L,0}=\Delta_{R,0}=\Delta_{0}=0.005$, $I_{c}=0.62\Delta_{0}e/\hbar$,
and $N_{k_{\parallel}}$ is the number of transverse momentum cuts in the
calculation. $T_{0}$, $A$, $\ell_{L}$, $\ell_{R}$, and $m$ are set to unity.

In Fig. 2 (c), we show $J(k_{\parallel})$, the contributions to the first
order Josephson current, for different conserved transverse momentum
$k_{\parallel}$. States near the equator of the Fermi surface of the doped
Weyl semimetal, with $k_{\parallel}\approx\tilde{k}_{F}$, contribute the most
to the Josephson tunneling, where the magnitude of
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ is at its maximum. As
$k_{\parallel}$ approaches zero, the contributions to the Josephson current
are suppressed as the magnitude of
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ reaches its minimum near its gap
nodes along the $k_{z}$ axis. For transverse momentum away from the Weyl Fermi
surfaces, $k_{\parallel}>\tilde{k}_{F}+{\Delta_{R,0}}/{v_{F}}$, its
contribution to the Josephson current vanishes due to vanishing low energy
density of states in $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole
superconductor.

Furthermore, in Fig. 2 (d), we show the total first order Josephson current
phase relation given in Eq. (31) numerically. It exhibits $2\pi$ periodicity
and nonvanishing critical current $I_{c}=0.62\Delta_{0}e/\hbar$. In contrast,
for a Josephson junction in the same geometry between an $s$-wave and a chiral
$p_{x}\pm ip_{y}$-wave superconductor, the first order Josephson current
vanishes when there is no spin-orbit coupling in the bulk nor at the
interface. Hence, though locally $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$
exhibits chiral $p_{x}\pm ip_{y}$ pairing symmetry, the Josephson current in
the above junction serves as an effective probe of the global angular momentum
of the pairing order $l_{z,\mathrm{glob}}=0$. We discuss later the case when
there exist spin-orbit interactions in Sec. VI.

####  IV.1.2 Method 2: Numerical results from tight-binding model

We further show the nonzero Josephson coupling in the same Josephson junction
between an $s$-wave superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor, as shown
in Fig. 2 (a), using a tight-binding model on a cubic lattice (see Appendix D
and E for details). As $k_{x}$ and $k_{y}$ are assumed to be conserved across
the junction, we take periodic boundary conditions in the $x$ and $y$
directions. We consider $2N_{z}$ sites along $z$ direction and take open
boundary conditions at the two ends at $n_{z}=1$ and $2N_{z}$. The junction
barrier is located between $n_{z}=N_{z}$ and $n_{z}=N_{z}+1$, in which only
spin-independent nearest neighbor hopping within the barrier is considered. At
given conserved momenta $k_{x}$ and $k_{y}$, we model the barrier along $z$
direction as

$\displaystyle
H_{\mathrm{link};k_{x},k_{y}}=t_{0}\sum_{\sigma=\uparrow,\downarrow}(c^{\dagger}_{k_{x},k_{y};\sigma,N_{z}+1}c_{k_{x},k_{y};\sigma,N_{z}}+\mathrm{h.c.})$
(33)

in which $t_{0}$ is the spin-independent hopping amplitude, and
$c_{k_{x},k_{y};\sigma,n_{z}}$ annihilates a particle at site $n_{z}$ with
momentum $k_{x}$ and $k_{y}$ and spin $\sigma$. If we increase the thickness
of the barrier to a few sites, there is no qualitative change in Josephson
current.

Figure 3: Numerical results for a tight-binding model of the Josephson
junction along $z$ direction between an $s$-wave superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor. (a)
Maximal cross section of bulk Fermi surface $\mathrm{FS}_{+}$ of the monopole
superconductor at $k_{z}=K_{0,z}$ (solid) and that of the Fermi surface
$\mathrm{FS}_{0}$ of the $s$-wave superconductor at $k_{z}=0$ (dashed). Here,
$\mu_{L}/t_{0}=3.88$ and $\mu_{R}/t_{0}=0.05$. (b) The Josephson energy phase
relation contributed by different momenta in (a),
$k_{\parallel,n}=(n-1)\tilde{k}_{F}/4$, $n=1,\cdots,5$, and
$\tilde{k}_{F}=0.2$. Here, the relative magnitude is given by
$\Delta^{(n)}_{0}=\Delta_{0}k_{\parallel,n}/\tilde{k}_{F}$, with
$\Delta_{0}/t=0.04$. Lowest energy states for a given radial cut correspond to
the cuts shown in (a) and are shown in color, while the gray bands represent
higher-energy contributions. (c) The total Josephson current phase relation
(red, solid) compared to that of a $z$ direction junction between a
$p_{x}-ip_{y}$ superconductor and an $s$-wave SC (black, dashed). (d) The
Fraunhofer pattern of a corner junction (inset) in the $xy$-plane of a
monopole superconductor (red, solid), compared to that of a $p_{x}+ip_{y}$
superconductor (black, dashed). Above, $I_{0}={e\Delta_{0}}/{\hbar}$, and the
system size is $2N_{z}=300$.

In Fig. 3, we present our numerical results of the Josephson current for the
junction between an $s$-wave superconductor and a
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor. Fermi
surfaces of the bulk tight-binding Hamiltonians at the two sides of the
junction are shown in Fig. 3 (a). Although the lattice has four-fold rotation
symmetry with respect to $z$ axis, the Fermi surfaces are nearly spherical in
the $k_{x}k_{y}$ plane. In Fig. 3 (b), contributions to the Josephson energy
phase relation are calculated for representative momenta
${\mathbf{k}}_{\parallel}$ at different radii
$k_{\parallel}=\sqrt{{{k}_{x}^{2}+{k}_{y}^{2}}}$. At a given value of
$k_{\parallel}$, different azimuthal angles, $\varphi_{\tilde{k}}$, contribute
nearly the same energy phase relation. Hence, the contributions from states at
two representative azimuthal angles, $\varphi_{\tilde{k}}=0$ and $\pi/4$,
shown in the inset of Fig. 3 (a), overlap with each other in Fig. 3 (b). At
$k_{\parallel}=k_{\parallel,1}=0$, the states close to north and south poles
of the Fermi surfaces contribute to the Josephson coupling. However, the
magnitude of the pairing gap of the monopole superconductor,
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$, vanishes at the poles of the
Weyl Fermi surfaces. As such, states at $k_{\parallel,1}$ correspond to
vanishingly small Josephson energy and therefore flat energy phase relation.
As $k_{\parallel}$ increases to nonzero values of $k_{\parallel,i}$ with
$i\geq 2$, the corresponding contribution to Josephson energy shows
significant $2\pi$ periodicity. Furthermore, the overall magnitude of
Josephson energy increases as the magnitude of pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ at $k_{\parallel,i}$ increases
until reaching its maximum at $k_{\parallel,5}\approx\tilde{k}_{F}$.

In addition to the net $2\pi$-periodic Josephson energy phase relation, there
are contributions from second order Josephson tunneling, where multiple
scatterings from the junction interface lead to $\pi$-periodic contributions
to the energy phase relation. Such higher order contributions are most
significant for states near the gap nodes of the
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole pairing order, for
example, at $k_{\parallel}=k_{\parallel,2}$ in Fig. 3 (b). We provide further
detailed explanation on this in Appendix B.

We calculate the total Josephson current at zero temperature,
$\left.I_{J}(\phi)\right|_{T=0}=(-{2e}/{\hbar})\sum_{\mathbf{k}_{\parallel}}{\partial
E_{J}(\mathbf{k}_{\parallel},\phi)}/{\partial\phi}.$ Here,
$E_{J}(\mathbf{k}_{\parallel},\phi)$ is the Josephson energy phase relation
for occupied states with transverse momentum $\mathbf{k}_{\parallel}$. In Fig.
3 (c), we show our numerical result of Josephson current with the summation
over transverse momentum $k_{\parallel}\leq\tilde{k}_{F}$ at representative
values, $k_{\parallel,i=1,\cdots,5}$, and taking account of the measure in
polar coordinates. For $k_{\parallel}>\tilde{k}_{F}+{\Delta_{R,0}}/{v_{F}}$,
there are no in-gap Andreev states contributing to the Josephson current as
discussed in Sec. IV.1.1. Moreover, comparing with Fig. 2 (c), the total
Josephson current phase relation in Fig. 3 (c) is modulated by $\pi$-periodic
contributions from higher order tunneling processes. However, the total
Josephson current, obtained by superimposing contributions from different
$k_{\parallel}$, is nonetheless $2\pi$-periodic, as described below.

The first order $2\pi$-periodic Josephson current between an $s$-wave
superconductor and $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole
superconductor is nonvanishing in this tight-binding model for the same
symmetry reasons discussed for the continuum model in Sec. IV.1.1. The tight-
binding model preserves the two parity-related Fermi surfaces which have
opposite Chern numbers. Although the Fermi surfaces are not perfectly
spherical, the topology of the effective pairing remains invariant and is
still classified by the pair monopole charge, $q_{p}=-1$. Additionally, the
global angular momentum of the pairing order is constrained to the $l_{z}=0$
channel, as both the kinetic and the proximitized pairing parts of the tight-
binding Hamiltonian are invariant under lattice rotations $R_{z}$. Therefore,
the effective pairing order likewise transforms according to conserved
$l_{z,\mathrm{glob}}=0$. However, compared to the effective pairing derived in
the continuum model in Eq. (24), there are small but nonzero contributions to
the pairing order from monopole harmonics in the $l=2$ partial wave channel as
a result of the lattice symmetry. Nonetheless, these contributions from higher
partial wave channels do not significantly affect the nonzero first order
Josephson coupling. Hence, because the tight binding model preserves the
symmetry and topology of the system, the essential features of the Josephson
current, for example its $2\pi$-periodicity, agree with the previous result
based on continuum model.

We additionally examine a corner junction between a
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor and
$s$-wave superconductor. The corner junction is symmetrically aligned with two
faces along the $x$ and $y$ directions, as shown schematically in the inset in
Fig. 3 (d). Because this geometry is sensitive to pairing phase difference
along the two orthogonal crystal axes [20], it can be used as an additional
probe of global angular momentum $l_{z,\mathrm{glob}}$. Based on the same
tight-binding model of $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole
superconductor, we study the resulting critical current in the presence of an
external magnetic field along the $z$ direction. As the external flux
$\Phi_{B}$ through the barrier of the corner junction varies, we obtain the
Fraunhofer pattern of the corner junction, as shown in Fig. 3 (e). The
Fraunhofer pattern is qualitatively the same as that of an SIS junction
between two $s$-wave superconductors, with nodes at integer flux quantum,
$\Phi_{0}=h/2e$. This is in contrast to, for example, a corner junction
between an $s$-wave and chiral $p$-wave superconductor, which gives rise to an
interference pattern that is asymmetric with respect to zero flux.

In alignment with the findings presented in the preceding section, in which
the Josephson current was analyzed through linear response theory, the results
from the tight-binding model further reinforce that, globally, the
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconducting pairing
order transforms according to $l_{z,\mathrm{glob}}=0$. Consequently, the
monopole superconductor has nonzero Josephson coupling with an $s$-wave
superconductor in a junction along the $z$ axis, leading to an overall
nonvanishing $2\pi$-periodic first order Josephson current. We next discuss
similar junction designs to probe monopole superconductors with nonzero global
angular momentum $z$ component, or when the pairing order is in a different
topological sector.

###  IV.2 Probing monopole superconducting orders with nonzero global angular
momenta $l_{z,\mathrm{glob}}$

We now demonstrate that the $z$ direction Josephson junction design introduced
in Sec. III.2 can be generalized to probe monopole pairing orders with nonzero
$z$-component global angular momentum. As an example, we consider the monopole
pairing order in the same topological sector as that in Sec. IV.1 given by the
pair monopole charge $q_{p}=-1$, but in a different angular momentum channel
with $l=l_{z}=2$. Following the symmetry argument, the first order Josephson
coupling is allowed only when coupled to another superconductor with pairing
order in the $l_{z}=2$ channel, for example, a $d_{x^{2}-y^{2}}+id_{xy}$ wave
superconductor, in an SIS junction with geometry preserving the rotational
symmetry about $z$ axis. To demonstrate the symmetry principle, here, we focus
on deriving the form factor
$\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}^{(q_{p}=-1)}$ in this system and
show that it leads to nonzero first order Josephson coupling. In Appendix F,
we also show nonvanishing first order Josephson current obtained numerically
from a tight-binding model of this system.

We first introduce the model of a monopole superconductor for which the
pairing order is in the topological sector defined by pair monopole charge
$q_{p}=-1$ yet now transforms under rotation $R_{z}$ according to global
angular momentum $l_{z}=2$. For simplicity, we start from the same band
Hamiltonian of the doped magnetic Weyl semimetal, as shown in Eq. (19). We now
consider the case of inter-Fermi surface $p_{x}+ip_{y}$ pairing in the
${|\uparrow\uparrow\rangle}$ channel, given by
$\Delta_{\mathrm{inter}}(\mathbf{k})=(\Delta_{0}/2\tilde{k}_{F})(k_{x}+ik_{y})(\sigma_{0}+\sigma_{z})$.
The band-projected superconducting order in Eq. (23) describing the low energy
pairing between helical Fermi surfaces enclosing the Weyl nodes at
$\pm\mathbf{K}_{0}$ is now given by

$\displaystyle\Delta_{\mathrm{MSC}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(q_{p}=-1,l_{z}=2)}}(\tilde{\mathbf{k}})$
$\displaystyle=\Delta_{0}\sqrt{\frac{4\pi}{5}}\mathcal{Y}_{-1;2,2}(\theta_{\tilde{k}},\varphi_{\tilde{k}})$
(34)
$\displaystyle=\Delta_{0}\cos^{2}\frac{\theta_{\tilde{k}}}{2}\sin\theta_{\tilde{k}}e^{i(2-\lambda)\varphi_{\tilde{k}}},$

where, in consistency with prior notation, $\lambda=\pm 1$ corresponds to a
gauge in which the Dirac string is located at
$\theta_{\tilde{k}}=\pi(1+\lambda)/2$. As the Chern numbers of the Fermi
surfaces that participate in the inter-Fermi surface pairing are unchanged,
the effective pairing order remains in the same topological sector
characterized by the pair monopole charge, $q_{p}=-1$. However, the pairing
order now transforms according to its global angular momentum,
$l_{z}=l_{z,\mathrm{glob}}=2$. The nonzero global angular momentum can be
attributed to the inter-Fermi surface
$(p_{x}+ip_{y}){|\uparrow\uparrow\rangle}$ pairing, which transforms according
to its total angular momentum $j_{z}=2$ in the presence of spin-orbit
coupling. When coupled to the spin-texture of the Weyl semimetal Fermi
surface, this induces the global $l_{z}=2$ winding of the monopole
superconducting order in momentum space.

Figure 4: Junction along $z$ direction between a $d_{x^{2}-y^{2}}+id_{xy}$
superconductor (left) and $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}$
monopole superconductor (right).

We now consider a Josephson junction along the $z$ direction between a
$d_{x^{2}-y^{2}}+id_{xy}$ superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}$ monopole superconductor, as shown
in Fig. 4. We employ the form factor for a junction between a $q_{p}=-1$
monopole superconductor and spin singlet superconductor,
$\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}^{(q_{p}=-1)}$, as shown in Eq.
(27). Only now,
$\Delta_{L,0}(\mathbf{k}^{\prime})=(\Delta_{L,0}/k_{F}^{2})(k^{\prime}_{x}+ik^{\prime}_{y})^{2}$
describes the $d_{x^{2}-y^{2}}+id_{xy}$ pairing order and
$\Delta_{R,\mathrm{MSC}}^{(q_{p}=-1)}=\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}(\tilde{\mathbf{k}})$
is the monopole superconducting order in Eq. (34). For simplicity, we consider
spin-indpendent tunneling that trivially respects the junction geometry, given
by constant $T_{0}(\mathbf{k},\mathbf{k})=T_{0}$ and
$\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=0$. The form factor for the
junction between a $d_{x^{2}-y^{2}}+id_{xy}$ superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}$ monopole superconductor
simplifies to

$\displaystyle\mathfrak{F}_{d+id,\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}(\mathbf{k},\mathbf{k}^{\prime})=\Delta_{L,0}\Delta_{R,0}T_{0}^{2}\frac{{k^{\prime}}^{2}}{k_{F}^{2}}\cos^{2}\frac{\theta_{\tilde{k}}}{2}\sin^{2}\theta_{\tilde{k}},$
(35)

with $\Delta_{R,0}$ being the amplitude of the inter-Fermi surface pairing.
Because both the $d_{x^{2}-y^{2}}+id_{xy}$ superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=2)}$ monopole superconductor transform
under rotation $R_{z}$ according to their global angular momentum $l_{z}=2$,
the form factor itself is invariant under $R_{z}$. Hence, per the same
symmetry arguments in Sec. IV.1.1, this leads to nonvanishing first order
Josephson coupling.

###  IV.3 Generalization to monopole pairing orders with higher monopole
charges

We further demonstrate the probe of global angular momentum for monopole
superconductors in different topological sectors characterized by higher
monopole charges, $|q_{p}|>1$. Although the topological sector is different,
the same symmetry principles to extract $l_{z,\mathrm{gob}}$ outlined in Sec.
III.2 still hold. We begin with the general model of the monopole
superconductor in Eq. (13) but now consider a multi-Weyl semimetal with inter-
Fermi surface pairing between two parity-related Fermi surfaces. For three-
dimensional Weyl semimetals on a lattice with $n$-fold rotational symmetry
($n=2,3,4,6$), there can exist single, double, or triple Weyl nodes [38, 39].
The lattice point group symmetry sets an upper bound to the chirality of the
Weyl node, with $|\nu|\leq[n/2]$, and similarly for the pair monopole charge,
$|q_{p}|\leq n/2$.

As an example, we study monopole pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=3,l_{z}=0)}$ that has pair monopole charge
$q_{p}=3$ and transforms according to global angular momentum $l_{z}=0$. Such
pairing order can arise in a doped triple Weyl semimetal with inter-Fermi
surface $s$-wave pairing. We consider a minimal model of a triple Weyl
semimetal on a three-dimensional hexagonal lattice, which hosts triple Weyl
nodes protected by $C_{6}$ symmetry (see Appendix E.2). At appropriate doping,
the system has two parity-related Fermi surfaces, $\mathrm{FS}_{1}$ and
$\mathrm{FS}_{2}$, with Chern numbers $\mathcal{C}_{1}=3$ and
$\mathcal{C}_{2}=-3$.

From the tight-binding model, we consider the continuum model of the band
Hamiltonians $\mathcal{H}_{\mathrm{kin},1}^{(\nu=-3)}$ and
$\mathcal{H}_{\mathrm{kin},2}^{(\nu=+3)}$ describing Fermi surfaces
$\mathrm{FS}_{1}$ and $\mathrm{FS}_{2}$ surrounding triple Weyl nodes at
$\pm\mathbf{K}_{0}=(0,0,\pm K_{0,z})$. The low-energy band Hamiltonian defined
near $\mathbf{K}_{0}$ is given by

$\displaystyle\mathcal{H}^{(\nu=-3)}_{\mathrm{kin},1}(\tilde{\mathbf{k}}+\mathbf{K}_{0})=-\mu+\mathbf{h}_{1}^{(\nu=-3)}(\tilde{\mathbf{k}})\cdot\bm{\sigma}$
(36) $\displaystyle=-\mu+\hbar v_{F}\tilde{k}_{z}\sigma_{z}+\hbar
v_{F}\left[(\tilde{k}_{x}+i\tilde{k}_{y})^{3}\sigma_{+}+\mathrm{h.c.}\right],$

in which $\sigma_{\pm}=\left(\sigma_{x}\pm i\sigma_{y}\right)/2$, and is
related to the other Fermi surface by parity,
$\mathcal{H}_{\mathrm{kin},1}^{(\nu=-3)}(\mathbf{k})=\sigma_{z}\mathcal{H}_{\mathrm{kin},2}^{(\nu=+3)}(-\mathbf{k})\sigma_{z}$.
The helical eigenstate at $\mathrm{FS}_{1}$ is
$\chi_{1,+}(\tilde{\mathbf{k}})=(\cos({\theta_{h_{1}}}/{2})e^{-i(\lambda-1)\varphi_{h_{1}}/2},\sin({\theta_{h_{1}}}/{2})e^{-i(\lambda+1)\varphi_{h_{1}}/2})$,
and at $\mathrm{FS}_{2}$, it is
$\chi_{2,+}(\tilde{\mathbf{k}})=\sigma_{z}\chi_{1,+}(-\tilde{\mathbf{k}})$.
Here, consistent with previous notation, $\lambda=\pm 1$ denotes the choice of
gauge, and $\theta_{h_{1}}$ and $\varphi_{h_{1}}$ are the polar and azimuthal
angles defined in $\mathbf{h}_{1}^{(\nu=-3)}$ spherical coordinates by
$\mathbf{h}_{1}^{(\nu=-3)}\equiv|h_{1}^{(\nu=-3)}|(\sin\theta_{h_{1}}\cos\varphi_{h_{1}},\sin\theta_{h_{1}}\sin\varphi_{h_{1}},\cos\theta_{h_{1}})$.
They are related to momentum space polar and azimuthal angles by
$\sin\theta_{h_{1}}=\sin^{3}\theta_{\tilde{k}}/{(\cos^{2}\theta_{\tilde{k}}+\sin^{6}\theta_{\tilde{k}})^{1/2}}$
and $\varphi_{h_{1}}=-3\varphi_{\tilde{k}}$. We study the case of inter-Fermi
surface $s$-wave pairing in Eq. (13), in which
$\Delta_{\mathrm{inter}}(\mathbf{k})=i\Delta_{0}\sigma_{y}$, with $\Delta_{0}$
being the pairing amplitude.

In the weak-coupling regime, we obtain the low energy band-projected pairing
order from Eq. (16) which takes the following form,

$\displaystyle\Delta_{\mathrm{MSC}}^{{(q_{p}=3,l_{z}=0)}}(\tilde{\mathbf{k}})$
$\displaystyle=-\Delta_{0}\sqrt{\frac{8\pi}{3}}\mathcal{Y}_{-1;1,0}(\theta_{h_{1}},\varphi_{h_{1}})$
(37)
$\displaystyle=-\Delta_{0}\frac{\sin^{3}\theta_{\tilde{k}}}{\sqrt{\cos^{2}\theta_{\tilde{k}}+\sin^{6}\theta_{\tilde{k}}}}e^{i3\lambda\varphi_{\tilde{k}}}.$

As the mapping from the spherical coordinates in $\mathbf{h}_{1}^{(\nu=-3)}$
space to those in $\tilde{\mathbf{k}}$ space is characterized by a winding
number of $-3$, the $q=-1$ monopole harmonic order defined here in
$\mathbf{h}_{1}^{(\nu=-3)}$ spherical coordinates implies that the pairing
order in momentum space is characterized by $q_{p}=3$ pair monopole charge,
which agrees with $q_{p}=(\mathcal{C}_{1}-\mathcal{C}_{2})/2=3$. As a result,
local angular momentum is given by $l_{z,\mathrm{loc}}=\pm 3$, while the
global angular momentum is $l_{z}=0$. Furthermore, the pair monopole charge
sets the lower bound of the partial wave channels, $l\geq 3$. Different
partial wave channels of $l\geq 3$ contribute to the pairing order,
$\Delta_{\mathrm{MSC}}^{{(q_{p}=3,l_{z}=0)}}(\tilde{\mathbf{k}})=\sum_{l\geq
3}\alpha_{l}\mathcal{Y}_{q_{p}=3;l;l_{z}=0}(\theta_{\tilde{k}},\varphi_{\tilde{k}})$
with decreasing values of superposition coefficients $\alpha_{l}$ as $l$
increases from $3$.

We now examine the Josephson coupling between the above $q_{p}=3$ monopole
superconductor and a spin singlet superconductor, following the methods
introduced in Sec. III.1. For comparison, we consider the case in which there
is bare spin-independent tunneling at the junction barrier. The form factor
for this system evaluates to

$\displaystyle\mathfrak{F}^{{(q_{p}=3)}}_{\mathrm{sing}-\mathrm{MSC}}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-4\pi{\Delta_{L,0}}(\mathbf{k}^{\prime})[{\Delta_{R,\mathrm{MSC}}^{{(q_{p}=3)}}(\tilde{\mathbf{k}})}]^{*}$
(38)
$\displaystyle\times{T_{0}(-\mathbf{k},-\mathbf{k}^{\prime})}{T_{0}(\mathbf{k},\mathbf{k}^{\prime})}\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{h_{1}})\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{h_{1}}),$

in which $\bm{\Omega}_{h_{1}}=(\theta_{h_{1}},\varphi_{h_{1}})$ is the solid
angle in $\mathrm{h}_{1}^{(\nu=3)}$-space and
$\Delta_{\mathrm{MSC}}^{(q_{p}=3}(\tilde{\mathbf{k}})$ is the effective
pairing order of the $q_{p}=3$ monopole superconductor. Similar to the case of
the $q_{p}=-1$ monopole superconductor, whose form factor is given in Eq.
(27), the half-integer monopole harmonics arise from the effective spin-
dependent tunneling arising from the spin-orbit coupling in the triple Weyl
semimetal. The above form factor is non-singular and can be used to determine
the first order Josephson coupling of this system.

Similar to Sec. IV.1, we further consider a simple case that the junction is
aligned along $z$ direction, the singlet superconductor is $s$-wave, and the
bare spin-independent tunneling is a constant,
$T_{0}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}$. Then, the form factor for this
junction simplifies to

$\mathfrak{F}_{s-\mathrm{MSC}}^{(q_{p}=3,l_{z}=0)}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=T_{0}^{2}\Delta_{L,0}\Delta_{R,0}\frac{\sin^{6}\theta_{\tilde{k}}}{\cos^{2}\theta_{\tilde{k}}+\sin^{6}\theta_{\tilde{k}}},$
(39)

in which $\Delta_{L,0}$ and $\Delta_{R,0}$ are the pairing amplitudes of the
$s$-wave superconductor and $\Delta_{\mathrm{MSC}}^{(q_{p}=3,l_{z}=0)}$
monopole superconductor respectively.
$\mathfrak{F}_{s-\mathrm{MSC}}^{(q_{p}=3,l_{z}=0)}$ itself is invariant under
rotation $R_{z}$ and parity, as both the $s$-wave superconductor and
$\Delta_{\mathrm{MSC}}^{(q_{p}=3,l_{z}=0)}$ monopole superconductor both form
irreps of $l_{z}=0$. As discussed in Sec. IV.1.1, this leads to nonvanishing
first order Josephson current. Hence, following the junction design principles
introduced in Sec. III.2, we may similarly design probes of monopole pairing
order in different topological sectors.

##  V Josephson junctions to probe $l_{z,\mathrm{loc}}$ of the monopole
superconducting order

We now explore how to design Josephson junctions to probe local angular
momentum $l_{z,\mathrm{loc}}$ of monopole superconducting order, following the
second junction design introduced in Sec. III.2. We consider two identical
monopole superconductors forming a Josephson junction oriented along, for
example, the $x$ direction, which is perpendicular to the high-symmetry
rotational axis along $z$. As discussed in Sec. III.2, the total Josephson
current can be regarded as a superposition of contributions from a family of
Josephson junctions between two 2D chiral superconductors, corresponding to
different values of conserved momentum $k_{z}$. It has been studied
extensively in, for example, Refs. 23, 29 that a planar junction between two
chiral superconductors can host $4\pi$-periodic Josephson current arising from
the tunneling of chiral boundary modes. Therefore, by studying the Josephson
energy phase relation in our system, we can obtain the information about local
symmetry of monopole pairing order. In the continuum model of a
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor, its
pairing order, described by Eq. (18), exhibits local chiral $p$-wave pairing
with $l_{z,\mathrm{loc}}=\mp 1$ for a given value of $\tilde{k}_{z}$,

$\Delta_{\mathrm{MSC}}^{{(q_{p}=-1,l_{z}=0)}}(\tilde{\mathbf{k}})=\tilde{\Delta}(\tilde{k}_{z})(\tilde{k}_{x}-i\lambda\tilde{k}_{y}).$
(40)

Here, the local pairing amplitude is modulated by $\tilde{k}_{z}$ as
$\tilde{\Delta}(\tilde{k}_{z})=-\Delta_{0}\sqrt{1-(\tilde{k}_{z}/\tilde{k})^{2}}$,
and $\lambda=\pm 1$ describes the local chirality of the pairing order, which
changes sign as $\tilde{k}_{z}$ varies from the north pole to south pole,
indicative of the nontrivial topological charge of the pairing order,
$q_{p}=-1$.

We next consider a Josephson junction along the $x$ direction between two
monopole superconductors with identical pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$, shown schematically in Fig. 5
(a). Periodic boundary conditions are taken in $y$ and $z$ directions so that
momenta $\tilde{k}_{y}$ and $\tilde{k}_{z}$ transverse to the junction barrier
are conserved. At $\tilde{k}_{y}=0$, the above monopole superconducting
pairing order can be reduced to a family of $p_{x}$-wave pairing orders
locally at different values of $\tilde{k}_{z}$ with pairing amplitudes
$\tilde{\Delta}(\tilde{k}_{z})$, each of which support a Majorana state
localized at the junction interface contributing $4\pi$-periodic Josephson
energy phase relation, $E_{J}\propto\cos(\phi/2)$. Therefore, we expect the
total Josephson current phase relation in our system to exhibit
$4\pi$-periodicity as a result of the local chiral $p$-wave pairing symmetry
of the $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole pairing order.

Figure 5:  (a) Single Josephson junction along the $x$ direction between two
identical monopole superconductors with pairing order
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ and the same crystal orientation.
(b) Same Josephson junction as in Fig. (a), only now that the superconductor
at the left (yellow) is rotated by an angle $\theta$ about the $z$ axis
pointing out of plane. (c) Bulk Fermi surface $\mathrm{FS}_{1}$ of the doped
Weyl semimetal on a cubic lattice. Representative momenta
$\tilde{\mathbf{k}}_{i=1,\cdots,5}$ are shown, with $\tilde{k}_{y}=0$ and
$\tilde{k}_{z}=-0.2,-0.1,0.0,0.1,0.15$, respectively. As a reference, we plot
the largest $\tilde{k}_{y}\tilde{k}_{z}$ cross sections ($\tilde{k}_{x}=0$) of
bulk Fermi surfaces, which are identical for the two sides of the junction.
(d) Josephson energy phase relations for states at the representative
transverse momenta (in color) and states at the same values of $\tilde{k}_{z}$
but $k_{y}=\pm 0.1,\pm 0.2$ (in gray). Here, the system size and other
parameters are the same as those taken in Fig. 3.

We study a tight-binding model of the Josephson junction between two identical
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductors on a
cubic lattice, similar to that discussed in Sec. IV.1.2, but with the junction
now aligned along the $x$ direction, as shown schematically in Fig. 5 (a). Due
to the mirror $M_{z}$ symmetry of the system, there are identical
contributions to the Josephson energy phase relation from $\mathrm{FS}_{1}$
and $\mathrm{FS}_{2}$. Thus, we first examine states at five representative
transverse momentum $\tilde{\mathbf{k}}_{i}=(0,\tilde{k}_{i})$ ($i=1,\dots,5$)
within $\mathrm{FS}_{1}$ in the $k_{y}k_{z}$ plane, as shown in Fig. 5 (c). At
each conserved $\tilde{\mathbf{k}}_{i}$, the effective 1D Josephson junction
along $x$ direction hosts midgap states localized at the junction interface,
which contribute to $4\pi$-periodic Josephson energy phase relation, as shown
in Fig. 5 (d). As the magnitude of $\tilde{k}_{z}$ increases towards the
locations of pairing gap nodes, the effective pairing amplitude,
$\tilde{\Delta}(\tilde{k}_{z})$, decreases. Furthermore, at phase difference
$\phi=\pi$, the junction interface plays the role of a domain wall, where the
pairing order changes its sign at the junction barrier. For each
representative transverse momentum at $\tilde{k}_{y}=0$, the pairing order is
odd under mirror $M_{x}$ reflection,
${\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}(\tilde{k}_{x},\tilde{k}_{y}=0,\tilde{k}_{z})}={-\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}(-\tilde{k}_{x},\tilde{k}_{y}=0,\tilde{k}_{z})}$.
Together with the BdG particle-hole symmetry, this system is analogous to the
Jackiw-Rebbi model [40] and exhibits zero energy Andreev states localized at
junction interface when $\phi=\pi$. When fermion parity is conserved at
$\phi=\pi$, the total energy phase relation is $4\pi$-periodic and leads to
Josephson current phase relation of the same periodicity, shown in Fig. 6 (a)
in red. As a result, when a nonzero external magnetic flux is varied through
the junction interface, the Josephson Fraunhofer pattern exhibits lifted nodes
at odd integer flux quanta, shown in red in Fig. 6 (b). The conserved fermion
parity can be achieved experimentally, for example, by taking a dynamic
measurement of the Josephson current while maintaining a fixed voltage to
prevent parity transitions [41, 42, 43].

Figure 6:  Josephson current of the junction shown in Fig. 5 without and with
magnetic flux through the junction interface. (a) The total Josephson current
phase relation and (b) the Fraunhofer pattern when conserving (red, solid) or
breaking (blue, dashed) fermion parity. Here, $I_{0}={e\Delta_{0}}/{\hbar}$,
$\Phi_{0}={h}/{2e}$, and $\Delta_{0}/t=0.04$.

The existence of $4\pi$-periodic Josephson energy phase relation in this
single Josephson junction oriented along the $x$ direction stems from the
local spatial symmetry property of the monopole pairing order, which contains
nonzero component of $p_{x}$-wave pairing, or say, nonzero components of
$l_{z,\mathrm{loc}}=\pm 1$. Furthermore, due to spin-orbit coupling from Weyl
band Hamiltonian, we examine the pairing order under combined spin and spatial
rotations. The $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole
superconductor is invariant under simultaneous spin and spatial rotations;
hence, the local $z$-component angular momentum, $l_{z,\mathrm{loc}}$, in
momentum space and of spin, $s_{z}$, of the pairing order satisfies
$l_{z,\mathrm{loc}}+s_{z}=0$. However, the existence of $4\pi$-periodic
Josephson current constrains $l_{z,\mathrm{loc}}\neq 0$, i.e. the monopole
superconducting order has the symmetry of a chiral superconducting order
locally in momentum space. Specifically, at conserved momentum
$\tilde{k}_{z}$, the above junction can only be described in analogy to that
between two superconductors with spin-polarized chiral $p$-wave pairing order,
in which the Cooper pair is locally described by
$(p_{x}+ip_{y})|\downarrow\downarrow\rangle$ or
$(p_{x}-ip_{y})|\uparrow\uparrow\rangle$ pairing, satisfying
$l_{z,\mathrm{loc}}=-s_{z}$. Later, in Sec. VI, we show that the
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconducting order is
indeed unique and different from these chiral superconductors globally.
Nonetheless, though the monopole superconducting order is inherently different
from any pairing based on spherical harmonic symmetry, at given
$\tilde{k}_{z}$ one can describe the symmetry of the monopole pairing order
locally in analogy to that of a 2D chiral superconducting pairing order. As
such, for given conserved $\tilde{k}_{z}$, the
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ superconductor can be described
locally in momentum space as having chiral $p$-wave symmetry.

If the symmetry of the spin-orbit coupled system is not known a priori, we
additionally examine a second Josephson junction to demonstrate the local
symmetry of pairing order. Consider the junction geometry in which the
monopole superconductor to the left of the junction barrier is rotated by an
angle $\theta$ about the $z$ axis, as shown in Fig. 5 (b). If the Josephson
energy phase relation is invariant under the rotation, exhibiting
$4\pi$-periodicity with no phase shift, the pairing order lies in the zero
total angular momentum sector with a nonzero orbital angular momentum
$l_{z,\mathrm{loc}}$. We check this numerically using the tight-binding model,
where the four-fold lattice rotation takes discrete values of $\theta=n\pi/2$
with $n\in\mathbb{Z}$, and obtain the same Josephson energy phase relation as
in Fig. 5 (d). Then, following the same reasoning as in the paragraph above,
this additional rotated junction demonstrates the local chiral $p$-wave
symmetry of the $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole pairing
order.

Generally, the monopole superconducting order is always nodal, and, locally in
momentum space, its pairing order can always be described by chiral spherical
harmonic pairing symmetry with angular momentum $l_{z,\mathrm{loc}}$ (see
Appendix C for details). We can consider the same Josephson junction geometry
to identify the local pairing symmetry of other monopole orders that have a
different pair monopole charge or a different angular momentum
$l_{z,\mathrm{glob}}$.

##  VI Josephson junctions to probe monopole superconducting order in the
presence of spin-orbit interactions

Lastly, in the presence of spin-orbit interactions at the junction interface
or spin-orbit coupling in the bulk, the spin and orbital angular momentum are
no longer conserved separately. The Josephson coupling is allowed when the
pairing orders at two sides of the junction share the same value of total
angular momentum along the common rotational symmetry axis following the
junction geometry. For example, consider an SIS Josephson junction along $z$
axis between an $s$-wave superconductor with pairing
$\Delta_{L,0}|\uparrow\downarrow-\downarrow\uparrow\rangle$ and a chiral
$p$-wave triplet superconductor with zero total angular momentum pairing, for
example
$\Delta_{R,0}[(p_{x}-ip_{y}){|\uparrow\uparrow\rangle}+(p_{x}+ip_{y}){|\downarrow\downarrow\rangle}]$.
In the presence of spin-orbit interactions, it shows nonzero first order
Josephson coupling, similar to that in the Josephson junction of the same
geometry between an $s$-wave superconductor and a monopole superconductor with
pairing order $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$. Consequently,
junction design (I) introduced in Sec. IV cannot distinguish
$\Delta^{(q_{p}=-1,l_{z}=0)}_{\mathrm{MSC}}$ monopole pairing order from
chiral $p$-wave pairing with $j_{z}=0$. Furthermore, junction design (II)
introduced in Sec. V similarly cannot distinguish them, as the Majorana
surface modes localized at the junction interface between chiral $p$-wave
superconductors also give rise $4\pi$-periodic Josephson current phase
relation. Hence, we must resort to a third junction design to discern the
$\mathrm{\Delta}_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ superconducting order in
the presence of spin-orbit interaction at the barrier.

As an example, we consider the effect of bare interface Rashba-type spin-orbit
interaction. Generally, when spin-orbit interaction preserves time-reversal
symmetry, the form factor between a spin singlet and a $q_{p}=-1$ monopole
superconductor takes the following form,

$\mathfrak{F}_{\mathrm{\mathrm{sing}-\mathrm{MSC}}}^{(q_{p}=-1)}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=\mathfrak{f}_{0}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})+\mathfrak{f}_{1}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})+\mathfrak{f}_{2}(\tilde{\mathbf{k}},\mathbf{k}^{\prime}),$
(41)

where

$\displaystyle\mathfrak{f}_{0}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-4\pi\Delta_{L,0}(\mathbf{k}^{\prime})[\Delta_{R,\mathrm{MSC}}^{(q_{p}=-1)}(\tilde{\mathbf{k}})]^{*}$
(42a) $\displaystyle\times
T_{0}^{2}(\mathbf{k},\mathbf{k}^{\prime})\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}}),$
$\displaystyle\mathfrak{f}_{1}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-4\pi\Delta_{L,0}(\mathbf{k}^{\prime})[\Delta_{R,\mathrm{MSC}}^{(q_{p}=-1)}(\tilde{\mathbf{k}})]^{*}T_{0}(\mathbf{k},\mathbf{k}^{\prime})$
(42b)
$\displaystyle\times\Big{\\{}\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},+\frac{1}{2}}^{2}(\bm{\Omega}_{\tilde{k}})[T_{x}(\mathbf{k},\mathbf{k}^{\prime})-iT_{y}(\mathbf{k},\mathbf{k}^{\prime})]$
$\displaystyle\hskip
10.00002pt+\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}^{2}(\bm{\Omega}_{\tilde{k}})[T_{x}(\mathbf{k},\mathbf{k}^{\prime})+iT_{y}(\mathbf{k},\mathbf{k}^{\prime})]\Big{\\}},$
and
$\displaystyle\mathfrak{f}_{2}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-4\pi\Delta_{L,0}(\mathbf{k}^{\prime})[\Delta_{R,\mathrm{MSC}}^{(q_{p}=-1)}(\tilde{\mathbf{k}})]^{*}$
(42c)
$\displaystyle\times|\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})|^{2}\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},-\frac{1}{2}}(\bm{\Omega}_{\tilde{k}})\mathcal{Y}_{-\frac{1}{2};\frac{1}{2},\frac{1}{2}}(\bm{\Omega}_{\tilde{k}}).$

Here, $\mathfrak{f}_{0}$, $\mathfrak{f}_{1}$, and $\mathfrak{f}_{2}$ denote
the contributions to the form factor that are independent of, linear in, and
quadratic in spin-orbit interaction amplitude, respectively. We take
$\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=-\mathbf{T}(-\mathbf{k},-\mathbf{k}^{\prime})$
and
$T_{0}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}(-\mathbf{k},-\mathbf{k}^{\prime})$
so that the tunneling Hamiltonian in Eq. (3) has time-reversal symmetry.

Now we consider a Josephson junction along $x$ direction between an $s$-wave
superconductor and the $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole
superconductor, as shown schematically in Fig. 7 (a). In the presence of
Rashba-type spin-orbit interaction at the junction interface,
$\mathbf{T}(\mathbf{k},\mathbf{k}^{\prime})=T_{\mathrm{SO}}\delta_{\mathbf{k}_{\parallel},\mathbf{k}^{\prime}_{\parallel}}\hat{\mathbf{x}}\times\hat{\mathbf{k}}$.
For simplicity, we take the spin-independent tunneling amplitude as a
constant, $T_{0}(\mathbf{k},\mathbf{k}^{\prime})=T_{0}$. Terms in the form
factor in Eq. (41) that contribute to nonvanishing first order Josephson
coupling are given by
$\mathfrak{f}_{0}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-T_{0}^{2}\Delta_{L,0}\Delta_{R,0}\sin^{2}\theta_{\tilde{k}}$
and
$\mathfrak{f}_{2}(\tilde{\mathbf{k}},\mathbf{k}^{\prime})=-T_{\mathrm{SO}}^{2}|\hat{\mathbf{x}}\times\hat{\mathbf{k}}|^{2}\Delta_{L,0}\Delta_{R,0}\sin^{2}\theta_{\tilde{k}}\delta_{\mathbf{k}_{\parallel},\mathbf{k}^{\prime}_{\parallel}},$
with $\Delta_{L,0}$ and $\Delta_{R,0}$ being the pairing amplitudes of the
$s$-wave and $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole
superconductor respectively. Specifically, contributions from
$\mathfrak{f}_{0}$ and $\mathfrak{f}_{2}$ are invariant under rotation $R_{x}$
and thus lead to nonzero first order Josephson current between the
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor and
$s$-wave superconductor in a junction along the $x$ direction.

In contrast, a
$(p_{x}-ip_{y}){|\uparrow\uparrow\rangle}+(p_{x}+ip_{y}){|\downarrow\downarrow\rangle}$
chiral $p$-wave superconductor has vanishing first order Josephson coupling
with an $s$-wave superconductor in the same junction geometry, shown
schematically in Fig. 7 (c), and with the same Rashba-type spin-orbit
interaction at the interface. The form factor in Eq. (10) simplifies to
$\mathfrak{F}_{\mathrm{sing}-\mathrm{trip}}\propto
T_{0}T_{\mathrm{SO}}\hat{x}\cdot(\hat{\mathbf{k}}\times\mathbf{d}_{R}^{*}(\mathbf{k}))$.
Here, the $d$-vector in consideration is
$\mathbf{d}_{R}(\mathbf{k})\propto(k_{y},-k_{x},0)$. The form factor reduces
to
$\mathfrak{F}_{\mathrm{sing}-\mathrm{trip}}(\mathbf{k},\mathbf{k}^{\prime})\propto
k_{x}k_{z}$, which is odd under mirror $M_{z}$ reflection, $k_{z}\mapsto-
k_{z}$. Therefore, the net first order Josephson current of this system is
zero. As such, the presence of first order Josephson current in this
additional junction design can be used to distinguish
$(p_{x}-ip_{y}){|\uparrow\uparrow\rangle}+(p_{x}+ip_{y}){|\downarrow\downarrow\rangle}$
chiral superconductor from the $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$
monopole superconductor. For other monopole superconducting pairing orders in
different $q_{p}$ or $l_{z}$ sectors, one can design similar junctions to
distinguish monopole superconducting order from chiral pairing orders in the
presence of bare spin-orbit interactions at the interface.

Figure 7:  Tight-binding results of Josephson current for a junction with an
$s$-wave superconductor in the presence of Rashba-type spin-orbit interaction
at the junction interface. (a) Schematic of single Josephson junction between
an $s$-wave superconductor and $\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$
monopole superconductor and along the $x$ direction. (b) $2\pi$-periodic
Josephson current phase relation (red, solid) for the junction shown in (a).
(c) Schematic of junction between an $s$-wave superconductor and
$(p_{x}-ip_{y})\left|\uparrow\uparrow\right\rangle+(p_{x}+ip_{y})\left|\downarrow\downarrow\right\rangle$
superconductor. (d) $\pi$-periodic Josephson current phase relation (black,
dashed) of the junction shown in (c). Results are computed using a system of
$2N_{x}=300$ sites and with the same parameters of the tight-binding band
Hamiltonians in Fig. 3. For simplicity, $t_{\mathrm{SO}}=t_{0}$ at the
interface.

We demonstrate the analytic results above using numerical calculation of the
Josephson current based on tight binding models. The system is analogous to
that introduced in Sec. IV.1.2, only now the junction is oriented along $x$
direction and includes spin-dependent tunneling at the interface. For given
conserved transverse momenta $k_{y}$ and $k_{z}$, we model the junction
barrier as

$\displaystyle H_{\mathrm{link}}(k_{y},k_{z})=$ (43)
$\displaystyle\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}c^{\dagger}_{k_{y},k_{z},\sigma,N_{x}+1}\left[\hat{t}(k_{y},k_{z})\right]_{\sigma,\sigma^{\prime}}c_{k_{y},k_{z},\sigma^{\prime},N_{x}}+\mathrm{h.c.}$

Above,
$\hat{t}(k_{y},k_{z})=t_{0}\sigma_{0}+\mathbf{t}(k_{y},k_{z})\cdot\bm{\sigma}$,
in which $t_{0}$ is the spin-independent hopping amplitude, and
$\mathbf{t}(k_{y},k_{z})=t_{\mathrm{SO}}(\sin k_{y}\hat{z}-\sin k_{z}\hat{y})$
is the Rashba-type spin-orbit coupling at the interface, with
$t_{\mathrm{SO}}$ being the spin-dependent hopping amplitude.

The Josephson current for the junction along $x$ direction between a
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ monopole superconductor and
$s$-wave superconductor is shown in Fig. 7 (b). In agreement with the above
arguments using the form factor
$\mathfrak{F}_{\mathrm{sing}-\mathrm{MSC}}^{(q_{p}=-1)}$, there is an overall
$2\pi$-periodic Josephson current corresponding to nonvanishing first order
Josephson coupling. In contrast, for a junction between an $s$-wave
superconductor and
$(p_{x}-ip_{y}){|\uparrow\uparrow\rangle}+(p_{x}+ip_{y})\left|\downarrow\downarrow\right\rangle$
chiral superconductor, there is no first order Josephson coupling, as shown in
Fig. 7 (d). Rather, the leading order Josephson coupling between the $s$-wave
and chiral $p$-wave superconductor is a second order process, corresponding to
$\pi$-periodic Josephson current. It is to be noted that because we have taken
the spin-independent and spin-dependent tunneling amplitudes to both be of
order $1$ in the tight-binding model, the magnitude of the Josephson critical
current for the two junctions is similar; however, this does not change the
qualitative aspects, specifically the periodicity. The numerical results are
in agreement with the above microscopic analysis and demonstrate that this
additional junction can be used to distinguish
$\Delta_{\mathrm{MSC}}^{(q_{p}=-1,l_{z}=0)}$ pairing order from other chiral
$p$-wave superconductors in the presence of spin-orbit interaction at the
interface. Similar results hold for a chiral
$(p_{x}+ip_{y})|\downarrow\downarrow\rangle$ superconductor,
$(p_{x}-ip_{y})\left|\uparrow\uparrow\right\rangle$ superconductor, or linear
superposition of the two. In other words, the monopole superconducting order
cannot be described as a simple linear superposition of spin-polarized chiral
superconducting pairing orders globally but rather exists in its own distinct
topological class.

##  VII Conclusion

We have proposed using a set of phase-sensitive Josephson junctions to probe
monopole superconducting orders and to distinguish them from other known
pairing orders with spherical harmonic symmetry. We develop the designs based
on two approaches. Firstly, we employed the linear response approach to derive
the first order Josephson current between two uniform superconductors. In
addition to accounting for the microscopic details of the tunneling processes,
including spin-orbit interactions at the junction interface and spin-orbit
coupling in the bulk superconductors, we emphasize that the symmetry of
superconducting pairing orders and the junction geometry is encoded in a form
factor $\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})$, which is nonsingular
and determines the symmetry selection rules of the first order Josephson
coupling. Complementing the first approach, we have also proposed symmetry and
topological principles in the design of phase-sensitive probes of monopole
superconducting order. We noted that the $\mathrm{U}(1)$ topological
obstruction of the pairing phase can be exhibited as a nonzero shift between
the global and local angular momentum of the monopole superconducting order.
Then, following symmetry principles, we proposed two classes of Josephson
junctions of monopole superconductors to probe the global and local angular
momentum of the superconducting order.

We considered the first class of Josephson junctions aligned along the
rotational axis ($z$ axis) between a monopole superconductor and other
superconductors with known pairing symmetry. As the first order
$2\pi$-periodic Josephson current is allowed only when superconducting orders
at both sides of the junction transform under rotation $R_{z}$ in the same
way, we can use this to probe $l_{z,\mathrm{glob}}$, the global angular
momentum of the monopole pairing order. We have studied Josephson currents in
this class of junctions between monopole superconductors with different pair
monopole charges and different $l_{z,\mathrm{glob}}$ and superconductors with
standard $s$\- and $d$-wave spherical harmonic symmetry, using both linear
response and tight-binding calculations. The second class of junctions,
designed to extract the local angular momentum $l_{z,\mathrm{loc}}$ of the
pairing order, consists of two identical monopole superconductors in a
junction aligned perpendicular to the rotational $z$ axis. The Josephson
current can be viewed as a superposition of contributions from different
values of conserved momentum $k_{z}$ local in momentum space. By analyzing the
existence of $4\pi$-periodic Josephson current at different orientations of
the monopole superconductor based on tight-binding models, we have identified
the local chiral symmetry of monopole pairing order characterized by angular
momentum $l_{z,\mathrm{loc}}$. The shift between $l_{z,\mathrm{glob}}$ and
$l_{z,\mathrm{loc}}$ indicates the nonzero pair monopole charge of the exotic
superconducting order. The underlying symmetry and topological principles of
this work provide useful guidance for the ongoing experimental efforts toward
realizing and identifying exotic monopole superconducting pairing orders in
doped topological semimetal materials.

## Acknowledgements

We acknowledge the support of the NSF CAREER Grant No. DMR-1848349. In
addition, YL acknowledges the Institute for Quantum Matter, an Energy Frontier
Research Center funded by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, under Award DESC0019331, for support during
the initial stage of this work. JYZ acknowledges partial support from the
Johns Hopkins University Theoretical Interdisciplinary Physics and Astronomy
Center.

## References

  * Li and Haldane [2018] Y. Li and F. D. M. Haldane, Phys. Rev. Lett. 120, 067003 (2018).
  * Murakami and Nagaosa [2003] S. Murakami and N. Nagaosa, Phys. Rev. Lett. 90, 057002 (2003).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Baidya and Vanderbilt [2020] S. Baidya and D. Vanderbilt, Phys. Rev. B 102, 165115 (2020).
  * Jeon  _et al._ [2021] K.-R. Jeon, B. K. Hazra, K. Cho, A. Chakraborty, J.-C. Jeon, H. Han, H. L. Meyerheim, T. Kontos, and S. S. P. Parkin, Nature Materials 20, 1358 (2021).
  * Xiao  _et al._ [2022] R. Xiao, J. T. Held, J. Rable, S. Ghosh, K. Wang, K. A. Mkhoyan, and N. Samarth, Phys. Rev. Materials 6, 024203 (2022).
  * Li [2020] Y. Li,  (2020), arXiv:2001.05984 .
  * Bobrow and Li [2022] E. Bobrow and Y. Li, (2022), 2204.04249v1 .
  * Read and Green [2000] N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * Volovik [2003] G. E. Volovik,  _The Universe in a Helium Droplet_ (Oxford University Press, 2003).
  * Wu and Yang [1976] T. T. Wu and C. N. Yang, Nucl. Phys. B 107, 365 (1976).
  * Haldane [1983] F. D. M. Haldane, Phys. Rev. Lett. 51, 605 (1983).
  * Goddard and Olive [1978] P. Goddard and D. I. Olive, Reports on Progress in Physics 41, 1357 (1978).
  * Atiyah and Hitchin [1985] M. F. Atiyah and N. J. Hitchin, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 315, 459 (1985).
  * Josephson [1962] B. Josephson, Physics Letters 1, 251 (1962).
  * Geshkenbein  _et al._ [1987] V. B. Geshkenbein, A. I. Larkin, and A. Barone, Physical Review B 36, 235 (1987).
  * Sigrist and M. Rice [1992] M. Sigrist and T. M. Rice, Journal of the Physical Society of Japan 61, 4283 (1992).
  * Wollman  _et al._ [1993] D. A. Wollman, D. J. Van Harlingen, W. C. Lee, D. M. Ginsberg, and A. J. Leggett, Phys. Rev. Lett. 71, 2134 (1993).
  * Tsuei  _et al._ [1994] C. C. Tsuei, J. R. Kirtley, C. C. Chi, L. S. Yu-Jahnes, A. Gupta, T. Shaw, J. Z. Sun, and M. B. Ketchen, Phys. Rev. Lett. 73, 593 (1994).
  * Van Harlingen [1995] D. J. Van Harlingen, Rev. Mod. Phys. 67, 515 (1995).
  * Tsuei and Kirtley [2000] C. C. Tsuei and J. R. Kirtley, Rev. Mod. Phys. 72, 969 (2000).
  * Kitaev [2001] A. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Kwon  _et al._ [2004] H.-J. Kwon, K. Sengupta, and V. M. Yakovenko, The European Physical Journal B 37, 349 (2004).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Chung and Zhang [2009] S. Chung and S. Zhang, Phys. Rev. Lett. 103, 235301 (2009).
  * Sau  _et al._ [2010] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. D. Sarma, Physical Review Letters 104 (2010), 10.1103/physrevlett.104.040502.
  * Alicea [2010] J. Alicea, Phys. Rev. B 81, 125318 (2010).
  * Grosfeld and Stern [2011] E. Grosfeld and A. Stern, Proceedings of the National Academy of Sciences 108, 11810 (2011).
  * Alicea [2012] J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
  * Ambegaokar and Baratoff [1963] V. Ambegaokar and A. Baratoff, Phys. Rev. Lett. 10, 486 (1963).
  * Sigrist and Ueda [1991] M. Sigrist and K. Ueda, Rev. Mod. Phys. 63, 239 (1991).
  * Bardeen [1961] J. Bardeen, Physical Review Letters 6, 57 (1961).
  * Cohen  _et al._ [1962] M. H. Cohen, L. M. Falicov, and J. C. Phillips, Physical Review Letters 8, 316 (1962).
  * Balian and Werthamer [1963] R. Balian and N. Werthamer, Phys. Rev. 131, 1553 (1963).
  * Leggett [1975] A. J. Leggett, Rev. Mod. Phys. 47, 331 (1975).
  * Geshkenbein and Larkin [1986] V. Geshkenbein and A. Larkin, JETP Lett. 43, 395 (1986).
  * Wu and Yang [1977] T. T. Wu and C. N. Yang, Phys. Rev. D 16, 1018 (1977).
  * Fang  _et al._ [2012] C. Fang, M. J. Gilbert, X. Dai, and B. A. Bernevig, Phys. Rev. Lett. 108, 266802 (2012).
  * Tsirkin  _et al._ [2017] S. S. Tsirkin, I. Souza, and D. Vanderbilt, Physical Review B 96 (2017), 10.1103/physrevb.96.045102.
  * Jackiw and Rebbi [1976] R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * Kurter  _et al._ [2015] C. Kurter, A. Finck, Y. S. Hor, and D. J. Van Harlingen, Nature Communications 6, 7130 (2015).
  * Abboud  _et al._ [2022] N. Abboud, V. Subramanyan, X.-Q. Sun, G. Yue, D. V. Harlingen, and S. Vishveshwara, Physical Review B 105 (2022), 10.1103/physrevb.105.214521.
  * Yue  _et al._ [2024] G. Yue, C. Zhang, E. D. Huemiller, J. H. Montone, G. R. Arias, D. G. Wild, J. Y. Zhang, D. R. Hamilton, X. Yuan, X. Yao, D. Jain, J. Moon, M. Salehi, N. Koirala, S. Oh, and D. J. Van Harlingen, Physical Review B 109 (2024), 10.1103/physrevb.109.094511.
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ (Springer US, Boston, MA, 2000).
  * Lüders and Usadel [1971] G. Lüders and K.-D. Usadel,  _The method of the correlation function in superconductivity theory_ (Springer, 1971).

##  Appendix A Linear response theory of Josephson coupling in presence of
spin-orbit interactions

We review the derivation of the first order Josephson current in an SIS
junction via linear response theory, accounting for a general form of spin-
orbit interaction in the bulk superconductors or at the junction interface. We
consider the Hamiltonian of the Josephson junction presented in Eq. (1), in
which two superconductors are connected by a weak link modelled by the
tunneling Hamiltonian in Eq. (3).

###  A.1 Evaluating first order Josephson coupling

The total electric tunneling current through the junction reads
$I(t)=|e|\langle\dot{N}_{L}(t)\rangle$, in which
$N_{L}=\sum_{\sigma,\mathbf{k}}c^{\dagger}_{L,\mathbf{k},\sigma}c_{L,\mathbf{k},\sigma}$
is the number of particles at the left side of the barrier. We evaluate using
linear response theory, in which $H_{\mathrm{link}}$ in Eq. (3) is treated
perturbatively in the thermal weight expansion to the first order. The first
order current is given by
$I(t)={-i(e/\hbar)\int_{-\infty}^{+\infty}\mathrm{d}t^{\prime}\Theta(t-t^{\prime})\langle[\dot{N}_{L}(t),H_{\mathrm{link}}(t^{\prime})]\rangle}$,
where $\Theta(t-t^{\prime})$ is the Heaviside function, and $\dot{N}_{L}(t)$
and $H_{\mathrm{link}}(t^{\prime})$ are in the interaction picture. It follows
that the current is a superposition of the normal current from quasiparticle
tunneling and the Josephson supercurrent from the tunneling of Cooper pairs
[44]. The Josephson tunneling current is given by the imaginary part of the
current-current correlation function,

$I_{J}(t,\phi)=-\frac{2e}{\hbar}\mathrm{Im}[e^{-2ieVt}\Phi_{\mathrm{ret}}(eV,\phi)],$
(44)

in which $eV=\mu_{L}-\mu_{R}$ is the voltage bias, and
$\phi=\phi_{R}-\phi_{L}$ is the overall $\mathrm{U}(1)$ phase difference of
the superconducting order parameters. Note that in the main text, we work at
zero bias, $eV=0$. Above, the correlation function is given by
$\Phi_{\mathrm{ret}}(t-t^{\prime},\phi)={-i\Theta(t-t^{\prime})}\langle[A(t),A(t^{\prime})]\rangle$,
in which

$A(t)=\sum_{\mathbf{k},\sigma;\mathbf{k}^{\prime},\sigma^{\prime}}T_{\sigma,\sigma^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})c^{\dagger}_{R,\mathbf{k},\sigma}(t)c_{L,\mathbf{k}^{\prime},\sigma^{\prime}}(t)$
(45)

is the current operator and $T(\mathbf{k},\mathbf{k}^{\prime})$ is the
tunneling matrix. It is assumed that the two states located at the left and
right sides satisfy
$\\{c_{L,\mathbf{k}^{\prime},\sigma^{\prime}},c_{R,\mathbf{k},\sigma}\\}=\\{c_{L,\mathbf{k}^{\prime},\sigma^{\prime}},c^{\dagger}_{R,\mathbf{k},\sigma}\\}=0$.
Although the above correlation function is derived from first order
perturbation theory, it corresponds to a second order process in
$T_{\sigma,\sigma^{\prime}}$, which can be regarded as the tunneling of a
Cooper pair.

We evaluate the corresponding correlation function in imaginary frequency
space,
$\Phi(i\omega,\phi)=-\int_{0}^{\beta}e^{i\omega\tau}\mathrm{d}\tau\langle\mathcal{T}_{\tau}A(\tau)A(0)\rangle$
and expand with Wick’s theorem:

$\displaystyle\Phi(i\omega,\phi)=-\sum_{\mathbf{k},\mathbf{k}^{\prime}}\sum_{\sigma,\sigma^{\prime},\rho,\rho^{\prime}}T_{\sigma,\sigma^{\prime}}(\mathbf{k},\mathbf{k}^{\prime})T_{\rho,\rho^{\prime}}(-\mathbf{k},-\mathbf{k}^{\prime})\times\frac{1}{\beta}\sum_{ip_{n}}\mathcal{F}_{R,\sigma,\rho}^{\dagger}(\mathbf{k},ip_{n})\mathcal{F}_{L,\sigma^{\prime},\rho^{\prime}}(\mathbf{k}^{\prime},ip_{n}-i\omega).$
(46)

Above, $p_{n}=\pi(2n+1)/\beta$ with $n\in\mathbb{Z}$ are the fermion Matsubara
frequencies, and
$\mathcal{F}_{L,\sigma^{\prime},\rho^{\prime}}(\mathbf{k},ip_{n}-i\omega)$ and
$\mathcal{F}^{\dagger}_{R,\sigma,\rho}(\mathbf{k},ip_{n})$ are the Fourier
transforms of the anomalous Green’s functions,
$\mathcal{F}_{L,\sigma^{\prime},\rho^{\prime}}(\mathbf{k},\tau)=\langle\mathcal{T}_{\tau}c_{L,-\mathbf{k},\rho^{\prime}}(\tau)c_{L,\mathbf{k},\sigma^{\prime}}(0)\rangle$
and
$\mathcal{F}^{\dagger}_{R,\sigma,\rho}(\mathbf{k},\tau)=\langle\mathcal{T}_{\tau}c^{\dagger}_{R,\mathbf{k},\sigma}(\tau)c^{\dagger}_{R,-\mathbf{k},\rho}(0)\rangle$,
with $\mathcal{T}_{\tau}$ being the time-ordering operator. It follows that
the correlation function can be expressed as [31]

$\Phi(i\omega,\phi)=-\frac{1}{\beta}\sum_{ip_{n}}\sum_{\mathbf{k},\mathbf{k}^{\prime}}\mathrm{Tr}\left[\mathcal{F}_{L}(\mathbf{k},ip_{n}-i\omega)T^{\mathrm{T}}(-\mathbf{k},-\mathbf{k}^{\prime})\mathcal{F}^{*}_{R}(\mathbf{k},ip_{n})T(\mathbf{k},\mathbf{k}^{\prime})\right],$
(47)

in which the trace is taken over any internal degrees of freedom (spin,
orbital, etc.).

Following Eq. (2), we consider the following mean-field BdG Hamiltonian in
$(\mathbf{c}_{\alpha,\mathbf{k}},\mathbf{c}^{\dagger}_{\alpha,-\mathbf{k}})$
Nambu basis,

$\mathcal{H}_{\mathrm{BdG},\alpha}(\mathbf{k})=\left(\begin{array}[]{cc}\mathcal{H}_{\mathrm{kin},\alpha}(\mathbf{k})&\Delta_{\alpha}(\mathbf{k})e^{i\phi_{\alpha}}\\\
\Delta^{\dagger}_{\alpha}(\mathbf{k})e^{-i\phi_{\alpha}}&-\mathcal{H}^{T}_{\mathrm{kin},\alpha}(-\mathbf{k})\end{array}\right),$
(48)

with $\alpha=L(R)$ corresponding to the left (right) side. Here,
$\mathcal{H}_{\mathrm{kin},\alpha}(\mathbf{k})$ is the band Hamiltonian,
$\Delta_{\alpha}(\mathbf{k})$ is the pairing matrix, and $\phi_{\alpha}$ is
the overall $\mathrm{U}(1)$ superconducting phase, taken to be uniform. The
anomalous Green’s functions can generally be expressed as [45]

$\mathcal{F}_{\alpha}(\mathbf{k},ip_{n})=-\mathcal{G}_{\mathrm{kin},\alpha}(\mathbf{k},ip_{n})\Delta_{\alpha}(\mathbf{k})e^{i\phi_{\alpha}}\left[1+\mathcal{G}^{\mathrm{T}}_{\mathrm{kin},\alpha}(-\mathbf{k},-ip_{n})\Delta^{\dagger}_{\alpha}(\mathbf{k})\mathcal{G}_{\mathrm{kin},\alpha}(\mathbf{k},ip_{n})\Delta_{\alpha}(\mathbf{k})\right]^{-1}\mathcal{G}_{\mathrm{kin},\alpha}^{\mathrm{T}}(-\mathbf{k},-ip_{n})$
(49)

in which
$\mathcal{G}_{\mathrm{kin},\alpha}(\mathbf{k},ip_{n})=[ip_{n}-\mathcal{H}_{\mathrm{kin}}(\mathbf{k})]^{-1}$
is the single-particle Green’s function of the band Hamiltonian in the absence
of superconductivity. We work in the weak-coupling regime and consider one of
the two following scenarios, under which the anomalous Green’s functions
greatly simplify.

In the first scenario, the band Hamiltonian is spin-independent and thus the
Green’s functions $\mathcal{G}_{\mathrm{kin},\alpha}$ are diagonal in spin
space. For simplicity, we consider the case in which the Hamiltonian is in the
spin-1/2 representation; though, generalization to systems with other internal
degrees of freedom is straightforward. Consider the pairing matrix given by
$\Delta_{\alpha}(\mathbf{k})=(d_{\alpha,0}(\mathbf{k})\sigma_{0}+\mathbf{d}_{\alpha}(\mathbf{k})\cdot\bm{\sigma})i\sigma_{y},$
in which $d_{\alpha,0}$ is the spin singlet pairing and $\mathbf{d}_{\alpha}$
is the $d$-vector corresponding to spin triplet pairing. In this scenario, the
anomalous Green’s functions simplify as

$\mathcal{F}_{\alpha}(\mathbf{k},ip_{n})=\frac{-\Big{[}f_{\alpha,0}(\mathbf{k},ip_{n})+\mathbf{f}_{\alpha}(\mathbf{k},ip_{n})\cdot\bm{\sigma}\Big{]}(i\sigma_{y})}{\Big{[}p_{n}^{2}+\xi_{\alpha,\mathbf{k}}^{2}+\frac{1}{2}\mathrm{Tr}(\Delta_{\alpha}^{\dagger}\Delta_{\alpha})\Big{]}^{2}-\Big{[}(\mathbf{d}^{*}_{\alpha}\cdot\mathbf{d}^{*}_{\alpha})^{2}(d_{\alpha,0}^{2})+2|\mathbf{d}_{\alpha}|^{2}|d_{\alpha,0}|^{2}+(d^{*}_{\alpha,0})^{2}(\mathbf{d}_{\alpha}\cdot\mathbf{d}_{\alpha})+|\mathbf{d}_{\alpha}\times\mathbf{d}_{\alpha}^{*}|^{2}\Big{]}},$
(50)

in which

$\displaystyle f_{\alpha,0}(\mathbf{k},ip_{n})$ $\displaystyle\equiv
d_{\alpha,0}\Big{[}p_{n}^{2}+\xi_{\alpha,\mathbf{k}}^{2}+\frac{1}{2}\mathrm{Tr}(\Delta_{\alpha}^{\dagger}\Delta_{\alpha})\Big{]}-\Big{[}|\mathbf{d}|^{2}d_{\alpha,0}+(\mathbf{d}_{\alpha}\cdot\mathbf{d}_{\alpha})d_{\alpha,0}^{*}\Big{]},$
(51a) $\displaystyle\mathbf{f}_{\alpha}(\mathbf{k},ip_{n})$
$\displaystyle\equiv\Big{[}p_{n}^{2}+\xi_{\alpha,\mathbf{k}}^{2}+\frac{1}{2}\mathrm{Tr}(\Delta_{\alpha}^{\dagger}\Delta_{\alpha})\Big{]}\mathbf{d}_{\alpha}-\Big{[}d_{\alpha,0}^{2}\mathbf{d}_{\alpha}^{*}+|d_{\alpha,0}|^{2}\mathbf{d}_{\alpha}+\mathbf{d}_{\alpha}\times(\mathbf{d}_{\alpha}\times\mathbf{d}_{\alpha}^{*})\Big{]}$
(51b)

include corrections from nonunitary pairing orders. Above,
$\xi_{\alpha,\mathbf{k}}=\xi_{\alpha,-\mathbf{k}}$ is the band energy
dispersion at the left ($\alpha=L$) and right ($\alpha=R$) side of the
junction. $d_{\alpha,0}$ and $\mathbf{d}_{\alpha}$ are both functions of
$\mathbf{k}$, which we have dropped for convenience of notation. When the
pairing is unitary
($\Delta_{\alpha}^{\dagger}\Delta_{\alpha}\propto\sigma_{0}$), the anomalous
Green’s functions take a simpler form,

$\mathcal{F}_{\alpha}(\mathbf{k},ip_{n})=-\frac{e^{i\phi_{\alpha}}\Delta_{\alpha}(\mathbf{k})}{p_{n}^{2}+E^{2}_{\alpha,\mathbf{k}}},$
(52)

in which
$E_{\alpha,\mathbf{k}}=\sqrt{\xi_{\alpha,\mathbf{k}}^{2}+\mathrm{Tr}[\Delta_{\alpha}^{\dagger}(\mathbf{k})\Delta_{\alpha}(\mathbf{k})]/2}$
is the unperturbed dispersion of the BdG quasiparticle.

In the second scenario, there are spin-dependent terms in the band
Hamiltonian. Working in the weak-coupling regime, we transform to a band-
diagonal representation and then project to the states at the Fermi surface
(see Appendix B). Then, the anomalous Green’s functions take an analogous
form,

$\displaystyle\mathcal{F}_{\alpha}^{(bp)}(\mathbf{k},ip_{n})$
$\displaystyle=-\mathcal{G}^{(bp)}_{\mathrm{kin},\alpha}(\mathbf{k},ip_{n})\Delta^{(bp)}_{\alpha}(\mathbf{k})e^{i\phi_{\alpha}}\left[1+\mathcal{G}^{(b)\mathrm{T}}_{\mathrm{kin},\alpha}(-\mathbf{k},-ip_{n})\Delta^{(b)\dagger}_{\alpha}(\mathbf{k})\mathcal{G}^{(bp)}_{\mathrm{kin},\alpha}(\mathbf{k},ip_{n})\Delta^{(bp)}_{\alpha}(\mathbf{k})\right]^{-1}\mathcal{G}^{(b)\mathrm{T}}_{\mathrm{kin},\alpha}(-\mathbf{k},-ip_{n})$
$\displaystyle=-\frac{e^{i\phi_{\alpha}}\Delta^{(bp)}_{\alpha}(\mathbf{k})}{p_{n}^{2}+E^{(bp)2}_{\alpha,\mathbf{k}}}.$
(53)

Here,
$\mathcal{G}^{(bp)}_{\mathrm{kin,\alpha}}(\mathbf{k},ip_{n})=P_{+}^{(b)}(ip_{n}-\Lambda(\mathbf{k}))^{-1}P_{+}^{(b)}$
is the Green’s function for the band Hamiltonian evaluated in the band-
diagonal representation and projected to the Fermi level, with $P_{+}^{(b)}$
being the projection operator and $\Lambda(\mathbf{k})$ the diagonal matrix of
eigenvalues. $\Delta_{\alpha}^{(bp)}(\mathbf{k})$ is the projected effective
pairing matrix in the band-diagonal representation, as shown in Eq. (61) of
Appendix B, and
$E^{(bp)}_{\alpha,\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{(bp)2}+\mathrm{Tr}[\Delta^{(bp)\dagger}_{\alpha}(\mathbf{k})\Delta^{(bp)}_{\alpha}(\mathbf{k})]}$
is the energy of the BdG quasiparticle after projection, with
$\xi^{(bp)}_{\alpha,\mathbf{k}}=\xi^{(bp)}_{\alpha,-\mathbf{k}}$ being the
band dispersion of the eigenstate at the Fermi level. The above form is
justified in the weak-coupling regime, as pairing only occurs between states
near the Fermi surface.

Under the outlined scenarios, the correlation function can simply be expressed
as

$\Phi(i\omega,\phi)=-e^{-i\phi}\sum_{\mathbf{k},\mathbf{k}^{\prime}}w(E_{L,\mathbf{k}^{\prime}},E_{R,\mathbf{k}};\beta,i\omega)\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime}),$
(54)

as shown in Eq. (6). As described in Sec. II, the above correlation function
can be considered in three distinct parts as follows:

(I) The first term, $e^{-i\phi}$, includes superconducting phase difference,
$\phi=\phi_{R}-\phi_{L}$, and corresponds to the first order $2\pi$-periodic
Josephson current. For higher-order contributions to the Josephson current, it
is necessary to consider terms arising from multiple scatterings at the
junction interface.

(II) The second term is the function,
$w(\mathbf{k},\mathbf{k}^{\prime};\beta,i\omega)$, which encodes the density
of states and thermal weighting. After performing the Matsubara frequency
summation over fermion frequencies $ip_{n}$, the function is given by

$\displaystyle
w(E_{L,\mathbf{k}^{\prime}},E_{R,\mathbf{k}};\beta,i\omega)=\frac{1}{4E_{R,\mathbf{k}}E_{L,\mathbf{k}^{\prime}}}$
(55)
$\displaystyle\times\Big{[}\Big{(}n_{F}(E_{L,\mathbf{k}^{\prime}})-n_{F}(E_{R,\mathbf{k}})\Big{)}\left(\frac{1}{E_{L,\mathbf{k}^{\prime}}-E_{R,\mathbf{k}}-i\omega}+\frac{1}{E_{L,\mathbf{k}^{\prime}}-E_{R,\mathbf{k}}+i\omega}\right)$
$\displaystyle\hskip
10.00002pt+\Big{(}1-n_{F}(E_{L,\mathbf{k}^{\prime}})-n_{F}(E_{R,\mathbf{k}})\Big{)}\left(\frac{1}{E_{L,\mathbf{k}^{\prime}}+E_{R,\mathbf{k}}+i\omega}+\frac{1}{E_{L,\mathbf{k}^{\prime}}+E_{R,\mathbf{k}}-i\omega}\right)\Big{]},$

in which $n_{F}$ is the Fermi-Dirac distribution. The final expression for the
retarded correlation function is obtained via analytic continuation,
$i\omega\rightarrow eV+i\delta$ as $\delta\rightarrow 0^{+}$.

(III) The last and most essential term is the form factor,
$\mathfrak{F}(\mathbf{k},\mathbf{k}^{\prime})=\mathrm{Tr}[{\Delta_{L}(\mathbf{k}^{\prime})}{T^{\mathrm{T}}(-\mathbf{k},-\mathbf{k}^{\prime})}{\Delta_{R}^{*}(\mathbf{k})}{T(\mathbf{k},\mathbf{k}^{\prime})}],$
which encodes the symmetry of the order parameters in addition to the