# To the question of the integration of Plebansky Equation

A. N. Leznov111e-mail<EMAIL_ADDRESS>Universidad Autonoma del Estado de
Morelos, CCICAp,Cuernavaca, Mexico

###### Abstract

It is shown that corresponding to Plebansky equation of symmetry posses the
infinite set of solutions, which we present in explicit form. This fact leads
to conclusion about possibility to find series solutions of the Plebansky
equation in analytic form. Some classes of explicit solution are presented.

##  1 Introduction

The goal of the present paper is to investigate Plebansky equation on the
subject of its integrability. The Plebansky equation is non homogeneous
complex Monge-Ampher equation

$Det_{n}(\phi_{x_{i},\bar{x}_{j}})=1$ (1)

in the case $n=2$. $\phi$ is unknown function of $2n$ independent variables
$x_{i},\bar{x}_{j}$.

Equation of Plebansky arise in self-dual gravity [1] and we will use it in the
form

$v_{y,\bar{y}}v_{z,\bar{z}}-v_{y,\bar{z}}v_{z,\bar{y}}=1$ (2)

which will call as a basic form of this equation.

At this place we would like to notice that the general solution of homogeneous
Monge-Ampher equation (with zero instead of the unity in the right side of
(1)), was recently found in implicit form in [3].

We remind to the reader some known facts from the theory of differential
equations [2]. With each system of differential equations it is connected
symmetry system of equations, which arise after differentiation of the initial
system by some parameter and denotation such obtained derivative $\dot{a}=A$
as a new unknown function. If linear system of such equations have exact
solution on the class of solutions of the initial system the last one is
exactly integrable. This means that it is possible to present its solution (in
explicit or implicit form) depending on necessary number of arbitrary
functions sufficient for the statement the problem of Cauchy-Kovalevskai. Such
solution is called as the general one. If there exist series solutions of the
symmetry equation then initial system poses the same number of exact
solutions, but each solution is not the general one.

In the present paper we would like to show that in the case of Plebansky
equation second possibility take place. It may be found infinite series
solutions of its symmetry equation but not the general one. Thus it is not
possible to find general solution of this equation in analytic form. But
infinite number of its exact solutions exist and the goal of the present paper
consists in description the ways for obtaining them.

##  2 Preliminary manipulations

Now we would like to rewrite (2) in equivalent form

$(v_{\bar{y}}v_{z,\bar{z}}-v_{\bar{z}}v_{z,\bar{y}})_{y}+(v_{y,\bar{y}}v_{\bar{z}}-v_{y,\bar{z}}v_{\bar{y}})_{z}=2$
(3)
$(v_{y}v_{z,\bar{z}}-v_{z}v_{\bar{z},y})_{\bar{y}}+(v_{y,\bar{y}}v_{z}-v_{\bar{y},z}v_{y})_{\bar{z}}=2$
(4)

The last equations (3) and (4) can be partially resolved as

$v_{\bar{y}}v_{z,\bar{z}}-v_{\bar{z}}v_{z,\bar{y}}=\theta_{z}+y,\quad
v_{y,\bar{y}}v_{\bar{z}}-v_{y,\bar{z}}v_{\bar{y}}=-\theta_{y}+z$
$v_{y}v_{z,\bar{z}}-v_{z}v_{\bar{z},y}=\theta_{\bar{z}}+\bar{y},,\quad
v_{y,\bar{y}}v_{z}-v_{\bar{y},z}v_{y}=-\theta_{\bar{y}}+\bar{z}$

It is important to notice that function $\theta$ satisfy the the symmetry
equation corresponding to Plebansky one [2].

Resolving the last equalities with respect to $v_{\bar{y}},v_{\bar{z}}$ (or to
$v_{y},v_{z}$), taking into account the Plebansky equation (2) for $v$ we
obtain

$v_{\bar{y}}=v_{y,\bar{y}}\theta_{z}-v_{z,\bar{y}}\theta_{y}+zv_{z\bar{y}}+yv_{y,\bar{y}},\quad
v_{\bar{z}}=v_{y,\bar{z}}\theta_{z}-v_{z,\bar{z}}\theta_{y}+zv_{z\bar{z}}+yv_{y,\bar{z}}$
(5)
$v_{y}=v_{y,\bar{y}}\theta_{\bar{z}}-v_{y,\bar{z}}\theta_{\bar{y}}+\bar{z}v_{y,\bar{z}}+\bar{y}v_{y,\bar{y}},\quad
v_{z}=v_{z,\bar{y}}\theta_{\bar{z}}-v_{z,\bar{z}}\theta_{\bar{y}}+\bar{z}v_{z\bar{z}}+\bar{y}v_{z,\bar{y}}$
(6)

Equating second mixed derivatives of $v$ with respect to bar variables (or not
to the bar ones), we conclude that function $\theta$ in both cases is the
solution of the symmetry equation corresponding to (2) namely

$v_{\bar{y},y}\theta_{z,\bar{z}}+\theta_{\bar{y},y}v_{z,\bar{z}}-v_{y,\bar{z}}\theta_{\bar{y},z}-\theta_{y,\bar{z}}v_{\bar{y},z}=0$
(7)

##  3 Recurrent formula for solution of the symmetry equation

Rewriting (7) in two equivalent forms

$v_{\bar{y},y}\theta_{z}-v_{\bar{y},z}\theta_{y}=\tilde{\theta}_{\bar{y}},\quad
v_{\bar{z},y}\theta_{z}-v_{\bar{z},z}\theta_{y}=\tilde{\theta}_{\bar{z}}$ (8)
$v_{\bar{y},y}\theta_{\bar{z}}-v_{\bar{z},y}\theta_{\bar{y}}=\tilde{\theta}_{y},\quad
v_{\bar{y},z}\theta_{\bar{z}}-v_{\bar{z},z}\theta_{\bar{y}}=\tilde{\theta}_{z}$
(9)

By the same way as it was done in the previous section it is possible to check
that $\tilde{\theta}$ functions from (8),(9) satisfy the symmetry equation
(7). Obvious solution of symmetry equation is derivatives of $v$ with respect
to 4 independent arguments of the problem $y,z,\bar{y},\bar{z}$. Thus with
help of (8),(9) it is possible to construct infinite serie solutions of
symmetry equation (compare with [LYM]).

###  3.1 Static case

In this subsection we would like explain the case when symmetry equation has a
exact solution and how this connected with the integrable property of the
initial system.

Let us consider ”time independent” configurations, when
$v=v(z+\bar{z},y,\bar{y})$. In this case the Plebanski equation (2) and main
symmetry equation (5) are reduced correspondingly

$v_{y,\bar{y}}v_{z,z}-v_{y,z}v_{z,\bar{y}}=1$
$v_{\bar{y},y}\theta_{z}-v_{\bar{y},z}\theta_{y}=\tilde{\theta}_{\bar{y}},\quad
v_{z,y}\theta_{z}-v_{z,z}\theta_{y}=\tilde{\theta}_{z}$

Let us sick solution of the symmetry equation in terms of three independent
coordinates $(v_{z},y,\bar{y})$. The last symmetry system in these coordinates
looks as

$\theta_{v_{z}}v_{z,\bar{y}}+\theta_{\bar{y}}=\tilde{\theta}_{v_{z}}-v_{z,\bar{y}}\tilde{\theta},\quad\theta_{v_{z}}v_{z,z}=-v_{z,z}\tilde{\theta}_{y}$

From which follows that functions $\tilde{\theta},\theta$ are connected by the
linear system of equations

$\tilde{\theta}_{\bar{y}}=\theta_{v_{z}},\quad\tilde{\theta}_{v_{z}}=-\theta_{y}$

which is equivalent to three dimensional Laplace equation. Thus in the static
case the symmetry equation possesses the general solution and thus in this
case Plebansky equation is exactly integrable.

Below we present its general solution in implicit form

$v_{\bar{y}}=L_{v_{z}}(v_{z},y,\bar{y})\quad z=L_{y}(v_{z},y,\bar{y})$ (10)

where function $L$ satisfy Laplace equation in three dimension
$L_{v_{z},v_{z}}+L_{y,\bar{y}}=0$. General solution of Laplace equation in
three dimension depend on two arbitrary functions of two arguments. And thus
constructed above is a general solution of Plebansky equation in static case.

This solution may be obtained directly. Indeed Plebansky equation in this case
may be rewritten as

$(v_{\bar{y}})_{y}(v_{z})_{z}-(v_{z})_{y}(v_{\bar{y}})_{z}=1$

and functions $(v_{\bar{y}},v_{z})$ may be considered as transformed impulse
and coordinate with unity Poisson bracket between them. This system is
resolved by canonical transformation with generating function
$G(v_{z},y,\bar{y})$ by usual formulae

$v_{\bar{y}}=G_{v_{z}}(v_{z},y,\bar{y}),\quad z=G_{y}(v_{z},y,\bar{y})$

From condition of equality of the second mixed derivatives of $v$ function we
conclude that $G$ satisfy Laplace equation in three dimension. This is exactly
(10) above.

##  4 General strategy and example explaining it

General (5) is the system of two equations on two unknown functions
$v,\theta$. As a consequence we know that $v$ function satisfy Plebansky
equation and $\theta$ function the symmetry one. But all solutions of symmetry
equation are enumerated in the previous section. And thus if we change
$\theta$ in (5) on one of solution of the previous section we will obtain the
self consistent system of two equations only on one unknown function $v$.
Solution of this system (if it will be possible to find it) will be solution
of Plebansky equation corresponding to such chose solution of symmetry
equation.

To show that such idea is not mean less we at first consider simple example of
(5) under the chose $\theta=v_{\bar{y}}$. But result of solution of Plebansky
equation will be absolutely non trivial one.

We rewrite (5)

$v_{\bar{y}}=(z\frac{\partial}{\partial z}+y\frac{\partial}{\partial
y})v_{\bar{y}},\quad v_{\bar{z}}=1+(z\frac{\partial}{\partial
z}+y\frac{\partial}{\partial y})v_{\bar{z}}$

Two independent ordinary differential equations have obvious solution

$v_{\bar{y}}=(yz)^{1\over 2}X_{\bar{y}}(d,\bar{z},\bar{y}),\quad
v_{\bar{z}}=1+(yz)^{1\over 2}X_{\bar{z}}(d,\bar{z},\bar{y}),\quad
d\equiv{z\over y}$

Further

$v_{\bar{y},y}={1\over 2}d^{1\over 2}X_{\bar{y}}-d^{3\over
2}X_{\bar{y},d},\quad v_{\bar{y},z}={1\over 2}d^{-1\over
2}X_{\bar{y}}+d^{1\over 2}X_{\bar{y},d}$ $v_{\bar{z},y}={1\over 2}d^{1\over
2}X_{\bar{z}}-d^{3\over 2}X_{\bar{z},d},\quad v_{\bar{z},z}={1\over
2}d^{-1\over 2}X_{\bar{z}}+d^{1\over 2}X_{\bar{z},d}$

And equation of Plebansky (2) takes the form

$X_{\bar{y}}X_{\bar{z},D}-X_{\bar{z}}X_{\bar{y},D}=1,\quad D\equiv\ln d$ (11)

Resolving of the last equation is connected with second order ordinary
differential equations of the form $X_{D,D}=F(X,D)$ ($F$ arbitrary functions
of its arguments). Indeed solution of this equation depends on two arbitrary
parameters $c_{1},c_{2}$ and may be represented as some function depending on
3 arguments $X=X(D,c_{1},c_{2})$. Let us differential equation for $X$
function argumentson parameters $c_{i}$. We have
$X_{D,D,c_{i}}=F_{X}(X,D)X_{c_{i}}$. From the last equality we conclude
$(X_{D,c_{1}}X_{c_{2}}-X_{D,c_{2}}X_{c_{1}})_{D}=0$. We will assume that this
value is equal to unity (this is always possible to do by corresponding
canonical transformation). Now we will consider $c_{1},c_{2}$ as arbitrary
functions of the arguments $\bar{y},\bar{z}$. Then equation (11) Looks as

$(X_{c_{1}}X_{c_{2},D}-X_{c_{2}}X_{c_{1},D})((c_{1})_{\bar{y}}(c_{2})_{\bar{z}}-(c_{1})_{\bar{z}}(c_{2})_{\bar{y}})=((c_{1})_{\bar{y}}(c_{2})_{\bar{z}}-(c_{1})_{\bar{z}}(c_{2})_{\bar{y}})=1$

Thus obtained solution of initial Plebansky equation depend on two arbitrary
function $F(X,D)$ defined arbitrary equation of the second order and
generating function of canonical transformation resolving the last equation
for $c_{1},c_{2}$ functions.

We present the second way for solution equation (11). Let us consider
$X_{D},X$ as canonical transformed coordinate and impulse variable and
$\bar{z},\bar{y}$ as the same initial ones. Then generating function of
canonical transformation $W=W(X,\bar{y},D)$ satisfy the equations

$X_{D}=W_{X}(X,\bar{y},D),\quad\bar{z}=W_{\bar{y}}(X,\bar{y},D)$

From the second equation we have $X_{D}=-{W_{\bar{y},D}\over W_{\bar{y},X}}$
and after substitution into the first equation we obtain

$(W_{D}+{W_{X}^{2}\over 2})_{\bar{y}}=0,\quad W_{D}+{W_{X}^{2}\over 2}=F(X,D)$

The second one is Hamilton-Jacobi equation of the particle motion in potential
field $V=F(X,D)$, It leads to second order differential equation considered
above. Hamilton-Jacobi equation reduce the number of independent variables in
generating function $W$ from 3 up to 2. And thus solution of Plebansky
equation is determined by two functions $W,F$ each one of two independent
variables.

##  5 Equation of Plebansky in non usual variables

