# Flat-band Majorana bound states in topological Josephson junctions

Daisuke Oshima1, Satoshi Ikegaya1, Andreas P. Schnyder2, and Yukio Tanaka1
1Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan
2Max-Planck-Institut für Festkörperforschung, Heisenbergstrasse 1, D-70569
Stuttgart, Germany

###### Abstract

Nodal topological superconductors characterized by $p_{x}$-wave pairing
symmetry host flat-band Majorana bound states causing drastic anomalies in
low-energy electromagnetic responses. Nevertheless, the study of flat-band
Majorana bound states has been at a standstill owing to a serious lack of
candidate materials for $p_{x}$-wave superconductors. In this paper, by
expanding a scheme of planar topological Josephson junctions, we propose a
promising device realizing an effective $p_{x}$-wave superconductor.
Specifically, we consider a three-dimensional Josephson junction consisting of
a thin-film semiconductor hosting a persistent spin-helix state and two
conventional $s$-wave superconductors. We analytically obtain a topological
phase diagram and numerically demonstrate the emergence of flat-band Majorana
bound states by calculating the local density of states.

##  I Introduction

Majorana bound states (MBSs) in topological superconductors (SCs) have been a
central research topic in the field of superconductivity kane_10 ; zhang_11 ;
tanaka_12r ; sato_17 . An important finding of the intensive research efforts
over the past decade is that MBSs have various forms depending on the pairing
symmetries of the parent topological SCs. In two dimensions, for instance, a
chiral $p$-wave SC with broken time-reversal symmetry (TRS) hosts a MBS moving
along the edge of the sample in one direction (i.e., chiral MBS) volovik_97 ;
green_00 ; furusaki_01 , a helical $p$-wave SC preserving TRS harbors a pair
of MBSs moving in opposite directions (i.e., helical MBSs) maiti_06 ;
schnyder_08 ; zhang_09 ; tanaka_09 , and a nodal $p_{x}$-wave SC hosts flat-
band MBSs buchholtz_81 ; nagai_86 ; jerome_80 ; sarma_01 ; tanuma_01 ; sato_11
, will be the focus of this Letter. A remarkable feature of the flat-band MBSs
is that they have an extensive degeneracy at the Fermi level, unlike other
forms of MBSs. Accordingly, the flat-band MBSs of the $p_{x}$-wave SC affect
low-energy electromagnetic responses drastically and cause various anomalies:
a prominent zero-bias conductance peak in normal-metal/SC junctions bruder_90
; hu_94 ; tanaka_95 ; tanaka_96 ; asano_04 ; tanaka_04 ; tanaka_05(1) ;
asano_07 ; ikegaya_15 ; ikegaya_16(1) , a fractional current-phase
relationship in two-dimensional Josephson junctions tanaka_96-(2) ; barash_96
; kwon_04 ; asano_06(1) ; asano_06(2) ; ikegaya_16(2) , and a distinct
paramagnetic Meissner response to magnetic fields higashitani_97 ; suzuki_14 ;
suzuki_15 .

A crucial problem in the research of flat-band MBSs is a serious lack of
candidate materials for $p_{x}$-wave SCs. To resolve this problem, several
theoretical models exhibiting effective $p_{x}$-wave superconductivity have
been proposed, e.g., semiconductor/SC heterostructures alicea_10 ; you_13 ;
ikegaya_21 ; ikegaya_18  and helical $p$-wave SCs under magnetic fields
ikegaya_18 ; law_13 ; rosenow_14 . However, these proposals involve conditions
that are difficult to achieve in experiments, such as the necessity of a high
Zeeman potential exceeding the pair potential of the parent SC and/or the
absence of Rashba-type spin-orbit coupling (SOC) potentials, while these
systems intrinsically have structural inversion asymmetry causing a finite
Rashba SOC ikegaya_15 ; law_13 .

Figure 1: Schematic of a Josephson junction consisting of two conventional
$s$-wave SCs and a thin-film semiconductor hosting a PSH state. A magnetic
field $B$ is applied in the $x$-direction.

In this Letter, to resolve the stalemate in this research field, we propose an
alternative route for creating an effective $p_{x}$-wave superconductor. We
specifically consider a three-dimensional Josephson junction illustrated in
Fig. 1, where a thin-film semiconductor is sandwiched by two conventional
$s$-wave SCs. In addition, we assume that the semiconductor hosts a persistent
spin-helix (PSH) state, which is studied intensively in the field of
spintronics bernevig_06 ; chang_06 ; schliemann_17 ; kohda_17 . This work is
inspired by the innovative studies on planar topological Josephson junctions
(TJJs), which is attracting much attention recently flensberg_17 ; halperin_17
; haim_19 ; setiawan_19 ; ikegaya_20(2) ; shabani_20 ; yin_21 ; nichele_19 ;
yacoby_19 ; shabani_21 . In a planar TJJ, effective topological
superconductivity is realized at the one-dimensional interfacial region of the
two-dimensional Josephson junction. In our setup, we obtain an effective
$p_{x}$-wave superconductor in the vicinity of the two-dimensional
semiconductor segment of the junction. Thanks to the nature of TJJs
flensberg_17 ; halperin_17 , our setup can enter the topologically nontrivial
phase supporting the flat-band MBSs by applying Zeeman potentials smaller than
the pair potentials in the superconducting segment. Moreover, in the vicinity
of the semiconductor segment, the local inversion symmetry along the
$z$-direction is preserved, and thus the Rashba SOC can be negligible.
Therefore, our proposed system is free from the two difficult conditions
required in previous proposals. In addition, it has been recently shown that
the PSH can be realized not only in III–V semiconductor compounds awschalom_09
; salis_12 ; kohda_12 ; salis_14 , but also in two-dimensional ferroelectric
materials such as group-IV monochalcogenide monolayers tsymbal_18 ; ishii_19 ;
ishii_19(2) ; ishii_19(2) ; picozzi_19 ; zhao_19 ; jin_20 ; nardelli_20 ;
chang_20 . Hence, this rapid progress in the field of PSH makes it very likely
that our proposal can be successfully fabricated in experiments. Consequently,
we propose a promising recipe for making an effective $p_{x}$-wave SC
harboring flat-band MBSs.

##  II Effective Hamiltonian and Topological Phase Diagram

We first derive an effective Hamiltonian that describes the essential
properties of the proposed system. Let us start with a simple Bogoliubov–de
Gennes (BdG) Hamiltonian:

$\displaystyle\begin{split}&\check{H}(\bm{r})=\check{H}_{0}(\bm{r})+\check{H}_{\lambda}(\bm{r}),\\\
&\check{H}_{0}(\bm{r})=\left[\begin{array}[]{cccc}\xi(\bm{r})+V\hat{\sigma}_{x}&\Delta(z)(i\hat{\sigma}_{y})\\\
-\Delta^{\ast}(z)(i\hat{\sigma}_{y})&-\xi(\bm{r})\hat{\sigma}_{0}-V\hat{\sigma}_{x}\\\
\end{array}\right],\\\
&\check{H}_{\lambda}(\bm{r})=\left[\begin{array}[]{cc}-i\lambda\delta(z)\partial_{x}\hat{\sigma}_{z}&0\\\
0&-i\lambda\delta(z)\partial_{x}\hat{\sigma}_{z}\\\ \end{array}\right],\\\
&\xi(\bm{r})=-\frac{\hbar^{2}}{2m}(\partial^{2}_{x}+\partial^{2}_{y}+\partial^{2}_{z})-\mu,\\\
&\Delta(z)=\left\\{\begin{array}[]{cl}\Delta e^{i\varphi/2}&\text{for }z\leq
0\\\ \Delta e^{-i\varphi/2}&\text{for }z>0\end{array}\right.,\end{split}$ (1)

where $m$, $\mu$, and $V$ represent the effective electron mass, the chemical
potential, and the Zeeman potential induced by the magnetic field along the
$x$-direction, respectively. The superconducting phase difference is given by
$\varphi$. The Pauli matrices in spin space are denoted by
$\hat{\bm{\sigma}}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})$. The
junction interface is located at $z=0$, and we apply periodic boundary
conditions in the $x$\- and $y$-directions. In this Hamiltonian, we do not
describe explicitly the two-dimensional semiconductor hosting a PSH state.
Instead, we introduce $H_{\lambda}$ describing a unidirectional SOC potential
generating the PSH bernevig_06 ; chang_06 ; schliemann_17 ; kohda_17  only at
$z=0$. The SOC potential as in $H_{\lambda}$ is realized in zinc-blende III–V
semiconductor quantum wells grown along the $[110]$ direction bernevig_06 ;
chang_06 ; salis_14 , in quantum wells with equal amplitudes of Rashba and
Dresselhaus SOC bernevig_06 ; chang_06 ; awschalom_09 ; salis_12 ; kohda_12 ,
and in ferroelectric thin-film materials tsymbal_18 ; ishii_19 ; ishii_19(2) ;
picozzi_19 ; zhao_19 ; jin_20 ; nardelli_20 ; chang_20 . In what follows, we
construct an effective low-energy Hamiltonian that focuses on the subgap
states localized in the vicinity of the junction interface. To simplify the
discussion, we treat $H_{\lambda}$ as a perturbation. Within zeroth-order
perturbation, we first solve the BdG equation for the subgap states,

$\displaystyle\begin{split}&\int
d\bm{r}\check{U}^{\dagger}_{\mathrm{SG}}(\bm{r})\check{H}_{0}(\bm{r})\check{U}_{\mathrm{SG}}(\bm{r})=\check{E}_{0},\\\
&\check{E}_{0}=\mathrm{diag}[\varepsilon-V,\varepsilon+V,-\varepsilon+V,-\varepsilon-V],\\\
&\check{U}_{\mathrm{SG}}(\bm{r})=\check{U}_{\mathrm{SG}}(z)\frac{e^{ik_{x}x}e^{ik_{y}y}}{\sqrt{L_{x}L_{y}}},\end{split}$
(2)

where the conditions,

$\displaystyle
0<\varepsilon<\Delta,\quad\check{U}_{\mathrm{SG}}(z=\pm\infty)=0,$ (3)

are satisfied. Here $k_{x(y)}$ and $L_{x(y)}$ represent the momentum and width
in the $x$-direction ($y$-direction), respectively. The explicit form of
$\check{U}_{\mathrm{SG}}(\bm{r})$ is given in the Supplemental Material (SM)
sup_mat . Then, on the basis of the first-order perturbation theory, we
construct an effective Hamiltonian for the subgap states,

$\displaystyle\begin{split}\check{H}_{\mathrm{eff}}&=\int
d\bm{r}\check{U}^{\dagger}_{\mathrm{SG}}(\bm{r})\check{H}(\bm{r})\check{U}_{\mathrm{SG}}(\bm{r})\\\
&=\left[\begin{array}[]{cccc}\varepsilon-V&\tilde{\lambda}k_{x}&\tilde{\Delta}k_{x}&0\\\
\tilde{\lambda}k_{x}&\varepsilon+V&0&-\tilde{\Delta}k_{x}\\\
\tilde{\Delta}k_{x}&0&-\varepsilon+V&\tilde{\lambda}k_{x}\\\
0&-\tilde{\Delta}k_{x}&\tilde{\lambda}k_{x}&-\varepsilon-V\\\
\end{array}\right],\end{split}$ (4)

with $\varepsilon=\Delta\sqrt{1-D\sin^{2}(\varphi/2)}$,
$\tilde{\lambda}=\lambda\sqrt{D(1-D)}F$, and
$\tilde{\Delta}=\lambda\cot(\varphi/2)F$, where $D=D(k_{x},k_{y},\varphi)$ and
$F=F(k_{x},k_{y},\varphi)$ are the even functions of $k_{x}$, $k_{y}$, and
$\varphi$. The explicit forms of $D(k_{x},k_{y},\varphi)$ and
$F(k_{x},k_{y},\varphi)$ are given in the SM sup_mat . The effective
Hamiltonian of $\check{H}_{\mathrm{eff}}$ is equivalent to the BdG Hamiltonian
of a two-dimensional $p_{x}$-wave SC with SOC and Zeeman potential. Namely, we
can reinterpret $\varepsilon$ as the kinetic energy of the quasi-particle
states confined in the two-dimensional space (i.e., the subgap states
localized at the junction interface), $\tilde{\lambda}$ as the effective SOC
potential, and most importantly, $\tilde{\Delta}$ as the effective
$p_{x}$-wave pair potential acting on the quasi-particle states. Consequently,
we can expect that the present junction hosts flat-band MBSs at the surfaces
perpendicular to the $x$-direction.

To clarify the emergence of the flat-band MBSs, we analyze the topological
property of $\check{H}_{\mathrm{eff}}$. By diagonalizing
$\check{H}_{\mathrm{eff}}$, we obtain the energy eigenvalues as

$\displaystyle
E_{s}=\pm\sqrt{\left(\varepsilon+s\sqrt{V^{2}+\tilde{\lambda}^{2}k_{x}^{2}}\right)^{2}+\tilde{\Delta}^{2}k_{x}^{2}},$
(5)

for $s=\pm$. As discussed in the SM sup_mat  in detail, we obtain the two gap
nodes when the condition $X<V^{2}<\Delta^{2}$ is satisfied, where the lower
bound $X$ is given explicitly in the SM and satisfies

$\displaystyle\lim_{\Delta/\mu\rightarrow 0}X=\Delta^{2}\cos^{2}(\varphi/2).$
(6)

