\section{How to Use this Template}

\section{Introduction}
The liquid argon time projection chamber (LArTPC) detector technology provides high resolution event images and excellent calorimetric resolution for particle identification. The charged particles produce ionization electrons and scintillation light when they traverse liquid argon. The ionization electrons are drifted towards the anode planes under the electric field. In the single-phase design, the moving electrons produce current on the TPC wires at the anode. The wire signal is amplified and digitized by the electronics and then read out by the data acquisition system. The goal of the detector calibration is to convert the raw signal to the original energy deposition by removing the detector and physics effects, which we will discuss in detail in this paper. The calorimetric information is crucial for particle identification in a LArTPC, such as separating minimum ionizing particles (muons and charged pions) from highly ionizing particles (kaons and protons) and separating electrons from photons, which is the base for the neutrino cross section and oscillation measurements and the searches for beyond standard model physics. This paper summarizes the common techniques used to calibrate the TPC signal. The calibration of photon detector signal is not discussed in this paper.

\section{TPC Signal formation}
\label{sec:sigform}
The electron-ion pairs ($e^{-}$, Ar$^{+}$) are produced from energy loss by charged particles in liquid argon through the ionization process:
\begin{equation}
  N_{i} = \frac{\Delta E}{W_{\textrm{ion}}},
\end{equation}
where $N_{i}$ is the number of electron-ion pairs, $\Delta E$ is the energy loss, and $W_{\textrm{ion}} = 23.6 \pm 0.3$ eV~\cite{Miyajima:1974zz} is the ionization work function.

Some of the ionization electrons are recombined with surrounding molecular argon ions to form the excimer Ar$^{*}_{2}$. The excimers from both argon excitation and electron-ion recombination undergo dissociative decay to their ground state by emitting a vacuum ultraviolet photon. The free electrons that escape electron-ion recombination are drifted towards the wire planes under the electric field. Several effects can affect the electron drift. The electrons can be attached to contaminants in the liquid argon such as oxygen and water, which causes an attenuation of the signal on the TPC wires. The electron cloud is smeared both in the longitudinal and transverse directions by the diffusion effect. For a LArTPC located on the surface, there is a large flux of cosmic ray muons in the detector volume. Because of this, there is a significant accumulation of slow-moving argon ions (space charge) inside the detector, which leads to a position-dependent distortion to the electric field. The distorted electric field changes the electron-ion recombination rate and the trajectory of the free drifting electrons.

In a LArTPC with one grid plane and three instrumented wire planes, the drifting electrons first pass the grid plane, then pass two induction planes and eventually are collected by the wires on the collection plane. The electrons produce bipolar signals on the induction wires and unipolar signals on the collection wires. The wire signals is amplified by a preamplifier and then digitized by an analog-to-digital converter (ADC).
\section{LArTPC calibration}

In order to measure the energy loss per unit length ($dE/dx$) for particle identification, we need to convert the measured ADC counts to the energy deposition. The detector calibration procedure needs to remove all the instrument and physics effects following the reverse order of the signal formation described in Sec.~\ref{sec:sigform}. We discuss different ways to  measure each effect and practical methods to remove them.

\subsection{Electronics response and field response}
\label{sec:electronics}
Many LArTPC experiments have an electronics calibration system that has the capability to inject a known charge into each of the amplifiers connected to the TPC wires. In the ProtoDUNE pulser calibration system~\cite{Abi:2020mwi}, the injected charge is controlled by a six-bit voltage digital-to-analog converter (DAC):
\begin{equation}
  Q = SQ_{s},
\end{equation}
where $S$ is the DAC setting (0 - 63) and the step charge $Q_{s}$ = 3.43~fC~$\sim$ 21,400 electrons, which is roughly equivalent to the charge deposition of a minimum ionizing particle traveling parallel to the wire plane and perpendicular to that plane's wire direction on a single wire. The charge calibration is expressed as a gain for each channel normalized such that the product of the gain and the integral of the ADC counts over the pulse (pulse area) in a collection channel gives the charge in the pulse, i.e.~$Q = g A$. Figure~\ref{fig:pulser} shows the response of the ProtoDUNE electronics to an input charge of DAC setting 3.
\begin{figure}[htp!]
  \centering
  \includegraphics[width=0.70\textwidth]{figures/pulser.pdf}
\caption{
Response of the ProtoDUNE electronics to an input charge of DAC setting 3. The red curve is a fit to a simulated function of electronics response. Plot courtesy of David Adams (BNL) and the DUNE collaboration.
}
\label{fig:pulser}
\end{figure}

Figure~\ref{fig:tpcarea} shows the measured pulse area vs. DAC setting and
the fit for a typical collection channel.
The response is fairly linear over the DAC setting range (-5,~20) with saturation setting in
outside this range. Typical track charge deposits are one to four times the step charge and this saturation is only an issue for very heavily ionizing tracks. The gain for this channel is $g = (21.4$~ke)/(909.4~(ADC~count)-tick) =~23.5~e/((ADC~count)-tick).
\begin{figure}[htp!]
  \centering
  \includegraphics[width=0.70\textwidth]{figures/area_ch2080.pdf}
\caption{
Measured pulse area vs.\ DAC setting for a typical collection channel.
The red line shows the fit used to extract the gain. This figure is from~\cite{Abi:2020mwi}.
}
\label{fig:tpcarea}
\end{figure}

Figure~\ref{fig:tpcgains} shows the distribution of these gains for all channels.
Channels flagged as bad or especially noisy in an independent hand scan are shown
separately.
The gains for the remaining good channels are contained in a narrow peak with an
RMS of 5.1\% reflecting channel-to-channel response variation in the ADCs
and gain and shaping time variations in the amplifiers.
\begin{figure}[htp!]
\centering
  \includegraphics[width=0.70\textwidth]{figures/gaindist.pdf}
\caption{
Distribution of fitted gains for good (blue) and bad/noisy (red) channels.
The legend indicates the number of channels in each category and gives the mean (23.4~e/(ADC~count)/tick)) and RMS/mean (5.1\%) for the good channels. This figure is from~\cite{Abi:2020mwi}.
}
\label{fig:tpcgains}
\end{figure}

Ionization electrons follow the electric field lines as they pass through the wire planes producing direct and induced signals on nearby wires. Based on the Shockley-Ramo theorem, the instantaneous electric current $i$ on a particular electrode (wire) which is held at constant voltage, is given by
\begin{equation}
  i = e\nabla\phi \cdot \vec{v}_e,
\end{equation}
where $e$ is the charge in motion, and $\vec{v}_e$ is the
charge velocity at a given location, which is determined by the wire bias voltage settings. The weighting potential $\phi$
of a selected electrode at a given location is determined by virtually
removing the charge and setting the potential of the selected
electrode to unity while grounding all other conductors. For ProtoDUNE, the drift electric field and the weighting potential are calculated with Garfield~\cite{garfield}, as shown in Figure~\ref{fig:garfield}. The weighting field determines the bipolar signal shape on the induction plane and the unipolar signal shape on the collection plane.
\begin{figure}[htp!]
\centering
  \includegraphics[width=0.70\textwidth, trim={0cm 8cm 3cm 4.5cm},clip]{figures/VContour_trimmed.pdf}
\caption{
Garfield simulation of electron drift paths (yellow lines) in a 2D ProtoDUNE-SP TPC scheme and the equal weighting potential lines (green) for a given wire in the first induction plane, where the latter is shown in percentage from 1\% to 45\%. This figure is from~\cite{Abi:2020mwi}
}
\label{fig:garfield}
\end{figure}

The electronics response and the field response are normally removed from the raw wire signal through deconvolution~\cite{Baller:2017ugz}. If we ignore the charge induction on the neighboring wires, the time sampled wire signal read out by the data acquisition system, denoted $W(t)$, is the convolution of the ionization charge approaching the wire plane, $Q(t)$, with the field response $F(t)$, and the electronics response $E(t)$, plus the noise $N(t)$:
\begin{equation}
  W(t) = Q(t)*F(t)*E(t) + N(t).
\end{equation}
The time-dependent ionization charge can be recovered by deconvolution:
\begin{equation}
  Q^\prime(t) = \mathcal F^{-1}  \Bigg[ \frac{\Phi(f) \cdot {\mathcal F}(W(t))}{ {\mathcal F}(F(t)) \cdot {\mathcal F}(E(t))} \Bigg]
\end{equation}
where $\mathcal F$ is the Fourier Transform and $\Phi (f)$ is a filter function that maximizes the signal to noise ratio. The integral of filter function is usually normalized to 1 so it does not affect the charge reconstruction. The integral of electronics response is set to the gain measured by the pulser calibration system. For the collection wires, the integral of field response is normalized to 1. Therefore, the deconvolved signal $Q(t)$ is the ionization charge at the wire. For the induction wires, the reconstruction of charge is more difficult because of the bipolar signal shape.

A more advanced signal reconstruction technique is a two-dimensional (2D) deconvolution involving both the time and wire dimensions. The 2D deconvolution takes into account charge induction on the neighboring wires, which gives more accurate ionization charge information. It was first developed for the MicroBooNE experiment~\cite{Adams:2018dra} and then successfully used by the ProtoDUNE experiment~\cite{Abi:2020mwi}.

\subsection{Space Charge Effects}
\label{sec:sce}
In a large LArTPC located on the surface, space charge effects (SCE), i.e. the accumulation of positive argon ions ($Ar_{2}^{+}$) produced by the cosmic rays, distort the electric field significantly. Due to the ion mobility ($\mu_{i}\sim$ 10$^{-3} $cm$^{2}$V$^{-1}s^{-1}$) much smaller than the free electron one ($\mu_{e} \sim$ 500 cm$^{2}$V$^{-1}$s$^{-1}$), positive ions can survive in the drift region of the TPC for several minutes before being neutralized on the cathode or on the field shaping electrodes. In Ref.~\cite{Palestini:1998an}, the authors calculated the electric field $\mathscr{E}$ and space charge $\rho^{+}$, which we briefly summarize here. In the simplified model, $\mathscr{E}$ and $\rho^{+}$ are determined by the Maxwell and charge continuity equations:
\begin{equation}
  \frac{\partial\mathscr{E}}{\partial x} = \frac{\rho^{+}}{\varepsilon_{r}\varepsilon_{0}},
  \label{eq:maxwell}
\end{equation}
\begin{equation}
  \frac{\partial\rho^{+}}{\partial t} + \frac{\partial(\rho^{+} v_{i})}{\partial x} = J,
  \label{eq:continuity}
\end{equation}
where $x$ is the drift coordinate normal to the wire planes ($x = 0$ and $x = D$ define the anode and the cathode positions, respectively, and $D$ is the distance between anode and cathode), $\varepsilon_{0} = 8.854$ pF/m is the permittivity of vacuum, $\varepsilon_{r} = 1.504$ is the relative permittivity (dielectric constant) of liquid argon. The parameter $J$ is the average injected charge, which depends on the cosmic ray rate and the electric field applied. The ICARUS experiment measured $J = (1.9\pm0.1)\times10^{-10}$ C m$^{-3}$ s$^{-1}$ at the nominal electric field of 500 V/cm~\cite{Antonello:2020qht}. The MicroBooNE experiment estimated $J = 1.6\times10^{-10}$ C m$^{-3}$ s$^{-1}$ at the lower nominal electric field of 273 V/cm~\cite{Abratenko:2020bbx}. The parameter $v_{i} = \mu_{i}\mathscr{E}$ is the ion drift velocity. There is an uncertainty in the value of the mobility of positive argon ions $\mu_{i}$ (see references in~\cite{Palestini:2020dhv}). The measured mobility depends on several experimental conditions, such as the impurity level and temperature. A recent ICARUS measurement suggested $\mu_{i}$ is consistent with $0.9\times10^{-3}$ cm$^{2}$V$^{-1}$s$^{-1}$. MicroBooNE showed a reasonable agreement between the model prediction and the measured electric field~\cite{Abratenko:2020bbx} if a $\mu_{i}$ of $1.5\times10^{-3}$ cm$^{2}$V$^{-1}$s$^{-1}$ as reported in~\cite{doi:10.1063/1.334552} is used.

Eq.~\ref{eq:maxwell} and eq.~\ref{eq:continuity} are a one-dimensional approximation, which assumes symmetry in y and z coordinates. Introducing the dimensionless variable $\alpha$:
\begin{equation}
  \alpha = \frac{D}{\mathscr{E}_{0}}\sqrt{\frac{J}{\varepsilon_{r}\varepsilon_{0}\mu_{i}}},
\end{equation}
where $\mathscr{E}_{0} = V/D$ is the nominal electric field in absence of space charge produced by the high voltage $V$, eq.~\ref{eq:maxwell} and eq.~\ref{eq:continuity} can be solved to give:
\begin{equation}
  \mathscr{E}(x) = \mathscr{E}_{0}\sqrt{\left(\frac{\mathscr{E}_{a}}{\mathscr{E}_{0}}\right)^2+\alpha^{2}\frac{x^{2}}{D^{2}}},
\end{equation}
\begin{equation}
  \rho^{+}(x) = \frac{Jx}{\mu_{i}\mathscr{E}(x)},
\end{equation}
where $\mathscr{E}_a$ denotes the field at anode and is determined by the integral:
\begin{equation}
  \int_{0}^{D} \mathscr{E}(x) \,dx = V.
\end{equation}

The space charge increases the electric field at the cathode ($\mathscr{E}_{a}$) and decreases the electric field at the anode ($\mathscr{E}_{c}$) as shown in Figure~\ref{fig:Er}. When $\alpha\geq2$, the electric field at anode drops to 0 and the TPC stops working.
\begin{figure}[htp!]
\centering
  \includegraphics[width=0.70\textwidth]{figures/Er.pdf}
