# Supplementary Material: The universal spectral form factor for many-body
localization

Abhishodh Prakash<EMAIL_ADDRESS>International Centre for
Theoretical Sciences (ICTS-TIFR), Tata Institute of Fundamental Research,
Shivakote, Hesaraghatta Hobli, Bengaluru 560089, India  J. H. Pixley
<EMAIL_ADDRESS>Department of Physics and Astronomy, Center for
Materials Theory, Rutgers University, Piscataway, NJ 08854 USA  Manas Kulkarni
<EMAIL_ADDRESS>International Centre for Theoretical Sciences
(ICTS-TIFR), Tata Institute of Fundamental Research, Shivakote, Hesaraghatta
Hobli, Bengaluru 560089, India

###### Contents

  1. I S1. SFF for Poisson numbers
  2. II S2. Various time-scales for the SFF
  3. III S3. SFF for completely uncorrelated numbers to Poisson numbers
    1. III.1 S3a. SFF for completely uncorrelated levels (random diagonal ensemble)
    2. III.2 S3b. From Random levels to Poisson process
  4. IV S4. Disorder independence of the SFF in the MBL phase
  5. V S5. Quantifying the SFF deviations from exact result
  6. VI S6. Connection to the Berry Tabor conjecture and the $N\rightarrow\infty$ limit
    1. VI.1 S6a. Review of the statement of the conjecture
    2. VI.2 S6b. Matching the $N\rightarrow\infty$ limit for the Poisson SFF
    3. VI.3 S6c. Matching the $N\rightarrow\infty$ limit for the Random SFF

##  I S1. SFF for Poisson numbers

We now provide details of the derivation of the expressions for the spectral
form factors (SFF) of Poisson levels presented in the main text. Consider
levels from a Poisson process $\\{E_{i}\\}$ generated by summing gaps
$\\{\delta_{i}\\}$ with an exponential i.i.d

$\displaystyle P(\\{\delta_{i}\\})$ $\displaystyle=$
$\displaystyle\prod_{i}\rho(\delta_{i}),$ (S1) $\displaystyle\rho(\delta)$
$\displaystyle=$
$\displaystyle\frac{1}{\mu}\exp\left(-\frac{\delta}{\mu}\right),$ (S2)
$\displaystyle E_{i}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{i}\delta_{k}.$ (S3)

To compute the SFF, we need the joint two-point distribution
$P(E_{n},n;E_{m},m)$ i.e. the probability that the $m^{th}$ eigenvalue is
$E_{m}$ and the $n^{th}$ eigenvalue is $E_{n}$. Assuming with no loss of
generality $m>n$, this can be obtained as follows

$\displaystyle P(E_{n},n;E_{m},m)=p(E_{n},n)~{}p(E_{m}-E_{n},m-n),$ (S4)

where $p(E_{k},k)$ is the well known Poisson distribution

$p(E_{k},k)=\frac{e^{-\frac{E_{k}}{\mu}}}{\mu(k-1)!}\left(\frac{E_{k}}{\mu}\right)^{k-1}.$
(S5)

We assume $E_{k}>0$ and $k=1,2,3,\ldots$. Let us sketch the derivation of the
above equation

$\displaystyle P(E_{n},n;E_{m},m)$
$\displaystyle=\int_{0}^{\infty}d\delta_{1}\rho(\delta_{1})\cdots\int_{0}^{\infty}d\delta_{N}\rho(\delta_{N})~{}\delta\left(E_{m}-\sum_{i=1}^{m}\delta_{i}\right)~{}\delta\left(E_{n}-\sum_{i=1}^{n}\delta_{i}\right)$
$\displaystyle=\int_{0}^{\infty}d\delta_{1}\rho(\delta_{1})\cdots\int_{0}^{\infty}d\delta_{N}\rho(\delta_{N})~{}\int_{-\infty}^{\infty}\frac{dk}{2\pi}~{}e^{ik\left(E_{m}-\sum_{i=1}^{m}\delta_{i}\right)}~{}\int_{-\infty}^{\infty}\frac{dk^{\prime}}{2\pi}~{}e^{ik^{\prime}\left(E_{n}-\sum_{i=1}^{n}\delta_{i}\right)}$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dk}{2\pi}\int_{-\infty}^{\infty}\frac{dk^{\prime}}{2\pi}~{}e^{i(kE_{m}+k^{\prime}E_{n})}\left[\prod_{l=1}^{n}\int_{0}^{\infty}d\delta_{l}~{}\rho(\delta_{l})~{}e^{-i(k+k^{\prime})\delta_{l}}\right]\left[\prod_{l={n+1}}^{m}\int_{0}^{\infty}d\delta_{l}~{}\rho(\delta_{l})~{}e^{-ik\delta_{l}}\right]$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dk}{2\pi}\int_{-\infty}^{\infty}\frac{dk^{\prime}}{2\pi}~{}e^{i(kE_{m}+k^{\prime}E_{n})}\left[\langle
e^{-i(k+k^{\prime})\delta}\rangle_{\rho}\right]^{n}~{}\left[\langle
e^{-ik\delta}\rangle_{\rho}\right]^{m-n}.$ (S6)

We now change variables $\\{k,k^{\prime}\\}\mapsto\\{k,l=k+k^{\prime}\\}$

$\displaystyle P(E_{n},n;E_{m},m)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dk}{2\pi}~{}e^{ik(E_{m}-E_{n})}\left[\langle
e^{-ik\delta}\rangle_{\rho}\right]^{m-n}\int_{-\infty}^{\infty}\frac{dl}{2\pi}~{}e^{ilE_{n}}\left[\langle
e^{-il\delta}\rangle_{\rho}\right]^{n}$
$\displaystyle=p(E_{n},n)~{}p(E_{m}-E_{n},m-n),$ (S7)

where

$\displaystyle\langle e^{-ik\delta}\rangle_{\rho}$
$\displaystyle=\int_{0}^{\infty}d\delta~{}\rho(\delta)~{}e^{-ik\delta}=\frac{1}{1+i\mu
k},$ (S8) $\displaystyle p(E_{n},n)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dk}{2\pi}~{}e^{ikE_{n}}\left[\langle
e^{-ik\delta}\rangle_{\rho}\right]^{n},$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dk}{2\pi}~{}\frac{e^{ikE_{n}}}{(1+i\mu
k)^{n}}=\frac{e^{-\frac{E_{n}}{\mu}}}{\mu(n-1)!}\left(\frac{E_{n}}{\mu}\right)^{n-1}.$
(S9)

In the last step, the integral is performed by closing the contour in the
upper-half complex plane to enclose the $n^{th}$ order pole at $k=i/\mu$ and
using Cauchy’s integral formula.

Let us now proceed to compute the SFF

$\displaystyle K(\tau,N)$ $\displaystyle=\sum_{m,n=1}^{N}\langle
e^{i\tau(E_{m}-E_{n})}\rangle_{P}=\sum_{m=n=1}^{N}1+\sum_{m\neq
n=1}^{N}\langle e^{i\tau(E_{m}-E_{n})}\rangle_{P}$
$\displaystyle=N+\sum_{n=1}^{N}\sum_{m=n+1}^{N}\left(\langle
e^{i\tau(E_{m}-E_{n})}\rangle_{P}+\langle
e^{-i\tau(E_{m}-E_{n})}\rangle_{P}\right),$ (S10)

where

$\displaystyle\langle e^{i\tau(E_{m}-E_{n})}\rangle_{P}$
$\displaystyle=\int_{0}^{\infty}dE_{m}~{}\int_{0}^{\infty}dE_{n}~{}P(E_{n},n;E_{m},m)~{}e^{i\tau(E_{m}-E_{n})}$
$\displaystyle=\int_{0}^{\infty}dE_{n}~{}p(E_{n},n)\int_{0}^{\infty}dE_{m}~{}p(E_{m}-E_{n},m-n)~{}e^{i\tau(E_{m}-E_{n})}.$
(S11)

Changing the integration variables
$\\{E_{m},E_{n}\\}\mapsto\\{E_{mn}=E_{m}-E_{n},E_{n}\\}$, we get

$\displaystyle\langle e^{i\tau(E_{m}-E_{n})}\rangle_{P}$
$\displaystyle=\int_{0}^{\infty}dE_{n}~{}p(E_{n},n)\int_{0}^{\infty}dE_{mn}~{}p(E_{mn},m-n)~{}e^{i\tau
E_{mn}}$ $\displaystyle=\int_{0}^{\infty}dE_{mn}~{}p(E_{mn},m-n)~{}e^{i\tau
E_{mn}}=\frac{1}{(1-i\mu\tau)^{m-n}},$ (S12)

where we have used $\int_{0}^{\infty}dE_{n}~{}p(E_{n},n)=1$. Substituting the
above expression in eq. S10 gives

$\displaystyle K(\tau,N)$
$\displaystyle=N+\sum_{n=1}^{N}\sum_{m=n+1}^{N}\left[\frac{1}{(1+i\mu\tau)^{m-n}}+\frac{1}{(1-i\mu\tau)^{m-n}}\right].$
(S13)

The summand of the above series depends only on the differences $k=m-n$. As a
result, we can change the summation over $m$ and $n$ in terms of $k$ as
follows

$\displaystyle
K(\tau,N)=N+\sum_{k=1}^{N-1}(N-k)\left[\frac{1}{(1+i\mu\tau)^{k}}+\frac{1}{(1-i\mu\tau)^{k}}\right].$
(S14)

The above series can be summed up easily using standard results from summing
the geometric series

$\displaystyle
S^{G}_{n}=\sum_{k=1}^{n}br^{k-1}=b+br+br^{2}+\ldots+br^{n-1}=\frac{b\left(1-r^{n}\right)}{\left(1-r\right),}$
(S15)

and the arithmetico-geometric series

$\displaystyle S_{n}$
$\displaystyle=\sum_{k=1}^{n}\left(a+(k-1)d\right)br^{k-1}=ab+(a+d)br+(a+2d)br^{2}+\ldots+(a+(n-1)d)br^{n-1}$
$\displaystyle=\frac{ab}{1-r}+\frac{dbr(1-r^{n})}{(1-r)^{2}}-\frac{(a+nd)br^{n}}{(1-r)},$
(S16)

with $a=d=1$, $b=r=\frac{1}{(1\pm i\mu\tau)}$ and $n=N-1$ to get the final
answer

$K(\tau,N)=N+\frac{2}{(\mu\tau)^{2}}-\frac{(1+i\mu\tau)^{1-N}+(1-i\mu\tau)^{1-N}}{(\mu\tau)^{2}}.$
(S17)

We now proceed to calculate the connected SFF

$K_{c}(\tau,N)=K(\tau,N)-|\sum_{m}\langle e^{i\tau E_{m}}\rangle_{P}|^{2},$
(S18)

where

$\displaystyle\langle e^{i\tau E_{m}}\rangle_{P}$
$\displaystyle=\int_{0}^{\infty}dE_{m}~{}p(E_{m},m)~{}e^{i\tau
E_{m}}=\frac{1}{(1-i\mu\tau)^{m}},$ (S19) $\displaystyle\sum_{m=1}^{N}\langle
e^{i\tau E_{m}}\rangle_{P}$
$\displaystyle=\sum_{m=1}^{N}\frac{1}{(1-i\mu\tau)^{m}}=\frac{(1-i\mu\tau)^{-N}-1}{i\mu\tau}.$
(S20)

Plugging this into eq. S18 and some simplification, we get the final answer