The two gap nodes are located at $(k_{x},k_{y})=(0,\pm k_{F}\sqrt{1+\zeta})$
with

$\displaystyle\zeta=\frac{1}{\mu}\left[V^{2}-\frac{\Delta^{2}}{2}\left\\{1+\cos^{2}\frac{\varphi}{2}\right\\}\right]\sqrt{\frac{1}{V^{2}-\Delta^{2}\cos^{2}\frac{\varphi}{2}}},$

where $k_{F}=\sqrt{2m\mu}/\hbar$. Moreover, the energy spectrum becomes
completely gapless at $\varphi=\pi$, because the effective pair potential
$\tilde{\Delta}\propto\cot(\phi/2)$ vanishes. Consequently, we obtain two
distinct nodal superconducting phases as summarized in Fig. 2(a): the nodal
phase I for $\varphi<\pi$ and the nodal phase II for $\varphi>\pi$.

Next, we discuss the symmetry properties of the present junction. The
effective Hamiltonian $\check{H}_{\mathrm{eff}}$ has TRS and particle-hole
symmetry (PHS) as,

$\displaystyle\check{T}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}(\bm{k})\check{T}_{\mathrm{eff}}^{-1}=\check{H}_{\mathrm{eff}}(-\bm{k}),\quad\check{T}_{\mathrm{eff}}=-i\hat{\sigma}_{z}\check{\tau}_{z}\mathcal{K},$
(7)
$\displaystyle\check{C}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}(\bm{k})\check{C}_{\mathrm{eff}}^{-1}=-\check{H}_{\mathrm{eff}}(-\bm{k}),\quad\check{C}_{\mathrm{eff}}=\check{\tau}_{x}\mathcal{K},$
(8)

where the Pauli matrices in Nambu space are given by
$\check{\bm{\tau}}=(\check{\tau}_{x},\check{\tau}_{y},\check{\tau}_{z})$, and
$\mathcal{K}$ is the complex-conjugation operator. By combining
$\check{T}_{\mathrm{eff}}$ and $\check{C}_{\mathrm{eff}}$, we obtain a chiral
symmetry (CS) as

$\displaystyle\check{S}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}(\bm{k})\check{S}_{\mathrm{eff}}^{-1}=-\check{H}_{\mathrm{eff}}(\bm{k}),\quad\check{S}_{\mathrm{eff}}=\check{T}_{\mathrm{eff}}\check{C}_{\mathrm{eff}}.$
(9)

The original Hamiltonian $\check{H}(\bm{r})$ also preserves TRS, PHS, and CS
as

$\displaystyle\check{T}\check{H}(\bm{r})\check{T}^{-1}=\check{H}(\bm{r}),\quad\check{T}=\check{M}_{xy}\check{T}_{-},$
(10)
$\displaystyle\check{C}\check{H}(\bm{r})\check{C}^{-1}=-\check{H}(\bm{r}),\quad\check{C}=\check{\tau}_{x}\mathcal{K},$
(11)
$\displaystyle\check{S}\check{H}(\bm{r})\check{S}^{-1}=-\check{H}(\bm{r}),\quad\check{S}=\check{T}\check{C},$
(12)

where $\check{M}_{xy}=-i\hat{\sigma}_{z}\check{\tau}_{z}\mathcal{R}_{z}$
denotes mirror reflection symmetry with respect to the $xy$-plane, and
$\check{T}_{-}=i\hat{\sigma}_{y}\mathcal{K}$ represents conventional TRS
obeying $\check{T}_{-}^{2}=-1$; $\mathcal{R}_{z}$ describes the reflection of
the spatial coordinate $z$ (i.e., $z\rightarrow-z$). Furthermore, the present
system preserves the inversion symmetry with respect to the $y$-axis:

$\displaystyle\check{R}_{k_{y}}\check{H}_{\mathrm{eff}}(\bm{k})\check{R}_{k_{y}}^{-1}=\check{H}_{\mathrm{eff}}(k_{x},-k_{y})=\check{H}_{\mathrm{eff}}(\bm{k}),$
(13)
$\displaystyle\check{R}_{y}\check{H}(\bm{r})\check{R}_{y}^{-1}=\check{H}(x,-y,z)=\check{H}(\bm{r}),$
(14)

where $\check{R}_{k_{y}}$ ($\check{R}_{y}$) inverts the sign of $k_{y}$ ($y$)
as $k_{y}\rightarrow-k_{y}$ ($y\rightarrow-y$). As a result, we additionally
obtain a modified TRS and a modified PHS as

$\displaystyle\begin{split}&\check{\underline{T}}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}(\bm{k})\check{\underline{T}}_{\mathrm{eff}}^{-1}=\check{H}_{\mathrm{eff}}(-k_{x},k_{y}),\\\
&\check{\underline{T}}_{\mathrm{eff}}=\check{R}_{k_{y}}\check{T}_{\mathrm{eff}},\end{split}$
(15)

and

$\displaystyle\begin{split}&\check{\underline{C}}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}(\bm{k})\check{\underline{C}}_{\mathrm{eff}}^{-1}=-\check{H}_{\mathrm{eff}}(-k_{x},k_{y}),\\\
&\check{\underline{C}}_{\mathrm{eff}}=\check{R}_{k_{y}}\check{C}_{\mathrm{eff}},\end{split}$
(16)

where these modified symmetries do not involve the sign change in $k_{y}$. The
original BdG Hamiltonian also preserves

$\displaystyle\check{\underline{T}}\check{H}(\bm{r})\check{\underline{T}}^{-1}=\check{H}(\bm{r}),\quad\check{\underline{T}}=\check{R}_{y}\check{T},$
(17)
$\displaystyle\check{\underline{C}}\check{H}(\bm{r})\check{\underline{C}}^{-1}=-\check{H}(\bm{r}),\quad\check{\underline{C}}=\check{R}_{y}\check{C}.$
(18)

The relevant CS operator satisfies

$\displaystyle\check{\underline{S}}_{\mathrm{eff}}=\check{\underline{T}}_{\mathrm{eff}}\check{\underline{C}}_{\mathrm{eff}}=\check{S}_{\mathrm{eff}},$
(19)
$\displaystyle\check{\underline{S}}=\check{\underline{T}}\;\check{\underline{C}}=\check{S}.$
(20)

To characterize the nodal phases of $\check{H}_{\mathrm{eff}}(k_{x},k_{y})$
topologically, we set $k_{y}$ to a fixed value, and study
$\check{H}_{\mathrm{eff}}$ as a function of only $k_{x}$. I.e, we interpret
$\check{H}_{\mathrm{eff}}(k_{x},k_{y})$ as a one parameter family of one-
dimensional Hamiltonians with momentum $k_{x}$ and parameter $k_{y}$ sato_11 .
Since $\check{H}_{\mathrm{eff}}(k_{x},k_{y})$ satisfies Eq. (15) and Eq. (16)
with $\check{\underline{T}}_{\mathrm{eff}}^{2}=+1$ and
$\check{\underline{C}}_{\mathrm{eff}}^{2}=+1$, we can classify
$\check{H}_{\mathrm{eff}}(k_{x},k_{y})$ into the BDI symmetry class in one-
dimension schnyder_08 . Therefore, the topological property of
$\check{H}_{\mathrm{eff}}(k_{x},k_{y})$ can be characterized by a one-
dimensional winding number schnyder_08 ; sato_11 ,

$\displaystyle w(k_{y})=\frac{i}{4\pi}\int
dk_{x}\mathrm{Tr}[\check{S}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}^{-1}\partial_{k_{x}}\check{H}_{\mathrm{eff}}].$
(21)

By computing this as in the SM sup_mat , we obtain the winding number for the
nodal phase I (phase II) as

