\section{Introduction}
Model Predictive Control (MPC) has become the industry standard for model-based optimal control, with widespread applications and technical developments \cite{qin2003survey}. The objective of an MPC is to solve a finite horizon optimisation problem with respect to some constraints \cite{fontes2018guaranteed,RaM09}. In chemical industries, MPC has seen great success due to its ability to deliver ``optimal'' performance and to explicitly handle constraints. Often, a trade-off is made between prediction horizon length and the use of detailed models (attracting increased computational burdens) against the acceptable performance of the controller.

Whilst many chemical processes can be modelled with a high level of accuracy by using more than one time scale \cite{chang1984multi,bailey2018biochemical}, the use of such models can be computationally challenging \cite{pinnamaraju2018empirical}. In some cases, this introduced complexity is necessary since controller design using a model in only one timescale may cause instability when applied to the physical system. It is usually also true for such systems that a long prediction horizon (e.g., in the case of MPC) is needed for effective control. However, due to computational limitations, long optimisation horizons combined with relatively short time intervals are generally impractical \cite{findeisen2016exploiting}.

To overcome these complexity issues, one approach is to modify the standard MPC formulation for the multi-timescale setting \cite{tan2016model}, by introducing timescale specific models (using non-uniform time intervals) and hence prediction horizons. In this approach, small time intervals (capturing fast process dynamics) and large time intervals (capturing slower dynamics) are used to ensure short term accuracy whilst facilitating long term predictions. As a result, the computational burden associated with detailed predictions over long horizons can be significantly reduced by using models with larger time intervals without sacrificing significant accuracy or introducing significant complexities. Traditional methods to ensure the closed-loop stability of MPC schemes involve using the cost as a Lyapunov function \cite{RaM09}, which relies on models with uniform time intervals, and hence a uniformly spaced optimisation horizon to show that the cost is decreasing.

If arbitrary cost functions are considered, this stability issue is further compounded and consequently requires a tailored stability condition \cite{ellis2017economic}. In addition, if the reference trajectory is considered time-varying, stability guarantees for an optimal controller require a condition that is also reference-independent, e.g., those based on incremental stability \cite{angeli2002lyapunov}. Introduced by \cite{lohmiller1998contraction}, contraction theory facilitates stability analysis and control of nonlinear systems with respect to time-varying (feasible) references without redesigning the control law \cite{manchester2017control}. Contraction theory is particularly attractive in that it can be used to analyse the incremental stability of nonlinear systems and simultaneously synthesise a tracking controller using control contraction metrics (see, e.g., \cite{manchester2017control}). This has motivated developments for contraction-based MPCs \cite{mccloy2021differential}, offering significant flexibility over Lyapunov-based alternatives.

Inspired by the work in \cite{tan2016model}, we propose a non-uniformly spaced optimisation horizon nonlinear MPC with timescale specific contraction-based stability constraints. Extending from \cite{tan2016model}, we relax the linear system requirements to consider discrete-time nonlinear systems and propose a novel \textit{modular} contraction-based stability constraint to ensure closed-loop stability across multiple time intervals (for multiple timescale specific models). Moreover, this construction is reference-flexible and scalable, permitting tracking of \textit{a priori} unknown time-varying operation targets, arbitrary economic cost functions, and inclusion of additional timescale models in the MPC without increasing the computational burden.

The remainder of this article is structured as follows, Section \ref{sec:pre} presents the prerequisite contraction theory tools, and Section \ref{sec:pro} formulates the multi-timescale MPC problem. Section \ref{sec:ctr} develops contraction-based controller structures for systems subject to multiple timescales and imposes them as constraints on a multi-timescale MPC. Section \ref{sec:exa} demonstrates the overall method via numerical simulation and Section \ref{sec:conclusion} concludes this article.

\section{Preliminaries} \label{sec:pre}
Firstly, consider the single timescale nonlinear system
\begin{equation}\label{equ:pre cer sys}
    x_{k+1} = f(x_k) + g(x_k)u_k,
\end{equation}
where state and control are $x_k \in \mathcal{X} \subseteq \mathbb{R}^n$ and $u_k \in \mathcal{U} \subseteq \mathbb{R}^m$. The corresponding differential system of \eqref{equ:pre cer sys} is
\begin{equation}\label{equ:pre cer dif sys}
    \delta_{x_{k+1}} = A\delta_{x_k} + B\delta_{u_k},
\end{equation}
where Jacobian matrices of $f$ and $g$ in \eqref{equ:pre cer sys} are defined as $A:=\frac{\partial (f(x_k) + g(x_k)u_k)}{\partial x_k}$ and $B:=\frac{\partial (f(x_k) + g(x_k)u_k)}{\partial u_k}$  respectively, $\delta_{u_k} := \frac{\partial u_k}{\partial s}$ and $\delta_{x_k}:=\frac{\partial x_k}{\partial s}$ are vectors in the tangent space $T_x\mathcal{U}$ at $u_k$ and tangent space $T_x\mathcal{X}$ at $x_k$ respectively, where $s$ parameterises a path, $c(s): [0,1] \rightarrow \mathcal{X}$ between two points such that $c(0) = x, c(1) = x^*  \in \mathcal{X}$ (see Fig. \ref{fig:rie_geo}). Consider a feedback control law for the differential dynamics \eqref{equ:pre cer dif sys} as
\begin{equation}\label{equ:pre dif fed}
    \delta_{u_k} = K(x_k) \delta_{x_k},
\end{equation}
where $K$ is a state dependent feedback gain. The contraction condition for a discrete-time nonlinear control affine system in \eqref{equ:pre cer sys} can be described as follows,
\begin{theorem}[\cite{wei2021contractionsyntehsis}]\label{thm:pre ctr con}
    For a discrete-time nonlinear system \eqref{equ:pre cer sys}, with differential dynamics \eqref{equ:pre cer dif sys} and differential state-feedback controller \eqref{equ:pre dif fed}, provided a uniformly metric, $\underline{m}I \leq M(x_k) \leq \overline{m}I$, exists satisfying,
    \begin{equation}\label{equ:pre ctr con}
        (A+BK)^\top M_+(A+BK) - (1-\beta)M < 0,
    \end{equation}
    where $M_+=M(x_{k+1})$, the closed-loop system is contracting for some constant $0 < \beta \leq 1$. Furthermore, the closed-loop system is exponentially incrementally stable,  i.e.,
    \begin{equation}
        \label{equ:pre exp sta}
        |x_k-x^*_k| \leq  R e^{-\lambda k} |x_0 - x^*_0|,
    \end{equation}
    for some constant $R$, convergence rate $\lambda$ and any feasible reference trajectory $(x^*, u^*)$, satisfying \eqref{equ:pre cer sys}, where $x_0$ is the initial state and $x_k$ is the state at time-step $k$.

\end{theorem}
If condition \eqref{equ:pre ctr con} holds for all feasible state values in a region, that region is called a contraction region. In Theorem \ref{thm:pre ctr con}, the metric, $M$, is used in describing Riemannian geometry, which we briefly introduce here. We define the Riemannian distance, $d(x,x^*)$, as (see, e.g., \cite{do1992riemannian})
\begin{equation}\label{equ:Riemannian distance and energy}
    \begin{aligned}
    d(x,x^*) = d(c) :=\int_0^1 \sqrt{\delta^\top_{c(s)}M(c(s))\delta_{c(s)}}ds,
    \end{aligned}
\end{equation}
where $\delta_{c(s)} := \frac{\partial c(s)}{\partial s}$.
The shortest path in Riemannian space, or \textit{geodesic}, between $x$ and $x^*$ is defined as
\begin{equation}\label{equ:geodesic}
    \gamma(s) :=\argmin_{c(s)} {d(x,x^*)}.
\end{equation}

Leveraging Riemannian tools, a feasible contracting tracking controller controller for \eqref{equ:pre cer sys}, is obtained by integrating the differential feedback law \eqref{equ:pre dif fed} along the geodesic, $\gamma(s)$ \eqref{equ:geodesic}, as
\begin{equation}\label{equ:pre ctl int}
    u_k = u^*_k + \int_0^1K(\gamma(s))\frac{\partial \gamma(s)}{\partial s}\,ds.
\end{equation}

In summary, a contracting tracking controller controller ensures that the length of the minimum path (i.e., geodesic) between any two trajectories (e.g., the plant state, $x$, and desired state, $x^*$, trajectories), with respect to the metric $M$, shrinks with time, i.e., provided that the contraction condition \eqref{equ:pre ctr con} holds for the discrete-time nonlinear system \eqref{equ:pre cer sys},
we can employ a stabilizing feedback controller \eqref{equ:pre ctl int} to ensure convergence to feasible operating targets.
\begin{figure} \label{fig:rie_geo}
    \begin{center}
        \includegraphics[width=\linewidth]{spath.pdf}
        \caption{State and reference trajectories with $s$-parameterised paths.}
    \end{center}
    \vspace{-2.5mm}\squeezeup
\end{figure}

\section{Problem Formulation} \label{sec:pro}

\subsection{Multiple Timescale Nonlinear Processes}
In many chemical processes, the timescales of the underlying dynamics are distributed over a very large range. As such, the process dynamics can be described as functions of appropriate time intervals. Herein, we consider the timescale dependent nonlinear process model
\begin{equation}\label{eq:sys_dist_nom}
    x_{k+1} = f(x_k,\tau) + B(\tau) u_k
\end{equation}
where $\tau$ denotes the discrete time interval. For processes described by \eqref{eq:sys_dist_nom}, the process dynamics vary based on the time interval in which it is viewed, e.g., fast dynamics over short intervals and slow dynamics over long time intervals. The primary goal is to drive \eqref{eq:sys_dist_nom} to feasible target trajectories $x^*_k$, described as \begin{equation}\label{eq:sys_ref_nom}
    x^*_{k+1} = f(x^*_k,\tau) + B(\tau) u^*_k,