$K_{c}(\tau,N)=N+\frac{1}{(\mu\tau)^{2}}-\frac{(1+(\mu\tau)^{2})^{-N}}{(\mu\tau)^{2}}-\frac{i}{\mu\tau}\left[(1+i\mu\tau)^{-N}-(1-i\mu\tau)^{-N}\right].$
(S21)

Figure S1 shows the SFF computed numerically for $N$ Poisson levels and the
excellent agreement with the above analytical result. Note that the expression
for the SFF was obtained by the authors of [1] as a special case of a more
general result for uncorrelated gap distributions.

|
---|---
Figure S1: The connected SFF (left) and SFF (right) for $N=500$ numerically
generated Poisson levels with unit mean level spacing averaged over 50,000
realizations. The analytic result is shown as the black line.

##  II S2. Various time-scales for the SFF

Let us now closely analyze the various time-scales of the SFF for Poisson
levels. Let us focus on the connected SFF for concreteness and clarity. We
also assume that $N\gg 1$ although finite. Let us start with the expression
for $K^{P}_{c}(\tau,N)$

$K^{P}_{c}(\tau,N)=N+\frac{1}{(\mu\tau)^{2}}-\frac{(1+(\mu\tau)^{2})^{-N}}{(\mu\tau)^{2}}-\frac{i}{\mu\tau}\left[(1+i\mu\tau)^{-N}-(1-i\mu\tau)^{-N}\right]$
(S22)

For large values of $\mu\tau\gg 1$, we have

$\displaystyle K^{late}_{c}(\tau,N)\approx N.$ (S23)

For $\frac{1}{\sqrt{N}}<\mu\tau<1$, we have the result stated in the main-text

$\displaystyle K^{inter}_{c}(\tau,N)\approx N+\frac{1}{(\mu\tau)^{2}}.$ (S24)

For $\mu\tau<\frac{1}{\sqrt{N}}$, we can get the form of the SFF by a careful
Taylor expansion as follows

$\displaystyle K^{early}_{c}(\tau,N)$
$\displaystyle=N+\frac{1}{(\mu\tau)^{2}}-\frac{(1+(\mu\tau)^{2})^{-N}}{(\mu\tau)^{2}}-2\left(1+(\mu\tau)^{2}\right)^{-\frac{N}{2}}\frac{\sin\left(N\tan^{-1}(\mu\tau)\right)}{(\mu\tau)}$
$\displaystyle\approx
N+\frac{1}{(\mu\tau)^{2}}-\frac{(1-N(\mu\tau)^{2})}{(\mu\tau)^{2}}-2\left(1-\frac{N}{2}(\mu\tau)^{2}\right)\frac{\sin\left(N\mu\tau\right)}{(\mu\tau)}$
$\displaystyle\approx
2N\left(1-\frac{\sin\left(N\mu\tau\right)}{(N\mu\tau)}\right).$ (S25)

In this small $\tau$ regime, oscillations with a period $T=2\pi/(N\mu)$
appear.

|
---|---
Figure S2:  The connected SFF with (left) the early and intermediate forms
separated and (right) the various $\tau$ regimes separated .

The early-$\tau$ SFF can be further separated into two regions. For
$\frac{1}{N}<\mu\tau<\frac{1}{\sqrt{N}}$, we have $K_{c}^{early}\approx 2N$.
For $\mu\tau<\frac{1}{N}$, the sine term can be Taylor expanded to give us
$K^{early}_{c}\approx\frac{N^{3}(\mu\tau)^{2}}{3}$. The quadratic and constant
pieces meet at $\mu\tau=\frac{\sqrt{6}}{N}$. Putting all these together, we
see that for finite but large N, the connected SFF can be divided into four
$\tau$ regimes as shown in fig. S2.

Note that similar time-scales as well as scaling forms were also obtained in
[1] although their notion of universality (independence of underlying gap
distribution) is different from ours. While we consider the early $\tau$
regime to be non-universal since it is affected by the global density of
states, [1] considers the same to be universal since the expressions assume a
form that is independent of the underlying distribution of the uncorrelated
gaps.

##  III S3. SFF for completely uncorrelated numbers to Poisson numbers

The often-quoted relationship between Poisson processes and uncorrelated
numbers is clarified in this section. First, we derive the SFF for
uncorrelated numbers and demonstrate that it is different from the form
obtained previously in eqs. S17 and S21. Next, we present the setting when the
SFF for uncorrelated numbers reduces to that of the Poisson process.

###  III.1 S3a. SFF for completely uncorrelated levels (random diagonal
ensemble)

Consider an $N\times N$ matrix with random entries along the diagonal drawn
from some independent and identical distribution (i.i.d) $p(E)$. The
eigenvalues $\\{E_{i}\\}$ of such a matrix is simply the diagonal entries
sorted. As can be easily seen from the formulas for the SFF,

$\displaystyle K(\tau,N)$
$\displaystyle=\langle\sum_{m,n=1}^{N}e^{i\tau(E_{m}-E_{n})}\rangle,$ (S26)
$\displaystyle K_{c}(\tau,N)$
$\displaystyle=\langle\sum_{m,n=1}^{N}e^{i\tau(E_{m}-E_{n})}\rangle-|\sum_{m=1}^{N}\langle
e^{i\tau E_{m}}\rangle|^{2},$ (S27)

the ordering of $\\{E_{i}\\}$ is irrelevant to compute the SFF. We need the
joint two-point distribution $P(E_{n},n;E_{m},m)$, i.e. probability that the
$m^{th}$ eigenvalue is $E_{m}$ and the $n^{th}$ eigenvalue is $E_{n}$. Since
we can ignore ordering, this is simply

$P(E_{n},n;E_{m},m)=p(E_{n})~{}p(E_{m}).$ (S28)

Let us start with the disconnected SFF for random levels, $K^{R}(\tau,N)$

$\displaystyle K^{R}(\tau,N)$ $\displaystyle=\sum_{m,n}\langle
e^{i\tau(E_{m}-E_{n})}\rangle=\sum_{m=n}1+\sum_{m\neq n}\langle e^{i\tau
E_{m}}e^{-i\tau E_{n}}\rangle_{P}=N+N(N-1)|\langle e^{i\tau
E}\rangle_{p}|^{2},$ (S29)

where we have used

$\displaystyle\langle e^{i\tau E_{m}}e^{-i\tau E_{n}}\rangle_{P}$
$\displaystyle=$ $\displaystyle\langle e^{i\tau E_{m}}\rangle_{p}\langle
e^{-i\tau E_{n}}\rangle_{p},$ (S30) $\displaystyle\text{and }\langle e^{i\tau
E_{m}}\rangle$ $\displaystyle=$ $\displaystyle\langle e^{i\tau
E_{n}}\rangle_{p}=\langle e^{i\tau
E}\rangle_{p}=\int_{0}^{2D}dE~{}p(E)~{}e^{i\tau E}.$ (S31)

We can also get the connected SFF easily

$K^{R}_{c}(\tau,N)=K(\tau,N)-|\sum_{m}\langle e^{i\tau
E_{m}}\rangle|^{2}=K(\tau,N)-N^{2}|\langle e^{i\tau
E}\rangle|^{2}=N(1-|\langle e^{i\tau E}\rangle|^{2}),$ (S32)

|
---|---
Figure S3: The connected (left) and disconnected (right) SFF for random levels
with unit mean level spacing. The analytic result is the black line for
$N=500$ levels, which we compare with the numerically computed SFF using
levels drawn from a uniform distribution that is averaged over 50,000
realizations.

As a concrete example, if the levels are drawn from a uniform distribution of
width $D$ i.e.

$\displaystyle p(E)$ $\displaystyle=\begin{cases}\frac{1}{D}~{}for~{}0\leq
E\leq D\\\ 0~{}~{}otherwise\end{cases},$ (S33)

the expressions for SFF are as follows

$\displaystyle K^{R}(\tau,N)$
$\displaystyle=N+N(N-1)\left|\frac{\sin(D\tau/2)}{(D\tau/2)}\right|^{2},$
(S34) $\displaystyle K^{R}_{c}(\tau,N)$
$\displaystyle=N\left(1-\left|\frac{\sin(D\tau/2)}{(D\tau/2)}\right|^{2}\right).$
(S35)

It is more convenient to express these results in terms of the number of
levels $N$ and the mean level spacing $\mu$ which are related as $D=\mu N$.

$\displaystyle K^{R}(\tau,N)$ $\displaystyle=N+N(N-1)\left|\frac{\sin(\mu
N\tau/2)}{(\mu N\tau/2)}\right|^{2},$ (S36) $\displaystyle K^{R}_{c}(\tau,N)$
$\displaystyle=N\left(1-\left|\frac{\sin(\mu N\tau/2)}{(\mu
N\tau/2)}\right|^{2}\right).$ (S37)

Figure S3 shows a comparison of the analytic result for the SFF for random
levels and with it, the numerically computed SFF for $N$ random numbers drawn
from a uniform distribution with width $D=N$ (to ensure unit mean level
spacing) and its excellent matching with the above analytical expression.
Observe that $K_{c}(\tau,N)$ in eq. S32 strictly satisfies the inequality.

$\displaystyle K^{R}_{c}(\tau,N)\leq N$ (S38)

Recall that for Poisson levels, for large $N$ and intermediate values of
$\tau$ we have

$\displaystyle K^{P}_{c}(\tau,N)\approx N+\frac{1}{(\mu\tau)^{2}}\geq N.$
(S39)

Thus, the sub leading universal $\frac{1}{(\mu\tau)^{2}}$ form, which violates
the inequality, is impossible to be present in $K^{R}_{c}(\tau,N)$. We
conclude that the SFF for completely uncorrelated random levels is different
from the SFF for Poisson numbers, $K^{P}_{c}(\tau,N)$!

###  III.2 S3b. From Random levels to Poisson process

As we have seen, the SFF for numbers drawn from a Poisson process is different
from that of uncorrelated levels. In this section, we review the relationship
between uncorrelated levels and Poisson processes. We specify the condition
under which the SFF of the former reduces to the latter.

Let us once again consider $N$ uncorrelated random levels drawn from any
distribution. Let $E_{1},E_{2},\ldots,E_{N}$ be the same random numbers after
sorting. The process of sorting introduces a weak correlation that can be
observed in local spectral statistics. For example, the mean level spacing
distribution,

$\langle
P(s,N)\rangle=\frac{1}{N}\sum_{i}\langle\delta(s-E_{i+1}+E_{i})\rangle,$ (S40)

approaches the Poisson one

$\langle P(s,N)\rangle\xrightarrow{N\gg
1}P(s)=\frac{1}{\mu}e^{-\frac{s}{\mu}},$ (S41)

where $\mu$ is inherited from the original distribution of the random numbers.
This also extends to the distribution of larger spacings i.e.

$\langle
P(k,s,N)\rangle=\frac{1}{N}\sum_{i}\langle\delta(s-E_{i+k}+E_{i})\rangle\xrightarrow{N\gg
k}P(k,s)=\frac{s^{k-1}}{\mu^{k}(k-1)!}e^{-\frac{s}{\mu}}.$ (S42)

|
---|---
Figure S4: $\langle P(k,s,N)\rangle$ for various values of $k$ numerically
computed for $N=500$ sorted random numbers drawn from a uniform distribution
with mean level spacing $\mu=1$ averaged over $50,000$ samples. It can be seen
that $\langle P(k,s,N)\rangle$ agrees with the Poisson distribution, $P(k,s)$
(dashed lines) for small values of $k$ but deviates increasingly when $k\sim
N$.

