11institutetext: INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5,
50125 Firenze, Italy
11email<EMAIL_ADDRESS>22institutetext: Laboratoire de Physique de
l’Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université,
Université de Paris, F-75005 Paris, France

# A new analytical model of the cosmic-ray energy flux for Galactic diffuse
radio emission

Andrea Bracco 1122 Marco Padovani 11 Daniele Galli 11

(Received 15 February 2024; accepted 29 February 2024)

Low-frequency radio observations of diffuse synchrotron radiation offer a
unique vantage point for investigating the intricate relationship between gas
and magnetic fields in the formation of structures within the Galaxy, spanning
from the diffuse interstellar medium (ISM) to star-forming regions.

Achieving this pivotal objective hinges on a comprehensive understanding of
cosmic-ray properties, which dictate the effective energy distribution of
relativistic electrons, primarily responsible for the observable synchrotron
radiation. Notably, cosmic-ray electrons (CRe) with energies between 100 MeV
and 10 GeV play a crucial role in determining the majority of the sky
brightness below the GHz range. However, their energy flux ($j_{e}$) remains
elusive due to solar modulation.

We propose deriving observational constraints on this energy gap of
interstellar CRe through the brightness temperature spectral index of low-
frequency radio emission, here denoted as $\beta_{\rm obs}$. We introduce a
new parametric analytical model that fits available data of $j_{e}$ in
accordance with the $\beta_{\rm obs}$ values measured in the literature
between 50 MHz to 1 GHz for diffuse emission in the Milky Way. Our model
allows to account for multiple observations considering magnetic-field
strengths consistent with existing measurements below 10 $\mu$G. We present a
first all-sky map of the average component of the magnetic field perpendicular
to the line of sight and validate our methodology against state-of-the art
numerical simulations of the diffuse ISM.

This research makes headway in modeling Galactic diffuse emission with a
practical parametric form. It provides essential insights in preparation for
the imminent arrival of the Square Kilometre Array.

###### Key Words.:

magnetic fields – radiation mechanisms: non-thermal – ISM: cosmic rays – radio
continuum: ISM

##  1 Introduction

Low-frequency radio emission observations are paving the way for a
comprehensive study of Galactic magnetic fields traced by synchrotron
radiation. The breakthrough in sensitivity of current and upcoming radio
telescopes, including the Low Frequency Array (LOFAR, van Haarlem et al.,
2013), the Long Wavelength Array (LWA, Dowell et al., 2017), the New Extension
in Nançay Upgrading LOFAR (NenuFAR, Zarka et al., 2012), the C-Band All-Sky
Survey (Jones et al., 2018), and the Square Kilometre Array (SKA, Dewdney et
al., 2009) along with its precursors, promises the most detailed multi-scale
description of the Galactic magnetic field. This description encompasses both
its topology and its strength, ranging from the diffuse interstellar medium
(ISM) to star-forming regions (Heald et al., 2020).

Achieving this paramount objective, however, relies on the thorough
characterization of cosmic-ray (CR) properties, acceleration, and propagation.
These factors determine the effective energy distribution of relativistic
electrons, which are primarily responsible for the observable synchrotron
radiation111Positrons and secondary electrons also contribute to synchrotron
radiation but at a few % level (Orlando, 2018; Ponnada et al., 2024) (e.g.,
Ginzburg & Syrovatskii, 1964; Padovani & Galli, 2018). In particular, low-
energy cosmic-ray electrons (CRe) between 100 MeV and 10 GeV are those
relevant for most of the radio emission detected below the GHz range (Padovani
et al., 2021, hereafter P21). Unfortunately, because of solar modulation, the
energy flux of these electrons ($j_{e}$) cannot be measured from near-earth
direct observations (e.g., Gabici, 2022). Hence, $j_{e}$ is usually
interpolated in the GeV window (e.g., Orlando, 2018; Padovani & Galli, 2018;
Padovani et al., 2018; Unger & Farrar, 2023; Bracco et al., 2023) between the
MeV range observed with the Voyager spacecrafts (Cummings et al., 2016; Stone
et al., 2019) and the hundred-of-GeV range measured by facilities including
Fermi LAT (Ackermann et al., 2010), Pamela (Adriani et al., 2011), and the
Alpha Magnetic Spectrometer (AMS, Aguilar et al., 2014).

Building upon the pioneering works of Rockstroh & Webber (1978), Strong &
Wolfendale (1978), and Strong et al. (2000), and following up the
investigation conducted by P21, in this paper we propose that observational
constraints on the missing energy window of interstellar CRe can be obtained
through the spectral index of low-frequency radio emission, denoted as
$\beta_{\rm obs}$. We introduce a new analytical-parametric model that
accurately fits $j_{e}$, considering the values of $\beta_{\rm obs}$ measured
in the literature between 45 MHz and 408 MHz for Galactic diffuse emission. We
use our models of $j_{e}$ to estimate the strength of the Galactic magnetic
field averaged along the line of sight. Our observational results are
discussed through a comparison with synthetic data generated from state-of-
the-art magneto-hydrodynamic (MHD) simulations of the diffuse ISM.

By incorporating the constraint of $\beta_{\rm obs}$ in the shape of $j_{e}$,
we demonstrate an enhanced ability to statistically account for numerous
observations between 50 MHz and 1 GHz. This improvement is achieved without
resorting to unrealistic values for the strength of the Galactic magnetic
field, which, as measured through Zeeman splitting in the diffuse ISM and
rotation measures of pulsars, averages less than 10 $\mu$G (Heiles & Troland,
2005; Sobey et al., 2019). This research represents a notable advancement in
the modeling of Galactic diffuse emission, offering crucial insights to
prepare for the arrival of the SKA.

The paper is organized as follows: in Sect. 2, we introduce the basic
formalism of synchrotron emission and review several observations of
$\beta_{\rm obs}$ at low frequencies reported in the literature; in Sect. 3,
we present our $j_{e}$ model; in Sect. 4, we derive an all-sky map of the
magnetic-field strength and validate our methodology with MHD simulations. The
paper concludes with a discussion (Sect. 5) and a summary (Sect. 6).

##  2 Observations of the synchrotron spectral index

In this section, we first introduce the formalism linking $j_{e}$ and
$\beta_{\rm obs}$ (Sect. 2.1) and then summarize the observed values of
$\beta_{\rm obs}$ below 1 GHz (Sect. 2.2). We provide a list of observational
results that does not aim to be comprehensive but rather representative of
what can be found in the literature.

###  2.1 From first principles to brightness temperature

Galactic synchrotron diffuse emission arises from the interaction of
interstellar magnetic fields with CRe (e.g., Ginzburg & Syrovatskii, 1964;
Rybicki & Lightman, 1979). These relativistic particles gyrate about the field
lines emitting non-thermal radiation dependent on their energy ($E$) and the
magnetic-field strength perpendicular to the sightline
(${B}_{\perp}=|\@vec{B}_{\perp}|$). The total synchrotron emissivity in units
of power per unit volume, frequency ($\nu$), and solid angle is the sum of the
emissivities linearly polarized along ($\varepsilon_{\nu,\parallel}$) and
across ($\varepsilon_{\nu,\perp}$) $\@vec{B}_{\perp}$:

$\displaystyle\varepsilon_{\nu,\|}(\@vec{r})$ $\displaystyle=$
$\displaystyle\int_{m_{e}c^{2}}^{\infty}\frac{j_{e}(E)}{\varv_{e}}P_{\nu,\|}^{\rm
em}(E,{B_{\perp}}(\@vec{r}))\,{\rm d}E$ (1)
$\displaystyle\varepsilon_{\nu,\perp}(\@vec{r})$ $\displaystyle=$
$\displaystyle\int_{m_{e}c^{2}}^{\infty}\frac{j_{e}(E)}{\varv_{e}}P_{\nu,\perp}^{\rm
em}(E,{B_{\perp}}(\@vec{r}))\,{\rm d}E,$

where $\varv_{e}$ is the CRe velocity, $m_{e}$ is their mass, $c$ is the speed
of light, $\@vec{r}$ is the position, and $P^{\rm
em}_{\nu,\perp\leavevmode\nobreak\ \rm{\rm or}\leavevmode\nobreak\ \parallel}$
is the power per unit frequency for the two polarizations (see also equation 2
in P21 and Longair 2011). By integrating Eq. 1 along the line of sight (LOS),
and accounting for the telescope beam ($\Theta_{\rm T}$), one obtains the
synchrotron brightness temperature as

$\frac{T_{\rm b}(\nu)}{[{\rm K}]}=\frac{2\ln{(2)}\leavevmode\nobreak\
c^{2}10^{-23}}{\pi k_{\rm B}\theta_{\rm T}^{2}\nu^{2}}\left[\Theta_{\rm
T}\circledast\int_{\rm LOS}(\varepsilon_{\nu,\|}+\varepsilon_{\nu,\perp}){\rm
d}r\right],$ (2)

where $k_{\rm B}$ is the Boltzmann constant and the $\circledast$-symbol
represents the convolution with a Gaussian-approximated $\Theta_{\rm T}$ given
its full width half maximum (FWHM), $\theta_{\rm T}$, in units of radians.
Observationally, $T_{\rm b}(\nu)$ mostly results in a power-law function of
$\nu$ with spectral index $\beta_{\rm obs}$ (e.g., Reich & Reich, 1988a, b;
Roger et al., 1999). From Eqs. 1 and 2, and as already shown by P21, the value
of $\beta_{\rm obs}$ is not expected to be constant but rather to vary
depending on the MHD properties of the intervening ISM and on $j_{e}$.
Notably, it has long been known that energy losses may contribute to modifying
the synchrotron spectral index, as illustrated by several observations at low
frequencies in the radio band (e.g., Bridle, 1967; Sironi, 1974; Webster,
1975) as well as by diffusion and convection models of CRe propagation (e.g.,
Strong, 1977).

In this work we exploit these dependencies and invert the argument using
ancillary observations of $\beta_{\rm obs}$ to constrain $j_{e}$ and
${B}_{\perp}$. In the following, we will distinguish the observed values of
the spectral index, $\beta_{\rm obs}$, from those derived from models, denoted
as $\beta_{\rm syn}$.

###  2.2 Observed spectral indices below 1 GHz

Figure 1: Normalized histograms of the observed spectral indices, $\beta_{\rm
obs}$, between 45 MHz and 408 MHz from G11 as a function of Galactic latitude.
The 50th, the 16th, and the 84th percentiles of the orange histogram are
marked in black-dashed and -solid lines, respectively.

Figure 2: Our models for $s_{\rm CRe}$ obtained with Eq. 3 and their impact on
the modelled synchrotron spectral index $\beta_{\rm syn}$. Colours correspond
to the same models in all panels. Left panel: $j_{e}$ from models A, B, C, and
D with parameters listed in Table LABEL:tab:models that fit the black data
points from Voyager ($E<10^{8}$ eV, Cummings et al., 2016; Stone et al., 2019)
and Fermi LAT ($E>10^{10}$ eV, Ackermann et al., 2010). Central panel: $s_{\rm
CRe}$ values corresponding to the same models as in the left panel. Right
panel: values of $\beta_{\rm syn}$ for the same models as in the other panels
with the observed spectral indices from G11 in hatches – the 16th, the 50th,
and the 84th percentiles are marked as in Fig. 1. Models for $s_{\rm CRe}$ and
$\beta_{\rm syn}$ using $j_{e}$ from Orlando (2018) and Padovani et al. (2018)
are also shown in dashed-purple and -teal lines in the central and right
panels, respectively.

Measurements of the synchrotron spectral index from diffuse emission in the
Galaxy generally vary between $-2$ and $-3$ depending on the observed
frequency range, the sky coverage, and the FWHM of the instrument. Most
$\beta_{\rm obs}$ at frequencies below 408 MHz have been recently measured
with single dipole antennas designed for low-frequency cosmological
experiments, such as the Experiment to Detect the Global EoR Signature (EDGES,
Bowman et al., 2008) and the Large-Aperture Experiment to Detect the Dark Ages
(LEDA, Price et al., 2018). These experiments measured $-2.59<\beta_{\rm
obs}<-2.54$ between 50 MHz and 100 MHz (Mozdzen et al., 2019),
$-2.62<\beta_{\rm obs}<-2.60$ between 90 MHz and 190 MHz (Mozdzen et al.,
2017), and $-2.55<\beta_{\rm obs}<-2.45$ between 50 MHz and 87 MHz (Spinelli
et al., 2021). While these measurements integrate the synchrotron flux over
the full accessible sky, higher angular resolution observations of $\beta_{\rm
obs}$ are also available. In particular, in this work we refer to one of the
all-sky $\beta_{\rm obs}$ maps of reference for low-frequency studies derived
by Guzmán et al. (2011), hereafter G11, between 45 MHz and 408 MHz using data
from Alvarez et al. (1997), Maeda et al. (1999), and Haslam et al.
(1982)222The map can be downloaded at the following link:
http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/525/A138. This map, with a FWHM
of 5 degrees, shows large variations of $\beta_{\rm obs}$ ranging between
$-2.7$ and $-2.1$. In Fig. 1, we show the normalized histograms of these
$\beta_{\rm obs}$ as a function of Galactic latitude ($b$). Most of the
shallow values of $\beta_{\rm obs}$ are confined to the Galactic plane,
$|b|<5^{\circ}$ (compare the blue and orange curves), possibly due to free-
free absorption of radio photons along the LOS (e.g., Stanislavsky et al.,
2023). However, the distribution of the values of $\beta_{\rm obs}$ at
intermediate to high Galactic latitude above $|b|>5^{\circ}$ is stable with a
rather normal shape (see the green and orange curves). The 50th, 16th, and
84th percentiles of the orange histogram are $-2.53$, $-2.57$, and $-2.47$,
respectively.

