\pdfoutput=1
\documentclass[aps,prb,reprint,superscriptaddress,nofootinbib]{revtex4-1}


\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{dcolumn}
\usepackage{dsfont}
\usepackage{latexsym}
\usepackage{rotating}
\usepackage{color}
\usepackage{latexsym}
\usepackage{bbm}
\usepackage{subfigure}
\usepackage{float}
\usepackage{epsfig}
\usepackage{psfrag}
\usepackage{natbib, hyperref}
\usepackage{bm}
\usepackage{amsthm}
\usepackage{eucal}
\usepackage{mathrsfs}
\usepackage{url}
\usepackage{braket}
\usepackage{ulem}
\usepackage{calligra}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{newunicodechar}
\usepackage{marvosym}
\hyphenation{im-pre-se}
\usepackage[absolute]{textpos}
\usepackage[cal=cm]{mathalfa}
\usepackage{soul}
\def\bibsection{\section*{\refname}} 

\newcommand{\footremember}[2]{%
    \footnote{#2}
    \newcounter{#1}
    \setcounter{#1}{\value{footnote}}%
}

\usepackage{color} 
\newcommand{\red}{\color[rgb]{0.8,0,0}}
\newcommand{\blue}{\color[rgb]{0,0,0.6}}
\newcommand{\black}{\color[rgb]{0,0,0.9}}

\usepackage{hyperref}
\hypersetup{
colorlinks=true,final=true,
        linkcolor=blue,
        citecolor=blue,
        filecolor=blue,
        urlcolor=blue,
}

\begin{document}
\setlength{\abovedisplayskip}{3pt}
\setlength{\belowdisplayskip}{3pt}


%--------------------------------------------------------------------------------------
\setcounter{figure}{0}
\renewcommand{\figurename}{{\bf Fig.}}
\renewcommand{\thefigure}{{\bf \arabic{figure}}}
%--------------------------------------------------------------------------------------


\title{Majorana corner states on the dice lattice}


\author{Narayan Mohanta}
\email{narayan.mohanta@ph.iitr.ac.in}
\affiliation{Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA}
\affiliation{Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA}
\affiliation{Department of Physics, Indian Institute of Technology Roorkee, Roorkee 247667, India}

\author{Rahul Soni}
\affiliation{Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA}
\affiliation{Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA}

\author{Satoshi Okamoto}
\affiliation{Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA}

\author{Elbio Dagotto}
\affiliation{Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA}
\affiliation{Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA}

\begin{abstract}
Lattice geometry continues providing exotic new topological phases in condensed matter physics. Exciting recent examples are the higher-order topological phases, manifesting via localized lower-dimensional boundary states. Moreover, flat electronic bands with a non-trivial topology arise in various lattices and can hold a finite superfluid density, bounded by the Chern number $C$. Here we provide a general route to analyze the topological properties of flat bands in the superconducting state. We argue that the topological superconductivity induced in a flat band with $C\!=\!n$ is of order $n$ and the associated boundary Majorana states have a complex structure when $n\!>\!1$. As example, we show that an attractive interaction in the dice lattice, that exhibits flat bands with $C\!=\!2$, gives rise to second-order topological superconductivity with mixed singlet-triplet pairing. The second-order nature of the topological superconducting phase is revealed by the zero-energy Majorana bound states at the lattice corners. These findings suggest that flat bands with a finite Chern number provide feasible platforms for inducing topological superconductivity of various orders. 
\end{abstract}


\maketitle

%\noindent  \textbf{\small INTRODUCTION}\\
\section*{Introduction}
Higher order topology in quantum matter has recently generated a flurry of activity in several broad areas, including the field of superconductivity~\cite{Schindler_SciAdv2018,PhysRevLett.119.246401,PhysRevLett.120.026801,PhysRevB.97.205136,PhysRevB.98.201114,PhysRevB.97.205135,PhysRevX.9.011012,PhysRevB.99.041301,Li_2DMater2021,Ghosh_PRB2021,PhysRevResearch.2.012060,PhysRevResearch.2.043300}. At the boundaries and vortex cores, 
a topological superconductor harbors Majorana quasiparticles, with potential value in the long-sought area of decoherence-free quantum computing~\cite{Kitaev_2001,RevModPhys.80.1083,Sarma2015,PhysRevX.6.031016,RevModPhys.87.137,Mohanta_EPL2014}. Topological superconductivity can be induced, for example, by a Rashba spin-orbit coupling together with a magnetic field~\cite{PhysRevLett.105.077001,Mourik_Science2012,Rokhinson2012,Deng_Science2016}, and also by a spatially-modulated spin texture in proximity to a conventional superconductor~\cite{Desjardins2019,PhysRevApplied.12.034048,Mohanta_CommPhys2021,Herbrych2021}. A $n^{\rm th}$-order topological superconductor in $d$ dimensions hosts $(d\!-\!n)$-dimensional Majorana states~\cite{Song_PRL2017}. The corner-localized Majorana bound states (MBS) in a second-order two-dimensional topological superconductor are particularly interesting because a two-dimensional array of corner MBS -- useful for demonstrating non-Abelian statistics -- can be easily achieved~\cite{PhysRevB.98.165144,PhysRevResearch.2.032068,PhysRevResearch.2.043025}. These corner MBS have been proposed in many platforms including a topological insulator in proximity to a $d$-wave or $s^{\pm}$-wave superconductor, and a Josephson junction bilayer~\cite{PhysRevLett.121.096803,PhysRevLett.121.186801,PhysRevLett.122.126402}.

Following the discovery of unconventional superconductivity in twisted-bilayer graphene~\cite{Cao2018}, a series of studies suggested the possibility of electronic pairing from repulsive interaction in materials with a high density of states at the Fermi level, such as in a flat electronic band, leading to superconductivity with a high critical temperature~\cite{PhysRevB.101.014501,PhysRevLett.126.027002,Volovik2016,Aoki2020,PhysRevB.106.125155}. When a flat band is topologically nontrivial, the topological invariant places a lower bound on the superfluid weight $D_s$ \textit{i.e.} $D_s \! \ge \!C$, where $C$ is the Chern number of the flat band~\cite{Peotta2015,PhysRevLett.124.167002}. In this case, near a band-inversion wavevector, the Berry phase can convert a repulsive interaction between two oppositely-moving electrons into an effective attraction. 

Therefore, an important question that emerges is whether topological superconductivity can be realized in a topologically non-trivial flat band. Here, we address this general question and argue that there is a connection between the Chern number of the topological flat band and the order of the induced topological superconductivity. As an example, we consider the topological flat bands with $C\!=\!2$ of the dice lattice~\cite{PhysRevB.84.241103,PhysRevB.102.045105} in the presence of an attractive interaction and show that second-order topological superconductivity is induced by populating these flat band at the Fermi level. A hallmark of the induced second-order topological superconducting phase is found via the zero-energy MBS localized at the lattice corners, as shown schematically in Fig.~\ref{dice_flatbands}{\bf {\small e}}.
%---------------------------------------------
\begin{figure*}[t]
\begin{center}
\vspace{0mm}
\epsfig{file=Fig1.pdf,trim=0.0in 0.0in 0.0in 0.0in,clip=false, width=180mm}
\caption{{\bf Majorana corner states and topological flat bands in the dice lattice.} {\bf {\small a}} Dice lattice with open boundary conditions (OBC) and the schematic description of the corner-bound Majorana states in the second-order topological superconducting phase. There are three inequivalent sites in the unit cell, shown by the dashed lines. The triangles denote three-coordination sites and the hexagrams denote six-coordination sites. The black arrows surrounding the six-coordination site (middle layer) represent the vectors of the Rashba spin-orbit coupling, with clockwise sense of rotation for the upper triangles (bottom  layer) and counter-clockwise for the lower triangles (top layer). {\bf {\small b}} Band dispersion of the dice lattice with periodic boundary conditions (PBC), a finite spin-orbit coupling and a magnetic field, showing the topological nearly-flat bands with Chern number $C\!=\!\pm2$, close to the Fermi level ($E\!=\!0$). The topological superconducting phase is obtained by populating the lower topological flat band at the Fermi level.}
\label{dice_flatbands}
\vspace{-0mm}
\end{center}
\end{figure*}
%---------------------------------------------

The bipartite nature of the dice lattice, which can be envisaged as two merged triangular lattices, protects two degenerate flat bands coexisting with four dispersive bands. Such a geometry can be realized using a few layers of transition-metal oxides, dichalcogenides, and graphene. In the simplest realization of the dice lattice involving three (111) layers of cubic transition-metal oxides, such as in a SrTiO$_3$/SrIrO$_3$/SrTiO$_3$ trilayer, the cubic symmetry is reduced to trigonal symmetry. The strong spin-orbit coupling from the Ir$^{4+}$ ion and the broken inversion symmetry produces a Rashba spin-orbit coupling. In the reduced $D_{3d}$ symmetry of the trilayer, the Rashba spin-orbit coupling vectors lie in the plane parallel to the trilayer and have opposite senses of rotation for the top and the bottom layers of the three-coordination sites, surrounding the middle layer of six-coordination sites. In the presence of this Rashba spin-orbit coupling, the flat bands become isolated from the dispersive bands. Repulsive interactions in the flat bands then spontaneously generate ferro/ferri-magnetic order on the Kramer's pair of flat bands~\cite{PhysRevB.84.241103,PhysRevB.102.045105}, especially when they are close to half filling, and split them into two nearly-flat bands with Chern number $C\!=\! \pm2$, as shown in Fig.~\ref{dice_flatbands}{\bf {\small b}}. 

Besides the topological origin of superconductivity in the flat bands, superconductivity in the transition-metal oxide trilayer can also be induced by doping SrTiO$_3$~\cite{PhysRevLett.12.474,PhysRev.163.380}. Our calculations reveal that the superconducting state realized in the  dice-lattice topological flat bands, exhibits both singlet and triplet pairings. Such a singlet-triplet mixing is allowed by broken inversion symmetry in this oxide trilayer. Using symmetry analysis, we observed that the possible nearest-neihgbor pairing channels allowed in this lattice geometry are $d_{xy}$, $d_{x^2-y^2}$, $p_x$ and $p_y$. The nearest-neighbor pairing amplitudes were found to be complex numbers, illustrating the chiral nature of the induced topological superconducting phase. Besides finding corner MBS, supporting the second-order nature of the topological superconducting phase, we perform a Chern number analysis and obtain phase diagrams for the topological phase.

The dice lattice is a special case of the commonly-known $\alpha \!-\! \tau_{_3}$ lattice which interpolates between the dice lattice ($\alpha \!=\!0$, pseudospin 1) and the honeycomb lattice  ($\alpha \!=\!1$, pseudospin 1/2)~\cite{PhysRevB.94.125435}. It is also possible to transform the honeycomb lattice into the dice lattice, and vice versa, in an experimentally-simulated ultracold atomic gas platform~\cite{PhysRevB.73.144511}. Also, the results presented here are relevant to possible topological superconducting phases in twisted bilayer/multilayer graphene~\cite{PhysRevLett.121.087001}.\\

%\noindent  \textbf{\small RESULTS}\\
\section*{Results}
\noindent  \textbf{\small Model and set up}\\
The electron pairing in the topological flat bands of the dice lattice, realizable in a transition-metal oxide trilayer as discussed above, can be described by the following tight-binding Hamiltonian
\begin{align}
{\cal H}\!=\!&-t\! \sum_{\langle i\alpha,j\beta \rangle,\sigma} \!(c_{i\alpha\sigma}^{\dagger}c_{j\beta\sigma}^{\phantom{\dagger}}+{\rm H.c.})
\!-\mu\! \sum_{i,\alpha,\sigma} \! c_{i\alpha\sigma}^{\dagger}c_{i\alpha\sigma}^{\phantom{\dagger}} \nonumber \\
&-\lambda \! \sum_{\langle i\alpha,j\beta \rangle,\sigma,\sigma^{\prime}} \!(i [\mathbf{\hat{D}}_{ij}\cdot \boldsymbol{\sigma}]_{_{\alpha \beta}}^{\sigma \sigma^{\prime}}c_{i\alpha\sigma}^{\dagger}c_{j\beta\sigma^{\prime}}^{\phantom{\dagger}}+{\rm H.c.})  \nonumber \\
&-B_z \! \sum_{i,\alpha,\sigma,\sigma^{\prime}} \!( [\sigma_z]^{\sigma \sigma^{\prime}}c_{i\alpha\sigma}^{\dagger}c^{\phantom{\dagger}}_{i\alpha\sigma^{\prime}}+{\rm H.c.}) \nonumber \\
&-U_0 \sum_{i,\alpha} \!n_{i\alpha\uparrow}n_{i\alpha\downarrow}-\sum_{i,j,\alpha,\beta,\sigma,\sigma^{\prime}}U_1 n_{i\alpha\sigma}n_{j\beta\sigma^{\prime}},
\label{Ham}
\end{align}
%---------------------------------------------
\begin{figure*}[t]
\begin{center}
\vspace{-0mm}
\epsfig{file=Fig2.pdf,trim=0.0in 0.0in 0.0in 0.0in,clip=false, width=180mm}
\caption{{\bf Emergence of the corner MBS.} {\bf {\small a}} Real-space quasiparticle spectrum with varying chemical potential $\mu$, revealing the range $-1.16t \lesssim \mu \lesssim -0.82t$ (green region), where two pairs of low-energy eigenstates come close to zero energy, while other eigenstates move to higher energies, creating a topological energy gap that protects the Majorana bound states. {\bf {\small b}}, {\bf {\small c}} Plots of the local density of states $\rho (\bf{r})$ (in arbitrary units) in the topological superconducting phase ($\mu \!=\!-t$) and in the trivial superconducting phase ($\mu \!=\!-0.2t$). Other parameters are $\lambda \!=\!0.1t$ and $B_z\!=\!0.1t$. In {\bf {\small b}}, the localization of the MBS at the lattice corners indicates the second-order nature of the realized topological superconducting phase.}
\label{mu_var}
\vspace{-0mm}
\end{center}
\end{figure*}
%---------------------------------------------

\noindent where $t$ is the electron hopping amplitude, $i$ and $j$ are indices of different unit cells, $\alpha$ and $\beta$ represent indices of the three inequivalent sites within a unit cell, $\sigma=\uparrow,\downarrow$ labels the electron spin projection along the $z$ axis, $\langle \rangle$ represents nearest-neighbor sites, $\mu$ is the chemical potential, $\lambda$ is the strength of the Rashba spin-orbit coupling, $\mathbf{\hat{D}}_{ij}$ is the unit vector between unit cells $i$ and $j$, $\boldsymbol{\sigma}$ represents the Pauli matrices, $B_z$ is the strength of the magnetization field, the last two terms represent the onsite and non-local density-density attractive correlations with $U_0$ and $U_1$ as the strengths of the interactions, respectively, and $n_{i\alpha \sigma}=c_{i\alpha\sigma}^{\dagger}c_{i\alpha\sigma}$ is the electron density at the unit cell $i$, site $\alpha$ and spin $\sigma$. The density-density correlations are treated at the mean-field level and they give pairing amplitudes in different channels after decomposition. The pairing amplitudes between electrons at the same and nearest-neighbor (NN) sites are defined by $\Delta_{i\alpha\uparrow \downarrow}^{\rm onsite, s}=-U_0\langle c_{i\alpha\uparrow} c_{i\alpha\downarrow} \rangle$ (onsite singlet), $\Delta_{ij\alpha \beta \uparrow \downarrow}^{\rm NN, s}=-U_1\langle c_{i\alpha\uparrow} c_{j\beta\downarrow} \rangle$ (NN singlet), and $\Delta_{ij \alpha \beta\sigma \sigma}^{\rm NN, t}=-U_1\langle c_{i\alpha\sigma} c_{j\beta\sigma} \rangle$ (NN triplet). We choose $U_0$ and $U_1$ to be of the same value $U$; a different choice does not change the conclusions presented here because the pairing amplitudes are calculated self-consistently. For this self-consistent determination of the pairing amplitudes, we solve the Bogoliubov-de Gennes (BdG) equations, derived by performing the unitary transformation $c_{i\alpha\sigma}=\sum_{n}u_{i\alpha\sigma}^n\gamma_{n}^{\phantom{\dagger}}+v_{i\alpha\sigma}^{n*}\gamma_{n}^{\dagger}$ on the Hamiltonian~(\ref{Ham}), where $\gamma_{n}^{\phantom{\dagger}}$ is a fermionic annihilation operator in the $n^{\rm th}$ eigenstate, $u_{i\alpha\sigma}^n$ and $v_{i\alpha\sigma}^n$ are respectively the quasi-particle and quasi-hole amplitudes (see Methods for details). We use, throughout,  lattice spacing $a\!=\!1$, hopping scale $t\!=\!1$, and attractive potential $U\!=\!2t$. The triplet pairing amplitude is dynamically generated in the presence of Rashba spin-orbit coupling (broken inversion symmetry) and magnetic field (broken time-reversal symmetry)~\cite{PhysRevLett.87.037004,PhysRevB.97.214507}. Odd-parity, equal-spin pairing in the triplet channel is considered because such pairing is favored by parity fluctuations in the presence of Rashba spin-orbit coupling and a time-reversal symmetry-breaking Zeeman exchange field~\cite{PhysRevLett.115.207002}. Alternate routes to obtain spin-triplet topological pairing in similar systems include forward electron-phonon scattering which also suggests a robust equal-spin pairing~\cite{2022arXiv220709443L}.

\noindent  \textbf{\small \\Corner-localized MBS}\\
We investigate the emergence of the zero-energy MBS by inspecting the quasiparticle spectrum, obtained by numerically solving the Hamiltonian~(\ref{Ham}) while varying the chemical potential $\mu$. This procedure is repeated for many values of $\lambda$ and $B_z$, to search for signatures of the MBS in the quasiparticle spectrum. As shown in Fig.~\ref{mu_var}{\bf {\small a}}, at $(\lambda, B_z)\!=\!(0.1t, 0.1t)$ and  within the range $-1.16t \lesssim \mu \lesssim -0.82t$, two pairs of lowest-energy quasiparticle states remain close to zero energy while other low-energy levels move away towards higher energies, thus creating an energy gap. This energy gap provides topological protection to the near-zero-energy MBS preventing them from hybridizing with the higher-energy ordinary quasiparticle states. To study the real-space localization of these zero-energy MBS in a two-dimensional lattice, in Fig.~\ref{mu_var}{\bf {\small b}}-{\bf {\small c}} we plot the local density of states, obtained via $\rho(\mathbf{r})\!=\! \sum_{\alpha,\sigma}(|u_{i\alpha\sigma}^n|^2+v_{i\alpha\sigma}^n|^2)$ with $n$ corresponding to the lowest-positive energy eigenstate. We use two values for the chemical potential: $\mu \!=\!-t$, where the zero-energy states appear with a topological energy gap, and $\mu \!=\!-0.2t$, where the lowest-energy states are away from zero energy. At $\mu \!=\!-t$, the lowest-energy eigenstate is localized at the four lattice corners, while at $\mu \!=\!-0.2t$ it is distributed inside the bulk. The corner-localized zero-energy states provide a strong indication of the appearance of the MBS, and hence of the induced second-order topological superconducting phase.\\

\noindent  \textbf{\small Pairing symmetry}\\
The dice lattice has sites with coordination numbers three or six, and this feature distinguishes it from the triangular and hexagonal lattices. As discussed below, the presence of these two types of sites also determines the pairing symmetry in the superconducting state. Furthermore, the Rashba spin-orbit coupling enforces its symmetry in the superconducting pairing. 
%---------------------------------------------
\begin{figure}[t]
\begin{center}
\vspace{-0mm}
\epsfig{file=Fig3.pdf,trim=0.0in 0.0in 0.0in 0.0in,clip=false, width=85mm}
\caption{{\bf Pairing symmetry in the dice lattice.} Possible pairing symmetries around a six-coordination site in the dice lattice with Rashba spin-orbit coupling and induced magnetization. $\Delta$ is the pairing amplitude used for the illustration of the amplitudes along different neighbors.}
\label{symmetry}
\vspace{-4mm}
\end{center}
\end{figure}
%---------------------------------------------
In this $D_{3d}$ crystal symmetry with broken both inversion symmetry (due to Rashba spin-orbit coupling) and time-reversal symmetry (due to the induced magnetization), the possible pairing symmetries are \{$d_{xy}$, $d_{x^2-y^2}$\} in the singlet channel, and \{$p_x$, $p_y$\} in the triplet channel, as shown in Fig.~\ref{symmetry}. Mixing of the singlet and triplet components is allowed by the broken structural inversion symmetry in the discussed oxide trilayers~\cite{PhysRevLett.87.037004}. Therefore, a linear combination of these four types of pairing symmetry is stabilized. Figure~\ref{pairings} shows the profiles of the pairing amplitudes on the dice lattice at the same set of parameters where the corner MBS are found. The imaginary components of the nearest-neighbor pairing amplitudes are nonzero, implying a chiral mixed-parity topological superconducting state. The real part of the onsite singlet pairing amplitude Re($\Delta_{i\alpha\uparrow \downarrow}^{\rm onsite, s}$), in Fig.~\ref{pairings}{\bf {\small a}} top panel, clearly reveals a difference between the three and six coordination sites. While the onsite singlet pairing amplitude Re($\Delta_{i\alpha\uparrow \downarrow}^{\rm onsite, s}$) at the six-coordination sites is slightly smaller than that at the three-coordination sites, the nearest-neighbor singlet pairing amplitude Re($\Delta_{ij\alpha \beta \uparrow \downarrow}^{\rm NN, s}$), displayed in Fig.~\ref{pairings}{\bf {\small b}} top panel, shows the opposite behavior.
%---------------------------------------------
\onecolumngrid
\begin{figure}[b]
\begin{minipage}{18cm}
\begin{center}
\vspace{2mm}
\epsfig{file=Fig4.pdf,trim=0.0in 0.0in 0.0in 0.0in,clip=false, width=180mm}
\caption{{\bf Real-space profile of singlet and triplet pairing amplitudes.} Superconducting pairing amplitudes for different pairing channels: {\bf {\small a}} onsite singlet,  {\bf {\small b}} nearest-neighbor singlet, and {\bf {\small c}} nearest-neighbor triplet. Top and bottom rows show the real and imaginary parts of the pairing amplitudes. Model parameters are $\lambda \!=\!0.1t$, $B_z\!=\!0.1t$, and $\mu \!=\!-t$. The results indicate that the topological superconducting phase exhibits both singlet and triplet pairing components. The pairing amplitudes in the nearest-neighbor channel are complex numbers, implying a chiral pairing symmetry. OBC are used. The slightly different boundary effects among horizontal and tilted edges, as in the red vs white colored edges in lower panel {\bf a}, originate in the different site endings along those edges, as can be seen in Fig.~{\bf 1a}.}
\label{pairings}
\vspace{10mm}
\end{center}
\end{minipage}
\end{figure}
\twocolumngrid
%---------------------------------------------
%---------------------------------------------
\begin{figure*}[t]
\begin{center}
\vspace{5mm}
\epsfig{file=Fig5.pdf,trim=0.0in 0.0in 0.0in 0.0in,clip=false, width=170mm}
\caption{{\bf Quasiparticle band dispersion and topological phase diagram.} {\bf {\small a}}-{\bf {\small b}}  Bulk quasiparticle bands, obtained with periodic boundary conditions, along the momentum path $\Gamma (0,0)-{\rm K}(\frac{4\pi}{3\sqrt{3}},0)-{\rm M}(\frac{\pi}{\sqrt{3}},\frac{\pi}{3})-\Gamma$, below energy $E\!=\!0$, in the trivial ($\mu \!=\! -0.2t$) and in the topological ($\mu \!=\! -t$) superconducting phases. In {\bf {\small b}}, one band (in thick green) is isolated and carries a Chern number $C\!=\!2$. {\bf {\small c}}-{\bf {\small d}} Quasiparticle spectrum, obtained in a cylindrical geometry, at the same sets of parameters as in {\bf {\small a}}, {\bf {\small b}}, respectively. The pair of lowest-energy edge states (in red) crosses in the case of the topological superconducting phase in {\bf {\small d}}, near $k_1\!=\!0.95\pi$. Parameters used are $\lambda \!=\!0.1t$, $B_z\!=\!0.1t$, $\Delta^{\rm s}_{1,3}\!=\!0.5t$, and $\Delta^{\rm s}_{2}\!=\!0.45t$. {\bf {\small e}}-{\bf {\small f}} Plots of the Chern number of the occupied bands in the $\mu \!-\!B_z$ plane at two Rashba coupling strengths $\lambda \!=\!0.1t$ and $\lambda \!=\!0.3t$.}
\label{kspace_bands}
\vspace{-0mm}
\end{center}
\end{figure*}
%---------------------------------------------
The slight variation in the pairing amplitudes near the corners and edges of the lattice is due to the open boundary conditions used in the calculations. The real part of the nearest-neighbor triplet pairing amplitude Re($\Delta_{ij \alpha \beta\sigma \sigma}^{\rm NN, t}$) vanishes inside the bulk, but remains finite at the boundaries along the $\pm x$ directions. Such an exact cancellation of the pairing amplitude inside the bulk implies a $p_x$-wave pairing component. The imaginary part of the nearest-neighbor triplet pairing amplitude Im($\Delta_{ij \alpha \beta\sigma \sigma}^{\rm NN, t}$) is finite in the bulk and it is different in magnitude and changes its sign at the six and three coordination sites. This implies a $p_y$-wave pairing component because only a partial cancellation of the pairing amplitudes is possible along the $y$ direction in the dice lattice due to the lattice asymmetry along the $x$ and $y$ directions.\\

\noindent  \textbf{\small Topological superconducting transition}\\
The transition to the second-order topological superconducting phase can be understood by inspecting the quasiparticle band dispersion in momentum space. The realized topological superconducting phase has broken time-reversal symmetry, but it respects particle-hole symmetry, and hence it belongs to the Altland-Zirnbauer symmetry class D and can be identified by evaluating the Chern number~\cite{PhysRevB.78.195125}. We calculate the topological invariant of the superconducting phase as the sum of the Chern numbers of the quasiparticle bands below $E\!=\!0$. Isolated quasiparticle bands in the dice lattice have Chern number $C\!=\!\pm2$. Therefore, the Chern number of the total occupied states becomes finite whenever an odd number of topological quasiparticle bands appear as isolated at a given set of parameters. The momentum-space quasiparticle bands are obtained by diagonalizing the following BdG Hamiltonian at 
wavevector $\mathbf{k}$
\begin{align}
{\cal H}_{\rm BdG}(\mathbf{k})\!=\!
\begin{pmatrix}
\Psi_{\mathbf{k}}^{\dagger} & \Psi_{-\mathbf{k}}
\end{pmatrix}
\begin{pmatrix}
H_e   & H_{\Delta} \\
H_{\Delta}^* & H_h
\end{pmatrix}
\begin{pmatrix}
\Psi_{\mathbf{k}} \\
\Psi_{-\mathbf{k}}^{\dagger}
\end{pmatrix}
\label{H_k}
\end{align}
where $H_e$, $H_h$ and $H_{\Delta}$ are the matrices representing, respectively, the electron, the hole, and the pairing sectors of the Hamiltonian. The basis wave vector for $H_e$ is $\Psi_{\mathbf{k}}=[c_{\mathbf{k}1\uparrow}, c_{\mathbf{k}2\uparrow}, c_{\mathbf{k}3\uparrow},c_{\mathbf{k}1\downarrow}, c_{\mathbf{k}2\downarrow}, c_{\mathbf{k}3\downarrow}]^T$ (see Methods section for more information). Figures~\ref{kspace_bands}{\small {\bf a}}-{\small {\bf b}} show the quasiparticle band dispersion at two values of the chemical potential $\mu \!=\! -0.2t$ and $\mu \!=\! -t$. In the latter, an isolated topological band with $C\!=\!2$ appears, resulting in a second-order topological phase. In Fig.~\ref{kspace_bands}{\small {\bf c}}-{\small {\bf d}}, we plot the energy spectra using now a cylindrical geometry, with periodic boundary condition along  the $x$ direction and open boundary condition along  the $y$ direction, for the same two values of $\mu$ as above. In the trivial superconducting phase, Fig.~\ref{kspace_bands}{\small {\bf c}}, an energy gap prevails in the spectrum, while in the case of the topological superconducting phase, Fig.~\ref{kspace_bands}{\small {\bf d}}, two pairs of doubly-degenerate low energy quasiparticle bands cross near $k_1\!=\!0.95\pi$. The Chern number, used to construct a phase diagram for the second-order topological superconducting phase in the $\mu \!-\!B_z$ plane, is shown in Fig.~\ref{kspace_bands}{\small {\bf e}}-{\small {\bf f}}, for two values of the Rashba spin-orbit coupling strengths $\lambda \!=\!0.1t$ and $\lambda \!=\!0.3t$. The induced topological phase appears within a range of $\mu$, and this range tends to decrease with increasing the Zeeman field or the Rashba coupling.

%\noindent  \textbf{\small \\DISCUSSION}\\
\section*{Discussion}
Analogies between the topological superconductivity in flat bands, as found here, and the quantum-Hall insulator/superconductor interfaces can be drawn. Theoretically, it is known that a quantum Hall state with Chern number $C\!=\!1$, in proximity to a fully gapped $s$-wave superconductor, generates a topological (first-order) superconducting phase~\cite{PhysRevB.82.184516,PhysRevResearch.2.023063}. Likewise, the fractionalized MBS, \textit{i.e.} some realizations of parafermions, have been proposed in fractional quantum Hall states when in proximity to an $s$-wave superconductor~\cite{Clarke2013,PhysRevX.4.031009,PhysRevResearch.2.013232}. In the $C\!=\!2$ flat band considered here with an attractive pairing, the topological superconductivity is of order two. These findings suggest a close connection between the Chern number $C$ of the normal state and the degree of the induced topological superconductivity. 

To summarize, we showed that topological flat bands in the dice lattice with attractive interaction among electrons harbor second-order topological superconductivity. A signature of this exotic topological phase is revealed by the presence of MBS at the lattice corners. We propose that topological flat-bands provide a promising platform for realizing a general $n^{\rm th}$-order topological superconductor, where $n$ is the Chern number of the flat band. Topological flat bands with higher Chern numbers are found not only in the dice lattice, but also in kagom{\'e} and Lieb lattices~\cite{PhysRevB.104.235115,Okamoto_CommPhys2022}. Other than the examples of a few-layer graphene and a transition-metal-oxide trilayer, another candidate compound is CsV$_3$Sb$_5$~\cite{PhysRevLett.125.247002, Zhao2021, HU2022495}, where lattice geometry, flat-band topology and superconductivity may also produce Majorana states such as those discussed here. Hence, we expect future research will unveil topological superconductivity in a variety of compounds that exhibit topological flat bands. Furthermore, the superconducting transition temperature is proportional to the density of states at the Fermi level which is large for these flat-band systems. Therefore, when looking forward, topological flat-bands with higher Chern numbers provide an opportunity to search for higher-order topological superconductivity at high temperatures.\\

%\noindent  \textbf{\small METHODS}\\
\section*{Methods}
\noindent  \textbf{\small Self-consistent BdG formalism\\}
The attractive interaction term in the Hamiltonian~(\ref{Ham}) is decoupled into different pairing channels (singlet and triplet, onsite and nearest-neighbor). The mean-field Hamiltonian is then diagonalized using the BdG transformation $c_{i\alpha\sigma}\!=\!\sum_{n}u_{i\alpha\sigma}^n\gamma_{n}+v_{i\alpha\sigma}^{n*}\gamma_{n}^{\dagger}$, where $\gamma_{n}$ is a fermionic annihilation operator at the $n^{\rm th}$ eigenstate, $u_{i\alpha\sigma}^n$ and $v_{i\alpha\sigma}^n$ are  the quasi-particle and quasi-hole amplitudes, respectively. The quasi-particle and quasi-hole amplitudes are obtained by solving the BdG equations $\sum_{j}{\cal H}_{ij}\psi_{j}^n \!=\!\epsilon_n\psi_{i}^n$, where $\psi_{i}^n\!=\![u_{i\alpha\uparrow}^n, u_{i1\downarrow}^n, v_{i\alpha\uparrow}^n, v_{i\alpha\downarrow}^n]^T$ with $u_{i\alpha\uparrow}^n=[u_{i1\uparrow}^n,u_{i2\uparrow}^n,u_{i3\uparrow}^n]$ and similarly for other components, while $\epsilon_n$ is the energy eigenvalue of the $n^{\rm th}$ eigenstate. Then, the pairing amplitudes are calculated via $\Delta_{i\alpha\uparrow \downarrow}^{\rm onsite, s}\!=\!-U_0\langle c_{i\alpha\uparrow} c_{i\alpha\downarrow} \rangle$ (onsite singlet), $\Delta_{ij\alpha \beta \uparrow \downarrow}^{\rm NN, s}\!=\!-U_1\langle c_{i\alpha\uparrow} c_{j\beta\downarrow} \rangle$ (NN singlet), and $\Delta_{ij \alpha \beta\sigma \sigma}^{\rm NN, t}\!=\!-U_1\langle c_{i\alpha\sigma} c_{j\beta\sigma} \rangle$ (NN triplet). The self-consistency iterations continue until all the pairing amplitudes converge at all lattice sites, within a tolerance of 10$^{-8}$.

\noindent  \textbf{\small \\Momentum-space Hamiltonian\\}
The Hamiltonian~(\ref{H_k}) is expressed in the basis $\Psi_{\mathbf{k}}=[c_{\mathbf{k}1\uparrow}, c_{\mathbf{k}2\uparrow}, c_{\mathbf{k}3\uparrow},c_{\mathbf{k}1\downarrow}, c_{\mathbf{k}2\downarrow}, c_{\mathbf{k}3\downarrow}]^T$, where 1, 2, 3 denote the three inequivalent sites within a unit cell, and is given by
\begin{align}
&{\cal H}_{e}(\mathbf{k})\!=\! \nonumber \\
&\begin{pmatrix}
-B_z-\mu  & -t\gamma_{\mathbf{k}}^* & 0 & 0 &  -i\lambda \gamma_{\mathbf{k+}}^* & 0\\
-t\gamma_{\mathbf{k}}  & -B_z-\mu & -t\gamma_{\mathbf{k}}^* & i\lambda \gamma_{\mathbf{k-}} & 0 &  i\lambda \gamma_{\mathbf{k+}}^*\\
0  & -t\gamma_{\mathbf{k}} & -B_z-\mu & 0 &  -i\lambda\gamma_{\mathbf{k-}} & 0\\
0  & -i\lambda\gamma_{\mathbf{k-}}^* & 0 & B_z-\mu &  -t\gamma_{\mathbf{k}}^* & 0\\
i\lambda \gamma_{\mathbf{k+}}  & 0 & i\lambda\gamma_{\mathbf{k-}}^* & -t\gamma_{\mathbf{k}} & B_z-\mu &-t\gamma_{\mathbf{k}}^*\\
0  & -i\lambda\gamma_{\mathbf{k+}} & 0 & 0 &  -t\gamma_{\mathbf{k}} & B_z-\mu
\end{pmatrix},
\end{align}
where $\gamma_{\mathbf{k}}\!=\!1\!+\!e^{ik_1}\!+\!e^{ik_2}$, $\gamma_{\mathbf{k\pm}}\!=\!1+e^{ik_1\pm2\pi/3}+e^{ik_2\pm4\pi/3}$, $k_{1,2}\!=\!\mathbf{k}.\cdot \mathbf{e}_{1,2}$, and $\mathbf{e}_{1,2}$ are lattice translational vectors. The topological flat bands in Fig.~\ref{dice_flatbands}{\small{\bf B}} are obtained by diagonalizing ${\cal H}_{e}(\mathbf{k})$ at $t\!=\!1$, $\mu \!=\!0$, $\lambda \!=\!0.3t$, and $B_z\!=\!0.4t$. The hole part of the Hamiltonian~(\ref{H_k}) is given by ${\cal H}_{h}(\mathbf{k})\!=\![-{\cal H}_{e}(-\mathbf{k})]^T$, and the pairing part is given by
\begin{align}
{\cal H}_{\Delta}(\mathbf{k})\!=\!
\begin{pmatrix}
0  &  0 & 0 & \Delta_{1}^{\rm s} & 0  & 0\\
0  & 0 & 0 & 0 & \Delta_{2}^{\rm s}  & 0\\
0  & 0 & 0 & 0 & 0 & \Delta_{3}^{\rm s} \\
-\Delta_{1}^{\rm s} & 0 & 0 & 0 & 0 & 0\\
0 & -\Delta_{2}^{\rm s}  & 0 & 0 & 0 & 0\\
0 & 0 & -\Delta_{3}^{\rm s} & 0 & 0 & 0 \\
\end{pmatrix},
\end{align}
where $\Delta_{\alpha}^{\rm s}$ ($\alpha \!=\!1,2,3$) represents the onsite singlet pairing amplitude at site index $\alpha$.

%\noindent {\bf \\Acknowledgements:} 
\section*{Acknowledgements}
This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

\section*{Data availability}
All data obtained from numerical calculations have been presented in the paper.

\section*{Code availability}
Simulation codes are available from the corresponding author upon reasonable request.

\section*{Author contributions}
NM planned the work, performed numerical calculations and wrote the manuscript with inputs from all coauthors. All authors contributed to the analysis of the results, data presentation and manuscript writing.

\section*{Competing interests}
The authors declare no competing interests.

\section*{Correspondence}
Correspondence should be addressed to NM.

\bibliographystyle{naturemag}
\bibliography{Ref}



\end{document} 

