\section{Introduction}
\label{sec:intro}

The discovery of electronic insulators with topologically nontrivial band structures has led to remarkable progress in understanding gapped quantum phases. The prediction and experimental discovery of topological insulators (TIs)~\cite{kanemele,bernevigzhang,fukanemele,moorebalents,roy2009,konig2007,hsieh2008} and topological superconductors (TSCs)~\cite{readgreen,ivanov,stoneroy,zhang2018} led to a classification of gapped phases of non-interacting fermions~\cite{ryu2010,kitaev2009}; this ten-fold way encodes whether a system may host topologically nontrivial phases given the spatial dimension and the symmetries under which it is invariant. The nontrivial band topology of electronic states is manifest in striking universal properties, including robust gapless modes confined to the sample boundary and quantized response coefficients~\cite{hasankanermp,qizhangrmp}.

These concepts were later extended to crystalline symmetries, such as reflection, inversion, or rotation. Gapped phases protected by these symmetries are called topological crystalline insulators (TCIs)~\cite{fuTCI,hsiehTCI,okada2013,sessi2016,ma2017} and include higher-order topological insulators (HOTIs)~\cite{schindlerHOTI,benalcazar2017,langbehnHOTI,songHOTI,khalaf2018}. Specifically, an $n^{th}$ order TI/TSC in $d$ spatial dimensions is gapped everywhere except on a $d-n$ dimensional surface. More generally, TI/TSCs and HOTIs are examples of symmetry-protected topological (SPT)~\cite{chen2013,guwen2014} and crystalline SPT (cSPT)~\cite{song2017,buildingblock,elsethorngren} phases respectively, whose classification also accounts for interactions. Such phases have a trivial gapped bulk but host boundary (or hinge/corner) modes protected against local, symmetry-preserving perturbations~\cite{senthilSPT}.

A defining feature of topological phases is their robustness against disorder: provided the spectral (or mobility) gap remains finite and the disorder respects the symmetry protecting the TI/TSC, quantized coefficients and gapless edge modes persist~\cite{mong2012,ringel2012,schubert,foster2014}. Despite disorder breaking the lattice symmetries protecting TCIs, their boundary modes can evade localization when the full ensemble of disorder configurations remains symmetric~\cite{fulga2014,prodan2015}. Traditionally, robustness of topological states is established by adding disorder to a clean system, thereby assuming an underlying periodic reference state. This approach, while efficacious, fails when no such structure exists \textit{i.e.,} for aperiodic systems, including amorphous, quasiperiodic, and fractal systems. Nonetheless, topological phenomena have been shown to exist in both amorphous~\cite{agarwala2017,shibaglass,mitchell2018,bourne2018,agarwala2019,ojanen2019} and quasiperiodic~\cite{kraus1,kraus2,loring2016,bandres2016,quasicrystalQSH,spinbott} systems.
\begin{figure}[b]
\centering
\begin{subfigure}[t]{0.23\textwidth}
\includegraphics[width=\textwidth]{lattice.jpg}
\caption{}
\label{fig:gasket}
\end{subfigure}\quad
\begin{subfigure}[t]{0.23\textwidth}
\includegraphics[width=\textwidth]{realspaceChernNumber.jpg}
\caption{}
\label{fig:partition}
\end{subfigure}
\caption{(a) SG with ``periodic'' boundary conditions, such that all sites have coordination number four. (b) Regions $A,B,C$ considered in real-space Chern number calculations.}
\end{figure}

That the topology of quantum states can be defined in the absence of spatial regularity over long distances opens the door to finding topological phases on fractal lattices, which lack a natural distinction between bulk and boundary, and whose (typically non-integer) Hausdorff dimensions differ from their topological dimensions. Interest in fractal structures, which have a rich history~\cite{recursion,banavar,gefen1,gefen2,rammal,alexander}, has been revived given experimental advances in creating and manipulating synthetic lattices with arbitrary structures, in both photonic and electronic systems \cite{kollarhyper,lieblattice,pband,collins2017,girovsky,drost}. In particular, fractal lattices have been fabricated using focused ion beam milling~\cite{herek}, molecular chains~\cite{shang,tait2015,zhang2016}, and scanning-tunneling-microscopy (STM) techniques~\cite{kempkes2018}, with theoretical studies primarily focusing on localization and transport phenomena~\cite{vanVeen1,vanVeen2,pal2012,kosior2017,bpal2018,genzor2019}.

 \begin{figure*}[t]
 \centering
 \includegraphics[width=0.8\textwidth]{mainfig2.jpg}
 \caption{Gapless topological phase of the chiral $p+ip$ superconductor on the SG, with $g=5,\Delta=1,t=0.5, \mu=0.5$. Energy spectrum and probability densities of eigenvectors at indicated energies are shown. Color scale indicates values of $x,y,z$ coordinates, with dot size indicating the magnitude of the probability density at that point.}
 \label{fig:spectrum}
 \end{figure*}

However, our understanding of the influence of self-similar geometry on the topological character of electronic states remains nascent, having received attention only recently~\cite{fractalized,neupert}. In this paper, we fill this lacuna by developing a general framework elucidating the fate of topological states on fractal lattices embedded in two dimensions (2D). Through this picture, we find that the nature of thermodynamic phases---gapped vs gapless---on fractal lattices depends crucially on the ratio of bulk to edge coordinated sites. Focusing on the chiral $p$-wave superconductor on the Sierpinski gasket, we show that qualitative features obtained through numerical diagonalization can be understood simply through our framework. Besides characterizing the two distinct phases of this model using various theoretical tools, we further corroborate our understanding by studying both the $p$-wave superconductor and an HOTI on the Sierpinski carpet.

\section{Model}
\label{sec:model}

\begin{figure*}[t]
 \centering
 \includegraphics[width=0.8\textwidth]{mainfigTrivial2.jpg}
 \caption{Gapped trivial phase of the $p+ip$ superconductor on the SG, with $g=5,\Delta=1,t=0.5, \mu=2$. Energy spectrum and probability densities of eigenvectors at indicated energies are shown. Color scale indicates values of $x,y,z$ coordinates, with dot size indicating the magnitude of the probability density at that point.}
 \label{fig:spectrumTrivial}
 \end{figure*}

