v#1

# General relativistic precession and the long-term stability of the solar
system

Garett Brown${}^{1,2\;\star}$ and Hanno Rein1,2,3
1 Department of Physical and Environmental Sciences, University of Toronto at
Scarborough, Toronto, Ontario M1C 1A4, Canada,
2 Department of Physics, University of Toronto, Toronto, Ontario, M5S 3H4,
Canada,
3 Department of Astronomy and Astrophysics, University of Toronto, Toronto,
Ontario, M5S 3H4, Canada
$\star\;$<EMAIL_ADDRESS>

(Draft version: )

###### Abstract

The long-term evolution of the solar system is chaotic. In some cases, chaotic
diffusion caused by an overlap of secular resonances can increase the
eccentricity of planets when they enter into a linear secular resonance,
driving the system to instability. Previous work has shown that including
general relativistic contributions to the planets’ precession frequency is
crucial when modelling the solar system. It reduces the probability that the
solar system destabilizes within 5 Gyr by a factor of 60. We run 1280
additional $N$-body simulations of the solar system spanning 12.5 Gyr where we
allow the GR precession rate to vary with time. We develop a simple, unified,
Fokker-Planck advection-diffusion model that can reproduce the instability
time of Mercury with, without, and with time-varying GR precession. We show
that while ignoring GR precession does move Mercury’s precession frequency
closer to a resonance with Jupiter, this alone does not explain the increased
instability rate. It is necessary that there is also a significant increase in
the rate of diffusion. We find that the system responds smoothly to a change
in the precession frequency: There is no critical GR precession frequency
below which the solar system becomes significantly more unstable. Our results
show that the long-term evolution of the solar system is well described with
an advection-diffusion model.

###### keywords:

methods: numerical — gravitation — planets and satellites: dynamical evolution
and stability

##  1 Introduction

For many centuries, the motion of the solar system has been a subject of
interest to many cultures. Since Newton (1687) formulated his universal law of
gravitation, there has been a formal and precise way to describe this motion.
The paradigm shift of general relativity (GR) introduced by Einstein (1915)
further increased our understanding and ability to model the motion of the
solar system precisely. Shortly afterwards, post-Newtonian corrections derived
from the equations of general relativity provided a simpler approach to more
accurately model the discrepancies previously found through observation
(Eddington, 1923). Today, a combination of these precise numerical models
together with modern computational resources have made extensive studies on
the long-term behaviour of the solar system feasible. Modern numerical studies
of the solar system that focus primarily on its dynamical stability have
considered the classical solar system, which assumes no additional forces or
effects beyond Newtonian gravity, or they consider a constant strength of
general relativity (Batygin & Laughlin, 2008; Laskar & Gastineau, 2009; Zeebe,
2015; Abbot et al., 2021, 2022). We consider in more detail here how the
dynamical stability of the solar system depends on the general relativistic
perihelion precession rate of the planets.

We do this by artificially adjusting the strength of the first order post-
Newtonian correction over time, decreasing it from the values observed today
to completely absent. General modified versions of gravity have been developed
to account for dynamical features seen in galactic and galactic cluster scales
and are proposed to mimic the effect of dark matter (Milgrom, 1983; Moffat &
Toth, 2008; Verlinde, 2017). However, some of these theories result in changes
to the perihelion precession rate of Mercury which are not compatible with
current observations (Chan & Lee, 2023). We investigate a decrease in the GR
precession rate not because we think this mimics these modified theories of
gravity, but because it offers a clean way to experiment with the dynamics of
the solar system.

By comparing to previous work, we determine the effect the perihelion
precession rate of Mercury from general relativistic corrections has on the
probability that the system goes unstable. The results of this work reconfirm
the work of Laskar & Gastineau (2009), showing that the stability of the solar
system is highly dependent on the presence of general relativity. What is new
in this paper is that we also smoothly vary the strength of general
relativistic corrections to the perihelion precession over the lifetime of the
solar system and consider the implications for the very long-term evolution of
the solar system. We develop a physical model based on a diffusion process
(Mogavero & Laskar, 2021) that can explain the observed instability rate in
our numerical experiments. Our results show that a change in general
relativistic corrections does not only change the planets’ perihelion
precession rate, but also the rate at which the precession rate diffuses with
time. We find that the dependence of the diffusion coefficient and the general
relativistic corrections is smooth. We find no evidence of any critical
strength of GR that is required for stability.

To do this, we begin with a discussion of Mercury’s path to instability and
the role of secular resonances in Section 2. In Section 3 we describe the
numerical setup for our long-term ensemble of integrations of the solar
system, how we modify the first order post-Newtonian corrections that effect
the perihelion precession from general relativity, and the results of the
integrations. The diffusion model we use to interpret our results is developed
and compared to previous work in Section 4. In Section 5 we compare the model
to our new $N$-body data. Finally, we close in Section 6 with the conclusions
and implications.

##  2 The path to instability

The secular evolution of the solar system was first described by Laplace and
Lagrange (Laplace, 1775, 1776; Lagrange, 1776, 1778, 1781, 1782, 1783b, 1783a,
1784). See Laskar (2013) for a comprehensive historical review. This
description of the solar system averages over the mean motions (phases) of the
planets. The ground-breaking linear expansion is a perturbation theory first
order in the masses with a disturbing function second order in the
eccentricities and inclinations. It exhibits no changes to the semi-major axis
of the planets and over time shows only small changes to the eccentricities
and inclinations, but not by an amount to allow for orbit crossing scenarios.
The solutions contain the perihelion precession of planetary orbits as well as
the precession of the ascending nodes of the planets, but fundamentally they
are perfectly periodic. This initially seemed to _prove_ that the solar system
is dynamically stable for all time, but these successive approximations did
not provide rigorous bounds to actually prove it analytically (Poincaré,
1898).

From this Laplace-Lagrange linear expansion, the fundamental eigenmodes for
the variations in the eccentricities and inclinations of the solar system can
be calculated (Murray & Dermott, 1999). When Le Verrier (1840, 1841) continued
the perturbative expansion to higher order, significant terms emerged which
make crucial corrections to the linear equations and showed that the linear
theory could not be used for an indefinite period of time. In spite of
Poincaré (1899) proving the impossibility of an analytical solution to multi-
planetary dynamics over an infinite time interval, quasi-periodic
approximations were still the most accurate models at the time. Further
mathematical developments by Kolmogorov (1954), Arnol’d (1963), and Möser
(1962) (KAM) showed that for systems with a few degrees of freedom and in
small regions around the initial conditions, the trajectories remain
constrained to quasi-periodic solutions. These KAM toridal regions are only
isolating for systems with 2 degrees of freedom. As such, for systems with
more degrees of freedom, these constrained regions overlap and allow
trajectories to pass into chaotic regions and effectively diffuse through
phase space (Laskar, 2013). With the aid of computer algebra systems, the
integration of high-order secular expansions by Laskar (1985, 1986, 1990)
showed that the inner solar system is indeed chaotic. Recent developments in
high-order expansions continue to reveal insights into the dynamical richness
of the solar system unachievable by $N$-body integrations alone (Mogavero &
Laskar, 2021, 2022; Hoang et al., 2022).

Over short timescales of $\sim 1\,\rm{Myr}$ the evolution of the secular
frequencies is regular and can be calculated accurately either with computer
aided expansion to higher order, or with direct $N$-body simulations. However,
on timescales greater than $\sim 50\,\rm{Myr}$, the changes in secular
frequencies are chaotic (Laskar, 1990; Lithwick & Wu, 2011). This chaos in the
solar system is largely driven by overlapping secular resonances (Laskar,
1989). Various resonances driving chaos in the inner solar system are at play,
for example the Earth-Mars secular resonance
$2(\varpi_{4}-\varpi_{3})-(\Omega_{4}-\Omega_{3})$ and the Mercury-Venus-
Jupiter secular resonance $(\varpi_{1}-\varpi_{5})-(\Omega_{1}-\Omega_{2})$
(Laskar et al., 1992; Sussman & Wisdom, 1992), where
$\varpi_{1},...,\varpi_{8}$ are the longitude of the perihelia and
$\Omega_{1},...,\Omega_{8}$ are the longitude of the ascending nodes of the
planets. Sometimes these secular resonances are also expressed in terms of the
eigenfrequencies $g_{1},...,g_{8},s_{1},...,s_{8}$ of the solar system
associated with the precession of the perihelia and ascending nodes. We follow
the traditional convention to associate $g_{1}$ with Mercury, $g_{2}$ with
Venus, and so forth, with our discussion focusing mainly on Mercury and
Jupiter, and therefore $g_{1}$ and $g_{5}$. Thus, the Earth-Mars resonance
mentioned above would be $2(g_{4}-g_{3})-(s_{4}-s_{3})$.

Instabilities from the $g_{1}-g_{5}$ resonance can be seen in long-term solar
system integration of $N$-body models (Laskar & Gastineau, 2009; Zeebe, 2015;
Abbot et al., 2021, 2022; Brown & Rein, 2022). As the solar system evolves,
the planets push each other into or out of secular resonances through
exchanges of angular momentum (Laskar, 2000; Zakamska & Tremaine, 2004).
Reduced models help us to understand the nature of this process, specifically
how the $g_{1}-g_{5}$ resonance affects the eccentricity pumping of Mercury,
and how Mercury falls into the resonance in the first place after a diffusive
walk through phase space (Lithwick & Wu, 2011; Batygin et al., 2015; Mogavero
& Laskar, 2021). Although the frequencies of the secular eigenmodes of the
solar system are fixed in the simplest model, they do in fact show small
variations. In particular, the secular frequencies corresponding to the inner
solar system change as a result of mutual interactions of the inner planets
and the interactions of the inner planets with the outer gas giants (Laskar,
1990). As we will show in Section 4, a simple diffusion process is often a
remarkably accurate model that can explain most of the results seen in
ensembles of direct $N$-body integrations.

##  3 _N_ -body simulations

###  3.1 Methods