Note that crucially, the above is valid only for $k\ll N$. For $k\sim N$ on
the other hand, the distribution of $\langle P(k,s,N)\rangle$ would depend on
the underlying i.i.d distribution of the uncorrelated numbers. These facts can
be easily checked numerically as shown in fig. S4. In the limit
$N\rightarrow\infty$, we can expect the Poisson distribution for any finite
$k$. We can now understand why the SFF for random levels is different from
that of numbers from a Poisson process. Given uncorrelated random levels, the
SFF is computed using _all_ differences in eigenvalues $E_{i}-E_{j}$. When
$|i-j|\ll N$, the levels are Poissonian but when $|i-j|\sim N$, they are not.
As a result, the total SFF, that sums over _all_ differences $E_{i}-E_{j}$
deviates for Random levels from the Poisson form. However, we can recover the
Poisson SFF in the following way. Consider an ensemble of $N_{R}$ uncorrelated
random numbers and sort them (eigenvalues of $N_{R}\times N_{R}$ random
diagonal ensemble). From each sample, let us then choose only $N$ levels (from
the middle, say) out of $N_{R}$ levels to compute the SFF. If $N\ll N_{R}$,
the levels resemble a Poisson process in the sense all differences
$E_{i}-E_{j}$ follow a Poisson distribution and we should recover the Poisson
SFF. As $N_{R}\rightarrow\infty$, we should get the Poisson SFF for any $N$.
This can indeed be checked to be true numerically and shown in fig. S5 for the
connected SFF where the evolution from random to Poisson SFF can be clearly
seen. Starting from $N=N_{R}$ where the SFF matches with the form eq. S37,
with increasing $N_{R}$, the SFF form moves closer to the Poisson one.

These results can be summarized by taking two different limits of the spectral
form factor, namely for $N_{R}\gg 1$ we have

$\lim_{N\rightarrow N_{R}}K(\tau,N)\rightarrow
K^{R}(\tau,N)=N+N(N-1)\left|\frac{\sin(D\tau/2)}{(D\tau/2)}\right|^{2}$ (S43)

where $N=N_{R}=D$ and

$\lim_{N/N_{R}\rightarrow 0}K(\tau,N)\rightarrow
K^{P}(\tau,N)=N+\frac{2}{(\mu\tau)^{2}}-\frac{(1+i\mu\tau)^{1-N}+(1-i\mu\tau)^{1-N}}{(\mu\tau)^{2}}.$
(S44)

As we show below in Sec. S6, in the thermodynamic limit, i.e.
$N\rightarrow\infty$ both of these limits reduce to the conjecture of Berry
and Tabor for integrable systems [2, 3], namely

$\lim_{N_{R}\rightarrow\infty}\lim_{N\rightarrow
N_{R}}\frac{1}{N}K(\tau,N)=\lim_{N\rightarrow\infty}\lim_{N_{R}\rightarrow\infty}\frac{1}{N}K(\tau,N)=1+2\pi\delta(\tau).$
(S45)

---
Figure S5: The evolution of the connected SFF from Random eq. S37 to Poisson
eq. S21. $N=1000$ levels chosen from an ensemble of sorted $N_{R}$ numbers
with drawn from a uniform distribution with unit mean level spacing and
averaged over $10^{5}$ disorder realizations. .

We conclude this section with an additional comment on taking the
thermodynamic limit. In order to extract the Poisson behavior as
$N\rightarrow\infty$, especially the sub leading $\frac{1}{(\mu\tau)^{2}}$
signature for $K_{c}(\tau,N)$ we must respect the order of limits: first we
should take $N_{R}\rightarrow\infty$ _and then_ take $N\rightarrow\infty$. In
fact, the procedure to extract Poisson statistics from model MBL Hamiltonians
also needs similar care. Oganesyan and Huse [4] point out that given a
disordered Hamiltonian with $L$ sites, the number of random numbers used to
describe the microscopic model are $\mathcal{O}(L)$. However, the number of
eigenvalues produced are $\mathcal{O}(2^{L})$. This means that there would be
significant correlations between them. However, if we consider $N$ eigenvalues
from deep within the spectrum, as we have in the main text, the Poisson nature
is revealed if $N\ll\mathcal{N}_{L}$ where $\mathcal{N}_{L}$ is the total size
of the Hilbert space. To observe Poisson statistics for arbitrary $N$, we need
to take $L\rightarrow\infty$ first.

##  IV S4. Disorder independence of the SFF in the MBL phase

We now focus on the disorder dependence of the SFF in the model of the
disordered spin-chain that we considered in the main text. We list it here
again for completeness

$\displaystyle H(W)$
$\displaystyle=\sum_{i}J_{1}(S^{x}_{i}S^{x}_{i+1}+S^{y}_{i}S^{y}_{i+1}+\Delta
S^{z}_{i}S^{z}_{i+1})+w_{i}S^{z}_{i}$
$\displaystyle+\sum_{i}J_{2}(S^{x}_{i}S^{x}_{i+2}+S^{y}_{i}S^{y}_{i+2}+\Delta
S^{z}_{i}S^{z}_{i+2}).$ (S46)

As described in the main text, by keeping track of the adjacent gap ratio $r$
that is defined in terms successive gaps $\delta_{i}=E_{i+1}-E_{i}$ of the
energy spectrum $\\{E_{i}\\}$ as

$r_{i}=\frac{\mathrm{min}(\delta_{i},\delta_{i+1})}{\mathrm{max}(\delta_{i},\delta_{i+1})},$
(S47)

we can determine that the model above has a thermal phase for weak disorder
($W<W_{c}$) and an MBL phase at strong disorder ($W>W_{c}$). The critical
disorder is roughly around $W_{c}\approx 7.3$ where the curves for different
system sizes appear to cross as shown in fig. S6.

---
Figure S6: Adjacent gap ratio for the disordered Hamiltonian eq. S46 for three
system sizes, showing a clear crossing in the vicinity of the MBL transition
(this is shown as the inset in fig. 1 of the main text).

In the main text, we focused on data that is deep in the MBL regime (i.e.
where $\langle r\rangle\approx 0.39$ is nicely Poisson), for a specific
disorder strength $W=25$ where we showed that the SFF converges well with the
analytical result. We now test the validity of this result _across_ the MBL
phase. Figure S7 shows that for a wide range of disorder strengths $W\geq 10$,
for small values of $N=20,40,80$ the SFF is remarkably convergent with each
other and, importantly, with the analytical result. For larger values of
$N=900$, the universal $1/(\mu\tau)^{2}$ behavior can be observed at
intermediate $\tau$ values. This further strengthens the claim that MBL is
indeed a robust phase.

|
---|---
Figure S7: For system size $L=14$, a comparison of (left) the SFF for various
disorder strengths in the MBL phase of the disordered Hamiltonian eq. S46 for
small values of $N=20,40,80$ and (right) the reduced SFF for $N=900$. Up to
numerical uncertainty, we find remarkable agreement with the analytic result
across all of our numerical data for $W\geq 10$.

In table 1, we list the number of disorder samples used for the various
numerical calculations in obtaining the results presented in the main text as
well as the supplementary materials.

H(W) | L=12 | L=14 | L=16 | L=18
---|---|---|---|---
W $<$ 25 | 10,000 | 10,000 | 10,000 | -
W = 25 | 50,000 | 50,000 | 90,000 | 17,500

$H_{lbit}$ | L=18 | L=20 | L=22 | L=24 | L=26
---|---|---|---|---|---
| 50,000 | 50,000 | 50,000 | 50,000 | 12,000

Table 1: Number of disorder samples analyzed in studying (left) the disordered
spin chain model eq. S46 and (right) the truncated $l$-bit model eq. S51

##  V S5. Quantifying the SFF deviations from exact result

At finite sizes, deviations from the exact result due the deviation of the
many-body density of states (DOS) from the Poisson form are inevitable at
early values of $\tau$. This is minimized by picking out $N$ eigenvalues
deeper in the middle of the spectrum where the DOS is closest to Poisson. To
examine this, we compute the many-body DOS defined as

$\rho(E)=\frac{1}{N}\sum_{i=1}^{N}\delta(E-E_{i}),$ (S48)

in comparison with the DOS of Poisson numbers. The latter can be computed
exactly as follows

$\displaystyle\langle\rho(E)\rangle$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\langle\delta(E-E_{k})\rangle=\frac{1}{N}\sum_{k=1}^{N}\int_{0}^{\infty}dEP(k,E_{k})~{}\delta(E-E_{k})$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\left(\frac{E}{\mu}\right)^{k-1}\frac{e^{-\frac{E}{\mu}}}{\mu(k-1)!}=\frac{\Gamma\left(N,E/\mu\right)}{\mu~{}N!},$
(S49)

where, $\Gamma(N,x)$ is the _incomplete Gamma function_ defined (for integer
$N$) as

$\displaystyle\Gamma(N,x)=\int_{x}^{\infty}dt~{}t^{N-1}e^{-t}=e^{-x}~{}(N-1)!\sum_{k=0}^{N}\frac{x^{k}}{k!}.$
(S50)

As shown in fig. S8, we observe that for $N$ corresponding to increasingly
smaller fractions of the total Hilbert space dimension in the $S^{z}=0$
sector, $\mathcal{N}_{L}=\binom{L}{L/2}$, the DOS increasingly approaches the
analytical form computed for Poisson numbers eq. S49 for both the disordered
spin chain model $H(W)$ of eq. S46 in the MBL phase ($W=25$) as well as the
truncated $l$-bit model eq. S51 (we list here again the $l$-bit model for
completeness).

$\displaystyle H_{lbit}$
$\displaystyle=\sum_{i=1}^{L}\sum_{a=1}^{10}J^{a}_{i}\kappa^{z}_{i}\kappa^{z}_{i+1}\ldots\kappa^{z}_{i+a}$
(S51)

Despite minimizing this deviation by selecting the eigenvalues from the middle
of the many-body spectrum, differences are seen at the edge of the selected
eigenvalues. The SFF results at small $\tau$ values are dominated by large
energy differences, which probe the edges of the many-body spectrum where
there are bound to be deviations at finite sizes. On the other hand, for
intermediate values of $\tau$, the SFF probes the eigenvalues at smaller
scales where the DOS of the Poisson numbers is flat and we should see good
agreement with theoretical expectations, including the $\frac{1}{\tau^{2}}$
behavior.

|  |
---|---|---
|  |
Figure S8: The many-body density of states $\rho(E)$ for the truncated $l$-bit
model eq. S51 (top) with the disordered Hamiltonian eq. S46 with $W=25$
(bottom) for different choices for eigenvalues $N$ in the $S^{z}=0$ sector,
$\mathcal{N}_{L}=\binom{L}{L/2}$ for various system sizes $L$, averaged over
10000 disorder samples. The black dashed line is the expectation for a Poisson
density of states given in Eq. (S49), which agrees well with our numerical
results in the limit $N/\mathcal{N}_{L}\rightarrow 0$.

The DOS curvature and thus the deviation from the analytical result are
expected to reduce by increasing system size ($L$ and thus the size of the
Hilbert space $\binom{L}{L/2}$) and reducing $N$. This is clearest in the
connected SFF as shown in fig. S9. Consistent with the expectations stated
above, we find the largest area of deviation occurs at short $\tau$ values,
where as larger $\tau$ values (but well before the plateau sets in) that are
sampling energy differences from the flat DOS, our analytic results match the
numerics nicely.

|
---|---
Figure S9: Comparison of the analytic result for the SFF of Poisson levels for
various values of $L$ and $N$ with the numerics for the $l$-bit model (left)
and the disordered Hamiltonian (right).

