# Counterfactual Outcome Prediction using Structured State Space Model

Vishal Purohit

###### Abstract

Counterfactual outcome prediction in longitudinal data has recently gained
attention due to its potential applications in healthcare and social sciences.
In this paper, we explore the use of the state space model, a popular sequence
model, for this task. Specifically, we compare the performance of two models:
Treatment Effect Neural Controlled Differential Equation (TE-CDE) and
structured state space model (S4Model). While TE-CDE uses controlled
differential equations to address time-dependent confounding, it suffers from
optimization issues and slow training. In contrast, S4Model is more efficient
at modeling long-range dependencies and easier to train. We evaluate the
models on a simulated lung tumor growth dataset and find that S4Model
outperforms TE-CDE with 1.63x reduction in per epoch training time and 10x
better normalized mean squared error. Additionally, S4Model is more stable
during training and less sensitive to weight initialization than TE-CDE. Our
results suggest that the state space model may be a promising approach for
counterfactual outcome prediction in longitudinal data, with S4Model offering
a more efficient and effective alternative to TE-CDE.

conterfactual outcomes, state space model, continuous time modeling

##  1 Introduction

Healthcare systems are dynamic and continuously changing in order to improve
care and treatments for patients. In order to accomplish this goal, healthcare
providers rely on accurate and dependable tools that can predict patient
outcomes and provide guidance for treatment decisions. Standard clinical risk
prediction models are commonly used to forecast an individual’s risk of an
outcome based on their observed characteristics. However, these models are
frequently created using data from a population in which patients follow a mix
of treatment strategies. This makes them unsuitable for informing treatment
decisions. Counterfactual prediction models offer an alternative approach that
estimates a person’s risk of an outcome if they were to follow a specific
treatment pattern, taking into account other patient characteristics that
predict the outcome. Despite their potential advantages, the development of
decision models that predict counterfactual outcomes presents several
challenges. These challenges include determining the appropriate methods for
identifying and adjusting for confounding factors, accounting for time-varying
treatments and outcomes, and incorporating complex interactions between
different patient characteristics. Additionally, the high dimensionality of
healthcare data and the need to balance model complexity with interpretability
and generalizability further complicates the creation of accurate and reliable
counterfactual prediction models in healthcare. Nonetheless, the development
of such models has the potential to greatly enhance healthcare decision-making
and improve patient outcomes.

When dealing with longitudinal data, the observed data can occur at irregular
intervals, adding an additional layer of complexity that must be accounted for
by the algorithm. In addition to the challenge of irregular time-series data,
we also face the problem of time-dependent confounding. Confounding variables
are factors that are associated with both the exposure and the outcome of
interest, and they can introduce bias into causal inference if not properly
accounted for. In the context of cancer treatment, the history of patients’
covariates and their response to past treatments can be important confounding
variables that affect the choice of future treatments. Failing to account for
these variables can lead to biased estimates of the causal effect of a
particular treatment on patient outcomes. Addressing this issue requires the
use of various statistical methods, such as propensity score matching or
instrumental variables, to adjust for confounding variables and obtain more
accurate estimates of causal effects. However, identifying the appropriate
confounding variables to include in the analysis and the optimal method for
adjusting for them can be challenging. Moreover, the availability and quality
of data for these variables may vary across patients, making it difficult to
achieve a balanced and fair comparison between treatment groups. Thus, it is
crucial to carefully consider the potential sources of bias and confounding in
longitudinal data analysis and apply appropriate methods to minimize their
impact on causal inference. Additionally, developing novel and innovative
statistical methods that can handle irregular longitudinal data and time-
dependent confounding is an important area of ongoing research in the field of
causal inference.

Specifically, time-dependent confounding and distribution shifts pose unique
challenges in causal inference over time, which are not typically encountered
in standard time-series analysis. Conventional time-series models do not
account for the bias introduced by time-varying confounders, and therefore may
not be appropriate for analyzing observational data in which there is a
potential for confounding. To address this issue, researchers may need to use
more sophisticated causal inference methods that can account for time-varying
confounding and distribution shift, such as time-varying treatment models or
g-methods. These methods allow for the estimation of the causal effect of a
treatment or intervention over time, while also accounting for the potential
biases introduced by confounding factors that change over time.

Our study differs from TE-CDE (Seedat et al., 2022), which uses a continuous-
time model for estimating treatment effects over time. We use a general
sequence-to-sequence model that does not rely on continuity. Despite this
difference, our model significantly outperforms TE-CDE, achieving up to 100
times better results in some cases. For critical applications such as
counterfactual outcome estimation in a medical setting, it is essential to
have a model that can provide reliable and accurate estimates of treatment
effects. Our approach may not have the same theoretical soundness as TE-CDE,
but it has practical advantages in terms of computational efficiency and
scalability. However, we need further investigation to understand the factors
that affect the performance gap between our model and TE-CDE, and to explore
the strengths and weaknesses of each approach in different scenarios. We also
aim to develop new and innovative models that can combine the best features of
both approaches.

##  2 Related Work

This paper primarily deals with counterfactual outcome estimation in
irregularly sampled settings with time-dependent confounding. We briefly
outline the details of the method that is closest to our problem.

###  2.1 Continuous-Time Modeling of Counterfactual Outcomes Using Neural
Controlled Differential Equations

The ability to estimate counterfactual outcomes over time has immense
potential to revolutionize personalized healthcare by enabling decision-makers
to answer ”what-if” questions. However, current causal inference approaches
are limited in their ability to handle irregularly sampled data, which is a
common occurrence in practice. To address this issue, (Seedat et al., 2022)
propose a novel approach called Treatment Effect Neural Controlled
Differential Equation (TE-CDE), which models the underlying continuous-time
process explicitly using controlled differential equations. This allows for
the potential outcomes to be evaluated at any time point, even in the presence
of irregularly sampled data. TE-CDE incorporates adversarial training to
account for time-dependent confounding, which is a crucial challenge in
longitudinal settings not encountered in conventional time-series analysis.
The performance is evaluated on data generated by a lung tumour growth model
that reflects a range of clinical scenarios with irregular sampling.

###  2.2 SyncTwin: Treatment Effect Estimation with Longitudinal Outcomes
(Qian et al., 2021) and SurvITE: Learning Heterogeneous Treatment Effects from
Time-to-Event Data (Curth et al., 2021)

Though not related directly to our problem, we briefly explain these two
papers that deal with treatment effect estimation. SyncTwin is proposed for
estimating causal treatment effects using electronic health records. SyncTwin
learns a patient-specific time-constant representation from pre-treatment
observations to issue counterfactual predictions. SyncTwin is demonstrated to
be usable in real-world EHR and successfully reproduced the findings of a
randomized controlled clinical trial using observational data. On the other
hand SurvITE discussed the problem of inferring treatment effects from time-
to-event data, which involves estimating treatment effects on both
instantaneous risk and survival probabilities.

##  3 Problem Formulation

We follow the same problem formulation as in work (Seedat et al., 2022),. For
completeness, we summarize the problem formulation. Let $n$ be the number of
individuals whose covariates are collected over time [0,T]. Let
$\textbf{X}:[0,T]\rightarrow\mathbb{R}^{d}$. The covariate path over time
$[0,T]$ of each patient is $x_{i}\in\mathbb{R}^{d}$. Let $A$ define a
treatment process over period of time $[0,T]$ where $a_{i}\in\\{0,1\\}$.The
output $Y$ is dependent on treatment $A$ and patient covariates $X$. In our
paper, we examine the tumour size growth hence $y_{i}\in\mathbb{R}$ and
$\textbf{Y}:[0,T]\rightarrow\mathbb{R}$. We explicitly define time
$\mathcal{T}$ to be a record of times at which the observation from the
healthcare worker is made. For $i^{th}$ patient the record of time is denoted
by $t_{i}\in\mathcal{T}$. If S = {$X_{s},A_{s},Y_{s}$} represents a set of
observations over a period of time $[0,s]$. In longitudinal studies, the
observed data for each patient may be collected at irregular intervals due to
various reasons, such as missed appointments, treatment schedules, or patient
and healthcare worker preferences. This irregularity adds an additional layer
of complexity to the analysis, as the algorithm must be able to handle missing
data and account for the variability in observation times. In our study, we
aim to provide estimates of treatment effects at a specific time point beyond
the initial observation interval [0, s], denoted as $t_{i}^{{}^{\prime}}$. To
accomplish this goal, we make certain assumptions and use a specific
formulation that is described in more detail in (Seedat et al., 2022). By
addressing the challenges posed by irregularly sampled longitudinal data, we
can obtain more accurate and reliable estimates of treatment effects, which
can help inform clinical decision-making and improve patient outcomes.

##  4 S4Model : Structured State Space Model

The state space model (SSM) is a foundational scientific model mainly used in
fields such as control theory and computational neuroscience. A state space
model is defined as follows,

$\displaystyle x^{{}^{\prime}}(t)=\textbf{A}x(t)+\textbf{B}u(t),$ (1)
$\displaystyle y(t)=\textbf{C}x(t)+\textbf{D}u(t).$ (2)

In the field of control theory, the state space model (SSM) is a mathematical
framework used to describe the behavior of dynamic systems. It consists of a
set of differential equations that relate the system’s state variables to its
inputs and outputs. The state vector $x(t)$ is an N-dimensional vector that
represents the internal state of the system at time t, while the output vector
$y(t)$ is a vector that represents the system’s measurable outputs. The input
signal $u(t)$ is a one-dimensional vector that represents the inputs to the
system.

The state, input, output, and feedforward matrices A, B, C, D are used to
define the relationships between the state, input, and output vectors. In
recent work, SSM has been treated as a black box representation of the
sequence model where matrices A, B, C, D are learned using gradient descent.
To encode long sequences, HiPPO theory of continuous time memorization has
been leveraged, which specifies a class of certain matrices
$\textbf{A}\in\mathbb{R}^{N\times N}$ that allows the state $x(t)$ to memorize
the history of the input $u(t)$.