$\displaystyle
w(k_{y})=\left\\{\begin{array}[]{cl}+(-)1&\text{for}\quad|k_{y}|<k_{F}\sqrt{1+\zeta}\\\
0&\text{otherwise}\end{array}\right.,$ (24)

where $w(k_{y})=0$ irrespective of $k_{y}$ for the fully gapped phase. Since
the winding number is nonzero in a finite range of $k_{y}$, according to the
bulk-boundary correspondence, we obtain flat-band MBSs at a surface
perpendicular to the $x$-direction sato_11 . Unless the gap of
$\tilde{\Delta}$ closes or the CS of $\check{S}_{\mathrm{eff}}$ (and therefore
$\check{S}$) is broken, the flat-band MBSs can exist stably due to the
topological protection. Therefore, we can expect that higher-order
perturbations of $\lambda$ do not disturb the flat-band MBSs.

Moreover, we can additionally define a $\mathbb{Z}_{2}$ topological number as
$\mathbb{Z}_{2}(k_{y})=(-1)^{w(k_{y})}$ halperin_17 . In terms of the
$\mathbb{Z}_{2}$ topological number, the flat-band MBSs in both nodal phases I
and II are characterized by $\mathbb{Z}_{2}(k_{y})=-1$. When some perturbation
breaks the modified TRS and the CS, the symmetry class of
$\check{H}_{\mathrm{eff}}(k_{x},k_{y})$ changes into class D schnyder_08 .
Nevertheless, as long as the modified PHS of
$\check{\underline{C}}_{\mathrm{eff}}$ (and therefore $\check{\underline{C}}$)
is preserved, the flat-band MBSs still exist, but are now characterized by the
$\mathbb{Z}_{2}$ topological number schnyder_08 ; halperin_17 , rather than
the winding number. Hence, perturbations breaking the modified TRS might alter
the topological phase boundaries, but do not remove the flat-band MBSs. When
some perturbation breaks the modified PHS and the modified TS, the symmetry
class of $\check{H}_{\mathrm{eff}}(k_{x},k_{y})$ changes into class AIII
schnyder_08 . In this case, the system still preserves CS of
$\check{S}_{\mathrm{eff}}$ (and therefore $\check{S}$), leading to flat-band
MBSs characterized by the winding number $w(k_{y})$. In summary, we can expect
the emergence of flat-band MBSs as long as the present junction preserves the
CS or the modified PHS.

Figure 2: (a) Phase diagram of the nodal structure. The nodal phases satisfy
the condition $\Delta^{2}\cos(\varphi/2)<V^{2}<\Delta^{2}$, which is obtained
in the limit of $\mu\gg\Delta$. The nodal phase I and II are characterized by
the positive and negative winding number, respectively. Simultaneously, the
both nodal phases are characterized by the nontrivial $\mathbb{Z}_{2}$
topological number. (b) Zero-energy LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as a
function of the Zeeman potential $V$ and the superconducting phase difference
$\varphi$. In (c) and (d), we show $\rho_{\rm edge}(k_{y};E)$ on a log-scale,
i.e., $\log_{10}\rho_{\rm edge}(k_{y};E)$, as a function of the energy $E$ and
the momentum $k_{y}$. We plot the result for (c) the topologically nontrivial
phase with $(V,\varphi)=(0.7\Delta,0.65\pi)$ and that for (d) the
topologically trivial phase with $(V,\varphi)=(0.25\Delta,0.25\pi)$,
respectively.

##  III Flat-band Majorana bound states

To check the emergence of the flat-band MBSs, we describe the present junction
using a tight-binding model and calculate numerically the local density of
states (LDOS) at the surface perpendicular to the $x$-direction. We denote the
lattice constant by $a_{0}$ and parameterize the lattice sites by
$\bm{r}=i_{x}\bm{x}+i_{y}\bm{y}+i_{z}\bm{z}$, with
$|\bm{x}|=|\bm{y}|=|\bm{z}|=a_{0}$. The upper (lower) SC occupies $1\leq
i_{z}\leq L_{S}$ ($-L_{S}\leq i_{z}\leq-1$), and the two-dimensional
semiconductor is located at $i_{z}=0$ (see Fig. 1). We apply a periodic (an
open) boundary condition in the $y$-direction ($z$-direction). We assume that
the system is semi-infinite in the $x$-direction, i.e., $0\leq
i_{x}\leq\infty$. For latter convenience, we apply a Fourier transformation on
the annihilation operator at $\bm{r}$ with spin $\sigma$ as
$c_{\bm{r},\sigma}=\sum_{k_{y}}e^{ik_{y}a_{0}i_{y}}c_{\tilde{\bm{r}},k_{y},\sigma}/\sqrt{L_{y}}$,
where $\tilde{\bm{r}}=i_{x}\bm{x}+i_{z}\bm{z}$. Furthermore, in what follows,
we omit the subscript of $k_{y}$ in the annihilation operator for notational
simplicity, i.e., $c_{\tilde{\bm{r}},\sigma}$ means
$c_{\tilde{\bm{r}},k_{y},\sigma}$. The BdG tight-binding Hamiltonian reads
$H=\sum_{k_{y}}H(k_{y})$ with

$\displaystyle H(k_{y})=$
$\displaystyle-t\sum_{\tilde{\bm{r}}}\sum_{\sigma=\uparrow,\downarrow}\left(c^{\dagger}_{\tilde{\bm{r}}+\bm{x},\sigma}c_{\tilde{\bm{r}},\sigma}+c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}}+\bm{x},\sigma}\right)$
$\displaystyle-t\sum_{\tilde{\bm{r}},\sigma}\left(c^{\dagger}_{\tilde{\bm{r}}+\bm{z},\sigma}c_{\tilde{\bm{r}},\sigma}+c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}}+\bm{z},\sigma}\right)$
$\displaystyle+\sum_{\tilde{\bm{r}},\sigma}\xi_{i_{z}}c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}},\sigma}+\sum_{\tilde{\bm{r}},\sigma,\sigma^{\prime}}V(\hat{\sigma}_{y})_{\sigma,\sigma^{\prime}}c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}},\sigma^{\prime}}$
$\displaystyle+\sum_{i_{x}}\sum_{i_{z}=1}^{L_{S}}\Delta
e^{i\varphi/2}(c^{\dagger}_{\tilde{\bm{r}},k_{y},\uparrow}c^{\dagger}_{\tilde{\bm{r}},-k_{y},\uparrow}+\mathrm{H.c.})$
$\displaystyle+\sum_{i_{x}}\sum_{i_{z}=-L_{S}}^{-1}\Delta
e^{-i\varphi/2}(c^{\dagger}_{\tilde{\bm{r}},k_{y},\uparrow}c^{\dagger}_{\tilde{\bm{r}},-k_{y},\uparrow}+\mathrm{H.c.})$
$\displaystyle+\sum_{\tilde{\bm{r}},\sigma,\sigma^{\prime}}\frac{i\lambda}{2}\delta_{i_{z},0}(\hat{\sigma}_{z})_{\sigma,\sigma^{\prime}}$
$\displaystyle\qquad\times\left(c^{\dagger}_{\tilde{\bm{r}}+\bm{x},\sigma}c_{\tilde{\bm{r}},\sigma^{\prime}}-c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}}+\bm{x},\sigma^{\prime}}\right),$
(25) $\displaystyle\xi_{i_{z}}=$ $\displaystyle\left\\{\begin{aligned}
4t-\mu_{N}-2t\cos(k_{y}a_{0})\quad&{\rm for}\quad i_{z}=0\\\
6t-\mu_{S}-2t\cos(k_{y}a_{0})\quad&{\rm otherwise}\end{aligned}\right.,\,$
(26)

where $t$ is the nearest-neighbor hopping integral, and $\mu_{S(N)}$ is the
chemical potential in the SC (semiconductor). From the BdG Hamiltonian in Eq.
(25), we calculate the Green’s function
$G_{k_{y}}(\tilde{\bm{r}},\tilde{\bm{r}}^{\prime};E)$ for each $k_{y}$ by
using recursive Green’s function techniques fisher_81 ; ando_91 . The LDOS is
evaluated by

$\displaystyle\rho_{k_{y}}(\tilde{\bm{r}};E)=-\frac{1}{\pi}\left\\{{\rm
Tr}\left[G_{k_{y}}(\tilde{\bm{r}},\tilde{\bm{r}};E+i\delta)\right]\right\\},$
(27)

where $\delta$ is a small positive value added to the energy $E$. We
specifically focus on the LDOS at the edge of the semiconductor segment (i.e.,
$i_{x}=i_{z}=0$):

$\displaystyle\rho_{\rm edge}(k_{y};E)=\rho_{k_{y}}(i_{x}=0,i_{z}=0;E),$ (28)
$\displaystyle\bar{\rho}_{\rm edge}(E)=\sum_{k_{y}}\rho_{\rm edge}(k_{y};E).$
(29)

Hereafter, we fix the parameters as $t=10\Delta$, $\mu_{S}=30\Delta$,
$\mu_{N}=10\Delta$, $\lambda=10\Delta$, $\delta=10^{-4}\Delta$ and $L_{S}=30$.
Figure 2(b) shows the zero-energy LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as a
function of $V$ and $\varphi$. We see that $\bar{\rho}_{\rm edge}(E=0)$
increases significantly in the parameter regions corresponding to the
topologically nontrivial phases in Fig. 2(a). In Figs. 2(c) and 2(d), we show
$\rho_{\rm edge}(k_{y};E)$ for (c) the topologically nontrivial phase with
$(V,\varphi)=(0.7\Delta,0.65\pi)$ and for (d) the topologically trivial phase
with $(V,\varphi)=(0.25\Delta,0.25\pi)$, respectively. We plot
$\log_{10}\rho_{\rm edge}(k_{y};E)$ instead of the raw data of $\rho_{\rm
edge}(k_{y};E)$. As shown in Fig. 2(c), the LDOS in the topologically
nontrivial phase shows the drastic enhancement at $E=0$ for the broad range of
$k_{y}$, while that in the topologically trivial phase has a gapped structure.
On the basis of both analytical and numerical results, we demonstrate the
emergence of the flat-band MBSs due to the effective $p_{x}$-wave
superconductivity.

Figure 3: The zero-energy LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as a function of
$V$ and $\varphi$. The results are normalized by the same unit as the one in
Fig. 2(b). In (a) and (b), we consider (a) $(t_{1},t_{2})=(-2\Delta,-2\Delta)$
and (b) $(t_{1},t_{2})=(-2\Delta,2\Delta)$, respectively. In (a), both CS in
Eq. (12) and the modified PHS in Eq. (18) are preserved. In (b), the modified
PHS is preserved, while the CS is broken. In (c) and (d), we introduce the
Rashba SOC potential with
$(\lambda^{\prime}_{1},\lambda^{\prime}_{2})=(2.5\Delta,-2.5\Delta)$ and (d)
$(\lambda^{\prime}_{1},\lambda^{\prime}_{2})=(3.5\Delta,-1.5\Delta)$,
respectively. In (c), the CS is preserved, while the modified PHS is broken.
In (d), both CS and modified PHS are broken. The dotted line in (c) indicates
$(V,\varphi)$ taking the local minima of $\bar{\rho}_{\rm edge}(E=0)$.

Finally, we discuss the stability of the flat-band MBSs against perturbations.
We first consider a perturbation of the form,

$\displaystyle H_{\delta t}=-\sum_{\tilde{\bm{r}},\sigma}$
$\displaystyle(t_{1}\delta_{i_{z},1}+t_{2}\delta_{i_{z},-1})$
$\displaystyle\times\left(c^{\dagger}_{\tilde{\bm{r}}+\bm{z},\sigma}c_{\tilde{\bm{r}},\sigma}+c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}}+\bm{z},\sigma}\right),$
(30)

which lowers the transparency at the superconductor/semiconductor interfaces.
When $t_{1}=t_{2}$, the system preserves both CS in Eq. (12) and modified PHS
in Eq. (18). The figure 3(a) shows $\bar{\rho}_{\rm edge}(E=0)$ in the chiral
symmetric case as a function of $V$ and $\varphi$, where we choose
$t_{1}=t_{2}=-2\Delta$. Since the CS and the modified PHS are preserved, the
flat-band MBSs can remain at zero-energy, and therefore we find a significant
enhancement in the LDOS. Even so, as compared with Fig. 2(b), we find that the
LDOS around $(V,\varphi)=(0,\pi)$ is suppressed. Furthermore, by calculating
$\rho_{\rm edge}(k_{y};E)$, we see that the induced gap (which corresponds to
$\tilde{\Delta}$ in the first-order perturbation theory) is suppressed by
nonzero $t_{1(2)}$ (see the SM sup_mat ). The shrinking of the topologically
nontrivial phase and the suppression of the induced gap owing to the reduction
of the junction transparency are also found in the planar TJJs halperin_17 .
When $t_{1}\neq t_{2}$, the system preserves the modified PHS, while the CS is
broken. In Fig. 3(b), we show $\bar{\rho}_{\rm edge}(E=0)$, where we choose
$t_{1}=-2\Delta$ and $t_{2}=2\Delta$. Even in the absence of the CS, since the
modified PHS is preserved, we still find a significant enhancement in the
LDOS. By calculating $\rho_{\rm edge}(k_{y};E)$, we also confirm that the
enhancement is indeed caused by the flat-band MBSs (see also SM sup_mat ). We
note that the flat-band MBSs with $t_{1}\neq t_{2}$ are no longer
characterized by the winding number in Eq. (21), but by the nontrivial
$\mathbb{Z}_{2}$ topological number. We additionally find that the two nodal
phases originally separated by $\varphi=\pi$ seem to merge into a single nodal
phase. Then, we consider another perturbation,

$\displaystyle H_{R}=$
$\displaystyle\sum_{\tilde{\bm{r}},\sigma,\sigma^{\prime}}\left\\{\lambda^{\prime}_{1}\sin(k_{y}a_{0})\delta_{i_{z},1}+\lambda^{\prime}_{2}\sin(k_{y}a_{0})\delta_{i_{z},-1}\right\\}$
$\displaystyle\qquad\times(\hat{\sigma}_{x})_{\sigma,\sigma^{\prime}}c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}},\sigma^{\prime}}$
$\displaystyle-$
$\displaystyle\sum_{\tilde{\bm{r}},\sigma,\sigma^{\prime}}\left\\{\frac{i\lambda^{\prime}_{1}}{2}\delta_{i_{z},1}+\frac{i\lambda^{\prime}_{2}}{2}\delta_{i_{z},-1}\right\\}$
$\displaystyle\qquad\times(\hat{\sigma}_{y})_{\sigma,\sigma^{\prime}}\left(c^{\dagger}_{\tilde{\bm{r}}+\bm{x},\sigma}c_{\tilde{\bm{r}},\sigma^{\prime}}-c^{\dagger}_{\tilde{\bm{r}},\sigma}c_{\tilde{\bm{r}}+\bm{x},\sigma^{\prime}}\right),$
(31)

which describes the Rashba SOC potentials at one layer above and below the
semiconductor segment. When $\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}$, the
junction preserves the CS, while the modified PHS is broken. In Fig. 3(c), we
plot $\bar{\rho}_{\rm edge}(E=0)$, where we choose
$\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}=2.5\Delta$. Since the CS is
preserved, we still find a significant enhancement in the LDOS due to the
flat-band MBSs. We also find that $\bar{\rho}_{\rm edge}(E=0)$ becomes
asymmetric with respect to $\varphi=\pi$, i.e., the topologically nontrivial
region for $\varphi<\pi$ shrinks. By calculating $\rho_{\rm edge}(k_{y};E)$,
we also confirm that the amplitude of the induced gap also becomes asymmetric
with respect to $\varphi=\pi$, i.e., the induced gap for $\varphi<\pi$ is
smaller than that for $\varphi>\pi$ (see SM sup_mat ). Strictly speaking, the
phase boundary between two distinct nodal phases slightly deviates from
$\varphi=\pi$; see dotted line in Fig. 3(c). When
$\lambda^{\prime}_{1}\neq-\lambda^{\prime}_{2}$, both CS and modified PHS are
broken. In Fig. 3(d), we show $\bar{\rho}_{\rm edge}(E=0)$ with
$\lambda^{\prime}_{1}=3.5\Delta$ and $\lambda^{\prime}_{2}=-1.5\Delta$. We see
that the zero-energy LDOS is strongly suppressed, which implies the absence of
the flat-band MBSs. More specifically, by calculating $\rho_{\rm
edge}(k_{y};E)$, we confirm that the Rashba SOC potential causes the finite
slope in the dispersion of flat-band MBSs (see SM sup_mat ). The flat-band
MBSs in other effective $p_{x}$-wave superconductors are also destroyed by the
Rashba SOC potential in a similar manner ikegaya_15 ; law_13 . However, we
emphasize that the present junction intrinsically preserves inversion symmetry
with respect to the $z$-direction, and thus the Rashba SOC potential with
$\lambda^{\prime}_{1}\neq-\lambda^{\prime}_{2}$ becomes negligibly small. In
short, although the perturbations shrink the parameter regions for the
topologically nontrivial phase, we can obtain the flat-band MBSs as long as
the Josephson junction preserves the CS or the modified PHS.

##  IV Discussion

Within our numerical calculations, we obtain a sizable induced gap
corresponding to almost $25\%$ of the parent pair potential $\Delta$, which is
comparable to the induced gap of the planar TJJ halperin_17 . When we employ
Nb ($\Delta\approx 1.5\mathrm{meV}$ sutton_62 ) for the superconducting
segment, for instance, the induced gap is up to $375\mu\mathrm{eV}$. The
induced gap becomes larger with a larger SOC potential sustaining the PSH
state and with higher transparency at the semiconductor/superconductor
interface. This property of the induced gap provides a guide for designing
material combinations in experiments. Although electronic correlations are
weak in the thin-film semiconductor, studying the instability of the flat-band
MBSs against interactions would be an intriguing future work sigrist_00 ;
lee_14 ; schnyder_16 .

In summary, we consider a three-dimensional Josephson junction consisting of a
semiconductor hosting a PSH state and a conventional $s$-wave superconductor
(see Fig. 1). We demonstrate the emergence of the flat-band MBSs due to the
effective $p_{x}$-wave superconductivity in the vicinity of the semiconductor
segment. The proposed system has three significant advantages as follows: (i)
the system can be fabricated by interfacing existing materials, (ii) the
system can enter the topologically nontrivial phase without large Zeeman
potentials exceeding the pair potentials, and (iii) the system has a symmetric
structure leading to a vanishing Rashba SOC potential. Our proposal therefore
represents a promising approach to realize the effective $p_{x}$-wave
superconductor hosting the flat-band MBSs.

###### Acknowledgements.

We are grateful to S. Tamura for the fruitful discussions. S.I. is supported
by a Grant-in-Aid for JSPS Fellows (JSPS KAKENHI Grant No. JP21J00041). Y. T.
acknowledges support from a Grant-in-Aid for Scientific Research B (KAKENHI
Grant No. JP18H01176 and No. JP20H01857) and for Scientific Research A
(KAKENHI Grant No. JP20H00131), and Japan RFBR Bilateral Join Research
Projects/Seminars number 19-52-50026. This work is also supported by the JSPS
Core-to-Core Program (No. JPJSCCA20170002). D.O. and S.I. contributed equally
to this work.

