# Critical-time metric for risk analysis against sharp input anomalies:
computation and application case study

Arthur Perodou  Christophe Combastel  Ali Zolghadri

###### Abstract

This paper investigates the critical-time criteria as a security metric for
controlled systems subject to sharp input anomalies (attack, fault),
characterized by having high impact in a reduced amount of time (e.g. denial-
of-service, attack by upper saturation). The critical-time is the maximal
time-horizon for which a system can be considered to be safe after the
occurrence of an anomaly. This metric is expected to be useful for risk
analysis and treatment (prevention, detection, mitigation). In this work, the
computational problem of the critical-time for uncertain linear systems and
several classes of sharp input anomalies, depending on the input channel and
the set of abnormal signal values, is formulated based on the quadratic
constraints (QC) framework, representing sets by the intersection of QC
inequalities and equalities. An iterative LMI-based algorithm is then proposed
to provide an under-estimate of the critical-time. Finally, the potential of
the critical-time as a metric for defense design is illustrated and discussed
on the quadruple-tank case study through different relevant scenarios.

###### keywords:

Secure Control , Cyber-Physical Systems , Critical-time , Sharp Anomalies ,
Quadratic Constraints

††journal:  Journal of XXX

[label1]organization=Ampere-lab, Univ. Lyon, Ecole Centrale de Lyon, INSA
Lyon, Univ. Claude Bernard Lyon 1, CNRS UMR 5005,city=Ecully, postcode=69130,
country=France

[label2]organization=IMS-lab, Univ. Bordeaux, CNRS UMR 5218,city=Talence,
postcode=33400, country=France

##  1 Introduction

The ever-increasing integration of information and communication technologies
into engineering systems has led to the emergence of the cyber-physical
systems (CPS) vision [1], where cyber and physical parts are jointly co-
designed. The Achilles heel of any CPS is widely recognized to be its
fragility, i.e. its vulnerability to adverse anomalous events such as physical
faults and cyber intrusions that may result in catastrophic cascading
consequences and unacceptable scenarios [2]. Although historically a topic of
the Computer Science community, these anomalous events have highlighted the
need to consider the interconnection between the cyber and physical parts, and
it is now recognized that defense design cannot be based solely on traditional
IT security methods [3]. Therefore, over the past decade, the System and
Control community has developed numerous methods to address the problem of
security in controlled systems [2, 4, 5].

A key point for CPS security is to develop metrics that can be used
efficiently for risk analysis and treatment (prevention, detection,
mitigation) [6, 4]. In that sense, several criteria have been developed in
order to quantify the physical impact of cyber-physical attacks. For instance,
the maximum impact of stealthy attacks on the quadratic objective function of
an optimal control problem is computed and used for the design of the
controller and the detector in [7, 8, 9]; the distance of the set of reachable
states by stealthy attacks to a safe/critical region is evaluated and used as
a prevention tool in [6, 10] or as criteria for controller and detector design
in [11]; the maximum number of attacks on actuators and sensors under which
the state can still be estimated is proposed jointly with an associated
mitigation scheme in [12]; a metric defined by the ratio of the energy
required by the attack to deviate the system from the state trajectory over
the control energy needed for recovery is reported in [13].

Recently, a new security metric named critical-time was proposed in [14]. The
critical-time is the maximal time-horizon for which a system is considered to
be safe after the occurrence of an anomaly (attack, fault), that is the system
is not in a critical state allowing to recover to a nominal mode or to
reconfigure for a safe stop. The interest of the critical-time is manifold.
First, this time-based metric can be usefully included in each stage of the
risk management framework. In the risk analysis stage, it identifies system
vulnerabilities and anomalies with the greatest impact, associated with a low
critical-time. In the risk treatment stage, this metric may be used to prevent
anomalies, associated then with high or even infinite critical-time [10], to
provide a design objective for the synthesis of anomaly detectors [15, 16] or
for the development of a mitigation mechanism, such as reconfigurable control
[17], dissipation block activation [18] or software rejuvenation [19].
Moreover, maximizing the critical-time will leave extra-time for defense
mechanisms, including human operators, to detect and mitigate anomalies. This
may be achieved by tuning remaining degrees of freedom, as in a rational
security allocation approach [4]. Furthermore, an originality is to consider
the more general concept of anomaly, instead of only attacks, as the impact of
certain attacks and faults can be modeled in a similar way (see for instance
[20]). This especially allows this metric to be integrated in a secure and
safe approach where both attacks and faults are considered. In particular, the
critical-time generalizes the detection and reconfiguration delays classically
used in the fault-management literature [15, 17].

In this paper, the critical-time metric for controlled systems subject to
sharp input anomalies, characterized by having high impact in a reduced amount
of time (e.g. denial-of-service, attack by upper saturation), is investigated.
The associated computational problem for an uncertain linear discrete-time
model is addressed by taking advantage from the quadratic constraints (QC)
framework, representing sets by the intersection of QC inequalities and
equalities. This problem is reformulated in terms of LMI (Linear Matrix
Inequality) [21] and an iterative LMI-based algorithm is proposed to provide
an under-estimate of the critical-time. Unlike [14], limited to the worst-case
scenario where the attacker has full control of the inputs, the class of sharp
input anomalies is significantly extended. By explicitly modeling the
controller, sharp anomalies that only impact some inputs, and not all of them,
can be considered. In addition, the inclusion of quadratic constraint
equalities (QCE) allows for considering specific sharp anomalies, such as
denial-of-service (DoS) or attack by upper saturation. Moreover, this paper
relaxes the assumption on the initial state which is assumed now to be
uncertain, leading to a more practical and robust result. Finally, the
potential of the critical-time metric for defense design is thoroughly
illustrated on the quadruple-tank case study [22] through different scenarios.
In particular, it is shown how the critical time can be maximized by properly
selecting the remaining degrees of freedom, leaving extra-time for defense
mechanisms.

The paper is organized as follows. In Section 2, the QC framework for convex
set representations is introduced. In Section 3, the critical-time computation
problem of a controlled system, with an uncertain discrete-time linear model,
subject to sharp input anomalies is formulated. Section 4 demonstrates how the
computation can be performed based on LMI optimization. Finally, the approach
is applied on the quadruple-tank case study in Section 5, while Section 6
provides some concluding remarks.

##### Notations

Lower (upper) case letters are used for vectors (matrices). $\mathbb{R}^{+}$
denotes the set of non-negative real numbers and $\mathbb{R}^{n\times m}$ the
set of real-valued matrices of size $n\times m$. $I_{n}$ and $0_{n\times m}$
are respectively the identity matrix of $\mathbb{R}^{n\times n}$ and the zero
matrix of $\mathbb{R}^{n\times m}$. The subscripts are omitted when obvious
from the context. $X^{T}$ stands for transpose of $X$ while $M>(\text{resp.
}\geq)~{}0$ denotes positive (semi-) definiteness. For the sake of brevity,
$\llbracket k_{1},k_{2}\rrbracket$ refers to the integers from $k_{1}$ to
$k_{2}$ $z_{\llbracket
k_{1},k_{2}\rrbracket}:=\begin{bmatrix}z_{k_{1}}^{T}&\ldots&z_{k_{2}}^{T}\end{bmatrix}^{T}$.
Bold characters denote decision variables in a design problem or optimization
variables in an optimization problem.

##  2 Convex Set Representation by Quadratic Constraints Intersection

This section introduces and motivates the choice of Quadratic Constraints (QC)
as a computational framework. First, QC can be used to represent a large class
of convex sets. In addition, this framework allows for considering the
intersection of these sets, without explicitly computing it. Finally, its link
with LMI optimization provides efficient numerical methods.

######  Definition 1.

Consider a quadratic function $\sigma_{S}(\cdot)$ defined by:

$\forall z\in\mathbb{R}^{n_{z}},\qquad\sigma_{S}(z):=\begin{bmatrix}z\\\
1\end{bmatrix}^{T}S\begin{bmatrix}z\\\ 1\end{bmatrix}$

where $S=S^{T}\in\mathbb{R}^{(n_{z}+1)\times(n_{z}+1)}$. A vector
$z\in\mathbb{R}^{n_{z}}$ is said to satisfy a QC associated with $\sigma_{S}$
if $\sigma_{S}(z)\geq 0$ holds. In addition, $z\in\mathbb{R}^{n_{z}}$ is said
to satisfy a QC equality (QCE) associated with $\sigma_{S}$ if
$\sigma_{S}(z)=0$ holds.

It can be inferred from Definition 1 that the most common classes of convex
sets classically used [23] can be exactly expressed in the QC framework. For
example, an ellipsoid and a degenerated ellipsoid, such as an halfspace, is
inherently represented as a single QC [24]. The set $\mathcal{E}$ defined by

$\mathcal{E}:=\left\\{z\in\mathbb{R}^{n_{z}}\left|\sigma_{P_{\mathcal{E}}}(z)\geq
0\right.\right\\},\quad P_{{\mathcal{E}}}:=\begin{bmatrix}Q&s\\\
s^{T}&r\end{bmatrix}$

leads to an ellipsoid when $Q<0$, a degenerated ellipsoid when $Q\leq 0$ and a
halfspace when $Q=0$. The intersection of QCs offers then the possibility to
represent the intersection of ellipsoids and degenerated ellipsoids. In
particular, a polytope in $H$-representation
$\mathcal{H}:=\left\\{z\in\mathbb{R}^{n_{z}}\left|\ Hz\leq h\right.\right\\}$,
$h\in\mathbb{R}^{n_{h}}$, is represented by the intersection of halfspaces:

$\mathcal{H}=\bigcap_{i=1}^{n_{h}}{\left\\{z\in\mathbb{R}^{n_{z}}\left|\
\sigma_{P_{\mathcal{H}_{i}}}(z)\geq 0\right.\right\\}},\quad
P_{\mathcal{H}_{i}}:=\begin{bmatrix}0&H_{i}^{T}\\\ H_{i}&2h_{i}\end{bmatrix}$

where $H_{i}$ the $i$th line of $H$ and $h_{i}$ the $i$th element of $h$.
Moreover, a QCE especially enables to represent convex sets defined by a
linear equality, such as a hyperplane
$\mathcal{HP}:=\left\\{z\in\mathbb{R}^{n_{z}}\left|\ a^{T}z=b\right.\right\\}$
where $a\in\mathbb{R}^{n_{z}}$ and ${b\in\mathbb{R}}$. Planes, that are the
intersection of hyperplanes, can then be represented by the intersection of
QCEs. Finally, considering the intersection of both QCs and QCEs significantly
extends the scope. For instance, a zonotope
$\mathcal{Z}=\left\\{z\in\mathbb{R}^{n_{z}}\left|z=c+G\lambda,\
\lambda\in[-1,1]^{n_{\lambda}}\right.\right\\}$ can be expressed as the
intersection of $n_{z}$ hyperplanes and $n_{\lambda}$ degenerated ellipsoids,
i.e. QCEs and QCs. Zonotope bundles, that are sets defined as the intersection
of zonotopes, are also exactly represented, as well as constrained zonotopes
${\mathcal{CZ}=\left\\{z\in\mathbb{R}^{n_{z}}\left|z=c+G\lambda,\
D\lambda+d=0,\ \lambda\in[-1,1]^{n_{\lambda}}\right.\right\\}}$ which result
from the intersection of a zonotope with hyperplanes. Up to the numerical
reliability of optimization solvers, QCs provide a useful design framework to
exactly define and operate over a wide class of convex sets.

An important result of the QC framework is the S-procedure lemma that provides
a sufficient condition, in the form of a feasibility optimization problem
under LMI constraints, to test the satisfaction of a QC by a set of vectors
defined by the intersection of QCs [21, Chap. 2]. The next lemma provides an
adapted version when a QC is checked over multiple QCs and QCEs (see [25,
Theorem 2.3.3]).

######  Lemma 1 (S-procedure for intersection of QCs and QCEs).

Consider $(N+1)$ QCs associated with $\sigma_{S}$, $\sigma_{P_{1}}$, $\ldots$,
$\sigma_{P_{N}}$, and $M$ QCEs associated with $\sigma_{Q_{1}}$, $\ldots$,
$\sigma_{Q_{M}}$.
Then $(ii)\Rightarrow(i)$.

  * $(i)$

The QC $\sigma_{S}(z)\geq 0$ holds for all $z\in\mathbb{R}^{n_{z}}$ such that

$\displaystyle\forall p\in\llbracket 1,N\rrbracket$
$\displaystyle\sigma_{P_{p}}(z)\geq 0$ $\displaystyle\forall q\in\llbracket
1,M\rrbracket$ $\displaystyle\sigma_{Q_{q}}(z)=0$

  * $(ii)$

$\exists\left\\{\boldsymbol{\tau_{p}}\in\mathbb{R}^{+}\right\\}_{p=1}^{N}$,
$\exists\left\\{\boldsymbol{\mu_{q}}\in\mathbb{R}\right\\}_{q=1}^{M}$, such
that

$S-\sum_{p=1}^{N}{\boldsymbol{\tau_{p}}P_{p}}-\sum_{q=1}^{M}{\boldsymbol{\mu_{q}}Q_{q}}\geq
0$ (1)

######  Remark 1.

The scalars $\boldsymbol{\tau_{p}}$ and $\boldsymbol{\mu_{q}}$, classically
called multipliers, are the optimization variables in the constraint (1),
leading then to a feasibility problem over LMI constraints.

##  3 Problem Setting

###  3.1 Nominal system description

Consider an uncertain linear discrete-time system $(\Sigma)$:

$(\Sigma)\left\\{\begin{aligned}
x_{k+1}&=A_{\Sigma}x_{k}+B_{\Sigma}u_{k}^{a}+W_{\Sigma}w_{k},\quad
u_{k}^{a}\in\mathcal{U}\\\
y_{k}&=C_{\Sigma}x_{k}+V_{\Sigma}v_{k}\end{aligned}\right.$

where $x_{k}\in\mathbb{R}^{n}$ is the system state evaluated at time step $k$,
$u_{k}^{a}\in\mathbb{R}^{m}$ the applied input and $y_{k}$ the measured
output. The applied input $u_{k}^{a}$ is constrained to belong to the input-
limitation set $\mathcal{U}\subseteq\mathbb{R}^{m}$. The vectors
$w_{k}\in\mathbb{R}^{n_{w}}$ and $v_{k}\in\mathbb{R}^{n_{v}}$ represent
additive state disturbances and measurement noises.

The system $(\Sigma)$ is interconnected, through a communication network, to a
controller $(\mathcal{K})$ modeled as:

$(\mathcal{K})\left\\{\begin{aligned} \xi_{k+1}&=A_{K}\xi_{k}+B_{K}y_{k}\\\
u_{k}&=C_{K}\xi_{k}+D_{K}y_{k}\end{aligned}\right.$

where $\xi_{k}\in\mathbb{R}^{l}$ is the controller state and
$u_{k}\in\mathbb{R}^{m}$ is the controller output. Under normal conditions,
$u_{k}^{a}=u_{k}$.

###  3.2 Sharp input anomalies

In this paper, sharp anomalies, producing severe impact in a reduced amount of
time, on inputs are considered. This especially includes anomalies such as
sudden actuator failures, DoS and attack by saturation. From a security point
of view, sharp attacks result typically from adversaries with little knowledge
of the overall system, and may be accomplished from different manners
(fragility of a communication network protocol, communication channel jamming,
direct attack on the controller, or even resulting from false-data
measurements) [26]. Unlike the more sophisticated class of stealthy attacks,
sharp anomalies may be detected by a passive detection scheme, for instance
based on an observer or using a $\chi^{2}$-detector [2, 4].

The anomalies impact on the applied input is modeled as

$u_{k}^{a}=\Gamma_{u}u_{k}+\Gamma_{a}a_{k}$

where $a_{k}\in\mathcal{A}$ is an unknown signal and
$\mathcal{A}\subseteq\mathbb{R}^{m_{a}}$ the set of possible values for
$a_{k}$, while $\Gamma_{u}\in\mathbb{R}^{m\times m}$ and
$\Gamma_{a}\in\mathbb{R}^{m\times m_{a}}$ denote the input channels on which
the anomalies appear, and are typically binary-valued $\\{0,1\\}$-matrices. A
large class of sharp input anomalies can be modeled by the triplet
$(\Gamma_{u},\Gamma_{a},\mathcal{A})$. For instances,
$(\Gamma_{u},\Gamma_{a},\mathcal{A})=(0,I,\mathcal{U})$ models the worst-case
scenario where an adversary may apply any admissible input,
$(\Gamma_{u},\Gamma_{a},\mathcal{A})=(0,I,0)$ may be associated with a DoS,
$(\Gamma_{u},\Gamma_{a})=(I,I)$ models additive input faults or data-injection
attacks, that are with saturation objective when $\mathcal{A}=\\{a=\alpha\cdot
u^{a}_{max}\\}$ where $u^{a}_{max}$ the maximum admissible input and
$\alpha\in\mathbb{R}$ a large number, while

$\Gamma_{u}=\begin{bmatrix}I&0&0\\\ 0&0_{1\times 1}&0\\\
0&0&I\end{bmatrix}\quad\text{and}\quad\Gamma_{a}=\begin{bmatrix}0&0&0\\\
0&1&0\\\ 0&0&0\end{bmatrix}$

model anomalies that appear on the $i^{th}$ input channel. Normal operating
conditions are obtained when $(\Gamma_{u},\Gamma_{a})=(I,0)$.

###  3.3 Controlled system subject to sharp input anomalies

Figure 1: Controlled system subject to sharp input anomalies

The dynamics of the interconnected system $(\Sigma-\mathcal{K})$ subject to
input anomalies (Fig. 1) is then given by:

$(\Sigma-\mathcal{K})\left\\{\begin{aligned}
\overline{x}_{k+1}&=A\overline{x}_{k}+W\overline{w}_{k}+B_{a}a_{k}\\\
u^{a}_{k}&=C_{u}\overline{x}_{k}+V_{u}\overline{w}_{k}+\Gamma_{a}a_{k},\quad
u^{a}_{k}\in\mathcal{U}\end{aligned}\right.$ (2)

where
$\overline{x}_{k}:=\begin{bmatrix}x_{k}^{T}&\xi_{k}^{T}\end{bmatrix}^{T}$,
$\overline{w}_{k}:=\begin{bmatrix}w_{k}^{T}&v_{k}^{T}\end{bmatrix}^{T}$ and

$A=\begin{bmatrix}A_{\Sigma}+B_{\Sigma}\Gamma_{u}D_{K}C_{\Sigma}&B_{\Sigma}\Gamma_{u}C_{K}\\\
B_{K}C_{\Sigma}&A_{K}\end{bmatrix},\
B_{a}=\begin{bmatrix}B_{\Sigma}\Gamma_{a}\\\ 0\end{bmatrix}$
$W=\begin{bmatrix}W_{\Sigma}&B_{\Sigma}\Gamma_{u}D_{K}V_{\Sigma}\\\
0&B_{K}V_{\Sigma}\end{bmatrix},\quad\begin{aligned}
C_{u}&=\Gamma_{u}\begin{bmatrix}D_{K}C_{\Sigma}&C_{K}\end{bmatrix}\\\
V_{u}&=\Gamma_{u}\begin{bmatrix}0&D_{K}V_{\Sigma}\end{bmatrix}\end{aligned}$

This interconnected system $(\Sigma-\mathcal{K})$ is provided with a safety
set $\overline{\mathcal{S}}\subseteq\mathbb{R}^{n+l}$, i.e. a set of states
$\overline{x}_{k}$ for which $(\Sigma-\mathcal{K})$ is considered to be safe
at each time step $k$. For the ease of computation, and without loss of
generality, it is assumed that the anomalies start from $k_{a}=0$. The initial
state $\overline{x}_{0}$ is assumed to belong to a known and safe set
$\overline{\mathcal{X}}_{0}\subseteq\overline{\mathcal{S}}$, and the additive
uncertainties $\overline{w}_{k}$ to a bounded set
$\overline{\mathcal{W}}\subseteq\mathbb{R}^{n_{w}+n_{v}}$.

The defense objective is to design off-line a strategy such that safe
operation is preserved after the occurrence of an anomaly. In this
perspective, an (under)-estimation of the critical-time, that is the maximum
time window allowed before the system leaves the safety set, appears as a
crucial information to be provided to the defense designer.

###  3.4 Problem formulation

The problem considered in this paper is the computation of the critical-time
for the interconnected system (2) subject to sharp input anomalies defined by
$(\Gamma_{u},\Gamma_{a},\mathcal{A})$, based on the QC framework. To this end,
the safety set $\overline{\mathcal{S}}$ is assumed to be exactly represented
or inner-approximated by the intersection of QCs (Definition 1):

$\bigcap_{s=1}^{n_{\overline{\mathcal{S}}}}{\overline{\mathcal{S}}_{s}}\subseteq\overline{\mathcal{S}},\qquad\overline{\mathcal{S}}_{s}:=\left\\{\overline{x}\in\mathbb{R}^{n+l}\left|\sigma_{\overline{S}_{s}}(\overline{x})\geq
0\right.\right\\}$ (3)

and the initial set $\overline{\mathcal{X}}_{0}$, the input-limitation set
$\mathcal{U}$, the deviation set $\overline{\mathcal{W}}$ and the anomaly-
impact set $\mathcal{A}$ are exactly represented or over-approximated by the
intersection of QCs and QCEs:

$\overline{\mathcal{X}}_{0}\subseteq\left\\{\cap_{g=1}^{n_{\overline{\mathcal{X}}_{0}}}{\
\overline{\mathcal{X}}_{g}}\right\\}\bigcap\left\\{\cap_{h=1}^{m_{\overline{\mathcal{X}}_{0}}}{\
\overline{\mathcal{X}}_{h}^{e}}\right\\}$ (4)
$\mathcal{U}\subseteq\left\\{\cap_{i=1}^{n_{\mathcal{U}}}{\
\mathcal{U}_{i}}\right\\}\bigcap\left\\{\cap_{j=1}^{m_{\mathcal{U}}}{\
\mathcal{U}_{j}^{e}}\right\\}$ (5)
$\overline{\mathcal{W}}\subseteq\left\\{\cap_{p=1}^{n_{\overline{\mathcal{W}}}}{\
\overline{\mathcal{W}}_{p}}\right\\}\bigcap\left\\{\cap_{q=1}^{m_{\overline{\mathcal{W}}}}{\
\overline{\mathcal{W}}_{q}^{e}}\right\\}$ (6)
$\mathcal{A}\subseteq\left\\{\cap_{r=1}^{n_{\mathcal{A}}}{\
\mathcal{A}_{r}}\right\\}\bigcap\left\\{\cap_{t=1}^{m_{\mathcal{A}}}{\
\mathcal{A}_{t}^{e}}\right\\}$ (7)

with
$\overline{\mathcal{X}}_{g}:=\\{\overline{x}_{0}~{}\in~{}\mathbb{R}^{n+l}\ |\
\sigma_{\overline{X}_{g}}(\overline{x}_{0})~{}\geq~{}0\\}$,
$\overline{\mathcal{X}}^{e}_{h}:=\\{\overline{x}_{0}~{}\in~{}\mathbb{R}^{n+l}\
|\ \sigma_{\overline{X}^{e}_{h}}(x_{0})~{}=~{}0\\}$,
$\mathcal{U}_{i}:=\\{u^{a}~{}\in~{}\mathbb{R}^{m}\ |\
\sigma_{U_{i}}(u^{a})~{}\geq~{}0\\}$,
$\mathcal{U}^{e}_{j}:=\\{u^{a}~{}\in~{}\mathbb{R}^{m}\ |\
\sigma_{U^{e}_{j}}(u^{a})~{}=~{}0\\}$,
$\overline{\mathcal{W}}_{p}:=\\{\overline{w}~{}\in~{}\mathbb{R}^{n_{v}+n_{w}}\
|\ \sigma_{\overline{W}_{p}}(\overline{w})~{}\geq~{}0\\}$,
$\overline{\mathcal{W}}^{e}_{q}:=\\{\overline{w}~{}\in~{}\mathbb{R}^{n_{v}+n_{w}}\
|\ \sigma_{\overline{W}^{e}_{q}}(w)~{}=~{}0\\}$,
$\mathcal{A}_{r}:=\\{a~{}\in~{}\mathbb{R}^{m_{a}}\ |\
\sigma_{A_{r}}(a)~{}\geq~{}0\\}$, and
${\mathcal{A}^{e}_{t}:=\\{a~{}\in~{}\mathbb{R}^{m_{a}}\ |\
\sigma_{A^{e}_{t}}(a)~{}=~{}0\\}}$.

######  Remark 2.

The introduction of QCEs makes possible to exactly represent anomalies impact
leading to constant $a_{k}$, such as in the cases of DoS or attack by upper
saturation.

The problem can now be stated as follows.

######  Problem 1.

Given the system (2) subject to input anomalies modeled by
$(\Gamma_{u},\Gamma_{a},\mathcal{A})$, Given a safety set
$\overline{\mathcal{S}}$, an initial set
$\overline{\mathcal{X}}_{0}\subseteq\overline{\mathcal{S}}$, an input-
limitation set $\mathcal{U}$ and a deviation set $\overline{W}$, Given the
approximation (3)-(7) in the QC framework,
Find the maximal time-index $\boldsymbol{k_{f}}$ such that system (2) is safe
over the time-window $\llbracket 0,\boldsymbol{k_{f}}\rrbracket$ for all
initial states of $\overline{\mathcal{X}}_{0}$, for all anomalies of
$\mathcal{A}$ resulting in admissible inputs of $\mathcal{U}$, and for all
elements of $\overline{W}$:

$\max_{\boldsymbol{k_{f}}\in\mathbb{N}}{\boldsymbol{k_{f}}}$

subject to $\forall\overline{x}_{0}\in\overline{\mathcal{X}}_{0}$, $\forall
k\in\llbracket 0,\boldsymbol{k_{f}}\rrbracket$, $\forall a_{k}\in\mathcal{A}$
such that $u_{k}^{a}\in\mathcal{U}$, $\forall\overline{w}_{k}\in\overline{W}$,
$\overline{x}_{k}\in\overline{S}$.

##  4 Critical-time computation

In this section, first the following problem is solved (Theorem 1): how to
check if the state $\overline{x}_{k_{f}}$ belongs to the safety set
$\overline{\mathcal{S}}$, over all initial states of
$\overline{\mathcal{X}}_{0}$, occurred anomalies in $\mathcal{A}$ resulting in
admissible inputs in $\mathcal{U}$, and uncertainties in
$\overline{\mathcal{W}}$, for a fixed time-index $k_{f}$. Then, an iterative
algorithm (Algorithm 1) is proposed to find the maximal $\boldsymbol{k_{f}}$,
providing an underestimate of the critical-time.

######  Theorem 1.

Let $(\Sigma-\mathcal{K})$ be the interconnected system modeled by (2). Assume
that the safety set $\overline{\mathcal{S}}$, the initial-state set
$\overline{\mathcal{X}}_{0}$, the input-limitation set $\mathcal{U}$, the
deviation set $\overline{\mathcal{W}}$ and the anomaly set $\mathcal{A}$ are
provided with QC-based approximations given by (3)-(7). Assume that the time-
index $k_{f}\in\mathbb{N}$ is fixed.
Then $(i)\Rightarrow(ii)$.

  * $(i)$

$\forall s\in\llbracket 1,n_{\overline{\mathcal{S}}}\rrbracket$,
$\exists\left\\{\boldsymbol{\phi_{s,g}}\in\mathbb{R}^{+}\right\\}_{g=1}^{n_{\overline{\mathcal{X}}_{0}}}$,
$\exists\left\\{\boldsymbol{\varphi_{s,h}}\in\mathbb{R}\right\\}_{h=1}^{m_{\overline{\mathcal{X}}_{0}}}$,
$\forall k\in\llbracket 0,k_{f}-1\rrbracket$,
$\exists\left\\{\boldsymbol{\pi_{s,k,i}}\in\mathbb{R}^{+}\right\\}_{i=1}^{n_{\mathcal{U}}}$,
$\exists\left\\{\boldsymbol{\varpi_{s,k,j}}\in\mathbb{R}\right\\}_{j=1}^{m_{\mathcal{U}}}$,
$\exists\left\\{\boldsymbol{\rho_{s,k,p}}\in\mathbb{R}^{+}\right\\}_{p=1}^{n_{\overline{\mathcal{W}}}}$,
$\exists\left\\{\boldsymbol{\varrho_{s,k,q}}\in\mathbb{R}\right\\}_{q=1}^{m_{\overline{\mathcal{W}}}}$,
$\exists\left\\{\boldsymbol{\sigma_{s,k,r}}\in\mathbb{R}^{+}\right\\}_{r=1}^{n_{\mathcal{A}}}$,
$\exists\left\\{\boldsymbol{\varsigma_{s,k,t}}\in\mathbb{R}\right\\}_{t=1}^{m_{\mathcal{A}}}$,
such that:

$\displaystyle
M_{\overline{x}_{k_{f}}}^{T}\overline{S}_{s}M_{\overline{x}_{k_{f}}}-M_{\overline{x}_{0}}^{T}\left(\sum_{g=1}^{n_{\overline{\mathcal{X}}_{0}}}{\boldsymbol{\phi_{s,g}}\overline{X}_{g}}+\sum_{h=1}^{m_{\overline{\mathcal{X}}_{0}}}{\boldsymbol{\varphi_{s,h}}\overline{X}_{h}^{e}}\right)M_{\overline{x}_{0}}$
$\displaystyle-$
$\displaystyle\sum_{k=0}^{k_{f}-1}{N_{u_{k}^{a}}^{T}\left(\sum_{i=1}^{n_{\mathcal{U}}}{\boldsymbol{\pi_{s,k,i}}U_{i}}+\sum_{j=1}^{m_{\mathcal{U}}}{\boldsymbol{\varpi_{s,k,j}}U_{j}^{e}}\right)N_{u_{k}^{a}}}$
$\displaystyle-$
$\displaystyle\sum_{k=0}^{k_{f}-1}{E_{k}^{T}\left(\sum_{p=1}^{n_{\overline{\mathcal{W}}}}{\boldsymbol{\rho_{s,k,p}}\overline{W}_{p}}+\sum_{q=1}^{m_{\overline{\mathcal{W}}}}{\boldsymbol{\varrho_{s,k,q}}\overline{W}_{q}^{e}}\right)E_{k}}$
$\displaystyle-$
$\displaystyle\sum_{k=0}^{k_{f}-1}{F_{k}^{T}\left(\sum_{r=1}^{n_{\mathcal{A}}}{\boldsymbol{\sigma_{s,k,r}}A_{r}}+\sum_{t=1}^{m_{\mathcal{A}}}{\boldsymbol{\varsigma_{s,k,t}}A_{t}^{e}}\right)F_{k}}\geq
0$ (8)

with

$\displaystyle M_{\overline{x}_{k}}$
$\displaystyle=\begin{bmatrix}A^{k}&P_{k}&Q_{k}&0\\\ 0&0&0&1\end{bmatrix}$
$\displaystyle P_{k}$
$\displaystyle=\begin{bmatrix}A^{k-1}W&\ldots&AW&W&0_{(n+l)\times(k_{f}-k)(n_{v}+n_{w})}\end{bmatrix}$
$\displaystyle Q_{k}$
$\displaystyle=\begin{bmatrix}A^{k-1}B_{a}&\ldots&AB_{a}&B_{a}&0_{(n+l)\times(k_{f}-k)m_{a}}\end{bmatrix}$
$\displaystyle N_{u_{k}^{a}}$
$\displaystyle=C_{u}M_{\overline{x}_{k}}+V_{u}E_{k}+\Gamma_{a}F_{k}$

where $P_{0}=0$ and $Q_{0}=0$, and $E_{k}$, $F_{k}$ are such that:

$E_{k}z_{k_{f}}=\begin{bmatrix}\overline{w}_{k}\\\ 1\end{bmatrix}\qquad
F_{k}z_{k_{f}}=\begin{bmatrix}a_{k}\\\ 1\end{bmatrix}$

where
$z_{k_{f}}^{T}:=\begin{bmatrix}\overline{x}_{0}^{T}&\overline{w}_{\llbracket
0,k_{f}-1\rrbracket}^{T}&a_{\llbracket
0,k_{f}-1\rrbracket}^{T}&1\end{bmatrix}^{T}$.

  * $(ii)$

The system is safe at time-index $k_{f}$: $\overline{x}_{k_{f}}\in\mathcal{S}$
for all initial states $\overline{x}_{0}\in\overline{\mathcal{X}}_{0}$, for
all anomalies $a_{k}\in\mathcal{A}$ resulting in admissible inputs
$u_{k}^{a}\in\mathcal{U}$, and for all deviation elements
$\overline{w}_{k}\in\overline{\mathcal{W}}$.

###### Proof.

Using the recursive equation (2), it comes

$M_{\overline{x}_{k_{f}}}z_{k_{f}}=\begin{bmatrix}\overline{x}_{k_{f}}\\\
1\end{bmatrix}\quad
M_{\overline{x}_{0}}z_{k_{f}}=\begin{bmatrix}\overline{x}_{0}\\\
1\end{bmatrix}\quad N_{u_{k}^{a}}z_{k_{f}}=\begin{bmatrix}u_{k}^{a}\\\
1\end{bmatrix}$

Then, by applying Lemma 1 to (8),
$\sigma_{\overline{S}_{s}}(\overline{x}_{f})\geq 0$ holds for all
$\overline{x}_{0}$, $\left\\{\overline{w}_{k}\right\\}_{k=0}^{k_{f}-1}$, and
$\left\\{u_{k}^{a}\right\\}_{k=0}^{k_{f}-1}$ such that:

$\displaystyle\forall g\in\llbracket
1,n_{\overline{\mathcal{X}}_{0}}\rrbracket,\
\sigma_{\overline{X}_{g}}(\overline{x}_{0})\geq 0\quad\forall h\in\llbracket
1,m_{\overline{\mathcal{X}}_{0}}\rrbracket,\
\sigma_{\overline{X}_{h}^{e}}(\overline{x}_{0})=0$ $\displaystyle\forall
i\in\llbracket 1,n_{\mathcal{U}}\rrbracket,\ \sigma_{U_{i}}(u_{k}^{a})\geq 0\
\quad\forall j\in\llbracket 1,m_{\mathcal{U}}\rrbracket,\
\sigma_{U_{j}^{e}}(u_{k}^{a})=0$ $\displaystyle\forall p\in\llbracket
1,n_{\overline{\mathcal{W}}}\rrbracket,\
\sigma_{\overline{W}_{p}}(\overline{w}_{k})\geq 0\quad\forall q\in\llbracket
1,m_{\overline{\mathcal{W}}}\rrbracket,\
\sigma_{\overline{W}_{q}^{e}}(\overline{w}_{k})=0$ $\displaystyle\forall
r\in\llbracket 1,n_{\mathcal{A}}\rrbracket,\ \sigma_{A_{r}}(a_{k})\geq 0\
\quad\forall t\in\llbracket 1,m_{\mathcal{A}}\rrbracket,\
\sigma_{A_{t}^{e}}(a_{k})=0$

Thus, for all
$\overline{x}_{0}\in\left\\{\cap_{g=1}^{n_{\overline{\mathcal{X}}_{0}}}{\
\overline{\mathcal{X}}_{g}}\right\\}\bigcap\left\\{\cap_{h=1}^{m_{\overline{\mathcal{X}}_{0}}}{\
\overline{\mathcal{X}}_{h}^{e}}\right\\}$, for all
$\left\\{a_{k}\right\\}_{k=0}^{k_{f}-1}$ such that
$a_{k}\in\left\\{\cap_{r=1}^{n_{\mathcal{A}}}{\
\mathcal{A}_{r}}\right\\}\bigcap\left\\{\cap_{t=1}^{m_{\mathcal{A}}}{\
\mathcal{A}_{t}^{e}}\right\\}$ and
$u_{k}^{a}\in\left\\{\cap_{i=1}^{n_{\mathcal{U}}}{\
\mathcal{U}_{i}}\right\\}\bigcap\left\\{\cap_{j=1}^{m_{\mathcal{U}}}{\
\mathcal{U}_{j}^{e}}\right\\}$, and for all
$\left\\{\overline{w}_{k}\right\\}_{k=0}^{k_{f}-1}$ such that
$\overline{w}_{k}\in\left\\{\cap_{p=1}^{n_{\overline{\mathcal{W}}}}{\
\overline{\mathcal{W}}_{p}}\right\\}\bigcap\left\\{\cap_{q=1}^{m_{\overline{\mathcal{W}}}}{\
\overline{\mathcal{W}}_{q}^{e}}\right\\}$, one has
$\overline{x}_{k_{f}}\in\bigcap_{s=1}^{n_{\overline{\mathcal{S}}}}{\overline{\mathcal{S}}_{s}}$.
The inclusions (3)-(7) lead to the conclusion. ∎

######  Remark 3.

Condition $(i)$ is a feasibility optimization problem with LMI constraints. It
is then solvable in polynomial-time. This computational advantage comes at the
price of the conservatism of condition $(i)$, that is only sufficient. This
results from the application of the S-procedure (see [24], and the references
therein, for a detailed discussion).

In order to compute the critical-time $k_{c}$, the following iterative LMI-
based algorithm (Algorithm 1) is proposed: based on the sufficient condition
provided by Theorem 1, the time index $k_{f}$ is incremented until no solution
is found to the LMI feasibility problem. An under-estimate $k_{f}^{*}$ of the
critical-time $k_{c}$ is so obtained $k_{f}^{*}\leq k_{c}$ and ensures that
the system is safe at least until the obtained result $k_{f}^{*}$.

1 $k_{f}\leftarrow 0$

2 $t\leftarrow true$

3 while  _$t$_ do

4 $k_{f}\leftarrow k_{f}+1$

5 $t\leftarrow isFeasible$( $\forall s\in\llbracket
1,n_{\overline{\mathcal{S}}}\rrbracket$,
$\exists\left\\{\boldsymbol{\phi_{s,g}}\in\mathbb{R}^{+}\right\\}_{g=1}^{n_{\overline{\mathcal{X}}_{0}}}$,
$\exists\left\\{\boldsymbol{\varphi_{s,h}}\in\mathbb{R}\right\\}_{h=1}^{m_{\overline{\mathcal{X}}_{0}}}$,
$\forall k\in\llbracket 0,k_{f}-1\rrbracket$,
$\exists\left\\{\boldsymbol{\pi_{s,k,i}}\in\mathbb{R}^{+}\right\\}_{i=1}^{n_{\mathcal{U}}}$,
$\exists\left\\{\boldsymbol{\varpi_{s,k,j}}\in\mathbb{R}\right\\}_{j=1}^{m_{\mathcal{U}}}$,
$\exists\left\\{\boldsymbol{\rho_{s,k,p}}\in\mathbb{R}^{+}\right\\}_{p=1}^{n_{\overline{\mathcal{W}}}}$,
$\exists\left\\{\boldsymbol{\varrho_{s,k,q}}\in\mathbb{R}\right\\}_{q=1}^{m_{\overline{\mathcal{W}}}}$,
$\exists\left\\{\boldsymbol{\sigma_{s,k,r}}\in\mathbb{R}^{+}\right\\}_{r=1}^{n_{\mathcal{A}}}$,
$\exists\left\\{\boldsymbol{\varsigma_{s,k,t}}\in\mathbb{R}\right\\}_{t=1}^{m_{\mathcal{A}}}$,
such that (8) holds. )

6 /* true if a solution is found /*

7$k_{f}\leftarrow k_{f}-1$

return $k_{f}$

Algorithm 1 Critical-time computation

##  5 Case study - the quadruple-tank system

Figure 2: The quadruple-tank system

The quadruple-tank system (Fig. 2) originally described in [22], is considered
to illustrate the benefits of the critical-time. This benchmark has been used
in the CPS security literature for simulating cyber-attacks, evaluating their
physical impacts and developing counter-measures (see [6, 27] for instance).
It consists of four tanks with bottom outlet holes, arranged such that tank 3
is above tank 1 and tank 4 above tank 2. The inputs are the voltages
$v_{1_{in}}$, $v_{2_{in}}$ applied to two pumps, and the outputs are the
voltages $v_{1_{out}}$, $v_{2_{out}}$ from level measurement devices of tanks
1 and 2. The input liquid from pump 1 (resp. pump 2) is split into tank $1$
and $4$ (tank $2$ and $3$) using a valve with coefficient
$\gamma_{1}~{}\in~{}(0,1)$ (resp. $\gamma_{2}\in(0,1)$). The system is modeled
as follows:

$\left\\{\begin{aligned}
\dot{h}(t)&=A_{P}\sqrt{h(t)}+B_{P}(\gamma_{1},\gamma_{2})v_{in}(t)\\\
v_{out}(t)&=C_{P}h(t)+d(t)\end{aligned}\right.$ (9)

where $h=\begin{bmatrix}h_{1}&h_{2}&h_{3}&h_{4}\end{bmatrix}^{T}$, $h_{i}$ the
liquid level of tank $i$,
$v_{in}=\begin{bmatrix}v_{1_{in}}&v_{2_{in}}\end{bmatrix}^{T}$,
$v_{out}=\begin{bmatrix}v_{1_{out}}&v_{2_{out}}\end{bmatrix}^{T}$ and $d$ the
sensor noise. The matrices $A_{P}$ and $C_{P}$ depend on constant physical
parameters of the system, while $B_{P}(\gamma_{1},\gamma_{2})$ is an affine
matrix-valued function. See [22] for numerical values.

Due to physical limitations of the pumps, the inputs $v_{1_{in}}$ and
$v_{2_{in}}$ belong to $[v_{min},v_{max}]=[0\text{ V},10\text{ V}]$. Moreover,
for safety reasons, the levels $h_{1}$ and $h_{2}$ are constrained to belong
to $[h_{12_{min}},h_{12_{max}}]=[9\text{ cm},20\text{ cm}]$ and the levels
$h_{3}$ and $h_{4}$ to $[h_{34_{min}},h_{34_{max}}]=[0.1\text{ cm},5\text{
cm}]$. In particular, this prevents the tanks to overflow or to be empty, and
keeps a minimum output flow for tanks $1$ and $2$.

The system is remotely controlled through a communication network by PI
controllers, that compute the control input $v_{c}$ from a reference signal
$r$ and measurement $v_{out}$ such as $v_{i_{c}}=K_{i}(1+\frac{1}{T_{i}\cdot
s})(r_{i}-v_{i_{out}})$ where $K_{1}=3$, $K_{2}=2.7$, $T_{1}=30s$,
$T_{2}=40s$. The control objective is to set $h_{1}$ and $h_{2}$ to
$h_{1}^{0}=12.4\text{ cm}$ and $h_{2}^{0}=12.7\text{ cm}$. Anomalies are
assumed to appear when the process is at a stationary operating point
$(h^{0},v_{in}^{0})$. It is known [22] that, given $(h_{1}^{0},h_{2}^{0})$,
the values of $h_{3}^{0}$, $h_{4}^{0}$ and $v_{in}^{0}$ are fully determined
by the valves coefficients $(\gamma_{1},\gamma_{2})$. The choice of
$(\gamma_{1},\gamma_{2})$ impacts then both the operating point and the system
dynamics. Thus, the aim is to study the influence of $(\gamma_{1},\gamma_{2})$
on the critical-time. In addition to evaluate the time allowed to defense
mechanism given a couple $(\gamma_{1},\gamma_{2})$, this provides a means to
select couples with high critical-time. To this end, the critical-time is
computed for a finite number of arbitrarily-chosen $\gamma_{1}\in[0.65,0.75]$,
$\gamma_{2}\in[0.55,0.65]$.

###  5.1 Uncertain linear discrete-time approximation

For each couple $(\gamma_{1},\gamma_{2})$, the corresponding model (9) is
linearized around the associated operating point $(h^{0},v_{in}^{0})$. Then,
using the zero-order hold method with a sampling period of $1\text{s}$, the
model (2) is obtained, where $\overline{w}_{k}$ represents the deviation
between (9) and the approximation, especially resulting from linearization-
and-discretization errors and sensor noise. Based on several simulations, the
set $\overline{\mathcal{W}}$ is over-approximated (see Fig. 3) by the
intersection of two degenerated ellipsoids $\overline{\mathcal{W}}_{1}$ and
$\overline{\mathcal{W}}_{2}$,

$\overline{W}_{1}=\begin{bmatrix}Q_{\overline{W}_{1}}&0&s_{\overline{W}_{1}}\\\
0&0&0\\\
s_{\overline{W}_{1}}^{T}&0&r_{\overline{W}_{1}}\end{bmatrix}\qquad\overline{W}_{2}=\begin{bmatrix}Q_{\overline{W}_{2}}&0&s_{\overline{W}_{2}}\\\
0&0&0\\\ s_{\overline{W}_{2}}^{T}&0&r_{\overline{W}_{2}}\end{bmatrix}$

where

$Q_{\overline{W}_{1}}=\begin{bmatrix}-0.567&0&-0.489&0\\\ 0&0&0&0\\\
-0.489&0&-0.448&0\\\ 0&0&0&0\end{bmatrix}\quad
s_{\overline{W}_{1}}=\begin{bmatrix}0.000263\\\ 0\\\ 0.000433\\\
0\end{bmatrix}\quad r_{\overline{W}_{1}}=-1.20\cdot 10^{-8}$
$Q_{\overline{W}_{2}}=\begin{bmatrix}0&0&0&0\\\ 0&-0.611&0&-0.480\\\
0&0&0&0\\\ 0&-0.480&0&-0.409\\\ \end{bmatrix}\quad
s_{\overline{W}_{2}}=\begin{bmatrix}0\\\ 0.000569\\\ 0\\\
0.000501\end{bmatrix}\quad r_{\overline{W}_{2}}=-3.70\cdot 10^{-7}$

and a hyper-rectangle

$\left\\{\overline{w}\in\mathbb{R}^{6},\
\overline{w}_{min}\leq\overline{w}\leq\overline{w}_{max}\right\\}$

with

$\displaystyle\overline{w}_{min}$
$\displaystyle=-\begin{bmatrix}0.013&0.0027&0&0&0.137&0.126\end{bmatrix}^{T}$
$\displaystyle\overline{w}_{max}$
$\displaystyle=\begin{bmatrix}0.00011&0.0013&0.014&0.0045&0.143&0.147\end{bmatrix}^{T}$

Similarly, the initial set $\overline{\mathcal{X}}_{0}$ is over-approximated
by the hyper-rectangle

$\left\\{\overline{x}_{0}\in\mathbb{R}^{6},\
\overline{x}_{0_{min}}\leq\overline{x}_{0}\leq\overline{x}_{0_{max}}\right\\}$

with

$\displaystyle\overline{x}_{0_{min}}$
$\displaystyle=-\begin{bmatrix}0.065&0.047&0.032&0.028&0.066&0.061\end{bmatrix}^{T}$
$\displaystyle\overline{x}_{0_{max}}$
$\displaystyle=\begin{bmatrix}0.066&0.056&0.044&0.024&0.068&0.077\end{bmatrix}^{T}$

Figure 3: Projection of the hyper-rectangle (cyan solid line) and a
degenerated ellipsoid (red solid line) in the plane
$(\overline{w}_{k}(1),\overline{w}_{k}(3))$. Blue crosses denote samples
coming from simulations. The intersection of the rectangle and the ellipse
includes the whole set of crosses.

###  5.2 Critical-time computation and usage

Algorithm 1 is now applied to estimate the critical-time in three scenarios.
The computation is performed using matlab R2022b, using the semidefinite mode
of the cvx package [28] with the solver mosek [29], on a computer with 4-core
3.00 GHz CPU and 32GB RAM.

####  5.2.1 Anomaly-specific scenario

A DoS anomaly, resulting for instance from an unreliable communication network
or a jamming cyber-attack, is first considered. When an anomaly is especially
risky, an anomaly-specific defense strategy may be designed in the treatment
stage. As a case in point, a DoS may be detected using acknowledgment-based
communication protocol such as TCP [30]. This anomaly is modeled as
$(\Gamma_{u},\Gamma_{a},\mathcal{A})=(0_{2},I_{2},\mathcal{U}_{\gamma})$ where

$\mathcal{U}_{\gamma}=\left\\{u\in\mathbb{R}^{2},\
u=\begin{bmatrix}u_{\gamma_{1}}\\\ u_{\gamma_{2}}\end{bmatrix},\
u_{\gamma_{i}}=0-v_{0_{i}}\right\\}$

recalling that the values of $v_{0}$ depends on $(\gamma_{1},\gamma_{2})$. The
result is displayed in Table 1. Similarly, an attack-by-upper saturation is
considered and the estimated critical-time is displayed in Table 2. Both
results are compared with the result obtained from simulation on the non-
linear model (9). It can be especially observed that results from Algorithm 1
and from simulation are close, suggesting a moderate conservatism of the
computing approach developed in this paper.

Table 1: Critical-time (s) - DoS (diff. simulation vs algorithm) $\gamma_{1}$ $\gamma_{2}$ | $0.55$ | $0.575$ | $0.60$ | $0.625$ | $0.65$
---|---|---|---|---|---
$0.65$ | $13$ | $12$ |  $10$ $(+1)$ |  $8$ $(+1)$ |  $6$ $(+1)$
$0.675$ |  $13$ $(+1)$ |  $12$ $(+1)$ |  $11$ $(+1)$ |  $10$ $(+1)$ |  $8$ $(+1)$
$0.70$ |  $14$ $(+1)$ | $13$ | $12$ |  $11$ $(+1)$ |  $9$ $(+2)$
$0.725$ |  $10$ $(+2)$ |  $13$ $(+1)$ |  $12$ $(+1)$ | $12$ |  $10$ $(+1)$
$0.75$ |  $7$ $(+2)$ |  $10$ $(+1)$ |  $11$ $(+2)$ | $12$ |  $11$ $(+1)$
Table 2: Critical-time (s) - Upper Saturation (diff. simulation vs algorithm) $\gamma_{1}$ $\gamma_{2}$ | $0.55$ | $0.575$ | $0.60$ | $0.625$ | $0.65$
---|---|---|---|---|---
$0.65$ | $9$ | $11$ | $10$ | $8$ | $7$
$0.675$ | $7$ | $10$ | $13$ |  $13$ $(+1)$ | $13$
$0.70$ | $6$ | $9$ | $12$ |  $12$ $(+1)$ | $13$
$0.725$ | $5$ | $8$ | $11$ | $12$ | $12$
$0.75$ | $4$ | $7$ |  $10$ $(+1)$ | $11$ |  $11$ $(+1)$

####  5.2.2 Worst-case scenario

Figure 4: Impact on $h_{3}$ and $h_{4}$ of a DoS (dashed lines) and an attack
by upper saturation (dotted lines) on the quadruple tank system originally at
an operating point (solid lines). Grey patches refer to unsafe zones.

In this scenario, the class of anomalies leading to any admissible input
$(\Gamma_{u},\Gamma_{a},\mathcal{A})=(0,I,\mathcal{U})$ is considered. This
worst-case scenario may typically happen when an adversary takes full control
of the inputs, either through the communication network or by an attack on the
controller. In that case, a difficulty is to design a defense strategy for
anomalies with opposite impact (see for instance Fig. 4). Indeed, one may
observe for instance in Table 1 and Table 2 that the couple
$(\gamma_{1},\gamma_{2})=(0.70,0.55)$ leads to the highest critical-time
($14~{}\text{s}$) for a DoS anomaly but low critical-time ($6~{}\text{s}$)
when considering an attack-by-upper saturation. While an intuitive solution
may be found on this example and these two particular anomalies, a suitable
trade-off is expected to be more challenging to be made for numerous and
diverse anomalies. On the other hand, the computation of the critical-time for
all admissible anomalies naturally incorporates this tradeoff, without needing
to simulate each element of $\mathcal{U}$.

The computed under-estimation of the critical-time in this scenario is
displayed in Table 3. In addition, the computational time needed by Algorithm
1 at each iteration is plotted on Fig. 5. As expected, the computational time
increases moderately with the time horizon $k$ and remains tractalbe.. The
execution time of Algorithm 1 is $7.3~{}s$ on average.

Table 3: Critical-time (s) - Worst-case scenario $\gamma_{1}$ $\gamma_{2}$ | $0.55$ | $0.575$ | $0.60$ | $0.625$ | $0.65$
---|---|---|---|---|---
$0.65$ | $9$ | $11$ | $10$ | $8$ | $6$
$0.675$ | $7$ | $10$ | $11$ | $10$ | $8$
$0.70$ | $6$ | $9$ | $12$ | $11$ | $9$
$0.725$ | $5$ | $8$ | $11$ | $12$ | $10$
$0.75$ | $4$ | $7$ | $10$ | $11$ | $11$
Figure 5: Computational time per iteration in the worst-case scenario

####  5.2.3 Channel-dependent scenario

The classes of anomalies leading to any admissible input for a single input
channel is now considered. These are modeled by

$(\Gamma_{u_{1}},\Gamma_{a_{1}},\mathcal{A}_{1})=\left(\begin{bmatrix}0&0\\\
0&1\end{bmatrix},\begin{bmatrix}1\\\ 0\end{bmatrix},\mathcal{U}_{1}\right)$

for the first channel $C_{1}$ and by

$(\Gamma_{u_{2}},\Gamma_{a_{2}},\mathcal{A}_{2})=\left(\begin{bmatrix}1&0\\\
0&0\end{bmatrix},\begin{bmatrix}0\\\ 1\end{bmatrix},\mathcal{U}_{2}\right)$

for the second channel $C_{2}$, where

$\mathcal{U}_{1}:=\left\\{u\in\mathbb{R},\ \begin{bmatrix}u\\\
0\end{bmatrix}\in\mathcal{U}\right\\}\quad\mathcal{U}_{2}:=\left\\{u\in\mathbb{R},\
\begin{bmatrix}0\\\ u\end{bmatrix}\in\mathcal{U}\right\\}$

The estimated critical-time by Algorithm 1 is displayed in Fig. 4.

Table 4: Critical-time (s) - Channel-dependent scenario C1 $|$ C2 $\gamma_{1}$ $\gamma_{2}$ | $0.55$ | $0.575$ | $0.60$ | $0.625$ | $0.65$
---|---|---|---|---|---
$0.65$ |  $14$ $|$ $9$ |  $11$ $|$ $12$ |  $10$ $|$ $10$ |  $8$ $|$ $8$ |  $7$ $|$ $6$
$0.675$ |  $14$ $|$ $7$ |  $14$ $|$ $10$ |  $13$ $|$ $12$ |  $12$ $|$ $10$ |  $11$ $|$ $8$
$0.70$ |  $13$ $|$ $6$ |  $13$ $|$ $9$ |  $13$ $|$ $12$ |  $12$ $|$ $11$ |  $12$ $|$ $9$
$0.725$ |  $10$ $|$ $5$ |  $12$ $|$ $8$ |  $13$ $|$ $11$ |  $13$ $|$ $12$ |  $12$ $|$ $10$
$0.75$ |  $7$ $|$ $4$ |  $10$ $|$ $7$ |  $11$ $|$ $10$ |  $12$ $|$ $14$ |  $12$ $|$ $11$

####  5.2.4 Interpretation

In Table 1-4, it can be observed that the choice of $(\gamma_{1},\gamma_{2})$
has an important impact on the critical-time, varying from $4\text{s}$ to
$12\text{s}$ in Table 3 for instance. The interest of the critical-time metric
is then considered to be at least two-fold. First, it provides a means to
evaluate the maximum time allocated to defense mechanisms to detect and
mitigate anomalies. Second, it may be used as a prevention tool, as its
maximization leads to extra-time allowed to defense mechanisms. Indeed,
maximizing the critical-time may be viewed as a rational security allocation
approach [4], where the aim is to tune degrees of freedom in the controlled
system in order to make the security problem easier rather than developing
complex defense strategies.

Finally, the estimated critical-time associated with channel-dependent sharp
anomalies is especially useful for the design of the channel-isolation step,
required for the selection of a suitable mitigation strategy, in a detection-
and-isolation scheme [15]. Moreover, from a risk analysis perspective, it
allows to identify which channel should be protected as a priority for
defender with limited resources [6]. For instance, channel C2 appears here
more vulnerable than C1 as it generally leads to lower critical-time (Table
4).

##  6 Conclusion

In this paper, the critical-time metric was investigated for controlled
systems subject to sharp input anomalies. This metric has several advantages
for off-line risk analysis as it allows for quantifying the time-window
available to defense mechanisms during which safe operation can be ensured.
Using the QC framework, an iterative LMI-based algorithm has been established
to compute the critical-time for an uncertain linear discrete-time model.
Several classes of sharp input anomalies can be considered, depending on the
input channel and the set of values the abnormal signal can take. In addition,
the potential of the critical-time for risk analysis and treatment has been
illustrated on the quadruple tank benchmark through the analysis of a
representative set of distinct sharp anomalies scenarios.

Future work includes extension to LPV models, the investigation of the
critical-time for stealthy attacks, and the integration of specific
communication models.

## References

  * [1] E. A. Lee, The Past, Present and Future of Cyber-Physical Systems: A Focus on Models, Sensors 15 (3) (2015) 4837–4869.
  * [2] S. M. Dibaji, M. Pirani, D. Flamholz, A. M. Annaswamy, K. H. Johansson, A. Chakrabortty, A Systems and Control Perspective of CPS Security, Annual Reviews in Control 47 (2019).  doi:https://doi.org/10.1016/j.arcontrol.2019.04.011.
  * [3] H. Sandberg, S. Amin, K. H. Johansson, Cyberphysical Security in Networked Control Systems: An Introduction to the Issue, IEEE Control Systems Magazine 35 (1) (2015) 20–23.  doi:10.1109/MCS.2014.2364708.
  * [4] M. S. Chong, H. Sandberg, A. Teixeira, A Tutorial Introduction to Security and Privacy for Cyber-Physical Systems, in: 2019 18th European Control Conference (ECC), 2019, pp. 968–978.  doi:10.23919/ECC.2019.8795652.
  * [5] D. Wang, Z. Wang, B. Shen, F. E. Alsaadi, T. Hayat, Recent Advances on Filtering and Control for Cyber-Physical Systems under Security and Resource Constraints, Journal of the Franklin Institute 353 (11) (2016) 2451–2466.
  * [6] A. Teixeira, K. C. Sou, H. Sandberg, K. H. Johansson, Secure Control Systems: A Quantitative Risk Management Approach, IEEE Control Systems Magazine 35 (1) (2015) 24–45.  doi:10.1109/MCS.2014.2364709.
  * [7] A. Teixeira, H. Sandberg, K. H. Johansson, Strategic Stealthy Attacks: The output-to-output l2-gain, in: 2015 54th IEEE Conference on Decision and Control (CDC), 2015, pp. 2582–2587.  doi:10.1109/CDC.2015.7402605.
  * [8] A. M. H. Teixeira, Security Metrics for Control Systems, in: Safety, Security and Privacy for Cyber-Physical Systems, Springer, 2021, pp. 99–121.
  * [9] Q. Zhu, T. Basar, Game-Theoretic Methods for Robustness, Security, and Resilience of Cyberphysical Control Systems, IEEE Control Systems Magazine 35 (1) (2015) 46–65.  doi:10.1109/MCS.2014.2364710.
  * [10] C. Escudero, P. Massioni, G. Scorletti, E. Zamaï, Security of Control Systems: Prevention of Aging Attacks by means of Convex Robust Simulation Forecasts, IFAC-PapersOnLine 53 (2020).
  * [11] C. Murguia, I. Shames, J. Ruths, D. Nešić, Security Metrics and Synthesis of Secure Control Systems, Automatica 115 (2020).
  * [12] H. Fawzi, P. Tabuada, S. Diggavi, Secure Estimation and Control for Cyber-Physical Systems Under Adversarial Attacks, IEEE Transactions on Automatic Control 59 (6) (2014) 1454–1467.  doi:10.1109/TAC.2014.2303233.
  * [13] J. Börner, F. Steinke, Measuring LTI System Resilience against Adversarial Disturbances based on Efficient Generalized Eigenvalue Computations, in: 2021 IEEE Conference on Decision and Control (CDC), 2021.
  * [14] A. Perodou, C. Combastel, A. Zolghadri, Critical-Time Analysis of Cyber-Physical Systems subject to Actuator Attacks and Faults, in: 2021 IEEE Conference on Decision and Control (CDC), 2021.
  * [15] I. Hwang, S. Kim, Y. Kim, C. E. Seah, A Survey of Fault Detection, Isolation, and Reconfiguration Methods, IEEE Transactions on Control Systems Technology 18 (3) (2010) 636–653.  doi:10.1109/TCST.2009.2026285.
  * [16] P. Griffioen, S. Weerakkody, B. Sinopoli, O. Ozel, Y. Mo, A Tutorial on Detecting Security Attacks on Cyber-Physical Systems, in: 18th European Control Conference (ECC), IEEE, 2019.
  * [17] Y. Zhang, J. Jiang, Bibliographical Review on Reconfigurable Fault-Tolerant Control Systems, Annual Reviews in Control 32 (2) (2008) 229–252.
  * [18] A. Perodou, C. Combastel, A. Zolghadri, Towards Anomaly-Tolerant Systems by Dissipation Block Synthesis, in: 2021 5th International Conference on Control and Fault-Tolerant Systems (SysTol), IEEE, 2021, pp. 19–24.
  * [19] R. Romagnoli, P. Griffioen, B. H. Krogh, B. Sinopoli, Software Rejuvenation under Persistent Attacks in Constrained Environments, IFAC-PapersOnLine 53 (2) (2020) 4088–4094.
  * [20] D. Ye, S. Luo, A Co-Design Methodology for Cyber-Physical Systems under Actuator Fault and Cyber Attack, Journal of the Franklin Institute 356 (4) (2019) 1856–1879.
  * [21] S. P. Boyd, L. El Ghaoui, E. Feron, V. Balakrishnan, Linear Matrix Inequalities in System and Control Theory, Vol. 15, SIAM, 1994.
  * [22] K. H. Johansson, The Quadruple-Tank Process: a Multivariable Laboratory Process with an Adjustable Zero, IEEE Transactions on Control Systems Technology 8 (3) (2000) 456–465.  doi:10.1109/87.845876.
  * [23] M. Althoff, G. Frehse, A. Girard, Set Propagation Techniques for Reachability Analysis, Annual Review of Control, Robotics, and Autonomous Systems 4 (2021) 369–395.
  * [24] U. T. Jönsson, A lecture on the S-procedure, Lecture Note at the Royal Institute of technology, Sweden 23 (2001) 34–36.
  * [25] G. Scorletti, Approche Unifiée de l’Analyse et de la Commande des Systemes par Formulation LMI, Ph.D. thesis (1997).
URL www.theses.fr/1997PA112160/document

  * [26] H. S. Sánchez, D. Rotondo, T. Escobet, V. Puig, J. Quevedo, Bibliographical Review on Cyber Attacks from a Control Oriented Perspective, Annual Reviews in Control 48 (2019) 103–128.
  * [27] H. S. Sanchez, D. Rotondo, T. Escobet, V. Puig, J. Saludes, J. Quevedo, Detection of Replay Attacks in Cyber-Physical Systems using a Frequency-based Signature, Journal of the Franklin Institute 356 (5) (2019) 2798–2824.
  * [28] M. Grant, S. Boyd, CVX: Matlab Software for Disciplined Convex Programming, version 2.2, http://cvxr.com/cvx (Mar. 2020).
  * [29] MOSEK ApS, The MOSEK optimization toolbox for MATLAB manual. Version 9.3 (2021).
  * [30] L. Schenato, B. Sinopoli, M. Franceschetti, K. Poolla, S. S. Sastry, Foundations of Control and Estimation Over Lossy Networks, Proceedings of the IEEE 95 (1) (2007) 163–187.  doi:10.1109/JPROC.2006.887306.


# Imprint of Early Dark Energy in Stochastic Gravitational Wave Background

Chia-Feng Chang<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of California, Riverside, CA 92521, USA

###### Abstract

Early dark energy that relieves Hubble tension leaves a fingerprint in the
primordial stochastic gravitational wave background that originates from
cosmic string network. We show that the signal is not only detectable with
future planned gravitational wave experiments, but also distinguishable from
other astrophysical and cosmological signals in the gravitational wave
frequency spectrum.

##  I Introduction

Gravitational waves (GWs) search technologies Abbott  _et al._ (2016a, b,
2017a, 2017b, 2017c) have strong ongoing research interest. Implications from
these developments focus not only on discovering new astrophysical objects
Abbott  _et al._ (2017d, e), but also probing gravitation physics Lombriser
and Lima (2017); Lombriser and Taylor (2016); Creminelli and Vernizzi (2017);
Ezquiaga and Zumalacárregui (2017); Baker  _et al._ (2017); Sakstein and Jain
(2017) and non-standard cosmologies Chang and Cui (2021, 2020a); Caldwell  _et
al._ (2019); Cui  _et al._ (2019, 2018); Battye and Shellard (1994); Blasi
_et al._ (2020); Gouttenoire  _et al._ (2020a); and GW observations has been
used to measure the current Hubble rate $H_{0}$ Abbott  _et al._ (2017f);
Mukherjee  _et al._ (2021); Gayathri  _et al._ (2020); Mukherjee  _et al._
(2020). Although not yet sufficiently precise to resolve the Hubble tension,
estimates will be continuously improved.

Increasingly attention is focusing on possible solutions of Hubble tension, a
Hubble rate discrepancy between local observations, such as supernovae signals
Riess  _et al._ (2016, 2018); Breuval  _et al._ (2020); Riess  _et al._
(2021); Soltis  _et al._ (2021); Freedman  _et al._ (2020), lensing time
delays Bonvin  _et al._ (2017); Birrer  _et al._ (2019); Denzel  _et al._
(2021); Yang  _et al._ (2020); Birrer  _et al._ (2020); Millon  _et al._
(2020); Baxter and Sherwin (2021) and non-local searches, e.g. the cosmic
microwave background (CMB) Aghanim  _et al._ (2020); Aylor  _et al._ (2017);
Aiola  _et al._ (2020); Choi  _et al._ (2020). Local measurements confirm
$H_{0}$ is approximately five sigma statistical significance above observation
in the CMB with assuming standard cosmological model ($\Lambda$CDM). This
strong disagreement has been closely examined using many statistical Feeney
_et al._ (2018); Zhang  _et al._ (2017); Cardona  _et al._ (2017); Bennett
_et al._ (2014) and measurement methods Pesce  _et al._ (2020); Wong  _et al._
(2020); Birrer  _et al._ (2020); Huang  _et al._ (2019); Freedman  _et al._
(2020, 2019), and re-examining potential technical issues Benedict  _et al._
(2007); Humphreys  _et al._ (2013); Efstathiou (2014); Rigault  _et al._
(2015); Kenworthy  _et al._ (2019); Spergel  _et al._ (2015) in both local and
CMB measurements (see reviews Knox and Millea (2020); Verde  _et al._ (2019);
Wong  _et al._ (2020); Aylor  _et al._ (2019)). Various evidence suggests the
discrepancy arises due to a currently unknown physical phenomenon outside
conventional $\Lambda$CDM predictions.

V. Poulin et al. Poulin  _et al._ (2019, 2018); Smith  _et al._ (2020a) showed
that early dark energy (EDE) behaves as a cosmological constant when comoving
scaling factor $a(t)$ is smaller than critical $a_{c}\sim 10^{-4}$, and is
then diluted faster or equal to radiation-like component to relieve the Hubble
tension. EDE contributes up to $20\%$ energy density of the universe (model-
dependent) at $a_{c}$, consequently accelerating the universe expansion and
hence slightly delaying the universe entering the matter-domination era. This
framework brings $H_{0}$ in estimated CMB to be consistent with local
measurements and also with measured high and low redshifts Agrawal  _et al._
(2019); Smith  _et al._ (2020b); Murgia  _et al._ (2021); Lin  _et al._
(2020).

We adopted the simplest EDE model, i.e., slow-rolling potential
$V(\phi)\propto\phi^{2n}$ with scalar field $\phi$ Copeland  _et al._ (2006);
Linder (2008); Chang and Cui (2020b). Effective mass for $\phi$ is lighter
than the Hubble rate on the early universe, and hence Hubble friction
overdamps scalar field motion and freezes it. Consequently, the scalar field
behaves as a subdominant cosmological constant until the Hubble rate decreases
to approximately scalar effective mass, namely, the driving force overcomes
the Hubble friction. Subsequently, the field starts oscillating as a fluid
with equation of state $w_{\phi}=(n-1)/(n+1)$.

Cosmic strings, one dimension long-lived topological defects, are stable and
predictable sources for stochastic GW background (SGWB) and hence ideal
creator sources for GW as a messenger that carries new signal from the early
universe Vilenkin (1985); Caldwell and Allen (1992); Hindmarsh and Kibble
(1995). Such stable objects arise from beyond standard model theories, such as
spontaneously broken $U(1)$ Kibble (1976); Nielsen and Olesen (1973);
Vachaspati and Vilenkin (1984); Vilenkin and Shellard (2000) or superstring
theories Copeland  _et al._ (2004); Dvali and Vilenkin (2004); Polchinski
(2004); Jackson  _et al._ (2005). The GW frequency spectrum formed by the
cosmic string network is an approximate plateau over a broad frequency range
with parameter $G\mu$ dependence, where $G$ is the Newtonian gravitational
constant and $\mu$ is string tension. GW experiments EPTA Lentati  _et al._
(2015) and PPTA Shannon  _et al._ (2015) provide strong bounds on
$G\mu\lesssim 2\times 10^{-11}$ Cui  _et al._ (2019), but this remains in
tension with the NANOGrav 12.5yr result Arzoumanian  _et al._ (2020)
$G\mu\in(2,30)\times 10^{-11}$ in 95% C.L. Ellis and Lewicki (2021); Blasi
_et al._ (2021). The tension may be due to the different noise analysis
methods Hazboun  _et al._ (2020).

This paper shows that EDE imprints a GW signal in the SGWB formed from the
cosmic string network that is distinguishable from other astrophysical and
cosmological signals in the GW frequency spectrum. Such a unique spectrum
could be detected in future GW experiments LISA Amaro-Seoane  _et al._ (2017);
Bartolo  _et al._ (2016); Caprini  _et al._ (2019) and SKA Janssen  _et al._
(2015). We also contour the SKA and LISA sensitivities on EDE parameter space.

Figure 1:  Signal difference ratio versus GW frequencies as defined in Eq.(6).
The black solid curve: $G\mu=10^{-12}$, $a_{c}=10^{-4.48}$, $n=2$ and $f_{\rm
EDE}=4\%$. According to black curve, others are changing one parameter on
each, e.g. red: change $a_{c}\to 10^{-3.57}$, green: $G\mu\to 10^{-11}$, and
blue: $n\to 3$. The $a_{c}$ correspond to highest and lowest values in
$68\%$C.L. CMB analysis Poulin  _et al._ (2019). LISA and SKA sensitivities
targeting $G\mu=10^{-12}$ cosmic string GW spectrum (see Fig. 2 in Cui  _et
al._ (2019)) show as gray and yellow area, respectively. The lighter yellow is
for $G\mu=10^{-11}$.

##  II Framework

Nambu-Goto string action is an effective model that describes cosmic string
network evolution Auclair  _et al._ (2020). These strings formed at time
$t_{F}$ when temperature cools to symmetry breaking scale of theory. Shortly
after, the defects behave as a scaling invariant network that includes a few
long (super-horizon length) strings and a collection of closed loops chopped
from long strings. GW loop emission dominates the network energy loss, and
loop number density $n_{o}(t_{i},t)$ is characterized by Vilenkin and Shellard
(2000)

$\displaystyle dn_{o}(t_{i},t)=\frac{0.1}{\alpha}\int_{t_{F}}^{t}C_{\rm
eff}(t_{i})\frac{dt_{i}}{t_{i}^{4}}\left(\frac{a(t_{i})}{a(t)}\right)^{3},$
(1)

where loops form at $t_{i}$ and subsequently continuously diluted until time
$t$, factor $0.1$ represents $90\%$ string energy release to loop kinetic
energy and subsequent diluted off, and parameter $\alpha=0.1$ Blanco-Pillado
_et al._ (2014); Blanco-Pillado and Olum (2017) characterize initial loop size
$\ell(t_{i})=\alpha t_{i}$. The appendix reviews the simulation calibration
parameter $C_{\hbox{\scriptsize{eff}}}$ Cui  _et al._ (2018, 2019), which
controls loop production. Created loops are shortened by radiating GWs

$\displaystyle\ell(t)=\alpha t_{i}-\Gamma
G\mu(t-t_{i}),\;\;\;\;\;\;\hbox{with}\;\;\;\;t\geq t_{i},$ (2)

where GW emission parameter $\Gamma=50$ Blanco-Pillado  _et al._ (2014);
Blanco-Pillado and Olum (2017); Vilenkin (1981); Turok (1984); Quashnock and
Spergel (1990). Loop emission characterizes in normal modes
$k\in\mathbb{Z}^{+}$ oscillation with emitted GW frequency,

$\displaystyle f=\frac{2k}{\ell}=\frac{a(t)}{a(t_{0})}\frac{2k}{\alpha
t_{i}-\Gamma G\mu(t-t_{i})},$ (3)

where emitted GWs redshift to today ($t_{0}$).

Figure 2:  Energy density fraction of early dark energy
$f_{\hbox{\scriptsize{EDE}}}$ versus string tension parameter $G\mu$ with peak
frequency $f_{p}$ on upper $x$-axis, and fixing early dark energy potential
exponent parameter $n=2$ and critical redshift $a_{c}=10^{-4.48}$. The colored
region shows influenced cosmic string GW background Signal-to-Noise Ratio
(SNR) with LISA 4 years nominal mission operating period, see Eq.(10). Figure
3:  GW spectrum difference ratio $\Delta\Omega_{\rm GW}/\Omega_{\rm GW}$ at
peak frequency $f_{p}$ with fixed $G\mu=10^{-12}$. The colored areas and the
black circled areas represent the 95% C.L. and 68% C.L. in CMB numerical
analysis, respectively. The CMB analysis is directly quoted from Poulin  _et
al._ (2019), they scan early dark energy (EDE) total energy fraction $1\%\leq
f_{\rm EDE}\leq 13\%$ and critical scaling factor $10^{-4.8}\leq a_{c}\leq
10^{-3.4}$ in different EDE potential $\phi^{2n}$ exponents
$n=\\{2,3,\infty\\}$. The dashed curves show the signal-noise-ratio
$\hbox{SNR}=1$ with LISA 4 years collection. The curves on the boundaries
between each color (e.g. between lighter green and green) are $a_{c}$
variations with fixed $f_{\rm EDE}$ as the red dashed line.

Thus, the SGWB frequency spectrum can be computed by summing all normal modes,

$\displaystyle\Omega_{\hbox{\scriptsize{GW}}}(f)=\frac{f}{\rho_{c}}\frac{d\rho_{GW}}{df}=\sum_{k}\Omega_{\hbox{\scriptsize{GW}}}^{(k)}(f),$
(4)

with

$\displaystyle\Omega_{\hbox{\scriptsize{GW}}}^{(k)}(f)=\frac{1}{\rho_{c}}\frac{2k}{f}\frac{F_{\alpha}\Gamma^{(k)}G\mu^{2}}{\alpha(\alpha+\Gamma
G\mu)}$ (5)
$\displaystyle\times\int_{t_{F}}^{t_{0}}dt\frac{C_{\hbox{\scriptsize{eff}}}(t_{i})}{t_{i}^{4}}\left(\frac{a(t)}{a(t_{0})}\right)^{5}\left(\frac{a(t_{i})}{a(t)}\right)^{3}\theta(t_{i}-t_{F})\theta(\ell(t)),$

where critical density $\rho_{c}=3H_{0}^{2}/8\pi G$, cusp dominates GW
emission $\Gamma^{(k)}=\Gamma/(3.6\,k^{4/3})$ Cui  _et al._ (2019), and we sum
over $k$ modes up to $k\leq 10^{5}$.

To visualize EDE influence on SGWB, we define signal difference to SGWB
spectrum ratio as,

$\displaystyle\frac{\Delta\Omega_{\rm GW}}{\Omega_{\rm
GW}}(f)\equiv\frac{\Delta\Omega_{\rm GW}(f)}{\Omega_{\rm
GW}^{\scriptsize{\Lambda\hbox{CDM}}}(f)}\equiv\frac{\Omega_{\rm GW}^{\rm
EDE}(f)-\Omega_{\rm GW}^{\scriptsize{\Lambda\hbox{CDM}}}(f)}{\Omega_{\rm
GW}^{\scriptsize{\Lambda\hbox{CDM}}}(f)},$ (6)

where superscripts imply GW frequency spectra on different cosmologies. Fig. 1
shows that EDE modifies the spectrum in two frequency regions: firstly, in
contrast to $\Lambda$CDM, diluting EDE delays entering of matter domination,
which slows down the redshift on GWs that emitted by the loops with loop
formation time $t_{i}\lesssim\Gamma G\mu t_{c}/\alpha$ where $t_{c}$ is the
universe age at $a_{c}$. This mechanism peaks the GW spectrum at frequency
$f\sim 10^{-4}$ Hz. The resulting spectrum typically increases at
characterized frequency $f_{p}$, which can be estimated from Eq.(3) as,

$\displaystyle f_{p}\sim\frac{2\epsilon}{\Gamma G\mu
t_{c}}\frac{a_{c}}{a(t_{0})},$ (7)

where $\epsilon\sim\mathcal{O}(1)$ is a numerical parameter. The shape of this
peak approximately estimate to

$\displaystyle\frac{\Delta\Omega_{\rm GW}}{\Omega_{\rm
GW}}(f)\propto\left\\{\begin{aligned}
&\left(\frac{f}{f_{p}}\right)^{-0.23},\;\;\;\;\;\hbox{for}\;\;f\geq f_{p},\\\
&\left(\frac{f}{f_{p}}\right)^{3},\;\;\;\;\;\;\;\;\;\;\;\hbox{for}\;\;f<f_{p},\end{aligned}\right.$
(8)

and the slower diluting EDE (lower $n$) more significantly delays the universe
entering the matter domination epoch while GWs experience a longer and slower
dilution period, hence increasing the signal difference.

Second, faster universe expansion at $t_{c}$ reduces loop chopping efficiency
i.e., reduces $C_{\rm eff}(t_{i}=t_{c})$ in Eq.(1), due to less frequent
intercommutation between strings. Such a mechanism implies a signal difference
dip at characteristic frequency $f_{d}\sim 10^{-9}$ to $10^{-11}\,$Hz as shown
in Fig. 1, which can be computed from Eq.(3) with GWs emission today,

$\displaystyle f_{d}\sim\frac{2}{\alpha t_{c}}.$ (9)

Shortly after $t_{c}$, more long strings have entered the horizon, loop
chopping efficiency turns back to about that for $\Lambda$CDM, and again, the
delayed entering of matter domination era. Therefore, loop number density
increases slightly, and redshift effects on GW slow. This mechanism increases
signal difference at frequencies below but close to $f_{d}$.

Universal fitting on a wide frequency spectrum not only determines string
tension parameter $G\mu$, but also addresses EDE parameters. $f_{\rm EDE}$
proportionally controls both signal difference peak and dip amplitudes,
whereas peak amplitude is sensitive to $n$ but dip amplitude is independent.
Therefore, we expect peak and dip amplitudes pin down $f_{\rm EDE}$ and $n$,
and $f_{d}$ and $f_{p}$ can be used to address $a_{c}$ (as shown in Fig. 1).

Modified GW frequency spectra present in LISA Amaro-Seoane  _et al._ (2017);
Bartolo  _et al._ (2016); Caprini  _et al._ (2019) and SKA Janssen  _et al._
(2015) frequency sensitivities (Fig. 1) and hence we focus on analyzing
signal-to-noise ratio (SNR) for the remainder of this paper.

##  III Signal-to-Noise Ratio

Early dark energy has relatively small energy contribution in the early
universe, hence its influence GW variation is small compared to string SGWB,
but could still be larger than noise background if string tension is
sufficiently large. For comparing signal difference and noise background, we
define SNR as $\\{i={\rm SKA},{\rm LISA}\\}$

$\displaystyle\hbox{SNR}=\sqrt{T_{i}\int^{f_{\rm max}}_{f_{\rm
min}}df\left(\frac{\Delta\Omega_{\rm GW}(f)}{\Omega^{i}_{n}(f)}\right)^{2}},$
(10)

where the SKA and LISA are calculated separately, frequency is integrated on
the experimental sensitivity region, $T_{i}$ is the observation period, and
$\Omega^{i}_{\rm n}(f)$ is an effective strain noise spectral spectrum for
LISA and SKA, respectively. SNR analysis has been widely studied Schmitz
(2021); Smith and Caldwell (2019); Thrane and Romano (2013); Romano and
Cornish (2017) and numerical details are provided in the appendix.

The peak frequency $f_{p}$ is within the LISA goal sensitivities
$20\mu\,$Hz$\,\leq f\leq 1\,$Hz Amaro-Seoane  _et al._ (2017) for $G\mu\leq
2\times 10^{-11}$ as presented in Fig. 2, and the peaky signal is more
significant than noise background for most of interesting parameter space.
More string energy causes stronger cosmic string SGWB, and hence the signal is
more significant for larger $G\mu$.

Fig. 3 shows signal difference ratio $\Delta\Omega_{\rm GW}/\Omega_{\rm GW}$
at peak frequency $f_{p}$. The colored areas were studied in CMB analysis
Poulin  _et al._ (2019) for relieving the Hubble tension. We explicitly show
that the signal difference could exceed the noise background (SNR$>1$) in the
4 years LISA collection, and therefore the EDE signal is detectable. As
discussed, lower $n$ has a larger signal difference due to the later universe
transition to matter-domination era.

Fig. 4 shows LISA and SKA detection sensitivities to EDE parameter space.
Signal difference can be larger than noise background over few frequency
decades, including most sensitive LISA regions, see the gentle slope
$f^{-0.23}$ over frequencies $f>f_{p}$ in Eq.(8). Thus, LISA can capture the
EDE signal even though $f_{p}$ is outside LISA sensitivities. Higher $a_{c}$
corresponds to lower SNR because $f_{p}$ moves away from LISA sensitivities
with increasing $a_{c}$.

Figure 4:  Energy density fraction of early dark energy $f_{\rm EDE}$ versus
critical scaling factor $a_{c}$ with marking SNR for LISA 4 years operating
and SKA 20 years observation period, respectively. The black areas present the
CMB analysis with $68\%\,$C.L. from Poulin  _et al._ (2019).

The pulsar timing detector SKA can capture (SNR$>1$) the signal difference
dip, as presented in Fig. 4. There are two falls on the SKA curves on the
left-panel: the fall on lower $a_{c}$ is caused by the dip in signal
difference; whereas the fall on higher $a_{c}$ is due to the signal difference
peak. If the NANOGrav signal is due to cosmic string i.e. $G\mu\geq 2\times
10^{-11}$ Ellis and Lewicki (2021); Blasi  _et al._ (2021), then SKA and LISA
should detect EDE signals in future observations. Scanning the $1\sigma$
regions, detectable EDE signal with $n=\\{2,3,\infty\\}$ requires
approximately at least $G\mu\gtrsim\\{4.4,6.7,31\\}\times 10^{-14}$ and
$G\mu\gtrsim\\{4.2,20,29\\}\times 10^{-13}$ for LISA and SKA, respectively.

##  IV Distinguish signal from other sources

The EDE signal in cosmic string SGWB is distinguishable from other possible
influences, such as string parameter variations, or sub-dominated SGWB from
astrophysical or cosmological objects. For example, kinks or kink-kink
collision modes Auclair  _et al._ (2020), or string parameters $\alpha$ and
$\Gamma$ variations Cui  _et al._ (2019) would universally influence GW
amplitude. Therefore, they cannot cause local amplitude modification as does
EDE. We numerically checked that EDE influenced SGWB difference is much
flatter than cosmic string SGWB with influence from astrophysical objects,
such as binary black hole merger Cholis (2017); Mandic  _et al._ (2016) and
inspiral Bonetti and Sesana (2020); Babak  _et al._ (2017); Amaro-Seoane  _et
al._ (2007). On the other hand, a sub-dominated GWs from cosmological
phenomena, e.g. flat frequency spectrum from inflation Turner  _et al._
(1993); Boyle and Steinhardt (2008); Caprini and Figueroa (2018); peaky
spectrum from domain wall Saikawa (2017); Hiramatsu  _et al._ (2014); or first
order phase transition dynamics Caprini  _et al._ (2009a); Alanne  _et al._
(2020); Schmitz (2021): sound wave Hindmarsh  _et al._ (2015), bubble wall
collision Huber and Konstandin (2008); Weir (2016) and magnetohydrodynamic
turbulence Caprini  _et al._ (2009b); Binetruy  _et al._ (2012), frequency
spectra are distinguishable to EDE as well. In particular, the dip structure
signal difference at $f_{d}$ cannot be caused by any other physical phenomena.

##  V Conclusion

This paper proposed that SGWB originating from cosmic string network can be
used to test the Hubble tension solution EDE. EDE behaves as dark energy in
the early universe, then begins to dilute at the critical scaling factor
$a_{c}\sim 10^{-4}$ with total energy density fraction $f_{\rm EDE}\gtrsim
1\%$. It influences cosmic string SGWB by accelerating and decelerating
universe expansion rate $\dot{a}(t)$ due to a dark energy like equation of
state $\omega_{\phi}=-1$ and diluting faster than or equal to radiation-like
component $\omega_{\phi}\geq 1/3$, respectively. The decelerated universe
expansion rate locally increases cosmic string GW frequency spectra with
magnitude $0.1\%$ to $1\%$ in the frequency range $10^{-5}$ to $10^{-3}\,$Hz
which is within LISA sensitivity. And on the other side, the acceleration
reduces spectra with magnitude $\sim 5\%$ in the lower frequency region $f\sim
10^{-9}$ Hz, within the SKA search region. We also showed that EDE influenced
signals were stronger than LISA and SKA experimental noise background, and
hence detectable. Such spectral shapes are distinguishable from other SGWB
sourced by cosmological or astrophysical objects. Thus, GW detection provides
a possibility to probe EDE.

We considered the simplest EDE effective model Poulin  _et al._ (2019), but
future study will extend to more complex models Poulin  _et al._ (2018);
Agrawal  _et al._ (2019); Smith  _et al._ (2020b); Murgia  _et al._ (2021);
Smith  _et al._ (2020a); Lin  _et al._ (2020) to see whether those can be
detected with this method.

###### Acknowledgements.

##  VI Acknowledgement

We thank Yanou Cui, Simeon Bird for valuable discussion. CC thanks Ming-Feng
Ho for helpful comments on the draft. The author also thanks Academia Sinica
for its hospitality. This project is in part supported by the Provost’s
Scholars for the Advancement of Physical Sciences fellowship.

###  VI.1 Appendix A: VOS Model

In this appendix, we review the velocity-dependent one-scale (VOS) model that
used to predict the string network evolution Martins and Shellard (1996a, b,
2002), and most of the content can also be found in Refs. Cui  _et al._ (2019,
2018); Avelino and Sousa (2012); Sousa and Avelino (2013); Gouttenoire  _et
al._ (2020b). The VOS model is used to describe the evolution of long Nambu-
Goto string network in terms of a mean string velocity Martins and Shellard
(1996a, b, 2002)

$\displaystyle\bar{v}=\sqrt{\frac{m}{2}\frac{k(\bar{v})}{\left[k(\bar{v})+\bar{c}\right]}\left(1-\frac{2}{m}\right)},$
(11)

and a characteristic length,

$\displaystyle\xi=\frac{m}{2}\sqrt{\frac{k(\bar{v})\left[k(\bar{v})+\bar{c}\right]}{2(m-2)}},$
(12)

as a fraction of the horizon, where the $m$ is exponent of scaling factor in
universe energy density $\rho\propto a^{-m}$, the chopping parameter
$\bar{c}=0.23$ Martins and Shellard (2002), and the ansatz function Martins
and Shellard (2002)

$\displaystyle
k(\bar{v})=\frac{2\sqrt{2}}{\pi}(1-\bar{v}^{2})\left(1+2\sqrt{2}\bar{v}^{3}\right)\frac{1-8\bar{v}^{6}}{1+8\bar{v}^{6}}.$
(13)

The long string energy density express

$\displaystyle\rho_{L}=\frac{\mu}{\left(\xi t\right)^{2}},$ (14)

and the intercommutation between strings chops them to loops as energy losing
as

$\displaystyle\frac{d\rho_{L}}{dt}=\bar{c}\bar{v}\frac{\rho_{L}}{\xi
t}=\bar{c}\bar{v}\frac{\mu}{\left(\xi t\right)^{3}}.$ (15)

Consequently, the loops number density at evolution time $t$ with particular
loop size $\ell(t_{i})=\alpha t_{i}$ that formed at time $t_{i}$ reads

$\displaystyle
dn_{o}(t_{i},t)=\frac{0.1}{\alpha}\int_{t_{F}}^{t}C_{\hbox{\scriptsize{eff}}}(t_{i})\frac{dt_{i}}{t_{i}^{4}}\left(\frac{a(t_{i})}{a(t)}\right)^{3},$
(16)

with

$\displaystyle C_{\rm eff}(t_{i})$
$\displaystyle=\frac{\bar{c}}{\gamma}\bar{v}\xi^{-3}$ (17)

where $C_{\rm eff}$ implies that the energy gain from string network, and
$a^{3}$ is due to number density dilution. The $\gamma=\sqrt{2}$ is loop
Lorentz boost Blanco-Pillado  _et al._ (2014); Blanco-Pillado and Olum (2017).
EDE locally influences the universe expansion rate around $a(t)\sim a_{c}$,
and therefore it would influence the $C_{\rm eff}$ on such a period. A smaller
$C_{\rm eff}$ implies a faster universe expansion rate, i.e. smaller $m$. It
is because of a less frequently intercommuation rate in a faster expanding
universe.

We present $C_{\rm eff}$ versus scaling factor $a$ in Fig. 5, the bumps on the
black curve are caused by the changing relativistic degree of freedom in early
universe. EDE locally influences $C_{\rm eff}$ around $a_{c}$ as the colored
ranges. Before $a_{c}$, EDE behaves as a cosmological constant that
accelerates the universe expansion rate, and therefore $C_{\rm eff}$ is
smaller than the one in $\Lambda$CDM around the $a_{c}$. Shortly after, EDE
starts diluting with a dilution rate that is faster or equals to radiation-
like components, then the universe expansion rate $m$ turns to slightly larger
than the $m$ in $\Lambda$CDM. Consequently, $C_{\rm eff}$ increases to
slightly higher than in $\Lambda$CDM. All curves converge when EDE energy
density isn’t comparable to other components.

Figure 5:  $C_{\rm eff}$ versus scaling factor $a$. The black curve is in
$\Lambda$CDM. The modified universes with the EDE are shown as
$a_{c}=10^{-4.48}$ and $a_{c}=10^{-3.57}$ as blue and red, respectively.

###  VI.2 Appendix B: Numerical Details

####  VI.2.1 Parameters

The cosmological parameters we used in numerical is as following: the scale
factor for Hubble expansion rate $h=0.71$ Poulin  _et al._ (2019),
pressureless matter density of the universe $\Omega_{m}=0.31$, dark energy
density of the $\Lambda$CDM universe $\Omega_{\Lambda}=0.69$, today
temperature $T_{0}=2.726\,$K, and relativistic degrees of freedom $g_{*}(T)$
is quoted from Saikawa and Shirai (2018).

####  VI.2.2 LISA noise background

To estimate the noise background in the LISA to a stochastic gravitational
wave background, we express the result that concluded in Smith and Caldwell
(2019). The signal-to-noise ratio is defined as $\\{i={\rm SKA},{\rm LISA}\\}$

$\displaystyle\hbox{SNR}=\sqrt{T_{i}\int^{f_{\rm max}}_{f_{\rm
min}}df\left(\frac{\Omega_{\rm GW}(f)}{\Omega^{i}_{n}(f)}\right)^{2}},$ (18)

where $\Omega_{\rm GW}(f)$ represents the GW signal spectrum, $T_{i}$ is
experiment operating period, and the integration range is the experiment
sensitive region that is $20\mu\,$Hz $\leq f\leq$ 1 Hz in LISA. The effective
strain noise spectral spectrum is expressed as

$\displaystyle\Omega_{n}^{i}(f)=S_{n}^{i}(f)\frac{2\pi^{2}f^{3}}{3H_{0}^{2}},$
(19)

with the effective noise power spectrum,

$\displaystyle S_{n}^{\rm LISA}(f)\simeq
2\sqrt{2}\frac{20}{3}\left[\frac{S_{I}(f)}{(2\pi
f)^{4}}+S_{II}\right]\left[1+\left(\frac{3f}{4f_{*}}\right)^{2}\right]$ (20)

where $f_{*}=c/(2\pi L)$ with light speed $c$ and $L=2.5\times 10^{6}$ km, and
$S_{II}=3.6\times 10^{-41}\,$Hz-1 is an optical path-length fluctuation. The
acceleration noise $S_{I}(f)$ reads

$\displaystyle S_{I}(f)=5.76\times
10^{-48}\left[1+\left(\frac{f_{1}}{f}\right)^{2}\right]\hbox{s}^{-4}\hbox{Hz}^{-1},$
(21)

where $f_{1}=0.4\,$mHz. As discussed in Smith and Caldwell (2019), the SNR in
Eq.(18) is an idealization, we ignore foreground contamination, non-
Gaussianity, data interruption, and other systematic issues.

####  VI.2.3 SKA noise background

The noise background of interferometer and pulsar timing experiments have been
nicely reviewed in the appendix of Schmitz (2021). As the result, the
effective strain noise spectral spectrum is defined as Eq.(19), and the
effective noise power spectrum is given Moore  _et al._ (2015a); Hazboun  _et
al._ (2019); Hellings and Downs (1983)

$\displaystyle S_{n}^{\rm SKA}(f)=12\pi^{2}f^{2}\left[\frac{2}{N_{\rm
SKA}(N_{\rm SKA}-1)}\right]^{1/2}\frac{D^{\rm SKA}_{n}}{\zeta_{\rm rms}},$
(22)

where $\zeta_{\rm rms}\simeq 0.147$, number of pulsars $N_{\rm SKA}=50$
Janssen  _et al._ (2015); Moore  _et al._ (2015b), operating time $T_{\rm
SKA}=20$yrs Schmitz (2021); Siemens  _et al._ (2013), and the 20 years timing
noise spectra $D^{\rm SKA}_{n}\simeq 1.1\times 10^{-9}\,$Hz-3 Schmitz (2021);
Thrane and Romano (2013). The sensitive frequency range is from the operating
period $f_{\rm min}=1/T_{\rm SKA}$ to the cadence of the timing observation,
$f_{\rm max}=1/T_{c}$ where $T_{c}=1\,$week Janssen  _et al._ (2015); Moore
_et al._ (2015b).

## References

  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, VIRGO), Phys. Rev. Lett. 118, 221101 (2017a), [Erratum: Phys.Rev.Lett. 121, 129901 (2018)], arXiv:1706.01812 [gr-qc] .
  * Abbott  _et al._ (2017b) B. . P. . Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 851, L35 (2017b), arXiv:1711.05578 [astro-ph.HE] .
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 141101 (2017c), arXiv:1709.09660 [gr-qc] .
  * Abbott  _et al._ (2017d) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017d), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2017e) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017e), arXiv:1710.05833 [astro-ph.HE] .
  * Lombriser and Lima (2017) L. Lombriser and N. A. Lima, Phys. Lett. B 765, 382 (2017), arXiv:1602.07670 [astro-ph.CO] .
  * Lombriser and Taylor (2016) L. Lombriser and A. Taylor, JCAP 03, 031 (2016), arXiv:1509.08458 [astro-ph.CO] .
  * Creminelli and Vernizzi (2017) P. Creminelli and F. Vernizzi, Phys. Rev. Lett. 119, 251302 (2017), arXiv:1710.05877 [astro-ph.CO] .
  * Ezquiaga and Zumalacárregui (2017) J. M. Ezquiaga and M. Zumalacárregui, Phys. Rev. Lett. 119, 251304 (2017), arXiv:1710.05901 [astro-ph.CO] .
  * Baker  _et al._ (2017) T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, Phys. Rev. Lett. 119, 251301 (2017), arXiv:1710.06394 [astro-ph.CO] .
  * Sakstein and Jain (2017) J. Sakstein and B. Jain, Phys. Rev. Lett. 119, 251303 (2017), arXiv:1710.05893 [astro-ph.CO] .
  * Chang and Cui (2021) C.-F. Chang and Y. Cui, (2021), arXiv:2106.09746 [hep-ph] .
  * Chang and Cui (2020a) C.-F. Chang and Y. Cui, Phys. Dark Univ. 29, 100604 (2020a), arXiv:1910.04781 [hep-ph] .
  * Caldwell  _et al._ (2019) R. R. Caldwell, T. L. Smith, and D. G. E. Walker, Phys. Rev. D 100, 043513 (2019), arXiv:1812.07577 [astro-ph.CO] .
  * Cui  _et al._ (2019) Y. Cui, M. Lewicki, D. E. Morrissey, and J. D. Wells, JHEP 01, 081 (2019), arXiv:1808.08968 [hep-ph] .
  * Cui  _et al._ (2018) Y. Cui, M. Lewicki, D. E. Morrissey, and J. D. Wells, Phys. Rev. D 97, 123505 (2018), arXiv:1711.03104 [hep-ph] .
  * Battye and Shellard (1994) R. A. Battye and E. P. S. Shellard, Phys. Rev. Lett. 73, 2954 (1994), [Erratum: Phys.Rev.Lett. 76, 2203–2204 (1996)], arXiv:astro-ph/9403018 .
  * Blasi  _et al._ (2020) S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Res. 2, 043321 (2020), arXiv:2004.02889 [hep-ph] .
  * Gouttenoire  _et al._ (2020a) Y. Gouttenoire, G. Servant, and P. Simakachorn, JCAP 07, 016 (2020a), arXiv:1912.03245 [hep-ph] .
  * Abbott  _et al._ (2017f) B. P. Abbott _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), Nature 551, 85 (2017f), arXiv:1710.05835 [astro-ph.CO] .
  * Mukherjee  _et al._ (2021) S. Mukherjee, G. Lavaux, F. R. Bouchet, J. Jasche, B. D. Wandelt, S. M. Nissanke, F. Leclercq, and K. Hotokezaka, Astron. Astrophys. 646, A65 (2021), arXiv:1909.08627 [astro-ph.CO] .
  * Gayathri  _et al._ (2020) V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. Lousto, and R. O’Shaughnessy, (2020), arXiv:2009.14247 [astro-ph.HE] .
  * Mukherjee  _et al._ (2020) S. Mukherjee, A. Ghosh, M. J. Graham, C. Karathanasis, M. M. Kasliwal, I. Magaña Hernandez, S. M. Nissanke, A. Silvestri, and B. D. Wandelt, (2020), arXiv:2009.14199 [astro-ph.CO] .
  * Riess  _et al._ (2016) A. G. Riess _et al._ , Astrophys. J. 826, 56 (2016), arXiv:1604.01424 [astro-ph.CO] .
  * Riess  _et al._ (2018) A. G. Riess _et al._ , Astrophys. J. 861, 126 (2018), arXiv:1804.10655 [astro-ph.CO] .
  * Breuval  _et al._ (2020) L. Breuval _et al._ , Astron. Astrophys. 643, A115 (2020), arXiv:2006.08763 [astro-ph.SR] .
  * Riess  _et al._ (2021) A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Astrophys. J. Lett. 908, L6 (2021), arXiv:2012.08534 [astro-ph.CO] .
  * Soltis  _et al._ (2021) J. Soltis, S. Casertano, and A. G. Riess, Astrophys. J. Lett. 908, L5 (2021), arXiv:2012.09196 [astro-ph.GA] .
  * Freedman  _et al._ (2020) W. L. Freedman, B. F. Madore, T. Hoyt, I. S. Jang, R. Beaton, M. G. Lee, A. Monson, J. Neeley, and J. Rich,  (2020), 10.3847/1538-4357/ab7339, arXiv:2002.01550 [astro-ph.GA] .
  * Bonvin  _et al._ (2017) V. Bonvin _et al._ , Mon. Not. Roy. Astron. Soc. 465, 4914 (2017), arXiv:1607.01790 [astro-ph.CO] .
  * Birrer  _et al._ (2019) S. Birrer _et al._ , Mon. Not. Roy. Astron. Soc. 484, 4726 (2019), arXiv:1809.01274 [astro-ph.CO] .
  * Denzel  _et al._ (2021) P. Denzel, J. P. Coles, P. Saha, and L. L. R. Williams, Mon. Not. Roy. Astron. Soc. 501, 784 (2021), arXiv:2007.14398 [astro-ph.CO] .
  * Yang  _et al._ (2020) T. Yang, S. Birrer, and B. Hu, Mon. Not. Roy. Astron. Soc. 497, L56 (2020), arXiv:2003.03277 [astro-ph.CO] .
  * Birrer  _et al._ (2020) S. Birrer _et al._ , Astron. Astrophys. 643, A165 (2020), arXiv:2007.02941 [astro-ph.CO] .
  * Millon  _et al._ (2020) M. Millon _et al._ , Astron. Astrophys. 639, A101 (2020), arXiv:1912.08027 [astro-ph.CO] .
  * Baxter and Sherwin (2021) E. J. Baxter and B. D. Sherwin, Mon. Not. Roy. Astron. Soc. 501, 1823 (2021), arXiv:2007.04007 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Aylor  _et al._ (2017) K. Aylor _et al._ (SPT), Astrophys. J. 850, 101 (2017), arXiv:1706.10286 [astro-ph.CO] .
  * Aiola  _et al._ (2020) S. Aiola _et al._ (ACT), JCAP 12, 047 (2020), arXiv:2007.07288 [astro-ph.CO] .
  * Choi  _et al._ (2020) S. K. Choi _et al._ (ACT), JCAP 12, 045 (2020), arXiv:2007.07289 [astro-ph.CO] .
  * Feeney  _et al._ (2018) S. M. Feeney, D. J. Mortlock, and N. Dalmasso, Mon. Not. Roy. Astron. Soc. 476, 3861 (2018), arXiv:1707.00007 [astro-ph.CO] .
  * Zhang  _et al._ (2017) B. R. Zhang, M. J. Childress, T. M. Davis, N. V. Karpenka, C. Lidman, B. P. Schmidt, and M. Smith, Mon. Not. Roy. Astron. Soc. 471, 2254 (2017), arXiv:1706.07573 [astro-ph.CO] .
  * Cardona  _et al._ (2017) W. Cardona, M. Kunz, and V. Pettorino, JCAP 03, 056 (2017), arXiv:1611.06088 [astro-ph.CO] .
  * Bennett  _et al._ (2014) C. L. Bennett, D. Larson, J. L. Weiland, and G. Hinshaw, Astrophys. J. 794, 135 (2014), arXiv:1406.1718 [astro-ph.CO] .
  * Pesce  _et al._ (2020) D. W. Pesce _et al._ , Astrophys. J. Lett. 891, L1 (2020), arXiv:2001.09213 [astro-ph.CO] .
  * Wong  _et al._ (2020) K. C. Wong _et al._ , Mon. Not. Roy. Astron. Soc. 498, 1420 (2020), arXiv:1907.04869 [astro-ph.CO] .
  * Huang  _et al._ (2019) C. D. Huang, A. G. Riess, W. Yuan, L. M. Macri, N. L. Zakamska, S. Casertano, P. A. Whitelock, S. L. Hoffmann, A. V. Filippenko, and D. Scolnic,  (2019), 10.3847/1538-4357/ab5dbd, arXiv:1908.10883 [astro-ph.CO] .
  * Freedman  _et al._ (2019) W. L. Freedman _et al._ ,  (2019), 10.3847/1538-4357/ab2f73, arXiv:1907.05922 [astro-ph.CO] .
  * Benedict  _et al._ (2007) G. F. Benedict, B. E. McArthur, M. W. Feast, T. G. Barnes, T. E. Harrison, R. J. Patterson, J. W. Menzies, J. L. Bean, and W. L. Freedman, Astron. J. 133, 1810 (2007), [Erratum: Astron.J. 133, 2980 (2007)], arXiv:astro-ph/0612465 .
  * Humphreys  _et al._ (2013) E. M. L. Humphreys, M. J. Reid, J. M. Moran, L. J. Greenhill, and A. L. Argon, Astrophys. J. 775, 13 (2013), arXiv:1307.6031 [astro-ph.CO] .
  * Efstathiou (2014) G. Efstathiou, Mon. Not. Roy. Astron. Soc. 440, 1138 (2014), arXiv:1311.3461 [astro-ph.CO] .
  * Rigault  _et al._ (2015) M. Rigault _et al._ , Astrophys. J. 802, 20 (2015), arXiv:1412.6501 [astro-ph.CO] .
  * Kenworthy  _et al._ (2019) W. D. Kenworthy, D. Scolnic, and A. Riess, Astrophys. J. 875, 145 (2019), arXiv:1901.08681 [astro-ph.CO] .
  * Spergel  _et al._ (2015) D. N. Spergel, R. Flauger, and R. Hložek, Phys. Rev. D 91, 023518 (2015), arXiv:1312.3313 [astro-ph.CO] .
  * Knox and Millea (2020) L. Knox and M. Millea, Phys. Rev. D 101, 043533 (2020), arXiv:1908.03663 [astro-ph.CO] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess, Nature Astron. 3, 891 (2019), arXiv:1907.10625 [astro-ph.CO] .
  * Aylor  _et al._ (2019) K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghunathan, and W. L. K. Wu, Astrophys. J. 874, 4 (2019), arXiv:1811.00537 [astro-ph.CO] .
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Poulin  _et al._ (2018) V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Phys. Rev. D 98, 083525 (2018), arXiv:1806.10608 [astro-ph.CO] .
  * Smith  _et al._ (2020a) T. L. Smith, V. Poulin, and M. A. Amin, Phys. Rev. D 101, 063523 (2020a), arXiv:1908.06995 [astro-ph.CO] .
  * Agrawal  _et al._ (2019) P. Agrawal, F.-Y. Cyr-Racine, D. Pinner, and L. Randall, (2019), arXiv:1904.01016 [astro-ph.CO] .
  * Smith  _et al._ (2020b) T. L. Smith, V. Poulin, J. L. Bernal, K. K. Boddy, M. Kamionkowski, and R. Murgia, (2020b), arXiv:2009.10740 [astro-ph.CO] .
  * Murgia  _et al._ (2021) R. Murgia, G. F. Abellán, and V. Poulin, Phys. Rev. D 103, 063502 (2021), arXiv:2009.10733 [astro-ph.CO] .
  * Lin  _et al._ (2020) M.-X. Lin, W. Hu, and M. Raveri, Phys. Rev. D 102, 123523 (2020), arXiv:2009.08974 [astro-ph.CO] .
  * Copeland  _et al._ (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Linder (2008) E. V. Linder, Gen. Rel. Grav. 40, 329 (2008), arXiv:0704.2064 [astro-ph] .
  * Chang and Cui (2020b) C.-F. Chang and Y. Cui, Phys. Rev. D 102, 015003 (2020b), arXiv:1911.11885 [hep-ph] .
  * Vilenkin (1985) A. Vilenkin, Phys. Rept. 121, 263 (1985).
  * Caldwell and Allen (1992) R. R. Caldwell and B. Allen, Phys. Rev. D 45, 3447 (1992).
  * Hindmarsh and Kibble (1995) M. B. Hindmarsh and T. W. B. Kibble, Rept. Prog. Phys. 58, 477 (1995), arXiv:hep-ph/9411342 .
  * Kibble (1976) T. W. B. Kibble, J. Phys. A9, 1387 (1976).
  * Nielsen and Olesen (1973) H. B. Nielsen and P. Olesen, Nucl. Phys. B61, 45 (1973).
  * Vachaspati and Vilenkin (1984) T. Vachaspati and A. Vilenkin, Phys. Rev. D 30, 2036 (1984).
  * Vilenkin and Shellard (2000) A. Vilenkin and E. P. S. Shellard,  _Cosmic Strings and Other Topological Defects_ (Cambridge University Press, 2000).
  * Copeland  _et al._ (2004) E. J. Copeland, R. C. Myers, and J. Polchinski, JHEP 06, 013 (2004), arXiv:hep-th/0312067 .
  * Dvali and Vilenkin (2004) G. Dvali and A. Vilenkin, JCAP 03, 010 (2004), arXiv:hep-th/0312007 .
  * Polchinski (2004) J. Polchinski, in  _NATO Advanced Study Institute and EC Summer School on String Theory: From Gauge Interactions to Cosmology_ (2004) arXiv:hep-th/0412244 .
  * Jackson  _et al._ (2005) M. G. Jackson, N. T. Jones, and J. Polchinski, JHEP 10, 013 (2005), arXiv:hep-th/0405229 .
  * Lentati  _et al._ (2015) L. Lentati _et al._ , Mon. Not. Roy. Astron. Soc. 453, 2576 (2015), arXiv:1504.03692 [astro-ph.CO] .
  * Shannon  _et al._ (2015) R. M. Shannon _et al._ , Science 349, 1522 (2015), arXiv:1509.07320 [astro-ph.CO] .
  * Arzoumanian  _et al._ (2020) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Ellis and Lewicki (2021) J. Ellis and M. Lewicki, Phys. Rev. Lett. 126, 041304 (2021), arXiv:2009.06555 [astro-ph.CO] .
  * Blasi  _et al._ (2021) S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Lett. 126, 041305 (2021), arXiv:2009.06607 [astro-ph.CO] .
  * Hazboun  _et al._ (2020) J. S. Hazboun, J. Simon, X. Siemens, and J. D. Romano, Astrophys. J. Lett. 905, L6 (2020), arXiv:2009.05143 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Bartolo  _et al._ (2016) N. Bartolo _et al._ , JCAP 12, 026 (2016), arXiv:1610.06481 [astro-ph.CO] .
  * Caprini  _et al._ (2019) C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, JCAP 11, 017 (2019), arXiv:1906.09244 [astro-ph.CO] .
  * Janssen  _et al._ (2015) G. Janssen _et al._ , PoS AASKA14, 037 (2015), arXiv:1501.00127 [astro-ph.IM] .
  * Auclair  _et al._ (2020) P. Auclair _et al._ , JCAP 04, 034 (2020), arXiv:1909.00819 [astro-ph.CO] .
  * Blanco-Pillado  _et al._ (2014) J. J. Blanco-Pillado, K. D. Olum, and B. Shlaer, Phys. Rev. D 89, 023512 (2014), arXiv:1309.6637 [astro-ph.CO] .
  * Blanco-Pillado and Olum (2017) J. J. Blanco-Pillado and K. D. Olum, Phys. Rev. D 96, 104046 (2017), arXiv:1709.02693 [astro-ph.CO] .
  * Vilenkin (1981) A. Vilenkin, Phys. Lett. B 107, 47 (1981).
  * Turok (1984) N. Turok, Nucl. Phys. B 242, 520 (1984).
  * Quashnock and Spergel (1990) J. M. Quashnock and D. N. Spergel, Phys. Rev. D 42, 2505 (1990).
  * Schmitz (2021) K. Schmitz, JHEP 01, 097 (2021), arXiv:2002.04615 [hep-ph] .
  * Smith and Caldwell (2019) T. L. Smith and R. Caldwell, Phys. Rev. D 100, 104055 (2019), arXiv:1908.00546 [astro-ph.CO] .
  * Thrane and Romano (2013) E. Thrane and J. D. Romano, Phys. Rev. D 88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Romano and Cornish (2017) J. D. Romano and N. J. Cornish, Living Rev. Rel. 20, 2 (2017), arXiv:1608.06889 [gr-qc] .
  * Cholis (2017) I. Cholis, JCAP 06, 037 (2017), arXiv:1609.03565 [astro-ph.HE] .
  * Mandic  _et al._ (2016) V. Mandic, S. Bird, and I. Cholis, Phys. Rev. Lett. 117, 201102 (2016), arXiv:1608.06699 [astro-ph.CO] .
  * Bonetti and Sesana (2020) M. Bonetti and A. Sesana, Phys. Rev. D 102, 103023 (2020), arXiv:2007.14403 [astro-ph.GA] .
  * Babak  _et al._ (2017) S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D 95, 103012 (2017), arXiv:1703.09722 [gr-qc] .
  * Amaro-Seoane  _et al._ (2007) P. Amaro-Seoane, J. R. Gair, M. Freitag, M. Coleman Miller, I. Mandel, C. J. Cutler, and S. Babak, Class. Quant. Grav. 24, R113 (2007), arXiv:astro-ph/0703495 .
  * Turner  _et al._ (1993) M. S. Turner, M. J. White, and J. E. Lidsey, Phys. Rev. D 48, 4613 (1993), arXiv:astro-ph/9306029 .
  * Boyle and Steinhardt (2008) L. A. Boyle and P. J. Steinhardt, Phys. Rev. D 77, 063504 (2008), arXiv:astro-ph/0512014 .
  * Caprini and Figueroa (2018) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Saikawa (2017) K. Saikawa, Universe 3, 40 (2017), arXiv:1703.02576 [hep-ph] .
  * Hiramatsu  _et al._ (2014) T. Hiramatsu, M. Kawasaki, and K. Saikawa, JCAP 02, 031 (2014), arXiv:1309.5001 [astro-ph.CO] .
  * Caprini  _et al._ (2009a) C. Caprini, R. Durrer, T. Konstandin, and G. Servant, Phys. Rev. D 79, 083519 (2009a), arXiv:0901.1661 [astro-ph.CO] .
  * Alanne  _et al._ (2020) T. Alanne, T. Hugle, M. Platscher, and K. Schmitz, JHEP 03, 004 (2020), arXiv:1909.11356 [hep-ph] .
  * Hindmarsh  _et al._ (2015) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 92, 123009 (2015), arXiv:1504.03291 [astro-ph.CO] .
  * Huber and Konstandin (2008) S. J. Huber and T. Konstandin, JCAP 09, 022 (2008), arXiv:0806.1828 [hep-ph] .
  * Weir (2016) D. J. Weir, Phys. Rev. D 93, 124037 (2016), arXiv:1604.08429 [astro-ph.CO] .
  * Caprini  _et al._ (2009b) C. Caprini, R. Durrer, and G. Servant, JCAP 12, 024 (2009b), arXiv:0909.0622 [astro-ph.CO] .
  * Binetruy  _et al._ (2012) P. Binetruy, A. Bohe, C. Caprini, and J.-F. Dufaux, JCAP 06, 027 (2012), arXiv:1201.0983 [gr-qc] .
  * Martins and Shellard (1996a) C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 53, 575 (1996a), arXiv:hep-ph/9507335 .
  * Martins and Shellard (1996b) C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 54, 2535 (1996b), arXiv:hep-ph/9602271 .
  * Martins and Shellard (2002) C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 65, 043514 (2002), arXiv:hep-ph/0003298 .
  * Avelino and Sousa (2012) P. P. Avelino and L. Sousa, Phys. Rev. D 85, 083525 (2012), arXiv:1202.6298 [astro-ph.CO] .
  * Sousa and Avelino (2013) L. Sousa and P. P. Avelino, Phys. Rev. D 88, 023516 (2013), arXiv:1304.2445 [astro-ph.CO] .
  * Gouttenoire  _et al._ (2020b) Y. Gouttenoire, G. Servant, and P. Simakachorn, JCAP 07, 032 (2020b), arXiv:1912.02569 [hep-ph] .
  * Saikawa and Shirai (2018) K. Saikawa and S. Shirai, JCAP 05, 035 (2018), arXiv:1803.01038 [hep-ph] .
  * Moore  _et al._ (2015a) C. J. Moore, S. R. Taylor, and J. R. Gair, Class. Quant. Grav. 32, 055004 (2015a), arXiv:1406.5199 [astro-ph.IM] .
  * Hazboun  _et al._ (2019) J. S. Hazboun, J. D. Romano, and T. L. Smith, Phys. Rev. D 100, 104028 (2019), arXiv:1907.04341 [gr-qc] .
  * Hellings and Downs (1983) R. w. Hellings and G. s. Downs, Astrophys. J. Lett. 265, L39 (1983).
  * Moore  _et al._ (2015b) C. J. Moore, R. H. Cole, and C. P. L. Berry, Class. Quant. Grav. 32, 015014 (2015b), arXiv:1408.0740 [gr-qc] .
  * Siemens  _et al._ (2013) X. Siemens, J. Ellis, F. Jenet, and J. D. Romano, Class. Quant. Grav. 30, 224015 (2013), arXiv:1305.3196 [astro-ph.IM] .