Linear State Space Models (LSSL) (Gu et al., 2021) have been studied
extensively in the literature and have been shown to have multiple
interpretations as an ordinary differential equation, a recurrent model, and a
convolution. This paper focuses on structured space models (Gu et al., 2022),
which are easier to train and offer superior performance in relation to other
sequence models. For more details, readers can refer to the rich literature on
state space models.

Under the discrete-time view of LSSL, given an input sequence
$(u_{1},u_{2},...)$ discretized by time step of $\Delta t$, the resulting
discretized version can be obtained with the help of bilinear transform that
converts matrix A to $\bar{\textbf{A}}$. The discrete SSM (assume
$\textbf{D}=0$ for simplicity) is given by,

$\displaystyle x_{k}=\bar{\textbf{A}}x_{k-1}+\bar{\textbf{B}}u_{k},$ (3)
$\displaystyle y_{k}={\textbf{C}}x_{k}.$ (4)

Now the SSM in Equation 3 and Equation 4 is a mapping between input
$u_{k}\rightarrow y_{k}$. The matrix
$\bar{\textbf{A}}=(\textbf{I}-\frac{\Delta}{2}.\textbf{A})^{-1}(\textbf{I}+\frac{\Delta}{2}.\textbf{A})$
and
$\bar{\textbf{B}}=(\textbf{I}-\frac{\Delta}{2}.\textbf{A})^{-1}.\Delta{\textbf{B}}$.
The fundamental bottleneck in computing the discrete-time SSM Equation 3 that
it involves repeated matrix multiplication by A. For example, computing
Equation 3 naively as in the LSSL involves L successive multiplications by A,
requiring $O(N^{2}L)$ operations and $O(NL)$ space. According to Lemma 3.2 in
(Gu et al., 2022), when A is structured then one can perform faster
computation. For completeness, the lemma is as follows,

######  Lemma 4.1.

(Diagonalization of A) Conjugation is an equivalence relation on SSMs
$(\textbf{A},\textbf{B},\textbf{C})\sim(\textbf{V}^{-1}\textbf{A}\textbf{V},\textbf{V}^{-1}\textbf{B},\textbf{C}\textbf{V})$
and resulting SSM is,
$\tilde{x}^{{}^{\prime}}=\textbf{V}^{-1}\textbf{A}\textbf{V}x_{k-1}+\textbf{V}^{-1}\textbf{B}u_{k}$
and $y_{k}=\textbf{C}\textbf{V}x_{k}$ with $x=\textbf{V}\tilde{x}$

The diagonalization of HiPPO matrix A is unstable due to large values in
matrix V rendering the diagonalization numerically infeasible. (Gu et al.,
2022) observed that HiPPO matrix can be decomposed as a sum of a normal and
low-rank matrix. With help of linear algebraic techniques (Gu et al., 2022)
showed the following,

######  Lemma 4.2.

(Normal Plus Low Rank represenation (NLPR) of A) All HiPPO matrices have an
NPLR representation
$\mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{*}-\mathbf{P}\mathbf{Q}^{T}=\mathbf{V}(\mathbf{\Lambda}-(\mathbf{V}^{*}\mathbf{P})(\mathbf{V}^{*}\mathbf{Q})^{*})\mathbf{V}^{*}$
for unitary $\mathbf{V}\in\mathbb{C}^{N\times N}$ , diagonal
$\mathbf{\Lambda}$, and low-rank factorization $\mathbf{P}$ ,
$\mathbf{Q}\in\mathbb{R}^{N\times r}$

Under the NLPR formulation over $\mathbb{R}$ the complexity of SSM reduces to
$O(N)$ operations where $N$ is the state size.

###  4.1 From S4 to Deep S4

In order to initialize an SSM with A HiPPO matrix, we need to specify the
matrix A with a specific initialization. This SSM is unitarily equivalent to
some $(\mathbf{\Lambda}-\mathbf{P}\mathbf{Q}^{*},\mathbf{B},\mathbf{C})$ for
some diagonal $\mathbf{\Lambda}$ and vectors
$\textbf{P},\textbf{Q},\textbf{B}$, $\textbf{C}\in\mathbb{C}^{N\times 1}$.
These matrices and vectors comprise S4’s 5N trainable parameters.

S4 is a sequence model that defines a map from
$\mathbb{R}^{L}\rightarrow\textbf{R}^{L}$, i.e., a 1-D sequence map. However,
deep neural networks (DNNs) typically operate on feature maps of size H
instead of 1. To handle multiple features, S4 defines c-independent copies of
itself and then mixes the H features using a position-wise linear layer. This
results in a total of $O(H^{2})+O(HN)$ parameters per layer. Nonlinear
activation functions are also inserted between these layers.

Overall, S4 defines a sequence-to-sequence map of shape (batch size, sequence
length, hidden dimension). It is worth noting that the core S4 module is a
linear transformation, but the addition of non-linear transformations through
the depth of the network makes the overall deep SSM non-linear. A diagram of
S4model is show in Figure 1 and respective hyper-parameters are shown in Table
2.

It is important to note that S4 is designed for structured space models, which
are easier to train and offer superior performance compared to other sequence
models. In previous work on Linear State Space Models (LSSL) (Gu et al.,
2021), they have been shown to have multiple interpretations, such as an
ordinary differential equation, a recurrent model, and as convolution. Readers
can refer to the rich literature on state space models (Gu et al., 2021, 2022)
for further details. Thus, S4 is a deep structured space model that leverages
the HiPPO theory of continuous time memorization to encode long sequences. By
defining H-independent copies of itself and using position-wise linear layers
to mix features, it can handle multiple features and achieve superior
performance compared to other sequence models.

###  4.2 Optimization Objective

In our problem formulation, we have have a deep S4 model learning from
$S=\\{X_{s},A_{s},Y_{s}\\}$ to predict $Y_{t^{{}^{\prime}}>s}$ and
$A_{t^{{}^{\prime}}>s}$. Let the $f_{\phi}:\mathbb{R}^{c\times
L}\rightarrow\mathbb{R}^{c\times L}$ represent the mapping function
represented by $c$ S4 models and $\phi$ represent all the learnable
parameters. let output $f_{\phi}$ be given by latent state $z$. This state is
$z$ and is mapped to the treatment outcome (a) and lung tumour size (y) using
two different linear transformations. The predicted outcome is over a period
of time $[t,t_{k}]$ with observation times $(t_{1},.....,t_{k})$. The mean
square error (MSE) of outcome prediction is defined as,

$\mathcal{L}^{(y)}=\frac{1}{k}\Sigma_{i=1}^{k}(y_{t_{i}}-\hat{y}_{t_{i}})^{2}.$
(5)

The treatment loss is calculated using cross-entropy loss, which is given by,