\end{equation}
where the sequence $(x^*_{k+1},x^*_k,u^*_k)$ is a solution to \eqref{eq:sys_dist_nom} (note that in the following sections, employing an MPC negates the need to solve for $u^*_k$).

\subsection{MPC using Multiple Timescale Models}
For processes described by \eqref{eq:sys_dist_nom} to ``optimally'' track time-varying operating targets, i.e. $x_k\rightarrow x_k^*$, the incorporation of a prediction based control strategy, such as MPC, is naturally befitting. In order to design an MPC, the optimisation horizon must be long enough to cover the dominant slow dynamics (with long time constant), with a time interval small enough to capture the fast dynamics. In general, this is relatively impractical due to computational limitations. One approach is to incorporate multiple timescale models, which allows the optimisation horizon to capture the slower process dynamics without significantly increasing computational complexity. The controller, $u_k$, can then be computed to additionally optimise system economy via the minimisation of an arbitrary cost function, $\ell$, i.e., an MPC solves the following multi-timescale (multi-model) optimisation problem,
\begin{equation}
\begin{aligned}\label{eq:EMPCini}
\quad \min_{\hat{u}} & \sum_{i=0}^{N} \ell (\hat{x}_i,\hat{u}_i,\tau_i), \\
\text{s.t.} \quad & \hat{x}_0 = x_k,  \\
\quad & \hat{x}_{i+1} = f(\hat{x}_i,\tau_i) + B(\tau_i)\hat{u}_i, \\
& \hat{u}_i \in \mathcal{U}, \quad \hat{x}_i \in \mathcal{X},
\end{aligned}
\end{equation}
where $N$ is the prediction horizon, $\hat{x}_i$, and $\hat{u}_i$ are the respective $i$-th step state and control predictions, computed using the time interval $\tau_i$ and corresponding modelled functions $f(\tau_i)$ and $B(\tau_i)$. $x_k$ is the current state measurement, $\mathcal{U},\mathcal{X}$ are the sets of permissible inputs and states respectively.

The time interval, $\tau_i \in \mathbb{R}^+$, is varied after $N_i$ steps during the optimisation horizon according to the process dynamics, to account for the multiple timescale nature of the process. This causes the discrete-time optimisation model of the process to vary (non-uniformly) with the time interval. For brevity, we consider $p$ possible time intervals and note that the time intervals $\tau_i$ in the optimisation horizon are unrelated to the physical sampling period (e.g., denoted by $\tau_\Delta$) between which process measurements are sampled and control actions are applied, which remains constant (i.e., equal to $\tau_\Delta$).

Assuming feasibility of the optimisation problem~\eqref{eq:EMPCini}, the optimal input trajectory, $\mathbf{\hat{u}}^{opt} = (\hat{u}^{opt}_{0}, \hat{u}^{opt}_{1}, \cdots, \hat{u}^{opt}_{N-1}) \in \mathcal{U}^{N}$, can be computed. The MPC is then implemented in receding horizon fashion, by applying the first control action $\hat{u}^{opt}_{0}$ to system \eqref{eq:sys_dist_nom} until the next step, $k+1$, i.e.,
\begin{equation}
\label{eq:EMPCu}
u(k) = \hat{u}^{opt}_{0} : [k,k+1).
\end{equation}

There are two key obstacles in finding solutions to the optimisation problem in \eqref{eq:EMPCini}:
\begin{enumerate}
    \item[i.] Assessing stability over non-uniform time intervals.
    \item[ii.] Ensuring convergence of solutions to time-varying reference trajectories under arbitrary cost functions.
\end{enumerate}
Traditional methods to ensure the stability of MPC schemes involve using the cost function as a Lyapunov function \cite{RaM09}. This approach relies on models with uniform time intervals, and hence a uniformly spaced optimisation horizon to show that the cost is decreasing. This problem is further compounded when arbitrary cost functions are considered, whereby the cost function can be utilised to represent the process economy. Optimal solutions can therefore be associated with state values that deviate from target operating points. Since arbitrary cost functions are considered, a stability condition is required \cite{ellis2017economic}. As the reference trajectory is considered time-varying, stability guarantees for an optimal controller require a condition that is also reference-independent (incremental stability).
As a consequence, the MPC \eqref{eq:EMPCini} requires a contraction based constraint or equivalent, as will be demonstrated in the following sections.

\subsection{Objective and Approach}
To achieve the objective of state trajectory tracking, i.e., $x_k \to x^*_k$, whilst optimising system economy, for multi-timescale nonlinear processes, we require an MPC with conditions that can ensure convergency along the full prediction horizon under each timescale. To this effect, we propose to:
\begin{enumerate}
    \item[i.] Determine timescale dependent tractable stabilising controller conditions using contraction metrics
    \item[ii.] Impose stability conditions as constraints on an MPC
\end{enumerate}
In this way, tractable offline stability conditions can be constructed and imposed online on an MPC, which has the capacity for incorporating multiple timescale models whilst simultaneously optimising an economic cost and ensuring tracking of operating targets.

\section{Contraction-constrained MPC using Multiple Timescales}\label{sec:ctr}
In the following section, the contraction theory framework is leveraged to derive stability conditions to ensure closed-loop objective tracking of \eqref{eq:sys_dist_nom} over specific timescales. These conditions are transformed into tractable synthesis conditions and eventually used to characterise the property of stabilising controllers, forming the necessary MPC stability constraints.

\subsection{Timescale Dependent Control Contraction Metrics for Multiple Timescale Nonlinear Systems}
Since varying the time interval, $\tau$, results in multiple, timescale specific models of the process, we rewrite the dynamics in \eqref{eq:sys_dist_nom} using the following convenient state model
\begin{equation}\label{eq:sys_nom_taus}
    x_{\tau,k+1} = f_\tau(x_{\tau,k}) + B_\tau u_{\tau,k},
\end{equation}
and analogously for the reference model
\begin{equation}\label{eq:sys_ref_taus}
    x^*_{\tau,k+1} = f_\tau(x^*_{\tau,k}) + B_\tau u^*_{\tau,k}.
\end{equation}
Naturally, the corresponding differential dynamics for \eqref{eq:sys_nom_taus} are also timescale specific, described by
\begin{equation}\label{eq:dsys_nom_taus}
    \delta_{x_{\tau,k+1}} = A_\tau\delta_{x_{\tau,k}} + B_\tau\delta_{u_{\tau,k}},
\end{equation}
where $A_\tau := \frac{\partial f_\tau}{\partial x}$, and we consider a suitable differential feedback controller
\begin{equation}\label{eq:dsys_fed_taus}
    \delta_{u_{\tau,k}} = K_\tau(x_{\tau,k}) \delta_{x_{\tau,k}}.
\end{equation}
We then have the following direct result from Theorem \ref{thm:pre ctr con}.
\begin{corollary}\label{cor:cc_taus}
For each time interval, $\tau$, given the discrete-time nonlinear system \eqref{eq:sys_nom_taus}, with differential dynamics \eqref{eq:dsys_nom_taus} and differential state-feedback controller \eqref{eq:dsys_fed_taus}, provided uniformly bounded metrics, $\underline{m}_\tau I \leq M_\tau(x_{\tau,k}) \leq \overline{m}_\tau I$, exist satisfying,
    \begin{equation}\label{eq:dist ccon}
        (A_\tau+B_\tau K_\tau)^\top M_{\tau,+}(A_\tau+B_\tau K\tau) - (1-\beta_\tau)M_\tau < 0,
    \end{equation}
    where $M_{\tau,+}=M_\tau(x_{\tau,k+1})$, the differential system is contracting with rate $0 < \beta_\tau \leq 1$ on discrete intervals of $\tau$. Furthermore, over such intervals, the nonlinear system is exponentially incrementally stabilisable.
\end{corollary}
\begin{proof}
Suppose there exists a smooth feedback controller, e.g., with general timescale dependent structure $u_{\tau,k} = k_\tau(x_{\tau,k})+\bar{u}_{\tau,k}$, with a corresponding differential form as in \eqref{eq:dsys_fed_taus}, i.e., $\frac{\partial k_\tau}{\partial x} = K_\tau$. Then, substituting \eqref{eq:dsys_fed_taus} into \eqref{eq:dsys_nom_taus} yields, the closed-loop differential system for \eqref{eq:sys_nom_taus} as
    \begin{equation}\label{eq:diffcl_d}
        \delta_{x_{\tau,k+1}} = \left(A_\tau(x_{\tau,k}) + B_\tau K_\tau(x_{\tau,k})\right)\delta_{x_{\tau,k}}.
    \end{equation}
Substituting the closed-loop differential dynamics \eqref{eq:diffcl_d} and timescale specific differential Lyapunov function,
\begin{equation}\label{equ:generalised distance}
    V_\tau(x_{\tau,k},\delta_{x_{\tau,k}}) = \delta_{x_{\tau,k}}^\top M_\tau(x_{\tau,k})\delta_{x_{\tau,k}},
\end{equation}
into the exponential discrete-time Lyapunov condition of
\begin{equation}
    \label{eqe:Vdotcondition}
    V_{\tau,k+1} - (1-\beta_\tau)V_{\tau,k} < 0,
    \end{equation}
provides the timescale specific contraction condition in \eqref{eq:dist ccon}. To see that satisfaction of this condition implies stabilisability over each time interval, note that \eqref{eq:dist ccon} is a timescale specific instance of \eqref{equ:pre ctr con} in Theorem \ref{thm:pre ctr con}.
\end{proof}
The condition in \eqref{eq:dist ccon} can then be transformed into a tractable synthesis condition by adapting the results of \cite{wei2021control}.
\begin{corollary}\label{cor:cc_tractable}
For each pair of matrix functions $(W_\tau(x_{\tau,k}),L_\tau(x_{\tau,k}))$ satisfying instances of the condition
\begin{equation}\label{eq:tractable_con}
 \begin{bmatrix}
            W_{\tau,+} & A_\tau W_\tau +B_\tau L_\tau \\
           (A_\tau W_\tau +B_\tau L_\tau)^\top  & (1-\beta_\tau) W_\tau
        \end{bmatrix} > 0,
