# Jupiter Atmospheric Models and Outer Boundary Conditions for Giant Planet
Evolutionary Calculations

Yi-Xian Chen Department of Astrophysical Sciences, Princeton University, 4 Ivy
Lane, Princeton, NJ 08544, USA  Adam Burrows Department of Astrophysical
Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA  Ankan
Sur Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane,
Princeton, NJ 08544, USA  Roberto Tejada Arevalo Department of Astrophysical
Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA

###### Abstract

We present updated atmospheric tables suitable for calculating the post-
formation evolution and cooling of Jupiter and Jupiter-like exoplanets. These
tables are generated using a 1D radiative transfer modeling code that
incorporates the latest opacities and realistic prescriptions for stellar
irradiation and ammonia clouds. To ensure the accuracy of our model
parameters, we calibrate them against the measured temperature structure and
geometric albedo spectrum of Jupiter, its effective temperature, and its
inferred internal temperature. As a test case, we calculate the cooling
history of Jupiter using an adiabatic and homogeneous interior and compare
with extant models now used to evolve Jupiter and the giant planets. We find
that our model reasonably matches Jupiter after evolving a hot-start initial
condition to the present age of the solar system, with a discrepancy in
brightness temperature/radius within two per cent. Our algorithm allows us to
customize for different cloud, irradiation, and metallicity parameters. This
class of boundary conditions can be used to study the evolution of solar-
system giant planets and exoplanets with more complicated interior structures
and non-adiabatic, inhomogeneous internal profiles.

Jupiter, Gas giants, Exoplanets, Planet Atmospheres, Planet Evolution

††software: CoolTLusty (Hubeny & Lanz, 1995; Sudarsky et al., 2000, 2003,
2005; Burrows et al., 2008). Relevant boundary condition entropy tables used
for generating our figures are presented at
https://doi.org/10.5281/zenodo.8297690.

##  1 Introduction

The cooling and contraction of the mostly convective interiors of giant
planets occur after their formation by core accretion (Pollack et al., 1996;
D’Angelo et al., 2003; Li et al., 2021) or disk instability (Boss, 1997), and
is regulated by the radiative properties of the thin layers of atmosphere on
the planets’ surface. To determine the external boundary conditions for time-
dependent evolutionary calculations, it is crucial to generate a set of model
atmospheres at various values of surface gravity and internal emission and/or
entropy. The interior models, which provide information on the thermodynamic
properties and composition of the planet, can therefore be connected over time
using boundary conditions obtained through interpolation from this grid of
models.

Pioneering investigations into the evolution of Jupiter and Saturn, such as
Graboske et al. (1975), Pollack et al. (1977), and Hubbard (1977), were
conducted before the discovery of hot Jupiters in 1995 (Mayor & Queloz, 1995).
Subsequently, new sets of models were developed to cover a wider range of
parameters motivated by the need to study exoplanets (Burrows et al., 1997,
2001), and include better prescriptions for irradiation (Sudarsky et al.,
2000; Burrows et al., 2006a), opacities (Sharp & Burrows, 2007), and cloud
formation (Burrows et al., 2006b). These improvements have resulted in more
refined atmospheric boundary conditions for planetary evolution, including
ones customized for specific extrasolar sources (e.g. Burrows et al., 2003) or
solar system giant planets (Fortney et al., 2011).

Emission spectra of gas giants produced in self-consistent atmosphere
calculations can be directly compared with observed spectra of giant planets
to constrain their properties (Sudarsky et al., 2003). The launch of JWST has
ushered in a new era of exoplanet atmospheric detection with unprecedented
precision (Carter et al., 2022; Miles et al., 2023), requiring more complete
archives of reference theoretical spectra. In this sense, modeling the outer
atmospheres of gas giants and generating tabulated spectra are critical for
two main reasons. First, theoretical models are necessary to interpret
important physical quantities, such as temperature, gravity, radius,
composition, and albedo, from spectra obtained through direct detection.
Second, theoretical spectral models inform observers concerning which regions
of the spectrum are most likely to yield the most insight into a planet’s
atmosphere. Moreover, spectrum calculations provide more well-calibrated
atmospheric boundary conditions, which can be combined with updated knowledge
of opacity and interior equations for evolutionary models that constrain the
age of directly imaged planets.

This paper is structured as follows: In §2, we describe the numerical methods
used to calculate our models, focusing on additions and alterations, such as
ammonia clouds and irradiation, and compare to previous models for gas giant
evolution (Burrows et al., 1997; Fortney et al., 2011). In §3, we calibrate a
fiducial set of numerical parameters, including irradiation redistribution and
cloud modal size, based not only on Jupiter’s metallicity ($\sim$3$\times$
solar), gravity, and measured effective temperature, but also on observed
albedo and temperature-pressure profiles. In §4, we calculate the cooling of
gas giants using a standard interior equation of state, and illustrate the
effects of cloud and irradiation parameters on the evolution of gas-giant
radii, the internal temperature ($T_{\rm int}$), and the effective temperature
($T_{\rm eff}$). Finally, in §5, we discuss our conclusions and their
implications for the next generation of giant-planet evolutionary models that
incorporate more sophisticated thermal and compositional profiles and energy
transport modalities.

##  2 Model Atmospheres

We generate model atmospheres using the 1D atmosphere and spectral code
CoolTLusty (Hubeny & Lanz, 1995; Sudarsky et al., 2003, 2005; Burrows et al.,
2008). This code self-consistently calculates the spectrum and structure of a
plane-parallel atmosphere with radiative transfer, given a detailed suite of
opacity data compiled into thermochemical equilibrium tables (Burrows & Sharp,
1999; Sharp & Burrows, 2007). The opacities are treated using line-by-line
sampling and have recently been significantly updated in Lacy & Burrows (2023)
(see their Appendix A). Their strategy for calculating absorption cross
sections is mainly guided by recent progress by the ExomolOP (Chubb et al.,
2021) and EXOPLINES (Gharib-Nezhad et al., 2021) collaborations. A publicly
available Fortran code, exocross
111https://exocross.readthedocs.io/en/latest/, was employed to compute
absorption cross sections over a grid of temperatures, pressures, and
wavenumbers. User-defined aspects of exocross absorption cross-section
calculations include line lists, line profiles with pressure-dependent
broadening if desired (the Voigt profile is usually applied), line-wing
cutoffs, and an optional line strength threshold. For each molecule, the line
lists recommended by the ExoMol team were adopted, with additional updates to
cover higher temperature and larger wavelength ranges. The most relevant
updated molecules include water, ammonia, methane, and molecular hydrogen. The
line list for molecular ${\rm CH}_{4}$ (Yurchenko & Tennyson, 2014) is
additionally supplemented at shorter wavelengths by cross sections inferred
from Jupiter’s spectrum in Karkoschka (1994).

In our atmospheric calculations, we use 100 atmospheric layers, and 5000
frequency points spaced evenly in log from 0.5 to 300 $\mu$m. We neglect
disequilibrium chemistry discussed in Lacy & Burrows (2023), but add
treatments for stellar irradiation (§2.2). For fiducial evolutionary
calculations of a gas giant in isolation, the metallicity ($Z$), the intrinsic
temperature ($T_{\rm int}$; or equivalently the internal flux) and outer
planet radius are needed to construct the thermal boundary condition and
calculate the specific entropy at the base of the atmosphere’s radiative zone.
Inverting this table, $T_{\rm int}$ becomes a function of the surface gravity
and entropy per baryon ($S$) in the deep interior, sometimes parameterized by
$T_{10}$, the temperature that the interior isentrope would have if
extrapolated to a pressure of 10 bars (Burrows et al., 1997; Hubbard et al.,
1999; Fortney et al., 2011). When there is no external stellar irradiation
source, the total effective temperature ($T_{\rm eff}$) is equal to $T_{\rm
int}$, which is why in studies of isolated gas giant evolution $T_{\rm eff}$
is often simply used in the place of $T_{\rm int}$ (e.g. Burrows et al.,
1997).

When stellar irradiation is taken into account, $T_{\rm eff}$ will differ from
$T_{\rm int}$ due to the fraction of absorbed stellar irradiation by the
planet. Though we do not use such a quantity, this fraction is often, though
crudely, identified with a $T_{\rm eq}$, such that $T^{4}_{\rm eff}$ is set
equal to $T^{4}_{\rm eq}$ \+ $T^{4}_{\rm int}$ (Saumon et al., 1996; Sudarsky
et al., 2000). In CoolTLusty, the wavelength-dependent geometric albedos are
calculated self-consistently and the absorbed stellar heat is intrinsically
accounted for (Sudarsky et al., 2003; Burrows et al., 2006a). Although only
$T_{\rm int}$ is relevant to the internal evolution, $T_{\rm eff}$ can be
obtained from the total emission and compared to measured values for Jupiter.
Hence, when we calculate evolutionary tracks for $T_{\rm int}$, we
simultaneously obtain the associated $T_{\rm eff}$ evolution, along with the
frequency-dependent spectra.

###  2.1 Cloud Models

Since we focus on gas giants with effective temperatures of $80$ to $200$ K,
ammonia clouds will form in the atmosphere and have an impact on their late-
time evolution. Such clouds are a much-studied feature of both Jupiter and
Saturn (Brooke et al., 1998; Sromovsky & Fry, 2018; de Pater et al., 2019).
Water clouds appear earlier when the atmospheric temperatures go below
$\sim$400 K, but are quickly buried deep early on. Young, more massive, giant
planets should evince water clouds that for them would need to be included in
evolutionary calculations. However, such is not the case for the current
Jupiter and Saturn and we ignore them here. Our treatment of cloud opacity is
the same as in Lacy & Burrows (2023), and we use the same cloud shape
parameters (the compact “E”-type) and supersaturation factor (0.01). We adopt
the Clausius-Clapeyron line for the base of the cloud. The cloud spatial
distributions we employ are elaborated in earlier works that published cloudy
atmospheric models using CoolTLUSTY (Burrows et al., 2006b; Madhusudhan et
al., 2011). We assume the cloud species to have a Deirmendjian size
distribution (Deirmendjian, 1964):

$n(a)\propto\left(\frac{a}{a_{0}}\right)^{6}\exp\left[-6\left(\frac{a}{a_{0}}\right)\right]\,,$
(1)

with a default modal size ($a_{0}$) on the order of microns. A Deirmendjian
particle size distribution reproduces that of the Earth’s clouds for a $a_{0}$
of 4$\mu$m. As we show below (§3), $a_{0}\sim 1\mu$m is consistent with the
Bond and geometric albedos of Jupiter. However, we also experiment with larger
particle sizes to explore the model dependencies. For a given particle size
distribution, the frequency-dependent absorption cross section
$\sigma(\nu,a_{0})$ is then calculated using Mie theory (Kerker, 1969) and
converted into an opacity per unit mass (Lacy & Burrows, 2023).

###  2.2 Irradiation