## References

  * (1) M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) Y. Tanaka, M. Sato and N. Nagaosa, Symmetry and Topology in Superconductors —Odd-Frequency Pairing and Edge States—, J. Phys. Soc. Jpn. 81, 011013 (2012).
  * (4) M. Sato and Y. Ando, Topological superconductors: a review, Rep. Prog. Phys. 80, 076501 (2017).
  * (5) G. E. Volovik, On edge states in superconductors with time inversion symmetry breaking, JETP Lett. 66, 522 (1997).
  * (6) N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum Hall effect, Phys .Rev. B 61, 10267 (2000).
  * (7) A. Furusaki, M. Matsumoto, and M. Sigrist, Spontaneous Hall effect in a chiral $p$-wave superconductor, Phys .Rev. B 64, 054514 (2001).
  * (8) K. Sengupta, R. Roy, and M. Maiti, Spin Hall effect in triplet chiral superconductors and graphene, Phys .Rev. B 74, 094505 (2006).
  * (9) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * (10) X.-L. Qi, T. L. Hughes, S. Raghu, and S.-C. Zhang, Time-Reversal-Invariant Topological Superconductors and Superfluids in Two and Three Dimensions, Phys .Rev. Lett. 102, 187001 (2009).
  * (11) Y. Tanaka, T. Yokoyama, A. V. Balatsky, and N. Nagaosa, Theory of topological spin current in noncentrosymmetric superconductors, Phys .Rev. B 79, 060505(R) (2009).
  * (12) L. J. Buchholtz and G. Zwicknagl, Identification of $p$-wave superconductors, Phys. Rev. B 23, 5788 (1981).
  * (13) J. Hara and K. Nagai, A Polar State in a Slab as a Soluble Model of $p$-Wave Fermi Superfluid in Finite Geometry, Prog. Theor. Phys. 76, 1237 (1986).
  * (14) D. Jérome, A. Mazaud, M. Ribault, and K. Bechgaard, Superconductivity in a synthetic organic conductor (TMTSF)2PF6, J. Phys. (France) Lett. 41, L95 (1980).
  * (15) K. Sengupta, I. Žutić, H.-J. Kwon, V. M. Yakovenko, and S. Das Sarma, Midgap edge states and pairing symmetry of quasi-one-dimensional organic superconductors, Phys. Rev. B 63, 144531 (2001).
  * (16) Y. Tanuma, K. Kuroki, Y. Tanaka, and S. Kashiwaya, Theoretical study of quasiparticle states near the surface of a quasi-one-dimensional organic superconductor (TMTSF)2PF6, Phys. Rev. B 64, 214510 (2001).
  * (17) M. Sato, Y. Tanaka, K. Yada, and T. Yokoyama, Topology of Andreev bound states with flat dispersion, Phys. Rev. B 83, 224511 (2011).
  * (18) C. Bruder, Andreev scattering in anisotropic superconductors, Phys .Rev. B 41, 4017 (1990).
  * (19) C.-R. Hu, Midgap surface states as a novel signature for $d_{x_{a}^{2}-x_{b}^{2}}$-wave superconductivity, Phys. Rev. Lett. 72, 1526 (1994).
  * (20) Y. Tanaka and S. Kashiwaya, Theory of Tunneling Spectroscopy of $d$-Wave Superconductors, Phys. Rev. Lett. 74, 3451 (1995).
  * (21) S. Kashiwaya, Y. Tanaka, M. Koyanagi, and K. Kajimura, Theory for tunneling spectroscopy of anisotropic superconductors, Phys .Rev. B 53, 2667 (1996).
  * (22) Y. Asano, Y. Tanaka, and S. Kashiwaya, Phenomenological theory of zero-energy Andreev resonant states, Phys. Rev. B 69, 134501 (2004).
  * (23) Y. Tanaka and S. Kashiwaya, Anomalous charge transport in triplet superconductor junctions, Phys. Rev. B 70, 012507 (2004).
  * (24) Y. Tanaka, S. Kashiwaya, and T. Yokoyama, Theory of enhanced proximity effect by midgap Andreev resonant state in diffusive normal-metal/triplet superconductor junctions, Phys. Rev. B 71, 094513 (2005).
  * (25) Y. Asano, Y. Tanaka, A. A. Golubov, and S. Kashiwaya, Conductance Spectroscopy of Spin-Triplet Superconductors, Phys. Rev. Lett. 99, 067005 (2007).
  * (26) S. Ikegaya, Y. Asano, and Y. Tanaka, Anomalous proximity effect and theoretical design for its realization, Phys. Rev. B 91, 174511 (2015).
  * (27) S. Ikegaya, S.-I. Suzuki, Y. Tanaka, and Y. Asano, Quantization of conductance minimum and index theorem, Phys. Rev. B 94, 054512 (2016).
  * (28) Y. Tanaka and S. Kashiwaya, Theory of the Josephson effect in $d$-wave superconductors, Phys. Rev. B 53, 11957(R) (1996).
  * (29) Y. S. Barash, H. Burkhardt and D. Rainer, Low-Temperature Anomaly in the Josephson Critical Current of Junctions in $d$-Wave Superconductors, Phys. Rev. Lett. 77, 4070 (1996).
  * (30) H. -J. Kwon, K. Sengupta, and V. M. Yakovenko, Fractional ac Josephson effect in $p$\- and $d$-wave superconductors, Eur. Phys. J. B 37, 349$-$361(2004).
  * (31) Y. Asano, Y. Tanaka, and S. Kashiwaya, Anomalous Josephson Effect in $p$-Wave Dirty Junctions, Phys. Rev. Lett. 96, 097007 (2006).
  * (32) Y. Asano, Y. Tanaka, T. Yokoyama, and S. Kashiwaya, Josephson current through superconductor/diffusive-normal-metal/superconductor junctions: Interference effects governed by pairing symmetry, Phys. Rev. B 74, 064507 (2006).
  * (33) S. Ikegaya and Y. Asano, Degeneracy of Majorana bound states and fractional Josephson effect in a dirty SNS junction, J. Phys.: Condens. Matter 28, 375702 (2016).
  * (34) S. Higashitani, Mechanism of Paramagnetic Meissner Effect in High-Temperature Superconductors, J. Phys. Soc. Jpn., 66, 2556 (1997).
  * (35) S.-I. Suzuki and Y. Asano, Paramagnetic instability of small topological superconductors, Phys. Rev. B 89, 184508 (2014).
  * (36) S.-I. Suzuki and Y. Asano, Effects of surface roughness on the paramagnetic response of small unconventional superconductors, Phys. Rev. B 91, 214510 (2015).
  * (37) J. Alicea, Majorana fermions in a tunable semiconductor device, Phys. Rev. B 81, 125318 (2010).
  * (38) J. You, C. H. Oh, and V. Vedral, Majorana fermions in $s$-wave noncentrosymmetric superconductor with Dresselhaus (110) spin-orbit coupling, Phys. Rev. B 87, 054501 (2013).
  * (39) S. Ikegaya, J. Lee, A. P. Schnyder, and Y. Asano, Strong anomalous proximity effect from spin-singlet superconductors, Phys. Rev. B 104, L020502 (2021).
  * (40) S. Ikegaya, S. Kobayashi, and Y. Asano, Symmetry conditions of a nodal superconductor for generating robust flat-band Andreev bound states at its dirty surface, Phys. Rev. B 97, 174501 (2018).
  * (41) C. L. M. Wong, J. Liu, K. T. Law, P. A. Lee, Majorana flat bands and unidirectional Majorana edge states in gapless topological superconductors, Phys. Rev. B 88, 060504(R) (2013).
  * (42) T. Hyart, A. R. Wright, and B. Rosenow, Zeeman-field-induced topological phase transitions in triplet superconductors, Phys. Rev. B 90, 064507 (2014).
  * (43) B. A. Bernevig, J. Orenstein, and S.-C. Zhang, Exact SU(2) Symmetry and Persistent Spin Helix in a Spin-Orbit Coupled System, Phys. Rev. Lett. 97, 236601 (2006).
  * (44) M.-H. Liu, K.-W. Chen, S.-H. Chen, and C.-R. Chang, Persistent spin helix in Rashba-Dresselhaus two-dimensional electron systems, Phys. Rev. B 74, 235322 (2006).
  * (45) J. Schliemann, Colloquium: Persistent spin textures in semiconductor nanostructures, Rev. Mod. Phys. 89, 011001(2017).
  * (46) M. Kohda and G. Salis, Physics and application of persistent spin helix state in semiconductor heterostructures, Semicond. Sci. Technol. 32, 073002 (2017).
  * (47) M. Hell, M. Leijnse, and K. Flensberg, Two-Dimensional Platform for Networks of Majorana Bound States, Phys. Rev. Lett. 118, 107701 (2017).
  * (48) F. Pientka, A. Keselman, E. Berg, A. Yacoby, A. Stern, and B. I. Halperin, Topological Superconductivity in a Planar Josephson Junction, Phys. Rev. X 7, 021032 (2017).
  * (49) A. Haim and A. Stern, Benefits of Weak Disorder in One-Dimensional Topological Superconductors, Phys. Rev. Lett. 122, 126801 (2019).
  * (50) F. Setiawan, A. Stern, and E. Berg, Topological superconductivity in planar Josephson junctions: Narrowing down to the nanowire limit, Phys. Rev. B 99, 220506(R) (2019).
  * (51) S. Ikegaya, S. Tamura, D. Manske, and Y. Tanaka, Anomalous proximity effect of planar topological Josephson junctions, Phys. Rev. B 102, 140505(R) (2020).
  * (52) T. Zhou, M. C. Dartiailh, W. Mayer, J. E. Han, A. Matos-Abiague, J. Shabani, and I. Žutić, Phase Control of Majorana Bound States in a Topological X Junction, Phys. Rev. Lett. 124, 137001 (2020).
  * (53) S.-K. Jian and S. Yin, Chiral topological superconductivity in Josephson junctions, Phys. Rev. B 103, 134514 (2021).
  * (54) A. Fornieri, A. M. Whiticar, F. Setiawan, E. Portolés, A. C. C. Drachmann, A. Keselman, S. Gronin, C. Thomas, T. Wang, R. Kallaher, G. C. Gardner, E. Berg, M. J. Manfra, A. Stern, C. M. Marcus, and F. Nichele, Evidence of topological superconductivity in planar Josephson junctions, Nature 569, 89-92 (2019).
  * (55) H. Ren, F. Pientka, S. Hart, A. T. Pierce, M. Kosowsky, L. Lunczer, R. Schlereth, B. Scharf, E. M. Hankiewicz, L. W. Molenkamp, B. I. Halperin, and A. Yacoby, Topological superconductivity in a phase-controlled Josephson junction, Nature 569, 93-98 (2019).
  * (56) M. C. Dartiailh, W. Mayer, J. Yuan, K. S. Wickramasinghe, A. Matos-Abiague, I. Žutić, and J. Shabani, Phase Signature of Topological Transition in Josephson Junctions, Phys. Rev. Lett. 126, 036802 (2021).
  * (57) J. D. Koralek, C. P. Weber, J. Orenstein, B. A. Bernevig, S.-C. Zhang, S. Mack, and D. D. Awschalom, Emergence of the persistent spin helix in semiconductor quantum wells, Nature 458, 610-613 (2009).
  * (58) M. P. Walser, C. Reichl, W. Wegscheider, and G. Salis, Direct mapping of the formation of a persistent spin helix, Nat. Phys. 8, 752-762 (2012).
  * (59) M. Kohda, V. Lechner, Y. Kunihashi, T. Dollinger, P. Olbrich, C. Schonhuber, I. Caspers, V. V. Bel’kov, L. E. Golub, D. Weiss, K. Richter, J. Nitta, and S. D. Ganichev, Gate-controlled persistent spin helix state in (In,Ga)As quantum wells, Phys. Rev. B 86, 081306(R) (2012).
  * (60) Y. S. Chen, S. Fält, W. Wegscheider, and G. Salis, Unidirectional spin-orbit interaction and spin-helix state in a (110)-oriented GaAs/(Al,Ga)As quantum well, Phys. Rev. B 90, 121304(R) (2014).
  * (61) L. L. Tao and E. Y. Tsymbal, Persistent spin texture enforced by symmetry, Nat. Commun. 9, 2763 (2018).
  * (62) M. A. Absor and F. Ishii, Doping-induced persistent spin helix with a large spin splitting in monolayer SnSe, Phys. Rev. B 99, 075136 (2019).
  * (63) M. A. Absor and F. Ishii, Intrinsic persistent spin helix state in two-dimensional group-IV monochalcogenide $MX$ monolayers ($M=$Sn or Ge and $X=$S, Se, or Te), Phys. Rev. B 100, 115104 (2019).
  * (64) C. Autieri, P. Barone, J. Sławińska, and S. Picozzi, Persistent spin helix in Rashba-Dresselhaus ferroelectric CsBiNb2O7, Phys. Rev. Materials 3, 084416 (2019).
  * (65) H. Ai, X. Ma, X. Shao, W. Li, and M. Zhao, Reversible out-of-plane spin texture in a two-dimensional ferroelectric material for persistent spin helix, Phys. Rev. Mater. 3, 054407 (2019).
  * (66) H. Lee, J. Im, and H. Jin, Emergence of the giant out-of-plane Rashba effect and tunable nanoscale persistent spin helix in ferroelectric SnTe thin films, Appl. Phys. Lett. 116, 022411 (2020).
  * (67) J. Sławińska, F. T. Cerasoli, P. Gopal, M. Costa, S. Curtarolo, and M. B. Nardelli, Ultrathin SnTe films as a route towards all-in-one spintronics devices, 2D Mater., 7, 025026 (2020).
  * (68) S. Barraza-Lopez, B. M. Fregoso, J. W. Villanova, S. S. P. Parkin, and K. Chang, Colloquium: Physical properties of group-IV monochalcogenide monolayers, Rev. Mod. Phys. 93, 011001 (2021).
  * (69) See Supplemental Material at XXX for demonstrating the detailed derivation of the effective Hamiltonian in Eq. (4). We also show $\rho_{\rm edge}(k_{y};E)$ in the presence of the perturbations.
  * (70) P. A. Lee and D. S. Fisher, Anderson Localization in Two Dimensions, Phys. Rev. Lett. 47, 882 (1981).
  * (71) T. Ando, Quantum point contacts in magnetic fields, Phys. Rev. B 44, 8017 (1991).
  * (72) P. Townsend and J. Sutton, Investigation by Electron Tunneling of the Superconducting Energy Gaps in Nb, Ta, Sn, and Pb, Phys. Rev. 128, 591 (1962).
  * (73) C. Honerkamp, K. Wakabayashi and M. Sigrist, Instabilities at [110] surfaces of $d_{x^{2}-y^{2}}$ superconductors Europhys. Lett. 50, 368 (2000).
  * (74) A. C. Potter and P. A. Lee, Edge Ferromagnetism from Majorana Flat Bands: Application to Split Tunneling-Conductance Peaks in High-$T_{c}$ Cuprate Superconductors, Phys. Rev. Lett. 112, 117002 (2014).
  * (75) J. S. Hofmann, F. F. Assaad, and A. P. Schnyder, Edge instabilities of topological superconductors, Phys. Rev. B 93, 201116(R) (2016).

