# Using polarons for sub-$\mathrm{nK}$ quantum non-demolition thermometry in a
Bose-Einstein condensate

Mohammad Mehboudi  ICFO—Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Departament de Física, Universitat Autònoma de Barcelona, E-08193 Bellaterra,
Spain  Aniello Lampo  ICFO—Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Christos Charalambous  ICFO—Institut de Ciències Fotòniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Luis A. Correa  School of Mathematical Sciences and Centre for the Mathematics
and Theoretical Physics of Quantum Non-Equilibrium Systems, The University of
Nottingham, University Park, Nottingham NG7 2RD, United Kingdom  Kavli
Institute for Theoretical Physics University of California, Santa Barbara, CA
93106  Miguel Ángel García-March  ICFO—Institut de Ciències Fotòniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  Maciej Lewenstein  ICFO—Institut de Ciències Fotòniques,
The Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  ICREA, Lluís Companys 23, E-08010 Barcelona, Spain

###### Abstract

We introduce a novel minimally-disturbing method for sub-$\mathrm{nK}$
thermometry in a Bose-Einstein condensate (BEC). Our technique is based on the
Bose-polaron model; namely, an impurity embedded in the BEC acts as the
thermometer. We propose to detect temperature fluctuations from measurements
of the position and momentum of the impurity. Crucially, these cause minimal
back-action on the BEC and hence, realize a non-demolition temperature
measurement. Following the paradigm of the emerging field of quantum
thermometry, we combine tools from quantum parameter estimation and the theory
of open quantum systems to solve the problem in full generality. We thus avoid
any simplification, such as demanding thermalization of the impurity atoms, or
imposing weak dissipative interactions with the BEC. Our method is illustrated
with realistic experimental parameters common in many labs, thus showing that
it can compete with state-of-the-art destructive techniques, even when the
estimates are built from the outcomes of accessible (sub-optimal) quadrature
measurements.

Introduction— The ongoing efforts in the development of quantum technologies
is strongly fuelled by their many anticipated practical applications Celi  _et
al._ (2017). In the process, we are already benefiting from striking
experimental advances and much deeper theoretical insights. In particular,
ultracold atomic gases are a key platform for quantum technologies due to
their potential for quantum simulation Bloch  _et al._ (2008); Lewenstein  _et
al._ (2012). Nonetheless, operating a quantum simulator requires very precise
tuning of the parameters of the experiment, so as to ensure that the simulated
system behaves as intended. In particular, a precise temperature control is
essential, for instance, for the reconstruction of the equation of state of
the system Bloch  _et al._ (2012).

In current experimental setups, the main thermometric techniques are based on
time-of-flight measurements either directly on the BEC Leanhardt  _et al._
(2003); Gati  _et al._ (2006a, b), or on impurities embedded in it Olf  _et
al._ (2015); Spiegelhalder  _et al._ (2009). In the former case, temperatures
of few $\mathrm{nK}$, or even sub-$\mathrm{nK}$ might be estimated
efficiently, although at the price of destroying the BEC. On the contrary, the
latter protocols are less destructive, albeit efficient at relatively “large”
temperatures of $\sim 100$ $\mathrm{nK}$. Interestingly, recent proposals have
discussed minimally disturbing interferometric setups in which the temperature
is mapped onto a relative phase on a probe Stace (2010); Martín-Martínez  _et
al._ (2013); Sabín  _et al._ (2014), however, the underlying models are very
simple.

An effective non-demolition thermometric technique in the sub-$\mathrm{nK}$
regime is thus still missing. Any such strategy should be build upon a
comprehensive theoretical description and be capable of informing the choice
of the most sensitive temperature-dependent quantities to be measured. Here,
we propose what is, to the best of our knowledge, the first experimentally
feasible quantum non-demolition technique to measure the temperature of a BEC
in the sub-$\mathrm{nK}$ domain. It is based on the Bose polaron problem,
i.e., interrogation of an impurity that is embedded in the condensate, while
causing minimal disturbance to the cold atomic gas. The impurity problem has
been intensively studied in the context of polaron physics in strongly-
interacting Fermi Schirotzek  _et al._ (2009); Kohstall  _et al._ (2012);
Koschorreck  _et al._ (2012); Massignan  _et al._ (2014); Schmidt  _et al._
(2016); Lan and Lobo (2014); Levinsen and Parish (2014); Schmidt  _et al._
(2012); Płodzień  _et al._ (2018) or Bose gases Rath and Schmidt (2013); Côté
_et al._ (2002); Massignan  _et al._ (2005); Cucchietti and Timmermans (2006);
Lampo  _et al._ (2017); Pastukhov (2017); Yoshida  _et al._ (2018); Guenther
_et al._ (2018); Lampo  _et al._ (2018), as well as in solid state physics
Landau and Pekar (1948); Devreese and Alexandrov (2009); Alexandrov and
Devreese (2009), and mathematical physics Lieb and Yamazaki (1958); Lieb and
Thomas (1997); Frank  _et al._ (2010); Anapolitanos and Landon (2013); Lim
_et al._ (2018). We specifically avoid any unjustified simplifications—such as
complete thermalization of the impurities at the BEC temperature—and
investigate the problem in its full generality. The usefulness of our proposed
technique is finally illustrated with typical experimental parameters.

In our analysis, we benefit from the toolbox of the emergent field of quantum
thermometry Mehboudi  _et al._ , which combines quantum estimation theory and
the theory of open quantum systems. This will allow us to compare the ultimate
precision bounds on temperature estimation with the thermal sensitivity of
concrete experimentally feasible measurements.

The model— Let us consider an impurity of mass $m_{\rm I}$ (acting as the
temperature probe) embedded in a BEC of atoms of mass $m_{\rm B}$, chemical
potential $\mu$, and interatomic coupling strength $g_{\rm B}$. The condensate
(which makes up the sample), is confined in a one-dimensional harmonic well of
frequency $\omega_{\rm B}$, leading to a parabolic Thomas-Fermi potential with
radius $R=\sqrt{2\mu/m_{\rm B}\omega^{2}_{\rm B}}$. In turn, the impurity is
trapped in a harmonic potential of frequency $\Omega$. Finally, the
interspecies coupling (i.e., the probe-sample interaction) is denoted by
$g_{\rm IB}$. Our aim is to estimate the temperature $T$ of the BEC as
precisely as possible, while diminishing the ensuing disturbance.

We note that the global probe-sample Hamiltonian can be thought-of as a
quantum Brownian motion model consisting of the following contributions:

$\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{\rm I}+\hat{H}_{\rm
B}+\hat{H}_{\rm int}$ $\displaystyle=\frac{{\hat{p}}^{2}}{2m_{\rm
I}}+\frac{m_{\rm
I}\Omega^{2}}{2}\hat{x}^{2}+\sum\nolimits_{k}E_{k}\hat{b}^{\dagger}_{k}\hat{b}_{k}+\sum\nolimits_{k}\hbar
g_{k}\,\hat{x}\left(\hat{b}_{k}+\hat{b}^{\dagger}_{k}\right).$ (1)

Here, $\hat{H}_{\rm I}$ stands for the free Hamiltonian of the impurity, while
the term $\hat{H}_{\rm B}$ represents the BEC and encompasses all the
interacting modes of the atomic gas. Finally, $\hat{H}_{\rm int}$ is the
interaction between the impurity and the gas. In the second line of Eq. (1),
however, we express the BEC degrees of freedom in terms of the operators
$\hat{b}_{k}$ and $\hat{b}_{k}^{\dagger}$, that stand for the annihilation and
creation operators of a Bogoliubov mode with energy $E_{k}$. This new
representation is nothing but the result of a standard Bogoliubov
transformation, which diagonalizes $\hat{H}_{\rm B}$ and maps it into a non-
interacting form. In this picture, the last term in the second line of (1)
accounts for the interactions between the impurity and the Bogoliubov modes.
The interaction strength between the $k$th Bogoliubov mode and the impurity is
given by $g_{k}$, (see Lampo  _et al._ (2018) for details). Note that such
interactions exhibit, in general, a non-linear dependence on the position of
the impurity. The linear form presented in Eq. (1) is valid only near the
center of the confining potential of the BEC, i.e., when $x/R\ll 1$. Of
course, this leads to constraints on the values of the system parameters that
have been discussed in Lampo  _et al._ (2018). We underline that the values of
the physical quantities considered in the following fulfill the conditions
associated to the linear approximation in the interaction Hamiltonian.

In general, the spectrum of a Hamiltonian like $\hat{H}$ may not be bounded
from below Caldeira and Leggett (1983a, b); Weiss (2008). To rule out this
eventuality, it is common practice to shift the frequency of the Brownian
particle as
$\Omega^{2}\mapsto\tilde{\Omega}^{2}\coloneqq\Omega^{2}+2\sum_{k}g_{k}^{2}/\omega_{k}^{2}$
to compensate for the distortion caused by the coupling to the bath Correa
_et al._ (2017); Hovhannisyan and Correa (2017). Here, however, we will avoid
adding any terms “by hand”, since Eq. (1) has been consistently derived from
first principles Lampo  _et al._ (2018). Instead, we will limit ourselves to
choose parameters which fulfill the conditions described in Lampo  _et al._
(2018) that guarantee the positivity of the spectra.

Treating the impurity as a Brownian particle coupled to a bath of Bogoliubov
modes, allows us to exploit well-established techniques from the theory of
open quantum systems. Specifically, the motion of the probe is described
exactly by the second-order differential equation Weiss (2008)

$\frac{d^{2}\hat{x}(t)}{dt^{2}}+\Omega^{2}\hat{x}(t)+\frac{d}{dt}\int^{t}_{0}ds\,\Gamma(t-s)\,\hat{x}(s)=\frac{\hat{B}(t)}{m_{\rm
I}}.$ (2)

This is the quantum counterpart of the Langevin equation, introduced in 1909
for (classical) Brownian motion Langevin (1908); Mazo (2002). The term
$\hat{B}(t)$ on the right-hand side of Eq. (2) reads

$\hat{B}(t)\coloneqq\sum_{k}\hbar
g_{k}\left(\hat{b}^{\dagger}_{k}e^{-i\omega_{k}t}+\hat{b}_{k}e^{i\omega_{k}t}\right),$
(3)

and plays the role of a stochastic driving force. Memory effects enter in the
dynamics through the damping kernel

$\Gamma(\tau)\coloneqq\frac{1}{m_{\rm
I}}\int^{\infty}_{0}d\omega\,\frac{J(\omega)}{\omega}\cos(\omega\tau),$ (4)

where $J(\omega)\coloneqq\sum_{k\neq 0}\hbar
g_{k}^{2}\delta(\omega-\omega_{k})$ is the spectral density. For our model,
this is given by

$\displaystyle J(\omega)=m_{\rm I}\gamma\,\frac{\omega^{4}}{\omega^{3}_{\rm
B}}\,\Theta\left(\omega-\omega_{\rm B}\right),$ (5)

where $\gamma\coloneqq(2g_{\rm B}/m_{\rm I}\omega_{\rm
B}R^{3})\times(\eta\,\mu/\hbar\omega_{\rm B})^{2}$ and $\eta\coloneqq g_{\rm
IB}/g_{\rm B}$ Lampo  _et al._ (2018). The Heaviside function $\Theta(\cdot)$
introduces an ultraviolet cutoff, which regularizes the diverging high-
frequency behavior. Importantly, the long-time dynamics of the impurity—our
main focus—does not depend on the details of the cutoff Lampo  _et al._ (2017)
but rather, on the low frequency behaviour of $J(\omega)$. This is dictated by
the exponent of $\omega$ in the pre-factor of Eq. (5)—i.e., the “Ohmicity”
parameter Correa  _et al._ (2017).

Eq. (2) can be solved via Laplace transform (that is,
$\mathcal{L}_{z}[f(t)]\coloneqq\int_{0}^{\infty}dt\,e^{-tz}f(t)$). In
particular, the steady-state variances in position and momentum are given by

$\displaystyle\langle\hat{x}^{2}\rangle=\frac{\hbar}{2\pi}\int^{+\omega_{\rm
B}}_{-\omega_{\rm B}}d\omega\,\coth{(\hbar\omega/2k_{\rm
B}T)}\,\tilde{\chi}^{\prime\prime}(\omega),$ (6a)
$\displaystyle\langle\hat{p}^{2}\rangle=\frac{\hbar m^{2}_{\rm
I}}{2\pi}\int^{+\omega_{\rm B}}_{-\omega_{\rm
B}}d\omega\,\omega^{2}\coth{(\hbar\omega/2k_{\rm
B}T)}\,\tilde{\chi}^{\prime\prime}(\omega),$ (6b)

where $\langle\cdots\rangle$ denotes steady-state averaging and
$\tilde{\chi}^{\prime\prime}(\omega)$ stands for the response function that
reads as