Let us introduce notations $R=\ln v_{\bar{y}}+\ln v_{\bar{z}},\Delta=\ln
v_{\bar{y}}-\ln v_{\bar{z}}$. Or in other words
$v_{\bar{y}}=\exp({R+\Delta\over 2}),v_{\bar{z}}=\exp({R-\Delta\over 2})$. In
these notations the pair of equations (LABEL:LS) take the form

$R_{y},\theta_{z}-R_{z}\theta_{y}+zR_{z}+yR_{y}=2,\quad\Delta_{y}\theta_{z}-\Delta_{z}\theta_{y}+z\Delta_{z}+y\Delta_{y}=0$
(12)

In the last equations let us pass from independent variables $y,z$ to
independent variables $\theta,d={z\over y}$. Corresponding necessary formulae
are presented below

$y=Y(\theta,d,\bar{y},\bar{z}),\quad 1=Y_{\theta}\theta_{y}-Y_{d}{z\over
y^{2}},\quad 0=Y_{\theta}\theta_{z}+Y_{d}{1\over y}$ (13)
$R_{y}==R_{\theta}\theta_{y}-R_{d}{z\over y^{2}},\quad
R_{z}=R_{\theta}\theta_{z}+R_{d}{1\over y}$

and the same formulae for derivatives of the $\Delta$ function. In all
relations above its necessary to keep in mind that all function under
consideration depend also on bar arguments $\bar{y},\bar{z}$. Such dependence
will be taken into account on some forward steps of calculations.

In variables $\theta,d$ the system (12) looks as

$(Y^{2})_{\theta}-Y^{2}R_{\theta}=-R_{d},\quad(Y^{2}e^{-R})_{\theta}=(e^{-R})_{d},\quad
Y^{2}={\Delta_{d}\over\Delta_{\theta}}$ (14)

To do the last equality more symmetrical to $y,z$ variables let us multiply
the last one on $d={z\over y}$. We have in a consequence

$dY^{2}=yz={d\Delta_{d}\over\Delta_{\theta}}={\Delta_{D}\over\Delta_{\theta}},\quad
D=\ln d$ (15)

Now we would like to satisfy equation of Plebansky

$v_{y,\bar{y}}=({R_{y}+\Delta_{y}\over 2})\exp({R+\Delta\over 2}),\quad
v_{z,\bar{y}}=({R_{z}+\Delta_{z}\over 2})\exp({R+\Delta\over 2})$
$v_{y,\bar{z}}=({R_{y}-\Delta_{y}\over 2})\exp({R-\Delta\over 2}),\quad
v_{z,\bar{z}}=({R_{z}-\Delta_{z}\over 2})\exp({R-\Delta\over 2})$

And thus the equation of Plebansky looks as

$e^{R}(R_{z}\Delta_{y}-R_{y}\Delta_{z})=2$

In the last equation let us pass to variables $\theta,d$ with the help of
above formulae. We obtain

$R_{d}\Delta_{\theta}-R_{\theta}\Delta_{d}=e^{-R}(Y^{2})_{\theta}$

Let us compare the last equation with obtained above (14) ones. As a direct
corollary we have

$e^{-R}=-\Delta_{\theta}$ (16)

The last relations have as a direct consequence two first relations (14).

In all calculations above no information about dependence of all functions
involved with respect to bar arguments was not used. Now let use condition of
equality of the second mixed derivatives of $v$ function with respect to
$\bar{y},\bar{z}$ arguments in notations introduced above.

$(\exp({R+\Delta\over 2}))_{\bar{z}}=(\exp({R-\Delta\over
2}))_{\bar{y}},\quad\exp\Delta(R_{\bar{z}}+\Delta_{\bar{z}})=(R_{\bar{y}}-\Delta_{\bar{y}})$

But functions $R,\Delta$ are connected by (16) and function $\Delta$ depends
in its turn on $\theta$ and thus terms with derivatives on bar arguments look
as (we present left hand side term)

$-(\ln\Delta_{\theta})_{\bar{z}}+\Delta_{\bar{z}}+(-(\ln\Delta_{\theta})_{\theta}+\Delta_{\theta})\theta_{\bar{z}}$
(17)

After some computations using (15) we come to equation,which function $\Delta$
satisfy

$e^{5\Delta}((e^{-\Delta})_{\theta,\bar{z}}(e^{-\Delta})_{\theta,d}-(e^{-\Delta})_{\theta,\theta}(e^{-\Delta})_{\bar{z},d})=(e^{\Delta})_{\theta,\bar{y}}(e^{\Delta})_{\theta,d}-(e^{\Delta})_{\theta,\theta}(e^{\Delta})_{\bar{y},d}$
(18)

In the case when symmetry function depends only from non bar variables
equation (18) looks much more simple (in (17)
$\theta_{\bar{z}}=\theta_{\bar{y}}=0$)

$e^{3\Delta}(e^{-\Delta})_{\theta,\bar{z}}=(e^{\Delta})_{\theta,\bar{y}}$

This equation is equivalent to equation considered in section 4. It will be
explained in one of sections below.

After introduction new function $\theta=\Theta(\Delta,d,\bar{y},\bar{z})$ it
looks as

$e^{\Delta}(\Theta_{\Delta,\Delta}\Theta_{d,\bar{z}}-\Theta_{\Delta,d}\Theta_{\Delta,\bar{z}}+\Theta_{\Delta}\Theta_{d,\bar{z}})=\Theta_{\Delta,\Delta}\Theta_{d,\bar{y}}-\Theta_{\Delta,d}\Theta_{\Delta,\bar{y}})-\Theta_{\Delta}\Theta_{d,\bar{y}}$
(19)

or in variables $d=D(\Delta,\theta,\bar{y},\bar{z})$ it looks as

$e^{\Delta}Det_{3}\pmatrix{D_{\bar{z}}&D_{\Delta}&D_{\theta}\cr
D_{\Delta,\bar{z}}&D_{\Delta,\Delta}+D_{\Delta}&D_{\Delta\theta}\cr
D_{\theta,\bar{z}}&D_{\theta,\Delta}&D_{\theta,\theta}\cr}=Det_{3}\pmatrix{D_{\bar{y}}&D_{\Delta}&D_{\theta}\cr
D_{\Delta,\bar{y}}&D_{\Delta,\Delta}-D_{\Delta}&D_{\Delta\theta}\cr
D_{\theta,\bar{y}}&D_{\theta,\Delta}&D_{\theta,\theta}\cr}$ (20)

Each of three equations above are equivalent to the initial Plebansky equation
(2).

###  5.1 Equation (19) in the integral motion form

Let us introduce the following notation the operators of the differentiations

$L^{\pm}=e^{{\Delta\over 2}}{\frac{\partial}{\partial\bar{z}}}\pm
e^{-{\Delta\over 2}}{\frac{\partial}{\partial\bar{y}}},\quad
L^{0}={\frac{\partial}{\partial\Delta}}$

with the obvious commutation relations

$[L^{0},L^{\pm}]={1\over 2}L^{\mp},\quad[L^{+},L^{-}]=0$

(these are commutation relation of algebra of two dimension plane). In these
notations equation (19) may be rewritten as

$\Theta_{\Delta,\Delta}(L^{-}\Theta)_{d}+2\Theta_{\Delta}(L^{-}\Theta_{d})_{\Delta}-2\Theta_{\Delta}(L^{-}\Theta_{\Delta})_{d}-\Theta_{\Delta,d}(L^{-}\Theta_{\Delta})=0$

Or after dividing on $(\Theta_{\Delta})^{{1\over 2}}$ and trivial regrouping
of the terms we obtain

$((\Theta_{\Delta})^{{1\over
2}}(L^{-}\Theta_{d}))_{\Delta}=((\Theta_{\Delta})^{{1\over
2}}(L^{-}\Theta_{\Delta}))_{d}$ (21)

Or (19) may be rewritten in integral of motion form. We remind the reader that
all equations (18),(19), (20) where obtained from the equality of the
derivatives on the bar variables.

##  6 Some examples of particular solutions

In this section we present some particular solutions of equations (19) or
equivalent to it (18), (20). From this consideration it will be clear that
these equation leads indeed to solution of initial Plebansky equation.

###  6.1 The case when symmetry function do not depend on bar variables

From explicit form of symmetry equation of the second section it is clear that
simplest obvious its solution is $\theta=\theta(y,z)$. Equation describing
this situation are the following ones

$e^{3\Delta}(e^{-\Delta})_{\theta,\bar{z}}=(e^{\Delta})_{\theta,\bar{y}},\quad{\Delta_{d}\over\Delta_{\theta}}=Y^{2}(\theta,D)=yz$

From the second equation that function $e^{\Delta}$ really is the function of
only three variables $\bar{y},\bar{z},g(\theta,D)$. Let us seek solution of
the first equation in a form

$e^{\Delta}=-{X_{\bar{y}}\over X_{\bar{z}}}$

where $X=X(\bar{y},\bar{z},g(\theta,D))$. After such substitution first
equation takes the form

$({X_{\bar{y}}\over
X_{\bar{z}}})^{3}({X_{\bar{y},g}X_{\bar{z}}-X_{\bar{z},g}X_{\bar{y}}\over
X^{2}_{\bar{y}}})_{\bar{z}}=({X_{\bar{y},g}X_{\bar{z}}-X_{\bar{z},g}X_{\bar{y}}\over
X^{2}_{\bar{z}}})_{\bar{y}}$

From the last equality it follows immediately

$X_{\bar{y},g}X_{\bar{z}}-X_{\bar{z},g}X_{\bar{y}}=F(X,g)\to 1$

where $F$ arbitrary function of its two arguments. At last by substitution
$X\to Y(X,g)$ it is possible to equate $F$ to unity and we come back to
equation considered in 4 section. Solution of Plebansky equation id the
following one

$v_{\bar{y}}=\exp({R+\Delta\over
2}=({\exp\Delta\over-\Delta_{\theta}})^{{1\over 2}}={g_{\theta}}^{-{1\over
2}}X_{\bar{y}},\quad v_{\bar{z}}={g_{\theta}}^{-{1\over 2}}X_{\bar{z}},\quad
v={g_{\theta}}^{-{1\over 2}}X$

###  6.2 Solution do not depend on one bar coordinate

Let us assume that solution of equation (19) does not depend on $\bar{z}$
coordinate. Then equation for $\Theta$ function takes the form

$e^{-\Delta}(e^{-\Delta}\Theta_{\Delta})_{\Delta}\Theta_{D,\bar{y}}-e^{-\Delta}\Theta_{\Delta,D}e^{-\Delta}\Theta_{\Delta,\bar{y}})=0$

This is exactly three dimensional subclass of complex Monge-Ampher (four
dimensional equation), solution of which in implicit form is known [3].

It may be expressed via (in terms of) $\psi(\Delta,d,\bar{y})$ function which
is solution of the equation

$e^{\Delta}+F_{\psi}(D,\psi)+\bar{F}_{\psi}(\bar{y},\psi)=0$ (22)

In the last equation for $\psi(\Delta,D,\bar{y})$ two arbitrary functions
$F(D,\psi),\bar{F}(\bar{y},\psi)$ assumed to be known. Of course solution of
this equation (in general case) may be obtained only in implicit form.

After this solution of the equation for $\Theta(D,\bar{y},\Delta)$ function is
resolves as follows

$\Theta_{\Delta}=e^{\Delta}\psi,\quad\Delta_{\Theta}=e^{-\Delta}\psi^{-1},\quad\Theta_{\bar{y}}=-\bar{F}_{\bar{y}},\quad\Theta_{D}=F_{D}$

And this is a general solution of this equation.

The solution of the initial equation of Plebansky in coordinates
$z,y,\bar{z},\bar{y}$ is given by equations

$v_{\bar{y}}=e^{\Delta}(\psi)^{{1\over 2}},\quad v_{\bar{z}}=(\psi)^{{1\over
2}},\quad yz+F_{D}(D,\psi)=0$ (23)

We present below direct checking of the formulae above. Equation of Plebansky
in its basic form after substitution (23) looks as

${1\over 2}(e^{\Delta}_{y}\psi_{z}-e^{\Delta}_{z}\psi_{y})=-{1\over
2}F_{\psi,D}{y\psi_{y}+z\psi_{z}\over yz}=1$

In the last step of calculation necessary use last formulae of (23)

$y+F_{DD}{1\over z}+F_{D,\psi}\psi_{z}=0,\quad z-F_{DD}{1\over
y}+F_{D,\psi}\psi_{y}=0,\quad D=\ln z-\ln y$

####  6.2.1 Solution of the example from a section 4

Let us sick solution of (19), (21) in a form (of course such form is a direct
consequence of calculations of the 4 section)

$\Theta=(e^{-\Delta}-r(D,\bar{y},\bar{z}))^{-1}$

All necessary derivatives are the following ones

$\Theta_{\Delta}=e^{-{1\over
2}\Delta}\Theta^{2},\quad\Theta_{D}=r_{D}\Theta^{2},\quad\Theta_{\bar{y}}=r_{\bar{y}}\Theta^{2},\quad\Theta_{\bar{z}}=r_{\bar{z}}\Theta^{2}$

After calculation of second order derivatives and substitution into (21) we
obtain equation for $r$ function

$r^{3}({1\over r})_{D,\bar{y}}=r_{D,\bar{z}}$

This equation is exactly integrable (author have not met it in literature
before) with solution

$r=-{X_{\bar{z}}\over X_{\bar{y}}}$

where $X(D,\bar{y},\bar{z})$ is exactly the function from the section 4.
Indeed $r_{D}=-{X_{D,\bar{z}}X_{\bar{y}}-X_{D,\bar{y}}X_{\bar{z}}\over
X^{2}_{\bar{y}}}={1\over X^{2}_{\bar{y}}}$.