\end{equation}
where $W_\tau = M_\tau^{-1}$, $L_\tau=K_\tau W_\tau$, $W_{\tau,+}=W_\tau(x_{\tau,k+1})$, $M_\tau$ is uniformly bounded as $\underline{m}_\tau I \leq M_\tau(x_{\tau,k}) \leq \overline{m} _\tau I$ and $\beta_\tau \in (0,1]$, $M_\tau$ is then a discrete-time-control contraction metric and $K_\tau$ is a contracting differential feedback gain for \eqref{eq:sys_nom_taus}--\eqref{eq:dsys_fed_taus} over intervals of $\tau$.
\end{corollary}
\begin{proof}
Define $W_\tau(x_{\tau,k}) := M_\tau^{-1}(x_{\tau,k})$, $W_\tau(x_{\tau,k+1}) := M_\tau^{-1}(x_{\tau,k+1})$, $L_\tau(x_{\tau,k}) := K_\tau(x_{\tau,k})W_\tau(x_{\tau,k})$, apply Schur's complement to \eqref{eq:dist ccon}, and left/right multiply by $\diag\{I,W_\tau\}$ (and its transpose) to obtain \eqref{eq:tractable_con}.
\end{proof}
\begin{remark}
Offline solution to \eqref{eq:tractable_con} can be  completed, e.g., via the methods in \cite{wei2021contractionsyntehsis,wei2021discrete}, which in turn implies problem feasibility (in the sense of stabilising controller existence), whereby one particular feasible controller (unconstrained) is given by the feedback control structure obtained by integrating \eqref{eq:dsys_fed_taus} along the geodesic $\gamma_\tau(x_{\tau,k},x^*_{\tau,k})$, i.e.,
\begin{equation}\label{equ:ctl int taus}
    u_{\tau,k} = u^*_{\tau,k} + \int_0^1 K_\tau(\gamma_\tau(s))\delta_{\gamma_{\tau,k}(s)} \,\,ds.
\end{equation}
To additionally ensure physical constraint satisfaction, a common approach is to search offline for less aggressive contraction rates, resulting in smaller control action magnitudes, and hence comply with the system constraints (see, e.g., \cite{mccloy2021differential}).
\end{remark}

\subsection{Contraction-based Stability Constraint}
Optimal solutions to the MPC problem \eqref{eq:EMPCini} can be associated with state values which deviate from target operating points. As a consequence, the MPC requires a contraction based constraint or equivalent to ensure reference tracking, as will be constructed in the following.
\begin{corollary}\label{cor:cc_predH}
For \eqref{eq:sys_nom_taus}, the contraction constraint
\begin{equation}\label{eq:shrinking_Rdist}
\begin{split}
&r(x_{\tau,k},x^*_{\tau,k},u_{\tau,k},\beta_\tau):= d\left(\gamma_\tau(x_{\tau,k+1},x^*_{\tau,k+1})\right) \\&\qquad - \left(1-\beta_\tau) \right)^{1/2} d\left(\gamma_\tau(x_{\tau,k},x^*_{\tau,k})\right) \leq 0,
\end{split}
\end{equation}
where the geodesic $\gamma_\tau$ is computed using \eqref{equ:geodesic} with respect to the timescale specific metric $M_\tau$, defines the set of contracting tracking controllers which ensure the system state, $x_{\tau,k}$, is exponentially convergent to the feasible trajectory, $x^*_{\tau,k}$.
\end{corollary}
\begin{proof}
Integrating \eqref{eqe:Vdotcondition} with \eqref{equ:generalised distance} along the geodesic $\gamma_\tau(x_{\tau,k}, x_{\tau,k}^*)$ yields (with a slight abuse of notation)
\begin{equation}\label{equ:energy bound2}
\begin{aligned}
        &\int_0^1 \delta_{\gamma_{\tau,k+1}(s)}^\top M_\tau(\gamma_{\tau,k+1}(s))\delta_{\gamma_{\tau,k+1}(s)} \,\,ds \\&\qquad\leq \int_0^1 (1-\beta_\tau) \delta_{\gamma_{\tau,k}(s)}^\top M_\tau(\gamma_{\tau,k}(s))\delta_{\gamma_{\tau,k}(s)} \,\,ds.
\end{aligned}
\end{equation}
Using \eqref{equ:Riemannian distance and energy}, \eqref{equ:geodesic}, \eqref{equ:energy bound2}, the property $d(\gamma)^2 \leq d(c)^2$, and taking square roots gives the condition in \eqref{eq:shrinking_Rdist}.
\end{proof}
\begin{remark}
The contraction constraint in \eqref{eq:shrinking_Rdist}, has a reference-independent structure and is a nonlinear function affine in the control action, whereby the information required for the computation of $d(\gamma_\tau)$ is available from the current state measurement and knowledge of the system model and reference dynamics. In addition, these conditions define a set (or the conditional property) of controllers which guarantee the system is contracting (under Corollary \ref{cor:cc_taus}), whereby \eqref{equ:ctl int taus} provides one particular solution.
\end{remark}

\subsection{Constructing the MPC}
The multistage optimisation cost in \eqref{eq:EMPCini} can be reparameterised using the more convenient form in \eqref{eq:sys_nom_taus} (cf. \eqref{eq:sys_dist_nom}) as
\begin{equation}\label{eq:multi_cost}
     \sum_{i=0}^{N} \ell (\hat{x}_i,\hat{u}_i,\tau_i) = \sum^p_{i=1}  \sum^{N_{i}}_{j = 0} \ell(\hat{x}_{i,j},\hat{u}_{i,j},\tau_i),
\end{equation}
whereby the prediction horizon of length $N$ has been broken into (non necessarily uniform) segments based on the desired number of predictive steps using each time interval, i.e.,
\begin{equation}\label{eq:tausteps}
    N = \sum_{i=1}^p N_i, \qquad \hat{k} = \sum_{i=1}^p N_i \cdot \frac{\tau_i}{\tau_\Delta},
\end{equation}
whereby $\hat{k}$ denotes the total number of predicted discrete time-steps forward for the actual system in \eqref{eq:sys_dist_nom}. For example, given that the control actions for \eqref{eq:sys_dist_nom} are implemented on intervals of $\tau_\Delta$, suppose that $p=2$, $\tau_1 = \tau_\Delta$, $\tau_2 = 5\cdot\tau_\Delta$ and $N_1=N_2=1$. Then, in this scenario, the total number of predicted discrete time-steps forward for \eqref{eq:sys_dist_nom} would be $\hat{k} = 6$. Relative to a uniform timescale model ($\hat{k} = N = 2$) this approach offers an increased number of prediction steps without introducing additional complexity (or, by choosing an equivalent $\hat{k}$, at reduced computational complexity).

Each subsequent predictive model (corresponding to a specific $\tau_i$) used in \eqref{eq:multi_cost} is initialised with the last predicted state of the previous model, i.e.,
\begin{equation}\label{eq:multi_stage_state}
    \hat{x}_{1,0} = x_k, \qquad \hat{x}_{i,0} = \hat{x}_{i-1,N_j}, \quad i=2,\cdots,p
\end{equation}
and hence each sequence of predicted control inputs is defined as the concatenation along each segment, i.e.,
\begin{equation}\label{eq:multi_stage_control}
\begin{split}
\mathbf{\hat{u}} &= (\hat{u}_0,\cdots,\hat{u}_{N-1}) \\ &=\left((\hat{u}_{1,0},\dots,\hat{u}_{1,N_1-1}),(\hat{u}_{2,0},\dots,\hat{u}_{2,N_2-1}),\right.\\&\left.\qquad\cdots,(\hat{u}_{p,0},\dots,\hat{u}_{p,N_p-1}) \right)
\end{split}
\end{equation}

From \eqref{eq:sys_nom_taus}, \eqref{eq:shrinking_Rdist}, \eqref{eq:multi_cost}--\eqref{eq:multi_stage_control} we can construct the contraction-constrained multi-timescale MPC optimisation problem:
\begin{equation}
\begin{aligned}\label{eq:EMPC}
\quad \min_{\mathbf{\hat{u}}} & \sum^p_{i=1}  \sum^{N_{i}}_{j = 0} \ell(\hat{x}_{i,j},\hat{u}_{i,j},\tau_i), \\
\text{s.t.} \quad & \hat{x}_{1,0} = x_k, \quad \hat{x}_{i,0} = \hat{x}_{i-1,N_j},  \\
\quad & \hat{x}_{i,j+1} = f_i(\hat{x}_{i,j}) + B_i\hat{u}_{i,j}, \\
& \hat{u}_i \in \mathcal{U}, \quad \hat{x}_i \in \mathcal{X}, \\ &r(\hat{x}_{i,j},x^*_{i,j},\hat{u}_{i,j},\beta_{\tau_i}) \leq 0
\end{aligned}
\end{equation}

For the MPC in \eqref{eq:EMPC}, we then have the following.
\begin{theorem}\label{thm:mainmpc}
    An MPC with solution to \eqref{eq:EMPC} achieves economic optimisation simultaneously with reference tracking for multiple timescale nonlinear systems of the form in \eqref{eq:sys_dist_nom}.