$\displaystyle\tilde{\chi}^{\prime\prime}(\omega)=\frac{1}{m_{\rm
I}}\frac{\omega\,\zeta(\omega)}{\left[\omega\,\zeta(\omega)\right]^{2}+\left[\Omega^{2}-\omega^{2}+\omega\,\theta(\omega)\right]^{2}}.$
(7)

$\zeta\left(\omega\right)$ and $\theta\left(\omega\right)$ in Eq. (7) are,
respectively, the real and imaginary parts of
$\mathcal{L}_{\tilde{z}}\left[\Gamma(t)\right]$ evaluated at
$\tilde{z}=-i\omega+0^{+}$. Eqs. (4)–(7) thus allow to determine the steady-
state covariances of the impurity as a function of the system parameters—in
particular, the temperature of the BEC and the dissipation strength $\gamma$
Correa  _et al._ (2017). Note that, since the Hamiltonian (1) is linear in the
quadratures, Eqs. (6) fully characterize the steady state of the impurity
(together with $\langle\hat{x}\rangle=\langle\hat{p}\rangle=0$).

Thermometric performance— The inherent errors from quantum measurements give
rise to statistical uncertainty on the temperature estimate. Quantum
estimation theory allows us to place fundamental limits on the “error bars” of
the final temperature reading, and even to rank the various temperature-
dependent properties of the probe according to their thermal sensitivity. For
instance, let us build our temperature estimate from a large set of $\nu$
outcomes of independent measurements of some impurity observable $\hat{O}$
111Note that we are not only assuming to work with a large dataset, but also
that the estimator which maps measurement outcomes to temperature estimates is
unbiased Braunstein and Caves (1994).. We stress that these are either
measurements performed on independent impurity atoms, or measurements on the
same probe, but paced so that the BEC-impurity composite has time to reset to
its stationary state every time. By mere propagation of errors, the
uncertainty of the temperature inferred from such dataset would read Tóth and
Apellaniz (2014); Braunstein and Caves (1994)

$\displaystyle\delta
T({\hat{O}})\coloneqq\frac{\Delta{\hat{O}}}{\sqrt{\nu~{}\chi_{T}^{2}({\hat{O}})}},$
(8)

where
$\Delta^{2}{\hat{O}}\coloneqq\langle{\hat{O}}^{2}\rangle-\langle\hat{O}\rangle^{2}$
stands for the variance of $\hat{O}$ calculated on the stationary marginal of
the impurity $\hat{\varrho}_{\rm I}(T)$, and
$\chi_{T}({\hat{O}})\coloneqq\partial_{\xi}\,{\rm tr}\,[\hat{\varrho}_{\rm
I}(\xi)~{}\hat{O}]_{\xi=T}$ represents its (static) temperature
susceptibility.

In order to assess the performance of $\hat{O}$, it is essential to know which
is the minimum possible uncertainty (i.e., $(\delta
T)_{\min}\coloneqq\inf_{\hat{O}}\delta T(\hat{O})$). To this end, we introduce
the symmetric logarithmic derivative (SLD) $\hat{\Lambda}_{T}$, implicitly
defined as

$\displaystyle{\hat{\Lambda}}_{T}\,\hat{\varrho}_{\rm I}(T)+\hat{\varrho}_{\rm
I}(T)\,{\hat{\Lambda}}_{T}\equiv 2\left.\partial_{\xi}\,{\hat{\varrho}_{\rm
I}(\xi)}\,\right|_{\xi=T}.$ (9)

Coming back to the definition of $\chi_{T}(\hat{O})$, we notice that
$\chi_{T}({\hat{O}})=\frac{1}{2}\langle\hat{O}\,{\hat{\Lambda}}_{T}+{\hat{\Lambda}}_{T}\,{\hat{O}}\rangle-\langle\hat{O}\rangle\langle\hat{\Lambda}_{T}\rangle$,
while $\chi_{T}(\hat{\Lambda}_{T})=\Delta^{2}\hat{\Lambda}_{T}$. Making use of
the fact that $\Delta\hat{O}\,\Delta\hat{\Lambda}_{T}\geq\chi_{T}(\hat{O})$
allows to turn Eq. (8) into

$\displaystyle\delta
T({\hat{O}})\geq\frac{1}{\sqrt{\nu~{}\Delta^{2}{\hat{\Lambda}}_{T}}}\coloneqq\frac{1}{\sqrt{\nu~{}{\cal
F}(T)}},$ (10)

where we have introduced the quantum Fisher information (QFI) ${\cal
F}(T)\coloneqq\Delta^{2}{\hat{\Lambda}}_{T}$. Eq. (10) is nothing but the
quantum Cramér-Rao bound Braunstein and Caves (1994), and sets the ultimate
lower limit on the statistical error. Furthermore, by simply replacing
${\hat{O}}$ by ${\hat{\Lambda}_{T}}$ in (8), we can see that this bound is
saturated by performing complete projective measurements onto the eigenbasis
of the SLD.

Results— Owing to the simplified Hamiltonian in Eq. (1) we can write the SLD
and the QFI for temperature estimation solely in terms of the variances in
Eqs. (6) Monras (2013); Correa  _et al._ (2017), i.e.,

$\displaystyle{\hat{\Lambda}}_{T}$
$\displaystyle=C_{x}\left({\hat{x}}^{2}-\langle{\hat{x}}^{2}\rangle\right)+C_{p}\left({\hat{p}}^{2}-\langle{\hat{p}}^{2}\rangle\right),$
(11a) $\displaystyle{\cal F}(T)$
$\displaystyle=2\,C_{x}^{2}\,\langle{\hat{x}}^{2}\rangle^{2}+2\,C_{p}^{2}\,\langle{\hat{p}}^{2}\rangle^{2}-\hbar^{2}\,C_{x}\,C_{p},$
(11b)

where the coefficient $C_{x}$ is given by

$\displaystyle C_{x}$
$\displaystyle=\frac{4\,\langle{\hat{p}}^{2}\rangle^{2}\,\chi_{T}({\hat{x}}^{2})+\hbar^{2}\,\chi_{T}({\hat{p}}^{2})}{8\,\langle{\hat{x}}^{2}\rangle^{2}\,\langle{\hat{p}}^{2}\rangle^{2}-\hbar^{4}/2},$
(12)

and $C_{p}$ can be obtained by simply exchanging $\hat{x}$ and $\hat{p}$. That
is, by repeatedly measuring the observable (11a) on the impurity, the
temperature of the BEC can be estimated with the minimum possible uncertainty.
We are now in the position to plug in realistic numbers into the exact steady-
state marginal for the probe and explore the thermal sensitivity of our non-
demolition thermometric protocol at ultra-low temperatures.

Figure 1: (color online) (black) Optimal relative error $(\delta T)_{\rm
min}/T$ ($\nu=1$) as a function of the temperature of the BEC in a logarithmic
scale. Specifically, we work with impurities of Yb in a sea of ultra-cold K.
The temperature range for the BEC is $200\,$\mathrm{pK}$\leq T\leq
2\,$\mathrm{nK}$$. The trapping frequency of the gas (with $N=5000$ atoms) was
set to $\omega_{\rm B}=2\pi\times 100\,$\mathrm{Hz}$$, while
$\Omega=2\pi\times 10\,$\mathrm{Hz}$$ ($g_{rmB}=3\times
10^{-39}\,$\mathrm{J}\text{\,}\mathrm{m}$$). Different probe-sample coupling
ratios $\eta=g_{\rm IB}/g_{\rm B}$ were considered, namely (solid) $\eta=1$,
(dashed) $\eta=3$, and (dotted) $\eta=6$. For comparison, we also depicted the
relative error of a fully thermalized impurity (i.e., $\eta\rightarrow 0$)
(dot-dashed red). Note that, for $\eta=1$, the relative error can be kept
below $\sim 14$% from only $\nu\sim 100$ measurements. This is quantitatively
close to state-of-the-art destructive experimental techniques. See text for
discussion.
---
Figure 2: (color online) (dashed blue) Relative error for the position
quadrature $\delta T(\hat{x}^{2})/T$ and (dotted red) the momentum quadrature
$\delta T(\hat{x}^{2})/T$ as a function of $T$ for (left panel) $\eta=1$ and
(right panel) $\eta=6$. All parameters are the same as in Fig. 1. The minimum
relative error (solid black) is superimposed for reference. Even though both
measurement schemes are sub-optimal, they still might allow to draw estimates
with relative errors as low as $18\%$ for $\nu\sim 400$. Note that for $T\geq
0.5$ nK and by using the same data size $\nu=400$, one can achieve a relative
error below $10\%$.

As an illustrative example, we will work with a BEC of K atoms containing Yb
impurities. The qualitative picture would remain essentially unaltered
regardless of the atomic species considered. In Fig. 1, we plot the optimal
relative error $\sqrt{\nu}\,(\delta
T)_{\min}/T=\big{(}T\sqrt{\mathcal{F}_{T}}\big{)}^{-1}$ for various probe-
sample coupling strengths and temperatures ranging from $200\,$\mathrm{pK}$$
to $2\,$\mathrm{nK}$$. Specifically, keeping the interatomic and interspecies
couplings comparable (i.e., $\eta=1$) would allow to achieve a relative error
below $14$% from as few as $100$ measurement outcomes. That is, polaron
thermometry outperforms the interferometric technique proposed in Ref. Sabín
_et al._ (2014) by an order of magnitude. More importantly, unlike state-of-
the-art experimental methods (e.g., Leanhardt  _et al._ (2003); Gati  _et al._
(2006b, a); Sherson  _et al._ (2010); Olf  _et al._ (2015)), ours is non-
destructive.

We note, however, that the stronger the probe sample interaction, the worse
the estimation. Likewise, it can be clearly seen that, for strong dissipation,
the impurity deviates significantly from a thermal state at the temperature of
the sample. The first observation seems to be in striking contradiction with
the main results of Correa  _et al._ (2017), where a substantial dissipation-
driven enhancement was reported at low temperatures. Note however, that the
temperature range considered in Fig. 1 does not qualify as “low”, according to
the criteria of Ref. Correa  _et al._ (2017), namely $T\ll\hbar\omega_{\rm
B}/k_{B}$ (here, $T\sim\hbar\omega_{\rm B}/k_{B}$). When it comes to the
second observation, it is worth highlighting that the divergence between the
exact steady state of the impurity and a fully thermalised probe can be
sizeable in the $\mathrm{pK}$ range. This only comes to reinforce the idea
that simple dissipation models, such as a Gorini-Kossakowski-Lindblad-
Sudarshan master equation Lindblad (1976); Gorini  _et al._ (1976) are not
suitable for this type of analysis.

Recall that the above discussion assumes that the optimal measurement of Eq.
(11a) can be implemented. In practice, however, such a mixture of covariances
with temperature-dependent coefficients might be difficult to realize; the
bare quadratures $\langle\hat{x}^{2}\rangle$ or $\langle\hat{p}^{2}\rangle$
being easier to measure. The relative error of estimates based on these is
benchmarked against the ultimate lower bound in Figs. 2. Note that, at
$\eta=1$, $\langle\hat{x}^{2}\rangle$ and $\langle\hat{p}^{2}\rangle$ perform
similarly, while at stronger coupling, the position quadrature becomes a
significantly better temperature estimator. Also, under stronger dissipation,
$\langle\hat{x}^{2}\rangle$ gets closer to the optimal setting. Importantly,
our approach remains practically useful regardless of the strict sub-
optimality of $\hat{x}^{2}$—temperature estimates with $\delta
T(\hat{x}^{2})/T<18\%$ (or in the domain $T\geq 0.5$ nK, with $\delta
T(\hat{x}^{2})/T<10\%$) can still be constructed from relatively small
datasets of $\nu\sim 400$.

Conclusions— We have shown that impurities immersed in a BEC can be exploited
as temperature sensors. The key features of such thermometric scheme are that
(i) the temperature is estimated by monitoring the impurity atoms only—the BEC
itself does not need to be measured destructively, (ii) it can compete with
state-of-the-art thermometric techniques in the sub-$\mathrm{nK}$ range, and
(iii) the underlying analysis does not assume thermalization of the impurity
at the temperature of the BEC, but rather takes fully into account the strong
correlations built up between probe and sample.

In particular, we considered a cold atomic gas and an impurity both
harmonically confined in 1D at different trapping frequencies. Assuming that
the impurity remains localised around the minimum of the potential, allowed us
to “linearize” the model. We obtained the exact stationary state of the
impurity from the corresponding quantum Langevin equation and, using standard
tools from quantum estimation theory, we could eventually calculate the
minimum possible statistical uncertainty for a temperature measurement. In
particular, owing to our analysis being exact, we could verify that the usual
assumption of full thermalization for the impurity at the temperature of the
sample overestimates the performance of the scheme for typical parameters in
the $\mathrm{pK}$–$\mathrm{nK}$ range.