We consider a 2D spinless chiral $p$-wave superconductor (symmetry class D~\cite{ryu2010}) within the Bogoliubov-deGennes (BdG) framework, with the mean field lattice BCS Hamiltonian:
\begin{equation}
\label{bdgHam}
\hat{H} = -t \sum_{\langle \bf{r}, \bf{r'} \rangle} \hat{c}_{\bf{r}}^\dagger \hat{c}_{\bf{r'}} - \mu \sum_{\bf{r}}  \hat{c}_{\bf{r}}^\dagger \hat{c}_{\bf{r}} + \sum_{\bm{r},m} [ \Delta_m \hat{c}_{\bf{r} + \bm{e}_m}^\dagger \hat{c}_{\bf{r}}^\dagger + \text{h.c.} ],
\end{equation}
where $\hat{c}_{\bf{r}}^\dagger,\hat{c}_{\bf{r}}$ satisfy fermionic anti-commutation relations $\{\hat{c}_{\bf{r}},\hat{c}_{\bf{r'}}^\dagger \} = \delta_{\bf{r},\bf{r'}}$, $t$ is the nearest-neighbor hopping, $\mu$
is the chemical potential, and we set the lattice spacing $a = 1$. Specifying to a triangular lattice~\cite{cheng2010}, the pairing term $\Delta_m = \Delta e^{i \pi m/3}$ is defined on the nearest-neighbor bonds corresponding to the three lattice vectors $\bf{e}_m$ with azimuthal angles $m \pi/3$ ($m=0,1,2$). We introduce the standard Bogoliubov transformation: $\hat{c}_{\bf{r}} = \sum_{\bf{r}} \left[u_{n,\bf{r}} \hat{\gamma}_n + v_{n,\bf{r}} \hat{\gamma}_n^\dagger \right]$, where $\hat{\gamma}_n$ is the Bogoliubov quasiparticle annihilation operator and $(u_{n,\bf{r}}, v_{n,\bf{r}})^T$ diagonalizes the BdG Hamiltonian~\eqref{bdgHam}, with eigenvalue $E_n$.

We study this model on a Sierpinski gasket (SG) with ``periodic'' boundary conditions (see Fig.~\ref{fig:gasket}) \textit{i.e.,} with four gaskets arranged on alternating faces of an octahedron, ensuring that all lattice sites are equally (four) coordinated. We construct a lattice regulated (with a smallest triangle) SG recursively, by adding sites/bonds to a gasket at generation $g$ to arrive at the $g+1$ SG. The largest lattice we can probe numerically has $g=6$, with the total number of sites $N \sim 3^{g+1}$ at generation $g$.

Setting $\Delta > 0$ and noting that the Hamiltonian Eq.~\eqref{bdgHam} admits a topological phase on a triangular lattice for $-6t < \mu < 2t$ (see Appendix \ref{app:triangle}), we find that this model admits topologically distinct phases even on the SG. The qualitative distinction between the two phases is illustrated in Fig.~\ref{fig:spectrum}, which shows the spectrum and states for $g = 5$. For $\mu>2t$ or $\mu<-6t$, we find a fully gapped ``trivial'' phase (see Fig. \ref{fig:spectrumTrivial}), where eigenstates are delocalized, thereby behaving as bulk states in ordinary gapped systems. In the thermodynamic ($g \to \infty$) limit, the spectrum is self-similar, with infinitely many gaps. In contrast, for $-6t<\mu<2t$ we find that the amplitude of the largest gap in the spectrum decays exponentially with increasing generation (see Appendix \ref{gap:scaling}), such that the spectrum is strictly \textit{gapless} in the $g\to \infty$ limit. Thus, this parameter range describes a qualitatively distinct phase with emergent \textit{continuous} scale invariance, unlike the trivial phase which only possesses \textit{discrete} scale invariance. Particle-hole symmetry is present in both phases. While the spectra are obtained by numerically diagonalizing the BdG Hamiltonian~\eqref{bdgHam}, these can in principle also be obtained recursively (see Appendix \ref{app:recursive} for details).

An intriguing feature of the gapless phase is the \textit{edge-like} nature of eigenstates: in Fig.~\ref{fig:spectrum}, we plot the electronic densities for representative states at the indicated energies, revealing states sharply localized on triangular motifs formed by sites of various generations \textit{i.e.,} localized around the inner edges (or holes) of the SG. While states closest to $E = 0$ are localized on the outer edges, corresponding to the earliest generations, there is a hierarchy of states localized on inner edges created at subsequent generations of the SG. In the thermodynamic limit, we expect that \textit{all} eigenstates in this phase will be sharply localized along edges. Remarkably, these localized states are also \textit{chiral}, with a wave-packet initialized on any inner edge propagating in the direction opposite to that of one initialized on the outermost edge (see Appendix \ref{app:chiral}).

Surprisingly, we find that the transition between the trivially gapped and the gapless phase coincides with the trivial $\leftrightarrow$ topological transition of Eq.~\eqref{bdgHam} on the triangular lattice. This observation hints that the model on the SG inherits its behavior from one defined on a triangular lattice. Indeed, we can regard the inner edges of the SG as holes in a triangular lattice, which, in the topological phase of Eq.~\eqref{bdgHam}, host gapless chiral Majorana modes propagating counter to the outermost edge state~\cite{aliceareview}. Since the number of these holes increases with $g$, there are infinitely many gapless modes in the spectrum as $g\to \infty$, resulting in a gapless spectrum. This physical picture suggests that the chiral eigenstates in the gapless phase are descended from Majorana edge modes of the $p+ip$ state on a triangular lattice. We hence dub this the gapless \textit{topological} phase on the SG.

\begin{figure*}[t]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\includegraphics[width=\textwidth]{chernNoTrivPhaseGen4tp5dp5mu2.jpg}
\caption{}
\label{fig:cherntriv}
\end{subfigure}\quad
\begin{subfigure}[t]{0.45\textwidth}
\includegraphics[width=\textwidth]{chernNoTopPhaseGen4tp5d1mup5.jpg}
\caption{}
\label{fig:chernTop}
\end{subfigure}
\caption{The real-space Chern number (black curve) as a function of Fermi energy $E_{f}$, with the corresponding spectrum shown in red in (a) the trivial phase ($\mu=2$), and (b) the topological phase ($\mu=0.5$). Here, $g=4$, $t=0.5$, and $\Delta=0.5$. }
\end{figure*}

\section{Diagnostics}
\label{sec:diagnostics}

Before building on this intuitive picture and showing that it generalizes to other fractal lattices, such as the Sierpinski carpet (SC), and other topological states, we further characterize the two distinct phases of the $p+ip$ superconductor on the SG using some standard diagnostics.

\subsection{Real-space Chern number}
\label{ssec:realSpaceChernNo}

Since our model lacks translation invariance, and only retains (discrete) scale invariance, we cannot use the momentum-space Chern number to characterize the topological and trivial phases of the $p+ip$ superconductor on the SG. Thus, we instead compute the real-space Chern number introduced in Ref.~[\onlinecite{kitaev2006}], which reduces to the momentum-space Chern number in the presence of translation invariance:
\begin{equation}
    \mathcal{C}= 12 \pi i \sum_{j\in A}\,\sum_{k\in B}\,\sum_{l\in C}(\mathcal{P}_{jk}\mathcal{P}_{kl}\mathcal{P}_{lj} - \mathcal{P}_{jl}\mathcal{P}_{lk}\mathcal{P}_{kj}),
\label{realspace}
\end{equation}
where $\mathcal{P}$ projects onto occupied states with respect to a given chemical potential, and $j,k,l$ are indices corresponding to three distinct neighboring regions $A,B,C$, arranged counter-clockwise (see Fig.~\ref{fig:partition}). In Eq. \ref{realspace}, $\mathcal{P}_{ij}$ is a $2\times 2$ matrix whose rows correspond to $c^{\dagger}_{i}, c_{i}$, and whose columns correspond to  $c_{j}, c^{\dagger}_{j}$. Retaining the site basis, we rotate only the $k$ /pseudospin basis. We then diagonalize the $2 \times 2$ matrix in the expression for $\mathcal{C}$ such that pseudospin is now a good quantum number, and then take the trace. With $\widetilde{\mathcal{P}}_{ij}$ representing the $2 \times 2$ block after diagonalization, the expression for the Chern number can be rewritten as:

\begin{equation}
    \mathcal{C}= 12 \pi i \sum_{j\in A}\,\sum_{k\in B}\,\sum_{l\in C}\textrm{Tr}(\widetilde{\mathcal{P}}_{jk}\widetilde{\mathcal{P}}_{kl}\widetilde{\mathcal{P}}_{lj} - \widetilde{\mathcal{P}}_{jl}\widetilde{\mathcal{P}}_{lk}\widetilde{\mathcal{P}}_{kj}),
\label{chernTrace}
\end{equation}

For $g=5$ in the trivial phase, we find that $\mathcal{C} = 0$ for all gapped regions of the spectrum (see Fig. \ref{fig:cherntriv}). We have checked that this quantization becomes independent of the specific choice of regions $A,B,C$ at large $g \geq 4$ \textit{i.e.,} in the limit when the number of sites in each region becomes large. In the thermodynamic ($g \to \infty$) limit, the spectrum within the trivial phase displays an infinite hierarchy of self-similar gaps, and we expect that $\mathcal{C}$ will vanish identically for each of the infinitely many gaps in the spectrum.

In contrast, within the topological phase the gapped regions of the SG scale to zero and have a trivially quantized Chern number. As can be seen in Fig.~\ref{fig:chernTop}, we find that indeed $\mathcal{C}=0$ within the finite-size gaps at finite $g$ in the topological phase. Nevertheless, similarly to previous works on topological amorphous superconductors~\cite{ojanen2019} and on the quantum Hall effect on fractal lattices~\cite{neupert}, we expect the Chern number to take a non-trivial quantized value within the gapless regions due to the presence of a mobility gap and the topological nature of the phase. While our numerics suggest that the Chern number tends towards a quantized non-zero value with increasing $g$ in regions corresponding to low but non-zero density of states, we are numerically limited to $g \leq 5$, for which finite-size effects obscure the expected quantization.

Thus, in the thermodynamic limit, the trivial phase will exhibit a strictly quantized $\mathcal{C} = 0$ within the infinitely many gaps in the spectrum; on the other hand, although the spectrum becomes gapless in the topological phase, we expect that $\mathcal{C}$ converges to a non-trivial quantized value as $g \to \infty$ in the gapless regions due to the presence of a mobility gap~\cite{neupert}. Verifying the latter requires investigating the model on a SG with large $g$, which is beyond our current numerical capabilities.

\begin{figure*}[t]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\includegraphics[width=\textwidth]{IPRtrivial.jpg}
\caption{}
\label{fig:iprtriv}
\end{subfigure}\quad
\begin{subfigure}[t]{0.45\textwidth}
\includegraphics[width=\textwidth]{IPRtopological.jpg}
\caption{}
\label{fig:iprtop}
\end{subfigure}
\caption{IPR$_n$ (using open boundary conditions on a single SG), with $g = 4, t = 0.5, \Delta = 0.5$. All states are delocalized in the (a) trivial phase ($\mu = 2$) while the (b) topological phase ($\mu = 0.5$) exhibits states localized around edges.}
\end{figure*}

\begin{figure}[b]
 \centering
 \includegraphics[width=0.4\textwidth]{g5IPRtopological.jpg}
 \caption{IPR$_{n}$ with $g=5, t=0.5, \Delta=0.5$. Comparison with Fig.~\ref{fig:iprtop} clearly shows that the number of localized states in the middle of the spectrum increases with $g$.}
 \label{fig:iprtopg5}
 \end{figure}

\subsection{Inverse participation ratio}
\label{ssec:ipr}

\begin{figure}[b]
 \centering
 \includegraphics[width=0.4\textwidth]{levelStatistics.jpg}
 \caption{Distribution of normalized energy level spacings with disorder $W$, with $g = 4 , t = 0.5, \Delta = 0.5$. Level statistics shown for weak ($W=2$) and strong ($W=8$) disorder for the trivial ($\mu = 2$) and topological ($\mu = 0.5$) phase.}
 \label{fig:levelstatistics}
 \end{figure}
Another useful diagnostic is the inverse-participation-ratio (IPR) of the $n^{th}$ eigenstate~\cite{thouless1974,franz1996},
\begin{equation}
\textrm{IPR}_n = \frac{\sum_{\bf{r}} \left( |u_{n,\bf{r}}|^4 + |v_{n,\bf{r}}|^4 \right)  } { \big[ \sum_{\bf{r}}  \left( |u_{n,\bf{r}}|^2 + |v_{n,\bf{r}}|^2 \right) \big]^2}
\end{equation}
which scales as $L^{-2}$ for extended states but remains finite for localized states even in the thermodynamic limit. In the trivial phase, \textit{all} eigenstates are delocalized (see Fig.~\ref{fig:iprtriv}), reflecting their bulk nature. Increasing $g$ suppresses the IPR values further towards zero. In the topological phase, the IPR values instead abruptly jump between $\sim 0$ and $\sim 1$, with the latter corresponding to eigenstates localized along the various edges (or holes) of the SG, as in Fig.~\ref{fig:spectrum}. The number of localized states increases with $g$ (see Figs. \ref{fig:iprtop} and \ref{fig:iprtopg5}), consistent with the physical picture discussed above: cutting out holes from the triangular lattice does not introduce any edge modes in the trivial phase, and all states remain extended. In the topological phase however, additional gapless edge modes are introduced, with the number of such modes increasing with $g$. This agrees with the numerical observation of localized states with IPR$_n \sim 1$ as shown in Fig.~\ref{fig:iprtop}.

\subsection{Level Statistics}

Since disorder provides an independent probe of topology, we add an onsite term $\sum_{\bf{r}} V_{\bf{r}} \hat{c}_{\bf{r}}^\dagger \hat{c}_{\bf{r}}$ to the Hamiltonian~\eqref{bdgHam}, with $V_{\bf{r}}$ drawn randomly from the uniform distribution $[-W/2,W/2]$. In Fig.~\ref{fig:levelstatistics}, we plot the energy-level spacing distributions, averaged over 500 disorder realizations, for weak and strong disorder in both phases. The normalized level spacing is given by $s = |E_n - E_{n+1}|/\delta(E_n)$, with $\delta(E_n)$ the mean-level spacing near energy $E_n$. In the trivial phase, the distribution is Poissonian at both weak and strong disorder, consistent with a \textit{localized} phase. The level spacings in the topological phase follow unitary Wigner-Dyson (GUE) statistics at weak disorder ($W=2$) and transition to Poisson at strong ($W=8$) disorder, with the transition\footnote{Since we expect regions with quantized $\mathcal{C} \neq 0$ in the gapless topological phase (due to a mobility gap), the transition likely occurs through a critical delocalized state, as is the case for insulating systems with a quantized $\mathcal{C} \neq 0$ at small disorder~\cite{prodan2010}.} to the Anderson insulator occurring at $W \sim 5$. Agreement with the Wigner surmise (for $\beta =2$) at weak disorder (see Fig.~\ref{fig:levelstatistics}) indicates that the gapless topological phase is a \textit{diffusive metal}~\cite{mirlin2000,chou2014}.

\section{Recursive Decimation}
\label{sec:decimation}

We propose a physical picture which elucidates how topological states on 2D fractal lattices inherit their behavior from a ``parent'' state on an underlying periodic lattice. Consider the BdG Hamiltonian~\eqref{bdgHam} on a triangular lattice with open boundary conditions, lattice spacing $a$, and size $L=2^p a$, as shown in Fig.~\ref{fig:decimation}. We define bulk and edge sites as those with coordination number six and four, respectively\footnote{Four copies of the lattice can be arranged as in Fig.~\ref{fig:gasket} such that the corner sites are also four coordinated.}. We now decimate sites and bonds recursively to generate the SG. At the $g^{th}$ step ($g\geq 1$), we eliminate all sites and bonds contained inside $3^{g-1}$ inverted triangles of length $L/2^g$, introducing an additional $3^{g-1}$ inner boundaries into the lattice. The procedure continues until $g=g_c$, with $L/2^{g_c} = 2a$ ($g_c = p-1$), at which stage a generation $g_c$ SG is produced: the ratio of bulk sites $n_B(g)$ to edge sites $n_E(g)$ vanishes identically when $g = g_c$ (see Appendix \ref{app:decimation} for details). This process is illustrated in Fig.~\ref{fig:decimation}.

\begin{figure}[b]
    \centering
    \includegraphics[width=0.5\textwidth]{decimation.jpg}
    \caption{Decimating a triangular lattice recursively to generate the SG. Blue (black) dots denote sites with boundary (bulk) coordination. Sites and bonds inside the red (green) triangle(s) are eliminated at the first (second) step.}
    \label{fig:decimation}
\end{figure}

Starting in the topological phase, where a chiral Majorana mode propagates clockwise along the outermost boundary, each subsequent iteration introduces additional physical boundaries into the lattice, each hosting a chiral Majorana mode propagating counter-clockwise~\cite{aliceareview}. In the thermodynamic limit $L/a \to \infty$, the decimation is repeated infinitely many times ($g_c \to \infty$) until only boundary sites are left and a chiral Majorana mode propagates along each of the infinitely many edges, resulting in a gapless spectrum. Thus, the decimation picture shows that the chiral eigenstates of the gapless topological phase are intimately linked to the Majorana edge modes of the underlying $p$-wave state. Further, the absence of \textit{any} bulk sites explains why all bulk features of the underlying model are washed out as $g\to\infty$, with the novel gapless state effectively described by a self-similar network of chiral 1D Majorana modes.

\begin{figure*}[t]
\centering
\begin{subfigure}[t]{0.45\textwidth}
\includegraphics[width=\textwidth]{cornermodespectrum.jpg}
\label{fig:hotiSpectrum}
\caption{}
\end{subfigure}\quad
\begin{subfigure}[t]{0.45\textwidth}
\includegraphics[width=0.8\textwidth]{cornerModes.jpg}
\caption{}
\label{fig:hotiCornerMode}
\end{subfigure}
\caption{(a) Spectrum and (b) non-propagating corner modes in an HOTI defined on the SC ($g=3, \gamma=0.5, \lambda=1$).}
\label{fig:corner}
\end{figure*}

Starting instead in the trivial phase, each iteration only introduces additional gaps as no edge modes appear. The self-similar arrangement of the gaps is a consequence of discrete scale invariance of the generated SG, and the trivial $\to$ topological transition on the SG can be understood as the proliferation of chiral Majorana modes which occurs during the transition on the underlying periodic lattice. Our analytic picture naturally accounts for the phase boundaries of Eq.~\eqref{bdgHam} on the SG matching those on the triangular lattice. We also expect that the gapless topological phase inherits the robustness of the edge modes against arbitrary local perturbations respecting the symmetry protecting the parent ($p+ip$) state.

To further test our decimation picture, we place the $p+ip$ Hamiltonian on the SC. This lattice can be constructed by recursively decimating a square lattice, on which a topological phase exists for $-t<\mu< t$. However, the ratio $ \lim_{g\to \infty} n_B(g)/n_E(g) \sim 5$ for the SC, resulting in more bulk than edge coordinated sites. Crucially, the distance between gapless edge modes appearing along inner boundaries at each step of the decimation process \textit{decreases} with each iteration, such that each Majorana edge mode on the recursively generated SC is separated from one with opposite chirality by $3a$. In the thermodynamic limit, these edges states back-scatter and hybridize, leading to a gapped spectrum; we thus expect that bulk features of the underlying state persist on the SC as $g \to \infty$ even in the topological phase. Results obtained by numerically diagonalizing the BdG Hamiltonian Eq.~\eqref{bdgHam} on the SC vindicate our prediction: we find a trivial ($\mathcal{C} = 0)$ and a \textit{gapped} topological (with quantized $\mathcal{C} = 1$) phase, with phase boundaries matching those of the model on the square lattice (see Appendix \ref{app:carpet} for details).

We posit that the above analysis readily generalizes to \textit{any} parent 2D topological state protected by internal symmetries: for parameters corresponding to the topological phase on a triangular lattice, the model will admit a gapless topological phase on the SG, whose physics is governed by that of the 1D gapless edge states of the parent state. On the SC, for parameters corresponding to the topological phase on the square lattice, the spectrum will remain gapped and exhibit a nontrivial quantized topological invariant. Thus, the nature of topological states on a given fractal lattice depends crucially on whether $\lim_{g\to\infty}n_B(g)/n_E(g)$ remains finite or vanishes, resulting in a gapped or gapless topological phase respectively. The results of Ref.~[\onlinecite{fractalized}], which studied the half-BHZ model~\cite{bernevigbook} on the SG and SC, are in excellent agreement with our conjecture and support the generality of our arguments.

\section{HOTI on the Sierpinski carpet}

Extending the above ideas to topological states protected by spatial symmetries requires more care, since we must ensure that no symmetries protecting the underlying state are broken at any step of the recursive decimation, in order to stay within the same phase. For instance, for a cSPT protected by $C_4$ rotation, we can start from a square lattice and recursively generate the SC through decimation, resulting in a gapped topological phase in the thermodynamic limit. To demonstrate the applicability of our general framework to this case, we have studied the paradigmatic four-band model of an HOTI, introduced in Ref.~[\onlinecite{benalcazar2017}], on the SC. The real space Hamiltonian on the square lattice is given by:
\begin{equation}
\label{cornerModeHamiltonian}
H = - \sum_{m,n} \left[ \lambda^{(1)}_{m,n} \hat{c}_{m+1,n}^\dagger \hat{c}_{m,n} +  \lambda^{(2)}_{m,n} \hat{c}_{m,n+1}^\dagger \hat{c}_{m,n} + \textrm{h.c.} \right] \nonumber,
\end{equation}
where $\hat{c}_{m,n}^\dagger, \hat{c}_{m,n}$ are fermionic creation/annihilation operators for site $(m,n)$ of the square lattice, and where
\begin{equation*}
2 \lambda^{(1);(2)}_{m,n} = \lambda (1 + (-1)^{m;n})  + \gamma (1 - (-1)^{m;n}) \, .
\end{equation*}
This model preserves $C_{4}$ rotation, time-reversal, and charge-conjugation symmetries, and presents localized corner modes when $|\gamma/\lambda|<1$. Starting from the topological phase on the square lattice, we recursively decimate the lattice to generate the SC. Each iteration creates additional inner boundaries, each hosting protected gapless corner modes since no symmetries are broken at any stage. Following our general arguments, we expect a gapped topological phase on the SC as $g\to\infty$, with modes localized along the corners of infinitely many inner edges. As shown in Fig.~\ref{fig:corner}, we indeed find a gapped spectrum and corner modes on \textit{all} inner boundaries. While we are numerically limited to $g=3$, we expect this behavior persists for larger generations. We also note that the topological nature of the SC HOTI is protected only in the presence of a particle-hole symmetry in addition to a $C_4$ symmetry: in the absence of particle-hole symmetry, the zero energy modes can be shifted around without breaking the $C_4$ symmetry~\cite{benalcazar2017}. Besides the generalization to spatial symmetries, this analysis indicates that HOTIs remain well-defined on fractal lattices as long as symmetries protecting the parent state remain unbroken.

\section{Conclusions}
\label{sec:conclusions}

In this paper, we have presented general principles which determine the fate of 2D topological states on some fractal lattices, with numerics supporting our analytic arguments. Our results strongly suggest that lattices such as the SG (SC) can support gapless (gapped) topological phases, whose properties derive from those of an underlying parent state. Understanding the role of interactions remains an important open question, as does extending these ideas to 3D topological phases on \textit{e.g.}, the Sierpinski prism, where novel behavior could result from the rich structure of surface states. A more thorough investigation of the gapless topological phase of the $p+ip$ superconductor on the SG is also warranted and could shed light on its low-energy effective field theory as well as the observed \textit{topological} metal-to-insulator transition. Finally, given the progress in fabricating fractal lattices~\cite{herek,shang,tait2015,zhang2016,kempkes2018} and in realizing HOTIs on a variety of platforms~\cite{imhof2018,garcia2018,peterson2018,schindler2018}, experimentally realizing corner modes on a fractal lattice could be within reach.

\begin{acknowledgments}
We are especially grateful to Sheng-Jie Huang and Rahul Nandkishore for discussions which inspired this work. We also acknowledge stimulating conversations and correspondence with Yang-Zhi Chou, Victor Gurarie, Michael Hermele, Sergej Moroz, Titus Neupert, and Michael Pretko. The work of SP is supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (BES) under Award number DE-SC0014415. AP is supported by a PCTS Fellowship at Princeton University.
\end{acknowledgments}


\section{Introduction} \label{sec:intro}

In essence, the matched filter (MF) is a linear low-pass filter that maximizes the signal to noise ratio of the detected signal. When the noise can be assumed as the realization of an isotropic Gaussian random field (GRF), this filter provides the greatest probability of detection for a fixed probability of false detection or false alarm (PFA) \citep{kay98}. This property makes the MF a very popular detection technique. However, MF requires the a
priori knowledge of the position of the signal within the GRF (e.g. an emission line in a spectrum or a point source in an
astronomical map). In most practical applications this is not the case. For this reason, the MF is used assuming that, if present, the position of a
signal corresponds to a peak in the matched filtered data. A detection can be assigned to the peaks exceeding a prefixed threshold. In other words, the detection procedure is not based on the entire area of the GRF but only on the subset of
the points (pixels in the case of discrete data) corresponding to the position of a peak. Recently \citet{vio16} and \citet{vio17} (hereafter, VA16 and VVA17) have shown that, when based on the standard but wrong assumption that the probability density function (PDF) of the amplitudes of the peaks of a GRF is a Gaussian, this approach may lead to a severe underestimation of the PFA.  Moreover, the PFA does not provide
the probability that a given detection is spurious but  the probability that a generic peak can exceed, by chance, a fixed threshold. The critical  point is that the threshold
has to depend on the number of peaks. The greater the number of peaks the higher the threshold. In statistics this is a well-known problem called the multiple comparisons, multiple testing or multiple hypotheses problem,
whereas in other fields it is known as the look-elsewhere effect. The majority of the
proposed solutions are essentially of a non-parametric type, such as the methods that control the family-wise error rate \citep{leh05} and the procedures to control  the false discovery rate \citep{ben95}, with the latter which has been proposed for the
astronomical applications \citep{mil01, hop02}. The main limitation of these methods is that, unlike the parametric approaches, they are not able to exploit all the available information. Parametric methods have been produced, but in general they are not
easy to use. An example is the procedure proposed by \citet{vit11}, based on the computation of the Euler
characteristic of a GRF (a quantity more difficult to compute than the peak amplitudes), which needs numerical simulations to fix the value of some fundamental parameters.
For this reason, VA16 and VVA17 have introduced an efficient parametric approach able to provide a reliable estimate,
called specific probability of false alarm (SPFA),
of the probability of a false detection. Such a quantity is computed on the basis of the PDF of the peak amplitudes of a smooth and isotropic  GRF\footnote{Actually, there are situations where the isotropy condition can be relaxed (see the last paragraph of Sect.~4 in Vio17).}. In a recent work, \citet{pav18} reach the same conclusions as VA16 and VVA17 but adopting a different approach based on the Gumbel distribution, which represents the asymptotic PDF of the extreme
(i.e. the greatest value) of an isotropic GRF.

In this paper we show that, although the two approaches produce almost the same results, the method based on the PDF of the peak amplitudes (PAM) is more flexible than the method  based on the Gumbel distribution (GDM)  and moreover it permits an easy check of the conditions of its applicability. Hence, the detection procedure is more effective when PAM is used.

In Sect.~\ref{sec:MF1} the main characteristics of MF are reviewed. The reason why it underestimates the PFA when used in the standard way is explained in Sect.~\ref{sec:PAM}. In the same section,
the method suggested by VA16 and VVA17
to correctly compute this quantity is illustrated  and the quantity SPFA introduced. The GDM approach is presented in Sect.~\ref{sec:GDM}.
Finally,  in Sects.~\ref{sec:comparison}-\ref{sec:ALMA} the PAM and GDM performances are tested on a set of simulated GRFs as well on an interferometric map obtained with the Atacama Large Millimeter/submillimeter Array (ALMA) and the final remarks are deferred to Sect.~\ref{sec:conclusions}.

\section{Matched filter} \label{sec:MF1}

Given a discrete observed signal $\boldsymbol{x}$, the model assumed in the MF approach is
$\boldsymbol{x} = \boldsymbol{s} + \boldsymbol{n}$, where $\boldsymbol{s}$ is the deterministic signal to detect and $\boldsymbol{n}$ a zero-mean GRF  with known covariance matrix,
\begin{equation} \label{eq:C}
\boldsymbol{C} = {\rm E}[\boldsymbol{n} \boldsymbol{n}^T].
\end{equation}
Here, symbols  ${\rm E}[.]$ and  $^T$ denote the expectation operator and the vector or matrix transpose, respectively.

Under these conditions,  according to the Neyman-Pearson theorem \citep{ney33, kay98}, a detection is claimed when
\begin{equation} \label{eq:test1}
\mathcal{T}(\boldsymbol{x}) = \boldsymbol{x}^T \boldsymbol{f} > \gamma,
\end{equation}
with $\gamma$ a real constant and
\begin{equation} \label{eq:mf}
\boldsymbol{f}_s = \boldsymbol{C}^{-1} \boldsymbol{s}.
\end{equation}
Here $\boldsymbol{f}_s$ represents the matched filter.
The main characteristic of the MF is that it maximizes the probability of detection under the constraint of a fixed PFA .

Matched filter  works properly under two assumptions. The first assumption is that the signal $\boldsymbol{s}$ is known. Actually, in practical applications only the template $\boldsymbol{g}$ of the signal $\boldsymbol{s} = a \boldsymbol{g}$ is available
but not its amplitude $a$. However, this does not represent a true problem since the  MF in the form
\begin{equation} \label{eq:mf2}
\boldsymbol{f}_g = \boldsymbol{C}^{-1} \boldsymbol{g}
\end{equation}
does not affect the PFA but only the probability of detection (VA16). The second assumption requires that $\boldsymbol{x}$ and $\boldsymbol{s}$ have the same size. This implicitly means that the position of $\boldsymbol{s}$ within
$\boldsymbol{x}$  is known. Problems come up when this last assumption is relaxed and the sizes of $\boldsymbol{s}$ are smaller than the sizes of $\boldsymbol{x}$ (e.g. an emission line in a spectrum).
If the amplitude $a$ is also unknown, the standard approach consists in cross-correlating $\boldsymbol{x}$ with the MF given by~Eq. \eqref{eq:mf2} obtaining signal $\boldsymbol{\mathfrak{x}}$ typically standardized to zero-mean and unit variance.
A detection is claimed when a peak in $\boldsymbol{\mathfrak{x}}$ exceeds a threshold set to $u$.
If the number of signals $\boldsymbol{s}$ present in $\boldsymbol{x}$ is also unknown, this procedure has to be applied to all the most significant peaks. It is a widespread practice that the corresponding PFA is given by
\begin{equation} \label{eq:fd1}
\alpha = \Phi_c(u),
\end{equation}
where  $\Phi_c(u) = 1 - \Phi(u),$ with $\Phi(u)$ the standard Gaussian cumulative distribution function (CDF).

\section{Detection procedure based on the PDF of the amplitude of the peaks} \label{sec:PAM}

In VA16 and VVA17  it has been shown that the computation of the PFA by means of Eq.~\eqref{eq:fd1} can lead to a severe underestimation of the latter quantity. This is because the PDF of the peaks of a GRF is not a Gaussian as
implicitly assumed in Eq.~\eqref{eq:fd1}. For this reason, the correct PFA has to be estimated by means of
\begin{equation} \label{eq:corra}
 \alpha = \Psi_c(u),
\end{equation}
where
\begin{equation}
\Psi_c(u)= 1 - \Psi(u),
\end{equation}
with
\begin{equation} \label{eq:Psi}
\Psi(u)=\int_{-\infty}^{u} \psi(z) dz,
\end{equation}
and
\begin{equation} \label{eq:pdf_z1}
\psi(z) = \frac{\sqrt{3 - \kappa^2}}{\sqrt{6 \pi}} {\rm e}^{-\frac{3 z^2}{2(3 - \kappa^2)}} + \frac{2 \kappa z \sqrt{\pi}}{\sqrt{6}} \phi(z) \Phi\left(\frac{\kappa z}{\sqrt{3 - \kappa^2}} \right)
\end{equation}
 for the 1D case, and
\begin{multline} \label{eq:pdf_z2}
\psi(z) = \sqrt{3} \kappa^2 (z^2-1) \phi(z) \Phi \left( \frac{\kappa z}{\sqrt{2 - \kappa^2}} \right) + \frac{\kappa z \sqrt{3 ( 2 - \kappa^2)}}{2 \pi} {\rm e}^{-\frac{z^2}{2 - \kappa^2}}\\
+\frac{\sqrt{6}}{\sqrt{\pi (3 - \kappa^2)}} {\rm{e}^{-\frac{3 z^2}{2 (3-\kappa^2)}}} \Phi\left( \frac{\kappa z}{\sqrt{(3 - \kappa^2) (2 - \kappa^2)}} \right)
\end{multline}
for the 2D case.
These expressions represent  the PDF of the local maxima of a zero-mean, unit-variance homogeneous GRF \citep{che15a, che15b} \footnote{For the three-dimensional case see  \citet{che15b}.}.

Here, $\kappa$ is a parameter given by
\begin{equation} \label{eq:kd}
\kappa = - \frac{\varrho'(0)}{\sqrt{\varrho''(0)}},
\end{equation}
where $\varrho'(0)$ and $\varrho''(0)$ are, respectively, the first and second derivative with respect to $r^2$ of the autocorrelation function $\varrho(r)$ of the GRF at $r=0$, with $r$ the inter-point distance of the random field.
As shown in VVA17, this parameter can be estimated by means of a maximum likelihood approach,
\begin{equation} \label{eq:ml}
\hat{\kappa} = \underset{\kappa }{\arg\max} \sum_{i=1}^{N_p} \log{\left(\psi(z_i; \kappa)\right)},
\end{equation}
where $\{ z_i \}$, $i=1,2,\ldots, N_p$, are the local maxima in $\boldsymbol{\mathfrak{x}}$  \footnote{We recall that the function ``$ \underset{x}{\arg\max}[ H(x)]$'' provides the value of $x$ for which the function $H(x)$ has the greatest value.}.
In this way, it is possible to avoid the use of Eq.~\eqref{eq:kd}, which requires knowledge of $\varrho(r)$.
The condition of validity of these expressions is that $\varrho(r)$ is differentiable at least six times with respect to $r$, but it is conjectured that four times should be sufficient (Cheng: private communication).
Especially in astronomy, this is a weak condition because due to the point spread function of the instruments, the MF often has a Gaussian-like shape. Hence, the resulting $\boldsymbol{\mathfrak{x}}$ are characterized by very smooth $\varrho(r)$.

It is necessary to stress that the PFA given by Eq.~\eqref{eq:corra} does not provide the probability $\alpha$ that a specific detection is spurious, but the probability that a generic peak due to the noise in $\boldsymbol{\mathfrak{x}}$
can exceed, by chance, the threshold $u$.
If the number of such peaks is $N_p$, then a number  $\alpha \times N_p$  among them  is expected to exceed the prefixed detection threshold.
As a consequence, in spite of a low PFA, the reliability of a detection could actually be small. The solution proposed in VVA17 consists of a
preselection based on the PFA and then in the computation of the specific probability of false alarm (SPFA) for each detection.
This quantity can be computed by means of the order statistics, in particular by exploiting  the statistical characteristics of the greatest value of a finite sample of identical and independently distributed (iid) random variable from a given
PDF  \citep{hog13}. Under the iid condition, the PDF $\upsilon(z_{\max})$ of the largest value among a set of $N_p$ peaks $\{ z_i \}$ is given by
\begin{equation} \label{eq:gz}
\upsilon(z_{\max}) = N_p \left[ \Psi(z_{\max}) \right]^{N_p-1} \psi(z_{\max}).
\end{equation}
Hence, the SPFA can be evaluated by means of
\begin{equation} \label{eq:intz}
\alpha = \int_{z_{\max}}^{\infty} \upsilon(z') dz'.
\end{equation}
The numerical evaluation of this integral does not present particular difficulties since
\begin{align}
\alpha & = N_p \int_{z_{\max}}^{\infty} \left[\Psi(z) \right]^{N_p-1} d\Psi(z); \\
          & = \left[ \Psi(z) \right]^{N_p} \Big|_{z_{\rm max}}^{\infty}; \\
          & = 1-\Upsilon(z_{\rm max}),
\end{align}
with
\begin{equation} \label{eq:Ups}
\Upsilon(z_{\rm max})=\left[ \Psi(z_{\rm max}) \right]^{N_p}.
\end{equation}
This procedure is implicitly based on the assumption that there is only a signal $\boldsymbol{s}$ in $\boldsymbol{x}$. If the actual number is unknown, it has to be cyclically applied to all the remaining most prominent peaks in order of decreasing amplitude and stopped
when the estimated PFA is greater than a prefixed $\alpha^*$. In principle, $N_p$ should be lowered by one unit after any detection. This is because $N_p$ represents the number of peaks due to the noise.
However, since in the practical applications $N_p$ is on the order of thousands if not tens of thousands, this step is irrelevant.

As a final note, it is necessary to stress that the peaks  of a isotropic GRF, with $\varrho(r)$ typical of many astronomical observations, usually have a spatial distribution different from the spatial pattern characteristic
of a complete spatial random point process (CSRPP).
This is visible in the top panels of  Fig.~\ref{fig:fig_distance} where the PDF of the nearest neighbor
distances of the peaks of a simulated GRF with $\varrho(r)$ given by a two-dimensional circular Gaussian (typical of many astronomical images) with dispersion $\sigma_G=3$ in pixel units is compared with the corresponding
PDF of a simulated CSRPP.  From this figure it is clear that the PDF related to the peaks lacks
small values. Hence, the two processes are different on short spatial scales. This means that the iid condition for the peak amplitudes is not necessarily valid.
However,  on greater spatial scales, when $\varrho(r)$ is narrower than the
area spanned by the data (a basic situation for the application of the MF), this condition can be expected to hold with good accuracy. The rationale is that two generic points of a GRF with a distance $r$ such that $\varrho(r) \approx 0$
are essentially independent. The same holds for two generic peaks. This is confirmed by the bottom-left panel of Fig.~\ref{fig:fig_distance}  where the sample pair correlation function\footnote{The pair correlation function $\rho(r)$
of the spatial distribution of a set of points is given by
$\rho(r) = K'(r)/2 \pi r$ with $K'(r)$ the derivative of the Ripley's $K$-function with respect to $r$. For a CSRPP it is $\rho(r) = 1$ independently of $r$ \citep{bad16}.}  $\rho(r)$
of the peaks indicates that, for $r \ge 7,$ their spatial distribution is compatible with a  CSRPP.
For comparison, the bottom-right panel displays $\rho(r)$ for a true CSRPP.
The different behavior on small scales makes the peaks of the GRF show a more uniform spatial distribution than the spatial distribution of the points of the CSRPP (see Fig.~\ref{fig:fig_maps}).
In conclusion,  since  $\varrho(7) \approx 0$, it means that, for $r \ge 7,$ the peak amplitudes can be considered iid. As a result, most of the peak amplitudes of a GRF with $\varrho(r)$ typical of many astronomical applications (essentially Gaussian-like) can be expected to be approximately  iid. Hence,  Eq.~\eqref{eq:gz} is still applicable but possibly with an
effective number smaller than $N_p$ \citep[see Sect. 6 in][]{may05}.
This last point is due to the dependence among a set of random variables, which lowers their number of degrees of freedom\footnote{The term degrees of freedom refers to the number of
items that can be freely varied in calculating a statistic without violating any constraints.}.

\section{Detection procedure based on the Gumbel distribution} \label{sec:GDM}

Independently of their PDF, the CDF of the greatest value $x_{\rm max}$ of a set of $N \longrightarrow \infty$ iid random variables $\boldsymbol{x}$ is given by
\begin{equation} \label{eq:gumb}
G(x_{\rm max})=\exp{\left[ -(1+\gamma_g y)^{1/\gamma_g} \right]}
,\end{equation}
where
\begin{equation}
y=\frac{x_{\rm max}-a}{b}
\end{equation}
with $a$ and $b$ the location and the scale parameter, respectively. Three cases are possible, according to  $\gamma_g < 0$, $\gamma_g=0$ and $\gamma_g > 0$ \citep{cas04}. Assuming that $x_{\rm max}$ corresponds to a peak,
that is, $x_{\rm max}=z_{\rm max}$, and in the case of an isotropic GRF, it has been shown \citep{col11}
that the CDF of $z_{\rm max}$ is given by Eq.~\eqref{eq:gumb} with $\gamma_g =0$,
\begin{equation} \label{eq:gumb0}
G(z_{\rm max})=\exp{\left[ - {\rm e}^{-y} \right]},
\end{equation}
which represents the Gumbel distribution. With a two-dimensional, zero-mean, unit-variance GRF and in the regime  of slight clustering of the peaks, Eq.~\eqref{eq:gumb0} can be approximated by
\begin{equation}
G(z_{\rm max}) \approx \exp{\left[- \frac{\gamma^2 A z_{\rm max} }{(2 \pi)^{3/2} R_*^2} {\rm e}^{-z_{\rm max}^2/2} \right]},
\end{equation}
where $A$ is the area of the GRF,
\begin{equation}
\gamma=\frac{\sigma_1^2}{\sigma_0 \sigma_2},
\end{equation}
and
\begin{equation}
R_*=\frac{\sigma_1}{\sigma_2} \sqrt{2}.
\end{equation}
The SPFA for a specific $z_{\rm max}$ is given by $\alpha=1-G(z_{\rm max})$.

If the GRF is obtained by filtering a noise process characterized by a scale-free power spectrum $P(k)$,
\begin{equation}
P(k) = B k^n,
\end{equation}
with a filter $w_l(k)$ ($l$ a scale parameter), it is
\begin{equation} \label{eq:sigma_i}
\sigma_j^2= \frac{1}{2 \pi} \int_0^{\infty} k^{2 j +1} P(k) w_l^2 (k) dk.
\end{equation}
Details of the three-dimensional case can be found in \citet{col11}.

\section{Comparison with simulated data} \label{sec:comparison}

Given the different derivation of PAM and GDM, it is necessary to compare their performances. With this aim, a numerical experiment was carried out that is based on $5 x 10^3$
numerical simulations of $500 \times 500$ pixels GRFs characterized by a  circular Gaussian autocorrelation function with dispersion $\sigma_G=3$.
This kind of GRF corresponds to what is obtained by filtering a two-dimensional white-noise with a circular Gaussian $w_l(x)$ with standard deviation $l=3/\sqrt{2}$. The resulting $G(z_{\rm max}) $ can be expressed in the analytical form
\begin{equation} \label{eq:GG}
G(z_{\rm max}) \approx \exp{\left[- \frac{r}{4 \sqrt{2 \pi}}  z_{\rm max} {\rm e}^{-z^2_{\rm max}/2}  \right]},
\end{equation}
where $N_*$ is the ratio of the total area of the GRF to the window area $\pi  l^2$ \citep{pav18}.

Concerning the PAM, $\kappa=1$ and $N_p=2552$ are the theoretical values obtained as explained in VVA17, whereas
the corresponding mean value of the maximum likelihood parameter $\hat{\kappa}$ and the number of peaks $N_p$ obtained from each of the simulated GRFs  are  $\bar{\kappa}=1.01 \pm 0.01$ and $\bar{N}_p = 2435 \pm 26$.
The left panel of Fig.~\ref{fig:fig_pdf_circ} shows the PDFs $\upsilon_{\rm Th}(z_{\rm max})$ and $\upsilon_{\rm Mean}(z_{\rm max})$
corresponding to the theoretical value $N_p$ and the mean
value $\bar{N}_p$, respectively. Both slightly differ from the histogram $H(z_{\rm max})$.
Given the small statistical fluctuation of $\bar{\kappa}$ and $\bar{N}_p$, such a discrepancy can be explained by the condition of complete spatial randomness assumed for the spatial distribution of the peaks.
In particular, the number of effectively iid peaks has to be expected to be smaller than both $N$ and $\bar{N}_p$. Indeed, a maximum likelihood fit $\upsilon_{\rm ML}(z_{\rm max})$  of the PDF given by Eq.~\eqref{eq:gz} to the set of simulated
$\{ z_{\rm max} \}$ with $N_p$ as free parameter provides $\tilde{N}_p=2210$.
The right panel of Fig.~\ref{fig:fig_pdf_circ}  shows that something similar holds for the GDM since the theoretical value $N_*=17684$ is greater than the number $\tilde{N}_*=15330$ corresponding to the maximum likelihood fit
$g_{\rm ML}(z_{\rm max})$.
Here, it is worth stressing that $\upsilon_{\rm Th}(z_{\rm max})$ and $g_{\rm Th}(z_{\rm max}),$ and hence the corresponding CDFs $\Upsilon_{\rm Th}(z_{\rm max})$ and $G_{\rm Th}(z_{\rm max})$ shown in
Fig.~\ref{fig:fig_cdf_circ} are almost identical. The same holds for
$\upsilon_{\rm ML}(z_{\rm max})$ and $g_{\rm ML}(z_{\rm max})$ and the corresponding CDFs $\Upsilon_{\rm ML}(z_{\rm max})$ and $G_{\rm ML}(z_{\rm max})$ again shown in Fig.~\ref{fig:fig_cdf_circ}.
As expected the two approaches are practically equivalent, given that both methods are essentially based on the same assumptions. Moreover, in Fig.~\ref{fig:fig_cdf_circ} all these CDFs appear close to the sample CDF. This indicates
that the value of the parameters $N_p$ and $N_*$ is not of critical importance.

A further comparison of the two methods concerns the computational burden. In this respect the GDM is superior
in situations where the autocorrelation function has circular symmetry (i.e., it is characterized by only a scale parameter), since the more expensive step consists in the computation of the quantities $\sigma_j$, $j=0,1,2$, which requires the numerical computation of three one-dimensional integrals. On the other hand, PAM requires the numerical computation of a number of integrals equal to the number of peaks selected for the detection. In general  the difference in computational
time is on the order of a few seconds even in the case of large images. The computational superiority of GDM vanishes in the presence of autocorrelation functions when they are not circularly symmetric. This is because the three one-dimensional integrals become
two-dimensional for a two-dimensional GRF and three-dimensional  for a three-dimensional GRF, which is computationally expensive. Finally, as shown in Eq.~\eqref{eq:sigma_i}, the GDM requires knowledge of the Fourier transform $w(k)$ of the filter used to
obtain the GRF. If this is not available, it must be evaluated numerically with an additional computational cost.
In addition, the PAM is fully automatic given that the values of the parameters $\kappa$ and $N_p$ are estimated only from the peak amplitudes. Hence, contrary to the
GDM, there is no necessity to write a  specific code for a given autorrelation function.

A final aspect to consider is the easiness with PAM to check if the PDF of the peak amplitudes corresponding to $\hat{\kappa}$ is compatible with the sample PDF.  This check is fundamental for the reliability of the detection. With GDM only the Gaussianity of the entries of the random field can be checked, hence the algorithm has to be used as a black box. In this respect, Fig.~\ref{fig:fig_pdf_circ_01_01} (to compare to Fig.~\ref{fig:fig_pdf_circ}) shows the result of a numerical experiment similar
to that presented above, but
where the dispersion of the circular Gaussian autocorrelation function is $\sigma_G=1$. It is clear that both PAM and GDM work really badly. The reason is that, again, both methods are developed in the context of a continuous GRF. With
$\sigma_G=1$ the effects of the discretization become important. This point can be observed in Fig.~\ref{fig:fig_check}, which compares $\psi(z)$ for a typical realization of the GRF  when $\sigma_G=3$ (left panel) to that when
$\sigma_G=1$ (right panel). As a one-sample Kolmogorov-Smirnvov test with a $99\%$ confidence level indicates, only when $\sigma_G=3$ is the sample distribution of the peak amplitudes compatible with the PDF expected for a GRF and hence only in this case can the two methods be safely used. The conclusion is that a blind application of both methods can lead to wrong results.

\section{Comparison with an ALMA map} \label{sec:ALMA}

The peak amplitudes method and GDM are here applied to the standardized zero-mean unit-variance interferometric ALMA map in Fig.~\ref{fig:fig_map} with the aim of detecting point sources.  The data are taken from the ALMA project 2012.1.00173.S, a mosaicing of the Hubble Ultra-Deep Field (HUDF) in continuum \citep{dun17}. The HUDF was observed using a 45-pointing mosaic, with each pointing separated by 0.8 times the antenna beamsize. The details of the observing and data reduction procedures can be found in Dunlop et al. (2017). For the purpose of this work, in order to simplify the analysis, we have cut the outer edges of the total image and produced a symmetric, square-shaped image of size $1075 \times 1075$ pixels
(corresponding to $110\times 110$ arcsec$^2$). On this inner part we applied our algorithm and found $11959$ peaks.

As explained in Vio16, in ALMA maps the point sources and the blob-shaped structures due to the noise have a similar aspect and the MF cannot be applied. Consequently,
the detection test becomes a thresholding test, where a peak in the map is assigned to a point source if it exceeds a given threshold.
Before applying the detection algorithm, it is necessary to check whether the map meets the requirements mentioned above. In particular, the Gaussianity and isotropy of the noise background, the compatibility of the sample  PDF of the peak amplitudes
with the theoretical PDF $\psi(z),$  and the fact that the spatial distribution of the peaks is compatible with a CSRPP.

The results of these checks are presented in  Fig.~\ref{fig:fig_test}. In particular, the top-left panel shows that, although similar to a Gaussian, the histogram of the entries of the map indicates a PDF slightly leptokurtik\footnote{A leptokurtik PDF has a shape more peaked than a Gaussian.}.
The excess of values close to the mean can be understood taking into account that in general a zero-mean, isotropic GFR with a smooth autocorrelation function shows a lack of entries with high absolute value\footnote{This can be understood taking into account that an isotropic GFR with a specific autocorrelation function can be obtained by filtering a Gaussian white-noise process by means of a linear filter.
In most of the astronomical applications, the filters necessary to obtain a specific autocorrelation function are of low-pass type (e.g., for a Gaussian autocorrelation function the filter is also a Gaussian). As a consequence, in the resulting GRF the greatest values of the white-noise process are smoothed out in favor of values closer to the mean. As a matter of fact, this effect becomes important only for the maps
with dimensions comparable to the extension of the area where the autocorrelation function is significantly different from zero, but this is not the case for the map under consideration.}.
More importantly, from the top-right panel of the same
figure the histogram of the peak amplitudes appears compatible with the $\psi(z)$ corresponding  to a maximum likelihood estimate $\hat{\kappa}=0.98$. This is confirmed by a Kolmogorov-Smirnov test at a confidence level of $99\%$,

The isotropy of the noise background is supported by the bottom-left panel of Fig.~\ref{fig:fig_test}, which displays the sample autocorrelation functions along the $X$ and $Y$
directions versus the corresponding slices of a two-dimensional circular Gaussian resulting from a least-squares fit. The good agreement indicates that the ALMA map is compatible with an isotropic Gaussian GRF
whose autocorrelation function is a two-dimensional circular Gaussian with dispersion $\sigma_G \approx 3$.

Finally, as is visible in the bottom-right panel of the same figure, the sample pair correlation function $\rho(r)$ is compatible with a CSRPP for $r \ge 7$.
Since $\varrho(7) \approx 0,$ this means that, as seen in Sect.~\ref{sec:comparison}, the peak amplitudes do not present a relevant mutual statistical dependence, hence most of them can be considered iid.

Both PAM and GDM produce six statistically significant  detections, which are labeled in Fig.~\ref{fig:fig_map} with a number in order of decreasing intensity and shown in more detail in Fig.~\ref{fig:fig_submaps}. The SPFAs coming from the two methods are very similar, since for PAM the values are $0$, $4.65\rm{E}-11$, $1.51\rm{E}-06$, $5.76\rm{E}-06$, $3.68\rm{E}-02$, $7.05\rm{E}-02$, whereas for the GDM they become $0$, $0$, $1.54\rm{E}-06$, $5.84\rm{E}-06$, $3.74\rm{E}-02$, $7.16\rm{E}-02$.
This result confirms the equivalence of the two methods as it concerns the detection performance. Table~\ref{tab:sources} reports the coordinates of the detected sources and their identification with the sources reported in \citet{dun17}.

\begin{table*}
\caption{Identification of the sources on the map in Fig.~\ref{fig:fig_map}.}
\label{tab:sources}
\hskip -0.7cm
\begin{tabular}{c c c c}
\hline\hline
\\
Source number  & identification & coordinates & \\
& & RA & DEC\\

\hline
\\

Source 1 &   UDF2 & 3h32m43.53s &  -27d46$^\prime$39.25$^{\prime\prime}$ \\
Source 2 &   UDF3 & 3h32m38.55s &  -27d46$^\prime$34.57$^{\prime\prime}$ \\
Source 3 &   UDF5 & 3h32m36.96s &  -27d47$^\prime$27.13$^{\prime\prime}$ \\
Source 4 &   UDF4 & 3h32m41.02s &  -27d46$^\prime$31.58$^{\prime\prime}$ \\
Source 5 &   UDF6 & 3h32m34.43s &  -27d46$^\prime$59.72$^{\prime\prime}$ \\
Source 6 &  NOT IDENTIFIED &3h32m38.66s &  -27d48$^\prime$06.12$^{\prime\prime}$ \\

\\
\hline

\end{tabular}
\end{table*}

\section{Conclusions}  \label{sec:conclusions}

In this paper the performances of two techniques, the  PDF of the peak amplitudes method (PAM) and the Gumbel distribution method (GDM), have been compared in the context of the detection of weak signal embedded in noise.
The two methods have been applied to simulated signals and to observations taken with the ALMA interferometer.
We have shown that the two approaches are almost perfectly equivalent in their detection capability, but PAM proves to be more flexible and, at the same time, allows for an easy control of the condition of applicability of the technique. Hence, it appears
more appropriate in the search for weak sources in observations dominated by noise.

\begin{acknowledgements}

\end{acknowledgements}


\section{Motivation}
Another consequence of the quick growth of interest in cryptocurrencies and the lack of proper system security audits were the many hacks of exchanges.
The hacks dispossessed them of hundreds of thousands of bitcoins.
Because the Bitcoin network offers pseudo-anonymity to its users, it prevents any linking between addresses and any personally identifiable information\footnote{This assumes that the user is versed in privacy and does not for example reuse addresses.
It also assumes the use of techniques such as CoinJoin to prevent addresses from being connected. These techniques add additional overhead and are too advanced for the average Bitcoin user so they received little to no traction}.
Usually, when some bigger hacks happen on popular cryptocurrency exchanges, there is nothing we can do but wait until the thieves want to use their Bitcoins to buy goods or exchange them.
Frequently, users with tainted Bitcoins use mixers to launder them.
A mixer works so that the tainted Bitcoins which are deposited, are put into reserve whereas the clean Bitcoins are paid out to the customer on a new address.
Doing that, we lose the trace of the original address as it can not be directly linked to the new one. \\
This alongside exchanges not verifying their users allows malevolent actors to get by unscathed. With the increasing regulation on cryptocurrencies, financial institutions working with such currencies will be required to do extensive due diligence.
This will be mandatory not only of the involved parties but also of the origin of the coins giving every deposit into their system a risk score. This will raise the awareness about which coins exchanges process to lower their exposure in case of, for example receiving a subpoena from a legal entity investigating, to return them - like credit card charge-backs.
\\
In this paper, we propose a new method of determining the risk factor of received coins from a specific address based on their history and interactions with the network called TaintRank.
By having such information, vendors or exchanges can act differently based on the TaintRank the depositing address has.
Using the TaintRank score, an exchange can require additional verification or a store can refuse participate in the purchase.
Such a ranking system is not expected to be bulletproof, so we should not make harsh actions solely based on TaintRank score a certain address has.

\section{Related work}
The area of cryptocurrency network and transaction analysis has received quite a lot of attention in recent years. Researchers have done extensive research on the analysis of the level of anonymity that the network provides \cite{reidAnonimity} or with the use of different third-party services called mixers or tumblers \cite{de2017analysis,van2018bitcoin,moser2013inquiry}. These services promise to obfuscate the real identities of the owners by randomly exchanging the coins between different users. Such services only disconnect addresses and make direct connections impossible but in reality, they cannot launder the coins in the true essence of the word: someone will always receive marked coins and we can track them even through such services. \\
There has been some research in risk analysis of the Bitcoin network such as in Spagnuolo et al.\cite{spagnuolo2014bitiodine} where they propose a system capable of clustering different addresses, classifying them based on the activity and providing insight into how they interact with the network.
Malovrh \cite{malovrh2018bitcoin} investigates possibilities of using supervised learning as a technique for detecting different anomalies in the Bitcoin network which include criminal activities.
DebtRank, proposed by Battiston et al. \cite{battiston2012debtrank}, evaluates financial institutions based on the institutions they are loaning or debiting from.
Meiklejohn et al. \cite{meiklejohn2013fistful} deals with understanding the traceability of Bitcoin flows. It gives us a good estimation of how we expect stolen Bitcoins will move and be spent or exchanged. It also touches on some limitations of mixing services, which fail if we are dealing with multiple thousand stolen Bitcoins.

\section{Methods}\label{methods}
Methods we tested are based on node importance. We first construct a directed network where each node represents a Bitcoin address and each directed link a weighted transaction between the two addresses.
Each of the nodes in the network receives its own TaintRank based on different parameters such as node degree and edge weight.
We use prior knowledge to determine the tainted nodes. This knowledge included addresses belonging to a thief who stole Bitcoins from an exchange.
Since the entire Bitcoin network is too large to handle, we use some timeline after a specific event happened. In our case, this was a theft of Bitcoins from an exchange.
We took the data for 1 month after the attack happened and determined the TaintRank score for nodes in the network.\\
We test multiple different approaches to spread the taint throughout the network. First, we test a method where each node spreads its TaintRank by a fixed amount no matter the link weight or other parameters.
While we would lose track of the specific address that used a mixer, the person behind it will have an equal TaintRank score on their new addresses afterwards.
The negative side of this approach is, that even if an address receives a tiny fraction of Bitcoin from a very tainted address, it itself becomes equally tainted.\\
The further approaches we test take other parameters into the equation.
We use distance from the original thief node, weight of the links, value of the node that is linking to us and other parameters to determine TaintRank of the nodes.
For example, we did not give the same TaintRank score to each of the outgoing transaction, but we set its size proportional to the weight of the transaction.
Doing that, we exclude some nodes not involved in illegal activities but would be tainted using the first method.
But with this, less strict method, we are also open to attacks, that could reduce the TaintRank score of criminal addresses.
\subsection{Hypothesis}
The methods we test differ whether they give out TaintRank proportional to the weight of the link or whether they always give out the same amount no matter the weight.
In the first scenario, we are expecting to be vulnerable to a kind of attack, that spreads Bitcoins into many small fractions, and then ultimately joins some of them together to get almost untainted Bitcoins in some new address.
The second approach would, taint each node that interacts with tainted node (receives Bitcoin from) by a fixed amount.
Because of that, a criminal who knows how the approach works might just send tiny fractions of Bitcoins to random addresses to taint them, with the same amount he is tainted.
Doing that, many regular users will receive a high TaintRank, resulting in more work for the companies trying to filter genuinely tainted addresses.
The perfect combination will most likely comprise some combination of both approaches, taking many parameters into consideration when deciding for the TaintRank score of a node.

\section{Selected data}
Since its inception, the Bitcoin network has grown in size exponentially. It has millions of addresses and hundreds of millions of transaction. This scale poses a great challenge when trying to analyze it. To make the problem more contained, we select a subset of the network which proves to be a good decision.
We extract 285,591 transactions in a period between 20.06.2011 to 20.7.2011 which amasses to a network with parameters that can be seen in table~\ref{tab:network_data}.

\begin{table}[H]
    \begin{center}
    \caption{Bitcoin network data for the period between 20.06.2011 to 20.7.2011 .}
        \begin{tabular}{ l | c | c }
        \textbf{nodes} & \textbf{links} & \textbf{$\avg{k}$}\\
         \hline
        505,473 & 1,143,444 & 4.524  \\
        \end{tabular}
    \label{tab:network_data}
    \end{center}
\end{table}
The reason we chose this time frame was that one of the bigger Bitcoin thefts happened while the network was still small. This was the Mass MyBitcoin Theft which happened between 20th and 21st of June 2011.
The theft is associated with the address $1MAazCWMydsQB5ynYXqSGQDjNQMN3HFmEu$ which we use as our starting point.
We extract the data using a Bitcoin database index into a public Google BigQuery database which allows us to search with ease for the data using the well-known SQL query language. This removes a lot of issues because we do not have to write a program which fetches transactions from a full-history Bitcoin node which are not public. Using a query in BigQuery we exported all the inputs and the outputs of transactions for that given period. This gives us 100MB JSON file which we further process into a usable graph representation in the edgelist format using a custom script written in Golang. After that, we use scripts written in Python to analyze the graphs and try applying our taint analysis methods.
While processing the data we thought of a heuristic approach which could reduce the size of the overall network. The transactions on the Bitcoin network are made up of different inputs which map to outputs of earlier transactions. To spend an output in a new transaction the owner of the output must provide his signature. Given that a transaction contains inputs from multiple different addresses and that the transaction can be signed only after it was created, we believe it is safe to assume that the one who constructed the transaction is also in control of all the private keys used to sign the transaction. The counter-example to this would be the CoinJoin \cite{maurer2017anonymous} principle. This effectively allows multiple disparate entities to include their unspent outputs into a common transaction. But it also requires a party to act as a synchronizer between all other parties. This creates a high level of complexity which leads us to believe that such principles are rarely used.
Treating addresses of inputs in a transaction as owned by a common entity allowed us to reduce the size of our selected the network by almost a half.

\subsection{Degree distribution}
Degree distribution in our network could tell us a lot about the nodes.
We can assume that nodes with an extremely high degree (hubs) are some sort of exchanges or perhaps mixers.
For example, some nodes are actually Bitcoin exchanges, and we should not spread TaintRank to them since we assume that exchanges would check the incoming addresses, and not accept funds from a (overly) tainted node.
Firstly, we plot the degree distribution for both in and out degrees, which can be seen in figure~\ref{fig:degree_dist}.
\begin{figure}[h!]
    \includegraphics[width=8cm]{degree_dist.png}
    \centering
    \caption{Degree distribution in the created network for both in and out degree.}
    \label{fig:degree_dist}
\end{figure}

As seen from the plot, both vaguely resemble a power-law distribution.
But considering the spike in the beginning of $p_k^{out}$ and some clearly visible spikes at node degree of 100 and 1000, we should not conclude that the network follows a power-law degree distribution.
For further examination, we look at five addresses with the highest node degree as seen in table~\ref{tab:top_node_deg}

\begin{table}[H]
    \begin{center}
    \caption{Top five Bitcoin addresses in our network that have the highest $p_k^{in}$ and $p_k^{out}$. In the last column there is the believed owner.}
        \begin{tabular}{ p{5cm} | p{0.7cm} | p{1.6cm} }
        \small
        \textbf{Address} & \textbf{Degree} & \textbf{Belongs to}\\

         \hline
         $p_k^{in}$\\
          \hline
          \footnotesize{1eHhgW6vquBYhwMPhQ668HPjxTtpvZGPC} & 2,283  & Mt.Gox\\
          \footnotesize{1PJnjo4n2Rt5jWTUrCRr4inK2XmFPXqFC7} & 1,776  & Slush's pool\\
          \footnotesize{1CxjYUcBFmBftskGFST9MXhyQMQUNTuBzi} & 1,239  & Unknown\\
          \footnotesize{1MAazCWMydsQB5ynYXqSGQDjNQMN3HFmEu} & 627  & MyBTC thief\\
          \footnotesize{1MbSn15MZWNbkNWF72KopyQenCV2zdcWvr} & 598  & Unknown\\
          \hline
         $p_k^{out}$\\
          \hline
          \footnotesize{1MaZAHzEFfinRJ2dwK6YtNDfvWMBkiAxDr} & 27,055  & Unknown\\
          \footnotesize{18qr2srETSvQq4kP7yBYRqQ4LzmjhtRmcD} & 26,975  & Unknown\\
          \footnotesize{1PJnjo4n2Rt5jWTUrCRr4inK2XmFPXqFC7} & 10,607  & Slush's pool\\
          \footnotesize{1AgwESN7RKNZtaqzbqu6kPg3RS6C2qCgHi} & 7,548  & Unknown\\
          \footnotesize{1AZUPm5PC5QguquNsBg7HhWUYz5dfm2nU9} & 7,527  & Unknown\\
          \hline
        \end{tabular}
    \label{tab:top_node_deg}
    \end{center}
\end{table}

We can see that the the highest out-degree is one magnitude higher than the highest in-degree, which tells us that most addresses spread the funds they are sending out to multiple destinations.
Based on the findings of these addresses, we could adjust our algorithms, to take into account what some of these nodes represent.
But the issue is that most of the addresses have an unknown owner.

\section{Results}
With the methods described in ~\ref{methods}, we test our approach on the network.
In the selected data, we know the address where the initial stolen funds went and we base our propagation on it.
Based on that, we start to propagate the taint throughout the network by recursively following all the out links.
The number of nodes which could be reached following the out links from the original address is extremely large (227,656), indicating that the funds have most likely made contact with an exchange of some sort.
In the table~\ref{tab:subgraph_data}, we can see the characteristics of this subgraph.
When comparing them to the original graph we see that the average node degree is smaller, but no fundamental differences.
\begin{table}[H]
    \begin{center}
    \caption{Subgraph network data for the period between 20.06.2011 to 20.7.2011 .}
        \begin{tabular}{ l | c | c }
        \textbf{nodes} & \textbf{links} & \textbf{$\avg{k}$}\\
         \hline

        227,656 & 379,738 & 3.336
        \end{tabular}
    \label{tab:subgraph_data}
    \end{center}
\end{table}

Degree distribution of the subgraph can be seen in figure~\ref{fig:subgraph_deg_dist}.
It seems to follow power-law degree distribution better than the original graph.
The in degree distribution is almost textbook example of power-law, whereas out degree has some spikes that could make us question it.

\begin{figure}[h!]
    \includegraphics[width=8cm]{subgraph_deg_dist.png}
    \centering
    \caption{Degree distribution in the subgraph created from the network for both in and out degree.}
    \label{fig:subgraph_deg_dist}
\end{figure}
To reach the last node starting from the tainted one, we have to make 1660 steps.
Based on this observation, we have more than enough nodes, that require different TaintRank scores, since we know they are not equally tainted.
With this, we hope to gain some insight into nodes that are truly tainted, and make sure that nodes that are not really involved with the thief do not get tainted (too much).
In the next approaches, we will be using this subgraph to calculate TaintRank of the nodes.

\subsection{Fixed TaintRank spread}
First method we test consists of applying the TaintRank of the original node to all its successors.
With this method, as predicted, the drawback was that now we are dealing with almost half the nodes of the original network (227,656) being equally tainted.
Because of that, we can not differ if some of the nodes are more tainted than the others.
The positive side is that we did not miss any node which received tainted funds, but we also included lots of nodes that should not be tainted.

\subsection{Link weight and node value}
In this method, we spread TaintRank based on the current node value and the weight of the out-link.
With this approach, we try to take into account the amount of tainted funds that are being sent and how tainted is the address, which is sending us these funds.
We start in the root node and calculate TaintRank for each node in an iterative fashion using the following equation where $V_j$ is the value (sum of all in-links) of node j, $\Gamma_{i}$ represents the in-neighborhood of node $i$ and $w(i,j)$ represents the weight of the edge connecting i and j.

\[t_i = \sum_{j \in \Gamma_{i}} \frac{t_j\times w(j,i)}{V_j}\]

Similarly we tested the approach, where we consider $V_j$ as the sum of all out-links instead of in-links.
Reasoning behind it is, that we look only at the funds we are spreading instead of the total funds we are have received.
Following the equation, each node received TaintRank from every one of its predecessors, based on the amount received via this link and the value of the predecessor node.
Results of both approaches can be seen in figure~\ref{fig:link_weight_value}.
As we can see, the approach where we calculated the values from out-links gives a higher TaintRank, since we are only considering spent(sent) funds.
The initial drop is also less drastic than the one with in-links.
With this method, we should consider as tainted roughly first 100 nodes, right until the lines start to have linear pace.

\begin{figure}[H]
    \includegraphics[width=8cm]{link_weight.png}
    \centering
    \caption{Using link weight and node value to determine TaintRank. Displayed are top 1000 highest TaintRank scores for both approaches at determining the value of the node. Y-axis is logarithmic for better separation.}
    \label{fig:link_weight_value}
\end{figure}

The benefits of this approach are that not all nodes get tainted equally but rather proportional to the amount received.
Negatives include a dropping TaintRank as we split the stolen funds into lots of smaller pieces and spend them separately.

\subsection{Distance from the original tainted node}
In this method we want to take into account distance a certain address has from the root (thief) node.
At start, we set the TaintRank of root node to one and then set the TaintRank for other reachable nodes using iterative fashion.
In the equation we are iterating over all in-links of the currently visited node i and calculating $t_i$ based on their TaintRank and distance of node $i$ from the root node:
\[t_i = \sum_{j \in \Gamma_{i}} \frac{t_j}{distance(i, root)}\]

The results for the top 1000 highest TaintRank scores can be seen in figure~\ref{fig:link_distance}.
The top 1000 nodes are not necessarily equal to the ones found in previous method.
With this method, we can capture when a thief is splitting Bitcoins into tiny fractions, but does not make much distance from the original address.
We can see that the curve looks stair-like, where each step represents that we are further from the root node.
Considering the closeness, we could say that first 100-200 nodes gained relatively large TaintRank and should be considered, whereas the further ones matter less.

\begin{figure}[H]
    \includegraphics[width=8cm]{dist_taint.png}
    \centering
    \caption{Using distance to determine TaintRank. Displayed are top 1000 highest TaintRank scores.}
    \label{fig:link_distance}
\end{figure}

Our reasoning for this method was, the further you would be from the tainted node, the less tainted you would get.
Of course this kind of spreading leaves us vulnerable when a thief intentionally makes a long chain before spending the stolen Bitcoins.
\subsection{Combined approach}
With the above, we have seen the shortcomings of each approach.
In order to make it harder for the thief to shake off the taintedness, we try combining both approaches.
Combined approach method calculates taintedness of a node based on both methods. The equation becomes as follows:

\[t_i = \frac{\sum_{j \in \Gamma_{i}} \frac{t_j}{distance(i, root)} + \frac{t_j\times w(j,i)}{V_j}}{2}\]

Of course, because we used both methods, we have to divide final result by two.
We have also tried applying max value from both approaches instead of average between the two, but it produced similar results (around 400 nodes, that would be considered).
In figure~\ref{fig:all_approaches} we can see the differences between the previous approaches and the final, combined one.
On the chart, we can see that the amount of considerably tainted nodes is biggest with the combined approach (around 400 nodes would be considered as tainted).
Reason behind it is that it can handle edge cases, where a specific method would fail.
An example of such edge case, where the weight method would fail is when the thief would spread Bitcoins into tiny fractions across different addresses.
Each would then have a small TaintRank value.
Vice versa can be observed in the distance method, where the thief could make a very long chain of transactions and TaintRank would drop independent of the weight of links.
With this approach, we deal with such cases, where one method fails but the other succeeds.
As with all such problems, even this approach would fail in an event where the thief knows how TaintRank is being calculated and would act to trick it.
\begin{figure}[H]
    \includegraphics[width=8cm]{all_approaches.png}
    \centering
    \caption{Comparison of different approaches at calculating TaintRank. Displayed are top 1000 nodes based on highest TaintRank.}
    \label{fig:all_approaches}
\end{figure}

\subsection{PageRank like propagation}
For this method we define our taint propagation similarly to how PageRank works.
First let's define $m_i$ as all in-edges for node $n_i$, $m_i'$ as all tainted in-edges for node $n_i$, and $k_i$ and $k_i'$ as the number tainted and all out-edges for node $n_i$.\\
We define an initial taint value for every node in the network based on the ratio of total tainted in-edges and all in-edges of a node. Labeling of whether an edge is tainted was done during the data preparation phase.
\[t_i = \frac{m_i'}{m_i}\]
To calculate the actual taint value we use a formula which calculates a nodes taint based on the taints of the neighbour nodes $\Gamma_i$ with whom node $n_i$ has an in-edge
\[t_i = \sum_{e \in \Gamma_i} \frac{t_e}{k_e'}\]
This measure spreads the taint of a node over all its tainted out edges. The distinction of why we use only the tainted out edges is because the node might have legitimate transactions which happened before the ones which moved tainted funds. Such a measure not only marks nodes dealing with a high number of tainted transactions but also propagates this high taint on the nodes that received them. \\
We applied this method on the two instances of both the clustered and non-clustered networks. We also changed the number of times we recalculate the taint on a network by orders of magnitude. A single iteration takes six seconds on a higher range laptop running an eight core i7 processor and 16 GB of RAM. We use a script written in Python to calculate the ranks.\\
\subsubsection{Taint distributions}
We run the method on a non-clustered and clustered graph with one and five iterations. After calculating the values, we plot the taint distributions of these graphs. Figures \ref{fig:pagerank_1} and \ref{fig:pagerank_5} represent the distributions for non-clustered graph, while \ref{fig:pagerank_cluster_1} and \ref{fig:pagerank_cluster_5} represent the clustered graph. What we can observe from these figures is that with the increase of iterations, the taint distribution seems to resemble power-law distribution. Such a distribution is the most expressed with the clustered five iterations graph taint distribution. \\
The reason for such a distribution is difficult for us to reason about and it requires further investigation. It might be a coincidence, based on the data sample we worked on or it might have some to us unknown underlying mechanics.

\subsubsection{Nodes with the highest taint rank}
We select the best five nodes for every iteration count according to their taint rank. We then try searching for these addresses on the internet to try and reason on why they have such high ranks. This is easy for the non-clustered graph because every node is its own address. But with the clustered graph this proves to be a challenge. Not only do we have to check all the addresses in a cluster, but also because we have to reason why the clustering method works.\\
The top node with the highest taint rank in non-clustered graph with one iteration is $1eHhgW6vquBYhwMPhQ668HPjxTtpvZGPC$. After performing search on the internet, we believe address belongs to an exchange. Search on the internet lead us to believe that the address is connected to an exchange called Mt. Gox. Other four addresses seem to be burner addresses with high amounts of incoming and outgoing bitcoins. They could potentially be used by one of the mixing services.\\
Searching for the top nodes with the highest taint rank in non-clustered graph with five iterations gave us no information about the addresses' main purpose. The top address from the single iteration result was not on the list anymore while the other four addresses moved up a spot. A new address took the fifth place. Most of the addresses again look like burner addresses. To determine whether this is correct, we'll have to find a way to classify them.

\begin{figure}[h!]
\includegraphics[width=8cm]{pagerank_1.png}
\centering
\caption{PageRank like propagation for non-clustered with one iteration}
\label{fig:pagerank_1}
\end{figure}
\begin{figure}[h!]
\includegraphics[width=8cm]{pagerank_5.png}
\centering
\caption{PageRank like propagation for non-clustered with five iterations}
\label{fig:pagerank_5}
\end{figure}

\begin{figure}[h!]
\includegraphics[width=8cm]{pagerank_clustered_1.png}
\centering
\caption{PageRank like propagation for clustered with one iteration}
\label{fig:pagerank_cluster_1}
\end{figure}
\begin{figure}[h!]
\includegraphics[width=8cm]{pagerank_clustered_5.png}
\centering
\caption{PageRank like propagation for cluster with five iterations}
\label{fig:pagerank_cluster_5}
\end{figure}

\section{Discussion}
We've developed five different methods and all of them give at least slightly different results. This makes the evaluation even more difficult. Knowing whether a specific method's results are relevant, without understanding the underlying socioeconomic mechanics, makes this problem even more difficult to tackle.
Results provided by the approaches make a lot of sense.
Since we were calculating TaintRank on 45 \% nodes of the entire network, they must not be all highely tainted - the vast majority should be very mildly tainted.
This was achieved using all of the approaches.
We have seen the pros and cons of each and tried some combinations of them.
Entire computations are computationally very efficient with time complexity for the first four iterative approaches of $O(n)$.
The assigned scores give us some insight on what nodes we would like to avoid as a Bitcoin exchange company or a merchant.
If we would process these ill-gotten Bitcoins as an exchange company or vendor we might face some liabilities.
Potentially we would have to refund the person that was hacked and deal with the loss ourselves.
Such company could also be treated as a suspect in the theft.
As stated before, we should not solely rely on this rank, but it can be a first step we take.
Potential use case for the TaintRank system would be performing such scores in real time.
After an exchange got hacked, they would pinpoint the original hacker node and TaintRank would be calculated for the future nodes, preventing others from accepting funds from tainted addresses.
Same could be applied to some nodes, that are known to participate in illicit activities.
In a perfect scenario, we could rely solely on the most basic method we tested.
Giving each node that is reachable from the thief node a TaintRank of one and assume complete knowledge from all traders - so legitimate people would not accept funds from the tainted person.
But this is not feasible, so other approaches we tested would probably be applied.

\section{Conclusion}
Analyzing the Bitcoin network is a big challenge just by itself.
The ever-expanding nature makes it very difficult to process the whole network fast enough for it to stay relevant with every new block.
With the increasing popularity the number of unique addresses and transactions between them grow exponentially.
But given enough resources and sufficiently powerful processing infrastructure, this isn't something we couldn't solve. The bigger issue is the social, non-technical aspect of the network.
The questions of how to validate the results of developed methods. Different methods return different results and comparing them without any prior knowledge is next to impossible.
We believe that the methods we developed are a great step in the right direction. They show us that a taint score can be calculated for addresses. Given enough context, it can prove to be useful in determining whether an address is trustworthy. But we also believe that these methods are far from being applicable in real world scenarios where the livelihood of a crypto-currency exchange is at stake.

\section{References}