$\mathcal{L}^{(a)}=-\frac{1}{k}\Sigma_{i=1}^{k}\bigg{(}a_{t_{i}}\log(\hat{a}_{t_{i}})+(1-a_{t_{i}}\log(1-\hat{a}_{t_{i}})\bigg{)}$
(6)

We use a similar formulation as (Seedat et al., 2022), to balance the
representations and minimized the following loss function,

$\mathcal{L}_{total}=\frac{1}{n}\Sigma_{i=1}^{n}\bigg{(}\mathcal{L}_{i}^{(y)}-\mu\mathcal{L}_{i}^{(a)}\bigg{)},$
(7)

where $\mu$ is the hyper-parameter controlling the trade-off between treatment
and outcome prediction. Unlike TE-CDE, we fix the value of $\mu=0.5$ implying
outcome prediction to be prioritized over treatment predictions.

##  5 Experiments

In this section, we evaluate the ability of our proposed method to estimate
counterfactual outcomes using observational data. In real-world scenarios, it
is often impossible to obtain counterfactual outcomes, hence we use synthetic
data to perform an empirical evaluation of our approach. To simulate the
synthetic data, we utilize a simulation environment based on a
Pharmacokinetic-Pharmacodynamic (PK-PD) model of lung cancer tumor growth
(Seedat et al., 2022) This simulation framework allows us to compute
counterfactual outcomes for any time point and treatment plan. However, to
mimic the nature of real-world data, we simulate the irregularly sampled data
using the Hawkes process.

Our evaluation consists of two main parts. Firstly, we evaluate the accuracy
of our proposed method in estimating counterfactual outcomes for a given
treatment plan. Secondly, we compare the performance of our approach with
other state-of-the-art methods like TE-CDE for counterfactual estimation. We
estimate the counterfactual outcomes using our proposed approach and calculate
the mean squared error (MSE) between the estimated and simulated outcomes.

###  5.1 Datasets and Architecture

The PK-PD model used in our study has been previously described in (Seedat et
al., 2022) Our simulation environment allows for the manipulation of three key
parameters governing the generation of synthetic data. These parameters,
$\gamma_{c}$ and $\gamma_{r}$, are responsible for controlling the degree of
time-dependent confounding, whereby higher values of $\gamma_{c,r}$ indicate
that treatment assignment is influenced by tumour size. Cancer staging details
and the use of the Hawkes process are elaborated in (Seedat et al., 2022).
Within the Hawkes process, $\kappa$ is used to control the intensity of data
sampling between cancer stages. Specifically, higher stages of cancer result
in more frequent sampling, thus reflecting some of the characteristics of
real-world data that our algorithm may encounter. By manipulating these
parameters within our simulation framework, we are able to generate synthetic
data with known counterfactual outcomes, which serves as a robust testing
ground for the evaluation of our proposed method.

The deep S4 model consists of 4 layers with the latent dimension of 256. The
TE-CDE model uses a latent dimension of 128 for each of the two networks used
for treatment and outcome prediction. Both models use a batch size 32 for
training. In total, we have 10000 training datapoint, 1000 validation data
points (used for hyper-parameter tuning) and 10,000 test data points. We use
the validation data to find hyperparameters like learning rate. The best
learning rate for S4Model is 0.0005 for parameters A, B, C and 0.00002 for
other parameters. We use learning rate of 0.00001 for TE-CDE as it was the
best and most stable one. Both TE-CDE and S4Model are trained for 50 epochs on
Nvidia A30 GPU with 24GB VRAM.

###  5.2 Main Experiements

In Table 1 we show the side-by-side comparison of TE-CDE and S4Model
performance on the test dataset. Here $\kappa=\\{1,5,10,15,20\\}$, which shows
difficult levels of sampling intensity of the Hawks process. Higher value of
kappa implies sampling intensity is greater for more severe cancer stages.
Different values of $\gamma_{r,c}=\\{2,4,6,8,10\\}$ control the time-dependent
confounding. Under sever time-confounding (higher values of $\gamma_{r,c})$
our method still outperforms TE-CDE by a significant margin. We observe that
TE-CDE performs poorly in almost all combinations of $kappa$ and
$\gamma_{r,c}$. As we see in some of the case the TE-CDE model does not even
train completely due to early stopping criteria imposed by the authors.
Further, in Figure 2 we show the training loss curves for various combinations
of $\kappa$ and $\gamma_{r,c}$. Clearly, we see the disadvantage of using
neural-controlled differential equations as they are unstable to train and
diverge. Neural ODEs are well known for slower training as well as poor
convergence properties. The neural-controlled differential equations also
suffer from such issues and this is concerning since they are used in critical
applications like counterfactual outcome estimation. Additionally, TE-CDE
requires continuous values which are obtained through linear interpolation and
this need not reflect the real-world data. On the other hand, S4model
converges in every training setting with and outperforms TE-CDE by 100x at
most. Their superior convergence makes them a reliable candidate for modelling
longitudinal data. S4Model does not operate on continuous data hence it does
not require any interpolation. This clear evidence suggests that indeed the
general-purpose sequence models can beat the performance of continuous-time
models like TE-CDE.

Table 1: Results on test dataset. Here TC denotes Time Confounding parameters
$\gamma_{r,c}$ and SI denotes sampling intensity $(\kappa)$. Here total loss
in cross-entropy loss + RMSE

Algorithm | TE-CDE | S4Model
---|---|---
Metric | Total Loss | Total Loss
TC ($\rightarrow$) SI ($\downarrow$) | 2 | 4 | 6 | 8 | 10 | 2 | 4 | 6 | 8 | 10
$\kappa$ = 1 | 12.2971 | 3.5115 | 20.7077 | 31.8526 | 11.3300 | 0.3407 | 0.7396 | 1.2541 | 1.3585 | 1.3773
$\kappa$ = 5 | 3.4005 | 51.0468 | 3.4945 | 4.4499 | 38.0621 | 0.3358 | 0.7348 | 1.2421 | 1.3577 | 1.3694
$\kappa$ = 10 | 2.8759 | 4.0433 | 41.5988 | 10.3160 | 22.8966 | 0.3377 | 0.7428 | 1.2540 | 1.3609 | 1.3726
$\kappa$ = 15 | 59.9164 | 3.2268 | 13.2311 | 13.2911 | 10.2217 | 0.3403 | 0.7354 | 1.2465 | 1.3583 | 1.3742
$\kappa$ = 20 | 5.1836 | 3.1370 | 31.1797 | 8.3600 | 25.7174 | 0.3387 | 0.7396 | 1.2501 | 1.3530 | 1.3762

Figure 1: Diagram of Deep S4Model. Table 2: Hyper parameter table

S4model | TE-CDE
---|---
Hyper-parameter name | Value | Hyper-parameter name | Value
No. Layers | 4 | No. Layers | NA
Latent dimension | 128 | Latent dimension | 256
Batch size | 32 | Batch size | 32
Learning rate | 0.0005 | Learning rate | 0.00001
Optimizer | Adam | Optimizer | SGD
Number of epochs | 50 | Number of epochs | 50
Number of training samples | 10000 | Number of training samples | 10000
Number of test samples | 1000 | Number of test samples | 1000

Figure 2: Total loss (cross entropy + MSE loss) v/s epochs for various
$\kappa=\\{1,5,10,15,20\\}$. The values of $\gamma_{r,c}=\\{2,4,6,8,10\\}$.
The training curves and respective configuration of $\gamma_{r,c}$ are shown
as the last two numerical values in the plot legend. Figure 3: Bar plots
showing treated v/s untreated RMSE calculated on the test dataset for TE-CDE
and S4Model for various combination of $\kappa$ and $\gamma_{r,c}$

###  5.3 Impact of Domain Adversarial Training

As highlighted in (Seedat et al., 2022), domain adversarial training has been
shown to significantly improve the performance of TE-CDE. To investigate the
effectiveness of this technique on our proposed method, S4Model, we conducted
experiments to compare its performance with and without adversarial training.
To do so, we used the loss function defined in Equation 7 and set $\mu=0$ for
both models. We evaluated the performance of S4Model and TE-CDE side-by-side
on a dataset generated with $\kappa=5$ for various levels of time-confounding.
Our results show that the overall loss is significantly lower with domain
adversarial training, and S4Model outperforms TE-CDE by a considerable margin.
This highlights the effectiveness of domain adversarial training in achieving
balanced representations, which can effectively remove bias introduced by
time-dependent confounders and enable reliable counterfactual estimates. Our
findings suggest that domain adversarial training is a useful technique for
improving the performance of causal inference models, and should be considered
when analyzing complex datasets.

Table 3: Adv v/s Non-Adv

Algorithm ($\rightarrow$) | TE-CDE | S4Model
---|---|---
Metric ($\rightarrow$) | Total Loss | Total Loss
Time-Confounding ($\downarrow$) | Balanced | Imbalanced | Balanced | Imbalanced
$\gamma_{c,r}$ = 2 | 3.4005 | 18.0014 | 0.3358 | 1.8282
$\gamma_{c,r}$ = 4 | 51.0468 | 3.3101 | 0.7348 | 1.9711
$\gamma_{c,r}$ = 6 | 3.4945 | 5.8647 | 1.2421 | 2.9998
$\gamma_{c,r}$ = 8 | 4.4499 | 4.8950 | 1.3577 | 2.9572
$\gamma_{c,r}$ = 10 | 38.0621 | 11.2443 | 1.3694 | 2.7580

Table 4:  VRAM, time/epoch and number parameters of TE-CDE and S4Model

Algorithm | Training VRAM | Time/epoch(seconds) | # of parameters
---|---|---|---
TE-CDE | 1037MiB | 7.2 | 2238149
S4Model | 1475MiB | 4.4 | 679487

Table 5: Results of ablation study on number of layers in S4Model

Number of layers | Treatment Loss | Outcome Prediction Loss | Total Loss
---|---|---|---
2 | 0.0240 | 0.3203 | 0.3443
4 | 0.0267 | 0.3178 | 0.3445
6 | 0.0257 | 0.3042 | 0.3299
8 | 0.0240 | 0.2967 | 0.3207
10 | 0.0257 | 0.2956 | 0.3213

###  5.4 Training Efficiency and Stability

Apart from the superior empirical performance of S4Model we examine the
several practical metrics that are important for real-world use. To this end
we examine the amount of VRAM consumption, time taken per epoch, the total
number of trainable parameters and inference speed. It is critical that the
deployed model is smaller in memory footprint as well as faster on inferring
from observations without compromising on the overall performance. In Table 4
we show the necessary metrics. This shows that S4Model has 10x fewer
parameters and is 100x faster in training as well as inference yet achieves
10x and in some cases 100x better performance than TE-CDE. As noted previously
TE-CDE model can diverge during training. Even after multiple restarts, the
training did not converge due to extremely large loss values leading to ’nan’
errors. Note that we re-run TE-CDE multiple times yet it is very sensitive to
initialization conditions. Hence, making TE-CDE quite impractical for real-
world use.

##  6 Abalation Study

We study the effect of increasing the number of layers in S4Model and evaluate
it’s performance. Generally, increasing the number of layers should give us
good performance improvements. We show the results in Table 5 total loss
value, treatment loss and outcome prediction loss using various S4Models with
layers 2, 4, 6, 8, 10. We observe that increasing the number of layers
decreases the overall loss values and does not overfit the training data.
Another ablation study examines the latent dimension used by S4Model. We
examine the empirical performance for the latent dimension of sizes 8, 16, 32,
128, 256 and show the results in Table 6. We observe that S4Model is more
expressive when the latent dimension is higher and benefits the performance.

Table 6: Results on ablation study of the latent dimension of S4Model

Latent dimension | Treatment Loss | Outcome Prediction Loss | Total Loss
---|---|---|---
8 | 0.0528 | 0.7607 | 0.8135
16 | 0.0343 | 0.6842 | 0.7185
32 | 0.0298 | 0.5447 | 0.5745
128 | 0.0241 | 0.3748 | 0.3989
256 | 0.0242 | 0.3181 | 0.3423

##  7 Conclusion and Future Scope

In this paper, we examined the use of the general purpose sequence model -
Structured state space model for counterfactual outcome estimation. We
compared our method against the recently proposed method TE-CDE. Empirically
our method showed superior results compared to TE-CDE on the lung cancer
dataset. Under various time-confounding and sampling intensity settings we
outperform TE-CDE with a lower memory footprint, 10x less number of parameters
and guaranteed convergence. In summary, this work highlights the empirical
performance, scalability and efficiency of state space models for
counterfactual outcome prediction. As part of future work, we would like to
examine the applicability of the state space model under a sparse data
setting. We look for inspiration from compressive sensing theory (Wu et al.,
2019) and sparse coding to effectively reconstruct the sequence of observation
from a small subset of observations.

## References

  * Curth et al. (2021) Curth, A., Lee, C., and van der Schaar, M.  SurvITE: Learning heterogeneous treatment effects from time-to-event data.  In Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , 2021.
  * Gu et al. (2021) Gu, A., Johnson, I., Goel, K., Saab, K. K., Dao, T., Rudra, A., and Re, C.  Combining recurrent, convolutional, and continuous-time models with linear state space layers.  In Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , 2021.  URL https://openreview.net/forum?id=yWd42CWN3c.
  * Gu et al. (2022) Gu, A., Goel, K., and Re, C.  Efficiently modeling long sequences with structured state spaces.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=uYLFoz1vlAC.
  * Qian et al. (2021) Qian, Z., Zhang, Y., Bica, I., Wood, A., and van der Schaar, M.  Synctwin: Treatment effect estimation with longitudinal outcomes.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Seedat et al. (2022) Seedat, N., Imrie, F., Bellot, A., Qian, Z., and van der Schaar, M.  Continuous-time modeling of counterfactual outcomes using neural controlled differential equations.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162, pp. 19497–19521, 2022.
  * Wu et al. (2019) Wu, Y., Rosca, M., and Lillicrap, T. P.  Deep compressed sensing.  In _Proceedings of the 39th International Conference on Machine Learning_ , 2019.


# Dark Energy Radiation

Kim V. Berghaus  Department of Physics & Astronomy, The Johns Hopkins
University, Baltimore, MD 21218, USA  C.N. Yang Institute for Theoretical
Physics, Stony Brook University, NY 11794, USA  Peter W. Graham  Stanford
Institute for Theoretical Physics, Department of Physics, Stanford University,
Stanford, CA 94305, USA  David E. Kaplan  Department of Physics & Astronomy,
The Johns Hopkins University, Baltimore, MD 21218, USA  Guy D. Moore  Institut
für Kernphysik, Technische Universität Darmstadt, Schlossgartenstraße 2,
D-64289 Darmstadt, Germany  Surjeet Rajendran  Department of Physics &
Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA

###### Abstract

We show that if dark energy evolves in time, its dynamical component could be
dominated by a bath of dark radiation. Within current constraints this
radiation could have up to $\sim 10^{4}$ times more energy density than the
cosmic microwave background. We demonstrate particular models in which a
rolling scalar field generates different forms of dark radiation such as
hidden photons, milli-charged particles and even Standard Model neutrinos. We
find the leading effect on the late-time cosmological expansion history
depends on a single parameter beyond $\Lambda$CDM, namely the temperature of
the dark radiation today. Cosmological observations of this modified expansion
rate could provide a striking signature of this scenario. The dark radiation
itself could even be directly detectable in laboratory experiments, suggesting
a broader experimental program into the nature of dark energy.

###### Contents

  1. 1 Introduction
  2. 2 Friction and Radiation
    1. 2.1 Non-Abelian Gauge Sector
      1. 2.1.1 Non-Abelian with Fermions
    2. 2.2 Dark Photons and Milli-Charged Particles
    3. 2.3 Neutrinos
  3. 3 Cosmological Signatures
  4. 4 Discussion and Conclusion

##  1 Introduction

A variety of cosmological measurements have established that the universe is
currently undergoing accelerated expansion driven by dark energy. This dark
energy dominates the energy density of the universe. It does not clump into
galaxies and its evolution is constrained to be slower than that of the other
components of the universe. Notably, these properties have been established by
cosmological measurements that are sensitive only to the gravitational effects
of dark energy. Much like the case of dark matter, it is likely that the
properties of dark energy can be more deeply probed if the dark energy has
non-gravitational interactions with the standard model. How can the dark
energy couple to the standard model and what kinds of signals can these
interactions lead to? A systematic identification of these could lead to a
robust program of laboratory based probes of dark energy, similar to the
effort that is currently underway to identify the nature of dark matter.

Unlike dark matter, such a program has not been widely pursued for dark energy
due to the widespread theoretical assumption that the dark energy is a
cosmological constant. Indeed, a cosmological constant is the simplest
proposal for dark energy and such a constant cannot be probed in laboratory
experiments. But this is simply theory prejudice - the experimental nature of
physics requires observational probes to establish if the dark energy is
actually a constant or if it has a dynamical component, e.g., a kinetic
energy. This kinetic energy can couple to the standard model through non-
gravitational interactions and could potentially be probed in the laboratory.
In fact, there are significant theoretical motivations to search for such a
kinetic energy. The cosmological constant is observationally known to be
orders of magnitude smaller than identifiable standard model contributions to
it and thus appears to require a massive fine tuning. At present, the only
known mechanisms (e.g., [1, 2, 3, 4, 5, 6]) to naturally explain this fine
tuning requires dynamical evolution, implying the existence of a dynamical
component of dark energy.

Given the isotropic nature of dark energy, the simplest model of dynamical
dark energy is a light scalar field that is slowly rolling down its potential.
While the mass of this field has to be smaller than the Hubble scale ($\sim
10^{-43}$ GeV) today, its kinetic energy could be significantly larger —
constrained to be no larger than a few percent of the dark energy density
today $\sim\text{meV}^{4}$. The coupling of this kinetic energy to the
polarization of light (see e.g. [7, 8, 9] and references) and the spin of
nucleons/electrons has been investigated (see e.g. [10, 11, 12, 13]). In
these, the dark energy provides a homogeneous background and the motion of the
test particle against the background gives rise to polarization rotation/spin
precession. The spatial/temporal coherence and the frequency of these signals
are all set by the Hubble scale today i.e. the age of the universe.
Experimentally, while the extended temporal coherence is useful to combat
sensor noise, the DC nature of the signal gives rise to challenging
systematics.

In this paper, we point out that the phenomenology of dynamical dark energy is
richer than simply the motion of test particles against a homogeneous
background, giving rise to qualitatively different experimental signatures.
Specifically, if the scalar field has technically natural interactions with a
dark Yang-Mills sector, these interactions can give rise to additional
friction/dissipation channels for the kinetic energy of the scalar field.
These dissipation mechanisms can convert the coherent and homogenous kinetic
energy of the field into hot radiation. The energy density in this radiation
can be as large as a few percent of $\text{meV}^{4}$. While this is
considerably larger than the energy density set by $n_{\text{eff}}$
measurements from the CMB (in fact, this energy density is larger than that of
the CMB today), there is no conflict with the CMB since this radiation occurs
from the conversion of dark energy today rather than in the early universe.

If the dark energy were to have this phenomenology, there are considerable
experimental implications. First, these dissipation mechanisms can remove a
significant chunk of the kinetic energy of the field leading to null results
in some of the DC experimental probes discussed above. In particular the spin
precession measurements in laboratories today would have a greatly reduced
signal. However, searches for polarization rotation of the CMB (cosmic
birefringence) would still be important as even a tiny integrated kinetic
energy would be detectable beyond astrophysical limits. Second, the conversion
of the kinetic energy into hot radiation could be detected indirectly through
cosmological measurements of the expansion rate. Third, this conversion could
be detected directly by searching for this radiation to interact with a
detector. This requires new kinds of experiments. The hot radiation is at a
high frequency $\sim$ THz ($\sim$ meV) where systematic backgrounds are under
better control than a DC experiment. But, the signal is not coherent and thus
requires enhanced sensing.

The properties of the standard model imply that this radiation is not in the
form of photons. The high mass of the electron relative to the dark energy
scale as well as astrophysical limits on the direct coupling of light scalar
fields to electromagnetism makes the conversion into photons inefficient. In
this paper, we point out that this conversion can be efficient through Yang-
Mills sectors, hidden/dark photons, milli-charged particles and right handed
neutrinos.

The rest of this paper is organized as follows. In section 2, we outline the
basic technique that allows the kinetic energy of dark energy to be converted
into radiation in a technically natural way. We then show how this radiation
can manifest itself in the form of pure Yang-Mills gauge bosons (section 2.1),
hidden photons (section 2.2), milli-charged particles (section 2.2) and right
handed neutrinos (section 2.3). Following this, in section 3 we discuss
cosmological probes of this scenario and conclude in section 4 where we
comment on various direct detection possibilities.

##  2 Friction and Radiation

The point of this section is to show that dark energy could easily have a
radiation component consisting of a variety of possible particle types. We
focus on models in which a rolling scalar field couples to such a light
sector. The rolling scalar populates the light sector, which in turn
thermalizes and creates friction for the rolling scalar as a form of back-
reaction.

We take arguably the simplest model of a rolling scalar field, one with a
linear potential $V(\varphi)=-C\varphi$. The equation of motion for such a
field will be

$\ddot{\varphi}+3H\dot{\varphi}+\Upsilon\dot{\varphi}=C$ (1)

Where $H$ is the Hubble rate and $\Upsilon$ is the (generally temperature-
dependent) source of friction that we will exploit in this paper. The over-dot
denotes time derivatives in standard FRW coordinates. We have chosen the
origin of $\varphi$ to be the point of zero vacuum energy. Observation
requires $\varphi<0$ today (positive vacuum energy) and we will look for
solutions where $\ddot{\varphi}<(3H+\Upsilon)\dot{\varphi}$. Unless otherwise
noted, we will assume $\Upsilon\gg 3H$. Thus dark energy density will consist
of vacuum energy plus a quasi-thermal bath and a near-negligible contribution
from $\varphi$’s kinetic energy.

In the rest of this section, we present minimal models that can generate
different particle contents in the background dark radiation.

###  2.1 Non-Abelian Gauge Sector

For a non-Abelian gauge theory with fermions and an axion-like coupling of
$\varphi$ to the gauge fields, we have the Lagrangian

$-{\cal
L}=\frac{1}{2g^{2}}\>\mathrm{Tr}\>G_{\mu\nu}G^{\mu\nu}+\bar{\psi}\left(\not{D}+m\right)\psi+\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi+\frac{\varphi}{f}\frac{\>\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}-C\varphi$
(2)

where we take a gauge group $SU(N_{c})$ and the $N_{f}$ fermions $\psi$ are in
some representation $R$ of the gauge group. Contributions to dynamical
friction on a pseudo-scalar from a QCD-like sector has been of significant
interest as it affects the calculation of the axion dark-matter abundance
[14].

First we investigate the simplest case, pure Yang-Mills (essentially the
$m\rightarrow\infty$ limit). The field equation of motion is

$\partial_{\mu}\partial^{\mu}\varphi=-C+\frac{\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}f}.$
(3)