\caption{
Normalized electric field at the anode $\mathscr{E}_{a}/\mathscr{E}_{0}$ and cathode $\mathscr{E}_{c}/\mathscr{E}_{0}$ as a function of the parameter $\alpha$, reproduced from Ref.~\cite{Palestini:1998an}.
}
\label{fig:Er}
\end{figure}

There are many other effects that could affect the electric field. If the liquid argon purity is low, the contribution of negative ions from electrons attaching to the electronegative contaminants such as oxygen and water is not negligible~\cite{Luo:2020itx}. The fluid flow caused by the recirculation of liquid argon through the filter materials changes the space charge distribution since the fluid flow velocity and the ion drift velocity are of the same order. More sophisticated calculation of the space charge effects is discussed in~\cite{Palestini:1998an, Palestini:2020dhv}.

We now compare the electric field distortion in the three LArTPCs located on the surface: ICARUS, MicroBooNE and ProtoDUNE. Table~\ref{tab:sce} shows the running conditions of the three experiments. The calculated electric field distortion is systematically higher than the measurements. The fact that MicroBooNE prefers a larger value of positive argon ion mobility $\mu_{i}$ could be due to the different running conditions (e.g. argon purity and electric field). A careful treatment of detector-specific SCE is important for calibration.
\begin{table}[htp!]
\caption{Running conditions of the ICARUS, MicroBooNE and ProtoDUNE experiments and calculated electric field distortion in comparison with the measurements.}
\centering
\begin{tabular}{c|ccc}
\toprule
& \textbf{ICARUS}	& \textbf{MicroBooNE} & \textbf{ProtoDUNE}\\
\midrule
HV (kV) & 75 & 70 & 180 \\
$D$ (m) & 1.5 & 2.56 & 3.6\\
$\mathscr{E}_{0}$ (V/cm) & 500 & 273 & 500 \\
$J$ ($10^{-10}$ C m$^{-3}$ s$^{-1}$) & 1.9 & 1.6 & 1.9 \\
$\mu_{i}$ ($10^{-3}$ cm$^{2}$V$^{-1}$s$^{-1}$) & 0.9 & 1.5 & 0.9 \\
$\alpha$ & 0.378 & 0.839 & 0.907\\
$\mathscr{E}_{c}/\mathscr{E}_{0} - 1$ (calculated) & 4.7\% & 21.6\% & 25.0\%\\
$\mathscr{E}_{c}/\mathscr{E}_{0} - 1$ (measured) & 4\%~\cite{Antonello:2020qht} & 16\%~\cite{Abratenko:2020bbx, Adams:2019qrr} & 19\% $\pm$ 4\%~\cite{Abi:2020mwi}\\
$1 - \mathscr{E}_{a}/\mathscr{E}_{0}$ (calculated) & 2.4\% & 11.9\% & 14.0\%\\
$1 - \mathscr{E}_{a}/\mathscr{E}_{0}$ (measured) & 2\%~\cite{Antonello:2020qht} & 8\%~\cite{Abratenko:2020bbx, Adams:2019qrr} & 11\% $\pm$ 2\%~\cite{Abi:2020mwi}\\
\bottomrule
\end{tabular}
\label{tab:sce}
\end{table}

Different schemes to measure and correct for the space charge effects were developed for the MicroBooNE~\cite{Abratenko:2020bbx, Adams:2019qrr} and ProtoDUNE~\cite{Abi:2020mwi} experiments. One first measures the spatial distortions using either the UV laser tracks~\cite{Adams:2019qrr} or the cosmic ray muon tracks~\cite{Abratenko:2020bbx, Abi:2020mwi}. Without the electric field distortion, the reconstructed laser track is straight in the liquid argon. By comparing the reconstructed track trajectory in the distorted electric field with the true laser trajectory, one can measure the spatial distortions. Even though the muon tracks are not straight because of multiple Coulomb scattering in the absence of space charge effects, one can still measure the average spatial distortion using many muon tracks. The most straightforward method is to use two tracks that nearly cross. The comparison of the true crossing point with the reconstructed crossing point gives the spatial distortion at that particular position. The crossing-track method was employed in Ref.~\cite{Abratenko:2020bbx}. It was not practical to use this method with MicroBooNE's laser system because the set of locations in the detector where two laser beams can nearly cross is limited due to the limited motion of the reflecting mirrors. Instead, the authors developed a iterative method to measure the spatial distortions using single laser tracks~\cite{Adams:2019qrr}. In the first ProtoDUNE SCE measurement, the spatial distortions were measured at the front, back, top and bottom surfaces of the TPC. The spatial distortions in the bulk of the TPC volume were obtained through interpolation of measurements at the surfaces and scaling of simulated spatial distortions~\cite{Abi:2020mwi}.

Once the spatial distortion map is determined throughout the TPC volume, the electric field distortions associated with SCE can be computed. First, the local electron drift velocity is calculated from the spatial distortion map~\cite{Abi:2020mwi}. Once the local drift velocity $v(x, y, z)$ is determined throughout the TPC volume, the electric field magnitude is obtained by using the relationship between the electric field and the drift velocity, which is a function of the liquid argon temperature.

Both the spatial distortion map and the electric field map are used to correct the reconstructed track trajectory and the calorimetric measurement. The SCE affects the charge and energy loss per unit length ($dQ/dx$ and $dE/dx$) measurements in two ways. First, the spatial distortion changes the effective pitch $dx$ between two adjacent track hits. One can correct for this effect by recalculating the pitch $dx$ after correcting the two adjacent track hits using the spatial distortion map. Secondly, the distorted electric field changes the electron-ion recombination rate. One can correct for this effect by using the measured local electric field when converting $dQ/dx$ to $dE/dx$ using the recombination equation. More details on the recombination calibration will be discussed in Sec~\ref{sec:recomb}.

\subsection{Electron Attachment to Impurities}
In the presence of electronegative impurities, the concentration of free electrons $Q$ in liquid argon decreases exponentially with the drift time according to~\cite{Marchionni:2013tfa}
\begin{equation}
  \frac{dQ}{dt} = - k_{s}[S]Q,
\end{equation}
which leads to
\begin{equation}
  Q = Q_{0}e^{-k_{s}[S]t} = Q_{0}e^{-t/\tau_{e}},
\end{equation}
where $[S]$ is the concentration of electronegative impurities, $k_{s}$ is the electron attachment rate constant and $\tau_{e} = 1/(k_{s}[S])$ is the drift electron lifetime.

The liquid argon received from the supplier typically has contaminants of water, oxygen and nitrogen at the parts per million (ppm) level each. Water and oxygen capture drifting electrons and the concentration of these contaminants needs to be reduced by a factor of at least $10^{4}$ and maintained at this level to allow operation of the TPC. The purification system normally contains two filters~\cite{Benetti:1993mf, Curioni:2009rt, Adamowski:2014daa, Acciarri:2016smi, Abi:2020mwi}. The first filter that contains absorbent molecular sieve removes water contamination but can also remove small amounts of nitrogen and oxygen. The second filter that contains activated-copper-coated granules removes oxygen, and to a lesser extent, water. The nitrogen contaminant cannot be effectively filtered so the ultimate nitrogen concentration is set by the quality of the delivered argon. However, the nitrogen contaminant mainly suppresses the emission of scintillation photons and absorbs them as they propagate through the liquid argon~\cite{Jones:2013bca} and its impact on the drifting electrons is negligible.

The attachment rate constants to oxygen have been measured as a function of the electric field strength~\cite{doi:10.1021/j100564a006}. For electric fields of less than a few 100 V/cm, the attachment constant is measured to be $k_{O_{2}} = 9\times10^{10}\ \textrm{M}^{-1}\textrm{s}^{-1} = 3.6\ \textrm{ppm}^{-1} \mu s^{-1}$, which corresponds to $\tau_{e}(\mu s)\approx300/[\textrm{O}_{2}] (ppb)$. For increasing electric fields, the attachment rate constant decreases, as observed in several other measurements~\cite{Swan:1963dla, Bettini:1990rj}. At an electric field of 1 kV/cm, $\tau_{e}(\mu s)\approx500/[\textrm{O}_{2}] (ppb)$.

No direct measurement of attachment rate for water in liquid argon is reported in literature. A direct relation between the water concentration in the vapor above the liquid argon and the electron drift lifetime in the liquid argon is observed: (Drift Lifetime)$\times$(Water Concentration) = Constant, as reported in Ref.\cite{Andrews:2009zza}. In the same measurement, the authors concluded the concentration of water in the liquid argon effectively limits the drift electron lifetime.

The drift electron lifetime can be measured using the purity monitor~\cite{Carugno:1990kd, Amerio:2004ze, Amoruso:2004ti, Adamowski:2014daa, MicroBooNE:2016isc}. A purity monitor is a double gridded ion chamber immersed in the liquid argon volume, as shown in Figure~\ref{fig:prm}. The fraction of electrons generated via the photoelectric effect by the Xe lamp at the cathode that arrive at the anode $Q_{A}/Q_{C}$ after the electron drift time, $t$, is a measure of the electronegative impurity concentration and can be interpreted as the electron lifetime, $\tau_{e}$, such that
\begin{equation}
  Q_{A}/Q_{C} = e^{-t/\tau_{e}}.
\end{equation}
The SCE does not affect the operation of a purity monitor because the electrons are produced on the surface of the photocathode by the Xe lamp, and no ions are produced.
\begin{figure}[htp!]
\centering
  \includegraphics[width=0.70\textwidth]{figures/prm_1.pdf}
\caption{
A drawing of a liquid argon purity monitor. This figure is from~\cite{Adamowski:2014daa}.
}
\label{fig:prm}
\end{figure}

The drift electron lifetime can also be measured using cosmic ray muons recorded inside the TPC~\cite{Amoruso:2004ti, Bromberg:2015bqs, Anderson:2012vc, MicroBooNE:2017hgx, Acciarri:2019wgd, Abi:2020mwi}. The charge deposition per unit length measured by the muons is a function of the electron drift time for a constant drift electron lifetime:
\begin{equation}
  dQ/dx = (dQ/dx)_{0}e^{-t/\tau_{e}},
\end{equation}
where $(dQ/dx)_{0}$ is measured at the anode. For a LArTPC located on the surface, the SCE changes the measured $dQ/dx$ as discussed in Sec.~\ref{sec:sce}. Without correcting for the SCE, the measured electron lifetime is larger than the true electron lifetime because of the charge squeezing effect due to the spatial distortion and the change to the electron recombination due to the electric field distortion. ProtoDUNE developed a method to measure the drift electron lifetime using tracks tagged by the Cosmic-ray Tagging system (CRT). The selected TPC tracks are parallel to wire planes and perpendicular to the collection wire direction. Only the track segment in the central part of the TPC is used to minimize the spatial distortion caused by the SCE in the direction transverse to the nominal electric field. The change to the $dQ/dx$ measurement caused by the electric field distortion through recombination is small (less than 2\%) and corrected for using the measured electric field map. Figure~\ref{fig:lifetime} shows the lifetime measured in ProtoDUNE using CRT tagged muon tracks for two runs at different purity levels.
\begin{figure}[htp!]
\centering
  \includegraphics[width=0.49\textwidth]{figures/fit5759.pdf}
  \includegraphics[width=0.49\textwidth]{figures/fit5841.pdf}
\caption{
Plot the most probable value of the $dQ/dx$ distribution in ProtoDUNE as a function of the hit time, fit to an exponential decay function during a period of lower purity (left: $\tau_{e}$ = 10.4 ms) and during a period of higher purity (right: $\tau_{e}$ = 89 ms). This figure is from~\cite{Abi:2020mwi}.
}
\label{fig:lifetime}
\end{figure}

The two methods to measure the drift electron lifetime are complementary to each other. The purity monitors provide instant monitoring of the liquid argon purity but they measure purity in only a few locations outside the TPC. The muon method provides a direct measurement of the argon purity inside the TPC. The two methods may give different results of electron lifetime. The purity monitor normally operates at a much lower electric field ($\sim$ 25  V/cm) compared with the electric field in the TPC volume (500 V/cm). The rate constant for the attachment of electrons to the oxygen contaminant depends on the electric field, which means the measured lifetime using muons is higher than the one measured with the purity monitor if all contaminants are oxygen. The argon purity inside and outside the TPC can be different due to the fluid flow. The purity monitor measurement can be calibrated to provide the electron lifetime information inside the TPC, which is useful for the calibration of a LArTPC located deep underground (such as DUNE) where the cosmic flux is highly reduced.

\subsection{Diffusion}

During the drift, the ionization electron cloud spreads, which is an effect known as diffusion. The diffusion of electrons in strong electric fields is generally not isotropic. In general, the diffusion in the direction of the drift field (longitudinal diffusion) is smaller than the diffusion in the direction transverse to the field (transverse diffusion). The longitudinal and transverse sizes of an electron cloud are given by:
\begin{equation}
  \sigma_{L}^{2}(t) = \left(\frac{2D_{L}}{v_{d}^{2}}\right)t,
\end{equation}
\begin{equation}
  \sigma_{T}^{2}(t) = \left(\frac{2D_{T}}{v_{d}^{2}}\right)t,