We treat irradiation by a sun-like star using the approach found in Burrows et
al. (2006a) and Burrows et al. (2008). The solar power contribution at 5.2AU
is calculated by intercepting the incident solar flux with an area of $\pi
R_{J}^{2}$, multiplying by $f=1/2$ to account for an average zenith angle of
$60^{\circ}$ (Appleby, 1986; Marley & McKay, 1999; Fortney et al., 2011). Note
that though Burrows et al. (2008, see their Appendix D) proved that one should
apply a factor of 2/3 in the limit of strongly irradiated Hot Jupiters, a
factor of 1/2 may be more appropriate in the case of gas giants with moderate
to large orbital distances. Therefore, the frequency-integrated stellar flux
expressed via the $H$-moment is

$H_{\rm ext}=f\left(\frac{R_{*}}{d}\right)^{2}\frac{\sigma}{4\pi}T_{*}^{4},$
(2)

where $R_{*}$ is stellar radius, $d$ is the star-planet distance, and $T_{*}$
is the effective temperature of the stellar surface.

Furthermore, we apply a redistribution parameter $P_{\rm irr}$ to represent
the fraction of the flux that is redistributed. Effectively, our treatment
removes a fraction of the irradiation $H_{\rm irr}=P_{\rm irr}H_{\rm ext}$
from the day side (and constitutes an energy input to the night side).
Although all profiles we calculate in this paper belong to the day side, it’s
worth noting that conventional setups with no redistribution ($P_{\rm irr}=0$)
imply zero redistributed heating to the night side.

To redistribute heat from the day side, we add the additional sink term $-D$
to the radiative transfer equation, expressed in terms of the integrated
column mass of the atmosphere $m$:

$D(m)=\frac{2H_{\rm irr}}{m_{1}-m_{0}}\frac{m_{1}-m}{m_{1}-m_{0}}$ (3)

where the parameters $m_{1}$ and $m_{0}$ are column masses corresponding to
the limiting pressures $P_{0}$ and $P_{1}$ (the sink term is set to zero
outside this mass/pressure range), such that $D(m)$ linearly decreases with
$m$, achieving a value of zero at the bottom of the redistribution zone
($P_{1}$). For more details, see Appendices A & B of Burrows et al. (2008).
Integrating $D(m)$ over the column mass confined between the two limiting
pressures gives $H_{\rm irr}$. We restrict the redistribution altitudes to be
between $P_{0}=0.05$ and $P_{1}=0.5$ bar as limiting pressures, guided by the
temperature structure of Jupiter.

As we show below, larger $P_{\rm irr}$ leads to a steeper temperature
inversion in the stratosphere, and we select a fiducial $P_{\rm irr}$ of
$0.15$ based on calibration with measured Jupiter temperature-pressure
profiles.

##  3 Calibration with Jupiter

###  3.1 Fitting the Atmospheric Temperature and Geometric Albedo Profiles

Here, we determine a set of fiducial parameters for the incorporation of the
effects of clouds and irradiation into our atmosphere boundary conditions by
constraining our models to fit observation data, specifically the measured
albedo and temperature-pressure (TP) profiles of Jupiter (Seiff et al., 1998;
Karkoschka, 1998). For this purpose, we fix $\log(g[{\rm cgs}])=3.4$ and the
metal abundance ($Z$) at 3.16 times solar (Niemann et al., 1998). There are
some uncertainties in the helium abundance measurements and the fraction may
also vary with time in more modern evolutionary calculations, but we choose
$Y=0.25$ as the fiducial value and have tested that all outcomes are quite
insensitive to atmospheric $Y$ values between 0.22 and 0.28. In models with
irradiation, we assume the central stellar source of irradiation has a
blackbody spectrum with a temperature of $T_{*}=5777$K, radius of $R_{\odot}$,
and is a distance of $d=5.2$AU from the planet. Furthermore, as has been
traditional in the literature, we employ the “pseudo” effective temperature
$\tilde{T}_{\rm eff}$ as a constraint to help determine consistent values of
$T_{\rm int}$ for each set of cloud or irradiation parameters. More
specifically, we use as a fiducial $\tilde{T}_{\rm eff}=125.57\pm 0.07$K,
based on Cassini CIRS and VIMS observations (Li et al., 2012). These authors
found that while infrared emissions measured using the CIRS (Composite
Infrared Spectrometer, $>7\mu$m) account for most of Jupiter’s emissions, the
emissions around 5$\mu$m measured in the VIMS band (Visual and Infrared
Mapping Spectrometer, 4.5-5.5 $\mu$m) have a non-negligible $\sim 1\%$
contribution to the total emission flux, whose specific fraction varies with
latitude. To be consistent in matching this constraint, our $\tilde{T}_{\rm
eff}$ is defined as the sum of the thermal emission larger than 4.5 microns,
to differentiate with ${T}_{\rm eff}$, the total integrated brightness
temperature in theoretical models. Generally, $\tilde{T}_{\rm eff}$ varies
within a factor of 2% from ${T}_{\rm eff}$ and does not constitute a major
difference from previous evolution models of irradiated Jupiters. Moreover, as
we will elaborate below, only the tabulated values of convective zone entropy
$S$ will determine the evolution track, and not $\tilde{T}_{\rm eff}$. For the
purpose of generating an atmospheric boundary condition calibrated on
measurements of Jupiter, we explore the parameter space of $P_{\rm irr}$ and
$a_{0}$ of models with $\tilde{T}_{\rm eff}\approx 125.57$K, and converge (see
below) on a best-fit model with $P_{\rm irr}=0.15$ and $a_{0}=1\mu$m.

During this process, the Temperature-Pressure(TP) profile measured by the
Galileo entry probe (Seiff et al., 1998), which boasts a characteristic
temperature inversion, places a tight constraint on $P_{\rm irr}$. This
temperature inversion in the “stratosphere” of Jupiter is also seen in Voyager
data (Lindal et al., 1981), and attributed to the interaction between alkanes
and stellar irradiation (Yelle et al., 2001). The lower stratosphere, where
the temperature profile is relatively smooth, is mainly heated by absorption
of methane of sunlight in the near-IR wavelengths. In our modeling, the
redistribution introduces extra effective cooling between $0.05$ and $0.5$
bars. In Figure 1, we plot the TP profiles of characteristic atmospheric
models against the Galileo entry probe data. For comparison, the blue line
corresponds to a model with neither a cloud nor irradiation. In this model,
$T_{\rm int}=T_{\rm eff}\approx\tilde{T}_{\rm eff}$ and no temperature
inversion is observed. The orange line shows that the irradiation-only model
with $P_{\rm irr}=0.15$ fits reasonably well with the entry probe
measurements, and we found that $T_{\rm int}$ needed to be $\sim$97 K to
satisfy the $\tilde{T}_{\rm eff}=125$ K constraint. This is slightly different
from the 99 K estimation of Fortney et al. (2011), which may be a consequence
of their neglect of redistribution. In Figure 1, we include three additional
models with ammonia clouds of modal particle size $a_{0}=1\mu$m and different
values of $P_{\rm irr}$. Generally, at fixed $\tilde{T}_{\rm eff}$, cloud
existence or particle size does not have a strong impact on the TP profile
since the stratosphere is not heated by clouds, as can be seen by comparing
$P_{\rm irr}=0.15$ models with and without ammonia clouds. However, the
temperature profile transition at the “tropopause” is quite sensitive to
$P_{\rm irr}$. While for $P_{\rm irr}=0.3$ the thermal inversion at $\sim 0.5$
bar is too steep, it becomes too smooth without any redistribution cooling,
and for $P_{\rm irr}=0.15$ the TP profile fits well for all cloud particle
sizes. Moreover, a value for $T(1{\rm bar})$ of $166\pm 1$K, consistent with
Galileo entry probe data (Seiff et al., 1998), is well reproduced and comports
with other previous theoretical models of Jupiter’s atmospheric thermal
profile (Guillot et al., 2004). Note that our models still neglect other
potentially important heating/cooling sources such as aerosols Zhang et al.
(2015), wave-breaking (Young et al., 2005), etc., but we emphasize that we
mainly aim to improve upon past published literature interested in boundary
conditions for evolutionary models.

To break the degeneracy in the selection of cloud modal particle size, the
geometric albedo spectrum serves as an additional useful constraint. Inspired
by the conclusion that 1-10 micron ${\rm NH}_{3}$ ice particles fit the
Infrared Space Observatory data (Brooke et al., 1998), we experimented with
clouds with modal particle size ($a_{0}$) of 1, 3, and 10 microns. We note
there are more recent studies, based on New Horizons LEISA (Linear Etalon
Imaging Spectral Array) data, arguing that the characteristic particle size is
more likely between $\sim$2 and 4 microns (Sromovsky & Fry, 2018), with a
relatively wide distribution, and also that ${\rm NH}_{4}{\rm SH}$ solids can
contribute to the clouds.

Exploring this relevant range of parameters, we compare our geometric albedo
spectrum with the high resolution sub-micron geometric albedo spectrum from
the European Southern Observatory (Karkoschka, 1998). Our results for
different models and particle sizes are shown in Figure 2, plotted against the
data of Karkoschka (1998). We note that for the 1-micron case our geometric
albedo is consistent with observations for $\gtrsim 0.5\mu m$. However, the
irradiation-only model (red solid line) seriously underestimates the geometric
albedo profile; even the overall shape cannot be matched. Note that we do not
attempt to model the extra scatterer/chromophores (Carlson et al., 2016) below
$0.5\mu$m wavelength, about which there are significant uncertainties (Lacy et
al., 2019). For cloud sizes of 3 and 10 microns, although the change in
particle size does not affect the temperature profile, the geometric albedo
spectrum itself is slightly below the observations.

In conclusion, we find that the general geometric albedo spectrum can be
reproduced only with the inclusion of clouds and that its values sensitively
increase with decreasing cloud particle size. We also find that for
$a_{0}=1\mu$m the albedo spectrum matches well the measured values in
Karkoschka (1998). For this set of parameters, the internal temperature
($T_{\rm int}$) at a $\tilde{T}_{\rm eff}$ of $\sim 125.57$ K is constrained
to be $\approx$102 K, in contrast with a $T_{\rm int}$ of $\approx 97$ K for
the no-cloud case.

###  3.2 A Discussion Concerning the Bond Albedo

In theory, one can also relate the scattered/reflected fraction of solar
irradiation with an effective average Bond albedo $A_{B}$, where the fraction
that contributes to $\sigma T_{\rm eff}^{4}$ is $1-A_{B}$. Li et al. (2018)
analyzed the Cassini VIMS and ISS (Imaging Science Subsystem) data to estimate
a Bond albedo of 0.503$\pm$0.012, which is significantly larger than the
previous adopted values around 0.343 (Conrath et al., 1989; Fortney & Hubbard,
2003). By applying the incident solar flux at Jupiter’s orbital radius, they
derive an internal flux of 7.485 $\pm$ 0.163 W/m2, consistent with Jupiter’s
total emitted power of 14.098$\pm$0.031 W/m2 (from the earlier measurement of
$\tilde{T}_{\rm eff}=125.57\pm 0.07$K in Li et al. (2012)) and corresponding
to a $T_{\rm int}$ of $\approx 107$ K. As we noted, our cloudless model with
or without redistribution, as well as the cloudless Jupiter model of Fortney
et al. (2011), constrained by $\tilde{T}_{\rm eff}\approx 125$ K, both give
$T_{\rm int}\lesssim 100$K, which is inconsistent with the latest measurements
of Jupiter’s heat balance. As a matter of fact, this discrepancy implies there
should be significant cloud effects at work scattering away a larger fraction
of incident irradiation.