The last term is related to the Chern-Simons number density. When a thermal
average is taken of this equation, the last term will contribute to rolling
friction due to real-time sphaleron processes when $\dot{\varphi}\neq 0$. This
contribution is

$\left\langle\frac{\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}\right\rangle=\frac{\Gamma_{\mathrm{sph}}}{2T}\left(\frac{\dot{\varphi}}{f}\right)$
(4)

where the sphaleron rate is of order

$\Gamma_{\mathrm{sph}}\sim N_{c}^{5}\alpha^{5}T^{4}$ (5)

up to logs of $\alpha$. The factor $1/2T$ is the relation between the free
diffusion rate for topological charge $\Gamma_{\mathrm{sph}}$ and the rate at
which a chemical potential drives topological transitions, (see Ref. [15]).
For a discussion of the sphaleron rate see [16, 17, 18], and for a general
discussion of chemical potentials, particle number violation, and sphaleron
processes, see [19, 20, 21]. Thus the friction coefficient in (1) is [22]
$\Upsilon\sim(N_{c}\alpha)^{5}T^{3}/f^{2}$. This friction beats Hubble
friction when $f\ll T\sqrt{T/H}(N_{c}\alpha)^{5/2}$.

Now we can estimate the temperature as a function of the remaining parameters.
The dark radiation density, $\rho_{DR}$, will lose energy density via redshift
from the expanding universe, and gain energy density from the friction’s back
reaction:

$\dot{\rho}_{DR}=-4H\rho_{DR}+\Upsilon\dot{\varphi}^{2}$ (6)

where $\rho_{DR}=(2N_{c}^{2}-1)(\pi^{2}/30)T^{4}$, accounting for the gluons
and the thermalized axions (thermal excitations of the inflaton field) which
accompany them. We can plug in for $\dot{\varphi}$ by solving for it in (1) in
the slow-roll ($\ddot{\varphi}\rightarrow 0$) and thermal friction
($\Upsilon\gg H$) limits, giving $\dot{\varphi}\simeq C/\Upsilon$. If we also
assume quasi-steady-state behavior in (6), $\dot{\rho}_{DR}\sim 0$, we can
solve for $T$:

$T\sim\left(\frac{C^{2}f^{2}}{HN_{c}^{7}\alpha^{5}}\right)^{1/7}.$ (7)

The temperature’s slow ($-1/7$ power) dependence on $H$ shows this is a
reasonable approximation. However, in Section 3, we solve the exact equations
of motion analytically and numerically.

Note, we can choose parameters such that $T\sim$ meV. For thermal friction
dominance, this amounts to $f\ll$ TeV. As this is all occurring in a hidden
sector it is not subject to collider bounds. When one includes couplings to
standard model particles, there exists a broad range of parameter space
allowed by current experiments.

Also note, the following assumptions are required to make our analysis valid.
We assume that $\dot{\phi}$ is slowly varying on the $1/\alpha^{2}T$ time
scale on which the thermal medium reacts to these effects. If this is not the
case, then the response is more complicated and is not determined by the
sphaleron rate. However, this is a trivial constraint. We also assume that the
$\dot{\varphi}G\tilde{G}/f$ term is smaller than the other terms in the YM
equations of motion so we can linearize in it in the usual fluctuation-
dissipation sense. This happens as long as $\dot{\varphi}/f<T$, which is again
trivial to satisfy. Finally, we take $\alpha^{2}T>H$ so the Hubble rate is
small compared to the rate of thermalization, which simply requires $\alpha$
to not be too small.

Lastly, note that this quasi-steady-state behavior is relatively insensitive
to initial conditions. Setting $\dot{\varphi}=0$ and $T=0$ at early times
still asymptotically generates this friction and heat bath today. A large $T$
at early times will red-shift away and is anyway limited by CMB sensitivity to
additional degrees of freedom. A large $\dot{\varphi}$ at early times
redshifts even faster than radiation and is then damped by thermal dynamics.
Thus, the solution we find is an attractor.

####  2.1.1 Non-Abelian with Fermions

The existence of fermions in this sector can suppress the friction as
sphaleron transitions can build up chiral charge and associated chemical
potential, which suppresses the sphaleron rate. A fermion mass violates the
chiral symmetry and thus reduces the chiral charge. These issues only occur if
the fermion masses are very light; we will see below that, if the fermion mass
$m$ grows larger than $\alpha^{2}T$, the sphaleron rate approaches that of the
pure YM case.

If the mass is zero, there are two conserved fermionic numbers, the left-
handed fundamental number $n_{F}$ (for fermions in the fundamental
representation) and the left-handed anti-fundamental number $n_{\bar{F}}$,
with associated currents $j^{\mu}_{F}$ and $j^{\mu}_{\bar{F}}$. If we think of
the two Weyl fermions as one Dirac fermion, $n_{F}$ is the net left-particle
number and $n_{\bar{F}}$ is minus the net right-particle number. We will
assume initially $n_{F}=n_{\bar{F}}=0$. Again in the absence of a mass, each
current obeys

$\partial_{\mu}j^{\mu}_{F}=\frac{\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}=\partial_{\mu}j^{\mu}_{\bar{F}}$
(8)

where the field strength trace is the topological charge density. For a
general representation $R$ we would find

$\partial_{\mu}j^{\mu}_{R}=2T_{R}\frac{\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}$
(9)

where $T_{R}=C_{R}d_{R}/d_{A}$ is the trace normalization of the
representation ($1/2$ for fundamental, $N_{c}$ for adjoint, etc), $C_{R}$ is
the Casimir, $d_{R}$ is the dimension, and $d_{A}=N_{c}^{2}-1$ is the
dimension of the adjoint representation. Assuming a statistically uniform
system, this reduces to equations for particle number, such as

$\partial_{t}n_{R}=2T_{R}\frac{\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}.$
(10)

If the particle number is otherwise conserved (or for small mass, long lived),
then the fermions otherwise equilibrate, and the particle number is related to
a chemical potential $\mu_{R}$ via

$\displaystyle n_{R}(\mu)$
$\displaystyle=d_{R}\int\frac{d^{3}p}{(2\pi)^{3}}\left(\frac{1}{e^{(p-\mu_{R})/T}+1}-\frac{1}{e^{(p+\mu_{R})/T}+1}\right)\simeq\frac{d_{R}\mu_{R}T^{2}}{12}\,,$
(11) $\displaystyle\mu_{R}(n)$ $\displaystyle=\frac{12n_{R}}{d_{R}T^{2}}\,.$

Here the two terms in the first line are the particles and the antiparticles,
and we have expanded to linear order in $\mu_{F}$.

Now consider the equations of motion for the Chern-Simons number itself.
Besides the driving force $\dot{\varphi}$ acting on it, there is also the
energy cost of producing particle number against a chemical potential:

$\left\langle\frac{\mathrm{Tr}\>G_{\mu\nu}\tilde{G}^{\mu\nu}}{16\pi^{2}}\right\rangle=\frac{\Gamma_{\mathrm{sph}}}{2T}\left(\frac{\dot{\varphi}}{f}-\sum_{R}2T_{R}\mu_{R}\right)$
(12)

where we showed the general situation for multiple fields in different
representations $R$. The factor $2T_{R}$ is the number of particles produced
by a sphaleron event, each of which costs a free energy $\mu_{R}$.

Taking the thermal average of equations of motion (1) and (10), and inserting
the relations for the chemical potential (11), we find the following combined
equations (for a massless fermion in representation $R$):

$\displaystyle-\ddot{\varphi}=-C+\frac{\Gamma_{\mathrm{sph}}}{2Tf}\left(\frac{\dot{\varphi}}{f}-\frac{24T_{R}n_{R}}{d_{R}T^{2}}\right),$
(13)
$\displaystyle\partial_{t}n_{R}=\frac{T_{R}\Gamma_{\mathrm{sph}}}{T}\left(\frac{\dot{\varphi}}{f}-\frac{24T_{R}n_{R}}{d_{R}T^{2}}\right).$
(14)

We see that there is no solution where $\ddot{\varphi}$ vanishes
asymptotically; instead, the particle number grows and $\varphi$ continues to
accelerate. As a side note, there is still effective dissipation even in this
exactly massless case. These equations (13) and (14) can be solved for an
asymptotic solution which still has energy being dumped from the scalar field
into the thermal bath. This rate of energy transfer is highly suppressed and
we will not focus on these solutions.

However, in the case of a small non-zero fermion mass ($m\neq 0$), there are
chirality-violating scattering processes, with a matrix element proportional
to the mass at leading order. This means that there is an additional term in
the particle number evolution equation:

$\partial_{t}n_{R}=\frac{T_{R}\Gamma_{\mathrm{sph}}}{T}\left(\frac{\dot{\varphi}}{f}-\frac{24T_{R}n_{R}}{d_{R}T^{2}}\right)-\kappa\frac{n_{R}N_{c}\alpha
m^{2}}{T^{2}}$ (15)

with $\kappa$ a coefficient which one could calculate in principle and should
be ${\cal O}(1)$ [23]. Now we see that $n_{R}$ can approach a constant value.
Setting its time dependence in (15) to zero and solving for $n_{R}$, one can
now generate an effective friction rate $\Gamma_{\text{eff}}$ in the $\varphi$
equation of motion

$-\ddot{\varphi}=-C+\Gamma_{\text{eff}}\,\dot{\varphi}$ (16)

with

$\Gamma_{\text{eff}}=\frac{\Gamma_{\mathrm{sph}}}{2Tf^{2}}\left(\frac{\Gamma_{ch}}{\Gamma_{ch}+\frac{24T_{R}^{2}}{d_{R}T^{3}}\Gamma_{\mathrm{sph}}}\right)$
(17)

where we defined the chirality-violation rate $\Gamma_{ch}\equiv\kappa
N_{c}\alpha m^{2}/T$. If the $\Gamma_{ch}$ term is larger than the sphaleron
term – i.e., if $m^{2}/T^{2}\gg(24T_{R}^{2}/d_{R}\kappa)(N_{c}\alpha)^{4}$ –
then the particle number relaxes primarily due to spin flips before it can
back-react on the sphaleron processes, and thus for a large enough mass, the
original friction is restored. For a finite mass, friction at some level is
maintained and steady-state solutions are possible. In this limit where
$\Gamma_{ch}$ is small, the thermal friction beats Hubble friction if
$\Gamma_{\text{eff}}\simeq(\kappa
d_{R}N_{c}\alpha)/(48T_{R}^{2})(Tm^{2}/f^{2})>H$, and thus it is possible to
have significant dark energy radiation for a wide range of masses.

