# Supplemental Material. Excitons and carriers in transient absorption and
time–resolved ARPES spectroscopy: an abinitio approach

D. Sangalli  Istituto di Struttura della Materia of the National Research
Council, Via Salaria Km 29.3, I-00016 Montelibretti, Italy

##  I Convergence for the absorption spectrum and the excitonic wave–functions

###  I.1 Computational details for equilibrium properties and band structure

The ground state properties of LiF are computed with an energy cutoff on the
wave–functions of 80 Ry and a k-points grid $6\times 6\times 6$. Then a non-
self consistent calculation is performed on three k-points grids, i.e
16x16x16, 24x24x24, and 32x32x32 with 100 bands.

Figure 1: (left panel) Band structure in the energy range $-21$ \- $21$ eV
which contains the first 5 bands of LiF. (central panel) zoom in the valence
band region, bands from 2 to 4. (right panel) zoom in the conduction band
region, band 5. Points correspond to calculated eigenvalues on the 16x16x16
grid, while lines are obtained via 3D Fourier interpolation.

In Fig. 1 the band structure obtained from a Fourier interpolation starting
from the 16x16x16 is shown. The first occupied valence band lies at
$\approx-20$ eV, then there are three which are degenerate at $\Gamma$ with
$\epsilon_{v\Gamma}=0.$ eV and span an energy range of $\approx 3$ eV. The
first conduction band determines the gap with $\epsilon_{c1\Gamma}=14.5$ eV
spans an energy range of $\approx 6$ eV. A shissor operator of 6.05 eV is
applied in the plot to shift the KS band gap. The second conduction band start
at $\epsilon_{c2\Gamma}=20.0$ eV.

###  I.2 Absorption spectrum

For the calculation of the BSE matrix, bands from 2 to 5 are then considered
(for the 16x16x16 we also verified that using bands from 1 to 6 no significant
changes are observed). Indeed the critical parameter for the convergence of
the spectrum is the sampling of the Brillouin zone, as discussed more in
details in the present section. The other parameters used for the BSE
calculations are: 965 Reciprocal lattice vectors for the exchange, 59
reciprocal lattice vectors and 100 bands for the screening, 59 reciprocal
lattice vectors for the screened electron–hole interaction. The BSE is then
solved using the Haydock solver to obtain the spectrum and the SLEPC library
to obtaine the first 500 poles of the BSE matrix.

Here I report the convergence of the equilibrium absorption spectra and
amplitudes considering a ${\mathbf{k}}$–grid up to $32\times 32\times 32$.
Solving the BSE on such grid is demanding but feasible (there are $32^{3}\cdot
3=98.304$ transition, and the BSE has 9.663.676.416 matrix elements). We are
not able to perform real–time simulations on the same grid, mostly because the
real-time code is not as much optimized as the BSE code.

Figure 2: Convergence of the absorption spectrum with k-grids 16x16x16
(orange), 24x24x24 (red) and 32x32x32 (blue). For each simulation the first
500 eigenvalues of the excitonic hamiltonian are reported both weighted by the
oscillator strengths (under the spectrum, positive $y$-axis) and as a list of
bars also including dark transitions ()negative $y$-axis).

The spectrum is fully converged with the $24\times 24\times 24$ grid. Already
the $16\times 16\times 16$ grid gives a very good result, with the energy of
the first two bright excitonic peaks within 50 meV from the denser grid and
few oscillations in the continuum. In the continuum convergence is achieved
using a smearing paramenter $\eta=0.1$ eV. This means that real-time
simulations would be converged in such energy range using the same smearing
parameter. In the limit of very small smearing parameter we need to look at
the BSE eigenvalues. At the BSE level these are reasonably converged with the
$24\times 24\times 24$ grid in the bound region of the spectrum, but still not
converged in the continuum. For example the region around $E=14.8$ eV which we
consider in the main text has a very dense spectrum of poles. For the
$24\times 24\times 24$ we obtain 79 poles in a region of 50 meV (14.775
eV$<E<$14.824 eV), i.e. an average of 3 poles each 2 meV. However only very
few poles have non negligible intensity in the absorption spectrum (see inset
of fig. 2). In the limit of infinite k-points sampling we expect the number of
poles to become infinite in order to capture the absorption in the continuum.
However it is not possible to reach such situation numerically.

###  I.3 Excitonic wave-function

v

Figure 3: (left panel) Representation of the excitonic amplitudes
$F^{\lambda}_{\eta}(\omega)$ defined in eq. (1) for four different states and
$\eta=0.1\ eV$, with the k-grid 24x24x24. The matching energy $E_{\lambda}$ is
represented as a vertical bar with the same color for each state. The QP gap
is represented by a black dashed line. (right panel) the excitonic amplitude
just below the QP gap is compared for two different k-grids: 24x24x24 (blue
continuous) and 32x32x32 (blue dashed). The result slowly tends to a delta
like function.

We now move to the analysis of the excitonic wave–functions
$A^{\lambda}_{cv{\mathbf{k}}}$. To this end I introduce the excitonic
amplitudes which are a ${\mathbf{k}}$-integrated representation of the
excitonic wave–function:

$F^{eh(\lambda)}_{\eta}(\omega)=\sum_{cv{\mathbf{k}}}\frac{|A^{\lambda}_{cv{\mathbf{k}}}|^{2}}{\omega-\Delta\epsilon_{cv{\mathbf{k}}}+i\eta}.$
(1)

In Fig. 3 we see that for the lowest energy exciton, $\lambda=1s$,
$F^{eh(1s)}_{\eta}(\omega)$ has a very broad distribution and involves
electron–hole pairs up to more than $18\ eV$, i.e. $6\ eV$ above the
$\omega_{1s}$. This is due to the fact that the excitonic wave–function is
quite localized in space and requires many elements in the
$\psi_{c{\mathbf{k}}}\psi^{*}_{v{\mathbf{k}}}$ basis set to be represented.
However there is no relation between the poles of $F^{eh(1s)}_{\eta}(\omega)$
and the peak position $E_{1s}$ (which is marked by a vertical line in the
figure). The comparison between the pole $E_{1s}$ and the pole of
$F^{eh(1s)}(\omega)$ is quite similar to the comparison of the (AR)PES signal
of the coherent exciton state and the dephased state discussed in the main
text. Indeed the sole difference is how the energy of the valence band enters:
$\Delta\epsilon_{cv{\mathbf{k}}}$ and $\omega_{\lambda}$ vs
$\epsilon_{c{\mathbf{k}}}$ with $\omega_{\lambda}+\epsilon_{v{\mathbf{k}}}$.

It is interesting to inspect higher energy states. Above the QP gap,
$F^{\lambda}_{\eta}(\omega)$ are quite localized in energy. This means that
such $A^{\lambda}_{cv{\mathbf{k}}}$ states are more “delta like” and
correspond to specific electron–hole pairs. In Fig. 3 the state with index
121, at energy slightly below the QP gap, and the state with index 500, at
energy above the QP gap are shown. For such states the maxima of
$F^{\lambda}_{\eta}(\omega)$ move in the same way as the energies
$E_{\lambda}$. Compared to the maxima of $F^{\lambda}_{\eta}(\omega)$,
$E_{\lambda}$ is just shifted at lower energy by a constant. If the excitonic
state is delta–like, the energy shift would be roughly due to the diagonal
terms of the eh–interaction entering the BSE for such specific transition:
$W_{cv,cv}=\int
d^{3}{\mathbf{k}}W_{cv{\mathbf{k}},cv{\mathbf{k}}}A^{\lambda}_{cv{\mathbf{k}}}$.
In the limit of
$A^{\lambda}_{cv{\mathbf{k}}}\rightarrow\delta({\mathbf{k}}-{\mathbf{k}}_{0})$
such energy shift should go to zero since
$W_{cv{\mathbf{k}},cv{\mathbf{k}}}\approx 1/N_{{\mathbf{k}}}$. However
$A^{\lambda}_{cv{\mathbf{k}}}\neq\delta({\mathbf{k}}-{\mathbf{k}}_{0})$ since
it must be ortogonal to $A^{1s}_{cv{\mathbf{k}}}$. A converged quantitive
value of such energy shift as a function of the k-grid would however be very
difficult to achieve and it is beyond the scope of the present work. The goal
of the present discussion is to underline why in practice it would be very
difficult to converge the TD-HSEX simulations with the laser pulse tuned in
the continuum. As a consequence it is better to run TD-IP simulation in such
case. In a following section, we further analyze such hypothesis, comparing
the TR-ARPES spectrum of an HSEX simulations with the laser pulse in the
continuum, with the IP result.

##  II Non-equilibrium states generated by the pump pulse

###  II.1 Computational details for the real-time simulations

The same parameters used for BSE (i.e. energy cutoff and number of bands) are
also used for the real time propagation, with a time step of 5 as.
$\epsilon_{thresh}=10^{-5}$ Hartree is chosen as the numerical tolerance for
considering two energy levels $n$ and $m$ degenerate, so that $\rho_{nm}$ is
not subject to dephasing (see main text). The same energy threshold is used
for dipoles (see sec. IV.1.1 in the present SI). The real-time propagation has
roughly the same computational cost of a BSE calculations. In particular the
variation of the self energy can be constructed writing
$\Delta\Sigma^{HSEX}=K^{HSEX}\,\Delta G^{<}$, with $K^{HSEX}$ the BSE kernel.
There is, however, an advantage here. For BSE, $K^{HSEX}$ needs to be
constructed in the full BZ, for real time propagation it needs to be
constructed only for an ‘extended IBZ’, i.e. the IBZ defined by the symmetry
operations which are not broken by the presence of the laser field. The real
time propagation can be seen as an iterative solver of the BSE. It is less
demanding than the exact diagionalization, and in a way more similar to
perturbative schemes. The drawback is that the number of iterations (i.e. time
steps) may become very high. This depends on the goal of the propagation. To
obtain the results described in the present manuscript only linear absorption
properties are required and the cost is not much more demanding compared to
standard BSE iterative solvers like Haydock. For the TR-ARPES spectral
function however it is required to construct $G^{<}(t,t’)$ within GKBA. The
implementation I did is based on storing to disk the density matrix history,
which requires a lot of disk space. Alternative implementations which avoid
this are however possible.

###  II.2 Perturbative analysis of the EOM for the one-body density matrix

In view of the analysis of the state generated by the laser pulse we expand
the EOM for the density matrix $\rho$ up to second order in the pump laser
pulse ${\mathbf{E}}(t)$.

$\rho_{nm{\mathbf{k}}}(t)=\delta_{nm}f^{eq}_{n{\mathbf{k}}}+\rho_{nm{\mathbf{k}}}^{(1)}(t)+\rho_{nm{\mathbf{k}}}^{(2)}(t)+\cdots\,.$
(2)

For times preceding the activation of the pump
${\rho_{nm{\mathbf{k}}}=\delta_{nm}f^{eq}_{n{\mathbf{k}}}}$, i.e. $\rho$ is
defined by the equilibrium occupations of states. Later, when discussing the
transient absorption signal, we will consider the first order response in the
probe starting from the state generated by the pump. We split the analysis in
the independent particles (IP) part and in Hartree plus Screened EXchange part
(HSEX). The goal of this expansion is to highlight the role of the intraband
dipoles ${\mathbf{d}}_{nn{\mathbf{k}}}$, that are difficult to define and
which we neglect in our numerical implementation.

####  II.2.1 IP approximation

I express the EOM at the IP level in the following form

$i\partial_{t}\rho_{nm{\mathbf{k}}}-(\Delta\epsilon_{nm{\mathbf{k}}}-i\eta)\rho_{nm{\mathbf{k}}}=F_{nm{\mathbf{k}}}[\rho]\,,$
(3)

The source term of Eq. (3) reads

$F_{nm{\mathbf{k}}}[\rho]={\mathbf{E}}(t)\cdot\sum_{p}\left({\mathbf{d}}_{np{\mathbf{k}}}\rho_{pm{\mathbf{k}}}-\rho_{np{\mathbf{k}}}{\mathbf{d}}_{pm{\mathbf{k}}}\right)\,,$
(4)

where the ${\mathbf{d}}_{nm{\mathbf{k}}}$ are the matrix elements of the
dipole operator in the KS basis. Expanding order by order equations (3) and
(4) provides a chain of differential equations for the
$\rho^{(i)}_{nm{\mathbf{k}}}$:

$i\partial_{t}\rho^{(i)}_{nm{\mathbf{k}}}-(\Delta\epsilon_{nm{\mathbf{k}}}-i\eta)\rho^{(i)}_{nm{\mathbf{k}}}=F^{(i)}_{nm{\mathbf{k}}}\,,$
(5)

which can be hierarchically solved starting from the lowest order. At first
order

$F^{(1)}_{nm{\mathbf{k}}}=(f^{eq}_{n{\mathbf{k}}}-f^{eq}_{m{\mathbf{k}}}){\mathbf{E}}(t)\cdot{\mathbf{d}}_{nm{\mathbf{k}}}\,,$
(6)