More recent all-sky maps at the degree-scale level, depicting $\beta_{\rm
obs}$ values from 35 MHz to 80 MHz, are available from the LWA (Dowell et al.,
2017). These maps exhibit good agreement with those presented by G11, albeit
displaying slightly shallower $\beta_{\rm obs}$ values at intermediate
Galactic latitudes. The authors did not extensively investigate these
differences, attributing them to (i) variations in frequency coverage and (ii)
residual systematic effects present in their data.

Remarkably, $\beta_{\rm obs}$ maps of external galaxies, such as NGC 5775 and
NGC 891, have been produced at low frequency with LOFAR and the Very Large
Array (Mulcahy et al., 2018; Heald et al., 2020). We notice that, in both
cases, the maps appear consistent with the $\beta_{\rm obs}$ map of G11
outside of their galactic halo.

At higher frequencies, above 408 MHz, the $\beta_{\rm obs}$ values for diffuse
Galactic emission typically converge towards $-3$. A recent example of this
comes from the analysis of MeerKAT data (with FWHM of approximately 1 degree)
between 971 MHz and 1075 MHz at intermediate Galactic latitudes, which yielded
an average $\beta_{\rm obs}$ value of $-2.76\pm 0.15$ (Irfan et al., 2022, and
references therein).

##  3 Modeling the CRe energy flux

The ability of models to reproduce the observed synchrotron sky brightness and
$\beta_{\rm obs}$ strongly relies on the choice of $j_{e}$ and its
approximations (see Eq. 1). As discussed in P21, most theoretical models and
numerical simulations of Galactic synchrotron emission assume fixed power laws
for the energy distribution of CRe depending on the frequency range (e.g., Sun
et al., 2008; Waelkens et al., 2009; Sathyanarayana Rao et al., 2017; Reissl
et al., 2019; Wang et al., 2020). However, this assumption may lead to
significant biases in the interpretation of synchrotron data, notably
concerning the spatial variations of the spectral index both at low and high
frequency. The energy dependence of the slope of $j_{e}$ between 100 MeV and
50 GeV should be taken into account.

Some effort has been done in this direction, as shown for instance in the
works by Padovani et al. (2018) and Orlando (2018). The former proposed an
analytical four-parameter model that reproduces the power-law behaviour of
$j_{e}$ measured at low and high energies (see Sect. 1). The latter is based
on a combination of multi-frequency observations, Voyager measurements, and
the CRe propagation code called GALPROP (e.g. Porter et al., 2022, for a
recent release). However, as shown with dashed lines in the right panel of
Fig. 2, both cases struggle in reproducing the range of $\beta_{\rm obs}$
described in Sect. 2.2 (see hatches from G11)333We note that Orlando (2018)
developed their $j_{e}$ model including the radio data employed by G11.
However, since their post-processing was not the same, in particular
concerning their methods of component separation (i.e., thermal and non-
thermal components), the respective $\beta_{\rm obs}$ are different. We choose
to consider the work of G11, as it is widely referred to in the radio
community at low frequency. It is beyond our scopes to discuss at length the
inconsistencies between these two works. without resorting to unlikely large
values for the average $B_{\perp}$ along the LOS (hereafter, $\langle
B_{\perp}\rangle$) of the order of 20 $\mu$G.

In this context, we propose a new analytical model for $j_{e}$ tailored to
capture its variation in the energy gap encompassing 1 GeV. We introduce a
four-parameter model that describes the logarithmic derivative of $j_{e}$ with
respect to $E$ ($s_{\rm CRe}$) as follows

$s_{\rm
CRe}=\frac{d\log_{10}{j_{e}}}{d\log_{10}{E}}=\frac{s_{1}-s_{2}}{1+(E/E_{0})^{\gamma}}+s_{2},$
(3)

where $\gamma=1/\log{(\Delta E/[\rm eV])}$, and $\Delta E$ represents the
energy range around $E_{0}$ at which the transition between the low-energy and
high-energy slopes, $s_{1}$ and $s_{2}$, occurs. For comparison, the $s_{\rm
CRe}$ curves from Orlando (2018) and Padovani et al. (2018) are displayed in
the central panel of Fig. 2.

Instead of fitting the full set of independent data at once (from those of
$j_{e}$ to the radio ones) to illustrate the impact of our parameters, in Fig.
2, we show four realizations of our model that correspond to models A to D in
Table LABEL:tab:models.

Table 1: Model parameters entering Eq. 3 and shown in Figs. 2 and 6. model | $s_{1}$ | $s_{2}$ | $E_{0}$ | $\gamma$ | region
---|---|---|---|---|---
|  |  | [GeV] |  |
A | $-1.8$ | $-3.1$ | $3.2$ | $45$ | diffuse ISM
B | $-1.8$ | $-3.2$ | $3.2$ | $2.2$ | diffuse ISM
C | $-1.5$ | $-3.1$ | $1$ | $1.1$ | diffuse ISM
D | $-1.3$ | $-3.1$ | $0.6$ | $14$ | diffuse ISM
E | $-1.2$ | $-3.1$ | $3.2$ | $45$ | Orion-Taurus ridge

From left to right, models are shown with the same color for $j_{e}$, $s_{\rm
CRe}$, and $\beta_{\rm syn}$, respectively. Despite all fitting the Voyager
and Fermi-LAT data in the left panel, they behave differently in the 1-GeV
energy window, as depicted in the central panel. This energy range determines
the observable brightness and the corresponding $\beta_{\rm syn}$, as
illustrated in the right panel. Model A, characterized by a sharp change in
slope of $j_{e}$ at about 3 GeV (not distinguishable from a step function in
$s_{\rm CRe}$), is the type of model that converges in both energy and
frequency for values of $\langle B_{\perp}\rangle$ below 10 $\mu$G. Our
approach is purely phenomenological and, on average, suggests the existence of
two power laws of CRe contributing to synchrotron emission.

Figure 3: Modeled synchrotron spectral indices, $\beta_{\rm syn}$, reproducing
various observational results of $\beta_{\rm obs}$ (in colours) with model A
(solid lines), Orlando (2018) (dashed lines), and Padovani et al. (2018)
(dotted lines). Data are the following: between 50 - 100 MHz (red), and 90 -
190 MHz (blue) using EDGES (Mozdzen et al., 2017, 2019); between 45 - 408 MHz
using a set of multiple observations (G11); between 971 - 1075 MHz using
MeerKAT (Irfan et al., 2022); between 50 and 87 MHz using LEDA (Spinelli et
al., 2021).

Using model A as a reference, in Fig. 3, we also compare our estimates of
$\beta_{\rm syn}$ with the models of Padovani et al. (2018) and Orlando (2018)
for a number of observations below 1 GHz (see Sect. 2.2). As for the range of
$\beta_{\rm obs}$ from G11, we show that model A better reproduces $\beta_{\rm
obs}$ both at low and high frequency for realistic values of $\langle
B_{\perp}\rangle<10\,\mu$G. We notice that most single-dipole-antenna
experiments produce $\beta_{\rm obs}$ values consistent with smaller $\langle
B_{\perp}\rangle$ than other instruments. This is possibly due to averaging
Galactic diffuse emission and $B_{\perp}$ over larger sky areas.

##  4 Estimate of the magnetic-field strength

Figure 4: Map of the plane-of-the-sky magnetic-field strength averaged along
the line of sight, $\langle B_{\perp}\rangle$, derived with $\beta_{\rm obs}$
from G11 (see orange histogram in Fig. 1) and assuming a spatially uniform
$j_{e}$ described by model A. We label one possible instrumental pattern and a
few known radio structures: the Orion-Taurus ridge (Ori-Tau ridge, Bracco et
al., 2023), the North Polar Spur (NPS), Loop III, and Loop IIIs (e.g., Vidal
et al., 2015). A Galactic coordinate grid centered in the Galactic center with
steps of $30^{\circ}$ in latitude and longitude is overlaid.

Assuming model A for $j_{e}$, in this section we derive a map of $\langle
B_{\perp}\rangle$ from the $\beta_{\rm obs}$ map of G11 (Sect. 4.1). We
compare these results with synthetic data of synchrotron diffuse emission
(Sect. 4.2).

###  4.1 All-sky $\langle B_{\perp}\rangle$ map from $\beta_{\rm obs}$ between
45 and 408 MHz

In the previous sections, we have demonstrated that model A is capable of
statistically reproducing the correct values of $\beta_{\rm obs}$ observed at
intermediate-to-high Galactic latitudes, both by fitting the available $j_{e}$
data and for values of $\langle B_{\perp}\rangle$ lower than 10 $\mu$G. With
this model, we utilize the corresponding functional form of $\beta_{\rm syn}$,
which is dependent on $\langle B_{\perp}\rangle$, as depicted in red in the
right panel of Fig. 2, to interpolate the $\beta_{\rm obs}$ map presented in
G11 and infer the respective map of $\langle B_{\perp}\rangle$.

The resulting all-sky map is displayed in Fig. 4, shown with a linear scale
between 1 $\mu$G and 30 $\mu$G. A circular region around the North Celestial
Pole is masked as in G11. We notice that most of the sky above $|b|>5^{\circ}$
is consistent with values of $\langle B_{\perp}\rangle$ below 15 $\mu$G, as
expected from Sect. 3. This suggests that, over most of the diffuse sky,
considering a uniform spatial distribution of $j_{e}$ to describe the energy
flux of CRe is a reasonable assumption.

Figure 5: Estimate of $\langle B_{\perp}\rangle$ using synthetic intensity
data at 45 and 408 MHz obtained with model A from one MHD simulation (Bellomi
et al., 2020) already implemented in P21 to study diffuse Galactic synchrotron
emission. The initial mean magnetic-field strength of the simulation is 4
$\mu$G. The value of $\langle B_{\perp}\rangle$ is estimated as in Fig. 4.
Left panel: the ratio between the estimated $\langle B_{\perp}\rangle$ and the
true median value of the plane-of-the-sky magnetic field strength per pixel,
$\tilde{B}_{\perp}$, is shown as a function of $\tilde{B}_{\perp}$ for the
three integration axes. Right panel: for the green case shown in the left
panel, given the injected $j_{e}$-shape of model A in the simulations,
$\langle B_{\perp}\rangle$ is also estimated using presently existing models,
i.e. Orlando (2018) (in purple) and Padovani et al. (2018) (in teal). In all
panels levels show the (0.5, 1, 1.5, 2)-$\sigma$ equivalent contours of the
corresponding density of points.

However, a few regions stand out of the map. On one hand some striped patterns
in the West may be associated with instrumental artifacts present in the G11
map. On the other hand, most of the notable structures correspond to well-
known arcs in the radio sky, such as the Orion-Taurus ridge (Ori-Tau ridge,
Planck 2015 results. XXV., 2016; Bracco et al., 2023), the North Polar Spur
(NPS), and Loops III-IIIs (e.g., Vidal et al., 2015). In Sect. 5, we discuss
possible reasons that could explain their special appearance on the map with
values of $\langle B_{\perp}\rangle$ significantly larger than 15 $\mu$G.
Finally, some residual contamination of free-free absorption in the proximity
of the Galactic plane could bias the estimate of $\langle B_{\perp}\rangle$
toward large values (e.g., between Galactic longitudes,
$-30^{\circ}<l<+30^{\circ}$).

###  4.2 Validation of the $\langle B_{\perp}\rangle$ map with synthetic data