###  2.2 Dark Photons and Milli-Charged Particles

There are a number of ways to make $U(1)$ gauge bosons (dark photons) part of
the dark energy thermal bath. The simplest way is to take the QCD-like model
above and charge the fermions under an unbroken $U(1)$. Take the new fine-
structure constant to be $\alpha^{\prime}$. In the thermal bath, the rate of
dark photon production should be at least
$\Gamma_{\gamma^{\prime}}\sim\alpha\alpha^{\prime}T$. As long as this is much
larger than the Hubble expansion rate, a thermal abundance of dark photons
would be generated.

One does not need fermions in the non-Abelian sector to thermalize dark
photons. Additional ways to thermalize light $U(1)$s is through higher-
dimensional operators e.g.,

${\cal
L}\supset\frac{\varphi}{f^{\prime}}F^{\prime\mu\nu}{\tilde{F}^{\prime}}_{\mu\nu}$
(18)

where $F^{\prime}$ is the field strength of a dark photon. Because the YM to
dark photon cross section now goes as $T^{2}/(f^{\prime}f)^{2}$, a much lower
$f^{\prime},f$ would be required to thermalize this sector, roughly
$f^{\prime}\sim f\lesssim 10$ keV for $T\sim$ meV. The only other requirement
is that the mass of the dark photon, $m_{A^{\prime}}$ satisfies
$T>m_{A^{\prime}}\gg H$ (we discuss the massless case below).

This thermalization process can be sped up if there are interactions in the
$U(1)$ sector. This is because when thermal friction dominates, the rolling
dark energy field is all but guaranteed to thermalize as well. Roughly, if the
rate $\Gamma_{GG\rightarrow G\varphi}\simeq\sigma n_{G}\sim(\alpha
T^{3}/f^{2})>H$, then $\varphi$ particles are also thermal. Once a small
amount of dark photons is produced via (18), a runaway process can produce a
thermal abundance in this sector. For example, a sector of dark photons
$A^{\prime}_{\mu}$ and dark-charged fermions $\chi$ will have scattering cross
sections $\sigma_{\varphi\chi\rightarrow A^{\prime}\chi}$ and $\sigma_{\varphi
A^{\prime}\rightarrow\chi\chi}$ of order $\alpha^{\prime}/f^{\prime 2}$. The
growth rate of the number density of hidden sector particles $X\equiv
A^{\prime},\chi$ would be

$\frac{dn_{X}}{dt}\sim n_{X}n_{\varphi}\sigma_{\varphi X\rightarrow XX}\sim
n_{X}\frac{\alpha^{\prime}T^{3}}{f^{\prime 2}}.$ (19)

If the characteristic rate $\frac{\alpha^{\prime}T^{3}}{f^{\prime 2}}\gg H$
(up to a log factor), a thermal abundance of dark photons and dark-charged
fermions should be produced. If the dark photons kinetic-mix with regular
photons, this would also constitute an abundance of milli-charged fermions.

For completeness, we note that in the case of an extremely light or massless
dark photon without charged fermions, an abundance of very long wavelength
modes is also generated by a tachyonic instability, and in fact becomes a new
source of friction for $\varphi$ [24]. Via the coupling (18), a population of
vectors will be produced when $\dot{\phi}/f>m_{A^{\prime}}$ due to an
instability in the mode equation for the vectors [24, 25, 26]:

$\ddot{A^{\prime}_{\pm}}+\left(k^{2}\pm
4k\frac{\dot{\varphi}}{f}+m_{A^{\prime}}^{2}\right)A^{\prime}_{\pm}=0$ (20)

For very light or massless dark photons, with no charged fermions, this effect
dominates the friction mechanism and generates a background of non-thermal,
very low frequency ($\sim\dot{\varphi}/f^{\prime}$) vectors. Using the work of
[25], the field velocity is roughly constant at $\dot{\varphi}=\xi
Hf^{\prime}$, with $\xi\sim{\cal O}(100)$, corresponding to dark photon
frequencies around a hundred times the Hubble rate. The energy density in the
dark photon approaches the vacuum energy $V$ when $C\sim V/(\xi f^{\prime})$
[25]. This friction will dominate over the thermal friction from the YM sector
unless $f^{\prime}\gg\sqrt{T/H}f$, a limit in which the dark photons don’t
thermalize anyway. On the other hand, a full analysis of this version of dark
energy radiation friction should run parallel to the inflation version of this
model analyzed in [25], and we suggest its cosmological analysis as dark
energy for future work.

###  2.3 Neutrinos

One very interesting possibility is the dark-energy thermalization of Standard
Model neutrinos. This could happen through the coupling of right-handed
neutrinos to the dark energy radiation. For example, if fermions in the
general Lagrangian (2) are in the adjoint representation of the non-Abelian
group, this would allow for a simple way to thermalize neutral fermions, $N$,
through the coupling

${\cal L}\supset\frac{1}{f_{N}}G^{a}_{\mu\nu}\psi^{a}\sigma^{\mu\nu}N$ (21)

If $N$ is also a right-handed neutrino via the coupling $yh\ell N$, then it
will thermalize the Standard Model neutrinos with masses $m_{\nu}\lesssim T$,
potentially making them hotter than the standard temperature prediction for
the CMB$\nu$ background. While thermalization at late times will be generated
as long as $f_{N}\ll$ TeV, this sector will be frozen out at early times and
could easily have played a negligible role during Big Bang Nucleosynthesis and
recombination.

##  3 Cosmological Signatures

Thermal friction has a unique cosmological signature through its modification
of the late-time behavior of the dark-energy equation of state $w(a)$ as a
function of the scale factor of the universe, $a$. The dynamics also gives
rise to late-time spatial inhomogeneity (see [27]) with small deviations from
quintessence predictions, though these do not produce sizeable observables.
However, the background level evolution has distinctive features from a
cosmological constant $w_{CC}(a)=-1$, as well as other scalar field models,
which may be observable with upcoming experiments. In this section we present
the derivation of the thermal friction dark energy equation of state
$w_{\text{DER}}(a)$.

The equations of motion of the scalar field-radiation system in equation (1)
and (6) as a function of the scale factor $a$ are:

$a^{2}H^{2}\varphi^{\prime\prime}+\left(4aH^{2}+a^{2}H^{\prime}H+aH\Upsilon\right)\varphi^{\prime}=C$
(22)
$\frac{4\pi^{2}}{30}g_{*}aHT^{\prime}T^{3}+\frac{4\pi^{2}}{30}g_{*}HT^{4}=\Upsilon
a^{2}H^{2}{\varphi^{\prime}}^{2}$ (23)

where primes denote derivatives with respect to $a$ and $g_{*}$ is the
effective number of degrees of freedom in the dark radiation sector. We will
take $a=a_{0}\equiv 1$ to denote the scale factor today. Assuming the friction
coefficient of a pure non-Abelian gauge group $\Upsilon=T^{3}/\tilde{f}^{2}$
(where any dependence on $\alpha$ and $N_{c}$ has been absorbed in
$\tilde{f}$), one can now solve this system of equations for $\varphi(a)$ and
$T(a)$ with appropriate boundary conditions. However, focusing on the regime
in which thermal friction is dominant ($H\ll T^{3}/\tilde{f}^{2}$) at all
relevant times simplifies equations (22) and (23) to:

$\varphi^{\prime}(a)\approx\frac{\tilde{f}^{2}C}{aH(a)T^{3}(a)}$ (24)
$T^{\prime}(a)+\frac{T(a)}{a}=\frac{15aH(a){\varphi^{\prime}(a)}^{2}}{2\pi^{2}g_{*}\tilde{f}^{2}}$
(25)

where $g_{*}$ denotes the relativistic degrees of freedom in the radiation
energy density $\rho_{\text{DR}}\equiv\frac{\pi^{2}}{30}g_{*}T^{4}$. We solve
equation (25) perturbatively by approximating the Hubble function by its
counterpart with a constant dark energy, $H(a)\approx H^{(0)}(a)\equiv
H_{0}\sqrt{\frac{\Omega_{m}}{a^{3}}+1-\Omega_{m}}$, where $\Omega_{m}$ is the
ratio of matter density to critical energy density today (we are assuming a
flat geometry). In doing so, we obtain the following analytic expression for
the temperature:

$T(a)=\left(\frac{\frac{15C^{2}\tilde{f}^{2}}{2\pi^{2}g_{*}}\sqrt{\frac{\Omega_{m}}{a^{3}}+1-\Omega_{m}}\left(8a^{3}(1-\Omega_{m})-11\Omega_{m}\left(1-{}_{2}F_{1}\left[1,\frac{4}{3};\frac{11}{6};-\frac{a^{3}(1-\Omega_{m})}{\Omega_{m}}\right]\right)\right)}{8a^{3}H_{0}(1-\Omega_{m})^{2}}\right)^{\frac{1}{7}}$
(26)

which we plug into equation (24) and numerically integrate to solve for
$\varphi(a)$. Here, ${}_{2}F_{1}[a,b;c;z]$ is the Gaussian hypergeometric
function. In this regime the kinetic energy contribution is highly suppressed.
Thus, in this limit $T(a)$ and $\varphi(a)$, and their values today
$T_{0}\equiv T(a=1)$ and $\varphi_{0}\equiv\varphi(a=1)$, quantify the dark-
energy equation of state $w_{\text{DER}}$ as a function of $a$:

$w_{\text{DER}}(a)\simeq-1\frac{-C\varphi(a)}{-C\varphi(a)+\frac{\pi^{2}}{30}g_{*}T^{4}(a)}+\frac{1}{3}\frac{\frac{\pi^{2}}{30}g_{*}T^{4}(a)}{-C\varphi(a)+\frac{\pi^{2}}{30}g_{*}T^{4}(a)}.$
(27)

