††thanks: Corresponding author

# Planetary seismology as a test of modified gravity proposals

Aleksander Kozak<EMAIL_ADDRESS>Institute of Theoretical
physics, University of Wroclaw, pl. Maxa Borna 9, 50-206 Wroclaw, Poland
Aneta Wojnar<EMAIL_ADDRESS>Departamento de Física Teórica & IPARCOS,
Universidad Complutense de Madrid, E-28040, Madrid, Spain

###### Abstract

We demonstrate that it is possible to test models of gravity, such as Palatini
$f(R)$ and Eddington-inspired Born-Infeld models, using seismic data from
Earth. By incorporating additional limitations on Earth’s moment of inertia
and mass given from observational data, the models’ parameters can be
restricted to a $2\sigma$ level of accuracy. Our novel tool provides that
$\beta\lesssim 10^{9}\text{m}^{2}$ for Palatini and $\epsilon\lesssim 4\cdot
10^{9}\text{m}^{2}$ for Eddington-inspired Born-Infeld gravity. We also
discuss further enhancements to the proposed method, aimed at imposing even
more stringent constraints on modified gravity proposals.

##  I Introduction

Numerous proposals have been put forward to extend General Relativity (GR) in
order to address the mysteries of the dark sector of the Universe huterer ;
Copeland:2006wr ; Nojiri:2006ri ; nojiri2 ; nojiri3 ; Capozziello:2007ec ;
Carroll:2004de ; cantata , as well as the discrepancy between observations of
the visible components of galaxies and their dynamical mass rubin , and the
existence of apparently ”too massive” compact objects lina ; as ; craw ; NSBH
; abotHBH ; sak3 . These modifications to GR appear inevitable in the
cosmological regime although they must be significantly suppressed in small-
scale systems such as compact objects and the Solar System. Moreover, the
recent analyses of gravitational parameter space diagram baker  indicates
untested regions where curvature values correspond to galaxies and
(sub-)stellar objects. It is speculated that those untested domains, which
separate small scale systems from the cosmological one, could shed light on
corrections to GR.

Understanding a need of various but complementary tests for the gravitational
proposals and to fill the mentioned gap, we proposed a new one, based on the
planetary seismology olek1 . So far, only seismic data from low-mass stars
(asteroseismology) bell ; casa  and from the Sun (helioseismology) saltas1 ;
saltas2  has been used to constrain fundamental theories. Growing number and
accuracy of observational data on astrophysical bodies has already helped to
constrain or even to rule out some of the existing theories of gravity; for
example, Multi-messenger Astronomy ab1 ; ab2 ; ez ; pat ; ab3  excluded those
models which predict that the speed of gravitational waves is different than
the speed of light baker2 ; ez2 ; crem ; sakjain ; cope ; boran ; gong .
Another example are soft equations of state which were also ruled out in the
framework of GR111However, in the framework of modified gravity proposals
those equations of state are still a viable description of the matter
properties inside neutron stars, see review  and references therein. because
of inability to provide high neutron stars’ masses mar ; shi ; ruiz ; rez ;
fon ; cro ; ab4 ; cou .

However, in comparison to compact objects and even stellar ones in which an
equation of state as well as atmospheric properties play a crucial role,
carrying at the same time large uncertainties related to their description
olek4 ; debora ; wojnar , Earth seismology provides data which conveys
information on the planet’s interior poirier ; prem ; kus ; ken ; ken2 ; iris
. As we will demonstrate, the seismic data, together with high accuracy of
Earth’s mass and moment of inertia measurements, are novel and remarkable tool
to constrain models of gravity, at the same time using a well-understood
physics such that one can avoid some of uncertainties related to model’s
assumptions.

Moreover, our knowledge on our planet’s interior has been increased
significantly in last years, not only thanks to the improvements of
seismographic tools prem ; mush ; frost ; step ; pham , but also due to
reaching Earth’s interior’s temperatures and pressures in laboratories. This
allowed to study iron’s properties222Iron and its compounds are main elements
in Earth’s inner and outer core. and behaviour in these regimes laser . On the
other hand, new neutrino telescopes will also provide the information on
density, composition, and abundances of light elements in the outer core in
the deepest part of our planet topography ; top ; top2 ; top3 , decreasing
even more uncertainties related to Earth’s core characteristics.

One can also be worried about the insignificance of the modified gravity
effects in planetary physics. Indeed, the impact on the densities and
thicknesses of the layers is small but still notable olek2 ; olek3 ; wojnar2 ;
wojnar3 , and in light of the mentioned facts that we possess much more
information with better accuracy about the Solar System planets kaula ; bill ;
folk ; konopliv ; smith ; folk2 , and in particular, about Earth333And also on
Mars when seismic data obtained from analyzing waves created by marsquakes,
thumps of meteorite impacts, surface vibrations caused by activity in Mars’
atmosphere, and by weather phenomena, e.g. dust storms nasa , are ready.
ziemia , we are able to use the available data to constrain theories. As we
will see, with our simplified approach, we can do it up to a $2\sigma$ level
of accuracy.

The paper is organized as follows. In the section II we introduce the basic
notions relating to Earth’s modelling. We also recall the main assumptions and
tools used in the Preliminary Reference Earth Model (PREM), widely used in
geology and (exo-)planet science. Ricci-based gravity and Earth model
resulting from this proposal of gravity are presented in the section III. Our
methodology and results are discussed in the section IV. In the section V we
draw our conclusions.

##  II PREM - seismology and gravity

In what follows, let us discuss the most used global seismological Earth
model, which is a starting point for the more accurate Earth’s kus ; ken ;
ken2 ; iris  and other terrestrial wep  and even gaseous seager
(exo-)planets’ models. Preliminary Reference Earth Model prem  is based on
velocity-depth profiles given by the travel-time distance curves for seismic
waves and on periods of free oscillations poirier ; bolt ; bullen . It
provides pressure, density and elastic moduli profiles as functions of depth,
and as it will be clear in the further part, the lacking element to determine
them is a hydrostatic equilibrium equation.

The PREM is a 1-dimensional model and it adopts the following assumptions:

  * •

There is no exchange of heat between different layers (adiabatic compression)
therefore there is no additional term in the hydrostatic equilibrium equation
related to the temperature variation.

  * •

The planet is a spherical-symmetric ball in hydrostatic equilibrium given by

$\frac{dP}{dr}=-4\pi G\rho r^{-2}\int^{r}_{0}\rho r^{2}dr=:-\rho g$ (1)

where the pressure $P$, density $\rho$ and acceleration of gravity $g$ are
functions of the radial coordinate $r$ (or depth).

  * •

The planet consists of radially symmetric shells with the given444In the
original PREM model, the central density is not the theory’s free parameter,
but can be obtained by integrating the relevant equations. Instead, one
assumes the value of the density below the crust, and all densities in the
outer layers are described by Birch’s law. In our approach, we use the PREM
value, as we do not indent to change the model in the outermost layer due to a
weak effect of modified gravity.  density jump between the inner and outer
core $\Delta\rho=600$, central density $\rho_{c}=13050$ and density at the
mantle’s base $\rho_{m}=5563$ (in kg/m3).

Moreover, the mass equation

$M=4\pi\int_{0}^{R}r^{2}\rho(r)dr$ (2)

and moment of inertia

$I=\frac{8}{3}\pi\int_{0}^{R}r^{4}\rho(r)dr,$ (3)

where $R$ is Earth’s radius, play a role of the constraints: they are given by
observations with a high accuracy luzum ; chen .

On the other hand, the outer layers’ density profile (for upper mantle) is
described by the empirical Birch law

$\rho=a+bv_{p},$ (4)

where $a$ and $b$ are parameters which depend on the material in the upper
mantle prem  (they depend on the mean atomic mass) while $v_{p}$ is the
longitudinal elastic wave. Together with the transverse elastic wave $v_{s}$,
we can define the seismic parameter $\Phi_{s}$ as555See the derivation in, for
example, poirier .

$\Phi_{s}=v_{p}^{2}-\frac{4}{3}v_{s}^{2}.$ (5)

It is related to the elastic properties of an isotropic material; more
specifically, it is related to the bulk modulus $K$ (also called
incompressibility)

$\Phi_{s}=\frac{K}{\rho}.$ (6)

By the definition of the bulk modulus

$K=\frac{dP}{d\mathrm{ln}\rho}$ (7)

we can also express the seismic parameter in terms of the material’s
properties

$\Phi_{s}=\frac{dP}{d\rho},$ (8)

and hence it carries the information on the equation of state. Using this
feature in (1)

$\frac{d\rho}{dr}=-\rho g\Phi_{s}^{-1},$ (9)

together with the mass equation (2) and seismic data, that is, the
longitudinal and transverse elastic waves $v_{p}$ and $v_{s}$ prem , one can
obtain the density profile, as presented in the figure 1.

Figure 1: [color online] The density profile given by the Preliminary
Reference Earth Model prem  in which Newtonian gravity is assumed. The
velocities’ plots are obtained from data without using any theory of gravity.
The primary waves are the same as the longitudinal waves, while the secondary
waves are transverse in nature. The units are in km/s for velocities, while
the densities are in kg/m3.

It should be clear now how the planets’ models depend on theory of gravity;
modifying the hydrostatic equilibrium (1) or mass equations (2), the density
profile will vary with respect to the one provided by the PREM. Let us also
notice that more realistic seismic models of Earth also used Newtonian
hydrostatic equilibrium equations iris , therefore they are also gravity model
dependent. However, equipped with gravity-independent Earth model topography ;
top ; top2 ; top3 , as discussed in olek1 , can be a powerful tool in the
nearest future to constrain gravitational proposals.

##  III Earth model in Ricci-based gravity

In what follows, let us discuss a class of metric-affine proposals of gravity,
the so-called Ricci-based gravity (RBG) theories (see, e.g. alfonso ). This
particular class’ action can be written as

$\mathcal{S}=\int
d^{4}x\sqrt{-g}\mathcal{L}_{G}(g_{\mu\nu},R_{\mu\nu})+\mathcal{S}_{m}(g_{\mu\nu},\psi_{m})\
,$ (10)

where $g$ is the determinant of the space-time metric $g_{\mu\nu}$ while
$R_{\mu\nu}$ is the symmetric Ricci tensor, independent of the metric as it
has been constructed with the affine connection
$\Gamma\equiv\Gamma_{\mu\nu}^{\lambda}$ only. Let us now define an object
${M^{\mu}}_{\nu}\equiv g^{\mu\alpha}R_{\alpha\nu}$ which we will use to
construct the gravitational Lagrangian $\mathcal{L}_{G}$ such that we will
deal with a scalar function being built with powers of traces of
${M^{\mu}}_{\nu}$.

