# Ultracentral heavy ion collisions, transverse momentum and the equation of
state

Govert Nijs  Theoretical Physics Department, CERN, CH-1211 Genève 23,
Switzerland  Wilke van der Schee  Theoretical Physics Department, CERN,
CH-1211 Genève 23, Switzerland  Institute for Theoretical Physics, Utrecht
University, 3584 CC Utrecht, The Netherlands  NIKHEF, Amsterdam, The
Netherlands

###### Abstract

Ultracentral heavy ion collisions due to their exceptionally large
multiplicity probe an interesting regime of quark-gluon plasma where the size
is (mostly) fixed and fluctuations in the initial condition dominate. Spurred
by the recent measurement of the CMS collaboration we investigate the driving
factors of the increase of transverse momentum, including a complete analysis
of the influence of the QCD equation of state. Particularly interesting is the
influence of the centrality selection as well as the initial energy
deposition.

††preprint: CERN-TH-2023-234

Introduction - The hydrodynamic description of the quark-gluon plasma (QGP)
formed in relativistic heavy ion collisions at Brookhaven National Laboratory
and the Large Hadron Collider at CERN is one of the great successes of our
understanding of many-body physics within Quantum Chromodynamics. This
understanding relies on a relatively advanced interplay between an initial
far-from-equilibrium stage, a viscous hydrodynamic stage, a switch to the
hadronic phase and finally a relatively intricate global analysis of among
else the anisotropies of the spectra of the detected particles Busza  _et al._
(2018); Romatschke and Romatschke (2019); Bernhard  _et al._ (2019).

Recently, a novel way to look at the hydrodynamic QGP phase has been suggested
by looking at the mean transverse momentum of the top 1% highest multiplicity
events as a function of that multiplicity Gardim  _et al._ (2020a, b) (see
also Samanta  _et al._ (2023a)). These head-on collisions have a (relatively)
fixed transverse size and multiplicity fluctuations are mostly driven by
initial entropy and temperature fluctuations. According to hydrodynamics the
multiplicity is directly related to the entropy Muller and Rajagopal (2005)
and the transverse momentum to the temperature Gardim  _et al._ (2020b). As
such their relation has direct information on the equation of state (EoS)
Gardim  _et al._ (2020a). Together with a later more refined hydrodynamic
prediction from _Trajectum_ Nijs and van der Schee (2022); Giacalone  _et al._
(2023) this has recently been measured by the CMS collaboration CMS (2023);
Bernardes (2023) (see also results by ATLAS ATLAS (2023)), the result of which
is shown in Fig. 1.

It is important that both the experimental and theoretical systematic
uncertainties here are extremely small, often below 0.1%. The reason is that
Fig. 1 shows a ratio with respect to ‘normal’ central collisions in the 0–5%
class and that in this ratio experimental or theoretical systematic
uncertainties mostly cancel. The remaining non-trivial slope is conjectured to
equal the speed of sound $c_{s}^{2}$ at the temperature at the start of
hydrodynamics Gardim  _et al._ (2020a).

Given the importance of a possible direct extraction of the EoS this Letter
will give a complete overview regarding this measurement, including its
relation to the speed of sound. We will find that the EoS is only a
contributing factor, with important contributions from both the initial state
and the exact way how the centrality classes are determined.

Figure 1:  We show the mean transverse momentum versus multiplicity for
ultracentral events (up to the 0–0.001% class) normalized by the central 0–5%
class for twenty likely parameter settings (grey). The blue band is averaged
from the grey lines (1$\sigma$) together with a fit to Eqn. (1). The green
line shows our most likely result (MAP). CMS preliminary results are shown in
red CMS (2023); Bernardes (2023). The characteristic upward slope can among
else be understood by the speed of sound at the initialization temperature
(see text) Gardim  _et al._ (2020a, b).

Varying the EoS - In this work we will use the _Trajectum_ framework, which
consists of a generalized TRENTo initial stage Moreland  _et al._ (2015); Nijs
and van der Schee (2023), a viscous hydrodynamic stage and conversion to
hadrons that are then evolved using SMASH Weil  _et al._ (2016); Oliinychenko
_et al._ (2020); Sjostrand  _et al._ (2008). A complete description can be
found in Nijs  _et al._ (2021); Nijs and van der Schee (2023); Giacalone  _et
al._ (2023), whereby we take the most recent parameter settings from Giacalone
_et al._ (2023). In this particular work we use SMASH with only resonance
decays without the collision term. In the appendix we verify that the
collision term does not affect the results presented.

In order to test the conjecture that Gardim  _et al._ (2020a)

$p_{T}\propto N_{\rm ch}^{c_{s}^{2}}$ (1)

the most obvious attempt is to vary the speed of sound and see how Eqn. (1) is
affected. In this work we always label fits to Eqn. (1) by $c_{s}^{2}$ even
though we will find that a one-to-one correspondence with the EoS is not
always guaranteed.

Varying the EoS in heavy ion collisions was first done in Sangaline and Pratt
(2016), where also a non-trivial posterior was obtained using a global
analysis that agreed with results from the lattice. Here, this is done in a
slightly different way. In _Trajectum_ we follow Huovinen and Petreczky
(2010); Bernhard (2018), which constructs a continuous hybrid EoS that
interpolates between the hadron resonance gas (HRG) at low temperatures to the
lattice result by HotQCD at high temperatures Bazavov  _et al._ (2014). This
is done in a somewhat intricate way, starting with the HotQCD 10-parameter fit
of the pressure $P$ as function of the temperature $T$,

$\frac{P}{T^{4}}=\frac{1}{2}\bigl{(}1+\tanh[c_{t}(t-t_{0})]\bigr{)}\biggl{(}\frac{p_{\text{id}}+a_{n}/t+\cdots+d_{n}/t^{4}}{1+a_{d}/t+\cdots+d_{d}/t^{4}}\biggr{)},$
(2)

where $t=T/T_{c}$, $T_{c}=154$ MeV and
$p_{\text{id}}=95\pi^{2}/180=\tfrac{\pi^{2}}{90}N_{\rm dof}$ is the weakly
coupled limit for massless 3-flavor QCD. The ten parameters $t_{0}$, $c_{t}$
and $a_{n}$ till $d_{d}$ are given in Bazavov  _et al._ (2014). Next, using
thermodynamic identities the trace anomaly ($e-3P$) is computed and smoothly
connected with the HRG trace anomaly. Finally, from a starting temperature of
50 MeV the EoS is computed by integrating the trace anomaly Bernhard (2018).
It is important that only for the (physical) parameters in Bazavov  _et al._
(2014) this gives back the same $P/T^{4}$ as in Eqn. (2), where in particular
the parameter $p_{\rm id}$ is a measure for the number of degrees of freedom.

In the current work we vary the EoS in two ways. First, we vary both the
parameters $p_{\rm id}$ and $a_{n}$ such that $P/T^{4}$ in the UV stays fixed.
This is motivated by our prior knowledge that at high temperatures the EoS is
well constrained by our knowledge of perturbative QCD. We note that this
automatically limits large variations in the speed of sound, since the
integral of the speed of sound determines the effective number of degrees of
freedom, and is thereby quite unlike what is done in Sangaline and Pratt
(2016). Secondly, more similar to Sangaline and Pratt (2016), we vary
exclusively $p_{\rm id}$ to see the effect of the number of degrees of freedom
on the speed of sound and the final hadronic observables.

Fig. 2 shows the resulting EoS for several parameter values (top), its
$\langle p_{T}\rangle$ as a function of centrality, and the fitted values
using Eqn. (1) together with the actual values of $c_{s}^{2}$ at temperatures
of $0.3\,$GeV and at the average ultracentral initialization temperature
$\langle\langle T\rangle\rangle=0.42\,$GeV. Here the first average is within
an event, weighted by the energy density, and the second average is done in
the 0–1% centrality bin. It is clear that varying the EoS has a small but
significant effect on the absolute value of $\langle p_{T}\rangle$, but only a
small effect on the slope (bottom-left). Nevertheless, fitting $c_{s}^{2}$ to
these variations we find within uncertainties that there is a decreasing trend
that is consistent with the actual variations of the $c_{s}^{2}$ from the EoS.
Ideally one would like larger variations of $c_{s}^{2}$, but as can be seen it
is challenging to do so while keeping a causal $e/T^{4}$ that approaches the
correct UV limit.

As shown in Giacalone  _et al._ (2023) the absolute change in $\langle
p_{T}\rangle$ and other observables can be compensated by variations in other
parameters such as the viscosities. As such Giacalone  _et al._ (2023) did not
obtain a posterior constraint on $a_{n}$, but some mild correlations exist
between other parameters. It would be interesting to have more precise and
direct knowledge of the EoS and Fig. 2 suggests that the ultracentral $\langle
p_{T}\rangle$ can do this, but unfortunately only when the underlying
precision is high.