The value of the Hubble constant today, $H^{2}_{0}$, and the temperature
today, $T_{0}$ (which could potentially be measured by direct detection
discussed in Section 4), determines the parameter combinations $C\varphi_{0}$
and $C\tilde{f}$, the only combinations appearing in the dark-energy equation
of state itself (up to the discrete $g_{*}$). Since the constraints on
$w_{\text{DER}}$ are inferred assuming a $\Lambda$CDM cosmology, a dedicated
analysis is required to evaluate its bounds assuming a dynamic dark energy
with thermal friction, which is an interesting prospect beyond the scope of
the discussion presented here. For illustrative purposes, we plot
$w_{\text{DER}}(a_{0})\leq-0.95$ in Figure 1 which corresponds to temperatures
of about $T\lesssim 0.8\,\text{meV}$, assuming $g_{*}=7$ (the case of $SU(2)$
with $\varphi$ also thermalized). Though there are three model parameters, the
slope $C$, the value of the scalar field today $\phi_{0}$, and $\tilde{f}$ the
friction strength parameter, there is no phenomenological sensitivity to all
three parameters in this regime. This is due to
$w_{\text{DER}}^{\prime}(a_{0})$ being dominated by
$T^{\prime}(a_{0})\propto\left(C\tilde{f}\right)^{\frac{2}{7}}$ rather than
$\varphi^{\prime}(a_{0})$ which does not lift the degeneracy between different
combinations of $C$ and $\tilde{f}$. Thus the shape of the curves
$w_{\text{DER}}(a,T_{0})$ is fully determined by $w(a_{0})$ as shown on the
left of Figure 1.

Figure 1: The dark energy equation of state $w$ as a function of scale factor
$a$ (redshift $z$) for thermal friction on the left and quintessence on the
right. The abundances are defined as:
$\Omega_{\varphi}=\frac{-C\varphi_{0}}{\rho_{c}}$,
$\Omega_{\text{DR}}=\frac{\frac{\pi^{2}}{30}g_{*}T^{4}_{0}}{\rho_{c}}$,
$\Omega_{\dot{\varphi}_{Q}}=\frac{\frac{1}{2}a_{0}^{2}H^{2}_{0}{\varphi^{\prime}(a_{0})}^{2}}{\rho_{c}}$
where $\rho_{c}$ is the critical density today. We assume $g_{*}=7$
corresponding to a minimal thermal bath (comprised of three gauge bosons with
two polarizations each and the axion with one scalar degree of freedom), and
the best fit Planck value $\Omega_{m}=0.315$ [28].

In the opposite limit $H\gg\Upsilon$ ($\tilde{f}\rightarrow\infty$), we
recover the familiar quintessence solution in which kinetic energy dominates
and dark radiation is highly suppressed. To solve for the quintessence field
$\varphi_{Q}$ we turn off the thermal friction and recover the equation of
motion of the regular quintessence scalar field whose terms were previously
sub-dominant:

$a^{2}H^{2}\varphi_{Q}^{\prime\prime}+\left(4aH^{2}+a^{2}H^{\prime}H\right)\varphi_{Q}^{\prime}=C$
(28)

and find an equation of state:

$w_{Q}(a)\simeq-1\frac{-C\varphi_{Q}(a)}{-C\varphi_{Q}(a)+\frac{1}{2}a^{2}H(a)^{2}{\varphi_{Q}^{\prime}(a)}^{2}}+\frac{\frac{1}{2}a^{2}H(a)^{2}{\varphi_{Q}^{\prime}(a)}^{2}}{-C\varphi_{Q}(a)+\frac{1}{2}a^{2}H(a)^{2}{\varphi_{Q}^{\prime}(a)}^{2}}$
(29)

Using the Ansatz $H^{(0)}(a)$ (see text above equation (26)) we obtain the
analytical solution for the quintessence field $\varphi_{Q}(a)$:

$\varphi_{Q}(a)=\varphi_{Q,0}+\frac{2C\left(\text{tanh}^{-1}\left(\frac{\sqrt{1-\Omega_{m}}}{\sqrt{\frac{\Omega_{m}}{a^{3}}+1-\Omega_{m}}}\right)-\text{arctanh}^{-1}\left(\sqrt{1-\Omega_{m}}\right)\right)}{9H^{2}_{0}(1-\Omega_{m})^{\frac{3}{2}}}$
(30)

where $\varphi_{Q,0}\equiv\varphi_{Q}(a_{0})$. To the right side of Figure 1
we show the characteristic shape of quintessence curves with varying kinetic
energy contributions. At early times $a<0.3$, quintessence closely resembles a
cosmological constant. Only recently did the field obtain a sizeable kinetic
energy component $\Omega_{\dot{\varphi}}$ whose total amount is smaller than
the amount of dark radiation in the thermal friction model that results in the
same value for $w(a_{0})$ today. On the contrary, the change of
$w_{\text{DER}}(a)$ in the thermal friction model is more gradual. Even at
early times there exists a sizeable deviation from a cosmological constant.
The slowly rising temperature increases this deviation which asympotatically
approaches a constant value. Quintessence shows a differing trend with the
deviation from a cosmological constant linearly growing. Quintessence also
requires much larger field values venturing into the trans-planckian field
regime in order to remain within $w(a_{0})<-0.95$, whereas thermal friction
has a wide set of viable parameters that are sub-Planckian. For example the
gray dashed line on the left of Figure 1 corresponds to
$\varphi(a_{0})\sim-0.01\left(\frac{\tilde{f}}{1\text{ GeV}}\right)M_{pl}$,
whereas the gray dashed line on the right of Figure 1 corresponds to
$\varphi_{Q}(a_{0})\sim-2M_{pl}$. This is due to the thermal friction being
much larger than Hubble friction which allows for much steeper slopes within
the experimental constraints leading to smaller field values. The perturbative
descriptions of the two regimes discussed in this section are in very good
agreement with exact numerical solutions of (22) and (23), with deviations
smaller than $0.02\%$ for all curves shown in Figure 1.

As we have seen above, in our scenario there is actually only a single
additional parameter beyond the usual $\Lambda$CDM parameters which will
determine our entire equation of state as a function of time, $w_{DER}(a)$, at
leading order. Further, this also determines the temperature of the dark
radiation today, which could be separately measured in a direct detection
experiment. If we find evidence that dark energy is not a cosmological
constant then we will immediately want to measure its variation as a function
of cosmic time. If the (even roughly) measured $w(a)$ curve fits with our
constrained family of predictions, it would be a striking success of the
model. Further, it would then make a sharp prediction for the temperature of
the dark radiation which would start a race to directly detect this radiation
in the lab.

##  4 Discussion and Conclusion

The mechanisms discussed above can produce a thermalized local energy density
in dark radiation that can be as large as a few percent of the dark energy
density $\sim\text{meV}^{4}$ today, much larger than the energy density in the
CMB. This radiation would have a potentially observable effect on the
expansion rate of the universe, distinguishable from just a rolling scalar
field (quintessence). Interestingly, even though most of the kinetic energy of
the dark energy field has gone into radiation, the usual cosmic birefringence
signal could still be observed at the same level. The radiation might also be
directly detectable. This radiation may be in the form of pure Yang Mills
gluons, hidden photons, right handed neutrinos or milli-charged particles. It
would be interesting to develop experimental methods to directly detect these
particles, similar to dark matter direct detection efforts.

Of these types of radiation, a pure Yang Mills sectors is the most difficult
to detect since it can only interact with the nucleons and electrons through
operators of high dimension. The other portals are more promising. There are
many experimental efforts [29, 30, 31, 32, 33, 34, 35] currently underway to
detect hidden photons. These experiments are focused on detecting cold dark
matter at lower frequencies, with the experiment utilizing the expected
narrow-band nature of the dark matter signal to resonantly boost sensitivity.
But, there are emerging attempts to detect heavier hidden photons in the meV
range. In this energy range, due to the difficulty of creating electromagnetic
resonators, dark matter direct detection efforts are also broadband devices.
Moreover, the fundamental sensing technology necessary to make progress in
this energy range is the technology necessary to detect single photons in the
THz band. Due to the practical importance of this technology, it is likely
that there will be continued investment in developing these sensors. This will
lead to direct improvement in our ability to search for hidden photon dark
radiation produced by dark energy. Similarly, experiments have also been
proposed to detect milli-charged dark matter. It would be interesting to
evaluate the reach of these experiments to relativistic milli-charged
particles. While relavitistic particles will spend less time transiting a
detector than non-relativistic dark matter, they also interact more strongly
with magnetic fields and this could potentially be used to create magnetic
focusing devices to boost the experimental reach. We leave a detailed
discussion of this topic for future work.

The case of right handed neutrinos is particularly interesting. At the energy
scale $\sim$ meV relevant to dark energy radiation, the right handed neutrino
easily mixes with the left handed neutrino and can thus interact with the
standard model through the weak interactions. The right handed neutrinos
produced from dark energy could thus have a much larger energy density (by
$\sim$ four orders of magnitude) than the standard model (left handed)
neutrino background. Due to strong mixing, these right handed neutrinos can
also be searched for using experiments such as PTOLEMY [36] that are dedicated
to probing the cosmic neutrino background. In these experiments, the dark
energy produced neutrinos may have an energy density 1000x larger than the
standard model neutrino background giving rise to an enhanced rate in the
experiment. Moreover, these neutrinos could also be at higher energies $\sim$
meV as opposed to $\sim 0.1$ meV expected from the cosmic neutrino background.
The increased energy should help in detection in an experiment like PTOLEMY
that searches for a kinematic end point.