To quantify the deviation, we consider the quantity $\Delta$ defined as the
root mean square difference of the numerical connected SFF from the analytical
value at various values of $\tau$ defined as

$\displaystyle\Delta=\sqrt{\frac{\sum_{i=1}^{N_{\tau}}(K_{C}(\tau_{i})-K^{P}_{C}(\tau_{i}))^{2}}{N_{\tau}}},$
(S52)

$\tau_{i}$ is chosen from the interval $\tau_{i}\in(10^{-5},15.0)$ (with mean-
level spacing set to 1) where the SFF is non-trivial and $N_{\tau}$ is the
total number of $\tau$ points sampled. As seen in fig. S10, we see $\Delta$
decrease as expected with increasing $L$ as well as decreasing $N$.

|
---|---
Figure S10: Root mean square deviation of the connected SFF from analytical
result for $l$-bit (left) and Hamiltonian (right).

##  VI S6. Connection to the Berry Tabor conjecture and the
$N\rightarrow\infty$ limit

###  VI.1 S6a. Review of the statement of the conjecture

Berry and Tabor [2, 3] were interested in the spectra of point particles which
have classically integrable geodesics. For examples, free particles in a
manifold $\mathcal{M}$ with positive curvature such as spheres and tori or
with boundary conditions such as circular or rectangular. The spectrum is
basically that of the Schrodinger operator which is

$\displaystyle-\frac{\hbar^{2}}{2m}\triangle\phi_{i}=\lambda_{i}\phi_{i}$
(S53)

where $\triangle$ is the Laplace-Beltrami operator on $\mathcal{M}$ and the
eigenvalues $\lambda_{i}$ are positive definite, countably infinite and
unbounded from above $0<\lambda_{1}<\lambda_{2}<\ldots\rightarrow\infty$. The
question that the authors were interested is whether the distribution of
spectral gaps assumes an asymptotic form. Specifically, if we define

$\displaystyle P(s,N)=\frac{1}{N}\sum_{i}\delta(s-\lambda_{i+1}+\lambda_{i}),$
(S54)

the question is whether the following limit exists and what its form is

$P(s,N)\xrightarrow{N\rightarrow\infty}P(s).$ (S55)

Berry and Tabor conjecture that the limit does exist for sufficiently generic
systems and corresponds to that of a Poisson process

$\displaystyle P(s)=\frac{1}{\mu}e^{-\frac{s}{\mu}}.$ (S56)

The Berry Tabor conjecture is sometimes also stated in terms of the pair
density correlator

$\displaystyle
R_{2}(s)=\frac{1}{N}\sum_{i,j=1}^{N}\delta(s-E_{j}+E_{i})\xrightarrow{N\rightarrow\infty}\delta(s)+1,$
(S57)

and its Fourier transform, the spectral form factor

$\displaystyle\frac{K(\tau,N)}{N}=\frac{1}{N}\sum_{m,n=1}^{N}\langle
e^{i\tau(E_{m}-E_{n})}\rangle_{P}\xrightarrow{N\rightarrow\infty}1+2\pi\delta(\tau).$
(S58)

###  VI.2 S6b. Matching the $N\rightarrow\infty$ limit for the Poisson SFF

A consequence of the $l-bit$ hypothesis states that the MBL spectrum, due to
emergent integrability should also satisfy the Berry-Tabor conjecture. This is
supported by various measures of spectral correlations and statistics being
reproduced by Poisson numbers. We now show that the spectral form factor of
eq. S17 reduces to eq. S58 in the appropriate limit. Let us start with eq. S17
(setting $\mu=1$ for convenience)

$\displaystyle K^{P}(\tau,N)$
$\displaystyle=N+\frac{2}{\tau^{2}}-\frac{(1+i\tau)^{1-N}+(1-i\tau)^{1-N}}{\tau^{2}}.$
(S59)

We wish to recover the limit for integrable models eq. S58

$\lim_{N\rightarrow\infty}\frac{K^{P}(\tau,N)}{N}=1+2\pi\delta(\tau).$ (S60)

This amounts to proving that

$\displaystyle\lim_{N\rightarrow\infty}\eta_{N}(\tau)$
$\displaystyle=2\pi\delta(\tau)$ (S61) $\displaystyle\text{ where,
}\eta_{N}(\tau)$
$\displaystyle=\frac{2}{N\tau^{2}}-\frac{(1+i\tau)^{1-N}+(1-i\tau)^{1-N}}{N\tau^{2}}$
(S62)

It is very clear from the plots of $\eta_{N}(\tau)$ fig. S11 that the
functional form looks increasingly like the delta function with increasing
$N$.

---
Figure S11: $\eta_{N}(\tau)$ for various $N$

To prove the identity eq. S58, we need to find the limiting form of
$K^{P}(\tau,N)$ when viewed as an integration measure i.e. we want to show
that, for $f(\tau)$ that is well-behaved on the real line
$\tau\in{\mathbb{R}}$, we need to show that

$\displaystyle\lim_{N\rightarrow\infty}\int_{-\infty}^{\infty}d\tau
f(\tau)\eta_{N}(\tau)=2\pi f(0).$ (S63)

To do this, we need to prove that $\eta_{N}(\tau)$ satisfies three conditions
[5]:

  1. 1.

Symmetry: $\eta_{N}(-\tau)=\eta_{N}(\tau)$

  2. 2.

Singularity: $\lim_{N\rightarrow\infty}\frac{\eta_{N}(\tau\neq
0)}{\lim_{\tau\rightarrow 0}\eta_{N}(\tau)}=0$

  3. 3.

Normalization: $\lim_{N\rightarrow\infty}\int d\tau\eta_{N}(\tau)=2\pi$

The first i.e. symmetry is verified easily from the form of $\eta_{N}(\tau)$
in eq. S62. Singularity, which ensures that the function is increasingly
peaked at $\tau=0$ is also verified in a straight forward manner:

$\displaystyle\lim_{\tau\rightarrow 0}\eta_{N}(\tau)$ $\displaystyle=N-1$
(S64) $\displaystyle\lim_{N\rightarrow\infty}\frac{\eta_{N}(\tau\neq
0)}{\lim_{\tau\rightarrow 0}\eta_{N}(\tau)}$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{2}{N(N-1)\tau^{2}}-\frac{(1+i\tau)^{1-N}+(1-i\tau)^{1-N}}{N(N-1)\tau^{2}}=0$
(S65)

The last i.e. normalization requires some work. We need to compute the area
under the curve $\eta_{N}(\tau)$. We evaluate the integral using tricks used
to evaluate _Dirichlet integrals_ starting by splitting it up as follows

$\displaystyle\int_{-\infty}^{\infty}d\tau\eta_{N}(\tau)$
$\displaystyle=\mathcal{P}\left(\int_{-\infty}^{\infty}d\tau
I_{+}(N,\tau)\right)+\mathcal{P}\left(\int_{-\infty}^{\infty}d\tau
I_{-}(N,\tau)\right),$ (S66) $\displaystyle\text{ where, }I_{\pm}(N,\tau)$
$\displaystyle=\frac{1-(1+i\tau)^{1-N}}{N\tau^{2}}.$ (S67)

Note that while $\eta_{N}(\tau)$ is well behaved on the real line, $I_{\pm}$
are individually singular at $\tau=0$. $\mathcal{P}$ refers to the principal
component integration defined as follows

$\mathcal{P}\left(\int_{-\infty}^{\infty}d\tau
I_{\pm}(N,\tau)\right)=\lim_{\epsilon\rightarrow
0}\left(\int_{-\infty}^{-\epsilon}d\tau
I_{\pm}(N,\tau)+\int_{\epsilon}^{\infty}d\tau I_{\pm}(N,\tau)\right).$ (S68)

Let us proceed to evaluate each piece of the integral starting with $I_{+}$.
Notice that $I_{+}$ is singular at two points on the complex plane- $\tau=0$
and $\tau=i$.

---
Figure S12: Contour $C$ to evaluate the $I_{+}$ integral. The poles of $I_{+}$
are marked as $*$

Consider a contour that avoids both these points as shown in fig. S12.
$I_{+}(N,\tau)$ is analytic at all points inside the contour and hence, by
Cauchy’s theorem, we get

$\displaystyle\varointclockwise_{C}dzI_{+}(N,z)=0.$ (S69)

We can split the integral on the above contour as

$\displaystyle\int_{-\infty}^{-\epsilon}d\tau
I_{+}(N,\tau)+\int_{\epsilon}^{\infty}d\tau
I_{+}(N,\tau)+\int_{C_{\epsilon}}dzI_{+}(N,z)+\int_{C_{\infty}}dzI_{+}(N,z)=0$
(S70)

The contribution of the integral along the arc at infinity is $0$ since the
integrand $I_{+}$ vanishes at infinity. If we take the limit of the radius of
the semicircle encircling the origin, $\epsilon\rightarrow 0$, we have the
following

$\displaystyle\mathcal{P}\left(\int_{-\infty}^{\infty}d\tau
I_{+}(N,\tau)\right)=-\lim_{\epsilon\rightarrow
0}\int_{C_{\epsilon}}dzI_{+}(N,z).$ (S71)

The RHS of the above equation is evaluated easily. Along the semi-circular arc
$C_{\epsilon}$, the complex numbers can be parametrized as
$z=\epsilon~{}e^{i\theta}$. Thus we get

$\displaystyle\mathcal{P}\left(\int_{-\infty}^{\infty}d\tau
I_{+}(N,\tau)\right)$ $\displaystyle=-\lim_{\epsilon\rightarrow
0}\int_{-\pi}^{0}id\theta~{}(\epsilon~{}e^{i\theta})\frac{1-(1+i\epsilon~{}e^{i\theta})^{1-N}}{N(\epsilon~{}e^{i\theta})^{2}}$
$\displaystyle\approx-\lim_{\epsilon\rightarrow
0}\int_{-\pi}^{0}id\theta~{}(\epsilon~{}e^{i\theta})\frac{1-(1+i(1-N)\epsilon~{}e^{i\theta})}{N(\epsilon~{}e^{i\theta})^{2}}=\pi\left(1-\frac{1}{N}\right).$

The above steps can be repeated for $I_{-}$ which has poles at $0,-i$ and can
be evaluated by closing the contour on the upper half complex plane to get

$\displaystyle\mathcal{P}\left(\int_{-\infty}^{\infty}d\tau
I_{-}(N,\tau)\right)$ $\displaystyle=\pi\left(1-\frac{1}{N}\right).$ (S72)

Putting all these together, we get

$\displaystyle\int_{-\infty}^{\infty}d\tau\eta_{N}(\tau)$
$\displaystyle=2\pi\left(1-\frac{1}{N}\right),$ (S73)
$\displaystyle\lim_{N\rightarrow\infty}\int_{-\infty}^{\infty}d\tau\eta_{N}(\tau)$
$\displaystyle=2\pi.$ (S74)

Using these, for any well-behaved test function $f(\tau)$, we can conclude

$\displaystyle\int_{-\infty}^{\infty}d\tau f(\tau)\eta_{N}(\tau)\approx
f(0)\int_{-\infty}^{\infty}d\tau\eta_{N}(\tau)+\mathcal{O}\left(\frac{1}{N}\right)=2\pi
f(0)+\mathcal{O}\left(\frac{1}{N}\right).$ (S75)

Taking the limit of $N\rightarrow\infty$ above, we get eq. S63.

