# Terminator Habitability: the Case for Limited Water Availability on M-dwarf
Planets

Ana H. Lobo University of California, Irvine
Department of Physics & Astronomy
4129 Frederick Reines Hall
Irvine, CA 92697-4575, USA  Aomawa L. Shields University of California, Irvine
Department of Physics & Astronomy
4129 Frederick Reines Hall
Irvine, CA 92697-4575, USA  NASA NExSS Virtual Planetary Laboratory, Seattle,
WA  Igor Z. Palubski  University of California, Irvine
Department of Physics & Astronomy
4129 Frederick Reines Hall
Irvine, CA 92697-4575, USA  Eric Wolf University of Colorado, Boulder
Laboratory for Atmospheric and Space Physics,
Department of Atmospheric and Oceanic Sciences,
Boulder, CO  NASA NExSS Virtual Planetary Laboratory, Seattle, WA  NASA GSFC
Sellers Exoplanet Environments Collaboration, Greenbelt, MD

###### Abstract

Rocky planets orbiting M-dwarf stars are among the most promising and abundant
astronomical targets for detecting habitable climates. Planets in the M-dwarf
habitable zone are likely synchronously rotating, such that we expect
significant day-night temperature differences, and potentially limited
fractional habitability. Previous studies have focused on scenarios where
fractional habitability is confined to the substellar or “eye” region, but in
this paper we explore the possibility of planets with terminator habitability,
defined by the existence of a habitable band at the transition between a
scorching dayside and a glacial nightside. Using a global climate model, we
show that for water-limited planets it is possible to have scorching
temperatures in the “eye” and freezing temperatures on the nightside, while
maintaining a temperate climate in the terminator region, due to a reduced
atmospheric energy transport. Whereas on water-rich planets, increasing
stellar flux leads to increased atmospheric energy transport and a reduction
in day-night temperature differences, such that the terminator does not remain
habitable once the dayside temperatures approach runaway or moist greenhouse
limits. We also show that, while water-abundant simulations may result in
larger fractional habitability, they are vulnerable to water loss through
cold-trapping on the nightside surface or atmospheric water vapor escape,
suggesting that even if planets were formed with abundant water, their
climates could become water-limited and subject to terminator habitability.

planets and satellites: terrestrial planets — planets and satellites:
atmospheres

##  1 Introduction

Rocky planets in the habitable zones of M-dwarf stars are currently among the
most promising astronomical targets for planetary climate characterization
studies and the detection of habitable surface climate. M dwarfs, which make
up $\sim$70% of all stars (Bochanski et al., 2010), are expected to have
abundant rocky Earth-sized planets (Mulders et al., 2015). Also, their reduced
luminosity results in a habitable zone that is only $\sim$0.2 au from the star
(Kasting et al., 1993; Kopparapu et al., 2013), which facilitates detection
(Nutzman & Charbonneau, 2008). While M-dwarf prospects for habitability have
been extensively debated (Scalo et al., 2007; Tarter et al., 2007; Shields et
al., 2016), over the course of the last decade many of the most promising
candidates for habitability have been discovered orbiting M dwarfs, such as
TRAPPIST-1 (Gillon et al., 2017), Proxima Centauri (Anglada-Escudé et al.,
2016), and TOI-700 (Gilbert et al., 2020). With the recent launch and
deployment of the James Webb Space Telescope, it may soon be possible to
characterize the atmospheres of terrestrial planets in the habitable zones of
M dwarfs (Kreidberg & Loeb, 2016; Morley et al., 2017; Fauchez et al., 2019).
In preparation for upcoming observations, it is increasingly important that we
understand the full range of possible M-dwarf planetary climates and their
prospects for habitability.

The habitable zone’s close proximity to the star implies that the planets are
subject to rapid tidal locking, and are likely to be synchronously rotating
such that they have a permanent dayside (Armitage, 2009; Barnes, 2017).
Without incoming radiation, the nightside climate is determined by atmospheric
and ocean energy transport, and could be subject to atmospheric collapse
(Joshi et al., 1997). Whereas, on the highly irradiated dayside, greenhouse
effects could be amplified with increased atmospheric water vapor content
(Ingersoll, 1969), potentially leading to a runaway state. Despite the
nightside’s inevitable top-of-atmosphere radiative deficit, if the day-to-
night energy transport is sufficiently intense, it is possible for a large
portions of a synchronously rotating planet to remain above freezing (e.g. Hu
& Yang, 2014; Genio et al., 2019; Yang et al., 2019b), without entering a
runaway greenhouse state. Even neglecting ocean transport, aquaplanet
simulations can achieve moderate temperatures on synchronously rotating
planets (e.g. Joshi, 2003; Merlis & Schneider, 2010; Edson et al., 2011). A
temperate climate is achievable even when accounting for the longer wavelength
stellar spectra of M dwarfs, which results in changes to key radiative
feedbacks (e.g. Shields et al., 2016; Kopparapu et al., 2016; Turbet et al.,
2016; Komacek & Abbot, 2019). Yet, while we can easily simulate an arbitrary
habitable M-dwarf planet, the range of habitable planetary configurations is
by no means fully constrained, especially when we consider both ocean-covered
and land planets.

Water is a minimum requirement for life as we know it. Therefore, we are
naturally most interested in the climate of water worlds. Note our Earth is
effectively a water world even with $\sim$30 percent of its surface covered by
continents. For simplicity, in exoplanet modeling studies these are often
assumed to be completely ocean-covered worlds, but that is not necessarily the
case, especially for M-dwarf planets. Determining the amount of water expected
on a rocky M-dwarf planet in the habitable zone is no simple matter. Among
other things, it depends on whether we expect the planets to have formed in-
situ or have migrated. Due to the low luminosity of M dwarfs, their ice lines
are relatively close to the habitable zone ($\sim$ 0.3 AU), such that
migration in the inner disk could facilitate water-rich habitable zone planets
(Ogihara & Ida, 2009). On the other hand, in situ formation of volatile rich
planets would be less likely (Lissauer, 2007; Raymond et al., 2007, 2022),
even when accounting for volatile delivery through the migration of icy
planetesimals (Ciesla et al., 2015). There is also a higher risk of water loss
for M-dwarf planets due to early periods of increased stellar flux and
energetic flare activity (Bolmont et al., 2017; Tian & Ida, 2015; Luger &
Barnes, 2015; Ramirez & Kaltenegger, 2014), as well as periods of intense
tidal heating prior to reaching spin-orbit resonance (Barnes et al., 2013).
For the purpose of this study, we are specifically interested in moist M-dwarf
planets, but it is possible that these could be predominantly water-limited
land planets rather than ocean-covered worlds. Therefore, our study compares
both both water-rich ocean-covered scenarios, and water-limited land planet
scenarios.

Figure 1:  Stellar spectra for AD Leonis (red), HD22049, a K-star (orange),
and the Sun (yellow). Blue lines show water absorption values using the
Lorentz profile, from the HITRAN Database.

The presence of water on the surface of a planet can have a wide range of
effects on the climate, for example, altering the surface heat capacity
(Cronin & Emanuel, 2013; Donohoe et al., 2014). Water also impacts the surface
albedo, whether due to low albedo liquid water or the ice-albedo feedback.
Water plays a key role in the radiative budget, influencing cloud formation
and structure (Stevens, 2005), which in turn affects planetary albedo (Donohoe
& Battisti, 2011) and greenhouse effect intensity. It can also alter the
atmospheric energy transport (Held & Soden, 2006), enhance local energy
storage (Donohoe et al., 2014; Lobo & Bordoni, 2020), and functions as an
important greenhouse gas (Held & Soden, 2000) with a strong positive climate
feedback. Due to water’s various climate feedbacks and its effect on the
atmospheric structure, the habitable zone of a water-limited Earth twin is
broader than that of an aquaplanet Earth (Abe et al., 2011). But, while
water’s impact on climate is well-understood for Earth, many of these
fundamental climate feedbacks behave differently on M-dwarf planets due to the
lower frequency of the stellar radiation.

We typically refer to the stellar radiation as the “shortwave” flux
(abbreviated SW), to distinguish it from a planet’s “longwave” thermal
emissions (LW). However, compared to other stars M-dwarfs have relatively long
wavelength emissions. M dwarf effective temperatures range from 2000 to 3800K,
such that their peak emissions are at near-IR wavelength (Fig. 1) and a large
fraction of the SW overlaps with absorption bands for CO2 and H2O (Kasting et
al., 1993; Selsis et al., 2007). Water’s absorption is particularly important,
with multiple absorption bands near the spectra peak. For reference, we plot
H2O absorption in Fig. 1, using values from the HITRAN database. This implies
that unlike Earth’s atmosphere, which is predominantly heated from below due
to SW absorption at the surface, a moist M-dwarf planetary atmosphere can be
heated at various levels. It is also important to note that ice and snow
albedos are also lower for IR radiation (Dunkle & Bevans, 1956), such that the
ice-albedo feedback is weaker. On M-dwarf planets the effects of longer-
wavelength radiation result in warmer climates than would be obtained for
equal stellar flux from a higher frequency emitter, such as G and F dwarfs
(Shields et al., 2013, 2014, 2016). Therefore, the climate response to
increased stellar radiation on M-dwarf planets is not necessarily similar to
Earth’s, and climate scenarios must be explored in the context of specific
stellar types.

In this paper, we explore climate at the inner edge of the M-dwarf habitable
zone, to determine how fractional habitability changes as dayside temperatures
start to exceed habitable limits. In particular, we seek to determine whether
it is possible to sustain a habitable band at the transition between a
scorching dayside and a glacial nightside, a scenario which we refer to as
terminator habitability, or whether fractional habitability becomes impossible
once the substellar region surface temperatures surpass habitable limits.
While there is no broad consensus in terms of a habitable temperature range,
here we use a relatively narrow definition of temperatures between 0°C and
50°C, to provide a more conservative measure and to better track how
fractional habitability changes (Shields et al., 2016). Therefore, we are
exploring whether planets with a medium or large day-to-night temperature
gradient can exist specifically near the inner edge of the habitable zone,
such that the dayside temperatures surpass 50°C, and the resultant
consequences in terms of climate and fractional habitability.

For this study, we consider Earth-like planets orbiting AD Leonis, which has a
spectral classification of M3.5V, such that synchronously rotating planets in
its habitable zone are in a slowly rotating regime (Haqq-Misra et al., 2018;
Noda et al., 2017). On these planets, the rossby deformation radius is smaller
than the planetary radius, resulting in a dynamical regime with strong
convective activity at the substellar point, a thermally direct day-to-night
circulation, and tend to have roughly symmetric terminator properties. For
dimmer M-dwarf stars, the habitable zone can include synchronously rotating
planets that are more rapidly rotating, which have distinct dynamical
properties, including a stronger mean zonal flow that produces a band-like
structure at the equator (e.g. Haqq-Misra et al., 2018). But here we focus on
dynamics relevant to the habitable zone of brighter M dwarfs, which represent
the majority of observed M-dwarf stars (Perryman et al., 1997) and host the
majority of exoplanets detected so far (according to the NASA Exoplanet
Archive).

We use a 3D global climate model to determine whether it is possible to
sustain a temperature gradient large enough for a terminator habitability
scenario, and explore the implications of terminator habitability for future
climate characterization studies. It is not our goal to precisely quantify the
habitable zone edge, given that its location is dependent on a large range of
properties, including planetary radius and surface gravity (Yang et al.,
2019a; Thomson & Vallis, 2019), among many other factors (see e.g. Meadows &
Barnes, 2018; Kasting et al., 1993; Kopparapu et al., 2013; Yang et al., 2013,
2014a), but rather to explore the mechanisms through which the atmosphere
responds to increased stellar flux, including changes in the radiative budget
and atmospheric energy transport, to determine the viability of these surface
climate configurations. We begin with a study of water-abundant aquaplanet
simulations (Sections 3.1), followed by a comparison with water-limited land
planet simulations (Section 3.2). Note that for simplicity, we will define
habitability based only on surface temperature for the majority of the paper,
but we discuss long term water availability and the implications for
habitability in Section 3.3 and 3.4.

##  2 Methods

In this study, we use ExoCAM (Wolf, 2021; Wolf et al., 2022), which is a
modified version of the Community Atmosphere Model (CAM4), with correlated-k
radiative transfer (ExoRT)111 During revisions, a minor interpolation error
was found and fixed in the ExoCAM ExoRT package (commit 0223865, Wolf, 2021).
We ran simulations, branched from the previous spun-up cases, using the
updated model. These confirmed that the bug did not impact our results. , and
a finite-volume dynamical core. A complete description of the code and its
lineage is available in Wolf et al. (2022) and references within. In order to
isolate the effects of instellation on the day-night temperature contrasts and
terminator habitability thresholds, our simulations are Earth-like in terms of
radius and gravity, with an atmospheric composition of 40Pa CO2, 0.17Pa CH4,
and the remainder N2, totaling 1 bar for the dry component of the atmosphere,
with variable H2O also included. Simulations were run with a horizontal
resolution of $4\degree\times 5\degree$, and 40 vertical levels. We use the
stellar spectra of AD Leonis (Reid & Gizis, 1997; Segura et al., 2005), shown
in Fig. 1, prescribing zero eccentricity and zero obliquity, such that there
is no seasonal cycle. We also calculate the two-band snow and ice albedos,
which we refer to as the visible ($0.25<\lambda<0.75\mu m$) and near-IR
($0.75<\lambda<2.5\mu m$), weighted by the stellar spectrum:

Table 1: Albedo Values | Visible | Near-IR
---|---|---
Snow | 0.97 | 0.48
Ice | 0.64 | 0.17

All simulations use a time step of 1800s and are run until they reach steady
state. Results are shown as a time average of the last 10 years. Two
simulations did not reach steady state (Aq25 and Aq25h), and instead failed
due to numeric instabilities as they entered a runway greenhouse state.
Results for Aq25 are briefly discussed using the average of the last stable
month of model time. Results are shown on a planetary grid where the
substellar point is at 0° longitude and 0° latitude. Figures with a vertical
profile use the CESM hybrid sigma-pressure coordinate, plotted on a scale of
0-1 where 1 is the simulation’s reference surface pressure (Po).

The aquaplanet simulations (Aq) have uniform 50m mixed layer depth and no
ocean energy transport. The Aq34 aquaplanet is located 0.15au from the star.
We use this planet as a starting point for our analysis because, for a planet
orbiting AD Leonis, this orbital distance results in an Earth-like solar
constant and a mostly temperate dayside climate. To explore the inner edge of
the habitable zone and the transition into runaway greenhouse states, we also
use simulations with smaller orbital radii. Per Kepler’s third law, a
reduction in orbital radii is accompanied by a change in orbital period, such
that

$P=\Big{(}a^{3}\frac{M_{\bigodot}}{M_{\bigstar}}\Big{)}^{\frac{1}{2}},$ (1)

where $P$ is the orbital period in years, $M_{\bigstar}$ is stellar mass,
$M_{\bigodot}$ is the Sun’s mass, and $a$ is the orbital radius in au.
Therefore, to reduce orbital radii, and increase stellar flux, we compare
simulations with orbital and rotational periods of 34, 30, and 25 days (Table
2). We also include aquaplanet simulations with atmospheric surface pressure
reduced to 0.5bar (Aqh).

Given our interest in the greenhouse effect, and the enhanced importance of
water vapor as an absorber for redder stars, we also use land planet
simulations (L) where the planet’s water is limited. These simulations are
initialized without soil or ground water, with the ground properties of sand
everywhere (no oceans). All water is initially homogeneously distributed in
the atmosphere, where we prescribe the specific humidity (Q). This
configuration uses the Community Land Model (CLM) (Lawrence et al., 2010),
such that over time, water can remain in the atmosphere, precipitate,
accumulate in the soil, or solidify as snow/ice. We emphasize that these
water-limited simulations are significantly different from the aquaplanet
simulations, because aquaplanet simulations have unlimited water available for
evaporation at each grid box not covered in ice.

Table 2: List of Simulations  Name | Type |  $P_{\text{orb}}$ (days) | a (au) | F (W/m2) | Po
---|---|---|---|---|---
Aq34 | Aquaplanet | 34 days | 0.154 | 1378 | 1.0 bar
Aq34h | Aquaplanet | 34 days | 0.154 | 1378 | 0.5 bar
Aq30 | Aquaplanet | 30 days | 0.141 | 1629 | 1.0 bar
Aq30h | Aquaplanet | 30 days | 0.141 | 1629 | 0.5 bar
Aq25 | Aquaplanet | 25 days | 0.125 | 2077 | 1.0 bar
Aq25h | Aquaplanet | 25 days | 0.125 | 2077 | 0.5 bar
L34Qe3 | Q=0.0010 | 34 days | 0.154 | 1378 | 1.0 bar
L34Qe4 | Q=0.0001 | 34 days | 0.154 | 1378 | 1.0 bar
L25Qe4 | Q=0.0001 | 25 days | 0.125 | 2077 | 1.0 bar

##  3 Results

###  3.1 Water-Abundant Worlds

We begin with a simple comparison of the water-abundant cases (Aq34, Aq30, and
Aq25), examining the maximum and minimum surface temperatures. By definition,
in order for terminator habitability to occur a planet must sustain large day-
night temperature gradients. However, the simulations show that there is a
significant reduction in the surface temperature range as the planet moves
closer to the star (Fig. 2a) in agreement with the findings of Yang et al.
(2019b), Haqq-Misra et al. (2018), and Noda et al. (2017). Comparing Aq34 and
Aq30, we note that Aq30 is on average 25K warmer, but the day-night
temperature contrast ($\Delta$T) is only 34K (Table 3). This is unfavorable
for terminator habitability because, as the planet’s mean temperature
increases, the nightside warms more than the dayside, and the day-night
contrast becomes small. Therefore, at the inner edge of the habitable zone, we
would expect planets to either be mostly habitable (e.g. Fig. 3b) or, if their
temperatures become sufficiently high, entirely uninhabitable due to a global
runaway greenhouse effect. In this section, we explore the mechanisms that
lead to small day-night temperature differences in the water-abundant planets.

Figure 2:  Minimum (blue) and maximum (red) surface temperature (a) and precipitable water (b) for 1bar (circle) and 0.5bar (h) simulations. L markers show values for L34Qe4 and L25Qe4 simulations. The shaded region highlights the temperature range between 0 and 50°C, where we might expect climate favorable to complex life forms. Subplot c shows substellar (orange) and antistellar (purple) TOA OLR values. Table 3: Planetary surface temperatures and fractional habitability (FH).  | Mean T | T${}_{\text{max}}$ | T${}_{\text{min}}$ |  $\Delta$T | FH 0 - 50°C
---|---|---|---|---|---
Aq34 | 259K | 291K | 236K | 54K | 32%
Aq34h | 257K | 297K | 232K | 66K | 30%
Aq30 | 283K | 301K | 268K | 33K | 79%
Aq30h | 282K | 302K | 264K | 38K | 75%
L34Qe3 | 279K | 370K | 248K | 122K | 24%
L34Qe4 | 247K | 350K | 205K | 145K | 24%
L25Qe4 | 295K | 404K | 254K | 150K | 16%

Increasing the mean planetary temperature can result in many changes to the
radiative budget, including competing effects. For example, increasing water
vapor contributes to cloud formation and an increase in planetary albedo,
which has a strong effect on the climate of M-dwarf planets (Yang et al.,
2013; Yang & Abbot, 2014). As can be noted in Fig. 4c, regions with abundant
clouds reflect a larger fraction of incoming SW, which tends to reduce warming
particularly in the substellar region. Meanwhile, increased water vapor also
leads to stronger LW absorption and a stronger greenhouse effect. But,
individually, neither of these effects supports a reduction in the day-night
temperature contrast. While, the presence of clouds would reduce dayside
warming, the cloud fraction is already high in Aq34, such that there is not a
significant difference between the top-of-atmosphere (TOA) albedo in Aq34 and
Aq30. Meanwhile, the increased dayside greenhouse effect on its own, in the
absence of any energy and moisture transport, would actually promote increased
day-night temperature contrasts. Clearly that is not what occurs, as can be
noted in the surface temperatures (Fig. 2a) and precipitable water (Fig. 2b).
To fully understand the relationship between the day and nightside climate, we
have to explore the relationship between the radiative budget and the
atmospheric energy transport.

Figure 3:  Surface temperatures ($\degree$C), plotted with substellar point at
center. Black lines indicate the terminator.

As can be observed for Aq34, in Fig. 5 (e,j,o), water vapor (Q) maximizes near
the substellar region, and has large values even in the mid-troposphere. The
water vapor leads to peaks in SW absorption in the mid and upper atmospheric
layers (Fig. 5c,h), which when combined with reflection from clouds, prevents
most of the SW from reaching the surface (b,g,i). The atmospheric substellar
absorption of M-dwarf radiation on an Earth-like planet is stronger than for F
or G-dwarf instellation (Shields et al., 2019), such that in this simulation
less than 50 W/m2 of SW reaches the surface. Even with 400 W/m2 of LW, the
substellar surface is a radiative local minimum (Fig. 5,d). Near-surface
downwards radiative flux maxima occur in the “mid-latitudes” (roughly 50° from
the substellar point), creating regions of strong evaporation. These are also
regions of substantial sensible energy and moisture flux divergence due to the
large-scale overturning circulation’s near-surface winds patterns, such that
the “mid-latitude” temperatures are more moderate. Meanwhile, the resulting
transport of sensible energy and moisture helps sustain the substellar
region’s temperature maxima, despite the weak radiative values.

Figure 4:  a) Near surface ($\sigma=0.9$) temperatures (°C) and wind patterns
(quivers). b) Mid-troposphere ($\sigma=0.5$) geopotential height (km) and wind
patterns. c) TOA albedo and vertically integrated total cloud fraction (blue
contours), d) Surface evaporation rates and precipitation reaching the surface
in yellow contours ($10^{5}$ kg/m2s)

.

Figure 5:  Radiative profiles for Aq34. Left column shows temperature lapse
rate (K/km) calculated with area weighted average of region a: 0-30°, f:
30-60°, k: 60-90° from substellar point. The remaining columns, from left to
right, show downward SW (W/m2); solar heating rate (K/day); Total downward
radiation (SW + LW); Specific Humidity (kg/kg). The cross sections were taken
at longitudes 0°, 60°, and 90°.

With little SW reaching the surface, and most stellar radiation being absorbed
in the mid and upper atmosphere, we would expect weak atmospheric lapse rates
in the aquaplanet simulations. Fig. 5a shows atmospheric lapse rates for Aq34,
averaged within 30° radius of the substellar point. The near-surface lapse is
roughly adiabatic (8.8 K/km), and decreases above the boundary layer in the
region of strong convection, as would be expected for a region with moist
convection. In the region near $\sigma=0.5$, where we note strong SW
absorption (Fig. 5c), the lapse rate decreases further, dropping below 4 K/km
at $\sigma=0.4$. This is significantly lower than the local saturated
adiabatic lapse rate of $\sim$7 K/km, creating a region of convective
stability that limits the vertical flow (Fig. 6b).

Overall, the atmospheric circulation follows a structure typical for a planet
in a slowly rotating regime (e.g. Merlis & Schneider, 2010; Showman et al.,
2013; Haqq-Misra et al., 2018). There is strong near-surface convergence at
the substellar region (Fig. 4a), which drives an upward flow. The flow then
diverges in the mid/upper troposphere (Fig. 4b), and sinks in the nightside of
the planet. But, due to the region of convective stability, the overturning
circulation structure is shallow (6c) and above the region of convective
stability, a second smaller cell forms. The divisions between the cell layers
are indicated with black markers in Fig. 6c. If we considered only the cell’s
appearance, and the relatively small mass transport, we might erroneously
assume that atmospheric transport plays a reduced role in the water-abundant
simulations.

Figure 6:  Comparison of aquaplanet (Aq34) and land planet (L34Qe3 and L34Qe4)
simulations with orbital period of 34 days. Left column shows atmospheric
lapse rate, averaged between 0°and 30° from the substellar point. The middle
column shows vertical velocity (Pa/s), and the boundary layer depth (dashed
line). The right column shows mass stream functions (kg/s), and black markers
indicate levels where substellar horizontal flow divergence switches sign. The
velocity and stream function cross sections are averaged zonally from -15° to
15° latitude. Vertical velocity values for panel b are multiplied by a factor
of 2 to allow for easier comparison.

While the atmospheric circulation intensity reflects wind speeds and mass
transport, it does not necessarily provide a good indication of the net energy
transport. For a planet at steady state, in the absence of ocean energy
transport, any TOA imbalance must be equivalent to the vertically integrated
atmospheric energy transport, such that the energy budget (e.g. Neelin, 2007)
for a synchronously rotating planet without a seasonal cycle, can be
simplified to:

$R_{\text{toa}}-F_{\text{sfc}}=\nabla\cdot\langle\overline{\mathbf{v}h}\rangle,$
(2)

where $R_{\text{toa}}$ is top-of-atmosphere radiative fluxes, and
$F_{\text{sfc}}$ the surface fluxes, which include radiative, sensible, and
latent heat flux at the surface. The moist static energy term is defined such
that $h=c_{p}T+L_{v}Q+gZ$, comprised of dry enthalpy $c_{p}T$, latent energy
$L_{v}Q$, and potential energy $gZ$. $\langle\,\cdot\,\rangle$ denotes a
vertical integral, and $\overline{(\,\cdot\,)}$ denotes a temporal long-term
average.