A discovery of meV scale radiation in an experiment would be proof that this
radiation was produced at late times in the universe well after the CMB. The
two contenders for the source of such a signal would be dark matter and dark
energy. While the meV scale would be indicative of dark energy, it will be
necessary to develop experimental methods to distinguish these two
possibilities. The first check would be to probe the equation of state of dark
energy through cosmological measurements. With the advent of gravitational
wave detectors that offers the tantalizing possibility of genuine standard
candles, it might be possible to extend the reach of such cosmological probes.
Further, once it is known that the radiation couples to the Standard Model,
then the scalar dark energy field giving rise to the radiation must couple as
well at some level. Interestingly, this rolling scalar field could then give a
detectable cosmic birefringence signal in the CMB, which would be another
indication for this scenario.

This paper has focused on the conversion of dark energy into thermal radiation
through technically natural couplings between a rolling field and a gauge
sector. These kinds of couplings can also source another class of signals.
Tachyonic instability [25] can result in the conversion of the kinetic energy
of the field into ultra-long ($\sim 100/H_{0}$) wavelength modes of a vector
field. Canonical examples of such vector fields are $B-L$ gauge bosons and a
hidden photon. Signals from the former can potentially be looked for using
torsion balances/accelerometers[37]. The latter case may be phenomenologically
interesting - given the plasma mass of the photon, the hidden photon field is
likely to manifest itself as a spatially coherent magnetic field in the
universe and it might be interesting to see if such a magnetic field could be
the source of the observed long range coherence in extra-galactic magnetic
fields [38].

The seriousness of the cosmological constant problem and the most reasonable
approach to its solution makes a strong case for dynamical dark energy. Dark
matter also presents a deep mystery about the contents of our universe. While
dark matter was discovered through its gravitational effects, theoretical
motivations suggested a robust experimental program to directly detect it in
the laboratory. Dark energy too was discovered through its gravitational
impact on our universe. In this paper we have argued that theoretical
considerations also motivate a broader experimental program to directly detect
a dynamical component of dark energy in the laboratory.

## Acknowledgments

DK and SR are supported in part by the NSF under grant PHY-1818899. SR
acknowledges travel support from the Gordon and Betty Moore Foundation and the
American Physical Society enabling him to visit Stanford to complete this
work. SR is also supported by the DoE under a QuantISED grant for MAGIS. KB
was supported by in part by the NSF under grant PHY-1818899 and PHYS-1915093,
and the NASA award NNX17AK38G. SR and PWG are also supported by the SQMS
quantum center. PWG acknowledges the support provided by NSF Grant
PHY-1720397, the Heising-Simons Foundation Grant 2018-0765, the DOE HEP
QuantISED award 100495, and the Gordon and Betty Moore Foundation Grant
GBMF7946.

## References

  * [1] L.F. Abbott.  A Mechanism for Reducing the Value of the Cosmological Constant.  Phys. Lett. B, 150:427–430, 1985.
  * [2] Tom Banks.  Relaxation of the Cosmological Constant.  Phys. Rev. Lett., 52:1461–1463, 1984.
  * [3] Paul J. Steinhardt and Neil Turok.  Cosmic evolution in a cyclic universe.  Phys. Rev. D, 65:126003, 2002.
  * [4] Lasma Alberte, Paolo Creminelli, Andrei Khmelnitsky, David Pirtskhalava, and Enrico Trincherini.  Relaxing the Cosmological Constant: a Proof of Concept.  JHEP, 12:022, 2016.
  * [5] Peter W. Graham, David E. Kaplan, and Surjeet Rajendran.  Relaxation of the Cosmological Constant.  Phys. Rev. D, 100(1):015048, 2019.
  * [6] Peter W. Graham, David E. Kaplan, and Surjeet Rajendran.  Born again universe.  Phys. Rev. D, 97(4):044003, 2018.
  * [7] Latham A. Boyle, Robert R. Caldwell, and Marc Kamionkowski.  Spintessence! New models for dark matter and dark energy.  Phys. Lett. B, 545:17–22, 2002.
  * [8] Sean M. Carroll, George B. Field, and Roman Jackiw.  Limits on a Lorentz and Parity Violating Modification of Electrodynamics.  Phys. Rev. D, 41:1231, 1990.
  * [9] Yuto Minami and Eiichiro Komatsu.  New Extraction of the Cosmic Birefringence from the Planck 2018 Polarization Data.  Phys. Rev. Lett., 125(22):221301, 2020.
  * [10] M. Pospelov and M. Romalis.  Lorentz invariance on trial.  Phys. Today, 57N7:40–46, 2004.
  * [11] Michael V. Romalis and Robert R. Caldwell.  Laboratory search for a quintessence field.  2 2013.
  * [12] J.M. Brown, S.J. Smullin, T.W. Kornack, and M.V. Romalis.  New limit on Lorentz and CPT-violating neutron spin interactions.  Phys. Rev. Lett., 105:151604, 2010.
  * [13] Peter W. Graham, Selcuk Haciomeroglu, David E. Kaplan, Zhanibek Omarov, Surjeet Rajendran, and Yannis K. Semertzidis.  Storage Ring Probes of Dark Matter and Dark Energy.  5 2020.
  * [14] Larry D. McLerran, Emil Mottola, and Mikhail E. Shaposhnikov.  Sphalerons and Axion Dynamics in High Temperature {QCD}.  Phys. Rev. D, 43:2027–2035, 1991.
  * [15] Guy D. Moore.  Motion of Chern-Simons number at high temperatures under a chemical potential.  Nucl. Phys. B, 480:657–688, 1996.
  * [16] Peter B. Arnold, Dam Son, and Laurence G. Yaffe.  The Hot baryon violation rate is ${\cal O}(\alpha_{w}^{5}T^{4})$.  Phys. Rev. D, 55:6264–6273, 1997.
  * [17] Dietrich Bödeker.  On the effective dynamics of soft nonAbelian gauge fields at finite temperature.  Phys. Lett. B, 426:351–360, 1998.
  * [18] Guy D. Moore and Marcus Tassler.  The Sphaleron Rate in SU(N) Gauge Theory.  JHEP, 02:105, 2011.
  * [19] V.A. Kuzmin, V.A. Rubakov, and M.E. Shaposhnikov.  On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe.  Phys. Lett. B, 155:36, 1985.
  * [20] G.F. Giudice and Mikhail E. Shaposhnikov.  Strong sphalerons and electroweak baryogenesis.  Phys. Lett. B, 326:118–124, 1994.
  * [21] V.A. Rubakov and M.E. Shaposhnikov.  Electroweak baryon number nonconservation in the early universe and in high-energy collisions.  Usp. Fiz. Nauk, 166:493–537, 1996.
  * [22] Kim V. Berghaus, Peter W. Graham, and David E. Kaplan.  Minimal Warm Inflation.  JCAP, 03:034, 2020.
  * [23] A. Boyarsky, V. Cheianov, O. Ruchayskiy, and O. Sobol.  Evolution of the primordial axial charge across cosmic times.  7 2020.
  * [24] Bruce A. Campbell, Nemanja Kaloper, Richard Madden, and Keith A. Olive.  Physical properties of four-dimensional superstring gravity black hole solutions.  Nucl. Phys. B, 399:137–168, 1993.
  * [25] Mohamed M. Anber and Lorenzo Sorbo.  Naturally inflating on steep potentials through electromagnetic dissipation.  Phys. Rev. D, 81:043534, 2010.
  * [26] Anson Hook and Gustavo Marques-Tavares.  Relaxation from particle production.  JHEP, 12:101, 2016.
  * [27] Kim V. Berghaus and Karwal Tanvi.  Can thermal friction ease the Hubble Tension and the Large-Scale-Structure Tension, in preparation.  2021\.
  * [28] N. Aghanim et al.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys., 641:A6, 2020.
  * [29] Saptarshi Chaudhuri, Peter W. Graham, Kent Irwin, Jeremy Mardon, Surjeet Rajendran, and Yue Zhao.  Radio for hidden-photon dark matter detection.  Phys. Rev. D, 92(7):075012, 2015.
  * [30] Philip C. Bunting, Giorgio Gratta, Tom Melia, and Surjeet Rajendran.  Magnetic Bubble Chambers and Sub-GeV Dark Matter Direct Detection.  Phys. Rev. D, 95(9):095001, 2017.
  * [31] Gianpalo Carosi, Raphael Cervantes, Seth Kimes, Parashar Mohapatra, Rich Ottens, and Gray Rybka.  Orpheus: Extending the ADMX QCD Dark-Matter Axion Search to Higher Masses.  Springer Proc. Phys., 245:169–175, 2020.
  * [32] Simon Knapen, Tongyan Lin, Matt Pyle, and Kathryn M. Zurek.  Detection of Light Dark Matter With Optical Phonons in Polar Materials.  Phys. Lett. B, 785:386–390, 2018.
  * [33] Rouven Essig, Jesús Pérez-Ríos, Harikrishnan Ramani, and Oren Slone.  Direct Detection of Spin-(In)dependent Nuclear Scattering of Sub-GeV Dark Matter Using Molecular Excitations.  Phys. Rev. Research., 1:033105, 2019.
  * [34] Asimina Arvanitaki, Savas Dimopoulos, and Ken Van Tilburg.  Resonant absorption of bosonic dark matter in molecules.  Phys. Rev. X, 8(4):041001, 2018.
  * [35] Hao Chen, Rupak Mahapatra, Glenn Agnolet, Michael Nippe, Minjie Lu, Philip C. Bunting, Tom Melia, Surjeet Rajendran, Giorgio Gratta, and Jeffrey Long.  Quantum Detection using Magnetic Avalanches in Single-Molecule Magnets.  2 2020.
  * [36] E. Baracchini et al.  PTOLEMY: A Proposal for Thermal Relic Detection of Massive Neutrinos and Directional Detection of MeV Dark Matter.  8 2018.
  * [37] Peter W. Graham, David E. Kaplan, Jeremy Mardon, Surjeet Rajendran, and William A. Terrano.  Dark Matter Direct Detection with Accelerometers.  Phys. Rev. D, 93(7):075029, 2016.
  * [38] Kohei Kamada, Yuhsin Tsai, and Tanmay Vachaspati.  Magnetic Field Transfer From A Hidden Sector.  Phys. Rev. D, 98:043501, 2018.