We showed that, with only $100$ measurements, the relative error can be kept
below $14\%$ for temperatures as low as $200\,$\mathrm{pK}$$. Importantly, we
could also show that feasible sub-optimal quadrature
measurements—specifically, $\hat{x}^{2}$—allow for similar performances with
limited resources (i.e., datasets of just few hundreds of independent
measurements). Interestingly, we found that increasing the probe-sample
coupling does not improve the sensitivity of the protocol in the temperature
range under study due to the comparatively low typical trapping frequencies
($60$–$70$ $\mathrm{Hz}$).

Even though we illustrate our results with Yb impurities in a cold gas of K
atoms, our approach is completely general and could be straightforwardly
applied without limitations to other atomic species and temperature ranges.
Similar results are also expected in the 2D and 3D cases. In particular, such
an extension is straightforward for homogeneous BECs, the same position
squeezing effects giving rise to the enhanced sensitivity of $\hat{x}^{2}$ are
known to occur Lampo  _et al._ (2017); the main difference would be a larger
Ohmicity in Eq. (5).

In order to bring these promising quantum non-demolition thermometers a step
closer to experimental demonstrations, it would be interesting to study how
the unavoidable non-linearities could affect our results. Exploring whether
the entanglement between two impurities embedded in the BEC—recently studied
in Charalambous  _et al._ (2018)—can be used to boost thermometric performance
also remains an open challenge.

###### Acknowledgements.

We acknowledge fruitful discussions with A. Sanpera and A. Acín. This work has
been funded by the Spanish Ministry MINECO (project QIBEQI FIS2016-80773-P,
National Plan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No.
SEV-2015-0522, FPI), European Social Fund, Fundació Privada Cellex,
Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341, SGR 1381 and
CERCA/Program), ERC AdG OSYRIS, EU FETPRO QUIC, the National Science Centre,
Poland-Symfonia Grant No. 2016/20/W/ST4/00314, ERC StG GQCOP No. 637352, the
COST Action MP1403: “Nanoscale quantum optics”, and US National Science
Foundation under Grant No. NSF PHY-1748958. Financial support through a
scholarship from Programa Màsters d’Excel·lència (Fundació Catalunya-La
Pedrera) is appreciated. LAC thanks the Kavli Institute for Theoretical
Physics for their warm hospitality during the program “Thermodynamics of
quantum systems: Measurement, engines, and control”.

## References

  * Celi  _et al._ (2017) A. Celi, A. Sanpera, V. Ahufinger, and M. Lewenstein, Phys. Scr. 92, 013003 (2017).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Lewenstein  _et al._ (2012) M. Lewenstein, A. Sanpera, and V. Ahufinger,  _Ultracold atoms in optical lattices_ , Vol. 143 (Oxford University Press Oxford, 2012).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbene, Nat. Phys. 8, 267 (2012).
  * Leanhardt  _et al._ (2003) A. Leanhardt, T. Pasquini, M. Saba, A. Schirotzek, Y. Shin, D. Kielpinski, D. Pritchard, and W. Ketterle, Science 301, 1513 (2003).
  * Gati  _et al._ (2006a) R. Gati, B. Hemmerling, J. Fölling, M. Albiez, and M. K. Oberthaler, Phys. Rev. Lett. 96, 130404 (2006a).
  * Gati  _et al._ (2006b) R. Gati, J. Esteve, B. Hemmerling, T. Ottenstein, J. Appmeier, A. Weller, and M. Oberthaler, New J. Phys. 8, 189 (2006b).
  * Olf  _et al._ (2015) R. Olf, F. Fang, G. E. Marti, A. MacRae, and D. M. Stamper-Kurn, Nat. Phys. 11, 720 (2015).
  * Spiegelhalder  _et al._ (2009) F. M. Spiegelhalder, A. Trenkwalder, D. Naik, G. Hendl, F. Schreck, and R. Grimm, Phys. Rev. Lett. 103, 223203 (2009).
  * Stace (2010) T. M. Stace, Phys. Rev. A 82, 011611 (2010).
  * Martín-Martínez  _et al._ (2013) E. Martín-Martínez, A. Dragan, R. B. Mann, and I. Fuentes, New Journal of Physics 15, 053036 (2013).
  * Sabín  _et al._ (2014) C. Sabín, A. White, L. Hackermuller, and I. Fuentes, Scientific reports 4, 6436 (2014).
  * Schirotzek  _et al._ (2009) A. Schirotzek, C.-H. Wu, A. Sommer, and M. W. Zwierlein, Phys. Rev. Lett. 102, 230402 (2009).
  * Kohstall  _et al._ (2012) C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. M. Bruun, F. Schreck, and R. Grimm, Nature 485, 615 (2012).
  * Koschorreck  _et al._ (2012) M. Koschorreck, D. Pertot, E. Vogt, B. Fröhlich, M. Feld, and M. Köhl, Nature 485, 619 (2012).
  * Massignan  _et al._ (2014) P. Massignan, M. Zaccanti, and G. M. Bruun, Reports on Progress in Physics 77, 034401 (2014).
  * Schmidt  _et al._ (2016) R. Schmidt, H. R. Sadeghpour, and E. Demler, Phys. Rev. Lett. 116, 105302 (2016).
  * Lan and Lobo (2014) Z. Lan and C. Lobo, J. Indian I. Sci. 94, 179 (2014).
  * Levinsen and Parish (2014) J. Levinsen and M. M. Parish, (2014), arXiv:1408.2737 .
  * Schmidt  _et al._ (2012) R. Schmidt, T. Enss, V. Pietilä, and E. Demler, Phys. Rev. A 85, 021602 (2012).
  * Płodzień  _et al._ (2018) M. Płodzień, R. Demkowicz-Dobrzański, and T. Sowiński, Phys. Rev. A 97, 063619 (2018).
  * Rath and Schmidt (2013) S. P. Rath and R. Schmidt, Phys. Rev. A 88, 053632 (2013).
  * Côté  _et al._ (2002) R. Côté, V. Kharchenko, and M. D. Lukin, Phys. Rev. Lett. 89, 093001 (2002).
  * Massignan  _et al._ (2005) P. Massignan, C. J. Pethick, and H. Smith, Phys. Rev. A 71, 023606 (2005).
  * Cucchietti and Timmermans (2006) F. M. Cucchietti and E. Timmermans, Phys. Rev. Lett. 96, 210401 (2006).
  * Lampo  _et al._ (2017) A. Lampo, S. H. Lim, M. Á. García-March, and M. Lewenstein, Quantum 1, 30 (2017).
  * Pastukhov (2017) V. Pastukhov, Phys. Rev. A 96, 043625 (2017).
  * Yoshida  _et al._ (2018) S. M. Yoshida, S. Endo, J. Levinsen, and M. M. Parish, Phys. Rev. X 8, 011024 (2018).
  * Guenther  _et al._ (2018) N.-E. Guenther, P. Massignan, M. Lewenstein, and G. M. Bruun, Phys. Rev. Lett. 120, 050405 (2018).
  * Lampo  _et al._ (2018) A. Lampo, C. Charalambous, M. Á. García-March, and M. Lewenstein, arXiv:1803.08946 (2018).
  * Landau and Pekar (1948) L. D. Landau and S. I. Pekar, Zh. Eksp. Teor. Fiz. (1948).
  * Devreese and Alexandrov (2009) J. T. Devreese and A. S. Alexandrov, Reports on Progress in Physics 72, 066501 (2009).
  * Alexandrov and Devreese (2009) A. Alexandrov and J. Devreese,  _Advances in Polaron Physics_, Springer Series in Solid-State Sciences (Springer, 2009).
  * Lieb and Yamazaki (1958) E. H. Lieb and K. Yamazaki, Phys. Rev. 111, 728 (1958).
  * Lieb and Thomas (1997) E. H. Lieb and L. E. Thomas, Comm. Math. Phys. 183, 511 (1997).
  * Frank  _et al._ (2010) R. L. Frank, E. H. Lieb, R. Seiringer, and L. E. Thomas, Phys. Rev. Lett. 104, 210402 (2010).
  * Anapolitanos and Landon (2013) I. Anapolitanos and B. Landon, Lett. Math. Phys. 103, 1347 (2013).
  * Lim  _et al._ (2018) S. H. Lim, J. Wehr, A. Lampo, M. Á. García-March, and M. Lewenstein, Journal of Statistical Physics 170, 351 (2018).
  * (39) M. Mehboudi, L. A. Correa, and A. Sanpera, (in preparation) .
  * Caldeira and Leggett (1983a) A. Caldeira and A. Leggett, Physica A: Statistical Mechanics and its Applications 121, 587 (1983a).
  * Caldeira and Leggett (1983b) A. Caldeira and A. Leggett, Annals of Physics 149, 374 (1983b).
  * Weiss (2008) U. Weiss,  _Quantum dissipative systems_ , Vol. 13 (World Scientific Pub Co Inc, 2008).
  * Correa  _et al._ (2017) L. A. Correa, M. Perarnau-Llobet, K. V. Hovhannisyan, S. Hernández-Santana, M. Mehboudi, and A. Sanpera, Phys. Rev. A 96, 062103 (2017).
  * Hovhannisyan and Correa (2017) K. V. Hovhannisyan and L. A. Correa, arXiv:1712.03088 (2017).
  * Langevin (1908) P. Langevin, C. R. Acad. Sci. Paris. 146 (1908).
  * Mazo (2002) R. M. Mazo,  _Brownian Motion: Fluctuations, Dynamics, and Applications_ (Clarendon Press, 2002).
  * Note (1) Note that we are not only assuming to work with a large dataset, but also that the estimator which maps measurement outcomes to temperature estimates is unbiased Braunstein and Caves (1994).
  * Tóth and Apellaniz (2014) G. Tóth and I. Apellaniz, J. Phys. A: Math. Theor. 47, 424006 (2014).
  * Braunstein and Caves (1994) S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 3439 (1994).
  * Monras (2013) A. Monras, arXiv:1303.3682 (2013).
  * Sherson  _et al._ (2010) J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Nature 467, 68 (2010).
  * Lindblad (1976) G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * Charalambous  _et al._ (2018) C. Charalambous, M. Garcia-March, A. Lampo, M. Mehboudi, and M. Lewenstein, arXiv preprint arXiv:1805.00709 (2018).


# Increasing the skill of short-term wind speed ensemble forecasts combining
forecasts and observations via a new dynamic calibration

Gabriele Casciaro<EMAIL_ADDRESS>Francesco Ferrari
<EMAIL_ADDRESS>Daniele Lagomarsino Oneto
<EMAIL_ADDRESS>Andrea Lira-Loarca
<EMAIL_ADDRESS>Andrea Mazzino<EMAIL_ADDRESS>
###### Abstract

All numerical weather prediction models used for the wind industry need to
produce their forecasts starting from the main synoptic hours 00, 06, 12, and
18 UTC, once the analysis becomes available. The six-hour latency time between
two consecutive model runs calls for strategies to fill the gap by providing
new accurate predictions having, at least, hourly frequency. This is done to
accommodate the request of frequent, accurate and fresh information from
traders and system regulators to continuously adapt their work strategies.
Here, we propose a strategy where quasi-real time observed wind speed and
weather model predictions are combined by means of a novel Ensemble Model
Output Statistics (EMOS) strategy. The success of our strategy is measured by
comparisons against observed wind speed from SYNOP stations over Italy in the
years 2018 and 2019.

###### keywords:

wind forecasting , probabilistic forecasting , dynamic forecast calibration ,
ensemble model output statistics , wind forecast based on real-time conditions
, Numerical Weather Prediction models

††journal: Energy

[DICCA]organization=DICCA, Department of Civil, Chemical and Environmental
Engineering. University of Genoa, addressline=Via Montallegro 1, city=Genoa,
postcode=16145, state=Genoa, country=Italy [INFN] organization=INFN, National
Institute of Nuclear Physics, Genoa section, addressline=Via Dodecaneso 33,
city=Genoa, postcode=16146, state=Genoa, country=Italy
[MALGA]organization=MALGA, Machine Learning Genoa Center. University of Genoa,
addressline=Via Dodecaneso 35, city=Genoa, postcode=16146, state=Genoa,
country=Italy

First attempt to couple in an efficient and economic way real-time data and
ensemble predictions.

First assessment of the added value of real-time observations in a wind
calibration.

Real-time data can be easily and economically ingested in an EMOS-based
calibration.

Ingestion of real-time data produces noticeable benefits vs. static
calibrations.

Real-time data provide added value to the whole wind predictive probability
density.

##  1 Introduction