We validate the methodology used to derive the $\langle B_{\perp}\rangle$ map
by using synthetic data of synchrotron emission at 45 MHz and 408 MHz obtained
from one state-of-the-art MHD numerical simulation of the multiphase and
magnetized ISM (Bellomi et al., 2020). In summary, this numerical model,
generated using the RAMSES code (Teyssier, 2002; Fromang et al., 2006),
simulates a 50-pc cube with 128 pixels per side of warm-diffuse gas (with
initial temperatures and volume densities of $T=8000$ K and $n_{\rm H}=1.5$
cm-3, respectively), where cold and dense structures ($T<100$ K and $n_{\rm
H}>100$ cm-3) form due to thermal instability and hydro-magnetic turbulence.
The initial magnetic field considered in the simulation aligns with typical
values in the diffuse ISM (4 $\mu$G) and is directed along one of the main
axes of the cube (hereafter, the $x$ axis).

This model has already been utilized in P21 to simulate radio emission below
200 MHz. We follow the same recipe and derive $T_{\rm b}$ by means of Eqs. 1
and 2, where $B_{\perp}$ varies depending on the choice of the main axis of
the cube considered for the integration of Eq. 2, and $j_{e}$ is described by
model A. Given the simulated $T_{\rm b}$ maps, we measure $\beta_{\rm syn}$
between 45 MHz and 408 MHz for each LOS as done with real observations. Since
we do not apply any Gaussian beam to the synthetic data, the choice of
$\theta_{\rm T}$ does not affect the estimate of $\beta_{\rm syn}$. From the
$\beta_{\rm syn}$ maps, we interpolate back the value of $\langle
B_{\perp}\rangle$ using model A (see Sect. 4.1) and compare it with the true
median value of the simulated magnetic field along each sight line
($\tilde{B}_{\perp}$).

In the left panel of Fig. 5, we show the ratio between $\langle
B_{\perp}\rangle$ and $\tilde{B}_{\perp}$ as a function of the latter one. We
explore three cases depending on the integration axis, namely one along (plane
of the sky – POS – $yz$) and two across (POS, $xy$ and $xz$) the mean
magnetic-field direction.

Our methodology generally results in an excess bias toward the actual value of
$\tilde{B}_{\perp}$, which becomes more pronounced, up to a factor of 2, as
the magnetic-field strength decreases. Consequently, the most severe bias
arises when integrating the cube along the mean direction of the magnetic
field. In the other two cases our estimates are highly consistent (within
1-$\sigma$) with the true value of the magnetic field.

In the right panel of Fig. 5, for the POS equal to the $xy$ plane, we test the
bias on the estimate of $\tilde{B}_{\perp}$ in the case of misinterpreting the
energy flux of CRe. While always considering model A in the synthetic data, we
interpolate back $\langle B_{\perp}\rangle$ referring to the models of
Padovani et al. (2018) and Orlando (2018) in teal and purple colours,
respectively. As expected from Fig. 2, both cases significantly overestimate
$\tilde{B}_{\perp}$ compared to model A, reaching excess biases of more than
one order of magnitude in the case of Padovani et al. (2018).

##  5 Discussion

We have introduced a new analytical model for $j_{e}$ that depends on a
minimal set of parameters, making it easily usable for modeling the diffuse
synchrotron emission of the Galaxy. Nevertheless, it is important to clarify
that our phenomenological model does not represent the injection spectrum of
CRe at the emission source; instead, it corresponds to the effective CRe
spectrum after propagation and attenuation processes.

We have shown that, over most of the sky, $j_{e}$ can be considered rather
spatially uniform (at least at the angular scales probed by this work). A set
of model parameters (see model A in Table LABEL:tab:models) can be chosen such
that the values of $\langle B_{\perp}\rangle$ needed to reproduce available
$\beta_{\rm obs}$ at low frequencies are broadly consistent with the magnetic-
field strength measured in the diffuse ISM below 10 $\mu$G (e.g., Heiles &
Troland, 2005; Sobey et al., 2019). This is displayed in the all-sky map of
$\langle B_{\perp}\rangle$ shown in Fig. 4. However, a few structures
coinciding with known arcs in the radio sky seem to distinguish themselves
from the average trend with anomalously large values of $\langle
B_{\perp}\rangle$. The origin of these structures is highly debated (e.g.,
Sofue, 2015; Kataoka et al., 2021; Panopoulou et al., 2021). They could
represent the relic of several supernova remnants, which may affect the energy
distribution of CRe in the presence of stronger and compressed magnetic fields
(e.g., Morlino & Celli, 2021). As an example, in Fig. 6, we show the same as
in Fig. 2 but for the $\beta_{\rm obs}$ values measured by Bracco et al.
(2023) in one radio arc, namely the Ori-Tau ridge. In this case, it is not
possible to reproduce the observations between $-2.3$ and $-2.2$ using model A
unless resorting to $\langle B_{\perp}\rangle$ values of hundreds of $\mu$G.
Conversely, the range of $\beta_{\rm obs}$ in the Ori-Tau ridge can be
accounted for with more standard values of $\langle B_{\perp}\rangle$ using
our model E in Table LABEL:tab:models. Model E adjusts the CRe population at
low energy (below 1 GeV), departing from the local ISM conditions probed by
the Voyager spacecraft.

Figure 6: Same as in the left and right panels of Fig. 2 but for the observed
values of $\beta_{\rm obs}$ in the Orion-Taurus ridge in hatches (Bracco et
al., 2023). Models A and E for $\beta_{\rm syn}$ (right) and $j_{e}$ (left)
are overlaid in red and cyan lines, respectively (see Table LABEL:tab:models).

This example illustrates the versatility of our model and its ability to be
tailored to various astrophysical scenarios. Such adaptability is crucial for
optimizing existing modelling of synchrotron emission (e.g., Bracco et al.,
2022), especially as new discoveries at low frequencies continue to emerge.
For instance, our model could help us guide the interpretation of recent and
unprecedented findings made with the LOFAR telescope, one of the main SKA
precursors below 200 MHz (e.g., Zaroubi et al., 2015; Van Eck et al., 2017;
Jelić et al., 2018; Van Eck et al., 2019; Bracco et al., 2020; Turić et al.,
2021). These observations highlighted the unexpected correlation between
synchrotron polarization, affected by Faraday rotation (Ferrière et al.,
2021), with tracers of the neutral interstellar gas, clearing the path to a
new perspective on structure formation in the multiphase and magnetized
Galactic ISM.

##  6 Summary

Cosmic-ray electrons (CRe) with energies around 1 GeV are expected to be the
primary contributors to the diffuse synchrotron radiation observed in the
Galaxy at radio frequencies below 1 GHz. However, determining their exact
energy flux ($j_{e}$) has proven challenging due to solar modulation occurring
between 100 MeV and 10 GeV. Previous attempts to reconcile the radio non-
thermal emission with $j_{e}$ within this energy range have succeeded at the
cost of requiring unexpectedly large values for the average Galactic magnetic-
field strength perpendicular to the sight-line ($\langle B_{\perp}\rangle$).

In this work, we have introduced a new analytical model for $j_{e}$, capable
of fitting ancillary data on the CRe energy flux (i.e., Voyager and Fermi LAT)
while meeting observational constraints on the spectral index of the
brightness temperature at low frequencies, denoted as $\beta_{\rm obs}$. The
model achieves this while maintaining values of $\langle B_{\perp}\rangle$
consistent with existing measurements in the diffuse ISM below 10 $\mu$G
(e.g., Crutcher, 2012). We used the $\beta_{\rm obs}$ map obtained between 45
MHz and 408 MHz by G11, found the best model that reproduces these
observations, and validated its performance against various measurements of
$\beta_{\rm obs}$ in the literature between 50 MHz and 1 GHz.

Our fiducial model (model A in Table LABEL:tab:models) allowed us to derive,
for the first time, an all-sky map of $\langle B_{\perp}\rangle$ that we
tested against synthetic data of synchrotron emission produced from state-of-
the-art MHD simulations of the multiphase and magnetized ISM.

We found that, up to angular scales of a degree level, a uniform spatial
distribution of $j_{e}$ can be assumed for diffuse emission at intermediate-
to-high Galactic latitudes. Only a few regions in the sky, which correspond to
known arc-like structures in the radio band, are at variance with the average
trend. We demonstrated the versatility of our model to account for these
anomalous regions in the case of the Orion-Taurus ridge (Bracco et al., 2023).

This work represents a useful step forward in modeling Galactic diffuse
emission at low frequency given the adaptability of our data-driven four-
parameter model of $j_{e}$. It marks an important milestone in the
interpretation of current low-frequency observations of the diffuse ISM, such
as those from the LOFAR telescope, and for future analyses with the SKA coming
online.

###### Acknowledgements.

We express our gratitude to the referee, Andy Strong, whose comments improved
the quality of this manuscript. We are grateful to Elena Bellomi, Benjamin
Godard, and Patrick Hennebelle for providing us with the MHD simulations that
were used to generate synthetic data of synchrotron emission. We are thankful
to Marta Spinelli for discussions on low-frequency observations of the
synchrotron spectral index. AB acknowledges financial support from the INAF
initiative “IAF Astronomy Fellowships in Italy” (grant name MEGASKAT). AB also
thanks Elio for setting the pace with his memorable notes ”Olélé moliba
makasi”. In the analysis we made use of astropy (Astropy Collaboration et al.,
2018), scipy (Virtanen et al., 2020), and numpy (Harris et al., 2020).