On the other hand, the matter action

$\mathcal{S}_{m}=\int d^{4}x\sqrt{-g}\mathcal{L}_{m}(g_{\mu\nu},\psi_{m})$
(11)

in this proposal is minimally coupled to the metric, where $\psi_{m}$ are the
matter fields. That is, the antisymmetric part of the connection (torsion) is
neglected as in the case of minimally coupled bosonic fields666Moreover, even
fermionic particles, for example degenerate matter, can be effectively
described by a fluid approach, given by, for example, the perfect fluid
energy-momentum tensor. the problem trivializes to the symmetric part only
alfonso . Similarly, we focus on the symmetric part of the Ricci tensor to
avoid potential ghost-like instabilities BeltranJimenez:2019acz ;
Jimenez:2020dpn . Such an approach allows us to consider many theories of
gravity such as GR, Palatini $f(R)$ gravity, Eddington-inspired Born-Infeld
(EiBI) gravity vollick  and many of its extensions BeltranJimenez:2017doy .

Theories of gravity which can be included in the above gravitational action
have a nice feature: although the field equations possess rather complicated
forms, they can be rewritten in a more convenient way (see e.g
BeltranJimenez:2017doy )

${G^{\mu}}_{\nu}(q)=\frac{\kappa}{|\hat{\Omega}|^{1/2}}\left({T^{\mu}}_{\nu}-\delta^{\mu}_{\nu}\left(\mathcal{L}_{G}+\frac{T}{2}\right)\right)\
,$ (12)

where $|\hat{\Omega}|$ is a determinant of the deformation matrix defined
below, $T=g^{\mu\nu}T_{\mu\nu}$ is the trace of the energy-momentum tensor of
the matter fields
$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\partial\mathcal{L}_{m}}{\partial
g^{\mu\nu}}$. The Einstein tensor ${G^{\mu}}_{\nu}(q)$ is associated to a
tensor $q_{\mu\nu}$ such that the connection $\Gamma$ results to be the Levi-
Civita connection of it:

$\nabla_{\mu}^{\Gamma}(\sqrt{-q}q^{\alpha\beta})=0.$ (13)

This tensor, for the given formalism, is related to the space-time metric
$g_{\mu\nu}$ via the relation

$q_{\mu\nu}=g_{\mu\alpha}{\Omega^{\alpha}}_{\nu},$ (14)

where ${\Omega^{\alpha}}_{\nu}$ is the deformation matrix which depends on the
particular theory, that is, gravitational Lagrangian $\mathcal{L}_{G}$
considered. It turns out that $\mathcal{L}_{G}$ can be expressed on-shell as a
function of the matter fields and the space-time metric $g_{\mu\nu}$, so the
deformation matrix ${\Omega^{\alpha}}_{\nu}$ is also a function of them. Its
determinant is denoted by vertical bars as in equation (12). Moreover, it is
obvious now from the form of (12) that RBG theories provide the second-order
field equations which in vacuum, that is, ${T_{\mu}}^{\nu}=0$, reduce to the
GR counterparts. This also means that there is no extra degree of freedom
propagating in these theories apart from the usual two polarization of the
gravitational field.

We are interested in the non-relativistic limit of the field equations (12).
In the case of Palatini $f(R)$ toniato  and EiBI banados ; pani  gravities the
Poisson equation can be written in the form

$\nabla^{2}\phi=\frac{\kappa}{2}\Big{(}\rho+\alpha\nabla^{2}\rho\Big{)}$ (15)

where $\phi$ is gravitational potential, $\kappa=8\pi G$ while $\alpha$ a
theory parameter, taking the forms $\alpha=2\beta$ for Palatini $f(R)$
($\beta$ is the parameter accompanying the quadratic term) and
$\alpha=\epsilon/2$ for EiBI ($\epsilon=1/M_{BI}$, with $M_{BI}$ being the
Born-Infeld mass)777The similarity of the Poisson equation in those two
gravity proposals is not a coincidence - the EiBI gravity in the 1st order
approximation reduces to Palatini gravity with the quadratic term pani2 . On
the other hand, only the quadratic term enters the non-relativistic equations
since further curvature scalar terms enter the equations on the sixth order
toniato ..

The non-relativistic hydrostatic equilibrium and mass equations are

$\displaystyle\frac{d\phi}{dr}$ $\displaystyle=-\rho^{-1}\frac{dP}{dr},$ (16)
$\displaystyle M$ $\displaystyle=\int
4\pi^{\prime}\tilde{r}^{2}\rho(\tilde{r})d\tilde{r},$ (17)

or, specifically,

$\displaystyle\frac{dP}{dr}$
$\displaystyle=-\rho\left(\frac{GM(r)}{r^{2}}+\beta\kappa\rho^{\prime}(r)\right)\,\,\text{for
Palatini }f(R),$ (18) $\displaystyle\frac{dP}{dr}$
$\displaystyle=-\rho\left(\frac{\kappa M(r)}{8\pi
r^{2}}+\frac{\kappa\epsilon}{4}\frac{d\rho}{dr}\right)\,\,\,\,\,\,\,\text{for
EiBI}$ (19)

which clearly differ from (1) because of the correction term in the Poisson
equation (15). We also notice that in the non-relativistic limit, the relation
between the parameters is $\epsilon=4\beta$ thus we will denote now the
modified gravity correction by the parameter $\beta$ in the further part of
the paper. Therefore, because of the further numerical analysis, we will
rewrite those equations as

$\frac{d\rho}{dr}=-\rho g_{\text{eff}}\Phi_{s}^{-1}$ (20)

where, analogously to (1) and (9), we have introduced

$g_{\text{eff}}=\frac{GM(r)}{r^{2}}+\beta\kappa\rho^{\prime}(r).$

There already exist bounds on both proposals. It was showed that in the case
of Palatini gravity, the value of $\beta$ is related to the curvature regime
olmo2005gravity . It is so because of the fact that the Palatini curvature
scalar is proportional to the density resulting that its value also does. On
the other hand, the analytical examination of the weak-field limit yields that
$|\beta|\lesssim 2\times 10^{8}\rm\,m^{2}$ olmo2005gravity  while further
studies demonstrated that the Solar System experiments cannot deliver bounds
on the parameters because of the microphysics uncertainties toniato .
Considering EiBI gravity, the parameter $\epsilon$ lies within $-6.1\cdot
10^{15}\leq\epsilon\leq 1.1\cdot 10^{16}$m2 at $5\sigma$ confidence level
according to the newest bound pritam . Other bound of those or higher orders
were obtained earlier in jana ; casa2 ; ave ; ban . Similarly to General
Relativity, none of the considered models is able to explain the galaxy
rotation curves alejandro , so there is no bound obtained with the galaxies
catalogues yet.

Since we do not have yet a gravity-independent model of Earth, in the further
part then we will assume that PREM is an accurate model of our planet such
that we will use it, together with the mass and moment of inertia constraints,
to constrain the Ricci-based gravity. As we will see, we will prove that the
seismic data of Earth can be indeed used to restrict models of gravity to a
$2\sigma$ level of accuracy at least.

##  IV Numerical appraoch and results

The data set that served as a basis for our calculations was obtained from
prem  and references therein; it features measured values of longitudinal and
transverse seismic waves velocities for corresponding depths. Based on this
data and having made assumptions about the values of the free parameters of
the PREM: the density jump between the inner and outer core, as well as the
density at the base of the lower mantle (we neglect the upper mantle and the
crust since we assume Birch’s law there), one can calculate the density
profiles and integrate the result to get the total mass and the polar moment
of inertia. These values were treated as a measure of the accuracy of the
model, and it was these numbers that we compared our results to. The PREM
predicts the value of Earth’s mass with a very good accuracy; the measured
value is $M_{\oplus}=(5.9722\pm 0.0006)\times 10^{24}\text{kg}$ luzum ,
whereas PREM gives $M_{\text{PREM}}=5.9721\times 10^{24}\text{kg}$. The
situation is a bit worse when one considers the polar moment of inertia: the
value determined experimentally is $I_{\oplus}=(8.01736\pm 0.00097)\times
10^{37}\text{kg m}^{2}$ chen , while PREM yields
$I_{\text{PREM}}=8.01897\times 10^{37}\text{kg m}^{2}$, which is $\sim
1.5\sigma$ away from the expected value. This clearly results from the fact
that PREM is a one-dimensional model, and it does not incorporate an important
aspect that Earth is not a perfect sphere. This, obviously, poses some
problems when trying to assess the ’goodness’ of alternative gravity models.
Therefore, in our calculations, we decided to try to obtain the value
determined using PREM, while keeping the uncertainties coming from
measurements. We assume that, given a more accurate, 3-dimensional Earth
model, one would be able to increase the accuracy of the predictions (i.e. get
closer to the mean value obtained empirically), but the errors coming from
experiments would remain the same. For this reason, we aim at the following
values with accompanying errors:

$\begin{split}&M_{\text{target}}=(5.9721\pm 0.0006)\times 10^{24}\text{kg},\\\
&I_{\text{target}}=(8.01897\pm 0.00097)\times 10^{37}\text{kg
m}^{2}.\end{split}$ (21)

We adopted a brute-force approach to the problem of determining the values of
the theory’s free parameters: intervals not only of possible values of both
the density jump and the density at the base of the mantle were assumed (being
reasonably close to the experimentally measured values), but also of the
central density. This value is determined by the PREM given the total Earth’s
mass and its polar moment of inertia. We decided to make it a free parameter
and to allow it to change within a certain interval. Overall, the calculations
were performed using a Python script with $(\rho_{c},\rho_{m},\Delta\rho)$
(central density, density at the base of the mantle, density jump between
inner and outer core) taken from the following sets:

$\begin{split}&\rho_{c}\in[13050,13150]\>\text{kg m}^{-3},\\\
&\rho_{m}\in[5500,5600]\>\text{kg m}^{-3},\\\
&\Delta\rho\in[590,740]\>\text{kg m}^{-3}.\end{split}$ (22)

The variability of these three parameters was introduced in order to address
the uncertainty in determining exact values of the densities inside Earth
krasnoshchekov . We performed the calculations for different values of
$\beta$, ranging from $\beta=10^{8}\text{m}^{2}$ up to $\beta=9\times
10^{11}\text{m}^{2}$. It must be stressed once again that the model we are
using is very simple, and therefore the main aim of our computations is to
check which parameters would be crucial to a more sophisticated analysis, and
also to assess the order of magnitude of $\beta$, at which the effects of MG
are still in agreement with the observationally determined constraints. The
results of these calculations are shown in Fig 2 and 3.