\end{theorem}
\begin{proof}
By considering a discrete number of time intervals, $\tau_i$, $i=1,\cdots,p$, \eqref{eq:sys_dist_nom} can be described (potentially exactly) by $p$ models of the form in \eqref{eq:sys_nom_taus}. For each $i$-th model, the pair $(M_{\tau_i},K_{\tau_i})$ are found satisfying \eqref{eq:dist ccon} (e.g., via solution to \eqref{eq:tractable_con}).
Reference convergence over each horizon $N_i$ is ensured since economic cost minimising solutions to \eqref{eq:EMPC} are constrained to the set of control actions which possess the contracting property under Corollary \ref{cor:cc_predH}. Since each subsequent predictive model (used over the horizon of length $N_i$) is initialised with the last predicted state of the previous model via \eqref{eq:multi_stage_state}, the resulting sequence of predicted states (governed by a predicted control sequence of the form in \eqref{eq:multi_stage_control}) reparameterised from \eqref{eq:sys_nom_taus} and concatenated along each segment
\begin{equation}\label{eq:multi_stage_predicted_state}
\begin{split}
\mathbf{\hat{x}} &=\left((\hat{x}_{1,0},\dots,\hat{x}_{1,N_1}),(\hat{x}_{2,0},\dots,\hat{x}_{2,N_2}),\right.\\&\left.\qquad\cdots,(\hat{x}_{p,0},\dots,\hat{x}_{p,N_p}) \right)\\
&= (\hat{x}_0,\cdots,\hat{x}_{N})
\end{split}
\end{equation}
converges to the concatenated reference sequence reparameterised from \eqref{eq:sys_ref_taus}
\begin{equation}\label{eq:multi_stage_reference}
\begin{split}
\mathbf{x}^* &= (x^*_0,\cdots,x^*_{N})\\
&=\left((x^*_{1,0},\dots,x^*_{1,N_1}),(x^*_{2,0},\dots,x^*_{2,N_2}),\right.\\&\left.\qquad\cdots,(x^*_{p,0},\dots,x^*_{p,N_p}) \right),
\end{split}
\end{equation}
i.e., $d(\hat{x},x^*) \to 0$ and hence $\hat{x} \to x^*$. Initialising the predicted state sequence using the current state measurement, i.e., $\hat{x}_{1,0} = x_k$, and implementing the controller in receding horizon fashion \eqref{eq:EMPCu}, ensures $x_k \to x^*_k$ in the absence of prediction errors.
\end{proof}
\begin{remark}
Prediction accuracy depends on choosing an adequate number of timescale specific models (many processes are naturally described by an infinite number of timescales). For a finite number of intervals and hence models, the effects of prediction errors on the reference tracking error can be attenuated by using the differential dissipativity techniques developed in \cite{wei2021contractionsyntehsis} and imposing additional requirements during metric synthesis (via a truncated incremental $\mathcal{L}_2$ gain).
\end{remark}

\subsection{Discussion}

The proposed MPC approach achieves economic optimisation with stability for multiple timescale nonlinear systems. In particular, this carries several advantages:
\begin{itemize}
    \item Economic reference tracking -- the optimisation problem permits arbitrary economic cost functions, with reference tracking ensured by contraction constraints.
    \item Reference flexibility -- the contraction (stability) constraint is structurally reference-independent, constructed offline and updated online as new state measurements and reference information is available.
    \item Reduced computational burden -- with respect to an equivalent single timescale approach, longer prediction horizons can be considered by using longer time intervals (and incorporating appropriate slower dynamic models).
    \item Scalability via \textit{modular} stability constraint construction -- contraction constraints can be systematically added to the optimisation problem as additional time intervals (and hence models) are considered, without increasing offline synthesis or online computational requirements for an equivalent prediction length (each constraint is formed and imposed independently).
\end{itemize}

In future work, natural extensions can be made to the existing framework to simultaneously incorporate reduced order models (e.g., by adjusting model granularity \cite{findeisen2016exploiting}). Feasibility analysis, extension to a more general class of nonlinear systems and explicit consideration for prediction / timescale discretisation errors are also part of future studies.

\section{Illustrative Example}\label{sec:exa}
In the following numerical example, we consider the coupled-tank process described in~\cite{mccloy2014set}, consisting of two cylindrical tanks positioned one directly above the other. Liquid flows from an upper primary tank through a lower pilot tank to a reservoir, and a pump is used to push liquid from the reservoir back up to the upper tank. The objective is to minimise the pump energy input whilst achieving a target liquid level in the primary tank (matched by the accessible secondary tank). The discrete-time dynamics of the liquid levels $x_1$ (upper tank) and $x_2$ (lower tank) are described as
\begin{equation} \label{eq:coupledTanks}
\begin{aligned}
x_{1,k+1} & = x_{1,k} -\tau\frac{\sigma_1}{\alpha_1} \sqrt{2g x_{1,k}} + \tau\frac{k_p}{\alpha_1} u_k,\\
x_{2,k+1} & = x_{2,k} +\tau\frac{\sigma_2}{\alpha_2} \sqrt{2g x_{1,k}} - \tau\frac{\sigma_2}{\alpha_2}
\sqrt{2g x_{2,k}},
\end{aligned}
\end{equation}
where $u_k$ is the voltage applied to the pump, and the parameters are as follows:
$\alpha_1=155$cm$^2$ and $\alpha_2=15.5$cm$^2$ are the cross-sectional areas of tanks 1 and 2; $\sigma_1=\sigma_2= 0.178$cm$^2$ are the cross-sectional areas of each tank's outflow
orifice; $k_p=13.2$cm$^3$/Vs is the gain of the pump; $g=980$cm/s$^2$; and $\tau_\Delta = 1s$ is the natural sampling period of the system. The liquid levels are constrained to the range $x_1,x_2 \in [0,10]$cm, and the pump voltage is limited to $u \in [0,10]$V.

To capture the relatively fast and slow dynamics of each tank, two timescales ($p=2$ models), $\tau_1 = 1s$ and $\tau_2 = 10s$, are considered. Control contraction metrics are then synthesised via the following Sum-of-Squares programming formulation (see \cite{wei2021control} for more detail),
\begin{equation}\label{min:sos}
    \begin{aligned}
        &\underset{l_{\tau,c},w_{\tau,c},r}{\min} r\\
        &s.t. \ \phi^\top  \Omega_\tau \phi -rI \in \Sigma(x_k,u_k,\phi), ~r \geq 0,
    \end{aligned}
\end{equation}
where $\Omega_\tau$ is formed by the left hand side of \eqref{eq:tractable_con}, $\Sigma(x_k,u_k,\phi)$ is a Sum-of-squares function of the arguments $x_k,u_k$ and $\phi$, $\phi$ is a vector of monomials, $l_{\tau,c},w_{\tau,c}$ are the polynomial coefficients of elements of $W_\tau$ and $L_\tau$ in \eqref{eq:tractable_con}, and $r$ is a small positive constant. Solution pairs $(W_1,L_1)$ and $(W_2,L_2)$ to the Sum-of-Squares programming problem in \eqref{min:sos} satisfy \eqref{eq:tractable_con} and consequently the contraction condition in \eqref{eq:dist ccon} using the definitions $W=M^{-1},L=KW$. Hence, the pairs $(M_1,K_1)$ and $(M_2,K_2)$ were obtained from solutions to \eqref{min:sos}, for the two respective time scale systems, $\tau = 1$ and $\tau = 10$, with contraction rates $\beta_1=\beta_2 = 0.3$.

The MPC is then formulated as in \eqref{eq:EMPC} with two contraction metrics $M_1,M_2$ and two prediction horizons $N_1 = 1$, $N_2 = 3$ (corresponding to the two time intervals $\tau_1$ and $\tau_2$). The resulting total number of predicted time steps was computed via \eqref{eq:tausteps} as $\hat{k} = N_1 \cdot \frac{\tau_1}{\tau_\Delta} + N_2 \cdot \frac{\tau_2}{\tau_\Delta} = 31$. An economic cost function was designed to reflect the control effort, i.e.,
$\ell = ||u_k||^2$,
subject to the time-varying operating targets: $(x_1^*,x_2^*)$=$(2.5,2.5)$, $(5,5)$ and $(7.5,7.5)$ on the respective intervals of  $[0,40)$, $[40,80)$, $[80,120]s$. Constraint feasibility of \eqref{eq:shrinking_Rdist}, under the contraction rates $\beta_1,\beta_2$, was verified via constraint sampling of $u \in \mathcal{U}$ using $(M_1,K_1)$, $(M_2,K_2)$, the controller \eqref{equ:ctl int taus} and model \eqref{eq:coupledTanks} for $x\in \mathcal{X}$.

The closed-loop system was then simulated with the MPC \eqref{eq:EMPC} computed for each $\tau_\Delta=1s$ interval. The average computation time, per simulated step, for the proposed MPC, 'MPC($\tau_1,\tau_2$)', is compared with an equivalent single timescale MPC, 'MPC($\tau_\Delta$)', which, as shown in Table \ref{tbl:compt}, indicates that the proposed method requires significantly less computational effort for an equivalent number of prediction steps on average (in fact, the single timescale method is infeasible under the specified control period of $1s$).
\begin{table}[h]\caption{Computation Time Comparison.}\label{tbl:compt}
\centering
\begin{tabular}{c|c|c}
    Method & Predicted Steps & Ave. Comp. Time \\
    \hline
    MPC($\tau_1,\tau_2$) & $\hat{k} = 31$ & $0.425s$\\
    MPC($\tau_\Delta$) & $\hat{k} = 31$ & $3.106s$\\
\end{tabular}
\end{table}

Finally, the simulation results shown in Fig. \ref{fig:sim}, illustrate that the proposed contraction constrained multi-timescale MPC is capable of driving the system to the operational targets whilst minimising the input energy as per Theorem \ref{thm:mainmpc}.