Global cumulative installations of onshore and offshore wind are expected to
exceed 1 TW before 2025 (Global Wind Energy Council, 2021). This means that
the contribution of wind power in power systems is becoming increasingly
important. The downside is that detailed schedule plans and reserve capacity
must be properly set by power system regulators (Impram et al., 2020) facing
the intrinsic problem of the highly intermittent nature of wind, making this
very hard to predict. The accuracy of wind forecasts thus becomes an issue of
paramount importance for the wind industry.
In a recent work by Casciaro et al. (2021), a novel accurate Ensemble Model
Output Statistics (EMOS) strategy for calibrating wind speed/power forecasts
from an Ensemble Prediction System (EPS) has been proposed and its superiority
when compared against more parsimonious strategies in the 0-48 h look-ahead
forecast horizon clearly emerged. However, because all global weather models
start their run from analysis corresponding to the main synoptic hours 00, 06,
12, and 18 UTC, weather predictions (of any forecast horizons) necessarily
remain frozen for six hours. This limitation is in sharp contrast with the
needs of power system regulators, as well as of traders for marketing wind
energy, who need to adapt their strategies hour after hour in a quasi-
continuous way. It is thus very important to propose accurate strategies which
give fresh information on the wind speed in a given location continuously
evolving between two consecutive main synoptic hours. Proposing a strategy
with such characteristics is the main aim of the present paper. In plain
words, we propose a novel, parsimonious, dynamic EMOS strategy where the
parameters entering the EMOS predictive probability density function now also
depends on real-time (or quasi-real time) observed wind speed data. Our
strategy combines two well-known advantages for wind prediction by physical
methods (i.e. based on numerical weather prediction, NWP, models) and
statistical methods (i.e. statistical models based on time-series of past
measured observables and/or observation-driven machine-learning (ML) based
methods). While the former methods are suitable for predicting wind (and, more
generally, any meteorological observable) in the forecast horizon larger than,
say, 6 hours, the latter strategies turn out to have greater skills for short-
term forecasting, say, less than a few hours ahead. For a review of those
methods and discussions of their skills, see, e.g., Soman et al. (2010) and
the references therein. We will show that combining quasi real-time observed
data and predictions via a suitable EMOS strategy provides an optimal assembly
strategy which outperforms the single aforementioned strategies when acting
separately. This turns out to be the case both in terms of point indices (the
standard NMAE and correlation coefficient) and via suitable statistical
indices to assess the whole probability density function of the (calibrated)
ensemble forecast error.
For our dynamic strategy to be used operatively, the observed wind speed on
the site of interest must be available in real/quasi-real time, other than as
a record of past observations for implementing the static calibration.

The paper is organized as follows. In Sec. 2.1 we introduce the wind speed
observations dataset used both to train our calibration algorithm and for its
testing. Sec. 2.2 presents the forecast data from the Ensemble Prediction
System (EPS) of the ECMWF. Sec. 3.1 provides a quick review of the standard
Ensemble Model Output Statistics (EMOS). A very recent generalization of
standard EMOS to account for nonlinearities is described in Sec. 3.2. The
algorithm we propose to ingest real-time observed data in a EMOS strategy is
described in Sec. 4.1. The new algorithm is tested in different respects in
Secs. 4.2 and 4.3. The final section is devoted to draw some conclusions and
perspectives.

##  2 Wind data

###  2.1 Observed data from SYNOP stations

From 2018 to 2019, SYNOP meteorological data were collected at 42 locations
across Italy. According to ICAO specifications, the SYNOP anemometers record
wind speed (knots) as an average over 10 minutes at a nominal measurement
height of 10 meters a.g.l. (ICAO, 2007) with hourly frequency. The stations
are arranged fairly uniformly across the peninsula, as shown in Figure 1. Data
have been split using 2018 as training set and 2019 as test set.

Figure 1: Location of a selected group of Italian SYNOP stations used in the
present study for forecast calibration and testing. Colors are coded according
to the orography elevation.

###  2.2 Forecast data from the ECMWF Ensemble Prediction System (EPS)

The Ensemble Prediction System (EPS) from the European Centre for medium-range
weather forecasts (ECMWF), is made up of 51 members: a control forecast with
no perturbations and 50 forecasts created by adding small perturbations to the
best-known initial condition, using a mathematical formulation based on
singular vector decomposition and stochastic parameterization to represent
model uncertainties (Buizza, 1995; Leutbecher and Palmer, 2008).
The EPS used in this study has spectral triangular truncation with a cubic-
octahedral grid Tco639 and 91 layers with a top of atmosphere pressure of 0.01
hPa (Buizza, 2018) and a resolution of about 18 km (Persson, 2001).
Although four runs per day are made available from the ECMWF, in the present
study the sole EPS runs starting from the 00 UTC analysis have been considered
in the first 48-hour look-ahead forecast horizon.
Although there are robust evidences (Molteni et al., 1996; Montani et al.,
2019) that observed mean and variance correlate well with the EPS ensemble
mean and variance, all EPS members frequently underestimate and are
underdispersive with respect to observations. Then EPS forecasts must be
properly calibrated, especially for close-to-surface meteorological data. The
surface wind speed is just a relevant example considered in the present paper.

##  3 The static calibration strategy

We provide in this section a brief review of the so-called Ensemble Model
Output Statistics (EMOS) through which a raw ensemble forecast can be
accurately transformed into a predictive probability density function, with a
simultaneous correction for biases (additive and multiplicative) and
dispersion errors (Gneiting et al., 2005; Thorarinsdottir and Gneiting, 2010).
A generalization recently proposed by Casciaro et al. (2021) accounting for
nonlinear relationships between predictands and both predictors, and other
weather observables used as conditioning variables will be also reviewed.
Both strategies share the common feature of being static, meaning that they do
not exploit real-time (or quasi-real time) observed data.

###  3.1 The standard Ensemble Model Output Statistics

The standard EMOS strategy, here denoted as EMOS0, is a generalization of
multiple linear regressions, also known as Model Output Statistics (MOS)
(Glahn and Lowry, 1972), which is commonly used to calibrate deterministic
forecasts.
Let us consider $M$ ensemble member forecasts, $X_{1},\cdots,X_{M}$, of a
univariate continuous, positive-defined variable $Y$, the meteorological
observable of interest, here the wind speed in one specific location and given
look-ahead forecast time.
The EMOS method uses a parametric distribution of the following general form:

$Y\mid X_{1},\cdots,X_{M}\thicksim f(Y\mid X_{1},\cdots,X_{M})$ (1)

where the left-hand side denotes the fact that the distribution is conditional
on the ensemble member forecasts.
Gneiting et al. (2006) proposed the truncated normal distribution (TN) as a
model for the wind speed conditional distribution $f$. An alternative
distribution, the log-normal (LN), has been proposed by Baran and Lerch (2015)
even if the Authors found that the TN-LN mixture model outperforms the
traditional TN. The gamma distribution was finally suggested by Scheuerer and
Möller (2015). All these options, and other variants, have been extensively
discussed and compared by Wilks (2018).
Among all possible choices, in the present study we selected the gamma
distribution, $\mathcal{G}(\mu,\sigma^{2})$, as recently exploited by Casciaro
et al. (2021) with excellent results. This distribution is fully described by
two parameters: the shape parameter $k$ and the scale parameter $\theta$.
Their expressions in terms of the mean $\mu$ and variance $\sigma^{2}$ read:
$k=\mu^{2}/\sigma^{2}$ and $\theta=\sigma^{2}/\mu$ with $\mu$ and $\sigma^{2}$
given by:

$\mu=a+b_{1}X_{1}+\cdots+b_{M}X_{M}$ (2)

$\sigma^{2}=c+dS^{2}.$ (3)

The coefficients $a,b_{1},\cdots,b_{M},c,d$ are non-negative parameters and
$S^{2}$ is the variance of the EPS, a quantitative measure of the ensemble
spread.
Gneiting et al. (2005) proposed a strategy based on the minimization of the
Continuous Ranked Probability Score (CRPS) to determine the EMOS0 free
parameters (Hersbach, 2000). In plain terms, the CRPS is defined as

$CRPS(F,Y)=\int_{-\infty}^{\infty}[F(x)-H(x-Y)]^{2}dx$ (4)

where $F$ is the cumulative probability function associated to the gamma
density function $\mathcal{G}$, $Y$ is the observation, and $H$ is the
Heaviside function, which returns 0 when $x<Y$ and 1 otherwise.
A closed form for the CRPS for the gamma distribution has been obtained by
Scheuerer and Möller (2015) making the minimization procedure easy and fast.
For an observation-forecast pair ($Y,\mathbf{X}$) it reads:

$crps=Y\left[2P\left(k,\frac{Y}{\theta}\right)-1\right]-k\theta\left[2P\left(k+1,\frac{Y}{\theta}\right)-1\right]-\frac{\theta}{\beta\left(\frac{1}{2},k\right)}$
(5)

with $Y$ being the observation, $P$ the incomplete gamma function (Abramowitz
and Stegun, 1948), and $\beta$ the beta function. The forecast vector
$\mathbf{X}=(X_{1},\cdots,X_{M})$ comes into the expression (5) via the
parameters $k$ and $\theta$. The quantity to be minimized in a training set
where both observations and forecasts are available is:

$CRPS=\frac{1}{N}\sum_{i=1}^{N}crps(\mathbf{X}_{i},Y_{i})$ (6)

with $i$ denoting the i-th observation-forecast pair and $N$ is the total
number of pairs in the training set, here corresponding to the whole year
2018..
The CRPS combines calibration and informativeness in one index, thus allowing
the evaluation of predictive performance that is based on the paradigm of
maximizing the sharpness of the predictive distributions subject to
calibration (Gneiting et al., 2007).

###  3.2 Accounting for non-linearities and small-scale dynamic effects: an
evolution of the standard EMOS strategy

Very recently (Casciaro et al., 2021) proposed a novel EMOS strategy and
tested it against field measurements. The new strategy turned out to be
largely superior with respect to the standard EMOS. The main strengths of the
approach, baptized EMOS+4r by Casciaro et al. (2021), is that nonlinear
features can be easily and economically accounted by conditioning
meteorological variables. The strategy also deals with the issue related to
the model grid-point best representing weather conditions observed at the
ground station. This issue, particularly severe in regions having relevant
small-scale orographic variations, arises one to the coarse spatial resolution
of the EPS (as well as of any general circulation model).
To summarize the main idea of the method, let us consider $X_{i}^{j}$ the i-th
ensemble member forecast on the j-th model grid-point ($j=1,\cdots,4$ spans
over the nearest model grid-points to the ground station) and $S^{2\>j}$ the
members variance on the j-th model grid-point. We also denote by
$Z_{1},\cdots,Z_{q}$ the $q$ categorical variables (e.g., the wind direction,
the hour of the day, the height of the boundary layer, among others reported
by Casciaro et al. (2021)) expected to be useful to disentangle the EPS
forecast error.
The EMOS+4r predictive distribution is exactly as in the standard EMOS0 apart
the key fact that now the free parameters are best-fitted for each combination
of classes’ levels, via a training set, by minimizing the CRPS.
In plain terms, mean and variance of the predictive distribution are now given
by:

$\mu=a(Z_{1},\cdots,Z_{q})+\sum_{i=1,j=1}^{M,4}b_{ij}(Z_{1},\cdots,Z_{q})X_{i}^{j}$
(7)
$\sigma^{2}=c(Z_{1},\cdots,Z_{q})+\sum_{j=1}^{4}d_{j}(Z_{1},\cdots,Z_{q})S^{2\>j}$
(8)

where $j$ spans over the 4 model grid points around the station.
The last step of the strategy is to perform a final EMOS0 downstream of
previous calibration steps with the aim of allowing a synchronization of the
forecast to the current climate trend. This aim is achieved in terms of a
rolling training over the past 40 days without conditioning.

##  4 The dynamic calibration

The main strengths of the EMOS+4r strategy can be summarized as follow: i) it
accounts for nonlinear dependencies between predictands and both predictors
and other weather observables used as conditioning variables; ii) it deals
with the issue of the model grid-point selection best representing the weather
at the ground station; iii) it allows the calibrated forecast to be aligned
with the current weather climate via a rolling training.
These remarkable strengths are accompanied by the weak point of being unable
to ingest real-time weather information. These these could be used latter on
to extrapolate ahead in time the model error which can be quantified at the
present time if observations are made available outside of the synoptic hours.
Once the model error is extrapolated to the future via some suitable data-
driven statistical methods, it can be exploited to further correct forecast
error in the nearest future. How much the correction will be effective far
from the hour at which the real-time observation is available is an issue to
be addressed in the following sections.

###  4.1 The proposed algorithm

This section illustrates how the potential benefit of real-time observations
can be optimally integrated with the EMOS+4r strategy giving birth to a new
dynamic calibration strategy which will be also tested against different
alternative approaches.

In order to ensure greater clarity and better understanding, we report in
Figure 2 the kernel of the calibration, i.e. a general EMOS having N input
predictors for the ensemble mean, K input predictors for the ensemble
variance, and q conditioning variables. The values of N, q, and K will vary
from case to case in the different steps of the calibration we will show in
the following.