is non vanishing only for transitions between valence and conduction states.
The solution of equation (5) at first order shows that only the matrix
elements $\rho^{(1)}_{cv{\mathbf{k}}}$ with
$\Delta\epsilon_{cv{\mathbf{k}}}\sim\omega_{pump}$ are activated by the pump
through a resonance mechanism. After an initial transient regime, these terms
exhibit an oscillatory behavior eventually suppressed by an exponential
factor, _i.e._
$\exp(-i(\Delta\epsilon_{cv{\mathbf{k}}}-i\eta_{cv{\mathbf{k}}})t)$ if
decoherence is accounted for. Going ahead in the hierarchy of equations (5),
the ${\rho^{(2)}_{nm{\mathbf{k}}}}$ terms are activated. Using the fact that
$\rho^{(1)}_{cc^{\prime}{\mathbf{k}}}=\rho^{(1)}_{vv^{\prime}{\mathbf{k}}}=0$
we can simplify the second order term. In the $cc^{\prime}$ and $vv^{\prime}$
channel the source term associated to transitions among conduction states
reads

$\displaystyle
F^{(2)}_{cc^{\prime}{\mathbf{k}}}={\mathbf{E}}(t)\cdot\sum_{v}\left[\mathbf{d}_{cv{\mathbf{k}}}\rho_{vc^{\prime}{\mathbf{k}}}^{(1)}-\rho_{cv{\mathbf{k}}}^{(1)}\mathbf{d}_{vc^{\prime}{\mathbf{k}}}\right]\,,$
(7) $\displaystyle
F^{(2)}_{vv^{\prime}{\mathbf{k}}}={\mathbf{E}}(t)\cdot\sum_{c}\left[\mathbf{d}_{vc{\mathbf{k}}}\rho_{cv^{\prime}{\mathbf{k}}}^{(1)}-\rho_{vc{\mathbf{k}}}^{(1)}\mathbf{d}_{cv^{\prime}{\mathbf{k}}}\right]\,.$
(8)

The $\rho^{(2)}_{cc^{\prime}{\mathbf{k}}}$
($\rho^{(2)}_{vv^{\prime}{\mathbf{k}}}$) terms are activated when both
$\Delta\epsilon_{cv{\mathbf{k}}}$ and
$\Delta\epsilon_{c^{\prime}v{\mathbf{k}}}$
($\Delta\epsilon_{cv^{\prime}{\mathbf{k}}}$) lay within
$\omega_{0}-\Delta_{0}$ and $\omega_{0}+\Delta_{0}$. In the $cv$ channel

$F^{(2)}_{cv{\mathbf{k}}}={\mathbf{E}}(t)\cdot\left[\sum_{c^{\prime}}\mathbf{d}_{cc^{\prime}{\mathbf{k}}}\rho_{c^{\prime}v{\mathbf{k}}}^{(1)}-\sum_{v^{\prime}}\rho_{cv^{\prime}{\mathbf{k}}}^{(1)}\mathbf{d}_{v^{\prime}v{\mathbf{k}}}\right]\,,$
(9)

which shows the importance of intra-band dipoles appearing already at second
order in the coupling with the external field. In the implementation of the
yambo_rt code, dipoles are set to zero whenever
$\Delta\epsilon_{nm{\mathbf{k}}}<\epsilon_{thresh}$ (and in particular for
$n=m$), due to the fact that the dipoles $\mathbf{d}_{nm{\mathbf{k}}}$ are ill
defined in such cases. More in general, neglecting all valence–valence and
conduction–conduction matrix elements we obtain the No–Intraband (NI) term

$F^{(2NI)}_{cv{\mathbf{k}}}=0$ (10)

and, due to boundary conditions, ${\rho^{(2)}_{cv{\mathbf{k}}}=0}$. Since we
are interested in the leading order terms in the pump pulse
$\rho_{cv{\mathbf{k}}}^{(1)}$ and $\rho^{(2)}_{cc^{\prime}{\mathbf{k}}}$ we
can ignore the role of intra-band dipoles for weak pumps.

####  II.2.2 HSEX approximation

We now want to check how much the update of the HSEX self-energy influences
the above analysis. $\Delta\Sigma^{HSEX}$ is linear in $\rho$ (since we do not
update the screening)

$\Delta\Sigma^{HSEX}[\rho(t)]_{nm{\mathbf{k}}}=K^{HSEX}_{nm{\mathbf{k}},\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}\rho_{\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}},$
(11)

where $K^{HSEX}=\partial\Sigma^{HSEX}/\partial G$ is the standard ai-BSE
kernel, however the indexes are, in general, not restricted to the $cv$
channel. $\Delta\Sigma^{HSEX}$ enters the EOM in a term of the form
$Q^{HSEX}=[\Delta\Sigma^{HSEX},\rho]$. The leading order term is different
from zero only in the cv channel

$Q^{(1)}_{cv{\mathbf{k}}}=(f^{eq}_{v{\mathbf{k}}}-f^{eq}_{c{\mathbf{k}}})K^{HSEX}_{cv{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}},$
(12)

This terms gives a correction into the EOM