We run direct $N$-body simulations of the solar system which we will later
compare to our advection-diffusion model. All of our simulations for the solar
system (the Sun and eight planets) use exactly the same initial conditions
which are taken from NASA JPL Horizons data at the J2000 epoch. Even though
all the initial conditions are the same, because the systems are chaotic and
we use slightly different general relativistic corrections in each simulation,
the simulations diverge quickly. This is effectively the same as varying the
initial conditions of one planet by a tiny amount (Laskar & Gastineau, 2009).
We integrate simulations forward in time using REBOUND (Rein & Liu, 2012) and
the Wisdom-Holman integrator WHFast (Wisdom & Holman, 1992; Rein & Tamayo,
2015) with symplectic correctors and the lazy implementation of the kernel
method, WHCKL (Wisdom et al., 1996; Rein et al., 2019a). The WHCKL integrator
is well suited to provide highly accurate results for secularly evolving
systems (Rein et al., 2019b). We used a fixed timestep of
$dt=\sqrt{11}\,\mathrm{days}\approx 3.317\,\mathrm{days}$.

The aim of our simulations is not to exactly reproduce the solar system, but
to have a well defined model that can reproduce the most important dynamical
effects. Specifically, we do not consider the stellar evolution of the Sun
even though we integrate beyond the end of its life on the main sequence. We
also neglect the effects resulting from the solar oblateness, moons,
asteroids, tides, and other non-gravitational effects.

###  3.2 Perihelion precession from general relativity

Of all the contributions to Mercury’s perihelion precession, the most
important are the gravitational interactions from the other solar bodies
contributing $g_{1,\rm planets}=5.323\,^{\prime\prime}\mathrm{yr}^{-1}$
(arcseconds per year), followed by the precession from GR (the gravitoelectric
effect) with $g_{1,\rm gr}=0.4298\,^{\prime\prime}\mathrm{yr}^{-1}$, followed
by other less significant contributors such as solar oblateness which provide
$g_{1,\rm{other}}=2.8\times 10^{-4}\,{}^{\prime\prime}\mathrm{yr}^{-1}$ (Park
et al., 2017). Thus, the impact of general relativity follows closest behind
Newtonian planet-planet interactions and is more than three orders of
magnitude more influential than any additional effects. As discussed in
Section 2, the long-term stability of solar system is connected to Mercury’s
perihelion precession rate and thus the general relativistic effects are very
important (Laskar & Gastineau, 2009).

In this paper, we consider an $N$-body Newtonian model of the solar system. To
have a perihelion precession that is consistent with general relativity, we
include an additional potential term in our force calculation:

$\Phi_{\rm gr}=\alpha(t,\tau)\frac{6G^{2}M^{2}}{c^{2}r^{2}}$ (1)

where $r$ is the distance of the planet to the Sun and $M$, $G$, $c$ are the
mass of the Sun, the gravitational constant, and the speed of light
respectively. The perihelion precession rate caused by $\Phi_{\rm gr}$ is

$\dot{\varpi}_{{}_{\mathrm{GR}}}=\alpha(t,\tau)\frac{3GMn}{ac^{2}(1-e^{2})}$
(2)

where $n=\sqrt{GM/a^{3}}$ is the mean motion of a planet with semi-major axis
$a$ and $e$ is the planet’s eccentricity.

The time dependent parameter $\alpha(t,\tau)$ in the above equations enables
us to experiment with different strengths for the general relativistic
corrections. With $\alpha=1$, we refer to this potential as the standard first
order post-Newtonian corrections. This assumes that the Sun is the only body
in solar system massive enough so that particles orbiting it experience
general relativistic precession. This commonly used model gets the general
relativistic precession frequency right, but at the expense of introducing an
error on the mean motions on the order of $\mathcal{O}(GM/ac^{2})$ (Nobili &
Roxburgh, 1986; Saha & Tremaine, 1994). For the discussion in this paper, the
accuracy in the mean motion does not affect the results because of the
relative importance of Mercury’s perihelion precession rate to its stability
compared to its mean motion (there are no mean motion resonances in the inner
solar system).

With $\alpha=0$, we recover the purely classical regime with no general
relativistic precession. Although this is clearly inconsistent with many
observations that confirm general relativity, it can act as a useful case to
study if we want to understand the dynamics of the solar system. In previous
work, Laskar & Gastineau (2009) integrated thousands of solar system
simulations for billions of years and showed that without general relativistic
corrections, 60 per cent of the solutions resulted in an unstable solar system
after 5 Gyr (defined as Mercury having an eccentricity beyond 0.9).
Conversely, when general relativistic corrections were included, only 1 per
cent of the solutions were unstable (Laskar & Gastineau, 2009). See also Fig.
2 below.

In this paper, we go beyond simply turning general relativistic precession on
or off and allow for a time dependent parameter $\alpha$ defined as

$\alpha(t,\tau)=\begin{cases}1&t<0\\\ 1-(t/\tau)&0\leq t\leq\tau\\\
0&t>\tau\end{cases}\,.$ (3)

At the beginning of a simulation, $t=0$, the full GR corrections are taken
into account, consistent with present day observations. After $t=\tau$ no more
general relativistic corrections are applied and the simulation evolves
classically. We include the general relativistic corrections using REBOUNDx
(Tamayo et al., 2020) and the gr_potential module, modified to allow for the
time varying form given in equation (1).

###  3.3 _N_ -body results

Figure 1: The instability time of simulations in our ensemble that went
unstable as a function of the control parameter $\tau$ which determines the
time when general relativistic corrections to the perihelion precession are
turned off. Simulations are significantly more likely to go unstable before
12.5 Gyr when no general relativistic corrections are present (upper left
portion). The four different colours show the binning made for comparison to
the model in Section 5 and Fig. 4.

In this section we present the stability results of an ensemble of 1280 long-
term integrations of the solar system. We decrease the coupling strength,
$\alpha(t,\tau)$, of the general relativistic corrections to perihelion
precession linearly from the current value, 1, to 0 as defined in equation
(3). For each simulation we assign a different value for $\tau$. We sample
$\tau$ log-uniformly between $0.1$ Gyrs and $12.5$ Gyrs. We then integrate the
simulations for $12.5$ Gyrs, or until a collision or escape event occurs. We
do not carry out any simulations beyond the first physical collision or
escape.

We find that $1072$ out of $1280$ or 83.8 per cent of the simulations end up
in an instability, all involving a Mercury-Venus close encounter or collision.
175 simulations go unstable while some fractional GR precession is still
present. We present the instability times with respect to $\tau$ in Fig. 1.
The instability time is shown on the $y$-axis, $\tau$ is shown on the
$x$-axis. The upper-most diagonal line is when $\alpha(t,\tau)=0$ at $t=\tau$.
The remaining diagonal lines indicate the times when $\alpha(t,\tau)$ reaches
$0.25,0.5,0.75,$ and $0.95$ (from top to bottom).

##  4 Advection-diffusion model

Figure 2: A comparison showing the implementation of our model against the
previous work of Laskar & Gastineau (2009). The mean and $2\sigma$
uncertainties of the model is in good agreement with these numerical results
in both cases with and without GR. The dashed curve shows the effect of
changing the initial starting location of $g_{1}$ by removing $g_{1,\rm gr}$
from the GR model without altering the diffusion coefficient. The fact that
this curve does not match the numerical experiments without GR reveals that
the initial proximity between $g_{1}$ and $g_{5}$ alone does not account for
the increase in instabilities.

###  4.1 Fokker-Planck equation

In this section, we develop a simple model that is physically motivated and
can explain results from our new $N$-body simulations as well as the results
from previous studies that have both included and excluded general
relativistic corrections.

Our model describes the evolution of Mercury’s precession frequency $g_{1}$ as
an advection-diffusion process. This can be seen as a natural extension of
model used by Mogavero & Laskar (2021) who use a stochastic Wiener process to
describe the evolution of Mercury’s precession frequency $g_{1}$ (see their
section 8.2).

We use the Fokker-Planck equation to describe the evolution of
$p(g_{1},t,\tau)$, the probability density of the secular frequency $g_{1}$ as
a function of time $t$ and our model parameter $\tau$:

$\frac{\partial p(g_{1},t,\tau)}{\partial t}=-\mu(t,\tau)\frac{\partial
p(g_{1},t,\tau)}{\partial
g_{1}}+D(t,\tau)\frac{\partial^{2}p(g_{1},t,\tau)}{\partial g_{1}^{2}}.$ (4)

Here $\mu(t,\tau)$ describes the advection and $D(t,\tau)$ the diffusion of
$g_{1}$. We use the current precession frequency of Mercury,
$g_{1,0}=5.577^{\prime\prime}\mathrm{yr}^{-1}$, as the initial condition:

$p(g_{1},0,\tau)=\delta(g_{1}-g_{1,0}).$ (5)

There is an upper boundary which acts as a reflecting barrier at
$g_{1,\mathrm{max}}=5.72367^{\prime\prime}\mathrm{yr}^{-1}$ (Mogavero &
Laskar, 2021). We also have a lower boundary where $g_{1}$ overlaps with
$g_{5}=4.257^{\prime\prime}\mathrm{yr}^{-1}$, the precession frequency of
Jupiter. When this lower boundary is hit, the system is assumed to go unstable
on a very short timescale because of the $g_{1}-g_{5}$ secular resonance.

###  4.2 Model with general relativistic precession

In the limit where general relativistic corrections are present during the
entire simulation, Mogavero & Laskar (2021) show that the results of high-
order secular simulations can be closely matched with a constant diffusion
coefficient $D_{\rm gr}$ and a vanishing convection term $\mu=0$.

A more intuitive understanding of these quantities can be achieved by
rewriting the diffusion coefficient $D$ and the initial precession frequency
$g_{1,0}$ into a diffusion timescale $T_{D}$ and a dimensionless parameter
$\beta$ that describes how close the initial condition for $g_{1}$ is to the
lower boundary $g_{5}$ where we expect an instability to occur. Making use of
$(g_{1,\mathrm{max}}-g_{5})$, the distance between the two boundaries, we can
define the characteristic diffusion timescale as

$T_{D}=\frac{(g_{1,\mathrm{max}}-g_{5})^{2}}{4D}.$ (6)