Figure 2: Schematic view of the EMOS kernel of the dynamic calibration:
$\mathbf{Z}=(Z_{1},\cdots,Z_{q})$ are $q$ categorical conditioning variables
entering in the model parameters $a,b,c$, and $d$;
$\mathbf{X}=(X_{1},\cdots,X_{N})$ are $N$ input predictors for the mean,
$\mu$, of the EMOS predictive distribution; $\mathbf{S}=(S_{1},\cdots,S_{K})$
are $K$ input predictors for the variance, $\sigma$, of the EMOS predictive
distribution.

The first step (here denoted by Step 0) is to apply the EMOS kernel to each of
the 4 nearest model grid-points to the SYNOP station without using
conditioning variables. This step is schematically depicted in Figure 3 for a
representative model grid-point.

Figure 3: Schematic representation of the Step 0 of our dynamic calibration in
one of the 4 nearest model grid-points to a station. Input:
$\mathbf{X}=(X_{1},\cdots,X_{M})$ are the $M$ members of the ensemble
forecasts; $\mathbf{S}=(\sigma_{raw})$ is the square root of the variance of
the ensemble member forecasts (K=1). Output: $\mu_{0}$ and $\sigma_{0}$ the
parameters of the calibrated predictive distribution.

By applying it to each of the four model grid-points, we end up with
$\mu_{0}^{1},\cdots,\mu_{0}^{4};\sigma^{1}_{0},\cdots,\sigma^{4}_{0}$, with
the upper indices being associate to the model grid-point among the four
nearest to the ground station. The purpose of this first step is simply to
reduce the forecast BIAS while maintaining unchanged the correlation between
forecasts and observations. The resulting calibrated parameters from Step 0
are used as input to the second step of our calibration, Step 1.
New players $(O_{1},\cdots,O_{P})$ come into Step 1 and involve $P$ observed
data available at a reference hour, say $h$ (i.e. $h$ is the last hour at
which observed data are available and usable). In way of example, $P=2$ with
$O_{1}$ being the persistence built from the wind speed known at hour $h$, and
$O_{2}$ is the diurnal-cycle-based persistence, i.e. a forecast built in terms
of the observations occurred in the past 24 hours. Figure 4 summarizes how
observations are combined with the results from Step 0 carried out for the
four closest-to-the-station model grid-points.

Figure 4: Schematic representation of Step 1 of our dynamic calibration
involving the information coming from Step 0. Input:
$\mathbf{X}=(\mu_{0}^{1},\cdots,\mu_{0}^{4},O_{1},\cdots,O_{P})$, where
$\mu_{0}^{1},\cdots,\mu_{0}^{4}$ come from Step 0 iterated for the four
closest-to-the-station model grid-points, and $(O_{1},\cdots,O_{P})$ are $P$
variables related to station observations available at hour $h$;
$\mathbf{S}=(\sigma_{0}^{1},\cdots,\sigma_{0}^{4})$ comes from Step 0.
Conditioning variables: $\mathbf{Z}=(H)$ is the hour of the day (hourly step
within the 24 hours). Output: $\mu_{1}$ and $\sigma_{1}$ the calibrated
parameters of the Step 1 predictive distribution.

From Step 1 we end up with the parameters $\mu_{1},\sigma_{1}$ of the
predictive distribution. This latter now becomes a function not only of the
ensemble members but also of suitable observations.
Having used the hour of the day as conditioning variables, the information on
how far one is from the last available observation (at hour $h$) indirectly
comes into the model.
We now need to join all previous steps in terms of suitable conditioning
variables. This is the Step 2 of our calibration as summarized in Fig. 5.

Figure 5: Schematic representation of Step 2 of our dynamic calibration.
Input:
$\mathbf{X}=(\mu_{0}^{1},\cdots,\mu_{0}^{4},O_{1},\cdots,O_{P},\mu_{1})$ where
$\mu_{0}^{i}$ ($i=1,\cdots,4$) comes from Step 0, $(O_{1},\cdots,O_{P})$ are
$P$ station observations available at hour $h$; $\mu_{1}$ comes from Step 1;
$\mathbf{S}=(\sigma_{0}^{1},\cdots,\sigma_{0}^{4},\sigma_{1})$ where
$\sigma_{0}^{i}$ comes from Step 0, $\sigma_{1}$ comes from Step 1.
Conditioning variables: $\mathbf{Z}=(Z_{1},\cdots,Z_{q})$, selected following
the same procedure as in Casciaro et al. (2021). Output: $\mu_{2}$ and
$\sigma_{2}$, the Step 2 calibrated parameters of the predictive gamma
distribution.

The result of this last step of calibration is the calibrated predictive
probability gamma function, $\mathcal{G}(\mu_{2},\sigma_{2})$, which selects
the best model grid-point performance, disentangles the model error in terms
of suitable conditioning variables and, finally, synchronizes the forecast
with the most recent observed wind speed at the station. We however found (see
Sec. 4.2 for a quantitative analysis) that one more step is needed in order to
extend ahead in time as much as possible the benefit of real-time
observations. This final step (Step 3) is sketched in Figure 6. It simply
joins the static calibration EMOS+4r by Casciaro et al. (2021), with
parameters $\mu_{s},\sigma_{s}$, (where the subscript ’s’ is meant to stress
out that such parameters are from a ‘static’ calibration) with those coming
from Step 2. The final result is the gamma probability density function,
$\mathcal{G}(\mu_{d},\sigma_{d})$, with $\mu_{d}$ and $\sigma_{d}$ coming from
Step 3 of Fig. 6.

Figure 6: Schematic representation of Step 3, the last step of our dynamic
calibration. Input: $\mathbf{X}=(\mu_{2},\mu_{s})$ and
$\mathbf{S}=(\sigma_{2},\sigma_{s})$ where $\mu_{s}$ and $\sigma_{s}$ come
from the EMOS+4r static calibration by Casciaro et al. (2021) and reviewed in
Sec. 3.2. The parameters $\mu_{2}$ and $\sigma_{2}$ come from Step 2 of our
calibration. Conditioning variables: $\mathbf{Z}=(H)$ as in Fig. 4. Output:
$\mu_{d}$ and $\sigma_{d}$ the parameters of the final calibrated gamma
predictive probability density function.

###  4.2 Quantifying the importance of the different steps of calibration

Let us start to quantify the added value brought by single steps of our
dynamic calibration by assuming that observed data are available at h = 9 UTC.
Here $P=1$ and $O_{1}$ is the persistence built from the wind speed known at
hour $h$. In order to quantify the benefit brought by Step 0 of our
calibration, we compare in Fig. 7 the results of our dynamic calibration
obtained with and without Step 0, respectively. This is done in terms of the
skill score (see A) where the skill of the calibrated forecast from the
complete calibration (i.e. starting from Step 0) is comparatively quantified
against the calibrated forecast from the dynamic calibration now starting
directly from Step 1 using the $EMOS_{+4r}$ mean and variance as input. As one
can see from the figure, apart the first hour of forecast, there is a clear
added value brought by Step 0, as quantified both in terms of the NMAE index
and in terms of the correlation coefficient.

Figure 7: Mean skill score (of the skill scores computed at each station) for
both the NMAE index and the correlation coefficient, as a function of the
look-ahead time, for the complete dynamic calibration (i.e. the one starting
from Step 0) where the observed data at h = 9 UTC are used. The reference
calibration against which the skill scores have been calculated is the same
complete dynamic calibration which starts from Step 1 using the $EMOS_{+4r}$
mean and variance as input. The abscissa refers to the 24-hour forecast
horizon starting from the hour h = 9 UTC where the observations are available,
so that look-ahead time equal to 1 corresponds the forecast time at the hour
h+1 UTC, and so on for the other abscissa values.

In Fig. 8 the importance of Step 2 is quantified similarly to what we did in
Fig. 7: for both the NMAE index and the correlation coefficient, the skill
score of the complete dynamic calibration is computed taking the forecast
calibrated from the dynamic calibration without Step 2 as a reference.

Figure 8: As in Fig. 7 but now the reference calibration is the one carried
out in terms of our dynamic calibration without Step 2.

The results from Fig. 8 clearly indicate the added value brought by Step 2,
especially in the first 5 hours of forecast.
We conclude by quantifying the role of Step 3. This is done in Fig. 9 where
the reference calibrated forecast is the one coming from the dynamic
calibration now without using Step 3. Also in this case the importance of the
selected Step 3 is evident.

Figure 9: As in Fig. 7 but now the reference calibration is the one carried
out in terms of our dynamic calibration without Step 3.

###  4.3 Assessing the dynamic calibration via statistical indices

Having justified the added value brought by single steps of our dynamic
calibration, we now pass to assess the performaces of the dynamic calibration
against different forecasts, both static, based on calibrated forecast from
NWP models, and statistical, built solely in terms of past observed data,
using persistence (i.e. the simplest way to build a forecast from a given
observation) and a ML-based prediction. Here we consider the simplest way to
produce a forecast from observed data: the wind speed known at hour $h$ is
maintained for all considered look-ahead forecast horizons. This is thus a
prediction based on persistence whose values are encoded in the variable
$O_{1}$ in Step 1 and 2 of our dynamic calibration with $P=1$.
Let us start by emphasizing the weakness of the static calibration
$EMOS_{+4r}$ in the first forecast hours, a fact that motived us to propose a
dynamic calibration. To do that, let us imagine to have real-time data
available at the hour h (h belongs to the first 24 hours from the 00 UTC) and
to have at our disposal the $EMOS_{+4r}$ static calibration starting from 00
UTC with a forecast horizons of 48 hours. The $EMOS_{+4r}$ skill score (for
both NMAE and correlation coefficient ${\cal C}$) is reported in Fig. 10 where
the abscissa refers to the 24-hour forecast horizon starting from the hour h
where the observations are available (so that an abscissa equal to one
corresponds to the forecast time at the hour h+1). Such observation is not
ingested by the static $EMOS_{+4r}$ calibration but is used here to build a
forecast based on persistence which serves as reference forecast to evaluate
the $EMOS_{+4r}$ skill score of Fig. 10. The obtained skill scores are
computed averaging over all stations and over all hours h in the first 24
hours starting from 00 UTC. The shaded areas represent the 50 $\%$ confidence
interval around the median thus providing us an idea on how the skill score
varies among the different stations and by varying the hour h at which data
are available. As one can clearly see from Fig. 10 persistence overcomes
$EMOS_{+4r}$ in the first forecast hour while rapidly deteriorating from the
second hour onwards. This holds true in terms of both NMAE and correlation
coefficient.

Figure 10: Skill score of NMAE and correlation coefficient of the $EMOS_{+4r}$
static calibrated forecast using the persistence at time h as a reference
forecast. Continuous lines: median of the skill scores computed by considering
all stations and all hours h within the first 24 hours from 00 UTC. Shaded
areas represent the 50$\%$ confidence interval of the skill scores, accounting
for the variability of the latter among stations and while varying the hour h
at ourly steps between 00 UTC and 24 UTC of the same day.

The performance of our dynamic calibration is shown in Fig. 11 where the
associated skill scores are shown by taking the persistence as a reference
forecast. As one can see, our strategy overcomes persistence from the first
hour onwards while maintaining the quality of the $EMOS_{+4r}$ static
calibration already detectable from Fig. 10 for sufficiently long forecast
horizons.

Figure 11: As in Fig. 10 but now the skill of our dynamic calibration is
assesed against persistence.

This last remark can be clearly detected from Fig. 12 where the skill score of
the dynamic calibration is now presented taking the static calibration from
the $EMOS_{+4r}$ as a reference. A remarkable added value brought from the
dynamic calibration clearly emerges in the first 6-7 hours, progressively
reducing as the forecast horizon increases. Interestingly, the benefit carried
by the ingestion of the observed data never disappears, even for the farest
look-ahead forecast times.

Figure 12: As in Fig. 10 but now the skill of our dynamic calibration is
assesed against the static $EMOS_{+4r}$ calibration.

Let us conclude the assessment in terms of point error indices by evaluating
the skill of our dynamic calibration against the prediction built from a ML-
based algorithm detailed in B.

Figure 13: As in Fig. 10 but now the skill of our dynamic calibration is
assessed against the ML observation-driven predictions. Also note that in this
comparison we have confined the attention to the sole case corresponding to h
= 9 UTC with a forecast horizon of 12 hours. The $50\%$ confidence interval
thus gives us an idea on the variability of the skill score from station to
station.

This latter only uses observed data available at hour h = 9 UTC with a
forecast horizon of 12 hours. Results are reported in Fig. 13 and show how the
dynamic calibration outperforms the observation-driven ML-based prediction.

We have till now considered point error indices. Let us now pass to assess the
predictive probability density function as a whole. This is done in Fig. 14 in
terms of the so-called reliability index ($\Delta$) proposed by Delle Monache
et al. (2006), and the sharpness $Sh50$ corresponding to the average width of
the central 50% prediction interval of a forecast probability distribution as
proposed by Gneiting et al. (2007). Details are reported in A.