As is the case with Earth’s mean meridional circulation, the near-surface
convergence of warm air (in this case at the substellar point) results in net
convergence of dry enthalpy (or sensible energy), and latent energy. The
higher altitude divergent flow (Fig. 4b), tends to transport colder and dryer
air, but with larger amounts of potential energy. On Aq34, both branches
transport a significant amount of sensible energy (Fig. 7), such that the
sensible energy transport partially cancels out when vertically integrated.
The latent energy transport is comparatively small, such that the potential
energy transport of the upper branch ultimately determines the sign of the
vertically integrated energy transport (Fig. 8a). The small value of the
latent energy flux divergence indicates that, while water clearly impacts the
overall climate outcome, the latent energy transport has only a small effect
on the temperature gradients. The energy flux divergence (positive values in
Fig. 7 and 8) indicates net energy transport away from the substellar point.

Figure 7: Energy flux divergence in substellar region, for Aq34 (solid lines)
and L34Qe4 (dashed lines). Energy fluxes are broken down into sensible (red),
latent (blue), and potential energy (green) components, averaged across area
within 30° latitude of the substellar point, and weighted by vertical density
profile to facilitate comparison of upper and lower level energy transport.
Markers show the substellar boundary layer height for Aq34 (blue) and L34Qe4
(orange).

As we increase the stellar flux (Aq30), the increased dayside SW ($\Delta
SW_{down}$), is slightly reduced due to reflection from the TOA albedo
($\Delta SW_{up}$), and balanced partially by increased LW ($\Delta LW$), and
mostly by increased energy transport divergence ($\Delta E$). This can be
noted in Fig. 8b, which compares Aq30 to Aq34. The net energy transport
remains similar in shape (Fig. 8a,c), but the net sensible energy transport is
slightly reduced in favor of latent and potential energy transport.

Many of the trends noted for Aq30 persist for Aq25. Simulations that do not
achieve steady state, as is the case for Aq25, must be interpreted with
caution. Prior to reaching numerical instability, day-to-night temperature
differences become nearly negligible (Fig. 2a). The energy budget is not
properly closed, as we might expect for a simulation not at steady state. But
we notice sensible energy transport becomes small, likely due to the reduced
temperature gradients, and the majority of the transport is done by potential
and latent energy transport, which increase to similar magnitudes.

In summary, we can attribute the reduction of day-night temperature
differences in our water abundant planets to an increase in net energy
transport, which diverges additional energy away from the dayside, and
converges additional energy in the nightside. In other words, for increased
stellar radiation, the transport increases, reducing the intensity of dayside
warming, and enhancing nightside warming. With higher surface and atmospheric
temperatures, as well as increased moisture flux convergence, nightside
atmospheric water vapor increases (Fig. 2b), allowing for additional nightside
warming due to an increase in the local greenhouse effect.

As previously mentioned, the atmospheric energy transport is responsible for
resolving radiative imbalances (Eq. 2), such that these changes in the
atmospheric transport can also be more simply interpreted in the context of
the top-of-atmosphere radiative budget. As described in Koll & Abbot (2016),
for planets in a weak-temperature-gradient regime, the circulation can be
thought of as a heat engine driven by dayside heating, and cooling due to
global thermal emissions to space. The aquaplanet simulations have high
atmospheric water vapor content, resulting in optically thick dayside
atmospheres, such that their dayside thermal emissions are determined by the
upper tropospheric temperature. In such cases, increasing the stellar flux has
a small impact on the outgoing longwave radiation (OLR) (Yang & Abbot, 2014),
as can be noted in Fig. 2c and Fig. 8b. Meanwhile, the colder nightside has
low humidity and a relatively transparent atmosphere (except for Aq25), such
that warming is closely tied to increased OLR. Therefore, in a moist and
opaque atmosphere, increasing the stellar radiation as we go from Aq34 to Aq30
results in a small dayside OLR increase, a larger nightside OLR increase, and
requires a corresponding increase in day-to-night energy transport. The extra
SW absorbed in Aq30 is primarily balanced by a larger nightside OLR, which is
only achievable through a strengthening of the atmospheric energy transport
and an overall reduction in day-to-night temperature gradients.

Given that atmospheric energy transport is essential for reducing the day-to-
night temperature contrast, we might expect to see larger gradients on planets
with thinner atmospheres. However, that is not observed in our aquaplanet
simulations where we half the background N2 values. As can be noted in Fig. 2a
and Fig. 3, the temperature range remains almost unchanged when the surface
atmospheric pressure is halved. Instead, we note overall increased wind
speeds, such that the energy transport terms remain roughly similar. These
results are in general agreement with the findings from Zhang & Yang (2020),
which showed that background gas pressure has a relatively small impact on the
inner edge of the habitable zone for synchronously rotating Earth-like
planets. They show that the effects are non-monotonic, and significantly less
important than changes in the cloud scheme (Bin et al., 2018) and rotation
rate (Kopparapu et al., 2017). Therefore, we would not expect that pressure
alone would determine the viability of terminator habitability for the Earth-
like planets we are considering.

Figure 8:  Breakdown of the planetary energy budget. The left column shows
changes in the global mean radiative budget relative to Aq34, such that
$\Delta X=X-X_{Aq34}$, where X is incoming stellar radiation $SW_{down}$
(yellow), TOA reflected SW ($SW_{up}$, orange), TOA outgoing LW (maroon), the
net TOA radiation ($R_{TOA}$, in light blue), and the net energy flux
divergence ($E=\nabla\cdot\langle\overline{\mathbf{v}h}\rangle$, dashed black
line). For each simulation, the actual energy flux divergence values are
plotted on the right, broken down into sensible (red), latent (blue), and
potential energy (green) components, vertically integrated and meridionally
averaged. The thin maroon line shows TOA upward LW values. We also plot net
atmospheric energy input ($R_{\text{toa}}-F_{\text{sfc}}$), as a dashed light
blue line, which overlaps with the net energy flux divergence, confirming a
closed budget.

###  3.2 Water-Limited Planets

Here we switch our focus to water-limited land planets, using ExoCAM with a
land planet setup, initialized with varying amounts of atmospheric water vapor
(L34Qe3 and L34Qe4). Once temperatures equilibrate, these simulations have 20%
and 1% of Earth’s atmospheric precipitable water respectively, with small
amounts of water also in the form of surface ice and soil water, concentrated
on the nightside. We will begin by considering the differences between Aq34
and L34Qe4.

With reduced atmospheric water vapor, we expect a reduction in cloud coverage,
especially for low-level clouds. This expectation is confirmed, and results in
a reduction of the dayside albedo, which becomes largely homogeneous with
values between 0.2 and 0.26 in the water-limited simulations. There is still
some condensation, precipitation, and high altitude cloud formation, occurring
near the substellar point. But, with re-evaporation in the atmosphere, surface
precipitation becomes negligible. These precipitation patterns, noted in all
of our land planet simulations, appear to follow the “transition regime”
described in Ding & Wordsworth (2021). Despite having the same stellar flux as
Aq34, L34Qe4 has a larger net SW due to reduced cloud coverage and lower
albedo (negative $\Delta SW_{up}$ in Fig. 8d). The resulting increase in
substellar net SW must result in one of two outcomes: either increased TOA LW
emissions, or increased day-to-night energy transport.

Near the substellar point, a much larger fraction of SW reaches the surface
(Fig. 9), with dayside surface temperatures reaching maxima of 355K (L34Qe4).
Though there is still water-vapor SW absorption, the atmosphere’s vertical
temperature structure more closely resembles a dry region on Earth, with
dayside lapse rates remaining steeper than 6K/km until $\sigma=0.25$ (Fig. 6).
Without the low region of convective stability noted in Aq34, the water-
limited simulations have deeper overturning cells, ascending at the substellar
point (Fig. 6,f,i).

Examining the energy flux divergence (Fig. 8e), we note a strengthening of the
individual energy transport terms, but a weakening of the net transport. The
deepening of the atmospheric circulation, combined with strong high altitude
winds, leads to significantly larger values of potential energy transport away
from the “eye”, as can be noted in Fig. 7 which shows the vertical profile of
the energy flux divergence averaged across an area within 30° of the
substellar point. Sensible energy fluxes are intensified in both the upper and
lower branches of the circulation, but the lower branch increases
significantly in both intensity and depth. The boundary layer in the water-
limited simulation is significantly deeper (Fig. 7), and the lower branch of
the circulation occupies its full depth, reaching up to $\sigma=0.4$. This
effect can also be noted in the upward shift of the cell’s core (Fig. 6). This
“bottom-heavy” cell structure achieves higher values of sensible heat
conversion at the substellar point, which largely cancels out the increased
potential energy divergence, such that the net day-to-night energy transport
is reduced relative to Aq34 despite the increased net SW. Comparing $\Delta
SW_{up}$ to the change in net energy flux divergence ($\Delta E$), we note
that the effects of the reduced atmospheric energy transport are of comparable
intensity to the effects of the reduced albedo.

Figure 9:  Same as top row of Fig. 5, with each row showing lapse rates
averaged over a region 0-30° from the substellar point and cross sections
taken at 0° longitude, for Aq30, L34Qe3, and L34Qe4. To facilitate comparison
across simulations, specific humidity is shown on a logarithmic scale.

With reduced albedo and reduced energy transport, the large TOA imbalance is
resolved by an increase in dayside TOA LW emissions (positive $\Delta LW$ Fig.
8d). It is interesting to note that in the water-limited cases, the substellar
region becomes a region of LW maxima, rather than a minima as was the case for
the aquaplanet simulations (Aq34 and Aq30). Meanwhile, the reduced nightside
energy convergence results in a reduction of nightside temperatures. Together,
these effects allow the water-limited case to achieve equilibrium with larger
day-night temperature gradients, spanning 139K (L34Qe4) and 122K (L34Qe3). If
we assume that surface temperatures surpassing 50°C are hostile to life, and
temperatures nearing 100°C would be uninhabitable, both L34Qe3 and L34Qe4
water-limited scenarios can be considered cases of terminator habitability.

As with the aquaplanets, we can re-frame these results in terms of the
atmospheric radiative imbalances to obtain a simpler intuition. On these land
planets, the atmosphere is dryer, and therefore less cloudy and opaque. This
allows the dayside to be closer to radiative equilibrium (Fig. 8d), and
therefore implies a reduced energy transport to the nightside. The water-
limited L34Qe4 absorbs a greater fraction of the incoming stellar radiation at
the surface than Aq34, achieving higher surface temperatures, and higher OLR,
thus requiring a weaker energy transport (F ig. 8d). In turn, the nightside
receives less energy, and is colder than the aquaplanet counterpart. This
results in significantly larger day-night temperature contrasts, such that a
temperate terminator becomes not only easily achievable, but also a likely
scenario for land planets in the habitable zone.

The simulations discussed so far illustrate the possibility of terminator
habitability as an alternate climate configuration within the habitable zone
(at orbital period of 34 days). However, it is also worth emphasizing that
water-limited planets can endure higher stellar fluxes without entering a
runaway greenhouse state. For example, while Aq25 and Aq25h become unstable
due to runaway warming, an equivalent water-limited simulation (L25Qe4)
achieves steady state and sustains a temperate terminator climate (Fig. 3f).
The increase in substellar SW relative to L34Qe4 is balanced in roughly equal
parts by an increase in LW and an increase in local energy flux divergence.
The resulting climate is on average warmer than L34Qe4, and the day-to-night
temperature range increases to 154K (Fig. 2). The L25Qe4 nightside
temperatures remain just below freezing, dayside temperatures reach scorching
highs over 400K, and the terminator remains within the 270-320K range. While
the fractional habitability (based on surface temperatures) of L25Qe4 may be
smaller than Aq34 and Aq30 (Table 3), the contrasting results of L25Qe4 and
Aq25 illustrate how we might expect to observe planets with diverse surface
climates and regions of temperate surface temperatures even at the inner edge
of the habitable zone.

###  3.3 Long-Term Stability

The water-limited simulations show that terminator habitability is possible.
However, based on the results discussed so far, we have not yet established
whether they are ideal observational target in the search for life beyond
Earth. After all, these simulations depict planets with harsh climates,
limited regions of habitability, and limited water availability, which may
seem unappealing when compared to the climates of ocean worlds obtained with
aquaplanet simulations. However, arid exoplanets offer observational
advantages, because their reduced cloud coverage could facilitate detection of
water (Ding & Wordsworth, 2022) by JWST, potentially making them more
practical near-future targets. Also, a closer look at the aquaplanet moisture
budget shows that we must be cautious of assuming long-term climate stability,
which is crucial for increasing the odds of developing life.