\begin{figure}\label{fig:sim}
\centering
        \includegraphics[width=\linewidth]{sim.pdf}
        \caption{Simulation results for control of \eqref{eq:coupledTanks} using \eqref{eq:EMPC}.}
            \vspace{-2.5mm}
\end{figure}

\section{Conclusion} \label{sec:conclusion}
A contraction constrained MPC approach was developed for discrete-time nonlinear processes with multiple timescale dynamics. Through the differential system framework, stability conditions were derived to ensure convergence of the resulting closed-loop to feasible time-varying references. These conditions were also generalised to characterise a set of stabilising controllers and imposed as stability constraints on an MPC. The MPC would then search for a cost minimising controller amongst those stabilising controllers whilst utilising multiple timescale models. The result was an MPC capable of ensuring convergency of multiple timescale nonlinear processes. The overall approach was effectively demonstrated via simulation.

\bibliographystyle{IEEEtran}


\section{Motivation}
In the 1930s Anderson and Neddermeyer noted the observation of unexpected charged long-lived particles (LLPs) `less massive than protons but more penetrating than electrons' \cite{Anderson:1936zz,Neddermeyer:1937md}.  This unexpected discovery of new physics led Rabi to ask `who ordered that?', a question unanswered to this day.   We may infer two lessons for modern particle physics.  The first is that it would not be without precedent if new physics emerged at high energies in the form of charged LLPs.  A second is that a new physics discovery need not conform to any theoretical preconceptions nor answer any particular outstanding theoretical question.  We may, again, be left asking `who ordered that?' for decades to come.

It goes without saying that we should keep our eyes open for new LLPs, wherever we can.  Indeed, there are a wide variety of LLP searches undertaken across a range of energy scales, reflecting the wide range of possibilities that could give rise to them (see e.g.~\cite{Fairbairn:2006gg}).  At the TeV-scale the LHC provides numerous opportunities for LLP discovery.  One is raised by searching for large ionisation energy loss gradients ($dE/dx$) in the tracker of the ATLAS and CMS detectors, see e.g.~\cite{CMS:2011arq,ATLAS:2011ghv,ATLAS:2012urj,CMS:2012wcg,ATLAS:2013cab,CMS:2013czn,ATLAS:2015wsk,ATLAS:2016tbt,CMS:2016kce}. This allows to efficiently distinguish possible new-physics signals from the SM background, which at large $dE/dx$ is small to none, since the known particles with long lifetime are produced relativistically at the LHC.

This can be understood quantitatively by recalling that, for ionising particles with electric charge $Q$ (in multiples of the electron charge) and speed $v$ (with $\beta = v/c$ and $\gamma = 1/\sqrt{1- \beta^2}$), the mean energy loss per distance travelled is given by the Bethe-Bloch relation which, up to density effect corrections, takes the form
\begin{eqnarray}
-\left\langle \frac{dE}{dx}\right\rangle = 4\pi m_e n_e r_e^2  Q^2 \left(-1 +\frac{2}{\beta^2} \ln \frac{\beta\gamma}{I_e} \right)
 \, . \label{eq:BetheBloch}
\end{eqnarray}
Here $m_e$, $n_e$ and $r_e$ are the electron mass, number density in the medium and classical radius, while $I_e$ is a coefficient related to the mean excitation energy of the medium and the maximum energy transfer in a single electron collision. To adhere with conventions commonly used in high-energy physics~\cite{PDG}, from now on we will refer to $dE/dx$ as the mass stopping power, which is defined as $-\langle dE/dx \rangle /\rho$, where $\rho$ is the mass density of the material. The minus sign in the definition ensures that the mass stopping power $dE/dx$ is positive.

For a particle of unit charge, a large $dE/dx$ signal is expected if the LLP, leaving the ionising track in the Inner Detector, has a relatively small $\beta \lesssim 0.7$. This is the typical working assumption of the collaborations in analysing their data (with exceptions\footnote{These works consider direct production of multi-charged particles. These particles, being relatively slow, would give rise to an ionisation signal much greater than considered here, and a different analysis strategy is typically required.}~\cite{ATLAS:2013cab,CMS:2013czn,CMS:2016kce}). However, a large $dE/dx$ signal is possible also in the alternative assumption of $\emph{fast}$ ionising particles with $Q>1$.  To this end, in this work we develop a strategy to re-interpret the $Q=1$ ATLAS analyses in terms of this hypothesis.

We do not hesitate to add that we are motivated by  the recent  ATLAS announcement of an
excess in the large $dE/dx>2.4 \, \mathrm{MeV g^{-1} cm^2}$ data~\cite{ATLAS:2022pib,ATLAStalk} given by $7$ events in a region with known background of $0.7 \pm 0.4$ events, corresponding to a local (global) significance of 3.6$\sigma$ (3.3$\sigma$), if interpreted as due to metastable gluinos.  Only time will tell whether this excess of events is due to new physics or not. However, since these events lie in a signal-dominated region and the statistical significance will be rapidly tested with new data, it is timely to assess what new physics such events could correspond to.

The ATLAS collaboration analyses the excess as being due to relatively slow ($\beta \approx 0.5$) charge $Q=1$ particles with mass in the $1.0$--$2.5 \text{ TeV}$ range, considering benchmark models such as $R$-hadrons formed by gluinos in Split-Supersymmetry and metastable charginos or sleptons. However, the time-of-flight determination of $\beta$ instead indicates that all excess events have $\beta \simeq 1$ at $95\%$ confidence level, with uncertainties of about $0.1$ at  2$\sigma$~\cite{ATLAS:2022pib}.  This measurement disfavours any interpretation of the ATLAS excess in terms of conventional LLPs models.  In this work we propose a new scenario of boosted LLPs, which is consistent with observations.  Specifically, we find that the $dE/dx$ excess and the time-of-flight measurements can be explained by boosted ($\beta \approx 1$) particles in the TeV range and with larger electric charge, here $Q=2$ for concreteness, see Fig.~\ref{fig:beta}.  Such boosted states may be produced from the decay of a parent resonance with mass in the range $M_P \approx 4$--$6 \text{ TeV}$, opening the door to a new class of new physics scenarios that may be discoverable at the LHC.

\begin{figure}[t]
\centering
\includegraphics[width=25em]{beta} \caption{Distribution of the $dE/dx$ signal as function of $\beta$ for unit and doubly charged particles. The continuous lines denote the most probable values \eqref{eq:dEdxMPV}.  The simulated points are spread around it by the Crystal Ball distribution given in Appendix~\ref{app:calibration}. The grey bands indicate the regions of the ATLAS excess events, in particular their ionisation  $dE/dx$  and $\beta \approx 1$, as suggested by time-of-flight measurements.
\label{fig:beta}}
\end{figure}

\section{The {\boldmath$dE/dx$} signal for {\boldmath$Q>1$}}\label{sec:dEdx}
We begin by summarising the analysis strategy of the ATLAS collaboration~\cite{ATLAS:2012urj,ATLAS:2015wsk,ATLAS:2016tbt,ATLAS:2022pib}. The most probable value (MPV) of the mass stopping power is calibrated by light SM particles ($\pi^\pm, K^\pm, p$) with $|Q|=1$ \cite{Liang:2014via} and fit using a three-parameter ($c_0,c_1,c_2$) functional form of a phenomenological Bethe-Bloch-like relation
\begin{eqnarray}
\frac{dE}{dx}\bigg|_{\rm MPV}\!\!\!\!\!\!\!\!\!\!(\beta \gamma) = \frac{1+(\beta\gamma)^2}{\beta\gamma} \Big[ c_0 + c_1 \log_{10}(\beta\gamma) + c_2 \log_{10}^2 (\beta\gamma) \Big] \;.  \label{eq:dEdxMPV}
\end{eqnarray}
For each event with measured momentum $p_{\rm m}$ and ionisation energy loss $dE/dx|_{\rm m}$, an effective mass $m_{dE/dx}$ is obtained by inverting \eqref{eq:dEdxMPV}
\begin{eqnarray}
\frac{dE}{dx}\bigg|_{\rm MPV}\!\!\!\!\!\!\!\!\!\!(p_{\rm m}/m_{dE/dx}) = \frac{dE}{dx}\bigg|_{\rm m} \;.
\end{eqnarray}
The ATLAS excess is in the region
\begin{eqnarray}
1.0 \text{ TeV} \lesssim m_{dE/dx} \lesssim 2.5 \text{ TeV}~~.
\end{eqnarray}
For unit charge particles the effective mass is expected to be close to the physical mass of the particle causing the signal, $m_{dE/dx} \approx m$. However we find that two factors mainly spread $m_{dE/dx}$ around the physical value: the $dE/dx$ distribution around the MPV has a significant width, see Fig.~\ref{fig:beta}, and the momentum resolution in the Inner Detector deteriorates at large values of $p$. This latter factor dominates at the values of interest $p \gtrsim \text{ TeV}$.

The unit-charge hypothesis is implicit in the ATLAS analysis, so a direct recast to $Q=2$ is not possible. Thus to analyse the possibility of different charges we develop a different strategy.

\begin{figure}
\centering
\subfloat[]{\includegraphics[width=25em]{benchmark_m}
\label{fig:benchmark_m}}  \vspace{1em} \\
\subfloat[]{\includegraphics[width=25em]{benchmark_pT}
\label{fig:benchmark_pT}}  \vspace{1em} \\
\subfloat[]{\includegraphics[width=25em]{benchmark_dEdx}
\label{benchmark_dEdx}}
\caption{{\rm (a)} The effective-mass histogram for the observed data~\cite{ATLAS:2022pib} (black dots), background distribution (red line, taken from~\cite{ATLAS:2022pib}) and background plus $Q=2$ signal model (green line). {\rm (b)} Same for $p_T$.  {\rm (c)} Same for $dE/dx$, with the background extracted from~\cite{ATLAStalk}.} \label{fig:trio}
\end{figure}