###  VI.3 S6c. Matching the $N\rightarrow\infty$ limit for the Random SFF

We now obtain the same $N\rightarrow\infty$ limit for the SFF of uniformly
distributed random numbers and show that eq. S36 reduces to eq. S58. Let us
start with eq S36 setting $\mu=1$ for convenience as before.

$\displaystyle K^{R}(\tau,N)$
$\displaystyle=N+N(N-1)\left|\frac{\sin(N\tau/2)}{(N\tau/2)}\right|^{2}.$
(S76)

We wish to recover the limit for integrable models eq. S58

$\lim_{N\rightarrow\infty}\frac{K^{R}(\tau,N)}{N}=1+2\pi\delta(\tau).$ (S77)

This amounts to proving that

$\displaystyle\lim_{N\rightarrow\infty}\zeta_{N}(\tau)$
$\displaystyle=2\pi\delta(\tau)$ (S78) $\displaystyle\text{ where,
}\zeta_{N}(\tau)$
$\displaystyle=(N-1)\left|\frac{\sin(N\tau/2)}{(N\tau/2)}\right|^{2}.$ (S79)

It is very clear from the plots of $\zeta_{N}(\tau)$ fig. S13 that the
functional form looks increasingly like the delta function with increasing
$N$.

---
Figure S13: $\zeta_{N}(\tau)$ for various $N$

Again, we need to show that the functional form of $\zeta_{N}$ approaches that
of the delta function when interpreted as an integration measure i.e. it
satisfies [5]:

  1. 1.

Symmetry: $\zeta_{N}(-\tau)=\zeta_{N}(\tau)$

  2. 2.

Singularity: $\lim_{N\rightarrow\infty}\frac{\zeta_{N}(\tau\neq
0)}{\lim_{\tau\rightarrow 0}\zeta_{N}(\tau)}=0$

  3. 3.

Normalization: $\lim_{N\rightarrow\infty}\int d\tau\zeta_{N}(\tau)=2\pi$

The first two are again easily verified. The third too is verified as follows.
First, let us change variables $\frac{\mu{\cal N}\tau}{2}=x$ to get

$\displaystyle\int_{-\infty}^{\infty}d\tau\zeta_{N}(\tau)$
$\displaystyle=(N-1)\int_{-\infty}^{\infty}d\tau\left(\frac{\sin(N\tau/2)}{(N\tau/2)}\right)^{2}=\frac{2(N-1)}{N}\int_{-\infty}^{\infty}d\tau\left(\frac{\sin
x}{x}\right)^{2}=\frac{2\pi(N-1)}{N}$ (S80)
$\displaystyle\implies\lim_{N\rightarrow\infty}\int d\tau\zeta_{N}(\tau)$
$\displaystyle=2\pi.$ (S81)

The integral $\int_{-\infty}^{\infty}d\tau\left(\frac{\sin
x}{x}\right)^{2}=\pi$ can be evaluated in the following way. Let us define

$\displaystyle I(b)$
$\displaystyle=\int_{-\infty}^{\infty}d\tau\left(\frac{\sin(bx)}{x}\right)^{2}$
(S82) $\displaystyle\frac{dI(b)}{db}$
$\displaystyle=\int_{-\infty}^{\infty}d\tau\left(\frac{\sin(2bx)}{x}\right)=\pi.$
(S83)

The result of $\frac{dI(b)}{db}$ is a standard integral that can be obtained
using the Dirichlet integral trick in the previous section. We now perform the
definite integral to get

$I(b)=\pi b+const.$ (S84)

Finally, the constant is evaluated to $0$ by observing $I(0)=0$ in eq. S82.
Finally, setting $b=1$, we get

$I(1)=\int_{-\infty}^{\infty}d\tau\left(\frac{\sin x}{x}\right)^{2}=\pi.$
(S85)

This completes the proof.

## References

  * Riser  _et al._ [2020] R. Riser, V. A. Osipov, and E. Kanzieper, Nonperturbative theory of power spectrum in complex systems, Annals of Physics 413, 168065 (2020).
  * Berry and Tabor [1977] M. V. Berry and M. Tabor, Level clustering in the regular spectrum, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 356, 375 (1977).
  * Marklof [2001] J. Marklof, The berry-tabor conjecture, in  _European Congress of Mathematics_ , edited by C. Casacuberta, R. M. Miró-Roig, J. Verdera, and S. Xambó-Descamps (Birkhäuser Basel, Basel, 2001) pp. 421–427.
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Balakrishnan [2003] V. Balakrishnan, All about the dirac delta function, Resonance 8, 48 (2003).


# Enhancing Medical Image Segmentation: Optimizing Cross-Entropy Weights and
Post-Processing with Autoencoders

Pranav Singh
Tandon School of Engineering
New York University
<EMAIL_ADDRESS>Luoyao Chen
Center for Data Science
New York University
<EMAIL_ADDRESS>Mei Chen
Center for Data Science
New York University
<EMAIL_ADDRESS>Jinqian Pan
Center for Data Science
New York University
<EMAIL_ADDRESS>Raviteja Chukkapalli
Courant Institute of Mathematical Sciences
New York University
<EMAIL_ADDRESS>Shravan Chaudhari
Courant Institute of Mathematical Sciences
New York University
<EMAIL_ADDRESS>Jacopo Cirrone
Center for Data Science
New York University
<EMAIL_ADDRESS>
###### Abstract

The task of medical image segmentation presents unique challenges,
necessitating both localized and holistic semantic understanding to accurately
delineate areas of interest, such as critical tissues or aberrant features.
This complexity is heightened in medical image segmentation due to the high
degree of inter-class similarities, intra-class variations, and possible image
obfuscation. The segmentation task further diversifies when considering the
study of histopathology slides for autoimmune diseases like dermatomyositis.
The analysis of cell inflammation and interaction in these cases has been less
studied due to constraints in data acquisition pipelines. Despite the
progressive strides in medical science, we lack a comprehensive collection of
autoimmune diseases. As autoimmune diseases globally escalate in prevalence
and exhibit associations with COVID-19, their study becomes increasingly
essential. While there is existing research that integrates artificial
intelligence in the analysis of various autoimmune diseases, the exploration
of dermatomyositis remains relatively underrepresented. In this paper, we
present a deep-learning approach tailored for Medical image segmentation. Our
proposed method outperforms the current state-of-the-art techniques by an
average of 12.26% for U-Net and 12.04% for U-Net++ across the ResNet family of
encoders on the dermatomyositis dataset. Furthermore, we probe the importance
of optimizing loss function weights and benchmark our methodology on three
challenging medical image segmentation tasks.

##  1 Introduction

The development of potent CAD (Computer Aided Diagnosis) strategies has been
aided by advances in computational power and image analysis algorithms over
the past decade. Medical imaging is fundamental to these CAD methods.
Obtaining accurate results from CAD techniques relies on acquiring high-
quality medical imaging and corresponding annotation. These CAD approaches
facilitate various tasks such as image classification, segmentation, spatial
mapping, and tracking. Out of these, medical image segmentation is a
particularly challenging task due to several complexities. For example, in
skin lesion image segmentation, there exists significant intra-class
variability and inter-class similarity. This issue is exacerbated by the
presence of obscuration and low contrast, which makes the task of separating
the affected area from the surrounding image more challenging. On the other
hand, sometimes the data required to segment is very complex, with multiple
fine-grained and hard-to-segment objects, for example, in the case of
histopathology data of dermatomyositis (a kind of autoimmune disease). We
provide a few examples of the large variability and low contrast in Figure 1,
obscuration in Figure 2, and multiple hard-to-segment small objects in Figure
3.

Figure 1: Samples from the dermofit dataset for segmentation, we observe a
large lesion color variability from left to right. Interestingly, the
background color remains relatively consistent throughout the samples. We also
observe a change in contrast from left to right. These confounding factors
make the segmentation of skin lesions difficult.

Figure 2: This figure contains samples from the ISIC 2017 dataset. The
ISIC-2017 dataset exhibits obscuration and significant intra-class
variability, resembling the Dermofit dataset (from Figure 1) regarding inter-
class similarity and low contrast (rightmost image).

Figure 3: Semantic segmentation task as defined in Section 3 with input image
on the left and the corresponding ground truth on the right. On top, we have a
sample image (on the left) and the corresponding ground truth (on the right)
from the dermatomyositis dataset. Similarly, a sample from the dermofit
dataset is in the middle, and a sample from the ISIC-2017 dataset is at the
bottom. Unlike a single blob in the skin lesion dataset samples from dermofit
and ISIC 2017, we observe that the histopathology whole slide image has many
more fine-grained objects with hard-to-segment boundaries. For all the images
on the right, the yellow area represents the region of interest (foreground),
and the rest is the area other than the region of interest (background).
Figure 4: DEDL Architecture[22] with the APP (Autoencoder Post Processing) for
dermatomyositis image segmentation as described in Section 2.1. We use this
architecture as our baseline and propose changes to it to improve by around
12.26% for U-Net and 12.04% for U-Net++ - this is a considerable improvement
over DEDL, as DEDL improved over the previous state-of-the-art approach Van
Buren _et al_. [26] by around 5% for segmentation.

In addition to these modality-specific complexities, medical imaging datasets
are considerably smaller than natural datasets. The main reason for this is
the significant expenses and time involved in gathering, annotating medical
datasets and privacy concerns. Medical imaging datasets can only be labeled by
highly specialized clinicians instead of the possibility of crowdsourced
labeling in the case of natural datasets. Privacy concerns pose significant
challenges in the open sourcing of medical datasets, particularly for rare or
emerging diseases. Medical datasets are typically restricted to institutional
use, even when made available [16].

Despite significant progress in medical science, some diseases have not yet
been fully comprehended [4]. Autoimmune diseases are a notable category in
this context. The lack of a comprehensive catalog of autoimmune diseases,
unlike other diseases, is attributed to the diverse nature of their onset and
progression [24]. There are still important research questions for autoimmune
diseases regarding environmental triggers, pathogenesis, cell inflammation,
and interaction. Currently, there are over 80 classified autoimmune diseases.
Immune-modulatory drugs are commonly employed for the treatment of autoimmune
diseases. However, these drugs have a wide range of effects and lack
specificity for autoimmune diseases. Unfortunately, their usage is often
linked to other infections and malignant diseases as undesirable side effects.
Patients often have limited or no response to these treatments due to the
variability within these disorders. So, there is a pressing need for more
advanced, fast, and accurate ways to find novel relationships and pathologies
that can lead to more effective treatments for autoimmune diseases. To
accomplish this, it is imperative to develop precise and adaptable techniques
for analyzing autoimmune diseases related medical images. Implementing AI-
based Computer-Aided Diagnosis (CAD) is a potential strategy for achieving
this objective. In contrast to other diseases, however, lacking a definitive
list and a limited understanding of autoimmune disorders presents a challenge.
Consequently, there are few established data collection mechanisms for
autoimmune diseases. These factors contribute to the paucity of research on
the intersection of autoimmune diseases and CAD approaches. Most extant
research in this field is either outdated or lacks open-source methodologies.
The study of autoimmune diseases is paramount due to their increasing
prevalence[5, 8, 11]. Autoimmune diseases impact a significant portion of the
global population, ranging from 5% to 8%. These conditions cause considerable
distress to patients and have been found to have connections with COVID-19,
the primary cause of the recent worldwide pandemic [21, 23]. To bridge the
divide, Van Buren _et al_. [26] and Singh & Cirrone [22] have made attempts.
The main focus of these studies is dermatomyositis. This rare autoimmune
disease has received limited attention at the intersection of medical imaging
and the application of AI (Artificial Intelligence) for medical image
analysis. With this paper,

  * •