Supplements Materials
“Flat-band Majorana states in topological Josephson junctions”
Daisuke Oshima1, Satoshi Ikegaya1, Andreas P. Schnyder2, and Yukio Tanaka1

1Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan
2Max-Planck-Institut für Festkörperforschung, Heisenbergstrasse 1, D-70569
Stuttgart, Germany

##  V Detail of the derivation of the effective Hamiltonian

In this section, we derive an effective Hamiltonian describing the subgap
states localized around the Josephson junction interface. Let us consider a
Bogoliubov–de Gennes (BdG) Hamiltonian in Eq. (1) of the main text:

$\displaystyle\check{H}(\bm{r})=\check{H}_{0}(\bm{r})+\check{H}_{\lambda}(\bm{r}),$
(32)
$\displaystyle\check{H}_{0}(\bm{r})=\left[\begin{array}[]{cccc}\xi(\bm{r})&V&0&\Delta(z)\\\
V&\xi(\bm{r})&-\Delta(z)&0\\\ 0&-\Delta^{\ast}(z)&-\xi(\bm{r})&-V\\\
\Delta^{\ast}(z)&0&-V&-\xi(\bm{r})\\\ \end{array}\right],$ (37)
$\displaystyle\check{H}_{\lambda}(\bm{r})=\left[\begin{array}[]{cccc}-i\lambda\delta(z)\partial_{x}&0&0&0\\\
0&+i\lambda\delta(z)\partial_{x}&0&0\\\
0&0&-i\lambda\delta(z)\partial_{x}&0\\\
0&0&0&+i\lambda\delta(z)\partial_{x}\\\ \end{array}\right],$ (42)
$\displaystyle\xi(\bm{r})=-\frac{\hbar^{2}}{2m}(\partial^{2}_{x}+\partial^{2}_{y}+\partial^{2}_{z})-\mu,$
(43) $\displaystyle\Delta(z)=\left\\{\begin{array}[]{cl}\Delta
e^{i\varphi/2}&\text{for }z\leq 0\\\ \Delta e^{-i\varphi/2}&\text{for
}z>0\end{array}\right.,$ (46)

where $m$, $\mu$, $V$, and $\lambda$ represent the effective math of an
electron, chemical potential, Zeeman potential, and the strength of the spin-
orbit coupling potential sustaining a persistent spin helix state,
respectively. To derive the effective Hamiltonian for the subgap states, we
treat $\check{H}_{\lambda}$ as a perturbation. Within the zeroth-order
perturbation theory, the energy eigenvalue and wave function of the subgap
state satisfy the BdG equation,

$\displaystyle\check{H}_{0}(\bm{r})\check{U}_{\mathrm{SG}}(\bm{r})=\check{U}_{\mathrm{SG}}(\bm{r})\check{E}_{0},$
(47)
$\displaystyle\check{E}_{0}=\mathrm{diag}[\varepsilon-V,\varepsilon+V,-\varepsilon+V,-\varepsilon-V],$
(48)

under the conditions,

$\displaystyle 0\leq\varepsilon<\Delta$ (49)

and

$\displaystyle\lim_{z\rightarrow\pm\infty}\check{U}_{\mathrm{SG}}(\bm{r})=0.$
(50)

As a preliminary step of obtaining $\varepsilon$ and
$\check{U}_{\mathrm{SG}}(\bm{r})$, we first solve the BdG equation,

$\displaystyle\left[\begin{array}[]{cc}\xi(\bm{r})&s\Delta(z)\\\
s\Delta^{\ast}(z)&-\xi(\bm{r})\end{array}\right]\left[\begin{array}[]{c}u_{s}(\bm{r})\\\
v_{s}(\bm{r})\end{array}\right]=\varepsilon\left[\begin{array}[]{c}u_{s}(\bm{r})\\\
v_{s}(\bm{r})\end{array}\right],$ (57)

for $s=\pm$. The localized wave function for $z\leq 0$ and that for $z>0$ are
represented as

$\displaystyle\psi_{b,s}(\bm{r})=\left[\begin{array}[]{c}u_{b,s}(\bm{r})\\\
v_{b,s}(\bm{r})\end{array}\right]=\left(\frac{A_{b,s}}{\sqrt{2}}\left[\begin{array}[]{c}e^{i\eta}\\\
se^{-i\varphi/2}\end{array}\right]e^{-ik_{+}z}+\frac{B_{b,s}}{\sqrt{2}}\left[\begin{array}[]{c}se^{i\varphi/2}\\\
e^{i\eta}\end{array}\right]e^{ik_{-}z}\right)\frac{e^{ik_{x}x}}{\sqrt{L_{x}}}\frac{e^{ik_{y}y}}{\sqrt{L_{y}}},$
(64)

and

$\displaystyle\psi_{t,s}(\bm{r})=\left[\begin{array}[]{c}u_{t,s}(\bm{r})\\\
v_{t,s}(\bm{r})\end{array}\right]=\left(\frac{A_{t,s}}{\sqrt{2}}\left[\begin{array}[]{c}e^{i\eta}\\\
se^{i\varphi/2}\end{array}\right]e^{ik_{+}z}+\frac{B_{t,s}}{\sqrt{2}}\left[\begin{array}[]{c}se^{-i\varphi/2}\\\
e^{i\eta}\end{array}\right]e^{-ik_{-}z}\right)\frac{e^{ik_{x}x}}{\sqrt{L_{x}}}\frac{e^{ik_{y}y}}{\sqrt{L_{y}}},$
(71)

respectively, where

$\displaystyle
k_{\pm}=\sqrt{\frac{2m}{\hbar^{2}}\sqrt{\mu^{2}_{z}+\Omega^{2}}}e^{\pm
iQ},\quad
2Q=\mathrm{arctan}\frac{\Omega}{\mu_{z}},\quad\mu_{z}=\mu-\frac{\hbar^{2}}{2m}(k_{x}^{2}+k_{y}^{2}),$
(72)
$\displaystyle\eta=\mathrm{arctan}\frac{\Omega}{\varepsilon},\quad\Omega=\sqrt{\Delta^{2}-\varepsilon^{2}},$
(73)

and we apply a periodic boundary condition in the $x-$ ($y$-) direction with
$L_{x(y)}$. Here, $L_{x(y)}$ is the length in the $x$\- ($y$-) direction. From
the boundary condition at the junction interface,

$\displaystyle\lim_{z\rightarrow 0}\psi_{b,s}(\bm{r})=\lim_{z\rightarrow
0}\psi_{t,s}(\bm{r}),$ (74)
$\displaystyle\partial_{z}\psi_{b,s}(\bm{r})|_{z=0}=\partial_{z}\psi_{t,s}(\bm{r})|_{z=0},$
(75)

and a normalization condition,

$\displaystyle\int dx\int
dy\left\\{\int^{0}_{-\infty}dz\;\psi^{\dagger}_{b,s}(\bm{r})\psi_{b,s}(\bm{r})+\int^{\infty}_{0}dz\;\psi^{\dagger}_{t,s}(\bm{r})\psi_{t,s}(\bm{r})\right\\}=1,$
(76)

we obtain the energy eigenvalue as

$\displaystyle\varepsilon=\Delta\sqrt{1-D\sin^{2}(\varphi/2)},$ (77)
$\displaystyle
D=\frac{\Delta^{2}\sin^{2}(\varphi/2)-\mu_{z}^{2}+\mu_{z}\sqrt{2\Delta^{2}\sin^{2}(\varphi/2)+\mu_{z}^{2}}}{2\Delta^{2}\sin^{2}(\varphi/2)},$
(78)

and the numerical coefficients in the wave function as

$\displaystyle A_{b,s}=\frac{-sie^{-iQ/2}f_{-}}{\sqrt{N}},\quad
B_{b,s}=\frac{e^{-i\varphi/2}e^{iQ/2}f_{+}}{\sqrt{N}},\quad
A_{t,s}=\frac{se^{-i\varphi/2}e^{-iQ/2}f_{+}}{\sqrt{N}},\quad
B_{t,s}=\frac{ie^{iQ/2}f_{-}}{\sqrt{N}},$ (79) $\displaystyle
f_{\pm}=\sqrt{\frac{1}{2}\left(1\pm\frac{\sqrt{D\sin^{2}(\varphi/2)}}{\sqrt{1-D\sin^{2}(\varphi/2)}}\cot(\varphi/2)\right)},\quad
N=\frac{1}{\kappa}D(3-2D),\quad\kappa=\sqrt{\frac{2m}{\hbar^{2}}\sqrt{\mu^{2}_{z}+\Omega^{2}}}\sqrt{1-D},$
(80)

where the condition of

$\displaystyle-\frac{1}{2}\Delta|\sin(\varphi/2)|<\mu_{z},$ (81)

must be satisfied to realize the bound states, (i.e.,
$\lim_{z\rightarrow\pm\infty}\psi_{t,s}(\bm{r})=0$). From particle-hole
symmetry of the superconducting system, we also obtain

$\displaystyle\left[\begin{array}[]{cc}\xi(\bm{r})&s\Delta(z)\\\
s\Delta^{\ast}(z)&-\xi(\bm{r})\end{array}\right]\left[\begin{array}[]{c}-v^{\ast}_{s}(\bm{r})\\\
u_{s}^{\ast}(\bm{r})\end{array}\right]=-\varepsilon\left[\begin{array}[]{c}-v^{\ast}_{s}(\bm{r})\\\
u_{s}^{\ast}(\bm{r})\end{array}\right].$ (88)

Finally, we obtain

$\displaystyle\check{U}_{\mathrm{SG}}(\bm{r})=\check{Q}\;\check{U}^{\prime}_{\mathrm{SG}}(\bm{r})\;\check{W},$
(89)
$\displaystyle\check{U}^{\prime}_{\mathrm{SG}}(\bm{r})=\left[\begin{array}[]{cccc}u_{+}(\bm{r})&0&0&-v_{+}^{\ast}(\bm{r})\\\
0&u_{-}(\bm{r})&-v_{-}^{\ast}(\bm{r})&0\\\
0&v_{-}(\bm{r})&u_{-}^{\ast}(\bm{r})&0\\\
v_{+}(\bm{r})&0&0&u_{+}^{\ast}(\bm{r})\\\ \end{array}\right],$ (94)
$\displaystyle\check{Q}=\left[\begin{array}[]{cc}\hat{Q}&0\\\
0&\hat{Q}\end{array}\right],\quad\hat{Q}=\frac{1}{\sqrt{2}}\left[\begin{array}[]{cc}1&1\\\
-1&1\end{array}\right],\quad\check{W}=\left[\begin{array}[]{cc}\hat{w}&0\\\
0&\hat{w}^{\ast}\end{array}\right],\quad\hat{w}=\left[\begin{array}[]{cc}ie^{-i\eta/2}e^{i\varphi/4}&0\\\
0&-ie^{-i\eta/2}e^{i\varphi/4}\end{array}\right],$ (103)

where we introduce $\check{W}$ for later convenience. From the first-order
perturbation theory, we construct the effective Hamiltonian for the subgap
states as,

$\displaystyle\check{H}_{\mathrm{eff}}$ $\displaystyle=\int
d\bm{r}\check{U}^{\dagger}_{\mathrm{SG}}(\bm{r})\check{H}(\bm{r})\check{U}_{\mathrm{SG}}(\bm{r})$
$\displaystyle=\int
d\bm{r}\check{U}^{\dagger}_{\mathrm{SG}}(\bm{r})\check{H}_{0}(\bm{r})\check{U}_{\mathrm{SG}}(\bm{r})+\int
d\bm{r}\check{U}^{\dagger}_{\mathrm{SG}}(\bm{r})\check{H}_{\lambda}(\bm{r})\check{U}_{\mathrm{SG}}(\bm{r})$
$\displaystyle=\left[\begin{array}[]{cccc}\varepsilon-V&\tilde{\lambda}k_{x}&\tilde{\Delta}k_{x}&0\\\
\tilde{\lambda}k_{x}&\varepsilon+V&0&-\tilde{\Delta}k_{x}\\\
\tilde{\Delta}k_{x}&0&-\varepsilon+V&\tilde{\lambda}k_{x}\\\
0&-\tilde{\Delta}k_{x}&\tilde{\lambda}k_{x}&-\varepsilon-V\\\
\end{array}\right],$ (108)

with

$\displaystyle\tilde{\lambda}=\lambda\sqrt{D(1-D)}F,$ (109)
$\displaystyle\tilde{\Delta}=\lambda\cot(\varphi/2)F,$ (110) $\displaystyle
F=\sqrt{\frac{2m}{\hbar^{2}}\sqrt{\mu^{2}_{z}+\Omega^{2}}}\frac{\sqrt{1-D}}{D(3-2D)}\frac{\sqrt{D\sin^{2}(\varphi/2)}}{\sqrt{1-D\sin^{2}(\varphi/2)}},$
(111)

where $\tilde{\lambda}$ and $\tilde{\Delta}$ are the even-functions with
respect to $k_{x}$ and $k_{y}$. As also discussed in the main text, we find
that the effective Hamiltonian of $\check{H}_{\mathrm{eff}}$ is equivalent to
the BdG Hamiltonian of a two-dimensional $p_{x}$-wave superconductor in the
presence of the spin-orbit coupling potential and the Zeeman potential.

Then, we analyze the topological property of $\check{H}_{\mathrm{eff}}$. The
energy eigenvalues of $\check{H}_{\mathrm{eff}}$ are given by

$\displaystyle
E_{s}=\pm\sqrt{\left(\varepsilon+s\sqrt{V^{2}+\tilde{\lambda}^{2}k_{x}^{2}}\right)^{2}+\tilde{\Delta}^{2}k_{x}^{2}},$
(112)

for $s=\pm$. While $E_{+}$ is fully gapped irrespective of
$\bm{k}_{\parallel}=(k_{x},k_{y})$, the branch of $E_{-}$ can have
superconducting gap nodes when

$\displaystyle\varepsilon^{2}=V^{2}+\tilde{\lambda}^{2}k_{x}^{2},$ (113)
$\displaystyle\tilde{\Delta}k_{x}=0,$ (114)

is satisfied. At $k_{x}=0$, this condition is simplified to

$\displaystyle\varepsilon^{2}=V^{2}.$ (115)

From Eq. (115), we find the two superconducting gap nodes when the condition

$\displaystyle X<V^{2}<\Delta^{2},$ (116) $\displaystyle
X=\frac{\mu^{2}+\Delta^{2}\left\\{1+\cos^{2}(\varphi/2)\right\\}-\mu\sqrt{\mu^{2}+2\Delta^{2}\sin^{2}(\varphi/2)}}{2},$
(117)

is satisfied; the two gap nodes are located at

$\displaystyle\bm{k}_{\parallel}=(0,\pm\tilde{k}_{F}),$ (118)
$\displaystyle\tilde{k}_{F}=k_{F}\sqrt{1+\zeta},$ (119)
$\displaystyle\zeta=\frac{1}{\mu}\left[V^{2}-\frac{\Delta^{2}}{2}\left\\{1+\cos^{2}(\varphi/2)\right\\}\right]\sqrt{\frac{1}{V^{2}-\Delta^{2}\cos^{2}(\varphi/2)}},$
(120)

with $k_{F}=\sqrt{2m\mu}/\hbar$. In the limit of $\mu\gg\Delta$, this
condition in Eq. (116) is reduced to

$\displaystyle\Delta^{2}\cos^{2}(\varphi/2)<V^{2}<\Delta^{2}.$ (121)

In addition, $E_{-}$ at $\varphi=\pi$ becomes gapless irrespective of
$\bm{k}_{\parallel}$ because the effective pair potential
$\tilde{\Delta}\propto\cot(\phi/2)$ vanishes. We have two distinct nodal
phases separated by $\varphi=\pi$: the nodal phase I for $\varphi<\pi$ and the
nodal phase II for $\varphi>\pi$, as shown in Fig. 2(a) of the main text. The
effective Hamiltonian of $\check{H}_{\mathrm{eff}}$ has chiral symmetry
defined by

$\displaystyle\check{S}_{\mathrm{eff}}\check{H}_{\mathrm{eff}}(\bm{k}_{\parallel})\check{S}_{\mathrm{eff}}^{-1}=-\check{H}_{\mathrm{eff}}(\bm{k}_{\parallel}),\quad\check{S}_{\mathrm{eff}}=\left[\begin{array}[]{cccc}0&0&-i&0\\\
0&0&0&i\\\ i&0&0&0\\\ 0&-i&0&0\\\ \end{array}\right],$ (126)

where the original BdG Hamiltonian $\check{H}(\bm{r})$ also preserves chiral
symmetry as

$\displaystyle\check{S}H(\bm{r})\check{S}^{-1}=-H(\bm{r}),$ (127)
$\displaystyle\check{S}=\check{M}_{xy}\check{T}\check{C},$ (128)
$\displaystyle\check{M}_{xy}=\left[\begin{array}[]{cccc}-i&0&0&0\\\ 0&i&0&0\\\
0&0&i&0\\\ 0&0&0&-i\\\
\end{array}\right]\mathcal{R}_{z},\quad\check{T}=\left[\begin{array}[]{cccc}0&1&0&0\\\
-1&0&0&0\\\ 0&0&0&1\\\ 0&0&-1&0\\\
\end{array}\right]\mathcal{K},\quad\check{C}=\left[\begin{array}[]{cccc}0&0&1&0\\\
0&0&0&1\\\ 1&0&0&0\\\ 0&1&0&0\\\ \end{array}\right]\mathcal{K},$ (141)

where $\check{M}_{xy}$, $\check{T}$, and $\check{C}$ represent mirror
reflection symmetry with respect to the $xy$-plane, time-reversal symmetry,
and particle-hole symmetry, respectively; $\mathcal{R}_{z}$ describes the
reflection of spatial coordinate z (i.e., $z\rightarrow-z$), and $\mathcal{K}$
is the complex-conjugation operator. We can characterize the topological
property of $\check{H}_{\rm eff}$ by using a one-dimensional winding number
sato_11 ,

$\displaystyle w(k_{y})=\frac{i}{4\pi}\int
dk_{x}\mathrm{Tr}[\check{S}_{\mathrm{eff}}\left\\{\check{H}_{\mathrm{eff}}(\bm{k}_{\parallel})\right\\}^{-1}\partial_{k_{x}}\check{H}_{\mathrm{eff}}(\bm{k}_{\parallel})].$
(142)

It is possible to compute the winding number more simply by employing the
following procedures. First, according to Eq. (115), the emergence/position of
the nodes is irrelevant to $\tilde{\lambda}$. Therefore, since the topological
number is unchanged except when the superconducting gap closes, we can
evaluate the winding number precisely with setting $\tilde{\lambda}\rightarrow
0$. In this limit, $\check{H}_{\mathrm{eff}}$ can be block diagonalized as
$\check{H}_{\mathrm{eff}}\rightarrow\mathrm{diag}[\hat{H}_{+},\hat{H}_{-}]$,
where

$\displaystyle\hat{H}_{s}=\left[\begin{array}[]{cc}\varepsilon+sV&-s\tilde{\Delta}k_{x}\\\
-s\tilde{\Delta}k_{x}&-\varepsilon-sV\end{array}\right]$ (145)

for $s=\pm$. The chiral symmetry of each block component is given by

$\displaystyle\hat{S}_{s}\hat{H}_{s}(\bm{k}_{\parallel})\hat{S}_{s}^{-1}=-\hat{H}_{s}(\bm{k}_{\parallel}),\quad\hat{S}_{s}=s\left[\begin{array}[]{cc}0&i\\\
-i&0\end{array}\right].$ (148)

Since the superconducting gap nodes are originated from the $\hat{H}_{-}$
component, we can assess the winding number solely from $\hat{H}_{-}$. As a
result, the winding number can be calculated by

$\displaystyle w(k_{y})$ $\displaystyle=\frac{i}{4\pi}\int
dk_{x}\mathrm{Tr}[\hat{S}_{-}\left\\{\hat{H}_{-}(\bm{k}_{\parallel})\right\\}^{-1}\partial_{k_{x}}\hat{H}_{-}(\bm{k}_{\parallel})].$
(149)

According to the strict definition for the winding number in continuum model,
the integration in Eq. (149) is performed from $k_{x}=-\infty$ to $\infty$,
where we need to apply an additional hypothetical regulation, i.e.,
$\tilde{\Delta}\rightarrow 0$ at $k_{x}=\pm\infty$, so that we obtain
$\hat{H}_{-}(k_{x}=\infty)=\hat{H}_{-}(k_{x}=-\infty)$ sato_11 . Even so, we
can practically compute the winding number by using a simplified formula,
which is obtained independently of the details of the treatment at
$k_{x}=\pm\infty$ sato_11 :

$\displaystyle
w(k_{y})=\frac{1}{2}\sum_{\varepsilon(\bm{k}_{\parallel})=V}\mathrm{sgn}[\tilde{\Delta}k_{x}]\mathrm{sgn}[\partial_{k_{x}}\varepsilon(\bm{k}_{\parallel})],$
(150)

where $\sum_{\varepsilon(\bm{k}_{\parallel})=V}$ represents the summation over
$k_{x}$ satisfying $\varepsilon(\bm{k}_{\parallel})=V$ with the fixed $k_{y}$.
Finally, we obtain

$\displaystyle
w(k_{y})=\left\\{\begin{array}[]{cl}+1&\text{for}\quad|k_{y}|<\tilde{k}_{F}\\\
0&\text{otherwise}\end{array}\right.,$ (153)

for the nodal phase I in Fig. 2(a) of the main text, and

$\displaystyle
w(k_{y})=\left\\{\begin{array}[]{cl}-1&\text{for}\quad|k_{y}|<\tilde{k}_{F}\\\
0&\text{otherwise}\end{array}\right.,$ (156)

for the nodal phase II in Fig. 2(a) of the main text, while $w(k_{y})=0$
irrespective of $k_{y}$ for the fully gapped phase.

##  VI Flat-band Majorana bound states in the presence of perturbations

In this section, we show $\rho_{\mathrm{edge}}(k_{y};E)$ in the presence of
the perturbations discussed in the main text. Hereafter, the results of
$\rho_{\mathrm{edge}}(k_{y};E)$ are normalized by the same unit as the one in
Fig. 2(c) of the main text, and we plot
$\log_{10}\rho_{\mathrm{edge}}(k_{y};E)$ instead of the raw data of
$\rho_{\mathrm{edge}}(k_{y};E)$. In addition, we show
$\bar{\rho}_{\mathrm{edge}}(E=0)$ normalized by the same unit as the one in
Fig. 2(b) of the main text.

Firstly, we consider the system preserving both chiral symmetry [see Eq. (12)
of the main text] and modified particle-hole symmetry [see Eq. (18) of the
main text], where the junction transparency is decreased by the perturbation
$H_{\delta t}$ with $t_{1}=t_{2}$ [see Eq (28) of the main text]. Figure 4(a)
shows the $\bar{\rho}_{\mathrm{edge}}(E=0)$ with $t_{1}=t_{2}=-2\Delta$ as a
function of $V$ and $\varphi$, which is equivalent to Fig. 3(a) in the main
text. In Figs. 4(b) and 4(c), we show $\rho_{\rm edge}(k_{y};E)$ as a function
of $k_{y}$ and $E$, where we choose (b) $(V,\varphi)=(0.25\Delta,0.9\pi)$ and
(c) $(V,\varphi)=(0.7\Delta,0.65\pi)$, respectively [see also the dots in
Figs. 4(a) ]. Since both chiral symmetry and modified particle-hole symmetry
are preserved, we can find the flat-band MBSs. Nevertheless, in the vicinity
of the topological phase boundary, as shown in Fig. 4(b), the flat-band MBSs
appear at several separated regions with respect to $k_{y}$, where we also
find more than two gap nodes; similar flat-band MBSs also appear in the other
effective $p_{x}$-wave superconductors ikegaya_21 ; ikegaya_18 ; law_13 ;
rosenow_14 . In Fig. 4(d), we plot $\log_{10}\rho_{\rm edge}(k_{y};E)$ for the
ideal junction (i.e., $t_{1}=t_{2}=0$), where we choose $(V,\varphi)$ as same
as that in Fig. 4(c) [i.e., $(V,\varphi)=(0.7\Delta,0.65\pi)$]. From the
comparison between Fig. 4(c) and Fig. 4(d), especially around $k_{y}=0$, we
see that the induced gap is suppressed by decreasing the junction
transparency. In Figs. 5(a) and 5(b), we consider the junction with decreasing
the transparency further, i.e., $t_{1}=t_{2}=-5\Delta$. In Fig. 5(a), we plot
$\bar{\rho}_{\mathrm{edge}}(E=0)$ in the range from $0$ to $0.4$, while that
in Fig. 4(a) is shown in the range from $0$ to $1$ [see the color bars in
Figs. 4(a) and 5(a)]. We see that $\bar{\rho}_{\mathrm{edge}}(E=0)$ is
suppressed by decreasing the junction transparency. Figure 5(b) shows
$\rho_{\rm edge}(k_{y};E)$ at $(V,\varphi)=(0.8\Delta,0.7\pi)$, which is
located at deep inside of the topologically nontrivial region. We see that the
suppression of the induced gap becomes more evident by decreasing the junction
transparency. Especially, the induced gap around $k_{y}=0$ is disturbed
devastatingly, and thus the flat-band MBSs around $k_{y}=0$ disappear. We note
that the reduction of the induced gap by decreasing the junction transparency
is also found in the planar TJJs halperin_17 .

Secondary, we consider the system with broken chiral symmetry by $t_{1}\neq
t_{2}$, where the modified particle-hole symmetry is still preserved. Figure
6(a) shows the $\bar{\rho}_{\mathrm{edge}}(E=0)$ with $t_{1}=2\Delta$ and
$t_{2}=-2\Delta$, which is equivalent to Fig. 3(b) in the main text. In Figs.
6(b) and 6(c), we show $\rho_{\rm edge}(k_{y};E)$, where we choose (b)
$(V,\varphi)=(0.7\Delta,0.85\pi)$ and (c) $(V,\varphi)=(0.7\Delta,\pi)$,
respectively. Although the chiral symmetry is broken, since the modified
particle-hole symmetry is still preserved, we clearly find the emergence of
the flat-band MBSs protected by the $\mathbb{Z}_{2}$ index. Moreover, as shown
in Figs. 6(c), we obtain the flat-band MBSs at $\varphi=\pi$. The results
imply that the two nodal phases originally separated by $\varphi=\pi$ merge
into a single nodal phase hosting the flat-band MBSs.

Thirdly, we consider the system in the presence of the Rashbe spin-orbit
coupling (SOC) potential [see Eq. (29) of the main text] with preserving the
chiral symmetry (i.e., $\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}$), while
the modified particle-hole symmetry is broken. Figure 7(a) shows the
$\bar{\rho}_{\mathrm{edge}}(E=0)$ with
$\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}=2.5\Delta$, which is equivalent to
Fig. 3(c) in the main text. In Figs. 7(b) and 7(c), we show $\rho_{\rm
edge}(k_{y};E)$ at (b) $(V,\varphi)=(0.7\Delta,0.7\pi)$ and (c)
$(V,\varphi)=(0.7\Delta,1.3\pi)$, respectively. Since the chiral symmetry is
preserved, we obtain the flat-band MBSs. We additionally find that the induced
gaps become asymmetric with respect to $\varphi=\pi$, i.e., the induced gap
for $\varphi<\pi$ [Fig. 7(b)] is slightly smaller than that for $\varphi>\pi$
[Fig. 7(c)]. Such tendency becomes more evident by increasing the Rashba SOC
potential. Figure 8(a) shows the $\bar{\rho}_{\mathrm{edge}}(E=0)$ with
$\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}=5\Delta$. In Figs. 8(b) and 8(c),
we show $\rho_{\rm edge}(k_{y};E)$ at (b) $(V,\varphi)=(0.65\Delta,0.7\pi)$
and (c) $(V,\varphi)=(0.75\Delta,1.25\pi)$, respectively. As shown in Fig.
8(a), $\bar{\rho}_{\mathrm{edge}}(E=0)$ for $\varphi<\pi$ is distinctly
smaller than that for $\varphi>\pi$. As shown in Figs. 8(b) and 8(c), although
we obtain the flat-band MBSs in both $\varphi<\pi$ and $\varphi>\pi$, the
induced gap for $\varphi<\pi$ is strongly suppressed; we additionally find
that the energy spectrum has more than two gap nodes. The decay length of the
flat-band MBS becomes longer as the energy gap becomes smaller, and thus the
amplitude of wave function of the flat-band MBSs at the surface (i.e.,
$i_{x}=0$) becomes smaller. This explains the suppression of surface LDOS,
$\bar{\rho}_{\rm edge}(E=0)$, for $\varphi<\pi$. Strictly speaking, the phase
boundary between two distinct nodal phases slightly deviates from
$\varphi=\pi$; see dotted lines in Fig. 7(a) and Fig. 8(a). Actually, as shown
in Fig. 7(d), we confirm that the flat-band MBSs appears event at
$\varphi=\pi$.