####  6.2.2 Once more possible particular solution

Let us sick solution of (19), (21) in a form

$\Theta=P(\bar{y},\bar{z},\Delta)+Q(d,\Delta)$

After substitution into corresponding equations we come to the linear equation
for determining of the $P_{\Delta}$ function

$(e^{-{\Delta\over 2}}{\frac{\partial}{\partial\bar{y}}}-e^{{\Delta\over
2}}{\frac{\partial}{\partial\bar{z}}})P_{\Delta}=0$

with the obvious solution

$P=\int^{\Delta}d\delta p(e^{{-\delta\over 2}}\bar{z}+e^{{\delta\over
2}}\bar{y},\delta)$

But really for solution of Plebansky equation it is necessary only derivative
of $P$ function with respect to the $\Delta$ argument. Indeed

$v_{\bar{y}}=({e^{\Delta}\over\Delta_{\theta}})^{{1\over 2}},\quad
v_{\bar{z}}=({e^{-\Delta}\over\Delta_{\theta}})^{{1\over
2}},\quad\Delta_{\theta}={1\over\Theta_{\Delta}}={1\over
P_{\Delta}+Q_{\Delta}}$

and at last connection between $\theta$ and $\Delta$ is given by relation
$yz={\Delta_{D}\over\Delta_{\theta}}=-\Theta_{D}=-Q_{D}$.

##  7 Discrete transformation

Let us rewrite main equations of preliminary section once more

$v_{\bar{y}}=v_{y,\bar{y}}\theta_{z}-v_{z,\bar{y}}\theta_{y}+zv_{z\bar{y}}+yv_{y,\bar{y}},\quad
v_{\bar{z}}=v_{y,\bar{z}}\theta_{z}-v_{z,\bar{z}}\theta_{y}+zv_{z\bar{z}}+yv_{y,\bar{z}}$
(24)

and pay attention that in connection with (8) these equations may be rewritten
as

$v_{\bar{y}}=\tilde{\theta}_{\bar{y}}+zv_{z,\bar{y}}+yv_{y,\bar{y}},\quad
v_{\bar{z}}=\tilde{\theta}_{\bar{z}}+zv_{z,\bar{z}}+yv_{y,\bar{z}}$

The last equation may resolved with respect $\tilde{\theta}$ function
$\tilde{\theta}=v-zv_{z}-yv_{y}$. And thus after substitution this expression
for $\tilde{v}$ function we obtain linear system of equations for $\tilde{v}$
in a form

$\tilde{v}_{\bar{y}}=(-(Ov_{z})\frac{\partial}{\partial
y}+(Ov_{y})\frac{\partial}{\partial
z}+O)\tilde{v}_{\bar{y}},\quad\tilde{v}_{\bar{y}}=(-(Ov_{z})\frac{\partial}{\partial
y}+(Ov_{y})\frac{\partial}{\partial z}+O)\tilde{v}_{\bar{z}}$ (25)

where operator $O\equiv(y\frac{\partial}{\partial y}+z\frac{\partial}{\partial
z})$. Equations (25) is linear system of equations for determining
$\tilde{v}_{\bar{y}},\tilde{v}_{\bar{z}}$ functions by known $v$ solution of
Plebansky equation.

At this moment it is unknown to the author the systematic way for resolving
(25).

##  8 Outlook

The main result of the present paper is a new approach to the problem of
Plebansky equation. This approach allows rewrite Plebansky equation in
coordinates involving symmetry function and find series solutions of this
equation. Solutions obtained in such way depends on two arbitrary functions
each of two variables. Thus this is not a general solution of Plebansky
equation which must depend on two arbitrary function each of three variables.
But investigation of symmetry equation show that general solution of Plebansky
equation it is not possible present in analytic form. Situation exactly the
same as for instance in the case of famous nonlinear one dimensional
Schredinger equation, where it is possible to find infinite series solution of
soliton like type but not general one.

Constructed in the present paper solutions are connected with ordinary
differential equation of the second order $X_{z,z}=F(X,z)$ and it arise very
interesting problem to understand what relation has this equation to group of
symmetry of Plebansky equation responsible for its integrable properties.

##  9 Appendix

In this Appendix we consider simplest example from section $(6.0.1)$ chose
arbitrary functions in such simple form that all calculation are possible to
do in explicit form. Let us chose $F={(d+\psi)^{2}\over
2},\bar{F}={(\bar{y}+\psi)^{2}\over 2}$. The main equation allow determine
$\psi$ and all other necessary variables in explicit form in usual for
Plebansky equation variables

$\psi=-{e^{\Delta}+\bar{y}+d\over 2},\quad
yz={\Delta_{d}\over\Delta_{\theta}}=-\Theta_{d}=-(d+\psi)={e^{\Delta}+\bar{y}-d\over
2},$ $e^{\Delta}=2yz+D-\bar{y},\quad\psi=-(yz+d),\quad D=\ln{z\over y}$

In connection with results of section $(6.0.1)$ solution of basic Plebansky
equation looks as

$v=[(2yz+\ln{z\over y})\bar{y}-{(\bar{y})^{2}\over 2}+\bar{z}](yz+\ln{z\over
y})^{1\over 2}$

## References

  * [1] J.F.Plebanski Nucl. Phys. B373 (1992) 214-232.,of J.Math.Phys 16,(1975), 2395.
  * [2] L.V.Ovsjanikov Group Analysis of differential equations., Acad.Press New-York (1992).
  * [3] Fairlie D.B. and A. N. Leznov J. Phys. A Volume 33, Number 25, 30 June 2000 (4657-4661


for all $q\geq 1$. ∎

###  4.3 Consequences for the Wasserstein convergence of persistence diagrams

As a simple consequence of 4.2 and Proposition 4.12, we obtain that for $i<m$,
the $p$-Wasserstein convergence of $(\mathrm{dgm}_{i}({\mathsf{A}}_{n}))$ to
$\mathrm{dgm}_{i}({\mathsf{M}})$ holds if and only if $p>m$, as well as
precise asymptotics for the total persistence of
$(\mathrm{dgm}_{i}({\mathsf{A}}_{n}))$.

######  Corollary 4.17.

Let $p\geq 1$ and let $0\leq i<d$ be an integer. Under the same assumptions as
in Proposition 4.12, the following holds:

  * •

If $p>m$, then
${\mathbb{E}}[\mathrm{OT}^{p}_{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))]\to
0$ as $n\to\infty$.

  * •

If $p=m$,
${\mathbb{E}}[\mathrm{OT}^{p}_{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))]\to\mathrm{Pers}_{p}(\mu_{\infty,i,m})\mathrm{Vol}({\mathsf{M}})$
as $n\to\infty$, where $\mathrm{Vol}({\mathsf{M}})$ is the volume of
${\mathsf{M}}$.

  * •

If $p<m$ and $i<m$, then
${\mathbb{E}}[\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))]\to+\infty$
as $n\to\infty$.

Furthermore, for all $\alpha>0$,
$\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}))$ is equal to

$\displaystyle\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))+n^{1-\frac{\alpha}{m}}\mathrm{Pers}_{\alpha}(\mu_{\infty,i,m})\int_{\mathsf{M}}f(x)^{1-\frac{\alpha}{m}}\mathrm{d}x+o_{L^{1}}(n^{1-\frac{\alpha}{m}})+O_{L^{1}}\Big{(}\left(\frac{\log
n}{n}\right)^{\frac{1}{m}}\Big{)}.$

As noted earlier, both $\mathrm{dgm}_{i}({\mathsf{A}}_{n})$ and
$\mathrm{dgm}_{i}({\mathsf{M}})$ are trivial if $i\geq d$.

This corollary, whose proof can be found in the Appendix, gives a precise
answer to the questions raised in the introduction. First, when
${\mathsf{A}}_{n}$ is a random subset of a $m$-dimensional generic manifold in
${\mathbb{R}}^{d}$, the $\alpha$-total persistence of
$\mathrm{dgm}_{i}({\mathsf{A}}_{n})$ is not only bounded for $\alpha>d$ (as
was shown by Cohen-Steiner & al. [CSEHM10]), but for all $\alpha\geq m$.
Moreover, the sequence $\mathrm{dgm}_{i}({\mathsf{A}}_{n})$ converges for the
$\mathrm{OT}_{p}$ distance if $p>m$. A curious phenomenon can be observed in
the case $p=m$: the sequence does not converge to
$\mathrm{dgm}_{i}({\mathsf{M}})$ as one would expect, but its distance to the
power $p$ to $\mathrm{dgm}_{i}({\mathsf{M}})$ converges to some constant–in
that case, the cost to the power $p$ of matching all the points in Region (1)
to the diagonal $\partial\Omega$ neither converges to $0$ nor diverges, but is
asymptotically equal to this constant.

Using these asymptotics for the total persistence, we obtain regularity
guarantees for a large family of feature maps, called linear feature maps,
which includes feature maps introduced in [CFL+14, AEK+17, RHBK15, CWRW15,
KHF16, STR+18, CCI+20]. Let $(V,\|\cdot\|)$ be a normed vector space, and let
$\phi:\Omega\to V$ be a continuous bounded map. For $\alpha\geq 0$, the linear
feature map $\Phi_{\alpha}$ associated to $\phi$ and defined on the space
$\mathcal{D}_{f}$ of PDs having a finite number of points is defined for all
$a\in\mathcal{D}_{f}$ by $\Phi_{\alpha}(a)=\sum_{u\in
a}\mathrm{pers}(u)^{\alpha}\phi(u)\in V$.

######  Corollary 4.18.

Let $\alpha\geq 1$ and let $0\leq i<d$ be an integer. Under the same
assumptions as in Proposition 4.12, it holds that
$\Phi_{\alpha}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}))$ converges in probability
to $\Phi_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))$ whenever $\alpha>m$.

The proof is given in the Appendix. Remark that other weighting schemes are
possible. For instance, [KFH18] argued for using linear feature maps of the
form

$\Phi_{\alpha}(a)=\sum_{u\in a}\arctan(\mathrm{pers}(u)^{\alpha})\phi(u).$

Similar results would hold for such feature maps, as the map
$u\mapsto\arctan(\mathrm{pers}(u)^{\alpha})\phi(u)/\mathrm{pers}(u)^{\alpha}$
is continuous and bounded whenever $\phi$ is.

##  5 Numerical experiments

Figure 4: Left: the Čech PD $\mathrm{dgm}_{1}({\mathsf{A}}_{n})$ of a sample
of $n=10^{4}$ points sampled on a generic torus, with points in Regions (1),
(2) and (3) highlighted in different colors. Right: the persistence images of
$\mathrm{dgm}_{1}({\mathsf{A}}_{n})$ with weight $\mathrm{pers}^{p}$ for
different values of $p$.

We illustrate our results with synthetic experiments. We create a generic
submanifold of dimension $m$ by applying a random diffeomorphism $\Psi$ to a
given $m$-dimensional submanifold ${\mathsf{M}}_{0}$ (e.g. a torus). We then
draw a sample of $n$ i.i.d. observations sampled according to the pushforward
$P$ of the uniform distribution on ${\mathsf{M}}$ by $\Psi$.

Figure 5: Plot in log-log scale of
$\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))$ as a function of
$n$ for points sampled on a circle, $i=0$ (left), points sampled on a torus,
$i=0$ (center), points sampled on a torus, $i=1$ (right). Dashed lines have
slopes equal to $1-p/m$.

Continuity of feature maps. As a first experiment, we test the continuity of a
feature map, the persistence image [AEK+17]. In Figure 4, we plot the
persistence image of $\mathrm{dgm}_{1}({\mathsf{A}}_{n})$ where
${\mathsf{A}}_{n}$ is a sample of $n=10^{4}$ points on a generic torus. We
observe that the map is discontinuous for $p<2$: the two points with large
persistence corresponding to the PD of the underlying torus are nonapparent in
the image. For $p>2$, the two points are apparent, and the contribution of
points with small persistence (close to the lower edge) has vanished. In the
limit situation $p=2$, we see the contribution of both points with large and
small persistence. This phenomenon suggests the following heuristics: when in
presence of multiple datasets on $m$-dimensional objects whose global
geometries need to be distinguished, feature maps with weights
$\mathrm{pers}^{p}$ with $p>m$ should be used; when the relevant information
is the underlying density of the datasets, the choice $p<m$ should be
preferred.

Convergence of total persistences. We verify the rate of convergence of the
total persistence predicted by 4.2. For values of $n$ ranging from $10^{2}$ to
$10^{4}$, we compute
$\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))$ in three
scenarios: points sampled on a circle for $i=0$, and points sampled on a torus
for $i=0$ and $i=1$. The correct rates of convergence are observed on a log-
log plot, see Figure 5. For $i=1$, we remark that the asymptotic regime starts
at larger values of $n$, above $n=10^{3}$.

Convergence of $\mu_{n,i}$. We sample $n$ points on a torus by uniformly
sampling the two angles $(\theta,\phi)$ parametrizing the torus. We obtain a
(nonnuniform) probability measure, having density $f$. We then compute, for
various values of $n$, the measure $\mu_{n,1}$. The measure is approximated by
kernel density estimation (see Figure 6). We approximate in a similar manner
the measure $\mu_{\infty,1,2}$ by sampling $n=10^{5}$ points on a square. We
then apply the change of variable formula (21) to compute the theoretical
limit $\mu_{f,1}$. The distance $\mathrm{OT}_{2}(\mu_{n,1},\mu_{f,1})$ is then
computed by approximating the measures on a grid: the distance converges to
$0$ as predicted by 4.2. See also Figure 6.