However, we caution that given our planar calculations we have no associated
phase integrals to directly provide the Bond albedo from the geometric albedo
spectrum. The spectrum-integrated Bond albedo can be calculated from a
spectroscopic and atmospheric model by integrating the monochromatic spherical
albedo, weighted by the stellar flux $F_{*}(\lambda)$:

$A_{B}=\frac{\int F_{*}(\lambda)A_{s}(\lambda)d\lambda}{\int
F_{*}(\lambda)}\,,$ (4)

where $A_{s}(\lambda)$ is the spherical albedo (monochromatic Bond albedo).
The spherical albedo and monochromatic Bond albedo are the same and given by
the product of $A_{G}(\lambda)$ and the phase integral ($q(\lambda)$). Due to
the stellar flux weighting, the optical wavelengths dominate the total Bond
albedo. However, $q(\lambda)$ cannot be self-consistently determined from a 1D
planar model, and the scattering of light over all phase angles must be
studied in 2D (Marley & McKay, 1999; Cahoy et al., 2010; Madhusudhan &
Burrows, 2012).

Li et al. (2018) measured the realistic $q(\lambda)$ for Jupiter, averaged
over all phase angles, to be $\approx 1.3$ in the optical wavelength range of
interest (see their Figure 3). This means that in realistic multidimensional
models, we should expect a smaller fraction of the irradiation to contribute
to the planet emissions calculated in 1D; hence, it’s natural that the
internal $T_{\rm int}$ is larger in their measurements than in our fiducial
model (which boasts a consistent geometric albedo spectrum but lacks the phase
integral). We note that, informed by their model, if we multiply our $A_{G}$
by their factor of 1.3 to mimic a spherical albedo, we obtain a consistent
frequency-integrated value for the Bond albedo of 0.5. However, since the
phase integral is measured only for Jupiter and subject to many uncertainties,
we make no attempt to reconcile this mismatch with artificial treatments that
could lead to significant confusion 222However, see the next section for tests
varying $f$ that might inform the general dependence of $T_{\rm int}$ on the
phase integral.. The scattering phase function can be treated realistically
only with multi-dimensional radiative transfer simulations. These
uncertainties and caveats not withstanding, we conclude that our choice of
fiducial parameters: $P_{\rm irr}=0.15$ and $a_{0}=1\mu$m (with both
irradiation and ammonia clouds included) suitably reproduces the observation
data and is sufficient to provide practical atmosphere boundary conditions for
the evolution of Jupiter-like planets. Since redistribution and clouds seem
necessary to reproduce temperature inversions and the geometric albedo
spectrum of Jupiter, the inclusion of these effects is an improvement over
earlier realizations of atmospheric boundary conditions for evolutionary
calculations (Burrows et al., 1997; Fortney et al., 2011).

Figure 1: Comparison of the temperature-pressure profiles of different Jupiter
models at $\tilde{T}_{\rm eff}=125$K. Only the isolated (no irradiation or
cloud) models show no temperature inversion. All other models are irradiated
at 5.2 A.U., either without clouds or with clouds of modal particle size of
$1\mu$m. We observe that $P_{\rm irr}$ is closely related to the temperature
valley at the inversion around 0.1 bar, and that $P_{\rm irr}=0.15$ best
reproduces the observational data. Once the redistribution factor is set, the
existence of clouds does not alter the temperature-pressure profile. Other
particle sizes also have TP profiles that nearly coincide with the green line.
The Galileo entry probe data from Seiff et al. (1998) is plotted in black.
Figure 2: Geometric albedo spectra for different irradiated atmospheric models
without clouds (red) or with clouds of different modal particle sizes, fixing
$P_{\rm irr}=0.15$. The black spectrum is from Karkoschka (1998).  Figure 3:
The specific entropy in the convective zone under the atmosphere’s thin
radiative layer, plotted as a function of gravity, and effective temperature
$\tilde{T}_{\rm eff}$ (red) or internal temperature $T_{\rm int}$ (blue) for
our fiducial model ($f=0.5,P_{\rm irr}=0.15,a_{0}=1\mu$m). At low
temperatures, $\tilde{T}_{\rm eff}$ starts to become significantly larger than
$T_{\rm int}$ due to the increasing contribution of stellar irradiation.

To demonstrate its structure, in Figure 3 we plot the atmospheric entropy per
baryon at depth as a function of either $T_{\rm int}$ or $\tilde{T}_{\rm eff}$
and surface gravity using our fiducial boundary table. The entropy surface in
the $T_{\rm int},g$ plane is shown in blue (with wireframe) and the entropy
surface in the $\tilde{T}_{\rm eff},g$ plane is shown in red. With a planetary
evolution code one follows a gas giant along entropy surfaces from high
entropy and low $g$ (upper left) to low entropy and high $g$ (lower right),
interpolating the table within grid points. It’s apparent that at high entropy
and temperature the surfaces converge, while at low temperature
$\tilde{T}_{\rm eff}$ starts to deviate from $T_{\rm int}$ due to the
increasing contribution of irradiation. In the remainder of this paper, we
will omit the tilde symbol in $\tilde{T}_{\rm eff}$ for simplicity.

##  4 Homogeneous Adiabatic Evolutionary Calculations

In traditional adiabatic cooling calculations, the internal flux and entropy
are related by energy conservation:

$\frac{dL}{dm}=-T\frac{dS}{dt},$ (5)

where $L$ is the luminosity, $S$ is the specific entropy per mass, and $t$ is
the time. At a given time, $S$ would be the same throughout an adiabatic
envelope and independent of radius. $m$ is the integrated mass ($\int^{r}4\pi
r^{\prime 2}\rho(r^{\prime})dr^{\prime}$) and $T$ is the local temperature,
both functions of shell radius $r$. From this equation we can relate the
timestep between two models $\Delta t$ that differ in entropy by $\Delta S$.
The luminosity ($L$) is derived physically from the planet interior and is
equal to $4\pi R^{2}\sigma T^{4}_{\rm int}$, where the value of $T_{\rm int}$
at a given $S$ and surface gravity is inverted from a model atmosphere table
(see Figure 3). $T_{\rm eff}$ at a given time is also inverted from the model
atmosphere table, but does not go into the calculation of the variation of
internal entropy and the evolution of the temperature profile with time.

Figure 4: Left panel: Evolution of internal brightness temperature $T_{\rm
int}$ and emission temperature $T_{\rm eff}$ for our fiducial Jupiter model
(blue) evolved in our adiabatic code, the Fortney et al. (2011) boundary
condition (taken from their Figure 7), and the Burrows et al. (1997) boundary
condition evolved in our adiabatic code. Right panel: the evolution of the
model Jupiter radius. An extra evolutionary track for zero core mass is
plotted with a dashed line to show that a solid core is needed to make the
planet more compact.

In Figure 4, we plot the evolutionary tracks of a “hot-start” Jupiter for our
fiducial boundary conditions, as well as those from the Burrows et al. (1997)
and Fortney et al. (2011) models. All models used the Saumon et al. (1995)
equation of state. The left panel shows $T_{\rm int},\tilde{T}_{\rm eff}$ and
the right panel shows the radius evolution. The Fortney et al. (2011)
evolution track for Jupiter is taken from their Figure 7, for which they
included a 10$M_{\odot}$ metal core and approximately 15$M_{\odot}$ of extra
metallicity in the form of water in the H/He adiabatic envelope. In their
paper, they also used their adiabatic evolution code to calculate the
evolution of Jupiter with Burrows et al. (1997) tables that are for
$Z=Z_{\odot}$ and do not include irradiation. They obtain a low $T_{\rm int}$
of $\approx 100$K at 4.56 Gyrs (see their Figure 7) 333Fortney et al. (2011)
created effective temperatures for the Burrows et al. (1997) boundary
condition table through post-processing of irradiation, assuming a constant
Bond albedo of 0.343. To avoid confusion, we do not adopt this tactic. Hence,
no corresponding $\tilde{T}_{\rm eff}$s are plotted for the Burrows et al.
(1997) calculation.. However, we are unable to reproduce this result from
Fortney et al. (2011) using the Burrows et al. (1997) boundary condition table
at our disposal for any equivalent core/heavy-element mass between $0$ to
$30M_{\oplus}$. Instead, in Figure 4, we show an evolution track produced by
our code with $25M_{\oplus}$. This model yields at 4.567 Gyrs a $T_{\rm int}$
$110$ K, which is larger than the $T_{\rm int}$ from Fortney et al. (2011)
models, but consistent with earlier calculations performed by CoolTLusty.

Nevertheless, with our fiducial table, and a core mass of 25$M_{\oplus}$, we
obtain $\tilde{T}_{\rm eff}\sim 127K$ at 4.56 Gyrs. While $T_{\rm int}$ in our
calculation is smaller than found using the Burrows et al. (1997) table, the
total emission or effective temperature is reasonably close to the Li et al.
(2018) estimation. Moreover, adding the fiducial core mass has a pronounced
influence on the radius evolution. In the right panel of Figure 4, we also
compare the evolution of planetary radius with (solid line) and without
(dashed line) the core mass for the fiducial boundary condition table.
Generically, with a metal core or equivalent heavy-element mass in the
envelope, the gas giant has a larger mean density and is more compact.
Hereafter, we present all results with a solid core of 25 $M_{\oplus}$.

###  4.1 Dependence on the f factor

Figure 5: Left panel: Jupiter’s temperature evolution comparing different
zenith angle factors $f$; Right panel: Jupiter’s temperature evolution
comparing different cloud particle size parameters $a_{0}$.

In the fiducial models, $f=1/2$ is adopted to approximate the average zenith
angle for a moderately irradiated planet. For strong irradiation, this factor
approaches $2/3$ (Burrows et al., 2008) so there is a degree of freedom here
not constrained by planar models. In the left panel of Figure 5, we compare
$T_{\rm int}$ and $T_{\rm eff}$ evolutionary tracks using $f=0.67$ and our
fiducial tables (keeping $a_{0}=1\mu$m). We find that the $T_{\rm int}$ and
$T_{\rm eff}$ values at the current age of Jupiter increase and decrease,
respectively, with decreasing $f$, such that low $f$ evolutionary curves are
generally anchored between those of high $f$, although there are some sharp
non-linear transitions at $T_{\rm int}\sim$150-200K due to the onset of
ammonia clouds. This is consistent with the expectation that these
temperatures should converge towards the isolated planet case with $T_{\rm
eff}=T_{\rm int}$ (similar to the Burrows et al. (1997) evolutionary curve)
when $f$ approaches zero (i.e., when irradiation is ignored). Interestingly, a
decrease in $f$ may also generally represent the effect of an extra scatterer
(§3.2). If the fact that $A_{B}>A_{G}$ is a consequence of $q>1$, then one can
decrease $f$ by a factor of $\approx(A_{B}-A_{G})/(1-A_{G})$ to account for
this reduction. This will result in a larger $T_{\rm int}$ at 4.56 Gyrs, in
the direction of being slightly more consistent with Li et al. (2018).
However, the inclusion of a “phase-integral correction” in $f$ is quite
artificial, and we certainly do not expect $q$ or this modification of the
scattering fraction to be constant in time, Nevertheless, this crude $f$
pseudo-dependence illuminates to zeroth-order the dependence of $T_{\rm eff}$
and $T_{\rm int}$ on the phase integral, such that for multi-dimensional
radiative transfer codes that do include this effect, we expect evolutionary
tracks to move slightly closer to the observations.