and the parameter111This parameter is equivalent to the parameter $\alpha$ in
Mogavero & Laskar (2021). $\beta$ that describes our initial conditions as:

$\beta=(g_{1,0}-g_{5})/(g_{1,\mathrm{max}}-g_{5}).$ (7)

Using $g_{5}=4.257^{\prime\prime}\mathrm{yr}^{-1}$ we have $\beta=0.9$.
Mogavero & Laskar (2021) show that a value of $T_{D,\rm
gr}=27.6\,\mathrm{Gyr}$ reproduces the statistical results from their
simulations (which include general relativistic corrections) remarkably well.
The survival fraction as a function of time corresponding to this model is the
uppermost line plotted in Fig. 2. $N$-body data (an ensemble of 2501
simulations) from the previous work by Laskar & Gastineau (2009) is overlaid
as red squares for comparison. It shows how this advection-diffusion model
captures the instability rates of Mercury in $N$-body solar system simulations
for the first 5 Gyr very well. Obtaining a meaningful number of instabilities
beyond the remaining lifetime of the sun is costly, but we show below that
this model can also be used to give accurate results well beyond 5 Gyr.

###  4.3 Model without general relativistic precession

We can repeat the analysis from the previous section for the case without
general relativistic precession. In our notation, this is the case where
$\tau\rightarrow 0$ and $\alpha=0$. The contribution from general relativity
to the precession frequency of Mercury is $g_{1,\rm
gr}=0.4298^{\prime\prime}\mathrm{yr}^{-1}$. Without general relativity, our
initial condition is therefore shifted to

$p(g_{1},0,0)=\delta(g_{1}-(g_{1,0}-g_{1,\rm gr})),$ (8)

or equivalently $\beta=0.607$. If we assume this is the only change to the
system, then we end up with a survival fraction shown as a dashed line in Fig.
2. Clearly, although the rate of instability is higher, this does not fit the
numerical experiments by Laskar & Gastineau (2009).

This leads to the conclusion that we not only need to change the initial
conditions, but also the diffusion coefficient. To find the diffusion
coefficient we can make use of the $N$-body results by Laskar & Gastineau
(2009) who find that 60 per cent of all simulations have gone unstable after 5
Gyr, using an ensemble of 201 simulations. Choosing $T_{D,\rm
nogr}=3.22\,\mathrm{Gyr}$ results in a model that fits those results well, as
can be seen in Fig. 2. While there are some discrepancies between the model
and $N$-body results, especially in the 1-2 Gyr regime, we attribute this
discrepancy to the relatively small number of simulations used in the
ensemble.

###  4.4 Model with time dependent general relativistic precession

We are now in a position to describe a model for our case where we have a time
dependent $\alpha(t,\tau)$ given by equation (3), leading to a time dependent
general relativistic precession. In this paper, $\tau$ is the model parameter
that describes the timescale on which the general relativistic precession
changes.

For our model, the initial conditions are the present day precession frequency
$g_{1,0}$, i.e.:

$p(g_{1},0,\tau)=\delta(g_{1}-g_{1,0}).$ (9)

The diffusion coefficient is no longer a constant, but a function of time

$D(t,\tau)=\begin{cases}\left(\sqrt{D_{\rm
gr}}\cdot\left(1-\frac{t}{\tau}\right)+\sqrt{D_{\rm
nogr}}\cdot\frac{t}{\tau}\right)^{2}&t\leq\tau\\\ D_{\rm
nogr}&t>\tau\end{cases}$ (10)

We chose this specific form for the diffusion coefficient because the
interpolation between $D_{\rm gr}$ and $D_{\rm nogr}$ corresponds to a simple
linear interpolation in the forcing strength222Note that $D\sim t_{c}\cdot
F^{2}$ where $t_{c}$ is the characteristic forcing timescale and $F$ is the
characteristic forcing strength in a random walk process, see e.g. Rein &
Papaloizou (2009)..

In addition, we also now include the advection term in the Fokker-Planck
equation. Specifically:

$\mu(t,\tau)=\begin{cases}g_{1,\rm gr}/\tau&t\leq\tau\\\ 0&t>\tau\end{cases}.$
(11)

If we ignore the diffusion, then this form of $\mu$ moves the precession
frequency from the initial value of $g_{1,0}$ to $g_{1,0}-g_{1,gr}$ by time
$t=\tau$, i.e. when $\alpha=0$ and there is no more general relativistic
precession.

Note that there are no free parameters in this extension of the model. All the
parameters, most importantly $D_{\rm gr}$ and $D_{\rm nogr}$ are already fixed
by matching the simulations of Laskar & Gastineau (2009).

##  5 Results

Figure 3: The survival fraction of the entire ensemble of simulations compared
to the survival given by the model described in Section 4. The mean and
$2\sigma$ uncertainties are shown for the model.

Figure 4: The survival fraction of the ensemble of $N$-body simulations
separated into four bins of $\tau$ as shown in Fig. 1 is given in black. The
mean survival of the model described in Section 4 is shown for each bin with
$2\sigma$ uncertainties. The dotted vertical lines indicate when all of the
$N$-body simulations in the bin have completely decoupled from general
relativistic corrections. Thus, from $t=0$ to the dotted line, GR corrections
are gradually being turned off.

The survival fraction for our entire ensemble compared to the model is
presented in Fig. 3. It shows the results from our ensemble of 1280 $N$-body
simulations compared against the advection-diffusion model with time dependent
coefficients. The mean and $2\sigma$ uncertainties expected for an ensemble
size of 1280 are shown for the model.

For Fig. 4 we split our ensemble of simulations into four equally sized groups
to show the dependence on $\tau$. The colours match those used in Fig. 1. Each
panel includes the observed numerical survival fraction in our $N$-body
simulations along with the expected mean survival fraction given by our model
and the corresponding $2\sigma$ confidence interval for an ensemble size of
320. The plots also include a dotted vertical line indicating when the last
simulation in the bin has reached $\alpha(t,\tau)=0$. A notable feature in
some of the panels is that the model produces instabilities slightly earlier
than in the $N$-body experiments (starting around 4.5 Gyr). Even so, the model
and the numerical simulations agree remarkably well. In particular, note that
our model is not fit to the data, but rather uses a physically motivated model
and parameters which are solely determined by the limiting cases done in
previous work by Laskar & Gastineau (2009) and Mogavero & Laskar (2021).

The far left panel in Fig. 4 shows the survival fraction of simulations after
removing the GR precession quickly333Quickly here means that the simulations
in this bin all have the general relativistic corrections removed on a time
scale shorter than the lower bound on the destabilization time of Mercury
given by Mogavero & Laskar (2021) $\sim 0.56$ Gyr.. The results match those by
Laskar & Gastineau (2009): $\sim 60$ per cent go unstable within 5 Gyr. Note
that removing GR precession from these integrations (even in the far future)
does not immediately cause instability. In fact, some simulations continue to
remain stable (without GR) for more than $12$ Gyr.

##  6 Conclusions

In this paper, we explored how sensitive the stability of the solar system is
to changes in the perihelion precession caused by general relativistic
corrections. Earlier work by Laskar & Gastineau (2009) has shown that turning
off general relativistic corrections can increase the fraction of instability
within 5 Gyr by a factor of 60. We develop a simple advection-diffusion model
that can explain the results of previous $N$-body experiments which either
include or exclude general relativistic precession, thus extending the model
of Mogavero & Laskar (2021). We then go beyond simply turning the GR
precession on or off and instead vary the rate of GR precession smoothly over
time. We do this not because we think this is something that will occur in the
real solar system, but because it offers a way to experiment with this
dynamical system in a controlled way.

We show that our advection-diffusion model can naturally explain the case of a
time varying general relativistic precession remarkably well without the
introduction of any new free parameters. We are also confident that this model
provides a sound extrapolation for the time to Mercury’s dynamical instability
beyond 5 Gyr (neglecting changes to the stellar life cycle of the Sun and
other external factors). This confidence is based on the work of others
mentioned throughout this paper as well as the addition of our $N$-body
ensemble extending all the way to 12.5 Gyr (to date the largest set of solar
system integrations extended to this time). The model shows strong agreement
to $N$-body data for so many different regimes of general relativistic
precession. Additionally, the model reproduces the statistical results across
different analytical and numerical methods, again, without the introduction of
any new free parameters.

Without general relativistic corrections, the precession frequency of Mercury,
$g_{1}$, is closer to that of Jupiter, $g_{5}$, by about
$0.4298^{\prime\prime}\mathrm{yr}^{-1}$. However, our results show that this
alone does not explain the increased rate of instability observed in $N$-body
experiments. We show that our advection-diffusion model can represent the
instability rate only if we significantly decrease the diffusion timescale (or
increase the diffusion coefficient). In other words, although the
$g_{1}-g_{5}$ resonance determines when Mercury’s eccentricity reaches a
critical value, the process to get to this resonance is much faster without GR
corrections.

This result might not be surprising, given that the solar system is a complex
chaotic system with way more degrees of freedom than our model (which just
allows for a diffusion in $g_{1}$). Our simulations with time-varying general
relativistic precession allow us to gain some more insight into this process.
As we’ve shown in Section 5, our model matches the $N$-body experiments very
well given there are no free parameters aside from those used to match the
limiting cases of GR being completely on or off. This leads to several
conclusions. First, there is no critical value of the general relativistic
precession rate that needs to be crossed in order to significantly change the
instability rate. Second, because this is such a smooth transition all the way
from one limit to the other, we can rule out a process that involves a
sweeping secular resonance. If this were not so smooth, we would see an
increase in the instability rate at some critical GR precession frequency.
This makes sense because the diffusion in $g_{1}$ is strong enough so that a
small change due to the GR precession has a negligent effect. This provides
additional evidence that the stability of the solar system is robust to
moderate changes to the secular system (see also Laskar, 1990; Brown & Rein,
2022). The fact that we are able to model the evolution so well with a simple
diffusion model also shows that current numerical results (Laskar & Gastineau,
2009; Zeebe, 2015; Abbot et al., 2021) are robust against small perturbations
whether they are physical or numerical (Abbot et al., 2022). We expect that
statistical results are in agreement as long as simulations resolve secular
frequencies accurately enough so that physical diffusion (not numerical
diffusion or advection) is the dominant driver of the instability.