We improve upon the existing state-of-the-art approach [22] for
dermatomyositis segmentation by an average of 12.26% for U-Net and 12.04% for
Unet++ in Section 5.1. Additionally, we benchmark our approach on two other
challenging skin-related datasets.

  * •

We study the impact of adding a post-processing autoencoder in addition to
U-Net and U-Net++ on three medical imaging datasets in Section 5.1.1.

  * •

We investigate the significance of cross-entropy loss function weights on
three challenging medical imaging datasets in Section 5.1.2.

Encoder | Technique | U-Net
---|---|---
Baseline(w/o APP) | w/ Relu APP | w/ Gelu APP
ResNet-18 | DEDL | 0.4347 | 0.4608 | 0.4788
Ours | 0.5618 | 0.5479 | 0.5582
ResNet-34 | DEDL | 0.4774 | 0.4467 | 0.4983
Ours | 0.5306 | 0.5571 | 0.5606
ResNet-50 | DEDL | 0.3798 | 0.4187 | 0.3827
Ours | 0.5556 | 0.5495 | 0.5597
ResNet-101 | DEDL | 0.3718 | 0.4074 | 0.4402
Ours | 0.5502 | 0.5678 | 0.5497
Table 1: Performance comparison of DEDL [22] and our approach on the dermatomyositis dataset for U-Net. We repeat all experiments five times with different seed values and report the IoU on the test in the 95% CI (confidence interval).  Encoder | Technique | U-Net++
---|---|---
Baseline(w/o APP) | w/ Relu APP | w/ Gelu APP
ResNet-18 | DEDL | 0.5274 | 0.4177 | 0.4707
Ours | 0.5622 | 0.5679 | 0.5683
ResNet-34 | DEDL | 0.3745 | 0.4535 | 0.4678
Ours | 0.5536 | 0.5685 | 0.5633
ResNet-50 | DEDL | 0.4236 | 0.4685 | 0.4422
Ours | 0.5742 | 0.5698 | 0.5514
ResNet-101 | DEDL | 0.4311 | 0.4265 | 0.4467
Ours | 0.57 | 0.5727 | 0.5692
Table 2: Similar to Table 1, in this table, we compare the performance
comparison of DEDL [22] and our approach on U-Net++. We report IoU scores
averaged over five seed values in the 95% confidence interval (CI) over the
test set of the Dermatomyositis dataset in this table.

##  2 Background

Dataset | ResNet | U-Net
---|---|---
Baseline(w/o APP) | w/ Relu APP | w/ Gelu APP
Dermofit | ResNet18 | 0.7388 | 0.7477 | 0.7467
ResNet34 | 0.7576 | 0.7633 | 0.7525
ResNet50 | 0.7364 | 0.7338 | 0.7401
ResNet101 | 0.7252 | 0.7213 | 0.7258
Dermatomyositis | ResNet18 | 0.5618 | 0.5479 | 0.5582
ResNet34 | 0.5306 | 0.5571 | 0.5606
ResNet50 | 0.5556 | 0.5495 | 0.5597
ResNet101 | 0.5502 | 0.5678 | 0.5497
ISIC2017 | ResNet18 | 0.6458 | 0.6252 | 0.6357
ResNet34 | 0.6518 | 0.6227 | 0.6306
ResNet50 | 0.605 | 0.5984 | 0.6207
ResNet101 | 0.6267 | 0.6325 | 0.5884
Table 3: In this table we present the IoU (Intersection over Union) on the test set averaged over five seed values (in 95% CI) for U-Net trained with our proposed technique as mentioned in Section 3. Dataset | ResNet | UNet++
---|---|---
Baseline(w/o APP) | w/ Relu APP | w/ Gelu APP
Dermofit | ResNet18 | 0.744 | 0.7408 | 0.7366
ResNet34 | 0.754 | 0.7553 | 0.7599
ResNet50 | 0.737 | 0.7408 | 0.7379
ResNet101 | 0.7232 | 0.7264 | 0.7229
Dermatomyositis | ResNet18 | 0.5622 | 0.5679 | 0.5683
ResNet34 | 0.5536 | 0.5685 | 0.5633
ResNet50 | 0.5742 | 0.5698 | 0.5514
ResNet101 | 0.57 | 0.5727 | 0.5692
ISIC2017 | ResNet18 | 0.6096 | 0.6232 | 0.6005
ResNet34 | 0.6583 | 0.6423 | 0.6548
ResNet50 | 0.6103 | 0.6355 | 0.619
ResNet101 | 0.6018 | 0.6164 | 0.6041
Table 4: Similar, to table 3 in this table we present IoU averaged over five seed values (in 95% CI) for U-Net++. Supervision Level | Dataset | Baseline(w/o APP) | w/ Relu APP | w/ GeLu APP
---|---|---|---|---
U-Net | Dermofit | 0.7395 | 0.7415 | 0.7413
Dermatomoyositis | 0.5496 | 0.5556 | 0.5551
ISIC2017 | 0.6323 | 0.6197 | 0.6189
U-Net++ | Dermofit | 0.7396 | 0.7408 | 0.7393
Dermatomoyositis | 0.565 | 0.5697 | 0.5630
ISIC2017 | 0.6200 | 0.62935 | 0.6196
Table 5: Mean IoU in the 95% confidence interval when averaged over the entire
ResNet family for U-Net and U-Net++ on the dermofit, dermatomyositis, and the
ISIC 2017 dataset. We observe that adding ReLU autoencoder as a post-
processing unit improves performance for U-Net++ and, in almost all cases, for
U-Net.

Medical image segmentation separates the region of interest, usually a lesion,
cells, or other anatomical region of interest, from the slide background.
Traditional segmentation processes use pixel-level classification to group
pixels into different categories; in the case of semantic segmentation, these
would be background and foreground. But with the maturity of Convolutional
Neural Networks(CNNs), Ronneberger _et al_. [20] introduced U-Net - an
autoencoder-based architecture for biomedical segmentation. The U-Net consists
of an encoder and a decoder architecture, where the encoder acts as a feature
extractor, and the decoder learns the mask by using the extracted features as
input. In addition, the decoder also incorporates the feature maps from the
encoder to improve scaling up the representation to the image mask; these
connections are called ”Skip-connections.” Following U-Net, a wealth of
architectures have spawned: U-Net++[27], DeepLab [1], DeepLabV3+ [2], and
Feature Pyramid networks (FPN) [12]. All of these architectures build on the
autoencoder architecture of U-Net with skip connections. To increase the
receptive field of these architectures, various techniques have been
introduced, such as dilated networks [18] and nesting architecture, as in the
case of U-Net++[27]. Despite these advancements and complex architectures,
U-Net remains the choice of architecture for medical image segmentation [17,
22, 26].

###  2.1 Application of Segmentation Techniques for Autoimmune diseases.

Stafford _et al_. [24] conducted a comprehensive survey to examine the
application of AI in the context of autoimmune diseases. They observed the
median size of autoimmune datasets is much smaller (99-540 samples per
dataset) as compared to datasets pertaining to other medical modalities. The
scarce available data poses a significant challenge in acquiring informative
priors for artificial intelligence-based CAD approaches on these datasets
resulting in sub-par performance. Furthermore, most methodologies for
analyzing these datasets are antiquated and lack open-source availability[22].
To overcome these shortcomings Van Buren _et al_. [26] proposed the use of
U-Net for segmentation of whole slide images of dermatomyositis histopathology
data and open-sourced their approach. Given the considerable size of the whole
slide images (WSI) at 1408 $\times$ 1876, a tiling approach was employed to
partition the WSI into smaller 256 $\times$ 256 images, with padding.

They also used a combination of Dice and Binary cross entropy loss to
attenuate the problem of pixel distribution imbalance between the area
surrounding the region of interest (background pixels) and the region of
interest (foreground pixels). Due to this imbalance, the segmentation
architecture tends to focus more on the area surrounding the region of
interest than the region of interest for segmentation if unattended. Singh &
Cirrone [22] further improved on this benchmark by using U-Net and introduced
an “Autoencoder Post Processing” (APP) technique. The APP consists of stacked
linear layers for the encoder and decoder. This makes the autoencoder much
simpler than the convolution and skip connection based U-Net and U-Net++
architectures. After obtaining the mask from a U-Net or U-Net++, it is passed
through the APP. Since, the autoencoder consists of only stacked linear layers
it creates a noised version of the segmentation output from U-Net and U-Net++.
A mean squared error is then calculated between the autoencoder’s output and
ground truth. During training, the model trained with the help of the MSE loss
(calculated between the autoencoder output and the ground truth) and the cross
entropy loss (calculated between the U-Net/U-Net++ output and the ground
truth). This helps the model learn a more diverse set of features.

The autoencoder is only used during the training process. Hence there is only
a marginal increase in training time while the inference time remains
constant. They studied ReLU and GELU as activation functions for the linear
layers of APP and found that ReLU activations work better than GELU. For their
choice of architecture, they used U-Net and U-Net++ (nested U-net) with
Squeeze and Excitation [10] in the decoder for channel-level attention. To
navigate the problem of pixel distribution imbalance between the area
surrounding the region of interest (background pixels) and the region of
interest (foreground pixels), they used pixel-distribution ratio weights in
the cross-entropy loss. Wherein the background pixels used the ratio of
background pixels to total pixels as weights, and similarly, the foreground
pixels used the ratio of foreground pixels to total pixels as weights. With
these changes, they were able to improve on state of the art on the
dermatomyositis segmentation task [26] by around 5%. They also suggested a
change in evaluation metric from pixel accuracy to IoU (Intersection over
Union) as pixel accuracy does not correctly represent the quality of the
learned mask as opposed to the ground truth mask. Our study builds upon the
foundation established by Singh & Cirrone [22] as a baseline. Our proposed
methodology demonstrates significant improvement in performance as compared to
the state-of-the-art approach [22]. We achieve an average improvement of
12.26% for U-Net and 12.04% for U-Net++, as elaborated in Section 5.1. Next,
we benchmark our methodology on two complex skin lesion datasets in Section
5.1.1. Furthermore, we investigate the impact of autoencoder for post-
processing in Section 5.1.1 and the significance of loss function weights in
Section 5.1.2.

##  3 Methodology