Some comments concerning the integration technique are in order. First of all,
the modified hydrostatic equilibrium equation features a derivative of the
inverse of the seismic parameter $\Phi_{S}$, depending on the seismic
velocities. For simplicity, we fitted a curve to the data points and obtained
up to the third-order polynomials, whose coefficients are in agreement with
analytical approximations done by Dziewoński et al. in prem . Secondly, the
integration itself was carried out employing the Euler method with
$(\rho_{c},\rho_{m})$ treated as initial conditions at corresponding
boundaries, and $\Delta\rho$ used to compute the density at the base of the
outer core.

We also computed the errors for a more clear situation of fixed values
$(\rho_{c},\rho_{m},\Delta\rho)$, being in good agreement with the assumed
target mass and moment of inertia (i.e. within $1\sigma$ from the mean
values), for different values of the parameter $\beta$ in order to see what
the order of magnitude, at which the deviations exceed $1\sigma$ and
$2\sigma$, turns out to be. The results of this simple assessment are shown in
Fig 4.

A quick analysis of Fig 4 suggests that the effects of modified gravity become
noticeable and start exceeding the $2\sigma$ accuracy level when $\beta$ is of
order $\sim 10^{9}\text{m}^{2}$, which is similar to or lower than the upper
limits placed on the value of these parameters in other works olmo2005gravity
; ban . As one can see, the mass of the planet seems to be slightly more
sensitive to the change of the parameter $\beta$ than the polar moment of
inertia. This is due to the fact that the outermost layers of the planet,
which are unmodified, give a more significant contribution to the overall
moment of inertia than they do in the case of the mass since to calculate the
moment of inertia one raises the radius to the fourth power, while the
computation of mass requires second power. However, note that even a small
change in the core’s structure and composition of the terrestrial planets can
have a significant effect on the moment of inertia p1 ; p2 ; p3 ; p4 ; p5 ; p6
; p7 ; p8 ; p9 .

Figures 2 and 3 reveal that the density jump between the inner and outer core,
combined with an increase in the central density, has a small impact on the
overall mass and moment of inertia - at least within the $1\sigma$ region,
while the value of the density at the base of the mantle has a dramatic effect
on those quantities. This is caused by the intuitive fact that the greater
value of the jump is compensated by an increase in the central density. One
can clearly observe that the bright-blue region is shifted towards bigger
values of the density jump, meaning that a more dense core corresponds to a
bigger density difference between the inner and outer core. However, the
mantle base density determines the density profile for the outer layers and
its increase cannot be compensated by a decrease in any other quantity’s value
(for a fixed value of $\beta$). Since the mass and the moment of inertia are
determined with very good accuracy, even a small change in the boundary
conditions (here, the density at the base of the mantle) has a measurable
impact on the empirically determined quantities. In Fig 3 the effect of
modified gravity becomes visible. Interestingly, the region becomes shifted
toward greater values of the inner-outer core density jump, and toward lower
values of the density at the base of the mantle. The admissible region of the
mantle density is in all cases rather small, approximately equal to $4\text{kg
m}^{-3}$, which most likely not only exceeds the accuracy of today’s
experimental techniques but also might be irrelevant when a non-isotropic
model would be taken into account (structural differences between various
parts of the base of the mantle, such as local inhomogeneities, shifts in the
boundary between layers, and so on, might be greater than the uncertainty
related to the choice of the theory of gravity).

(a)

(b)

Figure 2: [color online] $1\sigma$ confidence regions of the theory parameters
$(\rho_{c},\rho_{m},\Delta\rho)$ for different values of the $\beta$
parameter, being of order $10^{8}$ m2. The darker color corresponds to lower
values of the central density, while the brighter one - to higher. The range
of the central density is shown in the color bar below the figures. The units
are kg/m3.

(a)

(b)

Figure 3: [color online] $1\sigma$ confidence regions of the theory parameters
$(\rho_{c},\rho_{m},\Delta\rho)$ for different values of the $\beta$
parameter, being of order $10^{9}$ m2. The darker color corresponds to lower
values of the central density, while the brighter one - to higher. The range
of the central density is shown in the color bar below the figures. The units
are kg/m3. (a)  Figure 4: [color online] Relative absolute error for the mass
and the moment of inertia of Earth. Red dots represent errors for the moment
of inertia, while blue ones correspond to the mass. The dark green stripe
represents a $1\sigma$ region for both quantities, while the light green
denotes a $2\sigma$ region. The values of $(\rho_{m},\rho_{c},\Delta\rho)$
chosen for numerical calculations are
$(5563,13050,600)\text{kg}/\text{m}^{3}$, respectively.

##  V Summary and conclusions

The present study was designed to determine if planetary seismology can be
used to constrain models of gravity. Because of our previous works, in which
we have proposed the procedure olek1  performed here and in which we have
studied in detail the structure and density profile dependencies on modified
gravity olek2 ; olek3 , we were convinced that constraining gravitational
models would be possible. However, it was not sure if one can do that with a
simplified Earth modeling as demonstrated in this work because of the
observational uncertainties. As proved in the presented study, one can put an
upper bound on the value of the $\beta$ parameter, for which the deviations of
Earth’s mass and polar moment of inertia do not exceed $2\sigma$. This value
is, for fixed density parameters inside the planet, of order $\sim
10^{9}\text{m}^{2}$, that is, $\beta\lesssim 10^{9}\text{m}^{2}$ for Palatini
$f(R)$ gravity, and $\epsilon\lesssim 4\cdot 10^{9}\text{m}^{2}$ for EiBI
one888On the other hand, the strong regime studies indicate the order $\sim
10^{8}$ m2 eva  for Palatini gravity while for the EiBI one the upper limit of
the parameter is $\lesssim 6\times 10^{6}\text{m}^{2}$ avelino ..

The seismic data prem  which we have used in our analysis consists of the
longitudinal and transverse elastic waves’ velocities and the depth of their
propagation. Both elastic waves carry information about the matter properties
inside Earth which can be encoded with the help of the seismic parameter,
appearing in the (modified) hydrostatic equilibrium and Poisson equations.
Solving them with respect to various values of a theory’s parameter and
boundary/initial conditions (such as central and jump densities between
layers), we obtain density profiles, differing with respect to the Preliminary
Reference Earth Model prem  which is based on the Newtonian gravity. On the
other hand, the observational values of Earth’s mass and moment of inertia
strongly constrain density profiles, and, as a result, the theory parameter.
Even though the model was rather simplistic, it was still possible to
constrain the Palatini parameter $\beta$ (EiBI parameter $\epsilon$) using the
data collected here on Earth.

However, our approach carries with it various limitations, caused by
assumptions and simplifications. The first and most crucial one, already
discussed previously in the text, is the spherical symmetry. Earth is not a
perfect sphere; moreover, the moment of inertia is sensitive to rotation and
particular symmetry induced by it. It is already a problem with the PREM,
whose density profile does not produce a moment of inertia which is in good
agreement with the observational value within its accuracy. Apart from it, the
PREM and our models are one-dimensional models, in addition assuming spherical
layers. Taking into account their imperfections and varying density jumps will
also have a non-zero effect on the moment of inertia and mass. The adiabatic
compression, that is, assuming that the temperature does not vary with depth
is yet another improvement which we need to take into account in our future
work.

One may also be concerned about using Birch’s law for the outer layers in
modified gravity. As already mentioned, it is an empirical law whose
coefficients were obtained experimentally. Although gravity has something to
say about matter properties999Such as, for example, chemical potential kulikov
and chemical reaction rates lecca , specific heat quantum , Debye temperature
and crystallization processes kalita , or equation of state kim ; wojnar4 .,
in this case we can safely use it. However, this law should be determined
again if we deal with seismic data from Mars - note that the coefficients
depend on the kind of material the outer layers are made of, which differ for
each terrestrial planet.

In spite of its current limitations, the study certainly offers a great tool
for constraining theories of gravity. Equipped with a quite simple Earth
model, we have been able to constrain the most popular Ricci-based gravities
up to $2\sigma$. Improving the points mentioned above, the method will provide
even stronger constraints while future neutrino telescopes will allow to use
the planet’s model which is gravity-independent, reducing even more
uncertainties related to the physical properties of the inner core. Work is
currently underway in this direction.

## Acknowledgements

AW acknowledges financial support from MICINN (Spain) Ayuda Juan de la Cierva
- incorporacíon 2020 No. IJC2020-044751-I.
The authors would like to thank prof. Dariusz Prorok for clarification
concerning the statistical data analysis.