In summary, we reconfirm that the solar system’s evolution is well described
by a simple (advection-)diffusion process, even in the presence of other small
perturbations. Here we considered the perturbation to be the change in the
general relativistic precession frequency, but any other physical effect that
changes the planets’ precession frequencies will lead to a similar result. The
initial motivation of this project was to use the stability of the solar
system as a test-bed for additional physics that are typically not included in
simulations of the solar system. For example, had our results shown that
slowly changing the general relativistic precession frequency leads to
significant increase in the instability rate, then we could have placed a
limit on alternative theories of general relativity by noting that the solar
system has not gone unstable yet. Clearly this is not the case. For future
work, we propose further investigation into the expected time to instability
for other secularly evolving planetary systems and whether or not a one
dimensional advection-diffusion model is representative of more systems or
only the solar system (see e.g. Hussain & Tamayo, 2020).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## Code availability

A repository containing a portion of the data underlying this article and code
for running the simulations and generating the figures can be found at
github.com/zyrxvo/GR-and-Long-term-Stability.

## ORCID iDs

Garett Brown  ID https://orcid.org/0000-0002-9354-3551 Hanno Rein  ID
https://orcid.org/0000-0003-1927-731X

## Acknowledgments

We are very grateful to an anonymous referee for a helpful review that
improved the quality of this paper. We would like to thank Scott Tremaine, Sam
Hadden, Dang Pham, and Mykhaylo Plotnykov for useful discussions. This
research has been supported by the NSERC Discovery Grants RGPIN-2014-04553 and
RGPIN-2020-04513. This research was made possible by the open-source projects
Jupyter (Kluyver et al., 2016), iPython (Pérez & Granger, 2007), matplotlib
(Hunter, 2007; Droettboom et al., 2016), and GNU Parallel (Tange, 2022). This
research was enabled in part by support provided by Digital Research Alliance
of Canada (formerly Compute Canada) (alliancecan.ca). Computations were
performed on the Niagara supercomputer (Loken et al., 2010; Ponce et al.,
2019) at the SciNet HPC Consortium (scinethpc.ca). SciNet is funded by: the
Canada Foundation for Innovation; the Government of Ontario; Ontario Research
Fund - Research Excellence; and the University of Toronto.

## References

  * Abbot et al. (2022) Abbot, D. S., Hernandez, D. M., Hadden, S., Webber, R. J., Afentakis, G. P., & Weare, J. 2022, arXiv e-prints, arXiv:2212.14844
  * Abbot et al. (2021) Abbot, D. S., Webber, R. J., Hadden, S., Seligman, D., & Weare, J. 2021, ApJ, 923, 236
  * Arnol’d (1963) Arnol’d, V. I. 1963, Russian Mathematical Surveys, 18, 9
  * Batygin & Laughlin (2008) Batygin, K. & Laughlin, G. 2008, The Astrophysical Journal, 683, 1207
  * Batygin et al. (2015) Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120
  * Brown & Rein (2022) Brown, G. & Rein, H. 2022, Monthly Notices of the Royal Astronomical Society, 515, 5942
  * Chan & Lee (2023) Chan, M. H. & Lee, C. M. 2023, MNRAS, 518, 6238
  * Droettboom et al. (2016) Droettboom, M., Hunter, J., Caswell, T. A., Firing, E., Nielsen, J. H., Elson, P., Root, B., Dale, D., Lee, J.-J., Seppänen, J. K., McDougall, D., Straw, A., May, R., Varoquaux, N., Yu, T. S., Ma, E., Moad, C., Silvester, S., Gohlke, C., Würtz, P., Hisch, T., Ariza, F., Cimarron, Thomas, I., Evans, J., Ivanov, P., Whitaker, J., Hobson, P., mdehoon, & Giuca, M. 2016, matplotlib: matplotlib v1.5.1
  * Eddington (1923) Eddington, A. S. 1923, The mathematical theory of relativity (The University Press)
  * Einstein (1915) Einstein, A. 1915, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften, 844
  * Hoang et al. (2022) Hoang, N. H., Mogavero, F., & Laskar, J. 2022, MNRAS, 514, 1342
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Hussain & Tamayo (2020) Hussain, N. & Tamayo, D. 2020, MNRAS, 491, 5258
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., Granger, B., Bussonnier, M., Frederic, J., Kelley, K., Hamrick, J., Grout, J., Corlay, S., et al. 2016, Positioning and Power in Academic Publishing: Players, Agents and Agendas, 87
  * Kolmogorov (1954) Kolmogorov, A. N. 1954, in Dokl. Akad. Nauk SSSR, Vol. 98, 527–530
  * Lagrange (1776) Lagrange, J. 1776, Sur l’altération des moyens mouvements des planètes
  * Lagrange (1778) Lagrange, J. L. 1778, Recherches sur les équations séculaires des mouvemens des noeuds, et des inclinaisons des orbites des planètes (de l’Imprimerie Royale)
  * Lagrange (1781) —. 1781, Théorie des variations séculaires des éléments des planètes, Première partie (Nouveaux Mémoires de l’Académie des Sciences et Belles-Lettres de Berlin)
  * Lagrange (1782) —. 1782, Théorie des variations séculaires des éléments des planètes, Seconde partie contenant la détermination de ces variations pour chacune des plaètes principales (Nouveaux Mémoires de l’Académie des Sciences et Belles-Lettres de Berlin)
  * Lagrange (1783a) —. 1783a, Sur les variations séculaires des mouvements moyens des planètes (Nouveaux Mémoires de l’Académie des Sciences et Belles-Lettres de Berlin)
  * Lagrange (1783b) —. 1783b, Théorie des variations périodiques des mouvements des planètes, Première partie (Nouveaux Mémoires de l’Académie des Sciences et Belles-Lettres de Berlin)
  * Lagrange (1784) —. 1784, Théorie des variations périodiques des mouvements des planètes, Seconde partie (Nouveaux Mémoires de l’Académie des Sciences et Belles-Lettres de Berlin)
  * Laplace (1775) Laplace, P. S. 1775, Mémoire de l’Académie royale des sciences de Paris
  * Laplace (1776) —. 1776, Mémoire de l’Académie royale des sciences de Paris, 8, 199
  * Laskar (1985) Laskar, J. 1985, A&A, 144, 133
  * Laskar (1986) —. 1986, A&A, 157, 59
  * Laskar (1989) —. 1989, Nat, 338, 237
  * Laskar (1990) Laskar, J. 1990, Icarus, 88, 266
  * Laskar (2000) —. 2000, Physical Review Letters, 84, 3240
  * Laskar (2013) —. 2013, in Chaos (Springer), 239–270
  * Laskar & Gastineau (2009) Laskar, J. & Gastineau, M. 2009, Nat, 459, 817
  * Laskar et al. (1992) Laskar, J., Quinn, T., & Tremaine, S. 1992, Icarus, 95, 148
  * Le Verrier (1840) Le Verrier, U. 1840, J. Math. Pures Appli, 4, 220
  * Le Verrier (1841) —. 1841, J. Math. Pures Appli, 28
  * Lithwick & Wu (2011) Lithwick, Y. & Wu, Y. 2011, ApJ, 739, 31
  * Loken et al. (2010) Loken, C., Gruner, D., Groer, L., Peltier, R., Bunn, N., Craig, M., Henriques, T., Dempsey, J., Yu, C.-H., Chen, J., Dursi, L. J., Chong, J., Northrup, S., Pinto, J., Knecht, N., & Zon, R. V. 2010, Journal of Physics: Conference Series, 256, 012026
  * Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
  * Moffat & Toth (2008) Moffat, J. W. & Toth, V. T. 2008, ApJ, 680, 1158
  * Mogavero & Laskar (2021) Mogavero, F. & Laskar, J. 2021, A&A, 655, A1
  * Mogavero & Laskar (2022) —. 2022, arXiv e-prints, arXiv:2205.03298
  * Möser (1962) Möser, J. 1962, Nachr. Akad. Wiss. Göttingen, II, 1
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics (Cambridge university press)
  * Newton (1687) Newton, I. 1687, Philosophiae Naturalis Principia Mathematica (Societatis Regiae)
  * Nobili & Roxburgh (1986) Nobili, A. & Roxburgh, I. W. 1986, in IAU Symposium, Vol. 114, Relativity in Celestial Mechanics and Astrometry. High Precision Dynamical Theories and Observational Verifications, ed. J. Kovalevsky & V. A. Brumberg, 105–110
  * Park et al. (2017) Park, R. S., Folkner, W. M., Konopliv, A. S., Williams, J. G., Smith, D. E., & Zuber, M. T. 2017, The Astronomical Journal, 153, 121
  * Pérez & Granger (2007) Pérez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21
  * Poincaré (1898) Poincaré, H. 1898, Annuaire du Bureau des Longitudes pour l’an 1898
  * Poincaré (1899) —. 1899, Les méthodes nouvelles de la mécanique céleste, Vol. 3 (Gauthier-Villars et fils, imprimeurs-libraires)
  * Ponce et al. (2019) Ponce, M., van Zon, R., Northrup, S., Gruner, D., Chen, J., Ertinaz, F., Fedoseev, A., Groer, L., Mao, F., Mundim, B. C., et al. 2019, in Proceedings of the Practice and Experience in Advanced Research Computing on Rise of the Machines (learning), ACM, 34
  * Rein & Liu (2012) Rein, H. & Liu, S. F. 2012, A&A, 537, A128
  * Rein & Papaloizou (2009) Rein, H. & Papaloizou, J. C. B. 2009, A&A, 497, 595
  * Rein & Tamayo (2015) Rein, H. & Tamayo, D. 2015, MNRAS, 452, 376
  * Rein et al. (2019a) Rein, H., Tamayo, D., & Brown, G. 2019a, Monthly Notices of the Royal Astronomical Society
  * Rein et al. (2019b) —. 2019b, Monthly Notices of the Royal Astronomical Society, 489, 4632
  * Saha & Tremaine (1994) Saha, P. & Tremaine, S. 1994, AJ, 108, 1962
  * Sussman & Wisdom (1992) Sussman, G. J. & Wisdom, J. 1992, Science, 257, 56
  * Tamayo et al. (2020) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2020, Monthly Notices of the Royal Astronomical Society, 491, 2885
  * Tange (2022) Tange, O. 2022, GNU Parallel 20221022 (‘Nord Stream’), GNU Parallel is a general parallelizer to run multiple serial command line programs in parallel without changing them.
  * Verlinde (2017) Verlinde, E. 2017, SciPost Physics, 2, 016
  * Wisdom & Holman (1992) Wisdom, J. & Holman, M. 1992, AJ, 104, 2022
  * Wisdom et al. (1996) Wisdom, J., Holman, M., & Touma, J. 1996, Fields Institute Communications, Vol. 10, p. 217, 10, 217
  * Zakamska & Tremaine (2004) Zakamska, N. L. & Tremaine, S. 2004, The Astronomical Journal, 128, 869
  * Zeebe (2015) Zeebe, R. E. 2015, ApJ, 798, 8