Aquaplanet simulations provide the atmosphere with an unlimited water source.
On the planet’s dayside, not only is there water available everywhere, but
there is no depletion in regions of negative net precipitation, nor resistance
to evaporation as would occur from drying soil (e.g. van de Griend & Owe,
1994). For Aq34, which has dayside temperatures comparable to Earth’s, this
results in an Earth-like amount of atmospheric precipitable water ($\sim$70%
of Earth’s atmospheric water). The “eye” is a region of negative net
precipitation, and there is continuous transport of water vapor to the
nightside, resulting in a total of 1.5 x 1016 kg/s of snowfall. For scale,
this rate of snowfall would cold-trap a volume of water equivalent to Earth’s
oceans in 90 thousand years. The presence of a deep global ocean circulation
could, in some cases, help deglaciate part or all of the nightside ocean (Hu &
Yang, 2014). But in the absence of a global ocean with favorable properties,
such as large depth (Hu & Yang, 2014) and high salinity (Olson et al., 2020),
we can assume a limited return flow for water deposited on the nightside. Even
accounting for ice sheet dynamics, which would imply a down-gradient flow of
ice towards the dayside, it is possible that over time a majority of the water
from an ocean-covered Earth-like planet could become cold-trapped (Menou,
2013; Leconte et al., 2013).

While Aq34 could be at risk from cold-trapping, the risk decreases
significantly for planets closer to the star. The Aq30 nightside temperatures
are higher, such that the equatorial region of the nightside remains above
freezing. The snowfall is confined to the higher latitudes, and the global
snowfall rate decreases by a factor of 100 relative to Aq34. However, the
atmospheric precipitable water increases to approximately 3.5 times Earth’s
atmosphere. In particular, we note an increase in water vapor content above
the atmospheric cold trap, reaching values of 0.003 kg/kg ($\sim$0.005
mol/mol), which is two orders of magnitude larger than Aq34 values, and
reaches the limit proposed by Kasting et al. (1984) and confirmed in a series
of planetary studies (e.g. Kasting et al., 2015; Wolf & Toon, 2015). This
raises the possibility that such a planet would be vulnerable to continuous
water vapor loss, especially during periods of strong flare activity. These
losses are in addition to water losses resulting from enhanced heating during
the M-dwarf’s protracted pre-main sequence phase, which could impact all of
these scenarios and cause the loss of large water reservoirs (e.g. Bolmont et
al., 2017; Luger & Barnes, 2015). Therefore, we have abundant reason to
question water-rich M-dwarf scenarios.

Meanwhile for land planets L34Qe3, L34Qe4, and even L25Qe4, while there is
still risk of water loss in their early history, high altitude moisture values
are comparable to Aq34 and well below the Kasting et al. (1984) limit.
Similarly, despite their extremely low nightside temperatures, at equilibrium
there is no snow being deposited on the nightside. There is some snow and ice
accumulation in L34Qe3 while the model spins up, which remain trapped. But the
remaining moisture is subsequently retained in the atmosphere above the
planetary boundary layer, primarily on the dayside. We have confirmed that
simulations initialized with more water (e.g. Q = 0.01) have similar night-
side water trapping, and slowly tend toward the climate in L34Qe3, but would
require hundreds more simulated years to fully reach equilibrium.

Based on our land planet simulations, water-limited atmospheres appear to
reach a more stable configuration than their aquaplanet counterparts, being
less vulnerable to additional water vapor loss or nightside cold-trapping. It
also suggests that ocean world configurations subject to significant water
loss, without entering a runaway greenhouse state, could eventually reach a
stable configuration with terminator habitability. These considerations,
combined with the observational advantages of water-limited planets (Ding &
Wordsworth, 2022), suggests that despite their potentially limited or
localized fractional habitability, land planets will be important
observational targets in the coming years, and will play a prominent role in
the early stages of exoplanet climate characterization.

###  3.4 Terminator Water Availability

Thus far we have defined terminator habitability in terms of surface
temperature. However, water is also a crucial ingredient for life as we know
it, and could pose a problem for water-limited scenarios. In our simulations,
rain is strongly concentrated in the substellar region (for both aquaplanet
and land-planet simulations), and snow occurs across the nightside in some Aq
simulations. For the rotational configurations considered here, the terminator
never receives significant amounts of precipitation. Therefore, even though
terminator evaporation rates tend to be low, in our simulations it is a region
of weak negative net precipitation. In other words, based on the atmospheric
moisture budget, the terminator is a region that would tend to have an arid
climate. However, water availability on land planets is not exclusively
determined by the atmospheric moisture budget, but also depends on factors
such as groundwater transport and glacier behavior. These effects are in turn
dependent on many poorly constrained planetary properties such as orography,
soil structure, and geothermal heat flux. While we cannot quantify these
effects with our climate model, we can discuss them qualitatively to explore
how they would impact terminator water availability.

The role of glacier flow has been widely discussed both in the context of
snowball Earth (e.g. Goodman & Pierrehumbert, 2003; Pollard & Kasting, 2005)
and as a potential mechanism for releasing cold-trapped water on synchronously
rotating planets (e.g. Leconte et al., 2013). Ice sheets can deform under
their own weight, resulting in an ice flow down the ice thickness gradient. On
a synchronously rotating planet, this flow could push ice towards the
terminator, where the temperate climate would result in melting. The flow
velocity depends on, among other things, the ice geometry, planetary gravity,
and geothermal heat flux. For Super-Earths, Yang et al. (2014b) argue it could
largely prevent nightside cold-trapping on water abundant planets.

Generally, the ice flow would be slower on Earth-sized planets, especially
those with more moderate geothermal heat fluxes, but it could still provide a
large enough source of melt-water to create a moist climate in the terminator
region. For example, on Greenland and Antarctica flow speeds and discharge
into the ocean can reach rates over 1km/year (Rignot et al., 2011; Rignot &
Mouginot, 2012). Currently, model representations of glacier dynamics are
limited, and not fully suitable for planetary purposes. For CESM, even a
configuration with an evolving ice sheet, where glacier area and elevation are
adjusted so as to conserve mass and energy, would not capture the dynamics
relevant for the ice flow we consider here. But, if we make a conservative
estimate, considering only the slower flow rates observed at the Antarctic
glacier edges ($\sim 100$m/yr), and assuming only a shallow ice layer reaches
the terminator ($\leq$ 10m), that would imply water transport rates of order
10-2 kg/s. This exceeds the terminator evaporation rates in the Aq30 and Aq34
simulations (see Fig. 4d) by two orders of magnitude. Given that aquaplanet
simulations have infinite water availability, we can cautiously treat their
evaporation rates as an upper limit, such that the comparison of glacier flow
to Aq terminator evaporation rates indicates that glacier flow could sustain
surface water near the terminator. Depending on the exact rates, this could
result in swamp-like regions, oceans, or in smaller lakes and rivers in
regions of topographic lows. Even if flow rates were significantly lower than
expected and surface runoff resulted in a large evaporating area, soil
evaporative resistance would tend to retain some local moisture that could be
beneficial for the development of life in the terminator region.

Of course, not all planets would have nightside glaciers. The presence of
large ice deposits presumes a planet with abundant initial water, that was
either retained throughout the early period of intense stellar emissions, or
delivered later on. It also presumes that ice would be stable on the
nightside. For planets receiving sufficiently high instellation, the nightside
ice can sublimate (Ding & Wordsworth, 2020), such that water would likely only
be present in vapor form. Our arid land planet simulations are not ideal for
exploring this limit, but based on the values from Ding & Wordsworth (2020),
we would expect nightside sublimation to be important for planets in L25Qe4’s
orbit. Therefore, L25 planets would be more likely to have a temperate but dry
terminator than their L30 or L34 counterparts.

A less recognized, but also potentially important mechanism to consider for
planets with stable surface water, is groundwater transport. Faulk et al.
(2020) showed that surface and subsurface methane flows on Titan could help
shape the global hydrological cycle, transporting liquid from the low-latitude
highlands to the lower polar regions, resulting in a moist polar climate. For
synchronously rotating planets, Turbet et al. (2016) argues alignment of
large-scale gravitational anomalies and the star-planet axis is favored
(Wieczorek, 2007), such that the largest topographic basins would tend to be
near the substellar point or anti-stellar point, which has served as
motivation for studies of sustellar continents on ocean worlds (Lewis et al.,
2018). For a land planet, where precipitation is strongly concentrated at the
substellar region, a substellar topographic high could result in a
hydrological cycle that continuously transports water across a larger portion
of the dayside. Whereas a topographic low could further concentrate water near
the substellar point. Therefore, additional exploration of topography, runoff
and subsurface transport is necessary to constrain surface water availability
for a land planet, especially for those with significant large-scale
topography.

Based on the above discussion of water availability, we might expect that
some, but not all planets with temperate terminators would also have moist
terminator climates. We expect that future work, exploring a wider range of
surface configurations will help better constrain water availability for land
planets. Model advancements to include glacier dynamics into the water budget
would further improve our understanding of these water worlds. Finally, an
important caveat to this work is the fact that we have only considered planets
with an Earth-like atmospheric composition. Varying atmospheric mass and
composition would certainly impact these results. For example, Ding &
Wordsworth (2020) showed that nightside cold-trapping might be avoided on
planets with high CO2 values, due to increased nightside temperatures.
Constraining exoplanet CO2 levels remains a challenge, due to their dependence
on outgassing and weathering rates (Walker et al., 1981). Therefore, future
work may wish to consider both terminator habitability and high CO2 “eye”
habitability scenarios, while exploring additional factors such as the impact
of groundwater flow and long term climate sensitivity in the presence of
M-dwarf flare activity to quantify prospects for sustained habitability.

##  4 Conclusions

In this study we explored the possibility of terminator habitability, defined
as the existence of a habitable band at the transition between a scorching
dayside and a glacial nightside on a synchronously rotating Earth-like planet
orbiting an M dwarf. In particular, we explore the viability of this
configuration at the inner edge of the habitable zone to determine if it is
possible for dayside temperatures to exceed typical habitable temperatures,
approaching the runaway greenhouse limit, without leading to a planet-wide
runaway greenhouse state.

We find that a temperate terminator is not achievable with aquaplanet
simulations that seek to reproduce ocean-covered planets, but can easily occur
on water-limited land planets. On aquaplanets, increasing stellar flux leads
to reduced day-night temperature gradients, such that the planet would tend to
reach a homogeneous climate before the dayside reaches a runaway greenhouse
state, never passing through a terminator habitability state. Whereas on
water-limited planets, we find that large day-night temperature gradients are
easily achievable without entering a runaway greenhouse state. We also find
that the water-limited land planet configurations may be favorable in terms of
long term climate stability, with reduced risk of nightside water cold-
trapping or water vapor escape. We expect that water-limited terminator
habitability scenarios could be a stable configuration for ocean-covered
worlds after significant water loss, especially in the case of nightside cold-
trapping.

There are still many uncertainties regarding the water content of habitable
zone M-dwarf planets. Based on our current understanding, it is possible that
water-limited planets could be abundant and possibly more common than ocean-
covered worlds. Therefore, terminator habitability may represent a significant
fraction of habitable M-dwarf planets. Compared to the temperate climates
obtained with aquaplanets, terminator habitability does offer reduced
fractional habitability. Also, while achieving a temperate terminator is
relatively easy on water-limited planets, constraining water availability at
the terminator remains a challenge. Overall, the lack of abundant surface
water in these simulations could pose a challenge for life to arise under
these conditions, but mechanisms including glacier flow could allow for
sufficient surface water accumulation to sustain locally moist and temperate
climates at or near the terminator. We expect that future studies exploring a
broader range of land planet configurations, in particular those using future
generations of surface and ice models, will find a wide range of habitable
terminator scenarios in regimes intermediate to the water-limited and
aquaplanet cases considered here.

##  5 Acknowledgments