###  4.2 Dependence on Particle Size

In the right panel of Figure 5, we display the evolution of $T_{\rm int}$ and
$T_{\rm eff}$ for fiducial and $a_{0}=3\mu$m models (fixing $f=1/2$). Just as
with the dependence on smaller $f$, for smaller particle size, $T_{\rm int}$
at 4.56 Gyrs is larger, but $T_{\rm eff}$ is not necessarily also larger; this
results from the smaller contribution from irradiation due to larger albedo.
For larger particle sizes, the two temperatures tend to diverge, and the
contrast between $T_{\rm int}$ and $T_{\rm eff}$ becomes large. In the small-
particle-size limit, for a negligible irradiation contribution (when incident
radiation is completely scattered), we expect these temperatures to converge,
similar to what is found in the $f\rightarrow 0$ limit.

###  4.3 Dependence on Abundances

All of the models above are for $Y=0.25$ and $Z=3.16Z_{\odot}$ in the
atmosphere. We also compare evolutionary tracks with $Y=0.25,Z=10Z_{\odot}$
and $Y=0.22,Z=3.16Z_{\odot}$ and find that varying $Y$ in the atmosphere has
little or no effect on the planetary evolution. However, an increase in $Z$
leads to an increase in $T_{\rm int}$ and $T_{\rm eff}$ and slows planetary
cooling by raising the atmospheric opacity. This effect arises even at earlier
times, unlike the effect of ammonia clouds which appears only at late stages.
Note again that these abundance variations are in the atmosphere, and that our
internal adiabats still have $Y=0.27,Z=3.16Z_{\odot}$ fixed. The dependence of
gas-giant evolution on internal compositions and their profiles (e.g. effect
of helium depletion) is properly a subject of future work.

##  5 Conclusion

In this study, we develop a set of atmospheric models using the 1D atmosphere
and spectral code CoolTLusty. Our goal is to create state-of-the-art boundary
condition tables that could be used for studying the evolution of gas giant
planets. The atmospheric opacities employed were significantly updated by Lacy
& Burrows (2023) and we implement realistic treatments of clouds and
irradiation to calibrate our parameters with the observed temperature
structure and albedo spectrum of Jupiter. We simulated the internal evolution
of a “hot-start” Jupiter using this set of tables in the context of the
traditional adiabatic paradigm. This approach is being challenged by the new
Juno data (Wahl et al., 2017; Bolton et al., 2017), but the viability of these
atmospheric boundary conditions for any evolutionary model is not compromised.
We find that with reasonable irradiation and cloud parameters we obtain an
atmospheric boundary condition table that cools down a “hot-start” Jupiter to
close to its current measured thermal state with its measured geometric
albedos.

In addition to providing a useful atmosphere reference model for future giant-
planet cooling calculations, we explored the dependence on various physical
parameters. The average zenith angle parameter $f$ and cloud modal particle
size $a_{0}$ have particular sway over the influence of irradiation. If the
absorption fraction of irradiation is very high, then the effective
temperature $T_{\rm eff}$ may not be able to cool down to near the measured
value of $\sim$125 K at Jupiter’s current age. However, by using fiducial sets
of cloud parameters to ensure a reasonable geometric albedo spectrum, we
observed that $T_{\rm eff}$ cools down faster than observed using previous
approaches (Burrows et al., 1997; Fortney et al., 2011), better matching the
surface observations of Jupiter. The average zenith angle also has an impact,
and for a slight reduction in $f$, $T_{\rm eff}$ also cools down a bit faster.
In addition, the atmospheric metallicity affects the cooling rate, with a
higher $Z$ raising the atmosphere opacity and resulting in slower cooling. The
helium fraction $Y$ in the atmosphere has a minimal effect. For this study, we
fixed the stellar luminosity and implemented only ammonia clouds.
Nevertheless, we believe incorporating time-changing stellar insolation and
allowing for the formation of water clouds at higher atmospheric temperatures
should have a less significant effect than ammonia clouds, since a giant’s
initial evolution across this parameter space is much more rapid than late-
stage cooling.

Our collection of tables, calibrated using data from Jupiter and covering a
relatively broad range of parameters that includes varying $Z$, $a_{0}$, and
$f$, serves as a foundational resource for advancing our understanding of the
evolution of Jupiter-like exoplanets and for modeling their spectral
evolution. This new boundary-condition dataset is available to the scientific
community for future research into the evolution of both solar-system giants
and giant exoplanets. Nevertheless, the effects of the composition gradients
inferred from Juno data (Wahl et al., 2017; Bolton et al., 2017) and of helium
rain (Stevenson, 1975; Stevenson & Salpeter, 1977; Fortney & Hubbard, 2004;
Mankovich et al., 2016; Mankovich & Fortney, 2020), as well as updates to the
H/He equation of state (Nettelmann et al., 2012; Militzer & Hubbard, 2013;
Miguel et al., 2016; Howard et al., 2023; Howard & Guillot, 2023), demand a
more generalized view beyond the adiabatic paradigm. We used the latter here
merely to test our atmosphere calculations in light of previous work. In
summary, these new boundary tables and atmospheres are meant to support the
next generation of comprehensive giant planet models, which is already well
underway (Nettelmann et al., 2015; Püstow et al., 2016; Vazan et al., 2016;
Mankovich et al., 2016; Vazan et al., 2018; Mankovich & Fortney, 2020; Miguel
& Vazan, 2023).

We thank Brianna Lacy for updated opacity tables and helpful discussions.
Funding (or partial funding) for this research was provided by the Center for
Matter at Atomic Pressures (CMAP), a National Science Foundation (NSF) Physics
Frontier Center, under Award PHY-2020249. Any opinions, findings, conclusions
or recommendations expressed in this material are those of the author(s) and
do not necessarily reflect those of the National Science Foundation.