We simulate a set of events with true momentum $p$ and ionisation $dE/dx$ given by the distribution around the Bethe-Bloch curve with $Q=2$. This is obtained by the phenomenological relation~\eqref{eq:dEdxMPV}, multiplied by the nominal charge factor $Q^2$.\footnote{This is sufficient at our level of accuracy. Notice that a direct calibration of the Bethe-Bloch curve for $Q>1$ is not available; one could expect an experimental charge factor $Q_{\rm eff}^2$, with $Q_{\rm eff} \approx 2$. We checked that changing $Q_{\rm eff}$ from its nominal value 2 by $\mathcal{O}(10\%)$ does not affect our results significantly.} Then, we take into account that the momentum reconstruction algorithm from the tracker implicitly assumes $Q=1$, since the transverse momentum is obtained as $p_T = Q B \rho$, with $\rho$ being the radius of curvature of the track and $B$ the magnetic field. Therefore, the momentum assigned by the experiments to the event is $p_{\rm rec} = p/2$, half of the real one. In our simulations $p_{\rm rec}$ is then spread by the detector resolution curve reported in~\cite{Benekos:2003xra}.

Finally, the $Q=1$ ATLAS algorithm is used to generate the $m_{dE/dx}$ histograms, as described at the beginning of this section. Because of the mismatch of charge the effective $m_{dE/dx}$ does not peak around the physical mass $m$ of the particle, but this is not a problem for the analysis, since we find that acceptable signal models are obtained by this procedure, which effectively allows us to interpret possible signals in ATLAS data in terms of $Q=2$ particles.

\section{Results for the ATLAS excess}\label{sec:results}
We now study the excess recently reported by the ATLAS collaboration~\cite{ATLAS:2022pib,ATLAStalk}. To this end, we obtain the Bethe-Bloch curve~\eqref{eq:dEdxMPV} from the calibration data extracted from~\cite{ATLAS:2022pib}, with results reported in Appendix~\ref{app:calibration}. The distribution around the MPV is fitted from the $\pi^\pm, K^\pm, p$ data in~\cite{ATLAS:2022pib} by a one-sided Crystal Ball function. We validate our procedure by reproducing the $Q=1$ signal models considered by the collaboration (see Appendix~\ref{app:validation}).

We consider a model given by a parent resonance with mass $M_P$, decaying into two metastable $Q=2$ daughter particles with mass $m_d$, which give rise to the observed excess. Possible explanations for the microscopic origin of this scenario are discussed in the next section. The decay of the parent particle into two daughters is approximated as isotropic in the rest frame and gets longitudinally boosted in the lab frame, working at leading order\footnote{However, we find that the effect of the boost is small, being subdominant compared to the detector momentum resolution.}. We then generate histograms by applying the procedure described in the previous section, with cuts $dE/dx>2.4 \, \mathrm{MeV g^{-1} cm^2}$, $p_T > 120 \text{ GeV}$, and perform a profile-likelihood analysis with Poissonian likelihoods.

\begin{figure}[t]
\centering
\includegraphics[width=20em]{fit} \caption{Profile-likelihood fit of the $m_{dE/dx}$, $p_T$ and $dE/dx$ distributions. Contours show the 1$\sigma$ and 2$\sigma$ preferred regions. The star indicates the best-fit point located at $M_P = 5.2~ {\rm TeV}$ and $m_d = 650~{\rm GeV}$. The dashed lines denote the corresponding values of $\beta$ for a decay at rest of the parent resonance.
\label{fig:fit}}
\end{figure}