We start with Singh & Cirrone’s [22] approach on the dermatomyositis dataset.
They use U-Net and U-Net++ as the choice of segmentation architecture with
Squeeze and Excitation [10] in the decoder for channel-level attention.
Similar to previous studies, our work focuses on semantic segmentation, where
the goal is to categorize each pixel in an image into a class. For semantic
segmentation, these classes would be - a region of interest (for example,
cells in the case of the dermatomyositis dataset) and an area other than the
region of interest or background (region other than the cells). For the
encoder, we study the performance of the entire ResNet family of CNNs[9]. We
use an encoder depth of three, increasing the convolution filter size from
128, 256, to 512. We initialize the encoder with ImageNet pre-trained weights.
In the decoder part of U-Net and U-Net++, we use a convolution channel scheme
of 256, 128, and 64. For each decoder block, we also use batch normalization
as well as squeeze-and-excitation channel excitation after the convolutional
layer. As discussed in [22, 26] and Section 2.1, the Dermatomyositis whole
slide images contain a lot more pixels without cells (background) as opposed
to with cells (foreground). To attenuate this imbalance in the distribution of
pixels, we use Cross Distribution Weights (CDW) in the cross-entropy loss. Van
Buren _et al_. [26] used random weights. In contrast, Singh & Cirrone [22]
used a ratio of the pixel with cells to total pixels and a ratio of non-cell
pixels to total pixels as weights for foreground and background, respectively,
in the cross-entropy loss. We propose to swap the weights and instead use the
ratio of the number of pixels not containing the cell to the total number of
pixels as the weight for the foreground. Similarly, the weight of the
background is the ratio of the number of pixels containing cells
(foreground/object of interest) to the total number of pixels. This
alternative weight assignment method aims to enhance the foreground
representation. This intuition is very similar to focal loss [13], wherein the
misclassification of the minority class is penalized more than that of the
minority class. To ensure our results are statistically significant, we
conduct all experiments over five different seed values and report the mean
values in the 95% confidence interval (C.I) over the five runs. We present
these results in Tables 3 and 4. Based on our proposal, we observe an average
improvement of 12.26% for U-Net and 12.04% for U-Net++. We further discuss
these in Section 5.1. Additionally, we benchmark our approach and study the
impact of autoencoder post-processing on two additional challenging
dermatology-related datasets - ISIC 2017 and the dermofit dataset in Section
5.1.1. Both datasets are challenging due to large intra-class variations and
inter-class similarities as depicted in Figure 1 and obscuration in 2.
Finally, we study the impact of using mean-frequency weights and compare the
results with distribution-swapped weights for U-Net and U-Net++ over the
ResNet family of encoders and three datasets in Section 5.1.2.

Encoder | Technique | U-Net | U-Net++
---|---|---|---
Baseline* | w/ ReLU APP | w/ GELU APP | Baseline* | w/ ReLU APP | w/ GELU APP
ResNet-18 | CDW | 0.5618 | 0.5479 | 0.5582 | 0.5622 | 0.5679 | 0.5683
Mean Frequency | 0.5645 | 0.5592 | 0.5405 | 0.5852 | 0.5603 | 0.5814
ResNet-34 | CDW | 0.5306 | 0.5571 | 0.5606 | 0.5536 | 0.5685 | 0.5633
Mean Frequency | 0.5555 | 0.5551 | 0.5616 | 0.5729 | 0.5763 | 0.57
ResNet-50 | CDW | 0.5556 | 0.5495 | 0.5597 | 0.5742 | 0.5698 | 0.5514
Mean Frequency | 0.5512 | 0.5652 | 0.5585 | 0.57 | 0.5723 | 0.5929
ResNet-101 | CDW | 0.5502 | 0.5678 | 0.5497 | 0.57 | 0.5727 | 0.5692
Mean Frequency | 0.5506 | 0.5537 | 0.5596 | 0.5892 | 0.5678 | 0.5773
Table 6: This table displays the average Intersection over Union (IoU) values obtained from five separate runs on the dermatomyositis test set, with their corresponding 95% confidence intervals (CI). In this context, CDW refers to the utilization of cross-distribution weights in the calculation of cross-entropy loss. Here, Baseline* represents the use of the segmentation architecture with autoencoder for post-processing (APP). Encoder | Technique | U-Net | U-Net++
---|---|---|---
Baseline* | w/ ReLU APP | w/ GELU APP | Baseline* | w/ ReLU APP | w/ GELU APP
ResNet-18 | CDW | 0.7388 | 0.7477 | 0.7467 | 0.744 | 0.7408 | 0.7366
Mean Frequency | 0.75 | 0.7498 | 0.7377 | 0.7469 | 0.7413 | 0.7449
ResNet-34 | CDW | 0.7576 | 0.7633 | 0.7525 | 0.754 | 0.7553 | 0.7599
Mean Frequency | 0.7533 | 0.7633 | 0.7535 | 0.7602 | 0.7635 | 0.7547
ResNet-50 | CDW | 0.7364 | 0.7338 | 0.7401 | 0.737 | 0.7408 | 0.7379
Mean Frequency | 0.7379 | 0.731 | 0.7385 | 0.7362 | 0.7358 | 0.7411
ResNet-101 | CDW | 0.7252 | 0.7213 | 0.7258 | 0.7232 | 0.7264 | 0.7229
Mean Frequency | 0.7212 | 0.7242 | 0.7236 | 0.7156 | 0.7247 | 0.7234
Table 7: Similar to Table 6, in this table, we present the IoU averaged over five runs in 95% confidence interval on the dermofit test set for U-Net and U-Net++. Like Table 6, CDW represents the scenario when cross-distribution weights are used for the cross-entropy loss and Baseline* represents the use of the segmentation architecture with autoencoder for post-processing (APP). Encoder | Technique | U-Net | U-Net++
---|---|---|---
Baseline* | w/ ReLU APP | w/ GELU APP | Baseline* | w/ ReLU APP | w/ GELU APP
ResNet-18 | CDW | 0.6458 | 0.6252 | 0.6357 | 0.6096 | 0.6232 | 0.6005
Mean Frequency | 0.6257 | 0.6394 | 0.6307 | 0.6177 | 0.6198 | 0.6074
ResNet-34 | CDW | 0.6518 | 0.6227 | 0.6306 | 0.6583 | 0.6423 | 0.6548
Mean Frequency | 0.6314 | 0.6409 | 0.6322 | 0.6412 | 0.6454 | 0.6513
ResNet-50 | CDW | 0.605 | 0.5984 | 0.6207 | 0.6103 | 0.6355 | 0.619
Mean Frequency | 0.6396 | 0.6223 | 0.6337 | 0.6354 | 0.628 | 0.646
ResNet-101 | CDW | 0.6267 | 0.6325 | 0.5884 | 0.6018 | 0.6164 | 0.6041
Mean Frequency | 0.6137 | 0.6283 | 0.6175 | 0.6049 | 0.6112 | 0.6016
Table 8: In this table, we showcase the IoU results for U-Net and U-Net++ on the ISIC 2017 test set. The values represent the average IoU over five runs in the 95% confidence interval. CDW, in this context, refers to the utilization of cross-distribution weights for cross-entropy loss, as demonstrated in Tables 6 and 7. Baseline* represents the use of the segmentation architecture with autoencoder for post-processing. Dataset | CDW | Median Frequency
---|---|---
Dermofit | 0.7395 | 0.7406
DM* | 0.5496 | 0.5564
ISIC2017 | 0.6197 | 0.6276
Table 9: In this table, we present the average IoU (in the 95% confidence interval) over the ResNet family of encoders and the three paradigms (Baseline, with ReLU and GELU APP) for U-Net over the dermatomyositis, the dermofit, and ISIC 2017 test sets from Tables 6, 7 and 8 respectively. Here, DM* represents the Dermatomyositis dataset, and CDW represents Cross Distribution Weight. Dataset | CDW | Median Frequency
---|---|---
Dermofit | 0.7396 | 0.7397
DM* | 0.565 | 0.5793
ISIC2017 | 0.62 | 0.6248
Table 10: Similar to Table 9, In this table, we provide the average Intersection over Union (IoU) values (in 95% CI) for the ResNet family for U-Net++ architecture on three different test sets: Dermatomyositis, Dermofit, and ISIC 2017. Dataset | CDW | Median Frequency
---|---|---
Dermofit | [0.3037, 0.6963] | [0.7180, 1.6466]
DM* | [0.1479, 0.8521] | [0.5986, 3.0348]
ISIC2017 | [0.2020, 0.7980] | [0.6265, 2.4755]
Table 11: Cross-entropy weights used for experimentation as described in
Section 5.1.1. Here, DM* represents the Dermatomyositis dataset as described
in Section 4.1 and CDW represents Cross Distribution Weight.

##  4 Experimental Details

###  4.1 Datasets

We use a 70-10-20 split for the dermatomyositis and dermofit datasets. For the
ISIC-2017 dataset, we use the same splits as used in the 2017 ISIC
competition.

##### Dermatomyositis:

We use the same dataset as used in previous works on dermatomyositis
segmentation [22, 26]. To give an idea about the modality of the dataset, we
show a random sample from the test set in Figure 3. The Dermatomyositis
dataset is collected from 198 muscle biopsies collected from seven
dermatomyositis patients. These files are then stored in TIFF format. Each
TIFF image contains eight slides that indicate the presence or absence of
phenotypic markers by setting binary thresholds for each channel (1-DAPI,
2-CXCR3, 3-CD19, 4-CXCR5, 5-PD1, 6-CD4, 7-CD27, 8-Autofluorescence). For
segmentation, we used the DAPI-stained image. Each whole slide image was tiled
into 480x480. We further expand on this in Section 4.1.2. This is a
particularly challenging dataset due to the large number of fine-grained
objects (cells) to be segmented per image, as discussed in Section 1.

##### Dermofit [7]:

As shown in Figure 3, the Dermofit dataset contains 1300 skin lesion RGB
images. These data are taken with a high-quality SLR camera in controlled
(ring flash) indoor illumination. The Dermofit dataset contains ten
categories; each includes a different number of instances: Actinic Keratosis
(AK): 45, Basal Cell Carcinoma (BCC): 239, Melanocytic Nevus / Mole (ML): 331,
Squamous Cell Carcinoma (SCC) sample 88, Seborrhoeic Keratosis (SK): 257,
Intraepithelial carcinoma (IEC): 78, Pyogenic Granuloma (PYO): 24, Haemangioma
(VASC): 96, Dermatofibroma (DF): 65, Melanoma (MEL): 76. No two images in this
dataset are of the same size, as a preprocessing step we interpolate all
images to 480x480 and then resize to 224x224 to ensure uniformity.

##### ISIC Challenge 2017 Dataset, Lesion Segmentation Task [3]:

The International Skin Imaging Collaboration (ISIC) is a large publicly
accessible dataset. We show a sample from the test set in Figure 3. In our
case, we use the segmentation dataset from 2017 and use the original splits
wherein 2,000 images were used as training, 150 images as validation, and 600
images as the test set. The ISIC 2017 and the Dermofit datasets described
above are skin lesion datasets with high intra-class variability and inter-
class similarities with obscuration areas of interest, as discussed in Section
1.

####  4.1.1 Common Implementation Details

We implemented all models in Pytorch [19] using a single NVIDIA RTX-8000 GPU
with 64 GB RAM and 3 CPU cores. All models are trained with an Adam optimizer
with an initial learning rate (lr) of 3.6e-4 and a weight decay 1e-5. We use a
cosine annealing scheduler with a maximum of 50 iterations and a minimum
learning rate of 3.4e-4 to adjust the learning rate based on each epoch. We
train all architectures for 50 epochs with batch size 16, followed by testing
on a held-out set. We use IoU (Intersection over Union) as our evaluation
metric on the test set. This aligns with previous work by Singh & Cirrone
[22]. We repeat all experiments with five different seed values and report the
mean value in the 95% confidence interval in all tables.

####  4.1.2 Data-Preprocessing

Images of the dermatomyositis dataset have a uniform size of 1408 $\times$
1876; we tiled each image into 12 sub-images of size 480 $\times$ 480 inline
with previous work [22]. In contrast, the Dermofit and the ISIC2017 datasets
contain images of different sizes, i.e., no two images in the dataset are the
same size. Additionally, since the other two datasets (dermofit and ISIC-2017)
contain skin lesions, they have significantly denser and larger mask labels
than the dermatomyositis dataset. Thus, a different image preprocessing step
is applied to the latter two datasets: bilinear interpolation to 480 $\times$
480 followed by a resize to 224 $\times$224\. For augmentation, we use the
same set of augmentation as used in Singh & Cirrone’s work [22], along with
Red channel normalization or ”Rnorm” [25] for all of our experiments.