In Fig. 14 the reference forecast is the static $EMOS_{+4r}$ calibration. The
increased skills of the dynamic calibration overcoming the static calibration
especially in the first forecast hours is evident, thus completing our quality
assessment.
A short summary of the most relevant error indices considered in the present
study are reported in Tab. 1 for our dynamic calibration and, for comparison,
for the static calibration and persistence. Results are presented for three
forecast time intervals starting from the hour h at which observations are
available. For the sake of example, the interval 1-2 hours means that the
reported indices are the average of the corresponding indices at hours h+1 and
h+2. A further average is performed over all possible h between 1 and 24 hours
(hourly step) and over all stations. The table confirms all conclusions drawn
in the paper on the superiority of our dynamic calibration against the static
one and, even more, against persistence.

Figure 14: As in Fig. 12 but now the skill scores are computed for $\Delta$ and the sharpness parameter $Sh50$. For the skill score of $\Delta$ the confidence interval (not shown) varies between -20 and 20 $\%$. Table 1: A final summary of relevant indices used to assess the skill of the calibrations. | Persistence |  | Static
---
calibration
| Dynamic
---
calibration
NMAE | 0.32 | 0.34 | 0.30 | 1-2 h
$\cal{C}$ | 0.76 | 0.71 | 0.78
Sh50 | - | 3.5 | 3.1
NMAE | 0.43 | 0.35 | 0.33 | 3-4 h
$\cal{C}$ | 0.62 | 0.71 | 0.74
Sh50 | - | 3.5 | 3.3
NMAE | 0.51 | 0.35 | 0.34 | 5-6 h
$\cal{C}$ | 0.52 | 0.71 | 0.72
Sh50 | - | 3.5 | 3.4

##  5 Conclusions and perspectives

A novel EMOS strategy, baptized $EMOS_{+4r}$, has been recently proposed by
Casciaro et al. (2021). This strategy differs from the standard EMOS in that
the free parameters entering the predictive EMOS probability density function
are categorical functions, and not constant as in the standard EMOS. This
allows the introduction of nonlinearities in the calibration strategy. These
categorical functions must be best-fitted for each combinantion of classes’
levels, via a training set, by minimizing the CRPS.
If, on one hand, the new calibration strategy brings a relevant added value
with respect to the standard EMOS calibration, on the other hand persistence
built in terms of an observation made available at the hour h turns out to be
more accurate at the next hour h+1, rapidly degrading at the next look-ahead
forecast times. This simple observation suggested us that persistence could
have been favorably ingested in the static $EMOS_{+4r}$ calibration with the
following two main aims: i) taking advantage of the known information from
observed data at hour h in between two consecutive synoptic hours; ii)
maintaining, and possibly increasing, the well-established superiority of
$EMOS_{+4r}$ over persistence after the hour h+1. Transforming this idea into
a new dynamic calibration strategy taking advantage of real-time observed data
has been the main result of the present paper.
At the same time, the new calibration maintains the well-established
superiority of $EMOS_{+4r}$ in the far forecast horizons over purely data-
driven statistical forecasts, including complex ML-based predictions.
The higher quality of the wind speed forecast over forecasts by static
approaches not ingesting real-time observations, and remaining frozen between
two successive six-hour separated synoptic hours, is expected to generate a
relevant added value for the wind power forecast with many useful applications
for the whole wind industry.
To conclude with some perspectives, our new strategy paves the way to a myriad
of other applications of interest for the green energy market. What we did for
the wind speed can indeed be easily generalized to other green sources,
including solar radiation and wave generation (Lira Loarca et al., 2021;
Ferrari et al., 2020; Besio et al., 2016). The key ingredient for the strategy
to be implemented is to have at disposal observed data in real/quasi-real
time, a requirement getting easier every day.

##  6 Acknowledgments

G.C. has been funded by the Italian bank foundation “Fondazione Carige”. A.M.
acknowledges the funding from the Interreg Italia-Francia Marittimo SICOMAR+
project (grant number D36C17000120006) and from the Compagnia di San Paolo
(Project MINIERA No. I34I20000380007). A.L.L has been funded by the Interreg
Italia-Francia Marittimo SINAPSI Project (grant number D64I18000160007). We
thank the Aeronautica Militare - Servizio Meteorologico - for providing us
with the SYNOP data as well as data from the EPS forecasts. Discussions with
Lorenzo Rosasco, Agnese Seminara, and Alessandro Verri are warmly
acknowledged.

##  Appendix A Statistical indices

The Skill Score (SS) index (Wilks, 2011) is used here to make the comparison
between different calibration strategies as quantitative as possible. It
assesses the performance of a given calibration by comparing its associated
statistical error index against the one corresponding to a reference forecast.
Namely,

$SS=\frac{A-A_{ref}}{A_{opt}-A_{ref}}$ (9)

where $A$ is the value of a suitable error index associated to the calibrated
forecast, $A_{ref}$ is the same as $A$ but relative to a reference forecast.
Finally, $A_{opt}$ refers to the optimal index value. A perfect calibration
yields SS=1, corresponding to the upper bound of SS. Values of SS smaller than
one (including negative values) indicate that the calibrated forecast is less
accurate than the reference one.
As far as the error indices are concerned ($A$ in Eq. (9)), here we consider
the normalized mean absolute error (NMAE), the correlation coefficient (${\cal
C}$), the so-called reliability index ($\Delta$) proposed by Delle Monache et
al. (2006), and the sharpness $Sh50$ corresponding to the average width of the
central 50% prediction interval of a forecast probability distribution as
proposed by Gneiting et al. (2007).
The NMAE is defined as:

$NMAE=\frac{\sum_{n=1}^{N}\left|X_{n}-Y_{n}\right|}{\sum_{n=1}^{N}Y_{n}}$ (10)

where $Y_{n}$ is the n-th observation and $X_{n}$ is the corresponding n-th
forecast (here corresponding to the mean of the 50 EPS ensemble) and $N$ is
the number of observation-forecast pairs in a given test set.
The correlation coefficient, $\cal{C}$, is a measure of linear dependence
between two variables and ranges from -1 to 1, with 1 representing the highest
correlation and -1 representing the highest anti-correlation (Wilks, 2011).
Quantitatively (Lee Rodgers and Nicewander, 1988),

${\cal{C}}=\frac{\sum_{n=1}^{N}(X_{n}-\overline{X})(Y_{n}-\overline{Y})}{N\sigma_{X}\sigma_{Y}}$
(11)

with:

$\sigma_{X}=\sqrt{\frac{\sum_{n=1}^{N}(X_{n}-\overline{X})^{2}}{N}}$ (12)
$\sigma_{Y}=\sqrt{\frac{\sum_{n=1}^{N}(Y_{n}-\overline{Y})^{2}}{N}}$ (13)

where $\overline{X}$ and $\overline{Y}$ are the mean values of $X$ and $Y$.
The goal of the probabilistic forecast, according to Gneiting et al. (2007),
is to maximize the sharpness of the predictive distribution subject to
calibration. Anderson (1996) and Hamill and Colucci (1997) proposed the use of
verification rank (VR) histograms to assess the calibration of ensemble
forecasts. VR histograms show the distribution of the ranks when the ranks of
the observations are pooled within the ordered ensemble forecasts. In a
calibrated ensemble, the observations and ensemble predictions should be
interchangeable, resulting in a uniform VR histogram. The continuous analogue
of the VR histogram is the probability integral transform (PIT) histogram
(Dawid, 1984; Diebold et al., 1997; Gneiting et al., 2007). The PIT value is
determined by the value of the predictive cumulative distribution function at
the verifying observation. For calibrated forecasts, the empirical cumulative
distribution function of PIT values should converge to the uniform
distribution.
The reliability index $\Delta$ was proposed by Delle Monache et al. (2006) to
quantify the deviation of VR histograms from uniformity. To quantify the
deviation from uniformity in the PIT histograms, we use here the following
definition of $\Delta$:

$\Delta=\sum_{i=1}^{m}\left|f_{i}-\frac{1}{m}\right|$ (14)

where $m$ is the histogram number of classes, each with a relative frequency
of $1/m$, and $f_{i}$ is the observed relative frequency in class $i$.
This index ranges from 0 to $+\infty$, with the bound zero corresponding to
optimality.
Sharpness is a property of forecasts that refers to the concentration of
predictive distributions. The sharper the forecasts, the more concentrated the
predictive distributions are, and the sharper the better, subject to
calibration.
To assess sharpness we average the widths of the central 50% prediction
intervals of forecast probability distributions at all time instants. Namely,

$Sh_{P}=\frac{1}{T}\sum_{t=1}^{T}\left[Q\left(f(t),\frac{1}{2}+\frac{P}{2}\right)-Q\left(f(t),\frac{1}{2}-\frac{P}{2}\right)\right]$
(15)

where $T$ is the total number of instants, $Q$ is the quantile of the forecast
distribution, $f(t)$ is the predictive distribution at time $t$ and $P$ is the
probability interval (here 50%).

##  Appendix B Technical Note on Machine Learning Models used as Benchmark

As a further benchmark, beside the naive Persistence Method, we use a Machine
Learning approach which is based only on time series of measured data. In last
years Machine Learning algorithms have been proliferating in the field of
data-driven wind speed prediction, establishing itself as a competitive
alternative to state-of-the-art statistical models in the field of time series
forecasting in general (Parmezan et al., 2019). Here we adopt the strategy
proposed in Lagomarsino Oneto et al. (2022), where wind time series are used
to map the forecast task into Supervised Learning problems which they solve
with a nonlinear kernel method (Bishop, 2016). In Lagomarsino Oneto et al.
(2022) they find that, depending on location, direction can be useful to
increase the prediction accuracy. To exploit this information they also
consider the two Cartesian components of the wind vector $u_{t},v_{t}$ which
can be easily calculated from the original signals of wind speed $s_{t}$ and
direction $\theta_{t}$ as $v_{t}=s_{t}\sin{(\theta_{t})}$ and
$u_{t}=s_{t}\cos{(\theta_{t})}$ (being the angle $\theta_{t}$ measured in
radiants). Following their steps we define Machine Learning models which are
able to learn, for each time $t$, the relation between the future value of the
wind speed at time $t+h$ and the past $\mu$ measurements (we will also refer
to $\mu$ as memory), i.e.

$\widehat{s}_{t+h}=\mathcal{F}(\eta_{t-\mu+1},\dots,\eta_{t})$ (16)

where $\mathcal{F}$ denotes the desired relation, $\widehat{s}_{t+h}$ our
prediction at horizon $h$ and $\eta_{t}$ the values of the input variables at
time $t$.
Learning function $\mathcal{F}$ can be seen as a regression problem in the
context of Supervised Learning. To solve the problem, input-output couples
$(\bm{\mathrm{x}}_{t},y_{t})_{t=1}^{n}$ can be defined as

$\bm{\mathrm{x}}_{t}=\big{[}\eta_{t-\mu+1},\dots,\eta_{t}\big{]}\in\mathbb{R}^{d}$
(17) $y_{t}=s_{t+h}$ (18)

where $n$ is the number of provided samples, $d=\mu\times k$ and $k$ is the
number of variables used to define $\eta_{t}$. Denote by
$\bm{\mathrm{{X}}}\in\mathbb{R}^{n\times d}$ (for a certain $d$) the lag
matrix whose rows are $\bm{\mathrm{x}}_{t}$ and
$\bm{\mathrm{y}}\in\mathbb{R}^{n}$ be the vector of the outputs with elements
$y_{t}$. Linear regression assumes that future wind behavior depends linearly
on its past trends: it aims to find coefficients
$\bm{\mathrm{\beta}}\in\mathbb{R}^{d}$ such that the following error term is
minimized

$\frac{1}{n}\sum_{i=1}^{n}\|\bm{\mathrm{x}}_{i}\bm{\mathrm{\beta}}-y_{i}\|^{2}=\frac{1}{n}\|\bm{\mathrm{{X}}}\bm{\mathrm{\beta}}-\bm{\mathrm{y}}\|^{2}$
(19)

and hence the formula $\bm{\mathrm{x}}_{\mathrm{new}}\bm{\mathrm{\beta}}$ is
used to infer how the wind will evolve in the future. To allow our models to
capture more complex dependencies we perform Kernel Ridge Regression (KRR)
introducing a non-linear transformation of the features via the kernel
function $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ (Murphy, 2012).
In our experiments we use the Gaussian kernel which is defined as
$k(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}}_{j})=e^{-\|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}_{j}\|^{2}/(2\sigma^{2})}$.
The solution to the KRR problem yields an estimator $\widehat{f}$ which can be
used for inference

$\widehat{f}(\bm{\mathrm{x}}_{\mathrm{new}})=k(\bm{\mathrm{x}}_{\mathrm{new}},\bm{\mathrm{{X}}})(\bm{\mathrm{{K}}}+n\lambda\mathrm{I})^{-1}\bm{\mathrm{y}}$
(20)