Figure 2: (top) We show the energy density and speed of sound as a function of
temperature for the variations as described in the text. (bottom) The varying
EoS has a significant effect on $\langle p_{T}\rangle$, but less so on the
slope as in Fig. 1. (bottom-right) Fitting the $c_{s}^{2}$ using Eqn. (1) we
nevertheless see a decreasing trend that is consistent with the EoS trend
(top) at temperatures 0.3 and $0.42\,$GeV (shown in red). Figure 3: (left) As
in Fig. 1 we show the ultracentral $\langle p_{T}\rangle$ for variations of
the nucleon width $w$ (purple), the initial state fluctuations by
$\sigma_{\text{fluct}}$ (red) and the TRENTo parameter $q$ (green). The lines
extend to the 0.001% centrality bin, so that it is clear that more
fluctuations (red) lead to a higher ultracentral multiplicity and $\langle
p_{T}\rangle$, but do not change the slope. On the other hand and $q$ affects
the slope significantly, which can be explained by the fact that ultracentral
collisions in the case of large $q$ select a smaller initial QGP as measured
by the RMS of the energy density (right). Figure 4:  One of the central
results of this Letter is that the slope of the ultracentral $\langle
p_{T}\rangle$ depends strongly on the centrality definition. Perhaps
unsurprisingly, selecting centrality using a calorimeter (upper-left and
upper-middle) causes a selection effect that favors a larger slope towards
ultracentral collisions. This effect is stronger when defining centrality in
the same region as the $\langle p_{T}\rangle$ due to auto-correlation (upper-
left). Similar effects are present when using charged multiplicity, whereby
especially a $p_{T}$ cut on multiplicity favors higher $\langle p_{T}\rangle$
(right panels).

Initial state and centrality definitions - Apart from the EoS we varied
several other aspects that include the shear viscosity, bulk viscosity,
nucleon size, the hydro starting time, initial energy fluctuations and the
TRENTo initial state parameters $p$ and $q$. The latter parameters are used
for a phenomenological description of the initial transverse energy density
$e(x_{\perp})$ as Moreland  _et al._ (2015); Nijs and van der Schee (2023)

$e(x_{\perp})\propto\left(\frac{\mathcal{T}_{L}(x_{\perp})^{p}+\mathcal{T}_{R}(x_{\perp})^{p}}{2}\right)^{q/p},$
(3)

in terms of the participant thickness functions of the left- and right-moving
nuclei $\mathcal{T}_{L,R}$. The mean of $\mathcal{T}_{L,R}$ is proportional to
the nucleon density of the nuclei and its fluctuations are given by a gamma
distribution with parameter $\sigma_{\text{fluct}}$.

Fig. 3 shows the analogue of Fig. 1 while varying the nucleon width $w$, $q$
and $\sigma_{\text{fluct}}$, which are the parameters that most directly
influence the ultracentral ratios as shown. As anticipated in Gardim  _et al._
(2020a) increasing $\sigma_{\text{fluct}}$ leads to a higher ultracentral
$\langle p_{T}\rangle$ and multiplicity (all curves stop at the class
0–0.001%), but the slope is unaffected. Perhaps surprisingly, the $q$
parameter defined above affects the slope most significantly. The reason is
that a large $q$ in combination with the ultracentral selection favors a small
initial droplet (Fig. 3, right). This can be understood, as a large $q$ can
mimic binary scaling energy deposition, which as opposed to participant energy
deposition is most effective in depositing entropy when all nucleons are in a
small area. A constant area is a crucial assumption in Gardim  _et al._
(2020a), as a smaller area automatically implies a stronger radial flow and
higher $\langle p_{T}\rangle$. We find that a large nucleon width changes the
shape of the curve, but does not change the slope significantly.