## References

  * Ackermann et al. (2010) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2010, Phys. Rev. D, 82, 092004
  * Adriani et al. (2011) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, ApJ, 742, 102
  * Aguilar et al. (2014) Aguilar, M., Aisa, D., Alvino, A., et al. 2014, Phys. Rev. Lett., 113, 121102
  * Alvarez et al. (1997) Alvarez, H., Aparici, J., May, J., & Olmos, F. 1997, A&AS, 124, 205
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bellomi et al. (2020) Bellomi, E., Godard, B., Hennebelle, P., et al. 2020, A&A, 643, A36
  * Bowman et al. (2008) Bowman, J. D., Rogers, A. E. E., & Hewitt, J. N. 2008, ApJ, 676, 1
  * Bracco et al. (2020) Bracco, A., Jelić, V., Marchal, A., et al. 2020, A&A, 644, L3
  * Bracco et al. (2022) Bracco, A., Ntormousi, E., Jelić, V., et al. 2022, A&A, 663, A37
  * Bracco et al. (2023) Bracco, A., Padovani, M., & Soler, J. D. 2023, A&A, 677, L11
  * Bridle (1967) Bridle, A. H. 1967, MNRAS, 136, 219
  * Crutcher (2012) Crutcher, R. M. 2012, ARA&A, 50, 29
  * Cummings et al. (2016) Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18
  * Dewdney et al. (2009) Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proceedings, 97, 1482
  * Dowell et al. (2017) Dowell, J., Taylor, G. B., Schinzel, F. K., Kassim, N. E., & Stovall, K. 2017, MNRAS, 469, 4537
  * Ferrière et al. (2021) Ferrière, K., West, J. L., & Jaffe, T. R. 2021, MNRAS, 507, 4968
  * Fromang et al. (2006) Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
  * Gabici (2022) Gabici, S. 2022, A&A Rev., 30, 4
  * Ginzburg & Syrovatskii (1964) Ginzburg, V. L. & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays
  * Guzmán et al. (2011) Guzmán, A. E., May, J., Alvarez, H., & Maeda, K. 2011, A&A, 525, A138 (G11)
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Haslam et al. (1982) Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1
  * Heald et al. (2020) Heald, G., Mao, S. A., Vacca, V., et al. 2020, Galaxies, 8, 53
  * Heiles & Troland (2005) Heiles, C. & Troland, T. H. 2005, ApJ, 624, 773
  * Irfan et al. (2022) Irfan, M. O., Bull, P., Santos, M. G., et al. 2022, MNRAS, 509, 4923
  * Jelić et al. (2018) Jelić, V., Prelogović, D., Haverkorn, M., Remeijn, J., & Klindžić, D. 2018, A&A, 615, L3
  * Jones et al. (2018) Jones, M. E., Taylor, A. C., Aich, M., et al. 2018, MNRAS, 480, 3224
  * Kataoka et al. (2021) Kataoka, J., Yamamoto, M., Nakamura, Y., et al. 2021, ApJ, 908, 14
  * Longair (2011) Longair, M. S. 2011, High Energy Astrophysics
  * Maeda et al. (1999) Maeda, K., Alvarez, H., Aparici, J., May, J., & Reich, P. 1999, A&AS, 140, 145
  * Morlino & Celli (2021) Morlino, G. & Celli, S. 2021, MNRAS, 508, 6142
  * Mozdzen et al. (2017) Mozdzen, T. J., Bowman, J. D., Monsalve, R. A., & Rogers, A. E. E. 2017, MNRAS, 464, 4995
  * Mozdzen et al. (2019) Mozdzen, T. J., Mahesh, N., Monsalve, R. A., Rogers, A. E. E., & Bowman, J. D. 2019, MNRAS, 483, 4411
  * Mulcahy et al. (2018) Mulcahy, D. D., Horneffer, A., Beck, R., et al. 2018, A&A, 615, A98
  * Orlando (2018) Orlando, E. 2018, MNRAS, 475, 2724
  * Padovani et al. (2021) Padovani, M., Bracco, A., Jelić, V., Galli, D., & Bellomi, E. 2021, A&A, 651, A116 (P21)
  * Padovani & Galli (2018) Padovani, M. & Galli, D. 2018, A&A, 620, L4
  * Padovani et al. (2018) Padovani, M., Galli, D., Ivlev, A. V., Caselli, P., & Ferrara, A. 2018, A&A, 619, A144
  * Panopoulou et al. (2021) Panopoulou, G. V., Dickinson, C., Readhead, A. C. S., Pearson, T. J., & Peel, M. W. 2021, ApJ, 922, 210
  * Planck 2015 results. XXV. (2016) Planck 2015 results. XXV. 2016, A&A, 594, A25
  * Ponnada et al. (2024) Ponnada, S. B., Panopoulou, G. V., Butsky, I. S., et al. 2024, MNRAS, 527, 11707
  * Porter et al. (2022) Porter, T. A., Jóhannesson, G., & Moskalenko, I. V. 2022, ApJS, 262, 30
  * Price et al. (2018) Price, D. C., Greenhill, L. J., Fialkov, A., et al. 2018, MNRAS, 478, 4193
  * Reich & Reich (1988a) Reich, P. & Reich, W. 1988a, A&AS, 74, 7
  * Reich & Reich (1988b) Reich, P. & Reich, W. 1988b, A&A, 196, 211
  * Reissl et al. (2019) Reissl, S., Brauer, R., Klessen, R. S., & Pellegrini, E. W. 2019, ApJ, 885, 15
  * Rockstroh & Webber (1978) Rockstroh, J. M. & Webber, W. R. 1978, ApJ, 224, 677
  * Roger et al. (1999) Roger, R. S., Costain, C. H., Landecker, T. L., & Swerdlyk, C. M. 1999, A&AS, 137, 7
  * Rybicki & Lightman (1979) Rybicki, G. B. & Lightman, A. P. 1979, Radiative processes in astrophysics
  * Sathyanarayana Rao et al. (2017) Sathyanarayana Rao, M., Subrahmanyan, R., Udaya Shankar, N., & Chluba, J. 2017, ApJ, 840, 33
  * Sironi (1974) Sironi, G. 1974, MNRAS, 166, 345
  * Sobey et al. (2019) Sobey, C., Bilous, A. V., Grießmeier, J. M., et al. 2019, MNRAS, 484, 3646
  * Sofue (2015) Sofue, Y. 2015, MNRAS, 447, 3824
  * Spinelli et al. (2021) Spinelli, M., Bernardi, G., Garsden, H., et al. 2021, MNRAS, 505, 1575
  * Stanislavsky et al. (2023) Stanislavsky, L. A., Bubnov, I. N., Konovalenko, A. A., Stanislavsky, A. A., & Yerin, S. N. 2023, A&A, 670, A157
  * Stone et al. (2019) Stone, E. C., Cummings, A. C., Heikkila, B. C., & Lal, N. 2019, Nature Astronomy, 3, 1013
  * Strong (1977) Strong, A. W. 1977, MNRAS, 181, 311
  * Strong et al. (2000) Strong, A. W., Moskalenko, I. V., & Reimer, O. 2000, ApJ, 537, 763
  * Strong & Wolfendale (1978) Strong, A. W. & Wolfendale, A. W. 1978, Journal of Physics G Nuclear Physics, 4, 1793
  * Sun et al. (2008) Sun, X. H., Reich, W., Waelkens, A., & Enßlin, T. A. 2008, A&A, 477, 573
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Turić et al. (2021) Turić, L., Jelić, V., Jaspers, R., et al. 2021, A&A, 654, A5
  * Unger & Farrar (2023) Unger, M. & Farrar, G. R. 2023, arXiv e-prints, arXiv:2312.13273
  * Van Eck et al. (2019) Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2019, A&A, 623, A71
  * Van Eck et al. (2017) Van Eck, C. L., Haverkorn, M., Alves, M. I. R., et al. 2017, A&A, 597, A98
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * Vidal et al. (2015) Vidal, M., Dickinson, C., Davies, R. D., & Leahy, J. P. 2015, MNRAS, 452, 656
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Waelkens et al. (2009) Waelkens, A., Jaffe, T., Reinecke, M., Kitaura, F. S., & Enßlin, T. A. 2009, A&A, 495, 697
  * Wang et al. (2020) Wang, J., Jaffe, T. R., Enßlin, T. A., et al. 2020, ApJS, 247, 18
  * Webster (1975) Webster, A. 1975, MNRAS, 171, 243
  * Zarka et al. (2012) Zarka, P., Girard, J. N., Tagger, M., & Denis, L. 2012, in SF2A-2012: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. S. Boissier, P. de Laverny, N. Nardetto, R. Samadi, D. Valls-Gabaud, & H. Wozniak, 687–694
  * Zaroubi et al. (2015) Zaroubi, S., Jelic, V., de Bruyn, A. G., et al. 2015, MNRAS, 454, L46