This material is based upon work supported by the National Science Foundation
under Award 1753373, and by a Clare Boothe Luce Professorship supported by the
Henry Luce Foundation. This research was also performed as part of the NASA’s
Virtual Planetary Laboratory, supported by the National Aeronautics and Space
Administration through the NASA Astrobiology Institute under solicitation
NNH12ZDA002C and Cooperative Agreement Number NNA13AA93A, and by the NASA
Astrobiology Program under grant 80NSSC18K0829 as part of the Nexus for
Exoplanet System Science (NExSS) research coordination network. We would like
to acknowledge high-performance computing support from Cheyenne
(doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information
Systems Laboratory, sponsored by the National Science Foundation.

## References

  * Abe et al. (2011) Abe, Y., Abe-Ouchi, A., Sleep, N. H., & Zahnle, K. J. 2011, Astrobiology, 11, 443, doi: 10.1089/ast.2010.0545
  * Anglada-Escudé et al. (2016) Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437, doi: 10.1038/nature19106
  * Armitage (2009) Armitage, P. J. 2009, Astrophysics of Planet Formation (Cambridge University Press), doi: 10.1017/CBO9780511802225
  * Barnes (2017) Barnes, R. 2017, Celestial Mechanics and Dynamical Astronomy, 129, 509, doi: 10.1007/s10569-017-9783-7
  * Barnes et al. (2013) Barnes, R., Mullins, K., Goldblatt, C., et al. 2013, Astrobiology, 13, 225, doi: 10.1089/ast.2012.0851
  * Bin et al. (2018) Bin, J., Tian, F., & Liu, L. 2018, Earth and Planetary Science Letters, 492, 121, doi: 10.1016/j.epsl.2018.04.003
  * Bochanski et al. (2010) Bochanski, J. J., Hawley, S. L., Covey, K. R., et al. 2010, The Astronomical Journal, 139, 2679, doi: 10.1088/0004-6256/139/6/2679
  * Bolmont et al. (2017) Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, Monthly Notices of the Royal Astronomical Society, 464, 3728, doi: 10.1093/mnras/stw2578
  * Ciesla et al. (2015) Ciesla, F. J., Mulders, G. D., Pascucci, I., & Apai, D. 2015, The Astrophysical Journal, 804, 9, doi: 10.1088/0004-637X/804/1/9
  * Cronin & Emanuel (2013) Cronin, T. W., & Emanuel, K. A. 2013, Journal of Advances in Modeling Earth Systems, 5, 843, doi: 10.1002/jame.20049
  * Ding & Wordsworth (2020) Ding, F., & Wordsworth, R. D. 2020, The Astrophysical Journal, 891, L18, doi: 10.3847/2041-8213/ab77d1
  * Ding & Wordsworth (2021) —. 2021, The Planetary Science Journal, 2, 201, doi: 10.3847/PSJ/ac2236
  * Ding & Wordsworth (2022) —. 2022, The Astrophysical Journal Letters, 925, L8, doi: 10.3847/2041-8213/ac4a5d
  * Donohoe & Battisti (2011) Donohoe, A., & Battisti, D. S. 2011, Journal of Climate, 24, 4402, doi: 10.1175/2011JCLI3946.1
  * Donohoe et al. (2014) Donohoe, A., Frierson, D. M. W., & Battisti, D. S. 2014, Climate Dynamics, 43, 1041, doi: 10.1007/s00382-013-1843-4
  * Dunkle & Bevans (1956) Dunkle, R. V., & Bevans, J. T. 1956, Journal of Atmospheric Sciences, 13, 212 , doi: 10.1175/1520-0469(1956)013<0212:AAAOTS>2.0.CO;2
  * Edson et al. (2011) Edson, A., Lee, S., Bannon, P., Kasting, J. F., & Pollard, D. 2011, Icarus, 212, 1, doi: 10.1016/j.icarus.2010.11.023
  * Fauchez et al. (2019) Fauchez, T. J., Turbet, M., Villanueva, G. L., et al. 2019, The Astrophysical Journal, 887, 194, doi: 10.3847/1538-4357/ab5862
  * Faulk et al. (2020) Faulk, S. P., Lora, J. M., Mitchell, J. L., & Milly, P. C. D. 2020, Nature Astronomy, 4, 390, doi: 10.1038/s41550-019-0963-0
  * Genio et al. (2019) Genio, A. D. D., Way, M. J., Amundsen, D. S., et al. 2019, Astrobiology, 19, 99, doi: 10.1089/ast.2017.1760
  * Gilbert et al. (2020) Gilbert, E. A., Barclay, T., Schlieder, J. E., et al. 2020, The Astronomical Journal, 160, 116, doi: 10.3847/1538-3881/aba4b2
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456, doi: 10.1038/nature21360
  * Goodman & Pierrehumbert (2003) Goodman, J. C., & Pierrehumbert, R. T. 2003, Journal of Geophysical Research, 108, 3308, doi: 10.1029/2002JC001471
  * Haqq-Misra et al. (2018) Haqq-Misra, J., Wolf, E. T., Joshi, M., Zhang, X., & Kopparapu, R. K. 2018, The Astrophysical Journal, 852, 67, doi: 10.3847/1538-4357/aa9f1f
  * Held & Soden (2000) Held, I. M., & Soden, B. J. 2000, Annual Review of Energy and the Environment, 25, 441, doi: 10.1146/annurev.energy.25.1.441
  * Held & Soden (2006) —. 2006, Journal of Climate, 19, 5686 , doi: 10.1175/JCLI3990.1
  * Hu & Yang (2014) Hu, Y., & Yang, J. 2014, Proceedings of the National Academy of Sciences, 111, 629, doi: 10.1073/pnas.1315215111
  * Ingersoll (1969) Ingersoll, A. P. 1969, Journal of Atmospheric Sciences, 26, 1191 , doi: 10.1175/1520-0469(1969)026<1191:TRGAHO>2.0.CO;2
  * Joshi (2003) Joshi, M. 2003, Astrobiology, 3, 415, doi: 10.1089/153110703769016488
  * Joshi et al. (1997) Joshi, M., Haberle, R., & Reynolds, R. 1997, Icarus, 129, 450, doi: 10.1006/icar.1997.5793
  * Kasting et al. (2015) Kasting, J. F., Chen, H., & Kopparapu, R. K. 2015, The Astrophysical Journal, 813, L3, doi: 10.1088/2041-8205/813/1/L3
  * Kasting et al. (1984) Kasting, J. F., Pollack, J. B., & Ackerman, T. P. 1984, Icarus, 57, 335, doi: 10.1016/0019-1035(84)90122-2
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108, doi: 10.1006/icar.1993.1010
  * Koll & Abbot (2016) Koll, D. D. B., & Abbot, D. S. 2016, The Astrophysical Journal, 825, 99, doi: 10.3847/0004-637x/825/2/99
  * Komacek & Abbot (2019) Komacek, T. D., & Abbot, D. S. 2019, The Astrophysical Journal, 871, 245, doi: 10.3847/1538-4357/aafb33
  * Kopparapu et al. (2017) Kopparapu, R. K., Wolf, E. T., Arney, G., et al. 2017, The Astrophysical Journal, 845, 5, doi: 10.3847/1538-4357/aa7cf9
  * Kopparapu et al. (2016) Kopparapu, R. K., Wolf, E. T., Haqq-Misra, J., et al. 2016, The Astrophysical Journal, 819, 84, doi: 10.3847/0004-637X/819/1/84
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, The Astrophysical Journal, 765, 131, doi: 10.1088/0004-637X/765/2/131
  * Kreidberg & Loeb (2016) Kreidberg, L., & Loeb, A. 2016, The Astrophysical Journal, 832, L12, doi: 10.3847/2041-8205/832/1/L12
  * Lawrence et al. (2010) Lawrence, D. M., Bonan, G. B., Flanner, M. G., et al. 2010, NCAR Technical Note 478
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., et al. 2013, Astronomy & Astrophysics, 554, A69, doi: 10.1051/0004-6361/201321042
  * Lewis et al. (2018) Lewis, N. T., Lambert, F. H., Boutle, I. A., et al. 2018, The Astrophysical Journal, 854, 171, doi: 10.3847/1538-4357/aaad0a
  * Lissauer (2007) Lissauer, J. J. 2007, The Astrophysical Journal, 660, L149, doi: 10.1086/518121
  * Lobo & Bordoni (2020) Lobo, A. H., & Bordoni, S. 2020, Icarus, 340, doi: 10.1016/j.icarus.2019.113592
  * Luger & Barnes (2015) Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119, doi: 10.1089/ast.2014.1231
  * Meadows & Barnes (2018) Meadows, V. S., & Barnes, R. K. 2018, Handbook of Exoplanets, 2771, doi: 10.1007/978-3-319-55333-7_57
  * Menou (2013) Menou, K. 2013, The Astrophysical Journal, 774, 51, doi: 10.1088/0004-637X/774/1/51
  * Merlis & Schneider (2010) Merlis, T. M., & Schneider, T. 2010, Journal of Advances in Modeling Earth Systems, 2, n/a, doi: 10.3894/JAMES.2010.2.13
  * Morley et al. (2017) Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, The Astrophysical Journal, 850, 121, doi: 10.3847/1538-4357/aa927b
  * Mulders et al. (2015) Mulders, G. D., Pascucci, I., & Apai, D. 2015, The Astrophysical Journal, 814, 130, doi: 10.1088/0004-637X/814/2/130
  * Neelin (2007) Neelin, J. D. 2007, The Global Circulation of the Atmosphere, 267
  * Noda et al. (2017) Noda, S., Ishiwatari, M., Nakajima, K., et al. 2017, Icarus, 282, 1, doi: 10.1016/j.icarus.2016.09.004
  * Nutzman & Charbonneau (2008) Nutzman, P., & Charbonneau, D. 2008, Publications of the Astronomical Society of the Pacific, 120, 317, doi: 10.1086/533420
  * Ogihara & Ida (2009) Ogihara, M., & Ida, S. 2009, The Astrophysical Journal, 699, 824, doi: 10.1088/0004-637X/699/1/824
  * Olson et al. (2020) Olson, S. L., Jansen, M., & Abbot, D. S. 2020, The Astrophysical Journal, 895, 19, doi: 10.3847/1538-4357/ab88c9
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49
  * Pollard & Kasting (2005) Pollard, D., & Kasting, J. F. 2005, Journal of Geophysical Research: Oceans, 110, doi: https://doi.org/10.1029/2004JC002525
  * Ramirez & Kaltenegger (2014) Ramirez, R. M., & Kaltenegger, L. 2014, The Astrophysical Journal, 797, L25, doi: 10.1088/2041-8205/797/2/L25
  * Raymond et al. (2007) Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, The Astrophysical Journal, 669, 606, doi: 10.1086/521587
  * Raymond et al. (2022) Raymond, S. N., Izidoro, A., Bolmont, E., et al. 2022, Nature Astronomy, 6, 80, doi: 10.1038/s41550-021-01518-6
  * Reid & Gizis (1997) Reid, I. N., & Gizis, J. E. 1997, The Astronomical Journal, 113, 2246, doi: 10.1086/118436
  * Rignot & Mouginot (2012) Rignot, E., & Mouginot, J. 2012, Geophysical Research Letters, 39, n/a, doi: 10.1029/2012GL051634
  * Rignot et al. (2011) Rignot, E., Mouginot, J., & Scheuchl, B. 2011, Science, 333, 1427, doi: 10.1126/science.1208336
  * Scalo et al. (2007) Scalo, J., Kaltenegger, L., Segura, A., et al. 2007, Astrobiology, 7, 85, doi: 10.1089/ast.2006.0125
  * Segura et al. (2005) Segura, A., Kasting, J. F., Meadows, V., et al. 2005, Astrobiology, 5, 706, doi: 10.1089/ast.2005.5.706
  * Selsis et al. (2007) Selsis, F., Kasting, J. F., Levrard, B., et al. 2007, A&A, 476, 1373, doi: 10.1051/0004-6361:20078091
  * Shields et al. (2016) Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Physics Reports, 663, 1, doi: 10.1016/j.physrep.2016.10.003
  * Shields et al. (2014) Shields, A. L., Bitz, C. M., Meadows, V. S., Joshi, M. M., & Robinson, T. D. 2014, The Astrophysical Journal, 785, L9, doi: 10.1088/2041-8205/785/1/L9
  * Shields et al. (2019) Shields, A. L., Bitz, C. M., & Palubski, I. 2019, The Astrophysical Journal Letters, 884, L2, doi: 10.3847/2041-8213/ab44ce
  * Shields et al. (2013) Shields, A. L., Meadows, V. S., Bitz, C. M., et al. 2013, Astrobiology, 13, 715, doi: 10.1089/ast.2012.0961
  * Showman et al. (2013) Showman, A. P., Wordsworth, R. D., Merlis, T. M., & Kaspi, Y. 2013, Comparative Climatology of Terrestrial Planets, doi: 10.2458/azu_uapress_9780816530595-ch12
  * Stevens (2005) Stevens, B. 2005, Annual Review of Earth and Planetary Sciences, 33, 605, doi: 10.1146/annurev.earth.33.092203.122658
  * Tarter et al. (2007) Tarter, J. C., Backus, P. R., Mancinelli, R. L., et al. 2007, Astrobiology, 7, 30, doi: 10.1089/ast.2006.0124
  * Thomson & Vallis (2019) Thomson, S. I., & Vallis, G. K. 2019, Quarterly Journal of the Royal Meteorological Society, 145, 2627, doi: 10.1002/qj.3582
  * Tian & Ida (2015) Tian, F., & Ida, S. 2015, Nature Geoscience, 8, 177, doi: 10.1038/ngeo2372
  * Turbet et al. (2016) Turbet, M., Leconte, J., Selsis, F., et al. 2016, Astronomy & Astrophysics, 596, A112, doi: 10.1051/0004-6361/201629577
  * van de Griend & Owe (1994) van de Griend, A. A., & Owe, M. 1994, Water Resources Research, 30, 181, doi: 10.1029/93WR02747
  * Walker et al. (1981) Walker, J. C. G., Hays, P. B., & Kasting, J. F. 1981, Journal of Geophysical Research, 86, 9776, doi: 10.1029/JC086iC10p09776
  * Wieczorek (2007) Wieczorek, M. 2007, Treatise on Geophysics, 165, doi: 10.1016/B978-044452748-6.00156-5
  * Wolf (2021) Wolf, E. 2021, GitHub repository.  https://github.com/storyofthewolf/ExoCAM
  * Wolf et al. (2022) Wolf, E. T., Kopparapu, R., Haqq-Misra, J., & Fauchez, T. J. 2022, The Planetary Science Journal, 3, 7, doi: 10.3847/PSJ/ac3f3d
  * Wolf & Toon (2015) Wolf, E. T., & Toon, O. B. 2015, Journal of Geophysical Research: Atmospheres, 120, 5775, doi: 10.1002/2015JD023302
  * Yang et al. (2019a) Yang, H., Komacek, T. D., & Abbot, D. S. 2019a, The Astrophysical Journal, 876, L27, doi: 10.3847/2041-8213/ab1d60
  * Yang & Abbot (2014) Yang, J., & Abbot, D. S. 2014, The Astrophysical Journal, 784, 155, doi: 10.1088/0004-637X/784/2/155
  * Yang et al. (2019b) Yang, J., Abbot, D. S., Koll, D. D. B., Hu, Y., & Showman, A. P. 2019b, The Astrophysical Journal, 871, 29, doi: 10.3847/1538-4357/aaf1a8
  * Yang et al. (2014a) Yang, J., Boué, G., Fabrycky, D. C., & Abbot, D. S. 2014a, The Astrophysical Journal, 787, L2, doi: 10.1088/2041-8205/787/1/L2
  * Yang et al. (2013) Yang, J., Cowan, N. B., & Abbot, D. S. 2013, The Astrophysical Journal, 771, L45, doi: 10.1088/2041-8205/771/2/L45
  * Yang et al. (2014b) Yang, J., Liu, Y., Hu, Y., & Abbot, D. S. 2014b, The Astrophysical Journal, 796, L22, doi: 10.1088/2041-8205/796/2/L22
  * Zhang & Yang (2020) Zhang, Y., & Yang, J. 2020, The Astrophysical Journal, 901, L36, doi: 10.3847/2041-8213/abb87f


# On the Relationship Between Iterated Statistical Linearization and
Quasi–Newton Methods

Anton Kullberg, Martin A. Skoglund, Isaac Skog, ,
and Gustaf Hendeby, .
This work was partially supported by the Wallenberg AI, Autonomous Systems and
Software Program (WASP) funded by the Knut and Alice Wallenberg
Foundation.Anton Kullberg, Martin Skoglund, and Gustaf Hendeby are all with
Linköping University, Sweden. Martin A. Skoglund is also affiliated with
Eriksholm Research Centre, Snekkersten, Denmark. (email: {anton.kullberg,
martin.skoglund, gustaf.hendeby }@liu.se). Isaac Skog is with Uppsala
University, Sweden. (email: isaac.skog@angstrom.uu.se).

###### Abstract

This letter investigates relationships between iterated filtering algorithms
based on statistical linearization, such as the iterated unscented Kalman
filter (iukf), and filtering algorithms based on quasi–Newton (qn) methods,
such as the qn iterated extended Kalman filter (qn–iekf). Firstly, it is shown
that the iukf and the iterated posterior linearization filter (iplf) can be
viewed as qn algorithms, by finding a Hessian correction in the qn–iekf such
that the iplf iterate updates are identical to that of the qn–iekf. Secondly,
it is shown that the iplf/iukf update can be rewritten such that it is
approximately identical to the qn–iekf, albeit for an additional correction
term. This enables a richer understanding of the properties of iterated
filtering algorithms based on statistical linearization.

###### Index Terms:

Nonlinear filtering, Statistical linearization, Quasi–Newton

††publicationid: pubid:

##  I Introduction

State estimation in discrete–time state–space models with additive Gaussian
noise, i.e., models such as

$\displaystyle\mathbf{x}_{k+1}$
$\displaystyle=\mathbf{f}(\mathbf{x}_{k})+\mathbf{w}_{k}$
$\displaystyle\mathbf{w}_{k}$ $\displaystyle\stackrel{{\scriptstyle\text{\tiny
i.i.d.}}}{{\sim}}\mathcal{N}(\mathbf{0},\mathbf{Q})$ (1a)
$\displaystyle\mathbf{y}_{k}$
$\displaystyle=\mathbf{h}(\mathbf{x}_{k})+\mathbf{e}_{k}$
$\displaystyle\mathbf{e}_{k}$ $\displaystyle\stackrel{{\scriptstyle\text{\tiny
i.i.d.}}}{{\sim}}\mathcal{N}(\mathbf{0},\mathbf{R})$ (1b)

has been studied extensively for decades. Here,
$\mathbf{x}_{k},~{}\mathbf{y}_{k},~{}\mathbf{w}_{k},~{}\mathbf{e}_{k}$ are the
state, measurement, process noise, and measurement noise, respectively. The
filtering problem is then to compute the marginal distributions
$p(\mathbf{x}_{k}|\mathbf{y}_{1:k})$, given a sequence of measurements
$\mathbf{y}_{1:k}$. In the case $\mathbf{f}$ and $\mathbf{h}$ are linear
functions, this is analytically tractable and the solution is given by the
Kalman filter, which is the optimal estimator in the mean–squared error sense
[1]. For nonlinear state–space models, analytical solutions generally do not
exist and approximate inference techniques have been developed for these
cases. The extended Kalman filter (ekf), introduced alongside the original
Kalman filter, was a simple way of extending the Kalman filter to the
nonlinear case [1]. The original ekf linearizes the dynamical model
$\mathbf{f}$ and the measurement model $\mathbf{h}$ at the current state
estimate using a first–order Taylor expansion and then applies the standard
Kalman filter updates. Since the ekf was developed, it has been identified
that this choice of linearization point may be suboptimal. This has lead to
the development of the iterated extended Kalman filter (iekf) [2, 3, 4]. This
is a family of approximate inference techniques that attempt to find a better
linearization point for the measurement model, which boils down to iterating
the measurement update a number of times. This family of inference techniques
include the line–search iekf, quasi–Newton (qn) iekf (qn–iekf) etc., which are
commonly referred to as damped iekfs [5].

Alongside the development of the iekf, other strategies for nonlinear
filtering were developed to circumvent the need for analytical linearization.
Particularly, the unscented Kalman filter (ukf), based on the unscented
transform was developed as a competitive alternative to the ekf [6, 7]. The
ukf is essentially based on deterministically sampling the prior and
propagating a set of “sigma points” through the nonlinear function, whereafter
an approximate Gaussian distribution can be constructed based on the
transformed points. This is done in both the time update, as well as the
measurement update. The ukf has since been shown to be equivalent to
statistically linearizing the nonlinear models and then applying the standard
Kalman recursions [8]. Therefore, both the ukf and other deterministically
sampled sigma–point methods can be interpreted as linearization–based
nonlinear filtering algorithms. Similarly to the ekf, the ukf has also been
extended to the iterated ukf (iukf) [9], which mimics the development of the
iekf closely.

Recently, another strategy for linearization–based filtering was introduced as
yet another alternative to the ekf and the sigma–point based filters, namely
the iterated posterior linearization filter (iplf) [10]. The iplf is based on
the idea of linearizing around the _posterior_
$p(\mathbf{x}_{k}|\mathbf{y}_{1:k})$. As the posterior is not available, the
iplf constructs an approximate posterior
$q_{i}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$ and then iterates the measurement
update, each time performing statistical linearization around the current
approximate posterior $q_{i}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$. Similarly to
the iekf, the iplf has been shown to diverge for some particular problems [9],
which has led to the development of damped versions of the iplf [11].

To gain a deeper understanding of the properties of the iterated statistical
linearization filters, such as the iukf and iplf, we seek to connect this
family of methods to classical qn methods, such as the qn–iekf. To that end,
we show that the iplf and the iukf (and other iterated statistical
linearization based filters) can be interpreted both as: (i) an exact qn
method, and (ii) an approximate qn method. We begin with mathematical
preliminaries and then detail the connections between qn and the iplf and
iukf.

##  II Preliminaries

Mathematical preliminaries are restated here for completeness.

###  II-A Statistical Linearization

Given a nonlinear model

$\mathbf{z}=\mathbf{g}(\mathbf{x}),$

we wish to find an affine representation

$\mathbf{g}(\mathbf{x})\approx\mathbf{A}\mathbf{x}+\mathbf{b}+\eta,$ (2)

with $\eta\sim\mathcal{N}(\eta;\mathbf{0},\bm{\Omega})$. In this affine
representation, there are three free parameters, $\mathbf{A},\mathbf{b}$, and
$\bm{\Omega}$. Statistical linearization finds these parameters by linearizing
w.r.t. a distribution $p(\mathbf{x})$. Practically, one may think of this as
constructing an affine function that best fits a number of samples of
$p(\mathbf{x})$ transformed through $\mathbf{g}(\mathbf{x})$. Assuming that
$p(\mathbf{x})=\mathcal{N}(\mathbf{x};\hat{\mathbf{x}},\mathbf{P})$,
statistical linearization selects the affine parameters as

$\displaystyle\mathbf{A}$ $\displaystyle=\Psi^{\top}\mathbf{P}^{-1}$ (3a)
$\displaystyle\mathbf{b}$
$\displaystyle=\bar{\mathbf{z}}-\mathbf{A}\hat{\mathbf{x}}$ (3b)
$\displaystyle\bm{\Omega}$
$\displaystyle=\Phi-\mathbf{A}\mathbf{P}\mathbf{A}^{\top}$ (3c)
$\displaystyle\bar{\mathbf{z}}$
$\displaystyle=\mathbb{E}[\mathbf{g}(\mathbf{x})]$ (3d) $\displaystyle\Psi$
$\displaystyle=\mathbb{E}[(\mathbf{x}-\hat{\mathbf{x}})(\mathbf{g}(\mathbf{x})-\bar{\mathbf{z}})^{\top}]$
(3e) $\displaystyle\Phi$
$\displaystyle=\mathbb{E}[(\mathbf{g}(\mathbf{x})-\bar{\mathbf{z}})(\mathbf{g}(\mathbf{x})-\bar{\mathbf{z}})^{\top}],$
(3f)

where the expectations are taken w.r.t. $p(\mathbf{x})$. The major difference
from analytical linearization is that ${\bm{\Omega}\neq 0}$, which implies
that the error in the linearization is captured. Typically, the expectations
in 3 are not analytically tractable and thus, practically, one often resorts
to some numerical integration technique.

###  II-B Quasi–Newton Optimization

A general nonlinear least–squares minimization problem is given by

$\hat{\mathbf{x}}=\operatorname*{\arg\min}_{\mathbf{x}}V(\mathbf{x}),\quad
V(\mathbf{x})=\frac{1}{2}r(\mathbf{x})^{\top}r(\mathbf{x}).$ (4)

One particular family of methods for solving these problems, is the Newton
family. This family of methods essentially finds the minimizing argument of 4
by starting at an initial guess $\mathbf{x}_{0}$ and iterating

$\mathbf{x}_{i+1}=\mathbf{x}_{i}-(\nabla^{2}V(\mathbf{x}_{i}))^{-1}\nabla
V(\mathbf{x}_{i}).$ (5)

Here, $\nabla^{2}V(\mathbf{x}_{i})$ and $\nabla V(\mathbf{x}_{i})$ are the
Hessian and the gradient of $V$ evaluated at $(\mathbf{x}_{i})$, respectively.
Note that convergence of the iterates typically benefit from step–size
correction, see _e.g._ , [12]. For nonlinear least–squares problems, the
gradient and Hessian are given by

$\displaystyle\nabla V(\mathbf{x})$
$\displaystyle=\mathbf{J}^{\top}(\mathbf{x})r(\mathbf{x})$ (6a)
$\displaystyle\nabla^{2}V(\mathbf{x})$
$\displaystyle=\mathbf{J}^{\top}(\mathbf{x})\mathbf{J}(\mathbf{x})+\sum_{i=1}^{n_{r}}[r(\mathbf{x})]_{i}\nabla^{2}[r(\mathbf{x})]_{i}$
(6b) $\displaystyle\mathbf{J}(\mathbf{x})$
$\displaystyle=\frac{dr(\mathbf{s})}{d\mathbf{s}}\bigg{\lvert}_{\mathbf{s}=\mathbf{x}}$
(6c)

where $[r(\mathbf{x})]_{i}$ is the $i$th component of $r(\mathbf{x})$ and
$n_{r}$ is the dimension of $r(\mathbf{x})$. As the Hessian of the cost
function can be computationally expensive to evaluate, approximate versions of
Newton’s method have been developed. In particular, the Gauss–Newton method
approximates the Hessian as

$\nabla^{2}V(\mathbf{x})\approx\mathbf{J}(\mathbf{x})^{\top}\mathbf{J}(\mathbf{x}),$
(7)

thus only requiring first–order information. As such, it is a Quasi–Newton
(qn) method since it operates as a Newton method with an approximate Hessian.
This approximation may be bad far from the optimum, which may affect
convergence. A remedy is to either approximate the Hessian directly in some
other way, or by introducing a correction term to the Gauss–Newton approximate
Hessian as

$\nabla^{2}V(\mathbf{x})\approx\mathbf{J}(\mathbf{x})^{\top}\mathbf{J}(\mathbf{x})+\mathbf{T},$
(8)

where $\mathbf{T}$ is supposed to capture second–order information and can be
chosen in a variety of ways, see [12] for an overview.

###  II-C Quasi–Newton iekf

The cost function for the general iekf is of the form 4 with [5]

$r(\mathbf{x})=\begin{bmatrix}\mathbf{R}^{-1/2}(\mathbf{y}_{k}-\mathbf{h}(\mathbf{x}))\\\
\mathbf{P}_{k|k-1}^{-1/2}(\hat{\mathbf{x}}_{k|k-1}-\mathbf{x})\end{bmatrix}.$
(9)

Hence, the Jacobian $\mathbf{J}(\mathbf{x})$ is given by

$\mathbf{J}(\mathbf{x})=-\begin{bmatrix}\mathbf{R}^{-1/2}\mathbf{H}(\mathbf{x})\\\
\mathbf{P}_{k|k-1}^{-1/2}\end{bmatrix},~{}\text{where}~{}\mathbf{H}(\mathbf{x})=\frac{d\mathbf{h}(\mathbf{s})}{d\mathbf{s}}\bigg{\lvert}_{\mathbf{s}=\mathbf{x}}.$
(10)

Now, using 5 and 8, the qn–iekf iterate update is given by [5]

$\displaystyle\mathbf{x}_{i+1}$
$\displaystyle=\hat{\mathbf{x}}+\mathbf{K}^{q}_{i}(\mathbf{y}_{k}-\mathbf{h}_{i}-\mathbf{H}_{i}\tilde{\mathbf{x}}_{i})-\mathbf{S}^{q}_{i}\mathbf{T}_{i}\tilde{\mathbf{x}}_{i}$
(11a) $\displaystyle\mathbf{P}_{i+1}$
$\displaystyle=\mathbf{P}-\mathbf{P}\mathbf{H}_{i}^{\top}(\mathbf{H}_{i}\mathbf{P}\mathbf{H}_{i}^{\top}+\mathbf{R})^{-1}\mathbf{H}_{i}\mathbf{P}$
(11b) $\displaystyle\mathbf{S}^{q}_{i}$
$\displaystyle\triangleq\left(\mathbf{H}_{i}^{\top}\mathbf{R}_{i}^{-1}\mathbf{H}_{i}+\mathbf{P}^{-1}+\mathbf{T}_{i}\right)^{-1}$
(11c) $\displaystyle\mathbf{K}^{q}_{i}$
$\displaystyle\triangleq\mathbf{S}^{q}_{i}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1},$
(11d)