Figure 6: Left: Heatmap of $\mu_{n,1}$ for $n=5\cdot 10^{4}$ points sampled on
the torus with density $f$. Center: Heatmap of $\mu_{f,1}$. Right:
$\mathrm{OT}_{2}$ distance between $\mu_{n,1}$ and $\mu_{f,1}$ (normalized by
$\mathrm{OT}_{2}(0,\mu_{f,i})$) for $n$ ranging from $10^{2}$ to $10^{5}$.

##  6 Conclusion

Under the manifold hypothesis, we have greatly refined earlier work regarding
the persistent homology of subsamples of compact sets, with especially strong
results when the sampling is either random or well-behaved. In particular, we
have precisely described the PDs of such samplings, and provided new
convergence guarantees w.r.t. the $p$-Wassertein distances, as well as
detailed asymptotics for their total $\alpha$-persistence. The main
limitations of our work were the assumptions that the data is sampled from a
submanifold, and without any noise. Relaxing those assumptions, as well as
establishing similar guarantees for Vietoris-Rips complexes, could be the
subject of future research. We also plan on exploring the consequences of our
findings regarding the persistent homology dimension [AAF+20, SSDE20, BLGS21,
Sch21] of submanifolds.

##  Appendix A Proofs of additional lemmas in Section 3

The goal of this section is to prove the following lemma, used in the proof of
3.7.

######  Lemma A.1.

Let ${\mathsf{A}}\subset{\mathbb{R}}^{d}$ be a finite set, and let
$a<b\in{\mathbb{R}}\cup\\{-\infty,+\infty\\}$ and $i\geq 0$. Let

$\bigcup_{j=1}^{L}C_{j}=\mathrm{Crit}({\mathsf{A}})\cap
d_{\mathsf{A}}^{-1}[a,b]$

be a covering of the set of critical points of $d_{\mathsf{A}}$ whose critical
value belongs to $[a,b]$,777When $a=-\infty$, we commit a minor abuse of
notation by writing $[a,b]$ rather than $(a,b]$, and similarly when
$b=+\infty$. and let

$N_{j}:=\\#\left(\bigcup_{z\in C_{j}}\sigma_{\mathsf{A}}(z)\right)$

be the cardinality of the union of the projections of the critical points of
$C_{j}$. Then the number of points in $\mathrm{dgm}_{i}({\mathsf{A}})$ such
that at least one of their coordinates belongs to $[a,b]$ is upper-bounded by
$\sum_{j=1}^{L}\binom{N_{j}}{i+2}$, hence by $\sum_{j=1}^{L}N_{j}^{i+2}$.

We first require an intermediate result. Following [BE17, Definition 4.2], we
say that a finite set ${\mathsf{B}}\subset{\mathbb{R}}^{d}$ is in general
position if for every ${\mathsf{S}}\subset{\mathsf{B}}$ of at most $d+1$
points, ${\mathsf{S}}$ is affinely independent in ${\mathbb{R}}^{d}$ and no
point of ${\mathsf{B}}\backslash{\mathsf{S}}$ lies on the smallest
circumsphere of ${\mathsf{S}}$.

######  Lemma A.2.

Let ${\mathsf{A}}\subset{\mathbb{R}}^{d}$ be a finite set. For any
$\varepsilon>0$, there exists a finite set
${\mathsf{B}}\subset{\mathbb{R}}^{d}$ such that:

  * •

${\mathsf{B}}$ is in general position,

  * •

the distance function $d_{\mathsf{B}}$ is topologically Morse,

  * •

there exists a bijection $\phi:{\mathsf{B}}\rightarrow{\mathsf{A}}$ with
$d(\phi(b),b)\leq\varepsilon$ for all $b\in{\mathsf{B}}$, and

  * •

there exists a map
$\psi:\mathrm{Crit}({\mathsf{B}})\rightarrow\mathrm{Crit}({\mathsf{A}})$ such
that $|d_{\mathsf{B}}(z)-d_{\mathsf{A}}(\psi(z))|\leq\varepsilon$ and
$\phi(\sigma_{\mathsf{B}}(z))\subset\sigma_{\mathsf{A}}(\psi(z))$ for all
$z\in\mathrm{Crit}({\mathsf{B}})$.

###### Proof.

Let $N$ be the cardinality of ${\mathsf{A}}=\\{p_{1},\ldots,p_{N}\\}$. As
being in general position is a dense property, we can find for any $\delta>0$
a set ${\mathsf{B}}:=\\{q_{1},\ldots,q_{N}\\}\subset{\mathbb{R}}^{d}$ such
that $d(q_{i},p_{i})\leq\delta$ and such that the set ${\mathsf{B}}$ is in
general position in ${\mathbb{R}}^{d}$. If $\delta<\varepsilon$ is small
enough, this defines a bijection
$\phi:{\mathsf{B}}\rightarrow{\mathsf{A}},q_{i}\mapsto p_{i}$ that satisfies
the second condition of the statement. As ${\mathsf{B}}$ is in general
position, the distance function $d_{\mathsf{B}}$ is then topologically Morse
(see [GR97]).

Moreover, Theorem 3.4 from [CCSL06], combined with the lower semicontinuity of
the norm of the generalized gradient $\nabla d_{\mathsf{A}}$ and the
compactness of ${\mathsf{A}}$, implies that for any $\rho>0$, each critical
point of ${\mathsf{B}}$ must be at distance less than $\rho$ from one of the
critical points of ${\mathsf{A}}$ as long as $\delta>0$ is small enough. This
defines a map
$\psi:\mathrm{Crit}({\mathsf{B}})\rightarrow\mathrm{Crit}({\mathsf{A}})$. If
$\delta+\rho<\varepsilon$, then
$|d_{\mathsf{B}}(z)-d_{\mathsf{A}}(\psi(z))|\leq|d_{\mathsf{B}}(z)-d_{\mathsf{A}}(z)+d_{\mathsf{A}}(z)-d_{\mathsf{A}}(\psi(z))|\leq\delta+\rho<\varepsilon$
and the first condition on $\psi$ is satisfied.

Furthermore, for each of the finitely many critical point
$z^{\prime}\in\mathrm{Crit}({\mathsf{A}})$, there exists
$\lambda(z^{\prime})>0$ such that
$B(z^{\prime},d_{\mathsf{A}}(z^{\prime})+\lambda(z^{\prime}))\cap{\mathsf{A}}=\sigma_{\mathsf{A}}(z^{\prime})$.
Assume now that $\delta>0$ is small enough that
$\delta<\min_{z^{\prime}\in\mathrm{Crit}({\mathsf{A}})}\lambda(z^{\prime})/4$
and that each critical point $z$ of ${\mathsf{B}}$ is at distance less than
$\min_{z^{\prime}\in\mathrm{Crit}({\mathsf{A}})}\lambda(z^{\prime})/4$ from
$\psi(z)\in\mathrm{Crit}({\mathsf{A}})$. Then such a $z$ is necessarily such
that $\sigma_{\mathsf{B}}(z)\subset\phi^{-1}(\sigma_{\mathsf{A}}(\psi(z)))$.
Hence $\psi$ is as required for $\delta$ small enough, and the lemma is
proved. ∎

###### Proof of Lemma A.1.

Lemma A.2 states that for any $\delta>0$, there exists a point cloud
${\mathsf{B}}$ in general position, a bijection
$\phi:{\mathsf{B}}\rightarrow{\mathsf{A}}$ with $d(\phi(b),b)\leq\delta$ for
all $b\in{\mathsf{B}}$, and a map
$\psi:\mathrm{Crit}({\mathsf{B}})\rightarrow\mathrm{Crit}({\mathsf{A}})$ such
that $|d_{\mathsf{B}}(z)-d_{\mathsf{A}}(\psi(z))|\leq\delta$ and
$\phi(\sigma_{\mathsf{B}}(z))\subset\sigma_{\mathsf{A}}(\psi(z))$ for all
$z\in\mathrm{Crit}({\mathsf{B}})$ and such that $d_{\mathsf{B}}$ is
topologically Morse. Let $\rho>0$ be such that all critical values of
$d_{\mathsf{A}}$ in $[a-\rho,b+\rho]$ belong to $[a,b]$, and let
${\mathsf{B}}$ be as described for some $\delta<\rho/2$. Since the distance
function $d_{\mathsf{B}}$ is topologically Morse, and as explained in the
proof of Proposition 3.6, the critical points
$z\in\mathrm{Crit}({\mathsf{B}})$ are topological critical points and are in
bijection (via their critical values) with the non-zero coordinates of the
points in (the union over all degrees of) the Čech persistence diagrams of
${\mathsf{B}}$. In particular, the number of points in
$\mathrm{dgm}_{i}({\mathsf{B}})$ such that at least one of their coordinates
belongs to $[a-\rho/2,b+\rho/2]$ is upper bounded by the number of points
$z\in\mathrm{Crit}({\mathsf{B}})$ of topological Morse index $i$ or $i+1$ such
that $d_{\mathsf{B}}(z)\in[a-\rho/2,b+\rho/2]$. For any such $z$, the point
$\psi(z)\in\mathrm{Crit}({\mathsf{A}})$ is such that
$d_{{\mathsf{A}}_{n}}(\psi(z))\in[a,b]$, and as such belongs to $C_{j}$ for
some $j\in\\{1,\ldots,L\\}$. Furthermore,
$\phi(\sigma_{\mathsf{B}}(z))\subset\sigma_{{\mathsf{A}}}(\psi(z))\subset\bigcup_{z^{\prime}\in
C_{j}}\sigma_{{\mathsf{A}}}(z^{\prime})$. As ${\mathsf{B}}$ is in general
position, each critical point of ${\mathsf{B}}$ of index $i$ is uniquely
identified by its $i+1$ projections on ${\mathsf{B}}$ (respectively $i+2$
projections for points of index $i+1$), and this bounds the total number of
$z\in\mathrm{Crit}({\mathsf{B}})$ such that
$d_{\mathsf{B}}(z)\in[a-\rho/2,b+\rho/2]$ by
$\sum_{j=1}^{L}\binom{N_{j}}{i+2}$.

We have shown that $\mathrm{dgm}_{i}({\mathsf{B}})$ does not have more than
$\sum_{j=1}^{n}\binom{N_{j}}{i+2}$ points with at least one of their
coordinates greater $\in[a-\rho/2,b+\rho/2]$. By the Bottleneck Stability
Theorem, there is a surjection from this set to the set of points in
$\mathrm{dgm}_{i}({\mathsf{A}})$ with at least one of their coordinates in
$[a,b]$ as soon as $\delta>0$ is small enough. Hence there are at most
$\sum_{j=1}^{n}\binom{N_{j}}{i+2}$ points in $\mathrm{dgm}_{i}({\mathsf{A}})$
with at least one of their coordinates in $[a,b]$, as desired. ∎

##  Appendix B Proofs of additional lemmas in Section 4

In this section, we give additional lemmas needed to prove 4.2. First, we
state two simple lemmas on properties of probability measures and of samples
on ${\mathsf{M}}$.

######  Lemma B.1.

Let ${\mathsf{M}}$ be a compact submanifold with positive reach. Let $P$ be a
probability measure having a density $f$ on ${\mathsf{M}}$ satisfying
$f_{\min}\leq f\leq f_{\max}$ for two strictly positive constants $f_{\min}$,
$f_{\max}$. There exist constants $c=c(m),C=C(m)$ such that for all $0\leq
r\leq\tau({\mathsf{M}})/4$

$cf_{\min}r^{m}\leq P(\overline{B}(x,r))\leq Cf_{\max}r^{m}.$ (39)

Let ${\mathsf{A}}_{n}$ be a sample of $n$ i.i.d. observations of law $P$.
Then, there exists $C_{0}=C_{0}({\mathsf{M}})$ depending on ${\mathsf{M}}$
such that for all $x\in{\mathsf{M}}$ and all $r>0$,

${\mathbb{P}}(d(x,{\mathsf{A}}_{n})\geq r)\leq\exp(-nC_{0}f_{\min}r^{m}).$
(40)

###### Proof.

For the first statement, see [AL18, Proposition 31]. Let us prove the second
one. Remark that the probability is zero for $r>\mathrm{diam}({\mathsf{M}})$.
Hence, we can assume that $r\leq\mathrm{diam}({\mathsf{M}})$. When
$r\leq\tau({\mathsf{M}})/4$, it holds that

$\displaystyle{\mathbb{P}}(d(x,{\mathsf{A}}_{n})\geq
r)=(1-P(\overline{B}(x,r)))^{n}\leq\exp(-nc_{m}f_{\min}r^{m}).$

When $\tau({\mathsf{M}})/4\leq r\leq\mathrm{diam}({\mathsf{M}})$, we write

$\displaystyle{\mathbb{P}}(d(x,{\mathsf{A}}_{n})\geq r)$
$\displaystyle\leq{\mathbb{P}}(d(x,{\mathsf{A}}_{n})\geq\tau({\mathsf{M}})/4)\leq\exp(-ncf_{\min}(\tau({\mathsf{M}})/4)^{m})$
$\displaystyle\leq\exp(-ncf_{\min}\frac{(\tau({\mathsf{M}})/4)^{m}}{\mathrm{diam}({\mathsf{M}})^{m}}r^{m}).$

Hence, the result holds with
$C_{0}=c\min\left(1,\frac{(\tau({\mathsf{M}})/4)^{m}}{\mathrm{diam}({\mathsf{M}})^{m}}\right)$.
∎

See 4.6

###### Proof.

The bound ${\mathbb{P}}(d_{H}({\mathsf{A}}_{n},{\mathsf{M}})>r)$ is given in
[Aam17, Lemma III.23]. Furthermore, [Aam17, Lemma III.23] also states that for
any $q>0$, there exists $C(q)$ depending on $f_{\min}$ and $m$ such that, with
probability at least $1-n^{-q/m}$, $d_{H}({\mathsf{A}}_{n},{\mathsf{M}})\leq
C(q)\left(\frac{\ln n}{n}\right)^{1/m}$. In particular, we obtain that