$a_{1},a_{2}\in{\mathfrak{A}}$ ($a_{1},a_{2}\in{\mathfrak{D}}_{s}$) realize an
extended 2-type $(\beta,j)$ in ${\mathfrak{A}}$ (resp. ${\mathfrak{D}}_{s}$)
if (1) $a_{1},a_{2}$ realize 2-outer-type $\beta$, and (2) either $i=\bot$, or
$\mathfrak{f}_{j}(a_{1})=a_{2}$, or $\mathfrak{f}_{j}(a_{2})=a_{1}$. Denote by
$\mbox{\large\boldmath$\beta$}^{+}$ the set of all extended 2-types realized
in ${\mathfrak{A}}$. For each $i\in\mathcal{I}_{5}$ select a pattern function
$\mathfrak{pat}^{{\mathfrak{A}}}_{i}:\mbox{\large\boldmath$\beta$}^{+}\to
A^{3}$ such that
$\mathfrak{pat}^{{\mathfrak{A}}}_{i}(\beta,j)=(a_{1},a_{2},a_{3})$ if
$\\{a_{1},a_{2}\\}$ realizes $(\beta,j)$ in ${\mathfrak{A}}$, and
${\mathfrak{A}}\models\varphi_{i}(a_{1},a_{2},a_{3}).$

_Stage 3: Constructing ${\mathfrak{D}}_{s+1}$ from ${\mathfrak{D}}_{s}$._ Let
$b_{1},b_{2}\in{\mathfrak{D}}_{s}$ be a pair of distinct elements realizing an
extended 2-type $(\beta,j)$ for which (c) or (c’) fails for some
$i\in\mathcal{I}_{5}$. Let
$(a_{1},a_{2},a_{3})=\mathfrak{pat}^{{\mathfrak{A}}}_{i}(\beta,j)$. Set
$({\mathfrak{D}}_{s+1}\\!\\!\restriction\\!\\!\\{b_{1},b_{2},b_{3}\\})\leftarrow
pull({\mathfrak{A}},b_{1}\leftarrow a_{1},b_{2}\leftarrow
a_{2},b_{3}\leftarrow a_{3}),$ where $b_{3}\in B\times\\{s+1\\}$ is any
element of the freshly added copy of ${\mathfrak{B}}$ with the same 1-type as
$a_{3}$. After providing a witness for $b_{1},b_{2}$ this way, we need to
complete ${\mathfrak{D}}$, i.e., we need to specify the relations between all
pairs and triples of elements for which we have not done it yet. Take any
subset $S\subseteq D_{s+1}$ of size $2$ or $3$, and search for an appropriate
pattern in ${\mathfrak{A}}$, then pull it into
${\mathfrak{D}}_{s+1}\\!\\!\restriction\\!\\!S$. We will omit the details
here. Consult the Section 4.1 and the algorithm defined there for the very
similar line of the reasoning. Therefore,
${\mathfrak{D}}_{s+1}\models\varphi_{1}$ (Condition (a)), and
${\mathfrak{D}}_{s+1}\models\varphi_{i}(d,\mathfrak{f}_{j}(d),d^{\prime})$ for
each $j\in\mathcal{I}_{3}$ (Condition (b)).

_Stage 4: Finalization._ Conditions (a) and (b) are satisfied after every
step, and after at most $|\mbox{\large\boldmath$\beta$}^{+}|$ many steps (the
number $k$), which is polynomial in $|\mbox{\large\boldmath$\beta$}|$ and
$|\varphi|$, all the $2$-types will satisfy (c) and (c’). Thus, all the
conditions are met for ${\mathfrak{D}}_{k}$. The size of the domain of
${\mathfrak{D}}_{k}$ is at most exponential in $|\varphi|$, since
${\mathfrak{D}}_{k}$ consists in fact of $|\mbox{\large\boldmath$\beta$}^{+}|$
copies of ${\mathfrak{B}}$, and both ${\mathfrak{B}}$ and
$|\mbox{\large\boldmath$\beta$}^{+}|$ are exponentially bounded in $\varphi$.

Let us now turn to the right-to-left direction of Lemma 12.

Right-to-left: building a model from a witness of satisfiability

Let ${\mathfrak{D}}$ be a witnesses of satisfiability. If $\varphi$ does not
have any $\mathcal{I}_{5}$-conjuncts then ${\mathfrak{D}}$ is already a model
of $\varphi$. In the other case, we will show how to construct a new witness
of satisfiability ${\mathfrak{D}}^{\prime}$, but bigger. Iterating this
construction indefinitely, and taking the limit structure, we will obtain an
infinite model of ${\mathfrak{A}}$.

The construction of a new witnesses of satisfiability is very similar to the
process of extracting a witnesses of satisfiability from a model. Hence, we
will just sketch the construction, and stress the differences. The domain of
${\mathfrak{D}}^{\prime}$ is the set333By $\binom{A}{2}$, where $A$ is a set,
we mean the set of all $2$-element subsets of $A$.

$D^{\prime}=D\times\\{(\bot,\emptyset)\\}\cup
D\times\mathcal{I}_{5}\times{\binom{D}{2}},$

and let the structure on
${\mathfrak{D}}^{\prime}\\!\\!\restriction\\!\\!D\times\\{s\\}$ be isomorphic
to ${\mathfrak{D}}$ for each
$s\in\\{(\bot,\emptyset)\\}\cup(\mathcal{I}_{5}\times\binom{D}{2})$. Then we
take care of witnesses for all $\mathcal{I}_{5}$-conjuncts for all the pairs
of elements in ${\mathfrak{D}}_{\bot,\emptyset}$. The witness for a conjunct
$\varphi_{i}$ and a pair $\\{d_{1},d_{2}\\}$ is provided in
$D_{i,\\{d_{1},d_{2}\\}}$. The major difference is that this time we pull
outer-substructures not from an existing model of $\varphi$, but rather from
${\mathfrak{D}}$ itself; this is always possible due to Conditions (c) and
(c’). Then we complete the structure, similarly as in our previous
constructions, by setting the relations between all the pairs and the triples
of elements for which we have not done it yet. Again, in this completion
process the outer-substructures are pulled from ${\mathfrak{D}}$, and not from
an external structure. So after that all the $\mathcal{I}_{2}$\- and
$\mathcal{I}_{3}$-conjuncts are satisfied, and all the pairs of elements from
$D$ have witnesses for all $\mathcal{I}_{5}$-conjuncts. We claim that the
newly obtained structure is again a witness of satisfiability, with all
elements from $D$ having the required witnesses, but with a strictly larger
domain. This allows us to continue the process ad infinitum, with the natural
limit structure being a model of $\varphi$.

Having shown Lemma 12, and noting that checking if a given structure
${\mathfrak{D}}$ is a witness of satisfiability for $\varphi$ can be easily
done in time polynomial in $|D|$ and $|\varphi|$, we get:

######  Theorem 13.

The satisfiability problem for $\mbox{\rm FO}^{3}_{-}$ is NExpTime-complete.

Note that in this subsection we says nothing about finite satisfiability for
$\mbox{\rm FO}^{3}_{-}$ or about its finite model property. These is left for
the next subsection.

###  5.2 Finite model property via probabilistic method

Here we show show that $\mbox{\rm FO}^{3}_{-}$ has the finite model property.
As in Section 3.3 we employ the probabilistic method. This time we will use as
finite building blocks the structures ${\mathfrak{B}}\models\varphi_{1}$
constructed deterministically. Randomness will come to play when those blocks
are joined.

Let $\varphi=\varphi_{1}\wedge\varphi_{2}$ be a normal form $\mbox{\rm
FO}^{3}_{-}$-sentence and ${\mathfrak{A}}$ its model. Wlog. we assume that any
$1$-type realized in ${\mathfrak{A}}$ is realized by at least three distinct
elements (if it is not the case, we replace ${\mathfrak{A}}$ by the structure
${\mathfrak{A}}^{\prime}$ with domain $A\times{\mathbb{N}}$, as in the
beginning of Section 3.3). As previously, we pick for every
$j\in\mathcal{I}_{3}$ a witness function $\mathfrak{f}_{j}:A\rightarrow A$
such that for every $a\in A$ we have ${\mathfrak{A}}\models\forall
z\varphi_{j}(a,\mathfrak{f}_{j}(a),z)$. Let
$\mbox{\large\boldmath$\alpha$}{},\mbox{\large\boldmath$\beta$}{}$ and
$\mbox{\large\boldmath$\gamma$}{}^{*}$ be the set of $1$-types realized in
${\mathfrak{A}}$, the set of $2$-types realized in ${\mathfrak{A}}$, and,
resp., the set of $3$-outer-types realized in ${\mathfrak{A}}$.