$\Delta\epsilon_{cv{\mathbf{k}}}\rho^{(1)}_{cv{\mathbf{k}}}\rightarrow\left(\Delta\epsilon_{cv{\mathbf{k}}}\delta_{v,\tilde{v}}\delta_{c,\tilde{c}}\delta({\mathbf{k}}-\tilde{{\mathbf{k}}})+(f^{eq}_{v{\mathbf{k}}}-f^{eq}_{c{\mathbf{k}}})K^{HSEX}_{cv{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\right)\rho^{(1)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}$
(13)

and gives a renormalization of the energies from
$\Delta\epsilon_{cv{\mathbf{k}}}$ to the poles of the equilibrium BSE at zero
momentum $\omega_{\lambda\mathbf{0}}$. As a consequence the state generated by
the pump oscillates at the excitonic energies.

The second order term has the form

$Q^{(2)}_{nm{\mathbf{k}}}=K^{HSEX}_{n\tilde{l}{\mathbf{k}},\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{l}m{\mathbf{k}}}-\rho^{(1)}_{n\tilde{l}{\mathbf{k}}}K^{HSEX}_{\tilde{l}m{\mathbf{k}},\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}+(f^{eq}_{m{\mathbf{k}}}-f^{eq}_{n{\mathbf{k}}})K^{HSEX}_{nm{\mathbf{k}},\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}\rho^{(2)}_{\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}$
(14)

This couples the $cv$ channel with the $cc^{\prime}$ and $vv^{\prime}$
channels. We can write this term explicitly in the $cc^{\prime}$ and
$vv^{\prime}$ channels:

$\displaystyle Q^{(2)}_{cc^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle\left(K^{HSEX}_{c\tilde{v}^{\prime}{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{v}^{\prime}c^{\prime}{\mathbf{k}}}-\rho^{(1)}_{c\tilde{v}^{\prime}{\mathbf{k}}}K^{HSEX}_{\tilde{v}^{\prime}c^{\prime}{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\right)\rho^{(1)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}$
(15) $\displaystyle Q^{(2)}_{vv^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle\left(K^{HSEX}_{v\tilde{c}^{\prime}{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{c}^{\prime}v^{\prime}{\mathbf{k}}}-\rho^{(1)}_{v\tilde{c}^{\prime}{\mathbf{k}}}K^{HSEX}_{\tilde{c}^{\prime}v^{\prime}{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\right)\rho^{(1)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}$
(16)

where the last term of eq. (14) gives no contribution. Accordingly
$\rho^{(2)}_{cc^{\prime}{\mathbf{k}}}$ and
$\rho^{(2)}_{vv^{\prime}{\mathbf{k}}}$ become independent from the terms
$\rho^{(2)}_{cv{\mathbf{k}}}$. This means that, neglecting the intra-band
dipoles, $\rho^{(2)}_{cc^{\prime}{\mathbf{k}}}$ and
$\rho^{(2)}_{vv^{\prime}{\mathbf{k}}}$ remains exact up to second order also
at the TD-HSEX level. Moreover looking at these terms for the case
$c=c^{\prime}$ we understand why, at the HSEX level,
$N_{c}=\sum_{c}\rho^{(2)}_{cc{\mathbf{k}}}$ displays fast oscillations as long
as $\rho_{cv{\mathbf{k}}}\neq 0$, as shown in the manuscript.

###  II.3 Generated state in the continuum

Figure 4: Time dependent polarization obtained propagating the equation of
motion for $\rho(t)$, with the laser pulse in the continuum at the IP level
(left panel) and at the HSEX level (right panel).

We now analyze the state generated with the laser pulse frequency in the
continuum, comparing the result of the TD-IP propagation with a full TD-HSEX
propagation. As shown in Fig. 4 the total poplation display a similar
behavior, although the total number of carriers generated is higher (by
roughly a factor 2) in the TD-HSEX case. This is due to the fact that the TD-
HSEX (i.e. BSE) absorption is higher than the IP absorption in the same energy
range. We will analyze more in details the distribution of the population with
the TR-ARPES signal. Here we also observe that, at the IP level, the residual
polarization is lower since more frequencies exist at the IP level in the
energy range considered by the pump pulse.

##  III Time Resolved ARPES singal

###  III.1 Construction of the $G^{<}(t,t^{\prime})$ Green function

The GKBA expression for $G^{<}$ is

$\displaystyle
G^{<}_{cc^{\prime}{\mathbf{k}}}(t,t^{\prime})=\rho_{c\tilde{n}{\mathbf{k}}}(t)G^{(a)}_{\tilde{n}c^{\prime}{\mathbf{k}}}(t,t^{\prime})-G^{(r)}_{c\tilde{n}{\mathbf{k}}}(t,t^{\prime})\rho_{\tilde{n}c^{\prime}{\mathbf{k}}}(t^{\prime}).$
(17)

The retarded and the advanced propagator at the TD-HSEX level are given by

$\displaystyle G^{(r)}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\theta(t-t^{\prime})T\left[e^{-i\int_{t^{\prime}}^{t}h^{HSEX}[\rho(t)]dt}\right].$
(18)

The propagator in practice is constructed as

$\displaystyle G^{(r)}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\theta(t-t^{\prime})\prod_{n=0}^{N}e^{-ih^{HSEX}[\rho(t-n\,dt)]dt}$
(19)

with $t-N\,dt=t^{\prime}$. Starting from the time step $t=t^{\prime}$ where
$G^{(r)}_{nm{\mathbf{k}}}(t,t)=-i\delta_{nm}$ we have

$\displaystyle
G^{(r)}(t+dt,t^{\prime}-dt)=e^{-ih^{HSEX}[\rho(t+dt)]dt}G^{(r)}(t,t^{\prime})e^{-ih^{HSEX}[\rho(t^{\prime}-dt)]dt}$
(20)

In the actual implementation I replace

$\displaystyle h^{HSEX}[\rho(t+dt)]$ $\displaystyle\rightarrow$ $\displaystyle
0.5\,(h^{HSEX}[\rho(t+dt)]-h^{HSEX}[\rho(t)])$ (21) $\displaystyle
h^{HSEX}[\rho(t^{\prime}-dt)]$ $\displaystyle\rightarrow$ $\displaystyle
0.5\,(h^{HSEX}[\rho(t^{\prime}-dt)]-h^{HSEX}[\rho(t^{\prime})]).$ (22)

The exponential $e^{-ih^{HSEX}[\rho(t)]dt}$ is approximated with a 5th order
expansion in $dt$ and, for each term on the left and on the right, a matrix-
matrix multiplication in the $nm$ space is performed. This is repeated for
each ${\mathbf{k}}$-point for each time-step. The Hamiltonian is defined as
(we consider only the non overlapping regime, i.e. we probe time ranges where
the pump pulse is off)

$h^{HSEX}[\rho(t)]=h^{eq}+\Delta\Sigma^{HSEX}[\rho(t)].$ (23)

In general

$\Delta\Sigma^{HSEX}[\rho(t)]_{nm{\mathbf{k}}}=K^{HSEX}_{nm{\mathbf{k}},\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}(\rho_{\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}(t)-\rho^{eq}_{\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}).$
(24)

The sum over ${\tilde{n}\tilde{m}\tilde{{\mathbf{k}}}}$ indexes is implicit.

To analyze the expression, we write again the perturbative expansion (this is
not done in the implementation). As discussed previously
$\rho_{cc^{\prime}{\mathbf{k}}}(t)\approx\rho^{(2)}_{cc^{\prime}{\mathbf{k}}}(t)$
(in the $vv^{\prime}$ channel instead
$\rho_{vv^{\prime}{\mathbf{k}}}(t)\approx\delta_{v,v^{\prime}}f^{(eq)}_{v{\mathbf{k}}}$),
while $\rho_{cv{\mathbf{k}}}(t)\approx\rho_{cv{\mathbf{k}}}^{(1)}(t)$. For the
self-energy and we consider the leading order terms only. If we restrict to
the $cv$ channel,

$\Delta\Sigma^{HSEX}[\rho(t)]_{cv{\mathbf{k}}}\approx
K^{HSEX}_{cv{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\rho^{(1)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}.$
(25)

where $K^{HSEX}$ is the standard ai-BSE kernel. The diagonal terms
$\Sigma_{nn{\mathbf{k}}}$ can be absorbed in a renormalization of the energies
entering the $h^{eq}$ hamiltonian. The propagators have leading terms
$G^{(r/a)}_{nn{\mathbf{k}}}(t,t^{\prime})\approx
G^{(r/a),(eq)}_{nn{\mathbf{k}}}(t,t^{\prime})$ and
$G^{(r/a)}_{cv{\mathbf{k}}}(t,t^{\prime})\approx
G^{(r/a),(1)}_{cv{\mathbf{k}}}(t,t^{\prime})$. The full $G^{<}$ can be split
into two terms, one with contribution from the terms of the density matrix
$\rho_{cc^{\prime}{\mathbf{k}}}(t)$ and one with contribution from the off-
diagonal elements $\rho_{cv{\mathbf{k}}}(t)$ Interestingly both terms are of
second order in the pump pulse:

$\displaystyle G^{<,(2a)}_{cc^{\prime}{\mathbf{k}}}(t,t^{\prime})$
$\displaystyle=$
$\displaystyle\rho^{(2)}_{c\tilde{c}{\mathbf{k}}}(t)G^{(a),(eq)}_{\tilde{c}c^{\prime}{\mathbf{k}}}(t,t^{\prime})-G^{(r),(eq)}_{c\tilde{c}{\mathbf{k}}}(t,t^{\prime})\rho^{(2)}_{\tilde{c}c^{\prime}{\mathbf{k}}}(t^{\prime}),$
(26) $\displaystyle G^{<,(2b)}_{cc^{\prime}{\mathbf{k}}}(t,t^{\prime})$
$\displaystyle=$
$\displaystyle\rho^{(1)}_{c\tilde{v}{\mathbf{k}}}(t)G^{(a),(1)}_{\tilde{v}c^{\prime}{\mathbf{k}}}(t,t^{\prime})-G^{(r),(1)}_{c\tilde{v}{\mathbf{k}}}(t,t^{\prime})\rho^{(1)}_{\tilde{v}c^{\prime}{\mathbf{k}}}(t^{\prime}),$
(27)

The signal for the case $\rho_{cv{\mathbf{k}}}=0$ is fully determined by eq.
(26), while the complete signal has also contribution from eq. (27). It is
remarkable to notice that, in the coherent exciton case, the contribution eq.
(27) not only gives the excitonic ARPES singal in the band gap, but is also
responsible for exactly removing the qp ARPES singal from the conduction band
due to eq. (26).

Figure 5: The TR-ARPES spectra function after partial de-phasing of the off-
diagonal elements of the density matrix. The numbers reported in the PE and
ARPES panels represent the magnification factor needed for the signal in the
conduction band region to reach similar intensity as the signal from the
valence band.

Here we also show (see Fig. 5) the situation in which the matrix elements of
the density matrix are partially dephased. In such case the cancellation
between the two contributions to $G^{<}_{cc{\mathbf{k}}}$ is not exact anymore
and both the signal from the excitonic state and from the carriers are
visible. This is the intermediate situation between the coherent excitonic
case and the non coherent carriers population shown in the manuscript. It is
reported here to highlight the importance of the cancellation between such two
terms. As discussed in the main text the procedure of sending to zero the off-
diagonal elements of the density matrix changes the total energy of the system
and it is thus not physical. This is somehow evident also from the ARPES
spectral function since lower photon energy is required to extract electrons
after these elements are sent to zero.

####  III.1.1 Interpolation of $G^{<}(t,t^{\prime})$ in ${\mathbf{k}}$-space

Figure 6: Integrated and $k$-resoved spectral function for LiF after the
action of an ultra short laser pulse with frequency resonant to the excitonic
peaks. Two different k-point samplings are used $G^{<}(\omega)$. The numbers
reported in the PE and ARPES panels represent the magnification factor needed
for the signal in the conduction band region to reach similar intensity as the
signal from the valence band.

The $G^{<}_{cc{\mathbf{k}}}(\omega)$ function is constructed after the TD-HSEX
propagation on a regular grid in ${\mathbf{k}}$-space. To plot the signal on a
given path a Fourier interpolation is performed (this is possible since the
physical signal is phase independent). However the interpolation of a function
is highly non trivial. We solve this doing a frequency by frequency
interpolation. This is much less efficient since the resulting signal tends to
be stretched along the fixed frequency axis. In this case it works reasonably
well since a very fine ${\mathbf{k}}$-points sampling is needed to converge
the absorption spectrum.

In Fig. 6 we show the interpolated signal starting from the $16\times 16\times
16$ sampling (left panel) and from the $24\times 24\times 24$. The excitonic
signal is very smooth already with the smaller grid, while the signal from the
valence band from the $16^{3}$ grid suffers from the stretching above
mentioned.

###  III.2 TR-ARPES signal in the continuum

Figure 7: Integrated and $k$-resoved spectral function for LiF after the
action of an ultra short laser pulse in the continuum. The $k$-resoved signal
is obtained interpolating $G^{<}({\mathbf{k}},\omega)$ frequency by frequency
from a regular grid. The band structure is superimposed to the signal (red
lines). A replica of the valence band shifted by the excitonic energy is also
represented with white lines. The three panels represent TD-IP simultion (left
panel) TD-HSEX simulation (central panel), and TD-HSEX simulation with
dephasing (right panel).

In Fig. 7 we compare the TR-ARPES signal from a system for the case
$\omega_{P_{a}}=14.8$ eV at the TD-IP and the TD-HSEX level. For the TD-HSEX
level we consider both the the result of the simulation without dephasing and
the result after sending to zero the off-diagonal matrix-elements of the
density matrix. The obtained signal is almost identical in the three cases,
with two peaks located at the point where the up-shifted valence band crosses
the conduction band. As for the case with the laser pulse resonant with the
exciton, the TD-HSEX solution predicts that the signal is a replica of the
valence band weighted by the BSE eigenvector which has in principle no direct
relation to the conduction band. However, in this energy range, the BSE
eigenvectors in the continuum $A^{\lambda{\mathbf{q}}}$ are of almost delta
like nature for LiF, i.e.
$A^{\lambda{\mathbf{q}}}\approx\delta_{v,v_{P_{a}}}\delta_{c,c_{P_{a}}}\delta({\mathbf{k}}-{\mathbf{k}}_{P_{a}})$
and the expression of the ARPES spectral function reduces to the IP expression
with
$\epsilon_{c{\mathbf{k}}}=(\epsilon_{v{\mathbf{k}}-{\mathbf{q}}}+\omega_{\lambda{\mathbf{q}}})$.
Numerically this results can be approached only very slowly as discussed
above. The results from this section show that our approach is completely
general and it works well in the excitonic case, but also in the continuum
where however a simpler independent particles approach and interpretation in
terms of carriers moving on the band structure works fine.

There is a small shift in ${\mathbf{k}}$–space between the IP and HSEX
approach (panel $a$ and panel $b$) which turns into an energy shift once
dephasing is included in the TD-HSEX simulation. The source of the shifts in
Fig. 7 can be understood considering these points:

  * (i)

the energy of the photo–emitted signal must be the energy from the valence
band, plus the energy of the laser pump $\epsilon_{v{\mathbf{k}}}+\omega_{0}$
if we do not artificially change the energy of the system (panels $a$ and $b$)

  * (ii)

the energy of the photo–emitted electron must lay on the conduction band
structure when it depends only on $G^{<,2a}_{cc{\mathbf{k}}}(\omega)$ (panels
$a$ and $c$)

  * (iii)

the dephasing alters the energy of the system. The energy of the photo–emitted
electrons is shifted from $\epsilon_{v{\mathbf{k}}}+\omega_{0}$ with
$\omega_{0}\approx E^{\lambda}$, to $\epsilon_{c{\mathbf{k}}}$, i.e. the same
shift discussed in the analysis of $F^{\lambda}$. However it does not change
where the signal lies in ${\mathbf{k}}$–space, since it just kills
$G^{<,2b}_{cc{\mathbf{k}}}(\omega)$, thus the signal from TD-HSEX simulation
always lies in the same ${\mathbf{k}}$–region (panels $a$ and $b$).

We expect all these shifts to go to zero at convergence, and thus regard them
as non-physical. However this should be fully verified. The existence of this
shifts is the main reason why in the main text we prefer to present the signal
obtained from within the IP propagation for the case in the continuum.

Figure 8: Integrated and $k$-resoved spectral function for LiF after the
action of an ultra short laser pulse with frequency in the continuum (bottom
panels). Two different time windows are considered to reconstruct
$G^{<}(\omega)$.

We conclude our analysis looking at the relative intensity of the two peaks.
Despite we are in the continuum with a very high density of poles in our
simulation, only few are optically active (see Fig. 2). This is reflected in
the behavior of the polarization, as shown in Fig. 8. This points to the fact
that, at variance with the case where we pump resonant with the excitonic
frequency, the PE signal may depend on the time window considered to construct
it. To address this we consider the spectrum generated into two different
energy windows: $160-200fs$ and $200-250fs$. Both spectra consist of two peaks
located in the same point. However, by changing the time window used to
extract the spectral function, the relative intensity of the two peaks changes
(see fig. 8).

##  IV Transient Absorption signal

###  IV.1 Perturbative expansion in the probe

####  IV.1.1 IP approximation

We now proceed with the perturbarbative expansion defined in the previous
section. We are interested in the linear response to the probe
${\mathbf{e}}(t)$. Since the signal generated by the pump pulse has mixed
first and second order character, the TrAbs signal will have in general mixed
second and third order character. After the action of the pump pulse the
density matrix contains both diagonal and off-diagonal elements and the third
order source terms have the general expression. In the $cc^{\prime}$ and
$vv^{\prime}$ channels we have

$\displaystyle F^{(1p)}_{cc^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{v}\left(\mathbf{d}_{cv{\mathbf{k}}}\rho_{vc^{\prime}{\mathbf{k}}}^{(1+2)}-\rho_{cv{\mathbf{k}}}^{(1+2)}\mathbf{d}_{vc^{\prime}{\mathbf{k}}}\right)+\sum_{c^{\prime\prime}}\left(\mathbf{d}_{cc^{\prime\prime}{\mathbf{k}}}\rho_{c^{\prime\prime}c^{\prime}{\mathbf{k}}}^{(2)}-\rho_{cc^{\prime\prime}{\mathbf{k}}}^{(2)}\mathbf{d}_{c^{\prime\prime}c^{\prime}{\mathbf{k}}}\right)\Bigg{]}\,.$
(28) $\displaystyle F^{(1p)}_{vv^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{c}\left(\mathbf{d}_{vc{\mathbf{k}}}\rho_{cv^{\prime}{\mathbf{k}}}^{(1+2)}-\rho_{vc{\mathbf{k}}}^{(1+2)}\mathbf{d}_{cv^{\prime}{\mathbf{k}}}\right)+\sum_{v^{\prime\prime}}\left(\mathbf{d}_{vv^{\prime\prime}{\mathbf{k}}}\rho_{v^{\prime\prime}v^{\prime}{\mathbf{k}}}^{(0+2)}-\rho_{vv^{\prime\prime}{\mathbf{k}}}^{(0+2)}\mathbf{d}_{v^{\prime\prime}v^{\prime}{\mathbf{k}}}\right)\Bigg{]}\,.$
(29)

In the $cv$ channel

$F^{(1p)}_{cv{\mathbf{k}}}={\mathbf{e}}(t)\cdot\Bigg{[}\sum_{c^{\prime}}\left(\mathbf{d}_{cc^{\prime}{\mathbf{k}}}\rho_{c^{\prime}v{\mathbf{k}}}^{(1+2)}-\rho_{cc^{\prime}{\mathbf{k}}}^{(2)}\mathbf{d}_{c^{\prime}v{\mathbf{k}}}\right)+\sum_{v^{\prime}}\left(\mathbf{d}_{cv^{\prime}{\mathbf{k}}}\rho_{v^{\prime}v{\mathbf{k}}}^{(0+2)}-\rho_{cv^{\prime}{\mathbf{k}}}^{(1+2)}\mathbf{d}_{v^{\prime}v{\mathbf{k}}}\right)\Bigg{]}\,,$
(30)

In both channels intra-band terms are needed. Not only the coupling with the
external pulse needs to include them, but also the dipoles are needed in the
definition of the polarization
${\mathbf{P}}=\sum_{nm}{\mathbf{d}}_{nm{\mathbf{k}}}\rho_{nm{\mathbf{k}}}$
from which the TrAbs signal is defined.

In the continuum the oscillatory terms quickly dephase and one can neglect
them. This gives:

$\displaystyle F^{(1p-DE)}_{cc^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{c^{\prime\prime}}\left(\mathbf{d}_{cc^{\prime\prime}{\mathbf{k}}}\rho_{c^{\prime\prime}c^{\prime}{\mathbf{k}}}^{(2)}-\rho_{cc^{\prime\prime}{\mathbf{k}}}^{(2)}\mathbf{d}_{c^{\prime\prime}c^{\prime}{\mathbf{k}}}\right)\Bigg{]}\,.$
(31) $\displaystyle F^{(1p-DE)}_{vv^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{v^{\prime\prime}}\left(\mathbf{d}_{vv^{\prime\prime}{\mathbf{k}}}\rho_{v^{\prime\prime}v^{\prime}{\mathbf{k}}}^{(0+2)}-\rho_{vv^{\prime\prime}{\mathbf{k}}}^{(0+2)}\mathbf{d}_{v^{\prime\prime}v^{\prime}{\mathbf{k}}}\right)\Bigg{]}\,.$
(32) $\displaystyle F^{(1p-DE)}_{cv{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{c^{\prime}}\left(-\rho_{cc^{\prime}{\mathbf{k}}}^{(2)}\mathbf{d}_{c^{\prime}v{\mathbf{k}}}\right)+\sum_{v^{\prime}}\left(\mathbf{d}_{cv^{\prime}{\mathbf{k}}}\rho_{v^{\prime}v{\mathbf{k}}}^{(0+2)}\right)\Bigg{]}\,,$
(33)

Thus the source term for the density matrix dynamics, which in turn determines
the transient absorption signal, depends on non–equilibrium occupations only.
The zero order term $\rho_{v^{\prime}v{\mathbf{k}}}^{(0)}$ is removed by
subtracting the equilibrium contribution. The $cc^{\prime}$ and $vv^{\prime}$
terms still give contribution to the polarization. However neglecting the
dipoles whenever $\Delta\epsilon_{nm{\mathbf{k}}}<\epsilon_{thresh}$ just
removes the low frequency terms, since $\rho_{cc^{\prime}{\mathbf{k}}}\propto
e^{-i\Delta\epsilon{cc^{\prime}{\mathbf{k}}}}$ (and similarly
$\rho_{vv^{\prime}{\mathbf{k}}}\propto
e^{-i\Delta\epsilon{vv^{\prime}{\mathbf{k}}}}$) and the transient absorption
for $\omega>\Delta\epsilon_{thresh}$ remains unaffected.

In the situation where, instead, there is an isolated pole, the
$\rho_{cv{\mathbf{k}}}$ terms can survive for longer times and the complete
expression for the source terms need to be used. Neglecting the intra-band
dipoles (and using, accordingly, $\rho^{(2)}_{cv{\mathbf{k}}}=0$) we obtain

$\displaystyle F^{(1p-NI)}_{cc^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{v}\left(\mathbf{d}_{cv{\mathbf{k}}}\rho_{vc^{\prime}{\mathbf{k}}}^{(1)}-\rho_{cv{\mathbf{k}}}^{(1)}\mathbf{d}_{vc^{\prime}{\mathbf{k}}}\right)\Bigg{]}\,.$
(34) $\displaystyle F^{(1p-NI)}_{vv^{\prime}{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{c}\left(\mathbf{d}_{vc{\mathbf{k}}}\rho_{cv^{\prime}{\mathbf{k}}}^{(1)}-\rho_{vc{\mathbf{k}}}^{(1)}\mathbf{d}_{cv^{\prime}{\mathbf{k}}}\right)\Bigg{]}\,.$
(35) $\displaystyle F^{(1p-NI)}_{cv{\mathbf{k}}}$ $\displaystyle=$
$\displaystyle{\mathbf{e}}(t)\cdot\Bigg{[}\sum_{c^{\prime}}\left(-\rho_{cc^{\prime}{\mathbf{k}}}^{(2)}\mathbf{d}_{c^{\prime}v{\mathbf{k}}}\right)+\sum_{v^{\prime}}\left(\mathbf{d}_{cv^{\prime}{\mathbf{k}}}\rho_{v^{\prime}v{\mathbf{k}}}^{(0+2)}\right)\Bigg{]}\,,$
(36)

All leading order terms are included, thus we expect the intra-band dipoles
not to have an important role in the coupling with the external field. The
resulting source term have contributions from first order terms in the pump
pulse, at variance with the dephased case, where only second order terms are
present. The first order terms affects only the $cc^{\prime}$ and
$cv^{\prime}$ channels. However, this may lead to oscillations in the $cv$
energy range, since $\rho_{cv{\mathbf{k}}}$ enters in such source terms. Thus
intra-band dipoles may have an important role in the definition of the
polarization. This is why a complete discussion on TrAbs signal may need to
account for intra-band dipoles in the excitonic case.

####  IV.1.2 HSEX approximation

We are interested to the HSEX self-energy to linear order in the probe pulse,
with a first term coming from
$[\Delta\Sigma^{(1+2)},\rho^{(1p)}]_{nm{\mathbf{k}}}$ and the second from
$[\Delta\Sigma^{(1p)},\rho^{(0+1+2)}]_{nm{\mathbf{k}}}$

$G^{(1p)}_{nm{\mathbf{k}}}=K^{HSEX}_{nl^{\prime}{\mathbf{k}},n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(1+2)}_{n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(1p)}_{l^{\prime}m{\mathbf{k}}}-\rho^{(1p)}_{nl^{\prime}{\mathbf{k}}}K^{HSEX}_{l^{\prime}m{\mathbf{k}},n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(1+2)}_{n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}+\\\
+K^{HSEX}_{nl^{\prime}{\mathbf{k}},n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(1p)}_{n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(0+1+2)}_{l^{\prime}m{\mathbf{k}}}-\rho^{(0+1+2)}_{nl^{\prime}{\mathbf{k}}}K^{HSEX}_{l^{\prime}m{\mathbf{k}},n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(1p)}_{n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}$
(37)

where here $\rho^{(0+1+2)}$ is the density matrix up to second order in the
pump, while $\rho^{(1p)}$ is the density matrix also including the first order
correction due to the probe. Taking the $cv$ channel, the second term can be
re-written as
$\tilde{K}^{P}_{nm{\mathbf{k}},n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}\rho^{(1p)}_{n^{\prime}m^{\prime}{\mathbf{k}}^{\prime}}$.
The first term instead contains
$(\Delta\Sigma^{(1+2)}_{cc{\mathbf{k}}}-\Delta\Sigma^{(1+2)}_{vv{\mathbf{k}}})\rho^{(1p)}_{cv{\mathbf{k}}}\equiv\Delta\tilde{\Sigma}^{P}_{cv{\mathbf{k}}}\rho^{(1p)}_{cv{\mathbf{k}}}$
which can be interpreted as a renormalization of the QP band structure, plus
$\Delta\Sigma^{(1+2)}_{cv{\mathbf{k}}}(\rho^{(1p)}_{vv{\mathbf{k}}}-\rho^{(1p)}_{cc{\mathbf{k}}})$
which points to the correlation between diagonal and off diagonal terms and
the role of intraband dipoles. Neglecting this last term we have, w.r.t. the
equilibrium case the change induced by the pump

$\left(\Delta\epsilon_{cv{\mathbf{k}}}\delta_{v,\tilde{v}}\delta_{c,\tilde{c}}\delta({\mathbf{k}}-\tilde{{\mathbf{k}}})+(f^{eq}_{v{\mathbf{k}}}-f^{eq}_{c{\mathbf{k}}})K^{HSEX}_{cv{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\right)\rho^{(1p)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\rightarrow\left((\Delta\epsilon_{cv{\mathbf{k}}}+\Delta\tilde{\Sigma}^{P}_{cv{\mathbf{k}}})\delta_{v,\tilde{v}}\delta_{c,\tilde{c}}\delta({\mathbf{k}}-\tilde{{\mathbf{k}}})+\tilde{K}^{P}_{cv{\mathbf{k}},\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}\right)\rho^{(1p)}_{\tilde{c}\tilde{v}\tilde{{\mathbf{k}}}}$
(38)

###  IV.2 On the role of intraband dipoles

Figure 9: Transient Absorption for LiF shown for the system driven with a
pulse resonant to the excitonic peak. In all cases the states are generate
with a TD-HSEX propagation up to 160 fs. Starting from the generated density
the system is then probed with a delta like pulse. In the time range 160-260
fs, the propagation in the first and second panel from the left is done
including intra-band dipoles. In the first panel the spectrum is reconstructed
starting from the polarization, in the second starting from the current.
Finally in the third panel the propagation is done setting to zero the intra-
band dipoles, and then computing the spectrum from the polarization. The probe
induced polarization
${\mathbf{P}}^{(i)}(t)={\mathbf{P}}_{Pp}(t)-{\mathbf{P}}_{P}(t)$ is shown in
the inset. The blue shaded line marks the laser pump frequency, the black
dashed line the position of the quasi-particle gap. The extra insets show a
zoom of $\varepsilon^{(eq)}_{2}(\omega)$ (black line) and
$\varepsilon^{(i)}_{2}(\omega)$ around the two main excitonic peaks.

To verify the role of intraband dipoles, we compare in Fig. 9 the transient
absorption signal using three different schemes which deal in different ways
with the intra-band dipoles. The first panel shows the transient absorption
generated from the Fourier transform of ${{\mathbf{P}}_{Pp}-{\mathbf{P}}_{P}}$
including intra-band diploles. Since such dipoles are ill defined for small
frequencies the non-equilibrium signal in presence of the coherent excitonic
state has a non physical behavior. This appears as a $1/\omega$ decay in the
detected energy range of the plot, which is due to a non-physical pole
appearing at $\omega_{X}\approx 4.11$ eV (not shown) with huge intensity
($\epsilon(\omega_{X})\approx 200$). The existence of such pole is responsible
for the behavior of the polarization, which is dominated by a characteristic
frequency in the first panel. If instead the spectrum is reconstructed
starting from the current ${{\mathbf{J}}_{Pp}-{\mathbf{J}}_{P}}$ (central
panel), such pole disappears. The current based approach is defined in terms
of the expectation value of the velocity operator and the intra–band dipoles
are well defined in such case. However the velocity operator should in
principle contain also the commutator with the HSEX self-energy which is
neglected here. As a last step we compute again the transient absorption
starting from the ${{\mathbf{P}}_{Pp}-{\mathbf{P}}_{P}}$, but setting to zero
all intra-band dipoles. The signal is quite similar to the one obtained in the
current based approach. In particular in both cases the differences between
the coherent excitonic case and the non coherent populations discussed in the
main text remain the same. Further checks on the role of the intra–band
dipoles could be achieved via the NEQ Berry Phase formulation.

###  IV.3 Tr-Abs signal in the continuum

Figure 10: Transient Absorption for LiF shown for the system driven in the
continuum. In both cases the states are generate with a TD-IP propagation up
to 160 fs. Starting from the generated density the system is then probed with
a delta like pulse propagating either the TD-IP EOM (left panel) and the TD-
HSEX EOM (right panel). The response is completly different, showing that
excitonic effects cannot be neglected in the probing step. In both cases the
role of coherences is negligible.

In this section we discuss in detail the idea of using an IP propagation for
the generation of the coherent carriers state. To this end we compare the
TrAbs obtained with IP pump - IP probe, IP pump - HSEX probe, and HSEX pump -
HSEX probe. In the IP - IP case the signal is due to Pauli blocking or
bleaching of the transitions, mostly the one at 14.8 eV. Other transitions
which share either the conduction band or the valence band with the transition
at 14.8 eV are also slightly bleached. All this can be easely interpreted in
terms of $\chi^{IP}[f_{n{\mathbf{k}}}(\tau)](\omega)$ with
$f_{n{\mathbf{k}}}=\rho_{nn{\mathbf{k}}}$. The signal completly neglect
excitonic effects in the pump process. This is not a good approximation, as
discussed in the main text. Instead, the IP - HSEX signal and the HSEX - HSEX
signal are quite similar (the main difference is the overall magnitude),
confirming the the IP approximation for the generation of the state is a
reasonable approximation.

###  IV.4 On the structure of the response function

I discuss here the resonant and non resonat term of the response function and
the relation between these two terms and the absorption and emission channel.
I first start writing the IP response function for a cold semiconductor, i.e.
with $f_{v{\mathbf{k}}}=1$ and $f_{c{\mathbf{k}}}=0$ for valence and
conduction respectively. The expression of a general response function
contains only $v\rightarrow c$ transitions:

$\chi^{AB,eq}(\omega)=\sum_{cv{\mathbf{k}}}\frac{(\langle
v{\mathbf{k}}|\hat{A}|c{\mathbf{k}}\rangle\langle
c{\mathbf{k}}|\hat{B}|v{\mathbf{k}}\rangle)f_{v{\mathbf{k}}}(1-f_{c{\mathbf{k}}})}{\omega-(\epsilon_{c{\mathbf{k}}}-\epsilon_{v{\mathbf{k}}})+i\eta}-\frac{(\langle
v{\mathbf{k}}|\hat{B}|c{\mathbf{k}}\rangle\langle
c{\mathbf{k}}|\hat{A}|v{\mathbf{k}}\rangle)f_{v{\mathbf{k}}}(1-f_{c{\mathbf{k}}})}{\omega+(\epsilon_{c{\mathbf{k}}}-\epsilon_{v{\mathbf{k}}})+i\eta}.$
(39)

The first term, or resonant term, describes absorption and has poles for
$\omega>0$. The second term, or anti-resonant term, would describe emission,
but at equilibrium has poles only for $\omega<0$, i.e. no emission takes place
from the ground state.

If we now consider the situation where electrons are injected from valence to
conduction, i.e. $f_{v{\mathbf{k}}}<1$ and $f_{c{\mathbf{k}}}>0$ we add extra
terms which contain new $c\rightarrow v$ transitions. These can be constructed
by swapping $c\leftrightarrow v$ in the previous expression (I neglect for the
present discussion $c\rightarrow c^{\prime}$ and $v\rightarrow v^{\prime}$)
transitions

$\chi^{AB,emiss}(\omega)=\sum_{cv{\mathbf{k}}}\frac{(\langle
c{\mathbf{k}}|\hat{A}|v{\mathbf{k}}\rangle\langle
v{\mathbf{k}}|\hat{B}|c{\mathbf{k}}\rangle)f_{c{\mathbf{k}}}(1-f_{v{\mathbf{k}}})}{\omega-(\epsilon_{v{\mathbf{k}}}-\epsilon_{c{\mathbf{k}}})+i\eta}-\frac{(\langle
c{\mathbf{k}}|\hat{B}|v{\mathbf{k}}\rangle\langle
v{\mathbf{k}}|\hat{A}|c{\mathbf{k}}\rangle)f_{c{\mathbf{k}}}(1-f_{v{\mathbf{k}}})}{\omega+(\epsilon_{v{\mathbf{k}}}-\epsilon_{c{\mathbf{k}}})+i\eta}.$
(40)

In this new terms the resonant term has poles for $\omega<0$, i.e. there is no
extra contribution to the absorption from these transitions, while the anti-
resonant term has now poles for $\omega>0$, and indeed emission can take place
from these new term. Considering only the $\omega>0$ region we thus have, for
valence to conduction transitions

$\chi^{AB,\omega>0}(\omega)=\sum_{cv{\mathbf{k}}}\frac{(\langle
v{\mathbf{k}}|\hat{A}|c{\mathbf{k}}\rangle\langle
c{\mathbf{k}}|\hat{B}|v{\mathbf{k}}\rangle)f_{v{\mathbf{k}}}(1-f_{c{\mathbf{k}}})}{\omega-(\epsilon_{c{\mathbf{k}}}-\epsilon_{v{\mathbf{k}}})+i\eta}-\frac{(\langle
c{\mathbf{k}}|\hat{B}|v{\mathbf{k}}\rangle\langle
v{\mathbf{k}}|\hat{A}|c{\mathbf{k}}\rangle)f_{c{\mathbf{k}}}(1-f_{v{\mathbf{k}}})}{\omega-(\epsilon_{c{\mathbf{k}}}-\epsilon_{v{\mathbf{k}}})+i\eta}.$
(41)

We see that, compared to equilibrium where
$f_{v{\mathbf{k}}}(1-f_{c{\mathbf{k}}})=1$ we have a reduction of absorption
due to Pauli blocking, since now $f_{v{\mathbf{k}}}(1-f_{c{\mathbf{k}}})<1$
and, moreover a further reduction the peak intensity due to the second term
with $f_{c{\mathbf{k}}}(1-f_{v{\mathbf{k}}})>0$ which is multiplied by $-1$.
The overall weight of the transition if
$(f_{v{\mathbf{k}}}-f_{c{\mathbf{k}}})$, as it should.

The same kind of structure can be constructed from the BSE propagator and
leads the correlated version of eq. (41), i.e. eq. (6) in the main text.
However in the bosonic case, as discussed, the absorption is not reduced but
enhanced and this exactly compensate the emission term. Notice that, similarly
to the fact that here $c\rightarrow c^{\prime}$ and $v\rightarrow v^{\prime}$
transitions are neglected here, in the main text
$\omega_{\lambda{\mathbf{q}}}\rightarrow\omega_{\lambda{\mathbf{q}}^{\prime}}$
transitions are neglected.

Let me add few final comments here. The way the bosonic occupations enter eq.
(6) in the main text, with $N_{\lambda{\mathbf{q}}}+1$ in one term, and
$N_{\lambda{\mathbf{q}}}$ in the other, is the “standard expression”, which,
for the electron–magnetic field, is associated to spontaneous emission,
stimulated emission and absorption of photons. Here the names are swapped. We
call ”absorption term” the one which contains the $N_{\lambda{\mathbf{q}}}+1$
factor, since we keep referring to the absorption of photons. Indeed
$N_{\lambda{\mathbf{q}}}+1$ factor describes the creation of excitons, i.e. it
could be called exciton emission term, and similarly the second term could be
called exciton absorption term, since it describes the destruction of an
exciton. In the electron–magnetic case the +1 term is associated to a
“spontaneous” process since we start from the electronic system, which
interacts with the photon field, in an excited state. Similarly the emission
of excitons is spontaneous (i.e. it does not need the existence of a coherent
excitonic source to take place) if we start from the photonic system, which
interacts with the “exciton field”, in an excited state. In this jargon there
would be available photons (i.e. the laser pulse) which can “decay” by
creating excitons to conserve the total energy of the system. There is one
case however, where the existence of photon populations is not needed: when
$\omega_{\lambda}<0$. This is the case where the EI is spontaneously created.
In the present manuscript I’m considering the case with $\omega_{\lambda}>0$
and the NEQ-EI is created via the presence of a laser pulse.


# Spectral Form Factors of Topological Phases

Anurag Sarkar<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology, Kanpur 208016, India  Subrata Pachhal<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology, Kanpur 208016, India
Adhip Agarwala<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology, Kanpur 208016, India  Diptarka Das<EMAIL_ADDRESS>Department of
Physics, Indian Institute of Technology, Kanpur 208016, India

###### Abstract

Signatures of dynamical quantum phase transitions and chaos can be found in
the time evolution of generalized partition functions such as spectral form
factors (SFF) and Loschmidt echos. While a lot of work has focused on the
nature of such systems in a variety of strongly interacting quantum theories,
in this work, we study their behavior in short-range entangled topological
phases - particularly focusing on the role of symmetry protected topological
zero modes. We show, using both analytical and numerical methods, how the
existence of such zero modes in any representative system can mask the SFF
with large period (akin to generalized Rabi) oscillations hiding any behavior
arising from the bulk of the spectrum. Moreover, in a quenched disordered
system these zero modes fundamentally change the late time universal behavior
reflecting the chaotic signatures of the zero energy manifold. Our study
uncovers the rich physics underlying the interplay of chaotic signatures and
topological characteristics in a quantum system.

Introduction: Ideas of thermalization and chaos have become all pervasive in
multiple sub-disciplines of physics, including classical and quantum many body
systems, quantum field theory, gravity and fluidsSrednicki (1994); D’Alessio
_et al._ (2016); Borgonovi  _et al._ (2016); Das  _et al._ (2018); Pollack
_et al._ (2020); Murugan  _et al._ (2021). These questions become crucial to
understand the phases of interacting quantum systems when they are either
driven or are coupled to baths, however, methods of characterizing chaos are
few and limited. The behavior of spectral form factor (SFF) in such systems
has been particularly illuminating. Most interestingly, one finds that SFFs
follow universal features independent of underlying microscopic details. SFFs
and its variants such as Loschmidt echo del Campo  _et al._ (2017), Fisher
zeros Bena  _et al._ (2005); Harrow  _et al._ (2020) etc., which are
generalizations of the equilibrium thermal partition function, have been
vigorously investigated for various quantum many body systems ranging from
lattice models to black holes Cotler  _et al._ (2017a); Saad  _et al._ (2018);
Dyer and Gur-Ari (2017). It is known that interacting many-body chaotic
systems show a dip-linear ramp structure in the SFF, which saturates at late
timesCotler  _et al._ (2017b); Liu (2018). This behavior however has
intriguing micro-structures which depends on the underlying symmetries, nature
of interactions, dimensionality and localization properties del Campo  _et
al._ (2017). However, the behavior of SFF vis-a-vis topological properties of
a quantum Hamiltonian has been little explored.

In this paper, we investigate the behavior of SFF in systems where the
Hamiltonians describe a symmetry protected topological phase. Topological
phases of matter, unlike symmetry broken systems, identify different
thermodynamic phases with different topological invariants which are often
protected by some discrete symmetries Ludwig (2016); Chiu  _et al._ (2016);
Hasan and Kane (2010); Qi and Zhang (2011). A paradigmatic model for
investigating topological systems is the Su-Schrieffer-Heeger (SSH) model Su
_et al._ (1979) where the system hosts chiral symmetry protected edge modes in
the topological phase. In this work we show that in general, existence of such
zero dimensional topologically protected eigenstates can fundamentally
transform the SFF showing large time oscillations. We show both numerically
and analytically, that these are generalizations of Rabi oscillations, where
the SFF locks between a few states of the complete many-body spectrum. We
further show that this in fact is a property of even higher order topological
phases Benalcazar  _et al._ (2017), where a $d$ dimensional topological phase
hosts $(d-2)$ dimensional boundary modes, with characteristic features in the
SFF.

Figure 1: (a) A topological phase achieved in a Hamiltonian system (see eqn.
(3) for SSH Model) is tuned by a parameter $v/w$. This sets an energy scale
$T^{*}=\frac{1}{\beta^{*}}$ (shown by a dashed line) below which the system
shows prominent oscillations (see (b), with $L=60,v/w=0.5,\beta=20$) in its
SFF. The same system in the trivial region ($v/w=2,\beta=10$) shows no signs
of such oscillations (see (c)).

Given the quadratic nature of the Hamiltonians, our study also adds to the
emerging area of identifying signatures of one-body chaos. This has recently
been explored in SYK-2Winer  _et al._ (2020); Liao  _et al._ (2020) as well as
in strongly coupled free gauge theories Chen (2022). To find such signatures
of one-body chaos and its interplay with topological order, we study a variant
of SSH model where a subregion of the bulk is randomized with all-to-all
random hoppings while keeping the edge state structure protected. We find that
the SFF shows a dip followed by an early time oscillating exponential ramp
reminiscent of one-body chaotic signatures in SYK-2 Winer  _et al._ (2020).
This develops into an intermediate linear ramp, that plateaus at late times in
the trivial phase. In the topological phase a different picture emerges. While
in any representative configuration, at late times there are Rabi
oscillations; under ensemble averaging, these oscillations get destroyed given
the random phase lags between various Rabi modes. This reflects the non-self
averaging characteristic of the SFF Prange (1997) in the topological, yet
disordered systems. Interestingly this averaged asymptotic value in the
topological phase is different from that in the trivial phase. We show that it
depends on the random matrix properties of the zero energy manifold. We end
the paper with a perspective on how this physics of interplay between
topological features and chaotic signatures may be a rich playground to
uncover a host of new phenomena in both lattice and field theoretic quantum
many body systems.

Model and SFF: Given a set of single particle eigenvalues $\epsilon_{n}$ of a
quantum non-interacting fermionic system, the generalized partition function
is given by

$Z(\beta+it)=\prod_{n}\Big{(}1+\exp(-(\beta+it)\epsilon_{n})\Big{)}$ (1)

where $n\in{0,\ldots,L-1}$, $L$ is the system size, $\beta$ is inverse
temperature and $t$ is the real time. SFF is defined as

$\text{SFF}(\beta,t)\equiv{\mathcal{Z}}_{2}(\beta,t)=\frac{Z(\beta+it)Z(\beta-
it)}{Z(\beta)^{2}}$ (2)

Therefore for any non-interacting system the SFF can be straightforwardly
evaluated once the single particle energies are known. In particular for the
SSH model we study, the Hamiltonian is given by:

$H=\sum_{i}-v(c^{\dagger}_{iA}c^{\vphantom{\dagger}}_{iB}+\text{h.c.})-w(c^{\dagger}_{iB}c^{\vphantom{\dagger}}_{i+1,A}+\text{h.c.})$
(3)

where $c^{\dagger}_{i\alpha},c^{\vphantom{\dagger}}_{i,\alpha}$ are the
fermion creation and annihilation operators at site $i$ for the orbitals
$\alpha\equiv A,B$. $v$ and $w$ are the intra and inter-unit cell hopping
strength. The system hosts a topological phase when $|v/w|<1$, characterized
by a non-trivial winding number and has two close to zero energy modes in an
open system with a edge localization length $\xi=[\log(|w/v|)]^{-1}$ which
reflects extent of the edge states. The system has time-reversal and sub-
lattice symmetry restricting it to BDI symmetry class which realizes general
off-diagonal real matrices of the free fermion ten-fold classification Altland
and Zirnbauer (1997); Agarwala  _et al._ (2017), such that the system spectra
is always symmetric about energy $E=0$.

For a eigenspectrum which is symmetric about $E=0$ the SFF (see eqn. (2)) can
be written as

$\text{SFF}=\prod_{\epsilon_{n}>0}Z_{2}^{\epsilon_{n}}(\beta,t)=\prod_{\epsilon_{n}>0}\frac{\left(\cosh(\beta\epsilon_{n})+\cos(\epsilon_{n}t)\right)^{2}}{\left(1+\cosh(\beta\epsilon_{n})\right)^{2}}$
(4)

which implies that the real-time behavior of the system takes a highly
convoluted form dependent on the frequencies of the single-particle energies.
This is rather uninteresting as can be seen by the behavior of SFF for $v/w=2$
(trivial region) (See Fig. 1(c)). However an interesting behavior emerges when
boundary modes dominate, whose energy $\epsilon_{1}\rightarrow 0$
exponentially in system size. This results in a time scale
$\sim\pi/\epsilon_{1}$ where the SFF first goes to zero and then oscillates
with the same frequency, drowning all noisy behavior due to higher energy
modes which are exponentially killed due to a finite $\beta$ (see SFF in Fig.
1(b) for $v/w=0.5$ in the topological regime).

Rabi Oscillations: We now discuss that the long time oscillations of the SFF
are in fact generalized Rabi oscillations where the system locks between the
boundary zero modes. To see this consider a general many body wavefunction
which is an equal superposition of all the many body basis states

$|\Psi\rangle=\frac{1}{2^{L}}\sum_{\\{N_{n}\\}}|\\{N_{n}\\}\rangle$ (5)

where each of the $|\\{N_{n}\\}\rangle$ specify the Fock state representation
labelled by occupancies ($N_{n}=0,1$) of the single particle state
$|\psi_{n}\rangle$ with single particle eigen-energy $\epsilon_{n}$. When this
state is quenched with the Hamiltonian and and compared with the initial
state, this gives a fidelity of

$\langle\Psi|\Psi(t)\rangle\propto Z(it)$ (6)

Thus instances where $Z(t)\rightarrow 0$ are Rabi oscillations of a pure state
under time evolution. At $\beta=0$ the fidelity behavior will be
uncharacteristic because all $\epsilon_{n}$s will show convoluted
oscillations. However the interpretation is distinct at a finite temperature
where the initial state (see eqn. (5)) can be generalized to
$\propto|\Psi(\beta,t)\rangle=\sum_{\\{N_{n}\\}}\exp(-\beta
E_{\\{N_{n}\\}}/2)|\\{N_{n}\\}\rangle$. When evolved in time, the
corresponding fidelity is $Z(\beta+it)/Z(\beta)$ and therefore the SFF is just
its absolute value squared. Hence the SFF going to zero are the Rabi
oscillations of the state $|\Psi(\beta,t)\rangle$. It is insightful to note
the behaviour of this state at finite temperature. Given the states are
normalized by the thermal occupancy factors, the many body states with largest
weights are the ones where $N_{n}=1\leavevmode\nobreak\
\forall\leavevmode\nobreak\ \epsilon_{n}<0$. In a periodic system, this
isolates a single state; however, in an open one dimensional SSH chain, when
in topological phase – this results in four states which are exponentially
close in their many body energies. These essentially correspond to different
ways of occupying the (right and left) boundary modes
($\equiv|R\rangle,|L\rangle$). Thus

$\displaystyle|\Psi(\beta,t)\rangle$ $\displaystyle\sim$
$\displaystyle\frac{1}{2}\Big{(}|\circ\circ\rangle+e^{-(\frac{\beta}{2}+it)\epsilon_{a}}|*\circ\rangle$
$\displaystyle+$ $\displaystyle
e^{-(\frac{\beta}{2}+it)\epsilon_{b}}|\circ*\rangle+|**\rangle\Big{)}\otimes|N_{n}=1\leavevmode\nobreak\
\forall\leavevmode\nobreak\ \epsilon_{n}<0\rangle$

where $\epsilon_{a},\epsilon_{b}$ ($\epsilon_{b}=-\epsilon_{a}$) represent the
anti-bonding and bonding orbital combined out of the left and right edge
states $(\\{|b\rangle,|a\rangle\\}=\frac{1}{\sqrt{2}}(|R\rangle\pm|L\rangle)$
and $*(\circ)$ represents the their occupancies (vacancies) in
$|n_{a},n_{b}\rangle$ basis. Other states get damped by the bulk gap between
the valence and conduction bands $\sim\exp(-\beta\Delta_{g})$ where
$\Delta_{g}\sim|w-v|$. This introduces a temperature scale
$T^{*}\sim\frac{1}{\beta^{*}}=\Delta_{g}$, below which the Rabi oscillations
are strong, while for $T>\Delta_{g}$ these Rabi oscillations dissolve with the
bulk signatures. As is clear these Rabi oscillations have a frequency
determined by the edge mode energy scale given by
$\sim\frac{\pi}{\epsilon_{a}}$. Unlike a single qubit Rabi oscillation here
the SFF behaves as $\sim|(1+\cos(\epsilon_{a}t))|^{2}$, thus the rise from the
minima is $\propto t^{4}$ (see Fig. 1). In Fig. 2 we show the behavior of the
SFF both in the trivial, topological and critical point; showing that the
topological phase (b) has dominant oscillations. The behavior in $\\{v/w,t\\}$
plane shows that the time period of Rabi oscillations are in fact dependent on
$\epsilon_{a}$ as it matches the dashed analytical curves expressed as:
$t_{\texttt{Rabi}}=2.653(2k+1)\pi\exp(-0.483L\log(v/w))$. Here, $k$ are
positive integers and $L$ denotes the system size. (see Fig. 2(d))

Figure 2: Contour plot of SFF in the $\beta-t$ plane for an open boundary SSH
model (L = 60) in trivial regime (a) with $v=1,w=0.5$, in topological regime
(b) with $v=0.5,w=1$ and at the critical point (c) with $v=0.5,w=0.5$. (d)
Behavior of SFF in the $v/w-t$ plane for $\beta=30$. The dashed lines are
analytical curves denoting the zeros of the SFF (see text).

Higher Order Topological Phases: In order to further study this physics, and
its generalizations, we study SFF in a higher order topological phase where a
two dimensional system hosts topologically protected four corner modes on a
square lattice Benalcazar  _et al._ (2017). The Hamiltonian is given by
$H(k)=[\gamma+\lambda\cos(k_{x})]\Gamma_{4}+\lambda\sin(k_{x})\Gamma_{3}+[\gamma+\lambda\cos(k_{y})]\Gamma_{2}+\lambda\sin(k_{y})\Gamma_{1}+\delta\Gamma_{0}$
where $\Gamma_{0}=\tau_{3}\sigma_{0},\Gamma_{k}=-\tau_{2}\sigma_{k}$ and
$\Gamma_{4}=\tau_{1}\sigma_{0}$; where $\tau,\sigma$ are Pauli matrices. The
system hosts four corner modes when $|\gamma/\lambda|<1$. Evaluating the SFF,
we again find long time oscillations (see Fig. 3) where now the effective
state spans over $2^{4}$ states which can be counted as occupancies of four
boundary modes. The oscillations follow
$\sim|(3+4\cos(\epsilon_{a}t)+\cos(2\epsilon_{a}t))|^{2}$ where $\epsilon_{a}$
is the exponentially small energy scale close to zero. Interestingly the rise
form zero is now $\sim t^{8}$ reflecting the higher number of zero modes in
the system (see Fig. 3). In fact for a multipole topological
insulatorBenalcazar  _et al._ (2017) with general number of $2p$ zero modes,
SFF would the scale $\propto t^{4p}$, in the topological phase for
$T<\Delta_{g}$. This is one of the key results of this work.

While the discussion above has been in terms of single particle eigenstates,
in general, SFFs in many body settings is easily decomposed into a form

$\text{SFF}=\frac{1}{Z^{2}}\sum_{m,n}e^{-\beta(E_{m}+E_{n})}e^{it(E_{m}-E_{n})}$
(8)

where $E_{m},E_{n}$ are many-body eigenvalues; it may then appear that that
the time averaged behavior must approach $Z(2\beta)/Z^{2}$ on general grounds.
Our results point out that even in an otherwise dense spectra, existence of
topologically protected zero modes may render this asymptotic value
insignificant and mask it with generalized Rabi-oscillations. It is natural to
pose what is the fate of such Rabi oscillations in presence of disorder, and
that is the question we next turn to.

Figure 3: (a) One-particle energy spectrum of 2D HOTI model for parameter
values $\delta=0,\leavevmode\nobreak\ \gamma=0.5$ and $\lambda=1$ with open
boundary condition on a $12\times 12$ lattice. In the inset the four
topological zero energy modes are shown. (b) Local probability density of the
zero modes showing that they are localized at the corners of the lattice (c)
Plot of SFF in the trivial regime ($\delta=0,\gamma=1,\lambda=0.5$) of the
HOTI model, showing no distinct Rabi oscillations. (d) Generalized Rabi
oscillations in the topological regime for the HOTI system (blue)
($\delta=0,\leavevmode\nobreak\ \gamma=0.5,\lambda=1$) and SSH system (red)
($v=0.5154,w=1,L=24$). The rise of the SFF in the former is $\sim t^{8}$
compared to the latter.

Random Bulk Hamiltonians and topological zero modes: Motivated by random
matrix theory (RMT), where the symmetry properties of random dense Hamiltonian
matrices determine their chaotic signatures Dyson (1962); Bertini  _et al._
(2018); Chan  _et al._ (2018), we mark a finite central region (excluding
boundaries) in the bulk of SSH system (see eqn. (3)) as $\equiv{\cal R}$,
where all-to-all hopping disorder of the form
$\sum_{i,j}w_{ij}\,c^{\dagger}_{iA}c^{\vphantom{\dagger}}_{jB}+\text{h.c.}$ is
introduced. Here $w_{ij}$’s are chosen from a Gaussian orthogonal ensemble
with a scale parameter $\frac{\sigma}{\sqrt{N_{\cal R}}}$ where $N_{\cal R}$
are the number of sites in the region ${\cal R}$ and $\\{i,j\\}\in{\cal R}$.
Since the disorder respects the sublattice, ${\cal R}$ excludes the edges, and
$\sigma<|w-v|$, hence every disorder configuration retains topologically
protected zero modes (for more details on stability to disorder, see
Supplemental Material(SM)).

We find that the disorder averaged SFF, irrespective of the topological
character of the phase, has the following behavior as a function of time: (i)
a dip (ii) a exponential oscillation and then a linear ramp (iii) a late time
saturation. The three regions are shown in Fig. 4(a). At short times, the
exponentially oscillating behavior is characteristic of signatures of single
particle chaos in systems such as SYK-2 Winer  _et al._ (2020); Liao  _et al._
(2020) or in free strongly coupled gauge/string theories Chen (2022). This is
in contrast to many body chaos that has a distinct linear ramp right after the
dip Cotler  _et al._ (2017b); Liu (2018). To understand the underlying
physics, it is instructive to consider a further simplified two band model
where the single particle density of states comprises of two semicircular
lobes ($\equiv\rho_{\pm}(E)$) each with a bandwidth of $2a$ but with a band
gap of $\Delta_{g}$ (details in SM). Assuming that the eigenvalue
distributions contributing to any lobe (say $\rho_{+}(E)$) are governed by
$N\times N$ random matrix drawn from the gaussian unitary ensemble the SFF is
given by (see eqn. (4))

${\left\langle{\mathcal{Z}}_{2}^{R}(\beta,t)\right\rangle}=\bigg{\langle}\exp\left(N\int
dE\,\rho_{+}(E)\log Z_{2}^{E}(\beta,t)\right)\bigg{\rangle}$ (9)

The angular brackets indicate the RMT averaging where all moments of the
density of states $\rho_{+}(E)$ contribute. At early times, the disconnected
piece dominates leading to $\sim\exp\big{[}\left(N\int
dE\,\langle\rho_{+}(E)\rangle\log Z_{2}^{E}(\beta,t)\right)\big{]}$. Under a
high temperature expansion the system exhibits the exponential ramp imbued
with short-time oscillations. The $\beta\rightarrow 0$ limit is given by

$\displaystyle{\left\langle{\mathcal{Z}}_{2}^{R}(0,t)\right\rangle}\approx\frac{1}{16^{N}}\exp\left(8N\sum_{k=1}^{\infty}\frac{(-1)^{k+1}J_{1}(kat)\cos(ka_{\texttt{max}}t)}{k^{2}at}\right)$
(10)

where, $a_{\texttt{max}}=\Delta+a$, and in a microscopic model is decided by
the energy scale where the bulk DOS peaks. This result points out that the
location of the dip is set by high energy scale $\pi/a_{\texttt{max}}$, which
is both system size independent and is unaffected by its topological
properties (see Fig. 4). The initial exponential oscillations next lead to a
linear ramp, which occurs at $t\sim L^{2/3}$. This linear ramp is a result of
the leading connected piece i.e., $\rho_{c}^{(2)}(E_{1},E_{2})$, just as in
the case of many-body SFF Cotler  _et al._ (2017c). The linear ramp further
leads to a saturation, which can be remarkably different if the underlying
system is trivial or topological. More interestingly the saturation physics is
temperature dependent, as we now discuss.

At high temperatures ($T>|w-v|$), irrespective of the topological features SFF
saturates to a value dominated by the Kubo gap of the bulk spectra. This is
essentially where all connected spectral correlations are featureless. For the
minimal model, this predicts emergence of a plateau at a time $t\sim 4N\pi$
(see SM) and $\sim L$ for bulk disordered SSH system. However now as $T$ is
reduced below the bulk gap scale ($T<|w-v|$) (see Fig. 1), depending on
whether we are in the trivial or the topological phase another saturation
plateau appears.

Figure 4: (a) SFF in random SSH model ($L=60,\beta=0,\sigma=0.01$, 30 random
sites in the bulk) shows the ramp structure in early time for both topological
($v=0.5,w=1$) and trivial regimes ($v=1,w=0.5$). The green and orange dashed
vertical lines denote the dip-time and plateau-time calculated from the toy
model. (b) The same system at a higher $\beta$ ($=10$) shows a step-like
plateau in late times in topological regime, while the trivial SFF shows no
such feature. (c) At very high $\beta$ ($=50$), the early time ramp
disappears; the late-time plateau in the topological SFF saturates around
$3/8$. Here, $\Delta_{g}^{-1}\sim 2$ and $\epsilon_{0}^{-1}\sim 10^{7}$. All
the SFF values are averaged over 2000 random configurations. For further
details, see text.

The difference between the trivial and the topological phases of the
disordered model at low temperatures is due to the presence of exponentially
close to zero energy modes ($\sim\pm\epsilon_{0}$) which as we know from the
clean limit leads to large time Rabi oscillations (see Fig. 1). The effective
SFF can be captured in the form $\sim
Z_{2}^{\epsilon_{0}}(\beta,t){\mathcal{Z}}_{2}^{R}(\beta,t)$ (see eqn. (4))
where at low temperatures $Z_{2}^{\epsilon_{0}}$ dominates and thus any single
configuration, even with bulk disorder, will lead to large time oscillations
after a time scale $t^{*}\sim\frac{1}{\epsilon_{0}}$. However, when
configuration averaged, SFF shows an intriguing behavior as is shown in Fig. 4
(b) and (c). At low temperatures the system reaches a different saturation
value at later times. The central result of our next discussion is to show
that this new saturation value is not governed by large-$N$ RMT results, as is
common in chaotic systems, but is rather governed by a small-$N$ RMT chaos, as
we discuss below.

Every peripheral disorder realization in region ${\cal R}$ perturbs the edge
modes which are in turn protected by the $(w,v)$ scale of the underlying SSH
Hamiltonian. Therefore at low energies i.e. finite $\beta$, i.e.,
$\beta>\Delta_{g}^{-1}$, the SFF gets dominated by RMTs reflecting the
couplings within the zero energy manifold itself. For instance the effective
Hamiltonian in context of our example becomes
$H_{\text{eff}}=\epsilon_{0}\sigma_{z}+\lambda\sigma_{x}$, where $\lambda$ can
be effectively considered to be drawn from probability distribution
$P(\lambda)$. Thus the configuration averaged $Z_{2}^{\epsilon_{0}}(\beta,t)$
becomes

$\displaystyle\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z_{2}^{\epsilon_{0}}(\beta,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}$
$\displaystyle=\int_{-\infty}^{\infty}d\lambda\,P(\lambda)\frac{\left(\cosh(\beta{\epsilon}_{\lambda})+\cos({\epsilon}_{\lambda}t)\right)^{2}}{\left(1+\cosh(\beta{\epsilon}_{\lambda})\right)^{2}}$
(11)

where $\epsilon_{\lambda}=\sqrt{\epsilon_{0}^{2}+\lambda^{2}}$. For a Gaussian
distribution of $\epsilon_{\lambda}$ around $\epsilon_{0}$ with variance
$\sigma$,
$\lim_{t\rightarrow\infty}\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z_{2}^{\epsilon_{0}}(\beta,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}=\alpha_{\text{topo}}$
is $\sim 3/8$ when $\epsilon_{o}\ll T<|w-v|$. Thus this late time behavior in
the topological phase is fundamentally distinct from the trivial regime where
no such zero energy manifold exists and the SFF just saturates to $1$ (see
SM). The analytical predictions exactly match our numerical results as shown
in Fig. 4 (b) and (c). Thus SFF and its utility in identifying chaotic
signatures, are fundamentally dependent on the topological properties of a
quantum system.

Outlook: In this work we have shown that interplay of topological order has
important implications for both dynamics and for emergence of chaos in quantum
systems. Our analysis establishes that emergence of zero energy boundary modes
fundamentally change the late time behavior of SFF, a versatile tool that has
been critical to diagnose chaos and thermalization in host of systems. In
particular we have shown that systems with close to zero boundary modes can
exhibit generalized Rabi oscillations which are characterized by topological
properties. Our mapping of SFF to the fidelity of a generalized state makes
this particularly transparent. We further analysed the role of disorder in the
system and showed, using toy models, that the late time behavior of the system
in fact holds crucial information of the topological order. While our work has
investigated the role of non-interacting topology on the SFF, the natural
question to pose is the role of both interactions and long-range topological
order. Given the existence of a gap, we expect our results to hold
perturbatively under interactions, a systematic study of which, is an exciting
future prospect.

Acknowledgement. A.S and S.P acknowledge funding from IIT Kanpur. A.A.
acknowledges support from IIT Kanpur Initiation Grant IITK/PHY/2022010. D.D.,
& A.S. would like to acknowledge the support provided by the Max Planck
Partner Group grant MAXPLA/PHY/2018577. D.D. would also like to acknowledge
the support provided by the MATRICS grant SERB/PHY/2020334.

## References

  * Srednicki (1994) Mark Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
  * D’Alessio  _et al._ (2016) Luca D’Alessio, Yariv Kafri, Anatoli Polkovnikov, and Marcos Rigol, “From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics,” Advances in Physics 65, 239–362 (2016), https://doi.org/10.1080/00018732.2016.1198134 .
  * Borgonovi  _et al._ (2016) F. Borgonovi, F.M. Izrailev, L.F. Santos, and V.G. Zelevinsky, “Quantum chaos and thermalization in isolated systems of interacting particles,” Physics Reports 626, 1–58 (2016), quantum chaos and thermalization in isolated systems of interacting particles.
  * Das  _et al._ (2018) Avijit Das, Saurish Chakrabarty, Abhishek Dhar, Anupam Kundu, David A. Huse, Roderich Moessner, Samriddhi Sankar Ray, and Subhro Bhattacharjee, “Light-cone spreading of perturbations and the butterfly effect in a classical spin chain,” Phys. Rev. Lett. 121, 024101 (2018).
  * Pollack  _et al._ (2020) Jason Pollack, Moshe Rozali, James Sully, and David Wakeham, “Eigenstate thermalization and disorder averaging in gravity,” Phys. Rev. Lett. 125, 021601 (2020).
  * Murugan  _et al._ (2021) Sugan Durai Murugan, Dheeraj Kumar, Subhro Bhattacharjee, and Samriddhi Sankar Ray, “Many-body chaos in thermalized fluids,” Phys. Rev. Lett. 127, 124501 (2021).
  * del Campo  _et al._ (2017) A. del Campo, J. Molina-Vilaplana, and J. Sonner, “Scrambling the spectral form factor: unitarity constraints and exact results,” Phys. Rev. D 95, 126008 (2017), arXiv:1702.04350 [hep-th] .
  * Bena  _et al._ (2005) I Bena, M Droz, and A Lipkowski, “Statistical mechanics of equilibrium and nonequilibrium phase transitions: The yang-lee formalism,” International Journal of Modern Physics B 19, 4269–4329 (2005).
  * Harrow  _et al._ (2020) Aram W. Harrow, Saeed Mehraban, and Mehdi Soleimanifar, “Classical algorithms, correlation decay, and complex zeros of partition functions of quantum many-body systems,” Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing (2020), 10.1145/3357713.3384322.
  * Cotler  _et al._ (2017a) Jordan S. Cotler, Guy Gur-Ari, Masanori Hanada, Joseph Polchinski, Phil Saad, Stephen H. Shenker, Douglas Stanford, Alexandre Streicher, and Masaki Tezuka, “Black Holes and Random Matrices,” JHEP 05, 118 (2017a), [Erratum: JHEP 09, 002 (2018)], arXiv:1611.04650 [hep-th] .
  * Saad  _et al._ (2018) Phil Saad, Stephen H. Shenker, and Douglas Stanford, “A semiclassical ramp in SYK and in gravity,” (2018), arXiv:1806.06840 [hep-th] .
  * Dyer and Gur-Ari (2017) Ethan Dyer and Guy Gur-Ari, “2D CFT Partition Functions at Late Times,” JHEP 08, 075 (2017), arXiv:1611.04592 [hep-th] .
  * Cotler  _et al._ (2017b) Jordan Cotler, Nicholas Hunter-Jones, Junyu Liu, and Beni Yoshida, “Chaos, Complexity, and Random Matrices,” JHEP 11, 048 (2017b), arXiv:1706.05400 [hep-th] .
  * Liu (2018) Junyu Liu, “Spectral form factors and late time quantum chaos,” Phys. Rev. D 98, 086026 (2018), arXiv:1806.05316 [hep-th] .
  * Ludwig (2016) Andreas W W Ludwig, “Topological phases: classification of topological insulators and superconductors of non-interacting fermions, and beyond,” Physica Scripta 2016, 014001 (2016).
  * Chiu  _et al._ (2016) Ching-Kai Chiu, Jeffrey C. Y. Teo, Andreas P. Schnyder, and Shinsei Ryu, “Classification of topological quantum matter with symmetries,” Rev. Mod. Phys. 88, 035005 (2016).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium : Topological insulators,” Rev. Mod. Phys. 82, 3045–3067 (2010).
  * Qi and Zhang (2011) Xiao-Liang Qi and Shou-Cheng Zhang, “Topological insulators and superconductors,” Rev. Mod. Phys. 83, 1057–1110 (2011).
  * Su  _et al._ (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, “Solitons in polyacetylene,” Phys. Rev. Lett. 42, 1698–1701 (1979).
  * Benalcazar  _et al._ (2017) Wladimir A Benalcazar, B Andrei Bernevig, and Taylor L Hughes, “Quantized electric multipole insulators,” Science 357, 61–66 (2017).
  * Winer  _et al._ (2020) Michael Winer, Shao-Kai Jian, and Brian Swingle, “An exponential ramp in the quadratic Sachdev-Ye-Kitaev model,” Phys. Rev. Lett. 125, 250602 (2020), arXiv:2006.15152 [cond-mat.stat-mech] .
  * Liao  _et al._ (2020) Yunxiang Liao, Amit Vikram, and Victor Galitski, “Many-body level statistics of single-particle quantum chaos,” Phys. Rev. Lett. 125, 250601 (2020), arXiv:2005.08991 [cond-mat.stat-mech] .
  * Chen (2022) Yiming Chen, “Spectral form factor for free large N gauge theory and strings,” JHEP 06, 137 (2022), arXiv:2202.04741 [hep-th] .
  * Prange (1997) R. E. Prange, “The spectral form factor is not self-averaging,” Phys. Rev. Lett. 78, 2280–2283 (1997).
  * Altland and Zirnbauer (1997) Alexander Altland and Martin R. Zirnbauer, “Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures,” Phys. Rev. B 55, 1142–1161 (1997).
  * Agarwala  _et al._ (2017) Adhip Agarwala, Arijit Haldar, and Vijay B. Shenoy, “The tenfold way redux: Fermionic systems with n-body interactions,” Annals of Physics 385, 469 – 511 (2017).
  * Dyson (1962) F. J. Dyson, “Statistical theory of the energy levels of complex systems. I,” J. Math. Phys. 3, 140–156 (1962).
  * Bertini  _et al._ (2018) Bruno Bertini, Pavel Kos, and Toma ž Prosen, “Exact spectral form factor in a minimal model of many-body quantum chaos,” Phys. Rev. Lett. 121, 264101 (2018).
  * Chan  _et al._ (2018) Amos Chan, Andrea De Luca, and J. T. Chalker, “Spectral statistics in spatially extended chaotic quantum many-body systems,” Phys. Rev. Lett. 121, 060601 (2018).
  * Cotler  _et al._ (2017c) Jordan Cotler, Nicholas Hunter-Jones, Junyu Liu, and Beni Yoshida, “Chaos, Complexity, and Random Matrices,” JHEP 11, 048 (2017c), arXiv:1706.05400 [hep-th] .
  * Brezin and Zee (1993) E. Brezin and A. Zee, “Universality of the correlations between eigenvalues of large random matrices,” Nucl. Phys. B 402, 613–627 (1993).
  * Resta (1998) Raffaele Resta, “Quantum-mechanical position operator in extended systems,” Phys. Rev. Lett. 80, 1800–1803 (1998).

Supplemental Material to
Spectral Form Factors of Topological Phases

##  I Toy model

Here we first discuss the toy model which effectively captures the physics of
SSH model with bulk disorder. The toy model consists of $N$ one-particle
energy states taken from a GUE of random matrices which has semicircular
$\rho(E)$ of radius $a$, centered around $a_{\texttt{max}}=a+\Delta$ (with
$\Delta>0$) and a corresponding negative $E$ copy (see Fig. 5$(a)$). The
positive energy branch of density of states $\rho_{+}(E)$ is distributed
between $a_{1}=a_{\texttt{max}}-a$ to $a_{2}=a_{\texttt{max}}+a$ with a
spectral gap of $\Delta$. The ensemble averaged $\rho_{+}(E)$ is given by

$\displaystyle{\left\langle\rho_{+}(E)\right\rangle}=\frac{2}{a^{2}\pi}\sqrt{(E-a_{1})(a_{2}-E)}$
(12)

Figure 5: (a) Toy model density of states. (b) Toy model ensemble averaged SFF
for $a_{\texttt{max}}=1.07,\,a=0.6,\,N=40$. The green, red and magenta lines
denote $t_{\texttt{dip}},\leavevmode\nobreak\
t_{\texttt{crossover}},\leavevmode\nobreak\ t_{\texttt{plateau}}$
respectively.

The many body SFF when ensemble averaged is given by ${\left\langle
Z^{R}_{2}(\beta,t)\right\rangle}$ which is

$\begin{split}&{\left\langle{\mathcal{Z}}^{R}_{2}(\beta,t)\right\rangle}={\left\langle\exp\left(N\int\text{d}E\,\rho_{+}(E)\,\log
Z^{E}_{2}\right)\right\rangle}=\sum^{\infty}_{n=0}\frac{1}{n!}N^{n}\int\left(\prod_{i=1}^{n}\text{d}E_{i}\right){\left\langle\rho_{+}(E_{1})\ldots\rho_{+}(E_{n})\right\rangle}\prod_{j=1}^{n}\log
Z^{E_{j}}(\beta,t)\end{split}$ (13)

where a typical $n$-point density correlator has the form

$\begin{split}&{\left\langle\rho_{+}(E_{1})\ldots\rho_{+}(E_{n})\right\rangle}\equiv{\left\langle\rho^{(n)}\right\rangle}=\prod_{i=1}^{n}{\left\langle\rho^{(1)}_{i}\right\rangle}+\sum_{\\{n_{i},m\\}}A^{\\{n_{i},m\\}}\prod_{i=1}^{k}\rho_{c}^{(n_{i})}\prod_{j}^{m}{\left\langle\rho^{(1)}_{j}\right\rangle}\text{
with: }\sum_{\begin{subarray}{c}i=1\\\
n_{i}>1\end{subarray}}^{k}n_{i}+m=n;\leavevmode\nobreak\ \sum_{i}^{k}n_{i}\neq
0\end{split}$ (14)

Here, the first term denotes the disconnected piece, $A^{\\{n_{i},m\\}}$
stands for the coefficient arising from the permutations of the indices,
$\rho_{c}^{(n_{i})}$ denotes the completely connected piece of the
$n_{i}$-point correlation and $\rho^{(1)}_{j}\equiv\rho_{+}(E_{j})$. Using the
explicit form for the joint distributions, as given through determinant of the
kernel Brezin and Zee (1993) we find :

$\displaystyle\prod_{i=1}^{n}{\left\langle\rho^{(1)}_{i}\right\rangle}$
$\displaystyle={\cal
O}(1),\,\,\,\,\prod_{i=1}^{k}\rho_{c}^{(n_{i})}\prod_{j}^{m}{\left\langle\rho^{(1)}_{j}\right\rangle}={\cal
O}\left(N^{m-n}\right)={\cal O}\left(N^{-\sum_{i}n_{i}}\right)$ (15)

Since $\text{min}(n_{i})\geq 2$ the second term is always at least ${\cal
O}(1/N^{2})$ suppressed compared to the fully disconnected one. Of course,
after evaluating the integrals in (13) there are non-trivial time growths
associated with various pieces. At early times, $t\ll N$, we can ignore these
time dependences and the dominant contribution at large $N$, comes from the
fully disconnected DOS correlator. Therefore the sum over $n$ in (13) re-
exponentiates and we find the early time, ensemble averaged SFF to be:

$\begin{split}&{\left\langle{\mathcal{Z}}^{R}_{2}(\beta,t)\right\rangle}\simeq\exp\left(N\int\text{d}E\,{\left\langle\rho_{+}(E)\right\rangle}\,\log
Z^{E}_{2}\right)=\exp\left(N\int_{a_{1}}^{a_{2}}\text{d}E\,\frac{2}{a^{2}\pi}\sqrt{(E-a_{1})(a_{2}-E)}\,\log
Z^{E}_{2}\right)\end{split}$ (16)

At high temperatures, $\beta=0$, the integral can be evaluated exactly in
terms of Bessel function:

${\left\langle{\mathcal{Z}}^{R}_{2}(0,t)\right\rangle}\simeq\frac{1}{16^{N}}\exp\left(8N\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k^{2}at}J_{1}(kat)\cos(ka_{\texttt{max}}t)\right)$
(17)

In the above expression, the higher frequency oscillations
$\cos(ka_{\texttt{max}}t)$ are suppressed as $k^{-2}$. Thus the $k=1$
component dominates and produces the dip at the first oscillation minima:
$t_{\texttt{dip}}=\pi/a_{\texttt{max}}$.

To see the dynamics at an intermediate time, we look at the connected
component of ${\left\langle\rho_{+}(E_{1})\rho_{+}(E_{2})\right\rangle}$
denoted as $\rho_{c}^{(2)}(E_{1},E_{2})$. In terms of sine kernel Cotler  _et
al._ (2017b), from (13) one obtains:

$\begin{split}&\frac{1}{2!}N^{2}\int\text{d}E_{1}\text{d}E_{2}\rho_{c}^{(2)}(E_{1},E_{2})\prod_{j=1}^{2}\log
Z^{E_{j}}_{2}=-\frac{N^{2}}{2!}\int\text{d}E_{1}\text{d}E_{2}\frac{\sin^{2}\left(N\pi(E_{1}-E_{2})\right)}{N^{2}\pi^{2}(E_{1}-E_{2})^{2}}\prod_{j=1}^{2}\log
Z^{E_{j}}_{2}\\\
&=-\frac{2}{\pi}\left(N\pi\log^{2}(4)-\log(4)\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}\left(4N\pi-
kt\right)\cos(ka_{\texttt{max}}t)\Theta\left(\frac{4N\pi}{k}-t\right)+\right.\\\
&\left.\sum_{k,l=1}^{\infty}\frac{(-1)^{k+l}}{kl}\left[\left(4N\pi-(k+l)t\right)\cos(|k-l|a_{\texttt{max}}t)\Theta\left(\frac{4N\pi}{(k+l)}-t\right)+\left(4N\pi-|k-l|t\right)\cos((k+l)a_{\texttt{max}}t)\Theta\left(\frac{4N\pi}{|k-l|}-t\right)\right]\right)\end{split}$
(18)

We see the 2-point connected piece generates linear terms of $\left(4N\pi-
kt\right)$ upto $t=4N\pi/k$. Thus the $k=1$ component gives the longest
surviving linear piece with the plateau time $t_{\texttt{plateau}}=4N\pi$ (see
Fig. 2$(b)$). Comparing the disconnected and connected pieces for $n=2$, one
can estimate the time $t_{\texttt{crossover}}$ at which the connected piece
starts to dominate over the disconnected piece. Using the asymptotic form of
Bessel functions one finds,

$N^{2}\frac{J_{1}(t_{\texttt{crossover}})}{t_{\texttt{crossover}}}\sim
N\leavevmode\nobreak\ \leavevmode\nobreak\ \Rightarrow\leavevmode\nobreak\
\leavevmode\nobreak\ t_{\texttt{crossover}}\sim\mathcal{O}(N^{2/3})$ (19)

So one sees a linear ramp starting from $t_{\texttt{crossover}}$ up to
$t_{\texttt{plateau}}$, and then the SFF plateaus.

As we see from the $n=2$ case, the surviving piece in long-time average is
$\sim\mathcal{O}(N)$ for the connected piece, while the disconnected piece has
a long-time averaged value of $\sim\mathcal{O}(N^{2})$. Thus in all $n$-point
correlations the disconnected piece gives the dominant long-time averaged
contribution. From (17), we see in long-time average the oscillations die
down, thus the plateau value becomes:

${\left\langle{\mathcal{Z}}^{R}_{2}(0,t\to\infty)\right\rangle}\sim e^{-N}$
(20)

##  II Late time SFF dynamics due to fluctuating zero modes:

In the 1D SSH model, there are two zero modes with energy $\pm\epsilon_{0}$ at
the two ends of the lattice chain. When disorder is introduced in the bulk,
the zero modes also fluctuate. To model the spectrum of these two fluctuating
zero modes we consider a Gaussian energy distribution centered around
$\epsilon_{0}$ with the variance $\sigma$. Thus the ensemble averaged zero
energy contribution to the SFF is:

$\displaystyle\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z^{\epsilon_{0}}_{2}(\beta,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}$
$\displaystyle=\int\text{d}\epsilon\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(\epsilon-\epsilon_{0})^{2}}{2\sigma^{2}}}\left(\frac{\cosh(\beta\epsilon)+\cos(\epsilon
t)}{\cosh(\beta\epsilon)+1}\right)^{2}$ (21)

For $\beta=0$ (i.e, $\beta\ll\epsilon_{0}^{-1}$) we have,

$\begin{split}\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z^{\epsilon_{0}}_{2}(0,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}&=\frac{1}{\sigma\sqrt{2\pi}}\int\text{d}\epsilon
e^{-\frac{(\epsilon-\epsilon_{0})^{2}}{2\sigma^{2}}}\left(\frac{1+\cos(\epsilon
t)}{2}\right)^{2}=\frac{3}{8}+\frac{1}{2}e^{-t^{2}\sigma^{2}/2}\cos(\epsilon_{0}t)+\frac{1}{8}e^{-2t^{2}\sigma^{2}}\cos(2\epsilon_{0}t)\end{split}$
(22)

At large $t$ the oscillations will average to zero, thus leading to a plateau
value of $3/8=0.375$.

To calculate the dip point, we can find the first minima of
$\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z^{\epsilon_{0}}_{2}(0,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}$ (ignoring the
$\cos(2\epsilon_{0}t)$) term and arrive at the transcendental equation:

$\begin{split}\partial_{t}\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z^{\epsilon_{0}}_{2}(0,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}=0\leavevmode\nobreak\
\Rightarrow\leavevmode\nobreak\
\tan(\epsilon_{0}t)=-\frac{t\sigma^{2}}{\epsilon_{0}}\end{split}$ (23)

When $\sigma\ll\epsilon_{0}$, in the above expression
$\tan(\epsilon_{0}t_{dip})\to 0$, thus:
$t_{\text{dip}}\to\frac{\pi}{\epsilon_{0}}$. In this case, the exponential
$\sim e^{-t^{2}\sigma^{2}/2}$ decay does not sufficiently suppress the
$\cos(\epsilon_{0}t)$ oscillations, thus one sees an oscillatory dip-ramp and
plateau. For $\sigma\gg\epsilon_{0}$, $\tan(\epsilon_{0}t_{dip})\to-\infty$,
hence: $t_{\text{dip}}\to\frac{\pi}{2\epsilon_{0}}.$ In this case, the
exponential decay dominates. Thus, one only sees a dip and plateau, without
any ramp. In the actual random SSH model, one has both the random excited
states and fluctuating zero modes, i.e.,
${\left\langle{\mathcal{Z}}^{\texttt{Full}}_{2}(\beta,t)\right\rangle}={\left\langle{\mathcal{Z}}^{R}_{2}(\beta,t)\right\rangle}\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z^{\epsilon_{0}}_{2}(\beta,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}$. Thus when
$\epsilon_{0}^{-1}\gg\beta>\Delta^{-1}$, the early-time dip-ramp due to
${\left\langle{\mathcal{Z}}^{R}_{2}(\beta,t)\right\rangle}$ starts to get
suppressed, and one starts to see the late time SFF dynamics of
$\mathopen{\hbox{\set@color${\langle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\langle}$}}Z^{\epsilon_{0}}_{2}(\beta,t)\mathclose{\hbox{\set@color${\rangle}$}\mkern
2.0mu\kern-3.49998pt\leavevmode\hbox{\set@color${\rangle}$}}$ (see Fig. 4 (b)
and (c) in the main text).

Figure 6: (a) The polarization $\mathcal{P}$ shows the stability of the
topological phase up to a critical $\sigma$ for the SSH Hamiltonian with
$90\%$ bulk region having all-to-all random hopping. The stability can also be
seen from the standard deviation (SD) of the polarization, which is clearly
related to the closing of the bulk gap ($\Delta E$) in the system. (b) Contour
plot of the bulk gap with the disorder strength $\sigma$ and the fraction of
random sites $N_{\mathcal{R}}/L$. For both the plots, the system size is
$L=100$, $v=0.5$, $w=1$ and all the data are averaged over $1000$ random
configurations.

##  III Stability to Disorder

In any topological phase presence of a bulk gap and underlying symmetries
protects the boundary modes under perturbative disorder. In our work, directed
towards studying the role of topological boundary modes in SFF, here we show
that even in presence of all-to-all random hopping (drawn from a Gaussian
distribution with zero mean and variance $\sigma$) in a finite region $\cal
R$, the boundary modes remain stable until a finite $\sigma$. We study
polarization, which is the the real space representation of the topological
invariant (winding number) in $1$D systems Resta (1998). The polarization is
defined as

$\mathcal{P}=\frac{1}{2\pi}\text{Im}\bigg{\\{}\text{Tr}\Big{[}\ln\big{(}\hat{P}D\hat{P}\big{)}\Big{]}\bigg{\\}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{mod $1$},$ (24)

where $\hat{P}=\sum_{E_{n}\leq E_{F}}|\psi_{n}\rangle{\langle\psi_{n}|}$ is
the ground state projection operator with $|\psi_{n}\rangle$ being the single-
particle eigenstate corresponding to the energy eigenvalue $E_{n}$ and $E_{F}$
is the Fermi energy of the system. The position ($x_{i}$) of all the lattice
sites (for system size $L$ there are $N=L/2$ unit-cells) are compactified and
exponentiated to give the operator $D=\text{diag}\big{[}e^{2\pi
x_{i\alpha}/N}\big{]}$, where $x_{i\alpha}=(i-1)$ for $i^{th}$ unit cell and
both $\alpha=A,B$ sublattices. In absence of disorder, in the topological
regime of SSH Hamiltonian ($|v/w|<1$), the polarization $\mathcal{P}$ is
quantized to $0.5$, while it is zero for the trivial insulating
phase($|v/w|>1$). Futhermore, in Fig. 6(a) we plot the polarization (see eqn.
24) of the system in the topological regime ($v/w=0.5$) as a function of bulk
disorder $\sigma$ in $90\%$ $(N_{\cal R}/L=0.9)$ of bulk sites. The system
shows quantized polarization as well as a finite bulk gap $\Delta E$ upto a
finite $\sigma$ confirming that indeed the topological properties of the
system remain stable even when bulk disorder randomizes the bulk spectra.
Futhermore the stability increases as the number of bulk sites disordered is
reduced; see Fig. 6(b) where the bulk gap is plotted as a function of $N_{\cal
R}/L$ and $\sigma$. We mention in passing that while our focus here has been
on properties of SFF - this model has a set of interesting ground state
properties which we will elaborate in future work.