$\displaystyle{\mathbb{E}}[d_{H}({\mathsf{A}}_{n},{\mathsf{M}})^{q}]\leq
C(q)^{q}\left(\frac{\ln
n}{n}\right)^{q/m}+\mathrm{diam}({\mathsf{M}})^{q}n^{-q/m},$

proving the second claim of the lemma. ∎

Next, we state a simple geometric lemma.

######  Lemma B.2.

Let ${\mathsf{A}}\subset{\mathsf{M}}$, $z\in\mathrm{Crit}({\mathsf{A}})$ and
$x\in\sigma_{\mathsf{A}}(z)$. Then
$\|\pi_{x}^{\bot}(z-x)\|\leq\|z-x\|^{2}/\tau({\mathsf{M}})$. Furthermore, if
$d_{\mathsf{A}}(z)\leq\tau(M)/\sqrt{2}$, then
$\|\pi_{x}(z-x)\|\geq\|z-x\|/\sqrt{2}$.

###### Proof.

The point $z$ can be written as a convex combination
$z=\sum_{k}\lambda_{k}y_{k}$ where the points $y_{k}$ are in
${\mathsf{A}}\subset{\mathsf{M}}$. Then, using [Fed59, Theorem 4.18],

$\displaystyle\|\pi_{x}^{\bot}(z-x)\|$
$\displaystyle\leq\sum_{k}\lambda_{k}\|\pi_{x}^{\bot}(y_{k}-x)\|\leq\sum_{k}\lambda_{k}\frac{\|y_{k}-x\|^{2}}{2\tau({\mathsf{M}})}$
$\displaystyle=\sum_{k}\lambda_{k}\frac{\|y_{k}-z\|^{2}+\|z-x\|^{2}+2\langle
y_{k}-z,z-x\rangle}{2\tau({\mathsf{M}})}=\frac{\|z-x\|^{2}}{\tau({\mathsf{M}})},$

as $\|y_{k}-z\|^{2}=\|x-z\|^{2}$ and $\sum_{k}\lambda_{k}(y_{k}-z)=0$. The
second inequality $\|\pi_{x}(z-x)\|\geq\|z-x\|/\sqrt{2}$ from the statement
follows from the first one through a direct computation. ∎

At last, the following property of the multinomial distribution was used in
the proof of 4.2.

######  Lemma B.3.

Let $K\geq 1$ and $L=(L_{1},\ldots,L_{K},L_{K+1})$ be a random multinomial
variable of parameters $n$ and $\alpha_{1},\ldots,\alpha_{K},\alpha_{K+1}$.
Then there exists $C=C(K,l)$ such that

${\mathbb{E}}[(\sum_{i=1}^{K}L_{i})^{l}|L_{i}\geq 1\ \forall i=1,\ldots,K]\leq
C_{l}(1+(n\sum_{i=1}^{K}\alpha_{i})^{l}).$

###### Proof.

Let $X_{1},\ldots,X_{n}$ be i.i.d. categorical variables of parameters
$\alpha_{1},\ldots,\alpha_{K},\alpha_{K+1}$. Define $L_{k}(p)=\sum_{r\leq
p}\mathbf{1}_{X_{r}=k}$; then $(L_{1}(n),\ldots,L_{K}(n),L_{K+1}(n))$ has the
same distribution as $(L_{1},\ldots,L_{K},L_{K+1})$, and we identify the two
in our notation. For a fixed $n$, and for any injective function
$\iota:\\{1,\ldots,k\\}\rightarrow\\{1,\ldots,n\\}$, consider the event

$E_{\iota}:=\\{L_{1}(\iota(1))=1,\ldots,L_{1}(\iota(K))=1\\},$

i.e. $\iota(i)$ is the first appearance of $i$ among the variables
$X_{1},\ldots,X_{n}$. Note that $E:=\\{L_{1},\ldots,L_{K}\geq
1\\}=\bigsqcup_{\iota}E_{\iota}$ where the sum is taken over all such
injective functions. Then

${\mathbb{E}}[(\sum_{i=1}^{K}L_{i})^{l}|E]=\sum_{\iota}{\mathbb{P}}(E_{\iota}|E){\mathbb{E}}[(\sum_{i=1}^{K}L_{i})^{l}|E_{\iota}].$

Fix a function $\iota$, and assume without loss of generality that
$\iota(1)<\iota(2)<\ldots<\iota(K)$. Conditioned by $A_{\iota}$, the variable
$Y_{i}:=\sum_{r=\iota(i)+1}^{\iota(i+1)-1}\mathbf{1}_{X_{r}\neq K+1}$ is a
binomial variable of parameters $\iota(i+1)-\iota(i)-2\leq n$ and
$\frac{\sum_{i=1}^{i}\alpha_{i}}{\alpha_{1}+\ldots+\alpha_{i}+\alpha_{K+1}}\leq\sum_{i=1}^{K}\alpha_{i}$.
Hence ${\mathbb{E}}[Y_{i}^{l}|E_{\iota}]\leq
C_{1}(l)(n\sum_{i=1}^{K}\alpha_{i})^{l}$ using classical bounds on the $l$-th
moment of a binomial variable, and we see that

$\displaystyle{\mathbb{E}}[(\sum_{i=1}^{K}L_{i})^{l}|E_{\iota}]$
$\displaystyle={\mathbb{E}}[(\sum_{r=1}^{n}\mathbf{1}_{X_{r}\neq
K+1})^{l}|E_{\iota}]={\mathbb{E}}[(K+\sum_{i=0}^{K}\sum_{r=\iota(i)+1}^{\iota(i+1)-1}\mathbf{1}_{X_{r}\neq
K+1})^{l}|E_{\iota}]$ $\displaystyle\leq
C_{2}(K,l)(K^{l}+\sum_{i=0}^{K}{\mathbb{E}}[Y_{i}^{l}|E_{\iota}])\leq
C_{3}(K,l)(1+(n\sum_{i=1}^{K}\alpha_{i})^{l})$

where we write $\iota(0)=0$ and $\iota(K+1)=n$ to simplify notation. ∎

##  Appendix C Proofs of 4.17 and 4.18

We restate and prove two corollaries from Section 4.

See 4.17

###### Proof.

Let $\varepsilon_{n}=d_{H}({\mathsf{A}}_{n},{\mathsf{M}})$. Consider the event
$E=\\{\varepsilon_{n}<\tau(M)/2\\}$. According to 2.3, on $E$, the number of
points of $\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n})$ is bounded by some
constant $N_{0}$ depending only on ${\mathsf{M}}$.

  * •

Consider the optimal matching $\gamma_{n}$ given by Proposition 4.12. On the
event $E$, this optimal matching sends all the points of
$\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n})$ to the diagonal $\partial\Omega$.
Thus, we have (when $E$ is satisfied)

$\begin{split}&|\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))-\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))|\\\
&\leq
N_{0}\max_{u\in\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n})}\|u-\gamma_{n}(u)\|_{\infty}^{p}+\\#(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))\max_{u\in\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n})}\|u-\gamma_{n}(u)\|_{\infty}^{p}.\end{split}$
(41)

Furthermore, note that

$\displaystyle\max_{u\in\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n})}\|u-\gamma_{n}(u)\|_{\infty}^{p}$
$\displaystyle\leq\max\left(\max_{(u_{1},u_{2})\in\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n})}|u_{2}-\gamma_{n}(u)_{2}|^{p},\varepsilon_{n}^{p}\right)$
$\displaystyle\leq\max_{(u_{1},u_{2})\in\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n})}|u_{2}-\gamma_{n}(u)_{2}|^{p}+\varepsilon_{n}^{p}.$

We take the expectation and apply the Cauchy-Schwarz inequality to obtain
(together with 4.2, Proposition 4.12 and Lemma 4.6):

$\displaystyle{\mathbb{E}}[\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))\mathbf{1}\\{E\\}]$
$\displaystyle\leq n^{1-p/m}\mathrm{Pers}_{p}(\mu_{f,i})+o(1).$

When $E$ is not realized, we crudely bound
$\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))$
by considering the matching sending every point to the diagonal. The cost of
this matching is bounded by $C_{\mathsf{M}}(1+n^{i+1})$, where $n^{i+1}$ is an
upper bound on the number of points in $\mathrm{dgm}_{i}({\mathsf{A}}_{n})$
and $C_{\mathsf{M}}$ is some constant depending on ${\mathsf{M}}$. Hence,

$\displaystyle{\mathbb{E}}[\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))\mathbf{1}\\{E^{c}\\}]\leq
C_{\mathsf{M}}(1+n^{i+1}){\mathbb{P}}(\varepsilon_{n}>\tau({\mathsf{M}})/2)=o(1),$

according to Lemma 4.6. This proves the first bullet point.

  * •

The second bullet point is proved likewise from (41). Indeed, in that case, it
holds that

${\mathbb{E}}[|\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))-\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))|\mathbf{1}\\{E\\}]=o(1).$
(42)

Also, using the same crude bound, we obtain that

${\mathbb{E}}[|\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))-\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))|\mathbf{1}\\{E^{c}\\}]=o(1).$
(43)

So far, we have proved that
$\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))=\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))+o_{L^{1}}(1)$.
According to 4.2, it holds that
$\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))=\mathrm{Pers}_{p}(\mu_{f,i})+o_{L^{1}}(1)$
for $p=m$, proving the second bullet point.

  * •

For the third bullet point, we use that on the event $E$,

$\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))\geq\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n})).$
(44)

The latter is equal to
$n^{1-p/m}\mathrm{Pers}_{p}(\mu_{f,i})+o_{L^{1}}(n^{1-p/m})$, with
$\mathrm{Pers}_{p}(\mu_{f,i})>0$: indeed, the support of the measure
$\mu_{f,i}$ is nontrivial for $i<m$, see [HST18, GTT19]. Hence,

$\displaystyle{\mathbb{E}}[\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))]$
$\displaystyle\geq{\mathbb{E}}[\mathrm{OT}_{p}^{p}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))\mathbf{1}\\{E\\}]$
$\displaystyle\geq{\mathbb{E}}[\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))]-{\mathbb{E}}[\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))\mathbf{1}\\{E^{c}\\}].$

The second term goes to zero (use the crude bound
$\mathrm{Pers}_{p}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))\leq
C_{\mathsf{M}}(1+n^{i+1})$), while the first one diverges. This proves the
third bullet point.

At last, we prove the formula for the asymptotic expansion of

$\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}))=\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(1)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n})).$
(45)

The first term is equal to
$n^{1-\alpha/m}\mathrm{Pers}_{\alpha}(\mu_{\infty,i,m})\int_{\mathsf{M}}f(x)^{1-\frac{\alpha}{m}}\mathrm{d}x+o_{L^{1}}(n^{1-\alpha/m})$
according to 4.2. Remark that for $u,v\in\Omega$,
$|\mathrm{pers}_{\alpha}(u)-\mathrm{pers}_{\alpha}(v)|\leq
2\alpha(\mathrm{pers}_{\alpha}(u)+\mathrm{pers}_{\alpha}(v))\|u-v\|_{\infty}$.
Hence, on the event $E=\\{\varepsilon_{n}<\tau({\mathsf{M}})/2\\}$,

$\displaystyle|\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))-\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))|$
$\displaystyle\leq
2\alpha\varepsilon_{n}(\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))).$

We crudely bound the persistence of a point in
$\mathrm{dgm}_{i}({\mathsf{A}}_{n})$ by $R({\mathsf{M}})$, the radius of the
smallest enclosing ball of ${\mathsf{M}}$, to obtain that

$\displaystyle|\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))-\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))|$
$\displaystyle\leq
2\alpha\varepsilon_{n}(R({\mathsf{M}})^{\alpha}+1)(\\#(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))+\\#(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))).$

We have already established that
$\\#(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))=O_{L^{2}}(1)$ (actually, it is
larger than $N_{0}$ with only exponentially small probability). Furthermore,
Proposition 4.12 states that
$\\#(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))=O_{L^{2}}(1)$. Lemma 4.6 also
states that $\varepsilon_{n}=O_{L^{2}}(((\log n)/n)^{1/m})$. Hence,

$\displaystyle|\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(2)}({\mathsf{A}}_{n}))+\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}^{(3)}({\mathsf{A}}_{n}))-\mathrm{Pers}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))|=O_{L^{1}}(((\log
n)/n)^{1/m}).$

This concludes the proof. ∎

See 4.18

###### Proof.

According to [DL21, Proposition 5.1], the feature map $\Phi_{\alpha}$ is
continuous with respect to the $\mathrm{OT}_{\alpha}$ distance. But we have
shown in 4.17 that
${\mathbb{E}}[\mathrm{OT}_{\alpha}^{\alpha}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))]\to
0$ as $n\to\infty$. In particular,
$\mathrm{OT}_{\alpha}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}),\mathrm{dgm}_{i}({\mathsf{M}}))$
converges in probability to $0$. By continuity, so does
$\|\Phi_{\alpha}(\mathrm{dgm}_{i}({\mathsf{A}}_{n}))-\Phi_{\alpha}(\mathrm{dgm}_{i}({\mathsf{M}}))\|$.
∎