Lastly, we consider the Rashba SOC potential breaking both chiral symmetry and
modified particle-hole symmetry (i.e.,
$\lambda^{\prime}_{1}\neq-\lambda^{\prime}_{2}$). Figure 9(a) shows the
$\bar{\rho}_{\mathrm{edge}}(E=0)$ with $\lambda^{\prime}_{1}=3.5\Delta$ and
$\lambda^{\prime}_{2}=-1.5\Delta$, which is equivalent to Fig. 3(d) in the
main text. The strong suppression in $\bar{\rho}_{\mathrm{edge}}(E=0)$
suggests the absence of the flat-band MBSs. In Fig. 9(b), we show $\rho_{\rm
edge}(k_{y};E)$ at $(V,\varphi)=(0.5\Delta,0.8\pi)$. We find that the
asymmetric Rashba SOC potential causes the finite slope in the dispersion of
flat-band MBSs.

Figure 4: (a) Zero-energy LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as a function of
the Zeeman potential $V$ and the superconducting phase difference $\varphi$,
where we choose $t_{1}=t_{2}=-2\Delta$; the system preserves both chiral
symmetry and modified particle-hole symmetry. This figure is equivalent to
Fig. 3(a) in the main text. In (b) and (c), we show $\rho_{\rm edge}(k_{y};E)$
on a log-scale, i.e., $\log_{10}\rho_{\rm edge}(k_{y};E)$, as a function of
the energy $E$ and the momentum $k_{y}$, where we choose (b)
$(V,\varphi)=(0.25\Delta,0.9\pi)$ and (c) $(V,\varphi)=(0.7\Delta,0.65\pi)$,
respectively. In (d), for comparison, we plot $\log_{10}\rho_{\rm
edge}(k_{y};E)$ for the ideal junction (i.e., $t_{1}=t_{2}=0$), where we
choose $(V,\varphi)=(0.7\Delta,0.65\pi)$. Figure 5: (a) Zero-energy LDOS,
$\bar{\rho}_{\rm edge}(E=0)$, as a function of $V$ and $\varphi$, where the
system preserves both chiral symmetry and modified particle-hole symmetry with
$t_{1}=t_{2}=-5\Delta$. In (b), we show $\log_{10}\rho_{\rm edge}(k_{y};E)$ as
a function of $E$ and $k_{y}$, where we choose
$(V,\varphi)=(0.8\Delta,0.7\pi)$. Figure 6: (a) Zero-energy LDOS,
$\bar{\rho}_{\rm edge}(E=0)$, as a function of $V$ and $\varphi$, where we
consider the asymmetric junction transparency with $t_{1}=2\Delta$ and
$t_{2}=-2\Delta$; the chiral symmetry is broken, while the modified particle-
hole symmetry is still preserved. This figure is equivalent to Fig. 3(b) in
the main text. In (b) and (c), we show $\log_{10}\rho_{\rm edge}(k_{y};E)$ as
a function of $E$ and $k_{y}$, where we choose (b)
$(V,\varphi)=(0.7\Delta,0.85\pi)$ and (c) $(V,\varphi)=(0.7\Delta,\pi)$,
respectively. Figure 7: (a) Zero-energy LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as
a function of $V$ and $\varphi$, where we consider the Rashba SOC potential
with $\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}=2.5\Delta$; the modified
particle-hole symmetry is broken, while the chiral symmetry is still
preserved. The dotted line indicates $(V,\varphi)$ taking the local minima of
$\bar{\rho}_{\rm edge}(E=0)$. This figure is equivalent to Fig. 3(c) in the
main text. In (b)-(d), we show $\log_{10}\rho_{\rm edge}(k_{y};E)$ as a
function of $E$ and $k_{y}$, where we choose (b)
$(V,\varphi)=(0.7\Delta,0.7\pi)$, (c) $(V,\varphi)=(0.7\Delta,1.3\pi)$, and
(d) $(V,\varphi)=(0.9\Delta,\pi)$, respectively. Figure 8: (a) Zero-energy
LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as a function of $V$ and $\varphi$, where
we choose $\lambda^{\prime}_{1}=-\lambda^{\prime}_{2}=5\Delta$; the modified
particle-hole symmetry is broken, while the chiral symmetry is still
preserved. The dotted line indicates $(V,\varphi)$ taking the local minima of
$\bar{\rho}_{\rm edge}(E=0)$. In (b) and (c), we show $\log_{10}\rho_{\rm
edge}(k_{y};E)$ as a function of $E$ and $k_{y}$, where we choose (b)
$(V,\varphi)=(0.65\Delta,0.7\pi)$ and (c) $(V,\varphi)=(0.75\Delta,1.25\pi)$,
respectively. Figure 9: (a) Zero-energy LDOS, $\bar{\rho}_{\rm edge}(E=0)$, as
a function of $V$ and $\varphi$, where both chiral symmetry and modified
particle-hole symmetry are broken by the Rashba SOC potentials with
$\lambda^{\prime}_{1}=3.5\Delta$ and $\lambda^{\prime}_{2}=-1.5\Delta$. This
figure is equivalent to Fig. 3(d) in the main text. In (b), we show
$\log_{10}\rho_{\rm edge}(k_{y};E)$ as a function of $E$ and $k_{y}$, where we
choose $(V,\varphi)=(0.5\Delta,0.8\pi)$.