Another crucial aspect of the ultracentral $\langle p_{T}\rangle$ is the
(experimental) centrality selection (see also Nijs and van der Schee (2022);
Samanta  _et al._ (2023a, b)). Traditionally centrality has been regarded as a
measure for impact parameter, but the more modern view is event activity.
Event activity can however be estimated in several ways, each looking at
activity in different areas of phase space. The two most standard estimators
are forward multiplicity (for ALICE V0 this is given by V0A ($2.8\leq\eta\leq
5.1$) and V0C ($-3.7\leq\eta\leq-1.7$) or forward energy as measured by a
calorimeter. Alternatively it is possible to use mid-rapidity activity (see
e.g. Aad  _et al._ (2023)), which however has the feature that centrality
often self-correlates with the observable measured. Lastly, while
experimentally challenging, a lack of detected spectators in the ZDC detector
can also signal a more central collision. Within _Trajectum_ we can mimic all
these strategies and hence make an apples-to-apples comparison with the
measurement.111For mid-rapidity versus forward centrality selection it is
important that the hydrodynamics in this work is boost invariant. The
difference between the two hence only includes thermodynamic fluctuations at
freeze-out, but ignores fluctuations in the initial condition. This hence
serves as a lower bound on the difference in the selection.

Fig. 4 shows six different centrality definitions, varying between
calorimetric (transverse energy) or charged multiplicity and varying rapidity
from the forward region (avoiding auto-correlation with the $\langle
p_{T}\rangle$ measurement) to mid-rapidity. For the charged multiplicity it is
important that detectors especially at mid-rapidity often are limited in
$p_{T}$ range, which we hence also varied. The standard _Trajectum_ centrality
selection is done with charged multiplicity at mid-rapidity ($|\eta|<0.8$ and
$p_{T}>0.4\,$GeV), which is used in all other plots.222Experiments like ALICE
almost always use forward centrality to avoid auto-correlation; in _Trajectum_
this is however in general often computationally expensive, as a larger
rapidity range needs to be simulated. Clearly, selecting centrality
calorimetrically favors events with a larger $\langle p_{T}\rangle$ and hence
has a larger slope going towards the most ultracentral events. Importantly,
defining centrality with charged multiplicity of relatively high $p_{T}$
particles (e.g. $p_{T}>0.4\,$GeV) has a similar effect of favoring high
$\langle p_{T}\rangle$ events. It is for this reason that our ‘standard’
centrality selection works relatively well in Fig. 1 even though CMS used a
forward calorimetric centrality definition. Lastly, as expected especially for
the calorimetric definitions the auto-correlated definition (e.g. defining
centrality at mid-rapidity just like $\langle p_{T}\rangle$) has a stronger
effect on $\langle p_{T}\rangle$.

While perhaps less surprising than the centrality cuts, it is clear that the
ultracentral $\langle p_{T}\rangle$ slope also depends on the $p_{T}$ cut of
$\langle p_{T}\rangle$ itself. This is shown in the Appendix, but here we
point out that an even more accurate theory-experiment comparison can be made
by focusing on charged particles with $p_{T}>0.3\,$GeV, since CMS relies on
extrapolation for particles with lower $p_{T}$. The fact that the agreement is
slightly different from Fig. 1 indicates that the CMS extrapolation differs
from the _Trajectum_ spectrum. In the Appendix we also made a comparison with
the ATLAS result using $p_{T}>0.5\,$GeV.

Figure 5: We show _Trajectum_ predictions for xenon (top) and oxygen (bottom)
collisions at LHC energies analogous to Fig. 1.

Discussion - One of the curious features of Fig. 1 is the local minimum that
arises around 1% centrality. As far as we are aware there is no
straightforward physical interpretation since all initial quantities such as
average entropy density, average temperature, impact parameters, size of the
QGP are monotonic around 1% (see Appendix). Nevertheless, around 1% a
qualitative change occurs where the impact parameter stops varying.

At first the results of this work may seem disappointing. After all, the
ultracentral $\langle p_{T}\rangle$ may not be a direct measurement of the
speed of sound. We, however, think this measurement may potentially be even
more exciting. Hydrodynamics and the EoS are a crucial aspect of the increase
in the $\langle p_{T}\rangle$, but several other factors contribute, including
subtle size variations (depending on parameters, such as $q$ discussed before)
as well as experimental considerations that includes the centrality
definition. The fact that the final measurement with its precision due to
experimental and theoretical uncertainty cancellations was predicted shows the
importance and success of understanding all stages of heavy ion collisions at
high precision and in particular highlights a novel possibility to constrain
the initial energy deposition.

In the future it would be especially interesting to study the ultracentral
$\langle p_{T}\rangle$ in other systems, such as xenon (see also ATLAS (2023))
and oxygen. Predictions from _Trajectum_ are provided in Fig. 5, where it can
be seen that the xenon uncertainty is considerably wider (presumably due to
its non-trivial shape Bally  _et al._ (2022a, b)). Oxygen collisions fluctuate
more due to their size, but the slope is relatively robust and consistent with
PbPb. Seeing such a slope in Oxygen collisions at the LHC would be a unique
signature of hydrodynamics in small droplets of QGP.

Acknowledgements - We thank Cesar Bernardes, Somadutta Bhatta, Giuliano
Giacalone, Wei Li, Jean-Yves Ollitrault and Jurgen Schukraft for interesting
discussions.

## References

  * Busza  _et al._ (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Romatschke and Romatschke (2019) P. Romatschke and U. Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Bernhard  _et al._ (2019) J. E. Bernhard, J. S. Moreland, and S. A. Bass, Nature Phys. 15, 1113 (2019).
  * Gardim  _et al._ (2020a) F. G. Gardim, G. Giacalone, and J.-Y. Ollitrault, Phys. Lett. B 809, 135749 (2020a), arXiv:1909.11609 [nucl-th] .
  * Gardim  _et al._ (2020b) F. G. Gardim, G. Giacalone, M. Luzum, and J.-Y. Ollitrault, Nature Phys. 16, 615 (2020b), arXiv:1908.09728 [nucl-th] .
  * Samanta  _et al._ (2023a) R. Samanta, S. Bhatta, J. Jia, M. Luzum, and J.-Y. Ollitrault, (2023a), arXiv:2303.15323 [nucl-th] .
  * Muller and Rajagopal (2005) B. Muller and K. Rajagopal, Eur. Phys. J. C 43, 15 (2005), arXiv:hep-ph/0502174 .
  * Nijs and van der Schee (2022) G. Nijs and W. van der Schee, Phys. Rev. C 106, 044903 (2022), arXiv:2110.13153 [nucl-th] .
  * Giacalone  _et al._ (2023) G. Giacalone, G. Nijs, and W. van der Schee, Phys. Rev. Lett. 131, 202302 (2023), arXiv:2305.00015 [nucl-th] .
  * CMS (2023) CMS (CMS), “Extracting the speed of sound in the strongly interacting matter created in relativistic nuclear collisions,”  (2023), CMS-PAS-HIN-23-003 .
  * Bernardes (2023) C. Bernardes, “Extracting the speed of sound in the strongly interacting matter in pbpb collisions,”  (2023), quark Matter 2023.
  * ATLAS (2023) ATLAS (ATLAS), “Measurement of Transverse Momentum Fluctuations in Xe+Xe and Pb+Pb Collisions with ATLAS,”  (2023).
  * Moreland  _et al._ (2015) J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Weil  _et al._ (2016) J. Weil _et al._ (SMASH), Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .
  * Oliinychenko  _et al._ (2020) D. Oliinychenko, V. Steinberg, J. Weil, M. Kretz, J. Staudenmaier, S. Ryu, A. Schäfer, J. Rothermel, J. Mohs, F. Li, H. E. (Petersen), L. Pang, D. Mitrovic, A. Goldschmidt, L. Geiger, J.-B. Rose, J. Hammelmann, and L. Prinz, “smash-transport/smash: Smash-1.8,”  (2020).
  * Sjostrand  _et al._ (2008) T. Sjostrand, S. Mrenna, and P. Z. Skands, Comput. Phys. Commun. 178, 852 (2008), arXiv:0710.3820 [hep-ph] .
  * Nijs  _et al._ (2021) G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. C 103, 054909 (2021), arXiv:2010.15134 [nucl-th] .
  * Nijs and van der Schee (2023) G. Nijs and W. van der Schee, (2023), arXiv:2304.06191 [nucl-th] .
  * Sangaline and Pratt (2016) E. Sangaline and S. Pratt, Phys. Rev. C 93, 024908 (2016), arXiv:1508.07017 [nucl-th] .
  * Huovinen and Petreczky (2010) P. Huovinen and P. Petreczky, Nucl. Phys. A 837, 26 (2010), arXiv:0912.2541 [hep-ph] .
  * Bernhard (2018) J. E. Bernhard,  _Bayesian parameter estimation for relativistic heavy-ion collisions_ , Ph.D. thesis, Duke U. (2018), arXiv:1804.06469 [nucl-th] .
  * Bazavov  _et al._ (2014) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Samanta  _et al._ (2023b) R. Samanta, J. a. P. Picchetti, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 108, 024908 (2023b), arXiv:2306.09294 [nucl-th] .
  * Aad  _et al._ (2023) G. Aad _et al._ (ATLAS), Phys. Rev. C 107, 054910 (2023), arXiv:2205.00039 [nucl-ex] .
  * Bally  _et al._ (2022a) B. Bally, M. Bender, G. Giacalone, and V. Somà, Phys. Rev. Lett. 128, 082301 (2022a), arXiv:2108.09578 [nucl-th] .
  * Bally  _et al._ (2022b) B. Bally _et al._ , (2022b), arXiv:2209.11042 [nucl-ex] .
  * Garcia-Montero  _et al._ (2022) O. Garcia-Montero, J. Staudenmaier, A. Schäfer, J. M. Torres-Rincon, and H. Elfner, Phys. Rev. C 105, 064906 (2022), arXiv:2107.08812 [hep-ph] .

##  I Appendix

Figure 6: We show for the $\langle p_{T}\rangle$ (left) and charged
multiplicity (middle) as a function of centrality the ratio of the full SMASH
result with respect to the approximation of doing resonance decays only.
Within statistical uncertainties the ratio does not depend strongly on
centrality in the ultracentral regime, which is why in this work we used SMASH
in resonance decay mode only (see right figure for the full ratio).

###  I.1 Full SMASH versus resonance only

In this work we showed results using SMASH with resonance decays only, as
opposed to a full hadronic cascade. The obvious advantage is that in this way
the computation is up to twenty times faster since a single hydro event can
now easily be used to generate 50 hadronic events without increasing the
computational time too much. Especially with central collisions and a large
rapidity range this would be prohibitively expensive using the full SMASH
code.

In this short section about Fig. 6 we verify however that due to the ratios
taken this approximation does not affect any of our results (for this we ran a
single run with full SMASH with about 5M events). The reason is relatively
clear, since in the ratio with the 0–5% class many of the variations cancel
(Fig. 6 right). Nevertheless, even without such cancellation it can be seen
that the charged $\langle p_{T}\rangle$ within uncertainties is almost
unaffected by the hadronic cascade (Fig. 6 left). The multiplicity with the
full cascade is about 1% higher (presumably due to proton-antiproton
annihilation into five pions Garcia-Montero  _et al._ (2022)), but the
centrality dependence of this effect is very mild. We note that for more
complicated observables (such as proton elliptic flow) this can be more
complicated, which is one more reason why the ultracentral $\langle
p_{T}\rangle$ is such an insightful observable.

###  I.2 Geometry

For completeness we show the impact parameter $b$, the size of the initial
droplet RMS and the initial temperature as a function of centrality and
normalized to the 0–5% class in Fig. 7. Consistent with Gardim  _et al._
(2020a) we see that towards ultracentral events the impact parameter becomes
almost zero, the size becomes approximately constant and by construction (we
select high multiplicity or high entropy states) the temperature keeps on
increasing. Nevertheless, we see an indication that the size may be decreasing
towards the 0.001% class, which is related to the $q$ parameter as explained
in the main text. If $q$ is sizeable then it is beneficial for a larger
entropy to have participating nucleons distributed in a smaller space. Keeping
everything fixed a smaller droplet leads to a higher $\langle p_{T}\rangle$.
This can be understood both from a stronger radial flow (due to the larger
gradients), or by noting that at fixed total entropy the smaller droplet has a
larger temperature that in turn leads to a higher $\langle p_{T}\rangle$
Gardim  _et al._ (2020b).

Figure 7: We show the impact parameter $b$, the RMS size of the initial QGP
(middle) and the initial average temperature at the hydrodynamic
initialization (right) as a function of centrality. Note that as in Fig. 1 the
curves are a variation of all parameters, including $\tau_{\rm hyd}$. As such
the temperature will vary significantly among the parameters, but much less so
when taking the ratio with respect to the 0–5% class.

###  I.3 Transverse momentum cut variations

In Fig. 9 we vary the cuts used to obtain both $\langle p_{T}\rangle$ and the
charged multiplicity for our MAP setting. It can be seen that the inclusive
result in the main text (top-left) gives rise to the largest slope. Comparing
the various cuts it can be seen that both tightening the lower or upper cut
reduce the slope. This is unsurprising, since the larger $\langle
p_{T}\rangle$ is a result of shifting particles from lower to higher $p_{T}$
and as such a smaller interval will only capture part of this effect. While
Eqn. (1) is conjectured for the inclusive $\langle p_{T}\rangle$, it is
important that in experiment only particles above a minimum $p_{T}$ cut can be
(accurately) measured. For ALICE this is typically $0.2\,$GeV and for
CMS/ATLAS around $0.3\,$GeV. The lower range hence relies on theoretical
modelling and a more direct comparison would compare for e.g. the
$p_{T}>0.3\,$GeV range. Secondly, mostly due to viscous corrections at freeze-
out and hard QCD processes there is a considerable theoretical uncertainty for
high-$p_{T}$ particles; for that reason it could be beneficial to e.g.
restrict to $p_{T}<2.0\,$GeV.

In Fig. 8 (left) we compare the $p_{T}>0.3\,$GeV setting with the CMS
preliminary result Bernardes (2023). Consistent with the theoretical
predictions CMS also finds a smaller slope than Fig. 1. The CMS uncertainties
are much smaller, since no (unknown) extrapolation to $p_{T}>0$ is necessary,
but it is noteworthy that also the theoretical systematic uncertainties are
significantly smaller.

Lastly, in Fig. 8 (right) we compare to the similar results from the ATLAS
collaboration ATLAS (2023), which use a $0.5\,$GeV $p_{T}$ cut as well as a
multiplicity based centrality selection using the same cut at mid-rapidity. As
in Fig. 1 we find that _Trajectum_ is slightly above the experimental result.

Figure 8: (left) We show the equivalent figure as Fig. 1, but now only taking
charged particles with $p_{T}>0.3\,$GeV. This has the advantage that
experimentally no extrapolation to $p_{T}>0$ is required, and hence a full
theory to experiment comparison can be made. Perhaps surprisingly, also the
theoretical systematic uncertainty is significantly smaller. (right) We show
the equivalent figure as Fig. 1 for our MAP settings, but now only taking the
same cuts and centrality selection as the $N_{\rm ch}^{\rm rec}$ method of
ATLAS (2023). Similar to Fig. 1 we find a good agreement, albeit with a
slightly higher slope than found in experiment. Figure 9:  In the main text we
always show the inclusive $\langle p_{T}\rangle$ (top-left, _Trajectum_ does
not produce hadrons with $p_{T}$ larger than $10\,$GeV). In the bottom row we
vary the lower $p_{T}$ cut, whereas the upper row also varies the upper cut-
off.

###  I.4 Predictors

In Gardim  _et al._ (2020a) relatively simple arguments were used to obtain
Eqn. (1), which includes assuming that multiplicity is proportional to the
total initial entropy and that the $\langle p_{T}\rangle$ is proportional to
the initial average temperature. As such it is interesting to verify whether
these arguments are consistent with the full _Trajectum_ simulations as done
in this work. Some insight can be obtained by looking at Fig. 10, which shows
temperature versus entropy (left), temperature versus multiplicity (middle)
and finally transverse momentum versus multiplicity (right), all for the same
MAP settings.

The second and third slope are consistent, which could indicate that the
temperature may indeed be a good estimator for the transverse momentum, at
least for our ’standard’ centrality selection (see Fig. 4 in the main text).
Switching from entropy to multiplicity decreases the slope and increases the
fluctuations (note the x-axis shift), which at least partly can be attributed
to Poisson fluctuations (reference Gardim  _et al._ (2020a) estimates an
effect of 10% for the ATLAS acceptance, but this will be larger in our
simulation since we use a rapidity range of only $|\eta|<1.5$).

We note that to then relate the slope to the speed of sound one also needs a
constant initial size of the plasma, which is not entirely trivial (see Fig.
7).

Figure 10: We show temperature versus total entropy (left), temperature versus
multiplicity (middle) and $\langle p_{T}\rangle$ versus multiplicity (right).
If multiplicity is proportional to the total entropy and the mean tranverse
momentum to the temperature then all these curves should be equivalent. We see
that the slope in particular decreases when going from entropy to
multiplicity, which at least partly can be explained by Poisson fluctuations.
For larger centralities around 5% we also see a non-trivial slope in the
average temperature that is not present for the $\langle p_{T}\rangle$. In
that regime we hence do not exactly satisfy the linear relation between
$\langle p_{T}\rangle$ and temperature as proposed in Gardim  _et al._
(2020a).


# EnergyAnalyzer: Using Static WCET Analysis Techniques to Estimate the Energy
Consumption of Embedded Applications

Simon Wegener , Kris K. Nikov , Jose Nunez-Yanez
and Kerstin Eder S. Wegener is with AbsInt, Saarbrücken, Germany. (e-mail:
swegener@absint.com)K. Nikov was with the University of Bristol, UK. (e-mail:
krisknikov@gmail.com)J. Nunez-Yanez is with Linköping University, Sweden.
(e-mail: jose.nunez-yanez@liu.se)K . Eder is with the University of Bristol,
UK. (e-mail<EMAIL_ADDRESS>
###### Abstract

This paper presents EnergyAnalyzer, a code-level static analysis tool for
estimating the energy consumption of embedded software based on statically
predictable hardware events. The tool utilises techniques usually used for
worst-case execution time (WCET) analysis together with bespoke energy models
developed for two predictable architectures—the ARM Cortex-M0 and the Gaisler
LEON3—to perform energy usage analysis. EnergyAnalyzer has been applied in
various use cases, such as selecting candidates for an optimised convolutional
neural network, analysing the energy consumption of a camera pill prototype,
and analysing the energy consumption of satellite communications software. The
tool was developed as part of a larger project called TeamPlay, which aimed to
provide a toolchain for developing embedded applications where energy
properties are first-class citizens, allowing the developer to reflect
directly on these properties at the source code level. The analysis
capabilities of EnergyAnalyzer are validated across a large number of
benchmarks for the two target architectures and the results show that the
statically estimated energy consumption has, with a few exceptions, less than
1% difference compared to the underlying empirical energy models which have
been validated on real hardware.

##  1 Introduction

Safety-critical embedded systems are used in various domains such as
transportation, aerospace, medical devices, and industrial control systems.
These systems are designed to meet certain non-functional requirements, such
as timing or energy usage constraints, in addition to functional requirements.
The satisfaction of these non-functional requirements is essential for the
correct operation of the system and the safety of its users. Failure to meet
these requirements can result in catastrophic consequences, such as loss of
life or severe financial losses. Therefore, it is necessary to ensure that
these systems are designed and implemented with reliable guarantees for their
non-functional requirements.

For timing constraints, reliable guarantees can be obtained by using sound
timing analysis methods. Timing analysis is a technique used to analyse the
temporal behaviour of a system and predict the worst-case execution time
(WCET) of tasks and other timing properties. Accurately determining a bound
for the WCET of a task is essential for ensuring that a system meets its
timing constraints and avoiding potential hazards.

Energy consumption is another crucial non-functional requirement for embedded
systems. Energy usage constraints are becoming more and more important due to
the increasing use of battery-powered and energy-constrained devices.
Moreover, reducing energy consumption can increase the lifetime of the device
and reduce its operating costs. However, ensuring that a system meets its
energy usage constraints is a challenging task, as energy consumption is
highly dependent on the system’s workload, input data, and hardware
characteristics. In contrast to timing analysis, which has a well-established
theoretical foundation, creating an energy model that yields safe yet tight
bounds for energy consumption is almost impossible. There are two primary
reasons for this. First, energy consumption is measured in physical units
(Joule), whereas processor cycles are a logical unit of time. Moreover, the
amount of energy consumed by a processor is highly specific to the actual
device. Two processors from the same production batch may already show a small
difference in energy consumption. Additionally, the amount of energy consumed
by one and the same processor may increase over time when the silicon inside
the processor degrades [GVGGS18]. Second, the actual amount of energy consumed
depends on the switching activity in the processor, which is highly data-
dependent. Thus, creating an energy model requires measuring all possible
input combinations for each instruction, which is usually not feasible
[MKE18]. To address these limitations, several research works proposed using
empirical methods to characterise energy models. For example, Georgiou et al.
[GKCE17] suggest using pseudo-randomly created data to characterise an
Instruction Set Architecture (ISA) energy model. This approach reduces the
number of input combinations needed to create the energy model and allows for
faster evaluation of the model.

In recent years, researchers have proposed using event counters to create more
accurate energy models for predictable architectures. Event counters are
hardware components that count the number of times certain events occur during
execution, such as instructions executed or cache misses. By using event
counters to create an energy model, the model can accurately capture the
energy consumption of a more diverse set of programs. Furthermore, event
counters are available on many modern processors, which makes the proposed
energy models more accessible to developers. Pallister et al. [PKME17]
proposed an event counter-based method for data-dependent energy modelling,
which is a more accurate way of modelling energy consumption for systems that
process variable data sets. The proposed method identifies the relationship
between the input data and the processor’s energy consumption and uses this
relationship to create an energy model. The authors evaluated their method on
two different processors and found that it provided more accurate predictions
of energy consumption than previous methods.

EnergyAnalyzer, presented in this paper, is a novel tool for static energy
consumption analysis. It focuses on estimating the energy consumption of
programs running on two specific architectures: the Gaisler LEON3
microprocessor [Leo], a radiation-tolerant microprocessor commonly used in the
space communications sector, and the ARM Cortex-M0 microcontroller [Cor],
known for its ultra-low power capabilities. The key features of EnergyAnalyzer
include:

  1. 1.

Utilising standard techniques from WCET analysis, for static worst-case energy
consumption (WCEC) analysis.

  2. 2.

Incorporating accurate microarchitectural analysis and energy models [NGC+22,
NMW+22] for the two target architectures—ARM Cortex-M0 and Gaisler LEON3.

  3. 3.

When validated against model predictions using real-time samples, the static
analysis shows ¡1% difference in estimated energy for the vast majority of
tested benchmarks.

##  2 Related Work

Several studies have attempted to construct worst-case energy models capable
of capturing the WCEC at the ISA level [JML06, WDH+15]. In Jayaseelan et al.
[JML06], the authors bound the WCEC on a simulated processor by maximizing the
switching activity factor for each simulated component to obtain a WCEC cost
for each ISA instruction. Although this method retrieves the WCEC for all the
ISA instructions, it could result in significant overestimation because the
absolute worst-case on the hardware simulation used for the energy model’s
characterization phase might be infeasible to be triggered by any program on
the actual hardware implementation of the same architecture. Additionally, the
approach is not feasible on physical hardware because there is no practical
way of maximizing the switching activity on hardware. To construct an
equivalent ISA energy model for a fabricated processor, one would need to
exhaustively search all combinations of valid data for the operands of an
instruction, making it infeasible in most cases due to the huge space of
possible input data combinations for each ISA instruction. Although Wägemann
et al. [WDH+15] claimed to have constructed an energy model capable of
capturing the WCEC, they did not provide the specifics of the model’s
characterization. Nevertheless, they tested this energy model on a benchmark
and admitted that such an absolute energy model could lead to significant
overestimations, making the retrieved energy consumption bounds less useful.

Ideally, data-sensitive energy models would be created to capture the energy
cost of executing an instruction based on the circuit switching activity
caused by the operands used. Such models can potentially capture the WCEC of a
program without overestimating the absolute energy consumption models.
However, recent work has demonstrated that finding the data that triggers the
WCEC is an NP-hard problem and that no practical method can approximate tight
energy consumption upper bounds within any level of confidence [MKE18].
Therefore, the authors of Georgiou et al. [GKCE17] suggested using pseudo-
randomly created data to characterize an ISA energy model, as their empirical
evidence showed that such models tend to be close to the actual worst case.
Although this approach is expected to yield loose upper-bound energy
consumption estimations, their experimental results showed a low level of
underestimation of the WCEC (less than 4%) for the programs tested. Such
estimations can still provide valuable guidance to the application programmer
to compare coding styles or algorithms in terms of resource consumption.

This work aims to estimate the WCEC using a similar approach to provide soft
but tight energy bounds. Design decisions can be made based on empirical
investigations to determine the level of over-provisioning that ensures the
required level of dependability for the given application. This approach
avoids the huge over-estimations and associated over-engineering at the system
level that absolute worst-case energy models inherently exhibit.

##  3 Tool Overview

Over the last several years, a more or less standard architecture for static
timing-analysis tools has emerged [WEE+08, HWH95, Erm03, TFW00], which is also
implemented in AbsInt’s WCET analyser aiT. One can distinguish four major
building blocks:

  1. 1.

The decoder translates the executable to an internal form that is used by the
other parts (value analysis, microarchitectural analysis, etc). Architecture
specific patterns decide whether an instruction is a call, branch, return, or
just an ordinary instruction. This knowledge is used to form the basic blocks
of the control-flow graph (CFG). Then, the control-flow between the basic
blocks is reconstructed. In most cases, this is done completely automatically.
However, if a target of a call or branch cannot be statically resolved, then
the user needs to write some annotations to guide the control-flow
reconstruction.

  2. 2.

Afterwards, the value analysis determines safe approximations of the values of
processor registers and memory cells for every program point and execution
context. These approximations are used to determine bounds on the iteration
numbers of loops and information about the addresses of memory accesses. Value
analysis information is also used to identify conditions that are always true
or always false. Such knowledge is used to infer that certain program parts
are never executed and therefore do not contribute to the worst-case resource
consumption. Value analysis is again architecture-dependent.

  3. 3.

The microarchitectural analysis then determines upper bounds for the execution
times of basic blocks by performing an abstract interpretation of the program
execution on the particular architecture, taking into account its pipeline,
caches, memory buses, and attached peripheral devices. The microarchitectural
analysis is even more architecture-dependent than the decoder and value
analysis, as the specification of the ISA alone does not suffice to create an
abstract model of the hardware’s timing behaviour, but the particular
specifics of a particular processor implementing this specification must be
taken into account (e.g., cache size, buffers, pre-fetching, etc). The
microarchitectural analysis is usually a composition of both pipeline and
cache analysis.

  4. 4.

Using the results of the preceding analysis phases, the path analysis phase
searches for the worst-case execution path. The analysis translates the
control-flow graph with the basic block timing bounds determined by the
microarchitectural analysis and the loop bounds derived by the value analysis
into an Integer Linear Program (ILP). The solution of the ILP yields a worst-
case path together with a safe upper bound of the WCET. Path analysis is
generic, i.e., does not depend on the target architecture.

The structure of EnergyAnalyzer is similar to the structure of aiT. In fact,
both tools share most components. In particular, they both use the same
decoder for CFG reconstruction and the same value, loop, control-flow, and
path analyses. Only the microarchitectural analysis differs. aiT uses a
microarchitectural timing model to derive safe upper bounds of the WCET for
each instruction. In contrast, EnergyAnalyzer employs the microarchitectural
energy models presented in Section 4 to estimate the energy consumption of
each instruction. Note, however, that for EnergyAnalyzer, the results of the
path analysis have to be interpreted in a slightly different way, because the
microarchitectural energy model of EnergyAnalyzer only produces an estimate of
the WCEC. Hence, during path analysis, worst-case execution frequencies of
basic blocks are combined with an approximate energy model to produce a tight
estimate of the worst-case energy behaviour which is not necessarily an upper
bound.

##  4 Microarchitectural Energy Analysis

A key component of EnergyAnalyzer is the underlying use of accurate energy
models for the target microarchitectures. Several types of energy modelling
techniques were explored and it was decided that hardware event-based power
modelling is the most suitable approach, since it is a popular and accurate
technique used for both CPU and full system modelling. In their research,
Rodrigues et al. [RAKK13] conducted a systematic review of Performance
Monitoring Unit (PMU) events also referred to as Performance Monitoring
Counters (PMCs) commonly used in modern microprocessors. They showed the
effectiveness of these events in characterising and modelling dynamic power
consumption. Several other studies have also explored accurate power modelling
[NYL13, RPBA+14, WDH+17, SSEM21, NNY20].

The PMC-based energy consumption estimation models were obtained via Ordinary
Least Squares [KNNL05] linear regression analysis, where coefficients,
$\beta_{x}$, are determined for each counter, $C_{x}$, to predict the overall
energy cost, i.e., $E=\sum_{x}(\beta_{x}\times C_{x})+\alpha$, with $\alpha$
being the residual error term. The coefficients $\beta_{x}$ are the constants
in the energy model that are program independent while the counters $C_{x}$
are the variables that depend on the program and its input. For a specific
program with known counter values, the energy model can be used to estimate
the energy consumed during the program’s execution.

For static-analysis-based energy consumption estimation, the overall energy
consumption estimate of a piece of code is typically constructed from the
estimates of the ISA basic blocks of the program. Thus, a PMC-based energy
model can enable energy consumption estimation via static analysis only if the
counters used for the modelling and prediction can be statically predicted at
the ISA basic block level. In order to make the model scalable for block-level
static analysis, we have trained without using an intercept, so the residual
is absorbed into the other event weights. This means that at time 0 the energy
predicted is also 0J. We have also used a Non-Negative Least Squares (NNLS)
solver to guarantee positive weights for all the events in the final energy
model, thus always guaranteeing predictable positive energy consumption values
from the model at discrete time slices [LH95].

The accuracy of the model has been evaluated by using PMC data from a test set
with the generated model equations. The measured power or energy values are
then compared to the estimations obtained from the model. The percentage
difference or mean absolute percentage error (MAPE) between them can be used
as an objective metric to quantify model accuracy. Several different models
for each platform were identified and the best performing ones were integrated
into EnergyAnalyzer. Additional details on model generation and validation
techniques used for both target platforms can be found in the accompanying
papers [NMW+22, NGC+22].

###  4.1 ARM Cortex-M0 Setup and Models

The target platform on which the Cortex-M0 models were developed and validated
is the STM32F0-Discovery board, which features the STM32F051 microprocessor
[STM]. The platform does not feature an on-chip PMU. Thus, a special
methodology was developed to obtain the necessary PMC information, using an
extended version of the Thumbulator instruction set simulator [Thu]. The
target platform allows ten different configurations, depending on CPU
frequency, wait states for flash memory access, and whether instruction pre-
fetch is enabled or not. Further details on how the energy analysis models
have been generated including a breakdown of the available hardware
configurations and associated model weights and performance are presented in
Nikov et al. [NGC+22].

We selected the following energy consumption model of the ARM Cortex-M0, using
six statically predictable PMCs for integration into EnergyAnalyzer. The model
offers an estimation error to physical measurements, calculated as Mean
Absolute Percentage Error (MAPE), of 2.8%, for all the data points used for
training and validation. The resulting energy estimation is measured in nJ:

$\displaystyle E_{\text{Cortex-M0}}=$ $\displaystyle\;0.972565030\times
C_{\text{executed instructions without multiplications}}$ $\displaystyle+$
$\displaystyle\;0.652871770\times C_{\text{RAM data reads}}$ $\displaystyle+$
$\displaystyle\;1.031341343\times C_{\text{RAM writes}}$ $\displaystyle+$
$\displaystyle\;1.037625441\times C_{\text{Flash data reads}}$
$\displaystyle+$ $\displaystyle\;1.354953706\times C_{\text{taken branches}}$
$\displaystyle+$ $\displaystyle\;2.274650563\times C_{\text{multiplication
instructions}}$

###  4.2 Gaisler LEON3 Setup and Models

The LEON3 energy models were trained and validated on the GR712RC evaluation
board [GR7]. Similarly to the STM32F0-Discovery board, this platform also does
not feature a PMU. In order to get the PMC measurements for the models, a new,
dual-platform approach using a Kintex UltraScale FPGA board was developed. The
programmable platform was loaded with a synthesised version of the LEON3
coupled together with the LEON3 Statistics Unit (L3STAT [Cob19]). The results
were synchronised with sensor measurements from the GR712RC platform to obtain
the complete data set for model generation and validation. More details on how
the energy analysis models have been generated are presented in Nikov et al.
[NMW+22].

The models presented in that paper describe fine-grained power models which
are trained and validated on all available samples. The models integrated into
EnergyAnalyzer use the same methodology, but with one key difference: the
samples in the data set are aggregated for each benchmark to create code-
block-sized models, making them more coarse-grained and the NNLS solver is
used to generate positive model weights, as required by EnergyAnalyzer. Since
average power models would not be very helpful for this purpose, total energy
consumption is used instead. Additionally, the set of events that is supported
by EnergyAnalyzer is limited, which has resulted in the creation of bespoke
models for integration. A list of all the supported PMCs can be found in Table
1. A subset of them, labelled ISA+Cache and shown in Table 2, has been
selected because it can be statically predicted with high accuracy.

Separate models are generated using each of the PMC subsets. In addition to
using the bottom-up and top-down search algorithms, detailed in [NMW+22], the
relatively small PMC sets also allow for a full-exhaustive search to be used.
The resulting models are then compared against a model that uses all available
PMCs. Table 3 presents the results of the coarse-grained model generation and
evaluation. As expected, the models computed using the larger All Supported
PMC list perform better than the ones using the ISA+Cache list. However, it is
interesting to note that the three search algorithms exploring the ISA+Cache
PMCs all converge on the same model with the _STORE_ event present in all
models generated, regardless of PMC selection and search method. The reason
why the ISA+Cache models perform worse than the All Supported models is that
the _TIME_ event, which is the single best predictor of power/energy according
to previous work [NNYH15, NNY20], is not included in the list. The ISA+Cache
energy model, with a MAPE of $<$8.3%, was integrated into EnergyAnalyzer to
provide energy estimations for the LEON3 in J.

# | Counter | # | Counter | # | Counter
---|---|---|---|---|---
$C_{0}$ | TIME | $C_{5}$ | WBHOLD | $C_{14}$ | LDST
$C_{1}$ | ICMISS | $C_{7}$ | IINST | $C_{15}$ | LOAD
$C_{2}$ | ICHOLD | $C_{11}$ | BRANCH | $C_{16}$ | STORE
$C_{3}$ | DCMISS | $C_{12}$ | CALL |  |
$C_{4}$ | DCHOLD | $C_{13}$ | TYPE2 |  |

Table 1: All supported PMCs by EnergyAnalyzer.

# | Counter | # | Counter | # | Counter
---|---|---|---|---|---
$C_{1}$ | ICMISS | $C_{11}$ | BRANCH | $C_{14}$ | LDST
$C_{3}$ | DCMISS | $C_{12}$ | CALL | $C_{15}$ | LOAD
$C_{7}$ | IINST | $C_{13}$ | TYPE2 | $C_{16}$ | STORE

Table 2: ISA+Cache subset of events.

Model | Expression | MAPE[%]
---|---|---
Train | Test
| Energy [J]
---
All Supported
All Events
| $\text{E}=0.155261+2.94155\text{e-}08\times{C_{0}}$
---
$+2.5661\text{e-}09\times{C_{2}}+9.93453\text{e-}09\times{C_{5}}$
$+8.97535\text{e-}10\times{C_{12}}+3.21255\text{e-}09\times{C_{13}}$
$+6.14384\text{e-}09\times{C_{15}}+4.54827\text{e-}08\times{C_{16}}$
1.14 | 0.29
| Energy [J]
---
All Supported
Bottom-Up
| $\text{E}=0+3.19557\text{e-}08\times{C_{0}}$
---
$+5.79224\text{e-}08\times{C_{16}}$
1.20 | 1.38
| Energy [J]
---
All Supported
Top-Down
| $\text{E}=0.131077+3.13122\text{e-}08\times{C_{0}}$
---
$+9.17778\text{e-}09\times{C_{5}}+2.99043\text{e-}09\times{C_{15}}$
$+3.92999\text{e-}08\times{C_{16}}$
1.02 | 1.54
| Energy [J]
---
All Supported
Full-Exhaustive
| $\text{E}=0.131087+3.13122\text{e-}08\times{C_{0}}$
---
$+9.17779\text{e-}09\times{C_{5}}+2.99043\text{e-}09\times{C_{14}}$
$+3.63095\text{e-}08\times{C_{16}}$
1.02 | 1.54
| Energy [J]
---
IsaCache
All Events
| $\text{E}=0+1.18567\text{e-}06\times{C_{3}}$
---
$+5.9072\text{e-}07\times{C_{12}}+3.88949\text{e-}08\times{C_{13}}$
$+8.03337\text{e-}08\times{C_{14}}+6.89885\text{e-}08\times{C_{16}}$
8.38 | 24.03
| Energy [J]
---
IsaCache
Bottom-Up
| $\text{E}=0+3.93365\text{e-}08\times{C_{7}}$
---
$+1.87111\text{e-}07\times{C_{16}}$
5.84 | 8.24
| Energy [J]
---
IsaCache
Top-Down
| $\text{E}=0+3.93365\text{e-}08\times{C_{7}}$
---
$+1.87111\text{e-}07\times{C_{16}}$
5.84 | 8.24
| Energy [J]
---
IsaCache
Full-Exhaustive
| $\text{E}=0+3.93365\text{e-}08\times{C_{7}}$
---
$+1.87111\text{e-}07\times{C_{16}}$
5.84 | 8.24

Table 3: Coarse-grained Model Results

##  5 Evaluation

We integrated the energy models from Section 4 into EnergyAnalyzer for ARM
Cortex-M0 and EnergyAnalyzer for LEON3, respectively. We evaluated the
integration with the help of the BEEBS benchmark suite [PHB13]. The goal of
the evaluation is to determine how close the statically estimated energy
consumption for a given workload is to the model estimation. Not all of the
BEEBS benchmarks exercise the worst-case path through the program during
execution. Thus, a comparison of the results of the analysis and the actual
measurements would compare in two orthogonal dimensions. First, it would
compare the tightness of the model with respect to the actual hardware
measurements. Second, it would compare the exercised path with the worst-case
path. In order to fix the comparison to one degree of freedom, we compared the
energy estimates obtained from static analysis and those obtained from the
energy model based on the actual PMC measurements from the platforms. The
tightness of the models has already been demonstrated in Section 4.

In contrast to the safety-critical embedded hard real-time software that is
usually analysed with aiT, the BEEBS benchmarks also contain dynamic memory
management using malloc and free. We did not analyse these benchmarks because
the manual annotation effort to get tight results would be too high. Some of
the benchmarks contain computed calls via function pointers that cannot be
resolved automatically. In this case we manually annotated the call targets.
Moreover, we specified constant data in some cases.

For the LEON3, only a subset of the BEEBS benchmarks has been measured on the
hardware setup, because the execution time of some of the benchmarks is too
low to synchronise the FPGA and the ASIC (see Section 4 and accompanying paper
[NMW+22]).

###  5.1 EnergyAnalyzer for ARM Cortex-M0

For some benchmarks, the static analysis was not able to derive all loop
bounds automatically. In this case, we used Thumbulator to derive flow
constraints for the ILP-based path analysis. However, the benchmark might not
exercise the worst-case path, and thus, using the simulation trace might not
result in the worst-case amount of loop iterations for each loop in the
program. For one of the benchmarks— _wikisort_ —the simulation with
Thumbulator fails because the binary allocated only 4096 bytes of stack, but
the routine WikiSort already needs 4520 bytes of stack. This causes a stack
overflow. Hence, some function pointer variables are overwritten, and the
benchmark cannot be executed correctly. We thus excluded the benchmark from
the evaluation. Two of the benchmarks— _qsort_ and _select_ —contain out-of-
bounds accesses.

Table 4 shows the results of the evaluation of EnergyAnalyzer for ARM
Cortex-M0. For 43 benchmarks, the difference between the model and the static
analysis is less than one percent, i.e., the execution path exercised during
the simulator run is the worst-case path. Note that the analysis results
include the energy consumption of the execution of the main routine, which is
not included in the simulator result, which only contains the path between the
start trigger and the stop trigger. However, the contribution of this overhead
is less than one mJ and hence, negligible. For the other benchmarks, the
static analysis selected different paths as worst-case execution paths. The
maximal observed difference between the simulator run and the static analysis
is 109% for benchmark _nsichneu_ , which models a state machine with many
different execution paths, and the static analysis was not able to prune
infeasible paths. Since the path analysis is a worst-case analysis, it
maximises over the possible execution paths. Hence, the path analysis selects
the worst-case combination which differs significantly from the simulated
execution path.

EnergyAnalyzer allows to trade performance for precision by specifying how
many calling and loop contexts should be distinguished during analysis. We
used this feature to increase the analysis precision. The analysis of most
benchmarks takes less than four minutes to complete, with the exception of
five benchmarks (_rijndael_ , _cubic_ , _sqrt_ , _nbody_ , _picojpeg_), which
took between 4 and 59 minutes.

Benchmark | Analysis Result | Model Result | $\Delta$ | Note
---|---|---|---|---
aha-compress |  078.885mJ |  078.828mJ | ¡ 1 % |
aha-mont64 |  099.396mJ |  099.396mJ | ¡ 1 % |
bubblesort | 366.763mJ | 366.762mJ | ¡ 1 % |
cnt |  042.813mJ |  042.804mJ | ¡ 1 % |
compress |  027.895mJ |  027.895mJ | ¡ 1 % |
crc |  009.623mJ |  009.623mJ | ¡ 1 % |
cubic |  007.801J |  004.138J | 89 % | flow constraints
duff |  004.349mJ |  004.349mJ | ¡ 1 % |
edn | 302.762mJ | 302.762mJ | ¡ 1 % |
expint |  043.315mJ |  043.315mJ | ¡ 1 % |
fac |  002.934mJ |  002.904mJ | 1 % |
fasta |  029.383J |  021.100J | 39 % | flow constraints
fdct |  012.292mJ |  012.292mJ | ¡ 1 % |
fibcall |  001.493mJ |  001.493mJ | ¡ 1 % |
fir |  001.994J |  001.994J | ¡ 1 % |
frac |  001.183J |  001.183J | ¡ 1 % |
insertsort |  003.089mJ |  003.089mJ | ¡ 1 % |
janne_complex |  001.402mJ |  001.402mJ | ¡ 1 % |
jfdctint |  031.481mJ |  031.476mJ | ¡ 1 % |
lcdnum | 886.941uJ | 805.000mJ | 10 % |
levenshtein | 400.926mJ | 400.926mJ | ¡ 1 % |
ludcmp | 174.559mJ | 174.559mJ | ¡ 1 % |
matmult-float |  001.537J |  001.537J | ¡ 1 % |
matmult-int | 842.724mJ | 842.649mJ | ¡ 1 % |
minver | 131.316mJ |  084.348mJ | 56 % | flow constraints
nbody |  025.844J |  025.844J | ¡ 1 % |
ndes | 293.387mJ | 293.297mJ | ¡ 1 % |
nettle-arcfour | 105.880mJ | 105.880mJ | ¡ 1 % |
nettle-cast128 |  023.214mJ |  023.211mJ | ¡ 1 % |
nettle-des |  022.595mJ |  022.595mJ | ¡ 1 % |
nettle-md5 |  005.467mJ |  005.467mJ | ¡ 1 % |
nettle-sha256 |  050.507mJ |  050.507mJ | ¡ 1 % |
newlib-exp |  070.439mJ |  070.439mJ | ¡ 1 % |
newlib-log |  052.954mJ |  052.954mJ | ¡ 1 % |
newlib-sqrt |  010.289mJ |  010.289mJ | ¡ 1 % |
nsichneu |  061.017mJ |  029.185mJ | 109 % |
picojpeg |  004.885J |  004.885J | ¡ 1 % |
prime | 209.663mJ | 209.663mJ | ¡ 1 % |
qsort |  027.294mJ |  020.408mJ | 34 % | flow constraints
qurt | 139.891mJ | 139.890mJ | ¡ 1 % |
rijndael |  007.176J |  007.042J | 2 % |
sglib-arraybinsearch |  076.596mJ |  076.596mJ | ¡ 1 % |
sglib-arrayheapsort |  086.857mJ |  086.857mJ | ¡ 1 % |
sglib-arrayquicksort |  065.600mJ |  065.600mJ | ¡ 1 % |
sglib-queue | 126.250mJ | 126.250mJ | ¡ 1 % |
slre | 206.734mJ | 206.734mJ | ¡ 1 % |
sqrt |  011.529J |  011.529J | ¡ 1 % |
st |  004.142J |  002.945J | 41 % | flow constraints
statemate |  013.331mJ |  009.308mJ | 43 % |
stb_perlin |  005.145J |  005.145J | ¡ 1 % |
stringsearch1 |  046.362mJ |  046.362mJ | ¡ 1 % |
strstr |  005.480mJ |  005.480mJ | ¡ 1 % |
trio-snprintf | 105.378mJ |  065.427mJ | 61 % | flow constraints
trio-sscanf | 139.345mJ |  071.618mJ | 95 % | flow constraints
ud |  021.863mJ |  021.862mJ | ¡ 1 % |
whetstone |  022.533J |  016.687J | 35 % | flow constraints

Table 4: Evaluation of the integration of the energy model for the ARM
Cortex-M0 into static energy consumption analysis.

###  5.2 EnergyAnalyzer for LEON3

Some of the BEEBS benchmarks contain floating-point computations. However,
since the FPGA implementation of the LEON3 was built without a FPU, the
benchmarks cannot use floating-point instructions but must use a software
library that emulates these floating-point computations. One of the
benchmarks— _minver_ —computes a matrix multiplication using floating-point
numbers, where one of the matrices is never initialised. We performed both the
standard worst-case analysis for this benchmark and an analysis where we
assumed that the computations only process normalised IEEE 754 floating-point
numbers and zero. This reflects the “flush to zero” option present in many
architectures. The computed energy consumption estimate is then cut in half,
which shows that enabling “flush to zero” in software floating-point
computations can save a lot of energy.

The LEON3 implements the SPARCv8 ISA, which uses register windows for fast
context switches, and for providing hardware support for the call stack.
However, the number of register windows is limited. The particular LEON3 model
used for our experiments, available on the the GR712RC board and its FPGA
equivalent, has eight register windows. Due to the overlapping nature of the
register windows, and their use as a ring buffer, only seven are usable.
Hence, in case a program needs more than seven register windows, the processor
triggers software traps to handle the register window overflow (and
underflow). This happens for two of the benchmarks— _picojpeg_ and _slre_.
Hence, the processor needs to execute trap functions when it detects a
register window overflow or a register window underflow. This causes
additional energy consumption which must be taken into account during a
system-level energy analysis.

Table 5 shows the results of the evaluation of EnergyAnalyzer for LEON3. The
analysis of most benchmarks takes less than four minutes to complete, but for
three benchmarks— _matmult-float_ , _nbody_ , and _picojpeg_ —the analysis
duration was 50 minutes, 45 minutes, and 24 minutes, respectively.

Benchmark | Analysis Result | Model Result | $\Delta$ | Note
---|---|---|---|---
aha-compress |  011.004 J |  011.004 J | 0 % |
aha-mont64 |  007.499 J |  007.491 J | ¡ 1 % |
bubblesort |  003.898 J |  003.889 J | ¡ 1 % |
edn |  039.186 J |  039.186 J | 0 % |
fir | 159.469 J | 159.469 J | 0 % |
frac |  059.391 J |  059.339 J | ¡ 1 % |
levenshtein |  025.506 J |  025.491 J | ¡ 1 % |
ludcmp |  010.992 J |  010.814 J | 2 % |
matmult-float |  002.847 J |  002.822 J | 1 % |
minver |  014.372 J |  004.643 J | 210 % | worst-case
minver |  007.398 J |  004.643 J | 59 % | assumptions
nbody |  004.512 J |  004.496 J | ¡ 1 % |
ndes |  024.828 J |  024.467 J | 1 % |
nettle-aes |  019.401 J |  019.389 J | ¡ 1 % |
nettle-arcfour |  009.644 J |  009.639 J | ¡ 1 % |
nettle-sha256 |  002.763 J |  002.754 J | ¡ 1 % |
newlib-exp |  004.374 J |  004.319 J | 1 % |
newlib-log |  003.284 J |  003.252 J | 1 % |
picojpeg | 503.732 J | 503.918 J | ¡ -1 % | traps
prime |  003.670 J |  003.667 J | ¡ 1 % |
qurt |  008.001 J |  007.958 J | 1 % |
sglib-arraybinsearch |  006.283 J |  006.281 J | ¡ 1 % |
sglib-arrayheapsort |  013.066 J |  013.062 J | ¡ 1 % |
sglib-arrayquicksort |  013.066 J |  013.052 J | ¡ 1 % |
sglib-queue |  013.901 J |  013.900 J | ¡ 1 % |
slre |  014.988 J |  015.261 J | -2 % | traps

Table 5: Evaluation of the integration of the _ISA+Cache_ energy model for the
LEON3 into static energy consumption analysis.

##  6 Integration into TeamPlay Toolchain and Case Studies

EnergyAnalyzer for ARM Cortex-M0 and EnergyAnalyzer for LEON3 can be used as
standalone tools to estimate the energy consumption of embedded software. They
provide a rich and user-friendly graphical user interface to ease the analysis
process. However, they have been developed during the TeamPlay project as part
of a larger toolchain where they enable multi-criteria optimisation in a
compiler, contract-based programming, and energy-aware scheduling. In the
following, we present the integration of EnergyAnalyzer into the WCET-aware C
compiler WCC [FL10].

The mechanisms implemented to integrate EnergyAnalyzer within WCC mirror the
mechanisms in place to perform WCET analysis using AbsInt’s WCET analyser aiT.
XTC files [Abs] are used to call aiT and EnergyAnalyzer in batch mode (i.e.,
without graphical user interface). An XTC file specifies the binary to be
analyzed, the entry point, the path to an annotation file which contains
details about the target architecture configuration as well as user-provided
annotations like flow facts, and the path to an XML report file. This XML
output file is then parsed by WCC after invocation of EnergyAnalyzer to
extract the analysis results and import them into WCC’s Low-Level IR at
function and basic block level. This attached energy data can further be
exploited by WCC to perform various compiler-level energy-aware optimisations,
and thus, establishing a smooth flow between compiler-level energy analysis
and optimisation. WCC supports source-level flow facts utilizing ANSI C
pragmas. WCC uses these pragmas to direct integrated analysers during
analysis. In general, a user can annotate their code with loop bounds,
recursion depths, and execution frequency of an instruction relative to some
other instruction. These source-level pragmas are translated within WCC into
AIS2 annotations for aiT and EnergyAnalyzer.

EnergyAnalyzer has been applied to several use cases in the course of the
TeamPlay project. First, it has been used to select candidates from a set of
implementations of compute elements for an optimised convolutional neural
network (CNN). Acting as an evaluation guide, it helped decide which
optimisations should be considered for the final CNN implementation and which
showed unacceptable energy consumption and should not be used. Second, it has
been used to analyse the energy consumption of a camera pill prototype. The
addition of an encryption algorithm showed a significant impact on the energy
consumption of the camera pill that also varied depending on the type of
encryption algorithm, with SPECK being an order of magnitude more energy
efficient than AES and PRESENT. EnergyAnalyzer closely predicted the actual
energy usage that was physically measured on the system, and thus is a viable
method for estimating energy consumption of a system such as the camera pill.
Third, EnergyAnalyzer has been used to analyse the energy consumption of a
piece of satellite software. One of the main challenges of the space industry
is power consumption, as spacecrafts usually have limited access to power
sources. EnergyAnalyzer proved to be a useful tool thanks to its support of
the LEON3 on the GR712RC platform, which is the most common processor ASIC
used by European space companies. The results showed a precise prediction of
the execution energy for the different binaries. Another interesting feature
of EnergyAnalyzer was the result visualisation using call and control-flow
graphs, showing the energy consumption for each of the functions inside the
binary which could be used not only for predicting and minimising energy
consumption but also for qualifying code for space.

More details on the integration of EnergyAnalyzer in the TeamPlay toolchain
and the use cases on which the toolchain has been applied are presented in
Rouxel et al. [RBE+22].

##  7 Conclusion

This paper presents our work on static energy consumption analysis for
embedded systems. We created energy models for two predictable architectures:
the ARM Cortex-M0 and the Gaisler LEON3, both achieving estimation errors of
less than 10% when validated against our two target hardware platforms. Both
models are based on hardware event counters, which can be predicted by static
analysis. The models are integrated into EnergyAnalyzer, a novel tool for
code-level static energy consumption analysis. Our evaluation results show a
good accuracy of the energy consumption analysis. The tool provides a user-
friendly graphical interface to analyse energy consumption at different levels
of granularity, from the entire program to individual functions and basic
blocks. EnergyAnalyzer is part of the TeamPlay toolchain [RBE+22], where it
enables multi-criteria code optimisation during compilation.

We demonstrated the usefulness of our tools in several case studies including
a camera-pill designed for medical diagnosis and a space communications
platform. In each case, the tool provided precise predictions of the energy
consumption and helped to identify energy bottlenecks. In conclusion, our work
provides a useful approach to analyse energy consumption in embedded systems.
A potential avenue for extension is to include peripheral energy consumption
for system-level analysis [WDD+18]. We believe that our approach can help to
design more energy-efficient embedded systems and applications, which is a
crucial step towards sustainable development.

## Acknowledgement

This research was supported by the European Union’s Horizon 2020 Research and
Innovation Programme under grant agreement No. 779882, TeamPlay (Time, Energy
and security Analysis for Multi/Many-core heterogeneous PLAtforms).

## References

  * [Abs] XTC Language Specification, Version 2.7.  https://www.absint.com/xtc/xtc-specification.pdf.  Accessed: 2020-01-13.
  * [Cob19] Cobham Gaisler AB.  GRLIB IP Core User’s Manual, 2019.  https://www.gaisler.com/products/grlib/grip.pdf.
  * [Cor] The ARM Cortex-M0 processor.  https://developer.arm.com/products/processors/cortex-m/cortex-m0.  Accessed: 2018-08-14.
  * [Erm03] Andreas Ermedahl.  A Modular Tool Architecture for Worst-Case Execution Time Analysis.  PhD thesis, Uppsala University, Sweden, 2003.
  * [FL10] Heiko Falk and Paul Lokuciejewski.  A compiler framework for the reduction of worst-case execution times.  Real-Time Systems, 46(2):251–300, 2010.
  * [GKCE17] Kyriakos Georgiou, Steve Kerrison, Zbigniew Chamski, and Kerstin Eder.  Energy transparency for deeply embedded programs.  ACM Trans. Archit. Code Optim., 14(1), Mar 2017.
  * [GR7] GR712RC Dual-Core LEON3-FT Development board.  https://www.gaisler.com/index.php/products/boards/gr712rc-board.  Accessed: 2018-09-11.
  * [GVGGS18] Xinfei Guo, Vaibhav Verma, Patricia Gonzalez-Guerrero, and Mircea R Stan.  When “things” get older: exploring circuit aging in iot applications.  In 2018 19th International Symposium on Quality Electronic Design (ISQED), pages 296–301. IEEE, 2018.
  * [HWH95] Christopher A. Healy, David B. Whalley, and Marion G. Harmon.  Integrating the timing analysis of pipelining and instruction caching.  In 16th IEEE Real-Time Systems Symposium, Palazzo dei Congressi, Via Matteotti, 1, Pisa, Italy, December 4-7, 1995, Proceedings, pages 288–297. IEEE Computer Society, 1995.
  * [JML06] R. Jayaseelan, T. Mitra, and Xianfeng Li.  Estimating the Worst-Case Energy Consumption of Embedded Software.  In Real-Time and Embedded Technology and Applications Symposium, 2006\. Proceedings of the 12th IEEE, pages 81–90, April 2006.
  * [KNNL05] Michael H Kutner, Chris Nachtsheim, John Neter, and William Li.  Applied linear statistical models.  McGraw-Hill Irwin, 2005.
  * [Leo] LEON3FT Fault-tolerant processor.  https://www.gaisler.com/index.php/products/processors/leon3ft.  Accessed: 2018-09-11.
  * [LH95] Charles L Lawson and Richard J Hanson.  Solving least squares problems.  SIAM, 1995.
  * [MKE18] Jeremy Morse, Steve Kerrison, and Kerstin Eder.  On the limitations of analyzing worst-case dynamic energy of processing.  ACM Trans. Embed. Comput. Syst., 17(3):59:1–59:22, February 2018\.
  * [NGC+22] Kris Nikov, Kyriakos Georgiou, Zbigniew Chamski, Kerstin Eder, and Jose Nunez-Yanez.  Accurate energy modelling on the cortex-m0 processor for profiling and static analysis.  In 2022 29th IEEE International Conference on Electronics, Circuits and Systems (ICECS), pages 1–4. IEEE, 2022.
  * [NMW+22] Kris Nikov, Marcos Martinez, Simon Wegener, Jose Nunez-Yanez, Zbigniew Chamski, Kyriakos Georgiou, and Kerstin Eder.  Robust and accurate fine-grain power models for embedded systems with no on-chip pmu.  IEEE Embedded Systems Letters, 2022.
  * [NNY20] Krastin Nikov and Jose Nunez-Yanez.  Intra and inter-core power modelling for single-isa heterogeneous processors.  International Journal of Embedded Systems, 12(3):324–340, 2020\.
  * [NNYH15] Krastin Nikov, Jose L Nunez-Yanez, and Matthew Horsnell.  Evaluation of hybrid run-time power models for the arm big. little architecture.  In 2015 IEEE 13th International Conference on Embedded and Ubiquitous Computing, pages 205–210. IEEE, 2015.
  * [NYL13] Jose Nunez-Yanez and Geza Lore.  Enabling accurate modeling of power and energy consumption in an arm-based system-on-chip.  Microprocessors and Microsystems, 37(3):319–332, 2013.
  * [PHB13] James Pallister, Simon J. Hollis, and Jeremy Bennett.  BEEBS: open benchmarks for energy measurements on embedded platforms.  CoRR, abs/1308.5174, 2013.
  * [PKME17] James Pallister, Steve Kerrison, Jeremy Morse, and Kerstin Eder.  Data dependent energy modeling for worst case energy consumption analysis.  In Proceedings of the 20th International Workshop on Software and Compilers for Embedded Systems, SCOPES 2017, pages 51–59, New York, NY, USA, 2017. Association for Computing Machinery.
  * [RAKK13] Rance Rodrigues, Arunachalam Annamalai, Israel Koren, and Sandip Kundu.  A study on the use of performance counters to estimate power in microprocessors.  IEEE Transactions on Circuits and Systems II: Express Briefs, 60(12):882–886, 2013.
  * [RBE+22] Benjamin Rouxel, Christopher Mark Brown, Emad Ebeid, Kerstin Eder, Heiko Falk, Clemens Grelck, Jesper Holst, Shashank Jadhav, Yoann Marquer, Marcos Martinez De Alejandro, et al.  The teamplay project: analysing and optimising time, energy, and security for cyber-physical systems.  In Proceedings the 2023 Design, Automation & Test in Europe Conference & Exhibition (DATE 2023). EDAA/IEEE, 2022.
  * [RPBA+14] Santhosh Kumar Rethinagiri, Oscar Palomar, Rabie Ben Atitallah, Smail Niar, Osman Unsal, and Adrian Cristal Kestelman.  System-level power estimation tool for embedded processor based platforms.  In Proceedings of the 6th Workshop on Rapid Simulation and Performance Evaluation: Methods and Tools, pages 1–8, 2014.
  * [SSEM21] Adam Seewald, Ulrik Pagh Schultz, Emad Ebeid, and Henrik Skov Midtiby.  Coarse-grained computation-oriented energy modeling for heterogeneous parallel embedded systems.  International Journal of Parallel Programming, 49(2):136–157, 2021\.
  * [STM] STM32F030x4/x6/x8/xC and STM32F070x6/xB advanced ARM-based 32-bit MCUs - Reference Manual.  Accessed: 2020-12-28.
  * [TFW00] Henrik Theiling, Christian Ferdinand, and Reinhard Wilhelm.  Fast and precise WCET prediction by separated cache and path analyses.  Real-Time Systems, 18(2/3):157–179, 2000.
  * [Thu] The Thumbulator Git repository modified for collecting performance monitoring counters for the STM32F0-Discovery board.  branch: prefetch-model, Git-tag: teamplay-D4.5. Accessed: 2020-11-09.
  * [WDD+18] Peter Wägemann, Christian Dietrich, Tobias Distler, Peter Ulbrich, and Wolfgang Schröder-Preikschat.  Whole-system worst-case energy-consumption analysis for energy-constrained real-time systems.  Leibniz International Proceedings in Informatics, LIPIcs 106 (2018), 106:24, 2018.
  * [WDH+15] P. Wägemann, T. Distler, T. Hönig, H. Janker, R. Kapitza, and W. Schröder-Preikschat.  Worst-Case Energy Consumption Analysis for Energy-Constrained Embedded Systems.  In 2015 27th Euromicro Conference on Real-Time Systems (ECRTS), pages 105–114, July 2015.
  * [WDH+17] Matthew J Walker, Stephan Diestelhorst, Andreas Hansson, Anup K Das, Sheng Yang, Bashir M Al-hashimi, and Geoff V Merrett.  Accurate and Stable Run-Time Power Modeling for Mobile and Embedded CPUs.  Ieee Transactions on Computer Aided Design of Integrated Circuits and Systems, 36(1):1–14, 2017.
  * [WEE+08] Reinhard Wilhelm, Jakob Engblom, Andreas Ermedahl, Niklas Holsti, Stephan Thesing, David Whalley, Guillem Bernat, Christian Ferdinand, Reinhold Heckmann, Tulika Mitra, Frank Mueller, Isabelle Puaut, Peter Puschner, Jan Staschulat, and Per Stenström.  The worst-case execution-time problem — overview of methods and survey of tools.  ACM Transactions on Embedded Computing Systems, 7(3):36:1–36:53, May 2008.