with simplified notation
$\tilde{\mathbf{x}}_{i}=\hat{\mathbf{x}}-\mathbf{x}_{i}$,
$\mathbf{h}_{i}=\mathbf{h}(\mathbf{x}_{i})$,
$\mathbf{H}_{i}=\mathbf{H}(\mathbf{x}_{i})$, and
$\mathbf{P}=\mathbf{P}_{k|k-1}$. Further, $\mathbf{T}_{i}$ is the current
Hessian correction which can be chosen freely, see [5].

###  II-D Iterated Posterior Linearization Filter

The iplf was initially developed through minimizing the Kullback–Leibler (kl)
divergence between the true marginal posterior
$p(\mathbf{x}_{k}|\mathbf{y}_{1:k})$ and an approximation
$q(\mathbf{x}_{k}|\mathbf{y}_{1:k})$, _i.e._ ,

$q(\mathbf{x}_{k}|\mathbf{y}_{1:k})=\min_{q}D_{\textsc{kl}}\left(p\|q\right),$
(12)

where

$D_{\textsc{kl}}\left(p\|q\right)=\int
p(\mathbf{x}_{k}|\mathbf{y}_{1:k})\log\frac{p(\mathbf{x}_{k}|\mathbf{y}_{1:k})}{q(\mathbf{x}_{k}|\mathbf{y}_{1:k})}d\mathbf{x}_{k}.$