## References

  * (1) D. Huterer and M. S. Turner, Phys. Rev. D 60, 081301 (1999).
  * (2) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006).
  * (3) S. Nojiri and S. D. Odintsov, Int. J. Geom. Meth. Mod. Phys. 4, 115 (2007).
  * (4) S. Nojiri, S.D. Odintsov, V.K. Oikonomou, Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution, Physics Reports 692 (2017).
  * (5) S. Nojiri, S.D. Odintsov, Unified cosmic history in modified gravity: from $F(R)$ theory to Lorentz non-invariant models, Physics Reports 505 (2011).
  * (6) S. Capozziello and M. Francaviglia, Gen. Rel. Grav. 40, 357 (2008).
  * (7) S. M. Carroll, A. De Felice, V. Duvvuri, D. A. Easson, M. Trodden, and M. S. Turner, Phys. Rev. D 71, 063513 (2005).
  * (8) Saridakis, et al., Modified gravity and cosmology: an update by the CANTATA network. Springer 2021.
  * (9) V. C. Rubin, N. Thonnard, W. K. Ford, Jr., Astrophys. J. 238, 471 (1980).
  * (10) M. Linares, T. Shahbaz, and J. Casares, The Astrophysical Journal 859, 54 (2018).
  * (11) J. Antoniadis et al., Science 340, 6131 (2012).
  * (12) F. Crawford, M. S. E. Roberts, J. W. T. Hessels, S. M. Ransom, M. Livingstone, C. R. Tam and V. M. Kaspi, Astrophys. J. 652, 1499 (2006).
  * (13) R. Abbott et al, Astroph. J. 896, L44 (2020).
  * (14) R. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 125, 101102 (2020).
  * (15) J. Sakstein, et al., Phys. Rev. Lett. 125, 261105 (2020).
  * (16) R. Saito, D. Yamauchi, S. Mizuno, J. Gleyzes and D. Langlois, JCAP 06, 008 (2015).
  * (17) T. Baker, D. Psaltis, and C. Skordis, Astroph. J. 802, 63 (2015).
  * (18) A. Kozak, A. Wojnar, Phys. Rev. D 104, 084097 (2021).
  * (19) E. P. Bellinger,J. Christensen-Dalsgaard, Astroph. J. Lett., 887, L1 (2019).
  * (20) J. Casanellas, In EPJ Web of Conferences (Vol. 101, p. 01015). EDP Sciences (2015).
  * (21) I. D. Saltas, I. Lopes, Phys. Rev. Lett. 123, 091103 (2019).
  * (22) I.D. Saltas, J. Christensen-Dalsgaard, Astronomy & Astrophysics 667, A115 (2022).
  * (23) B. P. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Ad- hikari, V. B. Adya, et al., Phys. Rev. Lett. 119, 161101 (2017).
  * (24) B. P. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Ad- hikari, V. Adya, et al., Astroph. J. Lett. 848, L13 (2017).
  * (25) J. M. Ezquiaga and M. Zumalacarregui, Front. Astron. Space Sci. 5, 44 (2018).
  * (26) B. Patricelli el al., MNRAS 513, 4159 (2022).
  * (27) B. P. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Ad- hikari, V. B. Adya, et al., Phys. Rev. Lett. 123, 011102 (2019).
  * (28) T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, Phys. Rev. Lett. 119, 251301 (2017).
  * (29) J. M. Ezquiaga and M. Zumalacarregui, Phys. Rev. Lett. 119, 251304 (2017).
  * (30) P. Creminelli and F. Vernizz, Phys. Rev. Lett. 119, 251302 (2017).
  * (31) J. Sakstein and B. Jain, Phys. Rev. Lett. 119, 251303 (2017).
  * (32) E. J. Copeland et al., Phys. Rev. Lett. 122, 061301 (2019).
  * (33) S. Boran et al, Phys. Rev. D 97, 041501 (2018).
  * (34) Y. Gong et al, Phys. Rev. D, 97, 084040 (2018).
  * (35) G. J. Olmo, D. Rubiera-Garcia, and A. Wojnar, Physics Reports 876, 1 (2020).
  * (36) B. Margalit and B. Metzger, Astrophys. J. 850, L19 (2017).
  * (37) M. Shibata, et al., Phys. Rev. D 96, 123012 (2017).
  * (38) M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D 97, 021501 (2018).
  * (39) L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. 852, L25 (2018).
  * (40) E. Fonseca, T. T. Pennucci, J. A. Ellis, I. H. Stairs, et al., Astrophys. J. 832, 167 (2016).
  * (41) H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. De- morest, et al., Nature Astronomy 4, 72 (2020).
  * (42) B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, et al., Astrophys. J. Letters 892, L3 (2020).
  * (43) M.W. Coughlin, M. W., et al., MNRAS 480, 3871 (2018).
  * (44) A. Kozak, K. Soieva, A. Wojnar, arXiv:2205.12812.
  * (45) D. Aguiar Gomes, A. Wojnar, arXiv:2206.04464.
  * (46) A. Wojnar, Phys. Rev. D 107, 044025 (2023).
  * (47) J.-P. Poirier, Introduction to the Physics of Earth’s Interior, Cambridge University Press, 2000.
  * (48) A. M. Dziewonski, D. L. Anderso, Preliminary reference Earth model, Phys. Earth Plan. Int. 25, 297 (1981).
  * (49) B. Kustowski et al, J. Geophys. Res. Solid Earth 113.B6 (2008).
  * (50) B. L. N. Kennett, E. R. Engdahl, Geophys. J. Int. 105.2, 429 (1991).
  * (51) B. L. N. Kennett, E. R. Engdahl, R. Buland, Geophys. J. Int. 122.1, 108 (1995).
  * (52) https://ds.iris.edu/ds/products/emc-referencemodels
  * (53) R. Butler, S. Tsuboi, Phys. Earth Planet. 321, 106802 (2021).
  * (54) D. A. Frost, & B. Romanowicz, Phys. Earth Planet. Inter. 286, 101 (2019).
  * (55) J. Stephenson, H. Tkalcic, M. J. & Sambridge, Geophys. Res. Solid Earth 126, e2020JB020545 (2021).
  * (56) T.-S. Pham, & H. Tkalcic, Nature Commun. 14, 754 (2023).
  * (57) S. Merkel, et al., Phys. Rev. Lett. 127.20, 205501 (2021).
  * (58) W. Winter, Walter, Neutrino Geophysics: Proceedings of Neutrino Sciences 2005. Springer, New York, NY, 2006. 285-307.
  * (59) A. Donini, S. Palomares-Ruiz, J. Salvado, Nature Physics 15.1, 37 (2019).
  * (60) S. Bourret, V. Van Elewyck, EPJ Web of Conferences. Vol. 207. EDP Sciences (2019).
  * (61) V. Van Elewyck, J. Coelho, E. Kaminski, L. Maderer, Europhysics News 52.1, 19 (2021).
  * (62) A. Kozak, A. Wojnar, Int. J. Geom. Meth. Mod. Phys. 19, 2250157 (2022).
  * (63) A. Kozak, A. Wojnar, Universe 8, 3, (2021).
  * (64) A. Wojnar, Phys. Rev. D 104, 104058 (2021).
  * (65) A. Wojnar, Phys. Rev. D 105, 124053 (2022).
  * (66) W.M. Kaula, An Introduction to Planetary Physics: The Terrestrial Planets, Wiley (1968).
  * (67) B.G. Bills, J. Geophys. Res. Planets 104.E12, 30773 (1999).
  * (68) W. M. Folkner et al., Science (1999), 1749 (1997).
  * (69) A. S. Konopliv and W. L. Sjogren, Publication 95 - 3, Jet Propulsion Laboratory, California Institute of Technology (1995).
  * (70) D.E. Smith et al., J. Geophys. Res. 98, 20871 (1995).
  * (71) W. M. Folkner et al., J. Geophys. Res. 102, 4057 (1997).
  * (72) J. G. Williams, Astroph. J. 108, 711 (1994).
  * (73) https://mars.nasa.gov/insight/spacecraft/instruments/seis/
  * (74) S. P. Weppner et al, MNRAS 452.2, 1375 (2015).
  * (75) S. Seager, et al., Astroph. j. 669.2, 1279 (2007).
  * (76) B. A. Bolt, Inside Earth, W. H. Freeman & Co., San Francisco (1982).
  * (77) K. E. Bullen, & B. A. Bolt, An Introduction to the Theory of Seismology, Cambridge University Press, Cambridge (1985).
  * (78) B. Luzum et al., Celest. Mech. Dyn. Astron. 110, 293 (2011).
  * (79) W. Chen, J. Ray, W. B. Shen, and C. L. Huang, J. Geod. 89, 179 (2015).
  * (80) F. Birch, Geophys. J. R. astr. Soc. 4, 295 (1961b).
  * (81) F. Birch, Solids under Pressure , W. Paul & D. M. Warschauer eds., McGraw Hill, New York (1963).
  * (82) V.I. Alfonso, et al., Class. Quantum Grav. 34, 235003 (2017).
  * (83) J. Beltran Jimenez, A. Delhom, EPJC 79, 1 (2019).
  * (84) J. Beltran Jiménez, A. Delhom, EPJC 80, 585 (2020).
  * (85) D. N. Vollick, Phys. Rev. D 69, 064030 (2004).
  * (86) J. Beltran Jimenez, L. Heisenberg, G. J. Olmo, D. Rubiera-Garcia, Phys. Rept. 727, 1 (2018).
  * (87) J. D. Toniato, D. C. Rodrigues, A. Wojnar, Phys. Rev. D 101, 064050 (2020).
  * (88) M. Banados, P.G. Ferreira, Phys. Rev. Lett. 105, 011101 (2010).
  * (89) P. Pani, V. Cardoso, T. Delsate, Phys. Rev. Lett. 107, 031101 (2011).
  * (90) G. J. Olmo, Phys. Rev. Lett. 95, 261102 (2005).
  * (91) P. Banerjee, et al., Astroph. J. 924, 20 (2022).
  * (92) S. Jana et al., Phys. Rev. D 97, 084011 (2018).
  * (93) J. Casanellas et al., Astroph. J. 745, 15 (2012).
  * (94) P. P. Avelino, Phys. Rev. D 85, 104053 (2012).
  * (95) S. Banerjee et al., JCAP 10, 004 (2017).
  * (96) A. Hernandez-Arboleda, D. C. Rodrigues, A. Wojnar (2022), arXiv:2204.03762.
  * (97) D. N. Krasnoshchekov, V. M. Ovtchinnikov, Dokl. Earth Sci. 478, 219 (2018).
  * (98) P. Pani, T. P. Sotiriou, Phys. Rev. Lett. 109, 251102 (2012).
  * (99) J-L. Margot et al, J. Geophys. Res. Planets 117.E12 (2012).
  * (100) B. Brugger et al, European Planetary Science Congress, EPSC2018-404, 2018.
  * (101) G. Steinbrugge et al, Geophys. Res. Lett. 48.3, e2020GL089895 (2021).
  * (102) H. Harder, G. Schubert, Icarus 151.1, 118 (2001)
  * (103) T. Spohn, F. Sohl, K. Wieczerkowski, V. Conzelmann, Planet. Space Sci. 49, 1561 (2001).
  * (104) M. A. Riner et al, J. Geophys. Res. Planetss 113.E8 (2008).
  * (105) J-L. Margot et al, Nature Astronomy 5, 676 (2021).
  * (106) A. S. Konopliv, C. F. Yoder, Geophys. Res. Lett. 23.14, 1857 (1996).
  * (107) K. C. Condie, Plate tectonics and crustal evolution, 4th ed., Butterworth-Heinemann (1997).
  * (108) I. Kulikov, P. Pronin, Int. J. Theor. Phys. 34, 1843 (1995).
  * (109) R. Verma, P. Nandi, Gen. Relativ. Gravit. 51, 143 (2019).
  * (110) S. Kalita, L. Sarmah, A. Wojnar, Phys. Rev. D 107, 044072 (2023).
  * (111) H-C. Kim, Phys. Rev. D 89, 064001 (2014).
  * (112) A. Wojnar, Phys. Rev. D 107, 044025 (2023).
  * (113) P. Lecca. Journal of Physics: Conference Series, volume 2090, page 012034. IOP Publishing, 2021.
  * (114) E. Lope-Oter, A. Wojnar, in preparation.
  * (115) P. P. Avelino, Phys.Rev. D 85, 104053 (2012).