# Optimal Scheduling of Battery Storage Systems in the Swedish Multi-FCR
Market Incorporating Battery Degradation and Technical Requirements ††thanks:
The work presented in this paper is financially supported by the following
projects: (i) ECom4Future FIWARE Driven Energy Communities for the Future -
received funding from the Clean Energy Transition Partnership under Grant
Agreement No. CETP22FP-2023-00316, and (ii) Implementation of the Vehicle to
Grid Services in Sweden (PAVE) - received funding from Vinnova (Sweden’s
innovation agency) within the FFI program under the Authority’s Dnr
2023-00785.
The authors are with the Department of Electrical Engineering, Chalmers
University of Technology, 412 96 Gothenburg, Sweden. (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Nima Mirzaei Alavijeh, Rahmat Khezri, Senior Member, IEEE, Mohammadreza
Mazidi, David Steen, and Le Anh Tuan, Member, IEEE

###### Abstract

This paper develops a novel mixed-integer linear programming (MILP) model for
optimal participation of battery energy storage systems (BESSs) in the Swedish
frequency containment reserve (FCR) markets. The developed model aims to
maximize the battery owner’s potential profit by considering battery
degradation and participation in multiple FCR markets, i.e., FCR in normal
operation (FCR-N), and FCR in disturbances (FCR-D) for up- and down-
regulations. Accordingly, a precise formulation of a detailed battery
degradation model and adherence to the technical requirements of the Swedish
FCR markets are incorporated into the developed model. To achieve more
practical results, simulations are conducted based on one minute time step
realistic data for the whole year 2022. The results show a potential profit of
k€ 708 for a 1MW/1MWh BESS by participating in multi-FCR market. Analyzing the
impact of considering degradation in the optimization problem has shown that
the annual battery aging cost could decrease by 5%-29% without a significant
effect on profit. The proposed model can be practically used by flexibility
asset owners to achieve profitable and sustainable operation strategies that
reduce battery degradation.

###### Index Terms:

Ancillary service, battery energy storage system, battery degradation,
frequency containment reserve, technical market requirements, optimal
scheduling

## Nomenclature

###  -A Acronyms

BESS

Battery energy storage systems

DA

Day-ahead

EV

Electric vehicle

FCR

Frequency containment reserve

FCR-N

Frequency containment reserve in normal operation

FCR-D

Frequency containment reserve in disturbances

LER

Limited energy reservoirs

MILP

Mixed-integer linear programming

OM

Operation and maintenance

SoE

State of energy

SoC

State of charge

###  -B Indices and Sets

$d\in\mathcal{D}$

Index and set of days in a year

$t\in\mathcal{T}_{d}$

Index and set of time steps in a day

$h\in\mathcal{H}_{d}$

Index and set of hours in a day

$\Xi$

Set of decision variables

###  -C Parameters and Functions

$\mathcal{A}h_{t,d}$

Ah throughput of BESS at time step $t$ of day $d$ [Ah]

$\mathcal{B}^{\mathrm{BAT}}$

Net present value of battery [€]

$\mathcal{B}_{d}^{\mathrm{DEG}}$

Cost of battery degradation on day $d$ [€]

$\mathcal{B}_{t,d}^{\mathrm{CAL}}$, $\mathcal{B}_{t,d}^{\mathrm{CYC}}$

Calendar and cycle aging cost of battery degradation on time step $t$ of day
$d$ [€]

$C_{d}^{\mathrm{DA}}$

Cost of BESS in day-ahead spot market on day $d$ [€]

$C^{\mathrm{REP}}$

Cost of replacement of battery cells [€]

$C^{\mathrm{REP}}$

Cost of operation and maintenance of BESS [€]

$e_{t/h}^{\mathrm{DR,X}}$, $e_{t/h}^{\mathrm{UR,X}}$

Energy content: Per unit activated energy for down/up regulation for FCR-X
service during $t$ or $h$ [h]

$EOL$

The percentage of retained capacity at end of life of BESS [%]

$F_{d}$

Total profit of day $d$ [€]

$f_{t,d}$

Grid frequency at time step $t$ of day $d$ [Hz]

$f_{n}$

Nominal grid frequency [Hz]

$f_{\mathrm{min}}^{\mathrm{N/D}}$

Minimum grid frequency for FCR-N or FCR-D activation [Hz]

$f_{\mathrm{max}}^{\mathrm{N/D}}$

Grid frequency for full activation of FCR-N or FCR-D services [Hz]

$\mathcal{I}_{t,d}$

C-rate of BESS at time step $t$ of day $d$

$i$

Interest rate [%]

$L$

Life time of BESS [Years]

$\mathcal{K}_{t}$

Temperature at time step $t$ [K]

$\underline{P}$, $\overline{P}$

Minimum and maximum charger power (MW)

$\underline{p}^{\Theta,X}$

Minimum bid size for FCR-X market, X $\in\\{\mathrm{N,DU,DD}\\}$ [MW]

$r^{sv}$

Ratio of salvage cost to replacement cost

$R_{d}^{\mathrm{DA}}$

Revenue of BESS in day-ahead spot market on day $d$ [€]

$R_{d}^{\mathrm{FCR}}$

Revenue from FCR market participation on day $d$ [€]

$R_{h,d}^{\mathrm{X}}$

Revenue from FCR-X market on hour $h$ of day $d$, X
$\in\\{\mathrm{N,DU,DD}\\}$ [€]

$\underline{S}$, $\overline{S}$

Minimum and maximum state-of-energy

$\mathcal{S}_{t,d}$

State of energy of BESS in time step $t$ of day $d$ [MWh]

$\mathcal{U}_{t,d}^{\mathrm{CAL}}$, $\mathcal{U}_{t,d}^{\mathrm{CYC}}$

Calendar and cycle aging of BESS on time step $t$ of day $d$ [%]

$\eta^{\mathrm{CH}},\eta^{\mathrm{DS}}$

Charging, discharging efficiency

$\rho^{\mathrm{spot}}_{h}$

Day-ahead spot market price at hour $h$ [€/MWh]

$\rho^{\mathrm{tax}}$

Fixed electricity tax [€/MWh]

$\rho^{\mathrm{g}}$

Grid utilization tariff [€/MWh]

$\rho^{\mathrm{X}}_{h}$

Capacity reimbursement price for FCR-X market at hour $h$, X
$\in\\{\mathrm{N,DU,DD}\\}$ [€/MW]

$\rho_{h}^{\mathrm{UR}}$, $\rho_{h}^{\mathrm{DR}}$

Up/down regulation prices at hour $h$ [€/MWh]

$\Delta p_{t,d}$

Deviation from baseline power due to the activation FCR services [MW]

$\Delta t$

The length of the simulation time step [s]

###  -D Variables

$b_{t,d}^{\mathrm{CH/DS}}$

Binary variable for charge/discharge power of BESS in time step $t$ of day $d$

$b_{h,d}^{\mathrm{CH/DS},\mathrm{bl}}$

Binary variable for baseline charge/discharge power of BESS at hour $h$ of day
$d$

$b_{h,d}^{X}$

Binary variable for participation in FCR-X market at hour $h$ of day $d$, X
$\in\\{\mathrm{N,DU,DD}\\}$

$p_{h,d}^{\Theta,\mathrm{X}}$

Power bid on FCR-X market on hour $h$ of day $d$, X
$\in\\{\mathrm{N,DU,DD}\\}$ [MW]

$p_{t,d}^{\mathrm{CH/DS}}$

Charge/discharge power of BESS in time step $t$ of day $d$ [MW]

$p_{h,d}^{\mathrm{CH/DS},\mathrm{bl}}$

Baseline (reference) charge/discharge power of BESS at hour $h$ of day $d$
[MW]

$p_{t,d}^{\mathrm{ND}}$, $p_{t,d}^{\mathrm{NU}}$

Activated up or down regulation for FCR-N service by BESS in time step $t$ of
day $d$ [MW]

$p_{t,d}^{\mathrm{DD}}$, $p_{t,d}^{\mathrm{DD}}$

Activated up or down regulation for FCR-D services by BESS in time step $t$ of
day $d$ [MW]

##  I Introduction

###  I-A Motivation and Aim

Frequency containment reserve (FCR) markets are one of the market-based tools
used by the Nordic electricity system operators for regulating grid frequency.
FCR services comprises three products: FCR in normal operation (FCR-N) and FCR
in disturbances (FCR-D) for up and down regulations, to manage small and
larger frequency fluctuations, respectively [1].

A recent market trend is the increasing participation of distributed
flexibility resources. Flexibility asset owners involved in providing FCR must
maintain fast response times and engage in shorter, but more frequent
activation events. Accordingly, battery energy storage systems (BESSs) are
considered ideal FCR providers because, as inverter-coupled resources, they
can deliver fast power response despite their limited energy capacity [2]. It
was shown by Hu et al. [3] that BESSs generate the highest revenue by
providing various FCRs in most European electricity markets. However, to
ensure a profitable scheduling strategy for BESSs, it is crucial to consider
the market requirements and precisely model the costs associated with battery
degradation.

This paper aims to develop an optimal scheduling model for the efficient
participation of BESSs in the Swedish energy and FCR markets. The objective is
to maximize profit while considering battery degradation accurately and
meeting all technical requirements of the Swedish markets.

###  I-B Literature Review

By analyzing the literature on the application of BESSs in the frequency
regulation of power systems, two distinct modeling approaches can be
identified, i.e., control-based and market-based approaches. The first
approach focuses on the optimal control of BESSs to effectively regulate
frequency. In [4], a distributed control strategy was proposed for multiple
BESSs to regulate frequency in the power system with high penetration of
renewable generation. In [5], an online frequency regulation strategy based on
the Lyapunov optimization technique was proposed. In [6], it was proposed that
the frequency regulation signal can be divided into a slow component, directed
to synchronous generators, and a fast component, directed to BESSs. Authors of
[7] proposed a robust control strategy to manage distributed BESSs for
frequency regulation. While these papers focused on developing control
strategies for BESSs to regulate frequency, they did not address the
optimization of profits that can be obtained by participating in the ancillary
service market.

Figure 1: Motivations, contributions, and benefits of the conducted study.

The participation of BESSs in ancillary service markets is investigated in the
second approach. An example of such an approach is [8], where a bidding
strategy and online control methodology for BESSs was proposed to enable
participation in both the day-ahead energy market and the FCR-N market. While
the battery degradation was modeled, participation in the FCR-D up/down market
was not considered. Furthermore, the proposed model was nonlinear which adds
computational burden and makes finding the global optimum solution
challenging. In [9], a two-stage stochastic optimization model was developed
for the optimal bidding of BESSs in both energy and ancillary service markets.
While the model effectively addressed participation in the multi-FCR market,
it fell short in full accounting for market requirements. Additionally, the
accuracy of battery degradation modeling in the proposed framework was deemed
simplistic. In [10], electric vehicles (EVs) were aggregated as a BESS, and a
two-level optimization model was proposed to determine their optimal
participation in energy and ancillary service markets. In the proposed model,
sole participation in the FCR-N market was considered, while factors such as
battery degradation and market requirements were neglected. In [11], robust
and stochastic methods were used to model the participation of aggregated EVs
in both the energy and ancillary service markets. In the developed model, a
constant annual activation ratio was calculated for FCR using real historical
data and incorporated into the supply-demand balance equation as a random
variable. While their developed model considered multi-FCR market
participation, it was not developed with the technical requirements of Nordic
FCR markets, and the battery degradation is not a detailed model. The authors
of [12] have introduced an optimal model for sizing and scheduling BESSs to
maximize income from participating in the energy and FCR-N markets. However,
the proposed model employed a simplified degradation model for the battery and
overlooked the FCR market requirements.

To ensure a profitable strategy for BESSs, the battery degradation should be
precisely modeled. This is because, during FCR-N activation, BESS is used as a
power buffer, alternating between charging and discharging cycles. This can
potentially result in substantial energy exchanges between the BESS and the
grid, thereby accelerating battery degradation in addition to high conversion
losses [13]. The impact of providing FCRs on the degradation of BESSs has been
examined in [14]. The results indicate that BESSs have been subjected to many
cycles with a low average cycle depth providing FCRs, leading to a dominance
of calendar aging over cyclic aging. In [15], degradation has been identified
as the most critical factor affecting the profitability of BESSs. Although
energy exchange between the BESS and the grid is significantly lower in FCR-D
up/down markets compared to the FCR-N market, FCR-N’s bid price is higher than
the FCR-D up/down in most of the time. Therefore, to maximize the
profitability of BESSs, it is essential to consider participation in multi-FCR
markets while accounting for degradation. As demonstrated in [13], selecting
the best combination of FCRs in each hour of BESS scheduling can increase the
profitability by 22% compared to delivering pure FCR-N. However, modeling the
participation of BESSs in multi-FCR markets is a burdensome and demanding task
nowadays due to the complicated and essential requirements of FCR markets. As
investigated in [16], neglecting market requirements may lead to situations
where balancing the BESS is insufficient to provide the FCRs, resulting in
penalties.

After a detailed review of the literature and thorough focus on relevant
research works, owing to the best of our knowledge, it has been found that no
work has been dedicated to providing an optimal scheduling model for BESSs in
the multi-FCR market considering battery degradation and meeting all technical
requirements of the Swedish markets.

###  I-C Approach and Contributions

To fill the research gaps, this study develops a new optimal BESS scheduling
model for multi-market FCR participation by considering all technical
requirements of the Swedish market and comprehensive battery degradation. The
proposed model is formulated with respect to a potential stacked hourly
participation in FCR markets while optimizing the purchased power from the DA
spot market acting as an optimal reference power for delivering FCR services.
From a practical point of view, the developed model in this study can be
utilized by flexibility asset owners in obtaining more sustainable and yet
profitable operation strategies for BESSs. The main contributions of this
study can be summarized as follows:

  * •

Formulating a MILP model for optimal scheduling of BESSs in DA spot and multi-
FCRs markets which can be solved efficiently through optimization software.

  * •

Incorporating all technical requirements of the Swedish FCR market in the
developed model and solving the optimal scheduling problem based on one minute
realistic data for one year.

  * •

Applying a detailed model of battery degradation in the developed model to
ensure the profitability of the model.

The proposed model enables time-efficient analyses on how the maximum
potential profit and participation strategy of BESS are impacted by input
parameters and model settings such as the inclusion of battery degradation,
and the inclusion of FCR markets technical requirements. By proposing this
model, the authors are pursuing the following objectives:

  * •

Modeling charge and discharge powers as variables linked to the bid size and
the activation criteria, enabling degradation to be considered in the
optimization.

  * •

Including endurance and power requirements from the latest technical
requirements for the participation of limited energy reservoirs in FCR
markets.

  * •

Enabling the possibility of stacked participation in multiple FCR markets at
each hour.

  * •

Optimizing the reference power considering both the DA spot and FCR markets.

The motivations, contributions, and benefits of the developed model in this
study are presented in Figure 1.

###  I-D Paper Organization

The rest of the paper is laid out as follows. The frequency containment
reserve and market requirements are presented in Section II. The model
formulation considering the battery degradation is described in Section III.
Case study and input data are elaborated in Section IV. The results and
discussion are given in Section V and conclusions are provided in Section VI.

##  II Swedish Frequency Containment Reserves

In Sweden, there are three services for FCR market: FCR-N, FCR-D up, and FCR-D
down. There are differences between regulations and technical requirements for
these FCR services. The main aspects of various FCR services are summarized in
Table I. For a realistic model formulation, applying these aspects is
essential.

TABLE I: Comparison of market and technical requirements of FCR services [17, 18] |  FCR-N |  FCR-D up |  FCR-D down
---|---|---|---
Market |  |  |
Trading time |  Day-ahead |  Day-ahead |  Day-ahead
Remuneration |  Capacity and energy |  Capacity |  Capacity
Min bid size |  0.1 MW |  0.1 MW |  0.1 MW
Technical requirements |  |  |
Required power upwards |  $1.34\cdot p_{h}^{\mathrm{\Theta,N}}$ |  $p_{h}^{\mathrm{\Theta,DU}}$ |  $0.2\cdot p_{h}^{\mathrm{\Theta,DD}}$
Required power downwards |  $1.34\cdot p_{h}^{\mathrm{\Theta,N}}$ |  $0.2\cdot p_{h}^{\mathrm{\Theta,DU}}$ |  $p_{h}^{\mathrm{\Theta,DD}}$
Required endurance upwards |  1h |  $\frac{1}{3}$h |  0
Required endurance downwards |  1h |  0 |  $\frac{1}{3}$h
Activation |  Automatic (Figure 2) |  Automatic (Figure 2) |  Automatic (Figure 2)

From a market perspective, the difference between the FCR services concerns
their remuneration. All three services are remunerated for their cleared bid
capacity based on pay-as-cleared prices [18]. In addition, FCR-N is
remunerated also based on the activated energy. This is because of the
relatively higher energy throughput in FCR-N services. Up-regulation energy is
compensated with up-regulation prices, while the service provider is charged
with down-regulation prices for the down-regulation energy.

To calculate realistic bidding and profit, technical requirements from the FCR
markets need to be considered. ENTSO-e is the European association for the
cooperation of transmission system operators that has provided the technical
requirements [17] for limited energy reservoirs (LER) such as EVs and
batteries. The power requirements are formulated based on the accepted bid
$P_{h}^{\mathrm{\Theta,X}}$. The endurance requirements indicate how long the
service provider must be able to provide the accepted bid at a full activation
scenario.

The activation of FCR services is based on the specified droop curves shown in
Figure 2. These curves should be strictly followed for the accepted bid size.
There are exceptions for LERs in case their state of charge (SOC) goes below
or above certain limits specified in [17]. In case these limits are passed,
the resource should deviate from the droop curves to bring back the SOC to the
acceptable ranges. In such cases, the service provider will not be remunerated
since it could not provide the service. These exceptions are not considered in
this study and instead a safety margin is implemented for the SOC by its lower
and upper bounds.

Figure 2: FCR activation power based on grid frequency: (a) FCR-N, (b) FCR-D

Considering the market and technical requirements in the model is essential
for a realistic estimation of bidding behavior and profitability of services.
In Section III, a novel MILP formulation is proposed for implementing these
requirements.

##  III Model Formulation

The developed optimization model in this study is presented in this section.
The model formulation is presented in different subsections representing the
objective function, and constraints related to the battery’s physical limits,
FCR participation rules and requirements, and battery degradation model.

###  III-A Objective Function

The objective function of the BESS scheduling problem is to maximize the daily
profit given by:

$\displaystyle f$ $\displaystyle=\max_{\Xi}\quad F_{d}$ (1a)
$\displaystyle=\max_{\Xi}\quad(R_{d}^{\mathrm{DA}}+R_{d}^{\mathrm{FCR}}-C_{d}^{\mathrm{DA}}-C_{d}^{\mathrm{DEG}}),\quad\forall
d\in\mathcal{D},$
$\displaystyle\Xi=\\{p_{h,d}^{\mathrm{CH,bl}},p_{h,d}^{\mathrm{DS,bl}},p_{t,d}^{\mathrm{CH}},p_{t,d}^{\mathrm{DS}},$
(1b) $\displaystyle\quad\quad\quad
p_{h,d}^{\mathrm{\Theta,N}},p_{h,d}^{\mathrm{\Theta,DU}},p_{h,d}^{\mathrm{\Theta,DD}}\mid
h\in\mathcal{H}_{d},t\in\mathcal{T}_{d}\\},\quad\forall d\in\mathcal{D}$

where $\mathcal{D}$ is the set of all the days in a year, $\mathcal{T}_{d}$ is
the set of all the timestamps in day d with a length of lower than one hour,
$\mathcal{H}_{d}$ is the set of all hours in day $d$. $R_{d}^{\mathrm{DA}}$ is
the revenue from DA spot market through discharging, $R_{d}^{\mathrm{FCR}}$ is
the revenue from FCR participation, $C_{d}^{\mathrm{DA}}$ is the cost in DA
spot market through charging, and $C_{d}^{\mathrm{DEG}}$ is the cost for
battery degradation. The decision variables of the optimization problem are
the hourly purchased and sold power in the DA spot market
($p_{h,d}^{\mathrm{CH,bl}},p_{h,d}^{\mathrm{DS,bl}}$) that are acting as the
reference point (baseline) for the FCR services, charging and discharging
power of the battery ($p_{t,d}^{\mathrm{CH}}$, $p_{t,d}^{\mathrm{DS}}$), power
bid to participate in FCR-N, FCR-D up and FCR-D down
($p_{h,d}^{\Theta,\mathrm{N}}$, $p_{h,d}^{\Theta,\mathrm{DU}}$,
$p_{h,d}^{\Theta,\mathrm{DD}}$). To avoid repetition, the indices related to
the day are not written in the rest of the equations.

The DA spot market revenue is formulated based on hourly spot market price
($\rho_{h}^{\mathrm{spot}}$), and discharging power
($p_{h}^{\mathrm{DS,bl}}$).

$R^{\mathrm{DA}}=\sum_{h\in\mathcal{H}}{p_{h}^{\mathrm{DS,bl}}\cdot(\rho_{h}^{\mathrm{spot}}+\rho^{\mathrm{tax}})}$
(2)

The revenue for multi-market FCR participation contains three parts of
revenues for FCR-N, down-regulation FCR-D, and up-regulation FCR-D.

$R^{\mathrm{FCR}}=\sum_{h\in\mathcal{H}}{R_{h}^{\mathrm{N}}+R_{h}^{\mathrm{DD}}+R_{h}^{\mathrm{DU}}}$
(3)

The FCR-N revenue includes three terms. The first term is the capacity
reimbursement for power bid. The capacity is reimbursed according to the FCR-N
market price ($\rho_{h}^{\mathrm{N}}$) and the amount of power bid
($p_{h}^{\Theta,\mathrm{N}}$). The second and third terms are the energy
compensation for upregulation and downregulation during the FCR-N
participation. These terms are calculated based on the concept of energy
content introduced by [13]. Energy content represents the per unit activated
energy for up or down ward regulations within a specific timestep $t$ or a
specific hour $h$. It is pre-calculated based on the frequency and droop
settings. Multiplying energy content ($e_{h}^{\mathrm{UR/DR,N}}$), FCR-N power
bid, and up and down regulation prices ($\rho_{h}^{\mathrm{UR/DR}}$) at each
hour will provides the hourly energy compensation for FCR-N services

$R_{h}^{\mathrm{N}}=p_{h}^{\Theta,\mathrm{N}}\cdot\left(\rho_{h}^{\mathrm{N}}\
+\rho_{h}^{\mathrm{UR}}\cdot e_{h}^{\mathrm{UR,N}}-\rho_{h}^{\mathrm{DR}}\cdot
e_{h}^{\mathrm{DR,N}}\right),\quad\forall h$ (4)

The revenues for FCR-D down and up are only based on capacity reimbursement
according to FCR-D up and down prices and the power bid.

$R_{h}^{\mathrm{DD}}=\rho_{h}^{\mathrm{DD}}\
p_{h}^{\Theta,\mathrm{DD}},\quad\forall\ h$ (5)
$R_{h}^{\mathrm{DU}}=\rho_{h}^{\mathrm{DU}}\
p_{h}^{\Theta,\mathrm{DU}},\quad\forall\ h$ (6)

The expenditure for charging in the DA spot market is formulated based on the
hourly spot price, grid utilization cost ($\rho^{g}$), electricity tax
($\rho^{\mathrm{tax}}$), and charging power.

$C^{\mathrm{DA}}=\sum_{h\in\mathcal{H}}{p_{h}^{\mathrm{CH,bl}}\left(\rho_{h}^{\mathrm{spot}}+\rho^{g}+\rho^{\mathrm{tax}}\right)}$
(7)

The cost of battery degradation is the sum of cycle aging cost
($\mathcal{B}_{t}^{\mathrm{CYC}}$) and calendar aging cost
($\mathcal{B}_{t}^{\mathrm{CAL}}$):

$C^{\mathrm{DEG}}=\sum_{t\in
T}\left(\mathcal{B}_{t}^{\mathrm{CAL}}+\mathcal{B}_{t}^{\mathrm{CYC}}\right)\
$ (8)

###  III-B Battery Constraints

The battery model contains several constraints related to power and state of
energy (SoE) in the optimization model. The charging and discharging power are
either zero or between $\underline{P}$ and $\bar{P}$. In fact, the power
cannot exceed the charger’s maximum power limit at any time. On the other
hand, due to low efficiency in low power, the chargers do not allow charging
or discharging with a power lower than the minimum limit.

$p_{h}^{\mathrm{CH,bl}},\
p_{h}^{\mathrm{DS,bl}},P_{t}^{\mathrm{CH}},P_{t}^{\mathrm{DS}}\in\left\\{0\right\\}\cup\left[\underline{P},\overline{P}\,\right],\quad\forall
h$ (9)

To avoid simultaneous charging and discharge, and to implement (9), binary
variables ($b_{t}^{\mathrm{CH}}$, $b_{t}^{\mathrm{DS}}$,
$b_{h}^{\mathrm{CH,bl}}$, $b_{h}^{\mathrm{DS,bl}}$) are considered. The sum of
the first pair and the sum of the second of binary variables should be less
than or equal to one at their respective timestep.

$b_{t}^{\mathrm{CH/DS}}\cdot\underline{P}\leq p_{t}^{\mathrm{CH/DS}}\leq
b_{t}^{\mathrm{CH/DS}}\cdot\overline{P},\quad\forall t$ (10)
$b_{h}^{\mathrm{CH/DS,bl}}\cdot\underline{P}\leq p_{h}^{\mathrm{CH/DS,bl}}\leq
b_{h}^{\mathrm{CH/DS,bl}}\cdot\overline{P},\quad\forall h$ (11)
$b_{t}^{\mathrm{CH}}+b_{t}^{\mathrm{DS}}\leq 1,\quad\forall t$ (12)
$b_{h}^{\mathrm{CH,bl}}+b_{h}^{\mathrm{DS,bl}}\leq 1,\quad\forall h$ (13)

The SoE ($\mathcal{S}_{t}$) is formulated based on the SoE from the previous
timestamp ($\mathcal{S}_{t-1,d}$), charging and discharging baseline power,
and the up and down regulation energy corresponding to FCR services.

$\displaystyle\mathcal{S}_{t}=\mathcal{S}_{t-1}$
$\displaystyle+\left(P_{h}^{\mathrm{CH,bl}}\cdot\eta^{\mathrm{CH}}-P_{h}^{\mathrm{DS,bl}}/\eta^{\mathrm{DS}}\right)\cdot\Delta
t$ (14)
$\displaystyle+P_{h_{t}}^{\mathrm{\Theta,N}}\cdot\left(e_{t}^{\mathrm{DR,N}}-e_{t}^{\mathrm{UR,N}}\right)$
$\displaystyle+P_{h_{t}}^{\mathrm{\Theta,DD}}\cdot e_{t}^{\mathrm{DR,DD}}$
$\displaystyle-P_{h_{t}}^{\mathrm{\Theta,DU}}\cdot
e_{t}^{\mathrm{UR,DU}},\quad\forall t$
$\underline{\mathcal{S}}\leq\mathcal{S}_{t}\leq\overline{\mathcal{S}},\quad\forall
t$ (15)

where $\eta^{\mathrm{CH}}$ and $\eta^{\mathrm{DS}}$ are the charging and
discharging efficiencies, and $e_{t}^{\mathrm{UR/DR,X}}$ is the energy content
of up- or down-regulation activation for different FCR markets. Hour $h_{t}$
represents the corresponding hour to timestamp $t$.

###  III-C FCR participation constraints

Several constraints need to be considered to satisfy the technical
requirements and physical limitations related to stacked multi-FCR market
participation. The activated power for delivering FCR-N service is formulated
based on the frequency deviations of the power system. It follows the droop
control shown in Figure 2.

$\displaystyle
p_{t}^{\mathrm{ND}}=\begin{cases}p_{h}^{\mathrm{\Theta,N}},&f_{t}>f_{\max}^{\mathrm{N}}\\\
p_{h}^{\mathrm{\Theta,N}}\left(\frac{f_{t}-f_{n}}{f_{\max}^{\mathrm{N}}-f_{n}}\right),&f_{n}\leq
f_{t}\leq f_{\max}^{\mathrm{N}}\\\ 0,&f_{t}<f_{n}\end{cases},\quad\forall t,$
(16) $\displaystyle p_{t}^{\mathrm{NU}}=\begin{cases}0,&f_{t}>f_{n}\\\
p_{h}^{\mathrm{\Theta,N}}\left(\frac{f_{t}-f_{n}}{f_{\min}^{\mathrm{N}}-f_{n}}\right),&f_{\min}^{\mathrm{N}}\leq
f_{t}\leq f_{n}\\\
p_{h}^{\mathrm{\Theta,N}},&f_{t}<f_{\min}^{\mathrm{N}}\end{cases},\quad\forall
t,$ (17)

where $f_{t}$ is the measured frequency of the system in Hz in timestamp $t$,
$f_{n}$ is the nominal frequency (50 Hz), $f_{\max}^{\mathrm{N}}$ and
$f_{\min}^{\mathrm{N}}$ are the maximum and minimum ranges of frequency for
FCR-N participation, 50.1 Hz and 49.9 Hz, respectively.

The activated power for delivering FCR-D markets is also formulated based on
the frequency deviations of the power system according to Figure 2.

$\displaystyle
p_{t}^{\mathrm{DD}}=\begin{cases}p_{h}^{\mathrm{\Theta,DD}},&f_{t}>f_{\max}^{\mathrm{D}}\\\
p_{h}^{\mathrm{\Theta,DD}}\left(\frac{f_{t}-f_{\max}^{\mathrm{N}}}{f_{\max}^{\mathrm{D}}-f_{\max}^{\mathrm{N}}}\right),&f_{max}^{\mathrm{N}}\leq
f_{t}\leq f_{\max}^{\mathrm{D}}\\\
0,&f_{t}<f_{\max}^{\mathrm{N}}\end{cases},\quad\forall\ t,$ (18)
$\displaystyle
p_{t}^{\mathrm{DU}}=\begin{cases}0,&f_{t}>f_{\min}^{\mathrm{N}}\\\
p_{h}^{\mathrm{\Theta,DU}}\left(\frac{f_{t}-f_{\min}^{\mathrm{N}}}{f_{\min}^{\mathrm{D}}-f_{\min}^{\mathrm{N}}}\right),&f_{\min}^{\mathrm{D}}\leq
f_{t}\leq f_{\min}^{\mathrm{N}}\\\ p_{h}^{\mathrm{\Theta,DU}},&f_{t}<\
f_{\min}^{\mathrm{D}}\end{cases},\quad\forall\ t,$ (19)

where $f_{\max}^{\mathrm{D}}$ and $f_{\min}^{\mathrm{D}}$ are the maximum and
minimum ranges of frequency for FCR-D participation, 50.5 Hz and 49.5 Hz,
respectively.

The activated power for FCR services should be delivered with respect to the
original plan of the battery. This is shown in (20) where the original plan is
the baseline power. The baseline is the traded energy in the DA spot market in
the context of this model. The deviation from the baseline is shown by $\Delta
p_{t}$. A positive $\Delta p_{t}$ corresponds to down-regulation and negative
to up-regulation because load convention is used for the power.

$\displaystyle\Delta p_{t}$ $\displaystyle=p_{t}-p_{h_{t}}^{\mathrm{bl}}$ (20)
$\displaystyle=\left(p_{t}^{\mathrm{CH}}-p_{t}^{\mathrm{DS}}\right)-\left(p_{h_{t}}^{\mathrm{CH,bl}}-p_{h_{t}}^{\mathrm{DS,bl}}\right),\quad\forall
t$

When bidding in the FCR markets, $\Delta p_{t}$ should be strictly equal to
the specified power levels from the droop controls:

$\Delta
P_{t}=P_{t}^{\mathrm{ND}}+P_{t}^{\mathrm{DD}}-P_{t}^{\mathrm{NU}}-P_{t}^{\mathrm{DU}},\quad\forall
t.$ (21)

The lower and upper bounds related to the bids need to be defined. Power bids
for each market should not exceed double the maximum power of the battery. The
upper bound is double the maximum power because the battery can, for example,
go from full charge to full discharge power. The minimum power is also
considered for each market participation ($\underline{p}^{\mathrm{\Theta,N}}$,
$\underline{p}^{\mathrm{\Theta,DU}}$, and
$\underline{p}^{\mathrm{\Theta,DD}}$). Binary variables ($b_{t}^{\mathrm{N}}$,
$b_{t}^{\mathrm{DD}}$, and $b_{t}^{\mathrm{DU}}$) are added to implement the
minimum power bid:

$b_{h}^{\mathrm{N}}\cdot\underline{p}^{\mathrm{\Theta,N}}\leq
p_{h}^{\mathrm{\Theta,N}}\leq b_{h}^{\mathrm{N}}\cdot\overline{P},\quad\forall
h,$ (22) $b_{h}^{\mathrm{DD}}\cdot\underline{p}^{\mathrm{\Theta,DD}}\leq
p_{h}^{\mathrm{\Theta,DD}}\leq b_{h}^{\mathrm{DD}}\cdot
2\overline{P},\quad\forall h,$ (23)
$b_{h}^{\mathrm{DU}}\cdot\underline{p}^{\mathrm{\Theta,DU}}\leq
p_{h}^{\mathrm{\Theta,DU}}\leq b_{h}^{\mathrm{DU}}\cdot
2\overline{P},\quad\forall h.$ (24)

To have results closer to real life, technical requirements for the FCR
markets (Table I) need to be considered. The power requirements are
transferred to the physical limitations of the battery in Figure 3 using the
maximum charge and discharge capacity ($\overline{P}$) and the baseline
($p_{h}^{\mathrm{bl}}$) which is the reference point for providing the up- or
down-regulation. The figure can be formulated as constraints provided in (25)
and (26).

$1.34p_{h}^{\mathrm{\Theta,N}}+p_{h}^{\mathrm{\Theta,DU}}+0.2p_{h}^{\mathrm{\Theta,DD}}\leq\overline{P}+p_{h}^{\mathrm{bl}},\quad\forall
h$ (25)
$1.34p_{h}^{\mathrm{\Theta,N}}+p_{h}^{\mathrm{\Theta,DD}}+0.2p_{h}^{\mathrm{\Theta,DU}}\leq\overline{P}-p_{h}^{\mathrm{bl}},\quad\forall
h$ (26)

Figure 3: Limits on FCR capacity bids considering the baseline and ENTSOe
requirements on power. N: $P_{h}^{\mathrm{\Theta,N}}$, DU:
$P_{h}^{\mathrm{\Theta,DU}}$, DD: $P_{h}^{\mathrm{\Theta,DD}}$

The energy requirements (i.e., endurance) can be implemented for the worst-
case scenarios for SoC at each hour. The potential extreme scenarios are
visualized in Figure 4 for up- and downwards activation when the reference
power is at charging. A similar analysis can be done for a reference power at
discharging. These most critical points are marked with red circles and can
occur in the following scenarios:

  1. 1.

If the service is not activated at all and the battery is locked to follow its
baseline: equation (27)

  2. 2.

If both FCR-N and FCR-D services are activated for 20 minutes: equation (28)

  3. 3.

If FCR-N and FCR-D services are simultaneously activated for 20 minutes and
FCR-N for the remaining 40 minutes within the hour: equation (29)

$\underline{\mathcal{S}}\leq S_{h-\Delta
t}+P_{h}^{\mathrm{bl}}\leq\overline{\mathcal{S}},\quad\forall h$ (27)
$\displaystyle\underline{\mathcal{S}}\leq\mathcal{S}_{h-\Delta
t}+\left(p_{h}^{\mathrm{bl}}+p_{h}^{\mathrm{\Theta,N}}+p_{h}^{\mathrm{\Theta,DD}}\right)\cdot\frac{1}{3}\leq\overline{\mathcal{S}},\quad\forall
h$ (28a) $\displaystyle\underline{\mathcal{S}}\leq\mathcal{S}_{h-\Delta
t}+\left(p_{h}^{\mathrm{bl}}-p_{h}^{\mathrm{\Theta,N}}-p_{h}^{\mathrm{\Theta,DU}}\right)\cdot\frac{1}{3}\leq\overline{\mathcal{S}},\quad\forall
h$ (28b) $\displaystyle\underline{\mathcal{S}}\leq\mathcal{S}_{h-\Delta
t}+p_{h}^{\mathrm{bl}}+p_{h}^{\mathrm{\Theta,N}}+\frac{1}{3}\cdot
p_{h}^{\mathrm{\Theta,DD}}\leq\overline{\mathcal{S}},\quad\forall h$ (29a)
$\displaystyle\underline{\mathcal{S}}\leq\mathcal{S}_{h-\Delta
t}+p_{h}^{\mathrm{bl}}-p_{h}^{\mathrm{\Theta,N}}-\frac{1}{3}\cdot
p_{h}^{\mathrm{\Theta,DU}}\leq\overline{\mathcal{S}},\quad\forall h$ (29b)

Figure 4: Examples of critical activation cases for endurance requirement when
reference power is at charging. (a) downwards regulation, (b) upwards
regulation

###  III-D Battery degradation

The costs of cycle aging and calendar aging of the battery are calculated
based on the degradation of the battery due to cycle aging
($\mathcal{U}_{t,d}^{\mathrm{CYC}}$) and calendar aging
($\mathcal{U}_{t,d}^{\mathrm{CAL}}$), percentage of retained capacity at end-
of-life (EOL), and the net present value of the battery
($\mathcal{B}^{\mathrm{BAT}}$).

$\mathcal{B}_{t}^{\mathrm{CAL}}=\mathcal{B}^{\mathrm{BAT}}\frac{\mathcal{U}_{t}^{\mathrm{CAL}}(\%)}{100\%-\mathrm{EOL}(\%)},\quad\forall\
t,$ (30)
$\mathcal{B}_{t}^{\mathrm{CYC}}=\mathcal{B}^{\mathrm{BAT}}\frac{\mathcal{U}_{t,d}^{\mathrm{CYC}}(\%)}{100\%-\mathrm{EOL}(\%)},\quad\forall\
t.$ (31)

Typically, the percentage of retained capacity at the EOL of the battery is
considered as 80% [19]. The net present value of the battery is calculated
based on the lifetime of the battery (L), interest rate (i), cost of
replacement ($C^{\mathrm{REP}}$), cost of operation and maintenance
($C^{\mathrm{OM}}$), and the ratio of salvage cost to replacement cost
($r^{\mathrm{sv}}$).

$\mathcal{B}^{\mathrm{BAT}}=\left(1-r^{\mathrm{sv}}\right)\
C^{\mathrm{REP}}\frac{1}{\left(1+i\right)^{L}}+C^{\mathrm{OM}}\frac{\left(1+i\right)^{L}-1}{{\alpha\left(1+i\right)}^{L}}$
(32)

The battery degradation model is taken from [20] where an experimental model
was developed for a Lithium-Nickel-Manganese-Cobalt + Lithium-Manganese oxide
battery. The authors derived nonlinear empirical aging models for the same
battery where the calendar aging depends on SOC, temperature
($\mathcal{K}_{t}$), and battery age in days ($\mathcal{Q}$), and the cycle
aging depends on temperature, C-rate ($\mathcal{I}_{t}$), and Ah-throughput
(${\mathcal{A}h}_{t}$).

$\mathcal{U}_{t}^{\mathrm{CAL}}=\mathcal{F}\left(\mathcal{S}_{t},\mathcal{K}_{t},\mathcal{Q}\right),\quad\forall
t$ (33)
$\mathcal{U}_{t}^{\mathrm{CYC}}=\mathcal{F}\left(\mathcal{K}_{t},\mathcal{I}_{t},{\mathcal{A}h}_{t}\right),\quad\forall
t.$ (34)

The model of calendar aging is a nonlinear function:

$\mathcal{U}_{t}^{\mathrm{CAL}}=G\left(\mathcal{S}_{t}\right)\
e^{\left(\frac{E_{a}}{R\
\mathcal{K}_{t}}\right)}\mathcal{Q}^{0.5},\quad\forall t,$ (35) $\displaystyle
G\left(\mathcal{S}_{t}\right)=\begin{cases}a_{1}{\mathcal{S}_{t}}^{2}+a_{2}\mathcal{S}_{t}+a_{3},&0\leq\mathcal{S}_{t}\leq
0.5Q\\\
b_{1}{\mathcal{S}_{t}}^{2}+b_{2}\mathcal{S}_{t}+b_{3},&0.5Q<\mathcal{S}_{t}\leq
0.7Q\\\
c_{1}{\mathcal{S}_{t}}^{2}+c_{2}\mathcal{S}_{t}+c_{3},&0.7Q<\mathcal{S}_{t}\leq
1.0Q\end{cases},\quad\forall t.$ (36)

Equation (35) is included in the problem by piece-wise linearization using
binaries and auxiliary variables over the three SoC spans in (36).

The cycle aging model [20] is also non-linear:

$\mathcal{U}_{t}^{\mathrm{CYC}}=\left(q_{1}{\mathcal{K}_{t}}^{2}+q_{2}\mathcal{K}_{t}+q_{3}\right)e^{q_{4}\mathcal{I}_{t}}{Ah}_{t},\quad\forall
t.$ (37)

The cycle aging is linearized over the operation range of the battery, i.e.,
power: 0.1-1.0 MW, SoC: 10%-90%, and assuming a relatively constant
temperature of 20${}^{\circ}\mathrm{C}$ because stationary batteries are
commonly kept in containers with regulated temperatures and effective battery
management systems.

Detail of mathematical linearization for calendar aging and cycle aging models
is presented in [19] and not repeated in this paper. The utilized aging
parameters are presented in Table II.

##  IV Simulated Cases and Input Data

The impact of degradation is assessed for five market participation cases:

  * •

Case w/o FCR: do not participate in FCR,

  * •

Case FCR-N: participate only in FCR-N,

  * •

Case FCR-DU: participate only in FCR-D up,

  * •

Case FCR-DD: participate only in FCR-D down,

  * •

Case multi: multi-market participation is allowed at each hour.

The five cases are run twice, once when degradation cost is included in the
objective function, and once without degradation cost. For the latter, the
degradation cost is post-calculated and included in calculating the profit.
The optimization problem is run for one day at a time and for a full year
using a time-step of one minute.

###  IV-A Comparison metrics

The results are assessed from two perspectives: monetary, and operation
strategy. The monetary perspective compares the cases based on profit and
degradation costs. The operation strategy compares the distribution of SOC,
baseline power, and bid size. In addition, it compares the number of hours for
participating in different combinations of markets.

###  IV-B Input data

The input data consists of the input parameters for the battery whose
degradation model was formulated and the market-related data including
frequency data and prices for the DA market, FCR-N, FCR-D up, FCR-D down, and
imbalance. Battery input parameters are provided in Table II, the battery data
are taken from [19], prices are taken from [21, 22], and frequency data is
taken from [23]. Market zone SE3 and the year 2022 are used for prices and
frequency data.

TABLE II: Battery input parameters Parameter |  Value |  Parameter |  Value
---|---|---|---
SoC range |  10%-90% |  Battery capacity |  1 MWh
Dis/Charge efficiency |  93% |  Min, max dis/charge power |  0.0MW, 1.0MW
Battery OM cost |  $2\%\cdot 137$ k€/MWh |  Battery replacement cost |  137 k€/MWh
EOL |  80% |  Interest rate |  5%
Battery calendar lifetime |  10 years |  Ratio of salvage cost |  0.5
Temperature |  20 ∘ C |  |
Cycle aging parameters[20]
At 20${}^{\circ}\mathrm{C}$:
$q_{1}{\mathcal{K}_{t}}^{2}+q_{2}\mathcal{K}_{t}+q_{3}=0.0008$, $q4=0.3903$
Calendar aging parameters[20, 24]
$E_{a}=24.5\mathrm{kJmol^{-1}},R=8.314J\mathrm{mol^{-1}K^{-1}}$, $a_{1}=-1.1$,
$a_{2}=89.7$, $a_{3}=1224.6$, $b_{1}=10.3$,$b_{2}=-1083.6$, $b_{3}=31447$,
$c_{1}=2.6$, $c_{2}=-409.5$, $c_{3}=22035$

##  V Results and Discussion

Figure 5: The battery utilization in different cases visualized by the
histogram of battery SoE ($\mathcal{S}_{t}$) and power ($p_{t}$) as a
percentage of all time step in a year

Monetary results are presented in Table III. Profit for cases without
considering degradation in the objective function is obtained by post-
calculating the aging cost and subtracting it from the objective function
values. The results show that the multi-market case has a significantly larger
profit compared to other market participation cases while having the second
smallest degradation cost. Another observation is that considering degradation
in the optimization does not have a considerable impact ($\mathrm{\leq 1\%}$)
on the profit except in case w/o FCR. However, the total annual aging cost has
been reduced by 5%-29% when battery degradation is considered in the
optimization. This highlights that considering degradation in the objective
function can provide operation strategies that not only result in a similar
profit but also a longer lifetime for the battery. The utilization of the
battery is visualized in Figure 5 for clarifying the reduction of battery
degradation. The SoE and power data points show a shift towards smaller
absolute values when degradation is considered. This explains the reduction of
the aging costs when degradation is considered in the optimization.

TABLE III: Annual monetary results in k€ | W/o deg. in objective | With deg. in objective |  $\Delta$ Cost tot. age.
---|---|---|---
Case |  Profit |  Cost cal. age. |  Cost cyc. age. |  Cost tot. age. |  Profit |  Cost cal. age. |  Cost cyc. age. |  Cost tot. age.
w/o FCR |  30 |  5.7 |  3.9 |  9.6 |  32 |  4.2 |  3.4 |  7.6 |  -21%
FCR-N |  213 |  3.9 |  6.2 |  10.2 |  214 |  3.8 |  5.9 |  9.6 |  -5%
FCR-DU |  560 |  6.3 |  3.1 |  9.5 |  561 |  4.4 |  2.6 |  6.9 |  -27%
FCR-DD |  303 |  3.5 |  2.7 |  6.2 |  303 |  2.6 |  2.3 |  4.9 |  -21%
Multi |  706 |  4.5 |  3.0 |  7.5 |  708 |  3.1 |  2.2 |  5.4 |  -29%

The total capacity loss of the battery over the year,
$\sum_{d}\sum_{t}\mathcal{U}_{t,d}(\%)$, is presented in Figure 6. The largest
degradation is expected in the FCR-N cases while the lowest degradation is
expected for the FCR-D down cases. This is due to the large energy throughput
and cycling for FCR-N cases while the FCR-D down cases require both lower
cycling and SoC levels. In addition, FCR-N cases have the largest cycle aging
costs while FCR-D up cases show the largest calendar aging. It is also worth
noting that the multi-market participation cases have shown a lower
degradation compared to the w/o FCR because the battery is utilized at lower
SoC and power setpoints.

Figure 6: Distribution of annual battery capacity loss for the different
market participation cases

To better understand the difference in the operation strategies, the
distribution of SoE in the beginning of each hour ($\mathcal{S}_{h}$),
baseline power ($p^{bl}_{h}$), and bid size ($p^{\Theta}$) are presented in
Figure 7. The dashed lines show the first, second, and third quartiles. SoE at
the beginning of each hour is considered as a part of the strategy because it
is one of the main variables for satisfying the endurance requirements. For
multi-market cases, the sum of bids in all the FCR markets is presented in
Figure 7.

Figure 7: Operation strategy including the distribution of SoE at the
beginning of each hour ($\mathcal{S}_{h}$), baseline power($p^{bl}_{h}$), and
bid size ($p^{\Theta}_{h}$) for the simulated cases

Several general observations can be made regarding the most dominant values
for the strategy variables presented in Figure 7. SoE at the beginning of each
hour ($\mathcal{S}_{h}$) seems less discrete in the FCR-N case compared to the
other cases. The dominant baseline power for all the cases is zero MW while
the FCR-N case shows again a less discrete distribution. These two
observations can highlight the potential complexities of real-life planning of
an optimal FCR-N strategy compared to the other market participation cases
which most likely boils down to the large energy throughput and a more regular
activation of FCR-N services.

The dominant bid sizes for cases FCR-D up and FCR-D down is 1 MW (Figure 7).
The larger bids were possible at hours with a non-zero reference power. The
bids in the FCR-N case are limited to 0.4 MW due to considering technical
requirements and the battery capacity. FCR-N is a symmetrical service with a
required 1h endurance. Therefore, for example, if the SoE is at 0.5MWh, only
0.4MW can be provided for 1h in each direction considering the allowed SoC
range of 10%-90%. For case Multi, the sum of the bids in all FCR markets are
presented showing a dominant sum at 1.6 MW. This is the dominant optimal
bidding strategy which is a simultaneous bid in FCR-D up and down markets
(Table IV). The sum is limited to 1.6 MW because of the technical requirement
of power indicating a 20% power availability in the opposite direction. The
bid size results are yet another observation regarding the impact of
considering technical requirements and their importance in obtaining realistic
bidding and operation strategies.

TABLE IV: Choice of market as the number of hours in a year Case |  None |  N |  DU |  DD |  N+DU |  N+DD |  DU+DD |  All
---|---|---|---|---|---|---|---|---
Without deg. in objective function
w/o FCR |  8760 |  - |  - |  - |  - |  - |  - |  -
FCR-N |  892 |  7868 |  - |  - |  - |  - |  - |  -
FCR-DU |  73 |  - |  8687 |  - |  - |  - |  - |  -
FCR-DD |  279 |  - |  - |  8481 |  - |  - |  - |  -
Multi |  0 |  3 |  1741 |  145 |  164 |  76 |  6110 |  521
| With deg. in objective function
w/o FCR |  8760 |  - |  - |  - |  - |  - |  - |  -
FCR-N |  715 |  8045 |  - |  - |  - |  - |  - |  -
FCR-DU |  70 |  - |  8690 |  - |  - |  - |  - |  -
FCR-DD |  202 |  - |  - |  8558 |  - |  - |  - |  -
Multi |  0 |  2 |  1540 |  147 |  93 |  75 |  6382 |  521

The mix of chosen markets in each case is presented in Table IV. The inclusion
of degradation in the objective function has led to a larger number of market
participation hours. This is because the BESS can reduce calendar aging by
regulating SoC. The change in the strategy comprises a larger number of
participation hours but with a reduction in the number of hours that have the
largest bidding power. This effect can be seen more clearly for the largest
bids in cases FCR-D down and FCR-D up presented in Figure 7. The dominant bid
combination for case Multi is simultaneous FCR-DU and FCR-DD participation.
This dominant bid combination besides the SoE levels in Figure 5 can clarify
why the total battery degradation in Multi case is in between FCR-D up and
FCR-D down cases (Figure 6).

It is worth noting that the model should not be seen as a ready-for-market
bidding algorithm. The presented profit values are the theoretical maximum
profit given the perspective of an oracle. In real-life, these profits might
not be achieved to their full extent due to the lack of information and
uncertainties in the input parameters. However, the model and its results can
be interpreted as a benchmark indicating the maximum potential profit, and the
potential impact of considering battery degradation in the optimization.
Hence, the model can be used as an oracle model to evaluate and compare real-
life bidding models that takes into account uncertainties in input parameters.

##  VI Conclusions

In this study, a novel MILP formulation was presented for Swedish multi-FCR
participation including a detailed cycle and calendar aging besides the most
recent technical requirements. The model was used for a full-year case study
of the year 2022 to analyze the impact of degradation and technical
requirements.

The results have shown that although considering degradation in the
optimization does not have a significant impact on the profit, it can decrease
the aging by 5%-29% leading to a more sustainable utilization of the battery.
The results have also clearly demonstrated the impact and importance of
fulfilling the technical requirements of the Swedish market. In addition, the
highest profit could have been achieved in the multi-market participation that
was dominated by simultaneous bidding in FCR-D up and down markets. The
largest total degradation was for the dedicated FCR-N participation case while
the lowest degradation was shown to be for the dedicated FCR-D down
participation case.

The proposed model and the results can be utilized by flexibility asset owners
to obtain more sustainable and yet profitable operation strategies. In
addition, the formulation can be seen as an oracle model providing a maximum
potential profit. This can be utilized as a benchmark for evaluating bidding
models considering uncertainties.

## References

  * [1] N. Modig, R. Eriksson, P. Ruokolainen, J. N. Ødegård, S. Weizenegger, and T. D. Fechtenburg, “Overview of frequency control in the nordic power system,” _Nordic Analysis Group_ , 2022.
  * [2] L. Meng, J. Zafar, S. K. Khadem, A. Collinson, K. C. Murchie, F. Coffele, and G. M. Burt, “Fast frequency response from energy storage systems—a review of grid standards, projects and technical issues,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 2, pp. 1566–1581, 2020.
  * [3] Y. Hu, M. Armada, and M. J. Sánchez, “Potential utilization of battery energy storage systems (bess) in the major european electricity markets,” _Applied Energy_ , vol. 322, p. 119512, 2022.
  * [4] T. Zhao, A. Parisio, and J. V. Milanović, “Distributed control of battery energy storage systems for improved frequency regulation,” _IEEE Transactions on Power Systems_ , vol. 35, no. 5, pp. 3729–3738, 2020.
  * [5] Q. Ma, W. Wei, L. Wu, and S. Mei, “Life-aware operation of battery energy storage in frequency regulation,” _IEEE Transactions on Sustainable Energy_ , vol. 14, no. 3, pp. 1725–1736, 2023.
  * [6] E. N. Sofia Guzman, M. Arriaga, C. A. Cañizares, J. W. Simpson-Porco, D. Sohm, and K. Bhattacharya, “Regulation signal design and fast frequency control with energy storage systems,” _IEEE Transactions on Power Systems_ , vol. 37, no. 1, pp. 224–236, 2022.
  * [7] A. Oshnoei, M. Kheradmandi, and S. M. Muyeen, “Robust control scheme for distributed battery energy storage systems in load frequency control,” _IEEE Transactions on Power Systems_ , vol. 35, no. 6, pp. 4781–4791, 2020.
  * [8] P. Astero and C. Evens, “Optimum operation of battery storage system in frequency containment reserves markets,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 6, pp. 4906–4915, 2020.
  * [9] I. M. Casla, A. Khodadadi, and L. Söder, “Optimal day ahead planning and bidding strategy of battery storage unit participating in nordic frequency markets,” _IEEE Access_ , vol. 10, pp. 76 870–76 883, 2022.
  * [10] S. Gao, H. Li, J. Jurasz, and R. Dai, “Optimal charging of electric vehicle aggregations participating in energy and ancillary service markets,” _IEEE Journal of Emerging and Selected Topics in Industrial Electronics_ , vol. 3, no. 2, pp. 270–278, 2022.
  * [11] I. Pavić, H. Pandžić, and T. Capuder, “Electric vehicle aggregator as an automatic reserves provider under uncertain balancing energy procurement,” _IEEE Transactions on Power Systems_ , vol. 38, no. 1, pp. 396–410, 2023.
  * [12] R. K. Subroto, D. Gebbran, A. B. Moreno, and T. Dragičević, “Bess optimal sizing and scheduling for energy arbitrage and frequency containment reserve via dual-loop optimization,” in _2022 IEEE Transportation Electrification Conference & Expo (ITEC)_. IEEE, 2022, pp. 941–946.
  * [13] A. Thingvad, C. Ziras, G. Le Ray, J. Engelhardt, R. R. Mosbæk, and M. Marinelli, “Economic value of multi-market bidding in nordic frequency markets,” in _2022 International Conference on Renewable Energies and Smart Technologies (REST)_ , vol. 1. IEEE, 2022, pp. 1–5.
  * [14] K. Jacqué, L. Koltermann, J. Figgener, S. Zurmühlen, and D. U. Sauer, “The influence of frequency containment reserve on the cycles of a hybrid stationary large-scale storage system,” _Journal of Energy Storage_ , vol. 52, p. 105040, 2022.
  * [15] A. Krupp, R. Beckmann, P. Draheim, E. Meschede, E. Ferg, F. Schuldt, and C. Agert, “Operating strategy optimization considering battery aging for a sector coupling system providing frequency containment reserve,” _Journal of Energy Storage_ , vol. 68, p. 107787, 2023.
  * [16] L. Koltermann, M. C. Cortés, J. Figgener, S. Zurmühlen, and D. U. Sauer, “Power curves of megawatt-scale battery storage technologies for frequency regulation and energy trading,” _Applied Energy_ , vol. 347, p. 121428, 2023.
  * [17] ENTSOe, “Technical requirements for frequency containment reserve provision in the nordic synchronous area,” ENTSOe, Tech. Rep., 2023.
  * [18] Information on different ancillary services. [Online]. Available: https://www.svk.se/en/stakeholders-portal/electricity-market/provision-of-ancillary-services/information-on-different-ancillary-services/
  * [19] R. Khezri, D. Steen, E. Wikner, and L. A. Tuan, “Optimal v2g scheduling of an ev with calendar and cycle aging of battery: An milp approach,” _IEEE Transactions on Transportation Electrification_ , pp. 1–1, 2024.
  * [20] J. Wang, J. Purewal, P. Liu, J. Hicks-Garner, S. Soukazian, E. Sherman, A. Sorenson, L. Vu, H. Tataria, and M. W. Verbrugge, “Degradation of lithium ion batteries employing graphite negatives and nickel–cobalt–manganese oxide+ spinel manganese oxide positives: Part 1, aging mechanisms and life estimation,” _Journal of Power Sources_ , vol. 269, pp. 937–948, 2014.
  * [21] ENTSO-E Transparency Platform. [Online]. Available: https://transparency.entsoe.eu/
  * [22] eSett Open Data. [Online]. Available: https://opendata.esett.com/
  * [23] Fingrid Open Data. [Online]. Available: https://data.fingrid.fi/en
  * [24] L. Calearo, A. Thingvad, and M. Marinelli, “Modeling of battery electric vehicles for degradation studies,” in _2019 54th International Universities Power Engineering Conference (UPEC)_. IEEE, 2019, pp. 1–6.