In particular, the divergence is used to find an “optimal” affine
approximation of the observation model, subsequently used as the “surrogate”
measurement model. However, as the objective 12 is not analytically tractable,
and also requires access to the true posterior
$p(\mathbf{x}_{k}|\mathbf{y}_{1:k})$, the iplf approximately minimizes this by
starting at some initial approximate posterior
$q_{0}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$. The approximate posterior is then
iteratively refined until a stopping condition is met. Essentially, the
recursions are similar to an iekf using a statistically linearized model. By
adapting the notation of [13], the iterate update of the iplf can be written
as

$\displaystyle\mathbf{x}_{i+1}^{\text{\tiny{IPLF}}}$
$\displaystyle=\hat{\mathbf{x}}+\mathbf{K}_{i}(\mathbf{y}_{k}-\mathbf{h}_{i}-\mathbf{H}_{i}\tilde{\mathbf{x}}_{i})$
(13a) $\displaystyle\mathbf{P}_{i+1}$
$\displaystyle=\mathbf{P}-\mathbf{K}_{i}\mathbf{S}_{i}\mathbf{K}_{i}^{\top}$
(13b) $\displaystyle\mathbf{K}_{i}$
$\displaystyle\triangleq\mathbf{P}_{i}\mathbf{H}_{i}^{\top}\mathbf{S}_{i}$
(13c) $\displaystyle\mathbf{S}_{i}$
$\displaystyle\triangleq\left(\mathbf{H}_{i}\mathbf{P}_{i}\mathbf{H}_{i}^{\top}+\mathbf{R}+\bm{\Omega}_{i}\right)^{-1},$
(13d)

where $\bm{\Omega}_{i},~{}\mathbf{h}_{i}$ and $\mathbf{H}_{i}$ are found
through statistical linearization of $\mathbf{h}$. Note that
$\mathbf{h}_{i}=\mathbf{b}+\mathbf{A}\hat{\mathbf{x}}$ and
$\mathbf{H}_{i}=\mathbf{A}$, see Section II-A. The iterates are initialized as
$\mathbf{x}_{i}=\hat{\mathbf{x}},~{}\mathbf{P}_{i}=\mathbf{P}$ and the updates
are then iterated until the approximate posterior
$q_{i}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$ does not significantly change, as
measured by

$D_{\textsc{kl}}(q_{i+1}\|q_{i}).$

Note that by choosing the covariance update 13b as
$\mathbf{P}_{i+1}=\mathbf{P}$ until the last iteration, the iukf presented in
[9] is obtained. We will now consider two relationships between the iplf
update 13a and the qn–iekf update 11a.

##  III Exact Quasi–Newton

Next, we show that the iplf can be viewed as an _exact_ qn method, _i.e._ ,
that it corresponds to qn with a particular choice of Hessian correction
$\mathbf{T}$. More precisely, we find a Hessian correction $\mathbf{T}_{i}$
such that the qn–iekf update 11a is equal to the iplf update 13a.

First, let
$\epsilon_{i}\triangleq\mathbf{y}_{k}-\mathbf{h}_{i}-\mathbf{H}_{i}\tilde{\mathbf{x}}_{i}$.
Now, setting 11a equal to 13a yields

$\displaystyle\mathbf{K}^{q}_{i}\epsilon_{i}-\mathbf{S}^{q}_{i}\mathbf{T}_{i}\tilde{\mathbf{x}}_{i}=\mathbf{K}_{i}\epsilon_{i}\iff$
$\displaystyle\mathbf{S}^{q}_{i}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\epsilon_{i}-\mathbf{S}^{q}_{i}\mathbf{T}_{i}\tilde{\mathbf{x}}_{i}=\mathbf{K}_{i}\epsilon_{i}\iff$
$\displaystyle\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\epsilon_{i}-\mathbf{T}_{i}\tilde{\mathbf{x}}_{i}=(\mathbf{S}^{q}_{i})^{-1}\mathbf{K}_{i}\epsilon_{i}\iff$
$\displaystyle\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\epsilon_{i}-\mathbf{T}_{i}\tilde{\mathbf{x}}_{i}=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\mathbf{H}_{i}+\mathbf{P}^{-1}+\mathbf{T}_{i})\mathbf{K}_{i}\epsilon_{i}\iff$
$\displaystyle\underbrace{(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\\!-\\!\big{(}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\mathbf{H}_{i}\\!+\\!\mathbf{P}^{-1}\big{)}\mathbf{K}_{i})}_{(*)}\epsilon_{i}=\mathbf{T}_{i}(\tilde{\mathbf{x}}_{i}\\!+\\!\mathbf{K}_{i}\epsilon_{i}).$

Now, note that $(*)$ can be written as

$\displaystyle(*)$
$\displaystyle=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}-\big{(}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\mathbf{H}_{i}+\mathbf{P}^{-1}\big{)}\mathbf{K}_{i})$
$\displaystyle=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}-\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\mathbf{H}_{i}\mathbf{P}_{i}\mathbf{H}_{i}^{\top}\mathbf{S}_{i}-\mathbf{P}^{-1}\mathbf{P}_{i}\mathbf{H}_{i}^{\top}\mathbf{S}_{i})$
$\displaystyle=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\left(\mathbf{S}_{i}^{-1}-\mathbf{H}_{i}\mathbf{P}_{i}\mathbf{H}_{i}^{\top}\right)\mathbf{S}_{i}-\mathbf{P}^{-1}\mathbf{P}_{i}\mathbf{H}_{i}^{\top}\mathbf{S}_{i})$
$\displaystyle=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\bm{\Omega}_{i}+\big{(}\mathbf{I}-\mathbf{P}^{-1}\mathbf{P}_{i}\big{)}\mathbf{H}_{i}^{\top})\mathbf{S}_{i}.$

Thus, we have

$\mathbf{T}_{i}(\tilde{\mathbf{x}}_{i}\\!+\\!\mathbf{K}_{i}\epsilon_{i})=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\bm{\Omega}_{i}+\big{(}\mathbf{I}-\mathbf{P}^{-1}\mathbf{P}_{i}\big{)}\mathbf{H}_{i}^{\top})\mathbf{S}_{i}\epsilon_{i}.$

Letting

$\displaystyle\mathbf{s}_{i}$
$\displaystyle\triangleq\tilde{\mathbf{x}}_{i}+\mathbf{K}_{i}\epsilon_{i}=\mathbf{x}_{i+1}^{\text{\tiny{IPLF}}}-\mathbf{x}_{i}$
(14a) $\displaystyle\mathbf{p}_{i}$
$\displaystyle\triangleq(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\bm{\Omega}_{i}+\big{(}\mathbf{I}-\mathbf{P}^{-1}\mathbf{P}_{i}\big{)}\mathbf{H}_{i}^{\top})\mathbf{S}_{i}\epsilon_{i},$
(14b)

we have

$\mathbf{T}_{i}\mathbf{s}_{i}=\mathbf{p}_{i}.$ (15)

which is similar to the _secant equation_ , see _e.g._[12, p. 24]. Hence, we
can follow a similar reasoning and procedure to find a solution. That is, we
impose that $\mathbf{T}_{i}$ be symmetric and that it is “close” to
$\mathbf{T}_{i-1}$, in some sense. Thus, we find $\mathbf{T}_{i}$ by
minimizing

$\displaystyle\min_{\mathbf{T}}$
$\displaystyle\lVert\mathbf{T}-\mathbf{T}_{i-1}\rVert_{\mathbf{W}}$ (16)
subject to
$\displaystyle\mathbf{T}=\mathbf{T}^{\top},\quad\mathbf{T}\mathbf{s}_{i}=\mathbf{p}_{i}.$

Here,
$\lVert\mathbf{A}\rVert_{\mathbf{W}}=\lVert\mathbf{W}^{1/2}\mathbf{A}\mathbf{W}^{1/2}\rVert$
and $\mathbf{W}$ is a nonsingular symmetric matrix. Now, let
$\mathbf{s}_{i},\mathbf{p}_{i}$ be in $\mathbb{R}^{n}$. Then, [14, Theorem
7.3] states that for any $\mathbf{c}\in\mathbb{R}^{n}$ such that
$\mathbf{c}^{\top}\mathbf{s}_{i}>0$ and
$\mathbf{W}\mathbf{c}=\mathbf{W}^{-1}\mathbf{s}_{i}$, the solution to 16 is
given by

$\mathbf{T}_{i}=\mathbf{T}_{i-1}+\frac{(\mathbf{p}_{i}-\mathbf{T}_{i-1}\mathbf{s}_{i})\mathbf{c}^{\top}+\mathbf{c}(\mathbf{p}_{i}-\mathbf{T}_{i-1}\mathbf{s}_{i})^{\top}}{\mathbf{c}^{\top}\mathbf{s}_{i}}\\\
-\frac{(\mathbf{p}_{i}-\mathbf{T}_{i-1}\mathbf{s}_{i})^{\top}\mathbf{s}_{i}}{(\mathbf{c}^{\top}\mathbf{s}_{i})^{2}}\mathbf{c}\mathbf{c}^{\top}.$
(17)