## References

  * (1) M. Sato, Y. Tanaka, K. Yada, and T. Yokoyama, Phys. Rev. B 83, 224511 (2011).
  * (2) S. Ikegaya, J. Lee, A. P. Schnyder, and Y. Asano, Phys. Rev. B 104, L020502 (2021).
  * (3) S. Ikegaya, S. Kobayashi, and Y. Asano, Phys. Rev. B 97, 174501 (2018).
  * (4) C. L. M. Wong, J. Liu, K. T. Law, P. A. Lee, Phys. Rev. B 88, 060504(R) (2013).
  * (5) T. Hyart, A. R. Wright, and B. Rosenow, Phys. Rev. B 90, 064507 (2014).
  * (6) F. Pientka, A. Keselman, E. Berg, A. Yacoby, A. Stern, and B. I. Halperin, Phys. Rev. X 7, 021032 (2017).


# Muons for cultural heritage

M. Moussawi1,∗ A. Giammanco1 V. Kumar1 M. Lagrange1

(1Centre for Cosmology, Particle Physics and Phenomenology (CP3), Université
catholique de Louvain, Chemin du Cyclotron 2, B-1348 Louvain-la-Neuve, Belgium
)

###### Abstract

Non-destructive subsurface imaging methods based on the absorption or
scattering of photons or neutrons are becoming increasingly popular in
cultural asset conservation. However, these techniques are limited by physical
and practical issues: their penetration depth may be insufficient for large
items, and they usually necessitate transferring the objects of interest to
specialised laboratories. The latter issue is recently being addressed by the
development of portable sources, but artificial radiation can be harmful and
is thus subjected to strict regulation. Muons are elementary particles that
are abundantly and freely created in the atmosphere by cosmic-ray
interactions. Their absorption and scattering in matter are respectively
dependent on the density and elemental composition of the substance they
traverse, suggesting that they could be used for subsurface remote imaging.
This novel technique, dubbed ”muography,” has been used in applications
ranging from geophysics to archaeology, but has remained largely unexplored
for a wide range of cultural heritage objects that are small by muography
standards but whose size and density are too large for conventional imaging
methods. This document outlines the general arguments and some early
simulation studies that aim at exploring the low-size limit of muography and
its relevance for cultural heritage preservation.

Proceedings of the Muon4Future Conference, 29-31 mai. 2023 at Venice, Italy.
Submitted to Proceeding in Science.

CP3-23-50

##  1 Introduction

Imaging methods based on X-rays have been widely used in the context of
cultural heritage preservation [1] due to their ability to penetrate various
materials. However, X-ray imaging has limitations when dealing with large or
dense objects like compact stone or metal, since they do not penetrate deep
enough. Alternative radiation types, such as MeV-range X-rays and neutrons,
offer some improvement but face challenges in transporting valuable objects to
specialized imaging facilities due to size, weight, and preservation concerns.
Various portable setups, such as X-ray fluorescence analysis (XRF) [2, 3] and
portable X-ray computed tomography (CT) systems [4], are available for
cultural heritage studies, but they have limitations in depth penetration and
radiation hazards. Neutron sources [5] offer greater depth [6], but raise
concerns about material activation. A recent advancement using a portable
proton accelerator [7] shows promise but also suffers from radiation hazard
concerns. In contrast, muography [8], which utilizes muons ($\mu$), elementary
particles naturally generated by cosmic-ray interactions in the atmosphere,
represents a promising solution. Cosmogenic muons have remarkable penetrating
capabilities, making them ideal for sub-surface imaging in a variety of
contexts including, as we argue in this paper, cultural heritage applications.
This technique includes two main methods: scattering-based and absorption-
based, sketched in Fig.2 and Fig.2 respectively. Absorption-based muography
measures muons absorption rate within materials, providing insights into their
density and composition, while scattering-based muography utilizes the
diffusion of muons to discriminate elements in multi-material objects [9].