## References

  * Appleby (1986) Appleby, J. F. 1986, Icarus, 65, 383, doi: 10.1016/0019-1035(86)90145-4
  * Bolton et al. (2017) Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Science, 356, 821, doi: 10.1126/science.aal2108
  * Boss (1997) Boss, A. P. 1997, Science, 276, 1836, doi: 10.1126/science.276.5320.1836
  * Brooke et al. (1998) Brooke, T. Y., Knacke, R. F., Encrenaz, T., et al. 1998, Icarus, 136, 1, doi: 10.1006/icar.1998.6013
  * Burrows et al. (2008) Burrows, A., Budaj, J., & Hubeny, I. 2008, ApJ, 678, 1436, doi: 10.1086/533518
  * Burrows et al. (2001) Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, Reviews of Modern Physics, 73, 719, doi: 10.1103/RevModPhys.73.719
  * Burrows & Sharp (1999) Burrows, A., & Sharp, C. M. 1999, ApJ, 512, 843, doi: 10.1086/306811
  * Burrows et al. (2003) Burrows, A., Sudarsky, D., & Hubbard, W. B. 2003, ApJ, 594, 545, doi: 10.1086/376897
  * Burrows et al. (2006a) Burrows, A., Sudarsky, D., & Hubeny, I. 2006a, ApJ, 650, 1140, doi: 10.1086/507269
  * Burrows et al. (2006b) —. 2006b, ApJ, 640, 1063, doi: 10.1086/500293
  * Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856, doi: 10.1086/305002
  * Cahoy et al. (2010) Cahoy, K. L., Marley, M. S., & Fortney, J. J. 2010, ApJ, 724, 189, doi: 10.1088/0004-637X/724/1/189
  * Carlson et al. (2016) Carlson, R. W., Baines, K. H., Anderson, M. S., Filacchione, G., & Simon, A. A. 2016, Icarus, 274, 106, doi: 10.1016/j.icarus.2016.03.008
  * Carter et al. (2022) Carter, A. L., Hinkley, S., Kammerer, J., et al. 2022, arXiv e-prints, arXiv:2208.14990, doi: 10.48550/arXiv.2208.14990
  * Chubb et al. (2021) Chubb, K. L., Rocchetto, M., Yurchenko, S. N., et al. 2021, A&A, 646, A21, doi: 10.1051/0004-6361/202038350
  * Conrath et al. (1989) Conrath, B. J., Hanel, R. A., & Samuelson, R. E. 1989, in Origin and Evolution of Planetary and Satellite Atmospheres, ed. S. K. Atreya, J. B. Pollack, & M. S. Matthews (Tucson, AZ, University of Arizona Press), 513–538
  * D’Angelo et al. (2003) D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540, doi: 10.1086/367555
  * de Pater et al. (2019) de Pater, I., Sault, R. J., Moeckel, C., et al. 2019, AJ, 158, 139, doi: 10.3847/1538-3881/ab3643
  * Deirmendjian (1964) Deirmendjian, D. 1964, Appl. Opt., 3, 187, doi: 10.1364/AO.3.000187
  * Fortney & Hubbard (2003) Fortney, J. J., & Hubbard, W. B. 2003, Icarus, 164, 228, doi: 10.1016/S0019-1035(03)00130-1
  * Fortney & Hubbard (2004) —. 2004, ApJ, 608, 1039, doi: 10.1086/420765
  * Fortney et al. (2011) Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., & Marley, M. S. 2011, ApJ, 729, 32, doi: 10.1088/0004-637X/729/1/32
  * Gharib-Nezhad et al. (2021) Gharib-Nezhad, E., Iyer, A. R., Line, M. R., et al. 2021, ApJS, 254, 34, doi: 10.3847/1538-4365/abf504
  * Graboske et al. (1975) Graboske, H. C., J., Pollack, J. B., Grossman, A. S., & Olness, R. J. 1975, ApJ, 199, 265, doi: 10.1086/153689
  * Guillot et al. (2004) Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, in Jupiter. The Planet, Satellites and Magnetosphere, ed. F. Bagenal, T. E. Dowling, & W. B. McKinnon, Vol. 1 (Cambrige University Press), 35–57
  * Howard & Guillot (2023) Howard, S., & Guillot, T. 2023, A&A, 672, L1, doi: 10.1051/0004-6361/202244851
  * Howard et al. (2023) Howard, S., Guillot, T., Bazot, M., et al. 2023, A&A, 672, A33, doi: 10.1051/0004-6361/202245625
  * Hubbard (1977) Hubbard, W. B. 1977, Icarus, 30, 305, doi: 10.1016/0019-1035(77)90164-6
  * Hubbard et al. (1999) Hubbard, W. B., Guillot, T., Marley, M. S., et al. 1999, Planet. Space Sci., 47, 1175, doi: 10.1016/S0032-0633(99)00042-2
  * Hubeny & Lanz (1995) Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875, doi: 10.1086/175226
  * Karkoschka (1994) Karkoschka, E. 1994, Icarus, 111, 174, doi: 10.1006/icar.1994.1139
  * Karkoschka (1998) —. 1998, Icarus, 133, 134, doi: 10.1006/icar.1998.5913
  * Kerker (1969) Kerker, M. 1969, The scattering of light and other electromagnetic radiation. (Academic Press)
  * Lacy & Burrows (2023) Lacy, B., & Burrows, A. 2023, arXiv e-prints, arXiv:2303.16295, doi: 10.48550/arXiv.2303.16295
  * Lacy et al. (2019) Lacy, B., Shlivko, D., & Burrows, A. 2019, AJ, 157, 132, doi: 10.3847/1538-3881/ab0415
  * Li et al. (2012) Li, L., Baines, K. H., Smith, M. A., et al. 2012, Journal of Geophysical Research (Planets), 117, E11002, doi: 10.1029/2012JE004191
  * Li et al. (2018) Li, L., Jiang, X., West, R. A., et al. 2018, Nature Communications, 9, 3709, doi: 10.1038/s41467-018-06107-2
  * Li et al. (2021) Li, Y.-P., Chen, Y.-X., Lin, D. N. C., & Zhang, X. 2021, ApJ, 906, 52, doi: 10.3847/1538-4357/abc883
  * Lindal et al. (1981) Lindal, G. F., Wood, G. E., Levy, G. S., et al. 1981, J. Geophys. Res., 86, 8721, doi: 10.1029/JA086iA10p08721
  * Madhusudhan & Burrows (2012) Madhusudhan, N., & Burrows, A. 2012, ApJ, 747, 25, doi: 10.1088/0004-637X/747/1/25
  * Madhusudhan et al. (2011) Madhusudhan, N., Burrows, A., & Currie, T. 2011, ApJ, 737, 34, doi: 10.1088/0004-637X/737/1/34
  * Mankovich et al. (2016) Mankovich, C., Fortney, J. J., & Moore, K. L. 2016, ApJ, 832, 113, doi: 10.3847/0004-637X/832/2/113
  * Mankovich & Fortney (2020) Mankovich, C. R., & Fortney, J. J. 2020, ApJ, 889, 51, doi: 10.3847/1538-4357/ab6210
  * Marley & McKay (1999) Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268, doi: 10.1006/icar.1998.6071
  * Mayor & Queloz (1995) Mayor, M., & Queloz, D. 1995, Nature, 378, 355, doi: 10.1038/378355a0
  * Miguel et al. (2016) Miguel, Y., Guillot, T., & Fayon, L. 2016, A&A, 596, A114, doi: 10.1051/0004-6361/201629732
  * Miguel & Vazan (2023) Miguel, Y., & Vazan, A. 2023, Remote Sensing, 15, 681, doi: 10.3390/rs15030681
  * Miles et al. (2023) Miles, B. E., Biller, B. A., Patapis, P., et al. 2023, ApJ, 946, L6, doi: 10.3847/2041-8213/acb04a
  * Militzer & Hubbard (2013) Militzer, B., & Hubbard, W. B. 2013, ApJ, 774, 148, doi: 10.1088/0004-637X/774/2/148
  * Nettelmann et al. (2012) Nettelmann, N., Becker, A., Holst, B., & Redmer, R. 2012, ApJ, 750, 52, doi: 10.1088/0004-637X/750/1/52
  * Nettelmann et al. (2015) Nettelmann, N., Fortney, J. J., Moore, K., & Mankovich, C. 2015, MNRAS, 447, 3422, doi: 10.1093/mnras/stu2634
  * Niemann et al. (1998) Niemann, H. B., Atreya, S. K., Carignan, G. R., et al. 1998, J. Geophys. Res., 103, 22831, doi: 10.1029/98JE01050
  * Pollack et al. (1977) Pollack, J. B., Grossman, A. S., Moore, R., & Graboske, H. C., J. 1977, Icarus, 30, 111, doi: 10.1016/0019-1035(77)90126-9
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Püstow et al. (2016) Püstow, R., Nettelmann, N., Lorenzen, W., & Redmer, R. 2016, Icarus, 267, 323, doi: 10.1016/j.icarus.2015.12.009
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Saumon et al. (1996) Saumon, D., Hubbard, W. B., Burrows, A., et al. 1996, ApJ, 460, 993, doi: 10.1086/177027
  * Seiff et al. (1998) Seiff, A., Kirk, D. B., Knight, T. C. D., et al. 1998, J. Geophys. Res., 103, 22857, doi: 10.1029/98JE01766
  * Sharp & Burrows (2007) Sharp, C. M., & Burrows, A. 2007, ApJS, 168, 140, doi: 10.1086/508708
  * Sromovsky & Fry (2018) Sromovsky, L. A., & Fry, P. M. 2018, Icarus, 307, 347, doi: 10.1016/j.icarus.2017.10.037
  * Stevenson (1975) Stevenson, D. J. 1975, Phys. Rev. B, 12, 3999, doi: 10.1103/PhysRevB.12.3999
  * Stevenson & Salpeter (1977) Stevenson, D. J., & Salpeter, E. E. 1977, ApJS, 35, 239, doi: 10.1086/190479
  * Sudarsky et al. (2003) Sudarsky, D., Burrows, A., & Hubeny, I. 2003, ApJ, 588, 1121, doi: 10.1086/374331
  * Sudarsky et al. (2005) Sudarsky, D., Burrows, A., Hubeny, I., & Li, A. 2005, ApJ, 627, 520, doi: 10.1086/430206
  * Sudarsky et al. (2000) Sudarsky, D., Burrows, A., & Pinto, P. 2000, ApJ, 538, 885, doi: 10.1086/309160
  * Vazan et al. (2018) Vazan, A., Helled, R., & Guillot, T. 2018, A&A, 610, L14, doi: 10.1051/0004-6361/201732522
  * Vazan et al. (2016) Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118, doi: 10.3847/0004-637X/829/2/118
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649, doi: 10.1002/2017GL073160
  * Yelle et al. (2001) Yelle, R. V., Griffith, C. A., & Young, L. A. 2001, Icarus, 152, 331, doi: 10.1006/icar.2001.6640
  * Young et al. (2005) Young, L. A., Yelle, R. V., Young, R., Seiff, A., & Kirk, D. B. 2005, Icarus, 173, 185, doi: 10.1016/j.icarus.2004.07.020
  * Yurchenko & Tennyson (2014) Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649, doi: 10.1093/mnras/stu326
  * Zhang et al. (2015) Zhang, X., West, R. A., Irwin, P. G. J., Nixon, C. A., & Yung, Y. L. 2015, Nature Communications, 6, 10231, doi: 10.1038/ncomms10231


# WAX-ML: A Python library for machine learning and feedback loops on
streaming data

Emmanuel Sérié<EMAIL_ADDRESS>Capital Fund Management
23 rue de l’Université
75007 Paris, France

###### Abstract

Wax is what you put on a surfboard to avoid slipping. It is an essential tool
to go surfing …We introduce WAX-ML a research-oriented Python library
providing tools to design powerful machine learning algorithms and feedback
loops working on streaming data. It strives to complement JAX with tools
dedicated to time series. WAX-ML makes JAX-based programs easy to use for end-
users working with pandas and xarray for data manipulation. It provides a
simple mechanism for implementing feedback loops, allows the implementation of
online learning and reinforcement learning algorithms with functions, and
makes them easy to integrate by end-users working with the object-oriented
reinforcement learning framework from the Gym library. It is released with an
Apache open-source license on GitHub at https://github.com/eserie/wax-ml.

Keywords: machine-learning, streaming data, time series, feedback loops,
Python, NumPy, JAX, Gym

##  1 Introduction

Since the advent of machine learning and its considerable development over the
last 20 years many tools and computer programs have been developed in this
field. The Python programming language (see Van Rossum and Drake (2009)),
which is a very flexible and easy-to-use programming language, has contributed
greatly to making machine learning programs more accessible.

Since its creation, a whole system of tools has been built to allow the
implementation and access to efficient numerical computation methods. To begin
with, the NumPy library, created in 2005 (see Harris et al. (2020) for a
recent review), allows to perform numerical computations and exposes, with a
simple API, multi-dimensional numerical arrays, and a set of functions
implementing mathematical operations on these arrays.

Subsequently, multiple tools have been built in different and complementary
directions to go further in the type of transformations that can be performed.

The Scikit-learn library, released in 2010 (see Buitinck et al. (2013)), has
proposed a high-level object-oriented API for exposing a very large set of
machine learning methods for supervised and unsupervised problems and has
gathered a large international community around it.

On the other hand, the pandas library (see Wes McKinney (2010)), launched in
2009 , has built, on top of NumPy, a high-level powerful and flexible
dataframe API that is better suited for data analysis and manipulation tools.
In particular, it includes "time-series functionalities such as date range
generation and frequency conversion, moving window statistics, date shifting
and lagging, creation of domain-specific time offsets and joining of time-
series". It also strives to have highly performance-optimized implementations
with critical code paths written in Cython or C.

Some time-series oriented machine learning tools have been released, such as
Prophet (see Taylor and Letham (2017)) in 2017 which is dedicated to end-users
working on time series and wanting to do prediction with black box statistical
modeling techniques.

The TensorFlow library, released in 2015 (see Abadi et al. (2015)), proposed a
set of tools to implement deep neural networks and has since built a large
community of developers and users around it. The TensorFlow developers rebuilt
a tensor API in Python, aiming to be compatible with NumPy but not completely,
which made TensorFlow a separate ecosystem from that built around NumPy. In
particular, this separation sometimes makes it difficult to develop research
ideas in a short development cycle, often requiring reimplementation of NumPy
code in TensorFlow. On the other hand, its widespread adoption in the research
community and industry for the implementation of deep-learning algorithms has
led to the development of very powerful optimization systems such as XLA
(Accelerated Linear Algebra) (see xla ), a domain-specific linear algebra
compiler that can accelerate TensorFlow models with potentially no changes to
the source code, resulting in improvements in speed and memory usage (see mlp
(2020)).

Since then, other libraries have been built to allow the implementation of XLA
programs outside of TensorFlow with other systems such as JAX (see Bradbury et
al. (2018)), Julia (see Bezanson et al. (2017)), PyTorch (see Paszke et al.
(2019)), Nx (see Team ).

In particular, JAX, a library released in 2018 (see Bradbury et al. (2018)),
has initiated a path to solving the "ecosystem separation" problem explained
above, by providing a functional research-oriented API in Python for writing
XLA programs with a NumPy API and primitive transformations to perform just-
in-time XLA compilation, automatic differentiation, parallelization, and
vectorization.