## References

  * [AAF+20] Henry Adams, Manuchehr Aminian, Elin Farnell, Michael Kirby, Joshua Mirth, Rachel Neville, Chris Peterson, and Clayton Shonkwiler.  A fractal dimension for measures via persistent homology.  In Topological Data Analysis: The Abel Symposium 2018, pages 1–31. Springer, 2020.
  * [Aam17] Eddie Aamari.  Rates of convergence for geometric inference.  PhD thesis, 2017.
  * [ACC+22] Andrew Aukerman, Mathieu Carrière, Chao Chen, Kevin Gardner, Raúl Rabadán, and Rami Vanguri.  Persistent homology based characterization of the breast cancer immune microenvironment: A feasibility study.  Journal of Computational Geometry, 12(2):183–206, 2022.
  * [ACSD23] Charles Arnal, David Cohen-Steiner, and Vincent Divol.  Critical points of the distance function to a generic submanifold.  arXiv preprint arXiv:2312.13147, 2023.
  * [AEK+17] Henry Adams, Tegan Emerson, Michael Kirby, Rachel Neville, Chris Peterson, Patrick Shipman, Sofya Chepushtanova, Eric Hanson, Francis Motta, and Lori Ziegelmeier.  Persistence images: A stable vector representation of persistent homology.  Journal of Machine Learning Research, 18(8):1–35, 2017.
  * [AKF+22] Dominique Attali, Hana Dal Poz Kouřimská, Christopher Fillmore, Ishika Ghosh, André Lieutier, Elizabeth Stephenson, and Mathijs Wintraecken.  Tight Bounds for the Learning of Homotopy à la Niyogi, Smale, and Weinberger for Subsets of Euclidean Spaces and of Riemannian Manifolds.  https://arxiv.org/abs/2206.10485, 2022.
  * [AL18] Eddie Aamari and Clément Levrard.  Stability and minimax optimality of tangential Delaunay complexes for manifold reconstruction.  Discrete & Computational Geometry, 59:923–971, 2018.
  * [BA14] Omer Bobrowski and Robert Adler.  Distance functions, critical points, and the topology of random Cech complexes.  Homology, Homotopy & Applications, 16(2), 2014.
  * [BCR+22] Bradley CA Brown, Anthony L Caterini, Brendan Leigh Ross, Jesse C Cresswell, and Gabriel Loaiza-Ganem.  Verifying the union of manifolds hypothesis for image data.  arXiv preprint arXiv:2207.02862, 2022.
  * [BE17] Urich Bauer and Herbert Edelsbrunner.  The Morse theory of Cech and Delaunay complexes.  Transactions of the American Mathematical Society, 369(5):pp. 3741–3762, 2017.
  * [BH22] Magnus B Botnan and Christian Hirsch.  On the consistency and asymptotic normality of multiparameter persistent Betti numbers.  Journal of Applied and Computational Topology, pages 1–38, 2022\.
  * [BHO18] Mickaël Buchet, Yasuaki Hiraoka, and Ippei Obayashi.  Persistent homology and materials informatics.  Nanoinformatics, pages 75–95, 2018.
  * [BK18] Omer Bobrowski and Matthew Kahle.  Topology of random geometric complexes: a survey.  Journal of applied and Computational Topology, 1:331–364, 2018\.
  * [BKS17] Omer Bobrowski, Matthew Kahle, and Primoz Skraba.  Maximally persistent cycles in random geometric complexes.  The Annals of Applied Probability, pages 2032–2060, 2017.
  * [BLGS21] Tolga Birdal, Aaron Lou, Leonidas J Guibas, and Umut Simsekli.  Intrinsic dimension, persistent homology and generalization in neural networks.  Advances in Neural Information Processing Systems, 34:6776–6789, 2021.
  * [BW17] Omer Bobrowski and Shmuel Weinberger.  On the vanishing of homology in random Cech complexes.  Random Structures & Algorithms, 51(1):14–51, 2017.
  * [BW20] Peter Bubenik and Alexander Wagner.  Embeddings of persistence diagrams into Hilbert spaces.  Journal of Applied and Computational Topology, 4(3):339–351, 2020\.
  * [CCI+20] Mathieu Carrière, Frédéric Chazal, Yuichi Ike, Théo Lacombe, Martin Royer, and Yuhei Umeda.  Perslay: A neural network layer for persistence diagrams and new graph topological signatures.  In International Conference on Artificial Intelligence and Statistics, pages 2786–2796. PMLR, 2020.
  * [CCO17] Mathieu Carriere, Marco Cuturi, and Steve Oudot.  Sliced Wasserstein kernel for persistence diagrams.  In International conference on machine learning, pages 664–673. PMLR, 2017.
  * [CCSL06] Frédéric Chazal, David Cohen-Steiner, and André Lieutier.  A sampling theory for compact sets in Euclidean space.  Discrete & Computational Geometry, 41:461–479, 06 2006.
  * [CDSGO16] Frédéric Chazal, Vin De Silva, Marc Glisse, and Steve Oudot.  The structure and stability of persistence modules, volume 10.  Springer, 2016.
  * [CFL+14] Frédéric Chazal, Brittany Terese Fasy, Fabrizio Lecci, Alessandro Rinaldo, and Larry Wasserman.  Stochastic convergence of persistence landscapes and silhouettes.  In Proceedings of the thirtieth annual symposium on Computational geometry, pages 474–483, 2014.
  * [CPH21] Luigi Caputi, Anna Pidnebesna, and Jaroslav Hlinka.  Promises and pitfalls of topological data analysis for brain connectivity analysis.  NeuroImage, 238:118245, 2021.
  * [CR20] Mathieu Carrière and Raúl Rabadán.  Topological data analysis of single-cell hi-c contact maps.  In Topological Data Analysis: The Abel Symposium 2018, pages 147–162. Springer, 2020.
  * [CSEH05] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer.  Stability of persistence diagrams.  In Proceedings of the twenty-first annual symposium on Computational geometry, pages 263–271, 2005.
  * [CSEHM10] David Cohen-Steiner, Herbert Edelsbrunner, John Harer, and Yuriy Mileyko.  Lipschitz functions have $L_{p}$-stable persistence.  Foundations of Computational Mathematics, 10(2):127–139, Apr 2010\.
  * [CWRW15] Yen-Chi Chen, Daren Wang, Alessandro Rinaldo, and Larry Wasserman.  Statistical analysis of persistence intensity functions.  arXiv preprint arXiv:1510.02502, 2015.
  * [DC19] Vincent Divol and Frédéric Chazal.  The density of expected persistence diagrams and its kernel based estimation.  Journal of Computational Geometry, 10(2):127–153, 2019.
  * [Div21] Vincent Divol.  Minimax adaptive estimation in manifold inference.  Electronic Journal of Statistics, 15(2):5888–5932, 2021.
  * [DKR+18] Joshua A Dijksman, Lenka Kovalcinova, Jie Ren, Robert P Behringer, Miroslav Kramar, Konstantin Mischaikow, and Lou Kondic.  Characterizing granular networks using topological metrics.  Physical Review E, 97(4):042903, 2018.
  * [DL21] Vincent Divol and Théo Lacombe.  Understanding the topology and the geometry of the space of persistence diagrams via optimal partial transport.  Journal of Applied and Computational Topology, 5:1–53, 2021.
  * [DP19] Vincent Divol and Wolfgang Polonik.  On the choice of weight functions for linear representations of persistence diagrams.  Journal of Applied and Computational Topology, 3(3):249–283, 2019\.
  * [EH10] Herbert Edelsbrunner and John Harer.  Computational Topology: An Introduction.  01 2010.
  * [ELZ02] Edelsbrunner, Letscher, and Zomorodian.  Topological persistence and simplification.  Discrete & computational geometry, 28:511–533, 2002.
  * [Fed59] Herbert Federer.  Curvature measures.  Transactions of the American Mathematical Society, 93(3):418–491, 1959.
  * [FLR+14] Brittany Terese Fasy, Fabrizio Lecci, Alessandro Rinaldo, Larry Wasserman, Sivaraman Balakrishnan, and Aarti Singh.  Confidence sets for persistence diagrams.  The Annals of Statistics, 42(6):2301, 2014.
  * [GPPVW12] Christopher R Genovese, Marco Perone-Pacifico, Isabella Verdinelli, and Larry Wasserman.  Manifold estimation and singular deconvolution under Hausdorff loss.  The Annals of Statistics, 40(2):941–963, 2012.
  * [GR97] Vladimir Gershkovich and Hyam Rubinstein.  Morse theory for min-type functions.  The Asian Journal of Mathematics, 1, 01 1997.
  * [Gro93] Karsten Grove.  Critical point theory for distance functions.  Proceedings of Symposia in Pure Mathematics, pages 357–385, 1993\.
  * [GTT19] Akshay Goel, Khanh Duy Trinh, and Kenkichi Tsunoda.  Strong law of large numbers for Betti numbers in the thermodynamic regime.  Journal of Statistical Physics, 174(4):865–892, 2019.
  * [Hat00] Allen Hatcher.  Algebraic topology.  Cambridge Univ. Press, Cambridge, 2000.
  * [HKN19] Christoph D Hofer, Roland Kwitt, and Marc Niethammer.  Learning representations of persistence barcodes.  Journal of Machine Learning Research, 20(126):1–45, 2019.
  * [HL23] Olympio Hacquard and Vadim Lebovici.  Euler characteristic tools for topological data analysis.  arXiv preprint arXiv:2303.14040, 2023.
  * [HST18] Yasuaki Hiraoka, Tomoyuki Shirai, and Khanh Duy Trinh.  Limit theorems for persistence diagrams.  The Annals of Applied Probability, 28(5):2740–2780, 2018.
  * [KFH18] Genki Kusano, Kenji Fukumizu, and Yasuaki Hiraoka.  Kernel method for persistence diagrams via kernel embedding and weight factor.  Journal of Machine Learning Research, 18(189):1–41, 2018.
  * [KHF16] Genki Kusano, Yasuaki Hiraoka, and Kenji Fukumizu.  Persistence weighted gaussian kernel for topological data analysis.  In International conference on machine learning, pages 2004–2013. PMLR, 2016.
  * [KKP+16] L Kondic, M Kramár, Luis A Pugnaloni, C Manuel Carlevaro, and K Mischaikow.  Structure of force networks in tapped particulate systems of disks and pentagons. ii. persistence analysis.  Physical Review E, 93(6):062903, 2016.
  * [KKZ+20] Kwangho Kim, Jisu Kim, Manzil Zaheer, Joon Kim, Frederic Chazal, and Larry Wasserman.  Pllay: Efficient topological layer based on persistent landscapes.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 15965–15977. Curran Associates, Inc., 2020.
  * [KM13] Matthew Kahle and Elizabeth Meckes.  Limit theorems for Betti numbers of random simplicial complexes.  Homology, Homotopy and Applications, 15(1):343–374, 2013.
  * [Mil63] J.W. Milnor.  Morse Theory.  Annals of mathematics studies. Princeton University Press, 1963.
  * [Mor59] Marston Morse.  Topologically non-degenerate functions on a compact n-manifold m.  Journal d’Analyse Mathématique, 7(1):189–208, 1959.
  * [MV21] Atish Mitra and Žiga Virk.  The space of persistence diagrams on $n$ points coarsely embeds into Hilbert space.  Proceedings of the American Mathematical Society, 149(6):2693–2703, 2021.
  * [Owa22] Takashi Owada.  Convergence of persistence diagram in the sparse regime.  The Annals of Applied Probability, 32(6):4706–4736, 2022.
  * [RHBK15] Jan Reininghaus, Stefan Huber, Ulrich Bauer, and Roland Kwitt.  A stable multi-scale kernel for topological machine learning.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4741–4748, 2015.
  * [RYB+20] Bastian Rieck, Tristan Yates, Christian Bock, Karsten Borgwardt, Guy Wolf, Nicholas Turk-Browne, and Smita Krishnaswamy.  Uncovering the topology of time-varying fMRI data using cubical persistence.  Advances in neural information processing systems, 33:6900–6912, 2020.
  * [Sch21] Benjamin Schweinhart.  Persistent homology and the upper box dimension.  Discrete & Computational Geometry, 65(2):331–364, 2021.
  * [SL22] Yara Skaf and Reinhard Laubenbacher.  Topological data analysis in biomedicine: A review.  Journal of Biomedical Informatics, 130:104082, 2022.
  * [SPCGB19] Ann E Sizemore, Jennifer E Phillips-Cremins, Robert Ghrist, and Danielle S Bassett.  The importance of the whole: topological data analysis for the network neuroscientist.  Network Neuroscience, 3(3):656–673, 2019.
  * [SSDE20] Umut Simsekli, Ozan Sener, George Deligiannidis, and Murat A Erdogdu.  Hausdorff dimension, heavy tails, and generalization in neural networks.  Advances in Neural Information Processing Systems, 33:5138–5151, 2020.
  * [STR+17] Mohammad Saadatfar, Hiroshi Takeuchi, Vanessa Robins, Nicolas Francois, and Yasuaki Hiraoka.  Pore configuration landscape of granular crystallization.  Nature communications, 8(1):15082, 2017.
  * [STR+18] Anirudh Som, Kowshik Thopalli, Karthikeyan Natesan Ramamurthy, Vinay Venkataraman, Ankita Shukla, and Pavan Turaga.  Perturbation robust representations of topological persistence diagrams.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 617–635, 2018.
  * [SYM23] Anna Song, Ka Man Yim, and Anthea Monod.  Generalized morse theory of distance functions to surfaces for persistent homology.  arXiv preprint arXiv:2306.14716, 2023.
  * [VDVW97] Aad W Van Der Vaart and Jon A Wellner.  Weak convergence and empirical processes: with applications to statistics.  Springer New York, 1997.
  * [Ver18] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [Wag21] Alexander Wagner.  Nonembeddability of persistence diagrams with $p>2$ Wasserstein metric.  Proceedings of the American Mathematical Society, 149(6):2673–2677, 2021.
  * [WB19] Jonathan Weed and Francis Bach.  Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance.  Bernoulli, 25(4A):2620–2648, 2019.


