\section{Introduction}
	
	The excitonic insulator (EI) is an unconventional phase of matter in which excitons condense into a macroscopically phase-coherent state \cite{Jerome1967,Kohn1967,Kohn1967mbp,Halperin1968}. 
	%
	This phase lies at the transition between semiconductors and semimetals \cite{Mott1968, Imada1998, Bronold2006}, and can be understood in terms of Bose-Einstein condensation (BEC) \cite{BECPitaevskiiStringari} or the Bardeen--Cooper--Schrieffer (BCS) \cite{BCS1957} theory of superconductivity \cite{Bronold2006,Phan2010,Seki2011,Perfetto2019}.
	%
	The EI is expected to exist in equilibrium and possibly at room temperature, and thus provide unprecedented access for the study of the BCS-BEC crossover, the metal-insulator transition and other strongly-correlated electron systems, as well as  potential applications such as thermoelectrics \cite{Wu2014} and dynamical band engineering \cite{Mor2017, Okazaki2018}. 
	
	Although it was first hypothesized decades ago, recent years have seen renewed interest in the EI, with the theoretical and experimental identification of several candidate materials, including transition metal \revised{chalcogenides} such as \ce{Ta2NiSe5}{}  \cite{Wakisaka2009, Kaneko2012, Kaneko2013,Seki2014, Kim2016,Sugimoto2016,Sugimoto2016a,Lu2017,Larkin2017, Mor2017, Larkin2018,Seo2018, Werdehausen2018ScienceAdvances,Okazaki2018, Domon2018,Lee2019}, \emph{1T}-\ce{TiSe2} \cite{DiSalvo1976,Pillo2000, Cercellier2007, Monney2009,Rossnagel2011,Hellmann2012, Kogar2017} and $\ce{TmSe_{1-x}Te_x}$ \cite{Wachter2004,Hulsen2006,Wachter2013}, as well as low-dimensional systems such as graphene, quantum wells, and quantum Hall bilayers \cite{Eisenstein2014,Zhu1995,Naveh1996,Khveshchenko2001, Eisenstein2004,Aleiner2007,Wu2015,Du2017,Varsano2017}.
	
	Despite this influx of candidate materials, confirmation of the existence of the EI phase in bulk materials has remained elusive. This is due to measured EI signatures being either hard to quantify, or difficult to separate from those of a normal insulator \cite{counterflow}. This situation signals an urgent need for a general distinguishing experimental signature which is unique to the EI.
	%
	
	
	To address this issue, we argue that a key distinguishing feature of the ideal EI is that it exhibits a spontaneously broken symmetry. 
	%
	The symmetry in question is {$U(1)$}, corresponding to the independent conservation of charge in the conduction and valence bands \cite{Phan2010,Wu2015}. While in practice this symmetry is not exact and different bands are generally coupled, it nevertheless can be approximately conserved if tunnelling between the bands is heavily suppressed. This could be due to a myriad of causes: the bands existing in two spatially separated planes (as in bilayer and quantum-well systems \cite{Eisenstein2004}); the localization of band orbitals around different atoms in the unit cell of \ce{Ta2NiSe5}{} \cite{Kaneko2013}; and the indirectness of the band gap in \emph{1T}-\ce{TiSe2} \cite{Rossnagel2011,Kogar2017}. Upon formation of the condensate, this symmetry spontaneously breaks. {(There is a separate $U(1)$ symmetry reflecting total charge conservation, which remains unbroken.)} 
	%
	Breaking of such a continuous symmetry gives rise to gapless Nambu--Goldstone collective modes (CMs) \cite{Nambu1960, Goldstone1961, Goldstone1962}. 
	%
	If the symmetry is not exact, the CMs will be gapped. Strong symmetry breaking could push the CM gap above the particle-hole continuum. However, if the symmetry violation is small, then the CM branch will remain low-lying. Therefore we suggest that the observation of soft modes, below the particle-hole continuum, is characteristic of an EI, providing vestigial evidence of the breaking of an underlying (approximate) symmetry.
	%
	Since such collective modes can be taken as evidence for an underlying EI phase, we 
	study the properties of these modes.  We focus on their robustness to disorder scattering, and outline their experimental consequences.
	
	
	We consider a minimal model which exhibits an excitonic phase transition: \revised{an extended Falicov--Kimball model \cite{Falicov1969} is driven into the EI phase by strong on-site electron--electron repulsion \cite{Batista2002, Batista2004, Ejima2014, Golez2016, Murakami2017,Tanabe2018, Murakami2020}. } The CMs are found via the density-density response function  $\chi\left(k,\omega\right)$ (defined below) within the random-phase approximation (RPA) \cite{Bohm1953, Pines1953, Chen2017}. We subsequently demonstrate the robustness of the CMs by calculating their mean free paths in the presence of a  disorder potential. Disorder can couple to the excitonic insulator in a variety of different channels, which we describe. \revised{In this we generalize Zittartz's previous study \cite{Zittartz1967} of EI quasiparticle disorder scattering and extend it to the collective modes.} We show that at small wavevectors the excitonic modes can have mean free paths orders of magnitude longer than those of single quasiparticles. We show that this can be  attributed to the relation between the CMs and the condensate broken symmetry, as well as their non-zero group velocity at small $k$.
	The large CM mean free paths suggest they should have an appreciable signature in thermal transport, and could potentially be observed directly in spatially-resolved ultrafast pump-probe experiments.
	
	
	The rest of this paper is organized as follows: In Sec.~\ref{sec:Model} we  present our model and the formalism of our diagrammatic linear response treatment, and we subsequently employ the RPA to identify the collective modes. We extend the formalism to incorporate disorder in Sec.~\ref{sec:Disorder}, and compute the disorder-induced scattering rates of both single particles and collective modes. In Sec.~\ref{sec:Discussion} we discuss our results qualitatively, formulate an effective field theory whereby we can interpret all of the results in a compact general fashion, and make connections to possible experimental consequences. Technical details of our calculations are given in the Appendices.
	
	
	\section{Excitonic Condensation in The Extended Falicov--Kimball Model} \label{sec:Model}
	
	
	\begin{figure}[tb]
		\centering
	
		\includegraphics[width=\apsfigurewidth]{fig_1}
		\caption{The non-interacting (dashed) and interacting (solid) electron bands. \revised{The interaction separates the overlapping bare bands into a semiconducting configuration. The subsequent condensation of excitons} 
	
		is glimpsed in the modest flattening of the valence band, indicating the hybridization of the original bands --- i.e. a correlated coexistence of electrons with holes.}
		\label{fig:quasiparticlebands}
	\end{figure}
	
	\heading{Model}
	We consider a simple model that exhibits excitonic condensation \revised{\cite{Batista2002, Batista2004, Ejima2014, Golez2016, Murakami2017,Tanabe2018, Murakami2020}.  The extended Falicov--Kimball model consists of itinerant spinless conduction and valence electrons, with a repulsive interaction between colocated particles:} 

	\begin{equation} \label{eq:Hubbard_Hamiltonian}
	H = \!\!\sum_{k,\,\alpha=c,v}\!\!\left(\epsilon_\alpha(k)+\Delta_\alpha \right)\tilde{c}_{k,\alpha}^\d \tilde{c}^{\phantom{\dagger}}_{k,\alpha}
	+U_0 \sum_{i=1}^{N} n_{i,c} n_{i,v}\,.
	\end{equation}
	$\tilde{c}_{k,c \left(v\right)}  = \frac{1}{\sqrt{N}} \sum_m e^{-ikm} c_{m,c \left(v\right)}$ are the annihilation operators in the conduction (valence) band, $\epsilon_{\alpha}(k)$ the band dispersions and $\Delta_\alpha$ the band offsets. $U_0 > 0$ is the on-site repulsion which drives the formation of the condensate. $N$ is the number of sites.
	We set $\hbar = 1$ throughout. The lattice constant is $a=1$, so that the first Brillouin zone (BZ) is $k\in\left[-\pi,\pi\right]$.
	We use nearest-neighbour hopping bands $\epsilon_{c,v}\left(k\right) = \mp J_{c,v} \cos {k}$, and our unit of energy is the average transfer integral $J = \frac{1}{2}\left(J_c + J_v \right)$.
	
	\revised{Ultimately we will be interested in bulk 3d systems, where long-range excitonic order is possible. This allows us to rely on mean-field theory, apply the RPA, and neglect localization due to disorder, assuming it is weak. While the numerical results below are calculated in one dimension for practicality, the formalism is readily evaluated in higher dimensions and will produce qualitatively similar results. We will point out  differences due to dimensionality where they arise. From these one-dimensional results we will infer the properties of 3d systems. This will be supported by an effective field theory which we solve in three dimensions.}  
	
	
	$H$ has a manifest $U(1)$ symmetry, as it does not mix states in different bands. In the excitonic insulator phase this symmetry will be broken, and the ground state of the system will have a non-zero value for the complex-valued order parameter 
	\begin{equation}
	\phi_i =\braket{c^\d_{i,c} c^{\phantom{\dagger}}_{i,v}} \neq 0 \,.
	\end{equation}
	
	We determine the order parameter for this condensed ground state within the mean-field approximation, using standard methods as described in Appendix \ref{app:GroundStateProperties}.
	As a starting point for our model parameters we choose $\Delta_c = -0.55~J,\, \Delta_v = -2.45~J,$ $U = 3.0~J$ and $J_c=J_v=J$. These values \revised{produce a semiconducting normal phase}, match \cite{Murakami2017} angle-resolved photoelectron spectroscopy measurements performed on \ce{Ta2NiSe5}{} \cite{Mor2017}, and place the model in the BEC regime \cite{Phan2010} (see Appendix \ref{app:GroundStateProperties}). However, we modify the transfer integrals to explicitly break particle--hole symmetry, setting $J_v = 1.2~J$ and $J_c = 0.8~J$. This is because electron--hole symmetry will protect the CMs from certain disorder scattering mechanisms, but in general we do not expect such a fine-tuned symmetry. \revised{The bare and interaction-driven quasiparticle bands, $\epsilon_\alpha (k)$ and $E_\pm(k)$, respectively,} are depicted in Fig.~\ref{fig:quasiparticlebands}. The system is at half-filling if the chemical potential is fixed at $\mu=0$, and we henceforth suppress it.
	
	
	\heading{Linear response}
	We investigate the collective modes through their signatures in the system linear response \cite{Murakami2020, Murakami2020Phonons}. 
	%
	%
	The observables whose response we are interested in are the electron density operators. We define them in real space as the local spinor vectors
	\begin{equation}
	\hat{n}^\mu_i = \Psi_i^\d \sigma^\mu \Psi_i^{\phantom{\dagger}},\quad (\mu=0,\,x,\,y,\,z)
	\end{equation}
	with $ \sigma^0 = \mathds{1}$ the identity matrix and the rest Pauli's, and defining the vector of annihilation operators $\Psi_i~=~(c_{i,c},~c_{i,v})^{\rm T}$, with the corresponding Fourier transform $\tilde{\Psi}_k$. The mean fields correspond to the expectation values of $\hat{n}^\mu$, as $\braket{\phi} = \frac{1}{2} \braket{\hat{n}^x + i \hat{n}^y }$ and $\braket{n_{c,v}} = \frac{1}{2} \braket{\hat{n}^0 \pm \hat{n}^z}$. $\hat{n}^0$ is the total electron density operator.
	
	The susceptibility is given by the Kubo formula \cite[Ch. 7]{AltlandSimons2010}
	\begin{align} 
	\chi_{\mu\nu}\left(i,t\right) 
	& = -i\Theta\left(t\right) \frac{\braket{GS|\left[ \hat{n}_i^\mu\left(t\right) , \hat{n}_0^\nu \left(0\right) \right]|GS}}{\braket{GS|GS}}\,, \label{ChiWickTheoremNew} 
	\end{align}
	evaluated with respect to the (interacting) ground state (GS) of the system, yielding the retarded response function. 
	This susceptibility is related to the time-ordered response function, which admits a diagrammatic expansion via Wick's Theorem \cite{PeskinSchroeder1995}. This is written as \cite{FetterWalecka}
	\begin{align} \label{eq:ChiWickTheoremConnected} 
	\chi_{\mu\nu}\left(i,t\right) 
	& = -i \Braket{GS\left|\mathcal{T} \lbrace \hat{n}_i^\mu\left(t\right) \hat{n}_0^\nu \left(0\right) \rbrace \right|GS}_\text{conn.}^\text{R}
	 = -i \times  \,
	\left[\sigma^\mu
	\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
	\begin{feynman}
	\vertex (m) at (-1.5,0) ;
	\vertex (n) at (+1.5,0) ;
	\vertex (ru) at (+0.5, +0.5) ;
	\vertex (rd) at (+0.5, -0.5) ;
	\vertex (lu) at (-0.5, +0.5) ;
	\vertex (ld) at (-0.5, -0.5) ;
	\diagram* {
		(m) -- [fermion, bend left = 20] (lu) -- (ru) -- [fermion, bend left = 20] (n),
		(n) -- [fermion, bend left = 20] (rd) -- (ld) -- [fermion, bend left = 20] (m),
	};
	\draw [fill=gray] (ld) rectangle (ru);
	\end{feynman}
	\end{tikzpicture} \, \sigma^\nu \right]^\text{R} \,.
	\end{align}
	Here $\mathcal{T}$ is the time-ordering operator. The notation $\left[ \dots \right]^\text{R}$ signifies that the linear response function $\chi$ is obtained by retarding the quantity in brackets, i.e. pushing all its poles to the lower half of the complex frequency plane \cite[p. 174]{FetterWalecka}. Henceforth we carry out calculations in terms of time-ordered functions, and the retardation at the end of the procedure is implied. The diagrammatic expansion for $\chi$ is the sum of all connected diagrams with two external vertices, as shown.
	
	The collective modes are most easily identified in $k$-space. 
	Assuming translational invariance, one obtains  
	\begin{equation} \label{eq:Chi0KuboMomentumFrequencySpace}
	\chi_{\mu\nu} \! \left(k,\omega\right) 
	= \left[ \frac{1}{\omega+i\eta} \right]\!\! * \!\!
	\int_{-\infty}^{\infty}\!\!\!\!\!
	\Braket{\left[  \hat{n}_{k}^\mu\!\left(t \right),  
		\hat{n}_{-k}^\nu\!\left(0\right) \right]} e^{i\omega t} dt\,,
	\end{equation}
	where $\hat{n}_{k}^\mu = \frac{1}{\sqrt{N}} \sum_p \tilde{\Psi}^\d_{p-k}\sigma^\mu\tilde{\Psi}_p$ is the Fourier transform of $\hat{n}_i^\mu$ and $*$ denotes a convolution. With $\eta$ a positive infinitesimal, the first term is the Fourier transform of $-i\Theta\left(t\right)$, which introduces the poles of $\chi$. For a time-translation invariant system the second term is self-adjoint under the transposition $\mu \leftrightarrow \nu$, implying that the susceptibility has the analytical structure $\chi\left(\omega\right)=\sum_\text{poles}\frac{\text{Hermitian matrix}}{\text{simple pole}}$. We thus define the Hermitian spectral weight function $\mathcal{A}_{\mu\nu}(k,\omega)$ by
	\begin{equation}\label{eq:chiSpectralPoleDecomposition}
	\chi_{\mu\nu}(k,\omega)=
	\int_{-\infty}^{\infty}\frac{\mathcal{A}_{\mu\nu}(k,\omega')}{\omega^+-\omega'}d\omega'
	\end{equation}
	where $\omega^+ = \omega + i\eta$. 
	%
	The modes of the system are located where $\mathcal{A}_{\mu\nu}\left(k,\omega\right)$ is non-zero. In the excitonic insulating state we expect a continuous region of modes above the gap energy, corresponding to excitations of electrons from the valence to the conduction band. By contrast, the collective modes that lie below this particle-hole continuum will form a sharp line in the $\left(k,\omega\right)$ plane, with a dispersion relation $\omega=\omega_k$. As the collective modes are gapless, they should be separated from the continuum, at least for small momenta. 
	%
	While Eq.~\eqref{eq:chiSpectralPoleDecomposition} is suggestive of the Kramers-Kronig relation satisfied by the retarded $\chi$ \cite{AltlandSimons2010}, we remark that it does not define $\mathcal{A}_{\mu\nu}$ uniquely. Instead, leveraging the identity $\mathrm{Im}\left(\omega^+ - \omega'\right)^{-1} = -\pi \delta \left(\omega - \omega'\right)$ along with $\mathcal{A}$'s hermiticity, we invert this definition and find
	\begin{equation} \label{eq:SpectralWeightFunctionASolution}
	\mathcal{A}_{\mu\nu}
	= 	\frac{1}{2\pi i} \left[\left( \chi_{\mu\nu} \right)^\d - \chi_{\mu\nu}  \right]\,.
	\end{equation}
	
	\heading{Collective modes}
	The expansion \eqref{eq:ChiWickTheoremConnected} begins by evaluating the Kubo formula with respect to some initial ground state ansatz $\ket{0}$ which approximates the true interacting ground state. For this we take the ground state found within the mean-field approximation. This solution is worked out in Ref.~\cite{Murakami2017} and outlined in Appendix~\ref{app:GroundStateProperties}, where we also list the resulting Feynman diagrammatic rules. 
	We employ the random-phase approximation (RPA) \cite{Bohm1953, Pines1953, Chen2017}. Within the RPA, the interacting susceptibility is given by a Dyson equation \cite{FetterWalecka}, represented schematically by
	\begin{equation}
	\underset{\chi^\text{RPA}}{\includegraphics[raise=-3.5mm]{bubbleRPA.pdf}} = \underset{\chi^0}{\includegraphics[raise=-3.5mm]{bubble.pdf}} 
	+ \underset{\,\,\chi^0 \, * \,  U^\text{eff} \, * \, \chi^\text{RPA}}{\includegraphics[raise=-3.5mm]{bubble.pdf}\includegraphics[raise=-3.5mm]{bubbleRPA.pdf}}\,.
	\end{equation}
	%
	Here the clear bubble $\chi^0$ is the bare susceptibility, and $U^\text{eff}$ is a numerical vertex factor corresponding to the projection of the electron--electron interaction operator onto the sector or bubble diagrams.
	%
	The convolutions are turned into simple (matrix) products in momentum space, which we will use exclusively henceforth.
	%
	Isolating $\chi^\text{RPA}$, we obtain the familiar RPA result \cite[cf. Sec. 3]{Chen2017}
	\begin{equation} \label{eq:chiRPASolution}
	\chi^\text{RPA}_{\mu\nu}\left(k,\omega\right)
	=\chi^0_{\mu\nu}
	+\sum_{\lambda,\eta} \chi^0_{\mu\lambda} U^\text{eff}_{\lambda\eta}  \chi^\text{RPA}_{\eta\nu} 
	=\left[\frac{\chi^0}{1-U^\text{eff}\chi^0}\right]_{\mu\nu}\,.
	\end{equation}
	%
	Here the solution for $\chi^\text{RPA}$ is shorthand notation for matrix inversion.
	%
	Note that Eq.~\eqref{eq:chiRPASolution} is a $4\times4$ matrix equation;  $U^\mathrm{eff} = \left(U_0/2\right) \mathrm{diag}\left(1,-1,-1,-1\right)$ and $\chi^0$ are computed in Appendix \ref{app:GroundStateProperties}. For $k=0$ the components $\mu,\nu=x,\,y,\,z$ do not mix with the $\mu=\nu=0$ sector, as $\hat{n}^0_{k=0}$  is the total particle number and is conserved \cite[SI]{Murakami2017}. However, this does not hold for general $k$, and we must consider this equation in full \cite{particle-hole}.
	
	
	
	\begin{figure}[tb]
		\centering
	
		\includegraphics[width=\apsfigurewidth]{fig_2}
		\caption{Component $\mathcal{A}_{yy}$ of the spectral function (\ref{eq:SpectralWeightFunctionASolution}),  corresponding to order parameter oscillations in the ``phase direction''. The colour scale is logarithmic (in units $J^{-1}$), showing the singularities at the poles; the broadening regulator used here is \revised{ $\eta = 0.005~J$.} The system size is $N = 1024$ sites. The gapped two-particle excitation continuum is visible at high energies. In the symmetry-broken EI, a Goldstone branch of CMs appears, and considerable spectral weight is shifted to it from the continuum.}
		\label{fig:Chi0andChiRPA}
	\end{figure}
	
	
	
	\heading{Results}
	The RPA susceptibility is evaluated and plotted in Fig.~\ref{fig:Chi0andChiRPA}. The displayed quantity is $\mathcal{A}_{yy}$, which corresponds to the spectral density of modes that oscillate in the ``phase direction'' of the order parameter. A non-zero \revised{$\eta = 0.005~J$} is used to broaden mode poles. 
	%
	The bottom of the two-particle continuum is visible, occurring at the frequency of the inter-band gap seen in Fig.~\ref{fig:quasiparticlebands}. 
	%
	Additionally, a separated branch of modes $\omega_k$ appears at lower energies, revealing the expected Goldstone-like branch of CMs.

	At fixed $k$, this branch absorbs a considerable portion of the spectral weight which was once in the two-particle continuum in the bare $\chi^0$. 
	%
	The eigenvectors of $\mathcal{A}_{\mu\nu}(k,\omega_k)$ indicate the nature of the collective modes, and at small $k$ these are predominantly phase modes.
	
	
	\section{Disorder Scattering} \label{sec:Disorder}
	
	The focus of this paper is to determine the disorder-induced scattering rate of excitations. We add to the model an external disorder potential and a corresponding one-particle vertex,
	\begin{align} \label{eq:DisorderHamiltonian}
	H_\text{dis} & = \feynmandiagram [small, baseline={([yshift=-.5ex]current bounding box.center)}, layered layout, horizontal=a to b] {c -- [fermion] a -- [insertion={[style=red,size=5pt]0.5}]  b -- [fermion] d}; 
	
	= \sum_{\lambda} \sum_{i} V^\lambda_i \Psi^\d_i\sigma^\lambda\Psi_i^{\phantom{\dagger}} 
	= \frac{1}{\sqrt{N}} \sum_{\lambda} \sum_{k,q} \tilde{V}^\lambda_q \tilde{\Psi}^\d_{k+q} \sigma^\lambda \tilde{\Psi}_k^{\phantom{\dagger}} \,,
	\end{align}
	where $V^\lambda_i$ ($\lambda = 0,\,x,\,y,\,z$) is a real random disorder field with Fourier transform $\tilde{V}^\lambda_q$. We assume that $V^\lambda$ is normally distributed with zero mean and the following covariance relations
	represented by a ``disorder propagator''
	\begin{equation} \label{eq:DisorderPotentialCovarianceMatrix}
	\feynmandiagram [small, baseline={([yshift=-0.5ex]current bounding box.center)}, horizontal=a to b] {
		a -- [red, boson] b
	};=
	\left< V^\lambda_i V^{\lambda'}_j \right>  = K_{\left|i-j\right|}^{\lambda \lambda'} \Leftrightarrow
	\left< \tilde{V}^\lambda_q \tilde{V}^{\lambda'}_{-q} \right> =  \tilde{K}^{\lambda \lambda'}_q\,.
	\end{equation}
	
	For simplicity we henceforth assume the disorder potential is delta-correlated in real space, $\braket{V^\lambda_i V^{\lambda'}_j}~\propto~\delta_{ij}$, signifying an ensemble of point scatterers. Furthermore we take the different disorder channels to be uncorrelated.  We thus set $\tilde{K}^{\lambda\lambda'}_q = K_0 \delta_{\lambda \lambda'}$. 
	
	The four disorder channels can be classified with respect to the $U(1)$ symmetry of the Hamiltonian: Channels $\lambda = 0,\,z$ preserve the symmetry, while channels $\lambda = x,\, y$ violate it. We will show that the symmetry of the disorder has crucial consequences for the scattering of the CMs. Physically, channel $\lambda = 0$ corresponds to an external electrostatic potential.
	
	Disorder will induce a finite lifetime for all excitations, which at the level of RPA were infinitely long lived.
	This will manifest in the broadening of the poles of the single- and two-particle response functions. 
	We will obtain these induced scattering rates by computing the electron and CM (proper) self-energy \cite{FetterWalecka}. 
	Working to second order in $H_\mathrm{dis}$, this is equivalent to computing lifetimes with Fermi's Golden Rule (FGR) \cite{Sakurai2017}. However, the diagrammatic method is useful for computing the scattering rate of the CMs, whose wavefunctions are not found explicitly. 
		
	\revised{Our work extends that of Zittartz \cite{Zittartz1967}, who calculated the lifetime of the EI quasiparticle excitations in the presence of disorder. The present work builds upon it in three aspects: (i) We also compute the lifetime of the collective modes. (ii) We consider all possible disorder channels;  Zittartz's ``normal impurities'' correspond to our $\lambda=0$ disorder channel. (iii) We take into account the fact that the disorder potential is dressed by the condensate (see below). This turns out to be crucial in obtaining the correct lifetimes of the collective modes when scattered by symmetry-respecting disorder, such as Zittartz's.}
	
	
	\heading{Disorder dressing}
	The inclusion of disorder modifies the ground state of the system above which the collective modes are excited. This manifests as a change in the value of the mean fields that the electronic degrees of freedom are coupled to. This modulation corresponds to the condensate shifting to screen the disorder. Assuming the disorder is perturbative, this can be accounted for to leading order by using linear response, as $\delta \braket{n^\mu}\left(q\right) = \sum_\nu \chi^\mathrm{RPA}_{\mu\nu}\! \left(q,\omega = 0\right)\tilde{V}^\nu_q$. This implies a dressed disorder vertex
	\begin{equation}
	\tilde{\mathbb{V}}_q^\lambda = \tilde{V}_q^\lambda + \sum_{\mu,\nu} U^\mathrm{eff}_{\lambda \mu} \chi^\mathrm{RPA}_{\mu\nu} \!\left(q, \omega = 0\right)\tilde{V}_q^\nu\,,
	\end{equation}
	which we absorb into a dressed disorder propagator
	\begin{align}
	& \braket{\tilde{\mathbb{V}}_q^\lambda \tilde{\mathbb{V}}_{-q}^{\lambda'}} 
	= \tilde{\mathbb{K}}_q^{\lambda \lambda'}
	= \sum_{\mu,\nu} \left[\mathds{1} + U^\mathrm{eff}\chi^\mathrm{RPA}_{\omega=0}\right]_{\lambda \mu}\! \tilde{K}_q^{\mu\nu} \!\left[\mathds{1} + U^\mathrm{eff} \chi^\mathrm{RPA}_{\omega=0}\right]_{\nu\lambda'}^\d\,.
	\end{align}
	%
	Diagrammatically, this is equivalent to the screening
	\begin{align}
	& \feynmandiagram [small, baseline={([yshift=-0.5ex]current bounding box.center)}, horizontal=a to b] {
		a -- [red, thick, boson] b
	};
	= \disorderprop + \disorderprop \includegraphics[raise=-3.5mm]{bubbleRPA.pdf} \, +   \includegraphics[raise=-3.5mm]{bubbleRPA.pdf} \disorderprop + \includegraphics[raise=-3.5mm]{bubbleRPA.pdf} \disorderprop \includegraphics[raise=-3.5mm]{bubbleRPA.pdf} \,.
	\end{align}
	%
	Though the bare disorder channels could be uncorrelated (i.e. $\tilde{K}_q$ is diagonal), in general the channels become mixed in $\tilde{\mathbb{K}}_q$. However, it can be shown [Appendix \ref{app:GroundStateProperties}] that channel $\lambda = y$ remains uncorrelated with the rest.
	%
	The dressing of the disorder is a crucial effect if the disorder respects the {$U(1)$ symmetry of  the clean system} (shown below). 
	
	
	
	
	\subsection{Electron propagator}
	
	
	\begin{figure}[tb!]
		\centering
	
		\includegraphics[width=\apsfigurewidth]{fig_3}
		\caption{Single-particle scattering rates, for different scattering channels. 
			The mode energy is measured relative to the bottom of the conduction band \revised{$E_{k=0}$}.
			The peaks at the band edges are van Hove singularities in the one-dimensional density of states. 
			The features in channels $\lambda=x,z$ are explained in the text. The vertical dashed line marks the bottom of the electron--hole two-particle continuum, for reference.
			(Here $N = 1024$.)}
		\label{fig:SingleParticleDecayRates}
	\end{figure}
	
	
	\heading{Quasiparticle scattering rates}
	The dressed propagator can be determined by the usual Dyson equation   
	\begin{equation} \label{eq:SingleParticleDysonEquation}
	G = G^0 + G^0 \Sigma G 
	\,\Leftrightarrow\,
	G^{-1} = \left(G^0\right)^{-1} - \Sigma\,,
	\end{equation}
	with $\Sigma$ the proper self-energy insertion \cite{FetterWalecka}. Working to second order in $V^\lambda$, it is found by evaluating the diagram 
	\begin{equation} \label{eq:SingleParticleSelfEnergyDiagram}
	\Sigma\left(k,\omega\right) = -i \times
	\feynmandiagram [inline=(b.base),layered layout, horizontal=a to b] {
		a -- [fermion] b -- [boson, red, thick, half right] a
	};
	= 
	\frac{1}{N} \!\! \sum_{\lambda, \lambda',q} \!\! \tilde{\mathbb{K}}_q^{\lambda\lambda'} \sigma^{\lambda} G^0 \!\left(k + q,\omega\right) \sigma^{\lambda'}.
	\end{equation}
	As the disorder propagator is statistical rather than quantum-mechanical, tadpole diagrams are considered disconnected and cancel by the linked cluster theorem.
	
	Furthermore, $G^0\left(k,\omega\right)$ is diagonalized by a unitary matrix $R_k$ which rotates it to the basis of the hybridized quasiparticle excitations $\phi^\d_\pm\left(k\right)$ [Appendix~\ref{app:GroundStateProperties}], where it takes the form
	%
	\begin{equation} \label{eq:SingleParticleBareGDiagonalized}
	G^0\left(k,\omega\right) = R_k \begin{pmatrix}
	\omega^+ - E_+(k) & \\
	& \omega^-  - E_-(k)
	\end{pmatrix}^{-1} R_k^\d \,.
	\end{equation}
	%
	Thus, Eq.~\eqref{eq:SingleParticleDysonEquation} reads
	\begin{equation}
	G^{-1} = R_k \begin{pmatrix}
	\omega^+ - E_+  - \tilde{\Sigma}_{11} &  -\tilde{\Sigma}_{12} \\
	-\tilde{\Sigma}_{21} & \omega^-  - E_- - \tilde{\Sigma}_{22}
	\end{pmatrix} R_k^\d \,
	\end{equation}
	where $\tilde{\Sigma}\left(k, \omega\right) = R_k^\d \Sigma\left(k,\omega\right) R_k$ is  the self-energy ``rotated'' to the same eigenbasis. The poles of the dressed propagator are where it is singular, i.e. $\omega=E_k$ at which $\det{G^{-1}\left(k, E_k\right)}=0$. The off-diagonal elements of $\tilde{\Sigma}$ only come into play in post-leading contributions to the new $E_k$. Focusing henceforth on the positive energy poles, it is now shifted and broadened by 
	\begin{equation} \label{eq:SingleParticleSelfEnergyRate}
	\delta E_k = \mathrm{Re}\,{\tilde{\Sigma}_{11}\!\left(k,E_k\right)}, \quad 
	\Gamma^\text{QP}_k = - \mathrm{Im}\,{\tilde{\Sigma}_{11}\!\left(k,E_k\right)}.
	\end{equation} 
	
	The matrices $R_k$ are real, 
	%
	\begin{equation}
	R_k = \exp\left({\frac{i}{2}\sigma^y\theta_k}\right) = 
	\begin{pmatrix}
	\phantom{-}	\cos{\frac{\theta_k}{2}}	\,& \sin{\frac{\theta_k}{2}} \\ 
	-		\sin{\frac{\theta_k}{2}}	\,& \cos{\frac{\theta_k}{2}}
	\end{pmatrix}
	\end{equation}
	%
	where the mixing angle $\theta_k$ which hybridizes the original bands is the polar angle of the pseudomagnetic field $\vec{B}_k=\left(-\sin{\theta_k},0,\cos{\theta_k}\right)|B_k|$ obtained in the mean-field solution [Appendix~\ref{app:GroundStateProperties}]. They thus do not mix the real and imaginary parts of the self-energy.
	
	
	
	The self-energy diagram is computed in Appendix \ref{app:QuasiparticleSelfEnergy}. We find the scattering rate to be 
	\begin{equation} \label{eq:SingleParticleFermiGoldenRuleRateBothRates}
	\Gamma^{\mathrm{QP}}_k \!\!
	=\! \frac{\pi}{2}  g\!\left(E_k\right)  \revised{  \langle \Big| \! \sum_\lambda \! \tilde{\mathbb{V}}_{-2k}^{\lambda} \hat{B}_{k}^{\lambda} \Big|^2 \rangle  \!}
	= \! \frac{\pi}{2}  g\!\left(E_k\right) \! \sum_{\lambda,\lambda'} \! \hat{B}^\lambda_k  \tilde{\mathbb{K}}_{-2k}^{\lambda \lambda'} \hat{B}^{\revised{\lambda'}}_k,
	\end{equation}
	where $\hat{B}_k^0 = 1$. Here we have put the scattering rate into a form recognizable as FGR \cite{Sakurai2017}, and $g\left(E\right)$ is the single-particle density of states (DoS) per site. 
	The last inner product is the matrix element for elastic backward scattering with momentum transfer $q = -2k$. {Here we discarded the $q=0$ matrix element, as forward scattering does not affect transport.}
	Note that $\hat{B}_k^y=0$. As the $\tilde{\mathbb{K}}_q^{yy}$ element is not mixed with other elements in the screening process, we find that the $\lambda=y$ channel has a vanishing scattering rate.
	
	\revised{This result qualitatively reduces to that of Zittartz (cf. Eq.~(31) in Ref.~\cite{Zittartz1967}) if we only consider the electrostatic channel ($\lambda=\lambda'=0$), and neglect the disorder dressing (i.e. use $K^{00}$ instead of $\mathbb{K}^{00}$). However, the dressed $\lambda=0$ channel generally mixes with additional channels, giving rise to a momentum dependence other than through the density of states.}
	
	\heading{Results}
	Eq.~\eqref{eq:SingleParticleFermiGoldenRuleRateBothRates} is evaluated and the contribution of each channel $\lambda$ is plotted in Fig.~\ref{fig:SingleParticleDecayRates}. \revised{The quasiparticle energies are measured relative to the bottom of the (interacting) conduction band $E_{k=0}$.} The sharp peaks at the band edges are the one-dimensional van Hove singularities in the DoS.  
	It is apparent that the different scattering channels induce scattering rates of the same order of magnitude at most energies. The rate \eqref{eq:SingleParticleFermiGoldenRuleRateBothRates} can be interpreted as the inner product between the disorder potential and the pseudomagnetic field. The screening ``rotates'' the disorder channels, such that the potential may be orthogonal to the field at specific $k$ points, and this is the origin of the dips in channels $\lambda = x, z$. It is curious that this allows the backscattering rate to vanish. However, this is an artifact of one spatial dimension, as in higher dimensions angular summation over different momentum transfers will produce a non-zero rate.
	\subsection{Density-density response function}
	
	\heading{Diagrammatic expansion}
	The CM scattering rates may be found by resumming diagrammatic insertions to obtain the leading-order contributions to the exciton self-energy. For this to yield the scattering rates, we need to find an effective self-energy which satisfies a Dyson-type equation, that is schematically
	\begin{equation} \label{eq:ChiDysonExpansionSimple}
	\chi^\text{dis} = \chi^\text{RPA} + \chi^\text{RPA} * \Sigma_\text{eff} * \chi^\text{dis} 
	= \left[\left(\chi^\text{RPA}\right)^{-1} - \Sigma_\text{eff}\right]^{-1} \,.
	\end{equation}
	We enumerate all the RPA-type second-order disorder diagrams which enter expansion \eqref{eq:ChiWickTheoremConnected}, which are 
	\begin{widetext}
		\begin{equation} \label{eq:DisorderDiagrams}
		\Sigma \,=\, \includegraphics[raise=-3.5mm]{bubbleDisordered.pdf} \,=\, \bubbleClam \,+\, \bubblePurse \,+\, 
		\begin{tikzpicture}[baseline={([yshift=-.5ex]m.base)}]
		\begin{feynman}
		\vertex (m) at (0,0) ;
		
		\vertex (a) at (-4.0, 0);
		\vertex (b) at (+4.0,0);
		
		\vertex (u1) at (-0.5*4.0,4.0*0.414213*0.5);
		\vertex (u2) at (+0.5*4.0,4.0*0.414213*0.5);
		
		\diagram* {
			
				(a) -- [out = 45, in = 180] (u1) -- [out = 0, in = 135] (m) -- [out = -135, in = -45] (a),
				(m) -- [out = 45, in = 180] (u2) -- [out = 0, in = 135] (b) -- [out = -135, in = -45] (m)},
			(u1) -- [boson, red, thick, quarter left] (u2),
		};
		\end{feynman}
		\end{tikzpicture}
		\,+\,
		\begin{tikzpicture}[baseline={([yshift=-.5ex]m.base)}]
		\begin{feynman}
		\vertex (m) at (0,0) ;
		\vertex (n) at (+4.0,0) ;
		
		\vertex (a) at (-4.0, 0);
		\vertex (b) at (2*4.0,0);
		
		\vertex (u1) at (-0.5*4.0,4.0*0.414213*0.5);
		\vertex (u2) at (+1.5*4.0,4.0*0.414213*0.5);
		
		\diagram* {
			(m) -- [fermion, quarter left]  (n) -- [fermion, quarter left] (m),
			
				(a) -- [out = 45, in = 180] (u1) -- [out = 0, in = 135] (m) -- [out = -135, in = -45] (a),
				(n) -- [out = 45, in = 180] (u2) -- [out = 0, in = 135] (b) -- [out = -135, in = -45] (n)},
			(u1) -- [boson, red, thick, quarter left] (u2),
		};
		\end{feynman}
		\draw [pattern=north east lines] (m) arc [start angle=135, end angle=45, radius=4.0/sqrt(2)] arc [start angle=-45, end angle=-135, radius=4.0/sqrt(2)];
		
		\end{tikzpicture}\,.
		\end{equation}
	
		Here the direction of the fermion propagators is left unspecified to indicate that the disorder potential could couple to either the electron or hole propagator in each bubble; therefore, the diagrams have 1, 2, 4, and 4 realizations, respectively. The different diagrams in Eq.~\eqref{eq:DisorderDiagrams} are computed and their properties discussed in Appendix \ref{app:ExcitonSelfEnergyDiagrammatics}. 
		
		We can now systematically generate the diagrammatic expansion for the disordered RPA susceptibility $\chi^\text{dis}$ by concatenating arbitrary numbers of clean and disordered bubble diagrams. This is represented schematically by
	
		\begin{align} \label{eq:DisoerderResummation}
		i \chi^\text{dis} & \,=\, \includegraphics[raise=-3.5mm]{bubbleRPADisordered.pdf} \,=\, 
		\includegraphics[raise=-3.5mm]{bubbleRPA.pdf} \,+\, \includegraphics[raise=-3.5mm]{bubbleDisordered.pdf} \nonumber
		\,+\, \includegraphics[raise=-3.5mm]{bubbleRPA.pdf}\includegraphics[raise=-3.5mm]{bubbleDisordered.pdf} \,+\, \includegraphics[raise=-3.5mm]{bubbleDisordered.pdf}\includegraphics[raise=-3.5mm]{bubbleRPA.pdf} \,+\,
		\dots \nonumber \\ &
		= \includegraphics[raise=-3.5mm]{bubbleRPA.pdf} + \Big(1+\includegraphics[raise=-3.5mm]{bubbleRPA.pdf}\Big)\includegraphics[raise=-3.5mm]{bubbleDisordered.pdf}
		\sum_{n=0}^\infty \Bigg[\Big(1+\includegraphics[raise=-3.5mm]{bubbleRPA.pdf}\Big)\includegraphics[raise=-3.5mm]{bubbleDisordered.pdf}\Bigg]^n
		\Big(1+\includegraphics[raise=-3.5mm]{bubbleRPA.pdf}\Big) \nonumber \\ &
		= \Bigg[\includegraphics[raise=-3.5mm]{bubble.pdf} \,+\, \includegraphics[raise=-3.5mm]{bubbleDisordered.pdf} \Bigg] \Bigg[1 \,+\, \includegraphics[raise=-3.5mm]{bubbleRPADisordered.pdf}\Bigg]\,.
		\end{align}
	\end{widetext}
	Note that two consecutive clean RPA bubbles are not allowed, as this leads to double counting. Therefore, the Dyson relation in the last line displays the bare $\chi^0$ instead of $\chi^\text{RPA}$, and it cannot be used to extract the effective self-energy in Eq.~\eqref{eq:ChiDysonExpansionSimple}.
	
	Instead, we focus on the second line of Eq.~\eqref{eq:DisoerderResummation}, which reads
	\begin{align} \label{eq:DisorderSelfEnergyDefinition}
	\chi^\text{dis} & = \chi^\text{RPA} + \left(1 + \chi^\text{RPA} * U^\text{eff}\right) * \Sigma ~ *  
	\sum_{n=0}^\infty \left[ \left(  1 + U^\text{eff} \!*\!  \chi^\text{RPA} \right) 
	\! * \! U^\text{eff} \!*\! \Sigma  \right]^n \!\! * 
	\left(1 + U^\text{eff} \!*\! \chi^\text{RPA} \right).
	\end{align}
	Adding $\left(U^\text{eff}\right)^{-1}$ to both sides and regrouping products, we can rewrite

	\begin{equation} \label{eq:EffectivePropagator}
	\mathcal{G}^\text{dis} = 
	\mathcal{G}^\text{RPA} + 
	\mathcal{G}^\text{RPA} * \Sigma^\star * \mathcal{G}^\text{RPA}
	\end{equation}
	%
	with $\mathcal{G}^\text{dis, RPA} = \left(U^\text{eff}\right)^{-1} + \chi^\text{dis, RPA}$ and the self-energy insertion
	\begin{equation} \label{eq:ExcitonImproperSelfEnergy}
	\Sigma^\star = \left[U^\text{eff} * \Sigma * U^\text{eff} \right] * 
	\sum_{n=0}^{\infty} \left[ \mathcal{G}^\text{RPA}*
	\left[  U^\text{eff} * \Sigma * U^\text{eff}  \right]\right]^n\,.
	\end{equation}
	
	This solution illustrates the analogy between $\chi$ and the exciton propagator: 
	%
	From Eqs.~\eqref{eq:EffectivePropagator} and \eqref{eq:ExcitonImproperSelfEnergy} we identify $\mathcal{G}^\text{RPA,dis}$ as the ``bare'' and ``dressed'' effective exciton propagators, respectively, and
	$\Sigma^\star$ as the \emph{improper} self-energy insertion \cite[Sec. 9]{FetterWalecka}. 
	%
	The effective \emph{proper} self-energy is then read from the definition 
	$\Sigma^\star = \Sigma_\text{eff} * \sum_n \left[ \mathcal{G}^\text{RPA} * \Sigma_\text{eff} \right]^n$. 
	We thus conclude
	\begin{align}
	\mathcal{G}^\text{dis} & 
	= \mathcal{G}^\text{RPA} \!+ \mathcal{G}^\text{RPA} \!* \Sigma_\text{eff} * \mathcal{G}^\text{dis} 
	=\left[\left(\mathcal{G}^\text{RPA}\right)^{-1} \!- \Sigma_\text{eff}\right]^{-1}\!\!,
	\\
	\Sigma_\text{eff} & = U_\text{eff} * \Sigma * U_\text{eff}\,. \label{eq:ExcitonEffectiveSelfEnergy}
	\end{align}
	
	Note that $U^\text{eff}$ is regular, and therefore $\mathcal{G}$ has no other poles except those of $\chi$. Additionally, $U^\text{eff}$ is itself Hermitian, and therefore we can substitute $\mathcal{G}$ into Eq.~\eqref{eq:SpectralWeightFunctionASolution} instead of $\chi$ without affecting $\mathcal{A}$.
	
	\heading{Scattering rates}
	Given the effective excitonic self-energy, we now need to extract from it the collective mode scattering rates in a procedure analogous to the ``rotation'' used to extract the single-particle rates. 
	
	We recall that the collective modes are located where $\mathcal{A_{\mu\nu}}(k,\omega)$ is non-zero. 
	Specifically, each collective mode at $\left( k, \omega_k \right)$ is associated with an eigenvector of $\mathcal{A_{\mu\nu}}\left( k, \omega_k \right)$. Although it is a $4 \times 4$ matrix, only eigenvectors with non-zero eigenvalues represent physical modes. 
	Thus, we decompose the weight function as
	\begin{align} \label{eq:SpectralADecomposition}
	\mathcal{A} \!\left( k, \omega \right) & = 
	\delta \! \left( \omega - \omega_k \right) \mathcal{R}_k 
	\begin{pmatrix}
	\hspace{0.00cm} \mathcal{A}_1(k) \hspace{1.50cm} \\ 
	\hspace{0.50cm} \mathcal{A}_2(k) \hspace{1.00cm} \\ 
	\hspace{1.00cm} \mathcal{A}_3(k) \hspace{0.50cm} \\ 
	\hspace{1.50cm} \mathcal{A}_4(k) \hspace{0.00cm}
	\end{pmatrix} \mathcal{R}_k^\d  + \text{sum over continuum modes}\,,
	\end{align}
	where the diagonalizing matrices $\mathcal{R}_k$ are defined such that the real eigenvalues $\mathcal{A}_i\!\left(k\right)$ are ordered by decreasing magnitude. Furthermore, we only expect one soft mode, and thus that only $\mathcal{A}_1$ is non-zero. {In Appendix \ref{app:ExcitonSelfEnergyDiagrammatics} we show that $\mathcal{A}_1 \sim \omega_k^{-1}$ for small $k$.}
	%
	%
	Since the continuum modes are separated from the collective modes, we will ignore their contribution. Equally we ignore the sharp CM poles at negative energies and focus on $\omega > 0$.
	
	Plugging this into Eq.~\eqref{eq:chiSpectralPoleDecomposition}, close to a collective mode pole $\left(k, \omega_k\right)$ of the clean system we have	a pole of the form
	\begin{equation}
	\left[\mathcal{G}^\text{dis}( k, \omega )\right]^{-1}  \sim 
	\mathcal{A}_1^{-1} \!\left(k\right) \left[ \omega^+ - \omega_k - \mathcal{A}_1\!\left(k\right) \tilde{\Sigma}_\text{eff} \left( k, \omega_k \right) \right],
	\end{equation}
	again with a ``rotated'' self-energy $\tilde{\Sigma}_\text{eff} \left(k, \omega \right)= \big[\mathcal{R}^\d_k \Sigma_\text{eff} \left( k, \omega \right) \mathcal{R}^{\phantom{\d}}_k\big]_{11}$. The scattering rate is then identified as 
	\begin{equation} \label{eq:ExcitonDecayRate}
	\Gamma^\text{CM}_k = -\mathcal{A}_1 \! \left( k \right)   \mathrm{Im}\left[ \tilde{\Sigma}_\text{eff} \! \left( k, \omega_k \right)\right]\,.
	\end{equation}
	\revised{Here we evaluate the self-energy at the clean pole $\omega_k$, as the shift in the pole frequency due to $\mathrm{Re}[\tilde{\Sigma}_\mathrm{eff}]$ contributes a sub-leading correction. In this we assumed that the disorder is weak such that this shift is much smaller than $\omega_{k}$. Therefore, this treatment applies to modes at frequencies $\omega_{k}$ that are larger than the gap opened by introducing a symmetry-violating disorder potential.}
	
	\begin{figure}[tb]
		\centering
	
		\includegraphics[width=\apsfigurewidth]{fig_4}		
		\caption{Disorder-induced scattering rates of the collective modes, for different scattering channels. The scattering rate for channel $\sigma^y$ is zero at the resolution of our calculation. The dashed line marks the bottom of the two-particle continuum. Above the two-particle continuum gap, the calculation is impacted (and for channel $\sigma^0$, is dominated) by additional diagrams, which represent an exciton scattering into the two-particle continuum, i.e., disassociating; this contribution is omitted. The dips in channel $\lambda = x$ are analogous to the vanishing backwards scattering rate of the quasiparticles, cf. Fig.~\ref{fig:SingleParticleDecayRates} and surrounding discussion.}
		\label{fig:CollectiveModeDecayRates}
	\end{figure} 
	
	\heading{Low-energy regime}
	We will mostly be interested in the low-energy, long-wavelength dynamics of the CMs. If we restrict our treatment to energies below the particle--hole continuum, such that a CM cannot disassociate into an unbound electron--hole pair, the CM scattering rates can be considerably simplified. In Appendix \ref{app:ExcitonSelfEnergyDiagrammatics}, we show that in this regime only the last diagram in Eq.~\eqref{eq:DisorderDiagrams} affects the mode lifetime. Furthermore, this diagram can be evaluated explicitly at energies below the gap, and in analogy to Eq.~\eqref{eq:SingleParticleFermiGoldenRuleRateBothRates} we obtain a FGR-like expression for the CM scattering rate in terms of an effective CM vertex function [Eq.~\eqref{eq:CollectiveModeFermiGoldenRule1Dimension}]. This allows us to isolate the backwards scattering rate and discard the forward scattering rate, similar to our treatment of the quasiparticles. We refer to Appendix \ref{app:ExcitonSelfEnergyDiagrammatics} for the full calculation.
	
	\heading{Results} 
	We evaluate the scattering rates of the CMs numerically and plot them in Fig.~\ref{fig:CollectiveModeDecayRates}. The scattering rate in channel $\lambda = y$ is zero at the resolution of our calculation, and therefore is not shown --- this is elucidated in the next section. The dashed vertical line marks the bottom of the two-particle continuum, above which additional scattering processes become available (not shown here). Note the qualitative difference between channel $\lambda = x$, in which the scattering rate diverges at $\omega = 0$, despite a finite density of states, and channels $\lambda = 0 , z$, in which the rate vanishes at $\omega = 0$
	
	
	\section{Discussion} \label{sec:Discussion}
	
	
	\begin{figure}[bt]
		\centering
	
		\includegraphics[width=\apsfigurewidth]{fig_5} 
		\caption{Comparison between the mean free paths of the single-particle (dashed) and collective (solid) excitations, for small wavenumbers. The dotted grey lines mark the expected power-law scaling of the rate in each channel (see body of the text).		}
		\label{fig:FreePathComparison}
	\end{figure}
	
	
	
	
	\heading{Mean free paths at $k\rightarrow0$}
	To compare the transport properties of the collective modes and the quasiparticles, we consider each excitation's mean free path, defined by $l_\mathrm{mfp} = v_g/\Gamma$, with $v_g$ the group velocity and $\Gamma$ the disorder-induced scattering rate.	We focus on the dynamics at long wavelengths where clear qualitative difference arise. 
	
	Fig.~\ref{fig:FreePathComparison} shows the mean free paths of the quasiparticle excitations (dashed curves) and collective modes (solid curves) for small $k$. 
	For the quasiparticles, one sees that there are similar scattering rates for all disorder channels. This may be understood as follows: The single-particle free paths tend to zero when $k\rightarrow 0$, because of both a vanishing group velocity and a diverging DoS. If the scattering matrix element at $k = 0$ is non-zero, then at small momenta it can be treated as constant, and as in one dimension the DoS is inversely proportional to the group velocity, we find $l_\mathrm{mfp} \propto v_g^2 \sim k^2$.
	
	In contrast, with the collective modes we can identify three distinct behaviours at small $k$: (i) $l_\mathrm{mfp} \sim k^2$ in channel $\lambda = x$, similar to the quasiparticles; (ii) $l_\mathrm{mfp} \sim k^{-2}$ in channel $\lambda = z$; and (iii) $l_\mathrm{mfp} \sim k^{-6}$ in channel $\lambda = 0$. We stress that the latter two are obtained only if we include the effect of disorder screening, as discussed in Sec.~\ref{sec:Disorder}, showing that this is a crucial effect. As we now explain, this trimodality is related to the symmetry of each disorder channel: channel $x$ does not respect any of the symmetries of model \eqref{eq:Hubbard_Hamiltonian}, channel $z$ preserves its $U(1)$ condensate symmetry, and channel $\lambda=0$ also has manifest particle--hole symmetry. 
	
	\revised{We remark that these scaling laws are cut off at very small $k$, when the mode energy $\omega_k$ becomes comparable to a possible CM gap, whether intrinsic (in the case of an approximate symmetry) or induced by disorder.}
	
	\heading{Effective field theory}
	The properties that we have computed at long wavelengths can be interpreted within an effective field theory. This also allows the results to be extended to arbitrary dimension $D$.	For this purpose we consider a classical field theory for the order parameter of the condensate. To accommodate the three disorder channels discussed so far we must consider a vector order parameter; the natural field to examine is the pseudo-spin $\vec{n}_i = \frac{1}{2}\Psi^\d_i \vec{\sigma} \Psi^{\phantom{\d}}_i$ at each site.
	Taking a continuum limit, we write down the Hamiltonian density
	\begin{align}
	\mathcal{H} = & \mu n_z + \frac{1}{2} \Lambda n^2_z + \frac{1}{2} \alpha_{xy} \left[\left(\nabla n_x \right)^2 + \left(\nabla n_y \right)^2\right]  + \frac{1}{2} \alpha_z \left(\nabla n_z \right)^2  + \vec{V}\!\left(x\right)\cdot \vec{n}\,.
	\end{align}
	The symmetry of model \eqref{eq:Hubbard_Hamiltonian} is represented by $U(1)$ rotational symmetry around $n_z$.
	
	We may solve for the collective modes and compute the scattering rates of these modes in the face of a disorder potential; this is worked out in Appendix \ref{app:CollectiveModeFieldTheory}. There we show that the low-$k$ scaling of the rates in channels $x$ and $z$ are $k^{D-3}$ and $k^{D+1}$, agreeing with the results obtained above numerically for $D = 1$. Crucially, the qualitative difference between the channels is that at small $k$ the CMs are predominantly phase modes, which do not couple directly to the disorder channel $z$, as it respects the $U(1)$ symmetry, while it may couple to channel $x$ directly. Instead, the only mechanism which couples the phase to channel $z$ is the modulation of the ground state amplitude field configuration, which shifts to screen the disorder. 
	This mechanism also applies to channel $x$, but it is dominated by the direct coupling mechanism. Therefore, $V^z$  manifests as a spatial modulation of the coefficients in the CM wave equation without generating additional terms. Scattering in channel $\lambda = z$ is thus effectively the CM equivalent of Rayleigh scattering and scales as $\lambda^4$ in three dimensions
	
	The above field theory also predicts that channel $y$ does not induce any scattering, which is a result of arbitrarily fixing the symmetry-broken order parameter to be real. Physically, it reflects the surprising fact that unlike modulations in the stationary order parameter amplitude, modulations in its stationary phase do not scatter the collective modes. 
	
	Finally, let us address channel $\lambda=0$. While the phenomenological model above does not include this channel explicitly, we note that the same mechanism applies: As channel $\lambda = 0$ respects the condensate symmetry, it would not introduce any symmetry-violating terms in the field theory action. It may therefore do one of three things: (i) Introduce higher-order, symmetry-respecting terms, which we are not interested in, (ii) effectively contribute to $V^z$ (but crucially, not to $V^x$ or $V^y$), or (iii) spatially modulate the model parameters.  {As discussed above, the latter two options are equivalent, and thus yield the same scaling.} 	Therefore, the last ingredient required is to know how $V^0$ couples into the action. As $V^0$ corresponds to an electrostatic potential, we know on physical grounds that a constant potential would not affect the system, so $V^0$ may not appear directly and the lowest-order scalar to consider is $\nabla^2 V^0$. This gives us the additional $(k^2)^2$ factor between the scaling relation of channels $\lambda = z$ and $\lambda=0$. {We  demonstrate this explicitly within our lattice model in Appendix \ref{app:ExcitonSelfEnergyDiagrammatics}.}
	
	
	{This field-theoretical description allows us to readily study  the case when the $U(1)$ symmetry is only approximate. The CM at $k = 0$ then gains a gap of a magnitude $\delta$ that depends on the degree of symmetry breaking. The long-wavelength dynamics are therefore affected, the CM group velocity now vanishing at $k=0$. For cases in which the excitonic description of the insulator is a good starting point, this symmetry violation is weak such that the gap is smaller that the two-particle continuum gap, $\delta \ll \Delta$. Since disorder only induces elastic scattering, the properties of the collective modes that we have derived should all apply for CMs in the energy window  $\delta \lesssim \omega \lesssim \Delta$.}
	
	
	\heading{Experimental consequences} 
	We now turn to discuss some of the possible experimental  consequences of the low-energy collective modes of the excitonic insulator. 
	Insights from the field theory above allow us to deduce both thermodynamic and transport properties of the system at low temperatures.
	
	Firstly, the collective branch below the transition supplies additional degrees of freedom, which should add to the system specific heat at $T$ below the gap. The specific heat can be computed similarly to that of acoustic phonons, and at low temperatures we expect the usual result of Debye {\cite{Kittel}},  $C_V \sim \left(T/\Theta_\mathrm{CM}\right)^D$ where $D$ is the dimension of the system. Here $\Theta_\mathrm{CM}$ is an ``excitonic Debye temperature'' of the order of the CM band width. 
	%
	{The relative size of this contribution compared to that of acoustic phonons with (phonon) Debye temperature $\Theta_\mathrm{D}$ is therefore of the order of $\left(\Theta_\mathrm{D}/\Theta_\mathrm{CM}\right)^D$. It is therefore likely to be most apparent in situations in which  the excitonic ordering temperature, and therefore $\Theta_\mathrm{CM}$ are small. However, the contribution can be comparable to the observed anomalies in specific heat seen in a candidate high-transition-temperature excitonic insulator \cite{Lu2017}.}
	
	
	Secondly, the collective modes have potentially very striking consequences for the transport properties of the EI. We have shown that at large wavelengths, it is possible for CMs to have large mean-free paths, if disorder in the system is symmetry-preserving. 
	These large mean free paths, combined with the electronic-scale group velocity, imply that they could have a dominant signature in the material thermal conductivity $\kappa \propto C_V v_g l_\mathrm{mfp}$. At temperatures below the particle-hole gap,
	one expects the excitonic insulator not just to transport heat via phonons, but to have a large contribution from the collective modes (at least at temperatures above any symmetry-breaking gap).
	
	If the only scattering mechanism present is disorder, the qualitative differences in the scattering of the CMs by the various disorder channels manifest as a stark differences in the resulting thermal conductivity, which we compute in Appendix \ref{app:CollectiveModeThermodynamics}. For disorder that couples to the excitonic order parameter ($\lambda = x$), $l_\mathrm{mfp} \rightarrow 0 $ for vanishing $k$, and we obtain $\kappa \sim T^3$ in all dimensions. If only the symmetry-respecting channels are present ($\lambda =0,z$), then the scattering rate vanishes for $k \rightarrow 0$. Similar to the treatment of acoustic phonons \cite{Kittel}, we can truncate the mean free path with the system linear size $L$. Therefore, the conductivity follows the specific heat and $\kappa \sim L T^D $ at very low temperatures, before plateauing once all the system-size-limited CMs are saturated. 
	
	In practice, as the temperature is increased, CMs will collide with thermally excited particles such as electrons, phonons, and each other. The separation of energy scales between CMs and phonons could  give rise to an interesting signature. To demonstrate this, we consider the CMs to be coupled to optical phonons with Einstein temperature $\Theta_\mathrm{E}$. The phonons are slower degrees of freedom, and on electronic timescales could be considered as static disorder fields. If the CM--phonon coupling mechanism breaks the $U(1)$ symmetry, then this scattering process would dominate all others once a thermal phonon population is excited at temperatures above $\Theta_\mathrm{E}$. However, below this temperature, mean free paths are limited by the system size and the thermal conductivity would increase. A schematic plot in such a system is depicted in Fig.~\ref{fig:thermal_conductivity}. 
	
	\begin{figure}[tb]
		\centering
	
		\includegraphics[width=\apsfigurewidth]{fig_6}
		\caption{CM thermal conductivity as a function of temperature, when coupled to optical phonons. The vertical dashed line  marks the phonon Einstein temperature. Below it, phonon thermal excitation is suppressed, and the CM mean free path becomes dominated by the system size.}
		\label{fig:thermal_conductivity}
	\end{figure}
	
	{Measurements of thermal conductivities of EI candidates are scarce. However, in \ce{TmSe_{0.45}Te_{0.55}} the suspected EI transition can be driven by varying pressure \cite{Wachter2004}, admitting comparison between the normal and EI phases at fixed temperatures. Indeed, below $20~\mathrm{K}$ the EI phase exhibits thermal conductivity orders of magnitude larger than that in the normal phase. Additionally, in the EI phase the thermal conductivity increases with decreasing temperature.}	This anomalous behaviour is attributed to superfluidity in the condensate \cite{Wachter2004, Wachter2005}. 
	However, it might be explained by a scenario such as that in Fig.~\ref{fig:thermal_conductivity}, in which decreasing $T$ reduces the scattering rate, thereby increasing thermal conductivity without invoking superfluidity. 	
	
	Finally, we point out that the robust collective modes could have other striking properties in dynamical measurements. The EI transition in candidate materials is often accompanied by other transitions such as a charge density wave \cite{Kogar2017} or a structural phase transition \cite{Kaneko2013}. This implies that the condensate could be coupled to other excitations such as plasmons or phonons. These then provide a more direct probe to study the EI collective modes by hybridizing with them.
	%
	In particular, we consider the case where CMs hybridize with phonons. \revised{This case was explored extensively in a recent study \cite{Murakami2020Phonons}, where it was shown that hybridization with phonons would produce a CM signature in the system optical conductivity. Retaining our focus on transport properties, we note that} the hybridized mode would gain a group velocity on the electronic scale, allowing lattice distortions to propagate at velocities much faster than the normal speed of sound, and this could be measured in response to mechanical perturbations. 
	Beyond this, the robustness of the CMs to disorder scattering, compounded by their high velocity, could  impart a long ballistic propagation length to the hybridized mode. The phonon-CM hybrid may then be more mobile then the charge carriers in the material, leading to long-range lattice dynamics and correlations which do not appear in the normal phase. This could explain recent spatially-resolved pump-probe experiments on \ce{Ta2NiSe5}{} \cite{Andrich2020}.
	
	
	
	\heading{Summary}
	In this work we investigated the transport properties of collective modes in the excitonic insulator phase. We have developed a method to calculate the effects of disorder on the collective modes from a microscopic model. We have used this to demonstrate that the collective modes are robust against scattering by disorder which respects the symmetry that is spontaneously broken in the condensate. We showed that this is not a feature unique to our particular model, but is rather  a general property of the relationship between the collective modes and the condensate symmetry, by reproducing our results in a phenomenological effective field theory. We predict that collective modes could have striking consequences for the thermal conductivity, and could be probed by mechanical response. We suggest that this could explain recent spatio-temporal measurements in \ce{Ta2NiSe5}{}.
	
	We thank P. Andrich, H. Bretscher and A. Rao for stimulating discussions about their experimental work, and to
	D. Gole{\v{z}}, Y. Murakami, D. Bennett, O. Hart, M. McGinley, and A. Szabo for insightful remarks.
	This work was supported by EPSRC Grant No. EP/P034616/1 and by a Simons Investigator Award.
	B.~R.~ also gratefully acknowledges  the support of the Cambridge International Trust and Wolfson College, Cambridge.
	
	
	\bibliographystyle{apsrev4-1}
	