The approach followed by the JAX developers has been to propose some quite
low-level primitives working with well-identified programming constraints,
such as working with pure functions. This constraint seems to be a good
compromise between maintaining the flexibility of Python and having the
possibility to optimize numerical computations with systems like XLA. This
approach suggests writing numerical Python programs in a functional
programming style afterward. This is not very natural in a language like
Python which is at its core an object-oriented programming language. Once the
constraints given by the functional programming approach of JAX are
"digested", this can become a strength since functional programming allows one
to write notably modular, reusable, testable programs, which aim at making JAX
an efficient tool to implement research ideas.

This bet taken by JAX has somehow worked since a whole ecosystem has been
built around it with the advent of domain-specific libraries. For instance,
Optax (see Hessel et al. (2020)), "a gradient processing and optimization
library for JAX", Haiku (see Hennigan et al. (2020)), "a simple neural network
library for JAX", RLax (see Budden et al. (2020)), "a library built on top of
JAX that exposes useful building blocks for implementing reinforcement
learning agents", Flax (see Heek et al. (2020)), "a neural network library and
ecosystem for JAX designed for flexibility", …

With the launch of WAX-ML, we want to complement this growing ecosystem with
tools that make it easier to implement machine learning algorithms on
streaming data.

We want to push forward the efforts initiated by JAX developers to reconcile
the Python ecosystems around numerical computations by proposing tools
allowing the use of JAX programs on time-series data thanks to the use of the
pandas and xarray (see Hoyer and Hamman (2017)) libraries which are, for now,
reference libraries for the manipulation of labeled datasets and in particular
time-series data.

We even think that the TensorFlow and PyTorch ecosystem could also be made
compatible with programs implemented with JAX ecosystem tools thanks to
EagerPy (see Rauber et al. (2020)) which is a library that allows writing code
that works natively with PyTorch, TensorFlow, JAX, and NumPy. In this first
released version of WAX-ML, we use EagerPy to facilitate the conversion
between NumPy tensors (which are used by pandas and xarray) and JAX tensors.
We are also implementing some universal modules so that we can work with other
tensor libraries, but for now, we remain focused on running JAX programs on
streaming data.

##  2 WAX-ML goal

WAX-ML’s goal is to expose "traditional" algorithms that are often difficult
to find in standard Python ecosystem and are related to time-series and more
generally to streaming data.

It aims to make it easy to work with algorithms from very various
computational domains such as machine learning, online learning, reinforcement
learning, optimal control, time-series analysis, optimization, statistical
modeling.

For now, WAX-ML focuses on time-series algorithms as this is one of the areas
of machine learning that lacks the most dedicated tools. Working with time
series is notoriously known to be difficult and often requires very specific
algorithms (statistical modeling, filtering, optimal control).

Even though some of the modern machine learning methods such as RNN, LSTM, or
reinforcement learning can do an excellent job on some specific time-series
problems, most of the problems require using more traditional algorithms such
as linear and non-linear filters, FFT, the eigendecomposition of matrices
(e.g. Bouchaud et al. (2005)), principal component analysis (PCA) (e.g. Dong
and Qin (2018)), Riccati solvers for optimal control and filtering,…

By adopting a functional approach, inherited from JAX, WAX-ML aims to be an
efficient tool to combine modern machine learning approaches with more
traditional ones.

Some work has been done in this direction in (Lee-Thorp et al. (2021)) where
transformer encoder architectures are massively accelerated, with limited
accuracy costs, by replacing the self-attention sublayers with a standard,
non-parameterized Fast Fourier Transform (FFT).

WAX-ML may also be useful for developing research ideas in areas such as
online machine learning (see haz (2021); Hazan (2019)) and development of
control, reinforcement learning, and online optimization methods.

##  3 What does WAX-ML do?

Well, building WAX-ML, we had some pretty ambitious design and implementation
goals.

To do things right, we decided to start small and in an open-source design
from the beginning. We released it with an Apache open-source license (see
Sérié (2021)).

In this section, we give a quick overview of the features we have developed in
WAX-ML. They are described in more detail in section 3.7.

For now, WAX-ML contains transformation tools that we call "unroll"
transformations allowing us to apply any transformation, possibly stateful, on
sequential data. It generalizes the RNN architecture to any stateful
transformation allowing the implementation of any kind of "filter".

We have implemented a "stream" module, described in section 3.3, permitting us
to synchronize data streams with different time resolutions.

We have implemented some general pandas and xarray "accessors" permitting the
application of any JAX-functions on pandas and xarray data containers.

We have implemented a ready-to-use exponential moving average filter that we
exposed with two APIs. One for JAX users: as Haiku modules (see for instance
our `EWMA` module). A second one for pandas and xarray users: with drop-in
replacement of pandas `ewm` accessor.

We have implemented a simple module `OnlineSupervisedLearner` to implement
online learning algorithms for supervised machine learning problems.

We have implemented building blocks for designing feedback loops in
reinforcement learning, and have provided a module called `GymFeedback`
allowing the implementation of feedback loop as the introduced in the library
Gym (see Brockman et al. (2016)), and illustrated in Figure 4.

Finally, we have implemented some “universal” modules that can work with
TensorFlow, PyTorch, JAX, and NumPy tensors. At the moment, we have only
implemented a demonstration module for the exponential moving average that we
have called `EagerEWMA`.

###  3.1 Why use WAX-ML?

If you deal with time-series and are a pandas or xarray user, but you want to
use the impressive tools of the JAX ecosystem, then WAX-ML might be the right
tool for you, as it implements pandas and xarray accessors to apply JAX
functions.

If you are a user of JAX, you may be interested in adding WAX-ML to your
toolbox to address time-series problems.

###  3.2 Functional programming

In WAX-ML, we pursue a functional programming approach inherited from JAX.

In this sense, WAX-ML is not a framework, as most object-oriented libraries
offer. Instead, we implement "functions" that must be pure to exploit the JAX
ecosystem.

We use the "module" mechanism proposed by the Haiku library to easily generate
pure function pairs, called `init` and `apply` in Haiku, to implement programs
that require the management of parameters and/or state variables. In this way,
we can recover all the advantages of object-oriented programming but exposed
in the functional programming approach.

This approach gives a lot of freedom in the type of ideas that can be
implemented. For instance, JAX has been used recently to accelerate fluid
dynamics simulations (see Kochkov et al. (2021)) by two orders of magnitude.

WAX-ML does not want to reinvent the wheel by reimplementing every algorithm.
We want existing machine learning libraries to work well together while trying
to leverage their strength, which is easy to do with a functional programming
approach.

To demonstrate this, in the current version of WAX-ML, we have constructed
various examples, such as an exponential moving average (to serve as a toy
example), the implementation and calibration of an LSTM architecture with a
standard supervised machine learning workflow, and the implementation of
online learning and reinforcement learning architectures. They are treated
equally and laid out with flat organization in our sub-package `wax.modules`.

###  3.3 Synchronize streams