# Graph Neural Networks for Emulation of Finite-Element
Ice Dynamics in Greenland and Antarctic Ice Sheets

Younghyun Koo  Maryam Rahnemoonfar

###### Abstract

Although numerical models provide accurate solutions for ice sheet dynamics
based on physics laws, they accompany intensified computational demands to
solve partial differential equations. In recent years, convolutional neural
networks (CNNs) have been widely used as statistical emulators for those
numerical models. However, since CNNs operate on regular grids, they cannot
represent the refined meshes and computational efficiency of finite-element
numerical models. Therefore, instead of CNNs, this study adopts an equivariant
graph convolutional network (EGCN) as an emulator for the ice sheet dynamics
modeling. EGCN reproduces ice thickness and velocity changes in the Helheim
Glacier, Greenland, and Pine Island Glacier, Antarctica, with 260 times and 44
times faster computation time, respectively. Compared to the traditional CNN
and graph convolutional network, EGCN shows outstanding accuracy in thickness
prediction near fast ice streams by preserving the equivariance to the
translation and rotation of graphs.

Ice-sheet and Sea-level System Model (ISSM), finite element analysis (FEA),
surrogate model, graph convolutional network (GCN), graph attention network
(GAT), climate change, von Mises tensile stress

##  1 Introduction

As the mass losses in Greenland and Antarctic ice sheets have accelerated in
recent years (Otosaka et al., 2023), it becomes more important to accurately
model the ice sheets’ behavior and their interaction with climate change.
Although several numerical models have been proposed based on the physical
understanding of ice flow dynamics (e.g., Stokes equations), those numerical
models require intensive computation to solve complex partial differential
equations (PDEs). In order to emulate such computationally expensive numerical
models, computationally cheap machine learning models have been developed
leveraging the parallel processing capability of graphic processing units
(GPUs).

Most of them relied on convolutional neural network (CNN) architectures
(Jouvet et al., 2022; Jouvet & Cordonnier, 2023) that are only compatible with
regular Euclidean or grid-like structures, such as images. However, given that
finite-element numerical models are implemented on unstructured meshes rather
than regular grids, CNN cannot fully represent unstructured meshes of
numerical ice sheet models. Therefore, instead of traditional CNN
architecture, we choose graph neural network (GNN) as the backbone deep
learning architecture to emulate finite-element ice sheet models. Unlike CNN,
GNN can operate on any irregular non-Euclidean graph structures (i.e., any
data structures with nodes and edges) by updating node features iteratively
through message-passing processes between neighboring nodes (Zhang et al.,
2019). In particular, considering the dynamical behavior of ice sheets, we
adopt a special graph convolutional network (GCN) architecture designed to
preserve equivariance to rotation and translations of dynamic systems: so-
called equivariant graph convolutional network (EGCN) (Satorras et al., 2022).

In this study, we aim to develop GNN emulators, particularly EGCN, for the
Ice-sheet Sea-level System Model (ISSM) (Larour et al., 2012). The
computational efficiency of the ISSM numerical model results from adaptive
mesh refinement (AMR), which allocates computational resources depending on
the expected precision of ice velocity for individual finite elements (Larour
et al., 2012; dos Santos et al., 2021). We take the Helheim Glacier, Greenland
(Fig. 1a), and Pine Island Glacier (PIG), Antarctica (Fig. 1b), as our testing
sites because they are the representative ice sheets that have experienced
rapid acceleration and mass loss. Given that mass loss of Helheim Glacier is
primarily driven by calving (Choi et al., 2018; Cheng et al., 2022) and PIG is
driven by basal melting (Joughin et al., 2021; Jacobs et al., 2011), we
predict how the ice sheet dynamics would change by calving and melting
parameters.

We train GNN models using the simulation data acquired from the ISSM numerical
model and assess their fidelity and computational efficiency for modeling ice
dynamics in the Helheim Glacier and PIG. The main contributions of this
research are the following:

  * •

We adopt EGCN as a statistical emulator to reproduce the ice dynamics
simulated from ISSM. This study is the first application of EGCN on ice sheet
modeling.

  * •

We conduct extensive experiments for the Helheim Glacier and PIG to evaluate
the potential and superiority of EGCN architecture over traditional GCN and
CNN architectures in predicting ice dynamics.

Figure 1: (a) Location of the Helheim Glacier in Greenland; (b) Location of
the Pine Island Glacier (PIG) in Antarctica.

##  2 Methodology

The unstructured meshes of ISSM can be regarded as graph structures where the
connectivity between nodes is expressed via adjacency matrices. Additionally,
since the ice movements change the spatial domains of ice sheets, modeling of
ice sheet dynamics can be regarded as a dynamic graph problem. Hence, based on
the graph representation of the ISSM meshes, we adopt and test the EGCN, which
is specialized for the modeling of dynamically changing graph structures
(Satorras et al., 2022). We also test two baseline deep learning models:
normal GCN and CNN.

###  2.1 Equivariant Graph Convolutional Network (EGCN)

EGCN is designed to conserve equivariance to rotations, translations,
reflections, and permutations in a graph structure (Satorras et al., 2022):
this EGCN structure has shown greater generalizability to any graph structures
of dynamics systems. Let the $l$th graph convolutional layer receives a set of
node features
$\textbf{h}^{(l)}=\\{h_{1}^{(l)},h_{2}^{(l)},...,h_{N}^{(l)}\\}$,
$h_{i}^{(l)}\in\mathbb{R}^{F_{l}}$ as the input and produces a new set of node
features,
$\textbf{h}^{(l+1)}=\\{h_{1}^{(l+1)},h_{2}^{(l+1)},...,h_{N}^{(l+1)}\\}$,
$h_{i}^{(l+1)}\in\mathbb{R}^{F_{l+1}}$, for the next $l+1$th layer. $N$ is the
number of nodes; $F_{l}$ and $F_{l+1}$ is the number of features in each node
at $l$th layer and $l+1$ layer, respectively. Then, an equivariant graph
convolutional layer updates node features by using the following equations:

$m_{ij}=\phi_{e}(h_{i}^{(l)},h_{j}^{(l)},||x_{i}^{(l)}-x_{j}^{(l)}||^{2},a_{ij})$
(1) $x_{i}^{(l+1)}=x_{i}^{(l)}+C\sum_{j\neq
i}(x_{i}^{(l)}-x_{j}^{(l)})\phi_{x}(m_{ij})$ (2)
$m_{i}=x_{i}^{(l)}+C\sum_{j\neq i}m_{ij}$ (3)
$h_{i}^{(l+1)}=\phi_{h}(h_{i}^{(l)},m_{i})$ (4)

where $a_{ij}$ is the edge attributes, $x_{i}$ and $x_{j}$ are the coordinate
embeddings for node $i$ and $j$, respectively, and $C$ is a constant for
normalization computed as $1/|\mathcal{N}(i)|$. $\mathcal{N}(i)$ is the set of
neighbors of node $i$. For the edge attributes, we extract five attributes
from the connecting nodes: distance, surface slope, base slope, acceleration
of x-component velocity, and acceleration of y-component velocity. $\phi_{e}$,
$\phi_{x}$, and $\phi_{h}$ are the edge, position, and node operations,
respectively, which are approximated by single-layer MLPs with 128 hidden
features. Herein, we regard the x- and y-component ice velocities as the
displacement causing the coordinate changes of graphs (i.e., $x_{i}$ and
$x_{j}$), and ice thickness is represented by hidden features (i.e.,
$h_{i}^{(l)}$). Thus, the ice thickness is equivariant to the displacement
caused by ice flow. Our EGCN consists of one input layer, five equivariant
graph convolutional layers, and one output layer (Fig. 2).

Figure 2: Schematic illustration of GCN architectures.

###  2.2 Graph Convolutional Network (GCN)

We compare the EGCN with a normal GCN architecture proposed by Kipf & Welling
(2017). Our GCN consists of one input layer, five graph convolutional hidden
layers, and one output layer (Fig. 2). For each graph convolutional layer, the
number of features is set to 128. The weights of graph convolutional layers
are updated via the layer-wise propagation rule as follows:

$\begin{split}&h_{i}^{(l+1)}=\sigma(\sum_{j\in\mathcal{N}(i)}\frac{1}{c_{ij}}\textbf{W}^{(l)}h_{j}^{(l)})\\\
\end{split}$ (5)

where $c_{ij}$ is the product of the square root of node degrees (i.e.,
$c_{ij}=\sqrt{|\mathcal{N}(j|)}\sqrt{|\mathcal{N}(i)|}$), $\textbf{W}^{(l)}$
is a layer-specific trainable weight matrix
($\textbf{W}^{(l)}\in\mathbb{R}^{F_{l+1}\times F_{l}}$), and $\sigma(\cdot)$
is an activation function; we use the Leaky ReLU activation function of 0.01
negative slope in this study.

###  2.3 Convolutional Neural Network (CNN)

As a baseline model to be compared with GNN emulators, we train and test a
fully convolutional network (FCN), which has a similar architecture to Jouvet
et al. (2022). In the FCN architecture, the hidden graph convolutional layers
are replaced with convolutional layers; it consists of five hidden
convolutional layers. Each convolution has a 3$\times$3 kernel size and 128
features. The leaky ReLU activation function of 0.01 negative slope is applied
after each convolutional layer. Since the FCN takes regular grids as the input
and output, we interpolate all the irregular mesh construction of the ISSM
into a 1 km regular grid using bilinear interpolation.

##  3 Experiment for the Helheim Glacier

In this study, we apply the EGCN, GCN, and FCN emulators to predict the ice
sheet dynamics in the Helheim Glacier varying by calving threshold parameter
($\sigma_{\text{max}}$) of the von Mises (VM) calving law (Morlighem et al.,
2016).

###  3.1 Preparation of training data

To generate training datasets for the GNN emulators, we run transient
simulations of ice dynamics in the Helheim Glacier between 2007 and 2020. We
use the Shelfy-Stream Approximation (SSA) (MacAyeal, 1989) to explain ice flow
in the Helheim Glacier (Cheng et al., 2022; Choi et al., 2018). For the
initial condition of the model, we use the surface velocities from satellite
observations (Mouginot et al., 2017a, 2019), bed topography from BedMachine
Greenland v6 (Morlighem et al., 2017), surface mass balance (SMB) from the
Regional Atmosphere Model (Tedesco & Fettweis, 2020), and the ocean thermal
forcing from the reanalysis data (Wood et al., 2021). To examine the
sensitivity of ice dynamics to the calving parameter $\sigma_{\text{max}}$, we
implement the transient solutions for 7 different $\sigma_{\text{max}}$ values
(i.e., 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.0 MPa) based on the values
proposed by Choi et al. (2018).

We convert the finite-element simulation results into graph structures by
extracting adjacent matrices that represent the connectivity between nodes:
for a triangular mesh composed of three nodes, these nodes are connected to
each other. Our GNN takes 10 features of graph nodes as inputs, including
$\sigma_{\text{max}}$, time, SMB, initial x- and y-component ice velocities,
initial surface elevation, bed elevation, initial ice thickness, and initial
ice mask. Then, the output layer predicts 4 features of the graph nodes:
x-component ice velocity, y-component ice velocity, ice thickness, and ice
mask. All those input and output features are normalized into [-1, 1] for
stable learning using the nominal maximum and minimum values of those
variables.

A total of 1,827 graphs are divided into the train, validation, and test
datasets based on the $\sigma_{\text{max}}$ values to assess if our emulator
can be generalized for out-of-sample $\sigma_{\text{max}}$ values. The data
with $\sigma_{\text{max}}$ of 0.70, 0.80, 0.85, 0.90, and 1.0 MPa are used for
training and validation: we randomly split them into 70 % and 30 % for
training and testing, respectively. The rest of the data with
$\sigma_{\text{max}}$ 0.75 and 0.95 MPa are test datasets. Consequently, the
number of training, validation, and test datasets is 913, 392, and 522,
respectively. The model is optimized by Adam stochastic gradient descent
algorithm with the mean square error (MSE) loss function, 400 epochs, and
0.001 learning rate.

###  3.2 Results

Table 1: Test accuracy of ice velocity and thickness for deep learning emulators for the Helheim Glacier. All matrices are averaged for 2 testing $\sigma_{\text{max}}$ values: 75 and 95 MPa. Model | Ice velocity | Ice thickness
---|---|---
RMSE (m/year) | R | RMSE (m) | R
EGCN | 95.66 | 0.990 | 27.55 | 0.999
GCN | 104.65 | 0.989 | 48.22 | 0.998
FCN | 168.33 | 0.972 | 64.32 | 0.996

The overall test accuracy of GCN, EGCN, and FCN is shown in Table 1. All deep-
learning emulators show outstanding performances in ice velocity and thickness
prediction with R-values greater than 0.990 in most cases. EGCN shows the best
accuracy among them, with around 10 m/year of lower ice velocity RMSE and 20 m
of lower ice thickness RMSE than normal GCN. Whereas GCN uses the spatial
distance of neighboring nodes in determining their relative weights in the
propagation process (Eq. 5), EGCN uses the message passing from all nodes to
preserve the equivariance of the entire graph. The equivariance architecture
to the graph rotation and translation could guarantee more generalizability to
any graph conditions, leading to the improvement of model fidelity. In
particular, EGCN can potentially represent the ice thickness changes induced
by ice velocity because the spatial embeddings are used to predict hidden
feature embeddings (Eq. 1-4). Additionally, we highlight that GNN emulators
outperform traditional FCN. Since FCN uses the interpolated 1-km regular grid
for all points without any adjustment to ice velocity, FCN inhibits
delineating boundary conditions in a fine resolution, leading to significant
errors. Since ISSM and GNN emulators use irregular meshes, which assign a
finer resolution for fast ice and a coarser resolution for static ice, they
can describe the ice boundaries more accurately with fine resolutions.