##  5 Results and Discussion

###  5.1 Improvement over the current state-of-the-art for Dermatomyositis WSI
Segmentation [22]

Following the methodology (Section 3) and experimentation setup (Section 4),
we present the IoU averaged over five runs in the 95% confidence interval on
the test set for the Dermatomyositis dataset in Table 1 for U-Net and in Table
2 for U-Net++. We observe that our approach improves over Singh & Cirrone’s
approach (DEDL) [22] consistently over the entire ResNet family for baseline
as well as with APP (both ReLU and GELU based) for both U-Net and U-Net++.
When averaged over the ResNet family of encoders and the three paradigms
(baseline approach without using autoencoders, ReLU autoencoders, and GELU-
based autoencoders), we observed that our approach improves over the previous
state-of-art [22] for Dermatomyositis segmentation by 12.26% and 12.04% for
U-Net and U-Net++ respectively.

####  5.1.1 Impact of Incorporating Autoencoder Post-Processing.

As described in Section 2.1, Singh and Cirrone [22] introduced an ”Autoencoder
Post Processing” unit or APP after the main segmentation architecture. The
purpose of this autoencoder was to provide a more noised version of the
prediction from the U-Net or U-Net++. The mean square error loss between the
noised output and the ground truth mask, along with the weighted-cross entropy
loss between the output of the U-Net or U-Net++ and the ground truth, is
optimized during training. This is depicted in Figure 4. They studied the
impact of using APP with ReLU and GELU activations only on the Dermatomyositis
dataset. In this section, with our improved approach as presented in Section
3, we study the impact of adding APP on two additional challenging dermatology
datasets. We present the IoU over the test set in the 95% confidence interval
over the ResNet family of encoders in Tables 3 and 4 for ISIC 2017 and the
Dermofit dataset, respectively. Adding ReLu and GeLU-based APP improves
performance over the baseline architecture (with no APP) for U-Net and U-Net++
in most cases for the Dermofit. To better understand the result, we average
the IoU on the test set over the entire ResNet family for U-Net and U-Net++
and present the results in Table 5. From Table 5, we observe that the addition
of APP, especially ReLU-based APP, does improve performance over the baseline
(not using APP) in almost all cases for U-Net++ and U-Net. The addition of APP
is did not improve performance only in the case of the ISIC-2017 dataset for
U-Net.

####  5.1.2 Impact of Cross-entropy loss weights

In section 3, we explained our rationale for switching from distribution-based
weights to cross-distribution-based weights for the cross-entropy loss. In
this section, we study the impact of changing the cross-entropy weights from
cross-distribution-based weights to mean frequency weights[6]. The median
frequency weight received by each class is derived from the reciprocal of the
pixel ratio of a particular class, normalized by the median frequency of the
class[14, 15]. The median frequency and the cross-distribution weights,
calculated over our three datasets, are mentioned in Table 11. Mathematically,
median frequency weights ($w_{c}$) are defined as follows:
$w_{c}=\frac{\text{med\\_freq}}{n_{c}}$. Here, $n_{c}$ is the number of pixels
belonging to class $c$ in the training dataset, and med_freq is the median of
the frequency of pixels belonging to each class in the dataset.111In our case,
there are only two classes - foreground (region of interest) and background
(area other than the region of interest). Where $n_{c}=\frac{\alpha}{\beta}$,
here, $\alpha$ represents the number of pixels of a class, and $\beta$
represents the total number of pixels in images where the given class is
present. We compare the weights calculated by cross-distribution and median
frequency in Table 11. We provide the full comparative result of using cross-
distribution weights and mean frequency over the three datasets in the 95%
confidence interval averaged over five seed values Tables 6, 7 and 8 for the
Dermatomyositis, Dermofit, and the ISIC-2017 datasets, respectively.
Additionally, to summarize these results, we present the average over the
ResNet family and the three training paradigms (baseline without APP and APP
with ReLU and GELU layers) in Tables 9 and 10 for U-Net and U-Net++,
respectively. From these tables, we observe that median-frequency weights for
cross-entropy loss improve performance over cross-distribution weights,
although the improvement is marginal in almost all cases.

##  6 Conclusion

We observed that our approach of using Cross Distribution Weights (CDW)
improved segmentation performance over the previous state-of-the-art approach
for dermatomyositis segmentation [22] by 12.26% for U-Net and by 12.04% for
U-Net++ averaged over the ResNet family. Furthermore, adding APP (Autoencoder
Post Processing) improves segmentation performance marginally in the case of
dermatomyositis and dermofit datasets. In the case of the ISIC 2017 dataset,
the addition of APP is only useful in the case of U-Net++. We have open-
sourced our approach at https://github.com/pranavsinghps1/Enhancing-Medical-
Image-Segmentation. We hope that our study and open-sourced approach will
catalyze further research at the intersection of autoimmune diseases like
dermatomyositis and the application of AI as well as for other dermatology-
related datasets. This would help us better understand the immunology of
autoimmune diseases and answer some of the critical research questions to
develop improved healthcare solutions.222Potential negative societal impact:
Autoimmune diseases are extremely heterogeneous; the dermatomyositis dataset
used in our experiments is geographically restricted. Hence, this is a study
of a particular variant. This study might or not be generalizable for other
variants. Hence application on a wider scale for real-life scenarios should
only be trusted after clearance from the concerned health and safety governing
bodies.

##### Acknowledgements

We would like to thank NYU HPC team for assisting us with our computational
needs.

## References

  * [1] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  IEEE transactions on pattern analysis and machine intelligence, 40(4):834–848, 2017.
  * [2] Liang-Chieh Chen, Yukun Zhu, George Papandreou, Florian Schroff, and Hartwig Adam.  Encoder-decoder with atrous separable convolution for semantic image segmentation.  In Proceedings of the European conference on computer vision (ECCV), pages 801–818, 2018.
  * [3] Noel CF Codella, David Gutman, M Emre Celebi, Brian Helba, Michael A Marchetti, Stephen W Dusza, Aadi Kalloo, Konstantinos Liopyris, Nabin Mishra, Harald Kittler, et al.  Skin lesion analysis toward melanoma detection: A challenge at the 2017 international symposium on biomedical imaging (isbi), hosted by the international skin imaging collaboration (isic).  In 2018 IEEE 15th international symposium on biomedical imaging (ISBI 2018), pages 168–172. IEEE, 2018.
  * [4] Glinda S Cooper, Milele LK Bynum, and Emily C Somers.  Recent insights in the epidemiology of autoimmune diseases: improved prevalence estimates and understanding of clustering of diseases.  Journal of autoimmunity, 33(3-4):197–207, 2009.
  * [5] Michael Ehrenfeld, Angela Tincani, Laura Andreoli, Marco Cattalini, Assaf Greenbaum, Darja Kanduc, Jaume Alijotas-Reig, Vsevolod Zinserling, Natalia Semenova, Howard Amital, et al.  Covid-19 and autoimmunity.  Autoimmunity reviews, 19(8):102597, 2020.
  * [6] David Eigen and Rob Fergus.  Predicting depth, surface normals and semantic labels with a common multi-scale convolutional architecture.  CoRR, abs/1411.4734, 2014.
  * [7] Robert Fisher and Jonathan Rees.  Dermofit project datasets.  2017\.
  * [8] Caroline Galeotti and Jagadeesh Bayry.  Autoimmune and inflammatory diseases following covid-19.  Nature Reviews Rheumatology, 16(8):413–414, 2020.
  * [9] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition. corr abs/1512.03385 (2015), 2015.
  * [10] Jie Hu, Li Shen, and Gang Sun.  Squeeze-and-excitation networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7132–7141, 2018.
  * [11] Aaron Lerner, Patricia Jeremias, and Torsten Matthias.  The world incidence and prevalence of autoimmune diseases is increasing.  Int J Celiac Dis, 3(4):151–5, 2015.
  * [12] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2117–2125, 2017.
  * [13] Tsung-Yi Lin, Priya Goyal, Ross Girshick, Kaiming He, and Piotr Dollár.  Focal loss for dense object detection.  In Proceedings of the IEEE international conference on computer vision, pages 2980–2988, 2017.
  * [14] Yahui Liu, Jian Yao, Xiaohu Lu, Menghan Xia, Xingbo Wang, and Yuan Liu.  Roadnet: Learning to comprehensively analyze road networks in complex urban scenes from high-resolution remotely sensed images.  IEEE Transactions on Geoscience and Remote Sensing, 57(4):2043–2056, 2019.
  * [15] Yahui Liu, Jian Yao, Xiaohu Lu, Renping Xie, and Li Li.  Deepcrack: A deep hierarchical feature learning architecture for crack segmentation.  Neurocomputing, 338:139–153, 2019.
  * [16] Christos Matsoukas, Johan Fredin Haslum, Moein Sorkhei, Magnus Söderberg, and Kevin Smith.  What makes transfer learning work for medical images: feature reuse & other factors.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 9225–9234, 2022.
  * [17] Sachin Mehta, Ezgi Mercan, Jamen Bartlett, Donald Weaver, Joann Elmore, and Linda Shapiro.  Y-net: Joint segmentation and classification for diagnosis of breast biopsy images.  In International Conference on Medical image computing and computer-assisted intervention. Springer, 2018.
  * [18] Sachin Mehta, Mohammad Rastegari, Anat Caspi, Linda Shapiro, and Hannaneh Hajishirzi.  Espnet: Efficient spatial pyramid of dilated convolutions for semantic segmentation.  In Proceedings of the european conference on computer vision (ECCV), pages 552–568, 2018.
  * [19] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [20] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [21] Chetan Sharma and Jagadeesh Bayry.  High risk of autoimmune diseases after covid-19.  Nature Reviews Rheumatology, pages 1–2, 2023.
  * [22] Pranav Singh and Jacopo Cirrone.  A data-efficient deep learning framework for segmentation and classification of histopathology images.  In European Conference on Computer Vision, pages 385–405. Springer, 2022.
  * [23] Anuroop Sriram, Matthew Muckley, Koustuv Sinha, Farah Shamout, Joelle Pineau, Krzysztof J Geras, Lea Azour, Yindalon Aphinyanaphongs, Nafissa Yakubova, and William Moore.  Covid-19 prognosis via self-supervised representation learning and multi-image prediction.  arXiv preprint arXiv:2101.04909, 2021.
  * [24] IS Stafford, M Kellermann, E Mossotto, Robert Mark Beattie, Ben D MacArthur, and Sarah Ennis.  A systematic review of the applications of artificial intelligence and machine learning in autoimmune diseases.  NPJ digital medicine, 3(1):30, 2020.
  * [25] Dang N.H. Thanh, Uğur Erkan, V.B. Surya Prasath, Vivek Kumar, and Nguyen Ngoc Hien.  A skin lesion segmentation method for dermoscopic images based on adaptive thresholding with normalization of color models.  In 2019 6th International Conference on Electrical and Electronics Engineering (ICEEE), pages 116–120, 2019.
  * [26] Kayla Van Buren, Yi Li, Fanghao Zhong, Yuan Ding, Amrutesh Puranik, Cynthia A Loomis, Narges Razavian, and Timothy B Niewold.  Artificial intelligence and deep learning to map immune cell types in inflamed human tissue.  Journal of Immunological Methods, 505:113233, 2022.
  * [27] Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang.  Unet++: A nested u-net architecture for medical image segmentation.  In Deep learning in medical image analysis and multimodal learning for clinical decision support, pages 3–11. Springer, 2018.