\section{Projection, boundary conditions and physical states}
Any physical property only gets contributions from {\em physical states},
{\em i.e.} states that are not annihilated by the projection into the subspace with $D_{\vec{r} }= i$ for each $\vec{r}$. If one uses semi-open boundary conditions for Kitaev's honeycomb model, {\em i.e.} open in one direction 
and periodic in the other direction, with
the honeycomb lattice terminating in a pair of arm-chair edges in the open direction, the projection into the physical subspace only affects the occupation number of boundary
orbitals constructed from the ``dangling'' $b$ fermions~\cite{Zschocke_Vojta,Pedrocchi_Chesi_Loss}.  Here, ``dangling'' refers to the fact that these $b$ fermions do not appear in any term of the fermionized Hamiltonian. As a result, canonical fermion orbitals constructed entirely from these dangling $b$ fermions can have any occupation number in the ground state(s).

Indeed, this freedom reflects a ground state degeneracy of the fermionic Hamiltonian of Kitaev's honeycomb model. The projection operator constrains the parity of (the total occupation number of canonical orbitals constructed from) $c$-fermions relative to the parity of (the
total occupation number of canonical orbitals constructed from) $b$-fermions. Since these dangling $b$ orbitals have no role in the Hamiltonian, every $c$-fermion state has a physical version with fixed degeneracy equal to $2^{N_{\mathrm{dangling}} -1}$, where $N_{\mathrm{dangling}}$ denotes the number of canonical fermion orbitals constructed from dangling $b$ fermions.
Note that {\em internal edges} generated by the presence of vacancies also contribute to
$N_{\mathrm{dangling}}$ in systems diluted with nonmagnetic vacancies.
We have repeated this analysis for the SU(2) symmetric model studied here, and find that it remains valid in our case as well. Here, we
exploit this by choosing such semi-open boundary conditions which guarantee that
the projection operation does not change the value of $\chi(T)$ we obtain from the free-fermion spectrum.





\section{Numerical results}
We obtain the density of eigenvalues of ${\mathbf K}$ using multiprecision methods
used in earlier work~\cite{Sanyal_Damle_Motrunich_PRL}, focusing on samples with semi-open boundary conditions carefully
chosen to eliminate extraneous zero modes arising from the arm-chair edges. To achieve this, we use antiperiodic boundary conditions on the fermionic problem whenever the length of the arm-chair edge is a multiple of $4$, and periodic boundary conditions whenever it is a multiple of $2$ but not a multiple of $4$. We use a compensated vacancy
ensemble that ensures that each realization has the same number of vacancies on the two sublattices, and
forbids vacancies from being on the arm-chair edges (this eliminates some finite-size corrections). Additionally, we forbid vacancies from being adjacent to each other or to share a neighbour---this ensures that our results are dominated by samples without multiple disconnected components, unlike other work~\cite{Weik_etal_PRB} that has specifically focused on zero modes of the largest percolating component of a diluted graphene lattice.


In Fig.~\ref{Fig2Appendix,Fig3Appendix,Fig4Appendix,Fig5Appendix}, we display results for the integrated density of
nonzero energy states $N$ for four representative values of $n_v$. In each case, $N(\Gamma)$ fits extremely well to the crossover form advertised
earlier. By comparing results on two different vacancy ensembles with different
correlations between vacancies, we find that the crossover exponent $y$ and crossover scale $\epsilon_c$ are quasiuniversal, depending on $n_v$ mainly through the dependence
of $\rho_0$ on $n_v$. These results closely mirror earlier
results~\cite{Sanyal_Damle_Motrunich_PRL} on the tight-binding
model for graphene diluted with vacancies. The reason is clear: As argued in Ref~\cite{Sanyal_Motrunich_Damle_PRL}, the physics of these zero modes is expected to be largely insensitive to bond-disorder. Moreover, the crossover in the density of states is controlled by the density of these zero modes. Our results reflect this insensitivity, since the presence of $\pi$ flux on certain plaquettes is merely a specific kind of bond-disorder.
{\em And so on...
So this appendix will serve as a wrapper for all of Sambuddha's numerical results in the tarball I have attached to an earlier email. The hope is that this is eventually a supplemental material to a PRL. If that doesn't work out, this Appendix becomes a Computational Details section.}



\section{Infinite-disorder fixed point}


This form of the susceptibility is a consequence of an underlying infinite-disorder fixed point. Here we present an argument for this within the framework of the strong-disorder renormalization group approach developed in Ref.~\cite{Motrunich_Damle_Huse_spinchainPRB,Motrunich_Damle_Huse_GadeWegnerPRB,Motrunich_Damle_Huse_dirtySCPRB}.
Within this approach to the spectrum of the $f$ fermions, a  pair of low-energy wavefunctions is obtained at energies $\pm \Omega = 10^{-\Gamma} K_{\mathrm{av}}$ by decimating the strongest hopping matrix element in a renormalized hopping Hamiltonian with cutoff $\Omega = 10^{-\Gamma} K_{\mathrm{av}}$. At this value of cutoff, the renormalized random hopping problem has $N(\Gamma)$ surviving sites per unit area.  $N(\Gamma)$ thus gives the integrated density of states in the interval $0 < |\epsilon| < 10^{-\Gamma}K_{\mathrm{av}}$, thereby justifying the notation used. In this problem, let $N(\Gamma)d\Gamma/\Delta(\Gamma)$ be the number of sites per unit area which are coupled by effective hopping amplitudes of strength $\tilde{K}$ such that that $\log(K_{\mathrm{av}}/\tilde{K}) \in (\Gamma - d\Gamma, \Gamma)$. Then, this RG approach gives $|dN/d\Gamma| \sim N(\Gamma)/\Delta(\Gamma)$.

For this RG picture to be reliable and self-consistent, $\Delta(\Gamma)$, which represents the width of the distribution of hopping amplitudes in the renormalized problem with cutoff $\Omega = 10^{-\Gamma}K_{\mathrm{av}}$, must grow indefinitely in the low energy limit as $\Gamma \rightarrow \infty$. Conversely, if an exact computation yields $N(\Gamma) \times |dN/d\Gamma|^{-1}$ that grows indefinitely with $\Gamma$ in a bipartite random hopping problem, we are led to conclude that the low energy spectrum is controlled by an underlying strong disorder fixed point. Using this diagnostic, we see that the modified Gade-Wegner form of the density of states corresponds to distributions of renormalized hoppings widening with RG flow as $\Delta(\Gamma) \sim \Gamma^{1-1/x}$ with $x=3/2$, {\em implying that the form of the low temperature susceptibility obtained above is controlled by an infinite-disorder fixed point}.




\end{document}