where $\bm{\mathrm{{K}}}\in\mathbb{R}^{n\times n}$ is the kernel matrix with
values $\bm{\mathrm{{K}}}_{ij}=k(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}}_{j})$ and
$\lambda$ is a regularization parameter which ensures the problem is well-
posed and must be tuned in order to optimize the predicting performance of the
model.
Since we focus on predictions delivered daily at 9:00 a.m. then only inputs
vectors whose reference times $t$ correspond to that time of day have been
kept. Furthermore we aim to predict the wind speed from $1$ to $12$ hours
ahead with an hourly frequency and for each horizon we consider different
combinations of memory $\mu$ and input data $\eta_{t}$. We analyse all
memories $\mu$ up to $24$ hours and two different cases for $\eta_{t}$. In the
first case $\eta_{t}=s_{t}$, i.e. the input data is the wind speed itself; in
the second case $\eta_{t}=(u_{t},v_{t})$: the two components of the wind
vector, which encodes the wind speed but also the direction. A reduction of
the dimension of input vectors have also been evaluated performing a Principal
Component Analysis on the input samples and keeping subsets of principal
components with increasing importance. Finally for each pair of horizon and
location we select the setup, in terms of memory and input design (either with
or without dimensionality reduction), that maximises the Pearson Correlation
between predictions and actually observed values of the speed. Data from each
station have been split into disjoint training and test set as described in
section 3.1. Again following Lagomarsino Oneto et al. (2022) a static rather
than a rolling approach is used for training. We then further split the
training set to estimate model hyperparameters $\lambda$ and $\sigma$ using
five-fold cross-validation. We performed an extensive 2-step grid-search to
find the hyperparameters which minimize the $R^{2}$ score for each separate
location. In the first step a coarse-grained grid was used to identify the
area in which hyperparameters were acceptable, in the second step we zoomed in
around such area to further refine the hyperparameter values.

## References

  * Abramowitz and Stegun (1948) Abramowitz, M., Stegun, I. A., 1948. Handbook of mathematical functions with formulas, graphs, and mathematical tables. Vol. 55. US Government printing office.
  * Anderson (1996) Anderson, J. L., 1996. A method for producing and evaluating probabilistic forecasts from ensemble model integrations. Journal of Climate 9 (7), 1518–1530.
  * Baran and Lerch (2015) Baran, S., Lerch, S., 2015. Log-normal distribution based Ensemble Model Output Statistics models for probabilistic wind-speed forecasting. Quarterly Journal of the Royal Meteorological Society 141 (691), 2289–2299.
  * Besio et al. (2016) Besio, G., Mentaschi, L., Mazzino, A., 2016. Wave energy resource assessment in the Mediterranean Sea on the basis of a 35-year hindcast. Energy 94, 50–63.
  * Bishop (2016) Bishop, C. M., 2016. Pattern Recognition and Machine Learning. Springer New York, New York, NY, softcover reprint of the original 1st edition 2006 (corrected at 8th printing 2009) edition.
  * Buizza (1995) Buizza, R., 1995. Optimal perturbation time evolution and sensitivity of ensemble prediction to perturbation amplitude. Quarterly Journal of the Royal Meteorological Society 121 (527), 1705–1738.
  * Buizza (2018) Buizza, R., 2018. Ensemble forecasting and the need for calibration. In: Statistical Postprocessing of Ensemble Forecasts. Elsevier, pp. 15–48.
  * Casciaro et al. (2021) Casciaro, G., Ferrari, F., Mazzino, A., 2021. Novel strategies of Ensemble Model Output Statistics (EMOS) for calibrating wind speed/power forecasts. arXiv preprint arXiv:2108.12174.
  * Dawid (1984) Dawid, A. P., 1984. Present position and potential developments: Some personal views statistical theory the prequential approach. Journal of the Royal Statistical Society: Series A (General) 147 (2), 278–290.
  * Delle Monache et al. (2006) Delle Monache, L., Hacker, J. P., Zhou, Y., Deng, X., Stull, R. B., 2006\. Probabilistic aspects of meteorological and ozone regional ensemble forecasts. Journal of Geophysical Research: Atmospheres 111 (D24).
  * Diebold et al. (1997) Diebold, F. X., Gunther, T. A., Tay, A., 1997. Evaluating density forecasts.
  * Ferrari et al. (2020) Ferrari, F., Besio, G., Cassola, F., Mazzino, A., 2020. Optimized wind and wave energy resource assessment and offshore exploitability in the Mediterranean Sea. Energy 190, 116447.
  * Glahn and Lowry (1972) Glahn, H. R., Lowry, D. A., 1972. The use of model output statistics (MOS) in objective weather forecasting. Journal of Applied Meteorology and Climatology 11 (8), 1203–1211.
  * Global Wind Energy Council (2021) Global Wind Energy Council, 2021. GWEC, Global Wind Report 2021, 2021-07-25.
URL https://gwec.net/wp-content/uploads/2021/03/GWEC-Global-Wind-
Report-2021.pdf

  * Gneiting et al. (2007) Gneiting, T., Balabdaoui, F., Raftery, A. E., 2007. Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69 (2), 243–268.
  * Gneiting et al. (2006) Gneiting, T., Larson, K., Westrick, K., Genton, M. G., Aldrich, E., 2006\. Calibrated probabilistic forecasting at the stateline wind energy center: The regime-switching space–time method. Journal of the American Statistical Association 101 (475), 968–979.
  * Gneiting et al. (2005) Gneiting, T., Raftery, A. E., Westveld III, A. H., Goldman, T., 2005\. Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation. Monthly Weather Review 133 (5), 1098–1118.
  * Hamill and Colucci (1997) Hamill, T. M., Colucci, S. J., 1997. Verification of Eta–RSM short-range ensemble forecasts. Monthly Weather Review 125 (6), 1312–1327.
  * Hersbach (2000) Hersbach, H., 2000. Decomposition of the continuous ranked probability score for ensemble prediction systems. Weather and Forecasting 15 (5), 559–570.
  * ICAO (2007) ICAO, 2007. Meteorological Service for International Air Navigation: Annex 3 to the Convention on International Civil Aviation, 16th edn. International Civil Aviation Organization: Montreal, Canada.
  * Impram et al. (2020) Impram, S., Nese, S. V., Oral, B., 2020. Challenges of renewable energy penetration on power system flexibility: A survey. Energy Strategy Reviews 31, 100539.
  * Lagomarsino Oneto et al. (2022) Lagomarsino Oneto, D., Meanti, G., Pagliana, N., Verri, A., Mazzino, A., Rosasco, L., Seminara, A., 2022. Physics informed shallow machine learning for wind speed prediction. In preparation.
  * Lee Rodgers and Nicewander (1988) Lee Rodgers, J., Nicewander, W. A., 1988. Thirteen ways to look at the correlation coefficient. The American Statistician 42 (1), 59–66.
  * Leutbecher and Palmer (2008) Leutbecher, M., Palmer, T. N., 2008. Ensemble forecasting. Journal of Computational Physics 227 (7), 3515–3539.
  * Lira Loarca et al. (2021) Lira Loarca, A., Ferrari, F., Mazzino, A., Besio, G., 2021. Future wind and wave energy resources and exploitability in the Mediterranean Sea by 2100. Applied Energy 302, 117492.
  * Molteni et al. (1996) Molteni, F., Buizza, R., Palmer, T. N., Petroliagis, T., 1996. The ECMWF ensemble prediction system: Methodology and validation. Quarterly Journal of the Royal Meteorological Society 122 (529), 73–119.
  * Montani et al. (2019) Montani, A., Gascon, E., Hewson, T., Paccagnella, T., Pillosu, F., Pincini, G., 2019. Performance of ECMWF- and COSMO-based ensemble forecast systems for precipitation events over Italy. EMS Annual Meeting Abstracts Vol. 16, EMS2019-146.
  * Murphy (2012) Murphy, K. P., 2012. Machine learning: a probabilistic perspective. MIT press.
  * Parmezan et al. (2019) Parmezan, A. R. S., Souza, V. M., Batista, G. E., 2019. Evaluation of statistical and machine learning models for time series prediction: Identifying the state-of-the-art and the best conditions for the use of each model. Information Sciences 484, 302–337.
  * Persson (2001) Persson, A., 2001. User guide to ECMWF forecast products. ECMWF Meteorological Bulletin M3 2.
  * Scheuerer and Möller (2015) Scheuerer, M., Möller, D., 2015. Probabilistic wind speed forecasting on a grid based on ensemble model output statistics. The Annals of Applied Statistics 9 (3), 1328–1349.
  * Soman et al. (2010) Soman, S. S., Zareipour, H., Malik, O., Mandal, P., 2010. A review of wind power and wind speed forecasting methods with different time horizons. In: North American Power Symposium 2010. IEEE, pp. 1–8.
  * Thorarinsdottir and Gneiting (2010) Thorarinsdottir, T. L., Gneiting, T., 2010. Probabilistic forecasts of wind speed: Ensemble model output statistics by using heteroscedastic censored regression. Journal of the Royal Statistical Society: Series A (Statistics in Society) 173 (2), 371–388.
  * Wilks (2011) Wilks, D. S., 2011. Statistical methods in the atmospheric sciences. Vol. 100\. Academic Press.
  * Wilks (2018) Wilks, D. S., 2018. Univariate Ensemble Postprocessing. In: Statistical Postprocessing of Ensemble Forecasts. Elsevier, pp. 49–84.


# Mathematical Musings of a Urologist

John R. Akeroyd; Robert K. Powers; Ganesh Rao, M.D

(Date: August 1, 2023
)

###### Abstract.

The derivatives with respect to the variable $r$ of $\pi r^{2}$ and
$\frac{4}{3}\pi r^{3}$ are $2\pi r$ and $4\pi r^{2}$, respectively. This
relates, through the derivative, the area enclosed in a circle to the length
of that circle and, likewise, the volume of a sphere to the surface area of
that sphere. The reasons why this works are basic to a first course in
calculus. In this brief article, we expand on these ideas to shapes other than
circles and spheres. Our approach is with the first year calculus student in
mind.

##  1\. Introduction and the Reasons Why

Urologists are natuarally interested in surface area versus rate of change in
volume of spherically shaped objects. We review connections between these
things without pretending that anything we say in this brief article is
fundamentally new. In fact, numerous articles have been written on this topic
– we mention two: [2] and [3]. Our main result, Theorem 1.4, is accessible to
the first year calculus student. As mentioned in the abstract, the
circumference of a circle of radius $r$ is the derivative with respect to $r$
of the formula for the area enclosed in that circle, and likewise with the
surface area of a sphere of radius $r$ and the formula for the volume of the
region bounded by that sphere. In the case of the circle of radius $r$ (see
[5]),

$\pi
r^{2}=\int_{0}^{2\pi}\int_{0}^{r}\rho\,d\rho\,d\theta=\int_{0}^{r}2\pi\rho\,d\rho.$

And in the case of the sphere,

$\frac{4}{3}\pi
r^{3}=\int_{0}^{\pi}\int_{0}^{2\pi}\int_{0}^{r}\rho^{2}\sin{\varphi}\,d\rho\,d\theta\,d\varphi=\int_{0}^{r}4\pi\rho^{2}\,d\rho.$

In the first case, the integral in terms of polar coordinates gives the area
of the region enclosed by a circle of radius $r$ as an integral of the
circumference formula for a circle of radius $\rho$, $0<\rho<r$. And in the
second, the integral in terms of spherical coordinates gives the volume of the
region bounded by a sphere of radius $r$ as an integral of the surface area
formula for a sphere of radius $\rho$, $0<\rho<r$. So, it is as basic as the
Fundamental Theorem of Calculus, and certainly no accident, that the
derivatives of the area and volume formulae here are circumference and surface
area formulae, respectively. In the case of a square of sidelength $x$, the
area enclosed is $x^{2}$ and the perimeter is $4x$. Clearly, the derivative of
$x^{2}$ does not equal $4x$. Things work out nicely, though, if we express
terms differently: In a way that is more akin to polar and spherical
coordinates. Given a square $S$, let $r$ be the radius of the largest circle
that can be inscribed in $S$. Then the sidelength of $S$ is $2r$ and the area
of the region enclosed by $S$ is $4r^{2}$. So, the derivative with respect to
$r$ of the area of the region enclosed by $S$ is $8r$, which is the length of
$S$ as a curve. And if $T$ is an equilateral triangle and $r$ is the radius of
the largest circle that can be inscribed in $T$, then one-half the base of $T$
is $\sqrt{3}r$ and the height of $T$ is $3r$; and so the area of the region
enclosed by $T$ is $3\sqrt{3}r^{2}$. Whence, the derivative with respect to
$r$ of the area enclosed by $T$ is $6\sqrt{3}r$, which is the length of $T$ as
a curve. This method works for any regular $n$-gon. We pose this as an
exercise. The reader may consult [2] for this result.