##  4 Experiments for the Pine Island Glacier

We applied the same GCN, EGCN, and FCN architectures to the PIG in Antarctica.
Although the same input and output features are used for the PIG emulator,
only the calving parameter $\sigma_{\text{{max}}}$ is replaced with basal
melting rate because the ice dynamics in PIG are likely driven by melt-driven
thinning near the grounding line rather than calving.

###  4.1 Preparation of training data

Similar to the Helheim experiment, the training datasets are also collected
from the ISSM transient simulation for 20 years using the SSA. We implement
the ISSM simulations for different annual basal melting rate scenarios ranging
from 0 to 70 m/year for every 2 m/year. Consequently, we collect 8,640 graphs
from the experiments with 36 different melting rates. We divide the 8,640
graphs into train, validation, and test datasets based on the melting rate
values: melting rates of 10, 30, 50, and 70 m/year are used as validation
datasets, melting rates of 0, 20, 40, and 60 m/year as test datasets and the
rest are used as training datasets. The total graphs for training, validation,
and testing are 6,720, 960, and 960, respectively. We use the following
observation data as the initial ice conditions: ice velocity from the NASA
Making Earth System Data Records for Use in Research Environments (MEaSUREs)
(Mouginot et al., 2017b), 1 km Antarctic digital elevation model (DEM) (Bamber
et al., 2009), bedrock topography data of the Amundsen Sea continental shelf
(Nitsche et al., 2007), Antarctic surface temperature data (Comiso, 2000), and
Antarctic surface mass balance data (Vaughan et al., 1999).

###  4.2 Results

Table 2 summarizes the accuracy of ice velocity and thickness prediction from
deep learning emulators. All deep learning emulators show significant
performance in both ice velocity and thickness predictions with R greater than
0.997. GNN emulators, including EGCN and GCN, show better accuracy than FCN in
ice velocity prediction, and EGCN shows the best accuracy. EGCN shows a
remarkably low ice thickness RMSE compared to the other emulators, which
agrees with the results from the Helheim Glacier. Therefore, the equivariant
architecture of EGCN helps improve the predictability of ice thickness changes
caused by ice dynamics in the PIG as well.

Table 2: Test accuracy of ice velocity and thickness for deep learning emulators for the PIG. All matrices are averaged for 4 testing melting rates: 0, 20, 40, and 60 m/year. Model | Ice velocity | Ice thickness
---|---|---
RMSE (m/year) | R | RMSE (m) | R
EGCN | 55.29 | 0.997 | 14.75 | 0.999
GCN | 55.99 | 0.997 | 34.61 | 0.999
FCN | 56.69 | 0.997 | 20.18 | 0.999

##  5 Computational performance

We record and compare the time to generate the final transient simulations for
the Helheim Glacier and PIG experiments (Table 3). The computation time of
ISSM is the total elapsed time spent on a single node of the Texas Advanced
Computing Center (TACC) Frontera supercomputing cluster, which is equipped
with 56 cores of Intel 8280 Cascade Lake CPUs (192 GB memory). The computation
times of deep learning emulators are the total elapsed time on a CPU (Intel(R)
Core(TM) i7-11700F; 32 GB memory) and GPU (NVIDIA GeForce RTX 3070; 24GB
memory). In the Helheim experiment, GCN shows the most dramatic speed-up by
around 560 times faster computation time than ISSM, followed by EGCN (260
times speed-up) and FCN (250 times speed-up). In the PIG experiment, GCN and
EGCN show 50 times and 44 times faster computation time than ISSM,
respectively. FCN takes the longest time to reproduce the 20-year ice sheet
dynamics, with only 21 times speed-up compared to the ISSM simulation.

Table 3: Total computational time (in seconds) for producing ice sheet 20-year transient simulations in the PIG for 36 melting rates and training time to train deep learning emulators. Model | Helheim | PIG
---|---|---
CPU | GPU | CPU | GPU
ISSM | 6681.25 | - | 712.96 | -
GCN | 200.58 | 11.89 | 125.05 | 14.18
EGCN | 759.37 | 25.42 | 607.77 | 16.14
FCN | 851.71 | 27.01 | 1932.00 | 33.34

##  6 Conclusions

This study introduces an equivariant graph convolutional network (EGCN) as a
deep learning emulator to reproduce ice dynamics modeled by the Ice-sheet and
Sea-level System Model (ISSM) and compares its performance with traditional
graph convolutional network (GCN) and convolutional neural network (CNN). When
the EGCN is trained with transient simulations in the Helheim Glacier,
Greenland, and Pine Island Glacier (PIG), Antarctica, EGCN shows better
accuracy than the others in predicting ice velocity and thickness by
preserving the equivariance of graph structures. By reducing computational
time 50-500 times compared to the CPU-based numerical models, the EGCN and GCN
emulators will be promising tools to investigate the impacts of
parameterization on future ice behavior, which will contribute to the
improvement of the prediction accuracy of ice sheet mass loss and sea level
rise.

## References

  * Bamber et al. (2009) Bamber, J. L., Gomez-Dans, J. L., and Griggs, J. A.  A new 1 km digital elevation model of the antarctic derived from combined satellite radar and laser data – part 1: Data and methods.  _The Cryosphere_ , 3(1):101–111, 2009.  doi: 10.5194/tc-3-101-2009.  URL https://tc.copernicus.org/articles/3/101/2009/.
  * Cheng et al. (2022) Cheng, G., Morlighem, M., Mouginot, J., and Cheng, D.  Helheim glacier’s terminus position controls its seasonal and inter-annual ice flow variability.  _Geophysical Research Letters_ , 49(5):e2021GL097085, 2022.
  * Choi et al. (2018) Choi, Y., Morlighem, M., Wood, M., and Bondzio, J. H.  Comparison of four calving laws to model greenland outlet glaciers.  _The Cryosphere_ , 12(12):3735–3746, 2018.  doi: 10.5194/tc-12-3735-2018.  URL https://tc.copernicus.org/articles/12/3735/2018/.
  * Comiso (2000) Comiso, J. C.  Variability and trends in antarctic surface temperatures from in situ and satellite infrared measurements.  _Journal of Climate_ , 13(10):1674 – 1696, 2000.  doi: 10.1175/1520-0442(2000)013¡1674:VATIAS¿2.0.CO;2.  URL https://journals.ametsoc.org/view/journals/clim/13/10/1520-0442_2000_013_1674_vatias_2.0.co_2.xml.
  * dos Santos et al. (2021) dos Santos, T. D., Morlighem, M., and Seroussi, H.  Assessment of numerical schemes for transient, finite-element ice flow models using issm v4.18.  _Geoscientific Model Development_ , 14(5):2545–2573, 2021.  doi: 10.5194/gmd-14-2545-2021.  URL https://gmd.copernicus.org/articles/14/2545/2021/.
  * Jacobs et al. (2011) Jacobs, S. S., Jenkins, A., Giulivi, C. F., and Dutrieux, P.  Stronger ocean circulation and increased melting under pine island glacier ice shelf.  _Nature Geoscience_ , 4(8):519–523, 2011.
  * Joughin et al. (2021) Joughin, I., Shapero, D., Smith, B., Dutrieux, P., and Barham, M.  Ice-shelf retreat drives recent pine island glacier speedup.  _Science Advances_ , 7(24):eabg3080, 2021.  doi: 10.1126/sciadv.abg3080.  URL https://www.science.org/doi/abs/10.1126/sciadv.abg3080.
  * Jouvet & Cordonnier (2023) Jouvet, G. and Cordonnier, G.  Ice-flow model emulator based on physics-informed deep learning.  _Journal of Glaciology_ , pp. 1–15, 2023.  doi: 10.1017/jog.2023.73.
  * Jouvet et al. (2022) Jouvet, G., Cordonnier, G., Kim, B., Lüthi, M., Vieli, A., and Aschwanden, A.  Deep learning speeds up ice flow modelling by several orders of magnitude.  _Journal of Glaciology_ , 68(270):651–664, 2022.  doi: 10.1017/jog.2021.120.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks, 2017.
  * Larour et al. (2012) Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.  Continental scale, high order, high spatial resolution, ice sheet modeling using the ice sheet system model (issm).  _Journal of Geophysical Research: Earth Surface_ , 117(F1), 2012.  doi: https://doi.org/10.1029/2011JF002140.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JF002140.
  * MacAyeal (1989) MacAyeal, D. R.  Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream b, antarctica.  _Journal of Geophysical Research: Solid Earth_ , 94(B4):4071–4087, 1989.  doi: https://doi.org/10.1029/JB094iB04p04071.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB094iB04p04071.
  * Morlighem et al. (2016) Morlighem, M., Bondzio, J., Seroussi, H., Rignot, E., Larour, E., Humbert, A., and Rebuffi, S.  Modeling of store gletscher’s calving dynamics, west greenland, in response to ocean thermal forcing.  _Geophysical Research Letters_ , 43(6):2659–2666, 2016.  doi: https://doi.org/10.1002/2016GL067695.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016GL067695.
  * Morlighem et al. (2017) Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., O’Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.  Bedmachine v3: Complete bed topography and ocean bathymetry mapping of greenland from multibeam echo sounding combined with mass conservation.  _Geophysical Research Letters_ , 44(21):11,051–11,061, 2017.  doi: https://doi.org/10.1002/2017GL074954.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL074954.
  * Mouginot et al. (2017a) Mouginot, J., Rignot, E., Scheuchl, B., and Millan, R.  Comprehensive annual ice sheet velocity mapping using landsat-8, sentinel-1, and radarsat-2 data.  _Remote Sensing_ , 9(4):364, 2017a.
  * Mouginot et al. (2017b) Mouginot, J., Scheuchl, B., and Rignot, E.  Measures annual antarctic ice velocity maps, version 1, 2017b.  URL https://nsidc.org/data/NSIDC-0720/versions/1.
  * Mouginot et al. (2019) Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.  Forty-six years of greenland ice sheet mass balance from 1972 to 2018.  _Proceedings of the National Academy of Sciences_ , 116(19):9239–9244, 2019.  doi: 10.1073/pnas.1904242116.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1904242116.
  * Nitsche et al. (2007) Nitsche, F. O., Jacobs, S. S., Larter, R. D., and Gohl, K.  Bathymetry of the amundsen sea continental shelf: Implications for geology, oceanography, and glaciology.  _Geochemistry, Geophysics, Geosystems_ , 8(10), 2007.  doi: https://doi.org/10.1029/2007GC001694.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2007GC001694.
  * Otosaka et al. (2023) Otosaka, I. N., Shepherd, A., Ivins, E. R., Schlegel, N.-J., Amory, C., van den Broeke, M. R., Horwath, M., Joughin, I., King, M. D., Krinner, G., Nowicki, S., Payne, A. J., Rignot, E., Scambos, T., Simon, K. M., Smith, B. E., Sørensen, L. S., Velicogna, I., Whitehouse, P. L., A, G., Agosta, C., Ahlstrøm, A. P., Blazquez, A., Colgan, W., Engdahl, M. E., Fettweis, X., Forsberg, R., Gallée, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B. C., Harig, C., Helm, V., Khan, S. A., Kittel, C., Konrad, H., Langen, P. L., Lecavalier, B. S., Liang, C.-C., Loomis, B. D., McMillan, M., Melini, D., Mernild, S. H., Mottram, R., Mouginot, J., Nilsson, J., Noël, B., Pattle, M. E., Peltier, W. R., Pie, N., Roca, M., Sasgen, I., Save, H. V., Seo, K.-W., Scheuchl, B., Schrama, E. J. O., Schröder, L., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T. C., Vishwakarma, B. D., van Wessem, J. M., Wiese, D., van der Wal, W., and Wouters, B.  Mass balance of the greenland and antarctic ice sheets from 1992 to 2020.  _Earth System Science Data_ , 15(4):1597–1616, 2023.  doi: 10.5194/essd-15-1597-2023.  URL https://essd.copernicus.org/articles/15/1597/2023/.
  * Satorras et al. (2022) Satorras, V. G., Hoogeboom, E., and Welling, M.  E(n) equivariant graph neural networks, 2022.
  * Tedesco & Fettweis (2020) Tedesco, M. and Fettweis, X.  Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the greenland ice sheet.  _The Cryosphere_ , 14(4):1209–1223, 2020.
  * Vaughan et al. (1999) Vaughan, D. G., Bamber, J. L., Giovinetto, M., Russell, J., and Cooper, A. P. R.  Reassessment of net surface mass balance in antarctica.  _Journal of Climate_ , 12(4):933 – 946, 1999.  doi: 10.1175/1520-0442(1999)012¡0933:RONSMB¿2.0.CO;2.  URL https://journals.ametsoc.org/view/journals/clim/12/4/1520-0442_1999_012_0933_ronsmb_2.0.co_2.xml.
  * Wood et al. (2021) Wood, M., Rignot, E., Fenty, I., An, L., Bjørk, A., Broeke, M. v. d., Cai, C., Kane, E., Menemenlis, D., Millan, R., Morlighem, M., Mouginot, J., Noël, B., Scheuchl, B., Velicogna, I., Willis, J. K., and Zhang, H.  Ocean forcing drives glacier retreat in Greenland.  _Science Advances_ , 7(1):eaba7282, 2021.  ISSN 2375-2548.  doi: 10.1126/sciadv.aba7282.
  * Zhang et al. (2019) Zhang, S., Tong, H., Xu, J., and Maciejewski, R.  Graph convolutional networks: a comprehensive review.  _Computational Social Networks_ , 6(1):1–23, 2019.