On the singularity structure of a discrete modified-Korteweg-deVries equation

Basil Grammaticos Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France and Université de Paris, IJCLab, 91405 Orsay, France Thamizharasi
Tamizhmani SAS, Vellore Institute of Technology, Vellore - 632014, Tamil Nadu,
India Ralph Willox Graduate School of Mathematical Sciences, the University of
Tokyo, 3-8-1 Komaba, Meguro-ku, 153-8914 Tokyo, Japan

Abstract We study the singularities of a modified lattice Korteweg-deVries
(KdV) equation and show that it admits three families of singularities, with
analogous properties to those found in the lattice KdV equation. The first
family consists of localised singularities which can occupy an arbitrarily
large domain but which are, nevertheless, always confined. The second family
consists of one or more lines extending all the way from the south-west to the
north-east on the plane, involving a single finite value that depends on the
parameter that appears in the equation. We argue that the infinite extent of
this singularity is not incompatible with the confinement property or with the
integrability of the equation. The third family consists of horizontal strips
in which the product of values on vertically adjacent lattice sites is equal
to 1. In the case of the lattice KdV equation this type of singularity was
dubbed ‘taishi’. The taishi for the modified lattice KdV equation can interact
with singularities of the other two families, giving rise to very rich and
quite intricate singularity structures. Nonetheless, these interactions can be
described in a compact way through the formulation of a symbolic dynamics that
is similar to, but in a sense, simpler than that for the KdV case. We give an
interpretation of this symbolic dynamics in terms of a Box&Ball system related
to the ultradiscrete mKdV equation. These results show that taishi-type
singularities are not limited to the lattice KdV equation, but might very well
be a general feature of integrable lattice equations with deep connections to
other facets of their integrability.

Keywords: discrete integrable systems, lattice equations, singularity
structure, symbolic dynamics

1\. Introduction Singularities play a major role in physics and mathematics
[1]. As M. Kruskal used to point out, ‘singularities are where interesting
things happen’. In the domain of mathematical modelling of physical systems
the study of singularities turned out to be a most useful tool. Physical
phenomena can display, depending on conditions, either regular or irregular
behaviour. One thus expects that their mathematical models would exhibit the
same diversity of behaviour. This is indeed the case and the prognosis of
regular or irregular behaviour in dynamical systems is thus of particular
interest. The theory of chaos was in fact developed to describe a particular
type of irregular motion that precludes long-term predictions, often
exhibiting a crucial dependence on the details of the initial conditions.
Regular motion on the other hand is associated with predictability, allowing
for long-term predictions about the evolution of a given system.

Singularities are particularly effective in the prediction of the behaviour of
a given dynamical system [2]. In the case of differential equations, the
theory of Painlevé [3] links the singularity structure of a system to its
integrable (or non-integrable) character; integrability guaranteeing regular
and predictable behaviour. In the domain of discrete systems it turned out
that singularities are also important. The singularity confinement conjecture
[4] postulated that discrete equations integrable through spectral methods
must possess so-called confined singularities, where ‘confinement’ means that
any spontaneously appearing singularity again disappears after a finite number
of iteration steps.

The property of confinement of singularities was first observed [4] in a study
of the discrete (or lattice) KdV (d-KdV) equation [5]

$x_{m+1,n+1}=x_{m,n}+{1\over x_{m+1,n}}-{1\over x_{m,n+1}}.$ $None$

Figure 1 shows two examples of confined singularities for the d-KdV equation.

Figure 1

Although the d-KdV equation was instrumental in the discovery of the
singularity confinement property, the complete study of its singularities had
to wait for more than a quarter century. In fact, it was not until 2021 that
some of the present authors produced an exhaustive study [6, 7] of the
singularities of d-KdV and their interactions.

Among these singularities there are several that correspond to patterns of
infinite extent, as shown in Figures 2 and 3.

Figure 2

It is important to realize however that the existence of singularities of
infinite extent, such as those in Figure 2, does not invalidate the
confinement property as a requirement for integrability since it is impossible
to enter such a singularity coming from finite, regular, initial conditions.

The singularity shown in Figure 3 is of a completely new type, only discovered
recently [7]. Its origin becomes clear if we rewrite (1) as

${x_{m+1,n+1}x_{m+1,n}-1\over x_{m+1,n}}={x_{m,n}x_{m,n+1}-1\over x_{m,n+1}}.$
$None$

If at some point $(m,n)$ on the lattice the product $x_{m,n}x_{m,n+1}$ happens
to be equal to 1 then it will remain equal to 1 at the next point to the
right, i.e. at $(m+1,n)$. This singularity therefore consists in a horizontal
strip, running from west to east on the lattice, where at vertically adjacent
lattice points we have two values that are reciprocal. The corresponding
singularity pattern is shown in the upper part of Figure 3.

Figure 3

One can also imagine situations where more than one such strip exists, some of
them even being adjacent. In the lower part of Figure 3 we show such a
‘double’ strip, i.e. one where we have $x_{m,n}x_{m,n\pm 1}=1$. Let us remark
here that we could have reorganised equation (1) in a different way showing
the existence of vertical strips that satisfy the condition
$x_{m,n}x_{m+1,n}=-1$. In [7] we decided to dub the singularity associated to
a strip in which the product of two adjacent values is equal to 1 (or $-1$) a
‘taishi’. The term taishi derives from the Japanese naming of elementary
particles (which always have the suffix -shi) and the character for ‘band’ or
‘strip’ (‘tai’). Had the taishi been isolated singularities they would have
been of trifling importance. However, as we showed in [7] the taishi interact
with all the other singularities and these interactions lead to very
interesting results.

We shall illustrate this by an example, but before proceeding any further it
is useful to introduce a handy notation. In the following, when we write 0 or
$\infty$, what we really mean is that we introduce initial conditions in which
a small quantity appears, say $\epsilon$, for which we take the limit
$\epsilon\to 0$. Thus we take any zeros to be proportional to some power of
$\epsilon$, infinities to be proportional to some power of $1/\epsilon$ and
the values appearing in the taishi to be such that their product is equal to
$1+{\cal O}(\epsilon^{p})$ for some positive integer $p$. As explained in [7]
the unfettered richness of the singularities of d-KdV is only apparent if one
allows different singularities to correspond to different orders in
$\epsilon$. To simplify the notation we shall represent the different orders
in $\epsilon$ or $1/\epsilon$ as exponents $0^{1},0^{2},0^{3},\cdots$ and
$\infty^{1},\infty^{2},\infty^{3},\cdots$. In the case of the taishi, we shall
keep the $\epsilon$, and the notation $\epsilon^{p}$ means that the product
leading to the taishi is of order $1+{\cal O}(\epsilon^{p})$. In what follows
we shall refer to these exponents as ‘weights’.

In the light of this remark, the triple of lines in Figure 2 is to be
understood as representing the forward and backward evolution, according to
(1), of an initial condition that contains the singular part
$\\{\infty^{1},0^{1},\infty^{1}\\}$. (All other initial conditions are taken
to be generic and regular and can be given either on a staircase that includes
the aforementioned singular part or, for example, on the half-axes shown in
Figure 4).

One may thus wonder what would happen is one started from an initial condition
$\\{\infty^{1},0^{2},\infty^{1}\\}$ instead. This situation is depicted in
Figure 4. Looking at the patterns of 0 and $\infty$ in the (slightly greyed)
right-upper part of Figure 4, which depicts the forward evolution of the
initial data, one has the impression that a non-confined singularity pattern
arises from the leftmost infinity, whereas the middle zero and rightmost
infinity simply disappear. However this is not the case. Close inspection of
the values of $x_{m+3,n}$, $x_{m+3,n+1}$ and $x_{m+3,n+2}$ shows that their
products $x_{m+3,n}x_{m+3,n+1}$ and $x_{m+3,n+1}x_{m+3,n+2}$ are indeed of
order $1+{\cal O}(\epsilon)$ such that they take (finite) reciprocal values at
the limit $\epsilon\to 0$. This then indicates the existence of an outgoing
taishi. In order to find the incoming branch of this taishi we must iterate
backwards from the initial conditions.

Figure 4

This is shown in the lower left part of Figure 4. We find indeed that the
products $x_{m-1,n}x_{m-1,n-1}$ and $x_{m-1,n-1}x_{m-1,n-2}$ are again of
order $1+{\cal O}(\epsilon)$ and that $x_{m-1,n}$, $x_{m-1,n-1}$ and
$x_{m-1,n-2}$ take finite, reciprocal, values at the limit $\epsilon\to 0$.
Thus what at the outset looked like an ‘arbitrary’ initial condition, intended
to mimic the $\\{\infty^{1},0^{1},\infty^{1}\\}$ lines of Figure 2, actually
corresponds to an interaction of a line of $\infty^{2}$ with a double taishi
$\epsilon^{1},\epsilon^{1}$. The conclusion of this calculation then is that a
triple of infinite lines of the form $\\{\infty^{1},0^{2},\infty^{1}\\}$
cannot exist. However, this does not mean that infinite lines of the form
$\\{\infty^{p},0^{q},\infty^{r}\\}$ with $pqr\neq 1$ are impossible. It can be
verified that such lines do persist indefinitely if they satisfy the stability
condition $\min[p,q,r]=q$.

In [7] we presented a set of rules governing the dynamics of the interaction
of one or more taishi with a diagonal line of infinities for the d-KdV
equation. The basic assumption for this rule to work is that the line of
infinities has weight 1/2. The result of an interaction with a line of
infinities with general integer weight $q$ is thus obtained by repeating the
procedure below $2q$ times. Citing [7]:

Starting sufficiently low,

\- move upwards up to the first (i.e. lowest) horizontal strip in a taishi
with non-zero weight and subtract 1 from its weight and add 1 to the weight of
the strip just above it,

\- then move to the next horizontal strip with non-zero weight above those two
strips, subtract 1 from its weight and add 1 to the weight of the strip just
above it,

\- repeat the same procedure until there are no more taishi with strips with
non-zero weight left.

No rule for the interaction of the taishi with a line of zeros was provided in
[7]. However the answer is simple. In order to treat such lines of zeros it
suffices to follow the procedure described above while replacing ‘low’ by
‘high’, ‘upwards’ by ‘downwards’, ’subtract’ by ‘add’ and ‘above’ by ‘below’.
And of course one must keep in mind that an infinite line of zeros cannot
exist unless it is flanked on both sides by lines of infinities, as shown in
the bottom-right corner of Figure 2.