Our model provides an excellent fit to the excess, as shown in Fig.~\ref{fig:benchmark_m}, with an approximate local significance of about $ 4\sigma$ relative to a background-only scenario\footnote{Following ATLAS, local significances are estimated by fitting the $m_{dE/dx}$ histogram only, using Wilks' theorem.}.
Most importantly, in contrast to the benchmark models considered in~\cite{ATLAS:2022pib}, here the new-physics events have $\beta \approx 1$ by construction, due to the boost resulting from heavy parent decays. We have also checked that the excess is reproduced in the $p_T$ histogram (which is affected by the large $dE/dx$ cut), see Fig.~\ref{fig:benchmark_pT}, and that the $dE/dx$ distribution agrees with data (Fig.~\ref{benchmark_dEdx}).

Finally, we have studied the parameter space that can explain the excess, performing various parameter fits. The results are shown in Fig.~\ref{fig:fit}, obtained fitting all three $m_{dE/dx}$, $p_T$ and $dE/dx$ histograms. Since these are not independent, the confidence intervals are obtained by toy pseudo-experiments, as described in~Appendix~\ref{app:fits}.
All in all, the excess can be explained by charge-two particles in the mass range between hundreds of GeV and the TeV, produced boosted by a parent particle with $M_P \gtrsim 3.5 \text{ TeV}$.
The $m_{dE/dx}$ distribution alone is well reproduced for $M_P \gtrsim 2 \text{ TeV}$ (see Appendix~\ref{app:fits}), but the inclusion of the $p_T$ information shifts the preferred range of $M_P$ to larger values. Assuming that the excess persists, the local significance of the best-fit point of our model would reach about $6.2\sigma$ at the end of Run 3 (with $460 \units{fb}^{-1}$ of data).

Figure \ref{fig:fit} describes only the fit to experimental data, but carries no information about the physical production mechanism of the parent particle. As discussed in the next section, realistic models can efficiently produce resonances at the LHC only up to about 6~TeV and this explains why we have limited the vertical axis of Fig.~\ref{fig:fit}. The best-fit point, indicated by a star, lies within the interesting physical mass region.

\section{Microscopic physics} \label{sec:models}
Having demonstrated that the kinematic pattern of heavy resonance production followed by decay to doubly-charged LLPs provides a unique and interesting phenomenological scenario for $dE/dx$ searches, as well as a candidate explanation for the recently observed excess, it naturally follows to briefly explore the microscopic physics that could underlie such a scenario.

Even without specifying the details of the microscopic model, we can derive some general properties of the resonance and its phenomenological consequences from the basic features of our physical setup. We are considering a colourless parent particle $P$ with mass $M_P$ and spin $J_P$, which can decay into a pair of doubly-charged long-lived daughter particles with branching ratio $B_d \equiv {\rm BR}(P\to dd^c)$. Since $P$ is colour singlet, it can be resonantly produced at the LHC only in the partonic channels $i=gg,q{\bar q}$. We assume that $P$ is coupled to at least one of these possible initial states, with corresponding branching ratio $B_i \equiv {\rm BR}(P\to i)$.

The parent and daughter particles have rather characteristic properties and cannot be immediately embedded in conventional new-physics scenarios. If real, they are likely part of a richer structure and accompanied by other new particles\footnote{Incidentally, ref.~\cite{Akhmedov:2021qmr} made the interesting remark that a doubly-charged long-lived (or stable) scalar particle could combine with light nuclei and catalyse their fusion, even under low-temperature and low-density conditions.}. In this context, it is not surprising that the charge-two daughter is the first particle to be discovered. Indeed, charge-one daughters would have been missed by ATLAS, since they produce a four-times smaller ionisation energy loss (because $dE/dx \propto Q^2$), while neutral daughters do not generate ionisation tracks.

The total cross section for parent resonant production in $pp$ collisions with centre-of-mass energy $\sqrt{s}$ is, in narrow width approximation ({\it i.e.} $\Gamma_P \ll M_P$, where $\Gamma_P$ is the $P$ total decay width),
\begin{equation}}% can be used as {equation} or {eqnarray
\sigma_P =\frac{2J_P+1}{s} \, \frac{\Gamma_P}{M_P} \sum_i C_i B_i ~,
\end{equation}
\begin{equation}}% can be used as {equation} or {eqnarray
C_{gg} =\frac{\pi^2}{8} \int_\tau^1 \frac{dx}{x} \, f_g(x) f_g(\tau x) ~,
~~~
C_{q{\bar q}} =\frac{4\pi^2}{9} \int_\tau^1 \frac{dx}{x} \, \left[ f_q(x) f_{\bar q}(\tau x) + f_{\bar q}(x) f_{q}(\tau x) \right]~,
\label{eqCi}
\end{equation}
where $\tau = M_P^2/s$ and the parton distribution functions $f_{g,q,{\bar q}}$ are evaluated at $Q^2 = M_P^2$. The values of $C_i$, for characteristic values of $M_P$, are tabulated in table~\ref{tab:C}.

\begin{table}[h]
\begin{center}
\begin{tabular}{|c||c|c|c|c|c|c|}
\hline
$M_P$  {[TeV]} & $C_{gg}$ & $C_{u{\bar u}}$ & $C_{d{\bar d}}$ & $C_{s{\bar s}}$ & $C_{c{\bar c}}$ & $C_{b{\bar b}}$\\
\hline
3 & $3.2 \times 10^{-1}$ & 1.4 & $4.7 \times 10^{-1}$ & $ 1.3\times 10^{-2}$  & $ 4.5\times 10^{-3}$ & $1.7 \times 10^{-3}$\\
4 & $ 2.1 \times 10^{-2}$ & $1.6 \times 10^{-1}$ & $3.2 \times 10^{-2}$ &  $7.3 \times 10^{-4}$  & $ 2.7\times 10^{-4}$ & $9.1 \times 10^{-5}$\\
5 & $1.6 \times 10^{-3}$ & $2.0 \times 10^{-2}$ & $1.8 \times 10^{-3}$ &$4.2 \times 10^{-5}$  & $1.7 \times 10^{-5}$ & $5.6 \times 10^{-6}$ \\
6 & $1.2 \times 10^{-4}$ & $2.1 \times 10^{-3}$ & $8.3 \times 10^{-5}$ & $2.2 \times 10^{-6}$  & $1.1 \times 10^{-6}$ & $3.4 \times 10^{-7}$\\
 \hline
 \end{tabular}
 \caption{The $C_i$ coefficients defined in eq.~\eqref{eqCi}, evaluated using the PDFs MSTW2008NLO~\cite{Martin:2009iq}, for $\sqrt{s}= 13$~TeV and for relevant values of the parent mass $M_P$.}
 \label{tab:C}
 \end{center}
\end{table}

The number of events with anomalous ionising tracks from fast-moving daughter particles is
\begin{equation}}% can be used as {equation} or {eqnarray
N_{\rm ev} (pp \to P\to dd^c) = {\cal L}\, \epsilon\, \sigma_P \, B_d ~,
\end{equation}
where $\cal L$ is the integrated luminosity and $\epsilon$ is an efficiency factor, which we take to be 20\%. The requirement of reproducing the
best-fit signal of 5 events for ${\cal L}=139$~fb$^{-1}$ determines the combination $B_i\, B_d \, \Gamma_P/ M_P$. In table~\ref{tab:Lam} we show the prediction for a scalar resonance in the gluon channel $(i=gg ,\, J_P=0)$ and for a vector resonance in the quark channel $(i=q{\bar q}  ,\, J_P=1)$, where $B_{q{\bar q}}$ is the branching ratio into a single quark channel taking, for simplicity, a universal value of $B_{q{\bar q}}$ valid for all quark species.\footnote{For the general case of non-universal quark branching ratios, the extracted value of $B_{q{\bar q}}$ has to be interpreted as the weighted average $B_{q{\bar q}}= \sum_i C_i B_i / \sum_i C_i$, where the sum extends over the first five quark species.}

\begin{table}[h]
\begin{center}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{|c||c|c||c|c|}
\hline
{\footnotesize Resonance} & {\footnotesize Gluon channel} & {\footnotesize Quark channel} & {\footnotesize Scalar resonance} & {\footnotesize Vector resonance}
 \vspace{-3mm} \\
 {\footnotesize mass} &   {\footnotesize$(i=gg ,\, J_P=0)$} & {\footnotesize$(i=q{\bar q}  ,\, J_P=1)$} & {\footnotesize coupled to gluons} & {\footnotesize coupled to quarks} \\
 \hline
$M_P$ {[TeV]} & $B_{gg}\, B_d \, \Gamma_P/ M_P$ & $B_{q{\bar q}}\, B_d \, \Gamma_P/ M_P$ & $\Lambda_P/ \sqrt{B_d}$  {[TeV]} & $g_{Z'}\, |Q_q|\sqrt{B_d} $\\
\hline
3 & $2.4\times 10^{-4}$ & $1.4\times 10^{-5}$ & 12 & 0.013\\
4  & $3.7\times 10^{-3}$ & $1.3\times 10^{-4}$ &  4.0 & 0.041\\
5  & $4.9\times 10^{-2}$ & $1.2\times 10^{-3}$ &  1.4 & 0.12\\
6  & $6.5\times 10^{-1}$ & $1.2\times 10^{-2}$ &  0.4 & 0.39\\
 \hline
 \end{tabular}
 \renewcommand{\arraystretch}{1.0}
 \caption{The combinations $B_i\, B_d \, \Gamma_P/ M_P$ required to reproduce the ATLAS $dE/dx$ signal, in the case of gluon and quark channels. Also shown are the predictions for the combinations $\Lambda_P / \sqrt{B_d}$ (for the model with scalar resonance coupled to gluons) and  $g_{Z'} |Q_q|\sqrt{B_d} $ (for the model with vector resonance coupled to quarks).}
 \label{tab:Lam}
 \end{center}
\end{table}

The results in table~\ref{tab:Lam} show that, for the gluon channel, the narrow width approximation can hold for $M_P \lesssim 5$~TeV, but it deteriorates at larger masses where the parent's interactions start becoming non-perturbative. The situation is more favourable for the quark channel, where the narrow width approximation can be satisfied in a broader range of $M_P$, as long as both branching ratios $B_{q{\bar q}}$ and $B_d$ are not too small.

A robust and model-independent prediction of our setup with boosted LLPs is an excess of dijet events mediated by resonant parent production. This excess is simply correlated with the anomalous $dE/dx$ events and the ATLAS signal predicts a non-standard contribution to the total dijet cross section at $\sqrt{s} = 13$~TeV
\begin{equation}}% can be used as {equation} or {eqnarray
\sigma_P ({\rm dijet}) = 0.45~{\rm fb} \, \frac{B_{q{\bar q}}}{B_d} ~,
\label{sdijet}
\end{equation}
where we have included a 50\% efficiency factor and summed over five quark species in the final state, with a universal value of $B_{q{\bar q}}$.\footnote{For the general case of non-universal quark branching ratios, $B_{q{\bar q}}$ in eq.~\eqref{sdijet} has to be interpreted as the average $B_{q{\bar q}}= \sum_i  B_i / 5$, where the sum extends over the first five quark species.}

\begin{figure}[t]
\centering
\includegraphics[width=25em]{dijet} \caption{Lower bounds on the relative branching $B_d/B_i$, both for the gluon ($i=qq$) and quark ($i=q\bar q$) production channels, obtained from the CMS limits on dijet searches~\cite{CMS:2019gwf}. Other parameters have been fixed to reproduce the ATLAS $dE/dx$ signal.
\label{fig:dijet}}
\end{figure}

Using the LHC limits on resonant dijet production in gluon and quark channels~\cite{CMS:2019gwf}, we can extract the lower bounds on the ratios $B_d/B_{gg}$ and $B_d/B_{q \bar q}$ shown in Fig.~\ref{fig:dijet}. This figure shows that resonant dijet production at present gives only a mild constraint on the boosted LLPs interpretation of the $dE/dx$ excess.

Since the doubly-charged daughters must carry EW quantum numbers,
another generic phenomenological consequence of our setup is an irreducible Drell-Yan production of daughter pairs. This must be subdominant as compared to production via $P$ decay, because the former would give slower $d$ particles with much larger $dE/dx$.  We estimate that for $m_d \gtrsim 500 \text{ GeV}$ this effect is subdominant.

While excess of dijet events and EW pair-production of doubly-charged particles are unescapable and model-independent consequences of our interpretation of the ATLAS $dE/dx$ signal, other experimental signatures (such as dilepton or missing energy from parent decay, or contact interactions from virtual parent or daughter exchange) could be present in specific model realisations, as will be shown in the following.

We can now gain further insight on the microscopic structure of boosted LLPs by illustrating specific examples of models for the parent and daughter particles.

\subsection*{Scalar resonance coupled to gluons}
A simple microscopic model is given by a heavy singlet scalar $P$ coupled to gluons and doubly-charged daughter particles $d$ as
\begin{equation}}% can be used as {equation} or {eqnarray
\frac{\alpha_s}{\Lambda_P} \, P \, G^2_{\mu\nu} + \kappa\, P\, d^c d~,
\end{equation}
where $\Lambda_P$ is the scale of the dimension-five effective interaction and $\kappa$ is a coupling constant. We do not need to specify the daughter's spin. The model must also include a feeble interaction, possibly described by a higher-dimension effective operator, that allows for $d$ decay, making the daughter metastable.

The parent decay width into gluons is given by
\begin{equation}}% can be used as {equation} or {eqnarray
\frac{\Gamma_{gg}}{M_P}  = \frac{2 \alpha_s^2 M_P^2}{\pi \Lambda_P^2}~.
\end{equation}
The ATLAS $dE/dx$ signal predicts the effective scale of the model, through the combination $\Lambda_P / \sqrt{B_d}$, as shown in table~\ref{tab:Lam} for relevant values of $M_P$.

Since the ratio $M_P/\Lambda_P$ has dimensions of coupling (much like the commonly encountered combination $m_W/v$ in the SM), the model indicates a moderately strongly-coupled UV completion for $M_P \lesssim 5$~TeV, at least for not too small $B_d$. For larger $M_P$, the theory enters a strongly-coupled regime and any perturbative control is lost. Therefore, a generic consequence of this model is the likely existence of new coloured states not far beyond the mass scale of $M_P$, which should not exceed about 5 TeV.

Finally, we remark that our results are unchanged if the parent, instead of being scalar, is a pseudoscalar coupled to $G \widetilde{G}$, since the formul\ae~for the cross section and branching ratio remain the same.

\subsection*{Vector resonance coupled to quarks}
As an alternative microscopic model, one can take the parent resonance to be a $Z'$ boson of a $U(1)'$ gauge group under which at least the daughter particle and first-generation quarks are charged.\footnote{$Z'$ decays to pairs of unit-charge LLPs were discussed in \cite{Bauer:2009cc}.} Taking the simple case of a vector current with coupling constant $g_{Z'}$, the $Z'$ partial width of the decay into each particle pair $\psi$ is\footnote{The normalisation is chosen such that the gauge interaction of the fermionic current is $g_{Z'} Q_\psi Z'_\mu {\bar \psi} \gamma^\mu \psi$.}
\begin{equation}}% can be used as {equation} or {eqnarray
\frac{\Gamma (Z'\to {\bar \psi}\psi)}{ M_{Z'}} = \frac{{\cal N}_\psi \,  Q_\psi^2 \, g_{Z'}^2}{12\pi} ~,
\end{equation}
where $Q_\psi$ is the $\psi$ charge under the new $U(1)'$ gauge group and ${\cal N}_\psi$ is the number of effective species. Quarks correspond to ${\cal N}_q =3$, while the daughter particle  gives
\begin{equation}}% can be used as {equation} or {eqnarray
{\cal N}_d = \left\{
\begin{array}{ll}
\frac{N_d \, \beta (3-\beta^2)}{2} &{\rm for}~J_d=1/2 \\
\frac{N_d \, \beta^3}{2} &{\rm for}~J_d=0
\end{array} \right.
 ~,~~~ \beta = \sqrt{1-\frac{4m_d^2}{M_{Z'}^2}} ~,
\end{equation}
where $N_d$ is the daughter multiplicity.

The ATLAS $dE/dx$ signal gives a prediction for the gauge coupling $g_{Z'}$, up to a coefficient $|Q_q|\sqrt{B_d}$. The prediction is shown in table~\ref{tab:Lam}, under the simplifying assumption of a universal $U(1)'$ charge $Q_q$ for all quarks.\footnote{For non-universal quark charges, one can simply replace $Q_q$ with the weighted average $ \sum_i C_i Q_i / \sum_i C_i$, when it refers to the initial state, and with the average $ \sum_i Q_i / 5$, when it refers to the final state, where the sum extends over the first five quark species.}  As long as $B_d$ is not too small, the new gauge coupling constant $g_{Z'}$ is safely in the perturbative regime in the full range of relevant values of $M_P$.

\begin{figure}[t]
\centering
\includegraphics[width=25em]{dileptons} \caption{Lower bounds on the daughter $U(1)'$ charge in the combinations $\sqrt{{\cal N}_d}\, |Q_d/Q_q|$ (for the dijet channel) and $\sqrt{{\cal N}_d}\, |Q_d/Q_\ell |$ (for the dilepton channel) from CMS searches.
The value of the gauge coupling $g_{Z'}$ has been fixed to reproduce the ATLAS $dE/dx$ signal.
\label{fig:dileptons}}
\end{figure}

As discussed at the beginning of this section, we expect an irreducible contribution to dijet events. Moreover, depending on $U(1)'$ charge assignments, we can also expect new effects in dilepton events, if $Z'$ has a significant decay width into leptons. Interestingly, once the value of $g_{Z'}$ is fixed to reproduce the ATLAS $dE/dx$ excess, the predictions for the dijet and dilepton cross sections are fully determined by $M_P$ and the daughter effective charge $\sqrt{{\cal N}_d}\, Q_d$ in units of the quark charge $Q_q$ and lepton charge $Q_\ell$, respectively. The predictions are independent of $\Gamma_P$ and therefore are not affected by possible $Z'$ decay modes into other particles.

The present LHC limits on resonant contributions to dijet~\cite{CMS:2019gwf} and dilepton~\cite{CMS:2021ctt} events can be translated into lower bounds on $\sqrt{{\cal N}_d}\, |Q_d/Q_q|$  and $\sqrt{{\cal N}_d}\, |Q_d/Q_\ell |$, as shown in Fig.~\ref{fig:dileptons}. The figure shows that dijet limits are easily satisfied as long as the daughter charge is not much smaller than those of quarks. Dilepton searches provide stronger bounds on $Q_d$ with respect to the lepton charge. This may be taken as an indication that vector resonances with suppressed lepton couplings are favoured.

It is interesting to consider a $B$--$L$ gauge boson (such that $Q_q=1/3$ and $Q_\ell = -1$), where the predictions for dijet and dilepton resonant production are correlated. Figure~\ref{fig:dileptons} shows that dileptons are the most efficient channel to test a $B$--$L$ gauge boson. Present searches give the bound $\sqrt{{\cal N}_d}\, |Q_d| \gtrsim 5$, in the most relevant mass window $ M_P \approx 4$--$6~{\rm TeV}$, and therefore require a sufficiently large daughter charge and/or multiplicity. We remark that the lower bound on $\sqrt{{\cal N}_d}\, |Q_d|$ scales as the inverse square root of the experimental efficiency in the $dE/dx$ signal, which we estimated as 20\%. A precise assessment of the bound would require a detailed experimental analysis and is quite sensitive to future statistical improvements.

Depending on $U(1)'$ charge assignment and the spectrum of new particles accompanying the daughters, other phenomenological signatures are possible. An intriguing example is the $Z'$ effect which could explain the recent $W$-mass anomaly claimed by the CDF collaboration~\cite{CDF:2022hxs}. Such an explanation requires a non-vanishing Higgs charge $Q_h$ under $U(1)'$, suppressed lepton charges, and a gauge coupling $g_{Z'}|Q_h| \simeq M_{Z'}/8 \text{ TeV}$~\cite{Strumia:2022qkt}. Therefore, the same $Z'$ can simultaneously fit both the $dE/dx$ excess and the $M_W$ measurement by CDF  if $\sqrt{B_d} Q_q/Q_h = \{ 0.035, 0.082,0.19,0.52\}$ for $M_Z'/{\rm TeV} = \{3,4,5,6 \}$. The emerging picture shows a $Z'$ coupled with comparable strength to quarks, Higgs and daughters, while couplings to leptons must be relatively suppressed. This might be indicative of a vector resonance of a coloured strongly-coupled sector in the multi-TeV range interacting with the Higgs boson, as in composite Higgs models.

\subsection*{EW-charged resonances}

In the cases discussed above, the parent resonance is assumed to be an $SU(2)_L$ singlet. However, a heavy electroweak doublet parent $P$ would in principle work as well. An important qualitative difference arising in this case is that $P$ cannot decay into a pair of the same daughter particle. As a consequence, scenarios with electroweak-charged resonances would be favoured if future data show that tracks with large $dE/dx$ are never accompanied by another ionising track from the recoiling particle.

In the case of an electroweak-doublet parent, a natural choice would be for it to decay into an $SU(2)_L$ triplet with hypercharge one (containing the electric charge-two state) and a doublet, possibly identified with the SM Higgs boson. Notice that the unit charge component of the triplet would give a lower $dE/dx \approx 1 \, \mathrm{MeV g^{-1} cm^2}$ signal, hidden in the large background.

\subsection*{Coloured resonances}

Finally, we note that coloured parent resonances are more difficult to accommodate, if one insists to have decays into a pair of the same kind of daughter.  If the charge-two daughter were coloured it would have a significant QCD production with lower $\beta$, yielding a significant $dE/dx$ signal, close to the upper bound of the dynamic range of the ATLAS detector. Then the only possibility to explain the ATLAS excess is that the decay of the coloured resonance takes place into two different particles, and only the one with $Q \leq 1$ is coloured.

\section{Conclusions}
History has taught us to expect the unexpected in fundamental physics.  Not every discovery is foreseen, nor have they all provided the missing piece in an outstanding theoretical jigsaw puzzle.  This was true  for the archetypal LLP discovery of the muon.  In this discovery the muons were produced from the decays of heavier parent particles, the pions.  It is just a coincidence that the pion and muon masses are so close and, in principle, the parents could have been significantly heavier than the muons, boosting them in the parent rest-frame.

In this work we have considered whether history could repeat itself at the LHC by studying the phenomenology of boosted charged LLPs in $dE/dx$ searches.  We have shown that the phase space they occupy is distinct from commonly-considered scenarios where the LLPs are pair-produced in non-resonant processes. Plausible microscopic models of boosted LLPs exist, are consistent with present experimental limits, and can be searched for at future LHC runs.
A general feature, independent of the specific model realisation of boosted LLPs, is an additional irreducible signature that could be revealed at resonant dijet searches, produced by the coupling of the parent resonance with SM light quarks or gluons. Other, more model-dependent, signatures can be useful to obtain further confirmation of potential discoveries.

An exciting aspect of boosted LLPs is that, so far, they are the only known explanation for the recently reported $dE/dx$ excess by the ATLAS collaboration, consistent with the information from the time-of-flight measurement, suggesting that $\beta \approx 1$. We find that overall the quality of the fit provided by boosted LLPs for the excess is very good and suggests new particles in the TeV range. It is also quite interesting that the excess lies in a low-background region and therefore it can turn into a more-than-$5\sigma$ discovery at the LHC Run 3, if present observations indicate
a real new-physics phenomenon.

Whether this excess will evolve into a full-blown `who ordered that?' discovery will be a question of statistics, a question of systematics, and ultimately a question of corroborating results from CMS.  Nevertheless, even if the excess eventually evaporates, heavy boosted charged LLPs will remain an interesting item on the menu of unexpected discoveries, which should be investigated by LHC experimental collaborations.

\section*{Acknowledgments}
We thank Ismet Siral for clarifications on~\cite{ATLAStalk} as well as the rest of the  ATLAS $dE/dx$ analysis group for comments on the manuscript.

\section*{Note added}
After this paper appeared as a preprint, the ATLAS collaboration presented new results in the search for long-lived multi-charged particles~\cite{ATLAS:2022cob}. Unlike our work, they only consider direct production from Drell-Yan or photon fusion, which gives rise to multi-charged \emph{unboosted} particles that cannot explain the $dE/dx$ excess.
As we discussed in Sec.~\ref{sec:models}, this is however an irreducible complementary signal of our framework, for light-enough daughter particles. The ATLAS results for $Q=2$ exclude $m_D < 1.05\,\mathrm{TeV}$ at $95 \%$ C.L. for fermionic daughters while showing a mild excess, with  $4$ observed events in a region with $1.5$ expected background. For scalar daughters the corresponding limit inferred from the analysis in~\cite{Altakach:2022hgn} is $m_D \lesssim 700 \text{ GeV}$.

Recently, following previous hints~\cite{CMS:2019gwf, Dobrescu:2018psr}, the CMS collaboration has presented results for the dedicated analysis of heavy resonances decaying into two pairs of jets~\cite{CMS:2022usq}, a topology closely related to our explanation of the $dE/dx$ anomaly, excluding parent masses smaller than $7.6\text{ TeV}$ for a benchmark choice of couplings. The collaboration reports two anomalous events with parent mass $8$ TeV, decaying into a pair of $2$ TeV daughters. The former  could be potentially identified with our parent particle $P$, whereas the latter, being short-lived, cannot be our daughters $d$, but possibly EW partners of them.

\section*{Appendix}