Physics and Astrophysics of Planetary Systems

# The Atmospheres of Extrasolar Planets

Mark S. Marley  Mail Stop 245-3; NASA Ames Research Center; Moffett Field,
California 94035; USA

###### Abstract

The characteristics of irradiated solar system planetary atmospheres have been
studied for decades, consequently modern planetary science benefits from an
exhaustive body of ground- and space-based data. The study of extrasolar
planetary atmospheres, by contrast, is still in its infancy and currently
rests on a few score of datapoints, mostly of the transiting planets. This
short survey aims not to review this dynamic field but rather stresses the
importance of a few theoretical concepts and processes for our understanding
of exoplanet atmospheres. Topics covered include atmospheric structure and
dynamics, cloud processes and photochemistry of planetary atmospheres.
Influences on the albedos, spectra, and colors of extrasolar planets are
reviewed and caution is urged in the interpretation of exoplanet colors.

##  1 Introduction

The atmosphere controls a planet’s evolution through time and provides a
window into the chemical and physical conditions under which the planet
formed. Except for airless worlds, like Mercury, the atmosphere also mediates
the flow of information we receive about the nature of the planet. In this
contribution I aim to discuss a sampling of a few of the more important
concepts which relate to our understanding of extrasolar planetary
atmospheres. I will focus on topics that are important for understanding
spectral and photometric data from extrasolar planets, with an eye towards
providing examples from the solar system.

In some sense the trajectory of extrasolar planet science is recapitulating
the history of solar system planetary science. Three decades ago in the
University of Arizona Space Science series book “Jupiter” Wallace (1976)
reviewed the sparse data then available on Jupiter’s thermal emission spectrum
and presented models of the planet’s atmospheric structure. His effort to
piece together a consistent picture of the atmosphere from a handful of
thermal infrared photometric measurements, reflection spectra, and
transmission data from stellar occultations bears a striking resemblance to
current efforts to understand the transiting extrasolar planets.

Today, with our eyes now fixed on exoplanets, we are fortunate to have at hand
the insights gained from such work over the past half century of planetary
exploration. For a few transiting planets we again find ourselves in the
situation where we have a few broad band brightness temperature measurements
from which we are attempting to infer global atmospheric processes. As we
carry out our science we would do well to remember that this ground has been
crossed before: irradiated planetary atmospheres have been encountered before
the $21^{\rm st}$ century (a fact sometimes overlooked in the modern
literature). The excitement of our time, of course, is that the range of
parameter space, from planetary masses to the dizzying range of incident
radiation is far larger than previously encountered. Nevertheless we should
not be surprised that well understood processes familiar from solar system
planetary science, including photochemistry, hot stratospheres, and cloud
processes also play important roles in extrasolar planet atmospheres (Marley
1997).

To aid those attempting to traverse this new landscape of theory and data of
exoplanet atmospheres for the first time, I discuss in this chapter a few
elementary themes that are helpful for understanding the conceptual scenery.
While I will illustrate concepts by drawing from work in the field, I am not
in any sense attempting to provide a comprehensive review of exoplanet
science, as many such recent reviews are available (e.g., Charbonneau (2008),
Deming (2008), and Seager et al. (2008) review recent progress; Burrows et al.
(2001, 2006) and Barman et al. (2005) present useful surveys of much of the
theory of these objects; and Marley et al. (2007) review the pre-2007 field).
Instead I present a brief, illustrated (if somewhat idiosyncratic) guidebook
to a few of the important processes encountered at these new worlds. Hence I
will touch on atmospheric structure, cloud formation, the interpretation of
brightness temperature and the importance of photochemistry. I will mention
some of the processes contributing to stratospheric heating and conclude with
some brief comments about albedos and colors of planets.

##  2 Atmospheric Structure

The basic relationships governing a static atmosphere in hydrostatic
equilibrium are well known and presented in introductory textbooks, such as
Chamberlain & Hunten (1987). Such texts make for excellent introductory
reading as many of the fundamental relationships and processes governing
extrasolar planets are of course familiar from the solar system planets. One
particularly useful relation can be derived for an ideal gas atmosphere in
hydrostatic equilibrium. The number density of molecules per unit area in a
gaseous column above a specified height, $\cal N$, can be related to the local
number density of molecules, $n$, by

${\cal N}\approx{p(z)\over{g(z)m}}=n(z)H.$ $None$

Here the mean molecular weight is $m$, altitude is $z$ and $p,\,g$ and $H$ are
the pressure, gravity and the scale height; the equality is exact for an
isothermal atmosphere.