In particular, we choose $\mathbf{c}=\mathbf{s}_{i}$ which guarantees
$\mathbf{c}^{\top}\mathbf{s}_{i}=\mathbf{s}_{i}^{\top}\mathbf{s}_{i}>0$ as
long as $\mathbf{s}_{i}\neq 0$. Note that the case $\mathbf{s}_{i}=0$ is not
of interest, as this means that $\mathbf{x}_{i+1}=\mathbf{x}_{i}$ and hence,
fixed point convergence has been reached. The resulting recursion 17 is
generally known as the Powell symmetric Broyden (psb) update [14, p. 70].
However, even though the structure is the same as psb, the definition of
$\mathbf{p}_{i}$ is not. Thus, with these particular definitions, 17 does not
correspond exactly to psb.

Note that we may also follow the typical derivation of the
Davidon–Fletcher–Powell (dfp) update [12] by imposing that $\mathbf{T}_{i}$ be
positive definite. However, this results in the condition
$\mathbf{s}_{i}^{\top}\mathbf{p}_{i}>0$, which generally is not the case nor
is it trivial to guarantee.

Note that even though the derivation here focuses on the iplf, a similar
relationship holds for the ickf, iukf, and other iterated sigma–point filters
as well. Regardless, this shows that iterated filtering algorithms based on
statistical linearization can be viewed as qn methods with a particular choice
of Hessian correction $\mathbf{T}_{i}$. This can, for instance, lead to the
understanding of the convergence properties of the ickf, iukf, and iplf from a
qn perspective. Further, we can view the captured linearization error
$\bm{\Omega}$ as an adaptation of the Hessian approximation in order to
account for lost second order information. Next, we interpret the Hessian
correction 17 by considering the special case of the iukf and ickf, which
simplifies 14b.

###  III-A iukf/ickf interpretation

The iukf [9] essentially performs the updates 13a, but keeps
$\mathbf{P}_{i+1}=\mathbf{P}$ fixed until the last iteration. In this case,
14b simplifies to

$\mathbf{p}_{i}=\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\bm{\Omega}_{i}\mathbf{S}_{i}\epsilon_{i}.$
(18)

Henceforth, we assume that $\mathbf{s}_{i}$ is constant and given, such that
we can interpret the correction 17 only in terms of changes in 18. Inspecting
the components of 18, it is “weighted” with
$\mathbf{R}^{-1}\bm{\Omega}_{i}\mathbf{S}_{i}$. Hence, with decreasing
measurement uncertainty $\mathbf{R}$, the Hessian correction grows “larger”,
as the measurement model is more precise. Similarly, as the linearization
error $\bm{\Omega}_{i}$ grows, the Hessian correction also increases, which
makes sense as it indicates that the measurement function $\mathbf{h}$ is
highly nonlinear and the Hessian approximation 7 is most likely poor and needs
more correction. Also, note that if $\bm{\Omega}_{i}=\mathbf{0}$, _i.e._ , if
the model is completely linear at the current iterate, the Hessian is only
corrected according to the iterate difference $\mathbf{s}_{i}$. Lastly, as the
innovation covariance $\mathbf{S}_{i}^{-1}$ decreases, the correction grows,
essentially also indicating that the measurement carries a lot of information
that can be exploited. This interpretation should approximately hold for the
iplf as well, as usually $\mathbf{P}^{-1}\mathbf{P}_{i}\approx\mathbf{I}$.
However, a detailed analysis of the exact behavior of the iplf is non–trivial,
as $\mathbf{H}_{i},\bm{\Omega}_{i}$ and $\mathbf{S}_{i}$ all depend on the
previous iterate $\mathbf{P}_{i}$.

##  IV Approximate Quasi–Newton

As an alternative to the exact view in Section III, we can also view the iplf
as an _approximate_ Quasi–Newton method. Essentially, it boils down to
modifying 13a such that it approximately takes the form 11a.

Start with 13d and write

$\mathbf{S}_{i}=\left(\mathbf{H}_{i}\mathbf{P}_{i}\mathbf{H}_{i}^{\top}+\mathbf{R}+\bm{\Omega}_{i}\right)^{-1}=(\mathbf{R}+\bm{\Omega}_{i})^{-1}\big{(}\mathbf{I}-\\\
\mathbf{H}_{i}(\mathbf{H}_{i}^{\top}(\mathbf{R}+\bm{\Omega}_{i})^{-1}\mathbf{H}_{i}+\mathbf{P}_{i}^{-1})^{-1}\mathbf{H}_{i}^{\top}(\mathbf{R}+\bm{\Omega}_{i})^{-1}\big{)}.$

Now, note that

$(\mathbf{H}_{i}^{\top}(\mathbf{R}+\bm{\Omega}_{i})^{-1}\mathbf{H}_{i}+\mathbf{P}_{i}^{-1})^{-1}\\\
=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\mathbf{H}_{i}+\mathbf{P}_{i}^{-1}-\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}(\mathbf{R}^{-1}+\bm{\Omega}_{i}^{-1})^{-1}\mathbf{R}^{-1}\mathbf{H}_{i})^{-1}\\\
=(\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}\mathbf{H}_{i}+\mathbf{P}_{i}^{-1}+\tilde{\mathbf{T}}_{i})^{-1}\triangleq\tilde{\mathbf{S}}^{q}_{i},$

with

$\displaystyle\tilde{\mathbf{T}}_{i}$
$\displaystyle\triangleq-\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}(\mathbf{R}^{-1}+\bm{\Omega}_{i}^{-1})^{-1}\mathbf{R}^{-1}\mathbf{H}_{i}$
$\displaystyle=-\mathbf{H}_{i}^{\top}\left(\mathbf{R}^{-1}-\left(\mathbf{R}+\bm{\Omega}_{i}\right)^{-1}\right)\mathbf{H}_{i}.$

Plugging this into 13c yields

$\mathbf{P}_{i}\mathbf{H}_{i}^{\top}\left(\mathbf{I}-(\mathbf{R}+\bm{\Omega}_{i})^{-1}\mathbf{H}_{i}\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}\right)(\mathbf{R}+\bm{\Omega}_{i})^{-1}\\\
=\mathbf{P}_{i}\left((\tilde{\mathbf{S}}^{q}_{i})^{-1}-\mathbf{H}_{i}^{\top}(\mathbf{R}+\bm{\Omega}_{i})^{-1}\mathbf{H}_{i}\right)\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}(\mathbf{R}+\bm{\Omega}_{i})^{-1}\\\
=\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}(\mathbf{R}+\bm{\Omega}_{i})^{-1}\\\
=\underbrace{\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}}_{\triangleq\tilde{\mathbf{K}}^{q}_{i}}-\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}(\mathbf{R}^{-1}+\bm{\Omega}_{i}^{-1})^{-1}\mathbf{R}^{-1}.$

Hence, with
$\epsilon_{i}=\mathbf{y}_{k}-\mathbf{h}_{i}-\mathbf{H}_{i}\tilde{\mathbf{x}}_{i}$,
the iterate update becomes

$\mathbf{x}_{i+1}=\hat{\mathbf{x}}+\tilde{\mathbf{K}}^{q}_{i}\epsilon_{i}-\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}(\mathbf{R}^{-1}+\bm{\Omega}_{i}^{-1})^{-1}\mathbf{R}^{-1})\epsilon_{i}\\\
=\hat{\mathbf{x}}+\tilde{\mathbf{K}}^{q}_{i}\epsilon_{i}-\tilde{\mathbf{S}}^{q}_{i}\overbrace{(-\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}(\mathbf{R}^{-1}+\bm{\Omega}_{i}^{-1})^{-1}\mathbf{R}^{-1}\mathbf{H}_{i})}^{\tilde{\mathbf{T}}_{i}}\tilde{\mathbf{x}}_{i}\\\
-\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}\mathbf{R}^{-1}(\mathbf{R}^{-1}+\bm{\Omega}_{i}^{-1})^{-1}\mathbf{R}^{-1}(\mathbf{y}_{k}-\mathbf{h}_{i})\\\
=\hat{\mathbf{x}}+\tilde{\mathbf{K}}^{q}_{i}\epsilon_{i}-\tilde{\mathbf{S}}^{q}_{i}\tilde{\mathbf{T}}_{i}\tilde{\mathbf{x}}_{i}-\tilde{\mathbf{S}}^{q}_{i}\mathbf{H}_{i}^{\top}(\mathbf{R}^{-1}-(\mathbf{R}+\bm{\Omega}_{i})^{-1})(\mathbf{y}_{k}-\mathbf{h}_{i}).$

Hence, the iplf can be interpreted as performing modified Quasi-Newton with a
specific choice of Hessian correction $\mathbf{T}_{i}$ and an additional
correction term in the iterate update. The additional term essentially
corrects the iterate depending on how well the (nonlinear) observation model
captures the current observation at the iterate $\mathbf{x}_{i}$. Note that
both the Hessian correction and additional term are weighted by the “degree of
linearity” of the observation model at the specific iterate. Hence, if the
observation model is completely linear at the current iterate,
$\bm{\Omega}_{i}=\mathbf{0}$ which yields $\tilde{\mathbf{T}}_{i}=\mathbf{0}$.
Thus, there is no Hessian correction nor is there any additional correction
term. Hence, the iterate update reduces to the standard iekf update.

##  V Conclusion

In this letter, we have shown that iterated filtering algorithms based on
statistical linearization, such as the iukf and the iplf can be interpreted as
qn methods with a particular choice of Hessian correction. Through this
connection, we hope to enable a richer understanding of the properties of the
iukf, iplf, and other iterated sigma–point filters. In particular, it allows
us to view the use of statistical linearization as an adaptation of the
Hessian approximation in order to account for lost second order information.
It may also prove useful in determining the convergence differences between
the iukf and iplf as they differ only slightly from one another.

## References

  * [1] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” _J. of Basic Eng._ , vol. 82, no. 1, pp. 35–45, Mar. 1960.
  * [2] W. Denham and S. Pines, “Sequential estimation when measurement function nonlinearity is comparable to measurement error,” in _Guidance Control Conference_ , ser. Guidance, Navigation, and Control and Co-located Conferences, Aug. 1965.
  * [3] B. M. Bell and F. W. Cathey, “The iterated Kalman filter update as a Gauss-Newton method,” vol. 38, no. 2, pp. 294–297, Feb. 1993.
  * [4] P. J. G. Teunissen, “On the local convergence of the iterated extended Kalman filter,” in _Proceeding of the XX General Assembly of the IUGG, IAG Section IV_. Vienna: IAG, Aug. 1991, pp. 177–184.
  * [5] M. A. Skoglund, G. Hendeby, and D. Axehill, “Extended Kalman Filter Modifications Based on an Optimization view point,” in _2015 18th Int. Conf. on Inf. Fusion (FUSION)_ , Jul. 2015, pp. 1856–1861.
  * [6] S. Julier, J. Uhlmann, and H. Durrant-Whyte, “A New Approach for Filtering Nonlinear Systems,” in _Proc. of 1995 Am. Control Conf. - ACC’95_ , vol. 3, Seattle, WA, USA, 1995, pp. 1628–1632.
  * [7] S. Julier, “A General Method for Approximating Nonlinear Transformations of Probability Distributions,” Tech. Rep.
  * [8] T. Lefebvre, H. Bruyninckx, and J. De Schuller, “Comment on ”A new method for the nonlinear transformation of means and covariances in filters and estimators” [with authors’ reply],” _IEEE Trans. Automat. Contr._ , vol. 47, no. 8, pp. 1406–1409, Aug. 2002.
  * [9] M. A. Skoglund, F. Gustafsson, and G. Hendeby, “On Iterative Unscented Kalman Filter using Optimization,” in _2019 22th Int. Conf. on Inf. Fusion (FUSION)_ , Ottawa, ON, Canada, Jul. 2019, pp. 1–8.
  * [10] Á. F. García-Fernández, L. Svensson, and M. R. Morelande, “Iterated Statistical Linear Regression for Bayesian Updates,” in _17th Int. Conf. on Inf. Fusion (FUSION)_ , Salamanca, Spain, Jul. 2014, pp. 1–8.
  * [11] M. Raitoharju, L. Svensson, Á. F. García-Fernández, and R. Piché, “Damped Posterior Linearization Filter,” _IEEE Signal Processing Letters_ , vol. 25, no. 4, pp. 536–540, Apr. 2018.
  * [12] J. Nocedal and S. Wright, _Numerical Optimization (Springer Series in Operations Research and Financial Engineering)_ , 2nd ed. Springer, Jul. 2006.
  * [13] Á. F. García-Fernández, L. Svensson, M. R. Morelande, and S. Särkkä, “Posterior Linearization Filter: Principles and Implementation Using Sigma Points,” _IEEE Trans. Signal Process._ , vol. 63, no. 20, pp. 5561–5573, Oct. 2015.
  * [14] J. E. Dennis, Jr. and J. J. Moré, “Quasi-newton methods, motivation and theory,” _SIAM Review_ , vol. 19, no. 1, pp. 46–89, 1977.