In the present paper we study the singularity structure of a lattice equation
that is linked to the lattice KdV equation by a set of Miura transformations
and which can therefore be regarded as a discrete modified KdV (mKdV)
equation. We show that although there exist many similarities with the case of
the lattice KdV equation described above, there are also essential
differences. In particular the formulation of the symbolic dynamics for the
taishi interactions in the case of mKdV is much simpler, allowing us to
entertain hopes for a future rigorous explanation for their dynamics.

2\. A discrete modified KdV equation The discrete form of the KdV equation was
first derived by Hirota in [5]. In his monumental work [8] where he proposed
the discrete analogue of the KP equation (an equation known since as the
Hirota-Miwa [9] equation) he also derived a discrete analogue of the modified
KdV equation but only gave its bilinear form, i.e. in terms of tau functions.
A nonlinear form of a potential discrete mKdV equation was proposed by Capel
and collaborators in [10] and Hirota obtained in [11] the nonlinear form of
the mKdV he had previously derived. As some of the present authors showed in
[12] this equation in fact coincides with that of Capel et al., after some
elementary transformations.

In [13] Levi and Yamilov derived what at the time seemed to be two new partial
difference equations. They obtained these equations from the combination of
Miura transformations linking the differential-difference Volterra equation to
a modified version thereof. They constructed the Lax pairs for both equations
and concluded that these systems are integrable. The two equations of Levi and
Yamilov were analysed in [14] by some of the present authors. We started by
obtaining their elementary singularities and derived their bilinear forms. The
singularity patterns we obtained were all confined (thus strengthening the
argument in favour of integrability) but the bilinear forms turned out to
coincide with those already known for the discrete KdV and mKdV equations.
Proceeding to the continuum limit we showed that the two equations of Levi and
Yamilov are indeed discrete analogues of the KdV and mKdV equations. The fact
that different discrete equations can have the same continuum limit is not
astonishing. The discrete integrable domain is far richer than the continuous
one and when taking a continuum limit of an integrable discrete system to a
continuous system, there are simply not that many possibilities for the
ensuing continuous integrable system. Unfortunately, this also means that
naming discrete systems based on their continuum limit is not quite
appropriate, but a better option is not always available. In what follows we
shall work with the discrete modified KdV equation proposed by Levi and
Yamilov in [13] and, for the sake of brevity, we shall refer to it as the
d-mKdV equation (hoping that there will be no confusion with the Hirota-Capel
lattice equation).

The discrete modified KdV equation we shall work with has the form

$(1+v_{m,n}v_{m+1,n})(kv_{m+1,n+1}+v_{m,n+1}/k)=(1+v_{m,n+1}v_{m+1,n+1})(kv_{m,n}+v_{m+1,n}/k),$
$None$

where $k$ is a non-zero constant such that $k^{4}\neq 1$ (lest the equation
become degenerate, as pointed out in [13]). We remark that, given the form of
(3), the equation is invariant under a sign change of either $v$ or $k$. An
equivalent form of (3) can be obtained by rearranging the terms on both sides:

$(1-v_{m,n}v_{m,n+1})(kv_{m+1,n+1}-v_{m+1,n}/k)=(1-v_{m+1,n}v_{m+1,n+1})(kv_{m,n}-v_{m,n+1}/k).$
$None$

Equation (3) being a discrete version of the modified KdV equation, it is
expected to be related to the d-KdV equation by a Miura transformation.
Curiously, this relation was not given in [14] and we remedy this here. The
Miura relations have the form

$x_{m,n}x_{m+1,n}={(k-v_{m,n})(k+v_{m+1,n})\over(kv_{m,n}-1)(kv_{m+1,n}+1)},$
$None$

$x_{m,n}x_{m,n+1}={(k-v_{m,n})(k+v_{m,n+1})\over(kv_{m,n}+1)(kv_{m,n+1}-1)}.$
$None$

Eliminating $v$ between these two relations gives relation (1) for $x$ and,
similarly, the elimination of $x$ from (5) and (6) leads to the d-mKdV
equation (3) above. To obtain equation (4) by Miura transformation from the
d-KdV equation, it suffices to change each occurrence of $k$ (or $v$) in (5)
and (6) to $-k$ (or $-v$).

In the following, all results we shall present are obtained from the evolution
under equation (3) in the positive $m$ and $n$ directions – i.e. towards the
NE on the lattice, as indicated in Figure 5 – of specific sets of initial
conditions given on a staircase.

The elementary singularity of d-mKdV was already obtained in [14]. Whenever,
due to a special choice of initial conditions, $v$ takes the value $k$ at a
point $(m,n)$ on the lattice, the values $1/k$ and $-1/k$ appear at points
$(m,n+1)$ and $(m+1,n)$, respectively, whereas $v_{m+1,n+1}$ takes the value
$-k$ and the singularity disappears. It is confined. Figure 5 below shows such
a singularity pattern.

Figure 5

When one compares this singularity pattern with the simplest one for d-KdV,
presented in Figure 1, it is tempting to identify a 0 in the d-KdV case with a
value $\pm k$ for d-mKdV, and similarly an $\infty$ for d-KdV with a $\pm 1/k$
in the present case. The results presented in this paper sustain this analogy,
up to a point, but as we shall explain this is by no means a strict
equivalence.

In the light of the results of [7] and the comments just above, we must also
look for other possible singularities. It turns out that an infinite line of
values $1/k$, running from SW to NE, can exist. It plays the same role as a
line of infinities for d-KdV. Figure 6 shows such a singularity.

Figure 6

It goes without saying that, given the symmetry of (3), an infinite line of
values $-1/k$ is also a possibility. In fact, several infinite lines can
coexist, and can even be adjacent, without any restrictions on the plus or
minus sign associated to each line.

The interesting question at this point, given the capital role they play in
the case of the d-KdV equation, is whether taishi-like singularities also
exist for the d-mKdV equation (3). It turns out that this is indeed the case.
Already by direct inspection of (3) one sees that if $v_{m,n}v_{m+1,n}=-1$,
for some value of $n$, the same relation automatically holds also at $n+1$.
Similarly, from (4), if $v_{m,n}v_{m,n+1}=1$ for some $m$ the same is true at
$m+1$. (To be fair, the reason we chose to work with the Levi-Yamilov form of
d-mKdV was precisely the fact that the taishi condition can be assessed by
simple inspection of the equation). Figure 7 shows a vertical and a horizontal
taishi.

Figure 7

It is also interesting to point out that the taishi conditions in the present
case can be obtained from those for the d-KdV equation through the Miura
transformations. Assuming for instance that $x_{m,n}x_{m+1,n}=1$ in (5), one
readily obtains the condition $v_{m,n}v_{m+1,n}=1$, and similarly in the other
direction.

Having presented the elementary singularities of d-mKdV, we now turn to the
study of more complex singularity patterns for which we shall see that the
above ‘simple’ correspondence between the singularities for the d-KdV and
d-mKdV equations breaks down.

In what follows, when we say that a singularity occurs at $k$, we mean that
the corresponding calculations are actually performed with the value
$k+\gamma\epsilon$ where $\gamma$ is a factor specific to each singularity
(and similarly for singularities at $1/k$ as well as for the taishi
conditions).

3\. The interaction of singularities of the discrete mKdV equation Before
addressing the question of interaction of the taishi with the other
singularities of the discrete mKdV equation, it is interesting to show the
singularity pattern generated by two adjacent values $v_{m,n+1}$ and
$v_{m+1,n}$ that are equal to $k$. The resulting confined singularity pattern
is shown in Figure 8.

Figure 8

Comparing this figure to Figure 1, obtained in the case of KdV, we remark that
while the values on the outside corners of the singularity pattern are, as
expected, a value $k$ where we would have had a 0 for d-KdV and a value $1/k$
where an infinity would have appeared for d-KdV, the singularities appearing
in the interior corners of the cross figure take neither of these two special
values and the analogy between the singular values for these two systems
clearly breaks down in this case. Still, the pattern on the interior corners
does follow the inversions and changes of sign of the elementary singularity
of d-mKdV, just as is the case for the d-KdV. (Note that a pattern similar to
that of Figure 8 can also be obtained in the case of d-KdV, but only when the
‘epsilons’ in the two starting singularities are taken to be strictly
identical).

The simplest taishi interaction one can think of is one where the taishi
impinges upon an elementary singularity, either head-on or tangentially.
Figure 9, obtained for initial conditions on a staircase that includes the
values $1/d$, $d$ and $k$, shows a head-on collision.

Figure 9

We remark that the elementary singularity is not perturbed by its encounter
with the taishi. On the other hand if the collision is off-centre but still
touches the elementary singularity pattern tangentially on its north side, as
in the case shown in Figure 10, a mirror-image of the singularity is created
above (i.e. to the north of) the original singularity pattern. The taishi
however exits the interaction region completely unperturbed. Note that Figure
10 was obtained from initial conditions on a staircase that includes the
values $1/d$, $d$ and the value $k$ to the SE of the value $d$, but none of
the other singular values.

Figure 10

Similar behaviours in the case of the taishi interacting with the elementary
singularity were observed for the d-KdV equation.

Next we turn to the interaction of a taishi with an infinite line of values
$1/k$. The result of such an interaction is presented in Figure 11, which was
obtained from initial conditions on a staircase that includes the values
$1/f$, $f$ and one value $1/k$ either to the bottom or to the right of $f$.

Figure 11

We remark that the result of this interaction is a shift of both the taishi as
well as of the infinite line. This situation is reminiscent of what happens in
the case of KdV, but the details of the interaction are quite different. In
the case of the d-KdV equation the interaction of a weight 1 taishi with a
line of infinities with that same weight results in an outgoing taishi shifted
upwards by two positions, whereas we only observe a single shift in the
present case.

In order to analyse these interactions in more detail we must relinquish the
tacit assumption that the difference from the value 1 for the product of
adjacent values in the case of the taishi, and from the value $1/k$ for the
infinite line, are both of first order in the small parameter $\epsilon$. As
explained in the introduction, we must introduce weights for both the taishi
and the infinite line. For lack of a more convenient notation we shall use, in
the case of $1/k$, a subscript to indicate the order in $\epsilon$ of the
perturbation of the value $1/k$ for $v$ . For instance, when we use the symbol
$1/k_{2}$ we mean this particular $v$ takes the value $1/k+{\cal
O}(\epsilon^{2})$. When no subscript is given, this is to be understood as
simply meaning $1/k+{\cal O}(\epsilon)$.