Figure 1: Atmospheric temperature-pressure profiles for Uranus, Jupiter, and
Earth. In all three planets the temperature increases with depth below a few
hundred millibars. In each planet this atmospheric region–the
troposphere–transports heat by convection from the deep interior, in the case
of the giants, or the surface in the case of Earth. In each atmosphere the
temperature also rises at low pressure–the stratosphere–owing to the
absorption of a fraction of the incident ultraviolet light by photochemical
products (ozone in the case of Earth, various hydrocarbon products in the
giants’ atmospheres). Uranus and Jupiter data are from the Voyager Radio
Science occultation experiments and the Earth profile is from the 1976
Standard Atmosphere. All data are available on-line
(http://atmos.nmsu.edu/planetary_datasets/indextemppres.html).

Equation (1) is useful for understanding a variety of basic characteristics of
atmospheres. For example if a molecular absorber has a cross section for
interaction with radiation of $\sigma\,$($\rm cm^{2}$), then the column
optical depth above a given pressure surface is $\tau_{\rm col}\sim{\cal
N}\sigma$. All else being equal, $\tau_{\rm col}$ falls linearly with
pressure. In a typical planetary atmosphere energy is transported upwards by
convection until the atmosphere becomes optically thin to thermal radiation,
when $\tau_{\rm col}\sim 1$. Above this level the outgoing energy is
transported by radiation. The pressure level in the atmosphere of this
radiative-convective boundary clearly depends upon the composition of the
atmosphere, the opacity of the major atmospheric constituents, gravity and the
temperature. For example, in the Earth’s atmosphere the surface temperature is
about $290\,\rm K$. Only at a temperature of about $220\,\rm K$ near $100\,\rm
mb$ (at the tropopause) is the optical depth at the peak of the Planck
function low enough that the air can radiate efficiently to space (Figure 1).
At lower pressures still a gray atmosphere would reach a constant temperature,
known as the skin temperature $T_{0}=(1/2)^{1/4}T_{\rm eff}$, where $T_{\rm
eff}$ is the effective temperature.

Figure 2: Pressure-temperature profiles for $\sim 4.5\,\rm Gyr$-old Jupiter-
like planets from 0.02 to 10 AU (left to right) from a solar-type star. Thick
lines are convective regions while thin lines are radiative regions. The
profiles at 5 and 10 AU show deviations that arise from numerical noise in the
chemical equilibrium table near condensation points, but this has a negligible
effect on planetary evolution. Figure adapted from Fortney et al. (2007).

The column abundance of molecules above a given pressure level in the
atmosphere controls the level at which the atmosphere becomes optically thin
to outgoing thermal radiation, and hence the location of the tropopause. The
tropopause pressure will thus vary with gravity as well as with atmospheric
structure and composition. In a lower gravity atmosphere each molecule
“weighs” less, so the column number of molecules, $\cal N$ above a given
pressure level must be larger to compress the gas than in a higher gravity
atmosphere. Indeed Equation (1) tells us that the column optical depth is
inversely proportional to $g$, so the column number density above the 1 bar
surface on a planet with half Earth’s gravity and a similar atmosphere would
be twice that of our atmosphere. All else being equal, the tropopause would be
at a lower pressure. For gas giant planets, which tend to have roughly
constant radii regardless of mass, higher mass generally corresponds to higher
gravity and thus more transparent atmospheres to a given pressure level. Given
the actual complexity of varying depths of absorption of incident radiation
with variations in gravity and atmospheric composition, theoretical models,
which relate atmospheric temperature to pressure, are required to fully
recognize the subtleties of atmospheric structure. Such models are shown in
Figure 2.

##  3 Atmospheric Dynamics

The overall structure of an irradiated atmosphere depends both upon the depth
at which incident energy is absorbed and internal sources of energy. For a
giant planet, convection transports energy outwards from the planetary
interior, as the planet cools slowly over time. For planets more distant than
a few AU from their primary star most incident energy is absorbed fairly deep
in the atmosphere, below the depth at which the atmosphere becomes optically
thick in the thermal infrared (Fortney et al. 2007). This is because most
gasses are more transparent in the optical than in the infrared. As a result
the absorbed incident energy simply adds to the internal energy being
transported outwards by convection and the temperature profile resembles that
of Jupiter shown in Figure 1. Another consequence of this relatively deep
deposition of solar radiation is that the internal heatflow is preferentially
transported by convection to the poles of Jupiter yielding a relatively
isothermal planet at the radiative-convective boundary (Ingersoll & Porco
1978).

However, for those giant planets found closer to their primary stars, the
radiative-convective boundary is deeper but absorption of incident flux still
occurs at a similar altitude (to the extent that composition is unchanged).
Thus the large incident flux upon a hot Jupiter is absorbed above the
radiative-convective boundary (Figure 2). As a consequence an isothermal layer
appears between the top of the deep convective zone and the region of the
atmosphere in which incident flux is absorbed. In this case the deep internal
heatflow is distinct from the thermalized incident radiation and the global
temperature distribution is no longer relatively homogeneous and equator to
pole temperature gradients can be large.

The nature of the atmospheric circulation of giant planets, which is
ultimately driven by the various energy fluxes, depends as well on a variety
of influences, particularly including rotation rate and atmospheric scale
height. Showman et al. (2007) present a very useful introduction and review of
this topic. The atmospheric redistribution of energy by winds is
unquestionably of paramount importance for the hot Jupiters and the efficiency
of redistribution controls the global temperature map and consequently the
phase variation of thermal emission, which has been successfully measured for
multiple planets by Spitzer (notably by Knutson et al. 2007; Cowan & Agol give
a complete summary of observations through mid-2008). Since a planet’s thermal
emission can arise from different depths in the atmosphere at different
wavelengths (§5), any variation in redistribution efficiency with altitude
will manifest itself as differing thermal emission maps as a function of
wavelength (Burrows et al. 2008a; Fortney et al. 2008). Ultimately coupled
models of radiative transfer and dynamics, similar to terrestrial global
circulation models, will be required to understand all of the contributing
factors (Showman et al. 2008 present one such model and review recent progress
in the field).

##  4 Clouds

Clouds play an important role in planetary atmospheres. They scatter incident
light back to space and sequester condensed species from the overlying
atmosphere. Jupiter’s atmosphere provides a point of departure for
understanding the diversity of giant planet atmospheres that we expect to
encounter outside of the solar system (see Lodders 2006 for a full
discussion).

Figure 3: Cloud structure expected to be found in, from left to right,
Jupiter’s atmosphere, a Jupiter-like planet too warm for water clouds to form,
and a hot Jupiter. Note that the overall cloud structure is similar in each
panel, but cloud levels move upward as the atmosphere is warmed. Important
labeled condensates include perovskite ($\rm CaTiO_{3}$) which removes TiO
from the gas phase, the salts $\rm Na_{2}S$ and KCl which remove two important
alkali absorbers (Figure 7) from the gas, and water clouds which appear near
the top of the Jupiter-like structure on the left. Figure modified from
Lodders (2004).

To understand Jupiter’s cloud structure we might imagine an air parcel moving
upwards from the deep interior (Figure 3). We start at a temperature of about
$2000\rm\,K$ and slowly raise the gas parcel up; as the gas rises, it cools
adiabatically. In the deep atmosphere all gasses are well mixed. At the point
in the atmosphere where a condensible species’ saturation mixing ratio first
equals the vapor mixing ratio, a cloud base is found. The first constituents
to condense are refractory oxides such as perovskite and corundum, followed by
various magnesium silicates including enstatite and forsterite. Iron is
predicted to condense as a native metal. Although we cannot see this region of
Jupiter’s atmosphere, we know that iron clouds are there because hydrogen
sulfide gas is detected in Jupiter’s visible atmosphere111If iron grains were
distributed uniformly above the iron condensation layer, $\rm H_{2}S$ gas
would react to form iron sulfide, FeS, thus removing sulphur-bearing gasses
from the atmosphere (Lodders 2004). (Niemann et al. 1998). As we move upwards
in the atmosphere the temperature continues to fall and eventually water
clouds form, removing $\rm H_{2}O$ from the gas phase. Above the water clouds,
the atmosphere continues to cool until ammonia clouds form. It is the ammonia
clouds of Jupiter, dusted by various photochemical pollutants, that we see
reflecting sunlight back from the planet. In §9 we’ll consider how giant
planets somewhat different than Jupiter might appear.

Since clouds scatter and absorb incident stellar radiation as well as emergent
thermal radiation, they play a very important role in controlling the
appearance and thermal structure of a planet. However clouds are intrinsically
difficult to model from a priori physical considerations. The detailed
behavior of terrestrial cloud cover as a function of atmospheric temperature
is the leading source of uncertainty in global climate models, for example.
Although the chemistry is thought to be well understood, predicting cloud
behavior for extrasolar planets, including such issues as particle sizes,
vertical distribution, and any horizontal patchiness is difficult.
Nevertheless some efforts at cloud modeling have been made (see the recent
review by Helling et al. (2008)). The hot L-type ultracool dwarfs have
atmospheres with thick silicate and iron clouds. Accounting for the effects of
these clouds has proven challenging. Given the difficulty of cloud modeling in
general and our experience with brown dwarfs, it seems that model predictions
for the spectra of extrasolar planets must be regarded skeptically, at least
when clouds are expected.

##  5 Brightness Temperature

The spectrum of any planet is composed of two components: scattered incident
radiation from the planet’s star and thermal emitted flux from the planet. The
thermal flux represents both energy arising from processes interior to the
planet and re-radiated absorbed incident radiation. For solar system planets
these two components of the spectrum are usually well separated in wavelength,
but for the hottest exoplanets there can be substantial overlap between
thermal radiation and scattered incident light. For a planet, such as a
transiting planet, with a known radius, the thermal emission spectrum is often
equated for convenience, wavelength by wavelength, to the thermal emission
from a blackbody. For an isothermal solid sphere with emissivity unity the
observed spectrum would equal that of a blackbody with a fixed temperature.
However for a real planet the flux will differ from that of a blackbody with
the same radius and at each wavelength a “brightness” temperature
$T_{B}(\lambda)$ may be defined.

Figure 4: Contribution functions calculated for a cloudless hot-Jupiter
lacking a stratosphere. Contribution functions are calculated for various
Spitzer broadband filters (the unlabeled curves correspond to various IRAC
filters), K band, and the Kepler band at 450–900 nm (black solid curve). For
clarity some of the curves have been normalized to 0.5 or 0 .75 rather than
1\. Figure from Showman et al. (2008).

Although thermal emission data for transiting planets are often reported in
terms of of $T_{B}(\lambda)$, such data must be regarded with some care.
Except in special cases (e.g., an isothermal atmosphere) brightness
temperature is not a measure of physical temperature or effective temperature.
Rather it gives a weighted measure of atmospheric temperatures over a range of
pressures from which flux emerges from the planet. To see this it is useful to
consider the expression for the upwards or outgoing intensity measured from an
atmosphere as a function of frequency, $I_{\nu}(0,\mu)$, where $\mu$ is the
cosine of the angle from the vertical. Making use of Eq. (1):

$I_{\nu}(0,\mu)={{\sigma_{\nu}n(z_{0})}\over\mu}\int\limits_{0}^{\infty}B_{\nu}(T)\exp{\biggl{(}}-{{z-z_{0}}\over
H}-{\tau_{\nu}\over\mu}{\biggr{)}}dz.$ $None$

The exponential function in this equation, known as a weighting function (see
Chamberlain & Hunten 1987), describes the relative contribution to the
outgoing flux from each altitude $z$, defined here as the distance from the
peak of the function, $z_{0}$. Crudely the brightness temperature can indeed
be equated to the temperature at the peak of the contribution function with
height. More importantly, the emergent flux, expressed as a single brightness
temperature is actually measuring thermal emission from a range of altitudes
in the atmosphere around $z_{0}$ with varying physical temperatures. Since
$\tau_{\nu}$ can vary dramatically with frequency $\nu$, the contribution
function—and thus the brightness temperature—can be quite different at
different wavelengths. For example Figure 4 illustrates theoretical
contribution functions for several different commonly used Spitzer filters
applied to the atmosphere of a cloud-free hot Jupiter (Showman et al. 2008).
Note that for each bandpass the measured emitted flux emerges over different
vertical regions of the atmosphere with different temperatures. In regions of
low opacity one sees deeper into the atmosphere which, for a monotonically
increasing temperature profile with depth, means higher brightness
temperature. High opacity spectral regions correspond to lower pressures and
lower temperatures. However if there is an inverted temperature profile then a
situation could emerge where $T_{B}(\lambda_{1})>T_{B}(\lambda_{2})$ and
$T_{B}(\lambda_{2})<T_{B}(\lambda_{3})$ where
$\tau(\lambda_{1})>\tau(\lambda_{2})>\tau(\lambda_{3})$. This is a commonplace
occurrence in the atmospheres of solar system giants and the chromospheres of
stars.

Furthermore, because of limb effects, this range in pressures from which flux
emerges also varies over the disk. Thus even for a gray atmosphere, with
constant optical depth as a function of wavelength, the brightness temperature
is in general not equal to the effective temperature at all wavelengths.

For these reasons, while brightness temperatures are useful shorthands to
convey information about planetary spectra, they must be regarded with some
caution. For example $8\,\rm\mu m$ Spitzer observation of the hot Neptune GJ
436b (Deming et al. 2007, Demory et al. 2007) yield a brightness temperature
of $712\pm 36\,\rm K$ which is modestly above the predicted effective
temperature (Deming et al. 2007). Since we do not expect, in general for
$T_{B}=T_{\rm eff}$ the information content of this single datapoint is
limited. With atmosphere models and additional data points the value of each
brightness temperature measurement increases, as was the case with Wallace’s
1976 study of Jupiter cited in the introduction.

##  6 Photochemistry

Figure 5: Incident flux at the top of the atmospheres of two transiting
planets compared to that received by Jupiter. Vertical lines denote the
approximate maximum wavelengths at which various molecules can be dissociated
(after Marley et al. 2007). Given the large incident fluxes received by these
and other hot Jupiters, many of the most abundant atmospheric species will be
easily photolyzed, which will likely lead to a rich photochemistry.

Atmospheric molecules can be dissociated by the absorption of ultraviolet
light, a process that happens high in the atmosphere before most incident UV
light is scattered back to space. Photochemical products can then participate
in complex reaction chains, producing various molecular products. A familiar
example is atmospheric ozone in Earth’s stratosphere, which ultimately results
from the photodissociation of molecular oxygen. Yung & DeMore (1999) provide a
useful guide to these topics on many planets. Photochemical products can
themselves become important players in the atmospheric radiative transfer of
giant planets, as discussed in the next section.

Photochemistry has long been expected (Marley 1998, Liang et al. 2004) to be
important for hot Jupiter atmospheres and will likely be far more complex than
in the solar system. This is because molecular species that are condensed
below the jovian clouds (e.g., $\rm H_{2}O,\ H_{2}S,\ NH_{3}$) and thus
protected from photodissociation will be gaseous in such hot atmospheres. Some
of these species, such as $\rm H_{2}S$, are easily photodissociated (Figure
5), and will likely produce new or unexpected species. Sulfur and nitrogen
compounds, in particular, may be important players in hot Jupiter
photochemistry and perhaps haze production (Marley et al. 2007). While the
carbon photochemistry has been studied (Liang et al. 2004), preliminary work
on photochemistry in water-bearing $\rm H_{2}-He$ atmospheres suggests that
compounds including $\rm CO_{2}$, HCN, and $\rm C_{2}H_{6}$ will be present
well in excess of the abundance predicted by equilibrium chemistry (Troyer et
al. 2007). Photochemical products may play a role in the formation of hot
stratospheres (Marley 1998; Marley et al. 2007; Burrows et al. 2008a). This
area is certainly rich for further study.

Photochemistry can also be very important in terrestrial atmospheres as is the
case with Venus and Titan. An earthlike planet with a greater abundance of
methane could well be enveloped in a photochemical haze and not appear
anything like a ‘pale blue dot’, even if there are indeed underlying oceans
(Zahnle 2008). Such obscuration by hazes is a major concern if planets are to
be characterized solely by their colors (§9).

##  7 Temperature Inversions (Stratospheres)

A gray atmosphere becomes asymptotically isothermal at small optical depths,
reaching the skin temperature $T_{0}$ (§2). It is often the case that an
atmosphere is optically thick to incident radiation at some wavelengths at low
pressures where the atmosphere is simultaneously relatively transparent at
infrared wavelengths. In this case more incident energy may be absorbed than
can be emitted by isothermal atmosphere with temperature $T_{0}$. As a result
the atmospheric layer with strong absorption must, in the absence of other
energy transport mechanisms, heat up until the thermal emission from the layer
equals the absorbed incident flux. An inverted temperature structure with a
warm, radiative upper atmosphere overlying a cooler tropopause, such as that
shown for several planets’ atmospheres in Figure 1, is known as a
stratosphere. Almost all solar system planets with an atmosphere exhibit a
stratosphere. In Earth’s atmosphere ozone absorbs ultraviolet light, which
warms the stratosphere to $270\,\rm K$, about 50 Kelvin warmer than the
temperature at the top of the troposphere. Solar system giant planet
atmospheres are heated by UV absorption by a combination of methane and
hydrocarbon photochemical products, including $\rm C_{2}H_{2}$ and $\rm
C_{2}H_{6}$ and photochemically produced hazes (e.g., see the review for
Neptune by Bishop et al. 1995).

The atmosphere of Jupiter provides a specific example. Without an energy
source the planet’s middle atmosphere would be close to 104 K, (the skin
temperature for Jupiter with $T_{\rm eff}=124\,\rm K$), as seen above the
tropopause in Figure 1. In the region where most of the incident UV flux is
absorbed (near 10 mbar) there is little overlap between a 100 K Planck
function and the important thermal opacity sources, so little flux can be
emitted. Since the absorbed energy cannot be radiated away by a 100 K
atmosphere, the atmosphere warms and the Planck function moves to shorter
wavelengths. Eventually the blue side of the Planck function overlaps the
strong $\nu_{4}$ methane fundamental and the $\nu_{9}$ ethane bands at 7.7 and
$12.2\,\rm\mu m$, allowing the atmosphere to radiatively cool, balancing the
absorbed incident flux (Chamberlain & Hunten 1987). As in other solar system
giant planet atmospheres, these strong mid-infrared bands of ethane and
methane act as a thermostat, regulating the stratospheric temperatures.

Likewise in exoplanet atmospheres a balance must be struck between the
absorption of incident radiation and thermal emission. For hot Jupiters which
are so warm that even the most refractory Ti- and V-bearing compounds do not
condense, TiO and VO gas may be exceptionally important absorbers (Hubeny et
al. 2003; Burrows et al. 2007; Fortney et al. 2008). These gasses, while not
abundant, have extraordinarily large absorption cross sections across the
entire optical spectrum. When present these molecules can absorb much of the
remarkably high incident flux above $1\,\rm mbar$ where the atmosphere is
optically thin in the thermal infrared. The atmosphere therefore becomes very
hot, perhaps as hot as $1900\,\rm K$ (Knutson et al. 2008), hot enough for
emission by the near-infrared and optical bands of water, CO, and even TiO to
balance this influx of energy (Fortney et al. 2008). There is a hint from
transit spectra of HD209458b that TiO and VO are indeed present in the
atmosphere (Desert et al. 2008).

As in solar system atmospheres, photochemistry may also play an important role
in exoplanet stratospheres. Photochemical hazes or photochemical gaseous
products that absorb well in the optical and UV could also provide prodigious
energy sources for exoplanet stratospheric heating. Photochemical pathways
have not yet been fully explored for these planets (but see the work on the
carbon chemistry by Liang et al. 2004 and speculation by Burrows et al. 2007)
and the opacity spectra for many potential molecules are not well known. This
area remains ripe for further study.

##  8 Albedos

Albedo (from the Latin word for ‘white’) is a measure of the reflectivity of
an object. In planetary science one encounters a variety of albedos (from
single-scattering to geometric to Bond to spherical, to name a few) and care
must be taken to carefully define the term in use. From a planet-wide
perspective the albedo of most importance is the Bond albedo, $A$, the ratio
of incident energy reflected into all angles by a planet to the total incident
energy. The Bond albedo appears in the equation for the equilibrium
temperature of a rapidly rotating planet with radius $R$ receiving an incident
flux ${\cal F}$:

$4\pi R^{2}\sigma T_{\rm eq}^{4}=(1-A)\pi R^{2}{\cal F}$ $None$

The geometric albedo, $p_{\lambda}$, is defined as the ratio of a planet’s
reflectivity measured at zero phase angle (opposition) to that of a Lambert
disk of the same radius. Unlike the Bond albedo, the geometric albedo is a
function of wavelength and, because it is measured at opposition (when the
phase angle $\Phi=0$), does not require information on the dependence of
scattering with phase to measure. For a perfectly reflecting Lambert sphere
the geometric albedo is $2\over 3$; for a semi-infinite purely Rayleigh
scattering atmosphere it is $3\over 4$. Both such idealized, perfectly
scattering objects would have a Bond albedo of 1, but the Rayleigh atmosphere
sends more light directly back to the observer at zero phase angle and thus
has a higher geometric albedo. Observed geometric albedo spectra for Uranus
and Jupiter as well as a model spectrum for HD 209458b are shown in Figure 6.
Note that the hot Jupiter is quite dark beyond about $0.4\,\rm\mu m$ since
gaseous Na and K absorb most incident photons before they can be scattered
back to space (see also Burrows et al. 2008b). Uranus and Jupiter would
likewise be quite dark in the red if not for their cloud layers that scatter
red photons before they can be absorbed by methane (Figure 7; Marley et al.
1999).

Figure 6: Observed geometric albedo of Uranus and Jupiter (Karkoschka 1994),
compared to one possible model of HD 209458b which ignores thermal emission,
and the $1\sigma$ upper limit from MOST (Rowe et al. 2007). The geometric
albedo of a deep Rayleigh scattering atmosphere and a Lambertian sphere are
shown as well. Figure from Fortney et al. (in prep.).

The Bond albedo is related to the monochromatic geometric albedo $q_{\lambda}$
by

$A=\int\limits_{0}^{\infty}p_{\lambda}q_{\lambda}f_{\lambda}d\lambda\Bigr{/}\int\limits_{0}^{\infty}f_{\lambda}d\lambda$
$None$

where $q_{\lambda}$ is the monochromatic phase integral and $f_{\lambda}$ is
the incident monochromatic flux. The phase integral is a measure of the
angular distribution of scattered light. For a Lambert sphere and a Rayleigh
sphere $q_{\lambda}={3\over 2}$ and $4\over 3$ respectively. In the general
case, however, a planet is neither a perfectly Lamberian nor a Rayleigh
scatterer and the scattered light must be measured as a function of phase
angle $\Phi$ and $\lambda$. This is why the Bond albedos for solar system
giants could only be accurately determined after the Voyager spacecraft had
measured their brightness at many different phase angles so that $q_{\lambda}$
could be measured (e.g. Pollack et al. 1986 for Uranus, see also the review by
Conrath et al. 1989).

Without a measurement of $q_{\lambda}$ it is common to assume that a planet
scatters light isotropically following a Lambertian phase function. The ratio
$C$ of the brightness of a planet at a distance $a$ from its star as seen in
reflected light at an arbitrary phase angle $\Phi$ to the brightness of the
star can then be written as

$C_{\lambda}(\Phi)=q_{\lambda}(R/a)^{2}\Bigl{[}{{\sin(\Phi)+(\pi-\Phi)\cos(\Phi)}\over\pi}\Bigr{]}.$
$None$

At quadrature $\Phi={\pi\over 2}$ and $C_{\lambda}=q_{\lambda}(R/a)^{2}/\pi$.
Note that for the giant planets in Figure 6 the contrast would be far more
favorable in the visible than in the red or infrared bands.

Figure 7: Model geometric albedo spectra for a Jupiter-mass planet at various
distances from a solar type star. Dotted lines are models for the 3 and 2 AU
planets but without cloud opacity. While these spectra are cast as a function
of orbital radius, the same sequence would result for a giant planet kept at a
fixed distance and modeled with progressively younger ages, higher masses, or
earlier stellar spectral types, all of which would produce warmer effective
temperatures, all else being equal. For example an $8\,\rm M_{J}$ planet at an
age of about 1 Gyr would have a similar spectrum to the planet at 1.0 AU. A
$4\,\rm M_{J}$ at the same age would be similar to a cloudless planet (dotted
line) at 2 AU. Figure adapted from Fortney et al. (in prep.).

Because molecular bands tend to be stronger at longer wavelengths and because
Rayleigh scattering is more efficient in the blue, most solar system giant
planets have larger geometric albedos at shorter wavelengths than at longer
wavelengths. Given that the Bond albedo is weighted by the incident flux,
which varies with stellar type, the same planet with the same geometric albedo
spectrum will have a variety of different Bond albedos as the color of the
incident flux is changed. Around an early type star, which is brightest in the
blue, a Jupiter-like planet might have a large Bond albedo of 0.45. But around
an M star the same planet’s Bond albedo might be less than 0.1 (Marley et al.
1999). For this reason it is usually best to consider the geometric albedo
which does not overtly depend on the incident flux (other than through its
dependence on the temperature profile).

##  9 Colors of Planets

Early efforts to directly image extrasolar planets in reflected light, such as
by small or intermediate-sized coronagraphs, will likely first produce images
in a few broad passbands. Given such limited available data it is worthwhile
to consider what might be learned about the nature of the detected planets
from such data.

Figure 8: Observed and model broad-band colors of giant planets. J, S, U, & N
are the observed colors of the corresponding solar system giants. Labeled
crosses show model colors of giant planets at the given distances from a
solar-type primary star and a simple cloud model. Unlabeled crosses are
Jupiter and Saturn models with a different set of assumptions about the cloud
properties. Depending on the assumed effective temperature and other model
details (cloud structure, etc.), the color of a $1\,\rm M_{J}$ planet can vary
widely. Color alone is thus not a useful discriminant of mass.

Given an orbital separation from the primary star, a single photometric
detection, combined with an assumed phase function and bounds placed on the
geometric albedo, would allow a crude estimate of the planet’s size. Assuming
an upper limit of $p<0.75$ (the pure Rayleigh scattering limit) and a lower
limit $p>0.06$ (typical of low albedo asteroids (Dotto et al. 2002)), for
example, would result in an uncertainty in the radius inferred for a directly
imaged planet of a factor of 3.5. A bright planet with a radius slightly
larger than Earth’s could not be distinguished from a dark planet with
Neptune’s radius on the basis of brightness alone. If the planet were also
detected by other means, for example radial velocity or astrometric methods,
then the known mass would discriminate between these two extremes. Without
such a detection, however, the nature of the planet would have to be discerned
by spectroscopic or photometric methods. Even low resolution spectroscopy
likely will be beyond the reach of modest aperture space-based coronagraphic
telescopes. This means that planets will have to be characterized, at least
initially, by their broadband colors.

Indeed based on our experience in the solar system, broadband colors of giant
planets at first seem to be promising markers for discerning planet type
(Figure 8 for the giants). Uranus and Neptune are blue while Jupiter and
Saturn are red. Among the terrestrial planets the Moon, Mercury, and Mars are
red while the Earth is slightly blue. It has been suggested (e.g., Traub 2003)
that such color trends, when applied to exoplanets, may help identify the type
of planets detected by direct imaging. A planet might be imaged in a few broad
spectral filters and characterized by comparison to solar system planetary
colors. There are substantial difficulties with this approach, however, since
a single planet can have very different colors just depending on its
temperature and the range of plausible colors is certainly much larger than
that sampled by solar system planets. The reflectivity of giant planets in the
blue is strongly influenced by stratospheric hazes (e.g., see Baines &
Bergstralh (1986) for Uranus) while the brightness in the red depends upon
cloud properties. As a result both photochemical hazes and clouds (Figure 7)
can substantially alter the broadband color of a planet and are difficult to
model on an a priori basis. Also there will certainly be exoplanets with
characteristics quite different from solar system objects (warm Neptunes,
super-Earths, water worlds, etc.) and their colors may be surprising.

For example Figure 7 plots the geometric albedo spectra of a Jupiter mass
planet at various distances from its primary star and thus warmer (see Marley
et al. (1999), Sudarsky et al. (2003), and Burrows (2005) for more complete
discussions). An extrasolar planet in the same orbit could also be warmer than
Jupiter if it were younger, more massive,222Since giants more massive than
Jupiter take longer to cool, a four Jupiter mass planet that is 2 billion
years old, would have $T_{\rm eff}\sim 400\,\rm K$ and exhibit a spectrum
similar to one of the lower curves of Figure 7. or orbiting a brighter and
hotter star than our sun. As we imagine warming Jupiter (Figure 3), first the
ammonia clouds would evaporate, allowing ammonia to be well mixed throughout
the atmosphere and allowing us to see the water clouds, which underlie the
ammonia clouds. Since the water clouds are optically thick and good scatterers
in the optical, such a planet would likely be much brighter and whiter than
our current Jupiter. As we continue to warm the planet the atmosphere would
heat up and the cloud base would move progressively higher in the atmosphere.
Eventually, even the massive water clouds would evaporate, and suddenly we
would have a clear atmosphere. Jupiter would then appear very dark blue
(Figure 7), because red photons, which do not efficiently Rayleigh scatter,
would burrow down into the atmosphere, never to return. As the atmosphere
continues to warm, alkali metals, which at low temperatures are found as
chlorides such as KCl, would also evaporate, again altering the optical
spectrum of the planet (Figures 6 & 7) as they absorb strongly in the red
(Burrows et al. 2000). The important lesson from this thought experiment is
that the same planet can have vastly different spectra–and by extension colors
(Figure 8)–depending on its effective temperature. It would be a mistake to
presume that one could recognize a Jupiter-mass planet simply by its broadband
color. In this example Jupiter varies from red to white to blue over the span
of a few hundred Kelvin. Spectra, such as those shown in Figure 7, along with
inferences on age of the primary star (hence placing limits on the age of the
planet) would be the best means of discriminating the nature of the planet.

##  10 Conclusions

As we enter the age of the direct detection and characterization of extrasolar
planets it is important not to overlook the lessons learned from the past half
century of planetary exploration. Solar system planets give ample
demonstration of the importance of atmospheric dynamics, cloud and
condensation processes, and photochemistry in controlling the face of planets.
The great diversity of atmospheres seen in our own solar system, from the hazy
skies of Titan to the turbulent atmosphere of Jupiter to the blue vistas of
Earth, are emblematic of the diversity of processes that can affect properties
of planets. Exoplanets will exhibit even larger ranges of properties and we
should not be surprised by hazy yellow Earths, red, white, or blue Jupiters,
or other unexpected worlds. Certainly extrasolar planet atmospheres will be
influenced by some of the processes mentioned here. (Other yet-to-be-
discovered influences will undoubtedly be important as well.) Ultimately,
however, such discoveries will extend the journey of planetary exploration
begun by the Mariners and Voyagers of the last century out into the galaxy.

Acknowledgements: The author thanks Jonathan Fortney and Kevin Zahnle for
helpful comments on the manuscript and Thierry Montmerle and the entire Les
Houches staff for organizing an exceptionally productive and informative
program.

## References

  * [Baines & Bergstralh(1986)] Baines, K. H., & Bergstralh, J. T. 1986, Icarus, 65, 406
  * [Barman et al.(2005)] Barman, T. S., Hauschildt, P. H., & Allard, F. 2005, Ap. J., 632, 1132
  * [Bishop et al.(1995)] Bishop, J., Atreya, S. K., Romani, P. N., Orton, G. S., Sandel, B. R., & Yelle, R. V. 1995, Neptune and Triton, 427
  * [Burrows(2005)] Burrows, A. 2005, Nature, 433, 261
  * [Burrows et al.(2000)] Burrows, A., Marley, M. S., & Sharp, C. M. 2000, Ap. J., 531, 438
  * [Burrows et al.(2001)] Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, Reviews of Modern Physics, 73, 719
  * [Burrows et al.(2006)] Burrows, A., Sudarsky, D., & Hubeny, I. 2006, Ap. J., 650, 1140
  * [Burrows et al.(2007)] Burrows, A., Hubeny, I., Budaj, J., Knutson, H. A., & Charbonneau, D. 2007, Ap. J. Lett., 668, L171
  * [Burrows et al.(2008)] Burrows, A., Budaj, J., & Hubeny, I. 2008a, Ap. J., 678, 1436
  * [Burrows et al.(2008)] Burrows, A., Ibgui, L., & Hubeny, I. 2008b, Ap. J., 682, 1277
  * [1987] Chamberlain, J. & Hunten, D. 1987, Theory of Planetary Atmospheres (Orlando: Academic Press)
  * [Charbonneau(2008)] Charbonneau, D. 2008, ArXiv e-prints, 808, arXiv:0808.3007
  * [Conrath et al.(1989)] Conrath, B. J., Hanel, R. A., & Samuelson, R. E. 1989, Origin and Evolution of Planetary and Satellite Atmospheres, 513
  * [Cowan & Agol(2008)] Cowan, N. B., & Agol, E. 2008, ArXiv e-prints, 806, arXiv:0806.4606
  * [Deming(2008)] Deming, D. 2008, ArXiv e-prints, 808, arXiv:0808.1289
  * [Demory et al.(2007)] Demory, B.-O., et al. 2007, Astron. & Astrophys., 475, 1125
  * [Desert et al.(2008)] Desert, J. -., Vidal-Madjar, A., Lecavelier des Etangs, A., Sing, D., Ehrenreich, D., Hebrard, G., & Ferlet, R. 2008, ArXiv e-prints, 809, arXiv:0809.1865
  * [Dotto et al.(2002)] Dotto, E., Barucci, M. A., Müller, T. G., Brucato, J. R., Fulchignoni, M., Mennella, V., & Colangeli, L. 2002, A & A, 393, 1065
  * [Fortney et al.(2007)] Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, Ap. J., 659, 1661
  * [Fortney et al.(2008)] Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, Ap. J., 678, 1419
  * [Helling et al.(2008)] Helling, C., et al. 2008, MNRAS in press, arXiv:0809.3657
  * [Hubeny et al.(2003)] Hubeny, I., Burrows, A., & Sudarsky, D. 2003, Ap. J., 594, 1011
  * [Ingersoll & Porco(1978)] Ingersoll, A. P., & Porco, C. C. 1978, Icarus, 35, 27
  * [Karkoschka(1994)] Karkoschka, E. 1994, Icarus, 111, 174
  * [Knutson et al.(2007)] Knutson, H. A., et al. 2007, Nature, 447, 183
  * [Knutson et al.(2008)] Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, Ap. J., 673, 526
  * [Liang et al.(2004)] Liang, M.-C., Seager, S., Parkinson, C. D., Lee, A. Y.-T., & Yung, Y. L. 2004, Ap. J. Lett., 605, L61
  * [Lodders(2004)] Lodders, K. 2004, Science, 303, 323
  * [Marley(1998)] Marley, M. S. 1998, Brown Dwarfs and Extrasolar Planets, 134, 383
  * [Marley et al.(2007)] Marley, M. S., Fortney, J., Seager, S., & Barman, T. 2007, Protostars and Planets V, 733
  * [Niemann et al.(1998)] Niemann, H. B., et al. 1998, J. Geophys. Res., 103, 22831
  * [Pollack et al.(1986)] Pollack, J. B., Rages, K., Baines, K. H., Bergstralh, J. T., Wenkert, D., & Danielson, G. E. 1986, Icarus, 65, 442
  * [Rowe et al.(2007)] Rowe, J. F., et al. 2007, ArXiv e-prints, 711, arXiv:0711.4111
  * [Seager et al.(2008)] Seager, S., Deming, D., & Valenti, J. A. 2008, ArXiv e-prints, 808, arXiv:0808.1913
  * [Showman et al.(2007)] Showman, A. P., Menou, K., & Y-K. Cho, J. 2007, ArXiv e-prints, 710, arXiv:0710.2930
  * [Showman et al.(2008)] Showman, A. P., Fortney, J. J., Lian, Y., Marley, M. S., Freedman, R. S., Knutson, H. A., & Charbonneau, D. 2008, ArXiv e-prints, 809, arXiv:0809.2089
  * [Sudarsky et al.(2003)] Sudarsky, D., Burrows, A., & Hubeny, I. 2003, Ap. J., 588, 1121
  * [Traub(2003)] Traub, W. A. 2003, Scientific Frontiers in Research on Extrasolar Planets, 294, 595
  * [Troyer et al.(2007)] Troyer, J., Moses, J. I., Fegley, B., Lodders, K., Marley, M. S., & Fortney, J. J. 2007, Bulletin of the American Astronomical Society, 38, 450
  * [Wallace(1976)] Wallace, L. 1976, Jupiter, 284
  * [Yung & Demore(1999)] Yung, Y. L., & Demore, W. B. 1999, Photochemistry of planetary atmospheres / Yuk L. Yung, William B. DeMore. New York : Oxford University Press, 1999.
  * [Zahnle(2008)] Zahnle, K. 2008, Nature, 454, 41


# Functional data analysis: Application to the second and third wave of
COVID-19 pandemic in Poland

Patrycja Hęćka<EMAIL_ADDRESS>
###### Abstract

In this article we use the methods of functional data analysis to analyze the
number of positive tests, deaths, convalescents, hospitalized and intensive
care people during second and third wave of the COVID-19 pandemic in Poland.
For this purpose firstly we convert the data to smooth functions. Then we use
principal component analysis and multiple function-on-function linear
regression model to analyze waves of COVID-19 pandemic in Polish voivodeships.

###### keywords:

functional data analysis , COVID-19 , functional principal component analysis
, smooth functions , function-on-function regression

###### MSC:

[2020] 62R10 , 62P10

[inst1]organization=Department of Telecommunications and Teleinformatics,
Wrocław University of Science and Technology, addressline=Wybrzeże
Wyspiańskiego 27, city=Wrocław, postcode=50-370, country=Poland

##  1 Introduction

The SARS-CoV-2 virus has become a global problem since it was revealed at the
end of 2019 in China. On March 4 the first case was detected in Poland. In
this paper, we use the methods of functional data analysis to analyze the
number of hospitalized and intensive care people during the second and third
wave of the COVID-19 pandemic in selected, Polish voivodeships.

Determining the boundaries of individual coronavirus waves in Poland is
contractual. The first case of infection was found on March 4, 2020 in a
66-year-old man in hospital in Zielona Góra. The beginning of the first wave
in Poland is therefore assumed to be spring 2020. The second COVID-19 wave in
Poland lasted five months - from September to January. The peak of the second
wave fell in November with a record increase - 27,875 infections - November 7,
2020. The third wave began three months after the peak of the second wave. It
is considered that the beginning of the third wave is February 16, 2021. Her
record was on April 1 with 35,251 new cases of SARS-CoV-2 - the highest daily
number of infections since the beginning of the pandemic in Poland. However,
due to deficiencies in the published data to the public, we assume as the
second wave the period between October 23, 2020 and February 15, 2021, and as
the third wave the time from February 16 to July 5, 2021. The data we used has
been collected and published by Michał Rogalski ([1], contact:
<EMAIL_ADDRESS>see data source: http://bit.ly/covid19-poland,
accessed: March 1, 2022).

Figure (1) shows daily, discrete observations of the number of positive
COVID-19 test results, number of deaths, convalescents, hospitalized people
and people in serious condition.

Figure 1: Daily observations of the number of positive tests results, deaths,
convalescents, hospitalized people and people in critical condition since
October 23, 2020 to July 5, 2021 in all voivodeships in Poland.

In order to reduce the influence of the number of inhabitants of a given
voivodeship on the analysis the number of cases, hospitalized people,
convalescents, deaths, and people in a serious condition we divide by the
number of inhabitants of a given voivodeship. Number of inhabitants we
download from the website of the Statistical Information Centre (see: [2],
”Area and population in the territorial profile in 2021, Area, population
number and density, as of 1 January 2021”, date of publication: July 22, 2021,
accessed: March 7, 2022). Then the obtained numbers were multiplied by
100,000.

From the analysis of appropriately scaled values, it turns out that the
largest number of hospitalized and in serious condition people per 100,000
inhabitants between October 23, 2020 and July 5, 2021 was recorded in the
Świętokrzyskie Voivodeship. Most positive tests and deaths were noted in the
Kujawsko-Pomorskie Voivodeship. The largest number of convalescents was
recorded in the Warmińsko-Mazurskie Voivodeship. The least people were
hospitalized in the Wielkopolskie Voivodeship, people in a serious condition
in the Pomorskie Voivodeship, deaths in the Małopolskie Voivodeship, positive
tests in Podkarpackie, and convalescents in Podlaskie. For this reason for the
test sample of the model which will be analyzed in the next chapters we select
Świętokrzyskie, Wielkopolskie, Podkarpackie and Małopolskie voivodeships.

Since then, the variability over time of positive test results, deaths,
recoveries, hospitalized people and people in a serious condition will be
considered through functional variables respectively:
$X_{1}(t),X_{2}(t),X_{3}(t),Y_{1}(t),Y_{2}(t)$. Variables $X$ will be treated
as predictors, and $Y$ as response variables in functional regression models.
The observed data is the number of daily values of these five functional
variables for sixteen voivodeships in Poland.

Figure (2) shows development of the disease in selected, most diverse
voivodeships by discrete, daily, scaled observations, i.e. a set of curves
$\\{(x_{ij}(t),y_{ik}(t)):i=1,...,16;j=1,2,3;k=1,2\\}$.

Figure 2: Daily observations of the number of positive tests results, deaths,
convalescents, hospitalized people and people in critical condition since
October 23, 2020 to July 5, 2021 in selected voivodeships in Poland. The
number of people was divided by the number of inhabitants of a given
voivodeship and then multiplied by 100,000.

An important limitation on the FDA methods is that the functional variables
should be observed in the same interval. The classic solution to such a
problem is to transfer all curves to the same interval. In this work, we will
consider curves on the interval $T=[0,1]$. Hence, from this points $x_{ij}$
and $y_{ik}$ present the curves located at the interval $[0,1]$.

##  2 From functional data to smooth functions

There are a lot of methods of conversion data to smooth functions (see [3],
[5], [6]). Our goal is to use discrete data $y_{k}$, $k=1,...,m$ to estimate
the function $x$. The main procedure in statistics, mathematics, and
engineering for converting discrete data in smooth functions is the expansion
of the base.

Consider the vectors $x_{i}(t):i=1,...,n;t\in T=[0,1]$ and assume that the
observations $y_{ik}$ are available for each knots
$t_{i1},t_{i2},...,t_{im_{i}}\in T$. Then we can present each observation
$y_{ik}$ as follows:

$y_{ik}=x_{i}(t_{ik})+\epsilon_{ik},i=1,...,n;k=1,...,m_{i},$ (1)

where $\epsilon_{ik}$ is a noise contributes to the roughness of the analyzed
data.

Suppose that the sample curves belong to the finite dimensional space
generated by the set of basis functions $\\{\phi_{1}(t),...,\phi_{p}(t)\\}$.
Then we can represent each curve $x_{i}$ by a linear expansion of the form:

$x_{i}(t)=\sum_{j=1}^{p}\alpha_{ij}\phi_{j}(t),i=1,...,n.$ (2)

We define by $\alpha_{i}$ vectors of the basis coefficients:
$\alpha_{i}=(\alpha_{i1},...,\alpha_{ip})^{\prime}$, which lengths are equal
to p and can be estimated by different methods. One of the most popular method
is least squares method (see [3], section 4.2). The least squares estimators
of $\alpha_{i}$ are of the form
$\hat{\alpha}_{i}=(\Phi_{i}^{\prime}\Phi_{i})^{-1}\Phi_{i}^{\prime}y_{i},$
where $\Phi_{i}=(\phi_{j}(t_{ik}))_{m_{i}\times p},j=1,...,p;k=1,...,m_{i}$.

Smoothness of the function is controlled by the number of used basis
functions. The greater the number of basis functions p, the better the curve
fits to discrete points. The smaller the p, the smoother the curve is.
Decision about increasing or decreasing the number of basis functions is
related to achieving a compromise between bias and variance - ”bias-variance
trade-off” (see [3], section 4.5.1). One of the methods using to get a ”trade-
off” between variance and bias is minimalization of the mean-squared error
(MSE).

Figure (3) shows the fitting of 20 basis functions with equally spaced knots
in order to approximate seven curves showing the number of cases, deaths,
recoveries, hospitalized people and people in a serious condition because of
COVID-19. Coefficients for each functional form were fitted by least squares.
The number of basis functions was chosen so that the value of the mean squared
error was the smallest, but also to avoid overfitting of the model.

Figure 3: Fitted curves to daily observations of the number of positive tests
results, deaths, convalescents, hospitalized people and people in critical
condition since October 23, 2020 to July 5, 2021 in selected voivodeships in
Poland. 20 basis functions were used.

##  3 Multiple Function-on-Function Linear Model

In this section we show multiple function-on-function linear model which was
described in details e.g. by Acal, Escabias, Aguilera and Valderrama (see [4],
section 3.1) and by Xiong Cai, Liugen Xue and Jiguo Cao [7]. Such a model is
often used to characterize pandemic in Europe for example in France [10],
Italy [11] and Spain [4].

The MFFLR model let for estimation the functional response variable $Y$ from a
vector of $J$ functional predictor variables denoted by
$X=(X_{1},...,X_{J})^{\prime}.$ Consider a random sample from $(X,Y)$ denoted
by $(x_{i},y_{i}):i=1,...,n$ with $x_{i}=(x_{i1},x_{i2},...,x_{iJ})^{\prime}$.

Then we define the functional linear model as follows:

$y_{i}(t)=\alpha(t)+\sum_{j=1}^{J}\int_{T}x_{ij}(s)\beta_{j}(s,t)ds+\epsilon_{i}(t),i=1,...,n,$
(3)

where $\alpha(t)$ is the intercept function, $\beta_{j}(s,t)$ are coefficient
functions and $\epsilon_{i}(t)$ are independent functional errors.

Consider the decomposition of the principal components of the functional
response variable and the functional predictior variables, given by

$x_{ij}(t)=\overline{x}_{j}(t)+\sum_{l=1}^{n-1}\xi_{il}^{x_{j}}f_{l}^{x_{j}}(t),$
(4) $y_{i}(t)=\overline{y}(t)+\sum_{l=1}^{n-1}\xi_{il}^{y}f_{l}^{y}(t),$ (5)

where $\xi_{il}^{x_{j}}$ and $\xi_{il}^{y}$ are the principal components
scores. The eigenfunctions of the sample covariance $x_{ij}(t)$ and $y_{i}(t)$
are denoted by $f_{l}^{x_{j}}$ and $f_{l}^{y}$.

The principal components decomposition given by the formula (4) allows for the
transformation of the MFFLR model described by the formula (3) into a linear
regression model for each principal component response variable $Y$ from the
components of the principal functional predictors, given by the formula (6).

$\hat{\xi}_{ik}^{y}=\sum_{j=1}^{J}\sum_{l=1}^{n-1}b_{kl}^{x_{j}}\xi_{il}^{x_{j}}+\epsilon_{ik},i=1,...,n;k=1,...,n-1.$
(6)

The functional coefficients are given here by
$\beta_{j}(s,t)=\sum_{k=1}^{n-1}\sum_{l=1}^{n-1}b_{kl}^{x_{j}}f_{k}^{x_{j}}(s)f_{l}^{y}(t).$

Finally we get the following PC-MFFLR model for the functional response:

$\hat{y}_{i}(s)=\overline{y}(s)+\sum_{k=1}^{K}\hat{\xi}_{ik}^{y}f_{k}^{y}(s)=\overline{y}(s)+\sum_{k=1}^{K}(\sum_{j=1}^{J}\sum_{l\in
L_{kj}}\hat{b}_{kl}^{x_{j}}\xi_{il}^{x_{j}})f_{k}^{y}(s),$ (7)

where K is the number of principal components selected to the model, and
$\hat{b}_{kl}^{x_{j}}$ are linear least-squares estimators of the regression
coefficients $b_{kl}.$

Let’s assume that we have n completely observed curves for all variables and m
missing curves for the response variable. For missing response curves, the
parameters $b_{kl}$ in model (6) are estimated using complete n sample
response curves and predictors. Then the missing response curves
$y_{i}^{miss}(s):i=n+1,...,n+m$ are estimated by computing the principal
components scores of the predictors:
$\xi_{il}^{x_{j}}:i=n+1,...,n+m,l=1,...,n-1$ and inserting them into equation
(7). Then the estimated PC-MFFLR model can be used to predict the value of the
new response variable Y on the test sample.

We can solve the imputation problem by using the multiple function-on-function
linear regression model for each response variable $Y_{1}(t)$ (hospitalized
people) and $Y_{2}(t)$ (people in critical condition). Both functional
regression models are estimated from the full data of 12 voivodeships that
determine the training sample. Then predictions are made for the four,
selected voivodeships: Małopolskie, Wielkopolskie, Świętokrzyskie and
Podkarpackie.

##  4 Data analysis

In this section, we use Principal Component Analysis and a PC-MFFLR model in
order to analyze functional data and predict missing response curves. The
results were obtained with software R (packages ’fda’ [9], ’ggplot2’ [8]).

###  4.1 Principal components analysis for functional data

A lot of papers describe PCA in the functional context (for example see [3],
section 8).

Our first step is to estimate the principal functional components for each of
the five functional predictors. It turns out that the first principal
components explain respectively $46.93\%$, $42.27\%$, $38.78\%$, $41.17\%$,
$44.30\%$ of the variability of $X_{1}$, $X_{2}$, $X_{3}$, $Y_{1}$, $Y_{2}$
from the training sample. The second principal components explain: $26.59\%$,
$19.30\%$, $26.79\%$, $29.41\%$, $31.70\%$ of the variability. The third
principal components explain much less, respectively: $13.96\%$, $14.97\%$,
$18.91\%$, $15.45\%$, $12.71\%$.

####  4.1.1 Weight functions

Figure (4) presents the weight functions (harmonics) related to the first 3
principal components. The presented weight functions are coefficients that
enable the eigenvectors to be computed from the original basis.

Figure 4: Weight functions $f_{l}^{x_{j}};j=1,2,3$ and $f_{l}^{y_{k}};k=1,2$,
$l=1,2,3$ for the first three principal components.

For the number of hospitalized people, the first, second and third principal
components explain together about 86.03, and for the number of people in
critical condition 88.71 percentage of the variance. Other principal
components explain a small percentage of information. Graphs of weight
functions are difficult to interpret. In the next subsubsection we present
plots that can be helpful to analyze the functional principal components.
However, looking at the figure (4), we can gain some intuition about them.

The first principal component shows the general variability of the number of
hospitalized people depending on the wave of the pandemic. We can see that the
biggest difference between the second and the third wave in the number of
hospitalized people occurs at the second half of the analyzed time - the third
wave of the COVID-19 pandemic. We can see that the smallest differences
between voivodeships are during the second wave of the pandemic. This may
suggest that the prediction of the number of hospitalized people during the
third wave of COVID-19 may be the hardest. In between the second and the third
wave we see negative values of the coefficients, which tell us about a
decrease in the number of hospitalized people. Voivodeships for which the
value of $\xi_{i1}^{y_{1}}$ is high will have large differences between the
number of hospitalized people in the second and third wave. During the third
wave, the number of hospitalized people significantly exceed the number of
such people during the second wave of the pandemic. It turns out that the
highest value of this coefficient is achieved for the Łódzkie Voivodeship.

Since the weight function $f_{2}^{y_{k}};k=1,2$ must be orthogonal to
$f_{1}^{y_{k}}$, we cannot expect that the second principal component will
explain a greater percentage of the variance than the first principal
component. In the case of the number of hospitalized people, the second
principal component explain about 29.41 percentage of variability. We see that
it achieves positive values during all of the time of the second and third
waves of the pandemic. The second principal component can be interpreted as an
indicator of the number of people hospitalized during the entire analyzed
period.

The first principal component for the number of people in a critical condition
suggests that the difference in the number of people in a serious condition
between voivodeships was at a similar level for almost the entire second and
third waves of COVID-19. At this time, the values of the coefficients are
positive. The highest value of $\xi_{i1}^{y_{2}}$ is obtained for the
Kujawsko-Pomorskie Voivodeship.

The third and further principal components explain a much smaller proportion
of variance than the first two components. This is influenced by the fact that
they must be orthogonal to the first two principal components and also
orthogonal to each other. They are more difficult to interpret than first two
components.

Interpretation of weight functions is not always simple in the context of
functional PCA. In the subsubsection (4.1.2) we show more commonly used form
of presenting results.

####  4.1.2 Mean curve

A method that is helpful during the analysis of functional principal
components is to present the mean curve together with the functions obtained
by the addition and subtraction of the properly multiplied harmonic (weight)
functions of the principal components from the mean. Such a plot makes sense
because the principal components represent the variation around the mean.
Figures (5) and (6) presents the mean curves and the perturbations of the
sample mean curves obtained by adding and subtracting a multiple of weight
functions for hospitalized and in critical condition people.

Figure 5: Mean curve of hospitalized people with curves resulting from adding
(+) and subtracting (-) appropriately scaled harmonic coefficients from the
mean.

Figure 6: Mean curve of people in critical condition with curves resulting
from adding (+) and subtracting (-) appropriately scaled harmonic coefficients
from the mean.

Analyzing the figure (5), we can see that in terms of the number of
hospitalized people, the first principal component shows the differences
between the second and third wave of the pandemic, while the second principal
component shows the general number of hospitalized people during pandemic.

The figure (6) suggests that for the number of people in a serious condition
the situation was reversed. The first principal component focuses on the
general number of people in critical condition during pandemic, while the
second principal component shows the differences between the waves.

The next components explain a smaller percentage of variance. This can be see
by examining the plots of the fourth harmonic function, where the plots of (+)
and (-) often coincide with the plot of the mean. The smaller the percentage
of variance is explained by the principal component, the more the (+) and (-)
plots coincide with the mean function.

####  4.1.3 Plotting principal component scores

An important aspect of PCA is the examination of the scores of each curve on
each component (see [3], section 8.3.2).

Figures (7) and (8) show other graphs interesting during FPCA analysis. They
represent the values of the variables achieved on the first and second scores.

Figure 7: Values of the first and second scores for the number of hospitalized
people (voivodeships from the training sample).

Suggested by the conclusions from the analysis of the figure (5), on the right
side of the figure (7) there are voivodeships where there is the largest
difference between the second and third wave of the pandemic. On the left side
of the chart, those for which the difference in the number of hospitalized
people is not that big in both waves. At the top of the figure (7), there will
be voivodeships where the number of hospitalized people was the highest during
the second and third wave of the pandemic. Consistently, at the bottom of the
chart, we observe voivodeships with a small number of hospitalized people per
100,000 inhabitants. According to the figure (7), the largest number of
hospitalized people was recorded in Lubuskie, Podlaskie, Lubelskie and Łódzkie
voivodeships. Comparing with appropriately scaled, discrete, real
observations, the largest number of hospitalized people per 100,000
inhabitants was achieved in the following voivodeships: Łódzkie, Lubuskie and
Lubelskie. The voivodeships with the smallest number of hospitalized people,
i.e. Śląskie, Pomorskie and Opolskie, also coincide with the real data.

Hence, we can confirm the conclusion appearing during the analysis of figure
(4), that the first principal component shows the differences between the
second and third wave, and the second principal component is related to the
number of hospitalized people. To conclude in the upper right corner, we have
voivodeships with a large number of hospitalized cases and with large
differences between the second and third wave. Therefore, Lubuskie and Łódzkie
voivodeships had the highest number of hospitalized people per 100,000
inhabitants and additionally in these voivodeships both waves differed
significantly from each other.

We will now perform an analogous analysis for the figure (8), showing the
values of the first and second scores for the number of people in a critical
condition.

Figure 8: Values of the first and second scores for the number of people in
critical condition (voivodeships from the training sample).

In this case, at the top of the figure (8), we can see the voivodeships for
which the difference between the second and third wave of the pandemic is
high, and low at the bottom of the figure. On the left we have voivodeships
for which the number of people in a critical condition was low, and on the
right - high. Here the situation is the opposite compared to the analysis of
the scores for the number of hospitalized people.

According to the data presented on the figure (2), the largest number of
people in a serious condition per 100,000 inhabitants of a given voivodeship
was reached in the Kujawsko-Pomorskie, Lubuskie and Mazowieckie voivodeships,
and the smallest in: Opolskie, Pomorskie and Śląskie. Thus, we can see that
the analysis of figure (8) seems to be performed correctly.

Hence, the Kujawsko-Pomorskie and Lubuskie voivodeships had the largest number
of people in a critical condition, while the Podlaskie and Warmińsko-Mazurskie
voivodeships had the largest differences between the second and third wave of
the pandemic in the number of people in a serious condition. Similar analyzes
can be made for the number of positive test results, recoveries and deaths.

###  4.2 Function-on-Function Model

Let us consider a training sample composed of all voivodeships except
Świętokrzyskie, Małopolskie, Wielkopolskie and Podkarpackie. The listed
voivodeships will be a test sample and we will make predictions for them.

Equation (8) presents the reduction of the linear function-on-function to a
linear model for the first principal components in terms of the first
principal components of each predictor:

$\hat{\xi}^{y_{k}}_{i1}=\gamma_{0}+\xi_{i1}^{x_{1}}\gamma_{1}^{y_{k}}+\xi_{i1}^{x_{2}}\gamma_{2}^{y_{k}}+\xi_{i1}^{x_{3}}\gamma_{3}^{y_{k}}+\epsilon_{i}^{y_{k}},\textrm{
}k=1,2,i=1,...,16,$ (8)

where $\gamma_{0},\gamma_{1}^{y_{k}},\gamma_{2}^{y_{k}},\gamma_{3}^{y_{k}}$
are the appropriate coefficients obtained by fitting the linear model to the
data.

Based on such models, we estimate the first principal components $Y_{1}(t)$
and $Y_{2}(t)$ from the first principal components of the $X_{1}(t)$,
$X_{2}(t)$ and $X_{3}(t)$.

We predict $Y_{1}(t)$ and $Y_{2}(t)$ using the following equation (9):

$\hat{y}_{ik}(t)=\overline{y}_{k}(t)+\hat{\xi}_{i1}^{y_{k}}f_{1}^{y_{k}}(t),\textrm{
}k=1,2,i=1,...,16.$ (9)

In order to test the models on a training sample, we will use the mean squared
error.

Figure (9) shows the graphs of the observed curves, fitted in the section (2)
together with the predicted curves obtained by applying formula (9) for
several voivodeships selected from the training sample. In table (1) you can
see the values of the mean squared error calculated for all voivodeships from
the training sample.

Figure 9: Observed and predicted curves for selected voivodeships from the training sample. Table 1: Values of the mean squared error for $y_{i1}$ (the number of hospitalized people) and $y_{i2}$ (the number of people in a serious condition) for the voivodeships from the training sample. voivodeship |  MSE($y_{i1})$ |  MSE($y_{i2}$)
---|---|---
dolnośląskie |  4.087584 |  0.8424024
kujawsko-pomorskie |  5.801758 |  0.5078765
łódzkie |  8.146502 |  0.7999353
lubelskie |  8.434905 |  1.1796471
lubuskie |  10.79992 |  0.9807714
mazowieckie |  4.569325 |  0.6134904
opolskie |  9.170555 |  0.7959462
podlaskie |  9.612532 |  1.0092404
pomorskie |  7.526031 |  0.9541513
śląskie |  7.451072 |  1.1456427
warmińsko-mazurskie |  3.906237 |  1.0330436
zachodnio-pomorskie |  6.692964 |  0.9105502

In table (1) the lowest and highest MSE values are marked in bold,
respectively. We can see that the values of the mean squared error for the
number of people in a serious condition are much lower than the values of the
mean squared error for the number of hospitalized people. The lowest value of
$\textrm{MSE}(y_{i2})$ was obtained for the Kujawsko-Pomorskie Voivodeship,
and the highest for the Lubelskie Voivodeship. In turn, the value of the mean
squared error for $y_{i1}$ is the lowest for the Warmińsko-Mazurskie
Voivodeship and the highest for the Lubuskie Voivodeship.

Figure (10) presents predictions for voivodeships from the test sample, and
table (2) presents predicted and observed values during the first and last
days of the analyzed period of the pandemic. In order to compare the
predictions with the true observed data, the predictions in table (2) have
been properly scaled and present the number of hospitalized and seriously ill
people for voivodeship on a given day, and not the number of people per
100,000 inhabitants.

Analyzing the figure (10) and table (2) we can see that the predictions for
the Wielkopolskie Voivodeship turned out to be very close to the true values.
The worst predictions we get for Małopolskie Voivodeship. For example, on
November 1, 102 people in a serious condition were recorded in the
Wielkopolskie Voivodeship, and the model predicted 107 people; 217 were
observed in the Małopolskie Voivodeship, while the model predicted 126. On the
same day, 95 people were found in a serious condition in the Podkarpackie
Voivodeship, against the expected 76, and in the Świętokrzyskie Voivodeship,
55 against 41. In turn, the number of hospitalized people on November 1 in the
Małopolskie Voivodeship amounted to 2,366 people, and the model predicted
1,591 cases; in the Podkarpackie Voivodeship 976 against the predicted 977; in
the Świętokrzyskie Voivodeship 829 were observed against 589 predicted, and in
the Wielkopolskie Voivodeship, 1,292 against 1,700.

Analyzing the figure (10) and table (2) we can see that the closest results to
the observed values were predicted at the beginning and at the end of the
analyzed time. The largest difference between the observed and predicted
values is for the Małopolskie Voivodeship. Here the model sometimes predicts
almost 2 times lower values than observed.

Figure 10: Observed and predicted curves for the test sample. Table 2:
Predicted values for the number of hospitalized people and in the serious
condition compared to the observed values. Observations are denoted by ”obs”
and predictions by ”pred”.

Hospitalized people
---
|  małopolskie |  podkarpackie |  świętokrzyskie |  wielkopolskie
time |  obs/pred |  obs/pred |  obs/pred |  obs/pred
23.10 |  1671/1023 |  681/629 |  559/377 |  924/1088
24.10 |  1667/1086 |  750/668 |  581/401 |  929/1158
25.10 |  1763/1151 |  757/707 |  589/426 |  961/1229
26.10 |  1907/1216 |  796/747 |  704/450 |  1039/1300
27.10 |  1972/1281 |  819/786 |  742/474 |  1090/1370
28.10 |  2003/1345 |  813/826 |  750/498 |  1137/1439
29.10 |  2072/1408 |  848/865 |  792/522 |  1182/1507
30.10 |  2250/1471 |  883/903 |  789/545 |  1241/1574
31.10 |  2289/1532 |  917/941 |  779/567 |  1252/1638
1.11 |  2366/1591 |  976/977 |  829/589 |  1292/1700
… |  … |  … |  … |  …
26.06 |  72/106 |  31/66 |  69/39 |  81/111
27.06 |  72/100 |  26/62 |  61/37 |  79/105
28.06 |  66/95 |  25/58 |  57/34 |  82/99
29.06 |  63/89 |  20/55 |  54/32 |  82/93
30.06 |  57/82 |  15/51 |  51/30 |  63/86
1.07 |  53/76 |  13/47 |  47/27 |  42/79
2.07 |  51/69 |  13/43 |  39/25 |  44/71
3.07 |  53/62 |  13/39 |  36/22 |  44/62
4.07 |  48/54 |  15/34 |  32/19 |  44/53
5.07 |  41/46 |  14/29 |  33/15 |  39/43

People in critical condition
---
23.10 |  163/74 |  50/46 |  46/26 |  59/74
24.10 |  160/80 |  54/49 |  46/27 |  56/75
25.10 |  172/85 |  56/52 |  46/29 |  75/77
26.10 |  188/91 |  68/55 |  51/30 |  74/80
27.10 |  207/97 |  67/58 |  52/32 |  77/83
28.10 |  215/103 |  72/62 |  50/34 |  68/87
29.10 |  209/109 |  70/65 |  50/35 |  80/92
30.10 |  208/115 |  81/69 |  53/37 |  94/97
31.10 |  212/121 |  90/72 |  55/39 |  94/102
1.11 |  217/126 |  95/76 |  55/41 |  102/107
… |  … |  … |  … |  …
26.06 |  11/21 |  5/13 |  4/7 |  11/18
27.06 |  11/20 |  4/12 |  3/7 |  10/17
28.06 |  14/19 |  3/11 |  2/6 |  9/16
29.06 |  12/18 |  2/11 |  2/6 |  9/14
30.06 |  11/17 |  2/10 |  2/5 |  9/13
1.07 |  11/16 |  0/9 |  1/5 |  7/11
2.07 |  11/14 |  0/8 |  0/4 |  7/9
3.07 |  11/13 |  0/7 |  0/4 |  7/8
4.07 |  11/11 |  1/6 |  0/3 |  7/6
5.07 |  10/9 |  0/5 |  0/2 |  5/4

##  5 Conlusions

The COVID-19 pandemic has shocked the whole world. An epidemic of this scale
is relatively a new phenomenon, which is why it has attracted the attention of
a large number of analysts around the world. In the presented work, our target
was to fit a model based on functional analysis data to data of the second and
third wave of the COVID-19 pandemic in Poland. The aim of the model was to
predict the number of people in a serious condition and hospitalized people.
The following were used as predictors: the number of deaths, convalescents and
positive test results. Estimations of the parameters were made on the training
group consisting twelve voivodeships. Then we test the quality of these
parameters. Next, predictions were made for the voivodeships: Małopolskie,
Podkarpackie, Świętokrzyskie and Wielkopolskie.

Model predicts numbers well in most voivodeships. The biggest problem turned
out to be Małopolskie Voivodeship. The best predictions we observed for
Wielkopolskie and Świętokrzyskie voivodeships. The analysis of principal
components turned out to be interesting. The figures (7) and (8) showed that
during the second and third wave, most people were hospitalized in the
Lubuskie and Podlaskie voivodeships and the largest number of people in a
serious condition was in voivodeships Kujawsko-Pomorskie and Lubuskie. The
biggest difference between the second and third waves in the number of people
hospitalized was in the Łódzkie and Śląskie voivodeships. The biggest
difference between the second and third wave in the number of people in
serious condition was in Podlaskie and Warmińsko-Mazurskie voivodeships.

## Acknowledgements

The author is grateful to Prof. Krzysztof Topolski (University of Wrocław) for
helpful comments on the master’s thesis.

## Funding

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## Declaration of interest

None.

## References

  * [1] [dataset] The author of the dataset: Michał Rogalski,  COVID-19 w Polsce,  Dane zebrane na podstawie raportów podawanych przez Ministerstwo Zdrowia, danych z WSSE, PSSE, Urzędów Wojewódzkich, oraz tych uzyskanych w prośbach o dostęp do informacji publicznej, contact<EMAIL_ADDRESS><bit.ly/covid19-poland>, accessed: March 1, 2022.
  * [2] [dataset] Statistical Information Centre,  Area and population in the territorial profile in 2021,  Area, population number and density, as of 1 January 2021, date of publication: July 22, 2021, accessed: March 7, 2022, <https://stat.gov.pl/en/topics/population/population/area-and-population-in-the-territorial-profile-in-2021,4,15.html>.
  * [3] J. O. Ramsay, B. W. Silverman,  Functional Data Analysis,  Springer Series in Statistics, Springer, New York, NY, Edition Number 2, Published: 10 November 2010, https://doi.org/10.1007/b98888.
  * [4] Christian Acal, Manuel Escabias, Ana M. Aguilera, Mariano J. Valderrama,  COVID-19 Data Imputation by Multiple Function-on-Function Principal Component Regression,  Mathematics 2021, 9, 1237, Academic Editor: Jin-Ting Zhang, Published: 28 May 2021. https:// doi.org/10.3390/math9111237.
  * [5] Shahid Ullah, Caroline F Finch,  Applications of functional data analysis: A systematic review,  BMC Medical Research Methodology, 2013, <http://www.biomedcentral.com/1471- 2288/13/43>.
  * [6] De Boor, C.,  A practical guide to splines,  2001, revised edition, Berlin: Springer-Verlag.
  * [7] Xiong Cai, Liugen Xue and Jiguo Cao,  Variable Selection for Multiple Function-on-Function Linear Regression,  Statistica Sinica 32 (2022), 1435-1465, doi: https://doi.org/10.5705/ss.202020.0473.
  * [8] Wickham H,  ggplot2: Elegant Graphics for Data Analysis,  Springer-Verlag New York, ISBN 978-3-319-24277-4, <https://ggplot2.tidyverse.org>, R package version 3.3.0, published: 2022-05-03.
  * [9] J. O. Ramsay, Hadley Wickham, Spencer Graves, Giles Hooker,  fda: Functional Data Analysis in R,  R package version 5.1.7, published: 2022-04-26, <https://CRAN.R- project.org/package=fda>.
  * [10] Kayode Oshinubi, Firas Ibrahim, Mustapha Rachdi and Jacques Demongeot,  Functional data analysis: Application to daily observation of COVID-19 prevalence in France,  AIMS Mathematics, 7(4): 5347–5385, Published: 06 January 2022. DOI: 10.3934/math.2022298.
  * [11] Tobia Boschi, Jacopo Di Iorio, Lorenzo Testa, Marzia A. Cremona, Francesca Chiaromonte,  Functional data analysis characterizes the shapes of the first COVID-19 epidemic wave in Italy,  Sci Rep 11, 17054 (2021). https://doi.org/10.1038/s41598-021-95866-y.