Let us construct a model ${\mathfrak{B}}\models\varphi_{1}$, together with
functions $\mathfrak{pat}:B\rightarrow A$ and $\mathfrak{f}_{j}:B\rightarrow
B$, for $j=1,\ldots,|\mathcal{I}_{3}|$, the latter without fixed points, as in
Section 4.1. We may assume that the set of $1$-types realized in
${\mathfrak{B}}$ is $\alpha$ (cf. Section 5.1, Left-to-right, Stage 1).

We now define a sequence of (partly) random structures ${\mathfrak{C}}_{n}$,
for $n=1,2,\ldots$.

Stage 0: The domain. $C_{n}=B\times\\{1,\ldots,n\\}$. We make the
substructures ${\mathfrak{C}}_{n}\\!\\!\restriction\\!\\!B\times\\{m\\}$
isomorphic to ${\mathfrak{B}}$, for all $m=1,2,\ldots,n$, and set no other
connections in this stage. We naturally lift functions $\mathfrak{pat}$ and
the $\mathfrak{f}_{j}$ to ${\mathfrak{C}}_{n}$.

Stage 1: $2$-types. For every pair $c_{1},c_{2}\in C_{n}$ of distinct elements
for which ${\rm tp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2}})$ has not yet been
defined (note that $c_{1}$ and $c_{2}$ belong then to different copies of
${\mathfrak{B}}$ in this case) proceed as follows. Let $\alpha_{1}={\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{1}})$, $\alpha_{2}={\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{2}})$ and let
$\mbox{\large\boldmath$\beta$}_{0}$ be the set of $2$-types $\beta$ realized
in ${\mathfrak{A}}$ such that $\beta\\!\\!\restriction\\!\\!x_{1}=\alpha_{1}$
and $\beta\\!\\!\restriction\\!\\!x_{2}=\alpha_{2}$. Clearly
$\mbox{\large\boldmath$\beta$}_{0}$ is non-empty, since even if
$\alpha_{1}=\alpha_{2}$ we will have two distinct realizations of this
$1$-type by our assumption about ${\mathfrak{A}}$. Choose ${\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2}})$ at random from
$\mbox{\large\boldmath$\beta$}_{0}$, with uniform probability.

Stage 2a: $3$-outer-types. For every tuple $c_{1},c_{2},c_{3}\in C_{n}$ such
that $c_{2}=\mathfrak{f}_{j}(c_{1})$ for some $j$, and ${\rm
pretp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2},c_{3}})$ has not yet been defined
proceed as follows. Recall that by the construction of ${\mathfrak{B}}$,
$c_{2}$ may not be in the image of any $\mathfrak{f}_{k}$ for $j\not=k$ in
this case. Let $a_{1}=\mathfrak{pat}(c_{1})$ and
$a_{2}=\mathfrak{f}_{j}(a_{1})$, and $\alpha_{3}={\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{3}})$. Let
$\mbox{\large\boldmath$\gamma$}^{*}_{0}$ be the set of $3$-outer-types
realized in ${\mathfrak{A}}$ by tuples $a_{1},a_{2},a_{3}$, for all $a_{3}\in
A$ having $1$-type $\alpha_{3}$. Note that
$\mbox{\large\boldmath$\gamma$}^{*}_{0}$ is non-empty–even if $a_{1}$ and
$a_{2}$ have type $\alpha_{3}$ then at least one appropriate $a_{3}$,
different from $a_{1}$, $a_{2}$ exists since we assumed that every $1$-type
realized in ${\mathfrak{A}}$ is realized at least three times. Choose ${\rm
pretp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2},c_{3}})$ at random from
$\mbox{\large\boldmath$\gamma$}^{*}_{0}$, with uniform probability.

Stage 2b: $3$-outer-types contd. For every tuple $c_{1},c_{2},c_{3}\in C_{n}$
for which ${\rm pretp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2},c_{3}})$ has not
yet been defined proceed as follows. Let $\alpha_{1}={\rm
tp}^{{{\mathfrak{C_{n}}}}}({a_{1}})$, $\alpha_{2}={\rm
tp}^{{{\mathfrak{C_{n}}}}}({a_{2}}),\alpha_{3}={\rm
tp}^{{{\mathfrak{C_{n}}}}}({a_{3}})$ and let
$\mbox{\large\boldmath$\gamma$}^{*}_{0}$ be the set of $3$-outer-types
$\gamma^{*}$ realized in ${\mathfrak{A}}$ such that
$\gamma^{*}\\!\\!\restriction\\!\\!x_{1}=\alpha_{1}$ and
$\gamma^{*}\\!\\!\restriction\\!\\!x_{2}=\alpha_{2}$ and
$\gamma^{*}\\!\\!\restriction\\!\\!x_{3}=\alpha_{3}$. Again, using the
assumption about ${\mathfrak{A}}$ we note that
$\mbox{\large\boldmath$\gamma$}^{*}_{0}$ is non-empty. Choose ${\rm
pretp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2},c_{3}})$ at random from
$\mbox{\large\boldmath$\gamma$}^{*}_{0}$, with uniform probability.

Let us now estimate the probability that
${\mathfrak{C}}_{n}\not\models\varphi$. It is not difficult to observe that,
regardless of our random choices, ${\mathfrak{C}}_{n}\models\varphi_{1}$.
Indeed, the conjuncts indexed by $\mathcal{I}_{3}$ are respected since they
are respected in ${\mathfrak{B}}$ and we ensure that $\mathfrak{f}_{j}(c)$
remains an appropriate witness for $c$ in ${\mathfrak{C}}_{n}$, for every $c$
in Stage 2a; the conjuncts indexed by $\mathcal{I}_{2}$ are respected since
they are respected in ${\mathfrak{B}}$ and all outer-types defined in
${\mathfrak{C}}$ are copied from ${\mathfrak{A}}$ (Stages 2a, 2b). Hence, we
need to take into consideration only the satisfaction of $\varphi_{2}$. Let
$F^{i}_{c_{1},c_{2}}$, for a pair of distinct elements $c_{1},c_{2}\in C_{n}$
and $i\in\mathcal{I}_{5}$, denotes the event
“${\mathfrak{C}}_{n}\not\models\exists z\varphi_{i}(c_{1},c_{2},z)$”. We first
estimate its probability.