We first show an interaction of a taishi of weight 2 with an infinite line of
weight 1 (Figure 12).

Figure 12

We remark that the taishi is split into two parts, sharing equally the initial
weights, which is again quite different from the d-KdV case where the
impinging taishi would not split but would just be shifted upwards by one
position.

The situation when a taishi of weight 1 impinges upon an infinite line of
weight 2 is shown in Figure 13.

Figure 13

We could of course continue presenting more and more examples of taishi
interactions. However, the set of examples we have at our disposal already
suffices to formulate a symbolic dynamics that describes the interaction of a
taishi with an infinite line of $\pm 1/k$. This will be the subject of section
4.

Before proceeding any further, however, it is interesting to present yet one
more case where the analogy between the d-KdV and the d-mKdV equations breaks
down. As shown in Figure 2 the d-KdV equation possesses a singularity pattern
consisting in three infinite lines where a zero is interspersed between two
infinities. (And of course also more complicated patterns where infinities
alternate with zeros). As we saw in the previous section, there is a strong
resemblance between the d-KdV and d-mKdV singularities where, in the latter,
values $1/k$ and $k$ play the roles played by $\infty$ and 0 for the former
equation. One would thus be tempted to posit the existence of a triple of
infinite lines of singularities for d-mKdV, where a value of $k$ is
interspersed between two $1/k$ values. In order to investigate this
possibility we start with precisely such an initial condition, namely
$v_{m-1,n}=1/k,v_{m,n}=k$ and $v_{m+1,n}=1/k$. (The weights of all three are
taken to be equal to 1, i.e. $v_{m-1,n}=1/k+{\cal O}(\epsilon)$ etc.; we have
verified that choosing different weights does not alter the conclusions). The
result of the evolution of this specific initial condition is shown in the
Figure 14 below.

Figure 14

We observe that the initial condition $1/k,k,1/k$ does not survive. Looking at
the upper right (slightly greyed) corner of the figure, one sees that only a
singularity $1/k$ with weight 3 emerges, accompanied by a taishi of overall
weight 2. In fact, iterating backwards one sees that the initial condition we
introduced was just the result of the interaction of a taishi of weight 2 with
an infinite line $1/k$ of weight 3. This situation is of course strongly
reminiscent of the instability observed in the case of Figure 4, for d-KdV,
due to the choice of the weight of the zero in the middle of the initial
condition. In the case of the d-KdV equation however there exist infinitely
many choices of the weight of this zero that do lead to a stable singularity
pattern. It seems that this is not possible for the d-mKdV equation. The
analogy between the singularity patterns of the two equations is thus less
straightforward than one might have surmised.

4\. A symbolic dynamics for the taishi interactions In the previous section we
have presented several examples of interactions of a taishi with an infinite
line of $1/k$. While in these relatively simple examples it was easy to follow
the evolution of the taishi, this is by no means always the case. We therefore
need a rule that will allow us to calculate the shape of the taishi after it
emerges from the interaction region, as well as the shift of the infinite line
of $1/k$ due to the interaction. The answer to the latter question is simple:
every single interaction with a taishi – be it simple or double – results in
an upwards shift of the infinite line by two lattice steps.

The changes the interaction induces on the d-mKdV taishi themselves can be
described using symbolic dynamics similar to that for the d-KdV case, though
slightly simpler. In the case of the d-KdV equation the symbolic dynamics for
its taishi was formulated assuming a somewhat artifical weight of 1/2 for the
line of infinities, requiring extra intermediary steps which are now no longer
necessary. As we shall see, the rules for the symbolic dynamics for the d-mKdV
taishi can be formulated based on the natural assumption that the infinite
line it interacts with has weight 1. Keeping this essential difference in mind
we can now formulate the symbolic dynamics rules in exactly the same terms as
for the d-KdV equation:

Starting well-below the lowest taishi (in the $n$ direction),

\- move upwards up in $n$ to the first (i.e. lowest) horizontal strip in a
taishi with non-zero weight and subtract 1 from its weight and add 1 to the
weight of the strip just above it,

\- then move to the next horizontal strip with non-zero weight above those two
strips (in the $n$ direction), subtract 1 from its weight and add 1 to the
weight of the strip just above it,

\- repeat the same procedure until there are no more taishi with strips with
non-zero weight left.

This is the rule for the interaction of an arbitrary collection of taishi with
an infinite line of $1/k$ with weight 1. If the weight $q$ of the line of
$1/k$ is greater than 1, then the above procedure should simply be repeated
$q$ times.

The interaction rule for a collection of taishi with a weight 1 line of values
$1/k$ can also be rephrased in the language of so-called ‘Box&Ball systems’
[15]. Let us represent the weights of all the taishi, before interaction, as
entries in an infinite vector $U$ parallel to the $n$ axis on the lattice (the
vector $U$ only contains the weights of the taishi at each position in $n$,
before interaction, and neglects all weights associated to any infinite lines
the taishi will interact with) and think of the changes in weights described
by the above rule, as the result of an interaction with a carrier that moves
along the $n$ axis, from position $-\infty$ to $+\infty$, interacting with the
entries in $U$ at each successive position it passes through. This carrier has
maximal capacity 1, i.e. it can only take the values 0 and 1 (in which case we
say that it is empty or full, respectively) and it is always empty for
sufficiently negative positions in $n$. An empty carrier does not interact
with zero entries in the vector $U$ but if the carrier arrives empty at a
position in $U$ with a non-zero entry, then the value in $U$ at this position
decreases by one and the value of the carrier is updated to become 1. On the
other hand, a carrier that arrives full at a certain position decreases its
state to 0 and adds 1 to the value in $U$ at that position. Note that this
necessarily happens at each position just above a position where the carrier
was updated from empty to full. Hence, the above interaction rule can be
interpreted as the update rule for a Box&Ball system with carrier [15], where
$U$ represents an infinite array of boxes that can each contain an arbitrary
number of balls (the weights of the taishi) and where the carrier, $V$, has
capacity 1.

If we denote the updates of the values of the carrier by $\hat{V}$ and the
(local) update of the vector $U$ by $\bar{U}$, then the interaction rule at
each position $n$ in the vector $U$ can be expressed as

$V=1\quad\Rightarrow\quad\hat{V}=0,\quad\bar{U}=U+1,\hskip 42.67912pt$ $None$
$V=0\quad\Rightarrow\quad\hat{V}=\min[1,U],\quad\bar{U}=U+V-\hat{V},$ $None$

which can be summarized in the closed form

$\hat{V}=\min[U,1-V]\,,\quad\bar{U}=U+V-\hat{V},$ $None$

under the boundary conditions $V,U=0$ for $|n|\sim\infty$ (for a vector $U$
that only has non-negative integer entries). Readers familiar with Box&Ball
systems will immediately notice that the update rule (9), subject to said
boundary conditions, is in fact the same as that of the famous Takahashi-
Satsuma soliton cellular automaton [16]. We shall come back to this point in
section 5.

Using this rule, the dynamics corresponding to figures 11, 12 and 13 can now
be simply represented by the diagrams below:

$\matrix{0\cr 0\cr 0\cr\bf 1\cr 0\cr 0}\Longrightarrow\matrix{0\cr 0\cr\bf
1\cr 0\cr 0\cr 0}\hskip 71.13188pt\matrix{0\cr 0\cr 0\cr\bf 2\cr 0\cr
0}\Longrightarrow\matrix{0\cr 0\cr\bf 1\cr\bf 1\cr 0\cr 0}\hskip
71.13188pt\matrix{0\cr 0\cr 0\cr\bf 1\cr 0\cr 0}\Longrightarrow\matrix{0\cr
0\cr 1\cr 0\cr 0\cr 0}\Longrightarrow\matrix{0\cr\bf 1\cr 0\cr 0\cr 0\cr 0}$

in which only the (immediately relevant) entries of the vector $U$ are given,
and where numbers in boldface correspond to the weights in the taishi before
interaction (on the left) and after it has left the interaction region (on the
right).

The interaction shown in Figure 14 however requires some explanation. First,
let us consider the generic interaction of the taishi of Figure 14 with an
infinite line of weight 3, i.e. obtained for the evolution of initial
conditions (on a staircase) that include the values $1/f,f$ and $1/k^{3}$ (as
shown in Figure 15) but none of the other singular values.

Figure 15

The corresponding symbolic dynamics is given by:

$\matrix{0\cr 0\cr 0\cr\bf 1\cr\bf 1\cr 0}\Longrightarrow\matrix{0\cr 0\cr
0\cr 2\cr 0\cr 0}\Longrightarrow\matrix{0\cr 0\cr 1\cr 1\cr 0\cr
0}\Longrightarrow\matrix{0\cr 0\cr\bf 2\cr 0\cr 0\cr 0}$

where the numbers in boldface again represent the initial and final weights of
the impinging taishi. The interaction depicted above is of course clearly
different from that in Figure 14. However, one should not forget that the
initial conditions that lead to the interactions shown in Figure 14 were
chosen such as to have consecutive values $1/k,k$ and $1/k$ in the $m$
direction, instead of the values that would naturally arise from the taishi
interaction, namely $1/k,k_{2},g,h,\cdots$ (as in Figure 15). This means that
we have, in fact, enforced one extra interaction with what turns out to be a
diagonal with weight 1 at the point of interaction (as shown by the thick line
in Figure 16 below).

Figure 16

The presence of this (weight 1) line alters the symbolic dynamics to

$\matrix{0\cr 0\cr 0\cr\bf 1\cr\bf 1\cr 0}\Longrightarrow\matrix{0\cr 0\cr
0\cr 2\cr 0\cr 0}\Longrightarrow\matrix{0\cr 0\cr 1\cr 1\cr 0\cr
0}\Longrightarrow\matrix{0\cr 0\cr 2\cr 0\cr 0\cr
0}\Longrightarrow\matrix{0\cr\bf 1\cr\bf 1\cr 0\cr 0\cr 0}$

in perfect agreement with what was obtained in Figure 14.

Imposing special (singular) values $k$ or $1/k$ at specific lattice sites can
have a deep effect on the evolution, transforming the singularity patterns one
would otherwise obtain for less specific initial values. Suppose for instance
that in the cross pattern shown in Figure 8 we demand that the value $f$
appearing in the interior of the cross is exactly equal to $1/k$, in analogy
to what happens in the case of the d-KdV equation. The consequence of this
choice is a profound modification of the singularity pattern, as can be seen
in Figure 17.