Physicists have brought a solution to the synchronization problem called the
Poincaré–Einstein synchronization (see synchronisation ). In WAX-ML we have
implemented a similar mechanism by defining a "local time", borrowing Henri
Poincaré terminology, to denominate the timestamps of the stream (the "local
stream") in which the user wants to apply transformations and unravel all
other streams. The other streams, which we have called "secondary streams",
are pushed back in the local stream using embedding maps which specify how to
convert timestamps from a secondary stream into timestamps in the local
stream.

This synchronization mechanism permits to work with secondary streams having
timestamps at frequencies that can be lower or higher than the local stream.
The data from these secondary streams are represented in the "local stream"
either with the use of a forward filling mechanism for lower frequencies or
with a buffering mechanism for higher frequencies.

Note that this simple synchronization scheme assumes that the different event
streams have fixed latencies.

We have implemented a "data tracing" mechanism to optimize access to out-of-
sync streams. This mechanism works on in-memory data. We perform the first
pass on the data, without actually accessing it, and determine the indices
necessary to later access the data. Doing so we are vigilant to not let any
"future" information pass through and thus guaranty a data processing that
respects causality.

The buffering mechanism used in the case of higher frequencies works with a
fixed buffer size (see the WAX-ML module `wax.modules.Buffer`) to allow the
use of JAX / XLA optimizations and efficient processing.

We give simple usage examples in our documentation and in working example
shown in Figure 2.

###  3.4 Working with streaming data

WAX-ML may complement JAX ecosystem by adding support for streaming data.

To do this, we implement a unique data tracing mechanism that prepares for
fast access to in-memory data and allows the execution of JAX tractable
functions such as `jit`, `grad`, `vmap` or `pmap` (see section 3.3).

This mechanism is somewhat special in that it works with time-series data.

The `wax.stream.Stream` object implements this idea. It uses Python generators
to synchronize multiple streaming data streams with potentially different
temporal resolutions. It works on in-memory data stored in `xarray.Dataset`.

To work with real streaming data, it should be possible to implement a buffer
mechanism running on any Python generator and to use the synchronization and
data tracing mechanisms implemented in WAX-ML to apply JAX transformations on
batches of data stored in memory. We have not yet implemented such a mechanism
but put it on our enhancement proposal list.

###  3.5 Adding support for time types in JAX

At the moment `datetime64` and `string_` types are not supported in JAX.

WAX-ML add support for these NumPy types in JAX by implementing an encoding
scheme for `datetime64` relying on pairs of 32-bit integers similar to
`PRNGKey` implemented in JAX for pseudo-random number generation. We also
implement an encoding scheme for `string_` relying on `LabelEncoder` from the
Scikit-learn library.

By providing these two encoding schemes, WAX-ML makes it easy to use JAX
algorithms on data of these types.

Currently, the types of time offsets supported by WAX-ML are quite limited and
we would like to collaborate with the pandas, xarray, and Astropy (see Astropy
Collaboration et al. (2013)) teams to develop further the time manipulation
tools in WAX-ML.

###  3.6 pandas and xarray accessors

WAX-ML implements pandas and xarray accessors to ease the usage of machine-
learning algorithms implemented with JAX-functions within high-level data APIs
such as `DataFrame` and `Series` in pandas and `Dataset` and `DataArray` in
xarray.

To load the accessors, users should have to simply run
`register_wax_accessors()` function and then use the "one-liner" syntax
`<data-container>.stream(…).apply(…)`

###  3.7 Implemented modules

We have some modules (inherited from Haiku modules) ready to be used in
`wax.modules`. They can be considered as "building blocks" that can be reused
to build more advanced programs to run on streaming data. We have some
"fundamental" modules that are specific to time series management, the
‘Buffer‘ module which implements the buffering mechanism, the ‘UpdateOnEvent‘
module which allows one to "freeze" the computations of a program and to
update them on some events in the "local stream".

To illustrate the use of this module, we show in Figure 1 how it can be used
to compute the trailing open, high, and close quantities of temperatures
recorded during the day.

Figure 1: Computation of trailing Open-High-Low-Close temperatures.

We have a few more specific modules that aim to reproduce some of the logic
that pandas users may be familiar with, such as `Lag` to implement a delay on
the input data, `Diff` to compute differences between values over time,
`PctChange` to compute the relative difference between values over time,
`RollingMean` to compute the moving average over time., `EWMA`, `EWMVar`,
`EWMCov`, to compute the exponential moving average, variance, and covariance
of the input data.

Finally, we implement domain-specific modules for online learning and
reinforcement learning such as `OnlineSupervisedLearner` (see section 4.2) and
`GymFeedback` (see section 3.9).

For now, WAX-ML offers direct access to some modules through specific
accessors for pandas and xarray users. For instance, we have an implementation
of the "exponential moving average" directly accessible through the accessor
`<data-container>.ewm(...).mean()` which provides a drop-in replacement for
the exponential moving average of pandas.

We have implemented some working examples in our documentation. In Figure 2,
we show an illustrative example where we apply exponential moving average
filters with different time scales on a temperature data set (borrowed from
the xarray library) while synchronizing two data streams ("air" and "ground")
operating at different frequencies.

Figure 2: We apply exponential moving average filters with different time
scales while synchronizing two data streams ("air" and "ground") operating at
different frequencies.

###  3.8 Speed

The use of JAX allows for leveraging hardware accelerators that optimize
programs for the CPU, GPU, and TPU architectures.

We have performed few benchmarks with WAX-ML on the computation of exponential
moving averages on rather large dataframes (with 1 million rows and 1000
columns). For this size of dataframe, we observed a speedup of 3 compared to
native pandas implementation. This is not a very important factor but it shows
that even on such a simple algorithm already implemented in C, the XLA
compilation can optimize things. The details of this benchmark are available
in our documentation.

###  3.9 Feedback loops

Feedback is a fundamental notion in time-series analysis and has a wide
history (e.g. Feedback ). So, we believe it is important to be able to
implement them well in WAX-ML.

A fundamental piece in the implementation of feedback loops is the delay
operator. We implement it with the delay module `Lag` which is itself
implemented with the `Buffer` module, a module implementing the buffering
mechanism.

The linear state-space models used to model linear time-invariant systems in
signal theory are a well-known place where feedbacks are used to implement for
instance infinite impulse response filters. This is easily implemented with
the WAX-ML tools and will be implemented at a later time.

Another example is control theory or reinforcement learning. In reinforcement
learning setup, an agent and an environment interact with a feedback loop (see
Figure 4). This generally results in a non-trivial global dynamic. In WAX-ML,
we propose a simple module called `GymFeedBack` that allows the implementation
of reinforcement learning experiments. This is built from an agent and an
environment, both possibly having parameters and state (see Figure 3). The
agent is in charge of generating an action from observations. The environment
is in charge of calculating a reward associated with the agent’s action and
preparing the next observation from some "raw observations" and the agent’s
action, which it gives back to the agent.

A feedback instance `GymFeedback(agent, env)` is a function that processes the
"raw observations" and returns a reward as represented in the Figure 4.
Equivalently, we can describe the function `GymFeedback(agent, env)`, after
transformation by Haiku transformation, by a pair of pure functions `init` and
`apply` that we describe in Figure 5.

We have made concrete use of this feedback mechanism in the online learning
application of section 4.2.

Figure 3: Agent and Environment are simple functions (possibly with parameters
and state). Figure 4: Description of a Gym feedback loop. Figure 5:
Description of a Gym feedback loop in term of pure functions init and apply.

###  3.10 Compatibility with other reinforcement learning frameworks

In addition, to ensuring compatibility with other tools in the Gym ecosystem
(see Brockman et al. (2016)), we propose a transformation mechanism to
transform functions into standard stateful Python objects following the Gym
API for agents and environments as implemented in the `deluca` framework (see
LLC (2021)).

##  4 Applications

###  4.1 Reconstructing the light curve of stars with LSTM

To demonstrate more advanced usages of WAX-ML with machine learning tools on
time series we have reproduced the study (see Pere (2020)) on the
reconstruction of the light curve of stars. We used the LSTM architecture to
predict the observed light flux through time and use the Haiku implementation.
We have then reproduced Haiku’s case study for LSTM on this concrete
application on the light curve of stars. We show an illustrative plot of the
study in Figure 6.

Here WAX-ML is used to transition efficiently from time-series data with Nan
values stored in a dataframe to a "standard" deep learning workflow in Haiku.

The study is available in our documentation through a reproducible notebook.

As a disclaimer, we emphasize that although our study uses real data, the
results presented in it should not be taken as scientific knowledge, since
neither the results nor the data source has been peer-reviewed by
astrophysists. The purpose of the study was only to demonstrate how WAX-ML can
be used when applying a "standard" machine-learning workflow, here LSTM, to
analyze time series. However, we welcome further collaboration to continue the
study.

Figure 6: light curve of stars reconstruction with LSTM. First: training
curve. Second: auto-regressively predicted values and truth values.

###  4.2 Online linear regression in non-stationary environment

We performed a case study by implementing an online learning non-stationary
linear regression problem following online learning techniques exposed in
(Hazan (2019)).

We go there progressively by showing how a linear regression problem can be
cast into an online learning problem thanks to the `OnlineSupervisedLearner`
module developed in WAX-ML.

Then, to tackle a non-stationary linear regression problem (i.e. with a weight
that can vary in time) we reformulate the problem into a reinforcement
learning problem that we implement with the `GymFeedBack` module of WAX-ML.

We then define an "agent" and an "environment" using simple functions
implemented with modules as illustrated in Figure 3. We have the agent which
is responsible for learning the weights of its internal linear model and the
environment which is responsible for generating labels and evaluating the
agent’s reward metric. We finally assemble them in a gym feedback loop as in
Figure 4 and Figure 5.

We experiment with a non-stationary environment that flips the sign of the
linear regression parameters at a given time step, known only to the
environment. We show an illustrative plot of the final result of the study in
Figure 7. More details about the study can be found in our
documentation111https://wax-ml.readthedocs.io/en/latest/.

This example shows that it is quite simple to implement the online learning
task with WAX-ML tools. In particular, the functional programming approach
allowed us to compose and to reuse functions to design the final task.

Figure 7: Online linear regression in a non-stationary environment. Left: The
regret (cumulative sum of losses) first becomes concave, which means that the
agent "learns something". Then, the regret curve has a bump at step 2000 where
it becomes locally linear. It finally ends in a concave regime concave regime,
which means that the agent has adapted to the new regime. Right: We see that
the weights converge to the correct values in both regimes.

##  5 Future plans

###  5.1 Feedback loops and control theory

We would like to implement other types of feedback loops in WAX-ML. For
instance, those of the standard control theory toolboxes, such as those
implemented in the SLICOT (see Benner et al. (1999)) library.

Many algorithms in this space are absent from the Python ecosystem and we aim
to provide JAX-based implementations and expose them with a simple API.

An idiomatic example in this field is the Kalman filter, a now-standard
algorithm that dates back to the 1950s (see Kalman (1960)). After 30 years of
existence, the Python ecosystem has still not integrated this algorithm into
widely adopted libraries. Some implementations can be found in Python
libraries such as `python-control` (see pyt ), statsmodels (see Seabold and
Perktold (2010)), or SciPy (see Virtanen et al. (2020)). Also, some machine
learning libraries such as Scikit-learn (see Buitinck et al. (2013)) and River
(see Montiel et al. (2020)) have some closed and non-solved issues on this
subject. Why has the Kalman filter not found its place in these libraries? We
think it may be because they have an object-oriented API, which makes them
very well suited to the specific problems of modern machine learning but, on
the other hand, prevents them from accommodating additional features such as
Kalman filtering. We think the functional approach of WAX-ML, inherited from
JAX, could well help to integrate a Kalman filter implementation in a machine
learning ecosystem.

It turns out that Python code written with JAX is not very far from Fortran
(see Fortran ), a mathematical FORmula TRANslating system. It should therefore
be quite easy and natural to reimplement standard algorithms implemented in
Fortran, such as those in the SLICOT library with JAX. It seems that some
questions about the integration of Fortran into JAX have already been raised.
As noted in an issue222issue 3950 reported on JAX’s GitHub page, and it might
even be possible to simply wrap Fortran code in JAX, which would avoid a
painful rewriting process.

Along with the implementation of good old algorithms, we would like to
implement more recent ones from the online learning literature which somehow
revisit the filtering and control problems. In particular, we would like to
implement the online learning version of the ARMA model developed in Anava et
al. (2013) and some online-learning versions of control theory algorithms, an
approach called "the non-stochastic control problem", such as the linear
quadratic regulator developed in Hazan et al. (2020).

###  5.2 Optimization

The JAX ecosystem already has a library dedicated to optimization: Optax,
which we use in WAX-ML. We could complete it by offering other first-order
algorithms such as the Alternating Direction Multiplier Method (ADMM) (e.g.
Boyd et al. (2011)). One can find functional implementations of proximal
algorithms in libraries such as `proxmin` (see Melchior (2021)),
`ProximalOperators` (see Stella et al. (2021)), or COSMO (see Garstka et al.
(2019)), which could give good reference implementations to start the work.

Another type of work took place around automatic differentiation and
optimization. In (Agrawal et al. (2019)) the authors implement differentiable
layers based on convex optimization in the library `cvxpylayers`. They have
implemented a JAX API but, at the moment, they cannot use the `jit`
compilation of JAX yet. We would be interested in helping to solve this
issue333issue 103 JAX’s GitHub page.

Furthermore, in the recent paper by Blondel et al. (2021), the authors propose
a new efficient and modular implicit differentiation technique with a JAX-
based implementation that should lead to a new open-source optimization
library in the JAX ecosystem.

###  5.3 Other algorithms

The machine learning libraries SciPy (see Virtanen et al. (2020)), Scikit-
learn (see Buitinck et al. (2013)), River (see Montiel et al. (2020)), and ml-
numpy (see Bourgin (2021)) implement many "traditional" machine learning
algorithms that should provide an excellent basis for linking or
reimplementing in JAX. WAX-ML could help to build a repository for JAX
versions of these algorithms.

##  6 Conclusion

We think that the simple tools we propose in WAX-ML might be useful for
performing machine-learning on streaming data and implementing feedback loops.

Some work could be continued to further develop the JAX ecosystem and
implement some of the algorithms missing in the Python ecosystem in the wide
adopted libraries.

We hope that this software will be well received by the open-source community
and that its development will also be actively pursued.

##  7 Collaborations

The WAX-ML team is open to discussion and collaboration with contributors from
any field who are interested in using WAX-ML for their problems on streaming
data. We are looking for use cases around data streaming in audio processing,
natural language processing, astrophysics, biology, finance, engineering …

We believe that good software design, especially in the scientific domain,
requires practical use cases and that the more diversified these use cases
are, the more the developed functionalities will be guaranteed to be well
implemented.

By making this software public, we hope to find enthusiasts who aim to develop
WAX-ML further!

##  8 Acknowledgements

I thank my family for supporting me during this period of writing the first
version of WAX-ML.

I thank my brother, Julien Sérié, musician, videographer, and talented graphic
designer, who made the logo of the software.

I would like to thank my colleagues at Capital Fund Management who encouraged
me to pursue this project and gave me constructive feedback while reading
early drafts of this document. I thanks Richard Beneyton, Raphael Benichou,
Jean-Philippe Bouchaud, David Eng, Tifenn Juguet, Julien Lafaye, Laurent
Laloux, Eric Lebigot, Marc Potters, Lamine Souiki, Marc Wouts, and the Open-
source committee of CFM for their early support. I thank Richard Beneyton,
Gawain Bolton and Gilles Zerah for having accepted to review the paper.

## References

  * (1) Python Control Systems Library.  URL https://github.com/python-control/python-control.
  * (2) XLA: Accelerated Linear Algebra.  URL https://www.tensorflow.org/xla.
  * mlp (2020) MLPerf, 2020.  URL https://blog.tensorflow.org/2020/07/tensorflow-2-mlperf-submissions.html.
  * haz (2021) Google Princeton AI and Hazan Lab @ Princeton University, 2021.  URL https://www.minregret.com/research/.
  * Abadi et al. (2015) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, 2015\.  URL http://tensorflow.org/.  Software available from tensorflow.org.
  * Agrawal et al. (2019) Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and Zico Kolter.  Differentiable Convex Optimization Layers, 2019.
  * Anava et al. (2013) Oren Anava, Elad Hazan, Shie Mannor, and Ohad Shamir.  Online Learning for Time Series Prediction.  In Shai Shalev-Shwartz and Ingo Steinwart, editors, _Proceedings of the 26th Annual Conference on Learning Theory_ , volume 30 of _Proceedings of Machine Learning Research_ , pages 172–184, Princeton, NJ, USA, 12–14 Jun 2013. Pmlr.  URL http://proceedings.mlr.press/v30/Anava13.html.
  * Astropy Collaboration et al. (2013) Astropy Collaboration, T. P. Robitaille, E. J. Tollerud, P. Greenfield, M. Droettboom, E. Bray, T. Aldcroft, M. Davis, A. Ginsburg, A. M. Price-Whelan, W. E. Kerzendorf, A. Conley, N. Crighton, K. Barbary, D. Muna, H. Ferguson, F. Grollier, M. M. Parikh, P. H. Nair, H. M. Unther, C. Deil, J. Woillez, S. Conseil, R. Kramer, J. E. H. Turner, L. Singer, R. Fox, B. A. Weaver, V. Zabalza, Z. I. Edwards, K. Azalee Bostroem, D. J. Burke, A. R. Casey, S. M. Crawford, N. Dencheva, J. Ely, T. Jenness, K. Labrie, P. L. Lim, F. Pierfederici, A. Pontzen, A. Ptak, B. Refsdal, M. Servillat, and O. Streicher.  Astropy: A community Python package for astronomy.  _aap_ , 558:A33, October 2013.  doi: 10.1051/0004-6361/201322068.
  * Benner et al. (1999) Peter Benner, Volker Mehrmann, Vasile Sima, Sabine Van Huffel, and Andras Varga.  _SLICOT—A Subroutine Library in Systems and Control Theory_ , pages 499–539.  Birkhäuser Boston, Boston, MA, 1999.  ISBN 978-1-4612-0571-5.  doi: 10.1007/978-1-4612-0571-5_10.  URL https://doi.org/10.1007/978-1-4612-0571-5%5F10.
  * Bezanson et al. (2017) Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah.  Julia: A fresh approach to numerical computing.  _SIAM review_ , 59(1):65–98, 2017.  URL https://doi.org/10.1137/141000671.
  * Blondel et al. (2021) Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares-López, Fabian Pedregosa, and Jean-Philippe Vert.  Efficient and Modular Implicit Differentiation, 2021.
  * Bouchaud et al. (2005) Jean-Philippe Bouchaud, Laurent Laloux, M. Augusta Miceli, and Marc Potters.  Large dimension forecasting models and random singular value spectra, 2005.
  * Bourgin (2021) David Bourgin.  ml-numpy, 2021.  URL https://github.com/ddbourgin/numpy-ml.
  * Boyd et al. (2011) Stephen P. Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  _Foundations and Trends in Machine Learning_ , 3(1):1–122, 2011.  URL http://dblp.uni-trier.de/db/journals/ftml/ftml3.html#BoydPCPE11.
  * Bradbury et al. (2018) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018.  URL http://github.com/google/jax.
  * Brockman et al. (2016) Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.  OpenAI Gym, 2016.
  * Budden et al. (2020) David Budden, Matteo Hessel, John Quan, Steven Kapturowski, Kate Baumli, Surya Bhupatiraju, Aurelia Guy, and Michael King.  RLax: Reinforcement Learning in JAX, 2020.  URL http://github.com/deepmind/rlax.
  * Buitinck et al. (2013) Lars Buitinck, Gilles Louppe, Mathieu Blondel, Fabian Pedregosa, Andreas Mueller, Olivier Grisel, Vlad Niculae, Peter Prettenhofer, Alexandre Gramfort, Jaques Grobler, Robert Layton, Jake VanderPlas, Arnaud Joly, Brian Holt, and Gaël Varoquaux.  API design for machine learning software: experiences from the Scikit-learn project.  In _ECML PKDD Workshop: Languages for Data Mining and Machine Learning_ , pages 108–122, 2013.
  * Dong and Qin (2018) Yining Dong and S. Joe Qin.  A novel dynamic PCA algorithm for dynamic data modeling and process monitoring.  _Journal of Process Control_ , 67:1–11, 2018.  ISSN 0959-1524.  doi: https://doi.org/10.1016/j.jprocont.2017.05.002.  URL https://www.sciencedirect.com/science/article/pii/S095915241730094X.  Big Data: Data Science for Process Control and Operations.
  * (20) Feedback.  Feedback — Wikipedia, The Free Encyclopedia.  URL https://en.wikipedia.org/wiki/Feedback.  [Online; accessed 29-May-2021].
  * (21) Fortran.  Fortran — Wikipedia, The Free Encyclopedia.  URL https://fr.wikipedia.org/wiki/Fortran.  [Online; accessed 29-May-2021].
  * Garstka et al. (2019) Michael Garstka, Mark Cannon, and Paul Goulart.  COSMO: A conic operator splitting method for large convex problems.  In _European Control Conference_ , 2019.  doi: 10.23919/ecc.2019.8796161.  URL https://arxiv.org/abs/1901.10887.
  * Harris et al. (2020) Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  doi: 10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Hazan (2019) Elad Hazan.  Introduction to Online Convex Optimization, 2019.
  * Hazan et al. (2020) Elad Hazan, Sham M. Kakade, and Karan Singh.  The Nonstochastic Control Problem, 2020.
  * Heek et al. (2020) Jonathan Heek, Anselm Levskaya, Avital Oliver, Marvin Ritter, Bertrand Rondepierre, Andreas Steiner, and Marc van Zee.  Flax: A neural network library and ecosystem for JAX, 2020.  URL http://github.com/google/flax.
  * Hennigan et al. (2020) Tom Hennigan, Trevor Cai, Tamara Norman, and Igor Babuschkin.  Haiku: Sonnet for JAX, 2020.  URL http://github.com/deepmind/dm-haiku.
  * Hessel et al. (2020) Matteo Hessel, David Budden, Fabio Viola, Mihaela Rosca, Eren Sezener, and Tom Hennigan.  Optax: composable gradient transformation and optimisation, in JAX!, 2020.  URL http://github.com/deepmind/optax.
  * Hoyer and Hamman (2017) S. Hoyer and J. Hamman.  xarray: N-D labeled arrays and datasets in Python.  _Journal of Open Research Software_ , 5(1), 2017.  doi: 10.5334/jors.148.  URL http://doi.org/10.5334/jors.148.
  * Kalman (1960) Rudolph Emil Kalman.  A New Approach to Linear Filtering and Prediction Problems.  _Journal of Basic Engineering_ , 82(1):35, 1960\.  doi: 10.1115/1.3662552.  URL http://dx.doi.org/10.1115/1.3662552.
  * Kochkov et al. (2021) Dmitrii Kochkov, Jamie A. Smith, Ayya Alieva, Qing Wang, Michael P. Brenner, and Stephan Hoyer.  Machine learning accelerated computational fluid dynamics, 2021.
  * Lee-Thorp et al. (2021) James Lee-Thorp, Joshua Ainslie, Ilya Eckstein, and Santiago Ontanon.  FNet: Mixing Tokens with Fourier Transforms, 2021.
  * LLC (2021) Google LLC.  Performant, differentiable reinforcement learning, 2021.  URL https://github.com/google/deluca.
  * Melchior (2021) Peter Melchior.  Proximal Minimization, 2021.  URL https://github.com/pmelchior/proxmin.
  * Montiel et al. (2020) Jacob Montiel, Max Halford, Saulo Martiello Mastelini, Geoffrey Bolmier, Raphael Sourty, Robin Vaysse, Adil Zouitine, Heitor Murilo Gomes, Jesse Read, Talel Abdessalem, and Albert Bifet.  River: machine learning for streaming data in Python, 2020.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Pere (2020) Christophe Pere.  https://towardsdatascience.com/how-to-use-deep-learning-for-time-series-forecasting-3f8a399cf205, 2020\.
  * Rauber et al. (2020) Jonas Rauber, Matthias Bethge, and Wieland Brendel.  EagerPy: Writing Code That Works Natively with PyTorch, TensorFlow, JAX, and NumPy, 2020.  URL https://eagerpy.jonasrauber.de.
  * Seabold and Perktold (2010) Skipper Seabold and Josef Perktold.  statsmodels: Econometric and statistical modeling with python.  In _9th Python in Science Conference_ , 2010.
  * Stella et al. (2021) Lorenzo Stella, Niccolò Antonello, and Mattias Fält.  ProximalOperators.jl, 2021.  URL https://kul-forbes.github.io/ProximalOperators.jl/latest/.
  * (41) Einstein synchronisation.  Einstein synchronisation — Wikipedia, The Free Encyclopedia.  URL https://en.wikipedia.org/wiki/Einstein%5Fsynchronisation.  [Online; accessed 29-May-2021].
  * Sérié (2021) Emmanuel Sérié.  WAX-ML: A Python library for machine-learning and feedback loops on streaming data, 2021.  URL http://github.com/eserie/wax-ml.
  * Taylor and Letham (2017) S. Taylor and Benjamin Letham.  Forecasting at Scale.  _PeerJ Prepr._ , 5:e3190, 2017.
  * (44) Nx Team.  Nx is a multi-dimensional tensors library for Elixir with multi-staged compilation to the CPU/GPU.  URL https://github.com/elixir-nx/nx/tree/main/nx.
  * Van Rossum and Drake (2009) Guido Van Rossum and Fred L. Drake.  _Python 3 Reference Manual_.  CreateSpace, Scotts Valley, CA, 2009.  ISBN 1441412697.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Wes McKinney (2010) Wes McKinney.  Data Structures for Statistical Computing in Python.  In Stéfan van der Walt and Jarrod Millman, editors, _Proceedings of the 9th Python in Science Conference_ , pages 56–61, 2010.  doi: 10.25080/Majora-92bf1922-00a.