Figure 1: In scattering muography, the object under investigation is
”sandwiched” between muon trackers.

Figure 2: In absorption muography, muon trackers are downstream of the object
of interest; 3D imaging can be obtained by combining multiple viewpoints.

Muography has proven to be highly effective in investigating cultural heritage
sites. The ScanPyramids project, utilizing absorption-based muography, made
the headlines by revealing an unexpected low-density anomaly deeply inside
Khufu’s Great Pyramid in 2017 [10], and then precisely characterizing a
previously unknown corridor in 2023 [11], which has been confirmed by visual
inspection via an endoscope. In another example, density anomalies
(potentially posing safety hazards) have been found in a rampart of a
defensive wall of Xi’an (China) in 2022 [12]. Furthermore, scattering-based
muography has been proposed to search for iron chains within the brickwork of
the Florence cathedral’s dome in Italy [13], and a proof-of-principle test on
a mock-up wall was successfully conducted to demonstrate the conceptual
validity of the method.

While most examples so far are applications to very large volumes of interest,
this paper advocates for the adoption of portable and safe muography as a
promising imaging approach for cultural heritage studies in a regime that is
new for muography (relatively low size) while being beyond reach for methods
based on other radiation sources. A preliminary simulation study using Geant4
[14] illustrates the potential applications and limitations of muography.

Figure 3: Picture of the wooden statue at the Africa Museum of Tervuren. From
project TOCOWO (https://tocowo.ugent.be/).

Figure 4: Geant4 simulation setup; green and blue panels represent the
scattering and absorption setups, respectively.

##  2 Simulated case studies

Each of the two muography techniques has its own sensitivity, applicability,
and limits.

In absorption muography, a single muon tracker is able to measure the 2D
projection of matter density, and the combination of measurements from
different viewpoints can give a 3D density map. However, it provides no
material discrimination apart from density, and small-size or low-density
objects do not stop enough muons to provide sufficient contrast.

In scattering muography, at least two muon trackers are needed, upstream and
downstream of the object of interest, to reconstruct the $\mu$ trajectory
before and after passing through it. This method naturally yields 3D
information, and is sensitive to elemental composition because the width of
the scattering angle distribution is a function of atomic number Z. However,
it is impractical for human-sized statues, as the object of interest must fit
between the two trackers. Either the object of interest is moved inside the
set-up, or a rather complex installation of the detectors must be performed
around the object. Therefore, this method is appropriate for relatively low-
size objects.

We perform a Monte Carlo simulation using CRY [15] to generate muons and a
Geant4 [14] model of an African statue (Fig. 4) made of hardwood. This object,
40 cm tall, has been studied with X-rays but its size nears the limit of that
method, while it is very small by muography standards. To investigate the
potential of muography for material identification, we scale the statue’s size
by factors two and four and introduce hidden cylinders of different materials
within its internal structure, as summarized in Table 1. In this first
exploratory study, we model an ideal detector (i.e. with 100% efficiency and
perfect resolution) made up of nine planes, as illustrated in Fig. 4. Among
these planes, the six indicated in green surround the target and are used for
scattering reconstruction, while the remaining three (in blue) are used in the
absorption reconstruction study.

Scenario | Statue size [cm3] | Cylinder material | Cylinder radius [cm]
---|---|---|---
I (a) | $80\times 30\times 30$ | $/$ | $/$
I (b) | $80\times 30\times 30$ | Air | 5
I (c) | $80\times 30\times 30$ | Bronze bar | 5
II | $160\times 60\times 60$ | Bronze bar | 10
Table 1: Different simulation scenarios.

Scattering reconstruction

Scattering muography is based on the measurement of muon deflections when
passing through an object. The deflection angle is measured by extrapolating
the incoming and outgoing trajectories observed by the two trackers and
determining their point of closest approach (POCA). This approach relies on
the interpretation of a POCA point as the actual place where the muon had a
single high-energy elastic interaction with a nucleus, neglecting the
occurrence of other electromagnetic interactions along its trajectory, which
is rough approximation of reality but has proven to be effective in many
applications. Figure 5 shows the distribution of POCA points obtained in the
three simulated scenarios denoted as I (a, b, c) in Table 1. These plots are
based on 5 million muons, roughly corresponding to an acquisition time of
$\sim 8$ hours, and they show how challenging it is to find a cavity within
this kind of statue, as opposed to finding a high-density insertion.

Figure 5: Distribution of POCA points, projected to a 2D plane for clarity, in
three simulated scenarios: (left) actual wooden sculpture, (middle) with a
cylindrical cavity, (right) with a cylindrical bronze rod.

The output of the muon-scattering reconstruction algorithm is a 3D
distribution of POCA points, each associated to a scattering angle. Based on
those raw data, some clustering algorithms can be used in order to
discriminate between different material densities and elements. At present two
methods are employed to analyze the object’s content: DBSCAN [16] and a
Neighborhood Sum [17] algorithm. DBSCAN only relies on the density of the POCA
points, while the Neighborhood Sum method can consider both the density of
POCA points and the scattering angle of the tracks. We apply DBSCAN in two
steps, to first remove noise points and then separate the volumes
corresponding to different materials using tighter clustering criteria; the
result is shown in Fig. 7. We applied the Neighborhood Sum without (Fig. 7
left) and with (Fig. 7 right) considering the additional information from the
scattering angles, and we obtain in both cases a good discrimination of the
two materials. This discrimination is not as precise as DBSCAN, however this
method is more appropriate for scenarios with low exposure times, where POCA
points are scarce, and the quantitative results it provides are still
reliable.

Figure 6: Using DBScan clustering algorithm (blue: hidden copper bar, yellow:
wood statue).

Figure 7: Using Neighborhood Sum clustering algorithm (orange: hidden copper
bar, magenta: wood statue), taking as input only the POCA positions (left) or
the scattering angles (right).

Absorption reconstruction

With the scenario II described in Table 1, we explore a challenging regime in
which the statue is very big for scattering muography and very small for
absorption muography. It is customary in this method, when the volume of
interest is very distant from the detector (e.g. when imaging the summit of a
volcano), to approximate the latter with a point, meaning that only the zenith
and azimuth angles ($\theta,\phi$) are important while the entry point of the
muon in the detector is not. However, to study human-sized sculptures we have
in general the possibility to position the detectors very close to the statue,
in order to maximize the resolution within the object, and this approximation
is no longer valid. For this study we develop a custom back-projection
reconstruction algorithm inspired by the methods of Refs. [18, 19]. As
illustrated in Fig.9, we extrapolate each muon track onto a voxelized volume,
and we count the number of times a voxel is hit by this backprojected
trajectory. Figure 9, based on the equivalent of two hours of data
acquisition, shows the 3D transmission map slice by slice in the voxelized
volume after selecting only muons with $E<800$ MeV, assuming that the detector
setup also contains a way to discriminate the muons above and below this
energy threshold. Energy discrimination can be achieved cheaply by introducing
a passive absorber before the last detector layer, used as a veto for
energetic muons, or more precisely by combining absorption, scattering, time
of flight, or other variables.

Figure 8: Sketch for the back-projection algorithm.

Figure 9: Transmission map slice by slice.

##  3 Conclusion and prospects

In this paper, we outlined the strengths and limitations of muography for
cultural heritage applications, and detailed a preliminary simulation study
with both scattering and absorption muography for the imaging of statues with
different size and containing different hidden materials. The next step will
be a systematic comparison of several more scenarios, in terms of material and
size of the statue and of the hidden volume. We will also take into account
various realistic scenarios for the detector resolution and the geometry of
the setup (e.g. distances between planes), to identify the best trade-off
between cost and statistical identification power, in preparation for actual
measurements with test objects.

Muography is inexpensive and portable; thanks to the muon penetration power,
it is complementary to other imaging methods. Absorption and scattering have
complementary strengths and weaknesses, but some limitations are in common for
both: long acquisition times are necessary, due to the relatively low natural
rate, and muon direction and energy cannot be controlled.

At the workshop, we were invited to comment on what an artificial muon beam
could do for this kind of studies. One could, indeed, overcome the
aforementioned muography drawbacks by using a muon beam where both muon energy
and direction could be controlled. Even a modest precision in these two
variables and a modest beam luminosity, by accelerator standards, would allow
to do much better than with muons from cosmic rays. The beam energy could be
optimized a priori based on the size and on the main material of the object;
if the inner composition is completely unknown, particularly interesting is
the possibility to scan the same object with beams at various energies. We
would benefit most if a transportable artificial muon source became
accessible. However, this would bring radiological hazards, like all methods
based on an artificial particle source, because of the byproducts of the
collisions needed in order to produce muons and antimuons.

## Acknowledgements

We are indebted to Tim De Kock of the Antwerp Cultural Heritage Sciences
(ARCHES) department at the University of Antwerp, Judy De Roy and Sam Huysmans
of the Royal Institute for Cultural Heritage (KIK-IRPA), and Matthieu Boone of
the Ghent University Centre for X-ray Tomography at the University of Ghent,
for their guidance on the definition of the targets of interest for this
potential application of muography. We thank the Africa Museum of Tervuren and
the project TOCOWO (https://tocowo.ugent.be/) for the model of a wooden statue
(Figure 4). This work was partially supported by the Fonds de la Recherche
Scientifique - FNRS under Grants No. T.0099.19 and J.0070.21, and by the EU
Horizon 2020 Research and Innovation Programme under the Grant Agreements No.
822185 (“INTENSE”) and No. 101021812 (“SilentBorder”).

## References

  * [1] M. P. Morigi _et al._ , “Application of X-ray computed tomography to cultural heritage diagnostics,” _Applied Physics A_ , vol. 100, no. 3, pp. 653–661, 2010.
  * [2] A. Bezur _et al._ , _Handheld XRF in Cultural Heritage: A Practical Workbook for Conservators_. Getty Conservation Institute, 2020.
  * [3] S. Ridolfi, “Portable x-ray fluorescence spectrometry for the analyses of cultural heritage,” in _IOP Conference Series: Materials Science and Engineering_ , vol. 37. IOP Publishing, 2012, p. 012001.
  * [4] F. Albertin _et al._ , “X-Ray Computed Tomography In Situ: An Opportunity for Museums and Restoration Laboratories,” _Heritage_ , vol. 2, no. 3, pp. 2028–2038, 2019.
  * [5] P. Kerr _et al._ , “Neutron transmission imaging with a portable dt neutron generator,” _Radiation Detection Technology and Methods_ , pp. 1–10, 2022.
  * [6] N. Kardjilov _et al._ , “Advances in neutron imaging,” _Materials Today_ , vol. 21, no. 6, pp. 652–672, 2018.
  * [7] F. Taccetti _et al._ , “MACHINA, the Movable Accelerator for Cultural Heritage In-situ Non-destructive Analysis: project overview,” _Rendiconti Lincei. Scienze Fisiche e Naturali_ , pp. 1–19, 2023.
  * [8] L. Bonechi _et al._ , “Atmospheric muons as an imaging tool,” _Rev. Phys._ , vol. 5, p. 100038, 2020.
  * [9] K. N. Borozdin _et al._ , “Radiographic imaging with cosmic-ray muons,” _Nature_ , vol. 422, p. 277, 2003.
  * [10] K. Morishima _et al._ , “Discovery of a big void in Khufu’s Pyramid by observation of cosmic-ray muons,” _Nature_ , vol. 552, no. 7685, p. 386, 2017.
  * [11] S. Procureur _et al._ , “Precise characterization of a corridor-shaped structure in Khufu’s Pyramid by observation of cosmic-ray muons,” _Nature Communications_ , vol. 14, no. 1, p. 1144, 2023.
  * [12] G. Liu _et al._ , “High-precision muography in archaeogeophysics: A case study on Xi’an defensive walls,” _Journal of Applied Physics_ , vol. 133, no. 1, p. 014901, 2023.
  * [13] E. Guardincerri _et al._ , “Imaging the dome of Santa Maria del Fiore using cosmic rays,” _Phil. Trans. R. Soc. A_ , vol. 377, p. 0136, 2018.
  * [14] S. Agostinelli _et al._ , “GEANT4 — a simulation toolkit,” _Nucl. Inst. Meth. A_ , vol. 506, p. 250, 2003.
  * [15] C. Hagmann, D. Lange, and D. Wright, “Cosmic-ray shower generator (CRY) for Monte Carlo transport codes,” in _2007 IEEE Nuclear Science Symposium_ , vol. 2. IEEE, 2007, pp. 1143–1146.
  * [16] M. Ester _et al._ , “A density-based algorithm for discovering clusters in large spatial databases with noise,” in _KDD-96 proceedings_ , ser. 96, no. 34, 1996, pp. 226–231.
  * [17] V. Kumar, _Comparative study of gas detectors and their suitability for imaging_. Mumbai: HBNI, 2022. [Online]. Available: http://www.hbni.ac.in/phdthesis/phys/PHYS05201604011.pdf
  * [18] M. Stapleton _et al._ , “Angle statistics reconstruction: a robust reconstruction algorithm for muon scattering tomography,” _JINST_ , vol. 9, no. 11, p. P11019, 2014.
  * [19] L. Bonechi _et al._ , “A projective reconstruction method of underground or hidden structures using atmospheric muon absorption data,” _JINST_ , vol. 10, no. 02, p. P02003, 2015.