Take a pair of distinct elements $c_{1},c_{2}\in C_{n}$ and
$i\in\mathcal{I}_{5}$. There are three cases: (a) $c_{2}=f_{j}(c_{1})$ for
some $j$, (b) $c_{1}=f_{j}(c_{2})$ for some $j$, and (c) none of the above
holds. Let us go in details through case (a). Let
$a_{1}=\mathfrak{pat}(c_{1})$, $a_{2}=\mathfrak{f}_{j}(a_{1})$. In this case
$c_{1}$ and $c_{2}$ belong to the same copy of ${\mathfrak{B}}$, say to
${\mathfrak{B}}\times\\{k\\}$, and ${\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2}})$ was set in Stage 0 to be equal to
${\rm tp}^{{{\mathfrak{A}}}}({a_{1},a_{2}})$. If
${\mathfrak{A}}\models\varphi_{i}(a_{1},a_{2},a_{1})$ then
${\mathfrak{C}}_{n}\models\varphi_{i}(c_{1},c_{2},c_{1})$, and if
${\mathfrak{A}}\models\varphi_{i}(a_{1},a_{2},a_{2})$ then
${\mathfrak{C}}_{n}\models\varphi_{i}(c_{1},c_{2},c_{2})$, so the pair
$a_{1},a_{2}$ has a witness for the $\varphi_{i}$. Otherwise, there is an
element $a_{3}$, $a_{3}\not\in\\{a_{1},a_{2}\\}$ such that
${\mathfrak{A}}\models\varphi_{i}(a_{1},a_{2},a_{3})$. Let $\alpha_{3}={\rm
tp}^{{{\mathfrak{A}}}}({a_{3}})$. Consider now any element $c_{3}$ of type
$\alpha_{3}$ in ${\mathfrak{B}}\times(\\{1,\ldots k-1,k+1,\ldots n)$. During
the construction of ${\mathfrak{C}}_{n}$ we randomly chose ${\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{3}})$ and ${\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{2},c_{3}})$ (Stage 1) and ${\rm
pretp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2},c_{3}})$ (Stage 2a). We had at most
$|\mbox{\large\boldmath$\beta$}|$ choices for ${\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{3}})$ and at least one of them was ${\rm
tp}^{{{\mathfrak{A}}}}({a_{1},a_{3}})$. Analogously, for ${\rm
tp}^{{{\mathfrak{C}}_{n}}}({c_{2},c_{3}})$. We had at most
$|\mbox{\large\boldmath$\gamma$}^{*}|$ choices for ${\rm
pretp}^{{{\mathfrak{C}}_{n}}}({c_{1},c_{2},c_{3}})$ and at least one of them
was ${\rm pretp}^{{{\mathfrak{A}}}}({a_{1},a_{2},a_{3}})$. Hence the
probability that $\\{c_{1}\rightarrow a_{1},c_{2}\rightarrow
a_{2},c_{3}\rightarrow a_{3}\\}$ is a partial isomorphism and hence
${\mathfrak{C}}_{n}\models\varphi_{i}(c_{1},c_{2},c_{3})$ is at least
$\epsilon=1/(|\mbox{\large\boldmath$\beta$}|^{2}\cdot|\mbox{\large\boldmath$\gamma$}^{*}|)$.
As there are at least $n-1$ elements of type $\alpha_{3}$ in
${\mathfrak{B}}\times(\\{1,\ldots k-1,k+1,\ldots n)$ (at least one in each
copy of ${\mathfrak{B}}$), and the structures on $c_{1},c_{2},c_{3}$ and
$c_{1},c_{2},c_{3}^{\prime}$ for different such elements
$c_{3},c_{3}^{\prime}$ are chosen independently, it follows that ${\rm
Prob}(F^{i}_{c_{1},c_{2}})\leqslant(1-\epsilon)^{n-1}$.

Note that the event “${\mathfrak{C}}_{n}\not\models\varphi$” is equal to
$\bigcup_{c_{1}\not=c_{2}}\bigcup_{i\in\mathcal{I}_{5}}F^{i}_{c_{1},c_{2}}$.
Hence

${\rm Prob}({\mathfrak{C}}_{n}\not\models\varphi)={\rm
Prob}(\bigcup_{c_{1}\not=c_{2}}\bigcup_{i\in\mathcal{I}_{5}}F^{i}_{c_{1},c_{2}})\leqslant\sum_{c_{1}\not=c_{2}}\sum_{i\in\mathcal{I}_{5}}{\rm
Prob}(F^{i}_{c_{1},c_{2}})=n|B|(n|B|-1)|\mathcal{I}_{5}|(1-\epsilon)^{n-1}.$

The limit of the above estimation when $n$ approaches infinity is $0$ (note
that $\epsilon$, $|B|$ and $|\mathcal{I}_{5}|$ do not depend on $n$). In
particular for some $n\in{\mathbb{N}}$ we have that ${\rm
Prob}({\mathfrak{C}}_{n}\not\models\varphi)<1$ and thus one of the possible
(finite number) of choices of ${\mathfrak{C}}_{n}$ is indeed a finite model of
$\varphi$.

We can refine this reasoning and achieve an exponential bound on the size of a
minimal finite model by applying the following claim:

There exists an exponential function $f$ such that
$\text{Prob}({\mathfrak{C}}_{n}\not\models\varphi)<1$ for $n\geqslant
f(|\varphi|)$.

######  Proof 5.2.

Denote by $|\varphi|$ the size of $\varphi$ in any reasonable encoding. In
Section 4.1 we estimated the size of the domain of ${\mathfrak{B}}$ to be
exponential in $|\varphi|$. Hence, there exist a polynomial $p_{1}(x)$ such
that $|B|\leqslant 2^{p_{1}(|\varphi|)}$. Similarly, from Lemma 1 we get that
there exists polynomials $q_{1}(x)$ and $q_{2}(x)$ such that
$|\mbox{\large\boldmath$\beta$}|\leqslant 2^{q_{1}(|\varphi|)}$ and
$|\mbox{\large\boldmath$\gamma$}^{*}|\leqslant 2^{q_{2}(|\varphi|)}$. Wlog. we
can assume that $p_{1}(x),q_{1}(x),q_{2}(x)\geqslant 10x$.

Let us state the bound on the probability of the event
“${\mathfrak{C}}_{n+1}\not\models\varphi$” in terms of $p(x):=2p_{1}(x),$ and
$q(x):=2q_{1}(x)+q_{2}(x)$:

$\text{Prob}({\mathfrak{C}}_{n+1}\not\models\varphi)=(n+1)|B|((n+1)|B|-1)|\mathcal{I}_{5}|(1-\epsilon)^{n}\leqslant
10n^{2}2^{p(|\varphi|)}|\varphi|(1-\varepsilon(|\varphi|))^{n},$ (12)

where $\varepsilon(x)$ is defined as $\varepsilon(x):=1/{2^{q(x)}}$ and
therefore $\varepsilon(|\varphi|)\leqslant
1/(|\mbox{\large\boldmath$\beta$}|^{2}\cdot|\mbox{\large\boldmath$\gamma$}^{*}|)=\epsilon$.

We need to show that taking $n=f(|\varphi|)$ for some exponential function $f$
makes the right-hand side of (12) strictly less than one.444We did not tried
to optimize this bound for the sake of simplicity. Let $f(x)=10\ln(2)\cdot
2^{2q(x)}p(x)$. Then taking log on both sides of:

$10f(x)^{2}2^{p(x)}x(1-\varepsilon(x))^{f(x)}<1$

gives us:

$\log(10)+2\log(f(x))+p(x)+\log(x)+f(x)\log(1-\varepsilon(x))<0.$

Recall that the Mercator Series states that
$\ln(1-x)=-\sum_{i=1}^{\infty}{x^{i}/i}$ for all $x\in(-1,1)$. Hence we want
to see that:

$\log(10)+2\log(f(x))+p(x)+\log(x)<\frac{f(x)}{\ln(2)}\sum_{i=1}^{\infty}{\frac{\varepsilon(x)^{i}}{i}}.$

Since $f(x)$ and $\varepsilon(x)$ are always positive, the above inequality is
implied by:

$\log(10)+2\log(f(x))+p(x)+\log(x)<f(x)\varepsilon(x)/\ln(2).$

Expanding the definition of $\varepsilon(x)$ and $f(x)$:

$\log(10)+2\log(10\ln(2))+4q(x)+2\log(p(x))+p(x)+\log(x)<10\ln(2)\cdot
2^{2q(x)}p(x)\cdot 2^{-q(x)}/\ln(2).$

Finally we get:

$\log(10)+2\log(10\ln(2))+4q(x)+2\log(p(x))+p(x)+\log(x)<10\cdot
2^{q(x)}p(x),$

and it is easy to see that each term on the left-hand side is bounded by
$2^{q(x)}p(x)$ (possibly for a large enough x), so the inequality holds.

Recalling Lemma 11 we conclude.

######  Theorem 14.

$\mbox{\rm FO}^{3}_{-}$ has the finite (exponential) model property. The
finite satisfiability problem (=satisfiability problem) for $\mbox{\rm
FO}^{3}_{-}$ is NExpTime-complete.

##  6 Conclusions

We proposed an extension of the uniform-one dimensional fragment by mixed
blocks of quantifiers, $\mbox{\rm AUF}_{1}$. We identified its two
subfragments: $\mbox{\rm AUF}_{1}^{-}$ and $\mbox{\rm AUF}_{1}^{3}$ with
exponential model property and NExpTime-complete satisfiability problem. We
also extended $\mbox{\rm AUF}_{1}^{3}$ to a more expressive decidable fragment
of $\mbox{\rm FO}^{3}$. The main open question which is left is whether the
whole $\mbox{\rm AUF}_{1}$ is decidable.

Another question is what happens to our logics if a use of equality (free or
uniform) is allowed? We conclude this paper with some notes concerning this
question.

###  6.1 Notes on equality

All the results are obtained in this paper in the absence of equality. We have
two options of adding equality: first, we can treat the equality symbol as
relation symbols from the signature and allow only its uniform use, as in the
formula $\forall xy(P(x)\wedge P(y)\rightarrow x=y)$; second, we can allow for
free use of equality, as e.g. in the formula $\forall xyz(R(x,y,z)\rightarrow
x\not=y)$. We recall that $\mbox{\rm UF}_{1}$ with free use of equality is
still decidable and has the exponential model property [KK14].

Regarding the free use of equality consider the following $\mbox{\rm
AUF}_{1}^{3}$ formula:

$\exists xS(x)\wedge\forall x\exists y\forall z(\neg S(y)\wedge
R(x,y,z)\wedge(x=z\vee\neg R(z,y,x))$

It is satisfied in the model whose universe is the set of natural numbers, $S$
is true only at $0$ and $R(x,y,z)$ is true iff $y=x+1$. It is readily verified
that there are no finite models (every $x$ needs to take a fresh $y$ as a
witness since otherwise $R(z,y,x)$ would not hold for one of the earlier $z$,
$z\not=x$).

The above example can be simply adapted to the case of $\mbox{\rm
AUF}_{1}^{-}$ with free use of equality. We just add a dummy existentially
quantified variable $t$ and require it to be equal to the previous,
universally quantified variable $z$. To accommodate all the variables we
increase the arity of $R$ by $2$ (one can think that the first and the last
position of $R$ from the previous example have been doubled):

$\exists xS(x)\wedge\forall x\exists y\forall z\exists t.(t=z\wedge\neg
S(y)\wedge R(x,x,y,z,t)\wedge(x=z\vee\neg R(z,t,y,x,x))).$

These examples show that satisfiability and finite satisfiability are
different problems for $\mbox{\rm AUF}_{1}^{-}$/$\mbox{\rm AUF}_{1}^{3}$ with
free use of equality. We do not know if they are decidable, however. This
issue is left for further investigations.

As for the uniform use, we suspect that it does not spoil the finite model
property. Such uniform use allow one to say that some $1$-types are realized
at most once. We believe that treating the realizations of such types with
appropriate care (extending the approach with _kings_ and _court_ for
$\mbox{\rm FO}^{2}$ from [GKV97]) one can show that the finite model property
survives. We leave the details to be worked out.

It is clear that allowing equality in $\mbox{\rm FO}^{3}_{-}$ spoils the
decidability, as $\mbox{\rm FO}^{3}$ contains formulas of the form $\forall
x\forall y\exists z\varphi$, for quantifier-free $\varphi$, which form an
undecidable subclass of Gödel class with equality [Gol84].

## Acknowledgements

This work is supported by NCN grant No. 2021/41/B/ST6/00996. The authors thank
the reviewers of the conference papers [FK23] and [Kie23], on which this
article is based, for their helpful comments.

## References

  * [AvBN98] Hajnal Andréka, Johan van Benthem, and István Németi.  Modal languages and bounded fragments of predicate logic.  Journal of Philosophical Logic, 27:217–274, 1998.  doi:10.1023/A:1004275029985.
  * [BBC+16] Saguy Benaim, Michael Benedikt, Witold Charatonik, Emanuel Kieronski, Rastislav Lenhardt, Filip Mazowiecki, and James Worrell.  Complexity of two-variable logic on finite trees.  ACM Trans. Comput. Log., 17(4):32:1–32:38, 2016.  doi:10.1145/2996796.
  * [BDM+11] Mikołaj Bojańczyk, Claire David, Anca Muscholl, Thomas Schwentick, and Luc Segoufin.  Two-variable logic on data words.  ACM Trans. Comput. Log., 12(4):27, 2011.  doi:10.1145/1970398.1970403.
  * [BKP23] Bartosz Bednarczyk, Daumantas Kojelis, and Ian Pratt-Hartmann.  On the limits of decision: the adjacent fragment of first-order logic.  In 50th International Colloquium on Automata, Languages, and Programming, ICALP 2023, volume 261 of LIPIcs, pages 111:1–111:21, 2023\.  doi:10.4230/LIPICS.ICALP.2023.111.
  * [BM17] Simone Bova and Fabio Mogavero.  Herbrand property, finite quasi-herbrand models, and a chandra-merlin theorem for quantified conjunctive queries.  In 32nd Annual ACM/IEEE Symposium on Logic in Computer Science, LICS 2017, pages 1–12. IEEE Computer Society, 2017.  doi:10.1109/LICS.2017.8005073.
  * [BMP17] Pierre Bourhis, Michael Morak, and Andreas Pieris.  Making cross products and guarded ontology languages compatible.  In Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI 2017, pages 880–886, 2017.  doi:10.24963/IJCAI.2017/122.
  * [BMSS09] Mikołaj Bojanczyk, Anca Muscholl, Thomas Schwentick, and Luc Segoufin.  Two-variable logic on data trees and XML reasoning.  J. ACM, 56(3), 2009.  doi:10.1145/1516512.1516515.
  * [BtCS15] Vince Bárány, Balder ten Cate, and Luc Segoufin.  Guarded negation.  J. ACM, 62(3):22, 2015.  doi:10.1145/2701414.
  * [CW13] Witold Charatonik and Piotr Witkowski.  Two-variable logic with counting and trees.  In 28th Annual ACM/IEEE Symposium on Logic in Computer Science, LICS 2013, pages 73–82, 2013.  doi:10.1109/LICS.2013.12.
  * [CW15] Witold Charatonik and Piotr Witkowski.  Two-variable logic with counting and a linear order.  In 24th EACSL Annual Conference on Computer Science Logic, CSL 2015,, volume 41 of LIPIcs, pages 631–647, 2015.  doi:10.4230/LIPICS.CSL.2015.631.
  * [FK23] Oskar Fiuk and Emanuel Kieronski.  An excursion to the border of decidability: between two- and three-variable logic.  In LPAR 2023: Proceedings of 24th International Conference on Logic for Programming, Artificial Intelligence and Reasoning, volume 94 of EPiC Series in Computing, pages 205–223, 2023.  doi:10.29007/1XNS.
  * [FKM24] Oskar Fiuk, Emanuel Kieronski, and Vincent Michielini.  On the complexity of Maslov class K.  In In preparation, 2024.
  * [GKV97] Erich Grädel, Phokion Kolaitis, and Moshe Y. Vardi.  On the decision problem for two-variable first-order logic.  Bulletin of Symbolic Logic, 3(1):53–69, 1997.  doi:10.2307/421196.
  * [Gol84] Warren D. Goldfarb.  The unsolvability of the Gödel class with identity.  J. Symb. Logic, 49:1237–1252, 1984.  doi:10.2307/2274274.
  * [GOR97] Erich Grädel, Martin Otto, and Eric Rosen.  Two-variable logic with counting is decidable.  In 12th Annual IEEE Symposium on Logic in Computer Science, LICS 1997, pages 306–317. IEEE, 1997.  doi:10.1109/LICS.1997.614957.
  * [GS83] Yuri Gurevich and Saharon Shelah.  Random models and the Gödel case of the decision problem.  J. Symbolic Logic, 48(4):1120–1124, 1983.  doi:10.2307/2273674.
  * [Hal35] Philip Hall.  On representatives of subsets.  Journal of the London Mathematical Society, s1-10(1):26–30, 1935\.  doi:10.1112/jlms/s1-10.37.26.
  * [HK14] Lauri Hella and Antti Kuusisto.  One-dimensional fragment of first-order logic.  In Proceedings of Advances in Modal Logic, 2014, pages 274–293, 2014.
  * [HS99] Ullrich Hustadt and Renate Schmidt.  Maslov’s class K revisited.  In Automated Deduction — CADE-16, pages 172–186, 1999.  doi:10.1007/3-540-48660-7\\_12.
  * [Kaz06] Yevgeny Kazakov.  Saturation-based decision procedures for extensions of the guarded fragment.  PhD thesis, Universität des Saarlandes, Saarbrücken, Germany, 2006\.
  * [Kie05] Emanuel Kieroński.  Results on the guarded fragment with equivalence or transitive relations.  In 19th International Workshop, CSL 2005, 14th Annual Conference of the EACSL, volume 3634 of LNCS, pages 309–324. Springer, 2005.  doi:10.1007/11538363\\_22.
  * [Kie19] Emanuel Kieronski.  One-dimensional guarded fragments.  In 44th International Symposium on Mathematical Foundations of Computer Science, MFCS 2019, volume 138 of LIPIcs, pages 16:1–16:14, 2019.  doi:10.4230/LIPIcs.MFCS.2019.16.
  * [Kie23] Emanuel Kieronski.  A uniform one-dimensional fragment with alternation of quantifiers.  In Fourteenth International Symposium on Games, Automata, Logics, and Formal Verification, GandALF 2023, volume 390 of EPTCS, pages 1–15, 2023.  doi:10.4204/EPTCS.390.1.
  * [KK14] Emanuel Kieronski and Antti Kuusisto.  Complexity and expressivity of uniform one-dimensional fragment with equality.  In Mathematical Foundations of Computer Science 2014 - 39th International Symposium, MFCS 2014, Proceedings, Part I, volume 8634 of Lecture Notes in Computer Science, pages 365–376, 2014.  doi:10.1007/978-3-662-44522-8\\_31.
  * [KK15] Emanuel Kieronski and Antti Kuusisto.  Uniform one-dimensional fragments with one equivalence relation.  In 24th EACSL Annual Conference on Computer Science Logic, CSL 2015, volume 41 of LIPIcs, pages 597–615, 2015.  doi:10.4230/LIPICS.CSL.2015.597.
  * [KLPS20] Andreas Krebs, Kamal Lodaya, Paritosh K. Pandya, and Howard Straubing.  Two-variable logics with some betweenness relations: Expressiveness, satisfiability and membership.  Log. Methods Comput. Sci., 16(3), 2020.  doi:10.23638/LMCS-16(3:16)2020.
  * [KMPT14] Emanuel Kieronski, Jakub Michaliszyn, Ian Pratt-Hartmann, and Lidia Tendera.  Two-variable first-order logic with equivalence closure.  SIAM J. Comput., 43(3):1012–1063, 2014.  doi:10.1137/120900095.
  * [KMW62] A. S. Kahr, Edward F. Moore, and Hao Wang.  Entscheidungsproblem reduced to the $\forall\exists\forall$ case.  Proc. Nat. Acad. Sci. U.S.A., 48:365–377, 1962.
  * [KO12] Emanuel Kieroński and Martin Otto.  Small substructures and decidability issues for first-order logic with two variables.  Journal of Symbolic Logic, 77:729–765, 2012.  doi:10.2178/JSL/1344862160.
  * [KR21] Emanuel Kieronski and Sebastian Rudolph.  Finite model theory of the triguarded fragment and related logics.  In 36th Annual ACM/IEEE Symposium on Logic in Computer Science, LICS 2021, pages 1–13, 2021.  doi:10.1109/LICS52264.2021.9470734.
  * [Kuu16] Antti Kuusisto.  On the uniform one-dimensional fragment.  In 29th International Workshop on Description Logics, Cape Town, South Africa, April 22-25, 2016., 2016.
  * [Lew80] Harry R. Lewis.  Complexity results for classes of quantificational formulas.  Journal of Computer and System Sciences, 21(3):317 – 353, 1980\.  doi:10.1016/0022-0000(80)90027-6.
  * [Mas71] Sergei J. Maslov.  The inverse method for establishing deducibility for logical calculi.  The Calculi of Symbolic Logic I: Proceedings of the Steklov Institute of Mathematics, 98, 1971.
  * [Mor75] Michael Mortimer.  On languages with two variables.  Zeitschrift für Mathematische Logik und Grundlagen der Mathematik, 21:135–140, 1975.
  * [MP15] Fabio Mogavero and Giuseppe Perelli.  Binding forms in first-order logic.  In 24th EACSL Annual Conference on Computer Science Logic, CSL 2015, volume 41 of LIPIcs, pages 648–665, 2015.  doi:10.4230/LIPICS.CSL.2015.648.
  * [Pra10] Ian Pratt-Hartmann.  The two-variable fragment with counting revisited.  In Logic, Language, Information and Computation, 17th International Workshop, WoLLIC 2010, pages 42–54, 2010.  doi:10.1007/978-3-642-13824-9\\_4.
  * [Pra15] Ian Pratt-Hartmann.  The two-variable fragment with counting and equivalence.  Math. Log. Q., 61(6):474–515, 2015.  doi:10.1002/malq.201400102.
  * [Pra21] Ian Pratt-Hartmann.  Fluted logic with counting.  In 48th International Colloquium on Automata, Languages, and Programming, ICALP 2021, volume 198 of LIPIcs, pages 141:1–141:17, 2021\.  doi:10.4230/LIPIcs.ICALP.2021.141.
  * [PST97] Leszek Pacholski, Wiesław Szwast, and Lidia Tendera.  Complexity of two-variable logic with counting.  In 12th Annual IEEE Symposium on Logic in Computer Science, LICS 1997, pages 318–327. IEEE, 1997.  doi:10.1109/LICS.1997.614958.
  * [PST19] Ian Pratt-Hartmann, Wieslaw Szwast, and Lidia Tendera.  The fluted fragment revisited.  J. Symb. Log., 84(3):1020–1048, 2019.  doi:10.1017/JSL.2019.33.
  * [Qui69] Willard V. Quine.  On the limits of decision.  In Proceedings of the 14th International Congress of Philosophy, volume III, pages 57–62, 1969.
  * [RŠ18] Sebastian Rudolph and Mantas Šimkus.  The triguarded fragment of first-order logic.  In 22nd International Conference on Logic for Programming, Artificial Intelligence and Reasoning, LPAR-22, volume 57 of EPiC Series in Computing, pages 604–619, 2018.  doi:10.29007/M8TS.
  * [Sco62] Dana Scott.  A decision method for validity of sentences in two variables.  Journal Symbolic Logic, 27:477, 1962.
  * [tCS13] Balder ten Cate and Luc Segoufin.  Unary negation.  Logical Methods in Comp. Sc., 9(3), 2013.  doi:10.2168/LMCS-9(3:25)2013.