######  Exercise 1.1.

Let $W$ be a regular $n$-gon and let $r$ be the radius of the largest circle
that can be inscribed in $W$. Let $A(n,r)$ be the area of the region enclosed
by $W$ and let $L(n,r)$ be the length of $W$ as a curve. Show that the
derivative of $A(n,r)$ with respect to $r$ is $L(n,r)$.

These ideas extend to three dimensions in the cases of a cube and other solids
such as a tetrahedron, whose four faces are identical equilateral triangles.
In the three dimensional case, the radius of the largest inscribed sphere is
the variable $r$. For the cube, the sidelength is $2r$, the volume enclosed by
the cube is $8r^{3}$ and the surface area is $24r^{2}$, which is the
derivative with respect to $r$ of the volume formula here. The details for our
tetrahedron are more cumbersome, but the idea is the same. With these regular
curves in $\mathbb{R}^{2}$ and regular surfaces in $\mathbb{R}^{3}$ that are
topologically equivalent to a sphere, formulae for curve length, area
enclosed, surface area and volume enclosed are not difficult to calculate.
But, we really do not need regularity to makes things work. We first set our
terminology. For an $n$-gon, the edges are the line segments that join the
vertices and for a surface that is made up of planar regions, the planar
regions are called facets. So, a triangle has three edges and a tetrahedron
has four facets. The essential requirement that we need to make things work is
that there is a largest inscribed circle (or a largest inscribed sphere, in
the three dimensioinal case) and that this circle (or sphere) be tangent to
each edge (or each facet). The next two lemmas are straightforward, but they
get to the central issue here. We omit their proofs since they are well-known.
For a vector $\vec{x}=(x_{1},x_{2},...,x_{n})$ in $\mathbb{R}^{n}$, the
Euclidean norm of $\vec{x}$ is given by
$\|\vec{x}\|:=\sqrt{x_{1}^{2}+x_{2}^{2}+\cdot\cdot\cdot x_{n}^{2}}$.

######  Lemma 1.2.

Let $\mathcal{L}$ be a line in the Cartesian plane $\mathbb{R}^{2}$ that does
not pass through the origin and let $c$ be a positive real number, not equal
to $1$. Then $c\mathcal{L}:=\\{c\vec{x}:\vec{x}\in\mathcal{L}\\}$ is a line in
$\mathbb{R}^{2}$ that is parallel to $\mathcal{L}$, there exists $\vec{x_{0}}$
in $\mathcal{L}$ such that $\|\vec{x_{0}}\|\leq|\vec{x}\|$ for all
$\vec{x}\in\mathcal{L}$ and the distance between $\mathcal{L}$ and
$c\mathcal{L}$ is $|1-c|\|\vec{x_{0}}\|$.

######  Lemma 1.3.

Let $\mathcal{P}$ be a plane in $\mathbb{R}^{3}$ that does not pass through
the origin and let $c$ be a positive real number, not equal to $1$. Then
$c\mathcal{P}:=\\{c\vec{x}:\vec{x}\in\mathcal{P}\\}$ is a plane in
$\mathbb{R}^{3}$ that is parallel to $\mathcal{P}$, there exists $\vec{x_{0}}$
in $\mathcal{P}$ such that $\|\vec{x_{0}}\|\leq\|\vec{x}\|$ for all
$\vec{x}\in\mathcal{P}$ and the distance between $\mathcal{P}$ and
$c\mathcal{P}$ is $|1-c|\|\vec{x_{0}}\|$.

With regard to our lemmas, notice that the circle about the origin in
$\mathbb{R}^{2}$ of radius $\|\vec{x}_{0}\|$ lies tangent to the line
$\mathcal{L}$ and intersects $\mathcal{L}$ at $\vec{x_{0}}$. Moreover, any
circle about the origin of radius less than $\|\vec{x}_{0}\|$ has no
intersection with $\mathcal{L}$. And, similarly, the sphere about the origin
in $\mathbb{R}^{3}$ of radius $\|\vec{x}_{0}\|$ lies tangent to the plane
$\mathcal{P}$ and intersects $\mathcal{P}$ at $\vec{x_{0}}$. Moreover, any
sphere about the origin of radius less than $\|\vec{x}_{0}\|$ has no
intersection with $\mathcal{P}$. Let $W$ be an $n$-gon (not necessarily
regular) in $\mathbb{R}^{2}$ and suppose that there is a largest circle that
can be inscribed inside $W$ and that this circle lies tangent to each of the
edges of $W$. Without loss of generality, we may assume that this circle is
centered at the origin. For any positive real number $c$, let
$cW=\\{c\vec{x}:\vec{x}\in W\\}$. Then $cW$ is an $n$-gon that has edges
parallel to corresponding edges of $W$ and the distances between these
corresponding parallel edges are all the same, independent of the edge; cf.,
Lemma 1.2. Continuing with this $W$, the length of $W$ as a curve is the
product of the length of $\frac{1}{r}W$ with $r$ and the area enclosed by $W$
is the product of the area enclosed by $\frac{1}{r}W$ with $r^{2}$; where
$\frac{1}{r}W$ is the multiple of $W$ whose largest inscribed circle has
radius equal to $1$. So, with $W$ fixed, the length of the curve $W$ and the
area of the region enclosed by $W$ are functions of the variable $r$, alone.
Similar observations carry over to surfaces in $\mathbb{R}^{3}$ that are
topologically equivalent to a sphere and that are made up of finitely many
facets, each of which has piecewise smooth boundary.

######  Theorem 1.4.

Let $W$ be an $n$-gon in $\mathbb{R}^{2}$ and suppose that there is a largest
circle that can be inscribed inside $W$ and that this circle lies tangent to
each edge of $W$. Let $r$ be the radius of that circle, let $A(r)$ be the area
of the region enclosed by $W$ and let $L(r)$ be the length of $W$ as a curve.
Then the derivative of $A(r)$ with respect to $r$ is $L(r)$. Similarly, let
$\Omega$ be a surface in $\mathbb{R}^{3}$ that is topologically equivalent to
a sphere. Further, suppose that $\Omega$ has finitely many planar facets each
having piecewise smooth boundary and suppose that there is a largest sphere
that can be inscribed inside $\Omega$ and that this sphere lies tangent to
each facet of $\Omega$. Let $r$ be the radius of that sphere, let $V(r)$ be
the volume of the region bounded by $\Omega$ and let $\sigma(r)$ be the
surface area of $\Omega$. Then the derivative of $V(r)$ with respect to $r$ is
$\sigma(r)$.

###### Proof.

We establish the result in the cases that $W$ is a triangle and that $\Omega$
is a tetrahedron. These will be quite illustrative of how things go in
general.

Let $W$ be any triangle. Place a small circle inside $W$, near a corner of
$W$, so that it is tangent to the two edges of $W$ which form that corner.
Move the circle away from the corner and expand the circle so that it remains
tangent to these two edges until the circle meets the third edge of $W$. We
thus find a circle inscribed in $W$ that is tangent to all three edges of $W$;
which, by construction, is the largest inscribed circle in $W$. This is
sometimes called the inscribed circle, or incircle, of the triangle $W$. Let
$r$ be the radius of this circle. As before, we may assume that this circle is
centered at the origin in $\mathbb{R}^{2}$. For $0<s<r$, consider
$\frac{s}{r}W$, which is a triangle inside $W$ that has sides parallel to
corresponding sides of $W$ and has its largest inscribed circle equal to
$\\{\vec{x}:\|\vec{x}\|=s\\}$. Moreover, the gaps between the sides of $W$ and
the corresponding sides of $\frac{s}{r}W$ are all equal to $r-s$; see Figure.

$\swarrow$$\frac{s}{r}W$$\swarrow$$W$Figure

Let $A(r)$ be the area of the region enclosed by $W$ and let $A(s)$ be the
area of the region enclosed by $\frac{s}{r}W$. Let $L(r)$ be the length of the
curve $W$ and let $L(s)$ be the length of the curve $\frac{s}{r}W$. Notice
that

$L(s)[r-s]<A(r)-A(s)<L(r)[r-s]\,\,\,\mbox{and that}\,\,\,L(s)\rightarrow
L(r),\,\,\,\mbox{as}\,\,\,s\rightarrow r.$

Therefore,

$(A(r)-A(s))/(r-s)\rightarrow L(r),\,\,\,\mbox{as}\,\,\,s\rightarrow r.$

The same type of argument holds for $s>r$, and we find that the derivative of
$A(r)$ with respect to $r$ is $L(r)$.

We now let $\Omega$ be any tetrahedron. Our proof here mimicks the proof for
the triangle. Place a small sphere inside $\Omega$, near a corner of $\Omega$,
so that it is tangent to the three facets of $\Omega$ which form that corner.
Move the sphere away from the corner and expand the sphere so that it remains
tangent to these three facets until it meets the fourth facet of $\Omega$. In
this way we find an inscribed sphere inside $\Omega$ that is tangent to all
four facets of $\Omega$. By construction, this sphere is the largest inscribed
sphere in $\Omega$. We may assume that this sphere is centered at the origin
in $\mathbb{R}^{3}$ and we let $r$ be the radius of this sphere. Proceeding as
in the case of the triangle, but with $\Omega$ in place of $W$, let $V(r)$ be
the volume of the region enclosed by $\Omega$ and let $V(s)$ be the volume of
the region enclosed by $\frac{s}{r}\Omega$. Let $\sigma(r)$ be the surface
area of $\Omega$ and let $\sigma(s)$ be the surface area of
$\frac{s}{r}\Omega$. Then we have:

$\sigma(s)[r-s]<V(r)-V(s)<\sigma(r)[r-s]\,\,\,\mbox{and
that}\,\,\,\sigma(s)\rightarrow\sigma(r),\,\,\,\mbox{as}\,\,\,s\rightarrow r.$

Therefore,

$(V(r)-V(s))/(r-s)\rightarrow\sigma(r),\,\,\,\mbox{as}\,\,\,s\rightarrow r.$

The same type of argument holds for $s>r$, and we find that the derivative of
$V(r)$ with respect to $r$ is $\sigma(r)$.

∎

The proof of Theorem 1.4 outlines the argument in general and establishes the
result for any triangle (in the two dimensional case) and any tetrahedron (in
the three dimensional case). Unfortunately, the strategy of defining things in
terms of a largest inscribed circle or sphere does not work well in general.
It does not even work for nonsquare rectangles. For an approach that works in
greater generality one may consult [3] and articles that reference this paper.

There is another idea that generalizes to curves like rectangles, indeed to
most piecewise smooth Jordan curves (cf., [4] and [1, Section 4.2]), though
there is no nice parameter $r$, as we have above. Take a circle (sphere, in
the three dimensional case) of radius $\varepsilon>0$ and roll it around the
curve (or surface) so that its center traces out a curve (or surface) on the
inside that has gap, on normal lines, to the curve (or surface) in question,
equal to $\varepsilon$. In this way one can generate a family of curves (or
surfaces) which expand out to the curve or surface in question and one can
take a derivative (as $\varepsilon$ tends to zero) of the area (or volume) of
the region bounded by the curve (or surface). This certainly works for all
rectangles and many curves that are piecewise smooth. We outline this for a
rectangle with sidelengths $a$ and $b$. For small $\varepsilon>0$, roll a
circle on the inside of the rectangle, up against the curve, to trace out a
Jordan curve that is inside the rectangle. That Jordan curve, in this case, is
actually a rectangle of sidelengths $a-2\varepsilon$ and $b-2\varepsilon$;
provided $\varepsilon$ is sufficiently small. Taking the “derivative” of the
area,

$\lim_{\varepsilon\rightarrow
0}[ab-(a-2\varepsilon)(b-2\varepsilon)]/\varepsilon=2a+2b,$

we have the perimeter of our original rectangle. So our theme here prevails,
though this method is less natural than the earlier one.

Acknowledgement. We are grateful to John Duncan for comments that have helped
to improve the exposition of this paper.

## References

  * [1] H. Guggenheimer, Differential Geometry, Dover Publications, 1977.
  * [2] J. Miller, Differentiating area and volume, Two-year College Mathematics Journal 9:1 (1978) 47-48. Reprinted by the MAA in A Century of Calculus, Part II, 1992, pp. 134-135.
  * [3] J. Tong, Area and perimeter, volume and surface Area, College Math. J., Vol. 28, No. 1, January (1997), p. 57.
  * [4] J. Munkres, Topology, 2nd edition, Pearson, 2014.
  * [5] R. E. Williamson, H. F. Trotter. Multivariable Mathematics : linear algebra, calculus, differential equations. United Kingdom, Prentice-Hall, 1979. John R. Akeroyd
University of Arkansas
Robert K. Powers
Radford University
Ganesh Rao, M.D.
First Urology