The value $1/k$ we imposed instead of $f$ requires – thinking of the backward
evolution – the existence of a whole line of $1/k$ interacting with a taishi
of weight 1 in the SW direction. As a result, this taishi is shifted upwards
once and exits the interaction area as expected. The singularities around and
inside the cross follow the same pattern as in d-KdV, mutatis mutandis, but
the exiting taishi modifies the righthand side of the cross. So while a
parallel with d-KdV does exist up to a certain point, treating the value $k$
in d-mKdV as the strict analogue of the d-KdV zero does not seem to work,
apart from the simplest cases. This makes the study of the singularities of
d-mKdV all the more interesting.

Figure 17

5\. Conclusion and outlook In this paper we studied the singularities of the
discrete modified KdV equation and their interactions. This study was
motivated by an exhaustive investigation of the possible singularities for the
discrete KdV equation that we carried out recently [7] and which led to the
identification of three different types of singularities. The first type are
localised, ‘confined’, singularities. These were already identified in [4] and
their discovery led to the formulation of the hypothesis that integrable
systems possess confined singularities. (It was discovered later that this
property does not necessarily hold for systems integrable through
linearisation, although when the integrability is obtained by spectral methods
no violation of the singularity confinement property is known). The second
type consists of singularities with infinite extent, more precisely of
diagonals of infinities going in the south-west to north-east direction on the
lattice and, in some cases, of lines of infinities alternating with lines of
zeros. We must stress here again that the infinite extent of such a
singularity is not incompatible with integrability.

The most interesting singularities however are those of the third type, the
existence of which was completely unexpected. These singularities consist of
horizontal (or vertical) strips, single or multiple, in which the product of
two vertically (horizontally) adjacent values is 1 (or $-1$). The
particularity of this third type of singularity is that it can interact with
singularities of all the other types, leading to very rich singular behaviour
for the d-KdV equation. We decided to introduce a new name for this type of
singularity, opting for the term ‘taishi’, inspired by the Japanese word for
‘band’.

In view of these results the question that of course immediately arises is
whether these taishi are an anomaly. Is their existence specific to the d-KdV
equation or, might they be a feature that one should expect from other
integrable lattice equations as well? It was thus natural to focus on the best
candidate to study this problem, given its affinity to KdV: its modified
counterpart. For our study we chose a particular discrete form of the modified
KdV equation, first proposed by Levi and Yamilov [13]. The main advantage of
this choice was the fact that the existence of the taishi could be assessed at
a glance, given the form of the equation. It thus turned out that this d-mKdV
indeed possesses the same type of singularities as the d-KdV equation (with
the exception of the infinite lines of zeros flanked by infinite lines of
infinities which do not exist in the d-mKdV case). Moreover, the taishi still
interact with the other types of singularities, in ways reminiscent of the
d-KdV case. However, ‘reminiscent’ does not mean ‘identical’ and several key
aspects of the interactions for the d-mKdV singularities differ from those for
the d-KdV equation. In particular, the symbolic dynamics that describes the
interaction of a set of taishi with an inifinite line of singular values turns
out to be simpler than for the d-KdV equation.

As explained in section 4, the symbolic dynamics that describes the
interaction of a (finite) set of taishi with an infinite line of singular
values $\pm 1/k$ with weight $q$, is in fact equivalent to $q$ time steps in a
Box&Ball system with carrier [15], with infinite capacity for the boxes (the
number of balls in which represent the weigths of the taishi) and with a
carrier that can only carry 1 ball. This particular system is the special case
$L=+\infty$, $M=1$ of the system

$U_{j}^{t+1}=V_{j}^{t}+\max[0,U_{j}^{t}+V_{j}^{t}-M]-\max[0,U_{j}^{t}+V_{j}^{t}-L],$
$None$ $V_{j+1}^{t}=V_{j}^{t}+U_{j}^{t}-U_{j}^{t+1},$ $None$

where $U_{j}^{t}$ represents the number of balls contained at time $t$ in the
box with index $j$ and where $V_{j}^{t}$ represents the number of balls in the
carrier when, at time $t$, it arrives at box $j$. Under the boundary
conditions $U_{j}^{t},V_{j}^{t}=0$ for large enough $|j|$ (at all $t$) and for
non-negative integer initial conditions $U_{j}^{0}$, this system describes the
evolution of a finite number of balls ($\sum_{j=-\infty}^{+\infty}U_{j}^{0}$)
over an infinite array of boxes – each with fixed capacity $L$ – due to the
movement of a carrier with capacity $M$ in the (positive) $j$ direction [15].
A simple calculation shows that the particular case $L=+\infty$, $M=1$ can
also be written as (9), which under the same boundary conditions actually
corresponds to the rule that defines one time step in the Takahashi-Satsuma
soliton cellular automaton [16], but in this case for the carrier variable $V$
and not for the $U$ variable that represents the boxes. This correspondence
between two apparently very different systems is a simple consequence of the
self-duality of the system $(10,11)$ under the involution
$(U,V;L,M;t,j)\to(V,U;M,L;j,t)$.

More important however is that the system $(10,11)$ is, in fact, an
ultradiscrete (or tropical) version of a (so-called Harrison-type) Yang-Baxter
map given in [17], which is known to define the Hirota-Capel discrete version
of the modified KdV equation [18]. We therefore arrive at the striking
conclusion that the symbolic dynamics that describe the interaction of the
taishi with type-two singularities for our d-mKdV equation (and, by extension,
also for the d-KdV equation) is actually equivalent to a tropical Yang-Baxter
map (or combinatorial $R$-matrix [19]) that governs the algebraic
integrability of the ultradiscrete version of the mKdV equation (or, depending
on your point of view, of the quantum integrable system with the same
$A_{1}^{(1)}$ symmetry, obtained through crystallization [20]). The fact that
soliton interactions define Yang-Baxter maps has been known for quite a while
[21,22] but that interactions of singularities are similarly related is quite
remarkable.

The results presented in the present paper thus not only confirm that the
behaviour we observed in the KdV case was not a singular occurrence (no pun
intended), but also point at a possibly deep connection between analytic
aspects (singularities) and algebraic aspects (Yang-Baxter maps) of
integrability in lattice equations. Besides the exploration of this intriguing
connection, we are also naturally led to speculate upon other possible paths
for future studies. The first question is whether the singularity structures
observed here also exist in other integrable lattice equations, and under
which guises. After all, KdV and m-KdV are like siblings and observing similar
behaviour in both is not all that astonishing. But even restricting ourselves
to the KdV case, the fact that the symbolic dynamics of the d-mKdV
singularities are very simple raises the hope of possibly finding a rigorous
proof for their formulation. In both [7] and the present paper the symbolic
dynamics were deduced by the observation of a slew of cases, but a rigorous
proof might also bring us closer to understanding the link to Box&Ball systems
(possibly a very tall order). A different direction one could consider is that
of establishing singular rules for these equations. Singular rules have been
introduced in numerical analysis [23,24] and, in particular, in convergence
accelerators in order to deal with accidentally vanishing denominators. The
existence of such rules is intimately related to the singularity confinement
property. It would thus be interesting to study the possible existence of such
rules for the various singularities we have discovered so far.

Acknowledgements RW would like to acknowledge support from the Japan Society
for the Promotion of Science (JSPS), through JSPS grant number 18K03355. He
also wishes to thank K. Nakata for useful discussions on the Box&Ball system
of section 4.

References

[1] A. Goriely, ‘Integrability and Nonintegrability of Dynamical Systems’,
World Scientific, Singapore (2001).

[2] A. Ramani, B. Grammaticos and T. Bountis, Physics Reports 180 (1989) 159.

[3] M.J. Ablowitz, A. Ramani and H. Segur, Lett. Nuov. Cim. 23 (1978) 333.

[4] B. Grammaticos, A. Ramani and V. Papageorgiou, Phys. Rev. Lett. 67 (1991)
1825.

[5] R. Hirota, J. Phys. Soc. Japan 43 (1977) 1424.

[6] D. Um, R. Willox, B. Grammaticos and A. Ramani, J. Phys. A: Math. Theor.
53 (2020) 114001.

[7] D. Um, A. Ramani, B. Grammaticos, R. Willox and J. Satsuma, J. Phys. A:
Math. Theor. 54 (2021) 095201.

[8] R. Hirota, J. Phys. Soc. Japan 503 (1981) 785.

[9] T. Miwa, Proc. Japan Acad. 58 (1982) 9.

[10] H. Capel, F.W. Nijhoff and V. Papageorgiou Phys. Lett A 153 (1991) 377.

[11] R. Hirota, J. Phys. Soc. Japan 672 (1998) 234.

[12] B. Grammaticos, A. Ramani, C. Scimiterna and R. Willox, J. Phys. A: Math.
Theor. 44 (2011) 152004

[13] D. Levi and R.I. Yamilov, “On a nonlinear integrable difference equation
on the square 3D-inconsistent”, preprint (2009) arXiv:0902.2126 [nlin.SI].

[14] A. Ramani, B. Grammaticos, R. Willox and J. Satsuma, J. Phys. A: Math.
Theor. 42 (2009) FT 282002.

[15] D. Takahashi and J. Matsukidaira, J. Phys. A: Math. Gen. 30 (1997) L733.

[16] D. Takahashi and J. Satsuma, J. Phys. Soc. Japan 59 (1990) 3514.

[17] V. G. Papageorgiou, A. G. Tongas and A. P. Veselov, J. Math. Phys. 47
(2006), 083502.

[18] S. Kakei, J.J.C. Nimmo and R. Willox, Glasgow Math. J. 51A (2009) 107.

[19] G. Hatayama, K. Hikami, R. Inoue, A. Kuniba, T. Takagi and T. Tokihiro,
J. Math. Phys. 42 (2001) 274.

[20] R. Inoue, A. Kuniba and T. Takagi, J. Phys. A: Math. Theor. 45 (2012)
073001.

[21] A.P. Veselov, Phys. Lett. A 314 (2003) 214.

[22] V.M. Goncharenko and A.P. Veselov, “Yang-Baxter maps and matrix
solitons”, in ‘New Trends in Integrability and Partial Solvability’, A.B.
Shabat et al. (eds.), Springer, Dordrecht (2004) pp. 191.

[23] P. Wynn, BIT 3 (1963), 175.

[24] C. Brezinski, ‘Algorithmes d’accélération de la convergence: étude
numérique’, Technip, Paris (1978).