\end{equation}
where $D_{L}$ and $D_{T}$ are longitudinal and transverse diffusion coefficients, $t$ is the electron drift time, and $v_d$ is the electron drift velocity. The diffusion coefficients are given by the Einstein-Smoluchowski relation~\cite{https://doi.org/10.1002/andp.19053220806, https://doi.org/10.1002/andp.19063261405}:
\begin{equation}
  D = \frac{\varepsilon_{e}}{e}\mu_{e},
\end{equation}
where $\varepsilon_{e}$ is the mean electron energy, $e$ is the electron charge, and $\mu_{e}$ is the electron mobility, which is a function of the electric field.
At relatively low electric field, the electrons gain so little energy from the field between the elastic atomic collisions that they come to thermal equilibrium with the medium. In this case, we have $\varepsilon_{e} = k\cdot T = 0.0075$ eV for argon at the nomal boiling point ($T$ = 87.3 K), where $k$ is the Boltzmann constant, and $\mu_{e} \approx 518$ cm$^{2}$/V/s~\cite{Li:2015rqa}. The corresponding coefficient is $D \sim 3.9$ cm$^{2}$/s, which is a lower limit to the diffusion coefficient. At strong electric fields the mean longitudinal electron energy increases while the electron mobility decreases. Therefore, the longitudinal diffusion coefficient is nearly a constant up to $\mathscr{E} = 1$ kV/cm. The ratio of the longitudinal to transverse diffusion coefficient is expressed as~\cite{Li:2015rqa}:
\begin{equation}
  \frac{D_{L}}{D_{T}} = 1 +\frac{\mathscr{E}}{\mu_{e}}\frac{\partial\mu_{e}}{\partial\mathscr{E}}.
\end{equation}
At low electric fields, the mobility is close to a constant, and we have $D_{L}/D_{T}\approx1$.

The longitudinal diffusion broadens the signal waveforms as a function of drift time while the transverse diffusion smears the signal between neighboring wires. A good understanding of the diffusion process is important because it can influence the event images and the accuracy of the drift coordinate measurement.

The ICARUS collaboration reported measurements of the longitudinal diffusion effect at four different electric fields: 100, 150, 250, 350 V/cm using a three-ton LArTPC with a maximum drift distance of 42 cm~\cite{Cennini:1994ha}. The longitudinal diffusion coefficient $D_{L}$ was derived from the analysis of the rise time of the signal on the collection wires. The resulting coefficients at different fields are consistent within the errors, which gives an average measurement of $D_{L} = 4.8\pm0.2$ cm$^{2}$/s. Li et al. from BNL reported measurements of $D_{L}$ between 100 and 2000 V/cm using a laser-pulsed gold photocathode with drift distance ranging from 5 to 60 mm~\cite{Li:2015rqa}. The ICARUS results show a good agreement with the prediction of Atrazhev and Timoshkin~\cite{689434}, while the BNL results are systematically higher than both.

Both MicroBooNE and ProtoDUNE are measuring the longitudinal diffusion coefficient $D_{L}$ by measuring the squared time width of a signal pulse, $\sigma_{t}^{2}$ as a function of drift time $t$:
\begin{equation}
  \sigma_{t}^{2}(t) \approx \sigma_{0}^{2} + \left(\frac{2D_{L}}{v_{d}^{2}}\right)t,
\end{equation}
where $\sigma_{0}$ is the signal width at anode, mostly determined by the electronics response discussed in Sec.~\ref{sec:electronics}. The results are expected to be published soon.

The diffusion process does not affect the charge reconstruction if hits are reconstructed with a variable time width and if the hit integral rather than the peak amplitude is used to measure the deposited charge since drift electrons are not lost by this process.

\subsection{Recombination}
\label{sec:recomb}
Electrons emitted by ionization are thermalized by interactions with the surrounding medium after which time they may recombine with nearby ions~\cite{Onsager:1938zz, https://doi.org/10.1002/andp.19133471205}. Electron-ion recombination introduces a non-linear relationship between $dE/dx$ and $dQ/dx$. The recombination effect can be measured using the stopping particles (muons or protons), which cover a wide range of $dE/dx$ values. For each point on a stopping track, $dE/dx$ is calculated from the distance to the track end (residual range)~\cite{Acciarri:2013met}, and $dQ/dx$ is calculated by converting the measured ADC counts to the number of electrons (discussed in ~\ref{sec:electronics}). The recombination effect can be parameterized by two empirical models. One is the Birks' model developed by the ICARUS collaboration~\cite{Amoruso:2004dy}:
\begin{equation}
\frac{dQ}{dx}(e/cm)=\frac{A_B}{W_{\textrm{ion}}}\left(\frac{\frac{dE}{dx}}{1+\frac{k_{B}}{\rho\mathscr{E}}\frac{dE}{dx}}\right),
\label{eqn:Birks_de_dx}
\end{equation}
where $A_{B}$ and $k_{B}$ are the two model parameters, $W_{\textrm{ion}} = 23.6$ eV is the ionization work function of arogn, $\mathscr{E}$ is the drift electric field, and $\rho$ is the liquid argon density. The other model is the modified Box model developed by the ArgoNeuT collaboration~\cite{Acciarri:2013met}:
\begin{equation}
\label{eqn:dq_dx_modbox}
\frac{dQ}{dx}(e/cm) =\frac{\ln(\frac{dE}{dx}\frac{\beta^{\prime}}{\rho\mathscr{E}}+\alpha)}{\frac{\beta^{\prime}}{\rho\mathscr{E}} W_{\textrm{ion}}},
\end{equation}
where $\alpha$ and $\beta^{\prime}$ are the two model parameters, and the other parameters are the same as in the Birks' model. In the detector calibration, $dE/dx$ is calculated by inverses of Eqs.~\ref{eqn:Birks_de_dx} and~\ref{eqn:dq_dx_modbox}:
\begin{equation}
  \frac{dE}{dx} = \frac{dQ/dx}{A_{B}/W_{\textrm{ion}}-k_{B}(dQ/dx)/(\rho\mathscr{E})},
  \label{eqn:de_dx_birks}
\end{equation}
or
\begin{equation}
  \frac{dE}{dx} = \frac{\rho\mathscr{E}}{\beta^{\prime}}(\exp(\beta^{\prime}W_{\textrm{ion}}(dQ/dx)/(\rho\mathscr{E}))-\alpha).
  \label{eqn:de_dx_box}
\end{equation}
The modified Box model describes data more accurately at higher $dQ/dx$. For $dQ/dx > (A_{B}\rho\mathscr{E})/(W_{\textrm{ion}}k_{B})$, $dE/dx < 0$ according to Eq.~\ref{eqn:de_dx_birks}. The LArIAT experiment combines both models for different ionization densities in the Michel electron measurement~\cite{Foreman:2019dzm}.

The recombination effects have been measured by the ICARUS~\cite{Amoruso:2004dy}, ArgoNeuT~\cite{Acciarri:2013met} and MicroBooNE~\cite{Adams:2019ssg} experiments at different electric fields. Table~\ref{tab:recomb} summarizes the recombination model parameters measured by the three experiments. We would like to point out that the two models with the best fit parameters are only valid in a given range of electric fields and $dE/dx$ values. Note the $k_{B}$ and $\beta^{\prime}$ parameters measured by the MicroBooNE experiment are lower than the ones measured by the other two experiments. Furthermore, the MicroBooNE recombination measurement was performed using an old version of simulation and before the SCE calibration was available. More details can be found in Ref.~\cite{Adams:2019ssg}.
\begin{table}[htp!]
\caption{The Birks and modified Box model parameters measured by the ICARUS, ArgoNeUT and MicroBooNE experiments}
\centering
\begin{tabular}{c|ccc}
\toprule
& \textbf{ICARUS}	& \textbf{ArgoNeuT} & \textbf{MicroBooNE}\\
\midrule
$\mathscr{E}$ (V/cm) & 200, 350, 500 & 481 & 273 \\
Sample & Stopping muons & Stopping protons & Stopping protons \\
Birks $A_{B}$ & $0.800 \pm 0.003$ & $0.806 \pm 0.010$ & $0.816 \pm 0.012$ \\
Birks $k_{B}$ ((kV/cm)(g/cm$^{2}$)/MeV) & $0.0486 \pm 0.0006$ & $0.052 \pm 0.001$ & $0.045 \pm 0.001$ \\
Box $\alpha$ & - & $0.93 \pm 0.02$ & $0.92 \pm 0.02$ \\
Box $\beta^{\prime}$ ((kV/cm)(g/cm$^{2}$)/MeV) & - & $0.212 \pm 0.002$ & $0.184 \pm 0.002$ \\
\bottomrule
\end{tabular}
\label{tab:recomb}
\end{table}

\subsection{Stopping Muon Calibration}
There are generally two steps in converting the measured raw ADC counts to energy. The first step is to convert $\frac{dQ}{dx}$ (ADC/cm) to $\frac{dQ}{dx}$ ($e$/cm). This step employs the measured electronics gain and removes effects that cause non-uniformities in the detector response, including the SCE and the free electron attenuation. The second step converts $\frac{dQ}{dx}$ ($e$/cm) to $\frac{dE}{dx}$ (MeV/cm) using either the Birks or modified Box recombination model. However, there could be uncertainties in each step of calibration that affect the final calibrated $dE/dx$ quantity. One can introduce a calibration constant to correct the calibrated charge. This calibration constant is a scaling factor that applies to $\frac{dQ}{dx}$ ($e$/cm), which can be fine-tuned using the stopping muons. The kinetic energy of the muon at each track hit can be determined by the residual range. A portion of the track can be selected that corresponds closely to a minimum ionizing particle, for which $dE/dx$ is very well understood theoretically to better than 1\%. The calibration constant ensures the measured $dE/dx$ agrees with the prediction made by the Landau-Vavilov function~\cite{Bichsel:1988if} in the MIP region. Figure~\ref{fig:cali_plots} shows the measured $dE/dx$ values in MicroBooNE data after tuning the calibration constant, which show a good agreement with the prediction in the MIP region. This universal calibration constant can be used to calibrate the detector response to all particles.
\begin{figure}[!ht]
\centering
\includegraphics[width=0.7\textwidth]{figures/dEdxuB.pdf}
\caption{Comparison between the predicted and the measured most probable value $dE/dx$ for stopping muons in MicroBooNE data. This figure is from~\cite{Adams:2019ssg}.}
\label{fig:cali_plots}
\end{figure}

In principle, if the electronics gains are correctly measured by the pulser calibration system and all the effects that cause non-uniformities in the detector response are properly corrected for, the calibration constant derived using stopping muons should be exactly 1. The goal of the LArTPC calibration is to understand all the detector and physics effects in order to achieve an accurate and robust calorimetric measurement.
\section{Conclusions}
In this article, we review the major effects that affect the calorimetric reconstruction in a LArTPC: electronics and field responses, space charge effects, electron attachment to impurities, diffusion and recombination. We also summarize the general methods to measure those effects and the results from different LArTPC experiments. Precise measurements of those effects improve the understanding of the detector response and energy resolution, which is crucial for reach physics sensitivity in the current and future LArTPC experiments.

\vspace{6pt}

\funding{Work supported by the Fermi National Accelerator Laboratory, managed and operated by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.}

\acknowledgments{The author thanks Bruce Baller, David Caratelli, Flavio Cavanna, Tom Junk, and Stephen Pordes from Fermilab for their careful reading of the paper and for their helpful comments.}

\externalbibliography{yes}


\section{Introduction}\label{Introduction}
In this paper, we study a hyperbolic model
for the dynamics of relativistic viscous fluids
and confirm its stability and causality in terms of decay and subluminal speed of its Fourier-Laplace modes.
The result is complete in the sense that it covers modes of arbitrary wave number and with respect
to arbitrary Lorentz frames.
Concretely, we address the system
\begin{eqnarray}\label{Teq}
{\partial \over \partial x^\beta}\left(T^{\alpha\beta}
+\Delta T^{\alpha\beta}_\Box
\right)
=0
\end{eqnarray}
which was recently proposed
as a candidate for a relativistic version of Navier-Stokes \cite{FT1,FT3}.
While in \eqref{Teq}
\begin{eqnarray} \label{T}
T^{\alpha\beta}=\rho U^\alpha U^\beta+p \Pi^{\alpha\beta}
\end{eqnarray}
is the usual ideal-fluid energy-momentum tensor, with $\rho,p, U^\alpha$ the fluid's energy,
pressure and 4-velocity, and $\Pi^{\alpha\beta}=g^{\alpha\beta}+U^\alpha U^\beta$ the projector
on the orthogonal complement of $U^\alpha$, the dissipation tensor reads
\begin{eqnarray}\label{DeltaT}
\begin{aligned}
-\Delta T^{\alpha\beta}_\Box
&=&\eta \Pi^{\alpha\gamma}\Pi^{\beta\delta}
\left[\frac{\partial U_{\gamma}}{\partial x^{\delta}}
+\frac{\partial U_{\delta}}{\partial x^{\gamma}}
-\frac{2}{3}g_{\gamma\delta}\frac{\partial U^{\epsilon}}{\partial x^{\epsilon}}\right]
+\check\zeta \Pi^{\alpha\beta}\frac{\partial U^{\gamma}}{\partial x^{\gamma}}\\
&&+\sigma\left[U^{\alpha}U^{\beta}\frac{\partial U^{\gamma}}{\partial x^{\gamma}}
-\left(\Pi^{\alpha\gamma}U^{\beta}
+\Pi^{\beta\gamma}U^{\alpha}\right)
U^{\delta}\frac{\partial U_{\gamma}}{\partial x^{\delta}}\right]
\\
&&+
\chi
\bigg[
\bigg(U^\alpha\frac{\partial \theta}{\partial x_\beta}
+U^\beta\frac{\partial \theta}{\partial x_\alpha}\bigg)
-g^{\alpha\beta}U^\gamma\frac{\partial \theta}{\partial x^\gamma}
\bigg],
\end{aligned}
\end{eqnarray}
where
\begin{eqnarray}\label{choicetildezeta}
\sigma=((4/3)\eta+\zeta)/(1-c_s^2)-c_s^2\chi\theta
\quad\text{and}\quad
\check\zeta=\zeta+c_s^2\sigma-c_s^2(1-c_s^2)\chi\theta
\end{eqnarray}
with $\eta,\zeta,\chi$ the coefficients of shear viscosity, bulk viscosity, and thermal conductivity,
and $0<c_s<1$ the speed of sound.
We refer the interested reader to \cite{FT3} for a detailed derivation and justification of
equations \eqref{Teq} with \eqref{DeltaT}, \eqref{choicetildezeta},
and to \cite{FT1,FT2,F20} for brief summaries of the history and current developments of
the mathematical theory of relativistic viscous fluid flow.  As in \cite{FT3} we will assume that
the fluid is barotropic, i.e., pressure and sound speed are functions of the
energy,
\begin{eqnarray}\label{P}
p=P(\rho),\quad c_s^2=P'(\rho),
\end{eqnarray}
which implies that the local state of the fluid has four degrees of freedom (three for the velocity,
which is constrained by $U_\alpha U^\alpha=-1$, and one for the thermodynamics) and correspondingly
the conservation laws \eqref{Teq} for energy-momentum suffice to determine the dynamics.

We focus here on Fourier-Laplace modes. These are solutions of the form
\begin{eqnarray}\label{FLmode}
\hat\psi_\alpha
e^{\lambda t+i\boldsymbol\xi\cdot \boldsymbol x}
\quad (\lambda,\boldsymbol\xi)\in \mathbb C\times\mathbb R^3,
\end{eqnarray}
of the linearization of \eqref{Teq} at some homogeneous reference state, with $\boldsymbol x$
referring to the spatial variables $x^1,x^2,x^3$  and $t$ to time $x^0$.
In relativity, the notion of a Fourier-Laplace mode with
its distinction between spatial and temporal frequencies is relative to the choice of a Lorentz
frame.
A non-trivial mode, \eqref{FLmode} with $\boldsymbol\xi\neq 0$, is called
(a) decaying, neutral, or growing if
$\hbox{Re}\{\lambda\}<0,=0,$ or $>0$, respectively, and
(b) subluminal, luminal, or superluminal if $|\hbox{Im}\{\lambda\}|<|\boldsymbol\xi|,=|\boldsymbol\xi|$, or
$>|\boldsymbol\xi|$, respectively.

\begin{theo}\label{ThmD}
For any values
\begin{eqnarray}\label{etahatzeta}
\eta>0,\quad \check\zeta\ge -\frac13\eta,\quad \sigma=\frac43\eta+\check\zeta,\quad\chi\ge 0,
\end{eqnarray}
the viscous system \eqref{Teq} with \eqref{DeltaT}--\eqref{P}
is dissipative in the sense that all its non-trivial Fourier-Laplace modes with respect to any
Lorentz frame are decaying.
\end{theo}
\begin{theo}\label{ThmC}
For system \eqref{Teq} with \eqref{DeltaT}--\eqref{P}
all non-trivial Fourier-Laplace modes with respect to any Lorentz frame are subluminal.
\end{theo}

Theorem 1 is proved in Section 2, Theorem 2 in Sections 3 and 4.
In Section 5, we address a counterpart of Theorem 1 for non-barotropic fluids; the corresponding
model, taken from \cite{FT2}, is introduced there.

The interpretation of Theorem 1 is clear -- it gives linear stability of homogeneous states in
all reference frames.  Rigorous stability results at the fully nonlinear level have been
accomplished by  Sroczinski \cite{S1,S2}.

Theorem 2  is related to  causality. Regarding modes, causality means that
$\lim_{|\boldsymbol\xi|\to\infty}(|\lambda(\boldsymbol\xi)|/|\boldsymbol\xi|)$ $\le 1$, a property that for
the models considered in the present paper has been checked in \cite{FT2,FT3}.
Our theorem proves a finer property, namely the fact that all \emph{phase velocities},
with $0<|\boldsymbol\xi|<\infty$,  are less
than the speed of light. Recall the recurrent considerations in the literature on interesting physical
systems which do not have this property (see \cite{B} for an early documentation).
For our models, the property is \emph{plausible}, as
both the second-order operator (corresponding to ``$|\boldsymbol\xi|=\infty$'') and the first order operator
(corresponding to ``$|\boldsymbol\xi|= 0$'') have  (sub-/)luminal phase velocity. Our \emph{proof} suggests the possibility of a more general principle granting
such ``interpolation''
between hyperbolic operators of different orders,
but the potential identification of such principle is here left to future work.

\section{Modes and stability in arbitrary frames}
\setcounter{equation}{0}
We first note that a solution to the linearized PDE which is a Fourier-Laplace mode with respect
to a given
frame is generally not a Fourier-Laplace mode with respect to a different frame. We thus briefly
consider \emph{general modes}, i.e., solutions of the linearized
equations of the form
\begin{eqnarray}\label{mode}
{{\hat\psi_\alpha}}
e^{i\xi^\beta  x_\beta},
\quad \xi^\beta\in\mathbb C^{4},
\end{eqnarray}
with, in principle, arbitrary \emph{complex} wave vectors. This \emph{is} a covariant notion;
under frame changes, the amplitude $\hat\psi_\alpha$ and notably the real part $\text{\rm Re }\xi^\beta$
and the imaginary part $\text{\rm Im }\xi^\beta$ of the wave vector transform as 4-vectors.
A general mode \eqref{mode} is a Fourier-Laplace
mode with respect to \emph{some} Lorentz frame if and only if
\begin{eqnarray}    \label{FLM_invariant_form}
\text{\rm Im }\xi^\beta \text{ is timelike or zero,}
\end{eqnarray}
in the latter case it is \emph{neutral}.
 It is a Fourier-Laplace mode with respect to the fluid rest frame,
if $\text{\rm Im }\xi^\beta$ is a multiple of the future-pointing fluid $4$-velocity
$U^\beta$ of the homogeneous reference state. A Fourier-Laplace mode is decaying if
$\text{\rm Im }\xi^\beta$ is future-pointing, and growing if $\text{\rm Im }\xi^\beta$ is past-pointing.
A Fourier-Laplace mode is \emph{subluminal} or \emph{luminal} if
\begin{eqnarray}
\text{\rm Re }\xi^\beta \text{ is spacelike or lightlike}.
\end{eqnarray}

We will study solutions that are Fourier-Laplace modes with respect to \emph{some} frame
by looking at them in the fluid's rest frame. To prepare for that perspective,
consider now an arbitrary general mode from the point of view of its
representation in the fluid's rest frame, let this representation be given by
\begin{eqnarray}\label{genmodeinrestframe}
\hat\psi_\alpha
e^{\lambda t+i\boldsymbol\xi\cdot \boldsymbol x},
\quad (\lambda,\boldsymbol\xi)\in \mathbb C\times\mathbb C^3;
\end{eqnarray}
i.e., relative to \eqref{mode} we have written $\boldsymbol\xi$ for $(\xi^1,\xi^2,\xi^3)$ and $\lambda$
for $i\xi^0$;
 the only difference of \eqref{genmodeinrestframe} from \eqref{FLmode} is that
$\boldsymbol\xi$ now is complex.
Assume that with respect to a fixed ``rest frame``, the dispersion relation for a covariant linear
system of partial differential equations is given by
\begin{eqnarray}\label{drrf}
0=\Pi(\lambda,\boldsymbol\xi)
\end{eqnarray}
i.e., $\Pi$ results from the symbol of the differential operator describing the system
with respect to the rest frame.

Consider now a solution to the PDE system
which is a Fourier-Laplace mode
\begin{eqnarray}\label{tildemode}
\widetilde{\hat\psi_\alpha}
e^{\tilde\lambda \tilde t+i\tilde\boldsymbol\xi\cdot\tilde{\boldsymbol x}}
=
\widetilde{{\hat\psi}_\alpha}
e^{Re\{\tilde\lambda\}\tilde t}e^{i\left(Im\{\tilde\lambda\}t+\tilde\boldsymbol\xi\cdot \tilde{\boldsymbol x} \right)},
\quad (\tilde\lambda,\tilde\boldsymbol\xi)\in \mathbb C\times\mathbb R^3,
\end{eqnarray}
with respect to a frame which is not the rest frame, but differs from it by a Lorentz boost.
In order to nevertheless study it by means
of the rest-frame dispersion relation \eqref{drrf} we note (cf.\ also \cite{HL}) that
if \eqref{FLmode} represents the same mode with respect to the rest frame,  the complex wave
vectors
$(\tilde\lambda,i\tilde\boldsymbol\xi)$ and $(\lambda,i\boldsymbol\xi)$ are related
through the Lorentz transformation
\begin{eqnarray}
(\lambda,\boldsymbol\xi)=\Lambda_v(\tilde\lambda,\tilde\boldsymbol\xi),
\end{eqnarray}
from the non-rest frame to the rest frame, i.e.,
\begin{eqnarray}\label{Lorentz_var}
\begin{aligned}
\lambda&=\gamma\tilde\lambda+i\gamma v\tilde\xi_1 \\
\xi_1&=-i\gamma v\tilde\lambda+\gamma\tilde\xi_1 \\
\xi_2&=\tilde\xi_2\\
\xi_3&=\tilde\xi_3,
\end{aligned}
\end{eqnarray}
with some $$v\in(-1,1)\text{ and }\gamma=(1-v^2)^{-1/2}.$$
The Fourier-Laplace modes in the tilde
frame are thus controlled by the dispersion relation
\begin{eqnarray} \label{disp_rel}
0=\tilde\Pi_v(\tilde\lambda,\tilde\boldsymbol\xi)
\end{eqnarray}
with
$$
\tilde\Pi_v(\tilde\lambda,\tilde\boldsymbol\xi)
\equiv
\Pi(\Lambda_v(\tilde\lambda,\tilde\boldsymbol\xi))=\Pi(\lambda,\boldsymbol\xi).
$$

In our specific situation, the linearized equations of motion in the fluid's rest frame,
\begin{eqnarray}
\begin{aligned}\label{PDElin}
\frac{\rho+p}\theta\left\{c_s^{-2}
{\partial \theta\over \partial t}
+\theta\nabla\cdot{\bf u}
\right\}
+\left\{\chi\left(\frac{\partial^2\theta}{\partial t^2}-\Delta\theta\right)\right\}
&=0,\\
\frac{\rho+p}\theta\left\{
\theta{\partial {\bf u}\over \partial t}
+\nabla\theta
\right\}
+\left\{{\sigma}\frac{\partial^2{\bf u}}{\partial t^2}
-({\eta}\,\nabla\cdot{\cal S}{\bf u}+{\check\zeta}\,\nabla(\nabla\cdot {\bf u}))\right\}&=0
\end{aligned}
\end{eqnarray}
yield, after a simple scaling, the dispersion relation
\begin{eqnarray}\label{dr}
0=\Pi(\lambda,\boldsymbol\xi)\equiv \det M(\lambda,\boldsymbol\xi),
\end{eqnarray}
with
(cf.\ \cite{FT3})
$$
M(\lambda,\boldsymbol\xi)
=
\begin{pmatrix}
c_s^{-2}\lambda  & i\boldsymbol\xi^\top \\
i\boldsymbol\xi & \lambda I
\end{pmatrix}
+
\begin{pmatrix}
\chi(\lambda^2+\boldsymbol\xi^2)
   & 0 \\
0&  N(\lambda,\boldsymbol\xi)
\end{pmatrix},
$$
where
$$
N(\lambda,\boldsymbol\xi)=\left(\sigma\lambda^2+\eta\boldsymbol\xi^2\right)I
+\left(\check\zeta+{1\over 3}\eta\right)\boldsymbol\xi\bxi^\top,
$$
with $\boldsymbol\xi^2\equiv\boldsymbol\xi^\top\boldsymbol\xi=\xi_1^2+\xi_2^2+\xi_3^2,$
and the dispersion relation factors as
\begin{eqnarray} \label{drsplit}
\Pi(\lambda,\boldsymbol\xi)=\Pi^L(\lambda,\boldsymbol\xi) \big(\Pi^T(\lambda,\boldsymbol\xi) \big)^2
\end{eqnarray}
with
$$
\Pi^L(\lambda,\boldsymbol\xi)=\{c_s^{-2}\lambda +\chi (\lambda^2+\boldsymbol\xi^2)\}
\{\lambda + \sigma\left(\lambda^2+\boldsymbol\xi^2\right)\}
+ \boldsymbol\xi^2
$$
and
$$
\Pi^T(\lambda,\boldsymbol\xi)=\lambda+\sigma\lambda^2+\eta\boldsymbol\xi^2;
$$
to simplify notation, we will use the again scaled versions
\begin{eqnarray}
\Pi^L(\lambda,\boldsymbol\xi)=\{\lambda +\hat\chi (\lambda^2+\boldsymbol\xi^2)\}
\{\lambda +\left(\lambda^2+\boldsymbol\xi^2\right)\}
+ c_s^2 \boldsymbol\xi^2
\end{eqnarray}
and
\begin{eqnarray}
\Pi^T(\lambda,\boldsymbol\xi)
=\lambda+\hat\sigma\lambda^2+\boldsymbol\xi^2
\end{eqnarray}
with $\hat\chi=c_s^2\chi/\sigma\ge 0$ and $\hat\sigma=\sigma/\eta\ge 1$.

Correspondingly, the Lorentz transformed dispersion relation
\eqref{disp_rel} also factors, as
\begin{eqnarray}\label{drtildesplit}
\tilde\Pi_v(\tilde\lambda,\tilde\boldsymbol\xi) =
\tilde\Pi^L_v(\tilde\lambda,\tilde\boldsymbol\xi)  \big(\tilde\Pi^T_v(\tilde\lambda,\tilde\boldsymbol\xi)\big)^2,
\end{eqnarray}
with
\begin{eqnarray} \label{tildePiLtildePiT}
\tilde\Pi^L_v(\tilde\lambda,\tilde\boldsymbol\xi)
\equiv
\Pi^L(\lambda,\boldsymbol\xi),\quad
\tilde\Pi^T_v(\tilde\lambda,\tilde\boldsymbol\xi)
\equiv
\Pi^T(\lambda,\boldsymbol\xi).
\end{eqnarray}
Concretely, setting $r^2\equiv \xi_2^2 +\xi_3^2 = \tilde\xi_2^2 + \tilde\xi_3^2$
and substituting the identities
$$
\lambda^2 + \boldsymbol\xi^2 = \tilde\lambda^2 + \tilde\boldsymbol\xi^2
\hspace{1cm} \text{and} \hspace{1cm}
\boldsymbol\xi^2 = r^2 - \gamma^2 (\tilde\lambda v + i \tilde\xi_1)^2,
$$
into $\tilde\Pi^L$, we obtain
\begin{eqnarray}  \label{Pi_L_gen2}
\tilde\Pi^L(\tilde\lambda,\tilde\boldsymbol\xi)
& = & \hat\chi  \big(\tilde\lambda^2+\tilde\boldsymbol\xi^2\big)^2 + \gamma (\hat\chi +1) \big(\tilde\lambda + i v \tilde\xi_1\big)\big(\tilde\lambda^2+\tilde\boldsymbol\xi^2\big) + \gamma^2 \big(\tilde\lambda + i v \tilde\xi_1\big)^2  \cr  & &    +\; c_s^2\big( r^2 - \gamma^2 (\tilde\lambda v + i \tilde\xi_1)^2 \big),
\end{eqnarray}
and a straightforward computation gives
\begin{align} \label{Pi_T_gen2}
\tilde\Pi^T(\tilde\lambda,\tilde\boldsymbol\xi)
&=\gamma(\tilde\lambda+iv\tilde\xi_1) +\hat\sigma(\gamma\tilde\lambda+i\gamma v\tilde\xi_1)^2 + (-i\gamma v\tilde\lambda+\gamma\tilde\xi_1)^2+r^2.
\end{align}
We now turn to the proof of Theorem 1, which expresses  a general principle
to conclude stability in all frames from stability in one frame, based on the idea
to bound Re$(\lambda)$ away from 0 by showing that neutral modes do not exist.\footnote{This idea
seems to have been brought up first in \cite{FT1}, Lemmas 6.4 and 6.6.}
\begin{prop}
There does not exist a nontrivial Fourier-Laplace mode of \eqref{PDElin}
in any Lorentz frame, which is neutral.
\end{prop}
\begin{proof}
While the question of whether a given \emph{general} mode \eqref{mode} is a Fourier-Laplace mode
with respect to a Lorentz frame generically depends on the choice of the frame,
a  mode that is a \emph{neutral} Fourier-Laplace mode with respect to one Lorentz frame, obviously
is a neutral Fourier-Laplace mode with respect to any other Lorentz frame as well. Thus if a
neutral nontrivial Fourier-Laplace mode with respect to some frame existed, it would readily be
a neutral nontrivial Fourier-Laplace mode with respect to the fluid's rest frame. But it is
known from \cite{FT3} that with respect to the fluid's rest frame, the model under consideration
allows only decaying, thus no neutral nontrivial Fourier-Laplace modes.
\end{proof}

{ \it Proof of Theorem 1.}
Assume that with respect to some frame, there is a non-trivial
Fourier-Laplace mode that does not decay, briefly: for some $\tilde v\in(-1,1)\setminus\{0\}$
there exists
$\tilde\boldsymbol\xi\in\mathbb R^3\setminus\{0\}$ such that $\tilde\Pi_{\tilde v}(.,\tilde\boldsymbol\xi)$ has a zero
$\tilde\lambda$ with positive real part. However, the set of (complex) zeroes of a continuously
parameterized family of (complex) polynomials depends continuously on the parameter. Using this
for the family of polynomials $\tilde \Pi_v(.,\tilde\boldsymbol\xi)$ with $\tilde\boldsymbol\xi$ fixed and the
parameter $v$ varying down or up to $0$, noting that during that variation the following
Lemma 1 keeps the set of zeroes within a fixed ball of finite radius $R(|\tilde v|,\{\tilde\boldsymbol\xi\})$ -- so none
can escape or come in across the boundary of that ball --, and recalling that all roots of
$\tilde\Pi_0(.,\tilde\boldsymbol\xi)$ have strictly negative real part, we see that the
Intermediate Value Theorem implies that there exists a value $v^*\in (-1,1)$ for which
$\tilde\Pi_{v^*}(.,\tilde\boldsymbol\xi)$ has a root with vanishing real part. This would contradict
Proposition 1.  \hfill $\Box$

\begin{lemma}
For all $\bar v\in [0,1)$ and any compact set $X\subset \mathbb R^3$, there exists $R=R(X,\bar v)>0$ such
that if $\tilde\lambda\in\mathbb C$ satisfies
$
{\tilde\Pi}_v(\tilde\lambda,\tilde\boldsymbol\xi)=0$ with $|v|\le \bar v\text{ and }\tilde\boldsymbol\xi\in X,
$
then
$|\tilde\lambda|\le R.$ In particular, the set of zeros of $\tilde\Pi_v(.,\boldsymbol\xi)$ depends continuously
on the parameter $v$.
\end{lemma}
\begin{proof}
Fix $\bar v$ and $X$ as stated, and consider the family
$$
\tilde\Pi_v(.,\tilde\boldsymbol\xi),\quad
(v,\tilde\boldsymbol\xi)\in\bar X=[-\bar v, \bar v]\times X
$$
of polynomials in $\tilde\lambda$. The assertion is true if (a) all coefficients of these polynomials
satisfy a uniform bound from above and (b) the leading-order coefficient is bounded uniformly away
from $0$. But (a) is clear due to the compactness of $\bar X$; notably, the Lorentz factor $\gamma$ is
bounded by its value $(1-\bar v)^{-1/2}$ for $v=\bar v$. Property (b) comes from the
(though sharp, cf.\ \cite{FT1})
 causality at
$|\boldsymbol\xi|=\infty$: for the
polynomial $\Pi$ in $\lambda$ and $\xi_1$ that gives the
dispersion relation with respect to the rest frame (considering $\xi_2,\xi_3$ as
parameters), the coefficient of the highest power $\lambda^n$ of $\lambda$
has absolute value $a>0$ larger or equal to that of the highest power $\xi_1^n$ of $\xi_1$.
In the mixing of $\tilde\lambda$ and $\tilde\xi_1$ through the Lorentz transformation
\eqref{Lorentz_var}, the coefficient of $\tilde\lambda^n$ in $\tilde\Pi_v$ thus remains bounded
from below by $(1-\bar v^n)a>0$.
\end{proof}

\section{Subluminality in absence of heat conduction}

In this and the next section we show that our system does not allow for
luminal modes. As we know from \cite{FT3} that for $v=0$ all modes are subluminal,
the continuity of the solution set of \eqref{disp_rel} with respect to the parameter $v$
according to Lemma 1 implies that all modes are subluminal, for any $v\in (-1,1)$  -- the assertion of Theorem 2.

Details of our proof  depend on whether the fluid does allow for heat conduction
or not. In this section, we consider the easier case that it does not, i.e., we assume that
the heat conductivity vanishes, $\chi=0$. We prove the non-existence
of luminal modes in the longitudinal case and in the transverse case separately in two subsections.

\subsection{Longitudinal modes}   \label{Sec_long_modes_chi0}
\setcounter{equation}0
The purpose of this subsection is to prove that there are no longitudinal luminal modes:
\begin{prop}\label{Prop_nonlum_L}
For any $v\in(-1,1)$,
$\tilde\Pi^L_v(\alpha\pm i|\tilde\boldsymbol\xi|,\tilde\boldsymbol\xi)\neq 0$ for all
$\alpha<0$ and $\tilde\boldsymbol\xi\in\mathbb R^3\setminus\{0\}$.
\end{prop}

To that end we fix
\begin{eqnarray}\label{fixcv}
c_s^2\in (0,1)\quad\text{and}\quad v\in(-1,1)
\end{eqnarray}
and lead the assumption that for some $\tilde\boldsymbol\xi\in\mathbb R^3\setminus\{0\}$
and some $\tilde\alpha<0$
\begin{eqnarray}\label{widerspruchsannahme}
\tilde\Pi^L_v(\tilde\alpha+iw,\tilde\boldsymbol\xi)=0\quad\text{with }w=|\tilde\boldsymbol\xi|\text{ or }w=-|\tilde\boldsymbol\xi|
\end{eqnarray}
to a contradiction. The idea is to show that on the one hand the imaginary part of equation \eqref{widerspruchsannahme} forces $\tilde \alpha$ to lie in a certain interval $\mathcal A\subset\mathbb R$ on which, on the other hand, the real part of \eqref{widerspruchsannahme} cannot hold.

Writing
$$
a=\gamma^2(1-c_s^2),\quad b= \gamma^2(1-c_s^2v^2),
$$
and using $\gamma^2 -1 = \gamma^2 v^2$ and $\gamma^2 (c_s^2-v^2) - c_s^2 = -a v^2$,
we find from \eqref{Pi_L_gen2} that
\begin{align}
{\rm Im}(\tilde\Pi^L_v(\alpha+iw,\tilde\boldsymbol\xi))
&=   \alpha  \big(2a +  \gamma \alpha \big) v\tilde\xi_1 + \alpha w \big( 2b + 3\gamma\alpha \big) \label{IM_PiL} \\
{\rm Re}(\tilde\Pi^L_v(\alpha+iw,\tilde\boldsymbol\xi)) &=
-av^2 \tilde\xi_1^2
- 2 w \big( a + \gamma \alpha \big) v\tilde\xi_1 + b (\alpha^2- w^2)
+ \gamma \alpha (\alpha^2 - 2 w^2)  + c_s^2 w^2.      \label{RE_PiL}
\end{align}
Introducing
\begin{equation}  \label{0chi_p1&2}
p_1(\alpha)= 2 a + \gamma\alpha
\hspace{1cm} \text{and} \hspace{1cm}
p_2(\alpha)= 2 b + 3\gamma\alpha,
\end{equation}
we write the imaginary part of \eqref{widerspruchsannahme} as
\begin{eqnarray}  \label{0chi_xi_soln}
p_1(\tilde\alpha)v\tilde\xi_1 = - p_2(\tilde\alpha) w.
\end{eqnarray}
\begin{lemma} \label{0chi_adm-reg_lemma}
(a) For any solution of \eqref{widerspruchsannahme}, $\tilde\alpha$ lies in the
interval
$$
\mathcal A=[\min\{\alpha_-,\alpha_+\},\max\{\alpha_-,\alpha_+\}]
$$
with
\begin{eqnarray} \label{0chi_alpha_+-}
\alpha_\pm
= \ -2\gamma \frac{1\pm v}{3\pm v}(1\mp vc_s^2).
\end{eqnarray}
(b) $\mathcal{A}$ consists of a single point if and only if
either $v=0$ or
\begin{eqnarray} \label{0chi_cr}
c_s^2 = 2/(3-v^2),
\end{eqnarray}
and in this case
\begin{eqnarray}\label{3alphas}
\alpha_- = \alpha_+ = \frac{-2}{\gamma(3-v^2)} \equiv \alpha_*.
\end{eqnarray}
\end{lemma}

\begin{proof}
To prove part (a), observe that \eqref{0chi_xi_soln} and the inequality $0<|\tilde\xi_1| \leq |w|$ implies
\begin{eqnarray}  \label{0chi_admis_reg}
|p_2(\tilde\alpha)|\le |vp_1(\tilde\alpha)|.
\end{eqnarray}
Inequality \eqref{0chi_admis_reg} holds if and only if $\tilde\alpha$ lies between the zeros of the affine functions
$$
p_+ \equiv p_2 + v p_1\quad \text{and}\quad p_- \equiv p_2 - v p_1
$$
given by
\begin{equation}  \label{0chi_i+-}
p_\pm(\alpha)   =  (3 \pm v) \gamma\alpha + 2 (b\pm va ),
\end{equation}
and these zeros are $\alpha_+$ and $\alpha_-$.
Part (b) directly follows by equating $\alpha_+$ and $\alpha_-$, which implies that either \eqref{0chi_cr} holds or $v=0$.
\end{proof}

Next, multiplying the real part of the dispersion relation \eqref{widerspruchsannahme} by $p_1^2(\alpha)$
and using \eqref{0chi_xi_soln} to substitute for $v \tilde\xi_1 p_1(\tilde\alpha)$ gives
\begin{equation}   \label{0chi_A-B-eqn}
0= A(\tilde\alpha) w^2 + B(\tilde\alpha)p_1^2(\tilde\alpha),
\end{equation}
with
\begin{align} \label{0chi_A&B}
A(\alpha) &\equiv   -a p_2(\alpha)^2 + 2 \big( a + \gamma \alpha \big) p_1(\alpha)p_2(\alpha)
-  (2 \gamma \alpha  + a) p_1(\alpha)^2 ,  \cr
&= 4\gamma^3\alpha^3 + 4b \gamma^2\alpha^2 - 4a c_s^4   \cr
&=  4 \gamma^2 B(\alpha)   - 4a c_s^4    \\
B(\alpha) &\equiv    \big(  b  + \gamma\alpha \big) \alpha^2 ,
\end{align}
as an equation equivalent to the dispersion relation \eqref{widerspruchsannahme}.
The essential idea is now to control the signs of $A$ and $B$ in a way that shows that
\eqref{0chi_A-B-eqn} cannot hold.
This is the subject of the next two lemmas.

\begin{lemma} \label{0chi_lemma_A}
(i) The polynomial $A$ is positive on the interior of $\mathcal{A}$.\\
(ii) $A$ vanishes on $\partial\mathcal A$ if and only if $c_s^2=2/(3-v^2)$, (in which case $\mathcal{A}$ is a singleton).
\end{lemma}

\begin{proof}
We write $\tilde A(\alpha) \equiv  4v^2  A(\alpha)$ as the quadratic form
\begin{eqnarray} \label{0chi_tA}
\tilde A =   r_1 {p_+}^2     +  2 r_2 p_+ p_-    + r_3  {p_-}^2   ,
\end{eqnarray}
in $p_+,p_-$ with coefficients
\begin{eqnarray} \nonumber
r_1(\alpha) &=&   - (1-v) \big((1-v) a + 2\gamma \alpha\big)      \cr
r_2(\alpha)&=&      (1-v^2) a  +2\gamma \alpha  \cr
r_3(\alpha) &=&    - (1+v) \big((1+v) a + 2\gamma \alpha\big)   ,
\end{eqnarray}
by substituting $p_2= \frac12(p_+ + p_-)$ and $vp_1=\frac12(p_+ - p_-)$  into \eqref{0chi_A&B}, and using the identities $v^2 \gamma^2 = \gamma^2 -1$, $\gamma^{-2}=1-v^2$.

We now investigate whether $\tilde A$ vanishes at a boundary point of $\mathcal A$. It clearly does vanish wherever both $p_+$ and $p_-$ vanish simultaneously; this is only the case when $c_s^2=2/(3-v^2)$, at the point displayed in equation \eqref{3alphas}. Assume now that
\begin{eqnarray}
c_s^2\neq \frac{2}{3-v^2} \quad \text{and} \quad  v \neq 0.
\end{eqnarray}
We first note that neither $p_+$ and $r_3$ nor $p_-$ and $r_1$ have common zeros. For assuming, say,
that $p_+$ and $r_3$ did, i.e.,
\begin{equation}
\alpha_+ \equiv  -2\gamma \frac{1+v}{3+v}(1-vc_s^2) \ = \
\alpha_3 \equiv  \frac{c_s^2 - 1}{2 \gamma (1-v)}     ,
\end{equation}
or equivalently
\begin{eqnarray}
(3+v)(c_s^2 - 1)   =   -4 \gamma^2 (1- v^2)(1- v c_s^2) ,
\end{eqnarray}
yields the contradiction that $c_s^2 = -\frac13$.
It now follows from \eqref{0chi_tA} that $\tilde A$ does not vanish at the boundary points of $\mathcal
A$ unless they coincide. In the special case $v=0$,
when by \eqref{3alphas} $\mathcal{A}$ collapses to the single point $\alpha_\pm = - 2/3$,
we find that $A(-2/3)= 16/27 - 4ac_s^4$, which vanishes if and only if $c_s^2 = 2/3$,
i.e., when \eqref{0chi_cr} holds. This proves assertion (ii).

To show (i), note that the third order polynomial $A$  satisfies
\begin{eqnarray}
A(-c_s^2/\gamma)=0,\quad A(0)<0,\quad\text{and} \quad A(+\infty)>0
\end{eqnarray}
and thus has both a negative and a positive root, and all roots of $A$ are real. By this, since $A$ has no zeros on $\partial A$ for
\begin{eqnarray}\label{notsep}
c_s^2\neq \frac{2}{3-v^2},
\end{eqnarray}
it follows that $A$ is positive on $\mathcal A$ for \emph{all}\ values of the parameters
\eqref{fixcv} as soon as it is positive on  $\mathcal A$ for one pair $(c_s^2,v)$ with
$c_s^2<2/(3-v^2)$ as well as for one pair with $c_s^2>2/(3-v^2)$. Indeed, the limiting cases
$(0,0)$ and $(1,0)$ have this property.   To see this, recall that  for $v=0$,  $\mathcal A$ is the singleton $\{-2/3\}$ and
\begin{equation}
A(\alpha) =  4 \big( \alpha^3 +  \alpha^2 - (1-c_s^2) c_s^4 \big),
\end{equation}
and thus on $\mathcal A$, $A$ assumes the value
\begin{equation}
{A}(-2/3)  =   \frac{16}{27} - 4(1-c_s^2) c_s^4  =  \frac{16}{27}>0
\end{equation}
both for $c_s^2=0$ and $c_s^2=1$. This completes the proof.
\end{proof}

\begin{lemma}  \label{0chi_lemma_B}
The polynomial $B$ is positive on $\mathcal{A}$.
\end{lemma}
\begin{proof}
By \eqref{0chi_A&B}, $B$ differs from $A$ by addition with a positive constant,  $4\gamma^2 B(\alpha)   = A(\alpha) + 4a c_s^4$. Thus, since by Lemma \ref{0chi_lemma_A} $A$ is positive on $\mathcal{A}$, $B$ is positive on $\mathcal{A}$ as well.
\end{proof}

In view of equation \eqref{0chi_A-B-eqn}, Lemmas \ref{0chi_lemma_A} and \ref{0chi_lemma_B}
readily imply the assertion of Proposition 2 except in the case \eqref{0chi_cr}.
In this case, in which $\alpha_+ = \alpha_- = \alpha_*$ by Lemma \ref{0chi_adm-reg_lemma},
we find from \eqref{0chi_p1&2}
that $\alpha_*$ is also a common zero of $p_1$ and $p_2$.
This implies by \eqref{0chi_xi_soln} that the imaginary part of the dispersion relation
\eqref{widerspruchsannahme} vanishes for ${\tilde\alpha} = \alpha_*$ and we cannot exactly
argue as above.

\begin{lemma}\label{sepcaseforchi=0}
In case $c_s^2=2/(3-v^2)$, one finds ${\rm Re}(\tilde\Pi_v^L(\alpha_*+iw,\tilde\boldsymbol\xi))\neq 0$.
\end{lemma}

\begin{proof}
In this case, $p_1,p_2$ and $A$ have the common zero $\alpha_*$ and $\hat A\equiv A/p_1^2$
satisfies $\hat A(\alpha_*)>0$. Instead of \eqref{0chi_A-B-eqn} one uses
$$
0=\hat A(\alpha_*)w^2+B(\alpha_*)
$$
which is a contradiction by Lemma \ref{0chi_lemma_B}.
\end{proof}

\subsection{Transverse modes}
The purpose of this subsection is to prove that there are no transverse luminal modes:
\begin{prop}\label{Prop_nonlum_T}
$\tilde\Pi^T_v(\alpha\pm i|\tilde\boldsymbol\xi|,\tilde\boldsymbol\xi)\neq 0$ for all
$\alpha<0$ and $\tilde\boldsymbol\xi\in\mathbb R^3\setminus\{0\}$.
\end{prop}

\begin{proof}
We fix $c_s^2\in (0,1)\quad\text{and}\quad v\in(-1,1)$  and lead the assumption that for some $\tilde\boldsymbol\xi\in\mathbb R^3\setminus\{0\}$ and some $\tilde\alpha<0$
\begin{eqnarray}\label{wid_ann_T}
\tilde\Pi_v^T(\tilde\alpha+iw,\tilde\boldsymbol\xi)=0\quad\text{with} \quad w=|\tilde\boldsymbol\xi|\text{ or }w=-|\tilde\boldsymbol\xi|
\end{eqnarray}
to a contradiction, following the idea of proof for Proposition \ref{Prop_nonlum_L} above. That is, in a first step we show the imaginary part of equation \eqref{wid_ann_T} forces $\tilde \alpha$ to lie again in some range $\mathcal A\subset\mathbb R$, and in the second step we show that the real part of \eqref{wid_ann_T} cannot hold on $\mathcal A$.  To begin, we substitute the ansatz $\tilde\lambda = {\tilde\alpha} + i w$ into the expression for $\tilde\Pi^T$ in \eqref{Pi_T_gen2}, which gives
\begin{align} \nonumber
\tilde\Pi_v^T(\tilde\lambda,\tilde\boldsymbol\xi)

&= \gamma({\tilde\alpha}+i(w+v\tilde\xi_1))+\hat\sigma\gamma^2\big({\tilde\alpha}+i(w+v\tilde\xi_1)\big)^2    +   \gamma^2 \big(-i v{\tilde\alpha}+(vw+\tilde\xi_1)\big)^2+ w^2 - \tilde\xi_1^2.
\end{align}
From this we find that the real and imaginary parts of condition \eqref{wid_ann_T} are given by
\begin{align}
\label{ImPiT}
0&=\text{Im\ \!}\tilde\Pi_v^T (\tilde\alpha+iw,\tilde\boldsymbol\xi)=\gamma(w+v\tilde\xi_1)
+2\hat\sigma\gamma^2{\tilde\alpha}(w+v\tilde\xi_1)
-2\gamma^2 v{\tilde\alpha}(vw+\tilde\xi_1),\\
\label{RePiT}
0&=\text{Re\ \!}\tilde\Pi_v^T (\tilde\alpha+iw,\tilde\boldsymbol\xi)
 = \big( 1 + (\hat\sigma-v^2)\gamma{\tilde\alpha} \big) \gamma {\tilde\alpha} + \gamma^2 (1-\hat\sigma)\big(w^2 + 2 wv \tilde\xi_1 + v^2 \tilde\xi_1^2 \big).
\end{align}
For the first step, we write \eqref{ImPiT} as
\begin{equation}\label{k}
 0=q_1({\tilde\alpha}) v \tilde\xi_1+q_2({\tilde\alpha})w
\end{equation}
for
\begin{eqnarray}
q_1(\alpha) \equiv  1+2\gamma\alpha(\hat\sigma-1)
\hspace{1cm} \text{and} \hspace{1cm}
q_2(\alpha) \equiv   1+2\gamma\alpha(\hat\sigma-v^2).
\end{eqnarray}
Since $|\tilde\xi_1| \leq |w|$, equation \eqref{k} implies that
\begin{eqnarray}
{\tilde\alpha} \in   \mathcal A \equiv \{\alpha:\: |q_2(\alpha)|\le |vq_1(\alpha)|\}
\end{eqnarray}
with $\mathcal A$ the interval between the
zeros of
$$
q_\pm=q_2\pm vq_1= \big(2\gamma (\hat\sigma\mp v)\alpha +1 \big)(1\pm v),
$$
i.e.,
$$
\mathcal A=[\min\{\alpha_-\alpha_+\},\max\{\alpha_-\alpha_+\}]\quad\text{with }
\alpha_\pm = -\frac{1}{2\gamma (\hat\sigma \mp v)} <0.
$$

For the second step, we multiply \eqref{RePiT} by
$q_1^2$ and, after substitution of \eqref{k}, write the resulting expression as
\begin{equation}\label{Aw^2+B=0}
0=A({\tilde\alpha})w^2+B({\tilde\alpha})
\end{equation}
with
\begin{eqnarray} \nonumber
A(\alpha) &=& -4 (\hat\sigma-1) \alpha^2  \quad   \leq 0, \quad \text{as } \hat\sigma \geq 1,     \cr
B(\alpha)  &=& q_1(\alpha)^2 \big(1 + (\hat\sigma -v^2)\gamma \alpha \big)\gamma \alpha .
\end{eqnarray}
As the zeros $0$, $\alpha_0$ and $\alpha_1$ of $B$ lie outside of $\mathcal{A}$,
$$
\alpha_0 = -\frac{1}{\gamma (\hat\sigma-v^2)} < \alpha_\pm < \alpha_1 = -\frac{1}{2\gamma(\hat\sigma-1)} <0,
$$
equation \eqref{Aw^2+B=0} is violated on $\mathcal A$, implying the sought after contradiction.
\end{proof}

\section{ Subluminality in presence of heat conduction}
\setcounter{equation}0
We turn to fluids allowing for heat conduction, $\chi>0$.
As heat conduction does not influence linear transverse
waves -- formally: $\chi$ does not appear in $\Pi_v^T$ --, it is obvious that  Proposition
\ref{Prop_nonlum_T} keeps valid.
The proof of Theorem \ref{ThmC} will thus be completed once we have shown
\begin{prop} \label{nonlum_L_chi>0}
The assertion of Proposition  \ref{Prop_nonlum_L} remains true for $\chi>0$.
\end{prop}
The rest of the section is devoted to proving this. In principle we proceed as in section 3.1.
We begin by fixing the parameters
$$
c_s^2 \in (0,1),  \hspace{1cm} \hat\chi = c_s^2\chi/\sigma >0  \hspace{1cm} \text{and} \hspace{1cm} v\in (-1,1),
$$
assume for contradiction there exist ${\tilde\alpha}< 0$ and $\tilde\boldsymbol\xi \in \mathbb R^3\setminus \{0\}$
such that \eqref{widerspruchsannahme} holds, and show analogues of Lemmas \ref{0chi_adm-reg_lemma}
to \ref{sepcaseforchi=0}.
This time,
\begin{align}
{\rm Im}(\tilde\Pi_v^L(\alpha+iw,\boldsymbol\xi))
&=  4 \hat\chi  w  \alpha^3  + \gamma \alpha^2 (\hat\chi +1) \big( 3 w+ v \tilde\xi_1  \big)  + 2 \gamma^2 \alpha (w+ v \tilde\xi_1)  -     2 c_s^2  \gamma^2  \alpha v (wv+ \tilde\xi_1)   \label{Im_Pi}   \\
{\rm Re}(\tilde\Pi_v^L(\alpha+iw,\boldsymbol\xi)))
&=   \hat\chi (\alpha^2  -4 w^2)\alpha^2
+ \gamma \alpha (\hat\chi +1) \big(\alpha^2 -2w(w+ v \tilde\xi_1)\big)
+ \gamma^2 \big(\alpha^2 - (w+ v \tilde\xi_1)^2\big)  \cr
& \ + c_s^2\big\{ r^2 - \gamma^2 \big(\alpha^2 v^2 - (wv+ \tilde\xi_1)^2\big) \big\}.    \label{Re_Pi}
\end{align}
Note first that equation \eqref{Im_Pi} is equivalent to
\begin{equation}  \label{xi_general_chi}
v p_1({\tilde\alpha}) \tilde\xi_1 = - w  p_2({\tilde\alpha})
\end{equation}
with now
\begin{equation} \label{p1_&_p2}
p_1(\alpha) = (\hat\chi+1)\gamma \alpha + 2a
\hspace{1cm} \text{and} \hspace{1cm}
p_2(\alpha) =  4 \hat\chi\alpha^2 + 3(\hat\chi+1)\gamma \alpha + 2b,
\end{equation}
where $a = \gamma^2(1-c_s^2)$ and $b  = \gamma^2(1-c_s^2v^2)$ as before. The constraint
$w^2\ge\tilde\xi_1^2$ again implies that ${\tilde\alpha}$ must lie in the set
\begin{eqnarray} \label{admis_reg}
\mathcal{A} \equiv \{\alpha:|p_2(\alpha)|\le|vp_1(\alpha)|\},
\end{eqnarray}
whose boundary points are the roots of the polynomials
$p_\pm \equiv p_2 \pm v p_1$ given by
\begin{equation} \label{p+_&_p-}
p_\pm(\alpha) =4 \hat\chi\alpha^2 + (3\pm v)(\hat\chi+1)\gamma \alpha + 2(b\pm va).
\end{equation}

\begin{lemma}  \label{Lemma_admissibility-region}
\textbf{(i)} The two polynomials $p_+$ and $p_-$ each have two distinct zeros,
$\alpha_\pm^- < \alpha_\pm^+$. They have a common zero if and only if either $v=0$,
in which case they
have both zeros in common,  $\alpha_-^+=\alpha_+^+$ and $\alpha_-^-=\alpha _+^-$ , or
\begin{eqnarray} \label{cr}
0 = \Gamma(c_s^2,\hat\chi,v)   \equiv   8 \hat\chi a^2 + \gamma^2 (b-3a)(\hat\chi+1)^2,
\end{eqnarray}
in which case $\alpha_-^+=\alpha_+^+$, while $\alpha_-^-\neq\alpha _+^-$.  \\
\textbf{(ii)} $\mathcal{A} = \mathcal{A}^- \cup \mathcal{A}^+$ with two disjoint compact intervals
$\mathcal A^- < \mathcal A^+$ having $\partial \mathcal{A^\pm} = \{\alpha_-^\pm,\alpha_+^\pm\}$. \\
\textbf{(iii)}  If $v=0$, the two intervals $\mathcal A^-$ and $\mathcal A^+$ are singletons, i.e., as
$v\rightarrow 0$ each interval $\mathcal A^\pm$ collapses continuously to a single point.
These points are the two distinct zeros of $p_+=p_-$. \\
\textbf{(iv)} For $v\neq 0$, $\mathcal A^-$ always has non-empty interior, while $\mathcal A^+$
collapses to a single point if and only if \eqref{cr} holds.
\end{lemma}

\begin{proof}
Part (i) directly follows from the discriminant $\Delta$ of $p_\pm$,
\begin{equation} \nonumber
\Delta= \gamma^2(3\pm v)^2(\hat\chi+1)^2-32\hat\chi(b\pm va),
\end{equation}
being strictly positive. This can be seen by using $(\hat\chi+1)^2/{\hat\chi}\geq 4$
to estimate
$$
\frac{b\pm va }{\gamma^2 (3\pm v)^2}<\frac18\le \frac{(\hat\chi+1)^2}{32\hat\chi}.
$$
Thus $p^+$ and $p^-$ both have two distinct zeros.
Clearly, if $v=0$, then $p_+=p_-=p_2$ and thus $\alpha_-^+=\alpha_+^+$ and $\alpha_-^-=\alpha _+^-$.
Assume now that $v\neq 0$ and $p^+$ and $p^-$, thus equivalently $p_1$ and $p_2$,
have a common zero. But the only zero of $p_1$ is
\begin{equation}    \label{zero_p1}
\alpha_* \equiv -\frac{2a}{\gamma(\hat\chi+1)} =\alpha^+_-=\alpha^+_+,
\quad \quad \text{and } \quad
p_2(\alpha_*) = 2\frac{(\hat\chi+1)^2}{\gamma^6} \Gamma(c_s^2,\hat\chi,v).
\end{equation}
Thus for $v\neq 0$, $\alpha_-^+=\alpha_+^+$ if and only if
\eqref{cr} holds while always $\alpha_-^-\neq\alpha_+^-$, and in any case
\begin{eqnarray}\label{4alphas}
\alpha_-^-,\alpha_+^-<\alpha_-^+,\alpha_+^+.
\end{eqnarray}

To see part (ii), note first that
as $|p_2(\alpha)| > |vp_1(\alpha)|$ for sufficiently large $|\alpha|$,
$\mathcal{A}$ is bounded.
The assertion thus follows from the fact $\partial \mathcal A=\{\alpha_-^-,\alpha_+^-,\alpha_-^+,\alpha_+^+\}$
and the above inequalities and equalities holding between these four numbers.

Parts (iii) and (iv) are direct consequences of (i) and (ii).
\end{proof}

Using that the imaginary part \eqref{Im_Pi} of the dispersion relation
vanishes for $\alpha=\tilde\alpha$,
we evaluate its real part \eqref{Re_Pi}
as
\begin{equation} \label{Aw^2+B}
A({\tilde\alpha}) w^2 + B({\tilde\alpha}) p_1^2({\tilde\alpha})   =0,
\end{equation}
in terms of the polynomials
\begin{align}
A(\alpha) &\equiv  \ -  a\: p_2(\alpha)^2  + 2 \big(\gamma \alpha (\hat\chi +1) +a \big) \: p_2(\alpha) p_1(\alpha)  - \big(4 \hat\chi \alpha^2  + 2 \gamma \alpha (\hat\chi +1) + a\big) p_1(\alpha)^2   \cr
& = \ 4\gamma^2 B(\alpha) - k    \label{A} \\
B(\alpha) &\equiv  \  \big(\hat\chi \alpha^2  + (\hat\chi +1) \gamma\alpha + b \big) \alpha^2 \label{B}
\end{align}
where $k \equiv 4 a c_s^4 \big( (\hat\chi-1)^2 + 4\hat\chi c_s^2 \big)^{-1}>0$.
\begin{lemma}  \label{Lemma_A}
The polynomial $A$ has the following properties:\\
(i) $A$ has a common zero with $p_+$ or $p_-$ if and only if the coincidence relation \eqref{cr} holds.\\
(ii) $A < 0$ on $\mathcal{A}^-$ and $A> 0$ on $\mathcal{A}^+$ in the generic case that \eqref{cr} does not
hold.\\
(iii)If \eqref{cr} holds, $A < 0$ on $\mathcal{A}^-$ and $A=0$ on $\mathcal{A}^+=\{\alpha_-^+=\alpha_+^+\}$.
\end{lemma}

\begin{proof}
To prove (i), we express $A$ as a quadratic form in $p^\pm \equiv p_2 \pm v p_1$.
Writing \eqref{A} as
$ \label{A_quad-form_0}
A  = a_1 p_1^2 + 2 a_2 p_1 p_2 + a_3 p_2^2
$
with
$
a_1 = - \big(4 \hat\chi \alpha^2  + 2 (\hat\chi +1) \gamma \alpha + a\big),
a_2 =  (\hat\chi +1) \gamma \alpha + a,
a_3 =   -  a,
$
and substituting now $vp_1 = (p_+ - p_-)/2$ and $p_2 =(p_+ + p_-)/2$ gives
\begin{equation}   \label{A_quad_form}
4v^2 A = A^+ (p^+)^2 + 2 C p^+ p^- + A^- (p^-)^2,
\end{equation}
with
\begin{eqnarray}
A^\pm =  a_1 \pm 2v a_2 +v^2 a_3
\hspace{.8cm} \text{and} \hspace{.8cm}
C = v^2 a_3 - a_1 .
\end{eqnarray}
This directly shows that $A$ vanishes either at a common zero of $p^+$ and $p^-$ --  whose existence is equivalent to \eqref{cr} by Lemma \ref{Lemma_admissibility-region} (i) --, or at common zeros of $A^\pm$ and $p^\mp$. However, neither $A^+$ and $p^-$, nor $A^-$ and $p^+$ have common zeros, since the resultants
$$
\begin{aligned}
\text{res}(A^\pm,p^\mp)&=\left|
\begin{matrix}
-4\hat\chi&2(\pm v-1)(\hat\chi+1)\gamma&-a(1\mp v)^2&0\\
0&-4\hat\chi&2(\pm v-1)(\hat\chi+1)\gamma&-a(1\mp v)^2\\
4\hat\chi&(3\mp v)(\hat\chi+1)\gamma&2(b\mp va)&0\\
0&4\hat\chi&(3\mp v)(\hat\chi+1)\gamma&2(b\mp va)
\end{matrix}
\right|\\
&=-4\hat\chi\left((1+3c_s^2)\hat\chi^2-2(2c_s^4+c_s^2+1)\hat\chi+(1+3c_s^2)\right)
\end{aligned}
$$
do not vanish for $\hat\chi >0$. This proves (i).

We now prove part (ii). For this we first verify (ii) in the special case
that $\hat\chi=1,v=0,c_s^2\neq 1/2$. We then extend this to the general parameter case by a continuity argument. So assume
$\hat\chi=1,v=0,c_s^2\neq \frac12$.  In this case,
$$
A(\alpha)= 4  B(\alpha)-k\quad\text{with }
B(\alpha)=(\alpha^2+2\alpha+1)\alpha^2
\text{ and }
k=c_s^2(1-c_s^2),$$
and
$$
p^\pm(\alpha)   = 4\alpha^2 + 6 \alpha + 2=2(\alpha+1)(2\alpha+1).
$$
At the zeros of $p^\pm$, we find that
$$
A(-1)=-k< 0
\hspace{1cm} \text{and}  \hspace{1cm}
A(-1/2)=1/4-c_s^2(1-c_s^2)>0.

$$
Since for $v=0$ the intervals $\mathcal{A}^-$ and  $\mathcal{A}^+$ collapse to the two zeros of
$p^+\equiv p^-$, we showed that $A$ is strictly negative on $\mathcal{A}^-=\{-1\}$ and strictly
positive on $\mathcal{A}^+=\{-\frac12\}$. This proves (ii) in the special case $\hat\chi=1,v=0$
for $c_s^2\neq 1/2$.

For $(\hat\chi,v)=(1,0)$ and $c_s^2\neq1/2$, the shape of $A$  is thus as follows:
As $A$ is of degree 4 with $A(\pm\infty)>0$, negative on $\mathcal A^-$, positive on $\mathcal A^+$
and has $A(0)=-k<0$, it has four real zeros, namely
one zero to the left of $\mathcal{A}^-$, a second one between $\mathcal{A}^-$ and
$\mathcal{A}^+$, a third one between $\mathcal{A}^+$ and $0$, and a forth in $(0,\infty)$.
Now, according to (i), this whole situation, i.e., the signs of $A$ on $\mathcal A^\pm$ and
the number and positioning of real roots of $A$
relative to $\mathcal A^-,\mathcal A^+,\{0\}$
is robust against changes of the parameters $(c_s^2,\hat\chi,v)$ in
$X=(0,1)\times(0,\infty)\times(-1,1)$ as long as these avoid the
surface $\Gamma_0=\{\Gamma=0\}\subset X$. As $\Gamma_0$ intersects the line segment
$(0,1)\times\{(1,0)\}$ exactly at the point $(c_s^2,1,0)$ with $c_s^2=1/2$,
the situation, in particular the sign behavior of $A$ on $\mathcal A^\pm$,
thus continues to all parameter values $(c_s^2,\hat\chi,v)\in X\setminus \Gamma_0$,
in either of the two connected components $\Gamma_0$ partitions $X\setminus \Gamma_0$.

(iii) For parameter values with $\Gamma(c_s^2,\hat\chi,v)=0$, the interval $\mathcal{A}^+$ collapses to a
single point. By \eqref{A_quad_form}, $A$ vanishes at this point because it is a common zero of $p^-$
and $p^+$.
\end{proof}

\begin{lemma} \label{Lemma_B}
The polynomial $B$ has the following properties:\\
(i) $B$ has a common zero with $p^+$ or $p^-$ if and only if $\hat\chi =1$ and $v=0$.\\
(ii) $B \leq 0$ on $\mathcal{A}^-$ and $B >0$ on $\mathcal{A}^+$.
\end{lemma}

\begin{proof}
(i) By symmetry of $B$ and $p_\pm$ with respect of $\pm v$, it suffices to
consider only $p_+$,
\begin{equation}  \label{no_common_zeros_eqn1}
p_+(\alpha)  = 4 \hat\chi\alpha^2 + (3+ v)(\hat\chi+1) \gamma \alpha  + 2(b + va).
\end{equation}
As $B(\alpha)=\alpha^2 \tilde B(\alpha)$ with
\begin{eqnarray}\label{no_common_zeros_eqn2}
\tilde B(\alpha) = \hat\chi \alpha^2  + (\hat\chi +1) \gamma  \alpha   + b
\end{eqnarray}
and $p_+(0)\neq 0$, we are thus looking for common zeros
of $p_+$ and $\tilde B$.
Since
\begin{eqnarray} \label{no_common_zeros_eqn2a}
(p^+ - 4 \tilde B)(\alpha)
&=& (v-1)(\hat\chi+1) \gamma \alpha  - 2(b -va) ,\cr
(b p^+ - 2(b + va) \tilde B)(\alpha)
&=& 2(b -va)  \hat\chi\alpha^2  + \big(b + v(b - 2a) \big)(\hat\chi+1) \gamma \alpha,
\end{eqnarray}
the existence of a common zero would imply that
\begin{eqnarray} \nonumber
\left|\begin{array}{cc}
(v-1)(\hat\chi+1) \gamma & - 2(b -va) \cr  2(b -va)  \hat\chi  &   \big(b + v(b - 2a) \big)(\hat\chi+1) \gamma
\end{array}\right|=0
\end{eqnarray}
or, equivalently,
\begin{equation}  \label{no_common_zeros_eqn3}
 \frac{(\hat\chi+1)^2}{4\hat\chi}=  \frac{(1+c_s^2v)^2}{1+v(v+2)c_s^2}.
\end{equation}
Note now that the left hand side of \eqref{no_common_zeros_eqn3} is  $\ge 1$
with equality holding only for $\hat\chi=1$.
On the other hand, as
\begin{eqnarray} \nonumber
(1+ c_s^2v)^2 - \big(1+v(v+2)c_s^2\big) = v^2 c_s^2 (c_s^2 -1)\le 0\quad\text{with ``='' only for }v=0,
\end{eqnarray}
 the expression on the right hand side of \eqref{no_common_zeros_eqn3} is
$\le 1$ with equality holding only for $v=0$.

Regarding (ii), obviously $B>0$ on $\mathcal A^+$ since $4\gamma^2 B=A+k$ with $k>0$ by \eqref{A}.
We prove $B\le 0$ on $\mathcal A^-$ first in the special case $v=0$ and then treat the
general case $v\in (-1,1)$ by a continuity argument. Assume that $v=0$ and $\hat\chi\neq 1$
(as the case $v=0$ and $\hat\chi=1$ is already covered by assertion (i)).
By (iii) of Lemma \ref{Lemma_admissibility-region},
$\mathcal{A}^-$ consists of the single point $\alpha^-\equiv\alpha_-^-=\alpha_+^-$, and $\tilde B$, now given
by
\begin{equation} \nonumber
\tilde B (\alpha) = \hat\chi \alpha^2  +   (\hat\chi +1) \alpha   + 1,
\end{equation}
has the zeros
$-{1}/{\hat\chi}$ and $ -1$, at which
\begin{eqnarray} \label{B_sign_eqn1}
p_+\left(-\frac{1}{\hat\chi}\right) = \frac{1}{\hat\chi} -1
\hspace{1cm} \text{and} \hspace{1cm}
p_+(-1) =   \hat\chi -1.
\end{eqnarray}
In particular, $p_+$ has different signs at the two zeros of $\tilde B$. Correspondingly,
$\alpha^-$ lies between these zeros, and $\tilde B(\alpha^-)<0$.
This proves (ii) in the case $v=0$, and by continuity, this also holds true for small
variations of $v$ around $v=0$. By continuity, upon variation of $(v,\hat\chi)$ in the connected domain
$(-1,1)\times(0,\infty)\setminus \{(0,1)\}$, the only way for $B$ to become positive on $\mathcal A^-$
would be to share a zero with $p_-$ or $p_+$; this would contradict property (i).
\end{proof}

Lemmas \ref{Lemma_A} and \ref{Lemma_B} readily imply the assertion of Proposition
\ref{nonlum_L_chi>0} except in the case when \eqref{cr} holds; analogously to the
the situation for $\hat\chi=0$, this case is covered by the following observation,
which is proved in the same way as Lemma \ref{sepcaseforchi=0}.
\begin{lemma} \label{sepcaseforchi>0}
In case \eqref{cr} holds, one finds ${\rm Re}(\tilde\Pi_v^L(\alpha_*+iw,\tilde\boldsymbol\xi))\neq 0$.
\end{lemma}

\section{Mode stability for non-barotropic fluids}
\setcounter{equation}0
In this section, we consider modes in viscous, heat conductive flow
of general
non-barotropic fluids. A non-barotropic fluid is given by an equation of state
$p=p(\theta,\psi)$ that represents the
pressure as a function of temperature and thermodynamic potential; in terms of these the inviscid
parts of the energy-momentum tensor and the particle flux are
\begin{eqnarray}\label{TN}
T^{\alpha\beta}=\theta^3\frac{\partial p(\theta,\psi)}{\partial \theta}\psi^\alpha\psi^\beta
+p(\theta,\psi)g^{\alpha\beta},
\quad
N^\beta=\frac{\partial p(\theta,\psi)}{\partial \psi}\psi^\beta,
\end{eqnarray}
written in the Godunov variables $\psi^\alpha=u^\alpha/\theta$ and $\psi=(\rho+p)/(n\theta)-s$, where
$n,s$ denote the particle number density and specific entropy.
The Navier-Stokes-Fourier model proposed in \cite{FT2} is
\begin{eqnarray}\label{nsf}
\frac{\partial}{\partial x^{\beta}}\left(T^{\alpha\beta}+\Delta T^{\alpha\beta}_\Box\right)=0,
\quad
\frac{\partial}{\partial x^{\beta}}\left(N^\beta+\Delta N^\beta_\Box\right)=0
\end{eqnarray}
with $\Delta T^{\alpha\beta}_\Box$ from \eqref{DeltaT} and
\begin{eqnarray}\label{DeltaN}
-\Delta N^\beta
\equiv
\kappa g^{\beta\gamma}{\partial \psi\over \partial x^\gamma}
+\tilde\sigma\left(U^\beta \Pi^{\gamma\delta}-U^\delta\Pi^{\beta\gamma}\right)
{\partial U_\gamma\over \partial x^\delta},
\end{eqnarray}
where $\tilde\sigma=\sigma n/(\rho+p)$ and $\kappa$ is the coefficient of thermal conductivity.

Subsuming $\psi^\alpha$ and $\psi$ under $\psi^a, a=0,1,2,3,4,$ with $\psi^4=\psi$,
$T^{\alpha\beta}$ and $N^\beta$ as $T^{a\beta}$, and
$\Delta T^{\alpha\beta}_\Box$ and $\Delta N^\beta_\Box$ as $\Delta T^{a\beta}_\Box$,
the equations can be written concisely  as
\begin{eqnarray}\label{lin}
A^{a\beta g}\frac{\partial \psi_g}{\partial x^\beta}
=
B^{a\beta g\delta}_\Box\frac{\partial^2\psi_g}{\partial x^\beta\partial x^\delta};
\end{eqnarray}
with coefficients
$$
A^{a\beta g}=\frac{\partial T^{a\beta}}{\partial\psi_g},\quad a,g=0,1,2,3,4,
$$
representing the Euler part div$(T)$ and $B^{a\beta g\delta}_\Box$
the Navier-Stokes-Fourier part div$(-\Delta T_\Box)$.
With respect to the rest frame, the coefficients are given (cf.\ \cite{FT2,F20})
by the matrices
$$
A^{a0g}\equiv\begin{pmatrix}
\theta^3p_{\theta\theta}&0&\theta p_{\psi\theta}-p_\psi\\
0&\theta^2p_\theta I_3 &0\\
\theta p_{\psi\theta}-p_\psi&0&\theta^{-1}p_{\psi\psi}
\end{pmatrix},\quad
A^{ajg}\equiv
\begin{pmatrix}
0&\theta^2 p_\theta e_j^\top&0\\
\theta^2 p_\theta e_j&0&p_\psi e_j\\
0&p_\psi e_j^\top&0
\end{pmatrix},
$$
and
$$
B^{a0g0}_\Box=
\begin{pmatrix}
0&0&0\\
0&-\sigma I_3&0\\
0&0&-\kappa
\end{pmatrix},
\quad
B^{ajgk}_\Box=
\begin{pmatrix}
0&0&0\\
0&
\eta I_3
+(\frac13\eta+\check\zeta) e_j e_k^\top&0\\
0&0&\kappa
\end{pmatrix},
\quad B^{ajg0}_\Box=0_{5}.
$$
The purpose of this section is to prove
\begin{theo}
Assume the fluid has speed of sound $0<c_s<1$. Then for system \eqref{nsf} with
\eqref{DeltaT} and \eqref{DeltaN} and
\begin{eqnarray}\label{etahatzetakappa}
\eta>0,\quad \check\zeta\ge -\frac13\eta,\quad \sigma=\frac43\eta+\check\zeta,\quad\kappa\geq 0,
\end{eqnarray}
all non-trivial Fourier-Laplace modes with respect to any Lorentz frame are decaying.
\end{theo}

For the case that in addition to viscosity and heat conduction the dissipation also includes diffusion, an analogue of Theorem 3 is implicitly contained in Sec.\ 3 of  \cite{S2}. Note however that formally Theorem 3 does not readily follow from that analogue as the limit of vanishing diffusion is singular.

\begin{proof}
Using direct analogues of Proposition 1 and Lemma 1, we see that it is sufficient to
show the assertion for rest-frame modes. We can thus also assume w.\ l.\ o.\ g.\ that
the wave vector $\boldsymbol\xi$ is of the form $(\xi,0,0)$ with $\xi \ge 0$.

The dispersion relation $0=\Pi(\lambda,(\xi,0,0))$ factors again as $\Pi=\Pi^L(\Pi^T)^2$ with
longitudinal part
$$
\Pi^L(\lambda,\xi)=
\left|\begin{matrix}
\lambda a_1&i\xi a_3&\lambda a_2\\
i\xi a_3&\lambda a_3+\sigma(\lambda^2+\xi^2)&i\xi a_4\\
\lambda a_2&i\xi a_4&\lambda a_5+\kappa(\lambda^2+\xi^2)
\end{matrix}\right|
$$
and transverse part
$$
\Pi^T(\lambda,\xi)=a_3\lambda+\sigma\lambda^2+\eta\xi^2
$$
with
$
a_1=\theta^3p_{\theta\theta},\
a_2=\theta p_{\psi\theta}-p_\psi, \
a_3=\theta^2p_\theta, \
a_4=p_\psi, \
a_5=\theta^{-1}p_{\psi\psi}.
$
Regarding the transverse part, the assertion follows in exactly the same way as in the barotropic case in \cite{FT1,FT3}. To study the longitudinal part, we write
\begin{eqnarray} \nonumber
\Pi^L(\lambda,\xi) = \lambda ^3\: f_1 + \lambda^2 (\lambda^2 + \xi^2)\: f_2 + \lambda \xi^2\: f_3 + \xi^2 (\lambda^2 + \xi^2)\: f_4 + \lambda (\lambda^2 + \xi^2)^2 \: f_5
\end{eqnarray}
with
\begin{eqnarray}
f_1 &=& a_3 (a_1 a_5 - a_2^2) >0 ,    \cr
f_2 &=&  \sigma (a_1 a_5 - a_2^2) + \kappa a_3 a_1  >0  , \cr
f_3 &=&   a_3^2 a_5 - 2 a_3 a_2 a_4 + a_1 a_4^2 >0 ,  \cr
f_4 &=&   \kappa a_3^2 \geq 0 ,   \cr
f_5 &=&  \sigma \kappa a_1 \geq 0  ,
\end{eqnarray}
where
the signs of the coefficients are
determined as follows: The positive definiteness of $A^{a0g}$ implies positivity of the diagonal
entries $a_1, a_3, a_5 >0$ and of the determinant $a_1 a_5 - a_2^2>0$. This directly implies
$f_1>0$, $f_5\geq 0$ and $f_2 = (f_4 a_1 +\sigma f_1)/a_3>0$; clearly $f_4\geq 0$.
Lastly $f_3$ has the same sign as $f_1$ by the dispersion relation of the inviscid (Euler) part,
\begin{eqnarray} \label{dispersion_Euler}
\pi^L(\lambda,\xi)= \lambda ^3 f_1  + \lambda \xi^2 f_3 =0,
\end{eqnarray}
which tells us that ${f_3}/{f_1} = c_s^2$. Now, to prove decay, we assume for contradiction that
\begin{eqnarray} \label{contradiction_non-bar}
\Pi^L(\lambda,\xi)= 0 \hspace{1cm} \text{for} \hspace{.5cm} \lambda = i\beta,
\end{eqnarray}
and compute
\begin{eqnarray}
{\rm Re}(\Pi^L(i\beta,\xi))  &=&  (\xi^2 - \beta^2) \big( \xi^2 f_4 - \beta^2 f_2 \big),   \cr
{\rm Im}(\Pi^L(i\beta,\xi))  &=& \beta \Big( \xi^2 f_3 - \beta^2 f_1 + (\xi^2 - \beta^2)^2 f_5 \Big) .
\end{eqnarray}
Solving the real part of \eqref{contradiction_non-bar} hence implies either $\xi^2 = \beta^2$ or
$\beta^2 = \xi^2f_4/f_2$. In the first case, the imaginary part of
\eqref{contradiction_non-bar} would give $f_3=f_1$,
which contradicts $c_s^2<1$.  In the second case, the imaginary part of \eqref{contradiction_non-bar} would imply
\begin{eqnarray}
f_2 f_3 - f_4 f_1 + \xi^2 \big(1-  f_4/f_2 \big)^2 f_5 =0,
\end{eqnarray}
a contradiction since $f_2 f_3 - f_4 f_1 = ({1}/{a_3})\big(\sigma f_1 f_3 + f_4 (a_1 a_4 - a_2 a_3)^2\big) >0$.
We have thus proved that Re$(\lambda)\neq 0$ for all modes.

A trivial computation in the spirit of the proof of Lemma  6.6 in \cite{FT1} shows that Re$(\lambda)<0$
for small $\xi>0$. As Re$(\lambda)$ cannot reach $0$, the assertion follows by continuity.
\end{proof}